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 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 pplay(klon,klev), 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
     !AA         zx_alpha1(i) = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
     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
10	! From phylmd/cltracrn.F,v 1.2 2005/05/25 13:10:09
11
12	use indicesol, only: nbsrf
13	use dimphy, only: klon, klev
14	use YOMCST, only: RD, rg
15
16	IMPLICIT none
17	!======================================================================
18	! Auteur(s): Alex/LMD) date: fev 99
19	! inspire de clqh + clvent
20	! Objet: diffusion verticale de traceurs avec quantite de traceur ds
21	! le sol ( reservoir de sol de radon )
22	!
23	! note : pour l'instant le traceur dans le sol et le flux sont
24	! calcules mais ils ne servent que de diagnostiques
25	! seule la tendance sur le traceur est sortie (d_tr)
26	!======================================================================
27	! Arguments:
28	! itr--- -input-R- le type de traceur 1- Rn 2 - Pb
29	! dtime----input-R- intervalle du temps (en second)
30	! u1lay----input-R- vent u de la premiere couche (m/s)
31	! v1lay----input-R- vent v de la premiere couche (m/s)
32	! coef-----input-R- le coefficient d'echange (m**2/s) l>1
33	! paprs----input-R- pression a inter-couche (Pa)
34	! pplay----input-R- pression au milieu de couche (Pa)
35	! delp-----input-R- epaisseur de couche (Pa)
36	! ftsol----input-R- temperature du sol (en Kelvin)
37	! tr-------input-R- traceurs
38	! trs------input-R- traceurs dans le sol
39	! masktr---input-R- Masque reservoir de sol traceur (1 = reservoir)
40	! fshtr----input-R- Flux surfacique de production dans le sol
41	! tautr----input-R- Constante de decroissance du traceur
42	! vdeptr---input-R- Vitesse de depot sec dans la couche brownienne
43	! hsoltr---input-R- Epaisseur equivalente du reservoir de sol
44	! lat-----input-R- latitude en degree
45	! d_tr-----output-R- le changement de "tr"
46	! d_trs----output-R- le changement de "trs"
47	!======================================================================
48	!======================================================================
49	REAL 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 pplay(klon,klev), 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	!AA zx_alpha1(i) = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
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	!-------------------------
167	! Au dessus des continents
168	!-------------------------
169	! Le pb peut se deposer partout : vdeptr = 10-3 m/s
170	! Le Rn est traiter commme une couche Brownienne puisque vdeptr = 0.
171	!
172	DO i = 1, klon
173	!
174	IF ( NINT(masktr(i)) .EQ. 1 ) THEN
175	zx_trs(i) = local_trs(i)
176	zx_a = zx_trs(i) &
177	+fshtr(i)dtimerotrhi(i) &
178	+rotrhi(i)masktr(i)zx_coef(i,1)/RG &
179	(zx_ctr(i,1)(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2)) &
180	+zx_alpha2(i)*zx_ctr(i,2))
181	zx_b = 1. + rotrhi(i)masktr(i)zx_coef(i,1)/RG &
182	* (1.-zx_dtr(i,1) &
183	(zx_alpha1(i)+zx_alpha2(i)zx_dtr(i,2))) &
184	+ dtime / tautr &
185	!AA: Pour l'instant, pour aller vite, le depot sec est traite
186	! comme une decroissance
187	+ dtime * vdeptr / hsoltr
188	zx_trs(i) = zx_a / zx_b
189	local_trs(i) = zx_trs(i)
190	ENDIF
191	!
192	! Si on est entre 60N et 70N on divise par 2 l'emanation
193	!--------------------------------------------------------
194	!
195	IF &
196	( (itr.eq.1.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GE.60. &
197	.AND.lat(i).LE.70.) &
198	.OR. &
199	(itr.eq.2.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GE.60. &
200	.AND.lat(i).LE.70.) ) &
201	THEN
202	zx_trs(i) = local_trs(i)
203	zx_a = zx_trs(i) &
204	+(fshtr(i)/2.)dtimerotrhi(i) &
205	+rotrhi(i)masktr(i)zx_coef(i,1)/RG &
206	(zx_ctr(i,1)(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2)) &
207	+zx_alpha2(i)*zx_ctr(i,2))
208	zx_b = 1. + rotrhi(i)masktr(i)zx_coef(i,1)/RG &
209	* (1.-zx_dtr(i,1) &
210	(zx_alpha1(i)+zx_alpha2(i)zx_dtr(i,2))) &
211	+ dtime / tautr &
212	+ dtime * vdeptr / hsoltr
213	zx_trs(i) = zx_a / zx_b
214	local_trs(i) = zx_trs(i)
215	ENDIF
216	!
217	!----------------------------------------------
218	! Au dessus des oceans et aux hautes latitudes
219	!----------------------------------------------
220	!
221	! au dessous de -60S pas d'emission de radon au dessus
222	! des oceans et des continents
223	!---------------------------------------------------------------
224
225	IF ( (itr.EQ.1.AND.NINT(masktr(i)).EQ.0) &
226	.OR. &
227	(itr.EQ.1.AND.NINT(masktr(i)).EQ.1.AND.lat(i).LT.-60.)) &
228	THEN
229	zx_trs(i) = 0.
230	local_trs(i) = 0.
231	END IF
232
233	! au dessus de 70 N pas d'emission de radon au dessus
234	! des oceans et des continents
235	!--------------------------------------------------------------
236	IF ( (itr.EQ.1.AND.NINT(masktr(i)).EQ.0) &
237	.OR. &
238	(itr.EQ.1.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GT.70.)) &
239	THEN
240	zx_trs(i) = 0.
241	local_trs(i) = 0.
242	END IF
243
244	! Au dessus des oceans la source est nulle
245	!-----------------------------------------
246	!
247	IF (itr.eq.1.AND.NINT(masktr(i)).EQ.0) THEN
248	zx_trs(i) = 0.
249	local_trs(i) = 0.
250	END IF
251	!
252	ENDDO ! sur le i=1,klon
253	!
254	!======================================================================
255	!==== une fois on a zx_trs, on peut faire l'iteration ========
256	!
257	DO i = 1, klon
258	local_tr(i,1) = zx_ctr(i,1)+zx_dtr(i,1)*zx_trs(i)
259	ENDDO
260	DO l = 2, klev
261	DO i = 1, klon
262	local_tr(i,l) &
263	= zx_ctr(i,l) + zx_dtr(i,l)*local_tr(i,l-1)
264	ENDDO
265	ENDDO
266	!======================================================================
267	!== Calcul des tendances du traceur ds le sol et dans l'atmosphere
268	!
269	DO l = 1, klev
270	DO i = 1, klon
271	d_tr(i,l) = local_tr(i,l) - tr(i,l)
272	ENDDO
273	ENDDO
274	DO i = 1, klon
275	d_trs(i) = local_trs(i) - trs(i)
276	ENDDO
277
278	END SUBROUTINE cltracrn