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_algo_andreas.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_algo_andreas.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: 16.7 KB
Line 
1!!! TO DO: consistent psi_m and psi_h needed!!! For now is those of NCAR !!!
2!!
3MODULE sbcblk_algo_andreas
4   !!======================================================================
5   !!                   ***  MODULE  sbcblk_algo_andreas  ***
6   !! Computes:
7   !!   * bulk transfer coefficients C_D, C_E and C_H
8   !!   * air temp. and spec. hum. adjusted from zt (2m) to zu (10m) if needed
9   !!   * the effective bulk wind speed at 10m Ubzu
10   !!  according to Andreas et al. (2015)
11   !!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
12   !!       Andreas, E.L., Mahrt, L. and Vickers, D. (2015),
13   !!       An improved bulk air–sea surface flux algorithm,
14   !!       including spray‐mediated transfer.
15   !!       Q.J.R. Meteorol. Soc., 141: 642-654. doi:10.1002/qj.2424
16   !!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
17   !!
18   !!   * bulk transfer coefficients C_D, C_E and C_H
19   !!   * air temp. and spec. hum. adjusted from zt (2m) to zu (10m) if needed
20   !!   * the effective bulk wind speed at z=zu: Ubzu
21   !!   => all these are used in bulk formulas in sbcblk.F90
22   !!
23   !!    Using the bulk formulation/param. of Large & Yeager 2008
24   !!
25   !!       Routine turb_andreas maintained and developed in AeroBulk
26   !!                     (https://github.com/brodeau/aerobulk/)
27   !!
28   !! ** Author: L. Brodeau, August 2020 / AeroBulk (https://github.com/brodeau/aerobulk)
29   !!----------------------------------------------------------------------
30   !! History :  4.x  !  2020-08  (L.Brodeau)   Original code
31   !!----------------------------------------------------------------------
32
33   !!----------------------------------------------------------------------
34   !!----------------------------------------------------------------------
35   USE dom_oce         ! ocean space and time domain
36   USE phycst          ! physical constants
37   USE sbc_phy         ! all thermodynamics functions, rho_air, q_sat, etc... !LB
38
39   IMPLICIT NONE
40   PRIVATE
41
42   !! Important (Brodeau fix):
43   REAL(wp), PARAMETER :: rRi_max = 0.15_wp   ! Bulk Ri above which the algorithm fucks up!
44   !                                          ! (increasing (>0) Ri means that surface layer increasingly stable and/or wind increasingly weak)
45   REAL(wp), PARAMETER :: rCs_min = 0.35E-3_wp ! minimum value to tolarate for CE and CH ! Must be larger than "Cx_min" !!!
46
47   PUBLIC :: TURB_ANDREAS, psi_m_andreas, psi_h_andreas
48
49   !! * Substitutions
50#  include "do_loop_substitute.h90"
51
52   !!----------------------------------------------------------------------
53CONTAINS
54
55   SUBROUTINE turb_andreas( zt, zu, sst, t_zt, ssq, q_zt, U_zu, &
56      &                     Cd, Ch, Ce, t_zu, q_zu, Ubzu,       &
57      &                     nb_iter, CdN, ChN, CeN               )
58      !!----------------------------------------------------------------------------------
59      !!                      ***  ROUTINE  turb_andreas  ***
60      !!
61      !! ** Purpose :   Computes turbulent transfert coefficients of surface
62      !!                fluxes according to Large & Yeager (2004) and Large & Yeager (2008)
63      !!                If relevant (zt /= zu), adjust temperature and humidity from height zt to zu
64      !!                Returns the effective bulk wind speed at zu to be used in the bulk formulas
65      !!
66      !! INPUT :
67      !! -------
68      !!    *  zt   : height for temperature and spec. hum. of air            [m]
69      !!    *  zu   : height for wind speed (usually 10m)                     [m]
70      !!    *  sst  : bulk SST                                                [K]
71      !!    *  t_zt : potential air temperature at zt                         [K]
72      !!    *  ssq  : specific humidity at saturation at SST                  [kg/kg]
73      !!    *  q_zt : specific humidity of air at zt                          [kg/kg]
74      !!    *  U_zu : scalar wind speed at zu                                 [m/s]
75      !!
76      !! OUTPUT :
77      !! --------
78      !!    *  Cd     : drag coefficient
79      !!    *  Ch     : sensible heat coefficient
80      !!    *  Ce     : evaporation coefficient
81      !!    *  t_zu   : pot. air temperature adjusted at wind height zu       [K]
82      !!    *  q_zu   : specific humidity of air        //                    [kg/kg]
83      !!    *  Ubzu   : bulk wind speed at zu                                 [m/s]
84      !!
85      !!
86      !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/)
87      !!----------------------------------------------------------------------------------
88      REAL(wp), INTENT(in   )                     ::   zt       ! height for t_zt and q_zt                    [m]
89      REAL(wp), INTENT(in   )                     ::   zu       ! height for U_zu                             [m]
90      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   sst      ! sea surface temperature                [Kelvin]
91      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   t_zt     ! potential air temperature              [Kelvin]
92      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   ssq      ! sea surface specific humidity           [kg/kg]
93      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_zt     ! specific air humidity at zt             [kg/kg]
94      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   U_zu     ! relative wind module at zu                [m/s]
95      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cd       ! transfer coefficient for momentum         (tau)
96      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ch       ! transfer coefficient for sensible heat (Q_sens)
97      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ce       ! transfert coefficient for evaporation   (Q_lat)
98      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   t_zu     ! pot. air temp. adjusted at zu               [K]
99      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   q_zu     ! spec. humidity adjusted at zu           [kg/kg]
100      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ubzu    ! bulk wind speed at zu                     [m/s]
101      !
102      INTEGER , INTENT(in   ), OPTIONAL                     :: nb_iter  ! number of iterations
103      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   CdN
104      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   ChN
105      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   CeN
106      !
107      INTEGER :: nbit, jit                    ! iterations...
108      LOGICAL :: l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U
109      !!
110      REAL(wp), DIMENSION(jpi,jpj) ::   u_star, t_star, q_star
111      REAL(wp), DIMENSION(jpi,jpj) ::   z0       ! roughness length (momentum) [m]
112      REAL(wp), DIMENSION(jpi,jpj) ::   UN10     ! Neutral wind speed at zu [m/s]
113      REAL(wp), DIMENSION(jpi,jpj) ::   zeta_u   ! stability parameter at height zu
114      REAL(wp), DIMENSION(jpi,jpj) ::   ztmp0, ztmp1, ztmp2
115      REAL(wp), DIMENSION(jpi,jpj) ::   RiB       ! square root of Cd
116      !!
117      !!----------------------------------------------------------------------------------
118      nbit = nb_iter0
119      IF( PRESENT(nb_iter) ) nbit = nb_iter
120
121      l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp ) ! testing "zu == zt" is risky with double precision
122
123      Ubzu = MAX( 0.25_wp , U_zu )   !  relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s
124
125      !! First guess:
126      UN10 = Ubzu
127      Cd   = 1.1E-3_wp
128      Ch   = 1.1E-3_wp
129      Ce   = 1.1E-3_wp
130      t_zu = t_zt
131      q_zu = q_zt
132
133      !! First guess of turbulent scales for scalars:
134      ztmp0  = SQRT(Cd)
135      t_star = Ch/ztmp0*(t_zu - sst) ! theta*
136      q_star = Ce/ztmp0*(q_zu - ssq) ! q*
137
138      ! Bulk Richardson number:
139      RiB(:,:) = Ri_bulk( zu, sst, t_zu, ssq, q_zu, Ubzu )
140
141
142      !! ITERATION BLOCK
143      DO jit = 1, nbit
144
145         WHERE ( RiB < rRi_max )
146            !! Normal condition case:
147            u_star = U_STAR_ANDREAS( UN10 )
148         ELSEWHERE
149            !! Extremely stable + weak wind !!!
150            !!  => for we force u* to be consistent with minimum value for CD:
151            !!  (otherwize algorithm becomes nonsense...)
152            u_star = SQRT(Cx_min) * Ubzu     ! Cd does not go below Cx_min !
153         ENDWHERE
154
155         !! Stability parameter :
156         zeta_u = zu*One_on_L( t_zu, q_zu, u_star, t_star, q_star )   ! zu * 1/L
157
158         !! Drag coefficient:
159         ztmp0 = u_star/Ubzu
160
161         Cd = MAX( ztmp0*ztmp0 , Cx_min )
162
163         !! Roughness length:
164         z0 = MIN( z0_from_Cd( zu, Cd,  ppsi=psi_m_andreas(zeta_u) ) , z0_sea_max )
165
166         !! z0t and z0q, based on LKB, just like into COARE 2.5:
167         ztmp0 = z0 * u_star / visc_air(t_zu) ! Re_r
168         ztmp1 = z0tq_LKB( 1, ztmp0, z0 )     ! z0t
169         ztmp2 = z0tq_LKB( 2, ztmp0, z0 )     ! z0q
170
171         !! Turbulent scales at zu :
172         ztmp0 = psi_h_andreas(zeta_u)  ! lolo: zeta_u for scalars???
173         t_star  = (t_zu - sst)*vkarmn/(LOG(zu) - LOG(ztmp1) - ztmp0)  ! theta* (ztmp1 == z0t in rhs term)
174         q_star  = (q_zu - ssq)*vkarmn/(LOG(zu) - LOG(ztmp2) - ztmp0)  !   q*   (ztmp2 == z0q in rhs term)
175
176         IF( (.NOT. l_zt_equal_zu).AND.( jit > 1 ) ) THEN
177            !! Re-updating temperature and humidity at zu if zt /= zu:
178            ztmp0 = zeta_u/zu*zt   ! zeta_t
179            ztmp0 = LOG(zt/zu) + psi_h_andreas(zeta_u) - psi_h_andreas(ztmp0)
180            t_zu = t_zt - t_star/vkarmn*ztmp0
181            q_zu = q_zt - q_star/vkarmn*ztmp0
182            RiB  = Ri_bulk( zu, sst, t_zu, ssq, q_zu, Ubzu ) !LOLO
183         ENDIF
184
185         !! Update neutral-stability wind at zu:
186         UN10 = MAX( 0.1_wp , UN10_from_ustar( zu, Ubzu, u_star, psi_m_andreas(zeta_u) ) ) ! UN10
187
188      END DO !DO jit = 1, nbit
189
190      ! Compute transfer coefficients at zu:
191      ztmp0 = u_star/Ubzu
192
193      Cd = MAX( ztmp0*ztmp0 , Cx_min )   ! the earlier use of Cx_min on u* should make use of Cx_min here unnecessary!
194
195      ztmp1 = t_zu - sst ;  ztmp1 = SIGN( MAX(ABS(ztmp1),1.E-6_wp), ztmp1 )  ! dt_zu
196      ztmp2 = q_zu - ssq ;  ztmp2 = SIGN( MAX(ABS(ztmp2),1.E-9_wp), ztmp2 )  ! dq_zu
197      Ch   = MAX( ztmp0*t_star/ztmp1 , rCs_min )
198      Ce   = MAX( ztmp0*q_star/ztmp2 , rCs_min )
199
200      !! Neutral-stability coefficients:
201      ztmp0 = 1._wp/LOG(zu/z0)
202      ztmp1 = z0 * u_star / visc_air(t_zu)  ! Re_r
203     
204      IF(PRESENT(CdN)) CdN = vkarmn2*ztmp0*ztmp0
205      IF(PRESENT(ChN)) ChN = vkarmn2*ztmp0/LOG(zu/z0tq_LKB( 1, ztmp1, z0 ))
206      IF(PRESENT(CeN)) CeN = vkarmn2*ztmp0/LOG(zu/z0tq_LKB( 2, ztmp1, z0 ))
207     
208   END SUBROUTINE turb_andreas
209
210
211   FUNCTION U_STAR_ANDREAS( pun10 )
212      !!----------------------------------------------------------------------------------
213      !! Estimate of the friction velocity as a function of the neutral-stability wind
214      !! speed at at 10m
215      !!
216      !! Origin: Eq.(2.2) of Andreas et al. (2015)
217      !!
218      !! ** Author: L. Brodeau, April 2020 / AeroBulk (https://github.com/brodeau/aerobulk/)
219      !!----------------------------------------------------------------------------------
220      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pun10          !: neutral-stability scalar wind speed at 10m (m/s)
221      REAL(wp), DIMENSION(jpi,jpj)             :: u_star_andreas !: friction velocity    [m/s]
222      !
223      INTEGER  ::     ji, jj ! dummy loop indices
224      REAL(wp) :: za, zt, zw ! local scalars
225      !!----------------------------------------------------------------------------------
226      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
227            zw  = pun10(ji,jj)
228            za = zw - 8.271_wp
229            zt = za + SQRT( 0.12_wp*za*za + 0.181_wp )
230            u_star_andreas(ji,jj) =   0.239_wp + 0.0433_wp * zt
231      END_2D
232   END FUNCTION U_STAR_ANDREAS
233
234
235   FUNCTION psi_m_andreas( pzeta )
236      !!----------------------------------------------------------------------------------
237      !!      Universal profile stability function for momentum
238      !!  TO DO !!!!!!!!!!!!!!!!!!!!!
239      !! LOLO: paper says Paulson 1970 when unstable and Grachev et al 2007 for STABLE
240      !!
241      !! pzeta : stability paramenter, z/L where z is altitude measurement
242      !!         and L is M-O length
243      !!
244      !! ** Author: L. Brodeau, April 2020 / AeroBulk (https://github.com/brodeau/aerobulk/)
245      !!----------------------------------------------------------------------------------
246      REAL(wp), DIMENSION(jpi,jpj) :: psi_m_andreas
247      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
248      !
249      REAL(wp), PARAMETER :: zam  = 5._wp      ! a_m (just below Eq.(9b)
250      REAL(wp), PARAMETER :: zbm = zam/6.5_wp  ! b_m (just below Eq.(9b)
251      !
252      REAL(wp), PARAMETER :: z1o3 = 1._wp/3._wp
253      REAL(wp), PARAMETER :: zsr3 = SQRT(3._wp)
254      !
255      INTEGER  ::   ji, jj    ! dummy loop indices
256      REAL(wp) :: zta, zx2, zx, zpsi_unst, zbbm, zpsi_stab,  zstab   ! local scalars
257      !!----------------------------------------------------------------------------------
258      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
259            !
260            zta = MIN( pzeta(ji,jj) , 15._wp ) !! Very stable conditions (L positif and big!)
261            !
262            !! *** Unstable: Paulson (1970): #LOLO: DOUBLE CHECK IT IS PAULSON!!!!!
263            zx2 = SQRT( ABS(1._wp - 16._wp*zta) )  ! (1 - 16z)^0.5
264            zx2 = MAX( zx2 , 1._wp )
265            zx  = SQRT(zx2)                          ! (1 - 16z)^0.25
266            zpsi_unst = 2._wp*LOG(ABS( (1._wp + zx )*0.5_wp ))   &
267               &            + LOG(ABS( (1._wp + zx2)*0.5_wp ))   &
268               &          - 2._wp*ATAN(zx) + rpi*0.5_wp
269            !
270            !! *** Stable: Grachev et al 2007 (SHEBA) [Eq.(12) Grachev et al 2007]:
271            zx   = ABS(1._wp + zta)**z1o3
272            zbbm = ABS( (1._wp - zbm)/zbm )**z1o3 ! B_m
273            !
274            zpsi_stab = -3.*zam/zbm*(zx - 1._wp) + zam*zbbm/(2.*zbm) * ( &
275               &        2.*LOG(ABS( (   zx     +   zbbm         )/(1._wp        +   zbbm   ) )) &
276               &         - LOG(ABS( (zx*zx - zx*zbbm + zbbm*zbbm)/(1._wp - zbbm + zbbm*zbbm) )) &
277               & + 2.*zsr3*( ATAN( (2.*zx - zbbm)/(zsr3*zbbm) ) - ATAN( (2._wp - zbbm)/(zsr3*zbbm) ) ) )
278            !
279            !
280            zstab = 0.5_wp + SIGN(0.5_wp, zta) ! zta > 0 => zstab = 1
281            !
282            psi_m_andreas(ji,jj) =       zstab  * zpsi_stab &  ! (zta > 0) Stable
283               &              + (1._wp - zstab) * zpsi_unst    ! (zta < 0) Unstable
284            !
285      END_2D
286   END FUNCTION psi_m_andreas
287
288
289   FUNCTION psi_h_andreas( pzeta )
290      !!----------------------------------------------------------------------------------
291      !! Universal profile stability function for temperature and humidity
292      !!
293      !!    TO DO
294      !!       !! LOLO: paper says Paulson 1970 when unstable and Grachev et al 2007 for STABLE
295      !!
296      !! pzeta : stability paramenter, z/L where z is altitude measurement
297      !!         and L is M-O length
298      !!
299      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
300      !!----------------------------------------------------------------------------------
301      REAL(wp), DIMENSION(jpi,jpj) :: psi_h_andreas
302      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
303      !
304      REAL(wp), PARAMETER ::  zah = 5._wp       ! a_h (just below Eq.(9b)
305      REAL(wp), PARAMETER ::  zbh = 5._wp       ! b_h (just below Eq.(9b)
306      REAL(wp), PARAMETER ::  zch = 3._wp       ! c_h (just below Eq.(9b)
307      REAL(wp), PARAMETER :: zbbh = SQRT(5._wp) ! B_h (just below Eq.(13)
308      !
309      INTEGER  ::   ji, jj     ! dummy loop indices
310      REAL(wp) :: zta, zz, zx2, zpsi_unst, zpsi_stab, zstab  ! local scalars
311      !!----------------------------------------------------------------------------------
312      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
313            !
314            zta = MIN( pzeta(ji,jj) , 15._wp ) !! Very stable conditions (L positif and large!)
315            !
316            !! *** Unstable: Paulson (1970): #LOLO: DOUBLE CHECK IT IS PAULSON!!!!!
317            zx2 = SQRT( ABS(1._wp - 16._wp*zta) )  ! (1 -16z)^0.5
318            zx2 = MAX( zx2 , 1._wp )
319            zpsi_unst = 2._wp*LOG( 0.5_wp*(1._wp + zx2) )
320            !
321            !! *** Stable: Grachev et al 2007 (SHEBA) [Eq.(13) Grachev et al 2007]:
322            zz = 2.*zta + zch
323            zpsi_stab = - 0.5*zbh*LOG(ABS(1._wp + zch*zta + zta*zta)) &
324               &        +  (-zah/zbbh + 0.5*zbh*zch/zbbh)  &
325               &          *( LOG(ABS((zz  - zbbh)/(zz  + zbbh))) &
326               &           - LOG(ABS((zch - zbbh)/(zch + zbbh)))    )
327            !
328            zstab = 0.5_wp + SIGN(0.5_wp, zta) ! zta > 0 => zstab = 1
329            !
330            psi_h_andreas(ji,jj) =            zstab  * zpsi_stab &  ! (zta > 0) Stable
331               &                   + (1._wp - zstab) * zpsi_unst    ! (zta < 0) Unstable
332            !
333      END_2D
334   END FUNCTION psi_h_andreas
335
336   !!======================================================================
337END MODULE sbcblk_algo_andreas
Note: See TracBrowser for help on using the repository browser.