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_coare.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_coare.F90 @ 11209

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

LB: making more efficient and more centralized use of physical/thermodynamics functions defined into "OCE/SBC/sbcblk_phy.F90".
  • Property svn:keywords set to Id
File size: 18.4 KB
Line 
1MODULE sbcblk_algo_coare
2   !!======================================================================
3   !!                   ***  MODULE  sbcblk_algo_coare  ***
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 COARE v3, Fairall et al., 2003
11   !!
12   !!       Routine turb_coare maintained and developed in AeroBulk
13   !!                     (https://github.com/brodeau/aerobulk)
14   !!
15   !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk)
16   !!----------------------------------------------------------------------
17   !! History :  4.0  !  2016-02  (L.Brodeau)   Original code
18   !!----------------------------------------------------------------------
19
20   !!----------------------------------------------------------------------
21   !!   turb_coare  : 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   USE sbcblk_phy     !LB: all thermodynamics functions, rho_air, q_sat, etc... #LB
34   USE in_out_manager ! I/O manager
35   USE iom            ! I/O manager library
36   USE lib_mpp        ! distribued memory computing library
37   USE prtctl         ! Print control
38   USE lib_fortran    ! to use key_nosignedzero
39
40   IMPLICIT NONE
41   PRIVATE
42
43   PUBLIC ::   TURB_COARE   ! called by sbcblk.F90
44
45   !                                              !! COARE own values for given constants:
46   REAL(wp), PARAMETER ::   zi0     = 600._wp      ! scale height of the atmospheric boundary layer...
47   REAL(wp), PARAMETER ::   Beta0   =   1.250_wp   ! gustiness parameter
48
49   INTEGER , PARAMETER ::   nb_itt = 5             ! number of itterations
50
51   !!----------------------------------------------------------------------
52CONTAINS
53
54   SUBROUTINE turb_coare( 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_coare  ***
59      !!
60      !! ** Purpose :   Computes turbulent transfert coefficients of surface
61      !!                fluxes according to Fairall et al. (2003)
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 (generally 10m)                   [m]
70      !!    *  U_zu : scalar wind speed at 10m                                [m/s]
71      !!    *  sst  : SST                                                     [K]
72      !!    *  t_zt : potential air temperature at zt                         [K]
73      !!    *  ssq  : specific humidity at saturation at SST                  [kg/kg]
74      !!    *  q_zt : specific humidity of air at zt                          [kg/kg]
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 10m                                [m/s]
85      !!
86      !!
87      !! ** Author: L. Brodeau, june 2016 / 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 at 10m                          [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) ::  &
108         &  u_star, t_star, q_star, &
109         &  dt_zu, dq_zu,    &
110         &  znu_a,           & !: Nu_air, Viscosity of air
111         &  z0, z0t
112      REAL(wp), DIMENSION(jpi,jpj) ::   zeta_u        ! stability parameter at height zu
113      REAL(wp), DIMENSION(jpi,jpj) ::   ztmp0, ztmp1, ztmp2
114      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zeta_t        ! stability parameter at height zt
115      !!----------------------------------------------------------------------
116      !
117      l_zt_equal_zu = .FALSE.
118      IF( ABS(zu - zt) < 0.01 ) l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision
119
120      IF( .NOT. l_zt_equal_zu )  ALLOCATE( zeta_t(jpi,jpj) )
121
122      !! First guess of temperature and humidity at height zu:
123      t_zu = MAX( t_zt , 199.0_wp )   ! who knows what's given on masked-continental regions...
124      q_zu = MAX( q_zt , 1.e-6_wp )   !               "
125
126      !! Pot. temp. difference (and we don't want it to be 0!)
127      dt_zu = t_zu - sst ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
128      dq_zu = q_zu - ssq ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
129
130      znu_a = visc_air(t_zu) ! Air viscosity (m^2/s) at zt given from temperature in (K)
131
132      ztmp2 = 0.5*0.5  ! initial guess for wind gustiness contribution
133      U_blk = SQRT(U_zu*U_zu + ztmp2)
134
135      ztmp2   = 10000.     ! optimization: ztmp2 == 1/z0 (with z0 first guess == 0.0001)
136      ztmp0   = LOG(zu*ztmp2)
137      ztmp1   = LOG(10.*ztmp2)
138      u_star = 0.035*U_blk*ztmp1/ztmp0       ! (u* = 0.035*Un10)
139
140
141      z0     = alfa_charn(U_blk)*u_star*u_star/grav + 0.11*znu_a/u_star
142      z0     = MIN(ABS(z0), 0.001)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO
143      z0t    = 1. / ( 0.1*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) )
144      z0t    = MIN(ABS(z0t), 0.001)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO
145
146      ztmp2  = vkarmn/ztmp0
147      Cd     = ztmp2*ztmp2    ! first guess of Cd
148
149      ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd
150
151      ztmp2 = Ri_bulk( zu, sst, t_zu, ssq, q_zu, U_blk ) ! Bulk Richardson Number (BRN)
152      !ztmp2  = grav*zu*(dt_zu + rctv0*t_zu*dq_zu)/(t_zu*U_blk*U_blk)  !! Bulk Richardson number ;       !Ribcu = -zu/(zi0*0.004*Beta0**3) !! Saturation Rib, zi0 = tropicalbound. layer depth
153
154      !! First estimate of zeta_u, depending on the stability, ie sign of BRN (ztmp2):
155      ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 )
156      ztmp0 = ztmp0*ztmp2
157      zeta_u = (1.-ztmp1) * (ztmp0/(1.+ztmp2/(-zu/(zi0*0.004*Beta0**3)))) & !  BRN < 0
158         &  +     ztmp1   * (ztmp0*(1. + 27./9.*ztmp2/ztmp0))               !  BRN > 0
159      !#LOLO: should make sure that the "ztmp0" of "27./9.*ztmp2/ztmp0" is "ztmp0*ztmp2" and not "ztmp0==vkarmn*vkarmn/LOG(zt/z0t)/Cd" !
160
161      !! First guess M-O stability dependent scaling params.(u*,t*,q*) to estimate z0 and z/L
162      ztmp0  = vkarmn/(LOG(zu/z0t) - psi_h_coare(zeta_u))
163
164      u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0)  - psi_m_coare(zeta_u))
165      t_star = dt_zu*ztmp0
166      q_star = dq_zu*ztmp0
167
168      ! What's need to be done if zt /= zu:
169      IF( .NOT. l_zt_equal_zu ) THEN
170         !! First update of values at zu (or zt for wind)
171         zeta_t = zt*zeta_u/zu
172         ztmp0 = psi_h_coare(zeta_u) - psi_h_coare(zeta_t)
173         ztmp1 = LOG(zt/zu) + ztmp0
174         t_zu = t_zt - t_star/vkarmn*ztmp1
175         q_zu = q_zt - q_star/vkarmn*ztmp1
176         q_zu = (0.5 + SIGN(0.5_wp,q_zu))*q_zu !Makes it impossible to have negative humidity :
177         !
178         dt_zu = t_zu - sst  ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
179         dq_zu = q_zu - ssq  ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
180      END IF
181
182      !! ITERATION BLOCK
183      DO j_itt = 1, nb_itt
184
185         !!Inverse of Monin-Obukov length (1/L) :
186         ztmp0 = One_on_L(t_zu, q_zu, u_star, t_star, q_star)  ! 1/L == 1/[Monin-Obukhov length]
187         ztmp0 = SIGN( MIN(ABS(ztmp0),200._wp), ztmp0 ) ! (prevents FPE from stupid values from masked region later on...) !#LOLO
188
189         ztmp1 = u_star*u_star   ! u*^2
190
191         !! Update wind at 10m taking into acount convection-related wind gustiness:
192         ! Ug = Beta*w*  (Beta = 1.25, Fairall et al. 2003, Eq.8):
193         ztmp2 = Beta0*Beta0*ztmp1*(MAX(-zi0*ztmp0/vkarmn,0._wp))**(2./3.)   ! => ztmp2 == Ug^2
194         !!   ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before 600.
195         U_blk = MAX(sqrt(U_zu*U_zu + ztmp2), 0.2_wp)        ! include gustiness in bulk wind speed
196         ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0.
197
198         !! Updating Charnock parameter, increases with the wind (Fairall et al., 2003 p. 577-578)
199         ztmp2 = alfa_charn(U_blk)  ! alpha Charnock parameter
200
201         !! Roughness lengthes z0, z0t (z0q = z0t) :
202         z0   = ztmp2*ztmp1/grav + 0.11*znu_a/u_star ! Roughness length (eq.6)
203         ztmp1 = z0*u_star/znu_a                           ! Re_r: roughness Reynolds number
204         z0t  = MIN( 1.1E-4_wp , 5.5E-5_wp*ztmp1**(-0.6_wp) ) ! Scalar roughness for both theta and q (eq.28)
205
206         !! Stability parameters:
207         zeta_u = zu*ztmp0
208         zeta_u = SIGN( MIN(ABS(zeta_u),50.0_wp), zeta_u )
209         IF( .NOT. l_zt_equal_zu ) THEN
210            zeta_t = zt*ztmp0
211            zeta_t = SIGN( MIN(ABS(zeta_t),50.0_wp), zeta_t )
212         END IF
213
214         !! Turbulent scales at zu=10m :
215         ztmp0   = psi_h_coare(zeta_u)
216         ztmp1   = vkarmn/(LOG(zu) - LOG(z0t) - ztmp0)
217
218         t_star = dt_zu*ztmp1
219         q_star = dq_zu*ztmp1
220         u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0) - psi_m_coare(zeta_u))
221
222         IF( .NOT. l_zt_equal_zu ) THEN
223            ! What's need to be done if zt /= zu
224            !! Re-updating temperature and humidity at zu :
225            ztmp2 = ztmp0 - psi_h_coare(zeta_t)
226            ztmp1 = log(zt/zu) + ztmp2
227            t_zu = t_zt - t_star/vkarmn*ztmp1
228            q_zu = q_zt - q_star/vkarmn*ztmp1
229            dt_zu = t_zu - sst ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6), dt_zu )
230            dq_zu = q_zu - ssq ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9), dq_zu )
231         END IF
232
233      END DO
234      !
235      ! compute transfer coefficients at zu :
236      ztmp0 = u_star/U_blk
237      Cd   = ztmp0*ztmp0
238      Ch   = ztmp0*t_star/dt_zu
239      Ce   = ztmp0*q_star/dq_zu
240      !
241      ztmp1 = zu + z0
242      Cdn = vkarmn*vkarmn / (log(ztmp1/z0 )*log(ztmp1/z0 ))
243      Chn = vkarmn*vkarmn / (log(ztmp1/z0t)*log(ztmp1/z0t))
244      Cen = Chn
245      !
246      IF( .NOT. l_zt_equal_zu ) DEALLOCATE( zeta_t )
247      !
248   END SUBROUTINE turb_coare
249
250
251   FUNCTION alfa_charn( pwnd )
252      !!-------------------------------------------------------------------
253      !! Compute the Charnock parameter as a function of the wind speed
254      !!
255      !! (Fairall et al., 2003 p.577-578)
256      !!
257      !! Wind below 10 m/s :  alfa = 0.011
258      !! Wind between 10 and 18 m/s : linear increase from 0.011 to 0.018
259      !! Wind greater than 18 m/s :  alfa = 0.018
260      !!
261      !! Author: L. Brodeau, june 2016 / AeroBulk  (https://github.com/brodeau/aerobulk/)
262      !!-------------------------------------------------------------------
263      REAL(wp), DIMENSION(jpi,jpj) :: alfa_charn
264      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pwnd   ! wind speed
265      !
266      INTEGER  ::   ji, jj         ! dummy loop indices
267      REAL(wp) :: zw, zgt10, zgt18
268      !!-------------------------------------------------------------------
269      !
270      DO jj = 1, jpj
271         DO ji = 1, jpi
272            !
273            zw = pwnd(ji,jj)   ! wind speed
274            !
275            ! Charnock's constant, increases with the wind :
276            zgt10 = 0.5 + SIGN(0.5_wp,(zw - 10)) ! If zw<10. --> 0, else --> 1
277            zgt18 = 0.5 + SIGN(0.5_wp,(zw - 18.)) ! If zw<18. --> 0, else --> 1
278            !
279            alfa_charn(ji,jj) =  (1. - zgt10)*0.011    &    ! wind is lower than 10 m/s
280               &     + zgt10*((1. - zgt18)*(0.011 + (0.018 - 0.011) &
281               &      *(zw - 10.)/(18. - 10.)) + zgt18*( 0.018 ) )    ! Hare et al. (1999)
282            !
283         END DO
284      END DO
285      !
286   END FUNCTION alfa_charn
287
288   FUNCTION psi_m_coare( pzeta )
289      !!----------------------------------------------------------------------------------
290      !! ** Purpose: compute the universal profile stability function for momentum
291      !!             COARE 3.0, Fairall et al. 2003
292      !!             pzeta : stability paramenter, z/L where z is altitude
293      !!                     measurement and L is M-O length
294      !!       Stability function for wind speed and scalars matching Kansas and free
295      !!       convection forms with weighting f convective form, follows Fairall et
296      !!       al (1996) with profile constants from Grachev et al (2000) BLM stable
297      !!       form from Beljaars and Holtslag (1991)
298      !!
299      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
300      !!----------------------------------------------------------------------------------
301      REAL(wp), DIMENSION(jpi,jpj) :: psi_m_coare
302      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
303      !
304      INTEGER  ::   ji, jj    ! dummy loop indices
305      REAL(wp) :: zta, zphi_m, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
306      !!----------------------------------------------------------------------------------
307      !
308      DO jj = 1, jpj
309         DO ji = 1, jpi
310            !
311            zta = pzeta(ji,jj)
312            !
313            zphi_m = ABS(1. - 15.*zta)**.25    !!Kansas unstable
314            !
315            zpsi_k = 2.*LOG((1. + zphi_m)/2.) + LOG((1. + zphi_m*zphi_m)/2.)   &
316               & - 2.*ATAN(zphi_m) + 0.5*rpi
317            !
318            zphi_c = ABS(1. - 10.15*zta)**.3333                   !!Convective
319            !
320            zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
321               &     - 1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
322            !
323            zf = zta*zta
324            zf = zf/(1. + zf)
325            zc = MIN(50._wp, 0.35_wp*zta)
326            zstab = 0.5 + SIGN(0.5_wp, zta)
327            !
328            psi_m_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & ! (zta < 0)
329               &                -   zstab     * ( 1. + 1.*zta     &                ! (zta > 0)
330               &                         + 0.6667*(zta - 14.28)/EXP(zc) + 8.525 )   !     "
331            !
332         END DO
333      END DO
334      !
335   END FUNCTION psi_m_coare
336
337
338   FUNCTION psi_h_coare( pzeta )
339      !!---------------------------------------------------------------------
340      !! Universal profile stability function for temperature and humidity
341      !! COARE 3.0, Fairall et al. 2003
342      !!
343      !! pzeta : stability paramenter, z/L where z is altitude measurement
344      !!         and L is M-O length
345      !!
346      !! Stability function for wind speed and scalars matching Kansas and free
347      !! convection forms with weighting f convective form, follows Fairall et
348      !! al (1996) with profile constants from Grachev et al (2000) BLM stable
349      !! form from Beljaars and Holtslag (1991)
350      !!
351      !! Author: L. Brodeau, june 2016 / AeroBulk
352      !!         (https://github.com/brodeau/aerobulk/)
353      !!----------------------------------------------------------------
354      !!
355      REAL(wp), DIMENSION(jpi,jpj) :: psi_h_coare
356      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
357      !
358      INTEGER  ::   ji, jj     ! dummy loop indices
359      REAL(wp) :: zta, zphi_h, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
360      !
361      DO jj = 1, jpj
362         DO ji = 1, jpi
363            !
364            zta = pzeta(ji,jj)
365            !
366            zphi_h = (ABS(1. - 15.*zta))**.5  !! Kansas unstable   (zphi_h = zphi_m**2 when unstable, zphi_m when stable)
367            !
368            zpsi_k = 2.*LOG((1. + zphi_h)/2.)
369            !
370            zphi_c = (ABS(1. - 34.15*zta))**.3333   !! Convective
371            !
372            zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
373               &    -1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
374            !
375            zf = zta*zta
376            zf = zf/(1. + zf)
377            zc = MIN(50._wp,0.35_wp*zta)
378            zstab = 0.5 + SIGN(0.5_wp, zta)
379            !
380            psi_h_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) &
381               &                -   zstab     * ( (ABS(1. + 2.*zta/3.))**1.5     &
382               &                           + .6667*(zta - 14.28)/EXP(zc) + 8.525 )
383            !
384         END DO
385      END DO
386      !
387   END FUNCTION psi_h_coare
388
389   !!======================================================================
390END MODULE sbcblk_algo_coare
Note: See TracBrowser for help on using the repository browser.