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

Annotation of /trunk/Sources/phylmd/phytrac.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 191 - (hide annotations)
Mon May 9 19:56:28 2016 UTC (8 years ago) by guez
File size: 14879 byte(s)
Extracted the call to read_comdissnew out of conf_gcm.

Made ok_instan a variable of module clesphys, itau_phy a variable of
module phyetat0_m, nid_ins a variable of module ini_histins_m, itap a
variable of new module time_phylmdz, so that histwrite_phy can be
called from any procedure without the need to cascade those variables
into that procedure. Made itau_w a variable of module time_phylmdz so
that it is computed only once per time step of physics.

Extracted variables of module clesphys which were in namelist
conf_phys_nml into their own namelist, clesphys_nml, and created
procedure read_clesphys reading clesphys_nml, to avoid side effect.

No need for double precision in procedure getso4fromfile. Assume there
is a single variable for the whole year in the NetCDF file instead of
one variable per month.

Created generic procedure histwrite_phy and removed procedure
write_histins, following LMDZ. histwrite_phy has only two arguments,
can be called from anywhere, and should manage the logic of writing or
not writing into various history files with various operations. So the
test on ok_instan goes inside histwrite_phy.

Test for raz_date in phyetat0 instead of physiq to avoid side effect.

Created procedure increment_itap to avoid side effect.

Removed unnecessary differences between procedures readsulfate and
readsulfate_pi.

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

  ViewVC Help
Powered by ViewVC 1.1.21