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

Contents of /trunk/phylmd/phytrac.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 227 - (show annotations)
Thu Nov 2 15:47:03 2017 UTC (6 years, 6 months ago) by guez
Original Path: trunk/Sources/phylmd/phytrac.f
File size: 13258 byte(s)
Rename phisinit to phis in restart.nc: clearer, same name as Fortran variable.

In aaam_bud, use rlat and rlon from phyetat0_m instead of having these
module variables associated to actual arguments in physiq.

In clmain, too many wind variables make the procedure hard to
understand. Use yu(:knon, 1) and yv(:knon, 1) instead of u1lay(:knon)
and v1lay(:knon). Note that when yu(:knon, 1) and yv(:knon, 1) are
used as actual arguments, they are probably copied to new arrays since
the elements are not contiguous. Rename yu10m to wind10m because this
is the norm of wind vector, not its zonal component. Rename yustar to
ustar. Rename uzon and vmer to u1 and v1 since these are wind
components at first layer and u1 and v1 are the names of corresponding
dummy arguments in stdlevvar.

In clmain, rename yzlev to zlev.

In clmain, screenc, stdlevvar and coefcdrag, remove the code
corresponding to zxli true (not used in LMDZ either).

Subroutine ustarhb becomes a function. Simplifications using the fact
that zx_alf2 = 0 and zx_alf1 = 1 (discarding the possibility to change
this).

In procedure vdif_kcay, remove unused dummy argument plev. Remove
useless computations of sss and sssq.

In clouds_gno, exp(100.) would overflow in single precision. Set
maximum to exp(80.) instead.

In physiq, use u(:, 1) and v(:, 1) as arguments to phytrac instead of
creating ad hoc variables yu1 and yv1.

In stdlevvar, rename dummy argument u_10m to wind10m, following the
corresponding modification in clmain. Simplifications using the fact
that ok_pred = 0 and ok_corr = 1 (discarding the possibility to change
this).

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, fm_therm, entr_therm, yu1, &
12 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 dimens_m, 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(klon, llm) ! coeff melange couche limite
81 real fm_therm(klon, llm+1), entr_therm(klon, llm) ! thermiques
82 REAL, intent(in):: yu1(:), yv1(:) ! (klon) vent au premier niveau
83
84 ! Arguments n\'ecessaires pour les sources et puits de traceur :
85 real, intent(in):: ftsol(:, :) ! (klon, nbsrf) surface temperature (K)
86 real, intent(in):: pctsrf(klon, nbsrf) ! Pourcentage de sol f(nature du sol)
87
88 ! Lessivage pour le on-line
89 REAL, intent(in):: frac_impa(klon, llm) ! fraction d'aerosols impactes
90 REAL, intent(in):: frac_nucl(klon, llm) ! fraction d'aerosols nuclees
91
92 ! Kerry Emanuel
93 real, intent(in):: da(klon, llm), phi(klon, llm, llm), mp(klon, llm)
94 REAL, intent(in):: upwd(klon, llm) ! saturated updraft mass flux
95 REAL, intent(in):: dnwd(klon, llm) ! saturated downdraft mass flux
96
97 real, intent(inout):: tr_seri(:, :, :) ! (klon, llm, nqmx - 2)
98 ! (mass fractions of tracers, excluding water, at mid-layers)
99
100 real, intent(in):: zmasse(:, :) ! (klon, llm)
101 ! (column-density of mass of air in a cell, in kg m-2)
102
103 integer, intent(in):: ncid_startphy
104
105 ! Local:
106
107 integer nsplit
108
109 ! TRACEURS
110
111 ! Sources et puits des traceurs:
112
113 ! Pour l'instant seuls les cas du rn et du pb ont ete envisages.
114
115 REAL source(klon) ! a voir lorsque le flux est prescrit
116 !
117 ! Pour la source de radon et son reservoir de sol
118
119 REAL, save:: trs(klon, nqmx - 2) ! Concentration de traceur dans le sol
120
121 REAL masktr(klon, nqmx - 2) ! Masque reservoir de sol traceur
122 ! Masque de l'echange avec la surface
123 ! (1 = reservoir) ou (possible => 1)
124 SAVE masktr
125 REAL fshtr(klon, nqmx - 2) ! Flux surfacique dans le reservoir de sol
126 SAVE fshtr
127 REAL hsoltr(nqmx - 2) ! Epaisseur equivalente du reservoir de sol
128 SAVE hsoltr
129 REAL tautr(nqmx - 2) ! Constante de decroissance radioactive
130 SAVE tautr
131 REAL vdeptr(nqmx - 2) ! Vitesse de depot sec dans la couche Brownienne
132 SAVE vdeptr
133 REAL scavtr(nqmx - 2) ! Coefficient de lessivage
134 SAVE scavtr
135
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, 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