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

source: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_coare.F90 @ 7968

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