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

Last change on this file since 12015 was 12015, checked in by gsamson, 3 months ago

dev_ASINTER-01-05_merged: merge dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk branch @rev11988 with dev_r11265_ASINTER-01_Guillaume_ABL1D branch @rev11937 (tickets #2159 and #2131); ORCA2_ICE(_ABL) reproductibility OK

  • 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 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 = 5        ! 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 (usually 10m)                     [m]
67      !!    *  sst  : bulk SST                                                [K]
68      !!    *  t_zt : potential air temperature at zt                         [K]
69      !!    *  ssq  : specific humidity at saturation at SST                  [kg/kg]
70      !!    *  q_zt : specific humidity of air at zt                          [kg/kg]
71      !!    *  U_zu : scalar wind speed at zu                                 [m/s]
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 zu                                 [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 at zt             [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 speed at zu                     [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_wp )   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      !! ITERATION BLOCK
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         ENDIF
174
175         ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 eq. 9a)...
176         !   In very rare low-wind conditions, the old way of estimating the
177         !   neutral wind speed at 10m leads to a negative value that causes the code
178         !   to crash. To prevent this a threshold of 0.25m/s is imposed.
179         ztmp2 = psi_m(zeta_u)
180         IF( ln_cdgw ) THEN      ! surface wave case
181            stab = vkarmn / ( vkarmn / sqrt_Cd_n10 - ztmp2 )  ! (stab == SQRT(Cd))
182            Cd   = stab * stab
183            ztmp0 = (LOG(zu/10._wp) - zpsi_h_u) / vkarmn / sqrt_Cd_n10
184            ztmp2 = stab / sqrt_Cd_n10   ! (stab == SQRT(Cd))
185            ztmp1 = 1._wp + Chn * ztmp0     
186            Ch    = Chn * ztmp2 / ztmp1  ! L&Y 2004 eq. (10b)
187            ztmp1 = 1._wp + Cen * ztmp0
188            Ce    = Cen * ztmp2 / ztmp1  ! L&Y 2004 eq. (10c)
189
190         ELSE
191         ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 eq. 9a)...
192         !   In very rare low-wind conditions, the old way of estimating the
193         !   neutral wind speed at 10m leads to a negative value that causes the code
194         !   to crash. To prevent this a threshold of 0.25m/s is imposed.
195         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))
196         ztmp0 = cd_neutral_10m(ztmp0)                                               ! Cd_n10
197         Cdn(:,:) = ztmp0
198         sqrt_Cd_n10 = sqrt(ztmp0)
199
200         stab    = 0.5_wp + sign(0.5_wp,zeta_u)                        ! update stability
201         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)
202         Chn(:,:) = Cx_n10
203
204         !! Update of transfer coefficients:
205         ztmp1 = 1._wp + sqrt_Cd_n10/vkarmn*(LOG(zu/10._wp) - ztmp2)   ! L&Y 2004 eq. (10a) (ztmp2 == psi_m(zeta_u))
206         Cd      = ztmp0 / ( ztmp1*ztmp1 )
207         stab = SQRT( Cd ) ! Temporary array !!! (stab == SQRT(Cd))
208
209         ztmp0 = (LOG(zu/10._wp) - zpsi_h_u) / vkarmn / sqrt_Cd_n10
210         ztmp2 = stab / sqrt_Cd_n10   ! (stab == SQRT(Cd))
211         ztmp1 = 1._wp + Cx_n10*ztmp0    ! (Cx_n10 == Ch_n10)
212         Ch  = Cx_n10*ztmp2 / ztmp1   ! L&Y 2004 eq. (10b)
213
214         Cx_n10  = 1.e-3_wp * (34.6_wp * sqrt_Cd_n10)  ! L&Y 2004 eq. (6b)    ! Cx_n10 == Ce_n10
215         Cen(:,:) = Cx_n10
216         ztmp1 = 1._wp + Cx_n10*ztmp0
217         Ce  = Cx_n10*ztmp2 / ztmp1  ! L&Y 2004 eq. (10c)
218         ENDIF
219
220      END DO !DO j_itt = 1, nb_itt
221
222   END SUBROUTINE turb_ncar
223
224
225   FUNCTION cd_neutral_10m( pw10 )
226      !!----------------------------------------------------------------------------------     
227      !! Estimate of the neutral drag coefficient at 10m as a function
228      !! of neutral wind  speed at 10m
229      !!
230      !! Origin: Large & Yeager 2008 eq.(11a) and eq.(11b)
231      !!
232      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
233      !!----------------------------------------------------------------------------------
234      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pw10           ! scalar wind speed at 10m (m/s)
235      REAL(wp), DIMENSION(jpi,jpj)             :: cd_neutral_10m
236      !
237      INTEGER  ::     ji, jj     ! dummy loop indices
238      REAL(wp) :: zgt33, zw, zw6 ! local scalars
239      !!----------------------------------------------------------------------------------
240      !
241      DO jj = 1, jpj
242         DO ji = 1, jpi
243            !
244            zw  = pw10(ji,jj)
245            zw6 = zw*zw*zw
246            zw6 = zw6*zw6
247            !
248            ! When wind speed > 33 m/s => Cyclone conditions => special treatment
249            zgt33 = 0.5_wp + SIGN( 0.5_wp, (zw - 33._wp) )   ! If pw10 < 33. => 0, else => 1
250            !
251            cd_neutral_10m(ji,jj) = 1.e-3_wp * ( &
252               &       (1._wp - zgt33)*( 2.7_wp/zw + 0.142_wp + zw/13.09_wp - 3.14807E-10_wp*zw6) & ! wind <  33 m/s
253               &      +    zgt33   *      2.34_wp )                                                 ! wind >= 33 m/s
254            !
255            cd_neutral_10m(ji,jj) = MAX(cd_neutral_10m(ji,jj), 1.E-6_wp)
256            !
257         END DO
258      END DO
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 jj = 1, jpj
280         DO ji = 1, jpi
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 DO
291      END DO
292   END FUNCTION psi_m
293
294
295   FUNCTION psi_h( pzeta )
296      !!----------------------------------------------------------------------------------
297      !! Universal profile stability function for temperature and humidity
298      !!    !! Psis, L&Y 2004 eq. (8c), (8d), (8e)
299      !!
300      !! pzeta : stability paramenter, z/L where z is altitude measurement
301      !!         and L is M-O length
302      !!
303      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
304      !!----------------------------------------------------------------------------------
305      REAL(wp), DIMENSION(jpi,jpj) :: psi_h
306      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
307      !
308      INTEGER  ::   ji, jj     ! dummy loop indices
309      REAL(wp) :: zx2, zstab  ! local scalars
310      !!----------------------------------------------------------------------------------
311      !
312      DO jj = 1, jpj
313         DO ji = 1, jpi
314            zx2 = SQRT( ABS( 1._wp - 16._wp*pzeta(ji,jj) ) )
315            zx2 = MAX( zx2 , 1._wp )
316            zstab = 0.5_wp + SIGN( 0.5_wp , pzeta(ji,jj) )
317            !
318            psi_h(ji,jj) =         zstab  * (-5._wp*pzeta(ji,jj))        &  ! Stable
319               &           + (1._wp - zstab) * (2._wp*LOG( (1._wp + zx2)*0.5_wp ))   ! Unstable
320            !
321         END DO
322      END DO
323   END FUNCTION psi_h
324
325   !!======================================================================
326END MODULE sbcblk_algo_ncar
Note: See TracBrowser for help on using the repository browser.