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

Contents of /trunk/phylmd/physiq.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 103 - (show annotations)
Fri Aug 29 13:00:05 2014 UTC (9 years, 8 months ago) by guez
File size: 61734 byte(s)
Renamed module cvparam to cv_param. Deleted procedure
cv_param. Changed variables of module cv_param into parameters.

In procedures cv_driver, cv_uncompress and cv3_uncompress, removed
some arguments giving dimensions and used module variables klon and
klev instead.

In procedures gradiv2, laplacien_gam and laplacien, changed
declarations of local variables because klevel is not always klev.

Removed code for nudging surface pressure.

Removed arguments pim and pjm of tau2alpha. Added assignment of false
to variable first.

Replaced real argument del of procedures foeew and FOEDE by logical
argument.

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

  ViewVC Help
Powered by ViewVC 1.1.21