Context Navigation

physiq.F @ 888

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

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

Download in other formats:

Original Format