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

Last change on this file since 11329 was 11329, checked in by laurent, 18 months ago

LB: cosmetic syntax homogenization.

  • Property svn:keywords set to Id
File size: 21.9 KB
Line 
1MODULE sbcblk_algo_ecmwf
2   !!======================================================================
3   !!                   ***  MODULE  sbcblk_algo_ecmwf  ***
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 IFS of ECMWF (cycle 40r1)
11   !!         based on IFS doc (avaible online on the ECMWF's website)
12   !!
13   !!       Routine turb_ecmwf maintained and developed in AeroBulk
14   !!                     (https://github.com/brodeau/aerobulk)
15   !!
16   !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk)
17   !!----------------------------------------------------------------------
18   !! History :  4.0  !  2016-02  (L.Brodeau)   Original code
19   !!----------------------------------------------------------------------
20
21   !!----------------------------------------------------------------------
22   !!   turb_ecmwf  : computes the bulk turbulent transfer coefficients
23   !!                   adjusts t_air and q_air from zt to zu m
24   !!                   returns the effective bulk wind speed at 10m
25   !!----------------------------------------------------------------------
26   USE oce             ! ocean dynamics and tracers
27   USE dom_oce         ! ocean space and time domain
28   USE phycst          ! physical constants
29   USE iom             ! I/O manager library
30   USE lib_mpp         ! distribued memory computing library
31   USE in_out_manager  ! I/O manager
32   USE prtctl          ! Print control
33   USE sbcwave, ONLY   :  cdn_wave ! wave module
34#if defined key_si3 || defined key_cice
35   USE sbc_ice         ! Surface boundary condition: ice fields
36#endif
37   USE lib_fortran     ! to use key_nosignedzero
38
39   USE sbc_oce         ! Surface boundary condition: ocean fields
40   USE sbcblk_phy      ! all thermodynamics functions, rho_air, q_sat, etc... !LB
41   USE sbcblk_skin     ! cool-skin/warm layer scheme (CSWL_ECMWF) !LB
42
43   IMPLICIT NONE
44   PRIVATE
45
46   PUBLIC ::   TURB_ECMWF   ! called by sbcblk.F90
47
48   !                   !! ECMWF own values for given constants, taken form IFS documentation...
49   REAL(wp), PARAMETER ::   charn0 = 0.018    ! Charnock constant (pretty high value here !!!
50   !                                          !    =>  Usually 0.011 for moderate winds)
51   REAL(wp), PARAMETER ::   zi0     = 1000.   ! scale height of the atmospheric boundary layer...1
52   REAL(wp), PARAMETER ::   Beta0    = 1.     ! gustiness parameter ( = 1.25 in COAREv3)
53   REAL(wp), PARAMETER ::   alpha_M = 0.11    ! For roughness length (smooth surface term)
54   REAL(wp), PARAMETER ::   alpha_H = 0.40    ! (Chapter 3, p.34, IFS doc Cy31r1)
55   REAL(wp), PARAMETER ::   alpha_Q = 0.62    !
56
57   INTEGER , PARAMETER ::   nb_itt = 5        ! number of itterations
58
59   !!----------------------------------------------------------------------
60CONTAINS
61
62   SUBROUTINE turb_ecmwf( zt, zu, T_s, t_zt, q_s, q_zt, U_zu, &
63      &                      Cd, Ch, Ce, t_zu, q_zu, U_blk,      &
64      &                      Cdn, Chn, Cen,                      &
65      &                      Qsw, rad_lw, slp,                   &
66      &                      Tsk_b )
67      !!----------------------------------------------------------------------
68       !!                      ***  ROUTINE  turb_ecmwf  ***
69      !!
70      !! ** Purpose :   Computes turbulent transfert coefficients of surface
71      !!                fluxes according to IFS doc. (cycle 40)
72      !!                If relevant (zt /= zu), adjust temperature and humidity from height zt to zu
73      !!                Returns the effective bulk wind speed at zu to be used in the bulk formulas
74      !!
75      !!                Applies the cool-skin warm-layer correction of the SST to T_s
76      !!                if the net shortwave flux at the surface (Qsw), the downwelling longwave
77      !!                radiative fluxes at the surface (rad_lw), and the sea-leve pressure (slp)
78      !!                are provided as (optional) arguments!
79      !!
80      !! INPUT :
81      !! -------
82      !!    *  zt   : height for temperature and spec. hum. of air            [m]
83      !!    *  zu   : height for wind speed (usually 10m)                     [m]
84      !!    *  U_zu : scalar wind speed at zu                                 [m/s]
85      !!    *  t_zt : potential air temperature at zt                         [K]
86      !!    *  q_zt : specific humidity of air at zt                          [kg/kg]
87      !!
88      !! INPUT/OUTPUT:
89      !! -------------
90      !!    *  T_s  : always "bulk SST" as input                              [K]
91      !!              -> unchanged "bulk SST" as output if CSWL not used      [K]
92      !!              -> skin temperature as output if CSWL used              [K]
93      !!
94      !!    *  q_s  : SSQ aka saturation specific humidity at temp. T_s       [kg/kg]
95      !!              -> doesn't need to be given a value if skin temp computed (in case l_use_skin=True)
96      !!              -> MUST be given the correct value if not computing skint temp. (in case l_use_skin=False)
97      !!
98      !! OPTIONAL INPUT (will trigger l_use_skin=TRUE if present!):
99      !! ---------------
100      !!    *  Qsw    : net solar flux (after albedo) at the surface (>0)     [W/m^2]
101      !!    *  rad_lw : downwelling longwave radiation at the surface  (>0)   [W/m^2]
102      !!    *  slp    : sea-level pressure                                    [Pa]
103      !!    *  Tsk_b  : estimate of skin temperature at previous time-step    [K]
104      !!
105      !! OUTPUT :
106      !! --------
107      !!    *  Cd     : drag coefficient
108      !!    *  Ch     : sensible heat coefficient
109      !!    *  Ce     : evaporation coefficient
110      !!    *  t_zu   : pot. air temperature adjusted at wind height zu       [K]
111      !!    *  q_zu   : specific humidity of air        //                    [kg/kg]
112      !!    *  U_blk  : bulk wind speed at zu                                 [m/s]
113      !!
114      !!
115      !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/)
116      !!----------------------------------------------------------------------------------
117      REAL(wp), INTENT(in   )                     ::   zt       ! height for t_zt and q_zt                    [m]
118      REAL(wp), INTENT(in   )                     ::   zu       ! height for U_zu                             [m]
119      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) ::   T_s      ! sea surface temperature                [Kelvin]
120      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   t_zt     ! potential air temperature              [Kelvin]
121      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) ::   q_s      ! sea surface specific humidity           [kg/kg]
122      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_zt     ! specific air humidity at zt             [kg/kg]
123      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   U_zu     ! relative wind module at zu                [m/s]
124      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cd       ! transfer coefficient for momentum         (tau)
125      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ch       ! transfer coefficient for sensible heat (Q_sens)
126      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ce       ! transfert coefficient for evaporation   (Q_lat)
127      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   t_zu     ! pot. air temp. adjusted at zu               [K]
128      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   q_zu     ! spec. humidity adjusted at zu           [kg/kg]
129      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   U_blk    ! bulk wind speed at zu                     [m/s]
130      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients
131      !
132      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   Qsw      !             [W/m^2]
133      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   rad_lw   !             [W/m^2]
134      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   slp      !             [Pa]
135      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   Tsk_b    !             [Pa]
136      !
137      INTEGER :: j_itt
138      LOGICAL :: l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U
139      !
140      REAL(wp), DIMENSION(jpi,jpj) ::  &
141         &  u_star, t_star, q_star, &
142         &  dt_zu, dq_zu,    &
143         &  znu_a,           & !: Nu_air, Viscosity of air
144         &  Linv,            & !: 1/L (inverse of Monin Obukhov length...
145         &  z0, z0t, z0q
146      !
147      ! Cool skin:
148      LOGICAL :: l_use_skin = .FALSE.
149      REAL(wp), DIMENSION(jpi,jpj) :: zsst    ! to back up the initial bulk SST
150      !
151      REAL(wp), DIMENSION(jpi,jpj) ::   func_m, func_h
152      REAL(wp), DIMENSION(jpi,jpj) ::   ztmp0, ztmp1, ztmp2
153      !!----------------------------------------------------------------------------------
154
155      ! Cool skin ?
156      IF( PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp) ) THEN
157         l_use_skin = .TRUE.
158      END IF
159      IF (lwp) PRINT *, ' *** LOLO(sbcblk_algo_ecmwf.F90) => l_use_skin =', l_use_skin
160
161      ! Identical first gess as in COARE, with IFS parameter values though...
162      !
163      l_zt_equal_zu = .FALSE.
164      IF( ABS(zu - zt) < 0.01_wp )   l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision
165
166      !! Initialization for cool skin:
167      zsst   = T_s    ! save the bulk SST
168      IF( l_use_skin ) THEN
169         ! First guess for skin temperature:
170         IF( PRESENT(Tsk_b) ) THEN
171            T_s = Tsk_b
172         ELSE
173            T_s = T_s - 0.25     ! sst - 0.25
174         END IF
175         q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s
176      END IF
177
178      !! First guess of temperature and humidity at height zu:
179      t_zu = MAX( t_zt ,  180._wp )   ! who knows what's given on masked-continental regions...
180      q_zu = MAX( q_zt , 1.e-6_wp )   !               "
181
182      !! Pot. temp. difference (and we don't want it to be 0!)
183      dt_zu = t_zu - T_s ;   dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
184      dq_zu = q_zu - q_s ;   dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
185
186      znu_a = visc_air(t_zu) ! Air viscosity (m^2/s) at zt given from temperature in (K)
187
188      U_blk = SQRT(U_zu*U_zu + 0.5_wp*0.5_wp) ! initial guess for wind gustiness contribution
189
190      ztmp0   = LOG(    zu*10000._wp) ! optimization: 10000. == 1/z0 (with z0 first guess == 0.0001)
191      ztmp1   = LOG(10._wp*10000._wp) !       "                    "               "
192      u_star = 0.035_wp*U_blk*ztmp1/ztmp0       ! (u* = 0.035*Un10)
193
194      z0     = charn0*u_star*u_star/grav + 0.11_wp*znu_a/u_star
195      z0     = MIN(ABS(z0), 0.001_wp)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO
196      z0t    = 1._wp / ( 0.1_wp*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) )
197      z0t    = MIN(ABS(z0t), 0.001_wp)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO
198
199      ztmp2  = vkarmn/ztmp0
200      Cd     = ztmp2*ztmp2    ! first guess of Cd
201
202      ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd
203
204      ztmp2 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN)
205
206      !! First estimate of zeta_u, depending on the stability, ie sign of BRN (ztmp2):
207      ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 )
208      func_m = ztmp0*ztmp2 ! temporary array !!
209      func_h = (1._wp-ztmp1) * (func_m/(1._wp+ztmp2/(-zu/(zi0*0.004_wp*Beta0**3)))) & !  BRN < 0 ! temporary array !!! func_h == zeta_u
210         &  +     ztmp1   * (func_m*(1._wp + 27._wp/9._wp*ztmp2/func_m))              !  BRN > 0
211      !#LB: should make sure that the "func_m" of "27./9.*ztmp2/func_m" is "ztmp0*ztmp2" and not "ztmp0==vkarmn*vkarmn/LOG(zt/z0t)/Cd" !
212
213      !! First guess M-O stability dependent scaling params.(u*,t*,q*) to estimate z0 and z/L
214      ztmp0  = vkarmn/(LOG(zu/z0t) - psi_h_ecmwf(func_h))
215
216      u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0)  - psi_m_ecmwf(func_h))
217      t_star = dt_zu*ztmp0
218      q_star = dq_zu*ztmp0
219
220      ! What needs to be done if zt /= zu:
221      IF( .NOT. l_zt_equal_zu ) THEN
222         !! First update of values at zu (or zt for wind)
223         ztmp0 = psi_h_ecmwf(func_h) - psi_h_ecmwf(zt*func_h/zu)    ! zt*func_h/zu == zeta_t
224         ztmp1 = LOG(zt/zu) + ztmp0
225         t_zu = t_zt - t_star/vkarmn*ztmp1
226         q_zu = q_zt - q_star/vkarmn*ztmp1
227         q_zu = (0.5_wp + SIGN(0.5_wp,q_zu))*q_zu !Makes it impossible to have negative humidity :
228         !
229         dt_zu = t_zu - T_s  ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
230         dq_zu = q_zu - q_s  ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
231      END IF
232
233
234      !! => that was same first guess as in COARE...
235
236
237      !! First guess of inverse of Monin-Obukov length (1/L) :
238      Linv = One_on_L( t_zu, q_zu, u_star, t_star, q_star )
239
240      !! Functions such as  u* = U_blk*vkarmn/func_m
241      ztmp0 = zu*Linv
242      func_m = LOG(zu) - LOG(z0)  - psi_m_ecmwf(ztmp0) + psi_m_ecmwf( z0*Linv)
243      func_h = LOG(zu) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv)
244
245      !! ITERATION BLOCK
246      DO j_itt = 1, nb_itt
247
248         !! Bulk Richardson Number at z=zu (Eq. 3.25)
249         ztmp0 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN)
250
251         !! New estimate of the inverse of the Monin-Obukhon length (Linv == zeta/zu) :
252         Linv = ztmp0*func_m*func_m/func_h / zu     ! From Eq. 3.23, Chap.3.2.3, IFS doc - Cy40r1
253         !! Note: it is slightly different that the L we would get with the usual
254         Linv = SIGN( MIN(ABS(Linv),200._wp), Linv ) ! (prevent FPE from stupid values from masked region later on...) !#LOLO
255
256         !! Update func_m with new Linv:
257         func_m = LOG(zu) -LOG(z0) - psi_m_ecmwf(zu*Linv) + psi_m_ecmwf(z0*Linv) ! LB: should be "zu+z0" rather than "zu" alone, but z0 is tiny wrt zu!
258
259         !! Need to update roughness lengthes:
260         u_star = U_blk*vkarmn/func_m
261         ztmp2  = u_star*u_star
262         ztmp1  = znu_a/u_star
263         z0     = MIN( ABS( alpha_M*ztmp1 + charn0*ztmp2/grav ) , 0.001_wp)
264         z0t    = MIN( ABS( alpha_H*ztmp1                     ) , 0.001_wp)   ! eq.3.26, Chap.3, p.34, IFS doc - Cy31r1
265         z0q    = MIN( ABS( alpha_Q*ztmp1                     ) , 0.001_wp)
266
267         !! Update wind at 10m taking into acount convection-related wind gustiness:
268         !! => Chap. 3.2, IFS doc - Cy40r1, Eq.3.17 and Eq.3.18 + Eq.3.8
269         ! Only true when unstable (L<0) => when ztmp0 < 0 => - !!!
270         ztmp2 = ztmp2 * ( MAX(-zi0*Linv/vkarmn , 0._wp))**(2._wp/3._wp) ! => w*^2  (combining Eq. 3.8 and 3.18, hap.3, IFS doc - Cy31r1)
271         !! => equivalent using Beta=1 (gustiness parameter, 1.25 for COARE, also zi0=600 in COARE..)
272         U_blk = MAX( SQRT(U_zu*U_zu + ztmp2) , 0.2_wp )              ! eq.3.17, Chap.3, p.32, IFS doc - Cy31r1
273         ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0.
274
275
276         !! Need to update "theta" and "q" at zu in case they are given at different heights
277         !! as well the air-sea differences:
278         IF( .NOT. l_zt_equal_zu ) THEN
279            !! Arrays func_m and func_h are free for a while so using them as temporary arrays...
280            func_h = psi_h_ecmwf(zu*Linv) ! temporary array !!!
281            func_m = psi_h_ecmwf(zt*Linv) ! temporary array !!!
282
283            ztmp2  = psi_h_ecmwf(z0t*Linv)
284            ztmp0  = func_h - ztmp2
285            ztmp1  = vkarmn/(LOG(zu) - LOG(z0t) - ztmp0)
286            t_star = dt_zu*ztmp1
287            ztmp2  = ztmp0 - func_m + ztmp2
288            ztmp1  = LOG(zt/zu) + ztmp2
289            t_zu   = t_zt - t_star/vkarmn*ztmp1
290
291            ztmp2  = psi_h_ecmwf(z0q*Linv)
292            ztmp0  = func_h - ztmp2
293            ztmp1  = vkarmn/(LOG(zu) - LOG(z0q) - ztmp0)
294            q_star = dq_zu*ztmp1
295            ztmp2  = ztmp0 - func_m + ztmp2
296            ztmp1  = LOG(zt/zu) + ztmp2
297            q_zu   = q_zt - q_star/vkarmn*ztmp1
298         END IF
299
300         !! Updating because of updated z0 and z0t and new Linv...
301         ztmp0 = zu*Linv
302         func_m = log(zu) - LOG(z0 ) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf(z0 *Linv)
303         func_h = log(zu) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv)
304
305         !! SKIN related part
306         !! -----------------
307         IF( l_use_skin ) THEN
308            !! compute transfer coefficients at zu : lolo: verifier...
309            Ch = vkarmn*vkarmn/(func_m*func_h)
310            ztmp1 = LOG(zu) - LOG(z0q) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0q*Linv)   ! func_q
311            Ce = vkarmn*vkarmn/(func_m*ztmp1)
312            ! Non-Solar heat flux to the ocean:
313            ztmp1 = U_blk*MAX(rho_air(t_zu, q_zu, slp), 1._wp)     ! rho*U10
314            ztmp2 = T_s*T_s
315            ztmp1 = ztmp1 * ( Ce*L_vap(T_s)*(q_zu - q_s) + Ch*cp_air(q_zu)*(t_zu - T_s) ) & ! Total turb. heat flux
316               &       +      emiss_w*(rad_lw - stefan*ztmp2*ztmp2)                         ! Net longwave flux
317            !!         => "ztmp1" is the net non-solar surface heat flux !
318            !! Updating the values of the skin temperature T_s and q_s :
319            CALL CSWL_ECMWF( Qsw, ztmp1, u_star, zsst, T_s )
320            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
321         END IF
322
323         IF( (l_use_skin).OR.(.NOT. l_zt_equal_zu) ) THEN
324            dt_zu = t_zu - T_s ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
325            dq_zu = q_zu - q_s ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
326         END IF
327
328      END DO !DO j_itt = 1, nb_itt
329
330      Cd = vkarmn*vkarmn/(func_m*func_m)
331      Ch = vkarmn*vkarmn/(func_m*func_h)
332      ztmp2 = log(zu/z0q) - psi_h_ecmwf(zu*Linv) + psi_h_ecmwf(z0q*Linv)   ! func_q
333      Ce = vkarmn*vkarmn/(func_m*ztmp2)
334
335      Cdn = vkarmn*vkarmn / (log(zu/z0 )*log(zu/z0 ))
336      Chn = vkarmn*vkarmn / (log(zu/z0t)*log(zu/z0t))
337      Cen = vkarmn*vkarmn / (log(zu/z0q)*log(zu/z0q))
338
339   END SUBROUTINE TURB_ECMWF
340
341
342   FUNCTION psi_m_ecmwf( pzeta )
343      !!----------------------------------------------------------------------------------
344      !! Universal profile stability function for momentum
345      !!     ECMWF / as in IFS cy31r1 documentation, available online
346      !!     at ecmwf.int
347      !!
348      !! pzeta : stability paramenter, z/L where z is altitude measurement
349      !!         and L is M-O length
350      !!
351      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
352      !!----------------------------------------------------------------------------------
353      REAL(wp), DIMENSION(jpi,jpj) :: psi_m_ecmwf
354      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
355      !
356      INTEGER  ::   ji, jj    ! dummy loop indices
357      REAL(wp) :: zzeta, zx, ztmp, psi_unst, psi_stab, stab
358      !!----------------------------------------------------------------------------------
359      !
360      DO jj = 1, jpj
361         DO ji = 1, jpi
362            !
363            zzeta = MIN( pzeta(ji,jj) , 5._wp ) !! Very stable conditions (L positif and big!):
364            !
365            ! Unstable (Paulson 1970):
366            !   eq.3.20, Chap.3, p.33, IFS doc - Cy31r1
367            zx = SQRT(ABS(1._wp - 16._wp*zzeta))
368            ztmp = 1._wp + SQRT(zx)
369            ztmp = ztmp*ztmp
370            psi_unst = LOG( 0.125_wp*ztmp*(1._wp + zx) )   &
371               &       -2._wp*ATAN( SQRT(zx) ) + 0.5_wp*rpi
372            !
373            ! Unstable:
374            ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1
375            psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) &
376               &       - zzeta - 2._wp/3._wp*5._wp/0.35_wp
377            !
378            ! Combining:
379            stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1
380            !
381            psi_m_ecmwf(ji,jj) = (1._wp - stab) * psi_unst & ! (zzeta < 0) Unstable
382               &                +      stab  * psi_stab      ! (zzeta > 0) Stable
383            !
384         END DO
385      END DO
386      !
387   END FUNCTION psi_m_ecmwf
388
389
390   FUNCTION psi_h_ecmwf( pzeta )
391      !!----------------------------------------------------------------------------------
392      !! Universal profile stability function for temperature and humidity
393      !!     ECMWF / as in IFS cy31r1 documentation, available online
394      !!     at ecmwf.int
395      !!
396      !! pzeta : stability paramenter, z/L where z is altitude measurement
397      !!         and L is M-O length
398      !!
399      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
400      !!----------------------------------------------------------------------------------
401      REAL(wp), DIMENSION(jpi,jpj) :: psi_h_ecmwf
402      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
403      !
404      INTEGER  ::   ji, jj     ! dummy loop indices
405      REAL(wp) ::  zzeta, zx, psi_unst, psi_stab, stab
406      !!----------------------------------------------------------------------------------
407      !
408      DO jj = 1, jpj
409         DO ji = 1, jpi
410            !
411            zzeta = MIN(pzeta(ji,jj) , 5._wp)   ! Very stable conditions (L positif and big!):
412            !
413            zx  = ABS(1._wp - 16._wp*zzeta)**.25        ! this is actually (1/phi_m)**2  !!!
414            !                                     ! eq.3.19, Chap.3, p.33, IFS doc - Cy31r1
415            ! Unstable (Paulson 1970) :
416            psi_unst = 2._wp*LOG(0.5_wp*(1._wp + zx*zx))   ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1
417            !
418            ! Stable:
419            psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) & ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1
420               &       - ABS(1._wp + 2._wp/3._wp*zzeta)**1.5_wp - 2._wp/3._wp*5._wp/0.35_wp + 1._wp
421            ! LB: added ABS() to avoid NaN values when unstable, which contaminates the unstable solution...
422            !
423            stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1
424            !
425            !
426            psi_h_ecmwf(ji,jj) = (1._wp - stab) * psi_unst &   ! (zzeta < 0) Unstable
427               &                +    stab    * psi_stab        ! (zzeta > 0) Stable
428            !
429         END DO
430      END DO
431      !
432   END FUNCTION psi_h_ecmwf
433
434
435
436
437
438   !!======================================================================
439END MODULE sbcblk_algo_ecmwf
Note: See TracBrowser for help on using the repository browser.