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 @ 11772

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

LB: solid updates+improvements of cool-skin/warm-layer capabilty of COARE and ECMWF bulk algorithms!
File size: 16.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 scheme according to 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      !!------------------------------------------------------------------
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, zdelta, zfr
90      !!---------------------------------------------------------------------
91      DO jj = 1, jpj
92         DO ji = 1, jpi
93
94            zQabs  = MIN( -0.1_wp , pQnsol(ji,jj) ) ! first guess, we do not miss a lot assuming 0 solar flux absorbed in the tiny layer of thicknes$
95            !                                       ! also, we ONLY consider when the viscous layer is loosing heat to the atmosphere, we only deal with cool-skin! => hence the "MIN( -0$
96
97            zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pQlat(ji,jj), pustar(ji,jj) )
98
99            DO jc = 1, 4 ! because implicit in terms of zdelta...
100               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 !
101               zQabs  = MIN( -0.1_wp , pQnsol(ji,jj) + zfr*pQsw(ji,jj) ) ! Total cooling at the interface
102               zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pQlat(ji,jj), pustar(ji,jj) )
103            END DO
104
105            dT_cs(ji,jj) = MIN( zQabs*zdelta/rk0_w , 0._wp )   ! temperature increment
106
107         END DO
108      END DO
109
110   END SUBROUTINE CS_COARE
111
112
113   SUBROUTINE WL_COARE( pQsw, pQnsol, pTau, pSST, iwait )
114      !!---------------------------------------------------------------------
115      !!
116      !!  Warm-Layer scheme according to COARE 3.6 (Fairall et al, 2019)
117      !!     ------------------------------------------------------------------
118      !!
119      !!  **   INPUT:
120      !!     *pQsw*       surface net solar radiation into the ocean     [W/m^2] => >= 0 !
121      !!     *pQnsol*     surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 !
122      !!     *pTau*       surface wind stress                            [N/m^2]
123      !!     *pSST*       bulk SST  (taken at depth gdept_1d(1))         [K]
124      !!     *iwait*      if /= 0 then wait before updating accumulated fluxes, we are within a converging itteration loop...
125      !!---------------------------------------------------------------------
126      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pQsw     ! surface net solar radiation into the ocean [W/m^2]     => >= 0 !
127      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pQnsol   ! surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 !
128      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTau     ! wind stress [N/m^2]
129      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pSST     ! bulk SST at depth gdept_1d(1) [K]
130      INTEGER ,                     INTENT(in)  :: iwait    ! if /= 0 then wait before updating accumulated fluxes
131      !!
132      INTEGER :: ji,jj
133      !
134      REAL(wp) :: zdTwl, zHwl, zQabs, zfr
135      REAL(wp) :: zqac, ztac
136      REAL(wp) :: zalpha, zcd1, zcd2, flg
137      !!---------------------------------------------------------------------
138
139      REAL(wp) :: ztime, znoon, zmidn
140      INTEGER  :: jl
141
142      LOGICAL :: l_exit, l_destroy_wl
143
144      !! INITIALIZATION:
145      zQabs  = 0._wp       ! total heat flux absorped in warm layer
146      zfr    = zfr0        ! initial value of solar flux absorption !LOLO: save it and use previous value !!!
147
148      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!
149
150      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 ...
151
152      IF (lwp) THEN
153         WRITE(numout,*) ''
154         WRITE(numout,*) 'LOLO: sbcblk_skin_coare => nsec_day, ztime =', nsec_day, ztime
155      END IF
156     
157      DO jj = 1, jpj
158         DO ji = 1, jpi
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            zalpha = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!)
171
172            zcd1 = SQRT(2._wp*rich*rCp0_w/(zalpha*grav*rho0_w))        !mess-o-constants 1
173            zcd2 = SQRT(2._wp*zalpha*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            END IF
185
186            IF ( .NOT. l_exit ) THEN
187               !! Initial test on initial guess of absorbed heat flux in warm-layer:
188               zfr = 1._wp - ( 0.28*0.014*(1. - EXP(-zHwl/0.014)) + 0.27*0.357*(1. - EXP(-zHwl/0.357)) &
189                  &        + 0.45*12.82*(1-EXP(-zHwl/12.82)) ) / zHwl
190               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!!!
191               PRINT *, '#LBD:  Initial Qsw & Qnsol:', NINT(pQsw(ji,jj)), NINT(pQnsol(ji,jj))
192               PRINT *, '#LBD:       =>Qabs:', zQabs,' zfr=', zfr
193
194               IF ( (ABS(zdTwl) < 1.E-6_wp) .AND. (zQabs <= 0._wp) ) THEN
195                  ! We have not started to build a WL yet (dT==0) and there's no way it can occur now
196                  ! since zQabs <= 0._wp
197                  ! => no need to go further
198                  PRINT *, '#LBD: we have not started to to build a WL yet (dT==0)'
199                  PRINT *, '#LBD: and theres no way it can occur now since zQabs=', zQabs
200                  PRINT *, '#LBD: => leaving without changing anything...'
201                  l_exit = .TRUE.
202               END IF
203
204            END IF
205
206            ! Okay test on updated absorbed flux:
207            !LOLO: remove??? has a strong influence !!!
208            IF ( (.NOT.(l_exit)) .AND. (Qnt_ac(ji,jj) + zQabs*rdt <= 0._wp) ) THEN
209               PRINT *, '#LBD: Oh boy! Next Qnt_ac looking weak! =>', Qnt_ac(ji,jj) + zQabs*rdt
210               PRINT *, '#LBD:  => time to destroy the warm-layer!'
211               l_exit       = .TRUE.
212               l_destroy_wl = .TRUE.
213            END IF
214
215
216            IF ( .NOT. l_exit) THEN
217
218               ! Two possibilities at this point:
219               ! 1/ A warm layer already exists (dT>0) but it is cooling down because Qabs<0
220               ! 2/ Regardless of WL formed (dT==0 or dT>0), we are in the process to initiate one or warm further it !
221
222               PRINT *, '#LBD:======================================================'
223               PRINT *, '#LBD: WL action makes sense now! => zQabs,dT_wl=', REAL(zQabs,4), REAL(zdTwl,4)
224               PRINT *, '#LBD:======================================================'
225               PRINT *, '#LBD: current values for Qac and Tac=', REAL(Qnt_ac(ji,jj),4), REAL(Tau_ac(ji,jj),4)
226
227               ztac = Tau_ac(ji,jj) + MAX(.002_wp , pTau(ji,jj))*rdt      ! updated momentum integral
228               PRINT *, '#LBD: updated value for Tac=',  REAL(ztac,4)
229
230               !! We update the value of absorbtion and zQabs:
231               !! some part is useless if Qsw=0 !!!
232               DO jl = 1, 5
233                  zfr = 1. - ( 0.28*0.014*(1. - EXP(-zHwl/0.014)) + 0.27*0.357*(1. - EXP(-zHwl/0.357)) &
234                     &        + 0.45*12.82*(1-EXP(-zHwl/12.82)) ) / zHwl
235                  zQabs = zfr*pQsw(ji,jj) + pQnsol(ji,jj)
236                  zqac  = Qnt_ac(ji,jj) + zQabs*rdt ! updated heat absorbed
237                  IF ( zqac <= 0._wp ) EXIT
238                  zHwl = MAX( MIN( Hwl_max , zcd1*ztac/SQRT(zqac)) , 0.1_wp ) ! Warm-layer depth
239               END DO
240               PRINT *, '#LBD: updated absorption and WL depth=',  REAL(zfr,4), REAL(zHwl,4)
241               PRINT *, '#LBD: updated value for Qabs=',  REAL(zQabs,4), 'W/m2'
242               PRINT *, '#LBD: updated value for Qac =',  REAL(zqac,4), 'J'
243
244               IF ( zqac <= 0._wp ) THEN
245                  l_destroy_wl = .TRUE.
246                  l_exit       = .TRUE.
247               ELSE
248                  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 !
249                  PRINT *, '#LBD: updated preliminary value for dT_wl=',  REAL(zdTwl,4)
250                  ! Warm layer correction
251                  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)
252                  zdTwl = zdTwl * ( flg + (1._wp-flg)*gdept_1d(1)/zHwl )
253               END IF
254
255            END IF !IF ( .NOT. l_exit)
256
257            IF ( l_destroy_wl ) THEN
258               zdTwl = 0._wp
259               zfr   = 0.75_wp
260               zHwl  = Hwl_max
261               zqac  = 0._wp
262               ztac  = 0._wp
263            END IF
264
265            PRINT *, '#LBD: exit values for Qac & Tac:', REAL(zqac,4), REAL(ztac,4)
266
267            IF ( iwait == 0 ) THEN
268               !! Iteration loop within bulk algo is over, time to update what needs to be updated:
269               dT_wl(ji,jj)  = zdTwl
270               Hz_wl(ji,jj)  = zHwl
271               PRINT *, '#LBD: FINAL EXIT values for dT_wl & Hz_wl:', REAL(dT_wl(ji,jj),4), REAL(Hz_wl(ji,jj),4)
272               Qnt_ac(ji,jj) = zqac ! Updating Qnt_ac, heat integral
273               Tau_ac(ji,jj) = ztac
274               PRINT *, '#LBD: FINAL EXIT values for Qac & Tac:', REAL(Qnt_ac(ji,jj),4), REAL(Tau_ac(ji,jj),4)
275               PRINT *, '#LBD'
276            END IF
277
278         END DO
279      END DO
280
281   END SUBROUTINE WL_COARE
282
283
284
285
286   FUNCTION delta_skin_layer( palpha, pQabs, pQlat, pustar_a )
287      !!---------------------------------------------------------------------
288      !! Computes the thickness (m) of the viscous skin layer.
289      !! Based on Fairall et al., 1996
290      !!
291      !! Fairall, C. W., Bradley, E. F., Godfrey, J. S., Wick, G. A.,
292      !! Edson, J. B., and Young, G. S. ( 1996), Cool‐skin and warm‐layer
293      !! effects on sea surface temperature, J. Geophys. Res., 101( C1), 1295-1308,
294      !! doi:10.1029/95JC03190.
295      !!
296      !! L. Brodeau, october 2019
297      !!---------------------------------------------------------------------
298      REAL(wp)                :: delta_skin_layer
299      REAL(wp), INTENT(in)    :: palpha   ! thermal expansion coefficient of sea-water (SST accurate enough!)
300      REAL(wp), INTENT(in)    :: pQabs    ! < 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
301      REAL(wp), INTENT(in)    :: pQlat    ! latent heat flux [W/m^2]
302      REAL(wp), INTENT(in)    :: pustar_a ! friction velocity in the air (u*) [m/s]
303      !!---------------------------------------------------------------------
304      REAL(wp) :: zusw, zusw2, zlamb, zQb
305      !!---------------------------------------------------------------------
306
307      zQb = pQabs + 0.026*MIN(pQlat,0._wp)*rCp0_w/rLevap/palpha   ! LOLO: Double check sign + division by palpha !!! units are okay!
308
309      zusw  = MAX(pustar_a, 1.E-4_wp) * sq_radrw    ! u* in the water
310      zusw2 = zusw*zusw
311
312      zlamb = 6._wp*( 1._wp + (palpha*zcon0/(zusw2*zusw2)*zQb)**0.75 )**(-1./3.) ! see eq.(14) in Fairall et al., 1996
313
314      delta_skin_layer = zlamb*rnu0_w/zusw
315
316   END FUNCTION delta_skin_layer
317
318   !!======================================================================
319END MODULE sbcblk_skin_coare
Note: See TracBrowser for help on using the repository browser.