Context Navigation

physiq.F @ 956

Last change on this file since 956 was 956, checked in by lmdzadmin, 16 years ago

Nettoyage du controle des parametres physiques. FH

Les parametres cycle_diurne, soil_model, new_oliq, ok_orodr, ok_orolf, ok_limitvrai, nbapp_rad et iflag_con
sont maintenant geres par la physique uniquement.
ecritphy est elimine.
dimphy.F90 et clesphys.h ne sont plus utilises par le code dynamique.
Le test academique obtenu en compilant avec
makegcm -p nophys gcm
fonctionne. FH
IM

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

File size: 110.1 KB

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

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

Download in other formats:

Original Format