Context Navigation

physiq.F @ 776

Last change on this file since 776 was 776, checked in by Laurent Fairhead, 17 years ago
Suite du merge entre la version et la HEAD: quelques modifications de Yann sur le LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 117.1 KB

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

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

Download in other formats:

Original Format