Context Navigation

physiq.F @ 927

Last change on this file since 927 was 927, checked in by lmdzadmin, 17 years ago
Ajout variables zmax0, f0 dans le startphy.nc FH IM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 124.3 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 ( jj_nb, jjmp1, nqmax,
1879	& nlevSTD, clevSTD, nbteta, ctetaSTD,
1880	& dtime, presnivs, ok_veget, ocean, iflag_pbl,
1881	& ok_mensuel,ok_journe,ok_hf,ok_instan,
1882	& rlon, rlat, nid_files)
1883
1884
1885	#ifdef histISCCP
1886	#include "ini_histISCCP.h"
1887	#endif
1888
1889	#ifdef histmthNMC
1890	#include "ini_histmthNMC.h"
1891	#endif
1892
1893	#include "ini_histday_seri.h"
1894
1895	#include "ini_paramLMDZ_phy.h"
1896
1897	#endif
1898
1899	cXXXPB Positionner date0 pour initialisation de ORCHIDEE
1900	date0 = zjulian
1901	C date0 = day_ini
1902	WRITE(,) 'physiq date0 : ',date0
1903	c
1904	c
1905	c
1906	c Prescrire l'ozone dans l'atmosphere
1907	c
1908	c
1909	cc DO i = 1, klon
1910	cc DO k = 1, klev
1911	cc CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1912	cc ENDDO
1913	cc ENDDO
1914	c
1915	#ifdef INCA
1916	call VTe(VTphysiq)
1917	call VTb(VTinca)
1918	iii = MOD(NINT(xjour),360)
1919	calday = FLOAT(iii) + gmtime
1920	WRITE(lunout,*) 'initial time ', xjour, calday
1921	#ifdef INCAINFO
1922	WRITE(lunout,*) 'Appel CHEMINI ...'
1923	#endif
1924	CALL chemini(
1925	$ rg,
1926	$ ra,
1927	$ airephy,
1928	$ rlat,
1929	$ rlon,
1930	$ presnivs,
1931	$ calday,
1932	$ klon,
1933	$ nqmax,
1934	$ pdtphys,
1935	$ annee_ref,
1936	$ day_ini)
1937	#ifdef INCAINFO
1938	WRITE(lunout,*) 'OK.'
1939	#endif
1940	call VTe(VTinca)
1941	call VTb(VTphysiq)
1942	#endif
1943	c
1944	ENDIF
1945	c
1946	c ************** Fin de IF ( debut ) *************
1947	c
1948
1949	! Tendances bidons pour les processus qui n'affectent pas certaines
1950	! variables.
1951	du0(:,:)=0.
1952	dv0(:,:)=0.
1953	dq0(:,:)=0.
1954	dql0(:,:)=0.
1955	c
1956	c Mettre a zero des variables de sortie (pour securite)
1957	c
1958	DO i = 1, klon
1959	d_ps(i) = 0.0
1960	ENDDO
1961	DO k = 1, klev
1962	DO i = 1, klon
1963	d_t(i,k) = 0.0
1964	d_u(i,k) = 0.0
1965	d_v(i,k) = 0.0
1966	ENDDO
1967	ENDDO
1968	DO iq = 1, nqmax
1969	DO k = 1, klev
1970	DO i = 1, klon
1971	d_qx(i,k,iq) = 0.0
1972	ENDDO
1973	ENDDO
1974	ENDDO
1975	da(:,:)=0.
1976	mp(:,:)=0.
1977	phi(:,:,:)=0.
1978	c
1979	c Ne pas affecter les valeurs entrees de u, v, h, et q
1980	c
1981	DO k = 1, klev
1982	DO i = 1, klon
1983	t_seri(i,k) = t(i,k)
1984	u_seri(i,k) = u(i,k)
1985	v_seri(i,k) = v(i,k)
1986	q_seri(i,k) = qx(i,k,ivap)
1987	ql_seri(i,k) = qx(i,k,iliq)
1988	qs_seri(i,k) = 0.
1989	ENDDO
1990	ENDDO
1991	IF (nqmax.GE.3) THEN
1992	DO iq = 3, nqmax
1993	DO k = 1, klev
1994	DO i = 1, klon
1995	tr_seri(i,k,iq-2) = qx(i,k,iq)
1996	ENDDO
1997	ENDDO
1998	ENDDO
1999	ELSE
2000	DO k = 1, klev
2001	DO i = 1, klon
2002	tr_seri(i,k,1) = 0.0
2003	ENDDO
2004	ENDDO
2005	ENDIF
2006	C
2007	DO i = 1, klon
2008	ztsol(i) = 0.
2009	ENDDO
2010	DO nsrf = 1, nbsrf
2011	DO i = 1, klon
2012	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
2013	ENDDO
2014	ENDDO
2015	cIM
2016	IF (ip_ebil_phy.ge.1) THEN
2017	ztit='after dynamic'
2018	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
2019	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2020	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2021	C Comme les tendances de la physique sont ajoute dans la dynamique,
2022	C on devrait avoir que la variation d'entalpie par la dynamique
2023	C est egale a la variation de la physique au pas de temps precedent.
2024	C Donc la somme de ces 2 variations devrait etre nulle.
2025	call diagphy(airephy,ztit,ip_ebil_phy
2026	e , zero_v, zero_v, zero_v, zero_v, zero_v
2027	e , zero_v, zero_v, zero_v, ztsol
2028	e , d_h_vcol+d_h_vcol_phy, d_qt, 0.
2029	s , fs_bound, fq_bound )
2030	END IF
2031
2032	c Diagnostiquer la tendance dynamique
2033	c
2034	IF (ancien_ok) THEN
2035	DO k = 1, klev
2036	DO i = 1, klon
2037	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
2038	d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
2039	ENDDO
2040	ENDDO
2041	ELSE
2042	DO k = 1, klev
2043	DO i = 1, klon
2044	d_t_dyn(i,k) = 0.0
2045	d_q_dyn(i,k) = 0.0
2046	ENDDO
2047	ENDDO
2048	ancien_ok = .TRUE.
2049	ENDIF
2050	c
2051	c Ajouter le geopotentiel du sol:
2052	c
2053	DO k = 1, klev
2054	DO i = 1, klon
2055	zphi(i,k) = pphi(i,k) + pphis(i)
2056	ENDDO
2057	ENDDO
2058	c
2059	c Verifier les temperatures
2060	c
2061	cIM BEG
2062	IF (check) THEN
2063	amn=MIN(ftsol(1,is_ter),1000.)
2064	amx=MAX(ftsol(1,is_ter),-1000.)
2065	DO i=2, klon
2066	amn=MIN(ftsol(i,is_ter),amn)
2067	amx=MAX(ftsol(i,is_ter),amx)
2068	ENDDO
2069	c
2070	PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
2071	ENDIF !(check) THEN
2072	cIM END
2073	c
2074	CALL hgardfou(t_seri,ftsol,'debutphy')
2075	c
2076	cIM BEG
2077	IF (check) THEN
2078	amn=MIN(ftsol(1,is_ter),1000.)
2079	amx=MAX(ftsol(1,is_ter),-1000.)
2080	DO i=2, klon
2081	amn=MIN(ftsol(i,is_ter),amn)
2082	amx=MAX(ftsol(i,is_ter),amx)
2083	ENDDO
2084	c
2085	PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
2086	ENDIF !(check) THEN
2087	cIM END
2088	c
2089	c Incrementer le compteur de la physique
2090	c
2091	itap = itap + 1
2092	julien = MOD(NINT(xjour),360)
2093	if (julien .eq. 0) julien = 360
2094	c
2095	c Mettre en action les conditions aux limites (albedo, sst, etc.).
2096	c Prescrire l'ozone et calculer l'albedo sur l'ocean.
2097	c
2098	IF (MOD(itap-1,lmt_pas) .EQ. 0) THEN
2099	WRITE(lunout,*)' PHYS cond julien ',julien
2100	CALL ozonecm( FLOAT(julien), rlat, paprs, wo)
2101	ENDIF
2102	c
2103	c Re-evaporer l'eau liquide nuageuse
2104	c
2105	DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse
2106	DO i = 1, klon
2107	zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
2108	c zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
2109	zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
2110	zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
2111	zb = MAX(0.0,ql_seri(i,k))
2112	za = - MAX(0.0,ql_seri(i,k))
2113	. * (zlvdcp(1.-zdelta)+zlsdcpzdelta)
2114	t_seri(i,k) = t_seri(i,k) + za
2115	q_seri(i,k) = q_seri(i,k) + zb
2116	ql_seri(i,k) = 0.0
2117	d_t_eva(i,k) = za
2118	d_q_eva(i,k) = zb
2119	ENDDO
2120	ENDDO
2121	cIM
2122	IF (ip_ebil_phy.ge.2) THEN
2123	ztit='after reevap'
2124	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,1,dtime
2125	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2126	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2127	call diagphy(airephy,ztit,ip_ebil_phy
2128	e , zero_v, zero_v, zero_v, zero_v, zero_v
2129	e , zero_v, zero_v, zero_v, ztsol
2130	e , d_h_vcol, d_qt, d_ec
2131	s , fs_bound, fq_bound )
2132	C
2133	END IF
2134
2135	c
2136	c=========================================================================
2137	! Calculs de l'orbite.
2138	! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
2139	! doit donc etre placé avant radlwsw et pbl_surface
2140
2141	! choix entre calcul de la longitude solaire vraie ou valeur fixee a
2142	! solarlong0
2143
2144	if (solarlong0<-999.) then
2145	CALL orbite(FLOAT(julien),zlongi,dist)
2146	else
2147	zlongi=solarlong0 ! longitude solaire vraie
2148	dist=1. ! distance au soleil / moyenne
2149	endif
2150
2151	print*,'Longitude solaire ',zlongi,solarlong0
2152
2153	! Avec ou sans cycle diurne
2154	IF (cycle_diurne) THEN
2155	zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
2156	CALL zenang(zlongi,gmtime,zdtime,rlat,rlon,rmu0,fract)
2157	ELSE
2158	rmu0 = -999.999
2159	ENDIF
2160
2161	if (mydebug) then
2162	call writefield_phy('u_seri',u_seri,llm)
2163	call writefield_phy('v_seri',v_seri,llm)
2164	call writefield_phy('t_seri',t_seri,llm)
2165	call writefield_phy('q_seri',q_seri,llm)
2166	endif
2167
2168	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2169	c Appel au pbl_surface : Planetary Boudary Layer et Surface
2170	c Cela implique tous les interactions des sous-surfaces et la partie diffusion
2171	c turbulent du couche limit.
2172	c
2173	c Certains varibales de sorties de pbl_surface sont utiliser que pour
2174	c ecriture des fihiers hist_XXXX.nc, ces sont :
2175	c qsol, zq2m, s_pblh, s_lcl,
2176	c s_capCL, s_oliqCL, s_cteiCL,s_pblT,
2177	c s_therm, s_trmb1, s_trmb2, s_trmb3,
2178	c zxrugs, zu10m, zv10m, fder,
2179	c zxqsurf, rh2m, zxfluxu, zxfluxv,
2180	c frugs, agesno, fsollw, fsolsw,
2181	c d_ts, fevap, fluxlat, t2m,
2182	c wfbils, wfbilo, fluxt, fluxu, fluxv,
2183	c
2184	c Certains ne sont pas utiliser du tout :
2185	c dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
2186	c
2187
2188	CALL pbl_surface(
2189	e dtime, date0, itap, julien,
2190	e debut, lafin,
2191	e rlon, rlat, rugoro, rmu0,
2192	e rain_fall, snow_fall, solsw, sollw,
2193	e t_seri, q_seri, u_seri, v_seri,
2194	e pplay, paprs, pctsrf,
2195	+ ftsol, falb1, falb2, u10m, v10m,
2196	s sollwdown, cdragh, cdragm, yu1, yv1,
2197	s albsol1, albsol2, sens, evap,
2198	s zxtsol, zxfluxlat, zt2m, qsat2m,
2199	s d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf,
2200	s ycoefh, pctsrf_new,
2201	d qsol, zq2m, s_pblh, s_lcl,
2202	d s_capCL, s_oliqCL, s_cteiCL,s_pblT,
2203	d s_therm, s_trmb1, s_trmb2, s_trmb3,
2204	d zxrugs, zu10m, zv10m, fder,
2205	d zxqsurf, rh2m, zxfluxu, zxfluxv,
2206	d frugs, agesno, fsollw, fsolsw,
2207	d d_ts, fevap, fluxlat, t2m,
2208	d wfbils, wfbilo, fluxt, fluxu, fluxv,
2209	- dsens, devap, zxsnow,
2210	- zxfluxt, zxfluxq, q2m, fluxq, pbl_tke )
2211	c
2212	c
2213	pctsrf(:,:) = pctsrf_new(:,:)
2214
2215	!-----------------------------------------------------------------------------------------
2216	! ajout des tendances de la diffusion turbulente
2217	CALL add_phys_tend(d_u_vdf,d_v_vdf,d_t_vdf,d_q_vdf,dql0,'vdf')
2218	!-----------------------------------------------------------------------------------------
2219
2220	if (mydebug) then
2221	call writefield_phy('u_seri',u_seri,llm)
2222	call writefield_phy('v_seri',v_seri,llm)
2223	call writefield_phy('t_seri',t_seri,llm)
2224	call writefield_phy('q_seri',q_seri,llm)
2225	endif
2226
2227
2228	IF (ip_ebil_phy.ge.2) THEN
2229	ztit='after surface_main'
2230	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2231	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2232	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2233	call diagphy(airephy,ztit,ip_ebil_phy
2234	e , zero_v, zero_v, zero_v, zero_v, sens
2235	e , evap , zero_v, zero_v, ztsol
2236	e , d_h_vcol, d_qt, d_ec
2237	s , fs_bound, fq_bound )
2238	END IF
2239
2240	c =================================================================== c
2241	c Calcul de Qsat
2242
2243	DO k = 1, klev
2244	DO i = 1, klon
2245	zx_t = t_seri(i,k)
2246	IF (thermcep) THEN
2247	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2248	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2249	zx_qs = MIN(0.5,zx_qs)
2250	zcor = 1./(1.-retv*zx_qs)
2251	zx_qs = zx_qs*zcor
2252	ELSE
2253	IF (zx_t.LT.t_coup) THEN
2254	zx_qs = qsats(zx_t)/pplay(i,k)
2255	ELSE
2256	zx_qs = qsatl(zx_t)/pplay(i,k)
2257	ENDIF
2258	ENDIF
2259	zqsat(i,k)=zx_qs
2260	ENDDO
2261	ENDDO
2262
2263	if (1.eq.1) then
2264	write(lunout,*) 'L qsat (g/kg) avant clouds_gno'
2265	write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
2266	endif
2267	c
2268	c Appeler la convection (au choix)
2269	c
2270	DO k = 1, klev
2271	DO i = 1, klon
2272	conv_q(i,k) = d_q_dyn(i,k)
2273	. + d_q_vdf(i,k)/dtime
2274	conv_t(i,k) = d_t_dyn(i,k)
2275	. + d_t_vdf(i,k)/dtime
2276	ENDDO
2277	ENDDO
2278	IF (check) THEN
2279	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2280	WRITE(lunout,*) "avantcon=", za
2281	ENDIF
2282	zx_ajustq = .FALSE.
2283	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
2284	IF (zx_ajustq) THEN
2285	DO i = 1, klon
2286	z_avant(i) = 0.0
2287	ENDDO
2288	DO k = 1, klev
2289	DO i = 1, klon
2290	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k))
2291	. *(paprs(i,k)-paprs(i,k+1))/RG
2292	ENDDO
2293	ENDDO
2294	ENDIF
2295	IF (iflag_con.EQ.1) THEN
2296	stop'reactiver le call conlmd dans physiq.F'
2297	c CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
2298	c . d_t_con, d_q_con,
2299	c . rain_con, snow_con, ibas_con, itop_con)
2300	ELSE IF (iflag_con.EQ.2) THEN
2301	CALL conflx(dtime, paprs, pplay, t_seri, q_seri,
2302	e conv_t, conv_q, -evap, omega,
2303	s d_t_con, d_q_con, rain_con, snow_con,
2304	s pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
2305	s kcbot, kctop, kdtop, pmflxr, pmflxs)
2306	WHERE (rain_con < 0.) rain_con = 0.
2307	WHERE (snow_con < 0.) snow_con = 0.
2308	DO i = 1, klon
2309	ibas_con(i) = klev+1 - kcbot(i)
2310	itop_con(i) = klev+1 - kctop(i)
2311	ENDDO
2312	ELSE IF (iflag_con.GE.3) THEN
2313	c nb of tracers for the KE convection:
2314	c MAF la partie traceurs est faite dans phytrac
2315	c on met ntra=1 pour limiter les appels mais on peut
2316	c supprimer les calculs / ftra.
2317	ntra = 1
2318
2319	c=====================================================================================
2320	cajout pour la parametrisation des poches froides:
2321	ccalcul de t_wake et t_undi: si pas de poches froides, t_wake=t_undi=t_seri
2322	do k=1,klev
2323	do i=1,klon
2324	if (iflag_wake.eq.1) then
2325	t_wake(i,k) = t_seri(i,k)
2326	. +(1-wake_s(i))*wake_deltat(i,k)
2327	q_wake(i,k) = q_seri(i,k)
2328	. +(1-wake_s(i))*wake_deltaq(i,k)
2329	t_undi(i,k) = t_seri(i,k)
2330	. -wake_s(i)*wake_deltat(i,k)
2331	q_undi(i,k) = q_seri(i,k)
2332	. -wake_s(i)*wake_deltaq(i,k)
2333	else
2334	t_wake(i,k) = t_seri(i,k)
2335	q_wake(i,k) = q_seri(i,k)
2336	t_undi(i,k) = t_seri(i,k)
2337	q_undi(i,k) = q_seri(i,k)
2338	endif
2339	enddo
2340	enddo
2341
2342	cc-- Calcul de l'energie disponible ALE (J/kg) et de la puissance disponible ALP (W/m2)
2343	cc-- pour le soulevement des particules dans le modele convectif
2344	c
2345	do i = 1,klon
2346	ALE(i) = 0.
2347	ALP(i) = 0.
2348	enddo
2349	c
2350	ccalcul de ale_wake et alp_wake
2351	do i = 1,klon
2352	if (iflag_wake.eq.1) then
2353	ale_wake(i) = 0.5wake_cstar(i)*2
2354	alp_wake(i) = wake_fip(i)
2355	else
2356	ale_wake(i) = 0.
2357	alp_wake(i) = 0.
2358	endif
2359	enddo
2360	ccombinaison avec ale et alp de couche limite: constantes si pas de couplage, valeurs calculees
2361	cdans le thermique sinon
2362	if (iflag_coupl.eq.0) then
2363	print*,'ALE et ALP imposes'
2364	do i = 1,klon
2365	con ne couple que ale
2366	c ALE(i) = max(ale_wake(i),Ale_bl(i))
2367	ALE(i) = max(ale_wake(i),ale_bl_prescr)
2368	con ne couple que alp
2369	c ALP(i) = alp_wake(i) + Alp_bl(i)
2370	ALP(i) = alp_wake(i) + alp_bl_prescr
2371	enddo
2372	else
2373	print*,'ALE et ALP couples au thermique'
2374	do i = 1,klon
2375	ALE(i) = max(ale_wake(i),Ale_bl(i))
2376	ALP(i) = alp_wake(i) + Alp_bl(i)
2377	c write(20,*)'ALE',ALE(i),Ale_bl(i),ale_wake(i)
2378	c write(21,*)'ALP',ALP(i),Alp_bl(i),alp_wake(i)
2379	enddo
2380	endif
2381
2382	cfin calcul ale et alp
2383	c=================================================================================================
2384
2385
2386	c sb, oct02:
2387	c Schema de convection modularise et vectorise:
2388	c (driver commun aux versions 3 et 4)
2389	c
2390	IF (ok_cvl) THEN ! new driver for convectL
2391
2392	CALL concvl (iflag_con,iflag_clos,
2393	. dtime,paprs,pplay,t_undi,q_undi,
2394	. t_wake,q_wake,
2395	. u_seri,v_seri,tr_seri,nbtr,
2396	. ALE,ALP,
2397	. ema_work1,ema_work2,
2398	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2399	. rain_con, snow_con, ibas_con, itop_con, sigd,
2400	. upwd,dnwd,dnwd0,
2401	. Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl,
2402	. pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd,
2403	. pmflxr,pmflxs,da,phi,mp,
2404	. ftd,fqd,lalim_conv,wght_th)
2405
2406	cIM cf. FH
2407	clwcon0=qcondc
2408	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
2409
2410	ELSE ! ok_cvl
2411	c MAF conema3 ne contient pas les traceurs
2412	CALL conema3 (dtime,
2413	. paprs,pplay,t_seri,q_seri,
2414	. u_seri,v_seri,tr_seri,ntra,
2415	. ema_work1,ema_work2,
2416	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2417	. rain_con, snow_con, ibas_con, itop_con,
2418	. upwd,dnwd,dnwd0,bas,top,
2419	. Ma,cape,tvp,rflag,
2420	. pbase
2421	. ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr
2422	. ,clwcon0)
2423
2424	ENDIF ! ok_cvl
2425
2426	c
2427	c Correction precip
2428	rain_con = rain_con * cvl_corr
2429	snow_con = snow_con * cvl_corr
2430	c
2431
2432	IF (.NOT. ok_gust) THEN
2433	do i = 1, klon
2434	wd(i)=0.0
2435	enddo
2436	ENDIF
2437
2438	c =================================================================== c
2439	c Calcul des proprietes des nuages convectifs
2440	c
2441
2442	c calcul des proprietes des nuages convectifs
2443	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
2444	call clouds_gno
2445	s (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
2446
2447	c =================================================================== c
2448
2449	DO i = 1, klon
2450	ema_pcb(i) = pbase(i)
2451	ENDDO
2452	DO i = 1, klon
2453
2454	! L'idicage de itop_con peut cacher un pb potentiel
2455	! FH sous la dictee de JYG, CR
2456	ema_pct(i) = paprs(i,itop_con(i)+1)
2457
2458	if (itop_con(i).gt.klev-3) then
2459	print*,'La convection monte trop haut '
2460	print*,'itop_con(,',i,',)=',itop_con(i)
2461	endif
2462	ENDDO
2463	DO i = 1, klon
2464	ema_cbmf(i) = ema_workcbmf(i)
2465	ENDDO
2466	ELSE IF (iflag_con.eq.0) THEN
2467	write(lunout,*) 'On n appelle pas la convection'
2468	clwcon0=0.
2469	rnebcon0=0.
2470	d_t_con=0.
2471	d_q_con=0.
2472	d_u_con=0.
2473	d_v_con=0.
2474	rain_con=0.
2475	snow_con=0.
2476	bas=1
2477	top=1
2478	ELSE
2479	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
2480	CALL abort
2481	ENDIF
2482
2483	c CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
2484	c . d_u_con, d_v_con)
2485
2486	!-----------------------------------------------------------------------------------------
2487	! ajout des tendances de la diffusion turbulente
2488	CALL add_phys_tend(d_u_con,d_v_con,d_t_con,d_q_con,dql0,'con')
2489	!-----------------------------------------------------------------------------------------
2490
2491	if (mydebug) then
2492	call writefield_phy('u_seri',u_seri,llm)
2493	call writefield_phy('v_seri',v_seri,llm)
2494	call writefield_phy('t_seri',t_seri,llm)
2495	call writefield_phy('q_seri',q_seri,llm)
2496	endif
2497
2498	cIM
2499	IF (ip_ebil_phy.ge.2) THEN
2500	ztit='after convect'
2501	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2502	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2503	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2504	call diagphy(airephy,ztit,ip_ebil_phy
2505	e , zero_v, zero_v, zero_v, zero_v, zero_v
2506	e , zero_v, rain_con, snow_con, ztsol
2507	e , d_h_vcol, d_qt, d_ec
2508	s , fs_bound, fq_bound )
2509	END IF
2510	C
2511	IF (check) THEN
2512	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2513	WRITE(lunout,*)"aprescon=", za
2514	zx_t = 0.0
2515	za = 0.0
2516	DO i = 1, klon
2517	za = za + airephy(i)/FLOAT(klon)
2518	zx_t = zx_t + (rain_con(i)+
2519	. snow_con(i))*airephy(i)/FLOAT(klon)
2520	ENDDO
2521	zx_t = zx_t/za*dtime
2522	WRITE(lunout,*)"Precip=", zx_t
2523	ENDIF
2524	IF (zx_ajustq) THEN
2525	DO i = 1, klon
2526	z_apres(i) = 0.0
2527	ENDDO
2528	DO k = 1, klev
2529	DO i = 1, klon
2530	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k))
2531	. *(paprs(i,k)-paprs(i,k+1))/RG
2532	ENDDO
2533	ENDDO
2534	DO i = 1, klon
2535	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime)
2536	. /z_apres(i)
2537	ENDDO
2538	DO k = 1, klev
2539	DO i = 1, klon
2540	IF (z_factor(i).GT.(1.0+1.0E-08) .OR.
2541	. z_factor(i).LT.(1.0-1.0E-08)) THEN
2542	q_seri(i,k) = q_seri(i,k) * z_factor(i)
2543	ENDIF
2544	ENDDO
2545	ENDDO
2546	ENDIF
2547	zx_ajustq=.FALSE.
2548
2549	c
2550	c=============================================================================
2551	cRR:Evolution de la poche froide: on ne fait pas de separation wake/env
2552	cpour la couche limite diffuse pour l instant
2553	c
2554	if (iflag_wake.eq.1) then
2555	DO k=1,klev
2556	DO i=1,klon
2557	dt_dwn(i,k) = ftd(i,k)
2558	wdt_PBL(i,k) = 0.
2559	dq_dwn(i,k) = fqd(i,k)
2560	wdq_PBL(i,k) = 0.
2561	M_dwn(i,k) = dnwd0(i,k)
2562	M_up(i,k) = upwd(i,k)
2563	dt_a(i,k) = d_t_con(i,k)/dtime - ftd(i,k)
2564	udt_PBL(i,k) = 0.
2565	dq_a(i,k) = d_q_con(i,k)/dtime - fqd(i,k)
2566	udq_PBL(i,k) = 0.
2567	ENDDO
2568	ENDDO
2569	c
2570	ccalcul caracteristiques de la poche froide
2571	call calWAKE (paprs,pplay,dtime
2572	: ,t_seri,q_seri,omega,ibas_con
2573	: ,dt_dwn,dq_dwn,M_dwn,M_up
2574	: ,dt_a,dq_a,sigd
2575	: ,wdt_PBL,wdq_PBL
2576	: ,udt_PBL,udq_PBL
2577	o ,wake_deltat,wake_deltaq,wake_dth
2578	o ,wake_h,wake_s,wake_dens
2579	o ,wake_pe,wake_fip,wake_gfl
2580	o ,dt_wake,dq_wake
2581	o ,wake_k, t_undi,q_undi
2582	o ,wake_omgbdth,wake_dp_omgb
2583	o ,wake_dtKE,wake_dqKE
2584	o ,wake_dtPBL,wake_dqPBL
2585	o ,wake_omg,wake_dp_deltomg
2586	o ,wake_spread,wake_Cstar,wake_d_deltat_gw
2587	o ,wake_ddeltat,wake_ddeltaq)
2588	c
2589	!-----------------------------------------------------------------------------------------
2590	! ajout des tendances des poches froides
2591	! Faire rapidement disparaitre l'ancien dt_wake pour garder un d_t_wake
2592	! coherent avec les autres d_t_...
2593	d_t_wake(:,:)=dt_wake(:,:)*dtime
2594	d_q_wake(:,:)=dq_wake(:,:)*dtime
2595	CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,'wake')
2596	!-----------------------------------------------------------------------------------------
2597
2598	endif
2599	c print*,'apres callwake iflag_cldcon=', iflag_cldcon
2600	c
2601	c===================================================================
2602	c Convection seche (thermiques ou ajustement)
2603	c===================================================================
2604	c
2605	call stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri
2606	s ,seuil_inversion,weak_inversion,dthmin)
2607
2608
2609
2610	d_t_ajsb(:,:)=0.
2611	d_q_ajsb(:,:)=0.
2612	d_t_ajs(:,:)=0.
2613	d_u_ajs(:,:)=0.
2614	d_v_ajs(:,:)=0.
2615	d_q_ajs(:,:)=0.
2616	clwcon0th(:,:)=0.
2617	fm_therm(:,:)=0.
2618	entr_therm(:,:)=0.
2619	c
2620	IF(prt_level>9)WRITE(lunout,*)
2621	. 'AVANT LA CONVECTION SECHE , iflag_thermals='
2622	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2623	if(iflag_thermals.lt.0) then
2624	c Rien
2625	c ====
2626	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
2627
2628
2629	else
2630
2631	c Thermiques
2632	c ==========
2633	IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals='
2634	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2635	print*,'JUSTE AVANT , iflag_thermals='
2636	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2637
2638
2639	if (iflag_thermals.gt.1) then
2640	call calltherm(pdtphys
2641	s ,pplay,paprs,pphi,weak_inversion
2642	s ,u_seri,v_seri,t_seri,q_seri,zqsat,debut
2643	s ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
2644	s ,fm_therm,entr_therm,zqasc,clwcon0th,lmax_th,ratqscth
2645	s ,ratqsdiff,zqsatth
2646	con rajoute ale et alp, et les caracteristiques de la couche alim
2647	s ,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0)
2648	endif
2649
2650	! call calltherm(pdtphys
2651	! s ,pplay,paprs,pphi
2652	! s ,u_seri,v_seri,t_seri,q_seri
2653	! s ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
2654	! s ,fm_therm,entr_therm)
2655
2656	c Ajustement sec
2657	c ==============
2658
2659	! Dans le cas où on active les thermiques, on fait partir l'ajustement
2660	! a partir du sommet des thermiques.
2661	! Dans le cas contraire, on demarre au niveau 1.
2662
2663	if (iflag_thermals.ge.13.or.iflag_thermals.eq.0) then
2664
2665	if(iflag_thermals.eq.0) then
2666	IF(prt_level>9)WRITE(lunout,*)'ajsec'
2667	limbas(:)=1
2668	else
2669	limbas(:)=lmax_th(:)
2670	endif
2671
2672	! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
2673	! pour des test de convergence numerique.
2674	! Le nouveau ajsec est a priori mieux, meme pour le cas
2675	! iflag_thermals = 0 (l'ancienne version peut faire des tendances
2676	! non nulles numeriquement pour des mailles non concernees.
2677
2678	if (iflag_thermals.eq.0) then
2679	CALL ajsec_convV2(paprs, pplay, t_seri,q_seri
2680	s , d_t_ajsb, d_q_ajsb)
2681	else
2682	CALL ajsec(paprs, pplay, t_seri,q_seri,limbas
2683	s , d_t_ajsb, d_q_ajsb)
2684	endif
2685
2686	!-----------------------------------------------------------------------------------------
2687	! ajout des tendances de l'ajustement sec ou des thermiques
2688	CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,'ajsb')
2689	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
2690	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
2691
2692	!-----------------------------------------------------------------------------------------
2693
2694	endif
2695
2696	endif
2697	c
2698	c===================================================================
2699	cIM
2700	IF (ip_ebil_phy.ge.2) THEN
2701	ztit='after dry_adjust'
2702	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2703	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2704	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2705	END IF
2706
2707
2708	c-------------------------------------------------------------------------
2709	c Caclul des ratqs
2710	c-------------------------------------------------------------------------
2711
2712	c print*,'calcul des ratqs'
2713	c ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
2714	c ----------------
2715	c on ecrase le tableau ratqsc calcule par clouds_gno
2716	if (iflag_cldcon.eq.1) then
2717	do k=1,klev
2718	do i=1,klon
2719	if(ptconv(i,k)) then
2720	ratqsc(i,k)=ratqsbas
2721	s +fact_cldcon*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
2722	else
2723	ratqsc(i,k)=0.
2724	endif
2725	enddo
2726	enddo
2727	else if (iflag_cldcon.eq.4) then
2728	ptconvth(:,:)=.false.
2729	ratqsc(:,:)=0.
2730	print*,'avant clouds_gno thermique'
2731	call clouds_gno
2732	s (klon,klev,q_seri,zqsat,clwcon0th,ptconvth,ratqsc,rnebcon0th)
2733	print*,' CLOUDS_GNO OK'
2734	endif
2735
2736	c ratqs stables
2737	c -------------
2738
2739	if (iflag_ratqs.eq.0) then
2740
2741	! Le cas iflag_ratqs=0 correspond a la version IPCC 2005 du modele.
2742	do k=1,klev
2743	do i=1, klon
2744	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2745	s min((paprs(i,1)-pplay(i,k))/(paprs(i,1)-30000.),1.)
2746	enddo
2747	enddo
2748
2749	! Pour iflag_ratqs=1 ou 2, le ratqs est constant au dessus de
2750	! 300 hPa (ratqshaut), varie lineariement en fonction de la pression
2751	! entre 600 et 300 hPa et est soit constant (ratqsbas) pour iflag_ratqs=1
2752	! soit lineaire (entre 0 a la surface et ratqsbas) pour iflag_ratqs=2
2753	! Il s'agit de differents tests dans la phase de reglage du modele
2754	! avec thermiques.
2755
2756	else if (iflag_ratqs.eq.1) then
2757
2758	do k=1,klev
2759	do i=1, klon
2760	if (pplay(i,k).ge.60000.) then
2761	ratqss(i,k)=ratqsbas
2762	else if ((pplay(i,k).ge.30000.).and.
2763	s (pplay(i,k).lt.60000.)) then
2764	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2765	s (60000.-pplay(i,k))/(60000.-30000.)
2766	else
2767	ratqss(i,k)=ratqshaut
2768	endif
2769	enddo
2770	enddo
2771
2772	else
2773
2774	do k=1,klev
2775	do i=1, klon
2776	if (pplay(i,k).ge.60000.) then
2777	ratqss(i,k)=ratqsbas
2778	s *(paprs(i,1)-pplay(i,k))/(paprs(i,1)-60000.)
2779	else if ((pplay(i,k).ge.30000.).and.
2780	s (pplay(i,k).lt.60000.)) then
2781	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2782	s (60000.-pplay(i,k))/(60000.-30000.)
2783	else
2784	ratqss(i,k)=ratqshaut
2785	endif
2786	enddo
2787	enddo
2788	endif
2789
2790
2791
2792
2793	c ratqs final
2794	c -----------
2795
2796	if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2
2797	s .or.iflag_cldcon.eq.4) then
2798
2799	! On ajoute une constante au ratqsc*2 pour tenir compte de
2800	! fluctuations turbulentes de petite echelle
2801
2802	do k=1,klev
2803	do i=1,klon
2804	if ((fm_therm(i,k).gt.1.e-10)) then
2805	ratqsc(i,k)=sqrt(ratqsc(i,k)2+0.052)
2806	endif
2807	enddo
2808	enddo
2809
2810	! les ratqs sont une conbinaison de ratqss et ratqsc
2811	! 1800s, en dur pour le moment, est le temps de
2812	! relaxation des ratqs
2813
2814	facteur=exp(-pdtphys/1800.)
2815
2816	print*,'WARNING ratqs a revoir '
2817	ratqs(:,:)=ratqsc(:,:)(1.-facteur)+ratqs(:,:)facteur
2818	ratqs(:,:)=sqrt(ratqs(:,:)2+ratqss(:,:)2)
2819
2820	else
2821	! on ne prend que le ratqs stable pour fisrtilp
2822	ratqs(:,:)=ratqss(:,:)
2823	endif
2824
2825
2826	c
2827	c Appeler le processus de condensation a grande echelle
2828	c et le processus de precipitation
2829	c-------------------------------------------------------------------------
2830	CALL fisrtilp(dtime,paprs,pplay,
2831	. t_seri, q_seri,ptconv,ratqs,
2832	. d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq,
2833	. rain_lsc, snow_lsc,
2834	. pfrac_impa, pfrac_nucl, pfrac_1nucl,
2835	. frac_impa, frac_nucl,
2836	. prfl, psfl, rhcl)
2837
2838	WHERE (rain_lsc < 0) rain_lsc = 0.
2839	WHERE (snow_lsc < 0) snow_lsc = 0.
2840	!-----------------------------------------------------------------------------------------
2841	! ajout des tendances de la diffusion turbulente
2842	CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,'lsc')
2843	!-----------------------------------------------------------------------------------------
2844	DO k = 1, klev
2845	DO i = 1, klon
2846	cldfra(i,k) = rneb(i,k)
2847	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
2848	ENDDO
2849	ENDDO
2850	IF (check) THEN
2851	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2852	WRITE(lunout,*)"apresilp=", za
2853	zx_t = 0.0
2854	za = 0.0
2855	DO i = 1, klon
2856	za = za + airephy(i)/FLOAT(klon)
2857	zx_t = zx_t + (rain_lsc(i)
2858	. + snow_lsc(i))*airephy(i)/FLOAT(klon)
2859	ENDDO
2860	zx_t = zx_t/za*dtime
2861	WRITE(lunout,*)"Precip=", zx_t
2862	ENDIF
2863	cIM
2864	IF (ip_ebil_phy.ge.2) THEN
2865	ztit='after fisrt'
2866	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2867	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2868	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2869	call diagphy(airephy,ztit,ip_ebil_phy
2870	e , zero_v, zero_v, zero_v, zero_v, zero_v
2871	e , zero_v, rain_lsc, snow_lsc, ztsol
2872	e , d_h_vcol, d_qt, d_ec
2873	s , fs_bound, fq_bound )
2874	END IF
2875
2876	if (mydebug) then
2877	call writefield_phy('u_seri',u_seri,llm)
2878	call writefield_phy('v_seri',v_seri,llm)
2879	call writefield_phy('t_seri',t_seri,llm)
2880	call writefield_phy('q_seri',q_seri,llm)
2881	endif
2882
2883	c
2884	c-------------------------------------------------------------------
2885	c PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
2886	c-------------------------------------------------------------------
2887
2888	c 1. NUAGES CONVECTIFS
2889	c
2890	cIM cf FH
2891	c IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
2892	IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
2893	snow_tiedtke=0.
2894	c print*,'avant calcul de la pseudo precip '
2895	c print*,'iflag_cldcon',iflag_cldcon
2896	if (iflag_cldcon.eq.-1) then
2897	rain_tiedtke=rain_con
2898	else
2899	c print*,'calcul de la pseudo precip '
2900	rain_tiedtke=0.
2901	c print*,'calcul de la pseudo precip 0'
2902	do k=1,klev
2903	do i=1,klon
2904	if (d_q_con(i,k).lt.0.) then
2905	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys
2906	s *(paprs(i,k)-paprs(i,k+1))/rg
2907	endif
2908	enddo
2909	enddo
2910	endif
2911	c
2912	c call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
2913	c
2914
2915	c Nuages diagnostiques pour Tiedtke
2916	CALL diagcld1(paprs,pplay,
2917	cIM cf FH . rain_con,snow_con,ibas_con,itop_con,
2918	. rain_tiedtke,snow_tiedtke,ibas_con,itop_con,
2919	. diafra,dialiq)
2920	DO k = 1, klev
2921	DO i = 1, klon
2922	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2923	cldliq(i,k) = dialiq(i,k)
2924	cldfra(i,k) = diafra(i,k)
2925	ENDIF
2926	ENDDO
2927	ENDDO
2928
2929	ELSE IF (iflag_cldcon.ge.3) THEN
2930	c On prend pour les nuages convectifs le max du calcul de la
2931	c convection et du calcul du pas de temps precedent diminue d'un facteur
2932	c facttemps
2933	facteur = pdtphys *facttemps
2934	do k=1,klev
2935	do i=1,klon
2936	rnebcon(i,k)=rnebcon(i,k)*facteur
2937	if (rnebcon0(i,k)clwcon0(i,k).gt.rnebcon(i,k)clwcon(i,k))
2938	s then
2939	rnebcon(i,k)=rnebcon0(i,k)
2940	clwcon(i,k)=clwcon0(i,k)
2941	endif
2942	enddo
2943	enddo
2944
2945	c
2946	cjq - introduce the aerosol direct and first indirect radiative forcings
2947	cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
2948	IF (ok_ade.OR.ok_aie) THEN
2949	#if defined(CPP_COUPLE) && !defined(INCA_AER)
2950	! Get sulfate aerosol distribution
2951	CALL readsulfate(rjourvrai, debut, sulfate)
2952	CALL readsulfate_preind(rjourvrai, debut, sulfate_pi)
2953
2954	! Calculate aerosol optical properties (Olivier Boucher)
2955	CALL aeropt(pplay, paprs, t_seri, sulfate, rhcl,
2956	. tau_ae, piz_ae, cg_ae, aerindex)
2957	#endif
2958	#if !defined(CPP_COUPLE)
2959	! Get sulfate aerosol distribution
2960	CALL readsulfate(rjourvrai, debut, sulfate)
2961	CALL readsulfate_preind(rjourvrai, debut, sulfate_pi)
2962
2963	! Calculate aerosol optical properties (Olivier Boucher)
2964	CALL aeropt(pplay, paprs, t_seri, sulfate, rhcl,
2965	. tau_ae, piz_ae, cg_ae, aerindex)
2966	#endif
2967	cym
2968	ELSE
2969	tau_ae(:,:,:)=0.0
2970	piz_ae(:,:,:)=0.0
2971	cg_ae(:,:,:)=0.0
2972	cym
2973	ENDIF
2974
2975	c
2976	cIM calcul nuages par le simulateur ISCCP
2977	c
2978	#ifdef histISCCP
2979	IF (ok_isccp) THEN
2980	cIM appel simulateur toutes les NINT(freq_ISCCP/dtime) heures
2981	IF (MOD(itap,NINT(freq_ISCCP/dtime)).EQ.0) THEN
2982	#include "calcul_simulISCCP.h"
2983	ENDIF !(MOD(itap,NINT(freq_ISCCP/dtime))
2984	ENDIF !ok_isccp
2985	#endif
2986
2987	c On prend la somme des fractions nuageuses et des contenus en eau
2988	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
2989	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
2990
2991	ENDIF
2992
2993	c
2994	c 2. NUAGES STARTIFORMES
2995	c
2996	IF (ok_stratus) THEN
2997	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
2998	DO k = 1, klev
2999	DO i = 1, klon
3000	IF (diafra(i,k).GT.cldfra(i,k)) THEN
3001	cldliq(i,k) = dialiq(i,k)
3002	cldfra(i,k) = diafra(i,k)
3003	ENDIF
3004	ENDDO
3005	ENDDO
3006	ENDIF
3007	c
3008	c Precipitation totale
3009	c
3010	DO i = 1, klon
3011	rain_fall(i) = rain_con(i) + rain_lsc(i)
3012	snow_fall(i) = snow_con(i) + snow_lsc(i)
3013	ENDDO
3014	cIM
3015	IF (ip_ebil_phy.ge.2) THEN
3016	ztit="after diagcld"
3017	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
3018	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3019	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3020	END IF
3021	c
3022	c Calculer l'humidite relative pour diagnostique
3023	c
3024	DO k = 1, klev
3025	DO i = 1, klon
3026	zx_t = t_seri(i,k)
3027	IF (thermcep) THEN
3028	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
3029	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
3030	zx_qs = MIN(0.5,zx_qs)
3031	zcor = 1./(1.-retv*zx_qs)
3032	zx_qs = zx_qs*zcor
3033	ELSE
3034	IF (zx_t.LT.t_coup) THEN
3035	zx_qs = qsats(zx_t)/pplay(i,k)
3036	ELSE
3037	zx_qs = qsatl(zx_t)/pplay(i,k)
3038	ENDIF
3039	ENDIF
3040	zx_rh(i,k) = q_seri(i,k)/zx_qs
3041	zqsat(i,k)=zx_qs
3042	ENDDO
3043	ENDDO
3044
3045	cIM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
3046	c equivalente a 2m (tpote) pour diagnostique
3047	c
3048	DO i = 1, klon
3049	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
3050	IF (thermcep) THEN
3051	IF(zt2m(i).LT.RTT) then
3052	Lheat=RLSTT
3053	ELSE
3054	Lheat=RLVTT
3055	ENDIF
3056	ELSE
3057	IF (zt2m(i).LT.RTT) THEN
3058	Lheat=RLSTT
3059	ELSE
3060	Lheat=RLVTT
3061	ENDIF
3062	ENDIF
3063	tpote(i) = tpot(i)*
3064	. EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
3065	ENDDO
3066
3067
3068	#ifdef INCA
3069	call VTe(VTphysiq)
3070	call VTb(VTinca)
3071	calday = FLOAT(julien) + gmtime
3072
3073	#ifdef INCA_AER
3074	call AEROSOL_METEO_CALC(calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs
3075	& ,prfl,psfl,pctsrf,airephy,xjour,rlat,rlon,u10m,v10m)
3076	#endif
3077
3078	#ifdef INCAINFO
3079	WRITE(lunout,*)'Appel CHEMHOOK_BEGIN ...'
3080	#endif
3081
3082	zxsnow_dummy(:) = 0.0
3083
3084	CALL chemhook_begin (calday,
3085	#if defined(INCA) && !defined(INCA_CH4) && !defined(INCA_NMHC) && !defined(INCA_AER)
3086	$ julien,
3087	$ gmtime,
3088	#endif
3089	$ pctsrf(1,1),
3090	$ rlat,
3091	$ rlon,
3092	$ airephy,
3093	$ paprs,
3094	$ pplay,
3095	$ ycoefh,
3096	$ pphi,
3097	$ t_seri,
3098	$ u,
3099	$ v,
3100	$ wo,
3101	$ q_seri,
3102	$ zxtsol,
3103	$ zxsnow_dummy,
3104	$ solsw,
3105	$ albsol1,
3106	$ rain_fall,
3107	$ snow_fall,
3108	$ itop_con,
3109	$ ibas_con,
3110	$ cldfra,
3111	$ iim,
3112	$ jjm,
3113	#ifdef INCA_AER
3114	$ tr_seri,
3115	$ ftsol,
3116	$ paprs,
3117	$ cdragh,
3118	$ cdragm,
3119	$ pctsrf,
3120	$ pdtphys,
3121	$ itap)
3122	#else
3123	$ tr_seri)
3124	#endif
3125
3126
3127	#ifdef INCAINFO
3128	WRITE(lunout,*)'OK.'
3129	#endif
3130	call VTe(VTinca)
3131	call VTb(VTphysiq)
3132	#endif
3133	c
3134	c Calculer les parametres optiques des nuages et quelques
3135	c parametres pour diagnostiques:
3136	c
3137	if (ok_newmicro) then
3138	CALL newmicro (paprs, pplay,ok_newmicro,
3139	. t_seri, cldliq, cldfra, cldtau, cldemi,
3140	. cldh, cldl, cldm, cldt, cldq,
3141	. flwp, fiwp, flwc, fiwc,
3142	e ok_aie,
3143	#if defined(CPP_COUPLE) && defined(INCA_AER)
3144	e ccm(:,:,1), ccm(:,:,2),
3145	#else
3146	e sulfate, sulfate_pi,
3147	#endif
3148	e bl95_b0, bl95_b1,
3149	s cldtaupi, re, fl)
3150	else
3151	CALL nuage (paprs, pplay,
3152	. t_seri, cldliq, cldfra, cldtau, cldemi,
3153	. cldh, cldl, cldm, cldt, cldq,
3154	e ok_aie,
3155	e sulfate, sulfate_pi,
3156	e bl95_b0, bl95_b1,
3157	s cldtaupi, re, fl)
3158
3159	endif
3160	c
3161	c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
3162	c
3163	IF (MOD(itaprad,radpas).EQ.0) THEN
3164
3165	DO i = 1, klon
3166	albsol1(i) = falb1(i,is_oce) * pctsrf(i,is_oce)
3167	. + falb1(i,is_lic) * pctsrf(i,is_lic)
3168	. + falb1(i,is_ter) * pctsrf(i,is_ter)
3169	. + falb1(i,is_sic) * pctsrf(i,is_sic)
3170	albsol2(i) = falb2(i,is_oce) * pctsrf(i,is_oce)
3171	. + falb2(i,is_lic) * pctsrf(i,is_lic)
3172	. + falb2(i,is_ter) * pctsrf(i,is_ter)
3173	. + falb2(i,is_sic) * pctsrf(i,is_sic)
3174	ENDDO
3175
3176	if (mydebug) then
3177	call writefield_phy('u_seri',u_seri,llm)
3178	call writefield_phy('v_seri',v_seri,llm)
3179	call writefield_phy('t_seri',t_seri,llm)
3180	call writefield_phy('q_seri',q_seri,llm)
3181	endif
3182
3183
3184	#if defined(CPP_COUPLE) && defined(INCA_AER)
3185	CALL radlwsw_inca ! nouveau rayonnement (compatible Arpege-IFS)
3186	e (kdlon,kflev,dist, rmu0, fract,
3187	e paprs, pplay,zxtsol,albsol, albsollw, t_seri,q_seri,
3188	e wo,
3189	e cldfra, cldemi, cldtau,
3190	s heat,heat0,cool,cool0,radsol,albpla,
3191	s topsw,toplw,solsw,sollw,
3192	s sollwdown,
3193	s topsw0,toplw0,solsw0,sollw0,
3194	s lwdn0, lwdn, lwup0, lwup,
3195	s swdn0, swdn, swup0, swup,
3196	e ok_ade, ok_aie, ! new for aerosol radiative effects
3197	e tau_inca, piz_inca, cg_inca, ! ="=
3198	s topswad_inca, solswad_inca, ! ="=
3199	s topswad0_inca, solswad0_inca, ! ="=
3200	s topsw_inca, topsw0_inca,
3201	s solsw_inca, solsw0_inca,
3202	e cldtaupi, ! ="=
3203	s topswai_inca, solswai_inca) ! ="=
3204	#else
3205	CALL radlwsw ! nouveau rayonnement (compatible Arpege-IFS)
3206	e (dist, rmu0, fract,
3207	e paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri,
3208	e wo,
3209	e cldfra, cldemi, cldtau,
3210	s heat,heat0,cool,cool0,radsol,albpla,
3211	s topsw,toplw,solsw,sollw,
3212	s sollwdown,
3213	s topsw0,toplw0,solsw0,sollw0,
3214	s lwdn0, lwdn, lwup0, lwup,
3215	s swdn0, swdn, swup0, swup,
3216	e ok_ade, ok_aie, ! new for aerosol radiative effects
3217	e tau_ae, piz_ae, cg_ae, ! ="=
3218	s topswad, solswad, ! ="=
3219	e cldtaupi, ! ="=
3220	s topswai, solswai) ! ="=
3221	#endif
3222	itaprad = 0
3223	ENDIF
3224	itaprad = itaprad + 1
3225
3226	if (iflag_radia.eq.0) then
3227	print *,'--------------------------------------------------'
3228	print *,'>>>> ATTENTION rayonnement desactive pour ce cas'
3229	print *,'>>>> heat et cool mis a zero '
3230	print *,'--------------------------------------------------'
3231	heat=0.
3232	cool=0.
3233	endif
3234
3235
3236	c
3237	c Ajouter la tendance des rayonnements (tous les pas)
3238	c
3239	DO k = 1, klev
3240	DO i = 1, klon
3241	t_seri(i,k) = t_seri(i,k)
3242	. + (heat(i,k)-cool(i,k)) * dtime/86400.
3243	ENDDO
3244	ENDDO
3245	c
3246	if (mydebug) then
3247	call writefield_phy('u_seri',u_seri,llm)
3248	call writefield_phy('v_seri',v_seri,llm)
3249	call writefield_phy('t_seri',t_seri,llm)
3250	call writefield_phy('q_seri',q_seri,llm)
3251	endif
3252
3253	cIM
3254	IF (ip_ebil_phy.ge.2) THEN
3255	ztit='after rad'
3256	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
3257	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3258	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3259	call diagphy(airephy,ztit,ip_ebil_phy
3260	e , topsw, toplw, solsw, sollw, zero_v
3261	e , zero_v, zero_v, zero_v, ztsol
3262	e , d_h_vcol, d_qt, d_ec
3263	s , fs_bound, fq_bound )
3264	END IF
3265	c
3266	c
3267	c Calculer l'hydrologie de la surface
3268	c
3269	c CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
3270	c . agesno, ftsol,fqsurf,fsnow, ruis)
3271	c
3272
3273	c
3274	c Calculer le bilan du sol et la derive de temperature (couplage)
3275	c
3276	DO i = 1, klon
3277	c bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
3278	c a la demande de JLD
3279	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
3280	ENDDO
3281	c
3282	cmoddeblott(jan95)
3283	c Appeler le programme de parametrisation de l'orographie
3284	c a l'echelle sous-maille:
3285	c
3286	IF (ok_orodr) THEN
3287	c
3288	c selection des points pour lesquels le shema est actif:
3289	igwd=0
3290	DO i=1,klon
3291	itest(i)=0
3292	c IF ((zstd(i).gt.10.0)) THEN
3293	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
3294	itest(i)=1
3295	igwd=igwd+1
3296	idx(igwd)=i
3297	ENDIF
3298	ENDDO
3299	c igwdim=MAX(1,igwd)
3300	c
3301	CALL drag_noro(klon,klev,dtime,paprs,pplay,
3302	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
3303	e igwd,idx,itest,
3304	e t_seri, u_seri, v_seri,
3305	cIM 141004 s zulow, zvlow, zustr, zvstr,
3306	s zulow, zvlow, zustrdr, zvstrdr,
3307	s d_t_oro, d_u_oro, d_v_oro)
3308	c
3309	c ajout des tendances
3310	!-----------------------------------------------------------------------------------------
3311	! ajout des tendances de la trainee de l'orographie
3312	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,'oro')
3313	!-----------------------------------------------------------------------------------------
3314	c
3315	ENDIF ! fin de test sur ok_orodr
3316	c
3317	if (mydebug) then
3318	call writefield_phy('u_seri',u_seri,llm)
3319	call writefield_phy('v_seri',v_seri,llm)
3320	call writefield_phy('t_seri',t_seri,llm)
3321	call writefield_phy('q_seri',q_seri,llm)
3322	endif
3323
3324	IF (ok_orolf) THEN
3325	c
3326	c selection des points pour lesquels le shema est actif:
3327	igwd=0
3328	DO i=1,klon
3329	itest(i)=0
3330	IF ((zpic(i)-zmea(i)).GT.100.) THEN
3331	itest(i)=1
3332	igwd=igwd+1
3333	idx(igwd)=i
3334	ENDIF
3335	ENDDO
3336	c igwdim=MAX(1,igwd)
3337	c
3338	CALL lift_noro(klon,klev,dtime,paprs,pplay,
3339	e rlat,zmea,zstd,zpic,
3340	e itest,
3341	e t_seri, u_seri, v_seri,
3342	s zulow, zvlow, zustrli, zvstrli,
3343	s d_t_lif, d_u_lif, d_v_lif)
3344	c
3345	!-----------------------------------------------------------------------------------------
3346	! ajout des tendances de la portance de l'orographie
3347	CALL add_phys_tend(d_u_lif,d_v_lif,d_t_lif,dq0,dql0,'lif')
3348	!-----------------------------------------------------------------------------------------
3349	c
3350	ENDIF ! fin de test sur ok_orolf
3351	c
3352	cIM cf. FLott BEG
3353	C STRESS NECESSAIRES: TOUTE LA PHYSIQUE
3354
3355	if (mydebug) then
3356	call writefield_phy('u_seri',u_seri,llm)
3357	call writefield_phy('v_seri',v_seri,llm)
3358	call writefield_phy('t_seri',t_seri,llm)
3359	call writefield_phy('q_seri',q_seri,llm)
3360	endif
3361
3362	DO i = 1, klon
3363	zustrph(i)=0.
3364	zvstrph(i)=0.
3365	ENDDO
3366	DO k = 1, klev
3367	DO i = 1, klon
3368	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
3369	c (paprs(i,k)-paprs(i,k+1))/rg
3370	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
3371	c (paprs(i,k)-paprs(i,k+1))/rg
3372	ENDDO
3373	ENDDO
3374	c
3375	cIM calcul composantes axiales du moment angulaire et couple des montagnes
3376	c
3377	IF (is_sequential) THEN
3378
3379	CALL aaam_bud (27,klon,klev,rjourvrai,gmtime,
3380	C ra,rg,romega,
3381	C rlat,rlon,pphis,
3382	C zustrdr,zustrli,zustrph,
3383	C zvstrdr,zvstrli,zvstrph,
3384	C paprs,u,v,
3385	C aam, torsfc)
3386	ENDIF
3387	cIM cf. FLott END
3388	cIM
3389	IF (ip_ebil_phy.ge.2) THEN
3390	ztit='after orography'
3391	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
3392	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3393	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3394	END IF
3395	c
3396	c
3397	cAA
3398	cAA Installation de l'interface online-offline pour traceurs
3399	cAA
3400	c====================================================================
3401	c Calcul des tendances traceurs
3402	c====================================================================
3403	C
3404	call phytrac ( rnpb,
3405	I itap,
3406	I julien,
3407	I gmtime,
3408	I debut,
3409	I lafin,
3410	I nqmax-2,
3411	I nlon,
3412	I nlev,
3413	I dtime,
3414	I u,
3415	I v,
3416	I t,
3417	I paprs,
3418	I pplay,
3419	I pmfu,
3420	I pmfd,
3421	I pen_u,
3422	I pde_u,
3423	I pen_d,
3424	I pde_d,
3425	I ycoefh,
3426	I fm_therm,
3427	I entr_therm,
3428	I yu1,
3429	I yv1,
3430	I ftsol,
3431	I pctsrf,
3432	I rlat,
3433	I frac_impa,
3434	I frac_nucl,
3435	I rlon,
3436	I presnivs,
3437	I pphis,
3438	I pphi,
3439	I albsol1,
3440	I qx(1,1,1),
3441	I rhcl,
3442	I cldfra,
3443	I rneb,
3444	I diafra,
3445	I cldliq,
3446	I itop_con,
3447	I ibas_con,
3448	I pmflxr,
3449	I pmflxs,
3450	I prfl,
3451	I psfl,
3452	I da,
3453	I phi,
3454	I mp,
3455	I upwd,
3456	I dnwd,
3457	#ifdef INCA
3458	I flxmass_w,
3459	#if defined(INCA_AER) && defined(CPP_COUPLE)
3460	I tau_inca,
3461	I piz_inca,
3462	I cg_inca,
3463	I ccm,
3464	I rfname,
3465	#endif
3466	#endif
3467	O tr_seri)
3468
3469	IF (offline) THEN
3470
3471	print*,'Attention on met a 0 les thermiques pour phystoke'
3472	call phystokenc (
3473	I nlon,nlev,pdtphys,rlon,rlat,
3474	I t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
3475	I fm_therm,entr_therm,
3476	I ycoefh,yu1,yv1,ftsol,pctsrf,
3477	I frac_impa, frac_nucl,
3478	I pphis,airephy,dtime,itap)
3479
3480
3481	ENDIF
3482
3483	c
3484	c Calculer le transport de l'eau et de l'energie (diagnostique)
3485	c
3486	CALL transp (paprs,zxtsol,
3487	e t_seri, q_seri, u_seri, v_seri, zphi,
3488	s ve, vq, ue, uq)
3489	c
3490	cIM global posePB BEG
3491	IF(1.EQ.0) THEN
3492	c
3493	CALL transp_lay (paprs,zxtsol,
3494	e t_seri, q_seri, u_seri, v_seri, zphi,
3495	s ve_lay, vq_lay, ue_lay, uq_lay)
3496	c
3497	ENDIF !(1.EQ.0) THEN
3498	cIM global posePB END
3499	c Accumuler les variables a stocker dans les fichiers histoire:
3500	c
3501	c+jld ec_conser
3502	DO k = 1, klev
3503	DO i = 1, klon
3504	ZRCPD = RCPD(1.0+RVTMP2q_seri(i,k))
3505	d_t_ec(i,k)=0.5/ZRCPD
3506	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
3507	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
3508	d_t_ec(i,k) = d_t_ec(i,k)/dtime
3509	END DO
3510	END DO
3511	c-jld ec_conser
3512	cIM
3513	IF (ip_ebil_phy.ge.1) THEN
3514	ztit='after physic'
3515	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
3516	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3517	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3518	C Comme les tendances de la physique sont ajoute dans la dynamique,
3519	C on devrait avoir que la variation d'entalpie par la dynamique
3520	C est egale a la variation de la physique au pas de temps precedent.
3521	C Donc la somme de ces 2 variations devrait etre nulle.
3522	call diagphy(airephy,ztit,ip_ebil_phy
3523	e , topsw, toplw, solsw, sollw, sens
3524	e , evap, rain_fall, snow_fall, ztsol
3525	e , d_h_vcol, d_qt, d_ec
3526	s , fs_bound, fq_bound )
3527	C
3528	d_h_vcol_phy=d_h_vcol
3529	C
3530	END IF
3531	C
3532	c=======================================================================
3533	c SORTIES
3534	c=======================================================================
3535
3536	cIM Interpolation sur les niveaux de pression du NMC
3537	c -------------------------------------------------
3538	c
3539	#include "calcul_STDlev.h"
3540	c
3541	c slp sea level pressure
3542	slp(:) = paprs(:,1)exp(pphis(:)/(RDt_seri(:,1)))
3543	c
3544	ccc prw = eau precipitable
3545	DO i = 1, klon
3546	prw(i) = 0.
3547	DO k = 1, klev
3548	prw(i) = prw(i) +
3549	. q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
3550	ENDDO
3551	ENDDO
3552	c
3553	cIM initialisation + calculs divers diag AMIP2
3554	c
3555	#include "calcul_divers.h"
3556	c
3557	#ifdef INCA
3558	call VTe(VTphysiq)
3559	call VTb(VTinca)
3560	#ifdef INCAINFO
3561	WRITE(lunout,*)'Appel CHEMHOOK_END ...'
3562	#endif
3563	CALL chemhook_end (calday,
3564	$ dtime,
3565	$ pplay,
3566	$ t_seri,
3567	$ tr_seri,
3568	$ nbtr,
3569	$ paprs,
3570	$ q_seri,
3571	$ annee_ref,
3572	$ day_ini,
3573	$ airephy,
3574	#ifdef INCA_AER
3575	$ xjour,
3576	$ pphi,
3577	$ pphis,
3578	$ zx_rh)
3579	#else
3580	$ xjour)
3581	#endif
3582	#ifdef INCAINFO
3583	WRITE(lunout,*)'OK.'
3584	#endif
3585	call VTe(VTinca)
3586	call VTb(VTphysiq)
3587	#endif
3588
3589	c=============================================================
3590	c
3591	c Convertir les incrementations en tendances
3592	c
3593	if (mydebug) then
3594	call writefield_phy('u_seri',u_seri,llm)
3595	call writefield_phy('v_seri',v_seri,llm)
3596	call writefield_phy('t_seri',t_seri,llm)
3597	call writefield_phy('q_seri',q_seri,llm)
3598	endif
3599
3600	DO k = 1, klev
3601	DO i = 1, klon
3602	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
3603	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
3604	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
3605	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
3606	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
3607	ENDDO
3608	ENDDO
3609	c
3610	IF (nqmax.GE.3) THEN
3611	DO iq = 3, nqmax
3612	DO k = 1, klev
3613	DO i = 1, klon
3614	d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
3615	ENDDO
3616	ENDDO
3617	ENDDO
3618	ENDIF
3619	c
3620	cIM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
3621	cIM global posePB#include "write_bilKP_ins.h"
3622	cIM global posePB#include "write_bilKP_ave.h"
3623	c
3624	c Sauvegarder les valeurs de t et q a la fin de la physique:
3625	c
3626	DO k = 1, klev
3627	DO i = 1, klon
3628	t_ancien(i,k) = t_seri(i,k)
3629	q_ancien(i,k) = q_seri(i,k)
3630	ENDDO
3631	ENDDO
3632	c
3633	!==========================================================================
3634	! Sorties des tendances pour un point particulier
3635	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
3636	! pour le debug
3637	! La valeur de igout est attribuee plus haut dans le programme
3638	!==========================================================================
3639
3640	if (klon.eq.1) then
3641	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
3642	write(lunout,*)
3643	s 'nlon,nlev,nqmax,debut,lafin,rjourvrai,gmtime,pdtphys'
3644	write(lunout,*)
3645	s nlon,nlev,nqmax,debut,lafin,rjourvrai,gmtime,pdtphys
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