Context Navigation

physiq.F @ 889

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

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

Download in other formats:

Original Format