Context Navigation

physiq.F @ 5415

Last change on this file since 5415 was 705, checked in by Laurent Fairhead, 18 years ago
Pour la conservation du flux d'eau OM, JLD LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 89.1 KB

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

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

Download in other formats:

Original Format