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

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

Changeset 13024


Ignore:
Timestamp:
2020-06-03T16:26:23+02:00 (4 years ago)
Author:
rblod
Message:

First version of new nesting tools merged with domaincfg, see ticket #2129

Location:
utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg
Files:
28 added
25 edited

Legend:

Unmodified
Added
Removed

  • utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/namelist_cfg

    r12414 r13024  
    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.              ! 
     44/ 
     45!----------------------------------------------------------------------- 
     46&namcrs        !   Grid coarsening for dynamics output and/or 
     47               !   passive tracer coarsened online simulations 
     48!----------------------------------------------------------------------- 
     49/ 
     50!----------------------------------------------------------------------- 
     51&namtsd    !   data : Temperature  & Salinity 
     52!----------------------------------------------------------------------- 
     53/ 
     54!----------------------------------------------------------------------- 
     55&namsbc        !   Surface Boundary Condition (surface module) 
     56!----------------------------------------------------------------------- 
     57/ 
     58!----------------------------------------------------------------------- 
     59&namsbc_core   !   namsbc_core  CORE bulk formulae 
     60!----------------------------------------------------------------------- 
     61/ 
     62!----------------------------------------------------------------------- 
     63&namtra_qsr    !   penetrative solar radiation 
     64!----------------------------------------------------------------------- 
     65/ 
     66!----------------------------------------------------------------------- 
     67&namsbc_rnf    !   runoffs namelist surface boundary condition 
     68!----------------------------------------------------------------------- 
     69/ 
     70!----------------------------------------------------------------------- 
     71&namsbc_ssr    !   surface boundary condition : sea surface restoring 
     72!----------------------------------------------------------------------- 
     73/ 
     74!----------------------------------------------------------------------- 
     75&namsbc_alb    !   albedo parameters 
     76!----------------------------------------------------------------------- 
     77/ 
     78!----------------------------------------------------------------------- 
     79&namberg       !   iceberg parameters 
     80!----------------------------------------------------------------------- 
     81/ 
     82!----------------------------------------------------------------------- 
     83&namlbc        !   lateral momentum boundary condition 
     84!----------------------------------------------------------------------- 
     85/ 
     86!----------------------------------------------------------------------- 
     87&nambfr        !   bottom friction 
     88!----------------------------------------------------------------------- 
     89/ 
     90!----------------------------------------------------------------------- 
     91&nambbc        !   bottom temperature boundary condition                (default: NO) 
     92!----------------------------------------------------------------------- 
     93   ln_trabbc   = .true.    !  Apply a geothermal heating at the ocean bottom 
     94/ 
     95!----------------------------------------------------------------------- 
     96&nambbl        !   bottom boundary layer scheme 
     97!----------------------------------------------------------------------- 
    1798/ 
    1899!----------------------------------------------------------------------- 
    19100&namcfg        !   parameters of the configuration 
    20101!----------------------------------------------------------------------- 
     102   ! 
     103   ln_e3_dep   = .true.    ! =T : e3=dk[depth] in discret sens. 
     104   !                       !      ===>>> will become the only possibility in v4.0 
     105   !                       ! =F : e3 analytical derivative of depth function 
     106   !                       !      only there for backward compatibility test with v3.6 
     107   !                       ! 
     108   cp_cfg      =  "orca"   !  name of the configuration 
     109   jp_cfg      =       2   !  resolution of the configuration 
     110   jpidta      =     182   !  1st lateral dimension ( >= jpi ) 
     111   jpjdta      =     149   !  2nd    "         "    ( >= jpj ) 
     112   jpkdta      =      31   !  number of levels      ( >= jpk ) 
     113   jpiglo      =     182   !  1st dimension of global domain --> i =jpidta 
     114   jpjglo      =     149   !  2nd    -                  -    --> j  =jpjdta 
     115   jperio      =       4   !  lateral cond. type (between 0 and 6) 
     116   ln_use_jattr = .false.  !  use (T) the file attribute: open_ocean_jstart, if present 
     117                           !  in netcdf input files, as the start j-row for reading 
     118   ln_domclo = .false.     ! computation of closed sea masks (see namclo) 
    21119/ 
    22120!----------------------------------------------------------------------- 
    23121&namzgr        !   vertical coordinate                                  (default: NO selection) 
    24122!----------------------------------------------------------------------- 
     123!----------------------------------------------------------------------- 
     124   ln_zco      = .false.   !  z-coordinate - full    steps 
     125   ln_zps      = .true.   !  z-coordinate - partial steps 
     126   ln_sco      = .false.   !  s- or hybrid z-s-coordinate 
     127   ln_isfcav   = .false.   !  ice shelf cavity             (T: see namzgr_isf) 
     128   ln_linssh   = .true.   !  linear free surface 
    25129/ 
    26130!----------------------------------------------------------------------- 
    … …  
    39143&namlbc        !   lateral momentum boundary condition                  (default: NO selection) 
    40144!----------------------------------------------------------------------- 
     145jpni = 0 
     146jpnj=0 
    41147/ 
    42148!----------------------------------------------------------------------- 
    43149&namagrif      !  AGRIF zoom                                            ("key_agrif") 
     150ln_bry_south = .TRUE. 
    44151!----------------------------------------------------------------------- 
    45152/ 

  • utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/namelist_ref

    r12414 r13024  
    4040&namdom        !   space and time domain (bathymetry, mesh, timestep) 
    4141!----------------------------------------------------------------------- 
    42    nn_bathy    =    1      !  compute (=0) or read (=1) the bathymetry file 
     42   nn_bathy    =    1      !  compute analyticaly (=0) or read (=1) the bathymetry file 
     43                           !  or compute (2) from external bathymetry 
     44   nn_interp   =    1                          ! type of interpolation (nn_bathy =2) 
     45   cn_topo     =  'bathymetry_ORCA12_V3.3.nc'  ! external topo file (nn_bathy =2) 
     46   cn_bath     =  'Bathymetry'                 ! topo name in file  (nn_bathy =2) 
     47   cn_lon      =  'nav_lon'                    ! lon  name in file  (nn_bathy =2) 
     48   cn_lat      =  'nav_lat'                    ! lat  name in file  (nn_bathy =2) 
    4349   rn_bathy    =    0.     !  value of the bathymetry. if (=0) bottom flat at jpkm1 
    4450   nn_msh      =    0      !  create (=1) a mesh file or not (=0) 
    … …  
    7379   ppkth2      =       48.029893720000 ! 
    7480   ppacr2      =       13.000000000000 ! 
    75 / 
     81// 
    7682!----------------------------------------------------------------------- 
    7783&namcfg        !   parameters of the configuration 
    … …  
    8288   !                       ! =F : e3 analytical derivative of depth function 
    8389   !                       !      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 
     90   !                       ! 
     91   cp_cfg      =  "orca"   !  name of the configuration 
     92   jp_cfg      =       2   !  resolution of the configuration 
     93   jpidta      =     182   !  1st lateral dimension ( >= jpi ) 
     94   jpjdta      =     149   !  2nd    "         "    ( >= jpj ) 
     95   jpkdta      =      31   !  number of levels      ( >= jpk ) 
     96   jpiglo      =     182   !  1st dimension of global domain --> i =jpidta 
     97   jpjglo      =     149   !  2nd    -                  -    --> j  =jpjdta 
     98   jperio      =       4   !  lateral cond. type (between 0 and 6) 
    9799   ln_use_jattr = .false.  !  use (T) the file attribute: open_ocean_jstart, if present 
    98100                           !  in netcdf input files, as the start j-row for reading 
    … …  
    182184   rn_sponge_tra = 2880.   !  coefficient for tracer   sponge layer [m2/s] 
    183185   rn_sponge_dyn = 2880.   !  coefficient for dynamics sponge layer [m2/s] 
    184    ln_chk_bathy  = .false. !  =T  check the parent bathymetry 
     186   ln_chk_bathy  = .FALSE. ! 
     187   npt_connect   = 2 
     188   npt_copy      = 2 
     189   ln_bry_south = .TRUE. 
     190/ 
     191!----------------------------------------------------------------------- 
     192&nam_tide      !   tide parameters                                      ("key_tide") 
     193!----------------------------------------------------------------------- 
     194   ln_tide_pot = .true.    !  use tidal potential forcing 
     195   ln_tide_ramp= .false.   ! 
     196   rdttideramp =    0.     ! 
     197   clname(1)   = 'DUMMY'   !  name of constituent - all tidal components must be set in namelist_cfg 
    185198/ 
    186199!----------------------------------------------------------------------- 
    … …  
    210223   rn_ice_sal    = 10.        !  si3 only:      --   salinity           -- 
    211224   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 
    232225/ 
    233226!----------------------------------------------------------------------- 
    … …  
    238231   nn_buffer   =   0       !  size in bytes of exported buffer ('B' case), 0 no exportation 
    239232   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) 
     233   jpni        =   1       !  jpni   number of processors following i (set automatically if < 1) 
     234   jpnj        =   1       !  jpnj   number of processors following j (set automatically if < 1) 
    242235/ 
    243236!----------------------------------------------------------------------- 

  • utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/agrif_boundary_connections.F90

    r12414 r13024  
    1313call Agrif_Bc_variable(e3t_copy_id, procname = connect_e3t_copy) 
    1414 
    15 ! Agrif_UseSpecialValue = .TRUE. 
    16 ! Agrif_SpecialValue = 0. 
    17 ! call Agrif_Bc_variable(e3t_connect_id, procname = connect_e3t_connect) 
    18 ! Agrif_UseSpecialValue = .FALSE. 
    19  
     15Allocate(e3t_interp(jpi,jpj,jpk)) 
     16e3t_interp = -10. 
     17Agrif_UseSpecialValue = .TRUE. 
     18Agrif_SpecialValue = 0. 
     19call Agrif_Bc_variable(e3t_connect_id, procname = connect_e3t_connect) 
     20Agrif_UseSpecialValue = .FALSE. 
     21 
     22! Call Agrif_make_connection() 
    2023 
    2124      Agrif_SpecialValue    = 0. 
    2225      Agrif_UseSpecialValue = ln_spc_dyn 
    2326      ! 
    24       CALL Agrif_Bc_variable( e3u_id, procname=connect_e3u ) 
    25       CALL Agrif_Bc_variable( e3v_id, procname=connect_e3v ) 
     27!      CALL Agrif_Bc_variable( e3u_id, procname=connect_e3u ) 
     28!      CALL Agrif_Bc_variable( e3v_id, procname=connect_e3v ) 
    2629      ! 
    2730      Agrif_UseSpecialValue = .FALSE. 
    … …  
    8386      ELSE 
    8487         mbkt(i1:i2,j1:j2) = nint(ptab(i1:i2,j1:j2)) 
    85           
    8688         WHERE (mbkt(i1:i2,j1:j2)==0) 
    8789           ssmask(i1:i2,j1:j2) = 0. 
     90         ELSEWHERE 
     91           ssmask(i1:i2,j1:j2) = 1. 
    8892         END WHERE 
    8993            
    … …  
    136140         do ji=i1,i2 
    137141           bathy_local (ji,jj) = SUM ( e3t_0(ji,jj, 1:mbkt(ji,jj) ) ) * ssmask(ji,jj) 
    138            print *,'ji = ',ji,jj,bathy_local(ji,jj),ptab(ji,jj,jpk+1) 
    139          enddo 
    140          enddo 
    141           
    142          ! DO jk=1,jpk 
    143            ! DO jj=j1,j2 
    144              ! DO ji=i1,i2 
    145                  ! e3t_0(ji,jj,jk) = MAX(ptab(ji,jj,jk),MIN(e3zps_min, e3t_1d(jk)*e3zps_rat )) 
    146                  ! e3t_0(ji,jj,jk) = MIN(e3t_0(ji,jj,jk),e3t_1d(jk)) 
    147              ! ENDDO 
    148            ! ENDDO 
    149          ! ENDDO 
     142         enddo 
     143         enddo 
     144          
     145         DO jk=1,jpk 
     146           DO jj=j1,j2 
     147             DO ji=i1,i2 
     148             if (e3t_interp(ji,jj,jk) == -10) then ! the connection has not yet been done 
     149                 e3t_interp(ji,jj,jk) = MAX(ptab(ji,jj,jk),MIN(e3zps_min, e3t_1d(jk)*e3zps_rat )) 
     150                 e3t_interp(ji,jj,jk) = MIN(e3t_interp(ji,jj,jk),e3t_1d(jk)) 
     151                 e3t_0(ji,jj,jk) = ztabramp(ji,jj)*e3t_0(ji,jj,jk)+(1.-ztabramp(ji,jj))*e3t_interp(ji,jj,jk) 
     152             endif 
     153             ENDDO 
     154           ENDDO 
     155         ENDDO 
    150156      ENDIF 
    151157      ! 

  • utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/agrif_connection.F90

    r12414 r13024  
    77      !!---------------------------------------------------------------------- 
    88      INTEGER  ::   ji, jj, ind1, ind2 
    9       INTEGER  ::   ispongearea 
     9      INTEGER  ::   ispongearea, istart 
    1010      REAL(wp) ::   z1_spongearea 
    1111      !!---------------------------------------------------------------------- 
    … …  
    1515 
    1616         ALLOCATE(ztabramp(jpi,jpj)) 
    17          ispongearea  = 1 + npt_connect * Agrif_irhox() 
     17         ispongearea = 1 + npt_connect * Agrif_irhox() 
     18         istart = npt_copy * Agrif_irhox() 
    1819         z1_spongearea = 1._wp / REAL( ispongearea ) 
    1920          
    … …  
    2122 
    2223         ! --- West --- ! 
    23          IF( (nbondi == -1) .OR. (nbondi == 2) ) THEN 
    24             ind1 = 1+nbghostcells 
    25             ind2 = 1+nbghostcells + ispongearea  
     24         IF( ((nbondi == -1) .OR. (nbondi == 2) ).AND. .NOT. (jperio == 1 .OR. jperio == 4 .OR. jperio == 6)) THEN 
     25            ind1 = 1+nbghostcells + istart 
     26            ind2 = ind1 + ispongearea  
    2627            DO jj = 1, jpj 
    2728               DO ji = ind1, ind2                 
    … …  
    3233 
    3334         ! --- East --- ! 
    34          IF( (nbondi == 1) .OR. (nbondi == 2) ) THEN 
    35             ind1 = nlci - nbghostcells - ispongearea 
    36             ind2 = nlci - nbghostcells 
     35         IF( ((nbondi == 1) .OR. (nbondi == 2) ).AND. .NOT. (jperio == 1 .OR. jperio == 4 .OR. jperio == 6)) THEN 
     36            !ind1 = nlci - nbghostcells - ispongearea 
     37            ind2 = nlci - nbghostcells - istart 
     38            ind1 = ind2 -ispongearea 
     39            
     40             
    3741            DO jj = 1, jpj 
    3842               DO ji = ind1, ind2 
    … …  
    4347 
    4448         ! --- South --- ! 
    45          IF( (nbondj == -1) .OR. (nbondj == 2) ) THEN 
    46             ind1 = 1+nbghostcells 
    47             ind2 = 1+nbghostcells + ispongearea 
     49         IF(( (nbondj == -1) .OR. (nbondj == 2) ).AND.(ln_bry_south)) THEN 
     50            ! ind1 = 1+nbghostcells 
     51            ! ind2 = 1+nbghostcells + ispongearea 
     52            ind1 = 1+nbghostcells + istart 
     53            ind2 = ind1 + ispongearea  
    4854            DO jj = ind1, ind2  
    4955               DO ji = 1, jpi 
    … …  
    5561         ! --- North --- ! 
    5662         IF( (nbondj == 1) .OR. (nbondj == 2) ) THEN 
    57             ind1 = nlcj - nbghostcells - ispongearea 
    58             ind2 = nlcj - nbghostcells 
     63            ! ind1 = nlcj - nbghostcells - ispongearea 
     64            ! ind2 = nlcj - nbghostcells 
     65             
     66            ind2 = nlcj - nbghostcells - istart 
     67            ind1 = ind2 -ispongearea 
     68             
    5969            DO jj = ind1, ind2 
    6070               DO ji = 1, jpi 
    … …  
    6777   END SUBROUTINE Agrif_connection 
    6878    
     79   SUBROUTINE Agrif_make_connection 
     80   use dom_oce 
     81   use agrif_parameters 
     82      !!---------------------------------------------------------------------- 
     83      !!                 *** ROUTINE  Agrif_Sponge *** 
     84      !!---------------------------------------------------------------------- 
     85      INTEGER  ::   ji, jj, ind1, ind2 
     86      INTEGER  ::   ispongearea, istart 
     87      REAL(wp) ::   z1_spongearea 
     88      !!---------------------------------------------------------------------- 
     89      ! 
     90         ! Define ramp from boundaries towards domain interior at T-points 
     91         ! Store it in ztabramp 
     92 
     93         ispongearea = 1 + npt_connect * Agrif_irhox() 
     94         istart = npt_copy * Agrif_irhox() 
     95 
     96         ! --- West --- ! 
     97         IF( (nbondi == -1) .OR. (nbondi == 2) ) THEN 
     98             
     99            ind1 = 1+nbghostcells + istart 
     100            ind2 = ind1 + ispongearea 
     101            DO jk=1,jpk             
     102             DO jj = 1, jpj 
     103               DO ji = ind1, ind2 
     104                 
     105                ! print *,'VAL = ',ztabramp(ji,jj)*e3t_interp(ji,jj,jk)+(1.-ztabramp(ji,jj))*e3t_0(ji,jj,jk), & 
     106                ! e3t_0(ji,jj,jk) 
     107                  e3t_0(ji,jj,jk) = ztabramp(ji,jj)*e3t_interp(ji,jj,jk)+(1.-ztabramp(ji,jj))*e3t_0(ji,jj,jk) 
     108               ENDDO 
     109              ENDDO 
     110            ENDDO 
     111         ENDIF 
     112 
     113         ! --- East --- ! 
     114         IF( (nbondi == 1) .OR. (nbondi == 2) ) THEN 
     115            ind2 = nlci - nbghostcells - istart 
     116            ind1 = ind2 -ispongearea 
     117            DO jk=1,jpk     
     118            DO jj = 1, jpj 
     119               DO ji = ind1, ind2 
     120                  e3t_0(ji,jj,jk) = ztabramp(ji,jj)*e3t_interp(ji,jj,jk)+(1.-ztabramp(ji,jj))*e3t_0(ji,jj,jk) 
     121               ENDDO 
     122            ENDDO 
     123            ENDDO 
     124         ENDIF 
     125 
     126         ! --- South --- ! 
     127         IF( (nbondj == -1) .OR. (nbondj == 2) ) THEN 
     128            ind1 = 1+nbghostcells + istart 
     129            ind2 = ind1 + ispongearea  
     130            DO jk=1,jpk     
     131            DO jj = ind1, ind2  
     132               DO ji = 1, jpi 
     133                  e3t_0(ji,jj,jk) = ztabramp(ji,jj)*e3t_interp(ji,jj,jk)+(1.-ztabramp(ji,jj))*e3t_0(ji,jj,jk) 
     134               END DO 
     135            ENDDO 
     136            ENDDO 
     137         ENDIF 
     138 
     139         ! --- North --- ! 
     140         IF( (nbondj == 1) .OR. (nbondj == 2) ) THEN 
     141             
     142            ind2 = nlcj - nbghostcells - istart 
     143            ind1 = ind2 -ispongearea 
     144            DO jk=1,jpk                
     145            DO jj = ind1, ind2 
     146               DO ji = 1, jpi 
     147                  e3t_0(ji,jj,jk) = ztabramp(ji,jj)*e3t_interp(ji,jj,jk)+(1.-ztabramp(ji,jj))*e3t_0(ji,jj,jk) 
     148               END DO 
     149            ENDDO 
     150            ENDDO 
     151         ENDIF 
     152      ! 
     153      ! 
     154   END SUBROUTINE Agrif_make_connection 
     155    
    69156#else 
    70157subroutine agrif_connection_empty 

  • utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/agrif_domzgr.F90

    r12414 r13024  
    2323      ! 
    2424      !!----------------------------------------------------------------------   
     25      INTEGER :: ji,jj 
    2526      ! 
    2627         western_side  = (nb == 1).AND.(ndir == 1) 
    … …  
    3031      IF( before) THEN 
    3132         ptab(i1:i2,j1:j2) = bathy(i1:i2,j1:j2) 
     33 
     34         do jj=j1,j2 
     35            do ji=i1,i2 
     36         ptab(ji,jj) = SUM ( e3t_0(ji,jj, 1:mbkt(ji,jj) ) ) * ssmask(ji,jj) 
     37      enddo 
     38      enddo 
    3239      ELSE 
    3340         bathy(i1:i2,j1:j2)=ptab 

  • utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/agrif_parameters.F90

    r12414 r13024  
    1010INTEGER :: npt_copy 
    1111INTEGER :: npt_connect 
    12  
     12LOGICAL :: ln_bry_south = .TRUE. 
    1313REAL(wp), PUBLIC, ALLOCATABLE, SAVE        , DIMENSION(:,:) ::   ztabramp 
     14REAL(wp), PUBLIC, ALLOCATABLE, SAVE        , DIMENSION(:,:,:) ::   e3t_interp 
    1415 
    1516end module agrif_parameters 

  • utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/agrif_recompute_scalefactors.f90

    r12414 r13024  
    99      ! Scale factors and depth at U-, V-, UW and VW-points 
    1010      DO jk = 1, jpk                        ! initialisation to z-scale factors 
    11 !         e3u_0 (:,:,jk) = e3t_1d(jk) 
    12 !         e3v_0 (:,:,jk) = e3t_1d(jk) 
     11         e3u_0 (:,:,jk) = e3t_1d(jk) 
     12         e3v_0 (:,:,jk) = e3t_1d(jk) 
    1313         e3uw_0(:,:,jk) = e3w_1d(jk) 
    1414         e3vw_0(:,:,jk) = e3w_1d(jk) 
    … …  
    1818         DO jj = 1, jpjm1 
    1919            DO ji = 1, jpim1   ! vector opt. 
    20 !               e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) ) 
    21 !               e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) ) 
     20               e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) ) 
     21               e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) ) 
    2222               e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) ) 
    2323               e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) ) 
    … …  
    2727      IF ( ln_isfcav ) THEN 
    2828      ! (ISF) define e3uw (adapted for 2 cells in the water column) 
     29      print *,'NOT READY SINCE:' 
     30      print *,'MBATHY HAS NOT BEEN CORRECTED / UPDATED' 
     31      print *,'EVEN NOT COMPUTED IN THE CASE ln_read_cfg = .TRUE.' 
     32      STOP 
    2933         DO jj = 2, jpjm1  
    3034            DO ji = 2, jpim1   ! vector opt.  

  • utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/agrif_update.F90

    r12414 r13024  
    66 
    77if (agrif_root()) return 
     8 
    89call agrif_update_variable(bottom_level_id,locupdate=(/npt_copy,0/),procname = update_bottom_level) 
     10 
    911 
    1012      Agrif_UseSpecialValueInUpdate = .TRUE. 
    … …  
    1416      Agrif_UseSpecialValueInUpdate = .FALSE. 
    1517 
    16 call agrif_update_variable(e3u_id,procname = update_e3u) 
    17 call agrif_update_variable(e3v_id,procname = update_e3v) 
     18!call agrif_update_variable(e3u_id,procname = update_e3u) 
     19!call agrif_update_variable(e3v_id,procname = update_e3v) 
    1820       
    1921end subroutine agrif_update_all 
    … …  
    4648         WHERE (mbkt(i1:i2,j1:j2)==0) 
    4749           ssmask(i1:i2,j1:j2) = 0. 
     50         ELSEWHERE 
     51           ssmask(i1:i2,j1:j2) = 1. 
    4852         END WHERE 
    4953            
    … …  
    6973               DO jj=j1,j2 
    7074                  DO ji=i1,i2 
    71                    if (mbkt(ji,jj) < jk) then 
     75                   if (mbkt(ji,jj) <= jk) then 
    7276                     tabres(ji,jj,jk) = e3t_0(ji,jj,jk) 
    7377                   else 
    … …  
    8185               DO jj=j1,j2 
    8286                  DO ji=i1,i2 
    83                    if (mbkt(ji,jj) < jk) then 
     87                   if (mbkt(ji,jj) <= jk) then 
     88                     e3t_0(ji,jj,jk) = MAX(tabres(ji,jj,jk),MIN(e3zps_min,e3t_1d(jk)*e3zps_rat)) 
     89                   else 
    8490                     e3t_0(ji,jj,jk) = e3t_1d(jk) 
    85                    else 
    86                      e3t_0(ji,jj,jk) = MAX(tabres(ji,jj,jk),MIN(e3zps_min,e3t_1d(jk)*e3zps_rat)) 
    8791                   endif 
    8892                  END DO 
    … …  
    113117          do jj=j1,j2 
    114118          do ji=i1,i2 
    115            if (min(mbkt(ji,jj),mbkt(ji+1,jj))<jk) then 
     119           if (min(mbkt(ji,jj),mbkt(ji+1,jj))<=jk) then 
    116120            tabres(ji,jj,jk) = zrhoy * e2u(ji,jj) * MIN(e3zps_min,e3t_1d(jk)*e3zps_rat) 
    117121           else 
    … …  
    125129            DO jj=j1,j2 
    126130               DO ji=i1,i2 
    127                  if (min(mbkt(ji,jj),mbkt(ji+1,jj))<jk) then 
    128                    e3u_0(ji,jj,jk)=e3t_1d(jk) 
     131                 if (min(mbkt(ji,jj),mbkt(ji+1,jj))<=jk) then 
     132                   e3u_0(ji,jj,jk)=MAX(tabres(ji,jj,jk) / e2u(ji,jj),MIN(e3zps_min,e3t_1d(jk)*e3zps_rat)) 
    129133                 else 
    130                    e3u_0(ji,jj,jk) = MAX(tabres(ji,jj,jk) / e2u(ji,jj),MIN(e3zps_min,e3t_1d(jk)*e3zps_rat)) 
     134                   e3u_0(ji,jj,jk) = e3t_1d(jk) 
    131135                 endif 
    132136               END DO 
    … …  
    157161          do jj=j1,j2 
    158162          do ji=i1,i2 
    159            if (min(mbkt(ji,jj),mbkt(ji,jj+1))<jk) then 
     163           if (min(mbkt(ji,jj),mbkt(ji,jj+1))<=jk) then 
    160164            tabres(ji,jj,jk) = zrhox * e1v(ji,jj) * MIN(e3zps_min,e3t_1d(jk)*e3zps_rat) 
    161165           else 
    … …  
    169173            DO jj=j1,j2 
    170174               DO ji=i1,i2 
    171                  if (min(mbkt(ji,jj),mbkt(ji,jj+1))<jk) then 
    172                    e3v_0(ji,jj,jk)=e3t_1d(jk) 
     175                 if (min(mbkt(ji,jj),mbkt(ji,jj+1))<=jk) then 
     176                   e3v_0(ji,jj,jk)=MAX(tabres(ji,jj,jk) / e1v(ji,jj),MIN(e3zps_min,e3t_1d(jk)*e3zps_rat)) 
    173177                 else 
    174                    e3v_0(ji,jj,jk) = MAX(tabres(ji,jj,jk) / e1v(ji,jj),MIN(e3zps_min,e3t_1d(jk)*e3zps_rat)) 
     178                   e3v_0(ji,jj,jk) = e3t_1d(jk) 
    175179                 endif 
    176180               END DO 

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

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

  • utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/bilinear_interp.f90

    r12414 r13024  
     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_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/dom_oce.F90

    r12414 r13024  
    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 
    … …  
    279281#if defined key_agrif 
    280282   LOGICAL, PUBLIC, PARAMETER ::   lk_agrif = .TRUE.    !: agrif flag 
     283   LOGICAL, PUBLIC            ::   lk_south, lk_north, lk_west, lk_east !: Child grid boundaries (interpolation or not) 
    281284#else 
    282285   LOGICAL, PUBLIC, PARAMETER ::   lk_agrif = .FALSE.   !: agrif flag 
    … …  
    326329         &      ff_f (jpi,jpj) ,    ff_t (jpi,jpj)                                     , STAT=ierr(3) ) 
    327330         ! 
    328       ALLOCATE( gdept_0(jpi,jpj,jpk) , gdepw_0(jpi,jpj,jpk) , STAT=ierr(4) ) 
     331      ! ALLOCATE( gdept_0(jpi,jpj,jpk) , gdepw_0(jpi,jpj,jpk) , gde3w_0(jpi,jpj,jpk) ,      & 
     332      !    &      gdept_b(jpi,jpj,jpk) , gdepw_b(jpi,jpj,jpk) ,                             & 
     333      !    &      gdept_n(jpi,jpj,jpk) , gdepw_n(jpi,jpj,jpk) , gde3w_n(jpi,jpj,jpk) , STAT=ierr(4) ) 
     334 
     335         ALLOCATE( gdept_0(jpi,jpj,jpk) , gdepw_0(jpi,jpj,jpk), STAT=ierr(4) ) 
     336 
     337         ! 
     338      ! 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) ,   & 
     339      !    &      e3t_b(jpi,jpj,jpk) , e3u_b(jpi,jpj,jpk) , e3v_b(jpi,jpj,jpk) ,                      e3w_b(jpi,jpj,jpk) ,   & 
     340      !    &      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) ,   & 
     341      !    &      e3t_a(jpi,jpj,jpk) , e3u_a(jpi,jpj,jpk) , e3v_a(jpi,jpj,jpk) ,                                             & 
     342      !    !                                                          ! 
     343      !    &      e3uw_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk) ,         & 
     344      !    &      e3uw_b(jpi,jpj,jpk) , e3vw_b(jpi,jpj,jpk) ,         & 
     345      !    &      e3uw_n(jpi,jpj,jpk) , e3vw_n(jpi,jpj,jpk) ,     STAT=ierr(5) ) 
     346 
     347         ! 
     348    !  ALLOCATE( ht_0(jpi,jpj) , hu_0(jpi,jpj) , hv_0(jpi,jpj) ,                                           & 
     349    !     &                      hu_b(jpi,jpj) , hv_b(jpi,jpj) , r1_hu_b(jpi,jpj) , r1_hv_b(jpi,jpj) ,     & 
     350    !     &      ht_n(jpi,jpj) , hu_n(jpi,jpj) , hv_n(jpi,jpj) , r1_hu_n(jpi,jpj) , r1_hv_n(jpi,jpj) ,     & 
     351    !     &                      hu_a(jpi,jpj) , hv_a(jpi,jpj) , r1_hu_a(jpi,jpj) , r1_hv_a(jpi,jpj) , STAT=ierr(6)  ) 
    329352         ! 
    330353      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) )                        
     354         &      e3uw_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk) , STAT=ierr(5) ) 
    332355         ! 
    333356      ALLOCATE( gdept_1d(jpk) , e3tp (jpi,jpj), e3wp(jpi,jpj) ,gdepw_1d(jpk) , e3t_1d(jpk) , e3w_1d(jpk) , STAT=ierr(6) ) 
    334357         ! 
    335       ALLOCATE( bathy(jpi,jpj),mbathy(jpi,jpj), tmask_i(jpi,jpj) , tmask_h(jpi,jpj) ,                        &  
     358      ALLOCATE( bathy(jpi,jpj),mbathy(jpi,jpj), tmask_i(jpi,jpj) , tmask_h(jpi,jpj) ,                        & 
    336359         &      ssmask (jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) ,     & 
    337360         &      mbkt   (jpi,jpj) , mbku   (jpi,jpj) , mbkv   (jpi,jpj) , STAT=ierr(7) ) 
    … …  
    340363         &      risfdep(jpi,jpj) , mikv(jpi,jpj) , mikf(jpi,jpj) , STAT=ierr(8) ) 
    341364         ! 
    342       ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk) ,     &  
     365      ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk) ,     & 
    343366         &      vmask(jpi,jpj,jpk) , fmask(jpi,jpj,jpk) , wmask(jpi,jpj,jpk) , STAT=ierr(9) ) 
    344367         ! 
    … …  
    351374      ALLOCATE( msk_opnsea  (jpi,jpj), msk_csundef (jpi,jpj),                        & 
    352375         &      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       ! 
     376         &      msk_csgrpglo(jpi,jpj), msk_csgrprnf(jpi,jpj), msk_csgrpemp(jpi,jpj), STAT=ierr(11) )      ! 
    355377      dom_oce_alloc = MAXVAL(ierr) 
    356378      ! 

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

    r12870 r13024  
    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      !!---------------------------------------------------------------------- 
     69      ! 
     70     ! IF( nn_timing == 1 )   CALL timing_start('dom_init') 
    6471      ! 
    6572      IF(lwp) THEN 
    … …  
    7178      !                       !==  Reference coordinate system  ==! 
    7279      ! 
    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 
     80      CALL dom_nam               ! read namelist ( namrun, namdom ) 
     81                  !   CALL dom_clo               ! Closed seas and lake 
     82          
     83      CALL dom_hgr               ! Horizontal mesh 
     84      CALL dom_zgr( ik_top, ik_bot )  ! Vertical mesh and bathymetry 
     85      CALL dom_msk( ik_top, ik_bot )  ! Masks 
     86      ! 
    8287      ! 
    8388      CALL dom_ctl                  ! print extrema of masked scale factors 
    8489      !  
     90#if ! defined key_agrif 
    8591      CALL cfg_write                ! create the configuration file 
     92#endif 
    8693      ! 
    8794   END SUBROUTINE dom_init 
    … …  
    103110         &             nn_stock, nn_write , ln_mskland  , ln_clobber   , nn_chunksz, nn_euler  ,     & 
    104111         &             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 ,                                & 
     112  !    NAMELIST/namdom/ nn_bathy, cn_topo, cn_bath, cn_lon, cn_lat, nn_interp,                        & 
     113  !       &             rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin,                       & 
     114  !       &             rn_atfp , rn_rdt   , ln_crs      , jphgr_msh ,                                & 
     115  !       &             ppglam0, ppgphi0, ppe1_deg, ppe2_deg, ppe1_m, ppe2_m,                         & 
     116  !       &             ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh,                  & 
     117  !       &             ppa2, ppkth2, ppacr2 
     118 
     119      NAMELIST/namdom/ nn_bathy, cn_topo, cn_bath, cn_lon, cn_lat, rn_scale, nn_interp, & 
     120         &              rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin,            & 
     121         &             rn_atfp , rn_rdt   ,  ln_crs      , jphgr_msh ,                               & 
    108122         &             ppglam0, ppgphi0, ppe1_deg, ppe2_deg, ppe1_m, ppe2_m,                         & 
    109123         &             ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh,                  & 
    110124         &             ppa2, ppkth2, ppacr2 
    111  
    112  
    113125 
    114126      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    … …  
    181193 
    182194 
    183  
     195      !  
     196      cn_topo ='' 
     197      cn_bath ='' 
     198      cn_lon  ='' 
     199      cn_lat  ='' 
     200      rn_scale = 1. 
    184201 
    185202      REWIND( numnam_ref )              ! Namelist namdom in reference namelist : space & time domain (bathymetry, mesh, timestep) 
    … …  
    193210      IF(lwm) WRITE ( numond, namdom ) 
    194211      ! 
     212 
     213 
     214 
    195215      IF(lwp) THEN 
    196216         WRITE(numout,*) 
    … …  
    198218         WRITE(numout,*) '      flag read/compute bathymetry      nn_bathy     = ', nn_bathy 
    199219         IF( nn_bathy == 2 ) THEN 
    200             WRITE(numout,*) '      compute bathymetry from file      cn_topo      = ', cn_topo 
     220            WRITE(numout,*) '   compute bathymetry from file      cn_topo      = ' , cn_topo 
     221            WRITE(numout,*) '   bathymetry name in file           cn_bath      = ' , cn_bath 
     222            WRITE(numout,*) '   longitude name in file            cn_lon       = ' , cn_lon 
     223            WRITE(numout,*) '   latitude  name in file            cn_lat       = ' , cn_lat 
     224            WRITE(numout,*) '   bathmetry scale factor            rn_scale     = ' , rn_scale  
    201225         ENDIF    
    202226         WRITE(numout,*) '      Depth (if =0 bathy=jpkm1)         rn_bathy     = ', rn_bathy 
    … …  
    255279      !!---------------------------------------------------------------------- 
    256280      ! 
     281#undef CHECK_DOM 
     282#ifdef CHECK_DOM 
    257283      IF(lk_mpp) THEN 
    258284         CALL mpp_minloc( 'dom_ctl', e1t(:,:), tmask_i(:,:), ze1min, iloc ) 
    … …  
    292318         WRITE(numout,"(14x,'e2t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2min, iimi2, ijmi2 
    293319      ENDIF 
     320#endif 
    294321      ! 
    295322      ! check that all processes are still there... If some process have an error, 

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

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

  • utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/domclo.F90

    r12442 r13024  
    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_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/domhgr.F90

    r12414 r13024  
    434434      ! ------------------------------------------ 
    435435      ! The equator line must be the latitude coordinate axe 
    436 ! (PM) be carefull with nperio/jperio  
     436 !(PM) be carefull with nperio/jperio  
    437437      IF( jperio == 2 ) THEN 
    438438         znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / REAL( jpi ) 
    … …  
    455455      ! 
    456456      INTEGER ::   inum   ! temporary logical unit 
     457      CHARACTER(LEN=135) :: coordinate_filename 
    457458      !!---------------------------------------------------------------------- 
    458459      ! 
    … …  
    463464      ENDIF 
    464465      ! 
    465       CALL iom_open( 'coordinates', inum ) 
     466      IF (ln_read_cfg) THEN 
     467         coordinate_filename=TRIM(cn_domcfg) 
     468      ELSE 
     469         coordinate_filename='coordinates' 
     470      ENDIF 
     471      CALL iom_open( coordinate_filename, inum ) 
    466472      ! 
    467473      CALL iom_get( inum, jpdom_data, 'glamt', glamt, lrowattr=ln_use_jattr ) 

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

    r12414 r13024  
    2626   USE domisf         ! domain: ice shelf 
    2727   USE domwri         ! domain: write the meshmask file 
     28   USE usrdef_fmask   ! user defined fmask 
    2829   USE bdy_oce        ! open boundary 
    2930   ! 
    … …  
    5253CONTAINS 
    5354 
    54    SUBROUTINE dom_msk 
     55   SUBROUTINE dom_msk( k_top, k_bot ) 
    5556      !!--------------------------------------------------------------------- 
    5657      !!                 ***  ROUTINE dom_msk  *** 
    … …  
    6263      !!      and ko_bot, the indices of the fist and last ocean t-levels which  
    6364      !!      are either defined in usrdef_zgr or read in zgr_read. 
    64       !!                The velocity masks (umask, vmask)  
     65      !!                The velocity masks (umask, vmask, wmask, wumask, wvmask)  
    6566      !!      are deduced from a product of the two neighboring tmask. 
    6667      !!                The vorticity mask (fmask) is deduced from tmask taking 
    … …  
    7778      !!                due to cyclic or North Fold boundaries as well as MPP halos. 
    7879      !! 
    79       !! ** Action :   tmask, umask, vmask, wmask : land/ocean mask  
     80      !! ** Action :   tmask, umask, vmask, wmask, wumask, wvmask : land/ocean mask  
    8081      !!                         at t-, u-, v- w, wu-, and wv-points (=0. or 1.) 
    8182      !!               fmask   : land/ocean mask at f-point (=0., or =1., or  
    … …  
    8586      !!               ssmask , ssumask, ssvmask, ssfmask : 2D ocean mask 
    8687      !!---------------------------------------------------------------------- 
     88 
     89      INTEGER, DIMENSION(:,:), INTENT(in) ::   k_top, k_bot   ! first and last ocean level 
    8790      ! 
    8891      INTEGER  ::   ji, jj, jk     ! dummy loop indices 
    … …  
    133136     ! N.B. tmask has already the right boundary conditions since mbathy is ok 
    134137     ! 
     138!      tmask(:,:,:) = 0._wp 
     139!      DO jk = 1, jpk 
     140!         DO jj = 1, jpj 
     141!            DO ji = 1, jpi 
     142!               IF(      ( REAL( mbathy (ji,jj) - jk, wp ) + 0.1_wp >= 0._wp )         & 
     143!               &  .AND. ( REAL( misfdep(ji,jj) - jk, wp ) - 0.1_wp <= 0._wp ) ) THEN 
     144!                  tmask(ji,jj,jk) = 1._wp 
     145!               END IF   
     146!            END DO 
     147!         END DO 
     148!      END DO     
     149  
     150!      IF ( ln_isfsubgl ) CALL zgr_isf_subgl 
     151 
     152      !  Ocean/land mask at t-point  (computed from ko_top and ko_bot) 
     153      ! ---------------------------- 
     154      ! 
    135155      tmask(:,:,:) = 0._wp 
    136       DO jk = 1, jpk 
     156      IF( ln_read_cfg) THEN 
    137157         DO jj = 1, jpj 
    138158            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   
     159               iktop = k_top(ji,jj) 
     160               ikbot = k_bot(ji,jj) 
     161               IF( iktop /= 0 ) THEN       ! water in the column 
     162                  tmask(ji,jj,iktop:ikbot  ) = 1._wp 
     163               ENDIF 
     164            END DO   
     165         END DO   
     166         ELSE 
     167         DO jk = 1, jpk 
     168            DO jj = 1, jpj 
     169               DO ji = 1, jpi 
     170                  IF(      ( REAL( mbathy (ji,jj) - jk, wp ) + 0.1_wp >= 0._wp )         & 
     171                  &  .AND. ( REAL( misfdep(ji,jj) - jk, wp ) - 0.1_wp <= 0._wp ) ) THEN 
     172                     tmask(ji,jj,jk) = 1._wp 
     173                  END IF 
     174               END DO 
    143175            END DO 
    144176         END DO 
    145       END DO     
    146   
    147       IF ( ln_isfsubgl ) CALL zgr_isf_subgl 
     177         IF ( ln_isfsubgl ) CALL zgr_isf_subgl 
     178      ENDIF 
     179 
    148180 
    149181!SF  add here lbc_lnk: bug not still understood : cause now domain configuration is read ! 
    … …  
    272304#if defined key_agrif  
    273305            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  
     306               IF(lk_east) fmask(nlci-1 , :     ,jk) = 0.e0      ! east  
     307               IF(lk_west) fmask(1      , :     ,jk) = 0.e0      ! west  
     308               IF(lk_north) fmask(:      ,nlcj-1 ,jk) = 0.e0      ! north  
     309               IF(lk_south) fmask(:      ,1      ,jk) = 0.e0      ! south  
    278310            ENDIF  
    279311#endif  
    … …  
    287319         ! 
    288320      ENDIF 
    289       ! 
    290       ! write out mesh mask 
    291       IF ( nn_msh > 0 ) CALL dom_wri 
     321       
     322      ! User defined alteration of fmask (use to reduce ocean transport in specified straits) 
     323      ! --------------------------------  
     324      ! 
     325 
     326      CALL usr_def_fmask( cn_cfg, nn_cfg, fmask ) 
    292327      ! 
    293328   END SUBROUTINE dom_msk 

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

    r12414 r13024  
    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_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/domutl.F90

    r12414 r13024  
    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 
    … …  
    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 
    … …  
    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 

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

    r12414 r13024  
    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 
    … …  
    9293CONTAINS        
    9394 
    94    SUBROUTINE dom_zgr 
     95   SUBROUTINE dom_zgr( k_top, k_bot ) 
    9596      !!---------------------------------------------------------------------- 
    9697      !!                ***  ROUTINE dom_zgr  *** 
    … …  
    109110      !! ** Action  :   define gdep., e3., mbathy and bathy 
    110111      !!---------------------------------------------------------------------- 
     112      INTEGER, DIMENSION(:,:), INTENT(out) ::   k_top, k_bot   ! ocean first and last level indices 
     113      ! 
    111114      INTEGER ::   ioptio, ibat   ! local integer 
    112115      INTEGER ::   ios 
    113116      ! 
     117      INTEGER :: jk 
     118      REAL(wp) ::   zrefdep             ! depth of the reference level (~10m) 
     119 
     120 
    114121      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav, ln_linssh 
    115122      !!---------------------------------------------------------------------- 
    … …  
    124131902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in configuration namelist', lwp ) 
    125132      IF(lwm) WRITE ( numond, namzgr ) 
     133 
     134      IF(ln_read_cfg) THEN 
     135         IF(lwp) WRITE(numout,*) 
     136         IF(lwp) WRITE(numout,*) '   ==>>>   Read vertical mesh in ', TRIM( cn_domcfg ), ' file' 
     137         ! 
     138         CALL zgr_read   ( ln_zco  , ln_zps  , ln_sco, ln_isfcav,   & 
     139            &              gdept_1d, gdepw_1d, e3t_1d, e3w_1d   ,   &    ! 1D gridpoints depth 
     140            &              gdept_0 , gdepw_0                    ,   &    ! gridpoints depth 
     141            &              e3t_0   , e3u_0   , e3v_0 , e3f_0    ,   &    ! vertical scale factors 
     142            &              e3w_0   , e3uw_0  , e3vw_0           ,   &    ! vertical scale factors 
     143            &              k_top   , k_bot            )                  ! 1st & last ocean level 
     144            ! 
     145!!gm to be remove when removing the OLD definition of e3 scale factors so that gde3w disappears 
     146!      ! Compute gde3w_0 (vertical sum of e3w) 
     147!      gde3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 
     148!      DO jk = 2, jpk 
     149!         gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk) 
     150!      END DO 
     151      ! 
     152 
     153      !                                ! top/bottom ocean level indices for t-, u- and v-points (f-point also for top) 
     154      CALL zgr_top_bot( k_top, k_bot )      ! with a minimum value set to 1 
     155 
     156      !                                ! deepest/shallowest W level Above/Below ~10m 
     157!!gm BUG in s-coordinate this does not work! 
     158      zrefdep = 10._wp - 0.1_wp * MINVAL( e3w_1d )                   ! ref. depth with tolerance (10% of minimum layer thickness) 
     159      nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m 
     160      nla10 = nlb10 - 1                                              ! deepest    W level Above ~10m 
     161!!gm end bug 
     162 
     163      ENDIF         
    126164 
    127165      IF(lwp) THEN                     ! Control print 
    … …  
    150188      IF ( ln_sco .AND. ln_isfcav ) ioptio = ioptio + 1 
    151189      IF( ioptio > 0 )   CALL ctl_stop( ' Cavity not tested/compatible with full step (zco) and sigma (ln_sco) ' ) 
     190 
     191      IF(.NOT.ln_read_cfg) THEN 
    152192      ! 
    153193      ! Build the vertical coordinate system 
    … …  
    163203      ! ----------------------------------- 
    164204                          CALL zgr_bot_level    ! deepest ocean level for t-, u- and v-points 
     205                          k_bot = mbkt 
    165206                          CALL zgr_top_level    ! shallowest ocean level for T-, U-, V- points 
     207                          k_top = mikt 
     208 
     209      ENDIF 
    166210      ! 
    167211      IF( nprint == 1 .AND. lwp )   THEN 
    168          WRITE(numout,*) ' MIN val mbathy  ', MINVAL(  mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 
     212         WRITE(numout,*) ' MIN val k_top   ', MINVAL(   k_top(:,:) ), ' MAX ', MAXVAL( k_top(:,:) ) 
     213         WRITE(numout,*) ' MIN val k_bot   ', MINVAL(   k_bot(:,:) ), ' MAX ', MAXVAL( k_bot(:,:) ) 
    169214         WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ),   & 
    170215            &                          ' w ', MINVAL( gdepw_0(:,:,:) ) 
    … …  
    181226            &                          ' w ', MAXVAL(   e3w_0(:,:,:) ) 
    182227      ENDIF 
     228 
    183229      ! 
    184230   END SUBROUTINE dom_zgr 
    185231 
     232   SUBROUTINE zgr_read( ld_zco  , ld_zps  , ld_sco  , ld_isfcav,   &   ! type of vertical coordinate 
     233      &                 pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d  ,   &   ! 1D reference vertical coordinate 
     234      &                 pdept , pdepw ,                            &   ! 3D t & w-points depth 
     235      &                 pe3t  , pe3u  , pe3v   , pe3f ,            &   ! vertical scale factors 
     236      &                 pe3w  , pe3uw , pe3vw         ,            &   !     -      -      - 
     237      &                 k_top  , k_bot    )                            ! top & bottom ocean level 
     238      !!--------------------------------------------------------------------- 
     239      !!              ***  ROUTINE zgr_read  *** 
     240      !! 
     241      !! ** Purpose :   Read the vertical information in the domain configuration file 
     242      !! 
     243      !!---------------------------------------------------------------------- 
     244      LOGICAL                   , INTENT(out) ::   ld_zco, ld_zps, ld_sco      ! vertical coordinate flags 
     245      LOGICAL                   , INTENT(out) ::   ld_isfcav                   ! under iceshelf cavity flag 
     246      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pdept_1d, pdepw_1d          ! 1D grid-point depth       [m] 
     247      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pe3t_1d , pe3w_1d           ! 1D vertical scale factors [m] 
     248      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pdept, pdepw                ! grid-point depth          [m] 
     249      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pe3t , pe3u , pe3v , pe3f   ! vertical scale factors    [m] 
     250      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pe3w , pe3uw, pe3vw         !    -       -      - 
     251      INTEGER , DIMENSION(:,:)  , INTENT(out) ::   k_top , k_bot               ! first & last ocean level 
     252      ! 
     253      INTEGER  ::   jk     ! dummy loop index 
     254      INTEGER  ::   inum   ! local logical unit 
     255      REAL(WP) ::   z_zco, z_zps, z_sco, z_cav 
     256      REAL(wp), DIMENSION(jpi,jpj) ::   z2d   ! 2D workspace 
     257      !!---------------------------------------------------------------------- 
     258      ! 
     259      IF(lwp) THEN 
     260         WRITE(numout,*) 
     261         WRITE(numout,*) '   zgr_read : read the vertical coordinates in ', TRIM( cn_domcfg ), ' file' 
     262         WRITE(numout,*) '   ~~~~~~~~' 
     263      ENDIF 
     264      ! 
     265      CALL iom_open( cn_domcfg, inum ) 
     266      ! 
     267      !                          !* type of vertical coordinate 
     268      CALL iom_get( inum, 'ln_zco'   , z_zco ) 
     269      CALL iom_get( inum, 'ln_zps'   , z_zps ) 
     270      CALL iom_get( inum, 'ln_sco'   , z_sco ) 
     271      IF( z_zco == 0._wp ) THEN   ;   ld_zco = .false.   ;   ELSE   ;   ld_zco = .true.   ;   ENDIF 
     272      IF( z_zps == 0._wp ) THEN   ;   ld_zps = .false.   ;   ELSE   ;   ld_zps = .true.   ;   ENDIF 
     273      IF( z_sco == 0._wp ) THEN   ;   ld_sco = .false.   ;   ELSE   ;   ld_sco = .true.   ;   ENDIF 
     274      ! 
     275      !                          !* ocean cavities under iceshelves 
     276      CALL iom_get( inum, 'ln_isfcav', z_cav ) 
     277      IF( z_cav == 0._wp ) THEN   ;   ld_isfcav = .false.   ;   ELSE   ;   ld_isfcav = .true.   ;   ENDIF 
     278      ! 
     279      !                          !* vertical scale factors 
     280      CALL iom_get( inum, jpdom_unknown, 'e3t_1d'  , pe3t_1d  )                     ! 1D reference coordinate 
     281      CALL iom_get( inum, jpdom_unknown, 'e3w_1d'  , pe3w_1d  ) 
     282      ! 
     283      CALL iom_get( inum, jpdom_data, 'e3t_0'  , pe3t  , lrowattr=ln_use_jattr )    ! 3D coordinate 
     284      CALL iom_get( inum, jpdom_data, 'e3u_0'  , pe3u  , lrowattr=ln_use_jattr ) 
     285      CALL iom_get( inum, jpdom_data, 'e3v_0'  , pe3v  , lrowattr=ln_use_jattr ) 
     286      CALL iom_get( inum, jpdom_data, 'e3f_0'  , pe3f  , lrowattr=ln_use_jattr ) 
     287      CALL iom_get( inum, jpdom_data, 'e3w_0'  , pe3w  , lrowattr=ln_use_jattr ) 
     288      CALL iom_get( inum, jpdom_data, 'e3uw_0' , pe3uw , lrowattr=ln_use_jattr ) 
     289      CALL iom_get( inum, jpdom_data, 'e3vw_0' , pe3vw , lrowattr=ln_use_jattr ) 
     290      ! 
     291      !                          !* depths 
     292      !                                   !- old depth definition (obsolescent feature) 
     293      IF(  iom_varid( inum, 'gdept_1d', ldstop = .FALSE. ) > 0  .AND.  & 
     294         & iom_varid( inum, 'gdepw_1d', ldstop = .FALSE. ) > 0  .AND.  & 
     295         & iom_varid( inum, 'gdept_0' , ldstop = .FALSE. ) > 0  .AND.  & 
     296         & iom_varid( inum, 'gdepw_0' , ldstop = .FALSE. ) > 0    ) THEN 
     297         CALL ctl_warn( 'zgr_read : old definition of depths and scale factors used ', &  
     298            &           '           depths at t- and w-points read in the domain configuration file') 
     299         CALL iom_get( inum, jpdom_unknown, 'gdept_1d', pdept_1d )    
     300         CALL iom_get( inum, jpdom_unknown, 'gdepw_1d', pdepw_1d ) 
     301         CALL iom_get( inum, jpdom_data   , 'gdept_0' , pdept , lrowattr=ln_use_jattr ) 
     302         CALL iom_get( inum, jpdom_data   , 'gdepw_0' , pdepw , lrowattr=ln_use_jattr ) 
     303         ! 
     304      ELSE                                !- depths computed from e3. scale factors 
     305         CALL e3_to_depth( pe3t_1d, pe3w_1d, pdept_1d, pdepw_1d )    ! 1D reference depth 
     306         CALL e3_to_depth( pe3t   , pe3w   , pdept   , pdepw    )    ! 3D depths 
     307         IF(lwp) THEN 
     308            WRITE(numout,*) 
     309            WRITE(numout,*) '              Reference 1D z-coordinate depth and scale factors:' 
     310            WRITE(numout, "(9x,' level  gdept_1d  gdepw_1d  e3t_1d   e3w_1d  ')" ) 
     311            WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, pdept_1d(jk), pdepw_1d(jk), pe3t_1d(jk), pe3w_1d(jk), jk = 1, jpk ) 
     312         ENDIF 
     313      ENDIF 
     314      ! 
     315      !                          !* ocean top and bottom level 
     316      CALL iom_get( inum, jpdom_data, 'top_level'    , z2d  , lrowattr=ln_use_jattr )   ! 1st wet T-points (ISF) 
     317      k_top(:,:) = NINT( z2d(:,:) ) 
     318      CALL iom_get( inum, jpdom_data, 'bottom_level' , z2d  , lrowattr=ln_use_jattr )   ! last wet T-points 
     319      k_bot(:,:) = NINT( z2d(:,:) ) 
     320      ! 
     321      ! reference depth for negative bathy (wetting and drying only) 
     322      ! IF( ll_wd )  CALL iom_get( inum,  'rn_wd_ref_depth' , ssh_ref   ) 
     323      ! 
     324      CALL iom_close( inum ) 
     325      ! 
     326   END SUBROUTINE zgr_read 
     327 
     328   SUBROUTINE zgr_top_bot( k_top, k_bot ) 
     329      !!---------------------------------------------------------------------- 
     330      !!                    ***  ROUTINE zgr_top_bot  *** 
     331      !! 
     332      !! ** Purpose :   defines the vertical index of ocean bottom (mbk. arrays) 
     333      !! 
     334      !! ** Method  :   computes from k_top and k_bot with a minimum value of 1 over land 
     335      !! 
     336      !! ** Action  :   mikt, miku, mikv :   vertical indices of the shallowest  
     337      !!                                     ocean level at t-, u- & v-points 
     338      !!                                     (min value = 1) 
     339      !! ** Action  :   mbkt, mbku, mbkv :   vertical indices of the deeptest  
     340      !!                                     ocean level at t-, u- & v-points 
     341      !!                                     (min value = 1 over land) 
     342      !!---------------------------------------------------------------------- 
     343      INTEGER , DIMENSION(:,:), INTENT(in) ::   k_top, k_bot   ! top & bottom ocean level indices 
     344      ! 
     345      INTEGER ::   ji, jj   ! dummy loop indices 
     346      REAL(wp), DIMENSION(jpi,jpj) ::   zk   ! workspace 
     347      !!---------------------------------------------------------------------- 
     348      ! 
     349      IF(lwp) WRITE(numout,*) 
     350      IF(lwp) WRITE(numout,*) '    zgr_top_bot : ocean top and bottom k-index of T-, U-, V- and W-levels ' 
     351      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~' 
     352      ! 
     353      mikt(:,:) = MAX( k_top(:,:) , 1 )    ! top    ocean k-index of T-level (=1 over land) 
     354      ! 
     355      mbkt(:,:) = MAX( k_bot(:,:) , 1 )    ! bottom ocean k-index of T-level (=1 over land) 
     356  
     357      !                                    ! N.B.  top     k-index of W-level = mikt 
     358      !                                    !       bottom  k-index of W-level = mbkt+1 
     359      DO jj = 1, jpjm1 
     360         DO ji = 1, jpim1 
     361            miku(ji,jj) = MAX(  mikt(ji+1,jj  ) , mikt(ji,jj)  ) 
     362            mikv(ji,jj) = MAX(  mikt(ji  ,jj+1) , mikt(ji,jj)  ) 
     363            mikf(ji,jj) = MAX(  mikt(ji  ,jj+1) , mikt(ji,jj), mikt(ji+1,jj  ), mikt(ji+1,jj+1)  ) 
     364            ! 
     365            mbku(ji,jj) = MIN(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  ) 
     366            mbkv(ji,jj) = MIN(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  ) 
     367         END DO 
     368      END DO 
     369      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk  
     370      zk(:,:) = REAL( miku(:,:), wp )   ;   CALL lbc_lnk( 'domzgr', zk, 'U', 1. )   ;   miku(:,:) = MAX( NINT( zk(:,:) ), 1 ) 
     371      zk(:,:) = REAL( mikv(:,:), wp )   ;   CALL lbc_lnk( 'domzgr', zk, 'V', 1. )   ;   mikv(:,:) = MAX( NINT( zk(:,:) ), 1 ) 
     372      zk(:,:) = REAL( mikf(:,:), wp )   ;   CALL lbc_lnk( 'domzgr', zk, 'F', 1. )   ;   mikf(:,:) = MAX( NINT( zk(:,:) ), 1 ) 
     373      ! 
     374      zk(:,:) = REAL( mbku(:,:), wp )   ;   CALL lbc_lnk( 'domzgr', zk, 'U', 1. )   ;   mbku(:,:) = MAX( NINT( zk(:,:) ), 1 ) 
     375      zk(:,:) = REAL( mbkv(:,:), wp )   ;   CALL lbc_lnk( 'domzgr', zk, 'V', 1. )   ;   mbkv(:,:) = MAX( NINT( zk(:,:) ), 1 ) 
     376      ! 
     377   END SUBROUTINE zgr_top_bot 
    186378 
    187379   SUBROUTINE zgr_z 
    … …  
    665857      ENDIF 
    666858 
    667       zbathy(:,:) = FLOAT( mbathy(:,:) ) 
    668       CALL lbc_lnk( 'domzgr',zbathy, 'T', 1._wp ) 
    669       mbathy(:,:) = INT( zbathy(:,:) ) 
    670  
     859      IF( lk_mpp ) THEN 
     860         zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     861         CALL lbc_lnk( 'toto',zbathy, 'T', 1._wp ) 
     862         mbathy(:,:) = INT( zbathy(:,:) ) 
     863      ENDIF 
    671864      !                                          ! East-west cyclic boundary conditions 
    672865      IF( jperio == 0 ) THEN 
    … …  
    10611254   END SUBROUTINE zgr_zps 
    10621255 
     1256 
    10631257   SUBROUTINE zgr_sco 
    10641258      !!---------------------------------------------------------------------- 
    … …  
    12621456         END DO 
    12631457         ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero 
    1264          CALL lbc_lnk( 'domzgr',zenv, 'T', 1._wp, 'no0' ) 
     1458         CALL lbc_lnk( 'toto',zenv, 'T', 1._wp, 'no0' ) 
    12651459         !                                                  ! ================ ! 
    12661460      END DO                                                !     End loop     ! 

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

    r12414 r13024  
    530530      INTEGER               :: idim3                ! id of the third dimension 
    531531      ! 
    532       INTEGER ::   nldi_save, nlei_save    !:patch before we remove periodicity and close boundaries in output files 
    533       INTEGER ::   nldj_save, nlej_save    !: 
     532    !  INTEGER ::   nldi_save, nlei_save    !:patch before we remove periodicity and close boundaries in output files 
     533    !  INTEGER ::   nldj_save, nlej_save    !: 
    534534      !--------------------------------------------------------------------- 
    535535      ! 
    … …  
    540540      ! without this, issue in some model decomposition 
    541541      ! seb: patch before we remove periodicity and close boundaries in output files 
    542       nldi_save = nldi   ;   nlei_save = nlei 
    543       nldj_save = nldj   ;   nlej_save = nlej 
    544       IF( nimpp           ==      1 ) nldi = 1 
    545       IF( nimpp + jpi - 1 == jpiglo ) nlei = jpi 
    546       IF( njmpp           ==      1 ) nldj = 1 
    547       IF( njmpp + jpj - 1 == jpjglo ) nlej = jpj 
     542    !  nldi_save = nldi   ;   nlei_save = nlei 
     543    !  nldj_save = nldj   ;   nlej_save = nlej 
     544    !  IF( nimpp           ==      1 ) nldi = 1 
     545    !  IF( nimpp + jpi - 1 == jpiglo ) nlei = jpi 
     546    !  IF( njmpp           ==      1 ) nldj = 1 
     547    !  IF( njmpp + jpj - 1 == jpjglo ) nlej = jpj 
    548548      ! 
    549549      ! define dimension variables if it is not already done 
    … …  
    719719      ENDIF 
    720720      ! 
    721       nldi = nldi_save   ;   nlei = nlei_save 
    722       nldj = nldj_save   ;   nlej = nlej_save 
     721  !    nldi = nldi_save   ;   nlei = nlei_save 
     722  !    nldj = nldj_save   ;   nlej = nlej_save 
    723723      !      
    724724   END SUBROUTINE iom_nf90_rp0123d 

  • utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/lib_mpp.F90

    r12414 r13024  
    16621662      CHARACTER(len=*),DIMENSION(:) ::   ldtxt 
    16631663      CHARACTER(len=*) ::   ldname 
    1664       INTEGER ::   function_value 
    16651664      INTEGER ::   kumnam_ref, knumnam_cfg , kumond , kstop 
    16661665      IF( PRESENT( localComm ) ) mpi_comm_oce = localComm 

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

    r12414 r13024  
    1717#      define MPI_OPERATION mpi_maxloc 
    1818#      define LOC_OPERATION MAXLOC 
    19 #      define ERRVAL -HUGE 
    2019#   endif 
    2120#   if defined OPERATION_MINLOC 
    2221#      define MPI_OPERATION mpi_minloc 
    2322#      define LOC_OPERATION MINLOC 
    24 #      define ERRVAL HUGE 
    2523#   endif 
    2624 
    … …  
    4442      ! 
    4543      idim = SIZE(kindex) 
     44      ALLOCATE ( ilocs(idim) ) 
    4645      ! 
    47       IF ( ALL(MASK_IN(:,:,:) /= 1._wp) ) THEN 
    48          ! special case for land processors 
    49          zmin = ERRVAL(zmin) 
    50          index0 = 0 
    51       ELSE 
    52          ALLOCATE ( ilocs(idim) ) 
    53          ! 
    54          ilocs = LOC_OPERATION( ARRAY_IN(:,:,:) , mask= MASK_IN(:,:,:) == 1._wp ) 
    55          zmin  = ARRAY_IN(ilocs(1),ilocs(2),ilocs(3)) 
    56          ! 
    57          kindex(1) = mig( ilocs(1) ) 
     46      ilocs = LOC_OPERATION( ARRAY_IN(:,:,:) , mask= MASK_IN(:,:,:) == 1._wp ) 
     47      zmin  = ARRAY_IN(ilocs(1),ilocs(2),ilocs(3)) 
     48      ! 
     49      kindex(1) = ilocs(1) + nimpp - 1 
    5850#  if defined DIM_2d || defined DIM_3d    /* avoid warning when kindex has 1 element */ 
    59          kindex(2) = mjg( ilocs(2) ) 
     51      kindex(2) = ilocs(2) + njmpp - 1 
    6052#  endif 
    6153#  if defined DIM_3d                      /* avoid warning when kindex has 2 elements */ 
    62          kindex(3) = ilocs(3) 
     54      kindex(3) = ilocs(3) 
    6355#  endif 
    64          !  
    65          DEALLOCATE (ilocs) 
    66          ! 
    67          index0 = kindex(1)-1   ! 1d index starting at 0 
     56      !  
     57      DEALLOCATE (ilocs) 
     58      ! 
     59      index0 = kindex(1)-1   ! 1d index starting at 0 
    6860#  if defined DIM_2d || defined DIM_3d   /* avoid warning when kindex has 1 element */ 
    69          index0 = index0 + jpiglo * (kindex(2)-1) 
     61      index0 = index0 + jpiglo * (kindex(2)-1) 
    7062#  endif 
    7163#  if defined DIM_3d                     /* avoid warning when kindex has 2 elements */ 
    72          index0 = index0 + jpiglo * jpjglo * (kindex(3)-1) 
     64      index0 = index0 + jpiglo * jpjglo * (kindex(3)-1) 
    7365#  endif 
    74       END IF 
    7566      zain(1,:) = zmin 
    7667      zain(2,:) = REAL(index0, wp) 
    … …  
    10798#undef LOC_OPERATION 
    10899#undef INDEX_TYPE 
    109 #undef ERRVAL 

  • utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/mppini.F90

    r12414 r13024  
    9999         CALL ctl_stop( 'mpp_init: equality  jpni = jpnj = jpnij = 1 is not satisfied',   & 
    100100            &           'the domain is lay out for distributed memory computing!' ) 
     101 
     102#if defined key_agrif 
     103     IF (.not.agrif_root()) THEN 
     104        CALL agrif_nemo_init 
     105     ENDIF 
     106 
     107      IF( .NOT. Agrif_Root() ) THEN       ! AGRIF children: specific setting (cf. agrif_user.F90) 
     108         print *,'nbcellsx = ',nbcellsx,nbghostcells_x 
     109         print *,'nbcellsy = ',nbcellsy,nbghostcells_y_s,nbghostcells_y_n 
     110         IF( jpiglo /= nbcellsx + 2 + 2*nbghostcells_x ) THEN 
     111            IF(lwp) THEN 
     112               WRITE(numout,*) 
     113               WRITE(numout,*) 'jpiglo should be: ', nbcellsx + 2 + 2*nbghostcells_x 
     114            ENDIF         
     115            CALL ctl_stop( 'STOP', 'mpp_init: Agrif children requires jpiglo == nbcellsx + 2 + 2*nbghostcells_x' ) 
     116         ENDIF    
     117         IF( jpjglo /= nbcellsy + 2 + nbghostcells_y_s + nbghostcells_y_n ) THEN 
     118            IF(lwp) THEN 
     119               WRITE(numout,*) 
     120               WRITE(numout,*) 'jpjglo shoud be: ', nbcellsy + 2 + nbghostcells_y_s + nbghostcells_y_n 
     121            ENDIF         
     122            CALL ctl_stop( 'STOP', & 
     123                'mpp_init: Agrif children requires jpjglo == nbcellsy + 2 + nbghostcells_y_s + nbghostcells_y_n' ) 
     124         ENDIF    
     125         IF( ln_use_jattr )   CALL ctl_stop( 'STOP', 'mpp_init:Agrif children requires ln_use_jattr = .false. ' ) 
     126      ENDIF 
     127#endif 
    101128         ! 
    102129   END SUBROUTINE mpp_init 
    … …  
    274301            IF(lwp) THEN 
    275302               WRITE(numout,*) 
    276                WRITE(numout,*) 'jpiglo shoud be: ', nbcellsx + 2 + 2*nbghostcells 
     303               WRITE(numout,*) 'jpiglo should be: ', nbcellsx + 2 + 2*nbghostcells 
    277304            ENDIF         
    278305            CALL ctl_stop( 'STOP', 'mpp_init: Agrif children requires jpiglo == nbcellsx + 2 + 2*nbghostcells' ) 

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

    r12414 r13024  
    408408      !!---------------------------------------------------------------------- 
    409409      ! 
    410       ierr = dom_oce_alloc   ()          ! ocean domain 
     410      ierr = 0 
     411      ierr = ierr + dom_oce_alloc   ()          ! ocean domain 
    411412      ! 
    412413      CALL mpp_sum( 'nemogcm', ierr ) 
    … …  
    416417 
    417418 
     419   SUBROUTINE nemo_partition( num_pes ) 
     420      !!---------------------------------------------------------------------- 
     421      !!                 ***  ROUTINE nemo_partition  *** 
     422      !! 
     423      !! ** Purpose : 
     424      !! 
     425      !! ** Method  : 
     426      !!---------------------------------------------------------------------- 
     427      INTEGER, INTENT(in) ::   num_pes   ! The number of MPI processes we have 
     428      ! 
     429      INTEGER, PARAMETER :: nfactmax = 20 
     430      INTEGER :: nfact ! The no. of factors returned 
     431      INTEGER :: ierr  ! Error flag 
     432      INTEGER :: ji 
     433      INTEGER :: idiff, mindiff, imin ! For choosing pair of factors that are closest in value 
     434      INTEGER, DIMENSION(nfactmax) :: ifact ! Array of factors 
     435      !!---------------------------------------------------------------------- 
     436      ! 
     437      ierr = 0 
     438      ! 
     439      CALL factorise( ifact, nfactmax, nfact, num_pes, ierr ) 
     440      ! 
     441      IF( nfact <= 1 ) THEN 
     442         WRITE (numout, *) 'WARNING: factorisation of number of PEs failed' 
     443         WRITE (numout, *) '       : using grid of ',num_pes,' x 1' 
     444         jpnj = 1 
     445         jpni = num_pes 
     446      ELSE 
     447         ! Search through factors for the pair that are closest in value 
     448         mindiff = 1000000 
     449         imin    = 1 
     450         DO ji = 1, nfact-1, 2 
     451            idiff = ABS( ifact(ji) - ifact(ji+1) ) 
     452            IF( idiff < mindiff ) THEN 
     453               mindiff = idiff 
     454               imin = ji 
     455            ENDIF 
     456         END DO 
     457         jpnj = ifact(imin) 
     458         jpni = ifact(imin + 1) 
     459      ENDIF 
     460      ! 
     461      jpnij = jpni*jpnj 
     462      ! 
     463   END SUBROUTINE nemo_partition 
     464 
     465 
     466   SUBROUTINE factorise( kfax, kmaxfax, knfax, kn, kerr ) 
     467      !!---------------------------------------------------------------------- 
     468      !!                     ***  ROUTINE factorise  *** 
     469      !! 
     470      !! ** Purpose :   return the prime factors of n. 
     471      !!                knfax factors are returned in array kfax which is of 
     472      !!                maximum dimension kmaxfax. 
     473      !! ** Method  : 
     474      !!---------------------------------------------------------------------- 
     475      INTEGER                    , INTENT(in   ) ::   kn, kmaxfax 
     476      INTEGER                    , INTENT(  out) ::   kerr, knfax 
     477      INTEGER, DIMENSION(kmaxfax), INTENT(  out) ::   kfax 
     478      ! 
     479      INTEGER :: ifac, jl, inu 
     480      INTEGER, PARAMETER :: ntest = 14 
     481      INTEGER, DIMENSION(ntest) ::   ilfax 
     482      !!---------------------------------------------------------------------- 
     483      ! 
     484      ! lfax contains the set of allowed factors. 
     485      ilfax(:) = (/(2**jl,jl=ntest,1,-1)/) 
     486      ! 
     487      ! Clear the error flag and initialise output vars 
     488      kerr  = 0 
     489      kfax  = 1 
     490      knfax = 0 
     491      ! 
     492      ! Find the factors of n. 
     493      IF( kn == 1 )   GOTO 20 
     494 
     495      ! nu holds the unfactorised part of the number. 
     496      ! knfax holds the number of factors found. 
     497      ! l points to the allowed factor list. 
     498      ! ifac holds the current factor. 
     499      ! 
     500      inu   = kn 
     501      knfax = 0 
     502      ! 
     503      DO jl = ntest, 1, -1 
     504         ! 
     505         ifac = ilfax(jl) 
     506         IF( ifac > inu )   CYCLE 
     507 
     508         ! Test whether the factor will divide. 
     509 
     510         IF( MOD(inu,ifac) == 0 ) THEN 
     511            ! 
     512            knfax = knfax + 1            ! Add the factor to the list 
     513            IF( knfax > kmaxfax ) THEN 
     514               kerr = 6 
     515               write (*,*) 'FACTOR: insufficient space in factor array ', knfax 
     516               return 
     517            ENDIF 
     518            kfax(knfax) = ifac 
     519            ! Store the other factor that goes with this one 
     520            knfax = knfax + 1 
     521            kfax(knfax) = inu / ifac 
     522            !WRITE (*,*) 'ARPDBG, factors ',knfax-1,' & ',knfax,' are ', kfax(knfax-1),' and ',kfax(knfax) 
     523         ENDIF 
     524         ! 
     525      END DO 
     526      ! 
     527   20 CONTINUE      ! Label 20 is the exit point from the factor search loop. 
     528      ! 
     529   END SUBROUTINE factorise 
     530 
     531 
     532   SUBROUTINE nemo_northcomms 
     533      !!---------------------------------------------------------------------- 
     534      !!                     ***  ROUTINE  nemo_northcomms  *** 
     535      !! ** Purpose :   Setup for north fold exchanges with explicit  
     536      !!                point-to-point messaging 
     537      !! 
     538      !! ** Method :   Initialization of the northern neighbours lists. 
     539      !!---------------------------------------------------------------------- 
     540      !!    1.0  ! 2011-10  (A. C. Coward, NOCS & J. Donners, PRACE) 
     541      !!    2.0  ! 2013-06 Setup avoiding MPI communication (I. Epicoco, S. Mocavero, CMCC)  
     542      !!---------------------------------------------------------------------- 
     543      INTEGER  ::   sxM, dxM, sxT, dxT, jn 
     544      INTEGER  ::   njmppmax 
     545      !!---------------------------------------------------------------------- 
     546      ! 
     547      njmppmax = MAXVAL( njmppt ) 
     548      ! 
     549      !initializes the north-fold communication variables 
     550      isendto(:) = 0 
     551      nsndto     = 0 
     552      ! 
     553      !if I am a process in the north 
     554      IF ( njmpp == njmppmax ) THEN 
     555          !sxM is the first point (in the global domain) needed to compute the 
     556          !north-fold for the current process 
     557          sxM = jpiglo - nimppt(narea) - nlcit(narea) + 1 
     558          !dxM is the last point (in the global domain) needed to compute the 
     559          !north-fold for the current process 
     560          dxM = jpiglo - nimppt(narea) + 2 
     561 
     562          !loop over the other north-fold processes to find the processes 
     563          !managing the points belonging to the sxT-dxT range 
     564   
     565          DO jn = 1, jpni 
     566                !sxT is the first point (in the global domain) of the jn 
     567                !process 
     568                sxT = nfiimpp(jn, jpnj) 
     569                !dxT is the last point (in the global domain) of the jn 
     570                !process 
     571                dxT = nfiimpp(jn, jpnj) + nfilcit(jn, jpnj) - 1 
     572                IF ((sxM .gt. sxT) .AND. (sxM .lt. dxT)) THEN 
     573                   nsndto = nsndto + 1 
     574                     isendto(nsndto) = jn 
     575                ELSEIF ((sxM .le. sxT) .AND. (dxM .ge. dxT)) THEN 
     576                   nsndto = nsndto + 1 
     577                     isendto(nsndto) = jn 
     578                ELSEIF ((dxM .lt. dxT) .AND. (sxT .lt. dxM)) THEN 
     579                   nsndto = nsndto + 1 
     580                     isendto(nsndto) = jn 
     581                END IF 
     582          END DO 
     583          nfsloop = 1 
     584          nfeloop = nlci 
     585          DO jn = 2,jpni-1 
     586           IF(nfipproc(jn,jpnj) .eq. (narea - 1)) THEN 
     587              IF (nfipproc(jn - 1 ,jpnj) .eq. -1) THEN 
     588                 nfsloop = nldi 
     589              ENDIF 
     590              IF (nfipproc(jn + 1,jpnj) .eq. -1) THEN 
     591                 nfeloop = nlei 
     592              ENDIF 
     593           ENDIF 
     594        END DO 
     595 
     596      ENDIF 
     597#if defined key_mpp_mpi 
     598      l_north_nogather = .TRUE. 
     599#endif 
     600   END SUBROUTINE nemo_northcomms 
     601 
    418602 
    419603   !!====================================================================== 

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

    r12414 r13024  
    5454   ! global domain size for AGRIF     !!! * total AGRIF computational domain * 
    5555   INTEGER, PUBLIC            ::   nbug_in_agrif_conv_do_not_remove_or_modify = 1 - 1 
    56    INTEGER, PUBLIC, PARAMETER ::   nbghostcells = 3                             !: number of ghost cells 
    57    INTEGER, PUBLIC            ::   nbcellsx   ! = jpiglo - 2 - 2*nbghostcells   !: number of cells in i-direction 
    58    INTEGER, PUBLIC            ::   nbcellsy   ! = jpjglo - 2 - 2*nbghostcells   !: number of cells in j-direction 
     56   INTEGER, PUBLIC, PARAMETER ::   nbghostcells = 3 !: number of ghost cells: default value 
     57   INTEGER, PUBLIC            ::   nbghostcells_x   !: number of ghost cells in i-direction 
     58   INTEGER, PUBLIC            ::   nbghostcells_y_s   !: number of ghost cells in j-direction at south 
     59   INTEGER, PUBLIC            ::   nbghostcells_y_n   !: number of ghost cells in j-direction at north 
     60   INTEGER, PUBLIC            ::   nbcellsx   ! = jpiglo - 2 - 2*nbghostcells_x   !: number of cells in i-direction 
     61   INTEGER, PUBLIC            ::   nbcellsy   ! = jpjglo - 2 - 2*nbghostcells_y   !: number of cells in j-direction 
    5962 
    6063   ! local domain size                !!! * local computational domain * 
Note: See TracChangeset for help on using the changeset viewer.