Context Navigation

physiq.F @ 594

Last change on this file since 594 was 594, checked in by Laurent Fairhead, 19 years ago
Bugs introduit par la synchro avec LOOP je sais pas comment GG LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 90.9 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 pct_ocean(klon,nbregdyn)
376	REAL rlonPOS(klon)
377
378	c sorties ISCCP
379
380	logical ok_isccp
381	real ecrit_isccp
382	integer nid_isccp
383	save ok_isccp, ecrit_isccp, nid_isccp
384
385	#ifdef histISCCP
386	data ok_isccp/.true./
387	#else
388	data ok_isccp/.false./
389	#endif
390
391	c sorties statistiques regime dynamique
392	logical ok_regdyn
393	real ecrit_regdyn
394	integer nid_regdyn
395	save ok_regdyn, ecrit_regdyn, nid_regdyn
396
397	#ifdef histREGDYN
398	c data ok_regdyn,ecrit_regdyn/.true.,0.125/
399	c data ok_regdyn,ecrit_regdyn/.true.,1./
400	data ok_regdyn/.true./
401	#else
402	data ok_regdyn/.false./
403	#endif
404
405	REAL zx_tau(kmaxm1), zx_pc(lmaxm1), zx_o500(iwmax)
406	DATA zx_tau/0.0, 0.3, 1.3, 3.6, 9.4, 23., 60./
407	DATA zx_pc/50., 180., 310., 440., 560., 680., 800./
408
409	c cldtopres pression au sommet des nuages
410	REAL cldtopres(lmaxm1)
411	DATA cldtopres/50., 180., 310., 440., 560., 680., 800./
412
413	INTEGER komega, nhoriRD
414
415	c taulev: numero du niveau de tau dans les sorties ISCCP
416	CHARACTER *4 taulev(kmaxm1)
417	DATA taulev/'tau1','tau2','tau3','tau4','tau5','tau6','tau7'/
418
419	REAL zx_lonx7(iimx7), zx_latx7(jjmp1x7)
420	INTEGER nhorix7
421	cIM: region='3d' <==> sorties en global
422	CHARACTER*3 region
423	PARAMETER(region='3d')
424	c
425	logical ok_hf
426	real ecrit_hf
427	integer nid_hf, nid_hf3d
428	save ok_hf, ecrit_hf, nid_hf, nid_hf3d
429
430	c QUESTION : noms de variables ?
431
432	#ifdef histhf
433	data ok_hf,ecrit_hf/.true.,0.25/
434	#else
435	data ok_hf/.false./
436	#endif
437
438	INTEGER longcles
439	PARAMETER ( longcles = 20 )
440	REAL clesphy0( longcles )
441	c
442	c Variables quasi-arguments
443	c
444	REAL xjour
445	SAVE xjour
446	c
447	c
448	c Variables propres a la physique
449	c
450	REAL dtime
451	SAVE dtime ! pas temporel de la physique
452	c
453	INTEGER radpas
454	SAVE radpas ! frequence d'appel rayonnement
455	c
456	REAL radsol(klon)
457	SAVE radsol ! bilan radiatif au sol calcule par code radiatif
458	c
459	REAL rlat(klon)
460	SAVE rlat ! latitude pour chaque point
461	c
462	REAL rlon(klon)
463	SAVE rlon ! longitude pour chaque point
464	c
465	cc INTEGER iflag_con
466	cc SAVE iflag_con ! indicateur de la convection
467	c
468	INTEGER itap
469	SAVE itap ! compteur pour la physique
470	c
471	REAL co2_ppm_etat0
472	c
473	REAL solaire_etat0
474	c
475	real slp(klon) ! sea level pressure
476
477	REAL ftsol(klon,nbsrf)
478	SAVE ftsol ! temperature du sol
479	c
480	REAL ftsoil(klon,nsoilmx,nbsrf)
481	SAVE ftsoil ! temperature dans le sol
482	c
483	REAL fevap(klon,nbsrf)
484	SAVE fevap ! evaporation
485	REAL fluxlat(klon,nbsrf)
486	SAVE fluxlat
487	c
488	REAL deltat(klon)
489	SAVE deltat ! ecart avec la SST de reference
490	c
491	REAL fqsurf(klon,nbsrf)
492	SAVE fqsurf ! humidite de l'air au contact de la surface
493	c
494	REAL qsol(klon)
495	SAVE qsol ! hauteur d'eau dans le sol
496	c
497	REAL fsnow(klon,nbsrf)
498	SAVE fsnow ! epaisseur neigeuse
499	c
500	REAL falbe(klon,nbsrf)
501	SAVE falbe ! albedo par type de surface
502	REAL falblw(klon,nbsrf)
503	SAVE falblw ! albedo par type de surface
504
505	c
506	c
507	c Parametres de l'Orographie a l'Echelle Sous-Maille (OESM):
508	c
509	REAL zmea(klon)
510	SAVE zmea ! orographie moyenne
511	c
512	REAL zstd(klon)
513	SAVE zstd ! deviation standard de l'OESM
514	c
515	REAL zsig(klon)
516	SAVE zsig ! pente de l'OESM
517	c
518	REAL zgam(klon)
519	save zgam ! anisotropie de l'OESM
520	c
521	REAL zthe(klon)
522	SAVE zthe ! orientation de l'OESM
523	c
524	REAL zpic(klon)
525	SAVE zpic ! Maximum de l'OESM
526	c
527	REAL zval(klon)
528	SAVE zval ! Minimum de l'OESM
529	c
530	REAL rugoro(klon)
531	SAVE rugoro ! longueur de rugosite de l'OESM
532	c
533	REAL zulow(klon),zvlow(klon),zustr(klon), zvstr(klon)
534	c
535	REAL zuthe(klon),zvthe(klon)
536	SAVE zuthe
537	SAVE zvthe
538	INTEGER igwd,idx(klon),itest(klon)
539	c
540	REAL agesno(klon,nbsrf)
541	SAVE agesno ! age de la neige
542	c
543	REAL alb_neig(klon)
544	SAVE alb_neig ! albedo de la neige
545	c
546	REAL run_off_lic_0(klon)
547	SAVE run_off_lic_0
548	cKE43
549	c Variables liees a la convection de K. Emanuel (sb):
550	c
551	REAL ema_workcbmf(klon) ! cloud base mass flux
552	SAVE ema_workcbmf
553
554	REAL ema_cbmf(klon) ! cloud base mass flux
555	SAVE ema_cbmf
556
557	REAL ema_pcb(klon) ! cloud base pressure
558	SAVE ema_pcb
559
560	REAL ema_pct(klon) ! cloud top pressure
561	SAVE ema_pct
562
563	REAL bas, top ! cloud base and top levels
564	SAVE bas
565	SAVE top
566
567	REAL Ma(klon,klev) ! undilute upward mass flux
568	SAVE Ma
569	REAL qcondc(klon,klev) ! in-cld water content from convect
570	SAVE qcondc
571	REAL ema_work1(klon, klev), ema_work2(klon, klev)
572	SAVE ema_work1, ema_work2
573	REAL wdn(klon), tdn(klon), qdn(klon)
574
575	REAL wd(klon) ! sb
576	SAVE wd ! sb
577
578	c Variables locales pour la couche limite (al1):
579	c
580	cAl1 REAL pblh(klon) ! Hauteur de couche limite
581	cAl1 SAVE pblh
582	c34EK
583	c
584	c Variables locales:
585	c
586	REAL cdragh(klon) ! drag coefficient pour T and Q
587	REAL cdragm(klon) ! drag coefficient pour vent
588	cAA
589	cAA Pour phytrac
590	cAA
591	REAL ycoefh(klon,klev) ! coef d'echange pour phytrac
592	REAL yu1(klon) ! vents dans la premiere couche U
593	REAL yv1(klon) ! vents dans la premiere couche V
594	REAL ffonte(klon,nbsrf) !Flux thermique utilise pour fondre la neige
595	REAL fqcalving(klon,nbsrf) !Flux d'eau "perdue" par la surface
596	c !et necessaire pour limiter la
597	c !hauteur de neige, en kg/m2/s
598	REAL zxffonte(klon), zxfqcalving(klon)
599
600	c$$$ LOGICAL offline ! Controle du stockage ds "physique"
601	c$$$ PARAMETER (offline=.false.)
602	c$$$ INTEGER physid
603	REAL pfrac_impa(klon,klev)! Produits des coefs lessivage impaction
604	save pfrac_impa
605	REAL pfrac_nucl(klon,klev)! Produits des coefs lessivage nucleation
606	save pfrac_nucl
607	REAL pfrac_1nucl(klon,klev)! Produits des coefs lessi nucl (alpha = 1)
608	save pfrac_1nucl
609	REAL frac_impa(klon,klev) ! fractions d'aerosols lessivees (impaction)
610	REAL frac_nucl(klon,klev) ! idem (nucleation)
611	#ifdef INCA
612	INTEGER :: iii
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	topswai(:)=0.
1008	topswad(:)=0.
1009	solswai(:)=0.
1010	solswad(:)=0.
1011	!rv
1012
1013	cym
1014	wfbils(:,:)=0
1015	cym
1016	modname = 'physiq'
1017	IF (if_ebil.ge.1) THEN
1018	DO i=1,klon
1019	zero_v(i)=0.
1020	END DO
1021	END IF
1022	ok_sync=.TRUE.
1023	IF (nqmax .LT. 2) THEN
1024	abort_message = 'eaux vapeur et liquide sont indispensables'
1025	CALL abort_gcm (modname,abort_message,1)
1026	ENDIF
1027	IF (debut) THEN
1028	CALL suphec ! initialiser constantes et parametres phys.
1029	c
1030	cIM 050204 BEG
1031	DO i=1, klon
1032	nday_rain(i)=0.
1033	ENDDO
1034	cIM 050204 END
1035	c
1036	c======================================================================
1037	cIM BEG
1038	DO k=1, nlevENS
1039	DO l=1, nlevSTD
1040	c
1041	bb=clevSTD(l)
1042	c
1043	IF(l.GE.2) THEN
1044	aa=clevSTD(l)
1045	bb=aa(1:lnblnk1(aa))
1046	ENDIF
1047	c
1048	IF(bb.EQ.clev(k)) THEN
1049	c print*,'k=',k,'l=',l,'clev=',clev(k)
1050	indENS(k)=l
1051	c print*,'k=',k,'l=',l,'clev=',clev(k),'indENS=',indENS(k)
1052	ENDIF
1053	c
1054	ENDDO
1055	ENDDO
1056	c
1057	ENDIF !debut
1058	cIM END
1059	xjour = rjourvrai
1060	c
1061	c Si c'est le debut, il faut initialiser plusieurs choses
1062	c ********
1063	c
1064	IF (debut) THEN
1065	C
1066	IF (if_ebil.ge.1) d_h_vcol_phy=0.
1067	c
1068	c appel a la lecture du run.def physique
1069	c
1070	call conf_phys(ocean, ok_veget, ok_journe, ok_mensuel,
1071	. ok_instan, fact_cldcon, facttemps,ok_newmicro,
1072	. iflag_cldcon,ratqsbas,ratqshaut, if_ebil,
1073	. ok_ade, ok_aie,
1074	. bl95_b0, bl95_b1,
1075	. iflag_thermals,nsplit_thermals)
1076
1077	c
1078	c
1079	c Initialiser les compteurs:
1080	c
1081
1082	frugs = 0.
1083	itap = 0
1084	itaprad = 0
1085	CALL phyetat0 ("startphy.nc",dtime,co2_ppm_etat0,solaire_etat0,
1086	. rlat,rlon,pctsrf, ftsol,ftsoil,deltat,fqsurf,qsol,fsnow,
1087	. falbe, falblw, fevap, rain_fall,snow_fall,solsw, sollwdown,
1088	. dlw,radsol,frugs,agesno,clesphy0,
1089	. zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,tabcntr0,
1090	. t_ancien, q_ancien, ancien_ok, rnebcon, ratqs,clwcon,
1091	. run_off_lic_0)
1092
1093	c ATTENTION : il faudra a terme relire q2 dans l'etat initial
1094	q2(:,:,:)=1.e-8
1095	c
1096	radpas = NINT( 86400./dtime/nbapp_rad)
1097	c
1098	C on remet le calendrier a zero
1099	c
1100	IF (raz_date .eq. 1) THEN
1101	itau_phy = 0
1102	ENDIF
1103
1104	c
1105	CALL printflag( tabcntr0,radpas,ok_ocean,ok_oasis ,ok_journe,
1106	, ok_instan, ok_region )
1107	c
1108	IF (ABS(dtime-pdtphys).GT.0.001) THEN
1109	WRITE(lunout,*) 'Pas physique n est pas correct',dtime,
1110	. pdtphys
1111	abort_message='Pas physique n est pas correct '
1112	call abort_gcm(modname,abort_message,1)
1113	ENDIF
1114	IF (nlon .NE. klon) THEN
1115	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,
1116	. klon
1117	abort_message='nlon et klon ne sont pas coherents'
1118	call abort_gcm(modname,abort_message,1)
1119	ENDIF
1120	IF (nlev .NE. klev) THEN
1121	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,
1122	. klev
1123	abort_message='nlev et klev ne sont pas coherents'
1124	call abort_gcm(modname,abort_message,1)
1125	ENDIF
1126	c
1127	IF (dtime*FLOAT(radpas).GT.21600..AND.cycle_diurne) THEN
1128	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
1129	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
1130	abort_message='Nbre d appels au rayonnement insuffisant'
1131	call abort_gcm(modname,abort_message,1)
1132	ENDIF
1133	WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
1134	WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=",
1135	. ok_cvl
1136	c
1137	cKE43
1138	c Initialisation pour la convection de K.E. (sb):
1139	IF (iflag_con.GE.3) THEN
1140
1141	WRITE(lunout,)"** Convection de Kerry Emanuel 4.3 "
1142	WRITE(lunout,*)
1143	. "On va utiliser le melange convectif des traceurs qui"
1144	WRITE(lunout,*)"est calcule dans convect4.3"
1145	WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
1146
1147	DO i = 1, klon
1148	ema_cbmf(i) = 0.
1149	ema_pcb(i) = 0.
1150	ema_pct(i) = 0.
1151	ema_workcbmf(i) = 0.
1152	ENDDO
1153	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
1154	DO i = 1, klon
1155	ibas_con(i) = 1
1156	itop_con(i) = klev+1
1157	ENDDO
1158	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
1159
1160	ENDIF
1161
1162	c34EK
1163	IF (ok_orodr) THEN
1164	DO i=1,klon
1165	rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1166	ENDDO
1167	CALL SUGWD(klon,klev,paprs,pplay)
1168	DO i=1,klon
1169	zuthe(i)=0.
1170	zvthe(i)=0.
1171	if(zstd(i).gt.10.)then
1172	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
1173	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
1174	endif
1175	ENDDO
1176	ENDIF
1177	c
1178	c
1179	lmt_pas = NINT(86400./dtime * 1.0) ! tous les jours
1180	WRITE(lunout,*)'La frequence de lecture surface est de ',
1181	. lmt_pas
1182	c
1183	ecrit_mth = NINT(86400./dtime *ecritphy) ! tous les ecritphy jours
1184	IF (ok_mensuel) THEN
1185	WRITE(lunout,*)'La frequence de sortie mensuelle est de ',
1186	. ecrit_mth
1187	ENDIF
1188	ecrit_day = NINT(86400./dtime *1.0) ! tous les jours
1189	IF (ok_journe) THEN
1190	WRITE(lunout,*)'La frequence de sortie journaliere est de ',
1191	. ecrit_day
1192	ENDIF
1193	ccc ecrit_ins = NINT(86400./dtime *0.5) ! 2 fois par jour
1194	ccc ecrit_ins = NINT(86400./dtime *0.25) ! 4 fois par jour
1195	ecrit_ins = NINT(86400./dtime/48.) ! a chaque pas de temps ==> PB. dans time_counter pour 1mois
1196	ecrit_ins = NINT(86400./dtime/12.) ! toutes les deux heures
1197	IF (ok_instan) THEN
1198	WRITE(lunout,*)'La frequence de sortie instant. est de ',
1199	. ecrit_ins
1200	ENDIF
1201	ecrit_reg = NINT(86400./dtime *0.25) ! 4 fois par jour
1202	IF (ok_region) THEN
1203	WRITE(lunout,*)'La frequence de sortie region est de ',
1204	. ecrit_reg
1205	ENDIF
1206
1207	c
1208	c Initialiser le couplage si necessaire
1209	c
1210	npas = 0
1211	nexca = 0
1212	if (ocean == 'couple') then
1213	npas = itaufin/ iphysiq
1214	nexca = 86400 / dtime
1215	write(lunout,*)' ##### Ocean couple #####'
1216	write(lunout,*)' Valeurs des pas de temps'
1217	write(lunout,*)' npas = ', npas
1218	write(lunout,*)' nexca = ', nexca
1219	endif
1220	c
1221	c
1222	cIM
1223	capemaxcels = 't_max(X)'
1224	t2mincels = 't_min(X)'
1225	t2maxcels = 't_max(X)'
1226
1227	c
1228	c=============================================================
1229	c Initialisation des sorties
1230	c=============================================================
1231
1232	#ifdef CPP_IOIPSL
1233
1234	#ifdef histhf
1235	#include "ini_histhf.h"
1236	#endif
1237
1238	#ifdef histday
1239	#include "ini_histday.h"
1240	#endif
1241
1242	#ifdef histmth
1243	#include "ini_histmth.h"
1244	#endif
1245
1246	#ifdef histins
1247	#include "ini_histins.h"
1248	#endif
1249
1250	#ifdef histISCCP
1251	#include "ini_histISCCP.h"
1252	#endif
1253
1254	#ifdef histmthNMC
1255	#include "ini_histmthNMC.h"
1256	#endif
1257
1258	#ifdef histREGDYN
1259	#include "ini_histREGDYN.h"
1260	#endif
1261
1262	#ifdef histISCCP
1263	#include "ini_histISCCP.h"
1264	#endif
1265	#endif
1266
1267	cXXXPB Positionner date0 pour initialisation de ORCHIDEE
1268	date0 = zjulian
1269	C date0 = day_ini
1270	WRITE(,) 'physiq date0 : ',date0
1271	c
1272	c
1273	c
1274	c Prescrire l'ozone dans l'atmosphere
1275	c
1276	c
1277	cc DO i = 1, klon
1278	cc DO k = 1, klev
1279	cc CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1280	cc ENDDO
1281	cc ENDDO
1282	c
1283	#ifdef INCA
1284	iii = MOD(NINT(xjour),360)
1285	calday = FLOAT(iii) + gmtime
1286	WRITE(lunout,*) 'initial time ', xjour, calday
1287	#ifdef INCAINFO
1288	WRITE(lunout,*) 'Appel CHEMINI ...'
1289	#endif
1290	CALL chemini( rpi,
1291	$ rg,
1292	$ ra,
1293	$ airephy,
1294	$ rlat,
1295	$ rlon,
1296	$ presnivs,
1297	$ calday,
1298	$ tracnam,
1299	$ natsnam,
1300	c $ mxoutflds,
1301	c $ outinst,
1302	c $ outtimav,
1303	$ klon,
1304	$ nqmax,
1305	$ pdtphys,
1306	$ annee_ref,
1307	$ day_ini)
1308	#ifdef INCAINFO
1309	WRITE(lunout,*) 'OK.'
1310	#endif
1311	#endif
1312	c
1313	ENDIF
1314	c
1315	c ************** Fin de IF ( debut ) *************
1316	c
1317	c
1318	c Mettre a zero des variables de sortie (pour securite)
1319	c
1320	DO i = 1, klon
1321	d_ps(i) = 0.0
1322	ENDDO
1323	DO k = 1, klev
1324	DO i = 1, klon
1325	d_t(i,k) = 0.0
1326	d_u(i,k) = 0.0
1327	d_v(i,k) = 0.0
1328	ENDDO
1329	ENDDO
1330	DO iq = 1, nqmax
1331	DO k = 1, klev
1332	DO i = 1, klon
1333	d_qx(i,k,iq) = 0.0
1334	ENDDO
1335	ENDDO
1336	ENDDO
1337	c
1338	c Ne pas affecter les valeurs entrees de u, v, h, et q
1339	c
1340	DO k = 1, klev
1341	DO i = 1, klon
1342	t_seri(i,k) = t(i,k)
1343	u_seri(i,k) = u(i,k)
1344	v_seri(i,k) = v(i,k)
1345	q_seri(i,k) = qx(i,k,ivap)
1346	ql_seri(i,k) = qx(i,k,iliq)
1347	qs_seri(i,k) = 0.
1348	ENDDO
1349	ENDDO
1350	IF (nqmax.GE.3) THEN
1351	DO iq = 3, nqmax
1352	DO k = 1, klev
1353	DO i = 1, klon
1354	tr_seri(i,k,iq-2) = qx(i,k,iq)
1355	ENDDO
1356	ENDDO
1357	ENDDO
1358	ELSE
1359	DO k = 1, klev
1360	DO i = 1, klon
1361	tr_seri(i,k,1) = 0.0
1362	ENDDO
1363	ENDDO
1364	ENDIF
1365	C
1366	DO i = 1, klon
1367	ztsol(i) = 0.
1368	ENDDO
1369	DO nsrf = 1, nbsrf
1370	DO i = 1, klon
1371	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1372	ENDDO
1373	ENDDO
1374	C
1375	IF (if_ebil.ge.1) THEN
1376	ztit='after dynamic'
1377	CALL diagetpq(airephy,ztit,ip_ebil,1,1,dtime
1378	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1379	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1380	C Comme les tendances de la physique sont ajoute dans la dynamique,
1381	C on devrait avoir que la variation d'entalpie par la dynamique
1382	C est egale a la variation de la physique au pas de temps precedent.
1383	C Donc la somme de ces 2 variations devrait etre nulle.
1384	call diagphy(airephy,ztit,ip_ebil
1385	e , zero_v, zero_v, zero_v, zero_v, zero_v
1386	e , zero_v, zero_v, zero_v, ztsol
1387	e , d_h_vcol+d_h_vcol_phy, d_qt, 0.
1388	s , fs_bound, fq_bound )
1389	END IF
1390
1391	c Diagnostiquer la tendance dynamique
1392	c
1393	IF (ancien_ok) THEN
1394	DO k = 1, klev
1395	DO i = 1, klon
1396	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
1397	d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
1398	ENDDO
1399	ENDDO
1400	ELSE
1401	DO k = 1, klev
1402	DO i = 1, klon
1403	d_t_dyn(i,k) = 0.0
1404	d_q_dyn(i,k) = 0.0
1405	ENDDO
1406	ENDDO
1407	ancien_ok = .TRUE.
1408	ENDIF
1409	c
1410	c Ajouter le geopotentiel du sol:
1411	c
1412	DO k = 1, klev
1413	DO i = 1, klon
1414	zphi(i,k) = pphi(i,k) + pphis(i)
1415	ENDDO
1416	ENDDO
1417	c
1418	c Verifier les temperatures
1419	c
1420	CALL hgardfou(t_seri,ftsol,'debutphy')
1421	c
1422	c Incrementer le compteur de la physique
1423	c
1424	itap = itap + 1
1425	julien = MOD(NINT(xjour),360)
1426	if (julien .eq. 0) julien = 360
1427	c
1428	c Mettre en action les conditions aux limites (albedo, sst, etc.).
1429	c Prescrire l'ozone et calculer l'albedo sur l'ocean.
1430	c
1431	IF (MOD(itap-1,lmt_pas) .EQ. 0) THEN
1432	WRITE(lunout,*)' PHYS cond julien ',julien
1433	CALL ozonecm( FLOAT(julien), rlat, paprs, wo)
1434	ENDIF
1435	c
1436	c Re-evaporer l'eau liquide nuageuse
1437	c
1438	DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse
1439	DO i = 1, klon
1440	zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1441	c zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1442	zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1443	zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
1444	zb = MAX(0.0,ql_seri(i,k))
1445	za = - MAX(0.0,ql_seri(i,k))
1446	. * (zlvdcp(1.-zdelta)+zlsdcpzdelta)
1447	t_seri(i,k) = t_seri(i,k) + za
1448	q_seri(i,k) = q_seri(i,k) + zb
1449	ql_seri(i,k) = 0.0
1450	d_t_eva(i,k) = za
1451	d_q_eva(i,k) = zb
1452	ENDDO
1453	ENDDO
1454	c
1455	IF (if_ebil.ge.2) THEN
1456	ztit='after reevap'
1457	CALL diagetpq(airephy,ztit,ip_ebil,2,1,dtime
1458	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1459	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1460	call diagphy(airephy,ztit,ip_ebil
1461	e , zero_v, zero_v, zero_v, zero_v, zero_v
1462	e , zero_v, zero_v, zero_v, ztsol
1463	e , d_h_vcol, d_qt, d_ec
1464	s , fs_bound, fq_bound )
1465	C
1466	END IF
1467	C
1468	c
1469	c Appeler la diffusion verticale (programme de couche limite)
1470	c
1471	DO i = 1, klon
1472	c if (.not. ok_veget) then
1473	c frugs(i,is_ter) = SQRT(frugs(i,is_ter)2+rugoro(i)2)
1474	c endif
1475	c frugs(i,is_lic) = rugoro(i)
1476	c frugs(i,is_oce) = rugmer(i)
1477	c frugs(i,is_sic) = 0.001
1478	zxrugs(i) = 0.0
1479	ENDDO
1480	DO nsrf = 1, nbsrf
1481	DO i = 1, klon
1482	c frugs(i,nsrf) = MAX(frugs(i,nsrf),0.001)
1483	frugs(i,nsrf) = MAX(frugs(i,nsrf),0.000015)
1484	ENDDO
1485	ENDDO
1486	DO nsrf = 1, nbsrf
1487	DO i = 1, klon
1488	zxrugs(i) = zxrugs(i) + frugs(i,nsrf)*pctsrf(i,nsrf)
1489	ENDDO
1490	ENDDO
1491	c
1492	C calculs necessaires au calcul de l'albedo dans l'interface
1493	c
1494	CALL orbite(FLOAT(julien),zlongi,dist)
1495	IF (cycle_diurne) THEN
1496	zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
1497	CALL zenang(zlongi,gmtime,zdtime,rlat,rlon,rmu0,fract)
1498	ELSE
1499	rmu0 = -999.999
1500	ENDIF
1501	cIM BEG
1502	DO i=1, klon
1503	sunlit(i)=1
1504	IF(rmu0(i).EQ.0.) sunlit(i)=0
1505	nbsunlit(1,i)=FLOAT(sunlit(i))
1506	ENDDO
1507	cIM END
1508	C Calcul de l'abedo moyen par maille
1509	albsol(:)=0.
1510	albsollw(:)=0.
1511	DO nsrf = 1, nbsrf
1512	DO i = 1, klon
1513	albsol(i) = albsol(i) + falbe(i,nsrf) * pctsrf(i,nsrf)
1514	albsollw(i) = albsollw(i) + falblw(i,nsrf) * pctsrf(i,nsrf)
1515	ENDDO
1516	ENDDO
1517	C
1518	C Repartition sous maille des flux LW et SW
1519	C Modif OM+PASB+JLD
1520	C Repartition du longwave par sous-surface linearisee
1521	Cn
1522
1523	DO nsrf = 1, nbsrf
1524	DO i = 1, klon
1525	c$$$ fsollw(i,nsrf) = sollwdown(i) - RSIGMAftsol(i,nsrf)*4
1526	c$$$ fsollw(i,nsrf) = sollw(i)
1527	fsollw(i,nsrf) = sollw(i)
1528	$ + 4.0RSIGMAztsol(i)*3 (ztsol(i)-ftsol(i,nsrf))
1529	fsolsw(i,nsrf) = solsw(i)*(1.-falbe(i,nsrf))/(1.-albsol(i))
1530	ENDDO
1531	ENDDO
1532
1533	fder = dlw
1534
1535
1536	CALL clmain(dtime,itap,date0,pctsrf,pctsrf_new,
1537	e t_seri,q_seri,u_seri,v_seri,
1538	e julien, rmu0, co2_ppm,
1539	e ok_veget, ocean, npas, nexca, ftsol,
1540	$ soil_model,cdmmax, cdhmax,
1541	$ ksta, ksta_ter, ok_kzmin, ftsoil, qsol,
1542	$ paprs,pplay,radsol, fsnow,fqsurf,fevap,falbe,falblw,
1543	$ fluxlat,
1544	cIM cf. JLD e rain_fall, snow_fall, solsw, sollw, sollwdown, fder,
1545	e rain_fall, snow_fall, fsolsw, fsollw, sollwdown, fder,
1546	e rlon, rlat, cuphy, cvphy, frugs,
1547	e debut, lafin, agesno,rugoro ,
1548	s d_t_vdf,d_q_vdf,d_u_vdf,d_v_vdf,d_ts,
1549	s fluxt,fluxq,fluxu,fluxv,cdragh,cdragm,
1550	s q2,
1551	s dsens, devap,
1552	s ycoefh,yu1,yv1, t2m, q2m, u10m, v10m,
1553	s fqcalving, ffonte, run_off_lic_0)
1554	c
1555	CXXX PB
1556	CXXX Incrementation des flux
1557	CXXX
1558
1559	zxfluxt=0.
1560	zxfluxq=0.
1561	zxfluxu=0.
1562	zxfluxv=0.
1563	DO nsrf = 1, nbsrf
1564	DO k = 1, klev
1565	DO i = 1, klon
1566	zxfluxt(i,k) = zxfluxt(i,k) +
1567	$ fluxt(i,k,nsrf) * pctsrf( i, nsrf)
1568	zxfluxq(i,k) = zxfluxq(i,k) +
1569	$ fluxq(i,k,nsrf) * pctsrf( i, nsrf)
1570	zxfluxu(i,k) = zxfluxu(i,k) +
1571	$ fluxu(i,k,nsrf) * pctsrf( i, nsrf)
1572	zxfluxv(i,k) = zxfluxv(i,k) +
1573	$ fluxv(i,k,nsrf) * pctsrf( i, nsrf)
1574	END DO
1575	END DO
1576	END DO
1577	DO i = 1, klon
1578	sens(i) = - zxfluxt(i,1) ! flux de chaleur sensible au sol
1579	c evap(i) = - fluxq(i,1) ! flux d'evaporation au sol
1580	evap(i) = - zxfluxq(i,1) ! flux d'evaporation au sol
1581	fder(i) = dlw(i) + dsens(i) + devap(i)
1582	ENDDO
1583
1584
1585	DO k = 1, klev
1586	DO i = 1, klon
1587	t_seri(i,k) = t_seri(i,k) + d_t_vdf(i,k)
1588	q_seri(i,k) = q_seri(i,k) + d_q_vdf(i,k)
1589	u_seri(i,k) = u_seri(i,k) + d_u_vdf(i,k)
1590	v_seri(i,k) = v_seri(i,k) + d_v_vdf(i,k)
1591	ENDDO
1592	ENDDO
1593	c
1594	IF (if_ebil.ge.2) THEN
1595	ztit='after clmain'
1596	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
1597	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1598	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1599	call diagphy(airephy,ztit,ip_ebil
1600	e , zero_v, zero_v, zero_v, zero_v, sens
1601	e , evap , zero_v, zero_v, ztsol
1602	e , d_h_vcol, d_qt, d_ec
1603	s , fs_bound, fq_bound )
1604	END IF
1605	C
1606	c
1607	c Incrementer la temperature du sol
1608	c
1609	DO i = 1, klon
1610	zxtsol(i) = 0.0
1611	zxfluxlat(i) = 0.0
1612	c
1613	zt2m(i) = 0.0
1614	zq2m(i) = 0.0
1615	zu10m(i) = 0.0
1616	zv10m(i) = 0.0
1617	cIM cf JLD ??
1618	zxffonte(i) = 0.0
1619	zxfqcalving(i) = 0.0
1620	c
1621	IF ( abs( pctsrf(i, is_ter) + pctsrf(i, is_lic) +
1622	$ pctsrf(i, is_oce) + pctsrf(i, is_sic) - 1.) .GT. EPSFRA)
1623	$ THEN
1624	WRITE(,) 'physiq : pb sous surface au point ', i,
1625	$ pctsrf(i, 1 : nbsrf)
1626	ENDIF
1627	ENDDO
1628	DO nsrf = 1, nbsrf
1629	DO i = 1, klon
1630	c IF (pctsrf(i,nsrf) .GE. EPSFRA) THEN
1631	ftsol(i,nsrf) = ftsol(i,nsrf) + d_ts(i,nsrf)
1632	cIM cf. JLD
1633	wfbils(i,nsrf) = ( fsolsw(i,nsrf) + fsollw(i,nsrf)
1634	$ + fluxt(i,1,nsrf) + fluxlat(i,nsrf) ) * pctsrf(i,nsrf)
1635	zxtsol(i) = zxtsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1636	zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf)*pctsrf(i,nsrf)
1637	cccIM
1638	zt2m(i) = zt2m(i) + t2m(i,nsrf)*pctsrf(i,nsrf)
1639	zq2m(i) = zq2m(i) + q2m(i,nsrf)*pctsrf(i,nsrf)
1640	zu10m(i) = zu10m(i) + u10m(i,nsrf)*pctsrf(i,nsrf)
1641	zv10m(i) = zv10m(i) + v10m(i,nsrf)*pctsrf(i,nsrf)
1642	cIM cf JLD ??
1643	zxffonte(i) = zxffonte(i) + ffonte(i,nsrf)*pctsrf(i,nsrf)
1644	zxfqcalving(i) = zxfqcalving(i) +
1645	. fqcalving(i,nsrf)*pctsrf(i,nsrf)
1646	c ENDIF
1647	ENDDO
1648	ENDDO
1649
1650	c
1651	c Si une sous-fraction n'existe pas, elle prend la temp. moyenne
1652	c
1653	DO nsrf = 1, nbsrf
1654	DO i = 1, klon
1655	IF (pctsrf(i,nsrf) .LT. epsfra) ftsol(i,nsrf) = zxtsol(i)
1656	cccIM
1657	IF (pctsrf(i,nsrf) .LT. epsfra) t2m(i,nsrf) = zt2m(i)
1658	IF (pctsrf(i,nsrf) .LT. epsfra) q2m(i,nsrf) = zq2m(i)
1659	IF (pctsrf(i,nsrf) .LT. epsfra) u10m(i,nsrf) = zu10m(i)
1660	IF (pctsrf(i,nsrf) .LT. epsfra) v10m(i,nsrf) = zv10m(i)
1661	cIM cf JLD ??
1662	IF (pctsrf(i,nsrf) .LT. epsfra) ffonte(i,nsrf) = zxffonte(i)
1663	IF (pctsrf(i,nsrf) .LT. epsfra)
1664	. fqcalving(i,nsrf) = zxfqcalving(i)
1665	ENDDO
1666	ENDDO
1667	c
1668	c
1669	c Calculer la derive du flux infrarouge
1670	c
1671	cXXX DO nsrf = 1, nbsrf
1672	DO i = 1, klon
1673	cXXX IF (pctsrf(i,nsrf) .GE. EPSFRA) THEN
1674	dlw(i) = - 4.0RSIGMAzxtsol(i)**3
1675	cXXX . *(ftsol(i,nsrf)-zxtsol(i))
1676	cXXX . *pctsrf(i,nsrf)
1677	cXXX ENDIF
1678	cXXX ENDDO
1679	ENDDO
1680	c
1681	c Appeler la convection (au choix)
1682	c
1683	DO k = 1, klev
1684	DO i = 1, klon
1685	conv_q(i,k) = d_q_dyn(i,k)
1686	. + d_q_vdf(i,k)/dtime
1687	conv_t(i,k) = d_t_dyn(i,k)
1688	. + d_t_vdf(i,k)/dtime
1689	ENDDO
1690	ENDDO
1691	IF (check) THEN
1692	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
1693	WRITE(lunout,*) "avantcon=", za
1694	ENDIF
1695	zx_ajustq = .FALSE.
1696	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
1697	IF (zx_ajustq) THEN
1698	DO i = 1, klon
1699	z_avant(i) = 0.0
1700	ENDDO
1701	DO k = 1, klev
1702	DO i = 1, klon
1703	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k))
1704	. *(paprs(i,k)-paprs(i,k+1))/RG
1705	ENDDO
1706	ENDDO
1707	ENDIF
1708	IF (iflag_con.EQ.1) THEN
1709	stop'reactiver le call conlmd dans physiq.F'
1710	c CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
1711	c . d_t_con, d_q_con,
1712	c . rain_con, snow_con, ibas_con, itop_con)
1713	ELSE IF (iflag_con.EQ.2) THEN
1714	CALL conflx(dtime, paprs, pplay, t_seri, q_seri,
1715	e conv_t, conv_q, zxfluxq(1,1), omega,
1716	s d_t_con, d_q_con, rain_con, snow_con,
1717	s pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
1718	s kcbot, kctop, kdtop, pmflxr, pmflxs)
1719	WHERE (rain_con < 0.) rain_con = 0.
1720	WHERE (snow_con < 0.) snow_con = 0.
1721	DO i = 1, klon
1722	ibas_con(i) = klev+1 - kcbot(i)
1723	itop_con(i) = klev+1 - kctop(i)
1724	ENDDO
1725	ELSE IF (iflag_con.GE.3) THEN
1726	c nb of tracers for the KE convection:
1727	if (nqmax .GE. 4) then
1728	ntra = nbtr
1729	else
1730	ntra = 1
1731	endif
1732	c
1733	c sb, oct02:
1734	c Schema de convection modularise et vectorise:
1735	c (driver commun aux versions 3 et 4)
1736	c
1737	IF (ok_cvl) THEN ! new driver for convectL
1738
1739	CALL concvl (iflag_con,
1740	. dtime,paprs,pplay,t_seri,q_seri,
1741	. u_seri,v_seri,tr_seri,nbtr,
1742	. ema_work1,ema_work2,
1743	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
1744	. rain_con, snow_con, ibas_con, itop_con,
1745	. upwd,dnwd,dnwd0,
1746	. Ma,cape,tvp,iflagctrl,
1747	. pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd)
1748	cIM cf. FH
1749	clwcon0=qcondc
1750
1751	ELSE ! ok_cvl
1752
1753	CALL conema3 (dtime,
1754	. paprs,pplay,t_seri,q_seri,
1755	. u_seri,v_seri,tr_seri,nbtr,
1756	. ema_work1,ema_work2,
1757	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
1758	. rain_con, snow_con, ibas_con, itop_con,
1759	. upwd,dnwd,dnwd0,bas,top,
1760	. Ma,cape,tvp,rflag,
1761	. pbase
1762	. ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr
1763	. ,clwcon0)
1764
1765	ENDIF ! ok_cvl
1766
1767	IF (.NOT. ok_gust) THEN
1768	do i = 1, klon
1769	wd(i)=0.0
1770	enddo
1771	ENDIF
1772
1773	c =================================================================== c
1774	c Calcul des proprietes des nuages convectifs
1775	c
1776	DO k = 1, klev
1777	DO i = 1, klon
1778	zx_t = t_seri(i,k)
1779	IF (thermcep) THEN
1780	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
1781	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
1782	zx_qs = MIN(0.5,zx_qs)
1783	zcor = 1./(1.-retv*zx_qs)
1784	zx_qs = zx_qs*zcor
1785	ELSE
1786	IF (zx_t.LT.t_coup) THEN
1787	zx_qs = qsats(zx_t)/pplay(i,k)
1788	ELSE
1789	zx_qs = qsatl(zx_t)/pplay(i,k)
1790	ENDIF
1791	ENDIF
1792	zqsat(i,k)=zx_qs
1793	ENDDO
1794	ENDDO
1795
1796	c calcul des proprietes des nuages convectifs
1797	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
1798	call clouds_gno
1799	s (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
1800
1801	c =================================================================== c
1802
1803	DO i = 1, klon
1804	ema_pcb(i) = pbase(i)
1805	ENDDO
1806	DO i = 1, klon
1807	ema_pct(i) = paprs(i,itop_con(i))
1808	ENDDO
1809	DO i = 1, klon
1810	ema_cbmf(i) = ema_workcbmf(i)
1811	ENDDO
1812	ELSE
1813	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
1814	CALL abort
1815	ENDIF
1816
1817	c CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
1818	c . d_u_con, d_v_con)
1819
1820	DO k = 1, klev
1821	DO i = 1, klon
1822	t_seri(i,k) = t_seri(i,k) + d_t_con(i,k)
1823	q_seri(i,k) = q_seri(i,k) + d_q_con(i,k)
1824	u_seri(i,k) = u_seri(i,k) + d_u_con(i,k)
1825	v_seri(i,k) = v_seri(i,k) + d_v_con(i,k)
1826	ENDDO
1827	ENDDO
1828	c
1829	IF (if_ebil.ge.2) THEN
1830	ztit='after convect'
1831	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
1832	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1833	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1834	call diagphy(airephy,ztit,ip_ebil
1835	e , zero_v, zero_v, zero_v, zero_v, zero_v
1836	e , zero_v, rain_con, snow_con, ztsol
1837	e , d_h_vcol, d_qt, d_ec
1838	s , fs_bound, fq_bound )
1839	END IF
1840	C
1841	IF (check) THEN
1842	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
1843	WRITE(lunout,*)"aprescon=", za
1844	zx_t = 0.0
1845	za = 0.0
1846	DO i = 1, klon
1847	za = za + airephy(i)/FLOAT(klon)
1848	zx_t = zx_t + (rain_con(i)+
1849	. snow_con(i))*airephy(i)/FLOAT(klon)
1850	ENDDO
1851	zx_t = zx_t/za*dtime
1852	WRITE(lunout,*)"Precip=", zx_t
1853	ENDIF
1854	IF (zx_ajustq) THEN
1855	DO i = 1, klon
1856	z_apres(i) = 0.0
1857	ENDDO
1858	DO k = 1, klev
1859	DO i = 1, klon
1860	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k))
1861	. *(paprs(i,k)-paprs(i,k+1))/RG
1862	ENDDO
1863	ENDDO
1864	DO i = 1, klon
1865	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime)
1866	. /z_apres(i)
1867	ENDDO
1868	DO k = 1, klev
1869	DO i = 1, klon
1870	IF (z_factor(i).GT.(1.0+1.0E-08) .OR.
1871	. z_factor(i).LT.(1.0-1.0E-08)) THEN
1872	q_seri(i,k) = q_seri(i,k) * z_factor(i)
1873	ENDIF
1874	ENDDO
1875	ENDDO
1876	ENDIF
1877	zx_ajustq=.FALSE.
1878	c
1879	IF (nqmax.GT.2) THEN !--melange convectif de traceurs
1880	c
1881	IF (iflag_con .NE. 2 .AND. debut) THEN
1882	WRITE(lunout,*)'Pour l instant, seul conflx fonctionne ',
1883	$ 'avec traceurs', iflag_con
1884	WRITE(lunout,*)' Mettre iflag_con',
1885	$ ' = 2 dans run.def et repasser'
1886	c CALL abort
1887	ENDIF
1888	c
1889	ENDIF !--nqmax.GT.2
1890	c
1891	c Appeler l'ajustement sec
1892	c
1893	c===================================================================
1894	c Convection seche (thermiques ou ajustement)
1895	c===================================================================
1896	c
1897	d_t_ajs(:,:)=0.
1898	d_u_ajs(:,:)=0.
1899	d_v_ajs(:,:)=0.
1900	d_q_ajs(:,:)=0.
1901	fm_therm(:,:)=0.
1902	entr_therm(:,:)=0.
1903	c
1904	IF(prt_level>9)WRITE(lunout,*)
1905	. 'AVANT LA CONVECTION SECHE , iflag_thermals='
1906	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
1907	if(iflag_thermals.lt.0) then
1908	c Rien
1909	c ====
1910	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
1911	else if(iflag_thermals.eq.0) then
1912
1913	c Ajustement sec
1914	c ==============
1915	IF(prt_level>9)WRITE(lunout,*)'ajsec'
1916	CALL ajsec(paprs, pplay, t_seri,q_seri, d_t_ajs, d_q_ajs)
1917	t_seri(:,:) = t_seri(:,:) + d_t_ajs(:,:)
1918	q_seri(:,:) = q_seri(:,:) + d_q_ajs(:,:)
1919	else
1920	c Thermiques
1921	c ==========
1922	IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals='
1923	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
1924	call calltherm(pdtphys
1925	s ,pplay,paprs,pphi
1926	s ,u_seri,v_seri,t_seri,q_seri
1927	s ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
1928	s ,fm_therm,entr_therm)
1929	endif
1930	c
1931	c===================================================================
1932	c
1933	IF (if_ebil.ge.2) THEN
1934	ztit='after dry_adjust'
1935	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
1936	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1937	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1938	END IF
1939
1940
1941	c-------------------------------------------------------------------------
1942	c Caclul des ratqs
1943	c-------------------------------------------------------------------------
1944
1945	c print*,'calcul des ratqs'
1946	c ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
1947	c ----------------
1948	c on ecrase le tableau ratqsc calcule par clouds_gno
1949	if (iflag_cldcon.eq.1) then
1950	do k=1,klev
1951	do i=1,klon
1952	if(ptconv(i,k)) then
1953	ratqsc(i,k)=ratqsbas
1954	s +fact_cldcon*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
1955	else
1956	ratqsc(i,k)=0.
1957	endif
1958	enddo
1959	enddo
1960	endif
1961
1962	c ratqs stables
1963	c -------------
1964	do k=1,klev
1965	cIM RAJOUT boucle do=i
1966	do i=1, klon
1967	cIM ratqss(:,k)=ratqsbas+(ratqshaut-ratqsbas)*
1968	cIM s min((paprs(:,1)-pplay(:,k))/(paprs(:,1)-30000.),1.)
1969	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
1970	s min((paprs(i,1)-pplay(i,k))/(paprs(i,1)-30000.),1.)
1971	cIM print*,' IMratqs STABLE i, k',i,k,ratqss(i,k)
1972	enddo
1973	enddo
1974
1975
1976	c ratqs final
1977	c -----------
1978	if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2) then
1979	c les ratqs sont une conbinaison de ratqss et ratqsc
1980	c ratqs final
1981	c 1e4 (en gros 3 heures), en dur pour le moment, est le temps de
1982	c relaxation des ratqs
1983	c facttemps=exp(-pdtphys/1.e4)
1984	facteur=exp(-pdtphys*facttemps)
1985	ratqs(:,:)=max(ratqs(:,:)*facteur,ratqss(:,:))
1986	ratqs(:,:)=max(ratqs(:,:),ratqsc(:,:))
1987	c print*,'calcul des ratqs fini'
1988	else
1989	c on ne prend que le ratqs stable pour fisrtilp
1990	ratqs(:,:)=ratqss(:,:)
1991	endif
1992
1993
1994	c
1995	c Appeler le processus de condensation a grande echelle
1996	c et le processus de precipitation
1997	c-------------------------------------------------------------------------
1998	CALL fisrtilp(dtime,paprs,pplay,
1999	. t_seri, q_seri,ptconv,ratqs,
2000	. d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq,
2001	. rain_lsc, snow_lsc,
2002	. pfrac_impa, pfrac_nucl, pfrac_1nucl,
2003	. frac_impa, frac_nucl,
2004	. prfl, psfl, rhcl)
2005
2006	WHERE (rain_lsc < 0) rain_lsc = 0.
2007	WHERE (snow_lsc < 0) snow_lsc = 0.
2008	DO k = 1, klev
2009	DO i = 1, klon
2010	t_seri(i,k) = t_seri(i,k) + d_t_lsc(i,k)
2011	q_seri(i,k) = q_seri(i,k) + d_q_lsc(i,k)
2012	ql_seri(i,k) = ql_seri(i,k) + d_ql_lsc(i,k)
2013	cldfra(i,k) = rneb(i,k)
2014	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
2015	ENDDO
2016	ENDDO
2017	IF (check) THEN
2018	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2019	WRITE(lunout,*)"apresilp=", za
2020	zx_t = 0.0
2021	za = 0.0
2022	DO i = 1, klon
2023	za = za + airephy(i)/FLOAT(klon)
2024	zx_t = zx_t + (rain_lsc(i)
2025	. + snow_lsc(i))*airephy(i)/FLOAT(klon)
2026	ENDDO
2027	zx_t = zx_t/za*dtime
2028	WRITE(lunout,*)"Precip=", zx_t
2029	ENDIF
2030	c
2031	IF (if_ebil.ge.2) THEN
2032	ztit='after fisrt'
2033	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
2034	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2035	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2036	call diagphy(airephy,ztit,ip_ebil
2037	e , zero_v, zero_v, zero_v, zero_v, zero_v
2038	e , zero_v, rain_lsc, snow_lsc, ztsol
2039	e , d_h_vcol, d_qt, d_ec
2040	s , fs_bound, fq_bound )
2041	END IF
2042	c
2043	c-------------------------------------------------------------------
2044	c PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
2045	c-------------------------------------------------------------------
2046
2047	c 1. NUAGES CONVECTIFS
2048	c
2049	IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
2050
2051	c Nuages diagnostiques pour Tiedtke
2052	CALL diagcld1(paprs,pplay,
2053	. rain_con,snow_con,ibas_con,itop_con,
2054	. diafra,dialiq)
2055	DO k = 1, klev
2056	DO i = 1, klon
2057	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2058	cldliq(i,k) = dialiq(i,k)
2059	cldfra(i,k) = diafra(i,k)
2060	ENDIF
2061	ENDDO
2062	ENDDO
2063
2064	ELSE IF (iflag_cldcon.eq.3) THEN
2065	c On prend pour les nuages convectifs le max du calcul de la
2066	c convection et du calcul du pas de temps précédent diminué d'un facteur
2067	c facttemps
2068	c facttemps=pdtphys/1.e4
2069	facteur = pdtphys *facttemps
2070	do k=1,klev
2071	do i=1,klon
2072	rnebcon(i,k)=rnebcon(i,k)*facteur
2073	if (rnebcon0(i,k)clwcon0(i,k).gt.rnebcon(i,k)clwcon(i,k))
2074	s then
2075	rnebcon(i,k)=rnebcon0(i,k)
2076	clwcon(i,k)=clwcon0(i,k)
2077	endif
2078	enddo
2079	enddo
2080
2081	cIM calcul nuages par le simulateur ISCCP
2082	IF (ok_isccp) THEN
2083	cIM calcul tau. emi nuages convectifs
2084	convfra(:,:)=rnebcon(:,:)
2085	convliq(:,:)=rnebcon(:,:)*clwcon(:,:)
2086	CALL newmicro (paprs, pplay,ok_newmicro,
2087	. t_seri, convliq, convfra, dtau_c, dem_c,
2088	. cldh_c, cldl_c, cldm_c, cldt_c, cldq_c,
2089	. flwp_c, fiwp_c, flwc_c, fiwc_c,
2090	e ok_aie,
2091	e sulfate, sulfate_pi,
2092	e bl95_b0, bl95_b1,
2093	s cldtaupi, re, fl)
2094	c
2095	cIM calcul tau. emi nuages startiformes
2096	CALL newmicro (paprs, pplay,ok_newmicro,
2097	. t_seri, cldliq, cldfra, dtau_s, dem_s,
2098	. cldh_s, cldl_s, cldm_s, cldt_s, cldq_s,
2099	. flwp_s, fiwp_s, flwc_s, fiwc_s,
2100	e ok_aie,
2101	e sulfate, sulfate_pi,
2102	e bl95_b0, bl95_b1,
2103	s cldtaupi, re, fl)
2104	c
2105	cldtot(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
2106
2107	cIM inversion des niveaux de pression ==> de haut en bas
2108	CALL haut2bas(klon, klev, pplay, pfull)
2109	CALL haut2bas(klon, klev, q_seri, qv)
2110	CALL haut2bas(klon, klev, cldtot, cc)
2111	CALL haut2bas(klon, klev, rnebcon, conv)
2112	CALL haut2bas(klon, klev, dtau_s, dtau_sH2B)
2113	CALL haut2bas(klon, klev, dtau_c, dtau_cH2B)
2114	CALL haut2bas(klon, klev, t_seri, at)
2115	CALL haut2bas(klon, klev, dem_s, dem_sH2B)
2116	CALL haut2bas(klon, klev, dem_c, dem_cH2B)
2117	CALL haut2bas(klon, klevp1, paprs, phalf)
2118
2119	c open(99,file='tautab.bin',access='sequential',
2120	c $ form='unformatted',status='old')
2121	c read(99) tautab
2122
2123	cIM210503
2124	IF (debut) THEN
2125	open(99,file='tautab.formatted', FORM='FORMATTED')
2126	read(99,'(f30.20)') tautab
2127	close(99)
2128	c
2129	open(99,file='invtau.formatted',form='FORMATTED')
2130	read(99,'(i10)') invtau
2131	close(99)
2132	c
2133	cIM: calcul coordonnees regions pour statistiques distribution
2134	cIM: nuages en ftion du regime dynamique pour regions oceaniques
2135	IF (ok_regdyn) THEN !histREGDYN
2136	nsrf=3
2137	DO nreg=1, nbregdyn
2138	DO i=1, klon
2139
2140	c IF (debut) THEN
2141	IF(rlon(i).LT.0.) THEN
2142	rlonPOS(i)=rlon(i)+360.
2143	ELSE
2144	rlonPOS(i)=rlon(i)
2145	ENDIF
2146	c ENDIF
2147
2148	pct_ocean(i,nreg)=0
2149
2150	c test si c'est 1 point d'ocean
2151	IF(pctsrf(i,nsrf).EQ.1.) THEN
2152
2153	IF(nreg.EQ.1) THEN
2154
2155	c TROP
2156	IF(rlat(i).GE.-30.AND.rlat(i).LE.30.) THEN
2157	pct_ocean(i,nreg)=1
2158	ENDIF
2159
2160	c PACIFIQUE NORD
2161	ELSEIF(nreg.EQ.2) THEN
2162	IF(rlat(i).GE.40.AND.rlat(i).LE.60.) THEN
2163	IF(rlonPOS(i).GE.160..AND.rlonPOS(i).LE.235.) THEN
2164	pct_ocean(i,nreg)=1
2165	ENDIF
2166	ENDIF
2167	c CALIFORNIE ST-CU
2168	ELSEIF(nreg.EQ.3) THEN
2169	IF(rlonPOS(i).GE.220..AND.rlonPOS(i).LE.250.) 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 HAWAI
2175	ELSEIF(nreg.EQ.4) THEN
2176	IF(rlonPOS(i).GE.180..AND.rlonPOS(i).LE.220.) THEN
2177	IF(rlat(i).GE.15.AND.rlat(i).LE.35.) THEN
2178	pct_ocean(i,nreg)=1
2179	ENDIF
2180	ENDIF
2181	c WARM POOL
2182	ELSEIF(nreg.EQ.5) THEN
2183	IF(rlonPOS(i).GE.70..AND.rlonPOS(i).LE.150.) THEN
2184	IF(rlat(i).GE.-5.AND.rlat(i).LE.20.) THEN
2185	pct_ocean(i,nreg)=1
2186	ENDIF
2187	ENDIF
2188	ENDIF !nbregdyn
2189	c TROP
2190	c IF(rlat(i).GE.-30.AND.rlat(i).LE.30.) THEN
2191	c pct_ocean(i)=.TRUE.
2192	c WRITE(,) 'pct_ocean =',i, rlon(i), rlat(i)
2193	c ENDIF !lon
2194	c ENDIF !lat
2195
2196	ENDIF !pctsrf
2197	ENDDO !klon
2198	ENDDO !nbregdyn
2199	ENDIF !ok_regdyn
2200
2201	cIM somme de toutes les nhistoW BEG
2202	DO nreg = 1, nbregdyn
2203	DO k = 1, kmaxm1
2204	DO l = 1, lmaxm1
2205	DO iw = 1, iwmax
2206	nhistoWt(k,l,iw,nreg)=0.
2207	ENDDO !iw
2208	ENDDO !l
2209	ENDDO !k
2210	ENDDO !nreg
2211	cIM somme de toutes les nhistoW END
2212	ENDIF
2213	cIM: initialisation de seed
2214	DO i=1, klon
2215	seed(i)=i+100
2216	ENDDO
2217
2218	cIM: pas de debug, debugcol
2219	debug=0
2220	debugcol=0
2221	cIM260503
2222	c o500 ==> distribution nuage ftion du regime dynamique a 500 hPa
2223	DO k=1, klevm1
2224	kp1=k+1
2225	c PRINT*,'k, presnivs',k,presnivs(k), presnivs(kp1)
2226	if(presnivs(k).GT.50000.AND.presnivs(kp1).LT.50000.) THEN
2227	DO i=1, klon
2228	o500(i)=omega(i,k)*RDAY/100.
2229	c if(i.EQ.1) print*,' 500hPa lev',k,presnivs(k),presnivs(kp1)
2230	ENDDO
2231	GOTO 1000
2232	endif
2233	1000 continue
2234	ENDDO
2235
2236	CALL ISCCP_CLOUD_TYPES(
2237	& debug,
2238	& debugcol,
2239	& klon,
2240	& sunlit,
2241	& klev,
2242	& ncol,
2243	& seed,
2244	& pfull,
2245	& phalf,
2246	& qv, cc, conv, dtau_sH2B, dtau_cH2B,
2247	& top_height,
2248	& overlap,
2249	& tautab,
2250	& invtau,
2251	& ztsol,
2252	& emsfc_lw,
2253	& at, dem_sH2B, dem_cH2B,
2254	& fq_isccp,
2255	& totalcldarea,
2256	& meanptop,
2257	& meantaucld,
2258	& boxtau,
2259	& boxptop)
2260
2261
2262	c passage de la grille (klon,7,7) a (iim,jjmp1,7,7)
2263	DO l=1, lmaxm1
2264	DO k=1, kmaxm1
2265	DO i=1, iim
2266	fq4d(i,1,k,l)=fq_isccp(1,k,l)
2267	ENDDO
2268	DO j=2, jjm
2269	DO i=1, iim
2270	ig=i+1+(j-2)*iim
2271	fq4d(i,j,k,l)=fq_isccp(ig,k,l)
2272	ENDDO
2273	ENDDO
2274	DO i=1, iim
2275	fq4d(i,jjmp1,k,l)=fq_isccp(klon,k,l)
2276	ENDDO
2277	ENDDO
2278	ENDDO
2279	c
2280	DO l=1, lmaxm1
2281	DO k=1, kmaxm1
2282	DO j=1, jjmp1
2283	DO i=1, iim
2284	ni=(i-1)*lmaxm1+l
2285	nj=(j-1)*kmaxm1+k
2286	fq3d(ni,nj)=fq4d(i,j,k,l)
2287	ENDDO
2288	ENDDO
2289	ENDDO
2290	ENDDO
2291
2292	c
2293	c calculs statistiques distribution nuage ftion du regime dynamique
2294	c
2295	c Ce calcul doit etre fait a partir de valeurs mensuelles ??
2296	CALL histo_o500_pctau(nbregdyn,pct_ocean,o500,fq_isccp,
2297	&histoW,nhistoW)
2298	c
2299	c nhistoWt = somme de toutes les nhistoW
2300	DO nreg=1, nbregdyn
2301	DO k = 1, kmaxm1
2302	DO l = 1, lmaxm1
2303	DO iw = 1, iwmax
2304	nhistoWt(k,l,iw,nreg)=nhistoWt(k,l,iw,nreg)+
2305	& nhistoW(k,l,iw,nreg)
2306	ENDDO
2307	ENDDO
2308	ENDDO
2309	ENDDO
2310	c
2311	ENDIF !ok_isccp
2312
2313	c On prend la somme des fractions nuageuses et des contenus en eau
2314	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
2315	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
2316
2317	ENDIF
2318
2319	c
2320	c 2. NUAGES STARTIFORMES
2321	c
2322	IF (ok_stratus) THEN
2323	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
2324	DO k = 1, klev
2325	DO i = 1, klon
2326	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2327	cldliq(i,k) = dialiq(i,k)
2328	cldfra(i,k) = diafra(i,k)
2329	ENDIF
2330	ENDDO
2331	ENDDO
2332	ENDIF
2333	c
2334	c Precipitation totale
2335	c
2336	DO i = 1, klon
2337	rain_fall(i) = rain_con(i) + rain_lsc(i)
2338	snow_fall(i) = snow_con(i) + snow_lsc(i)
2339	ENDDO
2340	c
2341	IF (if_ebil.ge.2) THEN
2342	ztit="after diagcld"
2343	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
2344	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2345	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2346	END IF
2347	c
2348	c Calculer l'humidite relative pour diagnostique
2349	c
2350	DO k = 1, klev
2351	DO i = 1, klon
2352	zx_t = t_seri(i,k)
2353	IF (thermcep) THEN
2354	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2355	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2356	zx_qs = MIN(0.5,zx_qs)
2357	zcor = 1./(1.-retv*zx_qs)
2358	zx_qs = zx_qs*zcor
2359	ELSE
2360	IF (zx_t.LT.t_coup) THEN
2361	zx_qs = qsats(zx_t)/pplay(i,k)
2362	ELSE
2363	zx_qs = qsatl(zx_t)/pplay(i,k)
2364	ENDIF
2365	ENDIF
2366	zx_rh(i,k) = q_seri(i,k)/zx_qs
2367	zqsat(i,k)=zx_qs
2368	ENDDO
2369	ENDDO
2370	cjq - introduce the aerosol direct and first indirect radiative forcings
2371	cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
2372	IF (ok_ade.OR.ok_aie) THEN
2373	! Get sulfate aerosol distribution
2374	CALL readsulfate(rjourvrai, debut, sulfate)
2375	CALL readsulfate_preind(rjourvrai, debut, sulfate_pi)
2376
2377	! Calculate aerosol optical properties (Olivier Boucher)
2378	CALL aeropt(pplay, paprs, t_seri, sulfate, rhcl,
2379	. tau_ae, piz_ae, cg_ae, aerindex)
2380	cym
2381	ELSE
2382	tau_ae(:,:,:)=0.0
2383	piz_ae(:,:,:)=0.0
2384	cg_ae(:,:,:)=0.0
2385	cym
2386	ENDIF
2387
2388	#ifdef INCA
2389	calday = FLOAT(julien) + gmtime
2390
2391	#ifdef INCA_AER
2392	call AEROSOL_METEO_CALC(calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs,
2393	& prfl,psfl,pctsrf,airephy,xjour,rlat,rlon)
2394	#endif
2395
2396	#ifdef INCAINFO
2397	WRITE(lunout,*)'Appel CHEMHOOK_BEGIN ...'
2398	#endif
2399	CALL chemhook_begin (calday,
2400	$ pctsrf(1,1),
2401	$ rlat,
2402	$ rlon,
2403	$ airephy,
2404	$ paprs,
2405	$ pplay,
2406	$ ycoefh,
2407	$ pphi,
2408	$ t_seri,
2409	$ u,
2410	$ v,
2411	$ wo,
2412	$ q_seri,
2413	$ zxtsol,
2414	$ zxsnow,
2415	$ solsw,
2416	$ albsol,
2417	$ rain_fall,
2418	$ snow_fall,
2419	$ itop_con,
2420	$ ibas_con,
2421	$ cldfra,
2422	$ iim,
2423	$ jjm,
2424	$ tr_seri)
2425	#ifdef INCAINFO
2426	WRITE(lunout,*)'OK.'
2427	#endif
2428	#endif
2429	c
2430	c Calculer les parametres optiques des nuages et quelques
2431	c parametres pour diagnostiques:
2432	c
2433	if (ok_newmicro) then
2434	CALL newmicro (paprs, pplay,ok_newmicro,
2435	. t_seri, cldliq, cldfra, cldtau, cldemi,
2436	. cldh, cldl, cldm, cldt, cldq,
2437	. flwp, fiwp, flwc, fiwc,
2438	e ok_aie,
2439	e sulfate, sulfate_pi,
2440	e bl95_b0, bl95_b1,
2441	s cldtaupi, re, fl)
2442	else
2443	CALL nuage (paprs, pplay,
2444	. t_seri, cldliq, cldfra, cldtau, cldemi,
2445	. cldh, cldl, cldm, cldt, cldq,
2446	e ok_aie,
2447	e sulfate, sulfate_pi,
2448	e bl95_b0, bl95_b1,
2449	s cldtaupi, re, fl)
2450
2451	endif
2452	c
2453	c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
2454	c
2455	IF (MOD(itaprad,radpas).EQ.0) THEN
2456	DO i = 1, klon
2457	albsol(i) = falbe(i,is_oce) * pctsrf(i,is_oce)
2458	. + falbe(i,is_lic) * pctsrf(i,is_lic)
2459	. + falbe(i,is_ter) * pctsrf(i,is_ter)
2460	. + falbe(i,is_sic) * pctsrf(i,is_sic)
2461	albsollw(i) = falblw(i,is_oce) * pctsrf(i,is_oce)
2462	. + falblw(i,is_lic) * pctsrf(i,is_lic)
2463	. + falblw(i,is_ter) * pctsrf(i,is_ter)
2464	. + falblw(i,is_sic) * pctsrf(i,is_sic)
2465	ENDDO
2466	CALL radlwsw ! nouveau rayonnement (compatible Arpege-IFS)
2467	e (dist, rmu0, fract,
2468	e paprs, pplay,zxtsol,albsol, albsollw, t_seri,q_seri,
2469	e wo,
2470	e cldfra, cldemi, cldtau,
2471	s heat,heat0,cool,cool0,radsol,albpla,
2472	s topsw,toplw,solsw,sollw,
2473	s sollwdown,
2474	s topsw0,toplw0,solsw0,sollw0,
2475	s lwdn0, lwdn, lwup0, lwup,
2476	s swdn0, swdn, swup0, swup,
2477	e ok_ade, ok_aie, ! new for aerosol radiative effects
2478	e tau_ae, piz_ae, cg_ae, ! ="=
2479	s topswad, solswad, ! ="=
2480	e cldtaupi, ! ="=
2481	s topswai, solswai) ! ="=
2482	itaprad = 0
2483	ENDIF
2484	itaprad = itaprad + 1
2485	c
2486	c Ajouter la tendance des rayonnements (tous les pas)
2487	c
2488	DO k = 1, klev
2489	DO i = 1, klon
2490	t_seri(i,k) = t_seri(i,k)
2491	. + (heat(i,k)-cool(i,k)) * dtime/86400.
2492	ENDDO
2493	ENDDO
2494	c
2495	IF (if_ebil.ge.2) THEN
2496	ztit='after rad'
2497	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
2498	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2499	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2500	call diagphy(airephy,ztit,ip_ebil
2501	e , topsw, toplw, solsw, sollw, zero_v
2502	e , zero_v, zero_v, zero_v, ztsol
2503	e , d_h_vcol, d_qt, d_ec
2504	s , fs_bound, fq_bound )
2505	END IF
2506	c
2507	c
2508	c Calculer l'hydrologie de la surface
2509	c
2510	c CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
2511	c . agesno, ftsol,fqsurf,fsnow, ruis)
2512	c
2513	DO i = 1, klon
2514	zxqsurf(i) = 0.0
2515	zxsnow(i) = 0.0
2516	ENDDO
2517	DO nsrf = 1, nbsrf
2518	DO i = 1, klon
2519	zxqsurf(i) = zxqsurf(i) + fqsurf(i,nsrf)*pctsrf(i,nsrf)
2520	zxsnow(i) = zxsnow(i) + fsnow(i,nsrf)*pctsrf(i,nsrf)
2521	ENDDO
2522	ENDDO
2523	c
2524	c Si une sous-fraction n'existe pas, elle prend la valeur moyenne
2525	c
2526	cXXX DO nsrf = 1, nbsrf
2527	cXXX DO i = 1, klon
2528	cXXX IF (pctsrf(i,nsrf).LT.epsfra) THEN
2529	cXXX fqsurf(i,nsrf) = zxqsurf(i)
2530	cXXX fsnow(i,nsrf) = zxsnow(i)
2531	cXXX ENDIF
2532	cXXX ENDDO
2533	cXXX ENDDO
2534	c
2535	c Calculer le bilan du sol et la derive de temperature (couplage)
2536	c
2537	DO i = 1, klon
2538	c bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
2539	c a la demande de JLD
2540	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
2541	ENDDO
2542	c
2543	cmoddeblott(jan95)
2544	c Appeler le programme de parametrisation de l'orographie
2545	c a l'echelle sous-maille:
2546	c
2547	IF (ok_orodr) THEN
2548	c
2549	c selection des points pour lesquels le shema est actif:
2550	igwd=0
2551	DO i=1,klon
2552	itest(i)=0
2553	c IF ((zstd(i).gt.10.0)) THEN
2554	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
2555	itest(i)=1
2556	igwd=igwd+1
2557	idx(igwd)=i
2558	ENDIF
2559	ENDDO
2560	c igwdim=MAX(1,igwd)
2561	c
2562	CALL drag_noro(klon,klev,dtime,paprs,pplay,
2563	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
2564	e igwd,idx,itest,
2565	e t_seri, u_seri, v_seri,
2566	s zulow, zvlow, zustr, zvstr,
2567	s d_t_oro, d_u_oro, d_v_oro)
2568	c
2569	c ajout des tendances
2570	DO k = 1, klev
2571	DO i = 1, klon
2572	t_seri(i,k) = t_seri(i,k) + d_t_oro(i,k)
2573	u_seri(i,k) = u_seri(i,k) + d_u_oro(i,k)
2574	v_seri(i,k) = v_seri(i,k) + d_v_oro(i,k)
2575	ENDDO
2576	ENDDO
2577	c
2578	ENDIF ! fin de test sur ok_orodr
2579	c
2580	IF (ok_orolf) THEN
2581	c
2582	c selection des points pour lesquels le shema est actif:
2583	igwd=0
2584	DO i=1,klon
2585	itest(i)=0
2586	IF ((zpic(i)-zmea(i)).GT.100.) THEN
2587	itest(i)=1
2588	igwd=igwd+1
2589	idx(igwd)=i
2590	ENDIF
2591	ENDDO
2592	c igwdim=MAX(1,igwd)
2593	c
2594	CALL lift_noro(klon,klev,dtime,paprs,pplay,
2595	e rlat,zmea,zstd,zpic,
2596	e itest,
2597	e t_seri, u_seri, v_seri,
2598	s zulow, zvlow, zustr, zvstr,
2599	s d_t_lif, d_u_lif, d_v_lif)
2600	c
2601	c ajout des tendances
2602	DO k = 1, klev
2603	DO i = 1, klon
2604	t_seri(i,k) = t_seri(i,k) + d_t_lif(i,k)
2605	u_seri(i,k) = u_seri(i,k) + d_u_lif(i,k)
2606	v_seri(i,k) = v_seri(i,k) + d_v_lif(i,k)
2607	ENDDO
2608	ENDDO
2609	c
2610	ENDIF ! fin de test sur ok_orolf
2611	c
2612	IF (if_ebil.ge.2) THEN
2613	ztit='after orography'
2614	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
2615	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2616	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2617	END IF
2618	c
2619	c
2620	cAA
2621	cAA Installation de l'interface online-offline pour traceurs
2622	cAA
2623	c====================================================================
2624	c Calcul des tendances traceurs
2625	c====================================================================
2626	C Pascale : il faut quand meme apeller phytrac car il gere les sorties
2627	cKE43 des traceurs => il faut donc mettre des flags a .false.
2628	IF (iflag_con.GE.3) THEN
2629	c on ajoute les tendances calculees par KE43
2630	cXXX OM on onhibe la convection sur les traceurs
2631	DO iq=1, nqmax-2 ! Sandrine a -3 ???
2632	cXXX OM on inhibe la convection sur les traceur
2633	cXXX DO k = 1, nlev
2634	cXXX DO i = 1, klon
2635	cXXX tr_seri(i,k,iq) = tr_seri(i,k,iq) + d_tr(i,k,iq)
2636	cXXX ENDDO
2637	cXXX ENDDO
2638	WRITE(iqn,'(i2.2)') iq
2639	CALL minmaxqfi(tr_seri(1,1,iq),0.,1.e33,'couche lim iq='//iqn)
2640	ENDDO
2641	CMAF modif pour garder info du nombre de traceurs auxquels
2642	C la physique s'applique
2643	ELSE
2644	CMAF modif pour garder info du nombre de traceurs auxquels
2645	C la physique s'applique
2646	C
2647	call phytrac (rnpb,
2648	I itap, julien, gmtime,
2649	I debut,lafin,
2650	I nqmax-2,
2651	I nlon,nlev,dtime,
2652	I u,v,t,paprs,pplay,
2653	I pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
2654	I ycoefh,fm_therm,entr_therm,yu1,yv1,ftsol,
2655	I pctsrf,rlat,frac_impa, frac_nucl,
2656	I rlon,presnivs,pphis,pphi,
2657	I albsol,
2658	I qx(1,1,1), rhcl,
2659	I cldfra, rneb, diafra, cldliq, itop_con,
2660	I ibas_con,
2661	I pmflxr,pmflxs,prfl,psfl,
2662	#ifdef INCA_CH4
2663	I flxmass_w,
2664	#endif
2665	O tr_seri)
2666	ENDIF
2667
2668	IF (offline) THEN
2669
2670	print*,'Attention on met a 0 les thermiques pour phystoke'
2671	call phystokenc (
2672	I nlon,nlev,pdtphys,rlon,rlat,
2673	I t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
2674	I fm_therm,entr_therm,
2675	I ycoefh,yu1,yv1,ftsol,pctsrf,
2676	I frac_impa, frac_nucl,
2677	I pphis,airephy,dtime,itap)
2678
2679
2680	ENDIF
2681
2682	c
2683	c Calculer le transport de l'eau et de l'energie (diagnostique)
2684	c
2685	CALL transp (paprs,zxtsol,
2686	e t_seri, q_seri, u_seri, v_seri, zphi,
2687	s ve, vq, ue, uq)
2688	c
2689	c Accumuler les variables a stocker dans les fichiers histoire:
2690	c
2691	c
2692	c
2693	c+jld ec_conser
2694	DO k = 1, klev
2695	DO i = 1, klon
2696	ZRCPD = RCPD(1.0+RVTMP2q_seri(i,k))
2697	d_t_ec(i,k)=0.5/ZRCPD
2698	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
2699	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
2700	d_t_ec(i,k) = d_t_ec(i,k)/dtime
2701	END DO
2702	END DO
2703	c-jld ec_conser
2704	IF (if_ebil.ge.1) THEN
2705	ztit='after physic'
2706	CALL diagetpq(airephy,ztit,ip_ebil,1,1,dtime
2707	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2708	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2709	C Comme les tendances de la physique sont ajoute dans la dynamique,
2710	C on devrait avoir que la variation d'entalpie par la dynamique
2711	C est egale a la variation de la physique au pas de temps precedent.
2712	C Donc la somme de ces 2 variations devrait etre nulle.
2713	call diagphy(airephy,ztit,ip_ebil
2714	e , topsw, toplw, solsw, sollw, sens
2715	e , evap, rain_fall, snow_fall, ztsol
2716	e , d_h_vcol, d_qt, d_ec
2717	s , fs_bound, fq_bound )
2718	C
2719	d_h_vcol_phy=d_h_vcol
2720	C
2721	END IF
2722	C
2723	c=======================================================================
2724	c SORTIES
2725	c=======================================================================
2726
2727	c Interpollation sur quelques niveaux de pression
2728	c -----------------------------------------------
2729	c
2730	c on moyenne mensuellement les champs 3D et on les interpole sur les niveaux STD
2731	c if(itap.EQ.1.OR.itap.EQ.13.OR.itap.EQ.25.OR.itap.EQ.37) THEN
2732	c if(MOD(itap,12).EQ.1) THEN
2733	cIM 120304 END
2734	DO k=1, nlevSTD
2735	call plevel(klon,klev,.true.,pplay,rlevSTD(k),
2736	. t_seri,tlevSTD(:,k))
2737	call plevel(klon,klev,.false.,pplay,rlevSTD(k),
2738	. u_seri,ulevSTD(:,k))
2739	call plevel(klon,klev,.false.,pplay,rlevSTD(k),
2740	. v_seri,vlevSTD(:,k))
2741	call plevel(klon,klev,.false.,pplay,rlevSTD(k),
2742	. zphi,philevSTD(:,k))
2743	call plevel(klon,klev,.false.,pplay,rlevSTD(k),
2744	. qx(:,:,ivap),qlevSTD(:,k))
2745	call plevel(klon,klev,.false.,pplay,rlevSTD(k),
2746	. zx_rh,rhlevSTD(:,k))
2747	ENDDO !nlevSTD
2748	c ENSEMBLES BEG
2749	DO k=1, nlevENS
2750	cIM 170304
2751	tlev(:,k)=tlevSTD(:,indENS(k))
2752	ulev(:,k)=ulevSTD(:,indENS(k))
2753	vlev(:,k)=vlevSTD(:,indENS(k))
2754	philev(:,k)=philevSTD(:,indENS(k))
2755	qlev(:,k)=qlevSTD(:,indENS(k))
2756	rhlev(:,k)=rhlevSTD(:,indENS(k))
2757	c
2758	call plevel(klon,klevp1,.true.,paprs,rlevENS(k),
2759	. omega,wlev(:,k))
2760	c
2761	ENDDO !k=1, nlevENS
2762	cIM 100304 BEG
2763	cIM interpolation a chaque pas de temps du SWup(clr) et SWdn(clr) a 200 hPa
2764	call plevel(klon,klevp1,.true.,paprs,20000.,
2765	$ swdn0,SWdn200clr)
2766	call plevel(klon,klevp1,.false.,paprs,20000.,
2767	$ swdn,SWdn200)
2768	call plevel(klon,klevp1,.false.,paprs,20000.,
2769	$ swup0,SWup200clr)
2770	call plevel(klon,klevp1,.false.,paprs,20000.,
2771	$ swup,SWup200)
2772	c
2773	call plevel(klon,klevp1,.false.,paprs,20000.,
2774	$ lwdn0,LWdn200clr)
2775	call plevel(klon,klevp1,.false.,paprs,20000.,
2776	$ lwdn,LWdn200)
2777	call plevel(klon,klevp1,.false.,paprs,20000.,
2778	$ lwup0,LWup200clr)
2779	call plevel(klon,klevp1,.false.,paprs,20000.,
2780	$ lwup,LWup200)
2781	c
2782	cIM 100304 END
2783	c
2784	c ENSEMBLES END
2785	c
2786	c slp sea level pressure
2787	slp(:) = paprs(:,1)exp(pphis(:)/(RDt_seri(:,1)))
2788	c
2789	ccc prw = eau precipitable
2790	DO i = 1, klon
2791	prw(i) = 0.
2792	DO k = 1, klev
2793	prw(i) = prw(i) +
2794	. q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
2795	ENDDO
2796	ENDDO
2797	c
2798	cIM sorties bilans energie cinetique et potentielle MJO
2799	DO k = 1, klev
2800	DO i = 1, klon
2801	d_u_oli(i,k) = d_u_oro(i,k) + d_u_lif(i,k)
2802	d_v_oli(i,k) = d_v_oro(i,k) + d_v_lif(i,k)
2803	ENDDO
2804	ENDDO
2805	c
2806	IF (MOD(itap-1,lmt_pas) .EQ. 0) THEN
2807	cIM PRINT *,' PHYS cond julien ',julien
2808	c CALL ozonecm( FLOAT(julien), rlat, paprs, wo)
2809	DO i = 1, klon
2810	total_rain(i)=rain_fall(i)+snow_fall(i)
2811	IF(total_rain(i).GT.0.) nday_rain(i)=nday_rain(i)+1.
2812	ENDDO
2813	c
2814	ENDIF
2815	c surface terre
2816	IF (debut) THEN
2817	DO i=1, klon
2818	IF(pctsrf_new(i,is_ter).GT.0.) THEN
2819	paire_ter(i)=airephy(i)*pctsrf_new(i,is_ter)
2820	ENDIF
2821	ENDDO
2822	ENDIF
2823	cIM 050204 END
2824
2825	c=============================================================
2826	c
2827	c Convertir les incrementations en tendances
2828	c
2829	DO k = 1, klev
2830	DO i = 1, klon
2831	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
2832	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
2833	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
2834	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
2835	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
2836	ENDDO
2837	ENDDO
2838	c
2839	IF (nqmax.GE.3) THEN
2840	DO iq = 3, nqmax
2841	DO k = 1, klev
2842	DO i = 1, klon
2843	d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
2844	ENDDO
2845	ENDDO
2846	ENDDO
2847	ENDIF
2848	c
2849	c Sauvegarder les valeurs de t et q a la fin de la physique:
2850	c
2851	DO k = 1, klev
2852	DO i = 1, klon
2853	t_ancien(i,k) = t_seri(i,k)
2854	q_ancien(i,k) = q_seri(i,k)
2855	ENDDO
2856	ENDDO
2857	c
2858	c 22.03.04 BEG
2859	c=============================================================
2860	c Ecriture des sorties
2861	c=============================================================
2862	#ifdef CPP_IOIPSL
2863
2864	#ifdef histhf
2865	#include "write_histhf.h"
2866	#endif
2867
2868	#ifdef histday
2869	#include "write_histday.h"
2870	#endif
2871
2872	#ifdef histmth
2873	#include "write_histmth.h"
2874	#endif
2875
2876	#ifdef histins
2877	#include "write_histins.h"
2878	#endif
2879
2880	#ifdef histREGDYN
2881	#include "write_histREGDYN.h"
2882	#endif
2883
2884	#ifdef histISCCP
2885	#include "write_histISCCP.h"
2886	#endif
2887
2888	#ifdef histmthNMC
2889	#include "write_histmthNMC.h"
2890	#endif
2891
2892	#endif
2893
2894	#ifdef INCA
2895	#ifdef INCAINFO
2896	WRITE(lunout,*)'Appel CHEMHOOK_END ...'
2897	#endif
2898	CALL chemhook_end (calday,
2899	$ dtime,
2900	$ pplay,
2901	$ t_seri,
2902	$ tr_seri,
2903	$ nbtr,
2904	$ paprs,
2905	$ q_seri,
2906	$ annee_ref,
2907	$ day_ini,
2908	$ xjour)
2909	#ifdef INCAINFO
2910	WRITE(lunout,*)'OK.'
2911	#endif
2912	#endif
2913
2914	c 22.03.04 END
2915	c
2916	c====================================================================
2917	c Si c'est la fin, il faut conserver l'etat de redemarrage
2918	c====================================================================
2919	c
2920	IF (lafin) THEN
2921	itau_phy = itau_phy + itap
2922	ccc IF (ok_oasis) CALL quitcpl
2923	CALL phyredem ("restartphy.nc",dtime,radpas,
2924	. rlat, rlon, pctsrf, ftsol, ftsoil, deltat, fqsurf, qsol,
2925	. fsnow, falbe,falblw, fevap, rain_fall, snow_fall,
2926	. solsw, sollwdown,dlw,
2927	. radsol,frugs,agesno,
2928	. zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,
2929	. t_ancien, q_ancien, rnebcon, ratqs, clwcon,run_off_lic_0)
2930	ENDIF
2931
2932
2933	RETURN
2934	END
2935	FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
2936	IMPLICIT none
2937	c
2938	c Calculer et imprimer l'eau totale. A utiliser pour verifier
2939	c la conservation de l'eau
2940	c
2941	#include "YOMCST.h"
2942	INTEGER klon,klev
2943	REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
2944	REAL aire(klon)
2945	REAL qtotal, zx, qcheck
2946	INTEGER i, k
2947	c
2948	zx = 0.0
2949	DO i = 1, klon
2950	zx = zx + aire(i)
2951	ENDDO
2952	qtotal = 0.0
2953	DO k = 1, klev
2954	DO i = 1, klon
2955	qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i)
2956	. *(paprs(i,k)-paprs(i,k+1))/RG
2957	ENDDO
2958	ENDDO
2959	c
2960	qcheck = qtotal/zx
2961	c
2962	RETURN
2963	END
2964	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
2965	IMPLICIT none
2966	c
2967	c Tranformer une variable de la grille physique a
2968	c la grille d'ecriture
2969	c
2970	INTEGER nfield,nlon,iim,jjmp1, jjm
2971	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
2972	c
2973	INTEGER i, n, ig
2974	c
2975	jjm = jjmp1 - 1
2976	DO n = 1, nfield
2977	DO i=1,iim
2978	ecrit(i,n) = fi(1,n)
2979	ecrit(i+jjm*iim,n) = fi(nlon,n)
2980	ENDDO
2981	DO ig = 1, nlon - 2
2982	ecrit(iim+ig,n) = fi(1+ig,n)
2983	ENDDO
2984	ENDDO
2985	RETURN
2986	END

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

Download in other formats:

Original Format