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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 11395 for NEMO/branches – NEMO

Changeset 11395 for NEMO/branches


Ignore:
Timestamp:
2019-08-02T16:19:00+02:00 (5 years ago)
Author:
mathiot
Message:

ENHANCE-02_ISF_nemo : Initial commit isf simplification (add ISF directory, moved isf routine in and split isf cavity and isf parametrisation, ...) (ticket #2142)

Location:
NEMO/branches/2019/ENHANCE-02_ISF_nemo
Files:
14 added
4 deleted
25 edited
5 copied

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2019/ENHANCE-02_ISF_nemo/cfgs/SHARED/field_def_nemo-oce.xml

    r10824 r11395  
    264264 
    265265          <!-- * variable related to ice shelf forcing * --> 
    266           <field id="fwfisf"       long_name="Ice shelf melting"                             unit="kg/m2/s"  /> 
    267           <field id="fwfisf3d"     long_name="Ice shelf melting"                             unit="kg/m2/s"  grid_ref="grid_T_3D" /> 
    268           <field id="qlatisf"      long_name="Ice shelf latent heat flux"                    unit="W/m2"     /> 
    269           <field id="qlatisf3d"    long_name="Ice shelf latent heat flux"                    unit="W/m2"     grid_ref="grid_T_3D" /> 
    270           <field id="qhcisf"       long_name="Ice shelf heat content flux"                   unit="W/m2"     /> 
    271           <field id="qhcisf3d"     long_name="Ice shelf heat content flux"                   unit="W/m2"     grid_ref="grid_T_3D" /> 
    272           <field id="isfgammat"    long_name="transfert coefficient for isf (temperature) "  unit="m/s"      /> 
    273           <field id="isfgammas"    long_name="transfert coefficient for isf (salinity)    "  unit="m/s"      /> 
    274           <field id="stbl"         long_name="salinity in the Losh tbl                    "  unit="PSU"      /> 
    275           <field id="ttbl"         long_name="temperature in the Losh tbl                 "  unit="C"        /> 
    276           <field id="utbl"         long_name="zonal current in the Losh tbl at T point    "  unit="m/s"      /> 
    277           <field id="vtbl"         long_name="merid current in the Losh tbl at T point    "  unit="m/s"      /> 
    278           <field id="thermald"     long_name="thermal driving of ice shelf melting        "  unit="C"        /> 
    279           <field id="tfrz"         long_name="top freezing point (used to compute melt)   "  unit="C"        /> 
    280           <field id="tinsitu"      long_name="top insitu temperature (used to cmpt melt)  "  unit="C"        /> 
    281           <field id="ustar"        long_name="ustar at T point used in ice shelf melting  "  unit="m/s"      /> 
     266          <field id="fwfisf_cav"      long_name="Ice shelf melting"                             unit="kg/m2/s"  /> 
     267          <field id="fwfisf_par"      long_name="Ice shelf melting"                             unit="kg/m2/s"  /> 
     268          <field id="qoceisf_cav"     long_name="Ice shelf ocean  heat flux"                    unit="W/m2"     /> 
     269          <field id="qoceisf_par"     long_name="Ice shelf ocean  heat flux"                    unit="W/m2"     /> 
     270          <field id="qlatisf_cav"     long_name="Ice shelf latent heat flux"                    unit="W/m2"     /> 
     271          <field id="qlatisf_par"     long_name="Ice shelf latent heat flux"                    unit="W/m2"     /> 
     272          <field id="qhcisf_cav"      long_name="Ice shelf heat content flux"                   unit="W/m2"     /> 
     273          <field id="qhcisf_par"      long_name="Ice shelf heat content flux"                   unit="W/m2"     /> 
     274          <field id="fwfisf3d_cav"    long_name="Ice shelf melting"                             unit="kg/m2/s"  grid_ref="grid_T_3D" /> 
     275          <field id="fwfisf3d_par"    long_name="Ice shelf melting"                             unit="kg/m2/s"  grid_ref="grid_T_3D" /> 
     276          <field id="qoceisf3d_cav"   long_name="Ice shelf ocean  heat flux"                    unit="W/m2"     grid_ref="grid_T_3D" /> 
     277          <field id="qoceisf3d_par"   long_name="Ice shelf ocean  heat flux"                    unit="W/m2"     grid_ref="grid_T_3D" /> 
     278          <field id="qlatisf3d_cav"   long_name="Ice shelf latent heat flux"                    unit="W/m2"     grid_ref="grid_T_3D" /> 
     279          <field id="qlatisf3d_par"   long_name="Ice shelf latent heat flux"                    unit="W/m2"     grid_ref="grid_T_3D" /> 
     280          <field id="qhcisf3d_cav"    long_name="Ice shelf heat content flux"                   unit="W/m2"     grid_ref="grid_T_3D" /> 
     281          <field id="qhcisf3d_par"    long_name="Ice shelf heat content flux"                   unit="W/m2"     grid_ref="grid_T_3D" /> 
     282          <field id="ttbl_cav"        long_name="temperature in the tracer sample depth      "  unit="C"        /> 
     283          <field id="ttbl_par"        long_name="temperature in the tracer sample depth      "  unit="C"        /> 
     284          <field id="isfthermald_cav" long_name="thermal driving of ice shelf melting        "  unit="C"        /> 
     285          <field id="isfthermald_par" long_name="thermal driving of ice shelf melting        "  unit="C"        /> 
     286          <field id="isfgammat"       long_name="transfert coefficient for isf (temperature) "  unit="m/s"      /> 
     287          <field id="isfgammas"       long_name="transfert coefficient for isf (salinity)    "  unit="m/s"      /> 
     288          <field id="stbl"            long_name="salinity in the Losh tbl                    "  unit="PSU"      /> 
     289          <field id="utbl"            long_name="zonal current in the Losh tbl at T point    "  unit="m/s"      /> 
     290          <field id="vtbl"            long_name="merid current in the Losh tbl at T point    "  unit="m/s"      /> 
     291          <field id="isfustar"        long_name="ustar at T point used in ice shelf melting  "  unit="m/s"      /> 
     292          <field id="qconisf"         long_name="Conductive heat flux through the ice shelf"    unit="W/m2"     /> 
    282293 
    283294          <!-- *_oce variables available with ln_blk_clio or ln_blk_core --> 

  • NEMO/branches/2019/ENHANCE-02_ISF_nemo/cfgs/SHARED/namelist_ref

    r10808 r11395  
    5252      cn_ocerst_outdir = "."         !  directory in which to write output ocean restarts 
    5353   ln_iscpl    = .false.   !  cavity evolution forcing or coupling to ice sheet model 
    54    nn_istate   =       0   !  output the initial state (1) or not (0) 
     54   nn_istate   =       1   !  output the initial state (1) or not (0) 
    5555   ln_rst_list = .false.   !  output restarts at list of times using nn_stocklist (T) or at set frequency with nn_stock (F) 
    5656   nn_stock    =    5840   !  frequency of creation of a restart file (modulo referenced to 1) 
    … …  
    6868!----------------------------------------------------------------------- 
    6969   ln_linssh   = .false.   !  =T  linear free surface  ==>>  model level are fixed in time 
    70    rn_isfhmin  =    1.00   !  treshold [m] to discriminate grounding ice from floating ice 
    7170   ! 
    7271   rn_rdt      = 5400.     !  time step for the dynamics and tracer 
    … …  
    7574   ln_crs      = .false.   !  Logical switch for coarsening module      (T => fill namcrs) 
    7675   ! 
    77    ln_meshmask = .false.   !  =T create a mesh file 
     76   ln_meshmask = .true.   !  =T create a mesh file 
    7877/ 
    7978!----------------------------------------------------------------------- 
    … …  
    180179!!   namsbc_rnf      river runoffs                                      (ln_rnf     =T) 
    181180!!   namsbc_apr      Atmospheric Pressure                               (ln_apr_dyn =T) 
    182 !!   namsbc_isf      ice shelf melting/freezing                         (ln_isfcav  =T : read (ln_read_cfg=T) or set or usr_def_zgr ) 
    183181!!   namsbc_iscpl    coupling option between land ice model and ocean   (ln_isfcav  =T) 
    184182!!   namsbc_wave     external fields from wave model                    (ln_wave    =T) 
    … …  
    436434/ 
    437435!----------------------------------------------------------------------- 
    438 &namsbc_isf    !  Top boundary layer (ISF)                              (ln_isfcav =T : read (ln_read_cfg=T)  
    439 !-----------------------------------------------------------------------             or set or usr_def_zgr ) 
     436&namisf       !  Top boundary layer (ISF)                              (ln_isf = T .AND. ln_isfcav: read (ln_read_cfg=T)  
     437!-----------------------------------------------------------------------                         or set or usr_def_zgr ) 
    440438   !                 ! type of top boundary layer  
    441    nn_isf      = 1         !  ice shelf melting/freezing 
    442                            !  1 = presence of ISF   ;  2 = bg03 parametrisation  
    443                            !  3 = rnf file for ISF  ;  4 = ISF specified freshwater flux 
    444                            !  options 1 and 4 need ln_isfcav = .true. (domzgr) 
    445       !              !  nn_isf = 1 or 2 cases: 
     439   ln_isfcav_mlt = .false.    ! ice shelf melting into the cavity 
     440      cn_isfcav_mlt = '3eq'   ! ice shelf melting formulation (spe/2eq/3eq/oasis) 
     441      !                       ! spe = fwfisf is read from a forcing field 
     442      !                       ! 2eq = ISOMIP  like: 2 equations formulation (Hunter et al., 2006) 
     443      !                       ! 3eq = ISOMIP+ like: 3 equations formulation (Asay-Davis et al., 2015) 
     444      !                       ! oasis = fwfisf is given by oasis 
     445      !              !  cn_isfcav_mlt = 2eq or 3eq cases: 
     446      cn_gammablk = 'ad15'    ! scheme to compute gammat/s (spe,ad15,hj99) 
     447      !                       ! ad15 = velocity dependend Gamma (u* * gammat/s)  (Jenkins et al. 2010) 
     448      !                       ! hj99 = velocity and stability dependent Gamma    (Holland et al. 1999) 
    446449      rn_gammat0  = 1.e-4     ! gammat coefficient used in blk formula 
    447450      rn_gammas0  = 1.e-4     ! gammas coefficient used in blk formula 
    448       !              !  nn_isf = 1 or 4 cases: 
    449       rn_hisf_tbl =  30.      ! thickness of the top boundary layer    (Losh et al. 2008) 
     451      ! 
     452      rn_htbl    =  30.      ! thickness of the top boundary layer    (Losh et al. 2008) 
    450453      !                       ! 0 => thickness of the tbl = thickness of the first wet cell 
    451       !              ! nn_isf = 1 case 
    452       nn_isfblk   = 1         ! 1 ISOMIP  like: 2 equations formulation (Hunter et al., 2006) 
    453       !                       ! 2 ISOMIP+ like: 3 equations formulation (Asay-Davis et al., 2015) 
    454       nn_gammablk = 1         ! 0 = cst Gammat (= gammat/s) 
    455       !                       ! 1 = velocity dependend Gamma (u* * gammat/s)  (Jenkins et al. 2010) 
    456       !                       ! 2 = velocity and stability dependent Gamma    (Holland et al. 1999) 
    457  
    458    !___________!_____________!___________________!___________!_____________!_________!___________!__________!__________!_______________! 
    459    !           !  file name  ! frequency (hours) ! variable  ! time interp.!  clim   ! 'yearly'/ ! weights  ! rotation ! land/sea mask ! 
    460    !           !             !  (if <0  months)  !   name    !  (logical)  !  (T/F)  ! 'monthly' ! filename ! pairing  ! filename      ! 
    461 !* nn_isf = 4 case 
    462    sn_fwfisf   = 'rnfisf'    ,         -12       ,'sowflisf' ,  .false.    , .true.  , 'yearly'  ,    ''    ,   ''     ,    '' 
    463 !* nn_isf = 3 case 
    464    sn_rnfisf   = 'rnfisf'    ,         -12       ,'sofwfisf' ,  .false.    , .true.  , 'yearly'  ,    ''    ,   ''     ,    '' 
    465 !* nn_isf = 2 and 3 cases  
    466    sn_depmax_isf ='rnfisf'    ,         -12       ,'sozisfmax',  .false.    , .true.  , 'yearly'  ,    ''    ,   ''     ,    '' 
    467    sn_depmin_isf ='rnfisf'    ,         -12       ,'sozisfmin',  .false.    , .true.  , 'yearly'  ,    ''    ,   ''     ,    '' 
    468 !* nn_isf = 2 case 
    469    sn_Leff_isf = 'rnfisf'    ,         -12       ,'Leff'     ,  .false.    , .true.  , 'yearly'  ,    ''    ,   ''     ,    '' 
     454      ! 
     455      !* 'spe' and 'oasis' case 
     456      !___________!_____________!___________________!___________!_____________!_________!___________!__________!__________!_______________! 
     457      !           !  file name  ! frequency (hours) ! variable  ! time interp.!  clim   ! 'yearly'/ ! weights  ! rotation ! land/sea mask ! 
     458      !           !             !  (if <0  months)  !   name    !  (logical)  !  (T/F)  ! 'monthly' ! filename ! pairing  ! filename      ! 
     459      sn_isfcav_fwf = 'isfmlt_cav',      -12       , 'fwflisf' ,  .false.    , .true.  , 'yearly'  ,    ''    ,   ''     ,    '' 
     460      ! 
     461   ! 
     462   ln_isfpar_mlt = .false.   ! ice shelf melting parametrised 
     463      cn_isfpar_mlt = 'spe'  ! ice shelf melting parametrisation (spe/bg03/oasis) 
     464      ! 
     465      !* all cases 
     466      !___________!_____________!___________________!___________!_____________!_________!___________!__________!__________!_______________! 
     467      !           !  file name  ! frequency (hours) ! variable  ! time interp.!  clim   ! 'yearly'/ ! weights  ! rotation ! land/sea mask ! 
     468      !           !             !  (if <0  months)  !   name    !  (logical)  !  (T/F)  ! 'monthly' ! filename ! pairing  ! filename      ! 
     469      sn_isfpar_zmax = 'isfmlt_par',       0        ,'sozisfmax',  .false.    , .true.  , 'yearly'  ,    ''    ,   ''     ,    '' 
     470      sn_isfpar_zmin = 'isfmlt_par',       0        ,'sozisfmin',  .false.    , .true.  , 'yearly'  ,    ''    ,   ''     ,    '' 
     471      !* 'spe' and 'oasis' case 
     472      sn_isfpar_fwf = 'isfmlt_par' ,      -12       ,'fwfisf'   ,  .false.    , .true.  , 'yearly'  ,    ''    ,   ''     ,    '' 
     473      !* 'bg03' case 
     474      sn_isfpar_Leff = 'isfmlt_par',       0        ,'Leff'     ,  .false.    , .true.  , 'yearly'  ,    ''    ,   ''     ,    '' 
     475      ! 
     476   ! 
     477   !ln_iscpl = .false. 
     478   !   nn_drown    = 10        ! number of iteration of the extrapolation loop (fill the new wet cells) 
     479   !   ln_hsb      = .false.   ! activate conservation module (conservation exact after a time of rn_fiscpl) 
     480   !   nn_fiscpl   = 43800     ! (number of time step) conservation period (maybe should be fix to the coupling frequencey of restart frequency) 
    470481/ 
    471482!----------------------------------------------------------------------- 

  • NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/BDY/bdyvol.F90

    r10481 r11395  
    1616   USE dom_oce        ! ocean space and time domain  
    1717   USE phycst         ! physical constants 
    18    USE sbcisf         ! ice shelf 
     18   USE isf            ! ice shelf 
    1919   ! 
    2020   USE in_out_manager ! I/O manager 
    … …  
    7777      ! Calculate the cumulate surface Flux z_cflxemp (m3/s) over all the domain 
    7878      ! ----------------------------------------------------------------------- 
    79       IF ( kc == 1 ) z_cflxemp = glob_sum( 'bdyvol', ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) * bdytmask(:,:) * e1e2t(:,:)  ) / rau0 
     79      IF ( kc == 1 ) z_cflxemp = glob_sum( 'bdyvol', ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) ) * bdytmask(:,:) * e1e2t(:,:)  ) / rau0 
    8080 
    8181      ! Compute bdy surface each cycle if non linear free surface 

  • NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DIA/diahsb.F90

    r10425 r11395  
    1818   USE sbc_oce        ! surface thermohaline fluxes 
    1919   USE sbcrnf         ! river runoff 
    20    USE sbcisf         ! ice shelves 
     20   USE isf            ! ice shelves 
    2121   USE domvvl         ! vertical scale factors 
    2222   USE traqsr         ! penetrative solar radiation 
    … …  
    9191      ! 1 - Trends due to forcing ! 
    9292      ! ------------------------- ! 
    93       z_frc_trd_v = r1_rau0 * glob_sum( 'diahsb', - ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) * surf(:,:) )   ! volume fluxes 
     93      z_frc_trd_v = r1_rau0 * glob_sum( 'diahsb', - ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) ) * surf(:,:) )   ! volume fluxes 
    9494      z_frc_trd_t =           glob_sum( 'diahsb', sbc_tsc(:,:,jp_tem) * surf(:,:) )                       ! heat fluxes 
    9595      z_frc_trd_s =           glob_sum( 'diahsb', sbc_tsc(:,:,jp_sal) * surf(:,:) )                       ! salt fluxes 
    … …  
    9898      IF( ln_rnf_sal)   z_frc_trd_s = z_frc_trd_s + glob_sum( 'diahsb', rnf_tsc(:,:,jp_sal) * surf(:,:) ) 
    9999      !                    ! Add ice shelf heat & salt input 
    100       IF( ln_isf    )   z_frc_trd_t = z_frc_trd_t + glob_sum( 'diahsb', risf_tsc(:,:,jp_tem) * surf(:,:) ) 
     100      IF( ln_isf    )   z_frc_trd_t = z_frc_trd_t & 
     101         &                          + glob_sum( 'diahsb', ( risf_cav_tsc(:,:,jp_tem) + risf_par_tsc(:,:,jp_tem) ) * surf(:,:) ) 
    101102      !                    ! Add penetrative solar radiation 
    102103      IF( ln_traqsr )   z_frc_trd_t = z_frc_trd_t + r1_rau0_rcp * glob_sum( 'diahsb', qsr     (:,:) * surf(:,:) ) 

  • NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DIA/diawri.F90

    r10425 r11395  
    2626   !!---------------------------------------------------------------------- 
    2727   USE oce            ! ocean dynamics and tracers  
     28   USE isf 
    2829   USE dom_oce        ! ocean space and time domain 
    2930   USE phycst         ! physical constants 
    … …  
    904905      CALL iom_rstput( 0, 0, inum, 'vomecrty', vn                )    ! now j-velocity 
    905906      CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn                )    ! now k-velocity 
     907      CALL iom_rstput( 0, 0, inum, 'risfdep', risfdep          )    ! now k-velocity 
     908      CALL iom_rstput( 0, 0, inum, 'fwfisf_cav', fwfisf_cav          )    ! now k-velocity 
     909      CALL iom_rstput( 0, 0, inum, 'rhisf_cav_tbl', rhisf_tbl_cav    )    ! now k-velocity 
     910      CALL iom_rstput( 0, 0, inum, 'rfrac_cav_tbl', rfrac_tbl_cav    )    ! now k-velocity 
     911      CALL iom_rstput( 0, 0, inum, 'misfkb_cav', REAL(misfkb_cav,8)    )    ! now k-velocity 
     912      CALL iom_rstput( 0, 0, inum, 'misfkt_cav', REAL(misfkt_cav,8)    )    ! now k-velocity 
     913      CALL iom_rstput( 0, 0, inum, 'fwfisf_par', fwfisf_par          )    ! now k-velocity 
     914      CALL iom_rstput( 0, 0, inum, 'rhisf_par_tbl', rhisf_tbl_par    )    ! now k-velocity 
     915      CALL iom_rstput( 0, 0, inum, 'rfrac_par_tbl', rfrac_tbl_par    )    ! now k-velocity 
     916      CALL iom_rstput( 0, 0, inum, 'misfkb_par', REAL(misfkb_par,8)    )    ! now k-velocity 
     917      CALL iom_rstput( 0, 0, inum, 'misfkt_par', REAL(misfkt_par,8)    )    ! now k-velocity 
    906918      IF( ALLOCATED(ahtu) ) THEN 
    907919         CALL iom_rstput( 0, 0, inum,  'ahtu', ahtu              )    ! aht at u-point 

  • NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DOM/dom_oce.F90

    r10068 r11395  
    3232   LOGICAL , PUBLIC ::   ln_linssh      !: =T  linear free surface ==>> model level are fixed in time 
    3333   LOGICAL , PUBLIC ::   ln_meshmask    !: =T  create a mesh-mask file (mesh_mask.nc) 
    34    REAL(wp), PUBLIC ::   rn_isfhmin     !: threshold to discriminate grounded ice to floating ice 
    3534   REAL(wp), PUBLIC ::   rn_rdt         !: time step for the dynamics and tracer 
    3635   REAL(wp), PUBLIC ::   rn_atfp        !: asselin time filter parameter 
    … …  
    170169   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmask_h            !: internal domain T-point mask (Figure 8.5 NEMO book) 
    171170 
    172    INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   misfdep                 !: top first ocean level             (ISF) 
    173    INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mikt, miku, mikv, mikf  !: top first wet T-, U-, V-, F-level (ISF) 
    174    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   risfdep                 !: Iceshelf draft                    (ISF) 
     171   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mikt, miku, mikv, mikf  !: top first wet T-, U-, V-, F-level           (ISF) 
    175172 
    176173   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ssmask, ssumask, ssvmask             !: surface mask at T-,U-, V- and F-pts 
    … …  
    281278         &      mbkt   (jpi,jpj) , mbku   (jpi,jpj) , mbkv   (jpi,jpj) , STAT=ierr(9) ) 
    282279         ! 
    283       ALLOCATE( misfdep(jpi,jpj) , mikt(jpi,jpj) , miku(jpi,jpj) ,     & 
    284          &      risfdep(jpi,jpj) , mikv(jpi,jpj) , mikf(jpi,jpj) , STAT=ierr(10) ) 
     280      ALLOCATE( mikt(jpi,jpj), miku(jpi,jpj), mikv(jpi,jpj), mikf(jpi,jpj), STAT=ierr(10) ) 
    285281         ! 
    286282      ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk) ,     &  

  • NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DOM/domain.F90

    r10425 r11395  
    140140                                       ! Read in masks to define closed seas and lakes  
    141141      ! 
    142       DO jj = 1, jpj                   ! depth of the iceshelves 
    143          DO ji = 1, jpi 
    144             ik = mikt(ji,jj) 
    145             risfdep(ji,jj) = gdepw_0(ji,jj,ik) 
    146          END DO 
    147       END DO 
    148       ! 
    149142      ht_0(:,:) = 0._wp  ! Reference ocean thickness 
    150143      hu_0(:,:) = 0._wp 
    … …  
    293286         &             nn_stock, nn_write , ln_mskland  , ln_clobber   , nn_chunksz, nn_euler  ,     & 
    294287         &             ln_cfmeta, ln_iscpl, ln_xios_read, nn_wxios 
    295       NAMELIST/namdom/ ln_linssh, rn_isfhmin, rn_rdt, rn_atfp, ln_crs, ln_meshmask 
     288      NAMELIST/namdom/ ln_linssh, rn_rdt, rn_atfp, ln_crs, ln_meshmask 
    296289#if defined key_netcdf4 
    297290      NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip 
    … …  
    412405         WRITE(numout,*) '      linear free surface (=T)                ln_linssh   = ', ln_linssh 
    413406         WRITE(numout,*) '      create mesh/mask file                   ln_meshmask = ', ln_meshmask 
    414          WRITE(numout,*) '      treshold to open the isf cavity         rn_isfhmin  = ', rn_isfhmin, ' [m]' 
    415407         WRITE(numout,*) '      ocean time step                         rn_rdt      = ', rn_rdt 
    416408         WRITE(numout,*) '      asselin time filter parameter           rn_atfp     = ', rn_atfp 

  • NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DOM/domwri.F90

    r10425 r11395  
    1616   !!   dom_stiff      : diagnose maximum grid stiffness/hydrostatic consistency (s-coordinate) 
    1717   !!---------------------------------------------------------------------- 
     18   USE isf             ! ice shelf 
    1819   USE dom_oce         ! ocean space and time domain 
    1920   USE phycst ,   ONLY :   rsmall 
    … …  
    155156       
    156157      ! note that mbkt is set to 1 over land ==> use surface tmask 
    157       zprt(:,:) = ssmask(:,:) * REAL( mbkt(:,:) , wp ) 
     158      zprt(:,:) = REAL( mbkt(:,:) , wp ) 
    158159      CALL iom_rstput( 0, 0, inum, 'mbathy', zprt, ktype = jp_i4 )     !    ! nb of ocean T-points 
    159       zprt(:,:) = ssmask(:,:) * REAL( mikt(:,:) , wp ) 
     160      zprt(:,:) = REAL( mikt(:,:) , wp ) 
    160161      CALL iom_rstput( 0, 0, inum, 'misf', zprt, ktype = jp_i4 )       !    ! nb of ocean T-points 
    161       zprt(:,:) = ssmask(:,:) * REAL( risfdep(:,:) , wp ) 
    162       CALL iom_rstput( 0, 0, inum, 'isfdraft', zprt, ktype = jp_r8 )   !    ! nb of ocean T-points 
    163162      !                                                         ! vertical mesh 
    164163      CALL iom_rstput( 0, 0, inum, 'e3t_0', e3t_0, ktype = jp_r8  )    !    ! scale factors 

  • NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN/divhor.F90

    r10425 r11395  
    2222   USE sbc_oce, ONLY : ln_rnf, ln_isf ! surface boundary condition: ocean 
    2323   USE sbcrnf          ! river runoff  
    24    USE sbcisf          ! ice shelf 
     24   USE isfhdiv         ! ice shelf 
    2525   USE iscplhsb        ! ice sheet / ocean coupling 
    2626   USE iscplini        ! ice sheet / ocean coupling 
    … …  
    100100      !  
    101101#endif 
    102       IF( ln_isf )   CALL sbc_isf_div( hdivn )      !==  ice shelf  ==!   (update hdivn field) 
    103       ! 
    104       IF( ln_iscpl .AND. ln_hsb )   CALL iscpl_div( hdivn ) !==  ice sheet  ==!   (update hdivn field) 
     102      IF( ln_isf )   CALL isf_hdiv( hdivn )      !==  ice shelf  ==!   (update hdivn field) 
    105103      ! 
    106104      CALL lbc_lnk( 'divhor', hdivn, 'T', 1. )   !   (no sign change) 

  • NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN/dynhpg.F90

    r10491 r11395  
    3131   !!---------------------------------------------------------------------- 
    3232   USE oce             ! ocean dynamics and tracers 
     33   USE isf             ! ice shelf 
    3334   USE sbc_oce         ! surface variable (only for the flag with ice shelf) 
    3435   USE dom_oce         ! ocean space and time domain 

  • NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN/dynnxt.F90

    r10425 r11395  
    2929   USE sbc_oce        ! Surface boundary condition: ocean fields 
    3030   USE sbcrnf         ! river runoffs 
    31    USE sbcisf         ! ice shelf 
     31   USE isfnxt 
    3232   USE phycst         ! physical constants 
    3333   USE dynadv         ! dynamics: vector invariant versus flux form 
    … …  
    241241               ENDIF 
    242242            END IF 
    243  
    244             IF ( ln_isf ) THEN   ! if ice shelf melting 
    245                DO jk = 1, jpkm1 ! Deal with isf separetely, as can be through depth too 
    246                   DO jj = 1, jpj 
    247                      DO ji = 1, jpi 
    248                         IF( misfkt(ji,jj) <=jk .and. jk < misfkb(ji,jj)  ) THEN 
    249                            e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) & 
    250                                 &          * ( e3t_n(ji,jj,jk) * r1_hisf_tbl(ji,jj) ) * tmask(ji,jj,jk) 
    251                         ELSEIF ( jk==misfkb(ji,jj) ) THEN 
    252                            e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) & 
    253                                 &          * ( e3t_n(ji,jj,jk) * r1_hisf_tbl(ji,jj) ) * ralpha(ji,jj) * tmask(ji,jj,jk) 
    254                         ENDIF 
    255                      END DO 
    256                   END DO 
    257                END DO 
    258             END IF 
     243            ! 
     244            ! ice shelf melting 
     245            IF ( ln_isf ) CALL isf_dynnxt( zcoef ) 
    259246            ! 
    260247            IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity 

  • NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN/dynspg_ts.F90

    r10742 r11395  
    3333   USE zdf_oce         ! vertical physics: variables 
    3434   USE zdfdrg          ! vertical physics: top/bottom drag coef. 
    35    USE sbcisf          ! ice shelf variable (fwfisf) 
     35   USE isf             ! ice shelf variable (fwfisf) 
    3636   USE sbcapr          ! surface boundary condition: atmospheric pressure 
    3737   USE dynadv    , ONLY: ln_dynadv_vec 
    … …  
    632632      !                                         ! Surface net water flux and rivers 
    633633      IF (ln_bt_fw) THEN 
    634          zssh_frc(:,:) = r1_rau0 * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) 
     634         zssh_frc(:,:) = r1_rau0 * ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) ) 
    635635      ELSE 
    636636         zztmp = r1_rau0 * r1_2 
    637637         zssh_frc(:,:) = zztmp * (  emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:)   & 
    638                 &                 + fwfisf(:,:) + fwfisf_b(:,:)                     ) 
     638                &                 + fwfisf_cav(:,:) + fwfisf_cav_b(:,:)             & 
     639                &                 + fwfisf_par(:,:) + fwfisf_par_b(:,:)             ) 
    639640      ENDIF 
    640641      ! 

  • NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DYN/sshwzv.F90

    r11293 r11395  
    1717   !!---------------------------------------------------------------------- 
    1818   USE oce            ! ocean dynamics and tracers variables 
     19   USE isf            ! ice shelf 
    1920   USE dom_oce        ! ocean space and time domain variables  
    2021   USE sbc_oce        ! surface boundary condition: ocean 
    … …  
    256257            sshb(:,:) = sshb(:,:) - zcoef * (     emp_b(:,:) - emp   (:,:)   & 
    257258               &                             -    rnf_b(:,:) + rnf   (:,:)   & 
    258                &                             + fwfisf_b(:,:) - fwfisf(:,:)   ) * ssmask(:,:) 
     259               &                             + fwfisf_cav_b(:,:) - fwfisf_cav(:,:)   & 
     260               &                             + fwfisf_par_b(:,:) - fwfisf_par(:,:)   ) * ssmask(:,:) 
    259261         ENDIF 
    260262         sshn(:,:) = ssha(:,:)                              ! now <-- after 

  • NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/iscplini.F90

    r11310 r11395  
    2323   PUBLIC   iscpl_init       
    2424   PUBLIC   iscpl_alloc  
    25     
    2625   !                                 !!* namsbc_iscpl namelist * 
    27    LOGICAL , PUBLIC ::   ln_hsb      !: 
     26   LOGICAL , PUBLIC ::   ln_iscpl_hsb !: 
    2827   INTEGER , PUBLIC ::   nn_fiscpl    !: 
    2928   INTEGER , PUBLIC ::   nn_drown     !: 
    30     
     29   !  
    3130   INTEGER , PUBLIC ::   nstp_iscpl   !: 
    3231   REAL(wp), PUBLIC ::   rdt_iscpl    !:  
    … …  
    5756      !!---------------------------------------------------------------------- 
    5857      INTEGER ::   ios           ! Local integer output status for namelist read 
    59       NAMELIST/namsbc_iscpl/ nn_fiscpl, ln_hsb, nn_drown 
     58      NAMELIST/namsbc_iscpl/ nn_fiscpl, ln_iscpl_hsb, nn_drown 
    6059      !!---------------------------------------------------------------------- 
    6160      ! 
    6261      nn_fiscpl = 0 
    63       ln_hsb    = .FALSE. 
     62      ln_iscpl_hsb    = .FALSE. 
    6463      REWIND( numnam_ref )              ! Namelist namsbc_iscpl in reference namelist : Ice sheet coupling 
    6564      READ  ( numnam_ref, namsbc_iscpl, IOSTAT = ios, ERR = 901) 
    … …  
    7675         WRITE(numout,*) 'iscpl_rst:' 
    7776         WRITE(numout,*) '~~~~~~~~~' 
    78          WRITE(numout,*) ' coupling     flag (ln_iscpl )            = ', ln_iscpl 
    79          WRITE(numout,*) ' conservation flag (ln_hsb   )            = ', ln_hsb 
    80          WRITE(numout,*) ' nb of stp for cons (rn_fiscpl)           = ', nstp_iscpl 
     77         WRITE(numout,*) ' coupling     flag (ln_iscpl    )            = ', ln_iscpl 
     78         WRITE(numout,*) ' conservation flag (ln_iscpl_hsb)            = ', ln_iscpl_hsb 
     79         WRITE(numout,*) ' nb of stp for cons (rn_fiscpl  )            = ', nstp_iscpl 
    8180         IF (nstp_iscpl .NE. nn_fiscpl) WRITE(numout,*) 'W A R N I N G: nb of stp for cons has been modified & 
    8281            &                                           (larger than run length)' 

  • NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/iscplrst.F90

    r11310 r11395  
    7171      CALL iscpl_rst_interpol( ztmask_b, zumask_b, zvmask_b, zsmask_b, ze3t_b, ze3u_b, ze3v_b, zdepw_b ) 
    7272      ! 
    73       IF ( ln_hsb ) THEN      ! compute correction if conservation needed 
     73      IF ( ln_iscpl_hsb ) THEN      ! compute correction if conservation needed 
    7474         IF( iscpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'rst_iscpl : unable to allocate rst_iscpl arrays' ) 
    7575         CALL iscpl_cons(ztmask_b, zsmask_b, ze3t_b, htsc_iscpl, hdiv_iscpl, rdt_iscpl) 
    … …  
    7979      IF( ln_meshmask .AND. ln_iscpl )   CALL dom_wri 
    8080      ! 
    81       IF ( ln_hsb ) THEN 
     81      IF ( ln_iscpl_hsb ) THEN 
    8282         cfile='correction' 
    8383         cfile = TRIM( cfile ) 

  • NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfmlt.F90

    r11310 r11395  
    1 MODULE sbcisf 
     1MODULE isfmlt 
    22   !!====================================================================== 
    33   !!                       ***  MODULE  sbcisf  *** 
    4    !! Surface module :  update surface ocean boundary condition under ice 
    5    !!                   shelf 
     4   !! Surface module :  compute iceshelf melt and heat flux 
    65   !!====================================================================== 
    76   !! History :  3.2  !  2011-02  (C.Harris  ) Original code isf cav 
    87   !!            X.X  !  2006-02  (C. Wang   ) Original code bg03 
    98   !!            3.4  !  2013-03  (P. Mathiot) Merging + parametrization 
    10    !!---------------------------------------------------------------------- 
    11  
    12    !!---------------------------------------------------------------------- 
    13    !!   sbc_isf       : update sbc under ice shelf 
     9   !!            4.1  !  2019-09  (P. Mathiot) Split param/explicit ice shelf and re-organisation 
     10   !!---------------------------------------------------------------------- 
     11 
     12   !!---------------------------------------------------------------------- 
     13   !!   isfmlt       : compute iceshelf melt and heat flux 
    1414   !!---------------------------------------------------------------------- 
    1515   USE oce            ! ocean dynamics and tracers 
    … …  
    2525   USE lbclnk         ! 
    2626   USE lib_fortran    ! glob_sum 
     27   ! 
     28   USE isfrst         ! iceshelf restart 
     29   USE isftbl         ! ice shelf boundary layer 
     30   USE isfpar         ! ice shelf parametrisation 
     31   USE isfcav         ! ice shelf cavity 
     32   USE isf            ! isf variables 
    2733 
    2834   IMPLICIT NONE 
     35 
    2936   PRIVATE 
    3037 
    31    PUBLIC   sbc_isf, sbc_isf_init, sbc_isf_div, sbc_isf_alloc  ! routine called in sbcmod and divhor 
    32  
    33    ! public in order to be able to output then  
    34  
    35    REAL(wp), PUBLIC ::   rn_hisf_tbl                 !: thickness of top boundary layer [m] 
    36    INTEGER , PUBLIC ::   nn_isf                      !: flag to choose between explicit/param/specified   
    37    INTEGER , PUBLIC ::   nn_isfblk                   !: flag to choose the bulk formulation to compute the ice shelf melting 
    38    INTEGER , PUBLIC ::   nn_gammablk                 !: flag to choose how the exchange coefficient is computed 
    39    REAL(wp), PUBLIC ::   rn_gammat0                  !: temperature exchange coeficient [] 
    40    REAL(wp), PUBLIC ::   rn_gammas0                  !: salinity    exchange coeficient [] 
    41  
    42    INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   misfkt   , misfkb        !: Level of ice shelf base 
    43    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rzisf_tbl                !: depth of calving front (shallowest point) nn_isf ==2/3 
    44    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rhisf_tbl, rhisf_tbl_0   !: thickness of tbl  [m] 
    45    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   r1_hisf_tbl              !: 1/thickness of tbl 
    46    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ralpha                   !: proportion of bottom cell influenced by tbl  
    47    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   risfLeff                 !: effective length (Leff) BG03 nn_isf==2 
    48    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ttbl, stbl, utbl, vtbl   !: top boundary layer variable at T point 
    49    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qisf                     !: net heat flux from ice shelf      [W/m2] 
    50    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   risf_tsc_b, risf_tsc     !: before and now T & S isf contents [K.m/s & PSU.m/s]   
    51  
    52    LOGICAL, PUBLIC ::   l_isfcpl = .false.       !: isf recieved from oasis 
    53  
    54    REAL(wp), PUBLIC, SAVE ::   rcpisf   = 2000.0_wp     !: specific heat of ice shelf             [J/kg/K] 
    55    REAL(wp), PUBLIC, SAVE ::   rkappa   = 1.54e-6_wp    !: heat diffusivity through the ice-shelf [m2/s] 
    56    REAL(wp), PUBLIC, SAVE ::   rhoisf   = 920.0_wp      !: volumic mass of ice shelf              [kg/m3] 
    57    REAL(wp), PUBLIC, SAVE ::   tsurf    = -20.0_wp      !: air temperature on top of ice shelf    [C] 
    58    REAL(wp), PUBLIC, SAVE ::   rLfusisf = 0.334e6_wp    !: latent heat of fusion of ice shelf     [J/kg] 
    59  
    60 !: Variable used in fldread to read the forcing file (nn_isf == 4 .OR. nn_isf == 3) 
    61    CHARACTER(len=100), PUBLIC           :: cn_dirisf  = './' !: Root directory for location of ssr files 
    62    TYPE(FLD_N)       , PUBLIC           :: sn_fwfisf         !: information about the isf melting file to be read 
    63    TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_fwfisf 
    64    TYPE(FLD_N)       , PUBLIC           :: sn_rnfisf         !: information about the isf melting param.   file to be read 
    65    TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_rnfisf            
    66    TYPE(FLD_N)       , PUBLIC           :: sn_depmax_isf     !: information about the grounding line depth file to be read 
    67    TYPE(FLD_N)       , PUBLIC           :: sn_depmin_isf     !: information about the calving   line depth file to be read 
    68    TYPE(FLD_N)       , PUBLIC           :: sn_Leff_isf       !: information about the effective length     file to be read 
    69     
     38   PUBLIC   isf_mlt, isf_mlt_init  ! routine called in sbcmod and divhor 
     39 
    7040   !!---------------------------------------------------------------------- 
    7141   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    7545CONTAINS 
    7646  
    77   SUBROUTINE sbc_isf( kt ) 
     47  SUBROUTINE isf_mlt( kt ) 
    7848      !!--------------------------------------------------------------------- 
    79       !!                  ***  ROUTINE sbc_isf  *** 
    80       !! 
    81       !! ** Purpose : Compute Salt and Heat fluxes related to ice_shelf  
    82       !!              melting and freezing  
    83       !! 
    84       !! ** Method  :  4 parameterizations are available according to nn_isf  
    85       !!               nn_isf = 1 : Realistic ice_shelf formulation 
    86       !!                        2 : Beckmann & Goose parameterization 
    87       !!                        3 : Specified runoff in deptht (Mathiot & al. ) 
    88       !!                        4 : specified fwf and heat flux forcing beneath the ice shelf 
     49      !!                  ***  ROUTINE isf_mlt  *** 
     50      !! 
     51      !! ** Purpose :  
     52      !! 
     53      !! ** Method  :  
     54      !! 
    8955      !!---------------------------------------------------------------------- 
    9056      INTEGER, INTENT(in) ::   kt   ! ocean time step 
    91       ! 
    92       INTEGER ::   ji, jj, jk   ! loop index 
    93       INTEGER ::   ikt, ikb     ! local integers 
    94       REAL(wp), DIMENSION(jpi,jpj) ::   zt_frz, zdep   ! freezing temperature (zt_frz) at depth (zdep)  
    95       REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   zqhcisf2d 
    96       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   zfwfisf3d, zqhcisf3d, zqlatisf3d 
    9757      !!--------------------------------------------------------------------- 
    9858      ! 
    99       IF( MOD( kt-1, nn_fsbc) == 0 ) THEN    ! compute salt and heat flux 
    100          ! 
    101          SELECT CASE ( nn_isf ) 
    102          CASE ( 1 )    ! realistic ice shelf formulation 
    103             ! compute T/S/U/V for the top boundary layer 
    104             CALL sbc_isf_tbl(tsn(:,:,:,jp_tem),ttbl(:,:),'T') 
    105             CALL sbc_isf_tbl(tsn(:,:,:,jp_sal),stbl(:,:),'T') 
    106             CALL sbc_isf_tbl(un(:,:,:)        ,utbl(:,:),'U') 
    107             CALL sbc_isf_tbl(vn(:,:,:)        ,vtbl(:,:),'V') 
    108             ! iom print 
    109             CALL iom_put('ttbl',ttbl(:,:)) 
    110             CALL iom_put('stbl',stbl(:,:)) 
    111             CALL iom_put('utbl',utbl(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)) 
    112             CALL iom_put('vtbl',vtbl(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)) 
    113             ! compute fwf and heat flux 
    114             ! compute fwf and heat flux 
    115             IF( .NOT.l_isfcpl ) THEN    ;   CALL sbc_isf_cav (kt) 
    116             ELSE                        ;   qisf(:,:)  = fwfisf(:,:) * rLfusisf  ! heat        flux 
    117             ENDIF 
    118             ! 
    119          CASE ( 2 )    ! Beckmann and Goosse parametrisation  
    120             stbl(:,:)   = soce 
    121             CALL sbc_isf_bg03(kt) 
    122             ! 
    123          CASE ( 3 )    ! specified runoff in depth (Mathiot et al., XXXX in preparation) 
    124             ! specified runoff in depth (Mathiot et al., XXXX in preparation) 
    125             IF( .NOT.l_isfcpl ) THEN 
    126                CALL fld_read ( kt, nn_fsbc, sf_rnfisf   ) 
    127                fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1)         ! fresh water flux from the isf (fwfisf <0 mean melting)  
    128             ENDIF 
    129             qisf(:,:)   = fwfisf(:,:) * rLfusisf             ! heat flux 
    130             stbl(:,:)   = soce 
    131             ! 
    132          CASE ( 4 )    ! specified fwf and heat flux forcing beneath the ice shelf 
    133             !          ! specified fwf and heat flux forcing beneath the ice shelf 
    134             IF( .NOT.l_isfcpl ) THEN 
    135                CALL fld_read ( kt, nn_fsbc, sf_fwfisf   ) 
    136                !CALL fld_read ( kt, nn_fsbc, sf_qisf   ) 
    137                fwfisf(:,:) = -sf_fwfisf(1)%fnow(:,:,1)            ! fwf 
    138             ENDIF 
    139             qisf(:,:)   = fwfisf(:,:) * rLfusisf               ! heat flux 
    140             stbl(:,:)   = soce 
    141             ! 
    142          END SELECT 
    143  
    144          ! compute tsc due to isf 
    145          ! isf melting implemented as a volume flux and we assume that melt water is at 0 PSU. 
    146          ! WARNING water add at temp = 0C, need to add a correction term (fwfisf * tfreez / rau0). 
    147          ! compute freezing point beneath ice shelf (or top cell if nn_isf = 3) 
    148          DO jj = 1,jpj 
    149             DO ji = 1,jpi 
    150                zdep(ji,jj)=gdepw_n(ji,jj,misfkt(ji,jj)) 
    151             END DO 
    152          END DO 
    153          CALL eos_fzp( stbl(:,:), zt_frz(:,:), zdep(:,:) ) 
    154           
    155          risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp - fwfisf(:,:) * zt_frz(:,:) * r1_rau0 ! 
    156          risf_tsc(:,:,jp_sal) = 0.0_wp 
    157  
    158          ! lbclnk 
    159          CALL lbc_lnk_multi( 'sbcisf', risf_tsc(:,:,jp_tem), 'T', 1., risf_tsc(:,:,jp_sal), 'T', 1., fwfisf,'T', 1., qisf, 'T', 1.) 
    160          ! output 
    161          IF( iom_use('iceshelf_cea') )   CALL iom_put( 'iceshelf_cea', -fwfisf(:,:)                      )   ! isf mass flux 
    162          IF( iom_use('hflx_isf_cea') )   CALL iom_put( 'hflx_isf_cea', risf_tsc(:,:,jp_tem) * rau0 * rcp )   ! isf sensible+latent heat (W/m2) 
    163          IF( iom_use('qlatisf' ) )       CALL iom_put( 'qlatisf'     , qisf(:,:)                         )   ! isf latent heat 
    164          IF( iom_use('fwfisf'  ) )       CALL iom_put( 'fwfisf'      , fwfisf(:,:)                       )   ! isf mass flux (opposite sign) 
    165  
    166          ! Diagnostics 
    167          IF( iom_use('fwfisf3d') .OR. iom_use('qlatisf3d') .OR. iom_use('qhcisf3d') .OR. iom_use('qhcisf')) THEN 
    168             ALLOCATE( zfwfisf3d(jpi,jpj,jpk) , zqhcisf3d(jpi,jpj,jpk) , zqlatisf3d(jpi,jpj,jpk) ) 
    169             ALLOCATE( zqhcisf2d(jpi,jpj) ) 
    170             ! 
    171             zfwfisf3d (:,:,:) = 0._wp                         ! 3d ice shelf melting (kg/m2/s) 
    172             zqhcisf3d (:,:,:) = 0._wp                         ! 3d heat content flux (W/m2) 
    173             zqlatisf3d(:,:,:) = 0._wp                         ! 3d ice shelf melting latent heat flux (W/m2) 
    174             zqhcisf2d (:,:)   = fwfisf(:,:) * zt_frz * rcp    ! 2d heat content flux (W/m2) 
    175             ! 
    176             DO jj = 1,jpj 
    177                DO ji = 1,jpi 
    178                   ikt = misfkt(ji,jj) 
    179                   ikb = misfkb(ji,jj) 
    180                   DO jk = ikt, ikb - 1 
    181                      zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf   (ji,jj) * r1_hisf_tbl(ji,jj) * e3t_n(ji,jj,jk) 
    182                      zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj) * e3t_n(ji,jj,jk) 
    183                      zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf     (ji,jj) * r1_hisf_tbl(ji,jj) * e3t_n(ji,jj,jk) 
    184                   END DO 
    185                   zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf   (ji,jj) * r1_hisf_tbl(ji,jj)   &  
    186                      &                                                                   * ralpha(ji,jj) * e3t_n(ji,jj,jk) 
    187                   zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj)   &  
    188                      &                                                                   * ralpha(ji,jj) * e3t_n(ji,jj,jk) 
    189                   zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf     (ji,jj) * r1_hisf_tbl(ji,jj)   &   
    190                      &                                                                   * ralpha(ji,jj) * e3t_n(ji,jj,jk) 
    191                END DO 
    192             END DO 
    193             ! 
    194             CALL iom_put('fwfisf3d' , zfwfisf3d (:,:,:)) 
    195             CALL iom_put('qlatisf3d', zqlatisf3d(:,:,:)) 
    196             CALL iom_put('qhcisf3d' , zqhcisf3d (:,:,:)) 
    197             CALL iom_put('qhcisf'   , zqhcisf2d (:,:  )) 
    198             ! 
    199             DEALLOCATE( zfwfisf3d, zqhcisf3d, zqlatisf3d ) 
    200             DEALLOCATE( zqhcisf2d ) 
    201          ENDIF 
    202          ! 
    203       ENDIF 
    204  
    205       IF( kt == nit000 ) THEN                         !   set the forcing field at nit000 - 1    ! 
    206          IF( ln_rstart .AND.    &                     ! Restart: read in restart file 
    207             &   iom_varid( numror, 'fwf_isf_b', ldstop = .FALSE. ) > 0 ) THEN 
    208             IF(lwp) WRITE(numout,*) '          nit000-1 isf tracer content forcing fields read in the restart file' 
    209             CALL iom_get( numror, jpdom_autoglo, 'fwf_isf_b', fwfisf_b(:,:)         , ldxios = lrxios )   ! before salt content isf_tsc trend 
    210             CALL iom_get( numror, jpdom_autoglo, 'isf_sc_b' , risf_tsc_b(:,:,jp_sal), ldxios = lrxios )   ! before salt content isf_tsc trend 
    211             CALL iom_get( numror, jpdom_autoglo, 'isf_hc_b' , risf_tsc_b(:,:,jp_tem), ldxios = lrxios )   ! before salt content isf_tsc trend 
    212          ELSE 
    213             fwfisf_b(:,:)    = fwfisf(:,:) 
    214             risf_tsc_b(:,:,:)= risf_tsc(:,:,:) 
    215          ENDIF 
    216       ENDIF 
     59      IF ( ln_isfcav_mlt ) THEN 
     60         ! 
     61         ! before time step  
     62         IF ( kt /= nit000 ) THEN  
     63            risf_cav_tsc_b (:,:,:) = risf_cav_tsc (:,:,:) 
     64            fwfisf_cav_b(:,:)      = fwfisf_cav(:,:) 
     65         END IF 
     66         ! 
     67         ! compute tbl lvl/h 
     68         CALL isf_tbl_lvl(ht_n, e3t_n, misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav) 
     69         ! 
     70         ! compute ice shelf melt 
     71         CALL isf_cav( kt, risf_cav_tsc, fwfisf_cav) 
     72         ! 
     73         ! write restart variables (qoceisf, qhcisf, fwfisf for now and before) 
     74         IF (lrst_oce) CALL isfrst_write(kt, 'cav', risf_cav_tsc, fwfisf_cav) 
     75         ! 
     76      END IF 
    21777      !  
    218       IF( lrst_oce ) THEN 
    219          IF(lwp) WRITE(numout,*) 
    220          IF(lwp) WRITE(numout,*) 'sbc : isf surface tracer content forcing fields written in ocean restart file ',   & 
    221             &                    'at it= ', kt,' date= ', ndastp 
    222          IF(lwp) WRITE(numout,*) '~~~~' 
    223          IF( lwxios ) CALL iom_swap(      cwxios_context          ) 
    224          CALL iom_rstput( kt, nitrst, numrow, 'fwf_isf_b', fwfisf(:,:)         , ldxios = lwxios ) 
    225          CALL iom_rstput( kt, nitrst, numrow, 'isf_hc_b' , risf_tsc(:,:,jp_tem), ldxios = lwxios ) 
    226          CALL iom_rstput( kt, nitrst, numrow, 'isf_sc_b' , risf_tsc(:,:,jp_sal), ldxios = lwxios ) 
    227          IF( lwxios ) CALL iom_swap(      cxios_context          ) 
    228       ENDIF 
    229       ! 
    230    END SUBROUTINE sbc_isf 
    231  
    232  
    233    INTEGER FUNCTION sbc_isf_alloc() 
    234       !!---------------------------------------------------------------------- 
    235       !!               ***  FUNCTION sbc_isf_rnf_alloc  *** 
    236       !!---------------------------------------------------------------------- 
    237       sbc_isf_alloc = 0       ! set to zero if no array to be allocated 
    238       IF( .NOT. ALLOCATED( qisf ) ) THEN 
    239          ALLOCATE(  risf_tsc(jpi,jpj,jpts), risf_tsc_b(jpi,jpj,jpts), qisf(jpi,jpj)   , & 
    240                &    rhisf_tbl(jpi,jpj)    , r1_hisf_tbl(jpi,jpj), rzisf_tbl(jpi,jpj)  , & 
    241                &    ttbl(jpi,jpj)         , stbl(jpi,jpj)       , utbl(jpi,jpj)       , & 
    242                &    vtbl(jpi, jpj)        , risfLeff(jpi,jpj)   , rhisf_tbl_0(jpi,jpj), & 
    243                &    ralpha(jpi,jpj)       , misfkt(jpi,jpj)     , misfkb(jpi,jpj)     , & 
    244                &    STAT= sbc_isf_alloc ) 
    245          ! 
    246          CALL mpp_sum ( 'sbcisf', sbc_isf_alloc ) 
    247          IF( sbc_isf_alloc /= 0 )   CALL ctl_stop( 'STOP', 'sbc_isf_alloc: failed to allocate arrays.' ) 
    248          ! 
    249       ENDIF 
    250    END FUNCTION 
    251  
    252  
    253   SUBROUTINE sbc_isf_init 
     78      IF ( ln_isfpar_mlt ) THEN 
     79         ! 
     80         ! before time step  
     81         IF ( kt /= nit000 ) THEN  
     82            risf_par_tsc_b (:,:) = risf_par_tsc (:,:) 
     83            fwfisf_par_b(:,:)    = fwfisf_par(:,:) 
     84         END IF 
     85         ! 
     86         ! compute misfkb, rhisf_tbl, rfrac (deepest level, thickness, fraction of deepest cell affected by tbl) 
     87         CALL isf_tbl_lvl(ht_n, e3t_n, misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par) 
     88         ! 
     89         ! compute ice shelf melt 
     90         CALL isf_par( kt, risf_par_tsc, fwfisf_par) 
     91         ! 
     92         ! write restart variables (qoceisf, qhcisf, fwfisf for now and before) 
     93         IF (lrst_oce) CALL isfrst_write(kt, 'par', risf_par_tsc, fwfisf_par) 
     94         ! 
     95      END IF 
     96      ! 
     97   END SUBROUTINE isf_mlt 
     98 
     99   SUBROUTINE isf_mlt_init 
    254100      !!--------------------------------------------------------------------- 
    255       !!                  ***  ROUTINE sbc_isf_init  *** 
    256       !! 
    257       !! ** Purpose : Initialisation of variables for iceshelf fluxes formulation 
    258       !! 
    259       !! ** Method  :  4 parameterizations are available according to nn_isf  
    260       !!               nn_isf = 1 : Realistic ice_shelf formulation 
    261       !!                        2 : Beckmann & Goose parameterization 
    262       !!                        3 : Specified runoff in deptht (Mathiot & al. ) 
    263       !!                        4 : specified fwf and heat flux forcing beneath the ice shelf 
    264       !!---------------------------------------------------------------------- 
    265       INTEGER               :: ji, jj, jk           ! loop index 
    266       INTEGER               :: ik                   ! current level index 
    267       INTEGER               :: ikt, ikb             ! top and bottom level of the isf boundary layer 
     101      !!                  ***  ROUTINE isfmlt_init  *** 
     102      !! 
     103      !! ** Purpose : 
     104      !! 
     105      !! ** Method  : 
     106      !! 
     107      !!---------------------------------------------------------------------- 
    268108      INTEGER               :: inum, ierror 
    269109      INTEGER               :: ios                  ! Local integer output status for namelist read 
    270       REAL(wp)              :: zhk 
    271       CHARACTER(len=256)    :: cvarzisf, cvarhisf   ! name for isf file 
    272       CHARACTER(LEN=32 )    :: cvarLeff             ! variable name for efficient Length scale 
    273       !!---------------------------------------------------------------------- 
    274       NAMELIST/namsbc_isf/ nn_isfblk, rn_hisf_tbl, rn_gammat0, rn_gammas0, nn_gammablk, nn_isf, & 
    275                          & sn_fwfisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf 
     110      INTEGER               :: ikt, ikb 
     111      INTEGER               :: ji, jj 
     112      !!---------------------------------------------------------------------- 
     113      NAMELIST/namisf/ ln_isfcav_mlt, cn_isfcav_mlt, cn_gammablk, rn_gammat0, rn_gammas0, rn_htbl, sn_isfcav_fwf,  & 
     114         &             ln_isfpar_mlt, cn_isfpar_mlt, sn_isfpar_fwf, sn_isfpar_zmin, sn_isfpar_zmax, sn_isfpar_Leff 
    276115      !!---------------------------------------------------------------------- 
    277116 
    278117      REWIND( numnam_ref )              ! Namelist namsbc_rnf in reference namelist : Runoffs  
    279       READ  ( numnam_ref, namsbc_isf, IOSTAT = ios, ERR = 901) 
    280 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc_isf in reference namelist', lwp ) 
     118      READ  ( numnam_ref, namisf, IOSTAT = ios, ERR = 901) 
     119901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namisf in reference namelist', lwp ) 
    281120 
    282121      REWIND( numnam_cfg )              ! Namelist namsbc_rnf in configuration namelist : Runoffs 
    283       READ  ( numnam_cfg, namsbc_isf, IOSTAT = ios, ERR = 902 ) 
    284 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namsbc_isf in configuration namelist', lwp ) 
    285       IF(lwm) WRITE ( numond, namsbc_isf ) 
    286  
     122      READ  ( numnam_cfg, namisf, IOSTAT = ios, ERR = 902 ) 
     123902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namisf in configuration namelist', lwp ) 
     124      IF(lwm) WRITE ( numond, namisf ) 
     125      ! 
    287126      IF(lwp) WRITE(numout,*) 
    288       IF(lwp) WRITE(numout,*) 'sbc_isf_init : heat flux of the ice shelf' 
     127      IF(lwp) WRITE(numout,*) 'isf_init : ice shelf initialisation' 
    289128      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
    290       IF(lwp) WRITE(numout,*) '   Namelist namsbc_isf :' 
    291       IF(lwp) WRITE(numout,*) '      type ice shelf melting/freezing         nn_isf      = ', nn_isf 
    292       IF(lwp) WRITE(numout,*) '      bulk formulation (nn_isf=1 only)        nn_isfblk   = ', nn_isfblk 
    293       IF(lwp) WRITE(numout,*) '      thickness of the top boundary layer     rn_hisf_tbl = ', rn_hisf_tbl 
    294       IF(lwp) WRITE(numout,*) '      gamma formulation                       nn_gammablk = ', nn_gammablk  
    295       IF(lwp) WRITE(numout,*) '      gammat coefficient                      rn_gammat0  = ', rn_gammat0   
    296       IF(lwp) WRITE(numout,*) '      gammas coefficient                      rn_gammas0  = ', rn_gammas0   
    297       IF(lwp) WRITE(numout,*) '      top drag coef. used (from namdrg_top)   rn_Cd0      = ', r_Cdmin_top  
    298  
    299  
    300                            !  1 = presence of ISF    2 = bg03 parametrisation  
    301                            !  3 = rnf file for isf   4 = ISF fwf specified 
    302                            !  option 1 and 4 need ln_isfcav = .true. (domzgr) 
    303       ! 
    304       ! Allocate public variable 
    305       IF ( sbc_isf_alloc()  /= 0 )         CALL ctl_stop( 'STOP', 'sbc_isf : unable to allocate arrays' ) 
     129      IF(lwp) WRITE(numout,*) '   Namelist namisf :' 
     130      IF(lwp) WRITE(numout,*) '      melt inside the cavity                  ln_isfcav_mlt   = ', ln_isfcav_mlt 
     131      IF ( ln_isfcav ) THEN 
     132         IF(lwp) WRITE(numout,*) '         melt formulation                        cn_isfcav_mlt   = ', cn_isfcav_mlt 
     133         IF(lwp) WRITE(numout,*) '         thickness of the top boundary layer     rn_htbl     = ', rn_htbl 
     134         IF(lwp) WRITE(numout,*) '         gamma formulation                       cn_gammablk = ', cn_gammablk  
     135         IF ( TRIM(cn_gammablk) .NE. 'spe' ) THEN  
     136            IF(lwp) WRITE(numout,*) '            gammat coefficient                      rn_gammat0  = ', rn_gammat0   
     137            IF(lwp) WRITE(numout,*) '            gammas coefficient                      rn_gammas0  = ', rn_gammas0   
     138            IF(lwp) WRITE(numout,*) '            top drag coef. used (from namdrg_top)   rn_Cd0      = ', r_Cdmin_top  
     139         END IF 
     140      END IF 
     141      IF(lwp) WRITE(numout,*) '' 
     142      IF(lwp) WRITE(numout,*) '      ice shelf melt parametrisation          ln_isfpar_mlt    = ', ln_isfpar_mlt 
     143      IF ( ln_isfpar_mlt ) THEN 
     144         IF(lwp) WRITE(numout,*) '         isf parametrisation formulation         cn_isfpar_mlt   = ', cn_isfpar_mlt 
     145      END IF 
     146      IF(lwp) WRITE(numout,*) '' 
     147      ! 
     148      ! Allocate public array 
     149      CALL isf_alloc() 
    306150      ! 
    307151      ! initialisation 
    308       qisf    (:,:)    = 0._wp   ;   fwfisf  (:,:) = 0._wp 
    309       risf_tsc(:,:,:)  = 0._wp   ;   fwfisf_b(:,:) = 0._wp 
    310       ! 
    311       ! define isf tbl tickness, top and bottom indice 
    312       SELECT CASE ( nn_isf ) 
    313       CASE ( 1 )  
    314          IF(lwp) WRITE(numout,*) 
    315          IF(lwp) WRITE(numout,*) '      ==>>>   presence of under iceshelf seas (nn_isf = 1)' 
    316          rhisf_tbl(:,:) = rn_hisf_tbl 
    317          misfkt   (:,:) = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv 
    318          ! 
    319       CASE ( 2 , 3 ) 
    320          IF( .NOT.l_isfcpl ) THEN 
    321              ALLOCATE( sf_rnfisf(1), STAT=ierror ) 
    322              ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) ) 
    323              CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 
    324           ENDIF 
    325           !  read effective lenght (BG03) 
    326           IF( nn_isf == 2 ) THEN 
    327             IF(lwp) WRITE(numout,*) 
    328             IF(lwp) WRITE(numout,*) '      ==>>>   bg03 parametrisation (nn_isf = 2)' 
    329             CALL iom_open( sn_Leff_isf%clname, inum ) 
    330             cvarLeff = TRIM(sn_Leff_isf%clvar) 
    331             CALL iom_get( inum, jpdom_data, cvarLeff, risfLeff , 1) 
    332             CALL iom_close(inum) 
    333             ! 
    334             risfLeff = risfLeff*1000.0_wp           !: convertion in m 
    335          ELSE 
    336             IF(lwp) WRITE(numout,*) 
    337             IF(lwp) WRITE(numout,*) '      ==>>>   rnf file for isf (nn_isf = 3)' 
    338          ENDIF 
    339          ! read depth of the top and bottom of the isf top boundary layer (in this case, isf front depth and grounding line depth) 
    340          CALL iom_open( sn_depmax_isf%clname, inum ) 
    341          cvarhisf = TRIM(sn_depmax_isf%clvar) 
    342          CALL iom_get( inum, jpdom_data, cvarhisf, rhisf_tbl, 1) !: depth of deepest point of the ice shelf base 
    343          CALL iom_close(inum) 
    344          ! 
    345          CALL iom_open( sn_depmin_isf%clname, inum ) 
    346          cvarzisf = TRIM(sn_depmin_isf%clvar) 
    347          CALL iom_get( inum, jpdom_data, cvarzisf, rzisf_tbl, 1) !: depth of shallowest point of the ice shelves base 
    348          CALL iom_close(inum) 
    349          ! 
    350          rhisf_tbl(:,:) = rhisf_tbl(:,:) - rzisf_tbl(:,:)        !: tickness isf boundary layer 
    351  
    352          !! compute first level of the top boundary layer 
    353          DO ji = 1, jpi 
    354             DO jj = 1, jpj 
    355                 ik = 2 
    356 !!gm potential bug: use gdepw_0 not _n 
    357                 DO WHILE ( ik <= mbkt(ji,jj) .AND. gdepw_n(ji,jj,ik) < rzisf_tbl(ji,jj) ) ;  ik = ik + 1 ;  END DO 
    358                 misfkt(ji,jj) = ik-1 
    359             END DO 
    360          END DO 
    361          ! 
    362       CASE ( 4 )  
    363          IF(lwp) WRITE(numout,*) 
    364          IF(lwp) WRITE(numout,*) '      ==>>>   specified fresh water flux in ISF (nn_isf = 4)' 
    365          ! as in nn_isf == 1 
    366          rhisf_tbl(:,:) = rn_hisf_tbl 
    367          misfkt   (:,:) = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv 
    368          ! 
    369          ! load variable used in fldread (use for temporal interpolation of isf fwf forcing) 
    370          IF( .NOT.l_isfcpl ) THEN 
    371            ALLOCATE( sf_fwfisf(1), STAT=ierror ) 
    372            ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) ) 
    373            CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 
    374          ENDIF 
    375          ! 
    376       CASE DEFAULT 
    377          CALL ctl_stop( 'sbc_isf_init: wrong value of nn_isf' ) 
    378       END SELECT 
    379           
    380       rhisf_tbl_0(:,:) = rhisf_tbl(:,:) 
    381  
    382       ! compute bottom level of isf tbl and thickness of tbl below the ice shelf 
     152      ! cav 
     153      misfkt_cav(:,:)     = mikt(:,:) ; misfkb_cav(:,:)       = mikt(:,:) 
     154      fwfisf_cav(:,:)     = 0.0_wp    ; fwfisf_cav_b(:,:)     = 0.0_wp 
     155      rhisf_tbl_cav(:,:)  = 1e-20     ; rfrac_tbl_cav(:,:)    = 0.0_wp 
     156      risf_cav_tsc(:,:,:) = 0.0_wp    ; risf_cav_tsc_b(:,:,:) = 0.0_wp 
     157      ! 
     158      mskisf_cav(:,:) = (1._wp - tmask(:,:,1)) * ssmask(:,:) 
     159      ! 
     160      ! par 
     161      misfkt_par(:,:)     = 1         ; misfkb_par(:,:)       = 1          
     162      fwfisf_par(:,:)     = 0.0_wp    ; fwfisf_par_b(:,:)     = 0.0_wp 
     163      rhisf_tbl_par(:,:)  = 1e-20     ; rfrac_tbl_par(:,:)    = 0.0_wp 
     164      risf_par_tsc(:,:,:) = 0.0_wp    ; risf_par_tsc_b(:,:,:) = 0.0_wp 
     165      ! 
     166      mskisf_par(:,:) = 0 
     167      ! 
     168      ! initialisation useful variable 
     169      r1_Lfusisf =  1._wp / rLfusisf 
     170      ! 
    383171      DO jj = 1,jpj 
    384172         DO ji = 1,jpi 
    385             ikt = misfkt(ji,jj) 
    386             ikb = misfkt(ji,jj) 
    387             ! thickness of boundary layer at least the top level thickness 
    388             rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3t_n(ji,jj,ikt)) 
    389  
    390             ! determine the deepest level influenced by the boundary layer 
    391             DO jk = ikt+1, mbkt(ji,jj) 
    392                IF( (SUM(e3t_n(ji,jj,ikt:jk-1)) < rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) )   ikb = jk 
    393             END DO 
    394             rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness. 
    395             misfkb(ji,jj) = ikb                                                   ! last wet level of the tbl 
    396             r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj) 
    397  
    398             zhk           = SUM( e3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1 
    399             ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / e3t_n(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer 
     173            ikt = mikt(ji,jj) 
     174            ikb = mbkt(ji,jj) 
     175            bathy  (ji,jj) = gdepw_0(ji,jj,ikb+1) 
     176            risfdep(ji,jj) = gdepw_0(ji,jj,ikt  ) 
    400177         END DO 
    401178      END DO 
    402  
    403       IF( lwxios ) THEN 
    404           CALL iom_set_rstw_var_active('fwf_isf_b') 
    405           CALL iom_set_rstw_var_active('isf_hc_b') 
    406           CALL iom_set_rstw_var_active('isf_sc_b') 
    407       ENDIF 
    408  
    409  
    410   END SUBROUTINE sbc_isf_init 
    411  
    412  
    413   SUBROUTINE sbc_isf_bg03(kt) 
    414       !!--------------------------------------------------------------------- 
    415       !!                  ***  ROUTINE sbc_isf_bg03  *** 
    416       !! 
    417       !! ** Purpose : add net heat and fresh water flux from ice shelf melting 
    418       !!          into the adjacent ocean 
    419       !! 
    420       !! ** Method  :   See reference 
    421       !! 
    422       !! ** Reference : Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean 
    423       !!         interaction for climate models", Ocean Modelling 5(2003) 157-170. 
    424       !!         (hereafter BG) 
    425       !! History :  06-02  (C. Wang) Original code 
    426       !!---------------------------------------------------------------------- 
    427       INTEGER, INTENT ( in ) :: kt 
    428       ! 
    429       INTEGER  :: ji, jj, jk ! dummy loop index 
    430       INTEGER  :: ik         ! current level 
    431       REAL(wp) :: zt_sum     ! sum of the temperature between 200m and 600m 
    432       REAL(wp) :: zt_ave     ! averaged temperature between 200m and 600m 
    433       REAL(wp) :: zt_frz     ! freezing point temperature at depth z 
    434       REAL(wp) :: zpress     ! pressure to compute the freezing point in depth 
    435       !!---------------------------------------------------------------------- 
    436       ! 
    437       DO ji = 1, jpi 
    438          DO jj = 1, jpj 
    439             ik = misfkt(ji,jj) 
    440             !! Initialize arrays to 0 (each step) 
    441             zt_sum = 0.e0_wp 
    442             IF ( ik > 1 ) THEN 
    443                ! 1. -----------the average temperature between 200m and 600m --------------------- 
    444                DO jk = misfkt(ji,jj),misfkb(ji,jj) 
    445                   ! Calculate freezing temperature 
    446                   zpress = grav*rau0*gdept_n(ji,jj,ik)*1.e-04 
    447                   CALL eos_fzp(stbl(ji,jj), zt_frz, zpress)  
    448                   zt_sum = zt_sum + (tsn(ji,jj,jk,jp_tem)-zt_frz) * e3t_n(ji,jj,jk) * tmask(ji,jj,jk)  ! sum temp 
    449                END DO 
    450                zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value 
    451                ! 2. ------------Net heat flux and fresh water flux due to the ice shelf 
    452                ! For those corresponding to zonal boundary     
    453                qisf(ji,jj) = - rau0 * rcp * rn_gammat0 * risfLeff(ji,jj) * e1t(ji,jj) * zt_ave  & 
    454                            & * r1_e1e2t(ji,jj) * tmask(ji,jj,jk) 
    455               
    456                fwfisf(ji,jj) = qisf(ji,jj) / rLfusisf          !fresh water flux kg/(m2s)                   
    457                fwfisf(ji,jj) = fwfisf(ji,jj) * ( soce / stbl(ji,jj) ) 
    458                !add to salinity trend 
    459             ELSE 
    460                qisf(ji,jj) = 0._wp   ;   fwfisf(ji,jj) = 0._wp 
    461             END IF 
    462          END DO 
    463       END DO 
    464       ! 
    465   END SUBROUTINE sbc_isf_bg03 
    466  
    467  
    468   SUBROUTINE sbc_isf_cav( kt ) 
    469       !!--------------------------------------------------------------------- 
    470       !!                     ***  ROUTINE sbc_isf_cav  *** 
    471       !! 
    472       !! ** Purpose :   handle surface boundary condition under ice shelf 
    473       !! 
    474       !! ** Method  : - 
    475       !! 
    476       !! ** Action  :   utau, vtau : remain unchanged 
    477       !!                taum, wndm : remain unchanged 
    478       !!                qns        : update heat flux below ice shelf 
    479       !!                emp, emps  : update freshwater flux below ice shelf 
    480       !!--------------------------------------------------------------------- 
    481       INTEGER, INTENT(in) ::   kt   ! ocean time step 
    482       ! 
    483       INTEGER  ::   ji, jj     ! dummy loop indices 
    484       INTEGER  ::   nit 
    485       LOGICAL  ::   lit 
    486       REAL(wp) ::   zlamb1, zlamb2, zlamb3 
    487       REAL(wp) ::   zeps1,zeps2,zeps3,zeps4,zeps6,zeps7 
    488       REAL(wp) ::   zaqe,zbqe,zcqe,zaqer,zdis,zsfrz,zcfac 
    489       REAL(wp) ::   zeps = 1.e-20_wp         
    490       REAL(wp) ::   zerr 
    491       REAL(wp), DIMENSION(jpi,jpj) ::   zfrz 
    492       REAL(wp), DIMENSION(jpi,jpj) ::   zgammat, zgammas  
    493       REAL(wp), DIMENSION(jpi,jpj) ::   zfwflx, zhtflx, zhtflx_b 
    494       !!--------------------------------------------------------------------- 
    495       ! 
    496       ! coeficient for linearisation of potential tfreez 
    497       ! Crude approximation for pressure (but commonly used) 
    498       IF ( l_useCT ) THEN   ! linearisation from Jourdain et al. (2017) 
    499          zlamb1 =-0.0564_wp 
    500          zlamb2 = 0.0773_wp 
    501          zlamb3 =-7.8633e-8 * grav * rau0 
    502       ELSE                  ! linearisation from table 4 (Asay-Davis et al., 2015) 
    503          zlamb1 =-0.0573_wp 
    504          zlamb2 = 0.0832_wp 
    505          zlamb3 =-7.53e-8 * grav * rau0 
    506       ENDIF 
    507       ! 
    508       ! initialisation 
    509       zgammat(:,:) = rn_gammat0 ; zgammas (:,:) = rn_gammas0 
    510       zhtflx (:,:) = 0.0_wp     ; zhtflx_b(:,:) = 0.0_wp     
    511       zfwflx (:,:) = 0.0_wp 
    512  
    513       ! compute ice shelf melting 
    514       nit = 1 ; lit = .TRUE. 
    515       DO WHILE ( lit )    ! maybe just a constant number of iteration as in blk_core is fine 
    516          SELECT CASE ( nn_isfblk ) 
    517          CASE ( 1 )   !  ISOMIP formulation (2 equations) for volume flux (Hunter et al., 2006) 
    518             ! Calculate freezing temperature 
    519             CALL eos_fzp( stbl(:,:), zfrz(:,:), risfdep(:,:) ) 
    520  
    521             ! compute gammat every where (2d) 
    522             CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx) 
    523              
    524             ! compute upward heat flux zhtflx and upward water flux zwflx 
    525             DO jj = 1, jpj 
    526                DO ji = 1, jpi 
    527                   zhtflx(ji,jj) =   zgammat(ji,jj)*rcp*rau0*(ttbl(ji,jj)-zfrz(ji,jj)) 
    528                   zfwflx(ji,jj) = - zhtflx(ji,jj)/rLfusisf 
    529                END DO 
    530             END DO 
    531  
    532             ! Compute heat flux and upward fresh water flux 
    533             qisf  (:,:) = - zhtflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 
    534             fwfisf(:,:) =   zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 
    535  
    536          CASE ( 2 )  ! ISOMIP+ formulation (3 equations) for volume flux (Asay-Davis et al., 2015) 
    537             ! compute gammat every where (2d) 
    538             CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx) 
    539  
    540             ! compute upward heat flux zhtflx and upward water flux zwflx 
    541             ! Resolution of a 2d equation from equation 21, 22 and 23 to find Sb (Asay-Davis et al., 2015) 
    542             DO jj = 1, jpj 
    543                DO ji = 1, jpi 
    544                   ! compute coeficient to solve the 2nd order equation 
    545                   zeps1 = rcp*rau0*zgammat(ji,jj) 
    546                   zeps2 = rLfusisf*rau0*zgammas(ji,jj) 
    547                   zeps3 = rhoisf*rcpisf*rkappa/MAX(risfdep(ji,jj),zeps) 
    548                   zeps4 = zlamb2+zlamb3*risfdep(ji,jj) 
    549                   zeps6 = zeps4-ttbl(ji,jj) 
    550                   zeps7 = zeps4-tsurf 
    551                   zaqe  = zlamb1 * (zeps1 + zeps3) 
    552                   zaqer = 0.5_wp/MIN(zaqe,-zeps) 
    553                   zbqe  = zeps1*zeps6+zeps3*zeps7-zeps2 
    554                   zcqe  = zeps2*stbl(ji,jj) 
    555                   zdis  = zbqe*zbqe-4.0_wp*zaqe*zcqe                
    556  
    557                   ! Presumably zdis can never be negative because gammas is very small compared to gammat 
    558                   ! compute s freeze 
    559                   zsfrz=(-zbqe-SQRT(zdis))*zaqer 
    560                   IF ( zsfrz < 0.0_wp ) zsfrz=(-zbqe+SQRT(zdis))*zaqer 
    561  
    562                   ! compute t freeze (eq. 22) 
    563                   zfrz(ji,jj)=zeps4+zlamb1*zsfrz 
    564    
    565                   ! zfwflx is upward water flux 
    566                   ! zhtflx is upward heat flux (out of ocean) 
    567                   ! compute the upward water and heat flux (eq. 28 and eq. 29) 
    568                   zfwflx(ji,jj) = rau0 * zgammas(ji,jj) * (zsfrz-stbl(ji,jj)) / MAX(zsfrz,zeps) 
    569                   zhtflx(ji,jj) = zgammat(ji,jj) * rau0 * rcp * (ttbl(ji,jj) - zfrz(ji,jj) )  
    570                END DO 
    571             END DO 
    572  
    573             ! compute heat and water flux 
    574             qisf  (:,:) = - zhtflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 
    575             fwfisf(:,:) =   zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 
    576  
    577          END SELECT 
    578  
    579          ! define if we need to iterate (nn_gammablk 0/1 do not need iteration) 
    580          IF ( nn_gammablk <  2 ) THEN ; lit = .FALSE. 
    581          ELSE                            
    582             ! check total number of iteration 
    583             IF (nit >= 100) THEN ; CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' ) 
    584             ELSE                 ; nit = nit + 1 
    585             END IF 
    586  
    587             ! compute error between 2 iterations 
    588             ! if needed save gammat and compute zhtflx_b for next iteration 
    589             zerr = MAXVAL(ABS(zhtflx-zhtflx_b)) 
    590             IF ( zerr <= 0.01_wp ) THEN ; lit = .FALSE. 
    591             ELSE                        ; zhtflx_b(:,:) = zhtflx(:,:) 
    592             END IF 
    593          END IF 
    594       END DO 
    595       ! 
    596       CALL iom_put('isfgammat', zgammat) 
    597       CALL iom_put('isfgammas', zgammas) 
    598       !  
    599    END SUBROUTINE sbc_isf_cav 
    600  
    601  
    602    SUBROUTINE sbc_isf_gammats(pgt, pgs, pqhisf, pqwisf ) 
    603       !!---------------------------------------------------------------------- 
    604       !! ** Purpose    : compute the coefficient echange for heat flux 
    605       !! 
    606       !! ** Method     : gamma assume constant or depends of u* and stability 
    607       !! 
    608       !! ** References : Holland and Jenkins, 1999, JPO, p1787-1800, eq 14 
    609       !!                Jenkins et al., 2010, JPO, p2298-2312 
    610       !!--------------------------------------------------------------------- 
    611       REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pgt   , pgs      !  
    612       REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pqhisf, pqwisf   !  
    613       ! 
    614       INTEGER  :: ji, jj                     ! loop index 
    615       INTEGER  :: ikt                        ! local integer 
    616       REAL(wp) :: zdku, zdkv                 ! U, V shear  
    617       REAL(wp) :: zPr, zSc, zRc              ! Prandtl, Scmidth and Richardson number  
    618       REAL(wp) :: zmob, zmols                ! Monin Obukov length, coriolis factor at T point 
    619       REAL(wp) :: zbuofdep, zhnu             ! Bouyancy length scale, sublayer tickness 
    620       REAL(wp) :: zhmax                      ! limitation of mol 
    621       REAL(wp) :: zetastar                   ! stability parameter 
    622       REAL(wp) :: zgmolet, zgmoles, zgturb   ! contribution of modelecular sublayer and turbulence  
    623       REAL(wp) :: zcoef                      ! temporary coef 
    624       REAL(wp) :: zdep 
    625       REAL(wp) :: zeps = 1.0e-20_wp     
    626       REAL(wp), PARAMETER :: zxsiN = 0.052_wp   ! dimensionless constant 
    627       REAL(wp), PARAMETER :: znu   = 1.95e-6_wp ! kinamatic viscosity of sea water (m2.s-1) 
    628       REAL(wp), DIMENSION(2) :: zts, zab 
    629       REAL(wp), DIMENSION(jpi,jpj) :: zustar   ! U, V at T point and friction velocity 
    630       !!--------------------------------------------------------------------- 
    631       ! 
    632       SELECT CASE ( nn_gammablk ) 
    633       CASE ( 0 ) ! gamma is constant (specified in namelist) 
    634          !! ISOMIP formulation (Hunter et al, 2006) 
    635          pgt(:,:) = rn_gammat0 
    636          pgs(:,:) = rn_gammas0 
    637  
    638       CASE ( 1 ) ! gamma is assume to be proportional to u* 
    639          !! Jenkins et al., 2010, JPO, p2298-2312 
    640          !! Adopted by Asay-Davis et al. (2015) 
    641          !! compute ustar (eq. 24) 
    642 !!gm  NB  use pCdU here so that it will incorporate local boost of Cd0 and log layer case : 
    643 !!         zustar(:,:) = SQRT( rCdU_top(:,:) * SQRT(utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + r_ke0_top) ) 
    644 !! or better :  compute ustar in zdfdrg  and use it here as well as in TKE, GLS and Co 
    645 !! 
    646 !!     ===>>>>    GM  to be done this chrismas 
    647 !!          
    648 !!gm end 
    649          zustar(:,:) = SQRT( r_Cdmin_top * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + r_ke0_top) ) 
    650  
    651          !! Compute gammats 
    652          pgt(:,:) = zustar(:,:) * rn_gammat0 
    653          pgs(:,:) = zustar(:,:) * rn_gammas0 
    654        
    655       CASE ( 2 ) ! gamma depends of stability of boundary layer 
    656          !! Holland and Jenkins, 1999, JPO, p1787-1800, eq 14 
    657          !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO) 
    658          !! compute ustar 
    659          zustar(:,:) = SQRT( r_Cdmin_top * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + r_ke0_top) ) 
    660  
    661          !! compute Pr and Sc number (can be improved) 
    662          zPr =   13.8_wp 
    663          zSc = 2432.0_wp 
    664  
    665          !! compute gamma mole 
    666          zgmolet = 12.5_wp * zPr ** (2.0/3.0) - 6.0_wp 
    667          zgmoles = 12.5_wp * zSc ** (2.0/3.0) - 6.0_wp 
    668  
    669          !! compute gamma 
    670          DO ji = 2, jpi 
    671             DO jj = 2, jpj 
    672                ikt = mikt(ji,jj) 
    673  
    674                IF( zustar(ji,jj) == 0._wp ) THEN           ! only for kt = 1 I think 
    675                   pgt = rn_gammat0 
    676                   pgs = rn_gammas0 
    677                ELSE 
    678                   !! compute Rc number (as done in zdfric.F90) 
    679 !!gm better to do it like in the new zdfric.F90   i.e. avm weighted Ri computation 
    680 !!gm moreover, use Max(rn2,0) to take care of static instabilities.... 
    681                   zcoef = 0.5_wp / e3w_n(ji,jj,ikt+1) 
    682                   !                                            ! shear of horizontal velocity 
    683                   zdku = zcoef * (  un(ji-1,jj  ,ikt  ) + un(ji,jj,ikt  )  & 
    684                      &             -un(ji-1,jj  ,ikt+1) - un(ji,jj,ikt+1)  ) 
    685                   zdkv = zcoef * (  vn(ji  ,jj-1,ikt  ) + vn(ji,jj,ikt  )  & 
    686                      &             -vn(ji  ,jj-1,ikt+1) - vn(ji,jj,ikt+1)  ) 
    687                   !                                            ! richardson number (minimum value set to zero) 
    688                   zRc = rn2(ji,jj,ikt+1) / MAX( zdku*zdku + zdkv*zdkv, zeps ) 
    689  
    690                   !! compute bouyancy  
    691                   zts(jp_tem) = ttbl(ji,jj) 
    692                   zts(jp_sal) = stbl(ji,jj) 
    693                   zdep        = gdepw_n(ji,jj,ikt) 
    694                   ! 
    695                   CALL eos_rab( zts, zdep, zab ) 
    696                   ! 
    697                   !! compute length scale  
    698                   zbuofdep = grav * ( zab(jp_tem) * pqhisf(ji,jj) - zab(jp_sal) * pqwisf(ji,jj) )  !!!!!!!!!!!!!!!!!!!!!!!!!!!! 
    699  
    700                   !! compute Monin Obukov Length 
    701                   ! Maximum boundary layer depth 
    702                   zhmax = gdept_n(ji,jj,mbkt(ji,jj)) - gdepw_n(ji,jj,mikt(ji,jj)) -0.001_wp 
    703                   ! Compute Monin obukhov length scale at the surface and Ekman depth: 
    704                   zmob   = zustar(ji,jj) ** 3 / (vkarmn * (zbuofdep + zeps)) 
    705                   zmols  = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt) 
    706  
    707                   !! compute eta* (stability parameter) 
    708                   zetastar = 1._wp / ( SQRT(1._wp + MAX(zxsiN * zustar(ji,jj) / ( ABS(ff_f(ji,jj)) * zmols * zRc ), 0._wp))) 
    709  
    710                   !! compute the sublayer thickness 
    711                   zhnu = 5 * znu / zustar(ji,jj) 
    712  
    713                   !! compute gamma turb 
    714                   zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / ( ABS(ff_f(ji,jj)) * zhnu )) & 
    715                   &      + 1._wp / ( 2 * zxsiN * zetastar ) - 1._wp / vkarmn 
    716  
    717                   !! compute gammats 
    718                   pgt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet) 
    719                   pgs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles) 
    720                END IF 
    721             END DO 
    722          END DO 
    723          CALL lbc_lnk_multi( 'sbcisf', pgt, 'T', 1., pgs, 'T', 1.) 
    724       END SELECT 
    725       ! 
    726    END SUBROUTINE sbc_isf_gammats 
    727  
    728  
    729    SUBROUTINE sbc_isf_tbl( pvarin, pvarout, cd_ptin ) 
    730       !!---------------------------------------------------------------------- 
    731       !!                  ***  SUBROUTINE sbc_isf_tbl  *** 
    732       !! 
    733       !! ** Purpose : compute mean T/S/U/V in the boundary layer at T- point 
    734       !! 
    735       !!---------------------------------------------------------------------- 
    736       REAL(wp), DIMENSION(:,:,:), INTENT(in   ) :: pvarin 
    737       REAL(wp), DIMENSION(:,:)  , INTENT(  out) :: pvarout 
    738       CHARACTER(len=1),           INTENT(in   ) :: cd_ptin ! point of variable in/out 
    739       ! 
    740       INTEGER ::   ji, jj, jk                ! loop index 
    741       INTEGER ::   ikt, ikb                    ! top and bottom index of the tbl 
    742       REAL(wp) ::   ze3, zhk 
    743       REAL(wp), DIMENSION(jpi,jpj) :: zhisf_tbl ! thickness of the tbl 
    744       !!---------------------------------------------------------------------- 
    745        
    746       ! initialisation 
    747       pvarout(:,:)=0._wp 
    748     
    749       SELECT CASE ( cd_ptin ) 
    750       CASE ( 'U' ) ! compute U in the top boundary layer at T- point  
    751          DO jj = 1,jpj 
    752             DO ji = 1,jpi 
    753                ikt = miku(ji,jj) ; ikb = miku(ji,jj) 
    754                ! thickness of boundary layer at least the top level thickness 
    755                zhisf_tbl(ji,jj) = MAX( rhisf_tbl_0(ji,jj) , e3u_n(ji,jj,ikt) ) 
    756  
    757                ! determine the deepest level influenced by the boundary layer 
    758                DO jk = ikt+1, mbku(ji,jj) 
    759                   IF ( (SUM(e3u_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk 
    760                END DO 
    761                zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3u_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness. 
    762  
    763                ! level fully include in the ice shelf boundary layer 
    764                DO jk = ikt, ikb - 1 
    765                   ze3 = e3u_n(ji,jj,jk) 
    766                   pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3 
    767                END DO 
    768  
    769                ! level partially include in ice shelf boundary layer  
    770                zhk = SUM( e3u_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj) 
    771                pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk) 
    772             END DO 
    773          END DO 
    774          DO jj = 2, jpj 
    775             DO ji = 2, jpi 
    776 !!gm a wet-point only average should be used here !!! 
    777                pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji-1,jj)) 
    778             END DO 
    779          END DO 
    780          CALL lbc_lnk('sbcisf', pvarout,'T',-1.) 
    781        
    782       CASE ( 'V' ) ! compute V in the top boundary layer at T- point  
    783          DO jj = 1,jpj 
    784             DO ji = 1,jpi 
    785                ikt = mikv(ji,jj) ; ikb = mikv(ji,jj) 
    786                ! thickness of boundary layer at least the top level thickness 
    787                zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3v_n(ji,jj,ikt)) 
    788  
    789                ! determine the deepest level influenced by the boundary layer 
    790                DO jk = ikt+1, mbkv(ji,jj) 
    791                   IF ( (SUM(e3v_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk 
    792                END DO 
    793                zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3v_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness. 
    794  
    795                ! level fully include in the ice shelf boundary layer 
    796                DO jk = ikt, ikb - 1 
    797                   ze3 = e3v_n(ji,jj,jk) 
    798                   pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3 
    799                END DO 
    800  
    801                ! level partially include in ice shelf boundary layer  
    802                zhk = SUM( e3v_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj) 
    803                pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk) 
    804             END DO 
    805          END DO 
    806          DO jj = 2, jpj 
    807             DO ji = 2, jpi 
    808 !!gm a wet-point only average should be used here !!! 
    809                pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji,jj-1)) 
    810             END DO 
    811          END DO 
    812          CALL lbc_lnk('sbcisf', pvarout,'T',-1.) 
    813  
    814       CASE ( 'T' ) ! compute T in the top boundary layer at T- point  
    815          DO jj = 1,jpj 
    816             DO ji = 1,jpi 
    817                ikt = misfkt(ji,jj) 
    818                ikb = misfkb(ji,jj) 
    819  
    820                ! level fully include in the ice shelf boundary layer 
    821                DO jk = ikt, ikb - 1 
    822                   ze3 = e3t_n(ji,jj,jk) 
    823                   pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3 
    824                END DO 
    825  
    826                ! level partially include in ice shelf boundary layer  
    827                zhk = SUM( e3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) 
    828                pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk) 
    829             END DO 
    830          END DO 
    831       END SELECT 
    832       ! 
    833       ! mask mean tbl value 
    834       pvarout(:,:) = pvarout(:,:) * ssmask(:,:) 
    835       ! 
    836    END SUBROUTINE sbc_isf_tbl 
    837        
    838  
    839    SUBROUTINE sbc_isf_div( phdivn ) 
    840       !!---------------------------------------------------------------------- 
    841       !!                  ***  SUBROUTINE sbc_isf_div  *** 
    842       !!        
    843       !! ** Purpose :   update the horizontal divergence with the runoff inflow 
    844       !! 
    845       !! ** Method  :    
    846       !!                CAUTION : risf_tsc(:,:,jp_sal) is negative (outflow) increase the  
    847       !!                          divergence and expressed in m/s 
    848       !! 
    849       !! ** Action  :   phdivn   decreased by the runoff inflow 
    850       !!---------------------------------------------------------------------- 
    851       REAL(wp), DIMENSION(:,:,:), INTENT( inout ) ::   phdivn   ! horizontal divergence 
    852       !  
    853       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    854       INTEGER  ::   ikt, ikb  
    855       REAL(wp) ::   zhk 
    856       REAL(wp) ::   zfact     ! local scalar 
    857       !!---------------------------------------------------------------------- 
    858       ! 
    859       zfact   = 0.5_wp 
    860       ! 
    861       IF(.NOT.ln_linssh ) THEN     ! need to re compute level distribution of isf fresh water 
    862          DO jj = 1,jpj 
    863             DO ji = 1,jpi 
    864                ikt = misfkt(ji,jj) 
    865                ikb = misfkt(ji,jj) 
    866                ! thickness of boundary layer at least the top level thickness 
    867                rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3t_n(ji,jj,ikt)) 
    868  
    869                ! determine the deepest level influenced by the boundary layer 
    870                DO jk = ikt, mbkt(ji,jj) 
    871                   IF ( (SUM(e3t_n(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 
    872                END DO 
    873                rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness. 
    874                misfkb(ji,jj) = ikb                                                  ! last wet level of the tbl 
    875                r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj) 
    876  
    877                zhk           = SUM( e3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)  ! proportion of tbl cover by cell from ikt to ikb - 1 
    878                ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / e3t_n(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer 
    879             END DO 
    880          END DO 
    881       END IF  
    882       ! 
    883       !==   ice shelf melting distributed over several levels   ==! 
    884       DO jj = 1,jpj 
    885          DO ji = 1,jpi 
    886                ikt = misfkt(ji,jj) 
    887                ikb = misfkb(ji,jj) 
    888                ! level fully include in the ice shelf boundary layer 
    889                DO jk = ikt, ikb - 1 
    890                   phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + ( fwfisf(ji,jj) + fwfisf_b(ji,jj) ) & 
    891                     &              * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact 
    892                END DO 
    893                ! level partially include in ice shelf boundary layer  
    894                phdivn(ji,jj,ikb) = phdivn(ji,jj,ikb) + ( fwfisf(ji,jj) & 
    895                     &            + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj)  
    896          END DO 
    897       END DO 
    898       ! 
    899    END SUBROUTINE sbc_isf_div 
    900  
     179      ! 
     180      IF ( ln_isfcav_mlt ) THEN 
     181         ! 
     182         ! initialisation  of cav variable 
     183         CALL isf_cav_init() 
     184         ! 
     185         ! read cav variable from restart 
     186         IF ( ln_rstart ) CALL isfrst_read('cav', risf_cav_tsc, fwfisf_cav, risf_cav_tsc_b, fwfisf_cav_b) 
     187         ! 
     188      END IF 
     189      ! 
     190      IF ( ln_isfpar_mlt ) THEN 
     191         ! 
     192         ! initialisation  of par variable 
     193         CALL isf_par_init() 
     194         ! 
     195         ! read par variable from restart 
     196         IF ( ln_rstart ) CALL isfrst_read('par', risf_par_tsc, fwfisf_par, risf_par_tsc_b, fwfisf_par_b) 
     197         ! 
     198      END IF 
     199      ! 
     200  END SUBROUTINE isf_mlt_init 
    901201   !!====================================================================== 
    902 END MODULE sbcisf 
     202END MODULE isfmlt 

  • NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/LDF/ldfslp.F90

    r10425 r11395  
    2121   !!---------------------------------------------------------------------- 
    2222   USE oce            ! ocean dynamics and tracers 
     23   USE isf            ! ice shelf 
    2324   USE dom_oce        ! ocean space and time domain 
    2425!   USE ldfdyn         ! lateral diffusion: eddy viscosity coef. 

  • NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/SBC/sbc_oce.F90

    r10882 r11395  
    123123   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fmmflx            !: freshwater budget: freezing/melting          [Kg/m2/s] 
    124124   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rnf    , rnf_b    !: river runoff        [Kg/m2/s]   
    125    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fwfisf , fwfisf_b !: ice shelf melting   [Kg/m2/s]   
    126125   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fwficb , fwficb_b !: iceberg melting [Kg/m2/s]   
    127126 
    … …  
    174173         &      sfx    (jpi,jpj) , sfx_b(jpi,jpj) , emp_tot(jpi,jpj), fmmflx(jpi,jpj), STAT=ierr(2) ) 
    175174         ! 
    176       ALLOCATE( fwfisf  (jpi,jpj), rnf  (jpi,jpj) , sbc_tsc  (jpi,jpj,jpts) , qsr_hc  (jpi,jpj,jpk) ,  & 
    177          &      fwfisf_b(jpi,jpj), rnf_b(jpi,jpj) , sbc_tsc_b(jpi,jpj,jpts) , qsr_hc_b(jpi,jpj,jpk) ,  & 
     175      ALLOCATE( rnf  (jpi,jpj) , sbc_tsc  (jpi,jpj,jpts) , qsr_hc  (jpi,jpj,jpk) ,  & 
     176         &      rnf_b(jpi,jpj) , sbc_tsc_b(jpi,jpj,jpts) , qsr_hc_b(jpi,jpj,jpk) ,  & 
    178177         &      fwficb  (jpi,jpj), fwficb_b(jpi,jpj), STAT=ierr(3) ) 
    179178         ! 

  • NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/SBC/sbccpl.F90

    r10617 r11395  
    3636   USE eosbn2         !  
    3737   USE sbcrnf  , ONLY : l_rnfcpl 
    38    USE sbcisf  , ONLY : l_isfcpl 
     38   USE isf     , ONLY : l_isfcpl, fwfisf_cav, fwfisf_par 
    3939#if defined key_cice 
    4040   USE ice_domain_size, only: ncat 
    … …  
    478478         IF(lwp) WRITE(numout,*) 
    479479         IF(lwp) WRITE(numout,*) '   iceshelf received from oasis ' 
     480         CALL ctl_stop('STOP','not coded') 
    480481      ENDIF 
    481482      ! 
    … …  
    14051406             rnf(:,:)    = rnf(:,:) + fwficb(:,:)   ! iceberg added to runfofs 
    14061407         ENDIF 
    1407          IF( srcv(jpr_isf)%laction )  fwfisf(:,:) = - frcv(jpr_isf)%z3(:,:,1)  ! fresh water flux from the isf (fwfisf <0 mean melting)   
     1408         IF( srcv(jpr_isf)%laction )  fwfisf_par(:,:) = - frcv(jpr_isf)%z3(:,:,1)  ! fresh water flux from the isf (fwfisf <0 mean melting)   
    14081409         
    14091410         IF( ln_mixcpl ) THEN   ;   emp(:,:) = emp(:,:) * xcplmask(:,:,0) + zemp(:,:) * zmsk(:,:) 
    … …  
    17081709      ENDIF 
    17091710      IF( srcv(jpr_isf)%laction ) THEN   ! iceshelf (fwfisf <0 mean melting) 
    1710         fwfisf(:,:) = - frcv(jpr_isf)%z3(:,:,1)   
     1711        fwfisf_par(:,:) = - frcv(jpr_isf)%z3(:,:,1)   
    17111712      ENDIF 
    17121713 
    … …  
    17471748      ENDIF 
    17481749      IF( srcv(jpr_isf)%laction ) THEN   ! iceshelf (fwfisf <0 mean melting) 
    1749         fwfisf(:,:) = - frcv(jpr_isf)%z3(:,:,1) 
     1750        fwfisf_par(:,:) = - frcv(jpr_isf)%z3(:,:,1) 
    17501751      ENDIF 
    17511752      ! 

  • NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/SBC/sbcfwb.F90

    r10570 r11395  
    2020   USE phycst         ! physical constants 
    2121   USE sbcrnf         ! ocean runoffs 
    22    USE sbcisf         ! ice shelf melting contribution 
     22   USE isf            ! ice shelf melting contribution 
    2323   USE sbcssr         ! Sea-Surface damping terms 
    2424   ! 
    … …  
    104104         ! 
    105105         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN 
    106             y_fwfnow(1) = local_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) - snwice_fmass(:,:) ) ) 
     106            y_fwfnow(1) = local_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) - snwice_fmass(:,:) ) ) 
    107107            CALL mpp_delay_sum( 'sbcfwb', 'fwb', y_fwfnow(:), z_fwfprv(:), kt == nitend - nn_fsbc + 1 ) 
    108108            z_fwfprv(1) = z_fwfprv(1) / area 
    … …  
    159159            ztmsk_neg(:,:) = tmask_i(:,:) - ztmsk_pos(:,:) 
    160160            !                                                  ! fwf global mean (excluding ocean to ice/snow exchanges)  
    161             z_fwf     = glob_sum( 'sbcfwb', e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) - snwice_fmass(:,:) ) ) / area 
     161            z_fwf     = glob_sum( 'sbcfwb', e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) - snwice_fmass(:,:) ) ) / area 
    162162            !             
    163163            IF( z_fwf < 0._wp ) THEN         ! spread out over >0 erp area to increase evaporation 

  • NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/SBC/sbcmod.F90

    r10499 r11395  
    3737#endif 
    3838   USE sbcice_cice    ! surface boundary condition: CICE sea-ice model 
    39    USE sbcisf         ! surface boundary condition: ice-shelf 
    4039   USE sbccpl         ! surface boundary condition: coupled formulation 
    4140   USE cpl_oasis3     ! OASIS routines for coupling 
    … …  
    4342   USE sbcrnf         ! surface boundary condition: runoffs 
    4443   USE sbcapr         ! surface boundary condition: atmo pressure  
    45    USE sbcisf         ! surface boundary condition: ice shelf 
     44   USE isf            ! surface boundary condition: ice shelf 
    4645   USE sbcfwb         ! surface boundary condition: freshwater budget 
    4746   USE icbstp         ! Icebergs 
    … …  
    239238#endif 
    240239      ! 
    241       IF( .NOT.ln_isf ) THEN        !* No ice-shelf in the domain : allocate and set to zero 
    242          IF( sbc_isf_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_isf arrays' ) 
    243          fwfisf  (:,:)   = 0._wp   ;   risf_tsc  (:,:,:) = 0._wp 
    244          fwfisf_b(:,:)   = 0._wp   ;   risf_tsc_b(:,:,:) = 0._wp 
    245       END IF 
    246240      IF( nn_ice == 0 ) THEN        !* No sea-ice in the domain : ice fraction is always zero 
    247241         IF( nn_components /= jp_iam_opa )   fr_i(:,:) = 0._wp    ! except for OPA in SAS-OPA coupled case 
    … …  
    329323      IF( ln_ssr      )   CALL sbc_ssr_init            ! Sea-Surface Restoring initialization 
    330324      ! 
    331       IF( ln_isf      )   CALL sbc_isf_init            ! Compute iceshelves 
    332       ! 
    333325                          CALL sbc_rnf_init            ! Runof initialization 
    334326      ! 
    … …  
    397389            rnf_b    (:,:  ) = rnf    (:,:  ) 
    398390            rnf_tsc_b(:,:,:) = rnf_tsc(:,:,:) 
    399          ENDIF 
    400          IF( ln_isf )  THEN 
    401             fwfisf_b  (:,:  ) = fwfisf  (:,:  )                
    402             risf_tsc_b(:,:,:) = risf_tsc(:,:,:)               
    403391         ENDIF 
    404392        ! 
    … …  
    450438         IF( .NOT. ln_passive_mode ) CALL lbc_lnk( 'sbcmod', emp, 'T', 1. ) ! ensure restartability with icebergs 
    451439      ENDIF 
    452  
    453       IF( ln_isf         )   CALL sbc_isf( kt )                   ! compute iceshelves 
    454440 
    455441      IF( ln_rnf         )   CALL sbc_rnf( kt )                   ! add runoffs to fresh water fluxes 
    … …  
    554540      ! 
    555541      IF(ln_ctl) THEN         ! print mean trends (used for debugging) 
    556          CALL prt_ctl(tab2d_1=fr_i              , clinfo1=' fr_i    - : ', mask1=tmask ) 
    557          CALL prt_ctl(tab2d_1=(emp-rnf + fwfisf), clinfo1=' emp-rnf - : ', mask1=tmask ) 
    558          CALL prt_ctl(tab2d_1=(sfx-rnf + fwfisf), clinfo1=' sfx-rnf - : ', mask1=tmask ) 
     542         CALL prt_ctl(tab2d_1=fr_i             , clinfo1=' fr_i     - : ', mask1=tmask ) 
     543         CALL prt_ctl(tab2d_1=(emp-rnf)        , clinfo1=' emp-rnf - : ', mask1=tmask ) 
     544         CALL prt_ctl(tab2d_1=(sfx-rnf)        , clinfo1=' sfx-rnf - : ', mask1=tmask ) 
    559545         CALL prt_ctl(tab2d_1=qns              , clinfo1=' qns      - : ', mask1=tmask ) 
    560546         CALL prt_ctl(tab2d_1=qsr              , clinfo1=' qsr      - : ', mask1=tmask ) 

  • NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/SBC/sbcrnf.F90

    r10523 r11395  
    1919   USE phycst         ! physical constants 
    2020   USE sbc_oce        ! surface boundary condition variables 
    21    USE sbcisf         ! PM we could remove it I think 
     21   USE isf            ! ice shelf 
    2222   USE eosbn2         ! Equation Of State 
    2323   USE closea         ! closed seas 

  • NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/tranxt.F90

    r10425 r11395  
    2828   USE sbc_oce         ! surface boundary condition: ocean 
    2929   USE sbcrnf          ! river runoffs 
    30    USE sbcisf          ! ice shelf melting 
     30   USE isf             ! ice shelf melting 
    3131   USE zdf_oce         ! ocean vertical mixing 
    3232   USE domvvl          ! variable volume 
    … …  
    312312                  ztc_f  = ztc_n  + atfp * ztc_d 
    313313                  ! 
    314                   IF( jk == mikt(ji,jj) ) THEN           ! first level  
    315                      ze3t_f = ze3t_f - zfact2 * ( (emp_b(ji,jj)    - emp(ji,jj)   )  & 
    316                             &                   + (fwfisf_b(ji,jj) - fwfisf(ji,jj))  ) 
     314                  IF( jk == 1 ) THEN           ! first level  
     315                     ze3t_f = ze3t_f - zfact2 * ( emp_b(ji,jj) - emp(ji,jj)   ) 
    317316                     ztc_f  = ztc_f  - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) ) 
    318317                  ENDIF 
    319318                  IF( ln_rnf_depth ) THEN 
    320319                     ! Rivers are not just at the surface must go down to nk_rnf(ji,jj) 
    321                      IF( mikt(ji,jj) <=jk .and. jk <= nk_rnf(ji,jj)  ) THEN 
     320                     IF( jk <= nk_rnf(ji,jj)  ) THEN 
    322321                        ze3t_f = ze3t_f - zfact2 * ( - (rnf_b(ji,jj) - rnf(ji,jj)   )  ) & 
    323322                    &                            * ( e3t_n(ji,jj,jk) / h_rnf(ji,jj) )  
    324323                     ENDIF 
    325324                  ELSE 
    326                      IF( jk == mikt(ji,jj) ) THEN           ! first level  
     325                     IF( jk == 1 ) THEN           ! first level  
    327326                        ze3t_f = ze3t_f - zfact2 * ( - (rnf_b(ji,jj)    - rnf(ji,jj)   ) )  
    328327                     ENDIF 
    … …  
    341340                  ! ice shelf 
    342341                  IF( ll_isf ) THEN 
    343                      ! level fully include in the Losch_2008 ice shelf boundary layer 
    344                      IF ( jk >= misfkt(ji,jj) .AND. jk < misfkb(ji,jj) )                          & 
    345                         ztc_f  = ztc_f  - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) )  & 
    346                                &                 * e3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj) 
    347                      ! level partially include in Losch_2008 ice shelf boundary layer  
    348                      IF ( jk == misfkb(ji,jj) )                                                   & 
    349                         ztc_f  = ztc_f  - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) )  & 
    350                                &                 * e3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj) * ralpha(ji,jj) 
     342                     IF ( ln_isfcav_mlt ) THEN 
     343                        ! level fully include in the Losch_2008 ice shelf boundary layer 
     344                        IF ( jk >= misfkt_cav(ji,jj) .AND. jk < misfkb_cav(ji,jj) ) THEN 
     345                           ztc_f  = ztc_f  - zfact1 * ( risf_cav_tsc(ji,jj,jn) - risf_cav_tsc_b(ji,jj,jn) ) & 
     346                              &                     * e3t_n(ji,jj,jk) / rhisf_tbl_cav(ji,jj) 
     347                           ze3t_f = ze3t_f - zfact2 * ( fwfisf_cav_b(ji,jj) - fwfisf_cav(ji,jj) )           & 
     348                              &                     * e3t_n(ji,jj,jk) / rhisf_tbl_cav(ji,jj) 
     349                        END IF 
     350                        ! level partially include in Losch_2008 ice shelf boundary layer  
     351                        IF ( jk == misfkb_cav(ji,jj) ) THEN 
     352                           ztc_f  = ztc_f  - zfact1 * ( risf_cav_tsc(ji,jj,jn) - risf_cav_tsc_b(ji,jj,jn) )  & 
     353                                  &                 * e3t_n(ji,jj,jk) / rhisf_tbl_cav(ji,jj) * rfrac_tbl_cav(ji,jj) 
     354                           ze3t_f = ze3t_f - zfact2 * ( fwfisf_cav_b(ji,jj) - fwfisf_cav(ji,jj) )            & 
     355                              &                     * e3t_n(ji,jj,jk) / rhisf_tbl_cav(ji,jj) * rfrac_tbl_cav(ji,jj) 
     356                        END IF 
     357                     END IF 
     358                     IF ( ln_isfpar_mlt ) THEN 
     359                        ! level fully include in the Losch_2008 ice shelf boundary layer 
     360                        IF ( jk >= misfkt_par(ji,jj) .AND. jk < misfkb_par(ji,jj) ) THEN 
     361                           ztc_f  = ztc_f  - zfact1 * ( risf_par_tsc(ji,jj,jn) - risf_par_tsc_b(ji,jj,jn) )  & 
     362                                  &                 * e3t_n(ji,jj,jk) / rhisf_tbl_par(ji,jj) 
     363                           ze3t_f = ze3t_f - zfact2 * ( fwfisf_par_b(ji,jj) - fwfisf_par(ji,jj) )            & 
     364                              &                     * e3t_n(ji,jj,jk) / rhisf_tbl_par(ji,jj) 
     365                        END IF 
     366                        ! level partially include in Losch_2008 ice shelf boundary layer  
     367                        IF ( jk == misfkb_par(ji,jj) ) THEN 
     368                           ztc_f  = ztc_f  - zfact1 * ( risf_par_tsc(ji,jj,jn) - risf_par_tsc_b(ji,jj,jn) )  & 
     369                                  &                 * e3t_n(ji,jj,jk) / rhisf_tbl_par(ji,jj) * rfrac_tbl_par(ji,jj) 
     370                           ze3t_f = ze3t_f - zfact2 * ( fwfisf_par_b(ji,jj) - fwfisf_par(ji,jj) )            & 
     371                              &                     * e3t_n(ji,jj,jk) / rhisf_tbl_par(ji,jj) * rfrac_tbl_par(ji,jj) 
     372                        END IF 
     373                     END IF 
    351374                  END IF 
    352375                  ! 

  • NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/TRA/trasbc.F90

    r10499 r11395  
    2222   USE sbcmod         ! ln_rnf   
    2323   USE sbcrnf         ! River runoff   
    24    USE sbcisf         ! Ice shelf    
    25    USE iscplini       ! Ice sheet coupling 
    2624   USE traqsr         ! solar radiation penetration 
    2725   USE trd_oce        ! trends: ocean variables 
    … …  
    155153      ! 
    156154      !---------------------------------------- 
    157       !       Ice Shelf effects (ISF) 
    158       !     tbl treated as in Losh (2008) JGR 
    159       !---------------------------------------- 
    160       ! 
    161 !!gm BUG ?   Why no differences between non-linear and linear free surface ? 
    162 !!gm         probably taken into account in r1_hisf_tbl : to be verified 
    163       IF( ln_isf ) THEN 
    164          zfact = 0.5_wp 
    165          DO jj = 2, jpj 
    166             DO ji = fs_2, fs_jpim1 
    167                ! 
    168                ikt = misfkt(ji,jj) 
    169                ikb = misfkb(ji,jj) 
    170                ! 
    171                ! level fully include in the ice shelf boundary layer 
    172                ! sign - because fwf sign of evapo (rnf sign of precip) 
    173                DO jk = ikt, ikb - 1 
    174                ! compute trend 
    175                   tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                                & 
    176                      &           + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem) )             & 
    177                      &           * r1_hisf_tbl(ji,jj) 
    178                END DO 
    179     
    180                ! level partially include in ice shelf boundary layer  
    181                ! compute trend 
    182                tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem)                                                 & 
    183                   &              + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem) )             & 
    184                   &              * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) 
    185  
    186             END DO 
    187          END DO 
    188       END IF 
    189       ! 
    190       !---------------------------------------- 
    191155      !        River Runoff effects 
    192156      !---------------------------------------- 
    … …  
    242206#endif 
    243207      ! 
    244       !---------------------------------------- 
    245       !        Ice Sheet coupling imbalance correction to have conservation 
    246       !---------------------------------------- 
    247       ! 
    248       IF( ln_iscpl .AND. ln_hsb) THEN         ! input of heat and salt due to river runoff  
    249          DO jk = 1,jpk 
    250             DO jj = 2, jpj  
    251                DO ji = fs_2, fs_jpim1 
    252                   zdep = 1._wp / e3t_n(ji,jj,jk)  
    253                   tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) - htsc_iscpl(ji,jj,jk,jp_tem) * zdep 
    254                   tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) - htsc_iscpl(ji,jj,jk,jp_sal) * zdep   
    255                END DO   
    256             END DO   
    257          END DO 
    258       ENDIF 
    259  
    260208      IF( l_trdtra )   THEN                      ! save the horizontal diffusive trends for further diagnostics 
    261209         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 

  • NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ZDF/zdfmxl.F90

    r10425 r11395  
    1212   !!---------------------------------------------------------------------- 
    1313   USE oce            ! ocean dynamics and tracers variables 
     14   USE isf            ! ice shelf 
    1415   USE dom_oce        ! ocean space and time domain variables 
    1516   USE trc_oce  , ONLY: l_offline         ! ocean space and time domain variables 

  • NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/nemogcm.F90

    r10588 r11395  
    6060   USE diacfl         ! CFL diagnostics               (dia_cfl_init routine) 
    6161   USE step           ! NEMO time-stepping                 (stp     routine) 
     62   USE isfmlt         ! ice shelf                     (isf_mlt_init routine) 
    6263   USE icbini         ! handle bergs, initialisation 
    6364   USE icbstp         ! handle bergs, calving, themodynamics and transport 
    … …  
    452453      !                                      ! Active tracers 
    453454      IF( ln_traqsr    )   CALL tra_qsr_init      ! penetrative solar radiation qsr 
     455      IF (ln_isf       )   CALL isf_mlt_init      ! ice shelf 
    454456                           CALL tra_bbc_init      ! bottom heat flux 
    455457                           CALL tra_bbl_init      ! advective (and/or diffusive) bottom boundary layer scheme 

  • NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/step.F90

    r11287 r11395  
    113113      IF( ln_apr_dyn )   CALL sbc_apr ( kstp )                   ! atmospheric pressure (NB: call before bdy_dta which needs ssh_ib)  
    114114      IF( ln_bdy     )   CALL bdy_dta ( kstp, time_offset=+1 )   ! update dynamic & tracer data at open boundaries 
     115      IF( ln_isf     )   CALL isf_mlt ( kstp ) 
    115116                         CALL sbc     ( kstp )                   ! Sea Boundary Condition (including sea-ice) 
    116117 
    … …  
    240241                         CALL tra_sbc       ( kstp )  ! surface boundary condition 
    241242      IF( ln_traqsr  )   CALL tra_qsr       ( kstp )  ! penetrative solar radiation qsr 
     243      IF( ln_isf     )   CALL tra_isf       ( kstp )  ! ice shelf heat flux 
    242244      IF( ln_trabbc  )   CALL tra_bbc       ( kstp )  ! bottom heat flux 
    243245      IF( ln_trabbl  )   CALL tra_bbl       ( kstp )  ! advective (and/or diffusive) bottom boundary layer scheme 

  • NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/step_oce.F90

    r10068 r11395  
    2222   USE sbcwave         ! Wave intialisation 
    2323 
     24   USE isfmlt          ! ice shelf boundary condition     (isf_mlt routine) 
     25 
    2426   USE traqsr          ! solar radiation penetration      (tra_qsr routine) 
     27   USE traisf          ! ice shelf                        (tra_isf routine) 
    2528   USE trasbc          ! surface boundary condition       (tra_sbc routine) 
    2629   USE trabbc          ! bottom boundary condition        (tra_bbc routine) 
Note: See TracChangeset for help on using the changeset viewer.