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 )

  ! 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

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