source: LMDZ6/branches/IPSLCM6.0.14/libf/phylmd/stratosphere_mask.F90 @ 5446

Last change on this file since 5446 was 3120, checked in by Laurent Fairhead, 7 years ago

Inclusion of bugfix 3119 from trunk
OB

  • Property svn:keywords set to Id
File size: 7.8 KB
RevLine 
[2772]1!
2! $Id: stratosphere_mask.F90 3120 2017-12-09 17:49:01Z fhourdin $
3!
[2992]4SUBROUTINE stratosphere_mask(missing_val, t_seri, pplay, xlat)
[2527]5
6!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
7!
8! determination of tropopause height and temperature from gridded temperature data
9!
[3120]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)
[3120]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
[2527]31!
32!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
33
34USE dimphy
[2536]35USE phys_local_var_mod, ONLY: stratomask
[2992]36USE phys_local_var_mod, ONLY: p_tropopause, z_tropopause, t_tropopause
[2772]37USE print_control_mod, ONLY: lunout, prt_level
[2527]38
39IMPLICIT NONE
40
[2992]41REAL, INTENT(IN)                       :: missing_val ! missing value, also XIOS
[2527]42REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri  ! Temperature
43REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay   ! pression pour le mileu de chaque couche (en Pa)
44REAL,DIMENSION(klon),INTENT(IN)        :: xlat    ! latitudes pour chaque point
45
46REAL, PARAMETER                        :: plimu=45000.
47REAL, PARAMETER                        :: pliml=7500.
48REAL, PARAMETER                        :: gamma=-0.002
49LOGICAL, PARAMETER                     :: dofill=.true.
[2992]50REAL,DIMENSION(klon)                   :: tp
[2527]51REAL,DIMENSION(klev)                   :: t, p
[2992]52INTEGER                                :: i, k, ifil
[2527]53REAL                                   :: ptrp, ttrp, ztrp, psrf, zsrf, pi
54
55pi     = 4.*ATAN(1.)
56
[2992]57!--computing tropopause
[2527]58DO i=1,klon
59  DO k=1,klev
60    t(k)=t_seri(i,klev+1-k)
61    p(k)=pplay(i,klev+1-k)
62  ENDDO
63  psrf=pplay(i,1)
64  zsrf=0.0
[2992]65  call twmo(missing_val, klev, t, p, psrf, zsrf, plimu, pliml, gamma, ptrp, ttrp, ztrp)
[2527]66  tp(i)=ptrp
[2992]67  p_tropopause(i)=ptrp
68  z_tropopause(i)=ztrp
69  t_tropopause(i)=ttrp
[2527]70ENDDO
71
[2992]72!--filling holes in tp but not in p_tropopause
[2527]73IF (dofill) THEN
74  ifil=0
75  DO i=1,klon
[2992]76  IF (ABS(tp(i)/missing_val-1.0).LT.0.01) THEN
[2527]77    !set missing values to very simple profile (neighbour averaging too expensive in LMDZ)
78    tp(i)=50000.-20000.*cos(xlat(i)/360.*2.*pi)
79    ifil=ifil+1
80  ENDIF
81  ENDDO
82!
83ENDIF
84!
85DO i=1, klon
86DO k=1, klev
87  IF (pplay(i,k).LT.tp(i)) THEN
[2536]88    stratomask(i,k)=1.0
[2527]89  ELSE
[2536]90    stratomask(i,k)=0.0
[2527]91  ENDIF
92ENDDO
93ENDDO
94
[2992]95IF (ifil.GT.0 .AND. prt_level >5) THEN
[2772]96  write(lunout,*)'Tropopause: number of undetermined values =', ifil
[2527]97ENDIF
98
99RETURN
100END SUBROUTINE stratosphere_mask
101
102!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
103! twmo
104!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
105
[2992]106subroutine twmo(missing_val, level, t, p, ps, zs, plimu, pliml, gamma, ptrp, ttrp, ztrp)
[2527]107
[3120]108! reference: Reichler, T., M. Dameris, and R. Sausen (GRL, 10.1029/2003GL018240, 2003)
109
[2527]110implicit none
111
112include "YOMCST.h"
113
114integer,intent(in)              :: level
[2992]115real,intent(in)                 :: missing_val
[2527]116real,intent(in),dimension(level):: t, p
117real,intent(in)                 :: plimu, pliml, gamma, ps, zs
118real,intent(out)                :: ptrp, ttrp, ztrp
119
120real,parameter                  :: deltaz = 2000.0
121
122real     :: faktor
123real     :: pmk, pm, a, b, tm, dtdp, dtdz, dlnp, tdlnp
124real     :: ag, bg, ptph, ttph, a0, b0
125real     :: pm0, tm0, pmk0, dtdz0
126real     :: p2km, asum, aquer
127real     :: pmk2, pm2, a2, b2, tm2, dtdp2, dtdz2
128integer  :: icount, jj, j
129
[2992]130ptrp=missing_val
131ttrp=missing_val
132ztrp=missing_val
[2527]133
[3120]134faktor = -RG/RD
[2527]135
136do 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.eq.level)     go to 999                 ! no, start level, initialize first
149   if (dtdz.le.gamma)  go to 999                 ! no, dt/dz < -2 K/km
150   if (dtdz0.gt.gamma) go to 999                 ! no, dt/dz below > -2 K/km
151   if (pm.gt.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.lt.pliml) go to 999
165   if (ptph.gt.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.gt.ptph) go to 110                ! doesn't happen
178       if(pm2.lt.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.le.gamma) go to 999           ! dt/dz above < gamma
192
193110 continue
194    enddo                   ! test next level
195
196888 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
204jj = LEVEL                                       ! bottom
205do while ((P(jj).gt.PS) .or. (T(jj).lt.100))     ! T must be valid too
206  jj=jj-1
207enddo
208
209DLNP = log(PS/P(jj))                             ! from surface pressure
210TM = T(jj)                                       ! take TM of lowest level (better: extrapolate)
211TDLNP = TM*DLNP
212
213do while ( (JJ.ge.2) .and. (PTRP.lt.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
218enddo
219
220DLNP = log(P(jj)/PTRP)                           ! up to tropopause pressure
221TM = 0.5 * (T(jj) + TTRP)                        ! use TTRP to get TM of this level
222TDLNP = TDLNP + TM*DLNP
223
[3120]224ZTRP = ZS + TDLNP*RD/RG
[2527]225
[3120]226!!if (ZTRP .lt. 0) then
227!!  print*,ZTRP
228!!  print*,PS
229!!  print*,P
230!!  print*,T
231!!  print*,ZS
232!!  stop
233!!endif
[2527]234
235return
236
237999 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
245enddo
246
247! no tropopouse found
248return
249end subroutine twmo
Note: See TracBrowser for help on using the repository browser.