Context Navigation

physiq.F @ 723

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

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

Download in other formats:

Original Format