Changeset 13024


Ignore:
Timestamp:
2020-06-03T16:26:23+02:00 (4 months 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