Context Navigation

physiq.F90 @ 1998

Last change on this file since 1998 was 1998, checked in by fhourdin, 10 years ago

Reactivation du calcul d'un zmax continu pour les thermiques

(thermcell_height.F90, thermcell_plume.F90, thermcell_dry.F90)
ouvre la voie à la réactivation d'une fermeture humide des thermiques
iflag_thermals_closure=2
(conf_phys_m.F90, thermcell.h, thermcell_main.F90)

Modification liée à la conservation de l'eau

(add_phys_tend.F90, add_pbl_tend.F90, physiq.F90)

Modifications liées au déclenchement stochastique

possibilité de revenir à la Ale déterministe pour le criter ALE>|CIN| iflag_trig_bl=2, 1 par défaut)
possibilité d'activer une fermeture statistique où ALP est divisé par la probabilité de déclenchement iflag_clos_bl=1 (0 par défaut, ancienne option 1 passée en =2)

Modification de l'entrainemement dans la version "stratocumulus" du

modèle du thermique (quand iflag_thermals_ed=8).
(modifie thermcell_plume.F90)

Catherine, Jean-Yves et Frédéric

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 122.4 KB

Line
1	! $Id: physiq.F90 1998 2014-03-19 14:14:30Z fhourdin $
2	!#define IO_DEBUG
3
4	SUBROUTINE physiq (nlon,nlev, &
5	debut,lafin,jD_cur, jH_cur,pdtphys, &
6	paprs,pplay,pphi,pphis,presnivs,clesphy0, &
7	u,v,t,qx, &
8	flxmass_w, &
9	d_u, d_v, d_t, d_qx, d_ps &
10	, dudyn &
11	, PVteta)
12
13	USE ioipsl, only: histbeg, histvert, histdef, histend, histsync, &
14	histwrite, ju2ymds, ymds2ju, ioget_year_len
15	USE comgeomphy
16	USE phys_cal_mod
17	USE write_field_phy
18	USE dimphy
19	USE infotrac
20	USE mod_grid_phy_lmdz
21	USE mod_phys_lmdz_para
22	USE iophy
23	USE misc_mod, mydebug=>debug
24	USE vampir
25	USE pbl_surface_mod, ONLY : pbl_surface
26	USE change_srf_frac_mod
27	USE surface_data, ONLY : type_ocean, ok_veget, ok_snow
28	USE phys_local_var_mod ! Variables internes non sauvegardees de la physique
29	USE phys_state_var_mod ! Variables sauvegardees de la physique
30	USE phys_output_var_mod ! Variables pour les ecritures des sorties
31	USE phys_output_write_mod
32	USE fonte_neige_mod, ONLY : fonte_neige_get_vars
33	USE phys_output_mod
34	USE phys_output_ctrlout_mod
35	USE iophy
36	use open_climoz_m, only: open_climoz ! ozone climatology from a file
37	use regr_pr_av_m, only: regr_pr_av
38	use netcdf95, only: nf95_close
39	!IM for NMC files
40	! use netcdf, only: nf90_fill_real
41	use netcdf
42	use mod_phys_lmdz_mpi_data, only: is_mpi_root
43	USE aero_mod
44	use ozonecm_m, only: ozonecm ! ozone of J.-F. Royer
45	use conf_phys_m, only: conf_phys
46	use radlwsw_m, only: radlwsw
47	use phyaqua_mod, only: zenang_an
48	USE control_mod
49	#ifdef REPROBUS
50	USE CHEM_REP, ONLY : Init_chem_rep_xjour
51	#endif
52	USE indice_sol_mod
53	USE phytrac_mod, ONLY : phytrac
54
55	!IM stations CFMIP
56	USE CFMIP_point_locations
57	use FLOTT_GWD_rando_m, only: FLOTT_GWD_rando
58
59	IMPLICIT none
60	!>======================================================================
61	!!
62	!! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
63	!!
64	!! Objet: Moniteur general de la physique du modele
65	!!AA Modifications quant aux traceurs :
66	!!AA - uniformisation des parametrisations ds phytrac
67	!!AA - stockage des moyennes des champs necessaires
68	!!AA en mode traceur off-line
69	!!======================================================================
70	!! CLEFS CPP POUR LES IO
71	!! =====================
72	#define histNMC
73	!!======================================================================
74	!! modif ( P. Le Van , 12/10/98 )
75	!!
76	!! Arguments:
77	!!
78	!! nlon----input-I-nombre de points horizontaux
79	!! nlev----input-I-nombre de couches verticales, doit etre egale a klev
80	!! debut---input-L-variable logique indiquant le premier passage
81	!! lafin---input-L-variable logique indiquant le dernier passage
82	!! jD_cur -R-jour courant a l'appel de la physique (jour julien)
83	!! jH_cur -R-heure courante a l'appel de la physique (jour julien)
84	!! pdtphys-input-R-pas d'integration pour la physique (seconde)
85	!! paprs---input-R-pression pour chaque inter-couche (en Pa)
86	!! pplay---input-R-pression pour le mileu de chaque couche (en Pa)
87	!! pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
88	!! pphis---input-R-geopotentiel du sol
89	!! presnivs-input_R_pressions approximat. des milieux couches ( en PA)
90	!! u-------input-R-vitesse dans la direction X (de O a E) en m/s
91	!! v-------input-R-vitesse Y (de S a N) en m/s
92	!! t-------input-R-temperature (K)
93	!! qx------input-R-humidite specifique (kg/kg) et d'autres traceurs
94	!! d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
95	!! d_q_dyn-input-R-tendance dynamique pour "q" (kg/kg/s)
96	!! flxmass_w -input-R- flux de masse verticale
97	!! d_u-----output-R-tendance physique de "u" (m/s/s)
98	!! d_v-----output-R-tendance physique de "v" (m/s/s)
99	!! d_t-----output-R-tendance physique de "t" (K/s)
100	!! d_qx----output-R-tendance physique de "qx" (kg/kg/s)
101	!! d_ps----output-R-tendance physique de la pression au sol
102	!!IM
103	!! PVteta--output-R-vorticite potentielle a des thetas constantes
104	!!======================================================================
105	include "dimensions.h"
106	integer jjmp1
107	parameter (jjmp1=jjm+1-1/jjm)
108	integer iip1
109	parameter (iip1=iim+1)
110
111	include "regdim.h"
112	include "dimsoil.h"
113	include "clesphys.h"
114	include "temps.h"
115	include "iniprint.h"
116	include "thermcell.h"
117	!======================================================================
118	LOGICAL ok_cvl ! pour activer le nouveau driver pour convection KE
119	PARAMETER (ok_cvl=.TRUE.)
120	LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface
121	PARAMETER (ok_gust=.FALSE.)
122	integer iflag_radia ! active ou non le rayonnement (MPL)
123	save iflag_radia
124	!$OMP THREADPRIVATE(iflag_radia)
125	!======================================================================
126	LOGICAL check ! Verifier la conservation du modele en eau
127	PARAMETER (check=.FALSE.)
128	LOGICAL ok_stratus ! Ajouter artificiellement les stratus
129	PARAMETER (ok_stratus=.FALSE.)
130	!======================================================================
131	REAL amn, amx
132	INTEGER igout
133	!======================================================================
134	! Clef controlant l'activation du cycle diurne:
135	!cc LOGICAL cycle_diurne
136	!cc PARAMETER (cycle_diurne=.FALSE.)
137	!======================================================================
138	! Modele thermique du sol, a activer pour le cycle diurne:
139	!cc LOGICAL soil_model
140	!cc PARAMETER (soil_model=.FALSE.)
141	!======================================================================
142	! Dans les versions precedentes, l'eau liquide nuageuse utilisee dans
143	! le calcul du rayonnement est celle apres la precipitation des nuages.
144	! Si cette cle new_oliq est activee, ce sera une valeur moyenne entre
145	! la condensation et la precipitation. Cette cle augmente les impacts
146	! radiatifs des nuages.
147	!cc LOGICAL new_oliq
148	!cc PARAMETER (new_oliq=.FALSE.)
149	!======================================================================
150	! Clefs controlant deux parametrisations de l'orographie:
151	!c LOGICAL ok_orodr
152	!cc PARAMETER (ok_orodr=.FALSE.)
153	!cc LOGICAL ok_orolf
154	!cc PARAMETER (ok_orolf=.FALSE.)
155	!======================================================================
156	LOGICAL ok_journe ! sortir le fichier journalier
157	save ok_journe
158	!$OMP THREADPRIVATE(ok_journe)
159	!
160	LOGICAL ok_mensuel ! sortir le fichier mensuel
161	save ok_mensuel
162	!$OMP THREADPRIVATE(ok_mensuel)
163	!
164	LOGICAL ok_instan ! sortir le fichier instantane
165	save ok_instan
166	!$OMP THREADPRIVATE(ok_instan)
167	!
168	LOGICAL ok_LES ! sortir le fichier LES
169	save ok_LES
170	!$OMP THREADPRIVATE(ok_LES)
171	!
172	LOGICAL callstats ! sortir le fichier stats
173	save callstats
174	!$OMP THREADPRIVATE(callstats)
175	!
176	LOGICAL ok_region ! sortir le fichier regional
177	PARAMETER (ok_region=.FALSE.)
178	!======================================================================
179	real seuil_inversion
180	save seuil_inversion
181	!$OMP THREADPRIVATE(seuil_inversion)
182	integer iflag_ratqs
183	save iflag_ratqs
184	!$OMP THREADPRIVATE(iflag_ratqs)
185	real facteur
186
187	REAL wmax_th(klon)
188	REAL tau_overturning_th(klon)
189
190	integer lmax_th(klon)
191	integer limbas(klon)
192	real ratqscth(klon,klev)
193	real ratqsdiff(klon,klev)
194	real zqsatth(klon,klev)
195
196	!======================================================================
197	!
198	INTEGER ivap ! indice de traceurs pour vapeur d'eau
199	PARAMETER (ivap=1)
200	INTEGER iliq ! indice de traceurs pour eau liquide
201	PARAMETER (iliq=2)
202
203	!
204	!
205	! Variables argument:
206	!
207	INTEGER nlon
208	INTEGER nlev
209	REAL, intent(in):: jD_cur, jH_cur
210
211	REAL pdtphys
212	LOGICAL debut, lafin
213	REAL paprs(klon,klev+1)
214	REAL pplay(klon,klev)
215	REAL pphi(klon,klev)
216	REAL pphis(klon)
217	REAL presnivs(klev)
218	REAL znivsig(klev)
219	real pir
220
221	REAL u(klon,klev)
222	REAL v(klon,klev)
223	REAL t(klon,klev),thetal(klon,klev)
224	! thetal: ligne suivante a decommenter si vous avez les fichiers MPL 20130625
225	! fth_fonctions.F90 et parkind1.F90
226	! sinon thetal=theta
227	! REAL fth_thetae,fth_thetav,fth_thetal
228	REAL qx(klon,klev,nqtot)
229	REAL flxmass_w(klon,klev)
230	REAL d_u(klon,klev)
231	REAL d_v(klon,klev)
232	REAL d_t(klon,klev)
233	REAL d_qx(klon,klev,nqtot)
234	REAL d_ps(klon)
235	! Variables pour le transport convectif
236	real da(klon,klev),phi(klon,klev,klev),mp(klon,klev)
237	! Variables pour le lessivage convectif
238	! RomP >>>
239	real phi2(klon,klev,klev)
240	real d1a(klon,klev),dam(klon,klev)
241	real ev(klon,klev),ep(klon,klev)
242	real clw(klon,klev),elij(klon,klev,klev)
243	real epmlmMm(klon,klev,klev),eplaMm(klon,klev)
244	! RomP <<<
245	!IM definition dynamique o_trac dans phys_output_open
246	! type(ctrl_out) :: o_trac(nqtot)
247	!
248	!IM Amip2 PV a theta constante
249	!
250	INTEGER nbteta
251	PARAMETER(nbteta=3)
252	CHARACTER*3 ctetaSTD(nbteta)
253	DATA ctetaSTD/'350','380','405'/
254	SAVE ctetaSTD
255	!$OMP THREADPRIVATE(ctetaSTD)
256	REAL rtetaSTD(nbteta)
257	DATA rtetaSTD/350., 380., 405./
258	SAVE rtetaSTD
259	!$OMP THREADPRIVATE(rtetaSTD)
260	!
261	REAL PVteta(klon,nbteta)
262	!
263	!MI Amip2 PV a theta constante
264
265	!ym INTEGER klevp1, klevm1
266	!ym PARAMETER(klevp1=klev+1,klevm1=klev-1)
267	!ym include "raddim.h"
268	!
269	!
270	!IM Amip2
271	! variables a une pression donnee
272	!
273	include "declare_STDlev.h"
274	!
275	!
276	include "radopt.h"
277	!
278	!
279
280
281	INTEGER debug
282	INTEGER n
283	!ym INTEGER npoints
284	!ym PARAMETER(npoints=klon)
285	!
286	INTEGER nregISCtot
287	PARAMETER(nregISCtot=1)
288	!
289	! imin_debut, nbpti, jmin_debut, nbptj : parametres pour sorties sur 1 region rectangulaire
290	! y compris pour 1 point
291	! imin_debut : indice minimum de i; nbpti : nombre de points en direction i (longitude)
292	! jmin_debut : indice minimum de j; nbptj : nombre de points en direction j (latitude)
293	INTEGER imin_debut, nbpti
294	INTEGER jmin_debut, nbptj
295	!IM: region='3d' <==> sorties en global
296	CHARACTER*3 region
297	PARAMETER(region='3d')
298	logical ok_hf
299	!
300	save ok_hf
301	!$OMP THREADPRIVATE(ok_hf)
302
303	INTEGER longcles
304	PARAMETER ( longcles = 20 )
305	REAL clesphy0( longcles )
306	!
307	! Variables propres a la physique
308	INTEGER itap
309	SAVE itap ! compteur pour la physique
310	!$OMP THREADPRIVATE(itap)
311	!
312	REAL,save :: solarlong0
313	!$OMP THREADPRIVATE(solarlong0)
314
315	!
316	! Parametres de l'Orographie a l'Echelle Sous-Maille (OESM):
317	!
318	!IM 141004 REAL zulow(klon),zvlow(klon),zustr(klon), zvstr(klon)
319	REAL zulow(klon),zvlow(klon)
320	!
321	INTEGER igwd,idx(klon),itest(klon)
322	!
323	! REAL,allocatable,save :: run_off_lic_0(:)
324	!!$OMP THREADPRIVATE(run_off_lic_0)
325	!ym SAVE run_off_lic_0
326	!KE43
327	! Variables liees a la convection de K. Emanuel (sb):
328	!
329	REAL bas, top ! cloud base and top levels
330	SAVE bas
331	SAVE top
332	!$OMP THREADPRIVATE(bas, top)
333
334	!
335	!=================================================================================================
336	!CR04.12.07: on ajoute les nouvelles variables du nouveau schema de convection avec poches froides
337	! Variables li\'ees \`a la poche froide (jyg)
338
339	REAL mip(klon,klev) ! mass flux shed by the adiab ascent at each level
340	!
341	REAL wape_prescr, fip_prescr
342	INTEGER it_wape_prescr
343	SAVE wape_prescr, fip_prescr, it_wape_prescr
344	!$OMP THREADPRIVATE(wape_prescr, fip_prescr, it_wape_prescr)
345	!
346	! variables supplementaires de concvl
347	REAL Tconv(klon,klev)
348	REAL sij(klon,klev,klev)
349
350	real, save :: alp_bl_prescr=0.
351	real, save :: ale_bl_prescr=0.
352
353	real, save :: ale_max=1000.
354	real, save :: alp_max=2.
355
356	real, save :: wake_s_min_lsp=0.1
357
358	!$OMP THREADPRIVATE(alp_bl_prescr,ale_bl_prescr)
359	!$OMP THREADPRIVATE(ale_max,alp_max)
360	!$OMP THREADPRIVATE(wake_s_min_lsp)
361
362
363	real ok_wk_lsp(klon)
364
365	!RC
366	! Variables li\'ees \`a la poche froide (jyg et rr)
367	! Version diagnostique pour l'instant : pas de r\'etroaction sur la convection
368
369	REAL t_wake(klon,klev),q_wake(klon,klev) ! wake pour la convection
370
371	REAL wake_dth(klon,klev) ! wake : temp pot difference
372
373	REAL wake_d_deltat_gw(klon,klev)! wake : delta T tendency due to Gravity Wave (/s)
374	REAL wake_omgbdth(klon,klev) ! Wake : flux of Delta_Theta transported by LS omega
375	REAL wake_dp_omgb(klon,klev) ! Wake : vertical gradient of large scale omega
376	REAL wake_dtKE(klon,klev) ! Wake : differential heating (wake - unpertubed) CONV
377	REAL wake_dqKE(klon,klev) ! Wake : differential moistening (wake - unpertubed) CONV
378	REAL wake_dtPBL(klon,klev) ! Wake : differential heating (wake - unpertubed) PBL
379	REAL wake_dqPBL(klon,klev) ! Wake : differential moistening (wake - unpertubed) PBL
380	REAL wake_ddeltat(klon,klev),wake_ddeltaq(klon,klev)
381	REAL wake_dp_deltomg(klon,klev) ! Wake : gradient vertical de wake_omg
382	REAL wake_spread(klon,klev) ! spreading term in wake_delt
383	!
384	!pourquoi y'a pas de save??
385	!
386	INTEGER wake_k(klon) ! Wake sommet
387	!
388	REAL t_undi(klon,klev) ! temperature moyenne dans la zone non perturbee
389	REAL q_undi(klon,klev) ! humidite moyenne dans la zone non perturbee
390	!
391	!jyg
392	!cc REAL wake_pe(klon) ! Wake potential energy - WAPE
393
394	REAL wake_gfl(klon) ! Gust Front Length
395	REAL wake_dens(klon)
396	!
397	!
398	REAL dt_dwn(klon,klev)
399	REAL dq_dwn(klon,klev)
400	REAL wdt_PBL(klon,klev)
401	REAL udt_PBL(klon,klev)
402	REAL wdq_PBL(klon,klev)
403	REAL udq_PBL(klon,klev)
404	REAL M_dwn(klon,klev)
405	REAL M_up(klon,klev)
406	REAL dt_a(klon,klev)
407	REAL dq_a(klon,klev)
408	REAL, SAVE :: alp_offset
409	!$OMP THREADPRIVATE(alp_offset)
410
411	!
412	!RR:fin declarations poches froides
413	!=======================================================================================================
414
415	REAL ztv(klon,klev),ztva(klon,klev)
416	REAL zpspsk(klon,klev)
417	REAL ztla(klon,klev),zqla(klon,klev)
418	REAL zthl(klon,klev)
419
420	!cc nrlmd le 10/04/2012
421
422	!--------Stochastic Boundary Layer Triggering: ALE_BL--------
423	!---Propri\'et\'es du thermiques au LCL
424	real zlcl_th(klon) ! Altitude du LCL calcul\'e continument (pcon dans thermcell_main.F90)
425	real fraca0(klon) ! Fraction des thermiques au LCL
426	real w0(klon) ! Vitesse des thermiques au LCL
427	real w_conv(klon) ! Vitesse verticale de grande \'echelle au LCL
428	real tke0(klon,klev+1) ! TKE au début du pas de temps
429	real therm_tke_max0(klon) ! TKE dans les thermiques au LCL
430	real env_tke_max0(klon) ! TKE dans l'environnement au LCL
431
432	!---D\'eclenchement stochastique
433	integer :: tau_trig(klon)
434
435	!--------Statistical Boundary Layer Closure: ALP_BL--------
436	!---Profils de TKE dans et hors du thermique
437	real pbl_tke_input(klon,klev+1,nbsrf)
438	real therm_tke_max(klon,klev) ! Profil de TKE dans les thermiques
439	real env_tke_max(klon,klev) ! Profil de TKE dans l'environnement
440
441
442	!cc fin nrlmd le 10/04/2012
443
444	! Variables locales pour la couche limite (al1):
445	!
446	!Al1 REAL pblh(klon) ! Hauteur de couche limite
447	!Al1 SAVE pblh
448	!34EK
449	!
450	! Variables locales:
451	!
452	!AA
453	!AA Pour phytrac
454	REAL u1(klon) ! vents dans la premiere couche U
455	REAL v1(klon) ! vents dans la premiere couche V
456
457	!@$$ LOGICAL offline ! Controle du stockage ds "physique"
458	!@$$ PARAMETER (offline=.false.)
459	!@$$ INTEGER physid
460	REAL frac_impa(klon,klev) ! fractions d'aerosols lessivees (impaction)
461	REAL frac_nucl(klon,klev) ! idem (nucleation)
462	! RomP >>>
463	REAL beta_prec_fisrt(klon,klev) ! taux de conv de l'eau cond (fisrt)
464	! RomP <<<
465
466	REAL :: calday
467
468	!IM cf FH pour Tiedtke 080604
469	REAL rain_tiedtke(klon),snow_tiedtke(klon)
470	!
471	!IM 050204 END
472	REAL devap(klon) ! evaporation et sa derivee
473	REAL dsens(klon) ! chaleur sensible et sa derivee
474
475	!
476	! Conditions aux limites
477	!
478	!
479	REAL :: day_since_equinox
480	! Date de l'equinoxe de printemps
481	INTEGER, parameter :: mth_eq=3, day_eq=21
482	REAL :: jD_eq
483
484	LOGICAL, parameter :: new_orbit = .true.
485
486	!
487	INTEGER lmt_pas
488	SAVE lmt_pas ! frequence de mise a jour
489	!$OMP THREADPRIVATE(lmt_pas)
490	real zmasse(klon, llm),exner(klon, llm)
491	! (column-density of mass of air in a cell, in kg m-2)
492	real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
493
494	!IM sorties
495	REAL un_jour
496	PARAMETER(un_jour=86400.)
497	INTEGER itapm1 !pas de temps de la physique du(es) mois precedents
498	SAVE itapm1 !mis a jour le dernier pas de temps du mois en cours
499	!$OMP THREADPRIVATE(itapm1)
500	!======================================================================
501	!
502	! Declaration des procedures appelees
503	!
504	EXTERNAL angle ! calculer angle zenithal du soleil
505	EXTERNAL alboc ! calculer l'albedo sur ocean
506	EXTERNAL ajsec ! ajustement sec
507	EXTERNAL conlmd ! convection (schema LMD)
508	!KE43
509	EXTERNAL conema3 ! convect4.3
510	EXTERNAL fisrtilp ! schema de condensation a grande echelle (pluie)
511	!AA
512	EXTERNAL fisrtilp_tr ! schema de condensation a grande echelle (pluie)
513	! ! stockage des coefficients necessaires au
514	! ! lessivage OFF-LINE et ON-LINE
515	EXTERNAL hgardfou ! verifier les temperatures
516	EXTERNAL nuage ! calculer les proprietes radiatives
517	!C EXTERNAL o3cm ! initialiser l'ozone
518	EXTERNAL orbite ! calculer l'orbite terrestre
519	EXTERNAL phyetat0 ! lire l'etat initial de la physique
520	EXTERNAL phyredem ! ecrire l'etat de redemarrage de la physique
521	EXTERNAL suphel ! initialiser certaines constantes
522	EXTERNAL transp ! transport total de l'eau et de l'energie
523	EXTERNAL ecribina ! ecrire le fichier binaire global
524	EXTERNAL ecribins ! ecrire le fichier binaire global
525	EXTERNAL ecrirega ! ecrire le fichier binaire regional
526	EXTERNAL ecriregs ! ecrire le fichier binaire regional
527	!IM
528	EXTERNAL haut2bas !variables de haut en bas
529	EXTERNAL ini_undefSTD !initialise a 0 une variable a 1 niveau de pression
530	EXTERNAL undefSTD !somme les valeurs definies d'1 var a 1 niveau de pression
531	! EXTERNAL moy_undefSTD !moyenne d'1 var a 1 niveau de pression
532	! EXTERNAL moyglo_aire !moyenne globale d'1 var ponderee par l'aire de la maille (moyglo_pondaire)
533	! !par la masse/airetot (moyglo_pondaima) et la vraie masse (moyglo_pondmass)
534	!
535	! Variables locales
536	!
537	REAL rhcl(klon,klev) ! humiditi relative ciel clair
538	REAL dialiq(klon,klev) ! eau liquide nuageuse
539	REAL diafra(klon,klev) ! fraction nuageuse
540	REAL cldliq(klon,klev) ! eau liquide nuageuse
541	!
542	!XXX PB
543	REAL fluxq(klon,klev, nbsrf) ! flux turbulent d'humidite
544	!
545	REAL zxfluxt(klon, klev)
546	REAL zxfluxq(klon, klev)
547	REAL zxfluxu(klon, klev)
548	REAL zxfluxv(klon, klev)
549
550	! Le rayonnement n'est pas calcule tous les pas, il faut donc
551	! sauvegarder les sorties du rayonnement
552	!ym SAVE heat,cool,albpla,topsw,toplw,solsw,sollw,sollwdown
553	!ym SAVE sollwdownclr, toplwdown, toplwdownclr
554	!ym SAVE topsw0,toplw0,solsw0,sollw0, heat0, cool0
555	!
556	INTEGER itaprad
557	SAVE itaprad
558	!$OMP THREADPRIVATE(itaprad)
559	!
560	REAL conv_q(klon,klev) ! convergence de l'humidite (kg/kg/s)
561	REAL conv_t(klon,klev) ! convergence de la temperature(K/s)
562
563	!
564	! REAL zxsnow(klon)
565	REAL zxsnow_dummy(klon)
566	!
567	REAL dist, rmu0(klon), fract(klon)
568	REAL zdtime, zlongi
569	!
570	REAL qcheck
571	REAL z_avant(klon), z_apres(klon), z_factor(klon)
572	LOGICAL zx_ajustq
573	!
574	REAL za, zb
575	REAL zx_t, zx_qs, zdelta, zcor, zlvdcp, zlsdcp
576	real zqsat(klon,klev)
577	INTEGER i, k, iq, ig, j, nsrf, ll, l, iiq
578	REAL t_coup
579	PARAMETER (t_coup=234.0)
580
581	!ym A voir plus tard !!
582	!ym REAL zx_relief(iim,jjmp1)
583	!ym REAL zx_aire(iim,jjmp1)
584	!
585	! Grandeurs de sorties
586	REAL s_capCL(klon)
587	REAL s_oliqCL(klon), s_cteiCL(klon)
588	REAL s_trmb1(klon), s_trmb2(klon)
589	REAL s_trmb3(klon)
590	!KE43
591	! Variables locales pour la convection de K. Emanuel (sb):
592
593	REAL tvp(klon,klev) ! virtual temp of lifted parcel
594	CHARACTER*40 capemaxcels !max(CAPE)
595
596	REAL rflag(klon) ! flag fonctionnement de convect
597	INTEGER iflagctrl(klon) ! flag fonctionnement de convect
598
599	! -- convect43:
600	INTEGER ntra ! nb traceurs pour convect4.3
601	REAL dtvpdt1(klon,klev), dtvpdq1(klon,klev)
602	REAL dplcldt(klon), dplcldr(klon)
603	!? . condm_con(klon,klev),conda_con(klon,klev),
604	!? . mr_con(klon,klev),ep_con(klon,klev)
605	!? . ,sadiab(klon,klev),wadiab(klon,klev)
606	! --
607	!34EK
608	!
609	! Variables du changement
610	!
611	! con: convection
612	! lsc: condensation a grande echelle (Large-Scale-Condensation)
613	! ajs: ajustement sec
614	! eva: evaporation de l'eau liquide nuageuse
615	! vdf: couche limite (Vertical DiFfusion)
616
617	! tendance nulles
618	REAL, dimension(klon,klev):: du0, dv0, dt0, dq0, dql0
619
620	!
621	!********************************************************
622	! declarations
623
624	!********************************************************
625	!IM 081204 END
626	!
627	REAL pen_u(klon,klev), pen_d(klon,klev)
628	REAL pde_u(klon,klev), pde_d(klon,klev)
629	INTEGER kcbot(klon), kctop(klon), kdtop(klon)
630	!
631	REAL ratqsc(klon,klev)
632	real ratqsbas,ratqshaut,tau_ratqs
633	save ratqsbas,ratqshaut,tau_ratqs
634	!$OMP THREADPRIVATE(ratqsbas,ratqshaut,tau_ratqs)
635
636	! Parametres lies au nouveau schema de nuages (SB, PDF)
637	real fact_cldcon
638	real facttemps
639	logical ok_newmicro
640	save ok_newmicro
641	!$OMP THREADPRIVATE(ok_newmicro)
642	!real ref_liq_pi(klon,klev), ref_ice_pi(klon,klev)
643	save fact_cldcon,facttemps
644	!$OMP THREADPRIVATE(fact_cldcon,facttemps)
645
646	integer iflag_cldcon
647	save iflag_cldcon
648	!$OMP THREADPRIVATE(iflag_cldcon)
649	logical ptconv(klon,klev)
650	!IM cf. AM 081204 BEG
651	logical ptconvth(klon,klev)
652	!IM cf. AM 081204 END
653	!
654	! Variables liees a l'ecriture de la bande histoire physique
655	!
656	!======================================================================
657	!
658
659	!
660	integer itau_w ! pas de temps ecriture = itap + itau_phy
661	!
662	!
663	! Variables locales pour effectuer les appels en serie
664	!
665	!IM RH a 2m (la surface)
666	REAL Lheat
667
668	INTEGER length
669	PARAMETER ( length = 100 )
670	REAL tabcntr0( length )
671	!
672	INTEGER ndex2d(iim*jjmp1)
673	!IM
674	!
675	!IM AMIP2 BEG
676	REAL moyglo, mountor
677	!IM 141004 BEG
678	REAL zustrdr(klon), zvstrdr(klon)
679	REAL zustrli(klon), zvstrli(klon)
680	REAL zustrph(klon), zvstrph(klon)
681	REAL zustrhi(klon), zvstrhi(klon)
682	REAL aam, torsfc
683	!IM 141004 END
684	!IM 190504 BEG
685	INTEGER ij, imp1jmp1
686	PARAMETER(imp1jmp1=(iim+1)*jjmp1)
687	!ym A voir plus tard
688	REAL zx_tmp(imp1jmp1), airedyn(iim+1,jjmp1)
689	REAL padyn(iim+1,jjmp1,klev+1)
690	REAL dudyn(iim+1,jjmp1,klev)
691	REAL rlatdyn(iim+1,jjmp1)
692	!IM 190504 END
693	LOGICAL ok_msk
694	REAL msk(klon)
695	!IM
696	REAL airetot, pi
697	!ym A voir plus tard
698	!ym REAL zm_wo(jjmp1, klev)
699	!IM AMIP2 END
700	!
701	REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique
702	REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D
703	REAL zx_tmp_2d(iim,jjmp1)
704	REAL zx_lon(iim,jjmp1), zx_lat(iim,jjmp1)
705	!
706	INTEGER nid_day_seri, nid_ctesGCM
707	SAVE nid_day_seri, nid_ctesGCM
708	!$OMP THREADPRIVATE(nid_day_seri,nid_ctesGCM)
709	!
710	!IM 280405 BEG
711	! INTEGER nid_bilKPins, nid_bilKPave
712	! SAVE nid_bilKPins, nid_bilKPave
713	! !$OMP THREADPRIVATE(nid_bilKPins, nid_bilKPave)
714	!
715	REAL ve_lay(klon,klev) ! transport meri. de l'energie a chaque niveau vert.
716	REAL vq_lay(klon,klev) ! transport meri. de l'eau a chaque niveau vert.
717	REAL ue_lay(klon,klev) ! transport zonal de l'energie a chaque niveau vert.
718	REAL uq_lay(klon,klev) ! transport zonal de l'eau a chaque niveau vert.
719	!
720	INTEGER nhori, nvert
721	REAL zsto
722	REAL zstophy, zout
723
724	real zjulian
725	save zjulian
726	!$OMP THREADPRIVATE(zjulian)
727
728	character*20 modname
729	character*80 abort_message
730	logical ok_sync
731	real date0
732	integer idayref
733
734	! essai writephys
735	integer fid_day, fid_mth, fid_ins
736	parameter (fid_ins = 1, fid_day = 2, fid_mth = 3)
737	integer prof2d_on, prof3d_on, prof2d_av, prof3d_av
738	parameter (prof2d_on = 1, prof3d_on = 2, &
739	prof2d_av = 3, prof3d_av = 4)
740	! Variables liees au bilan d'energie et d'enthalpi
741	REAL ztsol(klon)
742	REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec
743	REAL d_h_vcol_phy
744	REAL fs_bound, fq_bound
745	SAVE d_h_vcol_phy
746	!$OMP THREADPRIVATE(d_h_vcol_phy)
747	REAL zero_v(klon)
748	CHARACTER*40 ztit
749	INTEGER ip_ebil ! PRINT level for energy conserv. diag.
750	SAVE ip_ebil
751	DATA ip_ebil/0/
752	!$OMP THREADPRIVATE(ip_ebil)
753	INTEGER if_ebil ! level for energy conserv. dignostics
754	SAVE if_ebil
755	!$OMP THREADPRIVATE(if_ebil)
756	REAL q2m(klon,nbsrf) ! humidite a 2m
757
758	!IM: t2m, q2m, ustar, u10m, v10m et t2mincels, t2maxcels
759	CHARACTER*40 t2mincels, t2maxcels !t2m min., t2m max
760	CHARACTER*40 tinst, tave, typeval
761	REAL cldtaupi(klon,klev) ! Cloud optical thickness for pre-industrial (pi) aerosols
762
763
764	! Aerosol optical properties
765	CHARACTER*4, DIMENSION(naero_grp) :: rfname
766	REAL, DIMENSION(klon,klev) :: mass_solu_aero ! total mass concentration for all soluble aerosols[ug/m3]
767	REAL, DIMENSION(klon,klev) :: mass_solu_aero_pi ! - " - (pre-industrial value)
768
769	! Parameters
770	LOGICAL ok_ade, ok_aie ! Apply aerosol (in)direct effects or not
771	LOGICAL ok_cdnc ! ok cloud droplet number concentration (O. Boucher 01-2013)
772	REAL bl95_b0, bl95_b1 ! Parameter in Boucher and Lohmann (1995)
773	SAVE ok_ade, ok_aie, ok_cdnc, bl95_b0, bl95_b1
774	!$OMP THREADPRIVATE(ok_ade, ok_aie, ok_cdnc, bl95_b0, bl95_b1)
775	LOGICAL, SAVE :: aerosol_couple ! true : calcul des aerosols dans INCA
776	! false : lecture des aerosol dans un fichier
777	!$OMP THREADPRIVATE(aerosol_couple)
778	INTEGER, SAVE :: flag_aerosol
779	!$OMP THREADPRIVATE(flag_aerosol)
780	LOGICAL, SAVE :: new_aod
781	!$OMP THREADPRIVATE(new_aod)
782	!
783	!--STRAT AEROSOL
784	LOGICAL, SAVE :: flag_aerosol_strat
785	!$OMP THREADPRIVATE(flag_aerosol_strat)
786	!c-fin STRAT AEROSOL
787	!
788	! Declaration des constantes et des fonctions thermodynamiques
789	!
790	LOGICAL,SAVE :: first=.true.
791	!$OMP THREADPRIVATE(first)
792
793	integer, save:: read_climoz ! read ozone climatology
794	! (let it keep the default OpenMP shared attribute)
795	! Allowed values are 0, 1 and 2
796	! 0: do not read an ozone climatology
797	! 1: read a single ozone climatology that will be used day and night
798	! 2: read two ozone climatologies, the average day and night
799	! climatology and the daylight climatology
800
801	integer, save:: ncid_climoz ! NetCDF file containing ozone climatologies
802	! (let it keep the default OpenMP shared attribute)
803
804	real, pointer, save:: press_climoz(:)
805	! (let it keep the default OpenMP shared attribute)
806	! edges of pressure intervals for ozone climatologies, in Pa, in strictly
807	! ascending order
808
809	integer, save:: co3i = 0
810	! time index in NetCDF file of current ozone fields
811	!$OMP THREADPRIVATE(co3i)
812
813	integer ro3i
814	! required time index in NetCDF file for the ozone fields, between 1
815	! and 360
816
817	INTEGER ierr
818	include "YOMCST.h"
819	include "YOETHF.h"
820	include "FCTTRE.h"
821	!IM 100106 BEG : pouvoir sortir les ctes de la physique
822	include "conema3.h"
823	include "fisrtilp.h"
824	include "nuage.h"
825	include "compbl.h"
826	!IM 100106 END : pouvoir sortir les ctes de la physique
827	!
828	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
829	! Declarations pour Simulateur COSP
830	!============================================================
831	real :: mr_ozone(klon,klev)
832
833	!IM sorties fichier 1D paramLMDZ_phy.nc
834	REAL :: zx_tmp_0d(1,1)
835	INTEGER, PARAMETER :: np=1
836	REAL,dimension(klon_glo) :: rlat_glo
837	REAL,dimension(klon_glo) :: rlon_glo
838	REAL gbils(1), gevap(1), gevapt(1), glat(1), gnet0(1), gnet(1)
839	REAL grain(1), gtsol(1), gt2m(1), gprw(1)
840
841	!IM stations CFMIP
842	INTEGER, SAVE :: nCFMIP
843	!$OMP THREADPRIVATE(nCFMIP)
844	INTEGER, PARAMETER :: npCFMIP=120
845	INTEGER, ALLOCATABLE, SAVE :: tabCFMIP(:)
846	REAL, ALLOCATABLE, SAVE :: lonCFMIP(:), latCFMIP(:)
847	!$OMP THREADPRIVATE(tabCFMIP, lonCFMIP, latCFMIP)
848	INTEGER, ALLOCATABLE, SAVE :: tabijGCM(:)
849	REAL, ALLOCATABLE, SAVE :: lonGCM(:), latGCM(:)
850	!$OMP THREADPRIVATE(tabijGCM, lonGCM, latGCM)
851	INTEGER, ALLOCATABLE, SAVE :: iGCM(:), jGCM(:)
852	!$OMP THREADPRIVATE(iGCM, jGCM)
853	logical, dimension(nfiles) :: phys_out_filestations
854	logical, parameter :: lNMC=.FALSE.
855
856	!IM betaCRF
857	REAL, SAVE :: pfree, beta_pbl, beta_free
858	!$OMP THREADPRIVATE(pfree, beta_pbl, beta_free)
859	REAL, SAVE :: lon1_beta, lon2_beta, lat1_beta, lat2_beta
860	!$OMP THREADPRIVATE(lon1_beta, lon2_beta, lat1_beta, lat2_beta)
861	LOGICAL, SAVE :: mskocean_beta
862	!$OMP THREADPRIVATE(mskocean_beta)
863	REAL, dimension(klon, klev) :: beta ! facteur sur cldtaurad et cldemirad pour evaluer les retros liees aux CRF
864	REAL, dimension(klon, klev) :: cldtaurad ! epaisseur optique pour radlwsw pour tester "CRF off"
865	REAL, dimension(klon, klev) :: cldtaupirad ! epaisseur optique pour radlwsw pour tester "CRF off"
866	REAL, dimension(klon, klev) :: cldemirad ! emissivite pour radlwsw pour tester "CRF off"
867	REAL, dimension(klon, klev) :: cldfrarad ! fraction nuageuse
868
869	INTEGER :: nbtr_tmp ! Number of tracer inside concvl
870	REAL, dimension(klon,klev) :: sh_in ! Specific humidity entering in phytrac
871	integer iostat
872
873	REAL zzz
874
875	!======================================================================
876	! Gestion calendrier : mise a jour du module phys_cal_mod
877	!
878	CALL phys_cal_update(jD_cur,jH_cur)
879
880	!======================================================================
881	! Ecriture eventuelle d'un profil verticale en entree de la physique.
882	! Utilise notamment en 1D mais peut etre active egalement en 3D
883	! en imposant la valeur de igout.
884	!======================================================================d
885	if (prt_level.ge.1) then
886	igout=klon/2+1/klon
887	write(lunout,*) 'DEBUT DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
888	write(lunout,*) &
889	'nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys'
890	write(lunout,*) &
891	nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys
892
893	write(lunout,*) 'paprs, play, phi, u, v, t'
894	do k=1,klev
895	write(lunout,*) paprs(igout,k),pplay(igout,k),pphi(igout,k), &
896	u(igout,k),v(igout,k),t(igout,k)
897	enddo
898	write(lunout,*) 'ovap (g/kg), oliq (g/kg)'
899	do k=1,klev
900	write(lunout,) qx(igout,k,1)1000,qx(igout,k,2)*1000.
901	enddo
902	endif
903
904	!======================================================================
905
906	if (first) then
907
908	!CR:nvelles variables convection/poches froides
909
910	print*, '================================================='
911	print*, 'Allocation des variables locales et sauvegardees'
912	call phys_local_var_init
913	!
914	pasphys=pdtphys
915	! appel a la lecture du run.def physique
916	call conf_phys(ok_journe, ok_mensuel, &
917	ok_instan, ok_hf, &
918	ok_LES, &
919	callstats, &
920	solarlong0,seuil_inversion, &
921	fact_cldcon, facttemps,ok_newmicro,iflag_radia, &
922	iflag_cldcon,iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, &
923	ok_ade, ok_aie, ok_cdnc, aerosol_couple, &
924	flag_aerosol, flag_aerosol_strat, new_aod, &
925	bl95_b0, bl95_b1, &
926	! nv flags pour la convection et les poches froides
927	read_climoz, &
928	alp_offset)
929	call phys_state_var_init(read_climoz)
930	call phys_output_var_init
931	print*, '================================================='
932	!
933	dnwd0=0.0
934	ftd=0.0
935	fqd=0.0
936	cin=0.
937	!ym Attention pbase pas initialise dans concvl !!!!
938	pbase=0
939	!IM 180608
940
941	itau_con=0
942	first=.false.
943
944	endif ! first
945
946	!ym => necessaire pour iflag_con != 2
947	pmfd(:,:) = 0.
948	pen_u(:,:) = 0.
949	pen_d(:,:) = 0.
950	pde_d(:,:) = 0.
951	pde_u(:,:) = 0.
952	aam=0.
953
954	torsfc=0.
955	forall (k=1: llm) zmasse(:, k) = (paprs(:, k)-paprs(:, k+1)) / rg
956
957
958
959	modname = 'physiq'
960	!IM
961	IF (ip_ebil_phy.ge.1) THEN
962	DO i=1,klon
963	zero_v(i)=0.
964	END DO
965	END IF
966	ok_sync=.TRUE.
967
968	IF (debut) THEN
969	CALL suphel ! initialiser constantes et parametres phys.
970	ENDIF
971
972	if(prt_level.ge.1) print*,'CONVERGENCE PHYSIQUE THERM 1 '
973
974
975	!======================================================================
976	! Gestion calendrier : mise a jour du module phys_cal_mod
977	!
978	! CALL phys_cal_update(jD_cur,jH_cur)
979
980	!
981	! Si c'est le debut, il faut initialiser plusieurs choses
982	! ********
983	!
984	IF (debut) THEN
985	!rv
986	!CRinitialisation de wght_th et lalim_conv pour la definition de la couche alimentation
987	!de la convection a partir des caracteristiques du thermique
988	wght_th(:,:)=1.
989	lalim_conv(:)=1
990	!RC
991	ustar(:,:)=0.
992	u10m(:,:)=0.
993	v10m(:,:)=0.
994	rain_con(:)=0.
995	snow_con(:)=0.
996	topswai(:)=0.
997	topswad(:)=0.
998	solswai(:)=0.
999	solswad(:)=0.
1000
1001	wmax_th(:)=0.
1002	tau_overturning_th(:)=0.
1003
1004	IF (type_trac == 'inca') THEN
1005	! jg : initialisation jusqu'au ces variables sont dans restart
1006	ccm(:,:,:) = 0.
1007	tau_aero(:,:,:,:) = 0.
1008	piz_aero(:,:,:,:) = 0.
1009	cg_aero(:,:,:,:) = 0.
1010	END IF
1011
1012	rnebcon0(:,:) = 0.0
1013	clwcon0(:,:) = 0.0
1014	rnebcon(:,:) = 0.0
1015	clwcon(:,:) = 0.0
1016
1017	!IM
1018	IF (ip_ebil_phy.ge.1) d_h_vcol_phy=0.
1019	!
1020	print*,'iflag_coupl,iflag_clos,iflag_wake', &
1021	iflag_coupl,iflag_clos,iflag_wake
1022	print*,'CYCLE_DIURNE', cycle_diurne
1023	!
1024	IF (iflag_con.EQ.2.AND.iflag_cldcon.GT.-1) THEN
1025	abort_message = 'Tiedtke needs iflag_cldcon=-2 or -1'
1026	CALL abort_gcm (modname,abort_message,1)
1027	ENDIF
1028	!
1029	!
1030	! Initialiser les compteurs:
1031	!
1032	itap = 0
1033	itaprad = 0
1034
1035	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1036	!! Un petit travail \`a faire ici.
1037	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1038
1039	if (iflag_pbl>1) then
1040	PRINT*, "Using method MELLOR&YAMADA"
1041	endif
1042
1043	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1044	! FH 2008/05/02 changement lie a la lecture de nbapp_rad dans phylmd plutot que
1045	! dyn3d
1046	! Attention : la version precedente n'etait pas tres propre.
1047	! Il se peut qu'il faille prendre une valeur differente de nbapp_rad
1048	! pour obtenir le meme resultat.
1049	dtime=pdtphys
1050	radpas = NINT( 86400./dtime/nbapp_rad)
1051	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1052
1053	CALL phyetat0 ("startphy.nc",clesphy0,tabcntr0)
1054	IF (klon_glo==1) THEN
1055	coefh=0. ; coefm=0. ; pbl_tke=0.
1056	coefh(:,2,:)=1.e-2 ; coefm(:,2,:)=1.e-2 ; pbl_tke(:,2,:)=1.e-2
1057	PRINT*,'FH WARNING : lignes a supprimer'
1058	ENDIF
1059	!IM begin
1060	print*,'physiq: clwcon rnebcon ratqs',clwcon(1,1),rnebcon(1,1) &
1061	,ratqs(1,1)
1062	!IM end
1063
1064
1065
1066	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1067	!
1068	! on remet le calendrier a zero
1069	!
1070	IF (raz_date .eq. 1) THEN
1071	itau_phy = 0
1072	ENDIF
1073
1074	!IM cf. AM 081204 BEG
1075	PRINT*,'cycle_diurne3 =',cycle_diurne
1076	!IM cf. AM 081204 END
1077	!
1078	CALL printflag( tabcntr0,radpas,ok_journe, &
1079	ok_instan, ok_region )
1080	!
1081	IF (ABS(dtime-pdtphys).GT.0.001) THEN
1082	WRITE(lunout,*) 'Pas physique n est pas correct',dtime, &
1083	pdtphys
1084	abort_message='Pas physique n est pas correct '
1085	! call abort_gcm(modname,abort_message,1)
1086	dtime=pdtphys
1087	ENDIF
1088	IF (nlon .NE. klon) THEN
1089	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon, &
1090	klon
1091	abort_message='nlon et klon ne sont pas coherents'
1092	call abort_gcm(modname,abort_message,1)
1093	ENDIF
1094	IF (nlev .NE. klev) THEN
1095	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev, &
1096	klev
1097	abort_message='nlev et klev ne sont pas coherents'
1098	call abort_gcm(modname,abort_message,1)
1099	ENDIF
1100	!
1101	IF (dtime*REAL(radpas).GT.21600..AND.cycle_diurne) THEN
1102	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
1103	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
1104	abort_message='Nbre d appels au rayonnement insuffisant'
1105	call abort_gcm(modname,abort_message,1)
1106	ENDIF
1107	WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
1108	WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=", &
1109	ok_cvl
1110	!
1111	!KE43
1112	! Initialisation pour la convection de K.E. (sb):
1113	IF (iflag_con.GE.3) THEN
1114
1115	WRITE(lunout,)"** Convection de Kerry Emanuel 4.3 "
1116	WRITE(lunout,*) &
1117	"On va utiliser le melange convectif des traceurs qui"
1118	WRITE(lunout,*)"est calcule dans convect4.3"
1119	WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
1120
1121	DO i = 1, klon
1122	ema_cbmf(i) = 0.
1123	ema_pcb(i) = 0.
1124	ema_pct(i) = 0.
1125	! ema_workcbmf(i) = 0.
1126	ENDDO
1127	!IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
1128	DO i = 1, klon
1129	ibas_con(i) = 1
1130	itop_con(i) = 1
1131	ENDDO
1132	!IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
1133	!===============================================================================
1134	!CR:04.12.07: initialisations poches froides
1135	! Controle de ALE et ALP pour la fermeture convective (jyg)
1136	if (iflag_wake>=1) then
1137	CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr &
1138	,alp_bl_prescr, ale_bl_prescr)
1139	! 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
1140	! print*,'apres ini_wake iflag_cldcon=', iflag_cldcon
1141	endif
1142
1143	do i = 1,klon
1144	Ale_bl(i)=0.
1145	Alp_bl(i)=0.
1146	enddo
1147
1148	!================================================================================
1149	!IM stations CFMIP
1150	nCFMIP=npCFMIP
1151	OPEN(98,file='npCFMIP_param.data',status='old', &
1152	form='formatted',iostat=iostat)
1153	if (iostat == 0) then
1154	READ(98,*,end=998) nCFMIP
1155	998 CONTINUE
1156	CLOSE(98)
1157	CONTINUE
1158	IF(nCFMIP.GT.npCFMIP) THEN
1159	print*,'nCFMIP > npCFMIP : augmenter npCFMIP et recompiler'
1160	call abort_gcm("physiq", "", 1)
1161	else
1162	print*,'physiq npCFMIP=',npCFMIP,'nCFMIP=',nCFMIP
1163	ENDIF
1164
1165	!
1166	ALLOCATE(tabCFMIP(nCFMIP))
1167	ALLOCATE(lonCFMIP(nCFMIP), latCFMIP(nCFMIP))
1168	ALLOCATE(tabijGCM(nCFMIP))
1169	ALLOCATE(lonGCM(nCFMIP), latGCM(nCFMIP))
1170	ALLOCATE(iGCM(nCFMIP), jGCM(nCFMIP))
1171	!
1172	! lecture des nCFMIP stations CFMIP, de leur numero
1173	! et des coordonnees geographiques lonCFMIP, latCFMIP
1174	!
1175	CALL read_CFMIP_point_locations(nCFMIP, tabCFMIP, &
1176	lonCFMIP, latCFMIP)
1177	!
1178	! identification des
1179	! 1) coordonnees lonGCM, latGCM des points CFMIP dans la grille de LMDZ
1180	! 2) indices points tabijGCM de la grille physique 1d sur klon points
1181	! 3) indices iGCM, jGCM de la grille physique 2d
1182	!
1183	CALL LMDZ_CFMIP_point_locations(nCFMIP, lonCFMIP, latCFMIP, &
1184	tabijGCM, lonGCM, latGCM, iGCM, jGCM)
1185	!
1186	else
1187	ALLOCATE(tabijGCM(0))
1188	ALLOCATE(lonGCM(0), latGCM(0))
1189	ALLOCATE(iGCM(0), jGCM(0))
1190	end if
1191	else
1192	ALLOCATE(tabijGCM(0))
1193	ALLOCATE(lonGCM(0), latGCM(0))
1194	ALLOCATE(iGCM(0), jGCM(0))
1195	ENDIF
1196
1197	DO i=1,klon
1198	rugoro(i) = f_rugoro * MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1199	ENDDO
1200
1201	!34EK
1202	IF (ok_orodr) THEN
1203
1204	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1205	! FH sans doute a enlever de finitivement ou, si on le garde, l'activer
1206	! justement quand ok_orodr = false.
1207	! ce rugoro est utilise par la couche limite et fait double emploi
1208	! avec les param\'etrisations sp\'ecifiques de Francois Lott.
1209	! DO i=1,klon
1210	! rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1211	! ENDDO
1212	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1213	IF (ok_strato) THEN
1214	CALL SUGWD_strato(klon,klev,paprs,pplay)
1215	ELSE
1216	CALL SUGWD(klon,klev,paprs,pplay)
1217	ENDIF
1218
1219	DO i=1,klon
1220	zuthe(i)=0.
1221	zvthe(i)=0.
1222	if(zstd(i).gt.10.)then
1223	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
1224	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
1225	endif
1226	ENDDO
1227	ENDIF
1228	!
1229	!
1230	lmt_pas = NINT(86400./dtime * 1.0) ! tous les jours
1231	WRITE(lunout,*)'La frequence de lecture surface est de ', &
1232	lmt_pas
1233	!
1234	capemaxcels = 't_max(X)'
1235	t2mincels = 't_min(X)'
1236	t2maxcels = 't_max(X)'
1237	tinst = 'inst(X)'
1238	tave = 'ave(X)'
1239	!IM cf. AM 081204 BEG
1240	write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
1241	!IM cf. AM 081204 END
1242	!
1243	!=============================================================
1244	! Initialisation des sorties
1245	!=============================================================
1246
1247	#ifdef CPP_IOIPSL
1248
1249	!$OMP MASTER
1250	call phys_output_open(rlon,rlat,nCFMIP,tabijGCM, &
1251	iGCM,jGCM,lonGCM,latGCM, &
1252	jjmp1,nlevSTD,clevSTD,rlevSTD, &
1253	nbteta, ctetaSTD, dtime,ok_veget, &
1254	type_ocean,iflag_pbl,ok_mensuel,ok_journe, &
1255	ok_hf,ok_instan,ok_LES,ok_ade,ok_aie, &
1256	read_climoz, phys_out_filestations, &
1257	new_aod, aerosol_couple, &
1258	flag_aerosol_strat, pdtphys, paprs, pphis, &
1259	pplay, lmax_th, ptconv, ptconvth, ivap, &
1260	d_t, qx, d_qx, zmasse, ok_sync)
1261	!$OMP END MASTER
1262	!$OMP BARRIER
1263
1264	freq_outNMC(1) = ecrit_files(7)
1265	freq_outNMC(2) = ecrit_files(8)
1266	freq_outNMC(3) = ecrit_files(9)
1267	WRITE(lunout,*)'OK freq_outNMC(1)=',freq_outNMC(1)
1268	WRITE(lunout,*)'OK freq_outNMC(2)=',freq_outNMC(2)
1269	WRITE(lunout,*)'OK freq_outNMC(3)=',freq_outNMC(3)
1270
1271	include "ini_histday_seri.h"
1272
1273	include "ini_paramLMDZ_phy.h"
1274
1275	#endif
1276	ecrit_reg = ecrit_reg * un_jour
1277	ecrit_tra = ecrit_tra * un_jour
1278
1279	!XXXPB Positionner date0 pour initialisation de ORCHIDEE
1280	date0 = jD_ref
1281	WRITE(,) 'physiq date0 : ',date0
1282	!
1283	!
1284	!
1285	! Prescrire l'ozone dans l'atmosphere
1286	!
1287	!
1288	!c DO i = 1, klon
1289	!c DO k = 1, klev
1290	!c CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1291	!c ENDDO
1292	!c ENDDO
1293	!
1294	IF (type_trac == 'inca') THEN
1295	#ifdef INCA
1296	CALL VTe(VTphysiq)
1297	CALL VTb(VTinca)
1298	calday = REAL(days_elapsed) + jH_cur
1299	WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
1300
1301	CALL chemini( &
1302	rg, &
1303	ra, &
1304	airephy, &
1305	rlat, &
1306	rlon, &
1307	presnivs, &
1308	calday, &
1309	klon, &
1310	nqtot, &
1311	pdtphys, &
1312	annee_ref, &
1313	day_ref, &
1314	itau_phy)
1315
1316	CALL VTe(VTinca)
1317	CALL VTb(VTphysiq)
1318	#endif
1319	END IF
1320	!
1321	!
1322	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1323	! Nouvelle initialisation pour le rayonnement RRTM
1324	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1325
1326	call iniradia(klon,klev,paprs(1,1:klev+1))
1327
1328	!$omp single
1329	if (read_climoz >= 1) then
1330	call open_climoz(ncid_climoz, press_climoz)
1331	END IF
1332	!$omp end single
1333	!
1334	!IM betaCRF
1335	pfree=70000. !Pa
1336	beta_pbl=1.
1337	beta_free=1.
1338	lon1_beta=-180.
1339	lon2_beta=+180.
1340	lat1_beta=90.
1341	lat2_beta=-90.
1342	mskocean_beta=.FALSE.
1343
1344	OPEN(99,file='beta_crf.data',status='old', &
1345	form='formatted',err=9999)
1346	READ(99,*,end=9998) pfree
1347	READ(99,*,end=9998) beta_pbl
1348	READ(99,*,end=9998) beta_free
1349	READ(99,*,end=9998) lon1_beta
1350	READ(99,*,end=9998) lon2_beta
1351	READ(99,*,end=9998) lat1_beta
1352	READ(99,*,end=9998) lat2_beta
1353	READ(99,*,end=9998) mskocean_beta
1354	9998 Continue
1355	CLOSE(99)
1356	9999 Continue
1357	WRITE(,)'pfree=',pfree
1358	WRITE(,)'beta_pbl=',beta_pbl
1359	WRITE(,)'beta_free=',beta_free
1360	WRITE(,)'lon1_beta=',lon1_beta
1361	WRITE(,)'lon2_beta=',lon2_beta
1362	WRITE(,)'lat1_beta=',lat1_beta
1363	WRITE(,)'lat2_beta=',lat2_beta
1364	WRITE(,)'mskocean_beta=',mskocean_beta
1365	ENDIF
1366	!
1367	! ************** Fin de IF ( debut ) *************
1368	!
1369	!
1370	! Incrementer le compteur de la physique
1371	!
1372	itap = itap + 1
1373	!
1374	!
1375	! Update fraction of the sub-surfaces (pctsrf) and
1376	! initialize, where a new fraction has appeared, all variables depending
1377	! on the surface fraction.
1378	!
1379	CALL change_srf_frac(itap, dtime, days_elapsed+1, &
1380	pctsrf, falb1, falb2, ftsol, ustar, u10m, v10m, pbl_tke)
1381
1382
1383	! Update time and other variables in Reprobus
1384	IF (type_trac == 'repr') THEN
1385	#ifdef REPROBUS
1386	CALL Init_chem_rep_xjour(jD_cur-jD_ref+day_ref)
1387	print*,'xjour equivalent rjourvrai',jD_cur-jD_ref+day_ref
1388	CALL Rtime(debut)
1389	#endif
1390	END IF
1391
1392
1393	! Tendances bidons pour les processus qui n'affectent pas certaines
1394	! variables.
1395	du0(:,:)=0.
1396	dv0(:,:)=0.
1397	dt0 = 0.
1398	dq0(:,:)=0.
1399	dql0(:,:)=0.
1400	!
1401	! Mettre a zero des variables de sortie (pour securite)
1402	!
1403	DO i = 1, klon
1404	d_ps(i) = 0.0
1405	ENDDO
1406	DO k = 1, klev
1407	DO i = 1, klon
1408	d_t(i,k) = 0.0
1409	d_u(i,k) = 0.0
1410	d_v(i,k) = 0.0
1411	ENDDO
1412	ENDDO
1413	DO iq = 1, nqtot
1414	DO k = 1, klev
1415	DO i = 1, klon
1416	d_qx(i,k,iq) = 0.0
1417	ENDDO
1418	ENDDO
1419	ENDDO
1420	da(:,:)=0.
1421	mp(:,:)=0.
1422	phi(:,:,:)=0.
1423	! RomP >>>
1424	phi2(:,:,:)=0.
1425	beta_prec_fisrt(:,:)=0.
1426	beta_prec(:,:)=0.
1427	epmlmMm(:,:,:)=0.
1428	eplaMm(:,:)=0.
1429	d1a(:,:)=0.
1430	dam(:,:)=0.
1431	pmflxr=0.
1432	pmflxs=0.
1433	! RomP <<<
1434
1435	!
1436	! Ne pas affecter les valeurs entrees de u, v, h, et q
1437	!
1438	DO k = 1, klev
1439	DO i = 1, klon
1440	t_seri(i,k) = t(i,k)
1441	u_seri(i,k) = u(i,k)
1442	v_seri(i,k) = v(i,k)
1443	q_seri(i,k) = qx(i,k,ivap)
1444	ql_seri(i,k) = qx(i,k,iliq)
1445	qs_seri(i,k) = 0.
1446	ENDDO
1447	ENDDO
1448	tke0(:,:)=pbl_tke(:,:,is_ave)
1449	IF (nqtot.GE.3) THEN
1450	DO iq = 3, nqtot
1451	DO k = 1, klev
1452	DO i = 1, klon
1453	tr_seri(i,k,iq-2) = qx(i,k,iq)
1454	ENDDO
1455	ENDDO
1456	ENDDO
1457	ELSE
1458	DO k = 1, klev
1459	DO i = 1, klon
1460	tr_seri(i,k,1) = 0.0
1461	ENDDO
1462	ENDDO
1463	ENDIF
1464	!
1465	DO i = 1, klon
1466	ztsol(i) = 0.
1467	ENDDO
1468	DO nsrf = 1, nbsrf
1469	DO i = 1, klon
1470	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1471	ENDDO
1472	ENDDO
1473	!IM
1474	IF (ip_ebil_phy.ge.1) THEN
1475	ztit='after dynamic'
1476	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime &
1477	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
1478	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1479	! Comme les tendances de la physique sont ajoute dans la dynamique,
1480	! on devrait avoir que la variation d'entalpie par la dynamique
1481	! est egale a la variation de la physique au pas de temps precedent.
1482	! Donc la somme de ces 2 variations devrait etre nulle.
1483	call diagphy(airephy,ztit,ip_ebil_phy &
1484	, zero_v, zero_v, zero_v, zero_v, zero_v &
1485	, zero_v, zero_v, zero_v, ztsol &
1486	, d_h_vcol+d_h_vcol_phy, d_qt, 0. &
1487	, fs_bound, fq_bound )
1488	END IF
1489
1490	! Diagnostiquer la tendance dynamique
1491	!
1492	IF (ancien_ok) THEN
1493	DO k = 1, klev
1494	DO i = 1, klon
1495	d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
1496	d_v_dyn(i,k) = (v_seri(i,k)-v_ancien(i,k))/dtime
1497	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
1498	d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
1499	ENDDO
1500	ENDDO
1501	!!! RomP >>> td dyn traceur
1502	IF (nqtot.GE.3) THEN
1503	DO iq = 3, nqtot
1504	DO k = 1, klev
1505	DO i = 1, klon
1506	d_tr_dyn(i,k,iq-2)= &
1507	(tr_seri(i,k,iq-2)-tr_ancien(i,k,iq-2))/dtime
1508	! iiq=niadv(iq)
1509	! print*,i,k," d_tr_dyn",d_tr_dyn(i,k,iq-2),"tra:",iq,tname(iiq)
1510	ENDDO
1511	ENDDO
1512	ENDDO
1513	ENDIF
1514	!!! RomP <<<
1515	ELSE
1516	DO k = 1, klev
1517	DO i = 1, klon
1518	d_u_dyn(i,k) = 0.0
1519	d_v_dyn(i,k) = 0.0
1520	d_t_dyn(i,k) = 0.0
1521	d_q_dyn(i,k) = 0.0
1522	ENDDO
1523	ENDDO
1524	!!! RomP >>> td dyn traceur
1525	IF (nqtot.GE.3) THEN
1526	DO iq = 3, nqtot
1527	DO k = 1, klev
1528	DO i = 1, klon
1529	d_tr_dyn(i,k,iq-2)= 0.0
1530	ENDDO
1531	ENDDO
1532	ENDDO
1533	ENDIF
1534	!!! RomP <<<
1535	ancien_ok = .TRUE.
1536	ENDIF
1537	!
1538	! Ajouter le geopotentiel du sol:
1539	!
1540	DO k = 1, klev
1541	DO i = 1, klon
1542	zphi(i,k) = pphi(i,k) + pphis(i)
1543	ENDDO
1544	ENDDO
1545	!
1546	! Verifier les temperatures
1547	!
1548	!IM BEG
1549	IF (check) THEN
1550	amn=MIN(ftsol(1,is_ter),1000.)
1551	amx=MAX(ftsol(1,is_ter),-1000.)
1552	DO i=2, klon
1553	amn=MIN(ftsol(i,is_ter),amn)
1554	amx=MAX(ftsol(i,is_ter),amx)
1555	ENDDO
1556	!
1557	PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
1558	ENDIF !(check) THEN
1559	!IM END
1560	!
1561	CALL hgardfou(t_seri,ftsol,'debutphy')
1562	!
1563	!IM BEG
1564	IF (check) THEN
1565	amn=MIN(ftsol(1,is_ter),1000.)
1566	amx=MAX(ftsol(1,is_ter),-1000.)
1567	DO i=2, klon
1568	amn=MIN(ftsol(i,is_ter),amn)
1569	amx=MAX(ftsol(i,is_ter),amx)
1570	ENDDO
1571	!
1572	PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
1573	ENDIF !(check) THEN
1574	!IM END
1575	!
1576	! Mettre en action les conditions aux limites (albedo, sst, etc.).
1577	! Prescrire l'ozone et calculer l'albedo sur l'ocean.
1578	!
1579	if (read_climoz >= 1) then
1580	! Ozone from a file
1581	! Update required ozone index:
1582	ro3i = int((days_elapsed + jh_cur - jh_1jan) / ioget_year_len(year_cur) &
1583	* 360.) + 1
1584	if (ro3i == 361) ro3i = 360
1585	! (This should never occur, except perhaps because of roundup
1586	! error. See documentation.)
1587	if (ro3i /= co3i) then
1588	! Update ozone field:
1589	if (read_climoz == 1) then
1590	call regr_pr_av(ncid_climoz, (/"tro3"/), julien=ro3i, &
1591	press_in_edg=press_climoz, paprs=paprs, v3=wo)
1592	else
1593	! read_climoz == 2
1594	call regr_pr_av(ncid_climoz, (/"tro3 ", "tro3_daylight"/), &
1595	julien=ro3i, press_in_edg=press_climoz, paprs=paprs, v3=wo)
1596	end if
1597	! Convert from mole fraction of ozone to column density of ozone in a
1598	! cell, in kDU:
1599	forall (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) * rmo3 / rmd &
1600	* zmasse / dobson_u / 1e3
1601	! (By regridding ozone values for LMDZ only once every 360th of
1602	! year, we have already neglected the variation of pressure in one
1603	! 360th of year. So do not recompute "wo" at each time step even if
1604	! "zmasse" changes a little.)
1605	co3i = ro3i
1606	end if
1607	ELSEIF (MOD(itap-1,lmt_pas) == 0) THEN
1608	! Once per day, update ozone from Royer:
1609
1610	IF (solarlong0<-999.) then
1611	! Generic case with evolvoing season
1612	zzz=real(days_elapsed+1)
1613	ELSE IF (abs(solarlong0-1000.)<1.e-4) then
1614	! Particular case with annual mean insolation
1615	zzz=real(90) ! could be revisited
1616	IF (read_climoz/=-1) THEN
1617	abort_message ='read_climoz=-1 is recommended when solarlong0=1000.'
1618	CALL abort_gcm (modname,abort_message,1)
1619	ENDIF
1620	ELSE
1621	! Case where the season is imposed with solarlong0
1622	zzz=real(90) ! could be revisited
1623	ENDIF
1624	wo(:,:,1)=ozonecm(rlat, paprs,read_climoz,rjour=zzz)
1625	ENDIF
1626	!
1627	! Re-evaporer l'eau liquide nuageuse
1628	!
1629	DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse
1630	DO i = 1, klon
1631	zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1632	!jyg<
1633	! Attention : Arnaud a propose des formules completement differentes
1634	! A verifier !!!
1635	zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1636	IF (iflag_ice_thermo .EQ. 0) THEN
1637	zlsdcp=zlvdcp
1638	ENDIF
1639	!>jyg
1640
1641	zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
1642	zb = MAX(0.0,ql_seri(i,k))
1643	za = - MAX(0.0,ql_seri(i,k)) &
1644	* (zlvdcp(1.-zdelta)+zlsdcpzdelta)
1645	t_seri(i,k) = t_seri(i,k) + za
1646	q_seri(i,k) = q_seri(i,k) + zb
1647	ql_seri(i,k) = 0.0
1648	d_t_eva(i,k) = za
1649	d_q_eva(i,k) = zb
1650	ENDDO
1651	ENDDO
1652	!IM
1653	IF (ip_ebil_phy.ge.2) THEN
1654	ztit='after reevap'
1655	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,1,dtime &
1656	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
1657	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1658	call diagphy(airephy,ztit,ip_ebil_phy &
1659	, zero_v, zero_v, zero_v, zero_v, zero_v &
1660	, zero_v, zero_v, zero_v, ztsol &
1661	, d_h_vcol, d_qt, d_ec &
1662	, fs_bound, fq_bound )
1663	!
1664	END IF
1665
1666	!
1667	!=========================================================================
1668	! Calculs de l'orbite.
1669	! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
1670	! doit donc etre plac\'e avant radlwsw et pbl_surface
1671
1672	!!! jyg 17 Sep 2010 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1673	call ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
1674	day_since_equinox = (jD_cur + jH_cur) - jD_eq
1675	!
1676	! choix entre calcul de la longitude solaire vraie ou valeur fixee a
1677	! solarlong0
1678	if (solarlong0<-999.) then
1679	if (new_orbit) then
1680	! calcul selon la routine utilisee pour les planetes
1681	call solarlong(day_since_equinox, zlongi, dist)
1682	else
1683	! calcul selon la routine utilisee pour l'AR4
1684	CALL orbite(REAL(days_elapsed+1),zlongi,dist)
1685	endif
1686	else
1687	zlongi=solarlong0 ! longitude solaire vraie
1688	dist=1. ! distance au soleil / moyenne
1689	endif
1690	if(prt_level.ge.1) &
1691	write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
1692
1693
1694	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1695	! Calcul de l'ensoleillement :
1696	! ============================
1697	! Pour une solarlong0=1000., on calcule un ensoleillement moyen sur
1698	! l'annee a partir d'une formule analytique.
1699	! Cet ensoleillement est sym\'etrique autour de l'\'equateur et
1700	! non nul aux poles.
1701	IF (abs(solarlong0-1000.)<1.e-4) then
1702	call zenang_an(cycle_diurne,jH_cur,rlat,rlon,rmu0,fract)
1703	ELSE
1704	! Avec ou sans cycle diurne
1705	IF (cycle_diurne) THEN
1706	zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s)
1707	CALL zenang(zlongi,jH_cur,zdtime,rlat,rlon,rmu0,fract)
1708	ELSE
1709	CALL angle(zlongi, rlat, fract, rmu0)
1710	ENDIF
1711	ENDIF
1712
1713	! AI Janv 2014
1714	do i = 1, klon
1715	if (fract(i).le.0.) then
1716	JrNt(i)=0.
1717	else
1718	JrNt(i)=1.
1719	endif
1720	enddo
1721
1722	if (mydebug) then
1723	call writefield_phy('u_seri',u_seri,llm)
1724	call writefield_phy('v_seri',v_seri,llm)
1725	call writefield_phy('t_seri',t_seri,llm)
1726	call writefield_phy('q_seri',q_seri,llm)
1727	endif
1728
1729	!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1730	! Appel au pbl_surface : Planetary Boudary Layer et Surface
1731	! Cela implique tous les interactions des sous-surfaces et la partie diffusion
1732	! turbulent du couche limit.
1733	!
1734	! Certains varibales de sorties de pbl_surface sont utiliser que pour
1735	! ecriture des fihiers hist_XXXX.nc, ces sont :
1736	! qsol, zq2m, s_pblh, s_lcl,
1737	! s_capCL, s_oliqCL, s_cteiCL,s_pblT,
1738	! s_therm, s_trmb1, s_trmb2, s_trmb3,
1739	! zxrugs, zu10m, zv10m, fder,
1740	! zxqsurf, rh2m, zxfluxu, zxfluxv,
1741	! frugs, agesno, fsollw, fsolsw,
1742	! d_ts, fevap, fluxlat, t2m,
1743	! wfbils, wfbilo, fluxt, fluxu, fluxv,
1744	!
1745	! Certains ne sont pas utiliser du tout :
1746	! dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
1747	!
1748
1749	! Calcul de l'humidite de saturation au niveau du sol
1750
1751
1752
1753	if (iflag_pbl/=0) then
1754
1755	CALL pbl_surface( &
1756	dtime, date0, itap, days_elapsed+1, &
1757	debut, lafin, &
1758	rlon, rlat, rugoro, rmu0, &
1759	zsig, sollwdown, pphi, cldt, &
1760	rain_fall, snow_fall, solsw, sollw, &
1761	t_seri, q_seri, u_seri, v_seri, &
1762	pplay, paprs, pctsrf, &
1763	ftsol,falb1,falb2,ustar,u10m,v10m,wstar, &
1764	sollwdown, cdragh, cdragm, u1, v1, &
1765	albsol1, albsol2, sens, evap, &
1766	albsol3_lic,runoff, snowhgt, qsnow, to_ice, sissnow, &
1767	zxtsol, zxfluxlat, zt2m, qsat2m, &
1768	d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf, d_t_diss, &
1769	coefh(1:klon,1:klev,1:nbsrf+1), coefm(1:klon,1:klev,1:nbsrf+1), &
1770	slab_wfbils, &
1771	qsol, zq2m, s_pblh, s_lcl, &
1772	s_capCL, s_oliqCL, s_cteiCL,s_pblT, &
1773	s_therm, s_trmb1, s_trmb2, s_trmb3, &
1774	zxrugs, zustar, zu10m, zv10m, fder, &
1775	zxqsurf, rh2m, zxfluxu, zxfluxv, &
1776	frugs, agesno, fsollw, fsolsw, &
1777	d_ts, fevap, fluxlat, t2m, &
1778	wfbils, wfbilo, fluxt, fluxu, fluxv, &
1779	dsens, devap, zxsnow, &
1780	zxfluxt, zxfluxq, q2m, fluxq, pbl_tke )
1781
1782
1783	!---------------------------------------------------------------------
1784	! ajout des tendances de la diffusion turbulente
1785	IF (klon_glo==1) THEN
1786	CALL add_pbl_tend &
1787	(d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,paprs,'vdf')
1788	ELSE
1789	CALL add_phys_tend &
1790	(d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,paprs,'vdf')
1791	ENDIF
1792	!--------------------------------------------------------------------
1793
1794	if (mydebug) then
1795	call writefield_phy('u_seri',u_seri,llm)
1796	call writefield_phy('v_seri',v_seri,llm)
1797	call writefield_phy('t_seri',t_seri,llm)
1798	call writefield_phy('q_seri',q_seri,llm)
1799	endif
1800
1801	CALL evappot(klon,nbsrf,ftsol,pplay(:,1),cdragh, &
1802	t_seri(:,1),q_seri(:,1),u_seri(:,1),v_seri(:,1),evap_pot)
1803
1804
1805	IF (ip_ebil_phy.ge.2) THEN
1806	ztit='after surface_main'
1807	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
1808	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
1809	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1810	call diagphy(airephy,ztit,ip_ebil_phy &
1811	, zero_v, zero_v, zero_v, zero_v, sens &
1812	, evap , zero_v, zero_v, ztsol &
1813	, d_h_vcol, d_qt, d_ec &
1814	, fs_bound, fq_bound )
1815	END IF
1816
1817	ENDIF
1818	! =================================================================== c
1819	! Calcul de Qsat
1820
1821	DO k = 1, klev
1822	DO i = 1, klon
1823	zx_t = t_seri(i,k)
1824	IF (thermcep) THEN
1825	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
1826	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
1827	zx_qs = MIN(0.5,zx_qs)
1828	zcor = 1./(1.-retv*zx_qs)
1829	zx_qs = zx_qs*zcor
1830	ELSE
1831	IF (zx_t.LT.t_coup) THEN
1832	zx_qs = qsats(zx_t)/pplay(i,k)
1833	ELSE
1834	zx_qs = qsatl(zx_t)/pplay(i,k)
1835	ENDIF
1836	ENDIF
1837	zqsat(i,k)=zx_qs
1838	ENDDO
1839	ENDDO
1840
1841	if (prt_level.ge.1) then
1842	write(lunout,*) 'L qsat (g/kg) avant clouds_gno'
1843	write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
1844	endif
1845	!
1846	! Appeler la convection (au choix)
1847	!
1848	DO k = 1, klev
1849	DO i = 1, klon
1850	conv_q(i,k) = d_q_dyn(i,k) &
1851	+ d_q_vdf(i,k)/dtime
1852	conv_t(i,k) = d_t_dyn(i,k) &
1853	+ d_t_vdf(i,k)/dtime
1854	ENDDO
1855	ENDDO
1856	IF (check) THEN
1857	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
1858	WRITE(lunout,*) "avantcon=", za
1859	ENDIF
1860	zx_ajustq = .FALSE.
1861	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
1862	IF (zx_ajustq) THEN
1863	DO i = 1, klon
1864	z_avant(i) = 0.0
1865	ENDDO
1866	DO k = 1, klev
1867	DO i = 1, klon
1868	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k)) &
1869	*(paprs(i,k)-paprs(i,k+1))/RG
1870	ENDDO
1871	ENDDO
1872	ENDIF
1873
1874	! Calcule de vitesse verticale a partir de flux de masse verticale
1875	DO k = 1, klev
1876	DO i = 1, klon
1877	omega(i,k) = RG*flxmass_w(i,k) / airephy(i)
1878	END DO
1879	END DO
1880	if (prt_level.ge.1) write(lunout,*) 'omega(igout, :) = ', &
1881	omega(igout, :)
1882
1883	IF (iflag_con.EQ.1) THEN
1884	abort_message ='reactiver le call conlmd dans physiq.F'
1885	CALL abort_gcm (modname,abort_message,1)
1886	! CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
1887	! . d_t_con, d_q_con,
1888	! . rain_con, snow_con, ibas_con, itop_con)
1889	ELSE IF (iflag_con.EQ.2) THEN
1890	CALL conflx(dtime, paprs, pplay, t_seri, q_seri, &
1891	conv_t, conv_q, -evap, omega, &
1892	d_t_con, d_q_con, rain_con, snow_con, &
1893	pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
1894	kcbot, kctop, kdtop, pmflxr, pmflxs)
1895	d_u_con = 0.
1896	d_v_con = 0.
1897
1898	WHERE (rain_con < 0.) rain_con = 0.
1899	WHERE (snow_con < 0.) snow_con = 0.
1900	DO i = 1, klon
1901	ibas_con(i) = klev+1 - kcbot(i)
1902	itop_con(i) = klev+1 - kctop(i)
1903	ENDDO
1904	ELSE IF (iflag_con.GE.3) THEN
1905	! nb of tracers for the KE convection:
1906	! MAF la partie traceurs est faite dans phytrac
1907	! on met ntra=1 pour limiter les appels mais on peut
1908	! supprimer les calculs / ftra.
1909	ntra = 1
1910
1911	!=====================================================================================
1912	!ajout pour la parametrisation des poches froides:
1913	!calcul de t_wake et t_undi: si pas de poches froides, t_wake=t_undi=t_seri
1914	do k=1,klev
1915	do i=1,klon
1916	if (iflag_wake>=1) then
1917	t_wake(i,k) = t_seri(i,k) &
1918	+(1-wake_s(i))*wake_deltat(i,k)
1919	q_wake(i,k) = q_seri(i,k) &
1920	+(1-wake_s(i))*wake_deltaq(i,k)
1921	t_undi(i,k) = t_seri(i,k) &
1922	-wake_s(i)*wake_deltat(i,k)
1923	q_undi(i,k) = q_seri(i,k) &
1924	-wake_s(i)*wake_deltaq(i,k)
1925	else
1926	t_wake(i,k) = t_seri(i,k)
1927	q_wake(i,k) = q_seri(i,k)
1928	t_undi(i,k) = t_seri(i,k)
1929	q_undi(i,k) = q_seri(i,k)
1930	endif
1931	enddo
1932	enddo
1933
1934	!c-- Calcul de l'energie disponible ALE (J/kg) et de la puissance disponible ALP (W/m2)
1935	!c-- pour le soulevement des particules dans le modele convectif
1936	!
1937	do i = 1,klon
1938	ALE(i) = 0.
1939	ALP(i) = 0.
1940	enddo
1941	!
1942	!calcul de ale_wake et alp_wake
1943	if (iflag_wake>=1) then
1944	if (itap .le. it_wape_prescr) then
1945	do i = 1,klon
1946	ale_wake(i) = wape_prescr
1947	alp_wake(i) = fip_prescr
1948	enddo
1949	else
1950	do i = 1,klon
1951	!jyg ALE=WAPE au lieu de ALE = 1/2 Cstar**2
1952	!cc ale_wake(i) = 0.5wake_cstar(i)*2
1953	ale_wake(i) = wake_pe(i)
1954	alp_wake(i) = wake_fip(i)
1955	enddo
1956	endif
1957	else
1958	do i = 1,klon
1959	ale_wake(i) = 0.
1960	alp_wake(i) = 0.
1961	enddo
1962	endif
1963	!combinaison avec ale et alp de couche limite: constantes si pas de couplage, valeurs calculees
1964	!dans le thermique sinon
1965	if (iflag_coupl.eq.0) then
1966	if (debut.and.prt_level.gt.9) &
1967	WRITE(lunout,*)'ALE et ALP imposes'
1968	do i = 1,klon
1969	!on ne couple que ale
1970	! ALE(i) = max(ale_wake(i),Ale_bl(i))
1971	ALE(i) = max(ale_wake(i),ale_bl_prescr)
1972	!on ne couple que alp
1973	! ALP(i) = alp_wake(i) + Alp_bl(i)
1974	ALP(i) = alp_wake(i) + alp_bl_prescr
1975	enddo
1976	else
1977	IF(prt_level>9)WRITE(lunout,*)'ALE et ALP couples au thermique'
1978	! do i = 1,klon
1979	! ALE(i) = max(ale_wake(i),Ale_bl(i))
1980	! avant ALP(i) = alp_wake(i) + Alp_bl(i)
1981	! ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
1982	! write(20,*)'ALE',ALE(i),Ale_bl(i),ale_wake(i)
1983	! write(21,*)'ALP',ALP(i),Alp_bl(i),alp_wake(i)
1984	! enddo
1985
1986	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1987	! Modif FH 2010/04/27. Sans doute temporaire.
1988	! Deux options pour le alp_offset : constant si >?? 0 ou proportionnel ??a
1989	! w si <0
1990	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1991	do i = 1,klon
1992	ALE(i) = max(ale_wake(i),Ale_bl(i))
1993	!cc nrlmd le 10/04/2012----------Stochastic triggering--------------
1994	if (iflag_trig_bl.ge.1) then
1995	ALE(i) = max(ale_wake(i),Ale_bl_trig(i))
1996	endif
1997	!cc fin nrlmd le 10/04/2012
1998	if (alp_offset>=0.) then
1999	ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
2000	else
2001	ALP(i)=alp_wake(i)+Alp_bl(i)+alp_offset*min(omega(i,6),0.)
2002	if (alp(i)<0.) then
2003	print*,'ALP ',alp(i),alp_wake(i) &
2004	,Alp_bl(i),alp_offset*min(omega(i,6),0.)
2005	endif
2006	endif
2007	enddo
2008	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2009
2010	endif
2011	do i=1,klon
2012	if (alp(i)>alp_max) then
2013	IF(prt_level>9)WRITE(lunout,*) &
2014	'WARNING SUPER ALP (seuil=',alp_max, &
2015	'): i, alp, alp_wake,ale',i,alp(i),alp_wake(i),ale(i)
2016	alp(i)=alp_max
2017	endif
2018	if (ale(i)>ale_max) then
2019	IF(prt_level>9)WRITE(lunout,*) &
2020	'WARNING SUPER ALE (seuil=',ale_max, &
2021	'): i, alp, alp_wake,ale',i,ale(i),ale_wake(i),alp(i)
2022	ale(i)=ale_max
2023	endif
2024	enddo
2025
2026	!fin calcul ale et alp
2027	!=================================================================================================
2028
2029
2030	! sb, oct02:
2031	! Schema de convection modularise et vectorise:
2032	! (driver commun aux versions 3 et 4)
2033	!
2034	IF (ok_cvl) THEN ! new driver for convectL
2035
2036	IF (type_trac == 'repr') THEN
2037	nbtr_tmp=ntra
2038	ELSE
2039	nbtr_tmp=nbtr
2040	END IF
2041	!jyg iflag_con est dans clesphys
2042	!c CALL concvl (iflag_con,iflag_clos,
2043	CALL concvl (iflag_clos, &
2044	dtime,paprs,pplay,t_undi,q_undi, &
2045	t_wake,q_wake,wake_s, &
2046	u_seri,v_seri,tr_seri,nbtr_tmp, &
2047	ALE,ALP, &
2048	sig1,w01, &
2049	d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
2050	rain_con, snow_con, ibas_con, itop_con, sigd, &
2051	ema_cbmf,plcl,plfc,wbeff,upwd,dnwd,dnwd0, &
2052	Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl, &
2053	pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd, &
2054	! RomP >>>
2055	!! . pmflxr,pmflxs,da,phi,mp,
2056	!! . ftd,fqd,lalim_conv,wght_th)
2057	pmflxr,pmflxs,da,phi,mp,phi2,d1a,dam,sij,clw,elij, &
2058	ftd,fqd,lalim_conv,wght_th, &
2059	ev, ep,epmlmMm,eplaMm, &
2060	wdtrainA,wdtrainM)
2061	! RomP <<<
2062
2063	!IM begin
2064	! print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
2065	! .dnwd0(1,1),ftd(1,1),fqd(1,1)
2066	!IM end
2067	!IM cf. FH
2068	clwcon0=qcondc
2069	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
2070
2071	do i = 1, klon
2072	if (iflagctrl(i).le.1) itau_con(i)=itau_con(i)+1
2073	enddo
2074
2075	ELSE ! ok_cvl
2076
2077	! MAF conema3 ne contient pas les traceurs
2078	CALL conema3 (dtime, &
2079	paprs,pplay,t_seri,q_seri, &
2080	u_seri,v_seri,tr_seri,ntra, &
2081	sig1,w01, &
2082	d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
2083	rain_con, snow_con, ibas_con, itop_con, &
2084	upwd,dnwd,dnwd0,bas,top, &
2085	Ma,cape,tvp,rflag, &
2086	pbase &
2087	,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr &
2088	,clwcon0)
2089
2090	ENDIF ! ok_cvl
2091
2092	!
2093	! Correction precip
2094	rain_con = rain_con * cvl_corr
2095	snow_con = snow_con * cvl_corr
2096	!
2097
2098	IF (.NOT. ok_gust) THEN
2099	do i = 1, klon
2100	wd(i)=0.0
2101	enddo
2102	ENDIF
2103
2104	! =================================================================== c
2105	! Calcul des proprietes des nuages convectifs
2106	!
2107
2108	! calcul des proprietes des nuages convectifs
2109	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
2110	call clouds_gno &
2111	(klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
2112
2113	! =================================================================== c
2114
2115	DO i = 1, klon
2116	itop_con(i) = min(max(itop_con(i),1),klev)
2117	ibas_con(i) = min(max(ibas_con(i),1),itop_con(i))
2118	ENDDO
2119
2120	DO i = 1, klon
2121	ema_pcb(i) = paprs(i,ibas_con(i))
2122	ENDDO
2123	DO i = 1, klon
2124	! L'idicage de itop_con peut cacher un pb potentiel
2125	! FH sous la dictee de JYG, CR
2126	ema_pct(i) = paprs(i,itop_con(i)+1)
2127
2128	if (itop_con(i).gt.klev-3) then
2129	if(prt_level >= 9) then
2130	write(lunout,*)'La convection monte trop haut '
2131	write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
2132	endif
2133	endif
2134	ENDDO
2135	ELSE IF (iflag_con.eq.0) THEN
2136	write(lunout,*) 'On n appelle pas la convection'
2137	clwcon0=0.
2138	rnebcon0=0.
2139	d_t_con=0.
2140	d_q_con=0.
2141	d_u_con=0.
2142	d_v_con=0.
2143	rain_con=0.
2144	snow_con=0.
2145	bas=1
2146	top=1
2147	ELSE
2148	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
2149	call abort_gcm("physiq", "", 1)
2150	ENDIF
2151
2152	! CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
2153	! . d_u_con, d_v_con)
2154
2155	!-----------------------------------------------------------------------------------------
2156	! ajout des tendances de la diffusion turbulente
2157	CALL add_phys_tend(d_u_con,d_v_con,d_t_con,d_q_con,dql0,paprs,'con')
2158	!-----------------------------------------------------------------------------------------
2159
2160	if (mydebug) then
2161	call writefield_phy('u_seri',u_seri,llm)
2162	call writefield_phy('v_seri',v_seri,llm)
2163	call writefield_phy('t_seri',t_seri,llm)
2164	call writefield_phy('q_seri',q_seri,llm)
2165	endif
2166
2167	!IM
2168	IF (ip_ebil_phy.ge.2) THEN
2169	ztit='after convect'
2170	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
2171	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
2172	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2173	call diagphy(airephy,ztit,ip_ebil_phy &
2174	, zero_v, zero_v, zero_v, zero_v, zero_v &
2175	, zero_v, rain_con, snow_con, ztsol &
2176	, d_h_vcol, d_qt, d_ec &
2177	, fs_bound, fq_bound )
2178	END IF
2179	!
2180	IF (check) THEN
2181	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2182	WRITE(lunout,*)"aprescon=", za
2183	zx_t = 0.0
2184	za = 0.0
2185	DO i = 1, klon
2186	za = za + airephy(i)/REAL(klon)
2187	zx_t = zx_t + (rain_con(i)+ &
2188	snow_con(i))*airephy(i)/REAL(klon)
2189	ENDDO
2190	zx_t = zx_t/za*dtime
2191	WRITE(lunout,*)"Precip=", zx_t
2192	ENDIF
2193	IF (zx_ajustq) THEN
2194	DO i = 1, klon
2195	z_apres(i) = 0.0
2196	ENDDO
2197	DO k = 1, klev
2198	DO i = 1, klon
2199	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k)) &
2200	*(paprs(i,k)-paprs(i,k+1))/RG
2201	ENDDO
2202	ENDDO
2203	DO i = 1, klon
2204	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime) &
2205	/z_apres(i)
2206	ENDDO
2207	DO k = 1, klev
2208	DO i = 1, klon
2209	IF (z_factor(i).GT.(1.0+1.0E-08) .OR. &
2210	z_factor(i).LT.(1.0-1.0E-08)) THEN
2211	q_seri(i,k) = q_seri(i,k) * z_factor(i)
2212	ENDIF
2213	ENDDO
2214	ENDDO
2215	ENDIF
2216	zx_ajustq=.FALSE.
2217
2218	!
2219	!=============================================================================
2220	!RR:Evolution de la poche froide: on ne fait pas de separation wake/env
2221	!pour la couche limite diffuse pour l instant
2222	!
2223	if (iflag_wake>=1) then
2224	DO k=1,klev
2225	DO i=1,klon
2226	dt_dwn(i,k) = ftd(i,k)
2227	wdt_PBL(i,k) = 0.
2228	dq_dwn(i,k) = fqd(i,k)
2229	wdq_PBL(i,k) = 0.
2230	M_dwn(i,k) = dnwd0(i,k)
2231	M_up(i,k) = upwd(i,k)
2232	dt_a(i,k) = d_t_con(i,k)/dtime - ftd(i,k)
2233	udt_PBL(i,k) = 0.
2234	dq_a(i,k) = d_q_con(i,k)/dtime - fqd(i,k)
2235	udq_PBL(i,k) = 0.
2236	ENDDO
2237	ENDDO
2238
2239	if (iflag_wake==2) then
2240	ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
2241	DO k = 1,klev
2242	dt_dwn(:,k)= dt_dwn(:,k)+ &
2243	ok_wk_lsp(:)*(d_t_eva(:,k)+d_t_lsc(:,k))/dtime
2244	dq_dwn(:,k)= dq_dwn(:,k)+ &
2245	ok_wk_lsp(:)*(d_q_eva(:,k)+d_q_lsc(:,k))/dtime
2246	ENDDO
2247	endif
2248	!
2249	!calcul caracteristiques de la poche froide
2250	call calWAKE (paprs,pplay,dtime &
2251	,t_seri,q_seri,omega &
2252	,dt_dwn,dq_dwn,M_dwn,M_up &
2253	,dt_a,dq_a,sigd &
2254	,wdt_PBL,wdq_PBL &
2255	,udt_PBL,udq_PBL &
2256	,wake_deltat,wake_deltaq,wake_dth &
2257	,wake_h,wake_s,wake_dens &
2258	,wake_pe,wake_fip,wake_gfl &
2259	,dt_wake,dq_wake &
2260	,wake_k, t_undi,q_undi &
2261	,wake_omgbdth,wake_dp_omgb &
2262	,wake_dtKE,wake_dqKE &
2263	,wake_dtPBL,wake_dqPBL &
2264	,wake_omg,wake_dp_deltomg &
2265	,wake_spread,wake_Cstar,wake_d_deltat_gw &
2266	,wake_ddeltat,wake_ddeltaq)
2267	!
2268	!-----------------------------------------------------------------------------------------
2269	! ajout des tendances des poches froides
2270	! Faire rapidement disparaitre l'ancien dt_wake pour garder un d_t_wake
2271	! coherent avec les autres d_t_...
2272	d_t_wake(:,:)=dt_wake(:,:)*dtime
2273	d_q_wake(:,:)=dq_wake(:,:)*dtime
2274	CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,paprs,'wake')
2275	!-----------------------------------------------------------------------------------------
2276
2277	endif
2278	!
2279	!===================================================================
2280	!JYG
2281	IF (ip_ebil_phy.ge.2) THEN
2282	ztit='after wake'
2283	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
2284	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
2285	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2286	call diagphy(airephy,ztit,ip_ebil_phy &
2287	, zero_v, zero_v, zero_v, zero_v, zero_v &
2288	, zero_v, zero_v, zero_v, ztsol &
2289	, d_h_vcol, d_qt, d_ec &
2290	, fs_bound, fq_bound )
2291	END IF
2292
2293	! print*,'apres callwake iflag_cldcon=', iflag_cldcon
2294	!
2295	!===================================================================
2296	! Convection seche (thermiques ou ajustement)
2297	!===================================================================
2298	!
2299	call stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri &
2300	,seuil_inversion,weak_inversion,dthmin)
2301
2302
2303
2304	d_t_ajsb(:,:)=0.
2305	d_q_ajsb(:,:)=0.
2306	d_t_ajs(:,:)=0.
2307	d_u_ajs(:,:)=0.
2308	d_v_ajs(:,:)=0.
2309	d_q_ajs(:,:)=0.
2310	clwcon0th(:,:)=0.
2311	!
2312	! fm_therm(:,:)=0.
2313	! entr_therm(:,:)=0.
2314	! detr_therm(:,:)=0.
2315	!
2316	IF(prt_level>9)WRITE(lunout,*) &
2317	'AVANT LA CONVECTION SECHE , iflag_thermals=' &
2318	,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2319	if(iflag_thermals<0) then
2320	! Rien
2321	! ====
2322	IF(prt_level>9)WRITE(lunout,*)'pas de convection seche'
2323
2324
2325	else
2326
2327	! Thermiques
2328	! ==========
2329	IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals=' &
2330	,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2331
2332
2333	!cc nrlmd le 10/04/2012
2334	DO k=1,klev+1
2335	DO i=1,klon
2336	pbl_tke_input(i,k,is_oce)=pbl_tke(i,k,is_oce)
2337	pbl_tke_input(i,k,is_ter)=pbl_tke(i,k,is_ter)
2338	pbl_tke_input(i,k,is_lic)=pbl_tke(i,k,is_lic)
2339	pbl_tke_input(i,k,is_sic)=pbl_tke(i,k,is_sic)
2340	ENDDO
2341	ENDDO
2342	!cc fin nrlmd le 10/04/2012
2343
2344	if (iflag_thermals>=1) then
2345	call calltherm(pdtphys &
2346	,pplay,paprs,pphi,weak_inversion &
2347	,u_seri,v_seri,t_seri,q_seri,zqsat,debut &
2348	,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs &
2349	,fm_therm,entr_therm,detr_therm &
2350	,zqasc,clwcon0th,lmax_th,ratqscth &
2351	,ratqsdiff,zqsatth &
2352	!on rajoute ale et alp, et les caracteristiques de la couche alim
2353	,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca &
2354	,ztv,zpspsk,ztla,zthl &
2355	!cc nrlmd le 10/04/2012
2356	,pbl_tke_input,pctsrf,omega,airephy &
2357	,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
2358	,n2,s2,ale_bl_stat &
2359	,therm_tke_max,env_tke_max &
2360	,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
2361	,alp_bl_conv,alp_bl_stat &
2362	!cc fin nrlmd le 10/04/2012
2363	,zqla,ztva )
2364
2365	!cc nrlmd le 10/04/2012
2366	!-----------Stochastic triggering-----------
2367	if (iflag_trig_bl.ge.1) then
2368	!
2369	IF (prt_level .GE. 10) THEN
2370	print *,'cin, ale_bl_stat, alp_bl_stat ', &
2371	cin, ale_bl_stat, alp_bl_stat
2372	ENDIF
2373
2374
2375	!----Initialisations
2376	do i=1,klon
2377	proba_notrig(i)=1.
2378	random_notrig(i)=1e6ale_bl_stat(i)-int(1e6ale_bl_stat(i))
2379	if ( ale_bl_trig(i) .lt. abs(cin(i))+1.e-10 ) then
2380	tau_trig(i)=tau_trig_shallow
2381	else
2382	tau_trig(i)=tau_trig_deep
2383	endif
2384	enddo
2385	!
2386	IF (prt_level .GE. 10) THEN
2387	print *,'random_notrig, tau_trig ', &
2388	random_notrig, tau_trig
2389	print *,'s_trig,s2,n2 ', &
2390	s_trig,s2,n2
2391	ENDIF
2392
2393	!Option pour re-activer l'ancien calcul de Ale_bl (iflag_trig_bl=2)
2394	IF (iflag_trig_bl.eq.1) then
2395
2396	!----Tirage al\'eatoire et calcul de ale_bl_trig
2397	do i=1,klon
2398	if ( (ale_bl_stat(i) .gt. abs(cin(i))+1.e-10) ) then
2399	proba_notrig(i)=(1.-exp(-s_trig/s2(i)))** &
2400	(n2(i)*dtime/tau_trig(i))
2401	! print *, 'proba_notrig(i) ',proba_notrig(i)
2402	if (random_notrig(i) .ge. proba_notrig(i)) then
2403	ale_bl_trig(i)=ale_bl_stat(i)
2404	else
2405	ale_bl_trig(i)=0.
2406	endif
2407	else
2408	proba_notrig(i)=1.
2409	random_notrig(i)=0.
2410	ale_bl_trig(i)=0.
2411	endif
2412	enddo
2413
2414	ELSE IF (iflag_trig_bl.eq.2) then
2415
2416	do i=1,klon
2417	if ( (Ale_bl(i) .gt. abs(cin(i))+1.e-10) ) then
2418	proba_notrig(i)=(1.-exp(-s_trig/s2(i)))** &
2419	(n2(i)*dtime/tau_trig(i))
2420	! print *, 'proba_notrig(i) ',proba_notrig(i)
2421	if (random_notrig(i) .ge. proba_notrig(i)) then
2422	ale_bl_trig(i)=Ale_bl(i)
2423	else
2424	ale_bl_trig(i)=0.
2425	endif
2426	else
2427	proba_notrig(i)=1.
2428	random_notrig(i)=0.
2429	ale_bl_trig(i)=0.
2430	endif
2431	enddo
2432
2433	ENDIF
2434
2435	!
2436	IF (prt_level .GE. 10) THEN
2437	print *,'proba_notrig, ale_bl_trig ', &
2438	proba_notrig, ale_bl_trig
2439	ENDIF
2440
2441	endif !(iflag_trig_bl)
2442
2443	!-----------Statistical closure-----------
2444	if (iflag_clos_bl.eq.1) then
2445
2446	do i=1,klon
2447	!CR: alp probabiliste
2448	if (ale_bl_trig(i).gt.0.) then
2449	alp_bl(i)=alp_bl(i)/(1.-min(proba_notrig(i),0.999))
2450	endif
2451	enddo
2452
2453	else if (iflag_clos_bl.eq.2) then
2454
2455	!CR: alp calculee dans thermcell_main
2456	do i=1,klon
2457	alp_bl(i)=alp_bl_stat(i)
2458	enddo
2459
2460	else
2461
2462	alp_bl_stat(:)=0.
2463
2464	endif !(iflag_clos_bl)
2465
2466	IF (prt_level .GE. 10) THEN
2467	print *,'ale_bl_trig, alp_bl_stat ',ale_bl_trig, alp_bl_stat
2468	ENDIF
2469
2470	!cc fin nrlmd le 10/04/2012
2471
2472	! ----------------------------------------------------------------------
2473	! Transport de la TKE par les panaches thermiques.
2474	! FH : 2010/02/01
2475	! if (iflag_pbl.eq.10) then
2476	! call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm,
2477	! s rg,paprs,pbl_tke)
2478	! endif
2479	! ----------------------------------------------------------------------
2480	!IM/FH: 2011/02/23
2481	! Couplage Thermiques/Emanuel seulement si T<0
2482	if (iflag_coupl==2) then
2483	print*,'Couplage Thermiques/Emanuel seulement si T<0'
2484	do i=1,klon
2485	if (t_seri(i,lmax_th(i))>273.) then
2486	Ale_bl(i)=0.
2487	endif
2488	enddo
2489	endif
2490
2491	do i=1,klon
2492	! zmax_th(i)=pphi(i,lmax_th(i))/rg
2493	!CR:04/05/12:correction calcul zmax
2494	zmax_th(i)=zmax0(i)
2495	enddo
2496
2497	endif
2498
2499
2500	! Ajustement sec
2501	! ==============
2502
2503	! Dans le cas o\`u on active les thermiques, on fait partir l'ajustement
2504	! a partir du sommet des thermiques.
2505	! Dans le cas contraire, on demarre au niveau 1.
2506
2507	if (iflag_thermals>=13.or.iflag_thermals<=0) then
2508
2509	if(iflag_thermals.eq.0) then
2510	IF(prt_level>9)WRITE(lunout,*)'ajsec'
2511	limbas(:)=1
2512	else
2513	limbas(:)=lmax_th(:)
2514	endif
2515
2516	! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
2517	! pour des test de convergence numerique.
2518	! Le nouveau ajsec est a priori mieux, meme pour le cas
2519	! iflag_thermals = 0 (l'ancienne version peut faire des tendances
2520	! non nulles numeriquement pour des mailles non concernees.
2521
2522	if (iflag_thermals==0) then
2523	! Calling adjustment alone (but not the thermal plume model)
2524	CALL ajsec_convV2(paprs, pplay, t_seri,q_seri &
2525	, d_t_ajsb, d_q_ajsb)
2526	else if (iflag_thermals>0) then
2527	! Calling adjustment above the top of thermal plumes
2528	CALL ajsec(paprs, pplay, t_seri,q_seri,limbas &
2529	, d_t_ajsb, d_q_ajsb)
2530	endif
2531
2532	!-----------------------------------------------------------------------------------------
2533	! ajout des tendances de l'ajustement sec ou des thermiques
2534	CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,paprs,'ajsb')
2535	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
2536	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
2537
2538	!-----------------------------------------------------------------------------------------
2539
2540	endif
2541
2542	endif
2543	!
2544	!===================================================================
2545	!IM
2546	IF (ip_ebil_phy.ge.2) THEN
2547	ztit='after dry_adjust'
2548	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
2549	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
2550	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2551	call diagphy(airephy,ztit,ip_ebil_phy &
2552	, zero_v, zero_v, zero_v, zero_v, zero_v &
2553	, zero_v, zero_v, zero_v, ztsol &
2554	, d_h_vcol, d_qt, d_ec &
2555	, fs_bound, fq_bound )
2556	END IF
2557
2558
2559	!-------------------------------------------------------------------------
2560	! Computation of ratqs, the width (normalized) of the subrid scale
2561	! water distribution
2562	CALL calcratqs(klon,klev,prt_level,lunout, &
2563	iflag_ratqs,iflag_con,iflag_cldcon,pdtphys, &
2564	ratqsbas,ratqshaut,tau_ratqs,fact_cldcon, &
2565	ptconv,ptconvth,clwcon0th, rnebcon0th, &
2566	paprs,pplay,q_seri,zqsat,fm_therm, &
2567	ratqs,ratqsc)
2568
2569
2570	!
2571	! Appeler le processus de condensation a grande echelle
2572	! et le processus de precipitation
2573	!-------------------------------------------------------------------------
2574	IF (prt_level .GE.10) THEN
2575	print *,' ->fisrtilp '
2576	ENDIF
2577	!-------------------------------------------------------------------------
2578	CALL fisrtilp(dtime,paprs,pplay, &
2579	t_seri, q_seri,ptconv,ratqs, &
2580	d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq, &
2581	rain_lsc, snow_lsc, &
2582	pfrac_impa, pfrac_nucl, pfrac_1nucl, &
2583	frac_impa, frac_nucl, beta_prec_fisrt, &
2584	prfl, psfl, rhcl, &
2585	zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cldcon, &
2586	iflag_ice_thermo)
2587
2588	WHERE (rain_lsc < 0) rain_lsc = 0.
2589	WHERE (snow_lsc < 0) snow_lsc = 0.
2590	!-----------------------------------------------------------------------------------------
2591	! ajout des tendances de la diffusion turbulente
2592	CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,paprs,'lsc')
2593	!-----------------------------------------------------------------------------------------
2594	DO k = 1, klev
2595	DO i = 1, klon
2596	cldfra(i,k) = rneb(i,k)
2597	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
2598	ENDDO
2599	ENDDO
2600	IF (check) THEN
2601	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2602	WRITE(lunout,*)"apresilp=", za
2603	zx_t = 0.0
2604	za = 0.0
2605	DO i = 1, klon
2606	za = za + airephy(i)/REAL(klon)
2607	zx_t = zx_t + (rain_lsc(i) &
2608	+ snow_lsc(i))*airephy(i)/REAL(klon)
2609	ENDDO
2610	zx_t = zx_t/za*dtime
2611	WRITE(lunout,*)"Precip=", zx_t
2612	ENDIF
2613	!IM
2614	IF (ip_ebil_phy.ge.2) THEN
2615	ztit='after fisrt'
2616	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
2617	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
2618	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2619	call diagphy(airephy,ztit,ip_ebil_phy &
2620	, zero_v, zero_v, zero_v, zero_v, zero_v &
2621	, zero_v, rain_lsc, snow_lsc, ztsol &
2622	, d_h_vcol, d_qt, d_ec &
2623	, fs_bound, fq_bound )
2624	END IF
2625
2626	if (mydebug) then
2627	call writefield_phy('u_seri',u_seri,llm)
2628	call writefield_phy('v_seri',v_seri,llm)
2629	call writefield_phy('t_seri',t_seri,llm)
2630	call writefield_phy('q_seri',q_seri,llm)
2631	endif
2632
2633	!
2634	!-------------------------------------------------------------------
2635	! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
2636	!-------------------------------------------------------------------
2637
2638	! 1. NUAGES CONVECTIFS
2639	!
2640	!IM cf FH
2641	! IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
2642	IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
2643	snow_tiedtke=0.
2644	! print*,'avant calcul de la pseudo precip '
2645	! print*,'iflag_cldcon',iflag_cldcon
2646	if (iflag_cldcon.eq.-1) then
2647	rain_tiedtke=rain_con
2648	else
2649	! print*,'calcul de la pseudo precip '
2650	rain_tiedtke=0.
2651	! print*,'calcul de la pseudo precip 0'
2652	do k=1,klev
2653	do i=1,klon
2654	if (d_q_con(i,k).lt.0.) then
2655	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys &
2656	*(paprs(i,k)-paprs(i,k+1))/rg
2657	endif
2658	enddo
2659	enddo
2660	endif
2661	!
2662	! call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
2663	!
2664
2665	! Nuages diagnostiques pour Tiedtke
2666	CALL diagcld1(paprs,pplay, &
2667	!IM cf FH . rain_con,snow_con,ibas_con,itop_con,
2668	rain_tiedtke,snow_tiedtke,ibas_con,itop_con, &
2669	diafra,dialiq)
2670	DO k = 1, klev
2671	DO i = 1, klon
2672	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2673	cldliq(i,k) = dialiq(i,k)
2674	cldfra(i,k) = diafra(i,k)
2675	ENDIF
2676	ENDDO
2677	ENDDO
2678
2679	ELSE IF (iflag_cldcon.ge.3) THEN
2680	! On prend pour les nuages convectifs le max du calcul de la
2681	! convection et du calcul du pas de temps precedent diminue d'un facteur
2682	! facttemps
2683	facteur = pdtphys *facttemps
2684	do k=1,klev
2685	do i=1,klon
2686	rnebcon(i,k)=rnebcon(i,k)*facteur
2687	if (rnebcon0(i,k)clwcon0(i,k).gt.rnebcon(i,k)clwcon(i,k)) &
2688	then
2689	rnebcon(i,k)=rnebcon0(i,k)
2690	clwcon(i,k)=clwcon0(i,k)
2691	endif
2692	enddo
2693	enddo
2694
2695	!
2696	!jq - introduce the aerosol direct and first indirect radiative forcings
2697	!jq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
2698	IF (flag_aerosol .gt. 0) THEN
2699	IF (.NOT. aerosol_couple) &
2700	CALL readaerosol_optic( &
2701	debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
2702	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
2703	mass_solu_aero, mass_solu_aero_pi, &
2704	tau_aero, piz_aero, cg_aero, &
2705	tausum_aero, tau3d_aero)
2706	ELSE
2707	tausum_aero(:,:,:) = 0.
2708	tau_aero(:,:,:,:) = 0.
2709	piz_aero(:,:,:,:) = 0.
2710	cg_aero(:,:,:,:) = 0.
2711	ENDIF
2712	!
2713	!--STRAT AEROSOL
2714	!--updates tausum_aero,tau_aero,piz_aero,cg_aero
2715	IF (flag_aerosol_strat) THEN
2716	PRINT *,'appel a readaerosolstrat', mth_cur
2717	CALL readaerosolstrato(debut)
2718	ENDIF
2719	!--fin STRAT AEROSOL
2720
2721
2722	! On prend la somme des fractions nuageuses et des contenus en eau
2723
2724	if (iflag_cldcon>=5) then
2725
2726	do k=1,klev
2727	ptconvth(:,k)=fm_therm(:,k+1)>0.
2728	enddo
2729
2730	if (iflag_coupl==4) then
2731
2732	! Dans le cas iflag_coupl==4, on prend la somme des convertures
2733	! convectives et lsc dans la partie des thermiques
2734	! Le controle par iflag_coupl est peut etre provisoire.
2735	do k=1,klev
2736	do i=1,klon
2737	if (ptconv(i,k).and.ptconvth(i,k)) then
2738	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
2739	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
2740	else if (ptconv(i,k)) then
2741	cldfra(i,k)=rnebcon(i,k)
2742	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
2743	endif
2744	enddo
2745	enddo
2746
2747	else if (iflag_coupl==5) then
2748	do k=1,klev
2749	do i=1,klon
2750	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
2751	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
2752	enddo
2753	enddo
2754
2755	else
2756
2757	! Si on est sur un point touche par la convection profonde et pas
2758	! par les thermiques, on prend la couverture nuageuse et l'eau nuageuse
2759	! de la convection profonde.
2760
2761	!IM/FH: 2011/02/23
2762	! definition des points sur lesquels ls thermiques sont actifs
2763
2764	do k=1,klev
2765	do i=1,klon
2766	if (ptconv(i,k).and. .not. ptconvth(i,k)) then
2767	cldfra(i,k)=rnebcon(i,k)
2768	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
2769	endif
2770	enddo
2771	enddo
2772
2773	endif
2774
2775	else
2776
2777	! Ancienne version
2778	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
2779	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
2780	endif
2781
2782	ENDIF
2783
2784	! plulsc(:)=0.
2785	! do k=1,klev,-1
2786	! do i=1,klon
2787	! zzz=prfl(:,k)+psfl(:,k)
2788	! if (.not.ptconvth.zzz.gt.0.)
2789	! enddo prfl, psfl,
2790	! enddo
2791	!
2792	! 2. NUAGES STARTIFORMES
2793	!
2794	IF (ok_stratus) THEN
2795	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
2796	DO k = 1, klev
2797	DO i = 1, klon
2798	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2799	cldliq(i,k) = dialiq(i,k)
2800	cldfra(i,k) = diafra(i,k)
2801	ENDIF
2802	ENDDO
2803	ENDDO
2804	ENDIF
2805	!
2806	! Precipitation totale
2807	!
2808	DO i = 1, klon
2809	rain_fall(i) = rain_con(i) + rain_lsc(i)
2810	snow_fall(i) = snow_con(i) + snow_lsc(i)
2811	ENDDO
2812	!IM
2813	IF (ip_ebil_phy.ge.2) THEN
2814	ztit="after diagcld"
2815	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
2816	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
2817	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2818	call diagphy(airephy,ztit,ip_ebil_phy &
2819	, zero_v, zero_v, zero_v, zero_v, zero_v &
2820	, zero_v, zero_v, zero_v, ztsol &
2821	, d_h_vcol, d_qt, d_ec &
2822	, fs_bound, fq_bound )
2823	END IF
2824	!
2825	! Calculer l'humidite relative pour diagnostique
2826	!
2827	DO k = 1, klev
2828	DO i = 1, klon
2829	zx_t = t_seri(i,k)
2830	IF (thermcep) THEN
2831	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2832	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2833	zx_qs = MIN(0.5,zx_qs)
2834	zcor = 1./(1.-retv*zx_qs)
2835	zx_qs = zx_qs*zcor
2836	ELSE
2837	IF (zx_t.LT.t_coup) THEN
2838	zx_qs = qsats(zx_t)/pplay(i,k)
2839	ELSE
2840	zx_qs = qsatl(zx_t)/pplay(i,k)
2841	ENDIF
2842	ENDIF
2843	zx_rh(i,k) = q_seri(i,k)/zx_qs
2844	zqsat(i,k)=zx_qs
2845	ENDDO
2846	ENDDO
2847
2848	!IM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
2849	! equivalente a 2m (tpote) pour diagnostique
2850	!
2851	DO i = 1, klon
2852	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
2853	IF (thermcep) THEN
2854	IF(zt2m(i).LT.RTT) then
2855	Lheat=RLSTT
2856	ELSE
2857	Lheat=RLVTT
2858	ENDIF
2859	ELSE
2860	IF (zt2m(i).LT.RTT) THEN
2861	Lheat=RLSTT
2862	ELSE
2863	Lheat=RLVTT
2864	ENDIF
2865	ENDIF
2866	tpote(i) = tpot(i)* &
2867	EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
2868	ENDDO
2869
2870	IF (type_trac == 'inca') THEN
2871	#ifdef INCA
2872	CALL VTe(VTphysiq)
2873	CALL VTb(VTinca)
2874	calday = REAL(days_elapsed + 1) + jH_cur
2875
2876	call chemtime(itap+itau_phy-1, date0, dtime)
2877	IF (config_inca == 'aero') THEN
2878	CALL AEROSOL_METEO_CALC( &
2879	calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs, &
2880	prfl,psfl,pctsrf,airephy,rlat,rlon,u10m,v10m)
2881	END IF
2882
2883	zxsnow_dummy(:) = 0.0
2884
2885	CALL chemhook_begin (calday, &
2886	days_elapsed+1, &
2887	jH_cur, &
2888	pctsrf(1,1), &
2889	rlat, &
2890	rlon, &
2891	airephy, &
2892	paprs, &
2893	pplay, &
2894	coefh(1:klon,1:klev,is_ave), &
2895	pphi, &
2896	t_seri, &
2897	u, &
2898	v, &
2899	wo(:, :, 1), &
2900	q_seri, &
2901	zxtsol, &
2902	zxsnow_dummy, &
2903	solsw, &
2904	albsol1, &
2905	rain_fall, &
2906	snow_fall, &
2907	itop_con, &
2908	ibas_con, &
2909	cldfra, &
2910	iim, &
2911	jjm, &
2912	tr_seri, &
2913	ftsol, &
2914	paprs, &
2915	cdragh, &
2916	cdragm, &
2917	pctsrf, &
2918	pdtphys, &
2919	itap)
2920
2921	CALL VTe(VTinca)
2922	CALL VTb(VTphysiq)
2923	#endif
2924	END IF !type_trac = inca
2925	!
2926	! Calculer les parametres optiques des nuages et quelques
2927	! parametres pour diagnostiques:
2928	!
2929
2930	IF (aerosol_couple) THEN
2931	mass_solu_aero(:,:) = ccm(:,:,1)
2932	mass_solu_aero_pi(:,:) = ccm(:,:,2)
2933	END IF
2934
2935	if (ok_newmicro) then
2936	CALL newmicro (ok_cdnc, bl95_b0, bl95_b1, &
2937	paprs, pplay, t_seri, cldliq, cldfra, &
2938	cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
2939	flwp, fiwp, flwc, fiwc, &
2940	mass_solu_aero, mass_solu_aero_pi, &
2941	cldtaupi, re, fl, ref_liq, ref_ice, &
2942	ref_liq_pi, ref_ice_pi)
2943	else
2944	CALL nuage (paprs, pplay, &
2945	t_seri, cldliq, cldfra, cldtau, cldemi, &
2946	cldh, cldl, cldm, cldt, cldq, &
2947	ok_aie, &
2948	mass_solu_aero, mass_solu_aero_pi, &
2949	bl95_b0, bl95_b1, &
2950	cldtaupi, re, fl)
2951	endif
2952	!
2953	!IM betaCRF
2954	!
2955	cldtaurad = cldtau
2956	cldtaupirad = cldtaupi
2957	cldemirad = cldemi
2958
2959	!
2960	if(lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND. &
2961	lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
2962	!
2963	! global
2964	!
2965	DO k=1, klev
2966	DO i=1, klon
2967	if (pplay(i,k).GE.pfree) THEN
2968	beta(i,k) = beta_pbl
2969	else
2970	beta(i,k) = beta_free
2971	endif
2972	if (mskocean_beta) THEN
2973	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
2974	endif
2975	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
2976	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
2977	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
2978	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
2979	ENDDO
2980	ENDDO
2981	!
2982	else
2983	!
2984	! regional
2985	!
2986	DO k=1, klev
2987	DO i=1,klon
2988	!
2989	if (rlon(i).ge.lon1_beta.AND.rlon(i).le.lon2_beta.AND. &
2990	rlat(i).le.lat1_beta.AND.rlat(i).ge.lat2_beta) THEN
2991	if (pplay(i,k).GE.pfree) THEN
2992	beta(i,k) = beta_pbl
2993	else
2994	beta(i,k) = beta_free
2995	endif
2996	if (mskocean_beta) THEN
2997	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
2998	endif
2999	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
3000	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
3001	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
3002	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
3003	endif
3004	!
3005	ENDDO
3006	ENDDO
3007	!
3008	endif
3009	!
3010	! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
3011	!
3012	IF (MOD(itaprad,radpas).EQ.0) THEN
3013
3014	DO i = 1, klon
3015	albsol1(i) = falb1(i,is_oce) * pctsrf(i,is_oce) &
3016	+ falb1(i,is_lic) * pctsrf(i,is_lic) &
3017	+ falb1(i,is_ter) * pctsrf(i,is_ter) &
3018	+ falb1(i,is_sic) * pctsrf(i,is_sic)
3019	albsol2(i) = falb2(i,is_oce) * pctsrf(i,is_oce) &
3020	+ falb2(i,is_lic) * pctsrf(i,is_lic) &
3021	+ falb2(i,is_ter) * pctsrf(i,is_ter) &
3022	+ falb2(i,is_sic) * pctsrf(i,is_sic)
3023	ENDDO
3024
3025	if (mydebug) then
3026	call writefield_phy('u_seri',u_seri,llm)
3027	call writefield_phy('v_seri',v_seri,llm)
3028	call writefield_phy('t_seri',t_seri,llm)
3029	call writefield_phy('q_seri',q_seri,llm)
3030	endif
3031
3032	IF (aerosol_couple) THEN
3033	#ifdef INCA
3034	CALL radlwsw_inca &
3035	(kdlon,kflev,dist, rmu0, fract, solaire, &
3036	paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri, &
3037	wo(:, :, 1), &
3038	cldfrarad, cldemirad, cldtaurad, &
3039	heat,heat0,cool,cool0,radsol,albpla, &
3040	topsw,toplw,solsw,sollw, &
3041	sollwdown, &
3042	topsw0,toplw0,solsw0,sollw0, &
3043	lwdn0, lwdn, lwup0, lwup, &
3044	swdn0, swdn, swup0, swup, &
3045	ok_ade, ok_aie, &
3046	tau_aero, piz_aero, cg_aero, &
3047	topswad_aero, solswad_aero, &
3048	topswad0_aero, solswad0_aero, &
3049	topsw_aero, topsw0_aero, &
3050	solsw_aero, solsw0_aero, &
3051	cldtaupirad, &
3052	topswai_aero, solswai_aero)
3053
3054	#endif
3055	ELSE
3056	!
3057	!IM calcul radiatif pour le cas actuel
3058	!
3059	RCO2 = RCO2_act
3060	RCH4 = RCH4_act
3061	RN2O = RN2O_act
3062	RCFC11 = RCFC11_act
3063	RCFC12 = RCFC12_act
3064	!
3065	IF (prt_level .GE.10) THEN
3066	print *,' ->radlwsw, number 1 '
3067	ENDIF
3068	!
3069	CALL radlwsw &
3070	(dist, rmu0, fract, &
3071	paprs, pplay,zxtsol,albsol1, albsol2, &
3072	t_seri,q_seri,wo, &
3073	cldfrarad, cldemirad, cldtaurad, &
3074	ok_ade.OR.flag_aerosol_strat, ok_aie, flag_aerosol, &
3075	flag_aerosol_strat, &
3076	tau_aero, piz_aero, cg_aero, &
3077	cldtaupirad,new_aod, &
3078	zqsat, flwc, fiwc, &
3079	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
3080	heat,heat0,cool,cool0,radsol,albpla, &
3081	topsw,toplw,solsw,sollw, &
3082	sollwdown, &
3083	topsw0,toplw0,solsw0,sollw0, &
3084	lwdn0, lwdn, lwup0, lwup, &
3085	swdn0, swdn, swup0, swup, &
3086	topswad_aero, solswad_aero, &
3087	topswai_aero, solswai_aero, &
3088	topswad0_aero, solswad0_aero, &
3089	topsw_aero, topsw0_aero, &
3090	solsw_aero, solsw0_aero, &
3091	topswcf_aero, solswcf_aero, &
3092	ZLWFT0_i, ZFLDN0, ZFLUP0, &
3093	ZSWFT0_i, ZFSDN0, ZFSUP0)
3094
3095	!
3096	!IM 2eme calcul radiatif pour le cas perturbe ou au moins un
3097	!IM des taux doit etre different du taux actuel
3098	!IM Par defaut on a les taux perturbes egaux aux taux actuels
3099	!
3100	if (ok_4xCO2atm) then
3101	if (RCO2_per.NE.RCO2_act.OR.RCH4_per.NE.RCH4_act.OR. &
3102	RN2O_per.NE.RN2O_act.OR.RCFC11_per.NE.RCFC11_act.OR. &
3103	RCFC12_per.NE.RCFC12_act) THEN
3104	!
3105	RCO2 = RCO2_per
3106	RCH4 = RCH4_per
3107	RN2O = RN2O_per
3108	RCFC11 = RCFC11_per
3109	RCFC12 = RCFC12_per
3110	!
3111	IF (prt_level .GE.10) THEN
3112	print *,' ->radlwsw, number 2 '
3113	ENDIF
3114	!
3115	CALL radlwsw &
3116	(dist, rmu0, fract, &
3117	paprs, pplay,zxtsol,albsol1, albsol2, &
3118	t_seri,q_seri,wo, &
3119	cldfra, cldemi, cldtau, &
3120	ok_ade.OR.flag_aerosol_strat, ok_aie, flag_aerosol, &
3121	flag_aerosol_strat, &
3122	tau_aero, piz_aero, cg_aero, &
3123	cldtaupi,new_aod, &
3124	zqsat, flwc, fiwc, &
3125	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
3126	heatp,heat0p,coolp,cool0p,radsolp,albplap, &
3127	topswp,toplwp,solswp,sollwp, &
3128	sollwdownp, &
3129	topsw0p,toplw0p,solsw0p,sollw0p, &
3130	lwdn0p, lwdnp, lwup0p, lwupp, &
3131	swdn0p, swdnp, swup0p, swupp, &
3132	topswad_aerop, solswad_aerop, &
3133	topswai_aerop, solswai_aerop, &
3134	topswad0_aerop, solswad0_aerop, &
3135	topsw_aerop, topsw0_aerop, &
3136	solsw_aerop, solsw0_aerop, &
3137	topswcf_aerop, solswcf_aerop, &
3138	ZLWFT0_i, ZFLDN0, ZFLUP0, &
3139	ZSWFT0_i, ZFSDN0, ZFSUP0)
3140	endif
3141	endif
3142	!
3143	ENDIF ! aerosol_couple
3144	itaprad = 0
3145	ENDIF ! MOD(itaprad,radpas)
3146	itaprad = itaprad + 1
3147
3148	IF (iflag_radia.eq.0) THEN
3149	IF (prt_level.ge.9) THEN
3150	PRINT *,'--------------------------------------------------'
3151	PRINT *,'>>>> ATTENTION rayonnement desactive pour ce cas'
3152	PRINT *,'>>>> heat et cool mis a zero '
3153	PRINT *,'--------------------------------------------------'
3154	END IF
3155	heat=0.
3156	cool=0.
3157	sollw=0. ! MPL 01032011
3158	solsw=0.
3159	radsol=0.
3160	swup=0. ! MPL 27102011 pour les fichiers AMMA_profiles et AMMA_scalars
3161	swup0=0.
3162	swdn=0.
3163	swdn0=0.
3164	lwup=0.
3165	lwup0=0.
3166	lwdn=0.
3167	lwdn0=0.
3168	END IF
3169
3170	!
3171	! Ajouter la tendance des rayonnements (tous les pas)
3172	!
3173	DO k = 1, klev
3174	DO i = 1, klon
3175	t_seri(i,k) = t_seri(i,k) &
3176	+ (heat(i,k)-cool(i,k)) * dtime/RDAY
3177	ENDDO
3178	ENDDO
3179	!
3180	if (mydebug) then
3181	call writefield_phy('u_seri',u_seri,llm)
3182	call writefield_phy('v_seri',v_seri,llm)
3183	call writefield_phy('t_seri',t_seri,llm)
3184	call writefield_phy('q_seri',q_seri,llm)
3185	endif
3186
3187	!IM
3188	IF (ip_ebil_phy.ge.2) THEN
3189	ztit='after rad'
3190	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
3191	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
3192	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3193	call diagphy(airephy,ztit,ip_ebil_phy &
3194	, topsw, toplw, solsw, sollw, zero_v &
3195	, zero_v, zero_v, zero_v, ztsol &
3196	, d_h_vcol, d_qt, d_ec &
3197	, fs_bound, fq_bound )
3198	END IF
3199	!
3200	!
3201	! Calculer l'hydrologie de la surface
3202	!
3203	! CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
3204	! . agesno, ftsol,fqsurf,fsnow, ruis)
3205	!
3206
3207	!
3208	! Calculer le bilan du sol et la derive de temperature (couplage)
3209	!
3210	DO i = 1, klon
3211	! bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
3212	! a la demande de JLD
3213	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
3214	ENDDO
3215	!
3216	!moddeblott(jan95)
3217	! Appeler le programme de parametrisation de l'orographie
3218	! a l'echelle sous-maille:
3219	!
3220	IF (prt_level .GE.10) THEN
3221	print *,' call orography ? ', ok_orodr
3222	ENDIF
3223	!
3224	IF (ok_orodr) THEN
3225	!
3226	! selection des points pour lesquels le shema est actif:
3227	igwd=0
3228	DO i=1,klon
3229	itest(i)=0
3230	! IF ((zstd(i).gt.10.0)) THEN
3231	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
3232	itest(i)=1
3233	igwd=igwd+1
3234	idx(igwd)=i
3235	ENDIF
3236	ENDDO
3237	! igwdim=MAX(1,igwd)
3238	!
3239	IF (ok_strato) THEN
3240
3241	CALL drag_noro_strato(klon,klev,dtime,paprs,pplay, &
3242	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
3243	igwd,idx,itest, &
3244	t_seri, u_seri, v_seri, &
3245	zulow, zvlow, zustrdr, zvstrdr, &
3246	d_t_oro, d_u_oro, d_v_oro)
3247
3248	ELSE
3249	CALL drag_noro(klon,klev,dtime,paprs,pplay, &
3250	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
3251	igwd,idx,itest, &
3252	t_seri, u_seri, v_seri, &
3253	zulow, zvlow, zustrdr, zvstrdr, &
3254	d_t_oro, d_u_oro, d_v_oro)
3255	ENDIF
3256	!
3257	! ajout des tendances
3258	!-----------------------------------------------------------------------------------------
3259	! ajout des tendances de la trainee de l'orographie
3260	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,paprs,'oro')
3261	!-----------------------------------------------------------------------------------------
3262	!
3263	ENDIF ! fin de test sur ok_orodr
3264	!
3265	if (mydebug) then
3266	call writefield_phy('u_seri',u_seri,llm)
3267	call writefield_phy('v_seri',v_seri,llm)
3268	call writefield_phy('t_seri',t_seri,llm)
3269	call writefield_phy('q_seri',q_seri,llm)
3270	endif
3271
3272	IF (ok_orolf) THEN
3273	!
3274	! selection des points pour lesquels le shema est actif:
3275	igwd=0
3276	DO i=1,klon
3277	itest(i)=0
3278	IF ((zpic(i)-zmea(i)).GT.100.) THEN
3279	itest(i)=1
3280	igwd=igwd+1
3281	idx(igwd)=i
3282	ENDIF
3283	ENDDO
3284	! igwdim=MAX(1,igwd)
3285	!
3286	IF (ok_strato) THEN
3287
3288	CALL lift_noro_strato(klon,klev,dtime,paprs,pplay, &
3289	rlat,zmea,zstd,zpic,zgam,zthe,zpic,zval, &
3290	igwd,idx,itest, &
3291	t_seri, u_seri, v_seri, &
3292	zulow, zvlow, zustrli, zvstrli, &
3293	d_t_lif, d_u_lif, d_v_lif )
3294
3295	ELSE
3296	CALL lift_noro(klon,klev,dtime,paprs,pplay, &
3297	rlat,zmea,zstd,zpic, &
3298	itest, &
3299	t_seri, u_seri, v_seri, &
3300	zulow, zvlow, zustrli, zvstrli, &
3301	d_t_lif, d_u_lif, d_v_lif)
3302	ENDIF
3303	!
3304	!-----------------------------------------------------------------------------------------
3305	! ajout des tendances de la portance de l'orographie
3306	CALL add_phys_tend(d_u_lif,d_v_lif,d_t_lif,dq0,dql0,paprs,'lif')
3307	!-----------------------------------------------------------------------------------------
3308	!
3309	ENDIF ! fin de test sur ok_orolf
3310	! HINES GWD PARAMETRIZATION
3311
3312	IF (ok_hines) then
3313
3314	CALL hines_gwd(klon,klev,dtime,paprs,pplay, &
3315	rlat,t_seri,u_seri,v_seri, &
3316	zustrhi,zvstrhi, &
3317	d_t_hin, d_u_hin, d_v_hin)
3318	!
3319	! ajout des tendances
3320	CALL add_phys_tend(d_u_hin,d_v_hin,d_t_hin,dq0,dql0,paprs,'hin')
3321
3322	ENDIF
3323
3324	if (ok_gwd_rando) then
3325	call FLOTT_GWD_rando(DTIME, pplay, t_seri, u_seri, v_seri, &
3326	rain_fall + snow_fall, zustr_gwd_rando, zvstr_gwd_rando, &
3327	du_gwd_rando, dv_gwd_rando)
3328	CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0,paprs, &
3329	'flott_gwd_rando')
3330	end if
3331
3332	! STRESS NECESSAIRES: TOUTE LA PHYSIQUE
3333
3334	if (mydebug) then
3335	call writefield_phy('u_seri',u_seri,llm)
3336	call writefield_phy('v_seri',v_seri,llm)
3337	call writefield_phy('t_seri',t_seri,llm)
3338	call writefield_phy('q_seri',q_seri,llm)
3339	endif
3340
3341	DO i = 1, klon
3342	zustrph(i)=0.
3343	zvstrph(i)=0.
3344	ENDDO
3345	DO k = 1, klev
3346	DO i = 1, klon
3347	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime* &
3348	(paprs(i,k)-paprs(i,k+1))/rg
3349	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime* &
3350	(paprs(i,k)-paprs(i,k+1))/rg
3351	ENDDO
3352	ENDDO
3353	!
3354	!IM calcul composantes axiales du moment angulaire et couple des montagnes
3355	!
3356	IF (is_sequential .and. ok_orodr) THEN
3357	CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur, &
3358	ra,rg,romega, &
3359	rlat,rlon,pphis, &
3360	zustrdr,zustrli,zustrph, &
3361	zvstrdr,zvstrli,zvstrph, &
3362	paprs,u,v, &
3363	aam, torsfc)
3364	ENDIF
3365	!IM cf. FLott END
3366	!IM
3367	IF (ip_ebil_phy.ge.2) THEN
3368	ztit='after orography'
3369	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
3370	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
3371	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3372	call diagphy(airephy,ztit,ip_ebil_phy &
3373	, zero_v, zero_v, zero_v, zero_v, zero_v &
3374	, zero_v, zero_v, zero_v, ztsol &
3375	, d_h_vcol, d_qt, d_ec &
3376	, fs_bound, fq_bound )
3377	END IF
3378	!
3379	!
3380	!====================================================================
3381	! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
3382	!====================================================================
3383	! Abderrahmane 24.08.09
3384
3385	IF (ok_cosp) THEN
3386	! adeclarer
3387	#ifdef CPP_COSP
3388	IF (MOD(itap,NINT(freq_cosp/dtime)).EQ.0) THEN
3389
3390	print*,'freq_cosp',freq_cosp
3391	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
3392	! print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
3393	! s ref_liq,ref_ice
3394	call phys_cosp(itap,dtime,freq_cosp, &
3395	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
3396	ecrit_mth,ecrit_day,ecrit_hf, &
3397	klon,klev,rlon,rlat,presnivs,overlap, &
3398	fract,ref_liq,ref_ice, &
3399	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
3400	zu10m,zv10m,pphis, &
3401	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
3402	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
3403	prfl(:,1:klev),psfl(:,1:klev), &
3404	pmflxr(:,1:klev),pmflxs(:,1:klev), &
3405	mr_ozone,cldtau, cldemi)
3406
3407	! L calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
3408	! L cfaddbze,clcalipso2,dbze,cltlidarradar,
3409	! M clMISR,
3410	! R clisccp2,boxtauisccp,boxptopisccp,tclisccp,ctpisccp,
3411	! I tauisccp,albisccp,meantbisccp,meantbclrisccp)
3412
3413	ENDIF
3414
3415	#endif
3416	ENDIF !ok_cosp
3417	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3418	!AA
3419	!AA Installation de l'interface online-offline pour traceurs
3420	!AA
3421	!====================================================================
3422	! Calcul des tendances traceurs
3423	!====================================================================
3424	!
3425
3426	IF (type_trac=='repr') THEN
3427	sh_in(:,:) = q_seri(:,:)
3428	ELSE
3429	sh_in(:,:) = qx(:,:,ivap)
3430	END IF
3431
3432	call phytrac ( &
3433	itap, days_elapsed+1, jH_cur, debut, &
3434	lafin, dtime, u, v, t, &
3435	paprs, pplay, pmfu, pmfd, &
3436	pen_u, pde_u, pen_d, pde_d, &
3437	cdragh, coefh(1:klon,1:klev,is_ave), fm_therm, entr_therm, &
3438	u1, v1, ftsol, pctsrf, &
3439	zustar, zu10m, zv10m, &
3440	wstar(:,is_ave), ale_bl, ale_wake, &
3441	rlat, rlon, &
3442	frac_impa,frac_nucl, beta_prec_fisrt,beta_prec, &
3443	presnivs, pphis, pphi, albsol1, &
3444	sh_in, rhcl, cldfra, rneb, &
3445	diafra, cldliq, itop_con, ibas_con, &
3446	pmflxr, pmflxs, prfl, psfl, &
3447	da, phi, mp, upwd, &
3448	phi2, d1a, dam, sij, & !<<RomP
3449	wdtrainA, wdtrainM, sigd, clw,elij, & !<<RomP
3450	ev, ep, epmlmMm, eplaMm, & !<<RomP
3451	dnwd, aerosol_couple, flxmass_w, &
3452	tau_aero, piz_aero, cg_aero, ccm, &
3453	rfname, &
3454	d_tr_dyn, & !<<RomP
3455	tr_seri)
3456
3457	IF (offline) THEN
3458
3459	IF (prt_level.ge.9) &
3460	print*,'Attention on met a 0 les thermiques pour phystoke'
3461	call phystokenc ( &
3462	nlon,klev,pdtphys,rlon,rlat, &
3463	t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
3464	fm_therm,entr_therm, &
3465	cdragh,coefh(1:klon,1:klev,is_ave),u1,v1,ftsol,pctsrf, &
3466	frac_impa, frac_nucl, &
3467	pphis,airephy,dtime,itap, &
3468	qx(:,:,ivap),da,phi,mp,upwd,dnwd)
3469
3470
3471	ENDIF
3472
3473	!
3474	! Calculer le transport de l'eau et de l'energie (diagnostique)
3475	!
3476	CALL transp (paprs,zxtsol, &
3477	t_seri, q_seri, u_seri, v_seri, zphi, &
3478	ve, vq, ue, uq)
3479	!
3480	!IM global posePB BEG
3481	IF(1.EQ.0) THEN
3482	!
3483	CALL transp_lay (paprs,zxtsol, &
3484	t_seri, q_seri, u_seri, v_seri, zphi, &
3485	ve_lay, vq_lay, ue_lay, uq_lay)
3486	!
3487	ENDIF !(1.EQ.0) THEN
3488	!IM global posePB END
3489	! Accumuler les variables a stocker dans les fichiers histoire:
3490	!
3491
3492	!================================================================
3493	! Conversion of kinetic and potential energy into heat, for
3494	! parameterisation of subgrid-scale motions
3495	!================================================================
3496
3497	d_t_ec(:,:)=0.
3498	forall (k=1: llm) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
3499	CALL ener_conserv(klon,klev,pdtphys,u,v,t,qx(:,:,ivap), &
3500	u_seri,v_seri,t_seri,q_seri,pbl_tke(:,:,is_ave)-tke0(:,:), &
3501	zmasse,exner,d_t_ec)
3502	t_seri(:,:)=t_seri(:,:)+d_t_ec(:,:)
3503
3504	!IM
3505	IF (ip_ebil_phy.ge.1) THEN
3506	ztit='after physic'
3507	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime &
3508	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
3509	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3510	! Comme les tendances de la physique sont ajoute dans la dynamique,
3511	! on devrait avoir que la variation d'entalpie par la dynamique
3512	! est egale a la variation de la physique au pas de temps precedent.
3513	! Donc la somme de ces 2 variations devrait etre nulle.
3514
3515	call diagphy(airephy,ztit,ip_ebil_phy &
3516	, topsw, toplw, solsw, sollw, sens &
3517	, evap, rain_fall, snow_fall, ztsol &
3518	, d_h_vcol, d_qt, d_ec &
3519	, fs_bound, fq_bound )
3520	!
3521	d_h_vcol_phy=d_h_vcol
3522	!
3523	END IF
3524	!
3525	!=======================================================================
3526	! SORTIES
3527	!=======================================================================
3528	!
3529	!IM initialisation + calculs divers diag AMIP2
3530	!
3531	include "calcul_divers.h"
3532	!
3533	!IM Interpolation sur les niveaux de pression du NMC
3534	! -------------------------------------------------
3535	!
3536	include "calcul_STDlev.h"
3537	!
3538	! slp sea level pressure
3539	slp(:) = paprs(:,1)exp(pphis(:)/(RDt_seri(:,1)))
3540	!
3541	!cc prw = eau precipitable
3542	DO i = 1, klon
3543	prw(i) = 0.
3544	DO k = 1, klev
3545	prw(i) = prw(i) + &
3546	q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
3547	ENDDO
3548	ENDDO
3549	!
3550	IF (type_trac == 'inca') THEN
3551	#ifdef INCA
3552	CALL VTe(VTphysiq)
3553	CALL VTb(VTinca)
3554
3555	CALL chemhook_end ( &
3556	dtime, &
3557	pplay, &
3558	t_seri, &
3559	tr_seri, &
3560	nbtr, &
3561	paprs, &
3562	q_seri, &
3563	airephy, &
3564	pphi, &
3565	pphis, &
3566	zx_rh)
3567
3568	CALL VTe(VTinca)
3569	CALL VTb(VTphysiq)
3570	#endif
3571	END IF
3572
3573
3574	!
3575	! Convertir les incrementations en tendances
3576	!
3577	IF (prt_level .GE.10) THEN
3578	print *,'Convertir les incrementations en tendances '
3579	ENDIF
3580	!
3581	if (mydebug) then
3582	call writefield_phy('u_seri',u_seri,llm)
3583	call writefield_phy('v_seri',v_seri,llm)
3584	call writefield_phy('t_seri',t_seri,llm)
3585	call writefield_phy('q_seri',q_seri,llm)
3586	endif
3587
3588	DO k = 1, klev
3589	DO i = 1, klon
3590	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
3591	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
3592	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
3593	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
3594	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
3595	ENDDO
3596	ENDDO
3597	!
3598	IF (nqtot.GE.3) THEN
3599	DO iq = 3, nqtot
3600	DO k = 1, klev
3601	DO i = 1, klon
3602	d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
3603	ENDDO
3604	ENDDO
3605	ENDDO
3606	ENDIF
3607	!
3608	!IM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
3609	!IM global posePB include "write_bilKP_ins.h"
3610	!IM global posePB include "write_bilKP_ave.h"
3611	!
3612
3613	! Sauvegarder les valeurs de t et q a la fin de la physique:
3614	!
3615	DO k = 1, klev
3616	DO i = 1, klon
3617	u_ancien(i,k) = u_seri(i,k)
3618	v_ancien(i,k) = v_seri(i,k)
3619	t_ancien(i,k) = t_seri(i,k)
3620	q_ancien(i,k) = q_seri(i,k)
3621	ENDDO
3622	ENDDO
3623
3624	!!! RomP >>>
3625	IF (nqtot.GE.3) THEN
3626	DO iq = 3, nqtot
3627	DO k = 1, klev
3628	DO i = 1, klon
3629	tr_ancien(i,k,iq-2) = tr_seri(i,k,iq-2)
3630	ENDDO
3631	ENDDO
3632	ENDDO
3633	ENDIF
3634	!!! RomP <<<
3635	!==========================================================================
3636	! Sorties des tendances pour un point particulier
3637	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
3638	! pour le debug
3639	! La valeur de igout est attribuee plus haut dans le programme
3640	!==========================================================================
3641
3642	if (prt_level.ge.1) then
3643	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
3644	write(lunout,*) &
3645	'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
3646	write(lunout,*) &
3647	nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys, &
3648	pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce), &
3649	pctsrf(igout,is_sic)
3650	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
3651	do k=1,klev
3652	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k), &
3653	d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k), &
3654	d_t_eva(igout,k)
3655	enddo
3656	write(lunout,*) 'cool,heat'
3657	do k=1,klev
3658	write(lunout,*) cool(igout,k),heat(igout,k)
3659	enddo
3660
3661	write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
3662	do k=1,klev
3663	write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k), &
3664	d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
3665	enddo
3666
3667	write(lunout,*) 'd_ps ',d_ps(igout)
3668	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
3669	do k=1,klev
3670	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k), &
3671	d_qx(igout,k,1),d_qx(igout,k,2)
3672	enddo
3673	endif
3674
3675	!==========================================================================
3676
3677	!============================================================
3678	! Calcul de la temperature potentielle
3679	!============================================================
3680	DO k = 1, klev
3681	DO i = 1, klon
3682	!JYG/IM theta en debut du pas de temps
3683	!JYG/IM theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
3684	!JYG/IM theta en fin de pas de temps de physique
3685	theta(i,k)=t_seri(i,k)(100000./pplay(i,k))*(RD/RCPD)
3686	! thetal: 2 lignes suivantes a decommenter si vous avez les fichiers MPL 20130625
3687	! fth_fonctions.F90 et parkind1.F90
3688	! sinon thetal=theta
3689	! thetal(i,k)=fth_thetal(pplay(i,k),t_seri(i,k),q_seri(i,k),
3690	! : ql_seri(i,k))
3691	thetal(i,k)=theta(i,k)
3692	ENDDO
3693	ENDDO
3694	!
3695
3696	! 22.03.04 BEG
3697	!=============================================================
3698	! Ecriture des sorties
3699	!=============================================================
3700	#ifdef CPP_IOIPSL
3701
3702	! Recupere des varibles calcule dans differents modules
3703	! pour ecriture dans histxxx.nc
3704
3705	! Get some variables from module fonte_neige_mod
3706	CALL fonte_neige_get_vars(pctsrf, &
3707	zxfqcalving, zxfqfonte, zxffonte)
3708
3709
3710
3711
3712	!=============================================================
3713	! Separation entre thermiques et non thermiques dans les sorties
3714	! de fisrtilp
3715	!=============================================================
3716
3717	if (iflag_thermals>=1) then
3718	d_t_lscth=0.
3719	d_t_lscst=0.
3720	d_q_lscth=0.
3721	d_q_lscst=0.
3722	do k=1,klev
3723	do i=1,klon
3724	if (ptconvth(i,k)) then
3725	d_t_lscth(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
3726	d_q_lscth(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
3727	else
3728	d_t_lscst(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
3729	d_q_lscst(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
3730	endif
3731	enddo
3732	enddo
3733
3734	do i=1,klon
3735	plul_st(i)=prfl(i,lmax_th(i)+1)+psfl(i,lmax_th(i)+1)
3736	plul_th(i)=prfl(i,1)+psfl(i,1)
3737	enddo
3738	endif
3739
3740
3741	!On effectue les sorties:
3742
3743	CALL phys_output_write(itap, pdtphys, paprs, pphis, &
3744	pplay, lmax_th, aerosol_couple, &
3745	ok_ade, ok_aie, ivap, new_aod, ok_sync, &
3746	ptconv, read_climoz, clevSTD, &
3747	ptconvth, d_t, qx, d_qx, zmasse, &
3748	flag_aerosol_strat)
3749
3750
3751
3752
3753	include "write_histday_seri.h"
3754
3755	include "write_paramLMDZ_phy.h"
3756
3757	#endif
3758
3759	! 22.03.04 END
3760	!
3761	!====================================================================
3762	! Si c'est la fin, il faut conserver l'etat de redemarrage
3763	!====================================================================
3764	!
3765
3766	! -----------------------------------------------------------------
3767	! WSTATS: Saving statistics
3768	! -----------------------------------------------------------------
3769	! ("stats" stores and accumulates 8 key variables in file "stats.nc"
3770	! which can later be used to make the statistic files of the run:
3771	! "stats") only possible in 3D runs !
3772
3773
3774	IF (callstats) THEN
3775
3776	call wstats(klon,o_psol%name,"Surface pressure","Pa" &
3777	,2,paprs(:,1))
3778	call wstats(klon,o_tsol%name,"Surface temperature","K", &
3779	2,zxtsol)
3780	zx_tmp_fi2d(:) = rain_fall(:) + snow_fall(:)
3781	call wstats(klon,o_precip%name,"Precip Totale liq+sol", &
3782	"kg/(s*m2)",2,zx_tmp_fi2d)
3783	zx_tmp_fi2d(:) = rain_lsc(:) + snow_lsc(:)
3784	call wstats(klon,o_plul%name,"Large-scale Precip", &
3785	"kg/(s*m2)",2,zx_tmp_fi2d)
3786	zx_tmp_fi2d(:) = rain_con(:) + snow_con(:)
3787	call wstats(klon,o_pluc%name,"Convective Precip", &
3788	"kg/(s*m2)",2,zx_tmp_fi2d)
3789	call wstats(klon,o_sols%name,"Solar rad. at surf.", &
3790	"W/m2",2,solsw)
3791	call wstats(klon,o_soll%name,"IR rad. at surf.", &
3792	"W/m2",2,sollw)
3793	zx_tmp_fi2d(:) = topsw(:)-toplw(:)
3794	call wstats(klon,o_nettop%name,"Net dn radiatif flux at TOA", &
3795	"W/m2",2,zx_tmp_fi2d)
3796
3797
3798
3799	call wstats(klon,o_temp%name,"Air temperature","K", &
3800	3,t_seri)
3801	call wstats(klon,o_vitu%name,"Zonal wind","m.s-1", &
3802	3,u_seri)
3803	call wstats(klon,o_vitv%name,"Meridional wind", &
3804	"m.s-1",3,v_seri)
3805	call wstats(klon,o_vitw%name,"Vertical wind", &
3806	"m.s-1",3,omega)
3807	call wstats(klon,o_ovap%name,"Specific humidity", "kg/kg", &
3808	3,q_seri)
3809
3810
3811
3812	IF(lafin) THEN
3813	write (,) "Writing stats..."
3814	call mkstats(ierr)
3815	ENDIF
3816
3817	ENDIF !if callstats
3818
3819	IF (lafin) THEN
3820	itau_phy = itau_phy + itap
3821	CALL phyredem ("restartphy.nc")
3822	! open(97,form="unformatted",file="finbin")
3823	! write(97) u_seri,v_seri,t_seri,q_seri
3824	! close(97)
3825	!$OMP MASTER
3826	if (read_climoz >= 1) then
3827	if (is_mpi_root) then
3828	call nf95_close(ncid_climoz)
3829	end if
3830	deallocate(press_climoz) ! pointer
3831	end if
3832	!$OMP END MASTER
3833	ENDIF
3834
3835	! first=.false.
3836
3837	RETURN
3838	END SUBROUTINE physiq
3839	FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
3840	IMPLICIT none
3841	!
3842	! Calculer et imprimer l'eau totale. A utiliser pour verifier
3843	! la conservation de l'eau
3844	!
3845	include "YOMCST.h"
3846	INTEGER klon,klev
3847	REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
3848	REAL aire(klon)
3849	REAL qtotal, zx, qcheck
3850	INTEGER i, k
3851	!
3852	zx = 0.0
3853	DO i = 1, klon
3854	zx = zx + aire(i)
3855	ENDDO
3856	qtotal = 0.0
3857	DO k = 1, klev
3858	DO i = 1, klon
3859	qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i) &
3860	*(paprs(i,k)-paprs(i,k+1))/RG
3861	ENDDO
3862	ENDDO
3863	!
3864	qcheck = qtotal/zx
3865	!
3866	RETURN
3867	END FUNCTION qcheck
3868	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
3869	IMPLICIT none
3870	!
3871	! Tranformer une variable de la grille physique a
3872	! la grille d'ecriture
3873	!
3874	INTEGER nfield,nlon,iim,jjmp1, jjm
3875	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
3876	!
3877	INTEGER i, n, ig
3878	!
3879	jjm = jjmp1 - 1
3880	DO n = 1, nfield
3881	DO i=1,iim
3882	ecrit(i,n) = fi(1,n)
3883	ecrit(i+jjm*iim,n) = fi(nlon,n)
3884	ENDDO
3885	DO ig = 1, nlon - 2
3886	ecrit(iim+ig,n) = fi(1+ig,n)
3887	ENDDO
3888	ENDDO
3889	RETURN
3890	END SUBROUTINE gr_fi_ecrit

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format