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

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