Context Navigation

physiq.F @ 782

Last change on this file since 782 was 782, checked in by Laurent Fairhead, 17 years ago
Adaptation du code a la nouvelle interface avec les surface de Josefine LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 103.9 KB

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

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

Download in other formats:

Original Format