source: LMDZ.3.3/trunk/libf/phylmd/phytrac.F @ 2

Last change on this file since 2 was 2, checked in by lmdz, 25 years ago

Initial revision

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 17.4 KB
Line 
1      SUBROUTINE phytrac (rnpb,
2     I                   debutphy,
3     I                   nqmax,
4     I                   nlon,nlev,pdtphys,
5     I                   t_seri,paprs,pplay,
6     I                   pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
7     I                   coefh,yu1,yv1,ftsol,pctsrf,xlat,
8     I                   frac_impa, frac_nucl,
9     I                   xlon,presnivs,paire,pphis,
10     O                   tr_seri)
11      USE ioipsl
12      IMPLICIT none
13c======================================================================
14c Auteur(s) FH
15c Objet: Moniteur general des tendances traceurs
16c
17cAA Remarques en vrac:
18cAA--------------------
19cAA 1/ le call phytrac se fait avec nqmax-2 donc nous avons bien
20cAA les vrais traceurs (nbtr) dans phytrac (pas la vapeur ni eau liquide)
21cAA 2/ Le choix du radon et du pb se fait juste avec un data
22cAA    (peu propre). Peut-etre pourrait-on prevoir dans l'avenir
23cAA    une variable "type de traceur"
24c======================================================================
25#include "YOMCST.h"
26#include "dimensions.h"
27#include "dimphy.h"
28#include "indicesol.h"
29#include "control.h"
30c======================================================================
31
32c Arguments:
33c
34c   EN ENTREE:
35c   ==========
36c
37c   divers:
38c   -------
39c
40      integer nlon  ! nombre de points horizontaux
41      integer nlev  ! nombre de couches verticales
42      integer nqmax ! nombre de traceurs auxquels on applique la physique
43      real pdtphys  ! pas d'integration pour la physique (seconde)
44      real t_seri(nlon,nlev) ! temperature
45      real tr_seri(nlon,nlev,nbtr) ! traceur 
46      real paprs(nlon,nlev+1)  ! pression pour chaque inter-couche (en Pa)
47      real pplay(nlon,nlev)  ! pression pour le mileu de chaque couche (en Pa)
48      real presnivs(klev) ! pressions approximat. des milieux couches ( en PA)
49      real znivsig(klev) ! niveaux sigma
50      real paire(klon)
51      real pphis(klon)
52      logical debutphy       ! le flag de l'initialisation de la physique
53      integer ll
54c
55cAA Rem : nbtr : nombre de vrais traceurs est defini dans dimphy.h
56c
57c   convection:
58c   -----------
59c
60      REAL pmfu(nlon,nlev)  ! flux de masse dans le panache montant
61      REAL pmfd(nlon,nlev)  ! flux de masse dans le panache descendant
62      REAL pen_u(nlon,nlev) ! flux entraine dans le panache montant
63      REAL pde_u(nlon,nlev) ! flux detraine dans le panache montant
64      REAL pen_d(nlon,nlev) ! flux entraine dans le panache descendant
65      REAL pde_d(nlon,nlev) ! flux detraine dans le panache descendant
66c
67c   Couche limite:
68c   --------------
69c
70      REAL coefh(nlon,nlev) ! coeff melange CL
71      REAL yu1(nlon)        ! vents au premier niveau
72      REAL yv1(nlon)        ! vents au premier niveau
73      REAL xlat(nlon)       ! latitudes pour chaque point
74      REAL xlon(nlon)       ! longitudes pour chaque point
75
76c
77c   Lessivage:
78c   ----------
79c
80c pour le ON-LINE
81c
82      REAL frac_impa(nlon,nlev)  ! fraction d'aerosols impactes
83      REAL frac_nucl(nlon,nlev)  ! fraction d'aerosols nuclees
84c
85cAA
86cAA Arguments necessaires pour les sources et puits de traceur:
87cAA ----------------
88cAA
89      real ftsol(nlon,nbsrf)  ! Temperature du sol (surf)(Kelvin)
90      real pctsrf(nlon,nbsrf) ! Pourcentage de sol f(nature du sol)
91
92cAA ----------------------------
93cAA  VARIABLES LOCALES TRACEURS
94cAA ----------------------------
95cAA
96cAA Sources et puits des traceurs:
97cAA ------------------------------
98cAA
99cAA Pour l'instant seuls les cas du rn et du pb ont ete envisages.
100
101      REAL source(klon)       ! a voir lorsque le flux est prescrit
102      REAL topuits(klev,nbtr)  ! a voir
103cAA
104cAA Pour la source de radon et son reservoir de sol
105cAA ................................................
106 
107      REAL trs(klon,nbtr)    ! Conc. radon ds le sol
108      SAVE trs
109
110      REAL masktr(klon,nbtr) ! Masque reservoir de sol traceur
111c                            Masque de l'echange avec la surface
112c                           (1 = reservoir) ou (possible => 1 )
113      SAVE masktr
114      REAL fshtr(klon,nbtr)  ! Flux surfacique dans le reservoir de sol
115      SAVE fshtr
116      REAL hsoltr(nbtr)      ! Epaisseur equivalente du reservoir de sol
117      SAVE hsoltr
118      REAL tautr(nbtr)       ! Constante de decroissance radioactive
119      SAVE tautr
120      REAL vdeptr(nbtr)      ! Vitesse de depot sec dans la couche Brownienne
121      SAVE vdeptr
122      REAL scavtr(nbtr)      ! Coefficient de lessivage
123      SAVE scavtr
124cAA
125      CHARACTER*1 itn
126C maf ioipsl
127      CHARACTER*2 str2
128      INTEGER nhori, nvert
129      REAL zsto, zout, zjulian
130      INTEGER nid_tra
131      SAVE nid_tra
132      INTEGER ndex(1)
133      REAL zx_tmp_2d(iim,jjm+1), zx_tmp_3d(iim,jjm+1,klev)
134      REAL zx_lon(iim,jjm+1), zx_lat(iim,jjm+1)
135      INTEGER itra
136      SAVE itra                   ! compteur pour la physique
137C
138C Variables liees a l'ecriture de la bande histoire : phytrac.nc
139c
140      INTEGER ecrit_tra
141      SAVE ecrit_tra   
142C
143C nature du traceur
144c
145      logical aerosol(nbtr)  ! Nature du traceur
146c                            ! aerosol(it) = true  => aerosol
147c                            ! aerosol(it) = false => gaz
148c                            ! nat_trac(it) = 1. aerosolc
149      logical clsol(nbtr)    ! clsol(it) = true => CL sol calculee
150      logical radio(nbtr)    ! radio(it)=true => decroisssance radioactive
151      save aerosol,clsol,radio
152C
153c======================================================================
154c
155c Declaration des procedures appelees
156c
157c--modif convection tiedtke
158      INTEGER i, k, it
159
160      REAL delp(klon,klev)
161c--end modif
162c
163c Variables liees a l'ecriture de la bande histoire physique
164c
165c Variables locales pour effectuer les appels en serie
166c----------------------------------------------------
167c
168      REAL d_tr(klon,klev), d_trs(klon) ! tendances de traceurs
169      REAL d_tr_cl(klon,klev) ! tendance de traceurs  couche limite
170      REAL d_tr_cv(klon,klev) ! tendance de traceurs  convection
171      REAL d_tr_dec(klon,klev,nbtr) ! la tendance de la decroissance
172c                                   ! radioactive du rn - > pb
173      REAL d_tr_lessi_impa(klon,klev,nbtr) ! la tendance du lessivage
174c                                          ! par impaction
175      REAL d_tr_lessi_nucl(klon,klev,nbtr) ! la tendance du lessivage
176c                                          ! par nucleation
177      REAL fluxrn(klon,klev)
178      REAL fluxpb(klon,klev)
179      REAL pbimpa(klon,klev)
180      REAL pbnucl(klon,klev)
181      REAL rn(klon,klev)
182      REAL pb(klon,klev)
183      REAL flestottr(klon,klev,nbtr) ! flux de lessivage
184c                                    ! dans chaque couche
185
186C
187      character*20 modname
188      character*80 abort_message
189c
190c   Controles
191c-------------
192      logical first,couchelimite,convection,lessivage,sorties,
193     s        rnpb,inirnpb
194      save first,couchelimite,convection,lessivage,sorties,
195     s     inirnpb
196      data first,couchelimite,convection,lessivage,sorties
197     s     /.true.,.true.,.true.,.true.,.true./
198c
199c======================================================================
200c         ecrit_tra = NINT(86400./pdtphys *1.0)  ! tous les jours
201         modname='phytrac'
202
203         print*,'DANS PHYTRAC debutphy=',debutphy
204
205         ecrit_tra = NINT(86400./pdtphys *ecritphy)   
206         if (debutphy) then
207
208         if(nbtr.lt.nqmax) then
209           print*,'NQMAX=',nqmax
210           print*,'NBTR=',nbtr
211           abort_message='See above'
212           call abort_gcm(modname,abort_message,1)
213         endif
214
215         inirnpb=rnpb
216         PRINT*, 'La frequence de sortie traceurs est  ', ecrit_tra
217         itra=0
218C         
219         CALL ymds2ju(anneeref, 1, 1, 0.0, zjulian)
220         zjulian = zjulian + dayref
221c
222         CALL gr_fi_ecrit(1,klon,iim,jjm+1,xlon,zx_lon)
223         DO i = 1, iim
224            zx_lon(i,1) = xlon(i+1)
225            zx_lon(i,jjm+1) = xlon(i+1)
226         ENDDO
227         DO ll=1,klev
228            znivsig(ll)=float(ll)
229         ENDDO
230         CALL gr_fi_ecrit(1,klon,iim,jjm+1,xlat,zx_lat)
231         CALL histbeg("histrac", iim,zx_lon, jjm+1,zx_lat,
232     .                 1,iim,1,jjm+1, 0, zjulian, pdtphys,
233     .                 nhori, nid_tra)
234         call histvert(nid_tra, 'sig_s', 'Niveaux sigma','-',
235     .              klev, znivsig, nvert)
236C
237C         CALL histvert(nid_tra, "presnivs", "Vertical levels", "mb",
238C     .                 klev, presnivs, nvert)
239         zsto = pdtphys
240         zout = pdtphys * FLOAT(ecrit_tra)
241c
242         CALL histdef(nid_tra, "phis", "Surface geop. height", "-",
243     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
244     .                "once",  zsto,zout)
245c
246         CALL histdef(nid_tra, "aire", "Grid area", "-",
247     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
248     .                "once",  zsto,zout)
249c
250         DO it=1,nqmax
251         IF (it.LE.99) THEN
252         WRITE(str2,'(i2.2)') it
253C champ 3D
254         CALL histdef(nid_tra, "tr"//str2, "Tracer No."//str2, "U/kga",
255     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
256     .                "ave(X)", zsto,zout)
257         IF (lessivage) THEN
258         CALL histdef(nid_tra, "fl"//str2, "Tracer Flux No."//str2,
259     .                "U/m2/s",iim,jjm+1,nhori, klev,1,klev,nvert, 32,
260     .                "ave(X)", zsto,zout)
261         ENDIF
262         ELSE
263         PRINT*, "Trop de traceurs"
264         CALL abort
265         ENDIF
266         ENDDO
267         CALL histend(nid_tra)
268         ndex(1) = 0
269c
270         i = NINT(zout/zsto)
271         CALL gr_fi_ecrit(1,klon,iim,jjm+1,pphis,zx_tmp_2d)
272         CALL histwrite(nid_tra,"phis",i,zx_tmp_2d,iim*(jjm+1),ndex)
273C
274         i = NINT(zout/zsto)
275         CALL gr_fi_ecrit(1,klon,iim,jjm+1,paire,zx_tmp_2d)
276         CALL histwrite(nid_tra,"aire",i,zx_tmp_2d,iim*(jjm+1),ndex)
277c
278c======================================================================
279c   Initialisation de certaines variables pour le Rn et le Pb
280c======================================================================
281
282c Initialisation du traceur dans le sol (couche limite radonique)
283c
284         print*,'valeur de debut dans phytrac :',debutphy
285         do it=1,nqmax
286            do i=1,klon
287               trs(i,it) = 0.
288            enddo
289         END DO
290c Initialisation de la fraction d'aerosols lessivee
291c
292         DO it=1,nqmax
293           DO k = 1, nlev
294              DO i = 1, klon
295                 d_tr_lessi_impa(i,k,it) =  0.0
296                 d_tr_lessi_nucl(i,k,it) =  0.0
297             ENDDO
298           ENDDO
299         ENDDO
300c
301c Initialisation de la nature des traceurs
302c
303         DO it = 1, nqmax
304            aerosol(it) = .FALSE.  ! Tous les traceurs sont des gaz par defaut
305            radio(it) = .FALSE.    ! Par defaut pas de passage par radiornpb
306            clsol(it) = .FALSE.  ! Par defaut couche limite avec flux prescrit
307         ENDDO
308c
309      ENDIF  ! fin debutphy
310c Initialisation du traceur dans le sol (couche limite radonique)
311      if(inirnpb) THEN
312c
313         radio(1)= .true.
314         radio(2)= .true.
315         clsol(1)= .true.
316         clsol(2)= .true.
317         aerosol(2) = .TRUE. ! le Pb est un aerosol
318c
319         call initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr
320     .                   ,vdeptr,scavtr)
321         inirnpb=.false.
322      endif
323c======================================================================
324c   Calcul de l'effet de la convection
325c======================================================================
326
327      if (convection) then
328
329      print*,'Pas de temps dans phytrac : ',pdtphys
330      DO it=1, nqmax
331      CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
332     .            pplay, paprs, tr_seri(1,1,it), d_tr_cv)
333      DO k = 1, nlev
334      DO i = 1, klon
335         tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cv(i,k)
336      ENDDO
337      ENDDO
338      WRITE(itn,'(i1)') it
339      CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it='//itn)
340      ENDDO
341      print*,'apres nflxtr'
342
343
344      endif ! convection
345
346c======================================================================
347c   Calcul de l'effet de la couche limite
348c======================================================================
349
350      print*,'avant couchelimite'
351      if (couchelimite) then
352
353      DO k = 1, nlev
354      DO i = 1, klon
355         delp(i,k) = paprs(i,k)-paprs(i,k+1)
356      ENDDO
357      ENDDO
358
359C maf modif pour tenir compte du cas rnpb + traceur
360      DO it=1, nqmax
361      print *,'it',it,clsol(it)
362      if (clsol(it)) then  ! couche limite avec quantite dans le sol calculee
363          CALL cltracrn(it, pdtphys, yu1, yv1,
364     e                    coefh,t_seri,ftsol,pctsrf,
365     e                    tr_seri(1,1,it),trs(1,it),
366     e                    paprs, pplay, delp,
367     e                    masktr(1,it),fshtr(1,it),hsoltr(it),
368     e                    tautr(it),vdeptr(it),
369     e                    xlat,
370     s                    d_tr_cl,d_trs)
371          DO k = 1, nlev
372            DO i = 1, klon
373              tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k)
374            ENDDO
375          ENDDO
376c
377c Traceur ds sol
378c
379          DO i = 1, klon
380            trs(i,it) = trs(i,it) + d_trs(i)
381          END DO
382C
383C maf provisoire suppression des prints
384C         WRITE(itn,'(i1)') it
385C         CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracrn it='//itn)
386      else ! couche limite avec flux prescrit
387Cmaf provisoire source / traceur a creer
388        DO i=1, klon
389          source(i) = 0.0 ! pas de source, pour l'instant
390        ENDDO
391C
392          CALL cltrac(pdtphys, coefh,t_seri,
393     s               tr_seri(1,1,it), source,
394     e               paprs, pplay, delp,
395     s               d_tr )
396            DO k = 1, nlev
397               DO i = 1, klon
398                  tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr(i,k)
399               ENDDO
400            ENDDO
401Cmaf provisoire suppression des prints
402Cmaf          WRITE(itn,'(i1)') it
403cmaf          CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracn it='//itn)
404      endif
405      ENDDO
406c
407      endif ! couche limite
408
409      print*,'apres couchelimite'
410
411c======================================================================
412c   Calcul de l'effet du puits radioactif
413c======================================================================
414
415C MAF il faudrait faire une modification pour passer dans radiornpb
416c si radio=true mais pour l'instant radiornpb propre au cas rnpb
417      if(rnpb) then
418        call radiornpb (tr_seri,pdtphys,tautr,d_tr_dec)
419C
420        DO it=1,nqmax
421            if(radio(it)) then
422            DO k = 1, nlev
423               DO i = 1, klon
424                  tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it)
425               ENDDO
426            ENDDO
427            WRITE(itn,'(i1)') it
428            CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'puits rn it='//itn)
429            endif
430        ENDDO
431c
432      endif ! rnpb decroissance  radioactive
433C
434c======================================================================
435c   Calcul de l'effet de la precipitation
436c======================================================================
437
438      IF (lessivage) THEN
439
440      print*,'avant lessivage'
441
442       DO it = 1, nqmax
443           DO k = 1, nlev
444              DO i = 1, klon
445               d_tr_lessi_nucl(i,k,it) = 0.0
446               d_tr_lessi_impa(i,k,it) = 0.0
447               flestottr(i,k,it) = 0.0
448              ENDDO
449           ENDDO
450       ENDDO
451c
452c tendance des aerosols nuclees et impactes
453c
454       DO it = 1, nqmax
455         IF (aerosol(it)) THEN
456           DO k = 1, nlev
457              DO i = 1, klon
458               d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) +
459     s                  ( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
460               d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) +
461     s                  ( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
462              ENDDO
463           ENDDO
464         ENDIF
465       ENDDO
466c
467c Mises a jour des traceurs + calcul des flux de lessivage
468c Mise a jour due a l'impaction et a la nucleation
469c
470       DO it = 1, nqmax
471         IF (aerosol(it)) THEN
472           DO k = 1, nlev
473              DO i = 1, klon
474               tr_seri(i,k,it) = tr_seri(i,k,it) *
475     s              ( frac_impa(i,k) + frac_nucl(i,k) - 1. )   
476              ENDDO
477           ENDDO
478         ENDIF
479       ENDDO
480c
481c Flux lessivage total
482c
483      DO it = 1, nqmax
484           DO k = 1, nlev
485            DO i = 1, klon
486               flestottr(i,k,it) = flestottr(i,k,it) -
487     s                   ( d_tr_lessi_nucl(i,k,it)   +
488     s                     d_tr_lessi_impa(i,k,it) ) *
489     s                   ( paprs(i,k)-paprs(i,k+1) ) /
490     s                   (RG * pdtphys)
491            ENDDO
492           ENDDO
493c
494Cmaf suppression provisoire
495Cmaf        WRITE(itn,'(i1)') it
496Cmaf    CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'tr(lessi) it='//itn)
497      ENDDO
498c
499      print*,'apres lessivage'
500      ENDIF
501Cc
502      DO k = 1, nlev
503         DO i = 1, klon
504            fluxrn(i,k) = flestottr(i,k,1)
505            fluxpb(i,k) = flestottr(i,k,2)
506            rn(i,k) = tr_seri(i,k,1)
507            pb(i,k) = tr_seri(i,k,2)
508            pbnucl(i,k)=d_tr_lessi_nucl(i,k,2)
509            pbimpa(i,k)=d_tr_lessi_impa(i,k,2)
510         ENDDO
511      ENDDO
512      itra=itra+1
513      ndex(1) = 0
514      DO it=1,nqmax
515      IF (it.LE.99) THEN
516       WRITE(str2,'(i2.2)') it
517C champs 3D
518       CALL gr_fi_ecrit(klev,klon,iim,jjm+1,tr_seri(1,1,it),zx_tmp_3d)
519       CALL histwrite(nid_tra,"tr"//str2,itra,zx_tmp_3d,
520     .                                   iim*(jjm+1)*klev,ndex)
521       IF (lessivage) THEN
522       CALL gr_fi_ecrit(klev,klon,iim,jjm+1,flestottr(1,1,it),zx_tmp_3d)
523       CALL histwrite(nid_tra,"fl"//str2,itra,zx_tmp_3d,
524     .                                   iim*(jjm+1)*klev,ndex)
525      ENDIF
526      ELSE
527         PRINT*, "Trop de traceurs"
528         CALL abort
529      ENDIF
530      ENDDO
531
532      RETURN
533      END
Note: See TracBrowser for help on using the repository browser.