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 branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_coare.F90 @ 9019

Last change on this file since 9019 was 9019, checked in by timgraham, 6 years ago

Merge of dev_CNRS_2017 into branch
File size: 20.6 KB
Line 
1MODULE sbcblk_algo_coare
2   !!======================================================================
3   !!                       ***  MODULE  sbcblk_algo_coare  ***
4   !! Computes turbulent components of surface fluxes
5   !!         according to Fairall et al. 2003 (COARE v3)
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, Fairall et al. 2003
13   !!
14   !!
15   !!       Routine turb_coare 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_coare  : 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_lim3 || defined key_cice
35   USE sbc_ice        ! Surface boundary condition: ice fields
36#endif
37   !
38   USE in_out_manager ! I/O manager
39   USE iom            ! I/O manager library
40   USE lib_mpp        ! distribued memory computing library
41   USE wrk_nemo       ! work arrays
42   USE timing         ! Timing
43   USE prtctl         ! Print control
44   USE lib_fortran    ! to use key_nosignedzero
45
46   IMPLICIT NONE
47   PRIVATE
48
49   PUBLIC ::   TURB_COARE   ! called by sbcblk.F90
50
51   !                                              !! COARE own values for given constants:
52   REAL(wp), PARAMETER ::   zi0     = 600._wp      ! scale height of the atmospheric boundary layer...
53   REAL(wp), PARAMETER ::   Beta0   =   1.250_wp   ! gustiness parameter
54   REAL(wp), PARAMETER ::   rctv0   =   0.608_wp   ! constant to obtain virtual temperature...
55
56   !!----------------------------------------------------------------------
57CONTAINS
58
59   SUBROUTINE turb_coare( zt, zu, sst, t_zt, ssq, q_zt, U_zu, &
60      &                   Cd, Ch, Ce, t_zu, q_zu, U_blk,      &
61      &                   Cdn, Chn, Cen                       )
62      !!----------------------------------------------------------------------
63      !!                      ***  ROUTINE  turb_coare  ***
64      !!
65      !!            2015: L. Brodeau (brodeau@gmail.com)
66      !!
67      !! ** Purpose :   Computes turbulent transfert coefficients of surface
68      !!                fluxes according to Fairall et al. (2003)
69      !!                If relevant (zt /= zu), adjust temperature and humidity from height zt to zu
70      !!
71      !! ** Method : Monin Obukhov Similarity Theory
72      !!----------------------------------------------------------------------
73      !!
74      !! INPUT :
75      !! -------
76      !!    *  zt   : height for temperature and spec. hum. of air            [m]
77      !!    *  zu   : height for wind speed (generally 10m)                   [m]
78      !!    *  U_zu : scalar wind speed at 10m                                [m/s]
79      !!    *  sst  : SST                                                     [K]
80      !!    *  t_zt : potential air temperature at zt                         [K]
81      !!    *  ssq  : specific humidity at saturation at SST                  [kg/kg]
82      !!    *  q_zt : specific humidity of air at zt                          [kg/kg]
83      !!
84      !!
85      !! OUTPUT :
86      !! --------
87      !!    *  Cd     : drag coefficient
88      !!    *  Ch     : sensible heat coefficient
89      !!    *  Ce     : evaporation coefficient
90      !!    *  t_zu   : pot. air temperature adjusted at wind height zu       [K]
91      !!    *  q_zu   : specific humidity of air        //                    [kg/kg]
92      !!    *  U_blk  : bulk wind at 10m                                      [m/s]
93      !!----------------------------------------------------------------------
94      REAL(wp), INTENT(in   )                     ::   zt       ! height for t_zt and q_zt                    [m]
95      REAL(wp), INTENT(in   )                     ::   zu       ! height for U_zu                             [m]
96      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   sst      ! sea surface temperature                [Kelvin]
97      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   t_zt     ! potential air temperature              [Kelvin]
98      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   ssq      ! sea surface specific humidity           [kg/kg]
99      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_zt     ! specific air humidity at zt             [kg/kg]
100      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   U_zu     ! relative wind module at zu                [m/s]
101      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cd       ! transfer coefficient for momentum         (tau)
102      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ch       ! transfer coefficient for sensible heat (Q_sens)
103      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ce       ! transfert coefficient for evaporation   (Q_lat)
104      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   t_zu     ! pot. air temp. adjusted at zu               [K]
105      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   q_zu     ! spec. humidity adjusted at zu           [kg/kg]
106      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   U_blk    ! bulk wind at 10m                          [m/s]
107      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients
108      !
109      INTEGER :: j_itt
110      LOGICAL ::   l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U
111      INTEGER , PARAMETER ::   nb_itt = 4       ! number of itterations
112
113      REAL(wp), DIMENSION(:,:), POINTER  ::  &
114         &  u_star, t_star, q_star, &
115         &  dt_zu, dq_zu,    &
116         &  znu_a,           & !: Nu_air, Viscosity of air
117         &  z0, z0t
118      REAL(wp), DIMENSION(:,:), POINTER ::   zeta_u        ! stability parameter at height zu
119      REAL(wp), DIMENSION(:,:), POINTER ::   zeta_t        ! stability parameter at height zt
120      REAL(wp), DIMENSION(:,:), POINTER ::   ztmp0, ztmp1, ztmp2
121      !!----------------------------------------------------------------------
122      !
123      IF( nn_timing == 1 )  CALL timing_start('turb_coare')
124
125      CALL wrk_alloc( jpi,jpj,   u_star, t_star, q_star, zeta_u, dt_zu, dq_zu)
126      CALL wrk_alloc( jpi,jpj,   znu_a, z0, z0t, ztmp0, ztmp1, ztmp2 )
127
128      l_zt_equal_zu = .FALSE.
129      IF( ABS(zu - zt) < 0.01 ) l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision
130
131      IF( .NOT. l_zt_equal_zu )   CALL wrk_alloc( jpi,jpj, zeta_t )
132
133      !! First guess of temperature and humidity at height zu:
134      t_zu = MAX(t_zt , 0.0)    ! who knows what's given on masked-continental regions...
135      q_zu = MAX(q_zt , 1.e-6)  !               "
136
137      !! Pot. temp. difference (and we don't want it to be 0!)
138      dt_zu = t_zu - sst ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6), dt_zu )
139      dq_zu = q_zu - ssq ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9), dq_zu )
140
141      znu_a = visc_air(t_zt) ! Air viscosity (m^2/s) at zt given from temperature in (K)
142
143      ztmp2 = 0.5*0.5  ! initial guess for wind gustiness contribution
144      U_blk = SQRT(U_zu*U_zu + ztmp2)
145
146      ztmp2   = 10000.     ! optimization: ztmp2 == 1/z0 (with z0 first guess == 0.0001)
147      ztmp0   = LOG(zu*ztmp2)
148      ztmp1   = LOG(10.*ztmp2)
149      u_star = 0.035*U_blk*ztmp1/ztmp0       ! (u* = 0.035*Un10)
150
151
152      z0     = alfa_charn(U_blk)*u_star*u_star/grav + 0.11*znu_a/u_star
153      z0t    = 0.1*EXP(vkarmn/(0.00115/(vkarmn/ztmp1)))   !  WARNING: 1/z0t !
154
155      ztmp2  = vkarmn/ztmp0
156      Cd     = ztmp2*ztmp2    ! first guess of Cd
157
158      ztmp0 = vkarmn*vkarmn/LOG(zt*z0t)/Cd
159
160      !Ribcu = -zu/(zi0*0.004*Beta0**3) !! Saturation Rib, zi0 = tropicalbound. layer depth
161      ztmp2  = grav*zu*(dt_zu + rctv0*t_zu*dq_zu)/(t_zu*U_blk*U_blk)  !! Ribu Bulk Richardson number
162      ztmp1 = 0.5 + sign(0.5 , ztmp2)
163      ztmp0 = ztmp0*ztmp2
164      !!             Ribu < 0                                 Ribu > 0   Beta = 1.25
165      zeta_u = (1.-ztmp1) * (ztmp0/(1.+ztmp2/(-zu/(zi0*0.004*Beta0**3)))) &
166         &  +     ztmp1   * (ztmp0*(1. + 27./9.*ztmp2/ztmp0))
167
168      !! First guess M-O stability dependent scaling params.(u*,t*,q*) to estimate z0 and z/L
169      ztmp0  =  vkarmn/(LOG(zu*z0t) - psi_h_coare(zeta_u))
170
171      u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0)  - psi_m_coare(zeta_u))
172      t_star = dt_zu*ztmp0
173      q_star = dq_zu*ztmp0
174
175      ! What's need to be done if zt /= zu:
176      IF( .NOT. l_zt_equal_zu ) THEN
177
178         zeta_t = zt*zeta_u/zu
179
180         !! First update of values at zu (or zt for wind)
181         ztmp0 = psi_h_coare(zeta_u) - psi_h_coare(zeta_t)
182         ztmp1 = log(zt/zu) + ztmp0
183         t_zu = t_zt - t_star/vkarmn*ztmp1
184         q_zu = q_zt - q_star/vkarmn*ztmp1
185         q_zu = (0.5 + sign(0.5,q_zu))*q_zu !Makes it impossible to have negative humidity :
186
187         dt_zu = t_zu - sst  ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6), dt_zu )
188         dq_zu = q_zu - ssq  ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9), dq_zu )
189
190      END IF
191
192      !! ITERATION BLOCK
193      DO j_itt = 1, nb_itt
194
195         !!Inverse of Monin-Obukov length (1/L) :
196         ztmp0 = One_on_L(t_zu, q_zu, u_star, t_star, q_star)  ! 1/L == 1/[Monin-Obukhov length]
197
198         ztmp1 = u_star*u_star   ! u*^2
199
200         !! Update wind at 10m taking into acount convection-related wind gustiness:
201         ! Ug = Beta*w*  (Beta = 1.25, Fairall et al. 2003, Eq.8):
202         ztmp2 = Beta0*Beta0*ztmp1*(MAX(-zi0*ztmp0/vkarmn,0.))**(2./3.)   ! => ztmp2 == Ug^2
203         !!   ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before 600.
204         U_blk = MAX(sqrt(U_zu*U_zu + ztmp2), 0.2)        ! include gustiness in bulk wind speed
205         ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0.
206
207         !! Updating Charnock parameter, increases with the wind (Fairall et al., 2003 p. 577-578)
208         ztmp2 = alfa_charn(U_blk)  ! alpha Charnock parameter
209
210         !! Roughness lengthes z0, z0t (z0q = z0t) :
211         z0   = ztmp2*ztmp1/grav + 0.11*znu_a/u_star ! Roughness length (eq.6)
212         ztmp1 = z0*u_star/znu_a                             ! Re_r: roughness Reynolds number
213         z0t  = min( 1.1E-4 , 5.5E-5*ztmp1**(-0.6) )         ! Scalar roughness for both theta and q (eq.28)
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      CALL wrk_dealloc( jpi,jpj, u_star, t_star, q_star, zeta_u, dt_zu, dq_zu )
254      CALL wrk_dealloc( jpi,jpj, znu_a, z0, z0t, ztmp0, ztmp1, ztmp2 )
255      IF( .NOT. l_zt_equal_zu ) CALL wrk_dealloc( jpi,jpj, zeta_t )
256
257      IF( nn_timing == 1 )  CALL timing_stop('turb_coare')
258
259   END SUBROUTINE turb_coare
260
261
262   FUNCTION alfa_charn( pwnd )
263      !!-------------------------------------------------------------------
264      !! Compute the Charnock parameter as a function of the wind speed
265      !!
266      !! (Fairall et al., 2003 p.577-578)
267      !!
268      !! Wind below 10 m/s :  alfa = 0.011
269      !! Wind between 10 and 18 m/s : linear increase from 0.011 to 0.018
270      !! Wind greater than 18 m/s :  alfa = 0.018
271      !!
272      !! Author: L. Brodeau, june 2016 / AeroBulk  (https://sourceforge.net/p/aerobulk)
273      !!-------------------------------------------------------------------
274      REAL(wp), DIMENSION(jpi,jpj) :: alfa_charn
275      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pwnd   ! wind speed
276      !
277      INTEGER  ::   ji, jj         ! dummy loop indices
278      REAL(wp) :: zw, zgt10, zgt18
279      !!-------------------------------------------------------------------
280      !
281      DO jj = 1, jpj
282         DO ji = 1, jpi
283            !
284            zw = pwnd(ji,jj)   ! wind speed
285            !
286            ! Charnock's constant, increases with the wind :
287            zgt10 = 0.5 + SIGN(0.5,(zw - 10.)) ! If zw<10. --> 0, else --> 1
288            zgt18 = 0.5 + SIGN(0.5,(zw - 18.)) ! If zw<18. --> 0, else --> 1
289            !
290            alfa_charn(ji,jj) =  (1. - zgt10)*0.011    &    ! wind is lower than 10 m/s
291               &     + zgt10*((1. - zgt18)*(0.011 + (0.018 - 0.011) &
292               &      *(zw - 10.)/(18. - 10.)) + zgt18*( 0.018 ) )    ! Hare et al. (1999)
293            !
294         END DO
295      END DO
296      !
297   END FUNCTION alfa_charn
298
299
300   FUNCTION One_on_L( ptha, pqa, pus, pts, pqs )
301      !!------------------------------------------------------------------------
302      !!
303      !! Evaluates the 1./(Monin Obukhov length) from air temperature and
304      !!  specific humidity, and frictional scales u*, t* and q*
305      !!
306      !! Author: L. Brodeau, june 2016 / AeroBulk
307      !!         (https://sourceforge.net/p/aerobulk)
308      !!------------------------------------------------------------------------
309      REAL(wp), DIMENSION(jpi,jpj)             :: One_on_L         !: 1./(Monin Obukhov length) [m^-1]
310      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptha,  &  !: average potetntial air temperature [K]
311         &                                        pqa,   &  !: average specific humidity of air   [kg/kg]
312         &                                      pus, pts, pqs   !: frictional velocity, temperature and humidity
313      !
314      INTEGER  ::   ji, jj         ! dummy loop indices
315      REAL(wp) ::     zqa          ! local scalar
316      !!-------------------------------------------------------------------
317      !
318      DO jj = 1, jpj
319         DO ji = 1, jpi
320            !
321            zqa = (1. + rctv0*pqa(ji,jj))
322            !
323            One_on_L(ji,jj) =  grav*vkarmn*(pts(ji,jj)*zqa + rctv0*ptha(ji,jj)*pqs(ji,jj)) &
324               &                      / ( pus(ji,jj)*pus(ji,jj) * ptha(ji,jj)*zqa )
325            !
326         END DO
327      END DO
328      !
329   END FUNCTION One_on_L
330
331
332   FUNCTION psi_m_coare( pzeta )
333      !!----------------------------------------------------------------------------------
334      !! ** Purpose: compute the universal profile stability function for momentum
335      !!             COARE 3.0, Fairall et al. 2003
336      !!             pzeta : stability paramenter, z/L where z is altitude
337      !!                     measurement and L is M-O length
338      !!       Stability function for wind speed and scalars matching Kansas and free
339      !!       convection forms with weighting f convective form, follows Fairall et
340      !!       al (1996) with profile constants from Grachev et al (2000) BLM stable
341      !!       form from Beljaars and Holtslag (1991)
342      !!
343      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
344      !!----------------------------------------------------------------------------------
345      REAL(wp), DIMENSION(jpi,jpj) :: psi_m_coare
346      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
347      !
348      INTEGER  ::   ji, jj    ! dummy loop indices
349      REAL(wp) :: zta, zphi_m, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
350      !!----------------------------------------------------------------------------------
351      !
352      DO jj = 1, jpj
353         DO ji = 1, jpi
354            !
355            zta = pzeta(ji,jj)
356            !
357            zphi_m = ABS(1. - 15.*zta)**.25    !!Kansas unstable
358            !
359            zpsi_k = 2.*LOG((1. + zphi_m)/2.) + LOG((1. + zphi_m*zphi_m)/2.)   &
360               & - 2.*ATAN(zphi_m) + 0.5*rpi
361            !
362            zphi_c = ABS(1. - 10.15*zta)**.3333                   !!Convective
363            !
364            zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
365               &     - 1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
366            !
367            zf = zta*zta
368            zf = zf/(1. + zf)
369            zc = MIN(50., 0.35*zta)
370            zstab = 0.5 + SIGN(0.5, zta)
371            !
372            psi_m_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & ! (zta < 0)
373               &                -   zstab     * ( 1. + 1.*zta     &                ! (zta > 0)
374               &                         + 0.6667*(zta - 14.28)/EXP(zc) + 8.525 )   !     "
375            !
376         END DO
377      END DO
378      !
379   END FUNCTION psi_m_coare
380
381
382   FUNCTION psi_h_coare( pzeta )
383      !!---------------------------------------------------------------------
384      !! Universal profile stability function for temperature and humidity
385      !! COARE 3.0, Fairall et al. 2003
386      !!
387      !! pzeta : stability paramenter, z/L where z is altitude measurement
388      !!         and L is M-O length
389      !!
390      !! Stability function for wind speed and scalars matching Kansas and free
391      !! convection forms with weighting f convective form, follows Fairall et
392      !! al (1996) with profile constants from Grachev et al (2000) BLM stable
393      !! form from Beljaars and Holtslag (1991)
394      !!
395      !! Author: L. Brodeau, june 2016 / AeroBulk
396      !!         (https://sourceforge.net/p/aerobulk)
397      !!----------------------------------------------------------------
398      !!
399      REAL(wp), DIMENSION(jpi,jpj) :: psi_h_coare
400      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
401      !
402      INTEGER  ::   ji, jj     ! dummy loop indices
403      REAL(wp) :: zta, zphi_h, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
404      !
405      DO jj = 1, jpj
406         DO ji = 1, jpi
407            !
408            zta = pzeta(ji,jj)
409            !
410            zphi_h = (ABS(1. - 15.*zta))**.5  !! Kansas unstable   (zphi_h = zphi_m**2 when unstable, zphi_m when stable)
411            !
412            zpsi_k = 2.*LOG((1. + zphi_h)/2.)
413            !
414            zphi_c = (ABS(1. - 34.15*zta))**.3333   !! Convective
415            !
416            zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
417               &    -1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
418            !
419            zf = zta*zta
420            zf = zf/(1. + zf)
421            zc = MIN(50.,0.35*zta)
422            zstab = 0.5 + SIGN(0.5, zta)
423            !
424            psi_h_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) &
425               &                -   zstab     * ( (ABS(1. + 2.*zta/3.))**1.5     &
426               &                           + .6667*(zta - 14.28)/EXP(zc) + 8.525 )
427            !
428         END DO
429      END DO
430      !
431   END FUNCTION psi_h_coare
432
433
434   FUNCTION visc_air( ptak )
435      !!---------------------------------------------------------------------
436      !! Air kinetic viscosity (m^2/s) given from temperature in degrees...
437      !!
438      !! Author: L. Brodeau, june 2016 / AeroBulk
439      !!         (https://sourceforge.net/p/aerobulk)
440      !!---------------------------------------------------------------------
441      !!
442      REAL(wp), DIMENSION(jpi,jpj) :: visc_air
443      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak  ! air temperature in (K)
444      !
445      INTEGER  ::   ji, jj         ! dummy loop indices
446      REAL(wp) ::   ztc, ztc2      ! local scalar
447      !
448      DO jj = 1, jpj
449         DO ji = 1, jpi
450            ztc  = ptak(ji,jj) - rt0   ! air temp, in deg. C
451            ztc2 = ztc*ztc
452            visc_air(ji,jj) = 1.326E-5*(1. + 6.542E-3*ztc + 8.301E-6*ztc2 - 4.84E-9*ztc2*ztc)
453         END DO
454      END DO
455      !
456   END FUNCTION visc_air
457
458
459END MODULE sbcblk_algo_coare
Note: See TracBrowser for help on using the repository browser.