trunk/phylmd/cltracrn.f

module cltracrn_m

  IMPLICIT none

contains 

  SUBROUTINE cltracrn(itr, dtime, u1lay, v1lay, coef, cdragh, 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: Alexandre ARMENGAUD
    ! 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

    INTEGER itr
    ! itr--- -input-R- le type de traceur 1- Rn 2 - Pb

    REAL, intent(in):: dtime
    ! dtime----input-R- intervalle de temps (en second)
    REAL, intent(in):: u1lay(klon), v1lay(klon) ! vent de la premiere
                                                ! couche (m/s)
    REAL coef(:, 2:) ! (klon, 2:klev)
    ! coef-----input-R- le coefficient d'echange (m**2/s) l>1
    real cdragh(:) ! klon
    REAL, intent(in):: t(klon, klev) ! temperature (K)
    real, intent(in):: ftsol(klon, nbsrf), pctsrf(klon, nbsrf) 
    ! ftsol----input-R- temperature du sol (en Kelvin)
    REAL, intent(in):: tr(klon, klev) ! traceur
    REAL, intent(in):: trs(:) ! (klon) traceur dans le sol
    REAL, intent(in):: paprs(klon, klev+1)
    ! paprs----input-R- pression a l'inter-couche (Pa)
    real, intent(in):: pplay(klon, klev)
    ! pplay----input-R- pression au milieu de couche (Pa)
    real delp(klon, klev)
    ! delp-----input-R- epaisseur de couche (Pa)
    REAL masktr(klon) 
    ! masktr---input-R- Masque reservoir de sol traceur (1 = reservoir)
    REAL fshtr(klon) 
    ! fshtr----input-R- Flux surfacique de production dans le sol
    REAL hsoltr
    ! hsoltr---input-R- Epaisseur equivalente du reservoir de sol
    REAL tautr
    ! tautr----input-R- Constante de decroissance du traceur

    REAL, intent(in):: vdeptr
    ! vitesse de d\'ep\^ot sec dans la couche brownienne

    REAL, intent(in):: lat(klon) 
    ! lat-----input-R- latitude en degree
    REAL d_tr(klon, klev)
    ! d_tr-----output-R- le changement de "tr"

    REAL, intent(out):: d_trs(:) ! (klon) (diagnostic) changement de "trs"

    ! Local:
    INTEGER i, k, 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_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

    local_trs = trs

    ! Attention si dans pbl_surface 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) = cdragh(i) &
            * (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_a = local_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 d\'ep\^ot sec est trait\'e
          ! comme une d\'ecroissance :
          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
          local_trs(i) = zx_a / zx_b
       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_a = local_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
          local_trs(i) = zx_a / zx_b
       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
          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
          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
          local_trs(i) = 0.
       END IF
    ENDDO ! sur le i=1, klon

    ! une fois qu'on a local_trs, on peut faire l'iteration

    DO i = 1, klon
       local_tr(i, 1) = zx_ctr(i, 1)+zx_dtr(i, 1)*local_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
    d_trs = local_trs - trs

  END SUBROUTINE cltracrn

end module cltracrn_m
1	module cltracrn_m
2
3	IMPLICIT none
4
5	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 pbl_surface 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)dtimerotrhi(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.)dtimerotrhi(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