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

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

Ignore:
Timestamp:
2018-08-02T15:25:27+02:00 (6 years ago)
Author:
clem
Message:

nesting tools are partly rewritten (mostly for create_coordinates and bathy) to get better functionality. Now you can use the nesting to either define an agrif zoom or a regional domain (for bdy purposes). Also, the nesting tools output a domain_cfg.nc that can be directly used in NEMO4 without the need of DOMAINcfg tool. The option of median average for bathymetry interpolation still does not work properly but it's not new

Location:
utils/tools/NESTING/src
Files:
10 edited

Legend:

Unmodified
Added
Removed
  • utils/tools/NESTING/src/agrif_connect_topo.f90

    r9753 r10025  
    5151    TYPE(Coordinates) :: Grid      
    5252    ! 
    53     ! 
    54     diff = 0 
    55     IF(MOD(rho,2) .EQ. 0) diff = 1 
    56     ! 
    57     indx = 1 + nbghostcellsfine + CEILING(irafx/2.0) + diff       
    58     indy = 1 + nbghostcellsfine + CEILING(irafy/2.0) + diff 
    59     bornex = 1+nbghostcellsfine + CEILING(irafx/2.0) + diff - irafx 
    60     borney = 1+nbghostcellsfine + CEILING(irafy/2.0) + diff - irafy 
    61     bornex2 = nxfin - (nbghostcellsfine) + irafx - CEILING(irafx/2.0)  
    62     borney2 = nyfin - (nbghostcellsfine) + irafy - CEILING(irafy/2.0)  
    63     !     
    64     ALLOCATE(bathy_fin_constant(bornex-FLOOR(irafx/2.0):bornex2+FLOOR(irafx/2.0), & 
    65          borney-FLOOR(irafy/2.0):borney2+FLOOR(irafy/2.0))) 
    66     ! 
    67     DO j = borney,borney2,irafy 
    68  
    69        jpt = jmin + 1 + nbghostcellsfine + (j-indy)/irafy 
    70        IF(j<=1) jpt = jmin + 1 
    71  
    72        DO i = bornex,bornex2,irafx 
    73  
    74           ipt = imin + 1 + nbghostcellsfine + (i-indx)/irafx 
    75           IF(i<=1) ipt = imin + 1 
    76           !        
    77           DO jj = j-FLOOR(irafy/2.0),j+FLOOR(irafy/2.0)-diff 
    78              DO ii = i-FLOOR(irafx/2.0),i+FLOOR(irafx/2.0)-diff 
    79  
    80                 bathy_fin_constant(ii,jj) = coarse_bathy(ipt,jpt) 
    81  
     53    IF( ln_agrif_domain ) THEN 
     54    ! 
     55       diff = 0 
     56       IF(MOD(rho,2) .EQ. 0) diff = 1 
     57       ! 
     58       indx = 1 + nbghostcellsfine + CEILING(irafx/2.0) + diff       
     59       indy = 1 + nbghostcellsfine + CEILING(irafy/2.0) + diff 
     60       bornex = 1+nbghostcellsfine + CEILING(irafx/2.0) + diff - irafx 
     61       borney = 1+nbghostcellsfine + CEILING(irafy/2.0) + diff - irafy 
     62       bornex2 = nxfin - (nbghostcellsfine) + irafx - CEILING(irafx/2.0)  
     63       borney2 = nyfin - (nbghostcellsfine) + irafy - CEILING(irafy/2.0)  
     64       !     
     65       ALLOCATE(bathy_fin_constant(bornex-FLOOR(irafx/2.0):bornex2+FLOOR(irafx/2.0), & 
     66          borney-FLOOR(irafy/2.0):borney2+FLOOR(irafy/2.0))) 
     67       ! 
     68       DO j = borney,borney2,irafy 
     69           
     70          jpt = jmin + 1 + nbghostcellsfine + (j-indy)/irafy 
     71          IF(j<=1) jpt = jmin + 1 
     72           
     73          DO i = bornex,bornex2,irafx 
     74              
     75             ipt = imin + 1 + nbghostcellsfine + (i-indx)/irafx 
     76             IF(i<=1) ipt = imin + 1 
     77             !        
     78             DO jj = j-FLOOR(irafy/2.0),j+FLOOR(irafy/2.0)-diff 
     79                DO ii = i-FLOOR(irafx/2.0),i+FLOOR(irafx/2.0)-diff 
     80                    
     81                   bathy_fin_constant(ii,jj) = coarse_bathy(ipt,jpt) 
     82                    
     83                END DO 
    8284             END DO 
     85              
    8386          END DO 
    84  
    85        END DO 
    86     END DO 
     87       END DO 
     88        
     89    ELSE 
     90 
     91       ALLOCATE(bathy_fin_constant(1:nxfin,1:nyfin)) 
     92 
     93       DO j = 1,nyfin-irafy+1,irafy 
     94          jpt = jmin + FLOOR( (j - 1.) / irafy ) 
     95          ! 
     96          DO i = 1,nxfin-irafx+1,irafx 
     97             ipt = imin + FLOOR( (i - 1.) / irafx ) 
     98             ! 
     99             bathy_fin_constant(i:i+irafx-1,j:j+irafy-1) = coarse_bathy(ipt,jpt) 
     100             ! 
     101          END DO 
     102       END DO 
     103  
     104    ENDIF 
    87105    ! 
    88106    ! 
  • utils/tools/NESTING/src/agrif_create_bathy.f90

    r9753 r10025  
    1313  USE agrif_partial_steps 
    1414  USE agrif_connect_topo 
     15  USE agrif_interpolation 
    1516  ! 
    1617  IMPLICIT NONE 
     
    4041  !       
    4142  CHARACTER(len=80) :: namelistname 
    42   CHARACTER*100 :: Childmeter_file,Childlevel_file,Child_coordinates,child_ps      
    43   LOGICAL,DIMENSION(:,:),POINTER :: masksrc => NULL()   
     43  CHARACTER*100     :: Childmeter_file,Childlevel_file,Child_coordinates,child_ps, Child_domcfg     
    4444  LOGICAL :: identical_grids      
    45   INTEGER,DIMENSION(:,:),ALLOCATABLE ::mask_oce,trouble_points 
    4645  INTEGER :: i,j,num_links,nb,nbadd,status,narg,iargc      
    47   INTEGER,DIMENSION(:),POINTER :: src_add,dst_add => NULL()  
    48   INTEGER :: numlatfine,numlonfine,numlat,numlon,pos,pos2 
    49   REAL*8,DIMENSION(:,:),POINTER :: matrix,interpdata => NULL()      
    50   REAL*8, DIMENSION(:,:),POINTER :: bathy_fin_constant => NULL()   
    51   REAL*8,DIMENSION(:,:),ALLOCATABLE :: bathy_test,vardep,glamhr,gphihr 
    52   REAL*8,DIMENSION(:),ALLOCATABLE :: vardep1d 
    53   REAL*8, DIMENSION(:,:),POINTER :: gdepw_ps_interp => NULL()  
    54   REAL*8, DIMENSION(:,:),POINTER :: save_gdepw,rx,ry,maskedtopo 
     46  INTEGER :: jpj,jpi,pos,pos2 
     47  LOGICAL,DIMENSION(:,:),POINTER     :: masksrc => NULL()   
     48  INTEGER,DIMENSION(:,:),ALLOCATABLE :: mask_oce,trouble_points 
     49  INTEGER,DIMENSION(:)  ,POINTER     :: src_add,dst_add => NULL()  
     50  REAL*8, DIMENSION(:,:),POINTER     :: matrix,interpdata => NULL()      
     51  REAL*8, DIMENSION(:,:),POINTER     :: bathy_fin_constant => NULL()   
     52  REAL*8, DIMENSION(:,:),ALLOCATABLE :: bathy_test,vardep,glamhr,gphihr 
     53  REAL*8, DIMENSION(:)  ,ALLOCATABLE :: vardep1d 
     54  REAL*8, DIMENSION(:,:),POINTER     :: gdepw_ps_interp => NULL()  
     55  REAL*8, DIMENSION(:,:),POINTER     :: save_gdepw,rx,ry,maskedtopo 
    5556  REAL*8  :: Cell_lonmin,Cell_lonmax,Cell_latmin,Cell_latmax,wghts 
    56   LOGICAL :: Pacifique=.FALSE. 
     57  LOGICAL :: Pacifique = .FALSE. 
    5758  INTEGER :: boundary,xpos,ypos,iimin,iimax,jjmax,jjmin 
    58   INTEGER :: nbloops,nxhr,nyhr,ji,jj,nbiter,nbloopmax 
     59  INTEGER :: nxhr,nyhr,ji,jj,nbiter 
    5960  INTEGER :: ipt,jpt,iloc,jloc 
    6061  INTEGER, DIMENSION(2) :: i_min,i_max,j_min,j_max 
     
    7071  ! 
    7172  ! read input file (namelist.input) 
    72   ! 
    7373  CALL read_namelist(namelistname) 
    7474  !       
    7575  ! define names of child grid files 
     76  CALL set_child_name(parent_coordinate_file,Child_coordinates)  
     77  IF( TRIM(parent_meshmask_file) .NE. '/NULL' )    CALL set_child_name(parent_meshmask_file,Childlevel_file)             
     78  IF( TRIM(parent_meshmask_file) .NE. '/NULL' )    CALL set_child_name(parent_domcfg_output,Child_domcfg)  
    7679  ! 
    77   CALL set_child_name(parent_coordinate_file,Child_coordinates)  
    78   IF( TRIM(parent_meshmask_file) .NE. '/NULL' ) & 
    79        CALL set_child_name(parent_meshmask_file,Childlevel_file)             
    80   ! 
    81   ! 
    82   ! 
    83   ! 
    84   ! 
    85   !------------------------------------------------------------------ 
    86   ! ****First option : no new topo file / no partial steps 
    87   ! interpolate levels directly from parent file 
    88   !------------------------------------------------------------------ 
    89   ! 
    90   ! 
    91   ! 
    92   ! 
    93   ! 
    94   ! 
    95   IF(.NOT.new_topo .AND. .NOT.partial_steps) THEN       
     80  !                                                   ! ------------------------------------------------------------------ 
     81  IF( .NOT.new_topo .AND. .NOT.partial_steps ) THEN   ! **** First option : no new topo file & no partial steps 
     82     !                                                !                 ==> interpolate levels directly from parent file 
     83     !                                                ! ------------------------------------------------------------------ 
     84     WRITE(*,*) '*** First option: no new topo file & no partial steps' 
    9685     !       
    97      WRITE(*,*) 'First option' 
    98      !read coarse grid bathymetry and coordinates file 
    99      !            
    100      WRITE(*,*) 'No new topo file ...' 
     86     ! read coarse grid bathymetry and coordinates 
    10187     status = Read_Coordinates(TRIM(parent_coordinate_file),G0)     
    10288     status = Read_bathy_level(TRIM(parent_meshmask_file),G0) 
    10389     !            
    104      IF( imax > SIZE(G0%glamt,1) .OR. jmax > SIZE(G0%glamt,2) .OR. & 
    105           imax <= imin .OR. jmax <= jmin ) THEN                    
     90     ! read fine grid coordinates 
     91     status = Read_Coordinates(TRIM(Child_coordinates),G1,pacifique) 
     92     ! 
     93     ! stop if error in size 
     94     IF( imax > SIZE(G0%glamt,1) .OR. jmax > SIZE(G0%glamt,2) .OR. imax <= imin .OR. jmax <= jmin ) THEN                    
    10695        WRITE(*,*) 'ERROR ***** bad child grid definition ...' 
    10796        WRITE(*,*) 'please check imin,jmin,imax,jmax,jpizoom,jpjzoom values'        
    108         WRITE(*,*) ' ' 
    10997        STOP 
    11098     ENDIF 
    111      ! 
    112      !read fine grid coordinates file 
    113      !       
    114      status = Read_Coordinates(TRIM(Child_coordinates),G1,pacifique) 
    115      !   
    11699     IF( SIZE(G1%nav_lon,1) .NE. nxfin .OR. SIZE(G1%nav_lon,2) .NE. nyfin ) THEN 
    117         ! 
    118100        WRITE(*,*) 'ERROR ***** bad child coordinates file ...' 
    119         WRITE(*,*) ' ' 
    120         WRITE(*,*) 'please check that child coordinates file ' 
    121         WRITE(*,*) 'has been created with the same namelist ' 
    122         WRITE(*,*) ' ' 
     101        WRITE(*,*) 'please check that child coordinates file has been created with the same namelist' 
    123102        STOP 
    124         ! 
    125      ENDIF 
    126      ! 
    127      numlat =  SIZE(G0%nav_lat,2) 
    128      numlon =  SIZE(G0%nav_lat,1)     
    129      numlatfine =  SIZE(G1%nav_lat,2) 
    130      numlonfine =  SIZE(G1%nav_lat,1)   
     103     ENDIF 
     104     ! 
     105     jpj = SIZE(G0%nav_lat,2) 
     106     jpi = SIZE(G0%nav_lat,1)     
    131107     !            
    132      ALLOCATE(masksrc(numlon,numlat)) 
    133      ! 
    134108     ! create logical array masksrc 
    135      ! 
    136      WHERE(G0%bathy_level.LE.0)  
    137         masksrc = .FALSE. 
    138      ELSEWHERE 
    139         masksrc = .TRUE. 
     109     ALLOCATE(masksrc(jpi,jpj)) 
     110     WHERE(G0%bathy_level.LE.0)   ;   masksrc = .FALSE.   ; 
     111     ELSEWHERE                    ;   masksrc = .TRUE.    ; 
    140112     END WHERE 
    141113 
     114     ! change longitudes (from -180:180 to 0:360) 
    142115     IF ( Pacifique ) THEN 
    143         WHERE(G0%nav_lon < 0.001)  
    144            G0%nav_lon = G0%nav_lon + 360. 
    145         END WHERE 
    146      ENDIF 
    147  
    148  
    149      !-----------------           
    150      ! compute remapping matrix thanks to SCRIP package 
    151      ! 
    152      ! remapping process 
    153      !               
     116        WHERE(G0%nav_lon < 0.001)   G0%nav_lon = G0%nav_lon + 360. 
     117     ENDIF 
     118 
     119     ! From levels to meters 
    154120     ALLOCATE(G1%bathy_meter(nxfin,nyfin)) 
    155121     CALL levels_to_meter(G0) 
    156122     !              
    157      !             Call levels_to_meter(G1) 
    158      !              
    159      CALL get_remap_matrix(G0%nav_lat,G1%nav_lat,   & 
    160           G0%nav_lon,G1%nav_lon,   & 
    161           masksrc,matrix,src_add,dst_add) 
    162      CALL make_remap(G0%bathy_meter,G1%bathy_meter,nxfin,nyfin, & 
    163           matrix,src_add,dst_add)   
    164      !                                   
    165      !             
     123     ! compute remapping matrix thanks to SCRIP package 
     124     CALL get_remap_matrix(G0%nav_lat,G1%nav_lat,G0%nav_lon,G1%nav_lon,masksrc,matrix,src_add,dst_add) 
     125     CALL make_remap(G0%bathy_meter,G1%bathy_meter,nxfin,nyfin,matrix,src_add,dst_add)   
    166126     DEALLOCATE(masksrc) 
    167      !-----------------                                       
    168127     !       
    169      !  
    170128     ! compute constant bathymetry for Parent-Child bathymetry connection 
    171      !               
    172129     CALL init_constant_bathy(G0%bathy_meter,bathy_fin_constant) 
    173130     ! 
    174      boundary = npt_copy*irafx + nbghostcellsfine + 1  
    175      ! 
    176      ! connection carried out by copying parent grid values for the fine points 
    177      ! corresponding to 3 coarse grid cells at the boundaries 
    178      !                  
     131     ! replace child bathymetry by parent bathymetry at the boundaries 
     132     IF( ln_agrif_domain ) THEN 
     133        boundary = npt_copy*irafx + nbghostcellsfine + 1  
     134     ELSE 
     135        boundary = npt_copy*irafx 
     136     ENDIF 
    179137     G1%bathy_meter(1:boundary,:) = bathy_fin_constant(1:boundary,:) 
    180138     G1%bathy_meter(:,1:boundary) = bathy_fin_constant(:,1:boundary) 
    181139     G1%bathy_meter(nxfin-boundary+1:nxfin,:) = bathy_fin_constant(nxfin-boundary+1:nxfin,:) 
    182140     G1%bathy_meter(:,nyfin-boundary+1:nyfin) = bathy_fin_constant(:,nyfin-boundary+1:nyfin) 
     141     DEALLOCATE(bathy_fin_constant) 
    183142     !                   
    184      CALL smooth_topo(G1%bathy_meter(boundary:nxfin-boundary+1,boundary:nyfin-boundary+1),nbiter) 
    185      !              
     143     ! bathymetry smoothing (everywhere except at the boundaries)  
     144     IF( smoothing ) THEN 
     145        IF( ln_agrif_domain ) THEN 
     146           CALL smooth_topo(G1%bathy_meter(boundary:nxfin-boundary+1,boundary:nyfin-boundary+1),nbiter) 
     147        ELSE 
     148           CALL smooth_topo(G1%bathy_meter(boundary+1:nxfin-boundary,boundary+1:nyfin-boundary),nbiter) 
     149        ENDIF 
     150     ENDIF 
     151     ! 
     152     ! From meters to levels 
    186153     CALL meter_to_levels(G1) 
    187      !              
    188      DEALLOCATE(bathy_fin_constant) 
    189      ! 
    190      ! 
    191      !------------------------------------------------------------------ 
    192      ! ****Second option : new topo file or/and partial steps      
    193      !------------------------------------------------------------------  
    194      ! 
    195      ! 
    196      ! 
    197      ! 
    198      ! 
    199      ! 
    200      ! 
    201      ! 
    202   ELSE 
    203      ! 
    204      WRITE(*,*) 'Second option : partial steps' 
     154     G1%bathy_level=NINT(G1%bathy_level) 
     155     !        
     156     ! Check errors in bathy 
     157     DO j=1,nyfin 
     158        DO i=1,nxfin 
     159           IF (G1%bathy_level(i,j).LT.0.) THEN 
     160              PRINT *,'error in ',i,j,G1%bathy_level(i,j) 
     161              STOP 
     162           ENDIF 
     163        ENDDO 
     164     ENDDO 
     165     WHERE ((G1%bathy_level.LT.3.).AND.(G1%bathy_level.GT.0.))   G1%bathy_level=3 
     166     ! 
     167     ! remove closed seas 
     168     IF (removeclosedseas) THEN 
     169        ALLOCATE(bathy_test(nxfin,nyfin)) 
     170 
     171        bathy_test=0. 
     172        WHERE (G1%bathy_level(1,:)     .GT.0.)   bathy_test(1,:)=1 
     173        WHERE (G1%bathy_level(nxfin,:) .GT.0.)   bathy_test(nxfin,:)=1 
     174        WHERE (G1%bathy_level(:,1)     .GT.0.)   bathy_test(:,1)=1 
     175        WHERE (G1%bathy_level(:,nyfin) .GT.0.)   bathy_test(:,nyfin)=1 
     176 
     177        nbadd = 1 
     178        DO WHILE (nbadd.NE.0) 
     179           nbadd = 0 
     180           DO j=2,nyfin-1 
     181              DO i=2,nxfin-1 
     182                 IF (G1%bathy_level(i,j).GT.0.) THEN 
     183                    IF (MAX(bathy_test(i,j+1),bathy_test(i,j-1),bathy_test(i-1,j),bathy_test(i+1,j)).EQ.1) THEN 
     184                       IF (bathy_test(i,j).NE.1.) nbadd = nbadd + 1 
     185                       bathy_test(i,j)=1. 
     186                    ENDIF 
     187                 ENDIF 
     188              ENDDO 
     189           ENDDO 
     190        ENDDO 
     191 
     192        WHERE (bathy_test.EQ.0.)   G1%bathy_level = 0. 
     193 
     194        DEALLOCATE(bathy_test)            
     195     ENDIF 
     196     ! 
     197     ! store interpolation result in output file 
     198     CALL levels_to_meter(G1)   ! From levels to meters 
     199     status = Write_Bathy_level(TRIM(Childlevel_file),G1) 
     200     status = write_domcfg(TRIM(Child_domcfg),G1) 
     201 
     202     WRITE(*,*) '****** Bathymetry successfully created if no new topo ******' 
     203     STOP 
     204     ! 
     205     !                                                ! ----------------------------------------------------- 
     206  ELSE                                                ! **** Second option : new topo file or partial steps      
     207     !                                                ! -----------------------------------------------------  
     208     WRITE(*,*) '***Second option : new topo or partial steps' 
     209 
    205210     ! read fine and coarse grids coordinates file 
    206      !        
    207211     status = Read_Coordinates(TRIM(parent_coordinate_file),G0) 
    208212     status = Read_Coordinates(TRIM(Child_coordinates),G1,Pacifique) 
    209      !                         
    210      IF( imax > SIZE(G0%nav_lon,1) .OR. jmax > SIZE(G0%nav_lon,2) .OR. & 
    211          imax <= imin .OR. jmax <= jmin ) THEN                     
     213     ! 
     214     ! check error in size 
     215     IF( imax > SIZE(G0%nav_lon,1) .OR. jmax > SIZE(G0%nav_lon,2) .OR. imax <= imin .OR. jmax <= jmin ) THEN                     
    212216        WRITE(*,*) 'ERROR ***** bad child grid definition ...' 
    213217        WRITE(*,*) 'please check imin,jmin,imax,jmax,jpizoom,jpjzoom values'        
    214         WRITE(*,*) ' ' 
    215218        STOP 
    216219     ENDIF 
    217      ! 
    218  
    219      !   
    220220     IF( SIZE(G1%nav_lon,1) .NE. nxfin .OR. SIZE(G1%nav_lon,2) .NE. nyfin ) THEN 
    221         ! 
    222221        WRITE(*,*) 'ERROR ***** bad child coordinates file ...' 
    223         WRITE(*,*) ' ' 
    224         WRITE(*,*) 'please check that child coordinates file ' 
    225         WRITE(*,*) 'has been created with the same namelist ' 
    226         WRITE(*,*) ' ' 
     222        WRITE(*,*) 'please check that child coordinates file has been created with the same namelist' 
    227223        STOP 
    228         ! 
    229      ENDIF 
    230      !       
    231      ! read coarse grid bathymetry file 
    232      !--- 
    233      IF( new_topo ) THEN 
    234         WRITE(*,*) 'use new topo file : ',TRIM(elevation_database) 
     224     ENDIF 
     225     ! 
     226     ! === From here on: G0 is the grid associated with the new topography (like gebco or etopo) === 
     227     ! 
     228     ! read bathymetry data set => G0%bathy_meter 
     229     IF( new_topo ) THEN                                       ! read G0%bathy_meter (on a reduced grid) and G1 coordinates 
    235230        DEALLOCATE( G0%nav_lon, G0%nav_lat ) 
    236231        status = Read_bathy_meter(TRIM(elevation_database),G0,G1,Pacifique) 
    237      ELSE 
    238         WRITE(*,*) 'no new topo file' 
    239         status = Read_Bathymeter(TRIM(parent_bathy_meter),G0) 
     232     ELSE                                                      ! read G0%bathy_meter (on the global grid) 
     233        status = Read_Bathymeter(TRIM(parent_bathy_meter),G0)    
    240234        IF(Pacifique) THEN 
    241            WHERE(G0%nav_lon < 0.0001)  
    242               G0%nav_lon = G0%nav_lon + 360. 
    243            END WHERE 
    244         ENDIF 
    245      ENDIF 
    246      !---             
    247      numlatfine =  SIZE(G1%nav_lat,2) 
    248      numlonfine =  SIZE(G1%nav_lat,1)  
    249      !   
     235           WHERE(G0%nav_lon < 0.0001)   G0%nav_lon = G0%nav_lon + 360. 
     236        ENDIF 
     237     ENDIF 
    250238     !                
     239     ! what type of interpolation for bathymetry 
     240     IF( type_bathy_interp == 0 ) THEN 
     241        WRITE(*,*) 'Interpolation of high resolution bathymetry on child grid: Arithmetic average ...' 
     242     ELSE IF( type_bathy_interp == 1 ) THEN 
     243        WRITE(*,*) 'Interpolation of high resolution bathymetry on child grid: Median average ...' 
     244     ELSE IF( type_bathy_interp == 2 ) THEN      
     245        WRITE(*,*) 'Interpolation of high resolution bathymetry on child grid: Bilinear interpolation ...' 
     246     ELSE      
     247        WRITE(*,*) 'bad value for type_bathy_interp variable ( must be 0, 1 or 2 )' 
     248        STOP  
     249     ENDIF 
     250     ! 
     251     ! 1st allocation of child grid bathy 
    251252     ALLOCATE(G1%bathy_meter(nxfin,nyfin)) 
    252253     G1%bathy_meter(:,:)=0.                        
    253  
    254      WRITE(*,*) 'Interpolation of high resolution bathymetry on child grid' 
    255  
    256      IF( type_bathy_interp == 0 ) THEN 
    257         WRITE(*,*) 'Arithmetic average ...' 
    258      ELSE IF( type_bathy_interp == 1 ) THEN 
    259         WRITE(*,*) 'Median average ...' 
    260      ELSE IF( type_bathy_interp == 2 ) THEN      
    261         WRITE(*,*) 'Bilinear interpolation ...' 
    262      ELSE      
    263         WRITE(*,*) 'bad value for type_bathy_interp variable ( must be 0,1 or 2 )' 
    264         STOP  
    265      ENDIF 
    266      ! 
    267      !************************************ 
    268      !MEDIAN AVERAGE or ARITHMETIC AVERAGE 
    269      !************************************ 
    270      ! 
    271      IF( type_bathy_interp == 0 .OR. type_bathy_interp == 1 ) THEN  
    272         ! 
     254     ! 
     255     ! --------------------------------------------------------------------------------- 
     256     ! ===                 Bathymetry of the fine grid (step1)                       === 
     257     ! --------------------------------------------------------------------------------- 
     258     ! ==> It gives G1%bathy_meter from G0%bathy_meter 
     259     ! --------------------------------------------------------------------------------- 
     260     !                                                               ! -----------------------------  
     261     IF( type_bathy_interp == 0 .OR. type_bathy_interp == 1 ) THEN   ! arithmetic or median averages 
     262        !                                                            ! -----------------------------  
    273263        ALLOCATE(trouble_points(nxfin,nyfin)) 
    274         trouble_points = 0 
    275         ! 
    276         !  POINT DETECTION 
    277         ! 
     264        trouble_points(:,:) = 0 
    278265        !                        
    279         DO jj = 2,numlatfine 
    280            DO ji = 2,numlonfine 
     266        DO jj = 2,nyfin 
     267           DO ji = 2,nxfin 
    281268              !        
    282               ! FINE GRID CELLS DEFINITION                
    283               ! 
     269              ! fine grid cell extension                
    284270              Cell_lonmin = MIN(G1%glamf(ji-1,jj-1),G1%glamf(ji,jj-1),G1%glamf(ji,jj),G1%glamf(ji-1,jj)) 
    285271              Cell_lonmax = MAX(G1%glamf(ji-1,jj-1),G1%glamf(ji,jj-1),G1%glamf(ji,jj),G1%glamf(ji-1,jj)) 
     
    287273              Cell_latmax = MAX(G1%gphif(ji-1,jj-1),G1%gphif(ji,jj-1),G1%gphif(ji,jj),G1%gphif(ji-1,jj))  
    288274              !                
    289               ! SEARCH FOR ALL POINTS CONTAINED IN THIS CELL 
    290               ! 
     275              ! look for points in G0 (bathy dataset) contained in the fine grid cells   
    291276              iimin = 1 
    292277              DO WHILE( G0%nav_lon(iimin,1) < Cell_lonmin )  
     
    302287              DO WHILE( G0%nav_lon(iimax,1) <= Cell_lonmax )  
    303288                 iimax = iimax + 1 
    304               iimax = MIN( iimax,SIZE(G0%bathy_meter,1)) 
     289                 iimax = MIN( iimax,SIZE(G0%bathy_meter,1)) 
    305290              ENDDO 
    306291              !                                
     
    308293              DO WHILE( G0%nav_lat(iimax,jjmax) <= Cell_latmax )  
    309294                 jjmax = jjmax + 1 
    310               jjmax = MIN( jjmax,SIZE(G0%bathy_meter,2)) 
    311  
     295                 jjmax = MIN( jjmax,SIZE(G0%bathy_meter,2)) 
    312296              ENDDO 
    313297              ! 
    314               iimax = iimax-1 
    315               jjmax = jjmax-1 
     298              IF( ln_agrif_domain ) THEN 
     299                 iimax = iimax-1 
     300                 jjmax = jjmax-1 
     301              ELSE 
     302                 iimax = MAX(iimin,iimax-1) 
     303                 jjmax = MAX(jjmin,jjmax-1) 
     304              ENDIF 
    316305              !                
    317306              iimin = MAX( iimin,1 ) 
     
    324313 
    325314              IF( nxhr == 0 .OR. nyhr == 0 ) THEN 
     315                 ! 
    326316                 trouble_points(ji,jj) = 1 
     317                 ! 
    327318              ELSE 
    328  
    329                  ALLOCATE( vardep(nxhr,nyhr) ) 
    330                  ALLOCATE( mask_oce(nxhr,nyhr) ) 
    331                  mask_oce = 0        
    332  
     319                 ! 
     320                 ALLOCATE( vardep(nxhr,nyhr), mask_oce(nxhr,nyhr) ) 
    333321                 vardep(:,:) = G0%bathy_meter(iimin:iimax,jjmin:jjmax) 
    334  
    335                  WHERE( vardep(:,:) .GT. 0. )  mask_oce = 1 
    336  
    337 !                 IF( SUM(mask_oce) == 0 ) THEN 
    338                  IF( SUM(mask_oce) < 0.5*(nxhr*nyhr) ) THEN 
     322                 ! 
     323                 WHERE( vardep(:,:) .GT. 0. )   ;   mask_oce = 1 ; 
     324                 ELSEWHERE                      ;   mask_oce = 0 ; 
     325                 ENDWHERE 
     326                 ! 
     327                 IF( SUM(mask_oce) < 0.5*(nxhr*nyhr) ) THEN ! if more than half of the points are on land then bathy fine = 0 
    339328                    G1%bathy_meter(ji,jj) = 0. 
    340329                 ELSE 
    341                     IF( type_bathy_interp == 0 ) THEN 
    342                        ! Arithmetic average                    
    343                        G1%bathy_meter(ji,jj) = SUM (vardep(:,:)*mask_oce(:,:))/SUM(mask_oce) 
    344                     ELSE 
    345                        ! Median average           
    346                        ! 
    347                        vardep(:,:) = vardep(:,:)*mask_oce(:,:) - 100000*(1-mask_oce(:,:)) 
     330                    IF( type_bathy_interp == 0 ) THEN     ! Arithmetic average 
     331                       G1%bathy_meter(ji,jj) = SUM( vardep(:,:) * mask_oce(:,:) ) / SUM( mask_oce(:,:) ) 
     332                    ELSE                                  ! Median average 
     333                       vardep(:,:) = vardep(:,:) * mask_oce(:,:) - 100000 * ( 1 - mask_oce(:,:) ) 
    348334                       ALLOCATE(vardep1d(nxhr*nyhr)) 
    349335                       vardep1d = RESHAPE(vardep,(/ nxhr*nyhr /) ) 
     
    351337                       ! 
    352338                       ! Calculate median 
    353                        ! 
    354339                       IF (MOD(SUM(mask_oce),2) .NE. 0) THEN 
    355                           G1%bathy_meter(ji,jj) = vardep1d( SUM(mask_oce)/2 + 1) 
     340                          G1%bathy_meter(ji,jj) = vardep1d( SUM(mask_oce)/2 + 1 ) 
    356341                       ELSE 
    357                           G1%bathy_meter(ji,jj) = ( vardep1d(SUM(mask_oce)/2) + vardep1d(SUM(mask_oce)/2+1) )/2.0 
     342                          G1%bathy_meter(ji,jj) = 0.5 * ( vardep1d(SUM(mask_oce)/2) + vardep1d(SUM(mask_oce)/2+1) ) 
    358343                       END IF 
    359344                       DEALLOCATE(vardep1d)    
     
    361346                 ENDIF 
    362347                 DEALLOCATE (mask_oce,vardep) 
    363  
     348                 ! 
    364349              ENDIF 
    365350           ENDDO 
     
    367352 
    368353        IF( SUM( trouble_points ) > 0 ) THEN 
    369            PRINT*,'too much empty cells, proceed to bilinear interpolation !!' 
     354           PRINT*,'too much empty cells, proceed to bilinear interpolation' 
    370355           type_bathy_interp = 2 
    371356        ENDIF 
    372  
    373      ENDIF 
    374  
    375      ! 
    376      ! create logical array masksrc 
    377      ! 
    378      IF( type_bathy_interp == 2) THEN  
    379         ! 
    380  
    381         !             
     357            
     358     ENDIF 
     359     !                                                       ! -----------------------------  
     360     IF( type_bathy_interp == 2) THEN                        ! Bilinear interpolation 
     361        !                                                    ! -----------------------------  
    382362        identical_grids = .FALSE. 
    383363 
    384         IF( SIZE(G0%nav_lat,1) == SIZE(G1%nav_lat,1)  .AND.   & 
    385              SIZE(G0%nav_lat,2) == SIZE(G1%nav_lat,2)  .AND.   & 
    386              SIZE(G0%nav_lon,1) == SIZE(G1%nav_lon,1)  .AND.   & 
    387              SIZE(G0%nav_lon,2) == SIZE(G1%nav_lon,2)   ) THEN 
     364        IF( SIZE(G0%nav_lat,1) == SIZE(G1%nav_lat,1) .AND. SIZE(G0%nav_lat,2) == SIZE(G1%nav_lat,2)  .AND.   & 
     365            SIZE(G0%nav_lon,1) == SIZE(G1%nav_lon,1) .AND. SIZE(G0%nav_lon,2) == SIZE(G1%nav_lon,2) ) THEN 
    388366           IF( MAXVAL( ABS(G0%nav_lat(:,:)- G1%nav_lat(:,:)) ) < 0.0001 .AND.   & 
    389                 MAXVAL( ABS(G0%nav_lon(:,:)- G1%nav_lon(:,:)) ) < 0.0001 ) THEN 
     367               MAXVAL( ABS(G0%nav_lon(:,:)- G1%nav_lon(:,:)) ) < 0.0001 ) THEN 
     368              PRINT*,'same grid between ',elevation_database,' and child domain'     
    390369              G1%bathy_meter = G0%bathy_meter  
    391               PRINT*,'same grid between ',elevation_database,' and child domain'     
    392370              identical_grids = .TRUE.                           
    393371           ENDIF 
    394372        ENDIF 
    395  
    396373 
    397374        IF( .NOT. identical_grids ) THEN  
     
    400377           ALLOCATE(bathy_test(nxfin,nyfin)) 
    401378           ! 
    402            !                    Where(G0%bathy_meter.le.0.00001)  
    403            !                   masksrc = .false. 
    404            !               ElseWhere 
    405            ! 
    406            masksrc = .TRUE. 
    407            ! 
    408            !               End where                        
     379           !Where(G0%bathy_meter.le.0.00001)  
     380           !  masksrc = .false. 
     381           !ElseWhere 
     382              masksrc = .TRUE. 
     383           !End where                        
    409384           !             
    410385           ! compute remapping matrix thanks to SCRIP package             
    411            !                                   
    412            CALL get_remap_matrix(G0%nav_lat,G1%nav_lat,   & 
    413                 G0%nav_lon,G1%nav_lon,   & 
    414                 masksrc,matrix,src_add,dst_add) 
    415            CALL make_remap(G0%bathy_meter,bathy_test,nxfin,nyfin, & 
    416                 matrix,src_add,dst_add)   
     386           CALL get_remap_matrix(G0%nav_lat,G1%nav_lat,G0%nav_lon,G1%nav_lon,masksrc,matrix,src_add,dst_add) 
     387           CALL make_remap(G0%bathy_meter,bathy_test,nxfin,nyfin,matrix,src_add,dst_add)   
    417388           !                                   
    418389           G1%bathy_meter = bathy_test                
     
    424395        !             
    425396     ENDIF 
    426      ! 
    427      ! 
    428      ! 
    429      !------------------------------------------------------------------------------------------ 
    430      ! ! ****Third  option : Partial Steps 
    431      ! The code includes the  
    432      ! option to include partial cells at the bottom  
    433      ! in order to better resolve topographic variations 
    434      !------------------------------------------------------------------------------------------ 
    435      ! 
    436      ! At this step bathymetry in meters has already been interpolated on fine grid 
    437      ! 
    438      !                    
    439      IF( partial_steps ) THEN                 
    440         !                   
     397     ! --- 
     398     ! At this stage bathymetry in meters has already been interpolated on fine grid 
     399     !                    => G1%bathy_meter(nxfin,nyfin) 
     400     ! 
     401     ! Also G0 was the grid from the new bathymetry data set (etopo, gebco...) and not the coarse grid 
     402     ! --- 
     403     ! 
     404     ! --------------------------------------------------------------------------------- 
     405     ! ===                 Bathymetry of the fine grid (step2)                       === 
     406     ! --------------------------------------------------------------------------------- 
     407     ! ==> It gives an update of G1%bathy_meter and G1%bathy_level 
     408     ! --------------------------------------------------------------------------------- 
     409     ! From here on: G0 is the coarse grid 
     410     ! 
     411     ! Coarse grid bathymetry : G0%bathy_meter (on the global grid) 
     412     IF( partial_steps ) THEN                      
    441413        status = Read_Bathymeter(TRIM(parent_bathy_meter),G0) 
    442         DEALLOCATE(G0%nav_lat,G0%nav_lon) 
    443         status = Read_coordinates(TRIM(parent_coordinate_file),G0) 
    444         !------------------------                   
    445         IF (.NOT.ASSOCIATED(G0%gdepw_ps)) & 
    446              ALLOCATE(G0%gdepw_ps(SIZE(G0%bathy_meter,1),SIZE(G0%bathy_meter,2))) 
    447         IF (.NOT.ASSOCIATED(G1%gdepw_ps)) & 
    448              ALLOCATE(G1%gdepw_ps(SIZE(G1%bathy_meter,1),SIZE(G1%bathy_meter,2)))                   
    449         IF (.NOT.ASSOCIATED(gdepw_ps_interp)) & 
    450              ALLOCATE(gdepw_ps_interp(SIZE(G1%bathy_meter,1),SIZE(G1%bathy_meter,2))) 
    451         ! 
    452         !                        
    453         WRITE(*,*) 'Coarse grid : ' 
    454         CALL get_partial_steps(G0)  
    455         WRITE(*,*) ' ' 
    456         WRITE(*,*) 'Fine grid : ' 
    457         CALL get_partial_steps(G1)                 ! compute gdepw_ps for G1 
    458         CALL bathymetry_control(G0%Bathy_level)    !     
    459         CALL Check_interp(G0,gdepw_ps_interp)      ! interpolation in connection zone (3 coarse grid cells) 
    460         ! 
    461         boundary = npt_copy*irafx + nbghostcellsfine + 1                      
    462 ! J chanut: exclude matching if no open boundaries 
    463         IF (.NOT.ASSOCIATED(G1%wgt)) & 
    464              ALLOCATE(G1%wgt(SIZE(G1%bathy_meter,1),SIZE(G1%bathy_meter,2))) 
    465         G1%wgt(:,:) = 0. 
    466         IF ((.NOT.ASSOCIATED(G0%wgt)).AND.bathy_update) THEN  
    467              ALLOCATE(G0%wgt(SIZE(G0%nav_lat,1),SIZE(G0%nav_lat,2))) 
    468              G0%wgt(:,:) = 0. 
    469         ENDIF 
    470  
    471         DO jj=1,nyfin 
    472           ! West 
    473           IF (gdepw_ps_interp(nbghostcellsfine+1,jj)>0.) THEN 
    474              G1%gdepw_ps(1:boundary,jj) = gdepw_ps_interp(1:boundary,jj)  
    475              G1%wgt(1:boundary,jj) = 1. 
    476           ELSE 
    477              G1%gdepw_ps(1:nbghostcellsfine+1,jj)=0.  
    478           ENDIF  
    479           ! East 
    480           IF (gdepw_ps_interp(nxfin-nbghostcellsfine,jj)>0.) THEN 
    481              G1%gdepw_ps(nxfin-boundary+1:nxfin,jj)=gdepw_ps_interp(nxfin-boundary+1:nxfin,jj) 
    482              G1%wgt(nxfin-boundary+1:nxfin,jj) = 1. 
    483           ELSE 
    484              G1%gdepw_ps(nxfin-nbghostcellsfine:nxfin,jj) = 0. 
    485           ENDIF 
     414     ELSE 
     415        status = Read_Bathymeter(TRIM(parent_meshmask_file),G0) 
     416     ENDIF 
     417      
     418     ! Coarse grid coordinatees : G0 coordinates 
     419     DEALLOCATE(G0%nav_lat,G0%nav_lon) 
     420     status = Read_coordinates(TRIM(parent_coordinate_file),G0) 
     421 
     422     ! allocate temporary arrays                   
     423     IF (.NOT.ASSOCIATED(G0%gdepw_ps))       ALLOCATE(G0%gdepw_ps    (SIZE(G0%bathy_meter,1),SIZE(G0%bathy_meter,2))) 
     424     IF (.NOT.ASSOCIATED(G1%gdepw_ps))       ALLOCATE(G1%gdepw_ps    (SIZE(G1%bathy_meter,1),SIZE(G1%bathy_meter,2)))                   
     425     IF (.NOT.ASSOCIATED(gdepw_ps_interp))   ALLOCATE(gdepw_ps_interp(SIZE(G1%bathy_meter,1),SIZE(G1%bathy_meter,2))) 
     426     !                        
     427     IF( ln_agrif_domain ) THEN 
     428        boundary = npt_copy*irafx + nbghostcellsfine + 1 
     429     ELSE 
     430        boundary = npt_copy*irafx 
     431     ENDIF 
     432     ! 
     433     ! compute G0%gdepw_ps and G1%gdepw_ps 
     434     CALL get_partial_steps(G0)  
     435     CALL get_partial_steps(G1) 
     436     CALL bathymetry_control(G0%Bathy_level) 
     437 
     438     ! --------------------------------------- 
     439     ! Bathymetry at the boundaries (npt_copy)                       
     440     ! --------------------------------------- 
     441     ! 1st step: interpolate coarse bathy on the fine grid (using partial steps or not) 
     442     IF( partial_steps .AND. ln_agrif_domain ) THEN                    
     443        CALL Check_interp(G0,gdepw_ps_interp) 
     444     ELSE 
     445        gdepw_ps_interp = 0. * G1%gdepw_ps 
     446        !!CALL agrif_interp(G0%gdepw_ps,gdepw_ps_interp,'T') 
     447        CALL init_constant_bathy(G0%gdepw_ps,gdepw_ps_interp) 
     448     ENDIF 
     449 
     450     IF (.NOT.ASSOCIATED(G1%wgt))   ALLOCATE(G1%wgt(SIZE(G1%bathy_meter,1),SIZE(G1%bathy_meter,2))) 
     451     G1%wgt(:,:) = 0. 
     452     IF ((.NOT.ASSOCIATED(G0%wgt)).AND.bathy_update) THEN  
     453        ALLOCATE(G0%wgt(SIZE(G0%nav_lat,1),SIZE(G0%nav_lat,2))) 
     454        G0%wgt(:,:) = 0. 
     455     ENDIF 
     456     ! 
     457     ! 2nd step: copy parent bathymetry at the boundaries 
     458     DO jj=1,nyfin   ! West and East 
     459        IF ( gdepw_ps_interp(nbghostcellsfine+1,jj) > 0. ) THEN 
     460           G1%gdepw_ps(1:boundary,jj) = gdepw_ps_interp(1:boundary,jj)  
     461           G1%wgt(1:boundary,jj) = 1. 
     462        ELSE 
     463           G1%gdepw_ps(1:nbghostcellsfine+1,jj)=0.  
     464        ENDIF 
     465        ! 
     466        IF ( gdepw_ps_interp(nxfin-nbghostcellsfine,jj) > 0.) THEN 
     467           G1%gdepw_ps(nxfin-boundary+1:nxfin,jj)=gdepw_ps_interp(nxfin-boundary+1:nxfin,jj) 
     468           G1%wgt(nxfin-boundary+1:nxfin,jj) = 1. 
     469        ELSE 
     470           G1%gdepw_ps(nxfin-nbghostcellsfine:nxfin,jj) = 0. 
     471        ENDIF 
     472     END DO 
     473     ! 
     474     DO ji=1,nxfin    ! South and North 
     475        IF (gdepw_ps_interp(ji,nbghostcellsfine+1)>0.) THEN 
     476           G1%gdepw_ps(ji,1:boundary) = gdepw_ps_interp(ji,1:boundary) 
     477           G1%wgt(ji,1:boundary) = 1. 
     478        ELSE 
     479           G1%gdepw_ps(ji,1:nbghostcellsfine+1)=0.  
     480        ENDIF 
     481        ! 
     482        IF (gdepw_ps_interp(ji,nyfin-nbghostcellsfine)>0.) THEN 
     483           G1%gdepw_ps(ji,nyfin-boundary+1:nyfin)=gdepw_ps_interp(ji,nyfin-boundary+1:nyfin) 
     484           G1%wgt(ji,nyfin-boundary+1:nyfin) = 1. 
     485        ELSE 
     486           G1%gdepw_ps(ji,nyfin-nbghostcellsfine:nyfin) = 0. 
     487        ENDIF 
     488     END DO 
     489     ! 
     490     !clem: recalculate interpolation everywhere before linear connection (useless to me) 
     491     IF( partial_steps .AND. ln_agrif_domain ) THEN                 
     492        gdepw_ps_interp = 0. 
     493        CALL Check_interp(G0,gdepw_ps_interp) 
     494     ENDIF 
     495     ! 
     496     ! ------------------------------------------------------- 
     497     ! Bathymetry between boundaries and interior (npt_connect)                  
     498     ! -------------------------------------------------------- 
     499     ! Make linear connection (on npt_connect*irafx points) between the boundaries and the interior 
     500     IF( ln_agrif_domain ) THEN 
     501        boundary = (npt_copy + npt_connect)*irafx + nbghostcellsfine + 1 
     502     ELSE 
     503        boundary = (npt_copy + npt_connect)*irafx 
     504     ENDIF 
     505 
     506     IF( npt_connect > 0 ) THEN 
     507        WRITE(*,*) ' linear connection on ',npt_connect,'coarse grid points' 
     508 
     509        wghts = 1. 
     510        DO ji = boundary - npt_connect*irafx + 1 , boundary 
     511           wghts = wghts - (1. / (npt_connect*irafx + 1. ) ) 
     512           DO jj=1,nyfin 
     513              IF (G1%gdepw_ps(nbghostcellsfine+1,jj) > 0.)       G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj))   
     514           END DO 
    486515        END DO 
    487         DO ji=1,nxfin 
    488           ! South  
    489           IF (gdepw_ps_interp(ji,nbghostcellsfine+1)>0.) THEN 
    490              G1%gdepw_ps(ji,1:boundary) = gdepw_ps_interp(ji,1:boundary) 
    491              G1%wgt(ji,1:boundary) = 1. 
    492           ELSE 
    493              G1%gdepw_ps(ji,1:nbghostcellsfine+1)=0.  
    494           ENDIF 
    495           ! North 
    496           IF (gdepw_ps_interp(ji,nyfin-nbghostcellsfine)>0.) THEN 
    497              G1%gdepw_ps(ji,nyfin-boundary+1:nyfin)=gdepw_ps_interp(ji,nyfin-boundary+1:nyfin) 
    498              G1%wgt(ji,nyfin-boundary+1:nyfin) = 1. 
    499           ELSE 
    500              G1%gdepw_ps(ji,nyfin-nbghostcellsfine:nyfin) = 0. 
    501           ENDIF 
     516 
     517        wghts = 1. 
     518        DO ji = nxfin - (boundary - npt_connect*irafx), nxfin - boundary +1 , -1 
     519           wghts = wghts - (1. / (npt_connect*irafx + 1. ) ) 
     520           DO jj=1,nyfin 
     521              IF (G1%gdepw_ps(nxfin-nbghostcellsfine,jj) > 0.)   G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj)) 
     522           END DO 
    502523        END DO 
    503          
    504 !        G1%gdepw_ps(1:boundary,:) = gdepw_ps_interp(1:boundary,:) 
    505 !        G1%gdepw_ps(:,1:boundary) = gdepw_ps_interp(:,1:boundary) 
    506 !        G1%gdepw_ps(nxfin-boundary+1:nxfin,:) = gdepw_ps_interp(nxfin-boundary+1:nxfin,:) 
    507 !        G1%gdepw_ps(:,nyfin-boundary+1:nyfin) = gdepw_ps_interp(:,nyfin-boundary+1:nyfin) 
    508         !                    
    509  
    510         IF(.NOT. smoothing) WRITE(*,*) 'No smoothing process only connection is carried out' 
    511         WRITE(*,*) ' linear connection on ',npt_connect,'coarse grid points' 
    512  
    513         !             
    514         gdepw_ps_interp = 0. 
    515         CALL Check_interp(G0,gdepw_ps_interp)      ! interpolation in connection zone (3 coarse grid cells) 
    516         ! 
    517         ! 
    518         ! 
    519         ! 
    520         !                    LINEAR CONNECTION 
    521         ! 
    522         ! 
    523         ! 
    524         ! 
    525         ! 
     524 
    526525        wghts = 1. 
    527         DO ji = boundary + 1 , boundary + npt_connect * irafx 
    528             wghts = wghts - (1. / (npt_connect*irafx - 1. ) ) 
    529             DO jj=1,nyfin 
    530                IF (G1%gdepw_ps(nbghostcellsfine+1,jj) > 0.) THEN 
    531                   G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj))   
    532                ENDIF 
    533             END DO 
    534         END DO  
    535        
     526        DO jj = boundary - npt_connect*irafy + 1 , boundary 
     527           wghts = wghts - (1. / (npt_connect*irafy + 1. ) ) 
     528           DO ji=1,nxfin 
     529              IF (G1%gdepw_ps(ji,nbghostcellsfine+1) > 0.)       G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj)) 
     530           END DO 
     531        END DO 
     532 
    536533        wghts = 1. 
    537         DO ji = nxfin - boundary , nxfin - boundary - npt_connect * irafx +1 , -1 
    538             wghts = wghts - (1. / (npt_connect*irafx - 1. ) ) 
    539             DO jj=1,nyfin 
    540                IF (G1%gdepw_ps(nxfin-nbghostcellsfine,jj) > 0.) THEN 
    541                   G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj)) 
    542                ENDIF 
    543             END DO 
    544         END DO   
    545  
    546         wghts = 1. 
    547         DO jj = boundary + 1 , boundary + npt_connect * irafy 
    548             wghts = wghts - (1. / (npt_connect*irafy - 1. ) ) 
    549             DO ji=1,nxfin 
    550                IF (G1%gdepw_ps(ji,nbghostcellsfine+1) > 0.) THEN 
    551                   G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj)) 
    552                ENDIF 
    553             END DO 
    554         END DO 
    555  
    556         wghts = 1. 
    557         DO jj = nyfin - boundary , nyfin - boundary - npt_connect * irafy +1, -1 
    558             wghts = wghts - (1. / (npt_connect*irafy - 1. ) ) 
    559             DO ji=1,nxfin 
    560                IF (G1%gdepw_ps(ji,nyfin-nbghostcellsfine) > 0.) THEN 
    561                   G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj)) 
    562                ENDIF 
    563             END DO 
     534        DO jj = nyfin - (boundary - npt_connect*irafy) , nyfin - boundary +1, -1 
     535           wghts = wghts - (1. / (npt_connect*irafy + 1. ) ) 
     536           DO ji=1,nxfin 
     537              IF (G1%gdepw_ps(ji,nyfin-nbghostcellsfine) > 0.)   G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj)) 
     538           END DO 
    564539        END DO 
    565540        IF (.NOT.identical_grids) THEN 
     
    567542        ENDIF 
    568543 
    569         G1%bathy_meter = G1%gdepw_ps 
    570         !                      
    571         !  
    572 ! Chanut: remove smoothing if child grid bathymetry is found to already exist 
    573         IF(smoothing.AND.(.NOT.identical_grids)) THEN  
    574  
    575            ! 
    576            ! Smoothing to connect the connection zone (3 + npt_connect coarse grid cells) and the interior domain 
    577            ! 
    578 ! Chanut: smoothing everywhere then discard result in connection zone 
    579            CALL smooth_topo(G1%gdepw_ps(1:nxfin,1:nyfin),nbiter) 
    580            WHERE (G1%wgt(:,:)==0.) G1%bathy_meter(:,:) = G1%gdepw_ps(:,:) 
    581 !           boundary = (npt_copy+npt_connect)*irafx + nbghostcellsfine + 1  
    582 !           CALL smooth_topo(G1%gdepw_ps(boundary:nxfin-boundary+1,boundary:nyfin-boundary+1),nbiter) 
    583 !           G1%bathy_meter = G1%gdepw_ps                          
    584         ENDIF 
    585  
    586  
    587         ! 
    588         !   
    589         ! Remove closed seas 
    590         !                             
    591         IF (removeclosedseas) THEN 
    592            ALLOCATE(bathy_test(nxfin,nyfin)) 
    593            bathy_test=0. 
    594            WHERE (G1%bathy_meter(1,:).GT.0.) 
    595               bathy_test(1,:)=1 
    596            END WHERE 
    597            WHERE (G1%bathy_meter(nxfin,:).GT.0.) 
    598               bathy_test(nxfin,:)=1 
    599            END WHERE 
    600            WHERE (G1%bathy_meter(:,1).GT.0.) 
    601               bathy_test(:,1)=1 
    602            END WHERE 
    603            WHERE (G1%bathy_meter(:,nyfin).GT.0.) 
    604               bathy_test(:,nyfin)=1 
    605            END WHERE 
    606            nbadd = 1 
    607            DO WHILE (nbadd.NE.0) 
    608               nbadd = 0 
    609               DO j=2,nyfin-1 
    610                  DO i=2,nxfin-1 
    611                     IF (G1%bathy_meter(i,j).GT.0.) THEN 
    612                        IF (MAX(bathy_test(i,j+1),bathy_test(i,j-1), & 
    613                             bathy_test(i-1,j),bathy_test(i+1,j)).EQ.1) THEN 
    614                           IF (bathy_test(i,j).NE.1.) nbadd = nbadd + 1 
    615                           bathy_test(i,j)=1. 
    616                        ENDIF 
    617  
     544     ENDIF 
     545 
     546     ! replace G1%bathy_meter by G1%gdepw_ps 
     547     G1%bathy_meter = G1%gdepw_ps 
     548     !                      
     549     ! -------------------- 
     550     ! Bathymetry smoothing                  
     551     ! -------------------- 
     552     IF( smoothing .AND. (.NOT.identical_grids) ) THEN 
     553        ! Chanut: smoothing everywhere then discard result in connection zone 
     554        CALL smooth_topo(G1%gdepw_ps(:,:),nbiter) 
     555        WHERE (G1%wgt(:,:)==0.) G1%bathy_meter(:,:) = G1%gdepw_ps(:,:) 
     556     ELSE 
     557        WRITE(*,*) 'No smoothing process only connection is carried out' 
     558     ENDIF 
     559     ! 
     560     ! ------------------ 
     561     ! Remove closed seas 
     562     ! ------------------ 
     563     IF (removeclosedseas) THEN 
     564        ALLOCATE(bathy_test(nxfin,nyfin)) 
     565        bathy_test=0. 
     566        WHERE (G1%bathy_meter(1,:)    .GT.0.)   bathy_test(1,:)=1 
     567        WHERE (G1%bathy_meter(nxfin,:).GT.0.)   bathy_test(nxfin,:)=1 
     568        WHERE (G1%bathy_meter(:,1)    .GT.0.)   bathy_test(:,1)=1 
     569        WHERE (G1%bathy_meter(:,nyfin).GT.0.)   bathy_test(:,nyfin)=1 
     570 
     571        nbadd = 1 
     572        DO WHILE (nbadd.NE.0) 
     573           nbadd = 0 
     574           DO j=2,nyfin-1 
     575              DO i=2,nxfin-1 
     576                 IF (G1%bathy_meter(i,j).GT.0.) THEN 
     577                    IF (MAX(bathy_test(i,j+1),bathy_test(i,j-1),bathy_test(i-1,j),bathy_test(i+1,j)).EQ.1) THEN 
     578                       IF (bathy_test(i,j).NE.1.) nbadd = nbadd + 1 
     579                       bathy_test(i,j)=1. 
    618580                    ENDIF 
    619                  ENDDO 
     581 
     582                 ENDIF 
    620583              ENDDO 
    621584           ENDDO 
    622            WHERE (bathy_test.EQ.0.) 
    623               G1%bathy_meter = 0. 
    624            END WHERE 
    625            DEALLOCATE(bathy_test) 
    626         ENDIF 
    627         ! 
    628         ! Chanut: Compute partial step bathymetry once more 
    629         CALL get_partial_steps(G1)                 ! compute gdepw_ps for G1 
    630  
    631         IF(bathy_update) CALL Update_Parent_Bathy( G0,G1 )  
    632         ! 
     585        ENDDO 
     586        WHERE (bathy_test.EQ.0.)   G1%bathy_meter = 0. 
     587        DEALLOCATE(bathy_test) 
     588     ENDIF 
     589     ! 
     590     !                              ! ---------------- 
     591     IF( partial_steps ) THEN       ! If partial steps 
     592        !                           ! ---------------- 
     593        ! 
     594        CALL get_partial_steps(G1)  ! recompute gdepw_ps for G1 
     595 
     596        IF(bathy_update)   CALL Update_Parent_Bathy( G0,G1 )  
    633597        CALL set_child_name(parent_bathy_meter,child_ps) 
    634         status = Write_Bathy_meter(TRIM(child_ps),G1)         
    635  
    636         IF(bathy_update) status = Write_Bathy_meter(TRIM(updated_parent_file),G0) 
     598         
     599                           status = Write_Bathy_meter(TRIM(child_ps),G1) 
     600        IF(bathy_update)   status = Write_Bathy_meter(TRIM(updated_parent_file),G0) 
     601        IF(bathy_update)   status = write_domcfg(TRIM(updated_parent_domcfg),G0) 
    637602 
    638603        CALL get_partial_steps(G1) 
     
    640605        G1%bathy_level=NINT(G1%bathy_level) 
    641606        ! 
    642         IF( TRIM(parent_meshmask_file) .NE. '/NULL' ) & 
    643              status = Write_Bathy_level(TRIM(Childlevel_file),G1) 
     607        ! store interpolation result in output file 
     608        CALL levels_to_meter(G1)   ! From levels to meters 
     609        IF( TRIM(parent_meshmask_file) .NE. '/NULL' )    status = Write_Bathy_level(TRIM(Childlevel_file),G1) 
     610        IF( TRIM(parent_meshmask_file) .NE. '/NULL' )    status = write_domcfg(TRIM(Child_domcfg),G1) 
    644611        ! 
    645612        WRITE(*,*) '****** Bathymetry successfully created for partial cells ******' 
    646         WRITE(*,*) ' ' 
    647         ! 
    648         STOP          
    649      ENDIF 
    650      !             
    651      !--------------------------------end if partial steps------------------------------------ 
    652      ! 
    653      ! 
    654      status = Read_bathy_level(TRIM(parent_meshmask_file),G0) 
    655      !             
    656      CALL levels_to_meter(G0) 
    657      !  
    658      ! compute constant bathymetry for Parent-Child bathymetry connection 
    659      !               
    660      WHERE( G0%bathy_meter < 0.000001 ) G0%bathy_meter = 0. 
    661  
    662      CALL init_constant_bathy(G0%bathy_meter,bathy_fin_constant) 
    663      ! 
    664      boundary = npt_copy*irafx + nbghostcellsfine + 1    
    665      !              
    666      G1%bathy_meter(1:boundary,:) = bathy_fin_constant(1:boundary,:) 
    667      G1%bathy_meter(:,1:boundary) = bathy_fin_constant(:,1:boundary) 
    668      G1%bathy_meter(nxfin-boundary+1:nxfin,:) = bathy_fin_constant(nxfin-boundary+1:nxfin,:) 
    669      G1%bathy_meter(:,nyfin-boundary+1:nyfin) = bathy_fin_constant(:,nyfin-boundary+1:nyfin) 
    670      ! 
    671      ! bathymetry smoothing 
    672      !                   
    673      CALL smooth_topo(G1%bathy_meter(boundary:nxfin-boundary+1,boundary:nyfin-boundary+1),nbiter) 
    674      ! 
    675      ! convert bathymetry from meters to levels 
    676      ! 
    677      CALL meter_to_levels(G1)  
    678      !            
    679      DEALLOCATE(G1%bathy_meter)            
    680      !               
     613        STOP 
     614        ! 
     615        !                           ! ---------------- 
     616     ELSE                           ! No partial steps 
     617        !                           ! ---------------- 
     618        ! 
     619        ! convert bathymetry from meters to levels 
     620        CALL meter_to_levels(G1)  
     621 
     622        IF(bathy_update)   CALL Update_Parent_Bathy( G0,G1 )  
     623        CALL set_child_name(parent_bathy_meter,child_ps) 
     624         
     625                           status = Write_Bathy_meter(TRIM(child_ps),G1)         
     626        IF(bathy_update)   status = Write_Bathy_meter(TRIM(updated_parent_file),G0) 
     627        IF(bathy_update)   status = write_domcfg(TRIM(updated_parent_domcfg),G0) 
     628        ! 
     629        ! store interpolation result in output file 
     630        CALL levels_to_meter(G1)   ! From levels to meters 
     631        status = Write_Bathy_level(TRIM(Childlevel_file),G1) 
     632        status = write_domcfg(TRIM(Child_domcfg),G1) 
     633        ! 
     634        WRITE(*,*) '****** Bathymetry successfully created for full cells ******' 
     635        STOP 
     636        !  
     637     ENDIF 
    681638  ENDIF 
    682639  ! 
    683   ! 
    684   ! make connection thanks to constant and interpolated bathymetry 
    685   ! 
    686   !       
    687   G1%bathy_level=NINT(G1%bathy_level) 
    688   !        
    689   DO j=1,nyfin 
    690      DO i=1,nxfin 
    691         IF (g1%bathy_level(i,j).LT.0.) THEN 
    692            PRINT *,'error in ',i,j,g1%bathy_level(i,j) 
    693         ENDIF 
    694      ENDDO 
    695   ENDDO 
    696   !        
    697   WHERE ((G1%bathy_level.LT.3.).AND.(G1%bathy_level.GT.0.)) 
    698      G1%bathy_level=3 
    699   END WHERE 
    700   ! 
    701   ! possibility to remove closed seas 
    702   !       
    703   IF (removeclosedseas) THEN 
    704      ALLOCATE(bathy_test(nxfin,nyfin)) 
    705  
    706      bathy_test=0. 
    707      WHERE (G1%bathy_level(1,:).GT.0.) 
    708         bathy_test(1,:)=1 
    709      END WHERE 
    710  
    711      WHERE (G1%bathy_level(nxfin,:).GT.0.) 
    712         bathy_test(nxfin,:)=1 
    713      END WHERE 
    714  
    715      WHERE (G1%bathy_level(:,1).GT.0.) 
    716         bathy_test(:,1)=1 
    717      END WHERE 
    718  
    719      WHERE (G1%bathy_level(:,nyfin).GT.0.) 
    720         bathy_test(:,nyfin)=1 
    721      END WHERE 
    722  
    723      nbadd = 1 
    724  
    725      DO WHILE (nbadd.NE.0) 
    726         nbadd = 0 
    727         DO j=2,nyfin-1 
    728            DO i=2,nxfin-1 
    729               IF (G1%bathy_level(i,j).GT.0.) THEN 
    730                  IF (MAX(bathy_test(i,j+1),bathy_test(i,j-1),bathy_test(i-1,j),bathy_test(i+1,j)).EQ.1) THEN 
    731                     IF (bathy_test(i,j).NE.1.) nbadd = nbadd + 1 
    732                     bathy_test(i,j)=1. 
    733                  ENDIF 
    734  
    735               ENDIF 
    736            ENDDO 
    737         ENDDO 
    738  
    739      ENDDO 
    740  
    741      WHERE (bathy_test.EQ.0.) 
    742         G1%bathy_level = 0. 
    743      END WHERE 
    744      DEALLOCATE(bathy_test)            
    745   ENDIF 
    746  
    747  
    748   ! 
    749   ! store interpolation result in output file 
    750   ! 
    751   status = Write_Bathy_level(TRIM(Childlevel_file),G1) 
    752  
    753   WRITE(*,*) '****** Bathymetry successfully created for full cells ******' 
    754   WRITE(*,*) ' ' 
    755  
    756   STOP 
    757640END PROGRAM create_bathy 
    758641 
  • utils/tools/NESTING/src/agrif_create_coordinates.f90

    r9632 r10025  
    3939  ! 
    4040  ! read input file (namelist.input) 
    41   ! 
    4241  CALL read_namelist(namelistname) 
    43   ! 
    44  
    4542  !       
    4643  ! read parent coodinates file 
    47   ! 
    4844  status = Read_Coordinates(TRIM(parent_coordinate_file),G0) 
    4945  !  
    5046  ! define name of child coordinate file 
    51   !      
    5247  CALL set_child_name(parent_coordinate_file,Child_filename)  
    5348  ! 
     
    6459  ! 
    6560  ! allocation of child grid elements 
    66   ! 
    6761  CALL agrif_grid_allocate(G1,nxfin,nyfin) 
    6862  ! 
    69   ! 
    7063  !  check potential longitude problems  
    71   !       
    7264  IF( G0%glamt(imin,jmin) > G0%glamt(imax,jmax) ) THEN            
    73      WRITE(*,*) '    ' 
    74      WHERE ( G0%glamt < 0 ) 
    75         G0%glamt = G0%glamt + 360. 
    76      END WHERE 
    77      WHERE ( G0%glamf < 0 ) 
    78         G0%glamf = G0%glamf + 360. 
    79      END WHERE 
    80   ENDIF 
    81  
     65     WHERE ( G0%glamt < 0 )   G0%glamt = G0%glamt + 360. 
     66     WHERE ( G0%glamf < 0 )   G0%glamf = G0%glamf + 360. 
     67  ENDIF 
    8268  !       
    8369  ! interpolation from parent grid to child grid for 
     
    8975  !    gphi = latitude 
    9076  ! 
    91  
    9277  !       
    93   !> M. Dunphy ticket 2082: 
    9478  CALL agrif_interp(G0%glamt,G1%glamt,'T') 
    9579  CALL agrif_interp(G0%glamf,G1%glamf,'F') 
    96 !  G1%glamu = G1%glamf 
    97 !  G1%glamv = G1%glamt 
    9880  CALL agrif_interp(G0%glamu,G1%glamu,'U') 
    9981  CALL agrif_interp(G0%glamv,G1%glamv,'V')  
     
    10183  CALL agrif_interp(G0%gphit,G1%gphit,'T') 
    10284  CALL agrif_interp(G0%gphif,G1%gphif,'F') 
    103 !  G1%gphiu = G1%gphit 
    104 !  G1%gphiv = G1%gphif 
    10585  CALL agrif_interp(G0%gphiu,G1%gphiu,'U') 
    10686  CALL agrif_interp(G0%gphiv,G1%gphiv,'V') 
    107   !<  M. Dunphy ticket 2082 
    10887  ! 
    10988  ! 
     
    11392  ! 
    11493  ! Compute scale factors e1 e2 
    115   ! 
    116   DO j=1,nyfin 
    117      DO i=2,nxfin 
    118         G1%e1t(i,j) = ra * rad * SQRT( (COS(rad*G1%gphit(i,j))*(G1%glamu(i,j)-G1%glamu(i-1,j)))**2 & 
    119              + (G1%gphiu(i,j)-G1%gphiu(i-1,j))**2) 
    120         G1%e1v(i,j) = ra * rad * SQRT( (COS(rad*G1%gphiv(i,j))*(G1%glamf(i,j)-G1%glamf(i-1,j)))**2 & 
    121              + (G1%gphif(i,j)-G1%gphif(i-1,j))**2)                                                              
    122      END DO 
    123   END DO 
    124   !       
    125   G1%e1t(1,:)=G1%e1t(2,:) 
    126   G1%e1v(1,:)=G1%e1v(2,:) 
    127   ! 
    128   DO j=1,nyfin 
    129      DO i=1,nxfin-1 
    130         G1%e1u(i,j) = ra * rad * SQRT( (COS(rad*G1%gphiu(i,j))*(G1%glamt(i+1,j)-G1%glamt(i,j)))**2 & 
    131              + (G1%gphit(i+1,j)-G1%gphit(i,j))**2) 
    132         G1%e1f(i,j) = ra * rad * SQRT( (COS(rad*G1%gphif(i,j))*(G1%glamv(i+1,j)-G1%glamv(i,j)))**2 & 
    133              + (G1%gphiv(i+1,j)-G1%gphiv(i,j))**2)                                                             
    134      END DO 
    135   END DO 
    136   ! 
    137   G1%e1u(nxfin,:)=G1%e1u(nxfin-1,:) 
    138   G1%e1f(nxfin,:)=G1%e1f(nxfin-1,:) 
    139   ! 
    140   DO j=2,nyfin 
    141      DO i=1,nxfin                                      
    142         G1%e2t(i,j) = ra * rad * SQRT( (COS(rad*G1%gphit(i,j))*(G1%glamv(i,j)-G1%glamv(i,j-1)))**2 & 
    143              + (G1%gphiv(i,j)-G1%gphiv(i,j-1))**2) 
    144         G1%e2u(i,j) = ra * rad * SQRT( (COS(rad*G1%gphiu(i,j))*(G1%glamf(i,j)-G1%glamf(i,j-1)))**2 & 
    145              + (G1%gphif(i,j)-G1%gphif(i,j-1))**2)                                                                
    146      END DO 
    147   END DO 
    148   ! 
    149   G1%e2t(:,1)=G1%e2t(:,2) 
    150   G1%e2u(:,1)=G1%e2u(:,2) 
    151   ! 
    152   DO j=1,nyfin-1 
    153      DO i=1,nxfin 
    154         G1%e2v(i,j) = ra * rad * SQRT( (COS(rad*G1%gphiv(i,j))*(G1%glamt(i,j+1)-G1%glamt(i,j)))**2 & 
    155              + (G1%gphit(i,j+1)-G1%gphit(i,j))**2) 
    156         G1%e2f(i,j) = ra * rad * SQRT( (COS(rad*G1%gphif(i,j))*(G1%glamu(i,j+1)-G1%glamu(i,j)))**2 & 
    157              + (G1%gphiu(i,j+1)-G1%gphiu(i,j))**2) 
    158      END DO 
    159   END DO 
    160   ! 
    161   G1%e2v(:,nyfin)=G1%e2v(:,nyfin-1) 
    162   G1%e2f(:,nyfin)=G1%e2f(:,nyfin-1) 
    163  
    164  
    165   CALL agrif_interp(G0%e1t,G1%e1t,'T') 
    166   G1%e1t = G1%e1t / REAL(irafx) 
    167   CALL agrif_interp(G0%e2t,G1%e2t,'T') 
    168   G1%e2t = G1%e2t / REAL(irafy) 
    169  
    170   CALL agrif_interp(G0%e1u,G1%e1u,'U') 
    171   G1%e1u = G1%e1u / REAL(irafx) 
    172   CALL agrif_interp(G0%e2u,G1%e2u,'U') 
    173   G1%e2u = G1%e2u / REAL(irafy)   
    174  
    175   CALL agrif_interp(G0%e1v,G1%e1v,'V') 
    176   G1%e1v = G1%e1v / REAL(irafx) 
    177   CALL agrif_interp(G0%e2v,G1%e2v,'V') 
    178   G1%e2v = G1%e2v / REAL(irafy)  
    179  
    180   CALL agrif_interp(G0%e1f,G1%e1f,'F') 
    181   G1%e1f = G1%e1f / REAL(irafx) 
    182   CALL agrif_interp(G0%e2f,G1%e2f,'F') 
    183   G1%e2f = G1%e2f / REAL(irafy)                
    184  
    185   ! 
    186   WHERE ( G1%glamt > 180 ) 
    187      G1%glamt = G1%glamt - 360. 
    188   END WHERE 
    189   WHERE ( G1%glamf > 180 ) 
    190      G1%glamf = G1%glamf - 360. 
    191   END WHERE 
    192   WHERE ( G1%glamu > 180 ) 
    193      G1%glamu = G1%glamu - 360. 
    194   END WHERE 
    195   WHERE ( G1%glamv > 180 ) 
    196      G1%glamv = G1%glamv - 360. 
    197   END WHERE 
    198   ! 
     94!  DO j=1,nyfin 
     95!     DO i=2,nxfin 
     96!        G1%e1t(i,j) = ra * rad * SQRT( (COS(rad*G1%gphit(i,j))*(G1%glamu(i,j)-G1%glamu(i-1,j)))**2 & 
     97!             + (G1%gphiu(i,j)-G1%gphiu(i-1,j))**2) 
     98!        G1%e1v(i,j) = ra * rad * SQRT( (COS(rad*G1%gphiv(i,j))*(G1%glamf(i,j)-G1%glamf(i-1,j)))**2 & 
     99!             + (G1%gphif(i,j)-G1%gphif(i-1,j))**2)                                                              
     100!     END DO 
     101!  END DO 
     102!  !       
     103!  G1%e1t(1,:)=G1%e1t(2,:) 
     104!  G1%e1v(1,:)=G1%e1v(2,:) 
     105!  ! 
     106!  DO j=1,nyfin 
     107!     DO i=1,nxfin-1 
     108!        G1%e1u(i,j) = ra * rad * SQRT( (COS(rad*G1%gphiu(i,j))*(G1%glamt(i+1,j)-G1%glamt(i,j)))**2 & 
     109!             + (G1%gphit(i+1,j)-G1%gphit(i,j))**2) 
     110!        G1%e1f(i,j) = ra * rad * SQRT( (COS(rad*G1%gphif(i,j))*(G1%glamv(i+1,j)-G1%glamv(i,j)))**2 & 
     111!             + (G1%gphiv(i+1,j)-G1%gphiv(i,j))**2)                                                             
     112!     END DO 
     113!  END DO 
     114!  ! 
     115!  G1%e1u(nxfin,:)=G1%e1u(nxfin-1,:) 
     116!  G1%e1f(nxfin,:)=G1%e1f(nxfin-1,:) 
     117!  ! 
     118!  DO j=2,nyfin 
     119!     DO i=1,nxfin                                      
     120!        G1%e2t(i,j) = ra * rad * SQRT( (COS(rad*G1%gphit(i,j))*(G1%glamv(i,j)-G1%glamv(i,j-1)))**2 & 
     121!             + (G1%gphiv(i,j)-G1%gphiv(i,j-1))**2) 
     122!        G1%e2u(i,j) = ra * rad * SQRT( (COS(rad*G1%gphiu(i,j))*(G1%glamf(i,j)-G1%glamf(i,j-1)))**2 & 
     123!             + (G1%gphif(i,j)-G1%gphif(i,j-1))**2)                                                                
     124!     END DO 
     125!  END DO 
     126!  ! 
     127!  G1%e2t(:,1)=G1%e2t(:,2) 
     128!  G1%e2u(:,1)=G1%e2u(:,2) 
     129!  ! 
     130!  DO j=1,nyfin-1 
     131!     DO i=1,nxfin 
     132!        G1%e2v(i,j) = ra * rad * SQRT( (COS(rad*G1%gphiv(i,j))*(G1%glamt(i,j+1)-G1%glamt(i,j)))**2 & 
     133!             + (G1%gphit(i,j+1)-G1%gphit(i,j))**2) 
     134!        G1%e2f(i,j) = ra * rad * SQRT( (COS(rad*G1%gphif(i,j))*(G1%glamu(i,j+1)-G1%glamu(i,j)))**2 & 
     135!             + (G1%gphiu(i,j+1)-G1%gphiu(i,j))**2) 
     136!     END DO 
     137!  END DO 
     138!  ! 
     139!  G1%e2v(:,nyfin)=G1%e2v(:,nyfin-1) 
     140!  G1%e2f(:,nyfin)=G1%e2f(:,nyfin-1) 
     141 
     142 
     143  CALL agrif_interp(G0%e1t,G1%e1t,'T')   ;   G1%e1t = G1%e1t / REAL(irafx) 
     144  CALL agrif_interp(G0%e2t,G1%e2t,'T')   ;   G1%e2t = G1%e2t / REAL(irafy) 
     145 
     146  CALL agrif_interp(G0%e1u,G1%e1u,'U')   ;   G1%e1u = G1%e1u / REAL(irafx) 
     147  CALL agrif_interp(G0%e2u,G1%e2u,'U')   ;   G1%e2u = G1%e2u / REAL(irafy)   
     148 
     149  CALL agrif_interp(G0%e1v,G1%e1v,'V')   ;   G1%e1v = G1%e1v / REAL(irafx) 
     150  CALL agrif_interp(G0%e2v,G1%e2v,'V')   ;   G1%e2v = G1%e2v / REAL(irafy)  
     151 
     152  CALL agrif_interp(G0%e1f,G1%e1f,'F')   ;   G1%e1f = G1%e1f / REAL(irafx) 
     153  CALL agrif_interp(G0%e2f,G1%e2f,'F')   ;   G1%e2f = G1%e2f / REAL(irafy) 
     154  ! 
     155  WHERE ( G1%glamt > 180 )   G1%glamt = G1%glamt - 360. 
     156  WHERE ( G1%glamf > 180 )   G1%glamf = G1%glamf - 360. 
     157  WHERE ( G1%glamu > 180 )   G1%glamu = G1%glamu - 360. 
     158  WHERE ( G1%glamv > 180 )   G1%glamv = G1%glamv - 360. 
    199159  ! 
    200160  G1%nav_lon=G1%glamt 
    201161  G1%nav_lat=G1%gphit  
    202162  !               
    203   ! 
    204163  ! Write interpolation result in child coodinates file 
    205   !       
    206164  status = Write_Coordinates(TRIM(Child_filename),G1) 
    207  
    208165  ! 
    209166  WRITE(*,*) 'Child domain position : ' 
    210   WRITE(*,*) 'latmin =',G1%gphit(3,3) 
    211   WRITE(*,*) 'latmax =',G1%gphit(nxfin-2,nyfin-2) 
    212   WRITE(*,*) 'lonmin =',G1%glamt(3,3) 
    213   WRITE(*,*) 'lonmax =',G1%glamt(nxfin-2,nyfin-2) 
     167  IF( ln_agrif_domain ) THEN 
     168     WRITE(*,*) 'latmin =',G1%gphit(3,3) 
     169     WRITE(*,*) 'latmax =',G1%gphit(nxfin-2,nyfin-2) 
     170     WRITE(*,*) 'lonmin =',G1%glamt(3,3) 
     171     WRITE(*,*) 'lonmax =',G1%glamt(nxfin-2,nyfin-2) 
     172  ELSE 
     173     WRITE(*,*) 'latmin =',G1%gphit(1,1) 
     174     WRITE(*,*) 'latmax =',G1%gphit(nxfin,nyfin) 
     175     WRITE(*,*) 'lonmin =',G1%glamt(1,1) 
     176     WRITE(*,*) 'lonmax =',G1%glamt(nxfin,nyfin) 
     177  ENDIF 
    214178  STOP 
    215179END PROGRAM create_coordinate 
  • utils/tools/NESTING/src/agrif_create_data.f90

    r2143 r10025  
    1010  !************************************************************************ 
    1111  !                           * 
    12   ! PROGRAM  CREATE_DATA                     * 
     12  ! PROGRAM  CREATE_DATA                       * 
    1313  !                           * 
    1414  ! program to implement data interpolation to generate        * 
    1515  ! child grid forcing files                 * 
    16   !                           *                          * 
     16  !                           * 
    1717  !Interpolation is carried out using bilinear interpolation      * 
    1818  !routine from SCRIP package                *      
    1919  !                           * 
    20   !http://climate.lanl.gov/Software/SCRIP/            *                          * 
     20  !http://climate.lanl.gov/Software/SCRIP/            * 
    2121  !************************************************************************ 
    2222  ! 
    23   INTEGER :: narg,iargc,i 
     23  INTEGER           :: narg, iargc, ji 
    2424  CHARACTER(len=80) :: namelistname 
    2525 
     
    3333 
    3434  ! read input file (namelist.input) 
    35   ! 
    3635  CALL read_namelist(namelistname) 
    37  
    38   i = 1 
    3936  ! 
    4037  ! Interpolate U grid  data 
    41   ! 
    42   DO WHILE( TRIM(U_Files(i)) .NE. '/NULL' ) 
    43      PRINT *,'Grid U forcing files = ',u_files(i) 
     38  ji = 1 
     39  DO WHILE( TRIM(U_Files(ji)) .NE. '/NULL' ) 
     40     PRINT *,'Grid U forcing files = ',u_files(ji) 
    4441     !         
    45      CALL Interp_Extrap_var(U_FILES(i), 'U')  
    46      i = i+1                 
     42     CALL Interp_Extrap_var(U_FILES(ji), 'U')  
     43     ji = ji+1                 
    4744     !              
    4845  END DO 
    4946 
    50   i = 1 
    5147  ! 
    5248  ! Interpolate V grid  data 
    53   ! 
    54   DO WHILE( TRIM(V_Files(i)) .NE. '/NULL' ) 
    55      PRINT *,'Grid V forcing files = ',v_files(i) 
     49  ji = 1 
     50  DO WHILE( TRIM(V_Files(ji)) .NE. '/NULL' ) 
     51     PRINT *,'Grid V forcing files = ',v_files(ji) 
    5652     !         
    57      CALL Interp_Extrap_var(V_FILES(i), 'V')  
    58      i = i+1                 
     53     CALL Interp_Extrap_var(V_FILES(ji), 'V')  
     54     ji = ji+1                 
    5955     !              
    6056  END DO 
    61  
    62   i = 1 
    6357  ! 
    6458  ! Interpolate flux data 
    65   ! 
    66   DO WHILE( TRIM(Flx_Files(i)) .NE. '/NULL' ) 
    67      PRINT *,'flxfiles = ',flx_files(i) 
     59  ji = 1 
     60  DO WHILE( TRIM(Flx_Files(ji)) .NE. '/NULL' ) 
     61     PRINT *,'flxfiles = ',flx_files(ji) 
    6862     !         
    69      CALL Interp_Extrap_var(FLX_FILES(i), 'T')  
    70      i = i+1                 
     63     CALL Interp_Extrap_var(FLX_FILES(ji), 'T')  
     64     ji = ji+1                 
    7165     !              
    7266  END DO 
  • utils/tools/NESTING/src/agrif_create_restart.f90

    r6140 r10025  
    389389           WRITE(*,*) 'copy time_steps' 
    390390           CALL Read_Ncdf_var('time_steps',TRIM(restart_file),tabtemp1DInt)  
    391            CALL Write_Ncdf_var('time_steps',TRIM(timedimname),Child_file,tabtemp1DInt) 
     391           CALL Write_Ncdf_var('time_steps',TRIM(timedimname),Child_file,tabtemp1DInt,'integer') 
    392392           CALL Copy_Ncdf_att('time_steps',TRIM(restart_file),Child_file)  
    393393           DEALLOCATE(tabtemp1DInt) 
  • utils/tools/NESTING/src/agrif_interpolation.f90

    r9166 r10025  
    5454    !  
    5555    ! 
    56     CALL agrif_base_interp2(tabin,tabout,imin,jmin,typevar) 
    57  
     56    IF( ln_agrif_domain ) THEN 
     57       CALL agrif_base_interp2(tabin,tabout,imin,jmin,typevar) 
     58    ELSE 
     59       CALL agrif_base_interp3(tabin,tabout,imin,jmin,typevar) 
     60    ENDIF 
     61     
    5862    ! 
    5963  END SUBROUTINE agrif_interp 
    6064  ! 
    61   !******************************************************** 
    62   !   subroutine agrif_base_interp2       * 
    63   !******************************************************** 
    6465  !       
    6566  SUBROUTINE agrif_base_interp2(tabin,tabout,i_min,j_min,typevar) 
     
    189190  END SUBROUTINE agrif_base_interp2 
    190191  ! 
    191   !******************************************************** 
    192   !   subroutine polint             * 
    193   !******************************************************** 
     192  !       
     193  SUBROUTINE agrif_base_interp3(tabin,tabout,i_min,j_min,typevar) 
     194    ! 
     195    IMPLICIT NONE 
     196    REAL*8,DIMENSION(:,:) :: tabin,tabout 
     197    INTEGER :: i_min,j_min 
     198    CHARACTER(*) :: typevar    
     199 
     200    INTEGER :: nxf,nyf,zx,zy 
     201    INTEGER :: ji,jj,jif,jjf,jic,jjc,jdecx,jdecy 
     202    REAL*8 :: Ax, Bx, Ay, By 
     203 
     204    nxf = SIZE(tabout,1) 
     205    nyf = SIZE(tabout,2) 
     206 
     207    SELECT CASE(typevar) 
     208    CASE('T') 
     209       IF(MOD(irafx,2)==1) THEN ! odd 
     210          zx = 1 ; zy = 1 ; jdecx = FLOOR(irafx/2.) ; jdecy = FLOOR(irafy/2.) 
     211       ELSE                     ! even 
     212          zx = 2 ; zy = 2 ; jdecx = FLOOR(irafx/2.) ; jdecy = FLOOR(irafy/2.) 
     213       ENDIF 
     214    CASE('U') 
     215       IF(MOD(irafx,2)==1) THEN ! odd 
     216          zx = 1 ; zy = 1 ; jdecx = irafx - 1 ; jdecy = FLOOR(irafy/2.) 
     217       ELSE                     ! even 
     218          zx = 1 ; zy = 2 ; jdecx = irafx - 1 ; jdecy = FLOOR(irafy/2.) 
     219       ENDIF 
     220    CASE('V') 
     221       IF(MOD(irafx,2)==1) THEN ! odd 
     222          zx = 1 ; zy = 1 ; jdecx = FLOOR(irafx/2.) ; jdecy = irafy - 1 
     223       ELSE                     ! even 
     224          zx = 2 ; zy = 1 ; jdecx = FLOOR(irafx/2.) ; jdecy = irafy - 1 
     225       ENDIF 
     226    CASE('F') 
     227       IF(MOD(irafx,2)==1) THEN ! odd 
     228          zx = 1 ; zy = 1 ; jdecx = irafx - 1 ; jdecy = irafy - 1 
     229       ELSE                     ! even 
     230          zx = 1 ; zy = 1 ; jdecx = irafx - 1 ; jdecy = irafy - 1 
     231       ENDIF 
     232    END SELECT 
     233 
     234     
     235    DO jj = 1, nyf 
     236 
     237       jjf = jj - jdecy 
     238       jjc = j_min + FLOOR((jjf-1.) / irafy)  
     239        
     240       DO ji = 1, nxf 
     241           
     242          jif = ji - jdecx  
     243          jic = i_min + FLOOR((jif-1.) / irafx)  
     244           
     245          Bx = MOD( zx*jif-1, zx*irafx ) / REAL(zx*irafx) 
     246          By = MOD( zy*jjf-1, zy*irafy ) / REAL(zy*irafy) 
     247          Ax = 1. - Bx 
     248          Ay = 1. - By 
     249 
     250          tabout(ji,jj) = ( Bx * tabin(jic+1,jjc  ) + Ax * tabin(jic,jjc  ) ) * Ay + & 
     251             &            ( Bx * tabin(jic+1,jjc+1) + Ax * tabin(jic,jjc+1) ) * By 
     252       END DO 
     253    END DO 
     254     
     255    ! 
     256  END SUBROUTINE agrif_base_interp3 
     257 
    194258  !        
    195259  SUBROUTINE polint(xin,valin,n,x,val) 
     
    238302  END SUBROUTINE polint 
    239303  ! 
    240   !******************************************************** 
    241   !   subroutine polcoe             * 
    242   !******************************************************** 
    243304  !       
    244305  SUBROUTINE polcoe(xin,valin,n,cof) 
     
    709770          CALL Read_Ncdf_var('time_steps',filename,tabtemp1DInt)  
    710771          !       print *,'timedeph = ',tabtemp1DInt 
    711           CALL Write_Ncdf_var('time_steps','time_counter',Child_file,tabtemp1DInt) 
     772          CALL Write_Ncdf_var('time_steps','time_counter',Child_file,tabtemp1DInt,'integer') 
    712773          CALL Copy_Ncdf_att('time_steps',filename,Child_file)  
    713774          DEALLOCATE(tabtemp1DInt) 
     
    9841045  END SUBROUTINE Interp_Extrap_var 
    9851046  ! 
    986   !************************************************************** 
    987   !   end subroutine Interp_var 
    988   !**************************************************************           
    9891047  !  
    9901048END MODULE agrif_interpolation 
  • utils/tools/NESTING/src/agrif_partial_steps.f90

    r9694 r10025  
    276276    REAL,DIMENSION(N) :: gdepw,e3t 
    277277    REAL :: za0,za1,za2,zsur,zacr,zacr2,zkth,zkth2,zmin,zmax,zdepth 
    278     INTEGER :: kbathy,jk,diff 
    279     INTEGER :: bornex,borney,bornex2,borney2 
     278    INTEGER :: kbathy,jk 
    280279    ! 
    281280    IF ( ( pa0 == 0 .OR. pa1 == 0 .OR. psur == 0 ) & 
     
    348347       e3t(N) = e3t(N-1) 
    349348    END IF 
    350     ! 
    351     diff = 0       
    352     IF ( MOD(irafx,2) .EQ. 0 ) diff = 1 
    353     !        
    354     bornex = nbghostcellsfine + 1 + CEILING(irafx/2.0) + diff - irafx 
    355     borney = nbghostcellsfine + 1 + CEILING(irafy/2.0) + diff - irafy 
    356     bornex2 = nxfin - nbghostcellsfine - irafx - CEILING(irafx/2.0)  
    357     borney2 = nyfin - nbghostcellsfine - irafy - CEILING(irafy/2.0)                       
    358     ! 
     349    !               
    359350    ! 
    360351    ! west boundary 
    361     ! 
    362  
    363     CALL correct_level( gdepwchild,ParentGrid,gdepw,e3t,1,2+nbghostcellsfine+(npt_copy+npt_connect)*irafx-1, & 
    364          1,nyfin) 
    365  
     352    IF( ln_agrif_domain ) THEN 
     353       CALL correct_level( gdepwchild,ParentGrid,gdepw,e3t,1,2+nbghostcellsfine+(npt_copy+npt_connect)*irafx-1,1,nyfin) 
     354    ELSE 
     355       CALL correct_level( gdepwchild,ParentGrid,gdepw,e3t,1,(npt_copy+npt_connect)*irafx,1,nyfin) 
     356    ENDIF 
    366357    ! 
    367358    ! east boundary 
    368     ! 
    369  
    370     CALL correct_level( gdepwchild,ParentGrid,gdepw,e3t,nxfin-1-nbghostcellsfine-((npt_copy+npt_connect)*irafx-1),nxfin, & 
    371          1,nyfin) 
    372  
     359    IF( ln_agrif_domain ) THEN 
     360       CALL correct_level( gdepwchild,ParentGrid,gdepw,e3t,nxfin-1-nbghostcellsfine-((npt_copy+npt_connect)*irafx-1),nxfin,1,nyfin) 
     361    ELSE 
     362       CALL correct_level( gdepwchild,ParentGrid,gdepw,e3t,nxfin-((npt_copy+npt_connect)*irafx+1),nxfin,1,nyfin) 
     363    ENDIF 
    373364    ! 
    374365    ! north boundary 
    375     ! 
    376  
    377     CALL correct_level( gdepwchild,ParentGrid,gdepw,e3t,1,nxfin, & 
    378          nyfin-1 - nbghostcellsfine -((npt_copy+npt_connect)*irafy-1),nyfin ) 
    379  
     366    IF( ln_agrif_domain ) THEN 
     367       CALL correct_level( gdepwchild,ParentGrid,gdepw,e3t,1,nxfin,nyfin-1-nbghostcellsfine-((npt_copy+npt_connect)*irafy-1),nyfin) 
     368    ELSE 
     369       CALL correct_level( gdepwchild,ParentGrid,gdepw,e3t,1,nxfin,nyfin-((npt_copy+npt_connect)*irafy+1),nyfin) 
     370    ENDIF 
    380371    ! 
    381372    ! south boundary 
    382     ! 
    383     CALL correct_level( gdepwchild,ParentGrid,gdepw,e3t,1,nxfin, & 
    384          1,2+nbghostcellsfine+(npt_copy+npt_connect)*irafy-1 ) 
    385  
     373    IF( ln_agrif_domain ) THEN 
     374       CALL correct_level( gdepwchild,ParentGrid,gdepw,e3t,1,nxfin,1,2+nbghostcellsfine+(npt_copy+npt_connect)*irafy-1) 
     375    ELSE 
     376       CALL correct_level( gdepwchild,ParentGrid,gdepw,e3t,1,nxfin,1,(npt_copy+npt_connect)*irafy) 
     377    ENDIF 
    386378    !        
    387379    ! 
     
    398390    INTEGER :: kbathy,jk,indx,indy,diff 
    399391    REAL :: xdiff 
    400     INTEGER :: i,j,ji,ij,ii,jj,jpt,ipt 
     392    INTEGER :: i,j,ji,ij,ii,jj,jpt,ipt,i1,i2,j1,j2,ii1,ii2,jj1,jj2 
    401393    REAL*8 :: slopex, slopey,val,tmp1,tmp2,tmp3,tmp4 
    402394    INTEGER :: parentbathy 
     
    471463    jend = pty + agrif_int((y-ymin-dyfin/2.)/dyfin)                
    472464 
    473  
    474     ALLOCATE(gdepwtemp(ibegin-irafx:iend+irafx,jbegin-irafy:jend+irafy)) 
    475     ALLOCATE(parentbathytab(ibegin-irafx:iend+irafx,jbegin-irafy:jend+irafy)) 
    476  
     465    IF( ln_agrif_domain ) THEN 
     466       ALLOCATE(gdepwtemp(ibegin-irafx:iend+irafx,jbegin-irafy:jend+irafy)) 
     467       ALLOCATE(parentbathytab(ibegin-irafx:iend+irafx,jbegin-irafy:jend+irafy)) 
     468 
     469       i1 = ibegin 
     470       i2 = iend 
     471       j1 = jbegin 
     472       j2 = jend 
     473        
     474       ii1 = -FLOOR(irafx/2.0)+diff 
     475       ii2 =  FLOOR(irafx/2.0) 
     476       jj1 = -FLOOR(irafy/2.0)+diff 
     477       jj2 =  FLOOR(irafy/2.0) 
     478    ELSE 
     479       ibegin = minboundx 
     480       jbegin = minboundy 
     481       iend   = maxboundx ! (npt_copy+npt_connect)*irafx 
     482       jend   = maxboundy ! (npt_copy+npt_connect)*irafy 
     483       ! 
     484       ipbegin = imin + (ibegin-1)/irafx 
     485       jpbegin = jmin + (jbegin-1)/irafy 
     486       ipend   = ipbegin + (npt_copy+npt_connect) - 1 
     487       jpend   = jpbegin + (npt_copy+npt_connect) - 1 
     488       ! 
     489       i1 = ibegin 
     490       i2 = iend 
     491       j1 = jbegin 
     492       j2 = jend 
     493        
     494       ii1 = 0 
     495       ii2 = irafx - 1 
     496       jj1 = 0 
     497       jj2 = irafy - 1 
     498       ! 
     499       ALLOCATE(gdepwtemp(ibegin:iend,jbegin:jend)) 
     500       ALLOCATE(parentbathytab(ibegin:iend,jbegin:jend)) 
     501 
     502    ENDIF 
     503     
    477504 
    478505    jpt=jpbegin 
     
    482509 
    483510 
    484        DO i=ibegin,iend,irafx 
     511       DO i=i1,i2,irafx 
    485512 
    486513 
     
    524551          ! interpolation on fine grid points (connection zone) 
    525552          ! 
    526           DO ii = i-FLOOR(irafx/2.0)+diff,i+FLOOR(irafx/2.0) 
    527              x = ii-i - xdiff/2. 
     553          DO ii = i+ii1,i+ii2 
     554!!             x = ii-i - xdiff/2. 
    528555!!             val = parentgrid%gdepw_ps(ipt,jpt)+slopex * x 
    529556!! chanut: uncomment this to get nearest neighbor interpolation 
     
    563590 
    564591          slopey = vanleer(gdepwtemp(i,j-irafy:j+irafy:irafy))/REAL(irafy) 
    565  
     592           
    566593          tmp1 = (maxdepth - gdepwtemp(i,j)) / REAL(irafy) 
    567594          tmp2 = (gdepwtemp(i,j) - mindepth) / REAL(irafy) 
     
    583610 
    584611 
    585           DO jj = j-FLOOR(irafy/2.0)+diff,j+FLOOR(irafy/2.0) 
    586              y = jj-j - xdiff/2. 
     612          DO jj = j+jj1,j+jj2 
     613!!             y = jj-j - xdiff/2. 
    587614!!             val = gdepwtemp(i,j) + slopey*y 
    588615!! chanut: uncomment this to get nearest neighbor interpolation 
  • utils/tools/NESTING/src/agrif_readwrite.f90

    r9628 r10025  
    4040    TYPE(Coordinates) :: Grid 
    4141    !       
    42     CALL Read_Ncdf_var('glamt',name,Grid%glamt) 
    43     CALL Read_Ncdf_var('glamu',name,Grid%glamu) 
    44     CALL Read_Ncdf_var('glamv',name,Grid%glamv) 
    45     CALL Read_Ncdf_var('glamf',name,Grid%glamf) 
    46     CALL Read_Ncdf_var('gphit',name,Grid%gphit) 
    47     CALL Read_Ncdf_var('gphiu',name,Grid%gphiu) 
    48     CALL Read_Ncdf_var('gphiv',name,Grid%gphiv) 
    49     CALL Read_Ncdf_var('gphif',name,Grid%gphif) 
    50     CALL Read_Ncdf_var('e1t',name,Grid%e1t) 
    51     CALL Read_Ncdf_var('e1u',name,Grid%e1u) 
    52     CALL Read_Ncdf_var('e1v',name,Grid%e1v) 
    53     CALL Read_Ncdf_var('e1f',name,Grid%e1f) 
    54     CALL Read_Ncdf_var('e2t',name,Grid%e2t) 
    55     CALL Read_Ncdf_var('e2u',name,Grid%e2u) 
    56     CALL Read_Ncdf_var('e2v',name,Grid%e2v) 
    57     CALL Read_Ncdf_var('e2f',name,Grid%e2f) 
    58     CALL Read_Ncdf_var('nav_lon',name,Grid%nav_lon) 
    59     CALL Read_Ncdf_var('nav_lat',name,Grid%nav_lat)        
     42    CALL read_ncdf_var('glamt',name,Grid%glamt) 
     43    CALL read_ncdf_var('glamu',name,Grid%glamu) 
     44    CALL read_ncdf_var('glamv',name,Grid%glamv) 
     45    CALL read_ncdf_var('glamf',name,Grid%glamf) 
     46    CALL read_ncdf_var('gphit',name,Grid%gphit) 
     47    CALL read_ncdf_var('gphiu',name,Grid%gphiu) 
     48    CALL read_ncdf_var('gphiv',name,Grid%gphiv) 
     49    CALL read_ncdf_var('gphif',name,Grid%gphif) 
     50    CALL read_ncdf_var('e1t',name,Grid%e1t) 
     51    CALL read_ncdf_var('e1u',name,Grid%e1u) 
     52    CALL read_ncdf_var('e1v',name,Grid%e1v) 
     53    CALL read_ncdf_var('e1f',name,Grid%e1f) 
     54    CALL read_ncdf_var('e2t',name,Grid%e2t) 
     55    CALL read_ncdf_var('e2u',name,Grid%e2u) 
     56    CALL read_ncdf_var('e2v',name,Grid%e2v) 
     57    CALL read_ncdf_var('e2f',name,Grid%e2f) 
     58    CALL read_ncdf_var('nav_lon',name,Grid%nav_lon) 
     59    CALL read_ncdf_var('nav_lat',name,Grid%nav_lat)        
    6060    !  
    6161    IF( PRESENT(Pacifique) )THEN 
     
    103103    TYPE(Coordinates) :: Grid 
    104104    !       
    105     CALL Read_Ncdf_var('glamt',name,Grid%glamt,strt,cnt) 
    106     CALL Read_Ncdf_var('glamu',name,Grid%glamu,strt,cnt) 
    107     CALL Read_Ncdf_var('glamv',name,Grid%glamv,strt,cnt) 
    108     CALL Read_Ncdf_var('glamf',name,Grid%glamf,strt,cnt) 
    109     CALL Read_Ncdf_var('gphit',name,Grid%gphit,strt,cnt) 
    110     CALL Read_Ncdf_var('gphiu',name,Grid%gphiu,strt,cnt) 
    111     CALL Read_Ncdf_var('gphiv',name,Grid%gphiv,strt,cnt) 
    112     CALL Read_Ncdf_var('gphif',name,Grid%gphif,strt,cnt) 
    113     CALL Read_Ncdf_var('e1t',name,Grid%e1t,strt,cnt) 
    114     CALL Read_Ncdf_var('e1u',name,Grid%e1u,strt,cnt) 
    115     CALL Read_Ncdf_var('e1v',name,Grid%e1v,strt,cnt) 
    116     CALL Read_Ncdf_var('e1f',name,Grid%e1f,strt,cnt) 
    117     CALL Read_Ncdf_var('e2t',name,Grid%e2t,strt,cnt) 
    118     CALL Read_Ncdf_var('e2u',name,Grid%e2u,strt,cnt) 
    119     CALL Read_Ncdf_var('e2v',name,Grid%e2v,strt,cnt) 
    120     CALL Read_Ncdf_var('e2f',name,Grid%e2f,strt,cnt) 
    121     CALL Read_Ncdf_var('nav_lon',name,Grid%nav_lon,strt,cnt) 
    122     CALL Read_Ncdf_var('nav_lat',name,Grid%nav_lat,strt,cnt)        
     105    CALL read_ncdf_var('glamt',name,Grid%glamt,strt,cnt) 
     106    CALL read_ncdf_var('glamu',name,Grid%glamu,strt,cnt) 
     107    CALL read_ncdf_var('glamv',name,Grid%glamv,strt,cnt) 
     108    CALL read_ncdf_var('glamf',name,Grid%glamf,strt,cnt) 
     109    CALL read_ncdf_var('gphit',name,Grid%gphit,strt,cnt) 
     110    CALL read_ncdf_var('gphiu',name,Grid%gphiu,strt,cnt) 
     111    CALL read_ncdf_var('gphiv',name,Grid%gphiv,strt,cnt) 
     112    CALL read_ncdf_var('gphif',name,Grid%gphif,strt,cnt) 
     113    CALL read_ncdf_var('e1t',name,Grid%e1t,strt,cnt) 
     114    CALL read_ncdf_var('e1u',name,Grid%e1u,strt,cnt) 
     115    CALL read_ncdf_var('e1v',name,Grid%e1v,strt,cnt) 
     116    CALL read_ncdf_var('e1f',name,Grid%e1f,strt,cnt) 
     117    CALL read_ncdf_var('e2t',name,Grid%e2t,strt,cnt) 
     118    CALL read_ncdf_var('e2u',name,Grid%e2u,strt,cnt) 
     119    CALL read_ncdf_var('e2v',name,Grid%e2v,strt,cnt) 
     120    CALL read_ncdf_var('e2f',name,Grid%e2f,strt,cnt) 
     121    CALL read_ncdf_var('nav_lon',name,Grid%nav_lon,strt,cnt) 
     122    CALL read_ncdf_var('nav_lat',name,Grid%nav_lat,strt,cnt)        
    123123    ! 
    124124    WRITE(*,*) ' ' 
     
    146146    !            
    147147    dimnames = (/ 'x','y' /) 
    148     CALL Write_Ncdf_dim(dimnames(1),name,nxfin) 
    149     CALL Write_Ncdf_dim(dimnames(2),name,nyfin) 
    150     !       
    151     CALL Write_Ncdf_var('nav_lon',dimnames,name,Grid%nav_lon,'float')       
    152     CALL Write_Ncdf_var('nav_lat',dimnames,name,Grid%nav_lat,'float') 
    153     ! 
    154     CALL Write_Ncdf_var('glamt',dimnames,name,Grid%glamt,'double') 
    155     CALL Write_Ncdf_var('glamu',dimnames,name,Grid%glamu,'double') 
    156     CALL Write_Ncdf_var('glamv',dimnames,name,Grid%glamv,'double') 
    157     CALL Write_Ncdf_var('glamf',dimnames,name,Grid%glamf,'double') 
    158     CALL Write_Ncdf_var('gphit',dimnames,name,Grid%gphit,'double') 
    159     CALL Write_Ncdf_var('gphiu',dimnames,name,Grid%gphiu,'double') 
    160     CALL Write_Ncdf_var('gphiv',dimnames,name,Grid%gphiv,'double') 
    161     CALL Write_Ncdf_var('gphif',dimnames,name,Grid%gphif,'double')       
    162     CALL Write_Ncdf_var('e1t',dimnames,name,Grid%e1t,'double')       
    163     CALL Write_Ncdf_var('e1u',dimnames,name,Grid%e1u,'double')      
    164     CALL Write_Ncdf_var('e1v',dimnames,name,Grid%e1v,'double')       
    165     CALL Write_Ncdf_var('e1f',dimnames,name,Grid%e1f,'double') 
    166     CALL Write_Ncdf_var('e2t',dimnames,name,Grid%e2t,'double') 
    167     CALL Write_Ncdf_var('e2u',dimnames,name,Grid%e2u,'double') 
    168     CALL Write_Ncdf_var('e2v',dimnames,name,Grid%e2v,'double') 
    169     CALL Write_Ncdf_var('e2f',dimnames,name,Grid%e2f,'double') 
    170     !       
    171     CALL Copy_Ncdf_att('nav_lon',TRIM(parent_coordinate_file),name,MINVAL(Grid%nav_lon),MAXVAL(Grid%nav_lon)) 
    172     CALL Copy_Ncdf_att('nav_lat',TRIM(parent_coordinate_file),name,MINVAL(Grid%nav_lat),MAXVAL(Grid%nav_lat)) 
    173     CALL Copy_Ncdf_att('glamt',TRIM(parent_coordinate_file),name) 
    174     CALL Copy_Ncdf_att('glamu',TRIM(parent_coordinate_file),name) 
    175     CALL Copy_Ncdf_att('glamv',TRIM(parent_coordinate_file),name) 
    176     CALL Copy_Ncdf_att('glamf',TRIM(parent_coordinate_file),name) 
    177     CALL Copy_Ncdf_att('gphit',TRIM(parent_coordinate_file),name) 
    178     CALL Copy_Ncdf_att('gphiu',TRIM(parent_coordinate_file),name) 
    179     CALL Copy_Ncdf_att('gphiv',TRIM(parent_coordinate_file),name) 
    180     CALL Copy_Ncdf_att('gphif',TRIM(parent_coordinate_file),name) 
    181     CALL Copy_Ncdf_att('e1t',TRIM(parent_coordinate_file),name) 
    182     CALL Copy_Ncdf_att('e1u',TRIM(parent_coordinate_file),name) 
    183     CALL Copy_Ncdf_att('e1v',TRIM(parent_coordinate_file),name) 
    184     CALL Copy_Ncdf_att('e1f',TRIM(parent_coordinate_file),name) 
    185     CALL Copy_Ncdf_att('e2t',TRIM(parent_coordinate_file),name) 
    186     CALL Copy_Ncdf_att('e2u',TRIM(parent_coordinate_file),name) 
    187     CALL Copy_Ncdf_att('e2v',TRIM(parent_coordinate_file),name) 
    188     CALL Copy_Ncdf_att('e2f',TRIM(parent_coordinate_file),name)             
     148    CALL write_ncdf_dim(dimnames(1),name,nxfin) 
     149    CALL write_ncdf_dim(dimnames(2),name,nyfin) 
     150    !       
     151    CALL write_ncdf_var('nav_lon',dimnames,name,Grid%nav_lon,'float')       
     152    CALL write_ncdf_var('nav_lat',dimnames,name,Grid%nav_lat,'float') 
     153    ! 
     154    CALL write_ncdf_var('glamt',dimnames,name,Grid%glamt,'double') 
     155    CALL write_ncdf_var('glamu',dimnames,name,Grid%glamu,'double') 
     156    CALL write_ncdf_var('glamv',dimnames,name,Grid%glamv,'double') 
     157    CALL write_ncdf_var('glamf',dimnames,name,Grid%glamf,'double') 
     158    CALL write_ncdf_var('gphit',dimnames,name,Grid%gphit,'double') 
     159    CALL write_ncdf_var('gphiu',dimnames,name,Grid%gphiu,'double') 
     160    CALL write_ncdf_var('gphiv',dimnames,name,Grid%gphiv,'double') 
     161    CALL write_ncdf_var('gphif',dimnames,name,Grid%gphif,'double')       
     162    CALL write_ncdf_var('e1t',dimnames,name,Grid%e1t,'double')       
     163    CALL write_ncdf_var('e1u',dimnames,name,Grid%e1u,'double')      
     164    CALL write_ncdf_var('e1v',dimnames,name,Grid%e1v,'double')       
     165    CALL write_ncdf_var('e1f',dimnames,name,Grid%e1f,'double') 
     166    CALL write_ncdf_var('e2t',dimnames,name,Grid%e2t,'double') 
     167    CALL write_ncdf_var('e2u',dimnames,name,Grid%e2u,'double') 
     168    CALL write_ncdf_var('e2v',dimnames,name,Grid%e2v,'double') 
     169    CALL write_ncdf_var('e2f',dimnames,name,Grid%e2f,'double') 
     170    !       
     171    CALL copy_ncdf_att('nav_lon',TRIM(parent_coordinate_file),name,MINVAL(Grid%nav_lon),MAXVAL(Grid%nav_lon)) 
     172    CALL copy_ncdf_att('nav_lat',TRIM(parent_coordinate_file),name,MINVAL(Grid%nav_lat),MAXVAL(Grid%nav_lat)) 
     173    CALL copy_ncdf_att('glamt',TRIM(parent_coordinate_file),name) 
     174    CALL copy_ncdf_att('glamu',TRIM(parent_coordinate_file),name) 
     175    CALL copy_ncdf_att('glamv',TRIM(parent_coordinate_file),name) 
     176    CALL copy_ncdf_att('glamf',TRIM(parent_coordinate_file),name) 
     177    CALL copy_ncdf_att('gphit',TRIM(parent_coordinate_file),name) 
     178    CALL copy_ncdf_att('gphiu',TRIM(parent_coordinate_file),name) 
     179    CALL copy_ncdf_att('gphiv',TRIM(parent_coordinate_file),name) 
     180    CALL copy_ncdf_att('gphif',TRIM(parent_coordinate_file),name) 
     181    CALL copy_ncdf_att('e1t',TRIM(parent_coordinate_file),name) 
     182    CALL copy_ncdf_att('e1u',TRIM(parent_coordinate_file),name) 
     183    CALL copy_ncdf_att('e1v',TRIM(parent_coordinate_file),name) 
     184    CALL copy_ncdf_att('e1f',TRIM(parent_coordinate_file),name) 
     185    CALL copy_ncdf_att('e2t',TRIM(parent_coordinate_file),name) 
     186    CALL copy_ncdf_att('e2u',TRIM(parent_coordinate_file),name) 
     187    CALL copy_ncdf_att('e2v',TRIM(parent_coordinate_file),name) 
     188    CALL copy_ncdf_att('e2f',TRIM(parent_coordinate_file),name)             
    189189    ! 
    190190    WRITE(*,*) ' ' 
     
    209209    TYPE(Coordinates) :: Grid 
    210210    ! 
    211     CALL Read_Ncdf_var('mbathy',name,Grid%Bathy_level)     
     211    CALL read_ncdf_var('mbathy',name,Grid%Bathy_level)     
    212212    ! 
    213213    WRITE(*,*) ' ' 
     
    232232    CHARACTER(len=1),DIMENSION(2) :: dimnames 
    233233    ! 
    234     status = nf90_create(name,NF90_NOCLOBBER,ncid) 
     234    status = nf90_create(name,NF90_WRITE,ncid) 
    235235    status = nf90_close(ncid)           
    236236    ! 
    237237    dimnames = (/ 'x','y' /) 
    238     CALL Write_Ncdf_dim(dimnames(1),name,nxfin) 
    239     CALL Write_Ncdf_dim(dimnames(2),name,nyfin) 
    240     !       
    241     CALL Write_Ncdf_var('nav_lon',dimnames,name,Grid%nav_lon    ,'float') 
    242     CALL Write_Ncdf_var('nav_lat',dimnames,name,Grid%nav_lat    ,'float') 
    243     CALL Write_Ncdf_var('mbathy' ,dimnames,name,Grid%bathy_level,'float') 
    244     ! 
    245     CALL Copy_Ncdf_att('nav_lon',TRIM(parent_meshmask_file),name,MINVAL(Grid%nav_lon),MAXVAL(Grid%nav_lon)) 
    246     CALL Copy_Ncdf_att('nav_lat',TRIM(parent_meshmask_file),name,MINVAL(Grid%nav_lat),MAXVAL(Grid%nav_lat)) 
    247     CALL Copy_Ncdf_att('mbathy' ,TRIM(parent_meshmask_file),name)        
     238    CALL write_ncdf_dim(dimnames(1),name,nxfin) 
     239    CALL write_ncdf_dim(dimnames(2),name,nyfin) 
     240    !       
     241    CALL write_ncdf_var('nav_lon',dimnames,name,Grid%nav_lon    ,'float') 
     242    CALL write_ncdf_var('nav_lat',dimnames,name,Grid%nav_lat    ,'float') 
     243    CALL write_ncdf_var('mbathy' ,dimnames,name,Grid%bathy_level,'float') 
     244    ! 
     245    CALL copy_ncdf_att('nav_lon',TRIM(parent_meshmask_file),name,MINVAL(Grid%nav_lon),MAXVAL(Grid%nav_lon)) 
     246    CALL copy_ncdf_att('nav_lat',TRIM(parent_meshmask_file),name,MINVAL(Grid%nav_lat),MAXVAL(Grid%nav_lat)) 
     247    CALL copy_ncdf_att('mbathy' ,TRIM(parent_meshmask_file),name)        
    248248    ! 
    249249    WRITE(*,*) ' ' 
     
    276276       WRITE(*,*) ' etopo format for external high resolution database  ' 
    277277       WRITE(*,*) '****' 
    278        CALL Read_Ncdf_var('lon',name,topo_lon) 
    279        CALL Read_Ncdf_var('lat',name,topo_lat) 
     278       CALL read_ncdf_var('lon',name,topo_lon) 
     279       CALL read_ncdf_var('lat',name,topo_lat) 
    280280    ELSE IF( Dims_Existence('x',name) .AND. Dims_Existence('y',name) ) THEN 
    281281       WRITE(*,*) '****' 
    282282       WRITE(*,*) ' OPA format for external high resolution database  ' 
    283283       WRITE(*,*) '****' 
    284        CALL Read_Ncdf_var('nav_lon',name,CoarseGrid%nav_lon) 
    285        CALL Read_Ncdf_var('nav_lat',name,CoarseGrid%nav_lat) 
    286        CALL Read_Ncdf_var(parent_batmet_name,name,CoarseGrid%Bathy_meter) 
     284       CALL read_ncdf_var('nav_lon',name,CoarseGrid%nav_lon) 
     285       CALL read_ncdf_var('nav_lat',name,CoarseGrid%nav_lat) 
     286       CALL read_ncdf_var(parent_batmet_name,name,CoarseGrid%Bathy_meter) 
    287287       !             
    288288       IF ( PRESENT(Pacifique) ) THEN 
     
    303303    ENDIF 
    304304    ! 
    305     IF( MAXVAL(ChildGrid%glamt) > 180 ) THEN                  
     305    IF( MAXVAL(ChildGrid%glamt) > 180. ) THEN                  
    306306       !           
    307        WHERE( topo_lon < 0 ) 
    308           topo_lon = topo_lon + 360. 
    309        END WHERE 
     307       WHERE( topo_lon < 0. )   topo_lon = topo_lon + 360. 
    310308       !           
    311309       i_min = MAXLOC(topo_lon,mask = topo_lon < MINVAL(ChildGrid%nav_lon)-zdel) 
     
    316314       tabdim1 = ( SIZE(topo_lon) - i_min(1) + 1 ) + i_max(1)                     
    317315       !           
    318        IF(j_min(1)-2 >= 1 .AND. j_max(1)+3 <= SIZE(topo_lat,1) ) THEN 
    319           j_min(1) = j_min(1)-2 
    320           j_max(1) = j_max(1)+3 
     316       IF( ln_agrif_domain ) THEN 
     317          IF(j_min(1)-2 >= 1 .AND. j_max(1)+3 <= SIZE(topo_lat,1) ) THEN 
     318             j_min(1) = j_min(1)-2 
     319             j_max(1) = j_max(1)+3 
     320          ENDIF 
    321321       ENDIF 
    322322       tabdim2 = j_max(1) - j_min(1) + 1 
     
    352352    ELSE 
    353353       ! 
     354       WHERE( topo_lon > 180. )   topo_lon = topo_lon - 360. 
     355       ! 
    354356       i_min = MAXLOC(topo_lon,mask = topo_lon < MINVAL(ChildGrid%nav_lon)-zdel) 
    355357       i_max = MINLOC(topo_lon,mask = topo_lon > MAXVAL(ChildGrid%nav_lon)+zdel) 
     
    357359       j_max = MINLOC(topo_lat,mask = topo_lat > MAXVAL(ChildGrid%nav_lat)+zdel) 
    358360       !       
    359        IF(i_min(1)-2 >= 1 .AND. i_max(1)+3 <= SIZE(topo_lon,1) ) THEN 
    360           i_min(1) = i_min(1)-2 
    361           i_max(1) = i_max(1)+3 
     361       IF( ln_agrif_domain ) THEN 
     362          IF(i_min(1)-2 >= 1 .AND. i_max(1)+3 <= SIZE(topo_lon,1) ) THEN 
     363             i_min(1) = i_min(1)-2 
     364             i_max(1) = i_max(1)+3 
     365          ENDIF 
    362366       ENDIF 
    363367       tabdim1 = i_max(1) - i_min(1) + 1 
    364368       ! 
    365        IF(j_min(1)-2 >= 1 .AND. j_max(1)+3 <= SIZE(topo_lat,1) ) THEN 
    366           j_min(1) = j_min(1)-2 
    367           j_max(1) = j_max(1)+3 
     369       IF( ln_agrif_domain ) THEN 
     370          IF(j_min(1)-2 >= 1 .AND. j_max(1)+3 <= SIZE(topo_lat,1) ) THEN 
     371             j_min(1) = j_min(1)-2 
     372             j_max(1) = j_max(1)+3 
     373          ENDIF 
    368374       ENDIF 
    369375       tabdim2 = j_max(1) - j_min(1) + 1 
     
    423429    TYPE(Coordinates) :: Grid 
    424430    ! 
    425     CALL Read_Ncdf_var(parent_batmet_name,name,Grid%Bathy_meter)     
     431    CALL read_ncdf_var(parent_batmet_name,name,Grid%Bathy_meter)     
    426432    ! 
    427433    WRITE(*,*) ' ' 
     
    454460    dimnames = (/ 'x','y' /) 
    455461 
    456     CALL Write_Ncdf_dim(dimnames(1),name,nx) 
    457     CALL Write_Ncdf_dim(dimnames(2),name,ny) 
     462    CALL write_ncdf_dim(dimnames(1),name,nx) 
     463    CALL write_ncdf_dim(dimnames(2),name,ny) 
    458464    !      
    459     CALL Write_Ncdf_var('nav_lon'         ,dimnames,name,Grid%nav_lon    ,'float') 
    460     CALL Write_Ncdf_var('nav_lat'         ,dimnames,name,Grid%nav_lat    ,'float') 
    461     CALL Write_Ncdf_var(parent_batmet_name,dimnames,name,Grid%bathy_meter,'float') 
    462     CALL Write_Ncdf_var('weight'          ,dimnames,name,Grid%wgt        ,'float') 
    463     ! 
    464     CALL Copy_Ncdf_att('nav_lon'         ,TRIM(parent_bathy_meter),name,MINVAL(Grid%nav_lon),MAXVAL(Grid%nav_lon)) 
    465     CALL Copy_Ncdf_att('nav_lat'         ,TRIM(parent_bathy_meter),name,MINVAL(Grid%nav_lat),MAXVAL(Grid%nav_lat))   
    466     CALL Copy_Ncdf_att(parent_batmet_name,TRIM(parent_bathy_meter),name) 
     465    CALL write_ncdf_var('nav_lon'         ,dimnames,name,Grid%nav_lon    ,'float') 
     466    CALL write_ncdf_var('nav_lat'         ,dimnames,name,Grid%nav_lat    ,'float') 
     467    CALL write_ncdf_var(parent_batmet_name,dimnames,name,Grid%bathy_meter,'float') 
     468    CALL write_ncdf_var('weight'          ,dimnames,name,Grid%wgt        ,'float') 
     469    ! 
     470    CALL copy_ncdf_att('nav_lon'         ,TRIM(parent_bathy_meter),name,MINVAL(Grid%nav_lon),MAXVAL(Grid%nav_lon)) 
     471    CALL copy_ncdf_att('nav_lat'         ,TRIM(parent_bathy_meter),name,MINVAL(Grid%nav_lat),MAXVAL(Grid%nav_lat))   
     472    CALL copy_ncdf_att(parent_batmet_name,TRIM(parent_bathy_meter),name) 
    467473    ! 
    468474    WRITE(*,*) ' ' 
     
    474480  END FUNCTION Write_Bathy_meter 
    475481  ! 
     482  !***************************************************** 
     483  !   function write_domcfg(name,Grid) 
     484  !***************************************************** 
     485 
     486  INTEGER FUNCTION write_domcfg(name,Grid) 
     487     !----------------------------------------- 
     488     ! It creates a domain_cfg.nc used in NEMO4 
     489     !-----------------------------------------     
     490     ! 
     491     USE io_netcdf 
     492     ! 
     493     CHARACTER(*) name 
     494     TYPE(Coordinates) :: Grid 
     495     ! 
     496     INTEGER :: status, ncid 
     497     INTEGER :: nx, ny, jk 
     498     INTEGER :: ln_sco, ln_isfcav, ln_zco, ln_zps, jperio 
     499     REAL*8  :: rpi, rad, rday, rsiyea, rsiday, omega 
     500     ! 
     501     CHARACTER(len=1),     DIMENSION(3)     ::   dimnames 
     502     REAL*8 ,              DIMENSION(N)     ::   e3t_1d, e3w_1d, gdept_1d 
     503     REAL*8 , ALLOCATABLE, DIMENSION(:,:)   ::   ff_t, ff_f 
     504     INTEGER, ALLOCATABLE, DIMENSION(:,:)   ::   bottom_level, top_level 
     505     REAL*8 , ALLOCATABLE, DIMENSION(:,:,:) ::   e3t_0, e3u_0, e3v_0, e3f_0, e3w_0, e3uw_0, e3vw_0 
     506     ! 
     507     ! size of the Grid 
     508     nx = SIZE(Grid%bathy_meter,1) 
     509     ny = SIZE(Grid%bathy_meter,2) 
     510 
     511     ! allocate needed arrays for domain_cfg 
     512     ALLOCATE( ff_t(nx,ny), ff_f(nx,ny) ) 
     513     ALLOCATE( bottom_level(nx,ny), top_level(nx,ny) ) 
     514     ALLOCATE( e3t_0(nx,ny,N), e3u_0 (nx,ny,N), e3v_0 (nx,ny,N), e3f_0(nx,ny,N),  & 
     515        &      e3w_0(nx,ny,N), e3uw_0(nx,ny,N), e3vw_0(nx,ny,N) ) 
     516 
     517     ! some physical parameters 
     518     rpi    = 3.141592653589793 
     519     rad    = 3.141592653589793 / 180. 
     520     rday   = 24.*60.*60. 
     521     rsiyea = 365.25 * rday * 2. * rpi / 6.283076 
     522     rsiday = rday / ( 1. + rday / rsiyea ) 
     523     omega  = 2. * rpi / rsiday 
     524 
     525     ! Coriolis 
     526     ff_f(:,:) = 2. * omega * SIN( rad * Grid%gphif(:,:) )     ! compute it on the sphere at f-point 
     527     ff_t(:,:) = 2. * omega * SIN( rad * Grid%gphit(:,:) )     !    -        -       -    at t-point 
     528 
     529     ! top/bottom levels 
     530     bottom_level(:,:) = Grid%bathy_level(:,:) 
     531     top_level(:,:)    = MIN( 1, bottom_level(:,:) ) 
     532 
     533     ! vertical scale factors 
     534     CALL zgr_z( e3t_1d, e3w_1d, gdept_1d )     
     535     DO jk = 1, N 
     536        e3t_0  (:,:,jk) = e3t_1d  (jk) 
     537        e3u_0  (:,:,jk) = e3t_1d  (jk) 
     538        e3v_0  (:,:,jk) = e3t_1d  (jk) 
     539        e3f_0  (:,:,jk) = e3t_1d  (jk) 
     540        e3w_0  (:,:,jk) = e3w_1d  (jk) 
     541        e3uw_0 (:,:,jk) = e3w_1d  (jk) 
     542        e3vw_0 (:,:,jk) = e3w_1d  (jk) 
     543     END DO 
     544 
     545     ! logicals and others 
     546     ln_sco = 0 
     547     ln_isfcav = 0 
     548     IF( partial_steps ) THEN 
     549        ln_zps = 1 
     550        ln_zco = 0 
     551     ELSE 
     552        ln_zps = 0 
     553        ln_zco = 1 
     554     ENDIF 
     555 
     556     ! closed domain (agrif) 
     557     jperio = 0 
     558 
     559     ! ------------------- 
     560     ! write domain_cfg.nc 
     561     ! ------------------- 
     562     status = nf90_create(name,NF90_WRITE,ncid) 
     563     status = nf90_close(ncid) 
     564     ! 
     565     ! dimensions 
     566     dimnames = (/'x','y','z'/) 
     567     CALL write_ncdf_dim(dimnames(1),name,nx) 
     568     CALL write_ncdf_dim(dimnames(2),name,ny) 
     569     CALL write_ncdf_dim(dimnames(3),name,N) 
     570     !       
     571     ! variables 
     572     CALL write_ncdf_var('nav_lon',dimnames(1:2),name,Grid%nav_lon,'float')       
     573     CALL write_ncdf_var('nav_lat',dimnames(1:2),name,Grid%nav_lat,'float') 
     574     CALL write_ncdf_var('nav_lev',dimnames(3)  ,name,gdept_1d    ,'float') 
     575     ! 
     576     CALL write_ncdf_var('jpiglo',name,nx    ,'integer') 
     577     CALL write_ncdf_var('jpjglo',name,ny    ,'integer') 
     578     CALL write_ncdf_var('jpkglo',name,N     ,'integer') 
     579     CALL write_ncdf_var('jperio',name,jperio,'integer') 
     580     ! 
     581     CALL write_ncdf_var('ln_zco'   ,name,ln_zco   ,'integer') 
     582     CALL write_ncdf_var('ln_zps'   ,name,ln_zps   ,'integer') 
     583     CALL write_ncdf_var('ln_sco'   ,name,ln_sco   ,'integer') 
     584     CALL write_ncdf_var('ln_isfcav',name,ln_isfcav,'integer') 
     585 
     586     CALL write_ncdf_var('glamt',dimnames(1:2),name,Grid%glamt,'double') 
     587     CALL write_ncdf_var('glamu',dimnames(1:2),name,Grid%glamu,'double') 
     588     CALL write_ncdf_var('glamv',dimnames(1:2),name,Grid%glamv,'double') 
     589     CALL write_ncdf_var('glamf',dimnames(1:2),name,Grid%glamf,'double') 
     590     CALL write_ncdf_var('gphit',dimnames(1:2),name,Grid%gphit,'double') 
     591     CALL write_ncdf_var('gphiu',dimnames(1:2),name,Grid%gphiu,'double') 
     592     CALL write_ncdf_var('gphiv',dimnames(1:2),name,Grid%gphiv,'double') 
     593     CALL write_ncdf_var('gphif',dimnames(1:2),name,Grid%gphif,'double')       
     594 
     595     CALL write_ncdf_var('e1t',dimnames(1:2),name,Grid%e1t,'double')       
     596     CALL write_ncdf_var('e1u',dimnames(1:2),name,Grid%e1u,'double')      
     597     CALL write_ncdf_var('e1v',dimnames(1:2),name,Grid%e1v,'double')       
     598     CALL write_ncdf_var('e1f',dimnames(1:2),name,Grid%e1f,'double') 
     599     CALL write_ncdf_var('e2t',dimnames(1:2),name,Grid%e2t,'double') 
     600     CALL write_ncdf_var('e2u',dimnames(1:2),name,Grid%e2u,'double') 
     601     CALL write_ncdf_var('e2v',dimnames(1:2),name,Grid%e2v,'double') 
     602     CALL write_ncdf_var('e2f',dimnames(1:2),name,Grid%e2f,'double') 
     603 
     604     CALL write_ncdf_var('ff_f',dimnames(1:2),name,ff_f,'double') 
     605     CALL write_ncdf_var('ff_t',dimnames(1:2),name,ff_t,'double') 
     606 
     607     CALL write_ncdf_var('e3t_1d',dimnames(3),name,e3t_1d,'double') 
     608     CALL write_ncdf_var('e3w_1d',dimnames(3),name,e3w_1d,'double') 
     609 
     610     CALL write_ncdf_var('e3t_0' ,dimnames(:),name,e3t_0 ,'double') 
     611     CALL write_ncdf_var('e3w_0' ,dimnames(:),name,e3w_0 ,'double') 
     612     CALL write_ncdf_var('e3u_0' ,dimnames(:),name,e3u_0 ,'double') 
     613     CALL write_ncdf_var('e3v_0' ,dimnames(:),name,e3v_0 ,'double') 
     614     CALL write_ncdf_var('e3f_0' ,dimnames(:),name,e3f_0 ,'double') 
     615     CALL write_ncdf_var('e3w_0' ,dimnames(:),name,e3w_0 ,'double') 
     616     CALL write_ncdf_var('e3uw_0',dimnames(:),name,e3uw_0,'double') 
     617     CALL write_ncdf_var('e3vw_0',dimnames(:),name,e3vw_0,'double') 
     618 
     619     CALL write_ncdf_var('bottom_level',dimnames(1:2),name,bottom_level,'integer') 
     620     CALL write_ncdf_var('top_level'   ,dimnames(1:2),name,top_level   ,'integer') 
     621 
     622     CALL write_ncdf_var('bathy_meter' ,dimnames(1:2),name,Grid%bathy_meter,'float') 
     623 
     624     ! some attributes       
     625     CALL copy_ncdf_att('nav_lon',TRIM(parent_coordinate_file),name,MINVAL(Grid%nav_lon),MAXVAL(Grid%nav_lon)) 
     626     CALL copy_ncdf_att('nav_lat',TRIM(parent_coordinate_file),name,MINVAL(Grid%nav_lat),MAXVAL(Grid%nav_lat)) 
     627 
     628     CALL copy_ncdf_att('glamt',TRIM(parent_coordinate_file),name) 
     629     CALL copy_ncdf_att('glamu',TRIM(parent_coordinate_file),name) 
     630     CALL copy_ncdf_att('glamv',TRIM(parent_coordinate_file),name) 
     631     CALL copy_ncdf_att('glamf',TRIM(parent_coordinate_file),name) 
     632     CALL copy_ncdf_att('gphit',TRIM(parent_coordinate_file),name) 
     633     CALL copy_ncdf_att('gphiu',TRIM(parent_coordinate_file),name) 
     634     CALL copy_ncdf_att('gphiv',TRIM(parent_coordinate_file),name) 
     635     CALL copy_ncdf_att('gphif',TRIM(parent_coordinate_file),name) 
     636 
     637     CALL copy_ncdf_att('e1t',TRIM(parent_coordinate_file),name) 
     638     CALL copy_ncdf_att('e1u',TRIM(parent_coordinate_file),name) 
     639     CALL copy_ncdf_att('e1v',TRIM(parent_coordinate_file),name) 
     640     CALL copy_ncdf_att('e1f',TRIM(parent_coordinate_file),name) 
     641     CALL copy_ncdf_att('e2t',TRIM(parent_coordinate_file),name) 
     642     CALL copy_ncdf_att('e2u',TRIM(parent_coordinate_file),name) 
     643     CALL copy_ncdf_att('e2v',TRIM(parent_coordinate_file),name) 
     644     CALL copy_ncdf_att('e2f',TRIM(parent_coordinate_file),name)             
     645     ! 
     646     ! control print 
     647     WRITE(*,*) ' ' 
     648     WRITE(*,*) 'Writing domcfg file: ',name 
     649     WRITE(*,*) ' ' 
     650     ! 
     651     DEALLOCATE( ff_t, ff_f ) 
     652     DEALLOCATE( bottom_level, top_level ) 
     653     DEALLOCATE( e3t_0, e3u_0, e3v_0, e3f_0, e3w_0, e3uw_0, e3vw_0 ) 
     654     !       
     655     write_domcfg = 1 
     656 
     657  END FUNCTION write_domcfg 
     658  ! 
     659  SUBROUTINE zgr_z(e3t_1d, e3w_1d, gdept_1d) 
     660     !!---------------------------------------------------------------------- 
     661     !!                   ***  ROUTINE zgr_z (from NEMO4) *** 
     662     !!                    
     663     !! ** Purpose :   set the depth of model levels and the resulting  
     664     !!      vertical scale factors. 
     665     !! 
     666     !! ** Method  :   z-coordinate system (use in all type of coordinate) 
     667     !!        The depth of model levels is defined from an analytical 
     668     !!      function the derivative of which gives the scale factors. 
     669     !!        both depth and scale factors only depend on k (1d arrays). 
     670     !!              w-level: gdepw_1d  = gdep(k) 
     671     !!                       e3w_1d(k) = dk(gdep)(k)     = e3(k) 
     672     !!              t-level: gdept_1d  = gdep(k+0.5) 
     673     !!                       e3t_1d(k) = dk(gdep)(k+0.5) = e3(k+0.5) 
     674     !! 
     675     !! ** Action  : - gdept_1d, gdepw_1d : depth of T- and W-point (m) 
     676     !!              - e3t_1d  , e3w_1d   : scale factors at T- and W-levels (m) 
     677     !! 
     678     !! Reference : Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766. 
     679     !!---------------------------------------------------------------------- 
     680     INTEGER  ::   jk                   ! dummy loop indices 
     681     INTEGER  ::   jpk 
     682     REAL*8 ::   zt, zw                 ! temporary scalars 
     683     REAL*8 ::   zsur, za0, za1, zkth   ! Values set from parameters in 
     684     REAL*8 ::   zacr, zdzmin, zhmax    ! par_CONFIG_Rxx.h90 
     685     REAL*8 ::   za2, zkth2, zacr2      ! Values for optional double tanh function set from parameters  
     686 
     687     REAL*8, DIMENSION(:), INTENT(inout) ::  e3t_1d, e3w_1d, gdept_1d 
     688     REAL*8, DIMENSION(N)                ::  gdepw_1d 
     689     !!---------------------------------------------------------------------- 
     690     ! 
     691     ! 
     692     ! Set variables from parameters 
     693     ! ------------------------------ 
     694     zkth = ppkth       ;   zacr = ppacr 
     695     zdzmin = ppdzmin   ;   zhmax = pphmax 
     696     zkth2 = ppkth2     ;   zacr2 = ppacr2   ! optional (ldbletanh=T) double tanh parameters 
     697 
     698     ! If pa1 and pa0 and psur are et to pp_to_be_computed 
     699     !  za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr 
     700    IF ( ( pa0 == 0 .OR. pa1 == 0 .OR. psur == 0 ) & 
     701         .AND. ppdzmin.NE.0 .AND. pphmax.NE.0 ) THEN 
     702        ! 
     703        za1  = (  ppdzmin - pphmax / FLOAT(N-1)  )                                                      & 
     704           & / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(N-1) * (  LOG( COSH( (N - ppkth) / ppacr) )      & 
     705           &                                                   - LOG( COSH( ( 1  - ppkth) / ppacr) )  )  ) 
     706        za0  = ppdzmin - za1 *              TANH( (1-ppkth) / ppacr ) 
     707        zsur =   - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr )  ) 
     708     ELSE 
     709        za1 = pa1 ;       za0 = pa0 ;          zsur = psur 
     710        za2 = pa2                            ! optional (ldbletanh=T) double tanh parameter 
     711     ENDIF 
     712 
     713     ! Reference z-coordinate (depth - scale factor at T- and W-points) 
     714     ! ====================== 
     715     IF( ppkth == 0. ) THEN            !  uniform vertical grid  
     716 
     717        za1 = zhmax / FLOAT(N-1) 
     718 
     719        DO jk = 1, N 
     720           zw = FLOAT( jk ) 
     721           zt = FLOAT( jk ) + 0.5 
     722           gdepw_1d(jk) = ( zw - 1 ) * za1 
     723           gdept_1d(jk) = ( zt - 1 ) * za1 
     724           e3w_1d  (jk) =  za1 
     725           e3t_1d  (jk) =  za1 
     726        END DO 
     727     ELSE                                ! Madec & Imbard 1996 function 
     728        IF( .NOT. ldbletanh ) THEN 
     729           DO jk = 1, N 
     730              zw = REAL( jk ) 
     731              zt = REAL( jk ) + 0.5 
     732              gdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) )  ) 
     733              gdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) )  ) 
     734              e3w_1d  (jk) =          za0      + za1        * TANH(       (zw-zkth) / zacr   ) 
     735              e3t_1d  (jk) =          za0      + za1        * TANH(       (zt-zkth) / zacr   ) 
     736           END DO 
     737        ELSE 
     738           DO jk = 1, N 
     739              zw = FLOAT( jk ) 
     740              zt = FLOAT( jk ) + 0.5 
     741              ! Double tanh function 
     742              gdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth ) / zacr  ) )    & 
     743                 &                             + za2 * zacr2* LOG ( COSH( (zw-zkth2) / zacr2 ) )  ) 
     744              gdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth ) / zacr  ) )    & 
     745                 &                             + za2 * zacr2* LOG ( COSH( (zt-zkth2) / zacr2 ) )  ) 
     746              e3w_1d  (jk) =          za0      + za1        * TANH(       (zw-zkth ) / zacr  )      & 
     747                 &                             + za2        * TANH(       (zw-zkth2) / zacr2 ) 
     748              e3t_1d  (jk) =          za0      + za1        * TANH(       (zt-zkth ) / zacr  )      & 
     749                 &                             + za2        * TANH(       (zt-zkth2) / zacr2 ) 
     750           END DO 
     751        ENDIF 
     752        gdepw_1d(1) = 0.                    ! force first w-level to be exactly at zero 
     753     ENDIF 
     754 
     755     IF ( ln_e3_dep ) THEN      ! e3. = dk[gdep]    
     756        ! 
     757        !==>>>   need to be like this to compute the pressure gradient with ISF.  
     758        !        If not, level beneath the ISF are not aligned (sum(e3t) /= depth) 
     759        !        define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively 
     760        ! 
     761        DO jk = 1, N-1 
     762           e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk) 
     763        END DO 
     764        e3t_1d(N) = e3t_1d(N-1)   ! we don't care because this level is masked in NEMO 
     765 
     766        DO jk = 2, N 
     767           e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1) 
     768        END DO 
     769        e3w_1d(1  ) = 2. * (gdept_1d(1) - gdepw_1d(1)) 
     770     END IF 
     771 
     772     ! 
     773  END SUBROUTINE zgr_z 
     774 
     775 
    476776  !***************************************************** 
    477777  !   function set_child_name(Parentname,Childname) 
     
    501801    !    
    502802  END SUBROUTINE set_child_name 
    503   ! 
    504   !***************************************************** 
    505   !   function set_child_name(Parentname,Childname) 
    506   !***************************************************** 
    507803  ! 
    508804  !***************************************************** 
     
    558854  END SUBROUTINE get_interptype 
    559855  !       
    560   !*****************************************************               
    561   !   end subroutine get_interptype 
    562   !***************************************************** 
    563   !            
    564856  !***************************************************** 
    565857  !   subroutine Init_mask(name,Grid) 
     
    576868    ! 
    577869    IF(jpiglo == 1 .AND. jpjglo == 1) THEN 
    578        CALL Read_Ncdf_var('Bathy_level',name,Grid%Bathy_level) 
     870       CALL read_ncdf_var('Bathy_level',name,Grid%Bathy_level) 
    579871    ELSE 
    580        CALL Read_Ncdf_var('Bathy_level',name,Grid%Bathy_level,(/jpizoom,jpjzoom/),(/jpiglo,jpjglo/) ) 
     872       CALL read_ncdf_var('Bathy_level',name,Grid%Bathy_level,(/jpizoom,jpjzoom/),(/jpiglo,jpjglo/) ) 
    581873    ENDIF 
    582874 
     
    655947  ! 
    656948  !***************************************************** 
    657   !   end subroutine Init_mask 
    658   !***************************************************** 
    659   !  
    660   !***************************************************** 
    661949  !   subroutine Init_Tmask(name,Grid) 
    662950  !***************************************************** 
     
    672960    ! 
    673961    IF(jpiglo == 1 .AND. jpjglo == 1) THEN 
    674        CALL Read_Ncdf_var('Bathy_level',name,Grid%Bathy_level) 
     962       CALL read_ncdf_var('Bathy_level',name,Grid%Bathy_level) 
    675963    ELSE 
    676        CALL Read_Ncdf_var('Bathy_level',name,Grid%Bathy_level,(/jpizoom,jpjzoom/),(/jpiglo,jpjglo/) ) 
     964       CALL read_ncdf_var('Bathy_level',name,Grid%Bathy_level,(/jpizoom,jpjzoom/),(/jpiglo,jpjglo/) ) 
    677965    ENDIF 
    678966    ! 
     
    7121000    TYPE(Coordinates) :: Grid 
    7131001    ! 
    714     CALL Read_Ncdf_var('Bathy_level',filename,Grid%Bathy_level) 
     1002    CALL read_ncdf_var('Bathy_level',filename,Grid%Bathy_level) 
    7151003    !  
    7161004    nx = SIZE(Grid%Bathy_level,1) 
     
    7641052  END SUBROUTINE get_mask 
    7651053  ! 
    766   !***************************************************** 
    767   !   end subroutine get_mask 
    768   !***************************************************** 
    769   !        
    7701054  ! 
    7711055  !***************************************************** 
     
    8031087    ! 
    8041088  END SUBROUTINE write_dimg_var 
    805 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!              
    806 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     1089 
    8071090END MODULE agrif_readwrite 
  • utils/tools/NESTING/src/agrif_types.f90

    r9768 r10025  
    1515  TYPE Coordinates 
    1616     ! 
    17      REAL*8, DIMENSION(:,:), POINTER :: nav_lon,nav_lat => NULL() 
    18      REAL*8, DIMENSION(:,:), POINTER :: glamv, glamu, glamt, glamf => NULL() 
    19      REAL*8, DIMENSION(:,:), POINTER :: gphit, gphiu, gphiv, gphif => NULL() 
    20      REAL*8, DIMENSION(:,:), POINTER :: e1t, e1u, e1v, e1f => NULL() 
    21      REAL*8, DIMENSION(:,:), POINTER :: e2t, e2u, e2v, e2f => NULL() 
    22      REAL*8, DIMENSION(:,:), POINTER :: bathy_level => NULL() 
    23      REAL*8, DIMENSION(:,:), POINTER :: bathy_meter => NULL() 
    24      REAL*8, DIMENSION(:,:), POINTER :: wgt => NULL() 
    25      REAL*8, DIMENSION(:,:,:),POINTER :: fmask,umask,vmask,tmask => NULL() 
    26      REAL*8, DIMENSION(:,:,:),POINTER :: e3t_ps,e3w_ps,gdept_ps,gdepwps => NULL() 
    27      REAL*8, DIMENSION(:,:),POINTER :: gdepw_ps => NULL() 
    28      REAL*8, DIMENSION(:), POINTER :: gdeptht => NULL() 
    29      INTEGER, DIMENSION(:) , POINTER :: time_steps => NULL() 
     17     REAL*8,  DIMENSION(:,:),  POINTER :: nav_lon, nav_lat => NULL() 
     18     REAL*8,  DIMENSION(:,:), POINTER :: glamv, glamu, glamt, glamf => NULL() 
     19     REAL*8,  DIMENSION(:,:), POINTER :: gphit, gphiu, gphiv, gphif => NULL() 
     20     REAL*8,  DIMENSION(:,:), POINTER :: e1t, e1u, e1v, e1f => NULL() 
     21     REAL*8,  DIMENSION(:,:), POINTER :: e2t, e2u, e2v, e2f => NULL() 
     22     REAL*8,  DIMENSION(:,:), POINTER :: bathy_level => NULL() 
     23     REAL*8,  DIMENSION(:,:), POINTER :: bathy_meter => NULL() 
     24     REAL*8,  DIMENSION(:,:), POINTER :: wgt => NULL() 
     25     REAL*8,  DIMENSION(:,:,:),POINTER :: fmask, umask, vmask, tmask => NULL() 
     26     REAL*8,  DIMENSION(:,:,:),POINTER :: e3t_ps, e3w_ps, gdept_ps, gdepwps => NULL() 
     27     REAL*8,  DIMENSION(:,:),  POINTER :: gdepw_ps => NULL() 
     28     REAL*8,  DIMENSION(:),    POINTER :: gdeptht => NULL() 
     29     INTEGER, DIMENSION(:) ,   POINTER :: time_steps => NULL() 
    3030     !      
    3131  END TYPE Coordinates 
     
    3333  ! 
    3434  ! 
    35   CHARACTER*8,DIMENSION(10) :: flxtab = (/'socliot1','socliot2','socliopl',& 
     35  CHARACTER*8,DIMENSION(10) :: flxtab = (/'socliot1','socliot2','socliopl', & 
    3636       'socliocl','socliohu','socliowi','soshfldo','sohefldo','sowaflup','sofbt   '/) 
    3737  ! 
     
    4141  !************************************************************** 
    4242  ! 
    43   INTEGER irafx,irafy 
    44   INTEGER nxfin,nyfin 
    45   INTEGER, PARAMETER :: nbghostcellsfine = 3  
    46   INTEGER, PARAMETER :: nbghostcellscoarse = 1 
    47   !       
    48   INTEGER imin,jmin,imax,jmax,rho,rhot 
    49   INTEGER shlat 
    50   INTEGER N,type_bathy_interp 
     43  INTEGER ::   irafx, irafy 
     44  INTEGER ::   nxfin, nyfin 
     45  !       
     46  INTEGER ::   nbghostcellsfine, imin, jmin, imax, jmax, rho, rhot 
     47  INTEGER ::   shlat 
     48  INTEGER ::   N, type_bathy_interp 
    5149  !  
    52   INTEGER jpizoom,jpjzoom,npt_connect,npt_copy 
    53   !       
    54   REAL*8 rn_hmin 
    55   REAL*8 ppkth2, ppacr2, ppkth,ppacr,ppdzmin,pphmax,smoothing_factor,e3zps_min,e3zps_rat 
    56   REAL*8 psur,pa0,pa1,pa2,adatrj 
     50  INTEGER ::   jpizoom, jpjzoom, npt_connect, npt_copy 
     51  !       
     52  REAL*8 ::   rn_hmin 
     53  REAL*8 ::   ppkth2, ppacr2, ppkth, ppacr, ppdzmin, pphmax, smoothing_factor, e3zps_min, e3zps_rat 
     54  REAL*8 ::   psur, pa0, pa1, pa2, adatrj 
    5755  !        
    58   LOGICAL ldbletanh,ln_e3_dep 
    59   LOGICAL partial_steps,smoothing,bathy_update 
    60   LOGICAL new_topo,removeclosedseas,dimg,iom_activated 
     56  LOGICAL ::   ldbletanh, ln_e3_dep 
     57  LOGICAL ::   partial_steps, smoothing, bathy_update 
     58  LOGICAL ::   new_topo, removeclosedseas, dimg, iom_activated 
     59  LOGICAL ::   ln_agrif_domain 
    6160  !        
    62   CHARACTER*100 parent_meshmask_file,elevation_database,parent_bathy_meter 
    63   CHARACTER*100 elevation_name,parent_batmet_name 
    64   CHARACTER*100 parent_coordinate_file,restart_file,updated_parent_file,restart_trc_file 
    65   CHARACTER*100 dimg_output_file,interp_type 
    66   !       
    67   CHARACTER(len=80),DIMENSION(20) :: flx_Files, u_files, v_files 
    68   CHARACTER(len=255),DIMENSION(20) :: VAR_INTERP 
     61  CHARACTER*100 ::   parent_meshmask_file, elevation_database, parent_bathy_meter, parent_domcfg_output 
     62  CHARACTER*100 ::   elevation_name, parent_batmet_name 
     63  CHARACTER*100 ::   parent_coordinate_file, restart_file, updated_parent_file, updated_parent_domcfg, restart_trc_file 
     64  CHARACTER*100 ::   dimg_output_file, interp_type 
     65  !       
     66  CHARACTER(len=80) , DIMENSION(20) :: flx_Files, u_files, v_files 
     67  CHARACTER(len=255), DIMENSION(20) :: VAR_INTERP 
    6968  ! 
    7069  NAMELIST /input_output/iom_activated 
    7170  ! 
    72   NAMELIST /coarse_grid_files/parent_coordinate_file,parent_meshmask_file 
    73   !       
    74   NAMELIST /bathymetry/new_topo,elevation_database,elevation_name,smoothing,smoothing_factor, & 
    75        npt_connect,npt_copy,removeclosedseas,type_bathy_interp,rn_hmin       
    76   !       
    77   NAMELIST /nesting/imin,imax,jmin,jmax,rho,rhot,bathy_update,updated_parent_file       
    78   ! 
    79   NAMELIST /vertical_grid/ppkth,ppacr,ppdzmin,pphmax,psur,pa0,pa1,N,ldbletanh,ln_e3_dep,pa2,ppkth2,ppacr2 
     71  NAMELIST /coarse_grid_files/parent_coordinate_file, parent_meshmask_file, parent_domcfg_output 
     72  !       
     73  NAMELIST /bathymetry/new_topo, elevation_database, elevation_name, smoothing, smoothing_factor, & 
     74                       ln_agrif_domain, npt_connect, npt_copy, removeclosedseas, type_bathy_interp, rn_hmin       
     75  !       
     76  NAMELIST /nesting/nbghostcellsfine, imin, imax, jmin, jmax, rho, rhot, bathy_update, updated_parent_file, updated_parent_domcfg       
     77  ! 
     78  NAMELIST /vertical_grid/ppkth, ppacr, ppdzmin, pphmax, psur, pa0, pa1, N, ldbletanh, ln_e3_dep, pa2, ppkth2, ppacr2 
    8079  !  
    81   NAMELIST /partial_cells/partial_steps,parent_bathy_meter,parent_batmet_name,e3zps_min,e3zps_rat       
    82   ! 
    83   NAMELIST /nemo_coarse_grid/ jpizoom,jpjzoom  
     80  NAMELIST /partial_cells/partial_steps, parent_bathy_meter, parent_batmet_name, e3zps_min, e3zps_rat       
     81  ! 
     82  NAMELIST /nemo_coarse_grid/ jpizoom, jpjzoom  
    8483  !           
    85   NAMELIST /forcing_files/ flx_files, u_files, v_files  
     84  NAMELIST /forcing_files/ flx_files,  u_files, v_files  
    8685  !            
    8786  NAMELIST /interp/ VAR_INTERP 
    8887  !       
    89   NAMELIST /restart/ restart_file,shlat,dimg,dimg_output_file,adatrj,interp_type  
    90  
    91   NAMELIST /restart_trc/ restart_trc_file,interp_type  
     88  NAMELIST /restart/ restart_file, shlat, dimg, dimg_output_file, adatrj, interp_type  
     89 
     90  NAMELIST /restart_trc/ restart_trc_file, interp_type  
    9291  ! 
    9392CONTAINS 
     
    102101  !******************************************************** 
    103102  !        
    104   SUBROUTINE agrif_grid_allocate(Grid,nx,ny) 
     103  SUBROUTINE agrif_grid_allocate(Grid, nx, ny) 
    105104    ! 
    106105    TYPE(Coordinates) :: Grid 
    107     INTEGER :: nx,ny 
     106    INTEGER :: nx, ny 
    108107    ! 
    109108    ALLOCATE(Grid%nav_lon(nx,ny),Grid%nav_lat(nx,ny)) 
     
    175174       jmax = jmax + jpjzoom - 1 
    176175       ! 
    177        nxfin = (imax-imin)*irafx+2*nbghostcellsfine+2 
    178        nyfin = (jmax-jmin)*irafy+2*nbghostcellsfine+2 
     176       IF( ln_agrif_domain ) THEN 
     177          nxfin = (imax-imin)*irafx+2*nbghostcellsfine+2 
     178          nyfin = (jmax-jmin)*irafy+2*nbghostcellsfine+2 
     179       ELSE 
     180          nbghostcellsfine = 0 
     181          nxfin = (imax-imin+1)*irafx 
     182          nyfin = (jmax-jmin+1)*irafy 
     183       ENDIF 
    179184       ! 
    180185    ELSE 
    181186       ! 
    182187       PRINT *,'namelist file ''',TRIM(namelistname),''' not found' 
    183        STOP   
     188       STOP  
    184189       ! 
    185190    END IF 
  • utils/tools/NESTING/src/io_netcdf.f90

    r2143 r10025  
    2525  USE agrif_types 
    2626  !       
    27   INTERFACE Read_Ncdf_var 
    28      MODULE PROCEDURE Read_Ncdf_var1d_Real,   & 
    29           Read_Ncdf_var2d_Real,   & 
    30           Read_Ncdf_var2d_Real_bis,   & 
    31           Read_Ncdf_var3d_Real,   & 
    32           Read_Ncdf_var4d_Real,   & 
    33           Read_Ncdf_var3d_Real_t, & 
    34           Read_Ncdf_var4d_Real_t, & 
    35           Read_Ncdf_var4d_Real_nt, & 
    36           Read_Ncdf_var1d_Int,    & 
    37           Read_Ncdf_var2d_Int,    & 
    38           Read_Ncdf_var3d_Int,    & 
    39           Read_Ncdf_var4d_Int,    & 
    40           Read_Ncdf_var0d_Int,    & 
    41           Read_Ncdf_var0d_Real 
     27  INTERFACE read_ncdf_var 
     28     MODULE PROCEDURE  & 
     29        read_ncdf_var0d_real, read_ncdf_var1d_real, read_ncdf_var2d_real  , read_ncdf_var2d_real_bis,                        & 
     30        read_ncdf_var3d_real, read_ncdf_var4d_real, read_ncdf_var3d_real_t, read_ncdf_var4d_real_t, read_ncdf_var4d_real_nt, & 
     31        read_ncdf_var0d_int,  read_ncdf_var1d_int , read_ncdf_var2d_int   , read_ncdf_var3d_int   , read_ncdf_var4d_int 
    4232  END INTERFACE 
    4333  ! 
    44   INTERFACE Write_Ncdf_var 
    45      MODULE PROCEDURE Write_Ncdf_var1d_Real,   & 
    46           Write_Ncdf_var2d_Real,   & 
    47           Write_Ncdf_var3d_Real,   & 
    48           Write_Ncdf_var4d_Real,   & 
    49           Write_Ncdf_var3d_Real_t, & 
    50           Write_Ncdf_var4d_Real_t, & 
    51           Write_Ncdf_var4d_Real_nt, &                            
    52           Write_Ncdf_var2d_Real_bis,&            
    53           Write_Ncdf_var1d_Int,    & 
    54           Write_Ncdf_var2d_Int,    & 
    55           Write_Ncdf_var3d_Int,    & 
    56           Write_Ncdf_var4d_Int,    & 
    57           Write_Ncdf_var0d_Real 
     34  INTERFACE write_ncdf_var 
     35     MODULE PROCEDURE & 
     36        write_ncdf_var0d_real, write_ncdf_var1d_real  , write_ncdf_var2d_real  , write_ncdf_var3d_real,    & 
     37        write_ncdf_var4d_real, write_ncdf_var3d_real_t, write_ncdf_var4d_real_t, write_ncdf_var4d_real_nt, &                            
     38        write_ncdf_var2d_real_bis ,                                                                        &             
     39        write_ncdf_var0d_int, write_ncdf_var1d_int, write_ncdf_var2d_int, write_ncdf_var3d_int, write_ncdf_var4d_int 
    5840  END INTERFACE 
    5941  !  
    60   INTERFACE Copy_Ncdf_att 
    61      MODULE PROCEDURE Copy_Ncdf_att_latlon,Copy_Ncdf_att_var 
     42  INTERFACE copy_ncdf_att 
     43     MODULE PROCEDURE copy_ncdf_att_latlon,copy_ncdf_att_var 
    6244  END INTERFACE 
    6345  ! 
     
    6547  ! 
    6648  !**************************************************************** 
    67   !   subroutine Read_Ncdf_dim               * 
     49  !   subroutine read_ncdf_dim               * 
    6850  !                        * 
    6951  ! subroutine to retrieve value of a given dimension       * 
     
    7557  !**************************************************************** 
    7658  ! 
    77   SUBROUTINE Read_Ncdf_dim(dimname,file,dimval)    
     59  SUBROUTINE read_ncdf_dim(dimname,file,dimval)    
    7860    !       
    7961    IMPLICIT NONE 
     
    9779    status = nf90_close(ncid) 
    9880    !       
    99   END SUBROUTINE Read_Ncdf_dim 
     81  END SUBROUTINE read_ncdf_dim 
    10082  ! 
    10183  !************************************************************** 
    102   ! end subroutine Read_Ncdf_dim 
     84  ! end subroutine read_ncdf_dim 
    10385  !**************************************************************       
    10486  ! 
    10587  !**************************************************************** 
    106   !   subroutine Write_Ncdf_dim              * 
     88  !   subroutine write_ncdf_dim              * 
    10789  !                        * 
    10890  ! subroutine to write a dimension in a given file      * 
     
    11496  !**************************************************************** 
    11597  ! 
    116   SUBROUTINE Write_Ncdf_dim(dimname,file,dimval)    
     98  SUBROUTINE write_ncdf_dim(dimname,file,dimval)    
    11799    !       
    118100    IMPLICIT NONE 
     
    143125    status = nf90_close(ncid) 
    144126    !       
    145   END SUBROUTINE Write_Ncdf_dim 
     127  END SUBROUTINE write_ncdf_dim 
    146128  ! 
    147129  !************************************************************** 
    148   ! end subroutine Write_Ncdf_dim 
     130  ! end subroutine write_ncdf_dim 
    149131  !**************************************************************             
    150132  !  
    151133  !**************************************************************** 
    152   !   subroutine Read_Ncdf_var               * 
     134  !   subroutine read_ncdf_var               * 
    153135  !                        * 
    154136  ! subroutine to retrieve values of a given variable       * 
     
    159141  !                        * 
    160142  !**************************************************************** 
    161   !           
    162   !************************************************************** 
    163   !   subroutine Read_Ncdf_var1d_real 
    164   !**************************************************************           
    165   !       
    166   SUBROUTINE Read_Ncdf_var1d_Real(varname,file,tabvar) 
     143  !       
     144  SUBROUTINE read_ncdf_var1d_real(varname,file,tabvar) 
    167145    !       
    168146    IMPLICIT NONE 
     
    202180    status = nf90_close(ncid) 
    203181    !  
    204   END SUBROUTINE Read_Ncdf_var1d_Real 
     182  END SUBROUTINE read_ncdf_var1d_real 
    205183  !            
    206   !************************************************************** 
    207   ! end subroutine Read_Ncdf_var1d_real 
    208   !**************************************************************           
    209   !             
    210   !            
    211   !************************************************************** 
    212   !   subroutine Read_Ncdf_var2d_real 
    213   !**************************************************************           
    214   !       
    215   SUBROUTINE Read_Ncdf_var2d_Real(varname,file,tabvar) 
     184  !       
     185  SUBROUTINE read_ncdf_var2d_real(varname,file,tabvar) 
    216186    !       
    217187    IMPLICIT NONE 
     
    250220    status = nf90_close(ncid) 
    251221    !      
    252   END SUBROUTINE Read_Ncdf_var2d_Real 
     222  END SUBROUTINE read_ncdf_var2d_real 
    253223  !            
    254   !************************************************************** 
    255   ! end subroutine Read_Ncdf_var2d_real 
    256   !**************************************************************           
    257   !     
    258  
    259   !            
    260   !************************************************************** 
    261   !   subroutine Read_Ncdf_var2d_real_bis 
    262   !**************************************************************           
    263   !       
    264   SUBROUTINE Read_Ncdf_var2d_Real_bis(varname,file,tabvar,strt,cnt) 
     224  !       
     225  SUBROUTINE read_ncdf_var2d_real_bis(varname,file,tabvar,strt,cnt) 
    265226    !       
    266227    IMPLICIT NONE 
     
    301262    status = nf90_close(ncid) 
    302263    !      
    303   END SUBROUTINE Read_Ncdf_var2d_Real_bis 
     264  END SUBROUTINE read_ncdf_var2d_real_bis 
    304265  !            
    305   !************************************************************** 
    306   ! end subroutine Read_Ncdf_var2d_real_bis 
    307   !**************************************************************           
    308   !   
    309  
    310  
    311   !            
    312   !************************************************************** 
    313   !   subroutine Read_Ncdf_var3d_real 
    314   !**************************************************************           
    315   !       
    316   SUBROUTINE Read_Ncdf_var3d_Real(varname,file,tabvar) 
     266  !       
     267  SUBROUTINE read_ncdf_var3d_real(varname,file,tabvar) 
    317268    !       
    318269    IMPLICIT NONE 
     
    358309    status = nf90_close(ncid) 
    359310    !      
    360   END SUBROUTINE Read_Ncdf_var3d_Real 
     311  END SUBROUTINE read_ncdf_var3d_real 
    361312  !            
    362   !************************************************************** 
    363   ! end subroutine Read_Ncdf_var3d_real 
    364   !**************************************************************           
    365   !  
    366   !            
    367   !************************************************************** 
    368   !   subroutine Read_Ncdf_var4d_real 
    369   !**************************************************************           
    370   !       
    371   SUBROUTINE Read_Ncdf_var4d_Real(varname,file,tabvar) 
     313  !       
     314  SUBROUTINE read_ncdf_var4d_real(varname,file,tabvar) 
    372315    !       
    373316    IMPLICIT NONE 
     
    410353    status = nf90_close(ncid) 
    411354    !      
    412   END SUBROUTINE Read_Ncdf_var4d_Real 
    413  
    414   SUBROUTINE Read_Ncdf_var0d_Real(varname,file,tabvar) 
     355  END SUBROUTINE read_ncdf_var4d_real 
     356 
     357  SUBROUTINE read_ncdf_var0d_real(varname,file,tabvar) 
    415358    !       
    416359    IMPLICIT NONE 
     
    438381    status = nf90_close(ncid) 
    439382    !      
    440   END SUBROUTINE Read_Ncdf_var0d_Real 
    441  
    442   SUBROUTINE Read_Ncdf_var0d_Int(varname,file,tabvar) 
     383  END SUBROUTINE read_ncdf_var0d_real 
     384 
     385  SUBROUTINE read_ncdf_var0d_int(varname,file,tabvar) 
    443386    !       
    444387    IMPLICIT NONE 
     
    466409    status = nf90_close(ncid) 
    467410    !      
    468   END SUBROUTINE Read_Ncdf_var0d_Int 
     411  END SUBROUTINE read_ncdf_var0d_int 
    469412  !            
    470   !************************************************************** 
    471   ! end subroutine Read_Ncdf_var4d_real 
    472   !**************************************************************           
    473   !              
    474   !            
    475   !************************************************************** 
    476   !   subroutine Read_Ncdf_var1d_int 
    477   !**************************************************************           
    478   !       
    479   SUBROUTINE Read_Ncdf_var1d_Int(varname,file,tabvar) 
     413  !       
     414  SUBROUTINE read_ncdf_var1d_int(varname,file,tabvar) 
    480415    !       
    481416    IMPLICIT NONE 
     
    515450    status = nf90_close(ncid) 
    516451    !       
    517   END SUBROUTINE Read_Ncdf_var1d_Int 
     452  END SUBROUTINE read_ncdf_var1d_int 
    518453  !            
    519   !************************************************************** 
    520   ! end subroutine Read_Ncdf_var1d_int 
    521   !**************************************************************           
    522   !             
    523   !            
    524   !************************************************************** 
    525   !   subroutine Read_Ncdf_var2d_int 
    526   !**************************************************************           
    527   !       
    528   SUBROUTINE Read_Ncdf_var2d_Int(varname,file,tabvar) 
     454  !       
     455  SUBROUTINE read_ncdf_var2d_int(varname,file,tabvar) 
    529456    !       
    530457    IMPLICIT NONE 
     
    563490    status = nf90_close(ncid) 
    564491    !       
    565   END SUBROUTINE Read_Ncdf_var2d_Int 
     492  END SUBROUTINE read_ncdf_var2d_int 
    566493  !            
    567   !************************************************************** 
    568   ! end subroutine Read_Ncdf_var2d_Int 
    569   !**************************************************************           
    570   !     
    571  
    572   !            
    573   !************************************************************** 
    574   !   subroutine Read_Ncdf_var3d_Int 
    575   !**************************************************************           
    576   !       
    577   SUBROUTINE Read_Ncdf_var3d_Int(varname,file,tabvar) 
     494  !       
     495  SUBROUTINE read_ncdf_var3d_int(varname,file,tabvar) 
    578496    !       
    579497    IMPLICIT NONE 
     
    615533    status = nf90_close(ncid) 
    616534    !       
    617   END SUBROUTINE Read_Ncdf_var3d_Int 
     535  END SUBROUTINE read_ncdf_var3d_int 
    618536  !            
    619   !************************************************************** 
    620   ! end subroutine Read_Ncdf_var3d_int 
    621   !**************************************************************   
    622   ! 
    623   ! 
    624   ! 
    625   !************************************************************** 
    626   !   subroutine Read_Ncdf_var4d_int 
    627   !**************************************************************           
    628   !       
    629   SUBROUTINE Read_Ncdf_var4d_Int(varname,file,tabvar) 
     537  !       
     538  SUBROUTINE read_ncdf_var4d_int(varname,file,tabvar) 
    630539    !       
    631540    IMPLICIT NONE 
     
    668577    status = nf90_close(ncid) 
    669578    !      
    670   END SUBROUTINE Read_Ncdf_var4d_Int 
     579  END SUBROUTINE read_ncdf_var4d_int 
    671580  !            
    672   !************************************************************** 
    673   ! end subroutine Read_Ncdf_var4d_real 
    674   !**************************************************************          
    675   !     
    676   ! 
    677   !**************************************************************** 
    678   !   subroutine Write_Ncdf_var              * 
     581  ! 
     582  !**************************************************************** 
     583  !   subroutine write_ncdf_var              * 
    679584  !                        * 
    680585  ! subroutine to write a variable in a given file       * 
     
    688593  !     
    689594  !       
    690   !************************************************************** 
    691   !   subroutine Write_Ncdf_var1d_real 
    692   !**************************************************************     
    693   !       
    694   SUBROUTINE Write_Ncdf_var1d_Real(varname,dimname,file,tabvar,typevar) 
     595  SUBROUTINE write_ncdf_var1d_real(varname,dimname,file,tabvar,typevar) 
    695596    !       
    696597    IMPLICIT NONE 
    697598    !        
    698599    CHARACTER(*),INTENT(in) :: varname,file,dimname,typevar 
    699     REAL*8, DIMENSION(:), POINTER :: tabvar 
     600    REAL*8, DIMENSION(:), INTENT(in) :: tabvar 
    700601    ! 
    701602    ! local variables 
     
    711612    ENDIF 
    712613    !      
    713     status = nf90_inq_dimid(ncid,dimname, dimid) 
     614    status = nf90_inq_dimid(ncid,dimname,dimid) 
    714615    status = nf90_inq_varid(ncid,varname,varid) 
    715616    status = nf90_redef(ncid) 
     
    725626    status = nf90_close(ncid) 
    726627    !  
    727   END SUBROUTINE Write_Ncdf_var1d_Real 
    728   !       
    729   !************************************************************** 
    730   !   end subroutine Write_Ncdf_var1d_real 
    731   !************************************************************** 
    732   !       
    733   !************************************************************** 
    734   !   subroutine Write_Ncdf_var2d_real 
    735   !**************************************************************     
    736   !       
    737   SUBROUTINE Write_Ncdf_var2d_Real_bis(varname,dimname,file,tabvar,nbdim,typevar) 
     628  END SUBROUTINE write_ncdf_var1d_real 
     629  !       
     630  !       
     631  SUBROUTINE write_ncdf_var2d_real_bis(varname,dimname,file,tabvar,nbdim,typevar) 
    738632    !       
    739633    IMPLICIT NONE 
     
    743637    CHARACTER(*), DIMENSION(4) :: dimname 
    744638    REAL*8, DIMENSION(:,:) :: tabvar 
    745     REAL*8, DIMENSION(:,:,:),POINTER :: tabtemp3d 
    746     REAL*8, DIMENSION(:,:,:,:),POINTER :: tabtemp4d 
     639    REAL*8, ALLOCATABLE, DIMENSION(:,:,:) :: tabtemp3d 
     640    REAL*8, ALLOCATABLE, DIMENSION(:,:,:,:) :: tabtemp4d 
    747641    ! 
    748642    ! local variables 
     
    796690    IF(nbdim==3) status = nf90_put_var(ncid,varid,tabtemp3d) 
    797691    !       
    798     IF(ASSOCIATED( tabtemp3d ) ) DEALLOCATE( tabtemp3d )    
    799     IF(ASSOCIATED( tabtemp4d ) ) DEALLOCATE( tabtemp4d )      
    800     ! 
    801     status = nf90_close(ncid) 
    802     !  
    803   END SUBROUTINE Write_Ncdf_var2d_Real_bis 
    804   !       
    805   !************************************************************** 
    806   !   end subroutine Write_Ncdf_var2d_real 
    807   !**************************************************************   
    808   !    
    809   !       
    810   !************************************************************** 
    811   !   subroutine Write_Ncdf_var2d_real 
    812   !**************************************************************     
    813   !       
    814   SUBROUTINE Write_Ncdf_var2d_Real(varname,dimname,file,tabvar,typevar) 
     692    IF(ALLOCATED( tabtemp3d ) ) DEALLOCATE( tabtemp3d )    
     693    IF(ALLOCATED( tabtemp4d ) ) DEALLOCATE( tabtemp4d )      
     694    ! 
     695    status = nf90_close(ncid) 
     696    !  
     697  END SUBROUTINE write_ncdf_var2d_real_bis 
     698  !       
     699  !       
     700  SUBROUTINE write_ncdf_var2d_real(varname,dimname,file,tabvar,typevar) 
    815701    !       
    816702    !      implicit none 
     
    818704    CHARACTER(*),INTENT(in) :: varname,file,typevar 
    819705    CHARACTER(*), DIMENSION(2) :: dimname 
    820     REAL*8, DIMENSION(:,:), POINTER :: tabvar 
     706    REAL*8, DIMENSION(:,:), INTENT(in) :: tabvar 
    821707    ! 
    822708    ! local variables 
     
    851737    status = nf90_close(ncid) 
    852738    !  
    853   END SUBROUTINE Write_Ncdf_var2d_Real 
    854   !       
    855   !************************************************************** 
    856   !   end subroutine Write_Ncdf_var2d_real 
    857   !**************************************************************   
    858   ! 
    859   !       
    860   !************************************************************** 
    861   !   subroutine Write_Ncdf_var3d_real 
    862   !************************************************************** 
    863   !       
    864   SUBROUTINE Write_Ncdf_var3d_Real(varname,dimname,file,tabvar,typevar) 
     739  END SUBROUTINE write_ncdf_var2d_real 
     740  !       
     741  !       
     742  SUBROUTINE write_ncdf_var3d_real(varname,dimname,file,tabvar,typevar) 
    865743    !       
    866744    IMPLICIT NONE 
     
    868746    CHARACTER(*),INTENT(in) :: varname,file,typevar 
    869747    CHARACTER(*),DIMENSION(3),INTENT(in) :: dimname 
    870     REAL*8, DIMENSION(:,:,:), POINTER :: tabvar 
     748    REAL*8, DIMENSION(:,:,:), INTENT(in) :: tabvar 
    871749    ! 
    872750    ! local variables 
     
    902780    status = nf90_close(ncid) 
    903781    !  
    904   END SUBROUTINE Write_Ncdf_var3d_Real 
    905   !       
    906   !************************************************************** 
    907   !   end subroutine Write_Ncdf_var3d_real 
    908   !************************************************************** 
    909   !  
    910   ! 
    911   !************************************************************** 
    912   !   subroutine Write_Ncdf_var4d_real 
    913   !************************************************************** 
    914   !       
    915   SUBROUTINE Write_Ncdf_var4d_Real(varname,dimname,file,tabvar,typevar) 
     782  END SUBROUTINE write_ncdf_var3d_real 
     783  !       
     784  !       
     785  SUBROUTINE write_ncdf_var4d_real(varname,dimname,file,tabvar,typevar) 
    916786    !       
    917787    IMPLICIT NONE 
     
    919789    CHARACTER(*),INTENT(in) :: varname,file,typevar 
    920790    CHARACTER(*),DIMENSION(4),INTENT(in) :: dimname 
    921     REAL*8, DIMENSION(:,:,:,:), POINTER :: tabvar 
     791    REAL*8, DIMENSION(:,:,:,:), INTENT(in) :: tabvar 
    922792    ! 
    923793    ! local variables 
     
    959829    status = nf90_close(ncid) 
    960830    !  
    961   END SUBROUTINE Write_Ncdf_var4d_Real 
    962   !       
    963   !************************************************************** 
    964   !   end subroutine Write_Ncdf_var4d_real 
    965   !************************************************************** 
    966   ! 
    967   ! 
    968   !   
    969   !************************************************************** 
    970   !   subroutine Write_Ncdf_var1d_Int 
    971   !**************************************************************     
    972   !       
    973   SUBROUTINE Write_Ncdf_var1d_Int(varname,dimname,file,tabvar) 
    974     !       
    975     IMPLICIT NONE 
    976     !        
    977     CHARACTER(*),INTENT(in) :: varname,file,dimname 
    978     INTEGER, DIMENSION(:), POINTER :: tabvar 
     831  END SUBROUTINE write_ncdf_var4d_real 
     832  !       
     833  !       
     834  SUBROUTINE write_ncdf_var1d_int(varname,dimname,file,tabvar,typevar) 
     835    !       
     836    IMPLICIT NONE 
     837    !        
     838    CHARACTER(*),INTENT(in) :: varname,file,dimname,typevar 
     839    INTEGER, DIMENSION(:), INTENT(in) :: tabvar 
    979840    ! 
    980841    ! local variables 
     
    1000861    status = nf90_close(ncid) 
    1001862    !  
    1002   END SUBROUTINE Write_Ncdf_var1d_Int 
    1003   !       
    1004   !************************************************************** 
    1005   !   end subroutine Write_Ncdf_var1d_Int 
    1006   !************************************************************** 
    1007   !   
    1008   !       
    1009   !************************************************************** 
    1010   !   subroutine Write_Ncdf_var2d_Int 
    1011   !**************************************************************     
    1012   !       
    1013   SUBROUTINE Write_Ncdf_var2d_Int(varname,dimname,file,tabvar) 
    1014     !       
    1015     IMPLICIT NONE 
    1016     !        
    1017     CHARACTER(*),INTENT(in) :: varname,file 
    1018     CHARACTER(*), DIMENSION(2) :: dimname 
    1019     INTEGER, DIMENSION(:,:), POINTER :: tabvar 
     863  END SUBROUTINE write_ncdf_var1d_int 
     864  !       
     865  !       
     866  SUBROUTINE write_ncdf_var2d_int(varname,dimname,file,tabvar,typevar) 
     867    !       
     868    IMPLICIT NONE 
     869    !        
     870    CHARACTER(*), INTENT(in) :: varname,file,typevar 
     871    CHARACTER(*), DIMENSION(2), INTENT(in) :: dimname 
     872    INTEGER, DIMENSION(:,:), INTENT(in) :: tabvar 
    1020873    ! 
    1021874    ! local variables 
     
    1031884    ENDIF 
    1032885    !      
    1033     status = nf90_inq_dimid(ncid,dimname(1), dimid1) 
    1034     status = nf90_inq_dimid(ncid,dimname(2), dimid2) 
     886    status = nf90_inq_dimid(ncid,dimname(1),dimid1) 
     887    status = nf90_inq_dimid(ncid,dimname(2),dimid2) 
    1035888    status = nf90_inq_varid(ncid,varname,varid) 
    1036889    status = nf90_redef(ncid) 
    1037     status = nf90_def_var(ncid,varname,nf90_int,     & 
    1038          (/dimid1,dimid2/),varid) 
     890    status = nf90_def_var(ncid,varname,nf90_int,(/dimid1,dimid2/),varid) 
    1039891    status = nf90_enddef(ncid) 
    1040892    status = nf90_put_var(ncid,varid,tabvar)      
     
    1042894    status = nf90_close(ncid) 
    1043895    !  
    1044   END SUBROUTINE Write_Ncdf_var2d_Int 
    1045   !       
    1046   !************************************************************** 
    1047   !   end subroutine Write_Ncdf_var2d_Int 
    1048   !**************************************************************   
    1049   !  
    1050   !       
    1051   !************************************************************** 
    1052   !   subroutine Write_Ncdf_var3d_Int 
    1053   !************************************************************** 
    1054   !       
    1055   SUBROUTINE Write_Ncdf_var3d_Int(varname,dimname,file,tabvar) 
    1056     !       
    1057     IMPLICIT NONE 
    1058     !        
    1059     CHARACTER(*),INTENT(in) :: varname,file 
     896  END SUBROUTINE write_ncdf_var2d_int 
     897  !       
     898  !       
     899  SUBROUTINE write_ncdf_var3d_int(varname,dimname,file,tabvar,typevar) 
     900    !       
     901    IMPLICIT NONE 
     902    !        
     903    CHARACTER(*),INTENT(in) :: varname,file,typevar 
    1060904    CHARACTER(*),DIMENSION(3),INTENT(in) :: dimname 
    1061     INTEGER, DIMENSION(:,:,:), POINTER :: tabvar 
     905    INTEGER, DIMENSION(:,:,:), INTENT(in) :: tabvar 
    1062906    ! 
    1063907    ! local variables 
     
    1085929    status = nf90_close(ncid) 
    1086930    !  
    1087   END SUBROUTINE Write_Ncdf_var3d_Int 
    1088   !       
    1089   !************************************************************** 
    1090   !   end subroutine Write_Ncdf_var3d_Int 
    1091   !************************************************************** 
    1092   !  
    1093   ! 
    1094   !************************************************************** 
    1095   !   subroutine Write_Ncdf_var4d_Int 
    1096   !************************************************************** 
    1097   !       
    1098   SUBROUTINE Write_Ncdf_var4d_Int(varname,dimname,file,tabvar) 
    1099     !       
    1100     IMPLICIT NONE 
    1101     !        
    1102     CHARACTER(*),INTENT(in) :: varname,file 
     931  END SUBROUTINE write_ncdf_var3d_int 
     932  !       
     933  !       
     934  SUBROUTINE write_ncdf_var4d_int(varname,dimname,file,tabvar,typevar) 
     935    !       
     936    IMPLICIT NONE 
     937    !        
     938    CHARACTER(*),INTENT(in) :: varname,file,typevar 
    1103939    CHARACTER(*),DIMENSION(4),INTENT(in) :: dimname 
    1104     INTEGER, DIMENSION(:,:,:,:), POINTER :: tabvar 
     940    INTEGER, DIMENSION(:,:,:,:), INTENT(in) :: tabvar 
    1105941    ! 
    1106942    ! local variables 
     
    1129965    status = nf90_close(ncid) 
    1130966    !  
    1131   END SUBROUTINE Write_Ncdf_var4d_Int 
    1132   !       
    1133   !************************************************************** 
    1134   !   end subroutine Write_Ncdf_var4d_Int 
    1135   !**************************************************************  
    1136   ! 
    1137   ! 
    1138   !**************************************************************** 
    1139   !   subroutine Read_Ncdf_var_t             * 
     967  END SUBROUTINE write_ncdf_var4d_int 
     968  !       
     969  ! 
     970  !**************************************************************** 
     971  !   subroutine read_ncdf_var_t             * 
    1140972  !                        * 
    1141973  ! subroutine to read a variable in a given file for time t   * 
     
    1148980  !**************************************************************** 
    1149981  ! 
    1150   ! 
    1151   !************************************************************** 
    1152   !   subroutine Read_Ncdf_var3d_real_t 
    1153   !**************************************************************           
    1154   !       
    1155   SUBROUTINE Read_Ncdf_var3d_Real_t(varname,file,tabvar,time) 
     982  !       
     983  SUBROUTINE read_ncdf_var3d_real_t(varname,file,tabvar,time) 
    1156984    !       
    1157985    USE agrif_types      
     
    12011029    status = nf90_close(ncid) 
    12021030    !      
    1203   END SUBROUTINE Read_Ncdf_var3d_Real_t 
     1031  END SUBROUTINE read_ncdf_var3d_real_t 
    12041032  !            
    1205   !************************************************************** 
    1206   ! end subroutine Read_Ncdf_var3d_real_t 
    1207   !**************************************************************           
    1208   ! 
    1209   ! 
    1210   ! 
    1211   !**************************************************************** 
    1212   !   subroutine Write_Ncdf_var_t               * 
     1033  ! 
     1034  !**************************************************************** 
     1035  !   subroutine write_ncdf_var_t               * 
    12131036  !                        * 
    12141037  ! subroutine to write a variable in a given file for time t  * 
     
    12231046  ! 
    12241047  !       
    1225   !************************************************************** 
    1226   !   subroutine Write_Ncdf_var3d_real_t 
    1227   !************************************************************** 
    1228   !       
    1229   SUBROUTINE Write_Ncdf_var3d_Real_t(varname,dimname,file,tabvar,time,typevar) 
     1048  SUBROUTINE write_ncdf_var3d_real_t(varname,dimname,file,tabvar,time,typevar) 
    12301049    !       
    12311050    IMPLICIT NONE 
     
    12341053    CHARACTER(*),DIMENSION(3),INTENT(in) :: dimname 
    12351054    INTEGER :: time 
    1236     REAL*8, DIMENSION(:,:,:), POINTER :: tabvar 
     1055    REAL*8, DIMENSION(:,:,:), INTENT(in) :: tabvar 
    12371056    ! 
    12381057    ! local variables 
     
    12801099    status = nf90_close(ncid) 
    12811100    !  
    1282   END SUBROUTINE Write_Ncdf_var3d_Real_t 
    1283   !       
    1284   !************************************************************** 
    1285   !   end subroutine Write_Ncdf_var3d_real_t 
    1286   !************************************************************** 
    1287   ! 
    1288   ! 
    1289   !**************************************************************** 
    1290   !   subroutine Read_Ncdf_var_t             * 
     1101  END SUBROUTINE write_ncdf_var3d_real_t 
     1102  !       
     1103  ! 
     1104  !**************************************************************** 
     1105  !   subroutine read_ncdf_var_t             * 
    12911106  !                        * 
    12921107  ! subroutine to read a variable in a given file for time t   * 
     
    13001115  !**************************************************************** 
    13011116  ! 
    1302   ! 
    1303   !************************************************************** 
    1304   !   subroutine Read_Ncdf_var4d_real_nt 
    1305   !**************************************************************           
    1306   !       
    1307   SUBROUTINE Read_Ncdf_var4d_Real_nt(varname,file,tabvar,time,level) 
     1117  !       
     1118  SUBROUTINE read_ncdf_var4d_real_nt(varname,file,tabvar,time,level) 
    13081119    !       
    13091120    USE agrif_types      
     
    13531164    status = nf90_close(ncid) 
    13541165    !      
    1355   END SUBROUTINE Read_Ncdf_var4d_Real_nt 
     1166  END SUBROUTINE read_ncdf_var4d_real_nt 
    13561167  !            
    1357   !************************************************************** 
    1358   ! end subroutine Read_Ncdf_var4d_real_nt 
    1359   !**************************************************************           
    1360   ! 
    1361   ! 
    1362   !************************************************************** 
    1363   !   subroutine Read_Ncdf_var4d_real_t 
    1364   !**************************************************************           
    1365   !       
    1366   SUBROUTINE Read_Ncdf_var4d_Real_t(varname,file,tabvar,time) 
     1168  !       
     1169  SUBROUTINE read_ncdf_var4d_real_t(varname,file,tabvar,time) 
    13671170    !       
    13681171    USE agrif_types      
     
    14051208    status = nf90_close(ncid) 
    14061209    !      
    1407   END SUBROUTINE Read_Ncdf_var4d_Real_t 
     1210  END SUBROUTINE read_ncdf_var4d_real_t 
    14081211  !            
    1409   !************************************************************** 
    1410   ! end subroutine Read_Ncdf_var4d_real_t 
    1411   !**************************************************************           
    1412   ! 
    1413   ! 
    1414   !**************************************************************** 
    1415   !   subroutine Write_Ncdf_var_t               * 
     1212  !**************************************************************** 
     1213  !   subroutine write_ncdf_var_t               * 
    14161214  !                        * 
    14171215  ! subroutine to write a variable in a given file for time t  * 
     
    14271225  ! 
    14281226  !       
    1429   !************************************************************** 
    1430   !   subroutine Write_Ncdf_var4d_real_t 
    1431   !************************************************************** 
    1432   !       
    1433   SUBROUTINE Write_Ncdf_var4d_Real_t(varname,dimname,file,tabvar,time,typevar) 
     1227  SUBROUTINE write_ncdf_var4d_real_t(varname,dimname,file,tabvar,time,typevar) 
    14341228    !       
    14351229    IMPLICIT NONE 
     
    14381232    CHARACTER(*),DIMENSION(4),INTENT(in) :: dimname 
    14391233    INTEGER :: time,level 
    1440     REAL*8, DIMENSION(:,:,:,:), POINTER :: tabvar 
     1234    REAL*8, DIMENSION(:,:,:,:), INTENT(in) :: tabvar 
    14411235    ! 
    14421236    ! local variables 
     
    14841278    status = nf90_close(ncid) 
    14851279    !  
    1486   END SUBROUTINE Write_Ncdf_var4d_Real_t 
    1487   !       
    1488   !************************************************************** 
    1489   !   end subroutine Write_Ncdf_var4d_real_t 
    1490   !************************************************************** 
    1491   ! 
    1492   !       
    1493   !************************************************************** 
    1494   !   subroutine Write_Ncdf_var4d_real_nt 
    1495   !************************************************************** 
    1496   !       
    1497   SUBROUTINE Write_Ncdf_var4d_Real_nt(varname,dimname,file,tabvar,time,level,typevar) 
     1280  END SUBROUTINE write_ncdf_var4d_real_t 
     1281  !       
     1282  !       
     1283  SUBROUTINE write_ncdf_var4d_real_nt(varname,dimname,file,tabvar,time,level,typevar) 
    14981284    !       
    14991285    IMPLICIT NONE 
     
    15021288    CHARACTER(*),DIMENSION(4),INTENT(in) :: dimname 
    15031289    INTEGER :: time,level 
    1504     REAL*8, DIMENSION(:,:,:,:), POINTER :: tabvar 
     1290    REAL*8, DIMENSION(:,:,:,:), INTENT(in) :: tabvar 
    15051291    ! 
    15061292    ! local variables 
     
    15481334    status = nf90_close(ncid) 
    15491335    !  
    1550   END SUBROUTINE Write_Ncdf_var4d_Real_nt 
    1551  
    1552   SUBROUTINE Write_Ncdf_var0d_Real(varname,file,tabvar,typevar) 
     1336  END SUBROUTINE write_ncdf_var4d_real_nt 
     1337 
     1338  SUBROUTINE write_ncdf_var0d_real(varname,file,tabvar,typevar) 
    15531339    !       
    15541340    IMPLICIT NONE 
     
    15931379    status = nf90_close(ncid) 
    15941380    !  
    1595   END SUBROUTINE Write_Ncdf_var0d_Real 
    1596   !       
    1597   !************************************************************** 
    1598   !   end subroutine Write_Ncdf_var4d_real_nt 
    1599   !************************************************************** 
    1600   ! 
    1601   !**************************************************************** 
    1602   !   subroutine Read_Ncdf_VarName           * 
     1381  END SUBROUTINE write_ncdf_var0d_real 
     1382 
     1383  SUBROUTINE write_ncdf_var0d_int(varname,file,tabvar,typevar) 
     1384    !       
     1385    IMPLICIT NONE 
     1386    !        
     1387    CHARACTER(*),INTENT(in) :: varname,file,typevar 
     1388    INTEGER :: tabvar 
     1389    ! 
     1390    ! local variables 
     1391    ! 
     1392    INTEGER :: status,ncid 
     1393    INTEGER :: varid 
     1394    ! 
     1395    status = nf90_open(file,NF90_WRITE,ncid) 
     1396    IF (status/=nf90_noerr) THEN 
     1397       WRITE(*,*)"unable to open netcdf file : ",file 
     1398       STOP 
     1399    ENDIF 
     1400    !      
     1401    status = nf90_redef(ncid) 
     1402    status = nf90_def_var(ncid,TRIM(varname),nf90_int,varid) 
     1403    status = nf90_enddef(ncid) 
     1404    status = nf90_put_var(ncid,varid,tabvar) 
     1405    !      
     1406    IF (status/=nf90_noerr) THEN     
     1407       WRITE(*,*)"unable to store variable ",varname, & 
     1408       " in file ",file 
     1409       STOP 
     1410    ENDIF 
     1411    status = nf90_close(ncid) 
     1412    !  
     1413  END SUBROUTINE write_ncdf_var0d_int 
     1414 
     1415  ! 
     1416  !**************************************************************** 
     1417  !   subroutine read_ncdf_VarName           * 
    16031418  !                        * 
    16041419  ! subroutine to retrieve of all variables        * 
     
    16101425  !**************************************************************** 
    16111426  ! 
    1612   !************************************************************** 
    1613   !   subroutine Read_Ncdf_VarName 
    1614   !**************************************************************           
    1615   ! 
    1616   SUBROUTINE Read_Ncdf_VarName(filename,tabvarname) 
     1427  ! 
     1428  SUBROUTINE read_ncdf_VarName(filename,tabvarname) 
    16171429    !       
    16181430    CHARACTER(*),INTENT(in) :: filename 
     
    16371449    END DO 
    16381450 
    1639   END SUBROUTINE Read_Ncdf_Varname 
    1640   ! 
    1641   !************************************************************** 
    1642   !   end subroutine Read_Ncdf_VarName 
    1643   !**************************************************************           
    1644   ! 
    1645   !************************************************************** 
    1646   !   subroutine Copy_Ncdf_att 
    1647   !**************************************************************           
    1648   ! 
    1649   SUBROUTINE Copy_Ncdf_att_var(varname,filein,fileout) 
     1451  END SUBROUTINE read_ncdf_Varname 
     1452  ! 
     1453  ! 
     1454  SUBROUTINE copy_ncdf_att_var(varname,filein,fileout) 
    16501455    !       
    16511456    CHARACTER(*),INTENT(in) :: filein,fileout 
     
    17031508    !      print *,'ici2' 
    17041509    !  
    1705   END SUBROUTINE Copy_Ncdf_att_var 
    1706   ! 
    1707   !************************************************************** 
    1708   !   end subroutine Copy_Ncdf_att 
    1709   !**************************************************************          
    1710   !************************************************************** 
    1711   !   subroutine Copy_Ncdf_att 
    1712   !**************************************************************           
    1713   ! 
    1714   SUBROUTINE Copy_Ncdf_att_latlon(varname,filein,fileout,min,max) 
     1510  END SUBROUTINE copy_ncdf_att_var 
     1511  ! 
     1512  ! 
     1513  SUBROUTINE copy_ncdf_att_latlon(varname,filein,fileout,min,max) 
    17151514    !       
    17161515    CHARACTER(*),INTENT(in) :: filein,fileout 
     
    17601559    status = nf90_close(ncid_in) 
    17611560    status = nf90_close(ncid_out) 
    1762   END SUBROUTINE Copy_Ncdf_att_latlon 
     1561  END SUBROUTINE copy_ncdf_att_latlon 
     1562 
    17631563  !************************************************************* 
    17641564  !************************************************************** 
    1765   !************************************************************** 
    1766   !   function Get_NbDims 
    1767   !**************************************************************           
    17681565  ! 
    17691566  INTEGER FUNCTION Get_NbDims( varname , filename ) 
     
    17841581  END FUNCTION Get_NbDims 
    17851582  ! 
    1786   !************************************************************* 
    1787   !************************************************************** 
    1788   !   function Get_NbDims 
    1789   !**************************************************************           
    17901583  ! 
    17911584  LOGICAL FUNCTION Dims_Existence( dimname , filename ) 
     
    18111604  END FUNCTION Dims_Existence 
    18121605  ! 
    1813   !************************************************************* 
    1814   !************************************************************** 
    18151606END MODULE io_netcdf 
Note: See TracChangeset for help on using the changeset viewer.