libf/phylmd/cltracrn.f90

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,v 1.2 2005/05/25 13:10:09

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

  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

END SUBROUTINE cltracrn
1	SUBROUTINE cltracrn( itr, dtime,u1lay, v1lay, &
2	coef,t,ftsol,pctsrf, &
3	tr,trs,paprs,pplay,delp, &
4	masktr,fshtr,hsoltr,tautr,vdeptr, &
5	lat, &
6	d_tr,d_trs )
7
8	! From phylmd/cltracrn.F,v 1.2 2005/05/25 13:10:09
9
10	use indicesol, only: nbsrf
11	use dimphy, only: klon, klev
12	use SUPHEC_M, only: RD, rg
13
14	IMPLICIT none
15	!======================================================================
16	! Auteur(s): Alex/LMD) date: fev 99
17	! inspire de clqh + clvent
18	! Objet: diffusion verticale de traceurs avec quantite de traceur ds
19	! le sol ( reservoir de sol de radon )
20	!
21	! note : pour l'instant le traceur dans le sol et le flux sont
22	! calcules mais ils ne servent que de diagnostiques
23	! seule la tendance sur le traceur est sortie (d_tr)
24	!======================================================================
25	! Arguments:
26	! itr--- -input-R- le type de traceur 1- Rn 2 - Pb
27	! dtime----input-R- intervalle du temps (en second)
28	! u1lay----input-R- vent u de la premiere couche (m/s)
29	! v1lay----input-R- vent v de la premiere couche (m/s)
30	! coef-----input-R- le coefficient d'echange (m**2/s) l>1
31	! paprs----input-R- pression a inter-couche (Pa)
32	! pplay----input-R- pression au milieu de couche (Pa)
33	! delp-----input-R- epaisseur de couche (Pa)
34	! ftsol----input-R- temperature du sol (en Kelvin)
35	! tr-------input-R- traceurs
36	! trs------input-R- traceurs dans le sol
37	! masktr---input-R- Masque reservoir de sol traceur (1 = reservoir)
38	! fshtr----input-R- Flux surfacique de production dans le sol
39	! tautr----input-R- Constante de decroissance du traceur
40	! vdeptr---input-R- Vitesse de depot sec dans la couche brownienne
41	! hsoltr---input-R- Epaisseur equivalente du reservoir de sol
42	! lat-----input-R- latitude en degree
43	! d_tr-----output-R- le changement de "tr"
44	! d_trs----output-R- le changement de "trs"
45	!======================================================================
46	!======================================================================
47	REAL, intent(in):: dtime
48	REAL u1lay(klon), v1lay(klon)
49	REAL coef(klon,klev)
50	REAL, intent(in):: t(klon,klev) ! temperature (K)
51	real ftsol(klon,nbsrf), pctsrf(klon,nbsrf)
52	REAL tr(klon,klev), trs(klon)
53	REAL, intent(in):: paprs(klon,klev+1)
54	real, intent(in):: pplay(klon,klev)
55	real delp(klon,klev)
56	REAL masktr(klon)
57	REAL fshtr(klon)
58	REAL hsoltr
59	REAL tautr
60	REAL vdeptr
61	REAL, intent(in):: lat(klon)
62	REAL d_tr(klon,klev)
63	!======================================================================
64	REAL d_trs(klon) ! (diagnostic) traceur ds le sol
65	!======================================================================
66	INTEGER i, k, itr, n, l
67	REAL rotrhi(klon)
68	REAL zx_coef(klon,klev)
69	REAL zx_buf(klon)
70	REAL zx_ctr(klon,klev)
71	REAL zx_dtr(klon,klev)
72	REAL zx_trs(klon)
73	REAL zx_a, zx_b
74
75	REAL local_tr(klon,klev)
76	REAL local_trs(klon)
77	REAL zts(klon)
78	REAL zx_alpha1(klon), zx_alpha2(klon)
79	!======================================================================
80	!AA Pour l'instant les 4 types de surface ne sont pas pris en compte
81	!AA On fabrique avec zts un champ de temperature de sol
82	!AA que le pondere par la fraction de nature de sol.
83	!
84	! print*,'PASSAGE DANS CLTRACRN'
85
86	DO i = 1,klon
87	zts(i) = 0.
88	ENDDO
89	!
90	DO n=1,nbsrf
91	DO i = 1,klon
92	zts(i) = zts(i) + ftsol(i,n)*pctsrf(i,n)
93	ENDDO
94	ENDDO
95	!
96	DO i = 1,klon
97	rotrhi(i) = RD * zts(i) / hsoltr
98	END DO
99	!
100	DO k = 1, klev
101	DO i = 1, klon
102	local_tr(i,k) = tr(i,k)
103	ENDDO
104	ENDDO
105	!
106	DO i = 1, klon
107	local_trs(i) = trs(i)
108	ENDDO
109	!======================================================================
110	!AA Attention si dans clmain zx_alf1(i) = 1.0
111	!AA Il doit y avoir coherence (dc 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) = coef(i,1) &
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	!-------------------------
165	! Au dessus des continents
166	!-------------------------
167	! Le pb peut se deposer partout : vdeptr = 10-3 m/s
168	! Le Rn est traiter commme une couche Brownienne puisque vdeptr = 0.
169	!
170	DO i = 1, klon
171	!
172	IF ( NINT(masktr(i)) .EQ. 1 ) THEN
173	zx_trs(i) = local_trs(i)
174	zx_a = zx_trs(i) &
175	+fshtr(i)dtimerotrhi(i) &
176	+rotrhi(i)masktr(i)zx_coef(i,1)/RG &
177	(zx_ctr(i,1)(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2)) &
178	+zx_alpha2(i)*zx_ctr(i,2))
179	zx_b = 1. + rotrhi(i)masktr(i)zx_coef(i,1)/RG &
180	* (1.-zx_dtr(i,1) &
181	(zx_alpha1(i)+zx_alpha2(i)zx_dtr(i,2))) &
182	+ dtime / tautr &
183	!AA: Pour l'instant, pour aller vite, le depot sec est traite
184	! comme une decroissance
185	+ dtime * vdeptr / hsoltr
186	zx_trs(i) = zx_a / zx_b
187	local_trs(i) = zx_trs(i)
188	ENDIF
189	!
190	! Si on est entre 60N et 70N on divise par 2 l'emanation
191	!--------------------------------------------------------
192	!
193	IF &
194	( (itr.eq.1.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GE.60. &
195	.AND.lat(i).LE.70.) &
196	.OR. &
197	(itr.eq.2.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GE.60. &
198	.AND.lat(i).LE.70.) ) &
199	THEN
200	zx_trs(i) = local_trs(i)
201	zx_a = zx_trs(i) &
202	+(fshtr(i)/2.)dtimerotrhi(i) &
203	+rotrhi(i)masktr(i)zx_coef(i,1)/RG &
204	(zx_ctr(i,1)(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2)) &
205	+zx_alpha2(i)*zx_ctr(i,2))
206	zx_b = 1. + rotrhi(i)masktr(i)zx_coef(i,1)/RG &
207	* (1.-zx_dtr(i,1) &
208	(zx_alpha1(i)+zx_alpha2(i)zx_dtr(i,2))) &
209	+ dtime / tautr &
210	+ dtime * vdeptr / hsoltr
211	zx_trs(i) = zx_a / zx_b
212	local_trs(i) = zx_trs(i)
213	ENDIF
214	!
215	!----------------------------------------------
216	! Au dessus des oceans et aux hautes latitudes
217	!----------------------------------------------
218	!
219	! au dessous de -60S pas d'emission de radon au dessus
220	! des oceans et des continents
221	!---------------------------------------------------------------
222
223	IF ( (itr.EQ.1.AND.NINT(masktr(i)).EQ.0) &
224	.OR. &
225	(itr.EQ.1.AND.NINT(masktr(i)).EQ.1.AND.lat(i).LT.-60.)) &
226	THEN
227	zx_trs(i) = 0.
228	local_trs(i) = 0.
229	END IF
230
231	! au dessus de 70 N pas d'emission de radon au dessus
232	! des oceans et des continents
233	!--------------------------------------------------------------
234	IF ( (itr.EQ.1.AND.NINT(masktr(i)).EQ.0) &
235	.OR. &
236	(itr.EQ.1.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GT.70.)) &
237	THEN
238	zx_trs(i) = 0.
239	local_trs(i) = 0.
240	END IF
241
242	! Au dessus des oceans la source est nulle
243	!-----------------------------------------
244	!
245	IF (itr.eq.1.AND.NINT(masktr(i)).EQ.0) THEN
246	zx_trs(i) = 0.
247	local_trs(i) = 0.
248	END IF
249	!
250	ENDDO ! sur le i=1,klon
251	!
252	!======================================================================
253	!==== une fois on a zx_trs, on peut faire l'iteration ========
254	!
255	DO i = 1, klon
256	local_tr(i,1) = zx_ctr(i,1)+zx_dtr(i,1)*zx_trs(i)
257	ENDDO
258	DO l = 2, klev
259	DO i = 1, klon
260	local_tr(i,l) &
261	= zx_ctr(i,l) + zx_dtr(i,l)*local_tr(i,l-1)
262	ENDDO
263	ENDDO
264	!======================================================================
265	!== Calcul des tendances du traceur ds le sol et dans l'atmosphere
266	!
267	DO l = 1, klev
268	DO i = 1, klon
269	d_tr(i,l) = local_tr(i,l) - tr(i,l)
270	ENDDO
271	ENDDO
272	DO i = 1, klon
273	d_trs(i) = local_trs(i) - trs(i)
274	ENDDO
275
276	END SUBROUTINE cltracrn