Context Navigation

physiq.F @ 2279

Last change on this file since 2279 was 710, checked in by Laurent Fairhead, 18 years ago
Integration des modifications d'OM et JLD pour une meilleure prise en compte de la conservation du flux d'eau LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 89.0 KB

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

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

Download in other formats:

Original Format