Context Navigation

physiq.F @ 687

Last change on this file since 687 was 687, checked in by lmdzadmin, 19 years ago
Ajout diagnostiques rh, tpot a 2m et wbilo + lecture flag ip_ebil_phy dans conf_phys.F90 IM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 102.0 KB

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

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

Download in other formats:

Original Format