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_coare3p0.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_coare3p0.F90 @ 11284

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

LB: gooodbye "coare3p5", hello "coare3p6" + cleaner air humidity handling + support for Relative humidity.
File size: 22.3 KB
Line 
1MODULE sbcblk_algo_coare3p0
2   !!======================================================================
3   !!                   ***  MODULE  sbcblk_algo_coare3p0  ***
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_coare3p0 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_coare3p0  : 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_COARE3P0   ! 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.25_wp   ! gustiness parameter
50
51   INTEGER , PARAMETER ::   nb_itt = 5             ! number of itterations
52
53   !!----------------------------------------------------------------------
54CONTAINS
55
56   SUBROUTINE turb_coare3p0( 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_coare3p0  ***
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_3p0(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_3p0(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.1E-4_wp , 5.5E-5_wp*ztmp1**(-0.6_wp) ) ! Scalar roughness for both theta and q (eq.28) #LOLO: some use 1.15 not 1.1 !!!
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               &       +      rad_lw - emiss_w*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_coare3p0
315
316
317   FUNCTION alfa_charn_3p0( pwnd )
318      !!-------------------------------------------------------------------
319      !! Compute the Charnock parameter as a function of the wind speed
320      !!
321      !! (Fairall et al., 2003 p.577-578)
322      !!
323      !! Wind below 10 m/s :  alfa = 0.011
324      !! Wind between 10 and 18 m/s : linear increase from 0.011 to 0.018
325      !! Wind greater than 18 m/s :  alfa = 0.018
326      !!
327      !! Author: L. Brodeau, June 2016 / AeroBulk  (https://github.com/brodeau/aerobulk/)
328      !!-------------------------------------------------------------------
329      REAL(wp), DIMENSION(jpi,jpj) :: alfa_charn_3p0
330      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pwnd   ! wind speed
331      !
332      INTEGER  ::   ji, jj         ! dummy loop indices
333      REAL(wp) :: zw, zgt10, zgt18
334      !!-------------------------------------------------------------------
335      !
336      DO jj = 1, jpj
337         DO ji = 1, jpi
338            !
339            zw = pwnd(ji,jj)   ! wind speed
340            !
341            ! Charnock's constant, increases with the wind :
342            zgt10 = 0.5 + SIGN(0.5_wp,(zw - 10)) ! If zw<10. --> 0, else --> 1
343            zgt18 = 0.5 + SIGN(0.5_wp,(zw - 18.)) ! If zw<18. --> 0, else --> 1
344            !
345            alfa_charn_3p0(ji,jj) =  (1. - zgt10)*0.011    &    ! wind is lower than 10 m/s
346               &     + zgt10*((1. - zgt18)*(0.011 + (0.018 - 0.011) &
347               &      *(zw - 10.)/(18. - 10.)) + zgt18*( 0.018 ) )    ! Hare et al. (1999)
348            !
349         END DO
350      END DO
351      !
352   END FUNCTION alfa_charn_3p0
353
354   FUNCTION psi_m_coare( pzeta )
355      !!----------------------------------------------------------------------------------
356      !! ** Purpose: compute the universal profile stability function for momentum
357      !!             COARE 3.0, Fairall et al. 2003
358      !!             pzeta : stability paramenter, z/L where z is altitude
359      !!                     measurement and L is M-O length
360      !!       Stability function for wind speed and scalars matching Kansas and free
361      !!       convection forms with weighting f convective form, follows Fairall et
362      !!       al (1996) with profile constants from Grachev et al (2000) BLM stable
363      !!       form from Beljaars and Holtslag (1991)
364      !!
365      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
366      !!----------------------------------------------------------------------------------
367      REAL(wp), DIMENSION(jpi,jpj) :: psi_m_coare
368      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
369      !
370      INTEGER  ::   ji, jj    ! dummy loop indices
371      REAL(wp) :: zta, zphi_m, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
372      !!----------------------------------------------------------------------------------
373      !
374      DO jj = 1, jpj
375         DO ji = 1, jpi
376            !
377            zta = pzeta(ji,jj)
378            !
379            zphi_m = ABS(1. - 15.*zta)**.25    !!Kansas unstable
380            !
381            zpsi_k = 2.*LOG((1. + zphi_m)/2.) + LOG((1. + zphi_m*zphi_m)/2.)   &
382               & - 2.*ATAN(zphi_m) + 0.5*rpi
383            !
384            zphi_c = ABS(1. - 10.15*zta)**.3333                   !!Convective
385            !
386            zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
387               &     - 1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
388            !
389            zf = zta*zta
390            zf = zf/(1. + zf)
391            zc = MIN(50._wp, 0.35_wp*zta)
392            zstab = 0.5 + SIGN(0.5_wp, zta)
393            !
394            psi_m_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & ! (zta < 0)
395               &                -   zstab     * ( 1. + 1.*zta     &                ! (zta > 0)
396               &                         + 0.6667*(zta - 14.28)/EXP(zc) + 8.525 )   !     "
397            !
398         END DO
399      END DO
400      !
401   END FUNCTION psi_m_coare
402
403
404   FUNCTION psi_h_coare( pzeta )
405      !!---------------------------------------------------------------------
406      !! Universal profile stability function for temperature and humidity
407      !! COARE 3.0, Fairall et al. 2003
408      !!
409      !! pzeta : stability paramenter, z/L where z is altitude measurement
410      !!         and L is M-O length
411      !!
412      !! Stability function for wind speed and scalars matching Kansas and free
413      !! convection forms with weighting f convective form, follows Fairall et
414      !! al (1996) with profile constants from Grachev et al (2000) BLM stable
415      !! form from Beljaars and Holtslag (1991)
416      !!
417      !! Author: L. Brodeau, June 2016 / AeroBulk
418      !!         (https://github.com/brodeau/aerobulk/)
419      !!----------------------------------------------------------------
420      !!
421      REAL(wp), DIMENSION(jpi,jpj) :: psi_h_coare
422      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
423      !
424      INTEGER  ::   ji, jj     ! dummy loop indices
425      REAL(wp) :: zta, zphi_h, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
426      !
427      DO jj = 1, jpj
428         DO ji = 1, jpi
429            !
430            zta = pzeta(ji,jj)
431            !
432            zphi_h = (ABS(1. - 15.*zta))**.5  !! Kansas unstable   (zphi_h = zphi_m**2 when unstable, zphi_m when stable)
433            !
434            zpsi_k = 2.*LOG((1. + zphi_h)/2.)
435            !
436            zphi_c = (ABS(1. - 34.15*zta))**.3333   !! Convective
437            !
438            zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
439               &    -1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
440            !
441            zf = zta*zta
442            zf = zf/(1. + zf)
443            zc = MIN(50._wp,0.35_wp*zta)
444            zstab = 0.5 + SIGN(0.5_wp, zta)
445            !
446            psi_h_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) &
447               &                -   zstab     * ( (ABS(1. + 2.*zta/3.))**1.5     &
448               &                           + .6667*(zta - 14.28)/EXP(zc) + 8.525 )
449            !
450         END DO
451      END DO
452      !
453   END FUNCTION psi_h_coare
454
455   !!======================================================================
456END MODULE sbcblk_algo_coare3p0
Note: See TracBrowser for help on using the repository browser.