Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

physiq.F

Last change on this file was 556, checked in by lmdzadmin, 20 years ago
Initialisations diverses necessaires au portage Prism AC LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 87.1 KB

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

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

Download in other formats:

Original Format