Context Navigation

physiq.F @ 634

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

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

Download in other formats:

Original Format