[13820] | 1 | MODULE sbcblk_algo_ice_an05 |
---|
| 2 | !!====================================================================== |
---|
| 3 | !! *** MODULE sbcblk_algo_ice_an05 *** |
---|
| 4 | !! Computes turbulent components of surface fluxes over sea-ice |
---|
| 5 | !! |
---|
| 6 | !! Andreas, E.L., Jordan, R.E. & Makshtas, A.P. Parameterizing turbulent exchange over sea ice: the ice station weddell results. |
---|
| 7 | !! Boundary-Layer Meteorology 114, 439–460 (2005). https://doi.org/10.1007/s10546-004-1414-7 |
---|
| 8 | !! |
---|
| 9 | !! * bulk transfer coefficients C_D, C_E and C_H |
---|
| 10 | !! * air temp. and spec. hum. adjusted from zt (usually 2m) to zu (usually 10m) if needed |
---|
| 11 | !! * the "effective" bulk wind speed at zu: Ub (including gustiness contribution in unstable conditions) |
---|
| 12 | !! => all these are used in bulk formulas in sbcblk.F90 |
---|
| 13 | !! |
---|
| 14 | !! Routine turb_ice_an05 maintained and developed in AeroBulk |
---|
| 15 | !! (https://github.com/brodeau/aerobulk/) |
---|
| 16 | !! |
---|
| 17 | !! Author: Laurent Brodeau, Summer 2020 |
---|
| 18 | !! |
---|
| 19 | !!---------------------------------------------------------------------- |
---|
| 20 | USE par_kind, ONLY: wp |
---|
[14021] | 21 | USE par_oce, ONLY: jpi, jpj, Nis0, Nie0, Njs0, Nje0, nn_hls, ntsi, ntsj, ntei, ntej |
---|
[13830] | 22 | USE lib_mpp, ONLY: ctl_stop ! distribued memory computing library |
---|
[13820] | 23 | USE phycst ! physical constants |
---|
| 24 | USE sbc_phy ! Catalog of functions for physical/meteorological parameters in the marine boundary layer |
---|
| 25 | |
---|
| 26 | IMPLICIT NONE |
---|
| 27 | PRIVATE |
---|
| 28 | |
---|
| 29 | PUBLIC :: turb_ice_an05 |
---|
| 30 | |
---|
| 31 | INTEGER , PARAMETER :: nbit = 8 ! number of itterations |
---|
| 32 | |
---|
[13830] | 33 | !! * Substitutions |
---|
| 34 | # include "do_loop_substitute.h90" |
---|
[13820] | 35 | !!---------------------------------------------------------------------- |
---|
| 36 | CONTAINS |
---|
| 37 | |
---|
[13830] | 38 | SUBROUTINE turb_ice_an05( zt, zu, Ts_i, t_zt, qs_i, q_zt, U_zu, & |
---|
| 39 | & Cd_i, Ch_i, Ce_i, t_zu_i, q_zu_i, & |
---|
[13820] | 40 | & CdN, ChN, CeN, xz0, xu_star, xL, xUN10 ) |
---|
| 41 | !!---------------------------------------------------------------------- |
---|
| 42 | !! *** ROUTINE turb_ice_an05 *** |
---|
| 43 | !! |
---|
| 44 | !! ** Purpose : Computes turbulent transfert coefficients of surface |
---|
| 45 | !! fluxes according to: |
---|
| 46 | !! Andreas, E.L., Jordan, R.E. & Makshtas, A.P. Parameterizing turbulent exchange over sea ice: the ice station weddell results. |
---|
| 47 | !! Boundary-Layer Meteorology 114, 439–460 (2005). https://doi.org/10.1007/s10546-004-1414-7 |
---|
| 48 | !! |
---|
| 49 | !! If relevant (zt /= zu), adjust temperature and humidity from height zt to zu |
---|
| 50 | !! Returns the effective bulk wind speed at zu to be used in the bulk formulas |
---|
| 51 | !! |
---|
| 52 | !! INPUT : |
---|
| 53 | !! ------- |
---|
| 54 | !! * zt : height for temperature and spec. hum. of air [m] |
---|
| 55 | !! * zu : height for wind speed (usually 10m) [m] |
---|
| 56 | !! * Ts_i : surface temperature of sea-ice [K] |
---|
| 57 | !! * t_zt : potential air temperature at zt [K] |
---|
| 58 | !! * qs_i : saturation specific humidity at temp. Ts_i over ice [kg/kg] |
---|
| 59 | !! * q_zt : specific humidity of air at zt [kg/kg] |
---|
| 60 | !! * U_zu : scalar wind speed at zu [m/s] |
---|
| 61 | !! |
---|
| 62 | !! OUTPUT : |
---|
| 63 | !! -------- |
---|
| 64 | !! * Cd_i : drag coefficient over sea-ice |
---|
| 65 | !! * Ch_i : sensible heat coefficient over sea-ice |
---|
| 66 | !! * Ce_i : sublimation coefficient over sea-ice |
---|
| 67 | !! * t_zu_i : pot. air temp. adjusted at zu over sea-ice [K] |
---|
| 68 | !! * q_zu_i : spec. hum. of air adjusted at zu over sea-ice [kg/kg] |
---|
| 69 | !! |
---|
| 70 | !! OPTIONAL OUTPUT: |
---|
| 71 | !! ---------------- |
---|
| 72 | !! * CdN : neutral-stability drag coefficient |
---|
| 73 | !! * ChN : neutral-stability sensible heat coefficient |
---|
| 74 | !! * CeN : neutral-stability evaporation coefficient |
---|
| 75 | !! * xz0 : return the aerodynamic roughness length (integration constant for wind stress) [m] |
---|
| 76 | !! * xu_star : return u* the friction velocity [m/s] |
---|
| 77 | !! * xL : return the Obukhov length [m] |
---|
| 78 | !! * xUN10 : neutral wind speed at 10m [m/s] |
---|
| 79 | !! |
---|
| 80 | !! ** Author: L. Brodeau, January 2020 / AeroBulk (https://github.com/brodeau/aerobulk/) |
---|
| 81 | !!---------------------------------------------------------------------------------- |
---|
| 82 | REAL(wp), INTENT(in ) :: zt ! height for t_zt and q_zt [m] |
---|
| 83 | REAL(wp), INTENT(in ) :: zu ! height for U_zu [m] |
---|
| 84 | REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: Ts_i ! ice surface temperature [Kelvin] |
---|
| 85 | REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: t_zt ! potential air temperature [Kelvin] |
---|
| 86 | REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: qs_i ! sat. spec. hum. at ice/air interface [kg/kg] |
---|
| 87 | REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: q_zt ! spec. air humidity at zt [kg/kg] |
---|
| 88 | REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: U_zu ! relative wind module at zu [m/s] |
---|
| 89 | REAL(wp), INTENT(out), DIMENSION(jpi,jpj) :: Cd_i ! drag coefficient over sea-ice |
---|
| 90 | REAL(wp), INTENT(out), DIMENSION(jpi,jpj) :: Ch_i ! transfert coefficient for heat over ice |
---|
| 91 | REAL(wp), INTENT(out), DIMENSION(jpi,jpj) :: Ce_i ! transfert coefficient for sublimation over ice |
---|
| 92 | REAL(wp), INTENT(out), DIMENSION(jpi,jpj) :: t_zu_i ! pot. air temp. adjusted at zu [K] |
---|
| 93 | REAL(wp), INTENT(out), DIMENSION(jpi,jpj) :: q_zu_i ! spec. humidity adjusted at zu [kg/kg] |
---|
| 94 | !!---------------------------------------------------------------------------------- |
---|
| 95 | REAL(wp), INTENT(out), DIMENSION(jpi,jpj), OPTIONAL :: CdN |
---|
| 96 | REAL(wp), INTENT(out), DIMENSION(jpi,jpj), OPTIONAL :: ChN |
---|
| 97 | REAL(wp), INTENT(out), DIMENSION(jpi,jpj), OPTIONAL :: CeN |
---|
| 98 | REAL(wp), INTENT(out), DIMENSION(jpi,jpj), OPTIONAL :: xz0 ! Aerodynamic roughness length [m] |
---|
| 99 | REAL(wp), INTENT(out), DIMENSION(jpi,jpj), OPTIONAL :: xu_star ! u*, friction velocity |
---|
| 100 | REAL(wp), INTENT(out), DIMENSION(jpi,jpj), OPTIONAL :: xL ! zeta (zu/L) |
---|
| 101 | REAL(wp), INTENT(out), DIMENSION(jpi,jpj), OPTIONAL :: xUN10 ! Neutral wind at zu |
---|
| 102 | !!---------------------------------------------------------------------------------- |
---|
| 103 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: Ubzu |
---|
[13830] | 104 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ztmp0, ztmp1, ztmp2 ! temporary stuff |
---|
| 105 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: z0, dt_zu, dq_zu |
---|
[13820] | 106 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: u_star, t_star, q_star |
---|
[13830] | 107 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: znu_a !: Nu_air = kinematic viscosity of air |
---|
| 108 | REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zeta_u, zeta_t ! stability parameter at height zu |
---|
[13820] | 109 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: z0tq |
---|
| 110 | !! |
---|
| 111 | INTEGER :: jit |
---|
| 112 | LOGICAL :: l_zt_equal_zu = .FALSE. ! if q and t are given at same height as U |
---|
| 113 | !! |
---|
| 114 | LOGICAL :: lreturn_cdn=.FALSE., lreturn_chn=.FALSE., lreturn_cen=.FALSE. |
---|
| 115 | LOGICAL :: lreturn_z0=.FALSE., lreturn_ustar=.FALSE., lreturn_L=.FALSE., lreturn_UN10=.FALSE. |
---|
| 116 | !! |
---|
| 117 | CHARACTER(len=40), PARAMETER :: crtnm = 'turb_ice_an05@sbcblk_algo_ice_an05.f90' |
---|
| 118 | !!---------------------------------------------------------------------------------- |
---|
[13830] | 119 | ALLOCATE ( Ubzu(jpi,jpj), u_star(jpi,jpj), t_star(jpi,jpj), q_star(jpi,jpj), & |
---|
| 120 | & zeta_u(jpi,jpj), dt_zu(jpi,jpj), dq_zu(jpi,jpj), & |
---|
| 121 | & znu_a(jpi,jpj), ztmp1(jpi,jpj), ztmp2(jpi,jpj), & |
---|
| 122 | & z0(jpi,jpj), z0tq(jpi,jpj,2), ztmp0(jpi,jpj) ) |
---|
[13820] | 123 | |
---|
[13830] | 124 | lreturn_cdn = PRESENT(CdN) |
---|
| 125 | lreturn_chn = PRESENT(ChN) |
---|
| 126 | lreturn_cen = PRESENT(CeN) |
---|
| 127 | lreturn_z0 = PRESENT(xz0) |
---|
| 128 | lreturn_ustar = PRESENT(xu_star) |
---|
| 129 | lreturn_L = PRESENT(xL) |
---|
| 130 | lreturn_UN10 = PRESENT(xUN10) |
---|
[13820] | 131 | |
---|
| 132 | l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp ) |
---|
| 133 | IF( .NOT. l_zt_equal_zu ) ALLOCATE( zeta_t(jpi,jpj) ) |
---|
| 134 | |
---|
| 135 | !! Scalar wind speed cannot be below 0.2 m/s |
---|
| 136 | Ubzu = MAX( U_zu, wspd_thrshld_ice ) |
---|
| 137 | |
---|
| 138 | !! First guess of temperature and humidity at height zu: |
---|
| 139 | t_zu_i = MAX( t_zt , 100._wp ) ! who knows what's given on masked-continental regions... |
---|
| 140 | q_zu_i = MAX( q_zt , 0.1e-6_wp ) ! " |
---|
| 141 | |
---|
| 142 | !! Air-Ice differences (and we don't want it to be 0!) |
---|
| 143 | dt_zu = t_zu_i - Ts_i ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) |
---|
| 144 | dq_zu = q_zu_i - qs_i ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) |
---|
| 145 | |
---|
| 146 | znu_a = visc_air(t_zu_i) ! Air viscosity (m^2/s) at zt given from temperature in (K) |
---|
| 147 | |
---|
| 148 | !! Very crude first guesses of z0: |
---|
| 149 | z0 = 8.0E-4_wp |
---|
| 150 | |
---|
| 151 | !! Crude first guess of turbulent scales |
---|
| 152 | u_star = 0.035_wp*Ubzu*LOG( 10._wp/z0 )/LOG( zu/z0 ) |
---|
| 153 | z0 = rough_leng_m( u_star , znu_a ) |
---|
| 154 | |
---|
| 155 | DO jit = 1, 2 |
---|
| 156 | u_star = MAX ( Ubzu*vkarmn/(LOG(zu) - LOG(z0)) , 1.E-9 ) |
---|
| 157 | z0 = rough_leng_m( u_star , znu_a ) |
---|
| 158 | END DO |
---|
| 159 | |
---|
| 160 | z0tq = rough_leng_tq( z0, u_star , znu_a ) |
---|
| 161 | t_star = dt_zu*vkarmn/(LOG(zu/z0tq(:,:,1))) |
---|
| 162 | q_star = dq_zu*vkarmn/(LOG(zu/z0tq(:,:,2))) |
---|
[13830] | 163 | |
---|
| 164 | |
---|
[13820] | 165 | !! ITERATION BLOCK |
---|
| 166 | DO jit = 1, nbit |
---|
| 167 | |
---|
| 168 | !!Inverse of Obukov length (1/L) : |
---|
| 169 | ztmp0 = One_on_L(t_zu_i, q_zu_i, u_star, t_star, q_star) ! 1/L == 1/[Obukhov length] |
---|
| 170 | ztmp0 = SIGN( MIN(ABS(ztmp0),200._wp), ztmp0 ) ! (prevents FPE from stupid values from masked region later on...) |
---|
| 171 | |
---|
| 172 | !! Stability parameters "zeta" : |
---|
| 173 | zeta_u = zu*ztmp0 |
---|
| 174 | zeta_u = SIGN( MIN(ABS(zeta_u),50.0_wp), zeta_u ) |
---|
| 175 | IF( .NOT. l_zt_equal_zu ) THEN |
---|
| 176 | zeta_t = zt*ztmp0 |
---|
| 177 | zeta_t = SIGN( MIN(ABS(zeta_t),50.0_wp), zeta_t ) |
---|
| 178 | END IF |
---|
| 179 | |
---|
| 180 | !! Roughness lengthes z0, z0t, & z0q : |
---|
| 181 | z0 = rough_leng_m ( u_star , znu_a ) |
---|
| 182 | z0tq = rough_leng_tq( z0, u_star , znu_a ) |
---|
| 183 | |
---|
| 184 | !! Turbulent scales at zu : |
---|
| 185 | ztmp0 = psi_h_ice(zeta_u) |
---|
| 186 | t_star = dt_zu*vkarmn/(LOG(zu) - LOG(z0tq(:,:,1)) - ztmp0) |
---|
| 187 | q_star = dq_zu*vkarmn/(LOG(zu) - LOG(z0tq(:,:,2)) - ztmp0) |
---|
| 188 | u_star = MAX( Ubzu*vkarmn/(LOG(zu) - LOG(z0(:,:)) - psi_m_ice(zeta_u)) , 1.E-9 ) |
---|
| 189 | |
---|
| 190 | IF( .NOT. l_zt_equal_zu ) THEN |
---|
| 191 | !! Re-updating temperature and humidity at zu if zt /= zu : |
---|
| 192 | ztmp1 = LOG(zt/zu) + ztmp0 - psi_h_ice(zeta_t) |
---|
| 193 | t_zu_i = t_zt - t_star/vkarmn*ztmp1 |
---|
| 194 | q_zu_i = q_zt - q_star/vkarmn*ztmp1 |
---|
| 195 | dt_zu = t_zu_i - Ts_i ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) |
---|
| 196 | dq_zu = q_zu_i - qs_i ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) |
---|
| 197 | END IF |
---|
| 198 | |
---|
| 199 | END DO !DO jit = 1, nbit |
---|
| 200 | |
---|
| 201 | ! compute transfer coefficients at zu : |
---|
| 202 | ztmp0 = u_star/Ubzu |
---|
| 203 | Cd_i = ztmp0*ztmp0 |
---|
| 204 | Ch_i = ztmp0*t_star/dt_zu |
---|
| 205 | Ce_i = ztmp0*q_star/dq_zu |
---|
| 206 | |
---|
| 207 | IF( lreturn_cdn .OR. lreturn_chn .OR. lreturn_cen ) ztmp0 = 1._wp/LOG( zu/z0(:,:) ) |
---|
| 208 | IF( lreturn_cdn ) CdN = vkarmn2*ztmp0*ztmp0 |
---|
| 209 | IF( lreturn_chn ) ChN = vkarmn2*ztmp0/LOG(zu/z0tq(:,:,1)) |
---|
| 210 | IF( lreturn_cen ) CeN = vkarmn2*ztmp0/LOG(zu/z0tq(:,:,2)) |
---|
| 211 | |
---|
| 212 | IF( lreturn_z0 ) xz0 = z0 |
---|
| 213 | IF( lreturn_ustar ) xu_star = u_star |
---|
| 214 | IF( lreturn_L ) xL = 1./One_on_L(t_zu_i, q_zu_i, u_star, t_star, q_star) |
---|
| 215 | IF( lreturn_UN10 ) xUN10 = u_star/vkarmn*LOG(10./z0) |
---|
| 216 | |
---|
[13830] | 217 | DEALLOCATE ( Ubzu, u_star, t_star, q_star, zeta_u, dt_zu, dq_zu, z0, z0tq, znu_a, ztmp0, ztmp1, ztmp2 ) |
---|
[13820] | 218 | IF( .NOT. l_zt_equal_zu ) DEALLOCATE( zeta_t ) |
---|
| 219 | |
---|
| 220 | END SUBROUTINE turb_ice_an05 |
---|
| 221 | |
---|
| 222 | |
---|
| 223 | |
---|
| 224 | FUNCTION rough_leng_m( pus , pnua ) |
---|
| 225 | !!---------------------------------------------------------------------------------- |
---|
| 226 | !! Computes the roughness length of sea-ice according to Andreas et al. 2005, (eq. 19) |
---|
| 227 | !! |
---|
| 228 | !! Author: L. Brodeau, January 2020 / AeroBulk (https://github.com/brodeau/aerobulk/) |
---|
| 229 | !!---------------------------------------------------------------------------------- |
---|
| 230 | REAL(wp), DIMENSION(jpi,jpj) :: rough_leng_m ! roughness length over sea-ice [m] |
---|
| 231 | REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pus ! u* = friction velocity [m/s] |
---|
| 232 | REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pnua ! kinematic viscosity of air [m^2/s] |
---|
| 233 | !! |
---|
| 234 | INTEGER :: ji, jj ! dummy loop indices |
---|
| 235 | REAL(wp) :: zus, zz |
---|
| 236 | !!---------------------------------------------------------------------------------- |
---|
| 237 | DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) |
---|
| 238 | zus = MAX( pus(ji,jj) , 1.E-9_wp ) |
---|
| 239 | |
---|
| 240 | zz = (zus - 0.18_wp) / 0.1_wp |
---|
| 241 | |
---|
| 242 | rough_leng_m(ji,jj) = 0.135*pnua(ji,jj)/zus + 0.035*zus*zus/grav*( 5.*EXP(-zz*zz) + 1._wp ) ! Eq.(19) Andreas et al., 2005 |
---|
| 243 | END_2D |
---|
| 244 | !! |
---|
| 245 | END FUNCTION rough_leng_m |
---|
| 246 | |
---|
| 247 | FUNCTION rough_leng_tq( pz0, pus , pnua ) |
---|
| 248 | !!---------------------------------------------------------------------------------- |
---|
| 249 | !! Computes the roughness length of sea-ice according to Andreas et al. 2005, (eq. 22) |
---|
| 250 | !! => which still relies on Andreas 1987 ! |
---|
| 251 | !! |
---|
| 252 | !! Author: L. Brodeau, January 2020 / AeroBulk (https://github.com/brodeau/aerobulk/) |
---|
| 253 | !!---------------------------------------------------------------------------------- |
---|
| 254 | REAL(wp), DIMENSION(jpi,jpj,2) :: rough_leng_tq ! temp.,hum. roughness lengthes over sea-ice [m] |
---|
| 255 | REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pz0 ! roughness length [m] |
---|
| 256 | REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pus ! u* = friction velocity [m/s] |
---|
| 257 | REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pnua ! kinematic viscosity of air [m^2/s] |
---|
| 258 | !! |
---|
| 259 | INTEGER :: ji, jj ! dummy loop indices |
---|
| 260 | REAL(wp) :: zz0, zus, zre, zsmoot, ztrans, zrough |
---|
| 261 | REAL(wp) :: zb0, zb1, zb2, zlog, zlog2, zlog_z0s_on_z0 |
---|
| 262 | !!---------------------------------------------------------------------------------- |
---|
| 263 | DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) |
---|
| 264 | zz0 = pz0(ji,jj) |
---|
| 265 | zus = MAX( pus(ji,jj) , 1.E-9_wp ) |
---|
| 266 | zre = MAX( zus*zz0/pnua(ji,jj) , 0._wp ) ! Roughness Reynolds number |
---|
| 267 | |
---|
| 268 | !! *** TABLE 1 of Andreas et al. 2005 *** |
---|
[15122] | 269 | zsmoot = 0._wp ; ztrans = 0._wp ; zrough = 0._wp |
---|
| 270 | IF ( zre <= 0.135_wp ) THEN ! Smooth flow condition (R* <= 0.135): |
---|
| 271 | zsmoot = 1._wp |
---|
| 272 | ELSEIF( zre < 2.5_wp ) THEN ! Transition (0.135 < R* < 2.5) |
---|
| 273 | ztrans = 1._wp |
---|
| 274 | ELSE ! Rough ( R* > 2.5) |
---|
| 275 | zrough = 1._wp |
---|
| 276 | ENDIF |
---|
| 277 | |
---|
[13820] | 278 | zlog = LOG(zre) |
---|
| 279 | zlog2 = zlog*zlog |
---|
| 280 | |
---|
| 281 | !! z0t: |
---|
| 282 | zb0 = zsmoot*1.25_wp + ztrans*0.149_wp + zrough*0.317_wp |
---|
| 283 | zb1 = - ztrans*0.550_wp - zrough*0.565_wp |
---|
| 284 | zb2 = - zrough*0.183_wp |
---|
| 285 | zlog_z0s_on_z0 = zb0 + zb1*zlog + zb2*zlog2 |
---|
| 286 | rough_leng_tq(ji,jj,1) = zz0 * EXP( zlog_z0s_on_z0 ) |
---|
| 287 | |
---|
| 288 | !! z0q: |
---|
| 289 | zb0 = zsmoot*1.61_wp + ztrans*0.351_wp + zrough*0.396_wp |
---|
| 290 | zb1 = - ztrans*0.628_wp - zrough*0.512_wp |
---|
| 291 | zb2 = - zrough*0.180_wp |
---|
| 292 | zlog = LOG(zre) |
---|
| 293 | zlog_z0s_on_z0 = zb0 + zb1*zlog + zb2*zlog2 |
---|
| 294 | rough_leng_tq(ji,jj,2) = zz0 * EXP( zlog_z0s_on_z0 ) |
---|
| 295 | |
---|
| 296 | END_2D |
---|
| 297 | !! |
---|
| 298 | END FUNCTION rough_leng_tq |
---|
| 299 | |
---|
| 300 | |
---|
| 301 | |
---|
| 302 | FUNCTION psi_m_ice( pzeta ) |
---|
| 303 | !!---------------------------------------------------------------------------------- |
---|
| 304 | !! ** Purpose: compute the universal profile stability function for momentum |
---|
| 305 | !! |
---|
| 306 | !! |
---|
| 307 | !! Andreas et al 2005 == Jordan et al. 1999 |
---|
| 308 | !! |
---|
| 309 | !! Psi: |
---|
| 310 | !! Unstable => Paulson 1970 |
---|
| 311 | !! Stable => Holtslag & De Bruin 1988 |
---|
| 312 | !! |
---|
| 313 | !! pzeta : stability paramenter, z/L where z is altitude |
---|
| 314 | !! measurement and L is M-O length |
---|
| 315 | !! |
---|
| 316 | !! ** Author: L. Brodeau, 2020 / AeroBulk (https://github.com/brodeau/aerobulk/) |
---|
| 317 | !!---------------------------------------------------------------------------------- |
---|
| 318 | REAL(wp), DIMENSION(jpi,jpj) :: psi_m_ice |
---|
| 319 | REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta |
---|
| 320 | ! |
---|
| 321 | INTEGER :: ji, jj ! dummy loop indices |
---|
| 322 | REAL(wp) :: zta, zx, zpsi_u, zpsi_s, zstab |
---|
| 323 | !!---------------------------------------------------------------------------------- |
---|
| 324 | DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) ! |
---|
| 325 | zta = pzeta(ji,jj) |
---|
| 326 | ! |
---|
| 327 | ! Unstable stratification: |
---|
| 328 | zx = ABS(1._wp - 16._wp*zta)**.25 ! (16 here, not 15!) |
---|
| 329 | |
---|
| 330 | zpsi_u = LOG( (1._wp + zx*zx)/2. ) & ! Eq.(30) Jordan et al. 1999 |
---|
| 331 | & + 2.*LOG( (1._wp + zx )/2. ) & |
---|
| 332 | & - 2.*ATAN( zx ) + 0.5*rpi |
---|
| 333 | |
---|
| 334 | ! Stable stratification: |
---|
| 335 | zpsi_s = - ( 0.7_wp*zta + 0.75_wp*(zta - 14.3_wp)*EXP( -0.35*zta) + 10.7_wp ) ! Eq.(33) Jordan et al. 1999 |
---|
| 336 | |
---|
| 337 | !! Combine: |
---|
| 338 | zstab = 0.5 + SIGN(0.5_wp, zta) |
---|
| 339 | psi_m_ice(ji,jj) = (1._wp - zstab) * zpsi_u & ! Unstable (zta < 0) |
---|
| 340 | & + zstab * zpsi_s ! Stable (zta > 0) |
---|
| 341 | ! |
---|
| 342 | END_2D |
---|
| 343 | END FUNCTION psi_m_ice |
---|
| 344 | |
---|
| 345 | |
---|
| 346 | FUNCTION psi_h_ice( pzeta ) |
---|
| 347 | !!---------------------------------------------------------------------------------- |
---|
| 348 | !! ** Purpose: compute the universal profile stability function for |
---|
| 349 | !! temperature and humidity |
---|
| 350 | !! |
---|
| 351 | !! |
---|
| 352 | !! Andreas et al 2005 == Jordan et al. 1999 |
---|
| 353 | !! |
---|
| 354 | !! Psi: |
---|
| 355 | !! Unstable => Paulson 1970 |
---|
| 356 | !! Stable => Holtslag & De Bruin 1988 |
---|
| 357 | !! |
---|
| 358 | !! pzeta : stability paramenter, z/L where z is altitude |
---|
| 359 | !! measurement and L is M-O length |
---|
| 360 | !! |
---|
| 361 | !! ** Author: L. Brodeau, 2020 / AeroBulk (https://github.com/brodeau/aerobulk/) |
---|
| 362 | !!---------------------------------------------------------------------------------- |
---|
| 363 | REAL(wp), DIMENSION(jpi,jpj) :: psi_h_ice |
---|
| 364 | REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta |
---|
| 365 | ! |
---|
| 366 | INTEGER :: ji, jj ! dummy loop indices |
---|
| 367 | REAL(wp) :: zta, zx, zpsi_u, zpsi_s, zstab |
---|
| 368 | !!---------------------------------------------------------------------------------- |
---|
| 369 | DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) ! |
---|
| 370 | zta = pzeta(ji,jj) |
---|
| 371 | ! |
---|
| 372 | ! Unstable stratification: |
---|
| 373 | zx = ABS(1._wp - 16._wp*zta)**.25 ! (16 here, not 15!) |
---|
| 374 | |
---|
| 375 | zpsi_u = 2.*LOG( (1._wp + zx*zx)/2. ) ! Eq.(31) Jordan et al. 1999 |
---|
| 376 | |
---|
| 377 | ! Stable stratification (identical to Psi_m!): |
---|
| 378 | zpsi_s = - ( 0.7_wp*zta + 0.75_wp*(zta - 14.3_wp)*EXP( -0.35*zta) + 10.7_wp ) ! Eq.(33) Jordan et al. 1999 |
---|
| 379 | |
---|
| 380 | !! Combine: |
---|
| 381 | zstab = 0.5 + SIGN(0.5_wp, zta) |
---|
| 382 | psi_h_ice(ji,jj) = (1._wp - zstab) * zpsi_u & ! Unstable (zta < 0) |
---|
| 383 | & + zstab * zpsi_s ! Stable (zta > 0) |
---|
| 384 | ! |
---|
| 385 | END_2D |
---|
| 386 | END FUNCTION psi_h_ice |
---|
| 387 | |
---|
| 388 | !!====================================================================== |
---|
| 389 | END MODULE sbcblk_algo_ice_an05 |
---|