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

Annotation of /trunk/phylmd/phytrac.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 61 - (hide annotations)
Fri Apr 20 14:58:43 2012 UTC (12 years, 1 month ago) by guez
Original Path: trunk/libf/phylmd/phytrac.f90
File size: 17381 byte(s)
No more included file in LMDZE, not even "netcdf.inc".

Created a variable containing the list of common source files in
GNUmakefile. So we now also see clearly files that are specific to
each program.

Split module "histcom". Assembled resulting files in directory
"Histcom".

Removed aliasing in calls to "laplacien".

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

  ViewVC Help
Powered by ViewVC 1.1.21