source: NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/SBC/sbcblk_skin_ecmwf.F90 @ 12154

Last change on this file since 12154 was 12154, checked in by cetlod, 11 months ago

commit

File size: 15.7 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 sbcblk_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
48   !! Cool-skin related parameters:
49   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: dT_cs !: dT due to cool-skin effect
50   !                                                      ! => temperature difference between air-sea interface (z=0)
51   !                                                      !    and right below viscous layer (z=delta)
52
53   !! Warm-layer related parameters:
54   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: dT_wl !: dT due to warm-layer effect
55   !                                                      ! => difference between "almost surface (right below
56   !                                                      !    viscous layer, z=delta) and depth of bulk SST (z=gdept_1d(1))
57   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: Hz_wl !: depth (aka thickness) of warm-layer [m]
58   !
59   REAL(wp), PARAMETER, PUBLIC :: rd0  = 3.    !: Depth scale [m] of warm layer, "d" in Eq.11 (Zeng & Beljaars 2005)
60   REAL(wp), PARAMETER         :: zRhoCp_w = rho0_w*rCp0_w
61   !
62   REAL(wp), PARAMETER         :: rNuwl0 = 0.5  !: Nu (exponent of temperature profile) Eq.11
63   !                                            !: (Zeng & Beljaars 2005) !: set to 0.5 instead of
64   !                                            !: 0.3 to respect a warming of +3 K in calm
65   !                                            !: condition for the insolation peak of +1000W/m^2
66   !!----------------------------------------------------------------------
67CONTAINS
68
69
70   SUBROUTINE CS_ECMWF( pQsw, pQnsol, pustar, pSST )
71      !!---------------------------------------------------------------------
72      !!
73      !! Cool-skin parameterization, based on Fairall et al., 1996:
74      !!
75      !!  - Zeng X., and A. Beljaars, 2005: A prognostic scheme of sea surface
76      !!    skin temperature for modeling and data assimilation. Geophysical
77      !!    Research Letters, 32 (14) , pp. 1-4.
78      !!
79      !!------------------------------------------------------------------
80      !!
81      !!  **   INPUT:
82      !!     *pQsw*       surface net solar radiation into the ocean     [W/m^2] => >= 0 !
83      !!     *pQnsol*     surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 !
84      !!     *pustar*     friction velocity u*                           [m/s]
85      !!     *pSST*       bulk SST (taken at depth gdept_1d(1))          [K]
86      !!------------------------------------------------------------------
87      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQsw   ! net solar a.k.a shortwave radiation into the ocean (after albedo) [W/m^2]
88      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQnsol ! non-solar heat flux to the ocean [W/m^2]
89      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pustar  ! friction velocity, temperature and humidity (u*,t*,q*)
90      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pSST ! bulk SST [K]
91      !!---------------------------------------------------------------------
92      INTEGER  :: ji, jj, jc
93      REAL(wp) :: zQabs, zdlt, zfr, zalfa, zus
94      !!---------------------------------------------------------------------
95      DO jj = 1, jpj
96         DO ji = 1, jpi
97
98            zQabs = pQnsol(ji,jj) ! first guess of heat flux absorbed within the viscous sublayer of thicknes delta,
99            !                     !   => we DO not miss a lot assuming 0 solar flux absorbed in the tiny layer of thicknes zdlt...
100
101            zalfa = alpha_sw(pSST(ji,jj)) ! (crude) thermal expansion coefficient of sea-water [1/K]
102            zus   = pustar(ji,jj)
103
104            zdlt = delta_skin_layer( zalfa, zQabs, zus )
105
106            DO jc = 1, 4 ! because implicit in terms of zdlt...
107               zfr = MAX( 0.065_wp + 11._wp*zdlt  &
108                  &       - 6.6E-5_wp/zdlt*(1._wp - EXP(-zdlt/8.E-4_wp)) &
109                  &      , 0.01_wp ) ! Solar absorption, Eq.(5) Zeng & Beljaars, 2005
110               !                     !  =>  (WARNING: 0.065 rather than 0.137 in Fairal et al. 1996)
111               zQabs = pQnsol(ji,jj) + zfr*pQsw(ji,jj)
112               zdlt = delta_skin_layer( zalfa, zQabs, zus )
113            END DO
114
115            dT_cs(ji,jj) = zQabs*zdlt/rk0_w   ! temperature increment, yes dT_cs can actually > 0, if Qabs > 0 (rare but possible!)
116
117         END DO
118      END DO
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 jj = 1, jpj
176         DO ji = 1, jpi
177
178            zHwl = Hz_wl(ji,jj) ! first guess for warm-layer depth (and unique..., less advanced than COARE3p6 !)
179            ! it is = rd0 (3m) in default Zeng & Beljaars case...
180
181            !! Previous value of dT / warm-layer, adapted to depth:
182            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$
183            ztcorr = flg + (1._wp - flg)*gdept_1d(1)/zHwl
184            zdTwl_b = MAX ( dT_wl(ji,jj) / ztcorr , 0._wp )
185            ! zdTwl is the difference between "almost surface (right below viscous layer) and bottom of WL (here zHwl)
186            ! pdT         "                          "                                    and depth of bulk SST (here gdept_1d(1))!
187            !! => but of course in general the bulk SST is taken shallower than zHwl !!! So correction less pronounced!
188            !! => so here since pdT is difference between surface and gdept_1d(1), need to increase fof zdTwl !
189
190            zalfa = alpha_sw( pSST(ji,jj) ) ! (crude) thermal expansion coefficient of sea-water [1/K] (SST accurate enough!)
191
192            ! *** zfr = Fraction of solar radiation absorbed in warm layer (-)
193            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
194
195            zQabs = zfr*pQsw(ji,jj) + pQnsol(ji,jj)       ! tot heat absorbed in warm layer
196
197            zusw  = MAX( pustar(ji,jj), 1.E-4_wp ) * sq_radrw    ! u* in the water
198            zusw2 = zusw*zusw
199
200            ! Langmuir:
201            IF( l_pustk_known ) THEN
202               zLa = SQRT(zusw/MAX(pustk(ji,jj),1.E-6))
203            ELSE
204               zla = 0.3_wp
205            ENDIF
206            zfLa = MAX( zla**(-2._wp/3._wp) , 1._wp )   ! Eq.(6)
207
208            zwf = 0.5_wp + SIGN(0.5_wp, zQabs)  ! zQabs > 0. => 1.  / zQabs < 0. => 0.
209
210            zcst1 = vkarmn*grav*zalfa
211
212            ! 1/L when zQabs > 0 :
213            zL2 = zcst1*zQabs / (zRhoCp_w*zusw2*zusw)
214
215            zcst2 = zcst1 / ( 5._wp*zHwl*zusw2 )  !OR: zcst2 = zcst1*rNuwl0 / ( 5._wp*zHwl*zusw2 ) ???
216
217            zcst0 = rdt * (rNuwl0 + 1._wp) / zHwl
218
219            ZA = zcst0 * zQabs / ( rNuwl0 * zRhoCp_w )
220
221            zcst3 = -zcst0 * vkarmn * zusw * zfLa
222
223            !! Sorry about all these constants ( constant w.r.t zdTwl), it's for
224            !! the sake of optimizations... So all these operations are not done
225            !! over and over within the iteration loop...
226
227            !! T R U L L Y   I M P L I C I T => needs itteration
228            !! => have to itterate just because the 1/(Monin-Obukhov length), zL1, uses zdTwl when zQabs < 0..
229            !!    (without this term otherwize the implicit analytical solution is straightforward...)
230            zdTwl_n = zdTwl_b
231            DO jc = 1, 10
232
233               zdTwl_n = 0.5_wp * ( zdTwl_n + zdTwl_b ) ! semi implicit, for faster convergence
234
235               ! 1/L when zdTwl > 0 .AND. zQabs < 0 :
236               zL1 =         SQRT( zdTwl_n * zcst2 ) ! / zusw !!! Or??? => vkarmn * SQRT( zdTwl_n*grav*zalfa/( 5._wp*zHwl ) ) / zusw
237
238               ! Stability parameter (z/L):
239               zeta =  (1._wp - zwf) * zHwl*zL1   +   zwf * zHwl*zL2
240
241               ZB = zcst3 / PHI(zeta)
242
243               zdTwl_n = MAX ( zdTwl_b + ZA + ZB*zdTwl_n , 0._wp )  ! Eq.(6)
244
245            END DO
246
247            !! Update:
248            dT_wl(ji,jj) = zdTwl_n * ztcorr
249
250         END DO
251      END DO
252
253   END SUBROUTINE WL_ECMWF
254
255
256   FUNCTION delta_skin_layer( palpha, pQd, pustar_a )
257      !!---------------------------------------------------------------------
258      !! Computes the thickness (m) of the viscous skin layer.
259      !! Based on Fairall et al., 1996
260      !!
261      !! Fairall, C. W., Bradley, E. F., Godfrey, J. S., Wick, G. A.,
262      !! Edson, J. B., and Young, G. S. ( 1996), Cool‐skin and warm‐layer
263      !! effects on sea surface temperature, J. Geophys. Res., 101( C1), 1295-1308,
264      !! doi:10.1029/95JC03190.
265      !!
266      !! L. Brodeau, october 2019
267      !!---------------------------------------------------------------------
268      REAL(wp)                :: delta_skin_layer
269      REAL(wp), INTENT(in)    :: palpha   ! thermal expansion coefficient of sea-water (SST accurate enough!)
270      REAL(wp), INTENT(in)    :: pQd ! < 0 !!! part of the net heat flux actually absorbed in the WL [W/m^2]
271      !                              !  => term "Q + Rs*fs" in eq.6 of Fairall et al. 1996
272      REAL(wp), INTENT(in)    :: pustar_a ! friction velocity in the air (u*) [m/s]
273      !!---------------------------------------------------------------------
274      REAL(wp) :: zusw, zusw2, zlamb, ztf, ztmp
275      !!---------------------------------------------------------------------
276      ztf = 0.5_wp + SIGN(0.5_wp, pQd)  ! Qabs < 0 => cooling of the viscous layer => ztf = 0 (regular case)
277      !                                 ! Qabs > 0 => warming of the viscous layer => ztf = 1
278      !                                 !    (ex: weak evaporation and strong positive sensible heat flux)
279      zusw  = MAX(pustar_a, 1.E-4_wp) * sq_radrw    ! u* in the water
280      zusw2 = zusw*zusw
281      !
282      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
283      !  => 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
284      !
285      ztmp = rnu0_w/zusw
286      delta_skin_layer = (1._wp-ztf) *     zlamb*ztmp           &  ! regular case, Qd < 0, see Eq.(12) in Fairall et al., 1996
287         &               +   ztf     * MIN(6._wp*ztmp , 0.007_wp)  ! when Qd > 0
288   END FUNCTION delta_skin_layer
289
290
291   FUNCTION PHI( pzeta)
292      !!---------------------------------------------------------------------
293      !!
294      !! Takaya et al., 2010
295      !!  Eq.(5)
296      !! L. Brodeau, october 2019
297      !!---------------------------------------------------------------------
298      REAL(wp)                :: PHI
299      REAL(wp), INTENT(in)    :: pzeta    ! stability parameter
300      !!---------------------------------------------------------------------
301      REAL(wp) :: ztf, zzt2
302      !!---------------------------------------------------------------------
303      zzt2 = pzeta*pzeta
304      ztf = 0.5_wp + SIGN(0.5_wp, pzeta)  ! zeta > 0 => ztf = 1
305      !                                   ! zeta < 0 => ztf = 0
306      PHI =      ztf     * ( 1. + (5.*pzeta + 4.*zzt2)/(1. + 3.*pzeta + 0.25*zzt2) ) &   ! zeta > 0
307         &  + (1. - ztf) * 1./SQRT( 1. - 16.*(-ABS(pzeta)) )                             ! zeta < 0
308   END FUNCTION PHI
309
310   !!======================================================================
311END MODULE sbcblk_skin_ecmwf
Note: See TracBrowser for help on using the repository browser.