Context Navigation

physiq.F @ 826

Last change on this file since 826 was 814, checked in by Laurent Fairhead, 18 years ago
Rajout de la physique utilisant les thermiques FH LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 106.7 KB

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

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

Download in other formats:

Original Format