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

source: branches/UKMO/dev_merge_2017_CICE_interface/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_coare.F90 @ 9819

Last change on this file since 9819 was 9125, checked in by timgraham, 6 years ago

Removed wrk_arrays from whole code. No change in SETTE results from this.

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