Changeset 13655 for NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbcblk_algo_coare3p0.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_coare3p0.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 !! … … 37 37 38 38 USE sbc_oce ! Surface boundary condition: ocean fields 39 USE sbc blk_phy! all thermodynamics functions, rho_air, q_sat, etc... !LB39 USE sbc_phy ! all thermodynamics functions, rho_air, q_sat, etc... !LB 40 40 USE sbcblk_skin_coare ! cool-skin/warm layer scheme (CSWL_ECMWF) !LB 41 41 … … 50 50 REAL(wp), PARAMETER :: zi0 = 600._wp ! scale height of the atmospheric boundary layer... 51 51 REAL(wp), PARAMETER :: Beta0 = 1.25_wp ! gustiness parameter 52 53 INTEGER , PARAMETER :: nb_itt = 10 ! number of itterations 52 REAL(wp), PARAMETER :: zeta_abs_max = 50._wp 54 53 55 54 !!---------------------------------------------------------------------- … … 90 89 91 90 SUBROUTINE turb_coare3p0( 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, &91 & Cd, Ch, Ce, t_zu, q_zu, Ubzu, & 92 & nb_iter, Cdn, Chn, Cen, & ! optional output 94 93 & Qsw, rad_lw, slp, pdT_cs, & ! optionals for cool-skin (and warm-layer) 95 & pdT_wl, pHz_wl ) 94 & pdT_wl, pHz_wl ) ! optionals for warm-layer only 96 95 !!---------------------------------------------------------------------- 97 96 !! *** ROUTINE turb_coare3p0 *** … … 147 146 !! * t_zu : pot. air temperature adjusted at wind height zu [K] 148 147 !! * q_zu : specific humidity of air // [kg/kg] 149 !! * U _blk: bulk wind speed at zu [m/s]148 !! * Ubzu : bulk wind speed at zu [m/s] 150 149 !! 151 150 !! … … 167 166 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: t_zu ! pot. air temp. adjusted at zu [K] 168 167 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 ! 168 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Ubzu ! bulk wind speed at zu [m/s] 169 ! 170 INTEGER , INTENT(in ), OPTIONAL :: nb_iter ! number of iterations 171 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: CdN 172 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: ChN 173 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: CeN 172 174 REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: Qsw ! [W/m^2] 173 175 REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: rad_lw ! [W/m^2] … … 177 179 REAL(wp), INTENT( out), OPTIONAL, DIMENSION(jpi,jpj) :: pHz_wl ! [m] 178 180 ! 179 INTEGER :: j_itt181 INTEGER :: nbit, jit 180 182 LOGICAL :: l_zt_equal_zu = .FALSE. ! if q and t are given at same height as U 181 183 ! … … 194 196 IF( kt == nit000 ) CALL SBCBLK_ALGO_COARE3P0_INIT(l_use_cs, l_use_wl) 195 197 198 nbit = nb_iter0 199 IF( PRESENT(nb_iter) ) nbit = nb_iter 200 196 201 l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp ) ! testing "zu == zt" is risky with double precision 197 202 IF( .NOT. l_zt_equal_zu ) ALLOCATE( zeta_t(jpi,jpj) ) … … 211 216 ENDIF 212 217 213 214 218 !! First guess of temperature and humidity at height zu: 215 219 t_zu = MAX( t_zt , 180._wp ) ! who knows what's given on masked-continental regions... … … 222 226 znu_a = visc_air(t_zu) ! Air viscosity (m^2/s) at zt given from temperature in (K) 223 227 224 U _blk= SQRT(U_zu*U_zu + 0.5_wp*0.5_wp) ! initial guess for wind gustiness contribution228 Ubzu = SQRT(U_zu*U_zu + 0.5_wp*0.5_wp) ! initial guess for wind gustiness contribution 225 229 226 230 ztmp0 = LOG( zu*10000._wp) ! optimization: 10000. == 1/z0 (with z0 first guess == 0.0001) 227 231 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_3p0(U_zu)*u_star*u_star/grav + 0.11_wp*znu_a/u_star232 u_star = 0.035_wp*Ubzu*ztmp1/ztmp0 ! (u* = 0.035*Un10) 233 234 z0 = charn_coare3p0(U_zu)*u_star*u_star/grav + 0.11_wp*znu_a/u_star 231 235 z0 = MIN( MAX(ABS(z0), 1.E-9) , 1._wp ) ! (prevents FPE from stupid values from masked region later on) 232 236 … … 234 238 z0t = MIN( MAX(ABS(z0t), 1.E-9) , 1._wp ) ! (prevents FPE from stupid values from masked region later on) 235 239 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)240 Cd = MAX( (vkarmn/ztmp0)**2 , Cx_min ) ! first guess of Cd 241 242 ztmp0 = vkarmn2/LOG(zt/z0t)/Cd 243 244 ztmp2 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, Ubzu ) ! Bulk Richardson Number (BRN) 241 245 242 246 !! First estimate of zeta_u, depending on the stability, ie sign of BRN (ztmp2): 243 247 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 248 zeta_u = (1._wp - ztmp1) * ztmp0*ztmp2 / (1._wp - ztmp2*zi0*0.004_wp*Beta0**3/zu) & ! BRN < 0 249 & + ztmp1 * ( ztmp0*ztmp2 + 27._wp/9._wp*ztmp2*ztmp2 ) ! BRN > 0 250 249 251 !! First guess M-O stability dependent scaling params.(u*,t*,q*) to estimate z0 and z/L 250 252 ztmp0 = vkarmn/(LOG(zu/z0t) - psi_h_coare(zeta_u)) 251 253 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)254 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 255 t_star = dt_zu*ztmp0 254 256 q_star = dq_zu*ztmp0 … … 269 271 270 272 !! 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...)273 DO jit = 1, nbit 274 275 !!Inverse of Obukov length (1/L) : 276 ztmp0 = One_on_L(t_zu, q_zu, u_star, t_star, q_star) ! 1/L == 1/[Obukhov length] 277 ztmp0 = SIGN( MIN(ABS(ztmp0),200._wp), ztmp0 ) ! 1/L (prevents FPE from stupid values from masked region later on...) 276 278 277 279 ztmp1 = u_star*u_star ! u*^2 … … 280 282 ztmp2 = Beta0*Beta0*ztmp1*(MAX(-zi0*ztmp0/vkarmn,0._wp))**(2._wp/3._wp) ! square of wind gustiness contribution, ztmp2 == Ug^2 281 283 !! ! 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.284 Ubzu = MAX(SQRT(U_zu*U_zu + ztmp2), 0.2_wp) ! include gustiness in bulk wind speed 285 ! => 0.2 prevents Ubzu to be 0 in stable case when U_zu=0. 284 286 285 287 !! Stability parameters: 286 288 zeta_u = zu*ztmp0 287 zeta_u = SIGN( MIN(ABS(zeta_u), 50.0_wp), zeta_u )289 zeta_u = SIGN( MIN(ABS(zeta_u),zeta_abs_max), zeta_u ) 288 290 IF( .NOT. l_zt_equal_zu ) THEN 289 291 zeta_t = zt*ztmp0 290 zeta_t = SIGN( MIN(ABS(zeta_t), 50.0_wp), zeta_t )292 zeta_t = SIGN( MIN(ABS(zeta_t),zeta_abs_max), zeta_t ) 291 293 ENDIF 292 294 … … 296 298 !! Roughness lengthes z0, z0t (z0q = z0t) : 297 299 ztmp2 = u_star/vkarmn*LOG(10./z0) ! Neutral wind speed at 10m 298 z0 = alfa_charn_3p0(ztmp2)*ztmp1/grav + 0.11_wp*znu_a/u_star ! Roughness length (eq.6) [ ztmp1==u*^2 ]300 z0 = charn_coare3p0(ztmp2)*ztmp1/grav + 0.11_wp*znu_a/u_star ! Roughness length (eq.6) [ ztmp1==u*^2 ] 299 301 z0 = MIN( MAX(ABS(z0), 1.E-9) , 1._wp ) ! (prevents FPE from stupid values from masked region later on) 300 302 … … 309 311 t_star = dt_zu*ztmp1 310 312 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)313 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 314 313 315 IF( .NOT. l_zt_equal_zu ) THEN … … 318 320 ENDIF 319 321 320 321 322 IF( l_use_cs ) THEN 322 323 !! Cool-skin contribution 323 324 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, &325 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 326 & ztmp1, zeta_u, Qlat=ztmp2) ! Qnsol -> ztmp1 / Tau -> zeta_u 326 327 … … 330 331 IF( l_use_wl ) T_s(:,:) = T_s(:,:) + dT_wl(:,:)*tmask(:,:,1) 331 332 q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 332 333 333 ENDIF 334 334 335 335 IF( l_use_wl ) THEN 336 336 !! 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, &337 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 338 & ztmp1, zeta_u) ! Qnsol -> ztmp1 / Tau -> zeta_u 339 339 !! 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) )340 CALL WL_COARE( Qsw, ztmp1, zeta_u, zsst, MOD(nbit,jit) ) 341 341 342 342 !! Updating T_s and q_s !!! … … 351 351 ENDIF 352 352 353 END DO !DO j _itt = 1, nb_itt353 END DO !DO jit = 1, nbit 354 354 355 355 ! 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 356 ztmp0 = u_star/Ubzu 357 Cd = MAX( ztmp0*ztmp0 , Cx_min ) 358 Ch = MAX( ztmp0*t_star/dt_zu , Cx_min ) 359 Ce = MAX( ztmp0*q_star/dq_zu , Cx_min ) 365 360 366 361 IF( .NOT. l_zt_equal_zu ) DEALLOCATE( zeta_t ) 362 363 IF(PRESENT(Cdn)) Cdn = MAX( vkarmn2 / (LOG(zu/z0 )*LOG(zu/z0 )) , Cx_min ) 364 IF(PRESENT(Chn)) Chn = MAX( vkarmn2 / (LOG(zu/z0t)*LOG(zu/z0t)) , Cx_min ) 365 IF(PRESENT(Cen)) Cen = MAX( vkarmn2 / (LOG(zu/z0t)*LOG(zu/z0t)) , Cx_min ) 367 366 368 367 IF( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = dT_cs … … 375 374 376 375 377 FUNCTION alfa_charn_3p0( pwnd )376 FUNCTION charn_coare3p0( pwnd ) 378 377 !!------------------------------------------------------------------- 379 378 !! Compute the Charnock parameter as a function of the wind speed … … 387 386 !! Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 388 387 !!------------------------------------------------------------------- 389 REAL(wp), DIMENSION(jpi,jpj) :: alfa_charn_3p0388 REAL(wp), DIMENSION(jpi,jpj) :: charn_coare3p0 390 389 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pwnd ! wind speed 391 390 ! … … 393 392 REAL(wp) :: zw, zgt10, zgt18 394 393 !!------------------------------------------------------------------- 395 !396 394 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 397 !398 zw = pwnd(ji,jj) ! wind speed399 !400 ! Charnock's constant, increases with the wind :401 zgt10 = 0.5 + SIGN(0.5_wp,(zw - 10)) ! If zw<10. --> 0, else --> 1402 zgt18 = 0.5 + SIGN(0.5_wp,(zw - 18.)) ! If zw<18. --> 0, else --> 1403 !404 alfa_charn_3p0(ji,jj) = (1. - zgt10)*0.011 & ! wind is lower than 10 m/s405 & + zgt10*((1. - zgt18)*(0.011 + (0.018 - 0.011) &406 & *(zw - 10.)/(18. - 10.)) + zgt18*( 0.018 ) ) ! Hare et al. (1999)407 !395 ! 396 zw = pwnd(ji,jj) ! wind speed 397 ! 398 ! Charnock's constant, increases with the wind : 399 zgt10 = 0.5 + SIGN(0.5_wp,(zw - 10)) ! If zw<10. --> 0, else --> 1 400 zgt18 = 0.5 + SIGN(0.5_wp,(zw - 18.)) ! If zw<18. --> 0, else --> 1 401 ! 402 charn_coare3p0(ji,jj) = (1. - zgt10)*0.011 & ! wind is lower than 10 m/s 403 & + zgt10*((1. - zgt18)*(0.011 + (0.018 - 0.011) & 404 & *(zw - 10.)/(18. - 10.)) + zgt18*( 0.018 ) ) ! Hare et al. (1999) 405 ! 408 406 END_2D 409 ! 410 END FUNCTION alfa_charn_3p0 407 END FUNCTION charn_coare3p0 411 408 412 409 FUNCTION psi_m_coare( pzeta ) … … 429 426 REAL(wp) :: zta, zphi_m, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab 430 427 !!---------------------------------------------------------------------------------- 431 !432 428 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 ! 429 ! 430 zta = pzeta(ji,jj) 431 ! 432 zphi_m = ABS(1. - 15.*zta)**.25 !!Kansas unstable 433 ! 434 zpsi_k = 2.*LOG((1. + zphi_m)/2.) + LOG((1. + zphi_m*zphi_m)/2.) & 435 & - 2.*ATAN(zphi_m) + 0.5*rpi 436 ! 437 zphi_c = ABS(1. - 10.15*zta)**.3333 !!Convective 438 ! 439 zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) & 440 & - 1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447 441 ! 442 zf = zta*zta 443 zf = zf/(1. + zf) 444 zc = MIN(50._wp, 0.35_wp*zta) 445 zstab = 0.5 + SIGN(0.5_wp, zta) 446 ! 447 psi_m_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & ! (zta < 0) 448 & - zstab * ( 1. + 1.*zta & ! (zta > 0) 449 & + 0.6667*(zta - 14.28)/EXP(zc) + 8.525 ) ! " 455 450 END_2D 456 !457 451 END FUNCTION psi_m_coare 458 452 … … 474 468 !! (https://github.com/brodeau/aerobulk/) 475 469 !!---------------------------------------------------------------- 476 !!477 470 REAL(wp), DIMENSION(jpi,jpj) :: psi_h_coare 478 471 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta … … 480 473 INTEGER :: ji, jj ! dummy loop indices 481 474 REAL(wp) :: zta, zphi_h, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab 482 ! 475 !!---------------------------------------------------------------- 483 476 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 ! 477 ! 478 zta = pzeta(ji,jj) 479 ! 480 zphi_h = (ABS(1. - 15.*zta))**.5 !! Kansas unstable (zphi_h = zphi_m**2 when unstable, zphi_m when stable) 481 ! 482 zpsi_k = 2.*LOG((1. + zphi_h)/2.) 483 ! 484 zphi_c = (ABS(1. - 34.15*zta))**.3333 !! Convective 485 ! 486 zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) & 487 & -1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447 488 ! 489 zf = zta*zta 490 zf = zf/(1. + zf) 491 zc = MIN(50._wp,0.35_wp*zta) 492 zstab = 0.5 + SIGN(0.5_wp, zta) 493 ! 494 psi_h_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & 495 & - zstab * ( (ABS(1. + 2.*zta/3.))**1.5 & 496 & + .6667*(zta - 14.28)/EXP(zc) + 8.525 ) 505 497 END_2D 506 !507 498 END FUNCTION psi_h_coare 508 499
Note: See TracChangeset
for help on using the changeset viewer.
