Context Navigation

physiq.F @ 942

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

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

Download in other formats:

Original Format