Context Navigation

physiq.F @ 943

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

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

Download in other formats:

Original Format