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

Last change on this file was 11963, checked in by laurent, 11 months ago

More accurate comments/info, better syntax, simplifications, etc

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