| 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 | USE lmdz_yomcst |
|---|
| 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 : id_rn(radon), id_pb(plomb) |
|---|
| 27 | ! dtime....input-R- intervalle du temps (en secondes) ~ pdtphys |
|---|
| 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 | ! cdrag....input-R- cdrag |
|---|
| 31 | ! coef.....input-R- le coefficient d'echange (m**2/s) l>1, valable uniquement pour k entre 2 et klev |
|---|
| 32 | ! t........input-R- temperature (K) |
|---|
| 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 | !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)) == 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==id_rn.AND.NINT(masktr(i))==1.AND.lat(i)>=60..AND.lat(i)<=70.).OR. & |
|---|
| 200 | (itr==id_pb.AND.NINT(masktr(i))==1.AND.lat(i)>=60..AND.lat(i)<=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==id_rn.AND.NINT(masktr(i))==0).OR. & |
|---|
| 226 | (itr==id_rn.AND.NINT(masktr(i))==1.AND.lat(i)<-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==id_rn.AND.NINT(masktr(i))==0).OR. & |
|---|
| 235 | (itr==id_rn.AND.NINT(masktr(i))==1.AND.lat(i)>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==id_rn.AND.NINT(masktr(i))==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 |
|---|
