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), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: delta_vl !: thickness of the surface viscous layer (in the water) right below the air-sea interface [m] !! Warm-layer related parameters: REAL(wp), PARAMETER, PUBLIC :: H_wl_max = 20._wp !: maximum depth of warm layer (adjustable) ! REAL(wp), PARAMETER :: rich = 0.65_wp !: critical Richardson number ! REAL(wp), PARAMETER :: Qabs_thr = 50._wp !: threshold for heat flux absorbed in WL REAL(wp), PARAMETER :: zfr0 = 0.5_wp !: initial value of solar flux absorption ! REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: Tau_ac ! time integral / accumulated momentum Tauxdt => [N.s/m^2] (reset to zero every midnight) REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: Qnt_ac ! time integral / accumulated heat stored by the warm layer Qxdt => [J/m^2] (reset to zero every midnight) REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: H_wl ! depth of warm-layer [m] !!---------------------------------------------------------------------- CONTAINS SUBROUTINE CS_COARE( pQsw, pQnsol, pustar, pSST, pQlat, pdT ) !!--------------------------------------------------------------------- !! !! Cool-Skin scheme according to Fairall et al. 1996, revisited for 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 ! !! *pustar* friction velocity u* [m/s] !! *pSST* bulk SST (taken at depth gdept_1d(1)) [K] !! *pQlat* surface latent heat flux [K] !! !! ** INPUT/OUTPUT: !! *pdT* : as input => previous estimate of dT cool-skin !! as output => new estimate of dT cool-skin !! !!------------------------------------------------------------------ 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] !! REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pdT ! dT due to cool-skin effect !!--------------------------------------------------------------------- INTEGER :: ji, jj ! dummy loop indices REAL(wp) :: zQnet, zQnsol, zlamb, zdelta, zalpha_w, zfr, & & zz1, zz2, zus, & & ztf !!--------------------------------------------------------------------- DO jj = 1, jpj DO ji = 1, jpi zalpha_w = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!) zQnsol = MAX( 1._wp , - pQnsol(ji,jj) ) ! Non-solar heat loss to the atmosphere zdelta = delta_vl(ji,jj) ! using last value of delta !! Fraction of the shortwave flux absorbed by the cool-skin sublayer: zfr = MAX( 0.137_wp + 11._wp*zdelta - 6.6E-5_wp/zdelta*(1._wp - EXP(-zdelta/8.E-4_wp)) , 0.01_wp ) ! 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 ! zQnet = MAX( 1._wp , zQnsol - zfr*pQsw(ji,jj) ) ! Total cooling at the interface ztf = 0.5 + SIGN(0.5_wp, zQnet) ! Qt > 0 => cooling of the layer => ztf = 1 ! Qt < 0 => warming of the layer => ztf = 0 !! Term alpha*Qb (Qb is the virtual surface cooling inc. buoyancy effect of salinity due to evap): zz1 = zalpha_w*zQnet - 0.026*MIN(pQlat(ji,jj),0._wp)*rCp0_w/rLevap ! alpha*(Eq.8) == alpha*Qb "-" because Qlat < 0 !! LB: this terms only makes sense if > 0 i.e. in the cooling case !! so similar to what's done in ECMWF: zz1 = MAX(0._wp , zz1) ! 1. instead of 0.1 though ZQ = MAX(1.0,-pQlw(ji,jj) - pQsen(ji,jj) - pQlat(ji,jj)) zus = MAX(pustar(ji,jj), 1.E-4_wp) ! Laurent: too low wind (u*) might cause problem in stable cases: zz2 = zus*zus * roadrw zz2 = zz2*zz2 zlamb = 6._wp*( 1._wp + (rcst_cs*zz1/zz2)**0.75 )**(-1./3.) ! Lambda (Eq.14) (Saunders) ! Updating molecular sublayer thickness (delta): zz2 = rnu0_w/(SQRT(roadrw)*zus) zdelta = ztf * zlamb*zz2 & ! Eq.12 (when alpha*Qb>0 / cooling of layer) & + (1._wp - ztf) * MIN(0.007_wp , 6._wp*zz2 ) ! Eq.12 (when alpha*Qb<0 / warming of layer) !LB: changed 0.01 to 0.007 !! Once again with the new 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) zQnet = MAX( 1._wp , zQnsol - zfr*pQsw(ji,jj) ) ! Total cooling at the interface !! Update! pdT(ji,jj) = MIN( - zQnet*zdelta/rk0_w , 0._wp ) ! temperature increment delta_vl(ji,jj) = zdelta END DO END DO END SUBROUTINE CS_COARE SUBROUTINE WL_COARE( pQsw, pQnsol, pTau, pSST, iwait, pdT, Hwl, mask_wl ) !!--------------------------------------------------------------------- !! !! 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... !! !! ** OUTPUT: !! *pdT* dT due to warm-layer effect => difference between "almost surface (right below viscous layer) and depth of bulk SST !!--------------------------------------------------------------------- !! !! ** OPTIONAL OUTPUT: !! *Hwl* depth of warm layer [m] !! *mask_wl* mask for possible existence of a warm-layer (1) or not (0) !! !!------------------------------------------------------------------ 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 REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pdT ! dT due to warm-layer effect => difference between "almost surface (right below viscous layer) and depth of bulk SST !! REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: Hwl ! depth of warm layer [m] INTEGER(1), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: mask_wl ! mask for possible existence of a warm-layer (1) or not (0) ! ! INTEGER :: ji,jj ! REAL(wp) :: zdT_wl, zQabs, zfr, zdz REAL(wp) :: zqac, ztac REAL(wp) :: zalpha_w, zcd1, zcd2, flg !!--------------------------------------------------------------------- REAL(wp) :: ztime, znoon, zmidn INTEGER :: jl !! INITIALIZATION: pdT(:,:) = 0._wp ! dT initially set to 0._wp zdT_wl = 0._wp ! total warming (amplitude) in warm layer zQabs = 0._wp ! total heat absorped in warm layer zfr = zfr0 ! initial value of solar flux absorption ztac = 0._wp zqac = 0._wp IF ( PRESENT(mask_wl) ) mask_wl(:,:) = 0 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 zdz = MAX( MIN(H_wl(ji,jj),H_wl_max) , 0.1_wp) ! depth of warm layer !***** variables for warm layer *** zalpha_w = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!) zcd1 = SQRT(2._wp*rich*rCp0_w/(zalpha_w*grav*rho0_w)) !mess-o-constants 1 zcd2 = SQRT(2._wp*zalpha_w*grav/(rich*rho0_w))/(rCp0_w**1.5) !mess-o-constants 2 !******************************************************** !**** Compute apply warm layer correction ************* !******************************************************** znoon = MOD( 0.5_wp*(rdawn_dcy(ji,jj)+rdusk_dcy(ji,jj)), 1._wp ) ! 0 rnoon*24 = UTC time of local noon zmidn = MOD( znoon-0.5_wp , 1._wp ) IF (lwp .AND. (ji==2) .AND. (jj==2) ) WRITE(numout,*) 'LOLO: rdawn, rdusk, znoon, zmidn =', rdawn_dcy(ji,jj), rdusk_dcy(ji,jj), znoon, zmidn ! When midnight is past and dawn is not there yet, do what they call the "midnight reset": IF ( (ztime >= zmidn).AND.(ztime <= rdawn_dcy(ji,jj)) ) THEN IF (lwp .AND. (ji==2) .AND. (jj==2) ) WRITE(numout,*) ' LOLO [WL_COARE] MIDNIGHT RESET !!!!, ztime =', ztime zdz = H_wl_max Tau_ac(ji,jj) = 0._wp Qnt_ac(ji,jj) = 0._wp END IF IF ( ztime > rdawn_dcy(ji,jj) ) THEN ! Dawn, a new day starts at current location IF (lwp .AND. (ji==2) .AND. (jj==2) ) WRITE(numout,*) ' LOLO [WL_COARE] WE DO WL !!!!, ztime =', ztime !************************************ !**** set warm layer constants *** !************************************ zQabs = zfr*pQsw(ji,jj) + pQnsol(ji,jj) ! tot heat absorbed in warm layer !PRINT *, ' [WL_COARE] rdt, pQsw, pQnsol, zQabs =', rdt, REAL(pQsw(ji,jj),4), REAL(pQnsol(ji,jj),4), REAL(zQabs,4) IF ( zQabs >= Qabs_thr ) THEN ! Check for threshold !PRINT *, ' [WL_COARE] Tau_ac, Qnt_ac =', REAL(Tau_ac(ji,jj),4), REAL(Qnt_ac(ji,jj),4) !Tau_ac(ji,jj) = Tau_ac(ji,jj) + MAX(.002_wp , pTau(ji,jj))*rdt ! momentum integral ztac = Tau_ac(ji,jj) + MAX(.002_wp , pTau(ji,jj))*rdt ! updated momentum integral IF ( Qnt_ac(ji,jj) + zQabs*rdt > 0._wp ) THEN !check threshold for warm layer existence !****************************************** ! Compute the absorption profile !****************************************** DO jl = 1, 5 !loop 5 times for zfr zfr = 1. - ( 0.28*0.014*(1. - EXP(-zdz/0.014)) + 0.27*0.357*(1. - EXP(-zdz/0.357)) & & + 0.45*12.82*(1-EXP(-zdz/12.82)) ) / zdz zqac = Qnt_ac(ji,jj) + (zfr*pQsw(ji,jj) + pQnsol(ji,jj))*rdt ! updated heat absorbed IF ( zqac > 1._wp ) zdz = MAX( MIN( H_wl_max , zcd1*ztac/SQRT(zqac)) , 0.1_wp ) ! Warm-layer depth END DO ELSE !*********************** ! Warm layer wiped out !*********************** zfr = 0.75 zdz = H_wl_max zqac = Qnt_ac(ji,jj) + (zfr*pQsw(ji,jj) + pQnsol(ji,jj))*rdt ! updated heat absorbed END IF !IF ( Qnt_ac(ji,jj) + zQabs*rdt > 0._wp ) IF ( zqac > 1._wp ) zdT_wl = zcd2*zqac**1.5/ztac * MAX(zqac/ABS(zqac),0._wp) !! => IF(zqac>0._wp): zdT_wl=zcd2*zqac**1.5/ztac ; ELSE: zdT_wl=0. / ! normally: zqac > 0 ! ! Warm layer correction flg = 0.5_wp + SIGN( 0.5_wp , gdept_1d(1)-zdz ) ! => 1 when gdept_1d(1)>zdz (pdT(ji,jj) = zdT_wl) | 0 when gdept_1d(1)= Qabs_thr ) END IF ! IF ( isd_sol >= 21600 ) THEN ! (21600 == 6am) IF ( iwait == 0 ) THEN IF ( (zQabs >= Qabs_thr).AND.(ztime > rdawn_dcy(ji,jj)) ) THEN !PRINT *, ' [WL_COARE] WE UPDATE ACCUMULATED FLUXES !!!' Qnt_ac(ji,jj) = zqac ! Updating Qnt_ac, heat integral Tau_ac(ji,jj) = ztac ! IF ( PRESENT(mask_wl) ) mask_wl(ji,jj) = 1 END IF END IF H_wl(ji,jj) = zdz IF ( PRESENT(Hwl) ) Hwl(ji,jj) = H_wl(ji,jj) END DO END DO END SUBROUTINE WL_COARE !!====================================================================== END MODULE sbcblk_skin_coare