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_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC – NEMO

source: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ncar.F90 @ 11291

Last change on this file since 11291 was 11221, checked in by laurent, 5 years ago

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