Context Navigation

physiq.F @ 486

Last change on this file since 486 was 486, checked in by lmdzadmin, 20 years ago
Phasage avec la version de Ionela IM/LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 80.1 KB

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

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

Download in other formats:

Original Format