Index: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbcblk.F90
===================================================================
--- NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbcblk.F90	(revision 11308)
+++ NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbcblk.F90	(revision 11309)
@@ -77,16 +77,47 @@
 #endif 
 
-  INTERFACE cp_air
-    MODULE PROCEDURE cp_air_0d, cp_air_2d
-  END INTERFACE
-
-  !!Lolo: should ultimately be moved in the module with all physical constants ?
-!!gm  : In principle, yes.
-   REAL(wp)        , PARAMETER ::   Cp_dry = 1005.0       !: Specic heat of dry air, constant pressure      [J/K/kg]
-   REAL(wp)        , PARAMETER ::   Cp_vap = 1860.0       !: Specic heat of water vapor, constant pressure  [J/K/kg]
-   REAL(wp), PUBLIC, PARAMETER ::   R_dry = 287.05_wp     !: Specific gas constant for dry air              [J/K/kg]
-   REAL(wp)        , PARAMETER ::   R_vap = 461.495_wp    !: Specific gas constant for water vapor          [J/K/kg]
-   REAL(wp)        , PARAMETER ::   reps0 = R_dry/R_vap   !: ratio of gas constant for dry air and water vapor => ~ 0.622
-   REAL(wp), PUBLIC, PARAMETER ::   rctv0 = R_vap/R_dry   !: for virtual temperature (== (1-eps)/eps) => ~ 0.608
+   INTERFACE cp_air
+      MODULE PROCEDURE cp_air_0d, cp_air_2d
+   END INTERFACE
+
+   !                                   !!* Namelist namsbc_blk : bulk parameters
+   LOGICAL          ::   ln_NCAR        ! "NCAR"      algorithm   (Large and Yeager 2008)
+   LOGICAL          ::   ln_COARE_3p0   ! "COARE 3.0" algorithm   (Fairall et al. 2003)
+   LOGICAL          ::   ln_COARE_3p5   ! "COARE 3.5" algorithm   (Edson et al. 2013)
+   LOGICAL          ::   ln_ECMWF       ! "ECMWF"     algorithm   (IFS cycle 31)
+   !                                    !
+   REAL(wp)         ::   rn_pfac        !  multiplication factor for precipitation
+   REAL(wp), PUBLIC ::   rn_efac        !: multiplication factor for evaporation
+   REAL(wp), PUBLIC ::   rn_vfac        !: multiplication factor for ice/ocean velocity in the calculation of wind stress
+   REAL(wp)         ::   rn_zqt         !  z(q,t) : height of humidity and temperature measurements
+   REAL(wp)         ::   rn_zu          !  z(u)   : height of wind measurements
+   !                                    !
+   LOGICAL          ::   ln_Cd_L12      ! ice-atm drag = F( ice concentration )                        (Lupkes et al. JGR2012)
+   LOGICAL          ::   ln_Cd_L15      ! ice-atm drag = F( ice concentration, atmospheric stability ) (Lupkes et al. JGR2015)
+
+   INTEGER  ::   nblk                   ! choice of the bulk algorithm
+   !                                       ! associated indices:
+   INTEGER, PARAMETER ::   np_NCAR      = 1   ! "NCAR" algorithm        (Large and Yeager 2008)
+   INTEGER, PARAMETER ::   np_COARE_3p0 = 2   ! "COARE 3.0" algorithm   (Fairall et al. 2003)
+   INTEGER, PARAMETER ::   np_COARE_3p5 = 3   ! "COARE 3.5" algorithm   (Edson et al. 2013)
+   INTEGER, PARAMETER ::   np_ECMWF     = 4   ! "ECMWF" algorithm       (IFS cycle 31)
+
+   !                                                      !!! air parameters
+   REAL(wp)        , PARAMETER ::   Cp_dry = 1005.0        !  Specic heat of dry air, constant pressure      [J/K/kg]
+   REAL(wp)        , PARAMETER ::   Cp_vap = 1860.0        !  Specic heat of water vapor, constant pressure  [J/K/kg]
+   REAL(wp), PUBLIC, PARAMETER ::   R_dry  = 287.05_wp     !: Specific gas constant for dry air              [J/K/kg]
+   REAL(wp)        , PARAMETER ::   R_vap  = 461.495_wp    !  Specific gas constant for water vapor          [J/K/kg]
+   REAL(wp)        , PARAMETER ::   reps0  = R_dry/R_vap   !  ratio of gas constant for dry air and water vapor => ~ 0.622
+   REAL(wp), PUBLIC, PARAMETER ::   rctv0  = R_vap/R_dry   !: for virtual temperature (== (1-eps)/eps) => ~ 0.608
+   !                                                      !!! Bulk parameters
+   REAL(wp)        , PARAMETER ::   cpa    = 1000.5        ! specific heat of air (only used for ice fluxes now...)
+   REAL(wp)        , PARAMETER ::   Ls     =    2.839e6    ! latent heat of sublimation
+   REAL(wp)        , PARAMETER ::   Stef   =    5.67e-8    ! Stefan Boltzmann constant
+   REAL(wp)        , PARAMETER ::   Cd_ice =    1.4e-3     ! transfer coefficient over ice
+   REAL(wp)        , PARAMETER ::   albo   =    0.066      ! ocean albedo assumed to be constant
+   !
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Cd_atm , Ch_atm , Ce_atm   ! transfer coeffs. for momentum, sensible heat, and evaporation
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Cdn_oce, Chn_oce, Cen_oce  ! neutral  coeffs (needed for Lupkes et al. 2015 bulk scheme)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_zu, q_zu                 ! air temp. and spec. hum. at wind speed height (needed for Lupkes 2015 bulk scheme)
 
    INTEGER , PUBLIC, PARAMETER ::   jpfld   =11   !: maximum number of files to read
@@ -103,41 +134,5 @@
    INTEGER , PUBLIC, PARAMETER ::   jp_hpgj =11   !: index of ABL geostrophic wind or hpg (j-component) (m/s) at T-point
 
-   TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) ::   sf   !  structure of input fields (file informations, fields read)
-
-   !                                             !!! Bulk parameters
-   REAL(wp), PARAMETER ::   cpa    = 1000.5         ! specific heat of air (only used for ice fluxes now...)
-   REAL(wp), PARAMETER ::   Ls     =    2.839e6     ! latent heat of sublimation
-   REAL(wp), PARAMETER ::   Stef   =    5.67e-8     ! Stefan Boltzmann constant
-   REAL(wp), PARAMETER ::   Cd_ice =    1.4e-3      ! transfer coefficient over ice
-   REAL(wp), PARAMETER ::   albo   =    0.066       ! ocean albedo assumed to be constant
-   !
-   !                           !!* Namelist namsbc_blk : bulk parameters
-   LOGICAL  ::   ln_NCAR        ! "NCAR"      algorithm   (Large and Yeager 2008)
-   LOGICAL  ::   ln_COARE_3p0   ! "COARE 3.0" algorithm   (Fairall et al. 2003)
-   LOGICAL  ::   ln_COARE_3p5   ! "COARE 3.5" algorithm   (Edson et al. 2013)
-   LOGICAL  ::   ln_ECMWF       ! "ECMWF"     algorithm   (IFS cycle 31)
-   !
-   REAL(wp)         ::   rn_pfac        ! multiplication factor for precipitation
-   REAL(wp), PUBLIC ::   rn_efac        !: multiplication factor for evaporation
-   REAL(wp), PUBLIC ::   rn_vfac        !: multiplication factor for ice/ocean velocity in the calculation of wind stress
-   REAL(wp) ::   rn_zqt         ! z(q,t) : height of humidity and temperature measurements
-   REAL(wp) ::   rn_zu          ! z(u)   : height of wind measurements
-!!gm ref namelist initialize it so remove the setting to false below
-   LOGICAL  ::   ln_Cd_L12 = .FALSE. !  Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2012)
-   LOGICAL  ::   ln_Cd_L15 = .FALSE. !  Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2015)
-   !
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Cd_atm                    ! transfer coefficient for momentum      (tau)
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Ch_atm                    ! transfer coefficient for sensible heat (Q_sens)
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Ce_atm                    ! tansfert coefficient for evaporation   (Q_lat)
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_zu                      ! air temperature at wind speed height (needed by Lupkes 2015 bulk scheme)
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   q_zu                      ! air spec. hum.  at wind speed height (needed by Lupkes 2015 bulk scheme)
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   cdn_oce, chn_oce, cen_oce ! needed by Lupkes 2015 bulk scheme
-
-   INTEGER  ::   nblk           ! choice of the bulk algorithm
-   !                            ! associated indices:
-   INTEGER, PARAMETER ::   np_NCAR      = 1   ! "NCAR" algorithm        (Large and Yeager 2008)
-   INTEGER, PARAMETER ::   np_COARE_3p0 = 2   ! "COARE 3.0" algorithm   (Fairall et al. 2003)
-   INTEGER, PARAMETER ::   np_COARE_3p5 = 3   ! "COARE 3.5" algorithm   (Edson et al. 2013)
-   INTEGER, PARAMETER ::   np_ECMWF     = 4   ! "ECMWF" algorithm       (IFS cycle 31)
+   TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) ::   sf   !  structure of input atmospheric fields (file informations, fields read)
 
    !! * Substitutions
@@ -155,5 +150,5 @@
       !!-------------------------------------------------------------------
       ALLOCATE( Cd_atm (jpi,jpj), Ch_atm (jpi,jpj), Ce_atm (jpi,jpj), t_zu(jpi,jpj), q_zu(jpi,jpj), &
-         &      cdn_oce(jpi,jpj), chn_oce(jpi,jpj), cen_oce(jpi,jpj), STAT=sbc_blk_alloc )
+         &      Cdn_oce(jpi,jpj), Chn_oce(jpi,jpj), Cen_oce(jpi,jpj), STAT=sbc_blk_alloc )
       !
       CALL mpp_sum ( 'sbcblk', sbc_blk_alloc )
@@ -171,5 +166,5 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER  ::   ifpr                  ! dummy loop indice and argument
+      INTEGER  ::   jfpr                  ! dummy loop indice and argument
       INTEGER  ::   ios, ierror, ioptio   ! Local integer
       !!
@@ -233,17 +228,21 @@
       CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_init', 'surface boundary condition -- bulk formulae', 'namsbc_blk' )
       !
-      DO ifpr= 1, jpfld
+      DO jfpr= 1, jpfld
          !
-         IF( TRIM(sf(ifpr)%clrootname) /= 'NOT USED' ) THEN
-            IF( ln_abl .AND. (     ifpr == jp_wndi .OR. ifpr == jp_wndj .OR. ifpr == jp_humi   &
-               &              .OR. ifpr == jp_hpgi .OR. ifpr == jp_hpgj .OR. ifpr == jp_tair ) ) THEN
-               ALLOCATE( sf(ifpr)%fnow(jpi,jpj,jpka) )
-               IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf(ifpr)%fdta(jpi,jpj,jpka,2) )
-            ELSE
-               ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) )
-               IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )
+         IF( TRIM(sf(jfpr)%clrootname) == 'NOT USED' ) THEN    !--  not used field  --!   (only now allocated and set to zero) 
+            ALLOCATE( sf(jfpr)%fnow(jpi,jpj,1) )
+            sf(jfpr)%fnow(:,:,1) = 0._wp                 
+         ELSE                                                  !-- used field  --!
+            IF(   ln_abl    .AND.                                                      &
+               &    ( jfpr == jp_wndi .OR. jfpr == jp_wndj .OR. jfpr == jp_humi .OR.   &
+               &      jfpr == jp_hpgi .OR. jfpr == jp_hpgj .OR. jfpr == jp_tair     )  ) THEN   ! ABL: some fields are 3D input
+               ALLOCATE( sf(jfpr)%fnow(jpi,jpj,jpka) )
+               IF( slf_i(jfpr)%ln_tint )   ALLOCATE( sf(jfpr)%fdta(jpi,jpj,jpka,2) )
+            ELSE                                                                                ! others or Bulk fields are 2D fiels
+               ALLOCATE( sf(jfpr)%fnow(jpi,jpj,1) )
+               IF( slf_i(jfpr)%ln_tint )   ALLOCATE( sf(jfpr)%fdta(jpi,jpj,1,2) )
             ENDIF
             !
-            IF( slf_i(ifpr)%freqh > 0. .AND. MOD( NINT(3600. * slf_i(ifpr)%freqh), nn_fsbc * NINT(rdt) ) /= 0 )   &
+            IF( slf_i(jfpr)%freqh > 0. .AND. MOD( NINT(3600. * slf_i(jfpr)%freqh), nn_fsbc * NINT(rdt) ) /= 0 )   &
                &  CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rdt*nn_fsbc is NOT a submultiple of atmospheric forcing frequency.',   &
                &                 '               This is not ideal. You should consider changing either rdt or nn_fsbc value...' )
@@ -251,28 +250,23 @@
       END DO
       !
-      IF ( ln_wave ) THEN
-      !Activated wave module but neither drag nor stokes drift activated
-         IF ( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) )   THEN
-            CALL ctl_stop( 'STOP',  'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauwoc=F, ln_stcor=F' )
-      !drag coefficient read from wave model definable only with mfs bulk formulae and core 
-         ELSEIF (ln_cdgw .AND. .NOT. ln_NCAR )       THEN       
-             CALL ctl_stop( 'drag coefficient read from wave model definable only with NCAR and CORE bulk formulae')
-         ELSEIF (ln_stcor .AND. .NOT. ln_sdw)                             THEN
-             CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T')
+      IF( ln_wave ) THEN      ! surface waves
+         IF( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) )   &   ! Activated wave module but neither drag nor stokes drift activated
+            &   CALL ctl_stop( 'sbc_blk_init: Ask for wave coupling but ln_cdgw=ln_sdw=ln_tauwoc=ln_stcor=F' )
+         IF( ln_cdgw .AND. .NOT.ln_NCAR )                                 &   ! drag coefficient read from wave model only with NCAR bulk formulae 
+            &   CALL ctl_stop( 'sbc_blk_init: drag coefficient read from wave model need NCAR bulk formulae')
+         IF( ln_stcor .AND. .NOT.ln_sdw )                                 &
+            CALL ctl_stop( 'sbc_blk_init: Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T')
+      ELSE
+         IF( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) THEN
+            CALL ctl_warn( 'sbc_blk_init: ln_wave=F, set all wave-related namelist parameter to FALSE')
+            ln_cdgw =.FALSE.   ;   ln_sdw =.FALSE.   ;   ln_tauwoc =.FALSE.   ;   ln_stcor =.FALSE.   
          ENDIF
-      ELSE
-      IF ( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor )                & 
-         &   CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ',    &
-         &                  'with drag coefficient (ln_cdgw =T) '  ,                        &
-         &                  'or Stokes Drift (ln_sdw=T) ' ,                                 &
-         &                  'or ocean stress modification due to waves (ln_tauwoc=T) ',      &  
-         &                  'or Stokes-Coriolis term (ln_stcori=T)'  )
       ENDIF 
       !
-      IF( ln_abl ) THEN      ! ABL: read 3D fields for wind, temperature and humidity     
+      IF( ln_abl ) THEN       ! ABL: read 3D fields for wind, temperature, humidity and pressure gradient    
          rn_zqt = ght_abl(2)          ! set the bulk altitude to ABL first level
          rn_zu  = ght_abl(2)
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) '   ABL formulation: overwrite rn_zqt & rn_zu with ABL first level altitude'
+         IF(lwp) WRITE(numout,*) '  ABL formulation:   overwrite rn_zqt & rn_zu with ABL first level altitude'
       ENDIF
       !           
@@ -348,5 +342,5 @@
       !
       !                                            ! compute the surface ocean fluxes using bulk formulea
-      IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
+      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
          !
          CALL blk_oce_1( kt, sf(jp_wndi)%fnow(:,:,1), sf(jp_wndj)%fnow(:,:,1),   &   !   <<= in
@@ -359,9 +353,8 @@
             &                sf(jp_snow)%fnow(:,:,1), sst_m,                     &   !   <<= in
             &                zsen, zevp )                                            !   <=> in out
-         
       ENDIF
          
 #if defined key_cice
-      IF( MOD( kt - 1, nn_fsbc ) == 0 )   THEN
+      IF( MOD( kt-1, nn_fsbc ) == 0 )   THEN
          qlw_ice(:,:,1)   = sf(jp_qlw )%fnow(:,:,1)
          IF( ln_dm2dc ) THEN ; qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) )
@@ -468,11 +461,11 @@
       !
       CASE( np_NCAR      )   ;   CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, pssq, phumi, wndm,   &  ! NCAR-COREv2
-         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
+         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, Cdn_oce, Chn_oce, Cen_oce )
       CASE( np_COARE_3p0 )   ;   CALL turb_coare   ( rn_zqt, rn_zu, zst, ztpot, pssq, phumi, wndm,   &  ! COARE v3.0
-         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
+         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, Cdn_oce, Chn_oce, Cen_oce )
       CASE( np_COARE_3p5 )   ;   CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, pssq, phumi, wndm,   &  ! COARE v3.5
-         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
+         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, Cdn_oce, Chn_oce, Cen_oce )
       CASE( np_ECMWF     )   ;   CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, pssq, phumi, wndm,   &  ! ECMWF
-         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
+         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, Cdn_oce, Chn_oce, Cen_oce )
       CASE DEFAULT
          CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' )
