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

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

Improving syntax consistency

  • 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   !!----------------------------------------------------------------------
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 :: 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 ecmwf_init(l_use_cs, l_use_wl)
64      !!---------------------------------------------------------------------
65      !!                  ***  FUNCTION 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( ' 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      END IF
83      IF ( l_use_cs ) THEN
84         ierr = 0
85         ALLOCATE ( dT_cs(jpi,jpj), STAT=ierr )
86         IF( ierr > 0 ) CALL ctl_stop( ' ECMWF_INIT => allocation of dT_cs failed!' )
87         dT_cs(:,:) = -0.25_wp  ! First guess of skin correction
88      END IF
89   END SUBROUTINE 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) ::  &
185         &  u_star, t_star, q_star, &
186         &  dt_zu, dq_zu,    &
187         &  znu_a,           & !: Nu_air, Viscosity of air
188         &  Linv,            & !: 1/L (inverse of Monin Obukhov length...
189         &  z0, z0t, z0q
190      !
191      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
192         &                zsst,   &  ! to back up the initial bulk SST
193         &                pdTc,   &  ! SST increment "dT" for cool-skin correction           [K]
194         &                pdTw       ! SST increment "dT" for warm layer correction          [K]
195      !
196      REAL(wp), DIMENSION(jpi,jpj) ::   func_m, func_h
197      REAL(wp), DIMENSION(jpi,jpj) ::   ztmp0, ztmp1, ztmp2
198      CHARACTER(len=40), PARAMETER :: crtnm = 'turb_ecmwf@sbcblk_algo_ecmwf.F90'
199      !!----------------------------------------------------------------------------------
200
201      IF ( kt == nit000 ) CALL ECMWF_INIT(l_use_cs, l_use_wl)
202
203      l_zt_equal_zu = .FALSE.
204      IF( ABS(zu - zt) < 0.01_wp )   l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision
205
206      !! Initializations for cool skin and warm layer:
207      IF ( l_use_cs .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 cool-skin param!' )
209
210      IF ( l_use_wl .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) &
211         &   CALL ctl_stop( '['//TRIM(crtnm)//'] => ' , 'you need to provide Qsw, rad_lw & slp to use warm-layer param!' )
212
213      IF ( l_use_cs .OR. l_use_wl ) THEN
214         ALLOCATE ( zsst(jpi,jpj) )
215         zsst = T_s ! backing up the bulk SST
216         IF( l_use_cs ) T_s = T_s - 0.25_wp   ! First guess of correction
217         q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s
218      END IF
219
220
221      ! Identical first gess as in COARE, with IFS parameter values though...
222      !
223      !! First guess of temperature and humidity at height zu:
224      t_zu = MAX( t_zt ,  180._wp )   ! who knows what's given on masked-continental regions...
225      q_zu = MAX( q_zt , 1.e-6_wp )   !               "
226
227      !! Pot. temp. difference (and we don't want it to be 0!)
228      dt_zu = t_zu - T_s ;   dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
229      dq_zu = q_zu - q_s ;   dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
230
231      znu_a = visc_air(t_zu) ! Air viscosity (m^2/s) at zt given from temperature in (K)
232
233      U_blk = SQRT(U_zu*U_zu + 0.5_wp*0.5_wp) ! initial guess for wind gustiness contribution
234
235      ztmp0   = LOG(    zu*10000._wp) ! optimization: 10000. == 1/z0 (with z0 first guess == 0.0001)
236      ztmp1   = LOG(10._wp*10000._wp) !       "                    "               "
237      u_star = 0.035_wp*U_blk*ztmp1/ztmp0       ! (u* = 0.035*Un10)
238
239      z0     = charn0*u_star*u_star/grav + 0.11_wp*znu_a/u_star
240      z0     = MIN( MAX(ABS(z0), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on)
241
242      z0t    = 1._wp / ( 0.1_wp*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) )
243      z0t    = MIN( MAX(ABS(z0t), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on)
244
245      Cd     = (vkarmn/ztmp0)**2    ! first guess of Cd
246
247      ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd
248
249      ztmp2 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN)
250
251      !! First estimate of zeta_u, depending on the stability, ie sign of BRN (ztmp2):
252      ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 )
253      func_m = ztmp0*ztmp2 ! temporary array !!
254      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
255         &  +     ztmp1   * (func_m*(1._wp + 27._wp/9._wp*ztmp2/func_m))              !  BRN > 0
256      !#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" !
257
258      !! First guess M-O stability dependent scaling params.(u*,t*,q*) to estimate z0 and z/L
259      ztmp0  = vkarmn/(LOG(zu/z0t) - psi_h_ecmwf(func_h))
260
261      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)
262      t_star = dt_zu*ztmp0
263      q_star = dq_zu*ztmp0
264
265      ! What needs to be done if zt /= zu:
266      IF( .NOT. l_zt_equal_zu ) THEN
267         !! First update of values at zu (or zt for wind)
268         ztmp0 = psi_h_ecmwf(func_h) - psi_h_ecmwf(zt*func_h/zu)    ! zt*func_h/zu == zeta_t
269         ztmp1 = LOG(zt/zu) + ztmp0
270         t_zu = t_zt - t_star/vkarmn*ztmp1
271         q_zu = q_zt - q_star/vkarmn*ztmp1
272         q_zu = (0.5_wp + SIGN(0.5_wp,q_zu))*q_zu !Makes it impossible to have negative humidity :
273         !
274         dt_zu = t_zu - T_s  ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
275         dq_zu = q_zu - q_s  ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
276      END IF
277
278
279      !! => that was same first guess as in COARE...
280
281
282      !! First guess of inverse of Monin-Obukov length (1/L) :
283      Linv = One_on_L( t_zu, q_zu, u_star, t_star, q_star )
284
285      !! Functions such as  u* = U_blk*vkarmn/func_m
286      ztmp0 = zu*Linv
287      func_m = LOG(zu) - LOG(z0)  - psi_m_ecmwf(ztmp0) + psi_m_ecmwf( z0*Linv)
288      func_h = LOG(zu) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv)
289
290      !! ITERATION BLOCK
291      DO j_itt = 1, nb_itt
292
293         !! Bulk Richardson Number at z=zu (Eq. 3.25)
294         ztmp0 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN)
295
296         !! New estimate of the inverse of the Monin-Obukhon length (Linv == zeta/zu) :
297         Linv = ztmp0*func_m*func_m/func_h / zu     ! From Eq. 3.23, Chap.3.2.3, IFS doc - Cy40r1
298         !! Note: it is slightly different that the L we would get with the usual
299         Linv = SIGN( MIN(ABS(Linv),200._wp), Linv ) ! (prevent FPE from stupid values from masked region later on...) !#LOLO
300
301         !! Update func_m with new Linv:
302         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!
303
304         !! Need to update roughness lengthes:
305         u_star = U_blk*vkarmn/func_m
306         ztmp2  = u_star*u_star
307         ztmp1  = znu_a/u_star
308         z0     = MIN( ABS( alpha_M*ztmp1 + charn0*ztmp2/grav ) , 0.001_wp)
309         z0t    = MIN( ABS( alpha_H*ztmp1                     ) , 0.001_wp)   ! eq.3.26, Chap.3, p.34, IFS doc - Cy31r1
310         z0q    = MIN( ABS( alpha_Q*ztmp1                     ) , 0.001_wp)
311
312         !! 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)
313         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)
314         !!   ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before zi0
315         U_blk = MAX(SQRT(U_zu*U_zu + ztmp2), 0.2_wp)        ! include gustiness in bulk wind speed
316         ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0.
317
318
319         !! Need to update "theta" and "q" at zu in case they are given at different heights
320         !! as well the air-sea differences:
321         IF( .NOT. l_zt_equal_zu ) THEN
322            !! Arrays func_m and func_h are free for a while so using them as temporary arrays...
323            func_h = psi_h_ecmwf(zu*Linv) ! temporary array !!!
324            func_m = psi_h_ecmwf(zt*Linv) ! temporary array !!!
325
326            ztmp2  = psi_h_ecmwf(z0t*Linv)
327            ztmp0  = func_h - ztmp2
328            ztmp1  = vkarmn/(LOG(zu) - LOG(z0t) - ztmp0)
329            t_star = dt_zu*ztmp1
330            ztmp2  = ztmp0 - func_m + ztmp2
331            ztmp1  = LOG(zt/zu) + ztmp2
332            t_zu   = t_zt - t_star/vkarmn*ztmp1
333
334            ztmp2  = psi_h_ecmwf(z0q*Linv)
335            ztmp0  = func_h - ztmp2
336            ztmp1  = vkarmn/(LOG(zu) - LOG(z0q) - ztmp0)
337            q_star = dq_zu*ztmp1
338            ztmp2  = ztmp0 - func_m + ztmp2
339            ztmp1  = LOG(zt/zu) + ztmp2
340            q_zu   = q_zt - q_star/vkarmn*ztmp1
341         END IF
342
343         !! Updating because of updated z0 and z0t and new Linv...
344         ztmp0 = zu*Linv
345         func_m = log(zu) - LOG(z0 ) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf(z0 *Linv)
346         func_h = log(zu) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv)
347
348
349         IF( l_use_cs ) THEN
350            !! Cool-skin contribution
351
352            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, &
353               &                   ztmp1, ztmp0,  Qlat=ztmp2)  ! Qnsol -> ztmp1 / Tau -> ztmp0
354
355            CALL CS_ECMWF( Qsw, ztmp1, u_star, zsst )  ! Qnsol -> ztmp1
356
357            T_s(:,:) = zsst(:,:) + dT_cs(:,:)*tmask(:,:,1)
358            IF( l_use_wl ) T_s(:,:) = T_s(:,:) + dT_wl(:,:)*tmask(:,:,1)
359            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:))
360
361         END IF
362
363         IF( l_use_wl ) THEN
364            !! Warm-layer contribution
365            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, &
366               &                   ztmp1, ztmp2)  ! Qnsol -> ztmp1 / Tau -> ztmp2
367            CALL WL_ECMWF( Qsw, ztmp1, u_star, zsst )
368            !! Updating T_s and q_s !!!
369            T_s(:,:) = zsst(:,:) + dT_wl(:,:)*tmask(:,:,1)
370            IF( l_use_cs ) T_s(:,:) = T_s(:,:) + dT_cs(:,:)*tmask(:,:,1)
371            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:))
372         END IF
373
374         IF( l_use_cs .OR. l_use_wl .OR. (.NOT. l_zt_equal_zu) ) THEN
375            dt_zu = t_zu - T_s ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
376            dq_zu = q_zu - q_s ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
377         END IF
378
379      END DO !DO j_itt = 1, nb_itt
380
381      Cd = vkarmn*vkarmn/(func_m*func_m)
382      Ch = vkarmn*vkarmn/(func_m*func_h)
383      ztmp2 = log(zu/z0q) - psi_h_ecmwf(zu*Linv) + psi_h_ecmwf(z0q*Linv)   ! func_q
384      Ce = vkarmn*vkarmn/(func_m*ztmp2)
385
386      Cdn = vkarmn*vkarmn / (log(zu/z0 )*log(zu/z0 ))
387      Chn = vkarmn*vkarmn / (log(zu/z0t)*log(zu/z0t))
388      Cen = vkarmn*vkarmn / (log(zu/z0q)*log(zu/z0q))
389
390      IF ( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = dT_cs
391      IF ( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = dT_wl
392      IF ( l_use_wl .AND. PRESENT(pHz_wl) ) pHz_wl = Hz_wl
393
394      IF ( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst )
395
396   END SUBROUTINE turb_ecmwf
397
398
399   FUNCTION psi_m_ecmwf( pzeta )
400      !!----------------------------------------------------------------------------------
401      !! Universal profile stability function for momentum
402      !!     ECMWF / as in IFS cy31r1 documentation, available online
403      !!     at ecmwf.int
404      !!
405      !! pzeta : stability paramenter, z/L where z is altitude measurement
406      !!         and L is M-O length
407      !!
408      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
409      !!----------------------------------------------------------------------------------
410      REAL(wp), DIMENSION(jpi,jpj) :: psi_m_ecmwf
411      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
412      !
413      INTEGER  ::   ji, jj    ! dummy loop indices
414      REAL(wp) :: zzeta, zx, ztmp, psi_unst, psi_stab, stab
415      !!----------------------------------------------------------------------------------
416      DO jj = 1, jpj
417         DO ji = 1, jpi
418            !
419            zzeta = MIN( pzeta(ji,jj) , 5._wp ) !! Very stable conditions (L positif and big!):
420            !
421            ! Unstable (Paulson 1970):
422            !   eq.3.20, Chap.3, p.33, IFS doc - Cy31r1
423            zx = SQRT(ABS(1._wp - 16._wp*zzeta))
424            ztmp = 1._wp + SQRT(zx)
425            ztmp = ztmp*ztmp
426            psi_unst = LOG( 0.125_wp*ztmp*(1._wp + zx) )   &
427               &       -2._wp*ATAN( SQRT(zx) ) + 0.5_wp*rpi
428            !
429            ! Unstable:
430            ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1
431            psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) &
432               &       - zzeta - 2._wp/3._wp*5._wp/0.35_wp
433            !
434            ! Combining:
435            stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1
436            !
437            psi_m_ecmwf(ji,jj) = (1._wp - stab) * psi_unst & ! (zzeta < 0) Unstable
438               &                +      stab  * psi_stab      ! (zzeta > 0) Stable
439            !
440         END DO
441      END DO
442   END FUNCTION psi_m_ecmwf
443
444
445   FUNCTION psi_h_ecmwf( pzeta )
446      !!----------------------------------------------------------------------------------
447      !! Universal profile stability function for temperature and humidity
448      !!     ECMWF / as in IFS cy31r1 documentation, available online
449      !!     at ecmwf.int
450      !!
451      !! pzeta : stability paramenter, z/L where z is altitude measurement
452      !!         and L is M-O length
453      !!
454      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
455      !!----------------------------------------------------------------------------------
456      REAL(wp), DIMENSION(jpi,jpj) :: psi_h_ecmwf
457      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
458      !
459      INTEGER  ::   ji, jj     ! dummy loop indices
460      REAL(wp) ::  zzeta, zx, psi_unst, psi_stab, stab
461      !!----------------------------------------------------------------------------------
462      !
463      DO jj = 1, jpj
464         DO ji = 1, jpi
465            !
466            zzeta = MIN(pzeta(ji,jj) , 5._wp)   ! Very stable conditions (L positif and big!):
467            !
468            zx  = ABS(1._wp - 16._wp*zzeta)**.25        ! this is actually (1/phi_m)**2  !!!
469            !                                     ! eq.3.19, Chap.3, p.33, IFS doc - Cy31r1
470            ! Unstable (Paulson 1970) :
471            psi_unst = 2._wp*LOG(0.5_wp*(1._wp + zx*zx))   ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1
472            !
473            ! Stable:
474            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
475               &       - ABS(1._wp + 2._wp/3._wp*zzeta)**1.5_wp - 2._wp/3._wp*5._wp/0.35_wp + 1._wp
476            ! LB: added ABS() to avoid NaN values when unstable, which contaminates the unstable solution...
477            !
478            stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1
479            !
480            !
481            psi_h_ecmwf(ji,jj) = (1._wp - stab) * psi_unst &   ! (zzeta < 0) Unstable
482               &                +    stab    * psi_stab        ! (zzeta > 0) Stable
483            !
484         END DO
485      END DO
486   END FUNCTION psi_h_ecmwf
487
488
489   !!======================================================================
490END MODULE sbcblk_algo_ecmwf
Note: See TracBrowser for help on using the repository browser.