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

Last change on this file since 11615 was 11615, checked in by laurent, 4 years ago

LB: new cool-skin and warm-layer parameterizations for ECMWF and COARE3.6, COARE3.0 uses same CSWL param as for COARE3.6.

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