Context Navigation

physiq.F @ 625

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

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

Download in other formats:

Original Format