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 warm-layer effect => difference between "almost surface (right below viscous layer) and depth of bulk SST |
---|
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 warm-layer effect => difference between "almost surface (right below viscous layer) and depth of bulk SST |
---|
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 |
---|