Context Navigation

physiq.F @ 559

Last change on this file since 559 was 559, checked in by lmdzadmin, 20 years ago
Initialisations diverses YM LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 90.8 KB

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

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

Download in other formats:

Original Format