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_ncar.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_ncar.F90 @ 13655

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

Commit all my dev of 2020!

  • Property svn:keywords set to Id
File size: 17.4 KB
Line 
1MODULE sbcblk_algo_ncar
2   !!======================================================================
3   !!                   ***  MODULE  sbcblk_algo_ncar  ***
4   !! Computes:
5   !!   * bulk transfer coefficients C_D, C_E and C_H
6   !!   * air temp. and spec. hum. adjusted from zt (2m) to zu (10m) if needed
7   !!   * the effective bulk wind speed at 10m Ubzu
8   !!   => all these are used in bulk formulas in sbcblk.F90
9   !!
10   !!    Using the bulk formulation/param. of Large & Yeager 2008
11   !!
12   !!       Routine turb_ncar maintained and developed in AeroBulk
13   !!                     (https://github.com/brodeau/aerobulk/)
14   !!
15   !!                         L. Brodeau, 2015
16   !!=====================================================================
17   !! History :  3.6  !  2016-02  (L.Brodeau) successor of old turb_ncar of former sbcblk_core.F90
18   !!----------------------------------------------------------------------
19
20   !!----------------------------------------------------------------------
21   !!   turb_ncar  : computes the bulk turbulent transfer coefficients
22   !!                   adjusts t_air and q_air from zt to zu m
23   !!                   returns the effective bulk wind speed at 10m
24   !!----------------------------------------------------------------------
25   USE dom_oce         ! ocean space and time domain
26   USE sbc_oce, ONLY: ln_cdgw
27   USE sbcwave, ONLY: cdn_wave ! wave module
28   USE phycst          ! physical constants
29   USE sbc_phy         ! all thermodynamics functions, rho_air, q_sat, etc... !LB
30
31   IMPLICIT NONE
32   PRIVATE
33
34   PUBLIC :: TURB_NCAR   ! called by sbcblk.F90
35
36   !! * Substitutions
37#  include "do_loop_substitute.h90"
38
39   !!----------------------------------------------------------------------
40CONTAINS
41
42   SUBROUTINE turb_ncar(    zt, zu, sst, t_zt, ssq, q_zt, U_zu, &
43      &                     Cd, Ch, Ce, t_zu, q_zu, Ubzu,       &
44      &                     nb_iter, CdN, ChN, CeN               )
45      !!----------------------------------------------------------------------------------
46      !!                      ***  ROUTINE  turb_ncar  ***
47      !!
48      !! ** Purpose :   Computes turbulent transfert coefficients of surface
49      !!                fluxes according to Large & Yeager (2004) and Large & Yeager (2008)
50      !!                If relevant (zt /= zu), adjust temperature and humidity from height zt to zu
51      !!                Returns the effective bulk wind speed at zu to be used in the bulk formulas
52      !!
53      !! INPUT :
54      !! -------
55      !!    *  zt   : height for temperature and spec. hum. of air            [m]
56      !!    *  zu   : height for wind speed (usually 10m)                     [m]
57      !!    *  sst  : bulk SST                                                [K]
58      !!    *  t_zt : potential air temperature at zt                         [K]
59      !!    *  ssq  : specific humidity at saturation at SST                  [kg/kg]
60      !!    *  q_zt : specific humidity of air at zt                          [kg/kg]
61      !!    *  U_zu : scalar wind speed at zu                                 [m/s]
62      !!
63      !!
64      !! OUTPUT :
65      !! --------
66      !!    *  Cd     : drag coefficient
67      !!    *  Ch     : sensible heat coefficient
68      !!    *  Ce     : evaporation coefficient
69      !!    *  t_zu   : pot. air temperature adjusted at wind height zu       [K]
70      !!    *  q_zu   : specific humidity of air        //                    [kg/kg]
71      !!    *  Ubzu   : bulk wind speed at zu                                 [m/s]
72      !!
73      !! OPTIONAL OUTPUT:
74      !! ----------------
75      !!    * CdN      : neutral-stability drag coefficient
76      !!    * ChN      : neutral-stability sensible heat coefficient
77      !!    * CeN      : neutral-stability evaporation coefficient
78      !!
79      !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/)
80      !!----------------------------------------------------------------------------------
81      REAL(wp), INTENT(in   )                     ::   zt       ! height for t_zt and q_zt                    [m]
82      REAL(wp), INTENT(in   )                     ::   zu       ! height for U_zu                             [m]
83      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   sst      ! sea surface temperature                [Kelvin]
84      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   t_zt     ! potential air temperature              [Kelvin]
85      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   ssq      ! sea surface specific humidity           [kg/kg]
86      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_zt     ! specific air humidity at zt             [kg/kg]
87      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   U_zu     ! relative wind module at zu                [m/s]
88      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cd       ! transfer coefficient for momentum         (tau)
89      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ch       ! transfer coefficient for sensible heat (Q_sens)
90      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ce       ! transfert coefficient for evaporation   (Q_lat)
91      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   t_zu     ! pot. air temp. adjusted at zu               [K]
92      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   q_zu     ! spec. humidity adjusted at zu           [kg/kg]
93      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ubzu    ! bulk wind speed at zu                     [m/s]
94      !
95      INTEGER , INTENT(in   ), OPTIONAL                     :: nb_iter  ! number of iterations
96      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   CdN
97      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   ChN
98      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   CeN
99      !
100      INTEGER :: nbit, jit                    ! iterations...
101      LOGICAL :: l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U
102      !
103      REAL(wp), DIMENSION(jpi,jpj) ::   zCdN, zCeN, zChN        ! 10m neutral latent/sensible coefficient
104      REAL(wp), DIMENSION(jpi,jpj) ::   zsqrt_Cd, zsqrt_CdN   ! root square of Cd and Cd_neutral
105      REAL(wp), DIMENSION(jpi,jpj) ::   zeta_u        ! stability parameter at height zu
106      REAL(wp), DIMENSION(jpi,jpj) ::   ztmp0, ztmp1, ztmp2
107      !!----------------------------------------------------------------------------------
108      nbit = nb_iter0
109      IF( PRESENT(nb_iter) ) nbit = nb_iter
110
111      l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp ) ! testing "zu == zt" is risky with double precision
112
113      Ubzu = MAX( 0.5_wp , U_zu )   !  relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s
114
115      !! First guess of stability:
116      ztmp0 = virt_temp(t_zt, q_zt) - virt_temp(sst, ssq) ! air-sea difference of virtual pot. temp. at zt
117      ztmp1 = 0.5_wp + SIGN(0.5_wp,ztmp0)                 ! ztmp1 = 1 if dTv > 0  => STABLE, 0 if unstable
118
119      !! Neutral coefficients at 10m:
120      IF( ln_cdgw ) THEN      ! wave drag case
121         cdn_wave(:,:) = cdn_wave(:,:) + rsmall * ( 1._wp - tmask(:,:,1) )
122         zCdN   (:,:) = cdn_wave(:,:)
123      ELSE
124      zCdN = cd_n10_ncar( Ubzu )
125      ENDIF
126
127      zsqrt_CdN = SQRT( zCdN )
128
129      !! Initializing transf. coeff. with their first guess neutral equivalents :
130      Cd = zCdN
131      Ce = ce_n10_ncar( zsqrt_CdN )
132      Ch = ch_n10_ncar( zsqrt_CdN , ztmp1 )   ! ztmp1 is stability (1/0)
133      zsqrt_Cd = zsqrt_CdN
134
135      IF( ln_cdgw ) THEN
136         zCeN = Ce
137         zChN = Ch
138      ENDIF
139
140      !! Initializing values at z_u with z_t values:
141      t_zu = MAX( t_zt ,  180._wp )   ! who knows what's given on masked-continental regions...
142      q_zu = MAX( q_zt , 1.e-6_wp )   !               "
143
144
145      !! ITERATION BLOCK
146      DO jit = 1, nbit
147         !
148         ztmp1 = t_zu - sst   ! Updating air/sea differences
149         ztmp2 = q_zu - ssq
150
151         ! Updating turbulent scales :   (L&Y 2004 Eq. (7))
152         ztmp0 = zsqrt_Cd*Ubzu       ! u*
153         ztmp1 = Ch/zsqrt_Cd*ztmp1    ! theta*
154         ztmp2 = Ce/zsqrt_Cd*ztmp2    ! q*
155
156         ! Estimate the inverse of Obukov length (1/L) at height zu:
157         ztmp0 = One_on_L( t_zu, q_zu, ztmp0, ztmp1, ztmp2 )
158
159         !! Stability parameters :
160         zeta_u   = zu*ztmp0
161         zeta_u   = sign( min(abs(zeta_u),10._wp), zeta_u )
162
163         !! Shifting temperature and humidity at zu (L&Y 2004 Eq. (9b-9c))
164         IF( .NOT. l_zt_equal_zu ) THEN
165            ztmp0 = zt*ztmp0 ! zeta_t !
166            ztmp0 = SIGN( MIN(ABS(ztmp0),10._wp), ztmp0 )  ! Temporaty array ztmp0 == zeta_t !!!
167            ztmp0 = LOG(zt/zu) + psi_h_ncar(zeta_u) - psi_h_ncar(ztmp0)                   ! ztmp0 just used as temp array again!
168            t_zu = t_zt - ztmp1/vkarmn*ztmp0    ! ztmp1 is still theta*  L&Y 2004 Eq. (9b)
169            !!
170            q_zu = q_zt - ztmp2/vkarmn*ztmp0    ! ztmp2 is still q*      L&Y 2004 Eq. (9c)
171            q_zu = MAX(0._wp, q_zu)
172         END IF
173
174         ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 Eq. 9a)...
175         !   In very rare low-wind conditions, the old way of estimating the
176         !   neutral wind speed at 10m leads to a negative value that causes the code
177         !   to crash. To prevent this a threshold of 0.25m/s is imposed.
178         ztmp2 = psi_m_ncar(zeta_u)
179         IF( ln_cdgw ) THEN      ! surface wave case
180            zsqrt_Cd = vkarmn / ( vkarmn / zsqrt_CdN - ztmp2 )
181            Cd   = zsqrt_Cd * zsqrt_Cd
182            ztmp0 = (LOG(zu/10._wp) - psi_h_ncar(zeta_u)) / vkarmn / zsqrt_CdN
183            ztmp2 = zsqrt_Cd / zsqrt_CdN
184            ztmp1 = 1._wp + zChN * ztmp0
185            Ch    = zChN * ztmp2 / ztmp1  ! L&Y 2004 eq. (10b)
186            ztmp1 = 1._wp + zCeN * ztmp0
187            Ce    = zCeN * ztmp2 / ztmp1  ! L&Y 2004 eq. (10c)
188
189         ELSE
190         ztmp0 = MAX( 0.25_wp , UN10_from_CD(zu, Ubzu, Cd, ppsi=ztmp2) ) ! U_n10 (ztmp2 == psi_m_ncar(zeta_u))
191
192         zCdN = cd_n10_ncar(ztmp0)
193         zsqrt_CdN = sqrt(zCdN)
194
195         !! Update of transfer coefficients:
196
197         !! C_D
198         ztmp1  = 1._wp + zsqrt_CdN/vkarmn*(LOG(zu/10._wp) - ztmp2)   ! L&Y 2004 Eq. (10a) (ztmp2 == psi_m(zeta_u))
199         Cd     = MAX( zCdN / ( ztmp1*ztmp1 ), Cx_min )
200
201         !! C_H and C_E
202         zsqrt_Cd = SQRT( Cd )
203         ztmp0 = ( LOG(zu/10._wp) - psi_h_ncar(zeta_u) ) / vkarmn / zsqrt_CdN
204         ztmp2 = zsqrt_Cd / zsqrt_CdN
205
206         ztmp1 = 0.5_wp + SIGN(0.5_wp,zeta_u)                                ! update stability
207         zChN  = 1.e-3_wp * zsqrt_CdN*(18._wp*ztmp1 + 32.7_wp*(1._wp - ztmp1))  ! L&Y 2004 eq. (6c-6d)
208         zCeN  = 1.e-3_wp * (34.6_wp * zsqrt_CdN)                             ! L&Y 2004 eq. (6b)
209
210         Ch    = MAX( zChN*ztmp2 / ( 1._wp + zChN*ztmp0 ) , Cx_min ) ! L&Y 2004 eq. (10b)
211         Ce    = MAX( zCeN*ztmp2 / ( 1._wp + zCeN*ztmp0 ) , Cx_min ) ! L&Y 2004 eq. (10c)
212
213         ENDIF
214         
215      END DO !DO jit = 1, nbit
216     
217      IF(PRESENT(CdN)) CdN(:,:) = zCdN(:,:)
218      IF(PRESENT(CeN)) CeN(:,:) = zCeN(:,:)
219      IF(PRESENT(ChN)) ChN(:,:) = zChN(:,:)
220
221   END SUBROUTINE turb_ncar
222
223
224   FUNCTION cd_n10_ncar( pw10 )
225      !!----------------------------------------------------------------------------------
226      !! Estimate of the neutral drag coefficient at 10m as a function
227      !! of neutral wind  speed at 10m
228      !!
229      !! Origin: Large & Yeager 2008, Eq. (11)
230      !!
231      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
232      !!----------------------------------------------------------------------------------
233      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pw10           ! scalar wind speed at 10m (m/s)
234      REAL(wp), DIMENSION(jpi,jpj)             :: cd_n10_ncar
235      !
236      INTEGER  ::     ji, jj     ! dummy loop indices
237      REAL(wp) :: zgt33, zw, zw6 ! local scalars
238      !!----------------------------------------------------------------------------------
239      !
240      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
241            !
242            zw  = pw10(ji,jj)
243            zw6 = zw*zw*zw
244            zw6 = zw6*zw6
245            !
246            ! When wind speed > 33 m/s => Cyclone conditions => special treatment
247            zgt33 = 0.5_wp + SIGN( 0.5_wp, (zw - 33._wp) )   ! If pw10 < 33. => 0, else => 1
248            !
249            cd_n10_ncar(ji,jj) = 1.e-3_wp * ( &
250               &       (1._wp - zgt33)*( 2.7_wp/zw + 0.142_wp + zw/13.09_wp - 3.14807E-10_wp*zw6) & ! wind <  33 m/s
251               &      +    zgt33   *      2.34_wp )                                                 ! wind >= 33 m/s
252            !
253            cd_n10_ncar(ji,jj) = MAX( cd_n10_ncar(ji,jj), Cx_min )
254            !
255      END_2D
256      !
257   END FUNCTION cd_n10_ncar
258
259
260   FUNCTION ch_n10_ncar( psqrtcdn10 , pstab )
261      !!----------------------------------------------------------------------------------
262      !! Estimate of the neutral heat transfer coefficient at 10m      !!
263      !! Origin: Large & Yeager 2008, Eq. (9) and (12)
264
265      !!----------------------------------------------------------------------------------
266      REAL(wp), DIMENSION(jpi,jpj)             :: ch_n10_ncar
267      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psqrtcdn10 ! sqrt( CdN10 )
268      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pstab      ! stable ABL => 1 / unstable ABL => 0
269      !!----------------------------------------------------------------------------------
270      IF( ANY(pstab < -0.00001) .OR. ANY(pstab >  1.00001) ) THEN
271         PRINT *, 'ERROR: ch_n10_ncar@mod_blk_ncar.f90: pstab ='
272         PRINT *, pstab
273         STOP
274      END IF
275      !
276      ch_n10_ncar = MAX( 1.e-3_wp * psqrtcdn10*( 18._wp*pstab + 32.7_wp*(1._wp - pstab) )  , Cx_min )   ! Eq. (9) & (12) Large & Yeager, 2008
277      !
278   END FUNCTION ch_n10_ncar
279
280   FUNCTION ce_n10_ncar( psqrtcdn10 )
281      !!----------------------------------------------------------------------------------
282      !! Estimate of the neutral heat transfer coefficient at 10m      !!
283      !! Origin: Large & Yeager 2008, Eq. (9) and (13)
284      !!----------------------------------------------------------------------------------
285      REAL(wp), DIMENSION(jpi,jpj)             :: ce_n10_ncar
286      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psqrtcdn10 ! sqrt( CdN10 )
287      !!----------------------------------------------------------------------------------
288      ce_n10_ncar = MAX( 1.e-3_wp * ( 34.6_wp * psqrtcdn10 ) , Cx_min )
289      !
290   END FUNCTION ce_n10_ncar
291
292
293   FUNCTION psi_m_ncar( pzeta )
294      !!----------------------------------------------------------------------------------
295      !! Universal profile stability function for momentum
296      !!    !! Psis, L&Y 2004, Eq. (8c), (8d), (8e)
297      !!
298      !! pzeta : stability paramenter, z/L where z is altitude measurement
299      !!         and L is M-O length
300      !!
301      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
302      !!----------------------------------------------------------------------------------
303      REAL(wp), DIMENSION(jpi,jpj) :: psi_m_ncar
304      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
305      !
306      INTEGER  ::   ji, jj    ! dummy loop indices
307      REAL(wp) :: zta, zx2, zx, zpsi_unst, zpsi_stab,  zstab   ! local scalars
308      !!----------------------------------------------------------------------------------
309      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
310            zta = pzeta(ji,jj)
311            !
312            zx2 = SQRT( ABS(1._wp - 16._wp*zta) )  ! (1 - 16z)^0.5
313            zx2 = MAX( zx2 , 1._wp )
314            zx  = SQRT(zx2)                          ! (1 - 16z)^0.25
315            zpsi_unst = 2._wp*LOG( (1._wp + zx )*0.5_wp )   &
316               &            + LOG( (1._wp + zx2)*0.5_wp )   &
317               &          - 2._wp*ATAN(zx) + rpi*0.5_wp
318            !
319            zpsi_stab = -5._wp*zta
320            !
321            zstab = 0.5_wp + SIGN(0.5_wp, zta) ! zta > 0 => zstab = 1
322            !
323            psi_m_ncar(ji,jj) =          zstab  * zpsi_stab &  ! (zta > 0) Stable
324               &              + (1._wp - zstab) * zpsi_unst    ! (zta < 0) Unstable
325            !
326            !
327      END_2D
328   END FUNCTION psi_m_ncar
329
330
331   FUNCTION psi_h_ncar( pzeta )
332      !!----------------------------------------------------------------------------------
333      !! Universal profile stability function for temperature and humidity
334      !!    !! Psis, L&Y 2004, Eq. (8c), (8d), (8e)
335      !!
336      !! pzeta : stability paramenter, z/L where z is altitude measurement
337      !!         and L is M-O length
338      !!
339      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
340      !!----------------------------------------------------------------------------------
341      REAL(wp), DIMENSION(jpi,jpj) :: psi_h_ncar
342      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
343      !
344      INTEGER  ::   ji, jj     ! dummy loop indices
345      REAL(wp) :: zta, zx2, zpsi_unst, zpsi_stab, zstab  ! local scalars
346      !!----------------------------------------------------------------------------------
347      !
348      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
349            !
350            zta = pzeta(ji,jj)
351            !
352            zx2 = SQRT( ABS(1._wp - 16._wp*zta) )  ! (1 -16z)^0.5
353            zx2 = MAX( zx2 , 1._wp )
354            zpsi_unst = 2._wp*LOG( 0.5_wp*(1._wp + zx2) )
355            !
356            zpsi_stab = -5._wp*zta
357            !
358            zstab = 0.5_wp + SIGN(0.5_wp, zta) ! zta > 0 => zstab = 1
359            !
360            psi_h_ncar(ji,jj) =          zstab  * zpsi_stab &  ! (zta > 0) Stable
361               &              + (1._wp - zstab) * zpsi_unst    ! (zta < 0) Unstable
362            !
363      END_2D
364   END FUNCTION psi_h_ncar
365
366   !!======================================================================
367END MODULE sbcblk_algo_ncar
Note: See TracBrowser for help on using the repository browser.