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

Last change on this file since 2534 was 2527, checked in by oboucher, 8 years ago

Adding a routine to compute a stratospheric mask
Will be needed later on for stratospheric aerosols

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