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 @ 11772

Last change on this file since 11772 was 11772, checked in by laurent, 5 years ago

LB: solid updates+improvements of cool-skin/warm-layer capabilty of COARE and ECMWF bulk algorithms!

  • 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 ) 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      !!
84      IF ( l_use_cs ) THEN
85         ierr = 0
86         ALLOCATE ( dT_cs(jpi,jpj), STAT=ierr )
87         !IF( ierr > 0 ) STOP ' ECMWF_INIT => allocation of dT_cs failed!'
88         dT_cs(:,:) = -0.25_wp  ! First guess of skin correction
89      END IF
90   END SUBROUTINE ecmwf_init
91
92
93
94   SUBROUTINE turb_ecmwf( kt, zt, zu, T_s, t_zt, q_s, q_zt, U_zu, l_use_cs, l_use_wl, &
95      &                      Cd, Ch, Ce, t_zu, q_zu, U_blk,                           &
96      &                      Cdn, Chn, Cen,                                           &
97      &                      Qsw, rad_lw, slp, pdT_cs,                                & ! optionals for cool-skin (and warm-layer)
98      &                      pdT_wl, pHz_wl )                                                 ! optionals for warm-layer only
99      !!----------------------------------------------------------------------
100      !!                      ***  ROUTINE  turb_ecmwf  ***
101      !!
102      !! ** Purpose :   Computes turbulent transfert coefficients of surface
103      !!                fluxes according to IFS doc. (cycle 45r1)
104      !!                If relevant (zt /= zu), adjust temperature and humidity from height zt to zu
105      !!                Returns the effective bulk wind speed at zu to be used in the bulk formulas
106      !!
107      !!                Applies the cool-skin warm-layer correction of the SST to T_s
108      !!                if the net shortwave flux at the surface (Qsw), the downwelling longwave
109      !!                radiative fluxes at the surface (rad_lw), and the sea-leve pressure (slp)
110      !!                are provided as (optional) arguments!
111      !!
112      !! INPUT :
113      !! -------
114      !!    *  kt   : current time step (starts at 1)
115      !!    *  zt   : height for temperature and spec. hum. of air            [m]
116      !!    *  zu   : height for wind speed (usually 10m)                     [m]
117      !!    *  t_zt : potential air temperature at zt                         [K]
118      !!    *  q_zt : specific humidity of air at zt                          [kg/kg]
119      !!    *  U_zu : scalar wind speed at zu                                 [m/s]
120      !!    * l_use_cs : use the cool-skin parameterization
121      !!    * l_use_wl : use the warm-layer parameterization
122      !!
123      !! INPUT/OUTPUT:
124      !! -------------
125      !!    *  T_s  : always "bulk SST" as input                              [K]
126      !!              -> unchanged "bulk SST" as output if CSWL not used      [K]
127      !!              -> skin temperature as output if CSWL used              [K]
128      !!
129      !!    *  q_s  : SSQ aka saturation specific humidity at temp. T_s       [kg/kg]
130      !!              -> doesn't need to be given a value if skin temp computed (in case l_use_cs=True or l_use_wl=True)
131      !!              -> MUST be given the correct value if not computing skint temp. (in case l_use_cs=False or l_use_wl=False)
132      !!
133      !! OPTIONAL INPUT:
134      !! ---------------
135      !!    *  Qsw    : net solar flux (after albedo) at the surface (>0)     [W/m^2]
136      !!    *  rad_lw : downwelling longwave radiation at the surface  (>0)   [W/m^2]
137      !!    *  slp    : sea-level pressure                                    [Pa]
138      !!
139      !! OPTIONAL OUTPUT:
140      !! ----------------
141      !!    * pdT_cs  : SST increment "dT" for cool-skin correction           [K]
142      !!    * pdT_wl  : SST increment "dT" for warm-layer correction          [K]
143      !!    * pHz_wl  : thickness of warm-layer                               [m]
144      !!
145      !! OUTPUT :
146      !! --------
147      !!    *  Cd     : drag coefficient
148      !!    *  Ch     : sensible heat coefficient
149      !!    *  Ce     : evaporation coefficient
150      !!    *  t_zu   : pot. air temperature adjusted at wind height zu       [K]
151      !!    *  q_zu   : specific humidity of air        //                    [kg/kg]
152      !!    *  U_blk  : bulk wind speed at zu                                 [m/s]
153      !!
154      !!
155      !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/)
156      !!----------------------------------------------------------------------------------
157      INTEGER,  INTENT(in   )                     ::   kt       ! current time step
158      REAL(wp), INTENT(in   )                     ::   zt       ! height for t_zt and q_zt                    [m]
159      REAL(wp), INTENT(in   )                     ::   zu       ! height for U_zu                             [m]
160      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) ::   T_s      ! sea surface temperature                [Kelvin]
161      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   t_zt     ! potential air temperature              [Kelvin]
162      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) ::   q_s      ! sea surface specific humidity           [kg/kg]
163      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_zt     ! specific air humidity at zt             [kg/kg]
164      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   U_zu     ! relative wind module at zu                [m/s]
165      LOGICAL , INTENT(in   )                     ::   l_use_cs ! use the cool-skin parameterization
166      LOGICAL , INTENT(in   )                     ::   l_use_wl ! use the warm-layer parameterization
167      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cd       ! transfer coefficient for momentum         (tau)
168      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ch       ! transfer coefficient for sensible heat (Q_sens)
169      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ce       ! transfert coefficient for evaporation   (Q_lat)
170      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   t_zu     ! pot. air temp. adjusted at zu               [K]
171      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   q_zu     ! spec. humidity adjusted at zu           [kg/kg]
172      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   U_blk    ! bulk wind speed at zu                     [m/s]
173      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients
174      !
175      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   Qsw      !             [W/m^2]
176      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   rad_lw   !             [W/m^2]
177      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   slp      !             [Pa]
178      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pdT_cs
179      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pdT_wl   !             [K]
180      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pHz_wl   !             [m]
181      !
182      INTEGER :: j_itt
183      LOGICAL :: l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U
184      !
185      REAL(wp), DIMENSION(jpi,jpj) ::  &
186         &  u_star, t_star, q_star, &
187         &  dt_zu, dq_zu,    &
188         &  znu_a,           & !: Nu_air, Viscosity of air
189         &  Linv,            & !: 1/L (inverse of Monin Obukhov length...
190         &  z0, z0t, z0q
191      !
192      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
193         &                zsst,   &  ! to back up the initial bulk SST
194         &                pdTc,   &  ! SST increment "dT" for cool-skin correction           [K]
195         &                pdTw       ! SST increment "dT" for warm layer correction          [K]
196      !
197      REAL(wp), DIMENSION(jpi,jpj) ::   func_m, func_h
198      REAL(wp), DIMENSION(jpi,jpj) ::   ztmp0, ztmp1, ztmp2
199      CHARACTER(len=40), PARAMETER :: crtnm = 'turb_ecmwf@sbcblk_algo_ecmwf.F90'
200      !!----------------------------------------------------------------------------------
201
202      IF ( kt == nit000 ) CALL ECMWF_INIT(l_use_cs, l_use_wl)
203
204      l_zt_equal_zu = .FALSE.
205      IF( ABS(zu - zt) < 0.01_wp )   l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision
206
207      !! Initializations for cool skin and warm layer:
208      IF ( l_use_cs ) THEN
209         IF( .NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp)) ) THEN
210            PRINT *, ' * PROBLEM ('//TRIM(crtnm)//'): you need to provide Qsw, rad_lw & slp to use cool-skin param!'; STOP
211         END IF
212      END IF
213
214      IF ( l_use_wl ) THEN
215         IF( .NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) THEN
216            PRINT *, ' * PROBLEM ('//TRIM(crtnm)//'): you need to provide Qsw, rad_lw & slp to use warm-layer param!'; STOP
217         END IF
218      END IF
219
220      IF ( l_use_cs .OR. l_use_wl ) THEN
221         ALLOCATE ( zsst(jpi,jpj) )
222         zsst = T_s ! backing up the bulk SST
223         IF( l_use_cs ) T_s = T_s - 0.25_wp   ! First guess of correction
224         q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s !LOLO WL too!!!
225      END IF
226
227
228      ! Identical first gess as in COARE, with IFS parameter values though...
229      !
230      !! First guess of temperature and humidity at height zu:
231      t_zu = MAX( t_zt ,  180._wp )   ! who knows what's given on masked-continental regions...
232      q_zu = MAX( q_zt , 1.e-6_wp )   !               "
233
234      !! Pot. temp. difference (and we don't want it to be 0!)
235      dt_zu = t_zu - T_s ;   dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
236      dq_zu = q_zu - q_s ;   dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
237
238      znu_a = visc_air(t_zu) ! Air viscosity (m^2/s) at zt given from temperature in (K)
239
240      U_blk = SQRT(U_zu*U_zu + 0.5_wp*0.5_wp) ! initial guess for wind gustiness contribution
241
242      ztmp0   = LOG(    zu*10000._wp) ! optimization: 10000. == 1/z0 (with z0 first guess == 0.0001)
243      ztmp1   = LOG(10._wp*10000._wp) !       "                    "               "
244      u_star = 0.035_wp*U_blk*ztmp1/ztmp0       ! (u* = 0.035*Un10)
245
246      z0     = charn0*u_star*u_star/grav + 0.11_wp*znu_a/u_star
247      z0     = MIN(ABS(z0), 0.001_wp)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO
248      z0t    = 1._wp / ( 0.1_wp*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) )
249      z0t    = MIN(ABS(z0t), 0.001_wp)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO
250
251      Cd     = (vkarmn/ztmp0)**2    ! first guess of Cd
252
253      ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd
254
255      ztmp2 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN)
256
257      !! First estimate of zeta_u, depending on the stability, ie sign of BRN (ztmp2):
258      ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 )
259      func_m = ztmp0*ztmp2 ! temporary array !!
260      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
261         &  +     ztmp1   * (func_m*(1._wp + 27._wp/9._wp*ztmp2/func_m))              !  BRN > 0
262      !#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" !
263
264      !! First guess M-O stability dependent scaling params.(u*,t*,q*) to estimate z0 and z/L
265      ztmp0  = vkarmn/(LOG(zu/z0t) - psi_h_ecmwf(func_h))
266
267      u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0)  - psi_m_ecmwf(func_h))
268      t_star = dt_zu*ztmp0
269      q_star = dq_zu*ztmp0
270
271      ! What needs to be done if zt /= zu:
272      IF( .NOT. l_zt_equal_zu ) THEN
273         !! First update of values at zu (or zt for wind)
274         ztmp0 = psi_h_ecmwf(func_h) - psi_h_ecmwf(zt*func_h/zu)    ! zt*func_h/zu == zeta_t
275         ztmp1 = LOG(zt/zu) + ztmp0
276         t_zu = t_zt - t_star/vkarmn*ztmp1
277         q_zu = q_zt - q_star/vkarmn*ztmp1
278         q_zu = (0.5_wp + SIGN(0.5_wp,q_zu))*q_zu !Makes it impossible to have negative humidity :
279         !
280         dt_zu = t_zu - T_s  ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
281         dq_zu = q_zu - q_s  ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
282      END IF
283
284
285      !! => that was same first guess as in COARE...
286
287
288      !! First guess of inverse of Monin-Obukov length (1/L) :
289      Linv = One_on_L( t_zu, q_zu, u_star, t_star, q_star )
290
291      !! Functions such as  u* = U_blk*vkarmn/func_m
292      ztmp0 = zu*Linv
293      func_m = LOG(zu) - LOG(z0)  - psi_m_ecmwf(ztmp0) + psi_m_ecmwf( z0*Linv)
294      func_h = LOG(zu) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv)
295
296      !! ITERATION BLOCK
297      DO j_itt = 1, nb_itt
298
299         !! Bulk Richardson Number at z=zu (Eq. 3.25)
300         ztmp0 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN)
301
302         !! New estimate of the inverse of the Monin-Obukhon length (Linv == zeta/zu) :
303         Linv = ztmp0*func_m*func_m/func_h / zu     ! From Eq. 3.23, Chap.3.2.3, IFS doc - Cy40r1
304         !! Note: it is slightly different that the L we would get with the usual
305         Linv = SIGN( MIN(ABS(Linv),200._wp), Linv ) ! (prevent FPE from stupid values from masked region later on...) !#LOLO
306
307         !! Update func_m with new Linv:
308         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!
309
310         !! Need to update roughness lengthes:
311         u_star = U_blk*vkarmn/func_m
312         ztmp2  = u_star*u_star
313         ztmp1  = znu_a/u_star
314         z0     = MIN( ABS( alpha_M*ztmp1 + charn0*ztmp2/grav ) , 0.001_wp)
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         END IF
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         END IF
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         END IF
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         END IF
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      !
423      DO jj = 1, jpj
424         DO ji = 1, jpi
425            !
426            zzeta = MIN( pzeta(ji,jj) , 5._wp ) !! Very stable conditions (L positif and big!):
427            !
428            ! Unstable (Paulson 1970):
429            !   eq.3.20, Chap.3, p.33, IFS doc - Cy31r1
430            zx = SQRT(ABS(1._wp - 16._wp*zzeta))
431            ztmp = 1._wp + SQRT(zx)
432            ztmp = ztmp*ztmp
433            psi_unst = LOG( 0.125_wp*ztmp*(1._wp + zx) )   &
434               &       -2._wp*ATAN( SQRT(zx) ) + 0.5_wp*rpi
435            !
436            ! Unstable:
437            ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1
438            psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) &
439               &       - zzeta - 2._wp/3._wp*5._wp/0.35_wp
440            !
441            ! Combining:
442            stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1
443            !
444            psi_m_ecmwf(ji,jj) = (1._wp - stab) * psi_unst & ! (zzeta < 0) Unstable
445               &                +      stab  * psi_stab      ! (zzeta > 0) Stable
446            !
447         END DO
448      END DO
449      !
450   END FUNCTION psi_m_ecmwf
451
452
453   FUNCTION psi_h_ecmwf( pzeta )
454      !!----------------------------------------------------------------------------------
455      !! Universal profile stability function for temperature and humidity
456      !!     ECMWF / as in IFS cy31r1 documentation, available online
457      !!     at ecmwf.int
458      !!
459      !! pzeta : stability paramenter, z/L where z is altitude measurement
460      !!         and L is M-O length
461      !!
462      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
463      !!----------------------------------------------------------------------------------
464      REAL(wp), DIMENSION(jpi,jpj) :: psi_h_ecmwf
465      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
466      !
467      INTEGER  ::   ji, jj     ! dummy loop indices
468      REAL(wp) ::  zzeta, zx, psi_unst, psi_stab, stab
469      !!----------------------------------------------------------------------------------
470      !
471      DO jj = 1, jpj
472         DO ji = 1, jpi
473            !
474            zzeta = MIN(pzeta(ji,jj) , 5._wp)   ! Very stable conditions (L positif and big!):
475            !
476            zx  = ABS(1._wp - 16._wp*zzeta)**.25        ! this is actually (1/phi_m)**2  !!!
477            !                                     ! eq.3.19, Chap.3, p.33, IFS doc - Cy31r1
478            ! Unstable (Paulson 1970) :
479            psi_unst = 2._wp*LOG(0.5_wp*(1._wp + zx*zx))   ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1
480            !
481            ! Stable:
482            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
483               &       - ABS(1._wp + 2._wp/3._wp*zzeta)**1.5_wp - 2._wp/3._wp*5._wp/0.35_wp + 1._wp
484            ! LB: added ABS() to avoid NaN values when unstable, which contaminates the unstable solution...
485            !
486            stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1
487            !
488            !
489            psi_h_ecmwf(ji,jj) = (1._wp - stab) * psi_unst &   ! (zzeta < 0) Unstable
490               &                +    stab    * psi_stab        ! (zzeta > 0) Stable
491            !
492         END DO
493      END DO
494      !
495   END FUNCTION psi_h_ecmwf
496
497
498   !!======================================================================
499END MODULE sbcblk_algo_ecmwf
Note: See TracBrowser for help on using the repository browser.