Changeset 14789 for NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU/src/OCE/SBC/sbcblk_algo_coare3p6.F90
- Timestamp:
- 2021-05-05T13:18:04+02:00 (3 years ago)
- Location:
- NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/SBC/sbcblk_algo_coare3p6.F90 (modified) (25 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev _r12970_AGRIF_CMEMSext/AGRIF5 ^/vendors/AGRIF/dev@HEAD ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 ^/vendors/PPR@HEAD ext/PPR 8 9 9 10 # SETTE 10 ^/utils/CI/sette@1 3559sette11 ^/utils/CI/sette@14244 sette
-
- Property svn:externals
-
NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU/src/OCE/SBC/sbcblk_algo_coare3p6.F90
r13460 r14789 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 !! … … 16 16 !!---------------------------------------------------------------------- 17 17 !! History : 4.0 ! 2016-02 (L.Brodeau) Original code 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 27 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 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 ! Catalog of functions for physical/meteorological parameters in the marine boundary layer 40 31 USE sbcblk_skin_coare ! cool-skin/warm layer scheme (CSWL_ECMWF) !LB 41 32 … … 50 41 REAL(wp), PARAMETER :: zi0 = 600._wp ! scale height of the atmospheric boundary layer... 51 42 REAL(wp), PARAMETER :: Beta0 = 1.2_wp ! gustiness parameter 52 53 INTEGER , PARAMETER :: nb_itt = 10 ! number of itterations 43 REAL(wp), PARAMETER :: zeta_abs_max = 50._wp 54 44 55 45 !!---------------------------------------------------------------------- … … 90 80 91 81 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, &82 & Cd, Ch, Ce, t_zu, q_zu, Ubzu, & 83 & nb_iter, Cdn, Chn, Cen, & ! optional output 94 84 & Qsw, rad_lw, slp, pdT_cs, & ! optionals for cool-skin (and warm-layer) 95 & pdT_wl, pHz_wl ) 85 & pdT_wl, pHz_wl ) ! optionals for warm-layer only 96 86 !!---------------------------------------------------------------------- 97 87 !! *** ROUTINE turb_coare3p6 *** … … 147 137 !! * t_zu : pot. air temperature adjusted at wind height zu [K] 148 138 !! * q_zu : specific humidity of air // [kg/kg] 149 !! * U _blk: bulk wind speed at zu [m/s]139 !! * Ubzu : bulk wind speed at zu [m/s] 150 140 !! 151 141 !! … … 167 157 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: t_zu ! pot. air temp. adjusted at zu [K] 168 158 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 ! 159 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Ubzu ! bulk wind speed at zu [m/s] 160 ! 161 INTEGER , INTENT(in ), OPTIONAL :: nb_iter ! number of iterations 162 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: CdN 163 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: ChN 164 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: CeN 172 165 REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: Qsw ! [W/m^2] 173 166 REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: rad_lw ! [W/m^2] … … 177 170 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: pHz_wl ! [m] 178 171 ! 179 INTEGER :: j_itt172 INTEGER :: nbit, jit 180 173 LOGICAL :: l_zt_equal_zu = .FALSE. ! if q and t are given at same height as U 181 174 ! … … 194 187 IF( kt == nit000 ) CALL SBCBLK_ALGO_COARE3P6_INIT(l_use_cs, l_use_wl) 195 188 189 nbit = nb_iter0 190 IF( PRESENT(nb_iter) ) nbit = nb_iter 191 196 192 l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp ) ! testing "zu == zt" is risky with double precision 197 193 IF( .NOT. l_zt_equal_zu ) ALLOCATE( zeta_t(jpi,jpj) ) … … 211 207 ENDIF 212 208 213 214 209 !! First guess of temperature and humidity at height zu: 215 210 t_zu = MAX( t_zt , 180._wp ) ! who knows what's given on masked-continental regions... … … 222 217 znu_a = visc_air(t_zu) ! Air viscosity (m^2/s) at zt given from temperature in (K) 223 218 224 U _blk= SQRT(U_zu*U_zu + 0.5_wp*0.5_wp) ! initial guess for wind gustiness contribution219 Ubzu = SQRT(U_zu*U_zu + 0.5_wp*0.5_wp) ! initial guess for wind gustiness contribution 225 220 226 221 ztmp0 = LOG( zu*10000._wp) ! optimization: 10000. == 1/z0 (with z0 first guess == 0.0001) 227 222 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_star223 u_star = 0.035_wp*Ubzu*ztmp1/ztmp0 ! (u* = 0.035*Un10) 224 225 z0 = charn_coare3p6(U_zu)*u_star*u_star/grav + 0.11_wp*znu_a/u_star 231 226 z0 = MIN( MAX(ABS(z0), 1.E-9) , 1._wp ) ! (prevents FPE from stupid values from masked region later on) 232 227 … … 234 229 z0t = MIN( MAX(ABS(z0t), 1.E-9) , 1._wp ) ! (prevents FPE from stupid values from masked region later on) 235 230 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)231 Cd = MAX( (vkarmn/ztmp0)**2 , Cx_min ) ! first guess of Cd 232 233 ztmp0 = vkarmn2/LOG(zt/z0t)/Cd 234 235 ztmp2 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, Ubzu ) ! Bulk Richardson Number (BRN) 241 236 242 237 !! First estimate of zeta_u, depending on the stability, ie sign of BRN (ztmp2): 243 238 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" ! 239 zeta_u = (1._wp - ztmp1) * ztmp0*ztmp2 / (1._wp - ztmp2*zi0*0.004_wp*Beta0**3/zu) & ! BRN < 0 240 & + ztmp1 * ( ztmp0*ztmp2 + 27._wp/9._wp*ztmp2*ztmp2 ) ! BRN > 0 248 241 249 242 !! First guess M-O stability dependent scaling params.(u*,t*,q*) to estimate z0 and z/L 250 243 ztmp0 = vkarmn/(LOG(zu/z0t) - psi_h_coare(zeta_u)) 251 244 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)245 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 246 t_star = dt_zu*ztmp0 254 247 q_star = dq_zu*ztmp0 … … 269 262 270 263 !! 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...)264 DO jit = 1, nbit 265 266 !!Inverse of Obukov length (1/L) : 267 ztmp0 = One_on_L(t_zu, q_zu, u_star, t_star, q_star) ! 1/L == 1/[Obukhov length] 268 ztmp0 = SIGN( MIN(ABS(ztmp0),200._wp), ztmp0 ) ! 1/L (prevents FPE from stupid values from masked region later on...) 276 269 277 270 ztmp1 = u_star*u_star ! u*^2 … … 280 273 ztmp2 = Beta0*Beta0*ztmp1*(MAX(-zi0*ztmp0/vkarmn,0._wp))**(2._wp/3._wp) ! square of wind gustiness contribution, ztmp2 == Ug^2 281 274 !! ! 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.275 Ubzu = MAX(SQRT(U_zu*U_zu + ztmp2), 0.2_wp) ! include gustiness in bulk wind speed 276 ! => 0.2 prevents Ubzu to be 0 in stable case when U_zu=0. 284 277 285 278 !! Stability parameters: 286 279 zeta_u = zu*ztmp0 287 zeta_u = SIGN( MIN(ABS(zeta_u), 50.0_wp), zeta_u )280 zeta_u = SIGN( MIN(ABS(zeta_u),zeta_abs_max), zeta_u ) 288 281 IF( .NOT. l_zt_equal_zu ) THEN 289 282 zeta_t = zt*ztmp0 290 zeta_t = SIGN( MIN(ABS(zeta_t), 50.0_wp), zeta_t )283 zeta_t = SIGN( MIN(ABS(zeta_t),zeta_abs_max), zeta_t ) 291 284 ENDIF 292 285 … … 296 289 !! Roughness lengthes z0, z0t (z0q = z0t) : 297 290 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 ]291 z0 = charn_coare3p6(ztmp2)*ztmp1/grav + 0.11_wp*znu_a/u_star ! Roughness length (eq.6) [ ztmp1==u*^2 ] 299 292 z0 = MIN( MAX(ABS(z0), 1.E-9) , 1._wp ) ! (prevents FPE from stupid values from masked region later on) 300 293 … … 309 302 t_star = dt_zu*ztmp1 310 303 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)304 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 305 313 306 IF( .NOT. l_zt_equal_zu ) THEN … … 318 311 ENDIF 319 312 320 321 313 IF( l_use_cs ) THEN 322 314 !! Cool-skin contribution 323 315 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, &316 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 317 & ztmp1, zeta_u, Qlat=ztmp2) ! Qnsol -> ztmp1 / Tau -> zeta_u 326 318 … … 330 322 IF( l_use_wl ) T_s(:,:) = T_s(:,:) + dT_wl(:,:)*tmask(:,:,1) 331 323 q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 332 333 324 ENDIF 334 325 335 326 IF( l_use_wl ) THEN 336 327 !! 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, &328 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 329 & ztmp1, zeta_u) ! Qnsol -> ztmp1 / Tau -> zeta_u 339 330 !! 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) )331 CALL WL_COARE( Qsw, ztmp1, zeta_u, zsst, MOD(nbit,jit) ) 341 332 342 333 !! Updating T_s and q_s !!! … … 351 342 ENDIF 352 343 353 END DO !DO j _itt = 1, nb_itt344 END DO !DO jit = 1, nbit 354 345 355 346 ! 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 347 ztmp0 = u_star/Ubzu 348 Cd = MAX( ztmp0*ztmp0 , Cx_min ) 349 Ch = MAX( ztmp0*t_star/dt_zu , Cx_min ) 350 Ce = MAX( ztmp0*q_star/dq_zu , Cx_min ) 365 351 366 352 IF( .NOT. l_zt_equal_zu ) DEALLOCATE( zeta_t ) 353 354 IF(PRESENT(Cdn)) Cdn = MAX( vkarmn2 / (LOG(zu/z0 )*LOG(zu/z0 )) , Cx_min ) 355 IF(PRESENT(Chn)) Chn = MAX( vkarmn2 / (LOG(zu/z0t)*LOG(zu/z0t)) , Cx_min ) 356 IF(PRESENT(Cen)) Cen = MAX( vkarmn2 / (LOG(zu/z0t)*LOG(zu/z0t)) , Cx_min ) 367 357 368 358 IF( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = dT_cs … … 375 365 376 366 377 FUNCTION alfa_charn_3p6( pwnd )367 FUNCTION charn_coare3p6( pwnd ) 378 368 !!------------------------------------------------------------------- 379 369 !! Computes the Charnock parameter as a function of the Neutral wind speed at 10m … … 383 373 !! Author: L. Brodeau, July 2019 / AeroBulk (https://github.com/brodeau/aerobulk/) 384 374 !!------------------------------------------------------------------- 385 REAL(wp), DIMENSION(jpi,jpj) :: alfa_charn_3p6375 REAL(wp), DIMENSION(jpi,jpj) :: charn_coare3p6 386 376 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pwnd ! neutral wind speed at 10m 387 377 ! 388 378 REAL(wp), PARAMETER :: charn0_max = 0.028 !: value above which the Charnock parameter levels off for winds > 18 m/s 389 379 !!------------------------------------------------------------------- 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 )380 charn_coare3p6 = MAX( MIN( 0.0017_wp*pwnd - 0.005_wp , charn0_max) , 0._wp ) 381 !! 382 END FUNCTION charn_coare3p6 383 384 FUNCTION charn_coare3p6_wave( pus, pwsh, pwps ) 395 385 !!------------------------------------------------------------------- 396 386 !! Computes the Charnock parameter as a function of wave information and u* … … 400 390 !! Author: L. Brodeau, October 2019 / AeroBulk (https://github.com/brodeau/aerobulk/) 401 391 !!------------------------------------------------------------------- 402 REAL(wp), DIMENSION(jpi,jpj) :: alfa_charn_3p6_wave392 REAL(wp), DIMENSION(jpi,jpj) :: charn_coare3p6_wave 403 393 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pus ! friction velocity [m/s] 404 394 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pwsh ! significant wave height [m] 405 395 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pwps ! phase speed of dominant waves [m/s] 406 396 !!------------------------------------------------------------------- 407 alfa_charn_3p6_wave = ( pwsh*0.2_wp*(pus/pwps)**2.2_wp ) * grav/(pus*pus)408 !! 409 END FUNCTION alfa_charn_3p6_wave397 charn_coare3p6_wave = ( pwsh*0.2_wp*(pus/pwps)**2.2_wp ) * grav/(pus*pus) 398 !! 399 END FUNCTION charn_coare3p6_wave 410 400 411 401 … … 429 419 REAL(wp) :: zta, zphi_m, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab 430 420 !!---------------------------------------------------------------------------------- 431 !432 421 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 ! 422 ! 423 zta = pzeta(ji,jj) 424 ! 425 zphi_m = ABS(1. - 15.*zta)**.25 !!Kansas unstable 426 ! 427 zpsi_k = 2.*LOG((1. + zphi_m)/2.) + LOG((1. + zphi_m*zphi_m)/2.) & 428 & - 2.*ATAN(zphi_m) + 0.5*rpi 429 ! 430 zphi_c = ABS(1. - 10.15*zta)**.3333 !!Convective 431 ! 432 zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) & 433 & - 1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447 434 ! 435 zf = zta*zta 436 zf = zf/(1. + zf) 437 zc = MIN(50._wp, 0.35_wp*zta) 438 zstab = 0.5 + SIGN(0.5_wp, zta) 439 ! 440 psi_m_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & ! (zta < 0) 441 & - zstab * ( 1. + 1.*zta & ! (zta > 0) 442 & + 0.6667*(zta - 14.28)/EXP(zc) + 8.525 ) ! " 455 443 END_2D 456 !457 444 END FUNCTION psi_m_coare 458 445 … … 474 461 !! (https://github.com/brodeau/aerobulk/) 475 462 !!---------------------------------------------------------------- 476 !!477 463 REAL(wp), DIMENSION(jpi,jpj) :: psi_h_coare 478 464 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta … … 480 466 INTEGER :: ji, jj ! dummy loop indices 481 467 REAL(wp) :: zta, zphi_h, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab 482 ! 468 !!---------------------------------------------------------------- 483 469 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 ! 470 ! 471 zta = pzeta(ji,jj) 472 ! 473 zphi_h = (ABS(1. - 15.*zta))**.5 !! Kansas unstable (zphi_h = zphi_m**2 when unstable, zphi_m when stable) 474 ! 475 zpsi_k = 2.*LOG((1. + zphi_h)/2.) 476 ! 477 zphi_c = (ABS(1. - 34.15*zta))**.3333 !! Convective 478 ! 479 zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) & 480 & -1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447 481 ! 482 zf = zta*zta 483 zf = zf/(1. + zf) 484 zc = MIN(50._wp,0.35_wp*zta) 485 zstab = 0.5 + SIGN(0.5_wp, zta) 486 ! 487 psi_h_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & 488 & - zstab * ( (ABS(1. + 2.*zta/3.))**1.5 & 489 & + .6667*(zta - 14.28)/EXP(zc) + 8.525 ) 505 490 END_2D 506 !507 491 END FUNCTION psi_h_coare 508 492
Note: See TracChangeset
for help on using the changeset viewer.
