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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 53 - (hide annotations)
Fri Oct 7 13:11:58 2011 UTC (12 years, 7 months ago) by guez
File size: 71771 byte(s)


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

  ViewVC Help
Powered by ViewVC 1.1.21