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

Contents of /trunk/phylmd/phytrac.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 266 - (show annotations)
Thu Apr 19 17:54:55 2018 UTC (5 years, 11 months ago) by guez
File size: 13289 byte(s)
Define macros of the preprocessor CPP_IIM, CPP_JJM, CPP_LLM so we can
control the resolution from the compilation command, and automate
compilation for several resolutions.

In module yoethf_m, transform variables into named constants. So we do
not need procedure yoethf any longer.

Bug fix in program test_inter_barxy, missing calls to fyhyp and fxhyp,
and definition of rlatu.

Remove variable iecri of module conf_gcm_m. The files dyn_hist*.nc are
written every time step. We are simplifying the output system, pending
replacement by a whole new system.

Modify possible value of vert_sampling from "param" to
"strato_custom", following LMDZ. Default values of corresponding
namelist variables are now the values used for LMDZ CMIP6.

1 module phytrac_m
2
3 IMPLICIT none
4
5 private
6 public phytrac
7
8 contains
9
10 SUBROUTINE phytrac(julien, gmtime, firstcal, lafin, pdtphys, t_seri, paprs, &
11 pplay, pmfu, pmfd, pde_u, pen_d, coefh, cdragh, fm_therm, entr_therm, &
12 yu1, yv1, ftsol, pctsrf, frac_impa, frac_nucl, da, phi, mp, upwd, dnwd, &
13 tr_seri, zmasse, ncid_startphy)
14
15 ! From phylmd/phytrac.F, version 1.15, 2006/02/21 08:08:30 (SVN
16 ! revision 679) and phylmd/write_histrac.h, version 1.9,
17 ! 2006/02/21 08:08:30
18
19 ! Authors: Fr\'ed\'eric Hourdin, Abderrahmane Idelkadi, Marie-Alice
20 ! Foujols
21
22 ! Objet : moniteur g\'en\'eral des tendances des traceurs
23
24 ! L'appel de "phytrac" se fait avec "nqmx - 2" donc nous avons
25 ! bien les vrais traceurs, sans la vapeur d'eau ni l'eau liquide.
26
27 ! Modifications pour les traceurs :
28 ! - uniformisation des param\'etrisations dans phytrac
29 ! - stockage des moyennes des champs n\'ecessaires en mode traceur off-line
30
31 use abort_gcm_m, only: abort_gcm
32 use clesphys2, only: conv_emanuel
33 use cltrac_m, only: cltrac
34 use cltracrn_m, only: cltracrn
35 USE conf_gcm_m, ONLY: lmt_pas
36 use ctherm, only: iflag_thermals
37 use cvltr_m, only: cvltr
38 use dimensions, only: llm, nqmx
39 use dimphy, only: klon
40 use histwrite_phy_m, only: histwrite_phy
41 use indicesol, only: nbsrf
42 use iniadvtrac_m, only: tname
43 use initrrnpb_m, only: initrrnpb
44 use minmaxqfi_m, only: minmaxqfi
45 use netcdf, only: NF90_FILL_float
46 use netcdf95, only: nf95_inq_varid, nf95_get_var, nf95_put_var
47 use nflxtr_m, only: nflxtr
48 use nr_util, only: assert
49 use o3_chem_m, only: o3_chem
50 use phyetat0_m, only: rlat
51 use phyredem0_m, only: ncid_restartphy
52 use press_coefoz_m, only: press_coefoz
53 use radiornpb_m, only: radiornpb
54 use regr_pr_comb_coefoz_m, only: regr_pr_comb_coefoz
55 use SUPHEC_M, only: rg
56 use time_phylmdz, only: itap
57
58 integer, intent(in):: julien !jour julien, 1 <= julien <= 360
59 real, intent(in):: gmtime ! heure de la journ\'ee en fraction de jour
60 logical, intent(in):: firstcal ! first call to "calfis"
61 logical, intent(in):: lafin ! fin de la physique
62 real, intent(in):: pdtphys ! pas d'integration pour la physique (s)
63 real, intent(in):: t_seri(klon, llm) ! temperature, in K
64
65 real, intent(in):: paprs(klon, llm+1)
66 ! (pression pour chaque inter-couche, en Pa)
67
68 real, intent(in):: pplay(klon, llm)
69 ! (pression pour le mileu de chaque couche, en Pa)
70
71 ! convection:
72
73 REAL, intent(in):: pmfu(klon, llm) ! flux de masse dans le panache montant
74
75 REAL, intent(in):: pmfd(klon, llm)
76 ! flux de masse dans le panache descendant
77
78 REAL pde_u(klon, llm) ! flux detraine dans le panache montant
79 REAL pen_d(klon, llm) ! flux entraine dans le panache descendant
80 REAL coefh(:, 2:) ! (klon, 2:llm) coeff melange couche limite
81 real cdragh(:) ! (klon)
82 real fm_therm(klon, llm+1), entr_therm(klon, llm) ! thermiques
83 REAL, intent(in):: yu1(:), yv1(:) ! (klon) vent au premier niveau
84
85 ! Arguments n\'ecessaires pour les sources et puits de traceur :
86 real, intent(in):: ftsol(:, :) ! (klon, nbsrf) surface temperature (K)
87 real, intent(in):: pctsrf(klon, nbsrf) ! Pourcentage de sol f(nature du sol)
88
89 ! Lessivage pour le on-line
90 REAL, intent(in):: frac_impa(klon, llm) ! fraction d'aerosols impactes
91 REAL, intent(in):: frac_nucl(klon, llm) ! fraction d'aerosols nuclees
92
93 ! Kerry Emanuel
94 real, intent(in):: da(klon, llm), phi(klon, llm, llm), mp(klon, llm)
95 REAL, intent(in):: upwd(klon, llm) ! saturated updraft mass flux
96 REAL, intent(in):: dnwd(klon, llm) ! saturated downdraft mass flux
97
98 real, intent(inout):: tr_seri(:, :, :) ! (klon, llm, nqmx - 2)
99 ! (mass fractions of tracers, excluding water, at mid-layers)
100
101 real, intent(in):: zmasse(:, :) ! (klon, llm)
102 ! (column-density of mass of air in a cell, in kg m-2)
103
104 integer, intent(in):: ncid_startphy
105
106 ! Local:
107
108 integer nsplit
109
110 ! TRACEURS
111
112 ! Sources et puits des traceurs:
113
114 ! Pour l'instant seuls les cas du rn et du pb ont ete envisages.
115
116 REAL source(klon) ! a voir lorsque le flux est prescrit
117 !
118 ! Pour la source de radon et son reservoir de sol
119
120 REAL, save:: trs(klon, nqmx - 2) ! Concentration de traceur dans le sol
121
122 REAL, save:: masktr(klon, nqmx - 2) ! Masque reservoir de sol traceur
123 ! Masque de l'echange avec la surface
124 ! (1 = reservoir) ou (possible => 1)
125
126 REAL, save:: fshtr(klon, nqmx - 2)
127 ! Flux surfacique dans le reservoir de sol
128
129 REAL, save:: hsoltr(nqmx - 2) ! Epaisseur equivalente du reservoir de sol
130 REAL, save:: tautr(nqmx - 2) ! constante de d\'ecroissance radioactive
131
132 REAL, save:: vdeptr(nqmx - 2)
133 ! Vitesse de depot sec dans la couche Brownienne
134
135 REAL, save:: scavtr(nqmx - 2) ! Coefficient de lessivage
136 CHARACTER itn
137
138 logical, save:: aerosol(nqmx - 2) ! Nature du traceur
139 ! ! aerosol(it) = true => aerosol
140 ! ! aerosol(it) = false => gaz
141
142 logical, save:: clsol(nqmx - 2) ! couche limite sol flux
143 ! calcul\'ee, sinon prescrit
144 logical, save:: radio(nqmx - 2) ! d\'ecroisssance radioactive
145
146 ! convection tiedtke
147 INTEGER i, k, it
148 REAL delp(klon, llm)
149
150 ! Variables liees a l'ecriture de la bande histoire physique
151
152 ! Variables locales pour effectuer les appels en serie
153
154 REAL d_tr(klon, llm), d_trs(klon) ! tendances de traceurs
155 REAL d_tr_cl(klon, llm, nqmx - 2) ! tendance de traceurs couche limite
156
157 REAL d_tr_cv(klon, llm, nqmx - 2)
158 ! tendance de traceurs conv pour chq traceur
159
160 REAL d_tr_th(klon, llm, nqmx - 2) ! la tendance des thermiques
161 REAL d_tr_dec(klon, llm, 2) ! la tendance de la decroissance
162 ! ! radioactive du rn - > pb
163
164 REAL d_tr_lessi_impa(klon, llm, nqmx - 2)
165 ! tendance du lessivage par impaction
166
167 REAL d_tr_lessi_nucl(klon, llm, nqmx - 2)
168 ! tendance du lessivage par nucleation
169
170 REAL flestottr(klon, llm, nqmx - 2) ! flux de lessivage dans chaque couche
171
172 real ztra_th(klon, llm)
173 integer isplit, varid
174
175 ! Controls:
176 logical:: convection = .true.
177
178 !--------------------------------------
179
180 call assert(shape(zmasse) == (/klon, llm/), "phytrac zmasse")
181 call assert(shape(tr_seri) == (/klon, llm, nqmx - 2/), "phytrac tr_seri")
182
183 if (firstcal) then
184 ! Initialisation de certaines variables pour le radon et le plomb
185 ! Initialisation du traceur dans le sol (couche limite radonique)
186 trs(:, 2:) = 0.
187
188 call nf95_inq_varid(ncid_startphy, "trs", varid)
189 call nf95_get_var(ncid_startphy, varid, trs(:, 1))
190 if (any(trs(:, 1) == NF90_FILL_float)) call abort_gcm("phytrac", &
191 "some missing values in trs(:, 1)")
192
193 ! Initialisation de la fraction d'aerosols lessivee
194
195 d_tr_lessi_impa = 0.
196 d_tr_lessi_nucl = 0.
197
198 ! Initialisation de la nature des traceurs
199
200 aerosol = .FALSE. ! Tous les traceurs sont des gaz par defaut
201 radio = .FALSE. ! par d\'efaut pas de passage par "radiornpb"
202
203 if (nqmx >= 5) then
204 call press_coefoz ! read input pressure levels for ozone coefficients
205 end if
206
207 ! Initialisation du traceur dans le sol (couche limite radonique)
208 radio(1)= .true.
209 radio(2)= .true.
210 clsol(:2)= .true.
211 clsol(3:)= .false.
212 aerosol(2) = .TRUE. ! le Pb est un aerosol
213 call initrrnpb(pctsrf, masktr, fshtr, hsoltr, tautr, vdeptr, scavtr)
214 endif
215
216 if (convection) then
217 ! Calcul de l'effet de la convection
218 DO it=1, nqmx - 2
219 if (conv_emanuel) then
220 call cvltr(pdtphys, da, phi, mp, paprs, tr_seri(:, :, it), upwd, &
221 dnwd, d_tr_cv(:, :, it))
222 else
223 CALL nflxtr(pdtphys, pmfu, pmfd, pde_u, pen_d, paprs, &
224 tr_seri(:, :, it), d_tr_cv(:, :, it))
225 endif
226
227 DO k = 1, llm
228 DO i = 1, klon
229 tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cv(i, k, it)
230 ENDDO
231 ENDDO
232 WRITE(unit=itn, fmt='(i1)') it
233 CALL minmaxqfi(tr_seri(:, :, it), 0., 1e33, &
234 'convection, tracer index = ' // itn)
235 ENDDO
236 endif
237
238 ! Calcul de l'effet des thermiques
239
240 do it=1, nqmx - 2
241 do k=1, llm
242 do i=1, klon
243 d_tr_th(i, k, it)=0.
244 tr_seri(i, k, it)=max(tr_seri(i, k, it), 0.)
245 tr_seri(i, k, it)=min(tr_seri(i, k, it), 1e10)
246 enddo
247 enddo
248 enddo
249
250 if (iflag_thermals > 0) then
251 nsplit=10
252 DO it=1, nqmx - 2
253 do isplit=1, nsplit
254 ! Thermiques
255 call dqthermcell(klon, llm, pdtphys/nsplit &
256 , fm_therm, entr_therm, zmasse &
257 , tr_seri(1:klon, 1:llm, it), d_tr, ztra_th)
258
259 do k=1, llm
260 do i=1, klon
261 d_tr(i, k)=pdtphys*d_tr(i, k)/nsplit
262 d_tr_th(i, k, it)=d_tr_th(i, k, it)+d_tr(i, k)
263 tr_seri(i, k, it)=max(tr_seri(i, k, it)+d_tr(i, k), 0.)
264 enddo
265 enddo
266 enddo
267 ENDDO
268 endif
269
270 ! Calcul de l'effet de la couche limite
271
272 DO k = 1, llm
273 DO i = 1, klon
274 delp(i, k) = paprs(i, k)-paprs(i, k+1)
275 ENDDO
276 ENDDO
277
278 ! MAF modif pour tenir compte du cas traceur
279 DO it=1, nqmx - 2
280 if (clsol(it)) then
281 ! couche limite avec quantite dans le sol calculee
282 CALL cltracrn(it, pdtphys, yu1, yv1, coefh, cdragh, t_seri, ftsol, &
283 pctsrf, tr_seri(:, :, it), trs(:, it), paprs, pplay, delp, &
284 masktr(1, it), fshtr(1, it), hsoltr(it), tautr(it), &
285 vdeptr(it), rlat, d_tr_cl(1, 1, it), d_trs)
286 DO k = 1, llm
287 DO i = 1, klon
288 tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cl(i, k, it)
289 ENDDO
290 ENDDO
291
292 trs(:, it) = trs(:, it) + d_trs
293 else
294 ! couche limite avec flux prescrit
295 !MAF provisoire source / traceur a creer
296 DO i=1, klon
297 source(i) = 0. ! pas de source, pour l'instant
298 ENDDO
299
300 CALL cltrac(pdtphys, coefh, t_seri, tr_seri(:, :, it), source, &
301 paprs, pplay, delp, d_tr_cl(1, 1, it))
302 DO k = 1, llm
303 DO i = 1, klon
304 tr_seri(i, k, it) = tr_seri(i, k, it) + d_tr_cl(i, k, it)
305 ENDDO
306 ENDDO
307 endif
308 ENDDO
309
310 ! Calcul de l'effet du puits radioactif
311
312 ! MAF il faudrait faire une modification pour passer dans radiornpb
313 ! si radio=true
314 d_tr_dec = radiornpb(tr_seri, pdtphys, tautr)
315 DO it = 1, nqmx - 2
316 if (radio(it)) then
317 tr_seri(:, :, it) = tr_seri(:, :, it) + d_tr_dec(:, :, it)
318 WRITE(unit=itn, fmt='(i1)') it
319 CALL minmaxqfi(tr_seri(:, :, it), 0., 1e33, 'puits rn it='//itn)
320 endif
321 ENDDO
322
323 if (nqmx >= 5) then
324 ! Ozone as a tracer:
325 if (mod(itap - 1, lmt_pas) == 0) then
326 ! Once per day, update the coefficients for ozone chemistry:
327 call regr_pr_comb_coefoz(julien, paprs, pplay)
328 end if
329 call o3_chem(julien, gmtime, t_seri, zmasse, pdtphys, tr_seri(:, :, 3))
330 end if
331
332 ! Calcul de l'effet de la precipitation
333
334 d_tr_lessi_nucl = 0.
335 d_tr_lessi_impa = 0.
336 flestottr = 0.
337
338 ! tendance des aerosols nuclees et impactes
339
340 DO it = 1, nqmx - 2
341 IF (aerosol(it)) THEN
342 DO k = 1, llm
343 DO i = 1, klon
344 d_tr_lessi_nucl(i, k, it) = d_tr_lessi_nucl(i, k, it) + &
345 (1 - frac_nucl(i, k))*tr_seri(i, k, it)
346 d_tr_lessi_impa(i, k, it) = d_tr_lessi_impa(i, k, it) + &
347 (1 - frac_impa(i, k))*tr_seri(i, k, it)
348 ENDDO
349 ENDDO
350 ENDIF
351 ENDDO
352
353 ! Mises a jour des traceurs + calcul des flux de lessivage
354 ! Mise a jour due a l'impaction et a la nucleation
355
356 DO it = 1, nqmx - 2
357 IF (aerosol(it)) THEN
358 DO k = 1, llm
359 DO i = 1, klon
360 tr_seri(i, k, it) = tr_seri(i, k, it) * frac_impa(i, k) &
361 * frac_nucl(i, k)
362 ENDDO
363 ENDDO
364 ENDIF
365 ENDDO
366
367 ! Flux lessivage total
368 DO it = 1, nqmx - 2
369 DO k = 1, llm
370 DO i = 1, klon
371 flestottr(i, k, it) = flestottr(i, k, it) &
372 - (d_tr_lessi_nucl(i, k, it) + d_tr_lessi_impa(i, k, it)) &
373 * (paprs(i, k)-paprs(i, k+1)) / (RG * pdtphys)
374 ENDDO
375 ENDDO
376 ENDDO
377
378 ! Ecriture des sorties
379 CALL histwrite_phy("zmasse", zmasse)
380 DO it=1, nqmx - 2
381 CALL histwrite_phy(tname(it+2), tr_seri(:, :, it))
382 CALL histwrite_phy("fl"//tname(it+2), flestottr(:, :, it))
383 CALL histwrite_phy("d_tr_th_"//tname(it+2), d_tr_th(:, :, it))
384 CALL histwrite_phy("d_tr_cv_"//tname(it+2), d_tr_cv(:, :, it))
385 CALL histwrite_phy("d_tr_cl_"//tname(it+2), d_tr_cl(:, :, it))
386 ENDDO
387
388 if (lafin) then
389 call nf95_inq_varid(ncid_restartphy, "trs", varid)
390 call nf95_put_var(ncid_restartphy, varid, trs(:, 1))
391 endif
392
393 END SUBROUTINE phytrac
394
395 end module phytrac_m

  ViewVC Help
Powered by ViewVC 1.1.21