Context Navigation

physiq.F @ 714

Last change on this file since 714 was 704, checked in by Laurent Fairhead, 19 years ago
Inclusion des modifs de Y. Meurdesoif pour la version V3 LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 117.2 KB

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

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

Download in other formats:

Original Format