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 NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC – NEMO

source: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare.F90 @ 11182

Last change on this file since 11182 was 11182, checked in by laurent, 5 years ago

LB: adding possibility to use dew-point temperature rather than specific humidity as "air humidity" => ln_humi_dpt@namsbc_blk
  • Property svn:keywords set to Id
File size: 19.8 KB
Line 
1MODULE sbcblk_algo_coare
2   !!======================================================================
3   !!                   ***  MODULE  sbcblk_algo_coare  ***
4   !! Computes:
5   !!   * bulk transfer coefficients C_D, C_E and C_H
6   !!   * air temp. and spec. hum. adjusted from zt (2m) to zu (10m) if needed
7   !!   * the effective bulk wind speed at 10m U_blk
8   !!   => all these are used in bulk formulas in sbcblk.F90
9   !!
10   !!    Using the bulk formulation/param. of COARE v3, Fairall et al., 2003
11   !!
12   !!       Routine turb_coare maintained and developed in AeroBulk
13   !!                     (https://github.com/brodeau/aerobulk)
14   !!
15   !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk)
16   !!----------------------------------------------------------------------
17   !! History :  4.0  !  2016-02  (L.Brodeau)   Original code
18   !!----------------------------------------------------------------------
19
20   !!----------------------------------------------------------------------
21   !!   turb_coare  : computes the bulk turbulent transfer coefficients
22   !!                   adjusts t_air and q_air from zt to zu m
23   !!                   returns the effective bulk wind speed at 10m
24   !!----------------------------------------------------------------------
25   USE oce             ! ocean dynamics and tracers
26   USE dom_oce         ! ocean space and time domain
27   USE phycst          ! physical constants
28   USE sbc_oce         ! Surface boundary condition: ocean fields
29   USE sbcwave, ONLY   :  cdn_wave ! wave module
30#if defined key_si3 || defined key_cice
31   USE sbc_ice         ! Surface boundary condition: ice fields
32#endif
33   USE sbcblk_phy     !LB: all thermodynamics functions, rho_air, q_sat, etc...
34   USE in_out_manager ! I/O manager
35   USE iom            ! I/O manager library
36   USE lib_mpp        ! distribued memory computing library
37   USE prtctl         ! Print control
38   USE lib_fortran    ! to use key_nosignedzero
39
40   IMPLICIT NONE
41   PRIVATE
42
43   PUBLIC ::   TURB_COARE   ! called by sbcblk.F90
44
45   !                                              !! COARE own values for given constants:
46   REAL(wp), PARAMETER ::   zi0     = 600._wp      ! scale height of the atmospheric boundary layer...
47   REAL(wp), PARAMETER ::   Beta0   =   1.250_wp   ! gustiness parameter
48
49   INTEGER , PARAMETER ::   nb_itt = 5             ! number of itterations
50
51   !!----------------------------------------------------------------------
52CONTAINS
53
54   SUBROUTINE turb_coare( zt, zu, sst, t_zt, ssq, q_zt, U_zu, &
55      &                   Cd, Ch, Ce, t_zu, q_zu, U_blk,      &
56      &                   Cdn, Chn, Cen                       )
57      !!----------------------------------------------------------------------
58      !!                      ***  ROUTINE  turb_coare  ***
59      !!
60      !! ** Purpose :   Computes turbulent transfert coefficients of surface
61      !!                fluxes according to Fairall et al. (2003)
62      !!                If relevant (zt /= zu), adjust temperature and humidity from height zt to zu
63      !!                Returns the effective bulk wind speed at 10m to be used in the bulk formulas
64      !!
65      !!
66      !! INPUT :
67      !! -------
68      !!    *  zt   : height for temperature and spec. hum. of air            [m]
69      !!    *  zu   : height for wind speed (generally 10m)                   [m]
70      !!    *  U_zu : scalar wind speed at 10m                                [m/s]
71      !!    *  sst  : SST                                                     [K]
72      !!    *  t_zt : potential air temperature at zt                         [K]
73      !!    *  ssq  : specific humidity at saturation at SST                  [kg/kg]
74      !!    *  q_zt : specific humidity of air at zt                          [kg/kg]
75      !!
76      !!
77      !! OUTPUT :
78      !! --------
79      !!    *  Cd     : drag coefficient
80      !!    *  Ch     : sensible heat coefficient
81      !!    *  Ce     : evaporation coefficient
82      !!    *  t_zu   : pot. air temperature adjusted at wind height zu       [K]
83      !!    *  q_zu   : specific humidity of air        //                    [kg/kg]
84      !!    *  U_blk  : bulk wind speed at 10m                                [m/s]
85      !!
86      !!
87      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk)
88      !!----------------------------------------------------------------------------------
89      REAL(wp), INTENT(in   )                     ::   zt       ! height for t_zt and q_zt                    [m]
90      REAL(wp), INTENT(in   )                     ::   zu       ! height for U_zu                             [m]
91      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   sst      ! sea surface temperature                [Kelvin]
92      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   t_zt     ! potential air temperature              [Kelvin]
93      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   ssq      ! sea surface specific humidity           [kg/kg]
94      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_zt     ! specific air humidity at zt             [kg/kg]
95      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   U_zu     ! relative wind module at zu                [m/s]
96      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cd       ! transfer coefficient for momentum         (tau)
97      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ch       ! transfer coefficient for sensible heat (Q_sens)
98      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ce       ! transfert coefficient for evaporation   (Q_lat)
99      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   t_zu     ! pot. air temp. adjusted at zu               [K]
100      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   q_zu     ! spec. humidity adjusted at zu           [kg/kg]
101      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   U_blk    ! bulk wind at 10m                          [m/s]
102      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients
103      !
104      INTEGER :: j_itt
105      LOGICAL :: l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U
106      !
107      REAL(wp), DIMENSION(jpi,jpj) ::  &
108         &  u_star, t_star, q_star, &
109         &  dt_zu, dq_zu,    &
110         &  znu_a,           & !: Nu_air, Viscosity of air
111         &  z0, z0t
112      REAL(wp), DIMENSION(jpi,jpj) ::   zeta_u        ! stability parameter at height zu
113      REAL(wp), DIMENSION(jpi,jpj) ::   ztmp0, ztmp1, ztmp2
114      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zeta_t        ! stability parameter at height zt
115      !!----------------------------------------------------------------------
116      !
117      l_zt_equal_zu = .FALSE.
118      IF( ABS(zu - zt) < 0.01 ) l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision
119
120      IF( .NOT. l_zt_equal_zu )  ALLOCATE( zeta_t(jpi,jpj) )
121
122      !! First guess of temperature and humidity at height zu:
123      t_zu = MAX( t_zt , 199.0_wp )   ! who knows what's given on masked-continental regions...
124      q_zu = MAX( q_zt , 1.e-6_wp )   !               "
125
126      !! Pot. temp. difference (and we don't want it to be 0!)
127      dt_zu = t_zu - sst ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
128      dq_zu = q_zu - ssq ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
129
130      znu_a = visc_air(t_zu) ! Air viscosity (m^2/s) at zt given from temperature in (K)
131
132      ztmp2 = 0.5*0.5  ! initial guess for wind gustiness contribution
133      U_blk = SQRT(U_zu*U_zu + ztmp2)
134
135      ztmp2   = 10000.     ! optimization: ztmp2 == 1/z0 (with z0 first guess == 0.0001)
136      ztmp0   = LOG(zu*ztmp2)
137      ztmp1   = LOG(10.*ztmp2)
138      u_star = 0.035*U_blk*ztmp1/ztmp0       ! (u* = 0.035*Un10)
139
140
141      z0     = alfa_charn(U_blk)*u_star*u_star/grav + 0.11*znu_a/u_star
142      z0     = MIN(ABS(z0), 0.001)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO
143      z0t    = 1. / ( 0.1*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) )
144      z0t    = MIN(ABS(z0t), 0.001)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO
145
146      ztmp2  = vkarmn/ztmp0
147      Cd     = ztmp2*ztmp2    ! first guess of Cd
148
149      ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd
150
151      ztmp2  = grav*zu*(dt_zu + rctv0*t_zu*dq_zu)/(t_zu*U_blk*U_blk)  !! Ribu Bulk Richardson number ;       !Ribcu = -zu/(zi0*0.004*Beta0**3) !! Saturation Rib, zi0 = tropicalbound. layer depth
152
153      !! First estimate of zeta_u, depending on the stability, ie sign of Ribu (ztmp2):
154      ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 )
155      ztmp0 = ztmp0*ztmp2
156      zeta_u = (1.-ztmp1) * (ztmp0/(1.+ztmp2/(-zu/(zi0*0.004*Beta0**3)))) & !  Ribu < 0
157         &  +     ztmp1   * (ztmp0*(1. + 27./9.*ztmp2/ztmp0))               !  Ribu > 0
158      !#LOLO: should make sure that the "ztmp0" of "27./9.*ztmp2/ztmp0" is "ztmp0*ztmp2" and not "ztmp0==vkarmn*vkarmn/LOG(zt/z0t)/Cd" !
159
160      !! First guess M-O stability dependent scaling params.(u*,t*,q*) to estimate z0 and z/L
161      ztmp0  = vkarmn/(LOG(zu/z0t) - psi_h_coare(zeta_u))
162
163      u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0)  - psi_m_coare(zeta_u))
164      t_star = dt_zu*ztmp0
165      q_star = dq_zu*ztmp0
166
167      ! What's need to be done if zt /= zu:
168      IF( .NOT. l_zt_equal_zu ) THEN
169         !! First update of values at zu (or zt for wind)
170         zeta_t = zt*zeta_u/zu
171         ztmp0 = psi_h_coare(zeta_u) - psi_h_coare(zeta_t)
172         ztmp1 = LOG(zt/zu) + ztmp0
173         t_zu = t_zt - t_star/vkarmn*ztmp1
174         q_zu = q_zt - q_star/vkarmn*ztmp1
175         q_zu = (0.5 + SIGN(0.5_wp,q_zu))*q_zu !Makes it impossible to have negative humidity :
176         !
177         dt_zu = t_zu - sst  ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
178         dq_zu = q_zu - ssq  ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
179      END IF
180
181      !! ITERATION BLOCK
182      DO j_itt = 1, nb_itt
183
184         !!Inverse of Monin-Obukov length (1/L) :
185         ztmp0 = One_on_L(t_zu, q_zu, u_star, t_star, q_star)  ! 1/L == 1/[Monin-Obukhov length]
186         ztmp0 = SIGN( MIN(ABS(ztmp0),200._wp), ztmp0 ) ! (prevents FPE from stupid values from masked region later on...) !#LOLO
187
188         ztmp1 = u_star*u_star   ! u*^2
189
190         !! Update wind at 10m taking into acount convection-related wind gustiness:
191         ! Ug = Beta*w*  (Beta = 1.25, Fairall et al. 2003, Eq.8):
192         ztmp2 = Beta0*Beta0*ztmp1*(MAX(-zi0*ztmp0/vkarmn,0._wp))**(2./3.)   ! => ztmp2 == Ug^2
193         !!   ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before 600.
194         U_blk = MAX(sqrt(U_zu*U_zu + ztmp2), 0.2_wp)        ! include gustiness in bulk wind speed
195         ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0.
196
197         !! Updating Charnock parameter, increases with the wind (Fairall et al., 2003 p. 577-578)
198         ztmp2 = alfa_charn(U_blk)  ! alpha Charnock parameter
199
200         !! Roughness lengthes z0, z0t (z0q = z0t) :
201         z0   = ztmp2*ztmp1/grav + 0.11*znu_a/u_star ! Roughness length (eq.6)
202         ztmp1 = z0*u_star/znu_a                           ! Re_r: roughness Reynolds number
203         z0t  = MIN( 1.1E-4_wp , 5.5E-5_wp*ztmp1**(-0.6_wp) ) ! Scalar roughness for both theta and q (eq.28)
204
205         !! Stability parameters:
206         zeta_u = zu*ztmp0
207         zeta_u = SIGN( MIN(ABS(zeta_u),50.0_wp), zeta_u )
208         IF( .NOT. l_zt_equal_zu ) THEN
209            zeta_t = zt*ztmp0
210            zeta_t = SIGN( MIN(ABS(zeta_t),50.0_wp), zeta_t )
211         END IF
212
213         !! Turbulent scales at zu=10m :
214         ztmp0   = psi_h_coare(zeta_u)
215         ztmp1   = vkarmn/(LOG(zu) - LOG(z0t) - ztmp0)
216
217         t_star = dt_zu*ztmp1
218         q_star = dq_zu*ztmp1
219         u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0) - psi_m_coare(zeta_u))
220
221         IF( .NOT. l_zt_equal_zu ) THEN
222            ! What's need to be done if zt /= zu
223            !! Re-updating temperature and humidity at zu :
224            ztmp2 = ztmp0 - psi_h_coare(zeta_t)
225            ztmp1 = log(zt/zu) + ztmp2
226            t_zu = t_zt - t_star/vkarmn*ztmp1
227            q_zu = q_zt - q_star/vkarmn*ztmp1
228            dt_zu = t_zu - sst ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6), dt_zu )
229            dq_zu = q_zu - ssq ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9), dq_zu )
230         END IF
231
232      END DO
233      !
234      ! compute transfer coefficients at zu :
235      ztmp0 = u_star/U_blk
236      Cd   = ztmp0*ztmp0
237      Ch   = ztmp0*t_star/dt_zu
238      Ce   = ztmp0*q_star/dq_zu
239      !
240      ztmp1 = zu + z0
241      Cdn = vkarmn*vkarmn / (log(ztmp1/z0 )*log(ztmp1/z0 ))
242      Chn = vkarmn*vkarmn / (log(ztmp1/z0t)*log(ztmp1/z0t))
243      Cen = Chn
244      !
245      IF( .NOT. l_zt_equal_zu ) DEALLOCATE( zeta_t )
246      !
247   END SUBROUTINE turb_coare
248
249
250   FUNCTION alfa_charn( pwnd )
251      !!-------------------------------------------------------------------
252      !! Compute the Charnock parameter as a function of the wind speed
253      !!
254      !! (Fairall et al., 2003 p.577-578)
255      !!
256      !! Wind below 10 m/s :  alfa = 0.011
257      !! Wind between 10 and 18 m/s : linear increase from 0.011 to 0.018
258      !! Wind greater than 18 m/s :  alfa = 0.018
259      !!
260      !! Author: L. Brodeau, june 2016 / AeroBulk  (https://github.com/brodeau/aerobulk/)
261      !!-------------------------------------------------------------------
262      REAL(wp), DIMENSION(jpi,jpj) :: alfa_charn
263      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pwnd   ! wind speed
264      !
265      INTEGER  ::   ji, jj         ! dummy loop indices
266      REAL(wp) :: zw, zgt10, zgt18
267      !!-------------------------------------------------------------------
268      !
269      DO jj = 1, jpj
270         DO ji = 1, jpi
271            !
272            zw = pwnd(ji,jj)   ! wind speed
273            !
274            ! Charnock's constant, increases with the wind :
275            zgt10 = 0.5 + SIGN(0.5_wp,(zw - 10)) ! If zw<10. --> 0, else --> 1
276            zgt18 = 0.5 + SIGN(0.5_wp,(zw - 18.)) ! If zw<18. --> 0, else --> 1
277            !
278            alfa_charn(ji,jj) =  (1. - zgt10)*0.011    &    ! wind is lower than 10 m/s
279               &     + zgt10*((1. - zgt18)*(0.011 + (0.018 - 0.011) &
280               &      *(zw - 10.)/(18. - 10.)) + zgt18*( 0.018 ) )    ! Hare et al. (1999)
281            !
282         END DO
283      END DO
284      !
285   END FUNCTION alfa_charn
286
287
288   FUNCTION One_on_L( ptha, pqa, pus, pts, pqs )
289      !!------------------------------------------------------------------------
290      !!
291      !! Evaluates the 1./(Monin Obukhov length) from air temperature and
292      !!  specific humidity, and frictional scales u*, t* and q*
293      !!
294      !! Author: L. Brodeau, june 2016 / AeroBulk
295      !!         (https://github.com/brodeau/aerobulk/)
296      !!------------------------------------------------------------------------
297      REAL(wp), DIMENSION(jpi,jpj)             :: One_on_L         !: 1./(Monin Obukhov length) [m^-1]
298      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptha,  &  !: average potetntial air temperature [K]
299         &                                        pqa,   &  !: average specific humidity of air   [kg/kg]
300         &                                      pus, pts, pqs   !: frictional velocity, temperature and humidity
301      !
302      INTEGER  ::   ji, jj         ! dummy loop indices
303      REAL(wp) ::     zqa          ! local scalar
304      !!-------------------------------------------------------------------
305      !
306      DO jj = 1, jpj
307         DO ji = 1, jpi
308            !
309            zqa = (1. + rctv0*pqa(ji,jj))
310            !
311            One_on_L(ji,jj) =  grav*vkarmn*(pts(ji,jj)*zqa + rctv0*ptha(ji,jj)*pqs(ji,jj)) &
312               &                      / ( pus(ji,jj)*pus(ji,jj) * ptha(ji,jj)*zqa )
313            !
314         END DO
315      END DO
316      !
317   END FUNCTION One_on_L
318
319
320   FUNCTION psi_m_coare( pzeta )
321      !!----------------------------------------------------------------------------------
322      !! ** Purpose: compute the universal profile stability function for momentum
323      !!             COARE 3.0, Fairall et al. 2003
324      !!             pzeta : stability paramenter, z/L where z is altitude
325      !!                     measurement and L is M-O length
326      !!       Stability function for wind speed and scalars matching Kansas and free
327      !!       convection forms with weighting f convective form, follows Fairall et
328      !!       al (1996) with profile constants from Grachev et al (2000) BLM stable
329      !!       form from Beljaars and Holtslag (1991)
330      !!
331      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
332      !!----------------------------------------------------------------------------------
333      REAL(wp), DIMENSION(jpi,jpj) :: psi_m_coare
334      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
335      !
336      INTEGER  ::   ji, jj    ! dummy loop indices
337      REAL(wp) :: zta, zphi_m, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
338      !!----------------------------------------------------------------------------------
339      !
340      DO jj = 1, jpj
341         DO ji = 1, jpi
342            !
343            zta = pzeta(ji,jj)
344            !
345            zphi_m = ABS(1. - 15.*zta)**.25    !!Kansas unstable
346            !
347            zpsi_k = 2.*LOG((1. + zphi_m)/2.) + LOG((1. + zphi_m*zphi_m)/2.)   &
348               & - 2.*ATAN(zphi_m) + 0.5*rpi
349            !
350            zphi_c = ABS(1. - 10.15*zta)**.3333                   !!Convective
351            !
352            zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
353               &     - 1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
354            !
355            zf = zta*zta
356            zf = zf/(1. + zf)
357            zc = MIN(50._wp, 0.35_wp*zta)
358            zstab = 0.5 + SIGN(0.5_wp, zta)
359            !
360            psi_m_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & ! (zta < 0)
361               &                -   zstab     * ( 1. + 1.*zta     &                ! (zta > 0)
362               &                         + 0.6667*(zta - 14.28)/EXP(zc) + 8.525 )   !     "
363            !
364         END DO
365      END DO
366      !
367   END FUNCTION psi_m_coare
368
369
370   FUNCTION psi_h_coare( pzeta )
371      !!---------------------------------------------------------------------
372      !! Universal profile stability function for temperature and humidity
373      !! COARE 3.0, Fairall et al. 2003
374      !!
375      !! pzeta : stability paramenter, z/L where z is altitude measurement
376      !!         and L is M-O length
377      !!
378      !! Stability function for wind speed and scalars matching Kansas and free
379      !! convection forms with weighting f convective form, follows Fairall et
380      !! al (1996) with profile constants from Grachev et al (2000) BLM stable
381      !! form from Beljaars and Holtslag (1991)
382      !!
383      !! Author: L. Brodeau, june 2016 / AeroBulk
384      !!         (https://github.com/brodeau/aerobulk/)
385      !!----------------------------------------------------------------
386      !!
387      REAL(wp), DIMENSION(jpi,jpj) :: psi_h_coare
388      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
389      !
390      INTEGER  ::   ji, jj     ! dummy loop indices
391      REAL(wp) :: zta, zphi_h, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
392      !
393      DO jj = 1, jpj
394         DO ji = 1, jpi
395            !
396            zta = pzeta(ji,jj)
397            !
398            zphi_h = (ABS(1. - 15.*zta))**.5  !! Kansas unstable   (zphi_h = zphi_m**2 when unstable, zphi_m when stable)
399            !
400            zpsi_k = 2.*LOG((1. + zphi_h)/2.)
401            !
402            zphi_c = (ABS(1. - 34.15*zta))**.3333   !! Convective
403            !
404            zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
405               &    -1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
406            !
407            zf = zta*zta
408            zf = zf/(1. + zf)
409            zc = MIN(50._wp,0.35_wp*zta)
410            zstab = 0.5 + SIGN(0.5_wp, zta)
411            !
412            psi_h_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) &
413               &                -   zstab     * ( (ABS(1. + 2.*zta/3.))**1.5     &
414               &                           + .6667*(zta - 14.28)/EXP(zc) + 8.525 )
415            !
416         END DO
417      END DO
418      !
419   END FUNCTION psi_h_coare
420
421   !!======================================================================
422END MODULE sbcblk_algo_coare
Note: See TracBrowser for help on using the repository browser.