Context Navigation

physiq.F @ 5379

Last change on this file since 5379 was 694, checked in by lmdzadmin, 19 years ago
Correction bogue : debordement de memoire si klev < nlevSTD(=17) PLV/IM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 102.1 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	#define histmthNMC
45	#define histISCCP
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	#ifdef histmthNMC
1155	REAL zx_tmp_NC(iim,jjmp1,nlevSTD)
1156	REAL zx_tmp_fiNC(klon,nlevSTD)
1157	#endif
1158	REAL*8 zx_tmp2_fi3d(klon,klev) ! variable temporaire pour champs 3D
1159	REAL zx_tmp_2d(iim,jjmp1), zx_tmp_3d(iim,jjmp1,klev)
1160	REAL zx_lon(iim,jjmp1), zx_lat(iim,jjmp1)
1161	c
1162	INTEGER nid_day, nid_mth, nid_ins, nid_nmc, nid_day_seri
1163	INTEGER nid_ctesGCM
1164	SAVE nid_day, nid_mth, nid_ins, nid_nmc, nid_day_seri
1165	SAVE nid_ctesGCM
1166	c
1167	cIM 280405 BEG
1168	INTEGER nid_bilKPins, nid_bilKPave
1169	SAVE nid_bilKPins, nid_bilKPave
1170	c
1171	REAL ve_lay(klon,klev) ! transport meri. de l'energie a chaque niveau vert.
1172	REAL vq_lay(klon,klev) ! transport meri. de l'eau a chaque niveau vert.
1173	REAL ue_lay(klon,klev) ! transport zonal de l'energie a chaque niveau vert.
1174	REAL uq_lay(klon,klev) ! transport zonal de l'eau a chaque niveau vert.
1175	c
1176	cIM 280405 END
1177	c
1178	INTEGER nhori, nvert, nvert1, nvert3
1179	REAL zsto, zsto1, zsto2
1180	REAL zstophy, zstorad, zstohf, zstoday, zstomth, zout
1181	REAL zcals(napisccp), zcalh(napisccp), zoutj(napisccp)
1182	REAL zout_isccp(napisccp)
1183	SAVE zcals, zcalh, zoutj, zout_isccp
1184
1185	real zjulian
1186	save zjulian
1187
1188	character*20 modname
1189	character*80 abort_message
1190	logical ok_sync
1191	real date0
1192	integer idayref
1193
1194	C essai writephys
1195	integer fid_day, fid_mth, fid_ins
1196	parameter (fid_ins = 1, fid_day = 2, fid_mth = 3)
1197	integer prof2d_on, prof3d_on, prof2d_av, prof3d_av
1198	parameter (prof2d_on = 1, prof3d_on = 2,
1199	. prof2d_av = 3, prof3d_av = 4)
1200	character*30 nom_fichier
1201	character*10 varname
1202	character*40 vartitle
1203	character*20 varunits
1204	C Variables liees au bilan d'energie et d'enthalpi
1205	REAL ztsol(klon)
1206	REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
1207	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
1208	SAVE h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
1209	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
1210	REAL d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
1211	REAL d_h_vcol_phy
1212	REAL fs_bound, fq_bound
1213	SAVE d_h_vcol_phy
1214	REAL zero_v(klon)
1215	CHARACTER*15 ztit
1216
1217	c+jld ec_conser
1218	REAL d_t_ec(klon,klev) ! tendance du a la conersion Ec -> E thermique
1219	REAL ZRCPD
1220	c-jld ec_conser
1221	cIM: t2m, q2m, u10m, v10m et t2mincels, t2maxcels
1222	REAL t2m(klon,nbsrf), q2m(klon,nbsrf) !temperature, humidite a 2m
1223	REAL u10m(klon,nbsrf), v10m(klon,nbsrf) !vents a 10m
1224	REAL zt2m(klon), zq2m(klon) !temp., hum. 2m moyenne s/ 1 maille
1225	REAL zu10m(klon), zv10m(klon) !vents a 10m moyennes s/1 maille
1226	CHARACTER*40 t2mincels, t2maxcels !t2m min., t2m max
1227	CHARACTER*40 tinst, tave, typeval
1228	cjq Aerosol effects (Johannes Quaas, 27/11/2003)
1229	REAL sulfate(klon, klev) ! SO4 aerosol concentration [ug/m3]
1230	REAL sulfate_pi(klon, klev) ! SO4 aerosol concentration [ug/m3] (pre-industrial value)
1231	SAVE sulfate_pi
1232
1233	REAL cldtaupi(klon,klev) ! Cloud optical thickness for pre-industrial (pi) aerosols
1234
1235	REAL re(klon, klev) ! Cloud droplet effective radius
1236	REAL fl(klon, klev) ! denominator of re
1237
1238	REAL re_top(klon), fl_top(klon) ! CDR at top of liquid water clouds
1239
1240	! Aerosol optical properties
1241	REAL tau_ae(klon,klev,2), piz_ae(klon,klev,2)
1242	REAL cg_ae(klon,klev,2)
1243
1244	REAL topswad(klon), solswad(klon) ! Aerosol direct effect.
1245	! ok_ade=T -ADE=topswad-topsw
1246
1247	REAL topswai(klon), solswai(klon) ! Aerosol indirect effect.
1248	! ok_aie=T ->
1249	! ok_ade=T -AIE=topswai-topswad
1250	! ok_ade=F -AIE=topswai-topsw
1251
1252	REAL aerindex(klon) ! POLDER aerosol index
1253
1254	! Parameters
1255	LOGICAL ok_ade, ok_aie ! Apply aerosol (in)direct effects or not
1256	REAL bl95_b0, bl95_b1 ! Parameter in Boucher and Lohmann (1995)
1257	cym
1258	SAVE ok_ade, ok_aie, bl95_b0, bl95_b1
1259	cym
1260	c Anne
1261	SAVE u10m
1262	SAVE v10m
1263	SAVE t2m
1264	SAVE q2m
1265	SAVE ffonte
1266	SAVE fqcalving
1267	SAVE piz_ae
1268	SAVE tau_ae
1269	SAVE cg_ae
1270	SAVE rain_con
1271	SAVE snow_con
1272	SAVE topswai
1273	SAVE topswad
1274	SAVE solswai
1275	SAVE solswad
1276	SAVE d_u_con
1277	SAVE d_v_con
1278	SAVE rnebcon0
1279	SAVE clwcon0
1280	SAVE paire_ter
1281	c SAVE nhistoW
1282	c SAVE histoW
1283	c SAVE anne 20/09/2005
1284	SAVE pblh
1285	SAVE plcl
1286	SAVE capCL
1287	SAVE oliqCL
1288	SAVE cteiCL
1289	SAVE pblt
1290	SAVE therm
1291	SAVE trmb1
1292	SAVE trmb2
1293	SAVE trmb3
1294
1295	c fin Anne
1296	cjq-end
1297	c
1298	c Declaration des constantes et des fonctions thermodynamiques
1299	c
1300	#include "YOMCST.h"
1301	#include "YOETHF.h"
1302	#include "FCTTRE.h"
1303	cIM 100106 BEG : pouvoir sortir les ctes de la physique
1304	#include "conema3.h"
1305	#include "fisrtilp.h"
1306	#include "nuage.h"
1307	#include "compbl.h"
1308	cIM 100106 END : pouvoir sortir les ctes de la physique
1309	c
1310	c======================================================================
1311	modname = 'physiq'
1312	cIM
1313	IF (ip_ebil_phy.ge.1) THEN
1314	DO i=1,klon
1315	zero_v(i)=0.
1316	END DO
1317	END IF
1318	ok_sync=.TRUE.
1319	IF (nqmax .LT. 2) THEN
1320	abort_message = 'eaux vapeur et liquide sont indispensables'
1321	CALL abort_gcm (modname,abort_message,1)
1322	ENDIF
1323	IF (debut) THEN
1324	CALL suphec ! initialiser constantes et parametres phys.
1325	ENDIF
1326
1327
1328	c======================================================================
1329	xjour = rjourvrai
1330	c
1331	c Si c'est le debut, il faut initialiser plusieurs choses
1332	c ********
1333	c
1334	IF (debut) THEN
1335	C
1336	!rv
1337	u10m(:,:)=0.
1338	v10m(:,:)=0.
1339	t2m(:,:)=0.
1340	q2m(:,:)=0.
1341	ffonte(:,:)=0.
1342	fqcalving(:,:)=0.
1343	piz_ae(:,:,:)=0.
1344	tau_ae(:,:,:)=0.
1345	cg_ae(:,:,:)=0.
1346	rain_con(:)=0.
1347	snow_con(:)=0.
1348	bl95_b0=0.
1349	bl95_b1=0.
1350	topswai(:)=0.
1351	topswad(:)=0.
1352	solswai(:)=0.
1353	solswad(:)=0.
1354	!rv
1355	! anne
1356	d_u_con(:,:) = 0.0
1357	d_v_con(:,:) = 0.0
1358	rnebcon0(:,:) = 0.0
1359	clwcon0(:,:) = 0.0
1360	rnebcon(:,:) = 0.0
1361	clwcon(:,:) = 0.0
1362	paire_ter(:) = 0.0
1363	c nhistoW(:,:,:,:) = 0.0
1364	c histoW(:,:,:,:) = 0.0
1365	! fin anne
1366	! Anne 12/09/2005
1367
1368	pblh(:,:) =0. ! Hauteur de couche limite
1369	plcl(:,:) =0. ! Niveau de condensation de la CLA
1370	capCL(:,:) =0. ! CAPE de couche limite
1371	oliqCL(:,:) =0. ! eau_liqu integree de couche limite
1372	cteiCL(:,:) =0. ! cloud top instab. crit. couche limite
1373	pblt(:,:) =0. ! T a la Hauteur de couche limite
1374	therm(:,:) =0.
1375	trmb1(:,:) =0. ! deep_cape
1376	trmb2(:,:) =0. ! inhibition
1377	trmb3(:,:) =0. ! Point Omega
1378	! fin Anne
1379
1380	cym
1381	wfbils(:,:)=0
1382	cym
1383	cIM
1384	IF (ip_ebil_phy.ge.1) d_h_vcol_phy=0.
1385	c
1386	c appel a la lecture du run.def physique
1387	c
1388	call conf_phys(ocean, ok_veget, ok_journe, ok_mensuel,
1389	. ok_instan, fact_cldcon, facttemps,ok_newmicro,
1390	cIM . iflag_cldcon,ratqsbas,ratqshaut, if_ebil,
1391	. iflag_cldcon,ratqsbas,ratqshaut,
1392	. ok_ade, ok_aie,
1393	. bl95_b0, bl95_b1,
1394	. iflag_thermals,nsplit_thermals)
1395
1396	c
1397	c
1398	c Initialiser les compteurs:
1399	c
1400
1401	frugs = 0.
1402	itap = 0
1403	itaprad = 0
1404	CALL phyetat0 ("startphy.nc",dtime,co2_ppm_etat0,solaire_etat0,
1405	. rlat,rlon,pctsrf, ftsol,ftsoil,
1406	cIM "slab" ocean
1407	. ocean, tslab,seaice,
1408	. fqsurf,qsol,fsnow,
1409	cIM 220306 . falbe, falblw, fevap, rain_fall,snow_fall,solsw, sollwdown,
1410	. falbe, falblw, fevap, rain_fall,snow_fall,solsw, sollw,
1411	. dlw,radsol,frugs,agesno,clesphy0,
1412	. zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,tabcntr0,
1413	. t_ancien, q_ancien, ancien_ok, rnebcon, ratqs,clwcon,
1414	. run_off_lic_0)
1415
1416	c ATTENTION : il faudra a terme relire q2 dans l'etat initial
1417	q2(:,:,:)=1.e-8
1418	c
1419	radpas = NINT( 86400./dtime/nbapp_rad)
1420	c
1421	C on remet le calendrier a zero
1422	c
1423	IF (raz_date .eq. 1) THEN
1424	itau_phy = 0
1425	ENDIF
1426
1427	cIM cf. AM 081204 BEG
1428	PRINT*,'cycle_diurne3 =',cycle_diurne
1429	cIM cf. AM 081204 END
1430	c
1431	IF(ocean.NE.'force ') THEN
1432	ok_ocean=.TRUE.
1433	ENDIF
1434	c
1435	CALL printflag( tabcntr0,radpas,ok_ocean,ok_oasis ,ok_journe,
1436	, ok_instan, ok_region )
1437	c
1438	IF (ABS(dtime-pdtphys).GT.0.001) THEN
1439	WRITE(lunout,*) 'Pas physique n est pas correct',dtime,
1440	. pdtphys
1441	abort_message='Pas physique n est pas correct '
1442	call abort_gcm(modname,abort_message,1)
1443	ENDIF
1444	IF (nlon .NE. klon) THEN
1445	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,
1446	. klon
1447	abort_message='nlon et klon ne sont pas coherents'
1448	call abort_gcm(modname,abort_message,1)
1449	ENDIF
1450	IF (nlev .NE. klev) THEN
1451	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,
1452	. klev
1453	abort_message='nlev et klev ne sont pas coherents'
1454	call abort_gcm(modname,abort_message,1)
1455	ENDIF
1456	c
1457	IF (dtime*FLOAT(radpas).GT.21600..AND.cycle_diurne) THEN
1458	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
1459	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
1460	abort_message='Nbre d appels au rayonnement insuffisant'
1461	call abort_gcm(modname,abort_message,1)
1462	ENDIF
1463	WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
1464	WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=",
1465	. ok_cvl
1466	c
1467	cKE43
1468	c Initialisation pour la convection de K.E. (sb):
1469	IF (iflag_con.GE.3) THEN
1470
1471	WRITE(lunout,)"** Convection de Kerry Emanuel 4.3 "
1472	WRITE(lunout,*)
1473	. "On va utiliser le melange convectif des traceurs qui"
1474	WRITE(lunout,*)"est calcule dans convect4.3"
1475	WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
1476
1477	DO i = 1, klon
1478	ema_cbmf(i) = 0.
1479	ema_pcb(i) = 0.
1480	ema_pct(i) = 0.
1481	ema_workcbmf(i) = 0.
1482	ENDDO
1483	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
1484	DO i = 1, klon
1485	ibas_con(i) = 1
1486	itop_con(i) = 1
1487	ENDDO
1488	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
1489
1490	ENDIF
1491
1492	c34EK
1493	IF (ok_orodr) THEN
1494	DO i=1,klon
1495	rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1496	ENDDO
1497	CALL SUGWD(klon,klev,paprs,pplay)
1498	DO i=1,klon
1499	zuthe(i)=0.
1500	zvthe(i)=0.
1501	if(zstd(i).gt.10.)then
1502	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
1503	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
1504	endif
1505	ENDDO
1506	ENDIF
1507	c
1508	c
1509	lmt_pas = NINT(86400./dtime * 1.0) ! tous les jours
1510	WRITE(lunout,*)'La frequence de lecture surface est de ',
1511	. lmt_pas
1512	c
1513	cIM200505 ecrit_mth = NINT(86400./dtime *ecritphy) ! tous les ecritphy jours
1514	c IF (ok_mensuel) THEN
1515	c WRITE(lunout,*)'La frequence de sortie mensuelle est de ',
1516	c . ecrit_mth
1517	c ENDIF
1518	c ecrit_day = NINT(86400./dtime *1.0) ! tous les jours
1519	c IF (ok_journe) THEN
1520	c WRITE(lunout,*)'La frequence de sortie journaliere est de ',
1521	c . ecrit_day
1522	c ENDIF
1523	cIM 130904 BEG
1524	cIM 080205 ecrit_hf = 86400./dtime *0.25 ! toutes les 6h
1525	cIM 170305
1526	c ecrit_hf = 86400./dtime/12. ! toutes les 2h
1527	cIM 230305
1528	cIM200505 ecrit_hf = 86400./dtime *0.25 ! toutes les 6h
1529	c
1530	cIM200505 ecrit_hf2mth = ecrit_day/ecrit_hf*30
1531	c
1532	cIM200505 IF (ok_journe) THEN
1533	cIM200505 WRITE(lunout,*)'La frequence de sortie hf est de ',
1534	cIM200505 . ecrit_hf
1535	cIM200505 ENDIF
1536	cIM 130904 END
1537	ccc ecrit_ins = NINT(86400./dtime *0.5) ! 2 fois par jour
1538	ccc ecrit_ins = NINT(86400./dtime *0.25) ! 4 fois par jour
1539	c ecrit_ins = NINT(86400./dtime/48.) ! a chaque pas de temps ==> PB. dans time_counter pour 1mois
1540	c ecrit_ins = NINT(86400./dtime/12.) ! toutes les deux heures
1541	cIM200505 ecrit_ins = NINT(86400./dtime/8.) ! toutes les trois heures
1542	cIM200505 IF (ok_instan) THEN
1543	cIM200505 WRITE(lunout,*)'La frequence de sortie instant. est de ',
1544	cIM200505 . ecrit_ins
1545	cIM200505 ENDIF
1546	cIM200505 ecrit_reg = NINT(86400./dtime *0.25) ! 4 fois par jour
1547	cIM200505 IF (ok_region) THEN
1548	cIM200505 WRITE(lunout,*)'La frequence de sortie region est de ',
1549	cIM200505 . ecrit_reg
1550	cIM200505 ENDIF
1551	cIM 030306 BEG
1552	cIM ecrit_hf2mth = nombre de pas de temps de calcul de hf par mois apres lequel on ecrit
1553	cIM : ne pas modifier ecrit_hf2mth
1554	c
1555	ecrit_hf2mth = 30*1/ecrit_hf
1556	c ecrit_ins, ecrit_tra en secondes, chaque pas de temps de la physique
1557	ecrit_ins = dtime
1558	ecrit_tra = dtime
1559	cIM on passe les frequences de jours en secondes : ecrit_ins, ecrit_hf, ecrit_day, ecrit_mth, ecrit_tra, ecrit_reg
1560	ecrit_hf = ecrit_hf * un_jour
1561	ecrit_day = ecrit_day * un_jour
1562	ecrit_mth = ecrit_mth * un_jour
1563	ecrit_reg = ecrit_reg * un_jour
1564	cIM 030306 END
1565	c
1566	c Initialiser le couplage si necessaire
1567	c
1568	npas = 0
1569	nexca = 0
1570	if (ocean == 'couple') then
1571	npas = itaufin/ iphysiq
1572	nexca = 86400 / dtime
1573	write(lunout,*)' ##### Ocean couple #####'
1574	write(lunout,*)' Valeurs des pas de temps'
1575	write(lunout,*)' npas = ', npas
1576	write(lunout,*)' nexca = ', nexca
1577	endif
1578	c
1579	capemaxcels = 't_max(X)'
1580	t2mincels = 't_min(X)'
1581	t2maxcels = 't_max(X)'
1582	tinst = 'inst(X)'
1583	tave = 'ave(X)'
1584	cIM cf. AM 081204 BEG
1585	write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
1586	cIM cf. AM 081204 END
1587	c
1588	c=============================================================
1589	c Initialisation des sorties
1590	c=============================================================
1591
1592	#ifdef CPP_IOIPSL
1593
1594	#ifdef histhf
1595	#include "ini_histhf.h"
1596	#endif
1597
1598	#ifdef histday
1599	#include "ini_histday.h"
1600	cIM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
1601	c#include "ini_bilKP_ins.h"
1602	c#include "ini_bilKP_ave.h"
1603	#endif
1604
1605	#ifdef histmth
1606	#include "ini_histmth.h"
1607	#endif
1608
1609	#ifdef histins
1610	#include "ini_histins.h"
1611	#endif
1612
1613	#ifdef histISCCP
1614	#include "ini_histISCCP.h"
1615	#endif
1616
1617	#undef histmthNMC
1618	#define histmthNMC
1619	#ifdef histmthNMC
1620	#include "ini_histmthNMC.h"
1621	#endif
1622
1623	#include "ini_histday_seri.h"
1624
1625	#include "ini_paramLMDZ_phy.h"
1626
1627	#endif
1628
1629	cXXXPB Positionner date0 pour initialisation de ORCHIDEE
1630	date0 = zjulian
1631	C date0 = day_ini
1632	WRITE(,) 'physiq date0 : ',date0
1633	c
1634	c
1635	c
1636	c Prescrire l'ozone dans l'atmosphere
1637	c
1638	c
1639	cc DO i = 1, klon
1640	cc DO k = 1, klev
1641	cc CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1642	cc ENDDO
1643	cc ENDDO
1644	c
1645	#ifdef INCA
1646	iii = MOD(NINT(xjour),360)
1647	calday = FLOAT(iii) + gmtime
1648	WRITE(lunout,*) 'initial time ', xjour, calday
1649	#ifdef INCAINFO
1650	WRITE(lunout,*) 'Appel CHEMINI ...'
1651	#endif
1652	CALL chemini( rpi,
1653	$ rg,
1654	$ ra,
1655	$ airephy,
1656	$ rlat,
1657	$ rlon,
1658	$ presnivs,
1659	$ calday,
1660	$ klon,
1661	$ nqmax,
1662	$ pdtphys,
1663	$ annee_ref,
1664	$ day_ini)
1665	#ifdef INCAINFO
1666	WRITE(lunout,*) 'OK.'
1667	#endif
1668	#endif
1669	c
1670	ENDIF
1671	c
1672	c ************** Fin de IF ( debut ) *************
1673	c
1674	c
1675	c Mettre a zero des variables de sortie (pour securite)
1676	c
1677	DO i = 1, klon
1678	d_ps(i) = 0.0
1679	ENDDO
1680	DO k = 1, klev
1681	DO i = 1, klon
1682	d_t(i,k) = 0.0
1683	d_u(i,k) = 0.0
1684	d_v(i,k) = 0.0
1685	ENDDO
1686	ENDDO
1687	DO iq = 1, nqmax
1688	DO k = 1, klev
1689	DO i = 1, klon
1690	d_qx(i,k,iq) = 0.0
1691	ENDDO
1692	ENDDO
1693	ENDDO
1694	da(:,:)=0.
1695	mp(:,:)=0.
1696	phi(:,:,:)=0.
1697	c
1698	c Ne pas affecter les valeurs entrees de u, v, h, et q
1699	c
1700	DO k = 1, klev
1701	DO i = 1, klon
1702	t_seri(i,k) = t(i,k)
1703	u_seri(i,k) = u(i,k)
1704	v_seri(i,k) = v(i,k)
1705	q_seri(i,k) = qx(i,k,ivap)
1706	ql_seri(i,k) = qx(i,k,iliq)
1707	qs_seri(i,k) = 0.
1708	ENDDO
1709	ENDDO
1710	IF (nqmax.GE.3) THEN
1711	DO iq = 3, nqmax
1712	DO k = 1, klev
1713	DO i = 1, klon
1714	tr_seri(i,k,iq-2) = qx(i,k,iq)
1715	ENDDO
1716	ENDDO
1717	ENDDO
1718	ELSE
1719	DO k = 1, klev
1720	DO i = 1, klon
1721	tr_seri(i,k,1) = 0.0
1722	ENDDO
1723	ENDDO
1724	ENDIF
1725	C
1726	DO i = 1, klon
1727	ztsol(i) = 0.
1728	ENDDO
1729	DO nsrf = 1, nbsrf
1730	DO i = 1, klon
1731	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1732	ENDDO
1733	ENDDO
1734	cIM
1735	IF (ip_ebil_phy.ge.1) THEN
1736	ztit='after dynamic'
1737	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
1738	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1739	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1740	C Comme les tendances de la physique sont ajoute dans la dynamique,
1741	C on devrait avoir que la variation d'entalpie par la dynamique
1742	C est egale a la variation de la physique au pas de temps precedent.
1743	C Donc la somme de ces 2 variations devrait etre nulle.
1744	call diagphy(airephy,ztit,ip_ebil_phy
1745	e , zero_v, zero_v, zero_v, zero_v, zero_v
1746	e , zero_v, zero_v, zero_v, ztsol
1747	e , d_h_vcol+d_h_vcol_phy, d_qt, 0.
1748	s , fs_bound, fq_bound )
1749	END IF
1750
1751	c Diagnostiquer la tendance dynamique
1752	c
1753	IF (ancien_ok) THEN
1754	DO k = 1, klev
1755	DO i = 1, klon
1756	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
1757	d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
1758	ENDDO
1759	ENDDO
1760	ELSE
1761	DO k = 1, klev
1762	DO i = 1, klon
1763	d_t_dyn(i,k) = 0.0
1764	d_q_dyn(i,k) = 0.0
1765	ENDDO
1766	ENDDO
1767	ancien_ok = .TRUE.
1768	ENDIF
1769	c
1770	c Ajouter le geopotentiel du sol:
1771	c
1772	DO k = 1, klev
1773	DO i = 1, klon
1774	zphi(i,k) = pphi(i,k) + pphis(i)
1775	ENDDO
1776	ENDDO
1777	c
1778	c Verifier les temperatures
1779	c
1780	cIM BEG
1781	IF (check) THEN
1782	amn=MIN(ftsol(1,is_ter),1000.)
1783	amx=MAX(ftsol(1,is_ter),-1000.)
1784	DO i=2, klon
1785	amn=MIN(ftsol(i,is_ter),amn)
1786	amx=MAX(ftsol(i,is_ter),amx)
1787	ENDDO
1788	c
1789	PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
1790	ENDIF !(check) THEN
1791	cIM END
1792	c
1793	CALL hgardfou(t_seri,ftsol,'debutphy')
1794	c
1795	cIM BEG
1796	IF (check) THEN
1797	amn=MIN(ftsol(1,is_ter),1000.)
1798	amx=MAX(ftsol(1,is_ter),-1000.)
1799	DO i=2, klon
1800	amn=MIN(ftsol(i,is_ter),amn)
1801	amx=MAX(ftsol(i,is_ter),amx)
1802	ENDDO
1803	c
1804	PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
1805	ENDIF !(check) THEN
1806	cIM END
1807	c
1808	c Incrementer le compteur de la physique
1809	c
1810	itap = itap + 1
1811	julien = MOD(NINT(xjour),360)
1812	if (julien .eq. 0) julien = 360
1813	c
1814	c Mettre en action les conditions aux limites (albedo, sst, etc.).
1815	c Prescrire l'ozone et calculer l'albedo sur l'ocean.
1816	c
1817	IF (MOD(itap-1,lmt_pas) .EQ. 0) THEN
1818	WRITE(lunout,*)' PHYS cond julien ',julien
1819	CALL ozonecm( FLOAT(julien), rlat, paprs, wo)
1820	ENDIF
1821	c
1822	c Re-evaporer l'eau liquide nuageuse
1823	c
1824	DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse
1825	DO i = 1, klon
1826	zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1827	c zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1828	zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1829	zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
1830	zb = MAX(0.0,ql_seri(i,k))
1831	za = - MAX(0.0,ql_seri(i,k))
1832	. * (zlvdcp(1.-zdelta)+zlsdcpzdelta)
1833	t_seri(i,k) = t_seri(i,k) + za
1834	q_seri(i,k) = q_seri(i,k) + zb
1835	ql_seri(i,k) = 0.0
1836	d_t_eva(i,k) = za
1837	d_q_eva(i,k) = zb
1838	ENDDO
1839	ENDDO
1840	cIM
1841	IF (ip_ebil_phy.ge.2) THEN
1842	ztit='after reevap'
1843	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,1,dtime
1844	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1845	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1846	call diagphy(airephy,ztit,ip_ebil_phy
1847	e , zero_v, zero_v, zero_v, zero_v, zero_v
1848	e , zero_v, zero_v, zero_v, ztsol
1849	e , d_h_vcol, d_qt, d_ec
1850	s , fs_bound, fq_bound )
1851	C
1852	END IF
1853	C
1854	c
1855	c Appeler la diffusion verticale (programme de couche limite)
1856	c
1857	DO i = 1, klon
1858	c if (.not. ok_veget) then
1859	c frugs(i,is_ter) = SQRT(frugs(i,is_ter)2+rugoro(i)2)
1860	c endif
1861	c frugs(i,is_lic) = rugoro(i)
1862	c frugs(i,is_oce) = rugmer(i)
1863	c frugs(i,is_sic) = 0.001
1864	zxrugs(i) = 0.0
1865	ENDDO
1866	DO nsrf = 1, nbsrf
1867	DO i = 1, klon
1868	c frugs(i,nsrf) = MAX(frugs(i,nsrf),0.001)
1869	frugs(i,nsrf) = MAX(frugs(i,nsrf),0.000015)
1870	ENDDO
1871	ENDDO
1872	DO nsrf = 1, nbsrf
1873	DO i = 1, klon
1874	zxrugs(i) = zxrugs(i) + frugs(i,nsrf)*pctsrf(i,nsrf)
1875	ENDDO
1876	ENDDO
1877	c
1878	C calculs necessaires au calcul de l'albedo dans l'interface
1879	c
1880	CALL orbite(FLOAT(julien),zlongi,dist)
1881	IF (cycle_diurne) THEN
1882	zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
1883	CALL zenang(zlongi,gmtime,zdtime,rlat,rlon,rmu0,fract)
1884	ELSE
1885	rmu0 = -999.999
1886	ENDIF
1887	c
1888	C Calcul de l'abedo moyen par maille
1889	albsol(:)=0.
1890	albsollw(:)=0.
1891	DO nsrf = 1, nbsrf
1892	DO i = 1, klon
1893	albsol(i) = albsol(i) + falbe(i,nsrf) * pctsrf(i,nsrf)
1894	albsollw(i) = albsollw(i) + falblw(i,nsrf) * pctsrf(i,nsrf)
1895	ENDDO
1896	ENDDO
1897	C
1898	C Repartition sous maille des flux LW et SW
1899	C Modif OM+PASB+JLD
1900	C Repartition du longwave par sous-surface linearisee
1901	Cn
1902
1903	DO nsrf = 1, nbsrf
1904	DO i = 1, klon
1905	c$$$ fsollw(i,nsrf) = sollwdown(i) - RSIGMAftsol(i,nsrf)*4
1906	c$$$ fsollw(i,nsrf) = sollw(i)
1907	fsollw(i,nsrf) = sollw(i)
1908	$ + 4.0RSIGMAztsol(i)*3 (ztsol(i)-ftsol(i,nsrf))
1909	fsolsw(i,nsrf) = solsw(i)*(1.-falbe(i,nsrf))/(1.-albsol(i))
1910	ENDDO
1911	ENDDO
1912
1913	fder = dlw
1914
1915	IF (check) THEN
1916	amn=MIN(tslab(1),1000.)
1917	amx=MAX(tslab(1),-1000.)
1918	DO i=2, klon
1919	amn=MIN(tslab(i),amn)
1920	amx=MAX(tslab(i),amx)
1921	ENDDO
1922	c
1923	PRINT*,' debut avant clqh min max tslab',amn,amx
1924	ENDIF !(check) THEN
1925	c
1926	CALL clmain(dtime,itap,date0,pctsrf,pctsrf_new,
1927	e t_seri,q_seri,u_seri,v_seri,
1928	e julien, rmu0, co2_ppm,
1929	e ok_veget, ocean, npas, nexca, ftsol,
1930	$ soil_model,cdmmax, cdhmax,
1931	$ ksta, ksta_ter, ok_kzmin, ftsoil, qsol,
1932	cIM BAD $ paprs,pplay,radsol, fsnow,fqsurf,fevap,falbe,falblw,
1933	$ paprs,pplay, fsnow,fqsurf,fevap,falbe,falblw,
1934	$ fluxlat,
1935	e rain_fall, snow_fall,
1936	e fsolsw, fsollw, sollwdown, fder,
1937	e rlon, rlat, cuphy, cvphy, frugs,
1938	e debut, lafin, agesno,rugoro ,
1939	s d_t_vdf,d_q_vdf,d_u_vdf,d_v_vdf,d_ts,
1940	s fluxt,fluxq,fluxu,fluxv,cdragh,cdragm,
1941	s q2,
1942	s dsens, devap,
1943	s ycoefh,yu1,yv1, t2m, q2m, u10m, v10m,
1944	s pblh,capCL,oliqCL,cteiCL,pblT,
1945	s therm,trmb1,trmb2,trmb3,plcl,
1946	s fqcalving, ffonte, run_off_lic_0,
1947	cIM "slab" ocean
1948	s fluxo, fluxg, tslab, seaice)
1949	c
1950	CXXX PB
1951	CXXX Incrementation des flux
1952	CXXX
1953
1954	zxfluxt=0.
1955	zxfluxq=0.
1956	zxfluxu=0.
1957	zxfluxv=0.
1958	DO nsrf = 1, nbsrf
1959	DO k = 1, klev
1960	DO i = 1, klon
1961	zxfluxt(i,k) = zxfluxt(i,k) +
1962	$ fluxt(i,k,nsrf) * pctsrf( i, nsrf)
1963	zxfluxq(i,k) = zxfluxq(i,k) +
1964	$ fluxq(i,k,nsrf) * pctsrf( i, nsrf)
1965	zxfluxu(i,k) = zxfluxu(i,k) +
1966	$ fluxu(i,k,nsrf) * pctsrf( i, nsrf)
1967	zxfluxv(i,k) = zxfluxv(i,k) +
1968	$ fluxv(i,k,nsrf) * pctsrf( i, nsrf)
1969	END DO
1970	END DO
1971	END DO
1972	DO i = 1, klon
1973	sens(i) = - zxfluxt(i,1) ! flux de chaleur sensible au sol
1974	c evap(i) = - fluxq(i,1) ! flux d'evaporation au sol
1975	evap(i) = - zxfluxq(i,1) ! flux d'evaporation au sol
1976	fder(i) = dlw(i) + dsens(i) + devap(i)
1977	ENDDO
1978
1979
1980	DO k = 1, klev
1981	DO i = 1, klon
1982	t_seri(i,k) = t_seri(i,k) + d_t_vdf(i,k)
1983	q_seri(i,k) = q_seri(i,k) + d_q_vdf(i,k)
1984	u_seri(i,k) = u_seri(i,k) + d_u_vdf(i,k)
1985	v_seri(i,k) = v_seri(i,k) + d_v_vdf(i,k)
1986	ENDDO
1987	ENDDO
1988	cIM
1989	IF (ip_ebil_phy.ge.2) THEN
1990	ztit='after clmain'
1991	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
1992	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1993	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1994	call diagphy(airephy,ztit,ip_ebil_phy
1995	e , zero_v, zero_v, zero_v, zero_v, sens
1996	e , evap , zero_v, zero_v, ztsol
1997	e , d_h_vcol, d_qt, d_ec
1998	s , fs_bound, fq_bound )
1999	END IF
2000	C
2001	c
2002	c Incrementer la temperature du sol
2003	c
2004	DO i = 1, klon
2005	zxtsol(i) = 0.0
2006	zxfluxlat(i) = 0.0
2007	c
2008	zt2m(i) = 0.0
2009	zq2m(i) = 0.0
2010	zu10m(i) = 0.0
2011	zv10m(i) = 0.0
2012	cIM cf JLD ??
2013	zxffonte(i) = 0.0
2014	zxfqcalving(i) = 0.0
2015	cIM cf. AM 081204 BEG
2016	c
2017	s_pblh(i) = 0.0
2018	s_lcl(i) = 0.0
2019	s_capCL(i) = 0.0
2020	s_oliqCL(i) = 0.0
2021	s_cteiCL(i) = 0.0
2022	s_pblT(i) = 0.0
2023	s_therm(i) = 0.0
2024	s_trmb1(i) = 0.0
2025	s_trmb2(i) = 0.0
2026	s_trmb3(i) = 0.0
2027	c
2028	IF ( abs( pctsrf(i, is_ter) + pctsrf(i, is_lic) +
2029	$ pctsrf(i, is_oce) + pctsrf(i, is_sic) - 1.) .GT. EPSFRA)
2030	$ THEN
2031	WRITE(,) 'physiq : pb sous surface au point ', i,
2032	$ pctsrf(i, 1 : nbsrf)
2033	ENDIF
2034	ENDDO
2035	DO nsrf = 1, nbsrf
2036	DO i = 1, klon
2037	c IF (pctsrf(i,nsrf) .GE. EPSFRA) THEN
2038	ftsol(i,nsrf) = ftsol(i,nsrf) + d_ts(i,nsrf)
2039	cIM cf. JLD
2040	wfbils(i,nsrf) = ( fsolsw(i,nsrf) + fsollw(i,nsrf)
2041	$ + fluxt(i,1,nsrf) + fluxlat(i,nsrf) ) * pctsrf(i,nsrf)
2042	cIM
2043	wfbilo(i,nsrf) = ( fevap(i,nsrf) -
2044	$ (rain_fall(i) + snow_fall(i)) ) * pctsrf(i,nsrf)
2045	zxtsol(i) = zxtsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
2046	zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf)*pctsrf(i,nsrf)
2047	cccIM
2048	zt2m(i) = zt2m(i) + t2m(i,nsrf)*pctsrf(i,nsrf)
2049	zq2m(i) = zq2m(i) + q2m(i,nsrf)*pctsrf(i,nsrf)
2050	zu10m(i) = zu10m(i) + u10m(i,nsrf)*pctsrf(i,nsrf)
2051	zv10m(i) = zv10m(i) + v10m(i,nsrf)*pctsrf(i,nsrf)
2052	cIM cf JLD ??
2053	zxffonte(i) = zxffonte(i) + ffonte(i,nsrf)*pctsrf(i,nsrf)
2054	zxfqcalving(i) = zxfqcalving(i) +
2055	. fqcalving(i,nsrf)*pctsrf(i,nsrf)
2056	cIM cf. AM 081204 BEG
2057	s_pblh(i) = s_pblh(i) + pblh(i,nsrf)*pctsrf(i,nsrf)
2058	s_lcl(i) = s_lcl(i) + plcl(i,nsrf)*pctsrf(i,nsrf)
2059	s_capCL(i) = s_capCL(i) + capCL(i,nsrf) *pctsrf(i,nsrf)
2060	s_oliqCL(i) = s_oliqCL(i) + oliqCL(i,nsrf) *pctsrf(i,nsrf)
2061	s_cteiCL(i) = s_cteiCL(i) + cteiCL(i,nsrf) *pctsrf(i,nsrf)
2062	s_pblT(i) = s_pblT(i) + pblT(i,nsrf) *pctsrf(i,nsrf)
2063	s_therm(i) = s_therm(i) + therm(i,nsrf) *pctsrf(i,nsrf)
2064	s_trmb1(i) = s_trmb1(i) + trmb1(i,nsrf) *pctsrf(i,nsrf)
2065	s_trmb2(i) = s_trmb2(i) + trmb2(i,nsrf) *pctsrf(i,nsrf)
2066	s_trmb3(i) = s_trmb3(i) + trmb3(i,nsrf) *pctsrf(i,nsrf)
2067	c ENDIF
2068	ENDDO
2069	ENDDO
2070
2071	IF (check) THEN
2072	amn=MIN(ftsol(1,is_ter),1000.)
2073	amx=MAX(ftsol(1,is_ter),-1000.)
2074	DO i=2, klon
2075	amn=MIN(ftsol(i,is_ter),amn)
2076	amx=MAX(ftsol(i,is_ter),amx)
2077	ENDDO
2078	c
2079	PRINT*,' debut apres d_ts min max ftsol',itap,amn,amx
2080	ENDIF !(check) THEN
2081	c
2082	c Si une sous-fraction n'existe pas, elle prend la temp. moyenne
2083	c
2084	DO nsrf = 1, nbsrf
2085	DO i = 1, klon
2086	IF (pctsrf(i,nsrf) .LT. epsfra.OR.t2m(i,nsrf).EQ.0.) THEN
2087	ftsol(i,nsrf) = zxtsol(i)
2088	t2m(i,nsrf) = zt2m(i)
2089	q2m(i,nsrf) = zq2m(i)
2090	u10m(i,nsrf) = zu10m(i)
2091	v10m(i,nsrf) = zv10m(i)
2092	ffonte(i,nsrf) = zxffonte(i)
2093	fqcalving(i,nsrf) = zxfqcalving(i)
2094	pblh(i,nsrf)=s_pblh(i)
2095	plcl(i,nsrf)=s_lcl(i)
2096	capCL(i,nsrf)=s_capCL(i)
2097	oliqCL(i,nsrf)=s_oliqCL(i)
2098	cteiCL(i,nsrf)=s_cteiCL(i)
2099	pblT(i,nsrf)=s_pblT(i)
2100	therm(i,nsrf)=s_therm(i)
2101	trmb1(i,nsrf)=s_trmb1(i)
2102	trmb2(i,nsrf)=s_trmb2(i)
2103	trmb3(i,nsrf)=s_trmb3(i)
2104	ENDIF
2105	ENDDO
2106	ENDDO
2107	c
2108	c Calculer la derive du flux infrarouge
2109	c
2110	cXXX DO nsrf = 1, nbsrf
2111	DO i = 1, klon
2112	cXXX IF (pctsrf(i,nsrf) .GE. EPSFRA) THEN
2113	dlw(i) = - 4.0RSIGMAzxtsol(i)**3
2114	cXXX . *(ftsol(i,nsrf)-zxtsol(i))
2115	cXXX . *pctsrf(i,nsrf)
2116	cXXX ENDIF
2117	cXXX ENDDO
2118	ENDDO
2119	c
2120	c Appeler la convection (au choix)
2121	c
2122	DO k = 1, klev
2123	DO i = 1, klon
2124	conv_q(i,k) = d_q_dyn(i,k)
2125	. + d_q_vdf(i,k)/dtime
2126	conv_t(i,k) = d_t_dyn(i,k)
2127	. + d_t_vdf(i,k)/dtime
2128	ENDDO
2129	ENDDO
2130	IF (check) THEN
2131	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2132	WRITE(lunout,*) "avantcon=", za
2133	ENDIF
2134	zx_ajustq = .FALSE.
2135	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
2136	IF (zx_ajustq) THEN
2137	DO i = 1, klon
2138	z_avant(i) = 0.0
2139	ENDDO
2140	DO k = 1, klev
2141	DO i = 1, klon
2142	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k))
2143	. *(paprs(i,k)-paprs(i,k+1))/RG
2144	ENDDO
2145	ENDDO
2146	ENDIF
2147	IF (iflag_con.EQ.1) THEN
2148	stop'reactiver le call conlmd dans physiq.F'
2149	c CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
2150	c . d_t_con, d_q_con,
2151	c . rain_con, snow_con, ibas_con, itop_con)
2152	ELSE IF (iflag_con.EQ.2) THEN
2153	CALL conflx(dtime, paprs, pplay, t_seri, q_seri,
2154	e conv_t, conv_q, zxfluxq(1,1), omega,
2155	s d_t_con, d_q_con, rain_con, snow_con,
2156	s pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
2157	s kcbot, kctop, kdtop, pmflxr, pmflxs)
2158	WHERE (rain_con < 0.) rain_con = 0.
2159	WHERE (snow_con < 0.) snow_con = 0.
2160	DO i = 1, klon
2161	ibas_con(i) = klev+1 - kcbot(i)
2162	itop_con(i) = klev+1 - kctop(i)
2163	ENDDO
2164	ELSE IF (iflag_con.GE.3) THEN
2165	c nb of tracers for the KE convection:
2166	c MAF la partie traceurs est faite dans phytrac
2167	c on met ntra=1 pour limiter les appels mais on peut
2168	c supprimer les calculs / ftra.
2169	ntra = 1
2170	c sb, oct02:
2171	c Schema de convection modularise et vectorise:
2172	c (driver commun aux versions 3 et 4)
2173	c
2174	IF (ok_cvl) THEN ! new driver for convectL
2175
2176	CALL concvl (iflag_con,
2177	. dtime,paprs,pplay,t_seri,q_seri,
2178	. u_seri,v_seri,tr_seri,ntra,
2179	. ema_work1,ema_work2,
2180	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2181	. rain_con, snow_con, ibas_con, itop_con,
2182	. upwd,dnwd,dnwd0,
2183	. Ma,cape,tvp,iflagctrl,
2184	. pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd,
2185	. pmflxr,pmflxs,
2186	. da,phi,mp)
2187
2188	cIM cf. FH
2189	clwcon0=qcondc
2190	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
2191
2192	ELSE ! ok_cvl
2193	c MAF conema3 ne contient pas les traceurs
2194	CALL conema3 (dtime,
2195	. paprs,pplay,t_seri,q_seri,
2196	. u_seri,v_seri,tr_seri,ntra,
2197	. ema_work1,ema_work2,
2198	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2199	. rain_con, snow_con, ibas_con, itop_con,
2200	. upwd,dnwd,dnwd0,bas,top,
2201	. Ma,cape,tvp,rflag,
2202	. pbase
2203	. ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr
2204	. ,clwcon0)
2205
2206	ENDIF ! ok_cvl
2207
2208	IF (.NOT. ok_gust) THEN
2209	do i = 1, klon
2210	wd(i)=0.0
2211	enddo
2212	ENDIF
2213
2214	c =================================================================== c
2215	c Calcul des proprietes des nuages convectifs
2216	c
2217	DO k = 1, klev
2218	DO i = 1, klon
2219	zx_t = t_seri(i,k)
2220	IF (thermcep) THEN
2221	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2222	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2223	zx_qs = MIN(0.5,zx_qs)
2224	zcor = 1./(1.-retv*zx_qs)
2225	zx_qs = zx_qs*zcor
2226	ELSE
2227	IF (zx_t.LT.t_coup) THEN
2228	zx_qs = qsats(zx_t)/pplay(i,k)
2229	ELSE
2230	zx_qs = qsatl(zx_t)/pplay(i,k)
2231	ENDIF
2232	ENDIF
2233	zqsat(i,k)=zx_qs
2234	ENDDO
2235	ENDDO
2236
2237	c calcul des proprietes des nuages convectifs
2238	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
2239	call clouds_gno
2240	s (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
2241
2242	c =================================================================== c
2243
2244	DO i = 1, klon
2245	ema_pcb(i) = pbase(i)
2246	ENDDO
2247	DO i = 1, klon
2248	ema_pct(i) = paprs(i,itop_con(i))
2249	ENDDO
2250	DO i = 1, klon
2251	ema_cbmf(i) = ema_workcbmf(i)
2252	ENDDO
2253	ELSE
2254	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
2255	CALL abort
2256	ENDIF
2257
2258	c CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
2259	c . d_u_con, d_v_con)
2260
2261	DO k = 1, klev
2262	DO i = 1, klon
2263	t_seri(i,k) = t_seri(i,k) + d_t_con(i,k)
2264	q_seri(i,k) = q_seri(i,k) + d_q_con(i,k)
2265	u_seri(i,k) = u_seri(i,k) + d_u_con(i,k)
2266	v_seri(i,k) = v_seri(i,k) + d_v_con(i,k)
2267	ENDDO
2268	ENDDO
2269	cIM
2270	IF (ip_ebil_phy.ge.2) THEN
2271	ztit='after convect'
2272	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2273	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2274	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2275	call diagphy(airephy,ztit,ip_ebil_phy
2276	e , zero_v, zero_v, zero_v, zero_v, zero_v
2277	e , zero_v, rain_con, snow_con, ztsol
2278	e , d_h_vcol, d_qt, d_ec
2279	s , fs_bound, fq_bound )
2280	END IF
2281	C
2282	IF (check) THEN
2283	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2284	WRITE(lunout,*)"aprescon=", za
2285	zx_t = 0.0
2286	za = 0.0
2287	DO i = 1, klon
2288	za = za + airephy(i)/FLOAT(klon)
2289	zx_t = zx_t + (rain_con(i)+
2290	. snow_con(i))*airephy(i)/FLOAT(klon)
2291	ENDDO
2292	zx_t = zx_t/za*dtime
2293	WRITE(lunout,*)"Precip=", zx_t
2294	ENDIF
2295	IF (zx_ajustq) THEN
2296	DO i = 1, klon
2297	z_apres(i) = 0.0
2298	ENDDO
2299	DO k = 1, klev
2300	DO i = 1, klon
2301	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k))
2302	. *(paprs(i,k)-paprs(i,k+1))/RG
2303	ENDDO
2304	ENDDO
2305	DO i = 1, klon
2306	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime)
2307	. /z_apres(i)
2308	ENDDO
2309	DO k = 1, klev
2310	DO i = 1, klon
2311	IF (z_factor(i).GT.(1.0+1.0E-08) .OR.
2312	. z_factor(i).LT.(1.0-1.0E-08)) THEN
2313	q_seri(i,k) = q_seri(i,k) * z_factor(i)
2314	ENDIF
2315	ENDDO
2316	ENDDO
2317	ENDIF
2318	zx_ajustq=.FALSE.
2319	c
2320	c===================================================================
2321	c Convection seche (thermiques ou ajustement)
2322	c===================================================================
2323	c
2324	d_t_ajs(:,:)=0.
2325	d_u_ajs(:,:)=0.
2326	d_v_ajs(:,:)=0.
2327	d_q_ajs(:,:)=0.
2328	fm_therm(:,:)=0.
2329	entr_therm(:,:)=0.
2330	c
2331	IF(prt_level>9)WRITE(lunout,*)
2332	. 'AVANT LA CONVECTION SECHE , iflag_thermals='
2333	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2334	if(iflag_thermals.lt.0) then
2335	c Rien
2336	c ====
2337	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
2338	else if(iflag_thermals.eq.0) then
2339
2340	c Ajustement sec
2341	c ==============
2342	IF(prt_level>9)WRITE(lunout,*)'ajsec'
2343	CALL ajsec(paprs, pplay, t_seri,q_seri, d_t_ajs, d_q_ajs)
2344	t_seri(:,:) = t_seri(:,:) + d_t_ajs(:,:)
2345	q_seri(:,:) = q_seri(:,:) + d_q_ajs(:,:)
2346	else
2347	c Thermiques
2348	c ==========
2349	IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals='
2350	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2351	call calltherm(pdtphys
2352	s ,pplay,paprs,pphi
2353	s ,u_seri,v_seri,t_seri,q_seri
2354	s ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
2355	s ,fm_therm,entr_therm)
2356	endif
2357	c
2358	c===================================================================
2359	cIM
2360	IF (ip_ebil_phy.ge.2) THEN
2361	ztit='after dry_adjust'
2362	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2363	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2364	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2365	END IF
2366
2367
2368	c-------------------------------------------------------------------------
2369	c Caclul des ratqs
2370	c-------------------------------------------------------------------------
2371
2372	c print*,'calcul des ratqs'
2373	c ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
2374	c ----------------
2375	c on ecrase le tableau ratqsc calcule par clouds_gno
2376	if (iflag_cldcon.eq.1) then
2377	do k=1,klev
2378	do i=1,klon
2379	if(ptconv(i,k)) then
2380	ratqsc(i,k)=ratqsbas
2381	s +fact_cldcon*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
2382	else
2383	ratqsc(i,k)=0.
2384	endif
2385	enddo
2386	enddo
2387	endif
2388
2389	c ratqs stables
2390	c -------------
2391	do k=1,klev
2392	cIM RAJOUT boucle do=i
2393	do i=1, klon
2394	cIM ratqss(:,k)=ratqsbas+(ratqshaut-ratqsbas)*
2395	cIM s min((paprs(:,1)-pplay(:,k))/(paprs(:,1)-30000.),1.)
2396	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2397	s min((paprs(i,1)-pplay(i,k))/(paprs(i,1)-30000.),1.)
2398	cIM print*,' IMratqs STABLE i, k',i,k,ratqss(i,k)
2399	enddo
2400	enddo
2401
2402
2403	c ratqs final
2404	c -----------
2405	if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2) then
2406	c les ratqs sont une conbinaison de ratqss et ratqsc
2407	c ratqs final
2408	c 1e4 (en gros 3 heures), en dur pour le moment, est le temps de
2409	c relaxation des ratqs
2410	c facttemps=exp(-pdtphys/1.e4)
2411	facteur=exp(-pdtphys*facttemps)
2412	ratqs(:,:)=max(ratqs(:,:)*facteur,ratqss(:,:))
2413	ratqs(:,:)=max(ratqs(:,:),ratqsc(:,:))
2414	c print*,'calcul des ratqs fini'
2415	else
2416	c on ne prend que le ratqs stable pour fisrtilp
2417	ratqs(:,:)=ratqss(:,:)
2418	endif
2419
2420
2421	c
2422	c Appeler le processus de condensation a grande echelle
2423	c et le processus de precipitation
2424	c-------------------------------------------------------------------------
2425	CALL fisrtilp(dtime,paprs,pplay,
2426	. t_seri, q_seri,ptconv,ratqs,
2427	. d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq,
2428	. rain_lsc, snow_lsc,
2429	. pfrac_impa, pfrac_nucl, pfrac_1nucl,
2430	. frac_impa, frac_nucl,
2431	. prfl, psfl, rhcl)
2432
2433	WHERE (rain_lsc < 0) rain_lsc = 0.
2434	WHERE (snow_lsc < 0) snow_lsc = 0.
2435	DO k = 1, klev
2436	DO i = 1, klon
2437	t_seri(i,k) = t_seri(i,k) + d_t_lsc(i,k)
2438	q_seri(i,k) = q_seri(i,k) + d_q_lsc(i,k)
2439	ql_seri(i,k) = ql_seri(i,k) + d_ql_lsc(i,k)
2440	cldfra(i,k) = rneb(i,k)
2441	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
2442	ENDDO
2443	ENDDO
2444	IF (check) THEN
2445	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2446	WRITE(lunout,*)"apresilp=", za
2447	zx_t = 0.0
2448	za = 0.0
2449	DO i = 1, klon
2450	za = za + airephy(i)/FLOAT(klon)
2451	zx_t = zx_t + (rain_lsc(i)
2452	. + snow_lsc(i))*airephy(i)/FLOAT(klon)
2453	ENDDO
2454	zx_t = zx_t/za*dtime
2455	WRITE(lunout,*)"Precip=", zx_t
2456	ENDIF
2457	cIM
2458	IF (ip_ebil_phy.ge.2) THEN
2459	ztit='after fisrt'
2460	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2461	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2462	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2463	call diagphy(airephy,ztit,ip_ebil_phy
2464	e , zero_v, zero_v, zero_v, zero_v, zero_v
2465	e , zero_v, rain_lsc, snow_lsc, ztsol
2466	e , d_h_vcol, d_qt, d_ec
2467	s , fs_bound, fq_bound )
2468	END IF
2469	c
2470	c-------------------------------------------------------------------
2471	c PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
2472	c-------------------------------------------------------------------
2473
2474	c 1. NUAGES CONVECTIFS
2475	c
2476	cIM cf FH
2477	c IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
2478	IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
2479	snow_tiedtke=0.
2480	c print*,'avant calcul de la pseudo precip '
2481	c print*,'iflag_cldcon',iflag_cldcon
2482	if (iflag_cldcon.eq.-1) then
2483	rain_tiedtke=rain_con
2484	else
2485	c print*,'calcul de la pseudo precip '
2486	rain_tiedtke=0.
2487	c print*,'calcul de la pseudo precip 0'
2488	do k=1,klev
2489	do i=1,klon
2490	if (d_q_con(i,k).lt.0.) then
2491	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys
2492	s *(paprs(i,k)-paprs(i,k+1))/rg
2493	endif
2494	enddo
2495	enddo
2496	endif
2497	c
2498	c call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
2499	c
2500
2501	c Nuages diagnostiques pour Tiedtke
2502	CALL diagcld1(paprs,pplay,
2503	cIM cf FH . rain_con,snow_con,ibas_con,itop_con,
2504	. rain_tiedtke,snow_tiedtke,ibas_con,itop_con,
2505	. diafra,dialiq)
2506	DO k = 1, klev
2507	DO i = 1, klon
2508	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2509	cldliq(i,k) = dialiq(i,k)
2510	cldfra(i,k) = diafra(i,k)
2511	ENDIF
2512	ENDDO
2513	ENDDO
2514
2515	ELSE IF (iflag_cldcon.eq.3) THEN
2516	c On prend pour les nuages convectifs le max du calcul de la
2517	c convection et du calcul du pas de temps précédent diminué d'un facteur
2518	c facttemps
2519	c facttemps=pdtphys/1.e4
2520	facteur = pdtphys *facttemps
2521	do k=1,klev
2522	do i=1,klon
2523	rnebcon(i,k)=rnebcon(i,k)*facteur
2524	if (rnebcon0(i,k)clwcon0(i,k).gt.rnebcon(i,k)clwcon(i,k))
2525	s then
2526	rnebcon(i,k)=rnebcon0(i,k)
2527	clwcon(i,k)=clwcon0(i,k)
2528	endif
2529	enddo
2530	enddo
2531
2532	c
2533	cIM calcul nuages par le simulateur ISCCP
2534	c
2535	IF (ok_isccp) THEN
2536	#include "calcul_simulISCCP.h"
2537	ENDIF !ok_isccp
2538
2539	c On prend la somme des fractions nuageuses et des contenus en eau
2540	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
2541	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
2542
2543	ENDIF
2544
2545	c
2546	c 2. NUAGES STARTIFORMES
2547	c
2548	IF (ok_stratus) THEN
2549	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
2550	DO k = 1, klev
2551	DO i = 1, klon
2552	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2553	cldliq(i,k) = dialiq(i,k)
2554	cldfra(i,k) = diafra(i,k)
2555	ENDIF
2556	ENDDO
2557	ENDDO
2558	ENDIF
2559	c
2560	c Precipitation totale
2561	c
2562	DO i = 1, klon
2563	rain_fall(i) = rain_con(i) + rain_lsc(i)
2564	snow_fall(i) = snow_con(i) + snow_lsc(i)
2565	ENDDO
2566	cIM
2567	IF (ip_ebil_phy.ge.2) THEN
2568	ztit="after diagcld"
2569	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2570	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2571	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2572	END IF
2573	c
2574	c Calculer l'humidite relative pour diagnostique
2575	c
2576	DO k = 1, klev
2577	DO i = 1, klon
2578	zx_t = t_seri(i,k)
2579	IF (thermcep) THEN
2580	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2581	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2582	zx_qs = MIN(0.5,zx_qs)
2583	zcor = 1./(1.-retv*zx_qs)
2584	zx_qs = zx_qs*zcor
2585	ELSE
2586	IF (zx_t.LT.t_coup) THEN
2587	zx_qs = qsats(zx_t)/pplay(i,k)
2588	ELSE
2589	zx_qs = qsatl(zx_t)/pplay(i,k)
2590	ENDIF
2591	ENDIF
2592	zx_rh(i,k) = q_seri(i,k)/zx_qs
2593	zqsat(i,k)=zx_qs
2594	ENDDO
2595	ENDDO
2596	c
2597	cIM Calculer l'humidite relative a 2m (rh2m) pour diagnostique
2598	cIM ajout dependance type surface
2599	DO i = 1, klon
2600	rh2m(i)=0.
2601	qsat2m(i)=0.
2602	DO nsrf=1, nbsrf
2603	zx_t1(i,nsrf) = t2m(i,nsrf)
2604	IF (thermcep) THEN
2605	zdelta1(i,nsrf) = MAX(0.,SIGN(1.,rtt-zx_t1(i,nsrf)))
2606	zx_qs1(i,nsrf) = r2es *
2607	$ FOEEW(zx_t1(i,nsrf),zdelta1(i,nsrf))/paprs(i,1)
2608	zx_qs1(i,nsrf) = MIN(0.5,zx_qs1(i,nsrf))
2609	zcor1(i,nsrf) = 1./(1.-retv*zx_qs1(i,nsrf))
2610	zx_qs1(i,nsrf) = zx_qs1(i,nsrf)*zcor1(i,nsrf)
2611	ELSE
2612	c
2613	IF (zx_t.LT.RTT) THEN
2614	zx_qs = qsats(zx_t)/paprs(i,1)
2615	ELSE
2616	zx_qs = qsatl(zx_t)/paprs(i,1)
2617	ENDIF
2618	ENDIF
2619	zx_rh2m(i,nsrf) = q2m(i,nsrf)/zx_qs1(i,nsrf)
2620	zx_qsat2m(i,nsrf)=zx_qs1(i,nsrf)
2621	rh2m(i) = rh2m(i)+zx_rh2m(i,nsrf)*pctsrf(i,nsrf)
2622	qsat2m(i)=qsat2m(i)+zx_qsat2m(i,nsrf)*pctsrf(i,nsrf)
2623	ENDDO !nsrf
2624	ENDDO
2625	c
2626	cIM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
2627	c equivalente a 2m (tpote) pour diagnostique
2628	c
2629	DO i = 1, klon
2630	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
2631	IF (thermcep) THEN
2632	IF(zt2m(i).LT.RTT) then
2633	Lheat=RLSTT
2634	ELSE
2635	Lheat=RLVTT
2636	ENDIF
2637	ELSE
2638	IF (zt2m(i).LT.RTT) THEN
2639	Lheat=RLSTT
2640	ELSE
2641	Lheat=RLVTT
2642	ENDIF
2643	ENDIF
2644	tpote(i) = tpot(i)*
2645	. EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
2646	ENDDO
2647	c
2648	cjq - introduce the aerosol direct and first indirect radiative forcings
2649	cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
2650	IF (ok_ade.OR.ok_aie) THEN
2651	! Get sulfate aerosol distribution
2652	CALL readsulfate(rjourvrai, debut, sulfate)
2653	CALL readsulfate_preind(rjourvrai, debut, sulfate_pi)
2654
2655	! Calculate aerosol optical properties (Olivier Boucher)
2656	CALL aeropt(pplay, paprs, t_seri, sulfate, rhcl,
2657	. tau_ae, piz_ae, cg_ae, aerindex)
2658	cym
2659	ELSE
2660	tau_ae(:,:,:)=0.0
2661	piz_ae(:,:,:)=0.0
2662	cg_ae(:,:,:)=0.0
2663	cym
2664	ENDIF
2665
2666	#ifdef INCA
2667	calday = FLOAT(julien) + gmtime
2668
2669	#ifdef INCA_AER
2670	call AEROSOL_METEO_CALC(calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs,
2671	& prfl,psfl,pctsrf(1,3),airephy,xjour,rlat,rlon)
2672	#endif
2673
2674	#ifdef INCAINFO
2675	WRITE(lunout,*)'Appel CHEMHOOK_BEGIN ...'
2676	#endif
2677
2678	CALL chemhook_begin (calday,
2679	#if defined(INCA) && !defined(INCA_CH4) && !defined(INCA_NMHC) && !defined(INCA_AER)
2680	$ julien,
2681	$ gmtime,
2682	#endif
2683	$ pctsrf(1,1),
2684	$ rlat,
2685	$ rlon,
2686	$ airephy,
2687	$ paprs,
2688	$ pplay,
2689	$ ycoefh,
2690	$ pphi,
2691	$ t_seri,
2692	$ u,
2693	$ v,
2694	$ wo,
2695	$ q_seri,
2696	$ zxtsol,
2697	$ zxsnow,
2698	$ solsw,
2699	$ albsol,
2700	$ rain_fall,
2701	$ snow_fall,
2702	$ itop_con,
2703	$ ibas_con,
2704	$ cldfra,
2705	$ iim,
2706	$ jjm,
2707	#ifdef INCA_AER
2708	$ tr_seri,
2709	$ ftsol,
2710	$ paprs,
2711	$ cdragh,
2712	$ cdragm,
2713	$ pctsrf,
2714	$ pdtphys,
2715	$ itap)
2716	#else
2717	$ tr_seri)
2718	#endif
2719
2720
2721	#ifdef INCAINFO
2722	WRITE(lunout,*)'OK.'
2723	#endif
2724	#endif
2725	c
2726	c Calculer les parametres optiques des nuages et quelques
2727	c parametres pour diagnostiques:
2728	c
2729	if (ok_newmicro) then
2730	CALL newmicro (paprs, pplay,ok_newmicro,
2731	. t_seri, cldliq, cldfra, cldtau, cldemi,
2732	. cldh, cldl, cldm, cldt, cldq,
2733	. flwp, fiwp, flwc, fiwc,
2734	e ok_aie,
2735	e sulfate, sulfate_pi,
2736	e bl95_b0, bl95_b1,
2737	s cldtaupi, re, fl)
2738	else
2739	CALL nuage (paprs, pplay,
2740	. t_seri, cldliq, cldfra, cldtau, cldemi,
2741	. cldh, cldl, cldm, cldt, cldq,
2742	e ok_aie,
2743	e sulfate, sulfate_pi,
2744	e bl95_b0, bl95_b1,
2745	s cldtaupi, re, fl)
2746
2747	endif
2748	c
2749	c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
2750	c
2751	IF (MOD(itaprad,radpas).EQ.0) THEN
2752	DO i = 1, klon
2753	albsol(i) = falbe(i,is_oce) * pctsrf(i,is_oce)
2754	. + falbe(i,is_lic) * pctsrf(i,is_lic)
2755	. + falbe(i,is_ter) * pctsrf(i,is_ter)
2756	. + falbe(i,is_sic) * pctsrf(i,is_sic)
2757	albsollw(i) = falblw(i,is_oce) * pctsrf(i,is_oce)
2758	. + falblw(i,is_lic) * pctsrf(i,is_lic)
2759	. + falblw(i,is_ter) * pctsrf(i,is_ter)
2760	. + falblw(i,is_sic) * pctsrf(i,is_sic)
2761	ENDDO
2762	CALL radlwsw ! nouveau rayonnement (compatible Arpege-IFS)
2763	e (dist, rmu0, fract,
2764	e paprs, pplay,zxtsol,albsol, albsollw, t_seri,q_seri,
2765	e wo,
2766	e cldfra, cldemi, cldtau,
2767	s heat,heat0,cool,cool0,radsol,albpla,
2768	s topsw,toplw,solsw,sollw,
2769	s sollwdown,
2770	s topsw0,toplw0,solsw0,sollw0,
2771	s lwdn0, lwdn, lwup0, lwup,
2772	s swdn0, swdn, swup0, swup,
2773	e ok_ade, ok_aie, ! new for aerosol radiative effects
2774	e tau_ae, piz_ae, cg_ae, ! ="=
2775	s topswad, solswad, ! ="=
2776	e cldtaupi, ! ="=
2777	s topswai, solswai) ! ="=
2778	itaprad = 0
2779	ENDIF
2780	itaprad = itaprad + 1
2781	c
2782	c Ajouter la tendance des rayonnements (tous les pas)
2783	c
2784	DO k = 1, klev
2785	DO i = 1, klon
2786	t_seri(i,k) = t_seri(i,k)
2787	. + (heat(i,k)-cool(i,k)) * dtime/86400.
2788	ENDDO
2789	ENDDO
2790	cIM
2791	IF (ip_ebil_phy.ge.2) THEN
2792	ztit='after rad'
2793	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2794	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2795	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2796	call diagphy(airephy,ztit,ip_ebil_phy
2797	e , topsw, toplw, solsw, sollw, zero_v
2798	e , zero_v, zero_v, zero_v, ztsol
2799	e , d_h_vcol, d_qt, d_ec
2800	s , fs_bound, fq_bound )
2801	END IF
2802	c
2803	c
2804	c Calculer l'hydrologie de la surface
2805	c
2806	c CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
2807	c . agesno, ftsol,fqsurf,fsnow, ruis)
2808	c
2809	DO i = 1, klon
2810	zxqsurf(i) = 0.0
2811	zxsnow(i) = 0.0
2812	ENDDO
2813	DO nsrf = 1, nbsrf
2814	DO i = 1, klon
2815	zxqsurf(i) = zxqsurf(i) + fqsurf(i,nsrf)*pctsrf(i,nsrf)
2816	zxsnow(i) = zxsnow(i) + fsnow(i,nsrf)*pctsrf(i,nsrf)
2817	ENDDO
2818	ENDDO
2819	c
2820	c Si une sous-fraction n'existe pas, elle prend la valeur moyenne
2821	c
2822	cXXX DO nsrf = 1, nbsrf
2823	cXXX DO i = 1, klon
2824	cXXX IF (pctsrf(i,nsrf).LT.epsfra) THEN
2825	cXXX fqsurf(i,nsrf) = zxqsurf(i)
2826	cXXX fsnow(i,nsrf) = zxsnow(i)
2827	cXXX ENDIF
2828	cXXX ENDDO
2829	cXXX ENDDO
2830	c
2831	c Calculer le bilan du sol et la derive de temperature (couplage)
2832	c
2833	DO i = 1, klon
2834	c bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
2835	c a la demande de JLD
2836	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
2837	ENDDO
2838	c
2839	cmoddeblott(jan95)
2840	c Appeler le programme de parametrisation de l'orographie
2841	c a l'echelle sous-maille:
2842	c
2843	IF (ok_orodr) THEN
2844	c
2845	c selection des points pour lesquels le shema est actif:
2846	igwd=0
2847	DO i=1,klon
2848	itest(i)=0
2849	c IF ((zstd(i).gt.10.0)) THEN
2850	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
2851	itest(i)=1
2852	igwd=igwd+1
2853	idx(igwd)=i
2854	ENDIF
2855	ENDDO
2856	c igwdim=MAX(1,igwd)
2857	c
2858	CALL drag_noro(klon,klev,dtime,paprs,pplay,
2859	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
2860	e igwd,idx,itest,
2861	e t_seri, u_seri, v_seri,
2862	cIM 141004 s zulow, zvlow, zustr, zvstr,
2863	s zulow, zvlow, zustrdr, zvstrdr,
2864	s d_t_oro, d_u_oro, d_v_oro)
2865	c
2866	c ajout des tendances
2867	DO k = 1, klev
2868	DO i = 1, klon
2869	t_seri(i,k) = t_seri(i,k) + d_t_oro(i,k)
2870	u_seri(i,k) = u_seri(i,k) + d_u_oro(i,k)
2871	v_seri(i,k) = v_seri(i,k) + d_v_oro(i,k)
2872	ENDDO
2873	ENDDO
2874	c
2875	ENDIF ! fin de test sur ok_orodr
2876	c
2877	IF (ok_orolf) THEN
2878	c
2879	c selection des points pour lesquels le shema est actif:
2880	igwd=0
2881	DO i=1,klon
2882	itest(i)=0
2883	IF ((zpic(i)-zmea(i)).GT.100.) THEN
2884	itest(i)=1
2885	igwd=igwd+1
2886	idx(igwd)=i
2887	ENDIF
2888	ENDDO
2889	c igwdim=MAX(1,igwd)
2890	c
2891	CALL lift_noro(klon,klev,dtime,paprs,pplay,
2892	e rlat,zmea,zstd,zpic,
2893	e itest,
2894	e t_seri, u_seri, v_seri,
2895	s zulow, zvlow, zustrli, zvstrli,
2896	s d_t_lif, d_u_lif, d_v_lif)
2897	c
2898	c ajout des tendances
2899	DO k = 1, klev
2900	DO i = 1, klon
2901	t_seri(i,k) = t_seri(i,k) + d_t_lif(i,k)
2902	u_seri(i,k) = u_seri(i,k) + d_u_lif(i,k)
2903	v_seri(i,k) = v_seri(i,k) + d_v_lif(i,k)
2904	ENDDO
2905	ENDDO
2906	c
2907	ENDIF ! fin de test sur ok_orolf
2908	c
2909	cIM cf. FLott BEG
2910	C STRESS NECESSAIRES: TOUTE LA PHYSIQUE
2911
2912	DO i = 1, klon
2913	zustrph(i)=0.
2914	zvstrph(i)=0.
2915	ENDDO
2916	DO k = 1, klev
2917	DO i = 1, klon
2918	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
2919	c (paprs(i,k)-paprs(i,k+1))/rg
2920	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
2921	c (paprs(i,k)-paprs(i,k+1))/rg
2922	ENDDO
2923	ENDDO
2924	c
2925	cIM calcul composantes axiales du moment angulaire et couple des montagnes
2926	c
2927	CALL aaam_bud (27,klon,klev,rjourvrai,gmtime,
2928	C ra,rg,romega,
2929	C rlat,rlon,pphis,
2930	C zustrdr,zustrli,zustrph,
2931	C zvstrdr,zvstrli,zvstrph,
2932	C paprs,u,v,
2933	C aam, torsfc)
2934	cIM cf. FLott END
2935	cIM
2936	IF (ip_ebil_phy.ge.2) THEN
2937	ztit='after orography'
2938	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
2939	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2940	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2941	END IF
2942	c
2943	c
2944	cAA
2945	cAA Installation de l'interface online-offline pour traceurs
2946	cAA
2947	c====================================================================
2948	c Calcul des tendances traceurs
2949	c====================================================================
2950	C
2951	call phytrac ( rnpb,
2952	I itap,
2953	I julien,
2954	I gmtime,
2955	I debut,
2956	I lafin,
2957	I nqmax-2,
2958	I nlon,
2959	I nlev,
2960	I dtime,
2961	I u,
2962	I v,
2963	I t,
2964	I paprs,
2965	I pplay,
2966	I pmfu,
2967	I pmfd,
2968	I pen_u,
2969	I pde_u,
2970	I pen_d,
2971	I pde_d,
2972	I ycoefh,
2973	I fm_therm,
2974	I entr_therm,
2975	I yu1,
2976	I yv1,
2977	I ftsol,
2978	I pctsrf,
2979	I rlat,
2980	I frac_impa,
2981	I frac_nucl,
2982	I rlon,
2983	I presnivs,
2984	I pphis,
2985	I pphi,
2986	I albsol,
2987	I qx(1,1,1),
2988	I rhcl,
2989	I cldfra,
2990	I rneb,
2991	I diafra,
2992	I cldliq,
2993	I itop_con,
2994	I ibas_con,
2995	I pmflxr,
2996	I pmflxs,
2997	I prfl,
2998	I psfl,
2999	I da,
3000	I phi,
3001	I mp,
3002	I upwd,
3003	I dnwd,
3004	#ifdef INCA
3005	I flxmass_w,
3006	#endif
3007	O tr_seri)
3008
3009	IF (offline) THEN
3010
3011	print*,'Attention on met a 0 les thermiques pour phystoke'
3012	call phystokenc (
3013	I nlon,nlev,pdtphys,rlon,rlat,
3014	I t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
3015	I fm_therm,entr_therm,
3016	I ycoefh,yu1,yv1,ftsol,pctsrf,
3017	I frac_impa, frac_nucl,
3018	I pphis,airephy,dtime,itap)
3019
3020
3021	ENDIF
3022
3023	c
3024	c Calculer le transport de l'eau et de l'energie (diagnostique)
3025	c
3026	CALL transp (paprs,zxtsol,
3027	e t_seri, q_seri, u_seri, v_seri, zphi,
3028	s ve, vq, ue, uq)
3029	c
3030	cIM global posePB BEG
3031	IF(1.EQ.0) THEN
3032	c
3033	CALL transp_lay (paprs,zxtsol,
3034	e t_seri, q_seri, u_seri, v_seri, zphi,
3035	s ve_lay, vq_lay, ue_lay, uq_lay)
3036	c
3037	ENDIF !(1.EQ.0) THEN
3038	cIM global posePB END
3039	c Accumuler les variables a stocker dans les fichiers histoire:
3040	c
3041	c+jld ec_conser
3042	DO k = 1, klev
3043	DO i = 1, klon
3044	ZRCPD = RCPD(1.0+RVTMP2q_seri(i,k))
3045	d_t_ec(i,k)=0.5/ZRCPD
3046	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
3047	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
3048	d_t_ec(i,k) = d_t_ec(i,k)/dtime
3049	END DO
3050	END DO
3051	c-jld ec_conser
3052	cIM
3053	IF (ip_ebil_phy.ge.1) THEN
3054	ztit='after physic'
3055	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
3056	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
3057	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3058	C Comme les tendances de la physique sont ajoute dans la dynamique,
3059	C on devrait avoir que la variation d'entalpie par la dynamique
3060	C est egale a la variation de la physique au pas de temps precedent.
3061	C Donc la somme de ces 2 variations devrait etre nulle.
3062	call diagphy(airephy,ztit,ip_ebil_phy
3063	e , topsw, toplw, solsw, sollw, sens
3064	e , evap, rain_fall, snow_fall, ztsol
3065	e , d_h_vcol, d_qt, d_ec
3066	s , fs_bound, fq_bound )
3067	C
3068	d_h_vcol_phy=d_h_vcol
3069	C
3070	END IF
3071	C
3072	c=======================================================================
3073	c SORTIES
3074	c=======================================================================
3075
3076	cIM Interpolation sur les niveaux de pression du NMC
3077	c -------------------------------------------------
3078	c
3079	#include "calcul_STDlev.h"
3080	c
3081	c slp sea level pressure
3082	slp(:) = paprs(:,1)exp(pphis(:)/(RDt_seri(:,1)))
3083	c
3084	ccc prw = eau precipitable
3085	DO i = 1, klon
3086	prw(i) = 0.
3087	DO k = 1, klev
3088	prw(i) = prw(i) +
3089	. q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
3090	ENDDO
3091	ENDDO
3092	c
3093	cIM initialisation + calculs divers diag AMIP2
3094	c
3095	#include "calcul_divers.h"
3096	c
3097	#ifdef INCA
3098	#ifdef INCAINFO
3099	WRITE(lunout,*)'Appel CHEMHOOK_END ...'
3100	#endif
3101	CALL chemhook_end (calday,
3102	$ dtime,
3103	$ pplay,
3104	$ t_seri,
3105	$ tr_seri,
3106	$ nbtr,
3107	$ paprs,
3108	#ifdef INCA_CH4
3109	$ q_seri,
3110	#endif
3111	$ annee_ref,
3112	$ day_ini,
3113	#ifdef INCA_AER
3114	$ xjour,
3115	$ pphi,
3116	$ pphis,
3117	$ zx_rh,
3118	$ qx(1,1,1))
3119	#else
3120	$ xjour)
3121	#endif
3122	#ifdef INCAINFO
3123	WRITE(lunout,*)'OK.'
3124	#endif
3125	#endif
3126
3127	c=============================================================
3128	c
3129	c Convertir les incrementations en tendances
3130	c
3131	DO k = 1, klev
3132	DO i = 1, klon
3133	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
3134	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
3135	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
3136	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
3137	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
3138	ENDDO
3139	ENDDO
3140	c
3141	IF (nqmax.GE.3) THEN
3142	DO iq = 3, nqmax
3143	DO k = 1, klev
3144	DO i = 1, klon
3145	d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
3146	ENDDO
3147	ENDDO
3148	ENDDO
3149	ENDIF
3150	c
3151	cIM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
3152	cIM global posePB#include "write_bilKP_ins.h"
3153	cIM global posePB#include "write_bilKP_ave.h"
3154	c
3155	c Sauvegarder les valeurs de t et q a la fin de la physique:
3156	c
3157	DO k = 1, klev
3158	DO i = 1, klon
3159	t_ancien(i,k) = t_seri(i,k)
3160	q_ancien(i,k) = q_seri(i,k)
3161	ENDDO
3162	ENDDO
3163	c
3164	c 22.03.04 BEG
3165	c=============================================================
3166	c Ecriture des sorties
3167	c=============================================================
3168	#ifdef CPP_IOIPSL
3169
3170	#ifdef histhf
3171	#include "write_histhf.h"
3172	#endif
3173
3174	#ifdef histday
3175	#include "write_histday.h"
3176	#endif
3177
3178	#ifdef histmth
3179	#include "write_histmth.h"
3180	#endif
3181
3182	#ifdef histins
3183	#include "write_histins.h"
3184	#endif
3185
3186	#ifdef histISCCP
3187	#include "write_histISCCP.h"
3188	#endif
3189
3190	#ifdef histmthNMC
3191	#include "write_histmthNMC.h"
3192	#endif
3193
3194	#include "write_histday_seri.h"
3195
3196	#include "write_paramLMDZ_phy.h"
3197
3198	#endif
3199
3200	c 22.03.04 END
3201	c
3202	c====================================================================
3203	c Si c'est la fin, il faut conserver l'etat de redemarrage
3204	c====================================================================
3205	c
3206	IF (lafin) THEN
3207	itau_phy = itau_phy + itap
3208	ccc IF (ok_oasis) CALL quitcpl
3209	CALL phyredem ("restartphy.nc",dtime,radpas,
3210	. rlat, rlon, pctsrf, ftsol, ftsoil,
3211	cIM "slab" ocean
3212	. tslab, seaice,
3213	. fqsurf, qsol,
3214	. fsnow, falbe,falblw, fevap, rain_fall, snow_fall,
3215	cIM . solsw, sollwdown,dlw,
3216	. solsw, sollw,dlw,
3217	. radsol,frugs,agesno,
3218	. zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,
3219	. t_ancien, q_ancien, rnebcon, ratqs, clwcon,run_off_lic_0)
3220	ENDIF
3221
3222
3223	RETURN
3224	END
3225	FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
3226	IMPLICIT none
3227	c
3228	c Calculer et imprimer l'eau totale. A utiliser pour verifier
3229	c la conservation de l'eau
3230	c
3231	#include "YOMCST.h"
3232	INTEGER klon,klev
3233	REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
3234	REAL aire(klon)
3235	REAL qtotal, zx, qcheck
3236	INTEGER i, k
3237	c
3238	zx = 0.0
3239	DO i = 1, klon
3240	zx = zx + aire(i)
3241	ENDDO
3242	qtotal = 0.0
3243	DO k = 1, klev
3244	DO i = 1, klon
3245	qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i)
3246	. *(paprs(i,k)-paprs(i,k+1))/RG
3247	ENDDO
3248	ENDDO
3249	c
3250	qcheck = qtotal/zx
3251	c
3252	RETURN
3253	END
3254	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
3255	IMPLICIT none
3256	c
3257	c Tranformer une variable de la grille physique a
3258	c la grille d'ecriture
3259	c
3260	INTEGER nfield,nlon,iim,jjmp1, jjm
3261	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
3262	c
3263	INTEGER i, n, ig
3264	c
3265	jjm = jjmp1 - 1
3266	DO n = 1, nfield
3267	DO i=1,iim
3268	ecrit(i,n) = fi(1,n)
3269	ecrit(i+jjm*iim,n) = fi(nlon,n)
3270	ENDDO
3271	DO ig = 1, nlon - 2
3272	ecrit(iim+ig,n) = fi(1+ig,n)
3273	ENDDO
3274	ENDDO
3275	RETURN
3276	END

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

Download in other formats:

Original Format