/[lmdze]/trunk/phylmd/physiq.f
ViewVC logotype

Annotation of /trunk/phylmd/physiq.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 23 - (hide annotations)
Mon Dec 14 15:25:16 2009 UTC (14 years, 5 months ago) by guez
Original Path: trunk/libf/phylmd/physiq.f90
File size: 75335 byte(s) 
Split "orografi.f": one file for each procedure. Put the created files
in new directory "Orography".

Removed argument "vcov" of procedure "sortvarc". Removed arguments
"itau" and "time" of procedure "caldyn0". Removed arguments "itau",
"time" and "vcov" of procedure "sortvarc0".

Removed argument "time" of procedure "dynredem1". Removed NetCDF
variable "temps" in files "start.nc" and "restart.nc", because its
value is always 0.

Removed argument "nq" of procedures "iniadvtrac" and "leapfrog". The
number of "tracers read in "traceur.def" must now be equal to "nqmx",
or "nqmx" must equal 4 if there is no file "traceur.def". Replaced
variable "nq" by constant "nqmx" in "leapfrog".

NetCDF variable for ozone field in "coefoz.nc" must now be called
"tro3" instead of "r".

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