Context Navigation

physiq.F @ 658

Last change on this file since 658 was 658, checked in by lmdzadmin, 19 years ago
frequence de sortie des traceurs (en jours) choisie dans physiq.def IM MAF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 96.1 KB

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

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

Download in other formats:

Original Format