source: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/SBC/sbcblk_algo_ecmwf.F90 @ 12991

Last change on this file since 12991 was 12991, checked in by emanuelaclementi, 5 months ago

Included wave-current processes following Couvelard et al 2019 and Gurvan code -ticket #2155

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