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/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC – NEMO

source: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_coare.F90 @ 11804

Last change on this file since 11804 was 11804, checked in by laurent, 5 years ago

Use of "TURB_FLUXES@…" inside sbcblk.F90, positive latent HF & positive cool-skin temperature now allowed!
File size: 15.6 KB
Line 
1MODULE sbcblk_skin_coare
2   !!======================================================================
3   !!                   ***  MODULE  sbcblk_skin_coare  ***
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 sbcdcy          !LOLO: to know hour of dawn and dusk: rdawn_dcy and rdusk_dcy (needed in WL_COARE)
29   
30   USE lib_mpp         ! distribued memory computing library
31   USE in_out_manager  ! I/O manager
32   USE lib_fortran     ! to use key_nosignedzero
33
34   IMPLICIT NONE
35   PRIVATE
36
37   PUBLIC :: CS_COARE, WL_COARE
38
39   !! 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)
42   REAL(wp), PARAMETER :: zcon0 = -16._wp * grav * rho0_w * rCp0_w * rnu0_w*rnu0_w*rnu0_w / ( rk0_w*rk0_w ) ! "-" because ocean convention: Qabs > 0 => gain of heat for ocean!
43   !!                             => see eq.(14) in Fairall et al. 1996   (eq.(6) of Zeng aand Beljaars is WRONG! (typo?)
44
45   !! Warm-layer related parameters:
46   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: &
47      &                        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))
48      &                        Hz_wl,     &  !: depth of warm-layer [m]
49      &                        Qnt_ac,    &  !: time integral / accumulated heat stored by the warm layer Qxdt => [J/m^2] (reset to zero every midnight)
50      &                        Tau_ac        !: time integral / accumulated momentum 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, revisited for COARE 3.6 (Fairall et al., 2019)
66      !!
67      !! Fairall, C. W., Bradley, E. F., Godfrey, J. S., Wick, G. A.,
68      !! Edson, J. B., and Young, G. S. ( 1996), Cool‐skin and warm‐layer
69      !! effects on sea surface temperature, J. Geophys. Res., 101( C1), 1295-1308,
70      !! doi:10.1029/95JC03190.
71      !!
72      !!------------------------------------------------------------------
73      !!
74      !!  **   INPUT:
75      !!     *pQsw*       surface net solar radiation into the ocean     [W/m^2] => >= 0 !
76      !!     *pQnsol*     surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 !
77      !!     *pustar*     friction velocity u*                           [m/s]
78      !!     *pSST*       bulk SST (taken at depth gdept_1d(1))          [K]
79      !!     *pQlat*      surface latent heat flux                       [K]
80      !!------------------------------------------------------------------
81      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQsw   ! net solar a.k.a shortwave radiation into the ocean (after albedo) [W/m^2]
82      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQnsol ! non-solar heat flux to the ocean [W/m^2]
83      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pustar  ! friction velocity, temperature and humidity (u*,t*,q*)
84      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pSST ! bulk SST [K]
85      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQlat  ! latent heat flux [W/m^2]
86      !!---------------------------------------------------------------------
87      INTEGER  :: ji, jj, jc
88      REAL(wp) :: zQabs, zdelta, zfr
89      !!---------------------------------------------------------------------
90      DO jj = 1, jpj
91         DO ji = 1, jpi
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 zdelta...
95
96            zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pQlat(ji,jj), pustar(ji,jj) )
97
98            DO jc = 1, 4 ! because implicit in terms of zdelta...
99               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 !
100               zQabs = pQnsol(ji,jj) + zfr*pQsw(ji,jj)
101               zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pQlat(ji,jj), pustar(ji,jj) )
102            END DO
103
104            dT_cs(ji,jj) = zQabs*zdelta/rk0_w   ! temperature increment, yes dT_cs can actually > 0, if Qabs > 0 (rare but possible!)
105
106         END DO
107      END DO
108
109   END SUBROUTINE CS_COARE
110
111
112   SUBROUTINE WL_COARE( pQsw, pQnsol, pTau, pSST, iwait )
113      !!---------------------------------------------------------------------
114      !!
115      !!  Warm-Layer scheme according to COARE 3.6 (Fairall et al, 2019)
116      !!     ------------------------------------------------------------------
117      !!
118      !!  **   INPUT:
119      !!     *pQsw*       surface net solar radiation into the ocean     [W/m^2] => >= 0 !
120      !!     *pQnsol*     surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 !
121      !!     *pTau*       surface wind stress                            [N/m^2]
122      !!     *pSST*       bulk SST  (taken at depth gdept_1d(1))         [K]
123      !!     *iwait*      if /= 0 then wait before updating accumulated fluxes, we are within a converging itteration loop...
124      !!---------------------------------------------------------------------
125      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pQsw     ! surface net solar radiation into the ocean [W/m^2]     => >= 0 !
126      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pQnsol   ! surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 !
127      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTau     ! wind stress [N/m^2]
128      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pSST     ! bulk SST at depth gdept_1d(1) [K]
129      INTEGER ,                     INTENT(in)  :: iwait    ! if /= 0 then wait before updating accumulated fluxes
130      !!
131      INTEGER :: ji,jj
132      !
133      REAL(wp) :: zdTwl, zHwl, zQabs, zfr
134      REAL(wp) :: zqac, ztac
135      REAL(wp) :: zalpha, zcd1, zcd2, flg
136      !!---------------------------------------------------------------------
137
138      REAL(wp) :: ztime, znoon, zmidn
139      INTEGER  :: jl
140
141      LOGICAL :: l_exit, l_destroy_wl
142
143      !! INITIALIZATION:
144      zQabs  = 0._wp       ! total heat flux absorped in warm layer
145      zfr    = zfr0        ! initial value of solar flux absorption !LOLO: save it and use previous value !!!
146
147      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
149      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 ...
150
151      IF (lwp) THEN
152         WRITE(numout,*) ''
153         WRITE(numout,*) 'LOLO: sbcblk_skin_coare => nsec_day, ztime =', nsec_day, ztime
154      END IF
155     
156      DO jj = 1, jpj
157         DO ji = 1, jpi
158
159            l_exit       = .FALSE.
160            l_destroy_wl = .FALSE.
161
162            zdTwl =  dT_wl(ji,jj)                          ! value of previous time step as first guess
163            zHwl  = MAX( MIN(Hz_wl(ji,jj),Hwl_max),0.1_wp) !   "                  "           "
164
165            zqac = Qnt_ac(ji,jj) ! previous time step Qnt_ac
166            ztac = Tau_ac(ji,jj)
167
168            !*****  variables for warm layer  ***
169            zalpha = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!)
170
171            zcd1 = SQRT(2._wp*rich*rCp0_w/(zalpha*grav*rho0_w))        !mess-o-constants 1
172            zcd2 = SQRT(2._wp*zalpha*grav/(rich*rho0_w))/(rCp0_w**1.5) !mess-o-constants 2
173
174           
175            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
176            zmidn = MOD(              znoon-0.5_wp                 , 1._wp )
177            zmidn = MOD( zmidn + 0.125_wp , 1._wp ) ! 3 hours past the local midnight
178
179            IF ( (ztime >= zmidn) .AND. (ztime < rdawn_dcy(ji,jj)) ) THEN
180               ! Dawn reset to 0!
181               l_exit       = .TRUE.
182               l_destroy_wl = .TRUE.
183            END IF
184
185            IF ( .NOT. l_exit ) THEN
186               !! Initial test on initial guess of absorbed heat flux in warm-layer:
187               zfr = 1._wp - ( 0.28*0.014*(1. - EXP(-zHwl/0.014)) + 0.27*0.357*(1. - EXP(-zHwl/0.357)) &
188                  &        + 0.45*12.82*(1-EXP(-zHwl/12.82)) ) / zHwl
189               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!!!
190
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               END IF
197
198            END IF
199
200            ! Okay test on updated absorbed flux:
201            !LOLO: 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            END IF
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                  zfr = 1. - ( 0.28*0.014*(1. - EXP(-zHwl/0.014)) + 0.27*0.357*(1. - EXP(-zHwl/0.357)) &
221                     &        + 0.45*12.82*(1-EXP(-zHwl/12.82)) ) / zHwl
222                  zQabs = zfr*pQsw(ji,jj) + pQnsol(ji,jj)
223                  zqac  = Qnt_ac(ji,jj) + zQabs*rdt ! updated heat absorbed
224                  IF ( zqac <= 0._wp ) EXIT
225                  zHwl = MAX( MIN( Hwl_max , zcd1*ztac/SQRT(zqac)) , 0.1_wp ) ! Warm-layer depth
226               END DO
227
228               IF ( zqac <= 0._wp ) THEN
229                  l_destroy_wl = .TRUE.
230                  l_exit       = .TRUE.
231               ELSE
232                  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 !
233                  !PRINT *, '#LBD: updated preliminary value for dT_wl=',  REAL(zdTwl,4)
234                  ! Warm layer correction
235                  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)
236                  zdTwl = zdTwl * ( flg + (1._wp-flg)*gdept_1d(1)/zHwl )
237               END IF
238
239            END IF !IF ( .NOT. l_exit)
240
241            IF ( l_destroy_wl ) THEN
242               zdTwl = 0._wp
243               zfr   = 0.75_wp
244               zHwl  = Hwl_max
245               zqac  = 0._wp
246               ztac  = 0._wp
247            END IF
248
249            IF ( iwait == 0 ) THEN
250               !! Iteration loop within bulk algo is over, time to update what needs to be updated:
251               dT_wl(ji,jj)  = zdTwl
252               Hz_wl(ji,jj)  = zHwl
253               Qnt_ac(ji,jj) = zqac ! Updating Qnt_ac, heat integral
254               Tau_ac(ji,jj) = ztac
255            END IF
256
257         END DO
258      END DO
259
260   END SUBROUTINE WL_COARE
261
262
263
264
265   FUNCTION delta_skin_layer( palpha, pQd, pQlat, pustar_a )
266      !!---------------------------------------------------------------------
267      !! Computes the thickness (m) of the viscous skin layer.
268      !! Based on Fairall et al., 1996
269      !!
270      !! Fairall, C. W., Bradley, E. F., Godfrey, J. S., Wick, G. A.,
271      !! Edson, J. B., and Young, G. S. ( 1996), Cool‐skin and warm‐layer
272      !! effects on sea surface temperature, J. Geophys. Res., 101( C1), 1295-1308,
273      !! doi:10.1029/95JC03190.
274      !!
275      !! L. Brodeau, october 2019
276      !!---------------------------------------------------------------------
277      REAL(wp)                :: delta_skin_layer
278      REAL(wp), INTENT(in)    :: palpha   ! thermal expansion coefficient of sea-water (SST accurate enough!)
279      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
280      REAL(wp), INTENT(in)    :: pQlat    ! latent heat flux [W/m^2]
281      REAL(wp), INTENT(in)    :: pustar_a ! friction velocity in the air (u*) [m/s]
282      !!---------------------------------------------------------------------
283      REAL(wp) :: zusw, zusw2, zlamb, zQd, ztf, ztmp
284      !!---------------------------------------------------------------------
285     
286      zQd = pQd + 0.026*MIN(pQlat,0._wp)*rCp0_w/rLevap/palpha   ! LOLO: Double check sign + division by palpha !!! units are okay!
287
288      ztf = 0.5_wp + SIGN(0.5_wp, zQd)  ! Qabs < 0 => cooling of the viscous layer => ztf = 0 (regular case)
289      !                                 ! Qabs > 0 => warming of the viscous layer => ztf = 1 (ex: weak evaporation and strong positive sensible heat flux)
290      !
291      zusw  = MAX(pustar_a, 1.E-4_wp) * sq_radrw    ! u* in the water
292      zusw2 = zusw*zusw
293      !
294      zlamb = 6._wp*( 1._wp + MAX(palpha*zcon0/(zusw2*zusw2)*zQd, 0._wp)**0.75 )**(-1./3.) ! see Eq.(14) in Fairall et al., 1996
295      !  => zlamb is not used when Qd > 0, and since zcon0 < 0, we just use this "MAX" to prevent FPE errors (something_negative)**0.75
296      !
297      ztmp = rnu0_w/zusw
298      delta_skin_layer = (1._wp-ztf) *     zlamb*ztmp           &  ! regular case, Qd < 0, see Eq.(12) in Fairall et al., 1996
299         &               +   ztf     * MIN(6._wp*ztmp , 0.007_wp)  ! when Qd > 0
300   END FUNCTION delta_skin_layer
301
302   !!======================================================================
303END MODULE sbcblk_skin_coare
Note: See TracBrowser for help on using the repository browser.