Changeset 11963 for NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_coare.F90
- Timestamp:
- 2019-11-26T12:08:01+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_coare.F90
r11962 r11963 2 2 !!====================================================================== 3 3 !! *** MODULE sbcblk_skin_coare *** 4 !! Computes:5 !! * the surface skin temperature (aka SSST) based on the cool-skin/warm-layer6 !! scheme used at ECMWF7 !! Using formulation/param. of COARE 3.6 (Fairall et al., 2019)8 4 !! 9 !! From routine turb_ecmwf maintained and developed in AeroBulk10 !! (https://github.com/brodeau/aerobulk)5 !! Module that gathers the cool-skin and warm-layer parameterization used 6 !! in the COARE family of bulk parameterizations. 11 7 !! 12 !! ** Author: L. Brodeau, September 2019 / AeroBulk (https://github.com/brodeau/aerobulk) 8 !! Based on the last update for version COARE 3.6 (Fairall et al., 2019) 9 !! 10 !! Module 'sbcblk_skin_coare' also maintained and developed in AeroBulk (as 11 !! 'mod_skin_coare') 12 !! (https://github.com/brodeau/aerobulk) !! 13 !! ** Author: L. Brodeau, November 2019 / AeroBulk (https://github.com/brodeau/aerobulk) 13 14 !!---------------------------------------------------------------------- 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 15 !! History : 4.x ! 2019-11 (L.Brodeau) Original code 20 16 !!---------------------------------------------------------------------- 21 17 USE oce ! ocean dynamics and tracers … … 24 20 USE sbc_oce ! Surface boundary condition: ocean fields 25 21 26 USE sbcblk_phy ! LOLO: all thermodynamics functions, rho_air, q_sat, etc...27 28 USE sbcdcy ! LOLO: to know hour of dawn and dusk: rdawn_dcy and rdusk_dcy (needed in WL_COARE)29 22 USE sbcblk_phy ! misc. functions for marine ABL physics (rho_air, q_sat, bulk_formula, etc) 23 24 USE sbcdcy !#LB: to know hour of dawn and dusk: rdawn_dcy and rdusk_dcy (needed in WL_COARE) 25 30 26 USE lib_mpp ! distribued memory computing library 31 27 USE in_out_manager ! I/O manager … … 38 34 39 35 !! Cool-skin related parameters: 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) 36 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: dT_cs !: dT due to cool-skin effect 37 ! ! => temperature difference between air-sea interface (z=0) 38 ! ! and right below viscous layer (z=delta) 42 39 43 40 !! Warm-layer related parameters: 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) 41 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: dT_wl !: dT due to warm-layer effect 42 ! ! => difference between "almost surface (right below 43 ! ! viscous layer, z=delta) and depth of bulk SST (z=gdept_1d(1)) 44 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: Hz_wl !: depth (aka thickness) of warm-layer [m] 45 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: Qnt_ac !: time integral / accumulated heat stored by the warm layer 46 ! ! Qxdt => [J/m^2] (reset to zero every midnight) 47 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: Tau_ac !: time integral / accumulated momentum 48 ! ! Tauxdt => [N.s/m^2] (reset to zero every midnight) 49 49 50 50 REAL(wp), PARAMETER, PUBLIC :: Hwl_max = 20._wp !: maximum depth of warm layer (adjustable) … … 85 85 !!--------------------------------------------------------------------- 86 86 INTEGER :: ji, jj, jc 87 REAL(wp) :: zQabs, zd elta, zfr87 REAL(wp) :: zQabs, zdlt, zfr, zalfa, zqlat, zus 88 88 !!--------------------------------------------------------------------- 89 89 DO jj = 1, jpj … … 91 91 92 92 zQabs = pQnsol(ji,jj) ! first guess of heat flux absorbed within the viscous sublayer of thicknes delta, 93 ! ! => we DO not miss a lot assuming 0 solar flux absorbed in the tiny layer of thicknes zdelta... 94 95 zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pQlat(ji,jj), pustar(ji,jj) ) 96 97 DO jc = 1, 4 ! because implicit in terms of zdelta... 98 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 ! 93 ! ! => we DO not miss a lot assuming 0 solar flux absorbed in the tiny layer of thicknes zdlt... 94 95 zalfa = alpha_sw(pSST(ji,jj)) ! (crude) thermal expansion coefficient of sea-water [1/K] 96 zqlat = pQlat(ji,jj) 97 zus = pustar(ji,jj) 98 99 100 zdlt = delta_skin_layer( zalfa, zQabs, zqlat, zus ) 101 102 DO jc = 1, 4 ! because implicit in terms of zdlt... 103 zfr = MAX( 0.137_wp + 11._wp*zdlt & 104 & - 6.6E-5_wp/zdlt*(1._wp - EXP(-zdlt/8.E-4_wp)) & 105 & , 0.01_wp ) ! Solar absorption, Eq.16 (Fairall al. 1996b) 106 ! !LB: why 0.065 and not 0.137 like in the paper??? Beljaars & Zeng use 0.065, not 0.137 ! 99 107 zQabs = pQnsol(ji,jj) + zfr*pQsw(ji,jj) 100 zd elta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pQlat(ji,jj), pustar(ji,jj))108 zdlt = delta_skin_layer( zalfa, zQabs, zqlat, zus ) 101 109 END DO 102 110 103 dT_cs(ji,jj) = zQabs*zd elta/rk0_w ! temperature increment, yes dT_cs can actually > 0, if Qabs > 0 (rare but possible!)111 dT_cs(ji,jj) = zQabs*zdlt/rk0_w ! temperature increment, yes dT_cs can actually > 0, if Qabs > 0 (rare but possible!) 104 112 105 113 END DO … … 132 140 REAL(wp) :: zdTwl, zHwl, zQabs, zfr 133 141 REAL(wp) :: zqac, ztac 134 REAL(wp) :: zal pha, zcd1, zcd2, flg142 REAL(wp) :: zalfa, zcd1, zcd2, flg 135 143 !!--------------------------------------------------------------------- 136 144 … … 142 150 !! INITIALIZATION: 143 151 zQabs = 0._wp ! total heat flux absorped in warm layer 144 zfr = zfr0 ! initial value of solar flux absorption ! LOLO: save it and use previous value !!!152 zfr = zfr0 ! initial value of solar flux absorption !#LB: save it and use previous value !!! 145 153 146 154 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! … … 148 156 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 ... 149 157 150 IF (lwp) THEN151 WRITE(numout,*) ''152 WRITE(numout,*) 'LOLO: sbcblk_skin_coare => nsec_day, ztime =', nsec_day, ztime153 END IF154 155 158 DO jj = 1, jpj 156 159 DO ji = 1, jpi … … 166 169 167 170 !***** variables for warm layer *** 168 zal pha = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water(SST accurate enough!)169 170 zcd1 = SQRT(2._wp*rich*rCp0_w/(zal pha*grav*rho0_w)) !mess-o-constants 1171 zcd2 = SQRT(2._wp*zal pha*grav/(rich*rho0_w))/(rCp0_w**1.5) !mess-o-constants 2172 173 171 zalfa = alpha_sw( pSST(ji,jj) ) ! (crude) thermal expansion coefficient of sea-water [1/K] (SST accurate enough!) 172 173 zcd1 = SQRT(2._wp*rich*rCp0_w/(zalfa*grav*rho0_w)) !mess-o-constants 1 174 zcd2 = SQRT(2._wp*zalfa*grav/(rich*rho0_w))/(rCp0_w**1.5) !mess-o-constants 2 175 176 174 177 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 175 178 zmidn = MOD( znoon-0.5_wp , 1._wp ) 176 zmidn = MOD( zmidn + 0.125_wp , 1._wp ) ! 3 hours past the local midnight 177 178 IF 179 zmidn = MOD( zmidn + 0.125_wp , 1._wp ) ! 3 hours past the local midnight 180 181 IF( (ztime >= zmidn) .AND. (ztime < rdawn_dcy(ji,jj)) ) THEN 179 182 ! Dawn reset to 0! 180 183 l_exit = .TRUE. 181 184 l_destroy_wl = .TRUE. 182 END 183 184 IF 185 ENDIF 186 187 IF( .NOT. l_exit ) THEN 185 188 !! Initial test on initial guess of absorbed heat flux in warm-layer: 186 zfr = 1._wp - ( 0.28*0.014*(1. - EXP(-zHwl/0.014)) + 0.27*0.357*(1. - EXP(-zHwl/0.357)) & 187 & + 0.45*12.82*(1-EXP(-zHwl/12.82)) ) / zHwl 188 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!!! 189 190 IF ( (ABS(zdTwl) < 1.E-6_wp) .AND. (zQabs <= 0._wp) ) THEN 189 zQabs = frac_solar_abs(zHwl)*pQsw(ji,jj) + pQnsol(ji,jj) ! first guess of tot. heat flux absorbed in warm layer 190 ! ! => #LB: depends of zfr, which is wild guess... Wrong!!! 191 IF( (ABS(zdTwl) < 1.E-6_wp) .AND. (zQabs <= 0._wp) ) THEN 191 192 ! We have not started to build a WL yet (dT==0) and there's no way it can occur now 192 193 ! since zQabs <= 0._wp 193 194 ! => no need to go further 194 195 l_exit = .TRUE. 195 END 196 197 END 196 ENDIF 197 198 ENDIF 198 199 199 200 ! Okay test on updated absorbed flux: 200 ! LOLO: remove??? has a strong influence !!!201 IF ( (.NOT.(l_exit)) .AND.(Qnt_ac(ji,jj) + zQabs*rdt <= 0._wp) ) THEN201 !#LB: remove??? has a strong influence !!! 202 IF( (.NOT. l_exit).AND.(Qnt_ac(ji,jj) + zQabs*rdt <= 0._wp) ) THEN 202 203 l_exit = .TRUE. 203 204 l_destroy_wl = .TRUE. 204 END 205 206 207 IF 205 ENDIF 206 207 208 IF( .NOT. l_exit) THEN 208 209 209 210 ! Two possibilities at this point: … … 217 218 !! some part is useless if Qsw=0 !!! 218 219 DO jl = 1, 5 219 zfr = 1. - ( 0.28*0.014*(1. - EXP(-zHwl/0.014)) + 0.27*0.357*(1. - EXP(-zHwl/0.357)) & 220 & + 0.45*12.82*(1-EXP(-zHwl/12.82)) ) / zHwl 221 zQabs = zfr*pQsw(ji,jj) + pQnsol(ji,jj) 220 zQabs = frac_solar_abs(zHwl)*pQsw(ji,jj) + pQnsol(ji,jj) 222 221 zqac = Qnt_ac(ji,jj) + zQabs*rdt ! updated heat absorbed 223 IF 222 IF( zqac <= 0._wp ) EXIT 224 223 zHwl = MAX( MIN( Hwl_max , zcd1*ztac/SQRT(zqac)) , 0.1_wp ) ! Warm-layer depth 225 224 END DO 226 225 227 IF 226 IF( zqac <= 0._wp ) THEN 228 227 l_destroy_wl = .TRUE. 229 228 l_exit = .TRUE. … … 234 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) 235 234 zdTwl = zdTwl * ( flg + (1._wp-flg)*gdept_1d(1)/zHwl ) 236 END 237 238 END IF !IF( .NOT. l_exit)239 240 IF 235 ENDIF 236 237 ENDIF !IF( .NOT. l_exit) 238 239 IF( l_destroy_wl ) THEN 241 240 zdTwl = 0._wp 242 241 zfr = 0.75_wp … … 244 243 zqac = 0._wp 245 244 ztac = 0._wp 246 END 247 248 IF 245 ENDIF 246 247 IF( iwait == 0 ) THEN 249 248 !! Iteration loop within bulk algo is over, time to update what needs to be updated: 250 249 dT_wl(ji,jj) = zdTwl … … 252 251 Qnt_ac(ji,jj) = zqac ! Updating Qnt_ac, heat integral 253 252 Tau_ac(ji,jj) = ztac 254 END 253 ENDIF 255 254 256 255 END DO … … 276 275 REAL(wp) :: delta_skin_layer 277 276 REAL(wp), INTENT(in) :: palpha ! thermal expansion coefficient of sea-water (SST accurate enough!) 278 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 277 REAL(wp), INTENT(in) :: pQd ! < 0 !!! part of the net heat flux actually absorbed in the WL [W/m^2] 278 ! ! => term "Q + Rs*fs" in eq.6 of Fairall et al. 1996 279 279 REAL(wp), INTENT(in) :: pQlat ! latent heat flux [W/m^2] 280 280 REAL(wp), INTENT(in) :: pustar_a ! friction velocity in the air (u*) [m/s] … … 282 282 REAL(wp) :: zusw, zusw2, zlamb, zQd, ztf, ztmp 283 283 !!--------------------------------------------------------------------- 284 285 zQd = pQd + 0.026*MIN(pQlat,0._wp)*rCp0_w/rLevap/palpha ! LOLO: Double check sign + division by palpha !!! units are okay!284 285 zQd = pQd + 0.026*MIN(pQlat,0._wp)*rCp0_w/rLevap/palpha ! #LB: Double check sign + division by palpha !!! units are okay! 286 286 287 287 ztf = 0.5_wp + SIGN(0.5_wp, zQd) ! Qabs < 0 => cooling of the viscous layer => ztf = 0 (regular case) … … 299 299 END FUNCTION delta_skin_layer 300 300 301 302 FUNCTION frac_solar_abs( pHwl ) 303 !!--------------------------------------------------------------------- 304 !! Fraction of solar heat flux absorbed inside warm layer 305 !!--------------------------------------------------------------------- 306 REAL(wp) :: frac_solar_abs 307 REAL(wp), INTENT(in) :: pHwl ! thickness of warm-layer [m] 308 !!--------------------------------------------------------------------- 309 frac_solar_abs = 1._wp - ( 0.28*0.014 *(1._wp - EXP(-pHwl/0.014)) & 310 & + 0.27*0.357*(1._wp - EXP(-pHwl/0.357)) & 311 & + 0.45*12.82*(1-EXP(-pHwl/12.82)) ) / pHwl 312 END FUNCTION frac_solar_abs 313 301 314 !!====================================================================== 302 315 END MODULE sbcblk_skin_coare
Note: See TracChangeset
for help on using the changeset viewer.
