source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_coare.F90 @ 8585

Last change on this file since 8585 was 8585, checked in by clem, 3 years ago

add option Lupkes2015 for ice-atm drag

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