Context Navigation

physiq.F @ 739

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

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

Download in other formats:

Original Format