Changeset 14072 for NEMO/trunk/src/OCE/SBC/sbcblk_algo_ncar.F90
- Timestamp:
- 2020-12-04T08:48:38+01:00 (3 years ago)
- File:
-
- 1 edited
-
NEMO/trunk/src/OCE/SBC/sbcblk_algo_ncar.F90 (modified) (11 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/OCE/SBC/sbcblk_algo_ncar.F90
r13460 r14072 5 5 !! * bulk transfer coefficients C_D, C_E and C_H 6 6 !! * air temp. and spec. hum. adjusted from zt (2m) to zu (10m) if needed 7 !! * the effective bulk wind speed at 10m U _blk7 !! * the effective bulk wind speed at 10m Ubzu 8 8 !! => all these are used in bulk formulas in sbcblk.F90 9 9 !! … … 16 16 !!===================================================================== 17 17 !! History : 3.6 ! 2016-02 (L.Brodeau) successor of old turb_ncar of former sbcblk_core.F90 18 !! 4.2 ! 2020-12 (L. Brodeau) Introduction of various air-ice bulk parameterizations + improvements 18 19 !!---------------------------------------------------------------------- 19 20 … … 23 24 !! returns the effective bulk wind speed at 10m 24 25 !!---------------------------------------------------------------------- 25 USE oce ! ocean dynamics and tracers26 26 USE dom_oce ! ocean space and time domain 27 USE sbc_oce, ONLY: ln_cdgw 28 USE sbcwave, ONLY: cdn_wave ! wave module 27 29 USE phycst ! physical constants 28 USE sbc_oce ! Surface boundary condition: ocean fields 29 USE sbcwave, ONLY : cdn_wave ! wave module 30 #if defined key_si3 || defined key_cice 31 USE sbc_ice ! Surface boundary condition: ice fields 32 #endif 33 ! 34 USE iom ! I/O manager library 35 USE lib_mpp ! distribued memory computing library 36 USE in_out_manager ! I/O manager 37 USE prtctl ! Print control 38 USE lib_fortran ! to use key_nosignedzero 39 40 USE sbcblk_phy ! all thermodynamics functions, rho_air, q_sat, etc... !LB 30 USE sbc_phy ! Catalog of functions for physical/meteorological parameters in the marine boundary layer 41 31 42 32 IMPLICIT NONE … … 45 35 PUBLIC :: TURB_NCAR ! called by sbcblk.F90 46 36 47 INTEGER , PARAMETER :: nb_itt = 5 ! number of itterations48 37 !! * Substitutions 49 38 # include "do_loop_substitute.h90" … … 52 41 CONTAINS 53 42 54 SUBROUTINE turb_ncar( zt, zu, sst, t_zt, ssq, q_zt, U_zu, &55 & Cd, Ch, Ce, t_zu, q_zu, U_blk,&56 & Cdn, Chn, Cen)43 SUBROUTINE turb_ncar( zt, zu, sst, t_zt, ssq, q_zt, U_zu, & 44 & Cd, Ch, Ce, t_zu, q_zu, Ubzu, & 45 & nb_iter, CdN, ChN, CeN ) 57 46 !!---------------------------------------------------------------------------------- 58 47 !! *** ROUTINE turb_ncar *** … … 61 50 !! fluxes according to Large & Yeager (2004) and Large & Yeager (2008) 62 51 !! If relevant (zt /= zu), adjust temperature and humidity from height zt to zu 63 !! Returns the effective bulk wind speed at 10m to be used in the bulk formulas 64 !! 52 !! Returns the effective bulk wind speed at zu to be used in the bulk formulas 65 53 !! 66 54 !! INPUT : … … 82 70 !! * t_zu : pot. air temperature adjusted at wind height zu [K] 83 71 !! * q_zu : specific humidity of air // [kg/kg] 84 !! * U_blk : bulk wind speed at zu [m/s] 85 !! 72 !! * Ubzu : bulk wind speed at zu [m/s] 73 !! 74 !! OPTIONAL OUTPUT: 75 !! ---------------- 76 !! * CdN : neutral-stability drag coefficient 77 !! * ChN : neutral-stability sensible heat coefficient 78 !! * CeN : neutral-stability evaporation coefficient 86 79 !! 87 80 !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/) … … 99 92 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: t_zu ! pot. air temp. adjusted at zu [K] 100 93 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: q_zu ! spec. humidity adjusted at zu [kg/kg] 101 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: U_blk ! bulk wind speed at zu [m/s] 102 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Cdn, Chn, Cen ! neutral transfer coefficients 103 ! 104 INTEGER :: j_itt 94 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Ubzu ! bulk wind speed at zu [m/s] 95 ! 96 INTEGER , INTENT(in ), OPTIONAL :: nb_iter ! number of iterations 97 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: CdN 98 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: ChN 99 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: CeN 100 ! 101 INTEGER :: nbit, jit ! iterations... 105 102 LOGICAL :: l_zt_equal_zu = .FALSE. ! if q and t are given at same height as U 106 103 ! 107 REAL(wp), DIMENSION(jpi,jpj) :: Cx_n10! 10m neutral latent/sensible coefficient108 REAL(wp), DIMENSION(jpi,jpj) :: sqrt_Cd_n10 ! root square of Cd_n10104 REAL(wp), DIMENSION(jpi,jpj) :: zCdN, zCeN, zChN ! 10m neutral latent/sensible coefficient 105 REAL(wp), DIMENSION(jpi,jpj) :: zsqrt_Cd, zsqrt_CdN ! root square of Cd and Cd_neutral 109 106 REAL(wp), DIMENSION(jpi,jpj) :: zeta_u ! stability parameter at height zu 110 REAL(wp), DIMENSION(jpi,jpj) :: zpsi_h_u111 107 REAL(wp), DIMENSION(jpi,jpj) :: ztmp0, ztmp1, ztmp2 112 REAL(wp), DIMENSION(jpi,jpj) :: stab ! stability test integer 113 !!---------------------------------------------------------------------------------- 108 !!---------------------------------------------------------------------------------- 109 nbit = nb_iter0 110 IF( PRESENT(nb_iter) ) nbit = nb_iter 111 114 112 l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp ) ! testing "zu == zt" is risky with double precision 115 113 116 U _blk= MAX( 0.5_wp , U_zu ) ! relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s114 Ubzu = MAX( 0.5_wp , U_zu ) ! relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s 117 115 118 116 !! First guess of stability: 119 117 ztmp0 = virt_temp(t_zt, q_zt) - virt_temp(sst, ssq) ! air-sea difference of virtual pot. temp. at zt 120 stab = 0.5_wp + sign(0.5_wp,ztmp0) ! stab= 1 if dTv > 0 => STABLE, 0 if unstable118 ztmp1 = 0.5_wp + SIGN(0.5_wp,ztmp0) ! ztmp1 = 1 if dTv > 0 => STABLE, 0 if unstable 121 119 122 120 !! Neutral coefficients at 10m: 123 121 IF( ln_cdgw ) THEN ! wave drag case 124 122 cdn_wave(:,:) = cdn_wave(:,:) + rsmall * ( 1._wp - tmask(:,:,1) ) 125 z tmp0(:,:) = cdn_wave(:,:)123 zCdN (:,:) = cdn_wave(:,:) 126 124 ELSE 127 z tmp0 = cd_neutral_10m( U_blk)125 zCdN = cd_n10_ncar( Ubzu ) 128 126 ENDIF 129 127 130 sqrt_Cd_n10 = SQRT( ztmp0)128 zsqrt_CdN = SQRT( zCdN ) 131 129 132 130 !! Initializing transf. coeff. with their first guess neutral equivalents : 133 Cd = z tmp0134 Ce = 1.e-3_wp*( 34.6_wp * sqrt_Cd_n10)135 Ch = 1.e-3_wp*sqrt_Cd_n10*(18._wp*stab + 32.7_wp*(1._wp - stab))136 stab = sqrt_Cd_n10 ! Temporaty array !!! stab == SQRT(Cd)137 131 Cd = zCdN 132 Ce = ce_n10_ncar( zsqrt_CdN ) 133 Ch = ch_n10_ncar( zsqrt_CdN , ztmp1 ) ! ztmp1 is stability (1/0) 134 zsqrt_Cd = zsqrt_CdN 135 138 136 IF( ln_cdgw ) THEN 139 Cen= Ce140 Chn= Ch137 zCeN = Ce 138 zChN = Ch 141 139 ENDIF 142 140 143 !! First guess of temperature and humidity at height zu:141 !! Initializing values at z_u with z_t values: 144 142 t_zu = MAX( t_zt , 180._wp ) ! who knows what's given on masked-continental regions... 145 143 q_zu = MAX( q_zt , 1.e-6_wp ) ! " 146 144 145 147 146 !! ITERATION BLOCK 148 DO j _itt = 1, nb_itt147 DO jit = 1, nbit 149 148 ! 150 149 ztmp1 = t_zu - sst ! Updating air/sea differences 151 150 ztmp2 = q_zu - ssq 152 151 153 ! Updating turbulent scales : (L&Y 2004 eq. (7))154 ztmp0 = stab*U_blk ! u* (stab == SQRT(Cd))155 ztmp1 = Ch/ stab*ztmp1 ! theta* (stab == SQRT(Cd))156 ztmp2 = Ce/ stab*ztmp2 ! q* (stab == SQRT(Cd))157 158 ! Estimate the inverse of Monin-Obukov length (1/L) at height zu:152 ! Updating turbulent scales : (L&Y 2004 Eq. (7)) 153 ztmp0 = zsqrt_Cd*Ubzu ! u* 154 ztmp1 = Ch/zsqrt_Cd*ztmp1 ! theta* 155 ztmp2 = Ce/zsqrt_Cd*ztmp2 ! q* 156 157 ! Estimate the inverse of Obukov length (1/L) at height zu: 159 158 ztmp0 = One_on_L( t_zu, q_zu, ztmp0, ztmp1, ztmp2 ) 160 159 161 160 !! Stability parameters : 162 161 zeta_u = zu*ztmp0 163 zeta_u = sign( min(abs(zeta_u),10._wp), zeta_u ) 164 zpsi_h_u = psi_h( zeta_u ) 165 166 !! Shifting temperature and humidity at zu (L&Y 2004 eq. (9b-9c)) 162 zeta_u = sign( min(abs(zeta_u),10._wp), zeta_u ) 163 164 !! Shifting temperature and humidity at zu (L&Y 2004 Eq. (9b-9c)) 167 165 IF( .NOT. l_zt_equal_zu ) THEN 168 !! Array 'stab' is free for the moment so using it to store 'zeta_t'169 stab = zt*ztmp0170 stab = SIGN( MIN(ABS(stab),10._wp), stab ) ! Temporaty array stab == zeta_t !!!171 stab = LOG(zt/zu) + zpsi_h_u - psi_h(stab) ! stab just used as temp array again!172 t_zu = t_zt - ztmp1/vkarmn*stab ! ztmp1 is still theta* L&Y 2004 eq.(9b)173 q_zu = q_zt - ztmp2/vkarmn* stab ! ztmp2 is still q* L&Y 2004 eq.(9c)174 q_zu = max(0._wp, q_zu)175 END IF176 177 ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 eq. 9a)...166 ztmp0 = zt*ztmp0 ! zeta_t ! 167 ztmp0 = SIGN( MIN(ABS(ztmp0),10._wp), ztmp0 ) ! Temporaty array ztmp0 == zeta_t !!! 168 ztmp0 = LOG(zt/zu) + psi_h_ncar(zeta_u) - psi_h_ncar(ztmp0) ! ztmp0 just used as temp array again! 169 t_zu = t_zt - ztmp1/vkarmn*ztmp0 ! ztmp1 is still theta* L&Y 2004 Eq. (9b) 170 !! 171 q_zu = q_zt - ztmp2/vkarmn*ztmp0 ! ztmp2 is still q* L&Y 2004 Eq. (9c) 172 q_zu = MAX(0._wp, q_zu) 173 END IF 174 175 ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 Eq. 9a)... 178 176 ! In very rare low-wind conditions, the old way of estimating the 179 177 ! neutral wind speed at 10m leads to a negative value that causes the code 180 178 ! to crash. To prevent this a threshold of 0.25m/s is imposed. 181 ztmp2 = psi_m (zeta_u)179 ztmp2 = psi_m_ncar(zeta_u) 182 180 IF( ln_cdgw ) THEN ! surface wave case 183 stab = vkarmn / ( vkarmn / sqrt_Cd_n10 - ztmp2 ) ! (stab == SQRT(Cd))184 Cd = stab * stab185 ztmp0 = (LOG(zu/10._wp) - zpsi_h_u) / vkarmn / sqrt_Cd_n10186 ztmp2 = stab / sqrt_Cd_n10 ! (stab == SQRT(Cd))187 ztmp1 = 1._wp + Chn * ztmp0188 Ch = Chn* ztmp2 / ztmp1 ! L&Y 2004 eq. (10b)189 ztmp1 = 1._wp + Cen* ztmp0190 Ce = Cen* ztmp2 / ztmp1 ! L&Y 2004 eq. (10c)181 zsqrt_Cd = vkarmn / ( vkarmn / zsqrt_CdN - ztmp2 ) 182 Cd = zsqrt_Cd * zsqrt_Cd 183 ztmp0 = (LOG(zu/10._wp) - psi_h_ncar(zeta_u)) / vkarmn / zsqrt_CdN 184 ztmp2 = zsqrt_Cd / zsqrt_CdN 185 ztmp1 = 1._wp + zChN * ztmp0 186 Ch = zChN * ztmp2 / ztmp1 ! L&Y 2004 eq. (10b) 187 ztmp1 = 1._wp + zCeN * ztmp0 188 Ce = zCeN * ztmp2 / ztmp1 ! L&Y 2004 eq. (10c) 191 189 192 190 ELSE 193 ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 eq. 9a)... 194 ! In very rare low-wind conditions, the old way of estimating the 195 ! neutral wind speed at 10m leads to a negative value that causes the code 196 ! to crash. To prevent this a threshold of 0.25m/s is imposed. 197 ztmp0 = MAX( 0.25_wp , U_blk/(1._wp + sqrt_Cd_n10/vkarmn*(LOG(zu/10._wp) - ztmp2)) ) ! U_n10 (ztmp2 == psi_m(zeta_u)) 198 ztmp0 = cd_neutral_10m(ztmp0) ! Cd_n10 199 Cdn(:,:) = ztmp0 200 sqrt_Cd_n10 = sqrt(ztmp0) 201 202 stab = 0.5_wp + sign(0.5_wp,zeta_u) ! update stability 203 Cx_n10 = 1.e-3_wp*sqrt_Cd_n10*(18._wp*stab + 32.7_wp*(1._wp - stab)) ! L&Y 2004 eq. (6c-6d) (Cx_n10 == Ch_n10) 204 Chn(:,:) = Cx_n10 191 ztmp0 = MAX( 0.25_wp , UN10_from_CD(zu, Ubzu, Cd, ppsi=ztmp2) ) ! U_n10 (ztmp2 == psi_m_ncar(zeta_u)) 192 193 zCdN = cd_n10_ncar(ztmp0) 194 zsqrt_CdN = sqrt(zCdN) 205 195 206 196 !! Update of transfer coefficients: 207 ztmp1 = 1._wp + sqrt_Cd_n10/vkarmn*(LOG(zu/10._wp) - ztmp2) ! L&Y 2004 eq. (10a) (ztmp2 == psi_m(zeta_u)) 208 Cd = ztmp0 / ( ztmp1*ztmp1 ) 209 stab = SQRT( Cd ) ! Temporary array !!! (stab == SQRT(Cd)) 210 211 ztmp0 = (LOG(zu/10._wp) - zpsi_h_u) / vkarmn / sqrt_Cd_n10 212 ztmp2 = stab / sqrt_Cd_n10 ! (stab == SQRT(Cd)) 213 ztmp1 = 1._wp + Cx_n10*ztmp0 ! (Cx_n10 == Ch_n10) 214 Ch = Cx_n10*ztmp2 / ztmp1 ! L&Y 2004 eq. (10b) 215 216 Cx_n10 = 1.e-3_wp * (34.6_wp * sqrt_Cd_n10) ! L&Y 2004 eq. (6b) ! Cx_n10 == Ce_n10 217 Cen(:,:) = Cx_n10 218 ztmp1 = 1._wp + Cx_n10*ztmp0 219 Ce = Cx_n10*ztmp2 / ztmp1 ! L&Y 2004 eq. (10c) 197 198 !! C_D 199 ztmp1 = 1._wp + zsqrt_CdN/vkarmn*(LOG(zu/10._wp) - ztmp2) ! L&Y 2004 Eq. (10a) (ztmp2 == psi_m(zeta_u)) 200 Cd = MAX( zCdN / ( ztmp1*ztmp1 ), Cx_min ) 201 202 !! C_H and C_E 203 zsqrt_Cd = SQRT( Cd ) 204 ztmp0 = ( LOG(zu/10._wp) - psi_h_ncar(zeta_u) ) / vkarmn / zsqrt_CdN 205 ztmp2 = zsqrt_Cd / zsqrt_CdN 206 207 ztmp1 = 0.5_wp + SIGN(0.5_wp,zeta_u) ! update stability 208 zChN = 1.e-3_wp * zsqrt_CdN*(18._wp*ztmp1 + 32.7_wp*(1._wp - ztmp1)) ! L&Y 2004 eq. (6c-6d) 209 zCeN = 1.e-3_wp * (34.6_wp * zsqrt_CdN) ! L&Y 2004 eq. (6b) 210 211 Ch = MAX( zChN*ztmp2 / ( 1._wp + zChN*ztmp0 ) , Cx_min ) ! L&Y 2004 eq. (10b) 212 Ce = MAX( zCeN*ztmp2 / ( 1._wp + zCeN*ztmp0 ) , Cx_min ) ! L&Y 2004 eq. (10c) 213 220 214 ENDIF 221 215 222 END DO !DO j_itt = 1, nb_itt 216 END DO !DO jit = 1, nbit 217 218 IF(PRESENT(CdN)) CdN(:,:) = zCdN(:,:) 219 IF(PRESENT(CeN)) CeN(:,:) = zCeN(:,:) 220 IF(PRESENT(ChN)) ChN(:,:) = zChN(:,:) 223 221 224 222 END SUBROUTINE turb_ncar 225 223 226 224 227 FUNCTION cd_n eutral_10m( pw10 )228 !!---------------------------------------------------------------------------------- 225 FUNCTION cd_n10_ncar( pw10 ) 226 !!---------------------------------------------------------------------------------- 229 227 !! Estimate of the neutral drag coefficient at 10m as a function 230 228 !! of neutral wind speed at 10m 231 229 !! 232 !! Origin: Large & Yeager 2008 eq.(11a) and eq.(11b)230 !! Origin: Large & Yeager 2008, Eq. (11) 233 231 !! 234 232 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 235 233 !!---------------------------------------------------------------------------------- 236 234 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pw10 ! scalar wind speed at 10m (m/s) 237 REAL(wp), DIMENSION(jpi,jpj) :: cd_n eutral_10m235 REAL(wp), DIMENSION(jpi,jpj) :: cd_n10_ncar 238 236 ! 239 237 INTEGER :: ji, jj ! dummy loop indices … … 242 240 ! 243 241 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 244 !245 zw = pw10(ji,jj)246 zw6 = zw*zw*zw247 zw6 = zw6*zw6248 !249 ! When wind speed > 33 m/s => Cyclone conditions => special treatment250 zgt33 = 0.5_wp + SIGN( 0.5_wp, (zw - 33._wp) ) ! If pw10 < 33. => 0, else => 1251 !252 cd_neutral_10m(ji,jj) = 1.e-3_wp * ( &253 & (1._wp - zgt33)*( 2.7_wp/zw + 0.142_wp + zw/13.09_wp - 3.14807E-10_wp*zw6) & ! wind < 33 m/s254 & + zgt33 * 2.34_wp ) ! wind >= 33 m/s255 !256 cd_neutral_10m(ji,jj) = MAX(cd_neutral_10m(ji,jj), 1.E-6_wp)257 !242 ! 243 zw = pw10(ji,jj) 244 zw6 = zw*zw*zw 245 zw6 = zw6*zw6 246 ! 247 ! When wind speed > 33 m/s => Cyclone conditions => special treatment 248 zgt33 = 0.5_wp + SIGN( 0.5_wp, (zw - 33._wp) ) ! If pw10 < 33. => 0, else => 1 249 ! 250 cd_n10_ncar(ji,jj) = 1.e-3_wp * ( & 251 & (1._wp - zgt33)*( 2.7_wp/zw + 0.142_wp + zw/13.09_wp - 3.14807E-10_wp*zw6) & ! wind < 33 m/s 252 & + zgt33 * 2.34_wp ) ! wind >= 33 m/s 253 ! 254 cd_n10_ncar(ji,jj) = MAX( cd_n10_ncar(ji,jj), Cx_min ) 255 ! 258 256 END_2D 259 257 ! 260 END FUNCTION cd_neutral_10m 261 262 263 FUNCTION psi_m( pzeta ) 258 END FUNCTION cd_n10_ncar 259 260 261 FUNCTION ch_n10_ncar( psqrtcdn10 , pstab ) 262 !!---------------------------------------------------------------------------------- 263 !! Estimate of the neutral heat transfer coefficient at 10m !! 264 !! Origin: Large & Yeager 2008, Eq. (9) and (12) 265 266 !!---------------------------------------------------------------------------------- 267 REAL(wp), DIMENSION(jpi,jpj) :: ch_n10_ncar 268 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psqrtcdn10 ! sqrt( CdN10 ) 269 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pstab ! stable ABL => 1 / unstable ABL => 0 270 !!---------------------------------------------------------------------------------- 271 IF( ANY(pstab < -0.00001) .OR. ANY(pstab > 1.00001) ) THEN 272 PRINT *, 'ERROR: ch_n10_ncar@mod_blk_ncar.f90: pstab =' 273 PRINT *, pstab 274 STOP 275 END IF 276 ! 277 ch_n10_ncar = MAX( 1.e-3_wp * psqrtcdn10*( 18._wp*pstab + 32.7_wp*(1._wp - pstab) ) , Cx_min ) ! Eq. (9) & (12) Large & Yeager, 2008 278 ! 279 END FUNCTION ch_n10_ncar 280 281 FUNCTION ce_n10_ncar( psqrtcdn10 ) 282 !!---------------------------------------------------------------------------------- 283 !! Estimate of the neutral heat transfer coefficient at 10m !! 284 !! Origin: Large & Yeager 2008, Eq. (9) and (13) 285 !!---------------------------------------------------------------------------------- 286 REAL(wp), DIMENSION(jpi,jpj) :: ce_n10_ncar 287 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psqrtcdn10 ! sqrt( CdN10 ) 288 !!---------------------------------------------------------------------------------- 289 ce_n10_ncar = MAX( 1.e-3_wp * ( 34.6_wp * psqrtcdn10 ) , Cx_min ) 290 ! 291 END FUNCTION ce_n10_ncar 292 293 294 FUNCTION psi_m_ncar( pzeta ) 264 295 !!---------------------------------------------------------------------------------- 265 296 !! Universal profile stability function for momentum 266 !! !! Psis, L&Y 2004 eq. (8c), (8d), (8e)297 !! !! Psis, L&Y 2004, Eq. (8c), (8d), (8e) 267 298 !! 268 299 !! pzeta : stability paramenter, z/L where z is altitude measurement … … 271 302 !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 272 303 !!---------------------------------------------------------------------------------- 273 REAL(wp), DIMENSION(jpi,jpj) :: psi_m 304 REAL(wp), DIMENSION(jpi,jpj) :: psi_m_ncar 274 305 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta 275 306 ! 276 307 INTEGER :: ji, jj ! dummy loop indices 277 REAL(wp) :: z x2, zx,zstab ! local scalars308 REAL(wp) :: zta, zx2, zx, zpsi_unst, zpsi_stab, zstab ! local scalars 278 309 !!---------------------------------------------------------------------------------- 279 310 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 280 zx2 = SQRT( ABS( 1._wp - 16._wp*pzeta(ji,jj) ) ) 281 zx2 = MAX( zx2 , 1._wp ) 282 zx = SQRT( zx2 ) 283 zstab = 0.5_wp + SIGN( 0.5_wp , pzeta(ji,jj) ) 284 ! 285 psi_m(ji,jj) = zstab * (-5._wp*pzeta(ji,jj)) & ! Stable 286 & + (1._wp - zstab) * (2._wp*LOG((1._wp + zx)*0.5_wp) & ! Unstable 287 & + LOG((1._wp + zx2)*0.5_wp) - 2._wp*ATAN(zx) + rpi*0.5_wp) ! " 288 ! 311 zta = pzeta(ji,jj) 312 ! 313 zx2 = SQRT( ABS(1._wp - 16._wp*zta) ) ! (1 - 16z)^0.5 314 zx2 = MAX( zx2 , 1._wp ) 315 zx = SQRT(zx2) ! (1 - 16z)^0.25 316 zpsi_unst = 2._wp*LOG( (1._wp + zx )*0.5_wp ) & 317 & + LOG( (1._wp + zx2)*0.5_wp ) & 318 & - 2._wp*ATAN(zx) + rpi*0.5_wp 319 ! 320 zpsi_stab = -5._wp*zta 321 ! 322 zstab = 0.5_wp + SIGN(0.5_wp, zta) ! zta > 0 => zstab = 1 323 ! 324 psi_m_ncar(ji,jj) = zstab * zpsi_stab & ! (zta > 0) Stable 325 & + (1._wp - zstab) * zpsi_unst ! (zta < 0) Unstable 326 ! 327 ! 289 328 END_2D 290 END FUNCTION psi_m 291 292 293 FUNCTION psi_h ( pzeta )329 END FUNCTION psi_m_ncar 330 331 332 FUNCTION psi_h_ncar( pzeta ) 294 333 !!---------------------------------------------------------------------------------- 295 334 !! Universal profile stability function for temperature and humidity 296 !! !! Psis, L&Y 2004 eq. (8c), (8d), (8e)335 !! !! Psis, L&Y 2004, Eq. (8c), (8d), (8e) 297 336 !! 298 337 !! pzeta : stability paramenter, z/L where z is altitude measurement … … 301 340 !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 302 341 !!---------------------------------------------------------------------------------- 303 REAL(wp), DIMENSION(jpi,jpj) :: psi_h 342 REAL(wp), DIMENSION(jpi,jpj) :: psi_h_ncar 304 343 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta 305 344 ! 306 345 INTEGER :: ji, jj ! dummy loop indices 307 REAL(wp) :: z x2, zstab ! local scalars346 REAL(wp) :: zta, zx2, zpsi_unst, zpsi_stab, zstab ! local scalars 308 347 !!---------------------------------------------------------------------------------- 309 348 ! 310 349 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 311 zx2 = SQRT( ABS( 1._wp - 16._wp*pzeta(ji,jj) ) ) 312 zx2 = MAX( zx2 , 1._wp ) 313 zstab = 0.5_wp + SIGN( 0.5_wp , pzeta(ji,jj) ) 314 ! 315 psi_h(ji,jj) = zstab * (-5._wp*pzeta(ji,jj)) & ! Stable 316 & + (1._wp - zstab) * (2._wp*LOG( (1._wp + zx2)*0.5_wp )) ! Unstable 317 ! 350 ! 351 zta = pzeta(ji,jj) 352 ! 353 zx2 = SQRT( ABS(1._wp - 16._wp*zta) ) ! (1 -16z)^0.5 354 zx2 = MAX( zx2 , 1._wp ) 355 zpsi_unst = 2._wp*LOG( 0.5_wp*(1._wp + zx2) ) 356 ! 357 zpsi_stab = -5._wp*zta 358 ! 359 zstab = 0.5_wp + SIGN(0.5_wp, zta) ! zta > 0 => zstab = 1 360 ! 361 psi_h_ncar(ji,jj) = zstab * zpsi_stab & ! (zta > 0) Stable 362 & + (1._wp - zstab) * zpsi_unst ! (zta < 0) Unstable 363 ! 318 364 END_2D 319 END FUNCTION psi_h 365 END FUNCTION psi_h_ncar 320 366 321 367 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
