New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
sbcblk_skin_coare.F90 in NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/SBC – NEMO

source: NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/SBC/sbcblk_skin_coare.F90 @ 12406

Last change on this file since 12406 was 12406, checked in by davestorkey, 4 years ago
  1. Rename the namelist timestep parameter: rn_rdt -> rn_Dt
  2. Use the namelist parameter instead of the non-DOCTOR parameter: rdt -> rn_Dt

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