Context Navigation

physiq.F @ 953

Last change on this file since 953 was 953, checked in by lmdzadmin, 16 years ago
On supprime variables pbase, ibas qui servent pas ; print sous clef JYG/FH IM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 109.9 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_AER
1077	#ifdef CPP_COUPLE
1078	! Aerosol optical properties by INCA model
1079	CHARACTER*4 :: rfname(9)
1080
1081	#endif
1082	#endif
1083	REAL aerindex(klon) ! POLDER aerosol index
1084
1085	! Parameters
1086	LOGICAL ok_ade, ok_aie ! Apply aerosol (in)direct effects or not
1087	REAL bl95_b0, bl95_b1 ! Parameter in Boucher and Lohmann (1995)
1088	cym
1089	SAVE ok_ade, ok_aie, bl95_b0, bl95_b1
1090	c$OMP THREADPRIVATE(ok_ade, ok_aie, bl95_b0, bl95_b1)
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_AER
1212	#ifdef CPP_COUPLE
1213	tau_inca(:,:,:,:) = 0.
1214	piz_inca(:,:,:,:) = 0.
1215	cg_inca(:,:,:,:) = 0.
1216	ccm(:,:,:) = 0.
1217	topswai_inca(:) = 0.
1218	topswad_inca(:) = 0.
1219	topswad0_inca(:) = 0.
1220	topsw_inca(:,:) = 0.
1221	topsw0_inca(:,:) = 0.
1222	solswai_inca(:) = 0.
1223	solswad_inca(:) = 0.
1224	solswad0_inca(:) = 0.
1225	solsw_inca(:,:) = 0.
1226	solsw0_inca(:,:) = 0.
1227	#endif
1228	#endif
1229	!rv
1230	!ACo
1231	rnebcon0(:,:) = 0.0
1232	clwcon0(:,:) = 0.0
1233	rnebcon(:,:) = 0.0
1234	clwcon(:,:) = 0.0
1235	! fin anne
1236
1237	cIM
1238	IF (ip_ebil_phy.ge.1) d_h_vcol_phy=0.
1239	c
1240	c appel a la lecture du run.def physique
1241	c
1242	call conf_phys(ocean, ok_veget, ok_journe, ok_mensuel,
1243	. ok_instan, ok_hf,
1244	. solarlong0,seuil_inversion,
1245	. fact_cldcon, facttemps,ok_newmicro,iflag_radia,
1246	. iflag_cldcon,iflag_ratqs,ratqsbas,ratqshaut,
1247	. ok_ade, ok_aie,
1248	. bl95_b0, bl95_b1,
1249	. iflag_thermals,nsplit_thermals,
1250	cnv flags pour la convection et les poches froides
1251	. iflag_coupl,iflag_clos,iflag_wake)
1252
1253	print*,'iflag_coupl,iflag_clos,iflag_wake',
1254	. iflag_coupl,iflag_clos,iflag_wake
1255
1256	c
1257	c
1258	c Initialiser les compteurs:
1259	c
1260	itap = 0
1261	itaprad = 0
1262
1263	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1264	!! Un petit travail à faire ici.
1265	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1266
1267	if (iflag_pbl>1) then
1268	PRINT*, "Using method MELLOR&YAMADA"
1269	endif
1270
1271	CALL phyetat0 ("startphy.nc",dtime,co2_ppm_etat0,solaire_etat0,
1272	. rlat,rlon,pctsrf, ftsol,
1273	. ocean, ok_veget,
1274	. falb1, falb2, rain_fall,snow_fall,
1275	. solsw, sollw,
1276	. radsol,clesphy0,
1277	. zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,tabcntr0,
1278	. t_ancien, q_ancien, ancien_ok, rnebcon, ratqs,clwcon,
1279	. pbl_tke, zmax0, f0, ema_work1, ema_work2)
1280
1281
1282
1283
1284	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1285
1286
1287	DO i=1,klon
1288	IF ( abs( pctsrf(i, is_ter) + pctsrf(i, is_lic) +
1289	$ pctsrf(i, is_oce) + pctsrf(i, is_sic) - 1.) .GT. EPSFRA)
1290	$ THEN
1291	WRITE(,)
1292	$ 'physiq apres lecture de restart: pb sous surface au point ',
1293	$ i, pctsrf(i, 1 : nbsrf)
1294	ENDIF
1295	ENDDO
1296
1297	radpas = NINT( 86400./dtime/nbapp_rad)
1298	c
1299	C on remet le calendrier a zero
1300	c
1301	IF (raz_date .eq. 1) THEN
1302	itau_phy = 0
1303	ENDIF
1304
1305	cIM cf. AM 081204 BEG
1306	PRINT*,'cycle_diurne3 =',cycle_diurne
1307	cIM cf. AM 081204 END
1308	c
1309	CALL printflag( tabcntr0,radpas,ok_journe,
1310	, ok_instan, ok_region )
1311	c
1312	IF (ABS(dtime-pdtphys).GT.0.001) THEN
1313	WRITE(lunout,*) 'Pas physique n est pas correct',dtime,
1314	. pdtphys
1315	abort_message='Pas physique n est pas correct '
1316	! call abort_gcm(modname,abort_message,1)
1317	dtime=pdtphys
1318	ENDIF
1319	IF (nlon .NE. klon) THEN
1320	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,
1321	. klon
1322	abort_message='nlon et klon ne sont pas coherents'
1323	call abort_gcm(modname,abort_message,1)
1324	ENDIF
1325	IF (nlev .NE. klev) THEN
1326	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,
1327	. klev
1328	abort_message='nlev et klev ne sont pas coherents'
1329	call abort_gcm(modname,abort_message,1)
1330	ENDIF
1331	c
1332	IF (dtime*FLOAT(radpas).GT.21600..AND.cycle_diurne) THEN
1333	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
1334	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
1335	abort_message='Nbre d appels au rayonnement insuffisant'
1336	call abort_gcm(modname,abort_message,1)
1337	ENDIF
1338	WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
1339	WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=",
1340	. ok_cvl
1341	c
1342	cKE43
1343	c Initialisation pour la convection de K.E. (sb):
1344	IF (iflag_con.GE.3) THEN
1345
1346	WRITE(lunout,)"** Convection de Kerry Emanuel 4.3 "
1347	WRITE(lunout,*)
1348	. "On va utiliser le melange convectif des traceurs qui"
1349	WRITE(lunout,*)"est calcule dans convect4.3"
1350	WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
1351
1352	DO i = 1, klon
1353	ema_cbmf(i) = 0.
1354	ema_pcb(i) = 0.
1355	ema_pct(i) = 0.
1356	ema_workcbmf(i) = 0.
1357	ENDDO
1358	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
1359	DO i = 1, klon
1360	ibas_con(i) = 1
1361	itop_con(i) = 1
1362	ENDDO
1363	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
1364	c===============================================================================
1365	cCR:04.12.07: initialisations poches froides
1366	c Controle de ALE et ALP pour la fermeture convective (jyg)
1367	if (iflag_wake.eq.1) then
1368	CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr
1369	s ,alp_bl_prescr, ale_bl_prescr)
1370	c 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
1371	c print*,'apres ini_wake iflag_cldcon=', iflag_cldcon
1372	endif
1373
1374	do i = 1,klon
1375	wake_s(i) = 0.
1376	wake_fip(i) = 0.
1377	wake_cstar(i) = 0.
1378	DO k=1,klev
1379	wake_deltat(i,k)=0.
1380	wake_deltaq(i,k)=0.
1381	ENDDO
1382	enddo
1383	c================================================================================
1384
1385	ENDIF
1386
1387	rugoro=0.
1388	c34EK
1389	IF (ok_orodr) THEN
1390
1391	rugoro=0.
1392
1393	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1394	! FH sans doute a enlever de finitivement ou, si on le garde, l'activer
1395	! justement quand ok_orodr = false.
1396	! ce rugoro est utilise par la couche limite et fait double emploi
1397	! avec les paramétrisations spécifiques de Francois Lott.
1398	! DO i=1,klon
1399	! rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1400	! ENDDO
1401	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1402
1403	CALL SUGWD(klon,klev,paprs,pplay)
1404	DO i=1,klon
1405	zuthe(i)=0.
1406	zvthe(i)=0.
1407	if(zstd(i).gt.10.)then
1408	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
1409	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
1410	endif
1411	ENDDO
1412	ENDIF
1413	c
1414	c
1415	lmt_pas = NINT(86400./dtime * 1.0) ! tous les jours
1416	WRITE(lunout,*)'La frequence de lecture surface est de ',
1417	. lmt_pas
1418	c
1419	cIM200505 ecrit_mth = NINT(86400./dtime *ecritphy) ! tous les ecritphy jours
1420	c IF (ok_mensuel) THEN
1421	c WRITE(lunout,*)'La frequence de sortie mensuelle est de ',
1422	c . ecrit_mth
1423	c ENDIF
1424	c ecrit_day = NINT(86400./dtime *1.0) ! tous les jours
1425	c IF (ok_journe) THEN
1426	c WRITE(lunout,*)'La frequence de sortie journaliere est de ',
1427	c . ecrit_day
1428	c ENDIF
1429	cIM 130904 BEG
1430	cIM 080205 ecrit_hf = 86400./dtime *0.25 ! toutes les 6h
1431	cIM 170305
1432	c ecrit_hf = 86400./dtime/12. ! toutes les 2h
1433	cIM 230305
1434	cIM200505 ecrit_hf = 86400./dtime *0.25 ! toutes les 6h
1435	c
1436	cIM200505 ecrit_hf2mth = ecrit_day/ecrit_hf*30
1437	c
1438	cIM200505 IF (ok_journe) THEN
1439	cIM200505 WRITE(lunout,*)'La frequence de sortie hf est de ',
1440	cIM200505 . ecrit_hf
1441	cIM200505 ENDIF
1442	cIM 130904 END
1443	ccc ecrit_ins = NINT(86400./dtime *0.5) ! 2 fois par jour
1444	ccc ecrit_ins = NINT(86400./dtime *0.25) ! 4 fois par jour
1445	c ecrit_ins = NINT(86400./dtime/48.) ! a chaque pas de temps ==> PB. dans time_counter pour 1mois
1446	c ecrit_ins = NINT(86400./dtime/12.) ! toutes les deux heures
1447	cIM200505 ecrit_ins = NINT(86400./dtime/8.) ! toutes les trois heures
1448	cIM200505 IF (ok_instan) THEN
1449	cIM200505 WRITE(lunout,*)'La frequence de sortie instant. est de ',
1450	cIM200505 . ecrit_ins
1451	cIM200505 ENDIF
1452	cIM200505 ecrit_reg = NINT(86400./dtime *0.25) ! 4 fois par jour
1453	cIM200505 IF (ok_region) THEN
1454	cIM200505 WRITE(lunout,*)'La frequence de sortie region est de ',
1455	cIM200505 . ecrit_reg
1456	cIM200505 ENDIF
1457	cIM 030306 BEG
1458	cIM ecrit_hf2mth = nombre de pas de temps de calcul de hf par mois apres lequel on ecrit
1459	cIM : ne pas modifier ecrit_hf2mth
1460	c
1461	cIM 250308bad guide ecrit_hf2mth = 30*1/ecrit_hf
1462	ecrit_hf2mth = ecrit_mth/ecrit_hf
1463	c ecrit_ins en secondes, chaque pas de temps de la physique
1464	ecrit_ins = dtime
1465	cIM on passe les frequences de jours en secondes : ecrit_ins, ecrit_hf, ecrit_day, ecrit_mth, ecrit_tra, ecrit_reg
1466	ecrit_hf = ecrit_hf * un_jour
1467	!IM
1468	IF(ecrit_day.LE.1.) THEN
1469	ecrit_day = ecrit_day * un_jour !en secondes
1470	ENDIF
1471	!IM
1472	ecrit_mth = ecrit_mth * un_jour
1473	ecrit_reg = ecrit_reg * un_jour
1474	ecrit_tra = ecrit_tra * un_jour
1475	ecrit_ISCCP = ecrit_ISCCP * un_jour
1476	c
1477	PRINT*,'physiq ecrit_ hf day mth reg tra ISCCP hf2mth',
1478	. ecrit_hf,ecrit_day,ecrit_mth,ecrit_reg,ecrit_tra,ecrit_ISCCP,
1479	. ecrit_hf2mth
1480	cIM 030306 END
1481
1482	capemaxcels = 't_max(X)'
1483	t2mincels = 't_min(X)'
1484	t2maxcels = 't_max(X)'
1485	tinst = 'inst(X)'
1486	tave = 'ave(X)'
1487	cIM cf. AM 081204 BEG
1488	write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
1489	cIM cf. AM 081204 END
1490	c
1491	c=============================================================
1492	c Initialisation des sorties
1493	c=============================================================
1494
1495	#ifdef CPP_IOIPSL
1496
1497	c Commente par abderrahmane 11 2 08
1498	c#ifdef histhf
1499	c#include "ini_histhf.h"
1500	c#endif
1501
1502	c#ifdef histday
1503	c#include "ini_histday.h"
1504	cIM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
1505	c#include "ini_bilKP_ins.h"
1506	c#include "ini_bilKP_ave.h"
1507	c#endif
1508
1509	c#ifdef histmth
1510	c#include "ini_histmth.h"
1511	c#endif
1512
1513	c#ifdef histins
1514	c#include "ini_histins.h"
1515	c#endif
1516
1517	call phys_output_open(jjmp1,nqmax,nlevSTD,clevSTD,nbteta,
1518	& ctetaSTD,dtime,presnivs,ok_veget,
1519	& ocean,iflag_pbl,ok_mensuel,ok_journe,
1520	& ok_hf,ok_instan,nid_files)
1521
1522	#ifdef histISCCP
1523	#include "ini_histISCCP.h"
1524	#endif
1525
1526	#ifdef histmthNMC
1527	#include "ini_histmthNMC.h"
1528	#endif
1529
1530	#include "ini_histday_seri.h"
1531
1532	#include "ini_paramLMDZ_phy.h"
1533
1534	#endif
1535
1536	cXXXPB Positionner date0 pour initialisation de ORCHIDEE
1537	date0 = zjulian
1538	C date0 = day_ini
1539	WRITE(,) 'physiq date0 : ',date0
1540	c
1541	c
1542	c
1543	c Prescrire l'ozone dans l'atmosphere
1544	c
1545	c
1546	cc DO i = 1, klon
1547	cc DO k = 1, klev
1548	cc CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1549	cc ENDDO
1550	cc ENDDO
1551	c
1552	#ifdef INCA
1553	call VTe(VTphysiq)
1554	call VTb(VTinca)
1555	iii = MOD(NINT(xjour),360)
1556	calday = FLOAT(iii) + gmtime
1557	WRITE(lunout,*) 'initial time ', xjour, calday
1558	#ifdef INCAINFO
1559	WRITE(lunout,*) 'Appel CHEMINI ...'
1560	#endif
1561	CALL chemini(
1562	$ rg,
1563	$ ra,
1564	$ airephy,
1565	$ rlat,
1566	$ rlon,
1567	$ presnivs,
1568	$ calday,
1569	$ klon,
1570	$ nqmax,
1571	$ pdtphys,
1572	$ annee_ref,
1573	$ day_ini)
1574	#ifdef INCAINFO
1575	WRITE(lunout,*) 'OK.'
1576	#endif
1577	call VTe(VTinca)
1578	call VTb(VTphysiq)
1579	#endif
1580	c
1581	ENDIF
1582	c
1583	c ************** Fin de IF ( debut ) *************
1584	c
1585
1586	! Tendances bidons pour les processus qui n'affectent pas certaines
1587	! variables.
1588	du0(:,:)=0.
1589	dv0(:,:)=0.
1590	dq0(:,:)=0.
1591	dql0(:,:)=0.
1592	c
1593	c Mettre a zero des variables de sortie (pour securite)
1594	c
1595	DO i = 1, klon
1596	d_ps(i) = 0.0
1597	ENDDO
1598	DO k = 1, klev
1599	DO i = 1, klon
1600	d_t(i,k) = 0.0
1601	d_u(i,k) = 0.0
1602	d_v(i,k) = 0.0
1603	ENDDO
1604	ENDDO
1605	DO iq = 1, nqmax
1606	DO k = 1, klev
1607	DO i = 1, klon
1608	d_qx(i,k,iq) = 0.0
1609	ENDDO
1610	ENDDO
1611	ENDDO
1612	da(:,:)=0.
1613	mp(:,:)=0.
1614	phi(:,:,:)=0.
1615	c
1616	c Ne pas affecter les valeurs entrees de u, v, h, et q
1617	c
1618	DO k = 1, klev
1619	DO i = 1, klon
1620	t_seri(i,k) = t(i,k)
1621	u_seri(i,k) = u(i,k)
1622	v_seri(i,k) = v(i,k)
1623	q_seri(i,k) = qx(i,k,ivap)
1624	ql_seri(i,k) = qx(i,k,iliq)
1625	qs_seri(i,k) = 0.
1626	ENDDO
1627	ENDDO
1628	IF (nqmax.GE.3) THEN
1629	DO iq = 3, nqmax
1630	DO k = 1, klev
1631	DO i = 1, klon
1632	tr_seri(i,k,iq-2) = qx(i,k,iq)
1633	ENDDO
1634	ENDDO
1635	ENDDO
1636	ELSE
1637	DO k = 1, klev
1638	DO i = 1, klon
1639	tr_seri(i,k,1) = 0.0
1640	ENDDO
1641	ENDDO
1642	ENDIF
1643	C
1644	DO i = 1, klon
1645	ztsol(i) = 0.
1646	ENDDO
1647	DO nsrf = 1, nbsrf
1648	DO i = 1, klon
1649	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1650	ENDDO
1651	ENDDO
1652	cIM
1653	IF (ip_ebil_phy.ge.1) THEN
1654	ztit='after dynamic'
1655	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
1656	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1657	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1658	C Comme les tendances de la physique sont ajoute dans la dynamique,
1659	C on devrait avoir que la variation d'entalpie par la dynamique
1660	C est egale a la variation de la physique au pas de temps precedent.
1661	C Donc la somme de ces 2 variations devrait etre nulle.
1662	call diagphy(airephy,ztit,ip_ebil_phy
1663	e , zero_v, zero_v, zero_v, zero_v, zero_v
1664	e , zero_v, zero_v, zero_v, ztsol
1665	e , d_h_vcol+d_h_vcol_phy, d_qt, 0.
1666	s , fs_bound, fq_bound )
1667	END IF
1668
1669	c Diagnostiquer la tendance dynamique
1670	c
1671	IF (ancien_ok) THEN
1672	DO k = 1, klev
1673	DO i = 1, klon
1674	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
1675	d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
1676	ENDDO
1677	ENDDO
1678	ELSE
1679	DO k = 1, klev
1680	DO i = 1, klon
1681	d_t_dyn(i,k) = 0.0
1682	d_q_dyn(i,k) = 0.0
1683	ENDDO
1684	ENDDO
1685	ancien_ok = .TRUE.
1686	ENDIF
1687	c
1688	c Ajouter le geopotentiel du sol:
1689	c
1690	DO k = 1, klev
1691	DO i = 1, klon
1692	zphi(i,k) = pphi(i,k) + pphis(i)
1693	ENDDO
1694	ENDDO
1695	c
1696	c Verifier les temperatures
1697	c
1698	cIM BEG
1699	IF (check) THEN
1700	amn=MIN(ftsol(1,is_ter),1000.)
1701	amx=MAX(ftsol(1,is_ter),-1000.)
1702	DO i=2, klon
1703	amn=MIN(ftsol(i,is_ter),amn)
1704	amx=MAX(ftsol(i,is_ter),amx)
1705	ENDDO
1706	c
1707	PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
1708	ENDIF !(check) THEN
1709	cIM END
1710	c
1711	CALL hgardfou(t_seri,ftsol,'debutphy')
1712	c
1713	cIM BEG
1714	IF (check) THEN
1715	amn=MIN(ftsol(1,is_ter),1000.)
1716	amx=MAX(ftsol(1,is_ter),-1000.)
1717	DO i=2, klon
1718	amn=MIN(ftsol(i,is_ter),amn)
1719	amx=MAX(ftsol(i,is_ter),amx)
1720	ENDDO
1721	c
1722	PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
1723	ENDIF !(check) THEN
1724	cIM END
1725	c
1726	c Incrementer le compteur de la physique
1727	c
1728	itap = itap + 1
1729	julien = MOD(NINT(xjour),360)
1730	if (julien .eq. 0) julien = 360
1731	c
1732	c Mettre en action les conditions aux limites (albedo, sst, etc.).
1733	c Prescrire l'ozone et calculer l'albedo sur l'ocean.
1734	c
1735	IF (MOD(itap-1,lmt_pas) .EQ. 0) THEN
1736	WRITE(lunout,*)' PHYS cond julien ',julien
1737	CALL ozonecm( FLOAT(julien), rlat, paprs, wo)
1738	ENDIF
1739	c
1740	c Re-evaporer l'eau liquide nuageuse
1741	c
1742	DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse
1743	DO i = 1, klon
1744	zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1745	c zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1746	zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1747	zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
1748	zb = MAX(0.0,ql_seri(i,k))
1749	za = - MAX(0.0,ql_seri(i,k))
1750	. * (zlvdcp(1.-zdelta)+zlsdcpzdelta)
1751	t_seri(i,k) = t_seri(i,k) + za
1752	q_seri(i,k) = q_seri(i,k) + zb
1753	ql_seri(i,k) = 0.0
1754	d_t_eva(i,k) = za
1755	d_q_eva(i,k) = zb
1756	ENDDO
1757	ENDDO
1758	cIM
1759	IF (ip_ebil_phy.ge.2) THEN
1760	ztit='after reevap'
1761	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,1,dtime
1762	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1763	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1764	call diagphy(airephy,ztit,ip_ebil_phy
1765	e , zero_v, zero_v, zero_v, zero_v, zero_v
1766	e , zero_v, zero_v, zero_v, ztsol
1767	e , d_h_vcol, d_qt, d_ec
1768	s , fs_bound, fq_bound )
1769	C
1770	END IF
1771
1772	c
1773	c=========================================================================
1774	! Calculs de l'orbite.
1775	! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
1776	! doit donc etre placé avant radlwsw et pbl_surface
1777
1778	! choix entre calcul de la longitude solaire vraie ou valeur fixee a
1779	! solarlong0
1780
1781	if (solarlong0<-999.) then
1782	CALL orbite(FLOAT(julien),zlongi,dist)
1783	else
1784	zlongi=solarlong0 ! longitude solaire vraie
1785	dist=1. ! distance au soleil / moyenne
1786	endif
1787
1788	if(prt_level.ge.1) print*,'Longitude solaire ',zlongi,solarlong0
1789
1790	! Avec ou sans cycle diurne
1791	IF (cycle_diurne) THEN
1792	zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
1793	CALL zenang(zlongi,gmtime,zdtime,rlat,rlon,rmu0,fract)
1794	ELSE
1795	rmu0 = -999.999
1796	ENDIF
1797
1798	if (mydebug) then
1799	call writefield_phy('u_seri',u_seri,llm)
1800	call writefield_phy('v_seri',v_seri,llm)
1801	call writefield_phy('t_seri',t_seri,llm)
1802	call writefield_phy('q_seri',q_seri,llm)
1803	endif
1804
1805	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1806	c Appel au pbl_surface : Planetary Boudary Layer et Surface
1807	c Cela implique tous les interactions des sous-surfaces et la partie diffusion
1808	c turbulent du couche limit.
1809	c
1810	c Certains varibales de sorties de pbl_surface sont utiliser que pour
1811	c ecriture des fihiers hist_XXXX.nc, ces sont :
1812	c qsol, zq2m, s_pblh, s_lcl,
1813	c s_capCL, s_oliqCL, s_cteiCL,s_pblT,
1814	c s_therm, s_trmb1, s_trmb2, s_trmb3,
1815	c zxrugs, zu10m, zv10m, fder,
1816	c zxqsurf, rh2m, zxfluxu, zxfluxv,
1817	c frugs, agesno, fsollw, fsolsw,
1818	c d_ts, fevap, fluxlat, t2m,
1819	c wfbils, wfbilo, fluxt, fluxu, fluxv,
1820	c
1821	c Certains ne sont pas utiliser du tout :
1822	c dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
1823	c
1824
1825	CALL pbl_surface(
1826	e dtime, date0, itap, julien,
1827	e debut, lafin,
1828	e rlon, rlat, rugoro, rmu0,
1829	e rain_fall, snow_fall, solsw, sollw,
1830	e t_seri, q_seri, u_seri, v_seri,
1831	e pplay, paprs, pctsrf,
1832	+ ftsol, falb1, falb2, u10m, v10m,
1833	s sollwdown, cdragh, cdragm, yu1, yv1,
1834	s albsol1, albsol2, sens, evap,
1835	s zxtsol, zxfluxlat, zt2m, qsat2m,
1836	s d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf,
1837	s ycoefh, pctsrf_new,
1838	d qsol, zq2m, s_pblh, s_lcl,
1839	d s_capCL, s_oliqCL, s_cteiCL,s_pblT,
1840	d s_therm, s_trmb1, s_trmb2, s_trmb3,
1841	d zxrugs, zu10m, zv10m, fder,
1842	d zxqsurf, rh2m, zxfluxu, zxfluxv,
1843	d frugs, agesno, fsollw, fsolsw,
1844	d d_ts, fevap, fluxlat, t2m,
1845	d wfbils, wfbilo, fluxt, fluxu, fluxv,
1846	- dsens, devap, zxsnow,
1847	- zxfluxt, zxfluxq, q2m, fluxq, pbl_tke )
1848	c
1849	c
1850	pctsrf(:,:) = pctsrf_new(:,:)
1851
1852	!-----------------------------------------------------------------------------------------
1853	! ajout des tendances de la diffusion turbulente
1854	CALL add_phys_tend(d_u_vdf,d_v_vdf,d_t_vdf,d_q_vdf,dql0,'vdf')
1855	!-----------------------------------------------------------------------------------------
1856
1857	if (mydebug) then
1858	call writefield_phy('u_seri',u_seri,llm)
1859	call writefield_phy('v_seri',v_seri,llm)
1860	call writefield_phy('t_seri',t_seri,llm)
1861	call writefield_phy('q_seri',q_seri,llm)
1862	endif
1863
1864
1865	IF (ip_ebil_phy.ge.2) THEN
1866	ztit='after surface_main'
1867	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
1868	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1869	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1870	call diagphy(airephy,ztit,ip_ebil_phy
1871	e , zero_v, zero_v, zero_v, zero_v, sens
1872	e , evap , zero_v, zero_v, ztsol
1873	e , d_h_vcol, d_qt, d_ec
1874	s , fs_bound, fq_bound )
1875	END IF
1876
1877	c =================================================================== c
1878	c Calcul de Qsat
1879
1880	DO k = 1, klev
1881	DO i = 1, klon
1882	zx_t = t_seri(i,k)
1883	IF (thermcep) THEN
1884	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
1885	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
1886	zx_qs = MIN(0.5,zx_qs)
1887	zcor = 1./(1.-retv*zx_qs)
1888	zx_qs = zx_qs*zcor
1889	ELSE
1890	IF (zx_t.LT.t_coup) THEN
1891	zx_qs = qsats(zx_t)/pplay(i,k)
1892	ELSE
1893	zx_qs = qsatl(zx_t)/pplay(i,k)
1894	ENDIF
1895	ENDIF
1896	zqsat(i,k)=zx_qs
1897	ENDDO
1898	ENDDO
1899
1900	if (prt_level.ge.1) then
1901	write(lunout,*) 'L qsat (g/kg) avant clouds_gno'
1902	write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
1903	endif
1904	c
1905	c Appeler la convection (au choix)
1906	c
1907	DO k = 1, klev
1908	DO i = 1, klon
1909	conv_q(i,k) = d_q_dyn(i,k)
1910	. + d_q_vdf(i,k)/dtime
1911	conv_t(i,k) = d_t_dyn(i,k)
1912	. + d_t_vdf(i,k)/dtime
1913	ENDDO
1914	ENDDO
1915	IF (check) THEN
1916	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
1917	WRITE(lunout,*) "avantcon=", za
1918	ENDIF
1919	zx_ajustq = .FALSE.
1920	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
1921	IF (zx_ajustq) THEN
1922	DO i = 1, klon
1923	z_avant(i) = 0.0
1924	ENDDO
1925	DO k = 1, klev
1926	DO i = 1, klon
1927	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k))
1928	. *(paprs(i,k)-paprs(i,k+1))/RG
1929	ENDDO
1930	ENDDO
1931	ENDIF
1932	IF (iflag_con.EQ.1) THEN
1933	stop'reactiver le call conlmd dans physiq.F'
1934	c CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
1935	c . d_t_con, d_q_con,
1936	c . rain_con, snow_con, ibas_con, itop_con)
1937	ELSE IF (iflag_con.EQ.2) THEN
1938	CALL conflx(dtime, paprs, pplay, t_seri, q_seri,
1939	e conv_t, conv_q, -evap, omega,
1940	s d_t_con, d_q_con, rain_con, snow_con,
1941	s pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
1942	s kcbot, kctop, kdtop, pmflxr, pmflxs)
1943	WHERE (rain_con < 0.) rain_con = 0.
1944	WHERE (snow_con < 0.) snow_con = 0.
1945	DO i = 1, klon
1946	ibas_con(i) = klev+1 - kcbot(i)
1947	itop_con(i) = klev+1 - kctop(i)
1948	ENDDO
1949	ELSE IF (iflag_con.GE.3) THEN
1950	c nb of tracers for the KE convection:
1951	c MAF la partie traceurs est faite dans phytrac
1952	c on met ntra=1 pour limiter les appels mais on peut
1953	c supprimer les calculs / ftra.
1954	ntra = 1
1955
1956	c=====================================================================================
1957	cajout pour la parametrisation des poches froides:
1958	ccalcul de t_wake et t_undi: si pas de poches froides, t_wake=t_undi=t_seri
1959	do k=1,klev
1960	do i=1,klon
1961	if (iflag_wake.eq.1) then
1962	t_wake(i,k) = t_seri(i,k)
1963	. +(1-wake_s(i))*wake_deltat(i,k)
1964	q_wake(i,k) = q_seri(i,k)
1965	. +(1-wake_s(i))*wake_deltaq(i,k)
1966	t_undi(i,k) = t_seri(i,k)
1967	. -wake_s(i)*wake_deltat(i,k)
1968	q_undi(i,k) = q_seri(i,k)
1969	. -wake_s(i)*wake_deltaq(i,k)
1970	else
1971	t_wake(i,k) = t_seri(i,k)
1972	q_wake(i,k) = q_seri(i,k)
1973	t_undi(i,k) = t_seri(i,k)
1974	q_undi(i,k) = q_seri(i,k)
1975	endif
1976	enddo
1977	enddo
1978
1979	cc-- Calcul de l'energie disponible ALE (J/kg) et de la puissance disponible ALP (W/m2)
1980	cc-- pour le soulevement des particules dans le modele convectif
1981	c
1982	do i = 1,klon
1983	ALE(i) = 0.
1984	ALP(i) = 0.
1985	enddo
1986	c
1987	ccalcul de ale_wake et alp_wake
1988	do i = 1,klon
1989	if (iflag_wake.eq.1) then
1990	ale_wake(i) = 0.5wake_cstar(i)*2
1991	alp_wake(i) = wake_fip(i)
1992	else
1993	ale_wake(i) = 0.
1994	alp_wake(i) = 0.
1995	endif
1996	enddo
1997	ccombinaison avec ale et alp de couche limite: constantes si pas de couplage, valeurs calculees
1998	cdans le thermique sinon
1999	if (iflag_coupl.eq.0) then
2000	if (debut) print*,'ALE et ALP imposes'
2001	do i = 1,klon
2002	con ne couple que ale
2003	c ALE(i) = max(ale_wake(i),Ale_bl(i))
2004	ALE(i) = max(ale_wake(i),ale_bl_prescr)
2005	con ne couple que alp
2006	c ALP(i) = alp_wake(i) + Alp_bl(i)
2007	ALP(i) = alp_wake(i) + alp_bl_prescr
2008	enddo
2009	else
2010	print*,'ALE et ALP couples au thermique'
2011	do i = 1,klon
2012	ALE(i) = max(ale_wake(i),Ale_bl(i))
2013	ALP(i) = alp_wake(i) + Alp_bl(i)
2014	c write(20,*)'ALE',ALE(i),Ale_bl(i),ale_wake(i)
2015	c write(21,*)'ALP',ALP(i),Alp_bl(i),alp_wake(i)
2016	enddo
2017	endif
2018
2019	cfin calcul ale et alp
2020	c=================================================================================================
2021
2022
2023	c sb, oct02:
2024	c Schema de convection modularise et vectorise:
2025	c (driver commun aux versions 3 et 4)
2026	c
2027	IF (ok_cvl) THEN ! new driver for convectL
2028
2029	CALL concvl (iflag_con,iflag_clos,
2030	. dtime,paprs,pplay,t_undi,q_undi,
2031	. t_wake,q_wake,
2032	. u_seri,v_seri,tr_seri,nbtr,
2033	. ALE,ALP,
2034	. ema_work1,ema_work2,
2035	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2036	. rain_con, snow_con, ibas_con, itop_con, sigd,
2037	. upwd,dnwd,dnwd0,
2038	. Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl,
2039	. pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd,
2040	. pmflxr,pmflxs,da,phi,mp,
2041	. ftd,fqd,lalim_conv,wght_th)
2042
2043	cIM cf. FH
2044	clwcon0=qcondc
2045	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
2046
2047	ELSE ! ok_cvl
2048	c MAF conema3 ne contient pas les traceurs
2049	CALL conema3 (dtime,
2050	. paprs,pplay,t_seri,q_seri,
2051	. u_seri,v_seri,tr_seri,ntra,
2052	. ema_work1,ema_work2,
2053	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2054	. rain_con, snow_con, ibas_con, itop_con,
2055	. upwd,dnwd,dnwd0,bas,top,
2056	. Ma,cape,tvp,rflag,
2057	. pbase
2058	. ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr
2059	. ,clwcon0)
2060
2061	ENDIF ! ok_cvl
2062
2063	c
2064	c Correction precip
2065	rain_con = rain_con * cvl_corr
2066	snow_con = snow_con * cvl_corr
2067	c
2068
2069	IF (.NOT. ok_gust) THEN
2070	do i = 1, klon
2071	wd(i)=0.0
2072	enddo
2073	ENDIF
2074
2075	c =================================================================== c
2076	c Calcul des proprietes des nuages convectifs
2077	c
2078
2079	c calcul des proprietes des nuages convectifs
2080	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
2081	call clouds_gno
2082	s (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
2083
2084	c =================================================================== c
2085
2086	DO i = 1, klon
2087	ema_pcb(i) = pbase(i)
2088	ENDDO
2089	DO i = 1, klon
2090
2091	! L'idicage de itop_con peut cacher un pb potentiel
2092	! FH sous la dictee de JYG, CR
2093	ema_pct(i) = paprs(i,itop_con(i)+1)
2094
2095	if (itop_con(i).gt.klev-3) then
2096	print*,'La convection monte trop haut '
2097	print*,'itop_con(,',i,',)=',itop_con(i)
2098	endif
2099	ENDDO
2100	DO i = 1, klon
2101	ema_cbmf(i) = ema_workcbmf(i)
2102	ENDDO
2103	ELSE IF (iflag_con.eq.0) THEN
2104	write(lunout,*) 'On n appelle pas la convection'
2105	clwcon0=0.
2106	rnebcon0=0.
2107	d_t_con=0.
2108	d_q_con=0.
2109	d_u_con=0.
2110	d_v_con=0.
2111	rain_con=0.
2112	snow_con=0.
2113	bas=1
2114	top=1
2115	ELSE
2116	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
2117	CALL abort
2118	ENDIF
2119
2120	c CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
2121	c . d_u_con, d_v_con)
2122
2123	!-----------------------------------------------------------------------------------------
2124	! ajout des tendances de la diffusion turbulente
2125	CALL add_phys_tend(d_u_con,d_v_con,d_t_con,d_q_con,dql0,'con')
2126	!-----------------------------------------------------------------------------------------
2127
2128	if (mydebug) then
2129	call writefield_phy('u_seri',u_seri,llm)
2130	call writefield_phy('v_seri',v_seri,llm)
2131	call writefield_phy('t_seri',t_seri,llm)
2132	call writefield_phy('q_seri',q_seri,llm)
2133	endif
2134
2135	cIM
2136	IF (ip_ebil_phy.ge.2) THEN
2137	ztit='after convect'
2138	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2139	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2140	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2141	call diagphy(airephy,ztit,ip_ebil_phy
2142	e , zero_v, zero_v, zero_v, zero_v, zero_v
2143	e , zero_v, rain_con, snow_con, ztsol
2144	e , d_h_vcol, d_qt, d_ec
2145	s , fs_bound, fq_bound )
2146	END IF
2147	C
2148	IF (check) THEN
2149	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2150	WRITE(lunout,*)"aprescon=", za
2151	zx_t = 0.0
2152	za = 0.0
2153	DO i = 1, klon
2154	za = za + airephy(i)/FLOAT(klon)
2155	zx_t = zx_t + (rain_con(i)+
2156	. snow_con(i))*airephy(i)/FLOAT(klon)
2157	ENDDO
2158	zx_t = zx_t/za*dtime
2159	WRITE(lunout,*)"Precip=", zx_t
2160	ENDIF
2161	IF (zx_ajustq) THEN
2162	DO i = 1, klon
2163	z_apres(i) = 0.0
2164	ENDDO
2165	DO k = 1, klev
2166	DO i = 1, klon
2167	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k))
2168	. *(paprs(i,k)-paprs(i,k+1))/RG
2169	ENDDO
2170	ENDDO
2171	DO i = 1, klon
2172	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime)
2173	. /z_apres(i)
2174	ENDDO
2175	DO k = 1, klev
2176	DO i = 1, klon
2177	IF (z_factor(i).GT.(1.0+1.0E-08) .OR.
2178	. z_factor(i).LT.(1.0-1.0E-08)) THEN
2179	q_seri(i,k) = q_seri(i,k) * z_factor(i)
2180	ENDIF
2181	ENDDO
2182	ENDDO
2183	ENDIF
2184	zx_ajustq=.FALSE.
2185
2186	c
2187	c=============================================================================
2188	cRR:Evolution de la poche froide: on ne fait pas de separation wake/env
2189	cpour la couche limite diffuse pour l instant
2190	c
2191	if (iflag_wake.eq.1) then
2192	DO k=1,klev
2193	DO i=1,klon
2194	dt_dwn(i,k) = ftd(i,k)
2195	wdt_PBL(i,k) = 0.
2196	dq_dwn(i,k) = fqd(i,k)
2197	wdq_PBL(i,k) = 0.
2198	M_dwn(i,k) = dnwd0(i,k)
2199	M_up(i,k) = upwd(i,k)
2200	dt_a(i,k) = d_t_con(i,k)/dtime - ftd(i,k)
2201	udt_PBL(i,k) = 0.
2202	dq_a(i,k) = d_q_con(i,k)/dtime - fqd(i,k)
2203	udq_PBL(i,k) = 0.
2204	ENDDO
2205	ENDDO
2206	c
2207	ccalcul caracteristiques de la poche froide
2208	call calWAKE (paprs,pplay,dtime
2209	: ,t_seri,q_seri,omega
2210	: ,dt_dwn,dq_dwn,M_dwn,M_up
2211	: ,dt_a,dq_a,sigd
2212	: ,wdt_PBL,wdq_PBL
2213	: ,udt_PBL,udq_PBL
2214	o ,wake_deltat,wake_deltaq,wake_dth
2215	o ,wake_h,wake_s,wake_dens
2216	o ,wake_pe,wake_fip,wake_gfl
2217	o ,dt_wake,dq_wake
2218	o ,wake_k, t_undi,q_undi
2219	o ,wake_omgbdth,wake_dp_omgb
2220	o ,wake_dtKE,wake_dqKE
2221	o ,wake_dtPBL,wake_dqPBL
2222	o ,wake_omg,wake_dp_deltomg
2223	o ,wake_spread,wake_Cstar,wake_d_deltat_gw
2224	o ,wake_ddeltat,wake_ddeltaq)
2225	c
2226	!-----------------------------------------------------------------------------------------
2227	! ajout des tendances des poches froides
2228	! Faire rapidement disparaitre l'ancien dt_wake pour garder un d_t_wake
2229	! coherent avec les autres d_t_...
2230	d_t_wake(:,:)=dt_wake(:,:)*dtime
2231	d_q_wake(:,:)=dq_wake(:,:)*dtime
2232	CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,'wake')
2233	!-----------------------------------------------------------------------------------------
2234
2235	endif
2236	c print*,'apres callwake iflag_cldcon=', iflag_cldcon
2237	c
2238	c===================================================================
2239	c Convection seche (thermiques ou ajustement)
2240	c===================================================================
2241	c
2242	call stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri
2243	s ,seuil_inversion,weak_inversion,dthmin)
2244
2245
2246
2247	d_t_ajsb(:,:)=0.
2248	d_q_ajsb(:,:)=0.
2249	d_t_ajs(:,:)=0.
2250	d_u_ajs(:,:)=0.
2251	d_v_ajs(:,:)=0.
2252	d_q_ajs(:,:)=0.
2253	clwcon0th(:,:)=0.
2254	c
2255	IF(prt_level>9)WRITE(lunout,*)
2256	. 'AVANT LA CONVECTION SECHE , iflag_thermals='
2257	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2258	if(iflag_thermals.lt.0) then
2259	c Rien
2260	c ====
2261	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
2262
2263
2264	else
2265
2266	c Thermiques
2267	c ==========
2268	IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals='
2269	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2270
2271
2272	if (iflag_thermals.gt.1) then
2273	call calltherm(pdtphys
2274	s ,pplay,paprs,pphi,weak_inversion
2275	s ,u_seri,v_seri,t_seri,q_seri,zqsat,debut
2276	s ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
2277	s ,fm_therm,entr_therm,zqasc,clwcon0th,lmax_th,ratqscth
2278	s ,ratqsdiff,zqsatth
2279	con rajoute ale et alp, et les caracteristiques de la couche alim
2280	s ,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0)
2281	endif
2282
2283	! call calltherm(pdtphys
2284	! s ,pplay,paprs,pphi
2285	! s ,u_seri,v_seri,t_seri,q_seri
2286	! s ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
2287	! s ,fm_therm,entr_therm)
2288
2289	c Ajustement sec
2290	c ==============
2291
2292	! Dans le cas où on active les thermiques, on fait partir l'ajustement
2293	! a partir du sommet des thermiques.
2294	! Dans le cas contraire, on demarre au niveau 1.
2295
2296	if (iflag_thermals.ge.13.or.iflag_thermals.eq.0) then
2297
2298	if(iflag_thermals.eq.0) then
2299	IF(prt_level>9)WRITE(lunout,*)'ajsec'
2300	limbas(:)=1
2301	else
2302	limbas(:)=lmax_th(:)
2303	endif
2304
2305	! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
2306	! pour des test de convergence numerique.
2307	! Le nouveau ajsec est a priori mieux, meme pour le cas
2308	! iflag_thermals = 0 (l'ancienne version peut faire des tendances
2309	! non nulles numeriquement pour des mailles non concernees.
2310
2311	if (iflag_thermals.eq.0) then
2312	CALL ajsec_convV2(paprs, pplay, t_seri,q_seri
2313	s , d_t_ajsb, d_q_ajsb)
2314	else
2315	CALL ajsec(paprs, pplay, t_seri,q_seri,limbas
2316	s , d_t_ajsb, d_q_ajsb)
2317	endif
2318
2319	!-----------------------------------------------------------------------------------------
2320	! ajout des tendances de l'ajustement sec ou des thermiques
2321	CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,'ajsb')
2322	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
2323	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
2324
2325	!-----------------------------------------------------------------------------------------
2326
2327	endif
2328
2329	endif
2330	c
2331	c===================================================================
2332	cIM
2333	IF (ip_ebil_phy.ge.2) THEN
2334	ztit='after dry_adjust'
2335	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2336	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2337	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2338	END IF
2339
2340
2341	c-------------------------------------------------------------------------
2342	c Caclul des ratqs
2343	c-------------------------------------------------------------------------
2344
2345	c print*,'calcul des ratqs'
2346	c ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
2347	c ----------------
2348	c on ecrase le tableau ratqsc calcule par clouds_gno
2349	if (iflag_cldcon.eq.1) then
2350	do k=1,klev
2351	do i=1,klon
2352	if(ptconv(i,k)) then
2353	ratqsc(i,k)=ratqsbas
2354	s +fact_cldcon*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
2355	else
2356	ratqsc(i,k)=0.
2357	endif
2358	enddo
2359	enddo
2360	else if (iflag_cldcon.eq.4) then
2361	ptconvth(:,:)=.false.
2362	ratqsc(:,:)=0.
2363	if(prt_level.ge.9) print*,'avant clouds_gno thermique'
2364	call clouds_gno
2365	s (klon,klev,q_seri,zqsat,clwcon0th,ptconvth,ratqsc,rnebcon0th)
2366	if(prt_level.ge.9) print*,' CLOUDS_GNO OK'
2367	endif
2368
2369	c ratqs stables
2370	c -------------
2371
2372	if (iflag_ratqs.eq.0) then
2373
2374	! Le cas iflag_ratqs=0 correspond a la version IPCC 2005 du modele.
2375	do k=1,klev
2376	do i=1, klon
2377	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2378	s min((paprs(i,1)-pplay(i,k))/(paprs(i,1)-30000.),1.)
2379	enddo
2380	enddo
2381
2382	! Pour iflag_ratqs=1 ou 2, le ratqs est constant au dessus de
2383	! 300 hPa (ratqshaut), varie lineariement en fonction de la pression
2384	! entre 600 et 300 hPa et est soit constant (ratqsbas) pour iflag_ratqs=1
2385	! soit lineaire (entre 0 a la surface et ratqsbas) pour iflag_ratqs=2
2386	! Il s'agit de differents tests dans la phase de reglage du modele
2387	! avec thermiques.
2388
2389	else if (iflag_ratqs.eq.1) then
2390
2391	do k=1,klev
2392	do i=1, klon
2393	if (pplay(i,k).ge.60000.) then
2394	ratqss(i,k)=ratqsbas
2395	else if ((pplay(i,k).ge.30000.).and.
2396	s (pplay(i,k).lt.60000.)) then
2397	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2398	s (60000.-pplay(i,k))/(60000.-30000.)
2399	else
2400	ratqss(i,k)=ratqshaut
2401	endif
2402	enddo
2403	enddo
2404
2405	else
2406
2407	do k=1,klev
2408	do i=1, klon
2409	if (pplay(i,k).ge.60000.) then
2410	ratqss(i,k)=ratqsbas
2411	s *(paprs(i,1)-pplay(i,k))/(paprs(i,1)-60000.)
2412	else if ((pplay(i,k).ge.30000.).and.
2413	s (pplay(i,k).lt.60000.)) then
2414	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2415	s (60000.-pplay(i,k))/(60000.-30000.)
2416	else
2417	ratqss(i,k)=ratqshaut
2418	endif
2419	enddo
2420	enddo
2421	endif
2422
2423
2424
2425
2426	c ratqs final
2427	c -----------
2428
2429	if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2
2430	s .or.iflag_cldcon.eq.4) then
2431
2432	! On ajoute une constante au ratqsc*2 pour tenir compte de
2433	! fluctuations turbulentes de petite echelle
2434
2435	do k=1,klev
2436	do i=1,klon
2437	if ((fm_therm(i,k).gt.1.e-10)) then
2438	ratqsc(i,k)=sqrt(ratqsc(i,k)2+0.052)
2439	endif
2440	enddo
2441	enddo
2442
2443	! les ratqs sont une conbinaison de ratqss et ratqsc
2444	! 1800s, en dur pour le moment, est le temps de
2445	! relaxation des ratqs
2446
2447	facteur=exp(-pdtphys/1800.)
2448
2449	print*,'WARNING ratqs a revoir '
2450	ratqs(:,:)=ratqsc(:,:)(1.-facteur)+ratqs(:,:)facteur
2451	ratqs(:,:)=sqrt(ratqs(:,:)2+ratqss(:,:)2)
2452
2453	else
2454	! on ne prend que le ratqs stable pour fisrtilp
2455	ratqs(:,:)=ratqss(:,:)
2456	endif
2457
2458
2459	c
2460	c Appeler le processus de condensation a grande echelle
2461	c et le processus de precipitation
2462	c-------------------------------------------------------------------------
2463	CALL fisrtilp(dtime,paprs,pplay,
2464	. t_seri, q_seri,ptconv,ratqs,
2465	. d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq,
2466	. rain_lsc, snow_lsc,
2467	. pfrac_impa, pfrac_nucl, pfrac_1nucl,
2468	. frac_impa, frac_nucl,
2469	. prfl, psfl, rhcl)
2470
2471	WHERE (rain_lsc < 0) rain_lsc = 0.
2472	WHERE (snow_lsc < 0) snow_lsc = 0.
2473	!-----------------------------------------------------------------------------------------
2474	! ajout des tendances de la diffusion turbulente
2475	CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,'lsc')
2476	!-----------------------------------------------------------------------------------------
2477	DO k = 1, klev
2478	DO i = 1, klon
2479	cldfra(i,k) = rneb(i,k)
2480	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
2481	ENDDO
2482	ENDDO
2483	IF (check) THEN
2484	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2485	WRITE(lunout,*)"apresilp=", za
2486	zx_t = 0.0
2487	za = 0.0
2488	DO i = 1, klon
2489	za = za + airephy(i)/FLOAT(klon)
2490	zx_t = zx_t + (rain_lsc(i)
2491	. + snow_lsc(i))*airephy(i)/FLOAT(klon)
2492	ENDDO
2493	zx_t = zx_t/za*dtime
2494	WRITE(lunout,*)"Precip=", zx_t
2495	ENDIF
2496	cIM
2497	IF (ip_ebil_phy.ge.2) THEN
2498	ztit='after fisrt'
2499	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2500	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2501	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2502	call diagphy(airephy,ztit,ip_ebil_phy
2503	e , zero_v, zero_v, zero_v, zero_v, zero_v
2504	e , zero_v, rain_lsc, snow_lsc, ztsol
2505	e , d_h_vcol, d_qt, d_ec
2506	s , fs_bound, fq_bound )
2507	END IF
2508
2509	if (mydebug) then
2510	call writefield_phy('u_seri',u_seri,llm)
2511	call writefield_phy('v_seri',v_seri,llm)
2512	call writefield_phy('t_seri',t_seri,llm)
2513	call writefield_phy('q_seri',q_seri,llm)
2514	endif
2515
2516	c
2517	c-------------------------------------------------------------------
2518	c PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
2519	c-------------------------------------------------------------------
2520
2521	c 1. NUAGES CONVECTIFS
2522	c
2523	cIM cf FH
2524	c IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
2525	IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
2526	snow_tiedtke=0.
2527	c print*,'avant calcul de la pseudo precip '
2528	c print*,'iflag_cldcon',iflag_cldcon
2529	if (iflag_cldcon.eq.-1) then
2530	rain_tiedtke=rain_con
2531	else
2532	c print*,'calcul de la pseudo precip '
2533	rain_tiedtke=0.
2534	c print*,'calcul de la pseudo precip 0'
2535	do k=1,klev
2536	do i=1,klon
2537	if (d_q_con(i,k).lt.0.) then
2538	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys
2539	s *(paprs(i,k)-paprs(i,k+1))/rg
2540	endif
2541	enddo
2542	enddo
2543	endif
2544	c
2545	c call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
2546	c
2547
2548	c Nuages diagnostiques pour Tiedtke
2549	CALL diagcld1(paprs,pplay,
2550	cIM cf FH . rain_con,snow_con,ibas_con,itop_con,
2551	. rain_tiedtke,snow_tiedtke,ibas_con,itop_con,
2552	. diafra,dialiq)
2553	DO k = 1, klev
2554	DO i = 1, klon
2555	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2556	cldliq(i,k) = dialiq(i,k)
2557	cldfra(i,k) = diafra(i,k)
2558	ENDIF
2559	ENDDO
2560	ENDDO
2561
2562	ELSE IF (iflag_cldcon.ge.3) THEN
2563	c On prend pour les nuages convectifs le max du calcul de la
2564	c convection et du calcul du pas de temps precedent diminue d'un facteur
2565	c facttemps
2566	facteur = pdtphys *facttemps
2567	do k=1,klev
2568	do i=1,klon
2569	rnebcon(i,k)=rnebcon(i,k)*facteur
2570	if (rnebcon0(i,k)clwcon0(i,k).gt.rnebcon(i,k)clwcon(i,k))
2571	s then
2572	rnebcon(i,k)=rnebcon0(i,k)
2573	clwcon(i,k)=clwcon0(i,k)
2574	endif
2575	enddo
2576	enddo
2577
2578	c
2579	cjq - introduce the aerosol direct and first indirect radiative forcings
2580	cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
2581	IF (ok_ade.OR.ok_aie) THEN
2582	#if defined(CPP_COUPLE) && !defined(INCA_AER)
2583	! Get sulfate aerosol distribution
2584	CALL readsulfate(rjourvrai, debut, sulfate)
2585	CALL readsulfate_preind(rjourvrai, debut, sulfate_pi)
2586
2587	! Calculate aerosol optical properties (Olivier Boucher)
2588	CALL aeropt(pplay, paprs, t_seri, sulfate, rhcl,
2589	. tau_ae, piz_ae, cg_ae, aerindex)
2590	#endif
2591	#if !defined(CPP_COUPLE)
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	cym
2601	ELSE
2602	tau_ae(:,:,:)=0.0
2603	piz_ae(:,:,:)=0.0
2604	cg_ae(:,:,:)=0.0
2605	cym
2606	ENDIF
2607
2608	c
2609	cIM calcul nuages par le simulateur ISCCP
2610	c
2611	#ifdef histISCCP
2612	IF (ok_isccp) THEN
2613	cIM appel simulateur toutes les NINT(freq_ISCCP/dtime) heures
2614	IF (MOD(itap,NINT(freq_ISCCP/dtime)).EQ.0) THEN
2615	#include "calcul_simulISCCP.h"
2616	ENDIF !(MOD(itap,NINT(freq_ISCCP/dtime))
2617	ENDIF !ok_isccp
2618	#endif
2619
2620	c On prend la somme des fractions nuageuses et des contenus en eau
2621	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
2622	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
2623
2624	ENDIF
2625
2626	c
2627	c 2. NUAGES STARTIFORMES
2628	c
2629	IF (ok_stratus) THEN
2630	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
2631	DO k = 1, klev
2632	DO i = 1, klon
2633	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2634	cldliq(i,k) = dialiq(i,k)
2635	cldfra(i,k) = diafra(i,k)
2636	ENDIF
2637	ENDDO
2638	ENDDO
2639	ENDIF
2640	c
2641	c Precipitation totale
2642	c
2643	DO i = 1, klon
2644	rain_fall(i) = rain_con(i) + rain_lsc(i)
2645	snow_fall(i) = snow_con(i) + snow_lsc(i)
2646	ENDDO
2647	cIM
2648	IF (ip_ebil_phy.ge.2) THEN
2649	ztit="after diagcld"
2650	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2651	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2652	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2653	END IF
2654	c
2655	c Calculer l'humidite relative pour diagnostique
2656	c
2657	DO k = 1, klev
2658	DO i = 1, klon
2659	zx_t = t_seri(i,k)
2660	IF (thermcep) THEN
2661	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2662	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2663	zx_qs = MIN(0.5,zx_qs)
2664	zcor = 1./(1.-retv*zx_qs)
2665	zx_qs = zx_qs*zcor
2666	ELSE
2667	IF (zx_t.LT.t_coup) THEN
2668	zx_qs = qsats(zx_t)/pplay(i,k)
2669	ELSE
2670	zx_qs = qsatl(zx_t)/pplay(i,k)
2671	ENDIF
2672	ENDIF
2673	zx_rh(i,k) = q_seri(i,k)/zx_qs
2674	zqsat(i,k)=zx_qs
2675	ENDDO
2676	ENDDO
2677
2678	cIM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
2679	c equivalente a 2m (tpote) pour diagnostique
2680	c
2681	DO i = 1, klon
2682	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
2683	IF (thermcep) THEN
2684	IF(zt2m(i).LT.RTT) then
2685	Lheat=RLSTT
2686	ELSE
2687	Lheat=RLVTT
2688	ENDIF
2689	ELSE
2690	IF (zt2m(i).LT.RTT) THEN
2691	Lheat=RLSTT
2692	ELSE
2693	Lheat=RLVTT
2694	ENDIF
2695	ENDIF
2696	tpote(i) = tpot(i)*
2697	. EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
2698	ENDDO
2699
2700
2701	#ifdef INCA
2702	call VTe(VTphysiq)
2703	call VTb(VTinca)
2704	calday = FLOAT(julien) + gmtime
2705
2706	#ifdef INCA_AER
2707	call AEROSOL_METEO_CALC(calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs
2708	& ,prfl,psfl,pctsrf,airephy,xjour,rlat,rlon,u10m,v10m)
2709	#endif
2710
2711	#ifdef INCAINFO
2712	WRITE(lunout,*)'Appel CHEMHOOK_BEGIN ...'
2713	#endif
2714
2715	zxsnow_dummy(:) = 0.0
2716
2717	CALL chemhook_begin (calday,
2718	#if defined(INCA) && !defined(INCA_CH4) && !defined(INCA_NMHC) && !defined(INCA_AER)
2719	$ julien,
2720	$ gmtime,
2721	#endif
2722	$ pctsrf(1,1),
2723	$ rlat,
2724	$ rlon,
2725	$ airephy,
2726	$ paprs,
2727	$ pplay,
2728	$ ycoefh,
2729	$ pphi,
2730	$ t_seri,
2731	$ u,
2732	$ v,
2733	$ wo,
2734	$ q_seri,
2735	$ zxtsol,
2736	$ zxsnow_dummy,
2737	$ solsw,
2738	$ albsol1,
2739	$ rain_fall,
2740	$ snow_fall,
2741	$ itop_con,
2742	$ ibas_con,
2743	$ cldfra,
2744	$ iim,
2745	$ jjm,
2746	#ifdef INCA_AER
2747	$ tr_seri,
2748	$ ftsol,
2749	$ paprs,
2750	$ cdragh,
2751	$ cdragm,
2752	$ pctsrf,
2753	$ pdtphys,
2754	$ itap)
2755	#else
2756	$ tr_seri)
2757	#endif
2758
2759
2760	#ifdef INCAINFO
2761	WRITE(lunout,*)'OK.'
2762	#endif
2763	call VTe(VTinca)
2764	call VTb(VTphysiq)
2765	#endif
2766	c
2767	c Calculer les parametres optiques des nuages et quelques
2768	c parametres pour diagnostiques:
2769	c
2770	if (ok_newmicro) then
2771	CALL newmicro (paprs, pplay,ok_newmicro,
2772	. t_seri, cldliq, cldfra, cldtau, cldemi,
2773	. cldh, cldl, cldm, cldt, cldq,
2774	. flwp, fiwp, flwc, fiwc,
2775	e ok_aie,
2776	#if defined(CPP_COUPLE) && defined(INCA_AER)
2777	e ccm(:,:,1), ccm(:,:,2),
2778	#else
2779	e sulfate, sulfate_pi,
2780	#endif
2781	e bl95_b0, bl95_b1,
2782	s cldtaupi, re, fl)
2783	else
2784	CALL nuage (paprs, pplay,
2785	. t_seri, cldliq, cldfra, cldtau, cldemi,
2786	. cldh, cldl, cldm, cldt, cldq,
2787	e ok_aie,
2788	e sulfate, sulfate_pi,
2789	e bl95_b0, bl95_b1,
2790	s cldtaupi, re, fl)
2791
2792	endif
2793	c
2794	c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
2795	c
2796	IF (MOD(itaprad,radpas).EQ.0) THEN
2797
2798	DO i = 1, klon
2799	albsol1(i) = falb1(i,is_oce) * pctsrf(i,is_oce)
2800	. + falb1(i,is_lic) * pctsrf(i,is_lic)
2801	. + falb1(i,is_ter) * pctsrf(i,is_ter)
2802	. + falb1(i,is_sic) * pctsrf(i,is_sic)
2803	albsol2(i) = falb2(i,is_oce) * pctsrf(i,is_oce)
2804	. + falb2(i,is_lic) * pctsrf(i,is_lic)
2805	. + falb2(i,is_ter) * pctsrf(i,is_ter)
2806	. + falb2(i,is_sic) * pctsrf(i,is_sic)
2807	ENDDO
2808
2809	if (mydebug) then
2810	call writefield_phy('u_seri',u_seri,llm)
2811	call writefield_phy('v_seri',v_seri,llm)
2812	call writefield_phy('t_seri',t_seri,llm)
2813	call writefield_phy('q_seri',q_seri,llm)
2814	endif
2815
2816
2817	#if defined(CPP_COUPLE) && defined(INCA_AER)
2818	CALL radlwsw_inca ! nouveau rayonnement (compatible Arpege-IFS)
2819	e (kdlon,kflev,dist, rmu0, fract,
2820	e paprs, pplay,zxtsol,albsol, albsollw, t_seri,q_seri,
2821	e wo,
2822	e cldfra, cldemi, cldtau,
2823	s heat,heat0,cool,cool0,radsol,albpla,
2824	s topsw,toplw,solsw,sollw,
2825	s sollwdown,
2826	s topsw0,toplw0,solsw0,sollw0,
2827	s lwdn0, lwdn, lwup0, lwup,
2828	s swdn0, swdn, swup0, swup,
2829	e ok_ade, ok_aie, ! new for aerosol radiative effects
2830	e tau_inca, piz_inca, cg_inca, ! ="=
2831	s topswad_inca, solswad_inca, ! ="=
2832	s topswad0_inca, solswad0_inca, ! ="=
2833	s topsw_inca, topsw0_inca,
2834	s solsw_inca, solsw0_inca,
2835	e cldtaupi, ! ="=
2836	s topswai_inca, solswai_inca) ! ="=
2837	#else
2838	CALL radlwsw ! nouveau rayonnement (compatible Arpege-IFS)
2839	e (dist, rmu0, fract,
2840	e paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri,
2841	e wo,
2842	e cldfra, cldemi, cldtau,
2843	s heat,heat0,cool,cool0,radsol,albpla,
2844	s topsw,toplw,solsw,sollw,
2845	s sollwdown,
2846	s topsw0,toplw0,solsw0,sollw0,
2847	s lwdn0, lwdn, lwup0, lwup,
2848	s swdn0, swdn, swup0, swup,
2849	e ok_ade, ok_aie, ! new for aerosol radiative effects
2850	e tau_ae, piz_ae, cg_ae, ! ="=
2851	s topswad, solswad, ! ="=
2852	e cldtaupi, ! ="=
2853	s topswai, solswai) ! ="=
2854	#endif
2855	itaprad = 0
2856	ENDIF
2857	itaprad = itaprad + 1
2858
2859	if (iflag_radia.eq.0) then
2860	print *,'--------------------------------------------------'
2861	print *,'>>>> ATTENTION rayonnement desactive pour ce cas'
2862	print *,'>>>> heat et cool mis a zero '
2863	print *,'--------------------------------------------------'
2864	heat=0.
2865	cool=0.
2866	endif
2867
2868
2869	c
2870	c Ajouter la tendance des rayonnements (tous les pas)
2871	c
2872	DO k = 1, klev
2873	DO i = 1, klon
2874	t_seri(i,k) = t_seri(i,k)
2875	. + (heat(i,k)-cool(i,k)) * dtime/86400.
2876	ENDDO
2877	ENDDO
2878	c
2879	if (mydebug) then
2880	call writefield_phy('u_seri',u_seri,llm)
2881	call writefield_phy('v_seri',v_seri,llm)
2882	call writefield_phy('t_seri',t_seri,llm)
2883	call writefield_phy('q_seri',q_seri,llm)
2884	endif
2885
2886	cIM
2887	IF (ip_ebil_phy.ge.2) THEN
2888	ztit='after rad'
2889	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2890	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2891	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2892	call diagphy(airephy,ztit,ip_ebil_phy
2893	e , topsw, toplw, solsw, sollw, zero_v
2894	e , zero_v, zero_v, zero_v, ztsol
2895	e , d_h_vcol, d_qt, d_ec
2896	s , fs_bound, fq_bound )
2897	END IF
2898	c
2899	c
2900	c Calculer l'hydrologie de la surface
2901	c
2902	c CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
2903	c . agesno, ftsol,fqsurf,fsnow, ruis)
2904	c
2905
2906	c
2907	c Calculer le bilan du sol et la derive de temperature (couplage)
2908	c
2909	DO i = 1, klon
2910	c bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
2911	c a la demande de JLD
2912	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
2913	ENDDO
2914	c
2915	cmoddeblott(jan95)
2916	c Appeler le programme de parametrisation de l'orographie
2917	c a l'echelle sous-maille:
2918	c
2919	IF (ok_orodr) THEN
2920	c
2921	c selection des points pour lesquels le shema est actif:
2922	igwd=0
2923	DO i=1,klon
2924	itest(i)=0
2925	c IF ((zstd(i).gt.10.0)) THEN
2926	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
2927	itest(i)=1
2928	igwd=igwd+1
2929	idx(igwd)=i
2930	ENDIF
2931	ENDDO
2932	c igwdim=MAX(1,igwd)
2933	c
2934	CALL drag_noro(klon,klev,dtime,paprs,pplay,
2935	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
2936	e igwd,idx,itest,
2937	e t_seri, u_seri, v_seri,
2938	cIM 141004 s zulow, zvlow, zustr, zvstr,
2939	s zulow, zvlow, zustrdr, zvstrdr,
2940	s d_t_oro, d_u_oro, d_v_oro)
2941	c
2942	c ajout des tendances
2943	!-----------------------------------------------------------------------------------------
2944	! ajout des tendances de la trainee de l'orographie
2945	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,'oro')
2946	!-----------------------------------------------------------------------------------------
2947	c
2948	ENDIF ! fin de test sur ok_orodr
2949	c
2950	if (mydebug) then
2951	call writefield_phy('u_seri',u_seri,llm)
2952	call writefield_phy('v_seri',v_seri,llm)
2953	call writefield_phy('t_seri',t_seri,llm)
2954	call writefield_phy('q_seri',q_seri,llm)
2955	endif
2956
2957	IF (ok_orolf) THEN
2958	c
2959	c selection des points pour lesquels le shema est actif:
2960	igwd=0
2961	DO i=1,klon
2962	itest(i)=0
2963	IF ((zpic(i)-zmea(i)).GT.100.) THEN
2964	itest(i)=1
2965	igwd=igwd+1
2966	idx(igwd)=i
2967	ENDIF
2968	ENDDO
2969	c igwdim=MAX(1,igwd)
2970	c
2971	CALL lift_noro(klon,klev,dtime,paprs,pplay,
2972	e rlat,zmea,zstd,zpic,
2973	e itest,
2974	e t_seri, u_seri, v_seri,
2975	s zulow, zvlow, zustrli, zvstrli,
2976	s d_t_lif, d_u_lif, d_v_lif)
2977	c
2978	!-----------------------------------------------------------------------------------------
2979	! ajout des tendances de la portance de l'orographie
2980	CALL add_phys_tend(d_u_lif,d_v_lif,d_t_lif,dq0,dql0,'lif')
2981	!-----------------------------------------------------------------------------------------
2982	c
2983	ENDIF ! fin de test sur ok_orolf
2984	c
2985	cIM cf. FLott BEG
2986	C STRESS NECESSAIRES: TOUTE LA PHYSIQUE
2987
2988	if (mydebug) then
2989	call writefield_phy('u_seri',u_seri,llm)
2990	call writefield_phy('v_seri',v_seri,llm)
2991	call writefield_phy('t_seri',t_seri,llm)
2992	call writefield_phy('q_seri',q_seri,llm)
2993	endif
2994
2995	DO i = 1, klon
2996	zustrph(i)=0.
2997	zvstrph(i)=0.
2998	ENDDO
2999	DO k = 1, klev
3000	DO i = 1, klon
3001	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
3002	c (paprs(i,k)-paprs(i,k+1))/rg
3003	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
3004	c (paprs(i,k)-paprs(i,k+1))/rg
3005	ENDDO
3006	ENDDO
3007	c
3008	cIM calcul composantes axiales du moment angulaire et couple des montagnes
3009	c
3010	IF (is_sequential) THEN
3011
3012	CALL aaam_bud (27,klon,klev,rjourvrai,gmtime,
3013	C ra,rg,romega,
3014	C rlat,rlon,pphis,
3015	C zustrdr,zustrli,zustrph,
3016	C zvstrdr,zvstrli,zvstrph,
3017	C paprs,u,v,
3018	C aam, torsfc)
3019	ENDIF
3020	cIM cf. FLott END
3021	cIM
3022	IF (ip_ebil_phy.ge.2) THEN
3023	ztit='after orography'
3024	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
3025	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3026	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3027	END IF
3028	c
3029	c
3030	cAA
3031	cAA Installation de l'interface online-offline pour traceurs
3032	cAA
3033	c====================================================================
3034	c Calcul des tendances traceurs
3035	c====================================================================
3036	C
3037	call phytrac ( rnpb,
3038	I itap,
3039	I julien,
3040	I gmtime,
3041	I debut,
3042	I lafin,
3043	I nqmax-2,
3044	I nlon,
3045	I nlev,
3046	I dtime,
3047	I u,
3048	I v,
3049	I t,
3050	I paprs,
3051	I pplay,
3052	I pmfu,
3053	I pmfd,
3054	I pen_u,
3055	I pde_u,
3056	I pen_d,
3057	I pde_d,
3058	I ycoefh,
3059	I fm_therm,
3060	I entr_therm,
3061	I yu1,
3062	I yv1,
3063	I ftsol,
3064	I pctsrf,
3065	I rlat,
3066	I frac_impa,
3067	I frac_nucl,
3068	I rlon,
3069	I presnivs,
3070	I pphis,
3071	I pphi,
3072	I albsol1,
3073	I qx(1,1,1),
3074	I rhcl,
3075	I cldfra,
3076	I rneb,
3077	I diafra,
3078	I cldliq,
3079	I itop_con,
3080	I ibas_con,
3081	I pmflxr,
3082	I pmflxs,
3083	I prfl,
3084	I psfl,
3085	I da,
3086	I phi,
3087	I mp,
3088	I upwd,
3089	I dnwd,
3090	#ifdef INCA
3091	I flxmass_w,
3092	#if defined(INCA_AER) && defined(CPP_COUPLE)
3093	I tau_inca,
3094	I piz_inca,
3095	I cg_inca,
3096	I ccm,
3097	I rfname,
3098	#endif
3099	#endif
3100	O tr_seri)
3101
3102	IF (offline) THEN
3103
3104	print*,'Attention on met a 0 les thermiques pour phystoke'
3105	call phystokenc (
3106	I nlon,nlev,pdtphys,rlon,rlat,
3107	I t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
3108	I fm_therm,entr_therm,
3109	I ycoefh,yu1,yv1,ftsol,pctsrf,
3110	I frac_impa, frac_nucl,
3111	I pphis,airephy,dtime,itap)
3112
3113
3114	ENDIF
3115
3116	c
3117	c Calculer le transport de l'eau et de l'energie (diagnostique)
3118	c
3119	CALL transp (paprs,zxtsol,
3120	e t_seri, q_seri, u_seri, v_seri, zphi,
3121	s ve, vq, ue, uq)
3122	c
3123	cIM global posePB BEG
3124	IF(1.EQ.0) THEN
3125	c
3126	CALL transp_lay (paprs,zxtsol,
3127	e t_seri, q_seri, u_seri, v_seri, zphi,
3128	s ve_lay, vq_lay, ue_lay, uq_lay)
3129	c
3130	ENDIF !(1.EQ.0) THEN
3131	cIM global posePB END
3132	c Accumuler les variables a stocker dans les fichiers histoire:
3133	c
3134	c+jld ec_conser
3135	DO k = 1, klev
3136	DO i = 1, klon
3137	ZRCPD = RCPD(1.0+RVTMP2q_seri(i,k))
3138	d_t_ec(i,k)=0.5/ZRCPD
3139	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
3140	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
3141	d_t_ec(i,k) = d_t_ec(i,k)/dtime
3142	END DO
3143	END DO
3144	c-jld ec_conser
3145	cIM
3146	IF (ip_ebil_phy.ge.1) THEN
3147	ztit='after physic'
3148	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
3149	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3150	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3151	C Comme les tendances de la physique sont ajoute dans la dynamique,
3152	C on devrait avoir que la variation d'entalpie par la dynamique
3153	C est egale a la variation de la physique au pas de temps precedent.
3154	C Donc la somme de ces 2 variations devrait etre nulle.
3155	call diagphy(airephy,ztit,ip_ebil_phy
3156	e , topsw, toplw, solsw, sollw, sens
3157	e , evap, rain_fall, snow_fall, ztsol
3158	e , d_h_vcol, d_qt, d_ec
3159	s , fs_bound, fq_bound )
3160	C
3161	d_h_vcol_phy=d_h_vcol
3162	C
3163	END IF
3164	C
3165	c=======================================================================
3166	c SORTIES
3167	c=======================================================================
3168
3169	cIM Interpolation sur les niveaux de pression du NMC
3170	c -------------------------------------------------
3171	c
3172	#include "calcul_STDlev.h"
3173	c
3174	c slp sea level pressure
3175	slp(:) = paprs(:,1)exp(pphis(:)/(RDt_seri(:,1)))
3176	c
3177	ccc prw = eau precipitable
3178	DO i = 1, klon
3179	prw(i) = 0.
3180	DO k = 1, klev
3181	prw(i) = prw(i) +
3182	. q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
3183	ENDDO
3184	ENDDO
3185	c
3186	cIM initialisation + calculs divers diag AMIP2
3187	c
3188	#include "calcul_divers.h"
3189	c
3190	#ifdef INCA
3191	call VTe(VTphysiq)
3192	call VTb(VTinca)
3193	#ifdef INCAINFO
3194	WRITE(lunout,*)'Appel CHEMHOOK_END ...'
3195	#endif
3196	CALL chemhook_end (calday,
3197	$ dtime,
3198	$ pplay,
3199	$ t_seri,
3200	$ tr_seri,
3201	$ nbtr,
3202	$ paprs,
3203	$ q_seri,
3204	$ annee_ref,
3205	$ day_ini,
3206	$ airephy,
3207	#ifdef INCA_AER
3208	$ xjour,
3209	$ pphi,
3210	$ pphis,
3211	$ zx_rh)
3212	#else
3213	$ xjour)
3214	#endif
3215	#ifdef INCAINFO
3216	WRITE(lunout,*)'OK.'
3217	#endif
3218	call VTe(VTinca)
3219	call VTb(VTphysiq)
3220	#endif
3221
3222	c=============================================================
3223	c
3224	c Convertir les incrementations en tendances
3225	c
3226	if (mydebug) then
3227	call writefield_phy('u_seri',u_seri,llm)
3228	call writefield_phy('v_seri',v_seri,llm)
3229	call writefield_phy('t_seri',t_seri,llm)
3230	call writefield_phy('q_seri',q_seri,llm)
3231	endif
3232
3233	DO k = 1, klev
3234	DO i = 1, klon
3235	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
3236	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
3237	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
3238	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
3239	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
3240	ENDDO
3241	ENDDO
3242	c
3243	IF (nqmax.GE.3) THEN
3244	DO iq = 3, nqmax
3245	DO k = 1, klev
3246	DO i = 1, klon
3247	d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
3248	ENDDO
3249	ENDDO
3250	ENDDO
3251	ENDIF
3252	c
3253	cIM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
3254	cIM global posePB#include "write_bilKP_ins.h"
3255	cIM global posePB#include "write_bilKP_ave.h"
3256	c
3257	c Sauvegarder les valeurs de t et q a la fin de la physique:
3258	c
3259	DO k = 1, klev
3260	DO i = 1, klon
3261	t_ancien(i,k) = t_seri(i,k)
3262	q_ancien(i,k) = q_seri(i,k)
3263	ENDDO
3264	ENDDO
3265	c
3266	!==========================================================================
3267	! Sorties des tendances pour un point particulier
3268	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
3269	! pour le debug
3270	! La valeur de igout est attribuee plus haut dans le programme
3271	!==========================================================================
3272
3273	if (prt_level.ge.1) then
3274	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
3275	write(lunout,*)
3276	s 'nlon,nlev,nqmax,debut,lafin,rjourvrai,gmtime,pdtphys pct tlos'
3277	write(lunout,*)
3278	s nlon,nlev,nqmax,debut,lafin,rjourvrai,gmtime,pdtphys,
3279	s pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce),
3280	s pctsrf(igout,is_sic)
3281	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
3282	do k=1,nlev
3283	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k),
3284	s d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k),
3285	s d_t_eva(igout,k)
3286	enddo
3287	write(lunout,*) 'cool,heat'
3288	do k=1,nlev
3289	write(lunout,*) cool(igout,k),heat(igout,k)
3290	enddo
3291
3292	write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
3293	do k=1,nlev
3294	write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k),
3295	s d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
3296	enddo
3297
3298	write(lunout,*) 'd_ps ',d_ps(igout)
3299	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
3300	do k=1,nlev
3301	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k),
3302	s d_qx(igout,k,1),d_qx(igout,k,2)
3303	enddo
3304	endif
3305
3306	!==========================================================================
3307
3308	c============================================================
3309	c Calcul de la temperature potentielle
3310	c============================================================
3311	DO k = 1, klev
3312	DO i = 1, klon
3313	theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
3314	ENDDO
3315	ENDDO
3316	c
3317
3318	c 22.03.04 BEG
3319	c=============================================================
3320	c Ecriture des sorties
3321	c=============================================================
3322	#ifdef CPP_IOIPSL
3323
3324	c Recupere des varibles calcule dans differents modules
3325	c pour ecriture dans histxxx.nc
3326
3327	! Get some variables from module fonte_neige_mod
3328	CALL fonte_neige_get_vars(pctsrf,
3329	. zxfqcalving, zxfqfonte, zxffonte)
3330
3331	IF (ocean == 'slab') THEN
3332	! Get some variables from module ocean_slab_mod
3333	CALL ocean_slab_get_vars(tslab, seaice, fluxo, fluxg)
3334	ELSEIF (ocean == 'couple') THEN
3335	! Get some variables from module ocean_cpl_mod
3336	CALL ocean_cpl_get_vars(fluxo, fluxg)
3337	ELSE
3338	! Get some variables from module ocean_forced_mod
3339	CALL ocean_forced_get_vars(fluxo, fluxg)
3340	ENDIF
3341
3342
3343	c Commente par abderrahmane le 11 2 08
3344	c#ifdef histhf
3345	c#include "write_histhf.h"
3346	c#endif
3347
3348	c#ifdef histday
3349	c#include "write_histday.h"
3350	c#endif
3351
3352	c#ifdef histmth
3353	c#include "write_histmth.h"
3354	c#endif
3355
3356	c#ifdef histins
3357	c#include "write_histins.h"
3358	c#endif
3359
3360	#include "phys_output_write.h"
3361
3362	#ifdef histISCCP
3363	#include "write_histISCCP.h"
3364	#endif
3365
3366	#ifdef histmthNMC
3367	#include "write_histmthNMC.h"
3368	#endif
3369
3370	#include "write_histday_seri.h"
3371
3372	#include "write_paramLMDZ_phy.h"
3373
3374	#endif
3375
3376	c 22.03.04 END
3377	c
3378	c====================================================================
3379	c Si c'est la fin, il faut conserver l'etat de redemarrage
3380	c====================================================================
3381	c
3382
3383
3384	IF (lafin) THEN
3385	itau_phy = itau_phy + itap
3386	CALL phyredem ("restartphy.nc",dtime,radpas,ocean,
3387	. rlat, rlon, pctsrf, ftsol,
3388	. falb1, falb2, rain_fall,
3389	. snow_fall,
3390	. solsw, sollw,
3391	. radsol,
3392	. zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,
3393	. t_ancien, q_ancien, rnebcon, ratqs, clwcon,
3394	. pbl_tke, zmax0, f0, ema_work1, ema_work2)
3395	open(97,form="unformatted",file="finbin")
3396	write(97) u_seri,v_seri,t_seri,q_seri
3397	close(97)
3398	ENDIF
3399
3400
3401	RETURN
3402	END
3403	FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
3404	IMPLICIT none
3405	c
3406	c Calculer et imprimer l'eau totale. A utiliser pour verifier
3407	c la conservation de l'eau
3408	c
3409	#include "YOMCST.h"
3410	INTEGER klon,klev
3411	REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
3412	REAL aire(klon)
3413	REAL qtotal, zx, qcheck
3414	INTEGER i, k
3415	c
3416	zx = 0.0
3417	DO i = 1, klon
3418	zx = zx + aire(i)
3419	ENDDO
3420	qtotal = 0.0
3421	DO k = 1, klev
3422	DO i = 1, klon
3423	qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i)
3424	. *(paprs(i,k)-paprs(i,k+1))/RG
3425	ENDDO
3426	ENDDO
3427	c
3428	qcheck = qtotal/zx
3429	c
3430	RETURN
3431	END
3432	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
3433	IMPLICIT none
3434	c
3435	c Tranformer une variable de la grille physique a
3436	c la grille d'ecriture
3437	c
3438	INTEGER nfield,nlon,iim,jjmp1, jjm
3439	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
3440	c
3441	INTEGER i, n, ig
3442	c
3443	jjm = jjmp1 - 1
3444	DO n = 1, nfield
3445	DO i=1,iim
3446	ecrit(i,n) = fi(1,n)
3447	ecrit(i+jjm*iim,n) = fi(nlon,n)
3448	ENDDO
3449	DO ig = 1, nlon - 2
3450	ecrit(iim+ig,n) = fi(1+ig,n)
3451	ENDDO
3452	ENDDO
3453	RETURN
3454	END

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

Download in other formats:

Original Format