Changeset 11209 for NEMO/branches/2019
- Timestamp:
- 2019-07-03T09:09:17+02:00 (5 years ago)
- Location:
- NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC
- Files:
-
- 5 edited
-
sbcblk_algo_coare.F90 (modified) (3 diffs)
-
sbcblk_algo_coare3p5.F90 (modified) (3 diffs)
-
sbcblk_algo_ecmwf.F90 (modified) (7 diffs)
-
sbcblk_algo_ncar.F90 (modified) (10 diffs)
-
sbcblk_phy.F90 (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare.F90
r11182 r11209 31 31 USE sbc_ice ! Surface boundary condition: ice fields 32 32 #endif 33 USE sbcblk_phy !LB: all thermodynamics functions, rho_air, q_sat, etc... 33 USE sbcblk_phy !LB: all thermodynamics functions, rho_air, q_sat, etc... #LB 34 34 USE in_out_manager ! I/O manager 35 35 USE iom ! I/O manager library … … 149 149 ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd 150 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): 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 153 154 !! First estimate of zeta_u, depending on the stability, ie sign of BRN (ztmp2): 154 155 ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 ) 155 156 ztmp0 = ztmp0*ztmp2 156 zeta_u = (1.-ztmp1) * (ztmp0/(1.+ztmp2/(-zu/(zi0*0.004*Beta0**3)))) & ! Ribu< 0157 & + ztmp1 * (ztmp0*(1. + 27./9.*ztmp2/ztmp0)) ! Ribu> 0157 zeta_u = (1.-ztmp1) * (ztmp0/(1.+ztmp2/(-zu/(zi0*0.004*Beta0**3)))) & ! BRN < 0 158 & + ztmp1 * (ztmp0*(1. + 27./9.*ztmp2/ztmp0)) ! BRN > 0 158 159 !#LOLO: should make sure that the "ztmp0" of "27./9.*ztmp2/ztmp0" is "ztmp0*ztmp2" and not "ztmp0==vkarmn*vkarmn/LOG(zt/z0t)/Cd" ! 159 160 … … 285 286 END FUNCTION alfa_charn 286 287 287 288 FUNCTION One_on_L( ptha, pqa, pus, pts, pqs )289 !!------------------------------------------------------------------------290 !!291 !! Evaluates the 1./(Monin Obukhov length) from air temperature and292 !! specific humidity, and frictional scales u*, t* and q*293 !!294 !! Author: L. Brodeau, june 2016 / AeroBulk295 !! (https://github.com/brodeau/aerobulk/)296 !!------------------------------------------------------------------------297 REAL(wp), DIMENSION(jpi,jpj) :: One_on_L !: 1./(Monin Obukhov length) [m^-1]298 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptha, & !: average potetntial air temperature [K]299 & pqa, & !: average specific humidity of air [kg/kg]300 & pus, pts, pqs !: frictional velocity, temperature and humidity301 !302 INTEGER :: ji, jj ! dummy loop indices303 REAL(wp) :: zqa ! local scalar304 !!-------------------------------------------------------------------305 !306 DO jj = 1, jpj307 DO ji = 1, jpi308 !309 zqa = (1. + rctv0*pqa(ji,jj))310 !311 One_on_L(ji,jj) = grav*vkarmn*(pts(ji,jj)*zqa + rctv0*ptha(ji,jj)*pqs(ji,jj)) &312 & / ( pus(ji,jj)*pus(ji,jj) * ptha(ji,jj)*zqa )313 !314 END DO315 END DO316 !317 END FUNCTION One_on_L318 319 320 288 FUNCTION psi_m_coare( pzeta ) 321 289 !!---------------------------------------------------------------------------------- -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p5.F90
r11111 r11209 35 35 USE sbc_ice ! Surface boundary condition: ice fields 36 36 #endif 37 !37 USE sbcblk_phy !LB: all thermodynamics functions, rho_air, q_sat, etc... #LB 38 38 USE iom ! I/O manager library 39 39 USE lib_mpp ! distribued memory computing library … … 255 255 END SUBROUTINE turb_coare3p5 256 256 257 258 259 FUNCTION One_on_L( ptha, pqa, pus, pts, pqs )260 !!------------------------------------------------------------------------261 !!262 !! Evaluates the 1./(Monin Obukhov length) from air temperature and263 !! specific humidity, and frictional scales u*, t* and q*264 !!265 !! Author: L. Brodeau, june 2016 / AeroBulk266 !! (https://github.com/brodeau/aerobulk/)267 !!------------------------------------------------------------------------268 REAL(wp), DIMENSION(jpi,jpj) :: One_on_L !: 1./(Monin Obukhov length) [m^-1]269 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptha, & !: average potetntial air temperature [K]270 & pqa, & !: average specific humidity of air [kg/kg]271 & pus, pts, pqs !: frictional velocity, temperature and humidity272 !273 INTEGER :: ji, jj ! dummy loop indices274 REAL(wp) :: zqa ! local scalar275 !276 DO jj = 1, jpj277 DO ji = 1, jpi278 !279 zqa = (1. + rctv0*pqa(ji,jj))280 !281 One_on_L(ji,jj) = grav*vkarmn*(pts(ji,jj)*zqa + rctv0*ptha(ji,jj)*pqs(ji,jj)) &282 & / ( pus(ji,jj)*pus(ji,jj) * ptha(ji,jj)*zqa )283 !284 END DO285 END DO286 !287 END FUNCTION One_on_L288 289 290 257 FUNCTION psi_m_coare( pzeta ) 291 258 !!---------------------------------------------------------------------------------- … … 389 356 END FUNCTION psi_h_coare 390 357 391 392 FUNCTION visc_air( ptak )393 !!---------------------------------------------------------------------394 !! Air kinetic viscosity (m^2/s) given from temperature in degrees...395 !!396 !! Author: L. Brodeau, june 2016 / AeroBulk397 !! (https://github.com/brodeau/aerobulk/)398 !!---------------------------------------------------------------------399 REAL(wp), DIMENSION(jpi,jpj) :: visc_air400 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K]401 !402 INTEGER :: ji, jj ! dummy loop indices403 REAL(wp) :: ztc, ztc2 ! local scalar404 !405 DO jj = 1, jpj406 DO ji = 1, jpi407 ztc = ptak(ji,jj) - rt0 ! air temp, in deg. C408 ztc2 = ztc*ztc409 visc_air(ji,jj) = 1.326E-5*(1. + 6.542E-3*ztc + 8.301E-6*ztc2 - 4.84E-9*ztc2*ztc)410 END DO411 END DO412 !413 END FUNCTION visc_air414 415 358 !!====================================================================== 416 359 END MODULE sbcblk_algo_coare3p5 -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ecmwf.F90
r11182 r11209 38 38 39 39 USE sbc_oce ! Surface boundary condition: ocean fields 40 USE sbcblk_phy !LB: all thermodynamics functions, rho_air, q_sat, etc...40 USE sbcblk_phy !LB: all thermodynamics functions, rho_air, q_sat, etc... #LB 41 41 USE sbcblk_skin ! cool-skin/warm layer scheme (CSWL_ECMWF) 42 42 … … 72 72 ! !: => ACCEPTABLE IN MOST CONDITIONS ! (UNLESS: sunny + very calm/low-wind conditions) 73 73 ! !: => Otherwize use "nb_itt_wl = 10" 74 75 76 77 78 74 !!---------------------------------------------------------------------- 79 75 CONTAINS … … 213 209 ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd 214 210 215 ztmp2 = Ri_bulk( zu, t_zu, dt_zu, q_zu, dq_zu, U_blk ) ! Ribu = Bulk Richardson number216 217 !! First estimate of zeta_u, depending on the stability, ie sign of Ribu(ztmp2):211 ztmp2 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN) 212 213 !! First estimate of zeta_u, depending on the stability, ie sign of BRN (ztmp2): 218 214 ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 ) 219 215 func_m = ztmp0*ztmp2 ! temporary array !! 220 func_h = (1.-ztmp1) * (func_m/(1._wp+ztmp2/(-zu/(zi0*0.004_wp*Beta0**3)))) & ! Ribu< 0 ! temporary array !!! func_h == zeta_u221 & + ztmp1 * (func_m*(1._wp + 27._wp/9._wp*ztmp2/func_m)) ! Ribu> 0216 func_h = (1.-ztmp1) * (func_m/(1._wp+ztmp2/(-zu/(zi0*0.004_wp*Beta0**3)))) & ! BRN < 0 ! temporary array !!! func_h == zeta_u 217 & + ztmp1 * (func_m*(1._wp + 27._wp/9._wp*ztmp2/func_m)) ! BRN > 0 222 218 !#LB: should make sure that the "func_m" of "27./9.*ztmp2/func_m" is "ztmp0*ztmp2" and not "ztmp0==vkarmn*vkarmn/LOG(zt/z0t)/Cd" ! 223 219 … … 247 243 248 244 !! First guess of inverse of Monin-Obukov length (1/L) : 249 ztmp0 = (1. + rctv0*q_zu) ! the factor to apply to temp. to get virt. temp... 250 Linv = grav*vkarmn*(t_star*ztmp0 + rctv0*t_zu*q_star) / MAX( u_star*u_star * t_zu*ztmp0 , 1.E-9 ) ! #LOLO 251 Linv = SIGN( MIN(ABS(Linv),200._wp), Linv ) ! (prevent FPE from stupid values from masked region later on...) !#LOLO 245 Linv = One_on_L( t_zu, q_zu, u_star, t_star, q_star ) 252 246 253 247 !! Functions such as u* = U_blk*vkarmn/func_m … … 260 254 261 255 !! Bulk Richardson Number at z=zu (Eq. 3.25) 262 ztmp0 = Ri_bulk( zu, t_zu, dt_zu, q_zu, dq_zu, U_blk)256 ztmp0 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN) 263 257 264 258 !! New estimate of the inverse of the Monin-Obukhon length (Linv == zeta/zu) : … … 278 272 z0t = MIN( ABS( alpha_H*ztmp1 ) , 0.001_wp) ! eq.3.26, Chap.3, p.34, IFS doc - Cy31r1 279 273 z0q = MIN( ABS( alpha_Q*ztmp1 ) , 0.001_wp) 280 274 281 275 !! Update wind at 10m taking into acount convection-related wind gustiness: 282 276 !! => Chap. 3.2, IFS doc - Cy40r1, Eq.3.17 and Eq.3.18 + Eq.3.8 … … 447 441 448 442 449 FUNCTION Ri_bulk( pz, ptz, pdt, pqz, pdq, pub )450 !!----------------------------------------------------------------------------------451 !! Bulk Richardson number (Eq. 3.25 IFS doc)452 !!453 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)454 !!----------------------------------------------------------------------------------455 REAL(wp), DIMENSION(jpi,jpj) :: Ri_bulk !456 !457 REAL(wp) , INTENT(in) :: pz ! height above the sea [m]458 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptz ! air temperature at pz m [K]459 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pdt ! ptz - sst [K]460 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqz ! air temperature at pz m [kg/kg]461 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pdq ! pqz - ssq [kg/kg]462 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pub ! bulk wind speed [m/s]463 !!----------------------------------------------------------------------------------464 !465 Ri_bulk = grav*pz/(pub*pub) &466 & * ( pdt/(ptz - 0.5_wp*(pdt + grav*pz/(rCp_dry + rCp_vap*pqz))) &467 & + rctv0*pdq )468 !469 END FUNCTION Ri_bulk470 443 471 444 -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ncar.F90
r11111 r11209 38 38 USE lib_fortran ! to use key_nosignedzero 39 39 40 USE sbcblk_phy !LB: all thermodynamics functions, rho_air, q_sat, etc... #LB 40 41 41 42 IMPLICIT NONE … … 93 94 !! * t_zu : pot. air temperature adjusted at wind height zu [K] 94 95 !! * q_zu : specific humidity of air // [kg/kg] 95 !! * U_blk : bulk wind at 10m[m/s]96 !! * U_blk : bulk wind speed at 10m [m/s] 96 97 !!---------------------------------------------------------------------------------- 97 98 REAL(wp), INTENT(in ) :: zt ! height for t_zt and q_zt [m] … … 125 126 IF( ABS(zu - zt) < 0.01 ) l_zt_equal_zu = .TRUE. ! testing "zu == zt" is risky with double precision 126 127 127 U_blk = MAX( 0.5 , U_zu ) ! relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s128 U_blk = MAX( 0.5_wp , U_zu ) ! relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s 128 129 129 130 !! First guess of stability: 130 ztmp0 = t_zt*(1. + rctv0*q_zt) - sst*(1. + rctv0*ssq) ! air-sea difference of virtual pot. temp. at zt131 stab = 0.5 + sign(0.5,ztmp0) ! stab = 1 if dTv > 0 => STABLE, 0 if unstable131 ztmp0 = virt_temp(t_zt, q_zt) - virt_temp(sst, ssq) ! air-sea difference of virtual pot. temp. at zt 132 stab = 0.5_wp + sign(0.5_wp,ztmp0) ! stab = 1 if dTv > 0 => STABLE, 0 if unstable 132 133 133 134 !! Neutral coefficients at 10m: … … 159 160 160 161 ! Updating turbulent scales : (L&Y 2004 eq. (7)) 161 ztmp1 = Ch/stab*ztmp1 ! theta* (stab == SQRT(Cd)) 162 ztmp2 = Ce/stab*ztmp2 ! q* (stab == SQRT(Cd)) 163 164 ztmp0 = 1. + rctv0*q_zu ! multiply this with t and you have the virtual temperature 162 ztmp0 = stab*U_blk ! u* (stab == SQRT(Cd)) 163 ztmp1 = Ch/stab*ztmp1 ! theta* (stab == SQRT(Cd)) 164 ztmp2 = Ce/stab*ztmp2 ! q* (stab == SQRT(Cd)) 165 165 166 166 ! Estimate the inverse of Monin-Obukov length (1/L) at height zu: 167 ztmp0 = (grav*vkarmn/(t_zu*ztmp0)*(ztmp1*ztmp0 + rctv0*t_zu*ztmp2)) / (Cd*U_blk*U_blk) 168 ! ( Cd*U_blk*U_blk is U*^2 at zu ) 169 167 ztmp0 = One_on_L( t_zu, q_zu, ztmp0, ztmp1, ztmp2 ) 168 170 169 !! Stability parameters : 171 zeta_u = zu*ztmp0 ; zeta_u = sign( min(abs(zeta_u),10.0), zeta_u ) 170 zeta_u = zu*ztmp0 171 zeta_u = sign( min(abs(zeta_u),10.0_wp), zeta_u ) 172 172 zpsi_h_u = psi_h( zeta_u ) 173 173 … … 175 175 IF( .NOT. l_zt_equal_zu ) THEN 176 176 !! Array 'stab' is free for the moment so using it to store 'zeta_t' 177 stab = zt*ztmp0 ; stab = SIGN( MIN(ABS(stab),10.0), stab ) ! Temporaty array stab == zeta_t !!! 177 stab = zt*ztmp0 178 stab = SIGN( MIN(ABS(stab),10.0_wp), stab ) ! Temporaty array stab == zeta_t !!! 178 179 stab = LOG(zt/zu) + zpsi_h_u - psi_h(stab) ! stab just used as temp array again! 179 180 t_zu = t_zt - ztmp1/vkarmn*stab ! ztmp1 is still theta* L&Y 2004 eq.(9b) 180 181 q_zu = q_zt - ztmp2/vkarmn*stab ! ztmp2 is still q* L&Y 2004 eq.(9c) 181 q_zu = max(0. , q_zu)182 q_zu = max(0._wp, q_zu) 182 183 END IF 183 184 … … 194 195 195 196 ELSE 196 197 198 199 200 ztmp0 = MAX( 0.25 , U_blk/(1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - ztmp2)) ) ! U_n10 (ztmp2 == psi_m(zeta_u))201 202 203 204 205 stab = 0.5 + sign(0.5,zeta_u)! update stability206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 197 ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 eq. 9a)... 198 ! In very rare low-wind conditions, the old way of estimating the 199 ! neutral wind speed at 10m leads to a negative value that causes the code 200 ! to crash. To prevent this a threshold of 0.25m/s is imposed. 201 ztmp0 = MAX( 0.25_wp , U_blk/(1._wp + sqrt_Cd_n10/vkarmn*(LOG(zu/10._wp) - ztmp2)) ) ! U_n10 (ztmp2 == psi_m(zeta_u)) 202 ztmp0 = cd_neutral_10m(ztmp0) ! Cd_n10 203 Cdn(:,:) = ztmp0 204 sqrt_Cd_n10 = sqrt(ztmp0) 205 206 stab = 0.5_wp + sign(0.5_wp,zeta_u) ! update stability 207 Cx_n10 = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1. - stab)) ! L&Y 2004 eq. (6c-6d) (Cx_n10 == Ch_n10) 208 Chn(:,:) = Cx_n10 209 210 !! Update of transfer coefficients: 211 ztmp1 = 1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - ztmp2) ! L&Y 2004 eq. (10a) (ztmp2 == psi_m(zeta_u)) 212 Cd = ztmp0 / ( ztmp1*ztmp1 ) 213 stab = SQRT( Cd ) ! Temporary array !!! (stab == SQRT(Cd)) 214 215 ztmp0 = (LOG(zu/10.) - zpsi_h_u) / vkarmn / sqrt_Cd_n10 216 ztmp2 = stab / sqrt_Cd_n10 ! (stab == SQRT(Cd)) 217 ztmp1 = 1. + Cx_n10*ztmp0 ! (Cx_n10 == Ch_n10) 218 Ch = Cx_n10*ztmp2 / ztmp1 ! L&Y 2004 eq. (10b) 219 220 Cx_n10 = 1.e-3 * (34.6 * sqrt_Cd_n10) ! L&Y 2004 eq. (6b) ! Cx_n10 == Ce_n10 221 Cen(:,:) = Cx_n10 222 ztmp1 = 1. + Cx_n10*ztmp0 223 Ce = Cx_n10*ztmp2 / ztmp1 ! L&Y 2004 eq. (10c) 224 ENDIF 224 225 ! 225 226 END DO … … 252 253 ! 253 254 ! When wind speed > 33 m/s => Cyclone conditions => special treatment 254 zgt33 = 0.5 + SIGN( 0.5, (zw - 33.) ) ! If pw10 < 33. => 0, else => 1255 zgt33 = 0.5_wp + SIGN( 0.5_wp, (zw - 33._wp) ) ! If pw10 < 33. => 0, else => 1 255 256 ! 256 257 cd_neutral_10m(ji,jj) = 1.e-3 * ( & … … 258 259 & + zgt33 * 2.34 ) ! wind >= 33 m/s 259 260 ! 260 cd_neutral_10m(ji,jj) = MAX(cd_neutral_10m(ji,jj), 1.E-6 )261 cd_neutral_10m(ji,jj) = MAX(cd_neutral_10m(ji,jj), 1.E-6_wp) 261 262 ! 262 263 END DO … … 285 286 DO jj = 1, jpj 286 287 DO ji = 1, jpi 287 zx2 = SQRT( ABS( 1. - 16.*pzeta(ji,jj) ) )288 zx2 = MAX ( zx2 , 1.)288 zx2 = SQRT( ABS( 1._wp - 16._wp*pzeta(ji,jj) ) ) 289 zx2 = MAX( zx2 , 1._wp ) 289 290 zx = SQRT( zx2 ) 290 zstab = 0.5 + SIGN( 0.5, pzeta(ji,jj) )291 ! 292 psi_m(ji,jj) = zstab * (-5. *pzeta(ji,jj)) & ! Stable293 & + (1. - zstab) * (2.*LOG((1. + zx)*0.5) & ! Unstable294 & + LOG((1. + zx2)*0.5) - 2.*ATAN(zx) + rpi*0.5) ! "291 zstab = 0.5_wp + SIGN( 0.5_wp , pzeta(ji,jj) ) 292 ! 293 psi_m(ji,jj) = zstab * (-5._wp*pzeta(ji,jj)) & ! Stable 294 & + (1._wp - zstab) * (2._wp*LOG((1._wp + zx)*0.5_wp) & ! Unstable 295 & + LOG((1._wp + zx2)*0.5_wp) - 2._wp*ATAN(zx) + rpi*0.5_wp) ! " 295 296 ! 296 297 END DO … … 319 320 DO jj = 1, jpj 320 321 DO ji = 1, jpi 321 zx2 = SQRT( ABS( 1. - 16.*pzeta(ji,jj) ) )322 zx2 = MAX ( zx2 , 1.)323 zstab = 0.5 + SIGN( 0.5, pzeta(ji,jj) )324 ! 325 psi_h(ji,jj) = zstab * (-5. *pzeta(ji,jj)) & ! Stable326 & + (1. - zstab) * (2.*LOG( (1. + zx2)*0.5)) ! Unstable322 zx2 = SQRT( ABS( 1._wp - 16._wp*pzeta(ji,jj) ) ) 323 zx2 = MAX( zx2 , 1._wp ) 324 zstab = 0.5_wp + SIGN( 0.5_wp , pzeta(ji,jj) ) 325 ! 326 psi_h(ji,jj) = zstab * (-5._wp*pzeta(ji,jj)) & ! Stable 327 & + (1._wp - zstab) * (2._wp*LOG( (1._wp + zx2)*0.5_wp )) ! Unstable 327 328 ! 328 329 END DO -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_phy.F90
r11182 r11209 8 8 !!---------------------------------------------------------------------- 9 9 10 !! virt_temp : virtual (aka sensible) temperature (potential or absolute) 10 11 !! rho_air : density of (moist) air (depends on T_air, q_air and SLP 12 !! visc_air : kinematic viscosity (aka Nu_air) of air from temperature 13 !! L_vap : latent heat of vaporization of water as a function of temperature 11 14 !! cp_air : specific heat of (moist) air (depends spec. hum. q_air) 15 !! gamma_moist : adiabatic lapse-rate of moist air 16 !! One_on_L : 1. / ( Monin-Obukhov length ) 17 !! Ri_bulk : bulk Richardson number aka BRN 12 18 !! q_sat : saturation humidity as a function of SLP and temperature 13 !! gamma_moist : 14 !! L_vap : latent heat of vaporization of water as a function of temperature 15 !! visc_air : kinematic viscosity (aka Nu_air) of air from temperature 16 19 17 20 USE dom_oce ! ocean space and time domain 18 21 USE phycst ! physical constants … … 21 24 PRIVATE 22 25 26 INTERFACE gamma_moist 27 MODULE PROCEDURE gamma_moist_vctr, gamma_moist_sclr 28 END INTERFACE gamma_moist 29 30 PUBLIC virt_temp 23 31 PUBLIC rho_air 32 PUBLIC visc_air 33 PUBLIC L_vap 24 34 PUBLIC cp_air 35 PUBLIC gamma_moist 36 PUBLIC One_on_L 37 PUBLIC Ri_bulk 25 38 PUBLIC q_sat 26 PUBLIC gamma_moist27 PUBLIC L_vap28 PUBLIC visc_air29 39 30 40 !!---------------------------------------------------------------------- … … 35 45 CONTAINS 36 46 47 FUNCTION virt_temp( pta, pqa ) 48 !!------------------------------------------------------------------------ 49 !! 50 !! Compute the (absolute/potential) virtual temperature, knowing the 51 !! (absolute/potential) temperature and specific humidity 52 !! 53 !! If input temperature is absolute then output vitual temperature is absolute 54 !! If input temperature is potential then output vitual temperature is potential 55 !! 56 !! Author: L. Brodeau, June 2019 / AeroBulk 57 !! (https://github.com/brodeau/aerobulk/) 58 !!------------------------------------------------------------------------ 59 REAL(wp), DIMENSION(jpi,jpj) :: virt_temp !: 1./(Monin Obukhov length) [m^-1] 60 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pta, & !: absolute or potetntial air temperature [K] 61 & pqa !: specific humidity of air [kg/kg] 62 !!------------------------------------------------------------------- 63 ! 64 virt_temp(:,:) = pta(:,:) * (1._wp + rctv0*pqa(:,:)) 65 !! 66 !! This is exactly the same sing that: 67 !! virt_temp = pta * ( pwa + reps0) / (reps0*(1.+pwa)) 68 !! with wpa (mixing ration) defined as : pwa = pqa/(1.-pqa) 69 ! 70 END FUNCTION virt_temp 71 37 72 FUNCTION rho_air( ptak, pqa, pslp ) 38 73 !!------------------------------------------------------------------------------- … … 41 76 !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere 42 77 !! 43 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)78 !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 44 79 !!------------------------------------------------------------------------------- 45 80 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K] … … 53 88 END FUNCTION rho_air 54 89 90 FUNCTION visc_air(ptak) 91 !!---------------------------------------------------------------------------------- 92 !! Air kinetic viscosity (m^2/s) given from temperature in degrees... 93 !! 94 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 95 !!---------------------------------------------------------------------------------- 96 REAL(wp), DIMENSION(jpi,jpj) :: visc_air ! 97 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature in (K) 98 ! 99 INTEGER :: ji, jj ! dummy loop indices 100 REAL(wp) :: ztc, ztc2 ! local scalar 101 !!---------------------------------------------------------------------------------- 102 ! 103 DO jj = 1, jpj 104 DO ji = 1, jpi 105 ztc = ptak(ji,jj) - rt0 ! air temp, in deg. C 106 ztc2 = ztc*ztc 107 visc_air(ji,jj) = 1.326e-5*(1. + 6.542E-3*ztc + 8.301e-6*ztc2 - 4.84e-9*ztc2*ztc) 108 END DO 109 END DO 110 ! 111 END FUNCTION visc_air 112 113 FUNCTION L_vap( psst ) 114 !!--------------------------------------------------------------------------------- 115 !! *** FUNCTION L_vap *** 116 !! 117 !! ** Purpose : Compute the latent heat of vaporization of water from temperature 118 !! 119 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 120 !!---------------------------------------------------------------------------------- 121 REAL(wp), DIMENSION(jpi,jpj) :: L_vap ! latent heat of vaporization [J/kg] 122 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psst ! water temperature [K] 123 !!---------------------------------------------------------------------------------- 124 ! 125 L_vap = ( 2.501 - 0.00237 * ( psst(:,:) - rt0) ) * 1.e6 126 ! 127 END FUNCTION L_vap 55 128 56 129 FUNCTION cp_air( pqa ) … … 69 142 ! 70 143 END FUNCTION cp_air 144 145 FUNCTION gamma_moist_vctr( ptak, pqa ) 146 !!---------------------------------------------------------------------------------- 147 !! *** FUNCTION gamma_moist_vctr *** 148 !! 149 !! ** Purpose : Compute the moist adiabatic lapse-rate. 150 !! => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate 151 !! => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html 152 !! 153 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 154 !!---------------------------------------------------------------------------------- 155 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K] 156 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! specific humidity [kg/kg] 157 REAL(wp), DIMENSION(jpi,jpj) :: gamma_moist_vctr ! moist adiabatic lapse-rate 158 ! 159 INTEGER :: ji, jj ! dummy loop indices 160 REAL(wp) :: zta, zqa, zwa, ziRT ! local scalar 161 !!---------------------------------------------------------------------------------- 162 ! 163 DO jj = 1, jpj 164 DO ji = 1, jpi 165 zta = MAX( ptak(ji,jj), 180._wp) ! prevents screw-up over masked regions where field == 0. 166 zqa = MAX( pqa(ji,jj), 1.E-6_wp) ! " " " 167 ! 168 zwa = zqa / (1. - zqa) ! w is mixing ratio w = q/(1-q) | q = w/(1+w) 169 ziRT = 1._wp/(R_dry*zta) ! 1/RT 170 gamma_moist_vctr(ji,jj) = grav * ( 1._wp + rLevap*zwa*ziRT ) / ( rCp_dry + rLevap*rLevap*zwa*reps0*ziRT/zta ) 171 END DO 172 END DO 173 ! 174 END FUNCTION gamma_moist_vctr 175 176 FUNCTION gamma_moist_sclr( ptak, pqa ) 177 !!---------------------------------------------------------------------------------- 178 !! ** Purpose : Compute the moist adiabatic lapse-rate. 179 !! => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate 180 !! => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html 181 !! 182 !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 183 !!---------------------------------------------------------------------------------- 184 REAL(wp) :: gamma_moist_sclr 185 REAL(wp), INTENT(in) :: ptak, pqa ! air temperature (K) and specific humidity (kg/kg) 186 ! 187 REAL(wp) :: zta, zqa, zwa, ziRT ! local scalar 188 !!---------------------------------------------------------------------------------- 189 zta = MAX( ptak, 180._wp) ! prevents screw-up over masked regions where field == 0. 190 zqa = MAX( pqa, 1.E-6_wp) ! " " " 191 !! 192 zwa = zqa / (1._wp - zqa) ! w is mixing ratio w = q/(1-q) | q = w/(1+w) 193 ziRT = 1._wp / (R_dry*zta) ! 1/RT 194 gamma_moist_sclr = grav * ( 1._wp + rLevap*zwa*ziRT ) / ( rCp_dry + rLevap*rLevap*zwa*reps0*ziRT/zta ) 195 !! 196 END FUNCTION gamma_moist_sclr 197 198 FUNCTION One_on_L( ptha, pqa, pus, pts, pqs ) 199 !!------------------------------------------------------------------------ 200 !! 201 !! Evaluates the 1./(Monin Obukhov length) from air temperature and 202 !! specific humidity, and frictional scales u*, t* and q* 203 !! 204 !! Author: L. Brodeau, June 2016 / AeroBulk 205 !! (https://github.com/brodeau/aerobulk/) 206 !!------------------------------------------------------------------------ 207 REAL(wp), DIMENSION(jpi,jpj) :: One_on_L !: 1./(Monin Obukhov length) [m^-1] 208 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptha, & !: average potetntial air temperature [K] 209 & pqa, & !: average specific humidity of air [kg/kg] 210 & pus, pts, pqs !: frictional velocity, temperature and humidity 211 ! 212 INTEGER :: ji, jj ! dummy loop indices 213 REAL(wp) :: zqa ! local scalar 214 !!------------------------------------------------------------------- 215 ! 216 DO jj = 1, jpj 217 DO ji = 1, jpi 218 ! 219 zqa = (1._wp + rctv0*pqa(ji,jj)) 220 ! 221 One_on_L(ji,jj) = grav*vkarmn*(pts(ji,jj) + rctv0*ptha(ji,jj)*pqs(ji,jj)) & 222 & / MAX( pus(ji,jj)*pus(ji,jj)*ptha(ji,jj)*zqa , 1.E-9_wp ) 223 ! 224 END DO 225 END DO 226 ! 227 One_on_L = SIGN( MIN(ABS(One_on_L),200._wp), One_on_L ) ! (prevent FPE from stupid values over masked regions...) 228 ! 229 END FUNCTION One_on_L 230 231 FUNCTION Ri_bulk( pz, psst, ptha, pssq, pqa, pub ) 232 !!---------------------------------------------------------------------------------- 233 !! Bulk Richardson number according to "wide-spread equation"... 234 !! 235 !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/) 236 !!---------------------------------------------------------------------------------- 237 REAL(wp), DIMENSION(jpi,jpj) :: Ri_bulk 238 REAL(wp) , INTENT(in) :: pz ! height above the sea (aka "delta z") [m] 239 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psst ! SST [K] 240 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptha ! pot. air temp. at height "pz" [K] 241 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pssq ! 0.98*q_sat(SST) [kg/kg] 242 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air spec. hum. at height "pz" [kg/kg] 243 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pub ! bulk wind speed [m/s] 244 ! 245 INTEGER :: ji, jj ! dummy loop indices 246 REAL(wp) :: zqa, zta, zgamma, zdth_v, ztv, zsstv ! local scalars 247 !!------------------------------------------------------------------- 248 ! 249 DO jj = 1, jpj 250 DO ji = 1, jpi 251 ! 252 zqa = 0.5_wp*(pqa(ji,jj)+pssq(ji,jj)) ! ~ mean q within the layer... 253 zta = 0.5_wp*( psst(ji,jj) + ptha(ji,jj) - gamma_moist(ptha(ji,jj),zqa)*pz ) ! ~ mean absolute temperature of air within the layer 254 zta = 0.5_wp*( psst(ji,jj) + ptha(ji,jj) - gamma_moist(zta, zqa)*pz ) ! ~ mean absolute temperature of air within the layer 255 zgamma = gamma_moist(zta, zqa) ! Adiabatic lapse-rate for moist air within the layer 256 ! 257 zsstv = psst(ji,jj)*(1._wp + rctv0*pssq(ji,jj)) ! absolute==potential virtual SST (absolute==potential because z=0!) 258 ! 259 zdth_v = ptha(ji,jj)*(1._wp + rctv0*pqa(ji,jj)) - zsstv ! air-sea delta of "virtual potential temperature" 260 ! 261 ztv = 0.5_wp*( zsstv + (ptha(ji,jj) - zgamma*pz)*(1._wp + rctv0*pqa(ji,jj)) ) ! ~ mean absolute virtual temp. within the layer 262 ! 263 Ri_bulk(ji,jj) = grav*zdth_v*pz / ( ztv*pub(ji,jj)*pub(ji,jj) ) ! the usual definition of Ri_bulk 264 ! 265 END DO 266 END DO 267 END FUNCTION Ri_bulk 71 268 72 269 … … 108 305 109 306 110 FUNCTION gamma_moist( ptak, pqa )111 !!----------------------------------------------------------------------------------112 !! *** FUNCTION gamma_moist ***113 !!114 !! ** Purpose : Compute the moist adiabatic lapse-rate.115 !! => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate116 !! => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html117 !!118 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)119 !!----------------------------------------------------------------------------------120 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K]121 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! specific humidity [kg/kg]122 REAL(wp), DIMENSION(jpi,jpj) :: gamma_moist ! moist adiabatic lapse-rate123 !124 INTEGER :: ji, jj ! dummy loop indices125 REAL(wp) :: zrv, ziRT ! local scalar126 !!----------------------------------------------------------------------------------127 !128 DO jj = 1, jpj129 DO ji = 1, jpi130 zrv = pqa(ji,jj) / (1. - pqa(ji,jj))131 ziRT = 1. / (R_dry*ptak(ji,jj)) ! 1/RT132 gamma_moist(ji,jj) = grav * ( 1. + rLevap*zrv*ziRT ) / ( rCp_dry + rLevap*rLevap*zrv*reps0*ziRT/ptak(ji,jj) )133 END DO134 END DO135 !136 END FUNCTION gamma_moist137 138 139 FUNCTION L_vap( psst )140 !!---------------------------------------------------------------------------------141 !! *** FUNCTION L_vap ***142 !!143 !! ** Purpose : Compute the latent heat of vaporization of water from temperature144 !!145 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)146 !!----------------------------------------------------------------------------------147 REAL(wp), DIMENSION(jpi,jpj) :: L_vap ! latent heat of vaporization [J/kg]148 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psst ! water temperature [K]149 !!----------------------------------------------------------------------------------150 !151 L_vap = ( 2.501 - 0.00237 * ( psst(:,:) - rt0) ) * 1.e6152 !153 END FUNCTION L_vap154 155 156 157 FUNCTION visc_air(ptak)158 !!----------------------------------------------------------------------------------159 !! Air kinetic viscosity (m^2/s) given from temperature in degrees...160 !!161 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)162 !!----------------------------------------------------------------------------------163 REAL(wp), DIMENSION(jpi,jpj) :: visc_air !164 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature in (K)165 !166 INTEGER :: ji, jj ! dummy loop indices167 REAL(wp) :: ztc, ztc2 ! local scalar168 !!----------------------------------------------------------------------------------169 !170 DO jj = 1, jpj171 DO ji = 1, jpi172 ztc = ptak(ji,jj) - rt0 ! air temp, in deg. C173 ztc2 = ztc*ztc174 visc_air(ji,jj) = 1.326e-5*(1. + 6.542E-3*ztc + 8.301e-6*ztc2 - 4.84e-9*ztc2*ztc)175 END DO176 END DO177 !178 END FUNCTION visc_air179 180 181 307 !!====================================================================== 182 308 END MODULE sbcblk_phy
Note: See TracChangeset
for help on using the changeset viewer.
