Context Navigation

physiq.F @ 766

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

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

Download in other formats:

Original Format