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

Diff of /trunk/Sources/phylmd/cltracrn.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/libf/phylmd/cltracrn.f90 revision 38 by guez, Thu Jan 6 17:52:19 2011 UTC trunk/Sources/phylmd/cltracrn.f revision 225 by guez, Mon Oct 16 12:35:41 2017 UTC
# Line 1  Line 1 
1  SUBROUTINE cltracrn( itr, dtime,u1lay, v1lay, &  module cltracrn_m
      coef,t,ftsol,pctsrf, &  
      tr,trs,paprs,pplay,delp, &  
      masktr,fshtr,hsoltr,tautr,vdeptr, &  
      lat, &  
      d_tr,d_trs )  
   
   ! This procedure is clean: no C preprocessor directive, no include line.  
   ! From phylmd/cltracrn.F,v 1.2 2005/05/25 13:10:09  
   
   use indicesol, only: nbsrf  
   use dimphy, only: klon, klev  
   use SUPHEC_M, only: RD, rg  
2    
3    IMPLICIT none    IMPLICIT none
   !======================================================================  
   ! Auteur(s): Alex/LMD) date:  fev 99  
   !            inspire de clqh + clvent  
   ! Objet: diffusion verticale de traceurs avec quantite de traceur ds  
   !        le sol ( reservoir de sol de radon )  
   !          
   ! note : pour l'instant le traceur dans le sol et le flux sont  
   !        calcules mais ils ne servent que de diagnostiques  
   !        seule la tendance sur le traceur est sortie (d_tr)  
   !======================================================================  
   ! Arguments:  
   ! itr---  -input-R- le type de traceur 1- Rn 2 - Pb  
   ! dtime----input-R- intervalle du temps (en second)  
   ! u1lay----input-R- vent u de la premiere couche (m/s)  
   ! v1lay----input-R- vent v de la premiere couche (m/s)  
   ! coef-----input-R- le coefficient d'echange (m**2/s) l>1  
   ! paprs----input-R- pression a inter-couche (Pa)  
   ! pplay----input-R- pression au milieu de couche (Pa)  
   ! delp-----input-R- epaisseur de couche (Pa)  
   ! ftsol----input-R- temperature du sol (en Kelvin)  
   ! tr-------input-R- traceurs  
   ! trs------input-R- traceurs dans le sol  
   ! masktr---input-R- Masque reservoir de sol traceur (1 = reservoir)  
   ! fshtr----input-R- Flux surfacique de production dans le sol  
   ! tautr----input-R- Constante de decroissance du traceur  
   ! vdeptr---input-R- Vitesse de depot sec dans la couche brownienne  
   ! hsoltr---input-R- Epaisseur equivalente du reservoir de sol  
   ! lat-----input-R- latitude en degree  
   ! d_tr-----output-R- le changement de "tr"  
   ! d_trs----output-R- le changement de "trs"  
   !======================================================================  
   !======================================================================  
   REAL, intent(in):: dtime  
   REAL u1lay(klon), v1lay(klon)  
   REAL coef(klon,klev)  
   REAL, intent(in):: t(klon,klev) ! temperature (K)  
   real ftsol(klon,nbsrf), pctsrf(klon,nbsrf)  
   REAL tr(klon,klev), trs(klon)  
   REAL, intent(in):: paprs(klon,klev+1)  
   real, intent(in):: pplay(klon,klev)  
   real delp(klon,klev)  
   REAL masktr(klon)  
   REAL fshtr(klon)  
   REAL hsoltr  
   REAL tautr  
   REAL vdeptr  
   REAL, intent(in):: lat(klon)    
   REAL d_tr(klon,klev)  
   !======================================================================  
   REAL d_trs(klon)         ! (diagnostic) traceur ds le sol  
   !======================================================================  
   INTEGER i, k, itr, n, l  
   REAL rotrhi(klon)  
   REAL zx_coef(klon,klev)  
   REAL zx_buf(klon)  
   REAL zx_ctr(klon,klev)  
   REAL zx_dtr(klon,klev)  
   REAL zx_trs(klon)  
   REAL zx_a, zx_b  
   
   REAL local_tr(klon,klev)  
   REAL local_trs(klon)  
   REAL zts(klon)  
   REAL zx_alpha1(klon), zx_alpha2(klon)  
   !======================================================================  
   !AA Pour l'instant les 4 types de surface ne sont pas pris en compte  
   !AA On fabrique avec zts un champ de temperature de sol    
   !AA que le pondere par la fraction de nature de sol.  
   !  
   !      print*,'PASSAGE DANS CLTRACRN'  
   
   DO i = 1,klon  
      zts(i) = 0.  
   ENDDO  
   !  
   DO n=1,nbsrf  
      DO i = 1,klon  
         zts(i) = zts(i) + ftsol(i,n)*pctsrf(i,n)  
      ENDDO  
   ENDDO  
   !  
   DO i = 1,klon  
      rotrhi(i) = RD * zts(i) / hsoltr  
   END DO  
   !  
   DO k = 1, klev  
      DO i = 1, klon  
         local_tr(i,k) = tr(i,k)  
      ENDDO  
   ENDDO  
   !  
   DO i = 1, klon  
      local_trs(i) = trs(i)  
   ENDDO  
   !======================================================================  
   !AA   Attention si dans clmain zx_alf1(i) = 1.0  
   !AA   Il doit y avoir coherence (dc la meme chose ici)  
   
   DO i = 1, klon  
      zx_alpha1(i) = 1.0  
      zx_alpha2(i) = 1.0 - zx_alpha1(i)  
   ENDDO  
   !======================================================================  
   DO i = 1, klon  
      zx_coef(i,1) = coef(i,1) &  
           * (1.0+SQRT(u1lay(i)**2+v1lay(i)**2)) &  
           * pplay(i,1)/(RD*t(i,1))  
      zx_coef(i,1) = zx_coef(i,1) * dtime*RG  
   ENDDO  
   !  
   DO k = 2, klev  
      DO i = 1, klon  
         zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k)) &  
              *(paprs(i,k)*2/(t(i,k)+t(i,k-1))/RD)**2  
         zx_coef(i,k) = zx_coef(i,k) * dtime*RG  
      ENDDO  
   ENDDO  
   !======================================================================  
   DO i = 1, klon  
      zx_buf(i)      = delp(i,klev) + zx_coef(i,klev)  
      zx_ctr(i,klev) = local_tr(i,klev)*delp(i,klev)/zx_buf(i)  
      zx_dtr(i,klev) = zx_coef(i,klev) / zx_buf(i)  
   ENDDO  
   !  
   DO l = klev-1, 2 , -1  
      DO i = 1, klon  
         zx_buf(i) = delp(i,l)+zx_coef(i,l) &  
              +zx_coef(i,l+1)*(1.-zx_dtr(i,l+1))  
         zx_ctr(i,l) = ( local_tr(i,l)*delp(i,l) &  
              + zx_coef(i,l+1)*zx_ctr(i,l+1) )/zx_buf(i)  
         zx_dtr(i,l) = zx_coef(i,l) / zx_buf(i)  
      ENDDO  
   ENDDO  
   !  
   DO i = 1, klon  
      zx_buf(i) = delp(i,1) + zx_coef(i,2)*(1.-zx_dtr(i,2)) &  
           + masktr(i) * zx_coef(i,1) &  
           *( zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2) )  
      zx_ctr(i,1) = ( local_tr(i,1)*delp(i,1) &  
           + zx_ctr(i,2) &  
           *(zx_coef(i,2) &  
           - masktr(i) * zx_coef(i,1) &  
           *zx_alpha2(i) ) ) / zx_buf(i)  
      zx_dtr(i,1) = masktr(i) * zx_coef(i,1) / zx_buf(i)  
   ENDDO  
   !======================================================================  
   ! Calculer d'abord local_trs nouvelle quantite dans le reservoir  
   ! de sol  
   !  
   !-------------------------  
   ! Au dessus des continents  
   !-------------------------  
   ! Le pb peut se deposer partout : vdeptr = 10-3 m/s  
   ! Le Rn est traiter commme une couche Brownienne puisque vdeptr = 0.  
   !  
   DO i = 1, klon  
      !  
      IF ( NINT(masktr(i)) .EQ. 1 ) THEN  
         zx_trs(i) = local_trs(i)  
         zx_a = zx_trs(i) &  
              +fshtr(i)*dtime*rotrhi(i) &  
              +rotrhi(i)*masktr(i)*zx_coef(i,1)/RG &  
              *(zx_ctr(i,1)*(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2)) &  
              +zx_alpha2(i)*zx_ctr(i,2))  
         zx_b = 1. + rotrhi(i)*masktr(i)*zx_coef(i,1)/RG &  
              * (1.-zx_dtr(i,1) &  
              *(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2))) &  
              + dtime / tautr &  
              !AA: Pour l'instant, pour aller vite, le depot sec est traite  
              ! comme une decroissance  
              + dtime * vdeptr / hsoltr  
         zx_trs(i) = zx_a / zx_b  
         local_trs(i) = zx_trs(i)  
      ENDIF  
      !  
      ! Si on est entre 60N et 70N on divise par 2 l'emanation  
      !--------------------------------------------------------  
      !  
      IF &  
           ( (itr.eq.1.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GE.60. &  
           .AND.lat(i).LE.70.) &  
           .OR. &  
           (itr.eq.2.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GE.60. &  
           .AND.lat(i).LE.70.) ) &  
           THEN  
         zx_trs(i) = local_trs(i)  
         zx_a = zx_trs(i) &  
              +(fshtr(i)/2.)*dtime*rotrhi(i) &  
              +rotrhi(i)*masktr(i)*zx_coef(i,1)/RG &  
              *(zx_ctr(i,1)*(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2)) &  
              +zx_alpha2(i)*zx_ctr(i,2))  
         zx_b = 1. + rotrhi(i)*masktr(i)*zx_coef(i,1)/RG &  
              * (1.-zx_dtr(i,1) &  
              *(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2))) &  
              + dtime / tautr &  
              + dtime * vdeptr / hsoltr  
         zx_trs(i) = zx_a / zx_b  
         local_trs(i) = zx_trs(i)  
      ENDIF  
      !  
      !----------------------------------------------  
      ! Au dessus des oceans et aux hautes latitudes  
      !----------------------------------------------  
      !  
      ! au dessous de -60S  pas d'emission de radon au dessus  
      ! des oceans et des continents  
      !---------------------------------------------------------------  
   
      IF ( (itr.EQ.1.AND.NINT(masktr(i)).EQ.0) &  
           .OR. &  
           (itr.EQ.1.AND.NINT(masktr(i)).EQ.1.AND.lat(i).LT.-60.)) &  
           THEN  
         zx_trs(i) = 0.  
         local_trs(i) = 0.  
      END IF  
   
      ! au dessus de 70 N pas d'emission de radon au dessus  
      ! des oceans et des continents  
      !--------------------------------------------------------------  
      IF ( (itr.EQ.1.AND.NINT(masktr(i)).EQ.0) &  
           .OR. &  
           (itr.EQ.1.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GT.70.)) &  
           THEN  
         zx_trs(i) = 0.  
         local_trs(i) = 0.  
      END IF  
   
      ! Au dessus des oceans la source est nulle  
      !-----------------------------------------  
      !  
      IF (itr.eq.1.AND.NINT(masktr(i)).EQ.0) THEN  
         zx_trs(i) = 0.  
         local_trs(i) = 0.  
      END IF  
      !  
   ENDDO    ! sur le i=1,klon  
   !  
   !======================================================================  
   !==== une fois on a zx_trs, on peut faire l'iteration ========  
   !  
   DO i = 1, klon  
      local_tr(i,1) = zx_ctr(i,1)+zx_dtr(i,1)*zx_trs(i)  
   ENDDO  
   DO l = 2, klev  
      DO i = 1, klon  
         local_tr(i,l) &  
              = zx_ctr(i,l) + zx_dtr(i,l)*local_tr(i,l-1)  
      ENDDO  
   ENDDO  
   !======================================================================  
   !== Calcul des tendances du traceur ds le sol et dans l'atmosphere  
   !  
   DO l = 1, klev  
      DO i = 1, klon  
         d_tr(i,l) = local_tr(i,l) - tr(i,l)  
      ENDDO  
   ENDDO  
   DO i = 1, klon  
      d_trs(i) = local_trs(i) - trs(i)  
   ENDDO  
4    
5  END SUBROUTINE cltracrn  contains
6    
7      SUBROUTINE cltracrn(itr, dtime, u1lay, v1lay, coef, t, ftsol, pctsrf, tr, &
8           trs, paprs, pplay, delp, masktr, fshtr, hsoltr, tautr, vdeptr, lat, &
9           d_tr, d_trs)
10    
11        ! From phylmd/cltracrn.F, version 1.2 2005/05/25 13:10:09
12    
13        ! Author: Alexandre ARMENGAUD
14        ! Date: February 1999
15        ! Inspiré de clqh et clvent
16    
17        ! Objet: diffusion verticale de traceurs avec quantité de traceur
18        ! dans le sol (réservoir de sol de radon)
19    
20        ! Note : pour l'instant le traceur dans le sol et le flux sont
21        ! calculés mais ils ne servent que de diagnostics. Seule la
22        ! tendance sur le traceur est sortie (d_tr).
23    
24        use indicesol, only: nbsrf
25        use dimphy, only: klon, klev
26        use SUPHEC_M, only: RD, rg
27    
28        INTEGER itr
29        ! itr--- -input-R- le type de traceur 1- Rn 2 - Pb
30    
31        REAL, intent(in):: dtime
32        ! dtime----input-R- intervalle de temps (en second)
33        REAL, intent(in):: u1lay(klon), v1lay(klon) ! vent de la premiere
34                                                    ! couche (m/s)
35        REAL coef(klon, klev)
36        ! coef-----input-R- le coefficient d'echange (m**2/s) l>1
37        REAL, intent(in):: t(klon, klev) ! temperature (K)
38        real, intent(in):: ftsol(klon, nbsrf), pctsrf(klon, nbsrf)
39        ! ftsol----input-R- temperature du sol (en Kelvin)
40        REAL, intent(in):: tr(klon, klev) ! traceur
41        REAL, intent(in):: trs(:) ! (klon) traceur dans le sol
42        REAL, intent(in):: paprs(klon, klev+1)
43        ! paprs----input-R- pression a l'inter-couche (Pa)
44        real, intent(in):: pplay(klon, klev)
45        ! pplay----input-R- pression au milieu de couche (Pa)
46        real delp(klon, klev)
47        ! delp-----input-R- epaisseur de couche (Pa)
48        REAL masktr(klon)
49        ! masktr---input-R- Masque reservoir de sol traceur (1 = reservoir)
50        REAL fshtr(klon)
51        ! fshtr----input-R- Flux surfacique de production dans le sol
52        REAL hsoltr
53        ! hsoltr---input-R- Epaisseur equivalente du reservoir de sol
54        REAL tautr
55        ! tautr----input-R- Constante de decroissance du traceur
56    
57        REAL, intent(in):: vdeptr
58        ! vitesse de d\'ep\^ot sec dans la couche brownienne
59    
60        REAL, intent(in):: lat(klon)
61        ! lat-----input-R- latitude en degree
62        REAL d_tr(klon, klev)
63        ! d_tr-----output-R- le changement de "tr"
64    
65        REAL, intent(out):: d_trs(:) ! (klon) (diagnostic) changement de "trs"
66    
67        ! Local:
68        INTEGER i, k, n, l
69        REAL rotrhi(klon)
70        REAL zx_coef(klon, klev)
71        REAL zx_buf(klon)
72        REAL zx_ctr(klon, klev)
73        REAL zx_dtr(klon, klev)
74        REAL zx_a, zx_b
75    
76        REAL local_tr(klon, klev)
77        REAL local_trs(klon)
78        REAL zts(klon)
79        REAL zx_alpha1(klon), zx_alpha2(klon)
80    
81        !------------------------------------------------------------------
82    
83        ! Pour l'instant les quatre types de surface ne sont pas pris en
84        ! compte. On fabrique avec zts un champ de température de sol que
85        ! l'on pondère par la fraction de sol.
86    
87        DO i = 1, klon
88           zts(i) = 0.
89        ENDDO
90    
91        DO n=1, nbsrf
92           DO i = 1, klon
93              zts(i) = zts(i) + ftsol(i, n)*pctsrf(i, n)
94           ENDDO
95        ENDDO
96    
97        DO i = 1, klon
98           rotrhi(i) = RD * zts(i) / hsoltr
99        END DO
100    
101        DO k = 1, klev
102           DO i = 1, klon
103              local_tr(i, k) = tr(i, k)
104           ENDDO
105        ENDDO
106    
107        local_trs = trs
108    
109        ! Attention si dans clmain zx_alf1(i) = 1.
110        ! Il doit y avoir coherence (donc la meme chose ici)
111    
112        DO i = 1, klon
113           zx_alpha1(i) = 1.0
114           zx_alpha2(i) = 1.0 - zx_alpha1(i)
115        ENDDO
116    
117        DO i = 1, klon
118           zx_coef(i, 1) = coef(i, 1) &
119                * (1.0+SQRT(u1lay(i)**2+v1lay(i)**2)) &
120                * pplay(i, 1)/(RD*t(i, 1))
121           zx_coef(i, 1) = zx_coef(i, 1) * dtime*RG
122        ENDDO
123    
124        DO k = 2, klev
125           DO i = 1, klon
126              zx_coef(i, k) = coef(i, k)*RG/(pplay(i, k-1)-pplay(i, k)) &
127                   *(paprs(i, k)*2/(t(i, k)+t(i, k-1))/RD)**2
128              zx_coef(i, k) = zx_coef(i, k) * dtime*RG
129           ENDDO
130        ENDDO
131    
132        DO i = 1, klon
133           zx_buf(i) = delp(i, klev) + zx_coef(i, klev)
134           zx_ctr(i, klev) = local_tr(i, klev)*delp(i, klev)/zx_buf(i)
135           zx_dtr(i, klev) = zx_coef(i, klev) / zx_buf(i)
136        ENDDO
137    
138        DO l = klev-1, 2 , -1
139           DO i = 1, klon
140              zx_buf(i) = delp(i, l)+zx_coef(i, l) &
141                   +zx_coef(i, l+1)*(1.-zx_dtr(i, l+1))
142              zx_ctr(i, l) = (local_tr(i, l)*delp(i, l) &
143                   + zx_coef(i, l+1)*zx_ctr(i, l+1))/zx_buf(i)
144              zx_dtr(i, l) = zx_coef(i, l) / zx_buf(i)
145           ENDDO
146        ENDDO
147    
148        DO i = 1, klon
149           zx_buf(i) = delp(i, 1) + zx_coef(i, 2)*(1.-zx_dtr(i, 2)) &
150                + masktr(i) * zx_coef(i, 1) &
151                *(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i, 2))
152           zx_ctr(i, 1) = (local_tr(i, 1)*delp(i, 1) &
153                + zx_ctr(i, 2) &
154                *(zx_coef(i, 2) &
155                - masktr(i) * zx_coef(i, 1) &
156                *zx_alpha2(i))) / zx_buf(i)
157           zx_dtr(i, 1) = masktr(i) * zx_coef(i, 1) / zx_buf(i)
158        ENDDO
159    
160        ! Calculer d'abord local_trs nouvelle quantite dans le reservoir
161        ! de sol
162    
163        ! Au dessus des continents
164        ! Le pb peut se deposer partout : vdeptr = 10-3 m/s
165        ! Le Rn est traité commme une couche Brownienne puisque vdeptr = 0.
166    
167        DO i = 1, klon
168           IF (NINT(masktr(i)) .EQ. 1) THEN
169              zx_a = local_trs(i) &
170                   +fshtr(i)*dtime*rotrhi(i) &
171                   +rotrhi(i)*masktr(i)*zx_coef(i, 1)/RG &
172                   *(zx_ctr(i, 1)*(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i, 2)) &
173                   +zx_alpha2(i)*zx_ctr(i, 2))
174              ! Pour l'instant, pour aller vite, le d\'ep\^ot sec est trait\'e
175              ! comme une d\'ecroissance :
176              zx_b = 1. + rotrhi(i)*masktr(i)*zx_coef(i, 1)/RG &
177                   * (1.-zx_dtr(i, 1) &
178                   *(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i, 2))) &
179                   + dtime / tautr &
180                   + dtime * vdeptr / hsoltr
181              local_trs(i) = zx_a / zx_b
182           ENDIF
183    
184           ! Si on est entre 60N et 70N on divise par 2 l'emanation
185    
186           IF ((itr.eq.1.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GE.60. &
187                .AND.lat(i).LE.70.) .OR. &
188                (itr.eq.2.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GE.60. &
189                .AND.lat(i).LE.70.)) THEN
190              zx_a = local_trs(i) &
191                   +(fshtr(i)/2.)*dtime*rotrhi(i) &
192                   +rotrhi(i)*masktr(i)*zx_coef(i, 1)/RG &
193                   *(zx_ctr(i, 1)*(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i, 2)) &
194                   +zx_alpha2(i)*zx_ctr(i, 2))
195              zx_b = 1. + rotrhi(i)*masktr(i)*zx_coef(i, 1)/RG &
196                   * (1.-zx_dtr(i, 1) &
197                   *(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i, 2))) &
198                   + dtime / tautr &
199                   + dtime * vdeptr / hsoltr
200              local_trs(i) = zx_a / zx_b
201           ENDIF
202    
203           ! Au dessus des oceans et aux hautes latitudes
204    
205           ! au dessous de -60S pas d'emission de radon au dessus
206           ! des oceans et des continents
207    
208           IF ((itr.EQ.1.AND.NINT(masktr(i)).EQ.0) .OR. &
209                (itr.EQ.1.AND.NINT(masktr(i)).EQ.1.AND.lat(i).LT.-60.)) THEN
210              local_trs(i) = 0.
211           END IF
212    
213           ! au dessus de 70 N pas d'emission de radon au dessus
214           ! des oceans et des continents
215    
216           IF ((itr.EQ.1.AND.NINT(masktr(i)).EQ.0) .OR. &
217                (itr.EQ.1.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GT.70.)) THEN
218              local_trs(i) = 0.
219           END IF
220    
221           ! Au dessus des oceans la source est nulle
222    
223           IF (itr.eq.1.AND.NINT(masktr(i)).EQ.0) THEN
224              local_trs(i) = 0.
225           END IF
226        ENDDO ! sur le i=1, klon
227    
228        ! une fois qu'on a local_trs, on peut faire l'iteration
229    
230        DO i = 1, klon
231           local_tr(i, 1) = zx_ctr(i, 1)+zx_dtr(i, 1)*local_trs(i)
232        ENDDO
233        DO l = 2, klev
234           DO i = 1, klon
235              local_tr(i, l) &
236                   = zx_ctr(i, l) + zx_dtr(i, l)*local_tr(i, l-1)
237           ENDDO
238        ENDDO
239    
240        ! Calcul des tendances du traceur dans le sol et dans l'atmosphere
241    
242        DO l = 1, klev
243           DO i = 1, klon
244              d_tr(i, l) = local_tr(i, l) - tr(i, l)
245           ENDDO
246        ENDDO
247        d_trs = local_trs - trs
248    
249      END SUBROUTINE cltracrn
250    
251    end module cltracrn_m
Legend:
Removed from v.38  
changed lines
  Added in v.225
  ViewVC Help
Powered by ViewVC 1.1.21