source: LMDZ5/trunk/libf/phylmd/stratosphere_mask.F90 @ 5135

Last change on this file since 5135 was 2992, checked in by oboucher, 7 years ago

Adding p, z, and t of tropopause in output files.
Stratosphere_mask is streamlined a bit.
This includes a treatment of missing values.

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