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/2019/dev_ASINTER-01-05_merged/src/OCE/SBC – NEMO

source: NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/SBC/sbcblk_algo_ncar.F90 @ 12098

Last change on this file since 12098 was 12063, checked in by gsamson, 4 years ago

dev_ASINTER-01-05_merged: update branch with dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk@r12061 and trunk@r12055 + bugfix for agrif compatibility in sbcblk: sette tests with ref configs ok except ABL restartability (under investigation) (tickets #2159 and #2131)

  • Property svn:keywords set to Id
File size: 15.8 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
49   !!----------------------------------------------------------------------
50CONTAINS
51
52   SUBROUTINE turb_ncar( zt, zu, sst, t_zt, ssq, q_zt, U_zu, &
53      &                  Cd, Ch, Ce, t_zu, q_zu, U_blk,      &
54      &                  Cdn, Chn, Cen                       )
55      !!----------------------------------------------------------------------------------
56      !!                      ***  ROUTINE  turb_ncar  ***
57      !!
58      !! ** Purpose :   Computes turbulent transfert coefficients of surface
59      !!                fluxes according to Large & Yeager (2004) and Large & Yeager (2008)
60      !!                If relevant (zt /= zu), adjust temperature and humidity from height zt to zu
61      !!                Returns the effective bulk wind speed at 10m to be used in the bulk formulas
62      !!
63      !!
64      !! INPUT :
65      !! -------
66      !!    *  zt   : height for temperature and spec. hum. of air            [m]
67      !!    *  zu   : height for wind speed (usually 10m)                     [m]
68      !!    *  sst  : bulk SST                                                [K]
69      !!    *  t_zt : potential air temperature at zt                         [K]
70      !!    *  ssq  : specific humidity at saturation at SST                  [kg/kg]
71      !!    *  q_zt : specific humidity of air at zt                          [kg/kg]
72      !!    *  U_zu : scalar wind speed at zu                                 [m/s]
73      !!
74      !!
75      !! OUTPUT :
76      !! --------
77      !!    *  Cd     : drag coefficient
78      !!    *  Ch     : sensible heat coefficient
79      !!    *  Ce     : evaporation coefficient
80      !!    *  t_zu   : pot. air temperature adjusted at wind height zu       [K]
81      !!    *  q_zu   : specific humidity of air        //                    [kg/kg]
82      !!    *  U_blk  : bulk wind speed at zu                                 [m/s]
83      !!
84      !!
85      !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/)
86      !!----------------------------------------------------------------------------------
87      REAL(wp), INTENT(in   )                     ::   zt       ! height for t_zt and q_zt                    [m]
88      REAL(wp), INTENT(in   )                     ::   zu       ! height for U_zu                             [m]
89      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   sst      ! sea surface temperature                [Kelvin]
90      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   t_zt     ! potential air temperature              [Kelvin]
91      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   ssq      ! sea surface specific humidity           [kg/kg]
92      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_zt     ! specific air humidity at zt             [kg/kg]
93      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   U_zu     ! relative wind module at zu                [m/s]
94      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cd       ! transfer coefficient for momentum         (tau)
95      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ch       ! transfer coefficient for sensible heat (Q_sens)
96      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ce       ! transfert coefficient for evaporation   (Q_lat)
97      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   t_zu     ! pot. air temp. adjusted at zu               [K]
98      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   q_zu     ! spec. humidity adjusted at zu           [kg/kg]
99      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   U_blk    ! bulk wind speed at zu                     [m/s]
100      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients
101      !
102      INTEGER :: j_itt
103      LOGICAL :: l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U
104      !
105      REAL(wp), DIMENSION(jpi,jpj) ::   Cx_n10        ! 10m neutral latent/sensible coefficient
106      REAL(wp), DIMENSION(jpi,jpj) ::   sqrt_Cd_n10   ! root square of Cd_n10
107      REAL(wp), DIMENSION(jpi,jpj) ::   zeta_u        ! stability parameter at height zu
108      REAL(wp), DIMENSION(jpi,jpj) ::   zpsi_h_u
109      REAL(wp), DIMENSION(jpi,jpj) ::   ztmp0, ztmp1, ztmp2
110      REAL(wp), DIMENSION(jpi,jpj) ::   stab          ! stability test integer
111      !!----------------------------------------------------------------------------------
112      !
113      l_zt_equal_zu = .FALSE.
114      IF( ABS(zu - zt) < 0.01_wp )   l_zt_equal_zu = .TRUE.    ! 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      !! Initializing values at z_u with z_t values:
144      t_zu = t_zt   ;   q_zu = q_zt
145
146      !! ITERATION BLOCK
147      DO j_itt = 1, nb_itt
148         !
149         ztmp1 = t_zu - sst   ! Updating air/sea differences
150         ztmp2 = q_zu - ssq
151
152         ! Updating turbulent scales :   (L&Y 2004 eq. (7))
153         ztmp0 = stab*U_blk       ! u*       (stab == SQRT(Cd))
154         ztmp1 = Ch/stab*ztmp1    ! theta*   (stab == SQRT(Cd))
155         ztmp2 = Ce/stab*ztmp2    ! q*       (stab == SQRT(Cd))
156
157         ! Estimate the inverse of Monin-Obukov length (1/L) at height zu:
158         ztmp0 = One_on_L( t_zu, q_zu, ztmp0, ztmp1, ztmp2 )
159         
160         !! Stability parameters :
161         zeta_u   = zu*ztmp0
162         zeta_u = sign( min(abs(zeta_u),10._wp), zeta_u )
163         zpsi_h_u = psi_h( zeta_u )
164
165         !! Shifting temperature and humidity at zu (L&Y 2004 eq. (9b-9c))
166         IF( .NOT. l_zt_equal_zu ) THEN
167            !! Array 'stab' is free for the moment so using it to store 'zeta_t'
168            stab = zt*ztmp0
169            stab = SIGN( MIN(ABS(stab),10._wp), stab )  ! Temporaty array stab == zeta_t !!!
170            stab = LOG(zt/zu) + zpsi_h_u - psi_h(stab)                   ! stab just used as temp array again!
171            t_zu = t_zt - ztmp1/vkarmn*stab    ! ztmp1 is still theta*  L&Y 2004 eq.(9b)
172            q_zu = q_zt - ztmp2/vkarmn*stab    ! ztmp2 is still q*      L&Y 2004 eq.(9c)
173            q_zu = max(0._wp, q_zu)
174         ENDIF
175
176         ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 eq. 9a)...
177         !   In very rare low-wind conditions, the old way of estimating the
178         !   neutral wind speed at 10m leads to a negative value that causes the code
179         !   to crash. To prevent this a threshold of 0.25m/s is imposed.
180         ztmp2 = psi_m(zeta_u)
181         IF( ln_cdgw ) THEN      ! surface wave case
182            stab = vkarmn / ( vkarmn / sqrt_Cd_n10 - ztmp2 )  ! (stab == SQRT(Cd))
183            Cd   = stab * stab
184            ztmp0 = (LOG(zu/10._wp) - zpsi_h_u) / vkarmn / sqrt_Cd_n10
185            ztmp2 = stab / sqrt_Cd_n10   ! (stab == SQRT(Cd))
186            ztmp1 = 1._wp + Chn * ztmp0     
187            Ch    = Chn * ztmp2 / ztmp1  ! L&Y 2004 eq. (10b)
188            ztmp1 = 1._wp + Cen * ztmp0
189            Ce    = Cen * ztmp2 / ztmp1  ! L&Y 2004 eq. (10c)
190
191         ELSE
192         ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 eq. 9a)...
193         !   In very rare low-wind conditions, the old way of estimating the
194         !   neutral wind speed at 10m leads to a negative value that causes the code
195         !   to crash. To prevent this a threshold of 0.25m/s is imposed.
196         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))
197         ztmp0 = cd_neutral_10m(ztmp0)                                               ! Cd_n10
198         Cdn(:,:) = ztmp0
199         sqrt_Cd_n10 = sqrt(ztmp0)
200
201         stab    = 0.5_wp + sign(0.5_wp,zeta_u)                        ! update stability
202         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)
203         Chn(:,:) = Cx_n10
204
205         !! Update of transfer coefficients:
206         ztmp1 = 1._wp + sqrt_Cd_n10/vkarmn*(LOG(zu/10._wp) - ztmp2)   ! L&Y 2004 eq. (10a) (ztmp2 == psi_m(zeta_u))
207         Cd      = ztmp0 / ( ztmp1*ztmp1 )
208         stab = SQRT( Cd ) ! Temporary array !!! (stab == SQRT(Cd))
209
210         ztmp0 = (LOG(zu/10._wp) - zpsi_h_u) / vkarmn / sqrt_Cd_n10
211         ztmp2 = stab / sqrt_Cd_n10   ! (stab == SQRT(Cd))
212         ztmp1 = 1._wp + Cx_n10*ztmp0    ! (Cx_n10 == Ch_n10)
213         Ch  = Cx_n10*ztmp2 / ztmp1   ! L&Y 2004 eq. (10b)
214
215         Cx_n10  = 1.e-3_wp * (34.6_wp * sqrt_Cd_n10)  ! L&Y 2004 eq. (6b)    ! Cx_n10 == Ce_n10
216         Cen(:,:) = Cx_n10
217         ztmp1 = 1._wp + Cx_n10*ztmp0
218         Ce  = Cx_n10*ztmp2 / ztmp1  ! L&Y 2004 eq. (10c)
219         ENDIF
220
221      END DO !DO j_itt = 1, nb_itt
222
223   END SUBROUTINE turb_ncar
224
225
226   FUNCTION cd_neutral_10m( pw10 )
227      !!----------------------------------------------------------------------------------     
228      !! Estimate of the neutral drag coefficient at 10m as a function
229      !! of neutral wind  speed at 10m
230      !!
231      !! Origin: Large & Yeager 2008 eq.(11a) and eq.(11b)
232      !!
233      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
234      !!----------------------------------------------------------------------------------
235      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pw10           ! scalar wind speed at 10m (m/s)
236      REAL(wp), DIMENSION(jpi,jpj)             :: cd_neutral_10m
237      !
238      INTEGER  ::     ji, jj     ! dummy loop indices
239      REAL(wp) :: zgt33, zw, zw6 ! local scalars
240      !!----------------------------------------------------------------------------------
241      !
242      DO jj = 1, jpj
243         DO ji = 1, jpi
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 DO
259      END DO
260      !
261   END FUNCTION cd_neutral_10m
262
263
264   FUNCTION psi_m( pzeta )
265      !!----------------------------------------------------------------------------------
266      !! Universal profile stability function for momentum
267      !!    !! Psis, L&Y 2004 eq. (8c), (8d), (8e)
268      !!
269      !! pzeta : stability paramenter, z/L where z is altitude measurement
270      !!         and L is M-O length
271      !!
272      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
273      !!----------------------------------------------------------------------------------
274      REAL(wp), DIMENSION(jpi,jpj) :: psi_m
275      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
276      !
277      INTEGER  ::   ji, jj    ! dummy loop indices
278      REAL(wp) :: zx2, zx, zstab   ! local scalars
279      !!----------------------------------------------------------------------------------
280      DO jj = 1, jpj
281         DO ji = 1, jpi
282            zx2 = SQRT( ABS( 1._wp - 16._wp*pzeta(ji,jj) ) )
283            zx2 = MAX( zx2 , 1._wp )
284            zx  = SQRT( zx2 )
285            zstab = 0.5_wp + SIGN( 0.5_wp , pzeta(ji,jj) )
286            !
287            psi_m(ji,jj) =        zstab  * (-5._wp*pzeta(ji,jj))       &          ! Stable
288               &          + (1._wp - zstab) * (2._wp*LOG((1._wp + zx)*0.5_wp)   &          ! Unstable
289               &               + LOG((1._wp + zx2)*0.5_wp) - 2._wp*ATAN(zx) + rpi*0.5_wp)  !    "
290            !
291         END DO
292      END DO
293   END FUNCTION psi_m
294
295
296   FUNCTION psi_h( pzeta )
297      !!----------------------------------------------------------------------------------
298      !! Universal profile stability function for temperature and humidity
299      !!    !! Psis, L&Y 2004 eq. (8c), (8d), (8e)
300      !!
301      !! pzeta : stability paramenter, z/L where z is altitude measurement
302      !!         and L is M-O length
303      !!
304      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
305      !!----------------------------------------------------------------------------------
306      REAL(wp), DIMENSION(jpi,jpj) :: psi_h
307      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
308      !
309      INTEGER  ::   ji, jj     ! dummy loop indices
310      REAL(wp) :: zx2, zstab  ! local scalars
311      !!----------------------------------------------------------------------------------
312      !
313      DO jj = 1, jpj
314         DO ji = 1, jpi
315            zx2 = SQRT( ABS( 1._wp - 16._wp*pzeta(ji,jj) ) )
316            zx2 = MAX( zx2 , 1._wp )
317            zstab = 0.5_wp + SIGN( 0.5_wp , pzeta(ji,jj) )
318            !
319            psi_h(ji,jj) =         zstab  * (-5._wp*pzeta(ji,jj))        &  ! Stable
320               &           + (1._wp - zstab) * (2._wp*LOG( (1._wp + zx2)*0.5_wp ))   ! Unstable
321            !
322         END DO
323      END DO
324   END FUNCTION psi_h
325
326   !!======================================================================
327END MODULE sbcblk_algo_ncar
Note: See TracBrowser for help on using the repository browser.