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

source: branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_coare3p5.F90 @ 8637

Last change on this file since 8637 was 8637, checked in by gm, 7 years ago

#1911 (ENHANCE-09): PART I.3 - phasing with updated branch dev_r8183_ICEMODEL revision 8626
File size: 19.9 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_lim3 || 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 wrk_nemo        ! work arrays
41   USE timing          ! Timing
42   USE in_out_manager  ! I/O manager
43   USE prtctl          ! Print control
44   USE lib_fortran     ! to use key_nosignedzero
45
46   IMPLICIT NONE
47   PRIVATE
48
49   PUBLIC ::   TURB_COARE3P5   ! called by sbcblk.F90
50
51   !                                   ! COARE own values for given constants:
52   REAL(wp), PARAMETER ::   charn0_max =   0.028   ! value above which the Charnock paramter levels off for winds > 18
53   REAL(wp), PARAMETER ::   zi0        = 600.      ! scale height of the atmospheric boundary layer...1
54   REAL(wp), PARAMETER ::   Beta0      =   1.25    ! gustiness parameter
55   REAL(wp), PARAMETER ::   rctv0      =   0.608   ! constant to obtain virtual temperature...
56
57   !!----------------------------------------------------------------------
58CONTAINS
59
60   SUBROUTINE turb_coare3p5( zt, zu, sst, t_zt, ssq, q_zt, U_zu,  &
61      &                      Cd, Ch, Ce, t_zu, q_zu, U_blk,       &
62      &                      Cdn, Chn, Cen                        )
63      !!----------------------------------------------------------------------------------
64      !!                      ***  ROUTINE  turb_coare3p5  ***
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      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
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      !!----------------------------------------------------------------------------------
95      REAL(wp), INTENT(in   )                     ::   zt       ! height for t_zt and q_zt                   [m]
96      REAL(wp), INTENT(in   )                     ::   zu       ! height for U_zu                              [m]
97      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   sst      ! sea surface temperature              [Kelvin]
98      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   t_zt     ! potential air temperature            [Kelvin]
99      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   ssq      ! sea surface specific humidity         [kg/kg]
100      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_zt     ! specific air humidity                 [kg/kg]
101      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   U_zu     ! relative wind module at zu            [m/s]
102      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cd       ! transfer coefficient for momentum         (tau)
103      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ch       ! transfer coefficient for sensible heat (Q_sens)
104      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ce       ! transfert coefficient for evaporation   (Q_lat)
105      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   t_zu     ! pot. air temp. adjusted at zu             [K]
106      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   q_zu     ! spec. humidity adjusted at zu             [kg/kg]
107      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   U_blk    ! bulk wind at 10m                          [m/s]
108      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients
109      !
110      INTEGER :: j_itt
111      LOGICAL ::   l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U
112      INTEGER , PARAMETER ::   nb_itt = 4       ! number of itterations
113      !
114      REAL(wp), DIMENSION(:,:), POINTER  ::  &
115         &  u_star, t_star, q_star, &
116         &  dt_zu, dq_zu,    &
117         &  znu_a,           & !: Nu_air, Viscosity of air
118         &  z0, z0t
119      REAL(wp), DIMENSION(:,:), POINTER ::   zeta_u        ! stability parameter at height zu
120      REAL(wp), DIMENSION(:,:), POINTER ::   zeta_t        ! stability parameter at height zt
121      REAL(wp), DIMENSION(:,:), POINTER ::   ztmp0, ztmp1, ztmp2
122      !!----------------------------------------------------------------------------------
123      !
124      IF( nn_timing == 1 )  CALL timing_start('turb_coare3p5')
125
126      CALL wrk_alloc( jpi,jpj, u_star, t_star, q_star, zeta_u, dt_zu, dq_zu)
127      CALL wrk_alloc( jpi,jpj, znu_a, z0, z0t, ztmp0, ztmp1, ztmp2 )
128
129      l_zt_equal_zu = .FALSE.
130      IF( ABS(zu - zt) < 0.01 ) l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision
131
132      IF( .NOT. l_zt_equal_zu )   CALL wrk_alloc( jpi,jpj, zeta_t )
133
134      !! First guess of temperature and humidity at height zu:
135      t_zu = MAX(t_zt , 0.0)    ! who knows what's given on masked-continental regions...
136      q_zu = MAX(q_zt , 1.E-6)  !               "
137
138      !! Pot. temp. difference (and we don't want it to be 0!)
139      dt_zu = t_zu - sst ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6), dt_zu )
140      dq_zu = q_zu - ssq ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9), dq_zu )
141
142      znu_a = visc_air(t_zt) ! Air viscosity (m^2/s) at zt given from temperature in (K)
143
144      ztmp2 = 0.5*0.5  ! initial guess for wind gustiness contribution
145      U_blk = SQRT(U_zu*U_zu + ztmp2)
146
147      ztmp2   = 10000.     ! optimization: ztmp2 == 1/z0 (with z0 first guess == 0.0001)
148      ztmp0   = LOG(zu*ztmp2)
149      ztmp1   = LOG(10.*ztmp2)
150      u_star = 0.035*U_blk*ztmp1/ztmp0       ! (u* = 0.035*Un10)
151
152      !! COARE 3.5 first guess of UN10 is U_zu
153      ztmp2 = MIN( 0.0017*U_zu - 0.005 , charn0_max)  ! alpha Charnock parameter (Eq. 13 Edson al. 2013)
154      ztmp2 = MAX( ztmp2 , 0. )                       ! alpha Charnock parameter (Eq. 13 Edson al. 2013)
155      z0     = ztmp2*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         !! COARE 3.5: Charnock parameter is computed from the neutral wind speed at 10m: Eq. 13 (Edson al. 2013)
211         ztmp2 = u_star/vkarmn*LOG(10./z0)   ! UN10 Neutral wind at 10m!
212         ztmp2 = MIN( 0.0017*ztmp2 - 0.005 , charn0_max)  ! alpha Charnock parameter (Eq. 13 Edson al. 2013)
213         ztmp2 = MAX( ztmp2 , 0. )
214
215         !! Roughness lengthes z0, z0t (z0q = z0t) :
216         z0   = ztmp2*ztmp1/grav + 0.11*znu_a/u_star ! Roughness length (eq.6)
217         ztmp1 = z0*u_star/znu_a                             ! Re_r: roughness Reynolds number
218         !z0t  = MIN( 1.1E-4 , 5.5E-5*ztmp1**(-0.6) )   ! COARE 3.0
219         !! Chris Fairall and Jim Edsson, private communication, March 2016 / COARE 3.5 :
220         z0t   = MIN( 1.6e-4 , 5.8E-5*ztmp1**(-0.72)) ! These thermal roughness lengths give Stanton and
221         !z0q = z0t                                   ! Dalton numbers that closely approximate COARE3.0
222
223         !! Stability parameters:
224         zeta_u = zu*ztmp0 ; zeta_u = sign( min(abs(zeta_u),50.0), zeta_u )
225         IF( .NOT. l_zt_equal_zu ) THEN
226            zeta_t = zt*ztmp0 ;  zeta_t = sign( min(abs(zeta_t),50.0), zeta_t )
227         END IF
228
229         !! Turbulent scales at zu=10m :
230         ztmp0   = psi_h_coare(zeta_u)
231         ztmp1   = vkarmn/(LOG(zu) -LOG(z0t) - ztmp0)
232
233         t_star = dt_zu*ztmp1
234         q_star = dq_zu*ztmp1
235         u_star = U_blk*vkarmn/(LOG(zu) -LOG(z0) - psi_m_coare(zeta_u))
236
237         IF( .NOT. l_zt_equal_zu ) THEN
238            ! What's need to be done if zt /= zu
239            !! Re-updating temperature and humidity at zu :
240            ztmp2 = ztmp0 - psi_h_coare(zeta_t)
241            ztmp1 = log(zt/zu) + ztmp2
242            t_zu = t_zt - t_star/vkarmn*ztmp1
243            q_zu = q_zt - q_star/vkarmn*ztmp1
244            dt_zu = t_zu - sst ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6), dt_zu )
245            dq_zu = q_zu - ssq ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9), dq_zu )
246         END IF
247
248      END DO
249      !
250      ! compute transfer coefficients at zu :
251      ztmp0 = u_star/U_blk
252      Cd   = ztmp0*ztmp0
253      Ch   = ztmp0*t_star/dt_zu
254      Ce   = ztmp0*q_star/dq_zu
255      !
256      ztmp1 = zu + z0
257      Cdn = vkarmn*vkarmn / (log(ztmp1/z0 )*log(ztmp1/z0 ))
258      Chn = vkarmn*vkarmn / (log(ztmp1/z0t)*log(ztmp1/z0t))
259      Cen = Chn
260      !
261      CALL wrk_dealloc( jpi,jpj, u_star, t_star, q_star, zeta_u, dt_zu, dq_zu )
262      CALL wrk_dealloc( jpi,jpj, znu_a, z0, z0t, ztmp0, ztmp1, ztmp2 )
263      IF( .NOT. l_zt_equal_zu ) CALL wrk_dealloc( jpi,jpj, zeta_t )
264
265      IF( nn_timing == 1 )  CALL timing_stop('turb_coare3p5')
266
267   END SUBROUTINE turb_coare3p5
268
269
270
271   FUNCTION One_on_L( ptha, pqa, pus, pts, pqs )
272      !!------------------------------------------------------------------------
273      !!
274      !! Evaluates the 1./(Monin Obukhov length) from air temperature and
275      !!  specific humidity, and frictional scales u*, t* and q*
276      !!
277      !! Author: L. Brodeau, june 2016 / AeroBulk
278      !!         (https://sourceforge.net/p/aerobulk)
279      !!------------------------------------------------------------------------
280      REAL(wp), DIMENSION(jpi,jpj)             :: One_on_L         !: 1./(Monin Obukhov length) [m^-1]
281      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptha,  &  !: average potetntial air temperature [K]
282         &                                        pqa,   &  !: average specific humidity of air   [kg/kg]
283         &                                      pus, pts, pqs   !: frictional velocity, temperature and humidity
284      !
285      INTEGER  ::   ji, jj         ! dummy loop indices
286      REAL(wp) ::     zqa          ! local scalar
287      !
288      DO jj = 1, jpj
289         DO ji = 1, jpi
290            !
291            zqa = (1. + rctv0*pqa(ji,jj))
292            !
293            One_on_L(ji,jj) =  grav*vkarmn*(pts(ji,jj)*zqa + rctv0*ptha(ji,jj)*pqs(ji,jj)) &
294               &                      / ( pus(ji,jj)*pus(ji,jj) * ptha(ji,jj)*zqa )
295            !
296         END DO
297      END DO
298      !
299   END FUNCTION One_on_L
300
301
302   FUNCTION psi_m_coare( pzeta )
303      !!----------------------------------------------------------------------------------
304      !! ** Purpose: compute the universal profile stability function for momentum
305      !!             COARE 3.0, Fairall et al. 2003
306      !!             pzeta : stability paramenter, z/L where z is altitude
307      !!                     measurement and L is M-O length
308      !!       Stability function for wind speed and scalars matching Kansas and free
309      !!       convection forms with weighting f convective form, follows Fairall et
310      !!       al (1996) with profile constants from Grachev et al (2000) BLM stable
311      !!       form from Beljaars and Holtslag (1991)
312      !!
313      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
314      !!----------------------------------------------------------------------------------
315      REAL(wp), DIMENSION(jpi,jpj) :: psi_m_coare
316      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
317      !
318      INTEGER  ::   ji, jj    ! dummy loop indices
319      REAL(wp) :: zta, zphi_m, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
320      !!----------------------------------------------------------------------------------
321      !
322      DO jj = 1, jpj
323         DO ji = 1, jpi
324            !
325            zta = pzeta(ji,jj)
326            !
327            zphi_m = ABS(1. - 15.*zta)**.25    !!Kansas unstable
328            !
329            zpsi_k = 2.*LOG((1. + zphi_m)/2.) + LOG((1. + zphi_m*zphi_m)/2.)   &
330               & - 2.*ATAN(zphi_m) + 0.5*rpi
331            !
332            zphi_c = ABS(1. - 10.15*zta)**.3333                   !!Convective
333            !
334            zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
335               &     - 1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
336            !
337            zf = zta*zta
338            zf = zf/(1. + zf)
339            zc = MIN(50., 0.35*zta)
340            zstab = 0.5 + SIGN(0.5, zta)
341            !
342            psi_m_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & ! (zta < 0)
343               &                -   zstab     * ( 1. + 1.*zta     &                ! (zta > 0)
344               &                         + 0.6667*(zta - 14.28)/EXP(zc) + 8.525 )   !     "
345            !
346         END DO
347      END DO
348      !
349   END FUNCTION psi_m_coare
350
351
352   FUNCTION psi_h_coare( pzeta )
353      !!---------------------------------------------------------------------
354      !! Universal profile stability function for temperature and humidity
355      !! COARE 3.0, Fairall et al. 2003
356      !!
357      !! pzeta : stability paramenter, z/L where z is altitude measurement
358      !!         and L is M-O length
359      !!
360      !! Stability function for wind speed and scalars matching Kansas and free
361      !! convection forms with weighting f convective form, follows Fairall et
362      !! al (1996) with profile constants from Grachev et al (2000) BLM stable
363      !! form from Beljaars and Holtslag (1991)
364      !!
365      !! Author: L. Brodeau, june 2016 / AeroBulk
366      !!         (https://sourceforge.net/p/aerobulk)
367      !!----------------------------------------------------------------
368      !!
369      REAL(wp), DIMENSION(jpi,jpj) :: psi_h_coare
370      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
371      !
372      INTEGER  ::   ji, jj     ! dummy loop indices
373      REAL(wp) :: zta, zphi_h, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
374      !
375      DO jj = 1, jpj
376         DO ji = 1, jpi
377            !
378            zta = pzeta(ji,jj)
379            !
380            zphi_h = (ABS(1. - 15.*zta))**.5  !! Kansas unstable   (zphi_h = zphi_m**2 when unstable, zphi_m when stable)
381            !
382            zpsi_k = 2.*LOG((1. + zphi_h)/2.)
383            !
384            zphi_c = (ABS(1. - 34.15*zta))**.3333   !! Convective
385            !
386            zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
387               &    -1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
388            !
389            zf = zta*zta
390            zf = zf/(1. + zf)
391            zc = MIN(50.,0.35*zta)
392            zstab = 0.5 + SIGN(0.5, zta)
393            !
394            psi_h_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) &
395               &                -   zstab     * ( (ABS(1. + 2.*zta/3.))**1.5     &
396               &                           + .6667*(zta - 14.28)/EXP(zc) + 8.525 )
397            !
398         END DO
399      END DO
400      !
401   END FUNCTION psi_h_coare
402
403
404   FUNCTION visc_air( ptak )
405      !!---------------------------------------------------------------------
406      !! Air kinetic viscosity (m^2/s) given from temperature in degrees...
407      !!
408      !! Author: L. Brodeau, june 2016 / AeroBulk
409      !!         (https://sourceforge.net/p/aerobulk)
410      !!---------------------------------------------------------------------
411      REAL(wp), DIMENSION(jpi,jpj)             ::   visc_air
412      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak       ! air temperature   [K]
413      !
414      INTEGER  ::   ji, jj         ! dummy loop indices
415      REAL(wp) ::   ztc, ztc2      ! local scalar
416      !
417      DO jj = 1, jpj
418         DO ji = 1, jpi
419            ztc  = ptak(ji,jj) - rt0   ! air temp, in deg. C
420            ztc2 = ztc*ztc
421            visc_air(ji,jj) = 1.326E-5*(1. + 6.542E-3*ztc + 8.301E-6*ztc2 - 4.84E-9*ztc2*ztc)
422         END DO
423      END DO
424      !
425   END FUNCTION visc_air
426
427   !!======================================================================
428END MODULE sbcblk_algo_coare3p5
Note: See TracBrowser for help on using the repository browser.