Changeset 11275


Ignore:
Timestamp:
2019-07-17T15:05:13+02:00 (13 months ago)
Author:
smasson
Message:

dev_r11265_ABL : add ABL interface, see #2131

Location:
NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D
Files:
1 added
13 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/AGRIF_DEMO/EXPREF/namelist_cfg

    r11267 r11275  
    108108   sn_snow     = 'ncar_precip.15JUNE2009_fill',   -1.        , 'SNOW'    ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , '' 
    109109   sn_slp      = 'slp.15JUNE2009_fill'        ,    6.        , 'SLP'     ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , '' 
    110    sn_tdif     = 'taudif_core'                ,   24.        , 'taudif'  ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , '' 
    111110/ 
    112111!----------------------------------------------------------------------- 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/C1D_PAPA/EXPREF/namelist_cfg

    r11267 r11275  
    159159   sn_snow     = 'forcing_C1D_PAPA'   ,  3.   , 'sososnow',   .false.    , .false. , 'yearly'  , ''  , ''   , '' 
    160160   sn_slp      = 'forcing_C1D_PAPA'   ,  3.   , 'somslpre',   .true.     , .false. , 'yearly'  , ''  , ''   , '' 
    161    sn_tdif     = 'forcing_C1D_PAPA'   , 24.   , 'taudif'  ,   .false.    , .false. , 'yearly'  , ''  , ''   , '' 
    162  
    163161/ 
    164162!----------------------------------------------------------------------- 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/ORCA2_ICE_PISCES/EXPREF/context_nemo.xml

    r10255 r11275  
    3232      <axis id="iax_20C" long_name="20 degC isotherm"   unit="degC"            /> 
    3333      <axis id="iax_28C" long_name="28 degC isotherm"   unit="degC"            /> 
     34      <!-- ABL vertical axis definition --> 
     35      <axis id="ght_abl" long_name="ABL Vertical T levels" unit="m" positive="up"   /> 
     36      <axis id="ghw_abl" long_name="ABL Vertical W levels" unit="m" positive="up"   /> 
    3437    </axis_definition> 
    3538  
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg

    r11267 r11275  
    118118   sn_snow     = 'ncar_precip.15JUNE2009_fill',   -1.        , 'SNOW'    ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , '' 
    119119   sn_slp      = 'slp.15JUNE2009_fill'        ,    6.        , 'SLP'     ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , '' 
    120    sn_tdif     = 'taudif_core'                ,   24.        , 'taudif'  ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , '' 
    121120/ 
    122121!----------------------------------------------------------------------- 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/ORCA2_SAS_ICE/EXPREF/namelist_cfg

    r11267 r11275  
    7979   sn_snow     = 'ncar_precip.15JUNE2009_fill',   -1.        , 'SNOW'    ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , '' 
    8080   sn_slp      = 'slp.15JUNE2009_fill'        ,    6.        , 'SLP'     ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , '' 
    81    sn_tdif     = 'taudif_core'                ,   24.        , 'taudif'  ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , '' 
    8281/ 
    8382!----------------------------------------------------------------------- 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/SHARED/grid_def_nemo.xml

    r10226 r11275  
    5858       </grid> 
    5959 
     60       <!-- ABL grid definition --> 
     61       <grid id="grid_TA_3D"> 
     62         <domain domain_ref="grid_T" /> 
     63         <axis  id="ght_abl" /> 
     64       </grid> 
     65       <grid id="grid_TA_2D"> 
     66         <domain domain_ref="grid_T" /> 
     67       </grid> 
     68       <grid id="grid_WA_3D"> 
     69         <domain domain_ref="grid_T" /> 
     70         <axis  id="ghw_abl" /> 
     71       </grid> 
     72       <!--  --> 
     73 
    6074    </grid_definition> 
    6175     
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/SHARED/namelist_ref

    r11268 r11275  
    195195   ln_flx      = .false.   !  flux formulation                          (T => fill namsbc_flx ) 
    196196   ln_blk      = .false.   !  Bulk formulation                          (T => fill namsbc_blk ) 
     197   ln_abl      = .false.   !  ABL  formulation                          (T => fill namsbc_abl ) 
    197198      !              ! Type of coupling (Ocean/Ice/Atmosphere) : 
    198199   ln_cpl      = .false.   !  atmosphere coupled   formulation          ( requires key_oasis3 ) 
     
    258259      ln_Cd_L12   = .false.   !  air-ice drags = F(ice concentration) (Lupkes et al. 2012) 
    259260      ln_Cd_L15   = .false.   !  air-ice drags = F(ice concentration) (Lupkes et al. 2015) 
    260       ln_taudif   = .false.   !  HF tau contribution: use "mean of stress module - module of the mean stress" data 
    261261      rn_pfac     = 1.        !  multiplicative factor for precipitation (total & snow) 
    262262      rn_efac     = 1.        !  multiplicative factor for evaporation (0. or 1.) 
     
    277277   sn_snow     = 'ncar_precip.15JUNE2009_fill',   -1.        , 'SNOW'    ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , '' 
    278278   sn_slp      = 'slp.15JUNE2009_fill'        ,    6.        , 'SLP'     ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , '' 
    279    sn_tdif     = 'taudif_core'                ,   24         , 'taudif'  ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , '' 
     279   sn_hpgi     = 'NOT USED'                   ,   24.        , 'xxx'  ,   .false.   , .true. , 'yearly'  , '' , ''       , '' 
     280   sn_hpgj     = 'NOT USED'                   ,   24.        , 'xxx'  ,   .false.   , .true. , 'yearly'  , '' , ''       , '' 
    280281/ 
    281282!----------------------------------------------------------------------- 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/SPITZ12/EXPREF/namelist_cfg

    r11268 r11275  
    107107   sn_snow     = 'MARv3.6-9km-Svalbard-2hourly_spitz' ,  2. ,  'snow'    ,   .true.    , .false. , 'yearly'  , 'weights_bilin', '' , '' 
    108108   sn_slp      = 'MARv3.6-9km-Svalbard-2hourly_spitz' ,  2. ,  'slp'     ,   .true.    , .false. , 'yearly'  , 'weights_bilin', '' , '' 
    109    sn_tdif     = 'MARv3.6-9km-Svalbard-2hourly_spitz' ,  2. ,  'tdif'    ,   .true.    , .false. , 'yearly'  , 'weights_bilin', '' , '' 
    110109/ 
    111110!----------------------------------------------------------------------- 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/IOM/iom.F90

    r11223 r11275  
    2929   USE lib_mpp           ! MPP library 
    3030#if defined key_iomput 
    31    USE sbc_oce  , ONLY :   nn_fsbc         ! ocean space and time domain 
    32    USE trc_oce  , ONLY :   nn_dttrc        !  !: frequency of step on passive tracers 
    33    USE icb_oce  , ONLY :   nclasses, class_num       !  !: iceberg classes 
     31   USE sbc_oce  , ONLY :   nn_fsbc, ght_abl, ghw_abl, e3w_abl, jpka, jpkam1 
     32   USE trc_oce  , ONLY :   nn_dttrc                  ! frequency of step on passive tracers 
     33   USE icb_oce  , ONLY :   nclasses, class_num       ! iceberg classes 
    3434#if defined key_si3 
    3535   USE ice      , ONLY :   jpl 
     
    112112      ! 
    113113      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zt_bnds, zw_bnds 
     114      REAL(wp), DIMENSION(2   ,jpkam1)      ::   za_bnds     ! ABL vertical boundaries 
    114115      LOGICAL ::   ll_tmppatch = .TRUE.    !: seb: patch before we remove periodicity 
    115116      INTEGER ::   nldi_save, nlei_save    !:      and close boundaries in output files 
     
    195196      ! vertical grid definition 
    196197      IF(.NOT.llrst_context) THEN 
    197           CALL iom_set_axis_attr( "deptht",  paxis = gdept_1d ) 
    198           CALL iom_set_axis_attr( "depthu",  paxis = gdept_1d ) 
    199           CALL iom_set_axis_attr( "depthv",  paxis = gdept_1d ) 
    200           CALL iom_set_axis_attr( "depthw",  paxis = gdepw_1d ) 
    201  
     198          CALL iom_set_axis_attr(  "deptht", paxis = gdept_1d ) 
     199          CALL iom_set_axis_attr(  "depthu", paxis = gdept_1d ) 
     200          CALL iom_set_axis_attr(  "depthv", paxis = gdept_1d ) 
     201          CALL iom_set_axis_attr(  "depthw", paxis = gdepw_1d ) 
     202 
     203          ! ABL 
     204          IF( .NOT. ALLOCATED(ght_abl) ) THEN   ! force definition for xml files (xios)  
     205             ALLOCATE( ght_abl(jpka), ghw_abl(jpka), e3t_abl(jpka), e3w_abl(jpka) )   ! default allocation needed by iom 
     206             ght_abl(:) = -1._wp   ;   ghw_abl(:) = -1._wp 
     207             e3t_abl(:) = -1._wp   ;   e3w_abl(:) = -1._wp 
     208          ENDIF 
     209          CALL iom_set_axis_attr( "ght_abl", ght_abl(2:jpka) ) 
     210          CALL iom_set_axis_attr( "ghw_abl", ghw_abl(2:jpka) ) 
     211           
    202212          ! Add vertical grid bounds 
    203213          jkmin = MIN(2,jpk)  ! in case jpk=1 (i.e. sas2D) 
     
    208218          zw_bnds(2,1:jpkm1  ) = gdepw_1d(jkmin:jpk) 
    209219          zw_bnds(2,jpk:     ) = gdepw_1d(jpk) + e3t_1d(jpk) 
    210           CALL iom_set_axis_attr( "deptht", bounds=zw_bnds ) 
    211           CALL iom_set_axis_attr( "depthu", bounds=zw_bnds ) 
    212           CALL iom_set_axis_attr( "depthv", bounds=zw_bnds ) 
    213           CALL iom_set_axis_attr( "depthw", bounds=zt_bnds ) 
     220          CALL iom_set_axis_attr(  "deptht", bounds=zw_bnds ) 
     221          CALL iom_set_axis_attr(  "depthu", bounds=zw_bnds ) 
     222          CALL iom_set_axis_attr(  "depthv", bounds=zw_bnds ) 
     223          CALL iom_set_axis_attr(  "depthw", bounds=zt_bnds ) 
     224 
     225          ! ABL 
     226          za_bnds(1,:) = ghw_abl(1:jpkam1) 
     227          za_bnds(2,:) = ghw_abl(2:jpka  ) 
     228          CALL iom_set_axis_attr( "ght_abl", bounds=za_bnds ) 
     229          za_bnds(1,:) = ght_abl(2:jpka  ) 
     230          za_bnds(2,:) = ght_abl(2:jpka  ) + e3w_abl(2:jpka) 
     231          CALL iom_set_axis_attr( "ghw_abl", bounds=za_bnds ) 
    214232          ! 
    215233# if defined key_floats 
     
    19431961      ! 
    19441962      INTEGER :: ji, jj, jn, ni, nj 
    1945       INTEGER :: icnr, jcnr                                    ! Offset such that the vertex coordinate (i+icnr,j+jcnr) 
    1946       !                                                        ! represents the bottom-left corner of cell (i,j) 
    1947       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: z_bnds      ! Lat/lon coordinates of the vertices of cell (i,j) 
    1948       REAL(wp), ALLOCATABLE, DIMENSION(:,:)     :: z_fld       ! Working array to determine where to rotate cells 
    1949       REAL(wp), ALLOCATABLE, DIMENSION(:,:)     :: z_rot       ! Lat/lon working array for rotation of cells 
    1950       !!---------------------------------------------------------------------- 
    1951       ! 
    1952       ALLOCATE( z_bnds(4,jpi,jpj,2), z_fld(jpi,jpj), z_rot(4,2)  ) 
     1963      INTEGER :: icnr, jcnr                             ! Offset such that the vertex coordinate (i+icnr,j+jcnr) 
     1964      !                                                 ! represents the bottom-left corner of cell (i,j) 
     1965      REAL(wp), DIMENSION(4,jpi,jpj,2) ::   z_bnds      ! Lat/lon coordinates of the vertices of cell (i,j) 
     1966      REAL(wp), DIMENSION(  jpi,jpj  ) ::   z_fld       ! Working array to determine where to rotate cells 
     1967      REAL(wp), DIMENSION(4        ,2) ::   z_rot       ! Lat/lon working array for rotation of cells 
     1968      !!---------------------------------------------------------------------- 
    19531969      ! 
    19541970      ! Offset of coordinate representing bottom-left corner 
     
    20312047          &                                    bounds_lon = RESHAPE(z_bnds(:,nldi:nlei,nldj:nlej,2),(/ 4,ni*nj /)), nvertex=4 ) 
    20322048      ! 
    2033       DEALLOCATE( z_bnds, z_fld, z_rot )  
    2034       ! 
    20352049   END SUBROUTINE set_grid_bounds 
    20362050 
     
    21152129      f_op%timestep = nn_fsbc  ;  f_of%timestep =  0  ; CALL iom_set_field_attr('SBC'             , freq_op=f_op, freq_offset=f_of) 
    21162130      f_op%timestep = nn_fsbc  ;  f_of%timestep =  0  ; CALL iom_set_field_attr('SBC_scalar'      , freq_op=f_op, freq_offset=f_of) 
     2131      f_op%timestep = nn_fsbc  ;  f_of%timestep =  0  ; CALL iom_set_field_attr('ABL'             , freq_op=f_op, freq_offset=f_of) 
    21172132      f_op%timestep = nn_dttrc ;  f_of%timestep =  0  ; CALL iom_set_field_attr('ptrc_T'          , freq_op=f_op, freq_offset=f_of) 
    21182133      f_op%timestep = nn_dttrc ;  f_of%timestep =  0  ; CALL iom_set_field_attr('diad_T'          , freq_op=f_op, freq_offset=f_of) 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbc_oce.F90

    r10882 r11275  
    1111   !!            4.0  ! 2012-05  (C. Rousset) add attenuation coef for use in ice model  
    1212   !!            4.0  ! 2016-06  (L. Brodeau) new unified bulk routine (based on AeroBulk) 
     13   !!            4.0  ! 2019-03  (F. Lemarié, G. Samson) add compatibility with ABL mode     
    1314   !!---------------------------------------------------------------------- 
    1415 
     
    3435   LOGICAL , PUBLIC ::   ln_flx         !: flux      formulation 
    3536   LOGICAL , PUBLIC ::   ln_blk         !: bulk formulation 
     37   LOGICAL , PUBLIC ::   ln_abl         !: Atmospheric boundary layer model 
    3638#if defined key_oasis3 
    3739   LOGICAL , PUBLIC ::   lk_oasis = .TRUE.  !: OASIS used 
     
    7779   INTEGER , PUBLIC, PARAMETER ::   jp_flx     = 2        !: flux                          formulation 
    7880   INTEGER , PUBLIC, PARAMETER ::   jp_blk     = 3        !: bulk                          formulation 
    79    INTEGER , PUBLIC, PARAMETER ::   jp_purecpl = 4        !: Pure ocean-atmosphere Coupled formulation 
    80    INTEGER , PUBLIC, PARAMETER ::   jp_none    = 5        !: for OPA when doing coupling via SAS module 
     81   INTEGER , PUBLIC, PARAMETER ::   jp_abl     = 4        !: Atmospheric boundary layer    formulation 
     82   INTEGER , PUBLIC, PARAMETER ::   jp_purecpl = 5        !: Pure ocean-atmosphere Coupled formulation 
     83   INTEGER , PUBLIC, PARAMETER ::   jp_none    = 6        !: for OPA when doing coupling via SAS module 
    8184    
    8285   !!---------------------------------------------------------------------- 
     
    107110   INTEGER , PUBLIC ::  ncpl_qsr_freq            !: qsr coupling frequency per days from atmosphere 
    108111   ! 
    109    LOGICAL , PUBLIC ::   lhftau = .FALSE.        !: HF tau used in TKE: mean(stress module) - module(mean stress) 
    110112   !!                                   !!   now    ! before   !! 
    111113   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   utau   , utau_b   !: sea surface i-stress (ocean referential)     [N/m2] 
    112114   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   vtau   , vtau_b   !: sea surface j-stress (ocean referential)     [N/m2] 
    113115   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   taum              !: module of sea surface stress (at T-point)    [N/m2]  
    114    !! wndm is used onmpute surface gases exchanges in ice-free ocean or leads 
     116   !! wndm is used compute surface gases exchanges in ice-free ocean or leads 
    115117   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wndm              !: wind speed module at T-point (=|U10m-Uoce|)  [m/s] 
    116118   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qsr               !: sea heat flux:     solar                     [W/m2] 
     
    136138   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   atm_co2           !: atmospheric pCO2                             [ppm] 
    137139   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xcplmask          !: coupling mask for ln_mixcpl (warning: allocated in sbccpl) 
     140 
     141   !!--------------------------------------------------------------------- 
     142   !! ABL Vertical Domain size   
     143   !!--------------------------------------------------------------------- 
     144   INTEGER , PUBLIC            ::   jpka   = 2     !: ABL number of vertical levels (default definition) 
     145   INTEGER , PUBLIC            ::   jpkam1 = 1     !: jpka-1 
     146   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)   ::   ght_abl, ghw_abl          !: ABL geopotential height (needed for iom) 
     147   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)   ::   e3t_abl, e3w_abl          !: ABL vertical scale factors (needed for iom) 
    138148 
    139149   !!---------------------------------------------------------------------- 
     
    182192         &      ssu_m  (jpi,jpj) , sst_m(jpi,jpj) , frq_m(jpi,jpj) ,      & 
    183193         &      ssv_m  (jpi,jpj) , sss_m(jpi,jpj) , ssh_m(jpi,jpj) , STAT=ierr(4) ) 
    184          ! 
     194      ! 
    185195      ALLOCATE( e3t_m(jpi,jpj) , STAT=ierr(5) ) 
    186          ! 
     196      ! 
    187197      sbc_oce_alloc = MAXVAL( ierr ) 
    188198      CALL mpp_sum ( 'sbc_oce', sbc_oce_alloc ) 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbcblk.F90

    r11267 r11275  
    1818   !!            4.0  !  2016-10  (G. Madec)  introduce a sbc_blk_init routine 
    1919   !!            4.0  !  2016-10  (M. Vancoppenolle)  Introduce conduction flux emulator (M. Vancoppenolle)  
     20   !!            4.0  !  2019-03  (F. Lemarié & G. Samson)  add ABL compatibility (ln_abl=TRUE) 
    2021   !!---------------------------------------------------------------------- 
    2122 
     
    2324   !!   sbc_blk_init  : initialisation of the chosen bulk formulation as ocean surface boundary condition 
    2425   !!   sbc_blk       : bulk formulation as ocean surface boundary condition 
    25    !!   blk_oce       : computes momentum, heat and freshwater fluxes over ocean 
     26   !!       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)  
    2628   !!   rho_air       : density of (moist) air (depends on T_air, q_air and SLP 
    2729   !!   cp_air        : specific heat of (moist) air (depends spec. hum. q_air) 
     
    6567   PUBLIC   sbc_blk_init  ! called in sbcmod 
    6668   PUBLIC   sbc_blk       ! called in sbcmod 
     69   PUBLIC   blk_oce_1     ! called in sbcabl 
     70   PUBLIC   blk_oce_2     ! called in sbcabl 
     71   PUBLIC   rho_air       ! called in ablmod 
     72   PUBLIC   cp_air        ! called in ablmod 
    6773#if defined key_si3 
    6874   PUBLIC   blk_ice_tau   ! routine called in icesbc 
     
    7177#endif  
    7278 
    73 !!Lolo: should ultimately be moved in the module with all physical constants ? 
     79  INTERFACE cp_air 
     80    MODULE PROCEDURE cp_air_0d, cp_air_2d 
     81  END INTERFACE 
     82 
     83  !!Lolo: should ultimately be moved in the module with all physical constants ? 
    7484!!gm  : In principle, yes. 
    7585   REAL(wp), PARAMETER ::   Cp_dry = 1005.0       !: Specic heat of dry air, constant pressure      [J/K/kg] 
     
    8090   REAL(wp), PARAMETER ::   rctv0 = R_vap/R_dry   !: for virtual temperature (== (1-eps)/eps) => ~ 0.608 
    8191 
    82    INTEGER , PARAMETER ::   jpfld   =10           ! maximum number of files to read 
    83    INTEGER , PARAMETER ::   jp_wndi = 1           ! index of 10m wind velocity (i-component) (m/s)    at T-point 
    84    INTEGER , PARAMETER ::   jp_wndj = 2           ! index of 10m wind velocity (j-component) (m/s)    at T-point 
    85    INTEGER , PARAMETER ::   jp_tair = 3           ! index of 10m air temperature             (Kelvin) 
    86    INTEGER , PARAMETER ::   jp_humi = 4           ! index of specific humidity               ( % ) 
    87    INTEGER , PARAMETER ::   jp_qsr  = 5           ! index of solar heat                      (W/m2) 
    88    INTEGER , PARAMETER ::   jp_qlw  = 6           ! index of Long wave                       (W/m2) 
    89    INTEGER , PARAMETER ::   jp_prec = 7           ! index of total precipitation (rain+snow) (Kg/m2/s) 
    90    INTEGER , PARAMETER ::   jp_snow = 8           ! index of snow (solid prcipitation)       (kg/m2/s) 
    91    INTEGER , PARAMETER ::   jp_slp  = 9           ! index of sea level pressure              (Pa) 
    92    INTEGER , PARAMETER ::   jp_tdif =10           ! index of tau diff associated to HF tau   (N/m2)   at T-point 
    93  
    94    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf   ! structure of input fields (file informations, fields read) 
     92   INTEGER , PUBLIC, PARAMETER ::   jpfld   =11   !: maximum number of files to read 
     93   INTEGER , PUBLIC, PARAMETER ::   jp_wndi = 1   !: index of 10m wind velocity (i-component) (m/s)    at T-point 
     94   INTEGER , PUBLIC, PARAMETER ::   jp_wndj = 2   !: index of 10m wind velocity (j-component) (m/s)    at T-point 
     95   INTEGER , PUBLIC, PARAMETER ::   jp_tair = 3   !: index of 10m air temperature             (Kelvin) 
     96   INTEGER , PUBLIC, PARAMETER ::   jp_humi = 4   !: index of specific humidity               ( % ) 
     97   INTEGER , PUBLIC, PARAMETER ::   jp_qsr  = 5   !: index of solar heat                      (W/m2) 
     98   INTEGER , PUBLIC, PARAMETER ::   jp_qlw  = 6   !: index of Long wave                       (W/m2) 
     99   INTEGER , PUBLIC, PARAMETER ::   jp_prec = 7   !: index of total precipitation (rain+snow) (Kg/m2/s) 
     100   INTEGER , PUBLIC, PARAMETER ::   jp_snow = 8   !: index of snow (solid prcipitation)       (kg/m2/s) 
     101   INTEGER , PUBLIC, PARAMETER ::   jp_slp  = 9   !: index of sea level pressure              (Pa) 
     102   INTEGER , PUBLIC, PARAMETER ::   jp_hpgi =10   !: index of ABL geostrophic wind or hpg (i-component) (m/s) at T-point 
     103   INTEGER , PUBLIC, PARAMETER ::   jp_hpgj =11   !: index of ABL geostrophic wind or hpg (j-component) (m/s) at T-point 
     104 
     105   TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) ::   sf   !  structure of input fields (file informations, fields read) 
    95106 
    96107   !                                             !!! Bulk parameters 
     
    107118   LOGICAL  ::   ln_ECMWF       ! "ECMWF"     algorithm   (IFS cycle 31) 
    108119   ! 
    109    LOGICAL  ::   ln_taudif      ! logical flag to use the "mean of stress module - module of mean stress" data 
    110120   REAL(wp) ::   rn_pfac        ! multiplication factor for precipitation 
    111121   REAL(wp) ::   rn_efac        ! multiplication factor for evaporation 
     
    161171      !! 
    162172      !!---------------------------------------------------------------------- 
    163       INTEGER  ::   ifpr, jfld            ! dummy loop indice and argument 
     173      INTEGER  ::   ifpr                  ! dummy loop indice and argument 
    164174      INTEGER  ::   ios, ierror, ioptio   ! Local integer 
    165175      !! 
     
    168178      TYPE(FLD_N) ::   sn_wndi, sn_wndj, sn_humi, sn_qsr       ! informations about the fields to be read 
    169179      TYPE(FLD_N) ::   sn_qlw , sn_tair, sn_prec, sn_snow      !       "                        " 
    170       TYPE(FLD_N) ::   sn_slp , sn_tdif                        !       "                        " 
     180      TYPE(FLD_N) ::   sn_slp , sn_hpgi, sn_hpgj               !       "                        " 
    171181      NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw ,                &   ! input fields 
    172          &                 sn_tair, sn_prec, sn_snow, sn_slp, sn_tdif,                & 
     182         &                 sn_tair, sn_prec, sn_snow, sn_slp, sn_hpgi, sn_hpgj,       & 
    173183         &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF,             &   ! bulk algorithm 
    174          &                 cn_dir , ln_taudif, rn_zqt, rn_zu,                         &  
     184         &                 cn_dir , rn_zqt, rn_zu,                                    &  
    175185         &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15 
    176186      !!--------------------------------------------------------------------- 
     
    214224      slf_i(jp_tair) = sn_tair   ;   slf_i(jp_humi) = sn_humi 
    215225      slf_i(jp_prec) = sn_prec   ;   slf_i(jp_snow) = sn_snow 
    216       slf_i(jp_slp)  = sn_slp    ;   slf_i(jp_tdif) = sn_tdif 
    217       ! 
    218       lhftau = ln_taudif                     !- add an extra field if HF stress is used 
    219       jfld = jpfld - COUNT( (/.NOT.lhftau/) ) 
     226      slf_i(jp_slp)  = sn_slp 
     227      slf_i(jp_hpgi) = sn_hpgi   ;   slf_i(jp_hpgj) = sn_hpgj   
    220228      ! 
    221229      !                                      !- allocate the bulk structure 
    222       ALLOCATE( sf(jfld), STAT=ierror ) 
     230      ALLOCATE( sf(jpfld), STAT=ierror ) 
    223231      IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_init: unable to allocate sf structure' ) 
    224       DO ifpr= 1, jfld 
    225          ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) ) 
    226          IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) ) 
    227          IF( slf_i(ifpr)%freqh > 0. .AND. MOD( NINT(3600. * slf_i(ifpr)%freqh), nn_fsbc * NINT(rdt) ) /= 0 )   & 
    228             &  CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rdt*nn_fsbc is NOT a submultiple of atmospheric forcing frequency.',   & 
    229             &                 '               This is not ideal. You should consider changing either rdt or nn_fsbc value...' ) 
    230  
    231       END DO 
    232232      !                                      !- fill the bulk structure with namelist informations 
    233233      CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_init', 'surface boundary condition -- bulk formulae', 'namsbc_blk' ) 
     234      ! 
     235      DO ifpr= 1, jpfld 
     236         ! 
     237         IF( TRIM(sf(ifpr)%clrootname) /= 'NOT USED' ) THEN 
     238            IF( ln_abl .AND. (     ifpr == jp_wndi .OR. ifpr == jp_wndj .OR. ifpr == jp_humi   & 
     239               &              .OR. ifpr == jp_hpgi .OR. ifpr == jp_hpgj .OR. ifpr == jp_tair ) ) THEN 
     240               ALLOCATE( sf(ifpr)%fnow(jpi,jpj,jpka) ) 
     241               IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf(ifpr)%fdta(jpi,jpj,jpka,2) ) 
     242            ELSE 
     243               ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) ) 
     244               IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) ) 
     245            ENDIF 
     246            ! 
     247            IF( slf_i(ifpr)%freqh > 0. .AND. MOD( NINT(3600. * slf_i(ifpr)%freqh), nn_fsbc * NINT(rdt) ) /= 0 )   & 
     248               &  CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rdt*nn_fsbc is NOT a submultiple of atmospheric forcing frequency.',   & 
     249               &                 '               This is not ideal. You should consider changing either rdt or nn_fsbc value...' ) 
     250         ENDIF 
     251      END DO 
    234252      ! 
    235253      IF ( ln_wave ) THEN 
     
    252270      ENDIF  
    253271      ! 
     272      IF( ln_abl ) THEN      ! ABL: read 3D fields for wind, temperature and humidity      
     273         rn_zqt = ght_abl(2)          ! set the bulk altitude to ABL first level 
     274         rn_zu  = ght_abl(2) 
     275         IF(lwp) WRITE(numout,*) 
     276         IF(lwp) WRITE(numout,*) '   ABL formulation: overwrite rn_zqt & rn_zu with ABL first level altitude' 
     277      ENDIF 
    254278      !            
    255279      IF(lwp) THEN                     !** Control print 
     
    261285         WRITE(numout,*) '      "COARE 3.5" algorithm   (Edson et al. 2013)         ln_COARE_3p5 = ', ln_COARE_3p0 
    262286         WRITE(numout,*) '      "ECMWF"     algorithm   (IFS cycle 31)              ln_ECMWF     = ', ln_ECMWF 
    263          WRITE(numout,*) '      add High freq.contribution to the stress module     ln_taudif    = ', ln_taudif 
    264287         WRITE(numout,*) '      Air temperature and humidity reference height (m)   rn_zqt       = ', rn_zqt 
    265288         WRITE(numout,*) '      Wind vector reference height (m)                    rn_zu        = ', rn_zu 
     
    294317      !!      the 10m wind velocity (i-component) (m/s)    at T-point 
    295318      !!      the 10m wind velocity (j-component) (m/s)    at T-point 
    296       !!      the 10m or 2m specific humidity     ( % ) 
     319      !!      the 10m or 2m specific humidity     (kg/kg) 
    297320      !!      the solar heat                      (W/m2) 
    298321      !!      the Long wave                       (W/m2) 
    299322      !!      the 10m or 2m air temperature       (Kelvin) 
    300323      !!      the total precipitation (rain+snow) (Kg/m2/s) 
    301       !!      the snow (solid prcipitation)       (kg/m2/s) 
    302       !!      the tau diff associated to HF tau   (N/m2)   at T-point   (ln_taudif=T) 
    303       !!              (2) CALL blk_oce 
     324      !!      the snow (solid precipitation)      (kg/m2/s) 
     325      !!      ABL dynamical forcing (i/j-components of either hpg or geostrophic winds)  
     326      !!              (2) CALL blk_oce_1 and blk_oce_2 
    304327      !! 
    305328      !!      C A U T I O N : never mask the surface stress fields 
     
    318341      !!---------------------------------------------------------------------- 
    319342      INTEGER, INTENT(in) ::   kt   ! ocean time step 
     343      !! 
     344      REAL(wp), DIMENSION(jpi,jpj) ::   zssq, zcd_du, zsen, zevp 
    320345      !!--------------------------------------------------------------------- 
    321346      ! 
     
    323348      ! 
    324349      !                                            ! compute the surface ocean fluxes using bulk formulea 
    325       IF( MOD( kt - 1, nn_fsbc ) == 0 )   CALL blk_oce( kt, sf, sst_m, ssu_m, ssv_m ) 
    326  
     350      IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 
     351         ! 
     352         CALL blk_oce_1( kt, sf(jp_wndi)%fnow(:,:,1), sf(jp_wndj)%fnow(:,:,1),   &   !   <<= in 
     353            &                sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1),   &   !   <<= in 
     354            &                sf(jp_slp )%fnow(:,:,1), sst_m, ssu_m, ssv_m,       &   !   <<= in 
     355            &                zssq, zcd_du, zsen, zevp )                              !   =>> out 
     356          
     357         CALL blk_oce_2( kt, sf(jp_tair)%fnow(:,:,1), sf(jp_qsr )%fnow(:,:,1),   &   !   <<= in 
     358            &                sf(jp_qlw )%fnow(:,:,1), sf(jp_prec)%fnow(:,:,1),   &   !   <<= in 
     359            &                sf(jp_snow)%fnow(:,:,1), sst_m,                     &   !   <<= in 
     360            &                zsen, zevp )                                            !   <=> in out 
     361          
     362      ENDIF 
     363          
    327364#if defined key_cice 
    328365      IF( MOD( kt - 1, nn_fsbc ) == 0 )   THEN 
     
    343380 
    344381 
    345    SUBROUTINE blk_oce( kt, sf, pst, pu, pv ) 
    346       !!--------------------------------------------------------------------- 
    347       !!                     ***  ROUTINE blk_oce  *** 
    348       !! 
    349       !! ** Purpose :   provide the momentum, heat and freshwater fluxes at 
    350       !!      the ocean surface at each time step 
    351       !! 
    352       !! ** Method  :   bulk formulea for the ocean using atmospheric 
    353       !!      fields read in sbc_read 
     382   SUBROUTINE blk_oce_1( kt, pwndi, pwndj, ptair, phumi, &  ! inp 
     383              &              pslp,  pst  , pu   , pv,    &  ! inp 
     384              &              pssq, pcd_du, psen, pevp    )  ! out 
     385      !!--------------------------------------------------------------------- 
     386      !!                     ***  ROUTINE blk_oce_1  *** 
     387      !! 
     388      !! ** Purpose :   Computes Cd x |U|, Ch x |U|, Ce x |U| for ABL integration 
     389      !! 
     390      !! ** Method  :   bulk formulae using atmospheric 
     391      !!                fields from the ABL model at previous time-step 
     392      !! 
     393      !! ** Outputs : - pssq    : surface humidity used to compute latent heat flux (kg/kg) 
     394      !!              - pcd_du  : Cd x |dU| at T-points  (m/s) 
     395      !!              - psen    : Ch x |dU| at T-points  (m/s) 
     396      !!              - pevp    : Ce x |dU| at T-points  (m/s) 
     397      !!--------------------------------------------------------------------- 
     398      INTEGER , INTENT(in   )                 ::   kt     ! time step index 
     399      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pwndi  ! atmospheric wind at U-point              [m/s] 
     400      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pwndj  ! atmospheric wind at V-point              [m/s] 
     401      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   phumi  ! specific humidity at T-points            [kg/kg] 
     402      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   ptair  ! potential temperature at T-points        [Kelvin] 
     403      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pslp   ! sea-level pressure                       [Pa] 
     404      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pst    ! surface temperature                      [Celcius] 
     405      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pu     ! surface current at U-point (i-component) [m/s] 
     406      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pv     ! surface current at V-point (j-component) [m/s] 
     407      REAL(wp), INTENT(  out), DIMENSION(:,:) ::   pssq   ! specific humidity at pst                 [kg/kg] 
     408      REAL(wp), INTENT(  out), DIMENSION(:,:) ::   pcd_du ! Cd x |dU| at T-points                    [m/s] 
     409      REAL(wp), INTENT(  out), DIMENSION(:,:) ::   psen   ! Ch x |dU| at T-points                    [m/s] 
     410      REAL(wp), INTENT(  out), DIMENSION(:,:) ::   pevp   ! Ce x |dU| at T-points                    [m/s] 
     411      ! 
     412      INTEGER  ::   ji, jj               ! dummy loop indices 
     413      REAL(wp) ::   zztmp                ! local variable 
     414      REAL(wp), DIMENSION(jpi,jpj) ::   zwnd_i, zwnd_j    ! wind speed components at T-point 
     415      REAL(wp), DIMENSION(jpi,jpj) ::   zst               ! surface temperature in Kelvin 
     416      REAL(wp), DIMENSION(jpi,jpj) ::   zU_zu             ! bulk wind speed at height zu  [m/s] 
     417      REAL(wp), DIMENSION(jpi,jpj) ::   ztpot             ! potential temperature of air at z=rn_zqt [K] 
     418      REAL(wp), DIMENSION(jpi,jpj) ::   zrhoa             ! density of air   [kg/m^3] 
     419      !!--------------------------------------------------------------------- 
     420      ! 
     421      ! local scalars ( place there for vector optimisation purposes) 
     422      zst(:,:) = pst(:,:) + rt0      ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K) 
     423 
     424      ! ----------------------------------------------------------------------------- ! 
     425      !      0   Wind components and module at T-point relative to the moving ocean   ! 
     426      ! ----------------------------------------------------------------------------- ! 
     427 
     428      ! ... components ( U10m - U_oce ) at T-point (unmasked) 
     429#if defined key_cyclone 
     430      zwnd_i(:,:) = 0._wp 
     431      zwnd_j(:,:) = 0._wp 
     432      CALL wnd_cyc( kt, zwnd_i, zwnd_j )    ! add analytical tropical cyclone (Vincent et al. JGR 2012) 
     433      DO jj = 2, jpjm1 
     434         DO ji = fs_2, fs_jpim1   ! vect. opt. 
     435            pwndi(ji,jj) = pwndi(ji,jj) + zwnd_i(ji,jj) 
     436            pwndj(ji,jj) = pwndj(ji,jj) + zwnd_j(ji,jj) 
     437         END DO 
     438      END DO 
     439#endif 
     440      DO jj = 2, jpjm1 
     441         DO ji = fs_2, fs_jpim1   ! vect. opt. 
     442            zwnd_i(ji,jj) = (  pwndi(ji,jj) - rn_vfac * 0.5 * ( pu(ji-1,jj  ) + pu(ji,jj) )  ) 
     443            zwnd_j(ji,jj) = (  pwndj(ji,jj) - rn_vfac * 0.5 * ( pv(ji  ,jj-1) + pv(ji,jj) )  ) 
     444         END DO 
     445      END DO 
     446      CALL lbc_lnk_multi( 'sbcblk', zwnd_i, 'T', -1., zwnd_j, 'T', -1. ) 
     447      ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked) 
     448      wndm(:,:) = SQRT(  zwnd_i(:,:) * zwnd_i(:,:)   & 
     449         &             + zwnd_j(:,:) * zwnd_j(:,:)  ) * tmask(:,:,1) 
     450 
     451      ! ----------------------------------------------------------------------------- ! 
     452      !     I    Turbulent FLUXES                                                    ! 
     453      ! ----------------------------------------------------------------------------- ! 
     454 
     455      ! ... specific humidity at SST and IST tmask( 
     456      pssq(:,:) = 0.98 * q_sat( zst(:,:), pslp(:,:) ) 
     457      ! 
     458      IF( ln_abl ) THEN 
     459         ztpot = ptair(:,:) 
     460      ELSE 
     461         ! Estimate of potential temperature at z=rn_zqt, based on adiabatic lapse-rate 
     462         !    (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2 
     463         !    (since reanalysis products provide T at z, not theta !) 
     464         ztpot = ptair(:,:) + gamma_moist( ptair(:,:), phumi(:,:) ) * rn_zqt 
     465      ENDIF 
     466 
     467      SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point 
     468      ! 
     469      CASE( np_NCAR      )   ;   CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, pssq, phumi, wndm,   &  ! NCAR-COREv2 
     470         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
     471      CASE( np_COARE_3p0 )   ;   CALL turb_coare   ( rn_zqt, rn_zu, zst, ztpot, pssq, phumi, wndm,   &  ! COARE v3.0 
     472         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
     473      CASE( np_COARE_3p5 )   ;   CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, pssq, phumi, wndm,   &  ! COARE v3.5 
     474         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
     475      CASE( np_ECMWF     )   ;   CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, pssq, phumi, wndm,   &  ! ECMWF 
     476         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
     477      CASE DEFAULT 
     478         CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' ) 
     479      END SELECT 
     480 
     481 
     482      IF ( ln_abl ) THEN         !==  ABL formulation  ==!   multiplication by rho_air and turbulent fluxes computation done in ablstp 
     483!! FL do we need this multiplication by tmask ... ??? 
     484         DO jj = 1, jpj 
     485            DO ji = 1, jpi 
     486               zztmp = zU_zu(ji,jj) * tmask(ji,jj,1) 
     487               wndm(ji,jj)   = zztmp                   ! Store zU_zu in wndm to compute ustar2 in ablmod  
     488               pcd_du(ji,jj) = zztmp * Cd_atm(ji,jj) 
     489               psen(ji,jj)   = zztmp * Ch_atm(ji,jj) 
     490               pevp(ji,jj)   = zztmp * Ce_atm(ji,jj)        
     491            END DO 
     492         END DO 
     493 
     494      ELSE                       !==  BLK formulation  ==!   turbulent fluxes computation 
     495          
     496         !                             ! Compute true air density : 
     497         IF( ABS(rn_zu - rn_zqt) > 0.01 ) THEN     ! At zu: (probably useless to remove zrho*grav*rn_zu from SLP...) 
     498            zrhoa(:,:) = rho_air( t_zu(:,:) , q_zu(:,:) , pslp(:,:) ) 
     499         ELSE                                      ! At zt: 
     500            zrhoa(:,:) = rho_air( ptair(:,:), phumi(:,:), pslp(:,:) ) 
     501         END IF 
     502          
     503         DO jj = 1, jpj 
     504            DO ji = 1, jpi 
     505               zztmp         = zrhoa(ji,jj) * zU_zu(ji,jj) 
     506!!gm to be done by someone: check the efficiency of the call of cp_air within do loops  
     507               psen  (ji,jj) = cp_air( q_zu(ji,jj) ) * zztmp * Ch_atm(ji,jj) * (  zst(ji,jj) - t_zu(ji,jj) ) 
     508               pevp  (ji,jj) = rn_efac*MAX( 0._wp,      zztmp * Ce_atm(ji,jj) * ( pssq(ji,jj) - q_zu(ji,jj) ) ) 
     509               zztmp         = zztmp * Cd_atm(ji,jj) 
     510               pcd_du(ji,jj) = zztmp 
     511               taum  (ji,jj) = zztmp * wndm  (ji,jj)  
     512               zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) 
     513               zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj) 
     514            END DO 
     515         END DO 
     516 
     517         CALL iom_put( "taum_oce", taum )   ! output wind stress module 
     518          
     519         ! ... utau, vtau at U- and V_points, resp. 
     520         !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines 
     521         !     Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves 
     522         DO jj = 1, jpjm1 
     523            DO ji = 1, fs_jpim1 
     524               utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj  ) ) & 
     525                  &          * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1)) 
     526               vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji  ,jj+1) ) & 
     527                  &          * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1)) 
     528            END DO 
     529         END DO 
     530         CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. ) 
     531 
     532         IF(ln_ctl) THEN 
     533            CALL prt_ctl( tab2d_1=wndm  , clinfo1=' blk_oce_1: wndm   : ') 
     534            CALL prt_ctl( tab2d_1=utau  , clinfo1=' blk_oce_1: utau   : ', mask1=umask,   & 
     535               &          tab2d_2=vtau  , clinfo2='            vtau   : ', mask2=vmask ) 
     536         ENDIF 
     537         ! 
     538      ENDIF 
     539      ! 
     540      IF(ln_ctl) THEN 
     541         CALL prt_ctl( tab2d_1=pevp  , clinfo1=' blk_oce_1: pevp   : ' ) 
     542         CALL prt_ctl( tab2d_1=psen  , clinfo1=' blk_oce_1: psen   : ' ) 
     543         CALL prt_ctl( tab2d_1=pssq  , clinfo1=' blk_oce_1: pssq   : ' ) 
     544      ENDIF 
     545      ! 
     546  END SUBROUTINE blk_oce_1 
     547 
     548 
     549   SUBROUTINE blk_oce_2( kt, ptair, pqsr, pqlw, pprec,   &   ! <<= in 
     550              &              psnow, pst , psen, pevp     )   ! <<= in 
     551      !!--------------------------------------------------------------------- 
     552      !!                     ***  ROUTINE blk_oce_2  *** 
     553      !! 
     554      !! ** Purpose :   finalize the momentum, heat and freshwater fluxes computation  
     555      !!                at the ocean surface at each time step knowing Cd, Ch, Ce and  
     556      !!                atmospheric variables (from ABL or external data) 
    354557      !! 
    355558      !! ** Outputs : - utau    : i-component of the stress at U-point  (N/m2) 
     
    360563      !!              - qns     : Non Solar heat flux over the ocean    (W/m2) 
    361564      !!              - emp     : evaporation minus precipitation       (kg/m2/s) 
    362       !! 
    363       !!  ** Nota  :   sf has to be a dummy argument for AGRIF on NEC 
    364       !!--------------------------------------------------------------------- 
    365       INTEGER  , INTENT(in   )                 ::   kt    ! time step index 
    366       TYPE(fld), INTENT(inout), DIMENSION(:)   ::   sf    ! input data 
    367       REAL(wp) , INTENT(in)   , DIMENSION(:,:) ::   pst   ! surface temperature                      [Celcius] 
    368       REAL(wp) , INTENT(in)   , DIMENSION(:,:) ::   pu    ! surface current at U-point (i-component) [m/s] 
    369       REAL(wp) , INTENT(in)   , DIMENSION(:,:) ::   pv    ! surface current at V-point (j-component) [m/s] 
     565      !!--------------------------------------------------------------------- 
     566      INTEGER , INTENT(in)                 ::   kt    ! time step index 
     567      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptair 
     568      REAL(wp), INTENT(in), DIMENSION(:,:) ::   pqsr 
     569      REAL(wp), INTENT(in), DIMENSION(:,:) ::   pqlw 
     570      REAL(wp), INTENT(in), DIMENSION(:,:) ::   pprec 
     571      REAL(wp), INTENT(in), DIMENSION(:,:) ::   psnow 
     572      REAL(wp), INTENT(in), DIMENSION(:,:) ::   pst   ! surface temperature                      [Celcius] 
     573      REAL(wp), INTENT(in), DIMENSION(:,:) ::   psen 
     574      REAL(wp), INTENT(in), DIMENSION(:,:) ::   pevp 
    370575      ! 
    371576      INTEGER  ::   ji, jj               ! dummy loop indices 
    372       REAL(wp) ::   zztmp                ! local variable 
    373       REAL(wp), DIMENSION(jpi,jpj) ::   zwnd_i, zwnd_j    ! wind speed components at T-point 
    374       REAL(wp), DIMENSION(jpi,jpj) ::   zsq               ! specific humidity at pst 
    375       REAL(wp), DIMENSION(jpi,jpj) ::   zqlw, zqsb        ! long wave and sensible heat fluxes 
    376       REAL(wp), DIMENSION(jpi,jpj) ::   zqla, zevap       ! latent heat fluxes and evaporation 
     577      REAL(wp) ::   zztmp,zz1,zz2,zz3    ! local variable 
     578      REAL(wp), DIMENSION(jpi,jpj) ::   zqlw              ! long wave and sensible heat fluxes 
     579      REAL(wp), DIMENSION(jpi,jpj) ::   zqla              ! latent heat fluxes and evaporation 
    377580      REAL(wp), DIMENSION(jpi,jpj) ::   zst               ! surface temperature in Kelvin 
    378       REAL(wp), DIMENSION(jpi,jpj) ::   zU_zu             ! bulk wind speed at height zu  [m/s] 
    379       REAL(wp), DIMENSION(jpi,jpj) ::   ztpot             ! potential temperature of air at z=rn_zqt [K] 
    380       REAL(wp), DIMENSION(jpi,jpj) ::   zrhoa             ! density of air   [kg/m^3] 
    381581      !!--------------------------------------------------------------------- 
    382582      ! 
    383583      ! local scalars ( place there for vector optimisation purposes) 
    384584      zst(:,:) = pst(:,:) + rt0      ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K) 
    385  
    386       ! ----------------------------------------------------------------------------- ! 
    387       !      0   Wind components and module at T-point relative to the moving ocean   ! 
    388       ! ----------------------------------------------------------------------------- ! 
    389  
    390       ! ... components ( U10m - U_oce ) at T-point (unmasked) 
    391 !!gm    move zwnd_i (_j) set to zero  inside the key_cyclone ??? 
    392       zwnd_i(:,:) = 0._wp 
    393       zwnd_j(:,:) = 0._wp 
    394 #if defined key_cyclone 
    395       CALL wnd_cyc( kt, zwnd_i, zwnd_j )    ! add analytical tropical cyclone (Vincent et al. JGR 2012) 
    396       DO jj = 2, jpjm1 
    397          DO ji = fs_2, fs_jpim1   ! vect. opt. 
    398             sf(jp_wndi)%fnow(ji,jj,1) = sf(jp_wndi)%fnow(ji,jj,1) + zwnd_i(ji,jj) 
    399             sf(jp_wndj)%fnow(ji,jj,1) = sf(jp_wndj)%fnow(ji,jj,1) + zwnd_j(ji,jj) 
    400          END DO 
    401       END DO 
    402 #endif 
    403       DO jj = 2, jpjm1 
    404          DO ji = fs_2, fs_jpim1   ! vect. opt. 
    405             zwnd_i(ji,jj) = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pu(ji-1,jj  ) + pu(ji,jj) )  ) 
    406             zwnd_j(ji,jj) = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pv(ji  ,jj-1) + pv(ji,jj) )  ) 
    407          END DO 
    408       END DO 
    409       CALL lbc_lnk_multi( 'sbcblk', zwnd_i, 'T', -1., zwnd_j, 'T', -1. ) 
    410       ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked) 
    411       wndm(:,:) = SQRT(  zwnd_i(:,:) * zwnd_i(:,:)   & 
    412          &             + zwnd_j(:,:) * zwnd_j(:,:)  ) * tmask(:,:,1) 
    413585 
    414586      ! ----------------------------------------------------------------------------- ! 
     
    417589 
    418590      ! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle                    ! Short Wave 
    419       zztmp = 1. - albo 
    420       IF( ln_dm2dc ) THEN   ;   qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1) 
    421       ELSE                  ;   qsr(:,:) = zztmp *          sf(jp_qsr)%fnow(:,:,1)   * tmask(:,:,1) 
    422       ENDIF 
    423  
    424       zqlw(:,:) = (  sf(jp_qlw)%fnow(:,:,1) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:)  ) * tmask(:,:,1)   ! Long  Wave 
    425  
    426       ! ----------------------------------------------------------------------------- ! 
    427       !     II    Turbulent FLUXES                                                    ! 
    428       ! ----------------------------------------------------------------------------- ! 
    429  
    430       ! ... specific humidity at SST and IST tmask( 
    431       zsq(:,:) = 0.98 * q_sat( zst(:,:), sf(jp_slp)%fnow(:,:,1) ) 
    432       !! 
    433       !! Estimate of potential temperature at z=rn_zqt, based on adiabatic lapse-rate 
    434       !!    (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2 
    435       !!    (since reanalysis products provide T at z, not theta !) 
    436       ztpot = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1) ) * rn_zqt 
    437  
    438       SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point 
    439       ! 
    440       CASE( np_NCAR      )   ;   CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! NCAR-COREv2 
    441          &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    442       CASE( np_COARE_3p0 )   ;   CALL turb_coare   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! COARE v3.0 
    443          &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    444       CASE( np_COARE_3p5 )   ;   CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! COARE v3.5 
    445          &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    446       CASE( np_ECMWF     )   ;   CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! ECMWF 
    447          &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    448       CASE DEFAULT 
    449          CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' ) 
    450       END SELECT 
    451  
    452       !                          ! Compute true air density : 
    453       IF( ABS(rn_zu - rn_zqt) > 0.01 ) THEN     ! At zu: (probably useless to remove zrho*grav*rn_zu from SLP...) 
    454          zrhoa(:,:) = rho_air( t_zu(:,:)              , q_zu(:,:)              , sf(jp_slp)%fnow(:,:,1) ) 
    455       ELSE                                      ! At zt: 
    456          zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 
    457       END IF 
    458  
    459 !!      CALL iom_put( "Cd_oce", Cd_atm)  ! output value of pure ocean-atm. transfer coef. 
    460 !!      CALL iom_put( "Ch_oce", Ch_atm)  ! output value of pure ocean-atm. transfer coef. 
    461  
    462       DO jj = 1, jpj             ! tau module, i and j component 
    463          DO ji = 1, jpi 
    464             zztmp = zrhoa(ji,jj)  * zU_zu(ji,jj) * Cd_atm(ji,jj)   ! using bulk wind speed 
    465             taum  (ji,jj) = zztmp * wndm  (ji,jj) 
    466             zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) 
    467             zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj) 
    468          END DO 
    469       END DO 
    470  
    471       !                          ! add the HF tau contribution to the wind stress module 
    472       IF( lhftau )   taum(:,:) = taum(:,:) + sf(jp_tdif)%fnow(:,:,1) 
    473  
    474       CALL iom_put( "taum_oce", taum )   ! output wind stress module 
    475  
    476       ! ... utau, vtau at U- and V_points, resp. 
    477       !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines 
    478       !     Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves 
    479       DO jj = 1, jpjm1 
    480          DO ji = 1, fs_jpim1 
    481             utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj  ) ) & 
    482                &          * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1)) 
    483             vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji  ,jj+1) ) & 
    484                &          * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1)) 
    485          END DO 
    486       END DO 
    487       CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. ) 
     591      zztmp = 1._wp - albo 
     592      IF( ln_dm2dc ) THEN   ;   qsr(:,:) = zztmp * sbc_dcy( pqsr(:,:) ) * tmask(:,:,1) 
     593      ELSE                  ;   qsr(:,:) = zztmp *          pqsr(:,:)   * tmask(:,:,1) 
     594      ENDIF 
     595 
     596      zqlw(:,:) = ( pqlw(:,:) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:)  ) * tmask(:,:,1)   ! Long  Wave 
    488597 
    489598      !  Turbulent fluxes over ocean 
    490599      ! ----------------------------- 
    491600 
    492       ! zqla used as temporary array, for rho*U (common term of bulk formulae): 
    493       zqla(:,:) = zrhoa(:,:) * zU_zu(:,:) * tmask(:,:,1) 
    494  
    495       IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN 
    496          !! q_air and t_air are given at 10m (wind reference height) 
    497          zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - sf(jp_humi)%fnow(:,:,1)) ) ! Evaporation, using bulk wind speed 
    498          zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - ztpot(:,:)             )   ! Sensible Heat, using bulk wind speed 
    499       ELSE 
    500          !! q_air and t_air are not given at 10m (wind reference height) 
    501          ! Values of temp. and hum. adjusted to height of wind during bulk algorithm iteration must be used!!! 
    502          zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - q_zu(:,:) ) ) ! Evaporation, using bulk wind speed 
    503          zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - t_zu(:,:) )   ! Sensible Heat, using bulk wind speed 
    504       ENDIF 
    505  
    506       zqla(:,:) = L_vap(zst(:,:)) * zevap(:,:)     ! Latent Heat flux 
    507  
     601      zqla(:,:) = L_vap( zst(:,:) ) * pevp(:,:)     ! Latent Heat flux 
    508602 
    509603      IF(ln_ctl) THEN 
    510          CALL prt_ctl( tab2d_1=zqla  , clinfo1=' blk_oce: zqla   : ', tab2d_2=Ce_atm , clinfo2=' Ce_oce  : ' ) 
    511          CALL prt_ctl( tab2d_1=zqsb  , clinfo1=' blk_oce: zqsb   : ', tab2d_2=Ch_atm , clinfo2=' Ch_oce  : ' ) 
    512          CALL prt_ctl( tab2d_1=zqlw  , clinfo1=' blk_oce: zqlw   : ', tab2d_2=qsr, clinfo2=' qsr : ' ) 
    513          CALL prt_ctl( tab2d_1=zsq   , clinfo1=' blk_oce: zsq    : ', tab2d_2=zst, clinfo2=' zst : ' ) 
    514          CALL prt_ctl( tab2d_1=utau  , clinfo1=' blk_oce: utau   : ', mask1=umask,   & 
    515             &          tab2d_2=vtau  , clinfo2=           ' vtau : ', mask2=vmask ) 
    516          CALL prt_ctl( tab2d_1=wndm  , clinfo1=' blk_oce: wndm   : ') 
    517          CALL prt_ctl( tab2d_1=zst   , clinfo1=' blk_oce: zst    : ') 
     604         CALL prt_ctl( tab2d_1=zqla  , clinfo1=' blk_oce_2: zqla   : ' ) 
     605         CALL prt_ctl( tab2d_1=zqlw  , clinfo1=' blk_oce_2: zqlw   : ', tab2d_2=qsr, clinfo2=' qsr : ' ) 
    518606      ENDIF 
    519607 
     
    522610      ! ----------------------------------------------------------------------------- ! 
    523611      ! 
    524       emp (:,:) = (  zevap(:,:)                                          &   ! mass flux (evap. - precip.) 
    525          &         - sf(jp_prec)%fnow(:,:,1) * rn_pfac  ) * tmask(:,:,1) 
    526       ! 
    527       qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar 
    528          &     - sf(jp_snow)%fnow(:,:,1) * rn_pfac * rLfus                        &   ! remove latent melting heat for solid precip 
    529          &     - zevap(:,:) * pst(:,:) * rcp                                      &   ! remove evap heat content at SST 
    530          &     + ( sf(jp_prec)%fnow(:,:,1) - sf(jp_snow)%fnow(:,:,1) ) * rn_pfac  &   ! add liquid precip heat content at Tair 
    531          &     * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp                          & 
    532          &     + sf(jp_snow)%fnow(:,:,1) * rn_pfac                                &   ! add solid  precip heat content at min(Tair,Tsnow) 
    533          &     * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi 
    534       qns(:,:) = qns(:,:) * tmask(:,:,1) 
     612      emp (:,:) = ( pevp(:,:) - pprec(:,:) * rn_pfac  ) * tmask(:,:,1)             ! mass flux (evap. - precip.) 
     613      ! 
     614      zz1 = rn_pfac * rLfus 
     615      zz2 = rn_pfac * rcp 
     616      zz3 = rn_pfac * rcpi 
     617      qns(:,:) = (  zqlw(:,:) - psen(:,:) - zqla(:,:)                          &   ! Downward Non Solar 
     618         &        - psnow(:,:) * zztmp                                         &   ! remove latent melting heat for solid precip 
     619         &        - pevp(:,:) * pst(:,:) * rcp                                 &   ! remove evap heat content at SST 
     620         &        + ( pprec(:,:) - psnow(:,:) ) * ( ptair(:,:) - rt0 ) * zz2   &   ! add liquid precip heat content at Tair 
     621         &        + psnow(:,:) * ( MIN( ptair(:,:), rt0 ) - rt0 ) * zz3        &   ! add solid  precip heat content at min(Tair,Tsnow) 
     622         &       ) * tmask(:,:,1) 
    535623      ! 
    536624#if defined key_si3 
    537       qns_oce(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                ! non solar without emp (only needed by SI3) 
     625      qns_oce(:,:) = zqlw(:,:) - psen(:,:) - zqla(:,:)                             ! non solar without emp (only needed by SI3) 
    538626      qsr_oce(:,:) = qsr(:,:) 
    539627#endif 
     
    541629      IF ( nn_ice == 0 ) THEN 
    542630         CALL iom_put( "qlw_oce" ,   zqlw )                 ! output downward longwave heat over the ocean 
    543          CALL iom_put( "qsb_oce" , - zqsb )                 ! output downward sensible heat over the ocean 
     631         CALL iom_put( "qsb_oce" , - psen )                 ! output downward sensible heat over the ocean 
    544632         CALL iom_put( "qla_oce" , - zqla )                 ! output downward latent   heat over the ocean 
    545          CALL iom_put( "qemp_oce",   qns-zqlw+zqsb+zqla )   ! output downward heat content of E-P over the ocean 
     633         CALL iom_put( "qemp_oce",   qns-zqlw+psen+zqla )   ! output downward heat content of E-P over the ocean 
    546634         CALL iom_put( "qns_oce" ,   qns  )                 ! output downward non solar heat over the ocean 
    547635         CALL iom_put( "qsr_oce" ,   qsr  )                 ! output downward solar heat over the ocean 
    548636         CALL iom_put( "qt_oce"  ,   qns+qsr )              ! output total downward heat over the ocean 
    549          tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s] 
    550          sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s] 
     637         tprecip(:,:) = pprec(:,:) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s] 
     638         sprecip(:,:) = psnow(:,:) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s] 
    551639         CALL iom_put( 'snowpre', sprecip )                 ! Snow 
    552640         CALL iom_put( 'precip' , tprecip )                 ! Total precipitation 
     
    554642      ! 
    555643      IF(ln_ctl) THEN 
    556          CALL prt_ctl(tab2d_1=zqsb , clinfo1=' blk_oce: zqsb   : ', tab2d_2=zqlw , clinfo2=' zqlw  : ') 
    557          CALL prt_ctl(tab2d_1=zqla , clinfo1=' blk_oce: zqla   : ', tab2d_2=qsr  , clinfo2=' qsr   : ') 
    558          CALL prt_ctl(tab2d_1=pst  , clinfo1=' blk_oce: pst    : ', tab2d_2=emp  , clinfo2=' emp   : ') 
    559          CALL prt_ctl(tab2d_1=utau , clinfo1=' blk_oce: utau   : ', mask1=umask,   & 
    560             &         tab2d_2=vtau , clinfo2=              ' vtau  : ' , mask2=vmask ) 
    561       ENDIF 
    562       ! 
    563    END SUBROUTINE blk_oce 
    564  
    565  
    566  
     644         CALL prt_ctl(tab2d_1=zqlw , clinfo1=' blk_oce_2: zqlw  : ') 
     645         CALL prt_ctl(tab2d_1=zqla , clinfo1=' blk_oce_2: zqla  : ', tab2d_2=qsr  , clinfo2=' qsr   : ') 
     646         CALL prt_ctl(tab2d_1=emp  , clinfo1=' blk_oce_2: emp   : ') 
     647      ENDIF 
     648      ! 
     649   END SUBROUTINE blk_oce_2 
     650 
     651    
    567652   FUNCTION rho_air( ptak, pqa, pslp ) 
    568653      !!------------------------------------------------------------------------------- 
     
    584669 
    585670 
    586    FUNCTION cp_air( pqa ) 
     671   FUNCTION cp_air_0d( pqa ) 
     672      !!------------------------------------------------------------------------------- 
     673      !!                           ***  FUNCTION cp_air  *** 
     674      !! 
     675      !! ** Purpose : Compute specific heat (Cp) of moist air 
     676      !! 
     677      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
     678      !!------------------------------------------------------------------------------- 
     679      REAL(wp), INTENT(in) ::   pqa      ! air specific humidity         [kg/kg] 
     680      REAL(wp)             ::   cp_air_0d! specific heat of moist air   [J/K/kg] 
     681      !!------------------------------------------------------------------------------- 
     682      ! 
     683      Cp_air_0d = Cp_dry + Cp_vap * pqa 
     684      ! 
     685   END FUNCTION cp_air_0d 
     686 
     687 
     688   FUNCTION cp_air_2d( pqa ) 
    587689      !!------------------------------------------------------------------------------- 
    588690      !!                           ***  FUNCTION cp_air  *** 
     
    593695      !!------------------------------------------------------------------------------- 
    594696      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa      ! air specific humidity         [kg/kg] 
    595       REAL(wp), DIMENSION(jpi,jpj)             ::   cp_air   ! specific heat of moist air   [J/K/kg] 
     697      REAL(wp), DIMENSION(jpi,jpj)             ::   cp_air_2d! specific heat of moist air   [J/K/kg] 
    596698      !!------------------------------------------------------------------------------- 
    597699      ! 
    598       Cp_air = Cp_dry + Cp_vap * pqa 
    599       ! 
    600    END FUNCTION cp_air 
     700      Cp_air_2d = Cp_dry + Cp_vap * pqa 
     701      ! 
     702   END FUNCTION cp_air_2d 
    601703 
    602704 
     
    714816      Ce_atm(:,:) = Cd_ice 
    715817 
    716       wndm_ice(:,:) = 0._wp      !!gm brutal.... 
    717  
    718818      ! ------------------------------------------------------------ ! 
    719819      !    Wind module relative to the moving ice ( U10m - U_ice )   ! 
     
    743843      ! Computing density of air! Way denser that 1.2 over sea-ice !!! 
    744844      zrhoa (:,:) =  rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1)) 
    745  
    746       !!gm brutal.... 
    747       utau_ice  (:,:) = 0._wp 
    748       vtau_ice  (:,:) = 0._wp 
    749       !!gm end 
    750845 
    751846      ! ------------------------------------------------------------ ! 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbccpl.F90

    r10617 r11275  
    533533      !                                                      ! ------------------------- ! 
    534534      srcv(jpr_taum)%clname = 'O_TauMod'   ;   IF( TRIM(sn_rcv_taumod%cldes) == 'coupled' )   srcv(jpr_taum)%laction = .TRUE. 
    535       lhftau = srcv(jpr_taum)%laction 
    536535      ! 
    537536      !                                                      ! ------------------------- ! 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbcmod.F90

    r10499 r11275  
    1515   !!            3.6  ! 2014-11  (P. Mathiot, C. Harris) add ice shelves melting 
    1616   !!            4.0  ! 2016-06  (L. Brodeau) new general bulk formulation 
     17   !!            4.0  ! 2019-03  (F. Lemarié & G. Samson)  add ABL compatibility (ln_abl=TRUE) 
    1718   !!---------------------------------------------------------------------- 
    1819 
     
    3233   USE sbcflx         ! surface boundary condition: flux formulation 
    3334   USE sbcblk         ! surface boundary condition: bulk formulation 
     35   USE sbcabl         ! atmospheric boundary layer 
    3436   USE sbcice_if      ! surface boundary condition: ice-if sea-ice model 
    3537#if defined key_si3 
     
    9294      !! 
    9395      NAMELIST/namsbc/ nn_fsbc  ,                                                    & 
    94          &             ln_usr   , ln_flx   , ln_blk       ,                          & 
     96         &             ln_usr   , ln_flx   , ln_blk   , ln_abl,                      & 
    9597         &             ln_cpl   , ln_mixcpl, nn_components,                          & 
    9698         &             nn_ice   , ln_ice_embd,                                       & 
    9799         &             ln_traqsr, ln_dm2dc ,                                         & 
    98100         &             ln_rnf   , nn_fwb   , ln_ssr   , ln_isf    , ln_apr_dyn ,     & 
    99          &             ln_wave  , ln_cdgw  , ln_sdw   , ln_tauwoc  , ln_stcor   ,     & 
     101         &             ln_wave  , ln_cdgw  , ln_sdw   , ln_tauwoc  , ln_stcor  ,     & 
    100102         &             ln_tauw  , nn_lsm, nn_sdrift 
    101103      !!---------------------------------------------------------------------- 
     
    137139         WRITE(numout,*) '         flux         formulation                   ln_flx        = ', ln_flx 
    138140         WRITE(numout,*) '         bulk         formulation                   ln_blk        = ', ln_blk 
     141         WRITE(numout,*) '         ABL          formulation                   ln_abl        = ', ln_abl 
    139142         WRITE(numout,*) '      Type of coupling (Ocean/Ice/Atmosphere) : ' 
    140143         WRITE(numout,*) '         ocean-atmosphere coupled formulation       ln_cpl        = ', ln_cpl 
     
    225228      CASE( 1 )                        !- Ice-cover climatology ("Ice-if" model)   
    226229      CASE( 2 )                        !- SI3  ice model 
     230         IF( .NOT.( ln_blk .OR. ln_cpl .OR. ln_abl ) )   & 
     231            &                   CALL ctl_stop( 'sbc_init : SI3 sea-ice model requires ln_blk or ln_cpl or ln_abl = T' ) 
    227232      CASE( 3 )                        !- CICE ice model 
    228          IF( .NOT.( ln_blk .OR. ln_cpl ) )   CALL ctl_stop( 'sbc_init : CICE sea-ice model requires ln_blk or ln_cpl = T' ) 
    229          IF( lk_agrif                    )   CALL ctl_stop( 'sbc_init : CICE sea-ice model not currently available with AGRIF' )  
     233         IF( .NOT.( ln_blk .OR. ln_cpl .OR. ln_abl ) )   & 
     234            &                   CALL ctl_stop( 'sbc_init : CICE sea-ice model requires ln_blk or ln_cpl or ln_abl = T' ) 
     235         IF( lk_agrif                                )   & 
     236            &                   CALL ctl_stop( 'sbc_init : CICE sea-ice model not currently available with AGRIF' )  
    230237      CASE DEFAULT                     !- not supported 
    231238      END SELECT 
     
    270277      IF( ln_flx          ) THEN   ;   nsbc = jp_flx     ; icpt = icpt + 1   ;   ENDIF       ! flux                 formulation 
    271278      IF( ln_blk          ) THEN   ;   nsbc = jp_blk     ; icpt = icpt + 1   ;   ENDIF       ! bulk                 formulation 
     279      IF( ln_abl          ) THEN   ;   nsbc = jp_abl     ; icpt = icpt + 1   ;   ENDIF       ! ABL                  formulation 
    272280      IF( ll_purecpl      ) THEN   ;   nsbc = jp_purecpl ; icpt = icpt + 1   ;   ENDIF       ! Pure Coupled         formulation 
    273281      IF( ll_opa          ) THEN   ;   nsbc = jp_none    ; icpt = icpt + 1   ;   ENDIF       ! opa coupling via SAS module 
     
    281289         CASE( jp_flx     )   ;   WRITE(numout,*) '   ==>>>   flux formulation' 
    282290         CASE( jp_blk     )   ;   WRITE(numout,*) '   ==>>>   bulk formulation' 
     291         CASE( jp_abl     )   ;   WRITE(numout,*) '   ==>>>   ABL  formulation' 
    283292         CASE( jp_purecpl )   ;   WRITE(numout,*) '   ==>>>   pure coupled formulation' 
    284293!!gm abusive use of jp_none ??   ===>>> need to be check and changed by adding a jp_sas parameter 
     
    327336      IF( ln_blk      )   CALL sbc_blk_init            ! bulk formulae initialization 
    328337 
     338      IF( ln_abl      )    CALL sbc_abl_init           ! Atmospheric Boundary Layer (ABL) 
     339 
    329340      IF( ln_ssr      )   CALL sbc_ssr_init            ! Sea-Surface Restoring initialization 
    330341      ! 
     
    355366         CALL iom_set_rstw_var_active('sfx_b') 
    356367      ENDIF 
    357  
     368      !  
    358369   END SUBROUTINE sbc_init 
    359370 
     
    424435         IF( ll_sas    )       CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice )   ! OPA-SAS coupling: SAS receiving fields from OPA 
    425436                               CALL sbc_blk       ( kt )                    ! bulk formulation for the ocean 
     437                               ! 
     438      CASE( jp_abl     ) 
     439         IF( ll_sas    )       CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice )   ! OPA-SAS coupling: SAS receiving fields from OPA 
     440                               CALL sbc_abl       ( kt )                    ! ABL  formulation for the ocean 
    426441                               ! 
    427442      CASE( jp_purecpl )   ;   CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice )   ! pure coupled formulation 
     
    498513            CALL iom_get( numror, jpdom_autoglo, 'utau_b', utau_b, ldxios = lrxios )   ! before i-stress  (U-point) 
    499514            CALL iom_get( numror, jpdom_autoglo, 'vtau_b', vtau_b, ldxios = lrxios )   ! before j-stress  (V-point) 
    500             CALL iom_get( numror, jpdom_autoglo, 'qns_b' , qns_b, ldxios = lrxios )   ! before non solar heat flux (T-point) 
     515            CALL iom_get( numror, jpdom_autoglo,  'qns_b',  qns_b, ldxios = lrxios )   ! before non solar heat flux (T-point) 
    501516            ! The 3D heat content due to qsr forcing is treated in traqsr 
    502517            ! CALL iom_get( numror, jpdom_autoglo, 'qsr_b' , qsr_b, ldxios = lrxios  ) ! before     solar heat flux (T-point) 
Note: See TracChangeset for help on using the changeset viewer.