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/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC – NEMO

source: NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbcblk_skin_ecmwf.F90 @ 13655

Last change on this file since 13655 was 13655, checked in by laurent, 3 years ago

Commit all my dev of 2020!

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