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

Ignore:
Timestamp:
2019-07-04T12:34:28+02:00 (5 years ago)
Author:
laurent
Message:

LB: inclusion of "cool-skin/warm-layer" support (a la ECMWF) into COARE 3.0 algorithm.

File:
1 edited

Legend:

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

    r11209 r11215  
    1414   !!                     (https://github.com/brodeau/aerobulk) 
    1515   !! 
    16    !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk) 
     16   !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk) 
    1717   !!---------------------------------------------------------------------- 
    1818   !! History :  4.0  !  2016-02  (L.Brodeau)   Original code 
     
    3838 
    3939   USE sbc_oce         ! Surface boundary condition: ocean fields 
    40    USE sbcblk_phy      !LB: all thermodynamics functions, rho_air, q_sat, etc... #LB 
    41    USE sbcblk_skin     ! cool-skin/warm layer scheme (CSWL_ECMWF) 
     40   USE sbcblk_phy      ! all thermodynamics functions, rho_air, q_sat, etc... !LB 
     41   USE sbcblk_skin     ! cool-skin/warm layer scheme (CSWL_ECMWF) !LB 
    4242 
    4343   IMPLICIT NONE 
     
    5757   INTEGER , PARAMETER ::   nb_itt = 5        ! number of itterations 
    5858 
    59    !! Cool-Skin / Warm-Layer related parameters: 
    60    REAL(wp),    PARAMETER :: & 
    61       &  rdt0    = 3600.*1.5, & !: time step 
    62       &  rd0     = 3. ,      &  !: Depth scale [m], "d" in Eq.11 (Zeng & Beljaars 2005) 
    63       &  rNu0    = 0.5          !: Nu (exponent of temperature profile) Eq.11 
    64    !                            !: (Zeng & Beljaars 2005) !: set to 0.5 instead of 
    65    !                            !: 0.3 to respect a warming of +3 K in calm 
    66    !                            !: condition for the insolation peak of +1000W/m^2 
    67    INTEGER,    PARAMETER :: & 
    68       &  nb_itt_wl = 10         !: number of sub-itterations for solving the differential equation in warm-layer part 
    69    !                            !:  => use "nb_itt_wl = 1" for No itteration! => way cheaper !!! 
    70    !                            !:    => assumes balance between the last 2 terms of Eq.11 (lhs of eq.11 = 0) 
    71    !                            !:    => in that case no need for sub-itterations ! 
    72    !                            !:    => ACCEPTABLE IN MOST CONDITIONS ! (UNLESS: sunny + very calm/low-wind conditions) 
    73    !                            !:  => Otherwize use "nb_itt_wl = 10" 
    7459   !!---------------------------------------------------------------------- 
    7560CONTAINS 
     
    123108      !! 
    124109      !! 
    125       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
     110      !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    126111      !!---------------------------------------------------------------------------------- 
    127112      REAL(wp), INTENT(in   )                     ::   zt       ! height for t_zt and q_zt                    [m] 
     
    157142      LOGICAL :: l_use_skin = .FALSE. 
    158143      REAL(wp), DIMENSION(jpi,jpj) :: zsst    ! to back up the initial bulk SST 
    159  
    160144      ! 
    161145      REAL(wp), DIMENSION(jpi,jpj) ::   func_m, func_h 
    162146      REAL(wp), DIMENSION(jpi,jpj) ::   ztmp0, ztmp1, ztmp2 
    163147      !!---------------------------------------------------------------------------------- 
    164       ! 
     148 
    165149      ! Cool skin ? 
    166150      IF( PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp) ) THEN 
     
    169153      IF (lwp) PRINT *, ' *** LOLO(sbcblk_algo_ecmwf.F90) => l_use_skin =', l_use_skin 
    170154 
    171       ! Identical first gess as in COARE, with IFS parameter values though 
     155      ! Identical first gess as in COARE, with IFS parameter values though... 
    172156      ! 
    173157      l_zt_equal_zu = .FALSE. 
     
    182166 
    183167      !! First guess of temperature and humidity at height zu: 
    184       t_zu = MAX( t_zt , 0.0_wp )   ! who knows what's given on masked-continental regions... 
    185       q_zu = MAX( q_zt , 1.e-6_wp)   !               " 
     168      t_zu = MAX( t_zt ,  180._wp )   ! who knows what's given on masked-continental regions... 
     169      q_zu = MAX( q_zt , 1.e-6_wp )   !               " 
    186170 
    187171      !! Pot. temp. difference (and we don't want it to be 0!) 
     
    189173      dq_zu = q_zu - q_s ;   dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 
    190174 
    191       znu_a = visc_air(t_zt) ! Air viscosity (m^2/s) at zt given from temperature in (K) 
     175      znu_a = visc_air(t_zu) ! Air viscosity (m^2/s) at zt given from temperature in (K) 
    192176 
    193177      ztmp2 = 0.5_wp*0.5_wp  ! initial guess for wind gustiness contribution 
     
    214198      ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 ) 
    215199      func_m = ztmp0*ztmp2 ! temporary array !! 
    216       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 
     200      func_h = (1._wp-ztmp1) * (func_m/(1._wp+ztmp2/(-zu/(zi0*0.004_wp*Beta0**3)))) & !  BRN < 0 ! temporary array !!! func_h == zeta_u 
     201         &  +     ztmp1   * (func_m*(1._wp + 27._wp/9._wp*ztmp2/func_m))              !  BRN > 0 
    218202      !#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" ! 
    219203 
     
    259243         Linv = ztmp0*func_m*func_m/func_h / zu     ! From Eq. 3.23, Chap.3.2.3, IFS doc - Cy40r1 
    260244         !! Note: it is slightly different that the L we would get with the usual 
    261          !! expression, as in coare algorithm or in 'mod_phy.f90' (One_on_L_MO()) 
    262245         Linv = SIGN( MIN(ABS(Linv),200._wp), Linv ) ! (prevent FPE from stupid values from masked region later on...) !#LOLO 
    263246 
     
    315298         IF( l_use_skin ) THEN 
    316299            !! compute transfer coefficients at zu : lolo: verifier... 
    317             Cd = vkarmn*vkarmn/(func_m*func_m) 
    318300            Ch = vkarmn*vkarmn/(func_m*func_h) 
    319301            ztmp1 = LOG(zu) - LOG(z0q) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0q*Linv)   ! func_q 
     
    357339      !!         and L is M-O length 
    358340      !! 
    359       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
     341      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    360342      !!---------------------------------------------------------------------------------- 
    361343      REAL(wp), DIMENSION(jpi,jpj) :: psi_m_ecmwf 
     
    405387      !!         and L is M-O length 
    406388      !! 
    407       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
     389      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    408390      !!---------------------------------------------------------------------------------- 
    409391      REAL(wp), DIMENSION(jpi,jpj) :: psi_h_ecmwf 
Note: See TracChangeset for help on using the changeset viewer.