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 11615 for NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_phy.F90 – NEMO

Ignore:
Timestamp:
2019-09-30T15:18:21+02:00 (5 years ago)
Author:
laurent
Message:

LB: new cool-skin and warm-layer parameterizations for ECMWF and COARE3.6, COARE3.0 uses same CSWL param as for COARE3.6.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r11284 r11615  
    2929   END INTERFACE gamma_moist 
    3030 
     31   INTERFACE e_sat 
     32      MODULE PROCEDURE e_sat_vctr, e_sat_sclr 
     33   END INTERFACE e_sat 
     34 
     35   INTERFACE L_vap 
     36      MODULE PROCEDURE L_vap_vctr, L_vap_sclr 
     37   END INTERFACE L_vap 
     38 
     39   INTERFACE rho_air 
     40      MODULE PROCEDURE rho_air_vctr, rho_air_sclr 
     41   END INTERFACE rho_air 
     42 
     43   INTERFACE cp_air 
     44      MODULE PROCEDURE cp_air_vctr, cp_air_sclr 
     45   END INTERFACE cp_air 
     46 
     47      INTERFACE alpha_sw 
     48      MODULE PROCEDURE alpha_sw_vctr, alpha_sw_sclr 
     49   END INTERFACE alpha_sw 
     50 
     51 
     52 
    3153   PUBLIC virt_temp 
    3254   PUBLIC rho_air 
    … …  
    3961   PUBLIC q_sat 
    4062   PUBLIC q_air_rh 
     63   !: 
     64   PUBLIC update_qnsol_tau 
     65   PUBLIC alpha_sw 
    4166 
    4267   !!---------------------------------------------------------------------- 
    … …  
    7297   END FUNCTION virt_temp 
    7398 
    74    FUNCTION rho_air( ptak, pqa, pslp ) 
    75       !!------------------------------------------------------------------------------- 
    76       !!                           ***  FUNCTION rho_air  *** 
     99   FUNCTION rho_air_vctr( ptak, pqa, pslp ) 
     100      !!------------------------------------------------------------------------------- 
     101      !!                           ***  FUNCTION rho_air_vctr  *** 
    77102      !! 
    78103      !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere 
    … …  
    83108      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa       ! air specific humidity   [kg/kg] 
    84109      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp      ! pressure in                [Pa] 
    85       REAL(wp), DIMENSION(jpi,jpj)             ::   rho_air   ! density of moist air   [kg/m^3] 
    86       !!------------------------------------------------------------------------------- 
    87       ! 
    88       rho_air = pslp / (  R_dry*ptak * ( 1._wp + rctv0*pqa )  ) 
    89       ! 
    90    END FUNCTION rho_air 
     110      REAL(wp), DIMENSION(jpi,jpj)             ::   rho_air_vctr   ! density of moist air   [kg/m^3] 
     111      !!------------------------------------------------------------------------------- 
     112      rho_air_vctr = MAX( pslp / (R_dry*ptak * ( 1._wp + rctv0*pqa )) , 0.8_wp ) 
     113   END FUNCTION rho_air_vctr 
     114 
     115   FUNCTION rho_air_sclr( ptak, pqa, pslp ) 
     116      !!------------------------------------------------------------------------------- 
     117      !!                           ***  FUNCTION rho_air_sclr  *** 
     118      !! 
     119      !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere 
     120      !! 
     121      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
     122      !!------------------------------------------------------------------------------- 
     123      REAL(wp), INTENT(in) :: ptak           ! air temperature             [K] 
     124      REAL(wp), INTENT(in) :: pqa            ! air specific humidity   [kg/kg] 
     125      REAL(wp), INTENT(in) :: pslp           ! pressure in                [Pa] 
     126      REAL(wp)             :: rho_air_sclr   ! density of moist air   [kg/m^3] 
     127      !!------------------------------------------------------------------------------- 
     128      rho_air_sclr = MAX( pslp / (R_dry*ptak * ( 1._wp + rctv0*pqa )) , 0.8_wp ) 
     129   END FUNCTION rho_air_sclr 
     130    
     131 
    91132 
    92133   FUNCTION visc_air(ptak) 
    … …  
    113154   END FUNCTION visc_air 
    114155 
    115    FUNCTION L_vap( psst ) 
     156   FUNCTION L_vap_vctr( psst ) 
    116157      !!--------------------------------------------------------------------------------- 
    117       !!                           ***  FUNCTION L_vap  *** 
    118       !! 
    119       !! ** Purpose : Compute the latent heat of vaporization of water out of temperature 
    120       !! 
    121       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    122       !!---------------------------------------------------------------------------------- 
    123       REAL(wp), DIMENSION(jpi,jpj)             ::   L_vap   ! latent heat of vaporization   [J/kg] 
     158      !!                           ***  FUNCTION L_vap_vctr  *** 
     159      !! 
     160      !! ** Purpose : Compute the latent heat of vaporization of water from temperature 
     161      !! 
     162      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
     163      !!---------------------------------------------------------------------------------- 
     164      REAL(wp), DIMENSION(jpi,jpj)             ::   L_vap_vctr   ! latent heat of vaporization   [J/kg] 
    124165      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   psst   ! water temperature                [K] 
    125166      !!---------------------------------------------------------------------------------- 
    126167      ! 
    127       L_vap = (  2.501 - 0.00237 * ( psst(:,:) - rt0)  ) * 1.e6 
    128       ! 
    129    END FUNCTION L_vap 
    130  
    131    FUNCTION cp_air( pqa ) 
    132       !!------------------------------------------------------------------------------- 
    133       !!                           ***  FUNCTION cp_air  *** 
     168      L_vap_vctr = (  2.501_wp - 0.00237_wp * ( psst(:,:) - rt0)  ) * 1.e6_wp 
     169      ! 
     170   END FUNCTION L_vap_vctr 
     171 
     172   FUNCTION L_vap_sclr( psst ) 
     173      !!--------------------------------------------------------------------------------- 
     174      !!                           ***  FUNCTION L_vap_sclr  *** 
     175      !! 
     176      !! ** Purpose : Compute the latent heat of vaporization of water from temperature 
     177      !! 
     178      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
     179      !!---------------------------------------------------------------------------------- 
     180      REAL(wp)             ::   L_vap_sclr   ! latent heat of vaporization   [J/kg] 
     181      REAL(wp), INTENT(in) ::   psst         ! water temperature                [K] 
     182      !!---------------------------------------------------------------------------------- 
     183      ! 
     184      L_vap_sclr = (  2.501_wp - 0.00237_wp * ( psst - rt0)  ) * 1.e6_wp 
     185      ! 
     186   END FUNCTION L_vap_sclr 
     187 
     188 
     189   FUNCTION cp_air_vctr( pqa ) 
     190      !!------------------------------------------------------------------------------- 
     191      !!                           ***  FUNCTION cp_air_vctr  *** 
    134192      !! 
    135193      !! ** Purpose : Compute specific heat (Cp) of moist air 
    … …  
    138196      !!------------------------------------------------------------------------------- 
    139197      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa      ! air specific humidity         [kg/kg] 
    140       REAL(wp), DIMENSION(jpi,jpj)             ::   cp_air   ! specific heat of moist air   [J/K/kg] 
    141       !!------------------------------------------------------------------------------- 
    142       ! 
    143       cp_air = rCp_dry + rCp_vap * pqa 
    144       ! 
    145    END FUNCTION cp_air 
     198      REAL(wp), DIMENSION(jpi,jpj)             ::   cp_air_vctr   ! specific heat of moist air   [J/K/kg] 
     199      !!------------------------------------------------------------------------------- 
     200      cp_air_vctr = rCp_dry + rCp_vap * pqa 
     201   END FUNCTION cp_air_vctr 
     202 
     203   FUNCTION cp_air_sclr( pqa ) 
     204      !!------------------------------------------------------------------------------- 
     205      !!                           ***  FUNCTION cp_air_sclr  *** 
     206      !! 
     207      !! ** Purpose : Compute specific heat (Cp) of moist air 
     208      !! 
     209      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
     210      !!------------------------------------------------------------------------------- 
     211      REAL(wp), INTENT(in) :: pqa           ! air specific humidity         [kg/kg] 
     212      REAL(wp)             :: cp_air_sclr   ! specific heat of moist air   [J/K/kg] 
     213      !!------------------------------------------------------------------------------- 
     214      cp_air_sclr = rCp_dry + rCp_vap * pqa 
     215   END FUNCTION cp_air_sclr 
     216 
     217 
     218 
     219 
    146220 
    147221   FUNCTION gamma_moist_vctr( ptak, pqa ) 
    … …  
    275349 
    276350 
    277    FUNCTION e_sat_sclr( ptak, pslp ) 
     351   FUNCTION e_sat_vctr(ptak) 
     352      !!************************************************** 
     353      !! ptak:     air temperature [K] 
     354      !! e_sat:  water vapor at saturation [Pa] 
     355      !! 
     356      !! Recommended by WMO 
     357      !! 
     358      !! Goff, J. A., 1957: Saturation pressure of water on the new kelvin 
     359      !! temperature scale. Transactions of the American society of heating 
     360      !! and ventilating engineers, 347–354. 
     361      !! 
     362      !! rt0 should be 273.16 (triple point of water) and not 273.15 like here 
     363      !!************************************************** 
     364 
     365      REAL(wp), DIMENSION(jpi,jpj)             :: e_sat_vctr !: vapour pressure at saturation  [Pa] 
     366      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak    !: temperature (K) 
     367 
     368      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ztmp 
     369 
     370      ALLOCATE ( ztmp(jpi,jpj) ) 
     371 
     372      ztmp(:,:) = rtt0/ptak(:,:) 
     373 
     374      e_sat_vctr = 100.*( 10.**(10.79574*(1. - ztmp) - 5.028*LOG10(ptak/rtt0)         & 
     375         &       + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak/rtt0 - 1.)) )   & 
     376         &       + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614) ) 
     377 
     378      DEALLOCATE ( ztmp ) 
     379 
     380   END FUNCTION e_sat_vctr 
     381 
     382 
     383   FUNCTION e_sat_sclr( ptak ) 
    278384      !!---------------------------------------------------------------------------------- 
    279385      !!                   ***  FUNCTION e_sat_sclr  *** 
    … …  
    287393      !!---------------------------------------------------------------------------------- 
    288394      REAL(wp), INTENT(in) ::   ptak    ! air temperature                  [K] 
    289       REAL(wp), INTENT(in) ::   pslp    ! sea level atmospheric pressure   [Pa] 
    290395      REAL(wp)             ::   e_sat_sclr   ! water vapor at saturation   [kg/kg] 
    291396      ! 
    … …  
    325430         DO ji = 1, jpi 
    326431            ! 
    327             ze_sat =  e_sat_sclr( ptak(ji,jj), pslp(ji,jj) ) 
     432            ze_sat =  e_sat_sclr( ptak(ji,jj) ) 
    328433            ! 
    329434            q_sat(ji,jj) = reps0 * ze_sat/( pslp(ji,jj) - (1._wp - reps0)*ze_sat ) 
    … …  
    351456      DO jj = 1, jpj 
    352457         DO ji = 1, jpi 
    353             ze = prha(ji,jj)*e_sat_sclr(ptak(ji,jj), pslp(ji,jj)) 
     458            ze = prha(ji,jj)*e_sat_sclr(ptak(ji,jj)) 
    354459            q_air_rh(ji,jj) = ze*reps0/(pslp(ji,jj) - (1. - reps0)*ze) 
    355460         END DO 
    … …  
    357462      ! 
    358463   END FUNCTION q_air_rh 
     464 
     465 
     466   SUBROUTINE UPDATE_QNSOL_TAU( pTs, pqs, pTa, pqa, pust, ptst, pqst, pUb, pslp, prlw, & 
     467      &                         pQns, pTau,  & 
     468      &                         Qlat) 
     469      !!---------------------------------------------------------------------------------- 
     470      !! Purpose: returns the non-solar heat flux to the ocean aka "Qlat + Qsen + Qlw" 
     471      !!          and the module of the wind stress => pTau = Tau 
     472      !! ** Author: L. Brodeau, Sept. 2019 / AeroBulk (https://github.com/brodeau/aerobulk/) 
     473      !!---------------------------------------------------------------------------------- 
     474      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTs  ! water temperature at the air-sea interface [K] 
     475      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pqs  ! satur. spec. hum. at T=pTs   [kg/kg] 
     476      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTa  ! air temperature at z=zu [K] 
     477      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pqa  ! specific humidity at z=zu [kg/kg] 
     478      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pust ! u* 
     479      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: ptst ! t* 
     480      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pqst ! q* 
     481      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pUb  ! bulk wind speed at z=zu [m/s] 
     482      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pslp ! sea-level atmospheric pressure [Pa] 
     483      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: prlw ! downwelling longwave radiative flux [W/m^2] 
     484      ! 
     485      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pQns ! non-solar heat flux to the ocean aka "Qlat + Qsen + Qlw" [W/m^2]] 
     486      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pTau ! module of the wind stress [N/m^2] 
     487      ! 
     488      REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(out) :: Qlat 
     489      ! 
     490      REAL(wp) :: zdt, zdq, zCd, zCh, zCe, zUrho, zTs2, zz0, & 
     491         &        zQlat, zQsen, zQlw 
     492      INTEGER  ::   ji, jj     ! dummy loop indices 
     493      !!----------------------------------------------------------------------------------      
     494      DO jj = 1, jpj 
     495         DO ji = 1, jpi 
     496 
     497            zdt = pTa(ji,jj) - pTs(ji,jj) ;  zdt = SIGN( MAX(ABS(zdt),1.E-6_wp), zdt ) 
     498            zdq = pqa(ji,jj) - pqs(ji,jj) ;  zdq = SIGN( MAX(ABS(zdq),1.E-9_wp), zdq ) 
     499            zz0 = pust(ji,jj)/pUb(ji,jj) 
     500            zCd = zz0*zz0 
     501            zCh = zz0*ptst(ji,jj)/zdt 
     502            zCe = zz0*pqst(ji,jj)/zdq 
     503 
     504            zUrho = pUb(ji,jj)*MAX(rho_air(pTa(ji,jj), pqa(ji,jj), pslp(ji,jj)), 1._wp)     ! rho*U10 
     505            zTs2  = pTs(ji,jj)*pTs(ji,jj) 
     506 
     507            ! Wind stress module: 
     508            pTau(ji,jj) = zCd*zUrho*pUb(ji,jj) ! lolo? 
     509 
     510            ! Non-Solar heat flux to the ocean: 
     511            zQlat = MIN ( zUrho*zCe*L_vap( pTs(ji,jj)) * zdq , 0._wp )  ! we do not want a Qlat > 0 ! 
     512            zQsen = zUrho*zCh*cp_air(pqa(ji,jj)) * zdt 
     513            zQlw  = emiss_w*(prlw(ji,jj) - stefan*zTs2*zTs2) ! Net longwave flux 
     514 
     515            pQns(ji,jj) = zQlat + zQsen + zQlw 
     516             
     517            IF ( PRESENT(Qlat) ) Qlat(ji,jj) = zQlat             
     518         END DO 
     519      END DO 
     520   END SUBROUTINE UPDATE_QNSOL_TAU 
     521 
     522 
     523   FUNCTION alpha_sw_vctr( psst ) 
     524      !!--------------------------------------------------------------------------------- 
     525      !!                           ***  FUNCTION alpha_sw_vctr  *** 
     526      !! 
     527      !! ** Purpose : ROUGH estimate of the thermal expansion coefficient of sea-water at the surface (P =~ 1010 hpa) 
     528      !! 
     529      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
     530      !!---------------------------------------------------------------------------------- 
     531      REAL(wp), DIMENSION(jpi,jpj)             ::   alpha_sw_vctr   ! latent heat of vaporization   [J/kg] 
     532      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   psst   ! water temperature                [K] 
     533      !!---------------------------------------------------------------------------------- 
     534      alpha_sw_vctr = 2.1e-5_wp * MAX(psst(:,:)-rt0 + 3.2_wp, 0._wp)**0.79 
     535   END FUNCTION alpha_sw_vctr 
     536 
     537   FUNCTION alpha_sw_sclr( psst ) 
     538      !!--------------------------------------------------------------------------------- 
     539      !!                           ***  FUNCTION alpha_sw_sclr  *** 
     540      !! 
     541      !! ** Purpose : ROUGH estimate of the thermal expansion coefficient of sea-water at the surface (P =~ 1010 hpa) 
     542      !! 
     543      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
     544      !!---------------------------------------------------------------------------------- 
     545      REAL(wp)             ::   alpha_sw_sclr   ! latent heat of vaporization   [J/kg] 
     546      REAL(wp), INTENT(in) ::   psst   ! sea-water temperature                   [K] 
     547      !!---------------------------------------------------------------------------------- 
     548      alpha_sw_sclr = 2.1e-5_wp * MAX(psst-rt0 + 3.2_wp, 0._wp)**0.79 
     549   END FUNCTION alpha_sw_sclr 
     550 
    359551 
    360552 
Note: See TracChangeset for help on using the changeset viewer.