source: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ecmwf.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.

  • 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 (usually 10m)                     [m]
83      !!    *  U_zu : scalar wind speed at zu                                 [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 zu                                 [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 speed at zu                     [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      U_blk = SQRT(U_zu*U_zu + 0.5_wp*0.5_wp) ! initial guess for wind gustiness contribution
188
189      ztmp0   = LOG(    zu*10000._wp) ! optimization: 10000. == 1/z0 (with z0 first guess == 0.0001)
190      ztmp1   = LOG(10._wp*10000._wp) !       "                    "               "
191      u_star = 0.035_wp*U_blk*ztmp1/ztmp0       ! (u* = 0.035*Un10)
192
193      z0     = charn0*u_star*u_star/grav + 0.11_wp*znu_a/u_star
194      z0     = MIN(ABS(z0), 0.001_wp)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO
195      z0t    = 1._wp / ( 0.1_wp*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) )
196      z0t    = MIN(ABS(z0t), 0.001_wp)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO
197
198      ztmp2  = vkarmn/ztmp0
199      Cd     = ztmp2*ztmp2    ! first guess of Cd
200
201      ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd
202
203      ztmp2 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN)
204
205      !! First estimate of zeta_u, depending on the stability, ie sign of BRN (ztmp2):
206      ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 )
207      func_m = ztmp0*ztmp2 ! temporary array !!
208      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
209         &  +     ztmp1   * (func_m*(1._wp + 27._wp/9._wp*ztmp2/func_m))              !  BRN > 0
210      !#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" !
211
212      !! First guess M-O stability dependent scaling params.(u*,t*,q*) to estimate z0 and z/L
213      ztmp0  = vkarmn/(LOG(zu/z0t) - psi_h_ecmwf(func_h))
214
215      u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0)  - psi_m_ecmwf(func_h))
216      t_star = dt_zu*ztmp0
217      q_star = dq_zu*ztmp0
218
219      ! What needs to be done if zt /= zu:
220      IF( .NOT. l_zt_equal_zu ) THEN
221         !! First update of values at zu (or zt for wind)
222         ztmp0 = psi_h_ecmwf(func_h) - psi_h_ecmwf(zt*func_h/zu)    ! zt*func_h/zu == zeta_t
223         ztmp1 = LOG(zt/zu) + ztmp0
224         t_zu = t_zt - t_star/vkarmn*ztmp1
225         q_zu = q_zt - q_star/vkarmn*ztmp1
226         q_zu = (0.5_wp + SIGN(0.5_wp,q_zu))*q_zu !Makes it impossible to have negative humidity :
227         !
228         dt_zu = t_zu - T_s  ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
229         dq_zu = q_zu - q_s  ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
230      END IF
231
232
233      !! => that was same first guess as in COARE...
234
235
236      !! First guess of inverse of Monin-Obukov length (1/L) :
237      Linv = One_on_L( t_zu, q_zu, u_star, t_star, q_star )
238
239      !! Functions such as  u* = U_blk*vkarmn/func_m
240      ztmp0 = zu*Linv
241      func_m = LOG(zu) - LOG(z0)  - psi_m_ecmwf(ztmp0) + psi_m_ecmwf( z0*Linv)
242      func_h = LOG(zu) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv)
243
244      !! ITERATION BLOCK
245      DO j_itt = 1, nb_itt
246
247         !! Bulk Richardson Number at z=zu (Eq. 3.25)
248         ztmp0 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN)
249
250         !! New estimate of the inverse of the Monin-Obukhon length (Linv == zeta/zu) :
251         Linv = ztmp0*func_m*func_m/func_h / zu     ! From Eq. 3.23, Chap.3.2.3, IFS doc - Cy40r1
252         !! Note: it is slightly different that the L we would get with the usual
253         Linv = SIGN( MIN(ABS(Linv),200._wp), Linv ) ! (prevent FPE from stupid values from masked region later on...) !#LOLO
254
255         !! Update func_m with new Linv:
256         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!
257
258         !! Need to update roughness lengthes:
259         u_star = U_blk*vkarmn/func_m
260         ztmp2  = u_star*u_star
261         ztmp1  = znu_a/u_star
262         z0     = MIN( ABS( alpha_M*ztmp1 + charn0*ztmp2/grav ) , 0.001_wp)
263         z0t    = MIN( ABS( alpha_H*ztmp1                     ) , 0.001_wp)   ! eq.3.26, Chap.3, p.34, IFS doc - Cy31r1
264         z0q    = MIN( ABS( alpha_Q*ztmp1                     ) , 0.001_wp)
265
266         !! Update wind at 10m taking into acount convection-related wind gustiness:
267         !! => Chap. 3.2, IFS doc - Cy40r1, Eq.3.17 and Eq.3.18 + Eq.3.8
268         ! Only true when unstable (L<0) => when ztmp0 < 0 => - !!!
269         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)
270         !! => equivalent using Beta=1 (gustiness parameter, 1.25 for COARE, also zi0=600 in COARE..)
271         U_blk = MAX( SQRT(U_zu*U_zu + ztmp2) , 0.2_wp )              ! eq.3.17, Chap.3, p.32, IFS doc - Cy31r1
272         ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0.
273
274
275         !! Need to update "theta" and "q" at zu in case they are given at different heights
276         !! as well the air-sea differences:
277         IF( .NOT. l_zt_equal_zu ) THEN
278            !! Arrays func_m and func_h are free for a while so using them as temporary arrays...
279            func_h = psi_h_ecmwf(zu*Linv) ! temporary array !!!
280            func_m = psi_h_ecmwf(zt*Linv) ! temporary array !!!
281
282            ztmp2  = psi_h_ecmwf(z0t*Linv)
283            ztmp0  = func_h - ztmp2
284            ztmp1  = vkarmn/(LOG(zu) - LOG(z0t) - ztmp0)
285            t_star = dt_zu*ztmp1
286            ztmp2  = ztmp0 - func_m + ztmp2
287            ztmp1  = LOG(zt/zu) + ztmp2
288            t_zu   = t_zt - t_star/vkarmn*ztmp1
289
290            ztmp2  = psi_h_ecmwf(z0q*Linv)
291            ztmp0  = func_h - ztmp2
292            ztmp1  = vkarmn/(LOG(zu) - LOG(z0q) - ztmp0)
293            q_star = dq_zu*ztmp1
294            ztmp2  = ztmp0 - func_m + ztmp2
295            ztmp1  = LOG(zt/zu) + ztmp2
296            q_zu   = q_zt - q_star/vkarmn*ztmp1
297         END IF
298
299         !! Updating because of updated z0 and z0t and new Linv...
300         ztmp0 = zu*Linv
301         func_m = log(zu) - LOG(z0 ) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf(z0 *Linv)
302         func_h = log(zu) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv)
303
304         !! SKIN related part
305         !! -----------------
306         IF( l_use_skin ) THEN
307            !! compute transfer coefficients at zu : lolo: verifier...
308            Ch = vkarmn*vkarmn/(func_m*func_h)
309            ztmp1 = LOG(zu) - LOG(z0q) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0q*Linv)   ! func_q
310            Ce = vkarmn*vkarmn/(func_m*ztmp1)
311            ! Non-Solar heat flux to the ocean:
312            ztmp1 = U_blk*MAX(rho_air(t_zu, q_zu, slp), 1._wp)     ! rho*U10
313            ztmp2 = T_s*T_s
314            ztmp1 = ztmp1 * ( Ce*L_vap(T_s)*(q_zu - q_s) + Ch*cp_air(q_zu)*(t_zu - T_s) ) & ! Total turb. heat flux
315               &       +      emiss_w*(rad_lw - stefan*ztmp2*ztmp2)                         ! Net longwave flux
316            !!         => "ztmp1" is the net non-solar surface heat flux !
317            !! Updating the values of the skin temperature T_s and q_s :
318            CALL CSWL_ECMWF( Qsw, ztmp1, u_star, zsst, T_s )
319            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
320         END IF
321
322         IF( (l_use_skin).OR.(.NOT. l_zt_equal_zu) ) THEN
323            dt_zu = t_zu - T_s ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
324            dq_zu = q_zu - q_s ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
325         END IF
326
327      END DO !DO j_itt = 1, nb_itt
328
329      Cd = vkarmn*vkarmn/(func_m*func_m)
330      Ch = vkarmn*vkarmn/(func_m*func_h)
331      ztmp2 = log(zu/z0q) - psi_h_ecmwf(zu*Linv) + psi_h_ecmwf(z0q*Linv)   ! func_q
332      Ce = vkarmn*vkarmn/(func_m*ztmp2)
333
334      Cdn = vkarmn*vkarmn / (log(zu/z0 )*log(zu/z0 ))
335      Chn = vkarmn*vkarmn / (log(zu/z0t)*log(zu/z0t))
336      Cen = vkarmn*vkarmn / (log(zu/z0q)*log(zu/z0q))
337
338   END SUBROUTINE TURB_ECMWF
339
340
341   FUNCTION psi_m_ecmwf( pzeta )
342      !!----------------------------------------------------------------------------------
343      !! Universal profile stability function for momentum
344      !!     ECMWF / as in IFS cy31r1 documentation, available online
345      !!     at ecmwf.int
346      !!
347      !! pzeta : stability paramenter, z/L where z is altitude measurement
348      !!         and L is M-O length
349      !!
350      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
351      !!----------------------------------------------------------------------------------
352      REAL(wp), DIMENSION(jpi,jpj) :: psi_m_ecmwf
353      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
354      !
355      INTEGER  ::   ji, jj    ! dummy loop indices
356      REAL(wp) :: zzeta, zx, ztmp, psi_unst, psi_stab, stab
357      !!----------------------------------------------------------------------------------
358      !
359      DO jj = 1, jpj
360         DO ji = 1, jpi
361            !
362            zzeta = MIN( pzeta(ji,jj) , 5._wp ) !! Very stable conditions (L positif and big!):
363            !
364            ! Unstable (Paulson 1970):
365            !   eq.3.20, Chap.3, p.33, IFS doc - Cy31r1
366            zx = SQRT(ABS(1._wp - 16._wp*zzeta))
367            ztmp = 1._wp + SQRT(zx)
368            ztmp = ztmp*ztmp
369            psi_unst = LOG( 0.125_wp*ztmp*(1._wp + zx) )   &
370               &       -2._wp*ATAN( SQRT(zx) ) + 0.5_wp*rpi
371            !
372            ! Unstable:
373            ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1
374            psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) &
375               &       - zzeta - 2._wp/3._wp*5._wp/0.35_wp
376            !
377            ! Combining:
378            stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1
379            !
380            psi_m_ecmwf(ji,jj) = (1._wp - stab) * psi_unst & ! (zzeta < 0) Unstable
381               &                +      stab  * psi_stab      ! (zzeta > 0) Stable
382            !
383         END DO
384      END DO
385      !
386   END FUNCTION psi_m_ecmwf
387
388
389   FUNCTION psi_h_ecmwf( pzeta )
390      !!----------------------------------------------------------------------------------
391      !! Universal profile stability function for temperature and humidity
392      !!     ECMWF / as in IFS cy31r1 documentation, available online
393      !!     at ecmwf.int
394      !!
395      !! pzeta : stability paramenter, z/L where z is altitude measurement
396      !!         and L is M-O length
397      !!
398      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
399      !!----------------------------------------------------------------------------------
400      REAL(wp), DIMENSION(jpi,jpj) :: psi_h_ecmwf
401      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
402      !
403      INTEGER  ::   ji, jj     ! dummy loop indices
404      REAL(wp) ::  zzeta, zx, psi_unst, psi_stab, stab
405      !!----------------------------------------------------------------------------------
406      !
407      DO jj = 1, jpj
408         DO ji = 1, jpi
409            !
410            zzeta = MIN(pzeta(ji,jj) , 5._wp)   ! Very stable conditions (L positif and big!):
411            !
412            zx  = ABS(1._wp - 16._wp*zzeta)**.25        ! this is actually (1/phi_m)**2  !!!
413            !                                     ! eq.3.19, Chap.3, p.33, IFS doc - Cy31r1
414            ! Unstable (Paulson 1970) :
415            psi_unst = 2._wp*LOG(0.5_wp*(1._wp + zx*zx))   ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1
416            !
417            ! Stable:
418            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
419               &       - ABS(1._wp + 2._wp/3._wp*zzeta)**1.5_wp - 2._wp/3._wp*5._wp/0.35_wp + 1._wp
420            ! LB: added ABS() to avoid NaN values when unstable, which contaminates the unstable solution...
421            !
422            stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1
423            !
424            !
425            psi_h_ecmwf(ji,jj) = (1._wp - stab) * psi_unst &   ! (zzeta < 0) Unstable
426               &                +    stab    * psi_stab        ! (zzeta > 0) Stable
427            !
428         END DO
429      END DO
430      !
431   END FUNCTION psi_h_ecmwf
432
433
434
435
436
437   !!======================================================================
438END MODULE sbcblk_algo_ecmwf
Note: See TracBrowser for help on using the repository browser.