| 1 | !$Id $ |
|---|
| 2 | |
|---|
| 3 | SUBROUTINE cltracrn( itr, dtime,u1lay, v1lay, & |
|---|
| 4 | cdrag,coef,t,ftsol,pctsrf, & |
|---|
| 5 | tr,trs,paprs,pplay,delp, & |
|---|
| 6 | masktr,fshtr,hsoltr,tautr,vdeptr, & |
|---|
| 7 | lat,d_tr,d_trs ) |
|---|
| 8 | |
|---|
| 9 | USE dimphy |
|---|
| 10 | USE traclmdz_mod, ONLY : id_rn, id_pb |
|---|
| 11 | USE indice_sol_mod |
|---|
| 12 | |
|---|
| 13 | IMPLICIT NONE |
|---|
| 14 | !====================================================================== |
|---|
| 15 | ! Auteur(s): Alex/LMD) date: fev 99 |
|---|
| 16 | ! inspire de clqh + clvent |
|---|
| 17 | ! Objet: diffusion verticale de traceurs avec quantite de traceur ds |
|---|
| 18 | ! le sol ( reservoir de sol de radon ) |
|---|
| 19 | ! |
|---|
| 20 | ! note : pour l'instant le traceur dans le sol et le flux sont |
|---|
| 21 | ! calcules mais ils ne servent que de diagnostiques |
|---|
| 22 | ! seule la tendance sur le traceur est sortie (d_tr) |
|---|
| 23 | !--------------------------------------------------------------------- |
|---|
| 24 | ! Arguments: |
|---|
| 25 | ! itr......input-R- le type de traceur : id_rn(radon), id_pb(plomb) |
|---|
| 26 | ! dtime....input-R- intervalle du temps (en secondes) ~ pdtphys |
|---|
| 27 | ! u1lay....input-R- vent u de la premiere couche (m/s) |
|---|
| 28 | ! v1lay....input-R- vent v de la premiere couche (m/s) |
|---|
| 29 | ! cdrag....input-R- cdrag |
|---|
| 30 | ! coef.....input-R- le coefficient d'echange (m**2/s) l>1, valable uniquement pour k entre 2 et klev |
|---|
| 31 | ! t........input-R- temperature (K) |
|---|
| 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 | include "YOMCST.h" |
|---|
| 48 | ! |
|---|
| 49 | !Entrees |
|---|
| 50 | INTEGER,INTENT(IN) :: itr |
|---|
| 51 | REAL,INTENT(IN) :: dtime |
|---|
| 52 | REAL,DIMENSION(klon),INTENT(IN) :: u1lay, v1lay |
|---|
| 53 | REAL,DIMENSION(klon),INTENT(IN) :: cdrag |
|---|
| 54 | REAL,DIMENSION(klon,klev),INTENT(IN) :: coef, t |
|---|
| 55 | REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: ftsol, pctsrf |
|---|
| 56 | REAL,DIMENSION(klon,klev),INTENT(IN) :: tr |
|---|
| 57 | REAL,DIMENSION(klon),INTENT(IN) :: trs |
|---|
| 58 | REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs |
|---|
| 59 | REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay, delp |
|---|
| 60 | REAL,DIMENSION(klon),INTENT(IN) :: masktr |
|---|
| 61 | REAL,DIMENSION(klon),INTENT(IN) :: fshtr |
|---|
| 62 | REAL,INTENT(IN) :: hsoltr |
|---|
| 63 | REAL,INTENT(IN) :: tautr |
|---|
| 64 | REAL,INTENT(IN) :: vdeptr |
|---|
| 65 | REAL,DIMENSION(klon),INTENT(IN) :: lat |
|---|
| 66 | |
|---|
| 67 | !Sorties |
|---|
| 68 | REAL,DIMENSION(klon,klev),INTENT(OUT) :: d_tr |
|---|
| 69 | REAL,DIMENSION(klon),INTENT(OUT) :: d_trs ! (diagnostic) traceur ds le sol |
|---|
| 70 | |
|---|
| 71 | !Locales |
|---|
| 72 | REAL,DIMENSION(klon,klev) :: flux_tr ! (diagnostic) flux de traceur |
|---|
| 73 | INTEGER :: i, k, n, l |
|---|
| 74 | REAL,DIMENSION(klon) :: rotrhi |
|---|
| 75 | REAL,DIMENSION(klon,klev) :: zx_coef |
|---|
| 76 | REAL,DIMENSION(klon) :: zx_buf |
|---|
| 77 | REAL,DIMENSION(klon,klev) :: zx_ctr |
|---|
| 78 | REAL,DIMENSION(klon,klev) :: zx_dtr |
|---|
| 79 | REAL,DIMENSION(klon) :: zx_trs |
|---|
| 80 | REAL :: zx_a, zx_b |
|---|
| 81 | |
|---|
| 82 | REAL,DIMENSION(klon,klev) :: local_tr |
|---|
| 83 | REAL,DIMENSION(klon) :: local_trs |
|---|
| 84 | REAL,DIMENSION(klon) :: zts ! champ de temperature du sol |
|---|
| 85 | REAL,DIMENSION(klon) :: zx_alpha1, zx_alpha2 |
|---|
| 86 | !====================================================================== |
|---|
| 87 | !AA Pour l'instant les 4 types de surface ne sont pas pris en compte |
|---|
| 88 | !AA On fabrique avec zts un champ de temperature de sol |
|---|
| 89 | !AA que le pondere par la fraction de nature de sol. |
|---|
| 90 | |
|---|
| 91 | DO i = 1,klon |
|---|
| 92 | zts(i) = 0. |
|---|
| 93 | ENDDO |
|---|
| 94 | |
|---|
| 95 | DO n=1,nbsrf |
|---|
| 96 | DO i = 1,klon |
|---|
| 97 | zts(i) = zts(i) + ftsol(i,n)*pctsrf(i,n) |
|---|
| 98 | ENDDO |
|---|
| 99 | ENDDO |
|---|
| 100 | |
|---|
| 101 | DO i = 1,klon |
|---|
| 102 | rotrhi(i) = RD * zts(i) / hsoltr |
|---|
| 103 | ENDDO |
|---|
| 104 | |
|---|
| 105 | DO k = 1, klev |
|---|
| 106 | DO i = 1, klon |
|---|
| 107 | local_tr(i,k) = tr(i,k) |
|---|
| 108 | ENDDO |
|---|
| 109 | ENDDO |
|---|
| 110 | |
|---|
| 111 | DO i = 1, klon |
|---|
| 112 | local_trs(i) = trs(i) |
|---|
| 113 | ENDDO |
|---|
| 114 | !====================================================================== |
|---|
| 115 | !AA Attention si dans clmain zx_alf1(i) = 1.0 |
|---|
| 116 | !AA Il doit y avoir coherence (dc la meme chose ici) |
|---|
| 117 | |
|---|
| 118 | DO i = 1, klon |
|---|
| 119 | !AA zx_alpha1(i) = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2)) |
|---|
| 120 | zx_alpha1(i) = 1.0 |
|---|
| 121 | zx_alpha2(i) = 1.0 - zx_alpha1(i) |
|---|
| 122 | ENDDO |
|---|
| 123 | !====================================================================== |
|---|
| 124 | DO i = 1, klon |
|---|
| 125 | zx_coef(i,1) = cdrag(i)*(1.0+SQRT(u1lay(i)**2+v1lay(i)**2)) & |
|---|
| 126 | *pplay(i,1)/(RD*t(i,1)) |
|---|
| 127 | zx_coef(i,1) = zx_coef(i,1) * dtime*RG |
|---|
| 128 | ENDDO |
|---|
| 129 | |
|---|
| 130 | DO k = 2, klev |
|---|
| 131 | DO i = 1, klon |
|---|
| 132 | zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k)) & |
|---|
| 133 | *(paprs(i,k)*2/(t(i,k)+t(i,k-1))/RD)**2 |
|---|
| 134 | zx_coef(i,k) = zx_coef(i,k) * dtime*RG |
|---|
| 135 | ENDDO |
|---|
| 136 | ENDDO |
|---|
| 137 | !====================================================================== |
|---|
| 138 | DO i = 1, klon |
|---|
| 139 | zx_buf(i) = delp(i,klev) + zx_coef(i,klev) |
|---|
| 140 | zx_ctr(i,klev) = local_tr(i,klev)*delp(i,klev)/zx_buf(i) |
|---|
| 141 | zx_dtr(i,klev) = zx_coef(i,klev) / zx_buf(i) |
|---|
| 142 | ENDDO |
|---|
| 143 | |
|---|
| 144 | DO l = klev-1, 2 , -1 |
|---|
| 145 | DO i = 1, klon |
|---|
| 146 | zx_buf(i) = delp(i,l)+zx_coef(i,l) & |
|---|
| 147 | +zx_coef(i,l+1)*(1.-zx_dtr(i,l+1)) |
|---|
| 148 | |
|---|
| 149 | zx_ctr(i,l) = ( local_tr(i,l)*delp(i,l) & |
|---|
| 150 | + zx_coef(i,l+1)*zx_ctr(i,l+1) )/zx_buf(i) |
|---|
| 151 | zx_dtr(i,l) = zx_coef(i,l) / zx_buf(i) |
|---|
| 152 | ENDDO |
|---|
| 153 | ENDDO |
|---|
| 154 | |
|---|
| 155 | DO i = 1, klon |
|---|
| 156 | zx_buf(i) = delp(i,1) + zx_coef(i,2)*(1.-zx_dtr(i,2)) & |
|---|
| 157 | + masktr(i) * zx_coef(i,1) & |
|---|
| 158 | *( zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2) ) |
|---|
| 159 | |
|---|
| 160 | zx_ctr(i,1) = ( local_tr(i,1)*delp(i,1) & |
|---|
| 161 | + zx_ctr(i,2) & |
|---|
| 162 | *(zx_coef(i,2) & |
|---|
| 163 | - masktr(i) * zx_coef(i,1) & |
|---|
| 164 | *zx_alpha2(i) ) ) / zx_buf(i) |
|---|
| 165 | zx_dtr(i,1) = masktr(i) * zx_coef(i,1) / zx_buf(i) |
|---|
| 166 | ENDDO |
|---|
| 167 | !====================================================================== |
|---|
| 168 | ! Calculer d'abord local_trs nouvelle quantite dans le reservoir |
|---|
| 169 | ! de sol |
|---|
| 170 | !===================================================================== |
|---|
| 171 | |
|---|
| 172 | DO i = 1, klon |
|---|
| 173 | !------------------------- |
|---|
| 174 | ! Au dessus des continents |
|---|
| 175 | !-- |
|---|
| 176 | ! Le pb peut se deposer partout : vdeptr = 10-3 m/s |
|---|
| 177 | ! Le Rn est traiter commme une couche Brownienne puisque vdeptr = 0. |
|---|
| 178 | !------------------------------------------------------------------- |
|---|
| 179 | IF ( NINT(masktr(i)) .EQ. 1 ) THEN |
|---|
| 180 | zx_trs(i) = local_trs(i) |
|---|
| 181 | zx_a = zx_trs(i) & |
|---|
| 182 | +fshtr(i)*dtime*rotrhi(i) & |
|---|
| 183 | +rotrhi(i)*masktr(i)*zx_coef(i,1)/RG & |
|---|
| 184 | *(zx_ctr(i,1)*(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2)) & |
|---|
| 185 | +zx_alpha2(i)*zx_ctr(i,2)) |
|---|
| 186 | ! Pour l'instant, pour aller vite, le depot sec est traite comme une decroissance |
|---|
| 187 | zx_b = 1. + rotrhi(i)*masktr(i)*zx_coef(i,1)/RG & |
|---|
| 188 | * (1.-zx_dtr(i,1) & |
|---|
| 189 | *(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2))) & |
|---|
| 190 | + dtime / tautr & |
|---|
| 191 | + dtime * vdeptr / hsoltr |
|---|
| 192 | zx_trs(i) = zx_a / zx_b |
|---|
| 193 | local_trs(i) = zx_trs(i) |
|---|
| 194 | ENDIF |
|---|
| 195 | !-------------------------------------------------------- |
|---|
| 196 | ! Si on est entre 60N et 70N on divise par 2 l'emanation |
|---|
| 197 | !-------------------------------------------------------- |
|---|
| 198 | |
|---|
| 199 | IF ( (itr.eq.id_rn.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GE.60..AND.lat(i).LE.70.).OR. & |
|---|
| 200 | (itr.eq.id_pb.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GE.60..AND.lat(i).LE.70.) ) THEN |
|---|
| 201 | zx_trs(i) = local_trs(i) |
|---|
| 202 | zx_a = zx_trs(i) & |
|---|
| 203 | +(fshtr(i)/2.)*dtime*rotrhi(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 | ! |
|---|
| 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 | ! |
|---|
| 214 | zx_trs(i) = zx_a / zx_b |
|---|
| 215 | local_trs(i) = zx_trs(i) |
|---|
| 216 | ENDIF |
|---|
| 217 | |
|---|
| 218 | !---------------------------------------------- |
|---|
| 219 | ! Au dessus des oceans et aux hautes latitudes |
|---|
| 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.id_rn.AND.NINT(masktr(i)).EQ.0).OR. & |
|---|
| 226 | (itr.EQ.id_rn.AND.NINT(masktr(i)).EQ.1.AND.lat(i).LT.-60.)) 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.id_rn.AND.NINT(masktr(i)).EQ.0).OR. & |
|---|
| 235 | (itr.EQ.id_rn.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GT.70.)) THEN |
|---|
| 236 | zx_trs(i) = 0. |
|---|
| 237 | local_trs(i) = 0. |
|---|
| 238 | END IF |
|---|
| 239 | !--------------------------------------------- |
|---|
| 240 | ! Au dessus des oceans la source est nulle |
|---|
| 241 | !-------------------------------------------- |
|---|
| 242 | |
|---|
| 243 | IF (itr.eq.id_rn.AND.NINT(masktr(i)).EQ.0) THEN |
|---|
| 244 | zx_trs(i) = 0. |
|---|
| 245 | local_trs(i) = 0. |
|---|
| 246 | END IF |
|---|
| 247 | |
|---|
| 248 | ENDDO ! sur le i=1,klon |
|---|
| 249 | ! |
|---|
| 250 | !====================================================================== |
|---|
| 251 | ! Une fois on a zx_trs, on peut faire l'iteration |
|---|
| 252 | !====================================================================== |
|---|
| 253 | |
|---|
| 254 | DO i = 1, klon |
|---|
| 255 | local_tr(i,1) = zx_ctr(i,1)+zx_dtr(i,1)*zx_trs(i) |
|---|
| 256 | ENDDO |
|---|
| 257 | DO l = 2, klev |
|---|
| 258 | DO i = 1, klon |
|---|
| 259 | local_tr(i,l) = zx_ctr(i,l) + zx_dtr(i,l)*local_tr(i,l-1) |
|---|
| 260 | ENDDO |
|---|
| 261 | ENDDO |
|---|
| 262 | !====================================================================== |
|---|
| 263 | ! Calcul du flux de traceur (flux_tr): UA/(m**2 s) |
|---|
| 264 | !====================================================================== |
|---|
| 265 | DO i = 1, klon |
|---|
| 266 | flux_tr(i,1) = masktr(i)*zx_coef(i,1)/RG & |
|---|
| 267 | * (zx_alpha1(i)*local_tr(i,1)+zx_alpha2(i)*local_tr(i,2) & |
|---|
| 268 | -zx_trs(i)) / dtime |
|---|
| 269 | ENDDO |
|---|
| 270 | DO l = 2, klev |
|---|
| 271 | DO i = 1, klon |
|---|
| 272 | flux_tr(i,l) = zx_coef(i,l)/RG & |
|---|
| 273 | * (local_tr(i,l)-local_tr(i,l-1)) / dtime |
|---|
| 274 | ENDDO |
|---|
| 275 | ENDDO |
|---|
| 276 | !====================================================================== |
|---|
| 277 | ! Calcul des tendances du traceur ds le sol et dans l'atmosphere |
|---|
| 278 | !====================================================================== |
|---|
| 279 | DO l = 1, klev |
|---|
| 280 | DO i = 1, klon |
|---|
| 281 | d_tr(i,l) = local_tr(i,l) - tr(i,l) |
|---|
| 282 | ENDDO |
|---|
| 283 | ENDDO |
|---|
| 284 | DO i = 1, klon |
|---|
| 285 | d_trs(i) = local_trs(i) - trs(i) |
|---|
| 286 | ENDDO |
|---|
| 287 | |
|---|
| 288 | END SUBROUTINE cltracrn |
|---|
