1 |
|
|
! |
2 |
|
|
! $Id: stratosphere_mask.F90 3123 2017-12-13 13:16:38Z oboucher $ |
3 |
|
|
! |
4 |
|
71640 |
SUBROUTINE stratosphere_mask(missing_val, pphis, t_seri, pplay, xlat) |
5 |
|
|
|
6 |
|
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
7 |
|
|
! |
8 |
|
|
! determination of tropopause height and temperature from gridded temperature data |
9 |
|
|
! |
10 |
|
|
! reference: Reichler, T., M. Dameris, and R. Sausen (GRL, 10.1029/2003GL018240, 2003) |
11 |
|
|
! modified: 6/28/06 tjr |
12 |
|
|
! adapted to LMDZ by C. Kleinschmitt (2016-02-15) |
13 |
|
|
! committed to LMDz by O. Boucher (2016) with a mistake |
14 |
|
|
! mistake corrected by O. Boucher (2017-12-11) |
15 |
|
|
! |
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 |
25 |
|
|
! |
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 |
31 |
|
|
! |
32 |
|
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
33 |
|
|
|
34 |
|
|
USE dimphy |
35 |
|
|
USE phys_local_var_mod, ONLY: stratomask |
36 |
|
|
USE phys_local_var_mod, ONLY: p_tropopause, z_tropopause, t_tropopause |
37 |
|
|
USE print_control_mod, ONLY: lunout, prt_level |
38 |
|
|
|
39 |
|
|
IMPLICIT NONE |
40 |
|
|
|
41 |
|
|
INCLUDE "YOMCST.h" |
42 |
|
|
|
43 |
|
|
REAL, INTENT(IN) :: missing_val ! missing value, also XIOS |
44 |
|
|
REAL,DIMENSION(klon),INTENT(IN) :: pphis ! Geopotentiel de surface |
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. |
53 |
|
144 |
REAL,DIMENSION(klon) :: tp |
54 |
|
72 |
REAL,DIMENSION(klev) :: t, p |
55 |
|
|
INTEGER :: i, k, ifil |
56 |
|
|
REAL :: ptrp, ttrp, ztrp, psrf, zsrf, pi |
57 |
|
|
|
58 |
|
|
pi = 4.*ATAN(1.) |
59 |
|
|
|
60 |
|
|
!--computing tropopause |
61 |
✓✓ |
71640 |
DO i=1,klon |
62 |
✓✓ |
2862720 |
DO k=1,klev |
63 |
|
2791152 |
t(k)=t_seri(i,klev+1-k) |
64 |
|
2862720 |
p(k)=pplay(i,klev+1-k) |
65 |
|
|
ENDDO |
66 |
|
71568 |
psrf=pplay(i,1) |
67 |
|
71568 |
zsrf=pphis(i)/RG !--altitude de la surface |
68 |
|
71568 |
call twmo(missing_val, klev, t, p, psrf, zsrf, plimu, pliml, gamma, ptrp, ttrp, ztrp) |
69 |
|
71568 |
tp(i)=ptrp |
70 |
|
71568 |
p_tropopause(i)=ptrp |
71 |
|
71568 |
z_tropopause(i)=ztrp |
72 |
|
71640 |
t_tropopause(i)=ttrp |
73 |
|
|
ENDDO |
74 |
|
|
|
75 |
|
|
!--filling holes in tp but not in p_tropopause |
76 |
|
|
IF (dofill) THEN |
77 |
|
72 |
ifil=0 |
78 |
✓✓ |
71640 |
DO i=1,klon |
79 |
✗✓ |
71640 |
IF (ABS(tp(i)/missing_val-1.0).LT.0.01) THEN |
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 |
✓✓ |
71640 |
DO i=1, klon |
89 |
✓✓ |
2862792 |
DO k=1, klev |
90 |
✓✓ |
2862720 |
IF (pplay(i,k).LT.tp(i)) THEN |
91 |
|
1539327 |
stratomask(i,k)=1.0 |
92 |
|
|
ELSE |
93 |
|
1251825 |
stratomask(i,k)=0.0 |
94 |
|
|
ENDIF |
95 |
|
|
ENDDO |
96 |
|
|
ENDDO |
97 |
|
|
|
98 |
✗✓✗✗
|
72 |
IF (ifil.GT.0 .AND. prt_level >5) THEN |
99 |
|
|
write(lunout,*)'Tropopause: number of undetermined values =', ifil |
100 |
|
|
ENDIF |
101 |
|
|
|
102 |
|
72 |
RETURN |
103 |
|
|
END SUBROUTINE stratosphere_mask |
104 |
|
|
|
105 |
|
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
106 |
|
|
! twmo |
107 |
|
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
108 |
|
|
|
109 |
|
1356228 |
subroutine twmo(missing_val, level, t, p, ps, zs, plimu, pliml, gamma, ptrp, ttrp, ztrp) |
110 |
|
|
|
111 |
|
|
! reference: Reichler, T., M. Dameris, and R. Sausen (GRL, 10.1029/2003GL018240, 2003) |
112 |
|
|
|
113 |
|
|
implicit none |
114 |
|
|
|
115 |
|
|
include "YOMCST.h" |
116 |
|
|
|
117 |
|
|
integer,intent(in) :: level |
118 |
|
|
real,intent(in) :: missing_val |
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 |
|
|
|
133 |
|
71568 |
ptrp=missing_val |
134 |
|
71568 |
ttrp=missing_val |
135 |
|
71568 |
ztrp=missing_val |
136 |
|
|
|
137 |
|
71568 |
faktor = -RG/RD |
138 |
|
|
|
139 |
✗✓ |
1284660 |
do j=level,2,-1 |
140 |
|
|
|
141 |
|
|
! dt/dz |
142 |
|
1284660 |
pmk= 0.5 * (p(j-1)**rkappa+p(j)**rkappa) ! p**k at half level |
143 |
|
1284660 |
pm = pmk**(1./rkappa) ! p at half level |
144 |
|
1284660 |
a = (t(j-1)-t(j))/(p(j-1)**rkappa-p(j)**rkappa) ! dT/dp^k |
145 |
|
1284660 |
b = t(j)-(a*p(j)**rkappa) |
146 |
|
1284660 |
tm = a * pmk + b ! T at half level |
147 |
|
1284660 |
dtdp = a * rkappa * pm**(rkappa-1.) ! dT/dp at half level |
148 |
|
1284660 |
dtdz = faktor*dtdp*pm/tm ! dT/dz at half level |
149 |
|
|
|
150 |
|
|
! dt/dz valid? |
151 |
✓✓ |
1284660 |
if (j.eq.level) go to 999 ! no, start level, initialize first |
152 |
✓✓ |
1213092 |
if (dtdz.le.gamma) go to 999 ! no, dt/dz < -2 K/km |
153 |
✓✓ |
143177 |
if (dtdz0.gt.gamma) go to 999 ! no, dt/dz below > -2 K/km |
154 |
✓✓ |
96873 |
if (pm.gt.plimu) go to 999 ! no, pm too low |
155 |
|
|
|
156 |
|
|
! dtdz is valid, calculate tropopause pressure ptph |
157 |
|
72442 |
ag = (dtdz-dtdz0) / (pmk-pmk0) |
158 |
|
72442 |
bg = dtdz0 - (ag * pmk0) |
159 |
|
72442 |
ptph = exp(log((gamma-bg)/ag)/rkappa) |
160 |
|
|
|
161 |
|
|
! calculate temperature at this ptph assuming linear gamma |
162 |
|
|
! interpolation |
163 |
|
|
ttph = tm0 |
164 |
|
72442 |
ttph = ttph - (bg * log(pm0) + ag * (pm0**rkappa) /rkappa) / faktor*t(j) |
165 |
|
72442 |
ttph = ttph + (bg * log(ptph) + ag * (ptph**rkappa)/rkappa) / faktor*t(j) |
166 |
|
|
|
167 |
✗✓ |
72442 |
if (ptph.lt.pliml) go to 999 |
168 |
✗✓ |
72442 |
if (ptph.gt.plimu) go to 999 |
169 |
|
|
|
170 |
|
|
! 2nd test: dtdz above 2 km must not exceed gamma |
171 |
|
72442 |
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 |
✓✗ |
221025 |
do jj=j,2,-1 |
177 |
|
|
|
178 |
|
221025 |
pmk2 = .5 * (p(jj-1)**rkappa+p(jj)**rkappa) ! p mean ^kappa |
179 |
|
221025 |
pm2 = pmk2**(1/rkappa) ! p mean |
180 |
✗✓ |
221025 |
if(pm2.gt.ptph) go to 110 ! doesn't happen |
181 |
✓✓ |
221025 |
if(pm2.lt.p2km) go to 888 ! ptropo is valid |
182 |
|
|
|
183 |
|
149457 |
a2 = (t(jj-1)-t(jj)) ! a |
184 |
|
149457 |
a2 = a2/(p(jj-1)**rkappa-p(jj)**rkappa) |
185 |
|
149457 |
b2 = t(jj)-(a2*p(jj)**rkappa) ! b |
186 |
|
149457 |
tm2 = a2 * pmk2 + b2 ! T mean |
187 |
|
149457 |
dtdp2 = a2 * rkappa * (pm2**(rkappa-1)) ! dt/dp |
188 |
|
149457 |
dtdz2 = faktor*dtdp2*pm2/tm2 |
189 |
|
149457 |
asum = asum+dtdz2 |
190 |
|
149457 |
icount = icount+1 |
191 |
|
149457 |
aquer = asum/float(icount) ! dt/dz mean |
192 |
|
|
|
193 |
|
|
! discard ptropo ? |
194 |
✓✓ |
149457 |
if (aquer.le.gamma) go to 999 ! dt/dz above < gamma |
195 |
|
|
|
196 |
|
71568 |
110 continue |
197 |
|
|
enddo ! test next level |
198 |
|
|
|
199 |
|
|
888 continue ! ptph is valid |
200 |
|
71568 |
ptrp = ptph |
201 |
|
71568 |
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 |
✓✗✗✓
|
71568 |
do while ((P(jj).gt.PS) .or. (T(jj).lt.100)) ! T must be valid too |
209 |
|
|
jj=jj-1 |
210 |
|
|
enddo |
211 |
|
|
|
212 |
|
71568 |
DLNP = log(PS/P(jj)) ! from surface pressure |
213 |
|
|
TM = T(jj) ! take TM of lowest level (better: extrapolate) |
214 |
|
71568 |
TDLNP = TM*DLNP |
215 |
|
|
|
216 |
✓✗✓✓
|
1251825 |
do while ( (JJ.ge.2) .and. (PTRP.lt.P(jj-1)) ) |
217 |
|
1180257 |
DLNP = log(P(jj)/P(jj-1)) |
218 |
|
1180257 |
TM = 0.5 * (T(jj) + T(jj-1)) |
219 |
|
1180257 |
TDLNP = TDLNP + TM*DLNP |
220 |
|
71568 |
JJ=JJ-1 |
221 |
|
|
enddo |
222 |
|
|
|
223 |
|
71568 |
DLNP = log(P(jj)/PTRP) ! up to tropopause pressure |
224 |
|
71568 |
TM = 0.5 * (T(jj) + TTRP) ! use TTRP to get TM of this level |
225 |
|
71568 |
TDLNP = TDLNP + TM*DLNP |
226 |
|
|
|
227 |
|
71568 |
ZTRP = ZS + TDLNP*RD/RG |
228 |
|
|
|
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 |
237 |
|
|
|
238 |
|
1284660 |
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 |