Context Navigation

physiq.F @ 879

Last change on this file since 879 was 879, checked in by Laurent Fairhead, 18 years ago
Suite de la bascule vers une physique avec thermiques, nouvelle convection, poche froide ... LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 126.4 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	c
690	REAL,allocatable,save :: falbe(:,:)
691	c$OMP THREADPRIVATE(falbe) ! albedo par type de surface
692	c
693	REAL,allocatable,save :: falblw(:,:)
694	c$OMP THREADPRIVATE(falblw) ! albedo par type de surface
695
696	c
697	c
698	c Parametres de l'Orographie a l'Echelle Sous-Maille (OESM):
699	c
700	REAL,allocatable,save :: zmea(:)
701	c$OMP THREADPRIVATE(zmea)
702	cym SAVE zmea ! orographie moyenne
703	c
704	REAL,allocatable,save :: zstd(:)
705	c$OMP THREADPRIVATE(zstd)
706	cym SAVE zstd ! deviation standard de l'OESM
707	c
708	REAL,allocatable,save :: zsig(:)
709	c$OMP THREADPRIVATE(zsig)
710	cym SAVE zsig ! pente de l'OESM
711	c
712	REAL,allocatable,save :: zgam(:)
713	c$OMP THREADPRIVATE(zgam)
714	cym save zgam ! anisotropie de l'OESM
715	c
716	REAL,allocatable,save :: zthe(:)
717	c$OMP THREADPRIVATE(zthe)
718	cym SAVE zthe ! orientation de l'OESM
719	c
720	REAL,allocatable,save :: zpic(:)
721	c$OMP THREADPRIVATE(zpic)
722	cym SAVE zpic ! Maximum de l'OESM
723	c
724	REAL,allocatable,save :: zval(:)
725	c$OMP THREADPRIVATE(zval)
726	cym SAVE zval ! Minimum de l'OESM
727	c
728	REAL,allocatable,save :: rugoro(:)
729	c$OMP THREADPRIVATE(rugoro)
730	cym SAVE rugoro ! longueur de rugosite de l'OESM
731	c
732	cIM 141004 REAL zulow(klon),zvlow(klon),zustr(klon), zvstr(klon)
733	REAL zulow(klon),zvlow(klon)
734	c
735	REAL,allocatable,save :: zuthe(:),zvthe(:)
736	c$OMP THREADPRIVATE(zuthe,zvthe)
737	cym SAVE zuthe
738	cym SAVE zvthe
739	INTEGER igwd,idx(klon),itest(klon)
740	c
741	REAL agesno(klon,nbsrf)
742	c
743	REAL,allocatable,save :: alb_neig(:)
744	c$OMP THREADPRIVATE(alb_neig)
745	cym SAVE alb_neig ! albedo de la neige
746	c
747	c REAL,allocatable,save :: run_off_lic_0(:)
748	cc$OMP THREADPRIVATE(run_off_lic_0)
749	cym SAVE run_off_lic_0
750	cKE43
751	c Variables liees a la convection de K. Emanuel (sb):
752	c
753	REAL,allocatable,save :: ema_workcbmf(:) ! cloud base mass flux
754	c$OMP THREADPRIVATE(ema_workcbmf)
755	cym SAVE ema_workcbmf
756
757	REAL,allocatable,save :: ema_cbmf(:) ! cloud base mass flux
758	c$OMP THREADPRIVATE(ema_cbmf)
759	cym SAVE ema_cbmf
760
761	REAL,allocatable,save :: ema_pcb(:) ! cloud base pressure
762	c$OMP THREADPRIVATE(ema_pcb)
763	cym SAVE ema_pcb
764
765	REAL,allocatable,save :: ema_pct(:) ! cloud top pressure
766	c$OMP THREADPRIVATE(ema_pct)
767	cym SAVE ema_pct
768
769	REAL bas, top ! cloud base and top levels
770	SAVE bas
771	SAVE top
772	c$OMP THREADPRIVATE(bas, top)
773
774	REAL,allocatable,save :: Ma(:,:) ! undilute upward mass flux
775	c$OMP THREADPRIVATE(Ma)
776	cym SAVE Ma
777	REAL,allocatable,save :: qcondc(:,:) ! in-cld water content from convect
778	c$OMP THREADPRIVATE(qcondc)
779	cym SAVE qcondc
780	REAL,allocatable,save :: ema_work1(:, :), ema_work2(:, :)
781	c$OMP THREADPRIVATE(ema_work1,ema_work2)
782	cym SAVE ema_work1, ema_work2
783	REAL wdn(klon), tdn(klon), qdn(klon)
784
785	REAL,allocatable,save :: wd(:) ! sb
786	c$OMP THREADPRIVATE(wd)
787	cym SAVE wd ! sb
788
789	REAL,allocatable,save :: sigd(:)
790	c
791	c=================================================================================================
792	cCR04.12.07: on ajoute les nouvelles variables du nouveau schema de convection avec poches froides
793	c Variables liées à la poche froide (jyg)
794
795	REAL mip(klon,klev) ! mass flux shed by the adiab ascent at each level
796	REAL Vprecip(klon,klev) ! precipitation vertical profile
797	REAL,allocatable,save :: cin(:) ! CIN
798	REAL,allocatable,save :: ftd(:,:) ! differential heating between wake and environment
799	REAL,allocatable,save :: fqd(:,:) ! differential moistening between wake and environment
800	c34EK
801	c -- Variables de controle de ALE et ALP
802	REAL,allocatable,save :: ALE(:) ! Energie disponible pour soulevement : utilisee par la
803	c convection d'Emanuel pour le declenchement et la regulation
804	REAL,allocatable,save :: ALP(:) ! Puissance disponible pour soulevement !
805	c
806	REAL wape_prescr, fip_prescr
807	INTEGER it_wape_prescr
808	SAVE wape_prescr, fip_prescr, it_wape_prescr
809	c
810	c variables supplementaires de concvl
811	REAL Tconv(klon,klev)
812	REAL ment(klon,klev,klev),sij(klon,klev,klev)
813	REAL dd_t(klon,klev),dd_q(klon,klev)
814	cnouvelles variables pour le couplage convection-couche limite
815	real,allocatable,save :: Ale_bl(:)
816	real,allocatable,save :: Alp_bl(:)
817	integer,allocatable,save :: lalim_conv(:)
818	real,allocatable,save :: wght_th(:,:)
819	real alp_bl_prescr
820	save alp_bl_prescr
821	real ale_bl_prescr
822	save ale_bl_prescr
823	real ale_wake(klon)
824	real alp_wake(klon)
825	cRC
826	c Variables liées à la poche froide (jyg et rr)
827	c Version diagnostique pour l'instant : pas de rétroaction sur la convection
828
829	REAL t_wake(klon,klev),q_wake(klon,klev) ! wake pour la convection
830
831	REAL,allocatable,save :: wake_deltat(:,:) ! Wake : ecart de temperature avec la
832	c zone non perturbee
833	REAL,allocatable,save :: wake_deltaq(:,:) ! Wake : ecart d'humidite avec la
834	c zone non perturbee
835	REAL wake_dth(klon,klev) ! wake : temp pot difference
836
837	REAL wake_d_deltat_gw(klon,klev)! wake : delta T tendency due to Gravity Wave (/s)
838	REAL wake_omgbdth(klon,klev) ! Wake : flux of Delta_Theta transported by LS omega
839	REAL wake_dp_omgb(klon,klev) ! Wake : vertical gradient of large scale omega
840	REAL wake_dtKE(klon,klev) ! Wake : differential heating (wake - unpertubed) CONV
841	REAL wake_dqKE(klon,klev) ! Wake : differential moistening (wake - unpertubed) CONV
842	REAL wake_dtPBL(klon,klev) ! Wake : differential heating (wake - unpertubed) PBL
843	REAL wake_dqPBL(klon,klev) ! Wake : differential moistening (wake - unpertubed) PBL
844	REAL wake_omg(klon,klev) ! Wake : velocity difference (wake - unpertubed)
845	REAL wake_ddeltat(klon,klev),wake_ddeltaq(klon,klev)
846	REAL wake_dp_deltomg(klon,klev) ! Wake : gradient vertical de wake_omg
847	REAL wake_spread(klon,klev) ! spreading term in wake_delt
848	REAL,allocatable,save :: wake_Cstar(:) ! vitesse d'etalement de la poche
849	cpourquoi y'a pas de save??
850	REAL wake_h(klon) ! Wake : hauteur de la poche froide
851	REAL,allocatable,save :: wake_s(:) ! Wake : fraction surfacique occupee par
852	c la poche froide
853	INTEGER wake_k(klon) ! Wake sommet
854	c
855	REAL t_undi(klon,klev) ! temperature moyenne dans la zone non perturbee
856	REAL q_undi(klon,klev) ! humidite moyenne dans la zone non perturbee
857	c
858	REAL wake_pe(klon) ! Wake potential energy - WAPE
859	REAL,allocatable,save :: wake_fip(:) ! Gust Front Impinging power - ALP
860
861	REAL wake_gfl(klon) ! Gust Front Length
862	REAL wake_dens(klon)
863	c
864	REAL,allocatable,save :: dt_wake(:,:)
865	REAL,allocatable,save :: dq_wake(:,:) ! LS tendencies due to wake
866	c
867	REAL dt_dwn(klon,klev)
868	REAL dq_dwn(klon,klev)
869	REAL wdt_PBL(klon,klev)
870	REAL udt_PBL(klon,klev)
871	REAL wdq_PBL(klon,klev)
872	REAL udq_PBL(klon,klev)
873	REAL M_dwn(klon,klev)
874	REAL M_up(klon,klev)
875	REAL dt_a(klon,klev)
876	REAL dq_a(klon,klev)
877	c
878	cRR:fin declarations poches froides
879	c=======================================================================================================
880
881	c Variables locales pour la couche limite (al1):
882	c
883	cAl1 REAL pblh(klon) ! Hauteur de couche limite
884	cAl1 SAVE pblh
885	c34EK
886	c
887	c Variables locales:
888	c
889	REAL cdragh(klon) ! drag coefficient pour T and Q
890	REAL cdragm(klon) ! drag coefficient pour vent
891	cAA
892	cAA Pour phytrac
893	cAA
894	REAL ycoefh(klon,klev) ! coef d'echange pour phytrac
895	REAL yu1(klon) ! vents dans la premiere couche U
896	REAL yv1(klon) ! vents dans la premiere couche V
897
898	REAL zxffonte(klon), zxfqcalving(klon),zxfqfonte(klon)
899
900	c@$$ LOGICAL offline ! Controle du stockage ds "physique"
901	c@$$ PARAMETER (offline=.false.)
902	c@$$ INTEGER physid
903	REAL,allocatable,save :: pfrac_impa(:,:)! Produits des coefs lessivage impaction
904	c$OMP THREADPRIVATE(pfrac_impa)
905	cym save pfrac_impa
906	REAL,allocatable,save :: pfrac_nucl(:,:)! Produits des coefs lessivage nucleation
907	c$OMP THREADPRIVATE(pfrac_nucl)
908	cym save pfrac_nucl
909	REAL,allocatable,save :: pfrac_1nucl(:,:)! Produits des coefs lessi nucl (alpha = 1)
910	c$OMP THREADPRIVATE(pfrac_1nucl)
911	cym save pfrac_1nucl
912	REAL frac_impa(klon,klev) ! fractions d'aerosols lessivees (impaction)
913	REAL frac_nucl(klon,klev) ! idem (nucleation)
914	#ifdef INCA
915	INTEGER :: iii
916	REAL :: calday
917	#endif
918
919	cAA
920	REAL,allocatable,save :: rain_fall(:) ! pluie
921	c$OMP THREADPRIVATE(rain_fall)
922	REAL,allocatable,save :: snow_fall(:) ! neige
923	c$OMP THREADPRIVATE(snow_fall)
924	cym save snow_fall, rain_fall
925
926	cIM cf FH pour Tiedtke 080604
927	REAL rain_tiedtke(klon),snow_tiedtke(klon)
928	c
929
930	REAL,allocatable,save :: total_rain(:), nday_rain(:)
931	c$OMP THREADPRIVATE(total_rain,nday_rain)
932	cym save total_rain, nday_rain
933	cIM 050204 END
934	REAL evap(klon), devap(klon) ! evaporation et sa derivee
935	REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
936
937	REAL bils(klon) ! bilan de chaleur au sol
938	REAL wfbilo(klon,nbsrf) ! bilan d'eau, pour chaque
939	C ! type de sous-surface et pondere par la fraction
940	REAL wfbils(klon,nbsrf) ! bilan de chaleur au sol, pour chaque
941	C ! type de sous-surface et pondere par la fraction
942	REAL fder(klon)
943	REAL ve(klon) ! integr. verticale du transport meri. de l'energie
944	REAL vq(klon) ! integr. verticale du transport meri. de l'eau
945	REAL ue(klon) ! integr. verticale du transport zonal de l'energie
946	REAL uq(klon) ! integr. verticale du transport zonal de l'eau
947	c
948	REAL frugs(klon,nbsrf)
949	REAL zxrugs(klon) ! longueur de rugosite
950	c
951	c Conditions aux limites
952	c
953	INTEGER julien
954	c
955	INTEGER lmt_pas
956	SAVE lmt_pas ! frequence de mise a jour
957	c$OMP THREADPRIVATE(lmt_pas)
958	REAL,allocatable,save :: pctsrf(:,:)
959	c$OMP THREADPRIVATE(pctsrf)
960	cIM
961	REAL pctsrf_new(klon,nbsrf) !pourcentage surfaces issus d'ORCHIDEE
962
963	cym REAL paire_ter(klon) !surfaces terre
964	REAL,allocatable,save :: paire_ter(:) !surfaces terre
965	c$OMP THREADPRIVATE(paire_ter)
966
967	cIM
968	cym SAVE pctsrf ! sous-fraction du sol
969	REAL,allocatable,save :: albsol(:)
970	c$OMP THREADPRIVATE(albsol)
971	cym SAVE albsol ! albedo du sol total
972	REAL,allocatable,save :: albsollw(:)
973	c$OMP THREADPRIVATE(albsollw)
974	cym SAVE albsollw ! albedo du sol total
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 (klon.eq.1) then
1481	print*,'WARNING !!!! omega=0'
1482	omega=0.
1483	igout=1
1484	write(lunout,*) 'DEBUT DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
1485	write(lunout,*)
1486	s 'nlon,nlev,nqmax,debut,lafin,rjourvrai,gmtime,pdtphys'
1487	write(lunout,*)
1488	s nlon,nlev,nqmax,debut,lafin,rjourvrai,gmtime,pdtphys
1489
1490	write(lunout,*) 'papers, play, phi, u, v, t, omega'
1491	do k=1,nlev
1492	write(lunout,*) paprs(igout,k),pplay(igout,k),pphi(igout,k),
1493	s u(igout,k),v(igout,k),t(igout,k),omega(igout,k)
1494	enddo
1495	write(lunout,*) 'ovap (g/kg), oliq (g/kg)'
1496	do k=1,nlev
1497	write(lunout,) qx(igout,k,1)1000,qx(igout,k,2)*1000.
1498	enddo
1499	endif
1500
1501	c======================================================================
1502
1503	cym => necessaire pour iflag_con != 2
1504	pmfd(:,:) = 0.
1505	pen_u(:,:) = 0.
1506	pen_d(:,:) = 0.
1507	pde_d(:,:) = 0.
1508	pde_u(:,:) = 0.
1509	aam=0.
1510	torsfc=0.
1511
1512	if (first) then
1513
1514	allocate( t_ancien(klon,klev), q_ancien(klon,klev))
1515	allocate( swdn0(klon,klevp1), swdn(klon,klevp1))
1516	allocate( swup0(klon,klevp1), swup(klon,klevp1))
1517	allocate( SWdn200clr(klon), SWdn200(klon))
1518	allocate( SWup200clr(klon), SWup200(klon))
1519	allocate( lwdn0(klon,klevp1), lwdn(klon,klevp1))
1520	allocate( lwup0(klon,klevp1), lwup(klon,klevp1))
1521	allocate( LWdn200clr(klon), LWdn200(klon))
1522	allocate( LWup200clr(klon), LWup200(klon))
1523	allocate( LWdnTOA(klon), LWdnTOAclr(klon))
1524	allocate( radsol(klon))
1525	allocate( rlat(klon))
1526	allocate( rlon(klon))
1527	allocate( ftsol(klon,nbsrf))
1528	allocate( deltat(klon))
1529	allocate( falbe(klon,nbsrf))
1530	allocate( falblw(klon,nbsrf))
1531	allocate( zmea(klon))
1532	allocate( zstd(klon))
1533	allocate( zsig(klon))
1534	allocate( zgam(klon))
1535	allocate( zthe(klon))
1536	allocate( zpic(klon))
1537	allocate( zval(klon))
1538	allocate( rugoro(klon))
1539	allocate( zuthe(klon),zvthe(klon))
1540	allocate( alb_neig(klon))
1541	allocate( ema_workcbmf(klon))
1542	cCR:nvelles variables convection/poches froides
1543	allocate( sigd(klon))
1544	allocate( cin(klon))
1545	allocate( ftd(klon,klev))
1546	allocate( fqd(klon,klev))
1547	allocate( ALE(klon))
1548	allocate( ALP(klon))
1549	allocate( Ale_bl(klon))
1550	allocate( Alp_bl(klon))
1551	allocate( lalim_conv(klon))
1552	allocate( wght_th(klon,klev))
1553	allocate( wake_deltat(klon,klev))
1554	allocate( wake_deltaq(klon,klev))
1555	allocate( wake_Cstar(klon))
1556	allocate( wake_s(klon))
1557	allocate( wake_fip(klon))
1558	allocate( dt_wake(klon,klev))
1559	allocate( dq_wake(klon,klev))
1560	cfinCR
1561	allocate( ema_cbmf(klon))
1562	allocate( ema_pcb(klon))
1563	allocate( ema_pct(klon))
1564	allocate( Ma(klon,klev) )
1565	allocate( qcondc(klon,klev))
1566	allocate( ema_work1(klon, klev), ema_work2(klon, klev))
1567	allocate( wd(klon) )
1568	allocate( pfrac_impa(klon,klev))
1569	allocate( pfrac_nucl(klon,klev))
1570	allocate( pfrac_1nucl(klon,klev))
1571	allocate( rain_fall(klon) )
1572	allocate( snow_fall(klon) )
1573	allocate( total_rain(klon), nday_rain(klon))
1574	allocate( pctsrf(klon,nbsrf))
1575	allocate( albsol(klon))
1576	allocate( albsollw(klon))
1577	allocate( wo(klon,klev))
1578	allocate( clwcon(klon,klev),rnebcon(klon,klev))
1579	allocate( heat(klon,klev) )
1580	allocate( heat0(klon,klev) )
1581	allocate( cool(klon,klev) )
1582	allocate( cool0(klon,klev) )
1583	allocate( topsw(klon), toplw(klon), solsw(klon), sollw(klon))
1584	allocate( sollwdown(klon) )
1585	allocate( sollwdownclr(klon) )
1586	allocate( toplwdown(klon) )
1587	allocate( toplwdownclr(klon) )
1588	allocate( topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon))
1589	allocate( albpla(klon))
1590	allocate( cape(klon) )
1591	allocate( pbase(klon) )
1592	allocate( bbase(klon) )
1593	allocate( ibas_con(klon), itop_con(klon))
1594	allocate( ratqs(klon,klev))
1595	allocate( sulfate_pi(klon, klev))
1596	allocate( paire_ter(klon))
1597	allocate(tsumSTD(klon,nlevSTD,nout))
1598	allocate(usumSTD(klon,nlevSTD,nout))
1599	allocate(vsumSTD(klon,nlevSTD,nout))
1600	allocate(wsumSTD(klon,nlevSTD,nout))
1601	allocate(phisumSTD(klon,nlevSTD,nout))
1602	allocate(qsumSTD(klon,nlevSTD,nout))
1603	allocate(rhsumSTD(klon,nlevSTD,nout))
1604	allocate(uvsumSTD(klon,nlevSTD,nout))
1605	allocate(vqsumSTD(klon,nlevSTD,nout))
1606	allocate(vTsumSTD(klon,nlevSTD,nout))
1607	allocate(wqsumSTD(klon,nlevSTD,nout))
1608	allocate( vphisumSTD(klon,nlevSTD,nout))
1609	allocate( wTsumSTD(klon,nlevSTD,nout))
1610	allocate( u2sumSTD(klon,nlevSTD,nout))
1611	allocate( v2sumSTD(klon,nlevSTD,nout))
1612	allocate( T2sumSTD(klon,nlevSTD,nout))
1613	allocate( seed_old(klon,napisccp))
1614	allocate( pct_ocean(klon,nbregdyn))
1615	allocate( rlonPOS(klon))
1616	allocate( newsst(klon))
1617	allocate( zqasc(klon,klev))
1618	allocate( rain_con(klon))
1619	allocate( u10m(klon,nbsrf), v10m(klon,nbsrf))
1620	allocate( topswad(klon), solswad(klon))
1621	allocate( topswai(klon), solswai(klon) )
1622	#ifdef INCA_AER
1623	#ifdef CPP_COUPLE
1624	allocate( topswad_inca(klon), solswad_inca(klon))
1625	allocate( topswad0_inca(klon), solswad0_inca(klon))
1626	allocate( topswai_inca(klon), solswai_inca(klon))
1627	allocate( topsw_inca(klon,9), solsw_inca(klon,9))
1628	allocate( topsw0_inca(klon,9), solsw0_inca(klon,9))
1629	allocate( tau_inca(klon,klev,9,2))
1630	allocate( piz_inca(klon,klev,9,2))
1631	allocate( cg_inca(klon,klev,9,2))
1632	allocate( ccm(klon,klev,2) )
1633	#endif
1634	#endif
1635	allocate( clwcon0(klon,klev),rnebcon0(klon,klev))
1636	allocate( clwcon0th(klon,klev),rnebcon0th(klon,klev))
1637	allocate( tau_ae(klon,klev,2), piz_ae(klon,klev,2))
1638	allocate( cg_ae(klon,klev,2))
1639	allocate( snow_con(klon))
1640	allocate( tnondef(klon,nlevSTD,nout))
1641	allocate( d_u_con(klon,klev),d_v_con(klon,klev))
1642
1643
1644	paire_ter(:)=0.
1645	clwcon(:,:)=0.
1646	rnebcon(:,:)=0.
1647	ratqs(:,:)=0.
1648	sollw(:)=0.
1649	ema_work1(:,:)=0.
1650	ema_work2(:,:)=0.
1651	cym Attention pbase pas initialise dans concvl !!!!
1652	pbase(:)=0
1653
1654	first=.false.
1655	endif
1656
1657
1658	modname = 'physiq'
1659	cIM
1660	IF (ip_ebil_phy.ge.1) THEN
1661	DO i=1,klon
1662	zero_v(i)=0.
1663	END DO
1664	END IF
1665	ok_sync=.TRUE.
1666	IF (nqmax .LT. 2) THEN
1667	abort_message = 'eaux vapeur et liquide sont indispensables'
1668	CALL abort_gcm (modname,abort_message,1)
1669	ENDIF
1670	IF (debut) THEN
1671	CALL suphel ! initialiser constantes et parametres phys.
1672	ENDIF
1673
1674	print*,'CONVERGENCE PHYSIQUE THERM 1 '
1675
1676
1677	c======================================================================
1678	xjour = rjourvrai
1679	c
1680	c Si c'est le debut, il faut initialiser plusieurs choses
1681	c ********
1682	c
1683	IF (debut) THEN
1684	C
1685	!rv
1686	cCRinitialisation de wght_th et lalim_conv pour la definition de la couche alimentation
1687	cde la convection a partir des caracteristiques du thermique
1688	wght_th(:,:)=1.
1689	lalim_conv(:)=1
1690	cRC
1691	u10m(:,:)=0.
1692	v10m(:,:)=0.
1693	piz_ae(:,:,:)=0.
1694	tau_ae(:,:,:)=0.
1695	cg_ae(:,:,:)=0.
1696	rain_con(:)=0.
1697	snow_con(:)=0.
1698	bl95_b0=0.
1699	bl95_b1=0.
1700	topswai(:)=0.
1701	topswad(:)=0.
1702	solswai(:)=0.
1703	solswad(:)=0.
1704	#ifdef INCA_AER
1705	#ifdef CPP_COUPLE
1706	tau_inca(:,:,:,:) = 0.
1707	piz_inca(:,:,:,:) = 0.
1708	cg_inca(:,:,:,:) = 0.
1709	ccm(:,:,:) = 0.
1710	topswai_inca(:) = 0.
1711	topswad_inca(:) = 0.
1712	topswad0_inca(:) = 0.
1713	topsw_inca(:,:) = 0.
1714	topsw0_inca(:,:) = 0.
1715	solswai_inca(:) = 0.
1716	solswad_inca(:) = 0.
1717	solswad0_inca(:) = 0.
1718	solsw_inca(:,:) = 0.
1719	solsw0_inca(:,:) = 0.
1720	#endif
1721	#endif
1722	!rv
1723	!ACo
1724	d_u_con(:,:) = 0.0
1725	d_v_con(:,:) = 0.0
1726	rnebcon0(:,:) = 0.0
1727	clwcon0(:,:) = 0.0
1728	rnebcon(:,:) = 0.0
1729	clwcon(:,:) = 0.0
1730	paire_ter(:) = 0.0
1731	c nhistoW(:,:,:,:) = 0.0
1732	c histoW(:,:,:,:) = 0.0
1733
1734	cIM
1735	IF (ip_ebil_phy.ge.1) d_h_vcol_phy=0.
1736	c
1737	c appel a la lecture du run.def physique
1738	c
1739	call conf_phys(ocean, ok_veget, ok_journe, ok_mensuel,
1740	. ok_instan, ok_hf, 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	. falbe, falblw, 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 calculs necessaires au calcul de l'albedo dans l'interface
2259	c
2260	CALL orbite(FLOAT(julien),zlongi,dist)
2261	IF (cycle_diurne) THEN
2262	zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
2263	CALL zenang(zlongi,gmtime,zdtime,rlat,rlon,rmu0,fract)
2264	ELSE
2265	rmu0 = -999.999
2266	ENDIF
2267
2268	if (mydebug) then
2269	call writefield_phy('u_seri',u_seri,llm)
2270	call writefield_phy('v_seri',v_seri,llm)
2271	call writefield_phy('t_seri',t_seri,llm)
2272	call writefield_phy('q_seri',q_seri,llm)
2273	endif
2274
2275	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2276	c Appel au pbl_surface : Planetary Boudary Layer et Surface
2277	c Cela implique tous les interactions des sous-surfaces et la partie diffusion
2278	c turbulent du couche limit.
2279	c
2280	c Certains varibales de sorties de pbl_surface sont utiliser que pour
2281	c ecriture des fihiers hist_XXXX.nc, ces sont :
2282	c qsol, zq2m, s_pblh, s_lcl,
2283	c s_capCL, s_oliqCL, s_cteiCL,s_pblT,
2284	c s_therm, s_trmb1, s_trmb2, s_trmb3,
2285	c zxrugs, zu10m, zv10m, fder,
2286	c zxqsurf, rh2m, zxfluxu, zxfluxv,
2287	c frugs, agesno, fsollw, fsolsw,
2288	c d_ts, fevap, fluxlat, t2m,
2289	c wfbils, wfbilo, fluxt, fluxu, fluxv,
2290	c
2291	c Certains ne sont pas utiliser du tout :
2292	c dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
2293	c
2294	CALL pbl_surface(
2295	e dtime, date0, itap, julien,
2296	e debut, lafin,
2297	e rlon, rlat, rugoro, rmu0,
2298	e rain_fall, snow_fall, solsw, sollw,
2299	e t_seri, q_seri, u_seri, v_seri,
2300	e pplay, paprs, pctsrf,
2301	+ ftsol, falbe, falblw, u10m, v10m,
2302	s sollwdown, cdragh, cdragm, yu1, yv1,
2303	s albsol, albsollw, sens, evap,
2304	s zxtsol, zxfluxlat, zt2m, qsat2m,
2305	s d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf,
2306	s ycoefh, pctsrf_new,
2307	d qsol, zq2m, s_pblh, s_lcl,
2308	d s_capCL, s_oliqCL, s_cteiCL,s_pblT,
2309	d s_therm, s_trmb1, s_trmb2, s_trmb3,
2310	d zxrugs, zu10m, zv10m, fder,
2311	d zxqsurf, rh2m, zxfluxu, zxfluxv,
2312	d frugs, agesno, fsollw, fsolsw,
2313	d d_ts, fevap, fluxlat, t2m,
2314	d wfbils, wfbilo, fluxt, fluxu, fluxv,
2315	- dsens, devap, zxsnow,
2316	- zxfluxt, zxfluxq, q2m, fluxq, pbl_tke )
2317	c
2318	c
2319	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2320
2321	pctsrf(:,:) = pctsrf_new(:,:)
2322
2323	DO k = 1, klev
2324	DO i = 1, klon
2325	t_seri(i,k) = t_seri(i,k) + d_t_vdf(i,k)
2326	q_seri(i,k) = q_seri(i,k) + d_q_vdf(i,k)
2327	u_seri(i,k) = u_seri(i,k) + d_u_vdf(i,k)
2328	v_seri(i,k) = v_seri(i,k) + d_v_vdf(i,k)
2329	ENDDO
2330	ENDDO
2331
2332	if (mydebug) then
2333	call writefield_phy('u_seri',u_seri,llm)
2334	call writefield_phy('v_seri',v_seri,llm)
2335	call writefield_phy('t_seri',t_seri,llm)
2336	call writefield_phy('q_seri',q_seri,llm)
2337	endif
2338
2339
2340	IF (ip_ebil_phy.ge.2) THEN
2341	ztit='after surface_main'
2342	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2343	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2344	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2345	call diagphy(airephy,ztit,ip_ebil_phy
2346	e , zero_v, zero_v, zero_v, zero_v, sens
2347	e , evap , zero_v, zero_v, ztsol
2348	e , d_h_vcol, d_qt, d_ec
2349	s , fs_bound, fq_bound )
2350	END IF
2351
2352	c
2353	c Appeler la convection (au choix)
2354	c
2355	DO k = 1, klev
2356	DO i = 1, klon
2357	conv_q(i,k) = d_q_dyn(i,k)
2358	. + d_q_vdf(i,k)/dtime
2359	conv_t(i,k) = d_t_dyn(i,k)
2360	. + d_t_vdf(i,k)/dtime
2361	ENDDO
2362	ENDDO
2363	IF (check) THEN
2364	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2365	WRITE(lunout,*) "avantcon=", za
2366	ENDIF
2367	zx_ajustq = .FALSE.
2368	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
2369	IF (zx_ajustq) THEN
2370	DO i = 1, klon
2371	z_avant(i) = 0.0
2372	ENDDO
2373	DO k = 1, klev
2374	DO i = 1, klon
2375	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k))
2376	. *(paprs(i,k)-paprs(i,k+1))/RG
2377	ENDDO
2378	ENDDO
2379	ENDIF
2380	IF (iflag_con.EQ.1) THEN
2381	stop'reactiver le call conlmd dans physiq.F'
2382	c CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
2383	c . d_t_con, d_q_con,
2384	c . rain_con, snow_con, ibas_con, itop_con)
2385	ELSE IF (iflag_con.EQ.2) THEN
2386	CALL conflx(dtime, paprs, pplay, t_seri, q_seri,
2387	e conv_t, conv_q, -evap, omega,
2388	s d_t_con, d_q_con, rain_con, snow_con,
2389	s pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
2390	s kcbot, kctop, kdtop, pmflxr, pmflxs)
2391	WHERE (rain_con < 0.) rain_con = 0.
2392	WHERE (snow_con < 0.) snow_con = 0.
2393	DO i = 1, klon
2394	ibas_con(i) = klev+1 - kcbot(i)
2395	itop_con(i) = klev+1 - kctop(i)
2396	ENDDO
2397	ELSE IF (iflag_con.GE.3) THEN
2398	c nb of tracers for the KE convection:
2399	c MAF la partie traceurs est faite dans phytrac
2400	c on met ntra=1 pour limiter les appels mais on peut
2401	c supprimer les calculs / ftra.
2402	ntra = 1
2403
2404	c=====================================================================================
2405	cajout pour la parametrisation des poches froides:
2406	ccalcul de t_wake et t_undi: si pas de poches froides, t_wake=t_undi=t_seri
2407	do k=1,klev
2408	do i=1,klon
2409	if (iflag_wake.eq.1) then
2410	t_wake(i,k) = t_seri(i,k)
2411	. +(1-wake_s(i))*wake_deltat(i,k)
2412	q_wake(i,k) = q_seri(i,k)
2413	. +(1-wake_s(i))*wake_deltaq(i,k)
2414	t_undi(i,k) = t_seri(i,k)
2415	. -wake_s(i)*wake_deltat(i,k)
2416	q_undi(i,k) = q_seri(i,k)
2417	. -wake_s(i)*wake_deltaq(i,k)
2418	else
2419	t_wake(i,k) = t_seri(i,k)
2420	q_wake(i,k) = q_seri(i,k)
2421	t_undi(i,k) = t_seri(i,k)
2422	q_undi(i,k) = q_seri(i,k)
2423	endif
2424	enddo
2425	enddo
2426
2427	cc-- Calcul de l'energie disponible ALE (J/kg) et de la puissance disponible ALP (W/m2)
2428	cc-- pour le soulevement des particules dans le modele convectif
2429	c
2430	do i = 1,klon
2431	ALE(i) = 0.
2432	ALP(i) = 0.
2433	enddo
2434	c
2435	ccalcul de ale_wake et alp_wake
2436	do i = 1,klon
2437	if (iflag_wake.eq.1) then
2438	ale_wake(i) = 0.5wake_cstar(i)*2
2439	alp_wake(i) = wake_fip(i)
2440	else
2441	ale_wake(i) = 0.
2442	alp_wake(i) = 0.
2443	endif
2444	enddo
2445	ccombinaison avec ale et alp de couche limite: constantes si pas de couplage, valeurs calculees
2446	cdans le thermique sinon
2447	if (iflag_coupl.eq.0) then
2448	print*,'ALE et ALP imposes'
2449	do i = 1,klon
2450	con ne couple que ale
2451	c ALE(i) = max(ale_wake(i),Ale_bl(i))
2452	ALE(i) = max(ale_wake(i),ale_bl_prescr)
2453	con ne couple que alp
2454	c ALP(i) = alp_wake(i) + Alp_bl(i)
2455	ALP(i) = alp_wake(i) + alp_bl_prescr
2456	enddo
2457	else
2458	print*,'ALE et ALP couples au thermique'
2459	do i = 1,klon
2460	ALE(i) = max(ale_wake(i),Ale_bl(i))
2461	ALP(i) = alp_wake(i) + Alp_bl(i)
2462	c write(20,*)'ALE',ALE(i),Ale_bl(i),ale_wake(i)
2463	c write(21,*)'ALP',ALP(i),Alp_bl(i),alp_wake(i)
2464	enddo
2465	endif
2466
2467	cfin calcul ale et alp
2468	c=================================================================================================
2469
2470
2471	c sb, oct02:
2472	c Schema de convection modularise et vectorise:
2473	c (driver commun aux versions 3 et 4)
2474	c
2475	IF (ok_cvl) THEN ! new driver for convectL
2476
2477	CALL concvl (iflag_con,iflag_clos,
2478	. dtime,paprs,pplay,t_undi,q_undi,
2479	. t_wake,q_wake,
2480	. u_seri,v_seri,tr_seri,nbtr,
2481	. ALE,ALP,
2482	. ema_work1,ema_work2,
2483	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2484	. rain_con, snow_con, ibas_con, itop_con, sigd,
2485	. upwd,dnwd,dnwd0,
2486	. Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl,
2487	. pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd,
2488	. pmflxr,pmflxs,da,phi,mp,
2489	. ftd,fqd,lalim_conv,wght_th)
2490
2491	cIM cf. FH
2492	clwcon0=qcondc
2493	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
2494
2495	ELSE ! ok_cvl
2496	c MAF conema3 ne contient pas les traceurs
2497	CALL conema3 (dtime,
2498	. paprs,pplay,t_seri,q_seri,
2499	. u_seri,v_seri,tr_seri,ntra,
2500	. ema_work1,ema_work2,
2501	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2502	. rain_con, snow_con, ibas_con, itop_con,
2503	. upwd,dnwd,dnwd0,bas,top,
2504	. Ma,cape,tvp,rflag,
2505	. pbase
2506	. ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr
2507	. ,clwcon0)
2508
2509	ENDIF ! ok_cvl
2510
2511	c
2512	c Correction precip
2513	rain_con = rain_con * cvl_corr
2514	snow_con = snow_con * cvl_corr
2515	c
2516
2517	IF (.NOT. ok_gust) THEN
2518	do i = 1, klon
2519	wd(i)=0.0
2520	enddo
2521	ENDIF
2522
2523	c =================================================================== c
2524	c Calcul des proprietes des nuages convectifs
2525	c
2526	DO k = 1, klev
2527	DO i = 1, klon
2528	zx_t = t_seri(i,k)
2529	IF (thermcep) THEN
2530	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2531	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2532	zx_qs = MIN(0.5,zx_qs)
2533	zcor = 1./(1.-retv*zx_qs)
2534	zx_qs = zx_qs*zcor
2535	ELSE
2536	IF (zx_t.LT.t_coup) THEN
2537	zx_qs = qsats(zx_t)/pplay(i,k)
2538	ELSE
2539	zx_qs = qsatl(zx_t)/pplay(i,k)
2540	ENDIF
2541	ENDIF
2542	zqsat(i,k)=zx_qs
2543	ENDDO
2544	ENDDO
2545
2546	c calcul des proprietes des nuages convectifs
2547	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
2548	call clouds_gno
2549	s (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
2550
2551	c =================================================================== c
2552
2553	DO i = 1, klon
2554	ema_pcb(i) = pbase(i)
2555	ENDDO
2556	DO i = 1, klon
2557
2558	! L'idicage de itop_con peut cacher un pb potentiel
2559	! FH sous la dictee de JYG, CR
2560	ema_pct(i) = paprs(i,itop_con(i)+1)
2561
2562	if (itop_con(i).gt.klev-3) then
2563	print*,'La convection monte trop haut '
2564	print*,'itop_con(,',i,',)=',itop_con(i)
2565	endif
2566	ENDDO
2567	DO i = 1, klon
2568	ema_cbmf(i) = ema_workcbmf(i)
2569	ENDDO
2570	ELSE
2571	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
2572	CALL abort
2573	ENDIF
2574
2575	c CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
2576	c . d_u_con, d_v_con)
2577
2578	DO k = 1, klev
2579	DO i = 1, klon
2580	t_seri(i,k) = t_seri(i,k) + d_t_con(i,k)
2581	q_seri(i,k) = q_seri(i,k) + d_q_con(i,k)
2582	u_seri(i,k) = u_seri(i,k) + d_u_con(i,k)
2583	v_seri(i,k) = v_seri(i,k) + d_v_con(i,k)
2584	ENDDO
2585	ENDDO
2586
2587	if (mydebug) then
2588	call writefield_phy('u_seri',u_seri,llm)
2589	call writefield_phy('v_seri',v_seri,llm)
2590	call writefield_phy('t_seri',t_seri,llm)
2591	call writefield_phy('q_seri',q_seri,llm)
2592	endif
2593
2594	cIM
2595	IF (ip_ebil_phy.ge.2) THEN
2596	ztit='after convect'
2597	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2598	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2599	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2600	call diagphy(airephy,ztit,ip_ebil_phy
2601	e , zero_v, zero_v, zero_v, zero_v, zero_v
2602	e , zero_v, rain_con, snow_con, ztsol
2603	e , d_h_vcol, d_qt, d_ec
2604	s , fs_bound, fq_bound )
2605	END IF
2606	C
2607	IF (check) THEN
2608	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2609	WRITE(lunout,*)"aprescon=", za
2610	zx_t = 0.0
2611	za = 0.0
2612	DO i = 1, klon
2613	za = za + airephy(i)/FLOAT(klon)
2614	zx_t = zx_t + (rain_con(i)+
2615	. snow_con(i))*airephy(i)/FLOAT(klon)
2616	ENDDO
2617	zx_t = zx_t/za*dtime
2618	WRITE(lunout,*)"Precip=", zx_t
2619	ENDIF
2620	IF (zx_ajustq) THEN
2621	DO i = 1, klon
2622	z_apres(i) = 0.0
2623	ENDDO
2624	DO k = 1, klev
2625	DO i = 1, klon
2626	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k))
2627	. *(paprs(i,k)-paprs(i,k+1))/RG
2628	ENDDO
2629	ENDDO
2630	DO i = 1, klon
2631	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime)
2632	. /z_apres(i)
2633	ENDDO
2634	DO k = 1, klev
2635	DO i = 1, klon
2636	IF (z_factor(i).GT.(1.0+1.0E-08) .OR.
2637	. z_factor(i).LT.(1.0-1.0E-08)) THEN
2638	q_seri(i,k) = q_seri(i,k) * z_factor(i)
2639	ENDIF
2640	ENDDO
2641	ENDDO
2642	ENDIF
2643	zx_ajustq=.FALSE.
2644
2645	c
2646	c=============================================================================
2647	cRR:Evolution de la poche froide: on ne fait pas de separation wake/env
2648	cpour la couche limite diffuse pour l instant
2649	c
2650	if (iflag_wake.eq.1) then
2651	DO k=1,klev
2652	DO i=1,klon
2653	dt_dwn(i,k) = ftd(i,k)
2654	wdt_PBL(i,k) = 0.
2655	dq_dwn(i,k) = fqd(i,k)
2656	wdq_PBL(i,k) = 0.
2657	M_dwn(i,k) = dnwd0(i,k)
2658	M_up(i,k) = upwd(i,k)
2659	dt_a(i,k) = d_t_con(i,k)/dtime - ftd(i,k)
2660	udt_PBL(i,k) = 0.
2661	dq_a(i,k) = d_q_con(i,k)/dtime - fqd(i,k)
2662	udq_PBL(i,k) = 0.
2663	ENDDO
2664	ENDDO
2665	c
2666	ccalcul caracteristiques de la poche froide
2667	call calWAKE (paprs,pplay,dtime
2668	: ,t_seri,q_seri,omega,ibas_con
2669	: ,dt_dwn,dq_dwn,M_dwn,M_up
2670	: ,dt_a,dq_a,sigd
2671	: ,wdt_PBL,wdq_PBL
2672	: ,udt_PBL,udq_PBL
2673	o ,wake_deltat,wake_deltaq,wake_dth
2674	o ,wake_h,wake_s,wake_dens
2675	o ,wake_pe,wake_fip,wake_gfl
2676	o ,dt_wake,dq_wake
2677	o ,wake_k, t_undi,q_undi
2678	o ,wake_omgbdth,wake_dp_omgb
2679	o ,wake_dtKE,wake_dqKE
2680	o ,wake_dtPBL,wake_dqPBL
2681	o ,wake_omg,wake_dp_deltomg
2682	o ,wake_spread,wake_Cstar,wake_d_deltat_gw
2683	o ,wake_ddeltat,wake_ddeltaq)
2684	c
2685	cIncrementation des tendances de la poche froide
2686	DO k = 1, klev
2687	DO i = 1, klon
2688	t_seri(i,k) = t_seri(i,k) + dt_wake(i,k)*dtime
2689	q_seri(i,k) = q_seri(i,k) + dq_wake(i,k)*dtime
2690	ENDDO
2691	ENDDO
2692
2693	endif
2694	c print*,'apres callwake iflag_cldcon=', iflag_cldcon
2695	c
2696	c===================================================================
2697	c Convection seche (thermiques ou ajustement)
2698	c===================================================================
2699	c
2700	call stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri
2701	s ,seuil_inversion,weak_inversion,dthmin)
2702
2703
2704
2705	d_t_ajsb(:,:)=0.
2706	d_q_ajsb(:,:)=0.
2707	d_t_ajs(:,:)=0.
2708	d_u_ajs(:,:)=0.
2709	d_v_ajs(:,:)=0.
2710	d_q_ajs(:,:)=0.
2711	clwcon0th(:,:)=0.
2712	fm_therm(:,:)=0.
2713	entr_therm(:,:)=0.
2714	c
2715	IF(prt_level>9)WRITE(lunout,*)
2716	. 'AVANT LA CONVECTION SECHE , iflag_thermals='
2717	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2718	if(iflag_thermals.lt.0) then
2719	c Rien
2720	c ====
2721	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
2722
2723
2724	else
2725
2726	c Thermiques
2727	c ==========
2728	IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals='
2729	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2730	print*,'JUSTE AVANT , iflag_thermals='
2731	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2732
2733
2734	if (iflag_thermals.gt.1) then
2735	call calltherm(pdtphys
2736	s ,pplay,paprs,pphi,weak_inversion
2737	s ,u_seri,v_seri,t_seri,q_seri,zqsat,debut
2738	s ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
2739	s ,fm_therm,entr_therm,zqasc,clwcon0th,lmax_th,ratqscth
2740	s ,ratqsdiff,zqsatth
2741	con rajoute ale et alp, et les caracteristiques de la couche alim
2742	s ,Ale_bl,Alp_bl,lalim_conv,wght_th)
2743	endif
2744
2745	! call calltherm(pdtphys
2746	! s ,pplay,paprs,pphi
2747	! s ,u_seri,v_seri,t_seri,q_seri
2748	! s ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
2749	! s ,fm_therm,entr_therm)
2750
2751	c Ajustement sec
2752	c ==============
2753
2754	! Dans le cas où on active les thermiques, on fait partir l'ajustement
2755	! a partir du sommet des thermiques.
2756	! Dans le cas contraire, on demarre au niveau 1.
2757
2758	if (iflag_thermals.ge.13.or.iflag_thermals.eq.0) then
2759
2760	if(iflag_thermals.eq.0) then
2761	IF(prt_level>9)WRITE(lunout,*)'ajsec'
2762	limbas(:)=1
2763	else
2764	limbas(:)=lmax_th(:)
2765	endif
2766
2767	! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
2768	! pour des test de convergence numerique.
2769	! Le nouveau ajsec est a priori mieux, meme pour le cas
2770	! iflag_thermals = 0 (l'ancienne version peut faire des tendances
2771	! non nulles numeriquement pour des mailles non concernees.
2772
2773	if (iflag_thermals.eq.0) then
2774	CALL ajsec_convV2(paprs, pplay, t_seri,q_seri
2775	s , d_t_ajsb, d_q_ajsb)
2776	else
2777	CALL ajsec(paprs, pplay, t_seri,q_seri,limbas
2778	s , d_t_ajsb, d_q_ajsb)
2779	endif
2780
2781	t_seri(:,:) = t_seri(:,:) + d_t_ajsb(:,:)
2782	q_seri(:,:) = q_seri(:,:) + d_q_ajsb(:,:)
2783	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
2784	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
2785
2786	endif
2787
2788	endif
2789	c
2790	c===================================================================
2791	cIM
2792	IF (ip_ebil_phy.ge.2) THEN
2793	ztit='after dry_adjust'
2794	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2795	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2796	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2797	END IF
2798
2799
2800	c-------------------------------------------------------------------------
2801	c Caclul des ratqs
2802	c-------------------------------------------------------------------------
2803
2804	c print*,'calcul des ratqs'
2805	c ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
2806	c ----------------
2807	c on ecrase le tableau ratqsc calcule par clouds_gno
2808	if (iflag_cldcon.eq.1) then
2809	do k=1,klev
2810	do i=1,klon
2811	if(ptconv(i,k)) then
2812	ratqsc(i,k)=ratqsbas
2813	s +fact_cldcon*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
2814	else
2815	ratqsc(i,k)=0.
2816	endif
2817	enddo
2818	enddo
2819	else if (iflag_cldcon.eq.4) then
2820	ptconvth(:,:)=.false.
2821	ratqsc(:,:)=0.
2822	print*,'avant clouds_gno thermique'
2823	call clouds_gno
2824	s (klon,klev,q_seri,zqsat,clwcon0th,ptconvth,ratqsc,rnebcon0th)
2825	print*,' CLOUDS_GNO OK'
2826	endif
2827
2828	c ratqs stables
2829	c -------------
2830
2831	if (iflag_ratqs.eq.0) then
2832
2833	! Le cas iflag_ratqs=0 correspond a la version IPCC 2005 du modele.
2834	do k=1,klev
2835	do i=1, klon
2836	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2837	s min((paprs(i,1)-pplay(i,k))/(paprs(i,1)-30000.),1.)
2838	enddo
2839	enddo
2840
2841	! Pour iflag_ratqs=1 ou 2, le ratqs est constant au dessus de
2842	! 300 hPa (ratqshaut), varie lineariement en fonction de la pression
2843	! entre 600 et 300 hPa et est soit constant (ratqsbas) pour iflag_ratqs=1
2844	! soit lineaire (entre 0 a la surface et ratqsbas) pour iflag_ratqs=2
2845	! Il s'agit de differents tests dans la phase de reglage du modele
2846	! avec thermiques.
2847
2848	else if (iflag_ratqs.eq.1) then
2849
2850	do k=1,klev
2851	do i=1, klon
2852	if (pplay(i,k).ge.60000.) then
2853	ratqss(i,k)=ratqsbas
2854	else if ((pplay(i,k).ge.30000.).and.
2855	s (pplay(i,k).lt.60000.)) then
2856	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2857	s (60000.-pplay(i,k))/(60000.-30000.)
2858	else
2859	ratqss(i,k)=ratqshaut
2860	endif
2861	enddo
2862	enddo
2863
2864	else
2865
2866	do k=1,klev
2867	do i=1, klon
2868	if (pplay(i,k).ge.60000.) then
2869	ratqss(i,k)=ratqsbas
2870	s *(paprs(i,1)-pplay(i,k))/(paprs(i,1)-60000.)
2871	else if ((pplay(i,k).ge.30000.).and.
2872	s (pplay(i,k).lt.60000.)) then
2873	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2874	s (60000.-pplay(i,k))/(60000.-30000.)
2875	else
2876	ratqss(i,k)=ratqshaut
2877	endif
2878	enddo
2879	enddo
2880	endif
2881
2882
2883
2884
2885	c ratqs final
2886	c -----------
2887
2888	if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2
2889	s .or.iflag_cldcon.eq.4) then
2890
2891	! On ajoute une constante au ratqsc*2 pour tenir compte de
2892	! fluctuations turbulentes de petite echelle
2893
2894	do k=1,klev
2895	do i=1,klon
2896	if ((fm_therm(i,k).gt.1.e-10)) then
2897	ratqsc(i,k)=sqrt(ratqsc(i,k)2+0.052)
2898	endif
2899	enddo
2900	enddo
2901
2902	! les ratqs sont une conbinaison de ratqss et ratqsc
2903	! 1800s, en dur pour le moment, est le temps de
2904	! relaxation des ratqs
2905
2906	facteur=exp(-pdtphys/1800.)
2907
2908	print*,'WARNING ratqs a revoir '
2909	ratqs(:,:)=ratqsc(:,:)(1.-facteur)+ratqs(:,:)facteur
2910	ratqs(:,:)=sqrt(ratqs(:,:)2+ratqss(:,:)2)
2911
2912	else
2913	! on ne prend que le ratqs stable pour fisrtilp
2914	ratqs(:,:)=ratqss(:,:)
2915	endif
2916
2917
2918	c
2919	c Appeler le processus de condensation a grande echelle
2920	c et le processus de precipitation
2921	c-------------------------------------------------------------------------
2922	CALL fisrtilp(dtime,paprs,pplay,
2923	. t_seri, q_seri,ptconv,ratqs,
2924	. d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq,
2925	. rain_lsc, snow_lsc,
2926	. pfrac_impa, pfrac_nucl, pfrac_1nucl,
2927	. frac_impa, frac_nucl,
2928	. prfl, psfl, rhcl)
2929
2930	WHERE (rain_lsc < 0) rain_lsc = 0.
2931	WHERE (snow_lsc < 0) snow_lsc = 0.
2932	DO k = 1, klev
2933	DO i = 1, klon
2934	t_seri(i,k) = t_seri(i,k) + d_t_lsc(i,k)
2935	q_seri(i,k) = q_seri(i,k) + d_q_lsc(i,k)
2936	ql_seri(i,k) = ql_seri(i,k) + d_ql_lsc(i,k)
2937	cldfra(i,k) = rneb(i,k)
2938	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
2939	ENDDO
2940	ENDDO
2941	IF (check) THEN
2942	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2943	WRITE(lunout,*)"apresilp=", za
2944	zx_t = 0.0
2945	za = 0.0
2946	DO i = 1, klon
2947	za = za + airephy(i)/FLOAT(klon)
2948	zx_t = zx_t + (rain_lsc(i)
2949	. + snow_lsc(i))*airephy(i)/FLOAT(klon)
2950	ENDDO
2951	zx_t = zx_t/za*dtime
2952	WRITE(lunout,*)"Precip=", zx_t
2953	ENDIF
2954	cIM
2955	IF (ip_ebil_phy.ge.2) THEN
2956	ztit='after fisrt'
2957	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2958	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2959	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2960	call diagphy(airephy,ztit,ip_ebil_phy
2961	e , zero_v, zero_v, zero_v, zero_v, zero_v
2962	e , zero_v, rain_lsc, snow_lsc, ztsol
2963	e , d_h_vcol, d_qt, d_ec
2964	s , fs_bound, fq_bound )
2965	END IF
2966
2967	if (mydebug) then
2968	call writefield_phy('u_seri',u_seri,llm)
2969	call writefield_phy('v_seri',v_seri,llm)
2970	call writefield_phy('t_seri',t_seri,llm)
2971	call writefield_phy('q_seri',q_seri,llm)
2972	endif
2973
2974	c
2975	c-------------------------------------------------------------------
2976	c PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
2977	c-------------------------------------------------------------------
2978
2979	c 1. NUAGES CONVECTIFS
2980	c
2981	cIM cf FH
2982	c IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
2983	IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
2984	snow_tiedtke=0.
2985	c print*,'avant calcul de la pseudo precip '
2986	c print*,'iflag_cldcon',iflag_cldcon
2987	if (iflag_cldcon.eq.-1) then
2988	rain_tiedtke=rain_con
2989	else
2990	c print*,'calcul de la pseudo precip '
2991	rain_tiedtke=0.
2992	c print*,'calcul de la pseudo precip 0'
2993	do k=1,klev
2994	do i=1,klon
2995	if (d_q_con(i,k).lt.0.) then
2996	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys
2997	s *(paprs(i,k)-paprs(i,k+1))/rg
2998	endif
2999	enddo
3000	enddo
3001	endif
3002	c
3003	c call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
3004	c
3005
3006	c Nuages diagnostiques pour Tiedtke
3007	CALL diagcld1(paprs,pplay,
3008	cIM cf FH . rain_con,snow_con,ibas_con,itop_con,
3009	. rain_tiedtke,snow_tiedtke,ibas_con,itop_con,
3010	. diafra,dialiq)
3011	DO k = 1, klev
3012	DO i = 1, klon
3013	IF (diafra(i,k).GT.cldfra(i,k)) THEN
3014	cldliq(i,k) = dialiq(i,k)
3015	cldfra(i,k) = diafra(i,k)
3016	ENDIF
3017	ENDDO
3018	ENDDO
3019
3020	ELSE IF (iflag_cldcon.ge.3) THEN
3021	c On prend pour les nuages convectifs le max du calcul de la
3022	c convection et du calcul du pas de temps precedent diminue d'un facteur
3023	c facttemps
3024	facteur = pdtphys *facttemps
3025	do k=1,klev
3026	do i=1,klon
3027	rnebcon(i,k)=rnebcon(i,k)*facteur
3028	if (rnebcon0(i,k)clwcon0(i,k).gt.rnebcon(i,k)clwcon(i,k))
3029	s then
3030	rnebcon(i,k)=rnebcon0(i,k)
3031	clwcon(i,k)=clwcon0(i,k)
3032	endif
3033	enddo
3034	enddo
3035
3036	c
3037	cjq - introduce the aerosol direct and first indirect radiative forcings
3038	cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
3039	IF (ok_ade.OR.ok_aie) THEN
3040	#if defined(CPP_COUPLE) && !defined(INCA_AER)
3041	! Get sulfate aerosol distribution
3042	CALL readsulfate(rjourvrai, debut, sulfate)
3043	CALL readsulfate_preind(rjourvrai, debut, sulfate_pi)
3044
3045	! Calculate aerosol optical properties (Olivier Boucher)
3046	CALL aeropt(pplay, paprs, t_seri, sulfate, rhcl,
3047	. tau_ae, piz_ae, cg_ae, aerindex)
3048	#endif
3049	#if !defined(CPP_COUPLE)
3050	! Get sulfate aerosol distribution
3051	CALL readsulfate(rjourvrai, debut, sulfate)
3052	CALL readsulfate_preind(rjourvrai, debut, sulfate_pi)
3053
3054	! Calculate aerosol optical properties (Olivier Boucher)
3055	CALL aeropt(pplay, paprs, t_seri, sulfate, rhcl,
3056	. tau_ae, piz_ae, cg_ae, aerindex)
3057	#endif
3058	cym
3059	ELSE
3060	tau_ae(:,:,:)=0.0
3061	piz_ae(:,:,:)=0.0
3062	cg_ae(:,:,:)=0.0
3063	cym
3064	ENDIF
3065
3066	c
3067	cIM calcul nuages par le simulateur ISCCP
3068	c
3069	#ifdef histISCCP
3070	IF (ok_isccp) THEN
3071	cIM appel simulateur toutes les NINT(freq_ISCCP/dtime) heures
3072	IF (MOD(itap,NINT(freq_ISCCP/dtime)).EQ.0) THEN
3073	#include "calcul_simulISCCP.h"
3074	ENDIF !(MOD(itap,NINT(freq_ISCCP/dtime))
3075	ENDIF !ok_isccp
3076	#endif
3077
3078	c On prend la somme des fractions nuageuses et des contenus en eau
3079	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
3080	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
3081
3082	ENDIF
3083
3084	c
3085	c 2. NUAGES STARTIFORMES
3086	c
3087	IF (ok_stratus) THEN
3088	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
3089	DO k = 1, klev
3090	DO i = 1, klon
3091	IF (diafra(i,k).GT.cldfra(i,k)) THEN
3092	cldliq(i,k) = dialiq(i,k)
3093	cldfra(i,k) = diafra(i,k)
3094	ENDIF
3095	ENDDO
3096	ENDDO
3097	ENDIF
3098	c
3099	c Precipitation totale
3100	c
3101	DO i = 1, klon
3102	rain_fall(i) = rain_con(i) + rain_lsc(i)
3103	snow_fall(i) = snow_con(i) + snow_lsc(i)
3104	ENDDO
3105	cIM
3106	IF (ip_ebil_phy.ge.2) THEN
3107	ztit="after diagcld"
3108	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
3109	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3110	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3111	END IF
3112	c
3113	c Calculer l'humidite relative pour diagnostique
3114	c
3115	DO k = 1, klev
3116	DO i = 1, klon
3117	zx_t = t_seri(i,k)
3118	IF (thermcep) THEN
3119	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
3120	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
3121	zx_qs = MIN(0.5,zx_qs)
3122	zcor = 1./(1.-retv*zx_qs)
3123	zx_qs = zx_qs*zcor
3124	ELSE
3125	IF (zx_t.LT.t_coup) THEN
3126	zx_qs = qsats(zx_t)/pplay(i,k)
3127	ELSE
3128	zx_qs = qsatl(zx_t)/pplay(i,k)
3129	ENDIF
3130	ENDIF
3131	zx_rh(i,k) = q_seri(i,k)/zx_qs
3132	zqsat(i,k)=zx_qs
3133	ENDDO
3134	ENDDO
3135
3136	cIM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
3137	c equivalente a 2m (tpote) pour diagnostique
3138	c
3139	DO i = 1, klon
3140	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
3141	IF (thermcep) THEN
3142	IF(zt2m(i).LT.RTT) then
3143	Lheat=RLSTT
3144	ELSE
3145	Lheat=RLVTT
3146	ENDIF
3147	ELSE
3148	IF (zt2m(i).LT.RTT) THEN
3149	Lheat=RLSTT
3150	ELSE
3151	Lheat=RLVTT
3152	ENDIF
3153	ENDIF
3154	tpote(i) = tpot(i)*
3155	. EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
3156	ENDDO
3157
3158
3159	#ifdef INCA
3160	call VTe(VTphysiq)
3161	call VTb(VTinca)
3162	calday = FLOAT(julien) + gmtime
3163
3164	#ifdef INCA_AER
3165	call AEROSOL_METEO_CALC(calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs
3166	& ,prfl,psfl,pctsrf,airephy,xjour,rlat,rlon,u10m,v10m)
3167	#endif
3168
3169	#ifdef INCAINFO
3170	WRITE(lunout,*)'Appel CHEMHOOK_BEGIN ...'
3171	#endif
3172
3173	zxsnow_dummy(:) = 0.0
3174
3175	CALL chemhook_begin (calday,
3176	#if defined(INCA) && !defined(INCA_CH4) && !defined(INCA_NMHC) && !defined(INCA_AER)
3177	$ julien,
3178	$ gmtime,
3179	#endif
3180	$ pctsrf(1,1),
3181	$ rlat,
3182	$ rlon,
3183	$ airephy,
3184	$ paprs,
3185	$ pplay,
3186	$ ycoefh,
3187	$ pphi,
3188	$ t_seri,
3189	$ u,
3190	$ v,
3191	$ wo,
3192	$ q_seri,
3193	$ zxtsol,
3194	$ zxsnow_dummy,
3195	$ solsw,
3196	$ albsol,
3197	$ rain_fall,
3198	$ snow_fall,
3199	$ itop_con,
3200	$ ibas_con,
3201	$ cldfra,
3202	$ iim,
3203	$ jjm,
3204	#ifdef INCA_AER
3205	$ tr_seri,
3206	$ ftsol,
3207	$ paprs,
3208	$ cdragh,
3209	$ cdragm,
3210	$ pctsrf,
3211	$ pdtphys,
3212	$ itap)
3213	#else
3214	$ tr_seri)
3215	#endif
3216
3217
3218	#ifdef INCAINFO
3219	WRITE(lunout,*)'OK.'
3220	#endif
3221	call VTe(VTinca)
3222	call VTb(VTphysiq)
3223	#endif
3224	c
3225	c Calculer les parametres optiques des nuages et quelques
3226	c parametres pour diagnostiques:
3227	c
3228	if (ok_newmicro) then
3229	CALL newmicro (paprs, pplay,ok_newmicro,
3230	. t_seri, cldliq, cldfra, cldtau, cldemi,
3231	. cldh, cldl, cldm, cldt, cldq,
3232	. flwp, fiwp, flwc, fiwc,
3233	e ok_aie,
3234	#if defined(CPP_COUPLE) && defined(INCA_AER)
3235	e ccm(:,:,1), ccm(:,:,2),
3236	#else
3237	e sulfate, sulfate_pi,
3238	#endif
3239	e bl95_b0, bl95_b1,
3240	s cldtaupi, re, fl)
3241	else
3242	CALL nuage (paprs, pplay,
3243	. t_seri, cldliq, cldfra, cldtau, cldemi,
3244	. cldh, cldl, cldm, cldt, cldq,
3245	e ok_aie,
3246	e sulfate, sulfate_pi,
3247	e bl95_b0, bl95_b1,
3248	s cldtaupi, re, fl)
3249
3250	endif
3251	c
3252	c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
3253	c
3254	IF (MOD(itaprad,radpas).EQ.0) THEN
3255
3256	DO i = 1, klon
3257	albsol(i) = falbe(i,is_oce) * pctsrf(i,is_oce)
3258	. + falbe(i,is_lic) * pctsrf(i,is_lic)
3259	. + falbe(i,is_ter) * pctsrf(i,is_ter)
3260	. + falbe(i,is_sic) * pctsrf(i,is_sic)
3261	albsollw(i) = falblw(i,is_oce) * pctsrf(i,is_oce)
3262	. + falblw(i,is_lic) * pctsrf(i,is_lic)
3263	. + falblw(i,is_ter) * pctsrf(i,is_ter)
3264	. + falblw(i,is_sic) * pctsrf(i,is_sic)
3265	ENDDO
3266
3267	if (mydebug) then
3268	call writefield_phy('u_seri',u_seri,llm)
3269	call writefield_phy('v_seri',v_seri,llm)
3270	call writefield_phy('t_seri',t_seri,llm)
3271	call writefield_phy('q_seri',q_seri,llm)
3272	endif
3273
3274
3275	#if defined(CPP_COUPLE) && defined(INCA_AER)
3276	CALL radlwsw_inca ! nouveau rayonnement (compatible Arpege-IFS)
3277	e (kdlon,kflev,dist, rmu0, fract,
3278	e paprs, pplay,zxtsol,albsol, albsollw, t_seri,q_seri,
3279	e wo,
3280	e cldfra, cldemi, cldtau,
3281	s heat,heat0,cool,cool0,radsol,albpla,
3282	s topsw,toplw,solsw,sollw,
3283	s sollwdown,
3284	s topsw0,toplw0,solsw0,sollw0,
3285	s lwdn0, lwdn, lwup0, lwup,
3286	s swdn0, swdn, swup0, swup,
3287	e ok_ade, ok_aie, ! new for aerosol radiative effects
3288	e tau_inca, piz_inca, cg_inca, ! ="=
3289	s topswad_inca, solswad_inca, ! ="=
3290	s topswad0_inca, solswad0_inca, ! ="=
3291	s topsw_inca, topsw0_inca,
3292	s solsw_inca, solsw0_inca,
3293	e cldtaupi, ! ="=
3294	s topswai_inca, solswai_inca) ! ="=
3295	#else
3296	CALL radlwsw ! nouveau rayonnement (compatible Arpege-IFS)
3297	e (dist, rmu0, fract,
3298	e paprs, pplay,zxtsol,albsol, albsollw, t_seri,q_seri,
3299	e wo,
3300	e cldfra, cldemi, cldtau,
3301	s heat,heat0,cool,cool0,radsol,albpla,
3302	s topsw,toplw,solsw,sollw,
3303	s sollwdown,
3304	s topsw0,toplw0,solsw0,sollw0,
3305	s lwdn0, lwdn, lwup0, lwup,
3306	s swdn0, swdn, swup0, swup,
3307	e ok_ade, ok_aie, ! new for aerosol radiative effects
3308	e tau_ae, piz_ae, cg_ae, ! ="=
3309	s topswad, solswad, ! ="=
3310	e cldtaupi, ! ="=
3311	s topswai, solswai) ! ="=
3312	#endif
3313	itaprad = 0
3314	ENDIF
3315	itaprad = itaprad + 1
3316
3317	if (iflag_radia.eq.0) then
3318	print *,'--------------------------------------------------'
3319	print *,'>>>> ATTENTION rayonnement desactive pour ce cas'
3320	print *,'>>>> heat et cool mis a zero '
3321	print *,'--------------------------------------------------'
3322	heat=0.
3323	cool=0.
3324	endif
3325
3326
3327	c
3328	c Ajouter la tendance des rayonnements (tous les pas)
3329	c
3330	DO k = 1, klev
3331	DO i = 1, klon
3332	t_seri(i,k) = t_seri(i,k)
3333	. + (heat(i,k)-cool(i,k)) * dtime/86400.
3334	ENDDO
3335	ENDDO
3336	c
3337	if (mydebug) then
3338	call writefield_phy('u_seri',u_seri,llm)
3339	call writefield_phy('v_seri',v_seri,llm)
3340	call writefield_phy('t_seri',t_seri,llm)
3341	call writefield_phy('q_seri',q_seri,llm)
3342	endif
3343
3344	cIM
3345	IF (ip_ebil_phy.ge.2) THEN
3346	ztit='after rad'
3347	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
3348	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3349	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3350	call diagphy(airephy,ztit,ip_ebil_phy
3351	e , topsw, toplw, solsw, sollw, zero_v
3352	e , zero_v, zero_v, zero_v, ztsol
3353	e , d_h_vcol, d_qt, d_ec
3354	s , fs_bound, fq_bound )
3355	END IF
3356	c
3357	c
3358	c Calculer l'hydrologie de la surface
3359	c
3360	c CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
3361	c . agesno, ftsol,fqsurf,fsnow, ruis)
3362	c
3363
3364	c
3365	c Calculer le bilan du sol et la derive de temperature (couplage)
3366	c
3367	DO i = 1, klon
3368	c bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
3369	c a la demande de JLD
3370	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
3371	ENDDO
3372	c
3373	cmoddeblott(jan95)
3374	c Appeler le programme de parametrisation de l'orographie
3375	c a l'echelle sous-maille:
3376	c
3377	IF (ok_orodr) THEN
3378	c
3379	c selection des points pour lesquels le shema est actif:
3380	igwd=0
3381	DO i=1,klon
3382	itest(i)=0
3383	c IF ((zstd(i).gt.10.0)) THEN
3384	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
3385	itest(i)=1
3386	igwd=igwd+1
3387	idx(igwd)=i
3388	ENDIF
3389	ENDDO
3390	c igwdim=MAX(1,igwd)
3391	c
3392	CALL drag_noro(klon,klev,dtime,paprs,pplay,
3393	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
3394	e igwd,idx,itest,
3395	e t_seri, u_seri, v_seri,
3396	cIM 141004 s zulow, zvlow, zustr, zvstr,
3397	s zulow, zvlow, zustrdr, zvstrdr,
3398	s d_t_oro, d_u_oro, d_v_oro)
3399	c
3400	c ajout des tendances
3401	DO k = 1, klev
3402	DO i = 1, klon
3403	t_seri(i,k) = t_seri(i,k) + d_t_oro(i,k)
3404	u_seri(i,k) = u_seri(i,k) + d_u_oro(i,k)
3405	v_seri(i,k) = v_seri(i,k) + d_v_oro(i,k)
3406	ENDDO
3407	ENDDO
3408	c
3409	ENDIF ! fin de test sur ok_orodr
3410	c
3411	if (mydebug) then
3412	call writefield_phy('u_seri',u_seri,llm)
3413	call writefield_phy('v_seri',v_seri,llm)
3414	call writefield_phy('t_seri',t_seri,llm)
3415	call writefield_phy('q_seri',q_seri,llm)
3416	endif
3417
3418	IF (ok_orolf) THEN
3419	c
3420	c selection des points pour lesquels le shema est actif:
3421	igwd=0
3422	DO i=1,klon
3423	itest(i)=0
3424	IF ((zpic(i)-zmea(i)).GT.100.) THEN
3425	itest(i)=1
3426	igwd=igwd+1
3427	idx(igwd)=i
3428	ENDIF
3429	ENDDO
3430	c igwdim=MAX(1,igwd)
3431	c
3432	CALL lift_noro(klon,klev,dtime,paprs,pplay,
3433	e rlat,zmea,zstd,zpic,
3434	e itest,
3435	e t_seri, u_seri, v_seri,
3436	s zulow, zvlow, zustrli, zvstrli,
3437	s d_t_lif, d_u_lif, d_v_lif)
3438	c
3439	c ajout des tendances
3440	DO k = 1, klev
3441	DO i = 1, klon
3442	t_seri(i,k) = t_seri(i,k) + d_t_lif(i,k)
3443	u_seri(i,k) = u_seri(i,k) + d_u_lif(i,k)
3444	v_seri(i,k) = v_seri(i,k) + d_v_lif(i,k)
3445	ENDDO
3446	ENDDO
3447	c
3448	ENDIF ! fin de test sur ok_orolf
3449	c
3450	cIM cf. FLott BEG
3451	C STRESS NECESSAIRES: TOUTE LA PHYSIQUE
3452
3453	if (mydebug) then
3454	call writefield_phy('u_seri',u_seri,llm)
3455	call writefield_phy('v_seri',v_seri,llm)
3456	call writefield_phy('t_seri',t_seri,llm)
3457	call writefield_phy('q_seri',q_seri,llm)
3458	endif
3459
3460	DO i = 1, klon
3461	zustrph(i)=0.
3462	zvstrph(i)=0.
3463	ENDDO
3464	DO k = 1, klev
3465	DO i = 1, klon
3466	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
3467	c (paprs(i,k)-paprs(i,k+1))/rg
3468	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
3469	c (paprs(i,k)-paprs(i,k+1))/rg
3470	ENDDO
3471	ENDDO
3472	c
3473	cIM calcul composantes axiales du moment angulaire et couple des montagnes
3474	c
3475	IF (is_sequential) THEN
3476
3477	CALL aaam_bud (27,klon,klev,rjourvrai,gmtime,
3478	C ra,rg,romega,
3479	C rlat,rlon,pphis,
3480	C zustrdr,zustrli,zustrph,
3481	C zvstrdr,zvstrli,zvstrph,
3482	C paprs,u,v,
3483	C aam, torsfc)
3484	ENDIF
3485	cIM cf. FLott END
3486	cIM
3487	IF (ip_ebil_phy.ge.2) THEN
3488	ztit='after orography'
3489	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
3490	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3491	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3492	END IF
3493	c
3494	c
3495	cAA
3496	cAA Installation de l'interface online-offline pour traceurs
3497	cAA
3498	c====================================================================
3499	c Calcul des tendances traceurs
3500	c====================================================================
3501	C
3502	call phytrac ( rnpb,
3503	I itap,
3504	I julien,
3505	I gmtime,
3506	I debut,
3507	I lafin,
3508	I nqmax-2,
3509	I nlon,
3510	I nlev,
3511	I dtime,
3512	I u,
3513	I v,
3514	I t,
3515	I paprs,
3516	I pplay,
3517	I pmfu,
3518	I pmfd,
3519	I pen_u,
3520	I pde_u,
3521	I pen_d,
3522	I pde_d,
3523	I ycoefh,
3524	I fm_therm,
3525	I entr_therm,
3526	I yu1,
3527	I yv1,
3528	I ftsol,
3529	I pctsrf,
3530	I rlat,
3531	I frac_impa,
3532	I frac_nucl,
3533	I rlon,
3534	I presnivs,
3535	I pphis,
3536	I pphi,
3537	I albsol,
3538	I qx(1,1,1),
3539	I rhcl,
3540	I cldfra,
3541	I rneb,
3542	I diafra,
3543	I cldliq,
3544	I itop_con,
3545	I ibas_con,
3546	I pmflxr,
3547	I pmflxs,
3548	I prfl,
3549	I psfl,
3550	I da,
3551	I phi,
3552	I mp,
3553	I upwd,
3554	I dnwd,
3555	#ifdef INCA
3556	I flxmass_w,
3557	#if defined(INCA_AER) && defined(CPP_COUPLE)
3558	I tau_inca,
3559	I piz_inca,
3560	I cg_inca,
3561	I ccm,
3562	I rfname,
3563	#endif
3564	#endif
3565	O tr_seri)
3566
3567	IF (offline) THEN
3568
3569	print*,'Attention on met a 0 les thermiques pour phystoke'
3570	call phystokenc (
3571	I nlon,nlev,pdtphys,rlon,rlat,
3572	I t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
3573	I fm_therm,entr_therm,
3574	I ycoefh,yu1,yv1,ftsol,pctsrf,
3575	I frac_impa, frac_nucl,
3576	I pphis,airephy,dtime,itap)
3577
3578
3579	ENDIF
3580
3581	c
3582	c Calculer le transport de l'eau et de l'energie (diagnostique)
3583	c
3584	CALL transp (paprs,zxtsol,
3585	e t_seri, q_seri, u_seri, v_seri, zphi,
3586	s ve, vq, ue, uq)
3587	c
3588	cIM global posePB BEG
3589	IF(1.EQ.0) THEN
3590	c
3591	CALL transp_lay (paprs,zxtsol,
3592	e t_seri, q_seri, u_seri, v_seri, zphi,
3593	s ve_lay, vq_lay, ue_lay, uq_lay)
3594	c
3595	ENDIF !(1.EQ.0) THEN
3596	cIM global posePB END
3597	c Accumuler les variables a stocker dans les fichiers histoire:
3598	c
3599	c+jld ec_conser
3600	DO k = 1, klev
3601	DO i = 1, klon
3602	ZRCPD = RCPD(1.0+RVTMP2q_seri(i,k))
3603	d_t_ec(i,k)=0.5/ZRCPD
3604	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
3605	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
3606	d_t_ec(i,k) = d_t_ec(i,k)/dtime
3607	END DO
3608	END DO
3609	c-jld ec_conser
3610	cIM
3611	IF (ip_ebil_phy.ge.1) THEN
3612	ztit='after physic'
3613	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
3614	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3615	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3616	C Comme les tendances de la physique sont ajoute dans la dynamique,
3617	C on devrait avoir que la variation d'entalpie par la dynamique
3618	C est egale a la variation de la physique au pas de temps precedent.
3619	C Donc la somme de ces 2 variations devrait etre nulle.
3620	call diagphy(airephy,ztit,ip_ebil_phy
3621	e , topsw, toplw, solsw, sollw, sens
3622	e , evap, rain_fall, snow_fall, ztsol
3623	e , d_h_vcol, d_qt, d_ec
3624	s , fs_bound, fq_bound )
3625	C
3626	d_h_vcol_phy=d_h_vcol
3627	C
3628	END IF
3629	C
3630	c=======================================================================
3631	c SORTIES
3632	c=======================================================================
3633
3634	cIM Interpolation sur les niveaux de pression du NMC
3635	c -------------------------------------------------
3636	c
3637	#include "calcul_STDlev.h"
3638	c
3639	c slp sea level pressure
3640	slp(:) = paprs(:,1)exp(pphis(:)/(RDt_seri(:,1)))
3641	c
3642	ccc prw = eau precipitable
3643	DO i = 1, klon
3644	prw(i) = 0.
3645	DO k = 1, klev
3646	prw(i) = prw(i) +
3647	. q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
3648	ENDDO
3649	ENDDO
3650	c
3651	cIM initialisation + calculs divers diag AMIP2
3652	c
3653	#include "calcul_divers.h"
3654	c
3655	#ifdef INCA
3656	call VTe(VTphysiq)
3657	call VTb(VTinca)
3658	#ifdef INCAINFO
3659	WRITE(lunout,*)'Appel CHEMHOOK_END ...'
3660	#endif
3661	CALL chemhook_end (calday,
3662	$ dtime,
3663	$ pplay,
3664	$ t_seri,
3665	$ tr_seri,
3666	$ nbtr,
3667	$ paprs,
3668	$ q_seri,
3669	$ annee_ref,
3670	$ day_ini,
3671	$ airephy,
3672	#ifdef INCA_AER
3673	$ xjour,
3674	$ pphi,
3675	$ pphis,
3676	$ zx_rh)
3677	#else
3678	$ xjour)
3679	#endif
3680	#ifdef INCAINFO
3681	WRITE(lunout,*)'OK.'
3682	#endif
3683	call VTe(VTinca)
3684	call VTb(VTphysiq)
3685	#endif
3686
3687	c=============================================================
3688	c
3689	c Convertir les incrementations en tendances
3690	c
3691	if (mydebug) then
3692	call writefield_phy('u_seri',u_seri,llm)
3693	call writefield_phy('v_seri',v_seri,llm)
3694	call writefield_phy('t_seri',t_seri,llm)
3695	call writefield_phy('q_seri',q_seri,llm)
3696	endif
3697
3698	DO k = 1, klev
3699	DO i = 1, klon
3700	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
3701	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
3702	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
3703	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
3704	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
3705	ENDDO
3706	ENDDO
3707	c
3708	IF (nqmax.GE.3) THEN
3709	DO iq = 3, nqmax
3710	DO k = 1, klev
3711	DO i = 1, klon
3712	d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
3713	ENDDO
3714	ENDDO
3715	ENDDO
3716	ENDIF
3717	c
3718	cIM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
3719	cIM global posePB#include "write_bilKP_ins.h"
3720	cIM global posePB#include "write_bilKP_ave.h"
3721	c
3722	c Sauvegarder les valeurs de t et q a la fin de la physique:
3723	c
3724	DO k = 1, klev
3725	DO i = 1, klon
3726	t_ancien(i,k) = t_seri(i,k)
3727	q_ancien(i,k) = q_seri(i,k)
3728	ENDDO
3729	ENDDO
3730	c
3731	!==========================================================================
3732	! Sorties des tendances pour un point particulier
3733	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
3734	! pour le debug
3735	! La valeur de igout est attribuee plus haut dans le programme
3736	!==========================================================================
3737
3738	if (klon.eq.1) then
3739	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
3740	write(lunout,*)
3741	s 'nlon,nlev,nqmax,debut,lafin,rjourvrai,gmtime,pdtphys'
3742	write(lunout,*)
3743	s nlon,nlev,nqmax,debut,lafin,rjourvrai,gmtime,pdtphys
3744	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
3745	do k=1,nlev
3746	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k),
3747	s d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k),
3748	s d_t_eva(igout,k)
3749	enddo
3750	write(lunout,*) 'cool,heat'
3751	do k=1,nlev
3752	write(lunout,*) cool(igout,k),heat(igout,k)
3753	enddo
3754
3755	write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
3756	do k=1,nlev
3757	write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k),
3758	s d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
3759	enddo
3760
3761	write(lunout,*) 'd_ps ',d_ps(igout)
3762	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
3763	do k=1,nlev
3764	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k),
3765	s d_qx(igout,k,1),d_qx(igout,k,2)
3766	enddo
3767	endif
3768
3769	!==========================================================================
3770
3771	c============================================================
3772	c Calcul de la temperature potentielle
3773	c============================================================
3774	DO k = 1, klev
3775	DO i = 1, klon
3776	theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
3777	ENDDO
3778	ENDDO
3779	c
3780
3781	c 22.03.04 BEG
3782	c=============================================================
3783	c Ecriture des sorties
3784	c=============================================================
3785	#ifdef CPP_IOIPSL
3786
3787	c Recupere des varibles calcule dans differents modules
3788	c pour ecriture dans histxxx.nc
3789
3790	! Get some variables from module mod_fonte_neige
3791	CALL fonte_neige_get_vars(pctsrf,
3792	. zxfqcalving, zxfqfonte, zxffonte)
3793
3794	IF (ocean == 'slab') THEN
3795	! Get some variables from module oceanslab
3796	CALL ocean_slab_get_vars(tslab, seaice, fluxo, fluxg)
3797	ELSEIF (ocean == 'couple') THEN
3798	! Get some variables from module oceancpl
3799	CALL ocean_cpl_get_vars(fluxo, fluxg)
3800	ELSE
3801	! Get some variables from module oceanforced
3802	CALL ocean_forced_get_vars(fluxo, fluxg)
3803	ENDIF
3804
3805	#ifdef histhf
3806	#include "write_histhf.h"
3807	#endif
3808
3809	#ifdef histday
3810	#include "write_histday.h"
3811	#endif
3812
3813	#ifdef histmth
3814	#include "write_histmth.h"
3815	#endif
3816
3817	#ifdef histins
3818	#include "write_histins.h"
3819	#endif
3820
3821	#ifdef histISCCP
3822	#include "write_histISCCP.h"
3823	#endif
3824
3825	#ifdef histmthNMC
3826	#include "write_histmthNMC.h"
3827	#endif
3828
3829	#include "write_histday_seri.h"
3830
3831	#include "write_paramLMDZ_phy.h"
3832
3833	#endif
3834
3835	c 22.03.04 END
3836	c
3837	c====================================================================
3838	c Si c'est la fin, il faut conserver l'etat de redemarrage
3839	c====================================================================
3840	c
3841
3842
3843	IF (lafin) THEN
3844	itau_phy = itau_phy + itap
3845	CALL phyredem ("restartphy.nc",dtime,radpas,ocean,
3846	. rlat, rlon, pctsrf, ftsol,
3847	. falbe,falblw, rain_fall,
3848	. snow_fall,
3849	. solsw, sollw,
3850	. radsol,
3851	. zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,
3852	. t_ancien, q_ancien, rnebcon, ratqs, clwcon,
3853	. pbl_tke)
3854	ENDIF
3855
3856
3857	RETURN
3858	END
3859	FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
3860	IMPLICIT none
3861	c
3862	c Calculer et imprimer l'eau totale. A utiliser pour verifier
3863	c la conservation de l'eau
3864	c
3865	#include "YOMCST.h"
3866	INTEGER klon,klev
3867	REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
3868	REAL aire(klon)
3869	REAL qtotal, zx, qcheck
3870	INTEGER i, k
3871	c
3872	zx = 0.0
3873	DO i = 1, klon
3874	zx = zx + aire(i)
3875	ENDDO
3876	qtotal = 0.0
3877	DO k = 1, klev
3878	DO i = 1, klon
3879	qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i)
3880	. *(paprs(i,k)-paprs(i,k+1))/RG
3881	ENDDO
3882	ENDDO
3883	c
3884	qcheck = qtotal/zx
3885	c
3886	RETURN
3887	END
3888	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
3889	IMPLICIT none
3890	c
3891	c Tranformer une variable de la grille physique a
3892	c la grille d'ecriture
3893	c
3894	INTEGER nfield,nlon,iim,jjmp1, jjm
3895	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
3896	c
3897	INTEGER i, n, ig
3898	c
3899	jjm = jjmp1 - 1
3900	DO n = 1, nfield
3901	DO i=1,iim
3902	ecrit(i,n) = fi(1,n)
3903	ecrit(i+jjm*iim,n) = fi(nlon,n)
3904	ENDDO
3905	DO ig = 1, nlon - 2
3906	ecrit(iim+ig,n) = fi(1+ig,n)
3907	ENDDO
3908	ENDDO
3909	RETURN
3910	END

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

Download in other formats:

Original Format