Context Navigation

physiq.F @ 852

Last change on this file since 852 was 852, checked in by Laurent Fairhead, 17 years ago
Mise a jour de la physique avec thermiques avec la version de FH d'aout 2007 LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 109.3 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	REAL,SAVE,ALLOCATABLE :: topswad(:), solswad(:) ! Aerosol direct effect.
1330	c$OMP THREADPRIVATE(topswad,solswad)
1331	! ok_ade=T -ADE=topswad-topsw
1332
1333	REAL,SAVE,ALLOCATABLE :: topswai(:), solswai(:) ! Aerosol indirect effect.
1334	c$OMP THREADPRIVATE(topswai,solswai)
1335	! ok_aie=T ->
1336	! ok_ade=T -AIE=topswai-topswad
1337	! ok_ade=F -AIE=topswai-topsw
1338
1339	REAL aerindex(klon) ! POLDER aerosol index
1340
1341	! Parameters
1342	LOGICAL ok_ade, ok_aie ! Apply aerosol (in)direct effects or not
1343	REAL bl95_b0, bl95_b1 ! Parameter in Boucher and Lohmann (1995)
1344	cym
1345	SAVE ok_ade, ok_aie, bl95_b0, bl95_b1
1346	c$OMP THREADPRIVATE(ok_ade, ok_aie, bl95_b0, bl95_b1)
1347
1348	c
1349	c Declaration des constantes et des fonctions thermodynamiques
1350	c
1351	LOGICAL,SAVE :: first=.true.
1352	c$OMP THREADPRIVATE(first)
1353	#include "YOMCST.h"
1354	#include "YOETHF.h"
1355	#include "FCTTRE.h"
1356	cIM 100106 BEG : pouvoir sortir les ctes de la physique
1357	#include "conema3.h"
1358	#include "fisrtilp.h"
1359	#include "nuage.h"
1360	#include "compbl.h"
1361	cIM 100106 END : pouvoir sortir les ctes de la physique
1362	c
1363	c======================================================================
1364
1365	cym => necessaire pour iflag_con != 2
1366	pmfd(:,:) = 0.
1367	pen_u(:,:) = 0.
1368	pen_d(:,:) = 0.
1369	pde_d(:,:) = 0.
1370	pde_u(:,:) = 0.
1371	aam=0.
1372	torsfc=0.
1373
1374	if (first) then
1375
1376	allocate( t_ancien(klon,klev), q_ancien(klon,klev))
1377	allocate( swdn0(klon,klevp1), swdn(klon,klevp1))
1378	allocate( swup0(klon,klevp1), swup(klon,klevp1))
1379	allocate( SWdn200clr(klon), SWdn200(klon))
1380	allocate( SWup200clr(klon), SWup200(klon))
1381	allocate( lwdn0(klon,klevp1), lwdn(klon,klevp1))
1382	allocate( lwup0(klon,klevp1), lwup(klon,klevp1))
1383	allocate( LWdn200clr(klon), LWdn200(klon))
1384	allocate( LWup200clr(klon), LWup200(klon))
1385	allocate( LWdnTOA(klon), LWdnTOAclr(klon))
1386	allocate( radsol(klon))
1387	allocate( rlat(klon))
1388	allocate( rlon(klon))
1389	allocate( ftsol(klon,nbsrf))
1390	allocate( deltat(klon))
1391	allocate( falbe(klon,nbsrf))
1392	allocate( falblw(klon,nbsrf))
1393	allocate( zmea(klon))
1394	allocate( zstd(klon))
1395	allocate( zsig(klon))
1396	allocate( zgam(klon))
1397	allocate( zthe(klon))
1398	allocate( zpic(klon))
1399	allocate( zval(klon))
1400	allocate( rugoro(klon))
1401	allocate( zuthe(klon),zvthe(klon))
1402	allocate( alb_neig(klon))
1403	allocate( ema_workcbmf(klon))
1404	allocate( ema_cbmf(klon))
1405	allocate( ema_pcb(klon))
1406	allocate( ema_pct(klon))
1407	allocate( Ma(klon,klev) )
1408	allocate( qcondc(klon,klev))
1409	allocate( ema_work1(klon, klev), ema_work2(klon, klev))
1410	allocate( wd(klon) )
1411	allocate( pfrac_impa(klon,klev))
1412	allocate( pfrac_nucl(klon,klev))
1413	allocate( pfrac_1nucl(klon,klev))
1414	allocate( rain_fall(klon) )
1415	allocate( snow_fall(klon) )
1416	allocate( total_rain(klon), nday_rain(klon))
1417	allocate( pctsrf(klon,nbsrf))
1418	allocate( albsol(klon))
1419	allocate( albsollw(klon))
1420	allocate( wo(klon,klev))
1421	allocate( clwcon(klon,klev),rnebcon(klon,klev))
1422	allocate( heat(klon,klev) )
1423	allocate( heat0(klon,klev) )
1424	allocate( cool(klon,klev) )
1425	allocate( cool0(klon,klev) )
1426	allocate( topsw(klon), toplw(klon), solsw(klon), sollw(klon))
1427	allocate( sollwdown(klon) )
1428	allocate( sollwdownclr(klon) )
1429	allocate( toplwdown(klon) )
1430	allocate( toplwdownclr(klon) )
1431	allocate( topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon))
1432	allocate( albpla(klon))
1433	allocate( cape(klon) )
1434	allocate( pbase(klon) )
1435	allocate( bbase(klon) )
1436	allocate( ibas_con(klon), itop_con(klon))
1437	allocate( ratqs(klon,klev))
1438	allocate( sulfate_pi(klon, klev))
1439	allocate( paire_ter(klon))
1440	allocate(tsumSTD(klon,nlevSTD,nout))
1441	allocate(usumSTD(klon,nlevSTD,nout))
1442	allocate(vsumSTD(klon,nlevSTD,nout))
1443	allocate(wsumSTD(klon,nlevSTD,nout))
1444	allocate(phisumSTD(klon,nlevSTD,nout))
1445	allocate(qsumSTD(klon,nlevSTD,nout))
1446	allocate(rhsumSTD(klon,nlevSTD,nout))
1447	allocate(uvsumSTD(klon,nlevSTD,nout))
1448	allocate(vqsumSTD(klon,nlevSTD,nout))
1449	allocate(vTsumSTD(klon,nlevSTD,nout))
1450	allocate(wqsumSTD(klon,nlevSTD,nout))
1451	allocate( vphisumSTD(klon,nlevSTD,nout))
1452	allocate( wTsumSTD(klon,nlevSTD,nout))
1453	allocate( u2sumSTD(klon,nlevSTD,nout))
1454	allocate( v2sumSTD(klon,nlevSTD,nout))
1455	allocate( T2sumSTD(klon,nlevSTD,nout))
1456	allocate( seed_old(klon,napisccp))
1457	allocate( pct_ocean(klon,nbregdyn))
1458	allocate( rlonPOS(klon))
1459	allocate( newsst(klon))
1460	allocate( zqasc(klon,klev))
1461	allocate( rain_con(klon))
1462	allocate( u10m(klon,nbsrf), v10m(klon,nbsrf))
1463	allocate( topswad(klon), solswad(klon))
1464	allocate( topswai(klon), solswai(klon) )
1465	allocate( clwcon0(klon,klev),rnebcon0(klon,klev))
1466	allocate( clwcon0th(klon,klev),rnebcon0th(klon,klev))
1467	allocate( tau_ae(klon,klev,2), piz_ae(klon,klev,2))
1468	allocate( cg_ae(klon,klev,2))
1469	allocate( snow_con(klon))
1470	allocate( tnondef(klon,nlevSTD,nout))
1471	allocate( d_u_con(klon,klev),d_v_con(klon,klev))
1472
1473
1474	paire_ter(:)=0.
1475	clwcon(:,:)=0.
1476	rnebcon(:,:)=0.
1477	ratqs(:,:)=0.
1478	sollw(:)=0.
1479	ema_work1(:,:)=0.
1480	ema_work2(:,:)=0.
1481	cym Attention pbase pas initialise dans concvl !!!!
1482	pbase(:)=0
1483
1484	first=.false.
1485	endif
1486
1487
1488	modname = 'physiq'
1489	cIM
1490	IF (ip_ebil_phy.ge.1) THEN
1491	DO i=1,klon
1492	zero_v(i)=0.
1493	END DO
1494	END IF
1495	ok_sync=.TRUE.
1496	IF (nqmax .LT. 2) THEN
1497	abort_message = 'eaux vapeur et liquide sont indispensables'
1498	CALL abort_gcm (modname,abort_message,1)
1499	ENDIF
1500	IF (debut) THEN
1501	CALL suphec ! initialiser constantes et parametres phys.
1502	ENDIF
1503
1504	print*,'CONVERGENCE PHYSIQUE THERM 1 '
1505
1506
1507	c======================================================================
1508	xjour = rjourvrai
1509	c
1510	c Si c'est le debut, il faut initialiser plusieurs choses
1511	c ********
1512	c
1513	IF (debut) THEN
1514	C
1515	!rv
1516	u10m(:,:)=0.
1517	v10m(:,:)=0.
1518	piz_ae(:,:,:)=0.
1519	tau_ae(:,:,:)=0.
1520	cg_ae(:,:,:)=0.
1521	rain_con(:)=0.
1522	snow_con(:)=0.
1523	bl95_b0=0.
1524	bl95_b1=0.
1525	topswai(:)=0.
1526	topswad(:)=0.
1527	solswai(:)=0.
1528	solswad(:)=0.
1529	!rv
1530	! anne
1531	d_u_con(:,:) = 0.0
1532	d_v_con(:,:) = 0.0
1533	rnebcon0(:,:) = 0.0
1534	clwcon0(:,:) = 0.0
1535	rnebcon(:,:) = 0.0
1536	clwcon(:,:) = 0.0
1537	paire_ter(:) = 0.0
1538	c nhistoW(:,:,:,:) = 0.0
1539	c histoW(:,:,:,:) = 0.0
1540	! fin anne
1541
1542	cIM
1543	IF (ip_ebil_phy.ge.1) d_h_vcol_phy=0.
1544	c
1545	c appel a la lecture du run.def physique
1546	c
1547	call conf_phys(ocean, ok_veget, ok_journe, ok_mensuel,
1548	. ok_instan, ok_hf, seuil_inversion,
1549	. fact_cldcon, facttemps,ok_newmicro,
1550	. iflag_cldcon,iflag_ratqs,ratqsbas,ratqshaut,
1551	. ok_ade, ok_aie,
1552	. bl95_b0, bl95_b1,
1553	. iflag_thermals,nsplit_thermals)
1554
1555	c
1556	c
1557	c Initialiser les compteurs:
1558	c
1559	itap = 0
1560	itaprad = 0
1561
1562	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1563	!! Un petit travail à faire ici.
1564	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1565
1566	if (iflag_pbl>1) then
1567	PRINT*, "Using method MELLOR&YAMADA"
1568	endif
1569	! NB! pbl_tke could/should be read and written from (re)startphy.nc
1570	ALLOCATE(pbl_tke(klon,klev+1,nbsrf))
1571	pbl_tke(:,:,:) = 1.e-8
1572
1573
1574	CALL phyetat0 ("startphy.nc",dtime,co2_ppm_etat0,solaire_etat0,
1575	. rlat,rlon,pctsrf, ftsol,
1576	. ocean, ok_veget,
1577	. falbe, falblw, rain_fall,snow_fall,
1578	. solsw, sollw,
1579	. radsol,clesphy0,
1580	. zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,tabcntr0,
1581	. t_ancien, q_ancien, ancien_ok, rnebcon, ratqs,clwcon,
1582	. pbl_tke)
1583
1584
1585
1586
1587	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1588
1589
1590	DO i=1,klon
1591	IF ( abs( pctsrf(i, is_ter) + pctsrf(i, is_lic) +
1592	$ pctsrf(i, is_oce) + pctsrf(i, is_sic) - 1.) .GT. EPSFRA)
1593	$ THEN
1594	WRITE(,)
1595	$ 'physiq apres lecture de restart: pb sous surface au point ',
1596	$ i, pctsrf(i, 1 : nbsrf)
1597	ENDIF
1598	ENDDO
1599
1600	radpas = NINT( 86400./dtime/nbapp_rad)
1601	c
1602	C on remet le calendrier a zero
1603	c
1604	IF (raz_date .eq. 1) THEN
1605	itau_phy = 0
1606	ENDIF
1607
1608	cIM cf. AM 081204 BEG
1609	PRINT*,'cycle_diurne3 =',cycle_diurne
1610	cIM cf. AM 081204 END
1611	c
1612	CALL printflag( tabcntr0,radpas,ok_journe,
1613	, ok_instan, ok_region )
1614	c
1615	IF (ABS(dtime-pdtphys).GT.0.001) THEN
1616	WRITE(lunout,*) 'Pas physique n est pas correct',dtime,
1617	. pdtphys
1618	abort_message='Pas physique n est pas correct '
1619	! call abort_gcm(modname,abort_message,1)
1620	dtime=pdtphys
1621	ENDIF
1622	IF (nlon .NE. klon) THEN
1623	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,
1624	. klon
1625	abort_message='nlon et klon ne sont pas coherents'
1626	call abort_gcm(modname,abort_message,1)
1627	ENDIF
1628	IF (nlev .NE. klev) THEN
1629	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,
1630	. klev
1631	abort_message='nlev et klev ne sont pas coherents'
1632	call abort_gcm(modname,abort_message,1)
1633	ENDIF
1634	c
1635	IF (dtime*FLOAT(radpas).GT.21600..AND.cycle_diurne) THEN
1636	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
1637	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
1638	abort_message='Nbre d appels au rayonnement insuffisant'
1639	call abort_gcm(modname,abort_message,1)
1640	ENDIF
1641	WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
1642	WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=",
1643	. ok_cvl
1644	c
1645	cKE43
1646	c Initialisation pour la convection de K.E. (sb):
1647	IF (iflag_con.GE.3) THEN
1648
1649	WRITE(lunout,)"** Convection de Kerry Emanuel 4.3 "
1650	WRITE(lunout,*)
1651	. "On va utiliser le melange convectif des traceurs qui"
1652	WRITE(lunout,*)"est calcule dans convect4.3"
1653	WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
1654
1655	DO i = 1, klon
1656	ema_cbmf(i) = 0.
1657	ema_pcb(i) = 0.
1658	ema_pct(i) = 0.
1659	ema_workcbmf(i) = 0.
1660	ENDDO
1661	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
1662	DO i = 1, klon
1663	ibas_con(i) = 1
1664	itop_con(i) = 1
1665	ENDDO
1666	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
1667
1668	ENDIF
1669
1670	rugoro=0.
1671	c34EK
1672	IF (ok_orodr) THEN
1673
1674	rugoro=0.
1675
1676	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1677	! FH sans doute a enlever de finitivement ou, si on le garde, l'activer
1678	! justement quand ok_orodr = false.
1679	! ce rugoro est utilise par la couche limite et fait double emploi
1680	! avec les paramétrisations spécifiques de Francois Lott.
1681	! DO i=1,klon
1682	! rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1683	! ENDDO
1684	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1685
1686	CALL SUGWD(klon,klev,paprs,pplay)
1687	DO i=1,klon
1688	zuthe(i)=0.
1689	zvthe(i)=0.
1690	if(zstd(i).gt.10.)then
1691	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
1692	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
1693	endif
1694	ENDDO
1695	ENDIF
1696	c
1697	c
1698	lmt_pas = NINT(86400./dtime * 1.0) ! tous les jours
1699	WRITE(lunout,*)'La frequence de lecture surface est de ',
1700	. lmt_pas
1701	c
1702	cIM200505 ecrit_mth = NINT(86400./dtime *ecritphy) ! tous les ecritphy jours
1703	c IF (ok_mensuel) THEN
1704	c WRITE(lunout,*)'La frequence de sortie mensuelle est de ',
1705	c . ecrit_mth
1706	c ENDIF
1707	c ecrit_day = NINT(86400./dtime *1.0) ! tous les jours
1708	c IF (ok_journe) THEN
1709	c WRITE(lunout,*)'La frequence de sortie journaliere est de ',
1710	c . ecrit_day
1711	c ENDIF
1712	cIM 130904 BEG
1713	cIM 080205 ecrit_hf = 86400./dtime *0.25 ! toutes les 6h
1714	cIM 170305
1715	c ecrit_hf = 86400./dtime/12. ! toutes les 2h
1716	cIM 230305
1717	cIM200505 ecrit_hf = 86400./dtime *0.25 ! toutes les 6h
1718	c
1719	cIM200505 ecrit_hf2mth = ecrit_day/ecrit_hf*30
1720	c
1721	cIM200505 IF (ok_journe) THEN
1722	cIM200505 WRITE(lunout,*)'La frequence de sortie hf est de ',
1723	cIM200505 . ecrit_hf
1724	cIM200505 ENDIF
1725	cIM 130904 END
1726	ccc ecrit_ins = NINT(86400./dtime *0.5) ! 2 fois par jour
1727	ccc ecrit_ins = NINT(86400./dtime *0.25) ! 4 fois par jour
1728	c ecrit_ins = NINT(86400./dtime/48.) ! a chaque pas de temps ==> PB. dans time_counter pour 1mois
1729	c ecrit_ins = NINT(86400./dtime/12.) ! toutes les deux heures
1730	cIM200505 ecrit_ins = NINT(86400./dtime/8.) ! toutes les trois heures
1731	cIM200505 IF (ok_instan) THEN
1732	cIM200505 WRITE(lunout,*)'La frequence de sortie instant. est de ',
1733	cIM200505 . ecrit_ins
1734	cIM200505 ENDIF
1735	cIM200505 ecrit_reg = NINT(86400./dtime *0.25) ! 4 fois par jour
1736	cIM200505 IF (ok_region) THEN
1737	cIM200505 WRITE(lunout,*)'La frequence de sortie region est de ',
1738	cIM200505 . ecrit_reg
1739	cIM200505 ENDIF
1740	cIM 030306 BEG
1741	cIM ecrit_hf2mth = nombre de pas de temps de calcul de hf par mois apres lequel on ecrit
1742	cIM : ne pas modifier ecrit_hf2mth
1743	c
1744	ecrit_hf2mth = 30*1/ecrit_hf
1745	c ecrit_ins en secondes, chaque pas de temps de la physique
1746	ecrit_ins = dtime
1747	cIM on passe les frequences de jours en secondes : ecrit_ins, ecrit_hf, ecrit_day, ecrit_mth, ecrit_tra, ecrit_reg
1748	ecrit_hf = ecrit_hf * un_jour
1749	!IM
1750	IF(ecrit_day.LE.1.) THEN
1751	ecrit_day = ecrit_day * un_jour !en secondes
1752	ENDIF
1753	!IM
1754	ecrit_mth = ecrit_mth * un_jour
1755	ecrit_reg = ecrit_reg * un_jour
1756	ecrit_tra = ecrit_tra * un_jour
1757	ecrit_ISCCP = ecrit_ISCCP * un_jour
1758	c
1759	PRINT*,'physiq ecrit_ hf day mth reg tra ISCCP',ecrit_hf,
1760	. ecrit_day,ecrit_mth,ecrit_reg,ecrit_tra,ecrit_ISCCP
1761	cIM 030306 END
1762
1763	capemaxcels = 't_max(X)'
1764	t2mincels = 't_min(X)'
1765	t2maxcels = 't_max(X)'
1766	tinst = 'inst(X)'
1767	tave = 'ave(X)'
1768	cIM cf. AM 081204 BEG
1769	write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
1770	cIM cf. AM 081204 END
1771	c
1772	c=============================================================
1773	c Initialisation des sorties
1774	c=============================================================
1775
1776	#ifdef CPP_IOIPSL
1777
1778	#ifdef histhf
1779	#include "ini_histhf.h"
1780	#endif
1781
1782	#ifdef histday
1783	#include "ini_histday.h"
1784	cIM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
1785	c#include "ini_bilKP_ins.h"
1786	c#include "ini_bilKP_ave.h"
1787	#endif
1788
1789	#ifdef histmth
1790	#include "ini_histmth.h"
1791	#endif
1792
1793	#ifdef histins
1794	#include "ini_histins.h"
1795	#endif
1796
1797	#ifdef histISCCP
1798	#include "ini_histISCCP.h"
1799	#endif
1800
1801	#ifdef histmthNMC
1802	#include "ini_histmthNMC.h"
1803	#endif
1804
1805	#include "ini_histday_seri.h"
1806
1807	#include "ini_paramLMDZ_phy.h"
1808
1809	#endif
1810
1811	cXXXPB Positionner date0 pour initialisation de ORCHIDEE
1812	date0 = zjulian
1813	C date0 = day_ini
1814	WRITE(,) 'physiq date0 : ',date0
1815	c
1816	c
1817	c
1818	c Prescrire l'ozone dans l'atmosphere
1819	c
1820	c
1821	cc DO i = 1, klon
1822	cc DO k = 1, klev
1823	cc CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1824	cc ENDDO
1825	cc ENDDO
1826	c
1827	#ifdef INCA
1828	call VTe(VTphysiq)
1829	call VTb(VTinca)
1830	iii = MOD(NINT(xjour),360)
1831	calday = FLOAT(iii) + gmtime
1832	WRITE(lunout,*) 'initial time ', xjour, calday
1833	#ifdef INCAINFO
1834	WRITE(lunout,*) 'Appel CHEMINI ...'
1835	#endif
1836	CALL chemini(
1837	$ rg,
1838	$ ra,
1839	$ airephy,
1840	$ rlat,
1841	$ rlon,
1842	$ presnivs,
1843	$ calday,
1844	$ klon,
1845	$ nqmax,
1846	$ pdtphys,
1847	$ annee_ref,
1848	$ day_ini)
1849	#ifdef INCAINFO
1850	WRITE(lunout,*) 'OK.'
1851	#endif
1852	call VTe(VTinca)
1853	call VTb(VTphysiq)
1854	#endif
1855	c
1856	ENDIF
1857	c
1858	c ************** Fin de IF ( debut ) *************
1859	c
1860	c
1861	c Mettre a zero des variables de sortie (pour securite)
1862	c
1863	DO i = 1, klon
1864	d_ps(i) = 0.0
1865	ENDDO
1866	DO k = 1, klev
1867	DO i = 1, klon
1868	d_t(i,k) = 0.0
1869	d_u(i,k) = 0.0
1870	d_v(i,k) = 0.0
1871	ENDDO
1872	ENDDO
1873	DO iq = 1, nqmax
1874	DO k = 1, klev
1875	DO i = 1, klon
1876	d_qx(i,k,iq) = 0.0
1877	ENDDO
1878	ENDDO
1879	ENDDO
1880	da(:,:)=0.
1881	mp(:,:)=0.
1882	phi(:,:,:)=0.
1883	c
1884	c Ne pas affecter les valeurs entrees de u, v, h, et q
1885	c
1886	DO k = 1, klev
1887	DO i = 1, klon
1888	t_seri(i,k) = t(i,k)
1889	u_seri(i,k) = u(i,k)
1890	v_seri(i,k) = v(i,k)
1891	q_seri(i,k) = qx(i,k,ivap)
1892	ql_seri(i,k) = qx(i,k,iliq)
1893	qs_seri(i,k) = 0.
1894	ENDDO
1895	ENDDO
1896	IF (nqmax.GE.3) THEN
1897	DO iq = 3, nqmax
1898	DO k = 1, klev
1899	DO i = 1, klon
1900	tr_seri(i,k,iq-2) = qx(i,k,iq)
1901	ENDDO
1902	ENDDO
1903	ENDDO
1904	ELSE
1905	DO k = 1, klev
1906	DO i = 1, klon
1907	tr_seri(i,k,1) = 0.0
1908	ENDDO
1909	ENDDO
1910	ENDIF
1911	C
1912	DO i = 1, klon
1913	ztsol(i) = 0.
1914	ENDDO
1915	DO nsrf = 1, nbsrf
1916	DO i = 1, klon
1917	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1918	ENDDO
1919	ENDDO
1920	cIM
1921	IF (ip_ebil_phy.ge.1) THEN
1922	ztit='after dynamic'
1923	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
1924	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1925	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1926	C Comme les tendances de la physique sont ajoute dans la dynamique,
1927	C on devrait avoir que la variation d'entalpie par la dynamique
1928	C est egale a la variation de la physique au pas de temps precedent.
1929	C Donc la somme de ces 2 variations devrait etre nulle.
1930	call diagphy(airephy,ztit,ip_ebil_phy
1931	e , zero_v, zero_v, zero_v, zero_v, zero_v
1932	e , zero_v, zero_v, zero_v, ztsol
1933	e , d_h_vcol+d_h_vcol_phy, d_qt, 0.
1934	s , fs_bound, fq_bound )
1935	END IF
1936
1937	c Diagnostiquer la tendance dynamique
1938	c
1939	IF (ancien_ok) THEN
1940	DO k = 1, klev
1941	DO i = 1, klon
1942	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
1943	d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
1944	ENDDO
1945	ENDDO
1946	ELSE
1947	DO k = 1, klev
1948	DO i = 1, klon
1949	d_t_dyn(i,k) = 0.0
1950	d_q_dyn(i,k) = 0.0
1951	ENDDO
1952	ENDDO
1953	ancien_ok = .TRUE.
1954	ENDIF
1955	c
1956	c Ajouter le geopotentiel du sol:
1957	c
1958	DO k = 1, klev
1959	DO i = 1, klon
1960	zphi(i,k) = pphi(i,k) + pphis(i)
1961	ENDDO
1962	ENDDO
1963	c
1964	c Verifier les temperatures
1965	c
1966	cIM BEG
1967	IF (check) THEN
1968	amn=MIN(ftsol(1,is_ter),1000.)
1969	amx=MAX(ftsol(1,is_ter),-1000.)
1970	DO i=2, klon
1971	amn=MIN(ftsol(i,is_ter),amn)
1972	amx=MAX(ftsol(i,is_ter),amx)
1973	ENDDO
1974	c
1975	PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
1976	ENDIF !(check) THEN
1977	cIM END
1978	c
1979	CALL hgardfou(t_seri,ftsol,'debutphy')
1980	c
1981	cIM BEG
1982	IF (check) THEN
1983	amn=MIN(ftsol(1,is_ter),1000.)
1984	amx=MAX(ftsol(1,is_ter),-1000.)
1985	DO i=2, klon
1986	amn=MIN(ftsol(i,is_ter),amn)
1987	amx=MAX(ftsol(i,is_ter),amx)
1988	ENDDO
1989	c
1990	PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
1991	ENDIF !(check) THEN
1992	cIM END
1993	c
1994	c Incrementer le compteur de la physique
1995	c
1996	itap = itap + 1
1997	julien = MOD(NINT(xjour),360)
1998	if (julien .eq. 0) julien = 360
1999	c
2000	c Mettre en action les conditions aux limites (albedo, sst, etc.).
2001	c Prescrire l'ozone et calculer l'albedo sur l'ocean.
2002	c
2003	IF (MOD(itap-1,lmt_pas) .EQ. 0) THEN
2004	WRITE(lunout,*)' PHYS cond julien ',julien
2005	CALL ozonecm( FLOAT(julien), rlat, paprs, wo)
2006	ENDIF
2007	c
2008	c Re-evaporer l'eau liquide nuageuse
2009	c
2010	DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse
2011	DO i = 1, klon
2012	zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
2013	c zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
2014	zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
2015	zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
2016	zb = MAX(0.0,ql_seri(i,k))
2017	za = - MAX(0.0,ql_seri(i,k))
2018	. * (zlvdcp(1.-zdelta)+zlsdcpzdelta)
2019	t_seri(i,k) = t_seri(i,k) + za
2020	q_seri(i,k) = q_seri(i,k) + zb
2021	ql_seri(i,k) = 0.0
2022	d_t_eva(i,k) = za
2023	d_q_eva(i,k) = zb
2024	ENDDO
2025	ENDDO
2026	cIM
2027	IF (ip_ebil_phy.ge.2) THEN
2028	ztit='after reevap'
2029	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,1,dtime
2030	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2031	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2032	call diagphy(airephy,ztit,ip_ebil_phy
2033	e , zero_v, zero_v, zero_v, zero_v, zero_v
2034	e , zero_v, zero_v, zero_v, ztsol
2035	e , d_h_vcol, d_qt, d_ec
2036	s , fs_bound, fq_bound )
2037	C
2038	END IF
2039
2040	c
2041	C calculs necessaires au calcul de l'albedo dans l'interface
2042	c
2043	CALL orbite(FLOAT(julien),zlongi,dist)
2044	IF (cycle_diurne) THEN
2045	zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
2046	CALL zenang(zlongi,gmtime,zdtime,rlat,rlon,rmu0,fract)
2047	ELSE
2048	rmu0 = -999.999
2049	ENDIF
2050
2051	if (mydebug) then
2052	call writefield_phy('u_seri',u_seri,llm)
2053	call writefield_phy('v_seri',v_seri,llm)
2054	call writefield_phy('t_seri',t_seri,llm)
2055	call writefield_phy('q_seri',q_seri,llm)
2056	endif
2057
2058	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2059	c Appel au pbl_surface : Planetary Boudary Layer et Surface
2060	c Cela implique tous les interactions des sous-surfaces et la partie diffusion
2061	c turbulent du couche limit.
2062	c
2063	c Certains varibales de sorties de pbl_surface sont utiliser que pour
2064	c ecriture des fihiers hist_XXXX.nc, ces sont :
2065	c qsol, zq2m, s_pblh, s_lcl,
2066	c s_capCL, s_oliqCL, s_cteiCL,s_pblT,
2067	c s_therm, s_trmb1, s_trmb2, s_trmb3,
2068	c zxrugs, zu10m, zv10m, fder,
2069	c zxqsurf, rh2m, zxfluxu, zxfluxv,
2070	c frugs, agesno, fsollw, fsolsw,
2071	c d_ts, fevap, fluxlat, t2m,
2072	c wfbils, wfbilo, fluxt, fluxu, fluxv,
2073	c
2074	c Certains ne sont pas utiliser du tout :
2075	c dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
2076	c
2077	CALL pbl_surface(
2078	e dtime, date0, itap, julien,
2079	e debut, lafin,
2080	e rlon, rlat, rugoro, rmu0,
2081	e rain_fall, snow_fall, solsw, sollw,
2082	e t_seri, q_seri, u_seri, v_seri,
2083	e pplay, paprs, pctsrf,
2084	+ ftsol, falbe, falblw, u10m, v10m,
2085	s sollwdown, cdragh, cdragm, yu1, yv1,
2086	s albsol, albsollw, sens, evap,
2087	s zxtsol, zxfluxlat, zt2m, qsat2m,
2088	s d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf,
2089	s ycoefh, pctsrf_new,
2090	d qsol, zq2m, s_pblh, s_lcl,
2091	d s_capCL, s_oliqCL, s_cteiCL,s_pblT,
2092	d s_therm, s_trmb1, s_trmb2, s_trmb3,
2093	d zxrugs, zu10m, zv10m, fder,
2094	d zxqsurf, rh2m, zxfluxu, zxfluxv,
2095	d frugs, agesno, fsollw, fsolsw,
2096	d d_ts, fevap, fluxlat, t2m,
2097	d wfbils, wfbilo, fluxt, fluxu, fluxv,
2098	- dsens, devap, zxsnow,
2099	- zxfluxt, zxfluxq, q2m, fluxq, pbl_tke )
2100	c
2101	c
2102	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2103
2104	pctsrf(:,:) = pctsrf_new(:,:)
2105
2106	DO k = 1, klev
2107	DO i = 1, klon
2108	t_seri(i,k) = t_seri(i,k) + d_t_vdf(i,k)
2109	q_seri(i,k) = q_seri(i,k) + d_q_vdf(i,k)
2110	u_seri(i,k) = u_seri(i,k) + d_u_vdf(i,k)
2111	v_seri(i,k) = v_seri(i,k) + d_v_vdf(i,k)
2112	ENDDO
2113	ENDDO
2114
2115	if (mydebug) then
2116	call writefield_phy('u_seri',u_seri,llm)
2117	call writefield_phy('v_seri',v_seri,llm)
2118	call writefield_phy('t_seri',t_seri,llm)
2119	call writefield_phy('q_seri',q_seri,llm)
2120	endif
2121
2122
2123	IF (ip_ebil_phy.ge.2) THEN
2124	ztit='after surface_main'
2125	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2126	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2127	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2128	call diagphy(airephy,ztit,ip_ebil_phy
2129	e , zero_v, zero_v, zero_v, zero_v, sens
2130	e , evap , zero_v, zero_v, ztsol
2131	e , d_h_vcol, d_qt, d_ec
2132	s , fs_bound, fq_bound )
2133	END IF
2134
2135	c
2136	c Appeler la convection (au choix)
2137	c
2138	DO k = 1, klev
2139	DO i = 1, klon
2140	conv_q(i,k) = d_q_dyn(i,k)
2141	. + d_q_vdf(i,k)/dtime
2142	conv_t(i,k) = d_t_dyn(i,k)
2143	. + d_t_vdf(i,k)/dtime
2144	ENDDO
2145	ENDDO
2146	IF (check) THEN
2147	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2148	WRITE(lunout,*) "avantcon=", za
2149	ENDIF
2150	zx_ajustq = .FALSE.
2151	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
2152	IF (zx_ajustq) THEN
2153	DO i = 1, klon
2154	z_avant(i) = 0.0
2155	ENDDO
2156	DO k = 1, klev
2157	DO i = 1, klon
2158	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k))
2159	. *(paprs(i,k)-paprs(i,k+1))/RG
2160	ENDDO
2161	ENDDO
2162	ENDIF
2163	IF (iflag_con.EQ.1) THEN
2164	stop'reactiver le call conlmd dans physiq.F'
2165	c CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
2166	c . d_t_con, d_q_con,
2167	c . rain_con, snow_con, ibas_con, itop_con)
2168	ELSE IF (iflag_con.EQ.2) THEN
2169	CALL conflx(dtime, paprs, pplay, t_seri, q_seri,
2170	e conv_t, conv_q, -evap, omega,
2171	s d_t_con, d_q_con, rain_con, snow_con,
2172	s pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
2173	s kcbot, kctop, kdtop, pmflxr, pmflxs)
2174	WHERE (rain_con < 0.) rain_con = 0.
2175	WHERE (snow_con < 0.) snow_con = 0.
2176	DO i = 1, klon
2177	ibas_con(i) = klev+1 - kcbot(i)
2178	itop_con(i) = klev+1 - kctop(i)
2179	ENDDO
2180	ELSE IF (iflag_con.GE.3) THEN
2181	c nb of tracers for the KE convection:
2182	c MAF la partie traceurs est faite dans phytrac
2183	c on met ntra=1 pour limiter les appels mais on peut
2184	c supprimer les calculs / ftra.
2185	ntra = 1
2186	c sb, oct02:
2187	c Schema de convection modularise et vectorise:
2188	c (driver commun aux versions 3 et 4)
2189	c
2190	IF (ok_cvl) THEN ! new driver for convectL
2191
2192	CALL concvl (iflag_con,
2193	. dtime,paprs,pplay,t_seri,q_seri,
2194	. u_seri,v_seri,tr_seri,ntra,
2195	. ema_work1,ema_work2,
2196	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2197	. rain_con, snow_con, ibas_con, itop_con,
2198	. upwd,dnwd,dnwd0,
2199	. Ma,cape,tvp,iflagctrl,
2200	. pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd,
2201	. pmflxr,pmflxs,
2202	. da,phi,mp)
2203
2204	cIM cf. FH
2205	clwcon0=qcondc
2206	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
2207
2208	ELSE ! ok_cvl
2209	c MAF conema3 ne contient pas les traceurs
2210	CALL conema3 (dtime,
2211	. paprs,pplay,t_seri,q_seri,
2212	. u_seri,v_seri,tr_seri,ntra,
2213	. ema_work1,ema_work2,
2214	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2215	. rain_con, snow_con, ibas_con, itop_con,
2216	. upwd,dnwd,dnwd0,bas,top,
2217	. Ma,cape,tvp,rflag,
2218	. pbase
2219	. ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr
2220	. ,clwcon0)
2221
2222	ENDIF ! ok_cvl
2223
2224	c
2225	c Correction precip
2226	rain_con = rain_con * cvl_corr
2227	snow_con = snow_con * cvl_corr
2228	c
2229
2230	IF (.NOT. ok_gust) THEN
2231	do i = 1, klon
2232	wd(i)=0.0
2233	enddo
2234	ENDIF
2235
2236	c =================================================================== c
2237	c Calcul des proprietes des nuages convectifs
2238	c
2239	DO k = 1, klev
2240	DO i = 1, klon
2241	zx_t = t_seri(i,k)
2242	IF (thermcep) THEN
2243	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2244	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2245	zx_qs = MIN(0.5,zx_qs)
2246	zcor = 1./(1.-retv*zx_qs)
2247	zx_qs = zx_qs*zcor
2248	ELSE
2249	IF (zx_t.LT.t_coup) THEN
2250	zx_qs = qsats(zx_t)/pplay(i,k)
2251	ELSE
2252	zx_qs = qsatl(zx_t)/pplay(i,k)
2253	ENDIF
2254	ENDIF
2255	zqsat(i,k)=zx_qs
2256	ENDDO
2257	ENDDO
2258
2259	c calcul des proprietes des nuages convectifs
2260	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
2261	call clouds_gno
2262	s (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
2263
2264	c =================================================================== c
2265
2266	DO i = 1, klon
2267	ema_pcb(i) = pbase(i)
2268	ENDDO
2269	DO i = 1, klon
2270	ema_pct(i) = paprs(i,itop_con(i))
2271	if (itop_con(i).gt.klev-3) then
2272	print*,'La convection monte trop haut '
2273	print*,'itop_con(,',i,',)=',itop_con(i)
2274	endif
2275	ENDDO
2276	DO i = 1, klon
2277	ema_cbmf(i) = ema_workcbmf(i)
2278	ENDDO
2279	ELSE
2280	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
2281	CALL abort
2282	ENDIF
2283
2284	c CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
2285	c . d_u_con, d_v_con)
2286
2287	DO k = 1, klev
2288	DO i = 1, klon
2289	t_seri(i,k) = t_seri(i,k) + d_t_con(i,k)
2290	q_seri(i,k) = q_seri(i,k) + d_q_con(i,k)
2291	u_seri(i,k) = u_seri(i,k) + d_u_con(i,k)
2292	v_seri(i,k) = v_seri(i,k) + d_v_con(i,k)
2293	ENDDO
2294	ENDDO
2295
2296	if (mydebug) then
2297	call writefield_phy('u_seri',u_seri,llm)
2298	call writefield_phy('v_seri',v_seri,llm)
2299	call writefield_phy('t_seri',t_seri,llm)
2300	call writefield_phy('q_seri',q_seri,llm)
2301	endif
2302
2303	cIM
2304	IF (ip_ebil_phy.ge.2) THEN
2305	ztit='after convect'
2306	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2307	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2308	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2309	call diagphy(airephy,ztit,ip_ebil_phy
2310	e , zero_v, zero_v, zero_v, zero_v, zero_v
2311	e , zero_v, rain_con, snow_con, ztsol
2312	e , d_h_vcol, d_qt, d_ec
2313	s , fs_bound, fq_bound )
2314	END IF
2315	C
2316	IF (check) THEN
2317	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2318	WRITE(lunout,*)"aprescon=", za
2319	zx_t = 0.0
2320	za = 0.0
2321	DO i = 1, klon
2322	za = za + airephy(i)/FLOAT(klon)
2323	zx_t = zx_t + (rain_con(i)+
2324	. snow_con(i))*airephy(i)/FLOAT(klon)
2325	ENDDO
2326	zx_t = zx_t/za*dtime
2327	WRITE(lunout,*)"Precip=", zx_t
2328	ENDIF
2329	IF (zx_ajustq) THEN
2330	DO i = 1, klon
2331	z_apres(i) = 0.0
2332	ENDDO
2333	DO k = 1, klev
2334	DO i = 1, klon
2335	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k))
2336	. *(paprs(i,k)-paprs(i,k+1))/RG
2337	ENDDO
2338	ENDDO
2339	DO i = 1, klon
2340	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime)
2341	. /z_apres(i)
2342	ENDDO
2343	DO k = 1, klev
2344	DO i = 1, klon
2345	IF (z_factor(i).GT.(1.0+1.0E-08) .OR.
2346	. z_factor(i).LT.(1.0-1.0E-08)) THEN
2347	q_seri(i,k) = q_seri(i,k) * z_factor(i)
2348	ENDIF
2349	ENDDO
2350	ENDDO
2351	ENDIF
2352	zx_ajustq=.FALSE.
2353	c
2354	c===================================================================
2355	c Convection seche (thermiques ou ajustement)
2356	c===================================================================
2357	c
2358	call stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri
2359	s ,seuil_inversion,weak_inversion,dthmin)
2360
2361
2362
2363	d_t_ajsb(:,:)=0.
2364	d_q_ajsb(:,:)=0.
2365	d_t_ajs(:,:)=0.
2366	d_u_ajs(:,:)=0.
2367	d_v_ajs(:,:)=0.
2368	d_q_ajs(:,:)=0.
2369	clwcon0th(:,:)=0.
2370	fm_therm(:,:)=0.
2371	entr_therm(:,:)=0.
2372	c
2373	IF(prt_level>9)WRITE(lunout,*)
2374	. 'AVANT LA CONVECTION SECHE , iflag_thermals='
2375	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2376	if(iflag_thermals.lt.0) then
2377	c Rien
2378	c ====
2379	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
2380
2381
2382	else
2383
2384	c Thermiques
2385	c ==========
2386	IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals='
2387	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2388	print*,'JUSTE AVANT , iflag_thermals='
2389	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2390
2391
2392	if (iflag_thermals.gt.1) then
2393	call calltherm(pdtphys
2394	s ,pplay,paprs,pphi,weak_inversion
2395	s ,u_seri,v_seri,t_seri,q_seri,zqsat,debut
2396	s ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
2397	s ,fm_therm,entr_therm,zqasc,clwcon0th,lmax_th,ratqscth,
2398	s ratqsdiff,zqsatth)
2399	endif
2400
2401	! call calltherm(pdtphys
2402	! s ,pplay,paprs,pphi
2403	! s ,u_seri,v_seri,t_seri,q_seri
2404	! s ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
2405	! s ,fm_therm,entr_therm)
2406
2407	c Ajustement sec
2408	c ==============
2409
2410	! Dans le cas où on active les thermiques, on fait partir l'ajustement
2411	! a partir du sommet des thermiques.
2412	! Dans le cas contraire, on demarre au niveau 1.
2413
2414	if (iflag_thermals.ge.13.or.iflag_thermals.eq.0) then
2415
2416	if(iflag_thermals.eq.0) then
2417	IF(prt_level>9)WRITE(lunout,*)'ajsec'
2418	limbas(:)=1
2419	else
2420	limbas(:)=lmax_th(:)
2421	endif
2422
2423	! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
2424	! pour des test de convergence numerique.
2425	! Le nouveau ajsec est a priori mieux, meme pour le cas
2426	! iflag_thermals = 0 (l'ancienne version peut faire des tendances
2427	! non nulles numeriquement pour des mailles non concernees.
2428
2429	if (iflag_thermals.eq.0) then
2430	CALL ajsec_convV2(paprs, pplay, t_seri,q_seri
2431	s , d_t_ajsb, d_q_ajsb)
2432	else
2433	CALL ajsec(paprs, pplay, t_seri,q_seri,limbas
2434	s , d_t_ajsb, d_q_ajsb)
2435	endif
2436
2437	t_seri(:,:) = t_seri(:,:) + d_t_ajsb(:,:)
2438	q_seri(:,:) = q_seri(:,:) + d_q_ajsb(:,:)
2439	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
2440	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
2441
2442	endif
2443
2444	endif
2445	c
2446	c===================================================================
2447	cIM
2448	IF (ip_ebil_phy.ge.2) THEN
2449	ztit='after dry_adjust'
2450	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2451	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2452	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2453	END IF
2454
2455
2456	c-------------------------------------------------------------------------
2457	c Caclul des ratqs
2458	c-------------------------------------------------------------------------
2459
2460	c print*,'calcul des ratqs'
2461	c ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
2462	c ----------------
2463	c on ecrase le tableau ratqsc calcule par clouds_gno
2464	if (iflag_cldcon.eq.1) then
2465	do k=1,klev
2466	do i=1,klon
2467	if(ptconv(i,k)) then
2468	ratqsc(i,k)=ratqsbas
2469	s +fact_cldcon*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
2470	else
2471	ratqsc(i,k)=0.
2472	endif
2473	enddo
2474	enddo
2475	else if (iflag_cldcon.eq.4) then
2476	ptconvth(:,:)=.false.
2477	ratqsc(:,:)=0.
2478	print*,'avant clouds_gno thermique'
2479	call clouds_gno
2480	s (klon,klev,q_seri,zqsat,clwcon0th,ptconvth,ratqsc,rnebcon0th)
2481	print*,' CLOUDS_GNO OK'
2482	endif
2483
2484	c ratqs stables
2485	c -------------
2486
2487	if (iflag_ratqs.eq.0) then
2488
2489	! Le cas iflag_ratqs=0 correspond a la version IPCC 2005 du modele.
2490	do k=1,klev
2491	do i=1, klon
2492	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2493	s min((paprs(i,1)-pplay(i,k))/(paprs(i,1)-30000.),1.)
2494	enddo
2495	enddo
2496
2497	! Pour iflag_ratqs=1 ou 2, le ratqs est constant au dessus de
2498	! 300 hPa (ratqshaut), varie lineariement en fonction de la pression
2499	! entre 600 et 300 hPa et est soit constant (ratqsbas) pour iflag_ratqs=1
2500	! soit lineaire (entre 0 a la surface et ratqsbas) pour iflag_ratqs=2
2501	! Il s'agit de differents tests dans la phase de reglage du modele
2502	! avec thermiques.
2503
2504	else if (iflag_ratqs.eq.1) then
2505
2506	do k=1,klev
2507	do i=1, klon
2508	if (pplay(i,k).ge.60000.) then
2509	ratqss(i,k)=ratqsbas
2510	else if ((pplay(i,k).ge.30000.).and.
2511	s (pplay(i,k).lt.60000.)) then
2512	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2513	s (60000.-pplay(i,k))/(60000.-30000.)
2514	else
2515	ratqss(i,k)=ratqshaut
2516	endif
2517	enddo
2518	enddo
2519
2520	else
2521
2522	do k=1,klev
2523	do i=1, klon
2524	if (pplay(i,k).ge.60000.) then
2525	ratqss(i,k)=ratqsbas
2526	s *(paprs(i,1)-pplay(i,k))/(paprs(i,1)-60000.)
2527	else if ((pplay(i,k).ge.30000.).and.
2528	s (pplay(i,k).lt.60000.)) then
2529	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2530	s (60000.-pplay(i,k))/(60000.-30000.)
2531	else
2532	ratqss(i,k)=ratqshaut
2533	endif
2534	enddo
2535	enddo
2536	endif
2537
2538
2539
2540
2541	c ratqs final
2542	c -----------
2543
2544	if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2
2545	s .or.iflag_cldcon.eq.4) then
2546
2547	! On ajoute une constante au ratqsc*2 pour tenir compte de
2548	! fluctuations turbulentes de petite echelle
2549
2550	do k=1,klev
2551	do i=1,klon
2552	if ((fm_therm(i,k).gt.1.e-10)) then
2553	ratqsc(i,k)=sqrt(ratqsc(i,k)2+0.052)
2554	endif
2555	enddo
2556	enddo
2557
2558	! les ratqs sont une conbinaison de ratqss et ratqsc
2559	! 1800s, en dur pour le moment, est le temps de
2560	! relaxation des ratqs
2561
2562	facteur=exp(-pdtphys/1800.)
2563
2564	print*,'WARNING ratqs a revoir '
2565	ratqs(:,:)=ratqsc(:,:)(1.-facteur)+ratqs(:,:)facteur
2566	ratqs(:,:)=sqrt(ratqs(:,:)2+ratqss(:,:)2)
2567
2568	else
2569	! on ne prend que le ratqs stable pour fisrtilp
2570	ratqs(:,:)=ratqss(:,:)
2571	endif
2572
2573
2574	c
2575	c Appeler le processus de condensation a grande echelle
2576	c et le processus de precipitation
2577	c-------------------------------------------------------------------------
2578	CALL fisrtilp(dtime,paprs,pplay,
2579	. t_seri, q_seri,ptconv,ratqs,
2580	. d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq,
2581	. rain_lsc, snow_lsc,
2582	. pfrac_impa, pfrac_nucl, pfrac_1nucl,
2583	. frac_impa, frac_nucl,
2584	. prfl, psfl, rhcl)
2585
2586	WHERE (rain_lsc < 0) rain_lsc = 0.
2587	WHERE (snow_lsc < 0) snow_lsc = 0.
2588	DO k = 1, klev
2589	DO i = 1, klon
2590	t_seri(i,k) = t_seri(i,k) + d_t_lsc(i,k)
2591	q_seri(i,k) = q_seri(i,k) + d_q_lsc(i,k)
2592	ql_seri(i,k) = ql_seri(i,k) + d_ql_lsc(i,k)
2593	cldfra(i,k) = rneb(i,k)
2594	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
2595	ENDDO
2596	ENDDO
2597	IF (check) THEN
2598	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2599	WRITE(lunout,*)"apresilp=", za
2600	zx_t = 0.0
2601	za = 0.0
2602	DO i = 1, klon
2603	za = za + airephy(i)/FLOAT(klon)
2604	zx_t = zx_t + (rain_lsc(i)
2605	. + snow_lsc(i))*airephy(i)/FLOAT(klon)
2606	ENDDO
2607	zx_t = zx_t/za*dtime
2608	WRITE(lunout,*)"Precip=", zx_t
2609	ENDIF
2610	cIM
2611	IF (ip_ebil_phy.ge.2) THEN
2612	ztit='after fisrt'
2613	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2614	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2615	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2616	call diagphy(airephy,ztit,ip_ebil_phy
2617	e , zero_v, zero_v, zero_v, zero_v, zero_v
2618	e , zero_v, rain_lsc, snow_lsc, ztsol
2619	e , d_h_vcol, d_qt, d_ec
2620	s , fs_bound, fq_bound )
2621	END IF
2622
2623	if (mydebug) then
2624	call writefield_phy('u_seri',u_seri,llm)
2625	call writefield_phy('v_seri',v_seri,llm)
2626	call writefield_phy('t_seri',t_seri,llm)
2627	call writefield_phy('q_seri',q_seri,llm)
2628	endif
2629
2630	c
2631	c-------------------------------------------------------------------
2632	c PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
2633	c-------------------------------------------------------------------
2634
2635	c 1. NUAGES CONVECTIFS
2636	c
2637	cIM cf FH
2638	c IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
2639	IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
2640	snow_tiedtke=0.
2641	c print*,'avant calcul de la pseudo precip '
2642	c print*,'iflag_cldcon',iflag_cldcon
2643	if (iflag_cldcon.eq.-1) then
2644	rain_tiedtke=rain_con
2645	else
2646	c print*,'calcul de la pseudo precip '
2647	rain_tiedtke=0.
2648	c print*,'calcul de la pseudo precip 0'
2649	do k=1,klev
2650	do i=1,klon
2651	if (d_q_con(i,k).lt.0.) then
2652	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys
2653	s *(paprs(i,k)-paprs(i,k+1))/rg
2654	endif
2655	enddo
2656	enddo
2657	endif
2658	c
2659	c call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
2660	c
2661
2662	c Nuages diagnostiques pour Tiedtke
2663	CALL diagcld1(paprs,pplay,
2664	cIM cf FH . rain_con,snow_con,ibas_con,itop_con,
2665	. rain_tiedtke,snow_tiedtke,ibas_con,itop_con,
2666	. diafra,dialiq)
2667	DO k = 1, klev
2668	DO i = 1, klon
2669	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2670	cldliq(i,k) = dialiq(i,k)
2671	cldfra(i,k) = diafra(i,k)
2672	ENDIF
2673	ENDDO
2674	ENDDO
2675
2676	ELSE IF (iflag_cldcon.ge.3) THEN
2677	c On prend pour les nuages convectifs le max du calcul de la
2678	c convection et du calcul du pas de temps precedent diminue d'un facteur
2679	c facttemps
2680	facteur = pdtphys *facttemps
2681	do k=1,klev
2682	do i=1,klon
2683	rnebcon(i,k)=rnebcon(i,k)*facteur
2684	if (rnebcon0(i,k)clwcon0(i,k).gt.rnebcon(i,k)clwcon(i,k))
2685	s then
2686	rnebcon(i,k)=rnebcon0(i,k)
2687	clwcon(i,k)=clwcon0(i,k)
2688	endif
2689	enddo
2690	enddo
2691
2692	c
2693	cjq - introduce the aerosol direct and first indirect radiative forcings
2694	cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
2695	IF (ok_ade.OR.ok_aie) THEN
2696	! Get sulfate aerosol distribution
2697	CALL readsulfate(rjourvrai, debut, sulfate)
2698	CALL readsulfate_preind(rjourvrai, debut, sulfate_pi)
2699
2700	! Calculate aerosol optical properties (Olivier Boucher)
2701	CALL aeropt(pplay, paprs, t_seri, sulfate, rhcl,
2702	. tau_ae, piz_ae, cg_ae, aerindex)
2703	cym
2704	ELSE
2705	tau_ae(:,:,:)=0.0
2706	piz_ae(:,:,:)=0.0
2707	cg_ae(:,:,:)=0.0
2708	cym
2709	ENDIF
2710
2711	c
2712	cIM calcul nuages par le simulateur ISCCP
2713	c
2714	#ifdef histISCCP
2715	IF (ok_isccp) THEN
2716	cIM appel simulateur toutes les NINT(freq_ISCCP/dtime) heures
2717	IF (MOD(itap,NINT(freq_ISCCP/dtime)).EQ.0) THEN
2718	#include "calcul_simulISCCP.h"
2719	ENDIF !(MOD(itap,NINT(freq_ISCCP/dtime))
2720	ENDIF !ok_isccp
2721	#endif
2722
2723	c On prend la somme des fractions nuageuses et des contenus en eau
2724	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
2725	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
2726
2727	ENDIF
2728
2729	c
2730	c 2. NUAGES STARTIFORMES
2731	c
2732	IF (ok_stratus) THEN
2733	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
2734	DO k = 1, klev
2735	DO i = 1, klon
2736	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2737	cldliq(i,k) = dialiq(i,k)
2738	cldfra(i,k) = diafra(i,k)
2739	ENDIF
2740	ENDDO
2741	ENDDO
2742	ENDIF
2743	c
2744	c Precipitation totale
2745	c
2746	DO i = 1, klon
2747	rain_fall(i) = rain_con(i) + rain_lsc(i)
2748	snow_fall(i) = snow_con(i) + snow_lsc(i)
2749	ENDDO
2750	cIM
2751	IF (ip_ebil_phy.ge.2) THEN
2752	ztit="after diagcld"
2753	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2754	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2755	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2756	END IF
2757	c
2758	c Calculer l'humidite relative pour diagnostique
2759	c
2760	DO k = 1, klev
2761	DO i = 1, klon
2762	zx_t = t_seri(i,k)
2763	IF (thermcep) THEN
2764	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2765	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2766	zx_qs = MIN(0.5,zx_qs)
2767	zcor = 1./(1.-retv*zx_qs)
2768	zx_qs = zx_qs*zcor
2769	ELSE
2770	IF (zx_t.LT.t_coup) THEN
2771	zx_qs = qsats(zx_t)/pplay(i,k)
2772	ELSE
2773	zx_qs = qsatl(zx_t)/pplay(i,k)
2774	ENDIF
2775	ENDIF
2776	zx_rh(i,k) = q_seri(i,k)/zx_qs
2777	zqsat(i,k)=zx_qs
2778	ENDDO
2779	ENDDO
2780
2781	cIM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
2782	c equivalente a 2m (tpote) pour diagnostique
2783	c
2784	DO i = 1, klon
2785	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
2786	IF (thermcep) THEN
2787	IF(zt2m(i).LT.RTT) then
2788	Lheat=RLSTT
2789	ELSE
2790	Lheat=RLVTT
2791	ENDIF
2792	ELSE
2793	IF (zt2m(i).LT.RTT) THEN
2794	Lheat=RLSTT
2795	ELSE
2796	Lheat=RLVTT
2797	ENDIF
2798	ENDIF
2799	tpote(i) = tpot(i)*
2800	. EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
2801	ENDDO
2802
2803
2804	#ifdef INCA
2805	call VTe(VTphysiq)
2806	call VTb(VTinca)
2807	calday = FLOAT(julien) + gmtime
2808
2809	#ifdef INCA_AER
2810	call AEROSOL_METEO_CALC(calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs
2811	& ,prfl,psfl,pctsrf,airephy,xjour,rlat,rlon,u10m,v10m)
2812	#endif
2813
2814	#ifdef INCAINFO
2815	WRITE(lunout,*)'Appel CHEMHOOK_BEGIN ...'
2816	#endif
2817
2818	zxsnow_dummy(:) = 0.0
2819
2820	CALL chemhook_begin (calday,
2821	#if defined(INCA) && !defined(INCA_CH4) && !defined(INCA_NMHC) && !defined(INCA_AER)
2822	$ julien,
2823	$ gmtime,
2824	#endif
2825	$ pctsrf(1,1),
2826	$ rlat,
2827	$ rlon,
2828	$ airephy,
2829	$ paprs,
2830	$ pplay,
2831	$ ycoefh,
2832	$ pphi,
2833	$ t_seri,
2834	$ u,
2835	$ v,
2836	$ wo,
2837	$ q_seri,
2838	$ zxtsol,
2839	$ zxsnow_dummy,
2840	$ solsw,
2841	$ albsol,
2842	$ rain_fall,
2843	$ snow_fall,
2844	$ itop_con,
2845	$ ibas_con,
2846	$ cldfra,
2847	$ iim,
2848	$ jjm,
2849	#ifdef INCA_AER
2850	$ tr_seri,
2851	$ ftsol,
2852	$ paprs,
2853	$ cdragh,
2854	$ cdragm,
2855	$ pctsrf,
2856	$ pdtphys,
2857	$ itap)
2858	#else
2859	$ tr_seri)
2860	#endif
2861
2862
2863	#ifdef INCAINFO
2864	WRITE(lunout,*)'OK.'
2865	#endif
2866	call VTe(VTinca)
2867	call VTb(VTphysiq)
2868	#endif
2869	c
2870	c Calculer les parametres optiques des nuages et quelques
2871	c parametres pour diagnostiques:
2872	c
2873	if (ok_newmicro) then
2874	CALL newmicro (paprs, pplay,ok_newmicro,
2875	. t_seri, cldliq, cldfra, cldtau, cldemi,
2876	. cldh, cldl, cldm, cldt, cldq,
2877	. flwp, fiwp, flwc, fiwc,
2878	e ok_aie,
2879	e sulfate, sulfate_pi,
2880	e bl95_b0, bl95_b1,
2881	s cldtaupi, re, fl)
2882	else
2883	CALL nuage (paprs, pplay,
2884	. t_seri, cldliq, cldfra, cldtau, cldemi,
2885	. cldh, cldl, cldm, cldt, cldq,
2886	e ok_aie,
2887	e sulfate, sulfate_pi,
2888	e bl95_b0, bl95_b1,
2889	s cldtaupi, re, fl)
2890
2891	endif
2892	c
2893	c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
2894	c
2895	IF (MOD(itaprad,radpas).EQ.0) THEN
2896
2897	DO i = 1, klon
2898	albsol(i) = falbe(i,is_oce) * pctsrf(i,is_oce)
2899	. + falbe(i,is_lic) * pctsrf(i,is_lic)
2900	. + falbe(i,is_ter) * pctsrf(i,is_ter)
2901	. + falbe(i,is_sic) * pctsrf(i,is_sic)
2902	albsollw(i) = falblw(i,is_oce) * pctsrf(i,is_oce)
2903	. + falblw(i,is_lic) * pctsrf(i,is_lic)
2904	. + falblw(i,is_ter) * pctsrf(i,is_ter)
2905	. + falblw(i,is_sic) * pctsrf(i,is_sic)
2906	ENDDO
2907
2908	if (mydebug) then
2909	call writefield_phy('u_seri',u_seri,llm)
2910	call writefield_phy('v_seri',v_seri,llm)
2911	call writefield_phy('t_seri',t_seri,llm)
2912	call writefield_phy('q_seri',q_seri,llm)
2913	endif
2914
2915	CALL radlwsw ! nouveau rayonnement (compatible Arpege-IFS)
2916	e (dist, rmu0, fract,
2917	e paprs, pplay,zxtsol,albsol, albsollw, t_seri,q_seri,
2918	e wo,
2919	e cldfra, cldemi, cldtau,
2920	s heat,heat0,cool,cool0,radsol,albpla,
2921	s topsw,toplw,solsw,sollw,
2922	s sollwdown,
2923	s topsw0,toplw0,solsw0,sollw0,
2924	s lwdn0, lwdn, lwup0, lwup,
2925	s swdn0, swdn, swup0, swup,
2926	e ok_ade, ok_aie, ! new for aerosol radiative effects
2927	e tau_ae, piz_ae, cg_ae, ! ="=
2928	s topswad, solswad, ! ="=
2929	e cldtaupi, ! ="=
2930	s topswai, solswai) ! ="=
2931	itaprad = 0
2932	ENDIF
2933	itaprad = itaprad + 1
2934	c
2935	c Ajouter la tendance des rayonnements (tous les pas)
2936	c
2937	DO k = 1, klev
2938	DO i = 1, klon
2939	t_seri(i,k) = t_seri(i,k)
2940	. + (heat(i,k)-cool(i,k)) * dtime/86400.
2941	ENDDO
2942	ENDDO
2943	c
2944	if (mydebug) then
2945	call writefield_phy('u_seri',u_seri,llm)
2946	call writefield_phy('v_seri',v_seri,llm)
2947	call writefield_phy('t_seri',t_seri,llm)
2948	call writefield_phy('q_seri',q_seri,llm)
2949	endif
2950
2951	cIM
2952	IF (ip_ebil_phy.ge.2) THEN
2953	ztit='after rad'
2954	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2955	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2956	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2957	call diagphy(airephy,ztit,ip_ebil_phy
2958	e , topsw, toplw, solsw, sollw, zero_v
2959	e , zero_v, zero_v, zero_v, ztsol
2960	e , d_h_vcol, d_qt, d_ec
2961	s , fs_bound, fq_bound )
2962	END IF
2963	c
2964	c
2965	c Calculer l'hydrologie de la surface
2966	c
2967	c CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
2968	c . agesno, ftsol,fqsurf,fsnow, ruis)
2969	c
2970
2971	c
2972	c Calculer le bilan du sol et la derive de temperature (couplage)
2973	c
2974	DO i = 1, klon
2975	c bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
2976	c a la demande de JLD
2977	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
2978	ENDDO
2979	c
2980	cmoddeblott(jan95)
2981	c Appeler le programme de parametrisation de l'orographie
2982	c a l'echelle sous-maille:
2983	c
2984	IF (ok_orodr) THEN
2985	c
2986	c selection des points pour lesquels le shema est actif:
2987	igwd=0
2988	DO i=1,klon
2989	itest(i)=0
2990	c IF ((zstd(i).gt.10.0)) THEN
2991	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
2992	itest(i)=1
2993	igwd=igwd+1
2994	idx(igwd)=i
2995	ENDIF
2996	ENDDO
2997	c igwdim=MAX(1,igwd)
2998	c
2999	CALL drag_noro(klon,klev,dtime,paprs,pplay,
3000	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
3001	e igwd,idx,itest,
3002	e t_seri, u_seri, v_seri,
3003	cIM 141004 s zulow, zvlow, zustr, zvstr,
3004	s zulow, zvlow, zustrdr, zvstrdr,
3005	s d_t_oro, d_u_oro, d_v_oro)
3006	c
3007	c ajout des tendances
3008	DO k = 1, klev
3009	DO i = 1, klon
3010	t_seri(i,k) = t_seri(i,k) + d_t_oro(i,k)
3011	u_seri(i,k) = u_seri(i,k) + d_u_oro(i,k)
3012	v_seri(i,k) = v_seri(i,k) + d_v_oro(i,k)
3013	ENDDO
3014	ENDDO
3015	c
3016	ENDIF ! fin de test sur ok_orodr
3017	c
3018	if (mydebug) then
3019	call writefield_phy('u_seri',u_seri,llm)
3020	call writefield_phy('v_seri',v_seri,llm)
3021	call writefield_phy('t_seri',t_seri,llm)
3022	call writefield_phy('q_seri',q_seri,llm)
3023	endif
3024
3025	IF (ok_orolf) THEN
3026	c
3027	c selection des points pour lesquels le shema est actif:
3028	igwd=0
3029	DO i=1,klon
3030	itest(i)=0
3031	IF ((zpic(i)-zmea(i)).GT.100.) THEN
3032	itest(i)=1
3033	igwd=igwd+1
3034	idx(igwd)=i
3035	ENDIF
3036	ENDDO
3037	c igwdim=MAX(1,igwd)
3038	c
3039	CALL lift_noro(klon,klev,dtime,paprs,pplay,
3040	e rlat,zmea,zstd,zpic,
3041	e itest,
3042	e t_seri, u_seri, v_seri,
3043	s zulow, zvlow, zustrli, zvstrli,
3044	s d_t_lif, d_u_lif, d_v_lif)
3045	c
3046	c ajout des tendances
3047	DO k = 1, klev
3048	DO i = 1, klon
3049	t_seri(i,k) = t_seri(i,k) + d_t_lif(i,k)
3050	u_seri(i,k) = u_seri(i,k) + d_u_lif(i,k)
3051	v_seri(i,k) = v_seri(i,k) + d_v_lif(i,k)
3052	ENDDO
3053	ENDDO
3054	c
3055	ENDIF ! fin de test sur ok_orolf
3056	c
3057	cIM cf. FLott BEG
3058	C STRESS NECESSAIRES: TOUTE LA PHYSIQUE
3059
3060	if (mydebug) then
3061	call writefield_phy('u_seri',u_seri,llm)
3062	call writefield_phy('v_seri',v_seri,llm)
3063	call writefield_phy('t_seri',t_seri,llm)
3064	call writefield_phy('q_seri',q_seri,llm)
3065	endif
3066
3067	DO i = 1, klon
3068	zustrph(i)=0.
3069	zvstrph(i)=0.
3070	ENDDO
3071	DO k = 1, klev
3072	DO i = 1, klon
3073	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
3074	c (paprs(i,k)-paprs(i,k+1))/rg
3075	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
3076	c (paprs(i,k)-paprs(i,k+1))/rg
3077	ENDDO
3078	ENDDO
3079	c
3080	cIM calcul composantes axiales du moment angulaire et couple des montagnes
3081	c
3082	IF (is_sequential) THEN
3083
3084	CALL aaam_bud (27,klon,klev,rjourvrai,gmtime,
3085	C ra,rg,romega,
3086	C rlat,rlon,pphis,
3087	C zustrdr,zustrli,zustrph,
3088	C zvstrdr,zvstrli,zvstrph,
3089	C paprs,u,v,
3090	C aam, torsfc)
3091	ENDIF
3092	cIM cf. FLott END
3093	cIM
3094	IF (ip_ebil_phy.ge.2) THEN
3095	ztit='after orography'
3096	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
3097	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3098	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3099	END IF
3100	c
3101	c
3102	cAA
3103	cAA Installation de l'interface online-offline pour traceurs
3104	cAA
3105	c====================================================================
3106	c Calcul des tendances traceurs
3107	c====================================================================
3108	C
3109	call phytrac ( rnpb,
3110	I itap,
3111	I julien,
3112	I gmtime,
3113	I debut,
3114	I lafin,
3115	I nqmax-2,
3116	I nlon,
3117	I nlev,
3118	I dtime,
3119	I u,
3120	I v,
3121	I t,
3122	I paprs,
3123	I pplay,
3124	I pmfu,
3125	I pmfd,
3126	I pen_u,
3127	I pde_u,
3128	I pen_d,
3129	I pde_d,
3130	I ycoefh,
3131	I fm_therm,
3132	I entr_therm,
3133	I yu1,
3134	I yv1,
3135	I ftsol,
3136	I pctsrf,
3137	I rlat,
3138	I frac_impa,
3139	I frac_nucl,
3140	I rlon,
3141	I presnivs,
3142	I pphis,
3143	I pphi,
3144	I albsol,
3145	I qx(1,1,1),
3146	I rhcl,
3147	I cldfra,
3148	I rneb,
3149	I diafra,
3150	I cldliq,
3151	I itop_con,
3152	I ibas_con,
3153	I pmflxr,
3154	I pmflxs,
3155	I prfl,
3156	I psfl,
3157	I da,
3158	I phi,
3159	I mp,
3160	I upwd,
3161	I dnwd,
3162	#ifdef INCA
3163	I flxmass_w,
3164	#endif
3165	O tr_seri)
3166
3167	IF (offline) THEN
3168
3169	print*,'Attention on met a 0 les thermiques pour phystoke'
3170	call phystokenc (
3171	I nlon,nlev,pdtphys,rlon,rlat,
3172	I t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
3173	I fm_therm,entr_therm,
3174	I ycoefh,yu1,yv1,ftsol,pctsrf,
3175	I frac_impa, frac_nucl,
3176	I pphis,airephy,dtime,itap)
3177
3178
3179	ENDIF
3180
3181	c
3182	c Calculer le transport de l'eau et de l'energie (diagnostique)
3183	c
3184	CALL transp (paprs,zxtsol,
3185	e t_seri, q_seri, u_seri, v_seri, zphi,
3186	s ve, vq, ue, uq)
3187	c
3188	cIM global posePB BEG
3189	IF(1.EQ.0) THEN
3190	c
3191	CALL transp_lay (paprs,zxtsol,
3192	e t_seri, q_seri, u_seri, v_seri, zphi,
3193	s ve_lay, vq_lay, ue_lay, uq_lay)
3194	c
3195	ENDIF !(1.EQ.0) THEN
3196	cIM global posePB END
3197	c Accumuler les variables a stocker dans les fichiers histoire:
3198	c
3199	c+jld ec_conser
3200	DO k = 1, klev
3201	DO i = 1, klon
3202	ZRCPD = RCPD(1.0+RVTMP2q_seri(i,k))
3203	d_t_ec(i,k)=0.5/ZRCPD
3204	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
3205	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
3206	d_t_ec(i,k) = d_t_ec(i,k)/dtime
3207	END DO
3208	END DO
3209	c-jld ec_conser
3210	cIM
3211	IF (ip_ebil_phy.ge.1) THEN
3212	ztit='after physic'
3213	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
3214	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3215	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3216	C Comme les tendances de la physique sont ajoute dans la dynamique,
3217	C on devrait avoir que la variation d'entalpie par la dynamique
3218	C est egale a la variation de la physique au pas de temps precedent.
3219	C Donc la somme de ces 2 variations devrait etre nulle.
3220	call diagphy(airephy,ztit,ip_ebil_phy
3221	e , topsw, toplw, solsw, sollw, sens
3222	e , evap, rain_fall, snow_fall, ztsol
3223	e , d_h_vcol, d_qt, d_ec
3224	s , fs_bound, fq_bound )
3225	C
3226	d_h_vcol_phy=d_h_vcol
3227	C
3228	END IF
3229	C
3230	c=======================================================================
3231	c SORTIES
3232	c=======================================================================
3233
3234	cIM Interpolation sur les niveaux de pression du NMC
3235	c -------------------------------------------------
3236	c
3237	#include "calcul_STDlev.h"
3238	c
3239	c slp sea level pressure
3240	slp(:) = paprs(:,1)exp(pphis(:)/(RDt_seri(:,1)))
3241	c
3242	ccc prw = eau precipitable
3243	DO i = 1, klon
3244	prw(i) = 0.
3245	DO k = 1, klev
3246	prw(i) = prw(i) +
3247	. q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
3248	ENDDO
3249	ENDDO
3250	c
3251	cIM initialisation + calculs divers diag AMIP2
3252	c
3253	#include "calcul_divers.h"
3254	c
3255	#ifdef INCA
3256	call VTe(VTphysiq)
3257	call VTb(VTinca)
3258	#ifdef INCAINFO
3259	WRITE(lunout,*)'Appel CHEMHOOK_END ...'
3260	#endif
3261	CALL chemhook_end (calday,
3262	$ dtime,
3263	$ pplay,
3264	$ t_seri,
3265	$ tr_seri,
3266	$ nbtr,
3267	$ paprs,
3268	$ q_seri,
3269	$ annee_ref,
3270	$ day_ini,
3271	$ airephy,
3272	#ifdef INCA_AER
3273	$ xjour,
3274	$ pphi,
3275	$ pphis,
3276	$ zx_rh)
3277	#else
3278	$ xjour)
3279	#endif
3280	#ifdef INCAINFO
3281	WRITE(lunout,*)'OK.'
3282	#endif
3283	call VTe(VTinca)
3284	call VTb(VTphysiq)
3285	#endif
3286
3287	c=============================================================
3288	c
3289	c Convertir les incrementations en tendances
3290	c
3291	if (mydebug) then
3292	call writefield_phy('u_seri',u_seri,llm)
3293	call writefield_phy('v_seri',v_seri,llm)
3294	call writefield_phy('t_seri',t_seri,llm)
3295	call writefield_phy('q_seri',q_seri,llm)
3296	endif
3297
3298	DO k = 1, klev
3299	DO i = 1, klon
3300	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
3301	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
3302	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
3303	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
3304	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
3305	ENDDO
3306	ENDDO
3307	c
3308	IF (nqmax.GE.3) THEN
3309	DO iq = 3, nqmax
3310	DO k = 1, klev
3311	DO i = 1, klon
3312	d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
3313	ENDDO
3314	ENDDO
3315	ENDDO
3316	ENDIF
3317	c
3318	cIM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
3319	cIM global posePB#include "write_bilKP_ins.h"
3320	cIM global posePB#include "write_bilKP_ave.h"
3321	c
3322	c Sauvegarder les valeurs de t et q a la fin de la physique:
3323	c
3324	DO k = 1, klev
3325	DO i = 1, klon
3326	t_ancien(i,k) = t_seri(i,k)
3327	q_ancien(i,k) = q_seri(i,k)
3328	ENDDO
3329	ENDDO
3330	c
3331	c 22.03.04 BEG
3332	c=============================================================
3333	c Ecriture des sorties
3334	c=============================================================
3335	#ifdef CPP_IOIPSL
3336
3337	c Recupere des varibles calcule dans differents modules
3338	c pour ecriture dans histxxx.nc
3339
3340	! Get some variables from module mod_fonte_neige
3341	CALL fonte_neige_get_vars(pctsrf,
3342	. zxfqcalving, zxfqfonte, zxffonte)
3343
3344	IF (ocean == 'slab') THEN
3345	! Get some variables from module oceanslab
3346	CALL ocean_slab_get_vars(tslab, seaice, fluxo, fluxg)
3347	ELSEIF (ocean == 'couple') THEN
3348	! Get some variables from module oceancpl
3349	CALL ocean_cpl_get_vars(fluxo, fluxg)
3350	ELSE
3351	! Get some variables from module oceanforced
3352	CALL ocean_forced_get_vars(fluxo, fluxg)
3353	ENDIF
3354
3355	#ifdef histhf
3356	#include "write_histhf.h"
3357	#endif
3358
3359	#ifdef histday
3360	#include "write_histday.h"
3361	#endif
3362
3363	#ifdef histmth
3364	#include "write_histmth.h"
3365	#endif
3366
3367	#ifdef histins
3368	#include "write_histins.h"
3369	#endif
3370
3371	#ifdef histISCCP
3372	#include "write_histISCCP.h"
3373	#endif
3374
3375	#ifdef histmthNMC
3376	#include "write_histmthNMC.h"
3377	#endif
3378
3379	#include "write_histday_seri.h"
3380
3381	#include "write_paramLMDZ_phy.h"
3382
3383	#endif
3384
3385	c 22.03.04 END
3386	c
3387	c====================================================================
3388	c Si c'est la fin, il faut conserver l'etat de redemarrage
3389	c====================================================================
3390	c
3391
3392
3393	IF (lafin) THEN
3394	itau_phy = itau_phy + itap
3395	CALL phyredem ("restartphy.nc",dtime,radpas,ocean,
3396	. rlat, rlon, pctsrf, ftsol,
3397	. falbe,falblw, rain_fall,
3398	. snow_fall,
3399	. solsw, sollw,
3400	. radsol,
3401	. zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,
3402	. t_ancien, q_ancien, rnebcon, ratqs, clwcon,
3403	. pbl_tke)
3404	ENDIF
3405
3406
3407	RETURN
3408	END
3409	FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
3410	IMPLICIT none
3411	c
3412	c Calculer et imprimer l'eau totale. A utiliser pour verifier
3413	c la conservation de l'eau
3414	c
3415	#include "YOMCST.h"
3416	INTEGER klon,klev
3417	REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
3418	REAL aire(klon)
3419	REAL qtotal, zx, qcheck
3420	INTEGER i, k
3421	c
3422	zx = 0.0
3423	DO i = 1, klon
3424	zx = zx + aire(i)
3425	ENDDO
3426	qtotal = 0.0
3427	DO k = 1, klev
3428	DO i = 1, klon
3429	qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i)
3430	. *(paprs(i,k)-paprs(i,k+1))/RG
3431	ENDDO
3432	ENDDO
3433	c
3434	qcheck = qtotal/zx
3435	c
3436	RETURN
3437	END
3438	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
3439	IMPLICIT none
3440	c
3441	c Tranformer une variable de la grille physique a
3442	c la grille d'ecriture
3443	c
3444	INTEGER nfield,nlon,iim,jjmp1, jjm
3445	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
3446	c
3447	INTEGER i, n, ig
3448	c
3449	jjm = jjmp1 - 1
3450	DO n = 1, nfield
3451	DO i=1,iim
3452	ecrit(i,n) = fi(1,n)
3453	ecrit(i+jjm*iim,n) = fi(nlon,n)
3454	ENDDO
3455	DO ig = 1, nlon - 2
3456	ecrit(iim+ig,n) = fi(1+ig,n)
3457	ENDDO
3458	ENDDO
3459	RETURN
3460	END

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

Download in other formats:

Original Format