Context Navigation

physiq.F @ 937

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

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

Download in other formats:

Original Format