Context Navigation

physiq.F @ 929

Last change on this file since 929 was 929, checked in by lmdzadmin, 17 years ago
Ecriture plus compacte des sorties et ajout variables LMDZ AI IM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 124.4 KB

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

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

Download in other formats:

Original Format