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 @ 11111

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

LB: skin temperature now used by ECMWF bulk algorithm, and can be xios-ed ; new module "SBC/sbcblk_skin.F90" ; all physical constants defined in SBC now moved to "DOM/phycst.F90"; all air-properties and other ABL functions gathered in a new module: "SBC/sbcblk_phymbl.F90".
  • Property svn:keywords set to Id
File size: 16.4 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
41   IMPLICIT NONE
42   PRIVATE
43
44   PUBLIC ::   TURB_NCAR   ! called by sbcblk.F90
45   
46   !!----------------------------------------------------------------------
47CONTAINS
48
49   SUBROUTINE turb_ncar( zt, zu, sst, t_zt, ssq, q_zt, U_zu, &
50      &                  Cd, Ch, Ce, t_zu, q_zu, U_blk,      &
51      &                  Cdn, Chn, Cen                       )
52      !!----------------------------------------------------------------------------------
53      !!                      ***  ROUTINE  turb_ncar  ***
54      !!
55      !! ** Purpose :   Computes turbulent transfert coefficients of surface
56      !!                fluxes according to Large & Yeager (2004) and Large & Yeager (2008)
57      !!                If relevant (zt /= zu), adjust temperature and humidity from height zt to zu
58      !!                Returns the effective bulk wind speed at 10m to be used in the bulk formulas
59      !!
60      !! ** Method : Monin Obukhov Similarity Theory
61      !!             + Large & Yeager (2004,2008) closure: CD_n10 = f(U_n10)
62      !!
63      !! ** References :   Large & Yeager, 2004 / Large & Yeager, 2008
64      !!
65      !! ** Last update: Laurent Brodeau, June 2014:
66      !!    - handles both cases zt=zu and zt/=zu
67      !!    - optimized: less 2D arrays allocated and less operations
68      !!    - better first guess of stability by checking air-sea difference of virtual temperature
69      !!       rather than temperature difference only...
70      !!    - added function "cd_neutral_10m" that uses the improved parametrization of
71      !!      Large & Yeager 2008. Drag-coefficient reduction for Cyclone conditions!
72      !!    - using code-wide physical constants defined into "phycst.mod" rather than redifining them
73      !!      => 'vkarmn' and 'grav'
74      !!
75      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
76      !!
77      !! INPUT :
78      !! -------
79      !!    *  zt   : height for temperature and spec. hum. of air            [m]
80      !!    *  zu   : height for wind speed (generally 10m)                   [m]
81      !!    *  U_zu : scalar wind speed at 10m                                [m/s]
82      !!    *  sst  : SST                                                     [K]
83      !!    *  t_zt : potential air temperature at zt                         [K]
84      !!    *  ssq  : specific humidity at saturation at SST                  [kg/kg]
85      !!    *  q_zt : specific humidity of air at zt                          [kg/kg]
86      !!
87      !!
88      !! OUTPUT :
89      !! --------
90      !!    *  Cd     : drag coefficient
91      !!    *  Ch     : sensible heat coefficient
92      !!    *  Ce     : evaporation coefficient
93      !!    *  t_zu   : pot. air temperature adjusted at wind height zu       [K]
94      !!    *  q_zu   : specific humidity of air        //                    [kg/kg]
95      !!    *  U_blk  : bulk wind at 10m                                      [m/s]
96      !!----------------------------------------------------------------------------------
97      REAL(wp), INTENT(in   )                     ::   zt       ! height for t_zt and q_zt                    [m]
98      REAL(wp), INTENT(in   )                     ::   zu       ! height for U_zu                             [m]
99      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   sst      ! sea surface temperature                [Kelvin]
100      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   t_zt     ! potential air temperature              [Kelvin]
101      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   ssq      ! sea surface specific humidity           [kg/kg]
102      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_zt     ! specific air humidity                   [kg/kg]
103      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   U_zu     ! relative wind module at zu                [m/s]
104      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cd       ! transfer coefficient for momentum         (tau)
105      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ch       ! transfer coefficient for sensible heat (Q_sens)
106      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ce       ! transfert coefficient for evaporation   (Q_lat)
107      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   t_zu     ! pot. air temp. adjusted at zu               [K]
108      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   q_zu     ! spec. humidity adjusted at zu           [kg/kg]
109      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   U_blk    ! bulk wind at 10m                          [m/s]
110      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients
111      !
112      INTEGER ::   j_itt
113      LOGICAL ::   l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U
114      INTEGER , PARAMETER ::   nb_itt = 4       ! number of itterations
115      !
116      REAL(wp), DIMENSION(jpi,jpj) ::   Cx_n10        ! 10m neutral latent/sensible coefficient
117      REAL(wp), DIMENSION(jpi,jpj) ::   sqrt_Cd_n10   ! root square of Cd_n10
118      REAL(wp), DIMENSION(jpi,jpj) ::   zeta_u        ! stability parameter at height zu
119      REAL(wp), DIMENSION(jpi,jpj) ::   zpsi_h_u
120      REAL(wp), DIMENSION(jpi,jpj) ::   ztmp0, ztmp1, ztmp2
121      REAL(wp), DIMENSION(jpi,jpj) ::   stab          ! stability test integer
122      !!----------------------------------------------------------------------------------
123      !
124      l_zt_equal_zu = .FALSE.
125      IF( ABS(zu - zt) < 0.01 ) l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision
126
127      U_blk = MAX( 0.5 , U_zu )   !  relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s
128
129      !! First guess of stability:
130      ztmp0 = t_zt*(1. + rctv0*q_zt) - sst*(1. + rctv0*ssq) ! air-sea difference of virtual pot. temp. at zt
131      stab  = 0.5 + sign(0.5,ztmp0)                           ! stab = 1 if dTv > 0  => STABLE, 0 if unstable
132
133      !! Neutral coefficients at 10m:
134      IF( ln_cdgw ) THEN      ! wave drag case
135         cdn_wave(:,:) = cdn_wave(:,:) + rsmall * ( 1._wp - tmask(:,:,1) )
136         ztmp0   (:,:) = cdn_wave(:,:)
137      ELSE
138         ztmp0 = cd_neutral_10m( U_blk )
139      ENDIF
140
141      sqrt_Cd_n10 = SQRT( ztmp0 )
142
143      !! Initializing transf. coeff. with their first guess neutral equivalents :
144      Cd = ztmp0
145      Ce = 1.e-3*( 34.6 * sqrt_Cd_n10 )
146      Ch = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1. - stab))
147      stab = sqrt_Cd_n10   ! Temporaty array !!! stab == SQRT(Cd)
148 
149      IF( ln_cdgw )   Cen = Ce  ; Chn = Ch
150
151      !! Initializing values at z_u with z_t values:
152      t_zu = t_zt   ;   q_zu = q_zt
153
154      !!  * Now starting iteration loop
155      DO j_itt=1, nb_itt
156         !
157         ztmp1 = t_zu - sst   ! Updating air/sea differences
158         ztmp2 = q_zu - ssq
159
160         ! Updating turbulent scales :   (L&Y 2004 eq. (7))
161         ztmp1  = Ch/stab*ztmp1    ! theta*   (stab == SQRT(Cd))
162         ztmp2  = Ce/stab*ztmp2    ! q*       (stab == SQRT(Cd))
163
164         ztmp0 = 1. + rctv0*q_zu      ! multiply this with t and you have the virtual temperature
165
166         ! Estimate the inverse of Monin-Obukov length (1/L) at height zu:
167         ztmp0 =  (grav*vkarmn/(t_zu*ztmp0)*(ztmp1*ztmp0 + rctv0*t_zu*ztmp2)) / (Cd*U_blk*U_blk)
168         !                                                      ( Cd*U_blk*U_blk is U*^2 at zu )
169
170         !! Stability parameters :
171         zeta_u   = zu*ztmp0   ;  zeta_u = sign( min(abs(zeta_u),10.0), zeta_u )
172         zpsi_h_u = psi_h( zeta_u )
173
174         !! Shifting temperature and humidity at zu (L&Y 2004 eq. (9b-9c))
175         IF( .NOT. l_zt_equal_zu ) THEN
176            !! Array 'stab' is free for the moment so using it to store 'zeta_t'
177            stab = zt*ztmp0 ;  stab = SIGN( MIN(ABS(stab),10.0), stab )  ! Temporaty array stab == zeta_t !!!
178            stab = LOG(zt/zu) + zpsi_h_u - psi_h(stab)                   ! stab just used as temp array again!
179            t_zu = t_zt - ztmp1/vkarmn*stab    ! ztmp1 is still theta*  L&Y 2004 eq.(9b)
180            q_zu = q_zt - ztmp2/vkarmn*stab    ! ztmp2 is still q*      L&Y 2004 eq.(9c)
181            q_zu = max(0., q_zu)
182         END IF
183
184         ztmp2 = psi_m(zeta_u)
185         IF( ln_cdgw ) THEN      ! surface wave case
186            stab = vkarmn / ( vkarmn / sqrt_Cd_n10 - ztmp2 )  ! (stab == SQRT(Cd))
187            Cd   = stab * stab
188            ztmp0 = (LOG(zu/10.) - zpsi_h_u) / vkarmn / sqrt_Cd_n10
189            ztmp2 = stab / sqrt_Cd_n10   ! (stab == SQRT(Cd))
190            ztmp1 = 1. + Chn * ztmp0     
191            Ch    = Chn * ztmp2 / ztmp1  ! L&Y 2004 eq. (10b)
192            ztmp1 = 1. + Cen * ztmp0
193            Ce    = Cen * ztmp2 / ztmp1  ! L&Y 2004 eq. (10c)
194
195         ELSE
196            ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 eq. 9a)...
197            !   In very rare low-wind conditions, the old way of estimating the
198            !   neutral wind speed at 10m leads to a negative value that causes the code
199            !   to crash. To prevent this a threshold of 0.25m/s is imposed.
200            ztmp0 = MAX( 0.25 , U_blk/(1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - ztmp2)) ) ! U_n10 (ztmp2 == psi_m(zeta_u))
201            ztmp0 = cd_neutral_10m(ztmp0)                                               ! Cd_n10
202            Cdn(:,:) = ztmp0
203            sqrt_Cd_n10 = sqrt(ztmp0)
204
205            stab    = 0.5 + sign(0.5,zeta_u)                           ! update stability
206            Cx_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1. - stab))  ! L&Y 2004 eq. (6c-6d)    (Cx_n10 == Ch_n10)
207            Chn(:,:) = Cx_n10
208
209            !! Update of transfer coefficients:
210            ztmp1 = 1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - ztmp2)   ! L&Y 2004 eq. (10a) (ztmp2 == psi_m(zeta_u))
211            Cd      = ztmp0 / ( ztmp1*ztmp1 )
212            stab = SQRT( Cd ) ! Temporary array !!! (stab == SQRT(Cd))
213
214            ztmp0 = (LOG(zu/10.) - zpsi_h_u) / vkarmn / sqrt_Cd_n10
215            ztmp2 = stab / sqrt_Cd_n10   ! (stab == SQRT(Cd))
216            ztmp1 = 1. + Cx_n10*ztmp0    ! (Cx_n10 == Ch_n10)
217            Ch  = Cx_n10*ztmp2 / ztmp1   ! L&Y 2004 eq. (10b)
218
219            Cx_n10  = 1.e-3 * (34.6 * sqrt_Cd_n10)  ! L&Y 2004 eq. (6b)    ! Cx_n10 == Ce_n10
220            Cen(:,:) = Cx_n10
221            ztmp1 = 1. + Cx_n10*ztmp0
222            Ce  = Cx_n10*ztmp2 / ztmp1  ! L&Y 2004 eq. (10c)
223            ENDIF
224         !
225      END DO
226      !
227   END SUBROUTINE turb_ncar
228
229
230   FUNCTION cd_neutral_10m( pw10 )
231      !!----------------------------------------------------------------------------------     
232      !! Estimate of the neutral drag coefficient at 10m as a function
233      !! of neutral wind  speed at 10m
234      !!
235      !! Origin: Large & Yeager 2008 eq.(11a) and eq.(11b)
236      !!
237      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
238      !!----------------------------------------------------------------------------------
239      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pw10           ! scalar wind speed at 10m (m/s)
240      REAL(wp), DIMENSION(jpi,jpj)             :: cd_neutral_10m
241      !
242      INTEGER  ::     ji, jj     ! dummy loop indices
243      REAL(wp) :: zgt33, zw, zw6 ! local scalars
244      !!----------------------------------------------------------------------------------
245      !
246      DO jj = 1, jpj
247         DO ji = 1, jpi
248            !
249            zw  = pw10(ji,jj)
250            zw6 = zw*zw*zw
251            zw6 = zw6*zw6
252            !
253            ! When wind speed > 33 m/s => Cyclone conditions => special treatment
254            zgt33 = 0.5 + SIGN( 0.5, (zw - 33.) )   ! If pw10 < 33. => 0, else => 1
255            !
256            cd_neutral_10m(ji,jj) = 1.e-3 * ( &
257               &       (1. - zgt33)*( 2.7/zw + 0.142 + zw/13.09 - 3.14807E-10*zw6) & ! wind <  33 m/s
258               &      +    zgt33   *      2.34 )                                     ! wind >= 33 m/s
259            !
260            cd_neutral_10m(ji,jj) = MAX(cd_neutral_10m(ji,jj), 1.E-6)
261            !
262         END DO
263      END DO
264      !
265   END FUNCTION cd_neutral_10m
266
267
268   FUNCTION psi_m( pzeta )
269      !!----------------------------------------------------------------------------------
270      !! Universal profile stability function for momentum
271      !!    !! Psis, L&Y 2004 eq. (8c), (8d), (8e)
272      !!     
273      !! pzet0 : stability paramenter, z/L where z is altitude measurement                                         
274      !!         and L is M-O length
275      !!
276      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
277      !!----------------------------------------------------------------------------------
278      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pzeta
279      REAL(wp), DIMENSION(jpi,jpj)             ::   psi_m
280      !
281      INTEGER  ::   ji, jj         ! dummy loop indices
282      REAL(wp) :: zx2, zx, zstab   ! local scalars
283      !!----------------------------------------------------------------------------------
284      !
285      DO jj = 1, jpj
286         DO ji = 1, jpi
287            zx2 = SQRT( ABS( 1. - 16.*pzeta(ji,jj) ) )
288            zx2 = MAX ( zx2 , 1. )
289            zx  = SQRT( zx2 )
290            zstab = 0.5 + SIGN( 0.5 , pzeta(ji,jj) )
291            !
292            psi_m(ji,jj) =        zstab  * (-5.*pzeta(ji,jj))       &          ! Stable
293               &          + (1. - zstab) * (2.*LOG((1. + zx)*0.5)   &          ! Unstable
294               &               + LOG((1. + zx2)*0.5) - 2.*ATAN(zx) + rpi*0.5)  !    "
295            !
296         END DO
297      END DO
298      !
299   END FUNCTION psi_m
300
301
302   FUNCTION psi_h( pzeta )
303      !!----------------------------------------------------------------------------------
304      !! Universal profile stability function for temperature and humidity
305      !!    !! Psis, L&Y 2004 eq. (8c), (8d), (8e)
306      !!
307      !! pzet0 : stability paramenter, z/L where z is altitude measurement                                         
308      !!         and L is M-O length
309      !!
310      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
311      !!----------------------------------------------------------------------------------
312      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
313      REAL(wp), DIMENSION(jpi,jpj)             :: psi_h
314      !
315      INTEGER  ::   ji, jj    ! dummy loop indices
316      REAL(wp) :: zx2, zstab  ! local scalars
317      !!----------------------------------------------------------------------------------
318      !
319      DO jj = 1, jpj
320         DO ji = 1, jpi
321            zx2 = SQRT( ABS( 1. - 16.*pzeta(ji,jj) ) )
322            zx2 = MAX ( zx2 , 1. )
323            zstab = 0.5 + SIGN( 0.5 , pzeta(ji,jj) )
324            !
325            psi_h(ji,jj) =         zstab  * (-5.*pzeta(ji,jj))        &  ! Stable
326               &           + (1. - zstab) * (2.*LOG( (1. + zx2)*0.5 ))   ! Unstable
327            !
328         END DO
329      END DO
330      !
331   END FUNCTION psi_h
332
333   !!======================================================================
334END MODULE sbcblk_algo_ncar
Note: See TracBrowser for help on using the repository browser.