Changeset 11284


Ignore:
Timestamp:
2019-07-18T11:35:28+02:00 (13 months ago)
Author:
laurent
Message:

LB: gooodbye "coare3p5", hello "coare3p6" + cleaner air humidity handling + support for Relative humidity.

Location:
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC
Files:
2 added
2 deleted
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90

    r11266 r11284  
    4646#endif 
    4747   USE sbcblk_algo_ncar     ! => turb_ncar     : NCAR - CORE (Large & Yeager, 2009) 
    48    USE sbcblk_algo_coare    ! => turb_coare    : COAREv3.0 (Fairall et al. 2003) 
    49    USE sbcblk_algo_coare3p5 ! => turb_coare3p5 : COAREv3.5 (Edson et al. 2013) 
     48   USE sbcblk_algo_coare3p0 ! => turb_coare3p0 : COAREv3.0 (Fairall et al. 2003) 
     49   USE sbcblk_algo_coare3p6 ! => turb_coare3p6 : COAREv3.6 (Fairall et al. 2018 + Edson et al. 2013) 
    5050   USE sbcblk_algo_ecmwf    ! => turb_ecmwf    : ECMWF (IFS cycle 31) 
    5151   ! 
     
    8787   LOGICAL  ::   ln_NCAR        ! "NCAR"      algorithm   (Large and Yeager 2008) 
    8888   LOGICAL  ::   ln_COARE_3p0   ! "COARE 3.0" algorithm   (Fairall et al. 2003) 
    89    LOGICAL  ::   ln_COARE_3p5   ! "COARE 3.5" algorithm   (Edson et al. 2013) 
     89   LOGICAL  ::   ln_COARE_3p6   ! "COARE 3.6" algorithm   (Edson et al. 2013) 
    9090   LOGICAL  ::   ln_ECMWF       ! "ECMWF"     algorithm   (IFS cycle 31) 
    9191   ! 
     
    109109   !LB: 
    110110   LOGICAL  ::   ln_skin        ! use the cool-skin / warm-layer parameterization (only available in ECMWF and COARE algorithms) !LB 
    111    LOGICAL  ::   ln_humi_dpt    ! humidity read in files ("sn_humi") is dew-point temperature if .true. (specific humidity espected if .false.) !LB 
     111   LOGICAL  ::   ln_humi_sph    ! humidity read in files ("sn_humi") is specific humidity [kg/kg] if .true. !LB 
     112   LOGICAL  ::   ln_humi_dpt    ! humidity read in files ("sn_humi") is dew-point temperature [K] if .true. !LB 
     113   LOGICAL  ::   ln_humi_rlh    ! humidity read in files ("sn_humi") is relative humidity     [%] if .true. !LB 
     114   ! 
     115   INTEGER  ::   nhumi          ! choice of the bulk algorithm 
     116   !                            ! associated indices: 
     117   INTEGER, PARAMETER :: np_humi_sph = 1 
     118   INTEGER, PARAMETER :: np_humi_dpt = 2 
     119   INTEGER, PARAMETER :: np_humi_rlh = 3 
    112120   !LB. 
    113121 
     
    116124   INTEGER, PARAMETER ::   np_NCAR      = 1   ! "NCAR" algorithm        (Large and Yeager 2008) 
    117125   INTEGER, PARAMETER ::   np_COARE_3p0 = 2   ! "COARE 3.0" algorithm   (Fairall et al. 2003) 
    118    INTEGER, PARAMETER ::   np_COARE_3p5 = 3   ! "COARE 3.5" algorithm   (Edson et al. 2013) 
     126   INTEGER, PARAMETER ::   np_COARE_3p6 = 3   ! "COARE 3.6" algorithm   (Edson et al. 2013) 
    119127   INTEGER, PARAMETER ::   np_ECMWF     = 4   ! "ECMWF" algorithm       (IFS cycle 31) 
    120128 
     
    172180      NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw ,                &   ! input fields 
    173181         &                 sn_tair, sn_prec, sn_snow, sn_slp, sn_tdif,                & 
    174          &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF,             &   ! bulk algorithm 
     182         &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p6, ln_ECMWF,             &   ! bulk algorithm 
    175183         &                 cn_dir , ln_taudif, rn_zqt, rn_zu,                         & 
    176184         &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15,           & 
    177          &                 ln_skin, ln_humi_dpt                                                        ! cool-skin / warm-layer param. !LB 
     185         &                 ln_skin, ln_humi_sph, ln_humi_dpt, ln_humi_rlh    ! cool-skin / warm-layer !LB 
    178186      !!--------------------------------------------------------------------- 
    179187      ! 
     
    201209         nblk =  np_COARE_3p0   ;   ioptio = ioptio + 1 
    202210      ENDIF 
    203       IF( ln_COARE_3p5 ) THEN 
    204          nblk =  np_COARE_3p5   ;   ioptio = ioptio + 1 
     211      IF( ln_COARE_3p6 ) THEN 
     212         nblk =  np_COARE_3p6   ;   ioptio = ioptio + 1 
    205213      ENDIF 
    206214      IF( ln_ECMWF     ) THEN 
    207215         nblk =  np_ECMWF       ;   ioptio = ioptio + 1 
    208216      ENDIF 
    209       ! 
     217      IF( ioptio /= 1 )   CALL ctl_stop( 'sbc_blk_init: Choose one and only one bulk algorithm' ) 
     218 
     219 
     220 
    210221      !LB: 
    211222      !                             !** initialization of the cool-skin / warm-layer parametrization 
     
    215226         IF( sbc_blk_cswl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard arrays' ) 
    216227      END IF 
     228      ! 
     229      ioptio = 0 
     230      IF( ln_humi_sph ) THEN 
     231         nhumi =  np_humi_sph    ;   ioptio = ioptio + 1 
     232      ENDIF 
     233      IF( ln_humi_dpt ) THEN 
     234         nhumi =  np_humi_dpt    ;   ioptio = ioptio + 1 
     235      ENDIF 
     236      IF( ln_humi_rlh ) THEN 
     237         nhumi =  np_humi_rlh    ;   ioptio = ioptio + 1 
     238      ENDIF 
     239      IF( ioptio /= 1 )   CALL ctl_stop( 'sbc_blk_init: Choose one and only one type of air humidity' ) 
    217240      !LB. 
    218241 
    219       IF( ioptio /= 1 )   CALL ctl_stop( 'sbc_blk_init: Choose one and only one bulk algorithm' ) 
    220242      ! 
    221243      IF( ln_dm2dc ) THEN                 !* check: diurnal cycle on Qsr 
     
    278300         WRITE(numout,*) '      "NCAR"      algorithm   (Large and Yeager 2008)     ln_NCAR      = ', ln_NCAR 
    279301         WRITE(numout,*) '      "COARE 3.0" algorithm   (Fairall et al. 2003)       ln_COARE_3p0 = ', ln_COARE_3p0 
    280          WRITE(numout,*) '      "COARE 3.5" algorithm   (Edson et al. 2013)         ln_COARE_3p5 = ', ln_COARE_3p5 
     302         WRITE(numout,*) '      "COARE 3.6" algorithm (Fairall 2018 + Edson al 2013)ln_COARE_3p6 = ', ln_COARE_3p6 
    281303         WRITE(numout,*) '      "ECMWF"     algorithm   (IFS cycle 31)              ln_ECMWF     = ', ln_ECMWF 
    282304         WRITE(numout,*) '      add High freq.contribution to the stress module     ln_taudif    = ', ln_taudif 
     
    294316         CASE( np_NCAR      )   ;   WRITE(numout,*) '   ==>>>   "NCAR" algorithm        (Large and Yeager 2008)' 
    295317         CASE( np_COARE_3p0 )   ;   WRITE(numout,*) '   ==>>>   "COARE 3.0" algorithm   (Fairall et al. 2003)' 
    296          CASE( np_COARE_3p5 )   ;   WRITE(numout,*) '   ==>>>   "COARE 3.5" algorithm   (Edson et al. 2013)' 
     318         CASE( np_COARE_3p6 )   ;   WRITE(numout,*) '   ==>>>   "COARE 3.6" algorithm (Fairall 2018+Edson et al. 2013)' 
    297319         CASE( np_ECMWF     )   ;   WRITE(numout,*) '   ==>>>   "ECMWF" algorithm       (IFS cycle 31)' 
    298320         END SELECT 
     
    300322         WRITE(numout,*) 
    301323         WRITE(numout,*) '      use cool-skin/warm-layer parameterization (SSST)     ln_skin     = ', ln_skin !LB 
     324         ! 
     325         !LB: 
    302326         WRITE(numout,*) 
    303          WRITE(numout,*) '      air humidity read in files is dew-point temperature? ln_humi_dpt = ', ln_humi_dpt !LB 
     327         SELECT CASE( nhumi )              !* Print the choice of air humidity 
     328         CASE( np_humi_sph )   ;   WRITE(numout,*) '   ==>>>   air humidity is SPECIFIC HUMIDITY     [kg/kg]' 
     329         CASE( np_humi_dpt )   ;   WRITE(numout,*) '   ==>>>   air humidity is DEW-POINT TEMPERATURE [K]' 
     330         CASE( np_humi_rlh )   ;   WRITE(numout,*) '   ==>>>   air humidity is RELATIVE HUMIDITY     [%]' 
     331         END SELECT 
     332         !LB. 
    304333         ! 
    305334      ENDIF 
     
    360389         tatm_ice(:,:)    = sf(jp_tair)%fnow(:,:,1) 
    361390         !LB: 
    362          IF ( ln_humi_dpt ) THEN 
    363             qatm_ice(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 
    364          ELSE 
    365             qatm_ice(:,:) = sf(jp_humi)%fnow(:,:,1) 
    366          END IF 
     391         SELECT CASE( nhumi ) 
     392         CASE( np_humi_sph ) 
     393            qatm_ice(:,:) =           sf(jp_humi)%fnow(:,:,1) 
     394         CASE( np_humi_dpt ) 
     395            qatm_ice(:,:) = q_sat(    sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 
     396         CASE( np_humi_rlh ) 
     397            qatm_ice(:,:) = q_air_rh( 0.01_wp*sf(jp_humi)%fnow(:,:,1), sf(jp_tair)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1)) !RH is % percent in file 
     398         END SELECT 
    367399         !LB. 
    368400         tprecip(:,:)     = sf(jp_prec)%fnow(:,:,1) * rn_pfac 
     
    471503 
    472504      !LB: 
    473       ! spec. hum. of air at "rn_zqt" m above thes sea: 
    474       IF ( ln_humi_dpt ) THEN 
    475          ! spec. hum. of air is deduced from dew-point temperature and SLP: q_air = q_sat( d_air, slp) 
     505      ! zqair = specific humidity of air at "rn_zqt" m above thes sea: 
     506      SELECT CASE( nhumi ) 
     507      CASE( np_humi_sph ) 
     508         zqair(:,:) =        sf(jp_humi)%fnow(:,:,1)      ! what we read in file is already a spec. humidity! 
     509      CASE( np_humi_dpt ) 
    476510         IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => computing q_air out of d_air and slp !' 
    477511         zqair(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 
    478       ELSE 
    479          zqair(:,:) = sf(jp_humi)%fnow(:,:,1) ! what we read in file is already a spec. humidity! 
    480       END IF 
     512      CASE( np_humi_rlh ) 
     513         IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => computing q_air out of RH, t_air and slp !' 
     514         zqair(:,:) = q_air_rh( 0.01_wp*sf(jp_humi)%fnow(:,:,1), sf(jp_tair)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 
     515      END SELECT 
    481516      !LB. 
    482517 
     
    485520      !!    (since reanalysis products provide T at z, not theta !) 
    486521      ztpot = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), zqair(:,:) ) * rn_zqt 
    487        
     522 
    488523 
    489524      IF ( ln_skin ) THEN 
    490           
     525 
    491526         SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point 
    492527 
    493528         CASE( np_COARE_3p0 ) 
    494             CALL turb_coare   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! COARE v3.0 
     529            CALL turb_coare3p0 ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! COARE v3.0 
     530               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 
     531               &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1), & 
     532               &                Tsk_b=tsk(:,:) ) 
     533 
     534         CASE( np_COARE_3p6 ) 
     535            CALL turb_coare3p6 ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! COARE v3.6 
    495536               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 
    496537               &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1), & 
     
    502543               &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1), & 
    503544               &                Tsk_b=tsk(:,:) ) 
    504              
     545 
    505546         CASE DEFAULT 
    506547            CALL ctl_stop( 'STOP', 'sbc_oce: unsuported bulk formula selection for "ln_skin==.true."' ) 
     
    523564      ELSE !IF ( ln_skin ) 
    524565 
    525           
     566 
    526567         SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point 
    527568            ! 
     
    530571 
    531572         CASE( np_COARE_3p0 ) 
    532             IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_coare" WITHOUT CSWL optional arrays!!!' 
    533             CALL turb_coare   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! COARE v3.0 
     573            IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_coare3p0" WITHOUT CSWL optional arrays!!!' 
     574            CALL turb_coare3p0   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! COARE v3.0 
    534575               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    535576 
    536          CASE( np_COARE_3p5 )   ;   CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! COARE v3.5 
     577         CASE( np_COARE_3p6 )   ;   CALL turb_coare3p6( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! COARE v3.6 
    537578            &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    538579 
     
    802843 
    803844      !LB: 
    804       ! spec. hum. of air at "rn_zqt" m above thes sea: 
    805       IF ( ln_humi_dpt ) THEN 
    806          ! spec. hum. of air is deduced from dew-point temperature and SLP: q_air = q_sat( d_air, slp) 
     845      SELECT CASE( nhumi ) 
     846      CASE( np_humi_sph ) 
     847         zqair(:,:) =        sf(jp_humi)%fnow(:,:,1)      ! what we read in file is already a spec. humidity! 
     848      CASE( np_humi_dpt ) 
    807849         IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => ICE !!! computing q_air out of d_air and slp !' 
    808850         zqair(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 
    809       ELSE 
    810          zqair(:,:) = sf(jp_humi)%fnow(:,:,1) ! what we read in file is already a spec. humidity! 
    811       END IF 
     851      CASE( np_humi_rlh ) 
     852         IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => ICE !!! computing q_air out of RH, t_air and slp !' 
     853         zqair(:,:) = q_air_rh( 0.01_wp*sf(jp_humi)%fnow(:,:,1), sf(jp_tair)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 
     854      END SELECT 
    812855      !LB. 
    813856 
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_phy.F90

    r11266 r11284  
    1717   !!   Ri_bulk       : bulk Richardson number aka BRN 
    1818   !!   q_sat         : saturation humidity as a function of SLP and temperature 
    19  
     19   !!   q_air_rh      : specific humidity as a function of RH, t_air and SLP 
     20    
    2021   USE dom_oce        ! ocean space and time domain 
    2122   USE phycst         ! physical constants 
     
    3738   PUBLIC Ri_bulk 
    3839   PUBLIC q_sat 
     40   PUBLIC q_air_rh 
    3941 
    4042   !!---------------------------------------------------------------------- 
     
    115117      !!                           ***  FUNCTION L_vap  *** 
    116118      !! 
    117       !! ** Purpose : Compute the latent heat of vaporization of water from temperature 
     119      !! ** Purpose : Compute the latent heat of vaporization of water out of temperature 
    118120      !! 
    119121      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
     
    273275 
    274276 
     277   FUNCTION e_sat_sclr( ptak, pslp ) 
     278      !!---------------------------------------------------------------------------------- 
     279      !!                   ***  FUNCTION e_sat_sclr  *** 
     280      !!                  < SCALAR argument version > 
     281      !! ** Purpose : water vapor at saturation in [Pa] 
     282      !!              Based on accurate estimate by Goff, 1957 
     283      !! 
     284      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
     285      !! 
     286      !!    Note: what rt0 should be here, is 273.16 (triple point of water) and not 273.15 like here 
     287      !!---------------------------------------------------------------------------------- 
     288      REAL(wp), INTENT(in) ::   ptak    ! air temperature                  [K] 
     289      REAL(wp), INTENT(in) ::   pslp    ! sea level atmospheric pressure   [Pa] 
     290      REAL(wp)             ::   e_sat_sclr   ! water vapor at saturation   [kg/kg] 
     291      ! 
     292      REAL(wp) ::   zta, ztmp   ! local scalar 
     293      !!---------------------------------------------------------------------------------- 
     294      ! 
     295      zta = MAX( ptak , 180._wp )   ! air temp., prevents fpe0 errors dute to unrealistically low values over masked regions... 
     296      ztmp = rt0 / zta 
     297      ! 
     298      ! Vapour pressure at saturation [Pa] : WMO, (Goff, 1957) 
     299      e_sat_sclr = 100.*( 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(zta/rt0)        & 
     300         &    + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(zta/rt0 - 1.)) )  & 
     301         &    + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614) ) 
     302      ! 
     303   END FUNCTION e_sat_sclr 
     304 
     305 
    275306   FUNCTION q_sat( ptak, pslp ) 
    276307      !!---------------------------------------------------------------------------------- 
     
    288319      ! 
    289320      INTEGER  ::   ji, jj         ! dummy loop indices 
    290       REAL(wp) ::   zta, ze_sat, ztmp   ! local scalar 
    291       !!---------------------------------------------------------------------------------- 
    292       ! 
    293       DO jj = 1, jpj 
    294          DO ji = 1, jpi 
    295             ! 
    296             zta = MAX( ptak(ji,jj) , 180._wp )   ! air temp., prevents fpe0 errors dute to unrealistically low values over masked regions... 
    297             ztmp = rt0 / zta 
    298             ! 
    299             ! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957) 
    300             ze_sat = 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(zta/rt0)        & 
    301                &    + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(zta/rt0 - 1.)) )  & 
    302                &    + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614  ) 
    303             ! 
    304             q_sat(ji,jj) = reps0 * ze_sat/( 0.01_wp*pslp(ji,jj) - (1._wp - reps0)*ze_sat )   ! 0.01 because SLP is in [Pa] 
     321      REAL(wp) ::   ze_sat   ! local scalar 
     322      !!---------------------------------------------------------------------------------- 
     323      ! 
     324      DO jj = 1, jpj 
     325         DO ji = 1, jpi 
     326            ! 
     327            ze_sat =  e_sat_sclr( ptak(ji,jj), pslp(ji,jj) ) 
     328            ! 
     329            q_sat(ji,jj) = reps0 * ze_sat/( pslp(ji,jj) - (1._wp - reps0)*ze_sat ) 
    305330            ! 
    306331         END DO 
     
    308333      ! 
    309334   END FUNCTION q_sat 
     335 
     336   FUNCTION q_air_rh(prha, ptak, pslp) 
     337      !!---------------------------------------------------------------------------------- 
     338      !! Specific humidity of air out of Relative Humidity 
     339      !! 
     340      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
     341      !!---------------------------------------------------------------------------------- 
     342      REAL(wp), DIMENSION(jpi,jpj)             :: q_air_rh 
     343      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: prha        !: relative humidity      [fraction, not %!!!] 
     344      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak        !: air temperature        [K] 
     345      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pslp        !: atmospheric pressure   [Pa] 
     346      ! 
     347      INTEGER  ::   ji, jj      ! dummy loop indices 
     348      REAL(wp) ::   ze      ! local scalar 
     349      !!---------------------------------------------------------------------------------- 
     350      ! 
     351      DO jj = 1, jpj 
     352         DO ji = 1, jpi 
     353            ze = prha(ji,jj)*e_sat_sclr(ptak(ji,jj), pslp(ji,jj)) 
     354            q_air_rh(ji,jj) = ze*reps0/(pslp(ji,jj) - (1. - reps0)*ze) 
     355         END DO 
     356      END DO 
     357      ! 
     358   END FUNCTION q_air_rh 
    310359 
    311360 
Note: See TracChangeset for help on using the changeset viewer.