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

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

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