Context Navigation

physiq.F @ 677

Last change on this file since 677 was 677, checked in by Laurent Fairhead, 18 years ago
Differentes corrections pour g95 LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 97.5 KB

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

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

Download in other formats:

Original Format