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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 38 - (hide annotations)
Thu Jan 6 17:52:19 2011 UTC (13 years, 4 months ago) by guez
File size: 75225 byte(s)
Extracted ASCII art from "inigeom" into a separate text file in the
documentation.

"test_disvert" now creates a separate file for layer thicknesses.

Moved variables from module "yomcst" to module "suphec_m" because this
is where those variables are defined. Kept in "yomcst" only parameters
of Earth orbit. Gave the attribute "parameter" to some variables of
module "suphec_m".

Variables of module "yoethf" were defined in procedure "suphec". Moved
these definitions to a new procedure "yoethf" in module "yoethf_m".

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

  ViewVC Help
Powered by ViewVC 1.1.21