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

Last change on this file since 11826 was 11826, checked in by laurent, 11 months ago

SBCBLK-related physical constants all moved into "sbcblk_phy.F90".

File size: 15.3 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
43   !! Warm-layer related parameters:
44   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: &
45      &                        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))
46      &                        Hz_wl,     &  !: depth of warm-layer [m]
47      &                        Qnt_ac,    &  !: time integral / accumulated heat stored by the warm layer Qxdt => [J/m^2] (reset to zero every midnight)
48      &                        Tau_ac        !: time integral / accumulated momentum 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, revisited for COARE 3.6 (Fairall et al., 2019)
64      !!
65      !! Fairall, C. W., Bradley, E. F., Godfrey, J. S., Wick, G. A.,
66      !! Edson, J. B., and Young, G. S. ( 1996), Cool‐skin and warm‐layer
67      !! effects on sea surface temperature, J. Geophys. Res., 101( C1), 1295-1308,
68      !! doi:10.1029/95JC03190.
69      !!
70      !!------------------------------------------------------------------
71      !!
72      !!  **   INPUT:
73      !!     *pQsw*       surface net solar radiation into the ocean     [W/m^2] => >= 0 !
74      !!     *pQnsol*     surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 !
75      !!     *pustar*     friction velocity u*                           [m/s]
76      !!     *pSST*       bulk SST (taken at depth gdept_1d(1))          [K]
77      !!     *pQlat*      surface latent heat flux                       [K]
78      !!------------------------------------------------------------------
79      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQsw   ! net solar a.k.a shortwave radiation into the ocean (after albedo) [W/m^2]
80      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQnsol ! non-solar heat flux to the ocean [W/m^2]
81      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pustar  ! friction velocity, temperature and humidity (u*,t*,q*)
82      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pSST ! bulk SST [K]
83      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQlat  ! latent heat flux [W/m^2]
84      !!---------------------------------------------------------------------
85      INTEGER  :: ji, jj, jc
86      REAL(wp) :: zQabs, zdelta, zfr
87      !!---------------------------------------------------------------------
88      DO jj = 1, jpj
89         DO ji = 1, jpi
90
91            zQabs = pQnsol(ji,jj) ! first guess of heat flux absorbed within the viscous sublayer of thicknes delta,
92            !                     !   => we DO not miss a lot assuming 0 solar flux absorbed in the tiny layer of thicknes zdelta...
93
94            zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pQlat(ji,jj), pustar(ji,jj) )
95
96            DO jc = 1, 4 ! because implicit in terms of zdelta...
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 ) ! 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 !
98               zQabs = pQnsol(ji,jj) + zfr*pQsw(ji,jj)
99               zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pQlat(ji,jj), pustar(ji,jj) )
100            END DO
101
102            dT_cs(ji,jj) = zQabs*zdelta/rk0_w   ! temperature increment, yes dT_cs can actually > 0, if Qabs > 0 (rare but possible!)
103
104         END DO
105      END DO
106
107   END SUBROUTINE CS_COARE
108
109
110   SUBROUTINE WL_COARE( pQsw, pQnsol, pTau, pSST, iwait )
111      !!---------------------------------------------------------------------
112      !!
113      !!  Warm-Layer scheme according to COARE 3.6 (Fairall et al, 2019)
114      !!     ------------------------------------------------------------------
115      !!
116      !!  **   INPUT:
117      !!     *pQsw*       surface net solar radiation into the ocean     [W/m^2] => >= 0 !
118      !!     *pQnsol*     surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 !
119      !!     *pTau*       surface wind stress                            [N/m^2]
120      !!     *pSST*       bulk SST  (taken at depth gdept_1d(1))         [K]
121      !!     *iwait*      if /= 0 then wait before updating accumulated fluxes, we are within a converging itteration loop...
122      !!---------------------------------------------------------------------
123      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pQsw     ! surface net solar radiation into the ocean [W/m^2]     => >= 0 !
124      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pQnsol   ! surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 !
125      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTau     ! wind stress [N/m^2]
126      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pSST     ! bulk SST at depth gdept_1d(1) [K]
127      INTEGER ,                     INTENT(in)  :: iwait    ! if /= 0 then wait before updating accumulated fluxes
128      !!
129      INTEGER :: ji,jj
130      !
131      REAL(wp) :: zdTwl, zHwl, zQabs, zfr
132      REAL(wp) :: zqac, ztac
133      REAL(wp) :: zalpha, zcd1, zcd2, flg
134      !!---------------------------------------------------------------------
135
136      REAL(wp) :: ztime, znoon, zmidn
137      INTEGER  :: jl
138
139      LOGICAL :: l_exit, l_destroy_wl
140
141      !! INITIALIZATION:
142      zQabs  = 0._wp       ! total heat flux absorped in warm layer
143      zfr    = zfr0        ! initial value of solar flux absorption !LOLO: save it and use previous value !!!
144
145      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!
146
147      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 ...
148
149      IF (lwp) THEN
150         WRITE(numout,*) ''
151         WRITE(numout,*) 'LOLO: sbcblk_skin_coare => nsec_day, ztime =', nsec_day, ztime
152      END IF
153     
154      DO jj = 1, jpj
155         DO ji = 1, jpi
156
157            l_exit       = .FALSE.
158            l_destroy_wl = .FALSE.
159
160            zdTwl =  dT_wl(ji,jj)                          ! value of previous time step as first guess
161            zHwl  = MAX( MIN(Hz_wl(ji,jj),Hwl_max),0.1_wp) !   "                  "           "
162
163            zqac = Qnt_ac(ji,jj) ! previous time step Qnt_ac
164            ztac = Tau_ac(ji,jj)
165
166            !*****  variables for warm layer  ***
167            zalpha = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!)
168
169            zcd1 = SQRT(2._wp*rich*rCp0_w/(zalpha*grav*rho0_w))        !mess-o-constants 1
170            zcd2 = SQRT(2._wp*zalpha*grav/(rich*rho0_w))/(rCp0_w**1.5) !mess-o-constants 2
171
172           
173            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
174            zmidn = MOD(              znoon-0.5_wp                 , 1._wp )
175            zmidn = MOD( zmidn + 0.125_wp , 1._wp ) ! 3 hours past the local midnight
176
177            IF ( (ztime >= zmidn) .AND. (ztime < rdawn_dcy(ji,jj)) ) THEN
178               ! Dawn reset to 0!
179               l_exit       = .TRUE.
180               l_destroy_wl = .TRUE.
181            END IF
182
183            IF ( .NOT. l_exit ) THEN
184               !! Initial test on initial guess of absorbed heat flux in warm-layer:
185               zfr = 1._wp - ( 0.28*0.014*(1. - EXP(-zHwl/0.014)) + 0.27*0.357*(1. - EXP(-zHwl/0.357)) &
186                  &        + 0.45*12.82*(1-EXP(-zHwl/12.82)) ) / zHwl
187               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!!!
188
189               IF ( (ABS(zdTwl) < 1.E-6_wp) .AND. (zQabs <= 0._wp) ) THEN
190                  ! We have not started to build a WL yet (dT==0) and there's no way it can occur now
191                  ! since zQabs <= 0._wp
192                  ! => no need to go further
193                  l_exit = .TRUE.
194               END IF
195
196            END IF
197
198            ! Okay test on updated absorbed flux:
199            !LOLO: remove??? has a strong influence !!!
200            IF ( (.NOT.(l_exit)) .AND. (Qnt_ac(ji,jj) + zQabs*rdt <= 0._wp) ) THEN
201               l_exit       = .TRUE.
202               l_destroy_wl = .TRUE.
203            END IF
204
205
206            IF ( .NOT. l_exit) THEN
207
208               ! Two possibilities at this point:
209               ! 1/ A warm layer already exists (dT>0) but it is cooling down because Qabs<0
210               ! 2/ Regardless of WL formed (dT==0 or dT>0), we are in the process to initiate one or warm further it !
211
212               ztac = Tau_ac(ji,jj) + MAX(.002_wp , pTau(ji,jj))*rdt      ! updated momentum integral
213               !PRINT *, '#LBD: updated value for Tac=',  REAL(ztac,4)
214
215               !! We update the value of absorbtion and zQabs:
216               !! some part is useless if Qsw=0 !!!
217               DO jl = 1, 5
218                  zfr = 1. - ( 0.28*0.014*(1. - EXP(-zHwl/0.014)) + 0.27*0.357*(1. - EXP(-zHwl/0.357)) &
219                     &        + 0.45*12.82*(1-EXP(-zHwl/12.82)) ) / zHwl
220                  zQabs = zfr*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               END IF
236
237            END IF !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            END IF
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            END IF
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] => term "Q + Rs*fs" in eq.6 of Fairall et al. 1996
278      REAL(wp), INTENT(in)    :: pQlat    ! latent heat flux [W/m^2]
279      REAL(wp), INTENT(in)    :: pustar_a ! friction velocity in the air (u*) [m/s]
280      !!---------------------------------------------------------------------
281      REAL(wp) :: zusw, zusw2, zlamb, zQd, ztf, ztmp
282      !!---------------------------------------------------------------------
283     
284      zQd = pQd + 0.026*MIN(pQlat,0._wp)*rCp0_w/rLevap/palpha   ! LOLO: Double check sign + division by palpha !!! units are okay!
285
286      ztf = 0.5_wp + SIGN(0.5_wp, zQd)  ! Qabs < 0 => cooling of the viscous layer => ztf = 0 (regular case)
287      !                                 ! Qabs > 0 => warming of the viscous layer => ztf = 1 (ex: weak evaporation and strong positive sensible heat flux)
288      !
289      zusw  = MAX(pustar_a, 1.E-4_wp) * sq_radrw    ! u* in the water
290      zusw2 = zusw*zusw
291      !
292      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
293      !  => 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
294      !
295      ztmp = rnu0_w/zusw
296      delta_skin_layer = (1._wp-ztf) *     zlamb*ztmp           &  ! regular case, Qd < 0, see Eq.(12) in Fairall et al., 1996
297         &               +   ztf     * MIN(6._wp*ztmp , 0.007_wp)  ! when Qd > 0
298   END FUNCTION delta_skin_layer
299
300   !!======================================================================
301END MODULE sbcblk_skin_coare
Note: See TracBrowser for help on using the repository browser.