Context Navigation

physiq.F @ 616

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

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

Download in other formats:

Original Format