| 1 | ! $Id: stratosphere_mask.F90 5144 2024-07-29 21:01:04Z abarral $ |
|---|
| 2 | |
|---|
| 3 | SUBROUTINE stratosphere_mask(missing_val, pphis, t_seri, pplay, xlat) |
|---|
| 4 | |
|---|
| 5 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 6 | |
|---|
| 7 | ! determination of tropopause height and temperature from gridded temperature data |
|---|
| 8 | |
|---|
| 9 | ! reference: Reichler, T., M. Dameris, and R. Sausen (GRL, 10.1029/2003GL018240, 2003) |
|---|
| 10 | ! modified: 6/28/06 tjr |
|---|
| 11 | ! adapted to LMDZ by C. Kleinschmitt (2016-02-15) |
|---|
| 12 | ! committed to LMDz by O. Boucher (2016) with a mistake |
|---|
| 13 | ! mistake corrected by O. Boucher (2017-12-11) |
|---|
| 14 | |
|---|
| 15 | ! input: temp(nlon,nlat,nlev) 3D-temperature field |
|---|
| 16 | ! ps(nlon,nlat) 2D-surface pressure field |
|---|
| 17 | ! zs(nlon,nlat) 2D-surface height |
|---|
| 18 | ! nlon grid points in x |
|---|
| 19 | ! nlat grid points in y |
|---|
| 20 | ! pfull(nlon,nlat,nlev) full pressure levels in Pa |
|---|
| 21 | ! plimu upper limit for tropopause pressure |
|---|
| 22 | ! pliml lower limit for tropopause pressure |
|---|
| 23 | ! gamma tropopause criterion, e.g. -0.002 K/m |
|---|
| 24 | |
|---|
| 25 | ! output: p_tropopause(klon) tropopause pressure in Pa with missing values |
|---|
| 26 | ! t_tropopause(klon) tropopause temperature in K with missing values |
|---|
| 27 | ! z_tropopause(klon) tropopause height in m with missing values |
|---|
| 28 | ! stratomask stratospheric mask withtout missing values |
|---|
| 29 | ! ifil # of undetermined values |
|---|
| 30 | |
|---|
| 31 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 32 | |
|---|
| 33 | USE dimphy |
|---|
| 34 | USE phys_local_var_mod, ONLY: stratomask |
|---|
| 35 | USE phys_local_var_mod, ONLY: p_tropopause, z_tropopause, t_tropopause |
|---|
| 36 | USE lmdz_print_control, ONLY: lunout, prt_level |
|---|
| 37 | USE lmdz_yomcst |
|---|
| 38 | |
|---|
| 39 | IMPLICIT NONE |
|---|
| 40 | |
|---|
| 41 | REAL, INTENT(IN) :: missing_val ! missing value, also XIOS |
|---|
| 42 | REAL, DIMENSION(klon), INTENT(IN) :: pphis ! Geopotentiel de surface |
|---|
| 43 | REAL, DIMENSION(klon, klev), INTENT(IN) :: t_seri ! Temperature |
|---|
| 44 | REAL, DIMENSION(klon, klev), INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa) |
|---|
| 45 | REAL, DIMENSION(klon), INTENT(IN) :: xlat ! latitudes pour chaque point |
|---|
| 46 | |
|---|
| 47 | REAL, PARAMETER :: plimu = 45000. |
|---|
| 48 | REAL, PARAMETER :: pliml = 7500. |
|---|
| 49 | REAL, PARAMETER :: gamma = -0.002 |
|---|
| 50 | LOGICAL, PARAMETER :: dofill = .TRUE. |
|---|
| 51 | REAL, DIMENSION(klon) :: tp |
|---|
| 52 | REAL, DIMENSION(klev) :: t, p |
|---|
| 53 | INTEGER :: i, k, ifil |
|---|
| 54 | REAL :: ptrp, ttrp, ztrp, psrf, zsrf, pi |
|---|
| 55 | |
|---|
| 56 | pi = 4. * ATAN(1.) |
|---|
| 57 | |
|---|
| 58 | !--computing tropopause |
|---|
| 59 | DO i = 1, klon |
|---|
| 60 | DO k = 1, klev |
|---|
| 61 | t(k) = t_seri(i, klev + 1 - k) |
|---|
| 62 | p(k) = pplay(i, klev + 1 - k) |
|---|
| 63 | ENDDO |
|---|
| 64 | psrf = pplay(i, 1) |
|---|
| 65 | zsrf = pphis(i) / RG !--altitude de la surface |
|---|
| 66 | CALL twmo(missing_val, klev, t, p, psrf, zsrf, plimu, pliml, gamma, ptrp, ttrp, ztrp) |
|---|
| 67 | tp(i) = ptrp |
|---|
| 68 | p_tropopause(i) = ptrp |
|---|
| 69 | z_tropopause(i) = ztrp |
|---|
| 70 | t_tropopause(i) = ttrp |
|---|
| 71 | ENDDO |
|---|
| 72 | |
|---|
| 73 | !--filling holes in tp but not in p_tropopause |
|---|
| 74 | IF (dofill) THEN |
|---|
| 75 | ifil = 0 |
|---|
| 76 | DO i = 1, klon |
|---|
| 77 | IF (ABS(tp(i) / missing_val - 1.0)<0.01) THEN |
|---|
| 78 | !set missing values to very simple profile (neighbour averaging too expensive in LMDZ) |
|---|
| 79 | tp(i) = 50000. - 20000. * cos(xlat(i) / 360. * 2. * pi) |
|---|
| 80 | ifil = ifil + 1 |
|---|
| 81 | ENDIF |
|---|
| 82 | ENDDO |
|---|
| 83 | |
|---|
| 84 | ENDIF |
|---|
| 85 | |
|---|
| 86 | DO i = 1, klon |
|---|
| 87 | DO k = 1, klev |
|---|
| 88 | IF (pplay(i, k)<tp(i)) THEN |
|---|
| 89 | stratomask(i, k) = 1.0 |
|---|
| 90 | ELSE |
|---|
| 91 | stratomask(i, k) = 0.0 |
|---|
| 92 | ENDIF |
|---|
| 93 | ENDDO |
|---|
| 94 | ENDDO |
|---|
| 95 | |
|---|
| 96 | IF (ifil>0 .AND. prt_level >5) THEN |
|---|
| 97 | WRITE(lunout, *)'Tropopause: number of undetermined values =', ifil |
|---|
| 98 | ENDIF |
|---|
| 99 | |
|---|
| 100 | END SUBROUTINE stratosphere_mask |
|---|
| 101 | |
|---|
| 102 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 103 | ! twmo |
|---|
| 104 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 105 | |
|---|
| 106 | SUBROUTINE twmo(missing_val, level, t, p, ps, zs, plimu, pliml, gamma, ptrp, ttrp, ztrp) |
|---|
| 107 | |
|---|
| 108 | ! reference: Reichler, T., M. Dameris, and R. Sausen (GRL, 10.1029/2003GL018240, 2003) |
|---|
| 109 | |
|---|
| 110 | USE lmdz_yomcst |
|---|
| 111 | |
|---|
| 112 | IMPLICIT NONE |
|---|
| 113 | |
|---|
| 114 | INTEGER, INTENT(IN) :: level |
|---|
| 115 | REAL, INTENT(IN) :: missing_val |
|---|
| 116 | REAL, INTENT(IN), DIMENSION(level) :: t, p |
|---|
| 117 | REAL, INTENT(IN) :: plimu, pliml, gamma, ps, zs |
|---|
| 118 | REAL, INTENT(OUT) :: ptrp, ttrp, ztrp |
|---|
| 119 | |
|---|
| 120 | REAL, parameter :: deltaz = 2000.0 |
|---|
| 121 | |
|---|
| 122 | real :: faktor |
|---|
| 123 | real :: pmk, pm, a, b, tm, dtdp, dtdz, dlnp, tdlnp |
|---|
| 124 | real :: ag, bg, ptph, ttph, a0, b0 |
|---|
| 125 | real :: pm0, tm0, pmk0, dtdz0 |
|---|
| 126 | real :: p2km, asum, aquer |
|---|
| 127 | real :: pmk2, pm2, a2, b2, tm2, dtdp2, dtdz2 |
|---|
| 128 | INTEGER :: icount, jj, j |
|---|
| 129 | |
|---|
| 130 | ptrp = missing_val |
|---|
| 131 | ttrp = missing_val |
|---|
| 132 | ztrp = missing_val |
|---|
| 133 | |
|---|
| 134 | faktor = -RG / RD |
|---|
| 135 | |
|---|
| 136 | DO j = level, 2, -1 |
|---|
| 137 | |
|---|
| 138 | ! dt/dz |
|---|
| 139 | pmk = 0.5 * (p(j - 1)**rkappa + p(j)**rkappa) ! p**k at half level |
|---|
| 140 | pm = pmk**(1. / rkappa) ! p at half level |
|---|
| 141 | a = (t(j - 1) - t(j)) / (p(j - 1)**rkappa - p(j)**rkappa) ! dT/dp^k |
|---|
| 142 | b = t(j) - (a * p(j)**rkappa) |
|---|
| 143 | tm = a * pmk + b ! T at half level |
|---|
| 144 | dtdp = a * rkappa * pm**(rkappa - 1.) ! dT/dp at half level |
|---|
| 145 | dtdz = faktor * dtdp * pm / tm ! dT/dz at half level |
|---|
| 146 | |
|---|
| 147 | ! dt/dz valid? |
|---|
| 148 | IF (j==level) go to 999 ! no, start level, initialize first |
|---|
| 149 | IF (dtdz<=gamma) go to 999 ! no, dt/dz < -2 K/km |
|---|
| 150 | IF (dtdz0>gamma) go to 999 ! no, dt/dz below > -2 K/km |
|---|
| 151 | IF (pm>plimu) go to 999 ! no, pm too low |
|---|
| 152 | |
|---|
| 153 | ! dtdz is valid, calculate tropopause pressure ptph |
|---|
| 154 | ag = (dtdz - dtdz0) / (pmk - pmk0) |
|---|
| 155 | bg = dtdz0 - (ag * pmk0) |
|---|
| 156 | ptph = exp(log((gamma - bg) / ag) / rkappa) |
|---|
| 157 | |
|---|
| 158 | ! calculate temperature at this ptph assuming linear gamma |
|---|
| 159 | ! interpolation |
|---|
| 160 | ttph = tm0 |
|---|
| 161 | ttph = ttph - (bg * log(pm0) + ag * (pm0**rkappa) / rkappa) / faktor * t(j) |
|---|
| 162 | ttph = ttph + (bg * log(ptph) + ag * (ptph**rkappa) / rkappa) / faktor * t(j) |
|---|
| 163 | |
|---|
| 164 | IF (ptph<pliml) go to 999 |
|---|
| 165 | IF (ptph>plimu) go to 999 |
|---|
| 166 | |
|---|
| 167 | ! 2nd test: dtdz above 2 km must not exceed gamma |
|---|
| 168 | p2km = ptph + deltaz * (pm / tm) * faktor ! p at ptph + 2km |
|---|
| 169 | asum = 0.0 ! dtdz above |
|---|
| 170 | icount = 0 ! number of levels above |
|---|
| 171 | |
|---|
| 172 | ! test until apm < p2km |
|---|
| 173 | do jj = j, 2, -1 |
|---|
| 174 | |
|---|
| 175 | pmk2 = .5 * (p(jj - 1)**rkappa + p(jj)**rkappa) ! p mean ^kappa |
|---|
| 176 | pm2 = pmk2**(1 / rkappa) ! p mean |
|---|
| 177 | IF(pm2>ptph) go to 110 ! doesn't happen |
|---|
| 178 | IF(pm2<p2km) go to 888 ! ptropo is valid |
|---|
| 179 | |
|---|
| 180 | a2 = (t(jj - 1) - t(jj)) ! a |
|---|
| 181 | a2 = a2 / (p(jj - 1)**rkappa - p(jj)**rkappa) |
|---|
| 182 | b2 = t(jj) - (a2 * p(jj)**rkappa) ! b |
|---|
| 183 | tm2 = a2 * pmk2 + b2 ! T mean |
|---|
| 184 | dtdp2 = a2 * rkappa * (pm2**(rkappa - 1)) ! dt/dp |
|---|
| 185 | dtdz2 = faktor * dtdp2 * pm2 / tm2 |
|---|
| 186 | asum = asum + dtdz2 |
|---|
| 187 | icount = icount + 1 |
|---|
| 188 | aquer = asum / float(icount) ! dt/dz mean |
|---|
| 189 | |
|---|
| 190 | ! discard ptropo ? |
|---|
| 191 | IF (aquer<=gamma) go to 999 ! dt/dz above < gamma |
|---|
| 192 | |
|---|
| 193 | 110 continue |
|---|
| 194 | enddo ! test next level |
|---|
| 195 | |
|---|
| 196 | 888 continue ! ptph is valid |
|---|
| 197 | ptrp = ptph |
|---|
| 198 | ttrp = ttph |
|---|
| 199 | |
|---|
| 200 | ! now calculate height of tropopause by integrating hypsometric equation |
|---|
| 201 | ! from ps to ptrp |
|---|
| 202 | ! linearly interpolate in p (results are identical to p^kappa) |
|---|
| 203 | |
|---|
| 204 | jj = LEVEL ! bottom |
|---|
| 205 | DO while ((P(jj)>PS) .OR. (T(jj)<100)) ! T must be valid too |
|---|
| 206 | jj = jj - 1 |
|---|
| 207 | END DO |
|---|
| 208 | |
|---|
| 209 | DLNP = log(PS / P(jj)) ! from surface pressure |
|---|
| 210 | TM = T(jj) ! take TM of lowest level (better: extrapolate) |
|---|
| 211 | TDLNP = TM * DLNP |
|---|
| 212 | |
|---|
| 213 | DO while ((JJ>=2) .AND. (PTRP<P(jj - 1))) |
|---|
| 214 | DLNP = log(P(jj) / P(jj - 1)) |
|---|
| 215 | TM = 0.5 * (T(jj) + T(jj - 1)) |
|---|
| 216 | TDLNP = TDLNP + TM * DLNP |
|---|
| 217 | JJ = JJ - 1 |
|---|
| 218 | END DO |
|---|
| 219 | |
|---|
| 220 | DLNP = log(P(jj) / PTRP) ! up to tropopause pressure |
|---|
| 221 | TM = 0.5 * (T(jj) + TTRP) ! use TTRP to get TM of this level |
|---|
| 222 | TDLNP = TDLNP + TM * DLNP |
|---|
| 223 | |
|---|
| 224 | ZTRP = ZS + TDLNP * RD / RG |
|---|
| 225 | |
|---|
| 226 | !!if (ZTRP .lt. 0) THEN |
|---|
| 227 | !! PRINT*,'ZTRP=',ZTRP |
|---|
| 228 | !! PRINT*,'PS=',PS |
|---|
| 229 | !! PRINT*,'P=',P |
|---|
| 230 | !! PRINT*,'T=',T |
|---|
| 231 | !! PRINT*,'ZS=',ZS |
|---|
| 232 | !! stop |
|---|
| 233 | !!endif |
|---|
| 234 | |
|---|
| 235 | return |
|---|
| 236 | |
|---|
| 237 | 999 continue ! continue search at next higher level |
|---|
| 238 | tm0 = tm |
|---|
| 239 | pm0 = pm |
|---|
| 240 | pmk0 = pmk |
|---|
| 241 | dtdz0 = dtdz |
|---|
| 242 | a0 = a |
|---|
| 243 | b0 = b |
|---|
| 244 | |
|---|
| 245 | END DO |
|---|
| 246 | |
|---|
| 247 | ! no tropopouse found |
|---|
| 248 | |
|---|
| 249 | END SUBROUTINE twmo |
|---|
