New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
sbcblk_algo_ecmwf.F90 in NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC – NEMO

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

Last change on this file since 11615 was 11615, checked in by laurent, 4 years ago

LB: new cool-skin and warm-layer parameterizations for ECMWF and COARE3.6, COARE3.0 uses same CSWL param as for COARE3.6.

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