MODULE sbcblk_phy !!====================================================================== !! *** MODULE sbcblk_phy *** !! A set of functions to compute air themodynamics parameters !! needed by Aerodynamic Bulk Formulas !!===================================================================== !! 4.0 ! 2019 L. Brodeau from AeroBulk package (https://github.com/brodeau/aerobulk/) !!---------------------------------------------------------------------- !! virt_temp : virtual (aka sensible) temperature (potential or absolute) !! rho_air : density of (moist) air (depends on T_air, q_air and SLP !! visc_air : kinematic viscosity (aka Nu_air) of air from temperature !! L_vap : latent heat of vaporization of water as a function of temperature !! cp_air : specific heat of (moist) air (depends spec. hum. q_air) !! gamma_moist : adiabatic lapse-rate of moist air !! One_on_L : 1. / ( Monin-Obukhov length ) !! Ri_bulk : bulk Richardson number aka BRN !! q_sat : saturation humidity as a function of SLP and temperature !! q_air_rh : specific humidity as a function of RH, t_air and SLP USE dom_oce ! ocean space and time domain USE phycst ! physical constants IMPLICIT NONE PRIVATE INTERFACE gamma_moist MODULE PROCEDURE gamma_moist_vctr, gamma_moist_sclr END INTERFACE gamma_moist PUBLIC virt_temp PUBLIC rho_air PUBLIC visc_air PUBLIC L_vap PUBLIC cp_air PUBLIC gamma_moist PUBLIC One_on_L PUBLIC Ri_bulk PUBLIC q_sat PUBLIC q_air_rh !!---------------------------------------------------------------------- !! NEMO/OCE 4.0 , NEMO Consortium (2018) !! $Id: sbcblk.F90 10535 2019-01-16 17:36:47Z clem $ !! Software governed by the CeCILL license (see ./LICENSE) !!---------------------------------------------------------------------- CONTAINS FUNCTION virt_temp( pta, pqa ) !!------------------------------------------------------------------------ !! !! Compute the (absolute/potential) virtual temperature, knowing the !! (absolute/potential) temperature and specific humidity !! !! If input temperature is absolute then output vitual temperature is absolute !! If input temperature is potential then output vitual temperature is potential !! !! Author: L. Brodeau, June 2019 / AeroBulk !! (https://github.com/brodeau/aerobulk/) !!------------------------------------------------------------------------ REAL(wp), DIMENSION(jpi,jpj) :: virt_temp !: 1./(Monin Obukhov length) [m^-1] REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pta, & !: absolute or potetntial air temperature [K] & pqa !: specific humidity of air [kg/kg] !!------------------------------------------------------------------- ! virt_temp(:,:) = pta(:,:) * (1._wp + rctv0*pqa(:,:)) !! !! This is exactly the same sing that: !! virt_temp = pta * ( pwa + reps0) / (reps0*(1.+pwa)) !! with wpa (mixing ration) defined as : pwa = pqa/(1.-pqa) ! END FUNCTION virt_temp FUNCTION rho_air( ptak, pqa, pslp ) !!------------------------------------------------------------------------------- !! *** FUNCTION rho_air *** !! !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere !! !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) !!------------------------------------------------------------------------------- REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K] REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air specific humidity [kg/kg] REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pslp ! pressure in [Pa] REAL(wp), DIMENSION(jpi,jpj) :: rho_air ! density of moist air [kg/m^3] !!------------------------------------------------------------------------------- ! rho_air = pslp / ( R_dry*ptak * ( 1._wp + rctv0*pqa ) ) ! END FUNCTION rho_air FUNCTION visc_air(ptak) !!---------------------------------------------------------------------------------- !! Air kinetic viscosity (m^2/s) given from temperature in degrees... !! !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) !!---------------------------------------------------------------------------------- REAL(wp), DIMENSION(jpi,jpj) :: visc_air ! REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature in (K) ! INTEGER :: ji, jj ! dummy loop indices REAL(wp) :: ztc, ztc2 ! local scalar !!---------------------------------------------------------------------------------- ! DO jj = 1, jpj DO ji = 1, jpi ztc = ptak(ji,jj) - rt0 ! air temp, in deg. C ztc2 = ztc*ztc visc_air(ji,jj) = 1.326e-5*(1. + 6.542E-3*ztc + 8.301e-6*ztc2 - 4.84e-9*ztc2*ztc) END DO END DO ! END FUNCTION visc_air FUNCTION L_vap( psst ) !!--------------------------------------------------------------------------------- !! *** FUNCTION L_vap *** !! !! ** Purpose : Compute the latent heat of vaporization of water out of temperature !! !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) !!---------------------------------------------------------------------------------- REAL(wp), DIMENSION(jpi,jpj) :: L_vap ! latent heat of vaporization [J/kg] REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psst ! water temperature [K] !!---------------------------------------------------------------------------------- ! L_vap = ( 2.501 - 0.00237 * ( psst(:,:) - rt0) ) * 1.e6 ! END FUNCTION L_vap FUNCTION cp_air( pqa ) !!------------------------------------------------------------------------------- !! *** FUNCTION cp_air *** !! !! ** Purpose : Compute specific heat (Cp) of moist air !! !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) !!------------------------------------------------------------------------------- REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air specific humidity [kg/kg] REAL(wp), DIMENSION(jpi,jpj) :: cp_air ! specific heat of moist air [J/K/kg] !!------------------------------------------------------------------------------- ! cp_air = rCp_dry + rCp_vap * pqa ! END FUNCTION cp_air FUNCTION gamma_moist_vctr( ptak, pqa ) !!---------------------------------------------------------------------------------- !! *** FUNCTION gamma_moist_vctr *** !! !! ** Purpose : Compute the moist adiabatic lapse-rate. !! => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate !! => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html !! !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) !!---------------------------------------------------------------------------------- REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K] REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! specific humidity [kg/kg] REAL(wp), DIMENSION(jpi,jpj) :: gamma_moist_vctr ! moist adiabatic lapse-rate ! INTEGER :: ji, jj ! dummy loop indices REAL(wp) :: zta, zqa, zwa, ziRT ! local scalar !!---------------------------------------------------------------------------------- ! DO jj = 1, jpj DO ji = 1, jpi zta = MAX( ptak(ji,jj), 180._wp) ! prevents screw-up over masked regions where field == 0. zqa = MAX( pqa(ji,jj), 1.E-6_wp) ! " " " ! zwa = zqa / (1. - zqa) ! w is mixing ratio w = q/(1-q) | q = w/(1+w) ziRT = 1._wp/(R_dry*zta) ! 1/RT gamma_moist_vctr(ji,jj) = grav * ( 1._wp + rLevap*zwa*ziRT ) / ( rCp_dry + rLevap*rLevap*zwa*reps0*ziRT/zta ) END DO END DO ! END FUNCTION gamma_moist_vctr FUNCTION gamma_moist_sclr( ptak, pqa ) !!---------------------------------------------------------------------------------- !! ** Purpose : Compute the moist adiabatic lapse-rate. !! => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate !! => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html !! !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) !!---------------------------------------------------------------------------------- REAL(wp) :: gamma_moist_sclr REAL(wp), INTENT(in) :: ptak, pqa ! air temperature (K) and specific humidity (kg/kg) ! REAL(wp) :: zta, zqa, zwa, ziRT ! local scalar !!---------------------------------------------------------------------------------- zta = MAX( ptak, 180._wp) ! prevents screw-up over masked regions where field == 0. zqa = MAX( pqa, 1.E-6_wp) ! " " " !! zwa = zqa / (1._wp - zqa) ! w is mixing ratio w = q/(1-q) | q = w/(1+w) ziRT = 1._wp / (R_dry*zta) ! 1/RT gamma_moist_sclr = grav * ( 1._wp + rLevap*zwa*ziRT ) / ( rCp_dry + rLevap*rLevap*zwa*reps0*ziRT/zta ) !! END FUNCTION gamma_moist_sclr FUNCTION One_on_L( ptha, pqa, pus, pts, pqs ) !!------------------------------------------------------------------------ !! !! Evaluates the 1./(Monin Obukhov length) from air temperature and !! specific humidity, and frictional scales u*, t* and q* !! !! Author: L. Brodeau, June 2016 / AeroBulk !! (https://github.com/brodeau/aerobulk/) !!------------------------------------------------------------------------ REAL(wp), DIMENSION(jpi,jpj) :: One_on_L !: 1./(Monin Obukhov length) [m^-1] REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptha, & !: average potetntial air temperature [K] & pqa, & !: average specific humidity of air [kg/kg] & pus, pts, pqs !: frictional velocity, temperature and humidity ! INTEGER :: ji, jj ! dummy loop indices REAL(wp) :: zqa ! local scalar !!------------------------------------------------------------------- ! DO jj = 1, jpj DO ji = 1, jpi ! zqa = (1._wp + rctv0*pqa(ji,jj)) ! ! The main concern is to know whether, the vertical turbulent flux of virtual temperature, < u' theta_v' > is estimated with: ! a/ -u* [ theta* (1 + 0.61 q) + 0.61 theta q* ] => this is the one that seems correct! chose this one! ! or ! b/ -u* [ theta* + 0.61 theta q* ] ! One_on_L(ji,jj) = grav*vkarmn*( pts(ji,jj)*zqa + rctv0*ptha(ji,jj)*pqs(ji,jj) ) & & / MAX( pus(ji,jj)*pus(ji,jj)*ptha(ji,jj)*zqa , 1.E-9_wp ) ! END DO END DO ! One_on_L = SIGN( MIN(ABS(One_on_L),200._wp), One_on_L ) ! (prevent FPE from stupid values over masked regions...) ! END FUNCTION One_on_L FUNCTION Ri_bulk( pz, psst, ptha, pssq, pqa, pub ) !!---------------------------------------------------------------------------------- !! Bulk Richardson number according to "wide-spread equation"... !! !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/) !!---------------------------------------------------------------------------------- REAL(wp), DIMENSION(jpi,jpj) :: Ri_bulk REAL(wp) , INTENT(in) :: pz ! height above the sea (aka "delta z") [m] REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psst ! SST [K] REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptha ! pot. air temp. at height "pz" [K] REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pssq ! 0.98*q_sat(SST) [kg/kg] REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air spec. hum. at height "pz" [kg/kg] REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pub ! bulk wind speed [m/s] ! INTEGER :: ji, jj ! dummy loop indices REAL(wp) :: zqa, zta, zgamma, zdth_v, ztv, zsstv ! local scalars !!------------------------------------------------------------------- ! DO jj = 1, jpj DO ji = 1, jpi ! zqa = 0.5_wp*(pqa(ji,jj)+pssq(ji,jj)) ! ~ mean q within the layer... 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 zta = 0.5_wp*( psst(ji,jj) + ptha(ji,jj) - gamma_moist(zta, zqa)*pz ) ! ~ mean absolute temperature of air within the layer zgamma = gamma_moist(zta, zqa) ! Adiabatic lapse-rate for moist air within the layer ! zsstv = psst(ji,jj)*(1._wp + rctv0*pssq(ji,jj)) ! absolute==potential virtual SST (absolute==potential because z=0!) ! zdth_v = ptha(ji,jj)*(1._wp + rctv0*pqa(ji,jj)) - zsstv ! air-sea delta of "virtual potential temperature" ! ztv = 0.5_wp*( zsstv + (ptha(ji,jj) - zgamma*pz)*(1._wp + rctv0*pqa(ji,jj)) ) ! ~ mean absolute virtual temp. within the layer ! Ri_bulk(ji,jj) = grav*zdth_v*pz / ( ztv*pub(ji,jj)*pub(ji,jj) ) ! the usual definition of Ri_bulk ! END DO END DO END FUNCTION Ri_bulk FUNCTION e_sat_sclr( ptak, pslp ) !!---------------------------------------------------------------------------------- !! *** FUNCTION e_sat_sclr *** !! < SCALAR argument version > !! ** Purpose : water vapor at saturation in [Pa] !! Based on accurate estimate by Goff, 1957 !! !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) !! !! Note: what rt0 should be here, is 273.16 (triple point of water) and not 273.15 like here !!---------------------------------------------------------------------------------- REAL(wp), INTENT(in) :: ptak ! air temperature [K] REAL(wp), INTENT(in) :: pslp ! sea level atmospheric pressure [Pa] REAL(wp) :: e_sat_sclr ! water vapor at saturation [kg/kg] ! REAL(wp) :: zta, ztmp ! local scalar !!---------------------------------------------------------------------------------- ! zta = MAX( ptak , 180._wp ) ! air temp., prevents fpe0 errors dute to unrealistically low values over masked regions... ztmp = rt0 / zta ! ! Vapour pressure at saturation [Pa] : WMO, (Goff, 1957) e_sat_sclr = 100.*( 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(zta/rt0) & & + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(zta/rt0 - 1.)) ) & & + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614) ) ! END FUNCTION e_sat_sclr FUNCTION q_sat( ptak, pslp ) !!---------------------------------------------------------------------------------- !! *** FUNCTION q_sat *** !! !! ** Purpose : Specific humidity at saturation in [kg/kg] !! Based on accurate estimate of "e_sat" !! aka saturation water vapor (Goff, 1957) !! !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) !!---------------------------------------------------------------------------------- REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K] REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pslp ! sea level atmospheric pressure [Pa] REAL(wp), DIMENSION(jpi,jpj) :: q_sat ! Specific humidity at saturation [kg/kg] ! INTEGER :: ji, jj ! dummy loop indices REAL(wp) :: ze_sat ! local scalar !!---------------------------------------------------------------------------------- ! DO jj = 1, jpj DO ji = 1, jpi ! ze_sat = e_sat_sclr( ptak(ji,jj), pslp(ji,jj) ) ! q_sat(ji,jj) = reps0 * ze_sat/( pslp(ji,jj) - (1._wp - reps0)*ze_sat ) ! END DO END DO ! END FUNCTION q_sat FUNCTION q_air_rh(prha, ptak, pslp) !!---------------------------------------------------------------------------------- !! Specific humidity of air out of Relative Humidity !! !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) !!---------------------------------------------------------------------------------- REAL(wp), DIMENSION(jpi,jpj) :: q_air_rh REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: prha !: relative humidity [fraction, not %!!!] REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak !: air temperature [K] REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pslp !: atmospheric pressure [Pa] ! INTEGER :: ji, jj ! dummy loop indices REAL(wp) :: ze ! local scalar !!---------------------------------------------------------------------------------- ! DO jj = 1, jpj DO ji = 1, jpi ze = prha(ji,jj)*e_sat_sclr(ptak(ji,jj), pslp(ji,jj)) q_air_rh(ji,jj) = ze*reps0/(pslp(ji,jj) - (1. - reps0)*ze) END DO END DO ! END FUNCTION q_air_rh !!====================================================================== END MODULE sbcblk_phy