[11615] | 1 | MODULE sbcblk_skin_coare3p6 |
---|
| 2 | !!====================================================================== |
---|
| 3 | !! *** MODULE sbcblk_skin_coare3p6 *** |
---|
| 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 | |
---|
| 28 | USE lib_mpp ! distribued memory computing library |
---|
| 29 | USE in_out_manager ! I/O manager |
---|
| 30 | USE lib_fortran ! to use key_nosignedzero |
---|
| 31 | |
---|
| 32 | IMPLICIT NONE |
---|
| 33 | PRIVATE |
---|
| 34 | |
---|
| 35 | PUBLIC :: CS_COARE3P6, WL_COARE3P6 |
---|
| 36 | |
---|
| 37 | !! Cool-skin related parameters: |
---|
| 38 | REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: delta_vl !: thickness of the surface viscous layer (in the water) right below the air-sea interface [m] |
---|
| 39 | |
---|
| 40 | !! Warm-layer related parameters: |
---|
| 41 | REAL(wp), PARAMETER, PUBLIC :: H_wl_max = 20._wp !: maximum depth of warm layer (adjustable) |
---|
| 42 | ! |
---|
| 43 | REAL(wp), PARAMETER :: rich = 0.65_wp !: critical Richardson number |
---|
| 44 | ! |
---|
| 45 | REAL(wp), PARAMETER :: Qabs_thr = 50._wp !: threshold for heat flux absorbed in WL |
---|
| 46 | REAL(wp), PARAMETER :: zfr0 = 0.5_wp !: initial value of solar flux absorption |
---|
| 47 | ! |
---|
| 48 | REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: Tau_ac ! time integral / accumulated momentum Tauxdt => [N.s/m^2] (reset to zero every midnight) |
---|
| 49 | 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) |
---|
| 50 | REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: H_wl ! depth of warm-layer [m] |
---|
| 51 | !!---------------------------------------------------------------------- |
---|
| 52 | CONTAINS |
---|
| 53 | |
---|
| 54 | |
---|
| 55 | SUBROUTINE CS_COARE3P6( pQsw, pQnsol, pustar, pSST, pQlat, pdT ) |
---|
| 56 | !!--------------------------------------------------------------------- |
---|
| 57 | !! |
---|
| 58 | !! Cool-Skin scheme according to Fairall et al. 1996, revisited for COARE 3.6 (Fairall et al., 2019) |
---|
| 59 | !! ------------------------------------------------------------------ |
---|
| 60 | !! |
---|
| 61 | !! ** INPUT: |
---|
| 62 | !! *pQsw* surface net solar radiation into the ocean [W/m^2] => >= 0 ! |
---|
| 63 | !! *pQnsol* surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 ! |
---|
| 64 | !! *pustar* friction velocity u* [m/s] |
---|
| 65 | !! *pSST* bulk SST (taken at depth gdept_1d(1)) [K] |
---|
| 66 | !! *pQlat* surface latent heat flux [K] |
---|
| 67 | !! |
---|
| 68 | !! ** INPUT/OUTPUT: |
---|
| 69 | !! *pdT* : as input => previous estimate of dT cool-skin |
---|
| 70 | !! as output => new estimate of dT cool-skin |
---|
| 71 | !! |
---|
| 72 | !!------------------------------------------------------------------ |
---|
| 73 | REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQsw ! net solar a.k.a shortwave radiation into the ocean (after albedo) [W/m^2] |
---|
| 74 | REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQnsol ! non-solar heat flux to the ocean [W/m^2] |
---|
| 75 | REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pustar ! friction velocity, temperature and humidity (u*,t*,q*) |
---|
| 76 | REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pSST ! bulk SST [K] |
---|
| 77 | REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQlat ! latent heat flux [W/m^2] |
---|
| 78 | !! |
---|
| 79 | REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pdT ! dT due to cool-skin effect |
---|
| 80 | !!--------------------------------------------------------------------- |
---|
| 81 | INTEGER :: ji, jj ! dummy loop indices |
---|
| 82 | REAL(wp) :: zQnet, zQnsol, zlamb, zdelta, zalpha_w, zfr, & |
---|
| 83 | & zz1, zz2, zus, & |
---|
| 84 | & ztf |
---|
| 85 | !!--------------------------------------------------------------------- |
---|
| 86 | |
---|
| 87 | DO jj = 1, jpj |
---|
| 88 | DO ji = 1, jpi |
---|
| 89 | |
---|
| 90 | zalpha_w = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!) |
---|
| 91 | |
---|
| 92 | zQnsol = MAX( 1._wp , - pQnsol(ji,jj) ) ! Non-solar heat loss to the atmosphere |
---|
| 93 | |
---|
| 94 | zdelta = delta_vl(ji,jj) ! using last value of delta |
---|
| 95 | |
---|
| 96 | !! Fraction of the shortwave flux absorbed by the cool-skin sublayer: |
---|
| 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 ) ! 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 | |
---|
| 99 | zQnet = MAX( 1._wp , zQnsol - zfr*pQsw(ji,jj) ) ! Total cooling at the interface |
---|
| 100 | |
---|
| 101 | ztf = 0.5 + SIGN(0.5_wp, zQnet) ! Qt > 0 => cooling of the layer => ztf = 1 |
---|
| 102 | ! Qt < 0 => warming of the layer => ztf = 0 |
---|
| 103 | |
---|
| 104 | !! Term alpha*Qb (Qb is the virtual surface cooling inc. buoyancy effect of salinity due to evap): |
---|
| 105 | zz1 = zalpha_w*zQnet - 0.026*MIN(pQlat(ji,jj),0._wp)*rCp0_w/rLevap ! alpha*(Eq.8) == alpha*Qb "-" because Qlat < 0 |
---|
| 106 | !! LB: this terms only makes sense if > 0 i.e. in the cooling case |
---|
| 107 | !! so similar to what's done in ECMWF: |
---|
| 108 | zz1 = MAX(0._wp , zz1) ! 1. instead of 0.1 though ZQ = MAX(1.0,-pQlw(ji,jj) - pQsen(ji,jj) - pQlat(ji,jj)) |
---|
| 109 | |
---|
| 110 | zus = MAX(pustar(ji,jj), 1.E-4_wp) ! Laurent: too low wind (u*) might cause problem in stable cases: |
---|
| 111 | zz2 = zus*zus * roadrw |
---|
| 112 | zz2 = zz2*zz2 |
---|
| 113 | zlamb = 6._wp*( 1._wp + (rcst_cs*zz1/zz2)**0.75 )**(-1./3.) ! Lambda (Eq.14) (Saunders) |
---|
| 114 | |
---|
| 115 | ! Updating molecular sublayer thickness (delta): |
---|
| 116 | zz2 = rnu0_w/(SQRT(roadrw)*zus) |
---|
| 117 | zdelta = ztf * zlamb*zz2 & ! Eq.12 (when alpha*Qb>0 / cooling of layer) |
---|
| 118 | & + (1._wp - ztf) * MIN(0.007_wp , 6._wp*zz2 ) ! Eq.12 (when alpha*Qb<0 / warming of layer) |
---|
| 119 | !LB: changed 0.01 to 0.007 |
---|
| 120 | |
---|
| 121 | !! Once again with the new zdelta: |
---|
| 122 | 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) |
---|
| 123 | zQnet = MAX( 1._wp , zQnsol - zfr*pQsw(ji,jj) ) ! Total cooling at the interface |
---|
| 124 | |
---|
| 125 | !! Update! |
---|
| 126 | pdT(ji,jj) = MIN( - zQnet*zdelta/rk0_w , 0._wp ) ! temperature increment |
---|
| 127 | delta_vl(ji,jj) = zdelta |
---|
| 128 | |
---|
| 129 | END DO |
---|
| 130 | END DO |
---|
| 131 | |
---|
| 132 | END SUBROUTINE CS_COARE3P6 |
---|
| 133 | |
---|
| 134 | |
---|
| 135 | |
---|
| 136 | |
---|
| 137 | SUBROUTINE WL_COARE3P6( kt, pQsw, pQnsol, pTau, pSST, plon, isd, iwait, pdT, & |
---|
| 138 | & Hwl, mask_wl ) |
---|
| 139 | !!--------------------------------------------------------------------- |
---|
| 140 | !! |
---|
| 141 | !! Warm-Layer scheme according to COARE 3.6 (Fairall et al, 2019) |
---|
| 142 | !! ------------------------------------------------------------------ |
---|
| 143 | !! |
---|
| 144 | !! ** INPUT: |
---|
| 145 | !! *pQsw* surface net solar radiation into the ocean [W/m^2] => >= 0 ! |
---|
| 146 | !! *pQnsol* surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 ! |
---|
| 147 | !! *pTau* surface wind stress [N/m^2] |
---|
| 148 | !! *pSST* bulk SST (taken at depth gdept_1d(1)) [K] |
---|
| 149 | !! *plon* longitude [deg.E] |
---|
| 150 | !! *isd* current UTC time, counted in second since 00h of the current day |
---|
| 151 | !! *iwait* if /= 0 then wait before updating accumulated fluxes, we are within a converging itteration loop... |
---|
| 152 | !! |
---|
| 153 | !! ** OUTPUT: |
---|
| 154 | !! *pdT* dT due to warming at depth of pSST such that SST_actual = pSST + pdT |
---|
| 155 | !!--------------------------------------------------------------------- |
---|
| 156 | !! |
---|
| 157 | !! ** OPTIONAL OUTPUT: |
---|
| 158 | !! *Hwl* depth of warm layer [m] |
---|
| 159 | !! *mask_wl* mask for possible existence of a warm-layer (1) or not (0) |
---|
| 160 | !! |
---|
| 161 | !!------------------------------------------------------------------ |
---|
| 162 | INTEGER , INTENT(in) :: kt |
---|
| 163 | REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQsw ! surface net solar radiation into the ocean [W/m^2] => >= 0 ! |
---|
| 164 | REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQnsol ! surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 ! |
---|
| 165 | REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pTau ! wind stress [N/m^2] |
---|
| 166 | REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pSST ! bulk SST at depth gdept_1d(1) [K] |
---|
| 167 | REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: plon ! longitude [deg.E] |
---|
| 168 | INTEGER , INTENT(in) :: isd ! current UTC time, counted in second since 00h of the current day |
---|
| 169 | INTEGER , INTENT(in) :: iwait ! if /= 0 then wait before updating accumulated fluxes |
---|
| 170 | REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pdT ! dT due to warming at depth of pSST such that pSST_actual = pSST + pdT |
---|
| 171 | !! |
---|
| 172 | REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: Hwl ! depth of warm layer [m] |
---|
| 173 | INTEGER(1), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: mask_wl ! mask for possible existence of a warm-layer (1) or not (0) |
---|
| 174 | ! |
---|
| 175 | ! |
---|
| 176 | INTEGER :: ji,jj |
---|
| 177 | ! |
---|
| 178 | REAL(wp) :: zdT_wl, zQabs, zfr, zdz |
---|
| 179 | REAL(wp) :: zqac, ztac |
---|
| 180 | REAL(wp) :: zalpha_w, zcd1, zcd2, flg |
---|
| 181 | CHARACTER(len=128) :: cf_tmp |
---|
| 182 | !!--------------------------------------------------------------------- |
---|
| 183 | |
---|
| 184 | REAL(wp) :: rlag_gw_h, & ! local solar time lag in hours / Greenwich meridian (lon==0) => ex: ~ -10.47 hours for Hawai |
---|
| 185 | & rhr_sol ! local solar time in hours since midnight |
---|
| 186 | |
---|
| 187 | INTEGER :: ilag_gw_s, & ! local solar time LAG in seconds / Greenwich meridian (lon==0) => ex: ~ INT( -10.47*3600. ) seconds for Hawai |
---|
| 188 | & isd_sol, & ! local solar time in number of seconds since local solar midnight |
---|
| 189 | & jl |
---|
| 190 | |
---|
| 191 | !! INITIALIZATION: |
---|
| 192 | pdT(:,:) = 0._wp ! dT initially set to 0._wp |
---|
| 193 | zdT_wl = 0._wp ! total warming (amplitude) in warm layer |
---|
| 194 | zQabs = 0._wp ! total heat absorped in warm layer |
---|
| 195 | zfr = zfr0 ! initial value of solar flux absorption |
---|
| 196 | ztac = 0._wp |
---|
| 197 | zqac = 0._wp |
---|
| 198 | IF ( PRESENT(mask_wl) ) mask_wl(:,:) = 0 |
---|
| 199 | |
---|
| 200 | DO jj = 1, jpj |
---|
| 201 | DO ji = 1, jpi |
---|
| 202 | |
---|
| 203 | zdz = MAX( MIN(H_wl(ji,jj),H_wl_max) , 0.1_wp) ! depth of warm layer |
---|
| 204 | |
---|
| 205 | !! Need to know the local solar time from longitude and isd: |
---|
| 206 | rlag_gw_h = -1._wp * MODULO( ( 360._wp - MODULO(plon(ji,jj),360._wp) ) / 15._wp , 24._wp ) |
---|
| 207 | rlag_gw_h = -1._wp * SIGN( MIN(ABS(rlag_gw_h) , ABS(MODULO(rlag_gw_h,24._wp))), rlag_gw_h + 12._wp ) |
---|
| 208 | ilag_gw_s = INT( rlag_gw_h*3600._wp ) |
---|
| 209 | isd_sol = MODULO( isd + ilag_gw_s , 24*3600 ) |
---|
| 210 | rhr_sol = REAL( isd_sol , wp) / 3600._wp |
---|
| 211 | |
---|
| 212 | !***** variables for warm layer *** |
---|
| 213 | zalpha_w = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!) |
---|
| 214 | |
---|
| 215 | zcd1 = SQRT(2._wp*rich*rCp0_w/(zalpha_w*grav*rho0_w)) !mess-o-constants 1 |
---|
| 216 | zcd2 = SQRT(2._wp*zalpha_w*grav/(rich*rho0_w))/(rCp0_w**1.5) !mess-o-constants 2 |
---|
| 217 | |
---|
| 218 | !******************************************************** |
---|
| 219 | !**** Compute apply warm layer correction ************* |
---|
| 220 | !******************************************************** |
---|
| 221 | |
---|
| 222 | !IF (isd_sol <= rdt ) THEN !re-zero at midnight ! LOLO improve: risky if real midnight (00:00:00) is not a time in vtime... |
---|
| 223 | IF ( (rhr_sol > 23.5_wp).OR.(rhr_sol < 4._wp) ) THEN |
---|
| 224 | !PRINT *, ' [WL_COARE3P6] MIDNIGHT RESET !!!!, isd_sol =>', isd_sol |
---|
| 225 | zdz = H_wl_max |
---|
| 226 | Tau_ac(ji,jj) = 0._wp |
---|
| 227 | Qnt_ac(ji,jj) = 0._wp |
---|
| 228 | END IF |
---|
| 229 | |
---|
| 230 | IF ( rhr_sol > 5._wp ) THEN ! ( 5am) |
---|
| 231 | |
---|
| 232 | !PRINT *, ' [WL_COARE3P6] WE DO WL !' |
---|
| 233 | !PRINT *, ' [WL_COARE3P6] isd_sol, pTau, pSST, pdT =', isd_sol, REAL(pTau(ji,jj),4), REAL(pSST(ji,jj),4), REAL(pdT(ji,jj),4) |
---|
| 234 | |
---|
| 235 | !************************************ |
---|
| 236 | !**** set warm layer constants *** |
---|
| 237 | !************************************ |
---|
| 238 | |
---|
| 239 | zQabs = zfr*pQsw(ji,jj) + pQnsol(ji,jj) ! tot heat absorbed in warm layer |
---|
| 240 | |
---|
| 241 | !PRINT *, ' [WL_COARE3P6] rdt, pQsw, pQnsol, zQabs =', rdt, REAL(pQsw(ji,jj),4), REAL(pQnsol(ji,jj),4), REAL(zQabs,4) |
---|
| 242 | |
---|
| 243 | IF ( zQabs >= Qabs_thr ) THEN ! Check for threshold |
---|
| 244 | |
---|
| 245 | !PRINT *, ' [WL_COARE3P6] Tau_ac, Qnt_ac =', REAL(Tau_ac(ji,jj),4), REAL(Qnt_ac(ji,jj),4) |
---|
| 246 | |
---|
| 247 | !Tau_ac(ji,jj) = Tau_ac(ji,jj) + MAX(.002_wp , pTau(ji,jj))*rdt ! momentum integral |
---|
| 248 | ztac = Tau_ac(ji,jj) + MAX(.002_wp , pTau(ji,jj))*rdt ! updated momentum integral |
---|
| 249 | |
---|
| 250 | IF ( Qnt_ac(ji,jj) + zQabs*rdt > 0._wp ) THEN !check threshold for warm layer existence |
---|
| 251 | !****************************************** |
---|
| 252 | ! Compute the absorption profile |
---|
| 253 | !****************************************** |
---|
| 254 | DO jl = 1, 5 !loop 5 times for zfr |
---|
| 255 | zfr = 1. - ( 0.28*0.014*(1. - EXP(-zdz/0.014)) + 0.27*0.357*(1. - EXP(-zdz/0.357)) & |
---|
| 256 | & + 0.45*12.82*(1-EXP(-zdz/12.82)) ) / zdz |
---|
| 257 | zqac = Qnt_ac(ji,jj) + (zfr*pQsw(ji,jj) + pQnsol(ji,jj))*rdt ! updated heat absorbed |
---|
| 258 | IF ( zqac > 1._wp ) zdz = MAX( MIN( H_wl_max , zcd1*ztac/SQRT(zqac)) , 0.1_wp ) ! Warm-layer depth |
---|
| 259 | END DO |
---|
| 260 | |
---|
| 261 | ELSE |
---|
| 262 | !*********************** |
---|
| 263 | ! Warm layer wiped out |
---|
| 264 | !*********************** |
---|
| 265 | zfr = 0.75 |
---|
| 266 | zdz = H_wl_max |
---|
| 267 | zqac = Qnt_ac(ji,jj) + (zfr*pQsw(ji,jj) + pQnsol(ji,jj))*rdt ! updated heat absorbed |
---|
| 268 | |
---|
| 269 | END IF !IF ( Qnt_ac(ji,jj) + zQabs*rdt > 0._wp ) |
---|
| 270 | |
---|
| 271 | 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 ! |
---|
| 272 | |
---|
| 273 | ! Warm layer correction |
---|
| 274 | 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)<zdz (pdT(ji,jj) = zdT_wl*gdept_1d(1)/zdz) |
---|
| 275 | pdT(ji,jj) = zdT_wl * ( flg + (1._wp-flg)*gdept_1d(1)/zdz ) |
---|
| 276 | |
---|
| 277 | END IF ! IF ( zQabs >= Qabs_thr ) |
---|
| 278 | |
---|
| 279 | END IF ! IF ( isd_sol >= 21600 ) THEN ! (21600 == 6am) |
---|
| 280 | |
---|
| 281 | IF ( iwait == 0 ) THEN |
---|
| 282 | IF ( (zQabs >= Qabs_thr).AND.(rhr_sol >= 5._wp) ) THEN |
---|
| 283 | !PRINT *, ' [WL_COARE3P6] WE UPDATE ACCUMULATED FLUXES !!!' |
---|
| 284 | Qnt_ac(ji,jj) = zqac ! Updating Qnt_ac, heat integral |
---|
| 285 | Tau_ac(ji,jj) = ztac ! |
---|
| 286 | IF ( PRESENT(mask_wl) ) mask_wl(ji,jj) = 1 |
---|
| 287 | END IF |
---|
| 288 | END IF |
---|
| 289 | |
---|
| 290 | H_wl(ji,jj) = zdz |
---|
| 291 | |
---|
| 292 | IF ( PRESENT(Hwl) ) Hwl(ji,jj) = H_wl(ji,jj) |
---|
| 293 | |
---|
| 294 | END DO |
---|
| 295 | END DO |
---|
| 296 | |
---|
| 297 | END SUBROUTINE WL_COARE3P6 |
---|
| 298 | |
---|
| 299 | |
---|
| 300 | !!====================================================================== |
---|
| 301 | END MODULE sbcblk_skin_coare3p6 |
---|