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 12101 for utils/tools_UKMO_MERGE_2019/DOMAINcfg/src – NEMO

Ignore:
Timestamp:
2019-12-06T18:45:39+01:00 (4 years ago)
Author:
mathiot
Message:

merge ENHANCE-03_domcfg and ENHANCE-02_ISF_nemo

Location:
utils/tools_UKMO_MERGE_2019/DOMAINcfg/src
Files:
41 deleted
3 edited
47 copied

Legend:

Unmodified
Added
Removed

  • utils/tools_UKMO_MERGE_2019/DOMAINcfg/src/agrif_user.F90

    r12100 r12101  
    1212      !!---------------------------------------------------------------------- 
    1313   USE Agrif_Util 
    14    USE oce  
    1514   USE dom_oce 
    1615   USE nemogcm 

  • utils/tools_UKMO_MERGE_2019/DOMAINcfg/src/dom_oce.F90

    r12100 r12101  
    3636   REAL(wp), PUBLIC ::   rn_e3zps_rat    !: minimum thickness ration for partial steps 
    3737   INTEGER , PUBLIC ::   nn_msh          !: = 1 create a mesh-mask file 
    38    INTEGER , PUBLIC ::   nn_closea       !: =0 suppress closed sea/lake from the ORCA domain or not (=1) 
     38 
     39   INTEGER , PUBLIC ::   ntopo           !: = 0/1 ,compute/read the bathymetry file 
     40   INTEGER, PUBLIC  ::   nperio          !: type of lateral boundary condition 
     41   REAL(wp), PUBLIC ::   e3zps_min       !: miminum thickness for partial steps (meters) 
     42   REAL(wp), PUBLIC ::   e3zps_rat       !: minimum thickness ration for partial steps 
    3943 
    4044   INTEGER, PUBLIC :: nn_interp 
    … …  
    5054   LOGICAL, PUBLIC ::   lzoom_n    =  .FALSE.   !: North zoom type flag 
    5155 
     56   LOGICAL, PUBLIC ::   ln_domclo  =  .FALSE. 
    5257 
    5358   INTEGER       ::   jphgr_msh          !: type of horizontal mesh 
    … …  
    202207   LOGICAL, PUBLIC ::   ln_sco       !: s-coordinate or hybrid z-s coordinate 
    203208   LOGICAL, PUBLIC ::   ln_isfcav    !: presence of ISF  
    204    !                                                        !  ref.   ! before  !   now   ! after  ! 
    205    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3t_0 ,   e3t_b ,   e3t_n ,  e3t_a   !: t- vert. scale factor [m] 
    206    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3u_0 ,   e3u_b ,   e3u_n ,  e3u_a   !: u- vert. scale factor [m] 
    207    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3v_0 ,   e3v_b ,   e3v_n ,  e3v_a   !: v- vert. scale factor [m] 
    208    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3f_0           ,   e3f_n            !: f- vert. scale factor [m] 
    209    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3w_0 ,   e3w_b ,   e3w_n            !: w- vert. scale factor [m] 
    210    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::    e3uw_0 ,  e3uw_b ,  e3uw_n            !: uw-vert. scale factor [m] 
    211    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::    e3vw_0 ,  e3vw_b ,  e3vw_n            !: vw-vert. scale factor [m] 
    212  
    213    !                                                        !  ref.   ! before  !   now   ! 
    214    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdept_0 , gdept_b , gdept_n   !: t- depth              [m] 
    215    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdepw_0 , gdepw_b , gdepw_n   !: w- depth              [m] 
    216    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gde3w_0           , gde3w_n   !: w- depth (sum of e3w) [m] 
    217     
    218    !                                                      !  ref. ! before  !   now   !  after  ! 
    219    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ht_0            ,    ht_n             !: t-depth              [m] 
    220    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hu_0  ,    hu_b ,    hu_n ,    hu_a   !: u-depth              [m] 
    221    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hv_0  ,    hv_b ,    hv_n ,    hv_a   !: v-depth              [m] 
    222    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::           r1_hu_b , r1_hu_n , r1_hu_a   !: inverse of u-depth [1/m] 
    223    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::           r1_hv_b , r1_hv_n , r1_hv_a   !: inverse of v-depth [1/m] 
    224  
     209   ! 
     210   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  e3t_0, e3u_0 , e3v_0 , e3f_0 !: t-,u-,v-,f-vert. scale factor [m] 
     211   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  e3w_0, e3uw_0, e3vw_0        !: w-,uw-,vw-vert. scale factor [m] 
     212   ! 
     213   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdept_0 !: t- depth              [m] 
     214   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdepw_0 !: w- depth              [m] 
     215   ! 
    225216   INTEGER, PUBLIC ::   nla10              !: deepest    W level Above  ~10m (nlb10 - 1) 
    226217   INTEGER, PUBLIC ::   nlb10              !: shallowest W level Bellow ~10m (nla10 + 1)  
    227  
     218   ! 
    228219   !! 1D reference  vertical coordinate 
    229220   !! =-----------------====------ 
    … …  
    251242   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ssmask, ssumask, ssvmask             !: surface mask at T-,U-, V- and F-pts 
    252243   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, fmask   !: land/ocean mask at T-, U-, V- and F-pts 
    253    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: wmask, wumask, wvmask        !: land/ocean mask at WT-, WU- and WV-pts 
     244   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: wmask                        !: land/ocean mask at W- pts                
     245 
     246   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: msk_opnsea  , msk_csundef                !: open ocean mask, closed sea mask (all of them) 
     247   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: msk_csglo   , msk_csrnf   , msk_csemp    !: closed sea masks 
     248   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: msk_csgrpglo, msk_csgrprnf, msk_csgrpemp !: closed sea masks 
    254249 
    255250   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   tpol, fpol          !: north fold mask (jperio= 3 or 4) 
    … …  
    310305   INTEGER FUNCTION dom_oce_alloc() 
    311306      !!---------------------------------------------------------------------- 
    312       INTEGER, DIMENSION(12) :: ierr 
     307      INTEGER, DIMENSION(11) :: ierr 
    313308      !!---------------------------------------------------------------------- 
    314309      ierr(:) = 0 
    … …  
    331326         &      ff_f (jpi,jpj) ,    ff_t (jpi,jpj)                                     , STAT=ierr(3) ) 
    332327         ! 
    333       ALLOCATE( gdept_0(jpi,jpj,jpk) , gdepw_0(jpi,jpj,jpk) , gde3w_0(jpi,jpj,jpk) ,      & 
    334          &      gdept_b(jpi,jpj,jpk) , gdepw_b(jpi,jpj,jpk) ,                             & 
    335          &      gdept_n(jpi,jpj,jpk) , gdepw_n(jpi,jpj,jpk) , gde3w_n(jpi,jpj,jpk) , STAT=ierr(4) ) 
    336          ! 
    337       ALLOCATE( e3t_0(jpi,jpj,jpk) , e3u_0(jpi,jpj,jpk) , e3v_0(jpi,jpj,jpk) , e3f_0(jpi,jpj,jpk) , e3w_0(jpi,jpj,jpk) ,   & 
    338          &      e3t_b(jpi,jpj,jpk) , e3u_b(jpi,jpj,jpk) , e3v_b(jpi,jpj,jpk) ,                      e3w_b(jpi,jpj,jpk) ,   &  
    339          &      e3t_n(jpi,jpj,jpk) , e3u_n(jpi,jpj,jpk) , e3v_n(jpi,jpj,jpk) , e3f_n(jpi,jpj,jpk) , e3w_n(jpi,jpj,jpk) ,   &  
    340          &      e3t_a(jpi,jpj,jpk) , e3u_a(jpi,jpj,jpk) , e3v_a(jpi,jpj,jpk) ,                                             & 
    341          !                                                          ! 
    342          &      e3uw_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk) ,         & 
    343          &      e3uw_b(jpi,jpj,jpk) , e3vw_b(jpi,jpj,jpk) ,         &                
    344          &      e3uw_n(jpi,jpj,jpk) , e3vw_n(jpi,jpj,jpk) ,     STAT=ierr(5) )                        
    345          ! 
    346       ALLOCATE( ht_0(jpi,jpj) , hu_0(jpi,jpj) , hv_0(jpi,jpj) ,                                           & 
    347          &                      hu_b(jpi,jpj) , hv_b(jpi,jpj) , r1_hu_b(jpi,jpj) , r1_hv_b(jpi,jpj) ,     & 
    348          &      ht_n(jpi,jpj) , hu_n(jpi,jpj) , hv_n(jpi,jpj) , r1_hu_n(jpi,jpj) , r1_hv_n(jpi,jpj) ,     & 
    349          &                      hu_a(jpi,jpj) , hv_a(jpi,jpj) , r1_hu_a(jpi,jpj) , r1_hv_a(jpi,jpj) , STAT=ierr(6)  ) 
    350          ! 
    351          ! 
    352       ALLOCATE( gdept_1d(jpk) , e3tp (jpi,jpj), e3wp(jpi,jpj) ,gdepw_1d(jpk) , e3t_1d(jpk) , e3w_1d(jpk) , STAT=ierr(7) ) 
     328      ALLOCATE( gdept_0(jpi,jpj,jpk) , gdepw_0(jpi,jpj,jpk) , STAT=ierr(4) ) 
     329         ! 
     330      ALLOCATE( e3t_0 (jpi,jpj,jpk) , e3u_0 (jpi,jpj,jpk) , e3v_0(jpi,jpj,jpk) , e3f_0(jpi,jpj,jpk) , e3w_0(jpi,jpj,jpk) ,   & 
     331         &      e3uw_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk) , STAT=ierr(5) )                        
     332         ! 
     333      ALLOCATE( gdept_1d(jpk) , e3tp (jpi,jpj), e3wp(jpi,jpj) ,gdepw_1d(jpk) , e3t_1d(jpk) , e3w_1d(jpk) , STAT=ierr(6) ) 
    353334         ! 
    354335      ALLOCATE( bathy(jpi,jpj),mbathy(jpi,jpj), tmask_i(jpi,jpj) , tmask_h(jpi,jpj) ,                        &  
    355336         &      ssmask (jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) ,     & 
    356          &      mbkt   (jpi,jpj) , mbku   (jpi,jpj) , mbkv   (jpi,jpj) , STAT=ierr(9) ) 
     337         &      mbkt   (jpi,jpj) , mbku   (jpi,jpj) , mbkv   (jpi,jpj) , STAT=ierr(7) ) 
    357338         ! 
    358339      ALLOCATE( misfdep(jpi,jpj) , mikt(jpi,jpj) , miku(jpi,jpj) ,     & 
    359          &      risfdep(jpi,jpj) , mikv(jpi,jpj) , mikf(jpi,jpj) , STAT=ierr(10) ) 
     340         &      risfdep(jpi,jpj) , mikv(jpi,jpj) , mikf(jpi,jpj) , STAT=ierr(8) ) 
    360341         ! 
    361342      ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk) ,     &  
    362          &      vmask(jpi,jpj,jpk) , fmask(jpi,jpj,jpk) , STAT=ierr(11) ) 
    363          ! 
    364       ALLOCATE( wmask(jpi,jpj,jpk) , wumask(jpi,jpj,jpk), wvmask(jpi,jpj,jpk) , STAT=ierr(12) ) 
    365  
     343         &      vmask(jpi,jpj,jpk) , fmask(jpi,jpj,jpk) , wmask(jpi,jpj,jpk) , STAT=ierr(9) ) 
     344         ! 
    366345      ALLOCATE( hbatv (jpi,jpj) , hbatf (jpi,jpj) ,     & 
    367346         &      hbatt (jpi,jpj) , hbatu (jpi,jpj) ,     & 
    368347         &      scosrf(jpi,jpj) , scobot(jpi,jpj) ,     & 
    369348         &      hifv  (jpi,jpj) , hiff  (jpi,jpj) ,     & 
    370          &      hift  (jpi,jpj) , hifu  (jpi,jpj) , STAT=ierr(8) ) 
     349         &      hift  (jpi,jpj) , hifu  (jpi,jpj) , STAT=ierr(10) ) 
     350 
     351      ALLOCATE( msk_opnsea  (jpi,jpj), msk_csundef (jpi,jpj),                        & 
     352         &      msk_csglo   (jpi,jpj), msk_csrnf   (jpi,jpj), msk_csemp   (jpi,jpj), & 
     353         &      msk_csgrpglo(jpi,jpj), msk_csgrprnf(jpi,jpj), msk_csgrpemp(jpi,jpj), STAT=ierr(11) ) 
    371354      ! 
    372355      dom_oce_alloc = MAXVAL(ierr) 

  • utils/tools_UKMO_MERGE_2019/DOMAINcfg/src/domain.F90

    r12100 r12101  
    2121   !!   dom_ctl        : control print for the ocean domain 
    2222   !!---------------------------------------------------------------------- 
    23    USE oce             ! ocean variables 
    2423   USE dom_oce         ! domain: ocean 
    2524   USE phycst          ! physical constants 
    26  !  USE closea          ! closed seas 
    2725   USE domhgr          ! domain: set the horizontal mesh 
    2826   USE domzgr          ! domain: set the vertical mesh 
    29  !  USE domstp          ! domain: set the time-step 
    3027   USE dommsk          ! domain: set the mask system 
    31    USE domwri          ! domain: write the meshmask file 
    32    USE domvvl          ! variable volume 
     28   USE domclo          ! domain: set closed sea mask 
    3329   ! 
     30   USE lib_mpp         ! 
    3431   USE in_out_manager  ! I/O manager 
    3532   USE iom             !  
    36    USE wrk_nemo        ! Memory Allocation 
    37    USE lib_mpp         ! distributed memory computing library 
    38    USE lbclnk          ! ocean lateral boundary condition (or mpp link) 
    39    USE timing          ! Timing 
    4033 
    4134   IMPLICIT NONE 
    … …  
    6962      !!              - 1D configuration, move Coriolis, u and v at T-point 
    7063      !!---------------------------------------------------------------------- 
    71       INTEGER ::   jk          ! dummy loop indices 
    72       INTEGER ::   iconf = 0   ! local integers 
    73       REAL(wp), POINTER, DIMENSION(:,:) ::   z1_hu_0, z1_hv_0 
    74       !!---------------------------------------------------------------------- 
    75       ! 
    76      ! IF( nn_timing == 1 )   CALL timing_start('dom_init') 
    7764      ! 
    7865      IF(lwp) THEN 
    … …  
    8471      !                       !==  Reference coordinate system  ==! 
    8572      ! 
    86                      CALL dom_nam               ! read namelist ( namrun, namdom ) 
    87                   !   CALL dom_clo               ! Closed seas and lake 
    88           
    89                      CALL dom_hgr               ! Horizontal mesh 
    90                      CALL dom_zgr               ! Vertical mesh and bathymetry 
    91                      CALL dom_msk               ! Masks 
    92       ! 
    93       ht_0(:,:) = e3t_0(:,:,1) * tmask(:,:,1)   ! Reference ocean thickness 
    94       hu_0(:,:) = e3u_0(:,:,1) * umask(:,:,1) 
    95       hv_0(:,:) = e3v_0(:,:,1) * vmask(:,:,1) 
    96       DO jk = 2, jpk 
    97          ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk) 
    98          hu_0(:,:) = hu_0(:,:) + e3u_0(:,:,jk) * umask(:,:,jk) 
    99          hv_0(:,:) = hv_0(:,:) + e3v_0(:,:,jk) * vmask(:,:,jk) 
    100       END DO 
    101       ! 
    102       !              !==  time varying part of coordinate system  ==! 
    103       ! 
    104       IF( ln_linssh ) THEN          ! Fix in time : set to the reference one for all 
    105          !       before        !          now          !       after         ! 
    106          ;  gdept_b = gdept_0  ;   gdept_n = gdept_0   !        ---          ! depth of grid-points 
    107          ;  gdepw_b = gdepw_0  ;   gdepw_n = gdepw_0   !        ---          ! 
    108          ;                     ;   gde3w_n = gde3w_0   !        ---          ! 
    109          !                                                                   
    110          ;    e3t_b =   e3t_0  ;     e3t_n =   e3t_0   ;   e3t_a =  e3t_0    ! scale factors 
    111          ;    e3u_b =   e3u_0  ;     e3u_n =   e3u_0   ;   e3u_a =  e3u_0    ! 
    112          ;    e3v_b =   e3v_0  ;     e3v_n =   e3v_0   ;   e3v_a =  e3v_0    ! 
    113          ;                     ;     e3f_n =   e3f_0   !        ---          ! 
    114          ;    e3w_b =   e3w_0  ;     e3w_n =   e3w_0   !        ---          ! 
    115          ;   e3uw_b =  e3uw_0  ;    e3uw_n =  e3uw_0   !        ---          ! 
    116          ;   e3vw_b =  e3vw_0  ;    e3vw_n =  e3vw_0   !        ---          ! 
    117          ! 
    118          CALL wrk_alloc( jpi,jpj,   z1_hu_0, z1_hv_0 ) 
    119          ! 
    120          z1_hu_0(:,:) = ssumask(:,:) / ( hu_0(:,:) + 1._wp - ssumask(:,:) )     ! _i mask due to ISF 
    121          z1_hv_0(:,:) = ssvmask(:,:) / ( hv_0(:,:) + 1._wp - ssvmask(:,:) ) 
    122          ! 
    123          !        before       !          now          !       after         ! 
    124          ;                     ;      ht_n =    ht_0   !                     ! water column thickness 
    125          ;     hu_b =    hu_0  ;      hu_n =    hu_0   ;    hu_a =    hu_0   !  
    126          ;     hv_b =    hv_0  ;      hv_n =    hv_0   ;    hv_a =    hv_0   ! 
    127          ;  r1_hu_b = z1_hu_0  ;   r1_hu_n = z1_hu_0   ; r1_hu_a = z1_hu_0   ! inverse of water column thickness 
    128          ;  r1_hv_b = z1_hv_0  ;   r1_hv_n = z1_hv_0   ; r1_hv_a = z1_hv_0   ! 
    129          ! 
    130          CALL wrk_dealloc( jpi,jpj,   z1_hu_0, z1_hv_0 ) 
    131          ! 
    132       ELSE                         ! time varying : initialize before/now/after variables 
    133          ! 
    134          CALL dom_vvl_init  
    135          ! 
    136       ENDIF 
    137       ! 
    138       CALL cfg_write         ! create the configuration file 
    139       ! 
    140     !  IF( nn_timing == 1 )   CALL timing_stop('dom_init') 
     73      CALL dom_nam                  ! read namelist ( namrun, namdom ) 
     74      ! 
     75      CALL dom_hgr                  ! Horizontal mesh 
     76      ! 
     77      CALL dom_zgr                  ! Vertical mesh and bathymetry 
     78      ! 
     79      CALL dom_msk                  ! compute mask (needed by write_cfg) 
     80      ! 
     81      IF ( ln_domclo ) CALL dom_clo ! Closed seas and lake 
     82      ! 
     83      CALL dom_ctl                  ! print extrema of masked scale factors 
     84      !  
     85      CALL cfg_write                ! create the configuration file 
    14186      ! 
    14287   END SUBROUTINE dom_init 
    143  
    14488 
    14589   SUBROUTINE dom_nam 
    … …  
    160104         &             ln_cfmeta, ln_iscpl 
    161105      NAMELIST/namdom/ nn_bathy, cn_topo, cn_bath, cn_lon, cn_lat, nn_interp,                        & 
    162          &             rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin, rn_isfhmin,           & 
    163          &             rn_atfp , rn_rdt   , nn_closea   , ln_crs      , jphgr_msh ,                  & 
     106         &             rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin,                       & 
     107         &             rn_atfp , rn_rdt   , ln_crs      , jphgr_msh ,                                & 
    164108         &             ppglam0, ppgphi0, ppe1_deg, ppe2_deg, ppe1_m, ppe2_m,                         & 
    165109         &             ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh,                  & 
    … …  
    185129         WRITE(numout,*) '~~~~~~~ ' 
    186130         WRITE(numout,*) '   Namelist namrun' 
    187          WRITE(numout,*) '      job number                      nn_no      = ', nn_no 
    188131         WRITE(numout,*) '      experiment name for output      cn_exp     = ', cn_exp 
    189          WRITE(numout,*) '      file prefix restart input       cn_ocerst_in= ', cn_ocerst_in 
    190          WRITE(numout,*) '      restart input directory         cn_ocerst_indir= ', cn_ocerst_indir 
    191          WRITE(numout,*) '      file prefix restart output      cn_ocerst_out= ', cn_ocerst_out 
    192          WRITE(numout,*) '      restart output directory        cn_ocerst_outdir= ', cn_ocerst_outdir 
    193          WRITE(numout,*) '      restart logical                 ln_rstart  = ', ln_rstart 
    194          WRITE(numout,*) '      start with forward time step    nn_euler   = ', nn_euler 
    195          WRITE(numout,*) '      control of time step            nn_rstctl  = ', nn_rstctl 
    196          WRITE(numout,*) '      number of the first time step   nn_it000   = ', nn_it000 
    197          WRITE(numout,*) '      number of the last time step    nn_itend   = ', nn_itend 
    198          WRITE(numout,*) '      initial calendar date aammjj    nn_date0   = ', nn_date0 
    199          WRITE(numout,*) '      initial time of day in hhmm     nn_time0   = ', nn_time0 
    200          WRITE(numout,*) '      leap year calendar (0/1)        nn_leapy   = ', nn_leapy 
    201          WRITE(numout,*) '      initial state output            nn_istate  = ', nn_istate 
    202          IF( ln_rst_list ) THEN 
    203             WRITE(numout,*) '      list of restart dump times      nn_stocklist   =', nn_stocklist 
    204          ELSE 
    205             WRITE(numout,*) '      frequency of restart file       nn_stock   = ', nn_stock 
    206          ENDIF 
    207          WRITE(numout,*) '      frequency of output file        nn_write   = ', nn_write 
    208132         WRITE(numout,*) '      mask land points                ln_mskland = ', ln_mskland 
    209133         WRITE(numout,*) '      additional CF standard metadata ln_cfmeta  = ', ln_cfmeta 
    210134         WRITE(numout,*) '      overwrite an existing file      ln_clobber = ', ln_clobber 
    211135         WRITE(numout,*) '      NetCDF chunksize (bytes)        nn_chunksz = ', nn_chunksz 
    212          WRITE(numout,*) '      IS coupling at the restart step ln_iscpl   = ', ln_iscpl 
    213136      ENDIF 
    214137 
    … …  
    280203         WRITE(numout,*) '      min depth of the ocean    (>0) or    rn_hmin   = ', rn_hmin 
    281204         WRITE(numout,*) '      min number of ocean level (<0)       ' 
    282          WRITE(numout,*) '      treshold to open the isf cavity   rn_isfhmin   = ', rn_isfhmin, ' (m)' 
    283205         WRITE(numout,*) '      minimum thickness of partial      rn_e3zps_min = ', rn_e3zps_min, ' (m)' 
    284206         WRITE(numout,*) '         step level                     rn_e3zps_rat = ', rn_e3zps_rat 
    … …  
    288210         WRITE(numout,*) '           = 2   mesh and mask             ' 
    289211         WRITE(numout,*) '           = 3   mesh_hgr, msh_zgr and mask' 
    290          WRITE(numout,*) '      ocean time step                       rn_rdt    = ', rn_rdt 
    291          WRITE(numout,*) '      asselin time filter parameter         rn_atfp   = ', rn_atfp 
    292          WRITE(numout,*) '      suppression of closed seas (=0)       nn_closea = ', nn_closea 
    293          WRITE(numout,*) '      online coarsening of dynamical fields ln_crs    = ', ln_crs 
    294212         WRITE(numout,*) '      type of horizontal mesh jphgr_msh           = ', jphgr_msh 
    295213         WRITE(numout,*) '      longitude of first raw and column T-point ppglam0 = ', ppglam0 
    … …  
    337255      !!---------------------------------------------------------------------- 
    338256      ! 
    339 #undef CHECK_DOM 
    340 #ifdef CHECK_DOM 
    341257      IF(lk_mpp) THEN 
    342          CALL mpp_minloc( e1t(:,:), tmask_i(:,:), ze1min, iimi1,ijmi1 ) 
    343          CALL mpp_minloc( e2t(:,:), tmask_i(:,:), ze2min, iimi2,ijmi2 ) 
    344          CALL mpp_maxloc( e1t(:,:), tmask_i(:,:), ze1max, iima1,ijma1 ) 
    345          CALL mpp_maxloc( e2t(:,:), tmask_i(:,:), ze2max, iima2,ijma2 ) 
     258         CALL mpp_minloc( 'dom_ctl', e1t(:,:), tmask_i(:,:), ze1min, iloc ) 
     259         iimi1 = iloc(1) ; ijmi1 = iloc(2) 
     260         CALL mpp_minloc( 'dom_ctl', e2t(:,:), tmask_i(:,:), ze2min, iloc ) 
     261         iimi2 = iloc(1) ; ijmi2 = iloc(2) 
     262         CALL mpp_maxloc( 'dom_ctl', e1t(:,:), tmask_i(:,:), ze1max, iloc ) 
     263         iima1 = iloc(1) ; ijma1 = iloc(2) 
     264         CALL mpp_maxloc( 'dom_ctl', e2t(:,:), tmask_i(:,:), ze2max, iloc ) 
     265         iima2 = iloc(1) ; ijma2 = iloc(2) 
    346266      ELSE 
    347267         ze1min = MINVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp )     
    … …  
    372292         WRITE(numout,"(14x,'e2t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2min, iimi2, ijmi2 
    373293      ENDIF 
    374 #endif 
    375294      ! 
    376295   END SUBROUTINE dom_ctl 
    … …  
    490409      CALL iom_rstput( 0, 0, inum, 'bottom_level' , REAL( mbkt, wp )*ssmask , ktype = jp_i4 )   ! nb of ocean T-points 
    491410      CALL iom_rstput( 0, 0, inum, 'top_level'    , REAL( mikt, wp )*ssmask , ktype = jp_i4 )   ! nb of ocean T-points (ISF) 
    492       DO jj = 1,jpj 
    493          DO ji = 1,jpi 
    494             z2d (ji,jj) = SUM ( e3t_0(ji,jj, 1:mbkt(ji,jj) ) ) * ssmask(ji,jj)  
    495          END DO 
    496       END DO 
    497       CALL iom_rstput( 0, 0, inum, 'bathy_metry'   , z2d , ktype = jp_r4 ) 
    498  
    499       ! 
    500       IF( ln_sco ) THEN             ! s-coordinate: store grid stiffness ratio  (Not required anyway) 
    501          CALL dom_stiff( z2d ) 
    502          CALL iom_rstput( 0, 0, inum, 'stiffness', z2d )        !    ! Max. grid stiffness ratio 
    503       ENDIF 
     411      CALL iom_rstput( 0, 0, inum, 'isf_draft'    , risfdep , ktype = jp_r8 ) 
     412      CALL iom_rstput( 0, 0, inum, 'bathy_metry'  , bathy   , ktype = jp_r8 ) 
     413      ! 
     414      !                              !== closed sea ==! 
     415      IF (ln_domclo) THEN 
     416         ! mask for the open sea 
     417         CALL iom_rstput( 0, 0, inum, 'mask_opensea' , msk_opnsea  , ktype = jp_i4 ) 
     418         ! mask for all the under closed sea 
     419         CALL iom_rstput( 0, 0, inum, 'mask_csundef' , msk_csundef , ktype = jp_i4 ) 
     420         ! mask for global, local net precip, local net precip and evaporation correction 
     421         CALL iom_rstput( 0, 0, inum, 'mask_csglo'   , msk_csglo   , ktype = jp_i4 ) 
     422         CALL iom_rstput( 0, 0, inum, 'mask_csemp'   , msk_csemp   , ktype = jp_i4 ) 
     423         CALL iom_rstput( 0, 0, inum, 'mask_csrnf'   , msk_csrnf   , ktype = jp_i4 ) 
     424         ! mask for the various river mouth (in case multiple lake in the same outlet) 
     425         CALL iom_rstput( 0, 0, inum, 'mask_csgrpglo', msk_csgrpglo, ktype = jp_i4 ) 
     426         CALL iom_rstput( 0, 0, inum, 'mask_csgrpemp', msk_csgrpemp, ktype = jp_i4 ) 
     427         CALL iom_rstput( 0, 0, inum, 'mask_csgrprnf', msk_csgrprnf, ktype = jp_i4 ) 
     428      END IF 
    504429      ! 
    505430      !                                ! ============================ 

  • utils/tools_UKMO_MERGE_2019/DOMAINcfg/src/dombat.F90

    r12100 r12101  
    11MODULE dombat 
    22 
    3    USE oce               ! ocean variables 
    43   USE dom_oce           ! ocean domain 
    54!   USE closea            ! closed seas 
    … …  
    98   USE lbclnk            ! ocean lateral boundary conditions (or mpp link) 
    109   USE lib_mpp           ! distributed memory computing library 
    11    USE wrk_nemo          ! Memory allocation 
    12    USE timing            ! Timing 
    1310   USE agrif_modutil 
    1411   USE bilinear_interp 

  • utils/tools_UKMO_MERGE_2019/DOMAINcfg/src/domcfg.f90

    r9598 r12101  
    1515   USE in_out_manager  ! I/O manager 
    1616   USE lib_mpp         ! distributed memory computing library 
    17    USE timing          ! Timing 
    1817 
    1918   IMPLICIT NONE 
    … …  
    3635      !! 
    3736      !!---------------------------------------------------------------------- 
    38       ! 
    39       IF( nn_timing == 1 )  CALL timing_start('dom_cfg') 
    4037      ! 
    4138      IF(lwp) THEN                   ! Control print 
    … …  
    6057      CALL dom_glo                   ! global domain versus zoom and/or local domain 
    6158      ! 
    62       IF( nn_timing == 1 )  CALL timing_stop('dom_cfg') 
    63       ! 
    6459   END SUBROUTINE dom_cfg 
    65  
    6660 
    6761   SUBROUTINE dom_glo 
    … …  
    6963      !!                     ***  ROUTINE dom_glo  *** 
    7064      !! 
    71       !! ** Purpose :   initialization for global domain, zoom and local domain 
     65      !! ** Purpose :   initialization of global domain <--> local domain indices 
    7266      !! 
    7367      !! ** Method  :    
    7468      !! 
    75       !! ** Action  : - mig  , mjg :  
    76       !!              - mi0  , mi1   : 
    77       !!              - mj0, , mj1   : 
     69      !! ** Action  : - mig , mjg : local  domain indices ==> global domain indices 
     70      !!              - mi0 , mi1 : global domain indices ==> local  domain indices 
     71      !!              - mj0,, mj1   (global point not in the local domain ==> mi0>mi1 and/or mj0>mj1) 
    7872      !!---------------------------------------------------------------------- 
    7973      INTEGER ::   ji, jj   ! dummy loop argument 
    8074      !!---------------------------------------------------------------------- 
    81       !                              ! recalculate jpizoom/jpjzoom given lat/lon 
    8275      ! 
    83       !                        ! ============== ! 
    84       !                        !  Local domain  !  
    85       !                        ! ============== ! 
    86       DO ji = 1, jpi                 ! local domain indices ==> data domain indices 
    87         mig(ji) = ji + jpizoom - 1 + nimpp - 1 
     76      DO ji = 1, jpi                 ! local domain indices ==> global domain indices 
     77        mig(ji) = ji + nimpp - 1 
    8878      END DO 
    8979      DO jj = 1, jpj 
    90         mjg(jj) = jj + jpjzoom - 1 + njmpp - 1 
     80        mjg(jj) = jj + njmpp - 1 
    9181      END DO 
    92       ! 
    93       !                              ! data domain indices ==> local domain indices 
     82      !                              ! global domain indices ==> local domain indices 
    9483      !                                   ! (return (m.0,m.1)=(1,0) if data domain gridpoint is to the west/south of the  
    95       !                                   !local domain, or (m.0,m.1)=(jp.+1,jp.) to the east/north of local domain.  
    96       DO ji = 1, jpidta 
    97         mi0(ji) = MAX( 1, MIN( ji - jpizoom + 1 - nimpp + 1, jpi+1 ) ) 
    98         mi1(ji) = MAX( 0, MIN( ji - jpizoom + 1 - nimpp + 1, jpi   ) ) 
     84      !                                   ! local domain, or (m.0,m.1)=(jp.+1,jp.) to the east/north of local domain.  
     85      DO ji = 1, jpiglo 
     86        mi0(ji) = MAX( 1 , MIN( ji - nimpp + 1, jpi+1 ) ) 
     87        mi1(ji) = MAX( 0 , MIN( ji - nimpp + 1, jpi   ) ) 
    9988      END DO 
    100       DO jj = 1, jpjdta 
    101         mj0(jj) = MAX( 1, MIN( jj - jpjzoom + 1 - njmpp + 1, jpj+1 ) ) 
    102         mj1(jj) = MAX( 0, MIN( jj - jpjzoom + 1 - njmpp + 1, jpj   ) ) 
     89      DO jj = 1, jpjglo 
     90        mj0(jj) = MAX( 1 , MIN( jj - njmpp + 1, jpj+1 ) ) 
     91        mj1(jj) = MAX( 0 , MIN( jj - njmpp + 1, jpj   ) ) 
    10392      END DO 
    10493      IF(lwp) THEN                   ! control print 
    10594         WRITE(numout,*) 
    106          WRITE(numout,*) 'dom_glo : domain: data / local ' 
     95         WRITE(numout,*) 'dom_glo : domain: global <<==>> local ' 
    10796         WRITE(numout,*) '~~~~~~~ ' 
    108          WRITE(numout,*) '          data input domain    : jpidta = ', jpidta,   & 
    109             &                                            ' jpjdta = ', jpjdta, ' jpkdta = ', jpkdta 
    110          WRITE(numout,*) '          global or zoom domain: jpiglo = ', jpiglo,   & 
    111             &                                            ' jpjglo = ', jpjglo, ' jpk    = ', jpk 
    112          WRITE(numout,*) '          local domain         : jpi    = ', jpi   ,   & 
    113             &                                            ' jpj    = ', jpj   , ' jpk    = ', jpk 
     97         WRITE(numout,*) '   global domain:   jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpkglo = ', jpkglo 
     98         WRITE(numout,*) '   local  domain:   jpi    = ', jpi   , ' jpj    = ', jpj   , ' jpk    = ', jpk 
    11499         WRITE(numout,*) 
    115          WRITE(numout,*) '          south-west indices    jpizoom = ', jpizoom,   & 
    116             &                                           ' jpjzoom = ', jpjzoom 
     100         WRITE(numout,*) '   conversion from local to global domain indices (and vise versa) done' 
    117101         IF( nn_print >= 1 ) THEN 
    118102            WRITE(numout,*) 
    119             WRITE(numout,*) '          conversion local  ==> data i-index domain' 
     103            WRITE(numout,*) '          conversion local  ==> global i-index domain (mig)' 
    120104            WRITE(numout,25)              (mig(ji),ji = 1,jpi) 
    121105            WRITE(numout,*) 
    122             WRITE(numout,*) '          conversion data  ==> local  i-index domain' 
    123             WRITE(numout,*) '             starting index' 
    124             WRITE(numout,25)              (mi0(ji),ji = 1,jpidta) 
    125             WRITE(numout,*) '             ending index' 
    126             WRITE(numout,25)              (mi1(ji),ji = 1,jpidta) 
     106            WRITE(numout,*) '          conversion global ==> local  i-index domain' 
     107            WRITE(numout,*) '             starting index (mi0)' 
     108            WRITE(numout,25)              (mi0(ji),ji = 1,jpiglo) 
     109            WRITE(numout,*) '             ending index (mi1)' 
     110            WRITE(numout,25)              (mi1(ji),ji = 1,jpiglo) 
    127111            WRITE(numout,*) 
    128             WRITE(numout,*) '          conversion local  ==> data j-index domain' 
     112            WRITE(numout,*) '          conversion local  ==> global j-index domain (mjg)' 
    129113            WRITE(numout,25)              (mjg(jj),jj = 1,jpj) 
    130114            WRITE(numout,*) 
    131             WRITE(numout,*) '          conversion data  ==> local j-index domain' 
    132             WRITE(numout,*) '             starting index' 
    133             WRITE(numout,25)              (mj0(jj),jj = 1,jpjdta) 
    134             WRITE(numout,*) '             ending index' 
    135             WRITE(numout,25)              (mj1(jj),jj = 1,jpjdta) 
     115            WRITE(numout,*) '          conversion global ==> local j-index domain' 
     116            WRITE(numout,*) '             starting index (mj0)' 
     117            WRITE(numout,25)              (mj0(jj),jj = 1,jpjglo) 
     118            WRITE(numout,*) '             ending index (mj1)' 
     119            WRITE(numout,25)              (mj1(jj),jj = 1,jpjglo) 
    136120         ENDIF 
    137121      ENDIF 
    138122 25   FORMAT( 100(10x,19i4,/) ) 
    139  
    140       !                        ! ============== ! 
    141       !                        !  Zoom domain   ! 
    142       !                        ! ============== ! 
    143       !                              ! zoom control 
    144       IF( jpiglo + jpizoom - 1  >  jpidta .OR.   & 
    145           jpjglo + jpjzoom - 1  >  jpjdta      ) & 
    146           &   CALL ctl_stop( ' global or zoom domain exceed the data domain ! ' ) 
    147  
    148       !                              ! set zoom flag 
    149       IF( jpiglo < jpidta .OR. jpjglo < jpjdta )   lzoom = .TRUE. 
    150  
    151       !                              ! set zoom type flags 
    152       IF( lzoom .AND. jpizoom /= 1 )   lzoom_w = .TRUE.                     !  
    153       IF( lzoom .AND. jpjzoom /= 1 )   lzoom_s = .TRUE. 
    154       IF( lzoom .AND. jpiglo + jpizoom -1 /= jpidta )   lzoom_e = .TRUE. 
    155       IF( lzoom .AND. jpjglo + jpjzoom -1 /= jpjdta )   lzoom_n = .TRUE. 
    156       IF(lwp) THEN 
    157          WRITE(numout,*) 
    158          WRITE(numout,*) '          zoom flags : ' 
    159          WRITE(numout,*) '             lzoom   = ', lzoom  , ' (T = zoom, F = global )' 
    160          WRITE(numout,*) '             lzoom_e = ', lzoom_e, ' (T = forced closed east  boundary)' 
    161          WRITE(numout,*) '             lzoom_w = ', lzoom_w, ' (T = forced closed west  boundary)' 
    162          WRITE(numout,*) '             lzoom_s = ', lzoom_s, ' (T = forced closed South boundary)' 
    163          WRITE(numout,*) '             lzoom_n = ', lzoom_n, ' (T = forced closed North boundary)' 
    164       ENDIF 
    165       IF(  ( lzoom_e .OR. lzoom_w )  .AND.  ( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 )  )   & 
    166            &   CALL ctl_stop( ' Your zoom choice is inconsistent with east-west cyclic boundary condition' ) 
    167       IF(  lzoom_n  .AND.  ( 3 <= jperio .AND. jperio <= 6 )  )   & 
    168            &   CALL ctl_stop( ' Your zoom choice is inconsistent with North fold boundary condition' ) 
    169  
    170       !                              ! Pre-defined arctic/antarctic zoom of ORCA configuration flag 
    171       IF( cp_cfg == "orca" ) THEN 
    172          SELECT CASE ( jp_cfg ) 
    173          CASE ( 2 )                               !  ORCA_R2 configuration 
    174             IF(  cp_cfz == "arctic"    .AND. jpiglo  == 142    .AND. jpjglo  ==  53 .AND.   & 
    175                & jpizoom ==  21    .AND. jpjzoom ==  97         )   THEN 
    176               IF(lwp) WRITE(numout,*) '          ORCA configuration: arctic zoom ' 
    177             ENDIF 
    178             IF(  cp_cfz == "antarctic" .AND. jpiglo  == jpidta .AND. jpjglo  ==  50 .AND.   & 
    179                & jpizoom ==   1    .AND. jpjzoom ==   1         )   THEN 
    180               IF(lwp) WRITE(numout,*) '          ORCA configuration: antarctic zoom ' 
    181             ENDIF 
    182             !                              
    183          CASE ( 05 )                              !  ORCA_R05 configuration 
    184             IF(    cp_cfz == "arctic"    .AND. jpiglo  == 562    .AND. jpjglo  == 202 .AND.   & 
    185                & jpizoom ==  81    .AND. jpjzoom == 301         )   THEN 
    186               IF(lwp) WRITE(numout,*) '          ORCA configuration: arctic zoom ' 
    187             ENDIF 
    188             IF(    cp_cfz == "antarctic" .AND. jpiglo  == jpidta .AND. jpjglo  == 187 .AND.   & 
    189                & jpizoom ==   1    .AND. jpjzoom ==   1         )   THEN 
    190               IF(lwp) WRITE(numout,*) '          ORCA configuration: antarctic zoom ' 
    191             ENDIF 
    192          END SELECT 
    193          ! 
    194       ENDIF 
    195123      ! 
    196124   END SUBROUTINE dom_glo 
    197  
    198125   !!====================================================================== 
    199126END MODULE domcfg 

  • utils/tools_UKMO_MERGE_2019/DOMAINcfg/src/domhgr.F90

    r12100 r12101  
    2828   USE in_out_manager ! I/O manager 
    2929   USE lib_mpp        ! MPP library 
    30    USE timing         ! Timing 
    3130 
    3231   IMPLICIT NONE 
    … …  
    111110      INTEGER  ::   ie1e2u_v             ! fag for u- & v-surface read in coordinate file or not 
    112111      !!---------------------------------------------------------------------- 
    113       ! 
    114   !    IF( nn_timing == 1 )  CALL timing_start('dom_hgr') 
    115112      ! 
    116113      IF(lwp) THEN 
    … …  
    437434      ! ------------------------------------------ 
    438435      ! The equator line must be the latitude coordinate axe 
    439  
    440 !      IF( nperio == 2 ) THEN 
    441 !         znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / REAL( jpi ) 
    442 !         IF( znorme > 1.e-13 ) CALL ctl_stop( ' ===>>>> : symmetrical condition: rerun with good equator line' ) 
    443 !      ENDIF 
    444       ! 
    445     !  IF( nn_timing == 1 )  CALL timing_stop('dom_hgr') 
     436! (PM) be carefull with nperio/jperio  
     437      IF( jperio == 2 ) THEN 
     438         znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / REAL( jpi ) 
     439         IF( znorme > 1.e-13 ) CALL ctl_stop( ' ===>>>> : symmetrical condition: rerun with good equator line' ) 
     440      ENDIF 
    446441      ! 
    447442   END SUBROUTINE dom_hgr 

  • utils/tools_UKMO_MERGE_2019/DOMAINcfg/src/dommsk.F90

    r12100 r12101  
    99   !!             -   ! 1996-05  (G. Madec)  mask computed from tmask 
    1010   !!            8.0  ! 1997-02  (G. Madec)  mesh information put in domhgr.F 
    11    !!            8.1  ! 1997-07  (G. Madec)  modification of mbathy and fmask 
     11   !!            8.1  ! 1997-07  (G. Madec)  modification of kbat and fmask 
    1212   !!             -   ! 1998-05  (G. Roullet)  free surface 
    1313   !!            8.2  ! 2000-03  (G. Madec)  no slip accurate 
    … …  
    1717   !!            3.2  ! 2009-07  (R. Benshila) Suppression of rigid-lid option 
    1818   !!            3.6  ! 2015-05  (P. Mathiot) ISF: add wmask,wumask and wvmask 
    19    !!---------------------------------------------------------------------- 
    20  
    21    !!---------------------------------------------------------------------- 
    22    !!   dom_msk        : compute land/ocean mask 
    23    !!---------------------------------------------------------------------- 
    24    USE oce             ! ocean dynamics and tracers 
    25    USE dom_oce         ! ocean space and time domain 
     19   !!            4.0  ! 2016-06  (G. Madec, S. Flavoni)  domain configuration / user defined interface 
     20   !!---------------------------------------------------------------------- 
     21 
     22   !!---------------------------------------------------------------------- 
     23   !!   dom_msk       : compute land/ocean mask 
     24   !!---------------------------------------------------------------------- 
     25   USE dom_oce        ! ocean space and time domain 
     26   USE domisf         ! domain: ice shelf 
     27   USE domwri         ! domain: write the meshmask file 
     28   USE bdy_oce        ! open boundary 
    2629   ! 
    27    USE in_out_manager  ! I/O manager 
    28    USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    29    USE lib_mpp         ! 
    30    USE wrk_nemo        ! Memory allocation 
    31    USE timing          ! Timing 
     30   USE in_out_manager ! I/O manager 
     31   USE iom            ! IOM library 
     32   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
     33   USE lib_mpp        ! Massively Parallel Processing library 
    3234 
    3335   IMPLICIT NONE 
    … …  
    4244 
    4345   !! * Substitutions 
    44    !!---------------------------------------------------------------------- 
    45    !!                   ***  vectopt_loop_substitute  *** 
    46    !!---------------------------------------------------------------------- 
    47    !! ** purpose :   substitute the inner loop start/end indices with CPP macro 
    48    !!                allow unrolling of do-loop (useful with vector processors) 
    49    !!---------------------------------------------------------------------- 
    50    !!---------------------------------------------------------------------- 
    51    !! NEMO/OPA 3.7 , NEMO Consortium (2014) 
    52    !! $Id: vectopt_loop_substitute.h90 4990 2014-12-15 16:42:49Z timgraham $  
    53    !! Software governed by the CeCILL licence (./LICENSE) 
    54    !!---------------------------------------------------------------------- 
    55    !!---------------------------------------------------------------------- 
    56    !! NEMO/OPA 3.2 , LODYC-IPSL  (2009) 
    57    !! $Id: dommsk.F90 6140 2015-12-21 11:35:23Z timgraham $  
    58    !! Software governed by the CeCILL licence     (./LICENSE) 
     46#  include "vectopt_loop_substitute.h90" 
     47   !!---------------------------------------------------------------------- 
     48   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     49   !! $Id: dommsk.F90 10425 2018-12-19 21:54:16Z smasson $  
     50   !! Software governed by the CeCILL license (see ./LICENSE) 
    5951   !!---------------------------------------------------------------------- 
    6052CONTAINS 
    … …  
    6759      !!      zontal velocity points (u & v), vorticity points (f) points. 
    6860      !! 
    69       !! ** Method  :   The ocean/land mask is computed from the basin bathy- 
    70       !!      metry in level (mbathy) which is defined or read in dommba. 
    71       !!      mbathy equals 0 over continental T-point  
    72       !!      and the number of ocean level over the ocean. 
    73       !! 
    74       !!      At a given position (ji,jj,jk) the ocean/land mask is given by: 
    75       !!      t-point : 0. IF mbathy( ji ,jj) =< 0 
    76       !!                1. IF mbathy( ji ,jj) >= jk 
    77       !!      u-point : 0. IF mbathy( ji ,jj)  or mbathy(ji+1, jj ) =< 0 
    78       !!                1. IF mbathy( ji ,jj) and mbathy(ji+1, jj ) >= jk. 
    79       !!      v-point : 0. IF mbathy( ji ,jj)  or mbathy( ji ,jj+1) =< 0 
    80       !!                1. IF mbathy( ji ,jj) and mbathy( ji ,jj+1) >= jk. 
    81       !!      f-point : 0. IF mbathy( ji ,jj)  or mbathy( ji ,jj+1) 
    82       !!                   or mbathy(ji+1,jj)  or mbathy(ji+1,jj+1) =< 0 
    83       !!                1. IF mbathy( ji ,jj) and mbathy( ji ,jj+1) 
    84       !!                  and mbathy(ji+1,jj) and mbathy(ji+1,jj+1) >= jk. 
    85       !!      tmask_i : interior ocean mask at t-point, i.e. excluding duplicated 
    86       !!                rows/lines due to cyclic or North Fold boundaries as well 
    87       !!                as MPP halos. 
    88       !! 
    89       !!        The lateral friction is set through the value of fmask along 
    90       !!      the coast and topography. This value is defined by rn_shlat, a 
    91       !!      namelist parameter: 
     61      !! ** Method  :   The ocean/land mask  at t-point is deduced from ko_top  
     62      !!      and ko_bot, the indices of the fist and last ocean t-levels which  
     63      !!      are either defined in usrdef_zgr or read in zgr_read. 
     64      !!                The velocity masks (umask, vmask)  
     65      !!      are deduced from a product of the two neighboring tmask. 
     66      !!                The vorticity mask (fmask) is deduced from tmask taking 
     67      !!      into account the choice of lateral boundary condition (rn_shlat) : 
    9268      !!         rn_shlat = 0, free slip  (no shear along the coast) 
    9369      !!         rn_shlat = 2, no slip  (specified zero velocity at the coast) 
    … …  
    9571      !!         2 < rn_shlat, strong slip        | in the lateral boundary layer 
    9672      !! 
    97       !!      N.B. If nperio not equal to 0, the land/ocean mask arrays 
    98       !!      are defined with the proper value at lateral domain boundaries. 
    99       !! 
    100       !!      In case of open boundaries (lk_bdy=T): 
    101       !!        - tmask is set to 1 on the points to be computed bay the open 
    102       !!          boundaries routines. 
    103       !! 
    104       !! ** Action :   tmask    : land/ocean mask at t-point (=0. or 1.) 
    105       !!               umask    : land/ocean mask at u-point (=0. or 1.) 
    106       !!               vmask    : land/ocean mask at v-point (=0. or 1.) 
    107       !!               fmask    : land/ocean mask at f-point (=0. or 1.) 
    108       !!                          =rn_shlat along lateral boundaries 
    109       !!               tmask_i  : interior ocean mask 
     73      !!      tmask_i : interior ocean mask at t-point, i.e. excluding duplicated 
     74      !!                rows/lines due to cyclic or North Fold boundaries as well 
     75      !!                as MPP halos. 
     76      !!      tmask_h : halo mask at t-point, i.e. excluding duplicated rows/lines 
     77      !!                due to cyclic or North Fold boundaries as well as MPP halos. 
     78      !! 
     79      !! ** Action :   tmask, umask, vmask, wmask : land/ocean mask  
     80      !!                         at t-, u-, v- w, wu-, and wv-points (=0. or 1.) 
     81      !!               fmask   : land/ocean mask at f-point (=0., or =1., or  
     82      !!                         =rn_shlat along lateral boundaries) 
     83      !!               tmask_i : interior ocean mask  
     84      !!               tmask_h : halo mask 
     85      !!               ssmask , ssumask, ssvmask, ssfmask : 2D ocean mask 
    11086      !!---------------------------------------------------------------------- 
    111       INTEGER  ::   ji, jj, jk               ! dummy loop indices 
    112       INTEGER  ::   iif, iil, ii0, ii1, ii   ! local integers 
    113       INTEGER  ::   ijf, ijl, ij0, ij1       !   -       - 
    114       INTEGER  ::   ios 
    115       INTEGER  ::   isrow                    ! index for ORCA1 starting row 
    116       INTEGER , POINTER, DIMENSION(:,:) ::  imsk 
    117       REAL(wp), POINTER, DIMENSION(:,:) ::  zwf 
     87      ! 
     88      INTEGER  ::   ji, jj, jk     ! dummy loop indices 
     89      INTEGER  ::   iif, iil       ! local integers 
     90      INTEGER  ::   ijf, ijl       !   -       - 
     91      INTEGER  ::   iktop, ikbot   !   -       - 
     92      INTEGER  ::   ios, inum 
     93      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zwf   ! 2D workspace 
    11894      !! 
    11995      NAMELIST/namlbc/ rn_shlat, ln_vorlat 
     96      NAMELIST/nambdy/ ln_bdy ,nb_bdy, ln_coords_file, cn_coords_file,         & 
     97         &             ln_mask_file, cn_mask_file, cn_dyn2d, nn_dyn2d_dta,     & 
     98         &             cn_dyn3d, nn_dyn3d_dta, cn_tra, nn_tra_dta,             & 
     99         &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, & 
     100         &             cn_ice, nn_ice_dta,                                     & 
     101         &             rn_ice_tem, rn_ice_sal, rn_ice_age,                     & 
     102         &             ln_vol, nn_volctl, nn_rimwidth, nb_jpk_bdy 
    120103      !!--------------------------------------------------------------------- 
    121       ! 
    122   !    IF( nn_timing == 1 )  CALL timing_start('dom_msk') 
    123       ! 
    124       CALL wrk_alloc( jpi, jpj, imsk ) 
    125       CALL wrk_alloc( jpi, jpj, zwf  ) 
    126104      ! 
    127105      REWIND( numnam_ref )              ! Namelist namlbc in reference namelist : Lateral momentum boundary condition 
    128106      READ  ( numnam_ref, namlbc, IOSTAT = ios, ERR = 901 ) 
    129 901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namlbc in reference namelist', lwp ) 
    130  
     107901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namlbc in reference namelist', lwp ) 
    131108      REWIND( numnam_cfg )              ! Namelist namlbc in configuration namelist : Lateral momentum boundary condition 
    132109      READ  ( numnam_cfg, namlbc, IOSTAT = ios, ERR = 902 ) 
    133 902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namlbc in configuration namelist', lwp ) 
     110902   IF( ios >  0 )  CALL ctl_nam ( ios , 'namlbc in configuration namelist', lwp ) 
    134111      IF(lwm) WRITE ( numond, namlbc ) 
    135112       
    … …  
    142119         WRITE(numout,*) '      consistency with analytical form   ln_vorlat = ',ln_vorlat  
    143120      ENDIF 
    144  
    145       IF     (      rn_shlat == 0.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  free-slip ' 
    146       ELSEIF (      rn_shlat == 2.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  no-slip ' 
    147       ELSEIF ( 0. < rn_shlat .AND. rn_shlat < 2. ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  partial-slip ' 
    148       ELSEIF ( 2. < rn_shlat                     ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  strong-slip ' 
     121      ! 
     122      IF(lwp) WRITE(numout,*) 
     123      IF     (      rn_shlat == 0.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ==>>>   ocean lateral  free-slip' 
     124      ELSEIF (      rn_shlat == 2.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ==>>>   ocean lateral  no-slip' 
     125      ELSEIF ( 0. < rn_shlat .AND. rn_shlat < 2. ) THEN   ;   IF(lwp) WRITE(numout,*) '   ==>>>   ocean lateral  partial-slip' 
     126      ELSEIF ( 2. < rn_shlat                     ) THEN   ;   IF(lwp) WRITE(numout,*) '   ==>>>   ocean lateral  strong-slip' 
    149127      ELSE 
    150          WRITE(ctmp1,*) ' rn_shlat is negative = ', rn_shlat 
    151          CALL ctl_stop( ctmp1 ) 
    152       ENDIF 
    153  
    154       ! 1. Ocean/land mask at t-point (computed from mbathy) 
    155       ! ----------------------------- 
    156       ! N.B. tmask has already the right boundary conditions since mbathy is ok 
    157       ! 
     128         CALL ctl_stop( 'dom_msk: wrong value for rn_shlat (i.e. a negalive value). We stop.' ) 
     129      ENDIF 
     130 
     131     ! 1. Ocean/land mask at t-point (computed from mbathy) 
     132     ! ----------------------------- 
     133     ! N.B. tmask has already the right boundary conditions since mbathy is ok 
     134     ! 
    158135      tmask(:,:,:) = 0._wp 
    159136      DO jk = 1, jpk 
    160137         DO jj = 1, jpj 
    161138            DO ji = 1, jpi 
    162                IF( REAL( mbathy(ji,jj) - jk, wp ) + 0.1_wp >= 0._wp )   tmask(ji,jj,jk) = 1._wp 
    163             END DO   
    164          END DO   
    165       END DO   
    166        
    167       ! (ISF) define barotropic mask and mask the ice shelf point 
    168       ssmask(:,:)=tmask(:,:,1) ! at this stage ice shelf is not masked 
    169        
    170       DO jk = 1, jpk 
    171          DO jj = 1, jpj 
    172             DO ji = 1, jpi 
    173                IF( REAL( misfdep(ji,jj) - jk, wp ) - 0.1_wp >= 0._wp )   THEN 
    174                   tmask(ji,jj,jk) = 0._wp 
    175                END IF 
    176             END DO   
    177          END DO   
    178       END DO   
    179  
    180       ! Interior domain mask (used for global sum) 
    181       ! -------------------- 
    182    !   tmask_i(:,:) = ssmask(:,:)            ! (ISH) tmask_i = 1 even on the ice shelf 
    183  
    184   !    tmask_h(:,:) = 1._wp                 ! 0 on the halo and 1 elsewhere 
    185   !    iif = jpreci                         ! ??? 
    186   !    iil = nlci - jpreci + 1 
    187   !    ijf = jprecj                         ! ??? 
    188   !    ijl = nlcj - jprecj + 1 
    189  
    190    !   tmask_h( 1 :iif,   :   ) = 0._wp      ! first columns 
    191    !   tmask_h(iil:jpi,   :   ) = 0._wp      ! last  columns (including mpp extra columns) 
    192    !   tmask_h(   :   , 1 :ijf) = 0._wp      ! first rows 
    193    !   tmask_h(   :   ,ijl:jpj) = 0._wp      ! last  rows (including mpp extra rows) 
    194  
    195       ! north fold mask 
    196       ! --------------- 
    197    !   tpol(1:jpiglo) = 1._wp  
    198    !   fpol(1:jpiglo) = 1._wp 
    199    !   IF( jperio == 3 .OR. jperio == 4 ) THEN      ! T-point pivot 
    200    !      tpol(jpiglo/2+1:jpiglo) = 0._wp 
    201    !      fpol(     1    :jpiglo) = 0._wp 
    202    !      IF( mjg(nlej) == jpjglo ) THEN                  ! only half of the nlcj-1 row 
    203    !         DO ji = iif+1, iil-1 
    204    !            tmask_h(ji,nlej-1) = tmask_h(ji,nlej-1) * tpol(mig(ji)) 
    205    !         END DO 
    206    !      ENDIF 
    207    !   ENDIF 
    208       
    209    !   tmask_i(:,:) = tmask_i(:,:) * tmask_h(:,:) 
    210  
    211   !    IF( jperio == 5 .OR. jperio == 6 ) THEN      ! F-point pivot 
    212   !       tpol(     1    :jpiglo) = 0._wp 
    213   !       fpol(jpiglo/2+1:jpiglo) = 0._wp 
    214   !    ENDIF 
    215  
    216       ! 2. Ocean/land mask at u-,  v-, and z-points (computed from tmask) 
    217       ! ------------------------------------------- 
     139               IF(      ( REAL( mbathy (ji,jj) - jk, wp ) + 0.1_wp >= 0._wp )         & 
     140               &  .AND. ( REAL( misfdep(ji,jj) - jk, wp ) - 0.1_wp <= 0._wp ) ) THEN 
     141                  tmask(ji,jj,jk) = 1._wp 
     142               END IF   
     143            END DO 
     144         END DO 
     145      END DO     
     146  
     147      IF ( ln_isfsubgl ) CALL zgr_isf_subgl 
     148 
     149!SF  add here lbc_lnk: bug not still understood : cause now domain configuration is read ! 
     150!!gm I don't understand why...   
     151      CALL lbc_lnk( 'dommsk', tmask  , 'T', 1._wp )      ! Lateral boundary conditions 
     152 
     153     ! Mask corrections for bdy (read in mppini2) 
     154      REWIND( numnam_ref )              ! Namelist nambdy in reference namelist :Unstructured open boundaries 
     155      READ  ( numnam_ref, nambdy, IOSTAT = ios, ERR = 903) 
     156903   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy in reference namelist', lwp ) 
     157      REWIND( numnam_cfg )              ! Namelist nambdy in configuration namelist :Unstructured open boundaries 
     158      READ  ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 904 ) 
     159904   IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy in configuration namelist', lwp ) 
     160      ! ------------------------ 
     161      IF ( ln_bdy .AND. ln_mask_file ) THEN 
     162         CALL iom_open( cn_mask_file, inum ) 
     163         CALL iom_get ( inum, jpdom_data, 'bdy_msk', bdytmask(:,:) ) 
     164         CALL iom_close( inum ) 
     165         DO jk = 1, jpkm1 
     166            DO jj = 1, jpj 
     167               DO ji = 1, jpi 
     168                  tmask(ji,jj,jk) = tmask(ji,jj,jk) * bdytmask(ji,jj) 
     169               END DO 
     170            END DO 
     171         END DO 
     172      ENDIF 
     173          
     174      !  Ocean/land mask at u-, v-, and f-points   (computed from tmask) 
     175      ! ---------------------------------------- 
     176      ! NB: at this point, fmask is designed for free slip lateral boundary condition 
    218177      DO jk = 1, jpk 
    219178         DO jj = 1, jpjm1 
    220             DO ji = 1, jpim1   ! vector loop 
     179            DO ji = 1, fs_jpim1   ! vector loop 
    221180               umask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk) 
    222181               vmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji  ,jj+1,jk) 
    … …  
    228187         END DO 
    229188      END DO 
    230       ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at least 1 wet cell at u point 
    231 !     DO jj = 1, jpjm1 
    232 !         DO ji = 1, jpim1   ! vector loop 
    233 !            ssumask(ji,jj)  = ssmask(ji,jj) * ssmask(ji+1,jj  )  * MIN(1._wp,SUM(umask(ji,jj,:))) 
    234 !            ssvmask(ji,jj)  = ssmask(ji,jj) * ssmask(ji  ,jj+1)  * MIN(1._wp,SUM(vmask(ji,jj,:))) 
    235 !!         END DO 
    236 !         DO ji = 1, jpim1      ! NO vector opt. 
    237 !            ssfmask(ji,jj) =  ssmask(ji,jj  ) * ssmask(ji+1,jj  )   & 
    238 !               &            * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:))) 
    239 !         END DO 
    240 !      END DO 
    241       CALL lbc_lnk( 'toto',umask  , 'U', 1._wp )      ! Lateral boundary conditions 
    242       CALL lbc_lnk( 'toto',vmask  , 'V', 1._wp ) 
    243       CALL lbc_lnk( 'toto',fmask  , 'F', 1._wp ) 
    244  !     CALL lbc_lnk( 'toto',ssumask, 'U', 1._wp )      ! Lateral boundary conditions 
    245  !     CALL lbc_lnk( 'toto',ssvmask, 'V', 1._wp ) 
    246  !     CALL lbc_lnk( 'toto',ssfmask, 'F', 1._wp ) 
    247  
    248       ! 3. Ocean/land mask at wu-, wv- and w points  
    249       !---------------------------------------------- 
     189      CALL lbc_lnk_multi( 'dommsk', umask, 'U', 1., vmask, 'V', 1., fmask, 'F', 1. )      ! Lateral boundary conditions 
     190  
     191      ! Ocean/land mask at wu-, wv- and w points    (computed from tmask) 
     192      !----------------------------------------- 
    250193      wmask (:,:,1) = tmask(:,:,1)     ! surface 
    251       wumask(:,:,1) = umask(:,:,1) 
    252       wvmask(:,:,1) = vmask(:,:,1) 
    253194      DO jk = 2, jpk                   ! interior values 
    254195         wmask (:,:,jk) = tmask(:,:,jk) * tmask(:,:,jk-1) 
    255          wumask(:,:,jk) = umask(:,:,jk) * umask(:,:,jk-1)    
    256          wvmask(:,:,jk) = vmask(:,:,jk) * vmask(:,:,jk-1) 
    257196      END DO 
    258197 
     198 
     199      ! Ocean/land column mask at t-, u-, and v-points   (i.e. at least 1 wet cell in the vertical) 
     200      ! ---------------------------------------------- 
     201      ssmask (:,:) = MAXVAL( tmask(:,:,:), DIM=3 ) 
     202      ssumask(:,:) = MAXVAL( umask(:,:,:), DIM=3 ) 
     203      ssvmask(:,:) = MAXVAL( vmask(:,:,:), DIM=3 ) 
     204 
     205      ! Interior domain mask  (used for global sum) 
     206      ! -------------------- 
     207      ! 
     208      iif = nn_hls   ;   iil = nlci - nn_hls + 1 
     209      ijf = nn_hls   ;   ijl = nlcj - nn_hls + 1 
     210      ! 
     211      !                          ! halo mask : 0 on the halo and 1 elsewhere 
     212      tmask_h(:,:) = 1._wp                   
     213      tmask_h( 1 :iif,   :   ) = 0._wp      ! first columns 
     214      tmask_h(iil:jpi,   :   ) = 0._wp      ! last  columns (including mpp extra columns) 
     215      tmask_h(   :   , 1 :ijf) = 0._wp      ! first rows 
     216      tmask_h(   :   ,ijl:jpj) = 0._wp      ! last  rows (including mpp extra rows) 
     217      ! 
     218      !                          ! north fold mask 
     219      tpol(1:jpiglo) = 1._wp  
     220      fpol(1:jpiglo) = 1._wp 
     221      IF( jperio == 3 .OR. jperio == 4 ) THEN      ! T-point pivot 
     222         tpol(jpiglo/2+1:jpiglo) = 0._wp 
     223         fpol(     1    :jpiglo) = 0._wp 
     224         IF( mjg(nlej) == jpjglo ) THEN                  ! only half of the nlcj-1 row for tmask_h 
     225            DO ji = iif+1, iil-1 
     226               tmask_h(ji,nlej-1) = tmask_h(ji,nlej-1) * tpol(mig(ji)) 
     227            END DO 
     228         ENDIF 
     229      ENDIF 
     230      ! 
     231      IF( jperio == 5 .OR. jperio == 6 ) THEN      ! F-point pivot 
     232         tpol(     1    :jpiglo) = 0._wp 
     233         fpol(jpiglo/2+1:jpiglo) = 0._wp 
     234      ENDIF 
     235      ! 
     236      !                          ! interior mask : 2D ocean mask x halo mask  
     237      tmask_i(:,:) = ssmask(:,:) * tmask_h(:,:) 
     238 
     239 
    259240      ! Lateral boundary conditions on velocity (modify fmask) 
    260       ! ---------------------------------------      
    261       DO jk = 1, jpk 
    262          zwf(:,:) = fmask(:,:,jk)          
    263          DO jj = 2, jpjm1 
    264             DO ji = 2, jpim1   ! vector opt. 
    265                IF( fmask(ji,jj,jk) == 0._wp ) THEN 
    266                   fmask(ji,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1),   & 
    267                      &                                           zwf(ji-1,jj), zwf(ji,jj-1)  )  ) 
    268                ENDIF 
    269             END DO 
    270          END DO 
    271          DO jj = 2, jpjm1 
    272             IF( fmask(1,jj,jk) == 0._wp ) THEN 
    273                fmask(1  ,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) ) 
    274             ENDIF 
    275             IF( fmask(jpi,jj,jk) == 0._wp ) THEN 
    276                fmask(jpi,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) ) 
    277             ENDIF 
    278          END DO          
    279          DO ji = 2, jpim1 
    280             IF( fmask(ji,1,jk) == 0._wp ) THEN 
    281                fmask(ji, 1 ,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) ) 
    282             ENDIF 
    283             IF( fmask(ji,jpj,jk) == 0._wp ) THEN 
    284                fmask(ji,jpj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) ) 
    285             ENDIF 
    286          END DO 
    287       END DO 
    288       ! 
    289       IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA_R2 configuration 
    290          !                                                 ! Increased lateral friction near of some straits 
    291          !                                ! Gibraltar strait  : partial slip (fmask=0.5) 
    292          ij0 = 101   ;   ij1 = 101 
    293          ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp 
    294          ij0 = 102   ;   ij1 = 102 
    295          ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp 
    296          ! 
    297          !                                ! Bab el Mandeb : partial slip (fmask=1) 
    298          ij0 =  87   ;   ij1 =  88 
    299          ii0 = 160   ;   ii1 = 160   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp 
    300          ij0 =  88   ;   ij1 =  88 
    301          ii0 = 159   ;   ii1 = 159   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp 
    302          ! 
    303          !                                ! Danish straits  : strong slip (fmask > 2) 
    304 ! We keep this as an example but it is instable in this case  
    305 !         ij0 = 115   ;   ij1 = 115 
    306 !         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp 
    307 !         ij0 = 116   ;   ij1 = 116 
    308 !         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp 
    309          ! 
    310       ENDIF 
    311       ! 
    312       IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN   ! ORCA R1 configuration 
    313          !                                                 ! Increased lateral friction near of some straits 
    314          ! This dirty section will be suppressed by simplification process: 
    315          ! all this will come back in input files 
    316          ! Currently these hard-wired indices relate to configuration with 
    317          ! extend grid (jpjglo=332) 
    318          ! 
    319          isrow = 332 - jpjglo 
    320          ! 
    321          IF(lwp) WRITE(numout,*) 
    322          IF(lwp) WRITE(numout,*) '   orca_r1: increase friction near the following straits : ' 
    323          IF(lwp) WRITE(numout,*) '      Gibraltar ' 
    324          ii0 = 282           ;   ii1 = 283        ! Gibraltar Strait  
    325          ij0 = 241 - isrow   ;   ij1 = 241 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp   
    326  
    327          IF(lwp) WRITE(numout,*) '      Bhosporus ' 
    328          ii0 = 314           ;   ii1 = 315        ! Bhosporus Strait  
    329          ij0 = 248 - isrow   ;   ij1 = 248 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp   
    330  
    331          IF(lwp) WRITE(numout,*) '      Makassar (Top) ' 
    332          ii0 =  48           ;   ii1 =  48        ! Makassar Strait (Top)  
    333          ij0 = 189 - isrow   ;   ij1 = 190 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp   
    334  
    335          IF(lwp) WRITE(numout,*) '      Lombok ' 
    336          ii0 =  44           ;   ii1 =  44        ! Lombok Strait  
    337          ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp   
    338  
    339          IF(lwp) WRITE(numout,*) '      Ombai ' 
    340          ii0 =  53           ;   ii1 =  53        ! Ombai Strait  
    341          ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp   
    342  
    343          IF(lwp) WRITE(numout,*) '      Timor Passage ' 
    344          ii0 =  56           ;   ii1 =  56        ! Timor Passage  
    345          ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp   
    346  
    347          IF(lwp) WRITE(numout,*) '      West Halmahera ' 
    348          ii0 =  58           ;   ii1 =  58        ! West Halmahera Strait  
    349          ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp   
    350  
    351          IF(lwp) WRITE(numout,*) '      East Halmahera ' 
    352          ii0 =  55           ;   ii1 =  55        ! East Halmahera Strait  
    353          ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp   
    354          ! 
    355       ENDIF 
    356       ! 
    357       CALL lbc_lnk( 'toto',fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask 
    358       ! 
    359       ! CAUTION : The fmask may be further modified in dyn_vor_init ( dynvor.F90 ) 
    360       ! 
    361       CALL wrk_dealloc( jpi, jpj, imsk ) 
    362       CALL wrk_dealloc( jpi, jpj, zwf  ) 
    363       ! 
    364   !    IF( nn_timing == 1 )  CALL timing_stop('dom_msk') 
     241      ! ---------------------------------------   
     242      IF( rn_shlat /= 0 ) THEN      ! Not free-slip lateral boundary condition 
     243         ! 
     244         ALLOCATE( zwf(jpi,jpj) ) 
     245         ! 
     246         DO jk = 1, jpk 
     247            zwf(:,:) = fmask(:,:,jk)          
     248            DO jj = 2, jpjm1 
     249               DO ji = fs_2, fs_jpim1   ! vector opt. 
     250                  IF( fmask(ji,jj,jk) == 0._wp ) THEN 
     251                     fmask(ji,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1),   & 
     252                        &                                           zwf(ji-1,jj), zwf(ji,jj-1)  )  ) 
     253                  ENDIF 
     254               END DO 
     255            END DO 
     256            DO jj = 2, jpjm1 
     257               IF( fmask(1,jj,jk) == 0._wp ) THEN 
     258                  fmask(1  ,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) ) 
     259               ENDIF 
     260               IF( fmask(jpi,jj,jk) == 0._wp ) THEN 
     261                  fmask(jpi,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) ) 
     262               ENDIF 
     263            END DO          
     264            DO ji = 2, jpim1 
     265               IF( fmask(ji,1,jk) == 0._wp ) THEN 
     266                  fmask(ji, 1 ,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) ) 
     267               ENDIF 
     268               IF( fmask(ji,jpj,jk) == 0._wp ) THEN 
     269                  fmask(ji,jpj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) ) 
     270               ENDIF 
     271            END DO 
     272#if defined key_agrif  
     273            IF( .NOT. AGRIF_Root() ) THEN  
     274               IF ((nbondi ==  1).OR.(nbondi == 2)) fmask(nlci-1 , :     ,jk) = 0.e0      ! east  
     275               IF ((nbondi == -1).OR.(nbondi == 2)) fmask(1      , :     ,jk) = 0.e0      ! west  
     276               IF ((nbondj ==  1).OR.(nbondj == 2)) fmask(:      ,nlcj-1 ,jk) = 0.e0      ! north  
     277               IF ((nbondj == -1).OR.(nbondj == 2)) fmask(:      ,1      ,jk) = 0.e0      ! south  
     278            ENDIF  
     279#endif  
     280         END DO 
     281         ! 
     282         DEALLOCATE( zwf ) 
     283         ! 
     284         CALL lbc_lnk( 'dommsk', fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask 
     285         ! 
     286         ! CAUTION : The fmask may be further modified in dyn_vor_init ( dynvor.F90 ) depending on ln_vorlat 
     287         ! 
     288      ENDIF 
     289      ! 
     290      ! write out mesh mask 
     291      IF ( nn_msh > 0 ) CALL dom_wri 
    365292      ! 
    366293   END SUBROUTINE dom_msk 

  • utils/tools_UKMO_MERGE_2019/DOMAINcfg/src/domngb.F90

    r12100 r12101  
    1111   !!---------------------------------------------------------------------- 
    1212   USE dom_oce        ! ocean space and time domain 
     13   USE phycst 
    1314   ! 
    1415   USE in_out_manager ! I/O manager 
    … …  
    1819   PRIVATE 
    1920 
    20    PUBLIC   dom_ngb   ! routine called in iom.F90 module 
     21   PUBLIC   dom_ngb   ! routine called in iom.F90 and domclo.F90 module 
     22   PUBLIC   dist 
    2123 
    2224   !!---------------------------------------------------------------------- 
    … …  
    2729CONTAINS 
    2830 
    29    SUBROUTINE dom_ngb( plon, plat, kii, kjj, cdgrid, kkk ) 
     31   SUBROUTINE dom_ngb( plon, plat, kii, kjj, rdist, cdgrid, kkk ) 
    3032      !!---------------------------------------------------------------------- 
    3133      !!                    ***  ROUTINE dom_ngb  *** 
    … …  
    3739      !!---------------------------------------------------------------------- 
    3840      REAL(wp)        , INTENT(in   ) ::   plon, plat   ! longitude,latitude of the point 
     41      REAL(wp)        , INTENT(  out) ::   rdist        ! distance between the located point and the source 
    3942      INTEGER         , INTENT(  out) ::   kii, kjj     ! i-,j-index of the closes grid point 
    4043      INTEGER         , INTENT(in   ), OPTIONAL :: kkk  ! k-index of the mask level used 
    … …  
    4346      INTEGER :: ik         ! working level 
    4447      INTEGER , DIMENSION(2) ::   iloc 
    45       REAL(wp)               ::   zlon, zmini 
    4648      REAL(wp), DIMENSION(jpi,jpj) ::   zglam, zgphi, zmask, zdist 
    4749      !!-------------------------------------------------------------------- 
    … …  
    5759      END SELECT 
    5860 
    59       zlon       = MOD( plon       + 720., 360. )                                     ! plon between    0 and 360 
    60       zglam(:,:) = MOD( zglam(:,:) + 720., 360. )                                     ! glam between    0 and 360 
    61       IF( zlon > 270. )   zlon = zlon - 360.                                          ! zlon between  -90 and 270 
    62       IF( zlon <  90. )   WHERE( zglam(:,:) > 180. ) zglam(:,:) = zglam(:,:) - 360.   ! glam between -180 and 180 
    63       zglam(:,:) = zglam(:,:) - zlon 
     61      zdist = dist(plon, plat, zglam, zgphi) 
    6462 
    65       zgphi(:,:) = zgphi(:,:) - plat 
    66       zdist(:,:) = zglam(:,:) * zglam(:,:) + zgphi(:,:) * zgphi(:,:) 
    67        
    6863      IF( lk_mpp ) THEN   
    69          CALL mpp_minloc( 'domngb', zdist(:,:), zmask, zmini, iloc) 
     64         CALL mpp_minloc( 'domngb', zdist(:,:), zmask, rdist, iloc) 
    7065         kii = iloc(1) ; kjj = iloc(2) 
    7166      ELSE 
    7267         iloc(:) = MINLOC( zdist(:,:), mask = zmask(:,:) == 1.e0 ) 
     68         rdist = MINVAL( zdist(:,:) ) 
    7369         kii = iloc(1) + nimpp - 1 
    7470         kjj = iloc(2) + njmpp - 1 
    … …  
    7773   END SUBROUTINE dom_ngb 
    7874 
     75   FUNCTION dist(plonsrc, platsrc, plontrg, plattrg) 
     76      REAL(wp),                     INTENT(in) :: plonsrc, platsrc ! lat/lon of the source point 
     77      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: plontrg, plattrg ! lat/lon of the target (2d in this case) 
     78 
     79      REAL(wp) :: zxs, zys, zzs 
     80      REAL(wp), DIMENSION(jpi,jpj) :: zxt, zyt, zzt 
     81 
     82      REAL(wp), DIMENSION(jpi,jpj) :: dist ! distance between src and trg 
     83 
     84      zxs = COS( rad * platsrc ) * COS( rad * plonsrc ) 
     85      zys = COS( rad * platsrc ) * SIN( rad * plonsrc ) 
     86      zzs = SIN( rad * platsrc ) 
     87 
     88      zxt = COS( rad * plattrg ) * COS( rad * plontrg ) 
     89      zyt = COS( rad * plattrg ) * SIN( rad * plontrg ) 
     90      zzt = SIN( rad * plattrg ) 
     91 
     92      dist(:,:) = ( zxs - zxt(:,:) )**2   & 
     93         &      + ( zys - zyt(:,:) )**2   & 
     94         &      + ( zzs - zzt(:,:) )**2 
     95 
     96      dist(:,:) = ra * SQRT( dist(:,:) ) 
     97 
     98   END FUNCTION dist 
     99 
    79100   !!====================================================================== 
    80101END MODULE domngb 

  • utils/tools_UKMO_MERGE_2019/DOMAINcfg/src/domwri.F90

    r12100 r12101  
    154154      CALL iom_rstput( 0, 0, inum, 'ff_t', ff_t, ktype = jp_r8 ) 
    155155       
    156       ! note that mbkt is set to 1 over land ==> use surface tmask 
     156      ! note that mbkt and mikt is set to 1 over land ==> use surface tmask 
    157157      zprt(:,:) = ssmask(:,:) * REAL( mbkt(:,:) , wp ) 
    158       CALL iom_rstput( 0, 0, inum, 'mbathy', zprt, ktype = jp_i4 )     !    ! nb of ocean T-points 
     158      CALL iom_rstput( 0, 0, inum, 'mbathy', zprt, ktype = jp_i4 )                              !    ! nb of ocean T-points 
     159      CALL iom_rstput( 0, 0, inum, 'bathy_metry', bathy(:,:)   * ssmask(:,:), ktype = jp_r8 )   !    ! bathymetry 
    159160      zprt(:,:) = ssmask(:,:) * REAL( mikt(:,:) , wp ) 
    160       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 
    163       !                                                         ! vertical mesh 
     161      CALL iom_rstput( 0, 0, inum, 'misf', zprt, ktype = jp_i4 )                                !    ! first wet level 
     162      CALL iom_rstput( 0, 0, inum, 'isfdraft'   , risfdep(:,:) * ssmask(:,:), ktype = jp_r8 )   !    ! ice shelf draft 
     163      zprt(:,:) = ssmask(:,:) * REAL( mbkt(:,:) - mikt(:,:) + 1, wp ) 
     164      CALL iom_rstput( 0, 0, inum, 'mhw',zprt, ktype = jp_i4 ) 
     165      CALL iom_rstput( 0, 0, inum, 'hw' ,(bathy-risfdep)*ssmask, ktype = jp_r8 ) 
     166 
     167      !                  ! vertical mesh 
    164168      CALL iom_rstput( 0, 0, inum, 'e3t_0', e3t_0, ktype = jp_r8  )    !    ! scale factors 
    165169      CALL iom_rstput( 0, 0, inum, 'e3u_0', e3u_0, ktype = jp_r8  ) 
    … …  
    177181      ENDIF 
    178182      ! 
    179    !   IF( ll_wd ) CALL iom_rstput( 0, 0, inum, 'ht_0'   , ht_0   , ktype = jp_r8 ) 
    180  
    181183      !                                     ! ============================ 
    182184      CALL iom_close( inum )                !        close the files  

  • utils/tools_UKMO_MERGE_2019/DOMAINcfg/src/domzgr.F90

    r12100 r12101  
    1717   !!            3.4  ! 2012-08  (J. Siddorn) added Siddorn and Furner stretching function 
    1818   !!            3.4  ! 2012-12  (R. Bourdalle-Badie and G. Reffray)  modify C1D case   
    19    !!            3.6  ! 2014-11  (P. Mathiot and C. Harris) add ice shelf capabilitye 
     19   !!            3.6  ! 2014-11  (P. Mathiot and C. Harris) add ice shelf capability 
    2020   !!            3.?  ! 2015-11  (H. Liu) Modifications for Wetting/Drying 
    2121   !!---------------------------------------------------------------------- 
    … …  
    3535   !!       fgamma       : Siddorn and Furner 2012 stretching function 
    3636   !!--------------------------------------------------------------------- 
    37    USE oce               ! ocean variables 
    3837   USE dom_oce           ! ocean domain 
    39 !   USE closea            ! closed seas 
    4038   ! 
    4139   USE in_out_manager    ! I/O manager 
    … …  
    4341   USE lbclnk            ! ocean lateral boundary conditions (or mpp link) 
    4442   USE lib_mpp           ! distributed memory computing library 
    45    USE wrk_nemo          ! Memory allocation 
    46    USE timing            ! Timing 
     43   USE lib_fortran 
    4744   USE dombat 
     45   USE domisf 
    4846 
    4947   IMPLICIT NONE 
    … …  
    5149 
    5250   PUBLIC   dom_zgr        ! called by dom_init.F90 
    53  
     51   ! 
    5452   !                              !!* Namelist namzgr_sco * 
    5553   LOGICAL  ::   ln_s_sh94         ! use hybrid s-sig Song and Haidvogel 1994 stretching function fssig1 (ln_sco=T) 
    … …  
    6058   REAL(wp) ::   rn_rmax           ! maximum cut-off r-value allowed (0<rn_rmax<1) 
    6159   REAL(wp) ::   rn_hc             ! Critical depth for transition from sigma to stretched coordinates 
    62    INTEGER , PUBLIC ::   ntopo           !: = 0/1 ,compute/read the bathymetry file 
    63    REAL(wp), PUBLIC ::   e3zps_min       !: miminum thickness for partial steps (meters) 
    64    REAL(wp), PUBLIC ::   e3zps_rat       !: minimum thickness ration for partial steps 
    65    INTEGER, PUBLIC ::   nperio            !: type of lateral boundary condition 
    6660 
    6761   ! Song and Haidvogel 1994 stretching parameters 
    … …  
    121115      !!---------------------------------------------------------------------- 
    122116      ! 
    123   !    IF( nn_timing == 1 )   CALL timing_start('dom_zgr') 
    124117      ! 
    125118      REWIND( numnam_ref )              ! Namelist namzgr in reference namelist : Vertical coordinate 
    … …  
    152145      IF( ioptio /= 1 )   CALL ctl_stop( ' none or several vertical coordinate options used' ) 
    153146      ! 
     147      IF ( ln_isfcav ) CALL zgr_isf_nam 
    154148      ioptio = 0 
    155149      IF ( ln_zco .AND. ln_isfcav ) ioptio = ioptio + 1 
    … …  
    164158      IF( ln_zps      )   CALL zgr_zps          ! Partial step z-coordinate 
    165159      IF( ln_sco      )   CALL zgr_sco          ! s-coordinate or hybrid z-s coordinate 
     160                          CALL zgr_bat_ctl      ! check bathymetry (mbathy) and suppress isolated ocean points 
    166161      ! 
    167162      ! final adjustment of mbathy & check  
    168163      ! ----------------------------------- 
    169       IF( lzoom       )   CALL zgr_bat_zoom     ! correct mbathy in case of zoom subdomain 
    170                           CALL zgr_bat_ctl      ! check bathymetry (mbathy) and suppress isolated ocean points 
    171164                          CALL zgr_bot_level    ! deepest ocean level for t-, u- and v-points 
    172165                          CALL zgr_top_level    ! shallowest ocean level for T-, U-, V- points 
    … …  
    175168         WRITE(numout,*) ' MIN val mbathy  ', MINVAL(  mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 
    176169         WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ),   & 
    177             &                          ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ', MINVAL( gde3w_0(:,:,:) ) 
     170            &                          ' w ', MINVAL( gdepw_0(:,:,:) ) 
    178171         WRITE(numout,*) ' MIN val e3    t ', MINVAL(   e3t_0(:,:,:) ), ' f ', MINVAL(   e3f_0(:,:,:) ),  & 
    179172            &                          ' u ', MINVAL(   e3u_0(:,:,:) ), ' u ', MINVAL(   e3v_0(:,:,:) ),  & 
    … …  
    182175 
    183176         WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ),   & 
    184             &                          ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ', MAXVAL( gde3w_0(:,:,:) ) 
     177            &                          ' w ', MAXVAL( gdepw_0(:,:,:) ) 
    185178         WRITE(numout,*) ' MAX val e3    t ', MAXVAL(   e3t_0(:,:,:) ), ' f ', MAXVAL(   e3f_0(:,:,:) ),  & 
    186179            &                          ' u ', MAXVAL(   e3u_0(:,:,:) ), ' u ', MAXVAL(   e3v_0(:,:,:) ),  & 
    … …  
    188181            &                          ' w ', MAXVAL(   e3w_0(:,:,:) ) 
    189182      ENDIF 
    190       ! 
    191     !  IF( nn_timing == 1 )  CALL timing_stop('dom_zgr') 
    192183      ! 
    193184   END SUBROUTINE dom_zgr 
    … …  
    222213      REAL(wp) ::   za2, zkth2, zacr2      ! Values for optional double tanh function set from parameters  
    223214      !!---------------------------------------------------------------------- 
    224       ! 
    225    !   IF( nn_timing == 1 )  CALL timing_start('zgr_z') 
    226215      ! 
    227216      ! Set variables from parameters 
    … …  
    322311      IF ( ln_isfcav .OR. ln_e3_dep ) THEN      ! e3. = dk[gdep]    
    323312         ! 
    324 !==>>>   need to be like this to compute the pressure gradient with ISF.  
    325 !        If not, level beneath the ISF are not aligned (sum(e3t) /= depth) 
    326 !        define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively 
    327 ! 
    328313         DO jk = 1, jpkm1 
    329314            e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk)  
    … …  
    354339         IF( gdepw_1d(jk) <  0._wp .OR. gdept_1d(jk) <  0._wp )   CALL ctl_stop( 'dom:zgr_z: gdepw_1d or gdept_1d < 0 ' ) 
    355340      END DO 
    356       ! 
    357    !   IF( nn_timing == 1 )  CALL timing_stop('zgr_z') 
    358341      ! 
    359342   END SUBROUTINE zgr_z 
    … …  
    401384      !!---------------------------------------------------------------------- 
    402385      ! 
    403    !   IF( nn_timing == 1 )  CALL timing_start('zgr_bat') 
    404       ! 
    405386      IF(lwp) WRITE(numout,*) 
    406387      IF(lwp) WRITE(numout,*) '    zgr_bat : defines level and meter bathymetry' 
    … …  
    411392         !                                            ! global domain level and meter bathymetry (idta,zdta) 
    412393         ! 
    413          ALLOCATE( idta(jpidta,jpjdta), STAT=ierror ) 
     394         ALLOCATE( idta(jpiglo,jpjglo), STAT=ierror ) 
    414395         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate idta array' ) 
    415          ALLOCATE( zdta(jpidta,jpjdta), STAT=ierror ) 
     396         ALLOCATE( zdta(jpiglo,jpjglo), STAT=ierror ) 
    416397         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate zdta array' ) 
    417398         ! 
    … …  
    439420            IF(lwp) WRITE(numout,*) 
    440421            IF(lwp) WRITE(numout,*) '         bathymetry field: flat basin with a bump' 
    441             ii_bump = jpidta / 2                           ! i-index of the bump center 
    442             ij_bump = jpjdta / 2                           ! j-index of the bump center 
     422            ii_bump = jpiglo / 2                           ! i-index of the bump center 
     423            ij_bump = jpjglo / 2                           ! j-index of the bump center 
    443424            r_bump  = 50000._wp                            ! bump radius (meters)        
    444425            h_bump  =  2700._wp                            ! bump height (meters) 
    … …  
    450431            IF(lwp) WRITE(numout,*) '            background ocean depth = ', h_oce  , ' meters' 
    451432            !                                         
    452             DO jj = 1, jpjdta                              ! zdta : 
    453                DO ji = 1, jpidta 
     433            DO jj = 1, jpjglo                              ! zdta : 
     434               DO ji = 1, jpiglo 
    454435                  zi = FLOAT( ji - ii_bump ) * ppe1_m / r_bump 
    455436                  zj = FLOAT( jj - ij_bump ) * ppe2_m / r_bump 
    … …  
    467448            ENDIF 
    468449         ENDIF 
     450         ! 
    469451         !                                            ! set GLOBAL boundary conditions  
    470          !                                            ! Caution : idta on the global domain: use of jperio, not nperio 
    471452         IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) THEN 
    472453            idta( :    , 1    ) = -1                ;      zdta( :    , 1    ) = -1._wp 
    473             idta( :    ,jpjdta) =  0                ;      zdta( :    ,jpjdta) =  0._wp 
     454            idta( :    ,jpjglo) =  0                ;      zdta( :    ,jpjglo) =  0._wp 
    474455         ELSEIF( jperio == 2 ) THEN 
    475456            idta( :    , 1    ) = idta( : ,  3  )   ;      zdta( :    , 1    ) = zdta( : ,  3  ) 
    476             idta( :    ,jpjdta) = 0                 ;      zdta( :    ,jpjdta) =  0._wp 
     457            idta( :    ,jpjglo) = 0                 ;      zdta( :    ,jpjglo) =  0._wp 
    477458            idta( 1    , :    ) = 0                 ;      zdta( 1    , :    ) =  0._wp 
    478             idta(jpidta, :    ) = 0                 ;      zdta(jpidta, :    ) =  0._wp 
     459            idta(jpiglo, :    ) = 0                 ;      zdta(jpiglo, :    ) =  0._wp 
    479460         ELSE 
    480461            ih = 0                                  ;      zh = 0._wp 
    481462            IF( ln_sco )   ih = jpkm1               ;      IF( ln_sco )   zh = h_oce 
    482463            idta( :    , 1    ) = ih                ;      zdta( :    , 1    ) =  zh 
    483             idta( :    ,jpjdta) = ih                ;      zdta( :    ,jpjdta) =  zh 
     464            idta( :    ,jpjglo) = ih                ;      zdta( :    ,jpjglo) =  zh 
    484465            idta( 1    , :    ) = ih                ;      zdta( 1    , :    ) =  zh 
    485             idta(jpidta, :    ) = ih                ;      zdta(jpidta, :    ) =  zh 
     466            idta(jpiglo, :    ) = ih                ;      zdta(jpiglo, :    ) =  zh 
    486467         ENDIF 
    487468 
    … …  
    497478         risfdep(:,:)=0.e0 
    498479         misfdep(:,:)=1 
    499          ! 
    500          ! (ISF) TODO build ice draft netcdf file for isomip and build the corresponding part of code 
    501          IF( cp_cfg == "isomip" .AND. ln_isfcav ) THEN  
    502             risfdep(:,:)=200.e0  
    503             misfdep(:,:)=1  
    504             ij0 = 1 ; ij1 = 40  
    505             DO jj = mj0(ij0), mj1(ij1)  
    506                risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp  
    507             END DO  
    508             WHERE( bathy(:,:) <= 0._wp )  risfdep(:,:) = 0._wp  
    509          !  
    510          ELSEIF ( cp_cfg == "isomip2" .AND. ln_isfcav ) THEN 
    511          !  
    512             risfdep(:,:)=0.e0 
    513             misfdep(:,:)=1 
    514             ij0 = 1 ; ij1 = 40 
    515             DO jj = mj0(ij0), mj1(ij1) 
    516                risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp 
    517             END DO 
    518             WHERE( bathy(:,:) <= 0._wp )  risfdep(:,:) = 0._wp 
    519          END IF 
    520480         ! 
    521481         DEALLOCATE( idta, zdta ) 
    … …  
    591551               CALL iom_get  ( inum, jpdom_data, 'isf_draft', risfdep ) 
    592552               CALL iom_close( inum ) 
    593                WHERE( bathy(:,:) <= 0._wp )  risfdep(:,:) = 0._wp 
    594  
    595                ! set grounded point to 0  
    596                ! (a treshold could be set here if needed, or set it offline based on the grounded fraction) 
    597                WHERE ( bathy(:,:) <= risfdep(:,:) + rn_isfhmin ) 
    598                   misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 
    599                   mbathy (:,:) = 0 ; bathy  (:,:) = 0._wp 
    600                END WHERE 
    601553            END IF 
    602554            !        
    … …  
    633585      ENDIF 
    634586      ! 
    635     !  IF( nn_closea == 0 )   CALL clo_bat( bathy, mbathy )    !==  NO closed seas or lakes  ==! 
    636       !                        
    637587      IF ( .not. ln_sco ) THEN                                !==  set a minimum depth  ==! 
    638588         IF( rn_hmin < 0._wp ) THEN    ;   ik = - INT( rn_hmin )                                      ! from a nb of level 
    639589         ELSE                          ;   ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 )  ! from a depth 
    640590         ENDIF 
    641          zhmin = gdepw_1d(ik+1)                                                         ! minimum depth = ik+1 w-levels  
     591         zhmin = gdepw_1d(ik+1)                                                        ! minimum depth = ik+1 w-levels  
    642592         WHERE( bathy(:,:) <= 0._wp )   ;   bathy(:,:) = 0._wp                         ! min=0     over the lands 
    643          ELSE WHERE                     ;   bathy(:,:) = MAX(  zhmin , bathy(:,:)  )   ! min=zhmin over the oceans 
     593         ELSE WHERE ( risfdep == 0._wp );   bathy(:,:) = MAX(  zhmin , bathy(:,:)  )   ! min=zhmin over the oceans 
    644594         END WHERE 
    645595         IF(lwp) write(numout,*) 'Minimum ocean depth: ', zhmin, ' minimum number of ocean levels : ', ik 
    646596      ENDIF 
    647597      ! 
    648    !   IF( nn_timing == 1 )  CALL timing_stop('zgr_bat') 
    649       ! 
    650598   END SUBROUTINE zgr_bat 
    651  
    652  
    653    SUBROUTINE zgr_bat_zoom 
    654       !!---------------------------------------------------------------------- 
    655       !!                    ***  ROUTINE zgr_bat_zoom  *** 
    656       !! 
    657       !! ** Purpose : - Close zoom domain boundary if necessary 
    658       !!              - Suppress Med Sea from ORCA R2 and R05 arctic zoom 
    659       !! 
    660       !! ** Method  :  
    661       !! 
    662       !! ** Action  : - update mbathy: level bathymetry (in level index) 
    663       !!---------------------------------------------------------------------- 
    664       INTEGER ::   ii0, ii1, ij0, ij1   ! temporary integers 
    665       !!---------------------------------------------------------------------- 
    666       ! 
    667       IF(lwp) WRITE(numout,*) 
    668       IF(lwp) WRITE(numout,*) '    zgr_bat_zoom : modify the level bathymetry for zoom domain' 
    669       IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~' 
    670       ! 
    671       ! Zoom domain 
    672       ! =========== 
    673       ! 
    674       ! Forced closed boundary if required 
    675       IF( lzoom_s )   mbathy(  : , mj0(jpjzoom):mj1(jpjzoom) )      = 0 
    676       IF( lzoom_w )   mbathy(      mi0(jpizoom):mi1(jpizoom) , :  ) = 0 
    677       IF( lzoom_e )   mbathy(      mi0(jpiglo+jpizoom-1):mi1(jpiglo+jpizoom-1) , :  ) = 0 
    678       IF( lzoom_n )   mbathy(  : , mj0(jpjglo+jpjzoom-1):mj1(jpjglo+jpjzoom-1) )      = 0 
    679       ! 
    680       ! Configuration specific domain modifications 
    681       ! (here, ORCA arctic configuration: suppress Med Sea) 
    682       IF( cp_cfg == "orca" .AND. cp_cfz == "arctic" ) THEN 
    683          SELECT CASE ( jp_cfg ) 
    684          !                                        ! ======================= 
    685          CASE ( 2 )                               !  ORCA_R2 configuration 
    686             !                                     ! ======================= 
    687             IF(lwp) WRITE(numout,*) '                   ORCA R2 arctic zoom: suppress the Med Sea' 
    688             ii0 = 141   ;   ii1 = 162      ! Sea box i,j indices 
    689             ij0 =  98   ;   ij1 = 110 
    690             !                                     ! ======================= 
    691          CASE ( 05 )                              !  ORCA_R05 configuration 
    692             !                                     ! ======================= 
    693             IF(lwp) WRITE(numout,*) '                   ORCA R05 arctic zoom: suppress the Med Sea' 
    694             ii0 = 563   ;   ii1 = 642      ! zero over the Med Sea boxe 
    695             ij0 = 314   ;   ij1 = 370  
    696          END SELECT 
    697          ! 
    698          mbathy( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0   ! zero over the Med Sea boxe 
    699          ! 
    700       ENDIF 
    701       ! 
    702    END SUBROUTINE zgr_bat_zoom 
    703  
    704599 
    705600   SUBROUTINE zgr_bat_ctl 
    … …  
    727622      INTEGER ::   ji, jj, jl                    ! dummy loop indices 
    728623      INTEGER ::   icompt, ibtest, ikmax         ! temporary integers 
    729       REAL(wp), POINTER, DIMENSION(:,:) ::  zbathy 
    730       !!---------------------------------------------------------------------- 
    731       ! 
    732   !    IF( nn_timing == 1 )  CALL timing_start('zgr_bat_ctl') 
    733       ! 
    734       CALL wrk_alloc( jpi, jpj, zbathy ) 
     624      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::  zbathy 
     625      !!---------------------------------------------------------------------- 
     626      ! 
     627      ALLOCATE(zbathy(jpi,jpj)) 
    735628      ! 
    736629      IF(lwp) WRITE(numout,*) 
    … …  
    743636      icompt = 0 
    744637      DO jl = 1, 2 
    745          IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN 
     638         IF( l_Iperio ) THEN 
    746639            mbathy( 1 ,:) = mbathy(jpim1,:)           ! local domain is cyclic east-west 
    747640            mbathy(jpi,:) = mbathy(  2  ,:) 
    748641         ENDIF 
     642         zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     643         CALL lbc_lnk( 'domzgr',zbathy, 'T', 1._wp ) 
     644         mbathy(:,:) = INT( zbathy(:,:) ) 
     645          
    749646         DO jj = 2, jpjm1 
    750647            DO ji = 2, jpim1 
    … …  
    760657         END DO 
    761658      END DO 
    762    !   IF( lk_mpp )   CALL mpp_sum( icompt ) 
     659 
     660      IF( lk_mpp )   CALL mpp_sum( 'domzgr', icompt ) 
    763661      IF( icompt == 0 ) THEN 
    764662         IF(lwp) WRITE(numout,*)'     no isolated ocean grid points' 
    … …  
    766664         IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points suppressed' 
    767665      ENDIF 
    768       IF( lk_mpp ) THEN 
    769          zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    770          CALL lbc_lnk( 'toto',zbathy, 'T', 1._wp ) 
    771          mbathy(:,:) = INT( zbathy(:,:) ) 
    772       ENDIF 
     666 
     667      zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     668      CALL lbc_lnk( 'domzgr',zbathy, 'T', 1._wp ) 
     669      mbathy(:,:) = INT( zbathy(:,:) ) 
     670 
    773671      !                                          ! East-west cyclic boundary conditions 
    774       IF( nperio == 0 ) THEN 
    775          IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west boundary: nperio = ', nperio 
     672      IF( jperio == 0 ) THEN 
     673         IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west boundary: jperio = ', jperio 
    776674         IF( lk_mpp ) THEN 
    777675            IF( nbondi == -1 .OR. nbondi == 2 ) THEN 
    … …  
    790688            ENDIF 
    791689         ENDIF 
    792       ELSEIF( nperio == 1 .OR. nperio == 4 .OR. nperio ==  6 ) THEN 
    793          IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions on mbathy: nperio = ', nperio 
     690      ELSEIF( l_Iperio ) THEN 
     691         IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions on mbathy: jperio = ', jperio 
    794692         mbathy( 1 ,:) = mbathy(jpim1,:) 
    795693         mbathy(jpi,:) = mbathy(  2  ,:) 
    796       ELSEIF( nperio == 2 ) THEN 
    797          IF(lwp) WRITE(numout,*) '   equatorial boundary conditions on mbathy: nperio = ', nperio 
     694      ELSEIF( jperio == 2 ) THEN 
     695         IF(lwp) WRITE(numout,*) '   equatorial boundary conditions on mbathy: jperio = ', jperio 
    798696      ELSE 
    799697         IF(lwp) WRITE(numout,*) '    e r r o r' 
    800          IF(lwp) WRITE(numout,*) '    parameter , nperio = ', nperio 
     698         IF(lwp) WRITE(numout,*) '    parameter , jperio = ', jperio 
    801699         !         STOP 'dom_mba' 
    802700      ENDIF 
     701 
    803702      !  Boundary condition on mbathy 
    804703      IF( .NOT.lk_mpp ) THEN  
    … …  
    806705         !   ... mono- or macro-tasking: T-point, >0, 2D array, no slab 
    807706         zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    808          CALL lbc_lnk( 'toto',zbathy, 'T', 1._wp ) 
     707         CALL lbc_lnk( 'domzgr',zbathy, 'T', 1._wp ) 
    809708         mbathy(:,:) = INT( zbathy(:,:) ) 
    810709      ENDIF 
     710 
    811711      ! Number of ocean level inferior or equal to jpkm1 
    812       ikmax = 0 
    813       DO jj = 1, jpj 
    814          DO ji = 1, jpi 
    815             ikmax = MAX( ikmax, mbathy(ji,jj) ) 
    816          END DO 
    817       END DO 
    818 !!gm  !!! test to do:   ikmax = MAX( mbathy(:,:) )   ??? 
     712      zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     713      ikmax = glob_max( 'domzgr', zbathy(:,:) ) 
     714 
    819715      IF( ikmax > jpkm1 ) THEN 
    820716         IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' >  jpk-1' 
    … …  
    825721      ENDIF 
    826722      ! 
    827       CALL wrk_dealloc( jpi, jpj, zbathy ) 
    828       ! 
    829    !!   IF( nn_timing == 1 )  CALL timing_stop('zgr_bat_ctl') 
     723      DEALLOCATE( zbathy ) 
    830724      ! 
    831725   END SUBROUTINE zgr_bat_ctl 
    … …  
    845739      !!---------------------------------------------------------------------- 
    846740      INTEGER ::   ji, jj   ! dummy loop indices 
    847       REAL(wp), POINTER, DIMENSION(:,:) ::  zmbk 
    848       !!---------------------------------------------------------------------- 
    849       ! 
    850    !   IF( nn_timing == 1 )  CALL timing_start('zgr_bot_level') 
    851       ! 
    852       CALL wrk_alloc( jpi, jpj, zmbk ) 
     741      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::  zmbk 
     742      !!---------------------------------------------------------------------- 
     743      ! 
     744      ALLOCATE( zmbk(jpi,jpj) ) 
    853745      ! 
    854746      IF(lwp) WRITE(numout,*) 
    … …  
    866758      END DO 
    867759      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk  
    868       zmbk(:,:) = REAL( mbku(:,:), wp )   ;   CALL lbc_lnk('toto',zmbk,'U',1.)   ;   mbku  (:,:) = MAX( INT( zmbk(:,:) ), 1 ) 
    869       zmbk(:,:) = REAL( mbkv(:,:), wp )   ;   CALL lbc_lnk('toto',zmbk,'V',1.)   ;   mbkv  (:,:) = MAX( INT( zmbk(:,:) ), 1 ) 
    870       ! 
    871       CALL wrk_dealloc( jpi, jpj, zmbk ) 
    872       ! 
    873    !   IF( nn_timing == 1 )  CALL timing_stop('zgr_bot_level') 
     760      zmbk(:,:) = REAL( mbku(:,:), wp )   ;   CALL lbc_lnk('domzgr',zmbk,'U',1.)   ;   mbku  (:,:) = MAX( INT( zmbk(:,:) ), 1 ) 
     761      zmbk(:,:) = REAL( mbkv(:,:), wp )   ;   CALL lbc_lnk('domzgr',zmbk,'V',1.)   ;   mbkv  (:,:) = MAX( INT( zmbk(:,:) ), 1 ) 
     762      ! 
     763      DEALLOCATE( zmbk ) 
    874764      ! 
    875765   END SUBROUTINE zgr_bot_level 
    … …  
    889779      !!---------------------------------------------------------------------- 
    890780      INTEGER ::   ji, jj   ! dummy loop indices 
    891       REAL(wp), POINTER, DIMENSION(:,:) ::  zmik 
    892       !!---------------------------------------------------------------------- 
    893       ! 
    894    !   IF( nn_timing == 1 )  CALL timing_start('zgr_top_level') 
    895       ! 
    896       CALL wrk_alloc( jpi, jpj, zmik ) 
     781      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::  zmik 
     782      !!---------------------------------------------------------------------- 
     783      ! 
     784      ALLOCATE( zmik(jpi,jpj) ) 
    897785      ! 
    898786      IF(lwp) WRITE(numout,*) 
    … …  
    911799 
    912800      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk  
    913       zmik(:,:) = REAL( miku(:,:), wp )   ;   CALL lbc_lnk('toto',zmik,'U',1.)   ;   miku  (:,:) = MAX( INT( zmik(:,:) ), 1 ) 
    914       zmik(:,:) = REAL( mikv(:,:), wp )   ;   CALL lbc_lnk('toto',zmik,'V',1.)   ;   mikv  (:,:) = MAX( INT( zmik(:,:) ), 1 ) 
    915       zmik(:,:) = REAL( mikf(:,:), wp )   ;   CALL lbc_lnk('toto',zmik,'F',1.)   ;   mikf  (:,:) = MAX( INT( zmik(:,:) ), 1 ) 
    916       ! 
    917       CALL wrk_dealloc( jpi, jpj, zmik ) 
    918       ! 
    919    !   IF( nn_timing == 1 )  CALL timing_stop('zgr_top_level') 
     801      zmik(:,:) = REAL( miku(:,:), wp )   ;   CALL lbc_lnk('domzgr',zmik,'U',1.)   ;   miku  (:,:) = MAX( INT( zmik(:,:) ), 1 ) 
     802      zmik(:,:) = REAL( mikv(:,:), wp )   ;   CALL lbc_lnk('domzgr',zmik,'V',1.)   ;   mikv  (:,:) = MAX( INT( zmik(:,:) ), 1 ) 
     803      zmik(:,:) = REAL( mikf(:,:), wp )   ;   CALL lbc_lnk('domzgr',zmik,'F',1.)   ;   mikf  (:,:) = MAX( INT( zmik(:,:) ), 1 ) 
     804      ! 
     805      DEALLOCATE( zmik ) 
    920806      ! 
    921807   END SUBROUTINE zgr_top_level 
    … …  
    932818      INTEGER  ::   jk 
    933819      !!---------------------------------------------------------------------- 
    934       ! 
    935     !  IF( nn_timing == 1 )  CALL timing_start('zgr_zco') 
    936820      ! 
    937821      DO jk = 1, jpk 
    938822         gdept_0(:,:,jk) = gdept_1d(jk) 
    939823         gdepw_0(:,:,jk) = gdepw_1d(jk) 
    940          gde3w_0(:,:,jk) = gdepw_1d(jk) 
    941824         e3t_0  (:,:,jk) = e3t_1d  (jk) 
    942825         e3u_0  (:,:,jk) = e3t_1d  (jk) 
    … …  
    947830         e3vw_0 (:,:,jk) = e3w_1d  (jk) 
    948831      END DO 
    949       ! 
    950    !   IF( nn_timing == 1 )  CALL timing_stop('zgr_zco') 
    951832      ! 
    952833   END SUBROUTINE zgr_zco 
    … …  
    1004885      REAL(wp) ::   zdiff            ! temporary scalar 
    1005886      REAL(wp) ::   zmax             ! temporary scalar 
    1006       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zprt 
     887      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  zprt 
    1007888      !!--------------------------------------------------------------------- 
    1008889      ! 
    1009    !   IF( nn_timing == 1 )  CALL timing_start('zgr_zps') 
    1010       ! 
    1011       CALL wrk_alloc( jpi,jpj,jpk,   zprt ) 
     890      ALLOCATE( zprt(jpi,jpj,jpk) ) 
    1012891      ! 
    1013892      IF(lwp) WRITE(numout,*) 
    … …  
    1015894      IF(lwp) WRITE(numout,*) '    ~~~~~~~ ' 
    1016895      IF(lwp) WRITE(numout,*) '              mbathy is recomputed : bathy_level file is NOT used' 
     896 
     897      ! compute position of the ice shelf grounding line 
     898      ! set bathy and isfdraft to 0 where grounded 
     899      IF ( ln_isfcav ) CALL zgr_isf_zspace 
    1017900 
    1018901      ! bathymetry in level (from bathy_meter) 
    … …  
    1033916      END DO 
    1034917 
     918      ! Check compatibility between bathy and iceshelf draft 
     919      ! insure at least 2 wet level on the vertical under an ice shelf 
     920      ! compute misfdep and adjust isf draft if needed 
     921      IF ( ln_isfcav ) CALL zgr_isf_kspace 
     922 
    1035923      ! Scale factors and depth at T- and W-points 
    1036924      DO jk = 1, jpk                        ! intitialization to the reference z-coordinate 
    … …  
    1041929      END DO 
    1042930       
    1043       ! Bathy, iceshelf draft, scale factor and depth at T- and W- points in case of isf 
    1044       IF ( ln_isfcav ) CALL zgr_isf 
    1045  
    1046931      ! Scale factors and depth at T- and W-points 
    1047       IF ( .NOT. ln_isfcav ) THEN 
    1048          DO jj = 1, jpj 
    1049             DO ji = 1, jpi 
    1050                ik = mbathy(ji,jj) 
    1051                IF( ik > 0 ) THEN               ! ocean point only 
    1052                   ! max ocean level case 
    1053                   IF( ik == jpkm1 ) THEN 
    1054                      zdepwp = bathy(ji,jj) 
    1055                      ze3tp  = bathy(ji,jj) - gdepw_1d(ik) 
    1056                      ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 
    1057                      e3t_0(ji,jj,ik  ) = ze3tp 
    1058                      e3t_0(ji,jj,ik+1) = ze3tp 
    1059                      e3w_0(ji,jj,ik  ) = ze3wp 
    1060                      e3w_0(ji,jj,ik+1) = ze3tp 
    1061                      gdepw_0(ji,jj,ik+1) = zdepwp 
    1062                      gdept_0(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp 
    1063                      gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 
    1064                      ! 
    1065                   ELSE                         ! standard case 
    1066                      IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN  ;   gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 
    1067                      ELSE                                       ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 
    1068                      ENDIF 
    1069    !gm Bug?  check the gdepw_1d 
    1070                      !       ... on ik 
    1071                      gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )   & 
    1072                         &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   & 
    1073                         &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) )) 
    1074                      e3t_0  (ji,jj,ik) = e3t_1d  (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   &  
    1075                         &                             / ( gdepw_1d(      ik+1) - gdepw_1d(ik) )  
    1076                      e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   & 
    1077                         &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 
    1078                      !       ... on ik+1 
    1079                      e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
    1080                      e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
    1081                      gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 
    1082                   ENDIF 
    1083                ENDIF 
    1084             END DO 
    1085          END DO 
    1086          ! 
    1087          it = 0 
    1088          DO jj = 1, jpj 
    1089             DO ji = 1, jpi 
    1090                ik = mbathy(ji,jj) 
    1091                IF( ik > 0 ) THEN               ! ocean point only 
    1092                   e3tp (ji,jj) = e3t_0(ji,jj,ik) 
    1093                   e3wp (ji,jj) = e3w_0(ji,jj,ik) 
    1094                   ! test 
    1095                   zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik  ) 
    1096                   IF( zdiff <= 0._wp .AND. lwp ) THEN  
    1097                      it = it + 1 
    1098                      WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj 
    1099                      WRITE(numout,*) ' bathy = ', bathy(ji,jj) 
    1100                      WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 
    1101                      WRITE(numout,*) ' e3tp    = ', e3t_0  (ji,jj,ik), ' e3wp    = ', e3w_0  (ji,jj,ik  ) 
    1102                   ENDIF 
    1103                ENDIF 
    1104             END DO 
    1105          END DO 
    1106       END IF 
    1107       ! 
    1108       ! Scale factors and depth at U-, V-, UW and VW-points 
    1109       DO jk = 1, jpk                        ! initialisation to z-scale factors 
    1110          e3u_0 (:,:,jk) = e3t_1d(jk) 
    1111          e3v_0 (:,:,jk) = e3t_1d(jk) 
    1112          e3uw_0(:,:,jk) = e3w_1d(jk) 
    1113          e3vw_0(:,:,jk) = e3w_1d(jk) 
    1114       END DO 
    1115  
    1116       DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors 
    1117          DO jj = 1, jpjm1 
    1118             DO ji = 1, jpim1   ! vector opt. 
    1119                e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) ) 
    1120                e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) ) 
    1121                e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) ) 
    1122                e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) ) 
    1123             END DO 
    1124          END DO 
    1125       END DO 
    1126       IF ( ln_isfcav ) THEN 
    1127       ! (ISF) define e3uw (adapted for 2 cells in the water column) 
    1128          DO jj = 2, jpjm1  
    1129             DO ji = 2, jpim1   ! vector opt.  
    1130                ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj)) 
    1131                ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj)) 
    1132                IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji+1,jj  ,ikb  ) ) & 
    1133                                        &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj  ,ikb-1) ) 
    1134                ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1)) 
    1135                ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1)) 
    1136                IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji  ,jj+1,ikb  ) ) & 
    1137                                        &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji  ,jj+1,ikb-1) ) 
    1138             END DO 
    1139          END DO 
    1140       END IF 
    1141  
    1142       CALL lbc_lnk('toto', e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk('toto', e3uw_0, 'U', 1._wp )   ! lateral boundary conditions 
    1143       CALL lbc_lnk( 'toto',e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk('toto', e3vw_0, 'V', 1._wp ) 
    1144       ! 
    1145  
    1146       DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
    1147          WHERE( e3u_0 (:,:,jk) == 0._wp )   e3u_0 (:,:,jk) = e3t_1d(jk) 
    1148          WHERE( e3v_0 (:,:,jk) == 0._wp )   e3v_0 (:,:,jk) = e3t_1d(jk) 
    1149          WHERE( e3uw_0(:,:,jk) == 0._wp )   e3uw_0(:,:,jk) = e3w_1d(jk) 
    1150          WHERE( e3vw_0(:,:,jk) == 0._wp )   e3vw_0(:,:,jk) = e3w_1d(jk) 
    1151       END DO 
    1152        
    1153       ! Scale factor at F-point 
    1154       DO jk = 1, jpk                        ! initialisation to z-scale factors 
    1155          e3f_0(:,:,jk) = e3t_1d(jk) 
    1156       END DO 
    1157       DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors 
    1158          DO jj = 1, jpjm1 
    1159             DO ji = 1, jpim1   ! vector opt. 
    1160                e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) ) 
    1161             END DO 
    1162          END DO 
    1163       END DO 
    1164       CALL lbc_lnk('toto', e3f_0, 'F', 1._wp )       ! Lateral boundary conditions 
    1165       ! 
    1166       DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
    1167          WHERE( e3f_0(:,:,jk) == 0._wp )   e3f_0(:,:,jk) = e3t_1d(jk) 
    1168       END DO 
    1169 !!gm  bug ? :  must be a do loop with mj0,mj1 
    1170       !  
    1171       e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2 
    1172       e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:)  
    1173       e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:)  
    1174       e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:)  
    1175       e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:)  
    1176  
    1177       ! Control of the sign 
    1178       IF( MINVAL( e3t_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t_0 <= 0' ) 
    1179       IF( MINVAL( e3w_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w_0 <= 0' ) 
    1180       IF( MINVAL( gdept_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdept_0 <  0' ) 
    1181       IF( MINVAL( gdepw_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw_0 <  0' ) 
    1182       
    1183       ! Compute gde3w_0 (vertical sum of e3w) 
    1184       IF ( ln_isfcav ) THEN ! if cavity 
    1185          WHERE( misfdep == 0 )   misfdep = 1 
    1186          DO jj = 1,jpj 
    1187             DO ji = 1,jpi 
    1188                gde3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1) 
    1189                DO jk = 2, misfdep(ji,jj) 
    1190                   gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
    1191                END DO 
    1192                IF( misfdep(ji,jj) >= 2 )   gde3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj)) 
    1193                DO jk = misfdep(ji,jj) + 1, jpk 
    1194                   gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
    1195                END DO 
    1196             END DO 
    1197          END DO 
    1198       ELSE ! no cavity 
    1199          gde3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 
    1200          DO jk = 2, jpk 
    1201             gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk) 
    1202          END DO 
    1203       END IF 
    1204       ! 
    1205       CALL wrk_dealloc( jpi,jpj,jpk,   zprt ) 
    1206       ! 
    1207    !   IF( nn_timing == 1 )  CALL timing_stop('zgr_zps') 
    1208       ! 
    1209    END SUBROUTINE zgr_zps 
    1210  
    1211  
    1212    SUBROUTINE zgr_isf 
    1213       !!---------------------------------------------------------------------- 
    1214       !!                    ***  ROUTINE zgr_isf  *** 
    1215       !!    
    1216       !! ** Purpose :   check the bathymetry in levels 
    1217       !!    
    1218       !! ** Method  :   THe water column have to contained at least 2 cells 
    1219       !!                Bathymetry and isfdraft are modified (dig/close) to respect 
    1220       !!                this criterion. 
    1221       !!    
    1222       !! ** Action  : - test compatibility between isfdraft and bathy  
    1223       !!              - bathy and isfdraft are modified 
    1224       !!---------------------------------------------------------------------- 
    1225       INTEGER  ::   ji, jj, jl, jk       ! dummy loop indices 
    1226       INTEGER  ::   ik, it               ! temporary integers 
    1227       INTEGER  ::   icompt, ibtest       ! (ISF) 
    1228       INTEGER  ::   ibtestim1, ibtestip1 ! (ISF) 
    1229       INTEGER  ::   ibtestjm1, ibtestjp1 ! (ISF) 
    1230       REAL(wp) ::   zdepth           ! Ajusted ocean depth to avoid too small e3t 
    1231       REAL(wp) ::   zmax             ! Maximum and minimum depth 
    1232       REAL(wp) ::   zbathydiff       ! isf temporary scalar 
    1233       REAL(wp) ::   zrisfdepdiff     ! isf temporary scalar 
    1234       REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points 
    1235       REAL(wp) ::   zdepwp           ! Ajusted ocean depth to avoid too small e3t 
    1236       REAL(wp) ::   zdiff            ! temporary scalar 
    1237       REAL(wp), POINTER, DIMENSION(:,:)   ::   zrisfdep, zbathy, zmask   ! 2D workspace (ISH) 
    1238       INTEGER , POINTER, DIMENSION(:,:)   ::   zmbathy, zmisfdep         ! 2D workspace (ISH) 
    1239       !!--------------------------------------------------------------------- 
    1240       ! 
    1241   !!    IF( nn_timing == 1 )   CALL timing_start('zgr_isf') 
    1242       ! 
    1243       CALL wrk_alloc( jpi,jpj,   zbathy, zmask, zrisfdep) 
    1244       CALL wrk_alloc( jpi,jpj,   zmisfdep, zmbathy ) 
    1245  
    1246       ! (ISF) compute misfdep 
    1247       WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) /= 0 ) ;   misfdep(:,:) = 1   ! open water : set misfdep to 1   
    1248       ELSEWHERE                      ;                         misfdep(:,:) = 2   ! iceshelf : initialize misfdep to second level  
    1249       END WHERE   
    1250  
    1251       ! Compute misfdep for ocean points (i.e. first wet level)  
    1252       ! find the first ocean level such that the first level thickness  
    1253       ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where  
    1254       ! e3t_0 is the reference level thickness  
    1255       DO jk = 2, jpkm1  
    1256          zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )  
    1257          WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth )   misfdep(:,:) = jk+1  
    1258       END DO  
    1259       WHERE ( 0._wp < risfdep(:,:) .AND. risfdep(:,:) <= e3t_1d(1) ) 
    1260          risfdep(:,:) = 0. ; misfdep(:,:) = 1 
    1261       END WHERE 
    1262  
    1263       ! remove very shallow ice shelf (less than ~ 10m if 75L) 
    1264       WHERE (risfdep(:,:) <= 10._wp .AND. misfdep(:,:) > 1) 
    1265          misfdep = 0; risfdep = 0.0_wp; 
    1266          mbathy  = 0; bathy   = 0.0_wp; 
    1267       END WHERE 
    1268       WHERE (bathy(:,:) <= 30.0_wp .AND. gphit < -60._wp) 
    1269          misfdep = 0; risfdep = 0.0_wp; 
    1270          mbathy  = 0; bathy   = 0.0_wp; 
    1271       END WHERE 
    1272   
    1273 ! basic check for the compatibility of bathy and risfdep. I think it should be offline because it is not perfect and cannot solved all the situation 
    1274       icompt = 0  
    1275 ! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together 
    1276       DO jl = 1, 10      
    1277          ! check at each iteration if isf is grounded or not (1cm treshold have to be update after first coupling experiments) 
    1278          WHERE (bathy(:,:) <= risfdep(:,:) + rn_isfhmin) 
    1279             misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 
    1280             mbathy (:,:) = 0 ; bathy  (:,:) = 0._wp 
    1281          END WHERE 
    1282          WHERE (mbathy(:,:) <= 0)  
    1283             misfdep(:,:) = 0; risfdep(:,:) = 0._wp  
    1284             mbathy (:,:) = 0; bathy  (:,:) = 0._wp 
    1285          END WHERE 
    1286          IF( lk_mpp ) THEN 
    1287             zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
    1288             CALL lbc_lnk( 'toto',zbathy, 'T', 1. ) 
    1289             misfdep(:,:) = INT( zbathy(:,:) ) 
    1290  
    1291             CALL lbc_lnk( 'toto',risfdep,'T', 1. ) 
    1292             CALL lbc_lnk( 'toto',bathy,  'T', 1. ) 
    1293  
    1294             zbathy(:,:)  = FLOAT( mbathy(:,:) ) 
    1295             CALL lbc_lnk( 'toto',zbathy, 'T', 1. ) 
    1296             mbathy(:,:)  = INT( zbathy(:,:) ) 
    1297          ENDIF 
    1298          IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN  
    1299             misfdep( 1 ,:) = misfdep(jpim1,:)            ! local domain is cyclic east-west  
    1300             misfdep(jpi,:) = misfdep(  2  ,:)  
    1301             mbathy( 1 ,:)  = mbathy(jpim1,:)             ! local domain is cyclic east-west 
    1302             mbathy(jpi,:)  = mbathy(  2  ,:) 
    1303          ENDIF 
    1304  
    1305          ! split last cell if possible (only where water column is 2 cell or less) 
    1306          ! if coupled to ice sheet, we do not modify the bathymetry (can be discuss). 
    1307          IF ( .NOT. ln_iscpl) THEN 
    1308             DO jk = jpkm1, 1, -1 
    1309                zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
    1310                WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy) 
    1311                   mbathy(:,:) = jk 
    1312                   bathy(:,:)  = zmax 
    1313                END WHERE 
    1314             END DO 
    1315          END IF 
    1316   
    1317          ! split top cell if possible (only where water column is 2 cell or less) 
    1318          DO jk = 2, jpkm1 
    1319             zmax = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
    1320             WHERE( gdepw_1d(jk+1) > risfdep(:,:) .AND. risfdep(:,:) >= zmax .AND. misfdep + 1 >= mbathy) 
    1321                misfdep(:,:) = jk 
    1322                risfdep(:,:) = zmax 
    1323             END WHERE 
    1324          END DO 
    1325  
    1326   
    1327  ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition 
    1328          DO jj = 1, jpj 
    1329             DO ji = 1, jpi 
    1330                ! find the minimum change option: 
    1331                ! test bathy 
    1332                IF (risfdep(ji,jj) > 1) THEN 
    1333                   IF ( .NOT. ln_iscpl ) THEN 
    1334                      zbathydiff  =ABS(bathy(ji,jj)   - (gdepw_1d(mbathy (ji,jj)+1) & 
    1335                          &            + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 
    1336                      zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj)  ) & 
    1337                          &            - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 
    1338                      IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) <  misfdep(ji,jj)) THEN 
    1339                         IF (zbathydiff <= zrisfdepdiff) THEN 
    1340                            bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat ) 
    1341                            mbathy(ji,jj)= mbathy(ji,jj) + 1 
    1342                         ELSE 
    1343                            risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ) 
    1344                            misfdep(ji,jj) = misfdep(ji,jj) - 1 
    1345                         END IF 
    1346                      ENDIF 
    1347                   ELSE 
    1348                      IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) <  misfdep(ji,jj)) THEN 
    1349                         risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ) 
    1350                         misfdep(ji,jj) = misfdep(ji,jj) - 1 
    1351                      END IF 
    1352                   END IF 
    1353                END IF 
    1354             END DO 
    1355          END DO 
    1356   
    1357          ! At least 2 levels for water thickness at T, U, and V point. 
    1358          DO jj = 1, jpj 
    1359             DO ji = 1, jpi 
    1360                ! find the minimum change option: 
    1361                ! test bathy 
    1362                IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN 
    1363                   IF ( .NOT. ln_iscpl ) THEN  
    1364                      zbathydiff  =ABS(bathy(ji,jj)   - ( gdepw_1d(mbathy (ji,jj)+1) & 
    1365                          &                             + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 
    1366                      zrisfdepdiff=ABS(risfdep(ji,jj) - ( gdepw_1d(misfdep(ji,jj)  ) &  
    1367                          &                             - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 
    1368                      IF (zbathydiff <= zrisfdepdiff) THEN 
    1369                         mbathy(ji,jj) = mbathy(ji,jj) + 1 
    1370                         bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 
    1371                      ELSE 
    1372                         misfdep(ji,jj)= misfdep(ji,jj) - 1 
    1373                         risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat ) 
    1374                      END IF 
    1375                   ELSE 
    1376                      misfdep(ji,jj)= misfdep(ji,jj) - 1 
    1377                      risfdep(ji,jj)= gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat ) 
    1378                   END IF 
    1379                ENDIF 
    1380             END DO 
    1381          END DO 
    1382   
    1383  ! point V mbathy(ji,jj) == misfdep(ji,jj+1)  
    1384          DO jj = 1, jpjm1 
    1385             DO ji = 1, jpim1 
    1386                IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN 
    1387                   IF ( .NOT. ln_iscpl ) THEN  
    1388                      zbathydiff  =ABS(bathy(ji,jj    ) - ( gdepw_1d(mbathy (ji,jj)+1) & 
    1389                           &                              + MIN( e3zps_min, e3t_1d(mbathy (ji,jj  )+1)*e3zps_rat ))) 
    1390                      zrisfdepdiff=ABS(risfdep(ji,jj+1) - ( gdepw_1d(misfdep(ji,jj+1)) & 
    1391                           &                              - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat ))) 
    1392                      IF (zbathydiff <= zrisfdepdiff) THEN 
    1393                         mbathy(ji,jj) = mbathy(ji,jj) + 1 
    1394                         bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj  )) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj   )+1)*e3zps_rat ) 
    1395                      ELSE 
    1396                         misfdep(ji,jj+1)  = misfdep(ji,jj+1) - 1 
    1397                         risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat ) 
    1398                      END IF 
    1399                   ELSE 
    1400                      misfdep(ji,jj+1)  = misfdep(ji,jj+1) - 1 
    1401                      risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat ) 
    1402                   END IF 
    1403                ENDIF 
    1404             END DO 
    1405          END DO 
    1406   
    1407          IF( lk_mpp ) THEN 
    1408             zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
    1409             CALL lbc_lnk( 'toto',zbathy, 'T', 1. ) 
    1410             misfdep(:,:) = INT( zbathy(:,:) ) 
    1411  
    1412             CALL lbc_lnk( 'toto',risfdep,'T', 1. ) 
    1413             CALL lbc_lnk( 'toto',bathy,  'T', 1. ) 
    1414  
    1415             zbathy(:,:)  = FLOAT( mbathy(:,:) ) 
    1416             CALL lbc_lnk( 'toto',zbathy, 'T', 1. ) 
    1417             mbathy(:,:)  = INT( zbathy(:,:) ) 
    1418          ENDIF 
    1419  ! point V misdep(ji,jj) == mbathy(ji,jj+1)  
    1420          DO jj = 1, jpjm1 
    1421             DO ji = 1, jpim1 
    1422                IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) > 1) THEN 
    1423                   IF ( .NOT. ln_iscpl ) THEN  
    1424                      zbathydiff  =ABS(  bathy(ji,jj+1) - ( gdepw_1d(mbathy (ji,jj+1)+1) & 
    1425                            &                             + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ))) 
    1426                      zrisfdepdiff=ABS(risfdep(ji,jj  ) - ( gdepw_1d(misfdep(ji,jj  )  ) & 
    1427                            &                             - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )-1)*e3zps_rat ))) 
    1428                      IF (zbathydiff <= zrisfdepdiff) THEN 
    1429                         mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1 
    1430                         bathy  (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1)  ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ) 
    1431                      ELSE 
    1432                         misfdep(ji,jj)   = misfdep(ji,jj) - 1 
    1433                         risfdep(ji,jj)   = gdepw_1d(misfdep(ji,jj  )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )  )*e3zps_rat ) 
    1434                      END IF 
    1435                   ELSE 
    1436                      misfdep(ji,jj)   = misfdep(ji,jj) - 1 
    1437                      risfdep(ji,jj)   = gdepw_1d(misfdep(ji,jj  )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )  )*e3zps_rat ) 
    1438                   END IF 
    1439                ENDIF 
    1440             END DO 
    1441          END DO 
    1442   
    1443   
    1444          IF( lk_mpp ) THEN  
    1445             zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
    1446             CALL lbc_lnk( 'toto',zbathy, 'T', 1. ) 
    1447             misfdep(:,:) = INT( zbathy(:,:) ) 
    1448  
    1449             CALL lbc_lnk( 'toto',risfdep,'T', 1. ) 
    1450             CALL lbc_lnk( 'toto',bathy,  'T', 1. ) 
    1451  
    1452             zbathy(:,:)  = FLOAT( mbathy(:,:) ) 
    1453             CALL lbc_lnk( 'toto',zbathy, 'T', 1. ) 
    1454             mbathy(:,:)  = INT( zbathy(:,:) ) 
    1455          ENDIF  
    1456   
    1457  ! point U mbathy(ji,jj) == misfdep(ji,jj+1)  
    1458          DO jj = 1, jpjm1 
    1459             DO ji = 1, jpim1 
    1460                IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN 
    1461                   IF ( .NOT. ln_iscpl ) THEN  
    1462                   zbathydiff  =ABS(  bathy(ji  ,jj) - ( gdepw_1d(mbathy (ji,jj)+1) & 
    1463                        &                              + MIN( e3zps_min, e3t_1d(mbathy (ji  ,jj)+1)*e3zps_rat ))) 
    1464                   zrisfdepdiff=ABS(risfdep(ji+1,jj) - ( gdepw_1d(misfdep(ji+1,jj)) & 
    1465                        &                              - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat ))) 
    1466                   IF (zbathydiff <= zrisfdepdiff) THEN 
    1467                      mbathy(ji,jj) = mbathy(ji,jj) + 1 
    1468                      bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 
    1469                   ELSE 
    1470                      misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1 
    1471                      risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat ) 
    1472                   END IF 
    1473                   ELSE 
    1474                      misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1 
    1475                      risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat ) 
    1476                   ENDIF 
    1477                ENDIF 
    1478             ENDDO 
    1479          ENDDO 
    1480   
    1481          IF( lk_mpp ) THEN  
    1482             zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
    1483             CALL lbc_lnk( 'toto',zbathy, 'T', 1. ) 
    1484             misfdep(:,:) = INT( zbathy(:,:) ) 
    1485  
    1486             CALL lbc_lnk( 'toto',risfdep,'T', 1. ) 
    1487             CALL lbc_lnk( 'toto',bathy,  'T', 1. ) 
    1488  
    1489             zbathy(:,:)  = FLOAT( mbathy(:,:) ) 
    1490             CALL lbc_lnk( 'toto',zbathy, 'T', 1. ) 
    1491             mbathy(:,:)  = INT( zbathy(:,:) ) 
    1492          ENDIF  
    1493   
    1494  ! point U misfdep(ji,jj) == bathy(ji,jj+1)  
    1495          DO jj = 1, jpjm1 
    1496             DO ji = 1, jpim1 
    1497                IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) > 1) THEN 
    1498                   IF ( .NOT. ln_iscpl ) THEN  
    1499                      zbathydiff  =ABS(  bathy(ji+1,jj) - ( gdepw_1d(mbathy (ji+1,jj)+1) & 
    1500                           &                              + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat ))) 
    1501                      zrisfdepdiff=ABS(risfdep(ji  ,jj) - ( gdepw_1d(misfdep(ji  ,jj)  ) & 
    1502                           &                              - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)-1)*e3zps_rat ))) 
    1503                      IF (zbathydiff <= zrisfdepdiff) THEN 
    1504                         mbathy(ji+1,jj)  = mbathy (ji+1,jj) + 1 
    1505                         bathy (ji+1,jj)  = gdepw_1d(mbathy (ji+1,jj)  ) + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat ) 
    1506                      ELSE 
    1507                         misfdep(ji,jj)   = misfdep(ji  ,jj) - 1 
    1508                         risfdep(ji,jj)   = gdepw_1d(misfdep(ji  ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)   )*e3zps_rat ) 
    1509                      END IF 
    1510                   ELSE 
    1511                      misfdep(ji,jj)   = misfdep(ji  ,jj) - 1 
    1512                      risfdep(ji,jj)   = gdepw_1d(misfdep(ji  ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)   )*e3zps_rat ) 
    1513                   ENDIF 
    1514                ENDIF 
    1515             ENDDO 
    1516          ENDDO 
    1517   
    1518          IF( lk_mpp ) THEN 
    1519             zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
    1520             CALL lbc_lnk( 'toto',zbathy, 'T', 1. ) 
    1521             misfdep(:,:) = INT( zbathy(:,:) ) 
    1522  
    1523             CALL lbc_lnk( 'toto',risfdep,'T', 1. ) 
    1524             CALL lbc_lnk( 'toto',bathy,  'T', 1. ) 
    1525  
    1526             zbathy(:,:)  = FLOAT( mbathy(:,:) ) 
    1527             CALL lbc_lnk( 'toto',zbathy, 'T', 1. ) 
    1528             mbathy(:,:)  = INT( zbathy(:,:) ) 
    1529          ENDIF 
    1530       END DO 
    1531       ! end dig bathy/ice shelf to be compatible 
    1532       ! now fill single point in "coastline" of ice shelf, bathy, hole, and test again one cell tickness 
    1533       DO jl = 1,20 
    1534   
    1535  ! remove single point "bay" on isf coast line in the ice shelf draft' 
    1536          DO jk = 2, jpk 
    1537             WHERE (misfdep==0) misfdep=jpk 
    1538             zmask=0._wp 
    1539             WHERE (misfdep <= jk) zmask=1 
    1540             DO jj = 2, jpjm1 
    1541                DO ji = 2, jpim1 
    1542                   IF (misfdep(ji,jj) == jk) THEN 
    1543                      ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 
    1544                      IF (ibtest <= 1) THEN 
    1545                         risfdep(ji,jj)=gdepw_1d(jk+1) ; misfdep(ji,jj)=jk+1 
    1546                         IF (misfdep(ji,jj) > mbathy(ji,jj)) misfdep(ji,jj) = jpk 
    1547                      END IF 
    1548                   END IF 
    1549                END DO 
    1550             END DO 
    1551          END DO 
    1552          WHERE (misfdep==jpk) 
    1553              misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp 
    1554          END WHERE 
    1555          IF( lk_mpp ) THEN 
    1556             zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
    1557             CALL lbc_lnk( 'toto',zbathy, 'T', 1. ) 
    1558             misfdep(:,:) = INT( zbathy(:,:) ) 
    1559  
    1560             CALL lbc_lnk( 'toto',risfdep,'T', 1. ) 
    1561             CALL lbc_lnk('toto', bathy,  'T', 1. ) 
    1562  
    1563             zbathy(:,:)  = FLOAT( mbathy(:,:) ) 
    1564             CALL lbc_lnk( 'toto',zbathy, 'T', 1. ) 
    1565             mbathy(:,:)  = INT( zbathy(:,:) ) 
    1566          ENDIF 
    1567   
    1568  ! remove single point "bay" on bathy coast line beneath an ice shelf' 
    1569          DO jk = jpk,1,-1 
    1570             zmask=0._wp 
    1571             WHERE (mbathy >= jk ) zmask=1 
    1572             DO jj = 2, jpjm1 
    1573                DO ji = 2, jpim1 
    1574                   IF (mbathy(ji,jj) == jk .AND. misfdep(ji,jj) >= 2) THEN 
    1575                      ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 
    1576                      IF (ibtest <= 1) THEN 
    1577                         bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-1 
    1578                         IF (misfdep(ji,jj) > mbathy(ji,jj)) mbathy(ji,jj) = 0 
    1579                      END IF 
    1580                   END IF 
    1581                END DO 
    1582             END DO 
    1583          END DO 
    1584          WHERE (mbathy==0) 
    1585              misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp 
    1586          END WHERE 
    1587          IF( lk_mpp ) THEN 
    1588             zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
    1589             CALL lbc_lnk( 'toto',zbathy, 'T', 1. ) 
    1590             misfdep(:,:) = INT( zbathy(:,:) ) 
    1591  
    1592             CALL lbc_lnk( 'toto',risfdep,'T', 1. ) 
    1593             CALL lbc_lnk( 'toto',bathy,  'T', 1. ) 
    1594  
    1595             zbathy(:,:)  = FLOAT( mbathy(:,:) ) 
    1596             CALL lbc_lnk( 'toto',zbathy, 'T', 1. ) 
    1597             mbathy(:,:)  = INT( zbathy(:,:) ) 
    1598          ENDIF 
    1599   
    1600  ! fill hole in ice shelf 
    1601          zmisfdep = misfdep 
    1602          zrisfdep = risfdep 
    1603          WHERE (zmisfdep <= 1._wp) zmisfdep=jpk 
    1604          DO jj = 2, jpjm1 
    1605             DO ji = 2, jpim1 
    1606                ibtestim1 = zmisfdep(ji-1,jj  ) ; ibtestip1 = zmisfdep(ji+1,jj  ) 
    1607                ibtestjm1 = zmisfdep(ji  ,jj-1) ; ibtestjp1 = zmisfdep(ji  ,jj+1) 
    1608                IF( zmisfdep(ji,jj) >= mbathy(ji-1,jj  ) ) ibtestim1 = jpk 
    1609                IF( zmisfdep(ji,jj) >= mbathy(ji+1,jj  ) ) ibtestip1 = jpk 
    1610                IF( zmisfdep(ji,jj) >= mbathy(ji  ,jj-1) ) ibtestjm1 = jpk 
    1611                IF( zmisfdep(ji,jj) >= mbathy(ji  ,jj+1) ) ibtestjp1 = jpk 
    1612                ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 
    1613                IF( ibtest == jpk .AND. misfdep(ji,jj) >= 2) THEN 
    1614                   mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp 
    1615                END IF 
    1616                IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) >= 2) THEN 
    1617                   misfdep(ji,jj) = ibtest 
    1618                   risfdep(ji,jj) = gdepw_1d(ibtest) 
    1619                ENDIF 
    1620             ENDDO 
    1621          ENDDO 
    1622   
    1623          IF( lk_mpp ) THEN  
    1624             zbathy(:,:)  = FLOAT( misfdep(:,:) )  
    1625             CALL lbc_lnk( 'toto',zbathy,  'T', 1. )  
    1626             misfdep(:,:) = INT( zbathy(:,:) )  
    1627  
    1628             CALL lbc_lnk( 'toto',risfdep, 'T', 1. )  
    1629             CALL lbc_lnk( 'toto',bathy,   'T', 1. ) 
    1630  
    1631             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1632             CALL lbc_lnk( 'toto',zbathy,  'T', 1. ) 
    1633             mbathy(:,:) = INT( zbathy(:,:) ) 
    1634          ENDIF  
    1635  ! 
    1636  !! fill hole in bathymetry 
    1637          zmbathy (:,:)=mbathy (:,:) 
    1638          DO jj = 2, jpjm1 
    1639             DO ji = 2, jpim1 
    1640                ibtestim1 = zmbathy(ji-1,jj  ) ; ibtestip1 = zmbathy(ji+1,jj  ) 
    1641                ibtestjm1 = zmbathy(ji  ,jj-1) ; ibtestjp1 = zmbathy(ji  ,jj+1) 
    1642                IF( zmbathy(ji,jj) <  misfdep(ji-1,jj  ) ) ibtestim1 = 0 
    1643                IF( zmbathy(ji,jj) <  misfdep(ji+1,jj  ) ) ibtestip1 = 0 
    1644                IF( zmbathy(ji,jj) <  misfdep(ji  ,jj-1) ) ibtestjm1 = 0 
    1645                IF( zmbathy(ji,jj) <  misfdep(ji  ,jj+1) ) ibtestjp1 = 0 
    1646                ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 
    1647                IF( ibtest == 0 .AND. misfdep(ji,jj) >= 2) THEN 
    1648                   mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ; 
    1649                END IF 
    1650                IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) >= 2) THEN 
    1651                   mbathy(ji,jj) = ibtest 
    1652                   bathy(ji,jj)  = gdepw_1d(ibtest+1)  
    1653                ENDIF 
    1654             END DO 
    1655          END DO 
    1656          IF( lk_mpp ) THEN  
    1657             zbathy(:,:)  = FLOAT( misfdep(:,:) )  
    1658             CALL lbc_lnk( 'toto',zbathy,  'T', 1. )  
    1659             misfdep(:,:) = INT( zbathy(:,:) )  
    1660  
    1661             CALL lbc_lnk( 'toto',risfdep, 'T', 1. )  
    1662             CALL lbc_lnk( 'toto',bathy,   'T', 1. ) 
    1663  
    1664             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1665             CALL lbc_lnk( 'toto',zbathy,  'T', 1. ) 
    1666             mbathy(:,:) = INT( zbathy(:,:) ) 
    1667          ENDIF  
    1668  ! if not compatible after all check (ie U point water column less than 2 cells), mask U 
    1669          DO jj = 1, jpjm1 
    1670             DO ji = 1, jpim1 
    1671                IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji+1,jj) >= 1) THEN 
    1672                   mbathy(ji,jj)  = mbathy(ji,jj) - 1 ; bathy(ji,jj)   = gdepw_1d(mbathy(ji,jj)+1) ; 
    1673                END IF 
    1674             END DO 
    1675          END DO 
    1676          IF( lk_mpp ) THEN  
    1677             zbathy(:,:)  = FLOAT( misfdep(:,:) )  
    1678             CALL lbc_lnk( 'toto',zbathy,  'T', 1. )  
    1679             misfdep(:,:) = INT( zbathy(:,:) )  
    1680  
    1681             CALL lbc_lnk('toto', risfdep, 'T', 1. )  
    1682             CALL lbc_lnk('toto', bathy,   'T', 1. ) 
    1683  
    1684             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1685             CALL lbc_lnk( 'toto',zbathy,  'T', 1. ) 
    1686             mbathy(:,:) = INT( zbathy(:,:) ) 
    1687          ENDIF  
    1688  ! if not compatible after all check (ie U point water column less than 2 cells), mask U 
    1689          DO jj = 1, jpjm1 
    1690             DO ji = 1, jpim1 
    1691                IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji+1,jj) >= 1) THEN 
    1692                   mbathy(ji+1,jj)  = mbathy(ji+1,jj) - 1;   bathy(ji+1,jj)   = gdepw_1d(mbathy(ji+1,jj)+1) ; 
    1693                END IF 
    1694             END DO 
    1695          END DO 
    1696          IF( lk_mpp ) THEN  
    1697             zbathy(:,:)  = FLOAT( misfdep(:,:) )  
    1698             CALL lbc_lnk('toto', zbathy, 'T', 1. )  
    1699             misfdep(:,:) = INT( zbathy(:,:) )  
    1700  
    1701             CALL lbc_lnk('toto', risfdep,'T', 1. )  
    1702             CALL lbc_lnk( 'toto',bathy,  'T', 1. ) 
    1703  
    1704             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1705             CALL lbc_lnk( 'toto',zbathy, 'T', 1. ) 
    1706             mbathy(:,:) = INT( zbathy(:,:) ) 
    1707          ENDIF  
    1708  ! if not compatible after all check (ie V point water column less than 2 cells), mask V 
    1709          DO jj = 1, jpjm1 
    1710             DO ji = 1, jpi 
    1711                IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji,jj+1) >= 1) THEN 
    1712                   mbathy(ji,jj)  = mbathy(ji,jj) - 1 ; bathy(ji,jj)   = gdepw_1d(mbathy(ji,jj)+1) ; 
    1713                END IF 
    1714             END DO 
    1715          END DO 
    1716          IF( lk_mpp ) THEN  
    1717             zbathy(:,:)  = FLOAT( misfdep(:,:) )  
    1718             CALL lbc_lnk( 'toto',zbathy, 'T', 1. )  
    1719             misfdep(:,:) = INT( zbathy(:,:) )  
    1720  
    1721             CALL lbc_lnk( 'toto',risfdep,'T', 1. )  
    1722             CALL lbc_lnk('toto', bathy,  'T', 1. ) 
    1723  
    1724             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1725             CALL lbc_lnk( 'toto',zbathy, 'T', 1. ) 
    1726             mbathy(:,:) = INT( zbathy(:,:) ) 
    1727          ENDIF  
    1728  ! if not compatible after all check (ie V point water column less than 2 cells), mask V 
    1729          DO jj = 1, jpjm1 
    1730             DO ji = 1, jpi 
    1731                IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji,jj+1) >= 1) THEN 
    1732                   mbathy(ji,jj+1)  = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1) = gdepw_1d(mbathy(ji,jj+1)+1) ; 
    1733                END IF 
    1734             END DO 
    1735          END DO 
    1736          IF( lk_mpp ) THEN  
    1737             zbathy(:,:)  = FLOAT( misfdep(:,:) )  
    1738             CALL lbc_lnk( 'toto',zbathy, 'T', 1. )  
    1739             misfdep(:,:) = INT( zbathy(:,:) )  
    1740  
    1741             CALL lbc_lnk( 'toto',risfdep,'T', 1. )  
    1742             CALL lbc_lnk( 'toto',bathy,  'T', 1. ) 
    1743  
    1744             zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    1745             CALL lbc_lnk( 'toto',zbathy, 'T', 1. ) 
    1746             mbathy(:,:) = INT( zbathy(:,:) ) 
    1747          ENDIF  
    1748  ! if not compatible after all check, mask T 
    1749          DO jj = 1, jpj 
    1750             DO ji = 1, jpi 
    1751                IF (mbathy(ji,jj) <= misfdep(ji,jj)) THEN 
    1752                   misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0._wp ; mbathy(ji,jj)  = 0 ; bathy(ji,jj)   = 0._wp ; 
    1753                END IF 
    1754             END DO 
    1755          END DO 
    1756   
    1757          WHERE (mbathy(:,:) == 1) 
    1758             mbathy = 0; bathy = 0.0_wp ; misfdep = 0 ; risfdep = 0.0_wp 
    1759          END WHERE 
    1760       END DO  
    1761 ! end check compatibility ice shelf/bathy 
    1762       ! remove very shallow ice shelf (less than ~ 10m if 75L) 
    1763       WHERE (risfdep(:,:) <= 10._wp) 
    1764          misfdep = 1; risfdep = 0.0_wp; 
    1765       END WHERE 
    1766  
    1767       IF( icompt == 0 ) THEN  
    1768          IF(lwp) WRITE(numout,*)'     no points with ice shelf too close to bathymetry'  
    1769       ELSE  
    1770          IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry'  
    1771       ENDIF  
    1772  
    1773       ! compute scale factor and depth at T- and W- points 
    1774932      DO jj = 1, jpj 
    1775933         DO ji = 1, jpi 
    … …  
    1793951                  ELSE                                       ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 
    1794952                  ENDIF 
    1795       !            gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 
    1796953!gm Bug?  check the gdepw_1d 
    1797954                  !       ... on ik 
    … …  
    1799956                     &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   & 
    1800957                     &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) )) 
    1801                   e3t_0  (ji,jj,ik  ) = gdepw_0(ji,jj,ik+1) - gdepw_1d(ik  ) 
    1802                   e3w_0  (ji,jj,ik  ) = gdept_0(ji,jj,ik  ) - gdept_1d(ik-1) 
     958                  e3t_0  (ji,jj,ik) = e3t_1d  (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   &  
     959                     &                             / ( gdepw_1d(      ik+1) - gdepw_1d(ik) )  
     960                  e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   & 
     961                     &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 
    1803962                  !       ... on ik+1 
    1804963                  e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
    1805964                  e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
     965                  gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 
    1806966               ENDIF 
    1807967            ENDIF 
    … …  
    1829989      END DO 
    1830990      ! 
    1831       ! (ISF) Definition of e3t, u, v, w for ISF case 
    1832       DO jj = 1, jpj  
    1833          DO ji = 1, jpi  
    1834             ik = misfdep(ji,jj)  
    1835             IF( ik > 1 ) THEN               ! ice shelf point only  
    1836                IF( risfdep(ji,jj) < gdepw_1d(ik) )  risfdep(ji,jj)= gdepw_1d(ik)  
    1837                gdepw_0(ji,jj,ik) = risfdep(ji,jj)  
    1838 !gm Bug?  check the gdepw_0  
    1839             !       ... on ik  
    1840                gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) )   &  
    1841                   &                               * ( gdepw_1d(ik+1) - gdept_1d(ik)      )   &  
    1842                   &                               / ( gdepw_1d(ik+1) - gdepw_1d(ik)      )  
    1843                e3t_0  (ji,jj,ik  ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik)  
    1844                e3w_0  (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 
    1845  
    1846                IF( ik + 1 == mbathy(ji,jj) ) THEN               ! ice shelf point only (2 cell water column)  
    1847                   e3w_0  (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik)  
    1848                ENDIF  
    1849             !       ... on ik / ik-1  
    1850                e3w_0  (ji,jj,ik  ) = e3t_0  (ji,jj,ik) !2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik))  
    1851                gdept_0(ji,jj,ik-1) = gdept_0(ji,jj,ik) - e3w_0(ji,jj,ik) 
    1852                e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 
    1853                e3w_0  (ji,jj,ik-1) = gdept_0(ji,jj,ik-1) - gdept_1d(ik-2) 
    1854                gdepw_0(ji,jj,ik-1) = gdepw_0(ji,jj,ik) - e3t_0(ji,jj,ik-1) 
    1855             ENDIF  
    1856          END DO  
    1857       END DO  
    1858     
    1859       it = 0  
    1860       DO jj = 1, jpj  
    1861          DO ji = 1, jpi  
    1862             ik = misfdep(ji,jj)  
    1863             IF( ik > 1 ) THEN               ! ice shelf point only  
    1864                e3tp (ji,jj) = e3t_0(ji,jj,ik  )  
    1865                e3wp (ji,jj) = e3w_0(ji,jj,ik+1 )  
    1866             ! test  
    1867                zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik  )  
    1868                IF( zdiff <= 0. .AND. lwp ) THEN   
    1869                   it = it + 1  
    1870                   WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj  
    1871                   WRITE(numout,*) ' risfdep = ', risfdep(ji,jj)  
    1872                   WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff  
    1873                   WRITE(numout,*) ' e3tp  = ', e3tp(ji,jj), ' e3wp  = ', e3wp(ji,jj)  
    1874                ENDIF  
    1875             ENDIF  
    1876          END DO  
    1877       END DO  
    1878  
    1879       CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep ) 
    1880       CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy ) 
    1881       ! 
    1882   !    IF( nn_timing == 1 )   CALL timing_stop('zgr_isf') 
    1883       !       
    1884    END SUBROUTINE zgr_isf 
    1885  
     991      ! compute top scale factor if ice shelf 
     992      IF (ln_isfcav) CALL zps_isf 
     993      ! 
     994      ! Scale factors and depth at U-, V-, UW and VW-points 
     995      DO jk = 1, jpk                        ! initialisation to z-scale factors 
     996         e3u_0 (:,:,jk) = e3t_1d(jk) 
     997         e3v_0 (:,:,jk) = e3t_1d(jk) 
     998         e3uw_0(:,:,jk) = e3w_1d(jk) 
     999         e3vw_0(:,:,jk) = e3w_1d(jk) 
     1000      END DO 
     1001 
     1002      DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors 
     1003         DO jj = 1, jpjm1 
     1004            DO ji = 1, jpim1   ! vector opt. 
     1005               e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) ) 
     1006               e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) ) 
     1007               e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) ) 
     1008               e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) ) 
     1009            END DO 
     1010         END DO 
     1011      END DO 
     1012 
     1013      ! update e3uw in case only 2 cells in the water column 
     1014      IF ( ln_isfcav ) CALL zps_isf_e3uv_w 
     1015      ! 
     1016      CALL lbc_lnk('domzgr', e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk('domzgr', e3uw_0, 'U', 1._wp )   ! lateral boundary conditions 
     1017      CALL lbc_lnk('domzgr', e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk('domzgr', e3vw_0, 'V', 1._wp ) 
     1018      ! 
     1019      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
     1020         WHERE( e3u_0 (:,:,jk) == 0._wp )   e3u_0 (:,:,jk) = e3t_1d(jk) 
     1021         WHERE( e3v_0 (:,:,jk) == 0._wp )   e3v_0 (:,:,jk) = e3t_1d(jk) 
     1022         WHERE( e3uw_0(:,:,jk) == 0._wp )   e3uw_0(:,:,jk) = e3w_1d(jk) 
     1023         WHERE( e3vw_0(:,:,jk) == 0._wp )   e3vw_0(:,:,jk) = e3w_1d(jk) 
     1024      END DO 
     1025       
     1026      ! Scale factor at F-point 
     1027      DO jk = 1, jpk                        ! initialisation to z-scale factors 
     1028         e3f_0(:,:,jk) = e3t_1d(jk) 
     1029      END DO 
     1030      DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors 
     1031         DO jj = 1, jpjm1 
     1032            DO ji = 1, jpim1   ! vector opt. 
     1033               e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) ) 
     1034            END DO 
     1035         END DO 
     1036      END DO 
     1037      CALL lbc_lnk('domzgr', e3f_0, 'F', 1._wp )       ! Lateral boundary conditions 
     1038      ! 
     1039      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
     1040         WHERE( e3f_0(:,:,jk) == 0._wp )   e3f_0(:,:,jk) = e3t_1d(jk) 
     1041      END DO 
     1042!!gm  bug ? :  must be a do loop with mj0,mj1 
     1043      !  
     1044      e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2 
     1045      e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:)  
     1046      e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:)  
     1047      e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:)  
     1048      e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:)  
     1049 
     1050      ! Control of the sign 
     1051      IF( MINVAL( e3t_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t_0 <= 0' ) 
     1052      IF( MINVAL( e3w_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w_0 <= 0' ) 
     1053      IF( MINVAL( gdept_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdept_0 <  0' ) 
     1054      IF( MINVAL( gdepw_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw_0 <  0' ) 
     1055      ! 
     1056      ! if in the future gde3w_0 need to be compute, use the function defined in NEMO 
     1057      ! for now gde3w_0 computation is removed as not an output of domcfg 
     1058 
     1059      DEALLOCATE( zprt ) 
     1060      ! 
     1061   END SUBROUTINE zgr_zps 
    18861062 
    18871063   SUBROUTINE zgr_sco 
    … …  
    19351111      REAL(wp) ::   zrfact 
    19361112      ! 
    1937       REAL(wp), POINTER, DIMENSION(:,:  ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2 
    1938       REAL(wp), POINTER, DIMENSION(:,:  ) :: zenv, ztmp, zmsk, zri, zrj, zhbat 
     1113      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2 
     1114      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) :: zenv, ztmp, zmsk, zri, zrj, zhbat 
    19391115      !! 
    19401116      NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, & 
    … …  
    19421118     !!---------------------------------------------------------------------- 
    19431119      ! 
    1944    !!   IF( nn_timing == 1 )  CALL timing_start('zgr_sco') 
    1945       ! 
    1946       CALL wrk_alloc( jpi,jpj,   zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 ) 
     1120      ALLOCATE( zenv(jpi,jpj), ztmp(jpi,jpj), zmsk(jpi,jpj), zri(jpi,jpj), zrj(jpi,jpj), zhbat(jpi,jpj) , ztmpi1(jpi,jpj), ztmpi2(jpi,jpj), ztmpj1(jpi,jpj), ztmpj2(jpi,jpj) ) 
    19471121      ! 
    19481122      REWIND( numnam_ref )              ! Namelist namzgr_sco in reference namelist : Sigma-stretching parameters 
    … …  
    20241198 
    20251199      ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero 
    2026       CALL lbc_lnk( 'toto',zenv, 'T', 1._wp, 'no0' ) 
     1200      CALL lbc_lnk( 'domzgr',zenv, 'T', 1._wp, 'no0' ) 
    20271201      !  
    20281202      ! smooth the bathymetry (if required) 
    … …  
    20881262         END DO 
    20891263         ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero 
    2090          CALL lbc_lnk( 'toto',zenv, 'T', 1._wp, 'no0' ) 
     1264         CALL lbc_lnk( 'domzgr',zenv, 'T', 1._wp, 'no0' ) 
    20911265         !                                                  ! ================ ! 
    20921266      END DO                                                !     End loop     ! 
    … …  
    21321306      ! Apply lateral boundary condition 
    21331307!!gm  ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL 
    2134       zhbat(:,:) = hbatu(:,:)   ;   CALL lbc_lnk('toto', hbatu, 'U', 1._wp ) 
     1308      zhbat(:,:) = hbatu(:,:)   ;   CALL lbc_lnk('domzgr', hbatu, 'U', 1._wp ) 
    21351309      DO jj = 1, jpj 
    21361310         DO ji = 1, jpi 
    … …  
    21421316         END DO 
    21431317      END DO 
    2144       zhbat(:,:) = hbatv(:,:)   ;   CALL lbc_lnk('toto', hbatv, 'V', 1._wp ) 
     1318      zhbat(:,:) = hbatv(:,:)   ;   CALL lbc_lnk('domzgr', hbatv, 'V', 1._wp ) 
    21451319      DO jj = 1, jpj 
    21461320         DO ji = 1, jpi 
    … …  
    21511325         END DO 
    21521326      END DO 
    2153       zhbat(:,:) = hbatf(:,:)   ;   CALL lbc_lnk('toto', hbatf, 'F', 1._wp ) 
     1327      zhbat(:,:) = hbatf(:,:)   ;   CALL lbc_lnk('domzgr', hbatf, 'F', 1._wp ) 
    21541328      DO jj = 1, jpj 
    21551329         DO ji = 1, jpi 
    … …  
    21991373      ENDIF  
    22001374 
    2201       CALL lbc_lnk( 'toto',e3t_0 , 'T', 1._wp ) 
    2202       CALL lbc_lnk( 'toto',e3u_0 , 'U', 1._wp ) 
    2203       CALL lbc_lnk( 'toto',e3v_0 , 'V', 1._wp ) 
    2204       CALL lbc_lnk( 'toto',e3f_0 , 'F', 1._wp ) 
    2205       CALL lbc_lnk( 'toto',e3w_0 , 'W', 1._wp ) 
    2206       CALL lbc_lnk( 'toto',e3uw_0, 'U', 1._wp ) 
    2207       CALL lbc_lnk('toto', e3vw_0, 'V', 1._wp ) 
     1375      CALL lbc_lnk( 'domzgr',e3t_0 , 'T', 1._wp ) 
     1376      CALL lbc_lnk( 'domzgr',e3u_0 , 'U', 1._wp ) 
     1377      CALL lbc_lnk( 'domzgr',e3v_0 , 'V', 1._wp ) 
     1378      CALL lbc_lnk( 'domzgr',e3f_0 , 'F', 1._wp ) 
     1379      CALL lbc_lnk( 'domzgr',e3w_0 , 'W', 1._wp ) 
     1380      CALL lbc_lnk( 'domzgr',e3uw_0, 'U', 1._wp ) 
     1381      CALL lbc_lnk('domzgr', e3vw_0, 'V', 1._wp ) 
    22081382      ! 
    22091383        WHERE( e3t_0 (:,:,:) == 0._wp )   e3t_0 (:,:,:) = 1._wp 
    … …  
    22141388        WHERE( e3uw_0(:,:,:) == 0._wp )   e3uw_0(:,:,:) = 1._wp 
    22151389        WHERE( e3vw_0(:,:,:) == 0._wp )   e3vw_0(:,:,:) = 1._wp 
    2216  
    2217  
    2218 !!gm   I don't like that HERE we are supposed to set the reference coordinate (i.e. _0 arrays) 
    2219 !!gm   and only that !!!!! 
    2220 !!gm   THIS should be removed from here ! 
    2221       gdept_n(:,:,:) = gdept_0(:,:,:) 
    2222       gdepw_n(:,:,:) = gdepw_0(:,:,:) 
    2223       gde3w_n(:,:,:) = gde3w_0(:,:,:) 
    2224       e3t_n  (:,:,:) = e3t_0  (:,:,:) 
    2225       e3u_n  (:,:,:) = e3u_0  (:,:,:) 
    2226       e3v_n  (:,:,:) = e3v_0  (:,:,:) 
    2227       e3f_n  (:,:,:) = e3f_0  (:,:,:) 
    2228       e3w_n  (:,:,:) = e3w_0  (:,:,:) 
    2229       e3uw_n (:,:,:) = e3uw_0 (:,:,:) 
    2230       e3vw_n (:,:,:) = e3vw_0 (:,:,:) 
    2231 !!gm and obviously in the following, use the _0 arrays until the end of this subroutine 
    2232 !! gm end 
    22331390!! 
    22341391      ! HYBRID :  
    … …  
    22361393         DO ji = 1, jpi 
    22371394            DO jk = 1, jpkm1 
    2238                IF( scobot(ji,jj) >= gdept_n(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk ) 
     1395               IF( scobot(ji,jj) >= gdept_0(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk ) 
    22391396            END DO 
    22401397         END DO 
    … …  
    22461403         WRITE(numout,*) ' MIN val mbathy  ', MINVAL( mbathy(:,:)    ), ' MAX ', MAXVAL( mbathy(:,:) ) 
    22471404         WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ),   & 
    2248             &                          ' w ', MINVAL( gdepw_0(:,:,:) ), '3w '  , MINVAL( gde3w_0(:,:,:) ) 
     1405            &                          ' w ', MINVAL( gdepw_0(:,:,:) ) 
    22491406         WRITE(numout,*) ' MIN val e3    t ', MINVAL( e3t_0  (:,:,:) ), ' f '  , MINVAL( e3f_0  (:,:,:) ),   & 
    22501407            &                          ' u ', MINVAL( e3u_0  (:,:,:) ), ' u '  , MINVAL( e3v_0  (:,:,:) ),   & 
    … …  
    22531410 
    22541411         WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ),   & 
    2255             &                          ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w '  , MAXVAL( gde3w_0(:,:,:) ) 
     1412            &                          ' w ', MAXVAL( gdepw_0(:,:,:) ) 
    22561413         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( e3t_0  (:,:,:) ), ' f '  , MAXVAL( e3f_0  (:,:,:) ),   & 
    22571414            &                          ' u ', MAXVAL( e3u_0  (:,:,:) ), ' u '  , MAXVAL( e3v_0  (:,:,:) ),   & 
    … …  
    22981455               DO jk = 1, mbathy(ji,jj) 
    22991456                 ! check coordinate is monotonically increasing 
    2300                  IF (e3w_n(ji,jj,jk) <= 0._wp .OR. e3t_n(ji,jj,jk) <= 0._wp ) THEN 
     1457                 IF (e3w_0(ji,jj,jk) <= 0._wp .OR. e3t_0(ji,jj,jk) <= 0._wp ) THEN 
    23011458                    WRITE(ctmp1,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk 
    23021459                    WRITE(numout,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk 
    2303                     WRITE(numout,*) 'e3w',e3w_n(ji,jj,:) 
    2304                     WRITE(numout,*) 'e3t',e3t_n(ji,jj,:) 
     1460                    WRITE(numout,*) 'e3w',e3w_0(ji,jj,:) 
     1461                    WRITE(numout,*) 'e3t',e3t_0(ji,jj,:) 
    23051462                    CALL ctl_stop( ctmp1 ) 
    23061463                 ENDIF 
    23071464                 ! and check it has never gone negative 
    2308                  IF( gdepw_n(ji,jj,jk) < 0._wp .OR. gdept_n(ji,jj,jk) < 0._wp ) THEN 
     1465                 IF( gdepw_0(ji,jj,jk) < 0._wp .OR. gdept_0(ji,jj,jk) < 0._wp ) THEN 
    23091466                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk 
    23101467                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw   or gdept   =< 0  at point (i,j,k)= ', ji, jj, jk 
    2311                     WRITE(numout,*) 'gdepw',gdepw_n(ji,jj,:) 
    2312                     WRITE(numout,*) 'gdept',gdept_n(ji,jj,:) 
     1468                    WRITE(numout,*) 'gdepw',gdepw_0(ji,jj,:) 
     1469                    WRITE(numout,*) 'gdept',gdept_0(ji,jj,:) 
    23131470                    CALL ctl_stop( ctmp1 ) 
    23141471                 ENDIF 
    23151472                 ! and check it never exceeds the total depth 
    2316                  IF( gdepw_n(ji,jj,jk) > hbatt(ji,jj) ) THEN 
     1473                 IF( gdepw_0(ji,jj,jk) > hbatt(ji,jj) ) THEN 
    23171474                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk 
    23181475                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk 
    2319                     WRITE(numout,*) 'gdepw',gdepw_n(ji,jj,:) 
     1476                    WRITE(numout,*) 'gdepw',gdepw_0(ji,jj,:) 
    23201477                    CALL ctl_stop( ctmp1 ) 
    23211478                 ENDIF 
    … …  
    23241481               DO jk = 1, mbathy(ji,jj)-1 
    23251482                 ! and check it never exceeds the total depth 
    2326                 IF( gdept_n(ji,jj,jk) > hbatt(ji,jj) ) THEN 
     1483                IF( gdept_0(ji,jj,jk) > hbatt(ji,jj) ) THEN 
    23271484                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk 
    23281485                    WRITE(numout,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk 
    2329                     WRITE(numout,*) 'gdept',gdept_n(ji,jj,:) 
     1486                    WRITE(numout,*) 'gdept',gdept_0(ji,jj,:) 
    23301487                    CALL ctl_stop( ctmp1 ) 
    23311488                 ENDIF 
    … …  
    23351492      END DO 
    23361493      ! 
    2337       CALL wrk_dealloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 ) 
    2338       ! 
    2339    !!!   IF( nn_timing == 1 )  CALL timing_stop('zgr_sco') 
     1494      DEALLOCATE( zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 ) 
    23401495      ! 
    23411496   END SUBROUTINE zgr_sco 
    … …  
    23581513      REAL(wp) ::   ztmpu1, ztmpv1, ztmpf1 
    23591514      ! 
    2360       REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 
    2361       REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3            
    2362       !!---------------------------------------------------------------------- 
    2363  
    2364       CALL wrk_alloc( jpi,jpj,jpk,   z_gsigw3, z_gsigt3, z_gsi3w3                                      ) 
    2365       CALL wrk_alloc( jpi,jpj,jpk,   z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 
    2366  
    2367       z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp   ;   z_gsi3w3  = 0._wp 
     1515      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3 
     1516      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3            
     1517      !!---------------------------------------------------------------------- 
     1518 
     1519      ALLOCATE( z_gsigw3 (jpi,jpj,jpk), z_gsigt3 (jpi,jpj,jpk) ) 
     1520      ALLOCATE( z_esigt3 (jpi,jpj,jpk), z_esigw3 (jpi,jpj,jpk), z_esigtu3(jpi,jpj,jpk), z_esigtv3(jpi,jpj,jpk) ) 
     1521      ALLOCATE( z_esigtf3(jpi,jpj,jpk), z_esigwu3(jpi,jpj,jpk), z_esigwv3(jpi,jpj,jpk) ) 
     1522 
     1523      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp 
    23681524      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp  
    23691525      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp 
    … …  
    23921548            z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) ) 
    23931549            ! 
    2394             ! Coefficients for vertical depth as the sum of e3w scale factors 
    2395             z_gsi3w3(ji,jj,1) = 0.5_wp * z_esigw3(ji,jj,1) 
    2396             DO jk = 2, jpk 
    2397                z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk) 
    2398             END DO 
    2399             ! 
    24001550            DO jk = 1, jpk 
    24011551               zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 
    … …  
    24031553               gdept_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft ) 
    24041554               gdepw_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw ) 
    2405                gde3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft ) 
    24061555            END DO 
    24071556           ! 
    … …  
    24481597      END DO 
    24491598      ! 
    2450       CALL wrk_dealloc( jpi,jpj,jpk,   z_gsigw3, z_gsigt3, z_gsi3w3                                      ) 
    2451       CALL wrk_dealloc( jpi,jpj,jpk,  z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 
     1599      DEALLOCATE( z_gsigw3, z_gsigt3                                                        ) 
     1600      DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 
    24521601      ! 
    24531602   END SUBROUTINE s_sh94 
    … …  
    24761625      REAL(wp) ::   ztmpu1, ztmpv1, ztmpf1 
    24771626      ! 
    2478       REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 
    2479       REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3            
    2480       !!---------------------------------------------------------------------- 
    2481       ! 
    2482       CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3                                      ) 
    2483       CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 
    2484  
    2485       z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp   ;   z_gsi3w3  = 0._wp 
     1627      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3 
     1628      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3            
     1629      !!---------------------------------------------------------------------- 
     1630      ! 
     1631      ALLOCATE( z_gsigw3 (jpi,jpj,jpk), z_gsigt3 (jpi,jpj,jpk) ) 
     1632      ALLOCATE( z_esigt3 (jpi,jpj,jpk), z_esigw3 (jpi,jpj,jpk), z_esigtu3(jpi,jpj,jpk), z_esigtv3(jpi,jpj,jpk)) 
     1633      ALLOCATE( z_esigtf3(jpi,jpj,jpk), z_esigwu3(jpi,jpj,jpk), z_esigwv3(jpi,jpj,jpk) ) 
     1634 
     1635      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp 
    24861636      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp  
    24871637      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp 
    … …  
    25351685          z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk)) 
    25361686 
    2537           ! Coefficients for vertical depth as the sum of e3w scale factors 
    2538           z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1) 
    2539           DO jk = 2, jpk 
    2540              z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk) 
    2541           END DO 
    2542  
    25431687          DO jk = 1, jpk 
    25441688             gdept_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk) 
    25451689             gdepw_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk) 
    2546              gde3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk) 
    25471690          END DO 
    25481691 
    … …  
    26081751      ENDDO 
    26091752      ! 
    2610       CALL lbc_lnk('toto',e3t_0 ,'T',1.) ; CALL lbc_lnk('toto',e3u_0 ,'T',1.) 
    2611       CALL lbc_lnk('toto',e3v_0 ,'T',1.) ; CALL lbc_lnk('toto',e3f_0 ,'T',1.) 
    2612       CALL lbc_lnk('toto',e3w_0 ,'T',1.) 
    2613       CALL lbc_lnk('toto',e3uw_0,'T',1.) ; CALL lbc_lnk('toto',e3vw_0,'T',1.) 
    2614       ! 
    2615       CALL wrk_dealloc( jpi,jpj,jpk,   z_gsigw3, z_gsigt3, z_gsi3w3                                      ) 
    2616       CALL wrk_dealloc( jpi,jpj,jpk,  z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 
     1753      CALL lbc_lnk('domzgr',e3t_0 ,'T',1.) ; CALL lbc_lnk('domzgr',e3u_0 ,'T',1.) 
     1754      CALL lbc_lnk('domzgr',e3v_0 ,'T',1.) ; CALL lbc_lnk('domzgr',e3f_0 ,'T',1.) 
     1755      CALL lbc_lnk('domzgr',e3w_0 ,'T',1.) 
     1756      CALL lbc_lnk('domzgr',e3uw_0,'T',1.) ; CALL lbc_lnk('domzgr',e3vw_0,'T',1.) 
     1757      ! 
     1758      DEALLOCATE( z_gsigw3, z_gsigt3                                                        ) 
     1759      DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 
    26171760      ! 
    26181761   END SUBROUTINE s_sf12 
    … …  
    26311774      INTEGER  ::   ji, jj, jk       ! dummy loop argument 
    26321775      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars 
    2633       REAL(wp), POINTER, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w 
    2634       REAL(wp), POINTER, DIMENSION(:) :: z_esigt, z_esigw 
    2635       !!---------------------------------------------------------------------- 
    2636  
    2637       CALL wrk_alloc( jpk,   z_gsigw, z_gsigt, z_gsi3w ) 
    2638       CALL wrk_alloc( jpk,   z_esigt, z_esigw ) 
    2639  
    2640       z_gsigw  = 0._wp   ;   z_gsigt  = 0._wp   ;   z_gsi3w  = 0._wp 
     1776      REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_gsigw, z_gsigt 
     1777      REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_esigt, z_esigw 
     1778      !!---------------------------------------------------------------------- 
     1779 
     1780      ALLOCATE( z_gsigw(jpk), z_gsigt(jpk) ) 
     1781      ALLOCATE( z_esigt(jpk), z_esigw(jpk) ) 
     1782 
     1783      z_gsigw  = 0._wp   ;   z_gsigt  = 0._wp 
    26411784      z_esigt  = 0._wp   ;   z_esigw  = 0._wp  
    26421785 
    … …  
    26571800      z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) ) 
    26581801      ! 
    2659       ! Coefficients for vertical depth as the sum of e3w scale factors 
    2660       z_gsi3w(1) = 0.5_wp * z_esigw(1) 
    2661       DO jk = 2, jpk 
    2662          z_gsi3w(jk) = z_gsi3w(jk-1) + z_esigw(jk) 
    2663       END DO 
    2664 !!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays) 
    26651802      DO jk = 1, jpk 
    26661803         zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 
    … …  
    26681805         gdept_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft ) 
    26691806         gdepw_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw ) 
    2670          gde3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft ) 
    2671       END DO 
    2672 !!gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) 
     1807      END DO 
     1808 
    26731809      DO jj = 1, jpj 
    26741810         DO ji = 1, jpi 
    … …  
    26861822      END DO 
    26871823      ! 
    2688       CALL wrk_dealloc( jpk,   z_gsigw, z_gsigt, z_gsi3w ) 
    2689       CALL wrk_dealloc( jpk,   z_esigt, z_esigw          ) 
     1824      DEALLOCATE( z_gsigw, z_gsigt ) 
     1825      DEALLOCATE( z_esigt, z_esigw ) 
    26901826      ! 
    26911827   END SUBROUTINE s_tanh 

  • utils/tools_UKMO_MERGE_2019/DOMAINcfg/src/ioipsl.f90

    r6951 r12101  
    66! See IOIPSL/IOIPSL_License_CeCILL.txt 
    77! 
    8   USE errioipsl  
     8  USE errioipsl    
     9  USE calendar    
    910  USE stringop 
    10   USE mathelp     
    11   USE getincom 
    12   USE calendar    
    1311  USE fliocom     
    14   USE flincom     
    15   USE histcom     
    16   USE restcom 
     12 
    1713END MODULE ioipsl 

  • utils/tools_UKMO_MERGE_2019/DOMAINcfg/src/iom.F90

    r12100 r12101  
    3535   USE domngb          ! ocean space and time domain 
    3636   USE phycst          ! physical constants 
    37    USE dianam          ! build name of file 
    3837   USE xios 
    3938# endif 
    … …  
    6463   PRIVATE iom_set_rst_context, iom_set_rstw_active, iom_set_rstr_active 
    6564# endif 
    66    PUBLIC iom_set_rstw_var_active, iom_set_rstw_core, iom_set_rst_vars 
     65   PUBLIC iom_set_rstw_var_active, iom_set_rst_vars 
    6766 
    6867   INTERFACE iom_get 
    … …  
    348347#endif 
    349348   END SUBROUTINE iom_set_rstr_active 
    350  
    351    SUBROUTINE iom_set_rstw_core(cdmdl) 
    352       !!--------------------------------------------------------------------- 
    353       !!                   ***  SUBROUTINE  iom_set_rstw_core  *** 
    354       !! 
    355       !! ** Purpose :  set variables which are always in restart file  
    356       !!--------------------------------------------------------------------- 
    357    CHARACTER (len=*), INTENT (IN) :: cdmdl ! model OPA or SAS 
    358    CHARACTER(LEN=256)             :: clinfo    ! info character 
    359 #if defined key_iomput 
    360    IF(cdmdl == "OPA") THEN 
    361 !from restart.F90 
    362    CALL iom_set_rstw_var_active("rdt") 
    363    IF ( .NOT. ln_diurnal_only ) THEN 
    364         CALL iom_set_rstw_var_active('ub'  ) 
    365         CALL iom_set_rstw_var_active('vb'  ) 
    366         CALL iom_set_rstw_var_active('tb'  ) 
    367         CALL iom_set_rstw_var_active('sb'  ) 
    368         CALL iom_set_rstw_var_active('sshb') 
    369         ! 
    370         CALL iom_set_rstw_var_active('un'  ) 
    371         CALL iom_set_rstw_var_active('vn'  ) 
    372         CALL iom_set_rstw_var_active('tn'  ) 
    373         CALL iom_set_rstw_var_active('sn'  ) 
    374         CALL iom_set_rstw_var_active('sshn') 
    375         CALL iom_set_rstw_var_active('rhop') 
    376      ! extra variable needed for the ice sheet coupling 
    377         IF ( ln_iscpl ) THEN 
    378              CALL iom_set_rstw_var_active('tmask') 
    379              CALL iom_set_rstw_var_active('umask') 
    380              CALL iom_set_rstw_var_active('vmask') 
    381              CALL iom_set_rstw_var_active('smask') 
    382              CALL iom_set_rstw_var_active('e3t_n') 
    383              CALL iom_set_rstw_var_active('e3u_n') 
    384              CALL iom_set_rstw_var_active('e3v_n') 
    385              CALL iom_set_rstw_var_active('gdepw_n') 
    386         END IF 
    387       ENDIF 
    388       IF(ln_diurnal) CALL iom_set_rstw_var_active('Dsst') 
    389 !from trasbc.F90 
    390          CALL iom_set_rstw_var_active('sbc_hc_b') 
    391          CALL iom_set_rstw_var_active('sbc_sc_b') 
    392    ENDIF 
    393 #else 
    394         clinfo = 'iom_set_rstw_core: key_iomput is needed to use XIOS restart read/write functionality' 
    395         CALL ctl_stop('STOP', TRIM(clinfo)) 
    396 #endif 
    397    END SUBROUTINE iom_set_rstw_core 
    398349 
    399350   SUBROUTINE iom_set_rst_vars(fields) 
    … …  
    662613      ! start halo size for x,y dimensions 
    663614      ! end halo size for x,y dimensions 
     615      ! 
     616      INTEGER ::   nldi_save, nlei_save    !:patch before we remove periodicity and close boundaries in output files 
     617      INTEGER ::   nldj_save, nlej_save    !: 
     618      ! 
    664619      !--------------------------------------------------------------------- 
    665620      ! Initializations and control 
    … …  
    668623      clinfo = '                    iom_open ~~~  ' 
    669624      istop = nstop 
     625 
     626      ! use patch to force the writing off periodicity and close boundaries 
     627      ! without this, issue in some model decomposition 
     628      ! seb: patch before we remove periodicity and close boundaries in output files 
     629      nldi_save = nldi   ;   nlei_save = nlei 
     630      nldj_save = nldj   ;   nlej_save = nlej 
     631      IF( nimpp           ==      1 ) nldi = 1 
     632      IF( nimpp + jpi - 1 == jpiglo ) nlei = jpi 
     633      IF( njmpp           ==      1 ) nldj = 1 
     634      IF( njmpp + jpj - 1 == jpjglo ) nlej = jpj 
     635 
    670636      ! if iom_open is called for the first time: initialize iom_file(:)%nfid to 0 
    671637      ! (could be done when defining iom_file in f95 but not in f90) 
    … …  
    694660      ! do we read the overlap  
    695661      ! ugly patch SM+JMM+RB to overwrite global definition in some cases 
    696       llnoov = (jpni * jpnj ) == jpnij .AND. .NOT. lk_agrif 
     662      !llnoov = (jpni * jpnj ) == jpnij .AND. .NOT. lk_agrif 
     663      ! for domain_cfg, force to read the full domain 
     664      llnoov = .FALSE. 
    697665      ! create the file name by added, if needed, TRIM(Agrif_CFixed()) and TRIM(clsuffix) 
    698666      ! ============= 
    … …  
    792760         CALL iom_nf90_open( clname, kiomid, llwrt, llok, idompar, kdlev = kdlev ) 
    793761      ENDIF 
     762 
     763      nldi = nldi_save   ;   nlei = nlei_save 
     764      nldj = nldj_save   ;   nlej = nlej_save 
    794765      ! 
    795766   END SUBROUTINE iom_open 
    … …  
    10831054         ! do we read the overlap  
    10841055         ! ugly patch SM+JMM+RB to overwrite global definition in some cases 
    1085          llnoov = (jpni * jpnj ) == jpnij .AND. .NOT. lk_agrif  
     1056         !  
     1057         !llnoov = (jpni * jpnj ) == jpnij .AND. .NOT. lk_agrif  
     1058         ! for domain_cfg tools force to read the full domain 
     1059         llnoov = .FALSE. 
    10861060         ! check kcount and kstart optionals parameters... 
    10871061         IF( PRESENT(kcount) .AND. (.NOT. PRESENT(kstart)) ) CALL ctl_stop(trim(clinfo), 'kcount present needs kstart present') 
    … …  
    12641238               ENDIF 
    12651239            ENDIF 
    1266        
     1240            WRITE(numout,*) 'istart icnt',istart, ' ', icnt  
     1241            WRITE(numout,*) ' idx i1, i2, j1, j2 ',ix1, ix2, iy1, iy2 
    12671242            CALL iom_nf90_get( kiomid, idvar, inbdim, istart, icnt, ix1, ix2, iy1, iy2, pv_r1d, pv_r2d, pv_r3d ) 
    12681243 

  • utils/tools_UKMO_MERGE_2019/DOMAINcfg/src/iom_nf90.F90

    r12100 r12101  
    529529      INTEGER               :: idlv                 ! local variable 
    530530      INTEGER               :: idim3                ! id of the third dimension 
     531      ! 
     532      INTEGER ::   nldi_save, nlei_save    !:patch before we remove periodicity and close boundaries in output files 
     533      INTEGER ::   nldj_save, nlej_save    !: 
    531534      !--------------------------------------------------------------------- 
    532535      ! 
    533536      clinfo = '          iom_nf90_rp0123d, file: '//TRIM(iom_file(kiomid)%name)//', var: '//TRIM(cdvar) 
    534537      if90id = iom_file(kiomid)%nfid 
     538      ! 
     539      ! use patch to force the writing off periodicity and close boundaries 
     540      ! without this, issue in some model decomposition 
     541      ! seb: patch before we remove periodicity and close boundaries in output files 
     542      nldi_save = nldi   ;   nlei_save = nlei 
     543      nldj_save = nldj   ;   nlej_save = nlej 
     544      IF( nimpp           ==      1 ) nldi = 1 
     545      IF( nimpp + jpi - 1 == jpiglo ) nlei = jpi 
     546      IF( njmpp           ==      1 ) nldj = 1 
     547      IF( njmpp + jpj - 1 == jpjglo ) nlej = jpj 
    535548      ! 
    536549      ! define dimension variables if it is not already done 
    … …  
    705718         IF(lwp) WRITE(numout,*) TRIM(clinfo)//' written ok' 
    706719      ENDIF 
     720      ! 
     721      nldi = nldi_save   ;   nlei = nlei_save 
     722      nldj = nldj_save   ;   nlej = nlej_save 
    707723      !      
    708724   END SUBROUTINE iom_nf90_rp0123d 

  • utils/tools_UKMO_MERGE_2019/DOMAINcfg/src/lbclnk.F90

    r12100 r12101  
    7171   !!   lbc_bdy_lnk   : set the lateral BDY boundary condition 
    7272   !!---------------------------------------------------------------------- 
    73    USE oce            ! ocean dynamics and tracers    
    7473   USE dom_oce        ! ocean space and time domain  
    7574   USE in_out_manager ! I/O manager 

  • utils/tools_UKMO_MERGE_2019/DOMAINcfg/src/mpp_loc_generic.h90

    r12100 r12101  
    1717#      define MPI_OPERATION mpi_maxloc 
    1818#      define LOC_OPERATION MAXLOC 
     19#      define ERRVAL -HUGE 
    1920#   endif 
    2021#   if defined OPERATION_MINLOC 
    2122#      define MPI_OPERATION mpi_minloc 
    2223#      define LOC_OPERATION MINLOC 
     24#      define ERRVAL HUGE 
    2325#   endif 
    2426 
    … …  
    4244      ! 
    4345      idim = SIZE(kindex) 
    44       ALLOCATE ( ilocs(idim) ) 
    4546      ! 
    46       ilocs = LOC_OPERATION( ARRAY_IN(:,:,:) , mask= MASK_IN(:,:,:) == 1._wp ) 
    47       zmin  = ARRAY_IN(ilocs(1),ilocs(2),ilocs(3)) 
    48       ! 
    49       kindex(1) = ilocs(1) + nimpp - 1 
     47      IF ( ALL(MASK_IN(:,:,:) /= 1._wp) ) THEN 
     48         ! special case for land processors 
     49         zmin = ERRVAL(zmin) 
     50         index0 = 0 
     51      ELSE 
     52         ALLOCATE ( ilocs(idim) ) 
     53         ! 
     54         ilocs = LOC_OPERATION( ARRAY_IN(:,:,:) , mask= MASK_IN(:,:,:) == 1._wp ) 
     55         zmin  = ARRAY_IN(ilocs(1),ilocs(2),ilocs(3)) 
     56         ! 
     57         kindex(1) = mig( ilocs(1) ) 
    5058#  if defined DIM_2d || defined DIM_3d    /* avoid warning when kindex has 1 element */ 
    51       kindex(2) = ilocs(2) + njmpp - 1 
     59         kindex(2) = mjg( ilocs(2) ) 
    5260#  endif 
    5361#  if defined DIM_3d                      /* avoid warning when kindex has 2 elements */ 
    54       kindex(3) = ilocs(3) 
     62         kindex(3) = ilocs(3) 
    5563#  endif 
    56       !  
    57       DEALLOCATE (ilocs) 
    58       ! 
    59       index0 = kindex(1)-1   ! 1d index starting at 0 
     64         !  
     65         DEALLOCATE (ilocs) 
     66         ! 
     67         index0 = kindex(1)-1   ! 1d index starting at 0 
    6068#  if defined DIM_2d || defined DIM_3d   /* avoid warning when kindex has 1 element */ 
    61       index0 = index0 + jpiglo * (kindex(2)-1) 
     69         index0 = index0 + jpiglo * (kindex(2)-1) 
    6270#  endif 
    6371#  if defined DIM_3d                     /* avoid warning when kindex has 2 elements */ 
    64       index0 = index0 + jpiglo * jpjglo * (kindex(3)-1) 
     72         index0 = index0 + jpiglo * jpjglo * (kindex(3)-1) 
    6573#  endif 
     74      END IF 
    6675      zain(1,:) = zmin 
    6776      zain(2,:) = REAL(index0, wp) 
    … …  
    98107#undef LOC_OPERATION 
    99108#undef INDEX_TYPE 
     109#undef ERRVAL 

  • utils/tools_UKMO_MERGE_2019/DOMAINcfg/src/nemogcm.F90

    r12100 r12101  
    4444   !!   factorise     : calculate the factors of the no. of MPI processes 
    4545   !!---------------------------------------------------------------------- 
    46    USE step_oce       ! module used in the ocean time stepping module (step.F90) 
     46   USE dom_oce        ! ocean space and time domain variables 
     47   USE in_out_manager ! I/O manager 
     48   USE iom            ! 
    4749   USE domcfg         ! domain configuration               (dom_cfg routine) 
    4850   USE mppini         ! shared/distributed memory setting (mpp_init routine) 
    … …  
    141143      INTEGER  ::   ios, ilocal_comm   ! local integers 
    142144      CHARACTER(len=120), DIMENSION(60) ::   cltxt, cltxt2, clnam 
    143        ! 
    144       NAMELIST/namctl/ ln_ctl   , sn_cfctl, nn_print,ln_timing 
     145      !! 
     146      NAMELIST/namctl/ ln_ctl   , sn_cfctl, nn_print, nn_ictls, nn_ictle,   & 
     147         &             nn_isplt , nn_jsplt, nn_jctls, nn_jctle,             & 
     148         &             ln_timing, ln_diacfl 
    145149      NAMELIST/namcfg/ ln_e3_dep,                                & 
    146150         &             cp_cfg, cp_cfz, jp_cfg, jpidta, jpjdta, jpkdta, jpiglo, jpjglo, & 
    147          &             jpizoom, jpjzoom, jperio, ln_use_jattr 
     151         &             jperio, ln_use_jattr, ln_domclo 
    148152      !!---------------------------------------------------------------------- 
    149153      ! 
    … …  
    154158      CALL ctl_opn( numnam_cfg, 'namelist_cfg', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE. ) 
    155159      ! 
    156       REWIND( numnam_ref )              ! Namelist namctl in reference namelist : Control prints & Benchmark 
     160      REWIND( numnam_ref )              ! Namelist namctl in reference namelist 
    157161      READ  ( numnam_ref, namctl, IOSTAT = ios, ERR = 901 ) 
    158 901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namctl in reference namelist', .TRUE. ) 
    159  
    160       REWIND( numnam_cfg )              ! Namelist namctl in confguration namelist : Control prints & Benchmark 
     162901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namctl in reference namelist', .TRUE. ) 
     163      REWIND( numnam_cfg )              ! Namelist namctl in confguration namelist 
    161164      READ  ( numnam_cfg, namctl, IOSTAT = ios, ERR = 902 ) 
    162       902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namctl in configuration namelist', .TRUE. ) 
    163  
    164       ! 
    165       REWIND( numnam_ref )              ! Namelist namcfg in reference namelist : Control prints & Benchmark 
     165902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namctl in configuration namelist', .TRUE. ) 
     166      ! 
     167      REWIND( numnam_ref )              ! Namelist namcfg in reference namelist 
    166168      READ  ( numnam_ref, namcfg, IOSTAT = ios, ERR = 903 ) 
    167 903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcfg in reference namelist', .TRUE. ) 
    168  
    169       REWIND( numnam_cfg )              ! Namelist namcfg in confguration namelist : Control prints & Benchmark 
     169903   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namcfg in reference namelist', .TRUE. ) 
     170      REWIND( numnam_cfg )              ! Namelist namcfg in confguration namelist 
    170171      READ  ( numnam_cfg, namcfg, IOSTAT = ios, ERR = 904 ) 
    171 904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcfg in configuration namelist', .TRUE. )    
     172904   IF( ios >  0 )  CALL ctl_nam ( ios , 'namcfg in configuration namelist', .TRUE. )    
    172173 
    173174! Force values for AGRIF zoom (cf. agrif_user.F90) 
    … …  
    265266                            CALL     dom_cfg    ! Domain configuration 
    266267                            CALL     dom_init   ! Domain 
    267       IF( ln_ctl        )   CALL prt_ctl_init   ! Print control 
    268268      ! 
    269269   END SUBROUTINE nemo_init 
    … …  
    315315      jsplt     = nn_jsplt 
    316316 
    317      !  IF( .NOT.ln_read_cfg )   ln_closea = .false.   ! dealing possible only with a domcfg file 
    318317      ! 
    319318      !                             ! Parameter control 
    … …  
    409408      !!---------------------------------------------------------------------- 
    410409      ! 
    411       ierr =        oce_alloc       ()          ! ocean 
    412       ierr = ierr + dom_oce_alloc   ()          ! ocean domain 
     410      ierr = dom_oce_alloc   ()          ! ocean domain 
    413411      ! 
    414412      CALL mpp_sum( 'nemogcm', ierr ) 

  • utils/tools_UKMO_MERGE_2019/DOMAINcfg/src/par_oce.f90

    r9598 r12101  
    1313   PUBLIC 
    1414 
     15   CHARACTER(lc) ::   cp_cfg           !: name of the configuration 
     16   CHARACTER(lc) ::   cp_cfz           !: name of the zoom of configuration 
     17   INTEGER       ::   jp_cfg           !: resolution of the configuration 
     18 
     19   ! data size                                       !!! * size of all input files * 
     20   INTEGER       ::   jpidta           !: 1st lateral dimension ( >= jpi ) 
     21   INTEGER       ::   jpjdta           !: 2nd    "         "    ( >= jpj ) 
     22   INTEGER       ::   jpkdta           !: number of levels      ( >= jpk ) 
     23   LOGICAL       ::   ln_e3_dep        ! e3. definition flag 
     24   REAL(wp)      ::   pp_not_used       = 999999._wp   !: vertical grid parameter 
     25   REAL(wp)      ::   pp_to_be_computed = 999999._wp   !:    -      -       - 
     26   !!---------------------------------------------------------------------- 
     27   !!                   namcfg namelist parameters 
     28   !!---------------------------------------------------------------------- 
     29   LOGICAL       ::   ln_read_cfg      !: (=T) read the domain configuration file or (=F) not 
     30   CHARACTER(lc) ::      cn_domcfg        !: filename the configuration file to be read 
     31   LOGICAL       ::   ln_write_cfg     !: (=T) create the domain configuration file 
     32   CHARACTER(lc) ::      cn_domcfg_out    !: filename the configuration file to be read 
     33   ! 
     34   LOGICAL       ::   ln_use_jattr     !: input file read offset 
     35   !                                   !  Use file global attribute: open_ocean_jstart to determine start j-row  
     36   !                                   !  when reading input from those netcdf files that have the  
     37   !                                   !  attribute defined. This is designed to enable input files associated  
     38   !                                   !  with the extended grids used in the under ice shelf configurations to  
     39   !                                   !  be used without redundant rows when the ice shelves are not in use. 
     40   !  
     41 
     42   !!--------------------------------------------------------------------- 
     43   !! Domain Matrix size  
     44   !!--------------------------------------------------------------------- 
     45   ! configuration name & resolution   (required only in ORCA family case) 
     46   CHARACTER(lc) ::   cn_cfg           !: name of the configuration 
     47   INTEGER       ::   nn_cfg           !: resolution of the configuration  
     48 
     49   ! global domain size               !!! * total computational domain * 
     50   INTEGER       ::   jpiglo           !: 1st dimension of global domain --> i-direction 
     51   INTEGER       ::   jpjglo           !: 2nd    -                  -    --> j-direction 
     52   INTEGER       ::   jpkglo           !: 3nd    -                  -    --> k levels 
     53 
     54   ! global domain size for AGRIF     !!! * total AGRIF computational domain * 
     55   INTEGER, PUBLIC            ::   nbug_in_agrif_conv_do_not_remove_or_modify = 1 - 1 
     56   INTEGER, PUBLIC, PARAMETER ::   nbghostcells = 3                             !: number of ghost cells 
     57   INTEGER, PUBLIC            ::   nbcellsx   ! = jpiglo - 2 - 2*nbghostcells   !: number of cells in i-direction 
     58   INTEGER, PUBLIC            ::   nbcellsy   ! = jpjglo - 2 - 2*nbghostcells   !: number of cells in j-direction 
     59 
     60   ! local domain size                !!! * local computational domain * 
     61   INTEGER, PUBLIC ::   jpi   !                                                    !: first  dimension 
     62   INTEGER, PUBLIC ::   jpj   !                                                    !: second dimension 
     63   INTEGER, PUBLIC ::   jpk   ! = jpkglo                                           !: third  dimension 
     64   INTEGER, PUBLIC ::   jpim1 ! = jpi-1                                            !: inner domain indices 
     65   INTEGER, PUBLIC ::   jpjm1 ! = jpj-1                                            !:   -     -      - 
     66   INTEGER, PUBLIC ::   jpkm1 ! = jpk-1                                            !:   -     -      - 
     67   INTEGER, PUBLIC ::   jpij  ! = jpi*jpj                                          !:  jpi x jpj 
     68   INTEGER, PUBLIC ::   jpimax! = ( jpiglo-2*nn_hls + (jpni-1) ) / jpni + 2*nn_hls !: maximum jpi 
     69   INTEGER, PUBLIC ::   jpjmax! = ( jpjglo-2*nn_hls + (jpnj-1) ) / jpnj + 2*nn_hls !: maximum jpj 
     70 
     71   !!--------------------------------------------------------------------- 
     72   !! Active tracer parameters 
     73   !!--------------------------------------------------------------------- 
     74   INTEGER, PUBLIC, PARAMETER ::   jpts   = 2    !: Number of active tracers (=2, i.e. T & S ) 
     75   INTEGER, PUBLIC, PARAMETER ::   jp_tem = 1    !: indice for temperature 
     76   INTEGER, PUBLIC, PARAMETER ::   jp_sal = 2    !: indice for salinity 
     77 
    1578   !!---------------------------------------------------------------------- 
    1679   !!   Domain decomposition 
    … …  
    2285   INTEGER, PUBLIC, PARAMETER ::   jpr2di = 0   !: number of columns for extra outer halo  
    2386   INTEGER, PUBLIC, PARAMETER ::   jpr2dj = 0   !: number of rows    for extra outer halo  
    24    INTEGER, PUBLIC, PARAMETER ::   jpreci = 1   !: number of columns for overlap  
    25    INTEGER, PUBLIC, PARAMETER ::   jprecj = 1   !: number of rows    for overlap  
    26  
    27    !!---------------------------------------------------------------------- 
    28    !!                   namcfg namelist parameters 
    29    !!---------------------------------------------------------------------- 
    30    ! 
    31    LOGICAL       ::   ln_e3_dep        ! e3. definition flag 
    32    ! 
    33    CHARACTER(lc) ::   cp_cfg           !: name of the configuration 
    34    CHARACTER(lc) ::   cp_cfz           !: name of the zoom of configuration 
    35    INTEGER       ::   jp_cfg           !: resolution of the configuration 
    36  
    37    ! data size                                       !!! * size of all input files * 
    38    INTEGER       ::   jpidta           !: 1st lateral dimension ( >= jpi ) 
    39    INTEGER       ::   jpjdta           !: 2nd    "         "    ( >= jpj ) 
    40    INTEGER       ::   jpkdta           !: number of levels      ( >= jpk ) 
    41  
    42    ! global or zoom domain size                      !!! * computational domain * 
    43    INTEGER       ::   jpiglo           !: 1st dimension of global domain --> i 
    44    INTEGER       ::   jpjglo           !: 2nd    -                  -    --> j 
    45  
    46    ! zoom starting position  
    47    INTEGER       ::   jpizoom          !: left bottom (i,j) indices of the zoom 
    48    INTEGER       ::   jpjzoom          !: in data domain indices 
    49  
    50    ! Domain characteristics 
    51    INTEGER       ::   jperio           !: lateral cond. type (between 0 and 6) 
    52    !                                       !  = 0 closed                 ;   = 1 cyclic East-West 
    53    !                                       !  = 2 equatorial symmetric   ;   = 3 North fold T-point pivot 
    54    !                                       !  = 4 cyclic East-West AND North fold T-point pivot 
    55    !                                       !  = 5 North fold F-point pivot 
    56    !                                       !  = 6 cyclic East-West AND North fold F-point pivot 
    57  
    58    ! Input file read offset 
    59    LOGICAL       ::   ln_use_jattr     !: Use file global attribute: open_ocean_jstart to determine start j-row  
    60                                            ! when reading input from those netcdf files that have the  
    61                                            ! attribute defined. This is designed to enable input files associated  
    62                                            ! with the extended grids used in the under ice shelf configurations to  
    63                                            ! be used without redundant rows when the ice shelves are not in use. 
    64  
    65    !!  Values set to pp_not_used indicates that this parameter is not used in THIS config. 
    66    !!  Values set to pp_to_be_computed  indicates that variables will be computed in domzgr 
    67    REAL(wp)      ::   pp_not_used       = 999999._wp   !: vertical grid parameter 
    68    REAL(wp)      ::   pp_to_be_computed = 999999._wp   !:    -      -       - 
    69  
    70  
    71  
    72  
    73    !!--------------------------------------------------------------------- 
    74    !! Active tracer parameters 
    75    !!--------------------------------------------------------------------- 
    76    INTEGER, PUBLIC, PARAMETER ::   jpts   = 2    !: Number of active tracers (=2, i.e. T & S ) 
    77    INTEGER, PUBLIC, PARAMETER ::   jp_tem = 1    !: indice for temperature 
    78    INTEGER, PUBLIC, PARAMETER ::   jp_sal = 2    !: indice for salinity 
    79  
    80    !!--------------------------------------------------------------------- 
    81    !! Domain Matrix size  (if AGRIF, they are not all parameters) 
    82    !!--------------------------------------------------------------------- 
    83  
    84  
    85  
    86  
    87  
    88  
    89    INTEGER, PUBLIC  ::   jpi   ! = ( jpiglo-2*jpreci + (jpni-1) ) / jpni + 2*jpreci   !: first  dimension 
    90    INTEGER, PUBLIC  ::   jpj   ! = ( jpjglo-2*jprecj + (jpnj-1) ) / jpnj + 2*jprecj   !: second dimension 
    91    INTEGER, PUBLIC  ::   jpk   ! = jpkdta 
    92    INTEGER, PUBLIC  ::   jpim1 ! = jpi-1                                            !: inner domain indices 
    93    INTEGER, PUBLIC  ::   jpjm1 ! = jpj-1                                            !:   -     -      - 
    94    INTEGER, PUBLIC  ::   jpkm1 ! = jpk-1                                            !:   -     -      - 
    95    INTEGER, PUBLIC  ::   jpij  ! = jpi*jpj                                          !:  jpi x jpj 
     87   INTEGER, PUBLIC, PARAMETER ::   nn_hls = 1   !: halo width (applies to both rows and columns) 
    9688 
    9789   !!---------------------------------------------------------------------- 
    9890   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    99    !! $Id: par_oce.F90 5836 2015-10-26 14:49:40Z cetlod $  
    100    !! Software governed by the CeCILL licence (./LICENSE) 
     91   !! $Id: par_oce.F90 10068 2018-08-28 14:09:04Z nicolasmartin $  
     92   !! Software governed by the CeCILL license (see ./LICENSE) 
    10193   !!====================================================================== 
    10294END MODULE par_oce 
Note: See TracChangeset for help on using the changeset viewer.