Context Navigation

physiq.F @ 726

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

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

Download in other formats:

Original Format