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_coare3p5.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_coare3p5.F90 @ 11111

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

LB: skin temperature now used by ECMWF bulk algorithm, and can be xios-ed ; new module "SBC/sbcblk_skin.F90" ; all physical constants defined in SBC now moved to "DOM/phycst.F90"; all air-properties and other ABL functions gathered in a new module: "SBC/sbcblk_phymbl.F90".
  • Property svn:keywords set to Id
File size: 19.3 KB
RevLine 
[6723]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
[11111]16   !!                     (https://github.com/brodeau/aerobulk/)
[6723]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
[9570]34#if defined key_si3 || defined key_cice
[6723]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
[6727]47   PUBLIC ::   TURB_COARE3P5   ! called by sbcblk.F90
[6723]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
[9019]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                        )
[6723]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      !!
[11111]69      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
[6723]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]
[9019]105      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients
[6723]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      !
[9125]111      REAL(wp), DIMENSION(jpi,jpj) ::  &
[6723]112         &  u_star, t_star, q_star, &
113         &  dt_zu, dq_zu,    &
114         &  znu_a,           & !: Nu_air, Viscosity of air
115         &  z0, z0t
[9125]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
[6723]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
[9125]124      IF( .NOT. l_zt_equal_zu )   ALLOCATE( zeta_t(jpi,jpj) )
[6723]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      !
[9019]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      !
[9125]253      IF( .NOT. l_zt_equal_zu ) DEALLOCATE( zeta_t )
[9124]254      !
[6723]255   END SUBROUTINE turb_coare3p5
256
257
258
259   FUNCTION One_on_L( ptha, pqa, pus, pts, pqs )
260      !!------------------------------------------------------------------------
261      !!
262      !! Evaluates the 1./(Monin Obukhov length) from air temperature and
263      !!  specific humidity, and frictional scales u*, t* and q*
264      !!
265      !! Author: L. Brodeau, june 2016 / AeroBulk
[11111]266      !!         (https://github.com/brodeau/aerobulk/)
[6723]267      !!------------------------------------------------------------------------
268      REAL(wp), DIMENSION(jpi,jpj)             :: One_on_L         !: 1./(Monin Obukhov length) [m^-1]
269      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptha,  &  !: average potetntial air temperature [K]
270         &                                        pqa,   &  !: average specific humidity of air   [kg/kg]
271         &                                      pus, pts, pqs   !: frictional velocity, temperature and humidity
272      !
273      INTEGER  ::   ji, jj         ! dummy loop indices
274      REAL(wp) ::     zqa          ! local scalar
275      !
276      DO jj = 1, jpj
277         DO ji = 1, jpi
278            !
279            zqa = (1. + rctv0*pqa(ji,jj))
280            !
281            One_on_L(ji,jj) =  grav*vkarmn*(pts(ji,jj)*zqa + rctv0*ptha(ji,jj)*pqs(ji,jj)) &
282               &                      / ( pus(ji,jj)*pus(ji,jj) * ptha(ji,jj)*zqa )
283            !
284         END DO
285      END DO
286      !
287   END FUNCTION One_on_L
288
289
290   FUNCTION psi_m_coare( pzeta )
291      !!----------------------------------------------------------------------------------
292      !! ** Purpose: compute the universal profile stability function for momentum
293      !!             COARE 3.0, Fairall et al. 2003
294      !!             pzeta : stability paramenter, z/L where z is altitude
295      !!                     measurement and L is M-O length
296      !!       Stability function for wind speed and scalars matching Kansas and free
297      !!       convection forms with weighting f convective form, follows Fairall et
298      !!       al (1996) with profile constants from Grachev et al (2000) BLM stable
299      !!       form from Beljaars and Holtslag (1991)
300      !!
[11111]301      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
[6723]302      !!----------------------------------------------------------------------------------
303      REAL(wp), DIMENSION(jpi,jpj) :: psi_m_coare
304      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
305      !
306      INTEGER  ::   ji, jj    ! dummy loop indices
307      REAL(wp) :: zta, zphi_m, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
308      !!----------------------------------------------------------------------------------
309      !
310      DO jj = 1, jpj
311         DO ji = 1, jpi
312            !
313            zta = pzeta(ji,jj)
314            !
315            zphi_m = ABS(1. - 15.*zta)**.25    !!Kansas unstable
316            !
317            zpsi_k = 2.*LOG((1. + zphi_m)/2.) + LOG((1. + zphi_m*zphi_m)/2.)   &
318               & - 2.*ATAN(zphi_m) + 0.5*rpi
319            !
320            zphi_c = ABS(1. - 10.15*zta)**.3333                   !!Convective
321            !
322            zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
323               &     - 1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
324            !
325            zf = zta*zta
326            zf = zf/(1. + zf)
327            zc = MIN(50., 0.35*zta)
328            zstab = 0.5 + SIGN(0.5, zta)
329            !
330            psi_m_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & ! (zta < 0)
331               &                -   zstab     * ( 1. + 1.*zta     &                ! (zta > 0)
332               &                         + 0.6667*(zta - 14.28)/EXP(zc) + 8.525 )   !     "
333            !
334         END DO
335      END DO
336      !
337   END FUNCTION psi_m_coare
338
339
340   FUNCTION psi_h_coare( pzeta )
341      !!---------------------------------------------------------------------
342      !! Universal profile stability function for temperature and humidity
343      !! COARE 3.0, Fairall et al. 2003
344      !!
345      !! pzeta : stability paramenter, z/L where z is altitude measurement
346      !!         and L is M-O length
347      !!
348      !! Stability function for wind speed and scalars matching Kansas and free
349      !! convection forms with weighting f convective form, follows Fairall et
350      !! al (1996) with profile constants from Grachev et al (2000) BLM stable
351      !! form from Beljaars and Holtslag (1991)
352      !!
353      !! Author: L. Brodeau, june 2016 / AeroBulk
[11111]354      !!         (https://github.com/brodeau/aerobulk/)
[6723]355      !!----------------------------------------------------------------
356      !!
357      REAL(wp), DIMENSION(jpi,jpj) :: psi_h_coare
358      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
359      !
360      INTEGER  ::   ji, jj     ! dummy loop indices
361      REAL(wp) :: zta, zphi_h, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
362      !
363      DO jj = 1, jpj
364         DO ji = 1, jpi
365            !
366            zta = pzeta(ji,jj)
367            !
368            zphi_h = (ABS(1. - 15.*zta))**.5  !! Kansas unstable   (zphi_h = zphi_m**2 when unstable, zphi_m when stable)
369            !
370            zpsi_k = 2.*LOG((1. + zphi_h)/2.)
371            !
372            zphi_c = (ABS(1. - 34.15*zta))**.3333   !! Convective
373            !
374            zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
375               &    -1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
376            !
377            zf = zta*zta
378            zf = zf/(1. + zf)
379            zc = MIN(50.,0.35*zta)
380            zstab = 0.5 + SIGN(0.5, zta)
381            !
382            psi_h_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) &
383               &                -   zstab     * ( (ABS(1. + 2.*zta/3.))**1.5     &
384               &                           + .6667*(zta - 14.28)/EXP(zc) + 8.525 )
385            !
386         END DO
387      END DO
388      !
389   END FUNCTION psi_h_coare
390
391
392   FUNCTION visc_air( ptak )
393      !!---------------------------------------------------------------------
394      !! Air kinetic viscosity (m^2/s) given from temperature in degrees...
395      !!
396      !! Author: L. Brodeau, june 2016 / AeroBulk
[11111]397      !!         (https://github.com/brodeau/aerobulk/)
[6723]398      !!---------------------------------------------------------------------
399      REAL(wp), DIMENSION(jpi,jpj)             ::   visc_air
400      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak       ! air temperature   [K]
401      !
402      INTEGER  ::   ji, jj         ! dummy loop indices
403      REAL(wp) ::   ztc, ztc2      ! local scalar
404      !
405      DO jj = 1, jpj
406         DO ji = 1, jpi
407            ztc  = ptak(ji,jj) - rt0   ! air temp, in deg. C
408            ztc2 = ztc*ztc
409            visc_air(ji,jj) = 1.326E-5*(1. + 6.542E-3*ztc + 8.301E-6*ztc2 - 4.84E-9*ztc2*ztc)
410         END DO
411      END DO
412      !
413   END FUNCTION visc_air
414
415   !!======================================================================
416END MODULE sbcblk_algo_coare3p5
Note: See TracBrowser for help on using the repository browser.