Context Navigation

physiq.F @ 541

Last change on this file since 541 was 541, checked in by lmdzadmin, 20 years ago

Convergence avec la version d'Olivia Coindreau incluant:

le offline
les thermiques
mellor & yamada dans la couche limite

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 90.2 KB

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

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

Download in other formats:

Original Format