Context Navigation

tropopause_m.f90 @ 5918

Last change on this file since 5918 was 5650, checked in by dcugnet, 7 months ago
Bug fix: Coriolis parameters should be computed for each task and not broadcasted from the master.
File size: 17.8 KB

Line
1	MODULE tropopause_m
2
3	IMPLICIT NONE
4
5	PRIVATE
6
7	PUBLIC :: dyn_tropopause
8
9	REAL, PARAMETER :: DynPTrMin = 8.E+3 !--- Dyn tropopause pressures < DynPTrMin are set to DynPTrMin (Pa)
10	REAL, PARAMETER :: DynPTrMax = 4.E+4 !--- Dyn tropopause pressures > DynPTrMax are set to DynPTrMax (Pa)
11	REAL, PARAMETER :: theta0 = 380. !--- Default threshold for theta-defined tropopause (K)
12	REAL, PARAMETER :: pVort0 = 2.0 !--- Default threshold for PV-defined tropopause (PVU)
13	REAL, PARAMETER :: sg0 = 0.75 !--- Bottom->top PV=pv0e search loop starts at sigma=sg0 level
14	INTEGER, PARAMETER :: nadj = 3 !--- Threshold must be exceeded on nadj adjacent levels
15	INTEGER, PARAMETER :: ns = 2 !--- Number of neighbours used each side for vertical smoothing
16
17	CONTAINS
18
19	!===============================================================================================================================
20	FUNCTION dyn_tropopause(t, ts, paprs, pplay, rot, itrop, thet0, potV0) RESULT(pTrop)
21	USE assert_m, ONLY: assert
22	USE assert_eq_m, ONLY: assert_eq
23	USE dimphy, ONLY: klon, klev
24	USE geometry_mod, ONLY: latitude
25	USE strings_mod, ONLY: maxlen
26	USE yomcst_mod_h, ONLY: ROMEGA, RKAPPA, RG
27	USE vertical_layers_mod, ONLY: aps, bps, preff
28	USE lmdz_reprobus_wrappers, ONLY: itroprep
29	USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_REPROBUS
30	USE mod_phys_lmdz_para, ONLY: is_master
31	USE mod_phys_lmdz_transfert_para, ONLY : bcast
32	IMPLICIT NONE
33	REAL :: pTrop(klon) !--- Pressure at dynamical tropopause (Pa)
34	REAL, INTENT(IN) :: t(:,:) !--- Temperature at layers centers (K)
35	REAL, INTENT(IN) :: ts(:) !--- Temperature on surface layer interface (K)
36	REAL, INTENT(IN) :: paprs(:,:) !--- Pressure at layers interfaces (Pa)
37	REAL, INTENT(IN) :: pplay(:,:) !--- Pressure at layers centers (Pa)
38	REAL, INTENT(IN) :: rot(:,:) !--- Relative vorticity at layers centers (s-1)
39	INTEGER, OPTIONAL, INTENT(OUT) :: itrop(klon) !--- Last tropospheric layer idx
40	REAL, OPTIONAL, INTENT(IN) :: thet0 !--- Potential temperature at the tropopause (tropical region) (K)
41	REAL, OPTIONAL, INTENT(IN) :: potV0 !--- Potential vorticity at the tropopause (rest of globe) (PVU)
42	!------------------------------------------------------------------------------------------------------------------------------
43	CHARACTER(LEN=maxlen) :: modname !--- Current routine name
44	REAL :: Temp_edg(klon,klev) !--- Regular temperature at layers interfaces (except last one)(K)
45	REAL :: potTemp_edg(klon,klev) !--- Potential temperature at layers interfaces (except last one)(K)
46	REAL :: potTemp_cen(klon,klev) !--- Potential temperature at layers centers (K)
47	REAL :: potVort_cen(klon,klev) !--- Potential vorticity at layers centers (K)
48	REAL :: p_th0(klon) !--- Pressures at theta=380K (Pa)
49	REAL :: p_pv0(klon) !--- Pressures at PV=2PVU (Pa)
50	REAL :: al, th0, pv0 !--- Interpolation coefficient + potential temp. and PV thresholds
51	INTEGER :: i, k, kb, kt, kp, ib, ie, nw, n
52	INTEGER :: ith(klon) !--- Indices of first TH=380K layers (top -> bottom search)
53	INTEGER :: ipv(klon) !--- Indices of first PV=2PVU layers (top -> bottom search)
54	INTEGER :: ipv0(klon) !--- Indices of first PV=2PVU layers (bottom -> top search)
55	INTEGER :: ncons(klon) !--- Number of consecutive matching values found in vertical loops
56	INTEGER :: itr(klon) !--- Index of last layer with a center pressure lower than pTrop
57	INTEGER :: co(2*ns+1) !--- Binomial coefficients used compute smoothing weights "w(:,:)"
58	INTEGER, SAVE :: k0 !--- Start index (sigma=sg0) for 2PVU bottom->top search loop
59	REAL, ALLOCATABLE, SAVE :: fac(:) !--- Coriolis parameter: 2ROMEGASIN(cells centers latitudes) (s-1)
60	REAL, ALLOCATABLE, SAVE :: w(:,:) !--- Coefficients for vertical smoothing froutine "smooth"
61	LOGICAL, SAVE :: lFirst = .TRUE.
62	!$OMP THREADPRIVATE(k0, fac, w, lFirst)
63	!------------------------------------------------------------------------------------------------------------------------------
64	modname = 'dyn_tropopause'
65	CALL assert(SIZE(t, DIM=1) == klon, TRIM(modname)//" t klon")
66	CALL assert(SIZE(t, DIM=2) == klev, TRIM(modname)//" t klev")
67	CALL assert(SIZE(ts, DIM=1) == klon, TRIM(modname)//" ts klon")
68	CALL assert(SHAPE(paprs) == [klon, klev+1], TRIM(modname)//" paprs shape")
69	CALL assert(SHAPE(pplay) == [klon, klev ], TRIM(modname)//" pplay shape")
70	CALL assert(SHAPE(rot) == [klon, klev ], TRIM(modname)//" rot shape")
71
72	!--- MODIFY THE THRESHOLDS FOR THE DYNAMICAL TROPOPAUSE DEFINITION IN CASE THE CORRESPONDING OPTIONAL ARGUMENTS ARE USED
73	th0 = theta0; IF(PRESENT(thet0)) th0 = thet0 !--- Potential temperature at the tropopause (tropical region) (K)
74	pv0 = pVort0; IF(PRESENT(potV0)) pv0 = potV0 !--- Potential vorticity at the tropopause (rest of globe) (PVU)
75
76	IF(lFirst) THEN
77	ALLOCATE(fac(klon), w(ns+1, ns+1))
78
79	!--- COMPUTE THE CORIOLIS PARAMETER FOR PV ALCULATION ROUTINE "potentialVorticity"
80	DO i = 1, klon
81	fac(i) = 2. * ROMEGA * SIN(latitude(i))
82	END DO
83	!$OMP BARRIER
84
85	IF(is_master) THEN
86
87	!--- GET THE INDEX "k0" OF THE FIRST LOWER INTERFACE LAYER WITH SIGMA COORDINATE LOWER THAN "sg0"
88	!--- NOTE: "k0" DEPENDS ON VERTICAL DISCRETIZATION ONLY (VIA HYBRID COEFFS aps, bps) AND IS NOT SIMULATION-DEPENDENT
89	DO k0 = 1, klev; IF( aps(k0) / preff + bps(k0) < sg0 ) EXIT; END DO !--- START INDEX FOR BOTTOM->TOP PV SEARCH LOOP
90
91	!--- COMPUTE THE WEIGHTS FOR THE VERTICAL SMOOTHING ROUTINE "smooth"
92	co(:) = 0; w(:, :) = 0.
93	co(1) = 1; w(1, 1) = 1.
94	DO i = 1, ns
95	co(2:2ns+1) = co(2:2ns+1) + co(1:2*ns) !--- C(n+1,p+1) = C(n,p+1) + C(n,p)
96	co(2:2ns+1) = co(2:2ns+1) + co(1:2*ns) !--- C(n+1,p+1) = C(n,p+1) + C(n,p) AGAIN
97	w(i+1, 1:i+1) = REAL(co(i+1:2i+1))/REAL(SUM(co(i+1:2i+1)))
98	END DO
99
100	lFirst=.FALSE.
101	END IF
102	CALL bcast(k0)
103	CALL bcast(w)
104	CALL bcast(lFirst)
105	END IF
106
107	!=== DETERMINE THE PRESSURE AT WHICH THETA = th0 ============================================================================
108	CALL potentialTemperature(pplay, t, potTemp_cen) !--- POTENTIAL TEMPERATURE @ LAYERS CENTERS
109
110	!--- INDEX OF FIRST LAYERS WITH THETA<380K @ CENTER ON "nadj" CONSECUTIVE LAYERS
111	CALL getLayerIdx(potTemp_cen, th0, -1, nadj, ith) !--- FROM TOP TO BOTTOM
112
113	CALL getPressure(potTemp_cen, th0, ith, pplay, paprs, p_th0) !--- PRESSURE @ THETA = th0 SURFACE
114
115	!=== DETERMINE THE PRESSURE AT WHICH PV = pv0 ===============================================================================
116	CALL cen2edg(t, ts, pplay, paprs(:,1:klev), temp_edg) !--- TEMP @ LAYERS INTERFACES (EXCEPT LAST ONE)
117
118	CALL potentialTemperature (paprs(:,1:klev), temp_edg, potTemp_edg) !--- TPOT @ LAYERS INTERFACES (EXCEPT LAST ONE)
119
120	CALL potentialVorticity(rot, potTemp_edg, paprs(:,1:klev), potVort_cen) !--- ERTEL POTENTIAL VORTICITY @ LAYERS CENTERS
121
122	!--- INDEX OF FIRST LAYERS WITH PV<=2PVU @ CENTER ON "nadj" CONSECUTIVE LAYERS
123	CALL getLayerIdx(potVort_cen, pv0, -1, nadj, ipv) !--- FROM TOP TO BOTTOM
124	CALL getLayerIdx(potVort_cen, pv0, k0, nadj, ipv0) !--- FROM LAYER @ sig=sig0 TO TOP
125	DO i = 1, klon; n = 0 !--- CHOOSE BETWEEN BOTTOM AND TOP INDEX
126	IF(ipv0(i) == k0-1 .OR. ipv0(i) > ipv(i)) CYCLE !--- ipv0 CAN'T BE USED
127	DO k = ipv0(i), ipv(i); IF(potVort_cen(i, k) > pv0) n = n+1; END DO !--- NUMBER OF POINTS WITH PV>2PVU
128	IF(2 * n >= ipv(i)-ipv0(i)+1) ipv(i) = ipv0(i) !--- MORE THAN 50% > pv0 => LOWER POINT KEPT
129	END DO
130
131	CALL getPressure(potVort_cen, pv0, ipv, pplay, paprs, p_pv0) !--- PRESSURE @ PV = pv0 SURFACE
132
133	!=== DETERMINE THE UNFILTERED DYNAMICAL TROPOPAUSE PRESSURE FIELD (LOWER POINT BETWEEN THETA=380K AND PV=2PVU) ==============
134	DO i = 1, klon
135	pTrop(i) = MAX(p_th0(i), p_pv0(i))
136	END DO
137
138	!=== FILTER THE PRESSURE FIELD: TOO HIGH AND TOO LOW VALUES ARE CLIPPED =====================================================
139	DO i = 1, klon
140	IF(pTrop(i) < DynPTrMin) pTrop(i) = DynPTrMin
141	IF(pTrop(i) > DynPTrMax) pTrop(i) = DynPTrMax
142	END DO
143
144	!=== LAST VERTICAL INDEX WITH A PRESSURE HIGHER THAN TROPOPAUSE PRESSURE ====================================================
145	IF(.NOT.(PRESENT(itrop) .OR. CPPKEY_REPROBUS)) RETURN
146	DO i = 1, klon
147	DO k = 1, klev
148	IF(pplay(i,k+1) <= pTrop(i)) EXIT
149	END DO
150	IF(PRESENT(itrop )) itrop(i) = k
151	IF(CPPKEY_REPROBUS) itroprep(i) = k
152	END DO
153
154	CONTAINS
155
156	!===============================================================================================================================
157	SUBROUTINE cen2edg(v_cen, v0_edg, p_cen, p_edg, v_edg)
158	IMPLICIT NONE
159	REAL, DIMENSION(klon, klev), INTENT(IN) :: v_cen, p_cen, p_edg
160	REAL, DIMENSION(klon), INTENT(IN) :: v0_edg
161	REAL, DIMENSION(klon, klev), INTENT(OUT) :: v_edg
162	INTEGER :: i, k
163	DO i = 1, klon
164	v_edg(i, 1) = v0_edg(i)
165	END DO
166	DO k = 2, klev
167	DO i = 1, klon
168	al = LOG(p_edg(i, k-1)/p_cen(i, k)) / LOG(p_cen(i, k-1)/p_cen(i, k)) !--- CENTER -> INTERFACE INTERPOLATION COEFF
169	v_edg(i, k) = v_cen(i, k-1) + al * (v_cen(i, k) - v_cen(i, k-1)) !--- FIELD AT LAYER INTERFACE
170	END DO
171	END DO
172	END SUBROUTINE cen2edg
173	!===============================================================================================================================
174	SUBROUTINE getPressure(v_cen, v0, ix, p_cen, p_int, pre_v0)
175	IMPLICIT NONE
176	REAL, INTENT(IN) :: v_cen(klon, klev), v0
177	INTEGER, INTENT(IN) :: ix(klon)
178	REAL, INTENT(IN) :: p_cen(klon, klev), p_int(klon, klev+1)
179	REAL, INTENT(OUT) :: pre_v0(klon)
180	REAL :: al
181	INTEGER :: i, k
182	DO i = 1, klon; k = ix(i)
183	IF(k == 0) THEN
184	pre_v0(i) = p_int(i,1)
185	ELSE IF(k == klev) THEN
186	pre_v0(i) = p_int(i,klev+1)
187	ELSE
188	al = (v0 - v_cen(i, k+1)) / (v_cen(i, k) - v_cen(i, k+1))
189	pre_v0(i) = p_cen(i, k+1) * (p_cen(i, k) / p_cen(i, k+1))**al
190	END IF
191	END DO
192	END SUBROUTINE getPressure
193	!===============================================================================================================================
194	SUBROUTINE getLayerIdx(v, v0, k0, nadj, ix)
195	! Purpose: Search for the index of the last layer ix(i) with a value v(i,k) lower than or equal to v0.
196	! At least nadj adjacent layers must satisfy the criterium (less - as much as possible - near top or bottom).
197	! The search is done from: * top to bottom if k0 < 0 (from k=klev to k=\|k0\|)
198	! * bottom to top if k0 > 0 (from k=k0 to k=klev)
199	! - nominal case: k0 <= ix(i) < klev
200	! - special case: ix(i) == klev: ALL(v(i,k0:klev) <= v0)
201	! - special case: ix(i) == \|k0\|-1: ALL(v(i,k0:klev) > v0)
202	IMPLICIT NONE
203	REAL, INTENT(IN) :: v(klon, klev), v0
204	INTEGER, INTENT(IN) :: k0, nadj
205	INTEGER, INTENT(OUT) :: ix(klon)
206	INTEGER :: i, k, nc(klon)
207	nc(:) = 0
208	ix(:) = 0
209	IF(k0 < 0) THEN
210	!=== SEARCH FROM TOP TO BOTTOM: klev -> -k0
211	!--- ix(i) depends on nc(i), the number of adjacent layers with v(i,:) <= v0 (k is the index of the last tested layer)
212	!--- * nc(i) == nadj nominal case: enough matching values => ix(i) = k+nadj-1 (\|k0\|+nadj-1 <= k <= klev-nadj+1)
213	!--- particular case: all values are matching => ix(i) = klev (k = klev-nadj+1)
214	!--- * 0 < nc(i) < nadj bottom reached: nc<nadj matching values => ix(i) = k+nc(i)-1 (k = \|k0\|)
215	!--- * nc(i) == 0 bottom reached: no matching values => ix(i) = k (k = \|k0\|-1)
216	!--- So ix(i) = MAX(k, k+nc(i)-1) fits for each case.
217	DO k = klev, -1, -k0
218	DO i = 1, klon
219	IF(ix(i) /= 0) CYCLE !--- ADEQUATE LAYER ALREADY FOUND
220	nc(i) = nc(i) + 1
221	IF(ABS(v(i, k)) > v0) nc(i) = 0
222	IF(nc(i) /= nadj) CYCLE !--- nc<nadj ADJACENT LAYERS WITH v<=v0 FOUND
223	ix(i) = 1 !--- FAKE /=0 VALUE TO SKIP FOLLOWING ITERATIONS
224	END DO
225	END DO
226	DO i = 1, klon
227	ix(i) = MAX(k, k+nc(i)-1) !--- INDEX OF LOWEST LAYER WITH v<=v0
228	END DO
229	ELSE
230	!=== SEARCH FROM BOTTOM TO TOP: k0 -> klev
231	!--- ix(i) depends on nc(i), the number of adjacent layers with v(i,:) > v0 (k is the index of the last tested layer)
232	!--- * nc(i) == nadj nominal case: enough matching values => ix(i) = k-nadj ( k0 +nadj-1 <= k <= klev-nadj+1)
233	!--- particular case: all values are matching => ix(i) = k0-1 (k = k0+nadj-1)
234	!--- * 0 < nc(i) < nadj top reached: nc<nadj matching values => ix(i) = k-nc(i) (k = klev)
235	!--- * nc(i) == 0 top reached: no matching values => ix(i) = k (k = klev)
236	!--- So ix(i) = k-nc(i) fits for each case.
237	DO k = k0, klev
238	DO i = 1, klon
239	IF(ix(i) /= 0) CYCLE !--- ADEQUATE LAYER ALREADY FOUND
240	nc(i) = nc(i) + 1
241	IF(ABS(v(i, k)) <= v0) nc(i) = 0
242	IF(nc(i) /= nadj) CYCLE !--- nc<nadj ADJACENT LAYERS WITH v<=v0 FOUND
243	ix(i) = 1 !--- FAKE /=0 VALUE TO SKIP FOLLOWING ITERATIONS
244	END DO
245	END DO
246	DO i = 1, klon
247	ix(i) = k-nc(i) !--- INDEX OF LOWEST LAYER WITH v<=v0
248	END DO
249	END IF
250	END SUBROUTINE getLayerIdx
251	!===============================================================================================================================
252	SUBROUTINE potentialTemperature(pre, temp, tPot)
253	IMPLICIT NONE
254	REAL, DIMENSION(:, :), INTENT(IN) :: pre, temp
255	REAL, DIMENSION(SIZE(pre, 1), SIZE(pre, 2)), INTENT(OUT) :: tPot
256	REAL, ALLOCATABLE :: tmp(:,:)
257	CHARACTER(LEN=maxlen) :: modname
258	INTEGER :: i, k, ni, nk
259	modname = 'potentialTemperature'
260	ni = SIZE(pre, 1)
261	nk = SIZE(pre, 2)
262	CALL assert(SIZE(temp, DIM=1) == ni, TRIM(modname)//" SIZE(temp,1) SIZE(pre,1)")
263	CALL assert(SIZE(temp, DIM=2) == nk, TRIM(modname)//" SIZE(temp,2) SIZE(pre,2)")
264	ALLOCATE(tmp(ni, nk))
265	DO k = 1, nk !--- COMPUTE RAW FIELD
266	DO i = 1, ni
267	tmp(i, k) = temp(i, k) * (100000. / pre(i, k))**RKAPPA
268	END DO
269	END DO
270	DO k = 2, nk !--- ENSURE GROWING FIELD WITH ALTITUDE
271	DO i = 1, ni
272	IF(tmp(i, k)< tmp(i, k-1)) tmp(i, k) = tmp(i, k-1) + 1.E-5
273	END DO
274	END DO
275	CALL smooth(tmp, tPot) !--- FILTER THE FIELD
276	END SUBROUTINE potentialTemperature
277	!===============================================================================================================================
278	SUBROUTINE potentialVorticity(rot_cen, th_int, pint, pVor_cen)
279	IMPLICIT NONE
280	REAL, DIMENSION(klon, klev), INTENT(IN) :: rot_cen, th_int, pint
281	REAL, DIMENSION(klon, klev), INTENT(OUT) :: pVor_cen
282	REAL :: tmp(klon, klev)
283	INTEGER :: i, k, kp
284	DO k = 1, klev-1 !--- COMPUTE RAW FIELD
285	DO i = 1, klon
286	tmp(i, k) = -1.E6 * RG * (rot_cen(i, k) + fac(i)) * (th_int(i, k+1)-th_int(i, k)) / (pint(i, k+1)-pint(i, k))
287	END DO
288	END DO
289	DO i = 1, klon
290	tmp(i, klev) = tmp(i, klev-1)
291	END DO
292	CALL smooth(tmp, pVor_cen) !--- FILTER THE FIELD
293	END SUBROUTINE potentialVorticity
294	!===============================================================================================================================
295	SUBROUTINE smooth(v, vs)
296	! Purpose: Vertical smoothing of each profile v(i,:) using 2*ns+1 centered binomial weights (+/- ns points).
297	! Note: For levels near the bottom (k <= ns) or the top (k > klev-ns), a narrower set of weights (n<ns) is used.
298	! => in particular, first and last levels are left untouched.
299	IMPLICIT NONE
300	REAL, INTENT(IN) :: v (klon, klev)
301	REAL, INTENT(OUT) :: vs(klon, klev)
302	INTEGER :: i, j, k
303	vs(:, :) = 0.
304	DO k = 1, klev
305	n = MIN(k-1, klev-k, ns)
306	DO j = k-n, k+n
307	DO i = 1, klon
308	vs(i, k) = vs(i, k) + v(i, j) * w(n+1, 1+ABS(j-k))
309	END DO
310	END DO
311	END DO
312	END SUBROUTINE smooth
313
314	END FUNCTION dyn_tropopause
315
316	END MODULE tropopause_m

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format