Context Navigation

physiq.F @ 878

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

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

Download in other formats:

Original Format