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/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/SBC – NEMO

source: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/SBC/sbcblk_algo_andreas.F90 @ 15540

Last change on this file since 15540 was 15540, checked in by sparonuz, 3 years ago

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