Context Navigation

physiq.F @ 5156

Last change on this file since 5156 was 1791, checked in by Ehouarn Millour, 11 years ago

Enrichissement et organisation en module de la structure ctrl_out des variables de sortie: ajout des champs description, unité et type_ecrit. Adaptation en conséquence des routines histdef et histwrite.
UG
...................................................

Improvement and transformation into a module of the ctrl_out structure describing output vars. New fields: description, unit and type_ecrit. Creation of new routines histwrite and histdef to take advantage of this structure.
UG

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 136.8 KB

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

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

Download in other formats:

Original Format