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

source: NEMO/trunk/src/OCE/SBC/sbcblk_algo_ncar.F90 @ 12377

Last change on this file since 12377 was 12377, checked in by acc, 4 years ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.
  • 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      !
115      l_zt_equal_zu = .FALSE.
116      IF( ABS(zu - zt) < 0.01_wp )   l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision
117
118      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
119
120      !! First guess of stability:
121      ztmp0 = virt_temp(t_zt, q_zt) - virt_temp(sst, ssq) ! air-sea difference of virtual pot. temp. at zt
122      stab  = 0.5_wp + sign(0.5_wp,ztmp0)                           ! stab = 1 if dTv > 0  => STABLE, 0 if unstable
123
124      !! Neutral coefficients at 10m:
125      IF( ln_cdgw ) THEN      ! wave drag case
126         cdn_wave(:,:) = cdn_wave(:,:) + rsmall * ( 1._wp - tmask(:,:,1) )
127         ztmp0   (:,:) = cdn_wave(:,:)
128      ELSE
129      ztmp0 = cd_neutral_10m( U_blk )
130      ENDIF
131
132      sqrt_Cd_n10 = SQRT( ztmp0 )
133
134      !! Initializing transf. coeff. with their first guess neutral equivalents :
135      Cd = ztmp0
136      Ce = 1.e-3_wp*( 34.6_wp * sqrt_Cd_n10 )
137      Ch = 1.e-3_wp*sqrt_Cd_n10*(18._wp*stab + 32.7_wp*(1._wp - stab))
138      stab = sqrt_Cd_n10   ! Temporaty array !!! stab == SQRT(Cd)
139 
140      IF( ln_cdgw ) THEN
141   Cen = Ce
142   Chn = Ch
143      ENDIF
144
145      !! Initializing values at z_u with z_t values:
146      t_zu = t_zt   ;   q_zu = q_zt
147
148      !! ITERATION BLOCK
149      DO j_itt = 1, nb_itt
150         !
151         ztmp1 = t_zu - sst   ! Updating air/sea differences
152         ztmp2 = q_zu - ssq
153
154         ! Updating turbulent scales :   (L&Y 2004 eq. (7))
155         ztmp0 = stab*U_blk       ! u*       (stab == SQRT(Cd))
156         ztmp1 = Ch/stab*ztmp1    ! theta*   (stab == SQRT(Cd))
157         ztmp2 = Ce/stab*ztmp2    ! q*       (stab == SQRT(Cd))
158
159         ! Estimate the inverse of Monin-Obukov length (1/L) at height zu:
160         ztmp0 = One_on_L( t_zu, q_zu, ztmp0, ztmp1, ztmp2 )
161         
162         !! Stability parameters :
163         zeta_u   = zu*ztmp0
164         zeta_u = sign( min(abs(zeta_u),10._wp), zeta_u )
165         zpsi_h_u = psi_h( zeta_u )
166
167         !! Shifting temperature and humidity at zu (L&Y 2004 eq. (9b-9c))
168         IF( .NOT. l_zt_equal_zu ) THEN
169            !! Array 'stab' is free for the moment so using it to store 'zeta_t'
170            stab = zt*ztmp0
171            stab = SIGN( MIN(ABS(stab),10._wp), stab )  ! Temporaty array stab == zeta_t !!!
172            stab = LOG(zt/zu) + zpsi_h_u - psi_h(stab)                   ! stab just used as temp array again!
173            t_zu = t_zt - ztmp1/vkarmn*stab    ! ztmp1 is still theta*  L&Y 2004 eq.(9b)
174            q_zu = q_zt - ztmp2/vkarmn*stab    ! ztmp2 is still q*      L&Y 2004 eq.(9c)
175            q_zu = max(0._wp, q_zu)
176         ENDIF
177
178         ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 eq. 9a)...
179         !   In very rare low-wind conditions, the old way of estimating the
180         !   neutral wind speed at 10m leads to a negative value that causes the code
181         !   to crash. To prevent this a threshold of 0.25m/s is imposed.
182         ztmp2 = psi_m(zeta_u)
183         IF( ln_cdgw ) THEN      ! surface wave case
184            stab = vkarmn / ( vkarmn / sqrt_Cd_n10 - ztmp2 )  ! (stab == SQRT(Cd))
185            Cd   = stab * stab
186            ztmp0 = (LOG(zu/10._wp) - zpsi_h_u) / vkarmn / sqrt_Cd_n10
187            ztmp2 = stab / sqrt_Cd_n10   ! (stab == SQRT(Cd))
188            ztmp1 = 1._wp + Chn * ztmp0     
189            Ch    = Chn * ztmp2 / ztmp1  ! L&Y 2004 eq. (10b)
190            ztmp1 = 1._wp + Cen * ztmp0
191            Ce    = Cen * ztmp2 / ztmp1  ! L&Y 2004 eq. (10c)
192
193         ELSE
194         ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 eq. 9a)...
195         !   In very rare low-wind conditions, the old way of estimating the
196         !   neutral wind speed at 10m leads to a negative value that causes the code
197         !   to crash. To prevent this a threshold of 0.25m/s is imposed.
198         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))
199         ztmp0 = cd_neutral_10m(ztmp0)                                               ! Cd_n10
200         Cdn(:,:) = ztmp0
201         sqrt_Cd_n10 = sqrt(ztmp0)
202
203         stab    = 0.5_wp + sign(0.5_wp,zeta_u)                        ! update stability
204         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)
205         Chn(:,:) = Cx_n10
206
207         !! Update of transfer coefficients:
208         ztmp1 = 1._wp + sqrt_Cd_n10/vkarmn*(LOG(zu/10._wp) - ztmp2)   ! L&Y 2004 eq. (10a) (ztmp2 == psi_m(zeta_u))
209         Cd      = ztmp0 / ( ztmp1*ztmp1 )
210         stab = SQRT( Cd ) ! Temporary array !!! (stab == SQRT(Cd))
211
212         ztmp0 = (LOG(zu/10._wp) - zpsi_h_u) / vkarmn / sqrt_Cd_n10
213         ztmp2 = stab / sqrt_Cd_n10   ! (stab == SQRT(Cd))
214         ztmp1 = 1._wp + Cx_n10*ztmp0    ! (Cx_n10 == Ch_n10)
215         Ch  = Cx_n10*ztmp2 / ztmp1   ! L&Y 2004 eq. (10b)
216
217         Cx_n10  = 1.e-3_wp * (34.6_wp * sqrt_Cd_n10)  ! L&Y 2004 eq. (6b)    ! Cx_n10 == Ce_n10
218         Cen(:,:) = Cx_n10
219         ztmp1 = 1._wp + Cx_n10*ztmp0
220         Ce  = Cx_n10*ztmp2 / ztmp1  ! L&Y 2004 eq. (10c)
221         ENDIF
222
223      END DO !DO j_itt = 1, nb_itt
224
225   END SUBROUTINE turb_ncar
226
227
228   FUNCTION cd_neutral_10m( pw10 )
229      !!----------------------------------------------------------------------------------     
230      !! Estimate of the neutral drag coefficient at 10m as a function
231      !! of neutral wind  speed at 10m
232      !!
233      !! Origin: Large & Yeager 2008 eq.(11a) and eq.(11b)
234      !!
235      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
236      !!----------------------------------------------------------------------------------
237      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pw10           ! scalar wind speed at 10m (m/s)
238      REAL(wp), DIMENSION(jpi,jpj)             :: cd_neutral_10m
239      !
240      INTEGER  ::     ji, jj     ! dummy loop indices
241      REAL(wp) :: zgt33, zw, zw6 ! local scalars
242      !!----------------------------------------------------------------------------------
243      !
244      DO_2D_11_11
245         !
246         zw  = pw10(ji,jj)
247         zw6 = zw*zw*zw
248         zw6 = zw6*zw6
249         !
250         ! When wind speed > 33 m/s => Cyclone conditions => special treatment
251         zgt33 = 0.5_wp + SIGN( 0.5_wp, (zw - 33._wp) )   ! If pw10 < 33. => 0, else => 1
252         !
253         cd_neutral_10m(ji,jj) = 1.e-3_wp * ( &
254            &       (1._wp - zgt33)*( 2.7_wp/zw + 0.142_wp + zw/13.09_wp - 3.14807E-10_wp*zw6) & ! wind <  33 m/s
255            &      +    zgt33   *      2.34_wp )                                                 ! wind >= 33 m/s
256         !
257         cd_neutral_10m(ji,jj) = MAX(cd_neutral_10m(ji,jj), 1.E-6_wp)
258         !
259      END_2D
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_2D_11_11
281         zx2 = SQRT( ABS( 1._wp - 16._wp*pzeta(ji,jj) ) )
282         zx2 = MAX( zx2 , 1._wp )
283         zx  = SQRT( zx2 )
284         zstab = 0.5_wp + SIGN( 0.5_wp , pzeta(ji,jj) )
285         !
286         psi_m(ji,jj) =        zstab  * (-5._wp*pzeta(ji,jj))       &          ! Stable
287            &          + (1._wp - zstab) * (2._wp*LOG((1._wp + zx)*0.5_wp)   &          ! Unstable
288            &               + LOG((1._wp + zx2)*0.5_wp) - 2._wp*ATAN(zx) + rpi*0.5_wp)  !    "
289         !
290      END_2D
291   END FUNCTION psi_m
292
293
294   FUNCTION psi_h( pzeta )
295      !!----------------------------------------------------------------------------------
296      !! Universal profile stability function for temperature and humidity
297      !!    !! Psis, L&Y 2004 eq. (8c), (8d), (8e)
298      !!
299      !! pzeta : stability paramenter, z/L where z is altitude measurement
300      !!         and L is M-O length
301      !!
302      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
303      !!----------------------------------------------------------------------------------
304      REAL(wp), DIMENSION(jpi,jpj) :: psi_h
305      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
306      !
307      INTEGER  ::   ji, jj     ! dummy loop indices
308      REAL(wp) :: zx2, zstab  ! local scalars
309      !!----------------------------------------------------------------------------------
310      !
311      DO_2D_11_11
312         zx2 = SQRT( ABS( 1._wp - 16._wp*pzeta(ji,jj) ) )
313         zx2 = MAX( zx2 , 1._wp )
314         zstab = 0.5_wp + SIGN( 0.5_wp , pzeta(ji,jj) )
315         !
316         psi_h(ji,jj) =         zstab  * (-5._wp*pzeta(ji,jj))        &  ! Stable
317            &           + (1._wp - zstab) * (2._wp*LOG( (1._wp + zx2)*0.5_wp ))   ! Unstable
318         !
319      END_2D
320   END FUNCTION psi_h
321
322   !!======================================================================
323END MODULE sbcblk_algo_ncar
Note: See TracBrowser for help on using the repository browser.