GCC Code Coverage Report

Directory: ./
File: phys/stratosphere_mask.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 80 85 94.1%
Branches: 36 48 75.0%
Line Branch Exec Source
1 !
2 ! $Id: stratosphere_mask.F90 3123 2017-12-13 13:16:38Z oboucher $
3 !
4 119400 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 240 REAL,DIMENSION(klon) :: tp
54 120 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
2/2
✓ Branch 0 taken 119280 times.
✓ Branch 1 taken 120 times.
 119400 DO i=1,klon
62
2/2
✓ Branch 0 taken 4651920 times.
✓ Branch 1 taken 119280 times.
 4771200 DO k=1,klev
63 4651920 t(k)=t_seri(i,klev+1-k)
64 4771200 p(k)=pplay(i,klev+1-k)
65 ENDDO
66 119280 psrf=pplay(i,1)
67 119280 zsrf=pphis(i)/RG !--altitude de la surface
68 119280 call twmo(missing_val, klev, t, p, psrf, zsrf, plimu, pliml, gamma, ptrp, ttrp, ztrp)
69 119280 tp(i)=ptrp
70 119280 p_tropopause(i)=ptrp
71 119280 z_tropopause(i)=ztrp
72 119400 t_tropopause(i)=ttrp
73 ENDDO
74
75 !--filling holes in tp but not in p_tropopause
76 IF (dofill) THEN
77 120 ifil=0
78
2/2
✓ Branch 0 taken 119280 times.
✓ Branch 1 taken 120 times.
 119400 DO i=1,klon
79
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 119280 times.
 119400 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
2/2
✓ Branch 0 taken 119280 times.
✓ Branch 1 taken 120 times.
 119400 DO i=1, klon
89
2/2
✓ Branch 0 taken 4651920 times.
✓ Branch 1 taken 119280 times.
 4771320 DO k=1, klev
90
2/2
✓ Branch 0 taken 2570659 times.
✓ Branch 1 taken 2081261 times.
 4771200 IF (pplay(i,k).LT.tp(i)) THEN
91 2570659 stratomask(i,k)=1.0
92 ELSE
93 2081261 stratomask(i,k)=0.0
94 ENDIF
95 ENDDO
96 ENDDO
97
98
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 120 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 120 IF (ifil.GT.0 .AND. prt_level >5) THEN
99 ✗ write(lunout,*)'Tropopause: number of undetermined values =', ifil
100 ENDIF
101
102 120 RETURN
103 END SUBROUTINE stratosphere_mask
104
105 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
106 ! twmo
107 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
108
109 2255666 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 119280 ptrp=missing_val
134 119280 ttrp=missing_val
135 119280 ztrp=missing_val
136
137 119280 faktor = -RG/RD
138
139
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2136386 times.
 2136386 do j=level,2,-1
140
141 ! dt/dz
142 2136386 pmk= 0.5 * (p(j-1)**rkappa+p(j)**rkappa) ! p**k at half level
143 2136386 pm = pmk**(1./rkappa) ! p at half level
144 2136386 a = (t(j-1)-t(j))/(p(j-1)**rkappa-p(j)**rkappa) ! dT/dp^k
145 2136386 b = t(j)-(a*p(j)**rkappa)
146 2136386 tm = a * pmk + b ! T at half level
147 2136386 dtdp = a * rkappa * pm**(rkappa-1.) ! dT/dp at half level
148 2136386 dtdz = faktor*dtdp*pm/tm ! dT/dz at half level
149
150 ! dt/dz valid?
151
2/2
✓ Branch 0 taken 119280 times.
✓ Branch 1 taken 2017106 times.
 2136386 if (j.eq.level) go to 999 ! no, start level, initialize first
152
2/2
✓ Branch 0 taken 1789875 times.
✓ Branch 1 taken 227231 times.
 2017106 if (dtdz.le.gamma) go to 999 ! no, dt/dz < -2 K/km
153
2/2
✓ Branch 0 taken 74659 times.
✓ Branch 1 taken 152572 times.
 227231 if (dtdz0.gt.gamma) go to 999 ! no, dt/dz below > -2 K/km
154
2/2
✓ Branch 0 taken 32064 times.
✓ Branch 1 taken 120508 times.
 152572 if (pm.gt.plimu) go to 999 ! no, pm too low
155
156 ! dtdz is valid, calculate tropopause pressure ptph
157 120508 ag = (dtdz-dtdz0) / (pmk-pmk0)
158 120508 bg = dtdz0 - (ag * pmk0)
159 120508 ptph = exp(log((gamma-bg)/ag)/rkappa)
160
161 ! calculate temperature at this ptph assuming linear gamma
162 ! interpolation
163 ttph = tm0
164 120508 ttph = ttph - (bg * log(pm0) + ag * (pm0**rkappa) /rkappa) / faktor*t(j)
165 120508 ttph = ttph + (bg * log(ptph) + ag * (ptph**rkappa)/rkappa) / faktor*t(j)
166
167
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 120508 times.
 120508 if (ptph.lt.pliml) go to 999
168
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 120508 times.
 120508 if (ptph.gt.plimu) go to 999
169
170 ! 2nd test: dtdz above 2 km must not exceed gamma
171 120508 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
1/2
✓ Branch 0 taken 367209 times.
✗ Branch 1 not taken.
 367209 do jj=j,2,-1
177
178 367209 pmk2 = .5 * (p(jj-1)**rkappa+p(jj)**rkappa) ! p mean ^kappa
179 367209 pm2 = pmk2**(1/rkappa) ! p mean
180
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 367209 times.
 367209 if(pm2.gt.ptph) go to 110 ! doesn't happen
181
2/2
✓ Branch 0 taken 119280 times.
✓ Branch 1 taken 247929 times.
 367209 if(pm2.lt.p2km) go to 888 ! ptropo is valid
182
183 247929 a2 = (t(jj-1)-t(jj)) ! a
184 247929 a2 = a2/(p(jj-1)**rkappa-p(jj)**rkappa)
185 247929 b2 = t(jj)-(a2*p(jj)**rkappa) ! b
186 247929 tm2 = a2 * pmk2 + b2 ! T mean
187 247929 dtdp2 = a2 * rkappa * (pm2**(rkappa-1)) ! dt/dp
188 247929 dtdz2 = faktor*dtdp2*pm2/tm2
189 247929 asum = asum+dtdz2
190 247929 icount = icount+1
191 247929 aquer = asum/float(icount) ! dt/dz mean
192
193 ! discard ptropo ?
194
2/2
✓ Branch 0 taken 246701 times.
✓ Branch 1 taken 1228 times.
 247929 if (aquer.le.gamma) go to 999 ! dt/dz above < gamma
195
196 119280 110 continue
197 enddo ! test next level
198
199 888 continue ! ptph is valid
200 119280 ptrp = ptph
201 119280 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
2/4
✓ Branch 0 taken 119280 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 119280 times.
 119280 do while ((P(jj).gt.PS) .or. (T(jj).lt.100)) ! T must be valid too
209 ✗ jj=jj-1
210 enddo
211
212 119280 DLNP = log(PS/P(jj)) ! from surface pressure
213 TM = T(jj) ! take TM of lowest level (better: extrapolate)
214 119280 TDLNP = TM*DLNP
215
216
3/4
✓ Branch 0 taken 2081261 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 119280 times.
✓ Branch 3 taken 1961981 times.
 2081261 do while ( (JJ.ge.2) .and. (PTRP.lt.P(jj-1)) )
217 1961981 DLNP = log(P(jj)/P(jj-1))
218 1961981 TM = 0.5 * (T(jj) + T(jj-1))
219 1961981 TDLNP = TDLNP + TM*DLNP
220 119280 JJ=JJ-1
221 enddo
222
223 119280 DLNP = log(P(jj)/PTRP) ! up to tropopause pressure
224 119280 TM = 0.5 * (T(jj) + TTRP) ! use TTRP to get TM of this level
225 119280 TDLNP = TDLNP + TM*DLNP
226
227 119280 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 2136386 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
253