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.
Changeset 12377 for NEMO/trunk/src/OCE/SBC/sbcblk_algo_ecmwf.F90 – NEMO

Ignore:
Timestamp:
2020-02-12T15:39:06+01:00 (4 years ago)
Author:
acc
Message:

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.

Location:
NEMO/trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • NEMO/trunk

    • Property svn:externals

      •  

        old new  
        33^/utils/build/mk@HEAD         mk 
        44^/utils/tools@HEAD            tools 
        5 ^/vendors/AGRIF/dev@HEAD      ext/AGRIF 
         5^/vendors/AGRIF/dev_r11615_ENHANCE-04_namelists_as_internalfiles_agrif@HEAD      ext/AGRIF 
        66^/vendors/FCM@HEAD            ext/FCM 
        77^/vendors/IOIPSL@HEAD         ext/IOIPSL 

  • NEMO/trunk/src/OCE/SBC/sbcblk_algo_ecmwf.F90

    r10069 r12377  
    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   !! * Substitutions 
     48#  include "do_loop_substitute.h90" 
     49 
     50   !! ECMWF own values for given constants, taken form IFS documentation... 
    5051   REAL(wp), PARAMETER ::   charn0 = 0.018    ! Charnock constant (pretty high value here !!! 
    5152   !                                          !    =>  Usually 0.011 for moderate winds) 
    5253   REAL(wp), PARAMETER ::   zi0     = 1000.   ! scale height of the atmospheric boundary layer...1 
    5354   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] 
    5755   REAL(wp), PARAMETER ::   alpha_M = 0.11    ! For roughness length (smooth surface term) 
    5856   REAL(wp), PARAMETER ::   alpha_H = 0.40    ! (Chapter 3, p.34, IFS doc Cy31r1) 
    5957   REAL(wp), PARAMETER ::   alpha_Q = 0.62    ! 
     58 
     59   INTEGER , PARAMETER ::   nb_itt = 10             ! number of itterations 
     60 
    6061   !!---------------------------------------------------------------------- 
    6162CONTAINS 
    6263 
    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 
     64 
     65   SUBROUTINE sbcblk_algo_ecmwf_init(l_use_cs, l_use_wl) 
     66      !!--------------------------------------------------------------------- 
     67      !!                  ***  FUNCTION sbcblk_algo_ecmwf_init  *** 
    7668      !! 
    7769      !! INPUT : 
    7870      !! ------- 
     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) 
    79116      !!    *  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] 
     117      !!    *  zu   : height for wind speed (usually 10m)                     [m] 
    83118      !!    *  t_zt : potential air temperature at zt                         [K] 
    84       !!    *  ssq  : specific humidity at saturation at SST                  [kg/kg] 
    85119      !!    *  q_zt : specific humidity of air at zt                          [kg/kg] 
    86       !! 
     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] 
    87145      !! 
    88146      !! OUTPUT : 
    … …  
    93151      !!    *  t_zu   : pot. air temperature adjusted at wind height zu       [K] 
    94152      !!    *  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       !!---------------------------------------------------------------------------------- 
     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 
    100159      REAL(wp), INTENT(in   )                     ::   zt       ! height for t_zt and q_zt                    [m] 
    101160      REAL(wp), INTENT(in   )                     ::   zu       ! height for U_zu                             [m] 
    102       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   sst      ! sea surface temperature                [Kelvin] 
     161      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) ::   T_s      ! sea surface temperature                [Kelvin] 
    103162      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] 
     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] 
    106165      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 
    107168      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cd       ! transfer coefficient for momentum         (tau) 
    108169      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ch       ! transfer coefficient for sensible heat (Q_sens) 
    … …  
    110171      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   t_zu     ! pot. air temp. adjusted at zu               [K] 
    111172      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] 
     173      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   U_blk    ! bulk wind speed at zu                     [m/s] 
    113174      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients 
    114175      ! 
     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      ! 
    115183      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       ! 
     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 
    130201      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  
     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      ! 
    134221      !! 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)   !               " 
     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 )   !               " 
    137224 
    138225      !! 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 ! 
     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) 
    155242 
    156243      Cd     = (vkarmn/ztmp0)**2    ! first guess of Cd 
    157244 
    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 ) 
     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 ) 
    164251      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)) 
     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" ! 
    168255 
    169256      !! 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)) 
     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) 
    173260      t_star = dt_zu*ztmp0 
    174261      q_star = dq_zu*ztmp0 
    175262 
    176       ! What's need to be done if zt /= zu: 
     263      ! What needs to be done if zt /= zu: 
    177264      IF( .NOT. l_zt_equal_zu ) THEN 
    178          ! 
    179265         !! First update of values at zu (or zt for wind) 
    180266         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 
     267         ztmp1 = LOG(zt/zu) + ztmp0 
    182268         t_zu = t_zt - t_star/vkarmn*ztmp1 
    183269         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 ) 
    188          ! 
     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 ) 
    189274      ENDIF 
    190275 
    … …  
    194279 
    195280      !! 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 ) 
     281      Linv = One_on_L( t_zu, q_zu, u_star, t_star, q_star ) 
    198282 
    199283      !! 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  
     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) 
    205287 
    206288      !! ITERATION BLOCK 
    207       !! *************** 
    208  
    209289      DO j_itt = 1, nb_itt 
    210290 
    211291         !! Bulk Richardson Number at z=zu (Eq. 3.25) 
    212          ztmp0 = Ri_bulk(zu, t_zu, dt_zu, q_zu, dq_zu, U_blk) 
     292         ztmp0 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN) 
    213293 
    214294         !! 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 
     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...) 
    216298 
    217299         !! 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) 
     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! 
    220301 
    221302         !! Need to update roughness lengthes: 
    … …  
    223304         ztmp2  = u_star*u_star 
    224305         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 
     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 
    234314         ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0. 
    235315 
    … …  
    238318         !! as well the air-sea differences: 
    239319         IF( .NOT. l_zt_equal_zu ) THEN 
    240  
    241320            !! 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 !!! 
     321            func_h = psi_h_ecmwf(zu*Linv) ! temporary array !!! 
     322            func_m = psi_h_ecmwf(zt*Linv) ! temporary array !!! 
    244323 
    245324            ztmp2  = psi_h_ecmwf(z0t*Linv) 
    246325            ztmp0  = func_h - ztmp2 
    247             ztmp1  = vkarmn/(LOG(zu+z0) - LOG(z0t) - ztmp0) 
     326            ztmp1  = vkarmn/(LOG(zu) - LOG(z0t) - ztmp0) 
    248327            t_star = dt_zu*ztmp1 
    249328            ztmp2  = ztmp0 - func_m + ztmp2 
    … …  
    253332            ztmp2  = psi_h_ecmwf(z0q*Linv) 
    254333            ztmp0  = func_h - ztmp2 
    255             ztmp1  = vkarmn/(LOG(zu+z0) - LOG(z0q) - ztmp0) 
     334            ztmp1  = vkarmn/(LOG(zu) - LOG(z0q) - ztmp0) 
    256335            q_star = dq_zu*ztmp1 
    257336            ztmp2  = ztmp0 - func_m + ztmp2 
    258             ztmp1  = log(zt/zu) + ztmp2 
     337            ztmp1  = LOG(zt/zu) + ztmp2 
    259338            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 
     339         ENDIF 
    265340 
    266341         !! 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 
     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 
    273378 
    274379      Cd = vkarmn*vkarmn/(func_m*func_m) 
    275380      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 
     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 
    285395 
    286396 
    … …  
    294404      !!         and L is M-O length 
    295405      !! 
    296       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
     406      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    297407      !!---------------------------------------------------------------------------------- 
    298408      REAL(wp), DIMENSION(jpi,jpj) :: psi_m_ecmwf 
    … …  
    302412      REAL(wp) :: zzeta, zx, ztmp, psi_unst, psi_stab, stab 
    303413      !!---------------------------------------------------------------------------------- 
    304       ! 
    305       DO jj = 1, jpj 
    306          DO ji = 1, jpi 
    307             ! 
    308             zzeta = MIN( pzeta(ji,jj) , 5. ) !! Very stable conditions (L positif and big!): 
    309             ! 
    310             ! Unstable (Paulson 1970): 
    311             !   eq.3.20, Chap.3, p.33, IFS doc - Cy31r1 
    312             zx = SQRT(ABS(1. - 16.*zzeta)) 
    313             ztmp = 1. + SQRT(zx) 
    314             ztmp = ztmp*ztmp 
    315             psi_unst = LOG( 0.125*ztmp*(1. + zx) )   & 
    316                &       -2.*ATAN( SQRT(zx) ) + 0.5*rpi 
    317             ! 
    318             ! Unstable: 
    319             ! 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 
    322             ! 
    323             ! 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 
    328             ! 
    329          END DO 
    330       END DO 
    331       ! 
     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 
    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 
    … …  
    351457      !!---------------------------------------------------------------------------------- 
    352458      ! 
    353       DO jj = 1, jpj 
    354          DO ji = 1, jpi 
    355             ! 
    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  !!! 
    359             !                                     ! eq.3.19, Chap.3, p.33, IFS doc - Cy31r1 
    360             ! Unstable (Paulson 1970) : 
    361             psi_unst = 2.*LOG(0.5*(1. + zx*zx))   ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1 
    362             ! 
    363             ! 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.  
    366             ! LB: added ABS() to avoid NaN values when unstable, which contaminates the unstable solution... 
    367             ! 
    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 
    373             ! 
    374          END DO 
    375       END DO 
    376       ! 
     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 
    377480   END FUNCTION psi_h_ecmwf 
    378481 
    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 
    425482 
    426483   !!====================================================================== 
Note: See TracChangeset for help on using the changeset viewer.