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

Diff of /trunk/phylmd/cltracrn.f

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

trunk/libf/phylmd/cltracrn.f90 revision 10 by guez, Fri Apr 18 14:45:53 2008 UTC trunk/phylmd/cltracrn.f revision 254 by guez, Mon Feb 5 10:39:38 2018 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 YOMCST, 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, cdragh, t, ftsol, &
8           pctsrf, tr, trs, paprs, pplay, delp, masktr, fshtr, hsoltr, tautr, &
9           vdeptr, lat, 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(:, 2:) ! (klon, 2:klev)
36        ! coef-----input-R- le coefficient d'echange (m**2/s) l>1
37        real cdragh(:) ! klon
38        REAL, intent(in):: t(klon, klev) ! temperature (K)
39        real, intent(in):: ftsol(klon, nbsrf), pctsrf(klon, nbsrf)
40        ! ftsol----input-R- temperature du sol (en Kelvin)
41        REAL, intent(in):: tr(klon, klev) ! traceur
42        REAL, intent(in):: trs(:) ! (klon) traceur dans le sol
43        REAL, intent(in):: paprs(klon, klev+1)
44        ! paprs----input-R- pression a l'inter-couche (Pa)
45        real, intent(in):: pplay(klon, klev)
46        ! pplay----input-R- pression au milieu de couche (Pa)
47        real delp(klon, klev)
48        ! delp-----input-R- epaisseur de couche (Pa)
49        REAL masktr(klon)
50        ! masktr---input-R- Masque reservoir de sol traceur (1 = reservoir)
51        REAL fshtr(klon)
52        ! fshtr----input-R- Flux surfacique de production dans le sol
53        REAL hsoltr
54        ! hsoltr---input-R- Epaisseur equivalente du reservoir de sol
55        REAL tautr
56        ! tautr----input-R- Constante de decroissance du traceur
57    
58        REAL, intent(in):: vdeptr
59        ! vitesse de d\'ep\^ot sec dans la couche brownienne
60    
61        REAL, intent(in):: lat(klon)
62        ! lat-----input-R- latitude en degree
63        REAL d_tr(klon, klev)
64        ! d_tr-----output-R- le changement de "tr"
65    
66        REAL, intent(out):: d_trs(:) ! (klon) (diagnostic) changement de "trs"
67    
68        ! Local:
69        INTEGER i, k, n, l
70        REAL rotrhi(klon)
71        REAL zx_coef(klon, klev)
72        REAL zx_buf(klon)
73        REAL zx_ctr(klon, klev)
74        REAL zx_dtr(klon, klev)
75        REAL zx_a, zx_b
76    
77        REAL local_tr(klon, klev)
78        REAL local_trs(klon)
79        REAL zts(klon)
80        REAL zx_alpha1(klon), zx_alpha2(klon)
81    
82        !------------------------------------------------------------------
83    
84        ! Pour l'instant les quatre types de surface ne sont pas pris en
85        ! compte. On fabrique avec zts un champ de température de sol que
86        ! l'on pondère par la fraction de sol.
87    
88        DO i = 1, klon
89           zts(i) = 0.
90        ENDDO
91    
92        DO n=1, nbsrf
93           DO i = 1, klon
94              zts(i) = zts(i) + ftsol(i, n)*pctsrf(i, n)
95           ENDDO
96        ENDDO
97    
98        DO i = 1, klon
99           rotrhi(i) = RD * zts(i) / hsoltr
100        END DO
101    
102        DO k = 1, klev
103           DO i = 1, klon
104              local_tr(i, k) = tr(i, k)
105           ENDDO
106        ENDDO
107    
108        local_trs = trs
109    
110        ! Attention si dans clmain zx_alf1(i) = 1.
111        ! Il doit y avoir coherence (donc la meme chose ici)
112    
113        DO i = 1, klon
114           zx_alpha1(i) = 1.0
115           zx_alpha2(i) = 1.0 - zx_alpha1(i)
116        ENDDO
117    
118        DO i = 1, klon
119           zx_coef(i, 1) = cdragh(i) &
120                * (1.0+SQRT(u1lay(i)**2+v1lay(i)**2)) &
121                * pplay(i, 1)/(RD*t(i, 1))
122           zx_coef(i, 1) = zx_coef(i, 1) * dtime*RG
123        ENDDO
124    
125        DO k = 2, klev
126           DO i = 1, klon
127              zx_coef(i, k) = coef(i, k)*RG/(pplay(i, k-1)-pplay(i, k)) &
128                   *(paprs(i, k)*2/(t(i, k)+t(i, k-1))/RD)**2
129              zx_coef(i, k) = zx_coef(i, k) * dtime*RG
130           ENDDO
131        ENDDO
132    
133        DO i = 1, klon
134           zx_buf(i) = delp(i, klev) + zx_coef(i, klev)
135           zx_ctr(i, klev) = local_tr(i, klev)*delp(i, klev)/zx_buf(i)
136           zx_dtr(i, klev) = zx_coef(i, klev) / zx_buf(i)
137        ENDDO
138    
139        DO l = klev-1, 2 , -1
140           DO i = 1, klon
141              zx_buf(i) = delp(i, l)+zx_coef(i, l) &
142                   +zx_coef(i, l+1)*(1.-zx_dtr(i, l+1))
143              zx_ctr(i, l) = (local_tr(i, l)*delp(i, l) &
144                   + zx_coef(i, l+1)*zx_ctr(i, l+1))/zx_buf(i)
145              zx_dtr(i, l) = zx_coef(i, l) / zx_buf(i)
146           ENDDO
147        ENDDO
148    
149        DO i = 1, klon
150           zx_buf(i) = delp(i, 1) + zx_coef(i, 2)*(1.-zx_dtr(i, 2)) &
151                + masktr(i) * zx_coef(i, 1) &
152                *(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i, 2))
153           zx_ctr(i, 1) = (local_tr(i, 1)*delp(i, 1) &
154                + zx_ctr(i, 2) &
155                *(zx_coef(i, 2) &
156                - masktr(i) * zx_coef(i, 1) &
157                *zx_alpha2(i))) / zx_buf(i)
158           zx_dtr(i, 1) = masktr(i) * zx_coef(i, 1) / zx_buf(i)
159        ENDDO
160    
161        ! Calculer d'abord local_trs nouvelle quantite dans le reservoir
162        ! de sol
163    
164        ! Au dessus des continents
165        ! Le pb peut se deposer partout : vdeptr = 10-3 m/s
166        ! Le Rn est traité commme une couche Brownienne puisque vdeptr = 0.
167    
168        DO i = 1, klon
169           IF (NINT(masktr(i)) .EQ. 1) THEN
170              zx_a = local_trs(i) &
171                   +fshtr(i)*dtime*rotrhi(i) &
172                   +rotrhi(i)*masktr(i)*zx_coef(i, 1)/RG &
173                   *(zx_ctr(i, 1)*(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i, 2)) &
174                   +zx_alpha2(i)*zx_ctr(i, 2))
175              ! Pour l'instant, pour aller vite, le d\'ep\^ot sec est trait\'e
176              ! comme une d\'ecroissance :
177              zx_b = 1. + rotrhi(i)*masktr(i)*zx_coef(i, 1)/RG &
178                   * (1.-zx_dtr(i, 1) &
179                   *(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i, 2))) &
180                   + dtime / tautr &
181                   + dtime * vdeptr / hsoltr
182              local_trs(i) = zx_a / zx_b
183           ENDIF
184    
185           ! Si on est entre 60N et 70N on divise par 2 l'emanation
186    
187           IF ((itr.eq.1.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GE.60. &
188                .AND.lat(i).LE.70.) .OR. &
189                (itr.eq.2.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GE.60. &
190                .AND.lat(i).LE.70.)) THEN
191              zx_a = local_trs(i) &
192                   +(fshtr(i)/2.)*dtime*rotrhi(i) &
193                   +rotrhi(i)*masktr(i)*zx_coef(i, 1)/RG &
194                   *(zx_ctr(i, 1)*(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i, 2)) &
195                   +zx_alpha2(i)*zx_ctr(i, 2))
196              zx_b = 1. + rotrhi(i)*masktr(i)*zx_coef(i, 1)/RG &
197                   * (1.-zx_dtr(i, 1) &
198                   *(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i, 2))) &
199                   + dtime / tautr &
200                   + dtime * vdeptr / hsoltr
201              local_trs(i) = zx_a / zx_b
202           ENDIF
203    
204           ! Au dessus des oceans et aux hautes latitudes
205    
206           ! au dessous de -60S pas d'emission de radon au dessus
207           ! des oceans et des continents
208    
209           IF ((itr.EQ.1.AND.NINT(masktr(i)).EQ.0) .OR. &
210                (itr.EQ.1.AND.NINT(masktr(i)).EQ.1.AND.lat(i).LT.-60.)) THEN
211              local_trs(i) = 0.
212           END IF
213    
214           ! au dessus de 70 N pas d'emission de radon au dessus
215           ! des oceans et des continents
216    
217           IF ((itr.EQ.1.AND.NINT(masktr(i)).EQ.0) .OR. &
218                (itr.EQ.1.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GT.70.)) THEN
219              local_trs(i) = 0.
220           END IF
221    
222           ! Au dessus des oceans la source est nulle
223    
224           IF (itr.eq.1.AND.NINT(masktr(i)).EQ.0) THEN
225              local_trs(i) = 0.
226           END IF
227        ENDDO ! sur le i=1, klon
228    
229        ! une fois qu'on a local_trs, on peut faire l'iteration
230    
231        DO i = 1, klon
232           local_tr(i, 1) = zx_ctr(i, 1)+zx_dtr(i, 1)*local_trs(i)
233        ENDDO
234        DO l = 2, klev
235           DO i = 1, klon
236              local_tr(i, l) &
237                   = zx_ctr(i, l) + zx_dtr(i, l)*local_tr(i, l-1)
238           ENDDO
239        ENDDO
240    
241        ! Calcul des tendances du traceur dans le sol et dans l'atmosphere
242    
243        DO l = 1, klev
244           DO i = 1, klon
245              d_tr(i, l) = local_tr(i, l) - tr(i, l)
246           ENDDO
247        ENDDO
248        d_trs = local_trs - trs
249    
250      END SUBROUTINE cltracrn
251    
252    end module cltracrn_m
Legend:
Removed from v.10  
changed lines
  Added in v.254
  ViewVC Help
Powered by ViewVC 1.1.21