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

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

LB: skin temperature now used by ECMWF bulk algorithm, and can be xios-ed ; new module "SBC/sbcblk_skin.F90" ; all physical constants defined in SBC now moved to "DOM/phycst.F90"; all air-properties and other ABL functions gathered in a new module: "SBC/sbcblk_phymbl.F90".

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