Context Navigation

physiq.F @ 524

Last change on this file since 524 was 524, checked in by lmdzadmin, 20 years ago
Initial revision
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 88.7 KB

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

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

Download in other formats:

Original Format