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

Contents of /trunk/phylmd/physiq.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 108 - (show annotations)
Tue Sep 16 14:00:41 2014 UTC (9 years, 7 months ago) by guez
File size: 61789 byte(s)
Imported writefield from LMDZ. Close at the end of gcm the files which
were created by writefiled (not done in LMDZ).

Removed procedures for the output of Grads files. Removed calls to
dump2d. In guide, replaced calls to wrgrads by calls to writefield.

In vlspltqs, removed redundant programming of saturation
pressure. Call foeew from module FCTTRE instead.

Bug fix in interpre: size of w exceeding size of correponding actual
argument wg in advtrac.

In leapfrog, call guide until the end of the run, instead of six hours
before the end.

Bug fix in readsulfate_preind: type of arguments.

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

  ViewVC Help
Powered by ViewVC 1.1.21