Changeset 11111 for NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare.F90
- Timestamp:
- 2019-06-14T15:55:28+02:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare.F90
r10069 r11111 1 1 MODULE sbcblk_algo_coare 2 2 !!====================================================================== 3 !! *** MODULE sbcblk_algo_coare *** 4 !! Computes turbulent components of surface fluxes 5 !! according to Fairall et al. 2003 (COARE v3) 6 !! 3 !! *** MODULE sbcblk_algo_coare *** 4 !! Computes: 7 5 !! * bulk transfer coefficients C_D, C_E and C_H 8 6 !! * air temp. and spec. hum. adjusted from zt (2m) to zu (10m) if needed … … 10 8 !! => all these are used in bulk formulas in sbcblk.F90 11 9 !! 12 !! Using the bulk formulation/param. of COARE v3, Fairall et al. 2003 13 !! 10 !! Using the bulk formulation/param. of COARE v3, Fairall et al., 2003 14 11 !! 15 12 !! Routine turb_coare maintained and developed in AeroBulk 16 !! (http ://aerobulk.sourceforge.net/)13 !! (https://github.com/brodeau/aerobulk) 17 14 !! 18 !! Author: Laurent Brodeau, 2016, brodeau@gmail.com 19 !! 20 !!====================================================================== 21 !! History : 3.6 ! 2016-02 (L.Brodeau) Original code 15 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk) 16 !!---------------------------------------------------------------------- 17 !! History : 4.0 ! 2016-02 (L.Brodeau) Original code 22 18 !!---------------------------------------------------------------------- 23 19 … … 27 23 !! returns the effective bulk wind speed at 10m 28 24 !!---------------------------------------------------------------------- 29 USE oce ! ocean dynamics and tracers30 USE dom_oce ! ocean space and time domain31 USE phycst ! physical constants32 USE sbc_oce ! Surface boundary condition: ocean fields33 USE sbcwave, ONLY : cdn_wave ! wave module25 USE oce ! ocean dynamics and tracers 26 USE dom_oce ! ocean space and time domain 27 USE phycst ! physical constants 28 USE sbc_oce ! Surface boundary condition: ocean fields 29 USE sbcwave, ONLY : cdn_wave ! wave module 34 30 #if defined key_si3 || defined key_cice 35 USE sbc_ice ! Surface boundary condition: ice fields31 USE sbc_ice ! Surface boundary condition: ice fields 36 32 #endif 37 !33 USE sbcblk_phymbl !LB: all thermodynamics functions, rho_air, q_sat, etc... 38 34 USE in_out_manager ! I/O manager 39 35 USE iom ! I/O manager library … … 50 46 REAL(wp), PARAMETER :: zi0 = 600._wp ! scale height of the atmospheric boundary layer... 51 47 REAL(wp), PARAMETER :: Beta0 = 1.250_wp ! gustiness parameter 52 REAL(wp), PARAMETER :: rctv0 = 0.608_wp ! constant to obtain virtual temperature... 48 49 INTEGER , PARAMETER :: nb_itt = 5 ! number of itterations 53 50 54 51 !!---------------------------------------------------------------------- … … 61 58 !! *** ROUTINE turb_coare *** 62 59 !! 63 !! 2015: L. Brodeau (brodeau@gmail.com)64 !!65 60 !! ** Purpose : Computes turbulent transfert coefficients of surface 66 61 !! fluxes according to Fairall et al. (2003) 67 62 !! If relevant (zt /= zu), adjust temperature and humidity from height zt to zu 68 !! 69 !! ** Method : Monin Obukhov Similarity Theory 70 !!---------------------------------------------------------------------- 63 !! Returns the effective bulk wind speed at 10m to be used in the bulk formulas 64 !! 71 65 !! 72 66 !! INPUT : … … 88 82 !! * t_zu : pot. air temperature adjusted at wind height zu [K] 89 83 !! * q_zu : specific humidity of air // [kg/kg] 90 !! * U_blk : bulk wind at 10m [m/s] 91 !!---------------------------------------------------------------------- 84 !! * U_blk : bulk wind speed at 10m [m/s] 85 !! 86 !! 87 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk) 88 !!---------------------------------------------------------------------------------- 92 89 REAL(wp), INTENT(in ) :: zt ! height for t_zt and q_zt [m] 93 90 REAL(wp), INTENT(in ) :: zu ! height for U_zu [m] … … 106 103 ! 107 104 INTEGER :: j_itt 108 LOGICAL :: l_zt_equal_zu = .FALSE. ! if q and t are given at same height as U 109 INTEGER , PARAMETER :: nb_itt = 4 ! number of itterations 105 LOGICAL :: l_zt_equal_zu = .FALSE. ! if q and t are given at same height as U 110 106 111 107 REAL(wp), DIMENSION(jpi,jpj) :: & … … 125 121 126 122 !! First guess of temperature and humidity at height zu: 127 t_zu = MAX( t_zt , 0.0)! who knows what's given on masked-continental regions...128 q_zu = MAX( q_zt , 1.e-6)! "123 t_zu = MAX( t_zt , 199.0_wp ) ! who knows what's given on masked-continental regions... 124 q_zu = MAX( q_zt , 1.e-6_wp ) ! " 129 125 130 126 !! Pot. temp. difference (and we don't want it to be 0!) 131 dt_zu = t_zu - sst ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6 ), dt_zu )132 dq_zu = q_zu - ssq ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9 ), dq_zu )133 134 znu_a = visc_air(t_z t) ! Air viscosity (m^2/s) at zt given from temperature in (K)127 dt_zu = t_zu - sst ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 128 dq_zu = q_zu - ssq ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 129 130 znu_a = visc_air(t_zu) ! Air viscosity (m^2/s) at zt given from temperature in (K) 135 131 136 132 ztmp2 = 0.5*0.5 ! initial guess for wind gustiness contribution … … 144 140 145 141 z0 = alfa_charn(U_blk)*u_star*u_star/grav + 0.11*znu_a/u_star 146 z0t = 0.1*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) ! WARNING: 1/z0t ! 142 z0 = MIN(ABS(z0), 0.001) ! (prevent FPE from stupid values from masked region later on...) !#LOLO 143 z0t = 1. / ( 0.1*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) ) 144 z0t = MIN(ABS(z0t), 0.001) ! (prevent FPE from stupid values from masked region later on...) !#LOLO 147 145 148 146 ztmp2 = vkarmn/ztmp0 149 147 Cd = ztmp2*ztmp2 ! first guess of Cd 150 148 151 ztmp0 = vkarmn*vkarmn/LOG(zt*z0t)/Cd 152 153 !Ribcu = -zu/(zi0*0.004*Beta0**3) !! Saturation Rib, zi0 = tropicalbound. layer depth 154 ztmp2 = grav*zu*(dt_zu + rctv0*t_zu*dq_zu)/(t_zu*U_blk*U_blk) !! Ribu Bulk Richardson number 155 ztmp1 = 0.5 + sign(0.5 , ztmp2) 149 ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd 150 151 ztmp2 = grav*zu*(dt_zu + rctv0*t_zu*dq_zu)/(t_zu*U_blk*U_blk) !! Ribu Bulk Richardson number ; !Ribcu = -zu/(zi0*0.004*Beta0**3) !! Saturation Rib, zi0 = tropicalbound. layer depth 152 153 !! First estimate of zeta_u, depending on the stability, ie sign of Ribu (ztmp2): 154 ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 ) 156 155 ztmp0 = ztmp0*ztmp2 157 !! Ribu < 0 Ribu > 0 Beta = 1.25158 zeta_u = (1.-ztmp1) * (ztmp0/(1.+ztmp2/(-zu/(zi0*0.004*Beta0**3)))) &159 & + ztmp1 * (ztmp0*(1. + 27./9.*ztmp2/ztmp0))156 zeta_u = (1.-ztmp1) * (ztmp0/(1.+ztmp2/(-zu/(zi0*0.004*Beta0**3)))) & ! Ribu < 0 157 & + ztmp1 * (ztmp0*(1. + 27./9.*ztmp2/ztmp0)) ! Ribu > 0 158 !#LOLO: should make sure that the "ztmp0" of "27./9.*ztmp2/ztmp0" is "ztmp0*ztmp2" and not "ztmp0==vkarmn*vkarmn/LOG(zt/z0t)/Cd" ! 160 159 161 160 !! First guess M-O stability dependent scaling params.(u*,t*,q*) to estimate z0 and z/L 162 ztmp0 = vkarmn/(LOG(zu*z0t) - psi_h_coare(zeta_u))161 ztmp0 = vkarmn/(LOG(zu/z0t) - psi_h_coare(zeta_u)) 163 162 164 163 u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0) - psi_m_coare(zeta_u)) … … 168 167 ! What's need to be done if zt /= zu: 169 168 IF( .NOT. l_zt_equal_zu ) THEN 170 169 !! First update of values at zu (or zt for wind) 171 170 zeta_t = zt*zeta_u/zu 172 173 !! First update of values at zu (or zt for wind)174 171 ztmp0 = psi_h_coare(zeta_u) - psi_h_coare(zeta_t) 175 ztmp1 = log(zt/zu) + ztmp0172 ztmp1 = LOG(zt/zu) + ztmp0 176 173 t_zu = t_zt - t_star/vkarmn*ztmp1 177 174 q_zu = q_zt - q_star/vkarmn*ztmp1 178 q_zu = (0.5 + sign(0.5,q_zu))*q_zu !Makes it impossible to have negative humidity : 179 180 dt_zu = t_zu - sst ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6), dt_zu ) 181 dq_zu = q_zu - ssq ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9), dq_zu ) 182 175 q_zu = (0.5 + SIGN(0.5_wp,q_zu))*q_zu !Makes it impossible to have negative humidity : 176 ! 177 dt_zu = t_zu - sst ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 178 dq_zu = q_zu - ssq ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 183 179 END IF 184 180 … … 188 184 !!Inverse of Monin-Obukov length (1/L) : 189 185 ztmp0 = One_on_L(t_zu, q_zu, u_star, t_star, q_star) ! 1/L == 1/[Monin-Obukhov length] 186 ztmp0 = SIGN( MIN(ABS(ztmp0),200._wp), ztmp0 ) ! (prevents FPE from stupid values from masked region later on...) !#LOLO 190 187 191 188 ztmp1 = u_star*u_star ! u*^2 … … 193 190 !! Update wind at 10m taking into acount convection-related wind gustiness: 194 191 ! Ug = Beta*w* (Beta = 1.25, Fairall et al. 2003, Eq.8): 195 ztmp2 = Beta0*Beta0*ztmp1*(MAX(-zi0*ztmp0/vkarmn,0. ))**(2./3.) ! => ztmp2 == Ug^2192 ztmp2 = Beta0*Beta0*ztmp1*(MAX(-zi0*ztmp0/vkarmn,0._wp))**(2./3.) ! => ztmp2 == Ug^2 196 193 !! ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before 600. 197 U_blk = MAX(sqrt(U_zu*U_zu + ztmp2), 0.2 ) ! include gustiness in bulk wind speed194 U_blk = MAX(sqrt(U_zu*U_zu + ztmp2), 0.2_wp) ! include gustiness in bulk wind speed 198 195 ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0. 199 196 … … 203 200 !! Roughness lengthes z0, z0t (z0q = z0t) : 204 201 z0 = ztmp2*ztmp1/grav + 0.11*znu_a/u_star ! Roughness length (eq.6) 205 ztmp1 = z0*u_star/znu_a 206 z0t = min( 1.1E-4 , 5.5E-5*ztmp1**(-0.6) )! Scalar roughness for both theta and q (eq.28)202 ztmp1 = z0*u_star/znu_a ! Re_r: roughness Reynolds number 203 z0t = MIN( 1.1E-4_wp , 5.5E-5_wp*ztmp1**(-0.6_wp) ) ! Scalar roughness for both theta and q (eq.28) 207 204 208 205 !! Stability parameters: 209 zeta_u = zu*ztmp0 ; zeta_u = sign( min(abs(zeta_u),50.0), zeta_u ) 206 zeta_u = zu*ztmp0 207 zeta_u = SIGN( MIN(ABS(zeta_u),50.0_wp), zeta_u ) 210 208 IF( .NOT. l_zt_equal_zu ) THEN 211 zeta_t = zt*ztmp0 ; zeta_t = sign( min(abs(zeta_t),50.0), zeta_t ) 209 zeta_t = zt*ztmp0 210 zeta_t = SIGN( MIN(ABS(zeta_t),50.0_wp), zeta_t ) 212 211 END IF 213 212 214 213 !! Turbulent scales at zu=10m : 215 214 ztmp0 = psi_h_coare(zeta_u) 216 ztmp1 = vkarmn/(LOG(zu) - LOG(z0t) - ztmp0)215 ztmp1 = vkarmn/(LOG(zu) - LOG(z0t) - ztmp0) 217 216 218 217 t_star = dt_zu*ztmp1 219 218 q_star = dq_zu*ztmp1 220 u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0) - psi_m_coare(zeta_u))219 u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0) - psi_m_coare(zeta_u)) 221 220 222 221 IF( .NOT. l_zt_equal_zu ) THEN … … 259 258 !! Wind greater than 18 m/s : alfa = 0.018 260 259 !! 261 !! Author: L. Brodeau, june 2016 / AeroBulk (https:// sourceforge.net/p/aerobulk)260 !! Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 262 261 !!------------------------------------------------------------------- 263 262 REAL(wp), DIMENSION(jpi,jpj) :: alfa_charn … … 274 273 ! 275 274 ! Charnock's constant, increases with the wind : 276 zgt10 = 0.5 + SIGN(0.5 ,(zw - 10.)) ! If zw<10. --> 0, else --> 1277 zgt18 = 0.5 + SIGN(0.5 ,(zw - 18.)) ! If zw<18. --> 0, else --> 1275 zgt10 = 0.5 + SIGN(0.5_wp,(zw - 10)) ! If zw<10. --> 0, else --> 1 276 zgt18 = 0.5 + SIGN(0.5_wp,(zw - 18.)) ! If zw<18. --> 0, else --> 1 278 277 ! 279 278 alfa_charn(ji,jj) = (1. - zgt10)*0.011 & ! wind is lower than 10 m/s … … 294 293 !! 295 294 !! Author: L. Brodeau, june 2016 / AeroBulk 296 !! (https:// sourceforge.net/p/aerobulk)295 !! (https://github.com/brodeau/aerobulk/) 297 296 !!------------------------------------------------------------------------ 298 297 REAL(wp), DIMENSION(jpi,jpj) :: One_on_L !: 1./(Monin Obukhov length) [m^-1] … … 330 329 !! form from Beljaars and Holtslag (1991) 331 330 !! 332 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https:// sourceforge.net/p/aerobulk)331 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 333 332 !!---------------------------------------------------------------------------------- 334 333 REAL(wp), DIMENSION(jpi,jpj) :: psi_m_coare … … 356 355 zf = zta*zta 357 356 zf = zf/(1. + zf) 358 zc = MIN(50. , 0.35*zta)359 zstab = 0.5 + SIGN(0.5 , zta)357 zc = MIN(50._wp, 0.35_wp*zta) 358 zstab = 0.5 + SIGN(0.5_wp, zta) 360 359 ! 361 360 psi_m_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & ! (zta < 0) … … 383 382 !! 384 383 !! Author: L. Brodeau, june 2016 / AeroBulk 385 !! (https:// sourceforge.net/p/aerobulk)384 !! (https://github.com/brodeau/aerobulk/) 386 385 !!---------------------------------------------------------------- 387 386 !! … … 408 407 zf = zta*zta 409 408 zf = zf/(1. + zf) 410 zc = MIN(50. ,0.35*zta)411 zstab = 0.5 + SIGN(0.5 , zta)409 zc = MIN(50._wp,0.35_wp*zta) 410 zstab = 0.5 + SIGN(0.5_wp, zta) 412 411 ! 413 412 psi_h_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & … … 420 419 END FUNCTION psi_h_coare 421 420 422 423 FUNCTION visc_air( ptak ) 424 !!--------------------------------------------------------------------- 425 !! Air kinetic viscosity (m^2/s) given from temperature in degrees... 426 !! 427 !! Author: L. Brodeau, june 2016 / AeroBulk 428 !! (https://sourceforge.net/p/aerobulk) 429 !!--------------------------------------------------------------------- 430 !! 431 REAL(wp), DIMENSION(jpi,jpj) :: visc_air 432 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature in (K) 433 ! 434 INTEGER :: ji, jj ! dummy loop indices 435 REAL(wp) :: ztc, ztc2 ! local scalar 436 ! 437 DO jj = 1, jpj 438 DO ji = 1, jpi 439 ztc = ptak(ji,jj) - rt0 ! air temp, in deg. C 440 ztc2 = ztc*ztc 441 visc_air(ji,jj) = 1.326E-5*(1. + 6.542E-3*ztc + 8.301E-6*ztc2 - 4.84E-9*ztc2*ztc) 442 END DO 443 END DO 444 ! 445 END FUNCTION visc_air 446 447 421 !!====================================================================== 448 422 END MODULE sbcblk_algo_coare
Note: See TracChangeset
for help on using the changeset viewer.
