Context Navigation

physiq.F @ 1114

Last change on this file since 1114 was 1114, checked in by jghattas, 15 years ago

Creation du module infotrac:

contient les variables de advtrac.h
contient la subroutine iniadvtrac renommer en infotrac_init
le nombre des traceurs est lu dans tracer.def en dynamique (ou par default ou recu par INCA)
ce module est utilise dans la dynamique et la physique
contient aussi la variable nbtr qui avant etait stockee dans dimphy

Le fichier advtrac.h n'existe plus.
La compilation ne prend plus en compte le nombre de traceur.

/JG

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

File size: 114.7 KB

Line
1	c
2	c#define IO_DEBUG
3
4	SUBROUTINE physiq (nlon,nlev,
5	. debut,lafin,rjourvrai,gmtime,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
14	USE comgeomphy
15	USE write_field_phy
16	USE dimphy
17	USE infotrac
18	USE mod_grid_phy_lmdz
19	USE mod_phys_lmdz_para
20	USE iophy
21	USE misc_mod, mydebug=>debug
22	USE vampir
23	USE pbl_surface_mod, ONLY : pbl_surface
24	USE change_srf_frac_mod
25	USE surface_data, ONLY : type_ocean, ok_veget
26	USE phys_local_var_mod ! Variables internes non sauvegardees de la physique
27	USE phys_state_var_mod ! Variables sauvegardees de la physique
28	USE fonte_neige_mod, ONLY : fonte_neige_get_vars
29	USE phys_output_mod
30
31	IMPLICIT none
32	c======================================================================
33	c
34	c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
35	c
36	c Objet: Moniteur general de la physique du modele
37	cAA Modifications quant aux traceurs :
38	cAA - uniformisation des parametrisations ds phytrac
39	cAA - stockage des moyennes des champs necessaires
40	cAA en mode traceur off-line
41	c======================================================================
42	c CLEFS CPP POUR LES IO
43	c =====================
44	c#define histmthNMC
45	c#define histISCCP
46	c======================================================================
47	c modif ( P. Le Van , 12/10/98 )
48	c
49	c Arguments:
50	c
51	c nlon----input-I-nombre de points horizontaux
52	c nlev----input-I-nombre de couches verticales
53	c debut---input-L-variable logique indiquant le premier passage
54	c lafin---input-L-variable logique indiquant le dernier passage
55	c rjour---input-R-numero du jour de l'experience
56	c gmtime--input-R-temps universel dans la journee (0 a 86400 s)
57	c pdtphys-input-R-pas d'integration pour la physique (seconde)
58	c paprs---input-R-pression pour chaque inter-couche (en Pa)
59	c pplay---input-R-pression pour le mileu de chaque couche (en Pa)
60	c pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
61	c pphis---input-R-geopotentiel du sol
62	c presnivs-input_R_pressions approximat. des milieux couches ( en PA)
63	c u-------input-R-vitesse dans la direction X (de O a E) en m/s
64	c v-------input-R-vitesse Y (de S a N) en m/s
65	c t-------input-R-temperature (K)
66	c qx------input-R-humidite specifique (kg/kg) et d'autres traceurs
67	c d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
68	c d_q_dyn-input-R-tendance dynamique pour "q" (kg/kg/s)
69	c flxmass_w -input-R- flux de masse verticale
70	c d_u-----output-R-tendance physique de "u" (m/s/s)
71	c d_v-----output-R-tendance physique de "v" (m/s/s)
72	c d_t-----output-R-tendance physique de "t" (K/s)
73	c d_qx----output-R-tendance physique de "qx" (kg/kg/s)
74	c d_ps----output-R-tendance physique de la pression au sol
75	cIM
76	c PVteta--output-R-vorticite potentielle a des thetas constantes
77	c======================================================================
78	#include "dimensions.h"
79	integer jjmp1
80	parameter (jjmp1=jjm+1-1/jjm)
81	integer iip1
82	parameter (iip1=iim+1)
83
84	#include "regdim.h"
85	#include "indicesol.h"
86	#include "dimsoil.h"
87	#include "clesphys.h"
88	#include "control.h"
89	#include "temps.h"
90	#include "iniprint.h"
91	#include "thermcell.h"
92	c======================================================================
93	LOGICAL ok_cvl ! pour activer le nouveau driver pour convection KE
94	PARAMETER (ok_cvl=.TRUE.)
95	LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface
96	PARAMETER (ok_gust=.FALSE.)
97	integer iflag_radia ! active ou non le rayonnement (MPL)
98	save iflag_radia
99	c$OMP THREADPRIVATE(iflag_radia)
100	c======================================================================
101	LOGICAL check ! Verifier la conservation du modele en eau
102	PARAMETER (check=.FALSE.)
103	LOGICAL ok_stratus ! Ajouter artificiellement les stratus
104	PARAMETER (ok_stratus=.FALSE.)
105	c======================================================================
106	LOGICAL, SAVE :: rnpb=.TRUE.
107	c$OMP THREADPRIVATE(rnpb)
108	REAL amn, amx
109	INTEGER igout
110	c======================================================================
111	c Clef controlant l'activation du cycle diurne:
112	ccc LOGICAL cycle_diurne
113	ccc PARAMETER (cycle_diurne=.FALSE.)
114	c======================================================================
115	c Modele thermique du sol, a activer pour le cycle diurne:
116	ccc LOGICAL soil_model
117	ccc PARAMETER (soil_model=.FALSE.)
118	c======================================================================
119	c Dans les versions precedentes, l'eau liquide nuageuse utilisee dans
120	c le calcul du rayonnement est celle apres la precipitation des nuages.
121	c Si cette cle new_oliq est activee, ce sera une valeur moyenne entre
122	c la condensation et la precipitation. Cette cle augmente les impacts
123	c radiatifs des nuages.
124	ccc LOGICAL new_oliq
125	ccc PARAMETER (new_oliq=.FALSE.)
126	c======================================================================
127	c Clefs controlant deux parametrisations de l'orographie:
128	cc LOGICAL ok_orodr
129	ccc PARAMETER (ok_orodr=.FALSE.)
130	ccc LOGICAL ok_orolf
131	ccc PARAMETER (ok_orolf=.FALSE.)
132	c======================================================================
133	LOGICAL ok_journe ! sortir le fichier journalier
134	save ok_journe
135	c$OMP THREADPRIVATE(ok_journe)
136	c
137	LOGICAL ok_mensuel ! sortir le fichier mensuel
138	save ok_mensuel
139	c$OMP THREADPRIVATE(ok_mensuel)
140	c
141	LOGICAL ok_instan ! sortir le fichier instantane
142	save ok_instan
143	c$OMP THREADPRIVATE(ok_instan)
144	c
145	LOGICAL ok_LES ! sortir le fichier LES
146	save ok_LES
147	c$OMP THREADPRIVATE(ok_LES)
148	c
149	LOGICAL ok_region ! sortir le fichier regional
150	PARAMETER (ok_region=.FALSE.)
151	c======================================================================
152	real weak_inversion(klon),dthmin(klon)
153	real seuil_inversion
154	save seuil_inversion
155	c$OMP THREADPRIVATE(seuil_inversion)
156	integer iflag_ratqs
157	save iflag_ratqs
158	c$OMP THREADPRIVATE(iflag_ratqs)
159	REAL lambda_th(klon,klev),zz,znum,zden
160	REAL wmax_th(klon)
161	REAL zmax_th(klon)
162	REAL tau_overturning_th(klon)
163
164	integer lmax_th(klon)
165	integer limbas(klon)
166	real ratqscth(klon,klev)
167	real ratqsdiff(klon,klev)
168	real zqsatth(klon,klev)
169
170	c======================================================================
171	c
172	INTEGER ivap ! indice de traceurs pour vapeur d'eau
173	PARAMETER (ivap=1)
174	INTEGER iliq ! indice de traceurs pour eau liquide
175	PARAMETER (iliq=2)
176
177	c
178	c
179	c Variables argument:
180	c
181	INTEGER nlon
182	INTEGER nlev
183	REAL rjourvrai
184	REAL gmtime
185	REAL pdtphys
186	LOGICAL debut, lafin
187	REAL paprs(klon,klev+1)
188	REAL pplay(klon,klev)
189	REAL pphi(klon,klev)
190	REAL pphis(klon)
191	REAL presnivs(klev)
192	REAL znivsig(klev)
193	real pir
194
195	REAL u(klon,klev)
196	REAL v(klon,klev)
197	REAL t(klon,klev),theta(klon,klev)
198	REAL qx(klon,klev,nqtot)
199	REAL flxmass_w(klon,klev)
200	REAL omega(klon,klev) ! vitesse verticale en Pa/s
201	REAL d_u(klon,klev)
202	REAL d_v(klon,klev)
203	REAL d_t(klon,klev)
204	REAL d_qx(klon,klev,nqtot)
205	REAL d_ps(klon)
206	real da(klon,klev),phi(klon,klev,klev),mp(klon,klev)
207	c
208	cIM Amip2 PV a theta constante
209	c
210	INTEGER nbteta
211	PARAMETER(nbteta=3)
212	CHARACTER*3 ctetaSTD(nbteta)
213	DATA ctetaSTD/'350','380','405'/
214	SAVE ctetaSTD
215	c$OMP THREADPRIVATE(ctetaSTD)
216	REAL rtetaSTD(nbteta)
217	DATA rtetaSTD/350., 380., 405./
218	SAVE rtetaSTD
219	c$OMP THREADPRIVATE(rtetaSTD)
220	c
221	REAL PVteta(klon,nbteta)
222	REAL zx_tmp_3dte(iim,jjmp1,nbteta)
223	c
224	cMI Amip2 PV a theta constante
225
226	cym INTEGER klevp1, klevm1
227	cym PARAMETER(klevp1=klev+1,klevm1=klev-1)
228	cym#include "raddim.h"
229	c
230	c
231	cIM Amip2
232	c variables a une pression donnee
233	c
234	real rlevSTD(nlevSTD)
235	DATA rlevSTD/100000., 92500., 85000., 70000.,
236	.60000., 50000., 40000., 30000., 25000., 20000.,
237	.15000., 10000., 7000., 5000., 3000., 2000., 1000./
238	SAVE rlevstd
239	c$OMP THREADPRIVATE(rlevstd)
240	CHARACTER*4 clevSTD(nlevSTD)
241	DATA clevSTD/'1000','925 ','850 ','700 ','600 ',
242	.'500 ','400 ','300 ','250 ','200 ','150 ','100 ',
243	.'70 ','50 ','30 ','20 ','10 '/
244	SAVE clevSTD
245	c$OMP THREADPRIVATE(clevSTD)
246	c
247	CHARACTER*4 bb2
248	CHARACTER*2 bb3
249	c
250	real tlevSTD(klon,nlevSTD), qlevSTD(klon,nlevSTD)
251	real rhlevSTD(klon,nlevSTD), philevSTD(klon,nlevSTD)
252	real ulevSTD(klon,nlevSTD), vlevSTD(klon,nlevSTD)
253	real wlevSTD(klon,nlevSTD)
254
255	real twriteSTD(klon,nlevSTD,nfiles)
256	real qwriteSTD(klon,nlevSTD,nfiles)
257	real rhwriteSTD(klon,nlevSTD,nfiles)
258	real phiwriteSTD(klon,nlevSTD,nfiles)
259	real uwriteSTD(klon,nlevSTD,nfiles)
260	real vwriteSTD(klon,nlevSTD,nfiles)
261	real wwriteSTD(klon,nlevSTD,nfiles)
262	c
263	c nout : niveau de output des variables a une pression donnee
264	logical oknondef(klon,nlevSTD,nout)
265	c
266	c les produits uvSTD, vqSTD, .., T2STD sont calcules
267	c a partir des valeurs instantannees toutes les 6 h
268	c qui sont moyennees sur le mois
269	c
270	real uvSTD(klon,nlevSTD)
271	real vqSTD(klon,nlevSTD)
272	real vTSTD(klon,nlevSTD)
273	real wqSTD(klon,nlevSTD)
274	c
275	real vphiSTD(klon,nlevSTD)
276	real wTSTD(klon,nlevSTD)
277	real u2STD(klon,nlevSTD)
278	real v2STD(klon,nlevSTD)
279	real T2STD(klon,nlevSTD)
280	c
281	#include "radepsi.h"
282	#include "radopt.h"
283	c
284	c
285	c prw: precipitable water
286	real prw(klon)
287
288	REAL convliq(klon,klev) ! eau liquide nuageuse convective
289	REAL convfra(klon,klev) ! fraction nuageuse convective
290
291	REAL cldl_c(klon),cldm_c(klon),cldh_c(klon) !nuages bas, moyen et haut
292	REAL cldt_c(klon),cldq_c(klon) !nuage total, eau liquide integree
293	REAL cldl_s(klon),cldm_s(klon),cldh_s(klon) !nuages bas, moyen et haut
294	REAL cldt_s(klon),cldq_s(klon) !nuage total, eau liquide integree
295
296	INTEGER linv, kp1
297	c flwp, fiwp = Liquid Water Path & Ice Water Path (kg/m2)
298	c flwc, fiwc = Liquid Water Content & Ice Water Content (kg/kg)
299	REAL flwp(klon), fiwp(klon)
300	REAL flwc(klon,klev), fiwc(klon,klev)
301	REAL flwp_c(klon), fiwp_c(klon)
302	REAL flwc_c(klon,klev), fiwc_c(klon,klev)
303	REAL flwp_s(klon), fiwp_s(klon)
304	REAL flwc_s(klon,klev), fiwc_s(klon,klev)
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 *27 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 quasi-arguments
526	c
527	REAL xjour
528	SAVE xjour
529	c$OMP THREADPRIVATE(xjour)
530	c
531	c
532	c Variables propres a la physique
533	c
534	c INTEGER radpas
535	c SAVE radpas ! frequence d'appel rayonnement
536	ccccccccc$OMP THREADPRIVATE(radpas)
537	c
538	cc INTEGER iflag_con
539	c
540	INTEGER itap
541	SAVE itap ! compteur pour la physique
542	c$OMP THREADPRIVATE(itap)
543	c
544	real slp(klon) ! sea level pressure
545	c
546	REAL fevap(klon,nbsrf)
547	REAL fluxlat(klon,nbsrf)
548	c
549	REAL qsol(klon)
550	REAL,save :: solarlong0
551	c$OMP THREADPRIVATE(solarlong0)
552
553	c
554	c Parametres de l'Orographie a l'Echelle Sous-Maille (OESM):
555	c
556	cIM 141004 REAL zulow(klon),zvlow(klon),zustr(klon), zvstr(klon)
557	REAL zulow(klon),zvlow(klon)
558	c
559	INTEGER igwd,idx(klon),itest(klon)
560	c
561	REAL agesno(klon,nbsrf)
562	c
563	c REAL,allocatable,save :: run_off_lic_0(:)
564	cc$OMP THREADPRIVATE(run_off_lic_0)
565	cym SAVE run_off_lic_0
566	cKE43
567	c Variables liees a la convection de K. Emanuel (sb):
568	c
569	REAL bas, top ! cloud base and top levels
570	SAVE bas
571	SAVE top
572	c$OMP THREADPRIVATE(bas, top)
573
574	REAL wdn(klon), tdn(klon), qdn(klon)
575	c
576	c=================================================================================================
577	cCR04.12.07: on ajoute les nouvelles variables du nouveau schema de convection avec poches froides
578	c Variables liées à la poche froide (jyg)
579
580	REAL mip(klon,klev) ! mass flux shed by the adiab ascent at each level
581	REAL Vprecip(klon,klev) ! precipitation vertical profile
582	c
583	REAL wape_prescr, fip_prescr
584	INTEGER it_wape_prescr
585	SAVE wape_prescr, fip_prescr, it_wape_prescr
586	c$OMP THREADPRIVATE(wape_prescr, fip_prescr, it_wape_prescr)
587	c
588	c variables supplementaires de concvl
589	REAL Tconv(klon,klev)
590	REAL ment(klon,klev,klev),sij(klon,klev,klev)
591	REAL dd_t(klon,klev),dd_q(klon,klev)
592
593	real, save :: alp_bl_prescr=0.
594	real, save :: ale_bl_prescr=0.
595
596	real, save :: ale_max=100.
597	real, save :: alp_max=2.
598
599	c$OMP THREADPRIVATE(alp_bl_prescr,ale_bl_prescr)
600	c$OMP THREADPRIVATE(ale_max,alp_max)
601
602	real ale_wake(klon)
603	real alp_wake(klon)
604	cRC
605	c Variables liées à la poche froide (jyg et rr)
606	c Version diagnostique pour l'instant : pas de rétroaction sur la convection
607
608	REAL t_wake(klon,klev),q_wake(klon,klev) ! wake pour la convection
609
610	REAL wake_dth(klon,klev) ! wake : temp pot difference
611
612	REAL wake_d_deltat_gw(klon,klev)! wake : delta T tendency due to Gravity Wave (/s)
613	REAL wake_omgbdth(klon,klev) ! Wake : flux of Delta_Theta transported by LS omega
614	REAL wake_dp_omgb(klon,klev) ! Wake : vertical gradient of large scale omega
615	REAL wake_dtKE(klon,klev) ! Wake : differential heating (wake - unpertubed) CONV
616	REAL wake_dqKE(klon,klev) ! Wake : differential moistening (wake - unpertubed) CONV
617	REAL wake_dtPBL(klon,klev) ! Wake : differential heating (wake - unpertubed) PBL
618	REAL wake_dqPBL(klon,klev) ! Wake : differential moistening (wake - unpertubed) PBL
619	REAL wake_omg(klon,klev) ! Wake : velocity difference (wake - unpertubed)
620	REAL wake_ddeltat(klon,klev),wake_ddeltaq(klon,klev)
621	REAL wake_dp_deltomg(klon,klev) ! Wake : gradient vertical de wake_omg
622	REAL wake_spread(klon,klev) ! spreading term in wake_delt
623	c
624	cpourquoi y'a pas de save??
625	REAL wake_h(klon) ! Wake : hauteur de la poche froide
626	c
627	INTEGER wake_k(klon) ! Wake sommet
628	c
629	REAL t_undi(klon,klev) ! temperature moyenne dans la zone non perturbee
630	REAL q_undi(klon,klev) ! humidite moyenne dans la zone non perturbee
631	c
632	REAL wake_pe(klon) ! Wake potential energy - WAPE
633
634	REAL wake_gfl(klon) ! Gust Front Length
635	REAL wake_dens(klon)
636	c
637	c
638	REAL dt_dwn(klon,klev)
639	REAL dq_dwn(klon,klev)
640	REAL wdt_PBL(klon,klev)
641	REAL udt_PBL(klon,klev)
642	REAL wdq_PBL(klon,klev)
643	REAL udq_PBL(klon,klev)
644	REAL M_dwn(klon,klev)
645	REAL M_up(klon,klev)
646	REAL dt_a(klon,klev)
647	REAL dq_a(klon,klev)
648	c
649	cRR:fin declarations poches froides
650	c=======================================================================================================
651
652	REAL zw2(klon,klev+1)
653	REAL fraca(klon,klev+1)
654
655	c Variables locales pour la couche limite (al1):
656	c
657	cAl1 REAL pblh(klon) ! Hauteur de couche limite
658	cAl1 SAVE pblh
659	c34EK
660	c
661	c Variables locales:
662	c
663	REAL cdragh(klon) ! drag coefficient pour T and Q
664	REAL cdragm(klon) ! drag coefficient pour vent
665	cAA
666	cAA Pour phytrac
667	cAA
668	REAL coefh(klon,klev) ! coef d'echange pour phytrac, valable pour 2<=k<=klev
669	REAL u1(klon) ! vents dans la premiere couche U
670	REAL v1(klon) ! vents dans la premiere couche V
671
672	REAL zxffonte(klon), zxfqcalving(klon),zxfqfonte(klon)
673
674	c@$$ LOGICAL offline ! Controle du stockage ds "physique"
675	c@$$ PARAMETER (offline=.false.)
676	c@$$ INTEGER physid
677	REAL frac_impa(klon,klev) ! fractions d'aerosols lessivees (impaction)
678	REAL frac_nucl(klon,klev) ! idem (nucleation)
679	INTEGER :: iii
680	REAL :: calday
681
682	cIM cf FH pour Tiedtke 080604
683	REAL rain_tiedtke(klon),snow_tiedtke(klon)
684	c
685	cIM 050204 END
686	REAL evap(klon), devap(klon) ! evaporation et sa derivee
687	REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
688
689	REAL bils(klon) ! bilan de chaleur au sol
690	REAL wfbilo(klon,nbsrf) ! bilan d'eau, pour chaque
691	C ! type de sous-surface et pondere par la fraction
692	REAL wfbils(klon,nbsrf) ! bilan de chaleur au sol, pour chaque
693	C ! type de sous-surface et pondere par la fraction
694	REAL slab_wfbils(klon) ! bilan de chaleur au sol pour le cas de slab, sur les points d'ocean
695
696	REAL fder(klon)
697	REAL ve(klon) ! integr. verticale du transport meri. de l'energie
698	REAL vq(klon) ! integr. verticale du transport meri. de l'eau
699	REAL ue(klon) ! integr. verticale du transport zonal de l'energie
700	REAL uq(klon) ! integr. verticale du transport zonal de l'eau
701	c
702	REAL frugs(klon,nbsrf)
703	REAL zxrugs(klon) ! longueur de rugosite
704	c
705	c Conditions aux limites
706	c
707	INTEGER julien
708	c
709	INTEGER lmt_pas
710	SAVE lmt_pas ! frequence de mise a jour
711	c$OMP THREADPRIVATE(lmt_pas)
712
713	cIM sorties
714	REAL un_jour
715	PARAMETER(un_jour=86400.)
716	c======================================================================
717	c
718	c Declaration des procedures appelees
719	c
720	EXTERNAL angle ! calculer angle zenithal du soleil
721	EXTERNAL alboc ! calculer l'albedo sur ocean
722	EXTERNAL ajsec ! ajustement sec
723	EXTERNAL conlmd ! convection (schema LMD)
724	cKE43
725	EXTERNAL conema3 ! convect4.3
726	EXTERNAL fisrtilp ! schema de condensation a grande echelle (pluie)
727	cAA
728	EXTERNAL fisrtilp_tr ! schema de condensation a grande echelle (pluie)
729	c ! stockage des coefficients necessaires au
730	c ! lessivage OFF-LINE et ON-LINE
731	EXTERNAL hgardfou ! verifier les temperatures
732	EXTERNAL nuage ! calculer les proprietes radiatives
733	EXTERNAL o3cm ! initialiser l'ozone
734	EXTERNAL orbite ! calculer l'orbite terrestre
735	EXTERNAL ozonecm ! prescrire l'ozone
736	EXTERNAL phyetat0 ! lire l'etat initial de la physique
737	EXTERNAL phyredem ! ecrire l'etat de redemarrage de la physique
738	EXTERNAL radlwsw ! rayonnements solaire et infrarouge
739	EXTERNAL suphel ! initialiser certaines constantes
740	EXTERNAL transp ! transport total de l'eau et de l'energie
741	EXTERNAL ecribina ! ecrire le fichier binaire global
742	EXTERNAL ecribins ! ecrire le fichier binaire global
743	EXTERNAL ecrirega ! ecrire le fichier binaire regional
744	EXTERNAL ecriregs ! ecrire le fichier binaire regional
745	cIM
746	EXTERNAL haut2bas !variables de haut en bas
747	INTEGER lnblnk1
748	EXTERNAL lnblnk1 !enleve les blancs a la fin d'une variable de type
749	!caracter
750	EXTERNAL ini_undefSTD !initialise a 0 une variable a 1 niveau de pression
751	EXTERNAL undefSTD !somme les valeurs definies d'1 var a 1 niveau de pression
752	c EXTERNAL moy_undefSTD !moyenne d'1 var a 1 niveau de pression
753	c EXTERNAL moyglo_aire !moyenne globale d'1 var ponderee par l'aire de la maille (moyglo_pondaire)
754	c !par la masse/airetot (moyglo_pondaima) et la vraie masse (moyglo_pondmass)
755	c
756	c Variables locales
757	c
758	REAL rhcl(klon,klev) ! humiditi relative ciel clair
759	REAL dialiq(klon,klev) ! eau liquide nuageuse
760	REAL diafra(klon,klev) ! fraction nuageuse
761	REAL cldliq(klon,klev) ! eau liquide nuageuse
762	REAL cldfra(klon,klev) ! fraction nuageuse
763	REAL cldtau(klon,klev) ! epaisseur optique
764	REAL cldemi(klon,klev) ! emissivite infrarouge
765	c
766	CXXX PB
767	REAL fluxq(klon,klev, nbsrf) ! flux turbulent d'humidite
768	REAL fluxt(klon,klev, nbsrf) ! flux turbulent de chaleur
769	REAL fluxu(klon,klev, nbsrf) ! flux turbulent de vitesse u
770	REAL fluxv(klon,klev, nbsrf) ! flux turbulent de vitesse v
771	c
772	REAL zxfluxt(klon, klev)
773	REAL zxfluxq(klon, klev)
774	REAL zxfluxu(klon, klev)
775	REAL zxfluxv(klon, klev)
776	CXXX
777	c
778	REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous surface
779	REAL fsolsw(klon, nbsrf) ! flux solaire absorb. pour chaque sous surface
780	c Le rayonnement n'est pas calcule tous les pas, il faut donc
781	c sauvegarder les sorties du rayonnement
782	cym SAVE heat,cool,albpla,topsw,toplw,solsw,sollw,sollwdown
783	cym SAVE sollwdownclr, toplwdown, toplwdownclr
784	cym SAVE topsw0,toplw0,solsw0,sollw0, heat0, cool0
785	c
786	INTEGER itaprad
787	SAVE itaprad
788	c$OMP THREADPRIVATE(itaprad)
789	c
790	REAL conv_q(klon,klev) ! convergence de l'humidite (kg/kg/s)
791	REAL conv_t(klon,klev) ! convergence de la temperature(K/s)
792	c
793	REAL cldl(klon),cldm(klon),cldh(klon) !nuages bas, moyen et haut
794	REAL cldt(klon),cldq(klon) !nuage total, eau liquide integree
795	c
796	REAL zxtsol(klon), zxqsurf(klon), zxsnow(klon), zxfluxlat(klon)
797	REAL zxsnow_dummy(klon)
798	c
799	REAL dist, rmu0(klon), fract(klon)
800	REAL zdtime, zlongi
801	c
802	CHARACTER*2 str2
803	CHARACTER*2 iqn
804	c
805	REAL qcheck
806	REAL z_avant(klon), z_apres(klon), z_factor(klon)
807	LOGICAL zx_ajustq
808	c
809	REAL za, zb
810	REAL zx_t, zx_qs, zdelta, zcor, zfra, zlvdcp, zlsdcp
811	real zqsat(klon,klev)
812	INTEGER i, k, iq, ig, j, nsrf, ll, l, iiq, iff
813	REAL t_coup
814	PARAMETER (t_coup=234.0)
815	c
816	REAL zphi(klon,klev)
817	cym A voir plus tard !!
818	cym REAL zx_relief(iim,jjmp1)
819	cym REAL zx_aire(iim,jjmp1)
820	c
821	c Grandeurs de sorties
822	REAL s_pblh(klon), s_lcl(klon), s_capCL(klon)
823	REAL s_oliqCL(klon), s_cteiCL(klon), s_pblt(klon)
824	REAL s_therm(klon), s_trmb1(klon), s_trmb2(klon)
825	REAL s_trmb3(klon)
826	cKE43
827	c Variables locales pour la convection de K. Emanuel (sb):
828	c
829	REAL upwd(klon,klev) ! saturated updraft mass flux
830	REAL dnwd(klon,klev) ! saturated downdraft mass flux
831	REAL dnwd0(klon,klev) ! unsaturated downdraft mass flux
832	REAL tvp(klon,klev) ! virtual temp of lifted parcel
833	CHARACTER*40 capemaxcels !max(CAPE)
834
835	REAL rflag(klon) ! flag fonctionnement de convect
836	INTEGER iflagctrl(klon) ! flag fonctionnement de convect
837	c -- convect43:
838	INTEGER ntra ! nb traceurs pour convect4.3
839	REAL pori_con(klon) ! pressure at the origin level of lifted parcel
840	REAL plcl_con(klon),dtma_con(klon),dtlcl_con(klon)
841	REAL dtvpdt1(klon,klev), dtvpdq1(klon,klev)
842	REAL dplcldt(klon), dplcldr(klon)
843	c? . condm_con(klon,klev),conda_con(klon,klev),
844	c? . mr_con(klon,klev),ep_con(klon,klev)
845	c? . ,sadiab(klon,klev),wadiab(klon,klev)
846	c --
847	c34EK
848	c
849	c Variables du changement
850	c
851	c con: convection
852	c lsc: condensation a grande echelle (Large-Scale-Condensation)
853	c ajs: ajustement sec
854	c eva: evaporation de l'eau liquide nuageuse
855	c vdf: couche limite (Vertical DiFfusion)
856	REAL rneb(klon,klev)
857
858	! tendance nulles
859	REAL du0(klon,klev),dv0(klon,klev),dq0(klon,klev),dql0(klon,klev)
860
861	c
862	*********************************************************
863	* declarations
864
865	*********************************************************
866	cIM 081204 END
867	c
868	REAL pmfu(klon,klev), pmfd(klon,klev)
869	REAL pen_u(klon,klev), pen_d(klon,klev)
870	REAL pde_u(klon,klev), pde_d(klon,klev)
871	INTEGER kcbot(klon), kctop(klon), kdtop(klon)
872	REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1)
873	REAL prfl(klon,klev+1), psfl(klon,klev+1)
874	c
875	REAL rain_lsc(klon)
876	REAL snow_lsc(klon)
877	c
878	REAL ratqss(klon,klev),ratqsc(klon,klev)
879	real ratqsbas,ratqshaut
880	save ratqsbas,ratqshaut
881	c$OMP THREADPRIVATE(ratqsbas,ratqshaut)
882	real zpt_conv(klon,klev)
883
884	c Parametres lies au nouveau schema de nuages (SB, PDF)
885	real fact_cldcon
886	real facttemps
887	logical ok_newmicro
888	save ok_newmicro
889	c$OMP THREADPRIVATE(ok_newmicro)
890	save fact_cldcon,facttemps
891	c$OMP THREADPRIVATE(fact_cldcon,facttemps)
892	real facteur
893
894	integer iflag_cldcon
895	save iflag_cldcon
896	c$OMP THREADPRIVATE(iflag_cldcon)
897	logical ptconv(klon,klev)
898	cIM cf. AM 081204 BEG
899	logical ptconvth(klon,klev)
900	cIM cf. AM 081204 END
901	c
902	c Variables liees a l'ecriture de la bande histoire physique
903	c
904	c======================================================================
905	c
906	cIM cf. AM 081204 BEG
907	c declarations pour sortir sur une sous-region
908	integer imin_ins,imax_ins,jmin_ins,jmax_ins
909	save imin_ins,imax_ins,jmin_ins,jmax_ins
910	c$OMP THREADPRIVATE(imin_ins,imax_ins,jmin_ins,jmax_ins)
911	c real lonmin_ins,lonmax_ins,latmin_ins
912	c s ,latmax_ins
913	c data lonmin_ins,lonmax_ins,latmin_ins
914	c s ,latmax_ins/
915	c valeurs initiales s -5.,20.,41.,55./
916	c s 100.,130.,-20.,20./
917	c s -180.,180.,-90.,90./
918	c======================================================================
919	cIM cf. AM 081204 END
920
921	c
922	integer itau_w ! pas de temps ecriture = itap + itau_phy
923	c
924	c
925	c Variables locales pour effectuer les appels en serie
926	c
927	REAL zx_rh(klon,klev)
928	cIM RH a 2m (la surface)
929	REAL rh2m(klon), qsat2m(klon)
930	REAL tpot(klon), tpote(klon)
931	REAL Lheat
932
933	INTEGER length
934	PARAMETER ( length = 100 )
935	REAL tabcntr0( length )
936	c
937	INTEGER ndex2d(iimjjmp1),ndex3d(iimjjmp1*klev)
938	cIM
939	INTEGER ndex2d1(iwmax)
940	c
941	cIM AMIP2 BEG
942	REAL moyglo, mountor
943	cIM 141004 BEG
944	REAL zustrdr(klon), zvstrdr(klon)
945	REAL zustrli(klon), zvstrli(klon)
946	REAL zustrph(klon), zvstrph(klon)
947	REAL zustrhi(klon), zvstrhi(klon)
948	REAL aam, torsfc
949	cIM 141004 END
950	cIM 190504 BEG
951	INTEGER ij, imp1jmp1
952	PARAMETER(imp1jmp1=(iim+1)*jjmp1)
953	cym A voir plus tard
954	REAL zx_tmp(imp1jmp1), airedyn(iim+1,jjmp1)
955	REAL padyn(iim+1,jjmp1,klev+1)
956	REAL dudyn(iim+1,jjmp1,klev)
957	REAL rlatdyn(iim+1,jjmp1)
958	cIM 190504 END
959	LOGICAL ok_msk
960	REAL msk(klon)
961	cIM
962	REAL airetot, pi
963	cym A voir plus tard
964	cym REAL zm_wo(jjmp1, klev)
965	cIM AMIP2 END
966	c
967	REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique
968	REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D
969	c#ifdef histmthNMC
970	cym A voir plus tard !!!!
971	cym REAL zx_tmp_NC(iim,jjmp1,nlevSTD)
972	REAL zx_tmp_fiNC(klon,nlevSTD)
973	c#endif
974	REAL*8 zx_tmp2_fi3d(klon,klev) ! variable temporaire pour champs 3D
975	REAL zx_tmp_2d(iim,jjmp1), zx_tmp_3d(iim,jjmp1,klev)
976	REAL zx_lon(iim,jjmp1), zx_lat(iim,jjmp1)
977	c
978	INTEGER nid_day, nid_mth, nid_ins, nid_nmc, nid_day_seri
979	INTEGER nid_ctesGCM
980	SAVE nid_day, nid_mth, nid_ins, nid_nmc, nid_day_seri
981	SAVE nid_ctesGCM
982	c$OMP THREADPRIVATE(nid_day, nid_mth, nid_ins, nid_nmc)
983	c$OMP THREADPRIVATE(nid_day_seri,nid_ctesGCM)
984	c
985	cIM 280405 BEG
986	INTEGER nid_bilKPins, nid_bilKPave
987	SAVE nid_bilKPins, nid_bilKPave
988	c$OMP THREADPRIVATE(nid_bilKPins, nid_bilKPave)
989	c
990	REAL ve_lay(klon,klev) ! transport meri. de l'energie a chaque niveau vert.
991	REAL vq_lay(klon,klev) ! transport meri. de l'eau a chaque niveau vert.
992	REAL ue_lay(klon,klev) ! transport zonal de l'energie a chaque niveau vert.
993	REAL uq_lay(klon,klev) ! transport zonal de l'eau a chaque niveau vert.
994	c
995	cIM 280405 END
996	c
997	INTEGER nhori, nvert, nvert1, nvert3
998	REAL zsto, zsto1, zsto2
999	REAL zstophy, zstorad, zstohf, zstoday, zstomth, zout
1000	REAL zcals(napisccp), zcalh(napisccp), zoutj(napisccp)
1001	REAL zout_isccp(napisccp)
1002	SAVE zcals, zcalh, zoutj, zout_isccp
1003	c$OMP THREADPRIVATE(zcals, zcalh, zoutj, zout_isccp)
1004
1005	real zjulian
1006	save zjulian
1007	c$OMP THREADPRIVATE(zjulian)
1008
1009	character*20 modname
1010	character*80 abort_message
1011	logical ok_sync
1012	real date0
1013	integer idayref
1014
1015	C essai writephys
1016	integer fid_day, fid_mth, fid_ins
1017	parameter (fid_ins = 1, fid_day = 2, fid_mth = 3)
1018	integer prof2d_on, prof3d_on, prof2d_av, prof3d_av
1019	parameter (prof2d_on = 1, prof3d_on = 2,
1020	. prof2d_av = 3, prof3d_av = 4)
1021	character*30 nom_fichier
1022	character*10 varname
1023	character*40 vartitle
1024	character*20 varunits
1025	C Variables liees au bilan d'energie et d'enthalpi
1026	REAL ztsol(klon)
1027	REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
1028	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
1029	SAVE h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
1030	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
1031	c$OMP THREADPRIVATE(h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot,
1032	c$OMP+ h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot)
1033	REAL d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
1034	REAL d_h_vcol_phy
1035	REAL fs_bound, fq_bound
1036	SAVE d_h_vcol_phy
1037	c$OMP THREADPRIVATE(d_h_vcol_phy)
1038	REAL zero_v(klon)
1039	CHARACTER*15 ztit
1040	INTEGER ip_ebil ! PRINT level for energy conserv. diag.
1041	SAVE ip_ebil
1042	DATA ip_ebil/0/
1043	c$OMP THREADPRIVATE(ip_ebil)
1044	INTEGER if_ebil ! level for energy conserv. dignostics
1045	SAVE if_ebil
1046	c$OMP THREADPRIVATE(if_ebil)
1047	c+jld ec_conser
1048	REAL ZRCPD
1049	c-jld ec_conser
1050	REAL t2m(klon,nbsrf) ! temperature a 2m
1051	REAL q2m(klon,nbsrf) ! humidite a 2m
1052
1053	cIM: t2m, q2m, u10m, v10m et t2mincels, t2maxcels
1054	REAL zt2m(klon), zq2m(klon) !temp., hum. 2m moyenne s/ 1 maille
1055	REAL zu10m(klon), zv10m(klon) !vents a 10m moyennes s/1 maille
1056	CHARACTER*40 t2mincels, t2maxcels !t2m min., t2m max
1057	CHARACTER*40 tinst, tave, typeval
1058	cjq Aerosol effects (Johannes Quaas, 27/11/2003)
1059	REAL sulfate(klon, klev) ! SO4 aerosol concentration [ug/m3]
1060
1061	REAL cldtaupi(klon,klev) ! Cloud optical thickness for pre-industrial (pi) aerosols
1062
1063	REAL re(klon, klev) ! Cloud droplet effective radius
1064	REAL fl(klon, klev) ! denominator of re
1065
1066	REAL re_top(klon), fl_top(klon) ! CDR at top of liquid water clouds
1067
1068	! Aerosol optical properties
1069
1070	! Aerosol optical properties by INCA model
1071	CHARACTER*4 :: rfname(9)
1072	REAL aerindex(klon) ! POLDER aerosol index
1073
1074	! Parameters
1075	LOGICAL ok_ade, ok_aie ! Apply aerosol (in)direct effects or not
1076	REAL bl95_b0, bl95_b1 ! Parameter in Boucher and Lohmann (1995)
1077	SAVE ok_ade, ok_aie, bl95_b0, bl95_b1
1078	c$OMP THREADPRIVATE(ok_ade, ok_aie, bl95_b0, bl95_b1)
1079	LOGICAL, SAVE :: aerosol_couple ! true : calcul des aerosols dans INCA
1080	! false : lecture des aerosol dans un fichier
1081	c$OMP THREADPRIVATE(aerosol_couple)
1082
1083	c
1084	c Declaration des constantes et des fonctions thermodynamiques
1085	c
1086	LOGICAL,SAVE :: first=.true.
1087	c$OMP THREADPRIVATE(first)
1088	#include "YOMCST.h"
1089	#include "YOETHF.h"
1090	#include "FCTTRE.h"
1091	cIM 100106 BEG : pouvoir sortir les ctes de la physique
1092	#include "conema3.h"
1093	#include "fisrtilp.h"
1094	#include "nuage.h"
1095	#include "compbl.h"
1096	cIM 100106 END : pouvoir sortir les ctes de la physique
1097	c
1098	c======================================================================
1099	! Ecriture eventuelle d'un profil verticale en entree de la physique.
1100	! Utilise notamment en 1D mais peut etre active egalement en 3D
1101	! en imposant la valeur de igout.
1102	c======================================================================
1103
1104	if (prt_level.ge.1) then
1105	igout=klon/2+1/klon
1106	write(lunout,*) 'DEBUT DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
1107	write(lunout,*)
1108	s 'nlon,nlev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys'
1109	write(lunout,*)
1110	s nlon,nlev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys
1111
1112	write(lunout,*) 'papers, play, phi, u, v, t, omega'
1113	do k=1,nlev
1114	write(lunout,*) paprs(igout,k),pplay(igout,k),pphi(igout,k),
1115	s u(igout,k),v(igout,k),t(igout,k),omega(igout,k)
1116	enddo
1117	write(lunout,*) 'ovap (g/kg), oliq (g/kg)'
1118	do k=1,nlev
1119	write(lunout,) qx(igout,k,1)1000,qx(igout,k,2)*1000.
1120	enddo
1121	endif
1122
1123	c======================================================================
1124
1125	cym => necessaire pour iflag_con != 2
1126	pmfd(:,:) = 0.
1127	pen_u(:,:) = 0.
1128	pen_d(:,:) = 0.
1129	pde_d(:,:) = 0.
1130	pde_u(:,:) = 0.
1131	aam=0.
1132
1133	torsfc=0.
1134
1135	if (first) then
1136
1137	cCR:nvelles variables convection/poches froides
1138
1139	print*, '================================================='
1140	print*, 'Allocation des variables locales et sauvegardees'
1141	call phys_local_var_init
1142	call phys_state_var_init
1143	print*, '================================================='
1144
1145	cIM beg
1146	dnwd0=0.0
1147	ftd=0.0
1148	fqd=0.0
1149	cin=0.
1150	cym Attention pbase pas initialise dans concvl !!!!
1151	pbase=0
1152	paire_ter(:)=0.
1153	cIM 180608
1154	c pmflxr=0.
1155	c pmflxs=0.
1156	first=.false.
1157
1158	endif ! fisrt
1159
1160	modname = 'physiq'
1161	cIM
1162	IF (ip_ebil_phy.ge.1) THEN
1163	DO i=1,klon
1164	zero_v(i)=0.
1165	END DO
1166	END IF
1167	ok_sync=.TRUE.
1168
1169	IF (debut) THEN
1170	CALL suphel ! initialiser constantes et parametres phys.
1171	ENDIF
1172
1173	if(prt_level.ge.1) print*,'CONVERGENCE PHYSIQUE THERM 1 '
1174
1175
1176	c======================================================================
1177	xjour = rjourvrai
1178	c
1179	c Si c'est le debut, il faut initialiser plusieurs choses
1180	c ********
1181	c
1182	IF (debut) THEN
1183	C
1184	!rv
1185	cCRinitialisation de wght_th et lalim_conv pour la definition de la couche alimentation
1186	cde la convection a partir des caracteristiques du thermique
1187	wght_th(:,:)=1.
1188	lalim_conv(:)=1
1189	cRC
1190	u10m(:,:)=0.
1191	v10m(:,:)=0.
1192	piz_ae(:,:,:)=0.
1193	tau_ae(:,:,:)=0.
1194	cg_ae(:,:,:)=0.
1195	rain_con(:)=0.
1196	snow_con(:)=0.
1197	bl95_b0=0.
1198	bl95_b1=0.
1199	topswai(:)=0.
1200	topswad(:)=0.
1201	solswai(:)=0.
1202	solswad(:)=0.
1203
1204	lambda_th(:,:)=0.
1205	wmax_th(:)=0.
1206	tau_overturning_th(:)=0.
1207	IF (config_inca /= 'none') THEN
1208	tau_inca(:,:,:,:) = 0.
1209	piz_inca(:,:,:,:) = 0.
1210	cg_inca(:,:,:,:) = 0.
1211	ccm(:,:,:) = 0.
1212	topswai_inca(:) = 0.
1213	topswad_inca(:) = 0.
1214	topswad0_inca(:) = 0.
1215	topsw_inca(:,:) = 0.
1216	topsw0_inca(:,:) = 0.
1217	solswai_inca(:) = 0.
1218	solswad_inca(:) = 0.
1219	solswad0_inca(:) = 0.
1220	solsw_inca(:,:) = 0.
1221	solsw0_inca(:,:) = 0.
1222	END IF
1223
1224	rnebcon0(:,:) = 0.0
1225	clwcon0(:,:) = 0.0
1226	rnebcon(:,:) = 0.0
1227	clwcon(:,:) = 0.0
1228
1229	cIM
1230	IF (ip_ebil_phy.ge.1) d_h_vcol_phy=0.
1231	c
1232	c appel a la lecture du run.def physique
1233	c
1234	call conf_phys(ok_journe, ok_mensuel,
1235	. ok_instan, ok_hf,
1236	. ok_LES,
1237	. solarlong0,seuil_inversion,
1238	. fact_cldcon, facttemps,ok_newmicro,iflag_radia,
1239	. iflag_cldcon,iflag_ratqs,ratqsbas,ratqshaut,
1240	. ok_ade, ok_aie, aerosol_couple,
1241	. bl95_b0, bl95_b1,
1242	. iflag_thermals,nsplit_thermals,tau_thermals,
1243	. iflag_thermals_ed,iflag_thermals_optflux,
1244	cnv flags pour la convection et les poches froides
1245	. iflag_coupl,iflag_clos,iflag_wake)
1246
1247	print*,'iflag_coupl,iflag_clos,iflag_wake',
1248	. iflag_coupl,iflag_clos,iflag_wake
1249	print*,'CYCLE_DIURNE', cycle_diurne
1250	c
1251	IF (iflag_con.EQ.2.AND.iflag_cldcon.GT.-1) THEN
1252	abort_message = 'Tiedtke needs iflag_cldcon=-2 or -1'
1253	CALL abort_gcm (modname,abort_message,1)
1254	ENDIF
1255	c
1256	IF(ok_isccp.AND.iflag_con.LE.2) THEN
1257	abort_message = 'ISCCP-like outputs may be available for KE
1258	.(iflag_con >= 3); for Tiedtke (iflag_con=-2) put ok_isccp=n'
1259	CALL abort_gcm (modname,abort_message,1)
1260	ENDIF
1261	c
1262	c Initialiser les compteurs:
1263	c
1264	itap = 0
1265	itaprad = 0
1266
1267	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1268	!! Un petit travail à faire ici.
1269	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1270
1271	if (iflag_pbl>1) then
1272	PRINT*, "Using method MELLOR&YAMADA"
1273	endif
1274
1275	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1276	! FH 2008/05/02 changement lie a la lecture de nbapp_rad dans phylmd plutot que
1277	! dyn3d
1278	! Attention : la version precedente n'etait pas tres propre.
1279	! Il se peut qu'il faille prendre une valeur differente de nbapp_rad
1280	! pour obtenir le meme resultat.
1281	dtime=pdtphys
1282	radpas = NINT( 86400./dtime/nbapp_rad)
1283	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1284
1285	CALL phyetat0 ("startphy.nc",clesphy0,tabcntr0)
1286	cIM begin
1287	print*,'physiq: clwcon rnebcon ratqs',clwcon(1,1),rnebcon(1,1)
1288	$,ratqs(1,1)
1289	cIM end
1290
1291
1292
1293	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1294	c
1295	C on remet le calendrier a zero
1296	c
1297	IF (raz_date .eq. 1) THEN
1298	itau_phy = 0
1299	ENDIF
1300
1301	cIM cf. AM 081204 BEG
1302	PRINT*,'cycle_diurne3 =',cycle_diurne
1303	cIM cf. AM 081204 END
1304	c
1305	CALL printflag( tabcntr0,radpas,ok_journe,
1306	, ok_instan, ok_region )
1307	c
1308	IF (ABS(dtime-pdtphys).GT.0.001) THEN
1309	WRITE(lunout,*) 'Pas physique n est pas correct',dtime,
1310	. pdtphys
1311	abort_message='Pas physique n est pas correct '
1312	! call abort_gcm(modname,abort_message,1)
1313	dtime=pdtphys
1314	ENDIF
1315	IF (nlon .NE. klon) THEN
1316	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,
1317	. klon
1318	abort_message='nlon et klon ne sont pas coherents'
1319	call abort_gcm(modname,abort_message,1)
1320	ENDIF
1321	IF (nlev .NE. klev) THEN
1322	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,
1323	. klev
1324	abort_message='nlev et klev ne sont pas coherents'
1325	call abort_gcm(modname,abort_message,1)
1326	ENDIF
1327	c
1328	IF (dtime*FLOAT(radpas).GT.21600..AND.cycle_diurne) THEN
1329	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
1330	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
1331	abort_message='Nbre d appels au rayonnement insuffisant'
1332	call abort_gcm(modname,abort_message,1)
1333	ENDIF
1334	WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
1335	WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=",
1336	. ok_cvl
1337	c
1338	cKE43
1339	c Initialisation pour la convection de K.E. (sb):
1340	IF (iflag_con.GE.3) THEN
1341
1342	WRITE(lunout,)"** Convection de Kerry Emanuel 4.3 "
1343	WRITE(lunout,*)
1344	. "On va utiliser le melange convectif des traceurs qui"
1345	WRITE(lunout,*)"est calcule dans convect4.3"
1346	WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
1347
1348	DO i = 1, klon
1349	ema_cbmf(i) = 0.
1350	ema_pcb(i) = 0.
1351	ema_pct(i) = 0.
1352	ema_workcbmf(i) = 0.
1353	ENDDO
1354	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
1355	DO i = 1, klon
1356	ibas_con(i) = 1
1357	itop_con(i) = 1
1358	ENDDO
1359	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
1360	c===============================================================================
1361	cCR:04.12.07: initialisations poches froides
1362	c Controle de ALE et ALP pour la fermeture convective (jyg)
1363	if (iflag_wake.eq.1) then
1364	CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr
1365	s ,alp_bl_prescr, ale_bl_prescr)
1366	c 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
1367	c print*,'apres ini_wake iflag_cldcon=', iflag_cldcon
1368	endif
1369
1370	do i = 1,klon
1371	Ale_bl(i)=0.
1372	Alp_bl(i)=0.
1373	enddo
1374
1375	c================================================================================
1376
1377	ENDIF
1378
1379	rugoro=0.
1380	c34EK
1381	IF (ok_orodr) THEN
1382
1383	rugoro=0.
1384
1385	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1386	! FH sans doute a enlever de finitivement ou, si on le garde, l'activer
1387	! justement quand ok_orodr = false.
1388	! ce rugoro est utilise par la couche limite et fait double emploi
1389	! avec les paramétrisations spécifiques de Francois Lott.
1390	! DO i=1,klon
1391	! rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1392	! ENDDO
1393	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1394	IF (ok_strato) THEN
1395	CALL SUGWD_strato(klon,klev,paprs,pplay)
1396	ELSE
1397	CALL SUGWD(klon,klev,paprs,pplay)
1398	ENDIF
1399
1400	DO i=1,klon
1401	zuthe(i)=0.
1402	zvthe(i)=0.
1403	if(zstd(i).gt.10.)then
1404	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
1405	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
1406	endif
1407	ENDDO
1408	ENDIF
1409	c
1410	c
1411	lmt_pas = NINT(86400./dtime * 1.0) ! tous les jours
1412	WRITE(lunout,*)'La frequence de lecture surface est de ',
1413	. lmt_pas
1414	c
1415	cIM200505 ecrit_mth = NINT(86400./dtime *ecritphy) ! tous les ecritphy jours
1416	c IF (ok_mensuel) THEN
1417	c WRITE(lunout,*)'La frequence de sortie mensuelle est de ',
1418	c . ecrit_mth
1419	c ENDIF
1420	c ecrit_day = NINT(86400./dtime *1.0) ! tous les jours
1421	c IF (ok_journe) THEN
1422	c WRITE(lunout,*)'La frequence de sortie journaliere est de ',
1423	c . ecrit_day
1424	c ENDIF
1425	cIM 130904 BEG
1426	cIM 080205 ecrit_hf = 86400./dtime *0.25 ! toutes les 6h
1427	cIM 170305
1428	c ecrit_hf = 86400./dtime/12. ! toutes les 2h
1429	cIM 230305
1430	cIM200505 ecrit_hf = 86400./dtime *0.25 ! toutes les 6h
1431	c
1432	cIM200505 ecrit_hf2mth = ecrit_day/ecrit_hf*30
1433	c
1434	cIM200505 IF (ok_journe) THEN
1435	cIM200505 WRITE(lunout,*)'La frequence de sortie hf est de ',
1436	cIM200505 . ecrit_hf
1437	cIM200505 ENDIF
1438	cIM 130904 END
1439	ccc ecrit_ins = NINT(86400./dtime *0.5) ! 2 fois par jour
1440	ccc ecrit_ins = NINT(86400./dtime *0.25) ! 4 fois par jour
1441	c ecrit_ins = NINT(86400./dtime/48.) ! a chaque pas de temps ==> PB. dans time_counter pour 1mois
1442	c ecrit_ins = NINT(86400./dtime/12.) ! toutes les deux heures
1443	cIM200505 ecrit_ins = NINT(86400./dtime/8.) ! toutes les trois heures
1444	cIM200505 IF (ok_instan) THEN
1445	cIM200505 WRITE(lunout,*)'La frequence de sortie instant. est de ',
1446	cIM200505 . ecrit_ins
1447	cIM200505 ENDIF
1448	cIM200505 ecrit_reg = NINT(86400./dtime *0.25) ! 4 fois par jour
1449	cIM200505 IF (ok_region) THEN
1450	cIM200505 WRITE(lunout,*)'La frequence de sortie region est de ',
1451	cIM200505 . ecrit_reg
1452	cIM200505 ENDIF
1453	cIM 030306 BEG
1454	cIM ecrit_hf2mth = nombre de pas de temps de calcul de hf par mois apres lequel on ecrit
1455	cIM : ne pas modifier ecrit_hf2mth
1456	c
1457	cIM 250308bad guide ecrit_hf2mth = 30*1/ecrit_hf
1458	ecrit_hf2mth = ecrit_mth/ecrit_hf
1459	c ecrit_ins en secondes, chaque pas de temps de la physique
1460	ecrit_ins = dtime
1461	cIM on passe les frequences de jours en secondes : ecrit_ins, ecrit_hf, ecrit_day, ecrit_mth, ecrit_tra, ecrit_reg
1462	ecrit_hf = ecrit_hf * un_jour
1463	!IM
1464	IF(ecrit_day.LE.1.) THEN
1465	ecrit_day = ecrit_day * un_jour !en secondes
1466	ENDIF
1467	!IM
1468	ecrit_mth = ecrit_mth * un_jour
1469	ecrit_reg = ecrit_reg * un_jour
1470	ecrit_tra = ecrit_tra * un_jour
1471	ecrit_ISCCP = ecrit_ISCCP * un_jour
1472	ecrit_LES = ecrit_LES * un_jour
1473	c
1474	PRINT*,'physiq ecrit_ hf day mth reg tra ISCCP hf2mth',
1475	. ecrit_hf,ecrit_day,ecrit_mth,ecrit_reg,ecrit_tra,ecrit_ISCCP,
1476	. ecrit_hf2mth
1477	cIM 030306 END
1478
1479	capemaxcels = 't_max(X)'
1480	t2mincels = 't_min(X)'
1481	t2maxcels = 't_max(X)'
1482	tinst = 'inst(X)'
1483	tave = 'ave(X)'
1484	cIM cf. AM 081204 BEG
1485	write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
1486	cIM cf. AM 081204 END
1487	c
1488	c=============================================================
1489	c Initialisation des sorties
1490	c=============================================================
1491
1492	#ifdef CPP_IOIPSL
1493
1494	c$OMP MASTER
1495	call phys_output_open(jjmp1,nlevSTD,clevSTD,nbteta,
1496	& ctetaSTD,dtime,presnivs,ok_veget,
1497	& type_ocean,iflag_pbl,ok_mensuel,ok_journe,
1498	& ok_hf,ok_instan,ok_LES,ok_ade,ok_aie)
1499	c$OMP END MASTER
1500	c$OMP BARRIER
1501
1502	#ifdef histISCCP
1503	#include "ini_histISCCP.h"
1504	#endif
1505
1506	#ifdef histmthNMC
1507	#include "ini_histmthNMC.h"
1508	#endif
1509
1510	#include "ini_histday_seri.h"
1511
1512	#include "ini_paramLMDZ_phy.h"
1513
1514	#endif
1515
1516	cXXXPB Positionner date0 pour initialisation de ORCHIDEE
1517	date0 = zjulian
1518	C date0 = day_ini
1519	WRITE(,) 'physiq date0 : ',date0
1520	c
1521	c
1522	c
1523	c Prescrire l'ozone dans l'atmosphere
1524	c
1525	c
1526	cc DO i = 1, klon
1527	cc DO k = 1, klev
1528	cc CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1529	cc ENDDO
1530	cc ENDDO
1531	c
1532	IF (config_inca /= 'none') THEN
1533	#ifdef INCA
1534	CALL VTe(VTphysiq)
1535	CALL VTb(VTinca)
1536	iii = MOD(NINT(xjour),360)
1537	calday = FLOAT(iii) + gmtime
1538	WRITE(lunout,*) 'initial time ', xjour, calday
1539
1540	CALL chemini(
1541	$ rg,
1542	$ ra,
1543	$ airephy,
1544	$ rlat,
1545	$ rlon,
1546	$ presnivs,
1547	$ calday,
1548	$ klon,
1549	$ nqtot,
1550	$ pdtphys,
1551	$ annee_ref,
1552	$ day_ini)
1553
1554	CALL VTe(VTinca)
1555	CALL VTb(VTphysiq)
1556	#endif
1557	END IF
1558	c
1559	c
1560	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1561	! Nouvelle initialisation pour le rayonnement RRTM
1562	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1563
1564	call iniradia(klon,klev,paprs(1,1:klev+1))
1565
1566	ENDIF
1567	!
1568	! ************** Fin de IF ( debut ) *************
1569	!
1570	!
1571	! Incrementer le compteur de la physique
1572	!
1573	itap = itap + 1
1574	julien = MOD(NINT(xjour),360)
1575	if (julien .eq. 0) julien = 360
1576
1577	!
1578	! Update fraction of the sub-surfaces (pctsrf) and
1579	! initialize, where a new fraction has appeared, all variables depending
1580	! on the surface fraction.
1581	!
1582	CALL change_srf_frac(itap, dtime, julien,
1583	* pctsrf, falb1, falb2, ftsol, u10m, v10m, pbl_tke)
1584
1585
1586	! Tendances bidons pour les processus qui n'affectent pas certaines
1587	! variables.
1588	du0(:,:)=0.
1589	dv0(:,:)=0.
1590	dq0(:,:)=0.
1591	dql0(:,:)=0.
1592	c
1593	c Mettre a zero des variables de sortie (pour securite)
1594	c
1595	DO i = 1, klon
1596	d_ps(i) = 0.0
1597	ENDDO
1598	DO k = 1, klev
1599	DO i = 1, klon
1600	d_t(i,k) = 0.0
1601	d_u(i,k) = 0.0
1602	d_v(i,k) = 0.0
1603	ENDDO
1604	ENDDO
1605	DO iq = 1, nqtot
1606	DO k = 1, klev
1607	DO i = 1, klon
1608	d_qx(i,k,iq) = 0.0
1609	ENDDO
1610	ENDDO
1611	ENDDO
1612	da(:,:)=0.
1613	mp(:,:)=0.
1614	phi(:,:,:)=0.
1615	c
1616	c Ne pas affecter les valeurs entrees de u, v, h, et q
1617	c
1618	DO k = 1, klev
1619	DO i = 1, klon
1620	t_seri(i,k) = t(i,k)
1621	u_seri(i,k) = u(i,k)
1622	v_seri(i,k) = v(i,k)
1623	q_seri(i,k) = qx(i,k,ivap)
1624	ql_seri(i,k) = qx(i,k,iliq)
1625	qs_seri(i,k) = 0.
1626	ENDDO
1627	ENDDO
1628	IF (nqtot.GE.3) THEN
1629	DO iq = 3, nqtot
1630	DO k = 1, klev
1631	DO i = 1, klon
1632	tr_seri(i,k,iq-2) = qx(i,k,iq)
1633	ENDDO
1634	ENDDO
1635	ENDDO
1636	ELSE
1637	DO k = 1, klev
1638	DO i = 1, klon
1639	tr_seri(i,k,1) = 0.0
1640	ENDDO
1641	ENDDO
1642	ENDIF
1643	C
1644	DO i = 1, klon
1645	ztsol(i) = 0.
1646	ENDDO
1647	DO nsrf = 1, nbsrf
1648	DO i = 1, klon
1649	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1650	ENDDO
1651	ENDDO
1652	cIM
1653	IF (ip_ebil_phy.ge.1) THEN
1654	ztit='after dynamic'
1655	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
1656	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1657	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1658	C Comme les tendances de la physique sont ajoute dans la dynamique,
1659	C on devrait avoir que la variation d'entalpie par la dynamique
1660	C est egale a la variation de la physique au pas de temps precedent.
1661	C Donc la somme de ces 2 variations devrait etre nulle.
1662	call diagphy(airephy,ztit,ip_ebil_phy
1663	e , zero_v, zero_v, zero_v, zero_v, zero_v
1664	e , zero_v, zero_v, zero_v, ztsol
1665	e , d_h_vcol+d_h_vcol_phy, d_qt, 0.
1666	s , fs_bound, fq_bound )
1667	END IF
1668
1669	c Diagnostiquer la tendance dynamique
1670	c
1671	IF (ancien_ok) THEN
1672	DO k = 1, klev
1673	DO i = 1, klon
1674	d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
1675	d_v_dyn(i,k) = (v_seri(i,k)-v_ancien(i,k))/dtime
1676	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
1677	d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
1678	ENDDO
1679	ENDDO
1680	ELSE
1681	DO k = 1, klev
1682	DO i = 1, klon
1683	d_u_dyn(i,k) = 0.0
1684	d_v_dyn(i,k) = 0.0
1685	d_t_dyn(i,k) = 0.0
1686	d_q_dyn(i,k) = 0.0
1687	ENDDO
1688	ENDDO
1689	ancien_ok = .TRUE.
1690	ENDIF
1691	c
1692	c Ajouter le geopotentiel du sol:
1693	c
1694	DO k = 1, klev
1695	DO i = 1, klon
1696	zphi(i,k) = pphi(i,k) + pphis(i)
1697	ENDDO
1698	ENDDO
1699	c
1700	c Verifier les temperatures
1701	c
1702	cIM BEG
1703	IF (check) THEN
1704	amn=MIN(ftsol(1,is_ter),1000.)
1705	amx=MAX(ftsol(1,is_ter),-1000.)
1706	DO i=2, klon
1707	amn=MIN(ftsol(i,is_ter),amn)
1708	amx=MAX(ftsol(i,is_ter),amx)
1709	ENDDO
1710	c
1711	PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
1712	ENDIF !(check) THEN
1713	cIM END
1714	c
1715	CALL hgardfou(t_seri,ftsol,'debutphy')
1716	c
1717	cIM BEG
1718	IF (check) THEN
1719	amn=MIN(ftsol(1,is_ter),1000.)
1720	amx=MAX(ftsol(1,is_ter),-1000.)
1721	DO i=2, klon
1722	amn=MIN(ftsol(i,is_ter),amn)
1723	amx=MAX(ftsol(i,is_ter),amx)
1724	ENDDO
1725	c
1726	PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
1727	ENDIF !(check) THEN
1728	cIM END
1729	c
1730	c Mettre en action les conditions aux limites (albedo, sst, etc.).
1731	c Prescrire l'ozone et calculer l'albedo sur l'ocean.
1732	c
1733	IF (MOD(itap-1,lmt_pas) .EQ. 0) THEN
1734	if(prt_level.ge.1) WRITE(lunout,*)' PHYS cond julien ',julien
1735	CALL ozonecm( FLOAT(julien), rlat, paprs, wo)
1736	ENDIF
1737	c
1738	c Re-evaporer l'eau liquide nuageuse
1739	c
1740	DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse
1741	DO i = 1, klon
1742	zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1743	c zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1744	zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1745	zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
1746	zb = MAX(0.0,ql_seri(i,k))
1747	za = - MAX(0.0,ql_seri(i,k))
1748	. * (zlvdcp(1.-zdelta)+zlsdcpzdelta)
1749	t_seri(i,k) = t_seri(i,k) + za
1750	q_seri(i,k) = q_seri(i,k) + zb
1751	ql_seri(i,k) = 0.0
1752	d_t_eva(i,k) = za
1753	d_q_eva(i,k) = zb
1754	ENDDO
1755	ENDDO
1756	cIM
1757	IF (ip_ebil_phy.ge.2) THEN
1758	ztit='after reevap'
1759	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,1,dtime
1760	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1761	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1762	call diagphy(airephy,ztit,ip_ebil_phy
1763	e , zero_v, zero_v, zero_v, zero_v, zero_v
1764	e , zero_v, zero_v, zero_v, ztsol
1765	e , d_h_vcol, d_qt, d_ec
1766	s , fs_bound, fq_bound )
1767	C
1768	END IF
1769
1770	c
1771	c=========================================================================
1772	! Calculs de l'orbite.
1773	! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
1774	! doit donc etre placé avant radlwsw et pbl_surface
1775
1776	! choix entre calcul de la longitude solaire vraie ou valeur fixee a
1777	! solarlong0
1778
1779	if (solarlong0<-999.) then
1780	CALL orbite(FLOAT(julien),zlongi,dist)
1781	else
1782	zlongi=solarlong0 ! longitude solaire vraie
1783	dist=1. ! distance au soleil / moyenne
1784	endif
1785
1786	if(prt_level.ge.1) print*,'Longitude solaire ',zlongi,solarlong0
1787
1788	! Avec ou sans cycle diurne
1789	IF (cycle_diurne) THEN
1790	zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
1791	CALL zenang(zlongi,gmtime,zdtime,rlat,rlon,rmu0,fract)
1792	ELSE
1793	CALL angle(zlongi, rlat, fract, rmu0)
1794	ENDIF
1795
1796	if (mydebug) then
1797	call writefield_phy('u_seri',u_seri,llm)
1798	call writefield_phy('v_seri',v_seri,llm)
1799	call writefield_phy('t_seri',t_seri,llm)
1800	call writefield_phy('q_seri',q_seri,llm)
1801	endif
1802
1803	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1804	c Appel au pbl_surface : Planetary Boudary Layer et Surface
1805	c Cela implique tous les interactions des sous-surfaces et la partie diffusion
1806	c turbulent du couche limit.
1807	c
1808	c Certains varibales de sorties de pbl_surface sont utiliser que pour
1809	c ecriture des fihiers hist_XXXX.nc, ces sont :
1810	c qsol, zq2m, s_pblh, s_lcl,
1811	c s_capCL, s_oliqCL, s_cteiCL,s_pblT,
1812	c s_therm, s_trmb1, s_trmb2, s_trmb3,
1813	c zxrugs, zu10m, zv10m, fder,
1814	c zxqsurf, rh2m, zxfluxu, zxfluxv,
1815	c frugs, agesno, fsollw, fsolsw,
1816	c d_ts, fevap, fluxlat, t2m,
1817	c wfbils, wfbilo, fluxt, fluxu, fluxv,
1818	c
1819	c Certains ne sont pas utiliser du tout :
1820	c dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
1821	c
1822
1823	CALL pbl_surface(
1824	e dtime, date0, itap, julien,
1825	e debut, lafin,
1826	e rlon, rlat, rugoro, rmu0,
1827	e rain_fall, snow_fall, solsw, sollw,
1828	e t_seri, q_seri, u_seri, v_seri,
1829	e pplay, paprs, pctsrf,
1830	+ ftsol, falb1, falb2, u10m, v10m,
1831	s sollwdown, cdragh, cdragm, u1, v1,
1832	s albsol1, albsol2, sens, evap,
1833	s zxtsol, zxfluxlat, zt2m, qsat2m,
1834	s d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf,
1835	s coefh, slab_wfbils,
1836	d qsol, zq2m, s_pblh, s_lcl,
1837	d s_capCL, s_oliqCL, s_cteiCL,s_pblT,
1838	d s_therm, s_trmb1, s_trmb2, s_trmb3,
1839	d zxrugs, zu10m, zv10m, fder,
1840	d zxqsurf, rh2m, zxfluxu, zxfluxv,
1841	d frugs, agesno, fsollw, fsolsw,
1842	d d_ts, fevap, fluxlat, t2m,
1843	d wfbils, wfbilo, fluxt, fluxu, fluxv,
1844	- dsens, devap, zxsnow,
1845	- zxfluxt, zxfluxq, q2m, fluxq, pbl_tke )
1846
1847
1848	!-----------------------------------------------------------------------------------------
1849	! ajout des tendances de la diffusion turbulente
1850	CALL add_phys_tend(d_u_vdf,d_v_vdf,d_t_vdf,d_q_vdf,dql0,'vdf')
1851	!-----------------------------------------------------------------------------------------
1852
1853	if (mydebug) then
1854	call writefield_phy('u_seri',u_seri,llm)
1855	call writefield_phy('v_seri',v_seri,llm)
1856	call writefield_phy('t_seri',t_seri,llm)
1857	call writefield_phy('q_seri',q_seri,llm)
1858	endif
1859
1860
1861	IF (ip_ebil_phy.ge.2) THEN
1862	ztit='after surface_main'
1863	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
1864	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1865	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1866	call diagphy(airephy,ztit,ip_ebil_phy
1867	e , zero_v, zero_v, zero_v, zero_v, sens
1868	e , evap , zero_v, zero_v, ztsol
1869	e , d_h_vcol, d_qt, d_ec
1870	s , fs_bound, fq_bound )
1871	END IF
1872
1873	c =================================================================== c
1874	c Calcul de Qsat
1875
1876	DO k = 1, klev
1877	DO i = 1, klon
1878	zx_t = t_seri(i,k)
1879	IF (thermcep) THEN
1880	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
1881	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
1882	zx_qs = MIN(0.5,zx_qs)
1883	zcor = 1./(1.-retv*zx_qs)
1884	zx_qs = zx_qs*zcor
1885	ELSE
1886	IF (zx_t.LT.t_coup) THEN
1887	zx_qs = qsats(zx_t)/pplay(i,k)
1888	ELSE
1889	zx_qs = qsatl(zx_t)/pplay(i,k)
1890	ENDIF
1891	ENDIF
1892	zqsat(i,k)=zx_qs
1893	ENDDO
1894	ENDDO
1895
1896	if (prt_level.ge.1) then
1897	write(lunout,*) 'L qsat (g/kg) avant clouds_gno'
1898	write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
1899	endif
1900	c
1901	c Appeler la convection (au choix)
1902	c
1903	DO k = 1, klev
1904	DO i = 1, klon
1905	conv_q(i,k) = d_q_dyn(i,k)
1906	. + d_q_vdf(i,k)/dtime
1907	conv_t(i,k) = d_t_dyn(i,k)
1908	. + d_t_vdf(i,k)/dtime
1909	ENDDO
1910	ENDDO
1911	IF (check) THEN
1912	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
1913	WRITE(lunout,*) "avantcon=", za
1914	ENDIF
1915	zx_ajustq = .FALSE.
1916	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
1917	IF (zx_ajustq) THEN
1918	DO i = 1, klon
1919	z_avant(i) = 0.0
1920	ENDDO
1921	DO k = 1, klev
1922	DO i = 1, klon
1923	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k))
1924	. *(paprs(i,k)-paprs(i,k+1))/RG
1925	ENDDO
1926	ENDDO
1927	ENDIF
1928
1929	c Calcule de vitesse verticale a partir de flux de masse verticale
1930	DO k = 1, klev
1931	DO i = 1, klon
1932	omega(i,k) = RG*flxmass_w(i,k) / airephy(i)
1933	END DO
1934	END DO
1935
1936	IF (iflag_con.EQ.1) THEN
1937	stop'reactiver le call conlmd dans physiq.F'
1938	c CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
1939	c . d_t_con, d_q_con,
1940	c . rain_con, snow_con, ibas_con, itop_con)
1941	ELSE IF (iflag_con.EQ.2) THEN
1942	CALL conflx(dtime, paprs, pplay, t_seri, q_seri,
1943	e conv_t, conv_q, -evap, omega,
1944	s d_t_con, d_q_con, rain_con, snow_con,
1945	s pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
1946	s kcbot, kctop, kdtop, pmflxr, pmflxs)
1947	d_u_con = 0.
1948	d_v_con = 0.
1949
1950	WHERE (rain_con < 0.) rain_con = 0.
1951	WHERE (snow_con < 0.) snow_con = 0.
1952	DO i = 1, klon
1953	ibas_con(i) = klev+1 - kcbot(i)
1954	itop_con(i) = klev+1 - kctop(i)
1955	ENDDO
1956	ELSE IF (iflag_con.GE.3) THEN
1957	c nb of tracers for the KE convection:
1958	c MAF la partie traceurs est faite dans phytrac
1959	c on met ntra=1 pour limiter les appels mais on peut
1960	c supprimer les calculs / ftra.
1961	ntra = 1
1962
1963	c=====================================================================================
1964	cajout pour la parametrisation des poches froides:
1965	ccalcul de t_wake et t_undi: si pas de poches froides, t_wake=t_undi=t_seri
1966	do k=1,klev
1967	do i=1,klon
1968	if (iflag_wake.eq.1) then
1969	t_wake(i,k) = t_seri(i,k)
1970	. +(1-wake_s(i))*wake_deltat(i,k)
1971	q_wake(i,k) = q_seri(i,k)
1972	. +(1-wake_s(i))*wake_deltaq(i,k)
1973	t_undi(i,k) = t_seri(i,k)
1974	. -wake_s(i)*wake_deltat(i,k)
1975	q_undi(i,k) = q_seri(i,k)
1976	. -wake_s(i)*wake_deltaq(i,k)
1977	else
1978	t_wake(i,k) = t_seri(i,k)
1979	q_wake(i,k) = q_seri(i,k)
1980	t_undi(i,k) = t_seri(i,k)
1981	q_undi(i,k) = q_seri(i,k)
1982	endif
1983	enddo
1984	enddo
1985
1986	cc-- Calcul de l'energie disponible ALE (J/kg) et de la puissance disponible ALP (W/m2)
1987	cc-- pour le soulevement des particules dans le modele convectif
1988	c
1989	do i = 1,klon
1990	ALE(i) = 0.
1991	ALP(i) = 0.
1992	enddo
1993	c
1994	ccalcul de ale_wake et alp_wake
1995	do i = 1,klon
1996	if (iflag_wake.eq.1) then
1997	ale_wake(i) = 0.5wake_cstar(i)*2
1998	alp_wake(i) = wake_fip(i)
1999	else
2000	ale_wake(i) = 0.
2001	alp_wake(i) = 0.
2002	endif
2003	enddo
2004	ccombinaison avec ale et alp de couche limite: constantes si pas de couplage, valeurs calculees
2005	cdans le thermique sinon
2006	if (iflag_coupl.eq.0) then
2007	if (debut) print*,'ALE et ALP imposes'
2008	do i = 1,klon
2009	con ne couple que ale
2010	c ALE(i) = max(ale_wake(i),Ale_bl(i))
2011	ALE(i) = max(ale_wake(i),ale_bl_prescr)
2012	con ne couple que alp
2013	c ALP(i) = alp_wake(i) + Alp_bl(i)
2014	ALP(i) = alp_wake(i) + alp_bl_prescr
2015	enddo
2016	else
2017	IF(prt_level>9)WRITE(lunout,*)'ALE et ALP couples au thermique'
2018	do i = 1,klon
2019	ALE(i) = max(ale_wake(i),Ale_bl(i))
2020	ALP(i) = alp_wake(i) + Alp_bl(i)
2021	c write(20,*)'ALE',ALE(i),Ale_bl(i),ale_wake(i)
2022	c write(21,*)'ALP',ALP(i),Alp_bl(i),alp_wake(i)
2023	enddo
2024	endif
2025	do i=1,klon
2026	if (alp(i)>alp_max) then
2027	print*,'WARNING SUPER ALP (seuil=',alp_max,
2028	, '): i, alp, alp_wake,ale',i,alp(i),alp_wake(i),ale(i)
2029	alp(i)=alp_max
2030	endif
2031	if (ale(i)>ale_max) then
2032	print*,'WARNING SUPER ALE (seuil=',ale_max,
2033	, '): i, alp, alp_wake,ale',i,ale(i),ale_wake(i),alp(i)
2034	ale(i)=ale_max
2035	endif
2036	enddo
2037
2038	cfin calcul ale et alp
2039	c=================================================================================================
2040
2041
2042	c sb, oct02:
2043	c Schema de convection modularise et vectorise:
2044	c (driver commun aux versions 3 et 4)
2045	c
2046	IF (ok_cvl) THEN ! new driver for convectL
2047
2048	CALL concvl (iflag_con,iflag_clos,
2049	. dtime,paprs,pplay,t_undi,q_undi,
2050	. t_wake,q_wake,
2051	. u_seri,v_seri,tr_seri,nbtr,
2052	. ALE,ALP,
2053	. ema_work1,ema_work2,
2054	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2055	. rain_con, snow_con, ibas_con, itop_con, sigd,
2056	. upwd,dnwd,dnwd0,
2057	. Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl,
2058	. pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd,
2059	. pmflxr,pmflxs,da,phi,mp,
2060	. ftd,fqd,lalim_conv,wght_th)
2061
2062	cIM begin
2063	c print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
2064	c .dnwd0(1,1),ftd(1,1),fqd(1,1)
2065	cIM end
2066	cIM cf. FH
2067	clwcon0=qcondc
2068	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
2069
2070	ELSE ! ok_cvl
2071	c MAF conema3 ne contient pas les traceurs
2072	CALL conema3 (dtime,
2073	. paprs,pplay,t_seri,q_seri,
2074	. u_seri,v_seri,tr_seri,ntra,
2075	. ema_work1,ema_work2,
2076	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2077	. rain_con, snow_con, ibas_con, itop_con,
2078	. upwd,dnwd,dnwd0,bas,top,
2079	. Ma,cape,tvp,rflag,
2080	. pbase
2081	. ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr
2082	. ,clwcon0)
2083
2084	ENDIF ! ok_cvl
2085
2086	c
2087	c Correction precip
2088	rain_con = rain_con * cvl_corr
2089	snow_con = snow_con * cvl_corr
2090	c
2091
2092	IF (.NOT. ok_gust) THEN
2093	do i = 1, klon
2094	wd(i)=0.0
2095	enddo
2096	ENDIF
2097
2098	c =================================================================== c
2099	c Calcul des proprietes des nuages convectifs
2100	c
2101
2102	c calcul des proprietes des nuages convectifs
2103	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
2104	call clouds_gno
2105	s (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
2106
2107	c =================================================================== c
2108
2109	DO i = 1, klon
2110	ema_pcb(i) = pbase(i)
2111	ENDDO
2112	DO i = 1, klon
2113
2114	! L'idicage de itop_con peut cacher un pb potentiel
2115	! FH sous la dictee de JYG, CR
2116	ema_pct(i) = paprs(i,itop_con(i)+1)
2117
2118	if (itop_con(i).gt.klev-3) then
2119	print*,'La convection monte trop haut '
2120	print*,'itop_con(,',i,',)=',itop_con(i)
2121	endif
2122	ENDDO
2123	DO i = 1, klon
2124	ema_cbmf(i) = ema_workcbmf(i)
2125	ENDDO
2126	ELSE IF (iflag_con.eq.0) THEN
2127	write(lunout,*) 'On n appelle pas la convection'
2128	clwcon0=0.
2129	rnebcon0=0.
2130	d_t_con=0.
2131	d_q_con=0.
2132	d_u_con=0.
2133	d_v_con=0.
2134	rain_con=0.
2135	snow_con=0.
2136	bas=1
2137	top=1
2138	ELSE
2139	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
2140	CALL abort
2141	ENDIF
2142
2143	c CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
2144	c . d_u_con, d_v_con)
2145
2146	!-----------------------------------------------------------------------------------------
2147	! ajout des tendances de la diffusion turbulente
2148	CALL add_phys_tend(d_u_con,d_v_con,d_t_con,d_q_con,dql0,'con')
2149	!-----------------------------------------------------------------------------------------
2150
2151	if (mydebug) then
2152	call writefield_phy('u_seri',u_seri,llm)
2153	call writefield_phy('v_seri',v_seri,llm)
2154	call writefield_phy('t_seri',t_seri,llm)
2155	call writefield_phy('q_seri',q_seri,llm)
2156	endif
2157
2158	cIM
2159	IF (ip_ebil_phy.ge.2) THEN
2160	ztit='after convect'
2161	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2162	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2163	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2164	call diagphy(airephy,ztit,ip_ebil_phy
2165	e , zero_v, zero_v, zero_v, zero_v, zero_v
2166	e , zero_v, rain_con, snow_con, ztsol
2167	e , d_h_vcol, d_qt, d_ec
2168	s , fs_bound, fq_bound )
2169	END IF
2170	C
2171	IF (check) THEN
2172	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2173	WRITE(lunout,*)"aprescon=", za
2174	zx_t = 0.0
2175	za = 0.0
2176	DO i = 1, klon
2177	za = za + airephy(i)/FLOAT(klon)
2178	zx_t = zx_t + (rain_con(i)+
2179	. snow_con(i))*airephy(i)/FLOAT(klon)
2180	ENDDO
2181	zx_t = zx_t/za*dtime
2182	WRITE(lunout,*)"Precip=", zx_t
2183	ENDIF
2184	IF (zx_ajustq) THEN
2185	DO i = 1, klon
2186	z_apres(i) = 0.0
2187	ENDDO
2188	DO k = 1, klev
2189	DO i = 1, klon
2190	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k))
2191	. *(paprs(i,k)-paprs(i,k+1))/RG
2192	ENDDO
2193	ENDDO
2194	DO i = 1, klon
2195	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime)
2196	. /z_apres(i)
2197	ENDDO
2198	DO k = 1, klev
2199	DO i = 1, klon
2200	IF (z_factor(i).GT.(1.0+1.0E-08) .OR.
2201	. z_factor(i).LT.(1.0-1.0E-08)) THEN
2202	q_seri(i,k) = q_seri(i,k) * z_factor(i)
2203	ENDIF
2204	ENDDO
2205	ENDDO
2206	ENDIF
2207	zx_ajustq=.FALSE.
2208
2209	c
2210	c=============================================================================
2211	cRR:Evolution de la poche froide: on ne fait pas de separation wake/env
2212	cpour la couche limite diffuse pour l instant
2213	c
2214	if (iflag_wake.eq.1) then
2215	DO k=1,klev
2216	DO i=1,klon
2217	dt_dwn(i,k) = ftd(i,k)
2218	wdt_PBL(i,k) = 0.
2219	dq_dwn(i,k) = fqd(i,k)
2220	wdq_PBL(i,k) = 0.
2221	M_dwn(i,k) = dnwd0(i,k)
2222	M_up(i,k) = upwd(i,k)
2223	dt_a(i,k) = d_t_con(i,k)/dtime - ftd(i,k)
2224	udt_PBL(i,k) = 0.
2225	dq_a(i,k) = d_q_con(i,k)/dtime - fqd(i,k)
2226	udq_PBL(i,k) = 0.
2227	ENDDO
2228	ENDDO
2229	c
2230	ccalcul caracteristiques de la poche froide
2231	call calWAKE (paprs,pplay,dtime
2232	: ,t_seri,q_seri,omega
2233	: ,dt_dwn,dq_dwn,M_dwn,M_up
2234	: ,dt_a,dq_a,sigd
2235	: ,wdt_PBL,wdq_PBL
2236	: ,udt_PBL,udq_PBL
2237	o ,wake_deltat,wake_deltaq,wake_dth
2238	o ,wake_h,wake_s,wake_dens
2239	o ,wake_pe,wake_fip,wake_gfl
2240	o ,dt_wake,dq_wake
2241	o ,wake_k, t_undi,q_undi
2242	o ,wake_omgbdth,wake_dp_omgb
2243	o ,wake_dtKE,wake_dqKE
2244	o ,wake_dtPBL,wake_dqPBL
2245	o ,wake_omg,wake_dp_deltomg
2246	o ,wake_spread,wake_Cstar,wake_d_deltat_gw
2247	o ,wake_ddeltat,wake_ddeltaq)
2248	c
2249	!-----------------------------------------------------------------------------------------
2250	! ajout des tendances des poches froides
2251	! Faire rapidement disparaitre l'ancien dt_wake pour garder un d_t_wake
2252	! coherent avec les autres d_t_...
2253	d_t_wake(:,:)=dt_wake(:,:)*dtime
2254	d_q_wake(:,:)=dq_wake(:,:)*dtime
2255	CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,'wake')
2256	!-----------------------------------------------------------------------------------------
2257
2258	endif
2259	c print*,'apres callwake iflag_cldcon=', iflag_cldcon
2260	c
2261	c===================================================================
2262	c Convection seche (thermiques ou ajustement)
2263	c===================================================================
2264	c
2265	call stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri
2266	s ,seuil_inversion,weak_inversion,dthmin)
2267
2268
2269
2270	d_t_ajsb(:,:)=0.
2271	d_q_ajsb(:,:)=0.
2272	d_t_ajs(:,:)=0.
2273	d_u_ajs(:,:)=0.
2274	d_v_ajs(:,:)=0.
2275	d_q_ajs(:,:)=0.
2276	clwcon0th(:,:)=0.
2277	c
2278	fm_therm(:,:)=0.
2279	entr_therm(:,:)=0.
2280	detr_therm(:,:)=0.
2281	c
2282	IF(prt_level>9)WRITE(lunout,*)
2283	. 'AVANT LA CONVECTION SECHE , iflag_thermals='
2284	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2285	if(iflag_thermals.lt.0) then
2286	c Rien
2287	c ====
2288	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
2289
2290
2291	else
2292
2293	c Thermiques
2294	c ==========
2295	IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals='
2296	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2297
2298
2299	if (iflag_thermals.gt.1) then
2300	call calltherm(pdtphys
2301	s ,pplay,paprs,pphi,weak_inversion
2302	s ,u_seri,v_seri,t_seri,q_seri,zqsat,debut
2303	s ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
2304	s ,fm_therm,entr_therm,detr_therm
2305	s ,zqasc,clwcon0th,lmax_th,ratqscth
2306	s ,ratqsdiff,zqsatth
2307	con rajoute ale et alp, et les caracteristiques de la couche alim
2308	s ,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca)
2309	endif
2310
2311
2312	c Ajustement sec
2313	c ==============
2314
2315	! Dans le cas où on active les thermiques, on fait partir l'ajustement
2316	! a partir du sommet des thermiques.
2317	! Dans le cas contraire, on demarre au niveau 1.
2318
2319	if (iflag_thermals.ge.13.or.iflag_thermals.eq.0) then
2320
2321	if(iflag_thermals.eq.0) then
2322	IF(prt_level>9)WRITE(lunout,*)'ajsec'
2323	limbas(:)=1
2324	else
2325	limbas(:)=lmax_th(:)
2326	endif
2327
2328	! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
2329	! pour des test de convergence numerique.
2330	! Le nouveau ajsec est a priori mieux, meme pour le cas
2331	! iflag_thermals = 0 (l'ancienne version peut faire des tendances
2332	! non nulles numeriquement pour des mailles non concernees.
2333
2334	if (iflag_thermals.eq.0) then
2335	CALL ajsec_convV2(paprs, pplay, t_seri,q_seri
2336	s , d_t_ajsb, d_q_ajsb)
2337	else
2338	CALL ajsec(paprs, pplay, t_seri,q_seri,limbas
2339	s , d_t_ajsb, d_q_ajsb)
2340	endif
2341
2342	!-----------------------------------------------------------------------------------------
2343	! ajout des tendances de l'ajustement sec ou des thermiques
2344	CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,'ajsb')
2345	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
2346	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
2347
2348	!-----------------------------------------------------------------------------------------
2349
2350	endif
2351
2352	endif
2353	c
2354	c===================================================================
2355	cIM
2356	IF (ip_ebil_phy.ge.2) THEN
2357	ztit='after dry_adjust'
2358	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2359	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2360	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2361	END IF
2362
2363
2364	c-------------------------------------------------------------------------
2365	c Caclul des ratqs
2366	c-------------------------------------------------------------------------
2367
2368	c print*,'calcul des ratqs'
2369	c ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
2370	c ----------------
2371	c on ecrase le tableau ratqsc calcule par clouds_gno
2372	if (iflag_cldcon.eq.1) then
2373	do k=1,klev
2374	do i=1,klon
2375	if(ptconv(i,k)) then
2376	ratqsc(i,k)=ratqsbas
2377	s +fact_cldcon*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
2378	else
2379	ratqsc(i,k)=0.
2380	endif
2381	enddo
2382	enddo
2383
2384	c-----------------------------------------------------------------------
2385	c par nversion de la fonction log normale
2386	c-----------------------------------------------------------------------
2387	else if (iflag_cldcon.eq.4) then
2388	ptconvth(:,:)=.false.
2389	ratqsc(:,:)=0.
2390	if(prt_level.ge.9) print*,'avant clouds_gno thermique'
2391	call clouds_gno
2392	s (klon,klev,q_seri,zqsat,clwcon0th,ptconvth,ratqsc,rnebcon0th)
2393	if(prt_level.ge.9) print*,' CLOUDS_GNO OK'
2394
2395	c-----------------------------------------------------------------------
2396	c par calcul direct de l'ecart-type
2397	c-----------------------------------------------------------------------
2398
2399	else if (iflag_cldcon>=5) then
2400	wmax_th(:)=0.
2401	zmax_th(:)=0.
2402	do k=1,klev
2403	do i=1,klon
2404	wmax_th(i)=max(wmax_th(i),zw2(i,k))
2405	if (detr_therm(i,k).gt.0.) zmax_th(i)=pphi(i,k)/rg
2406	enddo
2407	enddo
2408	tau_overturning_th(:)=zmax_th(:)/max(0.5*wmax_th(:),0.1)
2409	print*,'TAU TH OK ',tau_overturning_th(1),detr_therm(1,3)
2410
2411	c On impose que l'air autour de la fraction couverte par le thermique
2412	c plus son air detraine durant tau_overturning_th soit superieur
2413	c a 0.1 q_seri
2414	zz=0.1
2415	do k=1,klev
2416	do i=1,klon
2417	lambda_th(i,k)=0.5*(fraca(i,k)+fraca(i,k+1))+
2418	s tau_overturning_th(i)*detr_therm(i,k)
2419	s *rg/(paprs(i,k)-paprs(i,k+1))
2420	znum=(1.-zz)*q_seri(i,k)
2421	zden=zqasc(i,k)-zz*q_seri(i,k)
2422	if (znum-lambda_th(i,k)*zden<0.) lambda_th(i,k)=znum/zden
2423	lambda_th(i,k)=min(lambda_th(i,k),0.9)
2424	enddo
2425	enddo
2426
2427	if(iflag_cldcon==5) then
2428	do k=1,klev
2429	do i=1,klon
2430	ratqsc(i,k)=sqrt(lambda_th(i,k)/(1.-lambda_th(i,k)))*
2431	s abs((zqasc(i,k)-q_seri(i,k))/q_seri(i,k))
2432	enddo
2433	enddo
2434	else if(iflag_cldcon==6) then
2435	do k=1,klev
2436	do i=1,klon
2437	ratqsc(i,k)=sqrt(lambda_th(i,k))*
2438	s (zqasc(i,k)-q_seri(i,k))/q_seri(i,k)
2439	enddo
2440	enddo
2441	endif
2442
2443	endif
2444
2445	c ratqs stables
2446	c -------------
2447
2448	if (iflag_ratqs.eq.0) then
2449
2450	! Le cas iflag_ratqs=0 correspond a la version IPCC 2005 du modele.
2451	do k=1,klev
2452	do i=1, klon
2453	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2454	s min((paprs(i,1)-pplay(i,k))/(paprs(i,1)-30000.),1.)
2455	enddo
2456	enddo
2457
2458	! Pour iflag_ratqs=1 ou 2, le ratqs est constant au dessus de
2459	! 300 hPa (ratqshaut), varie lineariement en fonction de la pression
2460	! entre 600 et 300 hPa et est soit constant (ratqsbas) pour iflag_ratqs=1
2461	! soit lineaire (entre 0 a la surface et ratqsbas) pour iflag_ratqs=2
2462	! Il s'agit de differents tests dans la phase de reglage du modele
2463	! avec thermiques.
2464
2465	else if (iflag_ratqs.eq.1) then
2466
2467	do k=1,klev
2468	do i=1, klon
2469	if (pplay(i,k).ge.60000.) then
2470	ratqss(i,k)=ratqsbas
2471	else if ((pplay(i,k).ge.30000.).and.
2472	s (pplay(i,k).lt.60000.)) then
2473	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2474	s (60000.-pplay(i,k))/(60000.-30000.)
2475	else
2476	ratqss(i,k)=ratqshaut
2477	endif
2478	enddo
2479	enddo
2480
2481	else
2482
2483	do k=1,klev
2484	do i=1, klon
2485	if (pplay(i,k).ge.60000.) then
2486	ratqss(i,k)=ratqsbas
2487	s *(paprs(i,1)-pplay(i,k))/(paprs(i,1)-60000.)
2488	else if ((pplay(i,k).ge.30000.).and.
2489	s (pplay(i,k).lt.60000.)) then
2490	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2491	s (60000.-pplay(i,k))/(60000.-30000.)
2492	else
2493	ratqss(i,k)=ratqshaut
2494	endif
2495	enddo
2496	enddo
2497	endif
2498
2499
2500
2501
2502	c ratqs final
2503	c -----------
2504
2505	if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2
2506	s .or.iflag_cldcon.ge.4) then
2507
2508	! On ajoute une constante au ratqsc*2 pour tenir compte de
2509	! fluctuations turbulentes de petite echelle
2510
2511	do k=1,klev
2512	do i=1,klon
2513	if ((fm_therm(i,k).gt.1.e-10)) then
2514	ratqsc(i,k)=sqrt(ratqsc(i,k)2+0.052)
2515	endif
2516	enddo
2517	enddo
2518
2519	! les ratqs sont une conbinaison de ratqss et ratqsc
2520	! 1800s, en dur pour le moment, est le temps de
2521	! relaxation des ratqs
2522
2523	facteur=exp(-pdtphys/1800.)
2524
2525	print*,'WARNING ratqs a revoir '
2526	ratqs(:,:)=ratqsc(:,:)(1.-facteur)+ratqs(:,:)facteur
2527	ratqs(:,:)=sqrt(ratqs(:,:)2+ratqss(:,:)2)
2528
2529	else
2530	! on ne prend que le ratqs stable pour fisrtilp
2531	ratqs(:,:)=ratqss(:,:)
2532	endif
2533
2534
2535	c
2536	c Appeler le processus de condensation a grande echelle
2537	c et le processus de precipitation
2538	c-------------------------------------------------------------------------
2539	CALL fisrtilp(dtime,paprs,pplay,
2540	. t_seri, q_seri,ptconv,ratqs,
2541	. d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq,
2542	. rain_lsc, snow_lsc,
2543	. pfrac_impa, pfrac_nucl, pfrac_1nucl,
2544	. frac_impa, frac_nucl,
2545	. prfl, psfl, rhcl)
2546
2547	WHERE (rain_lsc < 0) rain_lsc = 0.
2548	WHERE (snow_lsc < 0) snow_lsc = 0.
2549	!-----------------------------------------------------------------------------------------
2550	! ajout des tendances de la diffusion turbulente
2551	CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,'lsc')
2552	!-----------------------------------------------------------------------------------------
2553	DO k = 1, klev
2554	DO i = 1, klon
2555	cldfra(i,k) = rneb(i,k)
2556	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
2557	ENDDO
2558	ENDDO
2559	IF (check) THEN
2560	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2561	WRITE(lunout,*)"apresilp=", za
2562	zx_t = 0.0
2563	za = 0.0
2564	DO i = 1, klon
2565	za = za + airephy(i)/FLOAT(klon)
2566	zx_t = zx_t + (rain_lsc(i)
2567	. + snow_lsc(i))*airephy(i)/FLOAT(klon)
2568	ENDDO
2569	zx_t = zx_t/za*dtime
2570	WRITE(lunout,*)"Precip=", zx_t
2571	ENDIF
2572	cIM
2573	IF (ip_ebil_phy.ge.2) THEN
2574	ztit='after fisrt'
2575	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2576	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2577	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2578	call diagphy(airephy,ztit,ip_ebil_phy
2579	e , zero_v, zero_v, zero_v, zero_v, zero_v
2580	e , zero_v, rain_lsc, snow_lsc, ztsol
2581	e , d_h_vcol, d_qt, d_ec
2582	s , fs_bound, fq_bound )
2583	END IF
2584
2585	if (mydebug) then
2586	call writefield_phy('u_seri',u_seri,llm)
2587	call writefield_phy('v_seri',v_seri,llm)
2588	call writefield_phy('t_seri',t_seri,llm)
2589	call writefield_phy('q_seri',q_seri,llm)
2590	endif
2591
2592	c
2593	c-------------------------------------------------------------------
2594	c PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
2595	c-------------------------------------------------------------------
2596
2597	c 1. NUAGES CONVECTIFS
2598	c
2599	cIM cf FH
2600	c IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
2601	IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
2602	snow_tiedtke=0.
2603	c print*,'avant calcul de la pseudo precip '
2604	c print*,'iflag_cldcon',iflag_cldcon
2605	if (iflag_cldcon.eq.-1) then
2606	rain_tiedtke=rain_con
2607	else
2608	c print*,'calcul de la pseudo precip '
2609	rain_tiedtke=0.
2610	c print*,'calcul de la pseudo precip 0'
2611	do k=1,klev
2612	do i=1,klon
2613	if (d_q_con(i,k).lt.0.) then
2614	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys
2615	s *(paprs(i,k)-paprs(i,k+1))/rg
2616	endif
2617	enddo
2618	enddo
2619	endif
2620	c
2621	c call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
2622	c
2623
2624	c Nuages diagnostiques pour Tiedtke
2625	CALL diagcld1(paprs,pplay,
2626	cIM cf FH . rain_con,snow_con,ibas_con,itop_con,
2627	. rain_tiedtke,snow_tiedtke,ibas_con,itop_con,
2628	. diafra,dialiq)
2629	DO k = 1, klev
2630	DO i = 1, klon
2631	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2632	cldliq(i,k) = dialiq(i,k)
2633	cldfra(i,k) = diafra(i,k)
2634	ENDIF
2635	ENDDO
2636	ENDDO
2637
2638	ELSE IF (iflag_cldcon.ge.3) THEN
2639	c On prend pour les nuages convectifs le max du calcul de la
2640	c convection et du calcul du pas de temps precedent diminue d'un facteur
2641	c facttemps
2642	facteur = pdtphys *facttemps
2643	do k=1,klev
2644	do i=1,klon
2645	rnebcon(i,k)=rnebcon(i,k)*facteur
2646	if (rnebcon0(i,k)clwcon0(i,k).gt.rnebcon(i,k)clwcon(i,k))
2647	s then
2648	rnebcon(i,k)=rnebcon0(i,k)
2649	clwcon(i,k)=clwcon0(i,k)
2650	endif
2651	enddo
2652	enddo
2653
2654	c
2655	cjq - introduce the aerosol direct and first indirect radiative forcings
2656	cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
2657	IF (ok_ade.OR.ok_aie) THEN
2658	IF ( .NOT. aerosol_couple ) THEN
2659	! Get sulfate aerosol distribution
2660	CALL readsulfate(rjourvrai, debut, sulfate)
2661	CALL readsulfate_preind(rjourvrai, debut, sulfate_pi)
2662
2663	! Calculate aerosol optical properties (Olivier Boucher)
2664	CALL aeropt(pplay, paprs, t_seri, sulfate, rhcl,
2665	. tau_ae, piz_ae, cg_ae, aerindex)
2666	ENDIF
2667	ELSE
2668	tau_ae(:,:,:)=0.0
2669	piz_ae(:,:,:)=0.0
2670	cg_ae(:,:,:)=0.0
2671	ENDIF
2672
2673	cIM calcul nuages par le simulateur ISCCP
2674	c
2675	#ifdef histISCCP
2676	IF (ok_isccp) THEN
2677	c
2678	cIM lecture invtau, tautab des fichiers formattes
2679	c
2680	IF (debut) THEN
2681	c$OMP MASTER
2682	c
2683	open(99,file='tautab.formatted', FORM='FORMATTED')
2684	read(99,'(f30.20)') tautab_omp
2685	close(99)
2686	c
2687	open(99,file='invtau.formatted',form='FORMATTED')
2688	read(99,'(i10)') invtau_omp
2689
2690	c print*,'calcul_simulISCCP invtau_omp',invtau_omp
2691	c write(6,'(a,8i10)') 'invtau_omp',(invtau_omp(i),i=1,100)
2692
2693	close(99)
2694	c$OMP END MASTER
2695	c$OMP BARRIER
2696	tautab=tautab_omp
2697	invtau=invtau_omp
2698	c
2699	ENDIF !debut
2700	c
2701	cIM appel simulateur toutes les NINT(freq_ISCCP/dtime) heures
2702	IF (MOD(itap,NINT(freq_ISCCP/dtime)).EQ.0) THEN
2703	#include "calcul_simulISCCP.h"
2704	ENDIF !(MOD(itap,NINT(freq_ISCCP/dtime))
2705	ENDIF !ok_isccp
2706	#endif
2707
2708	c On prend la somme des fractions nuageuses et des contenus en eau
2709	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
2710	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
2711
2712	ENDIF
2713	c
2714	c 2. NUAGES STARTIFORMES
2715	c
2716	IF (ok_stratus) THEN
2717	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
2718	DO k = 1, klev
2719	DO i = 1, klon
2720	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2721	cldliq(i,k) = dialiq(i,k)
2722	cldfra(i,k) = diafra(i,k)
2723	ENDIF
2724	ENDDO
2725	ENDDO
2726	ENDIF
2727	c
2728	c Precipitation totale
2729	c
2730	DO i = 1, klon
2731	rain_fall(i) = rain_con(i) + rain_lsc(i)
2732	snow_fall(i) = snow_con(i) + snow_lsc(i)
2733	ENDDO
2734	cIM
2735	IF (ip_ebil_phy.ge.2) THEN
2736	ztit="after diagcld"
2737	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2738	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2739	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2740	END IF
2741	c
2742	c Calculer l'humidite relative pour diagnostique
2743	c
2744	DO k = 1, klev
2745	DO i = 1, klon
2746	zx_t = t_seri(i,k)
2747	IF (thermcep) THEN
2748	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2749	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2750	zx_qs = MIN(0.5,zx_qs)
2751	zcor = 1./(1.-retv*zx_qs)
2752	zx_qs = zx_qs*zcor
2753	ELSE
2754	IF (zx_t.LT.t_coup) THEN
2755	zx_qs = qsats(zx_t)/pplay(i,k)
2756	ELSE
2757	zx_qs = qsatl(zx_t)/pplay(i,k)
2758	ENDIF
2759	ENDIF
2760	zx_rh(i,k) = q_seri(i,k)/zx_qs
2761	zqsat(i,k)=zx_qs
2762	ENDDO
2763	ENDDO
2764
2765	cIM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
2766	c equivalente a 2m (tpote) pour diagnostique
2767	c
2768	DO i = 1, klon
2769	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
2770	IF (thermcep) THEN
2771	IF(zt2m(i).LT.RTT) then
2772	Lheat=RLSTT
2773	ELSE
2774	Lheat=RLVTT
2775	ENDIF
2776	ELSE
2777	IF (zt2m(i).LT.RTT) THEN
2778	Lheat=RLSTT
2779	ELSE
2780	Lheat=RLVTT
2781	ENDIF
2782	ENDIF
2783	tpote(i) = tpot(i)*
2784	. EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
2785	ENDDO
2786
2787	IF (config_inca /= 'none') THEN
2788	#ifdef INCA
2789	CALL VTe(VTphysiq)
2790	CALL VTb(VTinca)
2791	calday = FLOAT(julien) + gmtime
2792
2793	IF (config_inca == 'aero') THEN
2794	CALL AEROSOL_METEO_CALC(
2795	$ calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs,
2796	$ prfl,psfl,pctsrf,airephy,xjour,rlat,rlon,u10m,v10m)
2797	END IF
2798
2799	zxsnow_dummy(:) = 0.0
2800
2801	CALL chemhook_begin (calday,
2802	$ julien,
2803	$ gmtime,
2804	$ pctsrf(1,1),
2805	$ rlat,
2806	$ rlon,
2807	$ airephy,
2808	$ paprs,
2809	$ pplay,
2810	$ coefh,
2811	$ pphi,
2812	$ t_seri,
2813	$ u,
2814	$ v,
2815	$ wo,
2816	$ q_seri,
2817	$ zxtsol,
2818	$ zxsnow_dummy,
2819	$ solsw,
2820	$ albsol1,
2821	$ rain_fall,
2822	$ snow_fall,
2823	$ itop_con,
2824	$ ibas_con,
2825	$ cldfra,
2826	$ iim,
2827	$ jjm,
2828	$ tr_seri,
2829	$ ftsol,
2830	$ paprs,
2831	$ cdragh,
2832	$ cdragm,
2833	$ pctsrf,
2834	$ pdtphys,
2835	$ itap)
2836
2837	CALL VTe(VTinca)
2838	CALL VTb(VTphysiq)
2839	#endif
2840	END IF !config_inca /= 'none'
2841	c
2842	c Calculer les parametres optiques des nuages et quelques
2843	c parametres pour diagnostiques:
2844	c
2845
2846	IF (aerosol_couple) THEN
2847	sulfate(:,:) = ccm(:,:,1)
2848	sulfate_pi(:,:) = ccm(:,:,2)
2849	ENDIF
2850
2851	if (ok_newmicro) then
2852	CALL newmicro (paprs, pplay,ok_newmicro,
2853	. t_seri, cldliq, cldfra, cldtau, cldemi,
2854	. cldh, cldl, cldm, cldt, cldq,
2855	. flwp, fiwp, flwc, fiwc,
2856	e ok_aie,
2857	e sulfate, sulfate_pi,
2858	e bl95_b0, bl95_b1,
2859	s cldtaupi, re, fl)
2860	else
2861	CALL nuage (paprs, pplay,
2862	. t_seri, cldliq, cldfra, cldtau, cldemi,
2863	. cldh, cldl, cldm, cldt, cldq,
2864	e ok_aie,
2865	e sulfate, sulfate_pi,
2866	e bl95_b0, bl95_b1,
2867	s cldtaupi, re, fl)
2868
2869	endif
2870	c
2871	c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
2872	c
2873	IF (MOD(itaprad,radpas).EQ.0) THEN
2874
2875	DO i = 1, klon
2876	albsol1(i) = falb1(i,is_oce) * pctsrf(i,is_oce)
2877	. + falb1(i,is_lic) * pctsrf(i,is_lic)
2878	. + falb1(i,is_ter) * pctsrf(i,is_ter)
2879	. + falb1(i,is_sic) * pctsrf(i,is_sic)
2880	albsol2(i) = falb2(i,is_oce) * pctsrf(i,is_oce)
2881	. + falb2(i,is_lic) * pctsrf(i,is_lic)
2882	. + falb2(i,is_ter) * pctsrf(i,is_ter)
2883	. + falb2(i,is_sic) * pctsrf(i,is_sic)
2884	ENDDO
2885
2886	if (mydebug) then
2887	call writefield_phy('u_seri',u_seri,llm)
2888	call writefield_phy('v_seri',v_seri,llm)
2889	call writefield_phy('t_seri',t_seri,llm)
2890	call writefield_phy('q_seri',q_seri,llm)
2891	endif
2892
2893	IF (aerosol_couple) THEN
2894	#ifdef INCA
2895	CALL radlwsw_inca
2896	e (kdlon,kflev,dist, rmu0, fract, solaire,
2897	e paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri,
2898	e wo,
2899	e cldfra, cldemi, cldtau,
2900	s heat,heat0,cool,cool0,radsol,albpla,
2901	s topsw,toplw,solsw,sollw,
2902	s sollwdown,
2903	s topsw0,toplw0,solsw0,sollw0,
2904	s lwdn0, lwdn, lwup0, lwup,
2905	s swdn0, swdn, swup0, swup,
2906	e ok_ade, ok_aie,
2907	e tau_inca, piz_inca, cg_inca,
2908	s topswad_inca, solswad_inca,
2909	s topswad0_inca, solswad0_inca,
2910	s topsw_inca, topsw0_inca,
2911	s solsw_inca, solsw0_inca,
2912	e cldtaupi,
2913	s topswai_inca, solswai_inca)
2914	#endif
2915	ELSE
2916	CALL radlwsw ! nouveau rayonnement (compatible Arpege-IFS)
2917	e (dist, rmu0, fract,
2918	e paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri,
2919	e wo,
2920	e cldfra, cldemi, cldtau,
2921	s heat,heat0,cool,cool0,radsol,albpla,
2922	s topsw,toplw,solsw,sollw,
2923	s sollwdown,
2924	s topsw0,toplw0,solsw0,sollw0,
2925	s lwdn0, lwdn, lwup0, lwup,
2926	s swdn0, swdn, swup0, swup,
2927	e ok_ade, ok_aie, ! new for aerosol radiative effects
2928	e tau_ae, piz_ae, cg_ae, ! ="=
2929	s topswad, solswad, ! ="=
2930	e cldtaupi, ! ="=
2931	s topswai, solswai,zqsat,flwc,fiwc) ! ="=
2932	ENDIF
2933	itaprad = 0
2934	ENDIF
2935	itaprad = itaprad + 1
2936
2937	if (iflag_radia.eq.0) then
2938	print *,'--------------------------------------------------'
2939	print *,'>>>> ATTENTION rayonnement desactive pour ce cas'
2940	print *,'>>>> heat et cool mis a zero '
2941	print *,'--------------------------------------------------'
2942	heat=0.
2943	cool=0.
2944	endif
2945
2946	c
2947	c Ajouter la tendance des rayonnements (tous les pas)
2948	c
2949	DO k = 1, klev
2950	DO i = 1, klon
2951	t_seri(i,k) = t_seri(i,k)
2952	. + (heat(i,k)-cool(i,k)) * dtime/RDAY
2953	ENDDO
2954	ENDDO
2955	c
2956	if (mydebug) then
2957	call writefield_phy('u_seri',u_seri,llm)
2958	call writefield_phy('v_seri',v_seri,llm)
2959	call writefield_phy('t_seri',t_seri,llm)
2960	call writefield_phy('q_seri',q_seri,llm)
2961	endif
2962
2963	cIM
2964	IF (ip_ebil_phy.ge.2) THEN
2965	ztit='after rad'
2966	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2967	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2968	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2969	call diagphy(airephy,ztit,ip_ebil_phy
2970	e , topsw, toplw, solsw, sollw, zero_v
2971	e , zero_v, zero_v, zero_v, ztsol
2972	e , d_h_vcol, d_qt, d_ec
2973	s , fs_bound, fq_bound )
2974	END IF
2975	c
2976	c
2977	c Calculer l'hydrologie de la surface
2978	c
2979	c CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
2980	c . agesno, ftsol,fqsurf,fsnow, ruis)
2981	c
2982
2983	c
2984	c Calculer le bilan du sol et la derive de temperature (couplage)
2985	c
2986	DO i = 1, klon
2987	c bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
2988	c a la demande de JLD
2989	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
2990	ENDDO
2991	c
2992	cmoddeblott(jan95)
2993	c Appeler le programme de parametrisation de l'orographie
2994	c a l'echelle sous-maille:
2995	c
2996	IF (ok_orodr) THEN
2997	c
2998	c selection des points pour lesquels le shema est actif:
2999	igwd=0
3000	DO i=1,klon
3001	itest(i)=0
3002	c IF ((zstd(i).gt.10.0)) THEN
3003	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
3004	itest(i)=1
3005	igwd=igwd+1
3006	idx(igwd)=i
3007	ENDIF
3008	ENDDO
3009	c igwdim=MAX(1,igwd)
3010	c
3011	IF (ok_strato) THEN
3012
3013	CALL drag_noro_strato(klon,klev,dtime,paprs,pplay,
3014	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
3015	e igwd,idx,itest,
3016	e t_seri, u_seri, v_seri,
3017	s zulow, zvlow, zustrdr, zvstrdr,
3018	s d_t_oro, d_u_oro, d_v_oro)
3019
3020	ELSE
3021	CALL drag_noro(klon,klev,dtime,paprs,pplay,
3022	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
3023	e igwd,idx,itest,
3024	e t_seri, u_seri, v_seri,
3025	s zulow, zvlow, zustrdr, zvstrdr,
3026	s d_t_oro, d_u_oro, d_v_oro)
3027	ENDIF
3028	c
3029	c ajout des tendances
3030	!-----------------------------------------------------------------------------------------
3031	! ajout des tendances de la trainee de l'orographie
3032	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,'oro')
3033	!-----------------------------------------------------------------------------------------
3034	c
3035	ENDIF ! fin de test sur ok_orodr
3036	c
3037	if (mydebug) then
3038	call writefield_phy('u_seri',u_seri,llm)
3039	call writefield_phy('v_seri',v_seri,llm)
3040	call writefield_phy('t_seri',t_seri,llm)
3041	call writefield_phy('q_seri',q_seri,llm)
3042	endif
3043
3044	IF (ok_orolf) THEN
3045	c
3046	c selection des points pour lesquels le shema est actif:
3047	igwd=0
3048	DO i=1,klon
3049	itest(i)=0
3050	IF ((zpic(i)-zmea(i)).GT.100.) THEN
3051	itest(i)=1
3052	igwd=igwd+1
3053	idx(igwd)=i
3054	ENDIF
3055	ENDDO
3056	c igwdim=MAX(1,igwd)
3057	c
3058	IF (ok_strato) THEN
3059
3060	CALL lift_noro_strato(klon,klev,dtime,paprs,pplay,
3061	e rlat,zmea,zstd,zpic,zgam,zthe,zpic,zval,
3062	e igwd,idx,itest,
3063	e t_seri, u_seri, v_seri,
3064	s zulow, zvlow, zustrli, zvstrli,
3065	s d_t_lif, d_u_lif, d_v_lif )
3066
3067	ELSE
3068	CALL lift_noro(klon,klev,dtime,paprs,pplay,
3069	e rlat,zmea,zstd,zpic,
3070	e itest,
3071	e t_seri, u_seri, v_seri,
3072	s zulow, zvlow, zustrli, zvstrli,
3073	s d_t_lif, d_u_lif, d_v_lif)
3074	ENDIF
3075	c
3076	!-----------------------------------------------------------------------------------------
3077	! ajout des tendances de la portance de l'orographie
3078	CALL add_phys_tend(d_u_lif,d_v_lif,d_t_lif,dq0,dql0,'lif')
3079	!-----------------------------------------------------------------------------------------
3080	c
3081	ENDIF ! fin de test sur ok_orolf
3082	C HINES GWD PARAMETRIZATION
3083
3084	IF (ok_hines) then
3085
3086	CALL hines_gwd(klon,klev,dtime,paprs,pplay,
3087	i rlat,t_seri,u_seri,v_seri,
3088	o zustrhi,zvstrhi,
3089	o d_t_hin, d_u_hin, d_v_hin)
3090	c
3091	c ajout des tendances
3092	CALL add_phys_tend(d_u_hin,d_v_hin,d_t_hin,dq0,dql0,'lif')
3093
3094	ENDIF
3095	c
3096
3097	c
3098	cIM cf. FLott BEG
3099	C STRESS NECESSAIRES: TOUTE LA PHYSIQUE
3100
3101	if (mydebug) then
3102	call writefield_phy('u_seri',u_seri,llm)
3103	call writefield_phy('v_seri',v_seri,llm)
3104	call writefield_phy('t_seri',t_seri,llm)
3105	call writefield_phy('q_seri',q_seri,llm)
3106	endif
3107
3108	DO i = 1, klon
3109	zustrph(i)=0.
3110	zvstrph(i)=0.
3111	ENDDO
3112	DO k = 1, klev
3113	DO i = 1, klon
3114	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
3115	c (paprs(i,k)-paprs(i,k+1))/rg
3116	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
3117	c (paprs(i,k)-paprs(i,k+1))/rg
3118	ENDDO
3119	ENDDO
3120	c
3121	cIM calcul composantes axiales du moment angulaire et couple des montagnes
3122	c
3123	IF (is_sequential) THEN
3124
3125	CALL aaam_bud (27,klon,klev,rjourvrai,gmtime,
3126	C ra,rg,romega,
3127	C rlat,rlon,pphis,
3128	C zustrdr,zustrli,zustrph,
3129	C zvstrdr,zvstrli,zvstrph,
3130	C paprs,u,v,
3131	C aam, torsfc)
3132	ENDIF
3133	cIM cf. FLott END
3134	cIM
3135	IF (ip_ebil_phy.ge.2) THEN
3136	ztit='after orography'
3137	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
3138	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3139	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3140	END IF
3141	c
3142	c
3143	cAA
3144	cAA Installation de l'interface online-offline pour traceurs
3145	cAA
3146	c====================================================================
3147	c Calcul des tendances traceurs
3148	c====================================================================
3149	C
3150	IF (config_inca /= 'none') rnpb=.FALSE.
3151
3152	call phytrac ( rnpb,
3153	I itap,
3154	I julien,
3155	I gmtime,
3156	I debut,
3157	I lafin,
3158	I nlon,
3159	I nlev,
3160	I dtime,
3161	I u,
3162	I v,
3163	I t,
3164	I paprs,
3165	I pplay,
3166	I pmfu,
3167	I pmfd,
3168	I pen_u,
3169	I pde_u,
3170	I pen_d,
3171	I pde_d,
3172	I cdragh,
3173	I coefh,
3174	I fm_therm,
3175	I entr_therm,
3176	I u1,
3177	I v1,
3178	I ftsol,
3179	I pctsrf,
3180	I rlat,
3181	I frac_impa,
3182	I frac_nucl,
3183	I rlon,
3184	I presnivs,
3185	I pphis,
3186	I pphi,
3187	I albsol1,
3188	I qx(1,1,1),
3189	I rhcl,
3190	I cldfra,
3191	I rneb,
3192	I diafra,
3193	I cldliq,
3194	I itop_con,
3195	I ibas_con,
3196	I pmflxr,
3197	I pmflxs,
3198	I prfl,
3199	I psfl,
3200	I da,
3201	I phi,
3202	I mp,
3203	I upwd,
3204	I dnwd,
3205	I aerosol_couple,
3206	I flxmass_w,
3207	I tau_inca,
3208	I piz_inca,
3209	I cg_inca,
3210	I ccm,
3211	I rfname,
3212	O tr_seri)
3213
3214	IF (offline) THEN
3215
3216	print*,'Attention on met a 0 les thermiques pour phystoke'
3217	call phystokenc (
3218	I nlon,nlev,pdtphys,rlon,rlat,
3219	I t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
3220	I fm_therm,entr_therm,
3221	I cdragh,coefh,u1,v1,ftsol,pctsrf,
3222	I frac_impa, frac_nucl,
3223	I pphis,airephy,dtime,itap)
3224
3225
3226	ENDIF
3227
3228	c
3229	c Calculer le transport de l'eau et de l'energie (diagnostique)
3230	c
3231	CALL transp (paprs,zxtsol,
3232	e t_seri, q_seri, u_seri, v_seri, zphi,
3233	s ve, vq, ue, uq)
3234	c
3235	cIM global posePB BEG
3236	IF(1.EQ.0) THEN
3237	c
3238	CALL transp_lay (paprs,zxtsol,
3239	e t_seri, q_seri, u_seri, v_seri, zphi,
3240	s ve_lay, vq_lay, ue_lay, uq_lay)
3241	c
3242	ENDIF !(1.EQ.0) THEN
3243	cIM global posePB END
3244	c Accumuler les variables a stocker dans les fichiers histoire:
3245	c
3246	c+jld ec_conser
3247	DO k = 1, klev
3248	DO i = 1, klon
3249	ZRCPD = RCPD(1.0+RVTMP2q_seri(i,k))
3250	d_t_ec(i,k)=0.5/ZRCPD
3251	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
3252	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
3253	d_t_ec(i,k) = d_t_ec(i,k)/dtime
3254	END DO
3255	END DO
3256	c-jld ec_conser
3257	cIM
3258	IF (ip_ebil_phy.ge.1) THEN
3259	ztit='after physic'
3260	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
3261	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3262	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3263	C Comme les tendances de la physique sont ajoute dans la dynamique,
3264	C on devrait avoir que la variation d'entalpie par la dynamique
3265	C est egale a la variation de la physique au pas de temps precedent.
3266	C Donc la somme de ces 2 variations devrait etre nulle.
3267	call diagphy(airephy,ztit,ip_ebil_phy
3268	e , topsw, toplw, solsw, sollw, sens
3269	e , evap, rain_fall, snow_fall, ztsol
3270	e , d_h_vcol, d_qt, d_ec
3271	s , fs_bound, fq_bound )
3272	C
3273	d_h_vcol_phy=d_h_vcol
3274	C
3275	END IF
3276	C
3277	c=======================================================================
3278	c SORTIES
3279	c=======================================================================
3280
3281	cIM Interpolation sur les niveaux de pression du NMC
3282	c -------------------------------------------------
3283	c
3284	#include "calcul_STDlev.h"
3285	twriteSTD(:,:,1)=tsumSTD(:,:,2)
3286	qwriteSTD(:,:,1)=qsumSTD(:,:,2)
3287	rhwriteSTD(:,:,1)=rhsumSTD(:,:,2)
3288	phiwriteSTD(:,:,1)=phisumSTD(:,:,2)
3289	uwriteSTD(:,:,1)=usumSTD(:,:,2)
3290	vwriteSTD(:,:,1)=vsumSTD(:,:,2)
3291	wwriteSTD(:,:,1)=wsumSTD(:,:,2)
3292
3293	twriteSTD(:,:,2)=tsumSTD(:,:,1)
3294	qwriteSTD(:,:,2)=qsumSTD(:,:,1)
3295	rhwriteSTD(:,:,2)=rhsumSTD(:,:,1)
3296	phiwriteSTD(:,:,2)=phisumSTD(:,:,1)
3297	uwriteSTD(:,:,2)=usumSTD(:,:,1)
3298	vwriteSTD(:,:,2)=vsumSTD(:,:,1)
3299	wwriteSTD(:,:,2)=wsumSTD(:,:,1)
3300
3301	twriteSTD(:,:,3)=tlevSTD(:,:)
3302	qwriteSTD(:,:,3)=qlevSTD(:,:)
3303	rhwriteSTD(:,:,3)=rhlevSTD(:,:)
3304	phiwriteSTD(:,:,3)=philevSTD(:,:)
3305	uwriteSTD(:,:,3)=ulevSTD(:,:)
3306	vwriteSTD(:,:,3)=vlevSTD(:,:)
3307	wwriteSTD(:,:,3)=wlevSTD(:,:)
3308
3309	twriteSTD(:,:,4)=tlevSTD(:,:)
3310	qwriteSTD(:,:,4)=qlevSTD(:,:)
3311	rhwriteSTD(:,:,4)=rhlevSTD(:,:)
3312	phiwriteSTD(:,:,4)=philevSTD(:,:)
3313	uwriteSTD(:,:,4)=ulevSTD(:,:)
3314	vwriteSTD(:,:,4)=vlevSTD(:,:)
3315	wwriteSTD(:,:,4)=wlevSTD(:,:)
3316	c
3317	c slp sea level pressure
3318	slp(:) = paprs(:,1)exp(pphis(:)/(RDt_seri(:,1)))
3319	c
3320	ccc prw = eau precipitable
3321	DO i = 1, klon
3322	prw(i) = 0.
3323	DO k = 1, klev
3324	prw(i) = prw(i) +
3325	. q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
3326	ENDDO
3327	ENDDO
3328	c
3329	cIM initialisation + calculs divers diag AMIP2
3330	c
3331	#include "calcul_divers.h"
3332	c
3333	IF (config_inca /= 'none') THEN
3334	#ifdef INCA
3335	CALL VTe(VTphysiq)
3336	CALL VTb(VTinca)
3337
3338	CALL chemhook_end (calday,
3339	$ dtime,
3340	$ pplay,
3341	$ t_seri,
3342	$ tr_seri,
3343	$ nbtr,
3344	$ paprs,
3345	$ q_seri,
3346	$ annee_ref,
3347	$ day_ini,
3348	$ airephy,
3349	$ xjour,
3350	$ pphi,
3351	$ pphis,
3352	$ zx_rh)
3353
3354	CALL VTe(VTinca)
3355	CALL VTb(VTphysiq)
3356	#endif
3357	END IF
3358
3359	c=============================================================
3360	c
3361	c Convertir les incrementations en tendances
3362	c
3363	if (mydebug) then
3364	call writefield_phy('u_seri',u_seri,llm)
3365	call writefield_phy('v_seri',v_seri,llm)
3366	call writefield_phy('t_seri',t_seri,llm)
3367	call writefield_phy('q_seri',q_seri,llm)
3368	endif
3369
3370	DO k = 1, klev
3371	DO i = 1, klon
3372	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
3373	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
3374	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
3375	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
3376	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
3377	ENDDO
3378	ENDDO
3379	c
3380	IF (nqtot.GE.3) THEN
3381	DO iq = 3, nqtot
3382	DO k = 1, klev
3383	DO i = 1, klon
3384	d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
3385	ENDDO
3386	ENDDO
3387	ENDDO
3388	ENDIF
3389	c
3390	cIM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
3391	cIM global posePB#include "write_bilKP_ins.h"
3392	cIM global posePB#include "write_bilKP_ave.h"
3393	c
3394	c Sauvegarder les valeurs de t et q a la fin de la physique:
3395	c
3396	DO k = 1, klev
3397	DO i = 1, klon
3398	u_ancien(i,k) = u_seri(i,k)
3399	v_ancien(i,k) = v_seri(i,k)
3400	t_ancien(i,k) = t_seri(i,k)
3401	q_ancien(i,k) = q_seri(i,k)
3402	ENDDO
3403	ENDDO
3404	c
3405	!==========================================================================
3406	! Sorties des tendances pour un point particulier
3407	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
3408	! pour le debug
3409	! La valeur de igout est attribuee plus haut dans le programme
3410	!==========================================================================
3411
3412	if (prt_level.ge.1) then
3413	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
3414	write(lunout,*)
3415	s 'nlon,nlev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys pct tlos'
3416	write(lunout,*)
3417	s nlon,nlev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys,
3418	s pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce),
3419	s pctsrf(igout,is_sic)
3420	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
3421	do k=1,nlev
3422	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k),
3423	s d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k),
3424	s d_t_eva(igout,k)
3425	enddo
3426	write(lunout,*) 'cool,heat'
3427	do k=1,nlev
3428	write(lunout,*) cool(igout,k),heat(igout,k)
3429	enddo
3430
3431	write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
3432	do k=1,nlev
3433	write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k),
3434	s d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
3435	enddo
3436
3437	write(lunout,*) 'd_ps ',d_ps(igout)
3438	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
3439	do k=1,nlev
3440	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k),
3441	s d_qx(igout,k,1),d_qx(igout,k,2)
3442	enddo
3443	endif
3444
3445	!==========================================================================
3446
3447	c============================================================
3448	c Calcul de la temperature potentielle
3449	c============================================================
3450	DO k = 1, klev
3451	DO i = 1, klon
3452	theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
3453	ENDDO
3454	ENDDO
3455	c
3456
3457	c 22.03.04 BEG
3458	c=============================================================
3459	c Ecriture des sorties
3460	c=============================================================
3461	#ifdef CPP_IOIPSL
3462
3463	c Recupere des varibles calcule dans differents modules
3464	c pour ecriture dans histxxx.nc
3465
3466	! Get some variables from module fonte_neige_mod
3467	CALL fonte_neige_get_vars(pctsrf,
3468	. zxfqcalving, zxfqfonte, zxffonte)
3469
3470
3471	#include "phys_output_write.h"
3472
3473	#ifdef histISCCP
3474	#include "write_histISCCP.h"
3475	#endif
3476
3477	#ifdef histmthNMC
3478	#include "write_histmthNMC.h"
3479	#endif
3480
3481	#include "write_histday_seri.h"
3482
3483	#include "write_paramLMDZ_phy.h"
3484
3485	#endif
3486
3487	c 22.03.04 END
3488	c
3489	c====================================================================
3490	c Si c'est la fin, il faut conserver l'etat de redemarrage
3491	c====================================================================
3492	c
3493
3494
3495	IF (lafin) THEN
3496	itau_phy = itau_phy + itap
3497	CALL phyredem ("restartphy.nc")
3498	! open(97,form="unformatted",file="finbin")
3499	! write(97) u_seri,v_seri,t_seri,q_seri
3500	! close(97)
3501	ENDIF
3502
3503
3504	RETURN
3505	END
3506	FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
3507	IMPLICIT none
3508	c
3509	c Calculer et imprimer l'eau totale. A utiliser pour verifier
3510	c la conservation de l'eau
3511	c
3512	#include "YOMCST.h"
3513	INTEGER klon,klev
3514	REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
3515	REAL aire(klon)
3516	REAL qtotal, zx, qcheck
3517	INTEGER i, k
3518	c
3519	zx = 0.0
3520	DO i = 1, klon
3521	zx = zx + aire(i)
3522	ENDDO
3523	qtotal = 0.0
3524	DO k = 1, klev
3525	DO i = 1, klon
3526	qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i)
3527	. *(paprs(i,k)-paprs(i,k+1))/RG
3528	ENDDO
3529	ENDDO
3530	c
3531	qcheck = qtotal/zx
3532	c
3533	RETURN
3534	END
3535	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
3536	IMPLICIT none
3537	c
3538	c Tranformer une variable de la grille physique a
3539	c la grille d'ecriture
3540	c
3541	INTEGER nfield,nlon,iim,jjmp1, jjm
3542	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
3543	c
3544	INTEGER i, n, ig
3545	c
3546	jjm = jjmp1 - 1
3547	DO n = 1, nfield
3548	DO i=1,iim
3549	ecrit(i,n) = fi(1,n)
3550	ecrit(i+jjm*iim,n) = fi(nlon,n)
3551	ENDDO
3552	DO ig = 1, nlon - 2
3553	ecrit(iim+ig,n) = fi(1+ig,n)
3554	ENDDO
3555	ENDDO
3556	RETURN
3557	END

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

Download in other formats:

Original Format