Context Navigation

physiq.F @ 904

Last change on this file since 904 was 904, checked in by Laurent Fairhead, 16 years ago
Debut de la reecriture de la routine physi.F. On commence par une extension de hgardfou dans le module add_phys_tend.F90 qui verifie que les tendances physiques calculees dans les differentes parametrisations sont dans des limites tolerables avant de les ajouter aux variables d'etat de la dynamique FH LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 128.7 KB

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

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

Download in other formats:

Original Format