Context Navigation

physiq.F @ 757

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

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

Download in other formats:

Original Format