Changeset 13655 for NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbcblk_algo_coare3p6.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_coare3p6.F90
r13460 r13655 7 7 !! * bulk transfer coefficients C_D, C_E and C_H 8 8 !! * air temp. and spec. hum. adjusted from zt (2m) to zu (10m) if needed 9 !! * the effective bulk wind speed at 10m U _blk9 !! * the effective bulk wind speed at 10m Ubzu 10 10 !! => all these are used in bulk formulas in sbcblk.F90 11 11 !! … … 23 23 !! returns the effective bulk wind speed at 10m 24 24 !!---------------------------------------------------------------------- 25 USE oce ! ocean dynamics and tracers26 25 USE dom_oce ! ocean space and time domain 27 26 USE phycst ! physical constants 28 USE iom ! I/O manager library 29 USE lib_mpp ! distribued memory computing library 30 USE in_out_manager ! I/O manager 31 USE prtctl ! Print control 32 USE sbcwave, ONLY : cdn_wave ! wave module 33 #if defined key_si3 || defined key_cice 34 USE sbc_ice ! Surface boundary condition: ice fields 35 #endif 36 USE lib_fortran ! to use key_nosignedzero 37 38 USE sbc_oce ! Surface boundary condition: ocean fields 39 USE sbcblk_phy ! all thermodynamics functions, rho_air, q_sat, etc... !LB 27 USE lib_mpp, ONLY: ctl_stop ! distribued memory computing library 28 USE in_out_manager, ONLY: nit000 ! I/O manager 29 USE sbc_phy ! all thermodynamics functions, rho_air, q_sat, etc... !LB 40 30 USE sbcblk_skin_coare ! cool-skin/warm layer scheme (CSWL_ECMWF) !LB 41 31 … … 50 40 REAL(wp), PARAMETER :: zi0 = 600._wp ! scale height of the atmospheric boundary layer... 51 41 REAL(wp), PARAMETER :: Beta0 = 1.2_wp ! gustiness parameter 52 53 INTEGER , PARAMETER :: nb_itt = 10 ! number of itterations 42 REAL(wp), PARAMETER :: zeta_abs_max = 50._wp 54 43 55 44 !!---------------------------------------------------------------------- … … 90 79 91 80 SUBROUTINE turb_coare3p6( kt, zt, zu, T_s, t_zt, q_s, q_zt, U_zu, l_use_cs, l_use_wl, & 92 & Cd, Ch, Ce, t_zu, q_zu, U _blk,&93 & Cdn, Chn, Cen, &81 & Cd, Ch, Ce, t_zu, q_zu, Ubzu, & 82 & nb_iter, Cdn, Chn, Cen, & ! optional output 94 83 & Qsw, rad_lw, slp, pdT_cs, & ! optionals for cool-skin (and warm-layer) 95 & pdT_wl, pHz_wl ) 84 & pdT_wl, pHz_wl ) ! optionals for warm-layer only 96 85 !!---------------------------------------------------------------------- 97 86 !! *** ROUTINE turb_coare3p6 *** … … 147 136 !! * t_zu : pot. air temperature adjusted at wind height zu [K] 148 137 !! * q_zu : specific humidity of air // [kg/kg] 149 !! * U _blk: bulk wind speed at zu [m/s]138 !! * Ubzu : bulk wind speed at zu [m/s] 150 139 !! 151 140 !! … … 167 156 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: t_zu ! pot. air temp. adjusted at zu [K] 168 157 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: q_zu ! spec. humidity adjusted at zu [kg/kg] 169 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: U_blk ! bulk wind speed at zu [m/s] 170 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Cdn, Chn, Cen ! neutral transfer coefficients 171 ! 158 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Ubzu ! bulk wind speed at zu [m/s] 159 ! 160 INTEGER , INTENT(in ), OPTIONAL :: nb_iter ! number of iterations 161 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: CdN 162 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: ChN 163 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: CeN 172 164 REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: Qsw ! [W/m^2] 173 165 REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: rad_lw ! [W/m^2] … … 177 169 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: pHz_wl ! [m] 178 170 ! 179 INTEGER :: j_itt171 INTEGER :: nbit, jit 180 172 LOGICAL :: l_zt_equal_zu = .FALSE. ! if q and t are given at same height as U 181 173 ! … … 194 186 IF( kt == nit000 ) CALL SBCBLK_ALGO_COARE3P6_INIT(l_use_cs, l_use_wl) 195 187 188 nbit = nb_iter0 189 IF( PRESENT(nb_iter) ) nbit = nb_iter 190 196 191 l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp ) ! testing "zu == zt" is risky with double precision 197 192 IF( .NOT. l_zt_equal_zu ) ALLOCATE( zeta_t(jpi,jpj) ) … … 211 206 ENDIF 212 207 213 214 208 !! First guess of temperature and humidity at height zu: 215 209 t_zu = MAX( t_zt , 180._wp ) ! who knows what's given on masked-continental regions... … … 222 216 znu_a = visc_air(t_zu) ! Air viscosity (m^2/s) at zt given from temperature in (K) 223 217 224 U _blk= SQRT(U_zu*U_zu + 0.5_wp*0.5_wp) ! initial guess for wind gustiness contribution218 Ubzu = SQRT(U_zu*U_zu + 0.5_wp*0.5_wp) ! initial guess for wind gustiness contribution 225 219 226 220 ztmp0 = LOG( zu*10000._wp) ! optimization: 10000. == 1/z0 (with z0 first guess == 0.0001) 227 221 ztmp1 = LOG(10._wp*10000._wp) ! " " " 228 u_star = 0.035_wp*U _blk*ztmp1/ztmp0 ! (u* = 0.035*Un10)229 230 z0 = alfa_charn_3p6(U_zu)*u_star*u_star/grav + 0.11_wp*znu_a/u_star222 u_star = 0.035_wp*Ubzu*ztmp1/ztmp0 ! (u* = 0.035*Un10) 223 224 z0 = charn_coare3p6(U_zu)*u_star*u_star/grav + 0.11_wp*znu_a/u_star 231 225 z0 = MIN( MAX(ABS(z0), 1.E-9) , 1._wp ) ! (prevents FPE from stupid values from masked region later on) 232 226 … … 234 228 z0t = MIN( MAX(ABS(z0t), 1.E-9) , 1._wp ) ! (prevents FPE from stupid values from masked region later on) 235 229 236 Cd = (vkarmn/ztmp0)**2! first guess of Cd237 238 ztmp0 = vkarmn *vkarmn/LOG(zt/z0t)/Cd239 240 ztmp2 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U _blk) ! Bulk Richardson Number (BRN)230 Cd = MAX( (vkarmn/ztmp0)**2 , Cx_min ) ! first guess of Cd 231 232 ztmp0 = vkarmn2/LOG(zt/z0t)/Cd 233 234 ztmp2 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, Ubzu ) ! Bulk Richardson Number (BRN) 241 235 242 236 !! First estimate of zeta_u, depending on the stability, ie sign of BRN (ztmp2): 243 237 ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 ) 244 ztmp0 = ztmp0*ztmp2 245 zeta_u = (1._wp-ztmp1) * (ztmp0/(1._wp+ztmp2/(-zu/(zi0*0.004_wp*Beta0**3)))) & ! BRN < 0 246 & + ztmp1 * (ztmp0*(1._wp + 27._wp/9._wp*ztmp2/ztmp0)) ! BRN > 0 247 !#LB: should make sure that the "ztmp0" of "27./9.*ztmp2/ztmp0" is "ztmp0*ztmp2" and not "ztmp0==vkarmn*vkarmn/LOG(zt/z0t)/Cd" ! 248 238 zeta_u = (1._wp - ztmp1) * ztmp0*ztmp2 / (1._wp - ztmp2*zi0*0.004_wp*Beta0**3/zu) & ! BRN < 0 239 & + ztmp1 * ( ztmp0*ztmp2 + 27._wp/9._wp*ztmp2*ztmp2 ) ! BRN > 0 240 249 241 !! First guess M-O stability dependent scaling params.(u*,t*,q*) to estimate z0 and z/L 250 242 ztmp0 = vkarmn/(LOG(zu/z0t) - psi_h_coare(zeta_u)) 251 243 252 u_star = MAX ( U _blk*vkarmn/(LOG(zu) - LOG(z0) - psi_m_coare(zeta_u)) , 1.E-9 ) ! (MAX => prevents FPE from stupid values from masked region later on)244 u_star = MAX ( Ubzu*vkarmn/(LOG(zu) - LOG(z0) - psi_m_coare(zeta_u)) , 1.E-9 ) ! (MAX => prevents FPE from stupid values from masked region later on) 253 245 t_star = dt_zu*ztmp0 254 246 q_star = dq_zu*ztmp0 … … 269 261 270 262 !! ITERATION BLOCK 271 DO j _itt = 1, nb_itt272 273 !!Inverse of Monin-Obukov length (1/L) :274 ztmp0 = One_on_L(t_zu, q_zu, u_star, t_star, q_star) ! 1/L == 1/[ Monin-Obukhov length]275 ztmp0 = SIGN( MIN(ABS(ztmp0),200._wp), ztmp0 ) ! (prevents FPE from stupid values from masked region later on...)263 DO jit = 1, nbit 264 265 !!Inverse of Obukov length (1/L) : 266 ztmp0 = One_on_L(t_zu, q_zu, u_star, t_star, q_star) ! 1/L == 1/[Obukhov length] 267 ztmp0 = SIGN( MIN(ABS(ztmp0),200._wp), ztmp0 ) ! 1/L (prevents FPE from stupid values from masked region later on...) 276 268 277 269 ztmp1 = u_star*u_star ! u*^2 … … 280 272 ztmp2 = Beta0*Beta0*ztmp1*(MAX(-zi0*ztmp0/vkarmn,0._wp))**(2._wp/3._wp) ! square of wind gustiness contribution, ztmp2 == Ug^2 281 273 !! ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before zi0 282 U _blk= MAX(SQRT(U_zu*U_zu + ztmp2), 0.2_wp) ! include gustiness in bulk wind speed283 ! => 0.2 prevents U _blkto be 0 in stable case when U_zu=0.274 Ubzu = MAX(SQRT(U_zu*U_zu + ztmp2), 0.2_wp) ! include gustiness in bulk wind speed 275 ! => 0.2 prevents Ubzu to be 0 in stable case when U_zu=0. 284 276 285 277 !! Stability parameters: 286 278 zeta_u = zu*ztmp0 287 zeta_u = SIGN( MIN(ABS(zeta_u), 50.0_wp), zeta_u )279 zeta_u = SIGN( MIN(ABS(zeta_u),zeta_abs_max), zeta_u ) 288 280 IF( .NOT. l_zt_equal_zu ) THEN 289 281 zeta_t = zt*ztmp0 290 zeta_t = SIGN( MIN(ABS(zeta_t), 50.0_wp), zeta_t )282 zeta_t = SIGN( MIN(ABS(zeta_t),zeta_abs_max), zeta_t ) 291 283 ENDIF 292 284 … … 296 288 !! Roughness lengthes z0, z0t (z0q = z0t) : 297 289 ztmp2 = u_star/vkarmn*LOG(10./z0) ! Neutral wind speed at 10m 298 z0 = alfa_charn_3p6(ztmp2)*ztmp1/grav + 0.11_wp*znu_a/u_star ! Roughness length (eq.6) [ ztmp1==u*^2 ]290 z0 = charn_coare3p6(ztmp2)*ztmp1/grav + 0.11_wp*znu_a/u_star ! Roughness length (eq.6) [ ztmp1==u*^2 ] 299 291 z0 = MIN( MAX(ABS(z0), 1.E-9) , 1._wp ) ! (prevents FPE from stupid values from masked region later on) 300 292 … … 309 301 t_star = dt_zu*ztmp1 310 302 q_star = dq_zu*ztmp1 311 u_star = MAX( U _blk*vkarmn/(LOG(zu) - LOG(z0) - psi_m_coare(zeta_u)) , 1.E-9 ) ! (MAX => prevents FPE from stupid values from masked region later on)303 u_star = MAX( Ubzu*vkarmn/(LOG(zu) - LOG(z0) - psi_m_coare(zeta_u)) , 1.E-9 ) ! (MAX => prevents FPE from stupid values from masked region later on) 312 304 313 305 IF( .NOT. l_zt_equal_zu ) THEN … … 318 310 ENDIF 319 311 320 321 312 IF( l_use_cs ) THEN 322 313 !! Cool-skin contribution 323 314 324 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, &315 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, & 325 316 & ztmp1, zeta_u, Qlat=ztmp2) ! Qnsol -> ztmp1 / Tau -> zeta_u 326 317 … … 330 321 IF( l_use_wl ) T_s(:,:) = T_s(:,:) + dT_wl(:,:)*tmask(:,:,1) 331 322 q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 332 333 323 ENDIF 334 324 335 325 IF( l_use_wl ) THEN 336 326 !! Warm-layer contribution 337 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, &327 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, & 338 328 & ztmp1, zeta_u) ! Qnsol -> ztmp1 / Tau -> zeta_u 339 329 !! In WL_COARE or , Tau_ac and Qnt_ac must be updated at the final itteration step => add a flag to do this! 340 CALL WL_COARE( Qsw, ztmp1, zeta_u, zsst, MOD(nb _itt,j_itt) )330 CALL WL_COARE( Qsw, ztmp1, zeta_u, zsst, MOD(nbit,jit) ) 341 331 342 332 !! Updating T_s and q_s !!! … … 351 341 ENDIF 352 342 353 END DO !DO j _itt = 1, nb_itt343 END DO !DO jit = 1, nbit 354 344 355 345 ! compute transfer coefficients at zu : 356 ztmp0 = u_star/U_blk 357 Cd = ztmp0*ztmp0 358 Ch = ztmp0*t_star/dt_zu 359 Ce = ztmp0*q_star/dq_zu 360 361 ztmp1 = zu + z0 362 Cdn = vkarmn*vkarmn / (log(ztmp1/z0 )*log(ztmp1/z0 )) 363 Chn = vkarmn*vkarmn / (log(ztmp1/z0t)*log(ztmp1/z0t)) 364 Cen = Chn 346 ztmp0 = u_star/Ubzu 347 Cd = MAX( ztmp0*ztmp0 , Cx_min ) 348 Ch = MAX( ztmp0*t_star/dt_zu , Cx_min ) 349 Ce = MAX( ztmp0*q_star/dq_zu , Cx_min ) 365 350 366 351 IF( .NOT. l_zt_equal_zu ) DEALLOCATE( zeta_t ) 352 353 IF(PRESENT(Cdn)) Cdn = MAX( vkarmn2 / (LOG(zu/z0 )*LOG(zu/z0 )) , Cx_min ) 354 IF(PRESENT(Chn)) Chn = MAX( vkarmn2 / (LOG(zu/z0t)*LOG(zu/z0t)) , Cx_min ) 355 IF(PRESENT(Cen)) Cen = MAX( vkarmn2 / (LOG(zu/z0t)*LOG(zu/z0t)) , Cx_min ) 367 356 368 357 IF( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = dT_cs … … 375 364 376 365 377 FUNCTION alfa_charn_3p6( pwnd )366 FUNCTION charn_coare3p6( pwnd ) 378 367 !!------------------------------------------------------------------- 379 368 !! Computes the Charnock parameter as a function of the Neutral wind speed at 10m … … 383 372 !! Author: L. Brodeau, July 2019 / AeroBulk (https://github.com/brodeau/aerobulk/) 384 373 !!------------------------------------------------------------------- 385 REAL(wp), DIMENSION(jpi,jpj) :: alfa_charn_3p6374 REAL(wp), DIMENSION(jpi,jpj) :: charn_coare3p6 386 375 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pwnd ! neutral wind speed at 10m 387 376 ! 388 377 REAL(wp), PARAMETER :: charn0_max = 0.028 !: value above which the Charnock parameter levels off for winds > 18 m/s 389 378 !!------------------------------------------------------------------- 390 alfa_charn_3p6 = MAX( MIN( 0.0017_wp*pwnd - 0.005_wp , charn0_max) , 0._wp )391 !! 392 END FUNCTION alfa_charn_3p6393 394 FUNCTION alfa_charn_3p6_wave( pus, pwsh, pwps )379 charn_coare3p6 = MAX( MIN( 0.0017_wp*pwnd - 0.005_wp , charn0_max) , 0._wp ) 380 !! 381 END FUNCTION charn_coare3p6 382 383 FUNCTION charn_coare3p6_wave( pus, pwsh, pwps ) 395 384 !!------------------------------------------------------------------- 396 385 !! Computes the Charnock parameter as a function of wave information and u* … … 400 389 !! Author: L. Brodeau, October 2019 / AeroBulk (https://github.com/brodeau/aerobulk/) 401 390 !!------------------------------------------------------------------- 402 REAL(wp), DIMENSION(jpi,jpj) :: alfa_charn_3p6_wave391 REAL(wp), DIMENSION(jpi,jpj) :: charn_coare3p6_wave 403 392 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pus ! friction velocity [m/s] 404 393 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pwsh ! significant wave height [m] 405 394 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pwps ! phase speed of dominant waves [m/s] 406 395 !!------------------------------------------------------------------- 407 alfa_charn_3p6_wave = ( pwsh*0.2_wp*(pus/pwps)**2.2_wp ) * grav/(pus*pus)408 !! 409 END FUNCTION alfa_charn_3p6_wave396 charn_coare3p6_wave = ( pwsh*0.2_wp*(pus/pwps)**2.2_wp ) * grav/(pus*pus) 397 !! 398 END FUNCTION charn_coare3p6_wave 410 399 411 400 … … 429 418 REAL(wp) :: zta, zphi_m, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab 430 419 !!---------------------------------------------------------------------------------- 431 !432 420 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 433 ! 434 zta = pzeta(ji,jj) 435 ! 436 zphi_m = ABS(1. - 15.*zta)**.25 !!Kansas unstable 437 ! 438 zpsi_k = 2.*LOG((1. + zphi_m)/2.) + LOG((1. + zphi_m*zphi_m)/2.) & 439 & - 2.*ATAN(zphi_m) + 0.5*rpi 440 ! 441 zphi_c = ABS(1. - 10.15*zta)**.3333 !!Convective 442 ! 443 zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) & 444 & - 1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447 445 ! 446 zf = zta*zta 447 zf = zf/(1. + zf) 448 zc = MIN(50._wp, 0.35_wp*zta) 449 zstab = 0.5 + SIGN(0.5_wp, zta) 450 ! 451 psi_m_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & ! (zta < 0) 452 & - zstab * ( 1. + 1.*zta & ! (zta > 0) 453 & + 0.6667*(zta - 14.28)/EXP(zc) + 8.525 ) ! " 454 ! 421 ! 422 zta = pzeta(ji,jj) 423 ! 424 zphi_m = ABS(1. - 15.*zta)**.25 !!Kansas unstable 425 ! 426 zpsi_k = 2.*LOG((1. + zphi_m)/2.) + LOG((1. + zphi_m*zphi_m)/2.) & 427 & - 2.*ATAN(zphi_m) + 0.5*rpi 428 ! 429 zphi_c = ABS(1. - 10.15*zta)**.3333 !!Convective 430 ! 431 zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) & 432 & - 1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447 433 ! 434 zf = zta*zta 435 zf = zf/(1. + zf) 436 zc = MIN(50._wp, 0.35_wp*zta) 437 zstab = 0.5 + SIGN(0.5_wp, zta) 438 ! 439 psi_m_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & ! (zta < 0) 440 & - zstab * ( 1. + 1.*zta & ! (zta > 0) 441 & + 0.6667*(zta - 14.28)/EXP(zc) + 8.525 ) ! " 455 442 END_2D 456 !457 443 END FUNCTION psi_m_coare 458 444 … … 474 460 !! (https://github.com/brodeau/aerobulk/) 475 461 !!---------------------------------------------------------------- 476 !!477 462 REAL(wp), DIMENSION(jpi,jpj) :: psi_h_coare 478 463 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta … … 480 465 INTEGER :: ji, jj ! dummy loop indices 481 466 REAL(wp) :: zta, zphi_h, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab 482 ! 467 !!---------------------------------------------------------------- 483 468 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 484 ! 485 zta = pzeta(ji,jj) 486 ! 487 zphi_h = (ABS(1. - 15.*zta))**.5 !! Kansas unstable (zphi_h = zphi_m**2 when unstable, zphi_m when stable) 488 ! 489 zpsi_k = 2.*LOG((1. + zphi_h)/2.) 490 ! 491 zphi_c = (ABS(1. - 34.15*zta))**.3333 !! Convective 492 ! 493 zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) & 494 & -1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447 495 ! 496 zf = zta*zta 497 zf = zf/(1. + zf) 498 zc = MIN(50._wp,0.35_wp*zta) 499 zstab = 0.5 + SIGN(0.5_wp, zta) 500 ! 501 psi_h_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & 502 & - zstab * ( (ABS(1. + 2.*zta/3.))**1.5 & 503 & + .6667*(zta - 14.28)/EXP(zc) + 8.525 ) 504 ! 469 ! 470 zta = pzeta(ji,jj) 471 ! 472 zphi_h = (ABS(1. - 15.*zta))**.5 !! Kansas unstable (zphi_h = zphi_m**2 when unstable, zphi_m when stable) 473 ! 474 zpsi_k = 2.*LOG((1. + zphi_h)/2.) 475 ! 476 zphi_c = (ABS(1. - 34.15*zta))**.3333 !! Convective 477 ! 478 zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) & 479 & -1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447 480 ! 481 zf = zta*zta 482 zf = zf/(1. + zf) 483 zc = MIN(50._wp,0.35_wp*zta) 484 zstab = 0.5 + SIGN(0.5_wp, zta) 485 ! 486 psi_h_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & 487 & - zstab * ( (ABS(1. + 2.*zta/3.))**1.5 & 488 & + .6667*(zta - 14.28)/EXP(zc) + 8.525 ) 505 489 END_2D 506 !507 490 END FUNCTION psi_h_coare 508 491
Note: See TracChangeset
for help on using the changeset viewer.