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) REAL(wp), PARAMETER :: zcon0 = -16._wp * grav * rho0_w * rCp0_w * rnu0_w*rnu0_w*rnu0_w / ( rk0_w*rk0_w ) ! "-" because ocean convention: Qabs > 0 => gain of heat for ocean! !! => see eq.(14) in Fairall et al. 1996 (eq.(6) of Zeng aand Beljaars is WRONG! (typo?) !! 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 scheme according to 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 = MIN( -0.1_wp , pQnsol(ji,jj) ) ! first guess, we do not miss a lot assuming 0 solar flux absorbed in the tiny layer of thicknes$ ! ! also, we ONLY consider when the viscous layer is loosing heat to the atmosphere, we only deal with cool-skin! => hence the "MIN( -0$ 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 = MIN( -0.1_wp , pQnsol(ji,jj) + zfr*pQsw(ji,jj) ) ! Total cooling at the interface zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pQlat(ji,jj), pustar(ji,jj) ) END DO dT_cs(ji,jj) = MIN( zQabs*zdelta/rk0_w , 0._wp ) ! temperature increment 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*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!!! PRINT *, '#LBD: Initial Qsw & Qnsol:', NINT(pQsw(ji,jj)), NINT(pQnsol(ji,jj)) PRINT *, '#LBD: =>Qabs:', zQabs,' zfr=', zfr 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 PRINT *, '#LBD: we have not started to to build a WL yet (dT==0)' PRINT *, '#LBD: and theres no way it can occur now since zQabs=', zQabs PRINT *, '#LBD: => leaving without changing anything...' 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 PRINT *, '#LBD: Oh boy! Next Qnt_ac looking weak! =>', Qnt_ac(ji,jj) + zQabs*rdt PRINT *, '#LBD: => time to destroy the warm-layer!' 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 ! PRINT *, '#LBD:======================================================' PRINT *, '#LBD: WL action makes sense now! => zQabs,dT_wl=', REAL(zQabs,4), REAL(zdTwl,4) PRINT *, '#LBD:======================================================' PRINT *, '#LBD: current values for Qac and Tac=', REAL(Qnt_ac(ji,jj),4), REAL(Tau_ac(ji,jj),4) 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 PRINT *, '#LBD: updated absorption and WL depth=', REAL(zfr,4), REAL(zHwl,4) PRINT *, '#LBD: updated value for Qabs=', REAL(zQabs,4), 'W/m2' PRINT *, '#LBD: updated value for Qac =', REAL(zqac,4), 'J' 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) 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, zQb !!--------------------------------------------------------------------- zQb = pQabs + 0.026*MIN(pQlat,0._wp)*rCp0_w/rLevap/palpha ! LOLO: Double check sign + division by palpha !!! units are okay! zusw = MAX(pustar_a, 1.E-4_wp) * sq_radrw ! u* in the water zusw2 = zusw*zusw zlamb = 6._wp*( 1._wp + (palpha*zcon0/(zusw2*zusw2)*zQb)**0.75 )**(-1./3.) ! see eq.(14) in Fairall et al., 1996 delta_skin_layer = zlamb*rnu0_w/zusw END FUNCTION delta_skin_layer !!====================================================================== END MODULE sbcblk_skin_coare