source: NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/SBC/sbcblk_skin_coare.F90 @ 12154

Last change on this file since 12154 was 12154, checked in by cetlod, 10 months ago

commit

File size: 16.0 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
35   !! Cool-skin related parameters:
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)
39
40   !! Warm-layer related parameters:
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
50   REAL(wp), PARAMETER, PUBLIC :: Hwl_max = 20._wp    !: maximum depth of warm layer (adjustable)
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   !!----------------------------------------------------------------------
57CONTAINS
58
59
60   SUBROUTINE CS_COARE( pQsw, pQnsol, pustar, pSST, pQlat )
61      !!---------------------------------------------------------------------
62      !!
63      !! Cool-skin parameterization, based on Fairall et al., 1996,
64      !! revisited for COARE 3.6 (Fairall et al., 2019)
65      !!
66      !! Fairall, C. W., Bradley, E. F., Godfrey, J. S., Wick, G. A.,
67      !! Edson, J. B., and Young, G. S. ( 1996), Cool‐skin and warm‐layer
68      !! effects on sea surface temperature, J. Geophys. Res., 101( C1), 1295-1308,
69      !! doi:10.1029/95JC03190.
70      !!
71      !!------------------------------------------------------------------
72      !!
73      !!  **   INPUT:
74      !!     *pQsw*       surface net solar radiation into the ocean     [W/m^2] => >= 0 !
75      !!     *pQnsol*     surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 !
76      !!     *pustar*     friction velocity u*                           [m/s]
77      !!     *pSST*       bulk SST (taken at depth gdept_1d(1))          [K]
78      !!     *pQlat*      surface latent heat flux                       [K]
79      !!------------------------------------------------------------------
80      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQsw   ! net solar a.k.a shortwave radiation into the ocean (after albedo) [W/m^2]
81      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQnsol ! non-solar heat flux to the ocean [W/m^2]
82      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pustar  ! friction velocity, temperature and humidity (u*,t*,q*)
83      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pSST ! bulk SST [K]
84      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQlat  ! latent heat flux [W/m^2]
85      !!---------------------------------------------------------------------
86      INTEGER  :: ji, jj, jc
87      REAL(wp) :: zQabs, zdlt, zfr, zalfa, zqlat, zus
88      !!---------------------------------------------------------------------
89      DO jj = 1, jpj
90         DO ji = 1, jpi
91
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 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 !
107               zQabs = pQnsol(ji,jj) + zfr*pQsw(ji,jj)
108               zdlt = delta_skin_layer( zalfa, zQabs, zqlat, zus )
109            END DO
110
111            dT_cs(ji,jj) = zQabs*zdlt/rk0_w   ! temperature increment, yes dT_cs can actually > 0, if Qabs > 0 (rare but possible!)
112
113         END DO
114      END DO
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 jj = 1, jpj
159         DO ji = 1, jpi
160
161            l_exit       = .FALSE.
162            l_destroy_wl = .FALSE.
163
164            zdTwl =  dT_wl(ji,jj)                          ! value of previous time step as first guess
165            zHwl  = MAX( MIN(Hz_wl(ji,jj),Hwl_max),0.1_wp) !   "                  "           "
166
167            zqac = Qnt_ac(ji,jj) ! previous time step Qnt_ac
168            ztac = Tau_ac(ji,jj)
169
170            !*****  variables for warm layer  ***
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
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
178            zmidn = MOD(              znoon-0.5_wp                 , 1._wp )
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
182               ! Dawn reset to 0!
183               l_exit       = .TRUE.
184               l_destroy_wl = .TRUE.
185            ENDIF
186
187            IF( .NOT. l_exit ) THEN
188               !! Initial test on initial guess of absorbed heat flux in warm-layer:
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
192                  ! We have not started to build a WL yet (dT==0) and there's no way it can occur now
193                  ! since zQabs <= 0._wp
194                  ! => no need to go further
195                  l_exit = .TRUE.
196               ENDIF
197
198            ENDIF
199
200            ! Okay test on updated absorbed flux:
201            !#LB: remove??? has a strong influence !!!
202            IF( (.NOT. l_exit).AND.(Qnt_ac(ji,jj) + zQabs*rdt <= 0._wp) ) THEN
203               l_exit       = .TRUE.
204               l_destroy_wl = .TRUE.
205            ENDIF
206
207
208            IF( .NOT. l_exit) THEN
209
210               ! Two possibilities at this point:
211               ! 1/ A warm layer already exists (dT>0) but it is cooling down because Qabs<0
212               ! 2/ Regardless of WL formed (dT==0 or dT>0), we are in the process to initiate one or warm further it !
213
214               ztac = Tau_ac(ji,jj) + MAX(.002_wp , pTau(ji,jj))*rdt      ! updated momentum integral
215               !PRINT *, '#LBD: updated value for Tac=',  REAL(ztac,4)
216
217               !! We update the value of absorbtion and zQabs:
218               !! some part is useless if Qsw=0 !!!
219               DO jl = 1, 5
220                  zQabs = frac_solar_abs(zHwl)*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
225
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 !
231                  !PRINT *, '#LBD: updated preliminary value for dT_wl=',  REAL(zdTwl,4)
232                  ! Warm layer correction
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               ENDIF
236
237            ENDIF !IF( .NOT. l_exit)
238
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
245            ENDIF
246
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            ENDIF
254
255         END DO
256      END DO
257
258   END SUBROUTINE WL_COARE
259
260
261
262
263   FUNCTION delta_skin_layer( palpha, pQd, pQlat, pustar_a )
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!)
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      REAL(wp), INTENT(in)    :: pQlat    ! latent heat flux [W/m^2]
280      REAL(wp), INTENT(in)    :: pustar_a ! friction velocity in the air (u*) [m/s]
281      !!---------------------------------------------------------------------
282      REAL(wp) :: zusw, zusw2, zlamb, zQd, ztf, ztmp
283      !!---------------------------------------------------------------------
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
287      ztf = 0.5_wp + SIGN(0.5_wp, zQd)  ! Qabs < 0 => cooling of the viscous layer => ztf = 0 (regular case)
288      !                                 ! Qabs > 0 => warming of the viscous layer => ztf = 1 (ex: weak evaporation and strong positive sensible heat flux)
289      !
290      zusw  = MAX(pustar_a, 1.E-4_wp) * sq_radrw    ! u* in the water
291      zusw2 = zusw*zusw
292      !
293      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
294      !  => 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
295      !
296      ztmp = rnu0_w/zusw
297      delta_skin_layer = (1._wp-ztf) *     zlamb*ztmp           &  ! regular case, Qd < 0, see Eq.(12) in Fairall et al., 1996
298         &               +   ztf     * MIN(6._wp*ztmp , 0.007_wp)  ! when Qd > 0
299   END FUNCTION delta_skin_layer
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
314   !!======================================================================
315END MODULE sbcblk_skin_coare
Note: See TracBrowser for help on using the repository browser.