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/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_coare3p5.F90 @ 8879

Last change on this file since 8879 was 8879, checked in by frrh, 6 years ago

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