Changeset 13655 for NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbcblk_algo_ecmwf.F90
- Timestamp:
- 2020-10-21T16:15:13+02:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbcblk_algo_ecmwf.F90
r13460 r13655 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 !! … … 24 24 !! returns the effective bulk wind speed at 10m 25 25 !!---------------------------------------------------------------------- 26 USE oce ! ocean dynamics and tracers27 26 USE dom_oce ! ocean space and time domain 28 27 USE phycst ! physical constants 29 USE iom ! I/O manager library 30 USE lib_mpp ! distribued memory computing library 31 USE in_out_manager ! I/O manager 32 USE prtctl ! Print control 33 USE sbcwave, ONLY : cdn_wave ! wave module 34 #if defined key_si3 || defined key_cice 35 USE sbc_ice ! Surface boundary condition: ice fields 36 #endif 37 USE lib_fortran ! to use key_nosignedzero 38 39 USE sbc_oce ! Surface boundary condition: ocean fields 40 USE sbcblk_phy ! all thermodynamics functions, rho_air, q_sat, etc... !LB 28 USE lib_mpp, ONLY: ctl_stop ! distribued memory computing library 29 USE in_out_manager, ONLY: nit000 ! I/O manager 30 USE sbc_phy ! all thermodynamics functions, rho_air, q_sat, etc... !LB 41 31 USE sbcblk_skin_ecmwf ! cool-skin/warm layer scheme !LB 42 32 … … 45 35 46 36 PUBLIC :: SBCBLK_ALGO_ECMWF_INIT, TURB_ECMWF 47 !! * Substitutions48 # include "do_loop_substitute.h90"49 37 50 38 !! ECMWF own values for given constants, taken form IFS documentation... 51 REAL(wp), PARAMETER :: charn0 = 0.018! Charnock constant (pretty high value here !!!39 REAL(wp), PARAMETER, PUBLIC :: charn0_ecmwf = 0.018_wp ! Charnock constant (pretty high value here !!! 52 40 ! ! => Usually 0.011 for moderate winds) 53 41 REAL(wp), PARAMETER :: zi0 = 1000. ! scale height of the atmospheric boundary layer...1 … … 57 45 REAL(wp), PARAMETER :: alpha_Q = 0.62 ! 58 46 59 INTEGER , PARAMETER :: nb_itt = 10 ! number of itterations 47 !! * Substitutions 48 # include "do_loop_substitute.h90" 60 49 61 50 !!---------------------------------------------------------------------- … … 94 83 95 84 SUBROUTINE turb_ecmwf( kt, zt, zu, T_s, t_zt, q_s, q_zt, U_zu, l_use_cs, l_use_wl, & 96 & Cd, Ch, Ce, t_zu, q_zu, U _blk,&97 & Cdn, Chn, Cen, &85 & Cd, Ch, Ce, t_zu, q_zu, Ubzu, & 86 & nb_iter, Cdn, Chn, Cen, & ! optional output 98 87 & Qsw, rad_lw, slp, pdT_cs, & ! optionals for cool-skin (and warm-layer) 99 88 & pdT_wl, pHz_wl ) ! optionals for warm-layer only … … 151 140 !! * t_zu : pot. air temperature adjusted at wind height zu [K] 152 141 !! * q_zu : specific humidity of air // [kg/kg] 153 !! * U _blk: bulk wind speed at zu [m/s]142 !! * Ubzu : bulk wind speed at zu [m/s] 154 143 !! 155 144 !! … … 171 160 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: t_zu ! pot. air temp. adjusted at zu [K] 172 161 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: q_zu ! spec. humidity adjusted at zu [kg/kg] 173 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: U_blk ! bulk wind speed at zu [m/s] 174 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Cdn, Chn, Cen ! neutral transfer coefficients 175 ! 162 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Ubzu ! bulk wind speed at zu [m/s] 163 ! 164 INTEGER , INTENT(in ), OPTIONAL :: nb_iter ! number of iterations 165 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: CdN 166 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: ChN 167 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: CeN 176 168 REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: Qsw ! [W/m^2] 177 169 REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: rad_lw ! [W/m^2] … … 181 173 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: pHz_wl ! [m] 182 174 ! 183 INTEGER :: j_itt175 INTEGER :: nbit, jit 184 176 LOGICAL :: l_zt_equal_zu = .FALSE. ! if q and t are given at same height as U 185 177 ! … … 197 189 !!---------------------------------------------------------------------------------- 198 190 IF( kt == nit000 ) CALL SBCBLK_ALGO_ECMWF_INIT(l_use_cs, l_use_wl) 191 192 nbit = nb_iter0 193 IF( PRESENT(nb_iter) ) nbit = nb_iter 199 194 200 195 l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp ) ! testing "zu == zt" is risky with double precision … … 227 222 znu_a = visc_air(t_zu) ! Air viscosity (m^2/s) at zt given from temperature in (K) 228 223 229 U _blk= SQRT(U_zu*U_zu + 0.5_wp*0.5_wp) ! initial guess for wind gustiness contribution224 Ubzu = SQRT(U_zu*U_zu + 0.5_wp*0.5_wp) ! initial guess for wind gustiness contribution 230 225 231 226 ztmp0 = LOG( zu*10000._wp) ! optimization: 10000. == 1/z0 (with z0 first guess == 0.0001) 232 227 ztmp1 = LOG(10._wp*10000._wp) ! " " " 233 u_star = 0.035_wp*U _blk*ztmp1/ztmp0 ! (u* = 0.035*Un10)234 235 z0 = charn0 *u_star*u_star/grav + 0.11_wp*znu_a/u_star228 u_star = 0.035_wp*Ubzu*ztmp1/ztmp0 ! (u* = 0.035*Un10) 229 230 z0 = charn0_ecmwf*u_star*u_star/grav + 0.11_wp*znu_a/u_star 236 231 z0 = MIN( MAX(ABS(z0), 1.E-9) , 1._wp ) ! (prevents FPE from stupid values from masked region later on) 237 232 … … 239 234 z0t = MIN( MAX(ABS(z0t), 1.E-9) , 1._wp ) ! (prevents FPE from stupid values from masked region later on) 240 235 241 Cd = (vkarmn/ztmp0)**2! first guess of Cd242 243 ztmp0 = vkarmn *vkarmn/LOG(zt/z0t)/Cd244 245 ztmp2 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U _blk) ! Bulk Richardson Number (BRN)236 Cd = MAX( (vkarmn/ztmp0)**2 , Cx_min ) ! first guess of Cd 237 238 ztmp0 = vkarmn2/LOG(zt/z0t)/Cd 239 240 ztmp2 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, Ubzu ) ! Bulk Richardson Number (BRN) 246 241 247 242 !! First estimate of zeta_u, depending on the stability, ie sign of BRN (ztmp2): 248 243 ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 ) 249 func_m = ztmp0*ztmp2 ! temporary array !! 250 func_h = (1._wp-ztmp1) * (func_m/(1._wp+ztmp2/(-zu/(zi0*0.004_wp*Beta0**3)))) & ! BRN < 0 ! temporary array !!! func_h == zeta_u 251 & + ztmp1 * (func_m*(1._wp + 27._wp/9._wp*ztmp2/func_m)) ! BRN > 0 252 !#LB: should make sure that the "func_m" of "27./9.*ztmp2/func_m" is "ztmp0*ztmp2" and not "ztmp0==vkarmn*vkarmn/LOG(zt/z0t)/Cd" ! 253 244 func_h = (1._wp - ztmp1) * ztmp0*ztmp2 / (1._wp - ztmp2*zi0*0.004_wp*Beta0**3/zu) & ! BRN < 0 245 & + ztmp1 * ( ztmp0*ztmp2 + 27._wp/9._wp*ztmp2*ztmp2 ) ! BRN > 0 246 254 247 !! First guess M-O stability dependent scaling params.(u*,t*,q*) to estimate z0 and z/L 255 248 ztmp0 = vkarmn/(LOG(zu/z0t) - psi_h_ecmwf(func_h)) 256 249 257 u_star = MAX ( U _blk*vkarmn/(LOG(zu) - LOG(z0) - psi_m_ecmwf(func_h)) , 1.E-9 ) ! (MAX => prevents FPE from stupid values from masked region later on)250 u_star = MAX ( Ubzu*vkarmn/(LOG(zu) - LOG(z0) - psi_m_ecmwf(func_h)) , 1.E-9 ) ! (MAX => prevents FPE from stupid values from masked region later on) 258 251 t_star = dt_zu*ztmp0 259 252 q_star = dq_zu*ztmp0 … … 276 269 277 270 278 !! First guess of inverse of Monin-Obukov length (1/L) :271 !! First guess of inverse of Obukov length (1/L) : 279 272 Linv = One_on_L( t_zu, q_zu, u_star, t_star, q_star ) 280 273 281 !! Functions such as u* = U _blk*vkarmn/func_m274 !! Functions such as u* = Ubzu*vkarmn/func_m 282 275 ztmp0 = zu*Linv 283 276 func_m = LOG(zu) - LOG(z0) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf( z0*Linv) … … 285 278 286 279 !! ITERATION BLOCK 287 DO j _itt = 1, nb_itt280 DO jit = 1, nbit 288 281 289 282 !! Bulk Richardson Number at z=zu (Eq. 3.25) 290 ztmp0 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U _blk) ! Bulk Richardson Number (BRN)291 292 !! New estimate of the inverse of the Monin-Obukhon length (Linv == zeta/zu) :283 ztmp0 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, Ubzu ) ! Bulk Richardson Number (BRN) 284 285 !! New estimate of the inverse of the Obukhon length (Linv == zeta/zu) : 293 286 Linv = ztmp0*func_m*func_m/func_h / zu ! From Eq. 3.23, Chap.3.2.3, IFS doc - Cy40r1 294 287 !! Note: it is slightly different that the L we would get with the usual … … 299 292 300 293 !! Need to update roughness lengthes: 301 u_star = U _blk*vkarmn/func_m294 u_star = Ubzu*vkarmn/func_m 302 295 ztmp2 = u_star*u_star 303 296 ztmp1 = znu_a/u_star 304 z0 = MIN( ABS( alpha_M*ztmp1 + charn0 *ztmp2/grav ) , 0.001_wp)305 z0t = MIN( ABS( alpha_H*ztmp1 ) , 0.001_wp) ! eq.3.26, Chap.3, p.34, IFS doc - Cy31r1306 z0q = MIN( ABS( alpha_Q*ztmp1 ) , 0.001_wp)297 z0 = MIN( ABS( alpha_M*ztmp1 + charn0_ecmwf*ztmp2/grav ) , 0.001_wp) 298 z0t = MIN( ABS( alpha_H*ztmp1 ) , 0.001_wp) ! eq.3.26, Chap.3, p.34, IFS doc - Cy31r1 299 z0q = MIN( ABS( alpha_Q*ztmp1 ) , 0.001_wp) 307 300 308 301 !! Update wind at zu with convection-related wind gustiness in unstable conditions (Chap. 3.2, IFS doc - Cy40r1, Eq.3.17 and Eq.3.18 + Eq.3.8) 309 302 ztmp2 = Beta0*Beta0*ztmp2*(MAX(-zi0*Linv/vkarmn,0._wp))**(2._wp/3._wp) ! square of wind gustiness contribution (combining Eq. 3.8 and 3.18, hap.3, IFS doc - Cy31r1) 310 303 !! ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before zi0 311 U _blk= MAX(SQRT(U_zu*U_zu + ztmp2), 0.2_wp) ! include gustiness in bulk wind speed312 ! => 0.2 prevents U _blkto be 0 in stable case when U_zu=0.304 Ubzu = MAX(SQRT(U_zu*U_zu + ztmp2), 0.2_wp) ! include gustiness in bulk wind speed 305 ! => 0.2 prevents Ubzu to be 0 in stable case when U_zu=0. 313 306 314 307 … … 346 339 !! Cool-skin contribution 347 340 348 CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, U _blk, slp, rad_lw, &341 CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, & 349 342 & ztmp1, ztmp0, Qlat=ztmp2) ! Qnsol -> ztmp1 / Tau -> ztmp0 350 343 … … 359 352 IF( l_use_wl ) THEN 360 353 !! Warm-layer contribution 361 CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, U _blk, slp, rad_lw, &354 CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, & 362 355 & ztmp1, ztmp2) ! Qnsol -> ztmp1 / Tau -> ztmp2 363 356 CALL WL_ECMWF( Qsw, ztmp1, u_star, zsst ) … … 373 366 ENDIF 374 367 375 END DO !DO j _itt = 1, nb_itt376 377 Cd = vkarmn*vkarmn/(func_m*func_m)378 Ch = vkarmn*vkarmn/(func_m*func_h)379 ztmp2 = log(zu/z0q) - psi_h_ecmwf(zu*Linv) + psi_h_ecmwf(z0q*Linv) ! func_q380 Ce = vkarmn*vkarmn/(func_m*ztmp2)381 382 Cdn = vkarmn*vkarmn / (log(zu/z0 )*log(zu/z0 ))383 Chn = vkarmn*vkarmn / (log(zu/z0t)*log(zu/z0t))384 Cen = vkarmn*vkarmn / (log(zu/z0q)*log(zu/z0q))368 END DO !DO jit = 1, nbit 369 370 Cd = MAX( vkarmn2/(func_m*func_m) , Cx_min ) 371 Ch = MAX( vkarmn2/(func_m*func_h) , Cx_min ) 372 ztmp2 = LOG(zu/z0q) - psi_h_ecmwf(zu*Linv) + psi_h_ecmwf(z0q*Linv) ! func_q 373 Ce = MAX( vkarmn2/(func_m*ztmp2) , Cx_min ) 374 375 IF(PRESENT(Cdn)) Cdn = MAX( vkarmn2 / (LOG(zu/z0 )*LOG(zu/z0 )) , Cx_min ) 376 IF(PRESENT(Chn)) Chn = MAX( vkarmn2 / (LOG(zu/z0t)*LOG(zu/z0t)) , Cx_min ) 377 IF(PRESENT(Cen)) Cen = MAX( vkarmn2 / (LOG(zu/z0q)*LOG(zu/z0q)) , Cx_min ) 385 378 386 379 IF( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = dT_cs … … 408 401 ! 409 402 INTEGER :: ji, jj ! dummy loop indices 410 REAL(wp) :: zzeta, zx, ztmp, psi_unst, psi_stab, stab 411 !!---------------------------------------------------------------------------------- 403 REAL(wp) :: zta, zx2, zx, ztmp, zpsi_unst, zpsi_stab, zstab, zc 404 !!---------------------------------------------------------------------------------- 405 zc = 5._wp/0.35_wp 412 406 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 413 ! 414 zzeta = MIN( pzeta(ji,jj) , 5._wp ) !! Very stable conditions (L positif and big!): 415 ! 416 ! Unstable (Paulson 1970): 417 ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1 418 zx = SQRT(ABS(1._wp - 16._wp*zzeta)) 419 ztmp = 1._wp + SQRT(zx) 420 ztmp = ztmp*ztmp 421 psi_unst = LOG( 0.125_wp*ztmp*(1._wp + zx) ) & 422 & -2._wp*ATAN( SQRT(zx) ) + 0.5_wp*rpi 423 ! 424 ! Unstable: 425 ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1 426 psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) & 427 & - zzeta - 2._wp/3._wp*5._wp/0.35_wp 428 ! 429 ! Combining: 430 stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1 431 ! 432 psi_m_ecmwf(ji,jj) = (1._wp - stab) * psi_unst & ! (zzeta < 0) Unstable 433 & + stab * psi_stab ! (zzeta > 0) Stable 434 ! 407 ! 408 zta = MIN( pzeta(ji,jj) , 5._wp ) !! Very stable conditions (L positif and big!): 409 410 ! *** Unstable (Paulson 1970) [eq.3.20, Chap.3, p.33, IFS doc - Cy31r1] : 411 zx2 = SQRT( ABS(1._wp - 16._wp*zta) ) ! (1 - 16z)^0.5 412 zx = SQRT(zx2) ! (1 - 16z)^0.25 413 ztmp = 1._wp + zx 414 zpsi_unst = LOG( 0.125_wp*ztmp*ztmp*(1._wp + zx2) ) - 2._wp*ATAN( zx ) + 0.5_wp*rpi 415 416 ! *** Stable [eq.3.22, Chap.3, p.33, IFS doc - Cy31r1] : 417 zpsi_stab = -2._wp/3._wp*(zta - zc)*EXP(-0.35_wp*zta) & 418 & - zta - 2._wp/3._wp*zc 419 ! 420 zstab = 0.5_wp + SIGN(0.5_wp, zta) ! zta > 0 => zstab = 1 421 ! 422 psi_m_ecmwf(ji,jj) = zstab * zpsi_stab & ! (zta > 0) Stable 423 & + (1._wp - zstab) * zpsi_unst ! (zta < 0) Unstable 424 ! 435 425 END_2D 436 426 END FUNCTION psi_m_ecmwf … … 452 442 ! 453 443 INTEGER :: ji, jj ! dummy loop indices 454 REAL(wp) :: zzeta, zx, psi_unst, psi_stab, stab 455 !!---------------------------------------------------------------------------------- 444 REAL(wp) :: zta, zx2, zpsi_unst, zpsi_stab, zstab, zc 445 !!---------------------------------------------------------------------------------- 446 zc = 5._wp/0.35_wp 456 447 ! 457 448 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 458 ! 459 zzeta = MIN(pzeta(ji,jj) , 5._wp) ! Very stable conditions (L positif and big!): 460 ! 461 zx = ABS(1._wp - 16._wp*zzeta)**.25 ! this is actually (1/phi_m)**2 !!! 462 ! ! eq.3.19, Chap.3, p.33, IFS doc - Cy31r1 463 ! Unstable (Paulson 1970) : 464 psi_unst = 2._wp*LOG(0.5_wp*(1._wp + zx*zx)) ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1 465 ! 466 ! Stable: 467 psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) & ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1 468 & - ABS(1._wp + 2._wp/3._wp*zzeta)**1.5_wp - 2._wp/3._wp*5._wp/0.35_wp + 1._wp 469 ! LB: added ABS() to avoid NaN values when unstable, which contaminates the unstable solution... 470 ! 471 stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1 472 ! 473 ! 474 psi_h_ecmwf(ji,jj) = (1._wp - stab) * psi_unst & ! (zzeta < 0) Unstable 475 & + stab * psi_stab ! (zzeta > 0) Stable 476 ! 449 ! 450 zta = MIN(pzeta(ji,jj) , 5._wp) ! Very stable conditions (L positif and big!): 451 ! 452 ! *** Unstable (Paulson 1970) [eq.3.20, Chap.3, p.33, IFS doc - Cy31r1] : 453 zx2 = SQRT( ABS(1._wp - 16._wp*zta) ) ! (1 -16z)^0.5 454 zpsi_unst = 2._wp*LOG( 0.5_wp*(1._wp + zx2) ) 455 ! 456 ! *** Stable [eq.3.22, Chap.3, p.33, IFS doc - Cy31r1] : 457 zpsi_stab = -2._wp/3._wp*(zta - zc)*EXP(-0.35_wp*zta) & 458 & - ABS(1._wp + 2._wp/3._wp*zta)**1.5_wp - 2._wp/3._wp*zc + 1._wp 459 ! 460 ! LB: added ABS() to avoid NaN values when unstable, which contaminates the unstable solution... 461 ! 462 zstab = 0.5_wp + SIGN(0.5_wp, zta) ! zta > 0 => zstab = 1 463 ! 464 psi_h_ecmwf(ji,jj) = zstab * zpsi_stab & ! (zta > 0) Stable 465 & + (1._wp - zstab) * zpsi_unst ! (zta < 0) Unstable 466 ! 477 467 END_2D 478 468 END FUNCTION psi_h_ecmwf
Note: See TracChangeset
for help on using the changeset viewer.