New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 11209 for NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_phy.F90 – NEMO

Ignore:
Timestamp:
2019-07-03T09:09:17+02:00 (5 years ago)
Author:
laurent
Message:

LB: making more efficient and more centralized use of physical/thermodynamics functions defined into "OCE/SBC/sbcblk_phy.F90".

File:
1 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_phy.F90

    r11182 r11209  
    88   !!---------------------------------------------------------------------- 
    99 
     10   !!   virt_temp     : virtual (aka sensible) temperature (potential or absolute) 
    1011   !!   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 
    1114   !!   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 
    1218   !!   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 
    1720   USE dom_oce        ! ocean space and time domain 
    1821   USE phycst         ! physical constants 
    … …  
    2124   PRIVATE 
    2225 
     26   INTERFACE gamma_moist 
     27      MODULE PROCEDURE gamma_moist_vctr, gamma_moist_sclr 
     28   END INTERFACE gamma_moist 
     29 
     30   PUBLIC virt_temp 
    2331   PUBLIC rho_air 
     32   PUBLIC visc_air 
     33   PUBLIC L_vap 
    2434   PUBLIC cp_air 
     35   PUBLIC gamma_moist 
     36   PUBLIC One_on_L 
     37   PUBLIC Ri_bulk 
    2538   PUBLIC q_sat 
    26    PUBLIC gamma_moist 
    27    PUBLIC L_vap 
    28    PUBLIC visc_air 
    2939 
    3040   !!---------------------------------------------------------------------- 
    … …  
    3545CONTAINS 
    3646 
     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 
    3772   FUNCTION rho_air( ptak, pqa, pslp ) 
    3873      !!------------------------------------------------------------------------------- 
    … …  
    4176      !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere 
    4277      !! 
    43       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
     78      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    4479      !!------------------------------------------------------------------------------- 
    4580      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak      ! air temperature             [K] 
    … …  
    5388   END FUNCTION rho_air 
    5489 
     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 
    55128 
    56129   FUNCTION cp_air( pqa ) 
    … …  
    69142      ! 
    70143   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 
    71268 
    72269 
    … …  
    108305 
    109306 
    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_rate 
    116       !!     => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html 
    117       !! 
    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-rate 
    123       ! 
    124       INTEGER  ::   ji, jj         ! dummy loop indices 
    125       REAL(wp) :: zrv, ziRT        ! local scalar 
    126       !!---------------------------------------------------------------------------------- 
    127       ! 
    128       DO jj = 1, jpj 
    129          DO ji = 1, jpi 
    130             zrv = pqa(ji,jj) / (1. - pqa(ji,jj)) 
    131             ziRT = 1. / (R_dry*ptak(ji,jj))    ! 1/RT 
    132             gamma_moist(ji,jj) = grav * ( 1. + rLevap*zrv*ziRT ) / ( rCp_dry + rLevap*rLevap*zrv*reps0*ziRT/ptak(ji,jj) ) 
    133          END DO 
    134       END DO 
    135       ! 
    136    END FUNCTION gamma_moist 
    137  
    138  
    139    FUNCTION L_vap( psst ) 
    140       !!--------------------------------------------------------------------------------- 
    141       !!                           ***  FUNCTION L_vap  *** 
    142       !! 
    143       !! ** Purpose : Compute the latent heat of vaporization of water from temperature 
    144       !! 
    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.e6 
    152       ! 
    153    END FUNCTION L_vap 
    154  
    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 indices 
    167       REAL(wp) ::   ztc, ztc2   ! local scalar 
    168       !!---------------------------------------------------------------------------------- 
    169       ! 
    170       DO jj = 1, jpj 
    171          DO ji = 1, jpi 
    172             ztc  = ptak(ji,jj) - rt0   ! air temp, in deg. C 
    173             ztc2 = ztc*ztc 
    174             visc_air(ji,jj) = 1.326e-5*(1. + 6.542E-3*ztc + 8.301e-6*ztc2 - 4.84e-9*ztc2*ztc) 
    175          END DO 
    176       END DO 
    177       ! 
    178    END FUNCTION visc_air 
    179     
    180  
    181307   !!====================================================================== 
    182308END MODULE sbcblk_phy 
Note: See TracChangeset for help on using the changeset viewer.