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

Last change on this file since 11291 was 11291, checked in by laurent, 14 months ago

LB: more appropriated use of constant "emiss_w": emissivity of sea water.

File size: 21.5 KB
Line 
1MODULE sbcblk_algo_coare3p6
2   !!======================================================================
3   !!                   ***  MODULE  sbcblk_algo_coare3p6  ***
4   !!
5   !!                After Fairall et al 2018 & Edson et al 2013
6   !! Computes:
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   !!       Routine turb_coare3p6 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_coare3p6  : 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_COARE3P6   ! 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.2_wp   ! gustiness parameter
50
51   INTEGER , PARAMETER ::   nb_itt = 5             ! number of itterations
52
53   !!----------------------------------------------------------------------
54CONTAINS
55
56   SUBROUTINE turb_coare3p6( 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_coare3p6  ***
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 (usually 10m)                     [m]
77      !!    *  U_zu : scalar wind speed at zu                                 [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 zu                                 [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 speed at zu                     [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      U_blk = SQRT(U_zu*U_zu + 0.5_wp*0.5_wp) ! initial guess for wind gustiness contribution
181
182      ztmp0   = LOG(    zu*10000._wp) ! optimization: 10000. == 1/z0 (with z0 first guess == 0.0001)
183      ztmp1   = LOG(10._wp*10000._wp) !       "                    "               "
184      u_star = 0.035_wp*U_blk*ztmp1/ztmp0       ! (u* = 0.035*Un10)
185
186      z0     = alfa_charn_3p6(U_zu)*u_star*u_star/grav + 0.11_wp*znu_a/u_star
187      z0     = MIN(ABS(z0), 0.001_wp)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO
188      z0t    = 1._wp / ( 0.1_wp*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) )
189      z0t    = MIN(ABS(z0t), 0.001_wp)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO
190
191      ztmp2  = vkarmn/ztmp0
192      Cd     = ztmp2*ztmp2    ! first guess of Cd
193
194      ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd
195
196      ztmp2 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN)
197
198      !! First estimate of zeta_u, depending on the stability, ie sign of BRN (ztmp2):
199      ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 )
200      ztmp0 = ztmp0*ztmp2
201      zeta_u = (1._wp-ztmp1) * (ztmp0/(1._wp+ztmp2/(-zu/(zi0*0.004_wp*Beta0**3)))) & !  BRN < 0
202         &  +     ztmp1   * (ztmp0*(1._wp + 27._wp/9._wp*ztmp2/ztmp0))               !  BRN > 0
203      !#LB: should make sure that the "ztmp0" of "27./9.*ztmp2/ztmp0" is "ztmp0*ztmp2" and not "ztmp0==vkarmn*vkarmn/LOG(zt/z0t)/Cd" !
204
205      !! First guess M-O stability dependent scaling params.(u*,t*,q*) to estimate z0 and z/L
206      ztmp0  = vkarmn/(LOG(zu/z0t) - psi_h_coare(zeta_u))
207
208      u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0)  - psi_m_coare(zeta_u))
209      t_star = dt_zu*ztmp0
210      q_star = dq_zu*ztmp0
211
212      ! What needs to be done if zt /= zu:
213      IF( .NOT. l_zt_equal_zu ) THEN
214         !! First update of values at zu (or zt for wind)
215         zeta_t = zt*zeta_u/zu
216         ztmp0 = psi_h_coare(zeta_u) - psi_h_coare(zeta_t)
217         ztmp1 = LOG(zt/zu) + ztmp0
218         t_zu = t_zt - t_star/vkarmn*ztmp1
219         q_zu = q_zt - q_star/vkarmn*ztmp1
220         q_zu = (0.5_wp + SIGN(0.5_wp,q_zu))*q_zu !Makes it impossible to have negative humidity :
221         !
222         dt_zu = t_zu - T_s  ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
223         dq_zu = q_zu - q_s  ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
224      END IF
225
226      !! ITERATION BLOCK
227      DO j_itt = 1, nb_itt
228
229         !!Inverse of Monin-Obukov length (1/L) :
230         ztmp0 = One_on_L(t_zu, q_zu, u_star, t_star, q_star)  ! 1/L == 1/[Monin-Obukhov length]
231         ztmp0 = SIGN( MIN(ABS(ztmp0),200._wp), ztmp0 ) ! (prevents FPE from stupid values from masked region later on...) !#LOLO
232
233         ztmp1 = u_star*u_star   ! u*^2
234
235         !! Update wind at zu taking into acount convection-related wind gustiness:
236         ! Ug = Beta*w*  (Beta = 1.25, Fairall et al. 2003, Eq.8):
237         ztmp2 = Beta0*Beta0*ztmp1*(MAX(-zi0*ztmp0/vkarmn,0._wp))**(2./3.) ! square of wind gustiness contribution, ztmp2 == Ug^2
238         !!   ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before 600.
239         U_blk = MAX(sqrt(U_zu*U_zu + ztmp2), 0.2_wp)        ! include gustiness in bulk wind speed
240         ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0.
241
242         !! Stability parameters:
243         zeta_u = zu*ztmp0
244         zeta_u = SIGN( MIN(ABS(zeta_u),50.0_wp), zeta_u )
245         IF( .NOT. l_zt_equal_zu ) THEN
246            zeta_t = zt*ztmp0
247            zeta_t = SIGN( MIN(ABS(zeta_t),50.0_wp), zeta_t )
248         END IF
249
250         !! Adjustment the wind at 10m (not needed in the current algo form):
251         !IF ( zu \= 10._wp ) U10 = U_zu + u_star/vkarmn*(LOG(10._wp/zu) - psi_m_coare(10._wp*ztmp0) + psi_m_coare(zeta_u))
252         
253         !! Roughness lengthes z0, z0t (z0q = z0t) :
254         ztmp2 = u_star/vkarmn*LOG(10./z0)                                 ! Neutral wind speed at 10m
255         z0    = alfa_charn_3p6(ztmp2)*ztmp1/grav + 0.11_wp*znu_a/u_star   ! Roughness length (eq.6)
256         ztmp1 = z0*u_star/znu_a                                           ! Re_r: roughness Reynolds number
257         z0t   = MIN( 1.6E-4_wp , 5.8E-5_wp*ztmp1**(-0.72_wp))   ! COARE 3.6
258
259         !! Turbulent scales at zu :
260         ztmp0   = psi_h_coare(zeta_u)
261         ztmp1   = vkarmn/(LOG(zu) - LOG(z0t) - ztmp0) ! #LOLO: in ztmp0, some use psi_h_coare(zeta_t) rather than psi_h_coare(zeta_t) ???
262
263         t_star = dt_zu*ztmp1
264         q_star = dq_zu*ztmp1
265         u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0) - psi_m_coare(zeta_u))
266
267         IF( .NOT. l_zt_equal_zu ) THEN
268            ! What's need to be done if zt /= zu
269            !! Re-updating temperature and humidity at zu :
270            ztmp2 = ztmp0 - psi_h_coare(zeta_t)
271            ztmp1 = log(zt/zu) + ztmp2
272            t_zu = t_zt - t_star/vkarmn*ztmp1
273            q_zu = q_zt - q_star/vkarmn*ztmp1
274         END IF
275
276         !! SKIN related part
277         !! -----------------
278         IF( l_use_skin ) THEN
279            !! compute transfer coefficients at zu : lolo: verifier...
280            ztmp0 = u_star/U_blk
281            Ch   = ztmp0*t_star/dt_zu
282            Ce   = ztmp0*q_star/dq_zu
283            ! Non-Solar heat flux to the ocean:
284            ztmp1 = U_blk*MAX(rho_air(t_zu, q_zu, slp), 1._wp)     ! rho*U10
285            ztmp2 = T_s*T_s
286            ztmp1 = ztmp1 * ( Ce*L_vap(T_s)*(q_zu - q_s) + Ch*cp_air(q_zu)*(t_zu - T_s) ) & ! Total turb. heat flux
287               &       +      emiss_w*(rad_lw - stefan*ztmp2*ztmp2)                         ! Net longwave flux
288            !!         => "ztmp1" is the net non-solar surface heat flux !
289            !! Updating the values of the skin temperature T_s and q_s :
290            CALL CSWL_ECMWF( Qsw, ztmp1, u_star, zsst, T_s ) ! yes ECMWF, because more advanced than COARE (warm-layer added!)
291            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
292         END IF
293
294         IF( (l_use_skin).OR.(.NOT. l_zt_equal_zu) ) THEN
295            dt_zu = t_zu - T_s ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
296            dq_zu = q_zu - q_s ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
297         END IF
298
299      END DO !DO j_itt = 1, nb_itt
300
301      ! compute transfer coefficients at zu :
302      ztmp0 = u_star/U_blk
303      Cd   = ztmp0*ztmp0
304      Ch   = ztmp0*t_star/dt_zu
305      Ce   = ztmp0*q_star/dq_zu
306      !
307      ztmp1 = zu + z0
308      Cdn = vkarmn*vkarmn / (log(ztmp1/z0 )*log(ztmp1/z0 ))
309      Chn = vkarmn*vkarmn / (log(ztmp1/z0t)*log(ztmp1/z0t))
310      Cen = Chn
311      !
312      IF( .NOT. l_zt_equal_zu ) DEALLOCATE( zeta_t )
313      !
314   END SUBROUTINE turb_coare3p6
315
316
317   FUNCTION alfa_charn_3p6( pwnd )
318      !!-------------------------------------------------------------------
319      !! Computes the Charnock parameter as a function of the Neutral wind speed at 10m
320      !!
321      !!  (Eq. 13 in Edson et al., 2013)
322      !!
323      !! Author: L. Brodeau, July 2019 / AeroBulk  (https://github.com/brodeau/aerobulk/)
324      !!-------------------------------------------------------------------
325      REAL(wp), DIMENSION(jpi,jpj) :: alfa_charn_3p6
326      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pwnd   ! neutral wind speed at 10m
327      !
328      REAL(wp), PARAMETER :: charn0_max = 0.028  !: value above which the Charnock parameter levels off for winds > 18 m/s
329      !!-------------------------------------------------------------------
330      !
331      alfa_charn_3p6 = MAX( MIN( 0.0017_wp*pwnd - 0.005_wp , charn0_max) , 0._wp )
332      !
333   END FUNCTION alfa_charn_3p6
334
335   FUNCTION psi_m_coare( pzeta )
336      !!----------------------------------------------------------------------------------
337      !! ** Purpose: compute the universal profile stability function for momentum
338      !!             COARE 3.0, Fairall et al. 2003
339      !!             pzeta : stability paramenter, z/L where z is altitude
340      !!                     measurement and L is M-O length
341      !!       Stability function for wind speed and scalars matching Kansas and free
342      !!       convection forms with weighting f convective form, follows Fairall et
343      !!       al (1996) with profile constants from Grachev et al (2000) BLM stable
344      !!       form from Beljaars and Holtslag (1991)
345      !!
346      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
347      !!----------------------------------------------------------------------------------
348      REAL(wp), DIMENSION(jpi,jpj) :: psi_m_coare
349      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
350      !
351      INTEGER  ::   ji, jj    ! dummy loop indices
352      REAL(wp) :: zta, zphi_m, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
353      !!----------------------------------------------------------------------------------
354      !
355      DO jj = 1, jpj
356         DO ji = 1, jpi
357            !
358            zta = pzeta(ji,jj)
359            !
360            zphi_m = ABS(1. - 15.*zta)**.25    !!Kansas unstable
361            !
362            zpsi_k = 2.*LOG((1. + zphi_m)/2.) + LOG((1. + zphi_m*zphi_m)/2.)   &
363               & - 2.*ATAN(zphi_m) + 0.5*rpi
364            !
365            zphi_c = ABS(1. - 10.15*zta)**.3333                   !!Convective
366            !
367            zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
368               &     - 1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
369            !
370            zf = zta*zta
371            zf = zf/(1. + zf)
372            zc = MIN(50._wp, 0.35_wp*zta)
373            zstab = 0.5 + SIGN(0.5_wp, zta)
374            !
375            psi_m_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & ! (zta < 0)
376               &                -   zstab     * ( 1. + 1.*zta     &                ! (zta > 0)
377               &                         + 0.6667*(zta - 14.28)/EXP(zc) + 8.525 )   !     "
378            !
379         END DO
380      END DO
381      !
382   END FUNCTION psi_m_coare
383
384
385   FUNCTION psi_h_coare( pzeta )
386      !!---------------------------------------------------------------------
387      !! Universal profile stability function for temperature and humidity
388      !! COARE 3.0, Fairall et al. 2003
389      !!
390      !! pzeta : stability paramenter, z/L where z is altitude measurement
391      !!         and L is M-O length
392      !!
393      !! Stability function for wind speed and scalars matching Kansas and free
394      !! convection forms with weighting f convective form, follows Fairall et
395      !! al (1996) with profile constants from Grachev et al (2000) BLM stable
396      !! form from Beljaars and Holtslag (1991)
397      !!
398      !! Author: L. Brodeau, June 2016 / AeroBulk
399      !!         (https://github.com/brodeau/aerobulk/)
400      !!----------------------------------------------------------------
401      !!
402      REAL(wp), DIMENSION(jpi,jpj) :: psi_h_coare
403      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
404      !
405      INTEGER  ::   ji, jj     ! dummy loop indices
406      REAL(wp) :: zta, zphi_h, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
407      !
408      DO jj = 1, jpj
409         DO ji = 1, jpi
410            !
411            zta = pzeta(ji,jj)
412            !
413            zphi_h = (ABS(1. - 15.*zta))**.5  !! Kansas unstable   (zphi_h = zphi_m**2 when unstable, zphi_m when stable)
414            !
415            zpsi_k = 2.*LOG((1. + zphi_h)/2.)
416            !
417            zphi_c = (ABS(1. - 34.15*zta))**.3333   !! Convective
418            !
419            zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
420               &    -1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
421            !
422            zf = zta*zta
423            zf = zf/(1. + zf)
424            zc = MIN(50._wp,0.35_wp*zta)
425            zstab = 0.5 + SIGN(0.5_wp, zta)
426            !
427            psi_h_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) &
428               &                -   zstab     * ( (ABS(1. + 2.*zta/3.))**1.5     &
429               &                           + .6667*(zta - 14.28)/EXP(zc) + 8.525 )
430            !
431         END DO
432      END DO
433      !
434   END FUNCTION psi_h_coare
435
436   !!======================================================================
437END MODULE sbcblk_algo_coare3p6
Note: See TracBrowser for help on using the repository browser.