Context Navigation

physiq.F @ 651

Last change on this file since 651 was 651, checked in by Laurent Fairhead, 19 years ago
On fait dépendre la lecture de tslab et seaice dans phyetat0 du flag ocean (on n'a pas besoin de les lire si on n'est pas en ocean slab) LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 96.2 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
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(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	#ifdef histISCCP
1523	#include "ini_histISCCP.h"
1524	#endif
1525
1526	c#undef histmthNMC
1527	c#define histmthNMC
1528	#ifdef histmthNMC
1529	#include "ini_histmthNMC.h"
1530	#endif
1531
1532	#endif
1533
1534	cXXXPB Positionner date0 pour initialisation de ORCHIDEE
1535	date0 = zjulian
1536	C date0 = day_ini
1537	WRITE(,) 'physiq date0 : ',date0
1538	c
1539	c
1540	c
1541	c Prescrire l'ozone dans l'atmosphere
1542	c
1543	c
1544	cc DO i = 1, klon
1545	cc DO k = 1, klev
1546	cc CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1547	cc ENDDO
1548	cc ENDDO
1549	c
1550	#ifdef INCA
1551	iii = MOD(NINT(xjour),360)
1552	calday = FLOAT(iii) + gmtime
1553	WRITE(lunout,*) 'initial time ', xjour, calday
1554	#ifdef INCAINFO
1555	WRITE(lunout,*) 'Appel CHEMINI ...'
1556	#endif
1557	CALL chemini( rpi,
1558	$ rg,
1559	$ ra,
1560	$ airephy,
1561	$ rlat,
1562	$ rlon,
1563	$ presnivs,
1564	$ calday,
1565	$ klon,
1566	$ nqmax,
1567	$ pdtphys,
1568	$ annee_ref,
1569	$ day_ini)
1570	#ifdef INCAINFO
1571	WRITE(lunout,*) 'OK.'
1572	#endif
1573	#endif
1574	c
1575	ENDIF
1576	c
1577	c ************** Fin de IF ( debut ) *************
1578	c
1579	c
1580	c Mettre a zero des variables de sortie (pour securite)
1581	c
1582	DO i = 1, klon
1583	d_ps(i) = 0.0
1584	ENDDO
1585	DO k = 1, klev
1586	DO i = 1, klon
1587	d_t(i,k) = 0.0
1588	d_u(i,k) = 0.0
1589	d_v(i,k) = 0.0
1590	ENDDO
1591	ENDDO
1592	DO iq = 1, nqmax
1593	DO k = 1, klev
1594	DO i = 1, klon
1595	d_qx(i,k,iq) = 0.0
1596	ENDDO
1597	ENDDO
1598	ENDDO
1599	c
1600	c Ne pas affecter les valeurs entrees de u, v, h, et q
1601	c
1602	DO k = 1, klev
1603	DO i = 1, klon
1604	t_seri(i,k) = t(i,k)
1605	u_seri(i,k) = u(i,k)
1606	v_seri(i,k) = v(i,k)
1607	q_seri(i,k) = qx(i,k,ivap)
1608	ql_seri(i,k) = qx(i,k,iliq)
1609	qs_seri(i,k) = 0.
1610	ENDDO
1611	ENDDO
1612	IF (nqmax.GE.3) THEN
1613	DO iq = 3, nqmax
1614	DO k = 1, klev
1615	DO i = 1, klon
1616	tr_seri(i,k,iq-2) = qx(i,k,iq)
1617	ENDDO
1618	ENDDO
1619	ENDDO
1620	ELSE
1621	DO k = 1, klev
1622	DO i = 1, klon
1623	tr_seri(i,k,1) = 0.0
1624	ENDDO
1625	ENDDO
1626	ENDIF
1627	C
1628	DO i = 1, klon
1629	ztsol(i) = 0.
1630	ENDDO
1631	DO nsrf = 1, nbsrf
1632	DO i = 1, klon
1633	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1634	ENDDO
1635	ENDDO
1636	C
1637	IF (if_ebil.ge.1) THEN
1638	ztit='after dynamic'
1639	CALL diagetpq(airephy,ztit,ip_ebil,1,1,dtime
1640	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1641	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1642	C Comme les tendances de la physique sont ajoute dans la dynamique,
1643	C on devrait avoir que la variation d'entalpie par la dynamique
1644	C est egale a la variation de la physique au pas de temps precedent.
1645	C Donc la somme de ces 2 variations devrait etre nulle.
1646	call diagphy(airephy,ztit,ip_ebil
1647	e , zero_v, zero_v, zero_v, zero_v, zero_v
1648	e , zero_v, zero_v, zero_v, ztsol
1649	e , d_h_vcol+d_h_vcol_phy, d_qt, 0.
1650	s , fs_bound, fq_bound )
1651	END IF
1652
1653	c Diagnostiquer la tendance dynamique
1654	c
1655	IF (ancien_ok) THEN
1656	DO k = 1, klev
1657	DO i = 1, klon
1658	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
1659	d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
1660	ENDDO
1661	ENDDO
1662	ELSE
1663	DO k = 1, klev
1664	DO i = 1, klon
1665	d_t_dyn(i,k) = 0.0
1666	d_q_dyn(i,k) = 0.0
1667	ENDDO
1668	ENDDO
1669	ancien_ok = .TRUE.
1670	ENDIF
1671	c
1672	c Ajouter le geopotentiel du sol:
1673	c
1674	DO k = 1, klev
1675	DO i = 1, klon
1676	zphi(i,k) = pphi(i,k) + pphis(i)
1677	ENDDO
1678	ENDDO
1679	c
1680	c Verifier les temperatures
1681	c
1682	CALL hgardfou(t_seri,ftsol,'debutphy')
1683	c
1684	c Incrementer le compteur de la physique
1685	c
1686	itap = itap + 1
1687	julien = MOD(NINT(xjour),360)
1688	if (julien .eq. 0) julien = 360
1689	c
1690	c Mettre en action les conditions aux limites (albedo, sst, etc.).
1691	c Prescrire l'ozone et calculer l'albedo sur l'ocean.
1692	c
1693	IF (MOD(itap-1,lmt_pas) .EQ. 0) THEN
1694	WRITE(lunout,*)' PHYS cond julien ',julien
1695	CALL ozonecm( FLOAT(julien), rlat, paprs, wo)
1696	ENDIF
1697	c
1698	c Re-evaporer l'eau liquide nuageuse
1699	c
1700	DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse
1701	DO i = 1, klon
1702	zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1703	c zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1704	zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1705	zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
1706	zb = MAX(0.0,ql_seri(i,k))
1707	za = - MAX(0.0,ql_seri(i,k))
1708	. * (zlvdcp(1.-zdelta)+zlsdcpzdelta)
1709	t_seri(i,k) = t_seri(i,k) + za
1710	q_seri(i,k) = q_seri(i,k) + zb
1711	ql_seri(i,k) = 0.0
1712	d_t_eva(i,k) = za
1713	d_q_eva(i,k) = zb
1714	ENDDO
1715	ENDDO
1716	c
1717	IF (if_ebil.ge.2) THEN
1718	ztit='after reevap'
1719	CALL diagetpq(airephy,ztit,ip_ebil,2,1,dtime
1720	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1721	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1722	call diagphy(airephy,ztit,ip_ebil
1723	e , zero_v, zero_v, zero_v, zero_v, zero_v
1724	e , zero_v, zero_v, zero_v, ztsol
1725	e , d_h_vcol, d_qt, d_ec
1726	s , fs_bound, fq_bound )
1727	C
1728	END IF
1729	C
1730	c
1731	c Appeler la diffusion verticale (programme de couche limite)
1732	c
1733	DO i = 1, klon
1734	c if (.not. ok_veget) then
1735	c frugs(i,is_ter) = SQRT(frugs(i,is_ter)2+rugoro(i)2)
1736	c endif
1737	c frugs(i,is_lic) = rugoro(i)
1738	c frugs(i,is_oce) = rugmer(i)
1739	c frugs(i,is_sic) = 0.001
1740	zxrugs(i) = 0.0
1741	ENDDO
1742	DO nsrf = 1, nbsrf
1743	DO i = 1, klon
1744	c frugs(i,nsrf) = MAX(frugs(i,nsrf),0.001)
1745	frugs(i,nsrf) = MAX(frugs(i,nsrf),0.000015)
1746	ENDDO
1747	ENDDO
1748	DO nsrf = 1, nbsrf
1749	DO i = 1, klon
1750	zxrugs(i) = zxrugs(i) + frugs(i,nsrf)*pctsrf(i,nsrf)
1751	ENDDO
1752	ENDDO
1753	c
1754	C calculs necessaires au calcul de l'albedo dans l'interface
1755	c
1756	CALL orbite(FLOAT(julien),zlongi,dist)
1757	IF (cycle_diurne) THEN
1758	zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
1759	CALL zenang(zlongi,gmtime,zdtime,rlat,rlon,rmu0,fract)
1760	ELSE
1761	rmu0 = -999.999
1762	ENDIF
1763	c
1764	C Calcul de l'abedo moyen par maille
1765	albsol(:)=0.
1766	albsollw(:)=0.
1767	DO nsrf = 1, nbsrf
1768	DO i = 1, klon
1769	albsol(i) = albsol(i) + falbe(i,nsrf) * pctsrf(i,nsrf)
1770	albsollw(i) = albsollw(i) + falblw(i,nsrf) * pctsrf(i,nsrf)
1771	ENDDO
1772	ENDDO
1773	C
1774	C Repartition sous maille des flux LW et SW
1775	C Modif OM+PASB+JLD
1776	C Repartition du longwave par sous-surface linearisee
1777	Cn
1778
1779	DO nsrf = 1, nbsrf
1780	DO i = 1, klon
1781	c$$$ fsollw(i,nsrf) = sollwdown(i) - RSIGMAftsol(i,nsrf)*4
1782	c$$$ fsollw(i,nsrf) = sollw(i)
1783	fsollw(i,nsrf) = sollw(i)
1784	$ + 4.0RSIGMAztsol(i)*3 (ztsol(i)-ftsol(i,nsrf))
1785	fsolsw(i,nsrf) = solsw(i)*(1.-falbe(i,nsrf))/(1.-albsol(i))
1786	ENDDO
1787	ENDDO
1788
1789	fder = dlw
1790
1791
1792	CALL clmain(dtime,itap,date0,pctsrf,pctsrf_new,
1793	e t_seri,q_seri,u_seri,v_seri,
1794	e julien, rmu0, co2_ppm,
1795	e ok_veget, ocean, npas, nexca, ftsol,
1796	$ soil_model,cdmmax, cdhmax,
1797	$ ksta, ksta_ter, ok_kzmin, ftsoil, qsol,
1798	cIM BAD $ paprs,pplay,radsol, fsnow,fqsurf,fevap,falbe,falblw,
1799	$ paprs,pplay, fsnow,fqsurf,fevap,falbe,falblw,
1800	$ fluxlat,
1801	e rain_fall, snow_fall,
1802	e fsolsw, fsollw, sollwdown, fder,
1803	e rlon, rlat, cuphy, cvphy, frugs,
1804	e debut, lafin, agesno,rugoro ,
1805	s d_t_vdf,d_q_vdf,d_u_vdf,d_v_vdf,d_ts,
1806	s fluxt,fluxq,fluxu,fluxv,cdragh,cdragm,
1807	s q2,
1808	s dsens, devap,
1809	s ycoefh,yu1,yv1, t2m, q2m, u10m, v10m,
1810	cIM cf. AM 081204 BEG
1811	s pblh,capCL,oliqCL,cteiCL,pblT,
1812	s therm,trmb1,trmb2,trmb3,plcl,
1813	cIM cf. AM 081204 END
1814	s fqcalving, ffonte, run_off_lic_0,
1815	cIM "slab" ocean
1816	s fluxo, fluxg, tslab, seaice)
1817	c
1818	CXXX PB
1819	CXXX Incrementation des flux
1820	CXXX
1821
1822	zxfluxt=0.
1823	zxfluxq=0.
1824	zxfluxu=0.
1825	zxfluxv=0.
1826	DO nsrf = 1, nbsrf
1827	DO k = 1, klev
1828	DO i = 1, klon
1829	zxfluxt(i,k) = zxfluxt(i,k) +
1830	$ fluxt(i,k,nsrf) * pctsrf( i, nsrf)
1831	zxfluxq(i,k) = zxfluxq(i,k) +
1832	$ fluxq(i,k,nsrf) * pctsrf( i, nsrf)
1833	zxfluxu(i,k) = zxfluxu(i,k) +
1834	$ fluxu(i,k,nsrf) * pctsrf( i, nsrf)
1835	zxfluxv(i,k) = zxfluxv(i,k) +
1836	$ fluxv(i,k,nsrf) * pctsrf( i, nsrf)
1837	END DO
1838	END DO
1839	END DO
1840	DO i = 1, klon
1841	sens(i) = - zxfluxt(i,1) ! flux de chaleur sensible au sol
1842	c evap(i) = - fluxq(i,1) ! flux d'evaporation au sol
1843	evap(i) = - zxfluxq(i,1) ! flux d'evaporation au sol
1844	fder(i) = dlw(i) + dsens(i) + devap(i)
1845	ENDDO
1846
1847
1848	DO k = 1, klev
1849	DO i = 1, klon
1850	t_seri(i,k) = t_seri(i,k) + d_t_vdf(i,k)
1851	q_seri(i,k) = q_seri(i,k) + d_q_vdf(i,k)
1852	u_seri(i,k) = u_seri(i,k) + d_u_vdf(i,k)
1853	v_seri(i,k) = v_seri(i,k) + d_v_vdf(i,k)
1854	ENDDO
1855	ENDDO
1856	c
1857	IF (if_ebil.ge.2) THEN
1858	ztit='after clmain'
1859	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
1860	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1861	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1862	call diagphy(airephy,ztit,ip_ebil
1863	e , zero_v, zero_v, zero_v, zero_v, sens
1864	e , evap , zero_v, zero_v, ztsol
1865	e , d_h_vcol, d_qt, d_ec
1866	s , fs_bound, fq_bound )
1867	END IF
1868	C
1869	c
1870	c Incrementer la temperature du sol
1871	c
1872	DO i = 1, klon
1873	zxtsol(i) = 0.0
1874	zxfluxlat(i) = 0.0
1875	c
1876	zt2m(i) = 0.0
1877	zq2m(i) = 0.0
1878	zu10m(i) = 0.0
1879	zv10m(i) = 0.0
1880	cIM cf JLD ??
1881	zxffonte(i) = 0.0
1882	zxfqcalving(i) = 0.0
1883	cIM cf. AM 081204 BEG
1884	c
1885	s_pblh(i) = 0.0
1886	s_lcl(i) = 0.0
1887	s_capCL(i) = 0.0
1888	s_oliqCL(i) = 0.0
1889	s_cteiCL(i) = 0.0
1890	s_pblT(i) = 0.0
1891	s_therm(i) = 0.0
1892	s_trmb1(i) = 0.0
1893	s_trmb2(i) = 0.0
1894	s_trmb3(i) = 0.0
1895	c
1896	IF ( abs( pctsrf(i, is_ter) + pctsrf(i, is_lic) +
1897	$ pctsrf(i, is_oce) + pctsrf(i, is_sic) - 1.) .GT. EPSFRA)
1898	$ THEN
1899	WRITE(,) 'physiq : pb sous surface au point ', i,
1900	$ pctsrf(i, 1 : nbsrf)
1901	ENDIF
1902	ENDDO
1903	DO nsrf = 1, nbsrf
1904	DO i = 1, klon
1905	c IF (pctsrf(i,nsrf) .GE. EPSFRA) THEN
1906	ftsol(i,nsrf) = ftsol(i,nsrf) + d_ts(i,nsrf)
1907	cIM cf. JLD
1908	wfbils(i,nsrf) = ( fsolsw(i,nsrf) + fsollw(i,nsrf)
1909	$ + fluxt(i,1,nsrf) + fluxlat(i,nsrf) ) * pctsrf(i,nsrf)
1910	zxtsol(i) = zxtsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1911	zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf)*pctsrf(i,nsrf)
1912	cccIM
1913	zt2m(i) = zt2m(i) + t2m(i,nsrf)*pctsrf(i,nsrf)
1914	zq2m(i) = zq2m(i) + q2m(i,nsrf)*pctsrf(i,nsrf)
1915	zu10m(i) = zu10m(i) + u10m(i,nsrf)*pctsrf(i,nsrf)
1916	zv10m(i) = zv10m(i) + v10m(i,nsrf)*pctsrf(i,nsrf)
1917	cIM cf JLD ??
1918	zxffonte(i) = zxffonte(i) + ffonte(i,nsrf)*pctsrf(i,nsrf)
1919	zxfqcalving(i) = zxfqcalving(i) +
1920	. fqcalving(i,nsrf)*pctsrf(i,nsrf)
1921	cIM cf. AM 081204 BEG
1922	s_pblh(i) = s_pblh(i) + pblh(i,nsrf)*pctsrf(i,nsrf)
1923	s_lcl(i) = s_lcl(i) + plcl(i,nsrf)*pctsrf(i,nsrf)
1924	s_capCL(i) = s_capCL(i) + capCL(i,nsrf) *pctsrf(i,nsrf)
1925	s_oliqCL(i) = s_oliqCL(i) + oliqCL(i,nsrf) *pctsrf(i,nsrf)
1926	s_cteiCL(i) = s_cteiCL(i) + cteiCL(i,nsrf) *pctsrf(i,nsrf)
1927	s_pblT(i) = s_pblT(i) + pblT(i,nsrf) *pctsrf(i,nsrf)
1928	s_therm(i) = s_therm(i) + therm(i,nsrf) *pctsrf(i,nsrf)
1929	s_trmb1(i) = s_trmb1(i) + trmb1(i,nsrf) *pctsrf(i,nsrf)
1930	s_trmb2(i) = s_trmb2(i) + trmb2(i,nsrf) *pctsrf(i,nsrf)
1931	s_trmb3(i) = s_trmb3(i) + trmb3(i,nsrf) *pctsrf(i,nsrf)
1932	c ENDIF
1933	ENDDO
1934	ENDDO
1935
1936	c
1937	c Si une sous-fraction n'existe pas, elle prend la temp. moyenne
1938	c
1939	DO nsrf = 1, nbsrf
1940	DO i = 1, klon
1941	IF (pctsrf(i,nsrf) .LT. epsfra) ftsol(i,nsrf) = zxtsol(i)
1942	cccIM
1943	IF (pctsrf(i,nsrf) .LT. epsfra) t2m(i,nsrf) = zt2m(i)
1944	IF (pctsrf(i,nsrf) .LT. epsfra) q2m(i,nsrf) = zq2m(i)
1945	IF (pctsrf(i,nsrf) .LT. epsfra) u10m(i,nsrf) = zu10m(i)
1946	IF (pctsrf(i,nsrf) .LT. epsfra) v10m(i,nsrf) = zv10m(i)
1947	cIM cf JLD ??
1948	IF (pctsrf(i,nsrf) .LT. epsfra) ffonte(i,nsrf) = zxffonte(i)
1949	IF (pctsrf(i,nsrf) .LT. epsfra)
1950	. fqcalving(i,nsrf) = zxfqcalving(i)
1951	cIM cf. AM 081204 BEG
1952	IF (pctsrf(i,nsrf) .LT. epsfra) pblh(i,nsrf)=s_pblh(i)
1953	IF (pctsrf(i,nsrf) .LT. epsfra) plcl(i,nsrf)=s_lcl(i)
1954	IF (pctsrf(i,nsrf) .LT. epsfra) capCL(i,nsrf)=s_capCL(i)
1955	IF (pctsrf(i,nsrf) .LT. epsfra) oliqCL(i,nsrf)=s_oliqCL(i)
1956	IF (pctsrf(i,nsrf) .LT. epsfra) cteiCL(i,nsrf)=s_cteiCL(i)
1957	IF (pctsrf(i,nsrf) .LT. epsfra) pblT(i,nsrf)=s_pblT(i)
1958	IF (pctsrf(i,nsrf) .LT. epsfra) therm(i,nsrf)=s_therm(i)
1959	IF (pctsrf(i,nsrf) .LT. epsfra) trmb1(i,nsrf)=s_trmb1(i)
1960	IF (pctsrf(i,nsrf) .LT. epsfra) trmb2(i,nsrf)=s_trmb2(i)
1961	IF (pctsrf(i,nsrf) .LT. epsfra) trmb3(i,nsrf)=s_trmb3(i)
1962	ENDDO
1963	ENDDO
1964	c
1965	c
1966	c Calculer la derive du flux infrarouge
1967	c
1968	cXXX DO nsrf = 1, nbsrf
1969	DO i = 1, klon
1970	cXXX IF (pctsrf(i,nsrf) .GE. EPSFRA) THEN
1971	dlw(i) = - 4.0RSIGMAzxtsol(i)**3
1972	cXXX . *(ftsol(i,nsrf)-zxtsol(i))
1973	cXXX . *pctsrf(i,nsrf)
1974	cXXX ENDIF
1975	cXXX ENDDO
1976	ENDDO
1977	c
1978	c Appeler la convection (au choix)
1979	c
1980	DO k = 1, klev
1981	DO i = 1, klon
1982	conv_q(i,k) = d_q_dyn(i,k)
1983	. + d_q_vdf(i,k)/dtime
1984	conv_t(i,k) = d_t_dyn(i,k)
1985	. + d_t_vdf(i,k)/dtime
1986	ENDDO
1987	ENDDO
1988	IF (check) THEN
1989	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
1990	WRITE(lunout,*) "avantcon=", za
1991	ENDIF
1992	zx_ajustq = .FALSE.
1993	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
1994	IF (zx_ajustq) THEN
1995	DO i = 1, klon
1996	z_avant(i) = 0.0
1997	ENDDO
1998	DO k = 1, klev
1999	DO i = 1, klon
2000	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k))
2001	. *(paprs(i,k)-paprs(i,k+1))/RG
2002	ENDDO
2003	ENDDO
2004	ENDIF
2005	IF (iflag_con.EQ.1) THEN
2006	stop'reactiver le call conlmd dans physiq.F'
2007	c CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
2008	c . d_t_con, d_q_con,
2009	c . rain_con, snow_con, ibas_con, itop_con)
2010	ELSE IF (iflag_con.EQ.2) THEN
2011	CALL conflx(dtime, paprs, pplay, t_seri, q_seri,
2012	e conv_t, conv_q, zxfluxq(1,1), omega,
2013	s d_t_con, d_q_con, rain_con, snow_con,
2014	s pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
2015	s kcbot, kctop, kdtop, pmflxr, pmflxs)
2016	WHERE (rain_con < 0.) rain_con = 0.
2017	WHERE (snow_con < 0.) snow_con = 0.
2018	DO i = 1, klon
2019	ibas_con(i) = klev+1 - kcbot(i)
2020	itop_con(i) = klev+1 - kctop(i)
2021	ENDDO
2022	ELSE IF (iflag_con.GE.3) THEN
2023	c nb of tracers for the KE convection:
2024	c MAF la partie traceurs est faite dans phytrac
2025	c on met ntra=1 pour limiter les appels mais on peut
2026	c supprimer les calculs / ftra.
2027	ntra = 1
2028	c sb, oct02:
2029	c Schema de convection modularise et vectorise:
2030	c (driver commun aux versions 3 et 4)
2031	c
2032	IF (ok_cvl) THEN ! new driver for convectL
2033
2034	CALL concvl (iflag_con,
2035	. dtime,paprs,pplay,t_seri,q_seri,
2036	. u_seri,v_seri,tr_seri,ntra,
2037	. ema_work1,ema_work2,
2038	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2039	. rain_con, snow_con, ibas_con, itop_con,
2040	. upwd,dnwd,dnwd0,
2041	. Ma,cape,tvp,iflagctrl,
2042	. pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd,
2043	. pmflxr,pmflxs,
2044	. da,phi,mp)
2045
2046	cIM cf. FH
2047	clwcon0=qcondc
2048	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
2049
2050	ELSE ! ok_cvl
2051	c MAF conema3 ne contient pas les traceurs
2052	CALL conema3 (dtime,
2053	. paprs,pplay,t_seri,q_seri,
2054	. u_seri,v_seri,tr_seri,ntra,
2055	. ema_work1,ema_work2,
2056	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
2057	. rain_con, snow_con, ibas_con, itop_con,
2058	. upwd,dnwd,dnwd0,bas,top,
2059	. Ma,cape,tvp,rflag,
2060	. pbase
2061	. ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr
2062	. ,clwcon0)
2063
2064	ENDIF ! ok_cvl
2065
2066	IF (.NOT. ok_gust) THEN
2067	do i = 1, klon
2068	wd(i)=0.0
2069	enddo
2070	ENDIF
2071
2072	c =================================================================== c
2073	c Calcul des proprietes des nuages convectifs
2074	c
2075	DO k = 1, klev
2076	DO i = 1, klon
2077	zx_t = t_seri(i,k)
2078	IF (thermcep) THEN
2079	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2080	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2081	zx_qs = MIN(0.5,zx_qs)
2082	zcor = 1./(1.-retv*zx_qs)
2083	zx_qs = zx_qs*zcor
2084	ELSE
2085	IF (zx_t.LT.t_coup) THEN
2086	zx_qs = qsats(zx_t)/pplay(i,k)
2087	ELSE
2088	zx_qs = qsatl(zx_t)/pplay(i,k)
2089	ENDIF
2090	ENDIF
2091	zqsat(i,k)=zx_qs
2092	ENDDO
2093	ENDDO
2094
2095	c calcul des proprietes des nuages convectifs
2096	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
2097	call clouds_gno
2098	s (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
2099
2100	c =================================================================== c
2101
2102	DO i = 1, klon
2103	ema_pcb(i) = pbase(i)
2104	ENDDO
2105	DO i = 1, klon
2106	ema_pct(i) = paprs(i,itop_con(i))
2107	ENDDO
2108	DO i = 1, klon
2109	ema_cbmf(i) = ema_workcbmf(i)
2110	ENDDO
2111	ELSE
2112	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
2113	CALL abort
2114	ENDIF
2115
2116	c CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
2117	c . d_u_con, d_v_con)
2118
2119	DO k = 1, klev
2120	DO i = 1, klon
2121	t_seri(i,k) = t_seri(i,k) + d_t_con(i,k)
2122	q_seri(i,k) = q_seri(i,k) + d_q_con(i,k)
2123	u_seri(i,k) = u_seri(i,k) + d_u_con(i,k)
2124	v_seri(i,k) = v_seri(i,k) + d_v_con(i,k)
2125	ENDDO
2126	ENDDO
2127	c
2128	IF (if_ebil.ge.2) THEN
2129	ztit='after convect'
2130	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
2131	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2132	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2133	call diagphy(airephy,ztit,ip_ebil
2134	e , zero_v, zero_v, zero_v, zero_v, zero_v
2135	e , zero_v, rain_con, snow_con, ztsol
2136	e , d_h_vcol, d_qt, d_ec
2137	s , fs_bound, fq_bound )
2138	END IF
2139	C
2140	IF (check) THEN
2141	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2142	WRITE(lunout,*)"aprescon=", za
2143	zx_t = 0.0
2144	za = 0.0
2145	DO i = 1, klon
2146	za = za + airephy(i)/FLOAT(klon)
2147	zx_t = zx_t + (rain_con(i)+
2148	. snow_con(i))*airephy(i)/FLOAT(klon)
2149	ENDDO
2150	zx_t = zx_t/za*dtime
2151	WRITE(lunout,*)"Precip=", zx_t
2152	ENDIF
2153	IF (zx_ajustq) THEN
2154	DO i = 1, klon
2155	z_apres(i) = 0.0
2156	ENDDO
2157	DO k = 1, klev
2158	DO i = 1, klon
2159	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k))
2160	. *(paprs(i,k)-paprs(i,k+1))/RG
2161	ENDDO
2162	ENDDO
2163	DO i = 1, klon
2164	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime)
2165	. /z_apres(i)
2166	ENDDO
2167	DO k = 1, klev
2168	DO i = 1, klon
2169	IF (z_factor(i).GT.(1.0+1.0E-08) .OR.
2170	. z_factor(i).LT.(1.0-1.0E-08)) THEN
2171	q_seri(i,k) = q_seri(i,k) * z_factor(i)
2172	ENDIF
2173	ENDDO
2174	ENDDO
2175	ENDIF
2176	zx_ajustq=.FALSE.
2177	c
2178	c===================================================================
2179	c Convection seche (thermiques ou ajustement)
2180	c===================================================================
2181	c
2182	d_t_ajs(:,:)=0.
2183	d_u_ajs(:,:)=0.
2184	d_v_ajs(:,:)=0.
2185	d_q_ajs(:,:)=0.
2186	fm_therm(:,:)=0.
2187	entr_therm(:,:)=0.
2188	c
2189	IF(prt_level>9)WRITE(lunout,*)
2190	. 'AVANT LA CONVECTION SECHE , iflag_thermals='
2191	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2192	if(iflag_thermals.lt.0) then
2193	c Rien
2194	c ====
2195	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
2196	else if(iflag_thermals.eq.0) then
2197
2198	c Ajustement sec
2199	c ==============
2200	IF(prt_level>9)WRITE(lunout,*)'ajsec'
2201	CALL ajsec(paprs, pplay, t_seri,q_seri, d_t_ajs, d_q_ajs)
2202	t_seri(:,:) = t_seri(:,:) + d_t_ajs(:,:)
2203	q_seri(:,:) = q_seri(:,:) + d_q_ajs(:,:)
2204	else
2205	c Thermiques
2206	c ==========
2207	IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals='
2208	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2209	call calltherm(pdtphys
2210	s ,pplay,paprs,pphi
2211	s ,u_seri,v_seri,t_seri,q_seri
2212	s ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
2213	s ,fm_therm,entr_therm)
2214	endif
2215	c
2216	c===================================================================
2217	c
2218	IF (if_ebil.ge.2) THEN
2219	ztit='after dry_adjust'
2220	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
2221	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2222	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2223	END IF
2224
2225
2226	c-------------------------------------------------------------------------
2227	c Caclul des ratqs
2228	c-------------------------------------------------------------------------
2229
2230	c print*,'calcul des ratqs'
2231	c ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
2232	c ----------------
2233	c on ecrase le tableau ratqsc calcule par clouds_gno
2234	if (iflag_cldcon.eq.1) then
2235	do k=1,klev
2236	do i=1,klon
2237	if(ptconv(i,k)) then
2238	ratqsc(i,k)=ratqsbas
2239	s +fact_cldcon*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
2240	else
2241	ratqsc(i,k)=0.
2242	endif
2243	enddo
2244	enddo
2245	endif
2246
2247	c ratqs stables
2248	c -------------
2249	do k=1,klev
2250	cIM RAJOUT boucle do=i
2251	do i=1, klon
2252	cIM ratqss(:,k)=ratqsbas+(ratqshaut-ratqsbas)*
2253	cIM s min((paprs(:,1)-pplay(:,k))/(paprs(:,1)-30000.),1.)
2254	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
2255	s min((paprs(i,1)-pplay(i,k))/(paprs(i,1)-30000.),1.)
2256	cIM print*,' IMratqs STABLE i, k',i,k,ratqss(i,k)
2257	enddo
2258	enddo
2259
2260
2261	c ratqs final
2262	c -----------
2263	if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2) then
2264	c les ratqs sont une conbinaison de ratqss et ratqsc
2265	c ratqs final
2266	c 1e4 (en gros 3 heures), en dur pour le moment, est le temps de
2267	c relaxation des ratqs
2268	c facttemps=exp(-pdtphys/1.e4)
2269	facteur=exp(-pdtphys*facttemps)
2270	ratqs(:,:)=max(ratqs(:,:)*facteur,ratqss(:,:))
2271	ratqs(:,:)=max(ratqs(:,:),ratqsc(:,:))
2272	c print*,'calcul des ratqs fini'
2273	else
2274	c on ne prend que le ratqs stable pour fisrtilp
2275	ratqs(:,:)=ratqss(:,:)
2276	endif
2277
2278
2279	c
2280	c Appeler le processus de condensation a grande echelle
2281	c et le processus de precipitation
2282	c-------------------------------------------------------------------------
2283	CALL fisrtilp(dtime,paprs,pplay,
2284	. t_seri, q_seri,ptconv,ratqs,
2285	. d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq,
2286	. rain_lsc, snow_lsc,
2287	. pfrac_impa, pfrac_nucl, pfrac_1nucl,
2288	. frac_impa, frac_nucl,
2289	. prfl, psfl, rhcl)
2290
2291	WHERE (rain_lsc < 0) rain_lsc = 0.
2292	WHERE (snow_lsc < 0) snow_lsc = 0.
2293	DO k = 1, klev
2294	DO i = 1, klon
2295	t_seri(i,k) = t_seri(i,k) + d_t_lsc(i,k)
2296	q_seri(i,k) = q_seri(i,k) + d_q_lsc(i,k)
2297	ql_seri(i,k) = ql_seri(i,k) + d_ql_lsc(i,k)
2298	cldfra(i,k) = rneb(i,k)
2299	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
2300	ENDDO
2301	ENDDO
2302	IF (check) THEN
2303	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2304	WRITE(lunout,*)"apresilp=", za
2305	zx_t = 0.0
2306	za = 0.0
2307	DO i = 1, klon
2308	za = za + airephy(i)/FLOAT(klon)
2309	zx_t = zx_t + (rain_lsc(i)
2310	. + snow_lsc(i))*airephy(i)/FLOAT(klon)
2311	ENDDO
2312	zx_t = zx_t/za*dtime
2313	WRITE(lunout,*)"Precip=", zx_t
2314	ENDIF
2315	c
2316	IF (if_ebil.ge.2) THEN
2317	ztit='after fisrt'
2318	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
2319	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2320	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2321	call diagphy(airephy,ztit,ip_ebil
2322	e , zero_v, zero_v, zero_v, zero_v, zero_v
2323	e , zero_v, rain_lsc, snow_lsc, ztsol
2324	e , d_h_vcol, d_qt, d_ec
2325	s , fs_bound, fq_bound )
2326	END IF
2327	c
2328	c-------------------------------------------------------------------
2329	c PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
2330	c-------------------------------------------------------------------
2331
2332	c 1. NUAGES CONVECTIFS
2333	c
2334	cIM cf FH
2335	c IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
2336	IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
2337	snow_tiedtke=0.
2338	c print*,'avant calcul de la pseudo precip '
2339	c print*,'iflag_cldcon',iflag_cldcon
2340	if (iflag_cldcon.eq.-1) then
2341	rain_tiedtke=rain_con
2342	else
2343	c print*,'calcul de la pseudo precip '
2344	rain_tiedtke=0.
2345	c print*,'calcul de la pseudo precip 0'
2346	do k=1,klev
2347	do i=1,klon
2348	if (d_q_con(i,k).lt.0.) then
2349	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys
2350	s *(paprs(i,k)-paprs(i,k+1))/rg
2351	endif
2352	enddo
2353	enddo
2354	endif
2355	c
2356	c call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
2357	c
2358
2359	c Nuages diagnostiques pour Tiedtke
2360	CALL diagcld1(paprs,pplay,
2361	cIM cf FH . rain_con,snow_con,ibas_con,itop_con,
2362	. rain_tiedtke,snow_tiedtke,ibas_con,itop_con,
2363	. diafra,dialiq)
2364	DO k = 1, klev
2365	DO i = 1, klon
2366	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2367	cldliq(i,k) = dialiq(i,k)
2368	cldfra(i,k) = diafra(i,k)
2369	ENDIF
2370	ENDDO
2371	ENDDO
2372
2373	ELSE IF (iflag_cldcon.eq.3) THEN
2374	c On prend pour les nuages convectifs le max du calcul de la
2375	c convection et du calcul du pas de temps précédent diminué d'un facteur
2376	c facttemps
2377	c facttemps=pdtphys/1.e4
2378	facteur = pdtphys *facttemps
2379	do k=1,klev
2380	do i=1,klon
2381	rnebcon(i,k)=rnebcon(i,k)*facteur
2382	if (rnebcon0(i,k)clwcon0(i,k).gt.rnebcon(i,k)clwcon(i,k))
2383	s then
2384	rnebcon(i,k)=rnebcon0(i,k)
2385	clwcon(i,k)=clwcon0(i,k)
2386	endif
2387	enddo
2388	enddo
2389
2390	c
2391	cIM calcul nuages par le simulateur ISCCP
2392	c
2393	IF (ok_isccp) THEN
2394	#include "calcul_simulISCCP.h"
2395	ENDIF !ok_isccp
2396
2397	c On prend la somme des fractions nuageuses et des contenus en eau
2398	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
2399	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
2400
2401	ENDIF
2402
2403	c
2404	c 2. NUAGES STARTIFORMES
2405	c
2406	IF (ok_stratus) THEN
2407	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
2408	DO k = 1, klev
2409	DO i = 1, klon
2410	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2411	cldliq(i,k) = dialiq(i,k)
2412	cldfra(i,k) = diafra(i,k)
2413	ENDIF
2414	ENDDO
2415	ENDDO
2416	ENDIF
2417	c
2418	c Precipitation totale
2419	c
2420	DO i = 1, klon
2421	rain_fall(i) = rain_con(i) + rain_lsc(i)
2422	snow_fall(i) = snow_con(i) + snow_lsc(i)
2423	ENDDO
2424	c
2425	IF (if_ebil.ge.2) THEN
2426	ztit="after diagcld"
2427	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
2428	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2429	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2430	END IF
2431	c
2432	c Calculer l'humidite relative pour diagnostique
2433	c
2434	DO k = 1, klev
2435	DO i = 1, klon
2436	zx_t = t_seri(i,k)
2437	IF (thermcep) THEN
2438	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2439	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2440	zx_qs = MIN(0.5,zx_qs)
2441	zcor = 1./(1.-retv*zx_qs)
2442	zx_qs = zx_qs*zcor
2443	ELSE
2444	IF (zx_t.LT.t_coup) THEN
2445	zx_qs = qsats(zx_t)/pplay(i,k)
2446	ELSE
2447	zx_qs = qsatl(zx_t)/pplay(i,k)
2448	ENDIF
2449	ENDIF
2450	zx_rh(i,k) = q_seri(i,k)/zx_qs
2451	zqsat(i,k)=zx_qs
2452	ENDDO
2453	ENDDO
2454	cjq - introduce the aerosol direct and first indirect radiative forcings
2455	cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
2456	IF (ok_ade.OR.ok_aie) THEN
2457	! Get sulfate aerosol distribution
2458	CALL readsulfate(rjourvrai, debut, sulfate)
2459	CALL readsulfate_preind(rjourvrai, debut, sulfate_pi)
2460
2461	! Calculate aerosol optical properties (Olivier Boucher)
2462	CALL aeropt(pplay, paprs, t_seri, sulfate, rhcl,
2463	. tau_ae, piz_ae, cg_ae, aerindex)
2464	cym
2465	ELSE
2466	tau_ae(:,:,:)=0.0
2467	piz_ae(:,:,:)=0.0
2468	cg_ae(:,:,:)=0.0
2469	cym
2470	ENDIF
2471
2472	#ifdef INCA
2473	calday = FLOAT(julien) + gmtime
2474
2475	#ifdef INCA_AER
2476	call AEROSOL_METEO_CALC(calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs,
2477	& prfl,psfl,pctsrf(1,3),airephy,xjour,rlat,rlon)
2478	#endif
2479
2480	#ifdef INCAINFO
2481	WRITE(lunout,*)'Appel CHEMHOOK_BEGIN ...'
2482	#endif
2483
2484	CALL chemhook_begin (calday,
2485	$ pctsrf(1,1),
2486	$ rlat,
2487	$ rlon,
2488	$ airephy,
2489	$ paprs,
2490	$ pplay,
2491	$ ycoefh,
2492	$ pphi,
2493	$ t_seri,
2494	$ u,
2495	$ v,
2496	$ wo,
2497	$ q_seri,
2498	$ zxtsol,
2499	$ zxsnow,
2500	$ solsw,
2501	$ albsol,
2502	$ rain_fall,
2503	$ snow_fall,
2504	$ itop_con,
2505	$ ibas_con,
2506	$ cldfra,
2507	$ iim,
2508	$ jjm,
2509	#ifdef INCA_AER
2510	$ tr_seri,
2511	$ ftsol,
2512	$ paprs,
2513	$ cdragh,
2514	$ cdragm,
2515	$ pctsrf,
2516	$ pdtphys,
2517	$ itap)
2518	#else
2519	$ tr_seri)
2520	#endif
2521
2522
2523	#ifdef INCAINFO
2524	WRITE(lunout,*)'OK.'
2525	#endif
2526	#endif
2527	c
2528	c Calculer les parametres optiques des nuages et quelques
2529	c parametres pour diagnostiques:
2530	c
2531	if (ok_newmicro) then
2532	CALL newmicro (paprs, pplay,ok_newmicro,
2533	. t_seri, cldliq, cldfra, cldtau, cldemi,
2534	. cldh, cldl, cldm, cldt, cldq,
2535	. flwp, fiwp, flwc, fiwc,
2536	e ok_aie,
2537	e sulfate, sulfate_pi,
2538	e bl95_b0, bl95_b1,
2539	s cldtaupi, re, fl)
2540	else
2541	CALL nuage (paprs, pplay,
2542	. t_seri, cldliq, cldfra, cldtau, cldemi,
2543	. cldh, cldl, cldm, cldt, cldq,
2544	e ok_aie,
2545	e sulfate, sulfate_pi,
2546	e bl95_b0, bl95_b1,
2547	s cldtaupi, re, fl)
2548
2549	endif
2550	c
2551	c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
2552	c
2553	IF (MOD(itaprad,radpas).EQ.0) THEN
2554	DO i = 1, klon
2555	albsol(i) = falbe(i,is_oce) * pctsrf(i,is_oce)
2556	. + falbe(i,is_lic) * pctsrf(i,is_lic)
2557	. + falbe(i,is_ter) * pctsrf(i,is_ter)
2558	. + falbe(i,is_sic) * pctsrf(i,is_sic)
2559	albsollw(i) = falblw(i,is_oce) * pctsrf(i,is_oce)
2560	. + falblw(i,is_lic) * pctsrf(i,is_lic)
2561	. + falblw(i,is_ter) * pctsrf(i,is_ter)
2562	. + falblw(i,is_sic) * pctsrf(i,is_sic)
2563	ENDDO
2564	CALL radlwsw ! nouveau rayonnement (compatible Arpege-IFS)
2565	e (dist, rmu0, fract,
2566	e paprs, pplay,zxtsol,albsol, albsollw, t_seri,q_seri,
2567	e wo,
2568	e cldfra, cldemi, cldtau,
2569	s heat,heat0,cool,cool0,radsol,albpla,
2570	s topsw,toplw,solsw,sollw,
2571	s sollwdown,
2572	s topsw0,toplw0,solsw0,sollw0,
2573	s lwdn0, lwdn, lwup0, lwup,
2574	s swdn0, swdn, swup0, swup,
2575	e ok_ade, ok_aie, ! new for aerosol radiative effects
2576	e tau_ae, piz_ae, cg_ae, ! ="=
2577	s topswad, solswad, ! ="=
2578	e cldtaupi, ! ="=
2579	s topswai, solswai) ! ="=
2580	itaprad = 0
2581	ENDIF
2582	itaprad = itaprad + 1
2583	c
2584	c Ajouter la tendance des rayonnements (tous les pas)
2585	c
2586	DO k = 1, klev
2587	DO i = 1, klon
2588	t_seri(i,k) = t_seri(i,k)
2589	. + (heat(i,k)-cool(i,k)) * dtime/86400.
2590	ENDDO
2591	ENDDO
2592	c
2593	IF (if_ebil.ge.2) THEN
2594	ztit='after rad'
2595	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
2596	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2597	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2598	call diagphy(airephy,ztit,ip_ebil
2599	e , topsw, toplw, solsw, sollw, zero_v
2600	e , zero_v, zero_v, zero_v, ztsol
2601	e , d_h_vcol, d_qt, d_ec
2602	s , fs_bound, fq_bound )
2603	END IF
2604	c
2605	c
2606	c Calculer l'hydrologie de la surface
2607	c
2608	c CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
2609	c . agesno, ftsol,fqsurf,fsnow, ruis)
2610	c
2611	DO i = 1, klon
2612	zxqsurf(i) = 0.0
2613	zxsnow(i) = 0.0
2614	ENDDO
2615	DO nsrf = 1, nbsrf
2616	DO i = 1, klon
2617	zxqsurf(i) = zxqsurf(i) + fqsurf(i,nsrf)*pctsrf(i,nsrf)
2618	zxsnow(i) = zxsnow(i) + fsnow(i,nsrf)*pctsrf(i,nsrf)
2619	ENDDO
2620	ENDDO
2621	c
2622	c Si une sous-fraction n'existe pas, elle prend la valeur moyenne
2623	c
2624	cXXX DO nsrf = 1, nbsrf
2625	cXXX DO i = 1, klon
2626	cXXX IF (pctsrf(i,nsrf).LT.epsfra) THEN
2627	cXXX fqsurf(i,nsrf) = zxqsurf(i)
2628	cXXX fsnow(i,nsrf) = zxsnow(i)
2629	cXXX ENDIF
2630	cXXX ENDDO
2631	cXXX ENDDO
2632	c
2633	c Calculer le bilan du sol et la derive de temperature (couplage)
2634	c
2635	DO i = 1, klon
2636	c bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
2637	c a la demande de JLD
2638	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
2639	ENDDO
2640	c
2641	cmoddeblott(jan95)
2642	c Appeler le programme de parametrisation de l'orographie
2643	c a l'echelle sous-maille:
2644	c
2645	IF (ok_orodr) THEN
2646	c
2647	c selection des points pour lesquels le shema est actif:
2648	igwd=0
2649	DO i=1,klon
2650	itest(i)=0
2651	c IF ((zstd(i).gt.10.0)) THEN
2652	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
2653	itest(i)=1
2654	igwd=igwd+1
2655	idx(igwd)=i
2656	ENDIF
2657	ENDDO
2658	c igwdim=MAX(1,igwd)
2659	c
2660	CALL drag_noro(klon,klev,dtime,paprs,pplay,
2661	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
2662	e igwd,idx,itest,
2663	e t_seri, u_seri, v_seri,
2664	cIM 141004 s zulow, zvlow, zustr, zvstr,
2665	s zulow, zvlow, zustrdr, zvstrdr,
2666	s d_t_oro, d_u_oro, d_v_oro)
2667	c
2668	c ajout des tendances
2669	DO k = 1, klev
2670	DO i = 1, klon
2671	t_seri(i,k) = t_seri(i,k) + d_t_oro(i,k)
2672	u_seri(i,k) = u_seri(i,k) + d_u_oro(i,k)
2673	v_seri(i,k) = v_seri(i,k) + d_v_oro(i,k)
2674	ENDDO
2675	ENDDO
2676	c
2677	ENDIF ! fin de test sur ok_orodr
2678	c
2679	IF (ok_orolf) THEN
2680	c
2681	c selection des points pour lesquels le shema est actif:
2682	igwd=0
2683	DO i=1,klon
2684	itest(i)=0
2685	IF ((zpic(i)-zmea(i)).GT.100.) THEN
2686	itest(i)=1
2687	igwd=igwd+1
2688	idx(igwd)=i
2689	ENDIF
2690	ENDDO
2691	c igwdim=MAX(1,igwd)
2692	c
2693	CALL lift_noro(klon,klev,dtime,paprs,pplay,
2694	e rlat,zmea,zstd,zpic,
2695	e itest,
2696	e t_seri, u_seri, v_seri,
2697	s zulow, zvlow, zustrli, zvstrli,
2698	s d_t_lif, d_u_lif, d_v_lif)
2699	c
2700	c ajout des tendances
2701	DO k = 1, klev
2702	DO i = 1, klon
2703	t_seri(i,k) = t_seri(i,k) + d_t_lif(i,k)
2704	u_seri(i,k) = u_seri(i,k) + d_u_lif(i,k)
2705	v_seri(i,k) = v_seri(i,k) + d_v_lif(i,k)
2706	ENDDO
2707	ENDDO
2708	c
2709	ENDIF ! fin de test sur ok_orolf
2710	c
2711	cIM cf. FLott BEG
2712	C STRESS NECESSAIRES: TOUTE LA PHYSIQUE
2713
2714	DO i = 1, klon
2715	zustrph(i)=0.
2716	zvstrph(i)=0.
2717	ENDDO
2718	DO k = 1, klev
2719	DO i = 1, klon
2720	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
2721	c (paprs(i,k)-paprs(i,k+1))/rg
2722	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
2723	c (paprs(i,k)-paprs(i,k+1))/rg
2724	ENDDO
2725	ENDDO
2726	c
2727	cIM calcul composantes axiales du moment angulaire et couple des montagnes
2728	c
2729	CALL aaam_bud (27,klon,klev,rjourvrai,gmtime,
2730	C ra,rg,romega,
2731	C rlat,rlon,pphis,
2732	C zustrdr,zustrli,zustrph,
2733	C zvstrdr,zvstrli,zvstrph,
2734	C paprs,u,v,
2735	C aam, torsfc)
2736	cIM cf. FLott END
2737	c
2738	IF (if_ebil.ge.2) THEN
2739	ztit='after orography'
2740	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
2741	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2742	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2743	END IF
2744	c
2745	c
2746	cAA
2747	cAA Installation de l'interface online-offline pour traceurs
2748	cAA
2749	c====================================================================
2750	c Calcul des tendances traceurs
2751	c====================================================================
2752	C
2753	call phytrac (iflag_con,
2754	I rnpb,
2755	I itap,
2756	I julien,
2757	I gmtime,
2758	I debut,
2759	I lafin,
2760	I nqmax-2,
2761	I nlon,
2762	I nlev,
2763	I dtime,
2764	I u,
2765	I v,
2766	I t,
2767	I paprs,
2768	I pplay,
2769	I pmfu,
2770	I pmfd,
2771	I pen_u,
2772	I pde_u,
2773	I pen_d,
2774	I pde_d,
2775	I ycoefh,
2776	I fm_therm,
2777	I entr_therm,
2778	I yu1,
2779	I yv1,
2780	I ftsol,
2781	I pctsrf,
2782	I rlat,
2783	I frac_impa,
2784	I frac_nucl,
2785	I rlon,
2786	I presnivs,
2787	I pphis,
2788	I pphi,
2789	I albsol,
2790	I qx(1,1,1),
2791	I rhcl,
2792	I cldfra,
2793	I rneb,
2794	I diafra,
2795	I cldliq,
2796	I itop_con,
2797	I ibas_con,
2798	I pmflxr,
2799	I pmflxs,
2800	I prfl,
2801	I psfl,
2802	I da,
2803	I phi,
2804	I mp,
2805	I upwd,
2806	I dnwd,
2807	#ifdef INCA
2808	I flxmass_w,
2809	#endif
2810	O tr_seri)
2811
2812	IF (offline) THEN
2813
2814	print*,'Attention on met a 0 les thermiques pour phystoke'
2815	call phystokenc (
2816	I nlon,nlev,pdtphys,rlon,rlat,
2817	I t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
2818	I fm_therm,entr_therm,
2819	I ycoefh,yu1,yv1,ftsol,pctsrf,
2820	I frac_impa, frac_nucl,
2821	I pphis,airephy,dtime,itap)
2822
2823
2824	ENDIF
2825
2826	c
2827	c Calculer le transport de l'eau et de l'energie (diagnostique)
2828	c
2829	CALL transp (paprs,zxtsol,
2830	e t_seri, q_seri, u_seri, v_seri, zphi,
2831	s ve, vq, ue, uq)
2832	c
2833	cIM diag. bilKP
2834	c
2835	CALL transp_lay (paprs,zxtsol,
2836	e t_seri, q_seri, u_seri, v_seri, zphi,
2837	s ve_lay, vq_lay, ue_lay, uq_lay)
2838	c
2839	c Accumuler les variables a stocker dans les fichiers histoire:
2840	c
2841	c+jld ec_conser
2842	DO k = 1, klev
2843	DO i = 1, klon
2844	ZRCPD = RCPD(1.0+RVTMP2q_seri(i,k))
2845	d_t_ec(i,k)=0.5/ZRCPD
2846	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
2847	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
2848	d_t_ec(i,k) = d_t_ec(i,k)/dtime
2849	END DO
2850	END DO
2851	c-jld ec_conser
2852	IF (if_ebil.ge.1) THEN
2853	ztit='after physic'
2854	CALL diagetpq(airephy,ztit,ip_ebil,1,1,dtime
2855	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2856	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2857	C Comme les tendances de la physique sont ajoute dans la dynamique,
2858	C on devrait avoir que la variation d'entalpie par la dynamique
2859	C est egale a la variation de la physique au pas de temps precedent.
2860	C Donc la somme de ces 2 variations devrait etre nulle.
2861	call diagphy(airephy,ztit,ip_ebil
2862	e , topsw, toplw, solsw, sollw, sens
2863	e , evap, rain_fall, snow_fall, ztsol
2864	e , d_h_vcol, d_qt, d_ec
2865	s , fs_bound, fq_bound )
2866	C
2867	d_h_vcol_phy=d_h_vcol
2868	C
2869	END IF
2870	C
2871	c=======================================================================
2872	c SORTIES
2873	c=======================================================================
2874
2875	cIM Interpolation sur les niveaux de pression du NMC
2876	c -------------------------------------------------
2877	c
2878	#include "calcul_STDlev.h"
2879	c
2880	c slp sea level pressure
2881	slp(:) = paprs(:,1)exp(pphis(:)/(RDt_seri(:,1)))
2882	c
2883	ccc prw = eau precipitable
2884	DO i = 1, klon
2885	prw(i) = 0.
2886	DO k = 1, klev
2887	prw(i) = prw(i) +
2888	. q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
2889	ENDDO
2890	ENDDO
2891	c
2892	cIM initialisation + calculs divers diag AMIP2
2893	c
2894	#include "calcul_divers.h"
2895	c
2896	c=============================================================
2897	c
2898	c Convertir les incrementations en tendances
2899	c
2900	DO k = 1, klev
2901	DO i = 1, klon
2902	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
2903	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
2904	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
2905	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
2906	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
2907	ENDDO
2908	ENDDO
2909	c
2910	IF (nqmax.GE.3) THEN
2911	DO iq = 3, nqmax
2912	DO k = 1, klev
2913	DO i = 1, klon
2914	d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
2915	ENDDO
2916	ENDDO
2917	ENDDO
2918	ENDIF
2919	c
2920	cIM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
2921	c#include "write_bilKP_ins.h"
2922	c#include "write_bilKP_ave.h"
2923	c
2924	c Sauvegarder les valeurs de t et q a la fin de la physique:
2925	c
2926	DO k = 1, klev
2927	DO i = 1, klon
2928	t_ancien(i,k) = t_seri(i,k)
2929	q_ancien(i,k) = q_seri(i,k)
2930	ENDDO
2931	ENDDO
2932	c
2933	c 22.03.04 BEG
2934	c=============================================================
2935	c Ecriture des sorties
2936	c=============================================================
2937	#ifdef CPP_IOIPSL
2938
2939	#ifdef histhf
2940	#include "write_histhf.h"
2941	#endif
2942
2943	#ifdef histday
2944	#include "write_histday.h"
2945	#include "write_histday_seri.h"
2946	#endif
2947
2948	#ifdef histmth
2949	#include "write_histmth.h"
2950	#endif
2951
2952	#ifdef histins
2953	#include "write_histins.h"
2954	#endif
2955
2956	#ifdef histREGDYN
2957	#include "write_histREGDYN.h"
2958	#endif
2959
2960	#ifdef histISCCP
2961	#include "write_histISCCP.h"
2962	#endif
2963
2964
2965	#ifdef histmthNMC
2966	#include "write_histmthNMC.h"
2967	#endif
2968
2969	#endif
2970
2971	#ifdef INCA
2972	#ifdef INCAINFO
2973	WRITE(lunout,*)'Appel CHEMHOOK_END ...'
2974	#endif
2975	CALL chemhook_end (calday,
2976	$ dtime,
2977	$ pplay,
2978	$ t_seri,
2979	$ tr_seri,
2980	$ nbtr,
2981	$ paprs,
2982	#ifdef INCA_CH4
2983	$ q_seri,
2984	#endif
2985	$ annee_ref,
2986	$ day_ini,
2987	#ifdef INCA_AER
2988	$ xjour,
2989	$ pphi,
2990	$ pphis,
2991	$ zx_rh,
2992	$ qx(1,1,1))
2993	#else
2994	$ xjour)
2995	#endif
2996	#ifdef INCAINFO
2997	WRITE(lunout,*)'OK.'
2998	#endif
2999	#endif
3000
3001	c 22.03.04 END
3002	c
3003	c====================================================================
3004	c Si c'est la fin, il faut conserver l'etat de redemarrage
3005	c====================================================================
3006	c
3007	IF (lafin) THEN
3008	itau_phy = itau_phy + itap
3009	ccc IF (ok_oasis) CALL quitcpl
3010	CALL phyredem ("restartphy.nc",dtime,radpas,
3011	. rlat, rlon, pctsrf, ftsol, ftsoil,
3012	cIM "slab" ocean
3013	. tslab, seaice,
3014	. fqsurf, qsol,
3015	. fsnow, falbe,falblw, fevap, rain_fall, snow_fall,
3016	. solsw, sollwdown,dlw,
3017	. radsol,frugs,agesno,
3018	. zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,
3019	. t_ancien, q_ancien, rnebcon, ratqs, clwcon,run_off_lic_0)
3020	ENDIF
3021
3022
3023	RETURN
3024	END
3025	FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
3026	IMPLICIT none
3027	c
3028	c Calculer et imprimer l'eau totale. A utiliser pour verifier
3029	c la conservation de l'eau
3030	c
3031	#include "YOMCST.h"
3032	INTEGER klon,klev
3033	REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
3034	REAL aire(klon)
3035	REAL qtotal, zx, qcheck
3036	INTEGER i, k
3037	c
3038	zx = 0.0
3039	DO i = 1, klon
3040	zx = zx + aire(i)
3041	ENDDO
3042	qtotal = 0.0
3043	DO k = 1, klev
3044	DO i = 1, klon
3045	qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i)
3046	. *(paprs(i,k)-paprs(i,k+1))/RG
3047	ENDDO
3048	ENDDO
3049	c
3050	qcheck = qtotal/zx
3051	c
3052	RETURN
3053	END
3054	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
3055	IMPLICIT none
3056	c
3057	c Tranformer une variable de la grille physique a
3058	c la grille d'ecriture
3059	c
3060	INTEGER nfield,nlon,iim,jjmp1, jjm
3061	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
3062	c
3063	INTEGER i, n, ig
3064	c
3065	jjm = jjmp1 - 1
3066	DO n = 1, nfield
3067	DO i=1,iim
3068	ecrit(i,n) = fi(1,n)
3069	ecrit(i+jjm*iim,n) = fi(nlon,n)
3070	ENDDO
3071	DO ig = 1, nlon - 2
3072	ecrit(iim+ig,n) = fi(1+ig,n)
3073	ENDDO
3074	ENDDO
3075	RETURN
3076	END

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

Download in other formats:

Original Format