Changeset 13204


Ignore:
Timestamp:
2020-07-02T10:38:35+02:00 (5 weeks ago)
Author:
smasson
Message:

tools: update with tools_dev_r12970_AGRIF_CMEMS

Location:
utils/tools
Files:
8 deleted
38 edited
28 copied

Legend:

Unmodified
Added
Removed
  • utils/tools/DOMAINcfg/cpp_DOMAINcfg.fcm

    r12442 r13204  
    1  bld::tool::fppkeys key_mpp_mpi 
     1 bld::tool::fppkeys 
  • utils/tools/DOMAINcfg/namelist_cfg

    r12414 r13204  
    1515&namdom        !   space and time domain (bathymetry, mesh, timestep) 
    1616!----------------------------------------------------------------------- 
     17   nn_bathy    =    1      !  compute analyticaly (=0) or read (=1) the bathymetry file 
     18                           !  or compute (2) from external bathymetry 
     19   nn_interp   =    1                          ! type of interpolation (nn_bathy =2)                        
     20   cn_topo     =  'bathymetry_ORCA12_V3.3.nc'  ! external topo file (nn_bathy =2) 
     21   cn_bath     =  'Bathymetry'                 ! topo name in file  (nn_bathy =2) 
     22   cn_lon      =  'nav_lon'                    ! lon  name in file  (nn_bathy =2) 
     23   cn_lat      =  'nav_lat'                    ! lat  name in file  (nn_bathy =2) 
     24   rn_scale    = 1 
     25   rn_bathy    =    0.     !  value of the bathymetry. if (=0) bottom flat at jpkm1 
     26   jphgr_msh   =       0               !  type of horizontal mesh 
     27   ppglam0     =  999999.0             !  longitude of first raw and column T-point (jphgr_msh = 1) 
     28   ppgphi0     =  999999.0             ! latitude  of first raw and column T-point (jphgr_msh = 1) 
     29   ppe1_deg    =  999999.0             !  zonal      grid-spacing (degrees) 
     30   ppe2_deg    =  999999.0             !  meridional grid-spacing (degrees) 
     31   ppe1_m      =  999999.0             !  zonal      grid-spacing (degrees) 
     32   ppe2_m      =  999999.0             !  meridional grid-spacing (degrees) 
     33   ppsur       =   -4762.96143546300   !  ORCA r4, r2 and r05 coefficients 
     34   ppa0        =     255.58049070440   ! (default coefficients) 
     35   ppa1        =     245.58132232490   ! 
     36   ppkth       =      21.43336197938   ! 
     37   ppacr       =       3.0             ! 
     38   ppdzmin     =  999999.              !  Minimum vertical spacing 
     39   pphmax      =  999999.              !  Maximum depth 
     40   ldbletanh   =  .FALSE.              !  Use/do not use double tanf function for vertical coordinates 
     41   ppa2        =  999999.              !  Double tanh function parameters 
     42   ppkth2      =  999999.              ! 
     43   ppacr2      =  999999.              ! 
    1744/ 
    1845!----------------------------------------------------------------------- 
    1946&namcfg        !   parameters of the configuration 
    2047!----------------------------------------------------------------------- 
     48   ! 
     49   ln_e3_dep   = .true.    ! =T : e3=dk[depth] in discret sens. 
     50   !                       !      ===>>> will become the only possibility in v4.0 
     51   !                       ! =F : e3 analytical derivative of depth function 
     52   !                       !      only there for backward compatibility test with v3.6 
     53   !                       ! 
     54   cp_cfg      =  "orca"   !  name of the configuration 
     55   jp_cfg      =       2   !  resolution of the configuration 
     56   jpidta      =     182   !  1st lateral dimension ( >= jpi ) 
     57   jpjdta      =     149   !  2nd    "         "    ( >= jpj ) 
     58   jpkdta      =      31   !  number of levels      ( >= jpk ) 
     59   jpiglo      =     182   !  1st dimension of global domain --> i =jpidta 
     60   jpjglo      =     149   !  2nd    -                  -    --> j  =jpjdta 
     61   jperio      =       4   !  lateral cond. type (between 0 and 6) 
     62   ln_use_jattr = .false.  !  use (T) the file attribute: open_ocean_jstart, if present 
     63                           !  in netcdf input files, as the start j-row for reading 
     64   ln_domclo = .false.     ! computation of closed sea masks (see namclo) 
    2165/ 
    2266!----------------------------------------------------------------------- 
    2367&namzgr        !   vertical coordinate                                  (default: NO selection) 
    2468!----------------------------------------------------------------------- 
     69!----------------------------------------------------------------------- 
     70   ln_zco      = .false.   !  z-coordinate - full    steps 
     71   ln_zps      = .true.   !  z-coordinate - partial steps 
     72   ln_sco      = .false.   !  s- or hybrid z-s-coordinate 
     73   ln_isfcav   = .false.   !  ice shelf cavity             (T: see namzgr_isf) 
    2574/ 
    2675!----------------------------------------------------------------------- 
     
    4594/ 
    4695!----------------------------------------------------------------------- 
    47 &nambdy        !  unstructured open boundaries                          (default: OFF) 
    48 !----------------------------------------------------------------------- 
    49 / 
    50 !----------------------------------------------------------------------- 
    5196&namnc4        !   netcdf4 chunking and compression settings            ("key_netcdf4") 
    5297!----------------------------------------------------------------------- 
     
    55100&nammpp        !   Massively Parallel Processing                        ("key_mpp_mpi") 
    56101!----------------------------------------------------------------------- 
     102jpni = 0 
     103jpnj=0 
    57104/ 
    58 !----------------------------------------------------------------------- 
    59 &namctl        !   Control prints                                       (default: OFF) 
    60 !----------------------------------------------------------------------- 
    61 / 
  • utils/tools/DOMAINcfg/namelist_ref

    r12414 r13204  
    99&namrun        !   parameters of the run 
    1010!----------------------------------------------------------------------- 
    11    nn_no       =       0   !  Assimilation cycle index 
    1211   cn_exp      =  "ORCA2"  !  experience name 
    1312   nn_it000    =       1   !  first time step 
     
    1615   nn_time0    =       0   !  initial time of day in hhmm 
    1716   nn_leapy    =       0   !  Leap year calendar (1) or not (0) 
    18    ln_rstart   = .false.   !  start from rest (F) or from a restart file (T) 
    19       nn_euler    =    1      !  = 0 : start with forward time step if ln_rstart=T 
    20       nn_rstctl   =    0      !  restart control ==> activated only if ln_rstart=T 
    21       !                          !    = 0 nn_date0 read in namelist ; nn_it000 : read in namelist 
    22       !                          !    = 1 nn_date0 read in namelist ; nn_it000 : check consistancy between namelist and restart 
    23       !                          !    = 2 nn_date0 read in restart  ; nn_it000 : check consistancy between namelist and restart 
    24       cn_ocerst_in    = "restart"   !  suffix of ocean restart name (input) 
    25       cn_ocerst_indir = "."         !  directory from which to read input ocean restarts 
    26       cn_ocerst_out   = "restart"   !  suffix of ocean restart name (output) 
    27       cn_ocerst_outdir = "."         !  directory in which to write output ocean restarts 
    28    ln_iscpl    = .false.   !  cavity evolution forcing or coupling to ice sheet model 
    29    nn_istate   =       0   !  output the initial state (1) or not (0) 
    30    ln_rst_list = .false.   !  output restarts at list of times using nn_stocklist (T) or at set frequency with nn_stock (F) 
    31    nn_stock    =    5840   !  frequency of creation of a restart file (modulo referenced to 1) 
    32    nn_stocklist = 0,0,0,0,0,0,0,0,0,0 ! List of timesteps when a restart file is to be written 
    33    nn_write    =    5475   !  frequency of write in the output file   (modulo referenced to nn_it000) 
    3417   ln_mskland  = .false.   !  mask land points in NetCDF outputs (costly: + ~15%) 
    35    ln_cfmeta   = .false.   !  output additional data to netCDF files required for compliance with the CF metadata standard 
    3618   ln_clobber  = .true.    !  clobber (overwrite) an existing file 
    3719   nn_chunksz  =       0   !  chunksize (bytes) for NetCDF file (works only with iom_nf90 routines) 
     20   ln_cfmeta   = .false.   !  output additional data to netCDF files required for compliance with the CF metadata standard 
     21   ln_iscpl    = .false.   !  cavity evolution forcing or coupling to ice sheet model 
    3822/ 
    3923!----------------------------------------------------------------------- 
    4024&namdom        !   space and time domain (bathymetry, mesh, timestep) 
    4125!----------------------------------------------------------------------- 
    42    nn_bathy    =    1      !  compute (=0) or read (=1) the bathymetry file 
     26   ln_read_cfg = .false.   !  Read from a domain_cfg file 
     27   nn_bathy    =    1      !  compute analyticaly (=0) or read (=1) the bathymetry file 
     28                           !  or compute (2) from external bathymetry 
     29   nn_interp   =    1                          ! type of interpolation (nn_bathy =2) 
     30   cn_domcfg   = ' '       ! Name of the domain_cfg input file 
     31   cn_topo     =  'bathymetry_ORCA12_V3.3.nc'  ! external topo file (nn_bathy =2) 
     32   cn_bath     =  'Bathymetry'                 ! topo name in file  (nn_bathy =2) 
     33   cn_lon      =  'nav_lon'                    ! lon  name in file  (nn_bathy =2) 
     34   cn_lat      =  'nav_lat'                    ! lat  name in file  (nn_bathy =2) 
    4335   rn_bathy    =    0.     !  value of the bathymetry. if (=0) bottom flat at jpkm1 
    4436   nn_msh      =    0      !  create (=1) a mesh file or not (=0) 
     
    8274   !                       ! =F : e3 analytical derivative of depth function 
    8375   !                       !      only there for backward compatibility test with v3.6 
    84    ! 
    85    cp_cfg      = "default" !  name of the configuration 
    86    cp_cfz      = "no zoom" !  name of the zoom of configuration 
    87    jp_cfg      =      0    !  resolution of the configuration 
    88    jpkdta      =     31    !  number of levels      ( >= jpk ) 
    89    jpiglo      =     10    !  1st dimension of global domain --> i =jpidta 
    90    jpjglo      =     12    !  2nd    -                  -    --> j =jpjdta 
    91    jperio      =      0    !  lateral cond. type (between 0 and 6) 
    92                                  !  = 0 closed                 ;   = 1 cyclic East-West 
    93                                  !  = 2 equatorial symmetric   ;   = 3 North fold T-point pivot 
    94                                  !  = 4 cyclic East-West AND North fold T-point pivot 
    95                                  !  = 5 North fold F-point pivot 
    96                                  !  = 6 cyclic East-West AND North fold F-point pivot 
     76   !                       ! 
     77   cp_cfg      =  "orca"   !  name of the configuration 
     78   jp_cfg      =       2   !  resolution of the configuration 
     79   jpidta      =     182   !  1st lateral dimension ( >= jpi ) 
     80   jpjdta      =     149   !  2nd    "         "    ( >= jpj ) 
     81   jpkdta      =      31   !  number of levels      ( >= jpk ) 
     82   jpiglo      =     182   !  1st dimension of global domain --> i =jpidta 
     83   jpjglo      =     149   !  2nd    -                  -    --> j  =jpjdta 
     84   jperio      =       4   !  lateral cond. type (between 0 and 6) 
    9785   ln_use_jattr = .false.  !  use (T) the file attribute: open_ocean_jstart, if present 
    9886                           !  in netcdf input files, as the start j-row for reading 
     
    10694   ln_sco      = .false.   !  s- or hybrid z-s-coordinate 
    10795   ln_isfcav   = .false.   !  ice shelf cavity             (T: see namzgr_isf) 
    108    ln_linssh   = .false.   !  linear free surface 
    10996/ 
    11097!----------------------------------------------------------------------- 
     
    182169   rn_sponge_tra = 2880.   !  coefficient for tracer   sponge layer [m2/s] 
    183170   rn_sponge_dyn = 2880.   !  coefficient for dynamics sponge layer [m2/s] 
    184    ln_chk_bathy  = .false. !  =T  check the parent bathymetry 
    185 / 
    186 !----------------------------------------------------------------------- 
    187 &nambdy        !  unstructured open boundaries                          (default: OFF) 
    188 !----------------------------------------------------------------------- 
    189    ln_bdy         = .false.   !  Use unstructured open boundaries 
    190    nb_bdy         = 0         !  number of open boundary sets 
    191    ln_coords_file = .true.    !  =T : read bdy coordinates from file 
    192       cn_coords_file = 'coordinates.bdy.nc'  !  bdy coordinates files 
    193    ln_mask_file   = .false.   !  =T : read mask from file 
    194       cn_mask_file = ''        !  name of mask file (if ln_mask_file=.TRUE.) 
    195    cn_dyn2d    = 'none'       ! 
    196    nn_dyn2d_dta   =  0        !  = 0, bdy data are equal to the initial state 
    197       !                       !  = 1, bdy data are read in 'bdydata   .nc' files 
    198       !                       !  = 2, use tidal harmonic forcing data from files 
    199       !                       !  = 3, use external data AND tidal harmonic forcing 
    200    cn_dyn3d      =  'none'    ! 
    201    nn_dyn3d_dta  =  0         !  = 0, bdy data are equal to the initial state 
    202    !                          !  = 1, bdy data are read in 'bdydata   .nc' files 
    203    cn_tra        =  'none'    ! 
    204    nn_tra_dta    =  0         !  = 0, bdy data are equal to the initial state 
    205    !                          !  = 1, bdy data are read in 'bdydata   .nc' files 
    206    cn_ice        =  'none'    ! 
    207    nn_ice_dta    =  0         !  = 0, bdy data are equal to the initial state 
    208    !                          !  = 1, bdy data are read in 'bdydata   .nc' files 
    209    rn_ice_tem    = 270.       !  si3 only: arbitrary temperature of incoming sea ice 
    210    rn_ice_sal    = 10.        !  si3 only:      --   salinity           -- 
    211    rn_ice_age    = 30.        !  si3 only:      --   age                -- 
    212    ! 
    213    ln_tra_dmp    =.false.     !  open boudaries conditions for tracers 
    214    ln_dyn3d_dmp  =.false.     !  open boundary condition for baroclinic velocities 
    215    rn_time_dmp   =  1.        !  Damping time scale in days 
    216    rn_time_dmp_out = 1.        !  Outflow damping time scale 
    217    nn_rimwidth   = 10         !  width of the relaxation zone 
    218    ln_vol        = .false.    !  total volume correction (see nn_volctl parameter) 
    219    nn_volctl     =  1         !  = 0, the total water flux across open boundaries is zero 
    220    nb_jpk_bdy    = -1         ! number of levels in the bdy data (set < 0 if consistent with planned run) 
    221 / 
    222 !----------------------------------------------------------------------- 
    223 &namnc4        !   netcdf4 chunking and compression settings            ("key_netcdf4") 
    224 !----------------------------------------------------------------------- 
    225    nn_nchunks_i =   4       !  number of chunks in i-dimension 
    226    nn_nchunks_j =   4       !  number of chunks in j-dimension 
    227    nn_nchunks_k =   31      !  number of chunks in k-dimension 
    228    !                       !  setting nn_nchunks_k = jpk will give a chunk size of 1 in the vertical which 
    229    !                       !  is optimal for postprocessing which works exclusively with horizontal slabs 
    230    ln_nc4zip   = .true.    !  (T) use netcdf4 chunking and compression 
    231    !                       !  (F) ignore chunking information and produce netcdf3-compatible files 
     171   ln_chk_bathy  = .FALSE. ! 
     172   npt_connect   = 2 
     173   npt_copy      = 2 
    232174/ 
    233175!----------------------------------------------------------------------- 
     
    238180   nn_buffer   =   0       !  size in bytes of exported buffer ('B' case), 0 no exportation 
    239181   ln_nnogather =  .true.  !  activate code to avoid mpi_allgather use at the northfold 
    240    jpni        =   0       !  jpni   number of processors following i (set automatically if < 1) 
    241    jpnj        =   0       !  jpnj   number of processors following j (set automatically if < 1) 
     182   jpni        =   1       !  jpni   number of processors following i (set automatically if < 1) 
     183   jpnj        =   1       !  jpnj   number of processors following j (set automatically if < 1) 
    242184/ 
    243 !----------------------------------------------------------------------- 
    244 &namctl        !   Control prints                                       (default: OFF) 
    245 !----------------------------------------------------------------------- 
    246    ln_ctl = .FALSE.                 ! Toggle all report printing on/off (T/F); Ignored if sn_cfctl%l_config is T 
    247      sn_cfctl%l_config = .TRUE.     ! IF .true. then control which reports are written with the following 
    248        sn_cfctl%l_runstat = .FALSE. ! switches and which areas produce reports with the proc integer settings. 
    249        sn_cfctl%l_trcstat = .FALSE. ! The default settings for the proc integers should ensure 
    250        sn_cfctl%l_oceout  = .FALSE. ! that  all areas report. 
    251        sn_cfctl%l_layout  = .FALSE. ! 
    252        sn_cfctl%l_mppout  = .FALSE. ! 
    253        sn_cfctl%l_mpptop  = .FALSE. ! 
    254        sn_cfctl%procmin   = 0       ! Minimum area number for reporting [default:0] 
    255        sn_cfctl%procmax   = 1000000 ! Maximum area number for reporting [default:1000000] 
    256        sn_cfctl%procincr  = 1       ! Increment for optional subsetting of areas [default:1] 
    257        sn_cfctl%ptimincr  = 1       ! Timestep increment for writing time step progress info 
    258    nn_print    =    0      !  level of print (0 no extra print) 
    259    nn_ictls    =    0      !  start i indice of control sum (use to compare mono versus 
    260    nn_ictle    =    0      !  end   i indice of control sum        multi processor runs 
    261    nn_jctls    =    0      !  start j indice of control               over a subdomain) 
    262    nn_jctle    =    0      !  end   j indice of control 
    263    nn_isplt    =    1      !  number of processors in i-direction 
    264    nn_jsplt    =    1      !  number of processors in j-direction 
    265    ln_timing   = .false.   !  timing by routine write out in timing.output file 
    266    ln_diacfl   = .false.   !  CFL diagnostics write out in cfl_diagnostics.ascii 
    267 / 
  • utils/tools/DOMAINcfg/src/agrif_domzgr.F90

    r12414 r13204  
     1MODULE agrif_domzgr 
     2 
     3   USE agrif_profiles 
     4   USE dom_oce 
     5 
     6   IMPLICIT NONE 
     7   PRIVATE 
     8   
     9   PUBLIC :: agrif_create_bathy_meter 
     10  
     11 
     12CONTAINS 
     13 
    114#if defined key_agrif 
    2 subroutine agrif_domzgr 
    3 end subroutine agrif_domzgr 
    415 
    5 subroutine agrif_create_bathy_meter 
    6 use agrif_profiles 
    7 external :: init_bathy 
     16   SUBROUTINE agrif_create_bathy_meter 
    817 
    9 call Agrif_Init_variable(bathy_id, procname = init_bathy) 
     18      CALL Agrif_Init_variable(bathy_id, procname = init_bathy) 
    1019 
    11 end subroutine agrif_create_bathy_meter 
     20   END SUBROUTINE agrif_create_bathy_meter 
    1221 
    13     SUBROUTINE init_bathy( ptab, i1, i2, j1, j2, before, nb,ndir) 
    14     use dom_oce 
     22   SUBROUTINE init_bathy( ptab, i1, i2, j1, j2, before, nb,ndir) 
    1523      !!---------------------------------------------------------------------- 
    1624      !!                  ***  ROUTINE interpsshn  *** 
     
    2028      LOGICAL                         , INTENT(in   ) ::   before 
    2129      INTEGER                         , INTENT(in   ) ::   nb , ndir 
    22       LOGICAL  ::   western_side, eastern_side,northern_side,southern_side 
    2330      ! 
    2431      !!----------------------------------------------------------------------   
    25       ! 
    26          western_side  = (nb == 1).AND.(ndir == 1) 
    27          eastern_side  = (nb == 1).AND.(ndir == 2) 
    28          southern_side = (nb == 2).AND.(ndir == 1) 
    29          northern_side = (nb == 2).AND.(ndir == 2) 
     32      INTEGER :: ji,jj 
     33 
    3034      IF( before) THEN 
    3135         ptab(i1:i2,j1:j2) = bathy(i1:i2,j1:j2) 
     36         DO jj=j1,j2 
     37            DO ji=i1,i2 
     38               ptab(ji,jj) = SUM( e3t_0(ji,jj, 1:mbkt(ji,jj) ) ) * ssmask(ji,jj) 
     39            END DO 
     40         END DO 
    3241      ELSE 
    3342         bathy(i1:i2,j1:j2)=ptab 
     
    3544      ! 
    3645   END SUBROUTINE init_bathy 
     46 
    3747#else 
    38 subroutine agrif_domzgr_empty 
    39 end subroutine agrif_domzgr_empty 
     48   SUBROUTINE agrif_create_bathy_meter 
     49   END SUBROUTINE agrif_create_bathy_meter 
    4050#endif 
     51END MODULE agrif_domzgr 
  • utils/tools/DOMAINcfg/src/agrif_parameters.F90

    r12414 r13204  
    1 #ifdef key_agrif 
    2 module agrif_parameters 
    3 USE par_kind 
     1MODULE agrif_parameters 
     2    
     3   USE par_kind 
    44 
    5 INTEGER :: nn_cln_update 
    6 LOGICAL :: ln_spc_dyn 
    7 REAL(wp) :: rn_sponge_tra 
    8 REAL(wp) :: rn_sponge_dyn 
    9 LOGICAL :: ln_chk_bathy 
    10 INTEGER :: npt_copy 
    11 INTEGER :: npt_connect 
     5   PUBLIC  
    126 
    13 REAL(wp), PUBLIC, ALLOCATABLE, SAVE        , DIMENSION(:,:) ::   ztabramp 
     7#if defined key_agrif 
    148 
    15 end module agrif_parameters 
    16 #else 
    17 subroutine agrif_parameters_empty 
    18 end subroutine agrif_parameters_empty 
     9   INTEGER :: nn_cln_update 
     10   LOGICAL :: ln_spc_dyn 
     11   REAL(wp) :: rn_sponge_tra 
     12   REAL(wp) :: rn_sponge_dyn 
     13   LOGICAL :: ln_chk_bathy 
     14   INTEGER :: npt_copy 
     15   INTEGER :: npt_connect 
     16   REAL(wp), PUBLIC, ALLOCATABLE, SAVE        , DIMENSION(:,:) ::   ztabramp 
     17   REAL(wp), PUBLIC, ALLOCATABLE, SAVE        , DIMENSION(:,:,:) ::   e3t_interp 
     18 
    1919#endif 
     20 
     21END MODULE agrif_parameters 
  • utils/tools/DOMAINcfg/src/agrif_profiles.F90

    r12414 r13204  
    1 module agrif_profiles 
     1MODULE agrif_profiles 
    22 
    3 integer :: glamt_id, glamu_id, glamv_id,glamf_id 
    4 integer :: gphit_id, gphiu_id, gphiv_id,gphif_id 
    5 integer :: e1t_id, e1u_id, e1v_id, e1f_id 
    6 integer :: e2t_id, e2u_id, e2v_id, e2f_id 
     3  PUBLIC 
     4 
     5#if defined key_agrif 
     6 
     7INTEGER :: glamt_id, glamu_id, glamv_id,glamf_id 
     8INTEGER :: gphit_id, gphiu_id, gphiv_id,gphif_id 
     9INTEGER :: e1t_id, e1u_id, e1v_id, e1f_id 
     10INTEGER :: e2t_id, e2u_id, e2v_id, e2f_id 
    711 
    812 
    9 integer :: bathy_id 
     13INTEGER :: bathy_id 
    1014 
    1115! Vertical scale factors 
    1216 
    13 integer :: e3t_id 
    14 integer :: e3t_copy_id 
    15 integer :: e3t_connect_id 
    16 integer :: e3u_id, e3v_id 
     17INTEGER :: e3t_id 
     18INTEGER :: e3t_copy_id 
     19INTEGER :: e3t_connect_id 
     20INTEGER :: e3u_id, e3v_id 
    1721 
    1822! Bottom level 
    19 integer :: bottom_level_id 
    20 end module agrif_profiles 
     23INTEGER :: bottom_level_id 
     24 
     25# endif 
     26END MODULE agrif_profiles 
  • utils/tools/DOMAINcfg/src/agrif_user.F90

    r12414 r13204  
    11#if defined key_agrif 
    2 subroutine agrif_initworkspace() 
     2   SUBROUTINE agrif_user() 
     3      !!---------------------------------------------------------------------- 
     4      !!                 *** ROUTINE agrif_user *** 
     5      !!---------------------------------------------------------------------- 
     6   END SUBROUTINE agrif_user 
     7 
     8   SUBROUTINE agrif_initworkspace() 
    39      !!---------------------------------------------------------------------- 
    410      !!                 *** ROUTINE Agrif_InitWorkspace *** 
    511      !!---------------------------------------------------------------------- 
    6 end subroutine agrif_initworkspace 
    7 SUBROUTINE Agrif_InitValues 
     12   END SUBROUTINE agrif_initworkspace 
     13 
     14   SUBROUTINE Agrif_InitValues 
    815      !!---------------------------------------------------------------------- 
    916      !!                 *** ROUTINE Agrif_InitValues *** 
     
    1118      !! ** Purpose :: Declaration of variables to be interpolated 
    1219      !!---------------------------------------------------------------------- 
    13    USE Agrif_Util 
    14    USE dom_oce 
    15    USE nemogcm 
    16    USE domain 
    17    !! 
    18    IMPLICIT NONE 
    19     
    20   
    21    CALL nemo_init       !* Initializations of each fine grid 
    22  
    23    CALL dom_nam 
    24    CALL cfg_write         ! create the configuration file 
    25     
    26 END SUBROUTINE Agrif_InitValues 
    27  
    28 SUBROUTINE Agrif_InitValues_cont 
    29  
    30 use dom_oce 
    31     integer :: irafx, irafy 
    32     logical :: ln_perio 
    33     integer nx,ny 
    34  
    35 irafx = agrif_irhox() 
    36 irafy = agrif_irhoy() 
    37  
    38 nx=nlci ; ny=nlcj 
     20      USE Agrif_Util 
     21      USE dom_oce 
     22      USE nemogcm 
     23      USE domain 
     24      !! 
     25      IMPLICIT NONE 
     26 
     27      ! No temporal refinement 
     28      CALL Agrif_Set_coeffreft(1) 
     29 
     30      CALL nemo_init       !* Initializations of each fine grid 
     31 
     32      CALL dom_nam 
     33 
     34   END SUBROUTINE Agrif_InitValues 
     35 
     36   SUBROUTINE Agrif_InitValues_cont 
     37      !!---------------------------------------------------------------------- 
     38      !!                 *** ROUTINE Agrif_InitValues_cont *** 
     39      !! 
     40      !! ** Purpose :: Initialisation of variables to be interpolated 
     41      !!---------------------------------------------------------------------- 
     42      USE dom_oce 
     43      USE lbclnk 
     44      !! 
     45      IMPLICIT NONE 
     46      ! 
     47      INTEGER :: nx, ny 
     48      INTEGER :: irafx, irafy 
     49      LOGICAL :: ln_perio 
     50      ! 
     51      irafx = agrif_irhox() 
     52      irafy = agrif_irhoy() 
     53 
     54      nx = nlci ; ny = nlcj 
    3955 
    4056   !       IF(jperio /=1 .AND. jperio/=4 .AND. jperio/=6 ) THEN 
     
    4460   !          nbghostcellsfine_tot_y= nbghostcellsfine_y +1 
    4561   !       ELSE 
    46    !         nx = (nbcellsx)+2*nbghostcellsfine_x  
     62   !         nx = (nbcellsx)+2*nbghostcellsfine_x 
    4763   !         ny = (nbcellsy)+2*nbghostcellsfine+2 
    4864   !         nbghostcellsfine_tot_x= 1 
     
    5369   !       nx = nbcellsx+irafx 
    5470   !       ny = nbcellsy+irafy 
    55         
    56   WRITE(*,*) ' ' 
    57   WRITE(*,*)'Size of the High resolution grid: ',nx,' x ',ny 
    58   WRITE(*,*) ' ' 
    59         
    60        call agrif_init_lonlat() 
    61        ln_perio=.FALSE.  
    62        if( jperio ==1 .OR. jperio==2 .OR. jperio==4) ln_perio=.TRUE. 
    63  
    64        where (glamt < -180) glamt = glamt +360. 
    65        if (ln_perio) then 
    66          glamt(1,:)=glamt(nx-1,:) 
    67          glamt(nx,:)=glamt(2,:) 
    68        endif 
    69   
    70        where (glamu < -180) glamu = glamu +360. 
    71        if (ln_perio) then 
    72          glamu(1,:)=glamu(nx-1,:) 
    73          glamu(nx,:)=glamu(2,:) 
    74        endif 
    75  
    76       where (glamv < -180) glamv = glamv +360. 
    77        if (ln_perio) then 
    78          glamv(1,:)=glamv(nx-1,:) 
    79          glamv(nx,:)=glamv(2,:) 
    80        endif 
    81  
    82       where (glamf < -180) glamf = glamf +360. 
    83        if (ln_perio) then 
    84          glamf(1,:)=glamf(nx-1,:) 
    85          glamf(nx,:)=glamf(2,:) 
    86        endif 
    87  
    88        call agrif_init_scales() 
    89  
    90          
    91 END SUBROUTINE Agrif_InitValues_cont   
    92  
    93  
    94 subroutine agrif_declare_var() 
    95 use par_oce 
    96 use dom_oce 
    97 use agrif_profiles 
    98 use agrif_parameters 
    99  
    100    IMPLICIT NONE 
    101     
    102 integer :: ind1, ind2, ind3 
    103 integer nx,ny 
    104 integer nbghostcellsfine_tot_x, nbghostcellsfine_tot_y 
    105 INTEGER :: irafx 
    106 !!---------------------------------------------------------------------- 
    107  
    108    ! 1. Declaration of the type of variable which have to be interpolated 
    109    !--------------------------------------------------------------------- 
    110  nx=nlci ; ny=nlcj 
    111  
    112 !ind2 = nbghostcellsfine_tot_x + 1 
    113 !ind3 = nbghostcellsfine_tot_y + 1 
    114 ind2 = 2 + nbghostcells 
    115 ind3 = ind2 
    116 nbghostcellsfine_tot_x=nbghostcells+1 
    117 nbghostcellsfine_tot_y=nbghostcells+1 
    118  
    119 irafx = Agrif_irhox() 
    120  
    121 CALL agrif_nemo_init  ! specific namelist part if needed 
    122  
    123 CALL agrif_declare_variable((/2,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),glamt_id) 
    124 CALL agrif_declare_variable((/1,2/),(/ind2-1,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),glamu_id) 
    125 CALL agrif_declare_variable((/2,1/),(/ind2,ind3-1/),(/'x','y'/),(/1,1/),(/nx,ny/),glamv_id) 
    126 CALL agrif_declare_variable((/1,1/),(/ind2-1,ind3-1/),(/'x','y'/),(/1,1/),(/nx,ny/),glamf_id) 
    127  
    128 CALL agrif_declare_variable((/2,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),gphit_id) 
    129 CALL agrif_declare_variable((/1,2/),(/ind2-1,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),gphiu_id) 
    130 CALL agrif_declare_variable((/2,1/),(/ind2,ind3-1/),(/'x','y'/),(/1,1/),(/nx,ny/),gphiv_id) 
    131 CALL agrif_declare_variable((/1,1/),(/ind2-1,ind3-1/),(/'x','y'/),(/1,1/),(/nx,ny/),gphif_id) 
    132  
    133 ! Horizontal scale factors 
    134  
    135 CALL agrif_declare_variable((/2,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),e1t_id) 
    136 CALL agrif_declare_variable((/1,2/),(/ind2-1,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),e1u_id) 
    137 CALL agrif_declare_variable((/2,1/),(/ind2,ind3-1/),(/'x','y'/),(/1,1/),(/nx,ny/),e1v_id) 
    138 CALL agrif_declare_variable((/1,1/),(/ind2-1,ind3-1/),(/'x','y'/),(/1,1/),(/nx,ny/),e1f_id) 
    139  
    140 CALL agrif_declare_variable((/2,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),e2t_id) 
    141 CALL agrif_declare_variable((/1,2/),(/ind2-1,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),e2u_id) 
    142 CALL agrif_declare_variable((/2,1/),(/ind2,ind3-1/),(/'x','y'/),(/1,1/),(/nx,ny/),e2v_id) 
    143 CALL agrif_declare_variable((/1,1/),(/ind2-1,ind3-1/),(/'x','y'/),(/1,1/),(/nx,ny/),e2f_id) 
    144  
    145 ! Bathymetry 
    146  
    147 CALL agrif_declare_variable((/2,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),bathy_id) 
    148  
    149 ! Vertical scale factors 
    150 CALL agrif_declare_variable((/2,2,0/),(/ind2,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/nx,ny,jpk/),e3t_id) 
    151 CALL agrif_declare_variable((/2,2,0/),(/ind2,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/nx,ny,jpk/),e3t_copy_id) 
    152 CALL agrif_declare_variable((/2,2,0/),(/ind2,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/nx,ny,jpk+1/),e3t_connect_id) 
    153  
    154 CALL agrif_declare_variable((/1,2,0/),(/ind2-1,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/nx,ny,jpk/),e3u_id) 
    155 CALL agrif_declare_variable((/2,1,0/),(/ind2,ind3-1,0/),(/'x','y','N'/),(/1,1,1/),(/nx,ny,jpk/),e3v_id) 
    156  
    157 ! Bottom level 
    158  
    159 CALL agrif_declare_variable((/2,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),bottom_level_id) 
    160  
    161 CALL Agrif_Set_bcinterp(glamt_id,interp=AGRIF_linear) 
    162 CALL Agrif_Set_interp(glamt_id,interp=AGRIF_linear) 
    163 CALL Agrif_Set_bc( glamt_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 
    164  
    165 CALL Agrif_Set_bcinterp(glamu_id,interp=AGRIF_linear) 
    166 CALL Agrif_Set_interp(glamu_id,interp=AGRIF_linear) 
    167 CALL Agrif_Set_bc( glamu_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 
    168  
    169 CALL Agrif_Set_bcinterp(glamv_id,interp=AGRIF_linear) 
    170 CALL Agrif_Set_interp(glamv_id,interp=AGRIF_linear) 
    171 CALL Agrif_Set_bc( glamv_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 
    172  
    173 CALL Agrif_Set_bcinterp(glamf_id,interp=AGRIF_linear) 
    174 CALL Agrif_Set_interp(glamf_id,interp=AGRIF_linear) 
    175 CALL Agrif_Set_bc( glamf_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 
    176  
    177 CALL Agrif_Set_bcinterp(gphit_id,interp=AGRIF_linear) 
    178 CALL Agrif_Set_interp(gphit_id,interp=AGRIF_linear) 
    179 CALL Agrif_Set_bc( gphit_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 
    180  
    181 CALL Agrif_Set_bcinterp(gphiu_id,interp=AGRIF_linear) 
    182 CALL Agrif_Set_interp(gphiu_id,interp=AGRIF_linear) 
    183 CALL Agrif_Set_bc( gphiu_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 
    184  
    185 CALL Agrif_Set_bcinterp(gphiv_id,interp=AGRIF_linear) 
    186 CALL Agrif_Set_interp(gphiv_id,interp=AGRIF_linear) 
    187 CALL Agrif_Set_bc( gphiv_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 
    188  
    189 CALL Agrif_Set_bcinterp(gphif_id,interp=AGRIF_linear) 
    190 CALL Agrif_Set_interp(gphif_id,interp=AGRIF_linear) 
    191 CALL Agrif_Set_bc( gphif_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 
    192  
    193 ! 
    194  
    195 CALL Agrif_Set_bcinterp(e1t_id,interp=AGRIF_ppm) 
    196 CALL Agrif_Set_interp(e1t_id,interp=AGRIF_ppm) 
    197 CALL Agrif_Set_bc( e1t_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 
    198  
    199 CALL Agrif_Set_bcinterp(e1u_id, interp1=Agrif_linear, interp2=AGRIF_ppm) 
    200 CALL Agrif_Set_interp(e1u_id, interp1=Agrif_linear, interp2=AGRIF_ppm) 
    201 CALL Agrif_Set_bc( e1u_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 
    202  
    203 CALL Agrif_Set_bcinterp(e1v_id,interp1=AGRIF_ppm, interp2=Agrif_linear) 
    204 CALL Agrif_Set_interp(e1v_id, interp1=AGRIF_ppm, interp2=Agrif_linear) 
    205 CALL Agrif_Set_bc( e1v_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 
    206  
    207 CALL Agrif_Set_bcinterp(e1f_id,interp=AGRIF_linear) 
    208 CALL Agrif_Set_interp(e1f_id,interp=AGRIF_linear) 
    209 CALL Agrif_Set_bc( e1f_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 
    210  
    211 CALL Agrif_Set_bcinterp(e2t_id,interp=AGRIF_ppm) 
    212 CALL Agrif_Set_interp(e2t_id,interp=AGRIF_ppm) 
    213 CALL Agrif_Set_bc( e2t_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 
    214  
    215 CALL Agrif_Set_bcinterp(e2u_id,interp1=Agrif_linear, interp2=AGRIF_ppm) 
    216 CALL Agrif_Set_interp(e2u_id,interp1=Agrif_linear, interp2=AGRIF_ppm) 
    217 CALL Agrif_Set_bc( e2u_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 
    218  
    219 CALL Agrif_Set_bcinterp(e2v_id,interp1=AGRIF_ppm, interp2=Agrif_linear) 
    220 CALL Agrif_Set_interp(e2v_id,interp1=AGRIF_ppm, interp2=Agrif_linear) 
    221 CALL Agrif_Set_bc( e2v_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 
    222  
    223 CALL Agrif_Set_bcinterp(e2f_id,interp=AGRIF_linear) 
    224 CALL Agrif_Set_interp(e2f_id,interp=AGRIF_linear) 
    225 CALL Agrif_Set_bc( e2f_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 
    226  
    227 CALL Agrif_Set_bcinterp(bathy_id,interp=AGRIF_linear) 
    228 CALL Agrif_Set_interp(bathy_id,interp=AGRIF_linear) 
    229 CALL Agrif_Set_bc( bathy_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 
    230  
    231 ! Vertical scale factors 
    232 CALL Agrif_Set_bcinterp(e3t_id,interp=AGRIF_ppm) 
    233 CALL Agrif_Set_interp(e3t_id,interp=AGRIF_ppm) 
    234 CALL Agrif_Set_bc( e3t_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 
    235 CALL Agrif_Set_Updatetype( e3t_id, update = AGRIF_Update_Average) 
    236  
    237 CALL Agrif_Set_bcinterp(e3t_copy_id,interp=AGRIF_constant) 
    238 CALL Agrif_Set_interp(e3t_copy_id,interp=AGRIF_constant) 
    239 CALL Agrif_Set_bc( e3t_copy_id, (/-npt_copy*irafx-1,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/)) 
    240  
    241 CALL Agrif_Set_bcinterp(e3t_connect_id,interp=AGRIF_ppm) 
    242 CALL Agrif_Set_interp(e3t_connect_id,interp=AGRIF_ppm) 
    243 CALL Agrif_Set_bc( e3t_connect_id, (/-(npt_copy+npt_connect)*irafx-1,-npt_copy*irafx-2/)) 
    244  
    245 CALL Agrif_Set_bcinterp(e3u_id, interp1=Agrif_linear, interp2=AGRIF_ppm) 
    246 CALL Agrif_Set_interp(e3u_id, interp1=Agrif_linear, interp2=AGRIF_ppm) 
    247 CALL Agrif_Set_bc( e3u_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 
    248 CALL Agrif_Set_Updatetype(e3u_id,update1 = Agrif_Update_Copy, update2 = Agrif_Update_Average) 
    249  
    250 CALL Agrif_Set_bcinterp(e3v_id,interp1=AGRIF_ppm, interp2=Agrif_linear) 
    251 CALL Agrif_Set_interp(e3v_id, interp1=AGRIF_ppm, interp2=Agrif_linear) 
    252 CALL Agrif_Set_bc( e3v_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 
    253 CALL Agrif_Set_Updatetype(e3v_id,update1 = Agrif_Update_Average, update2 = Agrif_Update_Copy) 
    254     
    255 ! Bottom level 
    256 CALL Agrif_Set_bcinterp(bottom_level_id,interp=AGRIF_constant) 
    257 CALL Agrif_Set_interp(bottom_level_id,interp=AGRIF_constant) 
    258 CALL Agrif_Set_bc( bottom_level_id, (/-npt_copy*irafx-1,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/)) 
    259 CALL Agrif_Set_Updatetype( bottom_level_id, update = AGRIF_Update_Max) 
    260  
    261 end subroutine agrif_declare_var 
    262  
    263  
    264 subroutine agrif_init_lonlat() 
    265 use agrif_profiles 
    266 use agrif_util 
    267 external :: init_glamt, init_glamu, init_glamv, init_glamf 
    268 external :: init_gphit, init_gphiu, init_gphiv, init_gphif 
    269  
    270 call Agrif_Init_variable(glamt_id, procname = init_glamt) 
    271 call Agrif_Init_variable(glamu_id, procname = init_glamu) 
    272 call Agrif_Init_variable(glamv_id, procname = init_glamv) 
    273 call Agrif_Init_variable(glamf_id, procname = init_glamf) 
    274  
    275 call Agrif_Init_variable(gphit_id, procname = init_gphit) 
    276 call Agrif_Init_variable(gphiu_id, procname = init_gphiu) 
    277 call Agrif_Init_variable(gphiv_id, procname = init_gphiv) 
    278 call Agrif_Init_variable(gphif_id, procname = init_gphif) 
    279  
    280 end subroutine agrif_init_lonlat 
    281  
    282 subroutine agrif_init_scales() 
    283 use agrif_profiles 
    284 use agrif_util 
    285 external :: init_e1t, init_e1u, init_e1v, init_e1f 
    286 external :: init_e2t, init_e2u, init_e2v, init_e2f 
    287  
    288 call Agrif_Init_variable(e1t_id, procname = init_e1t) 
    289 call Agrif_Init_variable(e1u_id, procname = init_e1u) 
    290 call Agrif_Init_variable(e1v_id, procname = init_e1v) 
    291 call Agrif_Init_variable(e1f_id, procname = init_e1f) 
    292  
    293 call Agrif_Init_variable(e2t_id, procname = init_e2t) 
    294 call Agrif_Init_variable(e2u_id, procname = init_e2u) 
    295 call Agrif_Init_variable(e2v_id, procname = init_e2v) 
    296 call Agrif_Init_variable(e2f_id, procname = init_e2f) 
    297  
    298 end subroutine agrif_init_scales 
    299  
    300  
    301  
    302    SUBROUTINE init_glamt( ptab, i1, i2, j1, j2, before, nb,ndir) 
    303    use dom_oce 
     71 
     72      WRITE(*,*) ' ' 
     73      WRITE(*,*)'Size of the High resolution grid: ',nx,' x ',ny 
     74      WRITE(*,*) ' ' 
     75 
     76      CALL agrif_init_lonlat() 
     77      ln_perio = .FALSE. 
     78      IF( jperio == 1 .OR. jperio == 2 .OR. jperio == 4 ) ln_perio=.TRUE. 
     79 
     80      WHERE (glamt < -180) glamt = glamt +360. 
     81      WHERE (glamt > +180) glamt = glamt -360. 
     82 
     83      CALL lbc_lnk( 'glamt', glamt, 'T', 1._wp) 
     84      CALL lbc_lnk( 'gphit', gphit, 'T', 1._wp) 
     85 
     86      WHERE (glamu < -180) glamu = glamu +360. 
     87      WHERE (glamu > +180) glamu = glamu -360. 
     88      CALL lbc_lnk( 'glamu', glamu, 'U', 1._wp) 
     89      CALL lbc_lnk( 'gphiu', gphiu, 'U', 1._wp) 
     90 
     91      WHERE (glamv < -180) glamv = glamv +360. 
     92      WHERE (glamv > +180) glamv = glamv -360. 
     93      CALL lbc_lnk( 'glamv', glamv, 'V', 1._wp) 
     94      CALL lbc_lnk( 'gphiv', gphiv, 'V', 1._wp) 
     95 
     96      WHERE (glamf < -180) glamf = glamf +360. 
     97      WHERE (glamf > +180) glamf = glamf -360. 
     98      CALL lbc_lnk( 'glamf', glamf, 'F', 1._wp) 
     99      CALL lbc_lnk( 'gphif', gphif, 'F', 1._wp) 
     100 
     101      ! Correct South and North 
     102      IF ((nbondj == -1).OR.(nbondj == 2)) THEN 
     103         glamt(:,1) = glamt(:,2) 
     104         gphit(:,1) = gphit(:,2) 
     105         glamu(:,1) = glamu(:,2) 
     106         gphiu(:,1) = gphiu(:,2) 
     107         glamv(:,1) = glamv(:,2) 
     108         gphiv(:,1) = gphiv(:,2) 
     109      ENDIF 
     110      IF ((nbondj == 1).OR.(nbondj == 2)) THEN 
     111         glamt(:,jpj) = glamt(:,jpj-1) 
     112         gphit(:,jpj) = gphit(:,jpj-1) 
     113         glamu(:,jpj) = glamu(:,jpj-1) 
     114         gphiu(:,jpj) = gphiu(:,jpj-1) 
     115         glamv(:,jpj) = glamv(:,jpj-1) 
     116         gphiv(:,jpj) = gphiv(:,jpj-1) 
     117         glamf(:,jpj) = glamf(:,jpj-1) 
     118         gphif(:,jpj) = gphif(:,jpj-1) 
     119      ENDIF 
     120 
     121      ! Correct West and East 
     122      IF( jperio /= 1 ) THEN 
     123         IF((nbondi == -1) .OR. (nbondi == 2) ) THEN 
     124            glamt(1,:) = glamt(2,:) 
     125            gphit(1,:) = gphit(2,:) 
     126            glamu(1,:) = glamu(2,:) 
     127            gphiu(1,:) = gphiu(2,:) 
     128            glamv(1,:) = glamv(2,:) 
     129            gphiv(1,:) = gphiv(2,:) 
     130         ENDIF 
     131         IF( (nbondi == 1) .OR. (nbondi == 2) ) THEN 
     132            glamt(jpi,:) = glamt(jpi-1,:) 
     133            gphit(jpi,:) = gphit(jpi-1,:) 
     134            glamu(jpi,:) = glamu(jpi-1,:) 
     135            gphiu(jpi,:) = gphiu(jpi-1,:) 
     136            glamv(jpi,:) = glamv(jpi-1,:) 
     137            gphiv(jpi,:) = gphiv(jpi-1,:) 
     138            glamf(jpi,:) = glamf(jpi-1,:) 
     139            gphif(jpi,:) = gphif(jpi-1,:) 
     140         ENDIF 
     141      ENDIF 
     142 
     143      CALL agrif_init_scales() 
     144 
     145      ! Correct South and North 
     146      IF( (nbondj == -1) .OR. (nbondj == 2) ) THEN 
     147         e1t(:,1) = e1t(:,2) 
     148         e2t(:,1) = e2t(:,2) 
     149         e1u(:,1) = e1u(:,2) 
     150         e2u(:,1) = e2u(:,2) 
     151         e1v(:,1) = e1v(:,2) 
     152         e2v(:,1) = e2v(:,2) 
     153      ENDIF 
     154      IF( (nbondj == 1) .OR. (nbondj == 2) ) THEN 
     155         e1t(:,jpj) = e1t(:,jpj-1) 
     156         e2t(:,jpj) = e2t(:,jpj-1) 
     157         e1u(:,jpj) = e1u(:,jpj-1) 
     158         e2u(:,jpj) = e2u(:,jpj-1) 
     159         e1v(:,jpj) = e1v(:,jpj-1) 
     160         e2v(:,jpj) = e2v(:,jpj-1) 
     161         e1f(:,jpj) = e1f(:,jpj-1) 
     162         e2f(:,jpj) = e2f(:,jpj-1) 
     163      ENDIF 
     164 
     165      ! Correct West and East 
     166      IF( jperio /= 1 ) THEN 
     167         IF( (nbondj == -1) .OR. (nbondj == 2) ) THEN 
     168            e1t(1,:) = e1t(2,:) 
     169            e2t(1,:) = e2t(2,:) 
     170            e1u(1,:) = e1u(2,:) 
     171            e2u(1,:) = e2u(2,:) 
     172            e1v(1,:) = e1v(2,:) 
     173            e2v(1,:) = e2v(2,:) 
     174         ENDIF 
     175         IF( (nbondj == 1) .OR. (nbondj == 2) ) THEN 
     176            e1t(jpi,:) = e1t(jpi-1,:) 
     177            e2t(jpi,:) = e2t(jpi-1,:) 
     178            e1u(jpi,:) = e1u(jpi-1,:) 
     179            e2u(jpi,:) = e2u(jpi-1,:) 
     180            e1v(jpi,:) = e1v(jpi-1,:) 
     181            e2v(jpi,:) = e2v(jpi-1,:) 
     182            e1f(jpi,:) = e1f(jpi-1,:) 
     183            e2f(jpi,:) = e2f(jpi-1,:) 
     184         ENDIF 
     185      ENDIF 
     186 
     187   END SUBROUTINE Agrif_InitValues_cont 
     188 
     189 
     190   SUBROUTINE agrif_declare_var() 
     191      !!---------------------------------------------------------------------- 
     192      !!                 *** ROUTINE Agrif_InitValues_cont *** 
     193      !! 
     194      !! ** Purpose :: Declaration of variables to be interpolated 
     195      !!---------------------------------------------------------------------- 
     196      USE par_oce 
     197      USE dom_oce 
     198      USE agrif_profiles 
     199      USE agrif_parameters 
     200 
     201      IMPLICIT NONE 
     202 
     203      INTEGER :: ind1, ind2, ind3 
     204      INTEGER :: nx, ny 
     205      INTEGER ::nbghostcellsfine_tot_x, nbghostcellsfine_tot_y 
     206      INTEGER :: irafx 
     207 
     208      EXTERNAL :: nemo_mapping 
     209 
     210      ! 1. Declaration of the type of variable which have to be interpolated 
     211      !--------------------------------------------------------------------- 
     212 
     213      nx=nlci ; ny=nlcj 
     214 
     215      ind2 = 2 + nbghostcells_x 
     216      ind3 = 2 + nbghostcells_y_s 
     217      nbghostcellsfine_tot_x=nbghostcells_x+1 
     218      nbghostcellsfine_tot_y=max(nbghostcells_y_s,nbghostcells_y_n)+1 
     219 
     220      irafx = Agrif_irhox() 
     221 
     222      ! In case of East-West periodicity, prevent AGRIF interpolation at east and west boundaries 
     223      ! The procnames will not be CALLed at these boundaries 
     224      if (jperio == 1) THEN 
     225        CALL Agrif_Set_NearCommonBorderX(.TRUE.) 
     226        CALL Agrif_Set_DistantCommonBorderX(.TRUE.) 
     227      endif 
     228      if (.not.lk_south) THEN 
     229        CALL Agrif_Set_NearCommonBorderY(.TRUE.) 
     230      endif 
     231 
     232      CALL agrif_declare_variable((/2,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),glamt_id) 
     233      CALL agrif_declare_variable((/1,2/),(/ind2-1,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),glamu_id) 
     234      CALL agrif_declare_variable((/2,1/),(/ind2,ind3-1/),(/'x','y'/),(/1,1/),(/nx,ny/),glamv_id) 
     235      CALL agrif_declare_variable((/1,1/),(/ind2-1,ind3-1/),(/'x','y'/),(/1,1/),(/nx,ny/),glamf_id) 
     236 
     237      CALL agrif_declare_variable((/2,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),gphit_id) 
     238      CALL agrif_declare_variable((/1,2/),(/ind2-1,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),gphiu_id) 
     239      CALL agrif_declare_variable((/2,1/),(/ind2,ind3-1/),(/'x','y'/),(/1,1/),(/nx,ny/),gphiv_id) 
     240      CALL agrif_declare_variable((/1,1/),(/ind2-1,ind3-1/),(/'x','y'/),(/1,1/),(/nx,ny/),gphif_id) 
     241 
     242      ! Horizontal scale factors 
     243 
     244      CALL agrif_declare_variable((/2,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),e1t_id) 
     245      CALL agrif_declare_variable((/1,2/),(/ind2-1,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),e1u_id) 
     246      CALL agrif_declare_variable((/2,1/),(/ind2,ind3-1/),(/'x','y'/),(/1,1/),(/nx,ny/),e1v_id) 
     247      CALL agrif_declare_variable((/1,1/),(/ind2-1,ind3-1/),(/'x','y'/),(/1,1/),(/nx,ny/),e1f_id) 
     248 
     249      CALL agrif_declare_variable((/2,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),e2t_id) 
     250      CALL agrif_declare_variable((/1,2/),(/ind2-1,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),e2u_id) 
     251      CALL agrif_declare_variable((/2,1/),(/ind2,ind3-1/),(/'x','y'/),(/1,1/),(/nx,ny/),e2v_id) 
     252      CALL agrif_declare_variable((/1,1/),(/ind2-1,ind3-1/),(/'x','y'/),(/1,1/),(/nx,ny/),e2f_id) 
     253 
     254      ! Bathymetry 
     255 
     256      CALL agrif_declare_variable((/2,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),bathy_id) 
     257 
     258      ! Vertical scale factors 
     259      CALL agrif_declare_variable((/2,2,0/),(/ind2,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/nx,ny,jpk/),e3t_id) 
     260      CALL agrif_declare_variable((/2,2,0/),(/ind2,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/nx,ny,jpk/),e3t_copy_id) 
     261      CALL agrif_declare_variable((/2,2,0/),(/ind2,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/nx,ny,jpk+1/),e3t_connect_id) 
     262 
     263      CALL agrif_declare_variable((/1,2,0/),(/ind2-1,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/nx,ny,jpk/),e3u_id) 
     264      CALL agrif_declare_variable((/2,1,0/),(/ind2,ind3-1,0/),(/'x','y','N'/),(/1,1,1/),(/nx,ny,jpk/),e3v_id) 
     265 
     266      ! Bottom level 
     267 
     268      CALL agrif_declare_variable((/2,2/),(/ind2,ind3/),(/'x','y'/),(/1,1/),(/nx,ny/),bottom_level_id) 
     269 
     270      CALL Agrif_Set_bcinterp(glamt_id,interp=AGRIF_linear) 
     271      CALL Agrif_Set_interp(glamt_id,interp=AGRIF_linear) 
     272      CALL Agrif_Set_bc( glamt_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 
     273 
     274      CALL Agrif_Set_bcinterp(glamu_id,interp=AGRIF_linear) 
     275      CALL Agrif_Set_interp(glamu_id,interp=AGRIF_linear) 
     276      CALL Agrif_Set_bc( glamu_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)/) ) 
     277 
     278      CALL Agrif_Set_bcinterp(glamv_id,interp=AGRIF_linear) 
     279      CALL Agrif_Set_interp(glamv_id,interp=AGRIF_linear) 
     280      CALL Agrif_Set_bc( glamv_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)/) ) 
     281 
     282      CALL Agrif_Set_bcinterp(glamf_id,interp=AGRIF_linear) 
     283      CALL Agrif_Set_interp(glamf_id,interp=AGRIF_linear) 
     284      CALL Agrif_Set_bc( glamf_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)/) ) 
     285 
     286      CALL Agrif_Set_bcinterp(gphit_id,interp=AGRIF_linear) 
     287      CALL Agrif_Set_interp(gphit_id,interp=AGRIF_linear) 
     288      CALL Agrif_Set_bc( gphit_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 
     289 
     290      CALL Agrif_Set_bcinterp(gphiu_id,interp=AGRIF_linear) 
     291      CALL Agrif_Set_interp(gphiu_id,interp=AGRIF_linear) 
     292      CALL Agrif_Set_bc( gphiu_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)/) ) 
     293 
     294      CALL Agrif_Set_bcinterp(gphiv_id,interp=AGRIF_linear) 
     295      CALL Agrif_Set_interp(gphiv_id,interp=AGRIF_linear) 
     296      CALL Agrif_Set_bc( gphiv_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)/) ) 
     297 
     298      CALL Agrif_Set_bcinterp(gphif_id,interp=AGRIF_linear) 
     299      CALL Agrif_Set_interp(gphif_id,interp=AGRIF_linear) 
     300      CALL Agrif_Set_bc( gphif_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)/) ) 
     301 
     302      ! 
     303 
     304      CALL Agrif_Set_bcinterp(e1t_id,interp=AGRIF_ppm) 
     305      CALL Agrif_Set_interp(e1t_id,interp=AGRIF_ppm) 
     306      CALL Agrif_Set_bc( e1t_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 
     307 
     308      CALL Agrif_Set_bcinterp(e1u_id, interp1=Agrif_linear, interp2=AGRIF_ppm) 
     309      CALL Agrif_Set_interp(e1u_id, interp1=Agrif_linear, interp2=AGRIF_ppm) 
     310      CALL Agrif_Set_bc( e1u_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)/) ) 
     311 
     312      CALL Agrif_Set_bcinterp(e1v_id,interp1=AGRIF_ppm, interp2=Agrif_linear) 
     313      CALL Agrif_Set_interp(e1v_id, interp1=AGRIF_ppm, interp2=Agrif_linear) 
     314      CALL Agrif_Set_bc( e1v_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)/) ) 
     315 
     316      CALL Agrif_Set_bcinterp(e1f_id,interp=AGRIF_linear) 
     317      CALL Agrif_Set_interp(e1f_id,interp=AGRIF_linear) 
     318      CALL Agrif_Set_bc( e1f_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)/) ) 
     319 
     320      CALL Agrif_Set_bcinterp(e2t_id,interp=AGRIF_ppm) 
     321      CALL Agrif_Set_interp(e2t_id,interp=AGRIF_ppm) 
     322      CALL Agrif_Set_bc( e2t_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 
     323 
     324      CALL Agrif_Set_bcinterp(e2u_id,interp1=Agrif_linear, interp2=AGRIF_ppm) 
     325      CALL Agrif_Set_interp(e2u_id,interp1=Agrif_linear, interp2=AGRIF_ppm) 
     326      CALL Agrif_Set_bc( e2u_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)/) ) 
     327 
     328      CALL Agrif_Set_bcinterp(e2v_id,interp1=AGRIF_ppm, interp2=Agrif_linear) 
     329      CALL Agrif_Set_interp(e2v_id,interp1=AGRIF_ppm, interp2=Agrif_linear) 
     330      CALL Agrif_Set_bc( e2v_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)/) ) 
     331 
     332      CALL Agrif_Set_bcinterp(e2f_id,interp=AGRIF_linear) 
     333      CALL Agrif_Set_interp(e2f_id,interp=AGRIF_linear) 
     334      CALL Agrif_Set_bc( e2f_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)/) ) 
     335 
     336      CALL Agrif_Set_bcinterp(bathy_id,interp=AGRIF_linear) 
     337      CALL Agrif_Set_interp(bathy_id,interp=AGRIF_linear) 
     338      CALL Agrif_Set_bc( bathy_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 
     339 
     340      ! Vertical scale factors 
     341      CALL Agrif_Set_bcinterp(e3t_id,interp=AGRIF_ppm) 
     342      CALL Agrif_Set_interp(e3t_id,interp=AGRIF_ppm) 
     343      CALL Agrif_Set_bc( e3t_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 
     344      CALL Agrif_Set_Updatetype( e3t_id, update = AGRIF_Update_Average) 
     345 
     346      CALL Agrif_Set_bcinterp(e3t_copy_id,interp=AGRIF_constant) 
     347      CALL Agrif_Set_interp(e3t_copy_id,interp=AGRIF_constant) 
     348      CALL Agrif_Set_bc( e3t_copy_id, (/-npt_copy*irafx,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/)) 
     349 
     350      CALL Agrif_Set_bcinterp(e3t_connect_id,interp=AGRIF_ppm) 
     351      CALL Agrif_Set_interp(e3t_connect_id,interp=AGRIF_ppm) 
     352      CALL Agrif_Set_bc( e3t_connect_id, (/-(npt_copy+npt_connect)*irafx-1,-npt_copy*irafx/)) 
     353 
     354      CALL Agrif_Set_bcinterp(e3u_id, interp1=Agrif_linear, interp2=AGRIF_ppm) 
     355      CALL Agrif_Set_interp(e3u_id, interp1=Agrif_linear, interp2=AGRIF_ppm) 
     356      CALL Agrif_Set_bc( e3u_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 
     357      CALL Agrif_Set_Updatetype(e3u_id,update1 = Agrif_Update_Copy, update2 = Agrif_Update_Average) 
     358 
     359      CALL Agrif_Set_bcinterp(e3v_id,interp1=AGRIF_ppm, interp2=Agrif_linear) 
     360      CALL Agrif_Set_interp(e3v_id, interp1=AGRIF_ppm, interp2=Agrif_linear) 
     361      CALL Agrif_Set_bc( e3v_id, (/0,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/) ) 
     362      CALL Agrif_Set_Updatetype(e3v_id,update1 = Agrif_Update_Average, update2 = Agrif_Update_Copy) 
     363 
     364      ! Bottom level 
     365      CALL Agrif_Set_bcinterp(bottom_level_id,interp=AGRIF_constant) 
     366      CALL Agrif_Set_interp(bottom_level_id,interp=AGRIF_constant) 
     367      CALL Agrif_Set_bc( bottom_level_id, (/-npt_copy*irafx,max(nbghostcellsfine_tot_x,nbghostcellsfine_tot_y)-1/)) 
     368      CALL Agrif_Set_Updatetype( bottom_level_id, update = AGRIF_Update_Max) 
     369 
     370      CALL Agrif_Set_ExternalMapping(nemo_mapping) 
     371 
     372   END SUBROUTINE agrif_declare_var 
     373 
     374   SUBROUTINE nemo_mapping(ndim,ptx,pty,bounds,bounds_chunks,correction_required,nb_chunks) 
     375      USE dom_oce 
     376      INTEGER :: ndim 
     377      INTEGER :: ptx, pty 
     378      INTEGER,DIMENSION(ndim,2,2) :: bounds 
     379      INTEGER,DIMENSION(:,:,:,:),allocatable :: bounds_chunks 
     380      LOGICAL,DIMENSION(:),allocatable :: correction_required 
     381      INTEGER :: nb_chunks 
     382      INTEGER :: i 
     383 
     384      IF (agrif_debug_interp) THEN 
     385         DO i = 1, ndim 
     386             print *,'direction = ',i,bounds(i,1,2),bounds(i,2,2) 
     387         END DO 
     388      ENDIF 
     389 
     390      IF( bounds(2,2,2) > jpjglo ) THEN 
     391         IF( bounds(2,1,2) <= jpjglo ) THEN 
     392            nb_chunks = 2 
     393            ALLOCATE(bounds_chunks(nb_chunks,ndim,2,2)) 
     394            ALLOCATE(correction_required(nb_chunks)) 
     395            DO i = 1, nb_chunks 
     396               bounds_chunks(i,:,:,:) = bounds 
     397            END DO 
     398 
     399         ! FIRST CHUNCK (for j<=jpjglo) 
     400 
     401            ! Original indices 
     402            bounds_chunks(1,1,1,1) = bounds(1,1,2) 
     403            bounds_chunks(1,1,2,1) = bounds(1,2,2) 
     404            bounds_chunks(1,2,1,1) = bounds(2,1,2) 
     405            bounds_chunks(1,2,2,1) = jpjglo 
     406 
     407            bounds_chunks(1,1,1,2) = bounds(1,1,2) 
     408            bounds_chunks(1,1,2,2) = bounds(1,2,2) 
     409            bounds_chunks(1,2,1,2) = bounds(2,1,2) 
     410            bounds_chunks(1,2,2,2) = jpjglo 
     411 
     412            ! Correction required or not 
     413            correction_required(1)=.FALSE. 
     414 
     415         ! SECOND CHUNCK (for j>jpjglo) 
     416 
     417            !Original indices 
     418            bounds_chunks(2,1,1,1) = bounds(1,1,2) 
     419            bounds_chunks(2,1,2,1) = bounds(1,2,2) 
     420            bounds_chunks(2,2,1,1) = jpjglo-2 
     421            bounds_chunks(2,2,2,1) = bounds(2,2,2) 
     422 
     423           ! Where to find them 
     424           ! We use the relation TAB(ji,jj)=TAB(jpiglo-ji+2,jpjglo-2-(jj-jpjglo)) 
     425 
     426            IF (ptx == 2) THEN ! T, V points 
     427               bounds_chunks(2,1,1,2) = jpiglo-bounds(1,2,2)+2 
     428               bounds_chunks(2,1,2,2) = jpiglo-bounds(1,1,2)+2 
     429            ELSE ! U, F points 
     430               bounds_chunks(2,1,1,2) = jpiglo-bounds(1,2,2)+1 
     431               bounds_chunks(2,1,2,2) = jpiglo-bounds(1,1,2)+1 
     432            ENDIF 
     433 
     434            IF (pty == 2) THEN ! T, U points 
     435               bounds_chunks(2,2,1,2) = jpjglo-2-(bounds(2,2,2) -jpjglo) 
     436               bounds_chunks(2,2,2,2) = jpjglo-2-(jpjglo-2      -jpjglo) 
     437            ELSE ! V, F points 
     438               bounds_chunks(2,2,1,2) = jpjglo-3-(bounds(2,2,2) -jpjglo) 
     439               bounds_chunks(2,2,2,2) = jpjglo-3-(jpjglo-2      -jpjglo) 
     440            ENDIF 
     441       
     442            ! Correction required or not 
     443            correction_required(2)=.TRUE. 
     444 
     445         ELSE 
     446            
     447            nb_chunks = 1 
     448            ALLOCATE(bounds_chunks(nb_chunks,ndim,2,2)) 
     449            ALLOCATE(correction_required(nb_chunks)) 
     450            DO i=1,nb_chunks 
     451                bounds_chunks(i,:,:,:) = bounds 
     452            END DO 
     453 
     454            bounds_chunks(1,1,1,1) = bounds(1,1,2) 
     455            bounds_chunks(1,1,2,1) = bounds(1,2,2) 
     456            bounds_chunks(1,2,1,1) = bounds(2,1,2) 
     457            bounds_chunks(1,2,2,1) = bounds(2,2,2) 
     458 
     459            bounds_chunks(1,1,1,2) = jpiglo-bounds(1,2,2)+2 
     460            bounds_chunks(1,1,2,2) = jpiglo-bounds(1,1,2)+2 
     461 
     462            bounds_chunks(1,2,1,2) = jpjglo-2-(bounds(2,2,2)-jpjglo) 
     463            bounds_chunks(1,2,2,2) = jpjglo-2-(bounds(2,1,2)-jpjglo) 
     464 
     465            IF (ptx == 2) THEN ! T, V points 
     466               bounds_chunks(1,1,1,2) = jpiglo-bounds(1,2,2)+2 
     467               bounds_chunks(1,1,2,2) = jpiglo-bounds(1,1,2)+2 
     468            ELSE ! U, F points 
     469               bounds_chunks(1,1,1,2) = jpiglo-bounds(1,2,2)+1 
     470               bounds_chunks(1,1,2,2) = jpiglo-bounds(1,1,2)+1 
     471            ENDIF 
     472 
     473            IF (pty == 2) THEN ! T, U points 
     474               bounds_chunks(1,2,1,2) = jpjglo-2-(bounds(2,2,2) -jpjglo) 
     475               bounds_chunks(1,2,2,2) = jpjglo-2-(bounds(2,1,2) -jpjglo) 
     476            ELSE ! V, F points 
     477               bounds_chunks(1,2,1,2) = jpjglo-3-(bounds(2,2,2) -jpjglo) 
     478               bounds_chunks(1,2,2,2) = jpjglo-3-(bounds(2,1,2) -jpjglo) 
     479            ENDIF 
     480 
     481            correction_required(1)=.TRUE. 
     482 
     483         ENDIF  ! bounds(2,1,2) <= jpjglo 
     484 
     485      ELSE IF (bounds(1,1,2) < 1) THEN 
     486          
     487         IF (bounds(1,2,2) > 0) THEN 
     488            nb_chunks = 2 
     489            ALLOCATE(correction_required(nb_chunks)) 
     490            correction_required=.FALSE. 
     491            ALLOCATE(bounds_chunks(nb_chunks,ndim,2,2)) 
     492            DO i=1,nb_chunks 
     493               bounds_chunks(i,:,:,:) = bounds 
     494            END DO 
     495 
     496            bounds_chunks(1,1,1,2) = bounds(1,1,2)+jpiglo-2 
     497            bounds_chunks(1,1,2,2) = 1+jpiglo-2 
     498 
     499            bounds_chunks(1,1,1,1) = bounds(1,1,2) 
     500            bounds_chunks(1,1,2,1) = 1 
     501 
     502            bounds_chunks(2,1,1,2) = 2 
     503            bounds_chunks(2,1,2,2) = bounds(1,2,2) 
     504 
     505            bounds_chunks(2,1,1,1) = 2 
     506            bounds_chunks(2,1,2,1) = bounds(1,2,2) 
     507         ELSE 
     508            nb_chunks = 1 
     509            ALLOCATE(correction_required(nb_chunks)) 
     510            correction_required=.FALSE. 
     511            ALLOCATE(bounds_chunks(nb_chunks,ndim,2,2)) 
     512            DO i=1,nb_chunks 
     513               bounds_chunks(i,:,:,:) = bounds 
     514            END DO 
     515            bounds_chunks(1,1,1,2) = bounds(1,1,2)+jpiglo-2 
     516            bounds_chunks(1,1,2,2) = bounds(1,2,2)+jpiglo-2 
     517 
     518            bounds_chunks(1,1,1,1) = bounds(1,1,2) 
     519            bounds_chunks(1,1,2,1) = bounds(1,2,2) 
     520         ENDIF 
     521       
     522      ELSE 
     523       
     524         nb_chunks=1 
     525         ALLOCATE(correction_required(nb_chunks)) 
     526         correction_required=.FALSE. 
     527         ALLOCATE(bounds_chunks(nb_chunks,ndim,2,2)) 
     528         DO i=1,nb_chunks 
     529            bounds_chunks(i,:,:,:) = bounds 
     530         END DO 
     531         bounds_chunks(1,1,1,2) = bounds(1,1,2) 
     532         bounds_chunks(1,1,2,2) = bounds(1,2,2) 
     533         bounds_chunks(1,2,1,2) = bounds(2,1,2) 
     534         bounds_chunks(1,2,2,2) = bounds(2,2,2) 
     535 
     536         bounds_chunks(1,1,1,1) = bounds(1,1,2) 
     537         bounds_chunks(1,1,2,1) = bounds(1,2,2) 
     538         bounds_chunks(1,2,1,1) = bounds(2,1,2) 
     539         bounds_chunks(1,2,2,1) = bounds(2,2,2) 
     540 
     541      ENDIF 
     542 
     543   END SUBROUTINE nemo_mapping 
     544 
     545   FUNCTION agrif_external_switch_index(ptx,pty,i1,isens) 
     546      USE dom_oce 
     547      INTEGER :: ptx, pty, i1, isens 
     548      INTEGER :: agrif_external_switch_index 
     549 
     550      IF( isens == 1 )  THEN 
     551         IF( ptx == 2 ) THEN ! T, V points 
     552            agrif_external_switch_index = jpiglo-i1+2 
     553         ELSE ! U, F points 
     554            agrif_external_switch_index = jpiglo-i1+1 
     555         ENDIF 
     556      ELSE IF (isens ==2) THEN 
     557         IF (pty == 2) THEN ! T, U points 
     558            agrif_external_switch_index = jpjglo-2-(i1 -jpjglo) 
     559         ELSE ! V, F points 
     560            agrif_external_switch_index = jpjglo-3-(i1 -jpjglo) 
     561         ENDIF 
     562      ENDIF 
     563 
     564   END FUNCTION agrif_external_switch_index 
     565 
     566   SUBROUTINE correct_field(tab2d,i1,i2,j1,j2) 
     567      USE dom_oce 
     568      INTEGER :: i1,i2,j1,j2 
     569      REAL,DIMENSION(i1:i2,j1:j2) :: tab2d 
     570 
     571      INTEGER :: i,j 
     572      REAL,DIMENSION(i1:i2,j1:j2) :: tab2dtemp 
     573 
     574      tab2dtemp = tab2d 
     575 
     576      DO j=j1,j2 
     577         DO i=i1,i2 
     578        tab2d(i,j)=tab2dtemp(i2-(i-i1),j2-(j-j1)) 
     579         END DO 
     580      ENDDO 
     581 
     582   END SUBROUTINE correct_field 
     583 
     584   SUBROUTINE agrif_init_lonlat() 
     585      USE agrif_profiles 
     586      USE agrif_util 
     587      USE dom_oce 
     588      
     589      EXTERNAL :: init_glamt, init_glamu, init_glamv, init_glamf 
     590      EXTERNAL :: init_gphit, init_gphiu, init_gphiv, init_gphif 
     591      EXTERNAL :: longitude_linear_interp 
     592 
     593      INTEGER :: ji,jj,i1,i2,j1,j2 
     594      REAL, DIMENSION(jpi,jpj) :: tab2dtemp 
     595      INTEGER :: ind2,ind3 
     596      INTEGER :: irhox, irhoy 
     597 
     598      irhox = agrif_irhox() 
     599      irhoy = agrif_irhoy() 
     600      CALL Agrif_Set_external_linear_interp(longitude_linear_interp) 
     601 
     602      CALL Agrif_Init_variable(glamt_id, procname = init_glamt) 
     603      CALL Agrif_Init_variable(glamu_id, procname = init_glamu) 
     604      CALL Agrif_Init_variable(glamv_id, procname = init_glamv) 
     605      CALL Agrif_Init_variable(glamf_id, procname = init_glamf) 
     606      CALL Agrif_UnSet_external_linear_interp() 
     607 
     608      CALL Agrif_Init_variable(gphit_id, procname = init_gphit) 
     609      CALL Agrif_Init_variable(gphiu_id, procname = init_gphiu) 
     610      CALL Agrif_Init_variable(gphiv_id, procname = init_gphiv) 
     611      CALL Agrif_Init_variable(gphif_id, procname = init_gphif) 
     612 
     613   END SUBROUTINE agrif_init_lonlat 
     614 
     615   REAL FUNCTION longitude_linear_interp(x1,x2,coeff) 
     616      REAL :: x1, x2, coeff 
     617      REAL :: val_interp 
     618 
     619      IF( (x1*x2 <= -50*50) ) THEN 
     620         IF( x1 < 0 ) THEN 
     621            val_interp = coeff *(x1+360.) + (1.-coeff) *x2 
     622         ELSE 
     623            val_interp = coeff *x1 + (1.-coeff) *(x2+360.) 
     624         ENDIF 
     625         IF ((val_interp) >=180.) val_interp = val_interp - 360. 
     626      ELSE 
     627         val_interp = coeff * x1 + (1.-coeff) * x2 
     628      ENDIF 
     629      longitude_linear_interp = val_interp 
     630 
     631   END FUNCTION longitude_linear_interp 
     632 
     633   SUBROUTINE agrif_init_scales() 
     634      USE agrif_profiles 
     635      USE agrif_util 
     636      USE dom_oce 
     637      USE lbclnk 
     638      LOGICAL :: ln_perio 
     639      INTEGER nx,ny 
     640 
     641      EXTERNAL :: init_e1t, init_e1u, init_e1v, init_e1f 
     642      EXTERNAL :: init_e2t, init_e2u, init_e2v, init_e2f 
     643 
     644      nx = nlci ; ny = nlcj 
     645 
     646      ln_perio=.FALSE. 
     647      if( jperio ==1 .OR. jperio==2 .OR. jperio==4) ln_perio=.TRUE. 
     648 
     649      CALL Agrif_Init_variable(e1t_id, procname = init_e1t) 
     650      CALL Agrif_Init_variable(e1u_id, procname = init_e1u) 
     651      CALL Agrif_Init_variable(e1v_id, procname = init_e1v) 
     652      CALL Agrif_Init_variable(e1f_id, procname = init_e1f) 
     653 
     654      CALL Agrif_Init_variable(e2t_id, procname = init_e2t) 
     655      CALL Agrif_Init_variable(e2u_id, procname = init_e2u) 
     656      CALL Agrif_Init_variable(e2v_id, procname = init_e2v) 
     657      CALL Agrif_Init_variable(e2f_id, procname = init_e2f) 
     658 
     659      CALL lbc_lnk( 'e1t', e1t, 'T', 1._wp) 
     660      CALL lbc_lnk( 'e2t', e2t, 'T', 1._wp) 
     661      CALL lbc_lnk( 'e1u', e1u, 'U', 1._wp) 
     662      CALL lbc_lnk( 'e2u', e2u, 'U', 1._wp) 
     663      CALL lbc_lnk( 'e1v', e1v, 'V', 1._wp) 
     664      CALL lbc_lnk( 'e2v', e2v, 'V', 1._wp) 
     665      CALL lbc_lnk( 'e1f', e1f, 'F', 1._wp) 
     666      CALL lbc_lnk( 'e2f', e2f, 'F', 1._wp) 
     667 
     668   END SUBROUTINE agrif_init_scales 
     669 
     670   SUBROUTINE init_glamt( ptab, i1, i2, j1, j2,  before, nb,ndir) 
     671      USE dom_oce 
    304672      !!---------------------------------------------------------------------- 
    305673      !!                  ***  ROUTINE interpsshn  *** 
    306       !!----------------------------------------------------------------------   
     674      !!---------------------------------------------------------------------- 
     675      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
     676      REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
     677      LOGICAL                         , INTENT(in   ) ::   before 
     678      INTEGER                         , INTENT(in   ) ::   nb , ndir 
     679 
     680      ! 
     681      !!---------------------------------------------------------------------- 
     682      ! 
     683      IF( before) THEN 
     684         ptab(i1:i2,j1:j2) = glamt(i1:i2,j1:j2) 
     685      ELSE 
     686          glamt(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) 
     687      ENDIF 
     688      ! 
     689   END SUBROUTINE init_glamt 
     690 
     691   SUBROUTINE init_glamu( ptab, i1, i2, j1, j2, before, nb,ndir) 
     692      USE dom_oce 
     693      !!---------------------------------------------------------------------- 
     694      !!                  ***  ROUTINE interpsshn  *** 
     695      !!---------------------------------------------------------------------- 
    307696      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
    308697      REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
     
    311700      LOGICAL  ::   western_side, eastern_side,northern_side,southern_side 
    312701      ! 
    313       !!----------------------------------------------------------------------   
    314       ! 
    315          western_side  = (nb == 1).AND.(ndir == 1) 
    316          eastern_side  = (nb == 1).AND.(ndir == 2) 
    317          southern_side = (nb == 2).AND.(ndir == 1) 
    318          northern_side = (nb == 2).AND.(ndir == 2) 
    319       IF( before) THEN 
    320          ptab(i1:i2,j1:j2) = glamt(i1:i2,j1:j2) 
    321       ELSE 
    322          glamt(i1:i2,j1:j2)=ptab 
    323       ENDIF 
    324       ! 
    325    END SUBROUTINE init_glamt 
    326  
    327     SUBROUTINE init_glamu( ptab, i1, i2, j1, j2, before, nb,ndir) 
    328     use dom_oce 
     702      !!---------------------------------------------------------------------- 
     703      ! 
     704      IF( before) THEN 
     705         ptab(i1:i2,j1:j2) = glamu(i1:i2,j1:j2) 
     706      ELSE 
     707          glamu(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) 
     708      ENDIF 
     709      ! 
     710   END SUBROUTINE init_glamu 
     711 
     712   SUBROUTINE init_glamv( ptab, i1, i2, j1, j2, before, nb,ndir) 
     713      USE dom_oce 
    329714      !!---------------------------------------------------------------------- 
    330715      !!                  ***  ROUTINE interpsshn  *** 
    331       !!----------------------------------------------------------------------   
    332       INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
    333       REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
    334       LOGICAL                         , INTENT(in   ) ::   before 
    335       INTEGER                         , INTENT(in   ) ::   nb , ndir 
    336       LOGICAL  ::   western_side, eastern_side,northern_side,southern_side 
    337       ! 
    338       !!----------------------------------------------------------------------   
    339       ! 
    340          western_side  = (nb == 1).AND.(ndir == 1) 
    341          eastern_side  = (nb == 1).AND.(ndir == 2) 
    342          southern_side = (nb == 2).AND.(ndir == 1) 
    343          northern_side = (nb == 2).AND.(ndir == 2) 
    344       IF( before) THEN 
    345          ptab(i1:i2,j1:j2) = glamu(i1:i2,j1:j2) 
    346       ELSE 
    347          glamu(i1:i2,j1:j2)=ptab 
    348       ENDIF 
    349       ! 
    350    END SUBROUTINE init_glamu 
    351  
    352    SUBROUTINE init_glamv( ptab, i1, i2, j1, j2, before, nb,ndir) 
    353    use dom_oce 
     716      !!---------------------------------------------------------------------- 
     717      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
     718      REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
     719      LOGICAL                         , INTENT(in   ) ::   before 
     720      INTEGER                         , INTENT(in   ) ::   nb , ndir 
     721      ! 
     722      !!---------------------------------------------------------------------- 
     723      ! 
     724      IF( before) THEN 
     725         ptab(i1:i2,j1:j2) = glamv(i1:i2,j1:j2) 
     726      ELSE 
     727          glamv(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) 
     728      ENDIF 
     729      ! 
     730   END SUBROUTINE init_glamv 
     731 
     732   SUBROUTINE init_glamf( ptab, i1, i2, j1, j2,  before, nb,ndir) 
     733      USE dom_oce 
     734      !!---------------------------------------------------------------------- 
     735      !!                  ***  ROUTINE init_glamf  *** 
     736      !!---------------------------------------------------------------------- 
     737      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
     738      REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
     739      LOGICAL                         , INTENT(in   ) ::   before 
     740      INTEGER                         , INTENT(in   ) ::   nb , ndir 
     741      ! 
     742      !!---------------------------------------------------------------------- 
     743      ! 
     744      IF( before) THEN 
     745         ptab(i1:i2,j1:j2) = glamf(i1:i2,j1:j2) 
     746      ELSE 
     747          glamf(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) 
     748      ENDIF 
     749      ! 
     750   END SUBROUTINE init_glamf 
     751 
     752   SUBROUTINE init_gphit( ptab, i1, i2, j1, j2, before, nb,ndir) 
     753      USE dom_oce 
     754      !!---------------------------------------------------------------------- 
     755      !!                  ***  ROUTINE init_gphit  *** 
     756      !!---------------------------------------------------------------------- 
     757      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
     758      REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
     759      LOGICAL                         , INTENT(in   ) ::   before 
     760      INTEGER                         , INTENT(in   ) ::   nb , ndir 
     761      ! 
     762      !!---------------------------------------------------------------------- 
     763      ! 
     764      IF( before) THEN 
     765         ptab(i1:i2,j1:j2) = gphit(i1:i2,j1:j2) 
     766      ELSE 
     767         gphit(i1:i2,j1:j2)=ptab 
     768      ENDIF 
     769      ! 
     770   END SUBROUTINE init_gphit 
     771 
     772   SUBROUTINE init_gphiu( ptab, i1, i2, j1, j2, before, nb,ndir) 
     773      USE dom_oce 
     774      !!---------------------------------------------------------------------- 
     775      !!                  ***  ROUTINE init_gphiu  *** 
     776      !!---------------------------------------------------------------------- 
     777      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
     778      REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
     779      LOGICAL                         , INTENT(in   ) ::   before 
     780      INTEGER                         , INTENT(in   ) ::   nb , ndir 
     781      ! 
     782      !!---------------------------------------------------------------------- 
     783      ! 
     784      IF( before) THEN 
     785         ptab(i1:i2,j1:j2) = gphiu(i1:i2,j1:j2) 
     786      ELSE 
     787         gphiu(i1:i2,j1:j2)=ptab 
     788      ENDIF 
     789      ! 
     790   END SUBROUTINE init_gphiu 
     791 
     792   SUBROUTINE init_gphiv( ptab, i1, i2, j1, j2, before, nb,ndir) 
     793      USE dom_oce 
     794      !!---------------------------------------------------------------------- 
     795      !!                  ***  ROUTINE init_gphiv  *** 
     796      !!---------------------------------------------------------------------- 
     797      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
     798      REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
     799      LOGICAL                         , INTENT(in   ) ::   before 
     800      INTEGER                         , INTENT(in   ) ::   nb , ndir 
     801      ! 
     802      !!---------------------------------------------------------------------- 
     803 
     804      IF( before) THEN 
     805         ptab(i1:i2,j1:j2) = gphiv(i1:i2,j1:j2) 
     806      ELSE 
     807         gphiv(i1:i2,j1:j2)=ptab 
     808      ENDIF 
     809      ! 
     810   END SUBROUTINE init_gphiv 
     811 
     812 
     813   SUBROUTINE init_gphif( ptab, i1, i2, j1, j2, before, nb,ndir) 
     814      USE dom_oce 
     815      !!---------------------------------------------------------------------- 
     816      !!                  ***  ROUTINE init_gphif  *** 
     817      !!---------------------------------------------------------------------- 
     818      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
     819      REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
     820      LOGICAL                         , INTENT(in   ) ::   before 
     821      INTEGER                         , INTENT(in   ) ::   nb , ndir 
     822      ! 
     823      !!---------------------------------------------------------------------- 
     824      ! 
     825      IF( before) THEN 
     826         ptab(i1:i2,j1:j2) = gphif(i1:i2,j1:j2) 
     827      ELSE 
     828         gphif(i1:i2,j1:j2)=ptab 
     829      ENDIF 
     830      ! 
     831   END SUBROUTINE init_gphif 
     832 
     833 
     834   SUBROUTINE init_e1t( ptab, i1, i2, j1, j2, before, nb,ndir) 
     835      USE dom_oce 
     836      USE agrif_parameters 
     837      !!---------------------------------------------------------------------- 
     838      !!                  ***  ROUTINE init_e1t  *** 
     839      !!---------------------------------------------------------------------- 
     840      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
     841      REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
     842      LOGICAL                         , INTENT(in   ) ::   before 
     843      INTEGER                         , INTENT(in   ) ::   nb , ndir 
     844      ! 
     845      !!---------------------------------------------------------------------- 
     846      ! 
     847      INTEGER :: jj 
     848 
     849      IF( before) THEN 
     850        ! May need to extend at south boundary 
     851          IF (j1<1) THEN 
     852            IF (.NOT.agrif_child(lk_south)) THEN 
     853              IF ((nbondj == -1).OR.(nbondj == 2)) THEN 
     854                DO jj=1,j2 
     855                  ptab(i1:i2,jj)=e1t(i1:i2,jj) 
     856                ENDDO 
     857                DO jj=j1,0 
     858                  ptab(i1:i2,jj)=e1t(i1:i2,1) 
     859                ENDDO 
     860              ENDIF 
     861            ELSE 
     862              stop "OUT OF BOUNDS" 
     863            ENDIF 
     864          ELSE 
     865             ptab(i1:i2,j1:j2) = e1t(i1:i2,j1:j2) 
     866          ENDIF 
     867      ELSE 
     868         e1t(i1:i2,j1:j2)=ptab/Agrif_rhoy() 
     869      ENDIF 
     870      ! 
     871   END SUBROUTINE init_e1t 
     872 
     873   SUBROUTINE init_e1u( ptab, i1, i2, j1, j2, before, nb,ndir) 
     874      USE dom_oce 
     875      USE agrif_parameters 
     876      !!---------------------------------------------------------------------- 
     877      !!                  ***  ROUTINE init_e1u  *** 
     878      !!---------------------------------------------------------------------- 
     879      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
     880      REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
     881      LOGICAL                         , INTENT(in   ) ::   before 
     882      INTEGER                         , INTENT(in   ) ::   nb , ndir 
     883      ! 
     884      !!---------------------------------------------------------------------- 
     885      ! 
     886      INTEGER :: jj 
     887 
     888      IF( before) THEN 
     889          IF (j1<1) THEN 
     890            IF (.NOT.agrif_child(lk_south)) THEN 
     891              IF ((nbondj == -1).OR.(nbondj == 2)) THEN 
     892                DO jj=1,j2 
     893                  ptab(i1:i2,jj)=e1u(i1:i2,jj) 
     894                ENDDO 
     895                DO jj=j1,0 
     896                  ptab(i1:i2,jj)=e1u(i1:i2,1) 
     897                ENDDO 
     898              ENDIF 
     899            ELSE 
     900              stop "OUT OF BOUNDS" 
     901            ENDIF 
     902          ELSE 
     903             ptab(i1:i2,j1:j2) = e1u(i1:i2,j1:j2) 
     904          ENDIF 
     905      ELSE 
     906         e1u(i1:i2,j1:j2)=ptab/Agrif_rhoy() 
     907      ENDIF 
     908      ! 
     909   END SUBROUTINE init_e1u 
     910 
     911   SUBROUTINE init_e1v( ptab, i1, i2, j1, j2, before, nb,ndir) 
     912      USE dom_oce 
     913      !!---------------------------------------------------------------------- 
     914      !!                  ***  ROUTINE init_e1v  *** 
     915      !!---------------------------------------------------------------------- 
     916      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
     917      REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
     918      LOGICAL                         , INTENT(in   ) ::   before 
     919      INTEGER                         , INTENT(in   ) ::   nb , ndir 
     920      ! 
     921      !!---------------------------------------------------------------------- 
     922      ! 
     923      IF( before) THEN 
     924         ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) 
     925      ELSE 
     926         e1v(i1:i2,j1:j2)=ptab/Agrif_rhoy() 
     927      ENDIF 
     928      ! 
     929   END SUBROUTINE init_e1v 
     930 
     931   SUBROUTINE init_e1f( ptab, i1, i2, j1, j2, before, nb,ndir) 
     932      USE dom_oce 
     933      !!---------------------------------------------------------------------- 
     934      !!                  ***  ROUTINE init_e1f  *** 
     935      !!---------------------------------------------------------------------- 
     936      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
     937      REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
     938      LOGICAL                         , INTENT(in   ) ::   before 
     939      INTEGER                         , INTENT(in   ) ::   nb , ndir 
     940      ! 
     941      !!---------------------------------------------------------------------- 
     942      ! 
     943      IF( before) THEN 
     944         ptab(i1:i2,j1:j2) = e1f(i1:i2,j1:j2) 
     945      ELSE 
     946         e1f(i1:i2,j1:j2)=ptab/Agrif_rhoy() 
     947      ENDIF 
     948      ! 
     949   END SUBROUTINE init_e1f 
     950 
     951   SUBROUTINE init_e2t( ptab, i1, i2, j1, j2, before, nb,ndir) 
     952      USE dom_oce 
     953      USE agrif_parameters 
     954      !!---------------------------------------------------------------------- 
     955      !!                  ***  ROUTINE init_e2t  *** 
     956      !!---------------------------------------------------------------------- 
     957      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
     958      REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
     959      LOGICAL                         , INTENT(in   ) ::   before 
     960      INTEGER                         , INTENT(in   ) ::   nb , ndir 
     961      ! 
     962      !!---------------------------------------------------------------------- 
     963      ! 
     964      INTEGER :: jj 
     965 
     966      IF( before) THEN 
     967          IF (j1<1) THEN 
     968            IF (.NOT.agrif_child(lk_south)) THEN 
     969              IF ((nbondj == -1).OR.(nbondj == 2)) THEN 
     970                DO jj=1,j2 
     971                  ptab(i1:i2,jj)=e2t(i1:i2,jj) 
     972                ENDDO 
     973                DO jj=j1,0 
     974                  ptab(i1:i2,jj)=e2t(i1:i2,1) 
     975                ENDDO 
     976              ENDIF 
     977            ELSE 
     978              stop "OUT OF BOUNDS" 
     979            ENDIF 
     980          ELSE 
     981             ptab(i1:i2,j1:j2) = e2t(i1:i2,j1:j2) 
     982          ENDIF 
     983      ELSE 
     984         e2t(i1:i2,j1:j2)=ptab/Agrif_rhoy() 
     985      ENDIF 
     986      ! 
     987   END SUBROUTINE init_e2t 
     988 
     989   SUBROUTINE init_e2u( ptab, i1, i2, j1, j2, before, nb,ndir) 
     990   USE dom_oce 
     991   USE agrif_parameters 
    354992      !!---------------------------------------------------------------------- 
    355993      !!                  ***  ROUTINE interpsshn  *** 
    356       !!----------------------------------------------------------------------   
    357       INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
    358       REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
    359       LOGICAL                         , INTENT(in   ) ::   before 
    360       INTEGER                         , INTENT(in   ) ::   nb , ndir 
    361       LOGICAL  ::   western_side, eastern_side,northern_side,southern_side 
    362       ! 
    363       !!----------------------------------------------------------------------   
    364       ! 
    365          western_side  = (nb == 1).AND.(ndir == 1) 
    366          eastern_side  = (nb == 1).AND.(ndir == 2) 
    367          southern_side = (nb == 2).AND.(ndir == 1) 
    368          northern_side = (nb == 2).AND.(ndir == 2) 
    369       IF( before) THEN 
    370          ptab(i1:i2,j1:j2) = glamv(i1:i2,j1:j2) 
    371       ELSE 
    372          glamv(i1:i2,j1:j2)=ptab 
    373       ENDIF 
    374       ! 
    375    END SUBROUTINE init_glamv 
    376  
    377    SUBROUTINE init_glamf( ptab, i1, i2, j1, j2, before, nb,ndir) 
    378    use dom_oce 
     994      !!---------------------------------------------------------------------- 
     995      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
     996      REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
     997      LOGICAL                         , INTENT(in   ) ::   before 
     998      INTEGER                         , INTENT(in   ) ::   nb , ndir 
     999      ! 
     1000      !!---------------------------------------------------------------------- 
     1001      ! 
     1002      INTEGER :: jj 
     1003 
     1004      IF( before) THEN 
     1005          IF (j1<1) THEN 
     1006            IF (.NOT.agrif_child(lk_south)) THEN 
     1007              IF ((nbondj == -1).OR.(nbondj == 2)) THEN 
     1008                DO jj=1,j2 
     1009                  ptab(i1:i2,jj)=e2u(i1:i2,jj) 
     1010                ENDDO 
     1011                DO jj=j1,0 
     1012                  ptab(i1:i2,jj)=e2u(i1:i2,1) 
     1013                ENDDO 
     1014              ENDIF 
     1015            ELSE 
     1016              stop "OUT OF BOUNDS" 
     1017            ENDIF 
     1018          ELSE 
     1019             ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) 
     1020          ENDIF 
     1021      ELSE 
     1022         e2u(i1:i2,j1:j2)=ptab/Agrif_rhoy() 
     1023      ENDIF 
     1024      ! 
     1025   END SUBROUTINE init_e2u 
     1026 
     1027   SUBROUTINE init_e2v( ptab, i1, i2, j1, j2, before, nb,ndir) 
     1028      USE dom_oce 
    3791029      !!---------------------------------------------------------------------- 
    3801030      !!                  ***  ROUTINE interpsshn  *** 
    381       !!----------------------------------------------------------------------   
    382       INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
    383       REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
    384       LOGICAL                         , INTENT(in   ) ::   before 
    385       INTEGER                         , INTENT(in   ) ::   nb , ndir 
    386       LOGICAL  ::   western_side, eastern_side,northern_side,southern_side 
    387       ! 
    388       !!----------------------------------------------------------------------   
    389       ! 
    390          western_side  = (nb == 1).AND.(ndir == 1) 
    391          eastern_side  = (nb == 1).AND.(ndir == 2) 
    392          southern_side = (nb == 2).AND.(ndir == 1) 
    393          northern_side = (nb == 2).AND.(ndir == 2) 
    394       IF( before) THEN 
    395          ptab(i1:i2,j1:j2) = glamf(i1:i2,j1:j2) 
    396       ELSE 
    397          glamf(i1:i2,j1:j2)=ptab 
    398       ENDIF 
    399       ! 
    400    END SUBROUTINE init_glamf 
    401  
    402    SUBROUTINE init_gphit( ptab, i1, i2, j1, j2, before, nb,ndir) 
    403    use dom_oce 
     1031      !!---------------------------------------------------------------------- 
     1032      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
     1033      REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
     1034      LOGICAL                         , INTENT(in   ) ::   before 
     1035      INTEGER                         , INTENT(in   ) ::   nb , ndir 
     1036      ! 
     1037      !!---------------------------------------------------------------------- 
     1038      IF( before) THEN 
     1039         ptab(i1:i2,j1:j2) = e2v(i1:i2,j1:j2) 
     1040      ELSE 
     1041         e2v(i1:i2,j1:j2)=ptab/Agrif_rhoy() 
     1042      ENDIF 
     1043      ! 
     1044   END SUBROUTINE init_e2v 
     1045 
     1046   SUBROUTINE init_e2f( ptab, i1, i2, j1, j2, before, nb,ndir) 
     1047   USE dom_oce 
    4041048      !!---------------------------------------------------------------------- 
    4051049      !!                  ***  ROUTINE interpsshn  *** 
    406       !!----------------------------------------------------------------------   
    407       INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
    408       REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
    409       LOGICAL                         , INTENT(in   ) ::   before 
    410       INTEGER                         , INTENT(in   ) ::   nb , ndir 
    411       LOGICAL  ::   western_side, eastern_side,northern_side,southern_side 
    412       ! 
    413       !!----------------------------------------------------------------------   
    414       ! 
    415          western_side  = (nb == 1).AND.(ndir == 1) 
    416          eastern_side  = (nb == 1).AND.(ndir == 2) 
    417          southern_side = (nb == 2).AND.(ndir == 1) 
    418          northern_side = (nb == 2).AND.(ndir == 2) 
    419       IF( before) THEN 
    420          ptab(i1:i2,j1:j2) = gphit(i1:i2,j1:j2) 
    421       ELSE 
    422          gphit(i1:i2,j1:j2)=ptab 
    423       ENDIF 
    424       ! 
    425    END SUBROUTINE init_gphit 
    426  
    427     SUBROUTINE init_gphiu( ptab, i1, i2, j1, j2, before, nb,ndir) 
    428     use dom_oce 
    429       !!---------------------------------------------------------------------- 
    430       !!                  ***  ROUTINE interpsshn  *** 
    431       !!----------------------------------------------------------------------   
    432       INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
    433       REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
    434       LOGICAL                         , INTENT(in   ) ::   before 
    435       INTEGER                         , INTENT(in   ) ::   nb , ndir 
    436       LOGICAL  ::   western_side, eastern_side,northern_side,southern_side 
    437       ! 
    438       !!----------------------------------------------------------------------   
    439       ! 
    440          western_side  = (nb == 1).AND.(ndir == 1) 
    441          eastern_side  = (nb == 1).AND.(ndir == 2) 
    442          southern_side = (nb == 2).AND.(ndir == 1) 
    443          northern_side = (nb == 2).AND.(ndir == 2) 
    444       IF( before) THEN 
    445          ptab(i1:i2,j1:j2) = gphiu(i1:i2,j1:j2) 
    446       ELSE 
    447          gphiu(i1:i2,j1:j2)=ptab 
    448       ENDIF 
    449       ! 
    450    END SUBROUTINE init_gphiu 
    451  
    452     SUBROUTINE init_gphiv( ptab, i1, i2, j1, j2, before, nb,ndir) 
    453     use dom_oce 
    454       !!---------------------------------------------------------------------- 
    455       !!                  ***  ROUTINE interpsshn  *** 
    456       !!----------------------------------------------------------------------   
    457       INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
    458       REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
    459       LOGICAL                         , INTENT(in   ) ::   before 
    460       INTEGER                         , INTENT(in   ) ::   nb , ndir 
    461       LOGICAL  ::   western_side, eastern_side,northern_side,southern_side 
    462       ! 
    463       !!----------------------------------------------------------------------   
    464       ! 
    465          western_side  = (nb == 1).AND.(ndir == 1) 
    466          eastern_side  = (nb == 1).AND.(ndir == 2) 
    467          southern_side = (nb == 2).AND.(ndir == 1) 
    468          northern_side = (nb == 2).AND.(ndir == 2) 
    469       IF( before) THEN 
    470          ptab(i1:i2,j1:j2) = gphiv(i1:i2,j1:j2) 
    471       ELSE 
    472          gphiv(i1:i2,j1:j2)=ptab 
    473       ENDIF 
    474       ! 
    475    END SUBROUTINE init_gphiv 
    476  
    477  
    478    SUBROUTINE init_gphif( ptab, i1, i2, j1, j2, before, nb,ndir) 
    479    use dom_oce 
    480       !!---------------------------------------------------------------------- 
    481       !!                  ***  ROUTINE interpsshn  *** 
    482       !!----------------------------------------------------------------------   
    483       INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
    484       REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
    485       LOGICAL                         , INTENT(in   ) ::   before 
    486       INTEGER                         , INTENT(in   ) ::   nb , ndir 
    487       LOGICAL  ::   western_side, eastern_side,northern_side,southern_side 
    488       ! 
    489       !!----------------------------------------------------------------------   
    490       ! 
    491          western_side  = (nb == 1).AND.(ndir == 1) 
    492          eastern_side  = (nb == 1).AND.(ndir == 2) 
    493          southern_side = (nb == 2).AND.(ndir == 1) 
    494          northern_side = (nb == 2).AND.(ndir == 2) 
    495       IF( before) THEN 
    496          ptab(i1:i2,j1:j2) = gphif(i1:i2,j1:j2) 
    497       ELSE 
    498          gphif(i1:i2,j1:j2)=ptab 
    499       ENDIF 
    500       ! 
    501    END SUBROUTINE init_gphif 
    502  
    503  
    504    SUBROUTINE init_e1t( ptab, i1, i2, j1, j2, before, nb,ndir) 
    505    use dom_oce 
    506       !!---------------------------------------------------------------------- 
    507       !!                  ***  ROUTINE interpsshn  *** 
    508       !!----------------------------------------------------------------------   
    509       INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
    510       REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
    511       LOGICAL                         , INTENT(in   ) ::   before 
    512       INTEGER                         , INTENT(in   ) ::   nb , ndir 
    513       LOGICAL  ::   western_side, eastern_side,northern_side,southern_side 
    514       ! 
    515       !!----------------------------------------------------------------------   
    516       ! 
    517          western_side  = (nb == 1).AND.(ndir == 1) 
    518          eastern_side  = (nb == 1).AND.(ndir == 2) 
    519          southern_side = (nb == 2).AND.(ndir == 1) 
    520          northern_side = (nb == 2).AND.(ndir == 2) 
    521       IF( before) THEN 
    522          ptab(i1:i2,j1:j2) = e1t(i1:i2,j1:j2) 
    523       ELSE 
    524          e1t(i1:i2,j1:j2)=ptab/Agrif_rhoy() 
    525       ENDIF 
    526       ! 
    527    END SUBROUTINE init_e1t 
    528  
    529    SUBROUTINE init_e1u( ptab, i1, i2, j1, j2, before, nb,ndir) 
    530    use dom_oce 
    531       !!---------------------------------------------------------------------- 
    532       !!                  ***  ROUTINE interpsshn  *** 
    533       !!----------------------------------------------------------------------   
    534       INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
    535       REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
    536       LOGICAL                         , INTENT(in   ) ::   before 
    537       INTEGER                         , INTENT(in   ) ::   nb , ndir 
    538       LOGICAL  ::   western_side, eastern_side,northern_side,southern_side 
    539       ! 
    540       !!----------------------------------------------------------------------   
    541       ! 
    542          western_side  = (nb == 1).AND.(ndir == 1) 
    543          eastern_side  = (nb == 1).AND.(ndir == 2) 
    544          southern_side = (nb == 2).AND.(ndir == 1) 
    545          northern_side = (nb == 2).AND.(ndir == 2) 
    546       IF( before) THEN 
    547          ptab(i1:i2,j1:j2) = e1u(i1:i2,j1:j2) 
    548       ELSE 
    549          e1u(i1:i2,j1:j2)=ptab/Agrif_rhoy() 
    550       ENDIF 
    551       ! 
    552    END SUBROUTINE init_e1u 
    553  
    554    SUBROUTINE init_e1v( ptab, i1, i2, j1, j2, before, nb,ndir) 
    555    use dom_oce 
    556       !!---------------------------------------------------------------------- 
    557       !!                  ***  ROUTINE interpsshn  *** 
    558       !!----------------------------------------------------------------------   
    559       INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
    560       REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
    561       LOGICAL                         , INTENT(in   ) ::   before 
    562       INTEGER                         , INTENT(in   ) ::   nb , ndir 
    563       LOGICAL  ::   western_side, eastern_side,northern_side,southern_side 
    564       ! 
    565       !!----------------------------------------------------------------------   
    566       ! 
    567          western_side  = (nb == 1).AND.(ndir == 1) 
    568          eastern_side  = (nb == 1).AND.(ndir == 2) 
    569          southern_side = (nb == 2).AND.(ndir == 1) 
    570          northern_side = (nb == 2).AND.(ndir == 2) 
    571       IF( before) THEN 
    572          ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) 
    573       ELSE 
    574          e1v(i1:i2,j1:j2)=ptab/Agrif_rhoy() 
    575       ENDIF 
    576       ! 
    577    END SUBROUTINE init_e1v 
    578  
    579    SUBROUTINE init_e1f( ptab, i1, i2, j1, j2, before, nb,ndir) 
    580    use dom_oce 
    581       !!---------------------------------------------------------------------- 
    582       !!                  ***  ROUTINE interpsshn  *** 
    583       !!----------------------------------------------------------------------   
    584       INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
    585       REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
    586       LOGICAL                         , INTENT(in   ) ::   before 
    587       INTEGER                         , INTENT(in   ) ::   nb , ndir 
    588       LOGICAL  ::   western_side, eastern_side,northern_side,southern_side 
    589       ! 
    590       !!----------------------------------------------------------------------   
    591       ! 
    592          western_side  = (nb == 1).AND.(ndir == 1) 
    593          eastern_side  = (nb == 1).AND.(ndir == 2) 
    594          southern_side = (nb == 2).AND.(ndir == 1) 
    595          northern_side = (nb == 2).AND.(ndir == 2) 
    596       IF( before) THEN 
    597          ptab(i1:i2,j1:j2) = e1f(i1:i2,j1:j2) 
    598       ELSE 
    599          e1f(i1:i2,j1:j2)=ptab/Agrif_rhoy() 
    600       ENDIF 
    601       ! 
    602    END SUBROUTINE init_e1f 
    603  
    604   SUBROUTINE init_e2t( ptab, i1, i2, j1, j2, before, nb,ndir) 
    605    use dom_oce 
    606       !!---------------------------------------------------------------------- 
    607       !!                  ***  ROUTINE interpsshn  *** 
    608       !!----------------------------------------------------------------------   
    609       INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
    610       REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
    611       LOGICAL                         , INTENT(in   ) ::   before 
    612       INTEGER                         , INTENT(in   ) ::   nb , ndir 
    613       LOGICAL  ::   western_side, eastern_side,northern_side,southern_side 
    614       ! 
    615       !!----------------------------------------------------------------------   
    616       ! 
    617          western_side  = (nb == 1).AND.(ndir == 1) 
    618          eastern_side  = (nb == 1).AND.(ndir == 2) 
    619          southern_side = (nb == 2).AND.(ndir == 1) 
    620          northern_side = (nb == 2).AND.(ndir == 2) 
    621       IF( before) THEN 
    622          ptab(i1:i2,j1:j2) = e2t(i1:i2,j1:j2) 
    623       ELSE 
    624          e2t(i1:i2,j1:j2)=ptab/Agrif_rhoy() 
    625       ENDIF 
    626       ! 
    627    END SUBROUTINE init_e2t 
    628  
    629    SUBROUTINE init_e2u( ptab, i1, i2, j1, j2, before, nb,ndir) 
    630    use dom_oce 
    631       !!---------------------------------------------------------------------- 
    632       !!                  ***  ROUTINE interpsshn  *** 
    633       !!----------------------------------------------------------------------   
    634       INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
    635       REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
    636       LOGICAL                         , INTENT(in   ) ::   before 
    637       INTEGER                         , INTENT(in   ) ::   nb , ndir 
    638       LOGICAL  ::   western_side, eastern_side,northern_side,southern_side 
    639       ! 
    640       !!----------------------------------------------------------------------   
    641       ! 
    642          western_side  = (nb == 1).AND.(ndir == 1) 
    643          eastern_side  = (nb == 1).AND.(ndir == 2) 
    644          southern_side = (nb == 2).AND.(ndir == 1) 
    645          northern_side = (nb == 2).AND.(ndir == 2) 
    646       IF( before) THEN 
    647          ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) 
    648       ELSE 
    649          e2u(i1:i2,j1:j2)=ptab/Agrif_rhoy() 
    650       ENDIF 
    651       ! 
    652    END SUBROUTINE init_e2u 
    653  
    654    SUBROUTINE init_e2v( ptab, i1, i2, j1, j2, before, nb,ndir) 
    655    use dom_oce 
    656       !!---------------------------------------------------------------------- 
    657       !!                  ***  ROUTINE interpsshn  *** 
    658       !!----------------------------------------------------------------------   
    659       INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
    660       REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
    661       LOGICAL                         , INTENT(in   ) ::   before 
    662       INTEGER                         , INTENT(in   ) ::   nb , ndir 
    663       LOGICAL  ::   western_side, eastern_side,northern_side,southern_side 
    664       ! 
    665       !!----------------------------------------------------------------------   
    666       ! 
    667          western_side  = (nb == 1).AND.(ndir == 1) 
    668          eastern_side  = (nb == 1).AND.(ndir == 2) 
    669          southern_side = (nb == 2).AND.(ndir == 1) 
    670          northern_side = (nb == 2).AND.(ndir == 2) 
    671       IF( before) THEN 
    672          ptab(i1:i2,j1:j2) = e2v(i1:i2,j1:j2) 
    673       ELSE 
    674          e2v(i1:i2,j1:j2)=ptab/Agrif_rhoy() 
    675       ENDIF 
    676       ! 
    677    END SUBROUTINE init_e2v 
    678  
    679    SUBROUTINE init_e2f( ptab, i1, i2, j1, j2, before, nb,ndir) 
    680    use dom_oce 
    681       !!---------------------------------------------------------------------- 
    682       !!                  ***  ROUTINE interpsshn  *** 
    683       !!----------------------------------------------------------------------   
    684       INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
    685       REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
    686       LOGICAL                         , INTENT(in   ) ::   before 
    687       INTEGER                         , INTENT(in   ) ::   nb , ndir 
    688       LOGICAL  ::   western_side, eastern_side,northern_side,southern_side 
    689       ! 
    690       !!----------------------------------------------------------------------   
    691       ! 
    692          western_side  = (nb == 1).AND.(ndir == 1) 
    693          eastern_side  = (nb == 1).AND.(ndir == 2) 
    694          southern_side = (nb == 2).AND.(ndir == 1) 
    695          northern_side = (nb == 2).AND.(ndir == 2) 
     1050      !!---------------------------------------------------------------------- 
     1051      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
     1052      REAL, DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab 
     1053      LOGICAL                         , INTENT(in   ) ::   before 
     1054      INTEGER                         , INTENT(in   ) ::   nb , ndir 
     1055      ! 
     1056      !!---------------------------------------------------------------------- 
     1057      ! 
    6961058      IF( before) THEN 
    6971059         ptab(i1:i2,j1:j2) = e2f(i1:i2,j1:j2) 
     
    7031065 
    7041066 
    705 SUBROUTINE agrif_nemo_init 
    706 USE agrif_parameters 
    707 USE in_out_manager 
    708 USE lib_mpp 
    709  
    710     
    711    !! 
    712    IMPLICIT NONE 
    713     
    714    INTEGER ::   ios 
    715     
    716    NAMELIST/namagrif/ nn_cln_update,ln_spc_dyn,rn_sponge_tra,rn_sponge_dyn,ln_chk_bathy,npt_connect, npt_copy 
     1067   SUBROUTINE agrif_nemo_init 
     1068      USE agrif_parameters 
     1069      USE dom_oce 
     1070      USE in_out_manager 
     1071      USE lib_mpp 
     1072      !! 
     1073      IMPLICIT NONE 
     1074 
     1075      INTEGER ::   ios 
     1076 
     1077      NAMELIST/namagrif/ nn_cln_update,ln_spc_dyn,rn_sponge_tra,rn_sponge_dyn,ln_chk_bathy,npt_connect,   & 
     1078      & npt_copy 
    7171079 
    7181080      REWIND( numnam_ref )              ! Namelist namagrif in reference namelist : nesting parameters 
     
    7301092         WRITE(numout,*) '~~~~~~~' 
    7311093         WRITE(numout,*) '   Namelist namagrif : set nesting parameters' 
    732          WRITE(numout,*) '      npt_copy     = ', npt_copy 
    733          WRITE(numout,*) '      npt_connect  = ', npt_connect 
    734       ENDIF 
    735        
    736 END SUBROUTINE agrif_nemo_init 
    737     
    738  
    739 SUBROUTINE Agrif_detect( kg, ksizex ) 
     1094         WRITE(numout,*) '      npt_copy             = ', npt_copy 
     1095         WRITE(numout,*) '      npt_connect          = ', npt_connect 
     1096      ENDIF 
     1097 
     1098   ! Set the number of ghost cells according to periodicity 
     1099 
     1100      nbghostcells_x = nbghostcells 
     1101      nbghostcells_y_s = nbghostcells 
     1102      nbghostcells_y_n = nbghostcells 
     1103 
     1104      lk_west  = .NOT. ( Agrif_Ix() == 1 ) 
     1105      lk_east  = .NOT. ( Agrif_Ix() + nbcellsx/AGRIF_Irhox() == Agrif_Parent(jpiglo) -1 ) 
     1106      lk_south = .NOT. ( Agrif_Iy() == 1 ) 
     1107      lk_north = .NOT. ( Agrif_Iy() + nbcellsy/AGRIF_Irhoy() == Agrif_Parent(jpjglo) -1 ) 
     1108 
     1109      IF (.not.agrif_root()) THEN 
     1110        IF (jperio == 1) THEN 
     1111          nbghostcells_x = 0 
     1112        ENDIF 
     1113        IF (.NOT.lk_south) THEN 
     1114          nbghostcells_y_s = 0 
     1115        ENDIF 
     1116      ENDIF 
     1117 
     1118   END SUBROUTINE agrif_nemo_init 
     1119 
     1120   SUBROUTINE Agrif_detect( kg, ksizex ) 
    7401121      !!---------------------------------------------------------------------- 
    7411122      !!                      *** ROUTINE Agrif_detect *** 
    7421123      !!---------------------------------------------------------------------- 
    743    INTEGER, DIMENSION(2) :: ksizex 
    744    INTEGER, DIMENSION(ksizex(1),ksizex(2)) :: kg  
    745       !!---------------------------------------------------------------------- 
    746    ! 
    747    RETURN 
    748    ! 
    749 END SUBROUTINE Agrif_detect 
    750 SUBROUTINE agrif_before_regridding 
    751 END SUBROUTINE agrif_before_regridding 
    752  
    753 SUBROUTINE Agrif_InvLoc( indloc, nprocloc, i, indglob ) 
    754       !!---------------------------------------------------------------------- 
    755       !!                     *** ROUTINE Agrif_InvLoc *** 
    756       !!---------------------------------------------------------------------- 
    757    USE dom_oce 
    758    !! 
    759    IMPLICIT NONE 
    760    ! 
    761    INTEGER :: indglob, indloc, nprocloc, i 
    762       !!---------------------------------------------------------------------- 
    763    ! 
    764    SELECT CASE( i ) 
    765    CASE(1)   ;   indglob = indloc + nimppt(nprocloc+1) - 1 
    766    CASE(2)   ;   indglob = indloc + njmppt(nprocloc+1) - 1 
    767    CASE DEFAULT 
    768       indglob = indloc 
    769    END SELECT 
    770    ! 
    771 END SUBROUTINE Agrif_InvLoc 
    772  
    773 SUBROUTINE Agrif_get_proc_info( imin, imax, jmin, jmax ) 
     1124      INTEGER, DIMENSION(2) :: ksizex 
     1125      INTEGER, DIMENSION(ksizex(1),ksizex(2)) :: kg 
     1126      !!---------------------------------------------------------------------- 
     1127      ! 
     1128      RETURN 
     1129      ! 
     1130   END SUBROUTINE Agrif_detect 
     1131 
     1132   SUBROUTINE agrif_before_regridding 
     1133   END SUBROUTINE agrif_before_regridding 
     1134 
     1135# if defined  key_mpp_mpi 
     1136   SUBROUTINE Agrif_InvLoc( indloc, nprocloc, i, indglob ) 
     1137         !!---------------------------------------------------------------------- 
     1138         !!                     *** ROUTINE Agrif_InvLoc *** 
     1139         !!---------------------------------------------------------------------- 
     1140      USE dom_oce 
     1141      !! 
     1142      IMPLICIT NONE 
     1143      ! 
     1144      INTEGER :: indglob, indloc, nprocloc, i 
     1145         !!---------------------------------------------------------------------- 
     1146      ! 
     1147      SELECT CASE( i ) 
     1148      CASE(1)   ;   indglob = indloc + nimppt(nprocloc+1) - 1 
     1149      CASE(2)   ;   indglob = indloc + njmppt(nprocloc+1) - 1 
     1150      CASE DEFAULT 
     1151         indglob = indloc 
     1152      END SELECT 
     1153      ! 
     1154   END SUBROUTINE Agrif_InvLoc 
     1155 
     1156   SUBROUTINE Agrif_get_proc_info( imin, imax, jmin, jmax ) 
    7741157      !!---------------------------------------------------------------------- 
    7751158      !!                 *** ROUTINE Agrif_get_proc_info *** 
    7761159      !!---------------------------------------------------------------------- 
    777    USE par_oce 
    778    USE dom_oce  
    779    !! 
    780    IMPLICIT NONE 
    781    ! 
    782    INTEGER, INTENT(out) :: imin, imax 
    783    INTEGER, INTENT(out) :: jmin, jmax 
    784       !!---------------------------------------------------------------------- 
    785    ! 
    786    imin = nimppt(Agrif_Procrank+1)  ! ????? 
    787    jmin = njmppt(Agrif_Procrank+1)  ! ????? 
    788    imax = imin + jpi - 1 
    789    jmax = jmin + jpj - 1 
    790    !  
    791 END SUBROUTINE Agrif_get_proc_info 
    792  
    793 SUBROUTINE Agrif_estimate_parallel_cost(imin, imax,jmin, jmax, nbprocs, grid_cost) 
     1160      USE par_oce 
     1161      USE dom_oce 
     1162      !! 
     1163      IMPLICIT NONE 
     1164      ! 
     1165      INTEGER, INTENT(out) :: imin, imax 
     1166      INTEGER, INTENT(out) :: jmin, jmax 
     1167         !!---------------------------------------------------------------------- 
     1168      ! 
     1169      imin = nimppt(Agrif_Procrank+1)  ! ????? 
     1170      jmin = njmppt(Agrif_Procrank+1)  ! ????? 
     1171      imax = imin + jpi - 1 
     1172      jmax = jmin + jpj - 1 
     1173      ! 
     1174   END SUBROUTINE Agrif_get_proc_info 
     1175 
     1176   SUBROUTINE Agrif_estimate_parallel_cost(imin, imax,jmin, jmax, nbprocs, grid_cost) 
    7941177      !!---------------------------------------------------------------------- 
    7951178      !!                 *** ROUTINE Agrif_estimate_parallel_cost *** 
    7961179      !!---------------------------------------------------------------------- 
    797    USE par_oce 
    798    !! 
    799    IMPLICIT NONE 
    800    ! 
    801    INTEGER,  INTENT(in)  :: imin, imax 
    802    INTEGER,  INTENT(in)  :: jmin, jmax 
    803    INTEGER,  INTENT(in)  :: nbprocs 
    804    REAL(wp), INTENT(out) :: grid_cost 
    805       !!---------------------------------------------------------------------- 
    806    ! 
    807    grid_cost = REAL(imax-imin+1,wp)*REAL(jmax-jmin+1,wp) / REAL(nbprocs,wp) 
    808    ! 
    809 END SUBROUTINE Agrif_estimate_parallel_cost 
    810  
     1180      USE par_oce 
     1181      !! 
     1182      IMPLICIT NONE 
     1183      ! 
     1184      INTEGER,  INTENT(in)  :: imin, imax 
     1185      INTEGER,  INTENT(in)  :: jmin, jmax 
     1186      INTEGER,  INTENT(in)  :: nbprocs 
     1187      REAL(wp), INTENT(out) :: grid_cost 
     1188         !!---------------------------------------------------------------------- 
     1189      ! 
     1190      grid_cost = REAL(imax-imin+1,wp)*REAL(jmax-jmin+1,wp) / REAL(nbprocs,wp) 
     1191      ! 
     1192   END SUBROUTINE Agrif_estimate_parallel_cost 
     1193# endif 
    8111194#endif 
  • utils/tools/DOMAINcfg/src/bilinear_interp.f90

    r12414 r13204  
     1! 
    12MODULE bilinear_interp 
    2  
    3 END MODULE 
     3  ! 
     4  USE agrif_modutil 
     5  ! 
     6  !************************************************************************ 
     7  !                           * 
     8  ! MODULE  BILINEAR INTERP                  * 
     9  !                           * 
     10  ! bilinear interpolation routines from SCRIP package         *      
     11  !                           * 
     12  !http://climate.lanl.gov/Software/SCRIP/            * 
     13  !                           * 
     14  !Bilinear remapping                     * 
     15  !                           * 
     16  !************************************************************************ 
     17  !      
     18  !----------------------------------------------------------------------- 
     19  IMPLICIT NONE 
     20 
     21  !----------------------------------------------------------------------- 
     22  !     variables that describe each grid 
     23  !----------------------------------------------------------------------- 
     24  ! 
     25  INTEGER :: grid1_size,grid2_size,grid1_rank, grid2_rank 
     26  !       
     27  INTEGER, DIMENSION(:), ALLOCATABLE :: grid1_dims, grid2_dims   
     28  !   
     29 
     30  !----------------------------------------------------------------------- 
     31  !     grid coordinates and masks 
     32  !----------------------------------------------------------------------- 
     33  ! 
     34  LOGICAL, DIMENSION(:), ALLOCATABLE :: grid1_mask,grid2_mask         
     35  ! each grid center in radians 
     36  REAL*8,DIMENSION(:), ALLOCATABLE :: & 
     37       grid1_center_lat,  & 
     38       grid1_center_lon,  &  
     39       grid2_center_lat,  & 
     40       grid2_center_lon,  & 
     41       grid1_frac,        & ! fractional area of grid cells 
     42       grid2_frac           ! participating in remapping 
     43  ! 
     44  ! lat/lon bounding box for use in restricting grid searches 
     45  ! 
     46  REAL*8,DIMENSION(:,:), ALLOCATABLE :: grid1_bound_box,grid2_bound_box    
     47  ! 
     48  !----------------------------------------------------------------------- 
     49  !     bins for restricting searches 
     50  !----------------------------------------------------------------------- 
     51  ! 
     52  ! num of bins for restricted srch 
     53  INTEGER, PARAMETER :: num_srch_bins = 90   
     54  ! 
     55  ! min,max adds for grid cells in this lat bin 
     56  ! 
     57  INTEGER,DIMENSION(:,:), ALLOCATABLE :: bin_addr1,bin_addr2  
     58  ! 
     59  ! min,max longitude for each search bin 
     60  ! 
     61  REAL*8, DIMENSION(:,:), ALLOCATABLE  :: bin_lats,bin_lons  
     62 
     63  REAL*8, PARAMETER :: zero   = 0.0,  & 
     64       one    = 1.0,  & 
     65       two    = 2.0,  & 
     66       three  = 3.0,  & 
     67       four   = 4.0,  & 
     68       five   = 5.0,  &  
     69       half   = 0.5,  & 
     70       quart  = 0.25, & 
     71       bignum = 1.e+20, & 
     72       tiny   = 1.e-14, & 
     73       pi     = 3.14159265359, & 
     74       pi2    = two*pi, & 
     75       pih    = half*pi         
     76 
     77  REAL*8, PARAMETER :: deg2rad = pi/180. 
     78  !  
     79  ! max iteration count for i,j iteration  
     80  !     
     81  INTEGER , PARAMETER :: max_iter = 100    
     82  ! 
     83  ! convergence criterion 
     84  ! 
     85  REAL*8, PARAMETER :: converge = 1.e-10 
     86 
     87 
     88  !   
     89  INTEGER, PARAMETER :: norm_opt_none    = 1 & 
     90       ,norm_opt_dstarea = 2 & 
     91       ,norm_opt_frcarea = 3 
     92  ! 
     93  INTEGER, PARAMETER :: map_type_conserv  = 1 & 
     94       ,map_type_bilinear = 2 & 
     95       ,map_type_bicubic  = 3 & 
     96       ,map_type_distwgt  = 4 
     97  ! 
     98  INTEGER :: max_links_map1  &  ! current size of link arrays 
     99       ,num_links_map1  &  ! actual number of links for remapping 
     100       ,max_links_map2  &  ! current size of link arrays 
     101       ,num_links_map2  &  ! actual number of links for remapping 
     102       ,num_maps        &  ! num of remappings for this grid pair 
     103       ,num_wts         &  ! num of weights used in remapping 
     104       ,map_type        &  ! identifier for remapping method 
     105       ,norm_opt        &  ! option for normalization (conserv only) 
     106       ,resize_increment ! default amount to increase array size 
     107 
     108  INTEGER , DIMENSION(:), ALLOCATABLE :: & 
     109       grid1_add_map1, &  ! grid1 address for each link in mapping 1 
     110       grid2_add_map1, &  ! grid2 address for each link in mapping 1 
     111       grid1_add_map2, &  ! grid1 address for each link in mapping 2 
     112       grid2_add_map2    ! grid2 address for each link in mapping 2 
     113 
     114  REAL*8, DIMENSION(:,:), ALLOCATABLE ::   & 
     115       wts_map1, &   ! map weights for each link (num_wts,max_links) 
     116       wts_map2     ! map weights for each link (num_wts,max_links) 
     117  ! 
     118CONTAINS  
     119  !      
     120!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     121  !************************************************************************ 
     122  !   SUBROUTINE GRID_INIT 
     123  !************************************************************************ 
     124!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     125  ! 
     126  SUBROUTINE get_remap_matrix(grid1_lat,grid2_lat,grid1_lon,grid2_lon,mask, & 
     127       remap_matrix,source_add,destination_add) 
     128    ! 
     129    !----------------------------------------------------------------------- 
     130    !this routine makes any necessary changes (e.g. for 0,2pi longitude range) 
     131    !----------------------------------------------------------------------- 
     132    ! 
     133    REAL*8,DIMENSION(:,:) :: grid1_lat,grid2_lat,grid1_lon,grid2_lon 
     134    LOGICAL,DIMENSION(:,:) :: mask 
     135    !       
     136    INTEGER,DIMENSION(:),ALLOCATABLE :: source_add,destination_add   
     137    REAL*8,DIMENSION(:,:), ALLOCATABLE :: remap_matrix        
     138    ! 
     139    !----------------------------------------------------------------------- 
     140    ! local variables 
     141    !----------------------------------------------------------------------- 
     142    ! 
     143    INTEGER :: n,nele,i,j,ip1,jp1,n_add,e_add,ne_add,nx,ny 
     144    INTEGER :: xpos,ypos 
     145    !           
     146    ! integer mask 
     147    !  
     148    INTEGER, DIMENSION(:), ALLOCATABLE :: imask  
     149    ! 
     150    ! lat/lon intervals for search bins 
     151    ! 
     152    REAL*8 :: dlat,dlon            
     153    !       
     154    ! temps for computing bounding boxes 
     155    ! 
     156    REAL*8, DIMENSION(4) :: tmp_lats, tmp_lons   
     157    ! 
     158    !      write(*,*)'proceed to Bilinear interpolation ...' 
     159    ! 
     160!    IF(ALLOCATED(wts_map1)) DEALLOCATE(wts_map1) 
     161!    IF(ALLOCATED(grid1_add_map1)) DEALLOCATE(grid1_add_map1) 
     162!    IF(ALLOCATED(grid2_add_map1)) DEALLOCATE(grid2_add_map1) 
     163 
     164 
     165    ! 
     166    ALLOCATE(grid1_dims(2),grid2_dims(2)) 
     167    !       
     168    grid1_dims(1) = SIZE(grid1_lat,2) 
     169    grid1_dims(2) = SIZE(grid1_lat,1) 
     170    grid2_dims(1) = SIZE(grid2_lat,2) 
     171    grid2_dims(2) = SIZE(grid2_lat,1) 
     172    grid1_size = SIZE(grid1_lat,2) * SIZE(grid1_lat,1) 
     173    grid2_size = SIZE(grid2_lat,2) * SIZE(grid2_lat,1)   
     174    !       
     175    !----------------------------------------------------------------------- 
     176    !     allocate grid coordinates/masks and read data 
     177    !----------------------------------------------------------------------- 
     178    !      
     179    ALLOCATE( grid2_mask(grid2_size),         & 
     180         grid1_bound_box (4,grid1_size), & 
     181         grid2_bound_box (4,grid2_size), & 
     182         grid1_frac      (grid1_size),   & 
     183         grid2_frac      (grid2_size)) 
     184    ALLOCATE(imask(grid1_size)) 
     185    !                  
     186    !       
     187    grid1_frac = zero 
     188    grid2_frac = zero 
     189 
     190    ! 
     191    ! 2D array -> 1D array 
     192    ! 
     193    ALLOCATE(grid1_center_lat(SIZE(grid1_lat,1)*SIZE(grid1_lat,2))) 
     194    CALL tab2Dto1D(grid1_lat,grid1_center_lat) 
     195 
     196    ALLOCATE(grid1_center_lon(SIZE(grid1_lon,1)*SIZE(grid1_lon,2))) 
     197    CALL tab2Dto1D(grid1_lon,grid1_center_lon) 
     198 
     199    ALLOCATE(grid2_center_lat(SIZE(grid2_lat,1)*SIZE(grid2_lat,2))) 
     200    CALL tab2Dto1D(grid2_lat,grid2_center_lat) 
     201 
     202    ALLOCATE(grid2_center_lon(SIZE(grid2_lon,1)*SIZE(grid2_lon,2)))       
     203    CALL tab2Dto1D(grid2_lon,grid2_center_lon)  
     204 
     205    ALLOCATE(grid1_mask(SIZE(grid1_lat,1)*SIZE(grid1_lat,2))) 
     206    CALL logtab2Dto1D(mask,grid1_mask) 
     207    !       
     208    !      Write(*,*) ,'grid1_mask = ',grid1_mask                  
     209    ! 
     210    ! degrees to radian 
     211    ! 
     212    grid1_center_lat = grid1_center_lat*deg2rad 
     213    grid1_center_lon = grid1_center_lon*deg2rad 
     214    grid2_center_lat = grid2_center_lat*deg2rad 
     215    grid2_center_lon = grid2_center_lon*deg2rad 
     216 
     217    !----------------------------------------------------------------------- 
     218    !     convert longitudes to 0,2pi interval 
     219    !----------------------------------------------------------------------- 
     220 
     221    WHERE (grid1_center_lon .GT. pi2)  grid1_center_lon =       & 
     222         grid1_center_lon - pi2 
     223    WHERE (grid1_center_lon .LT. zero) grid1_center_lon =       & 
     224         grid1_center_lon + pi2 
     225    WHERE (grid2_center_lon .GT. pi2)  grid2_center_lon =       & 
     226         grid2_center_lon - pi2 
     227    WHERE (grid2_center_lon .LT. zero) grid2_center_lon =       & 
     228         grid2_center_lon + pi2 
     229 
     230    !----------------------------------------------------------------------- 
     231    ! 
     232    !     make sure input latitude range is within the machine values 
     233    !     for +/- pi/2  
     234    ! 
     235    !----------------------------------------------------------------------- 
     236 
     237    WHERE (grid1_center_lat >  pih) grid1_center_lat =  pih 
     238    WHERE (grid1_center_lat < -pih) grid1_center_lat = -pih 
     239    WHERE (grid2_center_lat >  pih) grid2_center_lat =  pih 
     240    WHERE (grid2_center_lat < -pih) grid2_center_lat = -pih 
     241 
     242    !----------------------------------------------------------------------- 
     243    ! 
     244    !     compute bounding boxes for restricting future grid searches 
     245    ! 
     246    !----------------------------------------------------------------------- 
     247    ! 
     248    nx = grid1_dims(1) 
     249    ny = grid1_dims(2) 
     250 
     251    DO n=1,grid1_size 
     252 
     253       !*** find N,S and NE points to this grid point 
     254 
     255       j = (n - 1)/nx +1 
     256       i = n - (j-1)*nx 
     257 
     258       IF (i < nx) THEN 
     259          ip1 = i + 1 
     260       ELSE 
     261          !*** assume cyclic 
     262          ip1 = 1 
     263          !*** but if it is not, correct 
     264          e_add = (j - 1)*nx + ip1 
     265          IF (ABS(grid1_center_lat(e_add) -     & 
     266               grid1_center_lat(n   )) > pih) THEN 
     267             ip1 = i 
     268          ENDIF 
     269          ip1=nx 
     270       ENDIF 
     271 
     272       IF (j < ny) THEN 
     273          jp1 = j+1 
     274       ELSE 
     275          !*** assume cyclic 
     276          jp1 = 1 
     277          !*** but if it is not, correct 
     278          n_add = (jp1 - 1)*nx + i 
     279          IF (ABS(grid1_center_lat(n_add) -             & 
     280               grid1_center_lat(n   )) > pih) THEN 
     281             jp1 = j 
     282          ENDIF 
     283          jp1=ny 
     284       ENDIF 
     285 
     286       n_add = (jp1 - 1)*nx + i 
     287       e_add = (j - 1)*nx + ip1 
     288       ne_add = (jp1 - 1)*nx + ip1 
     289 
     290       !*** find N,S and NE lat/lon coords and check bounding box 
     291 
     292       tmp_lats(1) = grid1_center_lat(n) 
     293       tmp_lats(2) = grid1_center_lat(e_add) 
     294       tmp_lats(3) = grid1_center_lat(ne_add) 
     295       tmp_lats(4) = grid1_center_lat(n_add) 
     296 
     297       tmp_lons(1) = grid1_center_lon(n) 
     298       tmp_lons(2) = grid1_center_lon(e_add) 
     299       tmp_lons(3) = grid1_center_lon(ne_add) 
     300       tmp_lons(4) = grid1_center_lon(n_add) 
     301 
     302       grid1_bound_box(1,n) = MINVAL(tmp_lats) 
     303       grid1_bound_box(2,n) = MAXVAL(tmp_lats) 
     304 
     305       grid1_bound_box(3,n) = MINVAL(tmp_lons) 
     306       grid1_bound_box(4,n) = MAXVAL(tmp_lons) 
     307    END DO 
     308 
     309    nx = grid2_dims(1) 
     310    ny = grid2_dims(2) 
     311 
     312    DO n=1,grid2_size 
     313 
     314       !*** find N,S and NE points to this grid point 
     315 
     316       j = (n - 1)/nx +1 
     317       i = n - (j-1)*nx 
     318 
     319       IF (i < nx) THEN 
     320          ip1 = i + 1 
     321       ELSE 
     322          !*** assume cyclic 
     323          ip1 = 1 
     324          !*** but if it is not, correct 
     325          e_add = (j - 1)*nx + ip1 
     326          IF (ABS(grid2_center_lat(e_add) -  & 
     327               grid2_center_lat(n   )) > pih) THEN 
     328             ip1 = i 
     329          ENDIF 
     330       ENDIF 
     331 
     332       IF (j < ny) THEN 
     333          jp1 = j+1 
     334       ELSE 
     335          !*** assume cyclic 
     336          jp1 = 1 
     337          !*** but if it is not, correct 
     338          n_add = (jp1 - 1)*nx + i 
     339          IF (ABS(grid2_center_lat(n_add) -  & 
     340               grid2_center_lat(n   )) > pih) THEN 
     341             jp1 = j 
     342          ENDIF 
     343       ENDIF 
     344 
     345       n_add = (jp1 - 1)*nx + i 
     346       e_add = (j - 1)*nx + ip1 
     347       ne_add = (jp1 - 1)*nx + ip1 
     348 
     349       !*** find N,S and NE lat/lon coords and check bounding box 
     350 
     351       tmp_lats(1) = grid2_center_lat(n) 
     352       tmp_lats(2) = grid2_center_lat(e_add) 
     353       tmp_lats(3) = grid2_center_lat(ne_add) 
     354       tmp_lats(4) = grid2_center_lat(n_add) 
     355 
     356       tmp_lons(1) = grid2_center_lon(n) 
     357       tmp_lons(2) = grid2_center_lon(e_add) 
     358       tmp_lons(3) = grid2_center_lon(ne_add) 
     359       tmp_lons(4) = grid2_center_lon(n_add) 
     360 
     361       grid2_bound_box(1,n) = MINVAL(tmp_lats) 
     362       grid2_bound_box(2,n) = MAXVAL(tmp_lats) 
     363       grid2_bound_box(3,n) = MINVAL(tmp_lons) 
     364       grid2_bound_box(4,n) = MAXVAL(tmp_lons) 
     365    END DO 
     366    ! 
     367    ! 
     368    ! 
     369    WHERE (ABS(grid1_bound_box(4,:) - grid1_bound_box(3,:)) > pi) 
     370       grid1_bound_box(3,:) = zero 
     371       grid1_bound_box(4,:) = pi2 
     372    END WHERE 
     373 
     374    WHERE (ABS(grid2_bound_box(4,:) - grid2_bound_box(3,:)) > pi) 
     375       grid2_bound_box(3,:) = zero 
     376       grid2_bound_box(4,:) = pi2 
     377    END WHERE 
     378 
     379    !*** 
     380    !*** try to check for cells that overlap poles 
     381    !*** 
     382 
     383    WHERE (grid1_center_lat > grid1_bound_box(2,:)) & 
     384         grid1_bound_box(2,:) = pih 
     385 
     386    WHERE (grid1_center_lat < grid1_bound_box(1,:)) & 
     387         grid1_bound_box(1,:) = -pih 
     388 
     389    WHERE (grid2_center_lat > grid2_bound_box(2,:)) & 
     390         grid2_bound_box(2,:) = pih 
     391 
     392    WHERE (grid2_center_lat < grid2_bound_box(1,:)) & 
     393         grid2_bound_box(1,:) = -pih 
     394 
     395    !----------------------------------------------------------------------- 
     396    !     set up and assign address ranges to search bins in order to  
     397    !     further restrict later searches 
     398    !----------------------------------------------------------------------- 
     399 
     400    ALLOCATE(bin_addr1(2,num_srch_bins)) 
     401    ALLOCATE(bin_addr2(2,num_srch_bins)) 
     402    ALLOCATE(bin_lats (2,num_srch_bins)) 
     403    ALLOCATE(bin_lons (2,num_srch_bins)) 
     404 
     405    dlat = pi/num_srch_bins 
     406 
     407    DO n=1,num_srch_bins 
     408       bin_lats(1,n) = (n-1)*dlat - pih 
     409       bin_lats(2,n) =     n*dlat - pih 
     410       bin_lons(1,n) = zero 
     411       bin_lons(2,n) = pi2 
     412       bin_addr1(1,n) = grid1_size + 1 
     413       bin_addr1(2,n) = 0 
     414       bin_addr2(1,n) = grid2_size + 1 
     415       bin_addr2(2,n) = 0 
     416    END DO 
     417 
     418    DO nele=1,grid1_size 
     419       DO n=1,num_srch_bins 
     420          IF (grid1_bound_box(1,nele) <= bin_lats(2,n) .AND.   & 
     421               grid1_bound_box(2,nele) >= bin_lats(1,n)) THEN 
     422             bin_addr1(1,n) = MIN(nele,bin_addr1(1,n)) 
     423             bin_addr1(2,n) = MAX(nele,bin_addr1(2,n)) 
     424          ENDIF 
     425       END DO 
     426    END DO 
     427 
     428    DO nele=1,grid2_size 
     429       DO n=1,num_srch_bins 
     430          IF (grid2_bound_box(1,nele) <= bin_lats(2,n) .AND.    & 
     431               grid2_bound_box(2,nele) >= bin_lats(1,n)) THEN 
     432             bin_addr2(1,n) = MIN(nele,bin_addr2(1,n)) 
     433             bin_addr2(2,n) = MAX(nele,bin_addr2(2,n)) 
     434          ENDIF 
     435       END DO 
     436    END DO 
     437    !  
     438    CALL init_remap_vars  
     439    CALL remap_bilin 
     440 
     441    ALLOCATE(remap_matrix(SIZE(wts_map1,1),SIZE(wts_map1,2)), & 
     442         source_add(SIZE(grid1_add_map1)),       & 
     443         destination_add(SIZE(grid2_add_map1))) 
     444 
     445    DO j = 1,SIZE(wts_map1,2) 
     446       DO i = 1,SIZE(wts_map1,1) 
     447 
     448          remap_matrix(i,j) = wts_map1(i,j) 
     449 
     450       END DO 
     451    END DO 
     452 
     453 
     454    source_add(:) = grid1_add_map1(:) 
     455    destination_add(:) = grid2_add_map1(:)                
     456    ! 
     457    WHERE(destination_add == 0) 
     458       destination_add = 1 
     459    END WHERE 
     460 
     461    WHERE(source_add == 0) 
     462       source_add = 1 
     463    END WHERE 
     464    !                
     465    !DEALLOCATE(grid1_bound_box,grid2_bound_box,grid1_center_lat,grid1_center_lon) 
     466    !DEALLOCATE(grid2_center_lat,grid2_center_lon,grid2_add_map1,grid1_add_map1,wts_map1) 
     467    !DEALLOCATE(grid1_frac,grid2_frac,grid1_dims,grid2_dims,grid2_mask,imask) 
     468    !DEALLOCATE(bin_addr1,bin_addr2,bin_lats,bin_lons) 
     469    !DEALLOCATE(grid1_mask)    
     470    !  
     471    !----------------------------------------------------------------------- 
     472 
     473  END SUBROUTINE get_remap_matrix 
     474 
     475 
     476 
     477 
     478 
     479 
     480 
     481 
     482 
     483 
     484 
     485 
     486 
     487 
     488 
     489 
     490 
     491 
     492 
     493 
     494 
     495 
     496 
     497 
     498  !*********************************************************************** 
     499!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     500  !************************************************************************ 
     501  !   SUBROUTINE REMAP_BILINEAR 
     502  !************************************************************************ 
     503!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     504 
     505  SUBROUTINE remap_bilin 
     506 
     507    !----------------------------------------------------------------------- 
     508    !     this routine computes the weights for a bilinear interpolation. 
     509    !----------------------------------------------------------------------- 
     510 
     511    !----------------------------------------------------------------------- 
     512    !     local variables 
     513    !----------------------------------------------------------------------- 
     514 
     515    INTEGER :: n,icount,dst_add,iter,nmap,nbmasked     
     516    !         
     517    ! address for the four source points 
     518    ! 
     519    INTEGER, DIMENSION(4) :: src_add,involved_pts 
     520    INTEGER, DIMENSION(1) :: minlon 
     521    INTEGER, DIMENSION(1) :: minlat 
     522    REAL*8, DIMENSION(4) :: distx,disty 
     523    REAL*8 :: normalize 
     524    !                
     525    ! latitudes longitudes of four bilinear corners 
     526    ! 
     527    REAL*8, DIMENSION(4) :: src_lats,src_lons 
     528    ! 
     529    ! bilinear weights for four corners 
     530    !       
     531    REAL*8, DIMENSION(4) :: wgts             
     532    ! 
     533    REAL*8 :: & 
     534         plat, plon,       &  ! lat/lon coords of destination point 
     535         iguess, jguess,   &  ! current guess for bilinear coordinate 
     536         thguess, phguess, &  ! current guess for lat/lon coordinate 
     537         deli, delj,       &  ! corrections to i,j 
     538         dth1, dth2, dth3, &  ! some latitude  differences 
     539         dph1, dph2, dph3, &  ! some longitude differences 
     540         dthp, dphp,       &  ! difference between point and sw corner 
     541         mat1, mat2, mat3, mat4, &  ! matrix elements 
     542         determinant, &     ! matrix determinant 
     543         sum_wgts          ! sum of weights for normalization 
     544 
     545    INTEGER lastsrc_add   
     546    REAL*8,DIMENSION(:),ALLOCATABLE :: cos_grid1_center_lat    
     547    REAL*8,DIMENSION(:),ALLOCATABLE :: cos_grid1_center_lon 
     548    REAL*8,DIMENSION(:),ALLOCATABLE :: sin_grid1_center_lat 
     549    REAL*8,DIMENSION(:),ALLOCATABLE :: sin_grid1_center_lon 
     550    INTEGER :: i 
     551    !       
     552    grid2_mask = .TRUE. 
     553    !       
     554    !       
     555    nmap = 1 
     556    ! 
     557    !*** 
     558    !*** loop over destination grid  
     559    !*** 
     560    !      print*,'grid2_size =',grid2_size 
     561    !       
     562    lastsrc_add=1 
     563    ! 
     564    Allocate(cos_grid1_center_lat(size(grid1_center_lat))) 
     565    Allocate(sin_grid1_center_lat(size(grid1_center_lat))) 
     566    Allocate(cos_grid1_center_lon(size(grid1_center_lon))) 
     567    Allocate(sin_grid1_center_lon(size(grid1_center_lon))) 
     568 
     569    do i=1,size(grid1_center_lat) 
     570        cos_grid1_center_lat(i)=cos(grid1_center_lat(i)) 
     571        sin_grid1_center_lat(i)=sin(grid1_center_lat(i)) 
     572    ENDDO 
     573 
     574    do i=1,size(grid1_center_lon) 
     575        cos_grid1_center_lon(i)=cos(grid1_center_lon(i)) 
     576        sin_grid1_center_lon(i)=sin(grid1_center_lon(i)) 
     577    ENDDO 
     578 
     579    grid_loop1: DO dst_add = 1, grid2_size 
     580 
     581       IF (.NOT. grid2_mask(dst_add)) CYCLE grid_loop1 
     582       !         
     583       plat = grid2_center_lat(dst_add) 
     584       plon = grid2_center_lon(dst_add) 
     585       !*** 
     586       !*** find nearest square of grid points on source grid 
     587       !*** 
     588       CALL grid_search_bilin(src_add, src_lats, src_lons,          & 
     589            plat, plon, grid1_dims,               & 
     590            grid1_center_lat, cos_grid1_center_lat, sin_grid1_center_lat, & 
     591            grid1_center_lon, cos_grid1_center_lon, sin_grid1_center_lon,  &  
     592            grid1_bound_box, bin_addr1, bin_addr2,lastsrc_add)            
     593       !*** 
     594       !*** check to see if points are land points 
     595       !***  
     596       !    
     597       IF (src_add(1) > 0) THEN 
     598          !      
     599          DO n=1,4 
     600             !           if(.not. grid1_mask(src_add(n))) nbmasked = nbmasked + 1 
     601             IF(.NOT. grid1_mask(src_add(n))) src_add(1) = 0 
     602          END DO 
     603          ! 
     604       ENDIF 
     605       ! 
     606       ! 
     607 
     608       !*** 
     609       !*** if point found, find local i,j coordinates for weights 
     610       !*** 
     611       IF (src_add(1) > 0) THEN 
     612          grid2_frac(dst_add) = one 
     613          !*** 
     614          !*** iterate to find i,j for bilinear approximation 
     615          !*** 
     616          dth1 = src_lats(2) - src_lats(1) 
     617          dth2 = src_lats(4) - src_lats(1) 
     618          dth3 = src_lats(3) - src_lats(2) - dth2 
     619 
     620          dph1 = src_lons(2) - src_lons(1) 
     621          dph2 = src_lons(4) - src_lons(1) 
     622          dph3 = src_lons(3) - src_lons(2) 
     623 
     624          IF (dph1 >  three*pih) dph1 = dph1 - pi2 
     625          IF (dph2 >  three*pih) dph2 = dph2 - pi2 
     626          IF (dph3 >  three*pih) dph3 = dph3 - pi2 
     627          IF (dph1 < -three*pih) dph1 = dph1 + pi2 
     628          IF (dph2 < -three*pih) dph2 = dph2 + pi2 
     629          IF (dph3 < -three*pih) dph3 = dph3 + pi2 
     630 
     631          dph3 = dph3 - dph2 
     632 
     633          iguess = half 
     634          jguess = half 
     635 
     636          iter_loop1: DO iter=1,max_iter 
     637 
     638             dthp = plat - src_lats(1) - dth1*iguess -        & 
     639                  dth2*jguess - dth3*iguess*jguess 
     640             dphp = plon - src_lons(1) 
     641 
     642             IF (dphp >  three*pih) dphp = dphp - pi2 
     643             IF (dphp < -three*pih) dphp = dphp + pi2 
     644 
     645             dphp = dphp - dph1*iguess - dph2*jguess -        & 
     646                  dph3*iguess*jguess 
     647 
     648             mat1 = dth1 + dth3*jguess 
     649             mat2 = dth2 + dth3*iguess 
     650             mat3 = dph1 + dph3*jguess 
     651             mat4 = dph2 + dph3*iguess 
     652 
     653             determinant = mat1*mat4 - mat2*mat3 
     654 
     655             deli = (dthp*mat4 - mat2*dphp)/determinant 
     656             delj = (mat1*dphp - dthp*mat3)/determinant 
     657 
     658             IF (ABS(deli) < converge .AND.                   & 
     659                  ABS(delj) < converge) EXIT iter_loop1 
     660 
     661             iguess = iguess + deli 
     662             jguess = jguess + delj 
     663 
     664          END DO iter_loop1 
     665 
     666          IF (iter <= max_iter) THEN 
     667 
     668             !*** 
     669             !*** successfully found i,j - compute weights 
     670             !*** 
     671 
     672             wgts(1) = (one-iguess)*(one-jguess) 
     673             wgts(2) = iguess*(one-jguess) 
     674             wgts(3) = iguess*jguess 
     675             wgts(4) = (one-iguess)*jguess 
     676             !                
     677             ! 
     678             CALL store_link_bilin(dst_add, src_add, wgts, nmap) 
     679 
     680          ELSE 
     681             PRINT *,'Point coords: ',plat,plon 
     682             PRINT *,'Dest grid lats: ',src_lats 
     683             PRINT *,'Dest grid lons: ',src_lons 
     684             PRINT *,'Dest grid addresses: ',src_add 
     685             PRINT *,'Current i,j : ',iguess, jguess 
     686             STOP 'Iteration for i,j exceed max iteration count' 
     687          ENDIF 
     688          ! 
     689          !*** 
     690          !*** search for bilinear failed - use a distance-weighted 
     691          !*** average instead (this is typically near the pole) 
     692          !*** 
     693       ELSE IF (src_add(1) < 0) THEN 
     694 
     695          src_add = ABS(src_add) 
     696          icount = 0 
     697          DO n=1,4 
     698             !            
     699             IF (grid1_mask(src_add(n))) THEN 
     700                icount = icount + 1 
     701             ELSE 
     702                src_lats(n) = zero 
     703             ENDIF 
     704             !      
     705          END DO 
     706 
     707          IF (icount > 0) THEN 
     708             !    
     709             !*** renormalize weights 
     710             ! 
     711             sum_wgts = SUM(src_lats) 
     712             wgts(1) = src_lats(1)/sum_wgts 
     713             wgts(2) = src_lats(2)/sum_wgts 
     714             wgts(3) = src_lats(3)/sum_wgts 
     715             wgts(4) = src_lats(4)/sum_wgts 
     716             ! 
     717             grid2_frac(dst_add) = one 
     718             CALL store_link_bilin(dst_add, src_add, wgts, nmap) 
     719          ENDIF 
     720 
     721          ! 
     722       ENDIF 
     723    END DO grid_loop1 
     724    !       
     725    !      Call sort_add(grid2_add_map1, grid1_add_map1, wts_map1)    
     726    !                
     727    ! 
     728    !----------------------------------------------------------------------- 
     729 
     730    deallocate(cos_grid1_center_lat) 
     731    deallocate(cos_grid1_center_lon) 
     732    deallocate(sin_grid1_center_lat) 
     733    deallocate(sin_grid1_center_lon) 
     734  END SUBROUTINE remap_bilin 
     735 
     736 
     737 
     738 
     739 
     740 
     741 
     742 
     743 
     744 
     745 
     746 
     747 
     748 
     749 
     750 
     751 
     752 
     753  !*********************************************************************** 
     754!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     755  !************************************************************************ 
     756  !   SUBROUTINE GRID_SEARCH_BILIN 
     757  !************************************************************************ 
     758!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     759  ! 
     760  SUBROUTINE grid_search_bilin(src_add, src_lats, src_lons,   & 
     761       plat, plon, src_grid_dims,      & 
     762       src_center_lat, cos_src_center_lat, sin_src_center_lat, & 
     763       src_center_lon, cos_src_center_lon, sin_src_center_lon,&  
     764       src_grid_bound_box,             & 
     765       src_bin_add, dst_bin_add,lastsrc_add) 
     766 
     767    !----------------------------------------------------------------------- 
     768    ! 
     769    !     this routine finds the location of the search point plat, plon 
     770    !     in the source grid and returns the corners needed for a bilinear 
     771    !     interpolation. 
     772    ! 
     773    !----------------------------------------------------------------------- 
     774 
     775    !----------------------------------------------------------------------- 
     776    !     output variables 
     777    !----------------------------------------------------------------------- 
     778    ! 
     779    ! address of each corner point enclosing P 
     780    ! 
     781    INTEGER,DIMENSION(4) :: src_add   
     782    REAL*8,DIMENSION(4) :: src_lats,src_lons   
     783    !       
     784    !----------------------------------------------------------------------- 
     785    !     input variables 
     786    !----------------------------------------------------------------------- 
     787    !  
     788    ! latitude, longitude of the search point 
     789    ! 
     790    REAL*8, INTENT(in) :: plat,plon    
     791    ! 
     792    ! size of each src grid dimension 
     793    ! 
     794    INTEGER, DIMENSION(2), INTENT(in) :: src_grid_dims   
     795    ! 
     796    ! latitude, longitude of each src grid center 
     797    ! 
     798    REAL*8, DIMENSION(:), INTENT(in) :: src_center_lat,src_center_lon   
     799    REAL*8, DIMENSION(:), INTENT(in) :: cos_src_center_lat 
     800    REAL*8, DIMENSION(:), INTENT(in) :: sin_src_center_lat 
     801    REAL*8, DIMENSION(:), INTENT(in) :: cos_src_center_lon 
     802    REAL*8, DIMENSION(:), INTENT(in) :: sin_src_center_lon 
     803    ! 
     804    ! bound box for source grid 
     805    ! 
     806    REAL*8, DIMENSION(:,:), INTENT(in) :: src_grid_bound_box  
     807    ! 
     808    ! latitude bins for restricting searches 
     809    ! 
     810    INTEGER, DIMENSION(:,:), INTENT(in) ::src_bin_add,dst_bin_add  
     811 
     812    INTEGER,OPTIONAL :: lastsrc_add 
     813    INTEGER :: loopsrc,l1,l2 
     814    !       
     815    !----------------------------------------------------------------------- 
     816    !     local variables 
     817    !----------------------------------------------------------------------- 
     818    ! 
     819    INTEGER :: n,next_n,srch_add,nx, ny,min_add, max_add,  & 
     820         i, j, jp1, ip1, n_add, e_add, ne_add 
     821 
     822 
     823    REAL*8 ::  vec1_lat, vec1_lon,vec2_lat, vec2_lon, cross_product,  & 
     824         cross_product_last,coslat_dst, sinlat_dst, coslon_dst, & 
     825         sinlon_dst,dist_min, distance  
     826 
     827    !----------------------------------------------------------------------- 
     828    !     restrict search first using bins 
     829    !----------------------------------------------------------------------- 
     830 
     831    src_add = 0 
     832 
     833    min_add = SIZE(src_center_lat) 
     834    max_add = 1 
     835    DO n=1,num_srch_bins 
     836       IF (plat >= bin_lats(1,n) .AND. plat <= bin_lats(2,n) .AND. & 
     837            plon >= bin_lons(1,n) .AND. plon <= bin_lons(2,n)) THEN 
     838          min_add = MIN(min_add, src_bin_add(1,n)) 
     839          max_add = MAX(max_add, src_bin_add(2,n)) 
     840       ENDIF 
     841    END DO 
     842 
     843    !----------------------------------------------------------------------- 
     844    !     now perform a more detailed search  
     845    !----------------------------------------------------------------------- 
     846 
     847    nx = src_grid_dims(1) 
     848    ny = src_grid_dims(2) 
     849 
     850    loopsrc=0 
     851    DO WHILE (loopsrc <= max_add) 
     852 
     853 
     854       l1=MAX(min_add,lastsrc_add-loopsrc) 
     855       l2=MIN(max_add,lastsrc_add+loopsrc)       
     856 
     857       loopsrc = loopsrc+1 
     858 
     859       srch_loop: DO srch_add = l1,l2,MAX(l2-l1,1) 
     860 
     861          !*** first check bounding box 
     862 
     863          IF (plat <= src_grid_bound_box(2,srch_add) .AND. &  
     864               plat >= src_grid_bound_box(1,srch_add) .AND.  & 
     865               plon <= src_grid_bound_box(4,srch_add) .AND.  & 
     866               plon >= src_grid_bound_box(3,srch_add)) THEN 
     867             !*** 
     868             !*** we are within bounding box so get really serious 
     869             !*** 
     870             !*** determine neighbor addresses 
     871             ! 
     872             j = (srch_add - 1)/nx +1 
     873             i = srch_add - (j-1)*nx 
     874             ! 
     875             IF (i < nx) THEN 
     876                ip1 = i + 1 
     877             ELSE 
     878                ip1 = 1 
     879             ENDIF 
     880             ! 
     881             IF (j < ny) THEN 
     882                jp1 = j+1 
     883             ELSE 
     884                jp1 = 1 
     885             ENDIF 
     886             ! 
     887             n_add = (jp1 - 1)*nx + i 
     888             e_add = (j - 1)*nx + ip1 
     889             ne_add = (jp1 - 1)*nx + ip1 
     890             ! 
     891             src_lats(1) = src_center_lat(srch_add) 
     892             src_lats(2) = src_center_lat(e_add) 
     893             src_lats(3) = src_center_lat(ne_add) 
     894             src_lats(4) = src_center_lat(n_add) 
     895             ! 
     896             src_lons(1) = src_center_lon(srch_add) 
     897             src_lons(2) = src_center_lon(e_add) 
     898             src_lons(3) = src_center_lon(ne_add) 
     899             src_lons(4) = src_center_lon(n_add) 
     900             ! 
     901             !*** 
     902             !*** for consistency, we must make sure all lons are in 
     903             !*** same 2pi interval 
     904             !*** 
     905             ! 
     906             vec1_lon = src_lons(1) - plon 
     907             IF (vec1_lon >  pi) THEN 
     908                src_lons(1) = src_lons(1) - pi2 
     909             ELSE IF (vec1_lon < -pi) THEN 
     910                src_lons(1) = src_lons(1) + pi2 
     911             ENDIF 
     912             DO n=2,4 
     913                vec1_lon = src_lons(n) - src_lons(1) 
     914                IF (vec1_lon >  pi) THEN 
     915                   src_lons(n) = src_lons(n) - pi2 
     916                ELSE IF (vec1_lon < -pi) THEN 
     917                   src_lons(n) = src_lons(n) + pi2 
     918                ENDIF 
     919             END DO 
     920             ! 
     921             corner_loop: DO n=1,4 
     922                next_n = MOD(n,4) + 1 
     923                !*** 
     924                !*** here we take the cross product of the vector making  
     925                !*** up each box side with the vector formed by the vertex 
     926                !*** and search point.  if all the cross products are  
     927                !*** positive, the point is contained in the box. 
     928                !*** 
     929                vec1_lat = src_lats(next_n) - src_lats(n) 
     930                vec1_lon = src_lons(next_n) - src_lons(n) 
     931                vec2_lat = plat - src_lats(n) 
     932                vec2_lon = plon - src_lons(n) 
     933                !*** 
     934                !*** check for 0,2pi crossings 
     935                !*** 
     936                IF (vec1_lon >  three*pih) THEN 
     937                   vec1_lon = vec1_lon - pi2 
     938                ELSE IF (vec1_lon < -three*pih) THEN 
     939                   vec1_lon = vec1_lon + pi2 
     940                ENDIF 
     941                IF (vec2_lon >  three*pih) THEN 
     942                   vec2_lon = vec2_lon - pi2 
     943                ELSE IF (vec2_lon < -three*pih) THEN 
     944                   vec2_lon = vec2_lon + pi2 
     945                ENDIF 
     946                ! 
     947                cross_product = vec1_lon*vec2_lat - vec2_lon*vec1_lat 
     948                ! 
     949                !*** 
     950                !*** if cross product is less than zero, this cell 
     951                !*** doesn't work 
     952                !*** 
     953                IF (n == 1) cross_product_last = cross_product 
     954                IF (cross_product*cross_product_last < zero) & 
     955                     EXIT corner_loop 
     956                cross_product_last = cross_product 
     957                ! 
     958             END DO corner_loop 
     959             !*** 
     960             !*** if cross products all same sign, we found the location 
     961             !*** 
     962             IF (n > 4) THEN 
     963                src_add(1) = srch_add 
     964                src_add(2) = e_add 
     965                src_add(3) = ne_add 
     966                src_add(4) = n_add 
     967                ! 
     968                lastsrc_add = srch_add 
     969                RETURN 
     970             ENDIF 
     971             !*** 
     972             !*** otherwise move on to next cell 
     973             !*** 
     974          ENDIF !bounding box check 
     975       END DO srch_loop 
     976 
     977 
     978    ENDDO 
     979 
     980 
     981    !*** 
     982    !*** if no cell found, point is likely either in a box that 
     983    !*** straddles either pole or is outside the grid.  fall back 
     984    !*** to a distance-weighted average of the four closest 
     985    !*** points.  go ahead and compute weights here, but store 
     986    !*** in src_lats and return -add to prevent the parent 
     987    !*** routine from computing bilinear weights 
     988    !*** 
     989    !print *,'Could not find location for ',plat,plon 
     990    !print *,'Using nearest-neighbor average for this point' 
     991    ! 
     992    coslat_dst = COS(plat) 
     993    sinlat_dst = SIN(plat) 
     994    coslon_dst = COS(plon) 
     995    sinlon_dst = SIN(plon) 
     996    ! 
     997    dist_min = bignum 
     998    src_lats = bignum 
     999    DO srch_add = min_add,max_add 
     1000       ! distance = ACOS(coslat_dst*COS(src_center_lat(srch_add))*   & 
     1001       !      (coslon_dst*COS(src_center_lon(srch_add)) +   & 
     1002       !      sinlon_dst*SIN(src_center_lon(srch_add)))+   & 
     1003       !      sinlat_dst*SIN(src_center_lat(srch_add))) 
     1004 
     1005       distance = ACOS(coslat_dst*cos_src_center_lat(srch_add)*   & 
     1006            (coslon_dst*cos_src_center_lon(srch_add) +   & 
     1007            sinlon_dst*sin_src_center_lon(srch_add))+   & 
     1008            sinlat_dst*sin_src_center_lat(srch_add)) 
     1009 
     1010       IF (distance < dist_min) THEN 
     1011          sort_loop: DO n=1,4 
     1012             IF (distance < src_lats(n)) THEN 
     1013                DO i=4,n+1,-1 
     1014                   src_add (i) = src_add (i-1) 
     1015                   src_lats(i) = src_lats(i-1) 
     1016                END DO 
     1017                src_add (n) = -srch_add 
     1018                src_lats(n) = distance 
     1019                dist_min = src_lats(4) 
     1020                EXIT sort_loop 
     1021             ENDIF 
     1022          END DO sort_loop 
     1023       ENDIF 
     1024    END DO 
     1025    ! 
     1026    src_lons = one/(src_lats + tiny) 
     1027    distance = SUM(src_lons) 
     1028    src_lats = src_lons/distance 
     1029 
     1030    !----------------------------------------------------------------------- 
     1031 
     1032  END SUBROUTINE grid_search_bilin 
     1033 
     1034 
     1035  !*********************************************************************** 
     1036!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     1037  !************************************************************************ 
     1038  !   SUBROUTINE STORE_LINK_BILIN 
     1039  !************************************************************************ 
     1040!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     1041 
     1042  SUBROUTINE store_link_bilin(dst_add, src_add, weights, nmap) 
     1043 
     1044    !----------------------------------------------------------------------- 
     1045    !     this routine stores the address and weight for four links  
     1046    !     associated with one destination point in the appropriate address  
     1047    !     and weight arrays and resizes those arrays if necessary. 
     1048    !----------------------------------------------------------------------- 
     1049 
     1050    !----------------------------------------------------------------------- 
     1051    !     input variables 
     1052    !----------------------------------------------------------------------- 
     1053    ! 
     1054    INTEGER :: dst_add,nmap 
     1055    ! 
     1056    INTEGER, DIMENSION(4) :: src_add 
     1057    ! 
     1058    REAL*8, DIMENSION(4) :: weights  
     1059 
     1060    !----------------------------------------------------------------------- 
     1061    ! 
     1062    !     local variables 
     1063    ! 
     1064    !----------------------------------------------------------------------- 
     1065 
     1066    INTEGER :: n,num_links_old    
     1067 
     1068    !----------------------------------------------------------------------- 
     1069    !     increment number of links and check to see if remap arrays need 
     1070    !     to be increased to accomodate the new link.  then store the 
     1071    !     link. 
     1072    !----------------------------------------------------------------------- 
     1073 
     1074    num_links_old  = num_links_map1 
     1075    num_links_map1 = num_links_old + 4 
     1076 
     1077    IF (num_links_map1 > max_links_map1) & 
     1078         CALL resize_remap_vars(1,resize_increment) 
     1079 
     1080    DO n=1,4 
     1081       grid1_add_map1(num_links_old+n) = src_add(n) 
     1082       grid2_add_map1(num_links_old+n) = dst_add 
     1083       wts_map1    (1,num_links_old+n) = weights(n) 
     1084    END DO 
     1085 
     1086    !----------------------------------------------------------------------- 
     1087 
     1088  END SUBROUTINE store_link_bilin 
     1089 
     1090 
     1091 
     1092 
     1093 
     1094 
     1095 
     1096 
     1097!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     1098  !************************************************************************ 
     1099  !   SUBROUTINE INIT_REMAP_VARS 
     1100  !************************************************************************ 
     1101!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     1102  ! 
     1103  SUBROUTINE init_remap_vars 
     1104 
     1105    !----------------------------------------------------------------------- 
     1106    ! 
     1107    !     this routine initializes some variables and provides an initial 
     1108    !     allocation of arrays (fairly large so frequent resizing  
     1109    !     unnecessary). 
     1110    ! 
     1111    !----------------------------------------------------------------------- 
     1112 
     1113    !----------------------------------------------------------------------- 
     1114    !     determine the number of weights 
     1115    !----------------------------------------------------------------------- 
     1116    num_wts = 1     ! bilinear interpolation 
     1117    !----------------------------------------------------------------------- 
     1118    !     initialize num_links and set max_links to four times the largest  
     1119    !     of the destination grid sizes initially (can be changed later). 
     1120    !     set a default resize increment to increase the size of link 
     1121    !     arrays if the number of links exceeds the initial size   
     1122    !----------------------------------------------------------------------- 
     1123 
     1124    num_links_map1 = 0 
     1125    max_links_map1 = 4*grid2_size 
     1126    IF (num_maps > 1) THEN 
     1127       num_links_map2 = 0 
     1128       max_links_map1 = MAX(4*grid1_size,4*grid2_size) 
     1129       max_links_map2 = max_links_map1 
     1130    ENDIF 
     1131 
     1132    resize_increment = 0.1*MAX(grid1_size,grid2_size) 
     1133 
     1134    !----------------------------------------------------------------------- 
     1135    !     allocate address and weight arrays for mapping 1   
     1136    !----------------------------------------------------------------------- 
     1137 
     1138    ALLOCATE (grid1_add_map1(max_links_map1),    & 
     1139         grid2_add_map1(max_links_map1),    & 
     1140         wts_map1(num_wts, max_links_map1)) 
     1141 
     1142    grid1_add_map1 = 0. 
     1143    grid2_add_map1 = 0. 
     1144    wts_map1 = 0. 
     1145 
     1146    !----------------------------------------------------------------------- 
     1147 
     1148  END SUBROUTINE init_remap_vars 
     1149 
     1150 
     1151 
     1152 
     1153 
     1154 
     1155 
     1156 
     1157 
     1158 
     1159 
     1160 
     1161 
     1162 
     1163  !*********************************************************************** 
     1164!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     1165  !************************************************************************ 
     1166  !   SUBROUTINE RESIZE_REMAP_VAR 
     1167  !************************************************************************ 
     1168!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     1169 
     1170  SUBROUTINE resize_remap_vars(nmap, increment) 
     1171 
     1172    !----------------------------------------------------------------------- 
     1173    !     this routine resizes remapping arrays by increasing(decreasing) 
     1174    !     the max_links by increment 
     1175    !----------------------------------------------------------------------- 
     1176 
     1177    !----------------------------------------------------------------------- 
     1178    !     input variables 
     1179    !----------------------------------------------------------------------- 
     1180 
     1181    INTEGER ::     & 
     1182         nmap,      &     ! identifies which mapping array to resize 
     1183         increment       ! the number of links to add(subtract) to arrays 
     1184 
     1185    !----------------------------------------------------------------------- 
     1186    !     local variables 
     1187    !----------------------------------------------------------------------- 
     1188 
     1189    INTEGER ::    & 
     1190         ierr,    &  ! error flag 
     1191         mxlinks   ! size of link arrays 
     1192 
     1193    INTEGER, DIMENSION(:), ALLOCATABLE ::    & 
     1194         add1_tmp,   & ! temp array for resizing address arrays 
     1195         add2_tmp  ! temp array for resizing address arrays 
     1196    ! 
     1197    ! temp array for resizing weight arrays 
     1198    ! 
     1199    REAL*8, DIMENSION(:,:), ALLOCATABLE :: wts_tmp    
     1200    ! 
     1201    !----------------------------------------------------------------------- 
     1202    !*** 
     1203    !*** allocate temporaries to hold original values 
     1204    !*** 
     1205    mxlinks = SIZE(grid1_add_map1) 
     1206    ALLOCATE (add1_tmp(mxlinks), add2_tmp(mxlinks), & 
     1207         wts_tmp(num_wts,mxlinks)) 
     1208 
     1209    add1_tmp = grid1_add_map1 
     1210    add2_tmp = grid2_add_map1 
     1211    wts_tmp  = wts_map1 
     1212 
     1213    !*** 
     1214    !*** deallocate originals and increment max_links then 
     1215    !*** reallocate arrays at new size 
     1216    !*** 
     1217 
     1218    DEALLOCATE (grid1_add_map1, grid2_add_map1, wts_map1) 
     1219    max_links_map1 = mxlinks + increment 
     1220    ALLOCATE (grid1_add_map1(max_links_map1),    & 
     1221         grid2_add_map1(max_links_map1),    & 
     1222         wts_map1(num_wts,max_links_map1)) 
     1223    !*** 
     1224    !*** restore original values from temp arrays and 
     1225    !*** deallocate temps 
     1226    !*** 
     1227    mxlinks = MIN(mxlinks, max_links_map1) 
     1228    grid1_add_map1(1:mxlinks) = add1_tmp (1:mxlinks) 
     1229    grid2_add_map1(1:mxlinks) = add2_tmp (1:mxlinks) 
     1230    wts_map1    (:,1:mxlinks) = wts_tmp(:,1:mxlinks) 
     1231    DEALLOCATE(add1_tmp, add2_tmp, wts_tmp) 
     1232 
     1233    !----------------------------------------------------------------------- 
     1234    ! 
     1235  END SUBROUTINE resize_remap_vars 
     1236  ! 
     1237  !************************************************************************ 
     1238  ! 
     1239END MODULE bilinear_interp 
     1240 
     1241!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     1242 
  • utils/tools/DOMAINcfg/src/dom_oce.F90

    r12414 r13204  
    4141   REAL(wp), PUBLIC ::   e3zps_min       !: miminum thickness for partial steps (meters) 
    4242   REAL(wp), PUBLIC ::   e3zps_rat       !: minimum thickness ration for partial steps 
     43   INTEGER , PUBLIC ::   nn_closea       !: =0 suppress closed sea/lake from the ORCA domain or not (=1) 
    4344 
    4445   INTEGER, PUBLIC :: nn_interp 
    45    CHARACTER(LEN=132), PUBLIC :: cn_topo 
     46   CHARACTER(LEN=200), PUBLIC :: cn_topo 
    4647   CHARACTER(LEN=132), PUBLIC :: cn_bath 
    4748   CHARACTER(LEN=132), PUBLIC :: cn_lon 
    4849   CHARACTER(LEN=132), PUBLIC :: cn_lat 
     50   REAL(wp), PUBLIC :: rn_scale 
    4951 
    5052   LOGICAL, PUBLIC ::   lzoom      =  .FALSE.   !: zoom flag 
     
    9395   !! ---------------------------- 
    9496   !                                   !!* Namelist namdom : time & space domain * 
    95    LOGICAL , PUBLIC ::   ln_linssh      !: =T  linear free surface ==>> model level are fixed in time 
    9697   LOGICAL , PUBLIC ::   ln_meshmask    !: =T  create a mesh-mask file (mesh_mask.nc) 
    9798   REAL(wp), PUBLIC ::   rn_isfhmin     !: threshold to discriminate grounded ice to floating ice 
     
    144145   INTEGER             , PUBLIC ::   narea            !: number for local area 
    145146   INTEGER             , PUBLIC ::   nbondi, nbondj   !: mark of i- and j-direction local boundaries 
    146    INTEGER, ALLOCATABLE, PUBLIC ::   nbondi_bdy(:)    !: mark i-direction local boundaries for BDY open boundaries 
    147    INTEGER, ALLOCATABLE, PUBLIC ::   nbondj_bdy(:)    !: mark j-direction local boundaries for BDY open boundaries 
    148    INTEGER, ALLOCATABLE, PUBLIC ::   nbondi_bdy_b(:)  !: mark i-direction of neighbours local boundaries for BDY open boundaries   
    149    INTEGER, ALLOCATABLE, PUBLIC ::   nbondj_bdy_b(:)  !: mark j-direction of neighbours local boundaries for BDY open boundaries   
    150147 
    151148   INTEGER, PUBLIC ::   npolj             !: north fold mark (0, 3 or 4) 
     
    279276#if defined key_agrif 
    280277   LOGICAL, PUBLIC, PARAMETER ::   lk_agrif = .TRUE.    !: agrif flag 
     278   LOGICAL, PUBLIC            ::   lk_south, lk_north, lk_west, lk_east !: Child grid boundaries (interpolation or not) 
    281279#else 
    282280   LOGICAL, PUBLIC, PARAMETER ::   lk_agrif = .FALSE.   !: agrif flag 
     
    326324         &      ff_f (jpi,jpj) ,    ff_t (jpi,jpj)                                     , STAT=ierr(3) ) 
    327325         ! 
    328       ALLOCATE( gdept_0(jpi,jpj,jpk) , gdepw_0(jpi,jpj,jpk) , STAT=ierr(4) ) 
     326      ! ALLOCATE( gdept_0(jpi,jpj,jpk) , gdepw_0(jpi,jpj,jpk) , gde3w_0(jpi,jpj,jpk) ,      & 
     327      !    &      gdept_b(jpi,jpj,jpk) , gdepw_b(jpi,jpj,jpk) ,                             & 
     328      !    &      gdept_n(jpi,jpj,jpk) , gdepw_n(jpi,jpj,jpk) , gde3w_n(jpi,jpj,jpk) , STAT=ierr(4) ) 
     329 
     330         ALLOCATE( gdept_0(jpi,jpj,jpk) , gdepw_0(jpi,jpj,jpk), STAT=ierr(4) ) 
     331 
     332         ! 
     333      ! 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) ,   & 
     334      !    &      e3t_b(jpi,jpj,jpk) , e3u_b(jpi,jpj,jpk) , e3v_b(jpi,jpj,jpk) ,                      e3w_b(jpi,jpj,jpk) ,   & 
     335      !    &      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) ,   & 
     336      !    &      e3t_a(jpi,jpj,jpk) , e3u_a(jpi,jpj,jpk) , e3v_a(jpi,jpj,jpk) ,                                             & 
     337      !    !                                                          ! 
     338      !    &      e3uw_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk) ,         & 
     339      !    &      e3uw_b(jpi,jpj,jpk) , e3vw_b(jpi,jpj,jpk) ,         & 
     340      !    &      e3uw_n(jpi,jpj,jpk) , e3vw_n(jpi,jpj,jpk) ,     STAT=ierr(5) ) 
     341 
     342         ! 
     343    !  ALLOCATE( ht_0(jpi,jpj) , hu_0(jpi,jpj) , hv_0(jpi,jpj) ,                                           & 
     344    !     &                      hu_b(jpi,jpj) , hv_b(jpi,jpj) , r1_hu_b(jpi,jpj) , r1_hv_b(jpi,jpj) ,     & 
     345    !     &      ht_n(jpi,jpj) , hu_n(jpi,jpj) , hv_n(jpi,jpj) , r1_hu_n(jpi,jpj) , r1_hv_n(jpi,jpj) ,     & 
     346    !     &                      hu_a(jpi,jpj) , hv_a(jpi,jpj) , r1_hu_a(jpi,jpj) , r1_hv_a(jpi,jpj) , STAT=ierr(6)  ) 
    329347         ! 
    330348      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) )                        
     349         &      e3uw_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk) , STAT=ierr(5) ) 
    332350         ! 
    333351      ALLOCATE( gdept_1d(jpk) , e3tp (jpi,jpj), e3wp(jpi,jpj) ,gdepw_1d(jpk) , e3t_1d(jpk) , e3w_1d(jpk) , STAT=ierr(6) ) 
    334352         ! 
    335       ALLOCATE( bathy(jpi,jpj),mbathy(jpi,jpj), tmask_i(jpi,jpj) , tmask_h(jpi,jpj) ,                        &  
     353      ALLOCATE( bathy(jpi,jpj),mbathy(jpi,jpj), tmask_i(jpi,jpj) , tmask_h(jpi,jpj) ,                        & 
    336354         &      ssmask (jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) ,     & 
    337355         &      mbkt   (jpi,jpj) , mbku   (jpi,jpj) , mbkv   (jpi,jpj) , STAT=ierr(7) ) 
     
    340358         &      risfdep(jpi,jpj) , mikv(jpi,jpj) , mikf(jpi,jpj) , STAT=ierr(8) ) 
    341359         ! 
    342       ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk) ,     &  
     360      ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk) ,     & 
    343361         &      vmask(jpi,jpj,jpk) , fmask(jpi,jpj,jpk) , wmask(jpi,jpj,jpk) , STAT=ierr(9) ) 
    344362         ! 
     
    351369      ALLOCATE( msk_opnsea  (jpi,jpj), msk_csundef (jpi,jpj),                        & 
    352370         &      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) ) 
    354       ! 
     371         &      msk_csgrpglo(jpi,jpj), msk_csgrprnf(jpi,jpj), msk_csgrpemp(jpi,jpj), STAT=ierr(11) )      ! 
    355372      dom_oce_alloc = MAXVAL(ierr) 
    356373      ! 
  • utils/tools/DOMAINcfg/src/domain.F90

    r12870 r13204  
    6262      !!              - 1D configuration, move Coriolis, u and v at T-point 
    6363      !!---------------------------------------------------------------------- 
     64      INTEGER ::   jk          ! dummy loop indices 
     65      INTEGER ::   iconf = 0   ! local integers 
     66      REAL(wp), POINTER, DIMENSION(:,:) ::   z1_hu_0, z1_hv_0 
     67      INTEGER , DIMENSION(jpi,jpj) ::   ik_top , ik_bot       ! top and bottom ocean level 
     68      !!---------------------------------------------------------------------- 
    6469      ! 
    6570      IF(lwp) THEN 
     
    7176      !                       !==  Reference coordinate system  ==! 
    7277      ! 
    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 
     78      CALL dom_nam               ! read namelist ( namrun, namdom ) 
     79                  !   CALL dom_clo               ! Closed seas and lake 
     80          
     81      CALL dom_hgr               ! Horizontal mesh 
     82      CALL dom_zgr( ik_top, ik_bot )  ! Vertical mesh and bathymetry 
     83      CALL dom_msk( ik_top, ik_bot )  ! Masks 
     84      ! 
    8285      ! 
    8386      CALL dom_ctl                  ! print extrema of masked scale factors 
    8487      !  
     88#if ! defined key_agrif 
    8589      CALL cfg_write                ! create the configuration file 
     90#endif 
    8691      ! 
    8792   END SUBROUTINE dom_init 
     
    98103      !!---------------------------------------------------------------------- 
    99104      USE ioipsl 
    100       NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,                 & 
    101                        nn_no   , cn_exp   , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl ,     & 
    102          &             nn_it000, nn_itend , nn_date0    , nn_time0     , nn_leapy  , nn_istate ,     & 
    103          &             nn_stock, nn_write , ln_mskland  , ln_clobber   , nn_chunksz, nn_euler  ,     & 
     105      NAMELIST/namrun/ cn_exp   ,    &           
     106         &             nn_it000, nn_itend , nn_date0    , nn_time0     , nn_leapy  ,     & 
     107         &             ln_mskland  , ln_clobber   , nn_chunksz,     & 
    104108         &             ln_cfmeta, ln_iscpl 
    105       NAMELIST/namdom/ nn_bathy, cn_topo, cn_bath, cn_lon, cn_lat, nn_interp,                        & 
    106          &             rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin,                       & 
    107          &             rn_atfp , rn_rdt   , ln_crs      , jphgr_msh ,                                & 
     109 
     110      NAMELIST/namdom/ ln_read_cfg, nn_bathy, cn_domcfg, cn_topo, cn_bath, cn_lon, cn_lat, rn_scale, nn_interp, & 
     111         &              rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin,            & 
     112         &             rn_atfp , rn_rdt   ,  ln_crs      , jphgr_msh ,                               & 
    108113         &             ppglam0, ppgphi0, ppe1_deg, ppe2_deg, ppe1_m, ppe2_m,                         & 
    109114         &             ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh,                  & 
    110115         &             ppa2, ppkth2, ppacr2 
    111  
    112  
    113116 
    114117      INTEGER  ::   ios                 ! Local integer output status for namelist read 
     
    137140 
    138141      cexper = cn_exp 
    139       nrstdt = nn_rstctl 
    140142      nit000 = nn_it000 
    141143      nitend = nn_itend 
    142144      ndate0 = nn_date0 
    143145      nleapy = nn_leapy 
    144       ninist = nn_istate 
    145       nstock = nn_stock 
    146       nstocklist = nn_stocklist 
    147       nwrite = nn_write 
    148       neuler = nn_euler 
    149       IF ( neuler == 1 .AND. .NOT. ln_rstart ) THEN 
    150          WRITE(ctmp1,*) 'ln_rstart =.FALSE., nn_euler is forced to 0 ' 
    151          CALL ctl_warn( ctmp1 ) 
    152          neuler = 0 
    153       ENDIF 
    154  
    155       !                             ! control of output frequency 
    156       IF ( nstock == 0 .OR. nstock > nitend ) THEN 
    157          WRITE(ctmp1,*) 'nstock = ', nstock, ' it is forced to ', nitend 
    158          CALL ctl_warn( ctmp1 ) 
    159          nstock = nitend 
    160       ENDIF 
    161       IF ( nwrite == 0 ) THEN 
    162          WRITE(ctmp1,*) 'nwrite = ', nwrite, ' it is forced to ', nitend 
    163          CALL ctl_warn( ctmp1 ) 
    164          nwrite = nitend 
    165       ENDIF 
    166  
    167  
    168  
    169  
    170       SELECT CASE ( nleapy )        ! Choose calendar for IOIPSL 
    171       CASE (  1 )  
    172          CALL ioconf_calendar('gregorian') 
    173          IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "gregorian", i.e. leap year' 
    174       CASE (  0 ) 
    175          CALL ioconf_calendar('noleap') 
    176          IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "noleap", i.e. no leap year' 
    177       CASE ( 30 ) 
    178          CALL ioconf_calendar('360d') 
    179          IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "360d", i.e. 360 days in a year' 
    180       END SELECT 
    181  
    182  
    183  
     146 
     147      !  
     148      cn_topo ='' 
     149      cn_bath ='' 
     150      cn_lon  ='' 
     151      cn_lat  ='' 
     152      rn_scale = 1. 
    184153 
    185154      REWIND( numnam_ref )              ! Namelist namdom in reference namelist : space & time domain (bathymetry, mesh, timestep) 
     
    193162      IF(lwm) WRITE ( numond, namdom ) 
    194163      ! 
     164 
     165 
     166 
    195167      IF(lwp) THEN 
    196168         WRITE(numout,*) 
     
    198170         WRITE(numout,*) '      flag read/compute bathymetry      nn_bathy     = ', nn_bathy 
    199171         IF( nn_bathy == 2 ) THEN 
    200             WRITE(numout,*) '      compute bathymetry from file      cn_topo      = ', cn_topo 
     172            WRITE(numout,*) '   compute bathymetry from file      cn_topo      = ' , cn_topo 
     173            WRITE(numout,*) '   bathymetry name in file           cn_bath      = ' , cn_bath 
     174            WRITE(numout,*) '   longitude name in file            cn_lon       = ' , cn_lon 
     175            WRITE(numout,*) '   latitude  name in file            cn_lat       = ' , cn_lat 
     176            WRITE(numout,*) '   bathmetry scale factor            rn_scale     = ' , rn_scale  
    201177         ENDIF    
    202178         WRITE(numout,*) '      Depth (if =0 bathy=jpkm1)         rn_bathy     = ', rn_bathy 
     
    255231      !!---------------------------------------------------------------------- 
    256232      ! 
     233#undef CHECK_DOM 
     234#ifdef CHECK_DOM 
    257235      IF(lk_mpp) THEN 
    258236         CALL mpp_minloc( 'dom_ctl', e1t(:,:), tmask_i(:,:), ze1min, iloc ) 
     
    292270         WRITE(numout,"(14x,'e2t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2min, iimi2, ijmi2 
    293271      ENDIF 
     272#endif 
    294273      ! 
    295274      ! check that all processes are still there... If some process have an error, 
     
    426405      CALL iom_rstput( 0, 0, inum, 'top_level'    , REAL( mikt, wp )*ssmask , ktype = jp_i4 )   ! nb of ocean T-points (ISF) 
    427406      CALL iom_rstput( 0, 0, inum, 'isf_draft'    , risfdep , ktype = jp_r8 ) 
    428       CALL iom_rstput( 0, 0, inum, 'bathy_metry'  , bathy   , ktype = jp_r8 ) 
     407      DO jj = 1,jpj 
     408         DO ji = 1,jpi 
     409            z2d (ji,jj) = SUM ( e3t_0(ji,jj, 1:mbkt(ji,jj) ) ) * ssmask(ji,jj)  
     410         END DO 
     411      END DO 
     412      CALL iom_rstput( 0, 0, inum, 'bathy_metry'   , z2d , ktype = jp_r8 ) 
    429413      ! 
    430414      !                              !== closed sea ==! 
  • utils/tools/DOMAINcfg/src/dombat.F90

    r12414 r13204  
    22 
    33   USE dom_oce           ! ocean domain 
    4 !   USE closea            ! closed seas 
    54   ! 
    65   USE in_out_manager    ! I/O manager 
     
    87   USE lbclnk            ! ocean lateral boundary conditions (or mpp link) 
    98   USE lib_mpp           ! distributed memory computing library 
     9#if defined key_agrif 
    1010   USE agrif_modutil 
     11   USE agrif_parameters 
     12#endif    
    1113   USE bilinear_interp 
    1214 
     
    2022   SUBROUTINE dom_bat 
    2123 
    22       INTEGER :: inum, isize, jsize, id, ji, jj 
     24      INTEGER :: inum, id, ji, jj,ji1,jj1 
     25      INTEGER :: iimin,iimax,jjmin,jjmax 
    2326      INTEGER :: tabdim1, tabdim2, nxhr, nyhr, nxyhr 
    2427      INTEGER, DIMENSION(2) :: ddims 
    2528      INTEGER, DIMENSION(3) :: status 
    26       INTEGER, DIMENSION(:,:), ALLOCATABLE :: trouble_points ,  vardep,mask_oce 
    27       INTEGER :: iimin,iimax,jjmin,jjmax 
    2829      INTEGER, DIMENSION(1) :: i_min,i_max 
    29       INTEGER, DIMENSION(1) ::j_min,j_max 
    30       REAL(wp), DIMENSION(jpi,jpj) :: myglamf 
    31       INTEGER,DIMENSION(:)  ,POINTER     ::   src_add,dst_add => NULL() 
    32       REAL(wp), DIMENSION(:)  ,ALLOCATABLE ::   vardep1d, lon_new1D,lat_new1D  
    33       REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bathy_new, lat_new, lon_new, bathy_test 
    34       REAL(wp), DIMENSION(:,:), ALLOCATABLE :: coarselon, coarselat, coarsebathy 
    35       REAL(wp) :: Cell_lonmin, Cell_lonmax, Cell_latmin, Cell_latmax 
     30      INTEGER, DIMENSION(1) :: j_min,j_max 
     31      INTEGER, DIMENSION(:)  ,ALLOCATABLE     ::   src_add,dst_add  
     32      INTEGER, DIMENSION(:,:), ALLOCATABLE :: trouble_points , vardep,mask_oce 
     33      REAL(wp) ,DIMENSION(jpi,jpj):: Cell_lonmin, Cell_lonmax, Cell_latmin, Cell_latmax 
     34      REAL(wp) ::zdel 
     35      REAL(wp), DIMENSION(:)  , ALLOCATABLE ::   lon_new1D , lat_new1D, vardep1d 
     36      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   coarselon, coarselat, coarsebathy, bathy_test 
     37      REAL(wp), DIMENSION(:,:),ALLOCATABLE     ::   matrix,interpdata  
     38      LOGICAL :: lonlat_2D, ln_pacifiq 
    3639      LOGICAL :: identical_grids 
    3740      LOGICAL, DIMENSION(:,:), ALLOCATABLE     ::   masksrc 
    38       REAL*8, DIMENSION(:,:),POINTER     ::   matrix,interpdata => NULL()  
    39       LOGICAL :: lonlat_2D 
    40  
     41      REAL(wp), DIMENSION(jpi,jpj) :: zglamt, zgphi, zglamu, zglamv, zglamf 
     42      REAL(wp) :: zshift 
     43  
    4144      CHARACTER(32) :: bathyfile, bathyname, lonname,latname        
    42  
    43      lonlat_2D=.false. 
    4445 
    4546      bathyfile=TRIM(cn_topo) 
     
    4950    
    5051      CALL iom_open( bathyfile, inum, lagrif=.FALSE. ) 
     52       
     53      ! check if lon/lat are 2D arrays 
     54      id = iom_varid( inum, lonname, ddims ) 
     55      IF (ddims(2)==0) THEN 
     56         lonlat_2D = .FALSE. 
     57      ELSE 
     58         lonlat_2D = .TRUE. 
     59      ENDIF    
     60       
    5161      id = iom_varid( inum, bathyname, ddims ) 
    52        
     62      ln_pacifiq = .FALSE. 
     63      zglamt(:,:) = glamt(:,:) 
     64      zglamu(:,:) = glamu(:,:) 
     65      zglamv(:,:) = glamv(:,:) 
     66      zglamf(:,:) = glamf(:,:) 
     67 
     68      IF( glamt(1,1) .GT. glamt(jpi,jpj) ) ln_pacifiq =.false.  
     69 
     70      zshift = 0. 
     71      IF( ln_pacifiq  ) THEN 
     72         zshift = 0.!Abs(minval(glamt)) +0.1  
     73         WHERE ( glamt < 0 ) 
     74            zglamt = zglamt + zshift + 360. 
     75         END WHERE 
     76         WHERE ( glamu < 0 ) 
     77            zglamu = zglamu + zshift +360. 
     78         END WHERE 
     79         WHERE ( glamv < 0 ) 
     80            zglamv = zglamv + zshift +360. 
     81         END WHERE 
     82         WHERE ( glamf < 0 ) 
     83            zglamf = zglamf + zshift +360. 
     84         END WHERE 
     85      ENDIF 
     86 
    5387      status=-1 
    54       ALLOCATE(lon_new  (ddims(1),ddims(2)), STAT=status(1))  
    55       ALLOCATE(lat_new  (ddims(1),ddims(2)), STAT=status(2))  
    56       ALLOCATE(bathy_new(ddims(1),ddims(2)), STAT=status(3))  
    57       IF( sum(status) /= 0 )   CALL ctl_stop( 'STOP', 'dom_bat : unable to allocate arrays' ) 
    5888 
    5989      IF (lonlat_2D) THEN 
    60         CALL iom_get  ( inum, jpdom_unknown, lonname, lon_new ) 
    61         CALL iom_get  ( inum, jpdom_unknown, latname, lat_new ) 
     90         ! here implicitly it's old topo database (orca format) 
     91         ALLOCATE(coarselon  (ddims(1),ddims(2)), STAT=status(1))  
     92         ALLOCATE(coarselat  (ddims(1),ddims(2)), STAT=status(2))  
     93         ALLOCATE(coarsebathy(ddims(1),ddims(2)), STAT=status(3))  
     94         IF( sum(status) /= 0 )   CALL ctl_stop( 'STOP', 'dom_bat : unable to allocate arrays' ) 
     95         CALL iom_get  ( inum, jpdom_unknown, lonname, coarselon ) 
     96         CALL iom_get  ( inum, jpdom_unknown, latname, coarselat ) 
     97         CALL iom_get  ( inum, jpdom_unknown, bathyname, coarsebathy ) 
     98         CALL iom_close (inum) 
     99         IF( ln_pacifiq ) THEN 
     100       !     WHERE(coarselon < 0.00001)  
     101                coarselon = Coarselon + zshift 
     102       !      END WHERE 
     103         ENDIF      
     104         ! equivalent to new database 
    62105      ELSE 
    63         ALLOCATE(lon_new1D(ddims(1)), lat_new1D(ddims(2))) 
    64         CALL iom_get  ( inum, jpdom_unknown, lonname, lon_new1D ) 
    65         CALL iom_get  ( inum, jpdom_unknown, latname, lat_new1D ) 
    66         DO ji=1, ddims(1) 
    67            lon_new(ji,:)=lon_new1D(ji) 
    68         ENDDO               
    69         DO ji=1, ddims(2) 
    70            lat_new(:,ji)=lat_new1D(ji) 
    71         ENDDO               
     106         ALLOCATE(lon_new1D(ddims(1)), lat_new1D(ddims(2))) 
     107         CALL iom_get  ( inum, jpdom_unknown, lonname, lon_new1D ) 
     108         CALL iom_get  ( inum, jpdom_unknown, latname, lat_new1D ) 
     109         IF( ln_pacifiq ) THEN 
     110            WHERE(lon_new1D < 0.00001)  
     111                lon_new1D = lon_new1D +360.!zshift 
     112             END WHERE 
     113         ENDIF                 
     114         zdel =  0.00    
     115         IF( MAXVAL(zglamf) > 180 + zshift ) THEN   
     116            !           
     117         !   WHERE( lon_new1D < 0 ) 
     118         !       lon_new1D = lon_new1D + 360. 
     119         !   END WHERE 
     120            !      
     121            i_min = MAXLOC(lon_new1D,mask = lon_new1D < MINVAL(zglamf(1:jpi-1,1:jpj-1)) ) 
     122            i_max = MINLOC(lon_new1D,mask = lon_new1D > MAXVAL(zglamf(1:jpi-1,1:jpj-1)) )                    
     123            j_min = MAXLOC(lat_new1D,mask = lat_new1D < MINVAL( gphif(1:jpi-1,1:jpj-1)) ) 
     124            j_max = MINLOC(lat_new1D,mask = lat_new1D > MAXVAL( gphif(1:jpi-1,1:jpj-1)) ) 
     125            ! 
     126            tabdim1 = ( SIZE(lon_new1D) - i_min(1) + 1 ) + i_max(1)                    
     127            !           
     128            IF(j_min(1)-2 >= 1 .AND. j_max(1)+3 <= SIZE(lat_new1D,1) ) THEN 
     129               j_min(1) = j_min(1) - 2 
     130               j_max(1) = j_max(1)+ 3 
     131            ENDIF 
     132            tabdim2 = j_max(1) - j_min(1) + 1 
     133            ! 
     134            ALLOCATE(coarselon  (tabdim1,tabdim2), STAT=status(1)) 
     135            ALLOCATE(coarselat  (tabdim1,tabdim2), STAT=status(2)) 
     136            ALLOCATE(Coarsebathy(tabdim1,tabdim2), STAT=status(3))  
     137 
     138            IF( SUM(status) /= 0 ) CALL ctl_stop( 'STOP', 'dom_bat : unable to allocate arrays' )           
     139            ! 
     140            DO ji = 1,tabdim1 
     141               coarselat(ji,:) = lat_new1D(j_min(1):j_max(1)) 
     142            END DO 
     143            
     144            ! 
     145            DO jj = 1, tabdim2                                  
     146               coarselon(1:SIZE(lon_new1D)-i_min(1)+1      ,jj) = lon_new1D(i_min(1):SIZE(lon_new1D)) 
     147               coarselon(2+SIZE(lon_new1D)-i_min(1):tabdim1,jj) = lon_new1D(1:i_max(1))    
     148            END DO 
     149            !  
     150            CALL iom_get(inum, jpdom_unknown, bathyname,coarsebathy(1:SIZE(lon_new1D)-i_min(1)+1,:), & 
     151            kstart=(/i_min(1),j_min(1)/), kcount=(/SIZE(lon_new1D)-i_min(1)+1,tabdim2/))   ! +1? 
     152            ! 
     153            CALL iom_get(inum, jpdom_unknown, bathyname,coarsebathy(2+SIZE(lon_new1D)-i_min(1):tabdim1,:), & 
     154            kstart=(/1,j_min(1)/),kcount=(/i_max(1),tabdim2/))                 
     155            ! 
     156         ELSE 
     157          !  WHERE( lon_new1D > (180. + zshift) )   lon_new1D = lon_new1D - 360. 
     158            i_min = MAXLOC(lon_new1D,mask = lon_new1D < MINVAL(zglamf)-zdel) 
     159            i_max = MINLOC(lon_new1D,mask = lon_new1D > MAXVAL(zglamf)+zdel) 
     160            j_min = MAXLOC(lat_new1D,mask = lat_new1D < MINVAL(gphif)-zdel) 
     161            j_max = MINLOC(lat_new1D,mask = lat_new1D > MAXVAL(gphif)+zdel) 
     162          
     163            i_min(1)=max(i_min(1),1) 
     164            !       
     165            IF(i_min(1)-2 >= 1 .AND. i_max(1)+3 <= SIZE(lon_new1D,1) ) THEN 
     166               i_min(1) = i_min(1)-2 
     167               i_max(1) = i_max(1)+3 
     168            ENDIF 
     169            tabdim1 = i_max(1) - i_min(1) + 1 
     170            ! 
     171            IF(j_min(1)-2 >= 1 .AND. j_max(1)+3 <= SIZE(lat_new1D,1) ) THEN 
     172               j_min(1) = j_min(1)-2 
     173               j_max(1) = j_max(1)+3 
     174            ENDIF 
     175            tabdim2 = j_max(1) - j_min(1) + 1 
     176            ! 
     177            ALLOCATE(coarselon  (tabdim1,tabdim2), STAT=status(1))  
     178            ALLOCATE(coarselat  (tabdim1,tabdim2), STAT=status(2))  
     179            ALLOCATE(coarsebathy(tabdim1,tabdim2), STAT=status(3))  
     180 
     181            IF( SUM(status) /= 0 ) CALL ctl_stop( 'STOP', 'dom_bat : unable to allocate arrays' )   
     182            ! 
     183            DO jj = 1,tabdim2 
     184               coarselon(:,jj) = lon_new1D(i_min(1):i_max(1)) 
     185            END DO 
     186            !       
     187            DO ji = 1,tabdim1 
     188              coarselat(ji,:) = lat_new1D(j_min(1):j_max(1))  
     189            END DO 
     190            ! 
     191            CALL iom_get(inum,jpdom_unknown,bathyname,coarsebathy, & 
     192            &                  kstart=(/i_min(1),j_min(1)/),kcount=(/tabdim1,tabdim2/)) 
     193            ! 
     194         ENDIF   ! > 180 
     195     
     196         DEALLOCATE(lon_new1D) ; DEALLOCATE(lat_new1D) 
     197         CALL iom_close (inum) 
     198          
     199         coarsebathy = coarsebathy *rn_scale 
     200        ! reset land to 0  
     201         WHERE (coarsebathy < 0.)         
     202          coarsebathy=0. 
     203         ENDWHERE  
     204  
     205      ENDIF   ! external 
     206 
     207      IF(lwp) THEN 
     208         WRITE(numout,*) 'Interpolation of high resolution bathymetry on child grid' 
     209         IF( nn_interp == 0 ) THEN 
     210            WRITE(numout,*) 'Arithmetic average ...' 
     211         ELSE IF( nn_interp == 1 ) THEN 
     212            WRITE(numout,*) 'Median average ...' 
     213         ELSE IF( nn_interp == 2 ) THEN      
     214            WRITE(numout,*) 'Bilinear interpolation ...' 
     215         ELSE      
     216            WRITE(*,*) 'bad value for nn_interp variable ( must be 0,1 or 2 )' 
     217            STOP  
     218         ENDIF 
    72219      ENDIF   
    73       CALL iom_get  ( inum, jpdom_unknown, bathyname, bathy_new ) 
    74       CALL iom_close (inum) 
    75       WHERE (bathy_new > 0.)         
    76         bathy_new=0. 
    77       ENDWHERE  
    78       bathy_new=-bathy_new  
    79  
    80        ! Eventually add here a pre-selection of the area (coarselon/lat) 
    81        i_min=10  ; j_min=10 
    82        i_max= ddims(1)-10 ; j_max=ddims(2)-10  
    83  
    84        tabdim1 = i_max(1) -  i_min(1) + 1  
    85        tabdim2 = j_max(1) - j_min(1) + 1 
    86        ! 
    87  
    88        ALLOCATE(coarselon(tabdim1,tabdim2)) 
    89        ALLOCATE(coarselat(tabdim1,tabdim2)) 
    90        ALLOCATE(coarsebathy(tabdim1,tabdim2))           
    91        ! 
    92        WHERE( lon_new < 0. )   lon_new = lon_new + 360. 
    93        myglamf=glamf 
    94        WHERE( myglamf < 0. )   myglamf = myglamf + 360. 
    95  
    96        coarselat(:,:)   =  lat_new  (i_min(1):i_max(1), j_min(1):j_max(1)) 
    97        coarselon  (:,:) =  lon_new  (i_min(1):i_max(1), j_min(1):j_max(1)) 
    98        coarsebathy(:,:) =  bathy_new(i_min(1):i_max(1), j_min(1):j_max(1)) 
    99  
    100        IF( nn_interp == 0 .OR. nn_interp == 1 ) THEN   ! arithmetic or median averages 
    101         !                                                            ! -----------------------------  
    102         ALLOCATE(trouble_points(jpi,jpj)) 
    103         trouble_points(:,:) = 0 
    104         !                        
    105         DO jj = 2, jpj 
    106            DO ji = 2, jpi 
    107               !        
    108               ! fine grid cell extension                
    109               Cell_lonmin = MIN(myglamf(ji-1,jj-1),myglamf(ji,jj-1),myglamf(ji,jj),myglamf(ji-1,jj)) 
    110               Cell_lonmax = MAX(myglamf(ji-1,jj-1),myglamf(ji,jj-1),myglamf(ji,jj),myglamf(ji-1,jj)) 
    111               Cell_latmin = MIN(gphif(ji-1,jj-1),gphif(ji,jj-1),gphif(ji,jj),gphif(ji-1,jj)) 
    112               Cell_latmax = MAX(gphif(ji-1,jj-1),gphif(ji,jj-1),gphif(ji,jj),gphif(ji-1,jj))  
    113               !                
    114               ! look for points in G0 (bathy dataset) contained in the fine grid cells   
    115               iimin = 1 
    116               DO WHILE( coarselon(iimin,1) < Cell_lonmin )  
    117                  iimin = iimin + 1 
    118               ENDDO 
    119               !                
    120               jjmin = 1 
    121               DO WHILE( coarselat(iimin,jjmin) < Cell_latmin )  
    122                  jjmin = jjmin + 1 
    123               ENDDO 
    124              !                 
    125               iimax = iimin  
    126               DO WHILE( coarselon(iimax,1) <= Cell_lonmax )  
    127                  iimax = iimax + 1 
    128                  iimax = MIN( iimax,SIZE(coarsebathy,1)) 
    129               ENDDO 
    130               !                                
    131               jjmax = jjmin  
    132               DO WHILE( coarselat(iimax,jjmax) <= Cell_latmax )  
    133                  jjmax = jjmax + 1 
    134                  jjmax = MIN( jjmax,SIZE(coarsebathy,2)) 
    135               ENDDO 
    136               ! 
    137               IF( .NOT. Agrif_Root() ) THEN 
    138                  iimax = iimax-1 
    139                  jjmax = jjmax-1 
    140               ELSE 
    141                  iimax = MAX(iimin,iimax-1) 
    142                  jjmax = MAX(jjmin,jjmax-1) 
    143               ENDIF 
    144               !                
    145               iimin = MAX( iimin,1 ) 
    146               jjmin = MAX( jjmin,1 ) 
    147               iimax = MIN( iimax,SIZE(coarsebathy,1)) 
    148               jjmax = MIN( jjmax,SIZE(coarsebathy,2)) 
    149  
    150               nxhr = iimax - iimin + 1 
    151               nyhr = jjmax - jjmin + 1                     
    152  
    153                  
    154               IF( nxhr == 0 .OR. nyhr == 0 ) THEN 
    155                  ! 
    156                  trouble_points(ji,jj) = 1 
    157                  ! 
    158               ELSE 
    159                  ! 
    160                  ALLOCATE( vardep(nxhr,nyhr), mask_oce(nxhr,nyhr) ) 
    161                  vardep(:,:) = coarsebathy(iimin:iimax,jjmin:jjmax) 
    162                  ! 
    163                  WHERE( vardep(:,:) .GT. 0. )   ;   mask_oce = 1 ; 
    164                  ELSEWHERE                      ;   mask_oce = 0 ; 
    165                  ENDWHERE 
    166                  ! 
    167                  nxyhr = nxhr*nyhr 
    168                  IF( SUM(mask_oce) < 0.5*(nxyhr) ) THEN ! if more than half of the points are on land then bathy fine = 0 
    169                     bathy(ji,jj) = 0. 
    170                  ELSE 
    171                     IF( nn_interp == 0 ) THEN     ! Arithmetic average 
    172                       bathy(ji,jj) = SUM( vardep(:,:) * mask_oce(:,:) ) / SUM( mask_oce(:,:) ) 
    173                     ELSE                                  ! Median average 
    174                        ALLOCATE(vardep1d(nxyhr)) 
    175                        vardep1d = RESHAPE(vardep,(/ nxyhr /) ) 
    176                        !!CALL ssort(vardep1d,nxyhr) 
    177                        CALL quicksort(vardep1d,1,nxyhr) 
    178                        ! 
    179                        ! Calculate median 
    180                        IF (MOD(nxyhr,2) .NE. 0) THEN 
    181                           bathy(ji,jj) = vardep1d( nxyhr/2 + 1 ) 
    182                        ELSE 
    183                           bathy(ji,jj) = 0.5 * ( vardep1d(nxyhr/2) + vardep1d(nxyhr/2+1) ) 
    184                        END IF 
    185                        DEALLOCATE(vardep1d)    
    186                     ENDIF 
    187                  ENDIF 
    188                  DEALLOCATE (mask_oce,vardep) 
    189                  ! 
    190               ENDIF 
    191            ENDDO 
    192         ENDDO 
    193  
    194         IF( SUM( trouble_points ) > 0 ) THEN 
    195            PRINT*,'too much empty cells, proceed to bilinear interpolation' 
    196            nn_interp = 2 
    197            stop 
    198         ENDIF            
    199      ENDIF 
    200  
    201 #undef MYTEST 
    202 #ifdef MYTEST 
    203          IF( nn_interp == 2) THEN                        ! Bilinear interpolation 
    204         !                                                    ! -----------------------------  
    205         identical_grids = .FALSE. 
    206  
    207         IF( SIZE(coarselat,1) == jpi .AND. SIZE(coarselat,2) == jpj  .AND.   & 
    208             SIZE(coarselon,1) == jpj .AND. SIZE(coarselon,2) == jpj ) THEN 
    209            IF( MAXVAL( ABS(coarselat(:,:)- gphit(:,:)) ) < 0.0001 .AND.   & 
    210                MAXVAL( ABS(coarselon(:,:)- glamt(:,:)) ) < 0.0001 ) THEN 
    211               PRINT*,'same grid between ', cn_topo,' and child domain'     
    212               bathy = bathy_new 
    213               identical_grids = .TRUE.                           
    214            ENDIF 
    215         ENDIF 
    216  
    217         IF( .NOT. identical_grids ) THEN  
    218  
    219            ALLOCATE(masksrc(SIZE(bathy_new,1),SIZE(bathy_new,2))) 
    220            ALLOCATE(bathy_test(jpi,jpj)) 
    221            ! 
    222            !Where(G0%bathy_meter.le.0.00001)  
    223            !  masksrc = .false. 
    224            !ElseWhere 
    225               masksrc = .TRUE. 
    226            !End where                        
    227            !             
    228            ! compute remapping matrix thanks to SCRIP package             
    229            CALL get_remap_matrix(coarselat,gphit,coarselon,glamt,masksrc,matrix,src_add,dst_add) 
    230            CALL make_remap(bathy_new,bathy_test,jpi,jpj,matrix,src_add,dst_add)   
    231            !                                   
    232            bathy = bathy_test                
    233            !             
    234            DEALLOCATE(masksrc) 
    235            DEALLOCATE(bathy_test)  
    236  
    237         ENDIF 
     220      bathy(:,:) = 0. 
     221      ! 
     222      !------------------------------------ 
     223      !MEDIAN AVERAGE or ARITHMETIC AVERAGE 
     224      !------------------------------------ 
     225      ! 
     226      IF( nn_interp == 0 .OR. nn_interp == 1 ) THEN  
     227         ! 
     228         ALLOCATE(trouble_points(jpi,jpj)) 
     229         trouble_points = 0 
     230         ! 
     231         !  POINT DETECTION 
     232         ! 
     233         !                        
     234         DO jj = 1,jpj 
     235            DO ji = 1,jpi 
     236               !      
     237               ! FINE GRID CELLS DEFINITION   
     238               ji1=max(ji-1,1) 
     239               jj1=max(jj-1,1)              
     240              
     241               ! 
     242               Cell_lonmin(ji,jj) = MIN(zglamf(ji1,jj1),zglamf(ji,jj1),zglamf(ji,jj),zglamf(ji1,jj)) 
     243               Cell_lonmax(ji,jj) = MAX(zglamf(ji1,jj1),zglamf(ji,jj1),zglamf(ji,jj),zglamf(ji1,jj)) 
     244               Cell_latmin(ji,jj) = MIN(gphif(ji1,jj1),gphif(ji,jj1),gphif(ji,jj),gphif(ji1,jj)) 
     245               Cell_latmax(ji,jj) = MAX(gphif(ji1,jj1),gphif(ji,jj1),gphif(ji,jj),gphif(ji1,jj)) 
     246               IF( ABS(Cell_lonmax(ji,jj) - Cell_lonmin(ji,jj) ) > 180 ) THEN 
     247                    zdel = Cell_lonmax(ji,jj) 
     248                    Cell_lonmax(ji,jj) = Cell_lonmin(ji,jj) 
     249                    Cell_lonmin(ji,jj) = zdel-360 
     250               ENDIF 
     251               !                
     252               ! SEARCH FOR ALL POINTS CONTAINED IN THIS CELL 
     253               !     
     254     !      ENDDO 
     255     !    ENDDO    
     256      !   CALL lbc_lnk( 'dom_bat', Cell_lonmin, 'T', 1. ) 
     257      !   CALL lbc_lnk( 'dom_bat', Cell_lonmax, 'T', 1. ) 
     258      !   CALL lbc_lnk( 'dom_bat', Cell_latmin, 'T', 1. ) 
     259      !   CALL lbc_lnk( 'dom_bat', Cell_latmax, 'T', 1. ) 
     260 
     261 
     262      !   DO jj = 2,jpj 
     263      !      DO ji = 2,jpi 
     264               iimin = 1 
     265               DO WHILE( coarselon(iimin,1) < Cell_lonmin(ji,jj) )  
     266                  iimin = iimin + 1 
     267             !     IF ( iimin .LE. 1 ) THEN 
     268             !     iimin = 1 
     269             !     EXIT 
     270             !     ENDIF 
     271               ENDDO 
     272               !                
     273               jjmin = 1 
     274               DO WHILE( coarselat(iimin,jjmin) < Cell_latmin(ji,jj) )  
     275                  jjmin = jjmin + 1 
     276             !     IF ( iimin .LE. 1 ) THEN 
     277             !     iimin = 1 
     278             !     EXIT 
     279             !     ENDIF 
     280               ENDDO 
     281               jjmin=max(1,jjmin) 
     282               !                 
     283               iimax = iimin  
     284               DO WHILE( coarselon(iimax,1)<= Cell_lonmax(ji,jj) )  
     285                  iimax = iimax + 1 
     286                  IF ( iimax .GE. SIZE(coarsebathy,1) ) THEN 
     287                  iimax = MIN( iimax,SIZE(coarsebathy,1)) 
     288                  EXIT 
     289                ENDIF 
     290               ENDDO 
     291               !                                
     292               jjmax = jjmin  
     293               DO WHILE( coarselat(iimax,jjmax) <= Cell_latmax(ji,jj) )  
     294                  jjmax = jjmax + 1 
     295                  IF ( jjmax .GE. SIZE(coarsebathy,2) ) THEN 
     296                  jjmax = MIN( jjmax,SIZE(coarsebathy,2)) 
     297                  EXIT 
     298                ENDIF 
     299               ENDDO 
     300               ! 
     301            !   iimax = iimax-1 
     302            !   jjmax = jjmax-1 
     303               !                
     304               iimin = MAX( iimin,1 ) 
     305               jjmin = MAX( jjmin,1 ) 
     306               iimax = MIN( iimax,SIZE(coarsebathy,1)) 
     307               jjmax = MIN( jjmax,SIZE(coarsebathy,2)) 
     308 
     309               nxhr = iimax - iimin + 1 
     310               nyhr = jjmax - jjmin + 1    
     311 
     312 
     313   
     314               IF( nxhr == 0 .OR. nyhr == 0 ) THEN 
     315                  trouble_points(ji,jj) = 1 
     316               ELSE 
     317                  ALLOCATE( vardep(nxhr,nyhr) ) 
     318                  ALLOCATE( mask_oce(nxhr,nyhr) ) 
     319                  mask_oce = 0        
     320                  vardep(:,:) = coarsebathy(iimin:iimax,jjmin:jjmax) 
     321                  WHERE( vardep(:,:) .GT. 0. )  
     322                     mask_oce = 1 
     323                  ENDWHERE 
     324                  nxyhr = nxhr*nyhr 
     325!                 IF( SUM(mask_oce) == 0 ) THEN 
     326                   IF( SUM(mask_oce) < 0.5*(nxhr*nyhr) ) THEN 
     327                     bathy(ji,jj) = 0. 
     328                   ELSE 
     329                     IF( nn_interp == 0 ) THEN 
     330                        ! Arithmetic average                    
     331                        bathy(ji,jj) = SUM (vardep(:,:)*mask_oce(:,:))/SUM(mask_oce) 
     332                     ELSE 
     333                        ! Median average        
     334                        ! 
     335                        ALLOCATE(vardep1d(nxhr*nyhr)) 
     336                        vardep1d = RESHAPE(vardep,(/ nxhr*nyhr /) ) 
     337                       ! CALL ssort(vardep1d,nxhr*nyhr) 
     338                        CALL quicksort(vardep1d,1,nxhr*nyhr) 
     339                        ! 
     340                        ! Calculate median 
     341                        ! 
     342                        IF (MOD(SUM(mask_oce),2) .NE. 0) THEN 
     343                           bathy(ji,jj) = vardep1d( nxyhr/2 + 1 ) 
     344                        ELSE 
     345                           bathy(ji,jj) =0.5 * ( vardep1d(nxyhr/2) + vardep1d(nxyhr/2+1) ) 
     346                        END IF 
     347                        DEALLOCATE(vardep1d)    
     348                     ENDIF 
     349                  ENDIF 
     350                  DEALLOCATE (mask_oce,vardep) 
     351               ENDIF 
     352            ENDDO 
     353         ENDDO 
     354 
     355         IF( SUM( trouble_points ) > 0 ) THEN 
     356            CALL ctl_warn ('too much empty cells, proceed to bilinear interpolation') 
     357            nn_interp = 2 
     358         ENDIF 
     359      ENDIF 
     360 
     361     ! 
     362     ! create logical array masksrc 
     363     ! 
     364      IF( nn_interp == 2) THEN  
     365         ! 
     366         identical_grids = .FALSE. 
     367 
     368# ifdef key_agrif 
     369         IF( Agrif_Parent(jpi) == jpi  .AND.   & 
     370             Agrif_Parent(jpj) == jpj  ) THEN 
     371            IF( MAXVAL( ABS(coarselat(:,:)- gphit(:,:)) ) < 0.0001 .AND.   & 
     372                 MAXVAL( ABS(coarselon(:,:)- gphit(:,:)) ) < 0.0001 ) THEN 
     373               bathy(:,:)  = coarsebathy(:,:)  
     374                IF(lwp) WRITE(numout,*) 'same grid between ', bathyname,' and child domain'                    
     375               identical_grids = .TRUE.                           
     376            ENDIF 
     377         ENDIF 
     378# endif 
     379         IF( .NOT. identical_grids ) THEN   
     380            ALLOCATE(masksrc(SIZE(coarsebathy,1),SIZE(coarsebathy,2))) 
     381            ALLOCATE(bathy_test(jpi,jpj)) 
     382            ! 
     383            !                    Where(G0%bathy_meter.le.0.00001)  
     384            !                  masksrc = .false. 
     385            !              ElseWhere 
     386            ! 
     387            masksrc = .TRUE. 
     388            ! 
     389            !              End where                        
     390            !             
     391            ! compute remapping matrix thanks to SCRIP package             
     392            !                                   
     393            CALL get_remap_matrix(coarselat,gphit,   & 
     394                 coarselon,zglamt,   & 
     395                 masksrc,matrix,src_add,dst_add) 
     396            CALL make_remap(coarsebathy,bathy_test,jpi,jpj, & 
     397                 matrix,src_add,dst_add)   
     398            !                                   
     399            bathy= bathy_test                
     400            !             
     401            DEALLOCATE(masksrc) 
     402            DEALLOCATE(bathy_test)  
     403         ENDIF 
    238404        !             
    239      ENDIF 
    240 #endif  
     405      ENDIF 
     406      CALL lbc_lnk( 'dom_bat', bathy, 'T', 1. ) 
     407 
     408       ! Correct South and North 
     409#if defined key_agrif 
     410      IF( lk_south ) THEN   
     411         IF( (nbondj == -1).OR.(nbondj == 2) ) THEN 
     412           bathy(:,1)=bathy(:,2) 
     413         ENDIF 
     414      ELSE 
     415            bathy(:,1) = 0. 
     416      ENDIF 
     417#else 
     418      IF ((nbondj == -1).OR.(nbondj == 2)) THEN 
     419            bathy(:,1)=bathy(:,2) 
     420      ENDIF 
     421#endif 
     422 
     423      IF ((nbondj == 1).OR.(nbondj == 2)) THEN 
     424         bathy(:,jpj)=bathy(:,jpj-1) 
     425      ENDIF 
     426 
     427      ! Correct West and East 
     428      IF (jperio /=1) THEN 
     429         IF ((nbondi == -1).OR.(nbondi == 2)) THEN 
     430            bathy(1,:)=bathy(2,:) 
     431         ENDIF 
     432         IF ((nbondi == 1).OR.(nbondi == 2)) THEN 
     433         bathy(jpi,:)=bathy(jpi-1,:) 
     434         ENDIF 
     435      ENDIF 
     436 
     437 
    241438   END SUBROUTINE dom_bat 
    242439 
  • utils/tools/DOMAINcfg/src/dombat_util.f90

    r12414 r13204  
    1 ! 
    21MODULE agrif_modutil 
    32  ! 
  • utils/tools/DOMAINcfg/src/domcfg.f90

    r12414 r13204  
    9999         WRITE(numout,*) 
    100100         WRITE(numout,*) '   conversion from local to global domain indices (and vise versa) done' 
    101          IF( nn_print >= 1 ) THEN 
    102             WRITE(numout,*) 
    103             WRITE(numout,*) '          conversion local  ==> global i-index domain (mig)' 
    104             WRITE(numout,25)              (mig(ji),ji = 1,jpi) 
    105             WRITE(numout,*) 
    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) 
    111             WRITE(numout,*) 
    112             WRITE(numout,*) '          conversion local  ==> global j-index domain (mjg)' 
    113             WRITE(numout,25)              (mjg(jj),jj = 1,jpj) 
    114             WRITE(numout,*) 
    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) 
    120          ENDIF 
     101 
     102!            WRITE(numout,*) 
     103!            WRITE(numout,*) '          conversion local  ==> global i-index domain (mig)' 
     104!            WRITE(numout,25)              (mig(ji),ji = 1,jpi) 
     105!            WRITE(numout,*) 
     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) 
     111!            WRITE(numout,*) 
     112!            WRITE(numout,*) '          conversion local  ==> global j-index domain (mjg)' 
     113!            WRITE(numout,25)              (mjg(jj),jj = 1,jpj) 
     114!            WRITE(numout,*) 
     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) 
    121120      ENDIF 
    122121 25   FORMAT( 100(10x,19i4,/) ) 
  • utils/tools/DOMAINcfg/src/domclo.F90

    r12442 r13204  
    1212   USE domngb          ! closest point algorithm 
    1313   USE phycst          ! rpi, rad, ra 
    14    USE domutil         ! flood filling algorithm (fill_pool) 
     14   USE domutl          ! flood filling algorithm (fill_pool) 
    1515   USE in_out_manager  ! I/O manager 
    1616   USE lbclnk          ! lateral boundary condition - MPP exchanges 
     
    8282      LOGICAL :: lskip     ! flag in case lake seed on land or already filled (...) 
    8383 
     84      NAMELIST/namclo/ rn_lon_opnsea, rn_lat_opnsea, nn_closea, sn_lake 
     85 
    8486      !!---------------------------------------------------------------------- 
    8587      !! 0 : Read namelist for closed sea definition 
     
    8991      IF(lwp) WRITE(numout,*)'~~~~~~~' 
    9092 
    91       NAMELIST/namclo/ rn_lon_opnsea, rn_lat_opnsea, nn_closea, sn_lake 
    9293      !!--------------------------------------------------------------------- 
    9394       
  • utils/tools/DOMAINcfg/src/domhgr.F90

    r12414 r13204  
    286286         gphi0 = zphi1 + zsin_alpha * ze1deg * REAL( jpjglo-2 , wp ) 
    287287         ! 
    288          IF( nprint==1 .AND. lwp )   THEN 
    289             WRITE(numout,*) '          ze1', ze1, 'cosalpha', zcos_alpha, 'sinalpha', zsin_alpha 
    290             WRITE(numout,*) '          ze1deg', ze1deg, 'glam0', glam0, 'gphi0', gphi0 
    291          ENDIF 
     288    !        WRITE(numout,*) '          ze1', ze1, 'cosalpha', zcos_alpha, 'sinalpha', zsin_alpha 
     289    !        WRITE(numout,*) '          ze1deg', ze1deg, 'glam0', glam0, 'gphi0', gphi0 
    292290         ! 
    293291         DO jj = 1, jpj 
     
    343341      e1_e2v(:,:) = e1v(:,:) / e2v(:,:) 
    344342 
    345       IF( lwp .AND. nn_print >=1 .AND. .NOT.ln_rstart ) THEN      ! Control print : Grid informations (if not restart) 
     343      IF( lwp ) THEN      ! Control print : Grid informations (if not restart) 
    346344         WRITE(numout,*) 
    347345         WRITE(numout,*) '          longitude and e1 scale factors' 
     
    434432      ! ------------------------------------------ 
    435433      ! The equator line must be the latitude coordinate axe 
    436 ! (PM) be carefull with nperio/jperio  
     434 !(PM) be carefull with nperio/jperio  
    437435      IF( jperio == 2 ) THEN 
    438436         znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / REAL( jpi ) 
     
    455453      ! 
    456454      INTEGER ::   inum   ! temporary logical unit 
     455      CHARACTER(LEN=135) :: coordinate_filename 
    457456      !!---------------------------------------------------------------------- 
    458457      ! 
     
    463462      ENDIF 
    464463      ! 
    465       CALL iom_open( 'coordinates', inum ) 
     464      IF (ln_read_cfg) THEN 
     465         coordinate_filename=TRIM(cn_domcfg) 
     466      ELSE 
     467         coordinate_filename='coordinates' 
     468      ENDIF 
     469      CALL iom_open( coordinate_filename, inum ) 
    466470      ! 
    467471      CALL iom_get( inum, jpdom_data, 'glamt', glamt, lrowattr=ln_use_jattr ) 
  • utils/tools/DOMAINcfg/src/domisf.F90

    r12414 r13204  
    1414   !!--------------------------------------------------------------------- 
    1515   USE dom_oce 
    16    USE domutil           ! flood filling algorithm 
     16   USE domutl            ! flood filling algorithm 
    1717   USE domngb            ! find nearest neighbourg 
    1818   USE in_out_manager    ! I/O manager 
  • utils/tools/DOMAINcfg/src/dommsk.F90

    r12414 r13204  
    2626   USE domisf         ! domain: ice shelf 
    2727   USE domwri         ! domain: write the meshmask file 
    28    USE bdy_oce        ! open boundary 
     28   USE usrdef_fmask   ! user defined fmask 
    2929   ! 
    3030   USE in_out_manager ! I/O manager 
     
    5252CONTAINS 
    5353 
    54    SUBROUTINE dom_msk 
     54   SUBROUTINE dom_msk( k_top, k_bot ) 
    5555      !!--------------------------------------------------------------------- 
    5656      !!                 ***  ROUTINE dom_msk  *** 
     
    6262      !!      and ko_bot, the indices of the fist and last ocean t-levels which  
    6363      !!      are either defined in usrdef_zgr or read in zgr_read. 
    64       !!                The velocity masks (umask, vmask)  
     64      !!                The velocity masks (umask, vmask, wmask, wumask, wvmask)  
    6565      !!      are deduced from a product of the two neighboring tmask. 
    6666      !!                The vorticity mask (fmask) is deduced from tmask taking 
     
    7777      !!                due to cyclic or North Fold boundaries as well as MPP halos. 
    7878      !! 
    79       !! ** Action :   tmask, umask, vmask, wmask : land/ocean mask  
     79      !! ** Action :   tmask, umask, vmask, wmask, wumask, wvmask : land/ocean mask  
    8080      !!                         at t-, u-, v- w, wu-, and wv-points (=0. or 1.) 
    8181      !!               fmask   : land/ocean mask at f-point (=0., or =1., or  
     
    8585      !!               ssmask , ssumask, ssvmask, ssfmask : 2D ocean mask 
    8686      !!---------------------------------------------------------------------- 
     87 
     88      INTEGER, DIMENSION(:,:), INTENT(in) ::   k_top, k_bot   ! first and last ocean level 
    8789      ! 
    8890      INTEGER  ::   ji, jj, jk     ! dummy loop indices 
     
    9496      !! 
    9597      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 
    10398      !!--------------------------------------------------------------------- 
    10499      ! 
     
    133128     ! N.B. tmask has already the right boundary conditions since mbathy is ok 
    134129     ! 
     130!      tmask(:,:,:) = 0._wp 
     131!      DO jk = 1, jpk 
     132!         DO jj = 1, jpj 
     133!            DO ji = 1, jpi 
     134!               IF(      ( REAL( mbathy (ji,jj) - jk, wp ) + 0.1_wp >= 0._wp )         & 
     135!               &  .AND. ( REAL( misfdep(ji,jj) - jk, wp ) - 0.1_wp <= 0._wp ) ) THEN 
     136!                  tmask(ji,jj,jk) = 1._wp 
     137!               END IF   
     138!            END DO 
     139!         END DO 
     140!      END DO     
     141  
     142!      IF ( ln_isfsubgl ) CALL zgr_isf_subgl 
     143 
     144      !  Ocean/land mask at t-point  (computed from ko_top and ko_bot) 
     145      ! ---------------------------- 
     146      ! 
    135147      tmask(:,:,:) = 0._wp 
    136       DO jk = 1, jpk 
     148      IF( ln_read_cfg) THEN 
    137149         DO jj = 1, jpj 
    138150            DO ji = 1, jpi 
    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   
     151               iktop = k_top(ji,jj) 
     152               ikbot = k_bot(ji,jj) 
     153               IF( iktop /= 0 ) THEN       ! water in the column 
     154                  tmask(ji,jj,iktop:ikbot  ) = 1._wp 
     155               ENDIF 
     156            END DO   
     157         END DO   
     158         ELSE 
     159         DO jk = 1, jpk 
     160            DO jj = 1, jpj 
     161               DO ji = 1, jpi 
     162                  IF(      ( REAL( mbathy (ji,jj) - jk, wp ) + 0.1_wp >= 0._wp )         & 
     163                  &  .AND. ( REAL( misfdep(ji,jj) - jk, wp ) - 0.1_wp <= 0._wp ) ) THEN 
     164                     tmask(ji,jj,jk) = 1._wp 
     165                  END IF 
     166               END DO 
    143167            END DO 
    144168         END DO 
    145       END DO     
    146   
    147       IF ( ln_isfsubgl ) CALL zgr_isf_subgl 
     169         IF ( ln_isfsubgl ) CALL zgr_isf_subgl 
     170      ENDIF 
     171 
    148172 
    149173!SF  add here lbc_lnk: bug not still understood : cause now domain configuration is read ! 
     
    151175      CALL lbc_lnk( 'dommsk', tmask  , 'T', 1._wp )      ! Lateral boundary conditions 
    152176 
    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) 
    156 903   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 ) 
    159 904   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           
    174177      !  Ocean/land mask at u-, v-, and f-points   (computed from tmask) 
    175178      ! ---------------------------------------- 
     
    272275#if defined key_agrif  
    273276            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  
     277               IF(lk_east) fmask(nlci-1 , :     ,jk) = 0.e0      ! east  
     278               IF(lk_west) fmask(1      , :     ,jk) = 0.e0      ! west  
     279               IF(lk_north) fmask(:      ,nlcj-1 ,jk) = 0.e0      ! north  
     280               IF(lk_south) fmask(:      ,1      ,jk) = 0.e0      ! south  
    278281            ENDIF  
    279282#endif  
     
    287290         ! 
    288291      ENDIF 
    289       ! 
    290       ! write out mesh mask 
    291       IF ( nn_msh > 0 ) CALL dom_wri 
     292       
     293      ! User defined alteration of fmask (use to reduce ocean transport in specified straits) 
     294      ! --------------------------------  
     295      ! 
     296 
     297      CALL usr_def_fmask( cn_cfg, nn_cfg, fmask ) 
    292298      ! 
    293299   END SUBROUTINE dom_msk 
  • utils/tools/DOMAINcfg/src/domngb.F90

    r12414 r13204  
    1111   !!---------------------------------------------------------------------- 
    1212   USE dom_oce        ! ocean space and time domain 
    13    USE phycst 
    1413   ! 
    1514   USE in_out_manager ! I/O manager 
    1615   USE lib_mpp        ! for mppsum 
     16   USE phycst 
    1717 
    1818   IMPLICIT NONE 
     
    4646      INTEGER :: ik         ! working level 
    4747      INTEGER , DIMENSION(2) ::   iloc 
     48      REAL(wp)               ::   zlon, zmini 
    4849      REAL(wp), DIMENSION(jpi,jpj) ::   zglam, zgphi, zmask, zdist 
    4950      !!-------------------------------------------------------------------- 
     
    5960      END SELECT 
    6061 
    61       zdist = dist(plon, plat, zglam, zgphi) 
     62      zlon       = MOD( plon       + 720., 360. )                                     ! plon between    0 and 360 
     63      zglam(:,:) = MOD( zglam(:,:) + 720., 360. )                                     ! glam between    0 and 360 
     64      IF( zlon > 270. )   zlon = zlon - 360.                                          ! zlon between  -90 and 270 
     65      IF( zlon <  90. )   WHERE( zglam(:,:) > 180. ) zglam(:,:) = zglam(:,:) - 360.   ! glam between -180 and 180 
     66      zglam(:,:) = zglam(:,:) - zlon 
    6267 
    6368      IF( lk_mpp ) THEN   
  • utils/tools/DOMAINcfg/src/domutl.F90

    r12414 r13204  
    1 MODULE domutil 
     1MODULE domutl 
    22   !! 
    33   !!====================================================================== 
     
    5454      INTEGER :: iim1, ijm1, ikm1                       ! working integer 
    5555      INTEGER :: nseed                                  ! size of the stack 
    56       TYPE (idx), POINTER :: seed 
     56      TYPE (idx), POINTER :: seed => NULL() 
    5757      !!----------------------------------------------------------------------  
    5858      ! 
    5959      ! initialisation seed 
    60       NULLIFY(seed) 
     60  !    NULLIFY(seed) 
    6161      ! 
    6262      ! create the first seed and update nseed for all processors 
     
    7272         END IF 
    7373      END IF 
    74       nseed=SUM(rseedmap); IF( lk_mpp )   CALL mpp_sum('domutil', nseed )  ! nseed =0 means on land => WARNING later on 
     74      nseed=SUM(rseedmap); IF( lk_mpp )   CALL mpp_sum('domutl', nseed )  ! nseed =0 means on land => WARNING later on 
    7575      ! 
    7676      ! loop until the stack size is 0 or if the pool is larger than the critical size 
     
    106106            ! 
    107107            ! exchange seed 
    108             nseed=SUM(rseedmap); IF( lk_mpp )   CALL mpp_sum('domutil', nseed )  ! this is the sum of all the point check this iteration 
     108            nseed=SUM(rseedmap); IF( lk_mpp )   CALL mpp_sum('domutl', nseed )  ! this is the sum of all the point check this iteration 
    109109            ! 
    110110            rseedmap_b(:,:,:)=rseedmap(:,:,:) 
    111             CALL lbc_lnk('domutil', rseedmap, 'T', 1.) 
     111            CALL lbc_lnk('domutl', rseedmap, 'T', 1.) 
    112112            ! 
    113113            ! build new list of seed 
     
    150150      INTEGER :: iim1, ijm1                   ! working integer 
    151151      INTEGER :: nseed                        ! size of the stack 
    152       TYPE (idx), POINTER :: seed 
     152      TYPE (idx), POINTER :: seed => NULL() 
    153153      !!---------------------------------------------------------------------- 
    154154      ! 
    155155      ! initialisation seed 
    156       NULLIFY(seed) 
     156      !NULLIFY(seed) 
    157157      ! 
    158158      ! create the first seed and update nseed for all processors 
     
    168168         END IF 
    169169      END IF 
    170       nseed=SUM(rseedmap); IF( lk_mpp )   CALL mpp_sum('domutil', nseed )  ! nseed =0 means on land => WARNING later on 
     170      nseed=SUM(rseedmap); IF( lk_mpp )   CALL mpp_sum('domutl', nseed )  ! nseed =0 means on land => WARNING later on 
    171171      ! 
    172172      ! loop until the stack size is 0 or if the pool is larger than the critical size 
     
    195195            ! 
    196196            ! exchange seed 
    197             nseed=SUM(rseedmap); IF( lk_mpp )   CALL mpp_sum('domutil', nseed )  ! this is the sum of all the point check this iteration 
    198             ! 
    199             CALL lbc_lnk('domutil', rseedmap, 'T', 1.) 
     197            nseed=SUM(rseedmap); IF( lk_mpp )   CALL mpp_sum('domutl', nseed )  ! this is the sum of all the point check this iteration 
     198            ! 
     199            CALL lbc_lnk('domutl', rseedmap, 'T', 1.) 
    200200            ! 
    201201            ! build new list of seed 
     
    259259         zpt_tmp => pt_idx%next 
    260260         DEALLOCATE(pt_idx) 
    261          NULLIFY(pt_idx) 
     261         !NULLIFY(pt_idx) 
     262         pt_idx => NULL() 
    262263         pt_idx => zpt_tmp 
    263264      ELSE 
    264          NULLIFY(pt_idx) 
     265         !NULLIFY(pt_idx) 
     266         pt_idx => NULL() 
    265267      ENDIF 
    266268   END SUBROUTINE del_head_idx 
    267269   ! 
    268 END MODULE domutil 
     270END MODULE domutl 
  • utils/tools/DOMAINcfg/src/domzgr.F90

    r12414 r13204  
    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 capability 
     19   !!            3.6  ! 2014-11  (P. Mathiot and C. Harris) add ice shelf capabilitye 
    2020   !!            3.?  ! 2015-11  (H. Liu) Modifications for Wetting/Drying 
    2121   !!---------------------------------------------------------------------- 
     
    3636   !!--------------------------------------------------------------------- 
    3737   USE dom_oce           ! ocean domain 
     38   USE depth_e3       ! depth <=> e3 
    3839   ! 
    3940   USE in_out_manager    ! I/O manager 
     
    4445   USE dombat 
    4546   USE domisf 
     47   USE agrif_domzgr 
    4648 
    4749   IMPLICIT NONE 
     
    9294CONTAINS        
    9395 
    94    SUBROUTINE dom_zgr 
     96   SUBROUTINE dom_zgr( k_top, k_bot ) 
    9597      !!---------------------------------------------------------------------- 
    9698      !!                ***  ROUTINE dom_zgr  *** 
     
    109111      !! ** Action  :   define gdep., e3., mbathy and bathy 
    110112      !!---------------------------------------------------------------------- 
     113      INTEGER, DIMENSION(:,:), INTENT(out) ::   k_top, k_bot   ! ocean first and last level indices 
     114      ! 
    111115      INTEGER ::   ioptio, ibat   ! local integer 
    112116      INTEGER ::   ios 
    113117      ! 
    114       NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav, ln_linssh 
     118      INTEGER :: jk 
     119      REAL(wp) ::   zrefdep             ! depth of the reference level (~10m) 
     120 
     121 
     122      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav 
    115123      !!---------------------------------------------------------------------- 
    116124      ! 
     
    124132902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in configuration namelist', lwp ) 
    125133      IF(lwm) WRITE ( numond, namzgr ) 
     134 
     135      IF(ln_read_cfg) THEN 
     136         IF(lwp) WRITE(numout,*) 
     137         IF(lwp) WRITE(numout,*) '   ==>>>   Read vertical mesh in ', TRIM( cn_domcfg ), ' file' 
     138         ! 
     139         CALL zgr_read   ( ln_zco  , ln_zps  , ln_sco, ln_isfcav,   & 
     140            &              gdept_1d, gdepw_1d, e3t_1d, e3w_1d   ,   &    ! 1D gridpoints depth 
     141            &              gdept_0 , gdepw_0                    ,   &    ! gridpoints depth 
     142            &              e3t_0   , e3u_0   , e3v_0 , e3f_0    ,   &    ! vertical scale factors 
     143            &              e3w_0   , e3uw_0  , e3vw_0           ,   &    ! vertical scale factors 
     144            &              k_top   , k_bot            )                  ! 1st & last ocean level 
     145            ! 
     146!!gm to be remove when removing the OLD definition of e3 scale factors so that gde3w disappears 
     147!      ! Compute gde3w_0 (vertical sum of e3w) 
     148!      gde3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 
     149!      DO jk = 2, jpk 
     150!         gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk) 
     151!      END DO 
     152      ! 
     153 
     154      !                                ! top/bottom ocean level indices for t-, u- and v-points (f-point also for top) 
     155      CALL zgr_top_bot( k_top, k_bot )      ! with a minimum value set to 1 
     156 
     157      !                                ! deepest/shallowest W level Above/Below ~10m 
     158!!gm BUG in s-coordinate this does not work! 
     159      zrefdep = 10._wp - 0.1_wp * MINVAL( e3w_1d )                   ! ref. depth with tolerance (10% of minimum layer thickness) 
     160      nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m 
     161      nla10 = nlb10 - 1                                              ! deepest    W level Above ~10m 
     162!!gm end bug 
     163 
     164      ENDIF         
    126165 
    127166      IF(lwp) THEN                     ! Control print 
     
    134173         WRITE(numout,*) '      s- or hybrid z-s-coordinate    ln_sco    = ', ln_sco 
    135174         WRITE(numout,*) '      ice shelf cavities             ln_isfcav = ', ln_isfcav 
    136          WRITE(numout,*) '      linear free surface            ln_linssh = ', ln_linssh 
    137       ENDIF 
    138  
    139       IF( ln_linssh .AND. lwp) WRITE(numout,*) '   linear free surface: the vertical mesh does not change in time' 
     175      ENDIF 
    140176 
    141177      ioptio = 0                       ! Check Vertical coordinate options 
     
    150186      IF ( ln_sco .AND. ln_isfcav ) ioptio = ioptio + 1 
    151187      IF( ioptio > 0 )   CALL ctl_stop( ' Cavity not tested/compatible with full step (zco) and sigma (ln_sco) ' ) 
     188 
     189      IF(.NOT.ln_read_cfg) THEN 
    152190      ! 
    153191      ! Build the vertical coordinate system 
     
    163201      ! ----------------------------------- 
    164202                          CALL zgr_bot_level    ! deepest ocean level for t-, u- and v-points 
     203                          k_bot = mbkt 
    165204                          CALL zgr_top_level    ! shallowest ocean level for T-, U-, V- points 
    166       ! 
    167       IF( nprint == 1 .AND. lwp )   THEN 
    168          WRITE(numout,*) ' MIN val mbathy  ', MINVAL(  mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 
     205                          k_top = mikt 
     206                          WHERE( bathy(:,:) <= 0._wp )  k_top(:,:) = 0  ! set k_top to zero over land 
     207      ENDIF 
     208      ! 
     209      IF( lwp )   THEN 
     210         WRITE(numout,*) ' MIN val k_top   ', MINVAL(   k_top(:,:) ), ' MAX ', MAXVAL( k_top(:,:) ) 
     211         WRITE(numout,*) ' MIN val k_bot   ', MINVAL(   k_bot(:,:) ), ' MAX ', MAXVAL( k_bot(:,:) ) 
    169212         WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ),   & 
    170213            &                          ' w ', MINVAL( gdepw_0(:,:,:) ) 
     
    181224            &                          ' w ', MAXVAL(   e3w_0(:,:,:) ) 
    182225      ENDIF 
     226 
    183227      ! 
    184228   END SUBROUTINE dom_zgr 
    185229 
     230   SUBROUTINE zgr_read( ld_zco  , ld_zps  , ld_sco  , ld_isfcav,   &   ! type of vertical coordinate 
     231      &                 pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d  ,   &   ! 1D reference vertical coordinate 
     232      &                 pdept , pdepw ,                            &   ! 3D t & w-points depth 
     233      &                 pe3t  , pe3u  , pe3v   , pe3f ,            &   ! vertical scale factors 
     234      &                 pe3w  , pe3uw , pe3vw         ,            &   !     -      -      - 
     235      &                 k_top  , k_bot    )                            ! top & bottom ocean level 
     236      !!--------------------------------------------------------------------- 
     237      !!              ***  ROUTINE zgr_read  *** 
     238      !! 
     239      !! ** Purpose :   Read the vertical information in the domain configuration file 
     240      !! 
     241      !!---------------------------------------------------------------------- 
     242      LOGICAL                   , INTENT(out) ::   ld_zco, ld_zps, ld_sco      ! vertical coordinate flags 
     243      LOGICAL                   , INTENT(out) ::   ld_isfcav                   ! under iceshelf cavity flag 
     244      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pdept_1d, pdepw_1d          ! 1D grid-point depth       [m] 
     245      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pe3t_1d , pe3w_1d           ! 1D vertical scale factors [m] 
     246      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pdept, pdepw                ! grid-point depth          [m] 
     247      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pe3t , pe3u , pe3v , pe3f   ! vertical scale factors    [m] 
     248      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pe3w , pe3uw, pe3vw         !    -       -      - 
     249      INTEGER , DIMENSION(:,:)  , INTENT(out) ::   k_top , k_bot               ! first & last ocean level 
     250      ! 
     251      INTEGER  ::   jk     ! dummy loop index 
     252      INTEGER  ::   inum   ! local logical unit 
     253      REAL(WP) ::   z_zco, z_zps, z_sco, z_cav 
     254      REAL(wp), DIMENSION(jpi,jpj) ::   z2d   ! 2D workspace 
     255      !!---------------------------------------------------------------------- 
     256      ! 
     257      IF(lwp) THEN 
     258         WRITE(numout,*) 
     259         WRITE(numout,*) '   zgr_read : read the vertical coordinates in ', TRIM( cn_domcfg ), ' file' 
     260         WRITE(numout,*) '   ~~~~~~~~' 
     261      ENDIF 
     262      ! 
     263      CALL iom_open( cn_domcfg, inum ) 
     264      ! 
     265      !                          !* type of vertical coordinate 
     266      CALL iom_get( inum, 'ln_zco'   , z_zco ) 
     267      CALL iom_get( inum, 'ln_zps'   , z_zps ) 
     268      CALL iom_get( inum, 'ln_sco'   , z_sco ) 
     269      IF( z_zco == 0._wp ) THEN   ;   ld_zco = .false.   ;   ELSE   ;   ld_zco = .true.   ;   ENDIF 
     270      IF( z_zps == 0._wp ) THEN   ;   ld_zps = .false.   ;   ELSE   ;   ld_zps = .true.   ;   ENDIF 
     271      IF( z_sco == 0._wp ) THEN   ;   ld_sco = .false.   ;   ELSE   ;   ld_sco = .true.   ;   ENDIF 
     272      ! 
     273      !                          !* ocean cavities under iceshelves 
     274      CALL iom_get( inum, 'ln_isfcav', z_cav ) 
     275      IF( z_cav == 0._wp ) THEN   ;   ld_isfcav = .false.   ;   ELSE   ;   ld_isfcav = .true.   ;   ENDIF 
     276      ! 
     277      !                          !* vertical scale factors 
     278      CALL iom_get( inum, jpdom_unknown, 'e3t_1d'  , pe3t_1d  )                     ! 1D reference coordinate 
     279      CALL iom_get( inum, jpdom_unknown, 'e3w_1d'  , pe3w_1d  ) 
     280      ! 
     281      CALL iom_get( inum, jpdom_data, 'e3t_0'  , pe3t  , lrowattr=ln_use_jattr )    ! 3D coordinate 
     282      CALL iom_get( inum, jpdom_data, 'e3u_0'  , pe3u  , lrowattr=ln_use_jattr ) 
     283      CALL iom_get( inum, jpdom_data, 'e3v_0'  , pe3v  , lrowattr=ln_use_jattr ) 
     284      CALL iom_get( inum, jpdom_data, 'e3f_0'  , pe3f  , lrowattr=ln_use_jattr ) 
     285      CALL iom_get( inum, jpdom_data, 'e3w_0'  , pe3w  , lrowattr=ln_use_jattr ) 
     286      CALL iom_get( inum, jpdom_data, 'e3uw_0' , pe3uw , lrowattr=ln_use_jattr ) 
     287      CALL iom_get( inum, jpdom_data, 'e3vw_0' , pe3vw , lrowattr=ln_use_jattr ) 
     288      ! 
     289      !                          !* depths 
     290      !                                   !- old depth definition (obsolescent feature) 
     291      IF(  iom_varid( inum, 'gdept_1d', ldstop = .FALSE. ) > 0  .AND.  & 
     292         & iom_varid( inum, 'gdepw_1d', ldstop = .FALSE. ) > 0  .AND.  & 
     293         & iom_varid( inum, 'gdept_0' , ldstop = .FALSE. ) > 0  .AND.  & 
     294         & iom_varid( inum, 'gdepw_0' , ldstop = .FALSE. ) > 0    ) THEN 
     295         CALL ctl_warn( 'zgr_read : old definition of depths and scale factors used ', &  
     296            &           '           depths at t- and w-points read in the domain configuration file') 
     297         CALL iom_get( inum, jpdom_unknown, 'gdept_1d', pdept_1d )    
     298         CALL iom_get( inum, jpdom_unknown, 'gdepw_1d', pdepw_1d ) 
     299         CALL iom_get( inum, jpdom_data   , 'gdept_0' , pdept , lrowattr=ln_use_jattr ) 
     300         CALL iom_get( inum, jpdom_data   , 'gdepw_0' , pdepw , lrowattr=ln_use_jattr ) 
     301         ! 
     302      ELSE                                !- depths computed from e3. scale factors 
     303         CALL e3_to_depth( pe3t_1d, pe3w_1d, pdept_1d, pdepw_1d )    ! 1D reference depth 
     304         CALL e3_to_depth( pe3t   , pe3w   , pdept   , pdepw    )    ! 3D depths 
     305         IF(lwp) THEN 
     306            WRITE(numout,*) 
     307            WRITE(numout,*) '              Reference 1D z-coordinate depth and scale factors:' 
     308            WRITE(numout, "(9x,' level  gdept_1d  gdepw_1d  e3t_1d   e3w_1d  ')" ) 
     309            WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, pdept_1d(jk), pdepw_1d(jk), pe3t_1d(jk), pe3w_1d(jk), jk = 1, jpk ) 
     310         ENDIF 
     311      ENDIF 
     312      ! 
     313      !                          !* ocean top and bottom level 
     314      CALL iom_get( inum, jpdom_data, 'top_level'    , z2d  , lrowattr=ln_use_jattr )   ! 1st wet T-points (ISF) 
     315      k_top(:,:) = NINT( z2d(:,:) ) 
     316      CALL iom_get( inum, jpdom_data, 'bottom_level' , z2d  , lrowattr=ln_use_jattr )   ! last wet T-points 
     317      k_bot(:,:) = NINT( z2d(:,:) ) 
     318      ! 
     319      ! reference depth for negative bathy (wetting and drying only) 
     320      ! IF( ll_wd )  CALL iom_get( inum,  'rn_wd_ref_depth' , ssh_ref   ) 
     321      ! 
     322      CALL iom_close( inum ) 
     323      ! 
     324   END SUBROUTINE zgr_read 
     325 
     326   SUBROUTINE zgr_top_bot( k_top, k_bot ) 
     327      !!---------------------------------------------------------------------- 
     328      !!                    ***  ROUTINE zgr_top_bot  *** 
     329      !! 
     330      !! ** Purpose :   defines the vertical index of ocean bottom (mbk. arrays) 
     331      !! 
     332      !! ** Method  :   computes from k_top and k_bot with a minimum value of 1 over land 
     333      !! 
     334      !! ** Action  :   mikt, miku, mikv :   vertical indices of the shallowest  
     335      !!                                     ocean level at t-, u- & v-points 
     336      !!                                     (min value = 1) 
     337      !! ** Action  :   mbkt, mbku, mbkv :   vertical indices of the deeptest  
     338      !!                                     ocean level at t-, u- & v-points 
     339      !!                                     (min value = 1 over land) 
     340      !!---------------------------------------------------------------------- 
     341      INTEGER , DIMENSION(:,:), INTENT(in) ::   k_top, k_bot   ! top & bottom ocean level indices 
     342      ! 
     343      INTEGER ::   ji, jj   ! dummy loop indices 
     344      REAL(wp), DIMENSION(jpi,jpj) ::   zk   ! workspace 
     345      !!---------------------------------------------------------------------- 
     346      ! 
     347      IF(lwp) WRITE(numout,*) 
     348      IF(lwp) WRITE(numout,*) '    zgr_top_bot : ocean top and bottom k-index of T-, U-, V- and W-levels ' 
     349      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~' 
     350      ! 
     351      mikt(:,:) = MAX( k_top(:,:) , 1 )    ! top    ocean k-index of T-level (=1 over land) 
     352      ! 
     353      mbkt(:,:) = MAX( k_bot(:,:) , 1 )    ! bottom ocean k-index of T-level (=1 over land) 
     354  
     355      !                                    ! N.B.  top     k-index of W-level = mikt 
     356      !                                    !       bottom  k-index of W-level = mbkt+1 
     357      DO jj = 1, jpjm1 
     358         DO ji = 1, jpim1 
     359            miku(ji,jj) = MAX(  mikt(ji+1,jj  ) , mikt(ji,jj)  ) 
     360            mikv(ji,jj) = MAX(  mikt(ji  ,jj+1) , mikt(ji,jj)  ) 
     361            mikf(ji,jj) = MAX(  mikt(ji  ,jj+1) , mikt(ji,jj), mikt(ji+1,jj  ), mikt(ji+1,jj+1)  ) 
     362            ! 
     363            mbku(ji,jj) = MIN(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  ) 
     364            mbkv(ji,jj) = MIN(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  ) 
     365         END DO 
     366      END DO 
     367      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk  
     368      zk(:,:) = REAL( miku(:,:), wp )   ;   CALL lbc_lnk( 'domzgr', zk, 'U', 1. )   ;   miku(:,:) = MAX( NINT( zk(:,:) ), 1 ) 
     369      zk(:,:) = REAL( mikv(:,:), wp )   ;   CALL lbc_lnk( 'domzgr', zk, 'V', 1. )   ;   mikv(:,:) = MAX( NINT( zk(:,:) ), 1 ) 
     370      zk(:,:) = REAL( mikf(:,:), wp )   ;   CALL lbc_lnk( 'domzgr', zk, 'F', 1. )   ;   mikf(:,:) = MAX( NINT( zk(:,:) ), 1 ) 
     371      ! 
     372      zk(:,:) = REAL( mbku(:,:), wp )   ;   CALL lbc_lnk( 'domzgr', zk, 'U', 1. )   ;   mbku(:,:) = MAX( NINT( zk(:,:) ), 1 ) 
     373      zk(:,:) = REAL( mbkv(:,:), wp )   ;   CALL lbc_lnk( 'domzgr', zk, 'V', 1. )   ;   mbkv(:,:) = MAX( NINT( zk(:,:) ), 1 ) 
     374      ! 
     375   END SUBROUTINE zgr_top_bot 
    186376 
    187377   SUBROUTINE zgr_z 
     
    665855      ENDIF 
    666856 
    667       zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    668       CALL lbc_lnk( 'domzgr',zbathy, 'T', 1._wp ) 
    669       mbathy(:,:) = INT( zbathy(:,:) ) 
    670  
     857      IF( lk_mpp ) THEN 
     858         zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     859         CALL lbc_lnk( 'toto',zbathy, 'T', 1._wp ) 
     860         mbathy(:,:) = INT( zbathy(:,:) ) 
     861      ENDIF 
    671862      !                                          ! East-west cyclic boundary conditions 
    672863      IF( jperio == 0 ) THEN 
     
    10611252   END SUBROUTINE zgr_zps 
    10621253 
     1254 
    10631255   SUBROUTINE zgr_sco 
    10641256      !!---------------------------------------------------------------------- 
     
    12621454         END DO 
    12631455         ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero 
    1264          CALL lbc_lnk( 'domzgr',zenv, 'T', 1._wp, 'no0' ) 
     1456         CALL lbc_lnk( 'toto',zenv, 'T', 1._wp, 'no0' ) 
    12651457         !                                                  ! ================ ! 
    12661458      END DO                                                !     End loop     ! 
     
    13411533        hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) ) 
    13421534 
    1343       IF( nprint == 1 .AND. lwp )   THEN 
     1535      IF( lwp )   THEN 
    13441536         WRITE(numout,*) ' MAX val hif   t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ),  & 
    13451537            &                        ' u ',   MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) ) 
     
    13971589         END DO 
    13981590      END DO 
    1399       IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ),   & 
     1591      IF(lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ),   & 
    14001592         &                                                       ' MAX ', MAXVAL( mbathy(:,:) ) 
    14011593 
    1402       IF( nprint == 1  .AND. lwp )   THEN         ! min max values over the local domain 
     1594      IF( lwp )   THEN         ! min max values over the local domain 
    14031595         WRITE(numout,*) ' MIN val mbathy  ', MINVAL( mbathy(:,:)    ), ' MAX ', MAXVAL( mbathy(:,:) ) 
    14041596         WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ),   & 
     
    17881980        z_gsigt(jk) = -fssig( REAL(jk,wp)        ) 
    17891981      END DO 
    1790       IF( nprint == 1 .AND. lwp )   WRITE(numout,*) 'z_gsigw 1 jpk    ', z_gsigw(1), z_gsigw(jpk) 
     1982      IF( lwp )   WRITE(numout,*) 'z_gsigw 1 jpk    ', z_gsigw(1), z_gsigw(jpk) 
    17911983      ! 
    17921984      ! Coefficients for vertical scale factors at w-, t- levels 
  • utils/tools/DOMAINcfg/src/in_out_manager.F90

    r12414 r13204  
    2222   !!---------------------------------------------------------------------- 
    2323   CHARACTER(lc) ::   cn_exp           !: experiment name used for output filename 
    24    CHARACTER(lc) ::   cn_ocerst_in     !: suffix of ocean restart name (input) 
    25    CHARACTER(lc) ::   cn_ocerst_indir  !: restart input directory 
    26    CHARACTER(lc) ::   cn_ocerst_out    !: suffix of ocean restart name (output) 
    27    CHARACTER(lc) ::   cn_ocerst_outdir !: restart output directory 
    28    LOGICAL       ::   ln_rstart        !: start from (F) rest or (T) a restart file 
    29    LOGICAL       ::   ln_rst_list      !: output restarts at list of times (T) or by frequency (F) 
    30    INTEGER       ::   nn_rstctl        !: control of the time step (0, 1 or 2) 
    31    INTEGER       ::   nn_rstssh   = 0  !: hand made initilization of ssh or not (1/0) 
    3224   INTEGER       ::   nn_it000         !: index of the first time step 
    3325   INTEGER       ::   nn_itend         !: index of the last time step 
     
    3527   INTEGER       ::   nn_time0         !: initial time of day in hhmm 
    3628   INTEGER       ::   nn_leapy         !: Leap year calendar flag (0/1 or 30) 
    37    INTEGER       ::   nn_istate        !: initial state output flag (0/1) 
    38    INTEGER       ::   nn_write         !: model standard output frequency 
    39    INTEGER       ::   nn_stock         !: restart file frequency 
    40    INTEGER, DIMENSION(10) :: nn_stocklist  !: restart dump times 
    4129   LOGICAL       ::   ln_mskland       !: mask land points in NetCDF outputs (costly: + ~15%) 
    4230   LOGICAL       ::   ln_cfmeta        !: output additional data to netCDF files required for compliance with the CF metadata standard 
     
    4533   LOGICAL       ::   ln_xios_read     !: use xios to read single file restart 
    4634   INTEGER       ::   nn_wxios         !: write resart using xios 0 - no, 1 - single, 2 - multiple file output 
    47    INTEGER       ::   nn_no            !: Assimilation cycle 
    4835 
    4936#if defined key_netcdf4 
     
    7461 
    7562   CHARACTER(lc) ::   cexper                      !: experiment name used for output filename 
    76    INTEGER       ::   nrstdt                      !: control of the time step (0, 1 or 2) 
    7763   INTEGER       ::   nit000                      !: index of the first time step 
    7864   INTEGER       ::   nitend                      !: index of the last time step 
    7965   INTEGER       ::   ndate0                      !: initial calendar date aammjj 
    8066   INTEGER       ::   nleapy                      !: Leap year calendar flag (0/1 or 30) 
    81    INTEGER       ::   ninist                      !: initial state output flag (0/1) 
    82    INTEGER       ::   nwrite                      !: model standard output frequency 
    83    INTEGER       ::   nstock                      !: restart file frequency 
    84    INTEGER, DIMENSION(10) :: nstocklist           !: restart dump times 
    85  
    86    !!---------------------------------------------------------------------- 
    87    !! was in restart but moved here because of the OFF line... better solution should be found... 
    88    !!---------------------------------------------------------------------- 
    89    INTEGER ::   nitrst                !: time step at which restart file should be written 
    90    LOGICAL ::   lrst_oce              !: logical to control the oce restart write  
    91    LOGICAL ::   lrst_ice              !: logical to control the ice restart write  
    92    INTEGER ::   numror = 0            !: logical unit for ocean restart (read). Init to 0 is needed for SAS (in daymod.F90) 
    93    INTEGER ::   numrir                !: logical unit for ice   restart (read) 
    94    INTEGER ::   numrow                !: logical unit for ocean restart (write) 
    95    INTEGER ::   numriw                !: logical unit for ice   restart (write) 
    96    INTEGER ::   nrst_lst              !: number of restart to output next 
    97  
    98    !!---------------------------------------------------------------------- 
    99    !!                    output monitoring 
    100    !!---------------------------------------------------------------------- 
    101    LOGICAL ::   ln_ctl           !: run control for debugging 
    102    TYPE :: sn_ctl                !: optional use structure for finer control over output selection 
    103       LOGICAL :: l_config  = .FALSE.  !: activate/deactivate finer control 
    104                                       !  Note if l_config is True then ln_ctl is ignored. 
    105                                       !  Otherwise setting ln_ctl True is equivalent to setting 
    106                                       !  all the following logicals in this structure True 
    107       LOGICAL :: l_runstat = .FALSE.  !: Produce/do not produce run.stat file (T/F) 
    108       LOGICAL :: l_trcstat = .FALSE.  !: Produce/do not produce tracer.stat file (T/F) 
    109       LOGICAL :: l_oceout  = .FALSE.  !: Produce all ocean.outputs    (T) or just one (F) 
    110       LOGICAL :: l_layout  = .FALSE.  !: Produce all layout.dat files (T) or just one (F) 
    111       LOGICAL :: l_mppout  = .FALSE.  !: Produce/do not produce mpp.output_XXXX files (T/F) 
    112       LOGICAL :: l_mpptop  = .FALSE.  !: Produce/do not produce mpp.top.output_XXXX files (T/F) 
    113                                       !  Optional subsetting of processor report files 
    114                                       !  Default settings of 0/1000000/1 should ensure all areas report. 
    115                                       !  Set to a more restrictive range to select specific areas 
    116       INTEGER :: procmin   = 0        !: Minimum narea to output 
    117       INTEGER :: procmax   = 1000000  !: Maximum narea to output 
    118       INTEGER :: procincr  = 1        !: narea increment to output 
    119       INTEGER :: ptimincr  = 1        !: timestep increment to output (time.step and run.stat) 
    120    END TYPE sn_ctl 
    121  
    122    TYPE (sn_ctl) :: sn_cfctl     !: run control structure for selective output 
    123    LOGICAL ::   ln_timing        !: run control for timing 
    124    LOGICAL ::   ln_diacfl        !: flag whether to create CFL diagnostics 
    125    INTEGER ::   nn_print         !: level of print (0 no print) 
    126    INTEGER ::   nn_ictls         !: Start i indice for the SUM control 
    127    INTEGER ::   nn_ictle         !: End   i indice for the SUM control 
    128    INTEGER ::   nn_jctls         !: Start j indice for the SUM control 
    129    INTEGER ::   nn_jctle         !: End   j indice for the SUM control