[11626] | 1 | MODULE sbcblk_skin_coare |
---|
[11615] | 2 | !!====================================================================== |
---|
[11626] | 3 | !! *** MODULE sbcblk_skin_coare *** |
---|
[11615] | 4 | !! Computes: |
---|
| 5 | !! * the surface skin temperature (aka SSST) based on the cool-skin/warm-layer |
---|
| 6 | !! scheme used at ECMWF |
---|
| 7 | !! Using formulation/param. of COARE 3.6 (Fairall et al., 2019) |
---|
| 8 | !! |
---|
| 9 | !! From routine turb_ecmwf maintained and developed in AeroBulk |
---|
| 10 | !! (https://github.com/brodeau/aerobulk) |
---|
| 11 | !! |
---|
| 12 | !! ** Author: L. Brodeau, September 2019 / AeroBulk (https://github.com/brodeau/aerobulk) |
---|
| 13 | !!---------------------------------------------------------------------- |
---|
| 14 | !! History : 4.0 ! 2016-02 (L.Brodeau) Original code |
---|
| 15 | !!---------------------------------------------------------------------- |
---|
| 16 | |
---|
| 17 | !!---------------------------------------------------------------------- |
---|
| 18 | !! cswl_ecmwf : computes the surface skin temperature (aka SSST) |
---|
| 19 | !! based on the cool-skin/warm-layer scheme used at ECMWF |
---|
| 20 | !!---------------------------------------------------------------------- |
---|
| 21 | USE oce ! ocean dynamics and tracers |
---|
| 22 | USE dom_oce ! ocean space and time domain |
---|
| 23 | USE phycst ! physical constants |
---|
| 24 | USE sbc_oce ! Surface boundary condition: ocean fields |
---|
| 25 | |
---|
| 26 | USE sbcblk_phy !LOLO: all thermodynamics functions, rho_air, q_sat, etc... |
---|
| 27 | |
---|
[11672] | 28 | USE sbcdcy !LOLO: to know hour of dawn and dusk: rdawn_dcy and rdusk_dcy (needed in WL_COARE) |
---|
| 29 | |
---|
[11615] | 30 | USE lib_mpp ! distribued memory computing library |
---|
| 31 | USE in_out_manager ! I/O manager |
---|
| 32 | USE lib_fortran ! to use key_nosignedzero |
---|
| 33 | |
---|
| 34 | IMPLICIT NONE |
---|
| 35 | PRIVATE |
---|
| 36 | |
---|
[11626] | 37 | PUBLIC :: CS_COARE, WL_COARE |
---|
[11615] | 38 | |
---|
| 39 | !! Cool-skin related parameters: |
---|
[11772] | 40 | REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: & |
---|
| 41 | & dT_cs !: dT due to cool-skin effect => temperature difference between air-sea interface (z=0) and right below viscous layer (z=delta) |
---|
[11615] | 42 | |
---|
| 43 | !! Warm-layer related parameters: |
---|
[11772] | 44 | REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: & |
---|
| 45 | & 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)) |
---|
| 46 | & Hz_wl, & !: depth of warm-layer [m] |
---|
| 47 | & Qnt_ac, & !: time integral / accumulated heat stored by the warm layer Qxdt => [J/m^2] (reset to zero every midnight) |
---|
| 48 | & Tau_ac !: time integral / accumulated momentum Tauxdt => [N.s/m^2] (reset to zero every midnight) |
---|
| 49 | |
---|
| 50 | REAL(wp), PARAMETER, PUBLIC :: Hwl_max = 20._wp !: maximum depth of warm layer (adjustable) |
---|
[11615] | 51 | ! |
---|
| 52 | REAL(wp), PARAMETER :: rich = 0.65_wp !: critical Richardson number |
---|
| 53 | ! |
---|
| 54 | REAL(wp), PARAMETER :: zfr0 = 0.5_wp !: initial value of solar flux absorption |
---|
| 55 | ! |
---|
| 56 | !!---------------------------------------------------------------------- |
---|
| 57 | CONTAINS |
---|
| 58 | |
---|
| 59 | |
---|
[11772] | 60 | SUBROUTINE CS_COARE( pQsw, pQnsol, pustar, pSST, pQlat ) |
---|
[11615] | 61 | !!--------------------------------------------------------------------- |
---|
| 62 | !! |
---|
[11804] | 63 | !! Cool-skin parameterization, based on Fairall et al., 1996, revisited for COARE 3.6 (Fairall et al., 2019) |
---|
[11615] | 64 | !! |
---|
[11772] | 65 | !! Fairall, C. W., Bradley, E. F., Godfrey, J. S., Wick, G. A., |
---|
| 66 | !! Edson, J. B., and Young, G. S. ( 1996), Cool‐skin and warm‐layer |
---|
| 67 | !! effects on sea surface temperature, J. Geophys. Res., 101( C1), 1295-1308, |
---|
| 68 | !! doi:10.1029/95JC03190. |
---|
| 69 | !! |
---|
| 70 | !!------------------------------------------------------------------ |
---|
| 71 | !! |
---|
[11615] | 72 | !! ** INPUT: |
---|
| 73 | !! *pQsw* surface net solar radiation into the ocean [W/m^2] => >= 0 ! |
---|
| 74 | !! *pQnsol* surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 ! |
---|
| 75 | !! *pustar* friction velocity u* [m/s] |
---|
| 76 | !! *pSST* bulk SST (taken at depth gdept_1d(1)) [K] |
---|
| 77 | !! *pQlat* surface latent heat flux [K] |
---|
| 78 | !!------------------------------------------------------------------ |
---|
[11804] | 79 | REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQsw ! net solar a.k.a shortwave radiation into the ocean (after albedo) [W/m^2] |
---|
| 80 | REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQnsol ! non-solar heat flux to the ocean [W/m^2] |
---|
| 81 | REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pustar ! friction velocity, temperature and humidity (u*,t*,q*) |
---|
| 82 | REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pSST ! bulk SST [K] |
---|
| 83 | REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQlat ! latent heat flux [W/m^2] |
---|
[11615] | 84 | !!--------------------------------------------------------------------- |
---|
[11772] | 85 | INTEGER :: ji, jj, jc |
---|
| 86 | REAL(wp) :: zQabs, zdelta, zfr |
---|
[11615] | 87 | !!--------------------------------------------------------------------- |
---|
| 88 | DO jj = 1, jpj |
---|
| 89 | DO ji = 1, jpi |
---|
| 90 | |
---|
[11804] | 91 | zQabs = pQnsol(ji,jj) ! first guess of heat flux absorbed within the viscous sublayer of thicknes delta, |
---|
| 92 | ! ! => we DO not miss a lot assuming 0 solar flux absorbed in the tiny layer of thicknes zdelta... |
---|
[11615] | 93 | |
---|
[11772] | 94 | zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pQlat(ji,jj), pustar(ji,jj) ) |
---|
[11615] | 95 | |
---|
[11772] | 96 | DO jc = 1, 4 ! because implicit in terms of zdelta... |
---|
[11804] | 97 | 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 ! |
---|
| 98 | zQabs = pQnsol(ji,jj) + zfr*pQsw(ji,jj) |
---|
[11772] | 99 | zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pQlat(ji,jj), pustar(ji,jj) ) |
---|
| 100 | END DO |
---|
[11615] | 101 | |
---|
[11804] | 102 | dT_cs(ji,jj) = zQabs*zdelta/rk0_w ! temperature increment, yes dT_cs can actually > 0, if Qabs > 0 (rare but possible!) |
---|
[11615] | 103 | |
---|
| 104 | END DO |
---|
| 105 | END DO |
---|
| 106 | |
---|
[11626] | 107 | END SUBROUTINE CS_COARE |
---|
[11615] | 108 | |
---|
| 109 | |
---|
[11772] | 110 | SUBROUTINE WL_COARE( pQsw, pQnsol, pTau, pSST, iwait ) |
---|
[11615] | 111 | !!--------------------------------------------------------------------- |
---|
| 112 | !! |
---|
| 113 | !! Warm-Layer scheme according to COARE 3.6 (Fairall et al, 2019) |
---|
| 114 | !! ------------------------------------------------------------------ |
---|
| 115 | !! |
---|
| 116 | !! ** INPUT: |
---|
| 117 | !! *pQsw* surface net solar radiation into the ocean [W/m^2] => >= 0 ! |
---|
| 118 | !! *pQnsol* surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 ! |
---|
| 119 | !! *pTau* surface wind stress [N/m^2] |
---|
| 120 | !! *pSST* bulk SST (taken at depth gdept_1d(1)) [K] |
---|
| 121 | !! *iwait* if /= 0 then wait before updating accumulated fluxes, we are within a converging itteration loop... |
---|
| 122 | !!--------------------------------------------------------------------- |
---|
| 123 | REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQsw ! surface net solar radiation into the ocean [W/m^2] => >= 0 ! |
---|
| 124 | REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQnsol ! surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 ! |
---|
| 125 | REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pTau ! wind stress [N/m^2] |
---|
| 126 | REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pSST ! bulk SST at depth gdept_1d(1) [K] |
---|
| 127 | INTEGER , INTENT(in) :: iwait ! if /= 0 then wait before updating accumulated fluxes |
---|
| 128 | !! |
---|
| 129 | INTEGER :: ji,jj |
---|
| 130 | ! |
---|
[11772] | 131 | REAL(wp) :: zdTwl, zHwl, zQabs, zfr |
---|
[11615] | 132 | REAL(wp) :: zqac, ztac |
---|
[11772] | 133 | REAL(wp) :: zalpha, zcd1, zcd2, flg |
---|
[11615] | 134 | !!--------------------------------------------------------------------- |
---|
| 135 | |
---|
[11672] | 136 | REAL(wp) :: ztime, znoon, zmidn |
---|
| 137 | INTEGER :: jl |
---|
[11615] | 138 | |
---|
[11772] | 139 | LOGICAL :: l_exit, l_destroy_wl |
---|
| 140 | |
---|
[11615] | 141 | !! INITIALIZATION: |
---|
[11772] | 142 | zQabs = 0._wp ! total heat flux absorped in warm layer |
---|
| 143 | zfr = zfr0 ! initial value of solar flux absorption !LOLO: save it and use previous value !!! |
---|
[11615] | 144 | |
---|
[11672] | 145 | 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! |
---|
| 146 | |
---|
| 147 | 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 ... |
---|
| 148 | |
---|
| 149 | IF (lwp) THEN |
---|
| 150 | WRITE(numout,*) '' |
---|
| 151 | WRITE(numout,*) 'LOLO: sbcblk_skin_coare => nsec_day, ztime =', nsec_day, ztime |
---|
| 152 | END IF |
---|
| 153 | |
---|
[11615] | 154 | DO jj = 1, jpj |
---|
| 155 | DO ji = 1, jpi |
---|
| 156 | |
---|
[11772] | 157 | l_exit = .FALSE. |
---|
| 158 | l_destroy_wl = .FALSE. |
---|
[11615] | 159 | |
---|
[11772] | 160 | zdTwl = dT_wl(ji,jj) ! value of previous time step as first guess |
---|
| 161 | zHwl = MAX( MIN(Hz_wl(ji,jj),Hwl_max),0.1_wp) ! " " " |
---|
| 162 | |
---|
| 163 | zqac = Qnt_ac(ji,jj) ! previous time step Qnt_ac |
---|
| 164 | ztac = Tau_ac(ji,jj) |
---|
| 165 | |
---|
[11615] | 166 | !***** variables for warm layer *** |
---|
[11772] | 167 | zalpha = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!) |
---|
[11615] | 168 | |
---|
[11772] | 169 | zcd1 = SQRT(2._wp*rich*rCp0_w/(zalpha*grav*rho0_w)) !mess-o-constants 1 |
---|
| 170 | zcd2 = SQRT(2._wp*zalpha*grav/(rich*rho0_w))/(rCp0_w**1.5) !mess-o-constants 2 |
---|
[11615] | 171 | |
---|
| 172 | |
---|
[11672] | 173 | 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 |
---|
| 174 | zmidn = MOD( znoon-0.5_wp , 1._wp ) |
---|
[11772] | 175 | zmidn = MOD( zmidn + 0.125_wp , 1._wp ) ! 3 hours past the local midnight |
---|
[11672] | 176 | |
---|
[11772] | 177 | IF ( (ztime >= zmidn) .AND. (ztime < rdawn_dcy(ji,jj)) ) THEN |
---|
| 178 | ! Dawn reset to 0! |
---|
| 179 | l_exit = .TRUE. |
---|
| 180 | l_destroy_wl = .TRUE. |
---|
[11615] | 181 | END IF |
---|
| 182 | |
---|
[11772] | 183 | IF ( .NOT. l_exit ) THEN |
---|
| 184 | !! Initial test on initial guess of absorbed heat flux in warm-layer: |
---|
| 185 | zfr = 1._wp - ( 0.28*0.014*(1. - EXP(-zHwl/0.014)) + 0.27*0.357*(1. - EXP(-zHwl/0.357)) & |
---|
| 186 | & + 0.45*12.82*(1-EXP(-zHwl/12.82)) ) / zHwl |
---|
| 187 | 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!!! |
---|
[11615] | 188 | |
---|
[11772] | 189 | IF ( (ABS(zdTwl) < 1.E-6_wp) .AND. (zQabs <= 0._wp) ) THEN |
---|
| 190 | ! We have not started to build a WL yet (dT==0) and there's no way it can occur now |
---|
| 191 | ! since zQabs <= 0._wp |
---|
| 192 | ! => no need to go further |
---|
| 193 | l_exit = .TRUE. |
---|
| 194 | END IF |
---|
[11615] | 195 | |
---|
[11772] | 196 | END IF |
---|
[11615] | 197 | |
---|
[11772] | 198 | ! Okay test on updated absorbed flux: |
---|
| 199 | !LOLO: remove??? has a strong influence !!! |
---|
| 200 | IF ( (.NOT.(l_exit)) .AND. (Qnt_ac(ji,jj) + zQabs*rdt <= 0._wp) ) THEN |
---|
| 201 | l_exit = .TRUE. |
---|
| 202 | l_destroy_wl = .TRUE. |
---|
| 203 | END IF |
---|
[11615] | 204 | |
---|
| 205 | |
---|
[11772] | 206 | IF ( .NOT. l_exit) THEN |
---|
[11615] | 207 | |
---|
[11772] | 208 | ! Two possibilities at this point: |
---|
| 209 | ! 1/ A warm layer already exists (dT>0) but it is cooling down because Qabs<0 |
---|
| 210 | ! 2/ Regardless of WL formed (dT==0 or dT>0), we are in the process to initiate one or warm further it ! |
---|
[11615] | 211 | |
---|
[11772] | 212 | ztac = Tau_ac(ji,jj) + MAX(.002_wp , pTau(ji,jj))*rdt ! updated momentum integral |
---|
[11775] | 213 | !PRINT *, '#LBD: updated value for Tac=', REAL(ztac,4) |
---|
[11615] | 214 | |
---|
[11772] | 215 | !! We update the value of absorbtion and zQabs: |
---|
| 216 | !! some part is useless if Qsw=0 !!! |
---|
| 217 | DO jl = 1, 5 |
---|
| 218 | zfr = 1. - ( 0.28*0.014*(1. - EXP(-zHwl/0.014)) + 0.27*0.357*(1. - EXP(-zHwl/0.357)) & |
---|
| 219 | & + 0.45*12.82*(1-EXP(-zHwl/12.82)) ) / zHwl |
---|
| 220 | zQabs = zfr*pQsw(ji,jj) + pQnsol(ji,jj) |
---|
| 221 | zqac = Qnt_ac(ji,jj) + zQabs*rdt ! updated heat absorbed |
---|
| 222 | IF ( zqac <= 0._wp ) EXIT |
---|
| 223 | zHwl = MAX( MIN( Hwl_max , zcd1*ztac/SQRT(zqac)) , 0.1_wp ) ! Warm-layer depth |
---|
| 224 | END DO |
---|
[11615] | 225 | |
---|
[11772] | 226 | IF ( zqac <= 0._wp ) THEN |
---|
| 227 | l_destroy_wl = .TRUE. |
---|
| 228 | l_exit = .TRUE. |
---|
| 229 | ELSE |
---|
| 230 | 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 ! |
---|
[11775] | 231 | !PRINT *, '#LBD: updated preliminary value for dT_wl=', REAL(zdTwl,4) |
---|
[11615] | 232 | ! Warm layer correction |
---|
[11772] | 233 | 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) |
---|
| 234 | zdTwl = zdTwl * ( flg + (1._wp-flg)*gdept_1d(1)/zHwl ) |
---|
| 235 | END IF |
---|
[11615] | 236 | |
---|
[11772] | 237 | END IF !IF ( .NOT. l_exit) |
---|
[11615] | 238 | |
---|
[11772] | 239 | IF ( l_destroy_wl ) THEN |
---|
| 240 | zdTwl = 0._wp |
---|
| 241 | zfr = 0.75_wp |
---|
| 242 | zHwl = Hwl_max |
---|
| 243 | zqac = 0._wp |
---|
| 244 | ztac = 0._wp |
---|
[11615] | 245 | END IF |
---|
| 246 | |
---|
[11772] | 247 | IF ( iwait == 0 ) THEN |
---|
| 248 | !! Iteration loop within bulk algo is over, time to update what needs to be updated: |
---|
| 249 | dT_wl(ji,jj) = zdTwl |
---|
| 250 | Hz_wl(ji,jj) = zHwl |
---|
| 251 | Qnt_ac(ji,jj) = zqac ! Updating Qnt_ac, heat integral |
---|
| 252 | Tau_ac(ji,jj) = ztac |
---|
| 253 | END IF |
---|
[11615] | 254 | |
---|
| 255 | END DO |
---|
| 256 | END DO |
---|
[11772] | 257 | |
---|
[11626] | 258 | END SUBROUTINE WL_COARE |
---|
[11615] | 259 | |
---|
| 260 | |
---|
[11772] | 261 | |
---|
| 262 | |
---|
[11804] | 263 | FUNCTION delta_skin_layer( palpha, pQd, pQlat, pustar_a ) |
---|
[11772] | 264 | !!--------------------------------------------------------------------- |
---|
| 265 | !! Computes the thickness (m) of the viscous skin layer. |
---|
| 266 | !! Based on Fairall et al., 1996 |
---|
| 267 | !! |
---|
| 268 | !! Fairall, C. W., Bradley, E. F., Godfrey, J. S., Wick, G. A., |
---|
| 269 | !! Edson, J. B., and Young, G. S. ( 1996), Cool‐skin and warm‐layer |
---|
| 270 | !! effects on sea surface temperature, J. Geophys. Res., 101( C1), 1295-1308, |
---|
| 271 | !! doi:10.1029/95JC03190. |
---|
| 272 | !! |
---|
| 273 | !! L. Brodeau, october 2019 |
---|
| 274 | !!--------------------------------------------------------------------- |
---|
| 275 | REAL(wp) :: delta_skin_layer |
---|
| 276 | REAL(wp), INTENT(in) :: palpha ! thermal expansion coefficient of sea-water (SST accurate enough!) |
---|
[11804] | 277 | 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 |
---|
[11772] | 278 | REAL(wp), INTENT(in) :: pQlat ! latent heat flux [W/m^2] |
---|
| 279 | REAL(wp), INTENT(in) :: pustar_a ! friction velocity in the air (u*) [m/s] |
---|
| 280 | !!--------------------------------------------------------------------- |
---|
[11804] | 281 | REAL(wp) :: zusw, zusw2, zlamb, zQd, ztf, ztmp |
---|
[11772] | 282 | !!--------------------------------------------------------------------- |
---|
[11804] | 283 | |
---|
| 284 | zQd = pQd + 0.026*MIN(pQlat,0._wp)*rCp0_w/rLevap/palpha ! LOLO: Double check sign + division by palpha !!! units are okay! |
---|
[11772] | 285 | |
---|
[11804] | 286 | ztf = 0.5_wp + SIGN(0.5_wp, zQd) ! Qabs < 0 => cooling of the viscous layer => ztf = 0 (regular case) |
---|
| 287 | ! ! Qabs > 0 => warming of the viscous layer => ztf = 1 (ex: weak evaporation and strong positive sensible heat flux) |
---|
| 288 | ! |
---|
[11772] | 289 | zusw = MAX(pustar_a, 1.E-4_wp) * sq_radrw ! u* in the water |
---|
| 290 | zusw2 = zusw*zusw |
---|
[11804] | 291 | ! |
---|
[11826] | 292 | 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 |
---|
| 293 | ! => 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 |
---|
[11804] | 294 | ! |
---|
| 295 | ztmp = rnu0_w/zusw |
---|
| 296 | delta_skin_layer = (1._wp-ztf) * zlamb*ztmp & ! regular case, Qd < 0, see Eq.(12) in Fairall et al., 1996 |
---|
| 297 | & + ztf * MIN(6._wp*ztmp , 0.007_wp) ! when Qd > 0 |
---|
[11772] | 298 | END FUNCTION delta_skin_layer |
---|
| 299 | |
---|
[11615] | 300 | !!====================================================================== |
---|
[11626] | 301 | END MODULE sbcblk_skin_coare |
---|