Context Navigation

physiq.F @ 844

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

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

Download in other formats:

Original Format