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_coare3p6.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_coare3p6.F90 @ 11623

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

LB: syntax improved...
File size: 16.0 KB
Line 
1MODULE sbcblk_skin_coare3p6
2   !!======================================================================
3   !!                   ***  MODULE  sbcblk_skin_coare3p6  ***
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 lib_mpp         ! distribued memory computing library
29   USE in_out_manager  ! I/O manager
30   USE lib_fortran     ! to use key_nosignedzero
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC :: CS_COARE3P6, WL_COARE3P6
36
37   !! Cool-skin related parameters:
38   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: delta_vl  !: thickness of the surface viscous layer (in the water) right below the air-sea interface [m]
39
40   !! Warm-layer related parameters:
41   REAL(wp), PARAMETER, PUBLIC :: H_wl_max = 20._wp    !: maximum depth of warm layer (adjustable)
42   !
43   REAL(wp), PARAMETER :: rich   = 0.65_wp   !: critical Richardson number
44   !
45   REAL(wp), PARAMETER :: Qabs_thr = 50._wp  !: threshold for heat flux absorbed in WL
46   REAL(wp), PARAMETER :: zfr0   = 0.5_wp    !: initial value of solar flux absorption
47   !
48   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: Tau_ac  ! time integral / accumulated momentum Tauxdt => [N.s/m^2] (reset to zero every midnight)
49   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: Qnt_ac    ! time integral / accumulated heat stored by the warm layer Qxdt => [J/m^2] (reset to zero every midnight)
50   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: H_wl     ! depth of warm-layer [m]
51   !!----------------------------------------------------------------------
52CONTAINS
53
54
55   SUBROUTINE CS_COARE3P6( pQsw, pQnsol, pustar, pSST, pQlat,  pdT )
56      !!---------------------------------------------------------------------
57      !!
58      !!  Cool-Skin scheme according to Fairall et al. 1996, revisited for COARE 3.6 (Fairall et al., 2019)
59      !!     ------------------------------------------------------------------
60      !!
61      !!  **   INPUT:
62      !!     *pQsw*       surface net solar radiation into the ocean     [W/m^2] => >= 0 !
63      !!     *pQnsol*     surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 !
64      !!     *pustar*     friction velocity u*                           [m/s]
65      !!     *pSST*       bulk SST (taken at depth gdept_1d(1))          [K]
66      !!     *pQlat*      surface latent heat flux                       [K]
67      !!
68      !!  **  INPUT/OUTPUT:
69      !!     *pdT*  : as input =>  previous estimate of dT cool-skin
70      !!              as output =>  new estimate of dT cool-skin
71      !!
72      !!------------------------------------------------------------------
73      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)    :: pQsw   ! net solar a.k.a shortwave radiation into the ocean (after albedo) [W/m^2]
74      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)    :: pQnsol ! non-solar heat flux to the ocean [W/m^2]
75      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)    :: pustar  ! friction velocity, temperature and humidity (u*,t*,q*)
76      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)    :: pSST ! bulk SST [K]
77      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)    :: pQlat  ! latent heat flux [W/m^2]
78      !!
79      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pdT    ! dT due to cool-skin effect
80      !!---------------------------------------------------------------------
81      INTEGER  ::   ji, jj     ! dummy loop indices
82      REAL(wp) :: zQnet, zQnsol, zlamb, zdelta, zalpha_w, zfr, &
83         &        zz1, zz2, zus, &
84         &        ztf
85      !!---------------------------------------------------------------------
86
87      DO jj = 1, jpj
88         DO ji = 1, jpi
89
90            zalpha_w = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!)
91
92            zQnsol = MAX( 1._wp , - pQnsol(ji,jj) ) ! Non-solar heat loss to the atmosphere
93
94            zdelta = delta_vl(ji,jj)   ! using last value of delta
95
96            !! Fraction of the shortwave flux absorbed by the cool-skin sublayer:
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 ) ! 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
99            zQnet = MAX( 1._wp , zQnsol - zfr*pQsw(ji,jj) )  ! Total cooling at the interface
100
101            ztf = 0.5 + SIGN(0.5_wp, zQnet) ! Qt > 0 => cooling of the layer => ztf = 1
102            !                                 Qt < 0 => warming of the layer => ztf = 0
103
104            !! Term alpha*Qb (Qb is the virtual surface cooling inc. buoyancy effect of salinity due to evap):
105            zz1 = zalpha_w*zQnet - 0.026*MIN(pQlat(ji,jj),0._wp)*rCp0_w/rLevap  ! alpha*(Eq.8) == alpha*Qb "-" because Qlat < 0
106            !! LB: this terms only makes sense if > 0 i.e. in the cooling case
107            !! so similar to what's done in ECMWF:
108            zz1 = MAX(0._wp , zz1)    ! 1. instead of 0.1 though ZQ = MAX(1.0,-pQlw(ji,jj) - pQsen(ji,jj) - pQlat(ji,jj))
109
110            zus = MAX(pustar(ji,jj), 1.E-4_wp) ! Laurent: too low wind (u*) might cause problem in stable cases:
111            zz2 = zus*zus * roadrw
112            zz2 = zz2*zz2
113            zlamb =  6._wp*( 1._wp + (rcst_cs*zz1/zz2)**0.75 )**(-1./3.) ! Lambda (Eq.14) (Saunders)
114
115            ! Updating molecular sublayer thickness (delta):
116            zz2    = rnu0_w/(SQRT(roadrw)*zus)
117            zdelta =      ztf    *          zlamb*zz2   &  ! Eq.12 (when alpha*Qb>0 / cooling of layer)
118               &    + (1._wp - ztf) * MIN(0.007_wp , 6._wp*zz2 )    ! Eq.12 (when alpha*Qb<0 / warming of layer)
119            !LB: changed 0.01 to 0.007
120
121            !! Once again with the new zdelta:
122            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)
123            zQnet = MAX( 1._wp , zQnsol - zfr*pQsw(ji,jj) ) ! Total cooling at the interface
124
125            !! Update!
126            pdT(ji,jj) =  MIN( - zQnet*zdelta/rk0_w , 0._wp )   ! temperature increment
127            delta_vl(ji,jj) = zdelta
128
129         END DO
130      END DO
131
132   END SUBROUTINE CS_COARE3P6
133
134
135
136
137   SUBROUTINE WL_COARE3P6( kt,  pQsw, pQnsol, pTau, pSST, plon, isd, iwait,  pdT, &
138      &                    Hwl, mask_wl )
139      !!---------------------------------------------------------------------
140      !!
141      !!  Warm-Layer scheme according to COARE 3.6 (Fairall et al, 2019)
142      !!     ------------------------------------------------------------------
143      !!
144      !!  **   INPUT:
145      !!     *pQsw*       surface net solar radiation into the ocean     [W/m^2] => >= 0 !
146      !!     *pQnsol*     surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 !
147      !!     *pTau*       surface wind stress                            [N/m^2]
148      !!     *pSST*       bulk SST  (taken at depth gdept_1d(1))         [K]
149      !!     *plon*       longitude                                      [deg.E]
150      !!     *isd*        current UTC time, counted in second since 00h of the current day
151      !!     *iwait*      if /= 0 then wait before updating accumulated fluxes, we are within a converging itteration loop...
152      !!
153      !!  **   OUTPUT:
154      !!     *pdT*   dT due to warm-layer effect => difference between "almost surface (right below viscous layer) and depth of bulk SST
155      !!---------------------------------------------------------------------
156      !!
157      !!   ** OPTIONAL OUTPUT:
158      !!     *Hwl*        depth of warm layer [m]
159      !!     *mask_wl*    mask for possible existence of a warm-layer (1) or not (0)
160      !!
161      !!------------------------------------------------------------------
162      INTEGER ,                     INTENT(in)  :: kt
163      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pQsw     ! surface net solar radiation into the ocean [W/m^2]     => >= 0 !
164      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pQnsol   ! surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 !
165      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTau     ! wind stress [N/m^2]
166      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pSST     ! bulk SST at depth gdept_1d(1) [K]
167      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: plon     ! longitude [deg.E]
168      INTEGER ,                     INTENT(in)  :: isd      ! current UTC time, counted in second since 00h of the current day
169      INTEGER ,                     INTENT(in)  :: iwait    ! if /= 0 then wait before updating accumulated fluxes
170      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pdT      ! dT due to warm-layer effect => difference between "almost surface (right below viscous layer) and depth of bulk SST
171      !!
172      REAL(wp),   DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: Hwl     ! depth of warm layer [m]
173      INTEGER(1), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: mask_wl ! mask for possible existence of a warm-layer (1) or not (0)
174      !
175      !
176      INTEGER :: ji,jj
177      !
178      REAL(wp) :: zdT_wl, zQabs, zfr, zdz
179      REAL(wp) :: zqac, ztac
180      REAL(wp) :: zalpha_w, zcd1, zcd2, flg
181      CHARACTER(len=128) :: cf_tmp
182      !!---------------------------------------------------------------------
183
184      REAL(wp) :: rlag_gw_h, &  ! local solar time lag in hours   / Greenwich meridian (lon==0) => ex: ~ -10.47 hours for Hawai
185         &        rhr_sol       ! local solar time in hours since midnight
186
187      INTEGER  :: ilag_gw_s, &  ! local solar time LAG in seconds / Greenwich meridian (lon==0) => ex: ~ INT( -10.47*3600. ) seconds for Hawai
188         &        isd_sol,   &  ! local solar time in number of seconds since local solar midnight
189         &        jl
190
191      !! INITIALIZATION:
192      pdT(:,:) = 0._wp    ! dT initially set to 0._wp
193      zdT_wl = 0._wp       ! total warming (amplitude) in warm layer
194      zQabs  = 0._wp       ! total heat absorped in warm layer
195      zfr    = zfr0        ! initial value of solar flux absorption
196      ztac   = 0._wp
197      zqac   = 0._wp
198      IF ( PRESENT(mask_wl) ) mask_wl(:,:) = 0
199
200      DO jj = 1, jpj
201         DO ji = 1, jpi
202
203            zdz   = MAX( MIN(H_wl(ji,jj),H_wl_max) , 0.1_wp) ! depth of warm layer
204
205            !! Need to know the local solar time from longitude and isd:
206            rlag_gw_h = -1._wp * MODULO( ( 360._wp - MODULO(plon(ji,jj),360._wp) ) / 15._wp , 24._wp )
207            rlag_gw_h = -1._wp * SIGN( MIN(ABS(rlag_gw_h) , ABS(MODULO(rlag_gw_h,24._wp))), rlag_gw_h + 12._wp )
208            ilag_gw_s = INT( rlag_gw_h*3600._wp )
209            isd_sol = MODULO( isd + ilag_gw_s , 24*3600 )
210            rhr_sol = REAL( isd_sol , wp) / 3600._wp
211
212            !*****  variables for warm layer  ***
213            zalpha_w = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!)
214
215            zcd1 = SQRT(2._wp*rich*rCp0_w/(zalpha_w*grav*rho0_w))        !mess-o-constants 1
216            zcd2 = SQRT(2._wp*zalpha_w*grav/(rich*rho0_w))/(rCp0_w**1.5) !mess-o-constants 2
217
218            !********************************************************
219            !****  Compute apply warm layer  correction *************
220            !********************************************************
221           
222            !IF (isd_sol <= rdt ) THEN    !re-zero at midnight ! LOLO improve: risky if real midnight (00:00:00) is not a time in vtime...
223            IF ( (rhr_sol > 23.5_wp).OR.(rhr_sol < 4._wp) ) THEN
224               !PRINT *, '  [WL_COARE3P6] MIDNIGHT RESET !!!!, isd_sol =>', isd_sol
225               zdz           = H_wl_max
226               Tau_ac(ji,jj) = 0._wp
227               Qnt_ac(ji,jj) = 0._wp
228            END IF
229
230            IF ( rhr_sol > 5._wp ) THEN  ! ( 5am)
231
232               !PRINT *, '  [WL_COARE3P6] WE DO WL !'
233               !PRINT *, '  [WL_COARE3P6] isd_sol, pTau, pSST, pdT =', isd_sol, REAL(pTau(ji,jj),4), REAL(pSST(ji,jj),4), REAL(pdT(ji,jj),4)
234
235               !************************************
236               !****   set warm layer constants  ***
237               !************************************
238
239               zQabs = zfr*pQsw(ji,jj) + pQnsol(ji,jj)       ! tot heat absorbed in warm layer
240
241               !PRINT *, '  [WL_COARE3P6] rdt,  pQsw, pQnsol, zQabs =', rdt, REAL(pQsw(ji,jj),4), REAL(pQnsol(ji,jj),4), REAL(zQabs,4)
242
243               IF ( zQabs >= Qabs_thr ) THEN         ! Check for threshold
244
245                  !PRINT *, '  [WL_COARE3P6] Tau_ac, Qnt_ac =', REAL(Tau_ac(ji,jj),4), REAL(Qnt_ac(ji,jj),4)
246
247                  !Tau_ac(ji,jj) = Tau_ac(ji,jj) + MAX(.002_wp , pTau(ji,jj))*rdt      ! momentum integral
248                  ztac = Tau_ac(ji,jj) + MAX(.002_wp , pTau(ji,jj))*rdt      ! updated momentum integral
249
250                  IF ( Qnt_ac(ji,jj) + zQabs*rdt > 0._wp ) THEN         !check threshold for warm layer existence
251                     !******************************************
252                     ! Compute the absorption profile
253                     !******************************************
254                     DO jl = 1, 5                           !loop 5 times for zfr
255                        zfr = 1. - ( 0.28*0.014*(1. - EXP(-zdz/0.014)) + 0.27*0.357*(1. - EXP(-zdz/0.357)) &
256                           &        + 0.45*12.82*(1-EXP(-zdz/12.82)) ) / zdz
257                        zqac = Qnt_ac(ji,jj) + (zfr*pQsw(ji,jj) + pQnsol(ji,jj))*rdt ! updated heat absorbed
258                        IF ( zqac > 1._wp ) zdz = MAX( MIN( H_wl_max , zcd1*ztac/SQRT(zqac)) , 0.1_wp ) ! Warm-layer depth
259                     END DO
260
261                  ELSE
262                     !***********************
263                     ! Warm layer wiped out
264                     !***********************
265                     zfr  = 0.75
266                     zdz  = H_wl_max
267                     zqac = Qnt_ac(ji,jj) + (zfr*pQsw(ji,jj) + pQnsol(ji,jj))*rdt ! updated heat absorbed
268
269                  END IF !IF ( Qnt_ac(ji,jj) + zQabs*rdt > 0._wp )
270
271                  IF ( zqac > 1._wp ) zdT_wl = zcd2*zqac**1.5/ztac * MAX(zqac/ABS(zqac),0._wp)  !! => IF(zqac>0._wp): zdT_wl=zcd2*zqac**1.5/ztac ; ELSE: zdT_wl=0. / ! normally: zqac > 0 !
272                 
273                  ! Warm layer correction
274                  flg = 0.5_wp + SIGN( 0.5_wp , gdept_1d(1)-zdz )               ! => 1 when gdept_1d(1)>zdz (pdT(ji,jj) = zdT_wl) | 0 when gdept_1d(1)<zdz (pdT(ji,jj) = zdT_wl*gdept_1d(1)/zdz)
275                  pdT(ji,jj) = zdT_wl * ( flg + (1._wp-flg)*gdept_1d(1)/zdz )
276                 
277               END IF ! IF ( zQabs >= Qabs_thr )
278
279            END IF ! IF ( isd_sol >= 21600 ) THEN  ! (21600 == 6am)
280
281            IF ( iwait == 0 ) THEN
282               IF ( (zQabs >= Qabs_thr).AND.(rhr_sol >= 5._wp) ) THEN
283                  !PRINT *, '  [WL_COARE3P6] WE UPDATE ACCUMULATED FLUXES !!!'
284                  Qnt_ac(ji,jj) = zqac ! Updating Qnt_ac, heat integral
285                  Tau_ac(ji,jj) = ztac !
286                  IF ( PRESENT(mask_wl) ) mask_wl(ji,jj) = 1
287               END IF
288            END IF
289
290            H_wl(ji,jj) = zdz
291
292            IF ( PRESENT(Hwl) ) Hwl(ji,jj) = H_wl(ji,jj)
293
294         END DO
295      END DO
296     
297   END SUBROUTINE WL_COARE3P6
298
299
300   !!======================================================================
301END MODULE sbcblk_skin_coare3p6
Note: See TracBrowser for help on using the repository browser.