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

Last change on this file since 11209 was 11209, checked in by laurent, 15 months 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: 17.1 KB
Line 
1MODULE sbcblk_algo_coare3p5
2   !!======================================================================
3   !!                       ***  MODULE  sbcblk_algo_coare3p5  ***
4   !! Computes turbulent components of surface fluxes
5   !!         according to Edson et al. 2013 (COARE v3.5) /JPO
6   !!
7   !!   * bulk transfer coefficients C_D, C_E and C_H
8   !!   * air temp. and spec. hum. adjusted from zt (2m) to zu (10m) if needed
9   !!   * the effective bulk wind speed at 10m U_blk
10   !!   => all these are used in bulk formulas in sbcblk.F90
11   !!
12   !!    Using the bulk formulation/param. of COARE v3.5, Edson et al. 2013
13   !!
14   !!
15   !!       Routine turb_coare3p5 maintained and developed in AeroBulk
16   !!                     (https://github.com/brodeau/aerobulk/)
17   !!
18   !!            Author: Laurent Brodeau, 2016, brodeau@gmail.com
19   !!
20   !!======================================================================
21   !! History :  3.6  !  2016-02  (L.Brodeau)  Original code
22   !!----------------------------------------------------------------------
23
24   !!----------------------------------------------------------------------
25   !!   turb_coare3p5  : computes the bulk turbulent transfer coefficients
26   !!                   adjusts t_air and q_air from zt to zu m
27   !!                   returns the effective bulk wind speed at 10m
28   !!----------------------------------------------------------------------
29   USE oce             ! ocean dynamics and tracers
30   USE dom_oce         ! ocean space and time domain
31   USE phycst          ! physical constants
32   USE sbc_oce         ! Surface boundary condition: ocean fields
33   USE sbcwave, ONLY   :  cdn_wave ! wave module
34#if defined key_si3 || defined key_cice
35   USE sbc_ice         ! Surface boundary condition: ice fields
36#endif
37   USE sbcblk_phy     !LB: all thermodynamics functions, rho_air, q_sat, etc... #LB
38   USE iom             ! I/O manager library
39   USE lib_mpp         ! distribued memory computing library
40   USE in_out_manager  ! I/O manager
41   USE prtctl          ! Print control
42   USE lib_fortran     ! to use key_nosignedzero
43
44   IMPLICIT NONE
45   PRIVATE
46
47   PUBLIC ::   TURB_COARE3P5   ! called by sbcblk.F90
48
49   !                                   ! COARE own values for given constants:
50   REAL(wp), PARAMETER ::   charn0_max =   0.028   ! value above which the Charnock paramter levels off for winds > 18
51   REAL(wp), PARAMETER ::   zi0        = 600.      ! scale height of the atmospheric boundary layer...1
52   REAL(wp), PARAMETER ::   Beta0      =   1.25    ! gustiness parameter
53
54   !!----------------------------------------------------------------------
55CONTAINS
56
57   SUBROUTINE turb_coare3p5( zt, zu, sst, t_zt, ssq, q_zt, U_zu,  &
58      &                      Cd, Ch, Ce, t_zu, q_zu, U_blk,       &
59      &                      Cdn, Chn, Cen                        )
60      !!----------------------------------------------------------------------------------
61      !!                      ***  ROUTINE  turb_coare3p5  ***
62      !!
63      !! ** Purpose :   Computes turbulent transfert coefficients of surface
64      !!                fluxes according to Fairall et al. (2003)
65      !!                If relevant (zt /= zu), adjust temperature and humidity from height zt to zu
66      !!
67      !! ** Method : Monin Obukhov Similarity Theory
68      !!
69      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
70      !!
71      !! INPUT :
72      !! -------
73      !!    *  zt   : height for temperature and spec. hum. of air            [m]
74      !!    *  zu   : height for wind speed (generally 10m)                   [m]
75      !!    *  U_zu : scalar wind speed at 10m                                [m/s]
76      !!    *  sst  : SST                                                     [K]
77      !!    *  t_zt : potential air temperature at zt                         [K]
78      !!    *  ssq  : specific humidity at saturation at SST                  [kg/kg]
79      !!    *  q_zt : specific humidity of air at zt                          [kg/kg]
80      !!
81      !!
82      !! OUTPUT :
83      !! --------
84      !!    *  Cd     : drag coefficient
85      !!    *  Ch     : sensible heat coefficient
86      !!    *  Ce     : evaporation coefficient
87      !!    *  t_zu   : pot. air temperature adjusted at wind height zu       [K]
88      !!    *  q_zu   : specific humidity of air        //                    [kg/kg]
89      !!    *  U_blk  : bulk wind at 10m                                      [m/s]
90      !!
91      !!----------------------------------------------------------------------------------
92      REAL(wp), INTENT(in   )                     ::   zt       ! height for t_zt and q_zt                   [m]
93      REAL(wp), INTENT(in   )                     ::   zu       ! height for U_zu                              [m]
94      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   sst      ! sea surface temperature              [Kelvin]
95      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   t_zt     ! potential air temperature            [Kelvin]
96      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   ssq      ! sea surface specific humidity         [kg/kg]
97      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_zt     ! specific air humidity                 [kg/kg]
98      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   U_zu     ! relative wind module at zu            [m/s]
99      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cd       ! transfer coefficient for momentum         (tau)
100      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ch       ! transfer coefficient for sensible heat (Q_sens)
101      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ce       ! transfert coefficient for evaporation   (Q_lat)
102      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   t_zu     ! pot. air temp. adjusted at zu             [K]
103      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   q_zu     ! spec. humidity adjusted at zu             [kg/kg]
104      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   U_blk    ! bulk wind at 10m                          [m/s]
105      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients
106      !
107      INTEGER :: j_itt
108      LOGICAL ::   l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U
109      INTEGER , PARAMETER ::   nb_itt = 4       ! number of itterations
110      !
111      REAL(wp), DIMENSION(jpi,jpj) ::  &
112         &  u_star, t_star, q_star, &
113         &  dt_zu, dq_zu,    &
114         &  znu_a,           & !: Nu_air, Viscosity of air
115         &  z0, z0t
116      REAL(wp), DIMENSION(jpi,jpj) ::   zeta_u        ! stability parameter at height zu
117      REAL(wp), DIMENSION(jpi,jpj) ::   ztmp0, ztmp1, ztmp2
118      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zeta_t        ! stability parameter at height zt
119      !!----------------------------------------------------------------------------------
120      !
121      l_zt_equal_zu = .FALSE.
122      IF( ABS(zu - zt) < 0.01 ) l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision
123
124      IF( .NOT. l_zt_equal_zu )   ALLOCATE( zeta_t(jpi,jpj) )
125
126      !! First guess of temperature and humidity at height zu:
127      t_zu = MAX(t_zt , 0.0)    ! who knows what's given on masked-continental regions...
128      q_zu = MAX(q_zt , 1.E-6)  !               "
129
130      !! Pot. temp. difference (and we don't want it to be 0!)
131      dt_zu = t_zu - sst ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6), dt_zu )
132      dq_zu = q_zu - ssq ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9), dq_zu )
133
134      znu_a = visc_air(t_zt) ! Air viscosity (m^2/s) at zt given from temperature in (K)
135
136      ztmp2 = 0.5*0.5  ! initial guess for wind gustiness contribution
137      U_blk = SQRT(U_zu*U_zu + ztmp2)
138
139      ztmp2   = 10000.     ! optimization: ztmp2 == 1/z0 (with z0 first guess == 0.0001)
140      ztmp0   = LOG(zu*ztmp2)
141      ztmp1   = LOG(10.*ztmp2)
142      u_star = 0.035*U_blk*ztmp1/ztmp0       ! (u* = 0.035*Un10)
143
144      !! COARE 3.5 first guess of UN10 is U_zu
145      ztmp2 = MIN( 0.0017*U_zu - 0.005 , charn0_max)  ! alpha Charnock parameter (Eq. 13 Edson al. 2013)
146      ztmp2 = MAX( ztmp2 , 0. )                       ! alpha Charnock parameter (Eq. 13 Edson al. 2013)
147      z0     = ztmp2*u_star*u_star/grav + 0.11*znu_a/u_star
148      z0t    = 0.1*EXP(vkarmn/(0.00115/(vkarmn/ztmp1)))   !  WARNING: 1/z0t !
149
150      ztmp2  = vkarmn/ztmp0
151      Cd     = ztmp2*ztmp2    ! first guess of Cd
152
153      ztmp0 = vkarmn*vkarmn/LOG(zt*z0t)/Cd
154
155      !Ribcu = -zu/(zi0*0.004*Beta0**3) !! Saturation Rib, zi0 = tropicalbound. layer depth
156      ztmp2  = grav*zu*(dt_zu + rctv0*t_zu*dq_zu)/(t_zu*U_blk*U_blk)  !! Ribu Bulk Richardson number
157      ztmp1 = 0.5 + sign(0.5 , ztmp2)
158      ztmp0 = ztmp0*ztmp2
159      !!             Ribu < 0                                 Ribu > 0   Beta = 1.25
160      zeta_u = (1.-ztmp1) * (ztmp0/(1.+ztmp2/(-zu/(zi0*0.004*Beta0**3)))) &
161         &  +     ztmp1   * (ztmp0*(1. + 27./9.*ztmp2/ztmp0))
162
163      !! First guess M-O stability dependent scaling params.(u*,t*,q*) to estimate z0 and z/L
164      ztmp0  =  vkarmn/(LOG(zu*z0t) - psi_h_coare(zeta_u))
165
166      u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0)  - psi_m_coare(zeta_u))
167      t_star = dt_zu*ztmp0
168      q_star = dq_zu*ztmp0
169
170      ! What's need to be done if zt /= zu:
171      IF( .NOT. l_zt_equal_zu ) THEN
172
173         zeta_t = zt*zeta_u/zu
174
175         !! First update of values at zu (or zt for wind)
176         ztmp0 = psi_h_coare(zeta_u) - psi_h_coare(zeta_t)
177         ztmp1 = log(zt/zu) + ztmp0
178         t_zu = t_zt - t_star/vkarmn*ztmp1
179         q_zu = q_zt - q_star/vkarmn*ztmp1
180         q_zu = (0.5 + sign(0.5,q_zu))*q_zu !Makes it impossible to have negative humidity :
181
182         dt_zu = t_zu - sst  ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6), dt_zu )
183         dq_zu = q_zu - ssq  ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9), dq_zu )
184
185      END IF
186
187      !! ITERATION BLOCK
188      DO j_itt = 1, nb_itt
189
190         !!Inverse of Monin-Obukov length (1/L) :
191         ztmp0 = One_on_L(t_zu, q_zu, u_star, t_star, q_star)  ! 1/L == 1/[Monin-Obukhov length]
192
193         ztmp1 = u_star*u_star   ! u*^2
194
195         !! Update wind at 10m taking into acount convection-related wind gustiness:
196         ! Ug = Beta*w*  (Beta = 1.25, Fairall et al. 2003, Eq.8):
197         ztmp2 = Beta0*Beta0*ztmp1*(MAX(-zi0*ztmp0/vkarmn,0.))**(2./3.)   ! => ztmp2 == Ug^2
198         !!   ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before 600.
199         U_blk = MAX(sqrt(U_zu*U_zu + ztmp2), 0.2)        ! include gustiness in bulk wind speed
200         ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0.
201
202         !! COARE 3.5: Charnock parameter is computed from the neutral wind speed at 10m: Eq. 13 (Edson al. 2013)
203         ztmp2 = u_star/vkarmn*LOG(10./z0)   ! UN10 Neutral wind at 10m!
204         ztmp2 = MIN( 0.0017*ztmp2 - 0.005 , charn0_max)  ! alpha Charnock parameter (Eq. 13 Edson al. 2013)
205         ztmp2 = MAX( ztmp2 , 0. )
206
207         !! Roughness lengthes z0, z0t (z0q = z0t) :
208         z0   = ztmp2*ztmp1/grav + 0.11*znu_a/u_star ! Roughness length (eq.6)
209         ztmp1 = z0*u_star/znu_a                             ! Re_r: roughness Reynolds number
210         !z0t  = MIN( 1.1E-4 , 5.5E-5*ztmp1**(-0.6) )   ! COARE 3.0
211         !! Chris Fairall and Jim Edsson, private communication, March 2016 / COARE 3.5 :
212         z0t   = MIN( 1.6e-4 , 5.8E-5*ztmp1**(-0.72)) ! These thermal roughness lengths give Stanton and
213         !z0q = z0t                                   ! Dalton numbers that closely approximate COARE3.0
214
215         !! Stability parameters:
216         zeta_u = zu*ztmp0 ; zeta_u = sign( min(abs(zeta_u),50.0), zeta_u )
217         IF( .NOT. l_zt_equal_zu ) THEN
218            zeta_t = zt*ztmp0 ;  zeta_t = sign( min(abs(zeta_t),50.0), zeta_t )
219         END IF
220
221         !! Turbulent scales at zu=10m :
222         ztmp0   = psi_h_coare(zeta_u)
223         ztmp1   = vkarmn/(LOG(zu) -LOG(z0t) - ztmp0)
224
225         t_star = dt_zu*ztmp1
226         q_star = dq_zu*ztmp1
227         u_star = U_blk*vkarmn/(LOG(zu) -LOG(z0) - psi_m_coare(zeta_u))
228
229         IF( .NOT. l_zt_equal_zu ) THEN
230            ! What's need to be done if zt /= zu
231            !! Re-updating temperature and humidity at zu :
232            ztmp2 = ztmp0 - psi_h_coare(zeta_t)
233            ztmp1 = log(zt/zu) + ztmp2
234            t_zu = t_zt - t_star/vkarmn*ztmp1
235            q_zu = q_zt - q_star/vkarmn*ztmp1
236            dt_zu = t_zu - sst ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6), dt_zu )
237            dq_zu = q_zu - ssq ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9), dq_zu )
238         END IF
239
240      END DO
241      !
242      ! compute transfer coefficients at zu :
243      ztmp0 = u_star/U_blk
244      Cd   = ztmp0*ztmp0
245      Ch   = ztmp0*t_star/dt_zu
246      Ce   = ztmp0*q_star/dq_zu
247      !
248      ztmp1 = zu + z0
249      Cdn = vkarmn*vkarmn / (log(ztmp1/z0 )*log(ztmp1/z0 ))
250      Chn = vkarmn*vkarmn / (log(ztmp1/z0t)*log(ztmp1/z0t))
251      Cen = Chn
252      !
253      IF( .NOT. l_zt_equal_zu ) DEALLOCATE( zeta_t )
254      !
255   END SUBROUTINE turb_coare3p5
256
257   FUNCTION psi_m_coare( pzeta )
258      !!----------------------------------------------------------------------------------
259      !! ** Purpose: compute the universal profile stability function for momentum
260      !!             COARE 3.0, Fairall et al. 2003
261      !!             pzeta : stability paramenter, z/L where z is altitude
262      !!                     measurement and L is M-O length
263      !!       Stability function for wind speed and scalars matching Kansas and free
264      !!       convection forms with weighting f convective form, follows Fairall et
265      !!       al (1996) with profile constants from Grachev et al (2000) BLM stable
266      !!       form from Beljaars and Holtslag (1991)
267      !!
268      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
269      !!----------------------------------------------------------------------------------
270      REAL(wp), DIMENSION(jpi,jpj) :: psi_m_coare
271      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
272      !
273      INTEGER  ::   ji, jj    ! dummy loop indices
274      REAL(wp) :: zta, zphi_m, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
275      !!----------------------------------------------------------------------------------
276      !
277      DO jj = 1, jpj
278         DO ji = 1, jpi
279            !
280            zta = pzeta(ji,jj)
281            !
282            zphi_m = ABS(1. - 15.*zta)**.25    !!Kansas unstable
283            !
284            zpsi_k = 2.*LOG((1. + zphi_m)/2.) + LOG((1. + zphi_m*zphi_m)/2.)   &
285               & - 2.*ATAN(zphi_m) + 0.5*rpi
286            !
287            zphi_c = ABS(1. - 10.15*zta)**.3333                   !!Convective
288            !
289            zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
290               &     - 1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
291            !
292            zf = zta*zta
293            zf = zf/(1. + zf)
294            zc = MIN(50., 0.35*zta)
295            zstab = 0.5 + SIGN(0.5, zta)
296            !
297            psi_m_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & ! (zta < 0)
298               &                -   zstab     * ( 1. + 1.*zta     &                ! (zta > 0)
299               &                         + 0.6667*(zta - 14.28)/EXP(zc) + 8.525 )   !     "
300            !
301         END DO
302      END DO
303      !
304   END FUNCTION psi_m_coare
305
306
307   FUNCTION psi_h_coare( pzeta )
308      !!---------------------------------------------------------------------
309      !! Universal profile stability function for temperature and humidity
310      !! COARE 3.0, Fairall et al. 2003
311      !!
312      !! pzeta : stability paramenter, z/L where z is altitude measurement
313      !!         and L is M-O length
314      !!
315      !! Stability function for wind speed and scalars matching Kansas and free
316      !! convection forms with weighting f convective form, follows Fairall et
317      !! al (1996) with profile constants from Grachev et al (2000) BLM stable
318      !! form from Beljaars and Holtslag (1991)
319      !!
320      !! Author: L. Brodeau, june 2016 / AeroBulk
321      !!         (https://github.com/brodeau/aerobulk/)
322      !!----------------------------------------------------------------
323      !!
324      REAL(wp), DIMENSION(jpi,jpj) :: psi_h_coare
325      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
326      !
327      INTEGER  ::   ji, jj     ! dummy loop indices
328      REAL(wp) :: zta, zphi_h, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
329      !
330      DO jj = 1, jpj
331         DO ji = 1, jpi
332            !
333            zta = pzeta(ji,jj)
334            !
335            zphi_h = (ABS(1. - 15.*zta))**.5  !! Kansas unstable   (zphi_h = zphi_m**2 when unstable, zphi_m when stable)
336            !
337            zpsi_k = 2.*LOG((1. + zphi_h)/2.)
338            !
339            zphi_c = (ABS(1. - 34.15*zta))**.3333   !! Convective
340            !
341            zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
342               &    -1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
343            !
344            zf = zta*zta
345            zf = zf/(1. + zf)
346            zc = MIN(50.,0.35*zta)
347            zstab = 0.5 + SIGN(0.5, zta)
348            !
349            psi_h_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) &
350               &                -   zstab     * ( (ABS(1. + 2.*zta/3.))**1.5     &
351               &                           + .6667*(zta - 14.28)/EXP(zc) + 8.525 )
352            !
353         END DO
354      END DO
355      !
356   END FUNCTION psi_h_coare
357
358   !!======================================================================
359END MODULE sbcblk_algo_coare3p5
Note: See TracBrowser for help on using the repository browser.