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

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

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

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

Location:
NEMO/trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • NEMO/trunk

    • Property svn:externals

      •  

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

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

    r12276 r12377  
    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_1     ! routine called in icesbc 
     72   PUBLIC   blk_ice_2     ! 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 
    134 #  include "vectopt_loop_substitute.h90" 
     130#  include "do_loop_substitute.h90" 
    135131   !!---------------------------------------------------------------------- 
    136132   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    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   
    182       REWIND( numnam_ref )                !* Namelist namsbc_blk in reference namelist : bulk parameters 
     179      !                             !** read bulk namelist 
    183180      READ  ( numnam_ref, namsbc_blk, IOSTAT = ios, ERR = 901) 
    184181901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc_blk in reference namelist' ) 
    185182      ! 
    186       REWIND( numnam_cfg )                !* Namelist namsbc_blk in configuration namelist : bulk parameters 
    187183      READ  ( numnam_cfg, namsbc_blk, IOSTAT = ios, ERR = 902 ) 
    188184902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namsbc_blk in configuration namelist' ) 
    … …  
    192188      !                             !** initialization of the chosen bulk formulae (+ check) 
    193189      !                                   !* select the bulk chosen in the namelist and check the choice 
    194                                                                ioptio = 0 
    195       IF( ln_NCAR      ) THEN   ;   nblk =  np_NCAR        ;   ioptio = ioptio + 1   ;   ENDIF 
    196       IF( ln_COARE_3p0 ) THEN   ;   nblk =  np_COARE_3p0   ;   ioptio = ioptio + 1   ;   ENDIF 
    197       IF( ln_COARE_3p5 ) THEN   ;   nblk =  np_COARE_3p5   ;   ioptio = ioptio + 1   ;   ENDIF 
    198       IF( ln_ECMWF     ) THEN   ;   nblk =  np_ECMWF       ;   ioptio = ioptio + 1   ;   ENDIF 
    199       ! 
     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 
    200203      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' ) 
    201233      ! 
    202234      IF( ln_dm2dc ) THEN                 !* check: diurnal cycle on Qsr 
    203235         IF( sn_qsr%freqh /= 24. )   CALL ctl_stop( 'sbc_blk_init: ln_dm2dc=T only with daily short-wave input' ) 
    204          IF( sn_qsr%ln_tint ) THEN  
     236         IF( sn_qsr%ln_tint ) THEN 
    205237            CALL ctl_warn( 'sbc_blk_init: ln_dm2dc=T daily qsr time interpolation done by sbcdcy module',   & 
    206238               &           '              ==> We force time interpolation = .false. for qsr' ) 
    … …  
    210242      !                                   !* set the bulk structure 
    211243      !                                      !- store namelist information in an array 
     244      IF( ln_blk ) jpfld = 9 
     245      IF( ln_abl ) jpfld = 11 
     246      ALLOCATE( slf_i(jpfld) ) 
     247      ! 
    212248      slf_i(jp_wndi) = sn_wndi   ;   slf_i(jp_wndj) = sn_wndj 
    213249      slf_i(jp_qsr ) = sn_qsr    ;   slf_i(jp_qlw ) = sn_qlw 
    214250      slf_i(jp_tair) = sn_tair   ;   slf_i(jp_humi) = sn_humi 
    215251      slf_i(jp_prec) = sn_prec   ;   slf_i(jp_snow) = sn_snow 
    216       slf_i(jp_slp)  = sn_slp    ;   slf_i(jp_tdif) = sn_tdif 
    217       ! 
    218       lhftau = ln_taudif                     !- add an extra field if HF stress is used 
    219       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 
    220256      ! 
    221257      !                                      !- allocate the bulk structure 
    222       ALLOCATE( sf(jfld), STAT=ierror ) 
     258      ALLOCATE( sf(jpfld), STAT=ierror ) 
    223259      IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_init: unable to allocate sf structure' ) 
    224       DO ifpr= 1, jfld 
    225          ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) ) 
    226          IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) ) 
    227          IF( slf_i(ifpr)%freqh > 0. .AND. MOD( NINT(3600. * slf_i(ifpr)%freqh), nn_fsbc * NINT(rdt) ) /= 0 )   & 
    228             &  CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rdt*nn_fsbc is NOT a submultiple of atmospheric forcing frequency.',   & 
    229             &                 '               This is not ideal. You should consider changing either rdt or nn_fsbc value...' ) 
    230  
     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 
    231281      END DO 
    232282      !                                      !- fill the bulk structure with namelist informations 
    233283      CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_init', 'surface boundary condition -- bulk formulae', 'namsbc_blk' ) 
    234284      ! 
    235       IF ( ln_wave ) THEN 
    236       !Activated wave module but neither drag nor stokes drift activated 
    237          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 
    238288            CALL ctl_stop( 'STOP',  'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauwoc=F, ln_stcor=F' ) 
    239       !drag coefficient read from wave model definable only with mfs bulk formulae and core  
    240          ELSEIF (ln_cdgw .AND. .NOT. ln_NCAR )       THEN        
    241              CALL ctl_stop( 'drag coefficient read from wave model definable only with NCAR and CORE bulk formulae') 
    242          ELSEIF (ln_stcor .AND. .NOT. ln_sdw)                             THEN 
    243              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') 
    244294         ENDIF 
    245295      ELSE 
    246       IF ( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor )                &  
    247          &   CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ',    & 
    248          &                  'with drag coefficient (ln_cdgw =T) '  ,                        & 
    249          &                  'or Stokes Drift (ln_sdw=T) ' ,                                 & 
    250          &                  'or ocean stress modification due to waves (ln_tauwoc=T) ',      &   
    251          &                  'or Stokes-Coriolis term (ln_stcori=T)'  ) 
    252       ENDIF  
    253       ! 
    254       !            
     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      ! 
    255316      IF(lwp) THEN                     !** Control print 
    256317         ! 
    257          WRITE(numout,*)                  !* namelist  
     318         WRITE(numout,*)                  !* namelist 
    258319         WRITE(numout,*) '   Namelist namsbc_blk (other than data information):' 
    259320         WRITE(numout,*) '      "NCAR"      algorithm   (Large and Yeager 2008)     ln_NCAR      = ', ln_NCAR 
    260321         WRITE(numout,*) '      "COARE 3.0" algorithm   (Fairall et al. 2003)       ln_COARE_3p0 = ', ln_COARE_3p0 
    261          WRITE(numout,*) '      "COARE 3.5" algorithm   (Edson et al. 2013)         ln_COARE_3p5 = ', ln_COARE_3p0 
    262          WRITE(numout,*) '      "ECMWF"     algorithm   (IFS cycle 31)              ln_ECMWF     = ', ln_ECMWF 
    263          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 
    264324         WRITE(numout,*) '      Air temperature and humidity reference height (m)   rn_zqt       = ', rn_zqt 
    265325         WRITE(numout,*) '      Wind vector reference height (m)                    rn_zu        = ', rn_zu 
    … …  
    275335         CASE( np_NCAR      )   ;   WRITE(numout,*) '   ==>>>   "NCAR" algorithm        (Large and Yeager 2008)' 
    276336         CASE( np_COARE_3p0 )   ;   WRITE(numout,*) '   ==>>>   "COARE 3.0" algorithm   (Fairall et al. 2003)' 
    277          CASE( np_COARE_3p5 )   ;   WRITE(numout,*) '   ==>>>   "COARE 3.5" algorithm   (Edson et al. 2013)' 
    278          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)' 
    279339         END SELECT 
    280340         ! 
     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         ! 
    281352      ENDIF 
    282353      ! 
    … …  
    291362      !!              (momentum, heat, freshwater and runoff) 
    292363      !! 
    293       !! ** Method  : (1) READ each fluxes in NetCDF files: 
    294       !!      the 10m wind velocity (i-component) (m/s)    at T-point 
    295       !!      the 10m wind velocity (j-component) (m/s)    at T-point 
    296       !!      the 10m or 2m specific humidity     ( % ) 
    297       !!      the solar heat                      (W/m2) 
    298       !!      the Long wave                       (W/m2) 
    299       !!      the 10m or 2m air temperature       (Kelvin) 
    300       !!      the total precipitation (rain+snow) (Kg/m2/s) 
    301       !!      the snow (solid prcipitation)       (kg/m2/s) 
    302       !!      the tau diff associated to HF tau   (N/m2)   at T-point   (ln_taudif=T) 
    303       !!              (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 
    304376      !! 
    305377      !!      C A U T I O N : never mask the surface stress fields 
    … …  
    318390      !!---------------------------------------------------------------------- 
    319391      INTEGER, INTENT(in) ::   kt   ! ocean time step 
    320       !!--------------------------------------------------------------------- 
     392      !!---------------------------------------------------------------------- 
     393      REAL(wp), DIMENSION(jpi,jpj) ::   zssq, zcd_du, zsen, zevp 
     394      REAL(wp) :: ztmp 
     395      !!---------------------------------------------------------------------- 
    321396      ! 
    322397      CALL fld_read( kt, nn_fsbc, sf )             ! input fields provided at the current time-step 
    323       ! 
     398 
     399      ! Sanity/consistence test on humidity at first time step to detect potential screw-up: 
     400      IF( kt == nit000 ) THEN 
     401         IF(lwp) WRITE(numout,*) '' 
     402#if defined key_agrif 
     403         IF(lwp) 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               IF (lwp) WRITE(numout,'("   Mean humidity value found on proc #",i6.6," is: ",f10.5)') 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         IF(lwp) WRITE(numout,*) ' === Sanity/consistence test on air humidity sucessfuly passed! ===' 
     425#endif 
     426         IF(lwp) WRITE(numout,*) '' 
     427      END IF !IF( kt == nit000 ) 
    324428      !                                            ! compute the surface ocean fluxes using bulk formulea 
    325       IF( MOD( kt - 1, nn_fsbc ) == 0 )   CALL blk_oce( kt, sf, sst_m, ssu_m, ssv_m ) 
    326  
     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            &                tsk_m, 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), tsk_m,                     &   !   <<= in 
     439            &                zsen, zevp )                                            !   <=> in out 
     440      ENDIF 
     441      ! 
    327442#if defined key_cice 
    328443      IF( MOD( kt - 1, nn_fsbc ) == 0 )   THEN 
    329444         qlw_ice(:,:,1)   = sf(jp_qlw )%fnow(:,:,1) 
    330          IF( ln_dm2dc ) THEN ; qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) 
    331          ELSE                ; qsr_ice(:,:,1) =          sf(jp_qsr)%fnow(:,:,1)  
    332          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 
    333450         tatm_ice(:,:)    = sf(jp_tair)%fnow(:,:,1) 
    334          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 
    335461         tprecip(:,:)     = sf(jp_prec)%fnow(:,:,1) * rn_pfac 
    336462         sprecip(:,:)     = sf(jp_snow)%fnow(:,:,1) * rn_pfac 
    … …  
    343469 
    344470 
    345    SUBROUTINE blk_oce( kt, sf, pst, pu, pv ) 
    346       !!--------------------------------------------------------------------- 
    347       !!                     ***  ROUTINE blk_oce  *** 
    348       !! 
    349       !! ** Purpose :   provide the momentum, heat and freshwater fluxes at 
    350       !!      the ocean surface at each time step 
    351       !! 
    352       !! ** Method  :   bulk formulea for the ocean using atmospheric 
    353       !!      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      &                  ptsk, 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                      [Celsius] 
     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(:,:) ::   ptsk   ! skin temp. (or SST if CS & WL not used)  [Celsius] 
     502      REAL(wp), INTENT(  out), DIMENSION(:,:) ::   pssq   ! specific humidity at pst                 [kg/kg] 
     503      REAL(wp), INTENT(  out), DIMENSION(:,:) ::   pcd_du ! Cd x |dU| at T-points                    [m/s] 
     504      REAL(wp), INTENT(  out), DIMENSION(:,:) ::   psen   ! Ch x |dU| at T-points                    [m/s] 
     505      REAL(wp), INTENT(  out), DIMENSION(:,:) ::   pevp   ! Ce x |dU| at T-points                    [m/s] 
     506      ! 
     507      INTEGER  ::   ji, jj               ! dummy loop indices 
     508      REAL(wp) ::   zztmp                ! local variable 
     509      REAL(wp), DIMENSION(jpi,jpj) ::   zwnd_i, zwnd_j    ! wind speed components at T-point 
     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      !                           ! Temporary conversion from Celcius to Kelvin (and set minimum value far above 0 K) 
     522      ptsk(:,:) = pst(:,:) + rt0  ! by default: skin temperature = "bulk SST" (will remain this way if NCAR algorithm used!) 
     523 
     524      ! ----------------------------------------------------------------------------- ! 
     525      !      0   Wind components and module at T-point relative to the moving ocean   ! 
     526      ! ----------------------------------------------------------------------------- ! 
     527 
     528      ! ... components ( U10m - U_oce ) at T-point (unmasked) 
     529#if defined key_cyclone 
     530      zwnd_i(:,:) = 0._wp 
     531      zwnd_j(:,:) = 0._wp 
     532      CALL wnd_cyc( kt, zwnd_i, zwnd_j )    ! add analytical tropical cyclone (Vincent et al. JGR 2012) 
     533      DO_2D_00_00 
     534         pwndi(ji,jj) = pwndi(ji,jj) + zwnd_i(ji,jj) 
     535         pwndj(ji,jj) = pwndj(ji,jj) + zwnd_j(ji,jj) 
     536      END_2D 
     537#endif 
     538      DO_2D_00_00 
     539         zwnd_i(ji,jj) = (  pwndi(ji,jj) - rn_vfac * 0.5 * ( pu(ji-1,jj  ) + pu(ji,jj) )  ) 
     540         zwnd_j(ji,jj) = (  pwndj(ji,jj) - rn_vfac * 0.5 * ( pv(ji  ,jj-1) + pv(ji,jj) )  ) 
     541      END_2D 
     542      CALL lbc_lnk_multi( 'sbcblk', zwnd_i, 'T', -1., zwnd_j, 'T', -1. ) 
     543      ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked) 
     544      wndm(:,:) = SQRT(  zwnd_i(:,:) * zwnd_i(:,:)   & 
     545         &             + zwnd_j(:,:) * zwnd_j(:,:)  ) * tmask(:,:,1) 
     546 
     547      ! ----------------------------------------------------------------------------- ! 
     548      !      I   Solar FLUX                                                           ! 
     549      ! ----------------------------------------------------------------------------- ! 
     550 
     551      ! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle                    ! Short Wave 
     552      zztmp = 1. - albo 
     553      IF( ln_dm2dc ) THEN 
     554         qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1) 
     555      ELSE 
     556         qsr(:,:) = zztmp *          sf(jp_qsr)%fnow(:,:,1)   * tmask(:,:,1) 
     557      ENDIF 
     558 
     559 
     560      ! ----------------------------------------------------------------------------- ! 
     561      !     II   Turbulent FLUXES                                                     ! 
     562      ! ----------------------------------------------------------------------------- ! 
     563 
     564      ! specific humidity at SST 
     565      pssq(:,:) = rdct_qsat_salt * q_sat( ptsk(:,:), pslp(:,:) ) 
     566 
     567      IF( ln_skin_cs .OR. ln_skin_wl ) THEN 
     568         !! Backup "bulk SST" and associated spec. hum. 
     569         zztmp1(:,:) = ptsk(:,:) 
     570         zztmp2(:,:) = pssq(:,:) 
     571      ENDIF 
     572 
     573      ! specific humidity of air at "rn_zqt" m above the sea 
     574      SELECT CASE( nhumi ) 
     575      CASE( np_humi_sph ) 
     576         zqair(:,:) = phumi(:,:)      ! what we read in file is already a spec. humidity! 
     577      CASE( np_humi_dpt ) 
     578         !IF(lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of d_air and slp !' !LBrm 
     579         zqair(:,:) = q_sat( phumi(:,:), pslp(:,:) ) 
     580      CASE( np_humi_rlh ) 
     581         !IF(lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of RH, t_air and slp !' !LBrm 
     582         zqair(:,:) = q_air_rh( 0.01_wp*phumi(:,:), ptair(:,:), pslp(:,:) ) !LB: 0.01 => RH is % percent in file 
     583      END SELECT 
     584      ! 
     585      ! potential temperature of air at "rn_zqt" m above the sea 
     586      IF( ln_abl ) THEN 
     587         ztpot = ptair(:,:) 
     588      ELSE 
     589         ! Estimate of potential temperature at z=rn_zqt, based on adiabatic lapse-rate 
     590         !    (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2 
     591         !    (since reanalysis products provide T at z, not theta !) 
     592         !#LB: because AGRIF hates functions that return something else than a scalar, need to 
     593         !     use scalar version of gamma_moist() ... 
     594         DO_2D_11_11 
     595            ztpot(ji,jj) = ptair(ji,jj) + gamma_moist( ptair(ji,jj), zqair(ji,jj) ) * rn_zqt 
     596         END_2D 
     597      ENDIF 
     598 
     599 
     600 
     601      !! Time to call the user-selected bulk parameterization for 
     602      !!  ==  transfer coefficients  ==!   Cd, Ch, Ce at T-point, and more... 
     603      SELECT CASE( nblk ) 
     604 
     605      CASE( np_NCAR      ) 
     606         CALL turb_ncar    ( rn_zqt, rn_zu, ptsk, ztpot, pssq, zqair, wndm,                              & 
     607            &                zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
     608 
     609      CASE( np_COARE_3p0 ) 
     610         CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, ptsk, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 
     611            &                zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,   & 
     612            &                Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) ) 
     613 
     614      CASE( np_COARE_3p6 ) 
     615         CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, ptsk, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 
     616            &                zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,   & 
     617            &                Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) ) 
     618 
     619      CASE( np_ECMWF     ) 
     620         CALL turb_ecmwf   ( kt, rn_zqt, rn_zu, ptsk, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl,  & 
     621            &                zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,   & 
     622            &                Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) ) 
     623 
     624      CASE DEFAULT 
     625         CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' ) 
     626 
     627      END SELECT 
     628 
     629      IF( ln_skin_cs .OR. ln_skin_wl ) THEN 
     630         !! ptsk and pssq have been updated!!! 
     631         !! 
     632         !! In the presence of sea-ice we forget about the cool-skin/warm-layer update of ptsk and pssq: 
     633         WHERE ( fr_i(:,:) > 0.001_wp ) 
     634            ! sea-ice present, we forget about the update, using what we backed up before call to turb_*() 
     635            ptsk(:,:) = zztmp1(:,:) 
     636            pssq(:,:) = zztmp2(:,:) 
     637         END WHERE 
     638      END IF 
     639 
     640      !!      CALL iom_put( "Cd_oce", zcd_oce)  ! output value of pure ocean-atm. transfer coef. 
     641      !!      CALL iom_put( "Ch_oce", zch_oce)  ! output value of pure ocean-atm. transfer coef. 
     642 
     643      IF( ABS(rn_zu - rn_zqt) < 0.1_wp ) THEN 
     644         !! If zu == zt, then ensuring once for all that: 
     645         t_zu(:,:) = ztpot(:,:) 
     646         q_zu(:,:) = zqair(:,:) 
     647      ENDIF 
     648 
     649 
     650      !  Turbulent fluxes over ocean  => BULK_FORMULA @ sbcblk_phy.F90 
     651      ! ------------------------------------------------------------- 
     652 
     653      IF( ln_abl ) THEN         !==  ABL formulation  ==!   multiplication by rho_air and turbulent fluxes computation done in ablstp 
     654         !! FL do we need this multiplication by tmask ... ??? 
     655         DO_2D_11_11 
     656            zztmp = zU_zu(ji,jj) !* tmask(ji,jj,1) 
     657            wndm(ji,jj)   = zztmp                   ! Store zU_zu in wndm to compute ustar2 in ablmod 
     658            pcd_du(ji,jj) = zztmp * zcd_oce(ji,jj) 
     659            psen(ji,jj)   = zztmp * zch_oce(ji,jj) 
     660            pevp(ji,jj)   = zztmp * zce_oce(ji,jj) 
     661         END_2D 
     662      ELSE                      !==  BLK formulation  ==!   turbulent fluxes computation 
     663         CALL BULK_FORMULA( rn_zu, ptsk(:,:), pssq(:,:), t_zu(:,:), q_zu(:,:), & 
     664            &               zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:),         & 
     665            &               wndm(:,:), zU_zu(:,:), pslp(:,:),                 & 
     666            &               taum(:,:), psen(:,:), zqla(:,:),                  & 
     667            &               pEvap=pevp(:,:), prhoa=rhoa(:,:) ) 
     668 
     669         zqla(:,:) = zqla(:,:) * tmask(:,:,1) 
     670         psen(:,:) = psen(:,:) * tmask(:,:,1) 
     671         taum(:,:) = taum(:,:) * tmask(:,:,1) 
     672         pevp(:,:) = pevp(:,:) * tmask(:,:,1) 
     673 
     674         ! Tau i and j component on T-grid points, using array "zcd_oce" as a temporary array... 
     675         zcd_oce = 0._wp 
     676         WHERE ( wndm > 0._wp ) zcd_oce = taum / wndm 
     677         zwnd_i = zcd_oce * zwnd_i 
     678         zwnd_j = zcd_oce * zwnd_j 
     679 
     680         CALL iom_put( "taum_oce", taum )   ! output wind stress module 
     681 
     682         ! ... utau, vtau at U- and V_points, resp. 
     683         !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines 
     684         !     Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves 
     685         DO_2D_10_10 
     686            utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj  ) ) & 
     687               &          * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1)) 
     688            vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji  ,jj+1) ) & 
     689               &          * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1)) 
     690         END_2D 
     691         CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. ) 
     692 
     693         IF(sn_cfctl%l_prtctl) THEN 
     694            CALL prt_ctl( tab2d_1=wndm  , clinfo1=' blk_oce_1: wndm   : ') 
     695            CALL prt_ctl( tab2d_1=utau  , clinfo1=' blk_oce_1: utau   : ', mask1=umask,   & 
     696               &          tab2d_2=vtau  , clinfo2='            vtau   : ', mask2=vmask ) 
     697         ENDIF 
     698         ! 
     699      ENDIF !IF( ln_abl ) 
     700       
     701      ptsk(:,:) = ( ptsk(:,:) - rt0 ) * tmask(:,:,1)  ! Back to Celsius 
     702             
     703      IF( ln_skin_cs .OR. ln_skin_wl ) THEN 
     704         CALL iom_put( "t_skin" ,  ptsk        )  ! T_skin in Celsius 
     705         CALL iom_put( "dt_skin" , ptsk - pst  )  ! T_skin - SST temperature difference... 
     706      ENDIF 
     707 
     708      IF(sn_cfctl%l_prtctl) THEN 
     709         CALL prt_ctl( tab2d_1=pevp  , clinfo1=' blk_oce_1: pevp   : ' ) 
     710         CALL prt_ctl( tab2d_1=psen  , clinfo1=' blk_oce_1: psen   : ' ) 
     711         CALL prt_ctl( tab2d_1=pssq  , clinfo1=' blk_oce_1: pssq   : ' ) 
     712      ENDIF 
     713      ! 
     714   END SUBROUTINE blk_oce_1 
     715 
     716 
     717   SUBROUTINE blk_oce_2( ptair, pqsr, pqlw, pprec,   &   ! <<= in 
     718      &                  psnow, ptsk, psen, pevp     )   ! <<= in 
     719      !!--------------------------------------------------------------------- 
     720      !!                     ***  ROUTINE blk_oce_2  *** 
     721      !! 
     722      !! ** Purpose :   finalize the momentum, heat and freshwater fluxes computation 
     723      !!                at the ocean surface at each time step knowing Cd, Ch, Ce and 
     724      !!                atmospheric variables (from ABL or external data) 
    354725      !! 
    355726      !! ** Outputs : - utau    : i-component of the stress at U-point  (N/m2) 
    … …  
    360731      !!              - qns     : Non Solar heat flux over the ocean    (W/m2) 
    361732      !!              - emp     : evaporation minus precipitation       (kg/m2/s) 
    362       !! 
    363       !!  ** Nota  :   sf has to be a dummy argument for AGRIF on NEC 
    364       !!--------------------------------------------------------------------- 
    365       INTEGER  , INTENT(in   )                 ::   kt    ! time step index 
    366       TYPE(fld), INTENT(inout), DIMENSION(:)   ::   sf    ! input data 
    367       REAL(wp) , INTENT(in)   , DIMENSION(:,:) ::   pst   ! surface temperature                      [Celcius] 
    368       REAL(wp) , INTENT(in)   , DIMENSION(:,:) ::   pu    ! surface current at U-point (i-component) [m/s] 
    369       REAL(wp) , INTENT(in)   , DIMENSION(:,:) ::   pv    ! surface current at V-point (j-component) [m/s] 
     733      !!--------------------------------------------------------------------- 
     734      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptair 
     735      REAL(wp), INTENT(in), DIMENSION(:,:) ::   pqsr 
     736      REAL(wp), INTENT(in), DIMENSION(:,:) ::   pqlw 
     737      REAL(wp), INTENT(in), DIMENSION(:,:) ::   pprec 
     738      REAL(wp), INTENT(in), DIMENSION(:,:) ::   psnow 
     739      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptsk   ! SKIN surface temperature   [Celsius] 
     740      REAL(wp), INTENT(in), DIMENSION(:,:) ::   psen 
     741      REAL(wp), INTENT(in), DIMENSION(:,:) ::   pevp 
    370742      ! 
    371743      INTEGER  ::   ji, jj               ! dummy loop indices 
    372       REAL(wp) ::   zztmp                ! local variable 
    373       REAL(wp), DIMENSION(jpi,jpj) ::   zwnd_i, zwnd_j    ! wind speed components at T-point 
    374       REAL(wp), DIMENSION(jpi,jpj) ::   zsq               ! specific humidity at pst 
    375       REAL(wp), DIMENSION(jpi,jpj) ::   zqlw, zqsb        ! long wave and sensible heat fluxes 
    376       REAL(wp), DIMENSION(jpi,jpj) ::   zqla, zevap       ! latent heat fluxes and evaporation 
    377       REAL(wp), DIMENSION(jpi,jpj) ::   zst               ! surface temperature in Kelvin 
    378       REAL(wp), DIMENSION(jpi,jpj) ::   zU_zu             ! bulk wind speed at height zu  [m/s] 
    379       REAL(wp), DIMENSION(jpi,jpj) ::   ztpot             ! potential temperature of air at z=rn_zqt [K] 
    380       REAL(wp), DIMENSION(jpi,jpj) ::   zrhoa             ! density of air   [kg/m^3] 
     744      REAL(wp) ::   zztmp,zz1,zz2,zz3    ! local variable 
     745      REAL(wp), DIMENSION(jpi,jpj) ::   ztskk             ! skin temp. in Kelvin  
     746      REAL(wp), DIMENSION(jpi,jpj) ::   zqlw              ! long wave and sensible heat fluxes       
     747      REAL(wp), DIMENSION(jpi,jpj) ::   zqla              ! latent heat fluxes and evaporation 
    381748      !!--------------------------------------------------------------------- 
    382749      ! 
    383750      ! local scalars ( place there for vector optimisation purposes) 
    384       zst(:,:) = pst(:,:) + rt0      ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K) 
    385  
     751 
     752 
     753      ztskk(:,:) = ptsk(:,:) + rt0  ! => ptsk in Kelvin rather than Celsius 
     754       
    386755      ! ----------------------------------------------------------------------------- ! 
    387       !      0   Wind components and module at T-point relative to the moving ocean   ! 
     756      !     III    Net longwave radiative FLUX                                        ! 
    388757      ! ----------------------------------------------------------------------------- ! 
    389758 
    390       ! ... components ( U10m - U_oce ) at T-point (unmasked) 
    391 !!gm    move zwnd_i (_j) set to zero  inside the key_cyclone ??? 
    392       zwnd_i(:,:) = 0._wp 
    393       zwnd_j(:,:) = 0._wp 
    394 #if defined key_cyclone 
    395       CALL wnd_cyc( kt, zwnd_i, zwnd_j )    ! add analytical tropical cyclone (Vincent et al. JGR 2012) 
    396       DO jj = 2, jpjm1 
    397          DO ji = fs_2, fs_jpim1   ! vect. opt. 
    398             sf(jp_wndi)%fnow(ji,jj,1) = sf(jp_wndi)%fnow(ji,jj,1) + zwnd_i(ji,jj) 
    399             sf(jp_wndj)%fnow(ji,jj,1) = sf(jp_wndj)%fnow(ji,jj,1) + zwnd_j(ji,jj) 
    400          END DO 
    401       END DO 
    402 #endif 
    403       DO jj = 2, jpjm1 
    404          DO ji = fs_2, fs_jpim1   ! vect. opt. 
    405             zwnd_i(ji,jj) = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pu(ji-1,jj  ) + pu(ji,jj) )  ) 
    406             zwnd_j(ji,jj) = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pv(ji  ,jj-1) + pv(ji,jj) )  ) 
    407          END DO 
    408       END DO 
    409       CALL lbc_lnk_multi( 'sbcblk', zwnd_i, 'T', -1., zwnd_j, 'T', -1. ) 
    410       ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked) 
    411       wndm(:,:) = SQRT(  zwnd_i(:,:) * zwnd_i(:,:)   & 
    412          &             + zwnd_j(:,:) * zwnd_j(:,:)  ) * tmask(:,:,1) 
     759      !! LB: now moved after Turbulent fluxes because must use the skin temperature rather that the SST 
     760      !! (ztskk is skin temperature if ln_skin_cs==.TRUE. .OR. ln_skin_wl==.TRUE.) 
     761      zqlw(:,:) = emiss_w * ( pqlw(:,:) - stefan*ztskk(:,:)*ztskk(:,:)*ztskk(:,:)*ztskk(:,:) ) * tmask(:,:,1)   ! Net radiative longwave flux 
     762 
     763      !  Latent flux over ocean 
     764      ! ----------------------- 
     765 
     766      ! use scalar version of L_vap() for AGRIF compatibility 
     767      DO_2D_11_11 
     768         zqla(ji,jj) = - L_vap( ztskk(ji,jj) ) * pevp(ji,jj)    ! Latent Heat flux !!GS: possibility to add a global qla to avoid recomputation after abl update 
     769      END_2D 
     770 
     771      IF(sn_cfctl%l_prtctl) THEN 
     772         CALL prt_ctl( tab2d_1=zqla  , clinfo1=' blk_oce_2: zqla   : ' ) 
     773         CALL prt_ctl( tab2d_1=zqlw  , clinfo1=' blk_oce_2: zqlw   : ', tab2d_2=qsr, clinfo2=' qsr : ' ) 
     774 
     775      ENDIF 
    413776 
    414777      ! ----------------------------------------------------------------------------- ! 
    415       !      I   Radiative FLUXES                                                     ! 
     778      !     IV    Total FLUXES                                                       ! 
    416779      ! ----------------------------------------------------------------------------- ! 
    417  
    418       ! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle                    ! Short Wave 
    419       zztmp = 1. - albo 
    420       IF( ln_dm2dc ) THEN   ;   qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1) 
    421       ELSE                  ;   qsr(:,:) = zztmp *          sf(jp_qsr)%fnow(:,:,1)   * tmask(:,:,1) 
    422       ENDIF 
    423  
    424       zqlw(:,:) = (  sf(jp_qlw)%fnow(:,:,1) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:)  ) * tmask(:,:,1)   ! Long  Wave 
    425  
    426       ! ----------------------------------------------------------------------------- ! 
    427       !     II    Turbulent FLUXES                                                    ! 
    428       ! ----------------------------------------------------------------------------- ! 
    429  
    430       ! ... specific humidity at SST and IST tmask( 
    431       zsq(:,:) = 0.98 * q_sat( zst(:,:), sf(jp_slp)%fnow(:,:,1) ) 
    432       !! 
    433       !! Estimate of potential temperature at z=rn_zqt, based on adiabatic lapse-rate 
    434       !!    (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2 
    435       !!    (since reanalysis products provide T at z, not theta !) 
    436       ztpot = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1) ) * rn_zqt 
    437  
    438       SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point 
    439       ! 
    440       CASE( np_NCAR      )   ;   CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! NCAR-COREv2 
    441          &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    442       CASE( np_COARE_3p0 )   ;   CALL turb_coare   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! COARE v3.0 
    443          &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    444       CASE( np_COARE_3p5 )   ;   CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! COARE v3.5 
    445          &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    446       CASE( np_ECMWF     )   ;   CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! ECMWF 
    447          &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    448       CASE DEFAULT 
    449          CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' ) 
    450       END SELECT 
    451  
    452       !                          ! Compute true air density : 
    453       IF( ABS(rn_zu - rn_zqt) > 0.01 ) THEN     ! At zu: (probably useless to remove zrho*grav*rn_zu from SLP...) 
    454          zrhoa(:,:) = rho_air( t_zu(:,:)              , q_zu(:,:)              , sf(jp_slp)%fnow(:,:,1) ) 
    455       ELSE                                      ! At zt: 
    456          zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 
    457       END IF 
    458  
    459 !!      CALL iom_put( "Cd_oce", Cd_atm)  ! output value of pure ocean-atm. transfer coef. 
    460 !!      CALL iom_put( "Ch_oce", Ch_atm)  ! output value of pure ocean-atm. transfer coef. 
    461  
    462       DO jj = 1, jpj             ! tau module, i and j component 
    463          DO ji = 1, jpi 
    464             zztmp = zrhoa(ji,jj)  * zU_zu(ji,jj) * Cd_atm(ji,jj)   ! using bulk wind speed 
    465             taum  (ji,jj) = zztmp * wndm  (ji,jj) 
    466             zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) 
    467             zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj) 
    468          END DO 
    469       END DO 
    470  
    471       !                          ! add the HF tau contribution to the wind stress module 
    472       IF( lhftau )   taum(:,:) = taum(:,:) + sf(jp_tdif)%fnow(:,:,1) 
    473  
    474       CALL iom_put( "taum_oce", taum )   ! output wind stress module 
    475  
    476       ! ... utau, vtau at U- and V_points, resp. 
    477       !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines 
    478       !     Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves 
    479       DO jj = 1, jpjm1 
    480          DO ji = 1, fs_jpim1 
    481             utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj  ) ) & 
    482                &          * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1)) 
    483             vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji  ,jj+1) ) & 
    484                &          * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1)) 
    485          END DO 
    486       END DO 
    487       CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. ) 
    488  
    489       !  Turbulent fluxes over ocean 
    490       ! ----------------------------- 
    491  
    492       ! zqla used as temporary array, for rho*U (common term of bulk formulae): 
    493       zqla(:,:) = zrhoa(:,:) * zU_zu(:,:) * tmask(:,:,1) 
    494  
    495       IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN 
    496          !! q_air and t_air are given at 10m (wind reference height) 
    497          zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - sf(jp_humi)%fnow(:,:,1)) ) ! Evaporation, using bulk wind speed 
    498          zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - ztpot(:,:)             )   ! Sensible Heat, using bulk wind speed 
    499       ELSE 
    500          !! q_air and t_air are not given at 10m (wind reference height) 
    501          ! Values of temp. and hum. adjusted to height of wind during bulk algorithm iteration must be used!!! 
    502          zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - q_zu(:,:) ) ) ! Evaporation, using bulk wind speed 
    503          zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - t_zu(:,:) )   ! Sensible Heat, using bulk wind speed 
    504       ENDIF 
    505  
    506       zqla(:,:) = L_vap(zst(:,:)) * zevap(:,:)     ! Latent Heat flux 
    507  
    508  
    509       IF(ln_ctl) THEN 
    510          CALL prt_ctl( tab2d_1=zqla  , clinfo1=' blk_oce: zqla   : ', tab2d_2=Ce_atm , clinfo2=' Ce_oce  : ' ) 
    511          CALL prt_ctl( tab2d_1=zqsb  , clinfo1=' blk_oce: zqsb   : ', tab2d_2=Ch_atm , clinfo2=' Ch_oce  : ' ) 
    512          CALL prt_ctl( tab2d_1=zqlw  , clinfo1=' blk_oce: zqlw   : ', tab2d_2=qsr, clinfo2=' qsr : ' ) 
    513          CALL prt_ctl( tab2d_1=zsq   , clinfo1=' blk_oce: zsq    : ', tab2d_2=zst, clinfo2=' zst : ' ) 
    514          CALL prt_ctl( tab2d_1=utau  , clinfo1=' blk_oce: utau   : ', mask1=umask,   & 
    515             &          tab2d_2=vtau  , clinfo2=           ' vtau : ', mask2=vmask ) 
    516          CALL prt_ctl( tab2d_1=wndm  , clinfo1=' blk_oce: wndm   : ') 
    517          CALL prt_ctl( tab2d_1=zst   , clinfo1=' blk_oce: zst    : ') 
    518       ENDIF 
    519  
    520       ! ----------------------------------------------------------------------------- ! 
    521       !     III    Total FLUXES                                                       ! 
    522       ! ----------------------------------------------------------------------------- ! 
    523       ! 
    524       emp (:,:) = (  zevap(:,:)                                          &   ! mass flux (evap. - precip.) 
    525          &         - sf(jp_prec)%fnow(:,:,1) * rn_pfac  ) * tmask(:,:,1) 
    526       ! 
    527       qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar 
    528          &     - sf(jp_snow)%fnow(:,:,1) * rn_pfac * rLfus                        &   ! remove latent melting heat for solid precip 
    529          &     - zevap(:,:) * pst(:,:) * rcp                                      &   ! remove evap heat content at SST 
    530          &     + ( sf(jp_prec)%fnow(:,:,1) - sf(jp_snow)%fnow(:,:,1) ) * rn_pfac  &   ! add liquid precip heat content at Tair 
    531          &     * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp                          & 
    532          &     + sf(jp_snow)%fnow(:,:,1) * rn_pfac                                &   ! add solid  precip heat content at min(Tair,Tsnow) 
    533          &     * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi 
     780      ! 
     781      emp (:,:) = (  pevp(:,:)                                       &   ! mass flux (evap. - precip.) 
     782         &         - pprec(:,:) * rn_pfac  ) * tmask(:,:,1) 
     783      ! 
     784      qns(:,:) = zqlw(:,:) + psen(:,:) + zqla(:,:)                   &   ! Downward Non Solar 
     785         &     - psnow(:,:) * rn_pfac * rLfus                        &   ! remove latent melting heat for solid precip 
     786         &     - pevp(:,:) * ptsk(:,:) * rcp                         &   ! remove evap heat content at SST 
     787         &     + ( pprec(:,:) - psnow(:,:) ) * rn_pfac               &   ! add liquid precip heat content at Tair 
     788         &     * ( ptair(:,:) - rt0 ) * rcp                          & 
     789         &     + psnow(:,:) * rn_pfac                                &   ! add solid  precip heat content at min(Tair,Tsnow) 
     790         &     * ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi 
    534791      qns(:,:) = qns(:,:) * tmask(:,:,1) 
    535792      ! 
    536793#if defined key_si3 
    537       qns_oce(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                ! non solar without emp (only needed by SI3) 
     794      qns_oce(:,:) = zqlw(:,:) + psen(:,:) + zqla(:,:)                             ! non solar without emp (only needed by SI3) 
    538795      qsr_oce(:,:) = qsr(:,:) 
    539796#endif 
    540797      ! 
     798      CALL iom_put( "rho_air"  , rhoa*tmask(:,:,1) )       ! output air density [kg/m^3] 
     799      CALL iom_put( "evap_oce" , pevp )                    ! evaporation 
     800      CALL iom_put( "qlw_oce"  , zqlw )                    ! output downward longwave heat over the ocean 
     801      CALL iom_put( "qsb_oce"  , psen )                    ! output downward sensible heat over the ocean 
     802      CALL iom_put( "qla_oce"  , zqla )                    ! output downward latent   heat over the ocean 
     803      tprecip(:,:) = pprec(:,:) * rn_pfac * tmask(:,:,1)   ! output total precipitation [kg/m2/s] 
     804      sprecip(:,:) = psnow(:,:) * rn_pfac * tmask(:,:,1)   ! output solid precipitation [kg/m2/s] 
     805      CALL iom_put( 'snowpre', sprecip )                   ! Snow 
     806      CALL iom_put( 'precip' , tprecip )                   ! Total precipitation 
     807      ! 
    541808      IF ( nn_ice == 0 ) THEN 
    542          CALL iom_put( "qlw_oce" ,   zqlw )                 ! output downward longwave heat over the ocean 
    543          CALL iom_put( "qsb_oce" , - zqsb )                 ! output downward sensible heat over the ocean 
    544          CALL iom_put( "qla_oce" , - zqla )                 ! output downward latent   heat over the ocean 
    545          CALL iom_put( "qemp_oce",   qns-zqlw+zqsb+zqla )   ! output downward heat content of E-P over the ocean 
    546          CALL iom_put( "qns_oce" ,   qns  )                 ! output downward non solar heat over the ocean 
    547          CALL iom_put( "qsr_oce" ,   qsr  )                 ! output downward solar heat over the ocean 
    548          CALL iom_put( "qt_oce"  ,   qns+qsr )              ! output total downward heat over the ocean 
    549          tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s] 
    550          sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s] 
    551          CALL iom_put( 'snowpre', sprecip )                 ! Snow 
    552          CALL iom_put( 'precip' , tprecip )                 ! Total precipitation 
    553       ENDIF 
    554       ! 
    555       IF(ln_ctl) THEN 
    556          CALL prt_ctl(tab2d_1=zqsb , clinfo1=' blk_oce: zqsb   : ', tab2d_2=zqlw , clinfo2=' zqlw  : ') 
    557          CALL prt_ctl(tab2d_1=zqla , clinfo1=' blk_oce: zqla   : ', tab2d_2=qsr  , clinfo2=' qsr   : ') 
    558          CALL prt_ctl(tab2d_1=pst  , clinfo1=' blk_oce: pst    : ', tab2d_2=emp  , clinfo2=' emp   : ') 
    559          CALL prt_ctl(tab2d_1=utau , clinfo1=' blk_oce: utau   : ', mask1=umask,   & 
    560             &         tab2d_2=vtau , clinfo2=              ' vtau  : ' , mask2=vmask ) 
    561       ENDIF 
    562       ! 
    563    END SUBROUTINE blk_oce 
    564  
    565  
    566  
    567    FUNCTION rho_air( ptak, pqa, pslp ) 
    568       !!------------------------------------------------------------------------------- 
    569       !!                           ***  FUNCTION rho_air  *** 
    570       !! 
    571       !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere 
    572       !! 
    573       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)  
    574       !!------------------------------------------------------------------------------- 
    575       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak      ! air temperature             [K] 
    576       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa       ! air specific humidity   [kg/kg] 
    577       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp      ! pressure in                [Pa] 
    578       REAL(wp), DIMENSION(jpi,jpj)             ::   rho_air   ! density of moist air   [kg/m^3] 
    579       !!------------------------------------------------------------------------------- 
    580       ! 
    581       rho_air = pslp / (  R_dry*ptak * ( 1._wp + rctv0*pqa )  ) 
    582       ! 
    583    END FUNCTION rho_air 
    584  
    585  
    586    FUNCTION cp_air( pqa ) 
    587       !!------------------------------------------------------------------------------- 
    588       !!                           ***  FUNCTION cp_air  *** 
    589       !! 
    590       !! ** Purpose : Compute specific heat (Cp) of moist air 
    591       !! 
    592       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    593       !!------------------------------------------------------------------------------- 
    594       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa      ! air specific humidity         [kg/kg] 
    595       REAL(wp), DIMENSION(jpi,jpj)             ::   cp_air   ! specific heat of moist air   [J/K/kg] 
    596       !!------------------------------------------------------------------------------- 
    597       ! 
    598       Cp_air = Cp_dry + Cp_vap * pqa 
    599       ! 
    600    END FUNCTION cp_air 
    601  
    602  
    603    FUNCTION q_sat( ptak, pslp ) 
    604       !!---------------------------------------------------------------------------------- 
    605       !!                           ***  FUNCTION q_sat  *** 
    606       !! 
    607       !! ** Purpose : Specific humidity at saturation in [kg/kg] 
    608       !!              Based on accurate estimate of "e_sat" 
    609       !!              aka saturation water vapor (Goff, 1957) 
    610       !! 
    611       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    612       !!---------------------------------------------------------------------------------- 
    613       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak    ! air temperature                       [K] 
    614       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp    ! sea level atmospheric pressure       [Pa] 
    615       REAL(wp), DIMENSION(jpi,jpj)             ::   q_sat   ! Specific humidity at saturation   [kg/kg] 
    616       ! 
    617       INTEGER  ::   ji, jj         ! dummy loop indices 
    618       REAL(wp) ::   ze_sat, ztmp   ! local scalar 
    619       !!---------------------------------------------------------------------------------- 
    620       ! 
    621       DO jj = 1, jpj 
    622          DO ji = 1, jpi 
    623             ! 
    624             ztmp = rt0 / ptak(ji,jj) 
    625             ! 
    626             ! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957) 
    627             ze_sat = 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(ptak(ji,jj)/rt0)        & 
    628                &    + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak(ji,jj)/rt0 - 1.)) )  & 
    629                &    + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614  ) 
    630                ! 
    631             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] 
    632             ! 
    633          END DO 
    634       END DO 
    635       ! 
    636    END FUNCTION q_sat 
    637  
    638  
    639    FUNCTION gamma_moist( ptak, pqa ) 
    640       !!---------------------------------------------------------------------------------- 
    641       !!                           ***  FUNCTION gamma_moist  *** 
    642       !! 
    643       !! ** Purpose : Compute the moist adiabatic lapse-rate. 
    644       !!     => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate 
    645       !!     => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html 
    646       !! 
    647       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    648       !!---------------------------------------------------------------------------------- 
    649       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak          ! air temperature       [K] 
    650       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa           ! specific humidity [kg/kg] 
    651       REAL(wp), DIMENSION(jpi,jpj)             ::   gamma_moist   ! moist adiabatic lapse-rate 
    652       ! 
    653       INTEGER  ::   ji, jj         ! dummy loop indices 
    654       REAL(wp) :: zrv, ziRT        ! local scalar 
    655       !!---------------------------------------------------------------------------------- 
    656       ! 
    657       DO jj = 1, jpj 
    658          DO ji = 1, jpi 
    659             zrv = pqa(ji,jj) / (1. - pqa(ji,jj)) 
    660             ziRT = 1. / (R_dry*ptak(ji,jj))    ! 1/RT 
    661             gamma_moist(ji,jj) = grav * ( 1. + rLevap*zrv*ziRT ) / ( Cp_dry + rLevap*rLevap*zrv*reps0*ziRT/ptak(ji,jj) ) 
    662          END DO 
    663       END DO 
    664       ! 
    665    END FUNCTION gamma_moist 
    666  
    667  
    668    FUNCTION L_vap( psst ) 
    669       !!--------------------------------------------------------------------------------- 
    670       !!                           ***  FUNCTION L_vap  *** 
    671       !! 
    672       !! ** Purpose : Compute the latent heat of vaporization of water from temperature 
    673       !! 
    674       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    675       !!---------------------------------------------------------------------------------- 
    676       REAL(wp), DIMENSION(jpi,jpj)             ::   L_vap   ! latent heat of vaporization   [J/kg] 
    677       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   psst   ! water temperature                [K] 
    678       !!---------------------------------------------------------------------------------- 
    679       ! 
    680       L_vap = (  2.501 - 0.00237 * ( psst(:,:) - rt0)  ) * 1.e6 
    681       ! 
    682    END FUNCTION L_vap 
     809         CALL iom_put( "qemp_oce" , qns-zqlw-psen-zqla )   ! output downward heat content of E-P over the ocean 
     810         CALL iom_put( "qns_oce"  ,   qns  )               ! output downward non solar heat over the ocean 
     811         CALL iom_put( "qsr_oce"  ,   qsr  )               ! output downward solar heat over the ocean 
     812         CALL iom_put( "qt_oce"   ,   qns+qsr )            ! output total downward heat over the ocean 
     813      ENDIF 
     814      ! 
     815      IF(sn_cfctl%l_prtctl) THEN 
     816         CALL prt_ctl(tab2d_1=zqlw , clinfo1=' blk_oce_2: zqlw  : ') 
     817         CALL prt_ctl(tab2d_1=zqla , clinfo1=' blk_oce_2: zqla  : ', tab2d_2=qsr  , clinfo2=' qsr   : ') 
     818         CALL prt_ctl(tab2d_1=emp  , clinfo1=' blk_oce_2: emp   : ') 
     819      ENDIF 
     820      ! 
     821   END SUBROUTINE blk_oce_2 
     822 
    683823 
    684824#if defined key_si3 
    … …  
    686826   !!   'key_si3'                                       SI3 sea-ice model 
    687827   !!---------------------------------------------------------------------- 
    688    !!   blk_ice_tau : provide the air-ice stress 
    689    !!   blk_ice_flx : provide the heat and mass fluxes at air-ice interface 
     828   !!   blk_ice_1  : provide the air-ice stress 
     829   !!   blk_ice_2  : provide the heat and mass fluxes at air-ice interface 
    690830   !!   blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux) 
    691831   !!   Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag 
    692    !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag  
     832   !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag 
    693833   !!---------------------------------------------------------------------- 
    694834 
    695    SUBROUTINE blk_ice_tau 
    696       !!--------------------------------------------------------------------- 
    697       !!                     ***  ROUTINE blk_ice_tau  *** 
     835   SUBROUTINE blk_ice_1( pwndi, pwndj, ptair, phumi, pslp , puice, pvice, ptsui,  &   ! inputs 
     836      &                  putaui, pvtaui, pseni, pevpi, pssqi, pcd_dui             )   ! optional outputs 
     837      !!--------------------------------------------------------------------- 
     838      !!                     ***  ROUTINE blk_ice_1  *** 
    698839      !! 
    699840      !! ** Purpose :   provide the surface boundary condition over sea-ice 
    … …  
    703844      !!                NB: ice drag coefficient is assumed to be a constant 
    704845      !!--------------------------------------------------------------------- 
     846      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pslp    ! sea-level pressure [Pa] 
     847      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pwndi   ! atmospheric wind at T-point [m/s] 
     848      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pwndj   ! atmospheric wind at T-point [m/s] 
     849      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   ptair   ! atmospheric wind at T-point [m/s] 
     850      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   phumi   ! atmospheric wind at T-point [m/s] 
     851      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   puice   ! sea-ice velocity on I or C grid [m/s] 
     852      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pvice   ! " 
     853      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   ptsui   ! sea-ice surface temperature [K] 
     854      REAL(wp) , INTENT(  out), DIMENSION(:,:  ), OPTIONAL ::   putaui  ! if ln_blk 
     855      REAL(wp) , INTENT(  out), DIMENSION(:,:  ), OPTIONAL ::   pvtaui  ! if ln_blk 
     856      REAL(wp) , INTENT(  out), DIMENSION(:,:  ), OPTIONAL ::   pseni   ! if ln_abl 
     857      REAL(wp) , INTENT(  out), DIMENSION(:,:  ), OPTIONAL ::   pevpi   ! if ln_abl 
     858      REAL(wp) , INTENT(  out), DIMENSION(:,:  ), OPTIONAL ::   pssqi   ! if ln_abl 
     859      REAL(wp) , INTENT(  out), DIMENSION(:,:  ), OPTIONAL ::   pcd_dui ! if ln_abl 
     860      ! 
    705861      INTEGER  ::   ji, jj    ! dummy loop indices 
    706       REAL(wp) ::   zwndi_f , zwndj_f, zwnorm_f   ! relative wind module and components at F-point 
    707862      REAL(wp) ::   zwndi_t , zwndj_t             ! relative wind components at T-point 
    708       REAL(wp), DIMENSION(jpi,jpj) ::   zrhoa     ! transfer coefficient for momentum      (tau) 
    709       !!--------------------------------------------------------------------- 
    710       ! 
    711       ! set transfer coefficients to default sea-ice values 
    712       Cd_atm(:,:) = Cd_ice 
    713       Ch_atm(:,:) = Cd_ice 
    714       Ce_atm(:,:) = Cd_ice 
    715  
    716       wndm_ice(:,:) = 0._wp      !!gm brutal.... 
     863      REAL(wp) ::   zootm_su                      ! sea-ice surface mean temperature 
     864      REAL(wp) ::   zztmp1, zztmp2                ! temporary arrays 
     865      REAL(wp), DIMENSION(jpi,jpj) ::   zcd_dui   ! transfer coefficient for momentum      (tau) 
     866      !!--------------------------------------------------------------------- 
     867      ! 
    717868 
    718869      ! ------------------------------------------------------------ ! 
    … …  
    720871      ! ------------------------------------------------------------ ! 
    721872      ! C-grid ice dynamics :   U & V-points (same as ocean) 
    722       DO jj = 2, jpjm1 
    723          DO ji = fs_2, fs_jpim1   ! vect. opt. 
    724             zwndi_t = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) )  ) 
    725             zwndj_t = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji  ,jj-1) + v_ice(ji,jj) )  ) 
    726             wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 
    727          END DO 
    728       END DO 
     873      DO_2D_00_00 
     874         zwndi_t = (  pwndi(ji,jj) - rn_vfac * 0.5_wp * ( puice(ji-1,jj  ) + puice(ji,jj) )  ) 
     875         zwndj_t = (  pwndj(ji,jj) - rn_vfac * 0.5_wp * ( pvice(ji  ,jj-1) + pvice(ji,jj) )  ) 
     876         wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 
     877      END_2D 
    729878      CALL lbc_lnk( 'sbcblk', wndm_ice, 'T',  1. ) 
    730879      ! 
    731880      ! Make ice-atm. drag dependent on ice concentration 
    732881      IF    ( ln_Cd_L12 ) THEN   ! calculate new drag from Lupkes(2012) equations 
    733          CALL Cdn10_Lupkes2012( Cd_atm ) 
    734          Ch_atm(:,:) = Cd_atm(:,:)       ! momentum and heat transfer coef. are considered identical 
     882         CALL Cdn10_Lupkes2012( Cd_ice ) 
     883         Ch_ice(:,:) = Cd_ice(:,:)       ! momentum and heat transfer coef. are considered identical 
     884         Ce_ice(:,:) = Cd_ice(:,:) 
    735885      ELSEIF( ln_Cd_L15 ) THEN   ! calculate new drag from Lupkes(2015) equations 
    736          CALL Cdn10_Lupkes2015( Cd_atm, Ch_atm )  
    737       ENDIF 
    738  
    739 !!      CALL iom_put( "Cd_ice", Cd_atm)  ! output value of pure ice-atm. transfer coef. 
    740 !!      CALL iom_put( "Ch_ice", Ch_atm)  ! output value of pure ice-atm. transfer coef. 
     886         CALL Cdn10_Lupkes2015( ptsui, pslp, Cd_ice, Ch_ice ) 
     887         Ce_ice(:,:) = Ch_ice(:,:)       ! sensible and latent heat transfer coef. are considered identical 
     888      ENDIF 
     889 
     890      !! IF ( iom_use("Cd_ice") ) CALL iom_put("Cd_ice", Cd_ice)   ! output value of pure ice-atm. transfer coef. 
     891      !! IF ( iom_use("Ch_ice") ) CALL iom_put("Ch_ice", Ch_ice)   ! output value of pure ice-atm. transfer coef. 
    741892 
    742893      ! local scalars ( place there for vector optimisation purposes) 
    743       ! Computing density of air! Way denser that 1.2 over sea-ice !!! 
    744       zrhoa (:,:) =  rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1)) 
    745  
    746       !!gm brutal.... 
    747       utau_ice  (:,:) = 0._wp 
    748       vtau_ice  (:,:) = 0._wp 
    749       !!gm end 
    750  
    751       ! ------------------------------------------------------------ ! 
    752       !    Wind stress relative to the moving ice ( U10m - U_ice )   ! 
    753       ! ------------------------------------------------------------ ! 
    754       ! C-grid ice dynamics :   U & V-points (same as ocean) 
    755       DO jj = 2, jpjm1 
    756          DO ji = fs_2, fs_jpim1   ! vect. opt. 
    757             utau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )            & 
    758                &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) ) 
    759             vtau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )            & 
    760                &          * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) ) 
    761          END DO 
    762       END DO 
    763       CALL lbc_lnk_multi( 'sbcblk', utau_ice, 'U', -1., vtau_ice, 'V', -1. ) 
    764       ! 
    765       ! 
    766       IF(ln_ctl) THEN 
    767          CALL prt_ctl(tab2d_1=utau_ice  , clinfo1=' blk_ice: utau_ice : ', tab2d_2=vtau_ice  , clinfo2=' vtau_ice : ') 
    768          CALL prt_ctl(tab2d_1=wndm_ice  , clinfo1=' blk_ice: wndm_ice : ') 
    769       ENDIF 
    770       ! 
    771    END SUBROUTINE blk_ice_tau 
    772  
    773  
    774    SUBROUTINE blk_ice_flx( ptsu, phs, phi, palb ) 
    775       !!--------------------------------------------------------------------- 
    776       !!                     ***  ROUTINE blk_ice_flx  *** 
     894      !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) 
     895      zcd_dui(:,:) = wndm_ice(:,:) * Cd_ice(:,:) 
     896 
     897      IF( ln_blk ) THEN 
     898         ! ------------------------------------------------------------ ! 
     899         !    Wind stress relative to the moving ice ( U10m - U_ice )   ! 
     900         ! ------------------------------------------------------------ ! 
     901         ! C-grid ice dynamics :   U & V-points (same as ocean) 
     902         DO_2D_00_00 
     903            putaui(ji,jj) = 0.5_wp * (  rhoa(ji+1,jj) * zcd_dui(ji+1,jj)             & 
     904               &                      + rhoa(ji  ,jj) * zcd_dui(ji  ,jj)  )          & 
     905               &         * ( 0.5_wp * ( pwndi(ji+1,jj) + pwndi(ji,jj) ) - rn_vfac * puice(ji,jj) ) 
     906            pvtaui(ji,jj) = 0.5_wp * (  rhoa(ji,jj+1) * zcd_dui(ji,jj+1)             & 
     907               &                      + rhoa(ji,jj  ) * zcd_dui(ji,jj  )  )          & 
     908               &         * ( 0.5_wp * ( pwndj(ji,jj+1) + pwndj(ji,jj) ) - rn_vfac * pvice(ji,jj) ) 
     909         END_2D 
     910         CALL lbc_lnk_multi( 'sbcblk', putaui, 'U', -1., pvtaui, 'V', -1. ) 
     911         ! 
     912         IF(sn_cfctl%l_prtctl)  CALL prt_ctl( tab2d_1=putaui  , clinfo1=' blk_ice: putaui : '   & 
     913            &                               , tab2d_2=pvtaui  , clinfo2='          pvtaui : ' ) 
     914      ELSE 
     915         zztmp1 = 11637800.0_wp 
     916         zztmp2 =    -5897.8_wp 
     917         DO_2D_11_11 
     918            pcd_dui(ji,jj) = zcd_dui (ji,jj) 
     919            pseni  (ji,jj) = wndm_ice(ji,jj) * Ch_ice(ji,jj) 
     920            pevpi  (ji,jj) = wndm_ice(ji,jj) * Ce_ice(ji,jj) 
     921            zootm_su       = zztmp2 / ptsui(ji,jj)   ! ptsui is in K (it can't be zero ??) 
     922            pssqi  (ji,jj) = zztmp1 * EXP( zootm_su ) / rhoa(ji,jj) 
     923         END_2D 
     924      ENDIF 
     925      ! 
     926      IF(sn_cfctl%l_prtctl)  CALL prt_ctl(tab2d_1=wndm_ice  , clinfo1=' blk_ice: wndm_ice : ') 
     927      ! 
     928   END SUBROUTINE blk_ice_1 
     929 
     930 
     931   SUBROUTINE blk_ice_2( ptsu, phs, phi, palb, ptair, phumi, pslp, pqlw, pprec, psnow  ) 
     932      !!--------------------------------------------------------------------- 
     933      !!                     ***  ROUTINE blk_ice_2  *** 
    777934      !! 
    778935      !! ** Purpose :   provide the heat and mass fluxes at air-ice interface 
    … …  
    784941      !! caution : the net upward water flux has with mm/day unit 
    785942      !!--------------------------------------------------------------------- 
    786       REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   ptsu   ! sea ice surface temperature 
     943      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   ptsu   ! sea ice surface temperature [K] 
    787944      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phs    ! snow thickness 
    788945      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phi    ! ice thickness 
    789946      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   palb   ! ice albedo (all skies) 
     947      REAL(wp), DIMENSION(:,:  ), INTENT(in)  ::   ptair 
     948      REAL(wp), DIMENSION(:,:  ), INTENT(in)  ::   phumi 
     949      REAL(wp), DIMENSION(:,:  ), INTENT(in)  ::   pslp 
     950      REAL(wp), DIMENSION(:,:  ), INTENT(in)  ::   pqlw 
     951      REAL(wp), DIMENSION(:,:  ), INTENT(in)  ::   pprec 
     952      REAL(wp), DIMENSION(:,:  ), INTENT(in)  ::   psnow 
    790953      !! 
    791954      INTEGER  ::   ji, jj, jl               ! dummy loop indices 
    792955      REAL(wp) ::   zst3                     ! local variable 
    793956      REAL(wp) ::   zcoef_dqlw, zcoef_dqla   !   -      - 
    794       REAL(wp) ::   zztmp, z1_rLsub           !   -      - 
     957      REAL(wp) ::   zztmp, zztmp2, z1_rLsub  !   -      - 
    795958      REAL(wp) ::   zfr1, zfr2               ! local variables 
    796959      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_st         ! inverse of surface temperature 
    … …  
    800963      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_dqsb        ! sensible  heat sensitivity over ice 
    801964      REAL(wp), DIMENSION(jpi,jpj)     ::   zevap, zsnw   ! evaporation and snw distribution after wind blowing (SI3) 
    802       REAL(wp), DIMENSION(jpi,jpj)     ::   zrhoa 
     965      REAL(wp), DIMENSION(jpi,jpj)     ::   zqair         ! specific humidity of air at z=rn_zqt [kg/kg] !LB 
    803966      REAL(wp), DIMENSION(jpi,jpj)     ::   ztmp, ztmp2 
    804967      !!--------------------------------------------------------------------- 
    805968      ! 
    806       zcoef_dqlw = 4.0 * 0.95 * Stef             ! local scalars 
    807       zcoef_dqla = -Ls * 11637800. * (-5897.8) 
    808       ! 
    809       zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 
     969      zcoef_dqlw = 4._wp * 0.95_wp * stefan             ! local scalars 
     970      zcoef_dqla = -rLsub * 11637800._wp * (-5897.8_wp) !LB: BAD! 
     971      ! 
     972      SELECT CASE( nhumi ) 
     973      CASE( np_humi_sph ) 
     974         zqair(:,:) =  phumi(:,:)      ! what we read in file is already a spec. humidity! 
     975      CASE( np_humi_dpt ) 
     976         zqair(:,:) = q_sat( phumi(:,:), pslp ) 
     977      CASE( np_humi_rlh ) 
     978         zqair(:,:) = q_air_rh( 0.01_wp*phumi(:,:), ptair(:,:), pslp(:,:) ) !LB: 0.01 => RH is % percent in file 
     979      END SELECT 
    810980      ! 
    811981      zztmp = 1. / ( 1. - albo ) 
    812       WHERE( ptsu(:,:,:) /= 0._wp )   ;   z1_st(:,:,:) = 1._wp / ptsu(:,:,:) 
    813       ELSEWHERE                       ;   z1_st(:,:,:) = 0._wp 
     982      WHERE( ptsu(:,:,:) /= 0._wp ) 
     983         z1_st(:,:,:) = 1._wp / ptsu(:,:,:) 
     984      ELSEWHERE 
     985         z1_st(:,:,:) = 0._wp 
    814986      END WHERE 
    815987      !                                     ! ========================== ! 
    … …  
    825997               qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj) 
    826998               ! Long  Wave (lw) 
    827                z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) 
     999               z_qlw(ji,jj,jl) = 0.95 * ( pqlw(ji,jj) - stefan * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) 
    8281000               ! lw sensitivity 
    8291001               z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3 
    … …  
    8331005               ! ----------------------------! 
    8341006 
    835                ! ... turbulent heat fluxes with Ch_atm recalculated in blk_ice_tau 
     1007               ! ... turbulent heat fluxes with Ch_ice recalculated in blk_ice_1 
    8361008               ! Sensible Heat 
    837                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)) 
     1009               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)) 
    8381010               ! Latent Heat 
    839                qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls  * Ch_atm(ji,jj) * wndm_ice(ji,jj) *  & 
    840                   &                ( 11637800. * EXP( -5897.8 * z1_st(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1) ) ) 
     1011               zztmp2 = EXP( -5897.8 * z1_st(ji,jj,jl) ) 
     1012               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa(ji,jj) * rLsub  * Ce_ice(ji,jj) * wndm_ice(ji,jj) *  & 
     1013                  &                ( 11637800. * zztmp2 / rhoa(ji,jj) - zqair(ji,jj) ) ) 
    8411014               ! Latent heat sensitivity for ice (Dqla/Dt) 
    8421015               IF( qla_ice(ji,jj,jl) > 0._wp ) THEN 
    843                   dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * Ch_atm(ji,jj) * wndm_ice(ji,jj) *  & 
    844                      &                 z1_st(ji,jj,jl)*z1_st(ji,jj,jl) * EXP(-5897.8 * z1_st(ji,jj,jl)) 
     1016                  dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * Ce_ice(ji,jj) * wndm_ice(ji,jj) *  & 
     1017                     &                 z1_st(ji,jj,jl) * z1_st(ji,jj,jl) * zztmp2 
    8451018               ELSE 
    8461019                  dqla_ice(ji,jj,jl) = 0._wp 
    … …  
    8481021 
    8491022               ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice) 
    850                z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_atm(ji,jj) * wndm_ice(ji,jj) 
     1023               z_dqsb(ji,jj,jl) = rhoa(ji,jj) * rCp_air * Ch_ice(ji,jj) * wndm_ice(ji,jj) 
    8511024 
    8521025               ! ----------------------------! 
    … …  
    8631036      END DO 
    8641037      ! 
    865       tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac * tmask(:,:,1)  ! total precipitation [kg/m2/s] 
    866       sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac * tmask(:,:,1)  ! solid precipitation [kg/m2/s] 
    867       CALL iom_put( 'snowpre', sprecip )                    ! Snow precipitation 
    868       CALL iom_put( 'precip' , tprecip )                    ! Total precipitation 
     1038      tprecip(:,:) = pprec(:,:) * rn_pfac * tmask(:,:,1)  ! total precipitation [kg/m2/s] 
     1039      sprecip(:,:) = psnow(:,:) * rn_pfac * tmask(:,:,1)  ! solid precipitation [kg/m2/s] 
     1040      CALL iom_put( 'snowpre', sprecip )                  ! Snow precipitation 
     1041      CALL iom_put( 'precip' , tprecip )                  ! Total precipitation 
    8691042 
    8701043      ! --- evaporation --- ! 
    … …  
    8831056      ! --- heat flux associated with emp --- ! 
    8841057      qemp_oce(:,:) = - ( 1._wp - at_i_b(:,:) ) * zevap(:,:) * sst_m(:,:) * rcp                  & ! evap at sst 
    885          &          + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp  & ! liquid precip at Tair 
     1058         &          + ( tprecip(:,:) - sprecip(:,:) ) * ( ptair(:,:) - rt0 ) * rcp               & ! liquid precip at Tair 
    8861059         &          +   sprecip(:,:) * ( 1._wp - zsnw ) *                                        & ! solid precip at min(Tair,Tsnow) 
    887          &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus ) 
     1060         &              ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus ) 
    8881061      qemp_ice(:,:) =   sprecip(:,:) * zsnw *                                                    & ! solid precip (only) 
    889          &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus ) 
     1062         &              ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus ) 
    8901063 
    8911064      ! --- total solar and non solar fluxes --- ! 
    … …  
    8951068 
    8961069      ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- ! 
    897       qprec_ice(:,:) = rhos * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus ) 
     1070      qprec_ice(:,:) = rhos * ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus ) 
    8981071 
    8991072      ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) --- 
    9001073      DO jl = 1, jpl 
    9011074         qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * rcpi * tmask(:,:,1) ) 
    902          !                         ! But we do not have Tice => consider it at 0degC => evap=0  
     1075         !                         ! But we do not have Tice => consider it at 0degC => evap=0 
    9031076      END DO 
    9041077 
    … …  
    9071080      zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )            ! zfr2 such that zfr1 + zfr2 to equal 1 
    9081081      ! 
    909       WHERE    ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) <  0.1_wp )       ! linear decrease from hi=0 to 10cm   
     1082      WHERE    ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) <  0.1_wp )       ! linear decrease from hi=0 to 10cm 
    9101083         qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) ) 
    9111084      ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp )       ! constant (zfr1) when hi>10cm 
    9121085         qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * zfr1 
    9131086      ELSEWHERE                                                         ! zero when hs>0 
    914          qtr_ice_top(:,:,:) = 0._wp  
     1087         qtr_ice_top(:,:,:) = 0._wp 
    9151088      END WHERE 
    9161089      ! 
    9171090 
    9181091      IF( iom_use('evap_ao_cea') .OR. iom_use('hflx_evap_cea') ) THEN 
    919          ztmp(:,:) = zevap(:,:) * ( 1._wp - at_i_b(:,:) )  
    920          CALL iom_put( 'evap_ao_cea'  , ztmp(:,:) * tmask(:,:,1) )   ! ice-free oce evap (cell average) 
    921          CALL iom_put( 'hflx_evap_cea', ztmp(:,:) * sst_m(:,:) * rcp * tmask(:,:,1) )   ! heat flux from evap (cell average) 
     1092         ztmp(:,:) = zevap(:,:) * ( 1._wp - at_i_b(:,:) ) 
     1093         IF( iom_use('evap_ao_cea'  ) )  CALL iom_put( 'evap_ao_cea'  , ztmp(:,:) * tmask(:,:,1) )   ! ice-free oce evap (cell average) 
     1094         IF( iom_use('hflx_evap_cea') )  CALL iom_put( 'hflx_evap_cea', ztmp(:,:) * sst_m(:,:) * rcp * tmask(:,:,1) )   ! heat flux from evap (cell average) 
    9221095      ENDIF 
    9231096      IF( iom_use('hflx_rain_cea') ) THEN 
    9241097         ztmp(:,:) = rcp * ( SUM( (ptsu-rt0) * a_i_b, dim=3 ) + sst_m(:,:) * ( 1._wp - at_i_b(:,:) ) ) 
    925          CALL iom_put( 'hflx_rain_cea', ( tprecip(:,:) - sprecip(:,:) ) * ztmp(:,:) )   ! heat flux from rain (cell average) 
     1098         IF( iom_use('hflx_rain_cea') )  CALL iom_put( 'hflx_rain_cea', ( tprecip(:,:) - sprecip(:,:) ) * ztmp(:,:) )   ! heat flux from rain (cell average) 
    9261099      ENDIF 
    9271100      IF( iom_use('hflx_snow_cea') .OR. iom_use('hflx_snow_ao_cea') .OR. iom_use('hflx_snow_ai_cea')  )  THEN 
    928           WHERE( SUM( a_i_b, dim=3 ) > 1.e-10 ) ;   ztmp(:,:) = rcpi * SUM( (ptsu-rt0) * a_i_b, dim=3 ) / SUM( a_i_b, dim=3 ) 
    929           ELSEWHERE                             ;   ztmp(:,:) = rcp * sst_m(:,:)     
    930           ENDWHERE 
    931           ztmp2(:,:) = sprecip(:,:) * ( ztmp(:,:) - rLfus )  
    932           CALL iom_put('hflx_snow_cea'   , ztmp2(:,:) ) ! heat flux from snow (cell average) 
    933           CALL iom_put('hflx_snow_ao_cea', ztmp2(:,:) * ( 1._wp - zsnw(:,:) ) ) ! heat flux from snow (over ocean) 
    934           CALL iom_put('hflx_snow_ai_cea', ztmp2(:,:) *           zsnw(:,:)   ) ! heat flux from snow (over ice) 
    935       ENDIF 
    936       ! 
    937       IF(ln_ctl) THEN 
     1101         WHERE( SUM( a_i_b, dim=3 ) > 1.e-10 ) 
     1102            ztmp(:,:) = rcpi * SUM( (ptsu-rt0) * a_i_b, dim=3 ) / SUM( a_i_b, dim=3 ) 
     1103         ELSEWHERE 
     1104            ztmp(:,:) = rcp * sst_m(:,:) 
     1105         ENDWHERE 
     1106         ztmp2(:,:) = sprecip(:,:) * ( ztmp(:,:) - rLfus ) 
     1107         IF( iom_use('hflx_snow_cea')    ) CALL iom_put('hflx_snow_cea'   , ztmp2(:,:) ) ! heat flux from snow (cell average) 
     1108         IF( iom_use('hflx_snow_ao_cea') ) CALL iom_put('hflx_snow_ao_cea', ztmp2(:,:) * ( 1._wp - zsnw(:,:) ) ) ! heat flux from snow (over ocean) 
     1109         IF( iom_use('hflx_snow_ai_cea') ) CALL iom_put('hflx_snow_ai_cea', ztmp2(:,:) *           zsnw(:,:)   ) ! heat flux from snow (over ice) 
     1110      ENDIF 
     1111      ! 
     1112      IF(sn_cfctl%l_prtctl) THEN 
    9381113         CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice: qla_ice  : ', tab3d_2=z_qsb   , clinfo2=' z_qsb    : ', kdim=jpl) 
    9391114         CALL prt_ctl(tab3d_1=z_qlw   , clinfo1=' blk_ice: z_qlw    : ', tab3d_2=dqla_ice, clinfo2=' dqla_ice : ', kdim=jpl) 
    … …  
    9441119      ENDIF 
    9451120      ! 
    946    END SUBROUTINE blk_ice_flx 
    947     
     1121   END SUBROUTINE blk_ice_2 
     1122 
    9481123 
    9491124   SUBROUTINE blk_ice_qcn( ld_virtual_itd, ptsu, ptb, phs, phi ) 
    … …  
    9541129      !!                to force sea ice / snow thermodynamics 
    9551130      !!                in the case conduction flux is emulated 
    956       !!                 
     1131      !! 
    9571132      !! ** Method  :   compute surface energy balance assuming neglecting heat storage 
    9581133      !!                following the 0-layer Semtner (1976) approach 
    … …  
    9791154      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zgfac   ! enhanced conduction factor 
    9801155      !!--------------------------------------------------------------------- 
    981        
     1156 
    9821157      ! -------------------------------------! 
    9831158      !      I   Enhanced conduction factor  ! 
    … …  
    9871162      ! 
    9881163      zgfac(:,:,:) = 1._wp 
    989        
     1164 
    9901165      IF( ld_virtual_itd ) THEN 
    9911166         ! 
    … …  
    9931168         zfac2 = EXP(1._wp) * 0.5_wp * zepsilon 
    9941169         zfac3 = 2._wp / zepsilon 
    995          !    
    996          DO jl = 1, jpl                 
    997             DO jj = 1 , jpj 
    998                DO ji = 1, jpi 
    999                   zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcnd_i * phs(ji,jj,jl) ) * zfac                            ! Effective thickness 
    1000                   IF( zhe >=  zfac2 )   zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) ) ! Enhanced conduction factor 
    1001                END DO 
    1002             END DO 
     1170         ! 
     1171         DO jl = 1, jpl 
     1172            DO_2D_11_11 
     1173               zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcnd_i * phs(ji,jj,jl) ) * zfac                            ! Effective thickness 
     1174               IF( zhe >=  zfac2 )   zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) ) ! Enhanced conduction factor 
     1175            END_2D 
    10031176         END DO 
    1004          !       
    1005       ENDIF 
    1006        
     1177         ! 
     1178      ENDIF 
     1179 
    10071180      ! -------------------------------------------------------------! 
    10081181      !      II   Surface temperature and conduction flux            ! 
    … …  
    10121185      ! 
    10131186      DO jl = 1, jpl 
    1014          DO jj = 1 , jpj 
    1015             DO ji = 1, jpi 
    1016                !                     
    1017                zkeff_h = zfac * zgfac(ji,jj,jl) / &                                    ! Effective conductivity of the snow-ice system divided by thickness 
    1018                   &      ( rcnd_i * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) ) 
    1019                ztsu    = ptsu(ji,jj,jl)                                                ! Store current iteration temperature 
    1020                ztsu0   = ptsu(ji,jj,jl)                                                ! Store initial surface temperature 
    1021                zqa0    = qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) ! Net initial atmospheric heat flux 
    1022                ! 
    1023                DO iter = 1, nit     ! --- Iterative loop 
    1024                   zqc   = zkeff_h * ( ztsu - ptb(ji,jj) )                              ! Conduction heat flux through snow-ice system (>0 downwards) 
    1025                   zqnet = zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc  ! Surface energy budget 
    1026                   ztsu  = ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h )              ! Temperature update 
    1027                END DO 
    1028                ! 
    1029                ptsu   (ji,jj,jl) = MIN( rt0, ztsu ) 
    1030                qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) ) 
    1031                qns_ice(ji,jj,jl) = qns_ice(ji,jj,jl) + dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) 
    1032                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) )  & 
    1033                              &   * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) ) 
    1034  
    1035                ! --- Diagnose the heat loss due to changing non-solar flux (as in icethd_zdf_bl99) --- ! 
    1036                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)  
    1037  
     1187         DO_2D_11_11 
     1188            ! 
     1189            zkeff_h = zfac * zgfac(ji,jj,jl) / &                                    ! Effective conductivity of the snow-ice system divided by thickness 
     1190               &      ( rcnd_i * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) ) 
     1191            ztsu    = ptsu(ji,jj,jl)                                                ! Store current iteration temperature 
     1192            ztsu0   = ptsu(ji,jj,jl)                                                ! Store initial surface temperature 
     1193            zqa0    = qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) ! Net initial atmospheric heat flux 
     1194            ! 
     1195            DO iter = 1, nit     ! --- Iterative loop 
     1196               zqc   = zkeff_h * ( ztsu - ptb(ji,jj) )                              ! Conduction heat flux through snow-ice system (>0 downwards) 
     1197               zqnet = zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc  ! Surface energy budget 
     1198               ztsu  = ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h )              ! Temperature update 
    10381199            END DO 
    1039          END DO 
     1200            ! 
     1201            ptsu   (ji,jj,jl) = MIN( rt0, ztsu ) 
     1202            qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) ) 
     1203            qns_ice(ji,jj,jl) = qns_ice(ji,jj,jl) + dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) 
     1204            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) )  & 
     1205               &   * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) ) 
     1206 
     1207            ! --- Diagnose the heat loss due to changing non-solar flux (as in icethd_zdf_bl99) --- ! 
     1208            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) 
     1209 
     1210         END_2D 
    10401211         ! 
    1041       END DO  
    1042       !       
     1212      END DO 
     1213      ! 
    10431214   END SUBROUTINE blk_ice_qcn 
    1044     
    1045  
    1046    SUBROUTINE Cdn10_Lupkes2012( Cd ) 
     1215 
     1216 
     1217   SUBROUTINE Cdn10_Lupkes2012( pcd ) 
    10471218      !!---------------------------------------------------------------------- 
    10481219      !!                      ***  ROUTINE  Cdn10_Lupkes2012  *** 
    10491220      !! 
    1050       !! ** Purpose :    Recompute the neutral air-ice drag referenced at 10m  
     1221      !! ** Purpose :    Recompute the neutral air-ice drag referenced at 10m 
    10511222      !!                 to make it dependent on edges at leads, melt ponds and flows. 
    10521223      !!                 After some approximations, this can be resumed to a dependency 
    10531224      !!                 on ice concentration. 
    1054       !!                 
     1225      !! 
    10551226      !! ** Method :     The parameterization is taken from Lupkes et al. (2012) eq.(50) 
    10561227      !!                 with the highest level of approximation: level4, eq.(59) 
    … …  
    10641235      !! 
    10651236      !!                 This new drag has a parabolic shape (as a function of A) starting at 
    1066       !!                 Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5  
     1237      !!                 Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5 
    10671238      !!                 and going down to Cdi(say 1.4e-3) for A=1 
    10681239      !! 
    … …  
    10741245      !! 
    10751246      !!---------------------------------------------------------------------- 
    1076       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Cd 
     1247      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pcd 
    10771248      REAL(wp), PARAMETER ::   zCe   = 2.23e-03_wp 
    10781249      REAL(wp), PARAMETER ::   znu   = 1._wp 
    … …  
    10891260 
    10901261      ! ice-atm drag 
    1091       Cd(:,:) = Cd_ice +  &                                                         ! pure ice drag 
    1092          &      zCe    * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp)  ! change due to sea-ice morphology 
    1093        
     1262      pcd(:,:) = rCd_ice +  &                                                         ! pure ice drag 
     1263         &      zCe     * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp)  ! change due to sea-ice morphology 
     1264 
    10941265   END SUBROUTINE Cdn10_Lupkes2012 
    10951266 
    10961267 
    1097    SUBROUTINE Cdn10_Lupkes2015( Cd, Ch ) 
     1268   SUBROUTINE Cdn10_Lupkes2015( ptm_su, pslp, pcd, pch ) 
    10981269      !!---------------------------------------------------------------------- 
    10991270      !!                      ***  ROUTINE  Cdn10_Lupkes2015  *** 
    11001271      !! 
    11011272      !! ** pUrpose :    Alternative turbulent transfert coefficients formulation 
    1102       !!                 between sea-ice and atmosphere with distinct momentum  
    1103       !!                 and heat coefficients depending on sea-ice concentration  
     1273      !!                 between sea-ice and atmosphere with distinct momentum 
     1274      !!                 and heat coefficients depending on sea-ice concentration 
    11041275      !!                 and atmospheric stability (no meltponds effect for now). 
    1105       !!                 
     1276      !! 
    11061277      !! ** Method :     The parameterization is adapted from Lupkes et al. (2015) 
    11071278      !!                 and ECHAM6 atmospheric model. Compared to Lupkes2012 scheme, 
    11081279      !!                 it considers specific skin and form drags (Andreas et al. 2010) 
    1109       !!                 to compute neutral transfert coefficients for both heat and  
     1280      !!                 to compute neutral transfert coefficients for both heat and 
    11101281      !!                 momemtum fluxes. Atmospheric stability effect on transfert 
    11111282      !!                 coefficient is also taken into account following Louis (1979). 
    … …  
    11161287      !!---------------------------------------------------------------------- 
    11171288      ! 
    1118       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Cd 
    1119       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Ch 
    1120       REAL(wp), DIMENSION(jpi,jpj)            ::   ztm_su, zst, zqo_sat, zqi_sat 
     1289      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   ptm_su ! sea-ice surface temperature [K] 
     1290      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pslp   ! sea-level pressure [Pa] 
     1291      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pcd    ! momentum transfert coefficient 
     1292      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pch    ! heat transfert coefficient 
     1293      REAL(wp), DIMENSION(jpi,jpj)            ::   zst, zqo_sat, zqi_sat 
    11211294      ! 
    11221295      ! ECHAM6 constants 
    … …  
    11461319      !!---------------------------------------------------------------------- 
    11471320 
    1148       ! mean temperature 
    1149       WHERE( at_i_b(:,:) > 1.e-20 )   ;   ztm_su(:,:) = SUM( t_su(:,:,:) * a_i_b(:,:,:) , dim=3 ) / at_i_b(:,:) 
    1150       ELSEWHERE                       ;   ztm_su(:,:) = rt0 
    1151       ENDWHERE 
    1152        
    11531321      ! Momentum Neutral Transfert Coefficients (should be a constant) 
    11541322      zCdn_form_tmp = zce10 * ( LOG( 10._wp / z0_form_ice + 1._wp ) / LOG( rn_zu / z0_form_ice + 1._wp ) )**2   ! Eq. 40 
    11551323      zCdn_skin_ice = ( vkarmn                                      / LOG( rn_zu / z0_skin_ice + 1._wp ) )**2   ! Eq. 7 
    1156       zCdn_ice      = zCdn_skin_ice   ! Eq. 7 (cf Lupkes email for details) 
     1324      zCdn_ice      = zCdn_skin_ice   ! Eq. 7 
    11571325      !zCdn_ice     = 1.89e-3         ! old ECHAM5 value (cf Eq. 32) 
    11581326 
    11591327      ! Heat Neutral Transfert Coefficients 
    1160       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) 
    1161       
     1328      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 
     1329 
    11621330      ! Atmospheric and Surface Variables 
    11631331      zst(:,:)     = sst_m(:,:) + rt0                                        ! convert SST from Celcius to Kelvin 
    1164       zqo_sat(:,:) = 0.98_wp * q_sat( zst(:,:)   , sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ocean [kg/kg] 
    1165       zqi_sat(:,:) = 0.98_wp * q_sat( ztm_su(:,:), sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ice   [kg/kg] 
    1166       ! 
    1167       DO jj = 2, jpjm1           ! reduced loop is necessary for reproducibility 
    1168          DO ji = fs_2, fs_jpim1 
    1169             ! Virtual potential temperature [K] 
    1170             zthetav_os = zst(ji,jj)    * ( 1._wp + rctv0 * zqo_sat(ji,jj) )   ! over ocean 
    1171             zthetav_is = ztm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) )   ! ocean ice 
    1172             zthetav_zu = t_zu (ji,jj)  * ( 1._wp + rctv0 * q_zu(ji,jj)    )   ! at zu 
    1173              
    1174             ! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead) 
    1175             zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj)     )**2   ! over ocean 
    1176             zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2   ! over ice 
    1177              
    1178             ! Momentum and Heat Neutral Transfert Coefficients 
    1179             zCdn_form_ice = zCdn_form_tmp * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta  ! Eq. 40 
    1180             zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) )               ! Eq. 53  
    1181                         
    1182             ! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead) 
    1183             z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water 
    1184             z0i = z0_skin_ice                                             ! over ice (cf Lupkes email for details) 
    1185             IF( zrib_o <= 0._wp ) THEN 
    1186                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 
    1187                zfhw = ( 1._wp + ( zbetah * ( zthetav_os - zthetav_zu )**r1_3 / ( Chn_oce(ji,jj) * MAX(0.01, wndm(ji,jj)) )   &     ! Eq. 26 
    1188                   &             )**zgamma )**z1_gamma 
    1189             ELSE 
    1190                zfmw = 1._wp / ( 1._wp + zam * zrib_o / SQRT( 1._wp + zrib_o ) )   ! Eq. 12 
    1191                zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) )   ! Eq. 28 
    1192             ENDIF 
    1193              
    1194             IF( zrib_i <= 0._wp ) THEN 
    1195                zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp)))   ! Eq.  9 
    1196                zfhi = 1._wp - zah * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp)))   ! Eq. 25 
    1197             ELSE 
    1198                zfmi = 1._wp / ( 1._wp + zam * zrib_i / SQRT( 1._wp + zrib_i ) )   ! Eq. 11 
    1199                zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) )   ! Eq. 27 
    1200             ENDIF 
    1201              
    1202             ! Momentum Transfert Coefficients (Eq. 38) 
    1203             Cd(ji,jj) = zCdn_skin_ice *   zfmi +  & 
    1204                &        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) ) 
    1205              
    1206             ! Heat Transfert Coefficients (Eq. 49) 
    1207             Ch(ji,jj) = zChn_skin_ice *   zfhi +  & 
    1208                &        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) ) 
    1209             ! 
    1210          END DO 
    1211       END DO 
    1212       CALL lbc_lnk_multi( 'sbcblk', Cd, 'T',  1., Ch, 'T', 1. ) 
     1332      zqo_sat(:,:) = rdct_qsat_salt * q_sat( zst(:,:)   , pslp(:,:) )   ! saturation humidity over ocean [kg/kg] 
     1333      zqi_sat(:,:) =                  q_sat( ptm_su(:,:), pslp(:,:) )   ! saturation humidity over ice   [kg/kg] 
     1334      ! 
     1335      DO_2D_00_00 
     1336         ! Virtual potential temperature [K] 
     1337         zthetav_os = zst(ji,jj)    * ( 1._wp + rctv0 * zqo_sat(ji,jj) )   ! over ocean 
     1338         zthetav_is = ptm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) )   ! ocean ice 
     1339         zthetav_zu = t_zu (ji,jj)  * ( 1._wp + rctv0 * q_zu(ji,jj)    )   ! at zu 
     1340 
     1341         ! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead) 
     1342         zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj)     )**2   ! over ocean 
     1343         zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2   ! over ice 
     1344 
     1345         ! Momentum and Heat Neutral Transfert Coefficients 
     1346         zCdn_form_ice = zCdn_form_tmp * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta  ! Eq. 40 
     1347         zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) )               ! Eq. 53 
     1348 
     1349         ! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead ?) 
     1350         z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water 
     1351         z0i = z0_skin_ice                                             ! over ice 
     1352         IF( zrib_o <= 0._wp ) THEN 
     1353            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 
     1354            zfhw = ( 1._wp + ( zbetah * ( zthetav_os - zthetav_zu )**r1_3 / ( Chn_oce(ji,jj) * MAX(0.01, wndm(ji,jj)) )   &     ! Eq. 26 
     1355               &             )**zgamma )**z1_gamma 
     1356         ELSE 
     1357            zfmw = 1._wp / ( 1._wp + zam * zrib_o / SQRT( 1._wp + zrib_o ) )   ! Eq. 12 
     1358            zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) )   ! Eq. 28 
     1359         ENDIF 
     1360 
     1361         IF( zrib_i <= 0._wp ) THEN 
     1362            zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp)))   ! Eq.  9 
     1363            zfhi = 1._wp - zah * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp)))   ! Eq. 25 
     1364         ELSE 
     1365            zfmi = 1._wp / ( 1._wp + zam * zrib_i / SQRT( 1._wp + zrib_i ) )   ! Eq. 11 
     1366            zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) )   ! Eq. 27 
     1367         ENDIF 
     1368 
     1369         ! Momentum Transfert Coefficients (Eq. 38) 
     1370         pcd(ji,jj) = zCdn_skin_ice *   zfmi +  & 
     1371            &        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) ) 
     1372 
     1373         ! Heat Transfert Coefficients (Eq. 49) 
     1374         pch(ji,jj) = zChn_skin_ice *   zfhi +  & 
     1375            &        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) ) 
     1376         ! 
     1377      END_2D 
     1378      CALL lbc_lnk_multi( 'sbcblk', pcd, 'T',  1., pch, 'T', 1. ) 
    12131379      ! 
    12141380   END SUBROUTINE Cdn10_Lupkes2015 
Note: See TracChangeset for help on using the changeset viewer.