Context Navigation

physiq.F @ 671

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

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

Download in other formats:

Original Format