New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 12182 for NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/SBC/sbcblk.F90 – NEMO

Ignore:
Timestamp:
2019-12-11T12:38:43+01:00 (4 years ago)
Author:
davestorkey
Message:

2019/dev_r11943_MERGE_2019: Merge in dev_ASINTER-01-05_merge.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/SBC/sbcblk.F90

    r11960 r12182  
    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) 
     20   !!            4.0  !  2019-03  (F. Lemarié & G. Samson)  add ABL compatibility (ln_abl=TRUE) 
    2021   !!---------------------------------------------------------------------- 
    2122 
     
    2324   !!   sbc_blk_init  : initialisation of the chosen bulk formulation as ocean surface boundary condition 
    2425   !!   sbc_blk       : bulk formulation as ocean surface boundary condition 
    25    !!   blk_oce       : computes momentum, heat and freshwater fluxes over ocean 
    26    !!   rho_air       : density of (moist) air (depends on T_air, q_air and SLP 
    27    !!   cp_air        : specific heat of (moist) air (depends spec. hum. q_air) 
    28    !!   q_sat         : saturation humidity as a function of SLP and temperature 
    29    !!   L_vap         : latent heat of vaporization of water as a function of temperature 
    30    !!             sea-ice case only :  
    31    !!   blk_ice_tau   : provide the air-ice stress 
    32    !!   blk_ice_flx   : provide the heat and mass fluxes at air-ice interface 
     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 : 
     29   !!   blk_ice_1   : provide the air-ice stress 
     30   !!   blk_ice_2   : provide the heat and mass fluxes at air-ice interface 
    3331   !!   blk_ice_qcn   : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux) 
    3432   !!   Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag 
    35    !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag  
     33   !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag 
    3634   !!---------------------------------------------------------------------- 
    3735   USE oce            ! ocean dynamics and tracers 
     
    4644   USE lib_fortran    ! to use key_nosignedzero 
    4745#if defined key_si3 
    48    USE ice     , ONLY :   u_ice, v_ice, jpl, a_i_b, at_i_b, t_su, rn_cnd_s, hfx_err_dif 
     46   USE ice     , ONLY :   jpl, a_i_b, at_i_b, rn_cnd_s, hfx_err_dif 
    4947   USE icethd_dh      ! for CALL ice_thd_snwblow 
    5048#endif 
    51    USE sbcblk_algo_ncar     ! => turb_ncar     : NCAR - CORE (Large & Yeager, 2009)  
    52    USE sbcblk_algo_coare    ! => turb_coare    : COAREv3.0 (Fairall et al. 2003)  
    53    USE sbcblk_algo_coare3p5 ! => turb_coare3p5 : COAREv3.5 (Edson et al. 2013) 
    54    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) 
    5553   ! 
    5654   USE iom            ! I/O manager library 
     
    6058   USE prtctl         ! Print control 
    6159 
     60   USE sbcblk_phy     ! a catalog of functions for physical/meteorological parameters in the marine boundary layer, rho_air, q_sat, etc... 
     61 
     62 
    6263   IMPLICIT NONE 
    6364   PRIVATE 
     
    6566   PUBLIC   sbc_blk_init  ! called in sbcmod 
    6667   PUBLIC   sbc_blk       ! called in sbcmod 
     68   PUBLIC   blk_oce_1     ! called in sbcabl 
     69   PUBLIC   blk_oce_2     ! called in sbcabl 
    6770#if defined key_si3 
    68    PUBLIC   blk_ice_tau   ! routine called in icesbc 
    69    PUBLIC   blk_ice_flx   ! routine called in icesbc 
     71   PUBLIC   blk_ice_   ! routine called in icesbc 
     72   PUBLIC   blk_ice_   ! routine called in icesbc 
    7073   PUBLIC   blk_ice_qcn   ! routine called in icesbc 
    71 #endif  
    72  
    73 !!Lolo: should ultimately be moved in the module with all physical constants ? 
    74 !!gm  : In principle, yes. 
    75    REAL(wp), PARAMETER ::   Cp_dry = 1005.0       !: Specic heat of dry air, constant pressure      [J/K/kg] 
    76    REAL(wp), PARAMETER ::   Cp_vap = 1860.0       !: Specic heat of water vapor, constant pressure  [J/K/kg] 
    77    REAL(wp), PARAMETER ::   R_dry = 287.05_wp     !: Specific gas constant for dry air              [J/K/kg] 
    78    REAL(wp), PARAMETER ::   R_vap = 461.495_wp    !: Specific gas constant for water vapor          [J/K/kg] 
    79    REAL(wp), PARAMETER ::   reps0 = R_dry/R_vap   !: ratio of gas constant for dry air and water vapor => ~ 0.622 
    80    REAL(wp), PARAMETER ::   rctv0 = R_vap/R_dry   !: for virtual temperature (== (1-eps)/eps) => ~ 0.608 
    81  
    82    INTEGER , PARAMETER ::   jpfld   =10           ! maximum number of files to read 
    83    INTEGER , PARAMETER ::   jp_wndi = 1           ! index of 10m wind velocity (i-component) (m/s)    at T-point 
    84    INTEGER , PARAMETER ::   jp_wndj = 2           ! index of 10m wind velocity (j-component) (m/s)    at T-point 
    85    INTEGER , PARAMETER ::   jp_tair = 3           ! index of 10m air temperature             (Kelvin) 
    86    INTEGER , PARAMETER ::   jp_humi = 4           ! index of specific humidity               ( % ) 
    87    INTEGER , PARAMETER ::   jp_qsr  = 5           ! index of solar heat                      (W/m2) 
    88    INTEGER , PARAMETER ::   jp_qlw  = 6           ! index of Long wave                       (W/m2) 
    89    INTEGER , PARAMETER ::   jp_prec = 7           ! index of total precipitation (rain+snow) (Kg/m2/s) 
    90    INTEGER , PARAMETER ::   jp_snow = 8           ! index of snow (solid prcipitation)       (kg/m2/s) 
    91    INTEGER , PARAMETER ::   jp_slp  = 9           ! index of sea level pressure              (Pa) 
    92    INTEGER , PARAMETER ::   jp_tdif =10           ! index of tau diff associated to HF tau   (N/m2)   at T-point 
    93  
    94    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf   ! structure of input fields (file informations, fields read) 
    95  
    96    !                                             !!! Bulk parameters 
    97    REAL(wp), PARAMETER ::   cpa    = 1000.5         ! specific heat of air (only used for ice fluxes now...) 
    98    REAL(wp), PARAMETER ::   Ls     =    2.839e6     ! latent heat of sublimation 
    99    REAL(wp), PARAMETER ::   Stef   =    5.67e-8     ! Stefan Boltzmann constant 
    100    REAL(wp), PARAMETER ::   Cd_ice =    1.4e-3      ! transfer coefficient over ice 
    101    REAL(wp), PARAMETER ::   albo   =    0.066       ! ocean albedo assumed to be constant 
    102    ! 
     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 
    10391   !                           !!* Namelist namsbc_blk : bulk parameters 
    10492   LOGICAL  ::   ln_NCAR        ! "NCAR"      algorithm   (Large and Yeager 2008) 
    10593   LOGICAL  ::   ln_COARE_3p0   ! "COARE 3.0" algorithm   (Fairall et al. 2003) 
    106    LOGICAL  ::   ln_COARE_3p5   ! "COARE 3.5" algorithm   (Edson et al. 2013) 
    107    LOGICAL  ::   ln_ECMWF       ! "ECMWF"     algorithm   (IFS cycle 31) 
     94   LOGICAL  ::   ln_COARE_3p6   ! "COARE 3.6" algorithm   (Edson et al. 2013) 
     95   LOGICAL  ::   ln_ECMWF       ! "ECMWF"     algorithm   (IFS cycle 45r1) 
    10896   ! 
    109    LOGICAL  ::   ln_taudif      ! logical flag to use the "mean of stress module - module of mean stress" data 
    110    REAL(wp) ::   rn_pfac        ! multiplication factor for precipitation 
    111    REAL(wp) ::   rn_efac        ! multiplication factor for evaporation 
    112    REAL(wp) ::   rn_vfac        ! multiplication factor for ice/ocean velocity in the calculation of wind stress 
    113    REAL(wp) ::   rn_zqt         ! z(q,t) : height of humidity and temperature measurements 
    114    REAL(wp) ::   rn_zu          ! z(u)   : height of wind measurements 
    115 !!gm ref namelist initialize it so remove the setting to false below 
    116    LOGICAL  ::   ln_Cd_L12 = .FALSE. !  Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2012) 
    117    LOGICAL  ::   ln_Cd_L15 = .FALSE. !  Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2015) 
     97   LOGICAL  ::   ln_Cd_L12      ! ice-atm drag = F( ice concentration )                        (Lupkes et al. JGR2012) 
     98   LOGICAL  ::   ln_Cd_L15      ! ice-atm drag = F( ice concentration, atmospheric stability ) (Lupkes et al. JGR2015) 
    11899   ! 
    119    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Cd_atm                    ! transfer coefficient for momentum      (tau) 
    120    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Ch_atm                    ! transfer coefficient for sensible heat (Q_sens) 
    121    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Ce_atm                    ! tansfert coefficient for evaporation   (Q_lat) 
    122    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_zu                      ! air temperature at wind speed height (needed by Lupkes 2015 bulk scheme) 
    123    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   q_zu                      ! air spec. hum.  at wind speed height (needed by Lupkes 2015 bulk scheme) 
    124    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   cdn_oce, chn_oce, cen_oce ! needed by Lupkes 2015 bulk scheme 
     100   REAL(wp)         ::   rn_pfac   ! multiplication factor for precipitation 
     101   REAL(wp), PUBLIC ::   rn_efac   ! multiplication factor for evaporation 
     102   REAL(wp), PUBLIC ::   rn_vfac   ! multiplication factor for ice/ocean velocity in the calculation of wind stress 
     103   REAL(wp)         ::   rn_zqt    ! z(q,t) : height of humidity and temperature measurements 
     104   REAL(wp)         ::   rn_zu     ! z(u)   : height of wind measurements 
     105   ! 
     106   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   Cd_ice , Ch_ice , Ce_ice   ! transfert coefficients over ice 
     107   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   Cdn_oce, Chn_oce, Cen_oce  ! neutral coeffs over ocean (L15 bulk scheme) 
     108   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   t_zu, q_zu                 ! air temp. and spec. hum. at wind speed height (L15 bulk scheme) 
     109 
     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 
    125121 
    126122   INTEGER  ::   nblk           ! choice of the bulk algorithm 
     
    128124   INTEGER, PARAMETER ::   np_NCAR      = 1   ! "NCAR" algorithm        (Large and Yeager 2008) 
    129125   INTEGER, PARAMETER ::   np_COARE_3p0 = 2   ! "COARE 3.0" algorithm   (Fairall et al. 2003) 
    130    INTEGER, PARAMETER ::   np_COARE_3p5 = 3   ! "COARE 3.5" algorithm   (Edson et al. 2013) 
    131    INTEGER, PARAMETER ::   np_ECMWF     = 4   ! "ECMWF" algorithm       (IFS cycle 31) 
     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) 
    132128 
    133129   !! * Substitutions 
     
    144140      !!             ***  ROUTINE sbc_blk_alloc *** 
    145141      !!------------------------------------------------------------------- 
    146       ALLOCATE( Cd_atm (jpi,jpj), Ch_atm (jpi,jpj), Ce_atm (jpi,jpj), t_zu(jpi,jpj), q_zu(jpi,jpj), & 
    147          &      cdn_oce(jpi,jpj), chn_oce(jpi,jpj), cen_oce(jpi,jpj), STAT=sbc_blk_alloc ) 
     142      ALLOCATE( t_zu(jpi,jpj)   , q_zu(jpi,jpj)   ,                                      & 
     143         &      Cdn_oce(jpi,jpj), Chn_oce(jpi,jpj), Cen_oce(jpi,jpj),                    & 
     144         &      Cd_ice (jpi,jpj), Ch_ice (jpi,jpj), Ce_ice (jpi,jpj), STAT=sbc_blk_alloc ) 
    148145      ! 
    149146      CALL mpp_sum ( 'sbcblk', sbc_blk_alloc ) 
     
    158155      !! ** Purpose :   choose and initialize a bulk formulae formulation 
    159156      !! 
    160       !! ** Method  :  
     157      !! ** Method  : 
    161158      !! 
    162159      !!---------------------------------------------------------------------- 
    163       INTEGER  ::   ifpr, jfld            ! dummy loop indice and argument 
     160      INTEGER  ::   jfpr                  ! dummy loop indice and argument 
    164161      INTEGER  ::   ios, ierror, ioptio   ! Local integer 
    165162      !! 
    166163      CHARACTER(len=100)            ::   cn_dir                ! Root directory for location of atmospheric forcing files 
    167       TYPE(FLD_N), DIMENSION(jpfld) ::   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 
    168165      TYPE(FLD_N) ::   sn_wndi, sn_wndj, sn_humi, sn_qsr       ! informations about the fields to be read 
    169166      TYPE(FLD_N) ::   sn_qlw , sn_tair, sn_prec, sn_snow      !       "                        " 
    170       TYPE(FLD_N) ::   sn_slp , sn_tdif                        !       "                        " 
     167      TYPE(FLD_N) ::   sn_slp , sn_hpgi, sn_hpgj               !       "                        " 
    171168      NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw ,                &   ! input fields 
    172          &                 sn_tair, sn_prec, sn_snow, sn_slp, sn_tdif,                & 
    173          &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF,             &   ! bulk algorithm 
    174          &                 cn_dir , ln_taudif, rn_zqt, rn_zu,                         &  
    175          &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15 
     169         &                 sn_tair, sn_prec, sn_snow, sn_slp, sn_hpgi, sn_hpgj,       & 
     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 
    176174      !!--------------------------------------------------------------------- 
    177175      ! 
     
    179177      IF( sbc_blk_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard arrays' ) 
    180178      ! 
    181       !                             !** read bulk namelist   
     179      !                             !** read bulk namelist 
    182180      READ  ( numnam_ref, namsbc_blk, IOSTAT = ios, ERR = 901) 
    183181901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc_blk in reference namelist' ) 
     
    190188      !                             !** initialization of the chosen bulk formulae (+ check) 
    191189      !                                   !* select the bulk chosen in the namelist and check the choice 
    192                                                                ioptio = 0 
    193       IF( ln_NCAR      ) THEN   ;   nblk =  np_NCAR        ;   ioptio = ioptio + 1   ;   ENDIF 
    194       IF( ln_COARE_3p0 ) THEN   ;   nblk =  np_COARE_3p0   ;   ioptio = ioptio + 1   ;   ENDIF 
    195       IF( ln_COARE_3p5 ) THEN   ;   nblk =  np_COARE_3p5   ;   ioptio = ioptio + 1   ;   ENDIF 
    196       IF( ln_ECMWF     ) THEN   ;   nblk =  np_ECMWF       ;   ioptio = ioptio + 1   ;   ENDIF 
    197       ! 
     190      ioptio = 0 
     191      IF( ln_NCAR      ) THEN 
     192         nblk =  np_NCAR        ;   ioptio = ioptio + 1 
     193      ENDIF 
     194      IF( ln_COARE_3p0 ) THEN 
     195         nblk =  np_COARE_3p0   ;   ioptio = ioptio + 1 
     196      ENDIF 
     197      IF( ln_COARE_3p6 ) THEN 
     198         nblk =  np_COARE_3p6   ;   ioptio = ioptio + 1 
     199      ENDIF 
     200      IF( ln_ECMWF     ) THEN 
     201         nblk =  np_ECMWF       ;   ioptio = ioptio + 1 
     202      ENDIF 
    198203      IF( ioptio /= 1 )   CALL ctl_stop( 'sbc_blk_init: Choose one and only one bulk algorithm' ) 
     204 
     205      !                             !** initialization of the cool-skin / warm-layer parametrization 
     206      IF( ln_skin_cs .OR. ln_skin_wl ) THEN 
     207         !! Some namelist sanity tests: 
     208         IF( ln_NCAR )      & 
     209            & CALL ctl_stop( 'sbc_blk_init: Cool-skin/warm-layer param. not compatible with NCAR algorithm' ) 
     210         IF( nn_fsbc /= 1 ) & 
     211            & CALL ctl_stop( 'sbc_blk_init: Please set "nn_fsbc" to 1 when using cool-skin/warm-layer param.') 
     212      END IF 
     213 
     214      IF( ln_skin_wl ) THEN 
     215         !! Check if the frequency of downwelling solar flux input makes sense and if ln_dm2dc=T if it is daily! 
     216         IF( (sn_qsr%freqh  < 0.).OR.(sn_qsr%freqh  > 24.) ) & 
     217            & CALL ctl_stop( 'sbc_blk_init: Warm-layer param. (ln_skin_wl) not compatible with freq. of solar flux > daily' ) 
     218         IF( (sn_qsr%freqh == 24.).AND.(.NOT. ln_dm2dc) ) & 
     219            & CALL ctl_stop( 'sbc_blk_init: Please set ln_dm2dc=T for warm-layer param. (ln_skin_wl) to work properly' ) 
     220      END IF 
     221 
     222      ioptio = 0 
     223      IF( ln_humi_sph ) THEN 
     224         nhumi =  np_humi_sph    ;   ioptio = ioptio + 1 
     225      ENDIF 
     226      IF( ln_humi_dpt ) THEN 
     227         nhumi =  np_humi_dpt    ;   ioptio = ioptio + 1 
     228      ENDIF 
     229      IF( ln_humi_rlh ) THEN 
     230         nhumi =  np_humi_rlh    ;   ioptio = ioptio + 1 
     231      ENDIF 
     232      IF( ioptio /= 1 )   CALL ctl_stop( 'sbc_blk_init: Choose one and only one type of air humidity' ) 
    199233      ! 
    200234      IF( ln_dm2dc ) THEN                 !* check: diurnal cycle on Qsr 
    201235         IF( sn_qsr%freqh /= 24. )   CALL ctl_stop( 'sbc_blk_init: ln_dm2dc=T only with daily short-wave input' ) 
    202          IF( sn_qsr%ln_tint ) THEN  
     236         IF( sn_qsr%ln_tint ) THEN 
    203237            CALL ctl_warn( 'sbc_blk_init: ln_dm2dc=T daily qsr time interpolation done by sbcdcy module',   & 
    204238               &           '              ==> We force time interpolation = .false. for qsr' ) 
     
    208242      !                                   !* set the bulk structure 
    209243      !                                      !- store namelist information in an array 
     244      IF( ln_blk ) jpfld = 9 
     245      IF( ln_abl ) jpfld = 11 
     246      ALLOCATE( slf_i(jpfld) ) 
     247      ! 
    210248      slf_i(jp_wndi) = sn_wndi   ;   slf_i(jp_wndj) = sn_wndj 
    211249      slf_i(jp_qsr ) = sn_qsr    ;   slf_i(jp_qlw ) = sn_qlw 
    212250      slf_i(jp_tair) = sn_tair   ;   slf_i(jp_humi) = sn_humi 
    213251      slf_i(jp_prec) = sn_prec   ;   slf_i(jp_snow) = sn_snow 
    214       slf_i(jp_slp)  = sn_slp    ;   slf_i(jp_tdif) = sn_tdif 
    215       ! 
    216       lhftau = ln_taudif                     !- add an extra field if HF stress is used 
    217       jfld = jpfld - COUNT( (/.NOT.lhftau/) ) 
     252      slf_i(jp_slp ) = sn_slp 
     253      IF( ln_abl ) THEN 
     254         slf_i(jp_hpgi) = sn_hpgi   ;   slf_i(jp_hpgj) = sn_hpgj 
     255      END IF 
    218256      ! 
    219257      !                                      !- allocate the bulk structure 
    220       ALLOCATE( sf(jfld), STAT=ierror ) 
     258      ALLOCATE( sf(jpfld), STAT=ierror ) 
    221259      IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_init: unable to allocate sf structure' ) 
    222       DO ifpr= 1, jfld 
    223          ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) ) 
    224          IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) ) 
    225          IF( slf_i(ifpr)%freqh > 0. .AND. MOD( NINT(3600. * slf_i(ifpr)%freqh), nn_fsbc * NINT(rdt) ) /= 0 )   & 
    226             &  CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rdt*nn_fsbc is NOT a submultiple of atmospheric forcing frequency.',   & 
    227             &                 '               This is not ideal. You should consider changing either rdt or nn_fsbc value...' ) 
    228  
     260      ! 
     261      DO jfpr= 1, jpfld 
     262         ! 
     263         IF( TRIM(sf(jfpr)%clrootname) == 'NOT USED' ) THEN    !--  not used field  --!   (only now allocated and set to zero) 
     264            ALLOCATE( sf(jfpr)%fnow(jpi,jpj,1) ) 
     265            sf(jfpr)%fnow(:,:,1) = 0._wp 
     266         ELSE                                                  !-- used field  --! 
     267            IF(   ln_abl    .AND.                                                      & 
     268               &    ( jfpr == jp_wndi .OR. jfpr == jp_wndj .OR. jfpr == jp_humi .OR.   & 
     269               &      jfpr == jp_hpgi .OR. jfpr == jp_hpgj .OR. jfpr == jp_tair     )  ) THEN   ! ABL: some fields are 3D input 
     270               ALLOCATE( sf(jfpr)%fnow(jpi,jpj,jpka) ) 
     271               IF( slf_i(jfpr)%ln_tint )   ALLOCATE( sf(jfpr)%fdta(jpi,jpj,jpka,2) ) 
     272            ELSE                                                                                ! others or Bulk fields are 2D fiels 
     273               ALLOCATE( sf(jfpr)%fnow(jpi,jpj,1) ) 
     274               IF( slf_i(jfpr)%ln_tint )   ALLOCATE( sf(jfpr)%fdta(jpi,jpj,1,2) ) 
     275            ENDIF 
     276            ! 
     277            IF( slf_i(jfpr)%freqh > 0. .AND. MOD( NINT(3600. * slf_i(jfpr)%freqh), nn_fsbc * NINT(rdt) ) /= 0 )   & 
     278               &  CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rdt*nn_fsbc is NOT a submultiple of atmospheric forcing frequency.',   & 
     279               &                 '               This is not ideal. You should consider changing either rdt or nn_fsbc value...' ) 
     280         ENDIF 
    229281      END DO 
    230282      !                                      !- fill the bulk structure with namelist informations 
    231283      CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_init', 'surface boundary condition -- bulk formulae', 'namsbc_blk' ) 
    232284      ! 
    233       IF ( ln_wave ) THEN 
    234       !Activated wave module but neither drag nor stokes drift activated 
    235          IF ( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) )   THEN 
     285      IF( ln_wave ) THEN 
     286         !Activated wave module but neither drag nor stokes drift activated 
     287         IF( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) )   THEN 
    236288            CALL ctl_stop( 'STOP',  'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauwoc=F, ln_stcor=F' ) 
    237       !drag coefficient read from wave model definable only with mfs bulk formulae and core  
    238          ELSEIF (ln_cdgw .AND. .NOT. ln_NCAR )       THEN        
    239              CALL ctl_stop( 'drag coefficient read from wave model definable only with NCAR and CORE bulk formulae') 
    240          ELSEIF (ln_stcor .AND. .NOT. ln_sdw)                             THEN 
    241              CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T') 
     289            !drag coefficient read from wave model definable only with mfs bulk formulae and core 
     290         ELSEIF(ln_cdgw .AND. .NOT. ln_NCAR )       THEN 
     291            CALL ctl_stop( 'drag coefficient read from wave model definable only with NCAR and CORE bulk formulae') 
     292         ELSEIF(ln_stcor .AND. .NOT. ln_sdw)                             THEN 
     293            CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T') 
    242294         ENDIF 
    243295      ELSE 
    244       IF ( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor )                &  
    245          &   CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ',    & 
    246          &                  'with drag coefficient (ln_cdgw =T) '  ,                        & 
    247          &                  'or Stokes Drift (ln_sdw=T) ' ,                                 & 
    248          &                  'or ocean stress modification due to waves (ln_tauwoc=T) ',      &   
    249          &                  'or Stokes-Coriolis term (ln_stcori=T)'  ) 
    250       ENDIF  
    251       ! 
    252       !            
     296         IF( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor )                & 
     297            &   CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ',    & 
     298            &                  'with drag coefficient (ln_cdgw =T) '  ,                        & 
     299            &                  'or Stokes Drift (ln_sdw=T) ' ,                                 & 
     300            &                  'or ocean stress modification due to waves (ln_tauwoc=T) ',      & 
     301            &                  'or Stokes-Coriolis term (ln_stcori=T)'  ) 
     302      ENDIF 
     303      ! 
     304      IF( ln_abl ) THEN       ! ABL: read 3D fields for wind, temperature, humidity and pressure gradient 
     305         rn_zqt = ght_abl(2)          ! set the bulk altitude to ABL first level 
     306         rn_zu  = ght_abl(2) 
     307         IF(lwp) WRITE(numout,*) 
     308         IF(lwp) WRITE(numout,*) '   ABL formulation: overwrite rn_zqt & rn_zu with ABL first level altitude' 
     309      ENDIF 
     310      ! 
     311      ! set transfer coefficients to default sea-ice values 
     312      Cd_ice(:,:) = rCd_ice 
     313      Ch_ice(:,:) = rCd_ice 
     314      Ce_ice(:,:) = rCd_ice 
     315      ! 
    253316      IF(lwp) THEN                     !** Control print 
    254317         ! 
    255          WRITE(numout,*)                  !* namelist  
     318         WRITE(numout,*)                  !* namelist 
    256319         WRITE(numout,*) '   Namelist namsbc_blk (other than data information):' 
    257320         WRITE(numout,*) '      "NCAR"      algorithm   (Large and Yeager 2008)     ln_NCAR      = ', ln_NCAR 
    258321         WRITE(numout,*) '      "COARE 3.0" algorithm   (Fairall et al. 2003)       ln_COARE_3p0 = ', ln_COARE_3p0 
    259          WRITE(numout,*) '      "COARE 3.5" algorithm   (Edson et al. 2013)         ln_COARE_3p5 = ', ln_COARE_3p0 
    260          WRITE(numout,*) '      "ECMWF"     algorithm   (IFS cycle 31)              ln_ECMWF     = ', ln_ECMWF 
    261          WRITE(numout,*) '      add High freq.contribution to the stress module     ln_taudif    = ', ln_taudif 
     322         WRITE(numout,*) '      "COARE 3.6" algorithm (Fairall 2018 + Edson al 2013)ln_COARE_3p6 = ', ln_COARE_3p6 
     323         WRITE(numout,*) '      "ECMWF"     algorithm   (IFS cycle 45r1)            ln_ECMWF     = ', ln_ECMWF 
    262324         WRITE(numout,*) '      Air temperature and humidity reference height (m)   rn_zqt       = ', rn_zqt 
    263325         WRITE(numout,*) '      Wind vector reference height (m)                    rn_zu        = ', rn_zu 
     
    273335         CASE( np_NCAR      )   ;   WRITE(numout,*) '   ==>>>   "NCAR" algorithm        (Large and Yeager 2008)' 
    274336         CASE( np_COARE_3p0 )   ;   WRITE(numout,*) '   ==>>>   "COARE 3.0" algorithm   (Fairall et al. 2003)' 
    275          CASE( np_COARE_3p5 )   ;   WRITE(numout,*) '   ==>>>   "COARE 3.5" algorithm   (Edson et al. 2013)' 
    276          CASE( np_ECMWF     )   ;   WRITE(numout,*) '   ==>>>   "ECMWF" algorithm       (IFS cycle 31)' 
     337         CASE( np_COARE_3p6 )   ;   WRITE(numout,*) '   ==>>>   "COARE 3.6" algorithm (Fairall 2018+Edson et al. 2013)' 
     338         CASE( np_ECMWF     )   ;   WRITE(numout,*) '   ==>>>   "ECMWF" algorithm       (IFS cycle 45r1)' 
    277339         END SELECT 
    278340         ! 
     341         WRITE(numout,*) 
     342         WRITE(numout,*) '      use cool-skin  parameterization (SSST)  ln_skin_cs  = ', ln_skin_cs 
     343         WRITE(numout,*) '      use warm-layer parameterization (SSST)  ln_skin_wl  = ', ln_skin_wl 
     344         ! 
     345         WRITE(numout,*) 
     346         SELECT CASE( nhumi )              !* Print the choice of air humidity 
     347         CASE( np_humi_sph )   ;   WRITE(numout,*) '   ==>>>   air humidity is SPECIFIC HUMIDITY     [kg/kg]' 
     348         CASE( np_humi_dpt )   ;   WRITE(numout,*) '   ==>>>   air humidity is DEW-POINT TEMPERATURE [K]' 
     349         CASE( np_humi_rlh )   ;   WRITE(numout,*) '   ==>>>   air humidity is RELATIVE HUMIDITY     [%]' 
     350         END SELECT 
     351         ! 
    279352      ENDIF 
    280353      ! 
     
    289362      !!              (momentum, heat, freshwater and runoff) 
    290363      !! 
    291       !! ** Method  : (1) READ each fluxes in NetCDF files: 
    292       !!      the 10m wind velocity (i-component) (m/s)    at T-point 
    293       !!      the 10m wind velocity (j-component) (m/s)    at T-point 
    294       !!      the 10m or 2m specific humidity     ( % ) 
    295       !!      the solar heat                      (W/m2) 
    296       !!      the Long wave                       (W/m2) 
    297       !!      the 10m or 2m air temperature       (Kelvin) 
    298       !!      the total precipitation (rain+snow) (Kg/m2/s) 
    299       !!      the snow (solid prcipitation)       (kg/m2/s) 
    300       !!      the tau diff associated to HF tau   (N/m2)   at T-point   (ln_taudif=T) 
    301       !!              (2) CALL blk_oce 
     364      !! ** Method  : 
     365      !!              (1) READ each fluxes in NetCDF files: 
     366      !!      the wind velocity (i-component) at z=rn_zu  (m/s) at T-point 
     367      !!      the wind velocity (j-component) at z=rn_zu  (m/s) at T-point 
     368      !!      the specific humidity           at z=rn_zqt (kg/kg) 
     369      !!      the air temperature             at z=rn_zqt (Kelvin) 
     370      !!      the solar heat                              (W/m2) 
     371      !!      the Long wave                               (W/m2) 
     372      !!      the total precipitation (rain+snow)         (Kg/m2/s) 
     373      !!      the snow (solid precipitation)              (kg/m2/s) 
     374      !!      ABL dynamical forcing (i/j-components of either hpg or geostrophic winds) 
     375      !!              (2) CALL blk_oce_1 and blk_oce_2 
    302376      !! 
    303377      !!      C A U T I O N : never mask the surface stress fields 
     
    316390      !!---------------------------------------------------------------------- 
    317391      INTEGER, INTENT(in) ::   kt   ! ocean time step 
    318       !!--------------------------------------------------------------------- 
     392      !!---------------------------------------------------------------------- 
     393      REAL(wp), DIMENSION(jpi,jpj) ::   zssq, zcd_du, zsen, zevp 
     394      REAL(wp) :: ztmp 
     395      !!---------------------------------------------------------------------- 
    319396      ! 
    320397      CALL fld_read( kt, nn_fsbc, sf )             ! input fields provided at the current time-step 
    321       ! 
     398 
     399      ! Sanity/consistence test on humidity at first time step to detect potential screw-up: 
     400      IF( kt == nit000 ) THEN 
     401         WRITE(numout,*) '' 
     402#if defined key_agrif 
     403         WRITE(numout,*) ' === AGRIF => Sanity/consistence test on air humidity SKIPPED! :( ===' 
     404#else 
     405         ztmp = SUM(tmask(:,:,1)) ! number of ocean points on local proc domain 
     406         IF( ztmp > 8._wp ) THEN ! test only on proc domains with at least 8 ocean points! 
     407            ztmp = SUM(sf(jp_humi)%fnow(:,:,1)*tmask(:,:,1))/ztmp ! mean humidity over ocean on proc 
     408            SELECT CASE( nhumi ) 
     409            CASE( np_humi_sph ) ! specific humidity => expect: 0. <= something < 0.065 [kg/kg] (0.061 is saturation at 45degC !!!) 
     410               IF(  (ztmp < 0._wp) .OR. (ztmp > 0.065)  ) ztmp = -1._wp 
     411            CASE( np_humi_dpt ) ! dew-point temperature => expect: 110. <= something < 320. [K] 
     412               IF( (ztmp < 110._wp).OR.(ztmp > 320._wp) ) ztmp = -1._wp 
     413            CASE( np_humi_rlh ) ! relative humidity => expect: 0. <= something < 100. [%] 
     414               IF(  (ztmp < 0._wp) .OR.(ztmp > 100._wp) ) ztmp = -1._wp 
     415            END SELECT 
     416            IF(ztmp < 0._wp) THEN 
     417               WRITE(numout,'("   Mean humidity value found on proc #",i5.5," is: ",f)') narea, ztmp 
     418               CALL ctl_stop( 'STOP', 'Something is wrong with air humidity!!!', & 
     419                  &   ' ==> check the unit in your input files'       , & 
     420                  &   ' ==> check consistence of namelist choice: specific? relative? dew-point?', & 
     421                  &   ' ==> ln_humi_sph -> [kg/kg] | ln_humi_rlh -> [%] | ln_humi_dpt -> [K] !!!' ) 
     422            END IF 
     423         END IF 
     424         WRITE(numout,*) ' === Sanity/consistence test on air humidity sucessfuly passed! ===' 
     425#endif 
     426         WRITE(numout,*) '' 
     427      END IF !IF( kt == nit000 ) 
    322428      !                                            ! compute the surface ocean fluxes using bulk formulea 
    323       IF( MOD( kt - 1, nn_fsbc ) == 0 )   CALL blk_oce( kt, sf, sst_m, ssu_m, ssv_m ) 
    324  
     429      IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 
     430         CALL blk_oce_1( kt, sf(jp_wndi)%fnow(:,:,1), sf(jp_wndj)%fnow(:,:,1),   &   !   <<= in 
     431            &                sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1),   &   !   <<= in 
     432            &                sf(jp_slp )%fnow(:,:,1), sst_m, ssu_m, ssv_m,       &   !   <<= in 
     433            &                sf(jp_qsr )%fnow(:,:,1), sf(jp_qlw )%fnow(:,:,1),   &   !   <<= in (wl/cs) 
     434            &                zssq, zcd_du, zsen, zevp )                              !   =>> out 
     435 
     436         CALL blk_oce_2(     sf(jp_tair)%fnow(:,:,1), sf(jp_qsr )%fnow(:,:,1),   &   !   <<= in 
     437            &                sf(jp_qlw )%fnow(:,:,1), sf(jp_prec)%fnow(:,:,1),   &   !   <<= in 
     438            &                sf(jp_snow)%fnow(:,:,1), sst_m,                     &   !   <<= in 
     439            &                zsen, zevp )                                            !   <=> in out 
     440      ENDIF 
     441      ! 
    325442#if defined key_cice 
    326443      IF( MOD( kt - 1, nn_fsbc ) == 0 )   THEN 
    327444         qlw_ice(:,:,1)   = sf(jp_qlw )%fnow(:,:,1) 
    328          IF( ln_dm2dc ) THEN ; qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) 
    329          ELSE                ; qsr_ice(:,:,1) =          sf(jp_qsr)%fnow(:,:,1)  
    330          ENDIF  
     445         IF( ln_dm2dc ) THEN 
     446            qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) 
     447         ELSE 
     448            qsr_ice(:,:,1) =          sf(jp_qsr)%fnow(:,:,1) 
     449         ENDIF 
    331450         tatm_ice(:,:)    = sf(jp_tair)%fnow(:,:,1) 
    332          qatm_ice(:,:)    = sf(jp_humi)%fnow(:,:,1) 
     451 
     452         SELECT CASE( nhumi ) 
     453         CASE( np_humi_sph ) 
     454            qatm_ice(:,:) =           sf(jp_humi)%fnow(:,:,1) 
     455         CASE( np_humi_dpt ) 
     456            qatm_ice(:,:) = q_sat(    sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 
     457         CASE( np_humi_rlh ) 
     458            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 
     459         END SELECT 
     460 
    333461         tprecip(:,:)     = sf(jp_prec)%fnow(:,:,1) * rn_pfac 
    334462         sprecip(:,:)     = sf(jp_snow)%fnow(:,:,1) * rn_pfac 
     
    341469 
    342470 
    343    SUBROUTINE blk_oce( kt, sf, pst, pu, pv ) 
    344       !!--------------------------------------------------------------------- 
    345       !!                     ***  ROUTINE blk_oce  *** 
    346       !! 
    347       !! ** Purpose :   provide the momentum, heat and freshwater fluxes at 
    348       !!      the ocean surface at each time step 
    349       !! 
    350       !! ** Method  :   bulk formulea for the ocean using atmospheric 
    351       !!      fields read in sbc_read 
     471   SUBROUTINE blk_oce_1( kt, pwndi, pwndj , ptair, phumi, &  ! inp 
     472      &              pslp , pst   , pu   , pv,    &  ! inp 
     473      &              pqsr , pqlw  ,               &  ! inp 
     474      &              pssq , pcd_du, psen , pevp   )  ! out 
     475      !!--------------------------------------------------------------------- 
     476      !!                     ***  ROUTINE blk_oce_1  *** 
     477      !! 
     478      !! ** Purpose :   if ln_blk=T, computes surface momentum, heat and freshwater fluxes 
     479      !!                if ln_abl=T, computes Cd x |U|, Ch x |U|, Ce x |U| for ABL integration 
     480      !! 
     481      !! ** Method  :   bulk formulae using atmospheric fields from : 
     482      !!                if ln_blk=T, atmospheric fields read in sbc_read 
     483      !!                if ln_abl=T, the ABL model at previous time-step 
     484      !! 
     485      !! ** Outputs : - pssq    : surface humidity used to compute latent heat flux (kg/kg) 
     486      !!              - pcd_du  : Cd x |dU| at T-points  (m/s) 
     487      !!              - psen    : Ch x |dU| at T-points  (m/s) 
     488      !!              - pevp    : Ce x |dU| at T-points  (m/s) 
     489      !!--------------------------------------------------------------------- 
     490      INTEGER , INTENT(in   )                 ::   kt     ! time step index 
     491      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pwndi  ! atmospheric wind at U-point              [m/s] 
     492      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pwndj  ! atmospheric wind at V-point              [m/s] 
     493      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   phumi  ! specific humidity at T-points            [kg/kg] 
     494      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   ptair  ! potential temperature at T-points        [Kelvin] 
     495      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pslp   ! sea-level pressure                       [Pa] 
     496      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pst    ! surface temperature                      [Celcius] 
     497      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pu     ! surface current at U-point (i-component) [m/s] 
     498      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pv     ! surface current at V-point (j-component) [m/s] 
     499      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pqsr   ! 
     500      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pqlw   ! 
     501      REAL(wp), INTENT(  out), DIMENSION(:,:) ::   pssq   ! specific humidity at pst                 [kg/kg] 
     502      REAL(wp), INTENT(  out), DIMENSION(:,:) ::   pcd_du ! Cd x |dU| at T-points                    [m/s] 
     503      REAL(wp), INTENT(  out), DIMENSION(:,:) ::   psen   ! Ch x |dU| at T-points                    [m/s] 
     504      REAL(wp), INTENT(  out), DIMENSION(:,:) ::   pevp   ! Ce x |dU| at T-points                    [m/s] 
     505      ! 
     506      INTEGER  ::   ji, jj               ! dummy loop indices 
     507      REAL(wp) ::   zztmp                ! local variable 
     508      REAL(wp), DIMENSION(jpi,jpj) ::   zwnd_i, zwnd_j    ! wind speed components at T-point 
     509      REAL(wp), DIMENSION(jpi,jpj) ::   zst               ! surface temperature in Kelvin 
     510      REAL(wp), DIMENSION(jpi,jpj) ::   zU_zu             ! bulk wind speed at height zu  [m/s] 
     511      REAL(wp), DIMENSION(jpi,jpj) ::   ztpot             ! potential temperature of air at z=rn_zqt [K] 
     512      REAL(wp), DIMENSION(jpi,jpj) ::   zqair             ! specific humidity     of air at z=rn_zqt [kg/kg] 
     513      REAL(wp), DIMENSION(jpi,jpj) ::   zcd_oce           ! momentum transfert coefficient over ocean 
     514      REAL(wp), DIMENSION(jpi,jpj) ::   zch_oce           ! sensible heat transfert coefficient over ocean 
     515      REAL(wp), DIMENSION(jpi,jpj) ::   zce_oce           ! latent   heat transfert coefficient over ocean 
     516      REAL(wp), DIMENSION(jpi,jpj) ::   zqla              ! latent heat flux 
     517      REAL(wp), DIMENSION(jpi,jpj) ::   zztmp1, zztmp2 
     518      !!--------------------------------------------------------------------- 
     519      ! 
     520      ! local scalars ( place there for vector optimisation purposes) 
     521      zst(:,:) = pst(:,:) + rt0      ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K) 
     522 
     523      ! ----------------------------------------------------------------------------- ! 
     524      !      0   Wind components and module at T-point relative to the moving ocean   ! 
     525      ! ----------------------------------------------------------------------------- ! 
     526 
     527      ! ... components ( U10m - U_oce ) at T-point (unmasked) 
     528#if defined key_cyclone 
     529      zwnd_i(:,:) = 0._wp 
     530      zwnd_j(:,:) = 0._wp 
     531      CALL wnd_cyc( kt, zwnd_i, zwnd_j )    ! add analytical tropical cyclone (Vincent et al. JGR 2012) 
     532      DO jj = 2, jpjm1 
     533         DO ji = fs_2, fs_jpim1   ! vect. opt. 
     534            pwndi(ji,jj) = pwndi(ji,jj) + zwnd_i(ji,jj) 
     535            pwndj(ji,jj) = pwndj(ji,jj) + zwnd_j(ji,jj) 
     536         END DO 
     537      END DO 
     538#endif 
     539      DO jj = 2, jpjm1 
     540         DO ji = fs_2, fs_jpim1   ! vect. opt. 
     541            zwnd_i(ji,jj) = (  pwndi(ji,jj) - rn_vfac * 0.5 * ( pu(ji-1,jj  ) + pu(ji,jj) )  ) 
     542            zwnd_j(ji,jj) = (  pwndj(ji,jj) - rn_vfac * 0.5 * ( pv(ji  ,jj-1) + pv(ji,jj) )  ) 
     543         END DO 
     544      END DO 
     545      CALL lbc_lnk_multi( 'sbcblk', zwnd_i, 'T', -1., zwnd_j, 'T', -1. ) 
     546      ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked) 
     547      wndm(:,:) = SQRT(  zwnd_i(:,:) * zwnd_i(:,:)   & 
     548         &             + zwnd_j(:,:) * zwnd_j(:,:)  ) * tmask(:,:,1) 
     549 
     550      ! ----------------------------------------------------------------------------- ! 
     551      !      I   Solar FLUX                                                           ! 
     552      ! ----------------------------------------------------------------------------- ! 
     553 
     554      ! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle                    ! Short Wave 
     555      zztmp = 1. - albo 
     556      IF( ln_dm2dc ) THEN 
     557         qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1) 
     558      ELSE 
     559         qsr(:,:) = zztmp *          sf(jp_qsr)%fnow(:,:,1)   * tmask(:,:,1) 
     560      ENDIF 
     561 
     562 
     563      ! ----------------------------------------------------------------------------- ! 
     564      !     II   Turbulent FLUXES                                                     ! 
     565      ! ----------------------------------------------------------------------------- ! 
     566 
     567      ! specific humidity at SST 
     568      pssq(:,:) = rdct_qsat_salt * q_sat( zst(:,:), pslp(:,:) ) 
     569 
     570      IF( ln_skin_cs .OR. ln_skin_wl ) THEN 
     571         zztmp1(:,:) = zst(:,:) 
     572         zztmp2(:,:) = pssq(:,:) 
     573      ENDIF 
     574 
     575      ! specific humidity of air at "rn_zqt" m above the sea 
     576      SELECT CASE( nhumi ) 
     577      CASE( np_humi_sph ) 
     578         zqair(:,:) = phumi(:,:)      ! what we read in file is already a spec. humidity! 
     579      CASE( np_humi_dpt ) 
     580         !IF(lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of d_air and slp !' !LBrm 
     581         zqair(:,:) = q_sat( phumi(:,:), pslp(:,:) ) 
     582      CASE( np_humi_rlh ) 
     583         !IF(lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of RH, t_air and slp !' !LBrm 
     584         zqair(:,:) = q_air_rh( 0.01_wp*phumi(:,:), ptair(:,:), pslp(:,:) ) !LB: 0.01 => RH is % percent in file 
     585      END SELECT 
     586      ! 
     587      ! potential temperature of air at "rn_zqt" m above the sea 
     588      IF( ln_abl ) THEN 
     589         ztpot = ptair(:,:) 
     590      ELSE 
     591         ! Estimate of potential temperature at z=rn_zqt, based on adiabatic lapse-rate 
     592         !    (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2 
     593         !    (since reanalysis products provide T at z, not theta !) 
     594         !#LB: because AGRIF hates functions that return something else than a scalar, need to 
     595         !     use scalar version of gamma_moist() ... 
     596         DO jj = 1, jpj 
     597            DO ji = 1, jpi 
     598               ztpot(ji,jj) = ptair(ji,jj) + gamma_moist( ptair(ji,jj), zqair(ji,jj) ) * rn_zqt 
     599            END DO 
     600         END DO 
     601      ENDIF 
     602 
     603 
     604 
     605      !! Time to call the user-selected bulk parameterization for 
     606      !!  ==  transfer coefficients  ==!   Cd, Ch, Ce at T-point, and more... 
     607      SELECT CASE( nblk ) 
     608 
     609      CASE( np_NCAR      ) 
     610         CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm,                              & 
     611            &                zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
     612 
     613      CASE( np_COARE_3p0 ) 
     614         CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 
     615            &                zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,   & 
     616            &                Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) ) 
     617 
     618      CASE( np_COARE_3p6 ) 
     619         CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 
     620            &                zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,   & 
     621            &                Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) ) 
     622 
     623      CASE( np_ECMWF     ) 
     624         CALL turb_ecmwf   ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl,  & 
     625            &                zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,   & 
     626            &                Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) ) 
     627 
     628      CASE DEFAULT 
     629         CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' ) 
     630 
     631      END SELECT 
     632 
     633      IF( ln_skin_cs .OR. ln_skin_wl ) THEN 
     634         !! In the presence of sea-ice we forget about the cool-skin/warm-layer update of zst and pssq: 
     635         WHERE ( fr_i < 0.001_wp ) 
     636            ! zst and pssq have been updated by cool-skin/warm-layer scheme and we keep it!!! 
     637            zst(:,:)  =  zst(:,:)*tmask(:,:,1) 
     638            pssq(:,:) = pssq(:,:)*tmask(:,:,1) 
     639         ELSEWHERE 
     640            ! we forget about the update... 
     641            zst(:,:)  = zztmp1(:,:) !#LB: using what we backed up before skin-algo 
     642            pssq(:,:) = zztmp2(:,:) !#LB:  "   "   " 
     643         END WHERE 
     644      END IF 
     645 
     646      !!      CALL iom_put( "Cd_oce", zcd_oce)  ! output value of pure ocean-atm. transfer coef. 
     647      !!      CALL iom_put( "Ch_oce", zch_oce)  ! output value of pure ocean-atm. transfer coef. 
     648 
     649      IF( ABS(rn_zu - rn_zqt) < 0.1_wp ) THEN 
     650         !! If zu == zt, then ensuring once for all that: 
     651         t_zu(:,:) = ztpot(:,:) 
     652         q_zu(:,:) = zqair(:,:) 
     653      ENDIF 
     654 
     655 
     656      !  Turbulent fluxes over ocean  => BULK_FORMULA @ sbcblk_phy.F90 
     657      ! ------------------------------------------------------------- 
     658 
     659      IF( ln_abl ) THEN         !==  ABL formulation  ==!   multiplication by rho_air and turbulent fluxes computation done in ablstp 
     660         !! FL do we need this multiplication by tmask ... ??? 
     661         DO jj = 1, jpj 
     662            DO ji = 1, jpi 
     663               zztmp = zU_zu(ji,jj) !* tmask(ji,jj,1) 
     664               wndm(ji,jj)   = zztmp                   ! Store zU_zu in wndm to compute ustar2 in ablmod 
     665               pcd_du(ji,jj) = zztmp * zcd_oce(ji,jj) 
     666               psen(ji,jj)   = zztmp * zch_oce(ji,jj) 
     667               pevp(ji,jj)   = zztmp * zce_oce(ji,jj) 
     668            END DO 
     669         END DO 
     670      ELSE                      !==  BLK formulation  ==!   turbulent fluxes computation 
     671         CALL BULK_FORMULA( rn_zu, zst(:,:), pssq(:,:), t_zu(:,:), q_zu(:,:), & 
     672            &               zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:),         & 
     673            &               wndm(:,:), zU_zu(:,:), pslp(:,:),                 & 
     674            &               taum(:,:), psen(:,:), zqla(:,:),                  & 
     675            &               pEvap=pevp(:,:), prhoa=rhoa(:,:) ) 
     676 
     677         zqla(:,:) = zqla(:,:) * tmask(:,:,1) 
     678         psen(:,:) = psen(:,:) * tmask(:,:,1) 
     679         taum(:,:) = taum(:,:) * tmask(:,:,1) 
     680         pevp(:,:) = pevp(:,:) * tmask(:,:,1) 
     681 
     682         ! Tau i and j component on T-grid points, using array "zcd_oce" as a temporary array... 
     683         zcd_oce = 0._wp 
     684         WHERE ( wndm > 0._wp ) zcd_oce = taum / wndm 
     685         zwnd_i = zcd_oce * zwnd_i 
     686         zwnd_j = zcd_oce * zwnd_j 
     687 
     688         CALL iom_put( "taum_oce", taum )   ! output wind stress module 
     689 
     690         ! ... utau, vtau at U- and V_points, resp. 
     691         !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines 
     692         !     Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves 
     693         DO jj = 1, jpjm1 
     694            DO ji = 1, fs_jpim1 
     695               utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj  ) ) & 
     696                  &          * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1)) 
     697               vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji  ,jj+1) ) & 
     698                  &          * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1)) 
     699            END DO 
     700         END DO 
     701         CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. ) 
     702 
     703         IF(ln_ctl) THEN 
     704            CALL prt_ctl( tab2d_1=wndm  , clinfo1=' blk_oce_1: wndm   : ') 
     705            CALL prt_ctl( tab2d_1=utau  , clinfo1=' blk_oce_1: utau   : ', mask1=umask,   & 
     706               &          tab2d_2=vtau  , clinfo2='            vtau   : ', mask2=vmask ) 
     707         ENDIF 
     708         ! 
     709      ENDIF 
     710      ! 
     711      IF(ln_ctl) THEN 
     712         CALL prt_ctl( tab2d_1=pevp  , clinfo1=' blk_oce_1: pevp   : ' ) 
     713         CALL prt_ctl( tab2d_1=psen  , clinfo1=' blk_oce_1: psen   : ' ) 
     714         CALL prt_ctl( tab2d_1=pssq  , clinfo1=' blk_oce_1: pssq   : ' ) 
     715      ENDIF 
     716      ! 
     717   END SUBROUTINE blk_oce_1 
     718 
     719 
     720   SUBROUTINE blk_oce_2( ptair, pqsr, pqlw, pprec,   &   ! <<= in 
     721      &          psnow, pst , psen, pevp     )   ! <<= in 
     722      !!--------------------------------------------------------------------- 
     723      !!                     ***  ROUTINE blk_oce_2  *** 
     724      !! 
     725      !! ** Purpose :   finalize the momentum, heat and freshwater fluxes computation 
     726      !!                at the ocean surface at each time step knowing Cd, Ch, Ce and 
     727      !!                atmospheric variables (from ABL or external data) 
    352728      !! 
    353729      !! ** Outputs : - utau    : i-component of the stress at U-point  (N/m2) 
     
    358734      !!              - qns     : Non Solar heat flux over the ocean    (W/m2) 
    359735      !!              - emp     : evaporation minus precipitation       (kg/m2/s) 
    360       !! 
    361       !!  ** Nota  :   sf has to be a dummy argument for AGRIF on NEC 
    362       !!--------------------------------------------------------------------- 
    363       INTEGER  , INTENT(in   )                 ::   kt    ! time step index 
    364       TYPE(fld), INTENT(inout), DIMENSION(:)   ::   sf    ! input data 
    365       REAL(wp) , INTENT(in)   , DIMENSION(:,:) ::   pst   ! surface temperature                      [Celcius] 
    366       REAL(wp) , INTENT(in)   , DIMENSION(:,:) ::   pu    ! surface current at U-point (i-component) [m/s] 
    367       REAL(wp) , INTENT(in)   , DIMENSION(:,:) ::   pv    ! surface current at V-point (j-component) [m/s] 
     736      !!--------------------------------------------------------------------- 
     737      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptair 
     738      REAL(wp), INTENT(in), DIMENSION(:,:) ::   pqsr 
     739      REAL(wp), INTENT(in), DIMENSION(:,:) ::   pqlw 
     740      REAL(wp), INTENT(in), DIMENSION(:,:) ::   pprec 
     741      REAL(wp), INTENT(in), DIMENSION(:,:) ::   psnow 
     742      REAL(wp), INTENT(in), DIMENSION(:,:) ::   pst   ! surface temperature                      [Celcius] 
     743      REAL(wp), INTENT(in), DIMENSION(:,:) ::   psen 
     744      REAL(wp), INTENT(in), DIMENSION(:,:) ::   pevp 
    368745      ! 
    369746      INTEGER  ::   ji, jj               ! dummy loop indices 
    370       REAL(wp) ::   zztmp                ! local variable 
    371       REAL(wp), DIMENSION(jpi,jpj) ::   zwnd_i, zwnd_j    ! wind speed components at T-point 
    372       REAL(wp), DIMENSION(jpi,jpj) ::   zsq               ! specific humidity at pst 
    373       REAL(wp), DIMENSION(jpi,jpj) ::   zqlw, zqsb        ! long wave and sensible heat fluxes 
    374       REAL(wp), DIMENSION(jpi,jpj) ::   zqla, zevap       ! latent heat fluxes and evaporation 
     747      REAL(wp) ::   zztmp,zz1,zz2,zz3    ! local variable 
     748      REAL(wp), DIMENSION(jpi,jpj) ::   zqlw              ! long wave and sensible heat fluxes 
     749      REAL(wp), DIMENSION(jpi,jpj) ::   zqla              ! latent heat fluxes and evaporation 
    375750      REAL(wp), DIMENSION(jpi,jpj) ::   zst               ! surface temperature in Kelvin 
    376       REAL(wp), DIMENSION(jpi,jpj) ::   zU_zu             ! bulk wind speed at height zu  [m/s] 
    377       REAL(wp), DIMENSION(jpi,jpj) ::   ztpot             ! potential temperature of air at z=rn_zqt [K] 
    378       REAL(wp), DIMENSION(jpi,jpj) ::   zrhoa             ! density of air   [kg/m^3] 
    379751      !!--------------------------------------------------------------------- 
    380752      ! 
     
    382754      zst(:,:) = pst(:,:) + rt0      ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K) 
    383755 
     756 
    384757      ! ----------------------------------------------------------------------------- ! 
    385       !      0   Wind components and module at T-point relative to the moving ocean   ! 
     758      !     III    Net longwave radiative FLUX                                        ! 
    386759      ! ----------------------------------------------------------------------------- ! 
    387760 
    388       ! ... components ( U10m - U_oce ) at T-point (unmasked) 
    389 !!gm    move zwnd_i (_j) set to zero  inside the key_cyclone ??? 
    390       zwnd_i(:,:) = 0._wp 
    391       zwnd_j(:,:) = 0._wp 
    392 #if defined key_cyclone 
    393       CALL wnd_cyc( kt, zwnd_i, zwnd_j )    ! add analytical tropical cyclone (Vincent et al. JGR 2012) 
    394       DO jj = 2, jpjm1 
    395          DO ji = fs_2, fs_jpim1   ! vect. opt. 
    396             sf(jp_wndi)%fnow(ji,jj,1) = sf(jp_wndi)%fnow(ji,jj,1) + zwnd_i(ji,jj) 
    397             sf(jp_wndj)%fnow(ji,jj,1) = sf(jp_wndj)%fnow(ji,jj,1) + zwnd_j(ji,jj) 
    398          END DO 
    399       END DO 
    400 #endif 
    401       DO jj = 2, jpjm1 
    402          DO ji = fs_2, fs_jpim1   ! vect. opt. 
    403             zwnd_i(ji,jj) = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pu(ji-1,jj  ) + pu(ji,jj) )  ) 
    404             zwnd_j(ji,jj) = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pv(ji  ,jj-1) + pv(ji,jj) )  ) 
    405          END DO 
    406       END DO 
    407       CALL lbc_lnk_multi( 'sbcblk', zwnd_i, 'T', -1., zwnd_j, 'T', -1. ) 
    408       ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked) 
    409       wndm(:,:) = SQRT(  zwnd_i(:,:) * zwnd_i(:,:)   & 
    410          &             + zwnd_j(:,:) * zwnd_j(:,:)  ) * tmask(:,:,1) 
    411  
    412       ! ----------------------------------------------------------------------------- ! 
    413       !      I   Radiative FLUXES                                                     ! 
    414       ! ----------------------------------------------------------------------------- ! 
    415  
    416       ! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle                    ! Short Wave 
    417       zztmp = 1. - albo 
    418       IF( ln_dm2dc ) THEN   ;   qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1) 
    419       ELSE                  ;   qsr(:,:) = zztmp *          sf(jp_qsr)%fnow(:,:,1)   * tmask(:,:,1) 
    420       ENDIF 
    421  
    422       zqlw(:,:) = (  sf(jp_qlw)%fnow(:,:,1) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:)  ) * tmask(:,:,1)   ! Long  Wave 
    423  
    424       ! ----------------------------------------------------------------------------- ! 
    425       !     II    Turbulent FLUXES                                                    ! 
    426       ! ----------------------------------------------------------------------------- ! 
    427  
    428       ! ... specific humidity at SST and IST tmask( 
    429       zsq(:,:) = 0.98 * q_sat( zst(:,:), sf(jp_slp)%fnow(:,:,1) ) 
    430       !! 
    431       !! Estimate of potential temperature at z=rn_zqt, based on adiabatic lapse-rate 
    432       !!    (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2 
    433       !!    (since reanalysis products provide T at z, not theta !) 
    434       ztpot = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1) ) * rn_zqt 
    435  
    436       SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point 
    437       ! 
    438       CASE( np_NCAR      )   ;   CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! NCAR-COREv2 
    439          &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    440       CASE( np_COARE_3p0 )   ;   CALL turb_coare   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! COARE v3.0 
    441          &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    442       CASE( np_COARE_3p5 )   ;   CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! COARE v3.5 
    443          &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    444       CASE( np_ECMWF     )   ;   CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! ECMWF 
    445          &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    446       CASE DEFAULT 
    447          CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' ) 
    448       END SELECT 
    449  
    450       !                          ! Compute true air density : 
    451       IF( ABS(rn_zu - rn_zqt) > 0.01 ) THEN     ! At zu: (probably useless to remove zrho*grav*rn_zu from SLP...) 
    452          zrhoa(:,:) = rho_air( t_zu(:,:)              , q_zu(:,:)              , sf(jp_slp)%fnow(:,:,1) ) 
    453       ELSE                                      ! At zt: 
    454          zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 
    455       END IF 
    456  
    457 !!      CALL iom_put( "Cd_oce", Cd_atm)  ! output value of pure ocean-atm. transfer coef. 
    458 !!      CALL iom_put( "Ch_oce", Ch_atm)  ! output value of pure ocean-atm. transfer coef. 
    459  
    460       DO jj = 1, jpj             ! tau module, i and j component 
    461          DO ji = 1, jpi 
    462             zztmp = zrhoa(ji,jj)  * zU_zu(ji,jj) * Cd_atm(ji,jj)   ! using bulk wind speed 
    463             taum  (ji,jj) = zztmp * wndm  (ji,jj) 
    464             zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) 
    465             zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj) 
    466          END DO 
    467       END DO 
    468  
    469       !                          ! add the HF tau contribution to the wind stress module 
    470       IF( lhftau )   taum(:,:) = taum(:,:) + sf(jp_tdif)%fnow(:,:,1) 
    471  
    472       CALL iom_put( "taum_oce", taum )   ! output wind stress module 
    473  
    474       ! ... utau, vtau at U- and V_points, resp. 
    475       !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines 
    476       !     Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves 
    477       DO jj = 1, jpjm1 
    478          DO ji = 1, fs_jpim1 
    479             utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj  ) ) & 
    480                &          * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1)) 
    481             vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji  ,jj+1) ) & 
    482                &          * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1)) 
    483          END DO 
    484       END DO 
    485       CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. ) 
     761      !! LB: now moved after Turbulent fluxes because must use the skin temperature rather that the SST 
     762      !! (zst is skin temperature if ln_skin_cs==.TRUE. .OR. ln_skin_wl==.TRUE.) 
     763      zqlw(:,:) = emiss_w * ( pqlw(:,:) - stefan*zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:) ) * tmask(:,:,1)   ! Net radiative longwave flux 
    486764 
    487765      !  Turbulent fluxes over ocean 
    488766      ! ----------------------------- 
    489767 
    490       ! zqla used as temporary array, for rho*U (common term of bulk formulae): 
    491       zqla(:,:) = zrhoa(:,:) * zU_zu(:,:) * tmask(:,:,1) 
    492  
    493       IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN 
    494          !! q_air and t_air are given at 10m (wind reference height) 
    495          zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - sf(jp_humi)%fnow(:,:,1)) ) ! Evaporation, using bulk wind speed 
    496          zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - ztpot(:,:)             )   ! Sensible Heat, using bulk wind speed 
    497       ELSE 
    498          !! q_air and t_air are not given at 10m (wind reference height) 
    499          ! Values of temp. and hum. adjusted to height of wind during bulk algorithm iteration must be used!!! 
    500          zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - q_zu(:,:) ) ) ! Evaporation, using bulk wind speed 
    501          zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - t_zu(:,:) )   ! Sensible Heat, using bulk wind speed 
    502       ENDIF 
    503  
    504       zqla(:,:) = L_vap(zst(:,:)) * zevap(:,:)     ! Latent Heat flux 
    505  
     768      ! use scalar version of L_vap() for AGRIF compatibility 
     769      DO jj = 1, jpj 
     770         DO ji = 1, jpi 
     771            zqla(ji,jj) = -1._wp * L_vap( zst(ji,jj) ) * pevp(ji,jj)    ! Latent Heat flux !!GS: possibility to add a global qla to avoid recomputation after abl update 
     772         ENDDO 
     773      ENDDO 
    506774 
    507775      IF(ln_ctl) THEN 
    508          CALL prt_ctl( tab2d_1=zqla  , clinfo1=' blk_oce: zqla   : ', tab2d_2=Ce_atm , clinfo2=' Ce_oce  : ' ) 
    509          CALL prt_ctl( tab2d_1=zqsb  , clinfo1=' blk_oce: zqsb   : ', tab2d_2=Ch_atm , clinfo2=' Ch_oce  : ' ) 
    510          CALL prt_ctl( tab2d_1=zqlw  , clinfo1=' blk_oce: zqlw   : ', tab2d_2=qsr, clinfo2=' qsr : ' ) 
    511          CALL prt_ctl( tab2d_1=zsq   , clinfo1=' blk_oce: zsq    : ', tab2d_2=zst, clinfo2=' zst : ' ) 
    512          CALL prt_ctl( tab2d_1=utau  , clinfo1=' blk_oce: utau   : ', mask1=umask,   & 
    513             &          tab2d_2=vtau  , clinfo2=           ' vtau : ', mask2=vmask ) 
    514          CALL prt_ctl( tab2d_1=wndm  , clinfo1=' blk_oce: wndm   : ') 
    515          CALL prt_ctl( tab2d_1=zst   , clinfo1=' blk_oce: zst    : ') 
     776         CALL prt_ctl( tab2d_1=zqla  , clinfo1=' blk_oce_2: zqla   : ' ) 
     777         CALL prt_ctl( tab2d_1=zqlw  , clinfo1=' blk_oce_2: zqlw   : ', tab2d_2=qsr, clinfo2=' qsr : ' ) 
     778 
    516779      ENDIF 
    517780 
    518781      ! ----------------------------------------------------------------------------- ! 
    519       !     III    Total FLUXES                                                       ! 
     782      !     IV    Total FLUXES                                                       ! 
    520783      ! ----------------------------------------------------------------------------- ! 
    521784      ! 
    522       emp (:,:) = (  zevap(:,:)                                          &   ! mass flux (evap. - precip.) 
    523          &         - sf(jp_prec)%fnow(:,:,1) * rn_pfac  ) * tmask(:,:,1) 
    524       ! 
    525       qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar 
    526          &     - sf(jp_snow)%fnow(:,:,1) * rn_pfac * rLfus                        &   ! remove latent melting heat for solid precip 
    527          &     - zevap(:,:) * pst(:,:) * rcp                                      &   ! remove evap heat content at SST 
    528          &     + ( sf(jp_prec)%fnow(:,:,1) - sf(jp_snow)%fnow(:,:,1) ) * rn_pfac  &   ! add liquid precip heat content at Tair 
    529          &     * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp                          & 
    530          &     + sf(jp_snow)%fnow(:,:,1) * rn_pfac                                &   ! add solid  precip heat content at min(Tair,Tsnow) 
    531          &     * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi 
     785      emp (:,:) = (  pevp(:,:)                                       &   ! mass flux (evap. - precip.) 
     786         &         - pprec(:,:) * rn_pfac  ) * tmask(:,:,1) 
     787      ! 
     788      qns(:,:) = zqlw(:,:) + psen(:,:) + zqla(:,:)                   &   ! Downward Non Solar 
     789         &     - psnow(:,:) * rn_pfac * rLfus                        &   ! remove latent melting heat for solid precip 
     790         &     - pevp(:,:) * pst(:,:) * rcp                          &   ! remove evap heat content at SST !LB??? pst is Celsius !? 
     791         &     + ( pprec(:,:) - psnow(:,:) ) * rn_pfac               &   ! add liquid precip heat content at Tair 
     792         &     * ( ptair(:,:) - rt0 ) * rcp                          & 
     793         &     + psnow(:,:) * rn_pfac                                &   ! add solid  precip heat content at min(Tair,Tsnow) 
     794         &     * ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi 
    532795      qns(:,:) = qns(:,:) * tmask(:,:,1) 
    533796      ! 
    534797#if defined key_si3 
    535       qns_oce(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                ! non solar without emp (only needed by SI3) 
     798      qns_oce(:,:) = zqlw(:,:) + psen(:,:) + zqla(:,:)                             ! non solar without emp (only needed by SI3) 
    536799      qsr_oce(:,:) = qsr(:,:) 
    537800#endif 
    538801      ! 
     802      CALL iom_put( "rho_air"  , rhoa*tmask(:,:,1) )       ! output air density [kg/m^3] 
     803      CALL iom_put( "evap_oce" , pevp )                    ! evaporation 
     804      CALL iom_put( "qlw_oce"  , zqlw )                    ! output downward longwave heat over the ocean 
     805      CALL iom_put( "qsb_oce"  , psen )                    ! output downward sensible heat over the ocean 
     806      CALL iom_put( "qla_oce"  , zqla )                    ! output downward latent   heat over the ocean 
     807      tprecip(:,:) = pprec(:,:) * rn_pfac * tmask(:,:,1)   ! output total precipitation [kg/m2/s] 
     808      sprecip(:,:) = psnow(:,:) * rn_pfac * tmask(:,:,1)   ! output solid precipitation [kg/m2/s] 
     809      CALL iom_put( 'snowpre', sprecip )                   ! Snow 
     810      CALL iom_put( 'precip' , tprecip )                   ! Total precipitation 
     811      ! 
    539812      IF ( nn_ice == 0 ) THEN 
    540          CALL iom_put( "qlw_oce" ,   zqlw )                 ! output downward longwave heat over the ocean 
    541          CALL iom_put( "qsb_oce" , - zqsb )                 ! output downward sensible heat over the ocean 
    542          CALL iom_put( "qla_oce" , - zqla )                 ! output downward latent   heat over the ocean 
    543          CALL iom_put( "qemp_oce",   qns-zqlw+zqsb+zqla )   ! output downward heat content of E-P over the ocean 
    544          CALL iom_put( "qns_oce" ,   qns  )                 ! output downward non solar heat over the ocean 
    545          CALL iom_put( "qsr_oce" ,   qsr  )                 ! output downward solar heat over the ocean 
    546          CALL iom_put( "qt_oce"  ,   qns+qsr )              ! output total downward heat over the ocean 
    547          tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s] 
    548          sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s] 
    549          CALL iom_put( 'snowpre', sprecip )                 ! Snow 
    550          CALL iom_put( 'precip' , tprecip )                 ! Total precipitation 
     813         CALL iom_put( "qemp_oce" , qns-zqlw-psen-zqla )   ! output downward heat content of E-P over the ocean 
     814         CALL iom_put( "qns_oce"  ,   qns  )               ! output downward non solar heat over the ocean 
     815         CALL iom_put( "qsr_oce"  ,   qsr  )               ! output downward solar heat over the ocean 
     816         CALL iom_put( "qt_oce"   ,   qns+qsr )            ! output total downward heat over the ocean 
     817      ENDIF 
     818      ! 
     819      IF( ln_skin_cs .OR. ln_skin_wl ) THEN 
     820         CALL iom_put( "t_skin" ,  (zst - rt0) * tmask(:,:,1) )           ! T_skin in Celsius 
     821         CALL iom_put( "dt_skin" , (zst - pst - rt0) * tmask(:,:,1) )     ! T_skin - SST temperature difference... 
    551822      ENDIF 
    552823      ! 
    553824      IF(ln_ctl) THEN 
    554          CALL prt_ctl(tab2d_1=zqsb , clinfo1=' blk_oce: zqsb   : ', tab2d_2=zqlw , clinfo2=' zqlw  : ') 
    555          CALL prt_ctl(tab2d_1=zqla , clinfo1=' blk_oce: zqla   : ', tab2d_2=qsr  , clinfo2=' qsr   : ') 
    556          CALL prt_ctl(tab2d_1=pst  , clinfo1=' blk_oce: pst    : ', tab2d_2=emp  , clinfo2=' emp   : ') 
    557          CALL prt_ctl(tab2d_1=utau , clinfo1=' blk_oce: utau   : ', mask1=umask,   & 
    558             &         tab2d_2=vtau , clinfo2=              ' vtau  : ' , mask2=vmask ) 
    559       ENDIF 
    560       ! 
    561    END SUBROUTINE blk_oce 
    562  
    563  
    564  
    565    FUNCTION rho_air( ptak, pqa, pslp ) 
    566       !!------------------------------------------------------------------------------- 
    567       !!                           ***  FUNCTION rho_air  *** 
    568       !! 
    569       !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere 
    570       !! 
    571       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)  
    572       !!------------------------------------------------------------------------------- 
    573       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak      ! air temperature             [K] 
    574       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa       ! air specific humidity   [kg/kg] 
    575       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp      ! pressure in                [Pa] 
    576       REAL(wp), DIMENSION(jpi,jpj)             ::   rho_air   ! density of moist air   [kg/m^3] 
    577       !!------------------------------------------------------------------------------- 
    578       ! 
    579       rho_air = pslp / (  R_dry*ptak * ( 1._wp + rctv0*pqa )  ) 
    580       ! 
    581    END FUNCTION rho_air 
    582  
    583  
    584    FUNCTION cp_air( pqa ) 
    585       !!------------------------------------------------------------------------------- 
    586       !!                           ***  FUNCTION cp_air  *** 
    587       !! 
    588       !! ** Purpose : Compute specific heat (Cp) of moist air 
    589       !! 
    590       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    591       !!------------------------------------------------------------------------------- 
    592       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa      ! air specific humidity         [kg/kg] 
    593       REAL(wp), DIMENSION(jpi,jpj)             ::   cp_air   ! specific heat of moist air   [J/K/kg] 
    594       !!------------------------------------------------------------------------------- 
    595       ! 
    596       Cp_air = Cp_dry + Cp_vap * pqa 
    597       ! 
    598    END FUNCTION cp_air 
    599  
    600  
    601    FUNCTION q_sat( ptak, pslp ) 
    602       !!---------------------------------------------------------------------------------- 
    603       !!                           ***  FUNCTION q_sat  *** 
    604       !! 
    605       !! ** Purpose : Specific humidity at saturation in [kg/kg] 
    606       !!              Based on accurate estimate of "e_sat" 
    607       !!              aka saturation water vapor (Goff, 1957) 
    608       !! 
    609       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    610       !!---------------------------------------------------------------------------------- 
    611       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak    ! air temperature                       [K] 
    612       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp    ! sea level atmospheric pressure       [Pa] 
    613       REAL(wp), DIMENSION(jpi,jpj)             ::   q_sat   ! Specific humidity at saturation   [kg/kg] 
    614       ! 
    615       INTEGER  ::   ji, jj         ! dummy loop indices 
    616       REAL(wp) ::   ze_sat, ztmp   ! local scalar 
    617       !!---------------------------------------------------------------------------------- 
    618       ! 
    619       DO jj = 1, jpj 
    620          DO ji = 1, jpi 
    621             ! 
    622             ztmp = rt0 / ptak(ji,jj) 
    623             ! 
    624             ! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957) 
    625             ze_sat = 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(ptak(ji,jj)/rt0)        & 
    626                &    + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak(ji,jj)/rt0 - 1.)) )  & 
    627                &    + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614  ) 
    628                ! 
    629             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] 
    630             ! 
    631          END DO 
    632       END DO 
    633       ! 
    634    END FUNCTION q_sat 
    635  
    636  
    637    FUNCTION gamma_moist( ptak, pqa ) 
    638       !!---------------------------------------------------------------------------------- 
    639       !!                           ***  FUNCTION gamma_moist  *** 
    640       !! 
    641       !! ** Purpose : Compute the moist adiabatic lapse-rate. 
    642       !!     => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate 
    643       !!     => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html 
    644       !! 
    645       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    646       !!---------------------------------------------------------------------------------- 
    647       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak          ! air temperature       [K] 
    648       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa           ! specific humidity [kg/kg] 
    649       REAL(wp), DIMENSION(jpi,jpj)             ::   gamma_moist   ! moist adiabatic lapse-rate 
    650       ! 
    651       INTEGER  ::   ji, jj         ! dummy loop indices 
    652       REAL(wp) :: zrv, ziRT        ! local scalar 
    653       !!---------------------------------------------------------------------------------- 
    654       ! 
    655       DO jj = 1, jpj 
    656          DO ji = 1, jpi 
    657             zrv = pqa(ji,jj) / (1. - pqa(ji,jj)) 
    658             ziRT = 1. / (R_dry*ptak(ji,jj))    ! 1/RT 
    659             gamma_moist(ji,jj) = grav * ( 1. + rLevap*zrv*ziRT ) / ( Cp_dry + rLevap*rLevap*zrv*reps0*ziRT/ptak(ji,jj) ) 
    660          END DO 
    661       END DO 
    662       ! 
    663    END FUNCTION gamma_moist 
    664  
    665  
    666    FUNCTION L_vap( psst ) 
    667       !!--------------------------------------------------------------------------------- 
    668       !!                           ***  FUNCTION L_vap  *** 
    669       !! 
    670       !! ** Purpose : Compute the latent heat of vaporization of water from temperature 
    671       !! 
    672       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    673       !!---------------------------------------------------------------------------------- 
    674       REAL(wp), DIMENSION(jpi,jpj)             ::   L_vap   ! latent heat of vaporization   [J/kg] 
    675       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   psst   ! water temperature                [K] 
    676       !!---------------------------------------------------------------------------------- 
    677       ! 
    678       L_vap = (  2.501 - 0.00237 * ( psst(:,:) - rt0)  ) * 1.e6 
    679       ! 
    680    END FUNCTION L_vap 
     825         CALL prt_ctl(tab2d_1=zqlw , clinfo1=' blk_oce_2: zqlw  : ') 
     826         CALL prt_ctl(tab2d_1=zqla , clinfo1=' blk_oce_2: zqla  : ', tab2d_2=qsr  , clinfo2=' qsr   : ') 
     827         CALL prt_ctl(tab2d_1=emp  , clinfo1=' blk_oce_2: emp   : ') 
     828      ENDIF 
     829      ! 
     830   END SUBROUTINE blk_oce_2 
     831 
    681832 
    682833#if defined key_si3 
     
    684835   !!   'key_si3'                                       SI3 sea-ice model 
    685836   !!---------------------------------------------------------------------- 
    686    !!   blk_ice_tau : provide the air-ice stress 
    687    !!   blk_ice_flx : provide the heat and mass fluxes at air-ice interface 
     837   !!   blk_ice_ : provide the air-ice stress 
     838   !!   blk_ice_ : provide the heat and mass fluxes at air-ice interface 
    688839   !!   blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux) 
    689840   !!   Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag 
    690    !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag  
     841   !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag 
    691842   !!---------------------------------------------------------------------- 
    692843 
    693    SUBROUTINE blk_ice_tau 
    694       !!--------------------------------------------------------------------- 
    695       !!                     ***  ROUTINE blk_ice_tau  *** 
     844   SUBROUTINE blk_ice_1( pwndi, pwndj, ptair, phumi, pslp , puice, pvice, ptsui,  &   ! inputs 
     845      &                  putaui, pvtaui, pseni, pevpi, pssqi, pcd_dui             )   ! optional outputs 
     846      !!--------------------------------------------------------------------- 
     847      !!                     ***  ROUTINE blk_ice_1  *** 
    696848      !! 
    697849      !! ** Purpose :   provide the surface boundary condition over sea-ice 
     
    701853      !!                NB: ice drag coefficient is assumed to be a constant 
    702854      !!--------------------------------------------------------------------- 
     855      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pslp    ! sea-level pressure [Pa] 
     856      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pwndi   ! atmospheric wind at T-point [m/s] 
     857      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pwndj   ! atmospheric wind at T-point [m/s] 
     858      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   ptair   ! atmospheric wind at T-point [m/s] 
     859      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   phumi   ! atmospheric wind at T-point [m/s] 
     860      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   puice   ! sea-ice velocity on I or C grid [m/s] 
     861      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pvice   ! " 
     862      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   ptsui   ! sea-ice surface temperature [K] 
     863      REAL(wp) , INTENT(  out), DIMENSION(:,:  ), OPTIONAL ::   putaui  ! if ln_blk 
     864      REAL(wp) , INTENT(  out), DIMENSION(:,:  ), OPTIONAL ::   pvtaui  ! if ln_blk 
     865      REAL(wp) , INTENT(  out), DIMENSION(:,:  ), OPTIONAL ::   pseni   ! if ln_abl 
     866      REAL(wp) , INTENT(  out), DIMENSION(:,:  ), OPTIONAL ::   pevpi   ! if ln_abl 
     867      REAL(wp) , INTENT(  out), DIMENSION(:,:  ), OPTIONAL ::   pssqi   ! if ln_abl 
     868      REAL(wp) , INTENT(  out), DIMENSION(:,:  ), OPTIONAL ::   pcd_dui ! if ln_abl 
     869      ! 
    703870      INTEGER  ::   ji, jj    ! dummy loop indices 
    704       REAL(wp) ::   zwndi_f , zwndj_f, zwnorm_f   ! relative wind module and components at F-point 
    705871      REAL(wp) ::   zwndi_t , zwndj_t             ! relative wind components at T-point 
    706       REAL(wp), DIMENSION(jpi,jpj) ::   zrhoa     ! transfer coefficient for momentum      (tau) 
    707       !!--------------------------------------------------------------------- 
    708       ! 
    709       ! set transfer coefficients to default sea-ice values 
    710       Cd_atm(:,:) = Cd_ice 
    711       Ch_atm(:,:) = Cd_ice 
    712       Ce_atm(:,:) = Cd_ice 
    713  
    714       wndm_ice(:,:) = 0._wp      !!gm brutal.... 
     872      REAL(wp) ::   zootm_su                      ! sea-ice surface mean temperature 
     873      REAL(wp) ::   zztmp1, zztmp2                ! temporary arrays 
     874      REAL(wp), DIMENSION(jpi,jpj) ::   zcd_dui   ! transfer coefficient for momentum      (tau) 
     875      !!--------------------------------------------------------------------- 
     876      ! 
    715877 
    716878      ! ------------------------------------------------------------ ! 
     
    720882      DO jj = 2, jpjm1 
    721883         DO ji = fs_2, fs_jpim1   ! vect. opt. 
    722             zwndi_t = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) )  ) 
    723             zwndj_t = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji  ,jj-1) + v_ice(ji,jj) )  ) 
     884            zwndi_t = (  pwndi(ji,jj) - rn_vfac * 0.5_wp * ( puice(ji-1,jj  ) + puice(ji,jj) )  ) 
     885            zwndj_t = (  pwndj(ji,jj) - rn_vfac * 0.5_wp * ( pvice(ji  ,jj-1) + pvice(ji,jj) )  ) 
    724886            wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 
    725887         END DO 
     
    729891      ! Make ice-atm. drag dependent on ice concentration 
    730892      IF    ( ln_Cd_L12 ) THEN   ! calculate new drag from Lupkes(2012) equations 
    731          CALL Cdn10_Lupkes2012( Cd_atm ) 
    732          Ch_atm(:,:) = Cd_atm(:,:)       ! momentum and heat transfer coef. are considered identical 
     893         CALL Cdn10_Lupkes2012( Cd_ice ) 
     894         Ch_ice(:,:) = Cd_ice(:,:)       ! momentum and heat transfer coef. are considered identical 
     895         Ce_ice(:,:) = Cd_ice(:,:) 
    733896      ELSEIF( ln_Cd_L15 ) THEN   ! calculate new drag from Lupkes(2015) equations 
    734          CALL Cdn10_Lupkes2015( Cd_atm, Ch_atm )  
    735       ENDIF 
    736  
    737 !!      CALL iom_put( "Cd_ice", Cd_atm)  ! output value of pure ice-atm. transfer coef. 
    738 !!      CALL iom_put( "Ch_ice", Ch_atm)  ! output value of pure ice-atm. transfer coef. 
     897         CALL Cdn10_Lupkes2015( ptsui, pslp, Cd_ice, Ch_ice ) 
     898         Ce_ice(:,:) = Ch_ice(:,:)       ! sensible and latent heat transfer coef. are considered identical 
     899      ENDIF 
     900 
     901      !! IF ( iom_use("Cd_ice") ) CALL iom_put("Cd_ice", Cd_ice)   ! output value of pure ice-atm. transfer coef. 
     902      !! IF ( iom_use("Ch_ice") ) CALL iom_put("Ch_ice", Ch_ice)   ! output value of pure ice-atm. transfer coef. 
    739903 
    740904      ! local scalars ( place there for vector optimisation purposes) 
    741       ! Computing density of air! Way denser that 1.2 over sea-ice !!! 
    742       zrhoa (:,:) =  rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1)) 
    743  
    744       !!gm brutal.... 
    745       utau_ice  (:,:) = 0._wp 
    746       vtau_ice  (:,:) = 0._wp 
    747       !!gm end 
    748  
    749       ! ------------------------------------------------------------ ! 
    750       !    Wind stress relative to the moving ice ( U10m - U_ice )   ! 
    751       ! ------------------------------------------------------------ ! 
    752       ! C-grid ice dynamics :   U & V-points (same as ocean) 
    753       DO jj = 2, jpjm1 
    754          DO ji = fs_2, fs_jpim1   ! vect. opt. 
    755             utau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )            & 
    756                &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) ) 
    757             vtau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )            & 
    758                &          * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) ) 
     905      !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) 
     906      zcd_dui(:,:) = wndm_ice(:,:) * Cd_ice(:,:) 
     907 
     908      IF( ln_blk ) THEN 
     909         ! ------------------------------------------------------------ ! 
     910         !    Wind stress relative to the moving ice ( U10m - U_ice )   ! 
     911         ! ------------------------------------------------------------ ! 
     912         ! C-grid ice dynamics :   U & V-points (same as ocean) 
     913         DO jj = 2, jpjm1 
     914            DO ji = fs_2, fs_jpim1   ! vect. opt. 
     915               putaui(ji,jj) = 0.5_wp * (  rhoa(ji+1,jj) * zcd_dui(ji+1,jj)             & 
     916                  &                      + rhoa(ji  ,jj) * zcd_dui(ji  ,jj)  )          & 
     917                  &         * ( 0.5_wp * ( pwndi(ji+1,jj) + pwndi(ji,jj) ) - rn_vfac * puice(ji,jj) ) 
     918               pvtaui(ji,jj) = 0.5_wp * (  rhoa(ji,jj+1) * zcd_dui(ji,jj+1)             & 
     919                  &                      + rhoa(ji,jj  ) * zcd_dui(ji,jj  )  )          & 
     920                  &         * ( 0.5_wp * ( pwndj(ji,jj+1) + pwndj(ji,jj) ) - rn_vfac * pvice(ji,jj) ) 
     921            END DO 
    759922         END DO 
    760       END DO 
    761       CALL lbc_lnk_multi( 'sbcblk', utau_ice, 'U', -1., vtau_ice, 'V', -1. ) 
    762       ! 
    763       ! 
    764       IF(ln_ctl) THEN 
    765          CALL prt_ctl(tab2d_1=utau_ice  , clinfo1=' blk_ice: utau_ice : ', tab2d_2=vtau_ice  , clinfo2=' vtau_ice : ') 
    766          CALL prt_ctl(tab2d_1=wndm_ice  , clinfo1=' blk_ice: wndm_ice : ') 
    767       ENDIF 
    768       ! 
    769    END SUBROUTINE blk_ice_tau 
    770  
    771  
    772    SUBROUTINE blk_ice_flx( ptsu, phs, phi, palb ) 
    773       !!--------------------------------------------------------------------- 
    774       !!                     ***  ROUTINE blk_ice_flx  *** 
     923         CALL lbc_lnk_multi( 'sbcblk', putaui, 'U', -1., pvtaui, 'V', -1. ) 
     924         ! 
     925         IF(ln_ctl)   CALL prt_ctl( tab2d_1=putaui  , clinfo1=' blk_ice: putaui : '   & 
     926            &                     , tab2d_2=pvtaui  , clinfo2='          pvtaui : ' ) 
     927      ELSE 
     928         zztmp1 = 11637800.0_wp 
     929         zztmp2 =    -5897.8_wp 
     930         DO jj = 1, jpj 
     931            DO ji = 1, jpi 
     932               pcd_dui(ji,jj) = zcd_dui (ji,jj) 
     933               pseni  (ji,jj) = wndm_ice(ji,jj) * Ch_ice(ji,jj) 
     934               pevpi  (ji,jj) = wndm_ice(ji,jj) * Ce_ice(ji,jj) 
     935               zootm_su       = zztmp2 / ptsui(ji,jj)   ! ptsui is in K (it can't be zero ??) 
     936               pssqi  (ji,jj) = zztmp1 * EXP( zootm_su ) / rhoa(ji,jj) 
     937            END DO 
     938         END DO 
     939      ENDIF 
     940      ! 
     941      IF(ln_ctl)  CALL prt_ctl(tab2d_1=wndm_ice  , clinfo1=' blk_ice: wndm_ice : ') 
     942      ! 
     943   END SUBROUTINE blk_ice_1 
     944 
     945 
     946   SUBROUTINE blk_ice_2( ptsu, phs, phi, palb, ptair, phumi, pslp, pqlw, pprec, psnow  ) 
     947      !!--------------------------------------------------------------------- 
     948      !!                     ***  ROUTINE blk_ice_2  *** 
    775949      !! 
    776950      !! ** Purpose :   provide the heat and mass fluxes at air-ice interface 
     
    782956      !! caution : the net upward water flux has with mm/day unit 
    783957      !!--------------------------------------------------------------------- 
    784       REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   ptsu   ! sea ice surface temperature 
     958      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   ptsu   ! sea ice surface temperature [K] 
    785959      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phs    ! snow thickness 
    786960      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phi    ! ice thickness 
    787961      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   palb   ! ice albedo (all skies) 
     962      REAL(wp), DIMENSION(:,:  ), INTENT(in)  ::   ptair 
     963      REAL(wp), DIMENSION(:,:  ), INTENT(in)  ::   phumi 
     964      REAL(wp), DIMENSION(:,:  ), INTENT(in)  ::   pslp 
     965      REAL(wp), DIMENSION(:,:  ), INTENT(in)  ::   pqlw 
     966      REAL(wp), DIMENSION(:,:  ), INTENT(in)  ::   pprec 
     967      REAL(wp), DIMENSION(:,:  ), INTENT(in)  ::   psnow 
    788968      !! 
    789969      INTEGER  ::   ji, jj, jl               ! dummy loop indices 
    790970      REAL(wp) ::   zst3                     ! local variable 
    791971      REAL(wp) ::   zcoef_dqlw, zcoef_dqla   !   -      - 
    792       REAL(wp) ::   zztmp, z1_rLsub           !   -      - 
     972      REAL(wp) ::   zztmp, zztmp2, z1_rLsub  !   -      - 
    793973      REAL(wp) ::   zfr1, zfr2               ! local variables 
    794974      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_st         ! inverse of surface temperature 
     
    798978      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_dqsb        ! sensible  heat sensitivity over ice 
    799979      REAL(wp), DIMENSION(jpi,jpj)     ::   zevap, zsnw   ! evaporation and snw distribution after wind blowing (SI3) 
    800       REAL(wp), DIMENSION(jpi,jpj)     ::   zrhoa 
    801       !!--------------------------------------------------------------------- 
    802       ! 
    803       zcoef_dqlw = 4.0 * 0.95 * Stef             ! local scalars 
    804       zcoef_dqla = -Ls * 11637800. * (-5897.8) 
    805       ! 
    806       zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 
     980      REAL(wp), DIMENSION(jpi,jpj)     ::   zqair         ! specific humidity of air at z=rn_zqt [kg/kg] !LB 
     981      !!--------------------------------------------------------------------- 
     982      ! 
     983      zcoef_dqlw = 4._wp * 0.95_wp * stefan             ! local scalars 
     984      zcoef_dqla = -rLsub * 11637800._wp * (-5897.8_wp) !LB: BAD! 
     985      ! 
     986      SELECT CASE( nhumi ) 
     987      CASE( np_humi_sph ) 
     988         zqair(:,:) =  phumi(:,:)      ! what we read in file is already a spec. humidity! 
     989      CASE( np_humi_dpt ) 
     990         zqair(:,:) = q_sat( phumi(:,:), pslp ) 
     991      CASE( np_humi_rlh ) 
     992         zqair(:,:) = q_air_rh( 0.01_wp*phumi(:,:), ptair(:,:), pslp(:,:) ) !LB: 0.01 => RH is % percent in file 
     993      END SELECT 
    807994      ! 
    808995      zztmp = 1. / ( 1. - albo ) 
    809       WHERE( ptsu(:,:,:) /= 0._wp )   ;   z1_st(:,:,:) = 1._wp / ptsu(:,:,:) 
    810       ELSEWHERE                       ;   z1_st(:,:,:) = 0._wp 
     996      WHERE( ptsu(:,:,:) /= 0._wp ) 
     997         z1_st(:,:,:) = 1._wp / ptsu(:,:,:) 
     998      ELSEWHERE 
     999         z1_st(:,:,:) = 0._wp 
    8111000      END WHERE 
    8121001      !                                     ! ========================== ! 
     
    8221011               qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj) 
    8231012               ! Long  Wave (lw) 
    824                z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) 
     1013               z_qlw(ji,jj,jl) = 0.95 * ( pqlw(ji,jj) - stefan * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) 
    8251014               ! lw sensitivity 
    8261015               z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3 
     
    8301019               ! ----------------------------! 
    8311020 
    832                ! ... turbulent heat fluxes with Ch_atm recalculated in blk_ice_tau 
     1021               ! ... turbulent heat fluxes with Ch_ice recalculated in blk_ice_1 
    8331022               ! Sensible Heat 
    834                z_qsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_atm(ji,jj) * wndm_ice(ji,jj) * (ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1)) 
     1023               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)) 
    8351024               ! Latent Heat 
    836                qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls  * Ch_atm(ji,jj) * wndm_ice(ji,jj) *  & 
    837                   &                ( 11637800. * EXP( -5897.8 * z1_st(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1) ) ) 
     1025               zztmp2 = EXP( -5897.8 * z1_st(ji,jj,jl) ) 
     1026               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa(ji,jj) * rLsub  * Ce_ice(ji,jj) * wndm_ice(ji,jj) *  & 
     1027                  &                ( 11637800. * zztmp2 / rhoa(ji,jj) - zqair(ji,jj) ) ) 
    8381028               ! Latent heat sensitivity for ice (Dqla/Dt) 
    8391029               IF( qla_ice(ji,jj,jl) > 0._wp ) THEN 
    840                   dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * Ch_atm(ji,jj) * wndm_ice(ji,jj) *  & 
    841                      &                 z1_st(ji,jj,jl)*z1_st(ji,jj,jl) * EXP(-5897.8 * z1_st(ji,jj,jl)) 
     1030                  dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * Ce_ice(ji,jj) * wndm_ice(ji,jj) *  & 
     1031                     &                 z1_st(ji,jj,jl) * z1_st(ji,jj,jl) * zztmp2 
    8421032               ELSE 
    8431033                  dqla_ice(ji,jj,jl) = 0._wp 
     
    8451035 
    8461036               ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice) 
    847                z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_atm(ji,jj) * wndm_ice(ji,jj) 
     1037               z_dqsb(ji,jj,jl) = rhoa(ji,jj) * rCp_air * Ch_ice(ji,jj) * wndm_ice(ji,jj) 
    8481038 
    8491039               ! ----------------------------! 
     
    8601050      END DO 
    8611051      ! 
    862       tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac * tmask(:,:,1)  ! total precipitation [kg/m2/s] 
    863       sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac * tmask(:,:,1)  ! solid precipitation [kg/m2/s] 
    864       CALL iom_put( 'snowpre', sprecip )                    ! Snow precipitation 
    865       CALL iom_put( 'precip' , tprecip )                    ! Total precipitation 
     1052      tprecip(:,:) = pprec(:,:) * rn_pfac * tmask(:,:,1)  ! total precipitation [kg/m2/s] 
     1053      sprecip(:,:) = psnow(:,:) * rn_pfac * tmask(:,:,1)  ! solid precipitation [kg/m2/s] 
     1054      CALL iom_put( 'snowpre', sprecip )                  ! Snow precipitation 
     1055      CALL iom_put( 'precip' , tprecip )                  ! Total precipitation 
    8661056 
    8671057      ! --- evaporation --- ! 
     
    8801070      ! --- heat flux associated with emp --- ! 
    8811071      qemp_oce(:,:) = - ( 1._wp - at_i_b(:,:) ) * zevap(:,:) * sst_m(:,:) * rcp                  & ! evap at sst 
    882          &          + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp  & ! liquid precip at Tair 
     1072         &          + ( tprecip(:,:) - sprecip(:,:) ) * ( ptair(:,:) - rt0 ) * rcp               & ! liquid precip at Tair 
    8831073         &          +   sprecip(:,:) * ( 1._wp - zsnw ) *                                        & ! solid precip at min(Tair,Tsnow) 
    884          &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus ) 
     1074         &              ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus ) 
    8851075      qemp_ice(:,:) =   sprecip(:,:) * zsnw *                                                    & ! solid precip (only) 
    886          &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus ) 
     1076         &              ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus ) 
    8871077 
    8881078      ! --- total solar and non solar fluxes --- ! 
     
    8921082 
    8931083      ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- ! 
    894       qprec_ice(:,:) = rhos * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus ) 
     1084      qprec_ice(:,:) = rhos * ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus ) 
    8951085 
    8961086      ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) --- 
    8971087      DO jl = 1, jpl 
    8981088         qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * rcpi * tmask(:,:,1) ) 
    899          !                         ! But we do not have Tice => consider it at 0degC => evap=0  
     1089         !                         ! But we do not have Tice => consider it at 0degC => evap=0 
    9001090      END DO 
    9011091 
     
    9041094      zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )            ! zfr2 such that zfr1 + zfr2 to equal 1 
    9051095      ! 
    906       WHERE    ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) <  0.1_wp )       ! linear decrease from hi=0 to 10cm   
     1096      WHERE    ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) <  0.1_wp )       ! linear decrease from hi=0 to 10cm 
    9071097         qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) ) 
    9081098      ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp )       ! constant (zfr1) when hi>10cm 
    9091099         qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * zfr1 
    9101100      ELSEWHERE                                                         ! zero when hs>0 
    911          qtr_ice_top(:,:,:) = 0._wp  
     1101         qtr_ice_top(:,:,:) = 0._wp 
    9121102      END WHERE 
    9131103      ! 
     
    9211111      ENDIF 
    9221112      ! 
    923    END SUBROUTINE blk_ice_flx 
    924     
     1113   END SUBROUTINE blk_ice_2 
     1114 
    9251115 
    9261116   SUBROUTINE blk_ice_qcn( ld_virtual_itd, ptsu, ptb, phs, phi ) 
     
    9311121      !!                to force sea ice / snow thermodynamics 
    9321122      !!                in the case conduction flux is emulated 
    933       !!                 
     1123      !! 
    9341124      !! ** Method  :   compute surface energy balance assuming neglecting heat storage 
    9351125      !!                following the 0-layer Semtner (1976) approach 
     
    9561146      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zgfac   ! enhanced conduction factor 
    9571147      !!--------------------------------------------------------------------- 
    958        
     1148 
    9591149      ! -------------------------------------! 
    9601150      !      I   Enhanced conduction factor  ! 
     
    9641154      ! 
    9651155      zgfac(:,:,:) = 1._wp 
    966        
     1156 
    9671157      IF( ld_virtual_itd ) THEN 
    9681158         ! 
     
    9701160         zfac2 = EXP(1._wp) * 0.5_wp * zepsilon 
    9711161         zfac3 = 2._wp / zepsilon 
    972          !    
    973          DO jl = 1, jpl                 
     1162         ! 
     1163         DO jl = 1, jpl 
    9741164            DO jj = 1 , jpj 
    9751165               DO ji = 1, jpi 
     
    9791169            END DO 
    9801170         END DO 
    981          !       
    982       ENDIF 
    983        
     1171         ! 
     1172      ENDIF 
     1173 
    9841174      ! -------------------------------------------------------------! 
    9851175      !      II   Surface temperature and conduction flux            ! 
     
    9911181         DO jj = 1 , jpj 
    9921182            DO ji = 1, jpi 
    993                !                     
     1183               ! 
    9941184               zkeff_h = zfac * zgfac(ji,jj,jl) / &                                    ! Effective conductivity of the snow-ice system divided by thickness 
    9951185                  &      ( rcnd_i * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) ) 
     
    10081198               qns_ice(ji,jj,jl) = qns_ice(ji,jj,jl) + dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) 
    10091199               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) )  & 
    1010                              &   * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) ) 
     1200                  &   * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) ) 
    10111201 
    10121202               ! --- Diagnose the heat loss due to changing non-solar flux (as in icethd_zdf_bl99) --- ! 
    1013                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)  
     1203               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) 
    10141204 
    10151205            END DO 
    10161206         END DO 
    10171207         ! 
    1018       END DO  
    1019       !       
     1208      END DO 
     1209      ! 
    10201210   END SUBROUTINE blk_ice_qcn 
    1021     
    1022  
    1023    SUBROUTINE Cdn10_Lupkes2012( Cd ) 
     1211 
     1212 
     1213   SUBROUTINE Cdn10_Lupkes2012( pcd ) 
    10241214      !!---------------------------------------------------------------------- 
    10251215      !!                      ***  ROUTINE  Cdn10_Lupkes2012  *** 
    10261216      !! 
    1027       !! ** Purpose :    Recompute the neutral air-ice drag referenced at 10m  
     1217      !! ** Purpose :    Recompute the neutral air-ice drag referenced at 10m 
    10281218      !!                 to make it dependent on edges at leads, melt ponds and flows. 
    10291219      !!                 After some approximations, this can be resumed to a dependency 
    10301220      !!                 on ice concentration. 
    1031       !!                 
     1221      !! 
    10321222      !! ** Method :     The parameterization is taken from Lupkes et al. (2012) eq.(50) 
    10331223      !!                 with the highest level of approximation: level4, eq.(59) 
     
    10411231      !! 
    10421232      !!                 This new drag has a parabolic shape (as a function of A) starting at 
    1043       !!                 Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5  
     1233      !!                 Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5 
    10441234      !!                 and going down to Cdi(say 1.4e-3) for A=1 
    10451235      !! 
     
    10511241      !! 
    10521242      !!---------------------------------------------------------------------- 
    1053       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Cd 
     1243      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pcd 
    10541244      REAL(wp), PARAMETER ::   zCe   = 2.23e-03_wp 
    10551245      REAL(wp), PARAMETER ::   znu   = 1._wp 
     
    10661256 
    10671257      ! ice-atm drag 
    1068       Cd(:,:) = Cd_ice +  &                                                         ! pure ice drag 
    1069          &      zCe    * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp)  ! change due to sea-ice morphology 
    1070        
     1258      pcd(:,:) = rCd_ice +  &                                                         ! pure ice drag 
     1259         &      zCe     * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp)  ! change due to sea-ice morphology 
     1260 
    10711261   END SUBROUTINE Cdn10_Lupkes2012 
    10721262 
    10731263 
    1074    SUBROUTINE Cdn10_Lupkes2015( Cd, Ch ) 
     1264   SUBROUTINE Cdn10_Lupkes2015( ptm_su, pslp, pcd, pch ) 
    10751265      !!---------------------------------------------------------------------- 
    10761266      !!                      ***  ROUTINE  Cdn10_Lupkes2015  *** 
    10771267      !! 
    10781268      !! ** pUrpose :    Alternative turbulent transfert coefficients formulation 
    1079       !!                 between sea-ice and atmosphere with distinct momentum  
    1080       !!                 and heat coefficients depending on sea-ice concentration  
     1269      !!                 between sea-ice and atmosphere with distinct momentum 
     1270      !!                 and heat coefficients depending on sea-ice concentration 
    10811271      !!                 and atmospheric stability (no meltponds effect for now). 
    1082       !!                 
     1272      !! 
    10831273      !! ** Method :     The parameterization is adapted from Lupkes et al. (2015) 
    10841274      !!                 and ECHAM6 atmospheric model. Compared to Lupkes2012 scheme, 
    10851275      !!                 it considers specific skin and form drags (Andreas et al. 2010) 
    1086       !!                 to compute neutral transfert coefficients for both heat and  
     1276      !!                 to compute neutral transfert coefficients for both heat and 
    10871277      !!                 momemtum fluxes. Atmospheric stability effect on transfert 
    10881278      !!                 coefficient is also taken into account following Louis (1979). 
     
    10931283      !!---------------------------------------------------------------------- 
    10941284      ! 
    1095       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Cd 
    1096       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Ch 
    1097       REAL(wp), DIMENSION(jpi,jpj)            ::   ztm_su, zst, zqo_sat, zqi_sat 
     1285      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   ptm_su ! sea-ice surface temperature [K] 
     1286      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pslp   ! sea-level pressure [Pa] 
     1287      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pcd    ! momentum transfert coefficient 
     1288      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pch    ! heat transfert coefficient 
     1289      REAL(wp), DIMENSION(jpi,jpj)            ::   zst, zqo_sat, zqi_sat 
    10981290      ! 
    10991291      ! ECHAM6 constants 
     
    11231315      !!---------------------------------------------------------------------- 
    11241316 
    1125       ! mean temperature 
    1126       WHERE( at_i_b(:,:) > 1.e-20 )   ;   ztm_su(:,:) = SUM( t_su(:,:,:) * a_i_b(:,:,:) , dim=3 ) / at_i_b(:,:) 
    1127       ELSEWHERE                       ;   ztm_su(:,:) = rt0 
    1128       ENDWHERE 
    1129        
    11301317      ! Momentum Neutral Transfert Coefficients (should be a constant) 
    11311318      zCdn_form_tmp = zce10 * ( LOG( 10._wp / z0_form_ice + 1._wp ) / LOG( rn_zu / z0_form_ice + 1._wp ) )**2   ! Eq. 40 
    11321319      zCdn_skin_ice = ( vkarmn                                      / LOG( rn_zu / z0_skin_ice + 1._wp ) )**2   ! Eq. 7 
    1133       zCdn_ice      = zCdn_skin_ice   ! Eq. 7 (cf Lupkes email for details) 
     1320      zCdn_ice      = zCdn_skin_ice   ! Eq. 7 
    11341321      !zCdn_ice     = 1.89e-3         ! old ECHAM5 value (cf Eq. 32) 
    11351322 
    11361323      ! Heat Neutral Transfert Coefficients 
    1137       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 (cf Lupkes email for details) 
    1138       
     1324      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 
     1325 
    11391326      ! Atmospheric and Surface Variables 
    11401327      zst(:,:)     = sst_m(:,:) + rt0                                        ! convert SST from Celcius to Kelvin 
    1141       zqo_sat(:,:) = 0.98_wp * q_sat( zst(:,:)   , sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ocean [kg/kg] 
    1142       zqi_sat(:,:) = 0.98_wp * q_sat( ztm_su(:,:), sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ice   [kg/kg] 
     1328      zqo_sat(:,:) = rdct_qsat_salt * q_sat( zst(:,:)   , pslp(:,:) )   ! saturation humidity over ocean [kg/kg] 
     1329      zqi_sat(:,:) =                  q_sat( ptm_su(:,:), pslp(:,:) )   ! saturation humidity over ice   [kg/kg] 
    11431330      ! 
    11441331      DO jj = 2, jpjm1           ! reduced loop is necessary for reproducibility 
     
    11461333            ! Virtual potential temperature [K] 
    11471334            zthetav_os = zst(ji,jj)    * ( 1._wp + rctv0 * zqo_sat(ji,jj) )   ! over ocean 
    1148             zthetav_is = ztm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) )   ! ocean ice 
     1335            zthetav_is = ptm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) )   ! ocean ice 
    11491336            zthetav_zu = t_zu (ji,jj)  * ( 1._wp + rctv0 * q_zu(ji,jj)    )   ! at zu 
    1150              
     1337 
    11511338            ! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead) 
    11521339            zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj)     )**2   ! over ocean 
    11531340            zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2   ! over ice 
    1154              
     1341 
    11551342            ! Momentum and Heat Neutral Transfert Coefficients 
    11561343            zCdn_form_ice = zCdn_form_tmp * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta  ! Eq. 40 
    1157             zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) )               ! Eq. 53  
    1158                         
    1159             ! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead) 
     1344            zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) )               ! Eq. 53 
     1345 
     1346            ! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead ?) 
    11601347            z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water 
    1161             z0i = z0_skin_ice                                             ! over ice (cf Lupkes email for details) 
     1348            z0i = z0_skin_ice                                             ! over ice 
    11621349            IF( zrib_o <= 0._wp ) THEN 
    11631350               zfmw = 1._wp - zam * zrib_o / ( 1._wp + 3._wp * zc2 * Cdn_oce(ji,jj) * SQRT( -zrib_o * ( rn_zu / z0w + 1._wp ) ) )  ! Eq. 10 
     
    11681355               zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) )   ! Eq. 28 
    11691356            ENDIF 
    1170              
     1357 
    11711358            IF( zrib_i <= 0._wp ) THEN 
    11721359               zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp)))   ! Eq.  9 
     
    11761363               zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) )   ! Eq. 27 
    11771364            ENDIF 
    1178              
     1365 
    11791366            ! Momentum Transfert Coefficients (Eq. 38) 
    1180             Cd(ji,jj) = zCdn_skin_ice *   zfmi +  & 
     1367            pcd(ji,jj) = zCdn_skin_ice *   zfmi +  & 
    11811368               &        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) ) 
    1182              
     1369 
    11831370            ! Heat Transfert Coefficients (Eq. 49) 
    1184             Ch(ji,jj) = zChn_skin_ice *   zfhi +  & 
     1371            pch(ji,jj) = zChn_skin_ice *   zfhi +  & 
    11851372               &        zChn_form_ice * ( zfhi * at_i_b(ji,jj) + zfhw * ( 1._wp - at_i_b(ji,jj) ) ) / MAX( 1.e-06, at_i_b(ji,jj) ) 
    11861373            ! 
    11871374         END DO 
    11881375      END DO 
    1189       CALL lbc_lnk_multi( 'sbcblk', Cd, 'T',  1., Ch, 'T', 1. ) 
     1376      CALL lbc_lnk_multi( 'sbcblk', pcd, 'T',  1., pch, 'T', 1. ) 
    11901377      ! 
    11911378   END SUBROUTINE Cdn10_Lupkes2015 
Note: See TracChangeset for help on using the changeset viewer.