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

Last change on this file since 9124 was 9124, checked in by gm, 6 years ago

dev_merge_2017: ln_timing instead of nn_timing + restricted timing to nemo_init and routine called by step in OPA_SRC

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