Context Navigation

physiq.F @ 839

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

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

Download in other formats:

Original Format