Changeset 11215 for NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare.F90
- Timestamp:
- 2019-07-04T12:34: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
r11209 r11215 13 13 !! (https://github.com/brodeau/aerobulk) 14 14 !! 15 !! ** Author: L. Brodeau, june 2016/ AeroBulk (https://github.com/brodeau/aerobulk)15 !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk) 16 16 !!---------------------------------------------------------------------- 17 17 !! History : 4.0 ! 2016-02 (L.Brodeau) Original code … … 26 26 USE dom_oce ! ocean space and time domain 27 27 USE phycst ! physical constants 28 USE sbc_oce ! Surface boundary condition: ocean fields 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 29 32 USE sbcwave, ONLY : cdn_wave ! wave module 30 33 #if defined key_si3 || defined key_cice 31 34 USE sbc_ice ! Surface boundary condition: ice fields 32 35 #endif 33 USE sbcblk_phy !LB: all thermodynamics functions, rho_air, q_sat, etc... #LB 34 USE in_out_manager ! I/O manager 35 USE iom ! I/O manager library 36 USE lib_mpp ! distribued memory computing library 37 USE prtctl ! Print control 38 USE lib_fortran ! to use key_nosignedzero 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 40 USE sbcblk_skin ! cool-skin/warm layer scheme (CSWL_ECMWF) !LB 39 41 40 42 IMPLICIT NONE … … 52 54 CONTAINS 53 55 54 SUBROUTINE turb_coare( zt, zu, sst, t_zt, ssq, q_zt, U_zu, &56 SUBROUTINE turb_coare( zt, zu, T_s, t_zt, q_s, q_zt, U_zu, & 55 57 & Cd, Ch, Ce, t_zu, q_zu, U_blk, & 56 & Cdn, Chn, Cen ) 58 & Cdn, Chn, Cen, & 59 & Qsw, rad_lw, slp ) 57 60 !!---------------------------------------------------------------------- 58 61 !! *** ROUTINE turb_coare *** … … 63 66 !! Returns the effective bulk wind speed at 10m to be used in the bulk formulas 64 67 !! 68 !! Applies the cool-skin warm-layer correction of the SST to T_s 69 !! if the net shortwave flux at the surface (Qsw), the downwelling longwave 70 !! radiative fluxes at the surface (rad_lw), and the sea-leve pressure (slp) 71 !! are provided as (optional) arguments! 65 72 !! 66 73 !! INPUT : … … 69 76 !! * zu : height for wind speed (generally 10m) [m] 70 77 !! * U_zu : scalar wind speed at 10m [m/s] 71 !! * sst : SST [K]72 78 !! * t_zt : potential air temperature at zt [K] 73 !! * ssq : specific humidity at saturation at SST [kg/kg]74 79 !! * q_zt : specific humidity of air at zt [kg/kg] 75 80 !! 81 !! INPUT/OUTPUT: 82 !! ------------- 83 !! * T_s : SST or skin temperature [K] 84 !! * q_s : SSQ aka saturation specific humidity at temp. T_s [kg/kg] 85 !! -> doesn't need to be given a value if skin temp computed (in case l_use_skin=True) 86 !! -> MUST be given the correct value if not computing skint temp. (in case l_use_skin=False) 87 !! 88 !! OPTIONAL INPUT (will trigger l_use_skin=TRUE if present!): 89 !! --------------- 90 !! * Qsw : net solar flux (after albedo) at the surface (>0) [W/m^2] 91 !! * rad_lw : downwelling longwave radiation at the surface (>0) [W/m^2] 92 !! * slp : sea-level pressure [Pa] 76 93 !! 77 94 !! OUTPUT : … … 85 102 !! 86 103 !! 87 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk)104 !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/) 88 105 !!---------------------------------------------------------------------------------- 89 106 REAL(wp), INTENT(in ) :: zt ! height for t_zt and q_zt [m] 90 107 REAL(wp), INTENT(in ) :: zu ! height for U_zu [m] 91 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: sst! sea surface temperature [Kelvin]108 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) :: T_s ! sea surface temperature [Kelvin] 92 109 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: t_zt ! potential air temperature [Kelvin] 93 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: ssq! sea surface specific humidity [kg/kg]110 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) :: q_s ! sea surface specific humidity [kg/kg] 94 111 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: q_zt ! specific air humidity at zt [kg/kg] 95 112 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: U_zu ! relative wind module at zu [m/s] … … 102 119 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Cdn, Chn, Cen ! neutral transfer coefficients 103 120 ! 121 REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: Qsw ! [W/m^2] 122 REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: rad_lw ! [W/m^2] 123 REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: slp ! [Pa] 124 ! 104 125 INTEGER :: j_itt 105 126 LOGICAL :: l_zt_equal_zu = .FALSE. ! if q and t are given at same height as U … … 113 134 REAL(wp), DIMENSION(jpi,jpj) :: ztmp0, ztmp1, ztmp2 114 135 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zeta_t ! stability parameter at height zt 115 !!---------------------------------------------------------------------- 116 ! 136 ! 137 ! Cool skin: 138 LOGICAL :: l_use_skin = .FALSE. 139 REAL(wp), DIMENSION(jpi,jpj) :: zsst ! to back up the initial bulk SST 140 !!---------------------------------------------------------------------------------- 141 ! 142 ! Cool skin ? 143 IF( PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp) ) THEN 144 l_use_skin = .TRUE. 145 END IF 146 IF (lwp) PRINT *, ' *** LOLO(sbcblk_algo_ecmwf.F90) => l_use_skin =', l_use_skin 147 117 148 l_zt_equal_zu = .FALSE. 118 IF( ABS(zu - zt) < 0.01 )l_zt_equal_zu = .TRUE. ! testing "zu == zt" is risky with double precision149 IF( ABS(zu - zt) < 0.01_wp ) l_zt_equal_zu = .TRUE. ! testing "zu == zt" is risky with double precision 119 150 120 151 IF( .NOT. l_zt_equal_zu ) ALLOCATE( zeta_t(jpi,jpj) ) 121 152 153 !! Initialization for cool skin: 154 IF( l_use_skin ) THEN 155 zsst = T_s ! save the bulk SST 156 T_s = T_s - 0.25 ! First guess of correction 157 q_s = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s 158 END IF 159 122 160 !! First guess of temperature and humidity at height zu: 123 t_zu = MAX( t_zt , 199.0_wp ) ! who knows what's given on masked-continental regions...161 t_zu = MAX( t_zt , 180._wp ) ! who knows what's given on masked-continental regions... 124 162 q_zu = MAX( q_zt , 1.e-6_wp ) ! " 125 163 126 164 !! Pot. temp. difference (and we don't want it to be 0!) 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 )165 dt_zu = t_zu - T_s ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 166 dq_zu = q_zu - q_s ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 129 167 130 168 znu_a = visc_air(t_zu) ! Air viscosity (m^2/s) at zt given from temperature in (K) 131 169 132 ztmp2 = 0.5 *0.5! initial guess for wind gustiness contribution170 ztmp2 = 0.5_wp*0.5_wp ! initial guess for wind gustiness contribution 133 171 U_blk = SQRT(U_zu*U_zu + ztmp2) 134 172 135 ztmp2 = 10000. ! optimization: ztmp2 == 1/z0 (with z0 first guess == 0.0001)173 ztmp2 = 10000._wp ! optimization: ztmp2 == 1/z0 (with z0 first guess == 0.0001) 136 174 ztmp0 = LOG(zu*ztmp2) 137 175 ztmp1 = LOG(10.*ztmp2) 138 u_star = 0.035*U_blk*ztmp1/ztmp0 ! (u* = 0.035*Un10) 139 140 141 z0 = alfa_charn(U_blk)*u_star*u_star/grav + 0.11*znu_a/u_star 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 176 u_star = 0.035_wp*U_blk*ztmp1/ztmp0 ! (u* = 0.035*Un10) 177 178 z0 = alfa_charn(U_blk)*u_star*u_star/grav + 0.11_wp*znu_a/u_star 179 z0 = MIN(ABS(z0), 0.001_wp) ! (prevent FPE from stupid values from masked region later on...) !#LOLO 180 z0t = 1._wp / ( 0.1_wp*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) ) 181 z0t = MIN(ABS(z0t), 0.001_wp) ! (prevent FPE from stupid values from masked region later on...) !#LOLO 145 182 146 183 ztmp2 = vkarmn/ztmp0 … … 149 186 ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd 150 187 151 ztmp2 = Ri_bulk( zu, sst, t_zu, ssq, q_zu, U_blk ) ! Bulk Richardson Number (BRN) 152 !ztmp2 = grav*zu*(dt_zu + rctv0*t_zu*dq_zu)/(t_zu*U_blk*U_blk) !! Bulk Richardson number ; !Ribcu = -zu/(zi0*0.004*Beta0**3) !! Saturation Rib, zi0 = tropicalbound. layer depth 188 ztmp2 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN) 153 189 154 190 !! First estimate of zeta_u, depending on the stability, ie sign of BRN (ztmp2): 155 191 ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 ) 156 192 ztmp0 = ztmp0*ztmp2 157 zeta_u = (1. -ztmp1) * (ztmp0/(1.+ztmp2/(-zu/(zi0*0.004*Beta0**3)))) & ! BRN < 0158 & + ztmp1 * (ztmp0*(1. + 27./9.*ztmp2/ztmp0)) ! BRN > 0159 !#L OLO: should make sure that the "ztmp0" of "27./9.*ztmp2/ztmp0" is "ztmp0*ztmp2" and not "ztmp0==vkarmn*vkarmn/LOG(zt/z0t)/Cd" !193 zeta_u = (1._wp-ztmp1) * (ztmp0/(1._wp+ztmp2/(-zu/(zi0*0.004_wp*Beta0**3)))) & ! BRN < 0 194 & + ztmp1 * (ztmp0*(1._wp + 27._wp/9._wp*ztmp2/ztmp0)) ! BRN > 0 195 !#LB: should make sure that the "ztmp0" of "27./9.*ztmp2/ztmp0" is "ztmp0*ztmp2" and not "ztmp0==vkarmn*vkarmn/LOG(zt/z0t)/Cd" ! 160 196 161 197 !! First guess M-O stability dependent scaling params.(u*,t*,q*) to estimate z0 and z/L … … 174 210 t_zu = t_zt - t_star/vkarmn*ztmp1 175 211 q_zu = q_zt - q_star/vkarmn*ztmp1 176 q_zu = (0.5 + SIGN(0.5_wp,q_zu))*q_zu !Makes it impossible to have negative humidity :212 q_zu = (0.5_wp + SIGN(0.5_wp,q_zu))*q_zu !Makes it impossible to have negative humidity : 177 213 ! 178 dt_zu = t_zu - sst; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )179 dq_zu = q_zu - ssq; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )214 dt_zu = t_zu - T_s ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 215 dq_zu = q_zu - q_s ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 180 216 END IF 181 217 … … 200 236 201 237 !! Roughness lengthes z0, z0t (z0q = z0t) : 202 z0 = ztmp2*ztmp1/grav + 0.11*znu_a/u_star! Roughness length (eq.6)203 ztmp1 = z0*u_star/znu_a 238 z0 = ztmp2*ztmp1/grav + 0.11_wp*znu_a/u_star ! Roughness length (eq.6) 239 ztmp1 = z0*u_star/znu_a ! Re_r: roughness Reynolds number 204 240 z0t = MIN( 1.1E-4_wp , 5.5E-5_wp*ztmp1**(-0.6_wp) ) ! Scalar roughness for both theta and q (eq.28) 205 241 … … 227 263 t_zu = t_zt - t_star/vkarmn*ztmp1 228 264 q_zu = q_zt - q_star/vkarmn*ztmp1 229 dt_zu = t_zu - sst ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6), dt_zu )230 dq_zu = q_zu - ssq ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9), dq_zu )231 265 END IF 232 266 233 END DO 234 ! 267 !! SKIN related part 268 !! ----------------- 269 IF( l_use_skin ) THEN 270 !! compute transfer coefficients at zu : lolo: verifier... 271 ztmp0 = u_star/U_blk 272 Ch = ztmp0*t_star/dt_zu 273 Ce = ztmp0*q_star/dq_zu 274 ! Non-Solar heat flux to the ocean: 275 ztmp1 = U_blk*MAX(rho_air(t_zu, q_zu, slp), 1._wp) ! rho*U10 276 ztmp2 = T_s*T_s 277 ztmp1 = ztmp1 * ( Ce*rLevap*(q_zu - q_s) + Ch*rCp_dry*(t_zu - T_s) ) & ! Total turb. heat flux 278 & + (rad_lw - emiss_w*stefan*ztmp2*ztmp2) ! Net longwave flux 279 !! Updating the values of the skin temperature T_s and q_s : 280 CALL CSWL_ECMWF( Qsw, ztmp1, u_star, zsst, T_s ) ! yes ECMWF, because more advanced than COARE (warm-layer added!) 281 q_s = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! 200 -> just to avoid numerics problem on masked regions if silly values are given 282 END IF 283 284 IF( (l_use_skin).OR.(.NOT. l_zt_equal_zu) ) THEN 285 dt_zu = t_zu - T_s ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 286 dq_zu = q_zu - q_s ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 287 END IF 288 289 END DO !DO j_itt = 1, nb_itt 290 235 291 ! compute transfer coefficients at zu : 236 292 ztmp0 = u_star/U_blk … … 259 315 !! Wind greater than 18 m/s : alfa = 0.018 260 316 !! 261 !! Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)317 !! Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 262 318 !!------------------------------------------------------------------- 263 319 REAL(wp), DIMENSION(jpi,jpj) :: alfa_charn … … 297 353 !! form from Beljaars and Holtslag (1991) 298 354 !! 299 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)355 !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 300 356 !!---------------------------------------------------------------------------------- 301 357 REAL(wp), DIMENSION(jpi,jpj) :: psi_m_coare … … 349 405 !! form from Beljaars and Holtslag (1991) 350 406 !! 351 !! Author: L. Brodeau, june 2016 / AeroBulk407 !! Author: L. Brodeau, June 2016 / AeroBulk 352 408 !! (https://github.com/brodeau/aerobulk/) 353 409 !!----------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.