source: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_coare.F90 @ 11962

MODULE sbcblk_skin_coare
   !!======================================================================
   !!                   ***  MODULE  sbcblk_skin_coare  ***
   !! Computes:
   !!   * the surface skin temperature (aka SSST) based on the cool-skin/warm-layer
   !!     scheme used at ECMWF
   !!    Using formulation/param. of COARE 3.6 (Fairall et al., 2019)
   !!
   !!   From routine turb_ecmwf maintained and developed in AeroBulk
   !!            (https://github.com/brodeau/aerobulk)
   !!
   !! ** Author: L. Brodeau, September 2019 / AeroBulk (https://github.com/brodeau/aerobulk)
   !!----------------------------------------------------------------------
   !! History :  4.0  !  2016-02  (L.Brodeau)   Original code
   !!----------------------------------------------------------------------

   !!----------------------------------------------------------------------
   !!  cswl_ecmwf : computes the surface skin temperature (aka SSST)
   !!               based on the cool-skin/warm-layer scheme used at ECMWF
   !!----------------------------------------------------------------------
   USE oce             ! ocean dynamics and tracers
   USE dom_oce         ! ocean space and time domain
   USE phycst          ! physical constants
   USE sbc_oce         ! Surface boundary condition: ocean fields

   USE sbcblk_phy      !LOLO: all thermodynamics functions, rho_air, q_sat, etc...

   USE sbcdcy          !LOLO: to know hour of dawn and dusk: rdawn_dcy and rdusk_dcy (needed in WL_COARE)
   
   USE lib_mpp         ! distribued memory computing library
   USE in_out_manager  ! I/O manager
   USE lib_fortran     ! to use key_nosignedzero

   IMPLICIT NONE
   PRIVATE

   PUBLIC :: CS_COARE, WL_COARE

   !! Cool-skin related parameters:
   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: &
      &                        dT_cs         !: dT due to cool-skin effect => temperature difference between air-sea interface (z=0) and right below viscous layer (z=delta)

   !! Warm-layer related parameters:
   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: &
      &                        dT_wl,     &  !: dT due to warm-layer effect => difference between "almost surface (right below viscous layer, z=delta) and depth of bulk SST (z=gdept_1d(1))
      &                        Hz_wl,     &  !: depth of warm-layer [m]
      &                        Qnt_ac,    &  !: time integral / accumulated heat stored by the warm layer Qxdt => [J/m^2] (reset to zero every midnight)
      &                        Tau_ac        !: time integral / accumulated momentum Tauxdt => [N.s/m^2] (reset to zero every midnight)

   REAL(wp), PARAMETER, PUBLIC :: Hwl_max = 20._wp    !: maximum depth of warm layer (adjustable)
   !
   REAL(wp), PARAMETER :: rich   = 0.65_wp   !: critical Richardson number
   !
   REAL(wp), PARAMETER :: zfr0   = 0.5_wp    !: initial value of solar flux absorption
   !
   !!----------------------------------------------------------------------
CONTAINS


   SUBROUTINE CS_COARE( pQsw, pQnsol, pustar, pSST, pQlat )
      !!---------------------------------------------------------------------
      !!
      !! Cool-skin parameterization, based on Fairall et al., 1996,
      !! revisited for COARE 3.6 (Fairall et al., 2019)
      !!
      !! Fairall, C. W., Bradley, E. F., Godfrey, J. S., Wick, G. A.,
      !! Edson, J. B., and Young, G. S. ( 1996), Cool‐skin and warm‐layer
      !! effects on sea surface temperature, J. Geophys. Res., 101( C1), 1295-1308,
      !! doi:10.1029/95JC03190.
      !!
      !!------------------------------------------------------------------
      !!
      !!  **   INPUT:
      !!     *pQsw*       surface net solar radiation into the ocean     [W/m^2] => >= 0 !
      !!     *pQnsol*     surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 !
      !!     *pustar*     friction velocity u*                           [m/s]
      !!     *pSST*       bulk SST (taken at depth gdept_1d(1))          [K]
      !!     *pQlat*      surface latent heat flux                       [K]
      !!------------------------------------------------------------------
      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQsw   ! net solar a.k.a shortwave radiation into the ocean (after albedo) [W/m^2]
      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQnsol ! non-solar heat flux to the ocean [W/m^2]
      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pustar  ! friction velocity, temperature and humidity (u*,t*,q*)
      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pSST ! bulk SST [K]
      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQlat  ! latent heat flux [W/m^2]
      !!---------------------------------------------------------------------
      INTEGER  :: ji, jj, jc
      REAL(wp) :: zQabs, zdelta, zfr
      !!---------------------------------------------------------------------
      DO jj = 1, jpj
         DO ji = 1, jpi

            zQabs = pQnsol(ji,jj) ! first guess of heat flux absorbed within the viscous sublayer of thicknes delta,
            !                     !   => we DO not miss a lot assuming 0 solar flux absorbed in the tiny layer of thicknes zdelta...

            zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pQlat(ji,jj), pustar(ji,jj) )

            DO jc = 1, 4 ! because implicit in terms of zdelta...
               zfr = MAX( 0.137_wp + 11._wp*zdelta - 6.6E-5_wp/zdelta*(1._wp - EXP(-zdelta/8.E-4_wp)) , 0.01_wp ) ! Solar absorption, Eq.16 (Fairall al. 1996b) /  !LB: why 0.065 and not 0.137 like in the paper??? Beljaars & Zeng use 0.065, not 0.137 !
               zQabs = pQnsol(ji,jj) + zfr*pQsw(ji,jj)
               zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pQlat(ji,jj), pustar(ji,jj) )
            END DO

            dT_cs(ji,jj) = zQabs*zdelta/rk0_w   ! temperature increment, yes dT_cs can actually > 0, if Qabs > 0 (rare but possible!)

         END DO
      END DO

   END SUBROUTINE CS_COARE


   SUBROUTINE WL_COARE( pQsw, pQnsol, pTau, pSST, iwait )
      !!---------------------------------------------------------------------
      !!
      !!  Warm-Layer scheme according to COARE 3.6 (Fairall et al, 2019)
      !!     ------------------------------------------------------------------
      !!
      !!  **   INPUT:
      !!     *pQsw*       surface net solar radiation into the ocean     [W/m^2] => >= 0 !
      !!     *pQnsol*     surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 !
      !!     *pTau*       surface wind stress                            [N/m^2]
      !!     *pSST*       bulk SST  (taken at depth gdept_1d(1))         [K]
      !!     *iwait*      if /= 0 then wait before updating accumulated fluxes, we are within a converging itteration loop...
      !!---------------------------------------------------------------------
      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pQsw     ! surface net solar radiation into the ocean [W/m^2]     => >= 0 !
      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pQnsol   ! surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 !
      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTau     ! wind stress [N/m^2]
      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pSST     ! bulk SST at depth gdept_1d(1) [K]
      INTEGER ,                     INTENT(in)  :: iwait    ! if /= 0 then wait before updating accumulated fluxes
      !!
      INTEGER :: ji,jj
      !
      REAL(wp) :: zdTwl, zHwl, zQabs, zfr
      REAL(wp) :: zqac, ztac
      REAL(wp) :: zalpha, zcd1, zcd2, flg
      !!---------------------------------------------------------------------

      REAL(wp) :: ztime, znoon, zmidn
      INTEGER  :: jl

      LOGICAL :: l_exit, l_destroy_wl

      !! INITIALIZATION:
      zQabs  = 0._wp       ! total heat flux absorped in warm layer
      zfr    = zfr0        ! initial value of solar flux absorption !LOLO: save it and use previous value !!!

      IF( .NOT. ln_dm2dc ) CALL sbc_dcy_param() ! we need to call sbc_dcy_param (sbcdcy.F90) because rdawn_dcy and rdusk_dcy are unkonwn otherwize!

      ztime = REAL(nsec_day,wp)/(24._wp*3600._wp) ! time of current time step since 00:00 for current day (UTC) -> ztime = 0 -> 00:00 / ztime = 0.5 -> 12:00 ...

      IF (lwp) THEN
         WRITE(numout,*) ''
         WRITE(numout,*) 'LOLO: sbcblk_skin_coare => nsec_day, ztime =', nsec_day, ztime
      END IF
      
      DO jj = 1, jpj
         DO ji = 1, jpi

            l_exit       = .FALSE.
            l_destroy_wl = .FALSE.

            zdTwl =  dT_wl(ji,jj)                          ! value of previous time step as first guess
            zHwl  = MAX( MIN(Hz_wl(ji,jj),Hwl_max),0.1_wp) !   "                  "           "

            zqac = Qnt_ac(ji,jj) ! previous time step Qnt_ac
            ztac = Tau_ac(ji,jj)

            !*****  variables for warm layer  ***
            zalpha = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!)

            zcd1 = SQRT(2._wp*rich*rCp0_w/(zalpha*grav*rho0_w))        !mess-o-constants 1
            zcd2 = SQRT(2._wp*zalpha*grav/(rich*rho0_w))/(rCp0_w**1.5) !mess-o-constants 2

            
            znoon = MOD( 0.5_wp*(rdawn_dcy(ji,jj)+rdusk_dcy(ji,jj)), 1._wp )   ! 0<rnoon<1. => rnoon*24 = UTC time of local noon
            zmidn = MOD(              znoon-0.5_wp                 , 1._wp )
            zmidn = MOD( zmidn + 0.125_wp , 1._wp ) ! 3 hours past the local midnight 

            IF ( (ztime >= zmidn) .AND. (ztime < rdawn_dcy(ji,jj)) ) THEN
               ! Dawn reset to 0!
               l_exit       = .TRUE.
               l_destroy_wl = .TRUE.
            END IF

            IF ( .NOT. l_exit ) THEN
               !! Initial test on initial guess of absorbed heat flux in warm-layer:
               zfr = 1._wp - ( 0.28*0.014*(1. - EXP(-zHwl/0.014)) + 0.27*0.357*(1. - EXP(-zHwl/0.357)) &
                  &        + 0.45*12.82*(1-EXP(-zHwl/12.82)) ) / zHwl
               zQabs = zfr*pQsw(ji,jj) + pQnsol(ji,jj) ! first guess of tot. heat flux absorbed in warm layer !LOLO: depends of zfr, which is wild guess... Wrong!!!

               IF ( (ABS(zdTwl) < 1.E-6_wp) .AND. (zQabs <= 0._wp) ) THEN
                  ! We have not started to build a WL yet (dT==0) and there's no way it can occur now
                  ! since zQabs <= 0._wp
                  ! => no need to go further
                  l_exit = .TRUE.
               END IF

            END IF

            ! Okay test on updated absorbed flux:
            !LOLO: remove??? has a strong influence !!!
            IF ( (.NOT.(l_exit)) .AND. (Qnt_ac(ji,jj) + zQabs*rdt <= 0._wp) ) THEN
               l_exit       = .TRUE.
               l_destroy_wl = .TRUE.
            END IF


            IF ( .NOT. l_exit) THEN

               ! Two possibilities at this point:
               ! 1/ A warm layer already exists (dT>0) but it is cooling down because Qabs<0
               ! 2/ Regardless of WL formed (dT==0 or dT>0), we are in the process to initiate one or warm further it !

               ztac = Tau_ac(ji,jj) + MAX(.002_wp , pTau(ji,jj))*rdt      ! updated momentum integral
               !PRINT *, '#LBD: updated value for Tac=',  REAL(ztac,4)

               !! We update the value of absorbtion and zQabs:
               !! some part is useless if Qsw=0 !!!
               DO jl = 1, 5
                  zfr = 1. - ( 0.28*0.014*(1. - EXP(-zHwl/0.014)) + 0.27*0.357*(1. - EXP(-zHwl/0.357)) &
                     &        + 0.45*12.82*(1-EXP(-zHwl/12.82)) ) / zHwl
                  zQabs = zfr*pQsw(ji,jj) + pQnsol(ji,jj)
                  zqac  = Qnt_ac(ji,jj) + zQabs*rdt ! updated heat absorbed
                  IF ( zqac <= 0._wp ) EXIT
                  zHwl = MAX( MIN( Hwl_max , zcd1*ztac/SQRT(zqac)) , 0.1_wp ) ! Warm-layer depth
               END DO

               IF ( zqac <= 0._wp ) THEN
                  l_destroy_wl = .TRUE.
                  l_exit       = .TRUE.
               ELSE
                  zdTwl = zcd2*zqac**1.5/ztac * MAX(zqac/ABS(zqac),0._wp)  !! => IF(zqac>0._wp): zdTwl=zcd2*zqac**1.5/ztac ; ELSE: zdTwl=0. / ! normally: zqac > 0 !
                  !PRINT *, '#LBD: updated preliminary value for dT_wl=',  REAL(zdTwl,4)
                  ! Warm layer correction
                  flg = 0.5_wp + SIGN( 0.5_wp , gdept_1d(1)-zHwl )               ! => 1 when gdept_1d(1)>zHwl (zdTwl = zdTwl) | 0 when gdept_1d(1)<zHwl (zdTwl = zdTwl*gdept_1d(1)/zHwl)
                  zdTwl = zdTwl * ( flg + (1._wp-flg)*gdept_1d(1)/zHwl )
               END IF

            END IF !IF ( .NOT. l_exit)

            IF ( l_destroy_wl ) THEN
               zdTwl = 0._wp
               zfr   = 0.75_wp
               zHwl  = Hwl_max
               zqac  = 0._wp
               ztac  = 0._wp
            END IF

            IF ( iwait == 0 ) THEN
               !! Iteration loop within bulk algo is over, time to update what needs to be updated:
               dT_wl(ji,jj)  = zdTwl
               Hz_wl(ji,jj)  = zHwl
               Qnt_ac(ji,jj) = zqac ! Updating Qnt_ac, heat integral
               Tau_ac(ji,jj) = ztac
            END IF

         END DO
      END DO

   END SUBROUTINE WL_COARE




   FUNCTION delta_skin_layer( palpha, pQd, pQlat, pustar_a )
      !!---------------------------------------------------------------------
      !! Computes the thickness (m) of the viscous skin layer.
      !! Based on Fairall et al., 1996
      !!
      !! Fairall, C. W., Bradley, E. F., Godfrey, J. S., Wick, G. A.,
      !! Edson, J. B., and Young, G. S. ( 1996), Cool‐skin and warm‐layer
      !! effects on sea surface temperature, J. Geophys. Res., 101( C1), 1295-1308,
      !! doi:10.1029/95JC03190.
      !!
      !! L. Brodeau, october 2019
      !!---------------------------------------------------------------------
      REAL(wp)                :: delta_skin_layer
      REAL(wp), INTENT(in)    :: palpha   ! thermal expansion coefficient of sea-water (SST accurate enough!)
      REAL(wp), INTENT(in)    :: pQd    ! < 0 !!! part of the net heat flux actually absorbed in the WL [W/m^2] => term "Q + Rs*fs" in eq.6 of Fairall et al. 1996
      REAL(wp), INTENT(in)    :: pQlat    ! latent heat flux [W/m^2]
      REAL(wp), INTENT(in)    :: pustar_a ! friction velocity in the air (u*) [m/s]
      !!---------------------------------------------------------------------
      REAL(wp) :: zusw, zusw2, zlamb, zQd, ztf, ztmp
      !!---------------------------------------------------------------------
      
      zQd = pQd + 0.026*MIN(pQlat,0._wp)*rCp0_w/rLevap/palpha   ! LOLO: Double check sign + division by palpha !!! units are okay!

      ztf = 0.5_wp + SIGN(0.5_wp, zQd)  ! Qabs < 0 => cooling of the viscous layer => ztf = 0 (regular case)
      !                                 ! Qabs > 0 => warming of the viscous layer => ztf = 1 (ex: weak evaporation and strong positive sensible heat flux)
      !
      zusw  = MAX(pustar_a, 1.E-4_wp) * sq_radrw    ! u* in the water
      zusw2 = zusw*zusw
      !
      zlamb = 6._wp*( 1._wp + MAX(palpha*rcst_cs/(zusw2*zusw2)*zQd, 0._wp)**0.75 )**(-1./3.) ! see Eq.(14) in Fairall et al., 1996
      !  => zlamb is not used when Qd > 0, and since rcst_cs < 0, we just use this "MAX" to prevent FPE errors (something_negative)**0.75
      !
      ztmp = rnu0_w/zusw
      delta_skin_layer = (1._wp-ztf) *     zlamb*ztmp           &  ! regular case, Qd < 0, see Eq.(12) in Fairall et al., 1996
         &               +   ztf     * MIN(6._wp*ztmp , 0.007_wp)  ! when Qd > 0
   END FUNCTION delta_skin_layer

   !!======================================================================
END MODULE sbcblk_skin_coare