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 @ 11266

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

LB: CSWL param now uses NEMO time step (rdt) and previous-t value of t_skin as first guess.
  • Property svn:keywords set to Id
File size: 21.9 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 2019 / 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 iom             ! I/O manager library
29   USE lib_mpp         ! distribued memory computing library
30   USE in_out_manager  ! I/O manager
31   USE prtctl          ! Print control
32   USE sbcwave, ONLY   :  cdn_wave ! wave module
33#if defined key_si3 || defined key_cice
34   USE sbc_ice         ! Surface boundary condition: ice fields
35#endif
36   USE lib_fortran     ! to use key_nosignedzero
37
38   USE sbc_oce         ! Surface boundary condition: ocean fields
39   USE sbcblk_phy      ! all thermodynamics functions, rho_air, q_sat, etc... !LB
40   USE sbcblk_skin     ! cool-skin/warm layer scheme (CSWL_ECMWF) !LB
41
42   IMPLICIT NONE
43   PRIVATE
44
45   PUBLIC ::   TURB_COARE   ! called by sbcblk.F90
46
47   !                                              !! COARE own values for given constants:
48   REAL(wp), PARAMETER ::   zi0     = 600._wp      ! scale height of the atmospheric boundary layer...
49   REAL(wp), PARAMETER ::   Beta0   =   1.250_wp   ! gustiness parameter
50
51   INTEGER , PARAMETER ::   nb_itt = 5             ! number of itterations
52
53   !!----------------------------------------------------------------------
54CONTAINS
55
56   SUBROUTINE turb_coare( zt, zu, T_s, t_zt, q_s, q_zt, U_zu, &
57      &                   Cd, Ch, Ce, t_zu, q_zu, U_blk,      &
58      &                   Cdn, Chn, Cen,                      &
59      &                   Qsw, rad_lw, slp, Tsk_b             )
60      !!----------------------------------------------------------------------
61      !!                      ***  ROUTINE  turb_coare  ***
62      !!
63      !! ** Purpose :   Computes turbulent transfert coefficients of surface
64      !!                fluxes according to Fairall et al. (2003)
65      !!                If relevant (zt /= zu), adjust temperature and humidity from height zt to zu
66      !!                Returns the effective bulk wind speed at 10m to be used in the bulk formulas
67      !!
68      !!                Applies the cool-skin warm-layer correction of the SST to T_s
69      !!                if the net shortwave flux at the surface (Qsw), the downwelling longwave
70      !!                radiative fluxes at the surface (rad_lw), and the sea-leve pressure (slp)
71      !!                are provided as (optional) arguments!
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      !!    *  t_zt : potential air temperature at zt                         [K]
79      !!    *  q_zt : specific humidity of air at zt                          [kg/kg]
80      !!
81      !! INPUT/OUTPUT:
82      !! -------------
83      !!    *  T_s  : always "bulk SST" as input                              [K]
84      !!              -> unchanged "bulk SST" as output if CSWL not used      [K]
85      !!              -> skin temperature as output if CSWL used              [K]
86      !!
87      !!    *  q_s  : SSQ aka saturation specific humidity at temp. T_s       [kg/kg]
88      !!              -> doesn't need to be given a value if skin temp computed (in case l_use_skin=True)
89      !!              -> MUST be given the correct value if not computing skint temp. (in case l_use_skin=False)
90      !!
91      !! OPTIONAL INPUT (will trigger l_use_skin=TRUE if present!):
92      !! ---------------
93      !!    *  Qsw    : net solar flux (after albedo) at the surface (>0)     [W/m^2]
94      !!    *  rad_lw : downwelling longwave radiation at the surface  (>0)   [W/m^2]
95      !!    *  slp    : sea-level pressure                                    [Pa]
96      !!    *  Tsk_b  : estimate of skin temperature at previous time-step    [K]
97      !!
98      !! OUTPUT :
99      !! --------
100      !!    *  Cd     : drag coefficient
101      !!    *  Ch     : sensible heat coefficient
102      !!    *  Ce     : evaporation coefficient
103      !!    *  t_zu   : pot. air temperature adjusted at wind height zu       [K]
104      !!    *  q_zu   : specific humidity of air        //                    [kg/kg]
105      !!    *  U_blk  : bulk wind speed at 10m                                [m/s]
106      !!
107      !!
108      !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/)
109      !!----------------------------------------------------------------------------------
110      REAL(wp), INTENT(in   )                     ::   zt       ! height for t_zt and q_zt                    [m]
111      REAL(wp), INTENT(in   )                     ::   zu       ! height for U_zu                             [m]
112      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) ::   T_s      ! sea surface temperature                [Kelvin]
113      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   t_zt     ! potential air temperature              [Kelvin]
114      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) ::   q_s      ! sea surface specific humidity           [kg/kg]
115      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_zt     ! specific air humidity at zt             [kg/kg]
116      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   U_zu     ! relative wind module at zu                [m/s]
117      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cd       ! transfer coefficient for momentum         (tau)
118      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ch       ! transfer coefficient for sensible heat (Q_sens)
119      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ce       ! transfert coefficient for evaporation   (Q_lat)
120      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   t_zu     ! pot. air temp. adjusted at zu               [K]
121      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   q_zu     ! spec. humidity adjusted at zu           [kg/kg]
122      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   U_blk    ! bulk wind at 10m                          [m/s]
123      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients
124      !
125      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   Qsw      !             [W/m^2]
126      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   rad_lw   !             [W/m^2]
127      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   slp      !             [Pa]
128      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   Tsk_b    !             [Pa]
129      !
130      INTEGER :: j_itt
131      LOGICAL :: l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U
132      !
133      REAL(wp), DIMENSION(jpi,jpj) ::  &
134         &  u_star, t_star, q_star, &
135         &  dt_zu, dq_zu,    &
136         &  znu_a,           & !: Nu_air, Viscosity of air
137         &  z0, z0t
138      REAL(wp), DIMENSION(jpi,jpj) ::   zeta_u        ! stability parameter at height zu
139      REAL(wp), DIMENSION(jpi,jpj) ::   ztmp0, ztmp1, ztmp2
140      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zeta_t        ! stability parameter at height zt
141      !
142      ! Cool skin:
143      LOGICAL :: l_use_skin = .FALSE.
144      REAL(wp), DIMENSION(jpi,jpj) :: zsst    ! to back up the initial bulk SST
145      !!----------------------------------------------------------------------------------
146      !
147      ! Cool skin ?
148      IF( PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp) ) THEN
149         l_use_skin = .TRUE.
150      END IF
151      IF (lwp) PRINT *, ' *** LOLO(sbcblk_algo_ecmwf.F90) => l_use_skin =', l_use_skin
152
153      l_zt_equal_zu = .FALSE.
154      IF( ABS(zu - zt) < 0.01_wp )   l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision
155
156      IF( .NOT. l_zt_equal_zu )  ALLOCATE( zeta_t(jpi,jpj) )
157
158      !! Initialization for cool skin:
159      zsst   = T_s    ! save the bulk SST
160      IF( l_use_skin ) THEN
161         ! First guess for skin temperature:
162         IF( PRESENT(Tsk_b) ) THEN
163            T_s = Tsk_b
164         ELSE
165            T_s = T_s - 0.25     ! sst - 0.25
166         END IF
167         q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s
168      END IF
169
170      !! First guess of temperature and humidity at height zu:
171      t_zu = MAX( t_zt ,  180._wp )   ! who knows what's given on masked-continental regions...
172      q_zu = MAX( q_zt , 1.e-6_wp )   !               "
173
174      !! Pot. temp. difference (and we don't want it to be 0!)
175      dt_zu = t_zu - T_s ;   dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
176      dq_zu = q_zu - q_s ;   dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
177
178      znu_a = visc_air(t_zu) ! Air viscosity (m^2/s) at zt given from temperature in (K)
179
180      ztmp2 = 0.5_wp*0.5_wp  ! initial guess for wind gustiness contribution
181      U_blk = SQRT(U_zu*U_zu + ztmp2)
182
183      ztmp2   = 10000._wp     ! optimization: ztmp2 == 1/z0 (with z0 first guess == 0.0001)
184      ztmp0   = LOG(zu*ztmp2)
185      ztmp1   = LOG(10.*ztmp2)
186      u_star = 0.035_wp*U_blk*ztmp1/ztmp0       ! (u* = 0.035*Un10)
187
188      z0     = alfa_charn(U_blk)*u_star*u_star/grav + 0.11_wp*znu_a/u_star
189      z0     = MIN(ABS(z0), 0.001_wp)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO
190      z0t    = 1._wp / ( 0.1_wp*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) )
191      z0t    = MIN(ABS(z0t), 0.001_wp)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO
192
193      ztmp2  = vkarmn/ztmp0
194      Cd     = ztmp2*ztmp2    ! first guess of Cd
195
196      ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd
197
198      ztmp2 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN)
199
200      !! First estimate of zeta_u, depending on the stability, ie sign of BRN (ztmp2):
201      ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 )
202      ztmp0 = ztmp0*ztmp2
203      zeta_u = (1._wp-ztmp1) * (ztmp0/(1._wp+ztmp2/(-zu/(zi0*0.004_wp*Beta0**3)))) & !  BRN < 0
204         &  +     ztmp1   * (ztmp0*(1._wp + 27._wp/9._wp*ztmp2/ztmp0))               !  BRN > 0
205      !#LB: should make sure that the "ztmp0" of "27./9.*ztmp2/ztmp0" is "ztmp0*ztmp2" and not "ztmp0==vkarmn*vkarmn/LOG(zt/z0t)/Cd" !
206
207      !! First guess M-O stability dependent scaling params.(u*,t*,q*) to estimate z0 and z/L
208      ztmp0  = vkarmn/(LOG(zu/z0t) - psi_h_coare(zeta_u))
209
210      u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0)  - psi_m_coare(zeta_u))
211      t_star = dt_zu*ztmp0
212      q_star = dq_zu*ztmp0
213
214      ! What needs to be done if zt /= zu:
215      IF( .NOT. l_zt_equal_zu ) THEN
216         !! First update of values at zu (or zt for wind)
217         zeta_t = zt*zeta_u/zu
218         ztmp0 = psi_h_coare(zeta_u) - psi_h_coare(zeta_t)
219         ztmp1 = LOG(zt/zu) + ztmp0
220         t_zu = t_zt - t_star/vkarmn*ztmp1
221         q_zu = q_zt - q_star/vkarmn*ztmp1
222         q_zu = (0.5_wp + SIGN(0.5_wp,q_zu))*q_zu !Makes it impossible to have negative humidity :
223         !
224         dt_zu = t_zu - T_s  ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
225         dq_zu = q_zu - q_s  ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
226      END IF
227
228      !! ITERATION BLOCK
229      DO j_itt = 1, nb_itt
230
231         !!Inverse of Monin-Obukov length (1/L) :
232         ztmp0 = One_on_L(t_zu, q_zu, u_star, t_star, q_star)  ! 1/L == 1/[Monin-Obukhov length]
233         ztmp0 = SIGN( MIN(ABS(ztmp0),200._wp), ztmp0 ) ! (prevents FPE from stupid values from masked region later on...) !#LOLO
234
235         ztmp1 = u_star*u_star   ! u*^2
236
237         !! Update wind at 10m taking into acount convection-related wind gustiness:
238         ! Ug = Beta*w*  (Beta = 1.25, Fairall et al. 2003, Eq.8):
239         ztmp2 = Beta0*Beta0*ztmp1*(MAX(-zi0*ztmp0/vkarmn,0._wp))**(2./3.)   ! => ztmp2 == Ug^2
240         !!   ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before 600.
241         U_blk = MAX(sqrt(U_zu*U_zu + ztmp2), 0.2_wp)        ! include gustiness in bulk wind speed
242         ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0.
243
244         !! Updating Charnock parameter, increases with the wind (Fairall et al., 2003 p. 577-578)
245         ztmp2 = alfa_charn(U_blk)  ! alpha Charnock parameter
246
247         !! Roughness lengthes z0, z0t (z0q = z0t) :
248         z0    = ztmp2*ztmp1/grav + 0.11_wp*znu_a/u_star  ! Roughness length (eq.6)
249         ztmp1 = z0*u_star/znu_a                          ! Re_r: roughness Reynolds number
250         z0t  = MIN( 1.1E-4_wp , 5.5E-5_wp*ztmp1**(-0.6_wp) ) ! Scalar roughness for both theta and q (eq.28)
251
252         !! Stability parameters:
253         zeta_u = zu*ztmp0
254         zeta_u = SIGN( MIN(ABS(zeta_u),50.0_wp), zeta_u )
255         IF( .NOT. l_zt_equal_zu ) THEN
256            zeta_t = zt*ztmp0
257            zeta_t = SIGN( MIN(ABS(zeta_t),50.0_wp), zeta_t )
258         END IF
259
260         !! Turbulent scales at zu=10m :
261         ztmp0   = psi_h_coare(zeta_u)
262         ztmp1   = vkarmn/(LOG(zu) - LOG(z0t) - ztmp0)
263
264         t_star = dt_zu*ztmp1
265         q_star = dq_zu*ztmp1
266         u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0) - psi_m_coare(zeta_u))
267
268         IF( .NOT. l_zt_equal_zu ) THEN
269            ! What's need to be done if zt /= zu
270            !! Re-updating temperature and humidity at zu :
271            ztmp2 = ztmp0 - psi_h_coare(zeta_t)
272            ztmp1 = log(zt/zu) + ztmp2
273            t_zu = t_zt - t_star/vkarmn*ztmp1
274            q_zu = q_zt - q_star/vkarmn*ztmp1
275         END IF
276
277         !! SKIN related part
278         !! -----------------
279         IF( l_use_skin ) THEN
280            !! compute transfer coefficients at zu : lolo: verifier...
281            ztmp0 = u_star/U_blk
282            Ch   = ztmp0*t_star/dt_zu
283            Ce   = ztmp0*q_star/dq_zu
284            ! Non-Solar heat flux to the ocean:
285            ztmp1 = U_blk*MAX(rho_air(t_zu, q_zu, slp), 1._wp)     ! rho*U10
286            ztmp2 = T_s*T_s
287            ztmp1 = ztmp1 * ( Ce*L_vap(T_s)*(q_zu - q_s) + Ch*cp_air(q_zu)*(t_zu - T_s) ) & ! Total turb. heat flux
288               &       +      rad_lw - emiss_w*stefan*ztmp2*ztmp2                           ! Net longwave flux
289            !!         => "ztmp1" is the net non-solar surface heat flux !
290            !! Updating the values of the skin temperature T_s and q_s :
291            CALL CSWL_ECMWF( Qsw, ztmp1, u_star, zsst, T_s ) ! yes ECMWF, because more advanced than COARE (warm-layer added!)
292            q_s = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp)  ! 200 -> just to avoid numerics problem on masked regions if silly values are given
293         END IF
294
295         IF( (l_use_skin).OR.(.NOT. l_zt_equal_zu) ) THEN
296            dt_zu = t_zu - T_s ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
297            dq_zu = q_zu - q_s ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
298         END IF
299
300      END DO !DO j_itt = 1, nb_itt
301
302      ! compute transfer coefficients at zu :
303      ztmp0 = u_star/U_blk
304      Cd   = ztmp0*ztmp0
305      Ch   = ztmp0*t_star/dt_zu
306      Ce   = ztmp0*q_star/dq_zu
307      !
308      ztmp1 = zu + z0
309      Cdn = vkarmn*vkarmn / (log(ztmp1/z0 )*log(ztmp1/z0 ))
310      Chn = vkarmn*vkarmn / (log(ztmp1/z0t)*log(ztmp1/z0t))
311      Cen = Chn
312      !
313      IF( .NOT. l_zt_equal_zu ) DEALLOCATE( zeta_t )
314      !
315   END SUBROUTINE turb_coare
316
317
318   FUNCTION alfa_charn( pwnd )
319      !!-------------------------------------------------------------------
320      !! Compute the Charnock parameter as a function of the wind speed
321      !!
322      !! (Fairall et al., 2003 p.577-578)
323      !!
324      !! Wind below 10 m/s :  alfa = 0.011
325      !! Wind between 10 and 18 m/s : linear increase from 0.011 to 0.018
326      !! Wind greater than 18 m/s :  alfa = 0.018
327      !!
328      !! Author: L. Brodeau, June 2016 / AeroBulk  (https://github.com/brodeau/aerobulk/)
329      !!-------------------------------------------------------------------
330      REAL(wp), DIMENSION(jpi,jpj) :: alfa_charn
331      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pwnd   ! wind speed
332      !
333      INTEGER  ::   ji, jj         ! dummy loop indices
334      REAL(wp) :: zw, zgt10, zgt18
335      !!-------------------------------------------------------------------
336      !
337      DO jj = 1, jpj
338         DO ji = 1, jpi
339            !
340            zw = pwnd(ji,jj)   ! wind speed
341            !
342            ! Charnock's constant, increases with the wind :
343            zgt10 = 0.5 + SIGN(0.5_wp,(zw - 10)) ! If zw<10. --> 0, else --> 1
344            zgt18 = 0.5 + SIGN(0.5_wp,(zw - 18.)) ! If zw<18. --> 0, else --> 1
345            !
346            alfa_charn(ji,jj) =  (1. - zgt10)*0.011    &    ! wind is lower than 10 m/s
347               &     + zgt10*((1. - zgt18)*(0.011 + (0.018 - 0.011) &
348               &      *(zw - 10.)/(18. - 10.)) + zgt18*( 0.018 ) )    ! Hare et al. (1999)
349            !
350         END DO
351      END DO
352      !
353   END FUNCTION alfa_charn
354
355   FUNCTION psi_m_coare( pzeta )
356      !!----------------------------------------------------------------------------------
357      !! ** Purpose: compute the universal profile stability function for momentum
358      !!             COARE 3.0, Fairall et al. 2003
359      !!             pzeta : stability paramenter, z/L where z is altitude
360      !!                     measurement and L is M-O length
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 (https://github.com/brodeau/aerobulk/)
367      !!----------------------------------------------------------------------------------
368      REAL(wp), DIMENSION(jpi,jpj) :: psi_m_coare
369      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
370      !
371      INTEGER  ::   ji, jj    ! dummy loop indices
372      REAL(wp) :: zta, zphi_m, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
373      !!----------------------------------------------------------------------------------
374      !
375      DO jj = 1, jpj
376         DO ji = 1, jpi
377            !
378            zta = pzeta(ji,jj)
379            !
380            zphi_m = ABS(1. - 15.*zta)**.25    !!Kansas unstable
381            !
382            zpsi_k = 2.*LOG((1. + zphi_m)/2.) + LOG((1. + zphi_m*zphi_m)/2.)   &
383               & - 2.*ATAN(zphi_m) + 0.5*rpi
384            !
385            zphi_c = ABS(1. - 10.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._wp, 0.35_wp*zta)
393            zstab = 0.5 + SIGN(0.5_wp, zta)
394            !
395            psi_m_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & ! (zta < 0)
396               &                -   zstab     * ( 1. + 1.*zta     &                ! (zta > 0)
397               &                         + 0.6667*(zta - 14.28)/EXP(zc) + 8.525 )   !     "
398            !
399         END DO
400      END DO
401      !
402   END FUNCTION psi_m_coare
403
404
405   FUNCTION psi_h_coare( pzeta )
406      !!---------------------------------------------------------------------
407      !! Universal profile stability function for temperature and humidity
408      !! COARE 3.0, Fairall et al. 2003
409      !!
410      !! pzeta : stability paramenter, z/L where z is altitude measurement
411      !!         and L is M-O length
412      !!
413      !! Stability function for wind speed and scalars matching Kansas and free
414      !! convection forms with weighting f convective form, follows Fairall et
415      !! al (1996) with profile constants from Grachev et al (2000) BLM stable
416      !! form from Beljaars and Holtslag (1991)
417      !!
418      !! Author: L. Brodeau, June 2016 / AeroBulk
419      !!         (https://github.com/brodeau/aerobulk/)
420      !!----------------------------------------------------------------
421      !!
422      REAL(wp), DIMENSION(jpi,jpj) :: psi_h_coare
423      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
424      !
425      INTEGER  ::   ji, jj     ! dummy loop indices
426      REAL(wp) :: zta, zphi_h, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
427      !
428      DO jj = 1, jpj
429         DO ji = 1, jpi
430            !
431            zta = pzeta(ji,jj)
432            !
433            zphi_h = (ABS(1. - 15.*zta))**.5  !! Kansas unstable   (zphi_h = zphi_m**2 when unstable, zphi_m when stable)
434            !
435            zpsi_k = 2.*LOG((1. + zphi_h)/2.)
436            !
437            zphi_c = (ABS(1. - 34.15*zta))**.3333   !! Convective
438            !
439            zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
440               &    -1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
441            !
442            zf = zta*zta
443            zf = zf/(1. + zf)
444            zc = MIN(50._wp,0.35_wp*zta)
445            zstab = 0.5 + SIGN(0.5_wp, zta)
446            !
447            psi_h_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) &
448               &                -   zstab     * ( (ABS(1. + 2.*zta/3.))**1.5     &
449               &                           + .6667*(zta - 14.28)/EXP(zc) + 8.525 )
450            !
451         END DO
452      END DO
453      !
454   END FUNCTION psi_h_coare
455
456   !!======================================================================
457END MODULE sbcblk_algo_coare
Note: See TracBrowser for help on using the repository browser.