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/r12581_ticket2418/src/OCE/SBC – NEMO

source: NEMO/branches/2020/r12581_ticket2418/src/OCE/SBC/sbcblk_algo_ncar.F90 @ 12623

Last change on this file since 12623 was 12623, checked in by smasson, 4 years ago

r12581_ticket2418: merge with trunk@12622, see #2418
  • Property svn:keywords set to Id
File size: 15.7 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 U_blk
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 oce             ! ocean dynamics and tracers
26   USE dom_oce         ! ocean space and time domain
27   USE phycst          ! physical constants
28   USE sbc_oce         ! Surface boundary condition: ocean fields
29   USE sbcwave, ONLY   :  cdn_wave ! wave module
30#if defined key_si3 || defined key_cice
31   USE sbc_ice         ! Surface boundary condition: ice fields
32#endif
33   !
34   USE iom             ! I/O manager library
35   USE lib_mpp         ! distribued memory computing library
36   USE in_out_manager  ! I/O manager
37   USE prtctl          ! Print control
38   USE lib_fortran     ! to use key_nosignedzero
39
40   USE sbcblk_phy      ! all thermodynamics functions, rho_air, q_sat, etc... !LB
41
42   IMPLICIT NONE
43   PRIVATE
44
45   PUBLIC :: TURB_NCAR   ! called by sbcblk.F90
46
47   INTEGER , PARAMETER ::   nb_itt = 5        ! number of itterations
48   !! * Substitutions
49#  include "do_loop_substitute.h90"
50
51   !!----------------------------------------------------------------------
52CONTAINS
53
54   SUBROUTINE turb_ncar( zt, zu, sst, t_zt, ssq, q_zt, U_zu, &
55      &                  Cd, Ch, Ce, t_zu, q_zu, U_blk,      &
56      &                  Cdn, Chn, Cen                       )
57      !!----------------------------------------------------------------------------------
58      !!                      ***  ROUTINE  turb_ncar  ***
59      !!
60      !! ** Purpose :   Computes turbulent transfert coefficients of surface
61      !!                fluxes according to Large & Yeager (2004) and Large & Yeager (2008)
62      !!                If relevant (zt /= zu), adjust temperature and humidity from height zt to zu
63      !!                Returns the effective bulk wind speed at 10m to be used in the bulk formulas
64      !!
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      !!
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      !!    *  U_blk  : 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(wp), 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(wp), 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(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   U_zu     ! relative wind module at zu                [m/s]
96      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cd       ! transfer coefficient for momentum         (tau)
97      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ch       ! transfer coefficient for sensible heat (Q_sens)
98      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ce       ! transfert coefficient for evaporation   (Q_lat)
99      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   t_zu     ! pot. air temp. adjusted at zu               [K]
100      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   q_zu     ! spec. humidity adjusted at zu           [kg/kg]
101      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   U_blk    ! bulk wind speed at zu                     [m/s]
102      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients
103      !
104      INTEGER :: j_itt
105      LOGICAL :: l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U
106      !
107      REAL(wp), DIMENSION(jpi,jpj) ::   Cx_n10        ! 10m neutral latent/sensible coefficient
108      REAL(wp), DIMENSION(jpi,jpj) ::   sqrt_Cd_n10   ! root square of Cd_n10
109      REAL(wp), DIMENSION(jpi,jpj) ::   zeta_u        ! stability parameter at height zu
110      REAL(wp), DIMENSION(jpi,jpj) ::   zpsi_h_u
111      REAL(wp), DIMENSION(jpi,jpj) ::   ztmp0, ztmp1, ztmp2
112      REAL(wp), DIMENSION(jpi,jpj) ::   stab          ! stability test integer
113      !!----------------------------------------------------------------------------------
114      l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp ) ! testing "zu == zt" is risky with double precision
115
116      U_blk = MAX( 0.5_wp , U_zu )   !  relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s
117
118      !! First guess of stability:
119      ztmp0 = virt_temp(t_zt, q_zt) - virt_temp(sst, ssq) ! air-sea difference of virtual pot. temp. at zt
120      stab  = 0.5_wp + sign(0.5_wp,ztmp0)                           ! stab = 1 if dTv > 0  => STABLE, 0 if unstable
121
122      !! Neutral coefficients at 10m:
123      IF( ln_cdgw ) THEN      ! wave drag case
124         cdn_wave(:,:) = cdn_wave(:,:) + rsmall * ( 1._wp - tmask(:,:,1) )
125         ztmp0   (:,:) = cdn_wave(:,:)
126      ELSE
127      ztmp0 = cd_neutral_10m( U_blk )
128      ENDIF
129
130      sqrt_Cd_n10 = SQRT( ztmp0 )
131
132      !! Initializing transf. coeff. with their first guess neutral equivalents :
133      Cd = ztmp0
134      Ce = 1.e-3_wp*( 34.6_wp * sqrt_Cd_n10 )
135      Ch = 1.e-3_wp*sqrt_Cd_n10*(18._wp*stab + 32.7_wp*(1._wp - stab))
136      stab = sqrt_Cd_n10   ! Temporaty array !!! stab == SQRT(Cd)
137 
138      IF( ln_cdgw ) THEN
139   Cen = Ce
140   Chn = Ch
141      ENDIF
142
143      !! First guess of temperature and humidity at height zu:
144      t_zu = MAX( t_zt ,  180._wp )   ! who knows what's given on masked-continental regions...
145      q_zu = MAX( q_zt , 1.e-6_wp )   !               "
146
147      !! ITERATION BLOCK
148      DO j_itt = 1, nb_itt
149         !
150         ztmp1 = t_zu - sst   ! Updating air/sea differences
151         ztmp2 = q_zu - ssq
152
153         ! Updating turbulent scales :   (L&Y 2004 eq. (7))
154         ztmp0 = stab*U_blk       ! u*       (stab == SQRT(Cd))
155         ztmp1 = Ch/stab*ztmp1    ! theta*   (stab == SQRT(Cd))
156         ztmp2 = Ce/stab*ztmp2    ! q*       (stab == SQRT(Cd))
157
158         ! Estimate the inverse of Monin-Obukov length (1/L) at height zu:
159         ztmp0 = One_on_L( t_zu, q_zu, ztmp0, ztmp1, ztmp2 )
160         
161         !! Stability parameters :
162         zeta_u   = zu*ztmp0
163         zeta_u = sign( min(abs(zeta_u),10._wp), zeta_u )
164         zpsi_h_u = psi_h( zeta_u )
165
166         !! Shifting temperature and humidity at zu (L&Y 2004 eq. (9b-9c))
167         IF( .NOT. l_zt_equal_zu ) THEN
168            !! Array 'stab' is free for the moment so using it to store 'zeta_t'
169            stab = zt*ztmp0
170            stab = SIGN( MIN(ABS(stab),10._wp), stab )  ! Temporaty array stab == zeta_t !!!
171            stab = LOG(zt/zu) + zpsi_h_u - psi_h(stab)                   ! stab just used as temp array again!
172            t_zu = t_zt - ztmp1/vkarmn*stab    ! ztmp1 is still theta*  L&Y 2004 eq.(9b)
173            q_zu = q_zt - ztmp2/vkarmn*stab    ! ztmp2 is still q*      L&Y 2004 eq.(9c)
174            q_zu = max(0._wp, q_zu)
175         ENDIF
176
177         ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 eq. 9a)...
178         !   In very rare low-wind conditions, the old way of estimating the
179         !   neutral wind speed at 10m leads to a negative value that causes the code
180         !   to crash. To prevent this a threshold of 0.25m/s is imposed.
181         ztmp2 = psi_m(zeta_u)
182         IF( ln_cdgw ) THEN      ! surface wave case
183            stab = vkarmn / ( vkarmn / sqrt_Cd_n10 - ztmp2 )  ! (stab == SQRT(Cd))
184            Cd   = stab * stab
185            ztmp0 = (LOG(zu/10._wp) - zpsi_h_u) / vkarmn / sqrt_Cd_n10
186            ztmp2 = stab / sqrt_Cd_n10   ! (stab == SQRT(Cd))
187            ztmp1 = 1._wp + Chn * ztmp0     
188            Ch    = Chn * ztmp2 / ztmp1  ! L&Y 2004 eq. (10b)
189            ztmp1 = 1._wp + Cen * ztmp0
190            Ce    = Cen * ztmp2 / ztmp1  ! L&Y 2004 eq. (10c)
191
192         ELSE
193         ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 eq. 9a)...
194         !   In very rare low-wind conditions, the old way of estimating the
195         !   neutral wind speed at 10m leads to a negative value that causes the code
196         !   to crash. To prevent this a threshold of 0.25m/s is imposed.
197         ztmp0 = MAX( 0.25_wp , U_blk/(1._wp + sqrt_Cd_n10/vkarmn*(LOG(zu/10._wp) - ztmp2)) ) ! U_n10 (ztmp2 == psi_m(zeta_u))
198         ztmp0 = cd_neutral_10m(ztmp0)                                               ! Cd_n10
199         Cdn(:,:) = ztmp0
200         sqrt_Cd_n10 = sqrt(ztmp0)
201
202         stab    = 0.5_wp + sign(0.5_wp,zeta_u)                        ! update stability
203         Cx_n10  = 1.e-3_wp*sqrt_Cd_n10*(18._wp*stab + 32.7_wp*(1._wp - stab))  ! L&Y 2004 eq. (6c-6d)    (Cx_n10 == Ch_n10)
204         Chn(:,:) = Cx_n10
205
206         !! Update of transfer coefficients:
207         ztmp1 = 1._wp + sqrt_Cd_n10/vkarmn*(LOG(zu/10._wp) - ztmp2)   ! L&Y 2004 eq. (10a) (ztmp2 == psi_m(zeta_u))
208         Cd      = ztmp0 / ( ztmp1*ztmp1 )
209         stab = SQRT( Cd ) ! Temporary array !!! (stab == SQRT(Cd))
210
211         ztmp0 = (LOG(zu/10._wp) - zpsi_h_u) / vkarmn / sqrt_Cd_n10
212         ztmp2 = stab / sqrt_Cd_n10   ! (stab == SQRT(Cd))
213         ztmp1 = 1._wp + Cx_n10*ztmp0    ! (Cx_n10 == Ch_n10)
214         Ch  = Cx_n10*ztmp2 / ztmp1   ! L&Y 2004 eq. (10b)
215
216         Cx_n10  = 1.e-3_wp * (34.6_wp * sqrt_Cd_n10)  ! L&Y 2004 eq. (6b)    ! Cx_n10 == Ce_n10
217         Cen(:,:) = Cx_n10
218         ztmp1 = 1._wp + Cx_n10*ztmp0
219         Ce  = Cx_n10*ztmp2 / ztmp1  ! L&Y 2004 eq. (10c)
220         ENDIF
221
222      END DO !DO j_itt = 1, nb_itt
223
224   END SUBROUTINE turb_ncar
225
226
227   FUNCTION cd_neutral_10m( pw10 )
228      !!----------------------------------------------------------------------------------     
229      !! Estimate of the neutral drag coefficient at 10m as a function
230      !! of neutral wind  speed at 10m
231      !!
232      !! Origin: Large & Yeager 2008 eq.(11a) and eq.(11b)
233      !!
234      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
235      !!----------------------------------------------------------------------------------
236      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pw10           ! scalar wind speed at 10m (m/s)
237      REAL(wp), DIMENSION(jpi,jpj)             :: cd_neutral_10m
238      !
239      INTEGER  ::     ji, jj     ! dummy loop indices
240      REAL(wp) :: zgt33, zw, zw6 ! local scalars
241      !!----------------------------------------------------------------------------------
242      !
243      DO_2D_11_11
244         !
245         zw  = pw10(ji,jj)
246         zw6 = zw*zw*zw
247         zw6 = zw6*zw6
248         !
249         ! When wind speed > 33 m/s => Cyclone conditions => special treatment
250         zgt33 = 0.5_wp + SIGN( 0.5_wp, (zw - 33._wp) )   ! If pw10 < 33. => 0, else => 1
251         !
252         cd_neutral_10m(ji,jj) = 1.e-3_wp * ( &
253            &       (1._wp - zgt33)*( 2.7_wp/zw + 0.142_wp + zw/13.09_wp - 3.14807E-10_wp*zw6) & ! wind <  33 m/s
254            &      +    zgt33   *      2.34_wp )                                                 ! wind >= 33 m/s
255         !
256         cd_neutral_10m(ji,jj) = MAX(cd_neutral_10m(ji,jj), 1.E-6_wp)
257         !
258      END_2D
259      !
260   END FUNCTION cd_neutral_10m
261
262
263   FUNCTION psi_m( pzeta )
264      !!----------------------------------------------------------------------------------
265      !! Universal profile stability function for momentum
266      !!    !! Psis, L&Y 2004 eq. (8c), (8d), (8e)
267      !!
268      !! pzeta : stability paramenter, z/L where z is altitude measurement
269      !!         and L is M-O length
270      !!
271      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
272      !!----------------------------------------------------------------------------------
273      REAL(wp), DIMENSION(jpi,jpj) :: psi_m
274      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
275      !
276      INTEGER  ::   ji, jj    ! dummy loop indices
277      REAL(wp) :: zx2, zx, zstab   ! local scalars
278      !!----------------------------------------------------------------------------------
279      DO_2D_11_11
280         zx2 = SQRT( ABS( 1._wp - 16._wp*pzeta(ji,jj) ) )
281         zx2 = MAX( zx2 , 1._wp )
282         zx  = SQRT( zx2 )
283         zstab = 0.5_wp + SIGN( 0.5_wp , pzeta(ji,jj) )
284         !
285         psi_m(ji,jj) =        zstab  * (-5._wp*pzeta(ji,jj))       &          ! Stable
286            &          + (1._wp - zstab) * (2._wp*LOG((1._wp + zx)*0.5_wp)   &          ! Unstable
287            &               + LOG((1._wp + zx2)*0.5_wp) - 2._wp*ATAN(zx) + rpi*0.5_wp)  !    "
288         !
289      END_2D
290   END FUNCTION psi_m
291
292
293   FUNCTION psi_h( pzeta )
294      !!----------------------------------------------------------------------------------
295      !! Universal profile stability function for temperature and humidity
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_h
304      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
305      !
306      INTEGER  ::   ji, jj     ! dummy loop indices
307      REAL(wp) :: zx2, zstab  ! local scalars
308      !!----------------------------------------------------------------------------------
309      !
310      DO_2D_11_11
311         zx2 = SQRT( ABS( 1._wp - 16._wp*pzeta(ji,jj) ) )
312         zx2 = MAX( zx2 , 1._wp )
313         zstab = 0.5_wp + SIGN( 0.5_wp , pzeta(ji,jj) )
314         !
315         psi_h(ji,jj) =         zstab  * (-5._wp*pzeta(ji,jj))        &  ! Stable
316            &           + (1._wp - zstab) * (2._wp*LOG( (1._wp + zx2)*0.5_wp ))   ! Unstable
317         !
318      END_2D
319   END FUNCTION psi_h
320
321   !!======================================================================
322END MODULE sbcblk_algo_ncar
Note: See TracBrowser for help on using the repository browser.