trunk/phylmd/cltracrn.f

module cltracrn_m

  IMPLICIT none

contains 

  SUBROUTINE cltracrn(itr, dtime, u1lay, v1lay, coef, t, ftsol, pctsrf, tr, &
       trs, paprs, pplay, delp, masktr, fshtr, hsoltr, tautr, vdeptr, lat, &
       d_tr, d_trs)

    ! From phylmd/cltracrn.F, version 1.2 2005/05/25 13:10:09

    ! Author: Alex A.
    ! Date: february 1999
    ! Inspiré de clqh et clvent

    ! Objet: diffusion verticale de traceurs avec quantité de traceur
    ! dans le sol (réservoir de sol de radon)

    ! Note : pour l'instant le traceur dans le sol et le flux sont
    ! calculés mais ils ne servent que de diagnostics. Seule la
    ! tendance sur le traceur est sortie (d_tr).

    use indicesol, only: nbsrf
    use dimphy, only: klon, klev
    use SUPHEC_M, only: RD, rg

    ! Arguments:
    ! itr--- -input-R- le type de traceur 1- Rn 2 - Pb
    ! dtime----input-R- intervalle de 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 l'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)

    !------------------------------------------------------------------

    ! Pour l'instant les quatre types de surface ne sont pas pris en
    ! compte. On fabrique avec zts un champ de température de sol que
    ! l'on pondère par la fraction de sol.

    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

    ! Attention si dans clmain zx_alf1(i) = 1.
    ! Il doit y avoir coherence (donc 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 traité 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))
          ! Pour l'instant, pour aller vite, le depot sec est traite
          ! comme une decroissance :
          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

       ! 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 qu'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 dans 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

  END SUBROUTINE cltracrn

end module cltracrn_m
1	guez	78	module cltracrn_m
2	guez	3
3	guez	78	IMPLICIT none
4	guez	3
5	guez	78	contains
6	guez	3
7	guez	78	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	guez	3
11	guez	78	! From phylmd/cltracrn.F, version 1.2 2005/05/25 13:10:09
12	guez	3
13	guez	78	! Author: Alex A.
14			! Date: february 1999
15			! Inspiré de clqh et clvent
16	guez	3
17	guez	78	! Objet: diffusion verticale de traceurs avec quantité de traceur
18			! dans le sol (réservoir de sol de radon)
19	guez	3
20	guez	78	! 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	guez	3
24	guez	78	use indicesol, only: nbsrf
25			use dimphy, only: klon, klev
26			use SUPHEC_M, only: RD, rg
27	guez	3
28	guez	78	! Arguments:
29			! itr--- -input-R- le type de traceur 1- Rn 2 - Pb
30			! dtime----input-R- intervalle de temps (en second)
31			! u1lay----input-R- vent u de la premiere couche (m/s)
32			! v1lay----input-R- vent v de la premiere couche (m/s)
33			! coef-----input-R- le coefficient d'echange (m**2/s) l>1
34			! paprs----input-R- pression a l'inter-couche (Pa)
35			! pplay----input-R- pression au milieu de couche (Pa)
36			! delp-----input-R- epaisseur de couche (Pa)
37			! ftsol----input-R- temperature du sol (en Kelvin)
38			! tr-------input-R- traceurs
39			! trs------input-R- traceurs dans le sol
40			! masktr---input-R- Masque reservoir de sol traceur (1 = reservoir)
41			! fshtr----input-R- Flux surfacique de production dans le sol
42			! tautr----input-R- Constante de decroissance du traceur
43			! vdeptr---input-R- Vitesse de depot sec dans la couche brownienne
44			! hsoltr---input-R- Epaisseur equivalente du reservoir de sol
45			! lat-----input-R- latitude en degree
46			! d_tr-----output-R- le changement de "tr"
47			! d_trs----output-R- le changement de "trs"
48	guez	3
49	guez	78	REAL, intent(in):: dtime
50			REAL u1lay(klon), v1lay(klon)
51			REAL coef(klon, klev)
52			REAL, intent(in):: t(klon, klev) ! temperature (K)
53			real ftsol(klon, nbsrf), pctsrf(klon, nbsrf)
54			REAL tr(klon, klev), trs(klon)
55			REAL, intent(in):: paprs(klon, klev+1)
56			real, intent(in):: pplay(klon, klev)
57			real delp(klon, klev)
58			REAL masktr(klon)
59			REAL fshtr(klon)
60			REAL hsoltr
61			REAL tautr
62			REAL vdeptr
63			REAL, intent(in):: lat(klon)
64			REAL d_tr(klon, klev)
65
66			REAL d_trs(klon) ! (diagnostic) traceur ds le sol
67
68			INTEGER i, k, itr, 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_trs(klon)
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			DO i = 1, klon
109			local_trs(i) = trs(i)
110			ENDDO
111
112			! Attention si dans clmain zx_alf1(i) = 1.
113			! Il doit y avoir coherence (donc la meme chose ici)
114
115			DO i = 1, klon
116			zx_alpha1(i) = 1.0
117			zx_alpha2(i) = 1.0 - zx_alpha1(i)
118			ENDDO
119
120			DO i = 1, klon
121			zx_coef(i, 1) = coef(i, 1) &
122			* (1.0+SQRT(u1lay(i)2+v1lay(i)2)) &
123			* pplay(i, 1)/(RD*t(i, 1))
124			zx_coef(i, 1) = zx_coef(i, 1) * dtime*RG
125			ENDDO
126
127			DO k = 2, klev
128			DO i = 1, klon
129			zx_coef(i, k) = coef(i, k)*RG/(pplay(i, k-1)-pplay(i, k)) &
130			(paprs(i, k)2/(t(i, k)+t(i, k-1))/RD)**2
131			zx_coef(i, k) = zx_coef(i, k) * dtime*RG
132			ENDDO
133			ENDDO
134
135			DO i = 1, klon
136			zx_buf(i) = delp(i, klev) + zx_coef(i, klev)
137			zx_ctr(i, klev) = local_tr(i, klev)*delp(i, klev)/zx_buf(i)
138			zx_dtr(i, klev) = zx_coef(i, klev) / zx_buf(i)
139			ENDDO
140
141			DO l = klev-1, 2 , -1
142			DO i = 1, klon
143			zx_buf(i) = delp(i, l)+zx_coef(i, l) &
144			+zx_coef(i, l+1)*(1.-zx_dtr(i, l+1))
145			zx_ctr(i, l) = (local_tr(i, l)*delp(i, l) &
146			+ zx_coef(i, l+1)*zx_ctr(i, l+1))/zx_buf(i)
147			zx_dtr(i, l) = zx_coef(i, l) / zx_buf(i)
148			ENDDO
149			ENDDO
150
151			DO i = 1, klon
152			zx_buf(i) = delp(i, 1) + zx_coef(i, 2)*(1.-zx_dtr(i, 2)) &
153			+ masktr(i) * zx_coef(i, 1) &
154			(zx_alpha1(i)+zx_alpha2(i)zx_dtr(i, 2))
155			zx_ctr(i, 1) = (local_tr(i, 1)*delp(i, 1) &
156			+ zx_ctr(i, 2) &
157			*(zx_coef(i, 2) &
158			- masktr(i) * zx_coef(i, 1) &
159			*zx_alpha2(i))) / zx_buf(i)
160			zx_dtr(i, 1) = masktr(i) * zx_coef(i, 1) / zx_buf(i)
161			ENDDO
162
163			! Calculer d'abord local_trs nouvelle quantite dans le reservoir
164			! de sol
165
166			! Au dessus des continents
167			! Le pb peut se deposer partout : vdeptr = 10-3 m/s
168			! Le Rn est traité commme une couche Brownienne puisque vdeptr = 0.
169
170			DO i = 1, klon
171			IF (NINT(masktr(i)) .EQ. 1) THEN
172			zx_trs(i) = local_trs(i)
173			zx_a = zx_trs(i) &
174			+fshtr(i)dtimerotrhi(i) &
175			+rotrhi(i)masktr(i)zx_coef(i, 1)/RG &
176			(zx_ctr(i, 1)(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i, 2)) &
177			+zx_alpha2(i)*zx_ctr(i, 2))
178			! Pour l'instant, pour aller vite, le depot sec est traite
179			! comme une decroissance :
180			zx_b = 1. + rotrhi(i)masktr(i)zx_coef(i, 1)/RG &
181			* (1.-zx_dtr(i, 1) &
182			(zx_alpha1(i)+zx_alpha2(i)zx_dtr(i, 2))) &
183			+ dtime / tautr &
184			+ dtime * vdeptr / hsoltr
185			zx_trs(i) = zx_a / zx_b
186			local_trs(i) = zx_trs(i)
187			ENDIF
188
189			! Si on est entre 60N et 70N on divise par 2 l'emanation
190
191			IF ((itr.eq.1.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GE.60. &
192			.AND.lat(i).LE.70.) .OR. &
193			(itr.eq.2.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GE.60. &
194			.AND.lat(i).LE.70.)) THEN
195			zx_trs(i) = local_trs(i)
196			zx_a = zx_trs(i) &
197			+(fshtr(i)/2.)dtimerotrhi(i) &
198			+rotrhi(i)masktr(i)zx_coef(i, 1)/RG &
199			(zx_ctr(i, 1)(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i, 2)) &
200			+zx_alpha2(i)*zx_ctr(i, 2))
201			zx_b = 1. + rotrhi(i)masktr(i)zx_coef(i, 1)/RG &
202			* (1.-zx_dtr(i, 1) &
203			(zx_alpha1(i)+zx_alpha2(i)zx_dtr(i, 2))) &
204			+ dtime / tautr &
205			+ dtime * vdeptr / hsoltr
206			zx_trs(i) = zx_a / zx_b
207			local_trs(i) = zx_trs(i)
208			ENDIF
209
210			! Au dessus des oceans et aux hautes latitudes
211
212			! au dessous de -60S pas d'emission de radon au dessus
213			! des oceans et des continents
214
215			IF ((itr.EQ.1.AND.NINT(masktr(i)).EQ.0) .OR. &
216			(itr.EQ.1.AND.NINT(masktr(i)).EQ.1.AND.lat(i).LT.-60.)) THEN
217			zx_trs(i) = 0.
218			local_trs(i) = 0.
219			END IF
220
221			! au dessus de 70 N pas d'emission de radon au dessus
222			! des oceans et des continents
223
224			IF ((itr.EQ.1.AND.NINT(masktr(i)).EQ.0) .OR. &
225			(itr.EQ.1.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GT.70.)) THEN
226			zx_trs(i) = 0.
227			local_trs(i) = 0.
228			END IF
229
230			! Au dessus des oceans la source est nulle
231
232			IF (itr.eq.1.AND.NINT(masktr(i)).EQ.0) THEN
233			zx_trs(i) = 0.
234			local_trs(i) = 0.
235			END IF
236
237			ENDDO ! sur le i=1, klon
238
239			! une fois qu'on a zx_trs, on peut faire l'iteration
240
241			DO i = 1, klon
242			local_tr(i, 1) = zx_ctr(i, 1)+zx_dtr(i, 1)*zx_trs(i)
243			ENDDO
244			DO l = 2, klev
245			DO i = 1, klon
246			local_tr(i, l) &
247			= zx_ctr(i, l) + zx_dtr(i, l)*local_tr(i, l-1)
248			ENDDO
249			ENDDO
250
251			! Calcul des tendances du traceur dans le sol et dans l'atmosphere
252
253			DO l = 1, klev
254			DO i = 1, klon
255			d_tr(i, l) = local_tr(i, l) - tr(i, l)
256			ENDDO
257			ENDDO
258			DO i = 1, klon
259			d_trs(i) = local_trs(i) - trs(i)
260			ENDDO
261
262			END SUBROUTINE cltracrn
263
264			end module cltracrn_m