Context Navigation

physiq.F @ 515

Last change on this file since 515 was 515, checked in by lmdzadmin, 21 years ago
Dernier bug sur histNMC versus histmthNMC IM
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 83.8 KB

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

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

Download in other formats:

Original Format