source: NEMO/trunk/src/OCE/SBC/sbcblk_algo_coare3p5.F90 @ 10522

Last change on this file since 10522 was 10069, checked in by nicolasmartin, 2 years ago

Fix mistakes of previous commit on SVN keywords property

  • Property svn:keywords set to Id
File size: 19.4 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   !!                     (http://aerobulk.sourceforge.net/)
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   !
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   REAL(wp), PARAMETER ::   rctv0      =   0.608   ! constant to obtain virtual temperature...
54
55   !!----------------------------------------------------------------------
56CONTAINS
57
58   SUBROUTINE turb_coare3p5( zt, zu, sst, t_zt, ssq, q_zt, U_zu,  &
59      &                      Cd, Ch, Ce, t_zu, q_zu, U_blk,       &
60      &                      Cdn, Chn, Cen                        )
61      !!----------------------------------------------------------------------------------
62      !!                      ***  ROUTINE  turb_coare3p5  ***
63      !!
64      !! ** Purpose :   Computes turbulent transfert coefficients of surface
65      !!                fluxes according to Fairall et al. (2003)
66      !!                If relevant (zt /= zu), adjust temperature and humidity from height zt to zu
67      !!
68      !! ** Method : Monin Obukhov Similarity Theory
69      !!
70      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
71      !!
72      !! INPUT :
73      !! -------
74      !!    *  zt   : height for temperature and spec. hum. of air            [m]
75      !!    *  zu   : height for wind speed (generally 10m)                   [m]
76      !!    *  U_zu : scalar wind speed at 10m                                [m/s]
77      !!    *  sst  : SST                                                     [K]
78      !!    *  t_zt : potential air temperature at zt                         [K]
79      !!    *  ssq  : specific humidity at saturation at SST                  [kg/kg]
80      !!    *  q_zt : specific humidity of air at zt                          [kg/kg]
81      !!
82      !!
83      !! OUTPUT :
84      !! --------
85      !!    *  Cd     : drag coefficient
86      !!    *  Ch     : sensible heat coefficient
87      !!    *  Ce     : evaporation coefficient
88      !!    *  t_zu   : pot. air temperature adjusted at wind height zu       [K]
89      !!    *  q_zu   : specific humidity of air        //                    [kg/kg]
90      !!    *  U_blk  : bulk wind at 10m                                      [m/s]
91      !!
92      !!----------------------------------------------------------------------------------
93      REAL(wp), INTENT(in   )                     ::   zt       ! height for t_zt and q_zt                   [m]
94      REAL(wp), INTENT(in   )                     ::   zu       ! height for U_zu                              [m]
95      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   sst      ! sea surface temperature              [Kelvin]
96      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   t_zt     ! potential air temperature            [Kelvin]
97      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   ssq      ! sea surface specific humidity         [kg/kg]
98      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_zt     ! specific air humidity                 [kg/kg]
99      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   U_zu     ! relative wind module at zu            [m/s]
100      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cd       ! transfer coefficient for momentum         (tau)
101      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ch       ! transfer coefficient for sensible heat (Q_sens)
102      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ce       ! transfert coefficient for evaporation   (Q_lat)
103      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   t_zu     ! pot. air temp. adjusted at zu             [K]
104      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   q_zu     ! spec. humidity adjusted at zu             [kg/kg]
105      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   U_blk    ! bulk wind at 10m                          [m/s]
106      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients
107      !
108      INTEGER :: j_itt
109      LOGICAL ::   l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U
110      INTEGER , PARAMETER ::   nb_itt = 4       ! number of itterations
111      !
112      REAL(wp), DIMENSION(jpi,jpj) ::  &
113         &  u_star, t_star, q_star, &
114         &  dt_zu, dq_zu,    &
115         &  znu_a,           & !: Nu_air, Viscosity of air
116         &  z0, z0t
117      REAL(wp), DIMENSION(jpi,jpj) ::   zeta_u        ! stability parameter at height zu
118      REAL(wp), DIMENSION(jpi,jpj) ::   ztmp0, ztmp1, ztmp2
119      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zeta_t        ! stability parameter at height zt
120      !!----------------------------------------------------------------------------------
121      !
122      l_zt_equal_zu = .FALSE.
123      IF( ABS(zu - zt) < 0.01 ) l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision
124
125      IF( .NOT. l_zt_equal_zu )   ALLOCATE( zeta_t(jpi,jpj) )
126
127      !! First guess of temperature and humidity at height zu:
128      t_zu = MAX(t_zt , 0.0)    ! who knows what's given on masked-continental regions...
129      q_zu = MAX(q_zt , 1.E-6)  !               "
130
131      !! Pot. temp. difference (and we don't want it to be 0!)
132      dt_zu = t_zu - sst ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6), dt_zu )
133      dq_zu = q_zu - ssq ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9), dq_zu )
134
135      znu_a = visc_air(t_zt) ! Air viscosity (m^2/s) at zt given from temperature in (K)
136
137      ztmp2 = 0.5*0.5  ! initial guess for wind gustiness contribution
138      U_blk = SQRT(U_zu*U_zu + ztmp2)
139
140      ztmp2   = 10000.     ! optimization: ztmp2 == 1/z0 (with z0 first guess == 0.0001)
141      ztmp0   = LOG(zu*ztmp2)
142      ztmp1   = LOG(10.*ztmp2)
143      u_star = 0.035*U_blk*ztmp1/ztmp0       ! (u* = 0.035*Un10)
144
145      !! COARE 3.5 first guess of UN10 is U_zu
146      ztmp2 = MIN( 0.0017*U_zu - 0.005 , charn0_max)  ! alpha Charnock parameter (Eq. 13 Edson al. 2013)
147      ztmp2 = MAX( ztmp2 , 0. )                       ! alpha Charnock parameter (Eq. 13 Edson al. 2013)
148      z0     = ztmp2*u_star*u_star/grav + 0.11*znu_a/u_star
149      z0t    = 0.1*EXP(vkarmn/(0.00115/(vkarmn/ztmp1)))   !  WARNING: 1/z0t !
150
151      ztmp2  = vkarmn/ztmp0
152      Cd     = ztmp2*ztmp2    ! first guess of Cd
153
154      ztmp0 = vkarmn*vkarmn/LOG(zt*z0t)/Cd
155
156      !Ribcu = -zu/(zi0*0.004*Beta0**3) !! Saturation Rib, zi0 = tropicalbound. layer depth
157      ztmp2  = grav*zu*(dt_zu + rctv0*t_zu*dq_zu)/(t_zu*U_blk*U_blk)  !! Ribu Bulk Richardson number
158      ztmp1 = 0.5 + sign(0.5 , ztmp2)
159      ztmp0 = ztmp0*ztmp2
160      !!             Ribu < 0                                 Ribu > 0   Beta = 1.25
161      zeta_u = (1.-ztmp1) * (ztmp0/(1.+ztmp2/(-zu/(zi0*0.004*Beta0**3)))) &
162         &  +     ztmp1   * (ztmp0*(1. + 27./9.*ztmp2/ztmp0))
163
164      !! First guess M-O stability dependent scaling params.(u*,t*,q*) to estimate z0 and z/L
165      ztmp0  =  vkarmn/(LOG(zu*z0t) - psi_h_coare(zeta_u))
166
167      u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0)  - psi_m_coare(zeta_u))
168      t_star = dt_zu*ztmp0
169      q_star = dq_zu*ztmp0
170
171      ! What's need to be done if zt /= zu:
172      IF( .NOT. l_zt_equal_zu ) THEN
173
174         zeta_t = zt*zeta_u/zu
175
176         !! First update of values at zu (or zt for wind)
177         ztmp0 = psi_h_coare(zeta_u) - psi_h_coare(zeta_t)
178         ztmp1 = log(zt/zu) + ztmp0
179         t_zu = t_zt - t_star/vkarmn*ztmp1
180         q_zu = q_zt - q_star/vkarmn*ztmp1
181         q_zu = (0.5 + sign(0.5,q_zu))*q_zu !Makes it impossible to have negative humidity :
182
183         dt_zu = t_zu - sst  ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6), dt_zu )
184         dq_zu = q_zu - ssq  ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9), dq_zu )
185
186      END IF
187
188      !! ITERATION BLOCK
189      DO j_itt = 1, nb_itt
190
191         !!Inverse of Monin-Obukov length (1/L) :
192         ztmp0 = One_on_L(t_zu, q_zu, u_star, t_star, q_star)  ! 1/L == 1/[Monin-Obukhov length]
193
194         ztmp1 = u_star*u_star   ! u*^2
195
196         !! Update wind at 10m taking into acount convection-related wind gustiness:
197         ! Ug = Beta*w*  (Beta = 1.25, Fairall et al. 2003, Eq.8):
198         ztmp2 = Beta0*Beta0*ztmp1*(MAX(-zi0*ztmp0/vkarmn,0.))**(2./3.)   ! => ztmp2 == Ug^2
199         !!   ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before 600.
200         U_blk = MAX(sqrt(U_zu*U_zu + ztmp2), 0.2)        ! include gustiness in bulk wind speed
201         ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0.
202
203         !! COARE 3.5: Charnock parameter is computed from the neutral wind speed at 10m: Eq. 13 (Edson al. 2013)
204         ztmp2 = u_star/vkarmn*LOG(10./z0)   ! UN10 Neutral wind at 10m!
205         ztmp2 = MIN( 0.0017*ztmp2 - 0.005 , charn0_max)  ! alpha Charnock parameter (Eq. 13 Edson al. 2013)
206         ztmp2 = MAX( ztmp2 , 0. )
207
208         !! Roughness lengthes z0, z0t (z0q = z0t) :
209         z0   = ztmp2*ztmp1/grav + 0.11*znu_a/u_star ! Roughness length (eq.6)
210         ztmp1 = z0*u_star/znu_a                             ! Re_r: roughness Reynolds number
211         !z0t  = MIN( 1.1E-4 , 5.5E-5*ztmp1**(-0.6) )   ! COARE 3.0
212         !! Chris Fairall and Jim Edsson, private communication, March 2016 / COARE 3.5 :
213         z0t   = MIN( 1.6e-4 , 5.8E-5*ztmp1**(-0.72)) ! These thermal roughness lengths give Stanton and
214         !z0q = z0t                                   ! Dalton numbers that closely approximate COARE3.0
215
216         !! Stability parameters:
217         zeta_u = zu*ztmp0 ; zeta_u = sign( min(abs(zeta_u),50.0), zeta_u )
218         IF( .NOT. l_zt_equal_zu ) THEN
219            zeta_t = zt*ztmp0 ;  zeta_t = sign( min(abs(zeta_t),50.0), zeta_t )
220         END IF
221
222         !! Turbulent scales at zu=10m :
223         ztmp0   = psi_h_coare(zeta_u)
224         ztmp1   = vkarmn/(LOG(zu) -LOG(z0t) - ztmp0)
225
226         t_star = dt_zu*ztmp1
227         q_star = dq_zu*ztmp1
228         u_star = U_blk*vkarmn/(LOG(zu) -LOG(z0) - psi_m_coare(zeta_u))
229
230         IF( .NOT. l_zt_equal_zu ) THEN
231            ! What's need to be done if zt /= zu
232            !! Re-updating temperature and humidity at zu :
233            ztmp2 = ztmp0 - psi_h_coare(zeta_t)
234            ztmp1 = log(zt/zu) + ztmp2
235            t_zu = t_zt - t_star/vkarmn*ztmp1
236            q_zu = q_zt - q_star/vkarmn*ztmp1
237            dt_zu = t_zu - sst ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6), dt_zu )
238            dq_zu = q_zu - ssq ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9), dq_zu )
239         END IF
240
241      END DO
242      !
243      ! compute transfer coefficients at zu :
244      ztmp0 = u_star/U_blk
245      Cd   = ztmp0*ztmp0
246      Ch   = ztmp0*t_star/dt_zu
247      Ce   = ztmp0*q_star/dq_zu
248      !
249      ztmp1 = zu + z0
250      Cdn = vkarmn*vkarmn / (log(ztmp1/z0 )*log(ztmp1/z0 ))
251      Chn = vkarmn*vkarmn / (log(ztmp1/z0t)*log(ztmp1/z0t))
252      Cen = Chn
253      !
254      IF( .NOT. l_zt_equal_zu ) DEALLOCATE( zeta_t )
255      !
256   END SUBROUTINE turb_coare3p5
257
258
259
260   FUNCTION One_on_L( ptha, pqa, pus, pts, pqs )
261      !!------------------------------------------------------------------------
262      !!
263      !! Evaluates the 1./(Monin Obukhov length) from air temperature and
264      !!  specific humidity, and frictional scales u*, t* and q*
265      !!
266      !! Author: L. Brodeau, june 2016 / AeroBulk
267      !!         (https://sourceforge.net/p/aerobulk)
268      !!------------------------------------------------------------------------
269      REAL(wp), DIMENSION(jpi,jpj)             :: One_on_L         !: 1./(Monin Obukhov length) [m^-1]
270      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptha,  &  !: average potetntial air temperature [K]
271         &                                        pqa,   &  !: average specific humidity of air   [kg/kg]
272         &                                      pus, pts, pqs   !: frictional velocity, temperature and humidity
273      !
274      INTEGER  ::   ji, jj         ! dummy loop indices
275      REAL(wp) ::     zqa          ! local scalar
276      !
277      DO jj = 1, jpj
278         DO ji = 1, jpi
279            !
280            zqa = (1. + rctv0*pqa(ji,jj))
281            !
282            One_on_L(ji,jj) =  grav*vkarmn*(pts(ji,jj)*zqa + rctv0*ptha(ji,jj)*pqs(ji,jj)) &
283               &                      / ( pus(ji,jj)*pus(ji,jj) * ptha(ji,jj)*zqa )
284            !
285         END DO
286      END DO
287      !
288   END FUNCTION One_on_L
289
290
291   FUNCTION psi_m_coare( pzeta )
292      !!----------------------------------------------------------------------------------
293      !! ** Purpose: compute the universal profile stability function for momentum
294      !!             COARE 3.0, Fairall et al. 2003
295      !!             pzeta : stability paramenter, z/L where z is altitude
296      !!                     measurement and L is M-O length
297      !!       Stability function for wind speed and scalars matching Kansas and free
298      !!       convection forms with weighting f convective form, follows Fairall et
299      !!       al (1996) with profile constants from Grachev et al (2000) BLM stable
300      !!       form from Beljaars and Holtslag (1991)
301      !!
302      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
303      !!----------------------------------------------------------------------------------
304      REAL(wp), DIMENSION(jpi,jpj) :: psi_m_coare
305      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
306      !
307      INTEGER  ::   ji, jj    ! dummy loop indices
308      REAL(wp) :: zta, zphi_m, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
309      !!----------------------------------------------------------------------------------
310      !
311      DO jj = 1, jpj
312         DO ji = 1, jpi
313            !
314            zta = pzeta(ji,jj)
315            !
316            zphi_m = ABS(1. - 15.*zta)**.25    !!Kansas unstable
317            !
318            zpsi_k = 2.*LOG((1. + zphi_m)/2.) + LOG((1. + zphi_m*zphi_m)/2.)   &
319               & - 2.*ATAN(zphi_m) + 0.5*rpi
320            !
321            zphi_c = ABS(1. - 10.15*zta)**.3333                   !!Convective
322            !
323            zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
324               &     - 1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
325            !
326            zf = zta*zta
327            zf = zf/(1. + zf)
328            zc = MIN(50., 0.35*zta)
329            zstab = 0.5 + SIGN(0.5, zta)
330            !
331            psi_m_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & ! (zta < 0)
332               &                -   zstab     * ( 1. + 1.*zta     &                ! (zta > 0)
333               &                         + 0.6667*(zta - 14.28)/EXP(zc) + 8.525 )   !     "
334            !
335         END DO
336      END DO
337      !
338   END FUNCTION psi_m_coare
339
340
341   FUNCTION psi_h_coare( pzeta )
342      !!---------------------------------------------------------------------
343      !! Universal profile stability function for temperature and humidity
344      !! COARE 3.0, Fairall et al. 2003
345      !!
346      !! pzeta : stability paramenter, z/L where z is altitude measurement
347      !!         and L is M-O length
348      !!
349      !! Stability function for wind speed and scalars matching Kansas and free
350      !! convection forms with weighting f convective form, follows Fairall et
351      !! al (1996) with profile constants from Grachev et al (2000) BLM stable
352      !! form from Beljaars and Holtslag (1991)
353      !!
354      !! Author: L. Brodeau, june 2016 / AeroBulk
355      !!         (https://sourceforge.net/p/aerobulk)
356      !!----------------------------------------------------------------
357      !!
358      REAL(wp), DIMENSION(jpi,jpj) :: psi_h_coare
359      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
360      !
361      INTEGER  ::   ji, jj     ! dummy loop indices
362      REAL(wp) :: zta, zphi_h, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
363      !
364      DO jj = 1, jpj
365         DO ji = 1, jpi
366            !
367            zta = pzeta(ji,jj)
368            !
369            zphi_h = (ABS(1. - 15.*zta))**.5  !! Kansas unstable   (zphi_h = zphi_m**2 when unstable, zphi_m when stable)
370            !
371            zpsi_k = 2.*LOG((1. + zphi_h)/2.)
372            !
373            zphi_c = (ABS(1. - 34.15*zta))**.3333   !! Convective
374            !
375            zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
376               &    -1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
377            !
378            zf = zta*zta
379            zf = zf/(1. + zf)
380            zc = MIN(50.,0.35*zta)
381            zstab = 0.5 + SIGN(0.5, zta)
382            !
383            psi_h_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) &
384               &                -   zstab     * ( (ABS(1. + 2.*zta/3.))**1.5     &
385               &                           + .6667*(zta - 14.28)/EXP(zc) + 8.525 )
386            !
387         END DO
388      END DO
389      !
390   END FUNCTION psi_h_coare
391
392
393   FUNCTION visc_air( ptak )
394      !!---------------------------------------------------------------------
395      !! Air kinetic viscosity (m^2/s) given from temperature in degrees...
396      !!
397      !! Author: L. Brodeau, june 2016 / AeroBulk
398      !!         (https://sourceforge.net/p/aerobulk)
399      !!---------------------------------------------------------------------
400      REAL(wp), DIMENSION(jpi,jpj)             ::   visc_air
401      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak       ! air temperature   [K]
402      !
403      INTEGER  ::   ji, jj         ! dummy loop indices
404      REAL(wp) ::   ztc, ztc2      ! local scalar
405      !
406      DO jj = 1, jpj
407         DO ji = 1, jpi
408            ztc  = ptak(ji,jj) - rt0   ! air temp, in deg. C
409            ztc2 = ztc*ztc
410            visc_air(ji,jj) = 1.326E-5*(1. + 6.542E-3*ztc + 8.301E-6*ztc2 - 4.84E-9*ztc2*ztc)
411         END DO
412      END DO
413      !
414   END FUNCTION visc_air
415
416   !!======================================================================
417END MODULE sbcblk_algo_coare3p5
Note: See TracBrowser for help on using the repository browser.