GCC Code Coverage Report


Directory: ./
File: phys/tropopause_m.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 0 47 0.0%
Branches: 0 116 0.0%

Line Branch Exec Source
1 MODULE tropopause_m
2
3 IMPLICIT NONE
4 PRIVATE
5 PUBLIC :: dyn_tropopause
6
7 CONTAINS
8
9 !-------------------------------------------------------------------------------
10 !
11 FUNCTION dyn_tropopause(t, ts, paprs, pplay, rot, itrop, thet0, pvor0)
12 !
13 !-------------------------------------------------------------------------------
14 USE assert_m, ONLY: assert
15 USE assert_eq_m, ONLY: assert_eq
16 USE dimphy, ONLY: klon, klev
17 USE geometry_mod, ONLY: latitude_deg, longitude_deg
18 USE vertical_layers_mod, ONLY: aps, bps, preff
19
20 !-------------------------------------------------------------------------------
21 ! Arguments:
22 REAL :: dyn_tropopause(klon) !--- Pressure at tropopause
23 REAL, INTENT(IN) :: t(:,:) !--- Cells-centers temperature
24 REAL, INTENT(IN) :: ts(:) !--- Surface temperature
25 REAL, INTENT(IN) :: paprs(:,:) !--- Cells-edges pressure
26 REAL, INTENT(IN) :: pplay(:,:) !--- Cells-centers pressure
27 REAL, INTENT(IN) :: rot(:,:) !--- Cells-centers relative vorticity
28 INTEGER, INTENT(OUT), OPTIONAL :: itrop(klon) !--- Last tropospheric layer idx
29 REAL, INTENT(IN), OPTIONAL :: thet0, pvor0
30 !-------------------------------------------------------------------------------
31 ! Local variables:
32 include "YOMCST.h"
33 REAL, PARAMETER :: DynPTrMin =8.E+3 !--- Thresholds for minimum and maximum
34 REAL, PARAMETER :: DynPTrMax =4.E+4 ! dynamical tropopause pressure (Pa).
35 CHARACTER(LEN=80) :: sub
36 INTEGER :: i, k, kb, kt, kp, ib, ie, nw
37 REAL :: al, th0, pv0
38 REAL, DIMENSION(klon,klev) :: tpot_cen, tpot_edg, pvor_cen
39 REAL, PARAMETER :: sg0=0.75 !--- Start level for PV=cte search loop
40 INTEGER, PARAMETER :: nadj=3 !--- Adjacent levs nb for thresholds detection
41 REAL, PARAMETER :: w(5)=[0.1,0.25,0.3,0.25,0.1] !--- Vertical smoothing
42 INTEGER, SAVE :: k0
43 LOGICAL, SAVE :: first=.TRUE.
44 !$OMP THREADPRIVATE(k0,first)
45 !-------------------------------------------------------------------------------
46 sub='dyn_tropopause'
47 CALL assert(SIZE(t ,1)==klon, TRIM(sub)//" t klon")
48 CALL assert(SIZE(t ,2)==klev, TRIM(sub)//" t klev")
49 CALL assert(SIZE(ts,1)==klon, TRIM(sub)//" ts klon")
50 CALL assert(SHAPE(paprs)==[klon,klev+1],TRIM(sub)//" paprs shape")
51 CALL assert(SHAPE(pplay)==[klon,klev ],TRIM(sub)//" pplay shape")
52 CALL assert(SHAPE(rot) ==[klon,klev ],TRIM(sub)//" rot shape")
53
54 !--- DEFAULT THRESHOLDS
55 th0=380.; IF(PRESENT(thet0)) th0=thet0 !--- In kelvins
56 pv0= 2.; IF(PRESENT(pvor0)) pv0=pvor0 !--- In PVU
57 IF(first) THEN
58 DO k0=1,klev; IF(aps(k0)/preff+bps(k0)<sg0) EXIT; END DO; first=.FALSE.
59 END IF
60
61 !--- POTENTIAL TEMPERATURE AT CELLS CENTERS AND INTERFACES
62 DO i = 1,klon
63 tpot_cen(i,1) = t(i,1)*(preff/pplay(i,1))**RKAPPA
64 tpot_edg(i,1) = ts(i) *(preff/paprs(i,1))**RKAPPA
65 DO k=2,klev
66 al = LOG(pplay(i,k-1)/paprs(i,k))/LOG(pplay(i,k-1)/pplay(i,k))
67 tpot_cen(i,k) = t(i,k) *(preff/pplay(i,k))**RKAPPA
68 tpot_edg(i,k) = (t(i,k-1)+al*(t(i,k)-t(i,k-1)))*(preff/paprs(i,k))**RKAPPA
69 !--- FORCE QUANTITIES TO BE GROWING
70 IF(tpot_edg(i,k)<tpot_edg(i,k-1)) tpot_edg(i,k)=tpot_edg(i,k-1)+1.E-5
71 IF(tpot_cen(i,k)<tpot_cen(i,k-1)) tpot_cen(i,k)=tpot_cen(i,k-1)+1.E-5
72 END DO
73 !--- VERTICAL SMOOTHING
74 tpot_cen(i,:)=smooth(tpot_cen(i,:),w)
75 tpot_edg(i,:)=smooth(tpot_edg(i,:),w)
76 END DO
77
78 !--- ERTEL POTENTIAL VORTICITY AT CELLS CENTERS (except in top layer)
79 DO i = 1, klon
80 DO k= 1, klev-1
81 pvor_cen(i,k)=-1.E6*RG*(rot(i,k)+2.*ROMEGA*SIN(latitude_deg(i)*RPI/180.))&
82 * (tpot_edg(i,k+1)-tpot_edg(i,k)) / (paprs(i,k+1)-paprs(i,k))
83 END DO
84 !--- VERTICAL SMOOTHING
85 pvor_cen(i,1:klev-1)=smooth(pvor_cen(i,1:klev-1),w)
86 END DO
87
88 !--- LOCATE TROPOPAUSE: LOWEST POINT BETWEEN THETA=380K AND PV=2PVU SURFACES.
89 DO i = 1, klon
90 !--- UPPER TROPOPAUSE: |PV|=2PVU POINT STARTING FROM TOP
91 DO kt=klev-1,1,-1; IF(ALL(ABS(pvor_cen(i,kt-nadj:kt))<=pv0)) EXIT; END DO
92 !--- LOWER TROPOPAUSE: |PV|=2PVU POINT STARTING FROM BOTTOM
93 DO kb=k0,klev-1; IF(ALL(ABS(pvor_cen(i,kb:kb+nadj))> pv0)) EXIT; END DO; kb=kb-1
94 !--- ISO-THETA POINT: THETA=380K STARTING FROM TOP
95 DO kp=klev-1,1,-1; IF(ALL(ABS(tpot_cen(i,kp-nadj:kp))<=th0)) EXIT; END DO
96 !--- CHOOSE BETWEEN LOWER AND UPPER TROPOPAUSE
97 IF(2*COUNT(ABS(pvor_cen(i,kb:kt))>pv0)>kt-kb+1) kt=kb
98 !--- PV-DEFINED TROPOPAUSE
99 al = (ABS(pvor_cen(i,kt+1))-pv0)/ABS(pvor_cen(i,kt+1)-pvor_cen(i,kt))
100 dyn_tropopause(i) = pplay(i,kt+1)*(pplay(i,kt)/pplay(i,kt+1))**al
101 !--- THETA=380K IN THE TROPICAL REGION
102 al = (tpot_cen(i,kp+1)-th0)/(tpot_cen(i,kp+1)-tpot_cen(i,kp))
103 dyn_tropopause(i) = MAX( pplay(i,kp+1)*(pplay(i,kp)/pplay(i,kp+1))**al, &
104 dyn_tropopause(i) )
105 !--- UNREALISTIC VALUES DETECTION
106 IF(dyn_tropopause(i)<DynPTrMin.OR.dyn_tropopause(i)>DynPTrMax) THEN
107 dyn_tropopause(i)=MIN(MAX(dyn_tropopause(i),DynPTrMax),DynPTrMin)
108 DO kt=1,klev-1; IF(pplay(i,kt+1)>dyn_tropopause(i)) EXIT; END DO; kp=kt
109 END IF
110 !--- LAST TROPOSPHERIC LAYER INDEX NEEDED
111 IF(PRESENT(itrop)) itrop(i)=MAX(kt,kp)
112 END DO
113
114 END FUNCTION dyn_tropopause
115
116
117 !-------------------------------------------------------------------------------
118 !
119 FUNCTION smooth(v,w)
120 !
121 !-------------------------------------------------------------------------------
122 ! Arguments:
123 REAL, INTENT(IN) :: v(:), w(:)
124 REAL, DIMENSION(SIZE(v)) :: smooth
125 !-------------------------------------------------------------------------------
126 ! Local variables:
127 INTEGER :: nv, nw, k, kb, ke, lb, le
128 !-------------------------------------------------------------------------------
129 nv=SIZE(v); nw=(SIZE(w)-1)/2
130 DO k=1,nv
131 kb=MAX(k-nw,1 ); lb=MAX(2+nw -k,1)
132 ke=MIN(k+nw,nv); le=MIN(1+nw+nv-k,1+2*nw)
133 smooth(k)=SUM(v(kb:ke)*w(lb:le))/SUM(w(lb:le))
134 END DO
135
136 END FUNCTION smooth
137 !
138 !-------------------------------------------------------------------------------
139
140 END MODULE tropopause_m
141