source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_coare3p5.F90 @ 8184

Last change on this file since 8184 was 7646, checked in by timgraham, 4 years ago

Merge of dev_merge_2016 into trunk. UPDATE TO ARCHFILES NEEDED for XIOS2.
LIM_SRC_s/limrhg.F90 to follow in next commit due to change of kind (I'm unable to do it in this commit).
Merged using the following steps:

1) svn merge —reintegrate svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk .
2) Resolve minor conflicts in sette.sh and namelist_cfg for ORCA2LIM3 (due to a change in trunk after branch was created)
3) svn commit
4) svn switch svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
5) svn merge svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2016/dev_merge_2016 .
6) At this stage I checked out a clean copy of the branch to compare against what is about to be committed to the trunk.
6) svn commit #Commit code to the trunk

In this commit I have also reverted a change to Fcheck_archfile.sh which was causing problems on the Paris machine.

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