Changeset 5607
- Timestamp:
- Apr 9, 2025, 9:47:14 PM (2 months ago)
- File:
-
- 1 edited
-
LMDZ6/trunk/libf/phylmd/tropopause_m.f90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
TabularUnified LMDZ6/trunk/libf/phylmd/tropopause_m.f90 ¶
r5578 r5607 1 1 MODULE tropopause_m 2 2 3 USE yomcst_mod_h4 IMPLICIT NONE 3 IMPLICIT NONE 4 5 5 PRIVATE 6 6 7 PUBLIC :: dyn_tropopause 7 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 8 17 CONTAINS 9 18 10 !------------------------------------------------------------------------------- 11 ! 12 FUNCTION dyn_tropopause(t, ts, paprs, pplay, rot, itrop, thet0, pvor0) 13 ! 14 !------------------------------------------------------------------------------- 19 !=============================================================================================================================== 20 FUNCTION dyn_tropopause(t, ts, paprs, pplay, rot, itrop, thet0, potV0) RESULT(pTrop) 15 21 USE assert_m, ONLY: assert 16 22 USE assert_eq_m, ONLY: assert_eq 17 23 USE dimphy, ONLY: klon, klev 18 24 USE geometry_mod, ONLY: latitude_deg, longitude_deg 19 USE vertical_layers_mod, ONLY: aps, bps, preff 25 USE strings_mod, ONLY: maxlen 26 USE yomcst_mod_h, ONLY: ROMEGA, RPI, RKAPPA, RG 27 USE vertical_layers_mod, ONLY: aps, bps, preff 20 28 USE lmdz_reprobus_wrappers, ONLY: itroprep 21 USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_REPROBUS 22 USE print_control_mod, ONLY: lunout 23 24 !------------------------------------------------------------------------------- 25 ! Arguments: 26 REAL :: dyn_tropopause(klon) !--- Pressure at tropopause 27 REAL, INTENT(IN) :: t(:,:) !--- Cells-centers temperature 28 REAL, INTENT(IN) :: ts(:) !--- Surface temperature 29 REAL, INTENT(IN) :: paprs(:,:) !--- Cells-edges pressure 30 REAL, INTENT(IN) :: pplay(:,:) !--- Cells-centers pressure 31 REAL, INTENT(IN) :: rot(:,:) !--- Cells-centers relative vorticity 32 INTEGER, INTENT(OUT), OPTIONAL :: itrop(klon) !--- Last tropospheric layer idx 33 REAL, INTENT(IN), OPTIONAL :: thet0, pvor0 34 !------------------------------------------------------------------------------- 35 ! Local variables: 36 37 REAL, PARAMETER :: DynPTrMin =8.E+3 !--- Thresholds for minimum and maximum 38 REAL, PARAMETER :: DynPTrMax =4.E+4 ! dynamical tropopause pressure (Pa). 39 CHARACTER(LEN=80) :: sub 40 INTEGER :: i, k, kb, kt, kp, ib, ie, nw 41 REAL :: al, th0, pv0 42 REAL, DIMENSION(klon,klev) :: tpot_cen, tpot_edg, pvor_cen 43 REAL, PARAMETER :: sg0=0.75 !--- Start level for PV=cte search loop 44 INTEGER, PARAMETER :: nadj=3 !--- Adjacent levs nb for thresholds detection 45 REAL, PARAMETER :: w(5)=[0.1,0.25,0.3,0.25,0.1] !--- Vertical smoothing 46 INTEGER, SAVE :: k0 47 INTEGER :: savkt 48 LOGICAL, SAVE :: first=.TRUE. 49 !$OMP THREADPRIVATE(k0,first) 50 !------------------------------------------------------------------------------- 51 sub='dyn_tropopause' 52 CALL assert(SIZE(t ,1)==klon, TRIM(sub)//" t klon") 53 CALL assert(SIZE(t ,2)==klev, TRIM(sub)//" t klev") 54 CALL assert(SIZE(ts,1)==klon, TRIM(sub)//" ts klon") 55 CALL assert(SHAPE(paprs)==[klon,klev+1],TRIM(sub)//" paprs shape") 56 CALL assert(SHAPE(pplay)==[klon,klev ],TRIM(sub)//" pplay shape") 57 CALL assert(SHAPE(rot) ==[klon,klev ],TRIM(sub)//" rot shape") 58 59 !--- DEFAULT THRESHOLDS 60 th0=380.; IF(PRESENT(thet0)) th0=thet0 !--- In kelvins 61 pv0= 2.; IF(PRESENT(pvor0)) pv0=pvor0 !--- In PVU 62 IF(first) THEN 63 DO k0=1,klev; IF(aps(k0)/preff+bps(k0)<sg0) EXIT; END DO; first=.FALSE. 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: 2*ROMEGA*SIN(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 !$OMP BARRIER 79 IF(is_master) THEN 80 81 !--- GET THE INDEX "k0" OF THE FIRST LOWER INTERFACE LAYER WITH SIGMA COORDINATE LOWER THAN "sg0" 82 !--- NOTE: "k0" DEPENDS ON VERTICAL DISCRETIZATION ONLY (VIA HYBRID COEFFS aps, bps) AND IS NOT SIMULATION-DEPENDENT 83 DO k0 = 1, klev; IF( aps(k0) / preff + bps(k0) < sg0 ) EXIT; END DO !--- START INDEX FOR BOTTOM->TOP PV SEARCH LOOP 84 85 !--- COMPUTE THE CORIOLIS PARAMETER FOR PV ALCULATION ROUTINE "potentialVorticity" 86 DO i = 1, klon 87 fac(:) = 2. * ROMEGA * SIN(latitude_deg(:) * RPI/180.) 88 END DO 89 90 !--- COMPUTE THE WEIGHTS FOR THE VERTICAL SMOOTHING ROUTINE "smooth" 91 co(:) = 0; w(:, :) = 0. 92 co(1) = 1; w(1, 1) = 1. 93 DO i = 1, ns 94 co(2:2*ns+1) = co(2:2*ns+1) + co(1:2*ns) !--- C(n+1,p+1) = C(n,p+1) + C(n,p) 95 co(2:2*ns+1) = co(2:2*ns+1) + co(1:2*ns) !--- C(n+1,p+1) = C(n,p+1) + C(n,p) AGAIN 96 w(i+1, 1:i+1) = REAL(co(i+1:2*i+1))/REAL(SUM(co(i+1:2*i+1))) 97 END DO 98 99 lFirst=.FALSE. 100 END IF 101 CALL bcast(k0) 102 CALL bcast(fac) 103 CALL bcast(w) 104 CALL bcast(lFirst) 64 105 END IF 65 106 66 !--- POTENTIAL TEMPERATURE AT CELLS CENTERS AND INTERFACES 67 DO i = 1,klon 68 tpot_cen(i,1) = t(i,1)*(preff/pplay(i,1))**RKAPPA 69 tpot_edg(i,1) = ts(i) *(preff/paprs(i,1))**RKAPPA 70 DO k=2,klev 71 al = LOG(pplay(i,k-1)/paprs(i,k))/LOG(pplay(i,k-1)/pplay(i,k)) 72 tpot_cen(i,k) = t(i,k) *(preff/pplay(i,k))**RKAPPA 73 tpot_edg(i,k) = (t(i,k-1)+al*(t(i,k)-t(i,k-1)))*(preff/paprs(i,k))**RKAPPA 74 !--- FORCE QUANTITIES TO BE GROWING 75 IF(tpot_edg(i,k)<tpot_edg(i,k-1)) tpot_edg(i,k)=tpot_edg(i,k-1)+1.E-5 76 IF(tpot_cen(i,k)<tpot_cen(i,k-1)) tpot_cen(i,k)=tpot_cen(i,k-1)+1.E-5 77 END DO 78 !--- VERTICAL SMOOTHING 79 tpot_cen(i,:)=smooth(tpot_cen(i,:),w) 80 tpot_edg(i,:)=smooth(tpot_edg(i,:),w) 81 END DO 82 83 !--- ERTEL POTENTIAL VORTICITY AT CELLS CENTERS (except in top layer) 84 DO i = 1, klon 85 DO k= 1, klev-1 86 pvor_cen(i,k)=-1.E6*RG*(rot(i,k)+2.*ROMEGA*SIN(latitude_deg(i)*RPI/180.))& 87 * (tpot_edg(i,k+1)-tpot_edg(i,k)) / (paprs(i,k+1)-paprs(i,k)) 88 END DO 89 !--- VERTICAL SMOOTHING 90 pvor_cen(i,1:klev-1)=smooth(pvor_cen(i,1:klev-1),w) 91 END DO 92 93 !--- LOCATE TROPOPAUSE: LOWEST POINT BETWEEN THETA=380K AND PV=2PVU SURFACES. 94 DO i = 1, klon 95 !--- UPPER TROPOPAUSE: |PV|=2PVU POINT STARTING FROM TOP 96 ! DO kt=klev-1,1,-1 97 ! savkt = kt 98 ! IF (kt-nadj == 0) THEN 99 ! WRITE(lunout,*)'ABORT_PHYSIC tropopause_m kt= ',kt 100 ! call abort_physic("tropopause_m", " kt = nadj", 1) 101 ! ENDIF 102 ! IF(ALL(ABS(pvor_cen(i,kt-nadj:kt))<=pv0)) THEN 103 ! EXIT 104 ! ENDIF 105 ! END DO 106 DO kt=klev-1,nadj+1,-1; savkt = kt; IF(ALL(ABS(pvor_cen(i,kt-nadj:kt))<=pv0)) EXIT; END DO 107 kt = savkt 108 !--- LOWER TROPOPAUSE: |PV|=2PVU POINT STARTING FROM BOTTOM 109 DO kb=k0,klev-1; IF(ALL(ABS(pvor_cen(i,kb:kb+nadj))> pv0)) EXIT; END DO; kb=kb-1 110 !--- ISO-THETA POINT: THETA=380K STARTING FROM TOP 111 DO kp=klev-1,1,-1; IF(ALL(ABS(tpot_cen(i,kp-nadj:kp))<=th0)) EXIT; END DO 112 !--- CHOOSE BETWEEN LOWER AND UPPER TROPOPAUSE 113 IF(2*COUNT(ABS(pvor_cen(i,kb:kt))>pv0)>kt-kb+1) kt=kb 114 !--- PV-DEFINED TROPOPAUSE 115 al = (ABS(pvor_cen(i,kt+1))-pv0)/ABS(pvor_cen(i,kt+1)-pvor_cen(i,kt)) 116 dyn_tropopause(i) = pplay(i,kt+1)*(pplay(i,kt)/pplay(i,kt+1))**al 117 !--- THETA=380K IN THE TROPICAL REGION 118 al = (tpot_cen(i,kp+1)-th0)/(tpot_cen(i,kp+1)-tpot_cen(i,kp)) 119 dyn_tropopause(i) = MAX( pplay(i,kp+1)*(pplay(i,kp)/pplay(i,kp+1))**al, & 120 dyn_tropopause(i) ) 121 !--- UNREALISTIC VALUES DETECTION 122 IF(dyn_tropopause(i)<DynPTrMin.OR.dyn_tropopause(i)>DynPTrMax) THEN 123 dyn_tropopause(i)=MIN(MAX(dyn_tropopause(i),DynPTrMax),DynPTrMin) 124 DO kt=1,klev-1; IF(pplay(i,kt+1)>dyn_tropopause(i)) EXIT; END DO; kp=kt 125 END IF 126 IF (CPPKEY_REPROBUS) THEN 127 itroprep(i)=MAX(kt,kp) 128 END IF 129 !--- LAST TROPOSPHERIC LAYER INDEX NEEDED 130 IF(PRESENT(itrop)) itrop(i)=MAX(kt,kp) 131 END DO 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 132 313 133 314 END FUNCTION dyn_tropopause 134 315 135 136 !-------------------------------------------------------------------------------137 !138 FUNCTION smooth(v,w)139 !140 !-------------------------------------------------------------------------------141 ! Arguments:142 REAL, INTENT(IN) :: v(:), w(:)143 REAL, DIMENSION(SIZE(v)) :: smooth144 !-------------------------------------------------------------------------------145 ! Local variables:146 INTEGER :: nv, nw, k, kb, ke, lb, le147 !-------------------------------------------------------------------------------148 nv=SIZE(v); nw=(SIZE(w)-1)/2149 DO k=1,nv150 kb=MAX(k-nw,1 ); lb=MAX(2+nw -k,1)151 ke=MIN(k+nw,nv); le=MIN(1+nw+nv-k,1+2*nw)152 smooth(k)=SUM(v(kb:ke)*w(lb:le))/SUM(w(lb:le))153 END DO154 155 END FUNCTION smooth156 !157 !-------------------------------------------------------------------------------158 159 316 END MODULE tropopause_m
Note: See TracChangeset
for help on using the changeset viewer.
