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

Last change on this file since 11845 was 11826, checked in by laurent, 4 years ago

SBCBLK-related physical constants all moved into "sbcblk_phy.F90".
File size: 14.6 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
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   !
48   REAL(wp), PARAMETER, PUBLIC :: rd0  = 3.    !: Depth scale [m] of warm layer, "d" in Eq.11 (Zeng & Beljaars 2005)
49   REAL(wp), PARAMETER         :: zRhoCp_w = rho0_w*rCp0_w
50   !
51   REAL(wp), PARAMETER         :: rNuwl0 = 0.5  !: Nu (exponent of temperature profile) Eq.11
52   !                                            !: (Zeng & Beljaars 2005) !: set to 0.5 instead of
53   !                                            !: 0.3 to respect a warming of +3 K in calm
54   !                                            !: condition for the insolation peak of +1000W/m^2
55   !!----------------------------------------------------------------------
56CONTAINS
57
58
59   SUBROUTINE CS_ECMWF( pQsw, pQnsol, pustar, pSST )
60      !!---------------------------------------------------------------------
61      !!
62      !! Cool-skin parameterization, based on Fairall et al., 1996:
63      !!
64      !! Fairall, C. W., Bradley, E. F., Godfrey, J. S., Wick, G. A.,
65      !! Edson, J. B., and Young, G. S. ( 1996), Cool‐skin and warm‐layer
66      !! effects on sea surface temperature, J. Geophys. Res., 101( C1), 1295-1308,
67      !! doi:10.1029/95JC03190.
68      !!
69      !!------------------------------------------------------------------
70      !!
71      !!  **   INPUT:
72      !!     *pQsw*       surface net solar radiation into the ocean     [W/m^2] => >= 0 !
73      !!     *pQnsol*     surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 !
74      !!     *pustar*     friction velocity u*                           [m/s]
75      !!     *pSST*       bulk SST (taken at depth gdept_1d(1))          [K]
76      !!------------------------------------------------------------------
77      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQsw   ! net solar a.k.a shortwave radiation into the ocean (after albedo) [W/m^2]
78      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQnsol ! non-solar heat flux to the ocean [W/m^2]
79      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pustar  ! friction velocity, temperature and humidity (u*,t*,q*)
80      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pSST ! bulk SST [K]
81      !!---------------------------------------------------------------------
82      INTEGER  :: ji, jj, jc
83      REAL(wp) :: zQabs, zdelta, zfr
84      !!---------------------------------------------------------------------
85      DO jj = 1, jpj
86         DO ji = 1, jpi
87
88            zQabs = pQnsol(ji,jj) ! first guess of heat flux absorbed within the viscous sublayer of thicknes delta,
89            !                     !   => we DO not miss a lot assuming 0 solar flux absorbed in the tiny layer of thicknes zdelta...
90
91            zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pustar(ji,jj) )
92
93            DO jc = 1, 4 ! because implicit in terms of zdelta...
94               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
95               !              =>  (WARNING: 0.065 rather than 0.137 in Fairal et al. 1996)
96               zQabs = pQnsol(ji,jj) + zfr*pQsw(ji,jj)
97               zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pustar(ji,jj) )
98            END DO
99
100            dT_cs(ji,jj) = zQabs*zdelta/rk0_w   ! temperature increment, yes dT_cs can actually > 0, if Qabs > 0 (rare but possible!)
101
102         END DO
103      END DO
104
105   END SUBROUTINE CS_ECMWF
106
107
108
109   SUBROUTINE WL_ECMWF( pQsw, pQnsol, pustar, pSST,  pustk )
110      !!---------------------------------------------------------------------
111      !!
112      !!  Warm-Layer scheme according to Zeng & Beljaars, 2005 (GRL)
113      !!  " A prognostic scheme of sea surface skin temperature for modeling and data assimilation "
114      !!
115      !!  STIL NO PROGNOSTIC EQUATION FOR THE DEPTH OF THE WARM-LAYER!
116      !!
117      !!    As included in IFS Cy45r1   /  E.C.M.W.F.
118      !!     ------------------------------------------------------------------
119      !!
120      !!  **   INPUT:
121      !!     *pQsw*       surface net solar radiation into the ocean     [W/m^2] => >= 0 !
122      !!     *pQnsol*     surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 !
123      !!     *pustar*     friction velocity u*                           [m/s]
124      !!     *pSST*       bulk SST  (taken at depth gdept_1d(1))         [K]
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)  :: pustar   ! friction velocity [m/s]
129      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pSST     ! bulk SST at depth gdept_1d(1) [K]
130      !!
131      REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) :: pustk ! surface Stokes velocity [m/s]
132      !
133      INTEGER :: ji, jj, jc
134      !
135      REAL(wp) :: &
136         & zHwl,    &  !: thickness of the warm-layer [m]
137         & ztcorr,  &  !: correction of dT w.r.t measurement depth of bulk SST (first T-point)
138         & zalpha, & !: thermal expansion coefficient of sea-water [1/K]
139         & zdTwl_b, zdTwl_n, & ! temp. diff. between "almost surface (right below viscous layer) and bottom of WL
140         & zfr, zeta, &
141         & zusw, zusw2, &
142         & zLa, zfLa, &
143         & flg, zwf, zQabs, &
144         & ZA, ZB, zL1, zL2, &
145         &  zcst0, zcst1, zcst2, zcst3
146      !
147      LOGICAL :: l_pustk_known
148      !!---------------------------------------------------------------------
149
150      l_pustk_known = .FALSE.
151      IF ( PRESENT(pustk) ) l_pustk_known = .TRUE.
152
153      DO jj = 1, jpj
154         DO ji = 1, jpi
155
156            zHwl = Hz_wl(ji,jj) ! first guess for warm-layer depth (and unique..., less advanced than COARE3p6 !)
157            ! it is = rd0 (3m) in default Zeng & Beljaars case...
158
159            !! Previous value of dT / warm-layer, adapted to depth:
160            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$
161            ztcorr = flg + (1._wp - flg)*gdept_1d(1)/zHwl
162            zdTwl_b = MAX ( dT_wl(ji,jj) / ztcorr , 0._wp )
163            ! zdTwl is the difference between "almost surface (right below viscous layer) and bottom of WL (here zHwl)
164            ! pdT         "                          "                                    and depth of bulk SST (here gdept_1d(1))!
165            !! => but of course in general the bulk SST is taken shallower than zHwl !!! So correction less pronounced!
166            !! => so here since pdT is difference between surface and gdept_1d(1), need to increase fof zdTwl !
167
168            zalpha = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!)
169
170
171            ! *** zfr = Fraction of solar radiation absorbed in warm layer (-)
172            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
173
174            zQabs = zfr*pQsw(ji,jj) + pQnsol(ji,jj)       ! tot heat absorbed in warm layer
175
176            zusw  = MAX( pustar(ji,jj), 1.E-4_wp ) * sq_radrw    ! u* in the water
177            zusw2 = zusw*zusw
178
179            ! Langmuir:
180            IF ( l_pustk_known ) THEN
181               zLa = SQRT(zusw/MAX(pustk(ji,jj),1.E-6))
182            ELSE
183               zla = 0.3_wp
184            END IF
185            zfLa = MAX( zla**(-2._wp/3._wp) , 1._wp )   ! Eq.(6)
186
187            zwf = 0.5_wp + SIGN(0.5_wp, zQabs)  ! zQabs > 0. => 1.  / zQabs < 0. => 0.
188
189            zcst1 = vkarmn*grav*zalpha
190
191            ! 1/L when zQabs > 0 :
192            zL2 = zcst1*zQabs / (zRhoCp_w*zusw2*zusw)
193               
194            zcst2 = zcst1 / ( 5._wp*zHwl*zusw2 )  !OR: zcst2 = zcst1*rNuwl0 / ( 5._wp*zHwl*zusw2 ) ???
195
196            zcst0 = rdt * (rNuwl0 + 1._wp) / zHwl
197           
198            ZA = zcst0 * zQabs / ( rNuwl0 * zRhoCp_w )
199
200            zcst3 = -zcst0 * vkarmn * zusw * zfLa
201
202            !! Sorry about all these constants ( constant w.r.t zdTwl), it's for
203            !! the sake of optimizations... So all these operations are not done
204            !! over and over within the iteration loop...
205           
206            !! T R U L L Y   I M P L I C I T => needs itteration
207            !! => have to itterate just because the 1/(Monin-Obukhov length), zL1, uses zdTwl when zQabs < 0..
208            !!    (without this term otherwize the implicit analytical solution is straightforward...)
209            zdTwl_n = zdTwl_b
210            DO jc = 1, 10
211               
212               zdTwl_n = 0.5_wp * ( zdTwl_n + zdTwl_b ) ! semi implicit, for faster convergence
213               
214               ! 1/L when zdTwl > 0 .AND. zQabs < 0 :
215               zL1 =         SQRT( zdTwl_n * zcst2 ) ! / zusw !!! Or??? => vkarmn * SQRT( zdTwl_n*grav*zalpha/( 5._wp*zHwl ) ) / zusw
216               !zL1 = vkarmn*SQRT( zdTwl_n       *grav*zalpha        / ( 5._wp*zHwl ) ) / zusw   ! => vkarmn outside, not inside zcst1 (just for this particular line) ???
217               
218               ! Stability parameter (z/L):
219               zeta =  (1._wp - zwf) * zHwl*zL1   +   zwf * zHwl*zL2
220
221               ZB = zcst3 / PHI(zeta)
222
223               zdTwl_n = MAX ( zdTwl_b + ZA + ZB*zdTwl_n , 0._wp )  ! Eq.(6)
224
225            END DO
226           
227            !! Update:
228            dT_wl(ji,jj) = zdTwl_n * ztcorr
229           
230         END DO
231      END DO
232
233   END SUBROUTINE WL_ECMWF
234
235
236
237   FUNCTION delta_skin_layer( palpha, pQd, pustar_a )
238      !!---------------------------------------------------------------------
239      !! Computes the thickness (m) of the viscous skin layer.
240      !! Based on Fairall et al., 1996
241      !!
242      !! Fairall, C. W., Bradley, E. F., Godfrey, J. S., Wick, G. A.,
243      !! Edson, J. B., and Young, G. S. ( 1996), Cool‐skin and warm‐layer
244      !! effects on sea surface temperature, J. Geophys. Res., 101( C1), 1295-1308,
245      !! doi:10.1029/95JC03190.
246      !!
247      !! L. Brodeau, october 2019
248      !!---------------------------------------------------------------------
249      REAL(wp)                :: delta_skin_layer
250      REAL(wp), INTENT(in)    :: palpha   ! thermal expansion coefficient of sea-water (SST accurate enough!)
251      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
252      REAL(wp), INTENT(in)    :: pustar_a ! friction velocity in the air (u*) [m/s]
253      !!---------------------------------------------------------------------
254      REAL(wp) :: zusw, zusw2, zlamb, ztf, ztmp
255      !!---------------------------------------------------------------------
256      ztf = 0.5_wp + SIGN(0.5_wp, pQd)  ! Qabs < 0 => cooling of the viscous layer => ztf = 0 (regular case)
257      !                                 ! Qabs > 0 => warming of the viscous layer => ztf = 1 (ex: weak evaporation and strong positive sensible heat flux)
258      !
259      zusw  = MAX(pustar_a, 1.E-4_wp) * sq_radrw    ! u* in the water
260      zusw2 = zusw*zusw
261      !
262      zlamb = 6._wp*( 1._wp + MAX(palpha*rcst_cs/(zusw2*zusw2)*pQd, 0._wp)**0.75 )**(-1./3.) ! see Eq.(14) in Fairall et al., 1996
263      !  => 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
264      !
265      ztmp = rnu0_w/zusw
266      delta_skin_layer = (1._wp-ztf) *     zlamb*ztmp           &  ! regular case, Qd < 0, see Eq.(12) in Fairall et al., 1996
267         &               +   ztf     * MIN(6._wp*ztmp , 0.007_wp)  ! when Qd > 0
268   END FUNCTION delta_skin_layer
269
270
271   FUNCTION PHI( pzeta)
272      !!---------------------------------------------------------------------
273      !!
274      !! Takaya et al., 2010
275      !!  Eq.(5)
276      !! L. Brodeau, october 2019
277      !!---------------------------------------------------------------------
278      REAL(wp)                :: PHI
279      REAL(wp), INTENT(in)    :: pzeta    ! stability parameter
280      !!---------------------------------------------------------------------
281      REAL(wp) :: ztf, zzt2
282      !!---------------------------------------------------------------------
283      !
284      zzt2 = pzeta*pzeta
285      !
286      ztf = 0.5_wp + SIGN(0.5_wp, pzeta)  ! zeta > 0 => ztf = 1
287      !                                   ! zeta < 0 => ztf = 0
288      PHI =      ztf     * ( 1. + (5.*pzeta + 4.*zzt2)/(1. + 3.*pzeta + 0.25*zzt2) ) &   ! zeta > 0
289         &  + (1. - ztf) * 1./SQRT( 1. - 16.*(-ABS(pzeta)) )                             ! zeta < 0
290      !
291   END FUNCTION PHI
292
293   !!======================================================================
294END MODULE sbcblk_skin_ecmwf
Note: See TracBrowser for help on using the repository browser.