source: NEMO/trunk/src/OCE/SBC/sbcblk_algo_ecmwf.F90 @ 12377

Last change on this file since 12377 was 12377, checked in by acc, 6 months ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge —ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The —ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

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