Context Navigation

physiq.F @ 996

Last change on this file since 996 was 996, checked in by lsce, 16 years ago

Modifications liées au calcul des nouveau sous-fractions
Nettoyage de ocean slab : il reste uniquement la version avec glace de mer forcé
Nouveaux variables pour distiguer la version et type d'ocean : type_ocean=force/slab/couple, version_ocean=opa8/nemo pour couplé ou version_ocean=sicOBS pour slab

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

File size: 109.5 KB

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

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

Download in other formats:

Original Format