/[lmdze]/trunk/libf/phylmd/physiq.f90
ViewVC logotype

Annotation of /trunk/libf/phylmd/physiq.f90

Parent Directory Parent Directory | Revision Log Revision Log

Revision 15 - (hide annotations)
Fri Aug 1 15:24:12 2008 UTC (15 years, 9 months ago) by guez
File size: 78386 byte(s) 
-- Minor modification of input/output:

Added variable "Sigma_O3_Royer" to "histday.nc". "ecrit_day" is not
modified in "physiq". Removed variables "pyu1", "pyv1", "ftsol1",
"ftsol2", "ftsol3", "ftsol4", "psrf1", "psrf2", "psrf3", "psrf4"
"mfu", "mfd", "en_u", "en_d", "de_d", "de_u", "coefh" from
"histrac.nc".

Variable "raz_date" of module "conf_gcm_m" has logical type instead of
integer type.

-- Should not change any result at run time:

Modified calls to "IOIPSL_Lionel" procedures because the interfaces of
these procedures have been simplified.

Changed name of variable in module "start_init_orog_m": "masque" to
"mask".

Created a module containing procedure "phyredem".

Removed arguments "punjours", "pdayref" and "ptimestep" of procedure
"iniphysiq".

Renamed procedure "gr_phy_write" to "gr_phy_write_2d". Created
procedure "gr_phy_write_3d".

Removed procedures "ini_undefstd", "moy_undefSTD", "calcul_STDlev",
"calcul_divers".
1 guez 3 module physiq_m
2    
3     ! This module is clean: no C preprocessor directive, no include line.
4    
5     IMPLICIT none
6    
7     private
8     public physiq
9    
10     contains
11    
12 guez 12 SUBROUTINE physiq(nq, firstcal, lafin, rdayvrai, gmtime, pdtphys, paprs, &
13 guez 13 pplay, pphi, pphis, presnivs, u, v, t, qx, omega, d_u, d_v, &
14 guez 3 d_t, d_qx, d_ps, dudyn, PVteta)
15    
16     ! From phylmd/physiq.F, v 1.22 2006/02/20 09:38:28
17    
18     ! Author : Z.X. Li (LMD/CNRS), date: 1993/08/18
19    
20     ! Objet: Moniteur general de la physique du modele
21     !AA Modifications quant aux traceurs :
22     !AA - uniformisation des parametrisations ds phytrac
23     !AA - stockage des moyennes des champs necessaires
24     !AA en mode traceur off-line
25    
26     USE ioipsl, only: ymds2ju, histwrite, histsync
27     use dimens_m, only: jjm, iim, llm
28     use indicesol, only: nbsrf, is_ter, is_lic, is_sic, is_oce, &
29     clnsurf, epsfra
30     use dimphy, only: klon, nbtr
31     use conf_gcm_m, only: raz_date, offline, iphysiq
32     use dimsoil, only: nsoilmx
33     use temps, only: itau_phy, day_ref, annee_ref, itaufin
34 guez 15 use clesphys, only: ecrit_hf, ecrit_ins, ecrit_mth, &
35 guez 12 cdmmax, cdhmax, &
36     co2_ppm, ecrit_reg, ecrit_tra, ksta, ksta_ter, &
37     ok_kzmin
38     use clesphys2, only: iflag_con, ok_orolf, ok_orodr, nbapp_rad, &
39     cycle_diurne, new_oliq, soil_model
40     use iniprint, only: prt_level
41 guez 3 use abort_gcm_m, only: abort_gcm
42     use YOMCST, only: rcpd, rtt, rlvtt, rg, ra, rsigma, retv, romega
43     use comgeomphy
44     use ctherm
45     use phytrac_m, only: phytrac
46     use oasis_m
47     use radepsi
48     use radopt
49     use yoethf
50     use ini_hist, only: ini_histhf, ini_histday, ini_histins
51     use orbite_m, only: orbite, zenang
52     use phyetat0_m, only: phyetat0, rlat, rlon
53     use hgardfou_m, only: hgardfou
54     use conf_phys_m, only: conf_phys
55 guez 15 use phyredem_m, only: phyredem
56 guez 3
57     ! Declaration des constantes et des fonctions thermodynamiques :
58     use fcttre, only: thermcep, foeew, qsats, qsatl
59    
60     ! Variables argument:
61    
62     INTEGER nq ! input nombre de traceurs (y compris vapeur d'eau)
63 guez 15
64     REAL, intent(in):: rdayvrai
65     ! (elapsed time since January 1st 0h of the starting year, in days)
66    
67 guez 3 REAL, intent(in):: gmtime ! heure de la journée en fraction de jour
68 guez 12 REAL, intent(in):: pdtphys ! pas d'integration pour la physique (seconde)
69 guez 7 LOGICAL, intent(in):: firstcal ! first call to "calfis"
70 guez 3 logical, intent(in):: lafin ! dernier passage
71    
72     REAL, intent(in):: paprs(klon, llm+1)
73     ! (pression pour chaque inter-couche, en Pa)
74 guez 15
75 guez 10 REAL, intent(in):: pplay(klon, llm)
76 guez 3 ! (input pression pour le mileu de chaque couche (en Pa))
77    
78     REAL pphi(klon, llm)
79     ! (input geopotentiel de chaque couche (g z) (reference sol))
80    
81     REAL pphis(klon) ! input geopotentiel du sol
82    
83     REAL presnivs(llm)
84     ! (input pressions approximat. des milieux couches ( en PA))
85    
86     REAL u(klon, llm) ! input vitesse dans la direction X (de O a E) en m/s
87     REAL v(klon, llm) ! input vitesse Y (de S a N) en m/s
88     REAL t(klon, llm) ! input temperature (K)
89    
90 guez 12 REAL, intent(in):: qx(klon, llm, nq)
91     ! (humidite specifique (kg/kg) et fractions massiques des autres traceurs)
92 guez 3
93     REAL omega(klon, llm) ! input vitesse verticale en Pa/s
94     REAL d_u(klon, llm) ! output tendance physique de "u" (m/s/s)
95     REAL d_v(klon, llm) ! output tendance physique de "v" (m/s/s)
96     REAL d_t(klon, llm) ! output tendance physique de "t" (K/s)
97     REAL d_qx(klon, llm, nq) ! output tendance physique de "qx" (kg/kg/s)
98     REAL d_ps(klon) ! output tendance physique de la pression au sol
99    
100     INTEGER nbteta
101     PARAMETER(nbteta=3)
102    
103     REAL PVteta(klon, nbteta)
104     ! (output vorticite potentielle a des thetas constantes)
105    
106     LOGICAL ok_cvl ! pour activer le nouveau driver pour convection KE
107     PARAMETER (ok_cvl=.TRUE.)
108     LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface
109     PARAMETER (ok_gust=.FALSE.)
110    
111     LOGICAL check ! Verifier la conservation du modele en eau
112     PARAMETER (check=.FALSE.)
113     LOGICAL ok_stratus ! Ajouter artificiellement les stratus
114     PARAMETER (ok_stratus=.FALSE.)
115    
116     ! Parametres lies au coupleur OASIS:
117     INTEGER, SAVE :: npas, nexca
118     logical rnpb
119     parameter(rnpb=.true.)
120    
121 guez 12 character(len=6), save:: ocean
122     ! (type de modèle océan à utiliser: "force" ou "slab" mais pas "couple")
123    
124 guez 3 logical ok_ocean
125     SAVE ok_ocean
126    
127     !IM "slab" ocean
128     REAL tslab(klon) !Temperature du slab-ocean
129     SAVE tslab
130     REAL seaice(klon) !glace de mer (kg/m2)
131     SAVE seaice
132     REAL fluxo(klon) !flux turbulents ocean-glace de mer
133     REAL fluxg(klon) !flux turbulents ocean-atmosphere
134    
135     ! Modele thermique du sol, a activer pour le cycle diurne:
136     logical ok_veget
137     save ok_veget
138     LOGICAL ok_journe ! sortir le fichier journalier
139     save ok_journe
140    
141     LOGICAL ok_mensuel ! sortir le fichier mensuel
142    
143     LOGICAL ok_instan ! sortir le fichier instantane
144     save ok_instan
145    
146     LOGICAL ok_region ! sortir le fichier regional
147     PARAMETER (ok_region=.FALSE.)
148    
149     ! pour phsystoke avec thermiques
150     REAL fm_therm(klon, llm+1)
151     REAL entr_therm(klon, llm)
152     real q2(klon, llm+1, nbsrf)
153     save q2
154    
155     INTEGER ivap ! indice de traceurs pour vapeur d'eau
156     PARAMETER (ivap=1)
157     INTEGER iliq ! indice de traceurs pour eau liquide
158     PARAMETER (iliq=2)
159    
160     REAL t_ancien(klon, llm), q_ancien(klon, llm)
161     SAVE t_ancien, q_ancien
162     LOGICAL ancien_ok
163     SAVE ancien_ok
164    
165     REAL d_t_dyn(klon, llm) ! tendance dynamique pour "t" (K/s)
166     REAL d_q_dyn(klon, llm) ! tendance dynamique pour "q" (kg/kg/s)
167    
168     real da(klon, llm), phi(klon, llm, llm), mp(klon, llm)
169    
170     !IM Amip2 PV a theta constante
171    
172     CHARACTER(LEN=3) ctetaSTD(nbteta)
173     DATA ctetaSTD/'350', '380', '405'/
174     REAL rtetaSTD(nbteta)
175     DATA rtetaSTD/350., 380., 405./
176    
177     !MI Amip2 PV a theta constante
178    
179     INTEGER klevp1
180     PARAMETER(klevp1=llm+1)
181    
182     REAL swdn0(klon, klevp1), swdn(klon, klevp1)
183     REAL swup0(klon, klevp1), swup(klon, klevp1)
184     SAVE swdn0, swdn, swup0, swup
185    
186     REAL SWdn200clr(klon), SWdn200(klon)
187     REAL SWup200clr(klon), SWup200(klon)
188     SAVE SWdn200clr, SWdn200, SWup200clr, SWup200
189    
190     REAL lwdn0(klon, klevp1), lwdn(klon, klevp1)
191     REAL lwup0(klon, klevp1), lwup(klon, klevp1)
192     SAVE lwdn0, lwdn, lwup0, lwup
193    
194     REAL LWdn200clr(klon), LWdn200(klon)
195     REAL LWup200clr(klon), LWup200(klon)
196     SAVE LWdn200clr, LWdn200, LWup200clr, LWup200
197    
198     !IM Amip2
199     ! variables a une pression donnee
200    
201     integer nlevSTD
202     PARAMETER(nlevSTD=17)
203     real rlevSTD(nlevSTD)
204     DATA rlevSTD/100000., 92500., 85000., 70000., &
205     60000., 50000., 40000., 30000., 25000., 20000., &
206     15000., 10000., 7000., 5000., 3000., 2000., 1000./
207     CHARACTER(LEN=4) clevSTD(nlevSTD)
208     DATA clevSTD/'1000', '925 ', '850 ', '700 ', '600 ', &
209     '500 ', '400 ', '300 ', '250 ', '200 ', '150 ', '100 ', &
210     '70 ', '50 ', '30 ', '20 ', '10 '/
211    
212     real tlevSTD(klon, nlevSTD), qlevSTD(klon, nlevSTD)
213     real rhlevSTD(klon, nlevSTD), philevSTD(klon, nlevSTD)
214     real ulevSTD(klon, nlevSTD), vlevSTD(klon, nlevSTD)
215     real wlevSTD(klon, nlevSTD)
216    
217     ! nout : niveau de output des variables a une pression donnee
218     INTEGER nout
219     PARAMETER(nout=3) !nout=1 : day; =2 : mth; =3 : NMC
220    
221     REAL tsumSTD(klon, nlevSTD, nout)
222     REAL usumSTD(klon, nlevSTD, nout), vsumSTD(klon, nlevSTD, nout)
223     REAL wsumSTD(klon, nlevSTD, nout), phisumSTD(klon, nlevSTD, nout)
224     REAL qsumSTD(klon, nlevSTD, nout), rhsumSTD(klon, nlevSTD, nout)
225    
226     SAVE tsumSTD, usumSTD, vsumSTD, wsumSTD, phisumSTD, &
227     qsumSTD, rhsumSTD
228    
229     logical oknondef(klon, nlevSTD, nout)
230     real tnondef(klon, nlevSTD, nout)
231     save tnondef
232    
233     ! les produits uvSTD, vqSTD, .., T2STD sont calcules
234     ! a partir des valeurs instantannees toutes les 6 h
235     ! qui sont moyennees sur le mois
236    
237     real uvSTD(klon, nlevSTD)
238     real vqSTD(klon, nlevSTD)
239     real vTSTD(klon, nlevSTD)
240     real wqSTD(klon, nlevSTD)
241    
242     real uvsumSTD(klon, nlevSTD, nout)
243     real vqsumSTD(klon, nlevSTD, nout)
244     real vTsumSTD(klon, nlevSTD, nout)
245     real wqsumSTD(klon, nlevSTD, nout)
246    
247     real vphiSTD(klon, nlevSTD)
248     real wTSTD(klon, nlevSTD)
249     real u2STD(klon, nlevSTD)
250     real v2STD(klon, nlevSTD)
251     real T2STD(klon, nlevSTD)
252    
253     real vphisumSTD(klon, nlevSTD, nout)
254     real wTsumSTD(klon, nlevSTD, nout)
255     real u2sumSTD(klon, nlevSTD, nout)
256     real v2sumSTD(klon, nlevSTD, nout)
257     real T2sumSTD(klon, nlevSTD, nout)
258    
259     SAVE uvsumSTD, vqsumSTD, vTsumSTD, wqsumSTD
260     SAVE vphisumSTD, wTsumSTD, u2sumSTD, v2sumSTD, T2sumSTD
261     !MI Amip2
262    
263     ! prw: precipitable water
264     real prw(klon)
265    
266     ! flwp, fiwp = Liquid Water Path & Ice Water Path (kg/m2)
267     ! flwc, fiwc = Liquid Water Content & Ice Water Content (kg/kg)
268     REAL flwp(klon), fiwp(klon)
269     REAL flwc(klon, llm), fiwc(klon, llm)
270    
271     INTEGER l, kmax, lmax
272     PARAMETER(kmax=8, lmax=8)
273     INTEGER kmaxm1, lmaxm1
274     PARAMETER(kmaxm1=kmax-1, lmaxm1=lmax-1)
275    
276     REAL zx_tau(kmaxm1), zx_pc(lmaxm1)
277     DATA zx_tau/0.0, 0.3, 1.3, 3.6, 9.4, 23., 60./
278     DATA zx_pc/50., 180., 310., 440., 560., 680., 800./
279    
280     ! cldtopres pression au sommet des nuages
281     REAL cldtopres(lmaxm1)
282     DATA cldtopres/50., 180., 310., 440., 560., 680., 800./
283    
284     ! taulev: numero du niveau de tau dans les sorties ISCCP
285     CHARACTER(LEN=4) taulev(kmaxm1)
286    
287     DATA taulev/'tau0', 'tau1', 'tau2', 'tau3', 'tau4', 'tau5', 'tau6'/
288     CHARACTER(LEN=3) pclev(lmaxm1)
289     DATA pclev/'pc1', 'pc2', 'pc3', 'pc4', 'pc5', 'pc6', 'pc7'/
290    
291     CHARACTER(LEN=28) cnameisccp(lmaxm1, kmaxm1)
292     DATA cnameisccp/'pc< 50hPa, tau< 0.3', 'pc= 50-180hPa, tau< 0.3', &
293     'pc= 180-310hPa, tau< 0.3', 'pc= 310-440hPa, tau< 0.3', &
294     'pc= 440-560hPa, tau< 0.3', 'pc= 560-680hPa, tau< 0.3', &
295     'pc= 680-800hPa, tau< 0.3', 'pc< 50hPa, tau= 0.3-1.3', &
296     'pc= 50-180hPa, tau= 0.3-1.3', 'pc= 180-310hPa, tau= 0.3-1.3', &
297     'pc= 310-440hPa, tau= 0.3-1.3', 'pc= 440-560hPa, tau= 0.3-1.3', &
298     'pc= 560-680hPa, tau= 0.3-1.3', 'pc= 680-800hPa, tau= 0.3-1.3', &
299     'pc< 50hPa, tau= 1.3-3.6', 'pc= 50-180hPa, tau= 1.3-3.6', &
300     'pc= 180-310hPa, tau= 1.3-3.6', 'pc= 310-440hPa, tau= 1.3-3.6', &
301     'pc= 440-560hPa, tau= 1.3-3.6', 'pc= 560-680hPa, tau= 1.3-3.6', &
302     'pc= 680-800hPa, tau= 1.3-3.6', 'pc< 50hPa, tau= 3.6-9.4', &
303     'pc= 50-180hPa, tau= 3.6-9.4', 'pc= 180-310hPa, tau= 3.6-9.4', &
304     'pc= 310-440hPa, tau= 3.6-9.4', 'pc= 440-560hPa, tau= 3.6-9.4', &
305     'pc= 560-680hPa, tau= 3.6-9.4', 'pc= 680-800hPa, tau= 3.6-9.4', &
306     'pc< 50hPa, tau= 9.4-23', 'pc= 50-180hPa, tau= 9.4-23', &
307     'pc= 180-310hPa, tau= 9.4-23', 'pc= 310-440hPa, tau= 9.4-23', &
308     'pc= 440-560hPa, tau= 9.4-23', 'pc= 560-680hPa, tau= 9.4-23', &
309     'pc= 680-800hPa, tau= 9.4-23', 'pc< 50hPa, tau= 23-60', &
310     'pc= 50-180hPa, tau= 23-60', 'pc= 180-310hPa, tau= 23-60', &
311     'pc= 310-440hPa, tau= 23-60', 'pc= 440-560hPa, tau= 23-60', &
312     'pc= 560-680hPa, tau= 23-60', 'pc= 680-800hPa, tau= 23-60', &
313     'pc< 50hPa, tau> 60.', 'pc= 50-180hPa, tau> 60.', &
314     'pc= 180-310hPa, tau> 60.', 'pc= 310-440hPa, tau> 60.', &
315     'pc= 440-560hPa, tau> 60.', 'pc= 560-680hPa, tau> 60.', &
316     'pc= 680-800hPa, tau> 60.'/
317    
318     !IM ISCCP simulator v3.4
319    
320     integer nid_hf, nid_hf3d
321     save nid_hf, nid_hf3d
322    
323     INTEGER longcles
324     PARAMETER ( longcles = 20 )
325    
326     ! Variables propres a la physique
327    
328     INTEGER, save:: radpas
329     ! (Radiative transfer computations are made every "radpas" call to
330     ! "physiq".)
331    
332     REAL radsol(klon)
333     SAVE radsol ! bilan radiatif au sol calcule par code radiatif
334    
335 guez 7 INTEGER, SAVE:: itap ! number of calls to "physiq"
336 guez 3
337     REAL ftsol(klon, nbsrf)
338     SAVE ftsol ! temperature du sol
339    
340     REAL ftsoil(klon, nsoilmx, nbsrf)
341     SAVE ftsoil ! temperature dans le sol
342    
343     REAL fevap(klon, nbsrf)
344     SAVE fevap ! evaporation
345     REAL fluxlat(klon, nbsrf)
346     SAVE fluxlat
347    
348     REAL fqsurf(klon, nbsrf)
349     SAVE fqsurf ! humidite de l'air au contact de la surface
350    
351     REAL qsol(klon)
352     SAVE qsol ! hauteur d'eau dans le sol
353    
354     REAL fsnow(klon, nbsrf)
355     SAVE fsnow ! epaisseur neigeuse
356    
357     REAL falbe(klon, nbsrf)
358     SAVE falbe ! albedo par type de surface
359     REAL falblw(klon, nbsrf)
360     SAVE falblw ! albedo par type de surface
361    
362 guez 13 ! Paramètres de l'orographie à l'échelle sous-maille (OESM) :
363     REAL, save:: zmea(klon) ! orographie moyenne
364     REAL, save:: zstd(klon) ! deviation standard de l'OESM
365     REAL, save:: zsig(klon) ! pente de l'OESM
366     REAL, save:: zgam(klon) ! anisotropie de l'OESM
367     REAL, save:: zthe(klon) ! orientation de l'OESM
368     REAL, save:: zpic(klon) ! Maximum de l'OESM
369     REAL, save:: zval(klon) ! Minimum de l'OESM
370     REAL, save:: rugoro(klon) ! longueur de rugosite de l'OESM
371 guez 3
372     REAL zulow(klon), zvlow(klon)
373    
374     INTEGER igwd, idx(klon), itest(klon)
375    
376     REAL agesno(klon, nbsrf)
377     SAVE agesno ! age de la neige
378    
379     REAL run_off_lic_0(klon)
380     SAVE run_off_lic_0
381     !KE43
382     ! Variables liees a la convection de K. Emanuel (sb):
383    
384     REAL bas, top ! cloud base and top levels
385     SAVE bas
386     SAVE top
387    
388     REAL Ma(klon, llm) ! undilute upward mass flux
389     SAVE Ma
390     REAL qcondc(klon, llm) ! in-cld water content from convect
391     SAVE qcondc
392     REAL ema_work1(klon, llm), ema_work2(klon, llm)
393     SAVE ema_work1, ema_work2
394    
395     REAL wd(klon) ! sb
396     SAVE wd ! sb
397    
398     ! Variables locales pour la couche limite (al1):
399    
400     ! Variables locales:
401    
402     REAL cdragh(klon) ! drag coefficient pour T and Q
403     REAL cdragm(klon) ! drag coefficient pour vent
404    
405     !AA Pour phytrac
406     REAL ycoefh(klon, llm) ! coef d'echange pour phytrac
407     REAL yu1(klon) ! vents dans la premiere couche U
408     REAL yv1(klon) ! vents dans la premiere couche V
409     REAL ffonte(klon, nbsrf) !Flux thermique utilise pour fondre la neige
410     REAL fqcalving(klon, nbsrf) !Flux d'eau "perdue" par la surface
411     ! !et necessaire pour limiter la
412     ! !hauteur de neige, en kg/m2/s
413     REAL zxffonte(klon), zxfqcalving(klon)
414    
415     REAL pfrac_impa(klon, llm)! Produits des coefs lessivage impaction
416     save pfrac_impa
417     REAL pfrac_nucl(klon, llm)! Produits des coefs lessivage nucleation
418     save pfrac_nucl
419     REAL pfrac_1nucl(klon, llm)! Produits des coefs lessi nucl (alpha = 1)
420     save pfrac_1nucl
421     REAL frac_impa(klon, llm) ! fractions d'aerosols lessivees (impaction)
422     REAL frac_nucl(klon, llm) ! idem (nucleation)
423    
424     !AA
425     REAL rain_fall(klon) ! pluie
426     REAL snow_fall(klon) ! neige
427     save snow_fall, rain_fall
428     !IM cf FH pour Tiedtke 080604
429     REAL rain_tiedtke(klon), snow_tiedtke(klon)
430    
431     REAL evap(klon), devap(klon) ! evaporation et sa derivee
432     REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
433     REAL dlw(klon) ! derivee infra rouge
434     SAVE dlw
435     REAL bils(klon) ! bilan de chaleur au sol
436     REAL fder(klon) ! Derive de flux (sensible et latente)
437     save fder
438     REAL ve(klon) ! integr. verticale du transport meri. de l'energie
439     REAL vq(klon) ! integr. verticale du transport meri. de l'eau
440     REAL ue(klon) ! integr. verticale du transport zonal de l'energie
441     REAL uq(klon) ! integr. verticale du transport zonal de l'eau
442    
443     REAL frugs(klon, nbsrf) ! longueur de rugosite
444     save frugs
445     REAL zxrugs(klon) ! longueur de rugosite
446    
447     ! Conditions aux limites
448    
449     INTEGER julien
450    
451 guez 7 INTEGER, SAVE:: lmt_pas ! number of time steps of "physics" per day
452 guez 3 REAL pctsrf(klon, nbsrf)
453     !IM
454     REAL pctsrf_new(klon, nbsrf) !pourcentage surfaces issus d'ORCHIDEE
455    
456     SAVE pctsrf ! sous-fraction du sol
457     REAL albsol(klon)
458     SAVE albsol ! albedo du sol total
459     REAL albsollw(klon)
460     SAVE albsollw ! albedo du sol total
461    
462     REAL, SAVE:: wo(klon, llm) ! ozone
463    
464     ! Declaration des procedures appelees
465    
466     EXTERNAL alboc ! calculer l'albedo sur ocean
467     EXTERNAL ajsec ! ajustement sec
468     EXTERNAL clmain ! couche limite
469     !KE43
470     EXTERNAL conema3 ! convect4.3
471     EXTERNAL fisrtilp ! schema de condensation a grande echelle (pluie)
472     EXTERNAL nuage ! calculer les proprietes radiatives
473     EXTERNAL ozonecm ! prescrire l'ozone
474     EXTERNAL radlwsw ! rayonnements solaire et infrarouge
475     EXTERNAL transp ! transport total de l'eau et de l'energie
476    
477     EXTERNAL ini_undefSTD !initialise a 0 une variable a 1 niveau de pression
478    
479 guez 7 EXTERNAL undefSTD
480     ! (somme les valeurs definies d'1 var a 1 niveau de pression)
481    
482 guez 3 ! Variables locales
483    
484     real clwcon(klon, llm), rnebcon(klon, llm)
485     real clwcon0(klon, llm), rnebcon0(klon, llm)
486    
487     save rnebcon, clwcon
488    
489     REAL rhcl(klon, llm) ! humiditi relative ciel clair
490     REAL dialiq(klon, llm) ! eau liquide nuageuse
491     REAL diafra(klon, llm) ! fraction nuageuse
492     REAL cldliq(klon, llm) ! eau liquide nuageuse
493     REAL cldfra(klon, llm) ! fraction nuageuse
494     REAL cldtau(klon, llm) ! epaisseur optique
495     REAL cldemi(klon, llm) ! emissivite infrarouge
496    
497     REAL fluxq(klon, llm, nbsrf) ! flux turbulent d'humidite
498     REAL fluxt(klon, llm, nbsrf) ! flux turbulent de chaleur
499     REAL fluxu(klon, llm, nbsrf) ! flux turbulent de vitesse u
500     REAL fluxv(klon, llm, nbsrf) ! flux turbulent de vitesse v
501    
502     REAL zxfluxt(klon, llm)
503     REAL zxfluxq(klon, llm)
504     REAL zxfluxu(klon, llm)
505     REAL zxfluxv(klon, llm)
506    
507     REAL heat(klon, llm) ! chauffage solaire
508     REAL heat0(klon, llm) ! chauffage solaire ciel clair
509     REAL cool(klon, llm) ! refroidissement infrarouge
510     REAL cool0(klon, llm) ! refroidissement infrarouge ciel clair
511     REAL topsw(klon), toplw(klon), solsw(klon), sollw(klon)
512     real sollwdown(klon) ! downward LW flux at surface
513     REAL topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)
514     REAL albpla(klon)
515     REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous surface
516     REAL fsolsw(klon, nbsrf) ! flux solaire absorb. pour chaque sous surface
517     ! Le rayonnement n'est pas calcule tous les pas, il faut donc
518     ! sauvegarder les sorties du rayonnement
519     SAVE heat, cool, albpla, topsw, toplw, solsw, sollw, sollwdown
520     SAVE topsw0, toplw0, solsw0, sollw0, heat0, cool0
521    
522     INTEGER itaprad
523     SAVE itaprad
524    
525     REAL conv_q(klon, llm) ! convergence de l'humidite (kg/kg/s)
526     REAL conv_t(klon, llm) ! convergence de la temperature(K/s)
527    
528     REAL cldl(klon), cldm(klon), cldh(klon) !nuages bas, moyen et haut
529     REAL cldt(klon), cldq(klon) !nuage total, eau liquide integree
530    
531     REAL zxtsol(klon), zxqsurf(klon), zxsnow(klon), zxfluxlat(klon)
532    
533     REAL dist, rmu0(klon), fract(klon)
534     REAL zdtime ! pas de temps du rayonnement (s)
535     real zlongi
536    
537     REAL z_avant(klon), z_apres(klon), z_factor(klon)
538     LOGICAL zx_ajustq
539    
540     REAL za, zb
541     REAL zx_t, zx_qs, zdelta, zcor, zlvdcp, zlsdcp
542     real zqsat(klon, llm)
543     INTEGER i, k, iq, nsrf
544     REAL t_coup
545     PARAMETER (t_coup=234.0)
546    
547     REAL zphi(klon, llm)
548    
549     !IM cf. AM Variables locales pour la CLA (hbtm2)
550    
551     REAL pblh(klon, nbsrf) ! Hauteur de couche limite
552     REAL plcl(klon, nbsrf) ! Niveau de condensation de la CLA
553     REAL capCL(klon, nbsrf) ! CAPE de couche limite
554     REAL oliqCL(klon, nbsrf) ! eau_liqu integree de couche limite
555     REAL cteiCL(klon, nbsrf) ! cloud top instab. crit. couche limite
556     REAL pblt(klon, nbsrf) ! T a la Hauteur de couche limite
557     REAL therm(klon, nbsrf)
558     REAL trmb1(klon, nbsrf) ! deep_cape
559     REAL trmb2(klon, nbsrf) ! inhibition
560     REAL trmb3(klon, nbsrf) ! Point Omega
561     ! Grdeurs de sorties
562     REAL s_pblh(klon), s_lcl(klon), s_capCL(klon)
563     REAL s_oliqCL(klon), s_cteiCL(klon), s_pblt(klon)
564     REAL s_therm(klon), s_trmb1(klon), s_trmb2(klon)
565     REAL s_trmb3(klon)
566    
567     ! Variables locales pour la convection de K. Emanuel (sb):
568    
569     REAL upwd(klon, llm) ! saturated updraft mass flux
570     REAL dnwd(klon, llm) ! saturated downdraft mass flux
571     REAL dnwd0(klon, llm) ! unsaturated downdraft mass flux
572     REAL tvp(klon, llm) ! virtual temp of lifted parcel
573     REAL cape(klon) ! CAPE
574     SAVE cape
575    
576     REAL pbase(klon) ! cloud base pressure
577     SAVE pbase
578     REAL bbase(klon) ! cloud base buoyancy
579     SAVE bbase
580     REAL rflag(klon) ! flag fonctionnement de convect
581     INTEGER iflagctrl(klon) ! flag fonctionnement de convect
582     ! -- convect43:
583     INTEGER ntra ! nb traceurs pour convect4.3
584     REAL dtvpdt1(klon, llm), dtvpdq1(klon, llm)
585     REAL dplcldt(klon), dplcldr(klon)
586    
587     ! Variables du changement
588    
589     ! con: convection
590     ! lsc: condensation a grande echelle (Large-Scale-Condensation)
591     ! ajs: ajustement sec
592     ! eva: evaporation de l'eau liquide nuageuse
593     ! vdf: couche limite (Vertical DiFfusion)
594     REAL d_t_con(klon, llm), d_q_con(klon, llm)
595     REAL d_u_con(klon, llm), d_v_con(klon, llm)
596     REAL d_t_lsc(klon, llm), d_q_lsc(klon, llm), d_ql_lsc(klon, llm)
597     REAL d_t_ajs(klon, llm), d_q_ajs(klon, llm)
598     REAL d_u_ajs(klon, llm), d_v_ajs(klon, llm)
599     REAL rneb(klon, llm)
600    
601     REAL pmfu(klon, llm), pmfd(klon, llm)
602     REAL pen_u(klon, llm), pen_d(klon, llm)
603     REAL pde_u(klon, llm), pde_d(klon, llm)
604     INTEGER kcbot(klon), kctop(klon), kdtop(klon)
605     REAL pmflxr(klon, llm+1), pmflxs(klon, llm+1)
606     REAL prfl(klon, llm+1), psfl(klon, llm+1)
607    
608     INTEGER ibas_con(klon), itop_con(klon)
609    
610     SAVE ibas_con, itop_con
611    
612     REAL rain_con(klon), rain_lsc(klon)
613     REAL snow_con(klon), snow_lsc(klon)
614     REAL d_ts(klon, nbsrf)
615    
616     REAL d_u_vdf(klon, llm), d_v_vdf(klon, llm)
617     REAL d_t_vdf(klon, llm), d_q_vdf(klon, llm)
618    
619     REAL d_u_oro(klon, llm), d_v_oro(klon, llm)
620     REAL d_t_oro(klon, llm)
621     REAL d_u_lif(klon, llm), d_v_lif(klon, llm)
622     REAL d_t_lif(klon, llm)
623    
624     REAL ratqs(klon, llm), ratqss(klon, llm), ratqsc(klon, llm)
625     real ratqsbas, ratqshaut
626     save ratqsbas, ratqshaut, ratqs
627    
628     ! Parametres lies au nouveau schema de nuages (SB, PDF)
629     real fact_cldcon
630     real facttemps
631     logical ok_newmicro
632     save ok_newmicro
633     save fact_cldcon, facttemps
634     real facteur
635    
636     integer iflag_cldcon
637     save iflag_cldcon
638    
639     logical ptconv(klon, llm)
640    
641     ! Variables liees a l'ecriture de la bande histoire physique
642    
643     integer itau_w ! pas de temps ecriture = itap + itau_phy
644    
645     ! Variables locales pour effectuer les appels en serie
646    
647     REAL t_seri(klon, llm), q_seri(klon, llm)
648     REAL ql_seri(klon, llm), qs_seri(klon, llm)
649     REAL u_seri(klon, llm), v_seri(klon, llm)
650    
651     REAL tr_seri(klon, llm, nbtr)
652     REAL d_tr(klon, llm, nbtr)
653    
654     REAL zx_rh(klon, llm)
655     INTEGER ndex2d(iim*(jjm + 1)), ndex3d(iim*(jjm + 1)*llm)
656    
657     REAL zustrdr(klon), zvstrdr(klon)
658     REAL zustrli(klon), zvstrli(klon)
659     REAL zustrph(klon), zvstrph(klon)
660     REAL aam, torsfc
661    
662     REAL dudyn(iim+1, jjm + 1, llm)
663    
664     REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique
665     REAL zx_tmp_fi3d(klon, llm) ! variable temporaire pour champs 3D
666    
667     REAL zx_tmp_2d(iim, jjm + 1), zx_tmp_3d(iim, jjm + 1, llm)
668    
669 guez 15 INTEGER, SAVE:: nid_day, nid_ins
670 guez 3
671     REAL ve_lay(klon, llm) ! transport meri. de l'energie a chaque niveau vert.
672     REAL vq_lay(klon, llm) ! transport meri. de l'eau a chaque niveau vert.
673     REAL ue_lay(klon, llm) ! transport zonal de l'energie a chaque niveau vert.
674     REAL uq_lay(klon, llm) ! transport zonal de l'eau a chaque niveau vert.
675    
676     REAL zsto
677    
678     character(len=20) modname
679     character(len=80) abort_message
680     logical ok_sync
681     real date0
682    
683     ! Variables liees au bilan d'energie et d'enthalpi
684     REAL ztsol(klon)
685     REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec
686     REAL d_h_vcol_phy
687     REAL fs_bound, fq_bound
688     SAVE d_h_vcol_phy
689     REAL zero_v(klon)
690     CHARACTER(LEN=15) ztit
691     INTEGER ip_ebil ! PRINT level for energy conserv. diag.
692     SAVE ip_ebil
693     DATA ip_ebil/0/
694     INTEGER if_ebil ! level for energy conserv. dignostics
695     SAVE if_ebil
696     !+jld ec_conser
697     REAL d_t_ec(klon, llm) ! tendance du a la conersion Ec -> E thermique
698     REAL ZRCPD
699     !-jld ec_conser
700     !IM: t2m, q2m, u10m, v10m
701     REAL t2m(klon, nbsrf), q2m(klon, nbsrf) !temperature, humidite a 2m
702     REAL u10m(klon, nbsrf), v10m(klon, nbsrf) !vents a 10m
703     REAL zt2m(klon), zq2m(klon) !temp., hum. 2m moyenne s/ 1 maille
704     REAL zu10m(klon), zv10m(klon) !vents a 10m moyennes s/1 maille
705     !jq Aerosol effects (Johannes Quaas, 27/11/2003)
706     REAL sulfate(klon, llm) ! SO4 aerosol concentration [ug/m3]
707    
708     REAL sulfate_pi(klon, llm)
709     ! (SO4 aerosol concentration [ug/m3] (pre-industrial value))
710     SAVE sulfate_pi
711    
712     REAL cldtaupi(klon, llm)
713     ! (Cloud optical thickness for pre-industrial (pi) aerosols)
714    
715     REAL re(klon, llm) ! Cloud droplet effective radius
716     REAL fl(klon, llm) ! denominator of re
717    
718     ! Aerosol optical properties
719     REAL tau_ae(klon, llm, 2), piz_ae(klon, llm, 2)
720     REAL cg_ae(klon, llm, 2)
721    
722     REAL topswad(klon), solswad(klon) ! Aerosol direct effect.
723     ! ok_ade=T -ADE=topswad-topsw
724    
725     REAL topswai(klon), solswai(klon) ! Aerosol indirect effect.
726     ! ok_aie=T ->
727     ! ok_ade=T -AIE=topswai-topswad
728     ! ok_ade=F -AIE=topswai-topsw
729    
730     REAL aerindex(klon) ! POLDER aerosol index
731    
732     ! Parameters
733     LOGICAL ok_ade, ok_aie ! Apply aerosol (in)direct effects or not
734     REAL bl95_b0, bl95_b1 ! Parameter in Boucher and Lohmann (1995)
735    
736     SAVE ok_ade, ok_aie, bl95_b0, bl95_b1
737     SAVE u10m
738     SAVE v10m
739     SAVE t2m
740     SAVE q2m
741     SAVE ffonte
742     SAVE fqcalving
743     SAVE piz_ae
744     SAVE tau_ae
745     SAVE cg_ae
746     SAVE rain_con
747     SAVE snow_con
748     SAVE topswai
749     SAVE topswad
750     SAVE solswai
751     SAVE solswad
752     SAVE d_u_con
753     SAVE d_v_con
754     SAVE rnebcon0
755     SAVE clwcon0
756     SAVE pblh
757     SAVE plcl
758     SAVE capCL
759     SAVE oliqCL
760     SAVE cteiCL
761     SAVE pblt
762     SAVE therm
763     SAVE trmb1
764     SAVE trmb2
765     SAVE trmb3
766    
767     !----------------------------------------------------------------
768    
769     modname = 'physiq'
770     IF (if_ebil >= 1) THEN
771     DO i=1, klon
772     zero_v(i)=0.
773     END DO
774     END IF
775     ok_sync=.TRUE.
776 guez 7 IF (nq < 2) THEN
777 guez 3 abort_message = 'eaux vapeur et liquide sont indispensables'
778 guez 13 CALL abort_gcm(modname, abort_message, 1)
779 guez 3 ENDIF
780    
781 guez 7 test_firstcal: IF (firstcal) THEN
782 guez 3 ! initialiser
783 guez 13 u10m=0.
784     v10m=0.
785     t2m=0.
786     q2m=0.
787     ffonte=0.
788     fqcalving=0.
789 guez 3 piz_ae(:, :, :)=0.
790     tau_ae(:, :, :)=0.
791     cg_ae(:, :, :)=0.
792     rain_con(:)=0.
793     snow_con(:)=0.
794     bl95_b0=0.
795     bl95_b1=0.
796     topswai(:)=0.
797     topswad(:)=0.
798     solswai(:)=0.
799     solswad(:)=0.
800    
801 guez 13 d_u_con = 0.0
802     d_v_con = 0.0
803     rnebcon0 = 0.0
804     clwcon0 = 0.0
805     rnebcon = 0.0
806     clwcon = 0.0
807 guez 3
808 guez 13 pblh =0. ! Hauteur de couche limite
809     plcl =0. ! Niveau de condensation de la CLA
810     capCL =0. ! CAPE de couche limite
811     oliqCL =0. ! eau_liqu integree de couche limite
812     cteiCL =0. ! cloud top instab. crit. couche limite
813     pblt =0. ! T a la Hauteur de couche limite
814     therm =0.
815     trmb1 =0. ! deep_cape
816     trmb2 =0. ! inhibition
817     trmb3 =0. ! Point Omega
818 guez 3
819     IF (if_ebil >= 1) d_h_vcol_phy=0.
820    
821     ! appel a la lecture du run.def physique
822    
823     call conf_phys(ocean, ok_veget, ok_journe, ok_mensuel, &
824     ok_instan, fact_cldcon, facttemps, ok_newmicro, &
825     iflag_cldcon, ratqsbas, ratqshaut, if_ebil, &
826     ok_ade, ok_aie, &
827     bl95_b0, bl95_b1, &
828     iflag_thermals, nsplit_thermals)
829    
830     ! Initialiser les compteurs:
831    
832     frugs = 0.
833     itap = 0
834     itaprad = 0
835 guez 12 CALL phyetat0("startphy.nc", pctsrf, ftsol, ftsoil, ocean, tslab, &
836     seaice, fqsurf, qsol, fsnow, &
837 guez 3 falbe, falblw, fevap, rain_fall, snow_fall, solsw, sollwdown, &
838 guez 13 dlw, radsol, frugs, agesno, &
839     zmea, zstd, zsig, zgam, zthe, zpic, zval, &
840 guez 3 t_ancien, q_ancien, ancien_ok, rnebcon, ratqs, clwcon, &
841     run_off_lic_0)
842    
843     ! ATTENTION : il faudra a terme relire q2 dans l'etat initial
844     q2(:, :, :)=1.e-8
845    
846 guez 12 radpas = NINT( 86400. / pdtphys / nbapp_rad)
847 guez 3
848     ! on remet le calendrier a zero
849 guez 15 IF (raz_date) itau_phy = 0
850 guez 3
851 guez 12 PRINT *, 'cycle_diurne = ', cycle_diurne
852 guez 3
853     IF(ocean.NE.'force ') THEN
854     ok_ocean=.TRUE.
855     ENDIF
856    
857 guez 12 CALL printflag(radpas, ok_ocean, ok_oasis, ok_journe, ok_instan, &
858     ok_region)
859 guez 3
860 guez 12 IF (pdtphys*REAL(radpas).GT.21600..AND.cycle_diurne) THEN
861     print *,'Nbre d appels au rayonnement insuffisant'
862     print *,"Au minimum 4 appels par jour si cycle diurne"
863 guez 3 abort_message='Nbre d appels au rayonnement insuffisant'
864     call abort_gcm(modname, abort_message, 1)
865     ENDIF
866 guez 12 print *,"Clef pour la convection, iflag_con=", iflag_con
867     print *,"Clef pour le driver de la convection, ok_cvl=", &
868 guez 3 ok_cvl
869    
870     ! Initialisation pour la convection de K.E. (sb):
871     IF (iflag_con >= 3) THEN
872    
873 guez 12 print *,"*** Convection de Kerry Emanuel 4.3 "
874 guez 3
875     !IM15/11/02 rajout initialisation ibas_con, itop_con cf. SB =>BEG
876     DO i = 1, klon
877     ibas_con(i) = 1
878     itop_con(i) = 1
879     ENDDO
880     !IM15/11/02 rajout initialisation ibas_con, itop_con cf. SB =>END
881    
882     ENDIF
883    
884     IF (ok_orodr) THEN
885 guez 13 rugoro = MAX(1e-5, zstd * zsig / 2)
886 guez 3 CALL SUGWD(klon, llm, paprs, pplay)
887 guez 13 else
888     rugoro = 0.
889 guez 3 ENDIF
890    
891 guez 12 lmt_pas = NINT(86400. / pdtphys) ! tous les jours
892 guez 7 print *, 'Number of time steps of "physics" per day: ', lmt_pas
893 guez 3
894 guez 12 ecrit_ins = NINT(ecrit_ins/pdtphys)
895     ecrit_hf = NINT(ecrit_hf/pdtphys)
896     ecrit_mth = NINT(ecrit_mth/pdtphys)
897     ecrit_tra = NINT(86400.*ecrit_tra/pdtphys)
898     ecrit_reg = NINT(ecrit_reg/pdtphys)
899 guez 3
900     ! Initialiser le couplage si necessaire
901    
902     npas = 0
903     nexca = 0
904    
905 guez 12 print *,'AVANT HIST IFLAG_CON=', iflag_con
906 guez 3
907     ! Initialisation des sorties
908    
909 guez 12 call ini_histhf(pdtphys, presnivs, nid_hf, nid_hf3d)
910     call ini_histday(pdtphys, presnivs, ok_journe, nid_day)
911     call ini_histins(pdtphys, presnivs, ok_instan, nid_ins)
912 guez 3 CALL ymds2ju(annee_ref, 1, int(day_ref), 0., date0)
913     !XXXPB Positionner date0 pour initialisation de ORCHIDEE
914     WRITE(*, *) 'physiq date0 : ', date0
915 guez 7 ENDIF test_firstcal
916 guez 3
917     ! Mettre a zero des variables de sortie (pour securite)
918    
919     DO i = 1, klon
920     d_ps(i) = 0.0
921     ENDDO
922     DO k = 1, llm
923     DO i = 1, klon
924     d_t(i, k) = 0.0
925     d_u(i, k) = 0.0
926     d_v(i, k) = 0.0
927     ENDDO
928     ENDDO
929     DO iq = 1, nq
930     DO k = 1, llm
931     DO i = 1, klon
932     d_qx(i, k, iq) = 0.0
933     ENDDO
934     ENDDO
935     ENDDO
936 guez 13 da=0.
937     mp=0.
938 guez 3 phi(:, :, :)=0.
939    
940     ! Ne pas affecter les valeurs entrees de u, v, h, et q
941    
942     DO k = 1, llm
943     DO i = 1, klon
944     t_seri(i, k) = t(i, k)
945     u_seri(i, k) = u(i, k)
946     v_seri(i, k) = v(i, k)
947     q_seri(i, k) = qx(i, k, ivap)
948     ql_seri(i, k) = qx(i, k, iliq)
949     qs_seri(i, k) = 0.
950     ENDDO
951     ENDDO
952     IF (nq >= 3) THEN
953     tr_seri(:, :, :nq-2) = qx(:, :, 3:nq)
954     ELSE
955     tr_seri(:, :, 1) = 0.
956     ENDIF
957    
958     DO i = 1, klon
959     ztsol(i) = 0.
960     ENDDO
961     DO nsrf = 1, nbsrf
962     DO i = 1, klon
963     ztsol(i) = ztsol(i) + ftsol(i, nsrf)*pctsrf(i, nsrf)
964     ENDDO
965     ENDDO
966    
967     IF (if_ebil >= 1) THEN
968     ztit='after dynamic'
969 guez 12 CALL diagetpq(airephy, ztit, ip_ebil, 1, 1, pdtphys &
970 guez 10 , t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs &
971 guez 3 , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
972     ! Comme les tendances de la physique sont ajoute dans la dynamique,
973     ! on devrait avoir que la variation d'entalpie par la dynamique
974     ! est egale a la variation de la physique au pas de temps precedent.
975     ! Donc la somme de ces 2 variations devrait etre nulle.
976     call diagphy(airephy, ztit, ip_ebil &
977     , zero_v, zero_v, zero_v, zero_v, zero_v &
978     , zero_v, zero_v, zero_v, ztsol &
979     , d_h_vcol+d_h_vcol_phy, d_qt, 0. &
980     , fs_bound, fq_bound )
981     END IF
982    
983     ! Diagnostiquer la tendance dynamique
984    
985     IF (ancien_ok) THEN
986     DO k = 1, llm
987     DO i = 1, klon
988 guez 12 d_t_dyn(i, k) = (t_seri(i, k)-t_ancien(i, k))/pdtphys
989     d_q_dyn(i, k) = (q_seri(i, k)-q_ancien(i, k))/pdtphys
990 guez 3 ENDDO
991     ENDDO
992     ELSE
993     DO k = 1, llm
994     DO i = 1, klon
995     d_t_dyn(i, k) = 0.0
996     d_q_dyn(i, k) = 0.0
997     ENDDO
998     ENDDO
999     ancien_ok = .TRUE.
1000     ENDIF
1001    
1002     ! Ajouter le geopotentiel du sol:
1003    
1004     DO k = 1, llm
1005     DO i = 1, klon
1006     zphi(i, k) = pphi(i, k) + pphis(i)
1007     ENDDO
1008     ENDDO
1009    
1010     ! Verifier les temperatures
1011    
1012     CALL hgardfou(t_seri, ftsol)
1013    
1014     ! Incrementer le compteur de la physique
1015    
1016 guez 7 itap = itap + 1
1017     julien = MOD(NINT(rdayvrai), 360)
1018 guez 3 if (julien == 0) julien = 360
1019    
1020     ! Mettre en action les conditions aux limites (albedo, sst, etc.).
1021     ! Prescrire l'ozone et calculer l'albedo sur l'ocean.
1022    
1023     IF (MOD(itap - 1, lmt_pas) == 0) THEN
1024     CALL ozonecm(REAL(julien), rlat, paprs, wo)
1025     ENDIF
1026    
1027     ! Re-evaporer l'eau liquide nuageuse
1028    
1029     DO k = 1, llm ! re-evaporation de l'eau liquide nuageuse
1030     DO i = 1, klon
1031     zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i, k))
1032     zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i, k))
1033     zdelta = MAX(0., SIGN(1., RTT-t_seri(i, k)))
1034     zb = MAX(0.0, ql_seri(i, k))
1035     za = - MAX(0.0, ql_seri(i, k)) &
1036     * (zlvdcp*(1.-zdelta)+zlsdcp*zdelta)
1037     t_seri(i, k) = t_seri(i, k) + za
1038     q_seri(i, k) = q_seri(i, k) + zb
1039     ql_seri(i, k) = 0.0
1040     ENDDO
1041     ENDDO
1042    
1043     IF (if_ebil >= 2) THEN
1044     ztit='after reevap'
1045 guez 12 CALL diagetpq(airephy, ztit, ip_ebil, 2, 1, pdtphys &
1046 guez 10 , t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs &
1047 guez 3 , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1048     call diagphy(airephy, ztit, ip_ebil &
1049     , zero_v, zero_v, zero_v, zero_v, zero_v &
1050     , zero_v, zero_v, zero_v, ztsol &
1051     , d_h_vcol, d_qt, d_ec &
1052     , fs_bound, fq_bound )
1053    
1054     END IF
1055    
1056     ! Appeler la diffusion verticale (programme de couche limite)
1057    
1058     DO i = 1, klon
1059     zxrugs(i) = 0.0
1060     ENDDO
1061     DO nsrf = 1, nbsrf
1062     DO i = 1, klon
1063     frugs(i, nsrf) = MAX(frugs(i, nsrf), 0.000015)
1064     ENDDO
1065     ENDDO
1066     DO nsrf = 1, nbsrf
1067     DO i = 1, klon
1068     zxrugs(i) = zxrugs(i) + frugs(i, nsrf)*pctsrf(i, nsrf)
1069     ENDDO
1070     ENDDO
1071    
1072     ! calculs necessaires au calcul de l'albedo dans l'interface
1073    
1074     CALL orbite(REAL(julien), zlongi, dist)
1075     IF (cycle_diurne) THEN
1076 guez 12 zdtime = pdtphys * REAL(radpas)
1077 guez 3 CALL zenang(zlongi, gmtime, zdtime, rmu0, fract)
1078     ELSE
1079     rmu0 = -999.999
1080     ENDIF
1081    
1082     ! Calcul de l'abedo moyen par maille
1083     albsol(:)=0.
1084     albsollw(:)=0.
1085     DO nsrf = 1, nbsrf
1086     DO i = 1, klon
1087     albsol(i) = albsol(i) + falbe(i, nsrf) * pctsrf(i, nsrf)
1088     albsollw(i) = albsollw(i) + falblw(i, nsrf) * pctsrf(i, nsrf)
1089     ENDDO
1090     ENDDO
1091    
1092     ! Repartition sous maille des flux LW et SW
1093     ! Repartition du longwave par sous-surface linearisee
1094    
1095     DO nsrf = 1, nbsrf
1096     DO i = 1, klon
1097     fsollw(i, nsrf) = sollw(i) &
1098     + 4.0*RSIGMA*ztsol(i)**3 * (ztsol(i)-ftsol(i, nsrf))
1099     fsolsw(i, nsrf) = solsw(i)*(1.-falbe(i, nsrf))/(1.-albsol(i))
1100     ENDDO
1101     ENDDO
1102    
1103     fder = dlw
1104    
1105 guez 12 CALL clmain(pdtphys, itap, date0, pctsrf, pctsrf_new, &
1106 guez 3 t_seri, q_seri, u_seri, v_seri, &
1107     julien, rmu0, co2_ppm, &
1108     ok_veget, ocean, npas, nexca, ftsol, &
1109     soil_model, cdmmax, cdhmax, &
1110     ksta, ksta_ter, ok_kzmin, ftsoil, qsol, &
1111     paprs, pplay, fsnow, fqsurf, fevap, falbe, falblw, &
1112     fluxlat, rain_fall, snow_fall, &
1113     fsolsw, fsollw, sollwdown, fder, &
1114     rlon, rlat, cuphy, cvphy, frugs, &
1115 guez 7 firstcal, lafin, agesno, rugoro, &
1116 guez 3 d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf, d_ts, &
1117     fluxt, fluxq, fluxu, fluxv, cdragh, cdragm, &
1118     q2, dsens, devap, &
1119     ycoefh, yu1, yv1, t2m, q2m, u10m, v10m, &
1120     pblh, capCL, oliqCL, cteiCL, pblT, &
1121     therm, trmb1, trmb2, trmb3, plcl, &
1122     fqcalving, ffonte, run_off_lic_0, &
1123     fluxo, fluxg, tslab, seaice)
1124    
1125     !XXX Incrementation des flux
1126    
1127     zxfluxt=0.
1128     zxfluxq=0.
1129     zxfluxu=0.
1130     zxfluxv=0.
1131     DO nsrf = 1, nbsrf
1132     DO k = 1, llm
1133     DO i = 1, klon
1134     zxfluxt(i, k) = zxfluxt(i, k) + &
1135     fluxt(i, k, nsrf) * pctsrf( i, nsrf)
1136     zxfluxq(i, k) = zxfluxq(i, k) + &
1137     fluxq(i, k, nsrf) * pctsrf( i, nsrf)
1138     zxfluxu(i, k) = zxfluxu(i, k) + &
1139     fluxu(i, k, nsrf) * pctsrf( i, nsrf)
1140     zxfluxv(i, k) = zxfluxv(i, k) + &
1141     fluxv(i, k, nsrf) * pctsrf( i, nsrf)
1142     END DO
1143     END DO
1144     END DO
1145     DO i = 1, klon
1146     sens(i) = - zxfluxt(i, 1) ! flux de chaleur sensible au sol
1147     evap(i) = - zxfluxq(i, 1) ! flux d'evaporation au sol
1148     fder(i) = dlw(i) + dsens(i) + devap(i)
1149     ENDDO
1150    
1151     DO k = 1, llm
1152     DO i = 1, klon
1153     t_seri(i, k) = t_seri(i, k) + d_t_vdf(i, k)
1154     q_seri(i, k) = q_seri(i, k) + d_q_vdf(i, k)
1155     u_seri(i, k) = u_seri(i, k) + d_u_vdf(i, k)
1156     v_seri(i, k) = v_seri(i, k) + d_v_vdf(i, k)
1157     ENDDO
1158     ENDDO
1159    
1160     IF (if_ebil >= 2) THEN
1161     ztit='after clmain'
1162 guez 12 CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, pdtphys &
1163 guez 10 , t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs &
1164 guez 3 , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1165     call diagphy(airephy, ztit, ip_ebil &
1166     , zero_v, zero_v, zero_v, zero_v, sens &
1167     , evap, zero_v, zero_v, ztsol &
1168     , d_h_vcol, d_qt, d_ec &
1169     , fs_bound, fq_bound )
1170     END IF
1171    
1172     ! Incrementer la temperature du sol
1173    
1174     DO i = 1, klon
1175     zxtsol(i) = 0.0
1176     zxfluxlat(i) = 0.0
1177    
1178     zt2m(i) = 0.0
1179     zq2m(i) = 0.0
1180     zu10m(i) = 0.0
1181     zv10m(i) = 0.0
1182     zxffonte(i) = 0.0
1183     zxfqcalving(i) = 0.0
1184    
1185     s_pblh(i) = 0.0
1186     s_lcl(i) = 0.0
1187     s_capCL(i) = 0.0
1188     s_oliqCL(i) = 0.0
1189     s_cteiCL(i) = 0.0
1190     s_pblT(i) = 0.0
1191     s_therm(i) = 0.0
1192     s_trmb1(i) = 0.0
1193     s_trmb2(i) = 0.0
1194     s_trmb3(i) = 0.0
1195    
1196     IF ( abs( pctsrf(i, is_ter) + pctsrf(i, is_lic) + &
1197     pctsrf(i, is_oce) + pctsrf(i, is_sic) - 1.) .GT. EPSFRA) &
1198     THEN
1199     WRITE(*, *) 'physiq : pb sous surface au point ', i, &
1200     pctsrf(i, 1 : nbsrf)
1201     ENDIF
1202     ENDDO
1203     DO nsrf = 1, nbsrf
1204     DO i = 1, klon
1205     ftsol(i, nsrf) = ftsol(i, nsrf) + d_ts(i, nsrf)
1206     zxtsol(i) = zxtsol(i) + ftsol(i, nsrf)*pctsrf(i, nsrf)
1207     zxfluxlat(i) = zxfluxlat(i) + fluxlat(i, nsrf)*pctsrf(i, nsrf)
1208    
1209     zt2m(i) = zt2m(i) + t2m(i, nsrf)*pctsrf(i, nsrf)
1210     zq2m(i) = zq2m(i) + q2m(i, nsrf)*pctsrf(i, nsrf)
1211     zu10m(i) = zu10m(i) + u10m(i, nsrf)*pctsrf(i, nsrf)
1212     zv10m(i) = zv10m(i) + v10m(i, nsrf)*pctsrf(i, nsrf)
1213     zxffonte(i) = zxffonte(i) + ffonte(i, nsrf)*pctsrf(i, nsrf)
1214     zxfqcalving(i) = zxfqcalving(i) + &
1215     fqcalving(i, nsrf)*pctsrf(i, nsrf)
1216     s_pblh(i) = s_pblh(i) + pblh(i, nsrf)*pctsrf(i, nsrf)
1217     s_lcl(i) = s_lcl(i) + plcl(i, nsrf)*pctsrf(i, nsrf)
1218     s_capCL(i) = s_capCL(i) + capCL(i, nsrf) *pctsrf(i, nsrf)
1219     s_oliqCL(i) = s_oliqCL(i) + oliqCL(i, nsrf) *pctsrf(i, nsrf)
1220     s_cteiCL(i) = s_cteiCL(i) + cteiCL(i, nsrf) *pctsrf(i, nsrf)
1221     s_pblT(i) = s_pblT(i) + pblT(i, nsrf) *pctsrf(i, nsrf)
1222     s_therm(i) = s_therm(i) + therm(i, nsrf) *pctsrf(i, nsrf)
1223     s_trmb1(i) = s_trmb1(i) + trmb1(i, nsrf) *pctsrf(i, nsrf)
1224     s_trmb2(i) = s_trmb2(i) + trmb2(i, nsrf) *pctsrf(i, nsrf)
1225     s_trmb3(i) = s_trmb3(i) + trmb3(i, nsrf) *pctsrf(i, nsrf)
1226     ENDDO
1227     ENDDO
1228    
1229     ! Si une sous-fraction n'existe pas, elle prend la temp. moyenne
1230    
1231     DO nsrf = 1, nbsrf
1232     DO i = 1, klon
1233 guez 7 IF (pctsrf(i, nsrf) < epsfra) ftsol(i, nsrf) = zxtsol(i)
1234 guez 3
1235 guez 7 IF (pctsrf(i, nsrf) < epsfra) t2m(i, nsrf) = zt2m(i)
1236     IF (pctsrf(i, nsrf) < epsfra) q2m(i, nsrf) = zq2m(i)
1237     IF (pctsrf(i, nsrf) < epsfra) u10m(i, nsrf) = zu10m(i)
1238     IF (pctsrf(i, nsrf) < epsfra) v10m(i, nsrf) = zv10m(i)
1239     IF (pctsrf(i, nsrf) < epsfra) ffonte(i, nsrf) = zxffonte(i)
1240     IF (pctsrf(i, nsrf) < epsfra) &
1241 guez 3 fqcalving(i, nsrf) = zxfqcalving(i)
1242 guez 7 IF (pctsrf(i, nsrf) < epsfra) pblh(i, nsrf)=s_pblh(i)
1243     IF (pctsrf(i, nsrf) < epsfra) plcl(i, nsrf)=s_lcl(i)
1244     IF (pctsrf(i, nsrf) < epsfra) capCL(i, nsrf)=s_capCL(i)
1245     IF (pctsrf(i, nsrf) < epsfra) oliqCL(i, nsrf)=s_oliqCL(i)
1246     IF (pctsrf(i, nsrf) < epsfra) cteiCL(i, nsrf)=s_cteiCL(i)
1247     IF (pctsrf(i, nsrf) < epsfra) pblT(i, nsrf)=s_pblT(i)
1248     IF (pctsrf(i, nsrf) < epsfra) therm(i, nsrf)=s_therm(i)
1249     IF (pctsrf(i, nsrf) < epsfra) trmb1(i, nsrf)=s_trmb1(i)
1250     IF (pctsrf(i, nsrf) < epsfra) trmb2(i, nsrf)=s_trmb2(i)
1251     IF (pctsrf(i, nsrf) < epsfra) trmb3(i, nsrf)=s_trmb3(i)
1252 guez 3 ENDDO
1253     ENDDO
1254    
1255     ! Calculer la derive du flux infrarouge
1256    
1257     DO i = 1, klon
1258     dlw(i) = - 4.0*RSIGMA*zxtsol(i)**3
1259     ENDDO
1260    
1261     ! Appeler la convection (au choix)
1262    
1263     DO k = 1, llm
1264     DO i = 1, klon
1265     conv_q(i, k) = d_q_dyn(i, k) &
1266 guez 12 + d_q_vdf(i, k)/pdtphys
1267 guez 3 conv_t(i, k) = d_t_dyn(i, k) &
1268 guez 12 + d_t_vdf(i, k)/pdtphys
1269 guez 3 ENDDO
1270     ENDDO
1271     IF (check) THEN
1272     za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)
1273 guez 12 print *, "avantcon=", za
1274 guez 3 ENDIF
1275     zx_ajustq = .FALSE.
1276     IF (iflag_con == 2) zx_ajustq=.TRUE.
1277     IF (zx_ajustq) THEN
1278     DO i = 1, klon
1279     z_avant(i) = 0.0
1280     ENDDO
1281     DO k = 1, llm
1282     DO i = 1, klon
1283     z_avant(i) = z_avant(i) + (q_seri(i, k)+ql_seri(i, k)) &
1284     *(paprs(i, k)-paprs(i, k+1))/RG
1285     ENDDO
1286     ENDDO
1287     ENDIF
1288     IF (iflag_con == 1) THEN
1289     stop 'reactiver le call conlmd dans physiq.F'
1290     ELSE IF (iflag_con == 2) THEN
1291 guez 12 CALL conflx(pdtphys, paprs, pplay, t_seri, q_seri, &
1292 guez 3 conv_t, conv_q, zxfluxq(1, 1), omega, &
1293     d_t_con, d_q_con, rain_con, snow_con, &
1294     pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
1295     kcbot, kctop, kdtop, pmflxr, pmflxs)
1296     WHERE (rain_con < 0.) rain_con = 0.
1297     WHERE (snow_con < 0.) snow_con = 0.
1298     DO i = 1, klon
1299     ibas_con(i) = llm+1 - kcbot(i)
1300     itop_con(i) = llm+1 - kctop(i)
1301     ENDDO
1302     ELSE IF (iflag_con >= 3) THEN
1303     ! nb of tracers for the KE convection:
1304     ! MAF la partie traceurs est faite dans phytrac
1305     ! on met ntra=1 pour limiter les appels mais on peut
1306     ! supprimer les calculs / ftra.
1307     ntra = 1
1308     ! Schema de convection modularise et vectorise:
1309     ! (driver commun aux versions 3 et 4)
1310    
1311     IF (ok_cvl) THEN ! new driver for convectL
1312     CALL concvl (iflag_con, &
1313 guez 12 pdtphys, paprs, pplay, t_seri, q_seri, &
1314 guez 3 u_seri, v_seri, tr_seri, ntra, &
1315     ema_work1, ema_work2, &
1316     d_t_con, d_q_con, d_u_con, d_v_con, d_tr, &
1317     rain_con, snow_con, ibas_con, itop_con, &
1318     upwd, dnwd, dnwd0, &
1319     Ma, cape, tvp, iflagctrl, &
1320     pbase, bbase, dtvpdt1, dtvpdq1, dplcldt, dplcldr, qcondc, wd, &
1321     pmflxr, pmflxs, &
1322     da, phi, mp)
1323    
1324     clwcon0=qcondc
1325 guez 13 pmfu=upwd+dnwd
1326 guez 3 ELSE ! ok_cvl
1327     ! MAF conema3 ne contient pas les traceurs
1328 guez 12 CALL conema3 (pdtphys, &
1329 guez 3 paprs, pplay, t_seri, q_seri, &
1330     u_seri, v_seri, tr_seri, ntra, &
1331     ema_work1, ema_work2, &
1332     d_t_con, d_q_con, d_u_con, d_v_con, d_tr, &
1333     rain_con, snow_con, ibas_con, itop_con, &
1334     upwd, dnwd, dnwd0, bas, top, &
1335     Ma, cape, tvp, rflag, &
1336     pbase &
1337     , bbase, dtvpdt1, dtvpdq1, dplcldt, dplcldr &
1338     , clwcon0)
1339     ENDIF ! ok_cvl
1340    
1341     IF (.NOT. ok_gust) THEN
1342     do i = 1, klon
1343     wd(i)=0.0
1344     enddo
1345     ENDIF
1346    
1347     ! Calcul des proprietes des nuages convectifs
1348    
1349     DO k = 1, llm
1350     DO i = 1, klon
1351     zx_t = t_seri(i, k)
1352     IF (thermcep) THEN
1353     zdelta = MAX(0., SIGN(1., rtt-zx_t))
1354     zx_qs = r2es * FOEEW(zx_t, zdelta)/pplay(i, k)
1355     zx_qs = MIN(0.5, zx_qs)
1356     zcor = 1./(1.-retv*zx_qs)
1357     zx_qs = zx_qs*zcor
1358     ELSE
1359 guez 7 IF (zx_t < t_coup) THEN
1360 guez 3 zx_qs = qsats(zx_t)/pplay(i, k)
1361     ELSE
1362     zx_qs = qsatl(zx_t)/pplay(i, k)
1363     ENDIF
1364     ENDIF
1365     zqsat(i, k)=zx_qs
1366     ENDDO
1367     ENDDO
1368    
1369     ! calcul des proprietes des nuages convectifs
1370 guez 13 clwcon0=fact_cldcon*clwcon0
1371 guez 3 call clouds_gno &
1372     (klon, llm, q_seri, zqsat, clwcon0, ptconv, ratqsc, rnebcon0)
1373     ELSE
1374 guez 12 print *, "iflag_con non-prevu", iflag_con
1375 guez 3 stop 1
1376     ENDIF
1377    
1378     DO k = 1, llm
1379     DO i = 1, klon
1380     t_seri(i, k) = t_seri(i, k) + d_t_con(i, k)
1381     q_seri(i, k) = q_seri(i, k) + d_q_con(i, k)
1382     u_seri(i, k) = u_seri(i, k) + d_u_con(i, k)
1383     v_seri(i, k) = v_seri(i, k) + d_v_con(i, k)
1384     ENDDO
1385     ENDDO
1386    
1387     IF (if_ebil >= 2) THEN
1388     ztit='after convect'
1389 guez 12 CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, pdtphys &
1390 guez 10 , t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs &
1391 guez 3 , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1392     call diagphy(airephy, ztit, ip_ebil &
1393     , zero_v, zero_v, zero_v, zero_v, zero_v &
1394     , zero_v, rain_con, snow_con, ztsol &
1395     , d_h_vcol, d_qt, d_ec &
1396     , fs_bound, fq_bound )
1397     END IF
1398    
1399     IF (check) THEN
1400     za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)
1401 guez 12 print *,"aprescon=", za
1402 guez 3 zx_t = 0.0
1403     za = 0.0
1404     DO i = 1, klon
1405     za = za + airephy(i)/REAL(klon)
1406     zx_t = zx_t + (rain_con(i)+ &
1407     snow_con(i))*airephy(i)/REAL(klon)
1408     ENDDO
1409 guez 12 zx_t = zx_t/za*pdtphys
1410     print *,"Precip=", zx_t
1411 guez 3 ENDIF
1412     IF (zx_ajustq) THEN
1413     DO i = 1, klon
1414     z_apres(i) = 0.0
1415     ENDDO
1416     DO k = 1, llm
1417     DO i = 1, klon
1418     z_apres(i) = z_apres(i) + (q_seri(i, k)+ql_seri(i, k)) &
1419     *(paprs(i, k)-paprs(i, k+1))/RG
1420     ENDDO
1421     ENDDO
1422     DO i = 1, klon
1423 guez 12 z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*pdtphys) &
1424 guez 3 /z_apres(i)
1425     ENDDO
1426     DO k = 1, llm
1427     DO i = 1, klon
1428     IF (z_factor(i).GT.(1.0+1.0E-08) .OR. &
1429 guez 7 z_factor(i) < (1.0-1.0E-08)) THEN
1430 guez 3 q_seri(i, k) = q_seri(i, k) * z_factor(i)
1431     ENDIF
1432     ENDDO
1433     ENDDO
1434     ENDIF
1435     zx_ajustq=.FALSE.
1436    
1437     ! Convection seche (thermiques ou ajustement)
1438    
1439 guez 13 d_t_ajs=0.
1440     d_u_ajs=0.
1441     d_v_ajs=0.
1442     d_q_ajs=0.
1443     fm_therm=0.
1444     entr_therm=0.
1445 guez 3
1446 guez 12 IF(prt_level>9)print *, &
1447 guez 3 'AVANT LA CONVECTION SECHE, iflag_thermals=' &
1448     , iflag_thermals, ' nsplit_thermals=', nsplit_thermals
1449 guez 7 if(iflag_thermals < 0) then
1450 guez 3 ! Rien
1451 guez 12 IF(prt_level>9)print *,'pas de convection'
1452 guez 3 else if(iflag_thermals == 0) then
1453     ! Ajustement sec
1454 guez 12 IF(prt_level>9)print *,'ajsec'
1455 guez 3 CALL ajsec(paprs, pplay, t_seri, q_seri, d_t_ajs, d_q_ajs)
1456 guez 13 t_seri = t_seri + d_t_ajs
1457     q_seri = q_seri + d_q_ajs
1458 guez 3 else
1459     ! Thermiques
1460 guez 12 IF(prt_level>9)print *,'JUSTE AVANT, iflag_thermals=' &
1461 guez 3 , iflag_thermals, ' nsplit_thermals=', nsplit_thermals
1462     call calltherm(pdtphys &
1463     , pplay, paprs, pphi &
1464     , u_seri, v_seri, t_seri, q_seri &
1465     , d_u_ajs, d_v_ajs, d_t_ajs, d_q_ajs &
1466     , fm_therm, entr_therm)
1467     endif
1468    
1469     IF (if_ebil >= 2) THEN
1470     ztit='after dry_adjust'
1471 guez 12 CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, pdtphys &
1472 guez 10 , t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs &
1473 guez 3 , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1474     END IF
1475    
1476     ! Caclul des ratqs
1477    
1478     ! ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
1479     ! on ecrase le tableau ratqsc calcule par clouds_gno
1480     if (iflag_cldcon == 1) then
1481     do k=1, llm
1482     do i=1, klon
1483     if(ptconv(i, k)) then
1484     ratqsc(i, k)=ratqsbas &
1485     +fact_cldcon*(q_seri(i, 1)-q_seri(i, k))/q_seri(i, k)
1486     else
1487     ratqsc(i, k)=0.
1488     endif
1489     enddo
1490     enddo
1491     endif
1492    
1493     ! ratqs stables
1494     do k=1, llm
1495     do i=1, klon
1496     ratqss(i, k)=ratqsbas+(ratqshaut-ratqsbas)* &
1497     min((paprs(i, 1)-pplay(i, k))/(paprs(i, 1)-30000.), 1.)
1498     enddo
1499     enddo
1500    
1501     ! ratqs final
1502     if (iflag_cldcon == 1 .or.iflag_cldcon == 2) then
1503     ! les ratqs sont une conbinaison de ratqss et ratqsc
1504     ! ratqs final
1505     ! 1e4 (en gros 3 heures), en dur pour le moment, est le temps de
1506     ! relaxation des ratqs
1507     facteur=exp(-pdtphys*facttemps)
1508 guez 13 ratqs=max(ratqs*facteur, ratqss)
1509     ratqs=max(ratqs, ratqsc)
1510 guez 3 else
1511     ! on ne prend que le ratqs stable pour fisrtilp
1512 guez 13 ratqs=ratqss
1513 guez 3 endif
1514    
1515     ! Appeler le processus de condensation a grande echelle
1516     ! et le processus de precipitation
1517 guez 12 CALL fisrtilp(pdtphys, paprs, pplay, &
1518 guez 3 t_seri, q_seri, ptconv, ratqs, &
1519     d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq, &
1520     rain_lsc, snow_lsc, &
1521     pfrac_impa, pfrac_nucl, pfrac_1nucl, &
1522     frac_impa, frac_nucl, &
1523     prfl, psfl, rhcl)
1524    
1525     WHERE (rain_lsc < 0) rain_lsc = 0.
1526     WHERE (snow_lsc < 0) snow_lsc = 0.
1527     DO k = 1, llm
1528     DO i = 1, klon
1529     t_seri(i, k) = t_seri(i, k) + d_t_lsc(i, k)
1530     q_seri(i, k) = q_seri(i, k) + d_q_lsc(i, k)
1531     ql_seri(i, k) = ql_seri(i, k) + d_ql_lsc(i, k)
1532     cldfra(i, k) = rneb(i, k)
1533     IF (.NOT.new_oliq) cldliq(i, k) = ql_seri(i, k)
1534     ENDDO
1535     ENDDO
1536     IF (check) THEN
1537     za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)
1538 guez 12 print *,"apresilp=", za
1539 guez 3 zx_t = 0.0
1540     za = 0.0
1541     DO i = 1, klon
1542     za = za + airephy(i)/REAL(klon)
1543     zx_t = zx_t + (rain_lsc(i) &
1544     + snow_lsc(i))*airephy(i)/REAL(klon)
1545     ENDDO
1546 guez 12 zx_t = zx_t/za*pdtphys
1547     print *,"Precip=", zx_t
1548 guez 3 ENDIF
1549    
1550     IF (if_ebil >= 2) THEN
1551     ztit='after fisrt'
1552 guez 12 CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, pdtphys &
1553 guez 10 , t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs &
1554 guez 3 , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1555     call diagphy(airephy, ztit, ip_ebil &
1556     , zero_v, zero_v, zero_v, zero_v, zero_v &
1557     , zero_v, rain_lsc, snow_lsc, ztsol &
1558     , d_h_vcol, d_qt, d_ec &
1559     , fs_bound, fq_bound )
1560     END IF
1561    
1562     ! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
1563    
1564     ! 1. NUAGES CONVECTIFS
1565    
1566     IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
1567     snow_tiedtke=0.
1568     if (iflag_cldcon == -1) then
1569     rain_tiedtke=rain_con
1570     else
1571     rain_tiedtke=0.
1572     do k=1, llm
1573     do i=1, klon
1574 guez 7 if (d_q_con(i, k) < 0.) then
1575 guez 3 rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i, k)/pdtphys &
1576     *(paprs(i, k)-paprs(i, k+1))/rg
1577     endif
1578     enddo
1579     enddo
1580     endif
1581    
1582     ! Nuages diagnostiques pour Tiedtke
1583     CALL diagcld1(paprs, pplay, &
1584     rain_tiedtke, snow_tiedtke, ibas_con, itop_con, &
1585     diafra, dialiq)
1586     DO k = 1, llm
1587     DO i = 1, klon
1588     IF (diafra(i, k).GT.cldfra(i, k)) THEN
1589     cldliq(i, k) = dialiq(i, k)
1590     cldfra(i, k) = diafra(i, k)
1591     ENDIF
1592     ENDDO
1593     ENDDO
1594    
1595     ELSE IF (iflag_cldcon == 3) THEN
1596 guez 7 ! On prend pour les nuages convectifs le max du calcul de la
1597     ! convection et du calcul du pas de temps précédent diminué d'un facteur
1598     ! facttemps
1599 guez 3 facteur = pdtphys *facttemps
1600     do k=1, llm
1601     do i=1, klon
1602     rnebcon(i, k)=rnebcon(i, k)*facteur
1603     if (rnebcon0(i, k)*clwcon0(i, k).gt.rnebcon(i, k)*clwcon(i, k)) &
1604     then
1605     rnebcon(i, k)=rnebcon0(i, k)
1606     clwcon(i, k)=clwcon0(i, k)
1607     endif
1608     enddo
1609     enddo
1610    
1611     ! On prend la somme des fractions nuageuses et des contenus en eau
1612 guez 13 cldfra=min(max(cldfra, rnebcon), 1.)
1613     cldliq=cldliq+rnebcon*clwcon
1614 guez 3
1615     ENDIF
1616    
1617     ! 2. NUAGES STARTIFORMES
1618    
1619     IF (ok_stratus) THEN
1620     CALL diagcld2(paprs, pplay, t_seri, q_seri, diafra, dialiq)
1621     DO k = 1, llm
1622     DO i = 1, klon
1623     IF (diafra(i, k).GT.cldfra(i, k)) THEN
1624     cldliq(i, k) = dialiq(i, k)
1625     cldfra(i, k) = diafra(i, k)
1626     ENDIF
1627     ENDDO
1628     ENDDO
1629     ENDIF
1630    
1631     ! Precipitation totale
1632    
1633     DO i = 1, klon
1634     rain_fall(i) = rain_con(i) + rain_lsc(i)
1635     snow_fall(i) = snow_con(i) + snow_lsc(i)
1636     ENDDO
1637    
1638     IF (if_ebil >= 2) THEN
1639     ztit="after diagcld"
1640 guez 12 CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, pdtphys &
1641 guez 10 , t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs &
1642 guez 3 , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1643     END IF
1644    
1645     ! Calculer l'humidite relative pour diagnostique
1646    
1647     DO k = 1, llm
1648     DO i = 1, klon
1649     zx_t = t_seri(i, k)
1650     IF (thermcep) THEN
1651     zdelta = MAX(0., SIGN(1., rtt-zx_t))
1652     zx_qs = r2es * FOEEW(zx_t, zdelta)/pplay(i, k)
1653     zx_qs = MIN(0.5, zx_qs)
1654     zcor = 1./(1.-retv*zx_qs)
1655     zx_qs = zx_qs*zcor
1656     ELSE
1657 guez 7 IF (zx_t < t_coup) THEN
1658 guez 3 zx_qs = qsats(zx_t)/pplay(i, k)
1659     ELSE
1660     zx_qs = qsatl(zx_t)/pplay(i, k)
1661     ENDIF
1662     ENDIF
1663     zx_rh(i, k) = q_seri(i, k)/zx_qs
1664     zqsat(i, k)=zx_qs
1665     ENDDO
1666     ENDDO
1667     !jq - introduce the aerosol direct and first indirect radiative forcings
1668     !jq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
1669     IF (ok_ade.OR.ok_aie) THEN
1670     ! Get sulfate aerosol distribution
1671 guez 7 CALL readsulfate(rdayvrai, firstcal, sulfate)
1672     CALL readsulfate_preind(rdayvrai, firstcal, sulfate_pi)
1673 guez 3
1674     ! Calculate aerosol optical properties (Olivier Boucher)
1675     CALL aeropt(pplay, paprs, t_seri, sulfate, rhcl, &
1676     tau_ae, piz_ae, cg_ae, aerindex)
1677     ELSE
1678     tau_ae(:, :, :)=0.0
1679     piz_ae(:, :, :)=0.0
1680     cg_ae(:, :, :)=0.0
1681     ENDIF
1682    
1683     ! Calculer les parametres optiques des nuages et quelques
1684     ! parametres pour diagnostiques:
1685    
1686     if (ok_newmicro) then
1687     CALL newmicro (paprs, pplay, ok_newmicro, &
1688     t_seri, cldliq, cldfra, cldtau, cldemi, &
1689     cldh, cldl, cldm, cldt, cldq, &
1690     flwp, fiwp, flwc, fiwc, &
1691     ok_aie, &
1692     sulfate, sulfate_pi, &
1693     bl95_b0, bl95_b1, &
1694     cldtaupi, re, fl)
1695     else
1696     CALL nuage (paprs, pplay, &
1697     t_seri, cldliq, cldfra, cldtau, cldemi, &
1698     cldh, cldl, cldm, cldt, cldq, &
1699     ok_aie, &
1700     sulfate, sulfate_pi, &
1701     bl95_b0, bl95_b1, &
1702     cldtaupi, re, fl)
1703    
1704     endif
1705    
1706     ! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
1707    
1708     IF (MOD(itaprad, radpas) == 0) THEN
1709     DO i = 1, klon
1710     albsol(i) = falbe(i, is_oce) * pctsrf(i, is_oce) &
1711     + falbe(i, is_lic) * pctsrf(i, is_lic) &
1712     + falbe(i, is_ter) * pctsrf(i, is_ter) &
1713     + falbe(i, is_sic) * pctsrf(i, is_sic)
1714     albsollw(i) = falblw(i, is_oce) * pctsrf(i, is_oce) &
1715     + falblw(i, is_lic) * pctsrf(i, is_lic) &
1716     + falblw(i, is_ter) * pctsrf(i, is_ter) &
1717     + falblw(i, is_sic) * pctsrf(i, is_sic)
1718     ENDDO
1719     ! nouveau rayonnement (compatible Arpege-IFS):
1720     CALL radlwsw(dist, rmu0, fract, &
1721     paprs, pplay, zxtsol, albsol, albsollw, t_seri, q_seri, &
1722     wo, &
1723     cldfra, cldemi, cldtau, &
1724     heat, heat0, cool, cool0, radsol, albpla, &
1725     topsw, toplw, solsw, sollw, &
1726     sollwdown, &
1727     topsw0, toplw0, solsw0, sollw0, &
1728     lwdn0, lwdn, lwup0, lwup, &
1729     swdn0, swdn, swup0, swup, &
1730     ok_ade, ok_aie, & ! new for aerosol radiative effects
1731     tau_ae, piz_ae, cg_ae, &
1732     topswad, solswad, &
1733     cldtaupi, &
1734     topswai, solswai)
1735     itaprad = 0
1736     ENDIF
1737     itaprad = itaprad + 1
1738    
1739     ! Ajouter la tendance des rayonnements (tous les pas)
1740    
1741     DO k = 1, llm
1742     DO i = 1, klon
1743     t_seri(i, k) = t_seri(i, k) &
1744 guez 12 + (heat(i, k)-cool(i, k)) * pdtphys/86400.
1745 guez 3 ENDDO
1746     ENDDO
1747    
1748     IF (if_ebil >= 2) THEN
1749     ztit='after rad'
1750 guez 12 CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, pdtphys &
1751 guez 10 , t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs &
1752 guez 3 , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1753     call diagphy(airephy, ztit, ip_ebil &
1754     , topsw, toplw, solsw, sollw, zero_v &
1755     , zero_v, zero_v, zero_v, ztsol &
1756     , d_h_vcol, d_qt, d_ec &
1757     , fs_bound, fq_bound )
1758     END IF
1759    
1760     ! Calculer l'hydrologie de la surface
1761    
1762     DO i = 1, klon
1763     zxqsurf(i) = 0.0
1764     zxsnow(i) = 0.0
1765     ENDDO
1766     DO nsrf = 1, nbsrf
1767     DO i = 1, klon
1768     zxqsurf(i) = zxqsurf(i) + fqsurf(i, nsrf)*pctsrf(i, nsrf)
1769     zxsnow(i) = zxsnow(i) + fsnow(i, nsrf)*pctsrf(i, nsrf)
1770     ENDDO
1771     ENDDO
1772    
1773     ! Calculer le bilan du sol et la derive de temperature (couplage)
1774    
1775     DO i = 1, klon
1776     bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
1777     ENDDO
1778    
1779 guez 13 !mod deb lott(jan95)
1780 guez 3 ! Appeler le programme de parametrisation de l'orographie
1781     ! a l'echelle sous-maille:
1782    
1783     IF (ok_orodr) THEN
1784     ! selection des points pour lesquels le shema est actif:
1785     igwd=0
1786     DO i=1, klon
1787     itest(i)=0
1788     IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
1789     itest(i)=1
1790     igwd=igwd+1
1791     idx(igwd)=i
1792     ENDIF
1793     ENDDO
1794    
1795 guez 12 CALL drag_noro(klon, llm, pdtphys, paprs, pplay, &
1796 guez 3 zmea, zstd, zsig, zgam, zthe, zpic, zval, &
1797     igwd, idx, itest, &
1798     t_seri, u_seri, v_seri, &
1799     zulow, zvlow, zustrdr, zvstrdr, &
1800     d_t_oro, d_u_oro, d_v_oro)
1801    
1802     ! ajout des tendances
1803     DO k = 1, llm
1804     DO i = 1, klon
1805     t_seri(i, k) = t_seri(i, k) + d_t_oro(i, k)
1806     u_seri(i, k) = u_seri(i, k) + d_u_oro(i, k)
1807     v_seri(i, k) = v_seri(i, k) + d_v_oro(i, k)
1808     ENDDO
1809     ENDDO
1810 guez 13 ENDIF
1811 guez 3
1812     IF (ok_orolf) THEN
1813    
1814     ! selection des points pour lesquels le shema est actif:
1815     igwd=0
1816     DO i=1, klon
1817     itest(i)=0
1818     IF ((zpic(i)-zmea(i)).GT.100.) THEN
1819     itest(i)=1
1820     igwd=igwd+1
1821     idx(igwd)=i
1822     ENDIF
1823     ENDDO
1824    
1825 guez 12 CALL lift_noro(klon, llm, pdtphys, paprs, pplay, &
1826 guez 3 rlat, zmea, zstd, zpic, &
1827     itest, &
1828     t_seri, u_seri, v_seri, &
1829     zulow, zvlow, zustrli, zvstrli, &
1830     d_t_lif, d_u_lif, d_v_lif)
1831    
1832     ! ajout des tendances
1833     DO k = 1, llm
1834     DO i = 1, klon
1835     t_seri(i, k) = t_seri(i, k) + d_t_lif(i, k)
1836     u_seri(i, k) = u_seri(i, k) + d_u_lif(i, k)
1837     v_seri(i, k) = v_seri(i, k) + d_v_lif(i, k)
1838     ENDDO
1839     ENDDO
1840    
1841     ENDIF ! fin de test sur ok_orolf
1842    
1843     ! STRESS NECESSAIRES: TOUTE LA PHYSIQUE
1844    
1845     DO i = 1, klon
1846     zustrph(i)=0.
1847     zvstrph(i)=0.
1848     ENDDO
1849     DO k = 1, llm
1850     DO i = 1, klon
1851 guez 12 zustrph(i)=zustrph(i)+(u_seri(i, k)-u(i, k))/pdtphys* &
1852 guez 3 (paprs(i, k)-paprs(i, k+1))/rg
1853 guez 12 zvstrph(i)=zvstrph(i)+(v_seri(i, k)-v(i, k))/pdtphys* &
1854 guez 3 (paprs(i, k)-paprs(i, k+1))/rg
1855     ENDDO
1856     ENDDO
1857    
1858     !IM calcul composantes axiales du moment angulaire et couple des montagnes
1859    
1860 guez 7 CALL aaam_bud (27, klon, llm, gmtime, &
1861 guez 3 ra, rg, romega, &
1862     rlat, rlon, pphis, &
1863     zustrdr, zustrli, zustrph, &
1864     zvstrdr, zvstrli, zvstrph, &
1865     paprs, u, v, &
1866     aam, torsfc)
1867    
1868     IF (if_ebil >= 2) THEN
1869     ztit='after orography'
1870 guez 12 CALL diagetpq(airephy, ztit, ip_ebil, 2, 2, pdtphys &
1871 guez 10 , t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs &
1872 guez 3 , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1873     END IF
1874    
1875     !AA Installation de l'interface online-offline pour traceurs
1876    
1877     ! Calcul des tendances traceurs
1878    
1879 guez 7 call phytrac(rnpb, itap, lmt_pas, julien, gmtime, firstcal, lafin, nq-2, &
1880 guez 12 pdtphys, u, v, t, paprs, pplay, pmfu, pmfd, pen_u, pde_u, pen_d, &
1881     pde_d, ycoefh, fm_therm, entr_therm, yu1, yv1, ftsol, pctsrf, &
1882     frac_impa, frac_nucl, presnivs, pphis, pphi, albsol, rhcl, cldfra, &
1883     rneb, diafra, cldliq, itop_con, ibas_con, pmflxr, pmflxs, prfl, &
1884     psfl, da, phi, mp, upwd, dnwd, tr_seri)
1885 guez 3
1886     IF (offline) THEN
1887    
1888     print*, 'Attention on met a 0 les thermiques pour phystoke'
1889     call phystokenc(pdtphys, rlon, rlat, &
1890     t, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
1891     fm_therm, entr_therm, &
1892     ycoefh, yu1, yv1, ftsol, pctsrf, &
1893     frac_impa, frac_nucl, &
1894 guez 12 pphis, airephy, pdtphys, itap)
1895 guez 3
1896     ENDIF
1897    
1898     ! Calculer le transport de l'eau et de l'energie (diagnostique)
1899    
1900     CALL transp (paprs, zxtsol, &
1901     t_seri, q_seri, u_seri, v_seri, zphi, &
1902     ve, vq, ue, uq)
1903    
1904     !IM diag. bilKP
1905    
1906     CALL transp_lay (paprs, zxtsol, &
1907     t_seri, q_seri, u_seri, v_seri, zphi, &
1908     ve_lay, vq_lay, ue_lay, uq_lay)
1909    
1910     ! Accumuler les variables a stocker dans les fichiers histoire:
1911    
1912     !+jld ec_conser
1913     DO k = 1, llm
1914     DO i = 1, klon
1915     ZRCPD = RCPD*(1.0+RVTMP2*q_seri(i, k))
1916     d_t_ec(i, k)=0.5/ZRCPD &
1917     *(u(i, k)**2+v(i, k)**2-u_seri(i, k)**2-v_seri(i, k)**2)
1918     t_seri(i, k)=t_seri(i, k)+d_t_ec(i, k)
1919 guez 12 d_t_ec(i, k) = d_t_ec(i, k)/pdtphys
1920 guez 3 END DO
1921     END DO
1922     !-jld ec_conser
1923     IF (if_ebil >= 1) THEN
1924     ztit='after physic'
1925 guez 12 CALL diagetpq(airephy, ztit, ip_ebil, 1, 1, pdtphys &
1926 guez 10 , t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs &
1927 guez 3 , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1928     ! Comme les tendances de la physique sont ajoute dans la dynamique,
1929     ! on devrait avoir que la variation d'entalpie par la dynamique
1930     ! est egale a la variation de la physique au pas de temps precedent.
1931     ! Donc la somme de ces 2 variations devrait etre nulle.
1932     call diagphy(airephy, ztit, ip_ebil &
1933     , topsw, toplw, solsw, sollw, sens &
1934     , evap, rain_fall, snow_fall, ztsol &
1935     , d_h_vcol, d_qt, d_ec &
1936     , fs_bound, fq_bound )
1937    
1938     d_h_vcol_phy=d_h_vcol
1939    
1940     END IF
1941    
1942     ! SORTIES
1943    
1944     !cc prw = eau precipitable
1945     DO i = 1, klon
1946     prw(i) = 0.
1947     DO k = 1, llm
1948     prw(i) = prw(i) + &
1949     q_seri(i, k)*(paprs(i, k)-paprs(i, k+1))/RG
1950     ENDDO
1951     ENDDO
1952    
1953     ! Convertir les incrementations en tendances
1954    
1955     DO k = 1, llm
1956     DO i = 1, klon
1957 guez 12 d_u(i, k) = ( u_seri(i, k) - u(i, k) ) / pdtphys
1958     d_v(i, k) = ( v_seri(i, k) - v(i, k) ) / pdtphys
1959     d_t(i, k) = ( t_seri(i, k)-t(i, k) ) / pdtphys
1960     d_qx(i, k, ivap) = ( q_seri(i, k) - qx(i, k, ivap) ) / pdtphys
1961     d_qx(i, k, iliq) = ( ql_seri(i, k) - qx(i, k, iliq) ) / pdtphys
1962 guez 3 ENDDO
1963     ENDDO
1964    
1965     IF (nq >= 3) THEN
1966     DO iq = 3, nq
1967     DO k = 1, llm
1968     DO i = 1, klon
1969 guez 13 d_qx(i, k, iq) = (tr_seri(i, k, iq-2) - qx(i, k, iq)) / pdtphys
1970 guez 3 ENDDO
1971     ENDDO
1972     ENDDO
1973     ENDIF
1974    
1975     ! Sauvegarder les valeurs de t et q a la fin de la physique:
1976    
1977     DO k = 1, llm
1978     DO i = 1, klon
1979     t_ancien(i, k) = t_seri(i, k)
1980     q_ancien(i, k) = q_seri(i, k)
1981     ENDDO
1982     ENDDO
1983    
1984     ! Ecriture des sorties
1985    
1986     call write_histhf
1987     call write_histday
1988     call write_histins
1989    
1990     ! Si c'est la fin, il faut conserver l'etat de redemarrage
1991    
1992     IF (lafin) THEN
1993     itau_phy = itau_phy + itap
1994 guez 13 CALL phyredem("restartphy.nc", rlat, rlon, pctsrf, ftsol, &
1995 guez 12 ftsoil, tslab, seaice, fqsurf, qsol, &
1996 guez 3 fsnow, falbe, falblw, fevap, rain_fall, snow_fall, &
1997     solsw, sollwdown, dlw, &
1998     radsol, frugs, agesno, &
1999 guez 13 zmea, zstd, zsig, zgam, zthe, zpic, zval, &
2000 guez 3 t_ancien, q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0)
2001     ENDIF
2002    
2003     contains
2004    
2005 guez 15 subroutine write_histday
2006 guez 3
2007 guez 15 use grid_change, only: gr_phy_write_3d
2008 guez 3
2009 guez 15 !------------------------------------------------
2010 guez 3
2011     if (ok_journe) THEN
2012     ndex2d = 0
2013     ndex3d = 0
2014     itau_w = itau_phy + itap
2015 guez 15 call histwrite(nid_day, "Sigma_O3_Royer", itau_w, gr_phy_write_3d(wo))
2016 guez 3
2017     if (ok_sync) then
2018     call histsync(nid_day)
2019     endif
2020     ENDIF
2021    
2022     End subroutine write_histday
2023    
2024     !****************************
2025    
2026     subroutine write_histhf
2027    
2028     ! From phylmd/write_histhf.h, v 1.5 2005/05/25 13:10:09
2029    
2030     ndex2d = 0
2031     ndex3d = 0
2032    
2033     itau_w = itau_phy + itap
2034    
2035     call write_histhf3d
2036    
2037     IF (ok_sync) THEN
2038     call histsync(nid_hf)
2039     ENDIF
2040    
2041     end subroutine write_histhf
2042    
2043     !***************************************************************
2044    
2045     subroutine write_histins
2046    
2047     ! From phylmd/write_histins.h, v 1.2 2005/05/25 13:10:09
2048    
2049     real zout
2050    
2051     !--------------------------------------------------
2052    
2053     IF (ok_instan) THEN
2054    
2055     ndex2d = 0
2056     ndex3d = 0
2057    
2058     ! Champs 2D:
2059    
2060 guez 12 zsto = pdtphys * ecrit_ins
2061     zout = pdtphys * ecrit_ins
2062 guez 3 itau_w = itau_phy + itap
2063    
2064     i = NINT(zout/zsto)
2065     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), pphis, zx_tmp_2d)
2066 guez 15 CALL histwrite(nid_ins, "phis", itau_w, zx_tmp_2d)
2067 guez 3
2068     i = NINT(zout/zsto)
2069     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), airephy, zx_tmp_2d)
2070 guez 15 CALL histwrite(nid_ins, "aire", itau_w, zx_tmp_2d)
2071 guez 3
2072     DO i = 1, klon
2073     zx_tmp_fi2d(i) = paprs(i, 1)
2074     ENDDO
2075     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)
2076 guez 15 CALL histwrite(nid_ins, "psol", itau_w, zx_tmp_2d)
2077 guez 3
2078     DO i = 1, klon
2079     zx_tmp_fi2d(i) = rain_fall(i) + snow_fall(i)
2080     ENDDO
2081     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)
2082 guez 15 CALL histwrite(nid_ins, "precip", itau_w, zx_tmp_2d)
2083 guez 3
2084     DO i = 1, klon
2085     zx_tmp_fi2d(i) = rain_lsc(i) + snow_lsc(i)
2086     ENDDO
2087     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)
2088 guez 15 CALL histwrite(nid_ins, "plul", itau_w, zx_tmp_2d)
2089 guez 3
2090     DO i = 1, klon
2091     zx_tmp_fi2d(i) = rain_con(i) + snow_con(i)
2092     ENDDO
2093     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)
2094 guez 15 CALL histwrite(nid_ins, "pluc", itau_w, zx_tmp_2d)
2095 guez 3
2096     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zxtsol, zx_tmp_2d)
2097 guez 15 CALL histwrite(nid_ins, "tsol", itau_w, zx_tmp_2d)
2098 guez 3 !ccIM
2099     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zt2m, zx_tmp_2d)
2100 guez 15 CALL histwrite(nid_ins, "t2m", itau_w, zx_tmp_2d)
2101 guez 3
2102     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zq2m, zx_tmp_2d)
2103 guez 15 CALL histwrite(nid_ins, "q2m", itau_w, zx_tmp_2d)
2104 guez 3
2105     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zu10m, zx_tmp_2d)
2106 guez 15 CALL histwrite(nid_ins, "u10m", itau_w, zx_tmp_2d)
2107 guez 3
2108     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zv10m, zx_tmp_2d)
2109 guez 15 CALL histwrite(nid_ins, "v10m", itau_w, zx_tmp_2d)
2110 guez 3
2111     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), snow_fall, zx_tmp_2d)
2112 guez 15 CALL histwrite(nid_ins, "snow", itau_w, zx_tmp_2d)
2113 guez 3
2114     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), cdragm, zx_tmp_2d)
2115 guez 15 CALL histwrite(nid_ins, "cdrm", itau_w, zx_tmp_2d)
2116 guez 3
2117     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), cdragh, zx_tmp_2d)
2118 guez 15 CALL histwrite(nid_ins, "cdrh", itau_w, zx_tmp_2d)
2119 guez 3
2120     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), toplw, zx_tmp_2d)
2121 guez 15 CALL histwrite(nid_ins, "topl", itau_w, zx_tmp_2d)
2122 guez 3
2123     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), evap, zx_tmp_2d)
2124 guez 15 CALL histwrite(nid_ins, "evap", itau_w, zx_tmp_2d)
2125 guez 3
2126     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), solsw, zx_tmp_2d)
2127 guez 15 CALL histwrite(nid_ins, "sols", itau_w, zx_tmp_2d)
2128 guez 3
2129     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), sollw, zx_tmp_2d)
2130 guez 15 CALL histwrite(nid_ins, "soll", itau_w, zx_tmp_2d)
2131 guez 3
2132     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), sollwdown, zx_tmp_2d)
2133 guez 15 CALL histwrite(nid_ins, "solldown", itau_w, zx_tmp_2d)
2134 guez 3
2135     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), bils, zx_tmp_2d)
2136 guez 15 CALL histwrite(nid_ins, "bils", itau_w, zx_tmp_2d)
2137 guez 3
2138     zx_tmp_fi2d(1:klon)=-1*sens(1:klon)
2139     ! CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), sens, zx_tmp_2d)
2140     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)
2141 guez 15 CALL histwrite(nid_ins, "sens", itau_w, zx_tmp_2d)
2142 guez 3
2143     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), fder, zx_tmp_2d)
2144 guez 15 CALL histwrite(nid_ins, "fder", itau_w, zx_tmp_2d)
2145 guez 3
2146     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), d_ts(1, is_oce), zx_tmp_2d)
2147 guez 15 CALL histwrite(nid_ins, "dtsvdfo", itau_w, zx_tmp_2d)
2148 guez 3
2149     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), d_ts(1, is_ter), zx_tmp_2d)
2150 guez 15 CALL histwrite(nid_ins, "dtsvdft", itau_w, zx_tmp_2d)
2151 guez 3
2152     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), d_ts(1, is_lic), zx_tmp_2d)
2153 guez 15 CALL histwrite(nid_ins, "dtsvdfg", itau_w, zx_tmp_2d)
2154 guez 3
2155     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), d_ts(1, is_sic), zx_tmp_2d)
2156 guez 15 CALL histwrite(nid_ins, "dtsvdfi", itau_w, zx_tmp_2d)
2157 guez 3
2158     DO nsrf = 1, nbsrf
2159     !XXX
2160     zx_tmp_fi2d(1 : klon) = pctsrf( 1 : klon, nsrf)*100.
2161     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)
2162     CALL histwrite(nid_ins, "pourc_"//clnsurf(nsrf), itau_w, &
2163 guez 15 zx_tmp_2d)
2164 guez 3
2165     zx_tmp_fi2d(1 : klon) = pctsrf( 1 : klon, nsrf)
2166     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)
2167     CALL histwrite(nid_ins, "fract_"//clnsurf(nsrf), itau_w, &
2168 guez 15 zx_tmp_2d)
2169 guez 3
2170     zx_tmp_fi2d(1 : klon) = fluxt( 1 : klon, 1, nsrf)
2171     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)
2172     CALL histwrite(nid_ins, "sens_"//clnsurf(nsrf), itau_w, &
2173 guez 15 zx_tmp_2d)
2174 guez 3
2175     zx_tmp_fi2d(1 : klon) = fluxlat( 1 : klon, nsrf)
2176     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)
2177     CALL histwrite(nid_ins, "lat_"//clnsurf(nsrf), itau_w, &
2178 guez 15 zx_tmp_2d)
2179 guez 3
2180     zx_tmp_fi2d(1 : klon) = ftsol( 1 : klon, nsrf)
2181     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)
2182     CALL histwrite(nid_ins, "tsol_"//clnsurf(nsrf), itau_w, &
2183 guez 15 zx_tmp_2d)
2184 guez 3
2185     zx_tmp_fi2d(1 : klon) = fluxu( 1 : klon, 1, nsrf)
2186     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)
2187     CALL histwrite(nid_ins, "taux_"//clnsurf(nsrf), itau_w, &
2188 guez 15 zx_tmp_2d)
2189 guez 3
2190     zx_tmp_fi2d(1 : klon) = fluxv( 1 : klon, 1, nsrf)
2191     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)
2192     CALL histwrite(nid_ins, "tauy_"//clnsurf(nsrf), itau_w, &
2193 guez 15 zx_tmp_2d)
2194 guez 3
2195     zx_tmp_fi2d(1 : klon) = frugs( 1 : klon, nsrf)
2196     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)
2197     CALL histwrite(nid_ins, "rugs_"//clnsurf(nsrf), itau_w, &
2198 guez 15 zx_tmp_2d)
2199 guez 3
2200     zx_tmp_fi2d(1 : klon) = falbe( 1 : klon, nsrf)
2201     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zx_tmp_fi2d, zx_tmp_2d)
2202     CALL histwrite(nid_ins, "albe_"//clnsurf(nsrf), itau_w, &
2203 guez 15 zx_tmp_2d)
2204 guez 3
2205     END DO
2206     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), albsol, zx_tmp_2d)
2207 guez 15 CALL histwrite(nid_ins, "albs", itau_w, zx_tmp_2d)
2208 guez 3 CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), albsollw, zx_tmp_2d)
2209 guez 15 CALL histwrite(nid_ins, "albslw", itau_w, zx_tmp_2d)
2210 guez 3
2211     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), zxrugs, zx_tmp_2d)
2212 guez 15 CALL histwrite(nid_ins, "rugs", itau_w, zx_tmp_2d)
2213 guez 3
2214     !IM cf. AM 081204 BEG
2215    
2216     !HBTM2
2217    
2218     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), s_pblh, zx_tmp_2d)
2219 guez 15 CALL histwrite(nid_ins, "s_pblh", itau_w, zx_tmp_2d)
2220 guez 3
2221     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), s_pblt, zx_tmp_2d)
2222 guez 15 CALL histwrite(nid_ins, "s_pblt", itau_w, zx_tmp_2d)
2223 guez 3
2224     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), s_lcl, zx_tmp_2d)
2225 guez 15 CALL histwrite(nid_ins, "s_lcl", itau_w, zx_tmp_2d)
2226 guez 3
2227     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), s_capCL, zx_tmp_2d)
2228 guez 15 CALL histwrite(nid_ins, "s_capCL", itau_w, zx_tmp_2d)
2229 guez 3
2230     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), s_oliqCL, zx_tmp_2d)
2231 guez 15 CALL histwrite(nid_ins, "s_oliqCL", itau_w, zx_tmp_2d)
2232 guez 3
2233     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), s_cteiCL, zx_tmp_2d)
2234 guez 15 CALL histwrite(nid_ins, "s_cteiCL", itau_w, zx_tmp_2d)
2235 guez 3
2236     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), s_therm, zx_tmp_2d)
2237 guez 15 CALL histwrite(nid_ins, "s_therm", itau_w, zx_tmp_2d)
2238 guez 3
2239     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), s_trmb1, zx_tmp_2d)
2240 guez 15 CALL histwrite(nid_ins, "s_trmb1", itau_w, zx_tmp_2d)
2241 guez 3
2242     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), s_trmb2, zx_tmp_2d)
2243 guez 15 CALL histwrite(nid_ins, "s_trmb2", itau_w, zx_tmp_2d)
2244 guez 3
2245     CALL gr_fi_ecrit(1, klon, iim, (jjm + 1), s_trmb3, zx_tmp_2d)
2246 guez 15 CALL histwrite(nid_ins, "s_trmb3", itau_w, zx_tmp_2d)
2247 guez 3
2248     !IM cf. AM 081204 END
2249    
2250     ! Champs 3D:
2251    
2252     CALL gr_fi_ecrit(llm, klon, iim, (jjm + 1), t_seri, zx_tmp_3d)
2253 guez 15 CALL histwrite(nid_ins, "temp", itau_w, zx_tmp_3d)
2254 guez 3
2255     CALL gr_fi_ecrit(llm, klon, iim, (jjm + 1), u_seri, zx_tmp_3d)
2256 guez 15 CALL histwrite(nid_ins, "vitu", itau_w, zx_tmp_3d)
2257 guez 3
2258     CALL gr_fi_ecrit(llm, klon, iim, (jjm + 1), v_seri, zx_tmp_3d)
2259 guez 15 CALL histwrite(nid_ins, "vitv", itau_w, zx_tmp_3d)
2260 guez 3
2261     CALL gr_fi_ecrit(llm, klon, iim, (jjm + 1), zphi, zx_tmp_3d)
2262 guez 15 CALL histwrite(nid_ins, "geop", itau_w, zx_tmp_3d)
2263 guez 3
2264     CALL gr_fi_ecrit(llm, klon, iim, (jjm + 1), pplay, zx_tmp_3d)
2265 guez 15 CALL histwrite(nid_ins, "pres", itau_w, zx_tmp_3d)
2266 guez 3
2267     CALL gr_fi_ecrit(llm, klon, iim, (jjm + 1), d_t_vdf, zx_tmp_3d)
2268 guez 15 CALL histwrite(nid_ins, "dtvdf", itau_w, zx_tmp_3d)
2269 guez 3
2270     CALL gr_fi_ecrit(llm, klon, iim, (jjm + 1), d_q_vdf, zx_tmp_3d)
2271 guez 15 CALL histwrite(nid_ins, "dqvdf", itau_w, zx_tmp_3d)
2272 guez 3
2273     if (ok_sync) then
2274     call histsync(nid_ins)
2275     endif
2276     ENDIF
2277    
2278     end subroutine write_histins
2279    
2280     !****************************************************
2281    
2282     subroutine write_histhf3d
2283    
2284     ! From phylmd/write_histhf3d.h, v 1.2 2005/05/25 13:10:09
2285    
2286     ndex2d = 0
2287     ndex3d = 0
2288    
2289     itau_w = itau_phy + itap
2290    
2291     ! Champs 3D:
2292    
2293     CALL gr_fi_ecrit(llm, klon, iim, (jjm + 1), t_seri, zx_tmp_3d)
2294 guez 15 CALL histwrite(nid_hf3d, "temp", itau_w, zx_tmp_3d)
2295 guez 3
2296     CALL gr_fi_ecrit(llm, klon, iim, (jjm + 1), qx(1, 1, ivap), zx_tmp_3d)
2297 guez 15 CALL histwrite(nid_hf3d, "ovap", itau_w, zx_tmp_3d)
2298 guez 3
2299     CALL gr_fi_ecrit(llm, klon, iim, (jjm + 1), u_seri, zx_tmp_3d)
2300 guez 15 CALL histwrite(nid_hf3d, "vitu", itau_w, zx_tmp_3d)
2301 guez 3
2302     CALL gr_fi_ecrit(llm, klon, iim, (jjm + 1), v_seri, zx_tmp_3d)
2303 guez 15 CALL histwrite(nid_hf3d, "vitv", itau_w, zx_tmp_3d)
2304 guez 3
2305 guez 6 if (nbtr >= 3) then
2306     CALL gr_fi_ecrit(llm, klon, iim, (jjm + 1), tr_seri(1, 1, 3), &
2307     zx_tmp_3d)
2308 guez 15 CALL histwrite(nid_hf3d, "O3", itau_w, zx_tmp_3d)
2309 guez 6 end if
2310 guez 3
2311     if (ok_sync) then
2312     call histsync(nid_hf3d)
2313     endif
2314    
2315     end subroutine write_histhf3d
2316    
2317     END SUBROUTINE physiq
2318    
2319     !****************************************************
2320    
2321     FUNCTION qcheck(klon, klev, paprs, q, ql, aire)
2322    
2323     ! From phylmd/physiq.F, v 1.22 2006/02/20 09:38:28
2324    
2325     use YOMCST
2326     IMPLICIT none
2327    
2328     ! Calculer et imprimer l'eau totale. A utiliser pour verifier
2329     ! la conservation de l'eau
2330    
2331     INTEGER klon, klev
2332     REAL, intent(in):: paprs(klon, klev+1)
2333     real q(klon, klev), ql(klon, klev)
2334     REAL aire(klon)
2335     REAL qtotal, zx, qcheck
2336     INTEGER i, k
2337    
2338     zx = 0.0
2339     DO i = 1, klon
2340     zx = zx + aire(i)
2341     ENDDO
2342     qtotal = 0.0
2343     DO k = 1, klev
2344     DO i = 1, klon
2345     qtotal = qtotal + (q(i, k)+ql(i, k)) * aire(i) &
2346     *(paprs(i, k)-paprs(i, k+1))/RG
2347     ENDDO
2348     ENDDO
2349    
2350     qcheck = qtotal/zx
2351    
2352     END FUNCTION qcheck
2353    
2354     end module physiq_m
  ViewVC Help
Powered by ViewVC 1.1.21