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

Annotation of /trunk/phylmd/phytrac.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 227 - (hide 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 guez 3 module phytrac_m
2    
3     IMPLICIT none
4    
5     private
6     public phytrac
7    
8     contains
9    
10 guez 202 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 guez 3
15 guez 213 ! 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 guez 3
19 guez 90 ! Authors: Fr\'ed\'eric Hourdin, Abderrahmane Idelkadi, Marie-Alice
20 guez 213 ! Foujols
21 guez 78
22 guez 90 ! Objet : moniteur g\'en\'eral des tendances des traceurs
23 guez 3
24 guez 78 ! L'appel de "phytrac" se fait avec "nqmx - 2" donc nous avons
25 guez 98 ! bien les vrais traceurs, sans la vapeur d'eau ni l'eau liquide.
26 guez 3
27 guez 47 ! Modifications pour les traceurs :
28 guez 213 ! - uniformisation des param\'etrisations dans phytrac
29 guez 90 ! - stockage des moyennes des champs n\'ecessaires en mode traceur off-line
30 guez 47
31 guez 78 use abort_gcm_m, only: abort_gcm
32 guez 182 use clesphys2, only: conv_emanuel
33 guez 120 use cltrac_m, only: cltrac
34 guez 78 use cltracrn_m, only: cltracrn
35 guez 202 USE conf_gcm_m, ONLY: lmt_pas
36 guez 3 use ctherm, only: iflag_thermals
37 guez 120 use cvltr_m, only: cvltr
38 guez 98 use dimens_m, only: llm, nqmx
39     use dimphy, only: klon
40 guez 201 use histwrite_phy_m, only: histwrite_phy
41 guez 78 use indicesol, only: nbsrf
42 guez 201 use iniadvtrac_m, only: tname
43 guez 98 use initrrnpb_m, only: initrrnpb
44 guez 17 use minmaxqfi_m, only: minmaxqfi
45 guez 171 use netcdf, only: NF90_FILL_float
46 guez 157 use netcdf95, only: nf95_inq_varid, nf95_get_var, nf95_put_var
47 guez 78 use nflxtr_m, only: nflxtr
48 guez 36 use nr_util, only: assert
49 guez 78 use o3_chem_m, only: o3_chem
50     use phyetat0_m, only: rlat
51 guez 157 use phyredem0_m, only: ncid_restartphy
52 guez 17 use press_coefoz_m, only: press_coefoz
53 guez 78 use radiornpb_m, only: radiornpb
54     use regr_pr_comb_coefoz_m, only: regr_pr_comb_coefoz
55     use SUPHEC_M, only: rg
56 guez 191 use time_phylmdz, only: itap
57 guez 3
58     integer, intent(in):: julien !jour julien, 1 <= julien <= 360
59 guez 90 real, intent(in):: gmtime ! heure de la journ\'ee en fraction de jour
60 guez 78 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 guez 3
65     real, intent(in):: paprs(klon, llm+1)
66     ! (pression pour chaque inter-couche, en Pa)
67    
68 guez 10 real, intent(in):: pplay(klon, llm)
69     ! (pression pour le mileu de chaque couche, en Pa)
70    
71 guez 78 ! convection:
72 guez 3
73 guez 62 REAL, intent(in):: pmfu(klon, llm) ! flux de masse dans le panache montant
74 guez 72
75     REAL, intent(in):: pmfd(klon, llm)
76     ! flux de masse dans le panache descendant
77    
78 guez 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 guez 225 real fm_therm(klon, llm+1), entr_therm(klon, llm) ! thermiques
82 guez 227 REAL, intent(in):: yu1(:), yv1(:) ! (klon) vent au premier niveau
83 guez 3
84 guez 90 ! Arguments n\'ecessaires pour les sources et puits de traceur :
85 guez 222 real, intent(in):: ftsol(:, :) ! (klon, nbsrf) surface temperature (K)
86 guez 203 real, intent(in):: pctsrf(klon, nbsrf) ! Pourcentage de sol f(nature du sol)
87 guez 3
88 guez 78 ! Lessivage pour le on-line
89 guez 213 REAL, intent(in):: frac_impa(klon, llm) ! fraction d'aerosols impactes
90     REAL, intent(in):: frac_nucl(klon, llm) ! fraction d'aerosols nuclees
91 guez 3
92 guez 78 ! Kerry Emanuel
93 guez 99 real, intent(in):: da(klon, llm), phi(klon, llm, llm), mp(klon, llm)
94 guez 78 REAL, intent(in):: upwd(klon, llm) ! saturated updraft mass flux
95     REAL, intent(in):: dnwd(klon, llm) ! saturated downdraft mass flux
96 guez 3
97 guez 98 real, intent(inout):: tr_seri(:, :, :) ! (klon, llm, nqmx - 2)
98 guez 78 ! (mass fractions of tracers, excluding water, at mid-layers)
99 guez 3
100 guez 78 real, intent(in):: zmasse(:, :) ! (klon, llm)
101 guez 17 ! (column-density of mass of air in a cell, in kg m-2)
102 guez 3
103 guez 191 integer, intent(in):: ncid_startphy
104 guez 3
105 guez 157 ! Local:
106    
107 guez 17 integer nsplit
108    
109 guez 78 ! TRACEURS
110 guez 17
111 guez 3 ! Sources et puits des traceurs:
112    
113     ! Pour l'instant seuls les cas du rn et du pb ont ete envisages.
114    
115 guez 78 REAL source(klon) ! a voir lorsque le flux est prescrit
116 guez 3 !
117     ! Pour la source de radon et son reservoir de sol
118    
119 guez 156 REAL, save:: trs(klon, nqmx - 2) ! Concentration de traceur dans le sol
120 guez 3
121 guez 98 REAL masktr(klon, nqmx - 2) ! Masque reservoir de sol traceur
122 guez 78 ! Masque de l'echange avec la surface
123     ! (1 = reservoir) ou (possible => 1)
124 guez 3 SAVE masktr
125 guez 98 REAL fshtr(klon, nqmx - 2) ! Flux surfacique dans le reservoir de sol
126 guez 3 SAVE fshtr
127 guez 98 REAL hsoltr(nqmx - 2) ! Epaisseur equivalente du reservoir de sol
128 guez 3 SAVE hsoltr
129 guez 98 REAL tautr(nqmx - 2) ! Constante de decroissance radioactive
130 guez 3 SAVE tautr
131 guez 98 REAL vdeptr(nqmx - 2) ! Vitesse de depot sec dans la couche Brownienne
132 guez 3 SAVE vdeptr
133 guez 98 REAL scavtr(nqmx - 2) ! Coefficient de lessivage
134 guez 3 SAVE scavtr
135    
136     CHARACTER itn
137    
138 guez 225 logical, save:: aerosol(nqmx - 2) ! Nature du traceur
139 guez 78 ! ! aerosol(it) = true => aerosol
140     ! ! aerosol(it) = false => gaz
141 guez 3
142 guez 225 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 guez 3 ! 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 guez 98 REAL d_tr_cl(klon, llm, nqmx - 2) ! tendance de traceurs couche limite
156 guez 156
157     REAL d_tr_cv(klon, llm, nqmx - 2)
158     ! tendance de traceurs conv pour chq traceur
159    
160 guez 98 REAL d_tr_th(klon, llm, nqmx - 2) ! la tendance des thermiques
161 guez 3 REAL d_tr_dec(klon, llm, 2) ! la tendance de la decroissance
162 guez 78 ! ! radioactive du rn - > pb
163 guez 3
164 guez 213 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 guez 3 real ztra_th(klon, llm)
173 guez 157 integer isplit, varid
174 guez 3
175     ! Controls:
176     logical:: convection = .true.
177    
178     !--------------------------------------
179    
180 guez 18 call assert(shape(zmasse) == (/klon, llm/), "phytrac zmasse")
181 guez 98 call assert(shape(tr_seri) == (/klon, llm, nqmx - 2/), "phytrac tr_seri")
182 guez 3
183 guez 7 if (firstcal) then
184 guez 3 ! Initialisation de certaines variables pour le radon et le plomb
185     ! Initialisation du traceur dans le sol (couche limite radonique)
186 guez 157 trs(:, 2:) = 0.
187 guez 3
188 guez 157 call nf95_inq_varid(ncid_startphy, "trs", varid)
189     call nf95_get_var(ncid_startphy, varid, trs(:, 1))
190 guez 171 if (any(trs(:, 1) == NF90_FILL_float)) call abort_gcm("phytrac", &
191     "some missing values in trs(:, 1)")
192 guez 3
193     ! Initialisation de la fraction d'aerosols lessivee
194    
195 guez 98 d_tr_lessi_impa = 0.
196     d_tr_lessi_nucl = 0.
197 guez 3
198     ! Initialisation de la nature des traceurs
199    
200 guez 225 aerosol = .FALSE. ! Tous les traceurs sont des gaz par defaut
201     radio = .FALSE. ! par d\'efaut pas de passage par "radiornpb"
202 guez 17
203 guez 98 if (nqmx >= 5) then
204 guez 17 call press_coefoz ! read input pressure levels for ozone coefficients
205     end if
206 guez 3
207 guez 62 ! Initialisation du traceur dans le sol (couche limite radonique)
208 guez 3 radio(1)= .true.
209     radio(2)= .true.
210 guez 225 clsol(:2)= .true.
211     clsol(3:)= .false.
212 guez 3 aerosol(2) = .TRUE. ! le Pb est un aerosol
213 guez 98 call initrrnpb(pctsrf, masktr, fshtr, hsoltr, tautr, vdeptr, scavtr)
214 guez 3 endif
215    
216     if (convection) then
217 guez 62 ! Calcul de l'effet de la convection
218 guez 98 DO it=1, nqmx - 2
219 guez 182 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 guez 78 CALL nflxtr(pdtphys, pmfu, pmfd, pde_u, pen_d, paprs, &
224     tr_seri(:, :, it), d_tr_cv(:, :, it))
225 guez 3 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 guez 78 CALL minmaxqfi(tr_seri(:, :, it), 0., 1e33, &
234 guez 3 'convection, tracer index = ' // itn)
235     ENDDO
236     endif
237    
238     ! Calcul de l'effet des thermiques
239    
240 guez 98 do it=1, nqmx - 2
241 guez 3 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 guez 78 tr_seri(i, k, it)=min(tr_seri(i, k, it), 1e10)
246 guez 3 enddo
247     enddo
248     enddo
249    
250     if (iflag_thermals > 0) then
251     nsplit=10
252 guez 98 DO it=1, nqmx - 2
253 guez 3 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 guez 78 ! Calcul de l'effet de la couche limite
271 guez 3
272 guez 225 DO k = 1, llm
273     DO i = 1, klon
274     delp(i, k) = paprs(i, k)-paprs(i, k+1)
275 guez 3 ENDDO
276 guez 225 ENDDO
277 guez 3
278 guez 225 ! 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 guez 3 ENDDO
290 guez 225 ENDDO
291 guez 3
292 guez 225 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 guez 3
300 guez 225 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 guez 3 ENDDO
306 guez 225 ENDDO
307     endif
308     ENDDO
309 guez 3
310 guez 78 ! Calcul de l'effet du puits radioactif
311 guez 3
312     ! MAF il faudrait faire une modification pour passer dans radiornpb
313 guez 98 ! 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 guez 3
323 guez 98 if (nqmx >= 5) then
324 guez 6 ! Ozone as a tracer:
325 guez 7 if (mod(itap - 1, lmt_pas) == 0) then
326     ! Once per day, update the coefficients for ozone chemistry:
327 guez 162 call regr_pr_comb_coefoz(julien, paprs, pplay)
328 guez 7 end if
329 guez 6 call o3_chem(julien, gmtime, t_seri, zmasse, pdtphys, tr_seri(:, :, 3))
330     end if
331 guez 3
332     ! Calcul de l'effet de la precipitation
333    
334 guez 213 d_tr_lessi_nucl = 0.
335     d_tr_lessi_impa = 0.
336     flestottr = 0.
337 guez 3
338 guez 213 ! tendance des aerosols nuclees et impactes
339 guez 3
340 guez 213 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 guez 3 ENDDO
349 guez 213 ENDDO
350     ENDIF
351     ENDDO
352 guez 3
353 guez 213 ! Mises a jour des traceurs + calcul des flux de lessivage
354     ! Mise a jour due a l'impaction et a la nucleation
355 guez 3
356 guez 213 DO it = 1, nqmx - 2
357     IF (aerosol(it)) THEN
358 guez 3 DO k = 1, llm
359     DO i = 1, klon
360 guez 213 tr_seri(i, k, it) = tr_seri(i, k, it) * frac_impa(i, k) &
361     * frac_nucl(i, k)
362 guez 3 ENDDO
363     ENDDO
364 guez 213 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 guez 3 ENDDO
376 guez 213 ENDDO
377 guez 3
378 guez 78 ! Ecriture des sorties
379 guez 201 CALL histwrite_phy("zmasse", zmasse)
380     DO it=1, nqmx - 2
381     CALL histwrite_phy(tname(it+2), tr_seri(:, :, it))
382 guez 213 CALL histwrite_phy("fl"//tname(it+2), flestottr(:, :, it))
383 guez 201 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 guez 3
388     if (lafin) then
389 guez 157 call nf95_inq_varid(ncid_restartphy, "trs", varid)
390     call nf95_put_var(ncid_restartphy, varid, trs(:, 1))
391 guez 3 endif
392    
393     END SUBROUTINE phytrac
394    
395     end module phytrac_m

  ViewVC Help
Powered by ViewVC 1.1.21