Ignore:
Timestamp:
2019-12-10T15:44:23+01:00 (10 months ago)
Author:
cetlod
Message:

commit

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/SBC/sbcblk_algo_ecmwf.F90

    r10069 r12154  
    11MODULE sbcblk_algo_ecmwf 
    22   !!====================================================================== 
    3    !!                       ***  MODULE  sbcblk_algo_ecmwf  *** 
    4    !! Computes turbulent components of surface fluxes 
    5    !!         according to the method in IFS of the ECMWF model 
    6    !! 
     3   !!                   ***  MODULE  sbcblk_algo_ecmwf  *** 
     4   !! Computes: 
    75   !!   * bulk transfer coefficients C_D, C_E and C_H 
    86   !!   * air temp. and spec. hum. adjusted from zt (2m) to zu (10m) if needed 
     
    108   !!   => all these are used in bulk formulas in sbcblk.F90 
    119   !! 
    12    !!    Using the bulk formulation/param. of IFS of ECMWF (cycle 31r2) 
     10   !!    Using the bulk formulation/param. of IFS of ECMWF (cycle 40r1) 
    1311   !!         based on IFS doc (avaible online on the ECMWF's website) 
    1412   !! 
     13   !!       Routine turb_ecmwf maintained and developed in AeroBulk 
     14   !!                     (https://github.com/brodeau/aerobulk) 
    1515   !! 
    16    !!       Routine turb_ecmwf maintained and developed in AeroBulk 
    17    !!                     (http://aerobulk.sourceforge.net/) 
    18    !! 
    19    !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
     16   !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk) 
    2017   !!---------------------------------------------------------------------- 
    2118   !! History :  4.0  !  2016-02  (L.Brodeau)   Original code 
     
    4138 
    4239   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 
    4342 
    4443   IMPLICIT NONE 
    4544   PRIVATE 
    4645 
    47    PUBLIC ::   TURB_ECMWF   ! called by sbcblk.F90 
    48  
    49    !                   !! ECMWF own values for given constants, taken form IFS documentation... 
     46   PUBLIC :: SBCBLK_ALGO_ECMWF_INIT, TURB_ECMWF 
     47 
     48   !! ECMWF own values for given constants, taken form IFS documentation... 
    5049   REAL(wp), PARAMETER ::   charn0 = 0.018    ! Charnock constant (pretty high value here !!! 
    5150   !                                          !    =>  Usually 0.011 for moderate winds) 
    5251   REAL(wp), PARAMETER ::   zi0     = 1000.   ! scale height of the atmospheric boundary layer...1 
    5352   REAL(wp), PARAMETER ::   Beta0    = 1.     ! gustiness parameter ( = 1.25 in COAREv3) 
    54    REAL(wp), PARAMETER ::   rctv0    = 0.608  ! constant to obtain virtual temperature... 
    55    REAL(wp), PARAMETER ::   Cp_dry = 1005.0   ! Specic heat of dry air, constant pressure      [J/K/kg] 
    56    REAL(wp), PARAMETER ::   Cp_vap = 1860.0   ! Specic heat of water vapor, constant pressure  [J/K/kg] 
    5753   REAL(wp), PARAMETER ::   alpha_M = 0.11    ! For roughness length (smooth surface term) 
    5854   REAL(wp), PARAMETER ::   alpha_H = 0.40    ! (Chapter 3, p.34, IFS doc Cy31r1) 
    5955   REAL(wp), PARAMETER ::   alpha_Q = 0.62    ! 
     56 
     57   INTEGER , PARAMETER ::   nb_itt = 10             ! number of itterations 
     58 
    6059   !!---------------------------------------------------------------------- 
    6160CONTAINS 
    6261 
    63    SUBROUTINE TURB_ECMWF( zt, zu, sst, t_zt, ssq , q_zt , U_zu,   & 
    64       &                   Cd, Ch, Ce , t_zu, q_zu, U_blk,         & 
    65       &                   Cdn, Chn, Cen                           ) 
    66       !!---------------------------------------------------------------------------------- 
    67       !!                      ***  ROUTINE  turb_ecmwf  *** 
    68       !! 
    69       !!            2015: L. Brodeau (brodeau@gmail.com) 
    70       !! 
    71       !! ** Purpose :   Computes turbulent transfert coefficients of surface 
    72       !!                fluxes according to IFS doc. (cycle 31) 
    73       !!                If relevant (zt /= zu), adjust temperature and humidity from height zt to zu 
    74       !! 
    75       !! ** Method : Monin Obukhov Similarity Theory 
     62 
     63   SUBROUTINE sbcblk_algo_ecmwf_init(l_use_cs, l_use_wl) 
     64      !!--------------------------------------------------------------------- 
     65      !!                  ***  FUNCTION sbcblk_algo_ecmwf_init  *** 
    7666      !! 
    7767      !! INPUT : 
    7868      !! ------- 
     69      !!    * l_use_cs : use the cool-skin parameterization 
     70      !!    * l_use_wl : use the warm-layer parameterization 
     71      !!--------------------------------------------------------------------- 
     72      LOGICAL , INTENT(in) ::   l_use_cs ! use the cool-skin parameterization 
     73      LOGICAL , INTENT(in) ::   l_use_wl ! use the warm-layer parameterization 
     74      INTEGER :: ierr 
     75      !!--------------------------------------------------------------------- 
     76      IF( l_use_wl ) THEN 
     77         ierr = 0 
     78         ALLOCATE ( dT_wl(jpi,jpj), Hz_wl(jpi,jpj), STAT=ierr ) 
     79         IF( ierr > 0 ) CALL ctl_stop( ' SBCBLK_ALGO_ECMWF_INIT => allocation of dT_wl & Hz_wl failed!' ) 
     80         dT_wl(:,:)  = 0._wp 
     81         Hz_wl(:,:)  = rd0 ! (rd0, constant, = 3m is default for Zeng & Beljaars) 
     82      ENDIF 
     83      IF( l_use_cs ) THEN 
     84         ierr = 0 
     85         ALLOCATE ( dT_cs(jpi,jpj), STAT=ierr ) 
     86         IF( ierr > 0 ) CALL ctl_stop( ' SBCBLK_ALGO_ECMWF_INIT => allocation of dT_cs failed!' ) 
     87         dT_cs(:,:) = -0.25_wp  ! First guess of skin correction 
     88      ENDIF 
     89   END SUBROUTINE sbcblk_algo_ecmwf_init 
     90 
     91 
     92 
     93   SUBROUTINE turb_ecmwf( kt, zt, zu, T_s, t_zt, q_s, q_zt, U_zu, l_use_cs, l_use_wl, & 
     94      &                      Cd, Ch, Ce, t_zu, q_zu, U_blk,                           & 
     95      &                      Cdn, Chn, Cen,                                           & 
     96      &                      Qsw, rad_lw, slp, pdT_cs,                                & ! optionals for cool-skin (and warm-layer) 
     97      &                      pdT_wl, pHz_wl )                                           ! optionals for warm-layer only 
     98      !!---------------------------------------------------------------------- 
     99      !!                      ***  ROUTINE  turb_ecmwf  *** 
     100      !! 
     101      !! ** Purpose :   Computes turbulent transfert coefficients of surface 
     102      !!                fluxes according to IFS doc. (cycle 45r1) 
     103      !!                If relevant (zt /= zu), adjust temperature and humidity from height zt to zu 
     104      !!                Returns the effective bulk wind speed at zu to be used in the bulk formulas 
     105      !! 
     106      !!                Applies the cool-skin warm-layer correction of the SST to T_s 
     107      !!                if the net shortwave flux at the surface (Qsw), the downwelling longwave 
     108      !!                radiative fluxes at the surface (rad_lw), and the sea-leve pressure (slp) 
     109      !!                are provided as (optional) arguments! 
     110      !! 
     111      !! INPUT : 
     112      !! ------- 
     113      !!    *  kt   : current time step (starts at 1) 
    79114      !!    *  zt   : height for temperature and spec. hum. of air            [m] 
    80       !!    *  zu   : height for wind speed (generally 10m)                   [m] 
    81       !!    *  U_zu : scalar wind speed at 10m                                [m/s] 
    82       !!    *  sst  : SST                                                     [K] 
     115      !!    *  zu   : height for wind speed (usually 10m)                     [m] 
    83116      !!    *  t_zt : potential air temperature at zt                         [K] 
    84       !!    *  ssq  : specific humidity at saturation at SST                  [kg/kg] 
    85117      !!    *  q_zt : specific humidity of air at zt                          [kg/kg] 
    86       !! 
     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] 
    87143      !! 
    88144      !! OUTPUT : 
     
    93149      !!    *  t_zu   : pot. air temperature adjusted at wind height zu       [K] 
    94150      !!    *  q_zu   : specific humidity of air        //                    [kg/kg] 
    95       !!    *  U_blk  : bulk wind at 10m                                      [m/s] 
    96       !! 
    97       !! 
    98       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    99       !!---------------------------------------------------------------------------------- 
     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 
    100157      REAL(wp), INTENT(in   )                     ::   zt       ! height for t_zt and q_zt                    [m] 
    101158      REAL(wp), INTENT(in   )                     ::   zu       ! height for U_zu                             [m] 
    102       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   sst      ! sea surface temperature                [Kelvin] 
     159      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) ::   T_s      ! sea surface temperature                [Kelvin] 
    103160      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   t_zt     ! potential air temperature              [Kelvin] 
    104       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   ssq      ! sea surface specific humidity           [kg/kg] 
    105       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_zt     ! specific air humidity                   [kg/kg] 
     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] 
    106163      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 
    107166      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cd       ! transfer coefficient for momentum         (tau) 
    108167      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ch       ! transfer coefficient for sensible heat (Q_sens) 
     
    110169      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   t_zu     ! pot. air temp. adjusted at zu               [K] 
    111170      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   q_zu     ! spec. humidity adjusted at zu           [kg/kg] 
    112       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   U_blk    ! bulk wind at 10m                          [m/s] 
     171      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   U_blk    ! bulk wind speed at zu                     [m/s] 
    113172      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients 
    114173      ! 
     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      ! 
    115181      INTEGER :: j_itt 
    116       LOGICAL ::   l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U 
    117       INTEGER , PARAMETER ::   nb_itt = 4       ! number of itterations 
    118       ! 
    119       REAL(wp), DIMENSION(jpi,jpj) ::   u_star, t_star, q_star,   & 
    120          &  dt_zu, dq_zu,    & 
    121          &  znu_a,           & !: Nu_air, Viscosity of air 
    122          &  Linv,            & !: 1/L (inverse of Monin Obukhov length... 
    123          &  z0, z0t, z0q 
    124       REAL(wp), DIMENSION(jpi,jpj) ::   func_m, func_h 
    125       REAL(wp), DIMENSION(jpi,jpj) ::   ztmp0, ztmp1, ztmp2 
    126       !!---------------------------------------------------------------------------------- 
    127       ! 
    128       ! Identical first gess as in COARE, with IFS parameter values though 
    129       ! 
     182      LOGICAL :: l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U 
     183      ! 
     184      REAL(wp), DIMENSION(jpi,jpj) ::  u_star, t_star, q_star 
     185      REAL(wp), DIMENSION(jpi,jpj) :: dt_zu, dq_zu      
     186      REAL(wp), DIMENSION(jpi,jpj) :: znu_a !: Nu_air, Viscosity of air 
     187      REAL(wp), DIMENSION(jpi,jpj) :: Linv  !: 1/L (inverse of Monin Obukhov length... 
     188      REAL(wp), DIMENSION(jpi,jpj) :: z0, z0t, z0q 
     189      ! 
     190      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zsst  ! to back up the initial bulk SST 
     191      ! 
     192      REAL(wp), DIMENSION(jpi,jpj) :: func_m, func_h 
     193      REAL(wp), DIMENSION(jpi,jpj) :: ztmp0, ztmp1, ztmp2 
     194      CHARACTER(len=40), PARAMETER :: crtnm = 'turb_ecmwf@sbcblk_algo_ecmwf.F90' 
     195      !!---------------------------------------------------------------------------------- 
     196 
     197      IF( kt == nit000 ) CALL SBCBLK_ALGO_ECMWF_INIT(l_use_cs, l_use_wl) 
     198 
    130199      l_zt_equal_zu = .FALSE. 
    131       IF( ABS(zu - zt) < 0.01 )   l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision 
    132  
    133  
     200      IF( ABS(zu - zt) < 0.01_wp )   l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision 
     201 
     202      !! Initializations for cool skin and warm layer: 
     203      IF( l_use_cs .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) & 
     204         &   CALL ctl_stop( '['//TRIM(crtnm)//'] => ' , 'you need to provide Qsw, rad_lw & slp to use cool-skin param!' ) 
     205 
     206      IF( l_use_wl .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) & 
     207         &   CALL ctl_stop( '['//TRIM(crtnm)//'] => ' , 'you need to provide Qsw, rad_lw & slp to use warm-layer param!' ) 
     208 
     209      IF( l_use_cs .OR. l_use_wl ) THEN 
     210         ALLOCATE ( zsst(jpi,jpj) ) 
     211         zsst = T_s ! backing up the bulk SST 
     212         IF( l_use_cs ) T_s = T_s - 0.25_wp   ! First guess of correction 
     213         q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s 
     214      ENDIF 
     215 
     216 
     217      ! Identical first gess as in COARE, with IFS parameter values though... 
     218      ! 
    134219      !! First guess of temperature and humidity at height zu: 
    135       t_zu = MAX( t_zt , 0.0 )   ! who knows what's given on masked-continental regions... 
    136       q_zu = MAX( q_zt , 1.e-6)   !               " 
     220      t_zu = MAX( t_zt ,  180._wp )   ! who knows what's given on masked-continental regions... 
     221      q_zu = MAX( q_zt , 1.e-6_wp )   !               " 
    137222 
    138223      !! Pot. temp. difference (and we don't want it to be 0!) 
    139       dt_zu = t_zu - sst   ;   dt_zu = SIGN( MAX(ABS(dt_zu),1.e-6), dt_zu ) 
    140       dq_zu = q_zu - ssq   ;   dq_zu = SIGN( MAX(ABS(dq_zu),1.e-9), dq_zu ) 
    141  
    142       znu_a = visc_air(t_zt) ! Air viscosity (m^2/s) at zt given from temperature in (K) 
    143  
    144       ztmp2 = 0.5 * 0.5 ! initial guess for wind gustiness contribution 
    145       U_blk = SQRT(U_zu*U_zu + ztmp2) 
    146  
    147       ! z0     = 0.0001 
    148       ztmp2   = 10000.     ! optimization: ztmp2 == 1/z0 
    149       ztmp0   = LOG(zu*ztmp2) 
    150       ztmp1   = LOG(10.*ztmp2) 
    151       u_star = 0.035*U_blk*ztmp1/ztmp0       ! (u* = 0.035*Un10) 
    152  
    153       z0     = charn0*u_star*u_star/grav + 0.11*znu_a/u_star 
    154       z0t    = 0.1*EXP(vkarmn/(0.00115/(vkarmn/ztmp1)))   !  WARNING: 1/z0t ! 
     224      dt_zu = t_zu - T_s ;   dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 
     225      dq_zu = q_zu - q_s ;   dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 
     226 
     227      znu_a = visc_air(t_zu) ! Air viscosity (m^2/s) at zt given from temperature in (K) 
     228 
     229      U_blk = SQRT(U_zu*U_zu + 0.5_wp*0.5_wp) ! initial guess for wind gustiness contribution 
     230 
     231      ztmp0   = LOG(    zu*10000._wp) ! optimization: 10000. == 1/z0 (with z0 first guess == 0.0001) 
     232      ztmp1   = LOG(10._wp*10000._wp) !       "                    "               " 
     233      u_star = 0.035_wp*U_blk*ztmp1/ztmp0       ! (u* = 0.035*Un10) 
     234 
     235      z0     = charn0*u_star*u_star/grav + 0.11_wp*znu_a/u_star 
     236      z0     = MIN( MAX(ABS(z0), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on) 
     237 
     238      z0t    = 1._wp / ( 0.1_wp*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) ) 
     239      z0t    = MIN( MAX(ABS(z0t), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on) 
    155240 
    156241      Cd     = (vkarmn/ztmp0)**2    ! first guess of Cd 
    157242 
    158       ztmp0 = vkarmn*vkarmn/LOG(zt*z0t)/Cd 
    159  
    160       ztmp2 = Ri_bulk( zu, t_zu, dt_zu, q_zu, dq_zu, U_blk )   ! Ribu = Bulk Richardson number 
    161  
    162       !! First estimate of zeta_u, depending on the stability, ie sign of Ribu (ztmp2): 
    163       ztmp1 = 0.5 + SIGN( 0.5 , ztmp2 ) 
     243      ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd 
     244 
     245      ztmp2 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN) 
     246 
     247      !! First estimate of zeta_u, depending on the stability, ie sign of BRN (ztmp2): 
     248      ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 ) 
    164249      func_m = ztmp0*ztmp2 ! temporary array !! 
    165       !!             Ribu < 0                                 Ribu > 0   Beta = 1.25 
    166       func_h = (1.-ztmp1)*(func_m/(1.+ztmp2/(-zu/(zi0*0.004*Beta0**3)))) &  ! temporary array !!! func_h == zeta_u 
    167          &  +     ztmp1*(func_m*(1. + 27./9.*ztmp2/ztmp0)) 
     250      func_h = (1._wp-ztmp1) * (func_m/(1._wp+ztmp2/(-zu/(zi0*0.004_wp*Beta0**3)))) & !  BRN < 0 ! temporary array !!! func_h == zeta_u 
     251         &  +     ztmp1   * (func_m*(1._wp + 27._wp/9._wp*ztmp2/func_m))              !  BRN > 0 
     252      !#LB: should make sure that the "func_m" of "27./9.*ztmp2/func_m" is "ztmp0*ztmp2" and not "ztmp0==vkarmn*vkarmn/LOG(zt/z0t)/Cd" ! 
    168253 
    169254      !! First guess M-O stability dependent scaling params.(u*,t*,q*) to estimate z0 and z/L 
    170       ztmp0   =        vkarmn/(LOG(zu*z0t) - psi_h_ecmwf(func_h)) 
    171  
    172       u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0)  - psi_m_ecmwf(func_h)) 
     255      ztmp0  = vkarmn/(LOG(zu/z0t) - psi_h_ecmwf(func_h)) 
     256 
     257      u_star = MAX ( U_blk*vkarmn/(LOG(zu) - LOG(z0)  - psi_m_ecmwf(func_h)) , 1.E-9 )  !  (MAX => prevents FPE from stupid values from masked region later on) 
    173258      t_star = dt_zu*ztmp0 
    174259      q_star = dq_zu*ztmp0 
    175260 
    176       ! What's need to be done if zt /= zu: 
     261      ! What needs to be done if zt /= zu: 
    177262      IF( .NOT. l_zt_equal_zu ) THEN 
    178          ! 
    179263         !! First update of values at zu (or zt for wind) 
    180264         ztmp0 = psi_h_ecmwf(func_h) - psi_h_ecmwf(zt*func_h/zu)    ! zt*func_h/zu == zeta_t 
    181          ztmp1 = log(zt/zu) + ztmp0 
     265         ztmp1 = LOG(zt/zu) + ztmp0 
    182266         t_zu = t_zt - t_star/vkarmn*ztmp1 
    183267         q_zu = q_zt - q_star/vkarmn*ztmp1 
    184          q_zu = (0.5 + sign(0.5,q_zu))*q_zu !Makes it impossible to have negative humidity : 
    185  
    186          dt_zu = t_zu - sst  ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6), dt_zu ) 
    187          dq_zu = q_zu - ssq  ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9), dq_zu ) 
     268         q_zu = (0.5_wp + SIGN(0.5_wp,q_zu))*q_zu !Makes it impossible to have negative humidity : 
    188269         ! 
     270         dt_zu = t_zu - T_s  ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 
     271         dq_zu = q_zu - q_s  ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 
    189272      ENDIF 
    190273 
     
    194277 
    195278      !! First guess of inverse of Monin-Obukov length (1/L) : 
    196       ztmp0 = (1. + rctv0*q_zu)  ! the factor to apply to temp. to get virt. temp... 
    197       Linv  =  grav*vkarmn*(t_star*ztmp0 + rctv0*t_zu*q_star) / ( u_star*u_star * t_zu*ztmp0 ) 
     279      Linv = One_on_L( t_zu, q_zu, u_star, t_star, q_star ) 
    198280 
    199281      !! Functions such as  u* = U_blk*vkarmn/func_m 
    200       ztmp1 = zu + z0 
    201       ztmp0 = ztmp1*Linv 
    202       func_m = LOG(ztmp1) -LOG(z0) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf(z0*Linv) 
    203       func_h = LOG(ztmp1*z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(1./z0t*Linv) 
    204  
     282      ztmp0 = zu*Linv 
     283      func_m = LOG(zu) - LOG(z0)  - psi_m_ecmwf(ztmp0) + psi_m_ecmwf( z0*Linv) 
     284      func_h = LOG(zu) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv) 
    205285 
    206286      !! ITERATION BLOCK 
    207       !! *************** 
    208  
    209287      DO j_itt = 1, nb_itt 
    210288 
    211289         !! Bulk Richardson Number at z=zu (Eq. 3.25) 
    212          ztmp0 = Ri_bulk(zu, t_zu, dt_zu, q_zu, dq_zu, U_blk) 
     290         ztmp0 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN) 
    213291 
    214292         !! New estimate of the inverse of the Monin-Obukhon length (Linv == zeta/zu) : 
    215          Linv = ztmp0*func_m*func_m/func_h / zu     ! From Eq. 3.23, Chap.3, p.33, IFS doc - Cy31r1 
     293         Linv = ztmp0*func_m*func_m/func_h / zu     ! From Eq. 3.23, Chap.3.2.3, IFS doc - Cy40r1 
     294         !! Note: it is slightly different that the L we would get with the usual 
     295         Linv = SIGN( MIN(ABS(Linv),200._wp), Linv ) ! (prevent FPE from stupid values from masked region later on...) 
    216296 
    217297         !! Update func_m with new Linv: 
    218          ztmp1 = zu + z0 
    219          func_m = LOG(ztmp1) -LOG(z0) - psi_m_ecmwf(ztmp1*Linv) + psi_m_ecmwf(z0*Linv) 
     298         func_m = LOG(zu) -LOG(z0) - psi_m_ecmwf(zu*Linv) + psi_m_ecmwf(z0*Linv) ! LB: should be "zu+z0" rather than "zu" alone, but z0 is tiny wrt zu! 
    220299 
    221300         !! Need to update roughness lengthes: 
     
    223302         ztmp2  = u_star*u_star 
    224303         ztmp1  = znu_a/u_star 
    225          z0    = alpha_M*ztmp1 + charn0*ztmp2/grav 
    226          z0t    = alpha_H*ztmp1                              ! eq.3.26, Chap.3, p.34, IFS doc - Cy31r1 
    227          z0q    = alpha_Q*ztmp1 
    228  
    229          !! Update wind at 10m taking into acount convection-related wind gustiness: 
    230          ! Only true when unstable (L<0) => when ztmp0 < 0 => - !!! 
    231          ztmp2 = ztmp2 * (MAX(-zi0*Linv/vkarmn,0.))**(2./3.) ! => w*^2  (combining Eq. 3.8 and 3.18, hap.3, IFS doc - Cy31r1) 
    232          !! => equivalent using Beta=1 (gustiness parameter, 1.25 for COARE, also zi0=600 in COARE..) 
    233          U_blk = MAX(sqrt(U_zu*U_zu + ztmp2), 0.2)              ! eq.3.17, Chap.3, p.32, IFS doc - Cy31r1 
     304         z0     = MIN( ABS( alpha_M*ztmp1 + charn0*ztmp2/grav ) , 0.001_wp) 
     305         z0t    = MIN( ABS( alpha_H*ztmp1                     ) , 0.001_wp)   ! eq.3.26, Chap.3, p.34, IFS doc - Cy31r1 
     306         z0q    = MIN( ABS( alpha_Q*ztmp1                     ) , 0.001_wp) 
     307 
     308         !! Update wind at zu with convection-related wind gustiness in unstable conditions (Chap. 3.2, IFS doc - Cy40r1, Eq.3.17 and Eq.3.18 + Eq.3.8) 
     309         ztmp2 = Beta0*Beta0*ztmp2*(MAX(-zi0*Linv/vkarmn,0._wp))**(2._wp/3._wp) ! square of wind gustiness contribution  (combining Eq. 3.8 and 3.18, hap.3, IFS doc - Cy31r1) 
     310         !!   ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before zi0 
     311         U_blk = MAX(SQRT(U_zu*U_zu + ztmp2), 0.2_wp)        ! include gustiness in bulk wind speed 
    234312         ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0. 
    235313 
     
    238316         !! as well the air-sea differences: 
    239317         IF( .NOT. l_zt_equal_zu ) THEN 
    240  
    241318            !! Arrays func_m and func_h are free for a while so using them as temporary arrays... 
    242             func_h = psi_h_ecmwf((zu+z0)*Linv) ! temporary array !!! 
    243             func_m = psi_h_ecmwf((zt+z0)*Linv) ! temporary array !!! 
     319            func_h = psi_h_ecmwf(zu*Linv) ! temporary array !!! 
     320            func_m = psi_h_ecmwf(zt*Linv) ! temporary array !!! 
    244321 
    245322            ztmp2  = psi_h_ecmwf(z0t*Linv) 
    246323            ztmp0  = func_h - ztmp2 
    247             ztmp1  = vkarmn/(LOG(zu+z0) - LOG(z0t) - ztmp0) 
     324            ztmp1  = vkarmn/(LOG(zu) - LOG(z0t) - ztmp0) 
    248325            t_star = dt_zu*ztmp1 
    249326            ztmp2  = ztmp0 - func_m + ztmp2 
     
    253330            ztmp2  = psi_h_ecmwf(z0q*Linv) 
    254331            ztmp0  = func_h - ztmp2 
    255             ztmp1  = vkarmn/(LOG(zu+z0) - LOG(z0q) - ztmp0) 
     332            ztmp1  = vkarmn/(LOG(zu) - LOG(z0q) - ztmp0) 
    256333            q_star = dq_zu*ztmp1 
    257334            ztmp2  = ztmp0 - func_m + ztmp2 
    258             ztmp1  = log(zt/zu) + ztmp2 
     335            ztmp1  = LOG(zt/zu) + ztmp2 
    259336            q_zu   = q_zt - q_star/vkarmn*ztmp1 
    260  
    261             dt_zu = t_zu - sst ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6), dt_zu ) 
    262             dq_zu = q_zu - ssq ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9), dq_zu ) 
    263  
    264          END IF 
     337         ENDIF 
    265338 
    266339         !! Updating because of updated z0 and z0t and new Linv... 
    267          ztmp1 = zu + z0 
    268          ztmp0 = ztmp1*Linv 
    269          func_m = log(ztmp1) - LOG(z0 ) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf(z0 *Linv) 
    270          func_h = log(ztmp1) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv) 
    271  
    272       END DO 
     340         ztmp0 = zu*Linv 
     341         func_m = log(zu) - LOG(z0 ) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf(z0 *Linv) 
     342         func_h = log(zu) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv) 
     343 
     344 
     345         IF( l_use_cs ) THEN 
     346            !! Cool-skin contribution 
     347 
     348            CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, U_blk, slp, rad_lw, & 
     349               &                   ztmp1, ztmp0,  Qlat=ztmp2)  ! Qnsol -> ztmp1 / Tau -> ztmp0 
     350 
     351            CALL CS_ECMWF( Qsw, ztmp1, u_star, zsst )  ! Qnsol -> ztmp1 
     352 
     353            T_s(:,:) = zsst(:,:) + dT_cs(:,:)*tmask(:,:,1) 
     354            IF( l_use_wl ) T_s(:,:) = T_s(:,:) + dT_wl(:,:)*tmask(:,:,1) 
     355            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 
     356 
     357         ENDIF 
     358 
     359         IF( l_use_wl ) THEN 
     360            !! Warm-layer contribution 
     361            CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, U_blk, slp, rad_lw, & 
     362               &                   ztmp1, ztmp2)  ! Qnsol -> ztmp1 / Tau -> ztmp2 
     363            CALL WL_ECMWF( Qsw, ztmp1, u_star, zsst ) 
     364            !! Updating T_s and q_s !!! 
     365            T_s(:,:) = zsst(:,:) + dT_wl(:,:)*tmask(:,:,1) ! 
     366            IF( l_use_cs ) T_s(:,:) = T_s(:,:) + dT_cs(:,:)*tmask(:,:,1) 
     367            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 
     368         ENDIF 
     369 
     370         IF( l_use_cs .OR. l_use_wl .OR. (.NOT. l_zt_equal_zu) ) THEN 
     371            dt_zu = t_zu - T_s ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 
     372            dq_zu = q_zu - q_s ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 
     373         ENDIF 
     374 
     375      END DO !DO j_itt = 1, nb_itt 
    273376 
    274377      Cd = vkarmn*vkarmn/(func_m*func_m) 
    275378      Ch = vkarmn*vkarmn/(func_m*func_h) 
    276       ztmp1 = log((zu + z0)/z0q) - psi_h_ecmwf((zu + z0)*Linv) + psi_h_ecmwf(z0q*Linv)   ! func_q 
    277       Ce = vkarmn*vkarmn/(func_m*ztmp1) 
    278  
    279       ztmp1 = zu + z0 
    280       Cdn = vkarmn*vkarmn / (log(ztmp1/z0 )*log(ztmp1/z0 )) 
    281       Chn = vkarmn*vkarmn / (log(ztmp1/z0t)*log(ztmp1/z0t)) 
    282       Cen = vkarmn*vkarmn / (log(ztmp1/z0q)*log(ztmp1/z0q)) 
    283  
    284    END SUBROUTINE TURB_ECMWF 
     379      ztmp2 = log(zu/z0q) - psi_h_ecmwf(zu*Linv) + psi_h_ecmwf(z0q*Linv)   ! func_q 
     380      Ce = vkarmn*vkarmn/(func_m*ztmp2) 
     381 
     382      Cdn = vkarmn*vkarmn / (log(zu/z0 )*log(zu/z0 )) 
     383      Chn = vkarmn*vkarmn / (log(zu/z0t)*log(zu/z0t)) 
     384      Cen = vkarmn*vkarmn / (log(zu/z0q)*log(zu/z0q)) 
     385 
     386      IF( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = dT_cs 
     387      IF( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = dT_wl 
     388      IF( l_use_wl .AND. PRESENT(pHz_wl) ) pHz_wl = Hz_wl 
     389 
     390      IF( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst ) 
     391 
     392   END SUBROUTINE turb_ecmwf 
    285393 
    286394 
     
    294402      !!         and L is M-O length 
    295403      !! 
    296       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
     404      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    297405      !!---------------------------------------------------------------------------------- 
    298406      REAL(wp), DIMENSION(jpi,jpj) :: psi_m_ecmwf 
     
    302410      REAL(wp) :: zzeta, zx, ztmp, psi_unst, psi_stab, stab 
    303411      !!---------------------------------------------------------------------------------- 
    304       ! 
    305412      DO jj = 1, jpj 
    306413         DO ji = 1, jpi 
    307414            ! 
    308             zzeta = MIN( pzeta(ji,jj) , 5. ) !! Very stable conditions (L positif and big!): 
     415            zzeta = MIN( pzeta(ji,jj) , 5._wp ) !! Very stable conditions (L positif and big!): 
    309416            ! 
    310417            ! Unstable (Paulson 1970): 
    311418            !   eq.3.20, Chap.3, p.33, IFS doc - Cy31r1 
    312             zx = SQRT(ABS(1. - 16.*zzeta)) 
    313             ztmp = 1. + SQRT(zx) 
     419            zx = SQRT(ABS(1._wp - 16._wp*zzeta)) 
     420            ztmp = 1._wp + SQRT(zx) 
    314421            ztmp = ztmp*ztmp 
    315             psi_unst = LOG( 0.125*ztmp*(1. + zx) )   & 
    316                &       -2.*ATAN( SQRT(zx) ) + 0.5*rpi 
     422            psi_unst = LOG( 0.125_wp*ztmp*(1._wp + zx) )   & 
     423               &       -2._wp*ATAN( SQRT(zx) ) + 0.5_wp*rpi 
    317424            ! 
    318425            ! Unstable: 
    319426            ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1 
    320             psi_stab = -2./3.*(zzeta - 5./0.35)*EXP(-0.35*zzeta) & 
    321                &       - zzeta - 2./3.*5./0.35 
     427            psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) & 
     428               &       - zzeta - 2._wp/3._wp*5._wp/0.35_wp 
    322429            ! 
    323430            ! Combining: 
    324             stab = 0.5 + SIGN(0.5, zzeta) ! zzeta > 0 => stab = 1 
    325             ! 
    326             psi_m_ecmwf(ji,jj) = (1. - stab) * psi_unst & ! (zzeta < 0) Unstable 
    327                &                +      stab  * psi_stab   ! (zzeta > 0) Stable 
     431            stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1 
     432            ! 
     433            psi_m_ecmwf(ji,jj) = (1._wp - stab) * psi_unst & ! (zzeta < 0) Unstable 
     434               &                +      stab  * psi_stab      ! (zzeta > 0) Stable 
    328435            ! 
    329436         END DO 
    330437      END DO 
    331       ! 
    332438   END FUNCTION psi_m_ecmwf 
    333439 
    334     
     440 
    335441   FUNCTION psi_h_ecmwf( pzeta ) 
    336442      !!---------------------------------------------------------------------------------- 
     
    342448      !!         and L is M-O length 
    343449      !! 
    344       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
     450      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    345451      !!---------------------------------------------------------------------------------- 
    346452      REAL(wp), DIMENSION(jpi,jpj) :: psi_h_ecmwf 
     
    354460         DO ji = 1, jpi 
    355461            ! 
    356             zzeta = MIN(pzeta(ji,jj) , 5.)   ! Very stable conditions (L positif and big!): 
    357             ! 
    358             zx  = ABS(1. - 16.*zzeta)**.25        ! this is actually (1/phi_m)**2  !!! 
     462            zzeta = MIN(pzeta(ji,jj) , 5._wp)   ! Very stable conditions (L positif and big!): 
     463            ! 
     464            zx  = ABS(1._wp - 16._wp*zzeta)**.25        ! this is actually (1/phi_m)**2  !!! 
    359465            !                                     ! eq.3.19, Chap.3, p.33, IFS doc - Cy31r1 
    360466            ! Unstable (Paulson 1970) : 
    361             psi_unst = 2.*LOG(0.5*(1. + zx*zx))   ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1 
     467            psi_unst = 2._wp*LOG(0.5_wp*(1._wp + zx*zx))   ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1 
    362468            ! 
    363469            ! Stable: 
    364             psi_stab = -2./3.*(zzeta - 5./0.35)*EXP(-0.35*zzeta) & ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1 
    365                &       - ABS(1. + 2./3.*zzeta)**1.5 - 2./3.*5./0.35 + 1.  
     470            psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) & ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1 
     471               &       - ABS(1._wp + 2._wp/3._wp*zzeta)**1.5_wp - 2._wp/3._wp*5._wp/0.35_wp + 1._wp 
    366472            ! LB: added ABS() to avoid NaN values when unstable, which contaminates the unstable solution... 
    367473            ! 
    368             stab = 0.5 + SIGN(0.5, zzeta) ! zzeta > 0 => stab = 1 
    369             ! 
    370             ! 
    371             psi_h_ecmwf(ji,jj) = (1. - stab) * psi_unst &   ! (zzeta < 0) Unstable 
    372                &                +    stab    * psi_stab     ! (zzeta > 0) Stable 
     474            stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1 
     475            ! 
     476            ! 
     477            psi_h_ecmwf(ji,jj) = (1._wp - stab) * psi_unst &   ! (zzeta < 0) Unstable 
     478               &                +    stab    * psi_stab        ! (zzeta > 0) Stable 
    373479            ! 
    374480         END DO 
    375481      END DO 
    376       ! 
    377482   END FUNCTION psi_h_ecmwf 
    378483 
    379  
    380    FUNCTION Ri_bulk( pz, ptz, pdt, pqz, pdq, pub ) 
    381       !!---------------------------------------------------------------------------------- 
    382       !! Bulk Richardson number (Eq. 3.25 IFS doc) 
    383       !! 
    384       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    385       !!---------------------------------------------------------------------------------- 
    386       REAL(wp), DIMENSION(jpi,jpj) ::   Ri_bulk   ! 
    387       ! 
    388       REAL(wp)                    , INTENT(in) ::   pz    ! height above the sea        [m] 
    389       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptz   ! air temperature at pz m     [K] 
    390       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pdt   ! ptz - sst                   [K] 
    391       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqz   ! air temperature at pz m [kg/kg] 
    392       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pdq   ! pqz - ssq               [kg/kg] 
    393       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pub   ! bulk wind speed           [m/s] 
    394       !!---------------------------------------------------------------------------------- 
    395       ! 
    396       Ri_bulk =   grav*pz/(pub*pub)                                          & 
    397          &      * ( pdt/(ptz - 0.5_wp*(pdt + grav*pz/(Cp_dry+Cp_vap*pqz)))   & 
    398          &          + rctv0*pdq ) 
    399       ! 
    400    END FUNCTION Ri_bulk 
    401  
    402  
    403    FUNCTION visc_air(ptak) 
    404       !!---------------------------------------------------------------------------------- 
    405       !! Air kinetic viscosity (m^2/s) given from temperature in degrees... 
    406       !! 
    407       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    408       !!---------------------------------------------------------------------------------- 
    409       REAL(wp), DIMENSION(jpi,jpj)             ::   visc_air   ! 
    410       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak       ! air temperature in (K) 
    411       ! 
    412       INTEGER  ::   ji, jj      ! dummy loop indices 
    413       REAL(wp) ::   ztc, ztc2   ! local scalar 
    414       !!---------------------------------------------------------------------------------- 
    415       ! 
    416       DO jj = 1, jpj 
    417          DO ji = 1, jpi 
    418             ztc  = ptak(ji,jj) - rt0   ! air temp, in deg. C 
    419             ztc2 = ztc*ztc 
    420             visc_air(ji,jj) = 1.326e-5*(1. + 6.542E-3*ztc + 8.301e-6*ztc2 - 4.84e-9*ztc2*ztc) 
    421          END DO 
    422       END DO 
    423       ! 
    424    END FUNCTION visc_air 
    425484 
    426485   !!====================================================================== 
Note: See TracChangeset for help on using the changeset viewer.