Changeset 12015


Ignore:
Timestamp:
2019-11-29T16:59:07+01:00 (8 months ago)
Author:
gsamson
Message:

dev_ASINTER-01-05_merged: merge dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk branch @rev11988 with dev_r11265_ASINTER-01_Guillaume_ABL1D branch @rev11937 (tickets #2159 and #2131); ORCA2_ICE(_ABL) reproductibility OK

Location:
NEMO/branches/2019/dev_ASINTER-01-05_merged/src
Files:
5 added
2 deleted
8 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_ASINTER-01-05_merged/src/ABL/ablmod.F90

    r11937 r12015  
    1717   USE phycst         ! physical constants 
    1818   USE dom_oce, ONLY  : tmask   
    19    USE sbc_oce, ONLY  : ght_abl, ghw_abl, e3t_abl, e3w_abl, jpka, jpkam1 
    20    USE sbcblk         ! use some physical constants for flux computation 
     19   USE sbc_oce, ONLY  : ght_abl, ghw_abl, e3t_abl, e3w_abl, jpka, jpkam1, rhoa 
     20   USE sbcblk         ! use rn_?fac 
     21   USE sbcblk_phy     ! use some physical constants for flux computation 
    2122   ! 
    2223   USE prtctl         ! Print control                    (prt_ctl routine) 
     
    100101      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptauj_ice    ! ice-surface tauy stress (V-point)      
    101102#endif      
    102      ! 
    103       REAL(wp), DIMENSION(1:jpi,1:jpj   )        ::   zrhoa, zwnd_i, zwnd_j 
     103      ! 
     104      REAL(wp), DIMENSION(1:jpi,1:jpj   )        ::   zwnd_i, zwnd_j 
    104105      REAL(wp), DIMENSION(1:jpi,2:jpka  )        ::   zCF     
    105106      REAL(wp), DIMENSION(1:jpi,1:jpj,1:jpka)    ::   z_cft      !--FL--to be removed after the test phase    
     
    529530            ztemp             = tq_abl  ( ji, jj, 2, nt_a, jp_ta )  
    530531            zhumi             = tq_abl  ( ji, jj, 2, nt_a, jp_qa )  
    531             zcff              = pslp_dta( ji, jj ) /   &              !<-- At this point ztemp and zhumi should not be zero ... 
    532                &                        (  R_dry*ztemp * ( 1._wp + rctv0*zhumi )  ) 
     532            !zcff              = pslp_dta( ji, jj ) /   &              !<-- At this point ztemp and zhumi should not be zero ... 
     533            !   &                        (  R_dry*ztemp * ( 1._wp + rctv0*zhumi )  ) 
     534            zcff              = rho_air( ztemp, zhumi, pslp_dta( ji, jj ) ) 
    533535            psen ( ji, jj )   =      cp_air(zhumi) * zcff * psen(ji,jj) * ( psst(ji,jj) + rt0 - ztemp ) 
    534536            pevp ( ji, jj )   = rn_efac*MAX( 0._wp,  zcff * pevp(ji,jj) * ( pssq(ji,jj)       - zhumi ) ) 
    535             zrhoa( ji, jj )   = zcff               
     537            rhoa( ji, jj )   = zcff               
    536538         END DO 
    537539      END DO 
     
    551553            zcff          = SQRT(  zwnd_i(ji,jj) * zwnd_i(ji,jj)   & 
    552554               &                 + zwnd_j(ji,jj) * zwnd_j(ji,jj)  )  ! * msk_abl(ji,jj) 
    553             zztmp         = zrhoa(ji,jj) * pcd_du(ji,jj) 
     555            zztmp         = rhoa(ji,jj) * pcd_du(ji,jj) 
    554556             
    555557            pwndm (ji,jj) =         zcff 
     
    593595               zztmp2 = 0.5_wp * ( v_abl(ji,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) ) 
    594596             
    595                ptaui_ice(ji,jj) = 0.5_wp * (  zrhoa(ji+1,jj) * pCd_du_ice(ji+1,jj)             & 
    596                   &                      +    zrhoa(ji  ,jj) * pCd_du_ice(ji  ,jj)  )          & 
     597               ptaui_ice(ji,jj) = 0.5_wp * (  rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj)             & 
     598                  &                      +    rhoa(ji  ,jj) * pCd_du_ice(ji  ,jj)  )          & 
    597599                  &         * ( zztmp1 - rn_vfac * pssu_ice(ji,jj) ) 
    598                ptauj_ice(ji,jj) = 0.5_wp * (  zrhoa(ji,jj+1) * pCd_du_ice(ji,jj+1)             & 
    599                   &                      +    zrhoa(ji,jj  ) * pCd_du_ice(ji,jj  )  )          & 
     600               ptauj_ice(ji,jj) = 0.5_wp * (  rhoa(ji,jj+1) * pCd_du_ice(ji,jj+1)             & 
     601                  &                      +    rhoa(ji,jj  ) * pCd_du_ice(ji,jj  )  )          & 
    600602                  &         * ( zztmp2 - rn_vfac * pssv_ice(ji,jj) ) 
    601603            END DO 
  • NEMO/branches/2019/dev_ASINTER-01-05_merged/src/ABL/sbcabl.F90

    r11858 r12015  
    341341         &                tq_abl(:,:,2,nt_n,jp_ta), tq_abl(:,:,2,nt_n,jp_qa),   &   !   <<= in 
    342342         &                sf(jp_slp )%fnow(:,:,1) , sst_m, ssu_m, ssv_m     ,   &   !   <<= in 
     343         &                sf(jp_qsr )%fnow(:,:,1) , sf(jp_qlw )%fnow(:,:,1) ,   &   !   <<= in 
    343344         &                zssq, zcd_du, zsen, zevp                          )       !   =>> out 
    344345   
  • NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/SBC/sbc_oce.F90

    r11275 r12015  
    116116   !! wndm is used compute surface gases exchanges in ice-free ocean or leads 
    117117   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wndm              !: wind speed module at T-point (=|U10m-Uoce|)  [m/s] 
     118   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rhoa              !: air density at "rn_zu" m above the sea       [kg/m3] !LB 
    118119   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qsr               !: sea heat flux:     solar                     [W/m2] 
    119120   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qns    , qns_b    !: sea heat flux: non solar                     [W/m2] 
     
    159160   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   frq_m     !: mean (nn_fsbc time-step) fraction of solar net radiation absorbed in the 1st T level [-] 
    160161 
     162   !!---------------------------------------------------------------------- 
     163   !!                     Cool-skin/Warm-layer 
     164   !!---------------------------------------------------------------------- 
     165   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tsk       !: sea-surface skin temperature (used if ln_skin_cs==.true. .OR. ln_skin_wl==.true.)  [K] !LB 
     166 
     167    
    161168   !! * Substitutions 
    162169#  include "vectopt_loop_substitute.h90" 
     
    177184      ! 
    178185      ALLOCATE( utau(jpi,jpj) , utau_b(jpi,jpj) , taum(jpi,jpj) ,     & 
    179          &      vtau(jpi,jpj) , vtau_b(jpi,jpj) , wndm(jpi,jpj) , STAT=ierr(1) )  
     186         &      vtau(jpi,jpj) , vtau_b(jpi,jpj) , wndm(jpi,jpj) , rhoa(jpi,jpj) , STAT=ierr(1) )  
    180187         ! 
    181188      ALLOCATE( qns_tot(jpi,jpj) , qns  (jpi,jpj) , qns_b(jpi,jpj),        & 
  • NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/SBC/sbcblk.F90

    r11586 r12015  
    1515   !!            3.7  !  2014-06  (L. Brodeau)  simplification and optimization of CORE bulk 
    1616   !!            4.0  !  2016-06  (L. Brodeau)  sbcblk_core becomes sbcblk and is not restricted to the CORE algorithm anymore 
    17    !!                 !                        ==> based on AeroBulk (http://aerobulk.sourceforge.net/) 
     17   !!                 !                        ==> based on AeroBulk (https://github.com/brodeau/aerobulk/) 
    1818   !!            4.0  !  2016-10  (G. Madec)  introduce a sbc_blk_init routine 
    19    !!            4.0  !  2016-10  (M. Vancoppenolle)  Introduce conduction flux emulator (M. Vancoppenolle)  
     19   !!            4.0  !  2016-10  (M. Vancoppenolle)  Introduce conduction flux emulator (M. Vancoppenolle) 
    2020   !!            4.0  !  2019-03  (F. Lemarié & G. Samson)  add ABL compatibility (ln_abl=TRUE) 
    2121   !!---------------------------------------------------------------------- 
     
    2424   !!   sbc_blk_init  : initialisation of the chosen bulk formulation as ocean surface boundary condition 
    2525   !!   sbc_blk       : bulk formulation as ocean surface boundary condition 
    26    !!       blk_oce_1 : computes pieces of momentum, heat and freshwater fluxes over ocean for ABL model  (ln_abl=TRUE)   
    27    !!       blk_oce_2 : finalizes momentum, heat and freshwater fluxes computation over ocean after the ABL step  (ln_abl=TRUE)  
    28    !!   rho_air       : density of (moist) air (depends on T_air, q_air and SLP 
    29    !!   cp_air        : specific heat of (moist) air (depends spec. hum. q_air) 
    30    !!   q_sat         : saturation humidity as a function of SLP and temperature 
    31    !!   L_vap         : latent heat of vaporization of water as a function of temperature 
    32    !!             sea-ice case only :  
     26   !!   blk_oce_1     : computes pieces of momentum, heat and freshwater fluxes over ocean for ABL model  (ln_abl=TRUE)   
     27   !!   blk_oce_2     : finalizes momentum, heat and freshwater fluxes computation over ocean after the ABL step  (ln_abl=TRUE)  
     28   !!             sea-ice case only : 
    3329   !!   blk_ice_1   : provide the air-ice stress 
    3430   !!   blk_ice_2   : provide the heat and mass fluxes at air-ice interface 
    3531   !!   blk_ice_qcn   : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux) 
    3632   !!   Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag 
    37    !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag  
     33   !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag 
    3834   !!---------------------------------------------------------------------- 
    3935   USE oce            ! ocean dynamics and tracers 
     
    5147   USE icethd_dh      ! for CALL ice_thd_snwblow 
    5248#endif 
    53    USE sbcblk_algo_ncar     ! => turb_ncar     : NCAR - CORE (Large & Yeager, 2009)  
    54    USE sbcblk_algo_coare    ! => turb_coare    : COAREv3.0 (Fairall et al. 2003)  
    55    USE sbcblk_algo_coare3p5 ! => turb_coare3p5 : COAREv3.5 (Edson et al. 2013) 
    56    USE sbcblk_algo_ecmwf    ! => turb_ecmwf    : ECMWF (IFS cycle 31)  
     49   USE sbcblk_algo_ncar     ! => turb_ncar     : NCAR - CORE (Large & Yeager, 2009) 
     50   USE sbcblk_algo_coare3p0 ! => turb_coare3p0 : COAREv3.0 (Fairall et al. 2003) 
     51   USE sbcblk_algo_coare3p6 ! => turb_coare3p6 : COAREv3.6 (Fairall et al. 2018 + Edson et al. 2013) 
     52   USE sbcblk_algo_ecmwf    ! => turb_ecmwf    : ECMWF (IFS cycle 45r1) 
    5753   ! 
    5854   USE iom            ! I/O manager library 
     
    6258   USE prtctl         ! Print control 
    6359 
     60   USE sbcblk_phy     !LB: all thermodynamics functions in the marine boundary layer, rho_air, q_sat, etc... 
     61 
     62 
    6463   IMPLICIT NONE 
    6564   PRIVATE 
     
    6968   PUBLIC   blk_oce_1     ! called in sbcabl 
    7069   PUBLIC   blk_oce_2     ! called in sbcabl 
    71    PUBLIC   rho_air       ! called in ablmod 
    72    PUBLIC   cp_air        ! called in ablmod 
    7370#if defined key_si3 
    7471   PUBLIC   blk_ice_1     ! routine called in icesbc 
    7572   PUBLIC   blk_ice_2     ! routine called in icesbc 
    7673   PUBLIC   blk_ice_qcn   ! routine called in icesbc 
    77 #endif  
    78  
    79   INTERFACE cp_air 
    80     MODULE PROCEDURE cp_air_0d, cp_air_2d 
    81   END INTERFACE 
    82  
    83    !                                   !!* Namelist namsbc_blk : bulk parameters 
    84    LOGICAL          ::   ln_NCAR        ! "NCAR"      algorithm   (Large and Yeager 2008) 
    85    LOGICAL          ::   ln_COARE_3p0   ! "COARE 3.0" algorithm   (Fairall et al. 2003) 
    86    LOGICAL          ::   ln_COARE_3p5   ! "COARE 3.5" algorithm   (Edson et al. 2013) 
    87    LOGICAL          ::   ln_ECMWF       ! "ECMWF"     algorithm   (IFS cycle 31) 
    88    !                                    ! 
    89    REAL(wp)         ::   rn_pfac        !  multiplication factor for precipitation 
    90    REAL(wp), PUBLIC ::   rn_efac        !: multiplication factor for evaporation 
    91    REAL(wp), PUBLIC ::   rn_vfac        !: multiplication factor for ice/ocean velocity in the calculation of wind stress 
    92    REAL(wp)         ::   rn_zqt         !  z(q,t) : height of humidity and temperature measurements 
    93    REAL(wp)         ::   rn_zu          !  z(u)   : height of wind measurements 
    94    !                                    ! 
    95    LOGICAL          ::   ln_Cd_L12      ! ice-atm drag = F( ice concentration )                        (Lupkes et al. JGR2012) 
    96    LOGICAL          ::   ln_Cd_L15      ! ice-atm drag = F( ice concentration, atmospheric stability ) (Lupkes et al. JGR2015) 
    97  
    98    INTEGER  ::   nblk                   ! choice of the bulk algorithm 
    99    !                                       ! associated indices: 
    100    INTEGER, PARAMETER ::   np_NCAR      = 1   ! "NCAR" algorithm        (Large and Yeager 2008) 
    101    INTEGER, PARAMETER ::   np_COARE_3p0 = 2   ! "COARE 3.0" algorithm   (Fairall et al. 2003) 
    102    INTEGER, PARAMETER ::   np_COARE_3p5 = 3   ! "COARE 3.5" algorithm   (Edson et al. 2013) 
    103    INTEGER, PARAMETER ::   np_ECMWF     = 4   ! "ECMWF" algorithm       (IFS cycle 31) 
    104  
    105    !                                                      !!! air parameters 
    106    REAL(wp)        , PARAMETER ::   Cp_dry = 1005.0        !  Specic heat of dry air, constant pressure      [J/K/kg] 
    107    REAL(wp)        , PARAMETER ::   Cp_vap = 1860.0        !  Specic heat of water vapor, constant pressure  [J/K/kg] 
    108    REAL(wp), PUBLIC, PARAMETER ::   R_dry = 287.05_wp     !: Specific gas constant for dry air              [J/K/kg] 
    109    REAL(wp)        , PARAMETER ::   R_vap  = 461.495_wp    !  Specific gas constant for water vapor          [J/K/kg] 
    110    REAL(wp)        , PARAMETER ::   reps0  = R_dry/R_vap   !  ratio of gas constant for dry air and water vapor => ~ 0.622 
    111    REAL(wp), PUBLIC, PARAMETER ::   rctv0 = R_vap/R_dry   !: for virtual temperature (== (1-eps)/eps) => ~ 0.608 
    112    !                                                      !!! Bulk parameters 
    113    REAL(wp)        , PARAMETER ::   cpa    = 1000.5        ! specific heat of air (only used for ice fluxes now...) 
    114    REAL(wp)        , PARAMETER ::   Ls     =    2.839e6    ! latent heat of sublimation 
    115    REAL(wp)        , PARAMETER ::   Stef   =    5.67e-8    ! Stefan Boltzmann constant 
    116    REAL(wp)        , PARAMETER ::   rcdice =    1.4e-3     ! transfer coefficient over ice 
    117    REAL(wp)        , PARAMETER ::   albo   =    0.066      ! ocean albedo assumed to be constant 
     74#endif 
     75 
     76   INTEGER , PUBLIC            ::   jpfld         ! maximum number of files to read 
     77   INTEGER , PUBLIC, PARAMETER ::   jp_wndi = 1   ! index of 10m wind velocity (i-component) (m/s)    at T-point 
     78   INTEGER , PUBLIC, PARAMETER ::   jp_wndj = 2   ! index of 10m wind velocity (j-component) (m/s)    at T-point 
     79   INTEGER , PUBLIC, PARAMETER ::   jp_tair = 3   ! index of 10m air temperature             (Kelvin) 
     80   INTEGER , PUBLIC, PARAMETER ::   jp_humi = 4   ! index of specific humidity               ( % ) 
     81   INTEGER , PUBLIC, PARAMETER ::   jp_qsr  = 5   ! index of solar heat                      (W/m2) 
     82   INTEGER , PUBLIC, PARAMETER ::   jp_qlw  = 6   ! index of Long wave                       (W/m2) 
     83   INTEGER , PUBLIC, PARAMETER ::   jp_prec = 7   ! index of total precipitation (rain+snow) (Kg/m2/s) 
     84   INTEGER , PUBLIC, PARAMETER ::   jp_snow = 8   ! index of snow (solid prcipitation)       (kg/m2/s) 
     85   INTEGER , PUBLIC, PARAMETER ::   jp_slp  = 9   ! index of sea level pressure              (Pa) 
     86   INTEGER , PUBLIC, PARAMETER ::   jp_hpgi =10   ! index of ABL geostrophic wind or hpg (i-component) (m/s) at T-point 
     87   INTEGER , PUBLIC, PARAMETER ::   jp_hpgj =11   ! index of ABL geostrophic wind or hpg (j-component) (m/s) at T-point 
     88 
     89   TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) ::   sf   ! structure of input atmospheric fields (file informations, fields read) 
     90 
     91   !                           !!* Namelist namsbc_blk : bulk parameters 
     92   LOGICAL  ::   ln_NCAR        ! "NCAR"      algorithm   (Large and Yeager 2008) 
     93   LOGICAL  ::   ln_COARE_3p0   ! "COARE 3.0" algorithm   (Fairall et al. 2003) 
     94   LOGICAL  ::   ln_COARE_3p6   ! "COARE 3.6" algorithm   (Edson et al. 2013) 
     95   LOGICAL  ::   ln_ECMWF       ! "ECMWF"     algorithm   (IFS cycle 45r1) 
     96   ! 
     97   REAL(wp)         ::   rn_pfac   ! multiplication factor for precipitation 
     98   REAL(wp), PUBLIC ::   rn_efac   ! multiplication factor for evaporation 
     99   REAL(wp), PUBLIC ::   rn_vfac   ! multiplication factor for ice/ocean velocity in the calculation of wind stress 
     100   REAL(wp)         ::   rn_zqt    ! z(q,t) : height of humidity and temperature measurements 
     101   REAL(wp)         ::   rn_zu     ! z(u)   : height of wind measurements 
     102   ! 
     103   LOGICAL  ::   ln_Cd_L12      ! ice-atm drag = F( ice concentration )                        (Lupkes et al. JGR2012) 
     104   LOGICAL  ::   ln_Cd_L15      ! ice-atm drag = F( ice concentration, atmospheric stability ) (Lupkes et al. JGR2015) 
    118105   ! 
    119106   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   Cd_ice , Ch_ice , Ce_ice   ! transfert coefficients over ice 
     
    121108   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   t_zu, q_zu                 ! air temp. and spec. hum. at wind speed height (L15 bulk scheme) 
    122109 
    123    !INTEGER , PUBLIC, PARAMETER ::   jpfld   =11   !: maximum number of files to read 
    124    INTEGER , PUBLIC            ::   jpfld         !: maximum number of files to read 
    125    INTEGER , PUBLIC, PARAMETER ::   jp_wndi = 1   !: index of 10m wind velocity (i-component) (m/s)    at T-point 
    126    INTEGER , PUBLIC, PARAMETER ::   jp_wndj = 2   !: index of 10m wind velocity (j-component) (m/s)    at T-point 
    127    INTEGER , PUBLIC, PARAMETER ::   jp_tair = 3   !: index of 10m air temperature             (Kelvin) 
    128    INTEGER , PUBLIC, PARAMETER ::   jp_humi = 4   !: index of specific humidity               ( % ) 
    129    INTEGER , PUBLIC, PARAMETER ::   jp_qsr  = 5   !: index of solar heat                      (W/m2) 
    130    INTEGER , PUBLIC, PARAMETER ::   jp_qlw  = 6   !: index of Long wave                       (W/m2) 
    131    INTEGER , PUBLIC, PARAMETER ::   jp_prec = 7   !: index of total precipitation (rain+snow) (Kg/m2/s) 
    132    INTEGER , PUBLIC, PARAMETER ::   jp_snow = 8   !: index of snow (solid prcipitation)       (kg/m2/s) 
    133    INTEGER , PUBLIC, PARAMETER ::   jp_slp  = 9   !: index of sea level pressure              (Pa) 
    134    INTEGER , PUBLIC, PARAMETER ::   jp_hpgi =10   !: index of ABL geostrophic wind or hpg (i-component) (m/s) at T-point 
    135    INTEGER , PUBLIC, PARAMETER ::   jp_hpgj =11   !: index of ABL geostrophic wind or hpg (j-component) (m/s) at T-point 
    136  
    137    TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) ::   sf   !  structure of input atmospheric fields (file informations, fields read) 
     110   LOGICAL  ::   ln_skin_cs     ! use the cool-skin (only available in ECMWF and COARE algorithms) !LB 
     111   LOGICAL  ::   ln_skin_wl     ! use the warm-layer parameterization (only available in ECMWF and COARE algorithms) !LB 
     112   LOGICAL  ::   ln_humi_sph    ! humidity read in files ("sn_humi") is specific humidity [kg/kg] if .true. !LB 
     113   LOGICAL  ::   ln_humi_dpt    ! humidity read in files ("sn_humi") is dew-point temperature [K] if .true. !LB 
     114   LOGICAL  ::   ln_humi_rlh    ! humidity read in files ("sn_humi") is relative humidity     [%] if .true. !LB 
     115   ! 
     116   INTEGER  ::   nhumi          ! choice of the bulk algorithm 
     117   !                            ! associated indices: 
     118   INTEGER, PARAMETER :: np_humi_sph = 1 
     119   INTEGER, PARAMETER :: np_humi_dpt = 2 
     120   INTEGER, PARAMETER :: np_humi_rlh = 3 
     121 
     122   INTEGER  ::   nblk           ! choice of the bulk algorithm 
     123   !                            ! associated indices: 
     124   INTEGER, PARAMETER ::   np_NCAR      = 1   ! "NCAR" algorithm        (Large and Yeager 2008) 
     125   INTEGER, PARAMETER ::   np_COARE_3p0 = 2   ! "COARE 3.0" algorithm   (Fairall et al. 2003) 
     126   INTEGER, PARAMETER ::   np_COARE_3p6 = 3   ! "COARE 3.6" algorithm   (Edson et al. 2013) 
     127   INTEGER, PARAMETER ::   np_ECMWF     = 4   ! "ECMWF" algorithm       (IFS cycle 45r1) 
    138128 
    139129   !! * Substitutions 
     
    165155      !! ** Purpose :   choose and initialize a bulk formulae formulation 
    166156      !! 
    167       !! ** Method  :  
     157      !! ** Method  : 
    168158      !! 
    169159      !!---------------------------------------------------------------------- 
     
    172162      !! 
    173163      CHARACTER(len=100)            ::   cn_dir                ! Root directory for location of atmospheric forcing files 
    174       !TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i                 ! array of namelist informations on the fields to read 
    175       TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   slf_i                 ! array of namelist informations on the fields to read 
     164      TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   slf_i        ! array of namelist informations on the fields to read 
    176165      TYPE(FLD_N) ::   sn_wndi, sn_wndj, sn_humi, sn_qsr       ! informations about the fields to be read 
    177166      TYPE(FLD_N) ::   sn_qlw , sn_tair, sn_prec, sn_snow      !       "                        " 
     
    179168      NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw ,                &   ! input fields 
    180169         &                 sn_tair, sn_prec, sn_snow, sn_slp, sn_hpgi, sn_hpgj,       & 
    181          &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF,             &   ! bulk algorithm 
    182          &                 cn_dir , rn_zqt, rn_zu,                                    &  
    183          &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15 
     170         &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p6, ln_ECMWF,             &   ! bulk algorithm 
     171         &                 cn_dir , rn_zqt, rn_zu,                                    & 
     172         &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15,           & 
     173         &                 ln_skin_cs, ln_skin_wl, ln_humi_sph, ln_humi_dpt, ln_humi_rlh  ! cool-skin / warm-layer !LB 
    184174      !!--------------------------------------------------------------------- 
    185175      ! 
     
    187177      IF( sbc_blk_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard arrays' ) 
    188178      ! 
    189       !                             !** read bulk namelist   
     179      !                             !** read bulk namelist 
    190180      REWIND( numnam_ref )                !* Namelist namsbc_blk in reference namelist : bulk parameters 
    191181      READ  ( numnam_ref, namsbc_blk, IOSTAT = ios, ERR = 901) 
     
    200190      !                             !** initialization of the chosen bulk formulae (+ check) 
    201191      !                                   !* select the bulk chosen in the namelist and check the choice 
    202                                                                ioptio = 0 
    203       IF( ln_NCAR      ) THEN   ;   nblk =  np_NCAR        ;   ioptio = ioptio + 1   ;   ENDIF 
    204       IF( ln_COARE_3p0 ) THEN   ;   nblk =  np_COARE_3p0   ;   ioptio = ioptio + 1   ;   ENDIF 
    205       IF( ln_COARE_3p5 ) THEN   ;   nblk =  np_COARE_3p5   ;   ioptio = ioptio + 1   ;   ENDIF 
    206       IF( ln_ECMWF     ) THEN   ;   nblk =  np_ECMWF       ;   ioptio = ioptio + 1   ;   ENDIF 
    207       ! 
     192      ioptio = 0 
     193      IF( ln_NCAR      ) THEN 
     194         nblk =  np_NCAR        ;   ioptio = ioptio + 1 
     195      ENDIF 
     196      IF( ln_COARE_3p0 ) THEN 
     197         nblk =  np_COARE_3p0   ;   ioptio = ioptio + 1 
     198      ENDIF 
     199      IF( ln_COARE_3p6 ) THEN 
     200         nblk =  np_COARE_3p6   ;   ioptio = ioptio + 1 
     201      ENDIF 
     202      IF( ln_ECMWF     ) THEN 
     203         nblk =  np_ECMWF       ;   ioptio = ioptio + 1 
     204      ENDIF 
    208205      IF( ioptio /= 1 )   CALL ctl_stop( 'sbc_blk_init: Choose one and only one bulk algorithm' ) 
     206 
     207      !                             !** initialization of the cool-skin / warm-layer parametrization 
     208      IF( ln_NCAR .AND. (ln_skin_cs .OR. ln_skin_wl) ) & 
     209         & CALL ctl_stop( 'sbc_blk_init: Cool-skin/warm-layer param. not compatible with NCAR algorithm!' ) 
     210      ! 
     211      ioptio = 0 
     212      IF( ln_humi_sph ) THEN 
     213         nhumi =  np_humi_sph    ;   ioptio = ioptio + 1 
     214      ENDIF 
     215      IF( ln_humi_dpt ) THEN 
     216         nhumi =  np_humi_dpt    ;   ioptio = ioptio + 1 
     217      ENDIF 
     218      IF( ln_humi_rlh ) THEN 
     219         nhumi =  np_humi_rlh    ;   ioptio = ioptio + 1 
     220      ENDIF 
     221      IF( ioptio /= 1 )   CALL ctl_stop( 'sbc_blk_init: Choose one and only one type of air humidity' ) 
    209222      ! 
    210223      IF( ln_dm2dc ) THEN                 !* check: diurnal cycle on Qsr 
    211224         IF( sn_qsr%freqh /= 24. )   CALL ctl_stop( 'sbc_blk_init: ln_dm2dc=T only with daily short-wave input' ) 
    212          IF( sn_qsr%ln_tint ) THEN  
     225         IF( sn_qsr%ln_tint ) THEN 
    213226            CALL ctl_warn( 'sbc_blk_init: ln_dm2dc=T daily qsr time interpolation done by sbcdcy module',   & 
    214227               &           '              ==> We force time interpolation = .false. for qsr' ) 
     
    234247      ALLOCATE( sf(jpfld), STAT=ierror ) 
    235248      IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_init: unable to allocate sf structure' ) 
    236       !                                      !- fill the bulk structure with namelist informations 
    237       CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_init', 'surface boundary condition -- bulk formulae', 'namsbc_blk' ) 
    238249      ! 
    239250      DO jfpr= 1, jpfld 
     
    258269         ENDIF 
    259270      END DO 
    260       ! 
    261       IF( ln_wave ) THEN      ! surface waves 
    262          IF( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) )   &   ! Activated wave module but neither drag nor stokes drift activated 
    263             &   CALL ctl_stop( 'sbc_blk_init: Ask for wave coupling but ln_cdgw=ln_sdw=ln_tauwoc=ln_stcor=F' ) 
    264          IF( ln_cdgw .AND. .NOT.ln_NCAR )                                 &   ! drag coefficient read from wave model only with NCAR bulk formulae  
    265             &   CALL ctl_stop( 'sbc_blk_init: drag coefficient read from wave model need NCAR bulk formulae') 
    266          IF( ln_stcor .AND. .NOT.ln_sdw )                                 & 
    267             CALL ctl_stop( 'sbc_blk_init: Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T') 
     271      !                                      !- fill the bulk structure with namelist informations 
     272      CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_init', 'surface boundary condition -- bulk formulae', 'namsbc_blk' ) 
     273      ! 
     274      IF( ln_wave ) THEN 
     275         !Activated wave module but neither drag nor stokes drift activated 
     276         IF( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) )   THEN 
     277            CALL ctl_stop( 'STOP',  'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauwoc=F, ln_stcor=F' ) 
     278            !drag coefficient read from wave model definable only with mfs bulk formulae and core 
     279         ELSEIF(ln_cdgw .AND. .NOT. ln_NCAR )       THEN 
     280            CALL ctl_stop( 'drag coefficient read from wave model definable only with NCAR and CORE bulk formulae') 
     281         ELSEIF(ln_stcor .AND. .NOT. ln_sdw)                             THEN 
     282            CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T') 
     283         ENDIF 
    268284      ELSE 
    269          IF( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) THEN 
    270             CALL ctl_warn( 'sbc_blk_init: ln_wave=F, set all wave-related namelist parameter to FALSE') 
    271             ln_cdgw =.FALSE.   ;   ln_sdw =.FALSE.   ;   ln_tauwoc =.FALSE.   ;   ln_stcor =.FALSE.    
    272          ENDIF 
    273       ENDIF  
     285         IF( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor )                & 
     286            &   CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ',    & 
     287            &                  'with drag coefficient (ln_cdgw =T) '  ,                        & 
     288            &                  'or Stokes Drift (ln_sdw=T) ' ,                                 & 
     289            &                  'or ocean stress modification due to waves (ln_tauwoc=T) ',      & 
     290            &                  'or Stokes-Coriolis term (ln_stcori=T)'  ) 
     291      ENDIF 
    274292      ! 
    275293      IF( ln_abl ) THEN       ! ABL: read 3D fields for wind, temperature, humidity and pressure gradient     
     
    281299      ! 
    282300      ! set transfer coefficients to default sea-ice values 
    283       Cd_ice(:,:) = rcdice 
    284       Ch_ice(:,:) = rcdice 
    285       Ce_ice(:,:) = rcdice 
     301      Cd_ice(:,:) = rCd_ice 
     302      Ch_ice(:,:) = rCd_ice 
     303      Ce_ice(:,:) = rCd_ice 
    286304      ! 
    287305      IF(lwp) THEN                     !** Control print 
    288306         ! 
    289          WRITE(numout,*)                  !* namelist  
     307         WRITE(numout,*)                  !* namelist 
    290308         WRITE(numout,*) '   Namelist namsbc_blk (other than data information):' 
    291309         WRITE(numout,*) '      "NCAR"      algorithm   (Large and Yeager 2008)     ln_NCAR      = ', ln_NCAR 
    292310         WRITE(numout,*) '      "COARE 3.0" algorithm   (Fairall et al. 2003)       ln_COARE_3p0 = ', ln_COARE_3p0 
    293          WRITE(numout,*) '      "COARE 3.5" algorithm   (Edson et al. 2013)         ln_COARE_3p5 = ', ln_COARE_3p0 
    294          WRITE(numout,*) '      "ECMWF"     algorithm   (IFS cycle 31)              ln_ECMWF     = ', ln_ECMWF 
     311         WRITE(numout,*) '      "COARE 3.6" algorithm (Fairall 2018 + Edson al 2013)ln_COARE_3p6 = ', ln_COARE_3p6 
     312         WRITE(numout,*) '      "ECMWF"     algorithm   (IFS cycle 45r1)            ln_ECMWF     = ', ln_ECMWF 
    295313         WRITE(numout,*) '      Air temperature and humidity reference height (m)   rn_zqt       = ', rn_zqt 
    296314         WRITE(numout,*) '      Wind vector reference height (m)                    rn_zu        = ', rn_zu 
     
    306324         CASE( np_NCAR      )   ;   WRITE(numout,*) '   ==>>>   "NCAR" algorithm        (Large and Yeager 2008)' 
    307325         CASE( np_COARE_3p0 )   ;   WRITE(numout,*) '   ==>>>   "COARE 3.0" algorithm   (Fairall et al. 2003)' 
    308          CASE( np_COARE_3p5 )   ;   WRITE(numout,*) '   ==>>>   "COARE 3.5" algorithm   (Edson et al. 2013)' 
    309          CASE( np_ECMWF     )   ;   WRITE(numout,*) '   ==>>>   "ECMWF" algorithm       (IFS cycle 31)' 
     326         CASE( np_COARE_3p6 )   ;   WRITE(numout,*) '   ==>>>   "COARE 3.6" algorithm (Fairall 2018+Edson et al. 2013)' 
     327         CASE( np_ECMWF     )   ;   WRITE(numout,*) '   ==>>>   "ECMWF" algorithm       (IFS cycle 45r1)' 
     328         END SELECT 
     329         ! 
     330         WRITE(numout,*) 
     331         WRITE(numout,*) '      use cool-skin  parameterization (SSST)  ln_skin_cs  = ', ln_skin_cs 
     332         WRITE(numout,*) '      use warm-layer parameterization (SSST)  ln_skin_wl  = ', ln_skin_wl 
     333         ! 
     334         WRITE(numout,*) 
     335         SELECT CASE( nhumi )              !* Print the choice of air humidity 
     336         CASE( np_humi_sph )   ;   WRITE(numout,*) '   ==>>>   air humidity is SPECIFIC HUMIDITY     [kg/kg]' 
     337         CASE( np_humi_dpt )   ;   WRITE(numout,*) '   ==>>>   air humidity is DEW-POINT TEMPERATURE [K]' 
     338         CASE( np_humi_rlh )   ;   WRITE(numout,*) '   ==>>>   air humidity is RELATIVE HUMIDITY     [%]' 
    310339         END SELECT 
    311340         ! 
     
    355384      CALL fld_read( kt, nn_fsbc, sf )             ! input fields provided at the current time-step 
    356385      ! 
    357       !                                            ! compute the surface ocean fluxes using bulk formulea 
    358       IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN 
    359          ! 
     386      IF( kt == nit000 ) tsk(:,:) = sst_m(:,:)*tmask(:,:,1)  ! no previous estimate of skin temperature => using bulk SST 
     387      ! 
     388      IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 
    360389         CALL blk_oce_1( kt, sf(jp_wndi)%fnow(:,:,1), sf(jp_wndj)%fnow(:,:,1),   &   !   <<= in 
    361390            &                sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1),   &   !   <<= in 
    362391            &                sf(jp_slp )%fnow(:,:,1), sst_m, ssu_m, ssv_m,       &   !   <<= in 
     392            &                sf(jp_qsr )%fnow(:,:,1), sf(jp_qlw )%fnow(:,:,1),   &   !   <<= in (wl/cs) 
    363393            &                zssq, zcd_du, zsen, zevp )                              !   =>> out 
    364394          
     
    368398            &                zsen, zevp )                                            !   <=> in out 
    369399      ENDIF 
    370           
     400      ! 
    371401#if defined key_cice 
    372       IF( MOD( kt-1, nn_fsbc ) == 0 )   THEN 
     402      IF( MOD( kt - 1, nn_fsbc ) == 0 )   THEN 
    373403         qlw_ice(:,:,1)   = sf(jp_qlw )%fnow(:,:,1) 
    374          IF( ln_dm2dc ) THEN ; qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) 
    375          ELSE                ; qsr_ice(:,:,1) =          sf(jp_qsr)%fnow(:,:,1)  
    376          ENDIF  
     404         IF( ln_dm2dc ) THEN 
     405            qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) 
     406         ELSE 
     407            qsr_ice(:,:,1) =          sf(jp_qsr)%fnow(:,:,1) 
     408         ENDIF 
    377409         tatm_ice(:,:)    = sf(jp_tair)%fnow(:,:,1) 
    378          qatm_ice(:,:)    = sf(jp_humi)%fnow(:,:,1) 
     410 
     411         SELECT CASE( nhumi ) 
     412         CASE( np_humi_sph ) 
     413            qatm_ice(:,:) =           sf(jp_humi)%fnow(:,:,1) 
     414         CASE( np_humi_dpt ) 
     415            qatm_ice(:,:) = q_sat(    sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 
     416         CASE( np_humi_rlh ) 
     417            qatm_ice(:,:) = q_air_rh( 0.01_wp*sf(jp_humi)%fnow(:,:,1), sf(jp_tair)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1)) !LB: 0.01 => RH is % percent in file 
     418         END SELECT 
     419 
    379420         tprecip(:,:)     = sf(jp_prec)%fnow(:,:,1) * rn_pfac 
    380421         sprecip(:,:)     = sf(jp_snow)%fnow(:,:,1) * rn_pfac 
     
    387428 
    388429 
    389    SUBROUTINE blk_oce_1( kt, pwndi, pwndj, ptair, phumi, &  ! inp 
    390               &              pslp,  pst  , pu   , pv,    &  ! inp 
    391               &              pssq, pcd_du, psen, pevp    )  ! out 
     430   SUBROUTINE blk_oce_1( kt, pwndi, pwndj , ptair, phumi, &  ! inp 
     431              &              pslp , pst   , pu   , pv,    &  ! inp 
     432              &              pqsr , pqlw  ,               &  ! inp 
     433              &              pssq , pcd_du, psen , pevp   )  ! out 
    392434      !!--------------------------------------------------------------------- 
    393435      !!                     ***  ROUTINE blk_oce_1  *** 
    394436      !! 
    395       !! ** Purpose :   Computes Cd x |U|, Ch x |U|, Ce x |U| for ABL integration 
    396       !! 
    397       !! ** Method  :   bulk formulae using atmospheric 
    398       !!                fields from the ABL model at previous time-step 
     437      !! ** Purpose :   if ln_blk=T, computes surface momentum, heat and freshwater fluxes 
     438      !!                if ln_abl=T, computes Cd x |U|, Ch x |U|, Ce x |U| for ABL integration 
     439      !! 
     440      !! ** Method  :   bulk formulae using atmospheric fields from : 
     441      !!                if ln_blk=T, atmospheric fields read in sbc_read 
     442      !!                if ln_abl=T, the ABL model at previous time-step 
    399443      !! 
    400444      !! ** Outputs : - pssq    : surface humidity used to compute latent heat flux (kg/kg) 
     
    412456      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pu     ! surface current at U-point (i-component) [m/s] 
    413457      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pv     ! surface current at V-point (j-component) [m/s] 
     458      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pqsr   !  
     459      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pqlw   !  
    414460      REAL(wp), INTENT(  out), DIMENSION(:,:) ::   pssq   ! specific humidity at pst                 [kg/kg] 
    415461      REAL(wp), INTENT(  out), DIMENSION(:,:) ::   pcd_du ! Cd x |dU| at T-points                    [m/s] 
     
    423469      REAL(wp), DIMENSION(jpi,jpj) ::   zU_zu             ! bulk wind speed at height zu  [m/s] 
    424470      REAL(wp), DIMENSION(jpi,jpj) ::   ztpot             ! potential temperature of air at z=rn_zqt [K] 
    425       REAL(wp), DIMENSION(jpi,jpj) ::   zrhoa             ! density of air   [kg/m^3] 
     471      REAL(wp), DIMENSION(jpi,jpj) ::   zqair             ! specific humidity     of air at z=rn_zqt [kg/kg] 
    426472      REAL(wp), DIMENSION(jpi,jpj) ::   zcd_oce           ! momentum transfert coefficient over ocean 
    427473      REAL(wp), DIMENSION(jpi,jpj) ::   zch_oce           ! sensible heat transfert coefficient over ocean 
    428474      REAL(wp), DIMENSION(jpi,jpj) ::   zce_oce           ! latent   heat transfert coefficient over ocean 
     475      REAL(wp), DIMENSION(jpi,jpj) ::   zqla              ! latent heat flux 
     476      REAL(wp), DIMENSION(jpi,jpj) ::   zztmp1, zztmp2 
    429477      !!--------------------------------------------------------------------- 
    430478      ! 
     
    460508 
    461509      ! ----------------------------------------------------------------------------- ! 
    462       !     I    Turbulent FLUXES                                                    ! 
     510      !      I   Solar FLUX                                                           ! 
    463511      ! ----------------------------------------------------------------------------- ! 
    464512 
    465       ! ... specific humidity at SST and IST tmask( 
    466       pssq(:,:) = 0.98 * q_sat( zst(:,:), pslp(:,:) ) 
    467       ! 
     513      ! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle                    ! Short Wave 
     514      zztmp = 1. - albo 
     515      IF( ln_dm2dc ) THEN 
     516         qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1) 
     517      ELSE 
     518         qsr(:,:) = zztmp *          sf(jp_qsr)%fnow(:,:,1)   * tmask(:,:,1) 
     519      ENDIF 
     520 
     521 
     522      ! ----------------------------------------------------------------------------- ! 
     523      !     II   Turbulent FLUXES                                                     ! 
     524      ! ----------------------------------------------------------------------------- ! 
     525 
     526      ! specific humidity at SST 
     527      pssq(:,:) = rdct_qsat_salt * q_sat( zst(:,:), pslp(:,:) ) 
     528 
     529      IF( ln_skin_cs .OR. ln_skin_wl ) THEN 
     530         zztmp1(:,:) = zst(:,:) 
     531         zztmp2(:,:) = pssq(:,:) 
     532      ENDIF 
     533 
     534      ! specific humidity of air at "rn_zqt" m above the sea 
     535      SELECT CASE( nhumi ) 
     536      CASE( np_humi_sph ) 
     537         zqair(:,:) = phumi(:,:)      ! what we read in file is already a spec. humidity! 
     538      CASE( np_humi_dpt ) 
     539         !IF(lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of d_air and slp !' !LBrm 
     540         zqair(:,:) = q_sat( phumi(:,:), pslp(:,:) ) 
     541      CASE( np_humi_rlh ) 
     542         !IF(lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of RH, t_air and slp !' !LBrm 
     543         zqair(:,:) = q_air_rh( 0.01_wp*phumi(:,:), ptair(:,:), pslp(:,:) ) !LB: 0.01 => RH is % percent in file 
     544      END SELECT 
     545      ! 
     546      ! potential temperature of air at "rn_zqt" m above the sea 
    468547      IF( ln_abl ) THEN 
    469548         ztpot = ptair(:,:) 
     
    472551         !    (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2 
    473552         !    (since reanalysis products provide T at z, not theta !) 
    474          ztpot = ptair(:,:) + gamma_moist( ptair(:,:), phumi(:,:) ) * rn_zqt 
    475       ENDIF 
    476  
    477       SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point 
    478       ! 
    479       CASE( np_NCAR      )   ;   CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, pssq, phumi, wndm,   &  ! NCAR-COREv2 
    480          &                                           zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, Cdn_oce, Chn_oce, Cen_oce ) 
    481       CASE( np_COARE_3p0 )   ;   CALL turb_coare   ( rn_zqt, rn_zu, zst, ztpot, pssq, phumi, wndm,   &  ! COARE v3.0 
    482          &                                           zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, Cdn_oce, Chn_oce, Cen_oce ) 
    483       CASE( np_COARE_3p5 )   ;   CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, pssq, phumi, wndm,   &  ! COARE v3.5 
    484          &                                           zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, Cdn_oce, Chn_oce, Cen_oce ) 
    485       CASE( np_ECMWF     )   ;   CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, pssq, phumi, wndm,   &  ! ECMWF 
    486          &                                           zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, Cdn_oce, Chn_oce, Cen_oce ) 
    487       CASE DEFAULT 
    488          CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' ) 
    489       END SELECT 
    490  
    491  
    492       IF ( ln_abl ) THEN         !==  ABL formulation  ==!   multiplication by rho_air and turbulent fluxes computation done in ablstp 
    493 !! FL do we need this multiplication by tmask ... ??? 
    494          DO jj = 1, jpj 
    495             DO ji = 1, jpi 
    496                zztmp = zU_zu(ji,jj) !* tmask(ji,jj,1) 
    497                wndm(ji,jj)   = zztmp                   ! Store zU_zu in wndm to compute ustar2 in ablmod  
    498                pcd_du(ji,jj) = zztmp * zcd_oce(ji,jj) 
    499                psen(ji,jj)   = zztmp * zch_oce(ji,jj) 
    500                pevp(ji,jj)   = zztmp * zce_oce(ji,jj)        
    501             END DO 
    502          END DO 
    503  
    504       ELSE                       !==  BLK formulation  ==!   turbulent fluxes computation 
     553         ztpot = ptair(:,:) + gamma_moist( ptair(:,:), zqair(:,:) ) * rn_zqt 
     554      ENDIF 
     555 
     556 
     557 
     558      IF( ln_skin_cs .OR. ln_skin_wl ) THEN 
     559 
     560         SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point 
     561 
     562         CASE( np_COARE_3p0 ) 
     563            CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 
     564               &                zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,   & 
     565               &                Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) ) 
     566 
     567         CASE( np_COARE_3p6 ) 
     568            CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 
     569               &                zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,   & 
     570               &                Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) ) 
     571 
     572         CASE( np_ECMWF     ) 
     573            CALL turb_ecmwf   ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl,  & 
     574               &                zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,   & 
     575               &                Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) ) 
     576 
     577         CASE DEFAULT 
     578            CALL ctl_stop( 'STOP', 'sbc_oce: unsuported bulk formula selection for "ln_skin_*==.true."' ) 
     579         END SELECT 
     580 
     581         !! In the presence of sea-ice we forget about the cool-skin/warm-layer update of zst and pssq: 
     582         WHERE ( fr_i < 0.001_wp ) 
     583            ! zst and pssq have been updated by cool-skin/warm-layer scheme and we keep it!!! 
     584            zst(:,:)  =  zst(:,:)*tmask(:,:,1) 
     585            pssq(:,:) = pssq(:,:)*tmask(:,:,1) 
     586         ELSEWHERE 
     587            ! we forget about the update... 
     588            zst(:,:)  = zztmp1(:,:) !LB: using what we backed up before skin-algo 
     589            pssq(:,:) = zztmp2(:,:) !LB:  "   "   " 
     590         END WHERE 
     591 
     592         !LB: Update of tsk, the "official" array for skin temperature 
     593         tsk(:,:) = zst(:,:) 
     594 
     595      ELSE !IF( ln_skin_cs .OR. ln_skin_wl ) 
     596 
     597         SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point 
     598            ! 
     599         CASE( np_NCAR      ) 
     600            CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm,                                & 
     601               &                zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
     602 
     603         CASE( np_COARE_3p0 ) 
     604            CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl,   & 
     605               &                 zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
     606 
     607         CASE( np_COARE_3p6 ) 
     608            CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl,   & 
     609               &                 zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
     610 
     611         CASE( np_ECMWF     ) 
     612            CALL turb_ecmwf   ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl,    & 
     613               &                zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
     614 
     615         CASE DEFAULT 
     616            CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' ) 
     617         END SELECT 
     618 
     619      ENDIF ! IF( ln_skin_cs .OR. ln_skin_wl ) 
     620 
     621      !!      CALL iom_put( "Cd_oce", zcd_oce)  ! output value of pure ocean-atm. transfer coef. 
     622      !!      CALL iom_put( "Ch_oce", zch_oce)  ! output value of pure ocean-atm. transfer coef. 
     623       
     624      IF( ABS(rn_zu - rn_zqt) < 0.1_wp ) THEN 
     625         !! If zu == zt, then ensuring once for all that: 
     626         t_zu(:,:) = ztpot(:,:) 
     627         q_zu(:,:) = zqair(:,:) 
     628      ENDIF 
     629       
     630 
     631      !  Turbulent fluxes over ocean  => BULK_FORMULA @ sbcblk_phy.F90 
     632      ! ------------------------------------------------------------- 
     633       
     634      IF (ln_abl) THEN 
     635 
     636         CALL BULK_FORMULA( rn_zu, zst(:,:), pssq(:,:), t_zu(:,:), q_zu(:,:), & 
     637            &               zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:),         & 
     638            &               wndm(:,:), zU_zu(:,:), pslp(:,:),                 & 
     639            &               pcd_du(:,:), psen(:,:), pevp(:,:),                & 
     640            &               prhoa=rhoa(:,:) ) 
     641 
     642      ELSE 
     643 
     644         CALL BULK_FORMULA( rn_zu, zst(:,:), pssq(:,:), t_zu(:,:), q_zu(:,:), & 
     645            &               zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:),         & 
     646            &               wndm(:,:), zU_zu(:,:), pslp(:,:),                 & 
     647            &               taum(:,:), psen(:,:), zqla(:,:),                  & 
     648            &               pEvap=pevp(:,:), prhoa=rhoa(:,:) ) 
     649 
     650         zqla(:,:) = zqla(:,:) * tmask(:,:,1) 
     651         psen(:,:) = psen(:,:) * tmask(:,:,1) 
     652         taum(:,:) = taum(:,:) * tmask(:,:,1) 
     653         pevp(:,:) = pevp(:,:) * tmask(:,:,1) 
    505654          
    506          !                             ! Compute true air density : 
    507          IF( ABS(rn_zu - rn_zqt) > 0.01 ) THEN     ! At zu: (probably useless to remove zrho*grav*rn_zu from SLP...) 
    508             zrhoa(:,:) = rho_air( t_zu(:,:) , q_zu(:,:) , pslp(:,:) ) 
    509          ELSE                                      ! At zt: 
    510             zrhoa(:,:) = rho_air( ptair(:,:), phumi(:,:), pslp(:,:) ) 
    511          END IF 
    512           
    513          DO jj = 1, jpj 
    514             DO ji = 1, jpi 
    515                zztmp         = zrhoa(ji,jj) * zU_zu(ji,jj) 
    516 !!gm to be done by someone: check the efficiency of the call of cp_air within do loops  
    517                psen  (ji,jj) = cp_air( q_zu(ji,jj) ) * zztmp * zch_oce(ji,jj) * (  zst(ji,jj) - t_zu(ji,jj) ) 
    518                pevp  (ji,jj) = rn_efac*MAX( 0._wp,     zztmp * zce_oce(ji,jj) * ( pssq(ji,jj) - q_zu(ji,jj) ) ) 
    519                zztmp         = zztmp * zcd_oce(ji,jj) 
    520                pcd_du(ji,jj) = zztmp 
    521                taum  (ji,jj) = zztmp * wndm  (ji,jj)  
    522                zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) 
    523                zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj) 
    524             END DO 
    525          END DO 
     655         ! Tau i and j component on T-grid points, using array "zcd_oce" as a temporary array... 
     656         zcd_oce = 0._wp 
     657         WHERE ( wndm > 0._wp ) zcd_oce = taum / wndm 
     658         zwnd_i = zcd_oce * zwnd_i 
     659         zwnd_j = zcd_oce * zwnd_j       
    526660 
    527661         CALL iom_put( "taum_oce", taum )   ! output wind stress module 
    528           
     662    
    529663         ! ... utau, vtau at U- and V_points, resp. 
    530664         !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines 
     
    539673         END DO 
    540674         CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. ) 
    541  
    542675         IF(ln_ctl) THEN 
    543676            CALL prt_ctl( tab2d_1=wndm  , clinfo1=' blk_oce_1: wndm   : ') 
     
    554687      ENDIF 
    555688      ! 
    556   END SUBROUTINE blk_oce_1 
     689   END SUBROUTINE blk_oce_1 
    557690 
    558691 
     
    593726      zst(:,:) = pst(:,:) + rt0      ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K) 
    594727 
     728 
    595729      ! ----------------------------------------------------------------------------- ! 
    596       !      I   Radiative FLUXES                                                     ! 
     730      !     III    Net longwave radiative FLUX                                        ! 
    597731      ! ----------------------------------------------------------------------------- ! 
    598732 
    599       ! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle                    ! Short Wave 
    600       zztmp = 1._wp - albo 
    601       IF( ln_dm2dc ) THEN   ;   qsr(:,:) = zztmp * sbc_dcy( pqsr(:,:) ) * tmask(:,:,1) 
    602       ELSE                  ;   qsr(:,:) = zztmp *          pqsr(:,:)   * tmask(:,:,1) 
    603       ENDIF 
    604  
    605       zqlw(:,:) = ( pqlw(:,:) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:)  ) * tmask(:,:,1)   ! Long  Wave 
     733      !! LB: now moved after Turbulent fluxes because must use the skin temperature rather that the SST 
     734      !! (zst is skin temperature if ln_skin_cs==.TRUE. .OR. ln_skin_wl==.TRUE.) 
     735      zqlw(:,:) = emiss_w * ( pqlw(:,:) - stefan*zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:) ) * tmask(:,:,1)   ! Net radiative longwave flux 
    606736 
    607737      !  Turbulent fluxes over ocean 
    608738      ! ----------------------------- 
    609739 
    610       zqla(:,:) = L_vap( zst(:,:) ) * pevp(:,:)     ! Latent Heat flux 
     740      zqla(:,:) = L_vap( zst(:,:) ) * pevp(:,:)     ! Latent Heat flux !!GS: possibility to add a global qla to avoid recomputation after abl update 
    611741 
    612742      IF(ln_ctl) THEN 
    613743         CALL prt_ctl( tab2d_1=zqla  , clinfo1=' blk_oce_2: zqla   : ' ) 
    614744         CALL prt_ctl( tab2d_1=zqlw  , clinfo1=' blk_oce_2: zqlw   : ', tab2d_2=qsr, clinfo2=' qsr : ' ) 
     745 
    615746      ENDIF 
    616747 
    617748      ! ----------------------------------------------------------------------------- ! 
    618       !     III    Total FLUXES                                                       ! 
     749      !     IV    Total FLUXES                                                       ! 
    619750      ! ----------------------------------------------------------------------------- ! 
    620751      ! 
    621       emp (:,:) = ( pevp(:,:) - pprec(:,:) * rn_pfac  ) * tmask(:,:,1)             ! mass flux (evap. - precip.) 
    622       ! 
    623       zz1 = rn_pfac * rLfus 
    624       zz2 = rn_pfac * rcp 
    625       zz3 = rn_pfac * rcpi 
    626       zztmp = 1._wp - albo 
    627       qns(:,:) = (  zqlw(:,:) - psen(:,:) - zqla(:,:)                          &   ! Downward Non Solar 
    628          &        - psnow(:,:) * zztmp                                         &   ! remove latent melting heat for solid precip 
    629          &        - pevp(:,:) * pst(:,:) * rcp                                 &   ! remove evap heat content at SST 
    630          &        + ( pprec(:,:) - psnow(:,:) ) * ( ptair(:,:) - rt0 ) * zz2   &   ! add liquid precip heat content at Tair 
    631          &        + psnow(:,:) * ( MIN( ptair(:,:), rt0 ) - rt0 ) * zz3        &   ! add solid  precip heat content at min(Tair,Tsnow) 
    632          &       ) * tmask(:,:,1) 
     752      emp (:,:) = (  pevp(:,:)                                       &   ! mass flux (evap. - precip.) 
     753         &         - pprec(:,:) * rn_pfac  ) * tmask(:,:,1) 
     754      ! 
     755      qns(:,:) = zqlw(:,:) + psen(:,:) + zqla(:,:)                   &   ! Downward Non Solar 
     756         &     - psnow(:,:) * rn_pfac * rLfus                        &   ! remove latent melting heat for solid precip 
     757         &     - pevp(:,:) * pst(:,:) * rcp                          &   ! remove evap heat content at SST !LB??? pst is Celsius !? 
     758         &     + ( pprec(:,:) - psnow(:,:) ) * rn_pfac               &   ! add liquid precip heat content at Tair 
     759         &     * ( ptair(:,:) - rt0 ) * rcp                          & 
     760         &     + psnow(:,:) * rn_pfac                                &   ! add solid  precip heat content at min(Tair,Tsnow) 
     761         &     * ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi 
     762      qns(:,:) = qns(:,:) * tmask(:,:,1) 
    633763      ! 
    634764#if defined key_si3 
    635       qns_oce(:,:) = zqlw(:,:) - psen(:,:) - zqla(:,:)                             ! non solar without emp (only needed by SI3) 
     765      qns_oce(:,:) = zqlw(:,:) + psen(:,:) + zqla(:,:)                             ! non solar without emp (only needed by SI3) 
    636766      qsr_oce(:,:) = qsr(:,:) 
    637767#endif 
    638768      ! 
    639       IF ( nn_ice == 0 ) THEN 
    640          CALL iom_put( "qlw_oce" ,   zqlw )                 ! output downward longwave heat over the ocean 
    641          CALL iom_put( "qsb_oce" , - psen )                 ! output downward sensible heat over the ocean 
    642          CALL iom_put( "qla_oce" , - zqla )                 ! output downward latent   heat over the ocean 
    643          CALL iom_put( "qemp_oce",   qns-zqlw+psen+zqla )   ! output downward heat content of E-P over the ocean 
    644          CALL iom_put( "qns_oce" ,   qns  )                 ! output downward non solar heat over the ocean 
    645          CALL iom_put( "qsr_oce" ,   qsr  )                 ! output downward solar heat over the ocean 
    646          CALL iom_put( "qt_oce"  ,   qns+qsr )              ! output total downward heat over the ocean 
    647          tprecip(:,:) = pprec(:,:) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s] 
    648          sprecip(:,:) = psnow(:,:) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s] 
    649          CALL iom_put( 'snowpre', sprecip )                 ! Snow 
    650          CALL iom_put( 'precip' , tprecip )                 ! Total precipitation 
     769      CALL iom_put( "rho_air"  ,   rhoa )                 ! output air density (kg/m^3) !#LB 
     770      CALL iom_put( "qlw_oce"  ,   zqlw )                 ! output downward longwave heat over the ocean 
     771      CALL iom_put( "qsb_oce"  ,   psen )                 ! output downward sensible heat over the ocean 
     772      CALL iom_put( "qla_oce"  ,   zqla )                 ! output downward latent   heat over the ocean 
     773      CALL iom_put( "evap_oce" ,   pevp )                 ! evaporation 
     774      CALL iom_put( "qemp_oce" ,   qns-zqlw-psen-zqla )   ! output downward heat content of E-P over the ocean 
     775      CALL iom_put( "qns_oce"  ,   qns  )                 ! output downward non solar heat over the ocean 
     776      CALL iom_put( "qsr_oce"  ,   qsr  )                 ! output downward solar heat over the ocean 
     777      CALL iom_put( "qt_oce"   ,   qns+qsr )              ! output total downward heat over the ocean 
     778      tprecip(:,:) = pprec(:,:) * rn_pfac * tmask(:,:,1)  ! output total precipitation [kg/m2/s] 
     779      sprecip(:,:) = psnow(:,:) * rn_pfac * tmask(:,:,1)  ! output solid precipitation [kg/m2/s] 
     780      CALL iom_put( 'snowpre', sprecip )                  ! Snow 
     781      CALL iom_put( 'precip' , tprecip )                  ! Total precipitation 
     782      IF( ln_skin_cs .OR. ln_skin_wl ) THEN 
     783         CALL iom_put( "t_skin" ,  (zst - rt0) * tmask(:,:,1) )           ! T_skin in Celsius 
     784         CALL iom_put( "dt_skin" , (zst - pst - rt0) * tmask(:,:,1) )     ! T_skin - SST temperature difference... 
    651785      ENDIF 
    652786      ! 
     
    660794 
    661795    
    662    FUNCTION rho_air( ptak, pqa, pslp ) 
    663       !!------------------------------------------------------------------------------- 
    664       !!                           ***  FUNCTION rho_air  *** 
    665       !! 
    666       !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere 
    667       !! 
    668       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)  
    669       !!------------------------------------------------------------------------------- 
    670       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak      ! air temperature             [K] 
    671       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa       ! air specific humidity   [kg/kg] 
    672       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp      ! pressure in                [Pa] 
    673       REAL(wp), DIMENSION(jpi,jpj)             ::   rho_air   ! density of moist air   [kg/m^3] 
    674       !!------------------------------------------------------------------------------- 
    675       ! 
    676       rho_air = pslp / (  R_dry*ptak * ( 1._wp + rctv0*pqa )  ) 
    677       ! 
    678    END FUNCTION rho_air 
    679  
    680  
    681    FUNCTION cp_air_0d( pqa ) 
    682       !!------------------------------------------------------------------------------- 
    683       !!                           ***  FUNCTION cp_air  *** 
    684       !! 
    685       !! ** Purpose : Compute specific heat (Cp) of moist air 
    686       !! 
    687       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    688       !!------------------------------------------------------------------------------- 
    689       REAL(wp), INTENT(in) ::   pqa      ! air specific humidity         [kg/kg] 
    690       REAL(wp)             ::   cp_air_0d! specific heat of moist air   [J/K/kg] 
    691       !!------------------------------------------------------------------------------- 
    692       ! 
    693       Cp_air_0d = Cp_dry + Cp_vap * pqa 
    694       ! 
    695    END FUNCTION cp_air_0d 
    696  
    697  
    698    FUNCTION cp_air_2d( pqa ) 
    699       !!------------------------------------------------------------------------------- 
    700       !!                           ***  FUNCTION cp_air  *** 
    701       !! 
    702       !! ** Purpose : Compute specific heat (Cp) of moist air 
    703       !! 
    704       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    705       !!------------------------------------------------------------------------------- 
    706       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa      ! air specific humidity         [kg/kg] 
    707       REAL(wp), DIMENSION(jpi,jpj)             ::   cp_air_2d! specific heat of moist air   [J/K/kg] 
    708       !!------------------------------------------------------------------------------- 
    709       ! 
    710       Cp_air_2d = Cp_dry + Cp_vap * pqa 
    711       ! 
    712    END FUNCTION cp_air_2d 
    713  
    714  
    715    FUNCTION q_sat( ptak, pslp ) 
    716       !!---------------------------------------------------------------------------------- 
    717       !!                           ***  FUNCTION q_sat  *** 
    718       !! 
    719       !! ** Purpose : Specific humidity at saturation in [kg/kg] 
    720       !!              Based on accurate estimate of "e_sat" 
    721       !!              aka saturation water vapor (Goff, 1957) 
    722       !! 
    723       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    724       !!---------------------------------------------------------------------------------- 
    725       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak    ! air temperature                       [K] 
    726       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp    ! sea level atmospheric pressure       [Pa] 
    727       REAL(wp), DIMENSION(jpi,jpj)             ::   q_sat   ! Specific humidity at saturation   [kg/kg] 
    728       ! 
    729       INTEGER  ::   ji, jj         ! dummy loop indices 
    730       REAL(wp) ::   ze_sat, ztmp   ! local scalar 
    731       !!---------------------------------------------------------------------------------- 
    732       ! 
    733       DO jj = 1, jpj 
    734          DO ji = 1, jpi 
    735             ! 
    736             ztmp = rt0 / ptak(ji,jj) 
    737             ! 
    738             ! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957) 
    739             ze_sat = 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(ptak(ji,jj)/rt0)        & 
    740                &    + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak(ji,jj)/rt0 - 1.)) )  & 
    741                &    + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614  ) 
    742                ! 
    743             q_sat(ji,jj) = reps0 * ze_sat/( 0.01_wp*pslp(ji,jj) - (1._wp - reps0)*ze_sat )   ! 0.01 because SLP is in [Pa] 
    744             ! 
    745          END DO 
    746       END DO 
    747       ! 
    748    END FUNCTION q_sat 
    749  
    750  
    751    FUNCTION gamma_moist( ptak, pqa ) 
    752       !!---------------------------------------------------------------------------------- 
    753       !!                           ***  FUNCTION gamma_moist  *** 
    754       !! 
    755       !! ** Purpose : Compute the moist adiabatic lapse-rate. 
    756       !!     => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate 
    757       !!     => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html 
    758       !! 
    759       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    760       !!---------------------------------------------------------------------------------- 
    761       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak          ! air temperature       [K] 
    762       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa           ! specific humidity [kg/kg] 
    763       REAL(wp), DIMENSION(jpi,jpj)             ::   gamma_moist   ! moist adiabatic lapse-rate 
    764       ! 
    765       INTEGER  ::   ji, jj         ! dummy loop indices 
    766       REAL(wp) :: zrv, ziRT        ! local scalar 
    767       !!---------------------------------------------------------------------------------- 
    768       ! 
    769       DO jj = 1, jpj 
    770          DO ji = 1, jpi 
    771             zrv = pqa(ji,jj) / (1. - pqa(ji,jj)) 
    772             ziRT = 1. / (R_dry*ptak(ji,jj))    ! 1/RT 
    773             gamma_moist(ji,jj) = grav * ( 1. + rLevap*zrv*ziRT ) / ( Cp_dry + rLevap*rLevap*zrv*reps0*ziRT/ptak(ji,jj) ) 
    774          END DO 
    775       END DO 
    776       ! 
    777    END FUNCTION gamma_moist 
    778  
    779  
    780    FUNCTION L_vap( psst ) 
    781       !!--------------------------------------------------------------------------------- 
    782       !!                           ***  FUNCTION L_vap  *** 
    783       !! 
    784       !! ** Purpose : Compute the latent heat of vaporization of water from temperature 
    785       !! 
    786       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    787       !!---------------------------------------------------------------------------------- 
    788       REAL(wp), DIMENSION(jpi,jpj)             ::   L_vap   ! latent heat of vaporization   [J/kg] 
    789       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   psst   ! water temperature                [K] 
    790       !!---------------------------------------------------------------------------------- 
    791       ! 
    792       L_vap = (  2.501 - 0.00237 * ( psst(:,:) - rt0)  ) * 1.e6 
    793       ! 
    794    END FUNCTION L_vap 
    795  
    796796#if defined key_si3 
    797797   !!---------------------------------------------------------------------- 
     
    802802   !!   blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux) 
    803803   !!   Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag 
    804    !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag  
     804   !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag 
    805805   !!---------------------------------------------------------------------- 
    806806 
     
    835835      REAL(wp) ::   zootm_su                      ! sea-ice surface mean temperature 
    836836      REAL(wp) ::   zztmp1, zztmp2                ! temporary arrays 
    837       REAL(wp), DIMENSION(jpi,jpj) ::   zrhoa     ! transfer coefficient for momentum      (tau) 
    838837      REAL(wp), DIMENSION(jpi,jpj) ::   zcd_dui   ! transfer coefficient for momentum      (tau) 
    839838      !!--------------------------------------------------------------------- 
     839      ! 
    840840 
    841841      ! ------------------------------------------------------------ ! 
     
    858858         Ce_ice(:,:) = Cd_ice(:,:) 
    859859      ELSEIF( ln_Cd_L15 ) THEN   ! calculate new drag from Lupkes(2015) equations 
    860          CALL Cdn10_Lupkes2015( ptsui, Cd_ice, Ch_ice )  
     860         CALL Cdn10_Lupkes2015( ptsui, pslp, Cd_ice, Ch_ice )  
    861861         Ce_ice(:,:) = Ch_ice(:,:)       ! sensible and latent heat transfer coef. are considered identical 
    862862      ENDIF 
    863863 
    864       IF ( iom_use("Cd_ice") ) CALL iom_put("Cd_ice", Cd_ice)   ! output value of pure ice-atm. transfer coef. 
    865       IF ( iom_use("Ch_ice") ) CALL iom_put("Ch_ice", Ch_ice)   ! output value of pure ice-atm. transfer coef. 
     864      !! IF ( iom_use("Cd_ice") ) CALL iom_put("Cd_ice", Cd_ice)   ! output value of pure ice-atm. transfer coef. 
     865      !! IF ( iom_use("Ch_ice") ) CALL iom_put("Ch_ice", Ch_ice)   ! output value of pure ice-atm. transfer coef. 
    866866 
    867867      ! local scalars ( place there for vector optimisation purposes) 
    868       ! Computing density of air! Way denser that 1.2 over sea-ice !!! 
    869       zrhoa  (:,:) = rho_air( ptair(:,:), phumi(:,:), pslp(:,:) ) 
     868      IF (ln_abl) rhoa  (:,:)  = rho_air( ptair(:,:), phumi(:,:), pslp(:,:) ) !!GS: rhoa must be (re)computed here with ABL to avoid division by zero after (TBI) 
    870869      zcd_dui(:,:) = wndm_ice(:,:) * Cd_ice(:,:) 
    871870       
     
    877876         DO jj = 2, jpjm1 
    878877            DO ji = fs_2, fs_jpim1   ! vect. opt. 
    879                putaui(ji,jj) = 0.5_wp * (  zrhoa(ji+1,jj) * zcd_dui(ji+1,jj)             & 
    880                   &                      + zrhoa(ji  ,jj) * zcd_dui(ji  ,jj)  )          & 
     878               putaui(ji,jj) = 0.5_wp * (  rhoa(ji+1,jj) * zcd_dui(ji+1,jj)             & 
     879                  &                      + rhoa(ji  ,jj) * zcd_dui(ji  ,jj)  )          & 
    881880                  &         * ( 0.5_wp * ( pwndi(ji+1,jj) + pwndi(ji,jj) ) - rn_vfac * puice(ji,jj) ) 
    882                pvtaui(ji,jj) = 0.5_wp * (  zrhoa(ji,jj+1) * zcd_dui(ji,jj+1)             & 
    883                   &                      + zrhoa(ji,jj  ) * zcd_dui(ji,jj  )  )          & 
     881               pvtaui(ji,jj) = 0.5_wp * (  rhoa(ji,jj+1) * zcd_dui(ji,jj+1)             & 
     882                  &                      + rhoa(ji,jj  ) * zcd_dui(ji,jj  )  )          & 
    884883                  &         * ( 0.5_wp * ( pwndj(ji,jj+1) + pwndj(ji,jj) ) - rn_vfac * pvice(ji,jj) ) 
    885884            END DO 
     
    898897               pevpi  (ji,jj) = wndm_ice(ji,jj) * Ce_ice(ji,jj) 
    899898               zootm_su       = zztmp2 / ptsui(ji,jj)   ! ptsui is in K (it can't be zero ??) 
    900                pssqi  (ji,jj) = zztmp1 * EXP( zootm_su ) / zrhoa(ji,jj) 
     899               pssqi  (ji,jj) = zztmp1 * EXP( zootm_su ) / rhoa(ji,jj) 
    901900            END DO 
    902901         END DO 
     
    934933      REAL(wp) ::   zst3                     ! local variable 
    935934      REAL(wp) ::   zcoef_dqlw, zcoef_dqla   !   -      - 
    936       REAL(wp) ::   zztmp, zztmp2, z1_rLsub          !   -      - 
     935      REAL(wp) ::   zztmp, zztmp2, z1_rLsub  !   -      - 
    937936      REAL(wp) ::   zfr1, zfr2               ! local variables 
    938937      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_st         ! inverse of surface temperature 
     
    942941      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_dqsb        ! sensible  heat sensitivity over ice 
    943942      REAL(wp), DIMENSION(jpi,jpj)     ::   zevap, zsnw   ! evaporation and snw distribution after wind blowing (SI3) 
    944       REAL(wp), DIMENSION(jpi,jpj)     ::   zrhoa 
    945       !!--------------------------------------------------------------------- 
    946       ! 
    947       zcoef_dqlw = 4.0 * 0.95 * Stef             ! local scalars 
    948       zcoef_dqla = -Ls * 11637800. * (-5897.8) 
    949       ! 
    950       zrhoa(:,:) = rho_air( ptair(:,:), phumi(:,:), pslp(:,:) ) 
     943      REAL(wp), DIMENSION(jpi,jpj)     ::   zqair         ! specific humidity of air at z=rn_zqt [kg/kg] !LB 
     944      !!--------------------------------------------------------------------- 
     945      ! 
     946      zcoef_dqlw = 4._wp * 0.95_wp * stefan             ! local scalars 
     947      zcoef_dqla = -rLsub * 11637800._wp * (-5897.8_wp) !LB: BAD! 
     948      ! 
     949      SELECT CASE( nhumi ) 
     950      CASE( np_humi_sph ) 
     951         zqair(:,:) =  phumi(:,:)      ! what we read in file is already a spec. humidity! 
     952      CASE( np_humi_dpt ) 
     953         zqair(:,:) = q_sat( phumi(:,:), pslp ) 
     954      CASE( np_humi_rlh ) 
     955         zqair(:,:) = q_air_rh( 0.01_wp*phumi(:,:), ptair(:,:), pslp(:,:) ) !LB: 0.01 => RH is % percent in file 
     956      END SELECT 
    951957      ! 
    952958      zztmp = 1. / ( 1. - albo ) 
    953       WHERE( ptsu(:,:,:) /= 0._wp )   ;   z1_st(:,:,:) = 1._wp / ptsu(:,:,:) 
    954       ELSEWHERE                       ;   z1_st(:,:,:) = 0._wp 
     959      WHERE( ptsu(:,:,:) /= 0._wp ) 
     960         z1_st(:,:,:) = 1._wp / ptsu(:,:,:) 
     961      ELSEWHERE 
     962         z1_st(:,:,:) = 0._wp 
    955963      END WHERE 
    956964      !                                     ! ========================== ! 
     
    966974               qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj) 
    967975               ! Long  Wave (lw) 
    968                z_qlw(ji,jj,jl) = 0.95 * ( pqlw(ji,jj) - Stef * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) 
     976               z_qlw(ji,jj,jl) = 0.95 * ( pqlw(ji,jj) - stefan * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) 
    969977               ! lw sensitivity 
    970978               z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3 
     
    976984               ! ... turbulent heat fluxes with Ch_ice recalculated in blk_ice_1 
    977985               ! Sensible Heat 
    978                z_qsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_ice(ji,jj) * wndm_ice(ji,jj) * (ptsu(ji,jj,jl) - ptair(ji,jj)) 
     986               z_qsb(ji,jj,jl) = rhoa(ji,jj) * rCp_air * Ch_ice(ji,jj) * wndm_ice(ji,jj) * (ptsu(ji,jj,jl) - ptair(ji,jj)) 
    979987               ! Latent Heat 
    980988               zztmp2 = EXP( -5897.8 * z1_st(ji,jj,jl) ) 
    981                qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls  * Ce_ice(ji,jj) * wndm_ice(ji,jj) *  & 
    982                   &                ( 11637800. * zztmp2 / zrhoa(ji,jj) - phumi(ji,jj) ) ) 
     989               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa(ji,jj) * rLsub  * Ce_ice(ji,jj) * wndm_ice(ji,jj) *  & 
     990                  &                ( 11637800. * zztmp2 / rhoa(ji,jj) - zqair(ji,jj) ) ) 
    983991               ! Latent heat sensitivity for ice (Dqla/Dt) 
    984992               IF( qla_ice(ji,jj,jl) > 0._wp ) THEN 
     
    990998 
    991999               ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice) 
    992                z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_ice(ji,jj) * wndm_ice(ji,jj) 
    993                 
     1000               z_dqsb(ji,jj,jl) = rhoa(ji,jj) * rCp_air * Ch_ice(ji,jj) * wndm_ice(ji,jj) 
     1001 
    9941002               ! ----------------------------! 
    9951003               !     III    Total FLUXES     ! 
     
    10421050      DO jl = 1, jpl 
    10431051         qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * rcpi * tmask(:,:,1) ) 
    1044          !                         ! But we do not have Tice => consider it at 0degC => evap=0  
     1052         !                         ! But we do not have Tice => consider it at 0degC => evap=0 
    10451053      END DO 
    10461054 
     
    10491057      zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )            ! zfr2 such that zfr1 + zfr2 to equal 1 
    10501058      ! 
    1051       WHERE    ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) <  0.1_wp )       ! linear decrease from hi=0 to 10cm   
     1059      WHERE    ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) <  0.1_wp )       ! linear decrease from hi=0 to 10cm 
    10521060         qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) ) 
    10531061      ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp )       ! constant (zfr1) when hi>10cm 
    10541062         qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * zfr1 
    10551063      ELSEWHERE                                                         ! zero when hs>0 
    1056          qtr_ice_top(:,:,:) = 0._wp  
     1064         qtr_ice_top(:,:,:) = 0._wp 
    10571065      END WHERE 
    10581066      ! 
     
    10641072         CALL prt_ctl(tab3d_1=ptsu    , clinfo1=' blk_ice: ptsu     : ', tab3d_2=qns_ice , clinfo2=' qns_ice  : ', kdim=jpl) 
    10651073         CALL prt_ctl(tab2d_1=tprecip , clinfo1=' blk_ice: tprecip  : ', tab2d_2=sprecip , clinfo2=' sprecip  : ') 
    1066       END IF 
     1074      ENDIF 
    10671075      ! 
    10681076   END SUBROUTINE blk_ice_2 
    1069     
     1077 
    10701078 
    10711079   SUBROUTINE blk_ice_qcn( ld_virtual_itd, ptsu, ptb, phs, phi ) 
     
    10761084      !!                to force sea ice / snow thermodynamics 
    10771085      !!                in the case conduction flux is emulated 
    1078       !!                 
     1086      !! 
    10791087      !! ** Method  :   compute surface energy balance assuming neglecting heat storage 
    10801088      !!                following the 0-layer Semtner (1976) approach 
     
    11011109      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zgfac   ! enhanced conduction factor 
    11021110      !!--------------------------------------------------------------------- 
    1103        
     1111 
    11041112      ! -------------------------------------! 
    11051113      !      I   Enhanced conduction factor  ! 
     
    11091117      ! 
    11101118      zgfac(:,:,:) = 1._wp 
    1111        
     1119 
    11121120      IF( ld_virtual_itd ) THEN 
    11131121         ! 
     
    11151123         zfac2 = EXP(1._wp) * 0.5_wp * zepsilon 
    11161124         zfac3 = 2._wp / zepsilon 
    1117          !    
    1118          DO jl = 1, jpl                 
     1125         ! 
     1126         DO jl = 1, jpl 
    11191127            DO jj = 1 , jpj 
    11201128               DO ji = 1, jpi 
     
    11241132            END DO 
    11251133         END DO 
    1126          !       
    1127       ENDIF 
    1128        
     1134         ! 
     1135      ENDIF 
     1136 
    11291137      ! -------------------------------------------------------------! 
    11301138      !      II   Surface temperature and conduction flux            ! 
     
    11361144         DO jj = 1 , jpj 
    11371145            DO ji = 1, jpi 
    1138                !                     
     1146               ! 
    11391147               zkeff_h = zfac * zgfac(ji,jj,jl) / &                                    ! Effective conductivity of the snow-ice system divided by thickness 
    11401148                  &      ( rcnd_i * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) ) 
     
    11531161               qns_ice(ji,jj,jl) = qns_ice(ji,jj,jl) + dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) 
    11541162               qml_ice(ji,jj,jl) = ( qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) - qcn_ice(ji,jj,jl) )  & 
    1155                              &   * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) ) 
     1163                  &   * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) ) 
    11561164 
    11571165               ! --- Diagnose the heat loss due to changing non-solar flux (as in icethd_zdf_bl99) --- ! 
    1158                hfx_err_dif(ji,jj) = hfx_err_dif(ji,jj) - ( dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) ) * a_i_b(ji,jj,jl)  
     1166               hfx_err_dif(ji,jj) = hfx_err_dif(ji,jj) - ( dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) ) * a_i_b(ji,jj,jl) 
    11591167 
    11601168            END DO 
    11611169         END DO 
    11621170         ! 
    1163       END DO  
    1164       !       
     1171      END DO 
     1172      ! 
    11651173   END SUBROUTINE blk_ice_qcn 
    1166     
     1174 
    11671175 
    11681176   SUBROUTINE Cdn10_Lupkes2012( pcd ) 
     
    11701178      !!                      ***  ROUTINE  Cdn10_Lupkes2012  *** 
    11711179      !! 
    1172       !! ** Purpose :    Recompute the neutral air-ice drag referenced at 10m  
     1180      !! ** Purpose :    Recompute the neutral air-ice drag referenced at 10m 
    11731181      !!                 to make it dependent on edges at leads, melt ponds and flows. 
    11741182      !!                 After some approximations, this can be resumed to a dependency 
    11751183      !!                 on ice concentration. 
    1176       !!                 
     1184      !! 
    11771185      !! ** Method :     The parameterization is taken from Lupkes et al. (2012) eq.(50) 
    11781186      !!                 with the highest level of approximation: level4, eq.(59) 
     
    11861194      !! 
    11871195      !!                 This new drag has a parabolic shape (as a function of A) starting at 
    1188       !!                 Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5  
     1196      !!                 Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5 
    11891197      !!                 and going down to Cdi(say 1.4e-3) for A=1 
    11901198      !! 
     
    12111219 
    12121220      ! ice-atm drag 
    1213       pcd(:,:) = rcdice +  &                                                         ! pure ice drag 
     1221      pcd(:,:) = rCd_ice +  &                                                         ! pure ice drag 
    12141222          &      zCe     * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp)  ! change due to sea-ice morphology 
    1215        
     1223 
    12161224   END SUBROUTINE Cdn10_Lupkes2012 
    12171225 
    12181226 
    1219    SUBROUTINE Cdn10_Lupkes2015( ptm_su, pcd, pch ) 
     1227   SUBROUTINE Cdn10_Lupkes2015( ptm_su, pslp, pcd, pch ) 
    12201228      !!---------------------------------------------------------------------- 
    12211229      !!                      ***  ROUTINE  Cdn10_Lupkes2015  *** 
    12221230      !! 
    12231231      !! ** pUrpose :    Alternative turbulent transfert coefficients formulation 
    1224       !!                 between sea-ice and atmosphere with distinct momentum  
    1225       !!                 and heat coefficients depending on sea-ice concentration  
     1232      !!                 between sea-ice and atmosphere with distinct momentum 
     1233      !!                 and heat coefficients depending on sea-ice concentration 
    12261234      !!                 and atmospheric stability (no meltponds effect for now). 
    1227       !!                 
     1235      !! 
    12281236      !! ** Method :     The parameterization is adapted from Lupkes et al. (2015) 
    12291237      !!                 and ECHAM6 atmospheric model. Compared to Lupkes2012 scheme, 
    12301238      !!                 it considers specific skin and form drags (Andreas et al. 2010) 
    1231       !!                 to compute neutral transfert coefficients for both heat and  
     1239      !!                 to compute neutral transfert coefficients for both heat and 
    12321240      !!                 momemtum fluxes. Atmospheric stability effect on transfert 
    12331241      !!                 coefficient is also taken into account following Louis (1979). 
     
    12391247      ! 
    12401248      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   ptm_su ! sea-ice surface temperature [K] 
     1249      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pslp   ! sea-level pressure [Pa] 
    12411250      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pcd    ! momentum transfert coefficient 
    12421251      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pch    ! heat transfert coefficient 
     
    12671276      REAL(wp) ::   z0w, z0i, zfmi, zfmw, zfhi, zfhw 
    12681277      REAL(wp) ::   zCdn_form_tmp 
    1269       !!--------------------------------------------------------------------------- 
    1270       ! 
     1278      !!---------------------------------------------------------------------- 
     1279 
    12711280      ! Momentum Neutral Transfert Coefficients (should be a constant) 
    12721281      zCdn_form_tmp = zce10 * ( LOG( 10._wp / z0_form_ice + 1._wp ) / LOG( rn_zu / z0_form_ice + 1._wp ) )**2   ! Eq. 40 
     
    12771286      ! Heat Neutral Transfert Coefficients 
    12781287      zChn_skin_ice = vkarmn**2 / ( LOG( rn_zu / z0_ice + 1._wp ) * LOG( rn_zu * z1_alpha / z0_skin_ice + 1._wp ) )   ! Eq. 50 + Eq. 52 
    1279       
     1288 
    12801289      ! Atmospheric and Surface Variables 
    12811290      zst(:,:)     = sst_m(:,:) + rt0                                        ! convert SST from Celcius to Kelvin 
    1282       zqo_sat(:,:) = 0.98_wp * q_sat( zst(:,:)   , sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ocean [kg/kg] 
    1283       zqi_sat(:,:) = 0.98_wp * q_sat( ptm_su(:,:), sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ice   [kg/kg] 
    1284       ! 
    1285       DO jj = 2, jpjm1           ! reduced loop is necessary for reproductibility 
     1291      zqo_sat(:,:) = rdct_qsat_salt * q_sat( zst(:,:)   , pslp(:,:) )   ! saturation humidity over ocean [kg/kg] 
     1292      zqi_sat(:,:) =                  q_sat( ptm_su(:,:), pslp(:,:) )   ! saturation humidity over ice   [kg/kg] 
     1293      ! 
     1294      DO jj = 2, jpjm1           ! reduced loop is necessary for reproducibility 
    12861295         DO ji = fs_2, fs_jpim1 
    12871296            ! Virtual potential temperature [K] 
     
    12891298            zthetav_is = ptm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) )   ! ocean ice 
    12901299            zthetav_zu = t_zu (ji,jj)  * ( 1._wp + rctv0 * q_zu(ji,jj)    )   ! at zu 
    1291              
     1300 
    12921301            ! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead) 
    12931302            zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj)     )**2   ! over ocean 
    12941303            zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2   ! over ice 
    1295              
     1304 
    12961305            ! Momentum and Heat Neutral Transfert Coefficients 
    12971306            zCdn_form_ice = zCdn_form_tmp * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta  ! Eq. 40 
    1298             zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) )               ! Eq. 53  
    1299                         
    1300             ! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead) 
     1307            zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) )               ! Eq. 53 
     1308 
     1309            ! Momentum and Heat Stability functions (!!GS: possibility to use psi_m_ecmwf instead ?) 
    13011310            z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water 
    13021311            z0i = z0_skin_ice                                             ! over ice 
     
    13091318               zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) )   ! Eq. 28 
    13101319            ENDIF 
    1311              
     1320 
    13121321            IF( zrib_i <= 0._wp ) THEN 
    13131322               zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp)))   ! Eq.  9 
     
    13171326               zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) )   ! Eq. 27 
    13181327            ENDIF 
    1319              
     1328 
    13201329            ! Momentum Transfert Coefficients (Eq. 38) 
    13211330            pcd(ji,jj) = zCdn_skin_ice *   zfmi +  & 
    13221331                &        zCdn_form_ice * ( zfmi * at_i_b(ji,jj) + zfmw * ( 1._wp - at_i_b(ji,jj) ) ) / MAX( 1.e-06, at_i_b(ji,jj) ) 
    1323              
     1332 
    13241333            ! Heat Transfert Coefficients (Eq. 49) 
    13251334            pch(ji,jj) = zChn_skin_ice *   zfhi +  & 
  • NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/SBC/sbcblk_algo_ecmwf.F90

    r10069 r12015  
    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   !!====================================================================== 
  • NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/SBC/sbcblk_algo_ncar.F90

    r10190 r12015  
    1111   !! 
    1212   !!       Routine turb_ncar maintained and developed in AeroBulk 
    13    !!                     (http://aerobulk.sourceforge.net/) 
     13   !!                     (https://github.com/brodeau/aerobulk/) 
    1414   !! 
    1515   !!                         L. Brodeau, 2015 
     
    2626   USE dom_oce         ! ocean space and time domain 
    2727   USE phycst          ! physical constants 
    28    USE sbc_oce         ! Surface boundary condition: ocean fields 
     28   USE iom             ! I/O manager library 
     29   USE lib_mpp         ! distribued memory computing library 
     30   USE in_out_manager  ! I/O manager 
     31   USE prtctl          ! Print control 
    2932   USE sbcwave, ONLY   :  cdn_wave ! wave module 
    3033#if defined key_si3 || defined key_cice 
    3134   USE sbc_ice         ! Surface boundary condition: ice fields 
    3235#endif 
    33    ! 
    34    USE iom             ! I/O manager library 
    35    USE lib_mpp         ! distribued memory computing library 
    36    USE in_out_manager  ! I/O manager 
    37    USE prtctl          ! Print control 
    3836   USE lib_fortran     ! to use key_nosignedzero 
    3937 
     38   USE sbc_oce         ! Surface boundary condition: ocean fields 
     39   USE sbcblk_phy      ! all thermodynamics functions, rho_air, q_sat, etc... !LB 
    4040 
    4141   IMPLICIT NONE 
    4242   PRIVATE 
    4343 
    44    PUBLIC ::   TURB_NCAR   ! called by sbcblk.F90 
    45  
    46    !                              ! NCAR own values for given constants: 
    47    REAL(wp), PARAMETER ::   rctv0 = 0.608   ! constant to obtain virtual temperature... 
    48     
     44   PUBLIC :: TURB_NCAR   ! called by sbcblk.F90 
     45 
     46   INTEGER , PARAMETER ::   nb_itt = 5        ! number of itterations 
     47 
    4948   !!---------------------------------------------------------------------- 
    5049CONTAINS 
     
    5352      &                  Cd, Ch, Ce, t_zu, q_zu, U_blk,      & 
    5453      &                  Cdn, Chn, Cen                       ) 
    55       !!---------------------------------------------------------------------------------- 
     54      !!---------------------------------------------------------------------- 
    5655      !!                      ***  ROUTINE  turb_ncar  *** 
    5756      !! 
     
    6160      !!                Returns the effective bulk wind speed at 10m to be used in the bulk formulas 
    6261      !! 
    63       !! ** Method : Monin Obukhov Similarity Theory 
    64       !!             + Large & Yeager (2004,2008) closure: CD_n10 = f(U_n10) 
    65       !! 
    66       !! ** References :   Large & Yeager, 2004 / Large & Yeager, 2008 
    67       !! 
    68       !! ** Last update: Laurent Brodeau, June 2014: 
    69       !!    - handles both cases zt=zu and zt/=zu 
    70       !!    - optimized: less 2D arrays allocated and less operations 
    71       !!    - better first guess of stability by checking air-sea difference of virtual temperature 
    72       !!       rather than temperature difference only... 
    73       !!    - added function "cd_neutral_10m" that uses the improved parametrization of 
    74       !!      Large & Yeager 2008. Drag-coefficient reduction for Cyclone conditions! 
    75       !!    - using code-wide physical constants defined into "phycst.mod" rather than redifining them 
    76       !!      => 'vkarmn' and 'grav' 
    77       !! 
    78       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    7962      !! 
    8063      !! INPUT : 
    8164      !! ------- 
    8265      !!    *  zt   : height for temperature and spec. hum. of air            [m] 
    83       !!    *  zu   : height for wind speed (generally 10m)                   [m] 
    84       !!    *  U_zu : scalar wind speed at 10m                                [m/s] 
    85       !!    *  sst  : SST                                                     [K] 
     66      !!    *  zu   : height for wind speed (usually 10m)                     [m] 
     67      !!    *  sst  : bulk SST                                                [K] 
    8668      !!    *  t_zt : potential air temperature at zt                         [K] 
    8769      !!    *  ssq  : specific humidity at saturation at SST                  [kg/kg] 
    8870      !!    *  q_zt : specific humidity of air at zt                          [kg/kg] 
     71      !!    *  U_zu : scalar wind speed at zu                                 [m/s] 
    8972      !! 
    9073      !! 
     
    9679      !!    *  t_zu   : pot. air temperature adjusted at wind height zu       [K] 
    9780      !!    *  q_zu   : specific humidity of air        //                    [kg/kg] 
    98       !!    *  U_blk  : bulk wind at 10m                                      [m/s] 
     81      !!    *  U_blk  : bulk wind speed at zu                                 [m/s] 
     82      !! 
     83      !! 
     84      !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    9985      !!---------------------------------------------------------------------------------- 
    10086      REAL(wp), INTENT(in   )                     ::   zt       ! height for t_zt and q_zt                    [m] 
     
    10389      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   t_zt     ! potential air temperature              [Kelvin] 
    10490      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] 
     91      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_zt     ! specific air humidity at zt             [kg/kg] 
    10692      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   U_zu     ! relative wind module at zu                [m/s] 
    10793      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cd       ! transfer coefficient for momentum         (tau) 
     
    11096      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   t_zu     ! pot. air temp. adjusted at zu               [K] 
    11197      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] 
     98      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   U_blk    ! bulk wind speed at zu                     [m/s] 
    11399      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients 
    114100      ! 
    115       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 
     101      INTEGER :: j_itt 
     102      LOGICAL :: l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U 
    118103      ! 
    119104      REAL(wp), DIMENSION(jpi,jpj) ::   Cx_n10        ! 10m neutral latent/sensible coefficient 
     
    126111      ! 
    127112      l_zt_equal_zu = .FALSE. 
    128       IF( ABS(zu - zt) < 0.01 ) l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision 
    129  
    130       U_blk = MAX( 0.5 , U_zu )   !  relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s 
     113      IF( ABS(zu - zt) < 0.01_wp )  l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision 
     114 
     115      U_blk = MAX( 0.5_wp , U_zu )   !  relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s 
    131116 
    132117      !! First guess of stability: 
    133       ztmp0 = t_zt*(1. + rctv0*q_zt) - sst*(1. + rctv0*ssq) ! air-sea difference of virtual pot. temp. at zt 
    134       stab  = 0.5 + sign(0.5,ztmp0)                           ! stab = 1 if dTv > 0  => STABLE, 0 if unstable 
     118      ztmp0 = virt_temp(t_zt, q_zt) - virt_temp(sst, ssq) ! air-sea difference of virtual pot. temp. at zt 
     119      stab  = 0.5_wp + sign(0.5_wp,ztmp0)                           ! stab = 1 if dTv > 0  => STABLE, 0 if unstable 
    135120 
    136121      !! Neutral coefficients at 10m: 
     
    139124         ztmp0   (:,:) = cdn_wave(:,:) 
    140125      ELSE 
    141          ztmp0 = cd_neutral_10m( U_blk ) 
     126      ztmp0 = cd_neutral_10m( U_blk ) 
    142127      ENDIF 
    143128 
     
    146131      !! Initializing transf. coeff. with their first guess neutral equivalents : 
    147132      Cd = ztmp0 
    148       Ce = 1.e-3*( 34.6 * sqrt_Cd_n10 ) 
    149       Ch = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1. - stab)) 
     133      Ce = 1.e-3_wp*( 34.6_wp * sqrt_Cd_n10 ) 
     134      Ch = 1.e-3_wp*sqrt_Cd_n10*(18._wp*stab + 32.7_wp*(1._wp - stab)) 
    150135      stab = sqrt_Cd_n10   ! Temporaty array !!! stab == SQRT(Cd) 
    151136  
    152       IF( ln_cdgw )   Cen = Ce  ; Chn = Ch 
     137      IF( ln_cdgw ) THEN 
     138   Cen = Ce 
     139   Chn = Ch 
     140      ENDIF 
    153141 
    154142      !! Initializing values at z_u with z_t values: 
    155143      t_zu = t_zt   ;   q_zu = q_zt 
    156144 
    157       !!  * Now starting iteration loop 
    158       DO j_itt=1, nb_itt 
     145      !! ITERATION BLOCK 
     146      DO j_itt = 1, nb_itt 
    159147         ! 
    160148         ztmp1 = t_zu - sst   ! Updating air/sea differences 
     
    162150 
    163151         ! Updating turbulent scales :   (L&Y 2004 eq. (7)) 
    164          ztmp1  = Ch/stab*ztmp1    ! theta*   (stab == SQRT(Cd)) 
    165          ztmp2  = Ce/stab*ztmp2    ! q*       (stab == SQRT(Cd)) 
    166  
    167          ztmp0 = 1. + rctv0*q_zu      ! multiply this with t and you have the virtual temperature 
     152         ztmp0 = stab*U_blk       ! u*       (stab == SQRT(Cd)) 
     153         ztmp1 = Ch/stab*ztmp1    ! theta*   (stab == SQRT(Cd)) 
     154         ztmp2 = Ce/stab*ztmp2    ! q*       (stab == SQRT(Cd)) 
    168155 
    169156         ! Estimate the inverse of Monin-Obukov length (1/L) at height zu: 
    170          ztmp0 =  (grav*vkarmn/(t_zu*ztmp0)*(ztmp1*ztmp0 + rctv0*t_zu*ztmp2)) / (Cd*U_blk*U_blk) 
    171          !                                                      ( Cd*U_blk*U_blk is U*^2 at zu ) 
    172  
     157         ztmp0 = One_on_L( t_zu, q_zu, ztmp0, ztmp1, ztmp2 ) 
     158          
    173159         !! Stability parameters : 
    174          zeta_u   = zu*ztmp0   ;  zeta_u = sign( min(abs(zeta_u),10.0), zeta_u ) 
     160         zeta_u   = zu*ztmp0 
     161         zeta_u = sign( min(abs(zeta_u),10._wp), zeta_u ) 
    175162         zpsi_h_u = psi_h( zeta_u ) 
    176163 
     
    178165         IF( .NOT. l_zt_equal_zu ) THEN 
    179166            !! Array 'stab' is free for the moment so using it to store 'zeta_t' 
    180             stab = zt*ztmp0 ;  stab = SIGN( MIN(ABS(stab),10.0), stab )  ! Temporaty array stab == zeta_t !!! 
     167            stab = zt*ztmp0 
     168            stab = SIGN( MIN(ABS(stab),10._wp), stab )  ! Temporaty array stab == zeta_t !!! 
    181169            stab = LOG(zt/zu) + zpsi_h_u - psi_h(stab)                   ! stab just used as temp array again! 
    182170            t_zu = t_zt - ztmp1/vkarmn*stab    ! ztmp1 is still theta*  L&Y 2004 eq.(9b) 
    183171            q_zu = q_zt - ztmp2/vkarmn*stab    ! ztmp2 is still q*      L&Y 2004 eq.(9c) 
    184             q_zu = max(0., q_zu) 
    185          END IF 
    186  
     172            q_zu = max(0._wp, q_zu) 
     173         ENDIF 
     174 
     175         ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 eq. 9a)... 
     176         !   In very rare low-wind conditions, the old way of estimating the 
     177         !   neutral wind speed at 10m leads to a negative value that causes the code 
     178         !   to crash. To prevent this a threshold of 0.25m/s is imposed. 
    187179         ztmp2 = psi_m(zeta_u) 
    188180         IF( ln_cdgw ) THEN      ! surface wave case 
    189181            stab = vkarmn / ( vkarmn / sqrt_Cd_n10 - ztmp2 )  ! (stab == SQRT(Cd)) 
    190182            Cd   = stab * stab 
    191             ztmp0 = (LOG(zu/10.) - zpsi_h_u) / vkarmn / sqrt_Cd_n10 
     183            ztmp0 = (LOG(zu/10._wp) - zpsi_h_u) / vkarmn / sqrt_Cd_n10 
    192184            ztmp2 = stab / sqrt_Cd_n10   ! (stab == SQRT(Cd)) 
    193             ztmp1 = 1. + Chn * ztmp0      
     185            ztmp1 = 1._wp + Chn * ztmp0      
    194186            Ch    = Chn * ztmp2 / ztmp1  ! L&Y 2004 eq. (10b) 
    195             ztmp1 = 1. + Cen * ztmp0 
     187            ztmp1 = 1._wp + Cen * ztmp0 
    196188            Ce    = Cen * ztmp2 / ztmp1  ! L&Y 2004 eq. (10c) 
    197189 
    198190         ELSE 
    199             ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 eq. 9a)... 
    200             !   In very rare low-wind conditions, the old way of estimating the 
    201             !   neutral wind speed at 10m leads to a negative value that causes the code 
    202             !   to crash. To prevent this a threshold of 0.25m/s is imposed. 
    203             ztmp0 = MAX( 0.25 , U_blk/(1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - ztmp2)) ) ! U_n10 (ztmp2 == psi_m(zeta_u)) 
    204             ztmp0 = cd_neutral_10m(ztmp0)                                               ! Cd_n10 
    205             Cdn(:,:) = ztmp0 
    206             sqrt_Cd_n10 = sqrt(ztmp0) 
    207  
    208             stab    = 0.5 + sign(0.5,zeta_u)                           ! update stability 
    209             Cx_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1. - stab))  ! L&Y 2004 eq. (6c-6d)    (Cx_n10 == Ch_n10) 
    210             Chn(:,:) = Cx_n10 
    211  
    212             !! Update of transfer coefficients: 
    213             ztmp1 = 1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - ztmp2)   ! L&Y 2004 eq. (10a) (ztmp2 == psi_m(zeta_u)) 
    214             Cd      = ztmp0 / ( ztmp1*ztmp1 ) 
    215             stab = SQRT( Cd ) ! Temporary array !!! (stab == SQRT(Cd)) 
    216  
    217             ztmp0 = (LOG(zu/10.) - zpsi_h_u) / vkarmn / sqrt_Cd_n10 
    218             ztmp2 = stab / sqrt_Cd_n10   ! (stab == SQRT(Cd)) 
    219             ztmp1 = 1. + Cx_n10*ztmp0    ! (Cx_n10 == Ch_n10) 
    220             Ch  = Cx_n10*ztmp2 / ztmp1   ! L&Y 2004 eq. (10b) 
    221  
    222             Cx_n10  = 1.e-3 * (34.6 * sqrt_Cd_n10)  ! L&Y 2004 eq. (6b)    ! Cx_n10 == Ce_n10 
    223             Cen(:,:) = Cx_n10 
    224             ztmp1 = 1. + Cx_n10*ztmp0 
    225             Ce  = Cx_n10*ztmp2 / ztmp1  ! L&Y 2004 eq. (10c) 
    226             ENDIF 
    227          ! 
    228       END DO 
    229       ! 
     191         ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 eq. 9a)... 
     192         !   In very rare low-wind conditions, the old way of estimating the 
     193         !   neutral wind speed at 10m leads to a negative value that causes the code 
     194         !   to crash. To prevent this a threshold of 0.25m/s is imposed. 
     195         ztmp0 = MAX( 0.25_wp , U_blk/(1._wp + sqrt_Cd_n10/vkarmn*(LOG(zu/10._wp) - ztmp2)) ) ! U_n10 (ztmp2 == psi_m(zeta_u)) 
     196         ztmp0 = cd_neutral_10m(ztmp0)                                               ! Cd_n10 
     197         Cdn(:,:) = ztmp0 
     198         sqrt_Cd_n10 = sqrt(ztmp0) 
     199 
     200         stab    = 0.5_wp + sign(0.5_wp,zeta_u)                        ! update stability 
     201         Cx_n10  = 1.e-3_wp*sqrt_Cd_n10*(18._wp*stab + 32.7_wp*(1._wp - stab))  ! L&Y 2004 eq. (6c-6d)    (Cx_n10 == Ch_n10) 
     202         Chn(:,:) = Cx_n10 
     203 
     204         !! Update of transfer coefficients: 
     205         ztmp1 = 1._wp + sqrt_Cd_n10/vkarmn*(LOG(zu/10._wp) - ztmp2)   ! L&Y 2004 eq. (10a) (ztmp2 == psi_m(zeta_u)) 
     206         Cd      = ztmp0 / ( ztmp1*ztmp1 ) 
     207         stab = SQRT( Cd ) ! Temporary array !!! (stab == SQRT(Cd)) 
     208 
     209         ztmp0 = (LOG(zu/10._wp) - zpsi_h_u) / vkarmn / sqrt_Cd_n10 
     210         ztmp2 = stab / sqrt_Cd_n10   ! (stab == SQRT(Cd)) 
     211         ztmp1 = 1._wp + Cx_n10*ztmp0    ! (Cx_n10 == Ch_n10) 
     212         Ch  = Cx_n10*ztmp2 / ztmp1   ! L&Y 2004 eq. (10b) 
     213 
     214         Cx_n10  = 1.e-3_wp * (34.6_wp * sqrt_Cd_n10)  ! L&Y 2004 eq. (6b)    ! Cx_n10 == Ce_n10 
     215         Cen(:,:) = Cx_n10 
     216         ztmp1 = 1._wp + Cx_n10*ztmp0 
     217         Ce  = Cx_n10*ztmp2 / ztmp1  ! L&Y 2004 eq. (10c) 
     218         ENDIF 
     219 
     220      END DO !DO j_itt = 1, nb_itt 
     221 
    230222   END SUBROUTINE turb_ncar 
    231223 
     
    238230      !! Origin: Large & Yeager 2008 eq.(11a) and eq.(11b) 
    239231      !! 
    240       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
     232      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    241233      !!---------------------------------------------------------------------------------- 
    242234      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pw10           ! scalar wind speed at 10m (m/s) 
     
    255247            ! 
    256248            ! When wind speed > 33 m/s => Cyclone conditions => special treatment 
    257             zgt33 = 0.5 + SIGN( 0.5, (zw - 33.) )   ! If pw10 < 33. => 0, else => 1 
    258             ! 
    259             cd_neutral_10m(ji,jj) = 1.e-3 * ( & 
    260                &       (1. - zgt33)*( 2.7/zw + 0.142 + zw/13.09 - 3.14807E-10*zw6) & ! wind <  33 m/s 
    261                &      +    zgt33   *      2.34 )                                     ! wind >= 33 m/s 
    262             ! 
    263             cd_neutral_10m(ji,jj) = MAX(cd_neutral_10m(ji,jj), 1.E-6) 
     249            zgt33 = 0.5_wp + SIGN( 0.5_wp, (zw - 33._wp) )   ! If pw10 < 33. => 0, else => 1 
     250            ! 
     251            cd_neutral_10m(ji,jj) = 1.e-3_wp * ( & 
     252               &       (1._wp - zgt33)*( 2.7_wp/zw + 0.142_wp + zw/13.09_wp - 3.14807E-10_wp*zw6) & ! wind <  33 m/s 
     253               &      +    zgt33   *      2.34_wp )                                                 ! wind >= 33 m/s 
     254            ! 
     255            cd_neutral_10m(ji,jj) = MAX(cd_neutral_10m(ji,jj), 1.E-6_wp) 
    264256            ! 
    265257         END DO 
     
    273265      !! Universal profile stability function for momentum 
    274266      !!    !! Psis, L&Y 2004 eq. (8c), (8d), (8e) 
    275       !!      
    276       !! pzet0 : stability paramenter, z/L where z is altitude measurement                                           
     267      !! 
     268      !! pzeta : stability paramenter, z/L where z is altitude measurement 
    277269      !!         and L is M-O length 
    278270      !! 
    279       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    280       !!---------------------------------------------------------------------------------- 
    281       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pzeta 
    282       REAL(wp), DIMENSION(jpi,jpj)             ::   psi_m 
    283       ! 
    284       INTEGER  ::   ji, jj         ! dummy loop indices 
     271      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
     272      !!---------------------------------------------------------------------------------- 
     273      REAL(wp), DIMENSION(jpi,jpj) :: psi_m 
     274      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta 
     275      ! 
     276      INTEGER  ::   ji, jj    ! dummy loop indices 
    285277      REAL(wp) :: zx2, zx, zstab   ! local scalars 
    286278      !!---------------------------------------------------------------------------------- 
    287       ! 
    288279      DO jj = 1, jpj 
    289280         DO ji = 1, jpi 
    290             zx2 = SQRT( ABS( 1. - 16.*pzeta(ji,jj) ) ) 
    291             zx2 = MAX ( zx2 , 1. ) 
     281            zx2 = SQRT( ABS( 1._wp - 16._wp*pzeta(ji,jj) ) ) 
     282            zx2 = MAX( zx2 , 1._wp ) 
    292283            zx  = SQRT( zx2 ) 
    293             zstab = 0.5 + SIGN( 0.5 , pzeta(ji,jj) ) 
    294             ! 
    295             psi_m(ji,jj) =        zstab  * (-5.*pzeta(ji,jj))       &          ! Stable 
    296                &          + (1. - zstab) * (2.*LOG((1. + zx)*0.5)   &          ! Unstable 
    297                &               + LOG((1. + zx2)*0.5) - 2.*ATAN(zx) + rpi*0.5)  !    " 
     284            zstab = 0.5_wp + SIGN( 0.5_wp , pzeta(ji,jj) ) 
     285            ! 
     286            psi_m(ji,jj) =        zstab  * (-5._wp*pzeta(ji,jj))       &          ! Stable 
     287               &          + (1._wp - zstab) * (2._wp*LOG((1._wp + zx)*0.5_wp)   &          ! Unstable 
     288               &               + LOG((1._wp + zx2)*0.5_wp) - 2._wp*ATAN(zx) + rpi*0.5_wp)  !    " 
    298289            ! 
    299290         END DO 
    300291      END DO 
    301       ! 
    302292   END FUNCTION psi_m 
    303293 
     
    308298      !!    !! Psis, L&Y 2004 eq. (8c), (8d), (8e) 
    309299      !! 
    310       !! pzet0 : stability paramenter, z/L where z is altitude measurement                                           
     300      !! pzeta : stability paramenter, z/L where z is altitude measurement 
    311301      !!         and L is M-O length 
    312302      !! 
    313       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    314       !!---------------------------------------------------------------------------------- 
     303      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
     304      !!---------------------------------------------------------------------------------- 
     305      REAL(wp), DIMENSION(jpi,jpj) :: psi_h 
    315306      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta 
    316       REAL(wp), DIMENSION(jpi,jpj)             :: psi_h 
    317       ! 
    318       INTEGER  ::   ji, jj    ! dummy loop indices 
     307      ! 
     308      INTEGER  ::   ji, jj     ! dummy loop indices 
    319309      REAL(wp) :: zx2, zstab  ! local scalars 
    320310      !!---------------------------------------------------------------------------------- 
     
    322312      DO jj = 1, jpj 
    323313         DO ji = 1, jpi 
    324             zx2 = SQRT( ABS( 1. - 16.*pzeta(ji,jj) ) ) 
    325             zx2 = MAX ( zx2 , 1. ) 
    326             zstab = 0.5 + SIGN( 0.5 , pzeta(ji,jj) ) 
    327             ! 
    328             psi_h(ji,jj) =         zstab  * (-5.*pzeta(ji,jj))        &  ! Stable 
    329                &           + (1. - zstab) * (2.*LOG( (1. + zx2)*0.5 ))   ! Unstable 
     314            zx2 = SQRT( ABS( 1._wp - 16._wp*pzeta(ji,jj) ) ) 
     315            zx2 = MAX( zx2 , 1._wp ) 
     316            zstab = 0.5_wp + SIGN( 0.5_wp , pzeta(ji,jj) ) 
     317            ! 
     318            psi_h(ji,jj) =         zstab  * (-5._wp*pzeta(ji,jj))        &  ! Stable 
     319               &           + (1._wp - zstab) * (2._wp*LOG( (1._wp + zx2)*0.5_wp ))   ! Unstable 
    330320            ! 
    331321         END DO 
    332322      END DO 
    333       ! 
    334323   END FUNCTION psi_h 
    335324 
  • NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/SBC/sbcdcy.F90

    r10425 r12015  
    77   !!   NEMO    2.0  !  2006-02  (S. Masson, G. Madec)  adaptation to NEMO 
    88   !!           3.1  !  2009-07  (J.M. Molines)  adaptation to v3.1 
     9   !!           4.*  !  2019-10  (L. Brodeau)  nothing really new, but the routine 
     10   !!                ! "sbc_dcy_param" has been extracted from old function "sbc_dcy" 
     11   !!                ! => this allows the warm-layer param of COARE3* to know the time 
     12   !!                ! of dawn and dusk even if "ln_dm2dc=.false." (rdawn_dcy & rdusk_dcy 
     13   !!                ! are now public) 
    914   !!---------------------------------------------------------------------- 
    1015 
     
    2227   IMPLICIT NONE 
    2328   PRIVATE 
    24     
     29 
    2530   INTEGER, PUBLIC ::   nday_qsr   !: day when parameters were computed 
    26     
     31 
    2732   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   raa , rbb  , rcc  , rab     ! diurnal cycle parameters 
    28    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rtmd, rdawn, rdusk, rscal   !    -      -       - 
    29    
     33   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rtmd, rscal   !    -      -       - 
     34   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: rdawn_dcy, rdusk_dcy   !    -      -       - 
     35 
    3036   PUBLIC   sbc_dcy        ! routine called by sbc 
     37   PUBLIC   sbc_dcy_param  ! routine used here and called by warm-layer parameterization (sbcblk_skin_coare*) 
    3138 
    3239   !!---------------------------------------------------------------------- 
    3340   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    34    !! $Id$  
     41   !! $Id$ 
    3542   !! Software governed by the CeCILL license (see ./LICENSE) 
    3643   !!---------------------------------------------------------------------- 
    3744CONTAINS 
    3845 
    39       INTEGER FUNCTION sbc_dcy_alloc() 
    40          !!---------------------------------------------------------------------- 
    41          !!                ***  FUNCTION sbc_dcy_alloc  *** 
    42          !!---------------------------------------------------------------------- 
    43          ALLOCATE( raa (jpi,jpj) , rbb  (jpi,jpj) , rcc  (jpi,jpj) , rab  (jpi,jpj) ,     & 
    44             &      rtmd(jpi,jpj) , rdawn(jpi,jpj) , rdusk(jpi,jpj) , rscal(jpi,jpj) , STAT=sbc_dcy_alloc ) 
    45             ! 
    46          CALL mpp_sum ( 'sbcdcy', sbc_dcy_alloc ) 
    47          IF( sbc_dcy_alloc /= 0 )   CALL ctl_stop( 'STOP', 'sbc_dcy_alloc: failed to allocate arrays' ) 
    48       END FUNCTION sbc_dcy_alloc 
     46   INTEGER FUNCTION sbc_dcy_alloc() 
     47      !!---------------------------------------------------------------------- 
     48      !!                ***  FUNCTION sbc_dcy_alloc  *** 
     49      !!---------------------------------------------------------------------- 
     50      ALLOCATE( raa (jpi,jpj) , rbb  (jpi,jpj) , rcc  (jpi,jpj) , rab  (jpi,jpj) ,     & 
     51         &      rtmd(jpi,jpj) , rdawn_dcy(jpi,jpj) , rdusk_dcy(jpi,jpj) , rscal(jpi,jpj) , STAT=sbc_dcy_alloc ) 
     52      ! 
     53      CALL mpp_sum ( 'sbcdcy', sbc_dcy_alloc ) 
     54      IF( sbc_dcy_alloc /= 0 )   CALL ctl_stop( 'STOP', 'sbc_dcy_alloc: failed to allocate arrays' ) 
     55   END FUNCTION sbc_dcy_alloc 
    4956 
    5057 
     
    6067      !! 
    6168      !! reference  : Bernie, DJ, E Guilyardi, G Madec, JM Slingo, and SJ Woolnough, 2007 
    62       !!              Impact of resolving the diurnal cycle in an ocean--atmosphere GCM.  
     69      !!              Impact of resolving the diurnal cycle in an ocean--atmosphere GCM. 
    6370      !!              Part 1: a diurnally forced OGCM. Climate Dynamics 29:6, 575-590. 
    6471      !!---------------------------------------------------------------------- 
    6572      LOGICAL , OPTIONAL          , INTENT(in) ::   l_mask    ! use the routine for night mask computation 
    66       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqsrin    ! input daily QSR flux  
     73      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqsrin    ! input daily QSR flux 
    6774      REAL(wp), DIMENSION(jpi,jpj)             ::   zqsrout   ! output QSR flux with diurnal cycle 
    6875      !! 
    6976      INTEGER  ::   ji, jj                                       ! dummy loop indices 
    7077      INTEGER, DIMENSION(jpi,jpj) :: imask_night ! night mask 
    71       REAL(wp) ::   ztwopi, zinvtwopi, zconvrad  
    7278      REAL(wp) ::   zlo, zup, zlousd, zupusd 
    73       REAL(wp) ::   zdsws, zdecrad, ztx, zsin, zcos 
    74       REAL(wp) ::   ztmp, ztmp1, ztmp2, ztest 
     79      REAL(wp) ::   ztmp, ztmp1, ztmp2 
    7580      REAL(wp) ::   ztmpm, ztmpm1, ztmpm2 
    76       !---------------------------statement functions------------------------ 
    77       REAL(wp) ::   fintegral, pt1, pt2, paaa, pbbb, pccc        ! dummy statement function arguments 
    78       fintegral( pt1, pt2, paaa, pbbb, pccc ) =                         & 
    79          &   paaa * pt2 + zinvtwopi * pbbb * SIN(pccc + ztwopi * pt2)   & 
    80          & - paaa * pt1 - zinvtwopi * pbbb * SIN(pccc + ztwopi * pt1) 
    8181      !!--------------------------------------------------------------------- 
    8282      ! 
    8383      ! Initialization 
    8484      ! -------------- 
    85       ztwopi    = 2._wp * rpi 
    86       zinvtwopi = 1._wp / ztwopi 
    87       zconvrad  = ztwopi / 360._wp 
    88  
    8985      ! When are we during the day (from 0 to 1) 
    9086      zlo = ( REAL(nsec_day, wp) - 0.5_wp * rdt ) / rday 
    9187      zup = zlo + ( REAL(nn_fsbc, wp)     * rdt ) / rday 
    92       !                                           
    93       IF( nday_qsr == -1 ) THEN       ! first time step only   
     88      ! 
     89      IF( nday_qsr == -1 ) THEN       ! first time step only 
    9490         IF(lwp) THEN 
    9591            WRITE(numout,*) 
     
    9894            WRITE(numout,*) 
    9995         ENDIF 
     96      ENDIF 
     97 
     98      ! Setting parameters for each new day: 
     99      CALL sbc_dcy_param() 
     100 
     101      !CALL iom_put( "rdusk_dcy", rdusk_dcy(:,:)*tmask(:,:,1) ) !LB 
     102      !CALL iom_put( "rdawn_dcy", rdawn_dcy(:,:)*tmask(:,:,1) ) !LB 
     103      !CALL iom_put( "rscal_dcy", rscal(:,:)*tmask(:,:,1) ) !LB 
     104 
     105 
     106      !     3. update qsr with the diurnal cycle 
     107      !     ------------------------------------ 
     108 
     109      imask_night(:,:) = 0 
     110      DO jj = 1, jpj 
     111         DO ji = 1, jpi 
     112            ztmpm = 0._wp 
     113            IF( ABS(rab(ji,jj)) < 1. ) THEN         ! day duration is less than 24h 
     114               ! 
     115               IF( rdawn_dcy(ji,jj) < rdusk_dcy(ji,jj) ) THEN       ! day time in one part 
     116                  zlousd = MAX(zlo, rdawn_dcy(ji,jj)) 
     117                  zlousd = MIN(zlousd, zup) 
     118                  zupusd = MIN(zup, rdusk_dcy(ji,jj)) 
     119                  zupusd = MAX(zupusd, zlo) 
     120                  ztmp = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
     121                  zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj) 
     122                  ztmpm = zupusd - zlousd 
     123                  IF( ztmpm .EQ. 0 ) imask_night(ji,jj) = 1 
     124                  ! 
     125               ELSE                                         ! day time in two parts 
     126                  zlousd = MIN(zlo, rdusk_dcy(ji,jj)) 
     127                  zupusd = MIN(zup, rdusk_dcy(ji,jj)) 
     128                  ztmp1 = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
     129                  ztmpm1=zupusd-zlousd 
     130                  zlousd = MAX(zlo, rdawn_dcy(ji,jj)) 
     131                  zupusd = MAX(zup, rdawn_dcy(ji,jj)) 
     132                  ztmp2 = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
     133                  ztmpm2 =zupusd-zlousd 
     134                  ztmp = ztmp1 + ztmp2 
     135                  ztmpm = ztmpm1 + ztmpm2 
     136                  zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj) 
     137                  IF(ztmpm .EQ. 0.) imask_night(ji,jj) = 1 
     138               ENDIF 
     139            ELSE                                   ! 24h light or 24h night 
     140               ! 
     141               IF( raa(ji,jj) > rbb(ji,jj) ) THEN           ! 24h day 
     142                  ztmp = fintegral(zlo, zup, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
     143                  zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj) 
     144                  imask_night(ji,jj) = 0 
     145                  ! 
     146               ELSE                                         ! No day 
     147                  zqsrout(ji,jj) = 0.0_wp 
     148                  imask_night(ji,jj) = 1 
     149               ENDIF 
     150            ENDIF 
     151         END DO 
     152      END DO 
     153      ! 
     154      IF( PRESENT(l_mask) .AND. l_mask ) THEN 
     155         zqsrout(:,:) = float(imask_night(:,:)) 
     156      ENDIF 
     157      ! 
     158   END FUNCTION sbc_dcy 
     159 
     160 
     161   SUBROUTINE sbc_dcy_param( ) 
     162      !! 
     163      INTEGER  ::   ji, jj                                       ! dummy loop indices 
     164      !INTEGER, DIMENSION(jpi,jpj) :: imask_night ! night mask 
     165      REAL(wp) ::   zdsws, zdecrad, ztx, zsin, zcos 
     166      REAL(wp) ::   ztmp, ztest 
     167      !---------------------------statement functions------------------------ 
     168      ! 
     169      IF( nday_qsr == -1 ) THEN       ! first time step only 
    100170         ! allocate sbcdcy arrays 
    101171         IF( sbc_dcy_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_dcy_alloc : unable to allocate arrays' ) 
    102172         ! Compute rcc needed to compute the time integral of the diurnal cycle 
    103          rcc(:,:) = zconvrad * glamt(:,:) - rpi 
     173         rcc(:,:) = rad * glamt(:,:) - rpi 
    104174         ! time of midday 
    105175         rtmd(:,:) = 0.5_wp - glamt(:,:) / 360._wp 
     
    107177      ENDIF 
    108178 
    109       ! If this is a new day, we have to update the dawn, dusk and scaling function   
     179      ! If this is a new day, we have to update the dawn, dusk and scaling function 
    110180      !---------------------- 
    111      
    112       !     2.1 dawn and dusk   
    113  
    114       ! nday is the number of days since the beginning of the current month  
    115       IF( nday_qsr /= nday ) THEN  
     181 
     182      !     2.1 dawn and dusk 
     183 
     184      ! nday is the number of days since the beginning of the current month 
     185      IF( nday_qsr /= nday ) THEN 
    116186         ! save the day of the year and the daily mean of qsr 
    117          nday_qsr = nday  
    118          ! number of days since the previous winter solstice (supposed to be always 21 December)          
     187         nday_qsr = nday 
     188         ! number of days since the previous winter solstice (supposed to be always 21 December) 
    119189         zdsws = REAL(11 + nday_year, wp) 
    120190         ! declination of the earths orbit 
    121          zdecrad = (-23.5_wp * zconvrad) * COS( zdsws * ztwopi / REAL(nyear_len(1),wp) ) 
     191         zdecrad = (-23.5_wp * rad) * COS( zdsws * 2._wp*rpi / REAL(nyear_len(1),wp) ) 
    122192         ! Compute A and B needed to compute the time integral of the diurnal cycle 
    123193 
     
    125195         DO jj = 1, jpj 
    126196            DO ji = 1, jpi 
    127                ztmp = zconvrad * gphit(ji,jj) 
     197               ztmp = rad * gphit(ji,jj) 
    128198               raa(ji,jj) = SIN( ztmp ) * zsin 
    129199               rbb(ji,jj) = COS( ztmp ) * zcos 
    130             END DO   
    131          END DO   
     200            END DO 
     201         END DO 
    132202         ! Compute the time of dawn and dusk 
    133203 
    134          ! rab to test if the day time is equal to 0, less than 24h of full day         
     204         ! rab to test if the day time is equal to 0, less than 24h of full day 
    135205         rab(:,:) = -raa(:,:) / rbb(:,:) 
    136206         DO jj = 1, jpj 
    137207            DO ji = 1, jpi 
    138                IF ( ABS(rab(ji,jj)) < 1._wp ) THEN         ! day duration is less than 24h 
    139          ! When is it night? 
    140                   ztx = zinvtwopi * (ACOS(rab(ji,jj)) - rcc(ji,jj)) 
    141                   ztest = -rbb(ji,jj) * SIN( rcc(ji,jj) + ztwopi * ztx ) 
    142          ! is it dawn or dusk? 
    143                   IF ( ztest > 0._wp ) THEN 
    144                      rdawn(ji,jj) = ztx 
    145                      rdusk(ji,jj) = rtmd(ji,jj) + ( rtmd(ji,jj) - rdawn(ji,jj) ) 
     208               IF( ABS(rab(ji,jj)) < 1._wp ) THEN         ! day duration is less than 24h 
     209                  ! When is it night? 
     210                  ztx = 1._wp/(2._wp*rpi) * (ACOS(rab(ji,jj)) - rcc(ji,jj)) 
     211                  ztest = -rbb(ji,jj) * SIN( rcc(ji,jj) + 2._wp*rpi * ztx ) 
     212                  ! is it dawn or dusk? 
     213                  IF( ztest > 0._wp ) THEN 
     214                     rdawn_dcy(ji,jj) = ztx 
     215                     rdusk_dcy(ji,jj) = rtmd(ji,jj) + ( rtmd(ji,jj) - rdawn_dcy(ji,jj) ) 
    146216                  ELSE 
    147                      rdusk(ji,jj) = ztx 
    148                      rdawn(ji,jj) = rtmd(ji,jj) - ( rdusk(ji,jj) - rtmd(ji,jj) ) 
     217                     rdusk_dcy(ji,jj) = ztx 
     218                     rdawn_dcy(ji,jj) = rtmd(ji,jj) - ( rdusk_dcy(ji,jj) - rtmd(ji,jj) ) 
    149219                  ENDIF 
    150220               ELSE 
    151                   rdawn(ji,jj) = rtmd(ji,jj) + 0.5_wp 
    152                   rdusk(ji,jj) = rdawn(ji,jj) 
    153                ENDIF 
    154              END DO   
    155          END DO   
    156          rdawn(:,:) = MOD( (rdawn(:,:) + 1._wp), 1._wp ) 
    157          rdusk(:,:) = MOD( (rdusk(:,:) + 1._wp), 1._wp ) 
     221                  rdawn_dcy(ji,jj) = rtmd(ji,jj) + 0.5_wp 
     222                  rdusk_dcy(ji,jj) = rdawn_dcy(ji,jj) 
     223               ENDIF 
     224            END DO 
     225         END DO 
     226         rdawn_dcy(:,:) = MOD( (rdawn_dcy(:,:) + 1._wp), 1._wp ) 
     227         rdusk_dcy(:,:) = MOD( (rdusk_dcy(:,:) + 1._wp), 1._wp ) 
    158228         !     2.2 Compute the scaling function: 
    159229         !         S* = the inverse of the time integral of the diurnal cycle from dawn to dusk 
     
    162232         DO jj = 1, jpj 
    163233            DO ji = 1, jpi 
    164                IF ( ABS(rab(ji,jj)) < 1._wp ) THEN         ! day duration is less than 24h 
     234               IF( ABS(rab(ji,jj)) < 1._wp ) THEN         ! day duration is less than 24h 
    165235                  rscal(ji,jj) = 0.0_wp 
    166                   IF ( rdawn(ji,jj) < rdusk(ji,jj) ) THEN      ! day time in one part 
    167                      IF( (rdusk(ji,jj) - rdawn(ji,jj) ) .ge. 0.001_wp ) THEN 
    168                        rscal(ji,jj) = fintegral(rdawn(ji,jj), rdusk(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj))  
    169                        rscal(ji,jj) = 1._wp / rscal(ji,jj) 
     236                  IF( rdawn_dcy(ji,jj) < rdusk_dcy(ji,jj) ) THEN      ! day time in one part 
     237                     IF( (rdusk_dcy(ji,jj) - rdawn_dcy(ji,jj) ) .ge. 0.001_wp ) THEN 
     238                        rscal(ji,jj) = fintegral(rdawn_dcy(ji,jj), rdusk_dcy(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
     239                        rscal(ji,jj) = 1._wp / rscal(ji,jj) 
    170240                     ENDIF 
    171241                  ELSE                                         ! day time in two parts 
    172                      IF( (rdusk(ji,jj) + (1._wp - rdawn(ji,jj)) ) .ge. 0.001_wp ) THEN 
    173                        rscal(ji,jj) = fintegral(0._wp, rdusk(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj))   & 
    174                           &         + fintegral(rdawn(ji,jj), 1._wp, raa(ji,jj), rbb(ji,jj), rcc(ji,jj))  
    175                        rscal(ji,jj) = 1. / rscal(ji,jj) 
     242                     IF( (rdusk_dcy(ji,jj) + (1._wp - rdawn_dcy(ji,jj)) ) .ge. 0.001_wp ) THEN 
     243                        rscal(ji,jj) = fintegral(0._wp, rdusk_dcy(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj))   & 
     244                           &         + fintegral(rdawn_dcy(ji,jj), 1._wp, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
     245                        rscal(ji,jj) = 1. / rscal(ji,jj) 
    176246                     ENDIF 
    177247                  ENDIF 
    178248               ELSE 
    179                   IF ( raa(ji,jj) > rbb(ji,jj) ) THEN         ! 24h day 
    180                      rscal(ji,jj) = fintegral(0._wp, 1._wp, raa(ji,jj), rbb(ji,jj), rcc(ji,jj))  
     249                  IF( raa(ji,jj) > rbb(ji,jj) ) THEN         ! 24h day 
     250                     rscal(ji,jj) = fintegral(0._wp, 1._wp, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
    181251                     rscal(ji,jj) = 1._wp / rscal(ji,jj) 
    182252                  ELSE                                          ! No day 
     
    184254                  ENDIF 
    185255               ENDIF 
    186             END DO   
    187          END DO   
     256            END DO 
     257         END DO 
    188258         ! 
    189259         ztmp = rday / ( rdt * REAL(nn_fsbc, wp) ) 
    190260         rscal(:,:) = rscal(:,:) * ztmp 
    191261         ! 
    192       ENDIF  
    193          !     3. update qsr with the diurnal cycle 
    194          !     ------------------------------------ 
    195  
    196       imask_night(:,:) = 0 
    197       DO jj = 1, jpj 
    198          DO ji = 1, jpi 
    199             ztmpm = 0._wp 
    200             IF( ABS(rab(ji,jj)) < 1. ) THEN         ! day duration is less than 24h 
    201                ! 
    202                IF( rdawn(ji,jj) < rdusk(ji,jj) ) THEN       ! day time in one part 
    203                   zlousd = MAX(zlo, rdawn(ji,jj)) 
    204                   zlousd = MIN(zlousd, zup) 
    205                   zupusd = MIN(zup, rdusk(ji,jj)) 
    206                   zupusd = MAX(zupusd, zlo) 
    207                   ztmp = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj))  
    208                   zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj) 
    209                   ztmpm = zupusd - zlousd 
    210                   IF ( ztmpm .EQ. 0 ) imask_night(ji,jj) = 1 
    211                   ! 
    212                ELSE                                         ! day time in two parts 
    213                   zlousd = MIN(zlo, rdusk(ji,jj)) 
    214                   zupusd = MIN(zup, rdusk(ji,jj)) 
    215                   ztmp1 = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj))  
    216                   ztmpm1=zupusd-zlousd 
    217                   zlousd = MAX(zlo, rdawn(ji,jj)) 
    218                   zupusd = MAX(zup, rdawn(ji,jj)) 
    219                   ztmp2 = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj))  
    220                   ztmpm2 =zupusd-zlousd 
    221                   ztmp = ztmp1 + ztmp2 
    222                   ztmpm = ztmpm1 + ztmpm2 
    223                   zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj) 
    224                   IF (ztmpm .EQ. 0.) imask_night(ji,jj) = 1 
    225                ENDIF 
    226             ELSE                                   ! 24h light or 24h night 
    227                ! 
    228                IF( raa(ji,jj) > rbb(ji,jj) ) THEN           ! 24h day 
    229                   ztmp = fintegral(zlo, zup, raa(ji,jj), rbb(ji,jj), rcc(ji,jj))  
    230                   zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj) 
    231                   imask_night(ji,jj) = 0 
    232                   ! 
    233                ELSE                                         ! No day 
    234                   zqsrout(ji,jj) = 0.0_wp 
    235                   imask_night(ji,jj) = 1 
    236                ENDIF 
    237             ENDIF 
    238          END DO   
    239       END DO   
    240       ! 
    241       IF( PRESENT(l_mask) .AND. l_mask ) THEN 
    242          zqsrout(:,:) = float(imask_night(:,:)) 
    243       ENDIF 
    244       ! 
    245    END FUNCTION sbc_dcy 
     262      ENDIF !IF( nday_qsr /= nday ) 
     263      ! 
     264   END SUBROUTINE sbc_dcy_param 
     265 
     266 
     267   FUNCTION fintegral( pt1, pt2, paaa, pbbb, pccc ) 
     268      REAL(wp), INTENT(in) :: pt1, pt2, paaa, pbbb, pccc 
     269      REAL(wp) :: fintegral 
     270      fintegral =   paaa * pt2 + 1._wp/(2._wp*rpi) * pbbb * SIN(pccc + 2._wp*rpi*pt2)   & 
     271         &        - paaa * pt1 - 1._wp/(2._wp*rpi) * pbbb * SIN(pccc + 2._wp*rpi*pt1) 
     272   END FUNCTION fintegral 
    246273 
    247274   !!====================================================================== 
  • NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/SBC/sbcmod.F90

    r11874 r12015  
    128128      ENDIF 
    129129#else 
    130       IF( lk_si3  )   nn_ice      = 2 
     130      !IF( lk_si3  )   nn_ice      = 2 
    131131      IF( lk_cice )   nn_ice      = 3 
    132132#endif 
     
    261261 
    262262      !                          ! Choice of the Surface Boudary Condition (set nsbc) 
     263      nday_qsr = -1   ! allow initialization at the 1st call !LB: now warm-layer of COARE* calls "sbc_dcy_param" of sbcdcy.F90! 
    263264      IF( ln_dm2dc ) THEN           !* daily mean to diurnal cycle 
    264          nday_qsr = -1   ! allow initialization at the 1st call 
     265         !LB:nday_qsr = -1   ! allow initialization at the 1st call 
    265266         IF( .NOT.( ln_flx .OR. ln_blk .OR. ln_abl ) .AND. nn_components /= jp_iam_opa )   & 
    266267            &   CALL ctl_stop( 'qsr diurnal cycle from daily values requires flux, bulk or abl formulation' ) 
     
    370371         CALL iom_set_rstw_var_active('sfx_b') 
    371372      ENDIF 
    372       !  
     373 
    373374   END SUBROUTINE sbc_init 
    374375 
     
    409410         emp_b (:,:) = emp (:,:) 
    410411         sfx_b (:,:) = sfx (:,:) 
    411          IF ( ln_rnf ) THEN 
     412         IF( ln_rnf ) THEN 
    412413            rnf_b    (:,:  ) = rnf    (:,:  ) 
    413414            rnf_tsc_b(:,:,:) = rnf_tsc(:,:,:) 
     
    451452      IF( ln_mixcpl )          CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice )   ! forced-coupled mixed formulation after forcing 
    452453      ! 
    453       IF ( ln_wave .AND. (ln_tauwoc .OR. ln_tauw) ) CALL sbc_wstress( )      ! Wind stress provided by waves  
     454      IF( ln_wave .AND. (ln_tauwoc .OR. ln_tauw) ) CALL sbc_wstress( )      ! Wind stress provided by waves  
    454455      ! 
    455456      !                                            !==  Misc. Options  ==! 
     
    485486!!$!clem: it looks like it is necessary for the north fold (in certain circumstances). Don't know why. 
    486487!!$      CALL lbc_lnk( 'sbcmod', emp, 'T', 1. ) 
    487       IF ( ll_wd ) THEN     ! If near WAD point limit the flux for now 
     488      IF( ll_wd ) THEN     ! If near WAD point limit the flux for now 
    488489         zthscl = atanh(rn_wd_sbcfra)                     ! taper frac default is .999  
    489490         zwdht(:,:) = sshn(:,:) + ht_0(:,:) - rn_wdmin1   ! do this calc of water 
Note: See TracChangeset for help on using the changeset viewer.