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

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

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/SBC/sbcblk.F90

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