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

Last change on this file since 92 was 53, checked in by lmdzadmin, 25 years ago

Correction sur le calendrier. LF

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