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_ecmwf.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_ecmwf.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: 14.9 KB
Line 
1MODULE sbcblk_skin_ecmwf
2   !!======================================================================
3   !!                   ***  MODULE  sbcblk_skin_ecmwf  ***
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 IFS of ECMWF (cycle 40r1)
8   !!   based on IFS doc (avaible online on the ECMWF's website)
9   !!
10   !!   From routine turb_ecmwf maintained and developed in AeroBulk
11   !!            (https://github.com/brodeau/aerobulk)
12   !!
13   !! ** Author: L. Brodeau, September 2019 / AeroBulk (https://github.com/brodeau/aerobulk)
14   !!----------------------------------------------------------------------
15   !! History :  4.0  !  2016-02  (L.Brodeau)   Original code
16   !!----------------------------------------------------------------------
17
18   !!----------------------------------------------------------------------
19   !!  cswl_ecmwf : computes the surface skin temperature (aka SSST)
20   !!               based on the cool-skin/warm-layer scheme used at ECMWF
21   !!----------------------------------------------------------------------
22   USE oce             ! ocean dynamics and tracers
23   USE dom_oce         ! ocean space and time domain
24   USE phycst          ! physical constants
25   USE sbc_oce         ! Surface boundary condition: ocean fields
26
27   USE sbcblk_phy      !LOLO: all thermodynamics functions, rho_air, q_sat, etc...
28
29   USE lib_mpp         ! distribued memory computing library
30   USE in_out_manager  ! I/O manager
31   USE lib_fortran     ! to use key_nosignedzero
32
33
34   IMPLICIT NONE
35   PRIVATE
36
37   PUBLIC :: CS_ECMWF, WL_ECMWF
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   !
50   REAL(wp), PARAMETER, PUBLIC :: rd0  = 3.    !: Depth scale [m] of warm layer, "d" in Eq.11 (Zeng & Beljaars 2005)
51   REAL(wp), PARAMETER         :: zRhoCp_w = rho0_w*rCp0_w
52   !
53   REAL(wp), PARAMETER         :: rNuwl0 = 0.5  !: Nu (exponent of temperature profile) Eq.11
54   !                                            !: (Zeng & Beljaars 2005) !: set to 0.5 instead of
55   !                                            !: 0.3 to respect a warming of +3 K in calm
56   !                                            !: condition for the insolation peak of +1000W/m^2
57   !!----------------------------------------------------------------------
58CONTAINS
59
60
61   SUBROUTINE CS_ECMWF( pQsw, pQnsol, pustar, pSST )
62      !!---------------------------------------------------------------------
63      !!
64      !! Cool-skin parameterization, based on Fairall et al., 1996:
65      !!
66      !! Fairall, C. W., Bradley, E. F., Godfrey, J. S., Wick, G. A.,
67      !! Edson, J. B., and Young, G. S. ( 1996), Cool‐skin and warm‐layer
68      !! effects on sea surface temperature, J. Geophys. Res., 101( C1), 1295-1308,
69      !! doi:10.1029/95JC03190.
70      !!
71      !!------------------------------------------------------------------
72      !!
73      !!  **   INPUT:
74      !!     *pQsw*       surface net solar radiation into the ocean     [W/m^2] => >= 0 !
75      !!     *pQnsol*     surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 !
76      !!     *pustar*     friction velocity u*                           [m/s]
77      !!     *pSST*       bulk SST (taken at depth gdept_1d(1))          [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      !!---------------------------------------------------------------------
84      INTEGER  :: ji, jj, jc
85      REAL(wp) :: zQabs, zdelta, zfr
86      !!---------------------------------------------------------------------
87      DO jj = 1, jpj
88         DO ji = 1, jpi
89
90            zQabs = pQnsol(ji,jj) ! first guess of heat flux absorbed within the viscous sublayer of thicknes delta,
91            !                     !   => we DO not miss a lot assuming 0 solar flux absorbed in the tiny layer of thicknes zdelta...
92
93            zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pustar(ji,jj) )
94
95            DO jc = 1, 4 ! because implicit in terms of zdelta...
96               zfr = MAX( 0.065_wp + 11._wp*zdelta - 6.6E-5_wp/zdelta*(1._wp - EXP(-zdelta/8.E-4_wp)) , 0.01_wp ) ! Solar absorption, Eq.(5) Zeng & Beljaars, 2005
97               !              =>  (WARNING: 0.065 rather than 0.137 in Fairal et al. 1996)
98               zQabs = pQnsol(ji,jj) + zfr*pQsw(ji,jj)
99               zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, 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_ECMWF
108
109
110
111   SUBROUTINE WL_ECMWF( pQsw, pQnsol, pustar, pSST,  pustk )
112      !!---------------------------------------------------------------------
113      !!
114      !!  Warm-Layer scheme according to Zeng & Beljaars, 2005 (GRL)
115      !!  " A prognostic scheme of sea surface skin temperature for modeling and data assimilation "
116      !!
117      !!  STIL NO PROGNOSTIC EQUATION FOR THE DEPTH OF THE WARM-LAYER!
118      !!
119      !!    As included in IFS Cy45r1   /  E.C.M.W.F.
120      !!     ------------------------------------------------------------------
121      !!
122      !!  **   INPUT:
123      !!     *pQsw*       surface net solar radiation into the ocean     [W/m^2] => >= 0 !
124      !!     *pQnsol*     surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 !
125      !!     *pustar*     friction velocity u*                           [m/s]
126      !!     *pSST*       bulk SST  (taken at depth gdept_1d(1))         [K]
127      !!------------------------------------------------------------------
128      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pQsw     ! surface net solar radiation into the ocean [W/m^2]     => >= 0 !
129      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pQnsol   ! surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 !
130      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pustar   ! friction velocity [m/s]
131      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pSST     ! bulk SST at depth gdept_1d(1) [K]
132      !!
133      REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) :: pustk ! surface Stokes velocity [m/s]
134      !
135      INTEGER :: ji, jj, jc
136      !
137      REAL(wp) :: &
138         & zHwl,    &  !: thickness of the warm-layer [m]
139         & ztcorr,  &  !: correction of dT w.r.t measurement depth of bulk SST (first T-point)
140         & zalpha, & !: thermal expansion coefficient of sea-water [1/K]
141         & zdTwl_b, zdTwl_n, & ! temp. diff. between "almost surface (right below viscous layer) and bottom of WL
142         & zfr, zeta, &
143         & zusw, zusw2, &
144         & zLa, zfLa, &
145         & flg, zwf, zQabs, &
146         & ZA, ZB, zL1, zL2, &
147         &  zcst0, zcst1, zcst2, zcst3
148      !
149      LOGICAL :: l_pustk_known
150      !!---------------------------------------------------------------------
151
152      l_pustk_known = .FALSE.
153      IF ( PRESENT(pustk) ) l_pustk_known = .TRUE.
154
155      DO jj = 1, jpj
156         DO ji = 1, jpi
157
158            zHwl = Hz_wl(ji,jj) ! first guess for warm-layer depth (and unique..., less advanced than COARE3p6 !)
159            ! it is = rd0 (3m) in default Zeng & Beljaars case...
160
161            !! Previous value of dT / warm-layer, adapted to depth:
162            flg = 0.5_wp + SIGN( 0.5_wp , gdept_1d(1)-zHwl )               ! => 1 when gdept_1d(1)>zHwl (dT_wl(ji,jj) = zdTwl) | 0 when z_s$
163            ztcorr = flg + (1._wp - flg)*gdept_1d(1)/zHwl
164            zdTwl_b = MAX ( dT_wl(ji,jj) / ztcorr , 0._wp )
165            ! zdTwl is the difference between "almost surface (right below viscous layer) and bottom of WL (here zHwl)
166            ! pdT         "                          "                                    and depth of bulk SST (here gdept_1d(1))!
167            !! => but of course in general the bulk SST is taken shallower than zHwl !!! So correction less pronounced!
168            !! => so here since pdT is difference between surface and gdept_1d(1), need to increase fof zdTwl !
169
170            zalpha = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!)
171
172
173            ! *** zfr = Fraction of solar radiation absorbed in warm layer (-)
174            zfr = 1._wp - 0.28_wp*EXP(-71.5_wp*zHwl) - 0.27_wp*EXP(-2.8_wp*zHwl) - 0.45_wp*EXP(-0.07_wp*zHwl)  !: Eq. 8.157
175
176            zQabs = zfr*pQsw(ji,jj) + pQnsol(ji,jj)       ! tot heat absorbed in warm layer
177
178            zusw  = MAX( pustar(ji,jj), 1.E-4_wp ) * sq_radrw    ! u* in the water
179            zusw2 = zusw*zusw
180
181            ! Langmuir:
182            IF ( l_pustk_known ) THEN
183               zLa = SQRT(zusw/MAX(pustk(ji,jj),1.E-6))
184            ELSE
185               zla = 0.3_wp
186            END IF
187            zfLa = MAX( zla**(-2._wp/3._wp) , 1._wp )   ! Eq.(6)
188
189            zwf = 0.5_wp + SIGN(0.5_wp, zQabs)  ! zQabs > 0. => 1.  / zQabs < 0. => 0.
190
191            zcst1 = vkarmn*grav*zalpha
192
193            ! 1/L when zQabs > 0 :
194            zL2 = zcst1*zQabs / (zRhoCp_w*zusw2*zusw)
195               
196            zcst2 = zcst1 / ( 5._wp*zHwl*zusw2 )  !OR: zcst2 = zcst1*rNuwl0 / ( 5._wp*zHwl*zusw2 ) ???
197
198            zcst0 = rdt * (rNuwl0 + 1._wp) / zHwl
199           
200            ZA = zcst0 * zQabs / ( rNuwl0 * zRhoCp_w )
201
202            zcst3 = -zcst0 * vkarmn * zusw * zfLa
203
204            !! Sorry about all these constants ( constant w.r.t zdTwl), it's for
205            !! the sake of optimizations... So all these operations are not done
206            !! over and over within the iteration loop...
207           
208            !! T R U L L Y   I M P L I C I T => needs itteration
209            !! => have to itterate just because the 1/(Monin-Obukhov length), zL1, uses zdTwl when zQabs < 0..
210            !!    (without this term otherwize the implicit analytical solution is straightforward...)
211            zdTwl_n = zdTwl_b
212            DO jc = 1, 10
213               
214               zdTwl_n = 0.5_wp * ( zdTwl_n + zdTwl_b ) ! semi implicit, for faster convergence
215               
216               ! 1/L when zdTwl > 0 .AND. zQabs < 0 :
217               zL1 =         SQRT( zdTwl_n * zcst2 ) ! / zusw !!! Or??? => vkarmn * SQRT( zdTwl_n*grav*zalpha/( 5._wp*zHwl ) ) / zusw
218               !zL1 = vkarmn*SQRT( zdTwl_n       *grav*zalpha        / ( 5._wp*zHwl ) ) / zusw   ! => vkarmn outside, not inside zcst1 (just for this particular line) ???
219               
220               ! Stability parameter (z/L):
221               zeta =  (1._wp - zwf) * zHwl*zL1   +   zwf * zHwl*zL2
222
223               ZB = zcst3 / PHI(zeta)
224
225               zdTwl_n = MAX ( zdTwl_b + ZA + ZB*zdTwl_n , 0._wp )  ! Eq.(6)
226
227            END DO
228           
229            !! Update:
230            dT_wl(ji,jj) = zdTwl_n * ztcorr
231           
232         END DO
233      END DO
234
235   END SUBROUTINE WL_ECMWF
236
237
238
239   FUNCTION delta_skin_layer( palpha, pQd, pustar_a )
240      !!---------------------------------------------------------------------
241      !! Computes the thickness (m) of the viscous skin layer.
242      !! Based on Fairall et al., 1996
243      !!
244      !! Fairall, C. W., Bradley, E. F., Godfrey, J. S., Wick, G. A.,
245      !! Edson, J. B., and Young, G. S. ( 1996), Cool‐skin and warm‐layer
246      !! effects on sea surface temperature, J. Geophys. Res., 101( C1), 1295-1308,
247      !! doi:10.1029/95JC03190.
248      !!
249      !! L. Brodeau, october 2019
250      !!---------------------------------------------------------------------
251      REAL(wp)                :: delta_skin_layer
252      REAL(wp), INTENT(in)    :: palpha   ! thermal expansion coefficient of sea-water (SST accurate enough!)
253      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
254      REAL(wp), INTENT(in)    :: pustar_a ! friction velocity in the air (u*) [m/s]
255      !!---------------------------------------------------------------------
256      REAL(wp) :: zusw, zusw2, zlamb, ztf, ztmp
257      !!---------------------------------------------------------------------
258      ztf = 0.5_wp + SIGN(0.5_wp, pQd)  ! Qabs < 0 => cooling of the viscous layer => ztf = 0 (regular case)
259      !                                 ! Qabs > 0 => warming of the viscous layer => ztf = 1 (ex: weak evaporation and strong positive sensible heat flux)
260      !
261      zusw  = MAX(pustar_a, 1.E-4_wp) * sq_radrw    ! u* in the water
262      zusw2 = zusw*zusw
263      !
264      zlamb = 6._wp*( 1._wp + MAX(palpha*zcon0/(zusw2*zusw2)*pQd, 0._wp)**0.75 )**(-1./3.) ! see Eq.(14) in Fairall et al., 1996
265      !  => zlamb is not used when Qd > 0, and since zcon0 < 0, we just use this "MAX" to prevent FPE errors (something_negative)**0.75
266      !
267      ztmp = rnu0_w/zusw
268      delta_skin_layer = (1._wp-ztf) *     zlamb*ztmp           &  ! regular case, Qd < 0, see Eq.(12) in Fairall et al., 1996
269         &               +   ztf     * MIN(6._wp*ztmp , 0.007_wp)  ! when Qd > 0
270   END FUNCTION delta_skin_layer
271
272
273   FUNCTION PHI( pzeta)
274      !!---------------------------------------------------------------------
275      !!
276      !! Takaya et al., 2010
277      !!  Eq.(5)
278      !! L. Brodeau, october 2019
279      !!---------------------------------------------------------------------
280      REAL(wp)                :: PHI
281      REAL(wp), INTENT(in)    :: pzeta    ! stability parameter
282      !!---------------------------------------------------------------------
283      REAL(wp) :: ztf, zzt2
284      !!---------------------------------------------------------------------
285      !
286      zzt2 = pzeta*pzeta
287      !
288      ztf = 0.5_wp + SIGN(0.5_wp, pzeta)  ! zeta > 0 => ztf = 1
289      !                                   ! zeta < 0 => ztf = 0
290      PHI =      ztf     * ( 1. + (5.*pzeta + 4.*zzt2)/(1. + 3.*pzeta + 0.25*zzt2) ) &   ! zeta > 0
291         &  + (1. - ztf) * 1./SQRT( 1. - 16.*(-ABS(pzeta)) )                             ! zeta < 0
292      !
293   END FUNCTION PHI
294
295   !!======================================================================
296END MODULE sbcblk_skin_ecmwf
Note: See TracBrowser for help on using the repository browser.