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 – NEMO

Changeset 11215


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.

Location:
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r11182 r11215  
    270270         WRITE(numout,*) '      "NCAR"      algorithm   (Large and Yeager 2008)     ln_NCAR      = ', ln_NCAR 
    271271         WRITE(numout,*) '      "COARE 3.0" algorithm   (Fairall et al. 2003)       ln_COARE_3p0 = ', ln_COARE_3p0 
    272          WRITE(numout,*) '      "COARE 3.5" algorithm   (Edson et al. 2013)         ln_COARE_3p5 = ', ln_COARE_3p0 
     272         WRITE(numout,*) '      "COARE 3.5" algorithm   (Edson et al. 2013)         ln_COARE_3p5 = ', ln_COARE_3p5 
    273273         WRITE(numout,*) '      "ECMWF"     algorithm   (IFS cycle 31)              ln_ECMWF     = ', ln_ECMWF 
    274274         WRITE(numout,*) '      add High freq.contribution to the stress module     ln_taudif    = ', ln_taudif 
    … …  
    472472      CASE( np_NCAR      )   ;   CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! NCAR-COREv2 
    473473         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    474       CASE( np_COARE_3p0 )   ;   CALL turb_coare   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! COARE v3.0 
    475          &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
     474 
     475      !LB: Skin!!!          
     476      CASE( np_COARE_3p0 ) 
     477         IF ( ln_skin ) THEN             
     478            CALL turb_coare   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! COARE v3.0 
     479               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 
     480               &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1)) 
     481            zst(:,:) = zst(:,:)*tmask(:,:,1) 
     482            zsq(:,:) = zsq(:,:)*tmask(:,:,1) 
     483         ELSE 
     484            IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_coare" WITHOUT CSWL optional arrays!!!'             
     485            CALL turb_coare   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! COARE v3.0 
     486               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
     487         END IF 
     488             
    476489      CASE( np_COARE_3p5 )   ;   CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! COARE v3.5 
    477490         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    … …  
    480493      CASE( np_ECMWF     ) 
    481494         IF ( ln_skin ) THEN             
    482             !IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_ecmwf" WITH CSWL optional arrays!!!' 
    483             !IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => BEFORE ZST(40:50,30) =', zst(40:50,30) 
    484495            CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! ECMWF          
    485496               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 
    486497               &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1)) 
    487             !LB: "zst" and "zsq" have been updated with skin temp. !!! (from bulk SST)... 
    488498            zst(:,:) = zst(:,:)*tmask(:,:,1) 
    489499            zsq(:,:) = zsq(:,:)*tmask(:,:,1) 
    490             !IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => AFTER ZST(40:50,30) =', zst(40:50,30) 
    491500         ELSE 
    492501            IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_ecmwf" WITHOUT CSWL optional arrays!!!' 
    493502            CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! ECMWF          
    494                &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
     503               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    495504         END IF 
    496505         !LB. 

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

    r11209 r11215  
    1313   !!                     (https://github.com/brodeau/aerobulk) 
    1414   !! 
    15    !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk) 
     15   !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk) 
    1616   !!---------------------------------------------------------------------- 
    1717   !! History :  4.0  !  2016-02  (L.Brodeau)   Original code 
    … …  
    2626   USE dom_oce         ! ocean space and time domain 
    2727   USE phycst          ! physical constants 
    28    USE sbc_oce         ! Surface boundary condition: ocean fields 
     28   USE iom             ! I/O manager library 
     29   USE lib_mpp         ! distribued memory computing library 
     30   USE in_out_manager  ! I/O manager 
     31   USE prtctl          ! Print control 
    2932   USE sbcwave, ONLY   :  cdn_wave ! wave module 
    3033#if defined key_si3 || defined key_cice 
    3134   USE sbc_ice         ! Surface boundary condition: ice fields 
    3235#endif 
    33    USE sbcblk_phy     !LB: all thermodynamics functions, rho_air, q_sat, etc... #LB 
    34    USE in_out_manager ! I/O manager 
    35    USE iom            ! I/O manager library 
    36    USE lib_mpp        ! distribued memory computing library 
    37    USE prtctl         ! Print control 
    38    USE lib_fortran    ! to use key_nosignedzero 
     36   USE lib_fortran     ! to use key_nosignedzero 
     37 
     38   USE sbc_oce         ! Surface boundary condition: ocean fields 
     39   USE sbcblk_phy      ! all thermodynamics functions, rho_air, q_sat, etc... !LB 
     40   USE sbcblk_skin     ! cool-skin/warm layer scheme (CSWL_ECMWF) !LB 
    3941 
    4042   IMPLICIT NONE 
    … …  
    5254CONTAINS 
    5355 
    54    SUBROUTINE turb_coare( zt, zu, sst, t_zt, ssq, q_zt, U_zu, & 
     56   SUBROUTINE turb_coare( zt, zu, T_s, t_zt, q_s, q_zt, U_zu, & 
    5557      &                   Cd, Ch, Ce, t_zu, q_zu, U_blk,      & 
    56       &                   Cdn, Chn, Cen                       ) 
     58      &                   Cdn, Chn, Cen,                      & 
     59      &                   Qsw, rad_lw, slp                   ) 
    5760      !!---------------------------------------------------------------------- 
    5861      !!                      ***  ROUTINE  turb_coare  *** 
    … …  
    6366      !!                Returns the effective bulk wind speed at 10m to be used in the bulk formulas 
    6467      !! 
     68      !!                Applies the cool-skin warm-layer correction of the SST to T_s 
     69      !!                if the net shortwave flux at the surface (Qsw), the downwelling longwave 
     70      !!                radiative fluxes at the surface (rad_lw), and the sea-leve pressure (slp) 
     71      !!                are provided as (optional) arguments! 
    6572      !! 
    6673      !! INPUT : 
    … …  
    6976      !!    *  zu   : height for wind speed (generally 10m)                   [m] 
    7077      !!    *  U_zu : scalar wind speed at 10m                                [m/s] 
    71       !!    *  sst  : SST                                                     [K] 
    7278      !!    *  t_zt : potential air temperature at zt                         [K] 
    73       !!    *  ssq  : specific humidity at saturation at SST                  [kg/kg] 
    7479      !!    *  q_zt : specific humidity of air at zt                          [kg/kg] 
    7580      !! 
     81      !! INPUT/OUTPUT: 
     82      !! ------------- 
     83      !!    *  T_s  : SST or skin temperature                                 [K] 
     84      !!    *  q_s  : SSQ aka saturation specific humidity at temp. T_s       [kg/kg] 
     85      !!              -> doesn't need to be given a value if skin temp computed (in case l_use_skin=True) 
     86      !!              -> MUST be given the correct value if not computing skint temp. (in case l_use_skin=False) 
     87      !! 
     88      !! OPTIONAL INPUT (will trigger l_use_skin=TRUE if present!): 
     89      !! --------------- 
     90      !!    *  Qsw    : net solar flux (after albedo) at the surface (>0)     [W/m^2] 
     91      !!    *  rad_lw : downwelling longwave radiation at the surface  (>0)   [W/m^2] 
     92      !!    *  slp    : sea-level pressure                                    [Pa] 
    7693      !! 
    7794      !! OUTPUT : 
    … …  
    85102      !! 
    86103      !! 
    87       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk) 
     104      !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    88105      !!---------------------------------------------------------------------------------- 
    89106      REAL(wp), INTENT(in   )                     ::   zt       ! height for t_zt and q_zt                    [m] 
    90107      REAL(wp), INTENT(in   )                     ::   zu       ! height for U_zu                             [m] 
    91       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   sst      ! sea surface temperature                [Kelvin] 
     108      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) ::   T_s      ! sea surface temperature                [Kelvin] 
    92109      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   t_zt     ! potential air temperature              [Kelvin] 
    93       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   ssq      ! sea surface specific humidity           [kg/kg] 
     110      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) ::   q_s      ! sea surface specific humidity           [kg/kg] 
    94111      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_zt     ! specific air humidity at zt             [kg/kg] 
    95112      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   U_zu     ! relative wind module at zu                [m/s] 
    … …  
    102119      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients 
    103120      ! 
     121      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   Qsw      !             [W/m^2] 
     122      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   rad_lw   !             [W/m^2] 
     123      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   slp      !             [Pa] 
     124      ! 
    104125      INTEGER :: j_itt 
    105126      LOGICAL :: l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U 
    … …  
    113134      REAL(wp), DIMENSION(jpi,jpj) ::   ztmp0, ztmp1, ztmp2 
    114135      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zeta_t        ! stability parameter at height zt 
    115       !!---------------------------------------------------------------------- 
    116       ! 
     136      ! 
     137      ! Cool skin: 
     138      LOGICAL :: l_use_skin = .FALSE. 
     139      REAL(wp), DIMENSION(jpi,jpj) :: zsst    ! to back up the initial bulk SST 
     140      !!---------------------------------------------------------------------------------- 
     141      ! 
     142      ! Cool skin ? 
     143      IF( PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp) ) THEN 
     144         l_use_skin = .TRUE. 
     145      END IF 
     146      IF (lwp) PRINT *, ' *** LOLO(sbcblk_algo_ecmwf.F90) => l_use_skin =', l_use_skin 
     147 
    117148      l_zt_equal_zu = .FALSE. 
    118       IF( ABS(zu - zt) < 0.01 ) l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision 
     149      IF( ABS(zu - zt) < 0.01_wp )  l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision 
    119150 
    120151      IF( .NOT. l_zt_equal_zu )  ALLOCATE( zeta_t(jpi,jpj) ) 
    121152 
     153      !! Initialization for cool skin: 
     154      IF( l_use_skin ) THEN 
     155         zsst   = T_s    ! save the bulk SST 
     156         T_s    = T_s - 0.25                      ! First guess of correction 
     157         q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s 
     158      END IF 
     159 
    122160      !! First guess of temperature and humidity at height zu: 
    123       t_zu = MAX( t_zt , 199.0_wp )   ! who knows what's given on masked-continental regions... 
     161      t_zu = MAX( t_zt ,  180._wp )   ! who knows what's given on masked-continental regions... 
    124162      q_zu = MAX( q_zt , 1.e-6_wp )   !               " 
    125163 
    126164      !! Pot. temp. difference (and we don't want it to be 0!) 
    127       dt_zu = t_zu - sst ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 
    128       dq_zu = q_zu - ssq ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 
     165      dt_zu = t_zu - T_s ;   dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 
     166      dq_zu = q_zu - q_s ;   dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 
    129167 
    130168      znu_a = visc_air(t_zu) ! Air viscosity (m^2/s) at zt given from temperature in (K) 
    131169 
    132       ztmp2 = 0.5*0.5  ! initial guess for wind gustiness contribution 
     170      ztmp2 = 0.5_wp*0.5_wp  ! initial guess for wind gustiness contribution 
    133171      U_blk = SQRT(U_zu*U_zu + ztmp2) 
    134172 
    135       ztmp2   = 10000.     ! optimization: ztmp2 == 1/z0 (with z0 first guess == 0.0001) 
     173      ztmp2   = 10000._wp     ! optimization: ztmp2 == 1/z0 (with z0 first guess == 0.0001) 
    136174      ztmp0   = LOG(zu*ztmp2) 
    137175      ztmp1   = LOG(10.*ztmp2) 
    138       u_star = 0.035*U_blk*ztmp1/ztmp0       ! (u* = 0.035*Un10) 
    139  
    140  
    141       z0     = alfa_charn(U_blk)*u_star*u_star/grav + 0.11*znu_a/u_star 
    142       z0     = MIN(ABS(z0), 0.001)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO 
    143       z0t    = 1. / ( 0.1*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) ) 
    144       z0t    = MIN(ABS(z0t), 0.001)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO 
     176      u_star = 0.035_wp*U_blk*ztmp1/ztmp0       ! (u* = 0.035*Un10) 
     177 
     178      z0     = alfa_charn(U_blk)*u_star*u_star/grav + 0.11_wp*znu_a/u_star 
     179      z0     = MIN(ABS(z0), 0.001_wp)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO 
     180      z0t    = 1._wp / ( 0.1_wp*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) ) 
     181      z0t    = MIN(ABS(z0t), 0.001_wp)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO 
    145182 
    146183      ztmp2  = vkarmn/ztmp0 
    … …  
    149186      ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd 
    150187 
    151       ztmp2 = Ri_bulk( zu, sst, t_zu, ssq, q_zu, U_blk ) ! Bulk Richardson Number (BRN) 
    152       !ztmp2  = grav*zu*(dt_zu + rctv0*t_zu*dq_zu)/(t_zu*U_blk*U_blk)  !! Bulk Richardson number ;       !Ribcu = -zu/(zi0*0.004*Beta0**3) !! Saturation Rib, zi0 = tropicalbound. layer depth 
     188      ztmp2 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN) 
    153189 
    154190      !! First estimate of zeta_u, depending on the stability, ie sign of BRN (ztmp2): 
    155191      ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 ) 
    156192      ztmp0 = ztmp0*ztmp2 
    157       zeta_u = (1.-ztmp1) * (ztmp0/(1.+ztmp2/(-zu/(zi0*0.004*Beta0**3)))) & !  BRN < 0 
    158          &  +     ztmp1   * (ztmp0*(1. + 27./9.*ztmp2/ztmp0))               !  BRN > 0 
    159       !#LOLO: should make sure that the "ztmp0" of "27./9.*ztmp2/ztmp0" is "ztmp0*ztmp2" and not "ztmp0==vkarmn*vkarmn/LOG(zt/z0t)/Cd" ! 
     193      zeta_u = (1._wp-ztmp1) * (ztmp0/(1._wp+ztmp2/(-zu/(zi0*0.004_wp*Beta0**3)))) & !  BRN < 0 
     194         &  +     ztmp1   * (ztmp0*(1._wp + 27._wp/9._wp*ztmp2/ztmp0))               !  BRN > 0 
     195      !#LB: should make sure that the "ztmp0" of "27./9.*ztmp2/ztmp0" is "ztmp0*ztmp2" and not "ztmp0==vkarmn*vkarmn/LOG(zt/z0t)/Cd" ! 
    160196 
    161197      !! First guess M-O stability dependent scaling params.(u*,t*,q*) to estimate z0 and z/L 
    … …  
    174210         t_zu = t_zt - t_star/vkarmn*ztmp1 
    175211         q_zu = q_zt - q_star/vkarmn*ztmp1 
    176          q_zu = (0.5 + SIGN(0.5_wp,q_zu))*q_zu !Makes it impossible to have negative humidity : 
     212         q_zu = (0.5_wp + SIGN(0.5_wp,q_zu))*q_zu !Makes it impossible to have negative humidity : 
    177213         ! 
    178          dt_zu = t_zu - sst  ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 
    179          dq_zu = q_zu - ssq  ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 
     214         dt_zu = t_zu - T_s  ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 
     215         dq_zu = q_zu - q_s  ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 
    180216      END IF 
    181217 
    … …  
    200236 
    201237         !! Roughness lengthes z0, z0t (z0q = z0t) : 
    202          z0   = ztmp2*ztmp1/grav + 0.11*znu_a/u_star ! Roughness length (eq.6) 
    203          ztmp1 = z0*u_star/znu_a                           ! Re_r: roughness Reynolds number 
     238         z0    = ztmp2*ztmp1/grav + 0.11_wp*znu_a/u_star ! Roughness length (eq.6) 
     239         ztmp1 = z0*u_star/znu_a                          ! Re_r: roughness Reynolds number 
    204240         z0t  = MIN( 1.1E-4_wp , 5.5E-5_wp*ztmp1**(-0.6_wp) ) ! Scalar roughness for both theta and q (eq.28) 
    205241 
    … …  
    227263            t_zu = t_zt - t_star/vkarmn*ztmp1 
    228264            q_zu = q_zt - q_star/vkarmn*ztmp1 
    229             dt_zu = t_zu - sst ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6), dt_zu ) 
    230             dq_zu = q_zu - ssq ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9), dq_zu ) 
    231265         END IF 
    232266 
    233       END DO 
    234       ! 
     267         !! SKIN related part 
     268         !! ----------------- 
     269         IF( l_use_skin ) THEN 
     270            !! compute transfer coefficients at zu : lolo: verifier... 
     271            ztmp0 = u_star/U_blk 
     272            Ch   = ztmp0*t_star/dt_zu 
     273            Ce   = ztmp0*q_star/dq_zu 
     274            ! Non-Solar heat flux to the ocean: 
     275            ztmp1 = U_blk*MAX(rho_air(t_zu, q_zu, slp), 1._wp)     ! rho*U10 
     276            ztmp2 = T_s*T_s 
     277            ztmp1 = ztmp1 * ( Ce*rLevap*(q_zu - q_s) + Ch*rCp_dry*(t_zu - T_s) ) & ! Total turb. heat flux 
     278               &       +    (rad_lw - emiss_w*stefan*ztmp2*ztmp2)                  ! Net longwave flux 
     279            !! Updating the values of the skin temperature T_s and q_s : 
     280            CALL CSWL_ECMWF( Qsw, ztmp1, u_star, zsst, T_s ) ! yes ECMWF, because more advanced than COARE (warm-layer added!) 
     281            q_s = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp)  ! 200 -> just to avoid numerics problem on masked regions if silly values are given 
     282         END IF 
     283 
     284         IF( (l_use_skin).OR.(.NOT. l_zt_equal_zu) ) THEN 
     285            dt_zu = t_zu - T_s ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 
     286            dq_zu = q_zu - q_s ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 
     287         END IF 
     288 
     289      END DO !DO j_itt = 1, nb_itt 
     290 
    235291      ! compute transfer coefficients at zu : 
    236292      ztmp0 = u_star/U_blk 
    … …  
    259315      !! Wind greater than 18 m/s :  alfa = 0.018 
    260316      !! 
    261       !! Author: L. Brodeau, june 2016 / AeroBulk  (https://github.com/brodeau/aerobulk/) 
     317      !! Author: L. Brodeau, June 2016 / AeroBulk  (https://github.com/brodeau/aerobulk/) 
    262318      !!------------------------------------------------------------------- 
    263319      REAL(wp), DIMENSION(jpi,jpj) :: alfa_charn 
    … …  
    297353      !!       form from Beljaars and Holtslag (1991) 
    298354      !! 
    299       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
     355      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    300356      !!---------------------------------------------------------------------------------- 
    301357      REAL(wp), DIMENSION(jpi,jpj) :: psi_m_coare 
    … …  
    349405      !! form from Beljaars and Holtslag (1991) 
    350406      !! 
    351       !! Author: L. Brodeau, june 2016 / AeroBulk 
     407      !! Author: L. Brodeau, June 2016 / AeroBulk 
    352408      !!         (https://github.com/brodeau/aerobulk/) 
    353409      !!---------------------------------------------------------------- 

  • 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.