Context Navigation

physiq.F @ 1032

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

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

Download in other formats:

Original Format