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