Context Navigation

physiq.F @ 973

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

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

Download in other formats:

Original Format