Context Navigation

physiq.F @ 567

Last change on this file since 567 was 567, checked in by lmdzadmin, 19 years ago
Corrections de bugs LMDZ4/INCA MAF LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 90.9 KB

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

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

Download in other formats:

Original Format