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.
agrif_create_bathy.f90 in utils/tools/NESTING/src – NEMO

source: utils/tools/NESTING/src/agrif_create_bathy.f90 @ 10027

Last change on this file since 10027 was 10027, checked in by clem, 6 years ago

solve issues with median averages (greatly increase speed and sort out land issues)

  • Property svn:keywords set to Id
File size: 27.0 KB
RevLine 
[2136]1!************************************************************************
2! Fortran 95 OPA Nesting tools                  *
3!                          *
4!     Copyright (C) 2005 Florian Lemarié (Florian.Lemarie@imag.fr)   *
5!                        Laurent Debreu (Laurent.Debreu@imag.fr)  *
6!************************************************************************
7!
8PROGRAM create_bathy
9  !
10  USE NETCDF
11  USE bilinear_interp
12  USE agrif_readwrite
13  USE agrif_partial_steps
14  USE agrif_connect_topo
[10025]15  USE agrif_interpolation
[2136]16  !
17  IMPLICIT NONE
18  !
19  !************************************************************************
20  !                           *
[2455]21  ! PROGRAM  CREATE_BATHY                 *
[2136]22  !                           *
23  ! program to implement bathymetry interpolation to generate     *
24  ! child grid bathymetry file                  *
25  !                           *
26  ! various options :                     *
27  !                           *
28  ! 1- Interpolation directly from parent bathymetry file (z-coord)  *
[2455]29  ! 2- Use new topo file in meters (for example etopo)      *
[2136]30  !                           *
31  ! vertical coordinates permitted : z-coord and partial steps    *
32  ! sigma coordinates is not yet implemented          *
33  !                           *
34  !Interpolation is carried out using bilinear interpolation      *
35  !routine from SCRIP package or median average          *     
36  !                           *
37  !http://climate.lanl.gov/Software/SCRIP/            *
38  !************************************************************************
39  !
40  ! variables declaration
41  !     
[10027]42  CHARACTER(len=80) ::   namelistname
43  CHARACTER*100     ::   Childmeter_file,Childlevel_file,Child_coordinates,child_ps, Child_domcfg   
44  LOGICAL ::   identical_grids     
45  INTEGER ::   nbadd,status,narg,iargc     
46  INTEGER ::   jpj,jpi
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
53  REAL*8, DIMENSION(:)  ,ALLOCATABLE ::   vardep1d
54  REAL*8, DIMENSION(:,:),POINTER     ::   gdepw_ps_interp => NULL() 
55  REAL*8  ::   Cell_lonmin,Cell_lonmax,Cell_latmin,Cell_latmax,wghts
56  LOGICAL ::   Pacifique = .FALSE.
57  INTEGER ::   boundary,iimin,iimax,jjmax,jjmin
58  INTEGER ::   nxhr,nyhr,nxyhr,ji,jj,nbiter
[2136]59
60  TYPE(Coordinates) :: G0,G1 
61  !     
62  narg = iargc()     
63  IF (narg == 0) THEN
64     namelistname = 'namelist.input'
65  ELSE
66     CALL getarg(1,namelistname)
67  ENDIF
68  !
69  ! read input file (namelist.input)
70  CALL read_namelist(namelistname)
71  !     
72  ! define names of child grid files
73  CALL set_child_name(parent_coordinate_file,Child_coordinates) 
[10025]74  IF( TRIM(parent_meshmask_file) .NE. '/NULL' )    CALL set_child_name(parent_meshmask_file,Childlevel_file)           
75  IF( TRIM(parent_meshmask_file) .NE. '/NULL' )    CALL set_child_name(parent_domcfg_output,Child_domcfg) 
[2136]76  !
[10025]77  !                                                   ! ------------------------------------------------------------------
78  IF( .NOT.new_topo .AND. .NOT.partial_steps ) THEN   ! **** First option : no new topo file & no partial steps
79     !                                                !                 ==> interpolate levels directly from parent file
80     !                                                ! ------------------------------------------------------------------
81     WRITE(*,*) '*** First option: no new topo file & no partial steps'
[2136]82     !     
[10025]83     ! read coarse grid bathymetry and coordinates
[2136]84     status = Read_Coordinates(TRIM(parent_coordinate_file),G0)   
[2455]85     status = Read_bathy_level(TRIM(parent_meshmask_file),G0)
[2136]86     !           
[10025]87     ! read fine grid coordinates
88     status = Read_Coordinates(TRIM(Child_coordinates),G1,pacifique)
89     !
90     ! stop if error in size
91     IF( imax > SIZE(G0%glamt,1) .OR. jmax > SIZE(G0%glamt,2) .OR. imax <= imin .OR. jmax <= jmin ) THEN                   
[2136]92        WRITE(*,*) 'ERROR ***** bad child grid definition ...'
93        WRITE(*,*) 'please check imin,jmin,imax,jmax,jpizoom,jpjzoom values'       
94        STOP
95     ENDIF
96     IF( SIZE(G1%nav_lon,1) .NE. nxfin .OR. SIZE(G1%nav_lon,2) .NE. nyfin ) THEN
97        WRITE(*,*) 'ERROR ***** bad child coordinates file ...'
[10025]98        WRITE(*,*) 'please check that child coordinates file has been created with the same namelist'
[2136]99        STOP
100     ENDIF
101     !
[10025]102     jpj = SIZE(G0%nav_lat,2)
103     jpi = SIZE(G0%nav_lat,1)   
[2136]104     !           
105     ! create logical array masksrc
[10025]106     ALLOCATE(masksrc(jpi,jpj))
107     WHERE(G0%bathy_level.LE.0)   ;   masksrc = .FALSE.   ;
108     ELSEWHERE                    ;   masksrc = .TRUE.    ;
[2136]109     END WHERE
110
[10025]111     ! change longitudes (from -180:180 to 0:360)
[2136]112     IF ( Pacifique ) THEN
[10025]113        WHERE(G0%nav_lon < 0.001)   G0%nav_lon = G0%nav_lon + 360.
[2136]114     ENDIF
115
[10025]116     ! From levels to meters
[2136]117     ALLOCATE(G1%bathy_meter(nxfin,nyfin))
118     CALL levels_to_meter(G0)
119     !             
[10025]120     ! compute remapping matrix thanks to SCRIP package
121     CALL get_remap_matrix(G0%nav_lat,G1%nav_lat,G0%nav_lon,G1%nav_lon,masksrc,matrix,src_add,dst_add)
122     CALL make_remap(G0%bathy_meter,G1%bathy_meter,nxfin,nyfin,matrix,src_add,dst_add) 
[2136]123     DEALLOCATE(masksrc)
124     !     
125     ! compute constant bathymetry for Parent-Child bathymetry connection
126     CALL init_constant_bathy(G0%bathy_meter,bathy_fin_constant)
127     !
[10025]128     ! replace child bathymetry by parent bathymetry at the boundaries
129     IF( ln_agrif_domain ) THEN
130        boundary = npt_copy*irafx + nbghostcellsfine + 1 
131     ELSE
132        boundary = npt_copy*irafx
133     ENDIF
[2136]134     G1%bathy_meter(1:boundary,:) = bathy_fin_constant(1:boundary,:)
135     G1%bathy_meter(:,1:boundary) = bathy_fin_constant(:,1:boundary)
136     G1%bathy_meter(nxfin-boundary+1:nxfin,:) = bathy_fin_constant(nxfin-boundary+1:nxfin,:)
137     G1%bathy_meter(:,nyfin-boundary+1:nyfin) = bathy_fin_constant(:,nyfin-boundary+1:nyfin)
[10025]138     DEALLOCATE(bathy_fin_constant)
[2136]139     !                 
[10025]140     ! bathymetry smoothing (everywhere except at the boundaries)
141     IF( smoothing ) THEN
142        IF( ln_agrif_domain ) THEN
143           CALL smooth_topo(G1%bathy_meter(boundary:nxfin-boundary+1,boundary:nyfin-boundary+1),nbiter)
144        ELSE
145           CALL smooth_topo(G1%bathy_meter(boundary+1:nxfin-boundary,boundary+1:nyfin-boundary),nbiter)
146        ENDIF
147     ENDIF
148     !
149     ! From meters to levels
[2136]150     CALL meter_to_levels(G1)
[10025]151     G1%bathy_level=NINT(G1%bathy_level)
152     !       
153     ! Check errors in bathy
[10027]154     DO jj=1,nyfin
155        DO ji=1,nxfin
156           IF (G1%bathy_level(ji,jj).LT.0.) THEN
157              PRINT *,'error in ',ji,jj,G1%bathy_level(ji,jj)
[10025]158              STOP
159           ENDIF
160        ENDDO
161     ENDDO
162     WHERE ((G1%bathy_level.LT.3.).AND.(G1%bathy_level.GT.0.))   G1%bathy_level=3
[2136]163     !
[10025]164     ! remove closed seas
165     IF (removeclosedseas) THEN
166        ALLOCATE(bathy_test(nxfin,nyfin))
167
168        bathy_test=0.
169        WHERE (G1%bathy_level(1,:)     .GT.0.)   bathy_test(1,:)=1
170        WHERE (G1%bathy_level(nxfin,:) .GT.0.)   bathy_test(nxfin,:)=1
171        WHERE (G1%bathy_level(:,1)     .GT.0.)   bathy_test(:,1)=1
172        WHERE (G1%bathy_level(:,nyfin) .GT.0.)   bathy_test(:,nyfin)=1
173
174        nbadd = 1
175        DO WHILE (nbadd.NE.0)
176           nbadd = 0
[10027]177           DO jj=2,nyfin-1
178              DO ji=2,nxfin-1
179                 IF (G1%bathy_level(ji,jj).GT.0.) THEN
180                    IF (MAX(bathy_test(ji,jj+1),bathy_test(ji,jj-1),bathy_test(ji-1,jj),bathy_test(ji+1,jj)).EQ.1) THEN
181                       IF (bathy_test(ji,jj).NE.1.) nbadd = nbadd + 1
182                       bathy_test(ji,jj)=1.
[10025]183                    ENDIF
184                 ENDIF
185              ENDDO
186           ENDDO
187        ENDDO
188
189        WHERE (bathy_test.EQ.0.)   G1%bathy_level = 0.
190
191        DEALLOCATE(bathy_test)           
192     ENDIF
[2136]193     !
[10025]194     ! store interpolation result in output file
195     CALL levels_to_meter(G1)   ! From levels to meters
196     status = Write_Bathy_level(TRIM(Childlevel_file),G1)
197     status = write_domcfg(TRIM(Child_domcfg),G1)
198
199     WRITE(*,*) '****** Bathymetry successfully created if no new topo ******'
200     STOP
[2136]201     !
[10025]202     !                                                ! -----------------------------------------------------
203  ELSE                                                ! **** Second option : new topo file or partial steps     
204     !                                                ! -----------------------------------------------------
205     WRITE(*,*) '***Second option : new topo or partial steps'
206
[2136]207     ! read fine and coarse grids coordinates file
208     status = Read_Coordinates(TRIM(parent_coordinate_file),G0)
209     status = Read_Coordinates(TRIM(Child_coordinates),G1,Pacifique)
[10025]210     !
211     ! check error in size
212     IF( imax > SIZE(G0%nav_lon,1) .OR. jmax > SIZE(G0%nav_lon,2) .OR. imax <= imin .OR. jmax <= jmin ) THEN                   
[2136]213        WRITE(*,*) 'ERROR ***** bad child grid definition ...'
214        WRITE(*,*) 'please check imin,jmin,imax,jmax,jpizoom,jpjzoom values'       
215        STOP
216     ENDIF
217     IF( SIZE(G1%nav_lon,1) .NE. nxfin .OR. SIZE(G1%nav_lon,2) .NE. nyfin ) THEN
218        WRITE(*,*) 'ERROR ***** bad child coordinates file ...'
[10025]219        WRITE(*,*) 'please check that child coordinates file has been created with the same namelist'
[2136]220        STOP
221     ENDIF
[10025]222     !
223     ! === From here on: G0 is the grid associated with the new topography (like gebco or etopo) ===
224     !
225     ! read bathymetry data set => G0%bathy_meter
226     IF( new_topo ) THEN                                       ! read G0%bathy_meter (on a reduced grid) and G1 coordinates
[2136]227        DEALLOCATE( G0%nav_lon, G0%nav_lat )
228        status = Read_bathy_meter(TRIM(elevation_database),G0,G1,Pacifique)
[10025]229     ELSE                                                      ! read G0%bathy_meter (on the global grid)
230        status = Read_Bathymeter(TRIM(parent_bathy_meter),G0)   
[2136]231        IF(Pacifique) THEN
[10025]232           WHERE(G0%nav_lon < 0.0001)   G0%nav_lon = G0%nav_lon + 360.
[2136]233        ENDIF
234     ENDIF
235     !               
[10025]236     ! what type of interpolation for bathymetry
[2136]237     IF( type_bathy_interp == 0 ) THEN
[10025]238        WRITE(*,*) 'Interpolation of high resolution bathymetry on child grid: Arithmetic average ...'
[2136]239     ELSE IF( type_bathy_interp == 1 ) THEN
[10025]240        WRITE(*,*) 'Interpolation of high resolution bathymetry on child grid: Median average ...'
[2136]241     ELSE IF( type_bathy_interp == 2 ) THEN     
[10025]242        WRITE(*,*) 'Interpolation of high resolution bathymetry on child grid: Bilinear interpolation ...'
[2136]243     ELSE     
[10025]244        WRITE(*,*) 'bad value for type_bathy_interp variable ( must be 0, 1 or 2 )'
[2136]245        STOP
246     ENDIF
247     !
[10025]248     ! 1st allocation of child grid bathy
249     ALLOCATE(G1%bathy_meter(nxfin,nyfin))
250     G1%bathy_meter(:,:)=0.                       
[2136]251     !
[10025]252     ! ---------------------------------------------------------------------------------
253     ! ===                 Bathymetry of the fine grid (step1)                       ===
254     ! ---------------------------------------------------------------------------------
255     ! ==> It gives G1%bathy_meter from G0%bathy_meter
256     ! ---------------------------------------------------------------------------------
257     !                                                               ! -----------------------------
258     IF( type_bathy_interp == 0 .OR. type_bathy_interp == 1 ) THEN   ! arithmetic or median averages
259        !                                                            ! -----------------------------
[2136]260        ALLOCATE(trouble_points(nxfin,nyfin))
[10025]261        trouble_points(:,:) = 0
[2136]262        !                       
[10027]263        DO jj = 2, nyfin
264           DO ji = 2, nxfin
[2136]265              !       
[10025]266              ! fine grid cell extension               
[2136]267              Cell_lonmin = MIN(G1%glamf(ji-1,jj-1),G1%glamf(ji,jj-1),G1%glamf(ji,jj),G1%glamf(ji-1,jj))
268              Cell_lonmax = MAX(G1%glamf(ji-1,jj-1),G1%glamf(ji,jj-1),G1%glamf(ji,jj),G1%glamf(ji-1,jj))
269              Cell_latmin = MIN(G1%gphif(ji-1,jj-1),G1%gphif(ji,jj-1),G1%gphif(ji,jj),G1%gphif(ji-1,jj))
[9694]270              Cell_latmax = MAX(G1%gphif(ji-1,jj-1),G1%gphif(ji,jj-1),G1%gphif(ji,jj),G1%gphif(ji-1,jj)) 
[2136]271              !               
[10025]272              ! look for points in G0 (bathy dataset) contained in the fine grid cells 
[2136]273              iimin = 1
274              DO WHILE( G0%nav_lon(iimin,1) < Cell_lonmin ) 
275                 iimin = iimin + 1
276              ENDDO
277              !               
278              jjmin = 1
279              DO WHILE( G0%nav_lat(iimin,jjmin) < Cell_latmin ) 
280                 jjmin = jjmin + 1
281              ENDDO
282              !               
283              iimax = iimin 
284              DO WHILE( G0%nav_lon(iimax,1) <= Cell_lonmax ) 
285                 iimax = iimax + 1
[10025]286                 iimax = MIN( iimax,SIZE(G0%bathy_meter,1))
[2136]287              ENDDO
288              !                               
289              jjmax = jjmin 
290              DO WHILE( G0%nav_lat(iimax,jjmax) <= Cell_latmax ) 
291                 jjmax = jjmax + 1
[10025]292                 jjmax = MIN( jjmax,SIZE(G0%bathy_meter,2))
[2136]293              ENDDO
294              !
[10025]295              IF( ln_agrif_domain ) THEN
296                 iimax = iimax-1
297                 jjmax = jjmax-1
298              ELSE
299                 iimax = MAX(iimin,iimax-1)
300                 jjmax = MAX(jjmin,jjmax-1)
301              ENDIF
[2136]302              !               
303              iimin = MAX( iimin,1 )
304              jjmin = MAX( jjmin,1 )
305              iimax = MIN( iimax,SIZE(G0%bathy_meter,1))
306              jjmax = MIN( jjmax,SIZE(G0%bathy_meter,2))
307
308              nxhr = iimax - iimin + 1
309              nyhr = jjmax - jjmin + 1                   
310
311              IF( nxhr == 0 .OR. nyhr == 0 ) THEN
[10025]312                 !
[2136]313                 trouble_points(ji,jj) = 1
[10025]314                 !
[2136]315              ELSE
[10025]316                 !
317                 ALLOCATE( vardep(nxhr,nyhr), mask_oce(nxhr,nyhr) )
[2136]318                 vardep(:,:) = G0%bathy_meter(iimin:iimax,jjmin:jjmax)
[10025]319                 !
320                 WHERE( vardep(:,:) .GT. 0. )   ;   mask_oce = 1 ;
321                 ELSEWHERE                      ;   mask_oce = 0 ;
322                 ENDWHERE
323                 !
[10027]324                 nxyhr = nxhr*nyhr
325                 IF( SUM(mask_oce) < 0.5*(nxyhr) ) THEN ! if more than half of the points are on land then bathy fine = 0
[2136]326                    G1%bathy_meter(ji,jj) = 0.
327                 ELSE
[10025]328                    IF( type_bathy_interp == 0 ) THEN     ! Arithmetic average
329                       G1%bathy_meter(ji,jj) = SUM( vardep(:,:) * mask_oce(:,:) ) / SUM( mask_oce(:,:) )
330                    ELSE                                  ! Median average
[10027]331                       ALLOCATE(vardep1d(nxyhr))
332                       vardep1d = RESHAPE(vardep,(/ nxyhr /) )
333                       !!CALL ssort(vardep1d,nxyhr)
334                       CALL quicksort(vardep1d,1,nxyhr)
[2136]335                       !
336                       ! Calculate median
[10027]337                       IF (MOD(nxyhr,2) .NE. 0) THEN
338                          G1%bathy_meter(ji,jj) = vardep1d( nxyhr/2 + 1 )
[2136]339                       ELSE
[10027]340                          G1%bathy_meter(ji,jj) = 0.5 * ( vardep1d(nxyhr/2) + vardep1d(nxyhr/2+1) )
[2136]341                       END IF
[9694]342                       DEALLOCATE(vardep1d)   
[2136]343                    ENDIF
344                 ENDIF
345                 DEALLOCATE (mask_oce,vardep)
[10025]346                 !
[2136]347              ENDIF
348           ENDDO
349        ENDDO
350
351        IF( SUM( trouble_points ) > 0 ) THEN
[10025]352           PRINT*,'too much empty cells, proceed to bilinear interpolation'
[2136]353           type_bathy_interp = 2
354        ENDIF
[10025]355           
[2136]356     ENDIF
[10025]357     !                                                       ! -----------------------------
358     IF( type_bathy_interp == 2) THEN                        ! Bilinear interpolation
359        !                                                    ! -----------------------------
[2136]360        identical_grids = .FALSE.
361
[10025]362        IF( SIZE(G0%nav_lat,1) == SIZE(G1%nav_lat,1) .AND. SIZE(G0%nav_lat,2) == SIZE(G1%nav_lat,2)  .AND.   &
363            SIZE(G0%nav_lon,1) == SIZE(G1%nav_lon,1) .AND. SIZE(G0%nav_lon,2) == SIZE(G1%nav_lon,2) ) THEN
[2136]364           IF( MAXVAL( ABS(G0%nav_lat(:,:)- G1%nav_lat(:,:)) ) < 0.0001 .AND.   &
[10025]365               MAXVAL( ABS(G0%nav_lon(:,:)- G1%nav_lon(:,:)) ) < 0.0001 ) THEN
366              PRINT*,'same grid between ',elevation_database,' and child domain'   
[2136]367              G1%bathy_meter = G0%bathy_meter 
368              identical_grids = .TRUE.                         
369           ENDIF
370        ENDIF
371
372        IF( .NOT. identical_grids ) THEN
373
374           ALLOCATE(masksrc(SIZE(G0%bathy_meter,1),SIZE(G0%bathy_meter,2)))
375           ALLOCATE(bathy_test(nxfin,nyfin))
376           !
[10025]377           !Where(G0%bathy_meter.le.0.00001)
378           !  masksrc = .false.
379           !ElseWhere
380              masksrc = .TRUE.
381           !End where                       
[2136]382           !           
383           ! compute remapping matrix thanks to SCRIP package           
[10025]384           CALL get_remap_matrix(G0%nav_lat,G1%nav_lat,G0%nav_lon,G1%nav_lon,masksrc,matrix,src_add,dst_add)
385           CALL make_remap(G0%bathy_meter,bathy_test,nxfin,nyfin,matrix,src_add,dst_add) 
[2136]386           !                                 
387           G1%bathy_meter = bathy_test               
388           !           
389           DEALLOCATE(masksrc)
390           DEALLOCATE(bathy_test) 
391
392        ENDIF
393        !           
394     ENDIF
[10025]395     ! ---
396     ! At this stage bathymetry in meters has already been interpolated on fine grid
397     !                    => G1%bathy_meter(nxfin,nyfin)
[2136]398     !
[10025]399     ! Also G0 was the grid from the new bathymetry data set (etopo, gebco...) and not the coarse grid
400     ! ---
[2136]401     !
[10025]402     ! ---------------------------------------------------------------------------------
403     ! ===                 Bathymetry of the fine grid (step2)                       ===
404     ! ---------------------------------------------------------------------------------
405     ! ==> It gives an update of G1%bathy_meter and G1%bathy_level
406     ! ---------------------------------------------------------------------------------
407     ! From here on: G0 is the coarse grid
[2136]408     !
[10025]409     ! Coarse grid bathymetry : G0%bathy_meter (on the global grid)
410     IF( partial_steps ) THEN                     
411        status = Read_Bathymeter(TRIM(parent_bathy_meter),G0)
412     ELSE
413        status = Read_Bathymeter(TRIM(parent_meshmask_file),G0)
414     ENDIF
415     
416     ! Coarse grid coordinatees : G0 coordinates
417     DEALLOCATE(G0%nav_lat,G0%nav_lon)
418     status = Read_coordinates(TRIM(parent_coordinate_file),G0)
419
420     ! allocate temporary arrays                 
421     IF (.NOT.ASSOCIATED(G0%gdepw_ps))       ALLOCATE(G0%gdepw_ps    (SIZE(G0%bathy_meter,1),SIZE(G0%bathy_meter,2)))
422     IF (.NOT.ASSOCIATED(G1%gdepw_ps))       ALLOCATE(G1%gdepw_ps    (SIZE(G1%bathy_meter,1),SIZE(G1%bathy_meter,2)))                 
423     IF (.NOT.ASSOCIATED(gdepw_ps_interp))   ALLOCATE(gdepw_ps_interp(SIZE(G1%bathy_meter,1),SIZE(G1%bathy_meter,2)))
424     !                       
425     IF( ln_agrif_domain ) THEN
426        boundary = npt_copy*irafx + nbghostcellsfine + 1
427     ELSE
428        boundary = npt_copy*irafx
429     ENDIF
[2136]430     !
[10025]431     ! compute G0%gdepw_ps and G1%gdepw_ps
432     CALL get_partial_steps(G0) 
433     CALL get_partial_steps(G1)
434     CALL bathymetry_control(G0%Bathy_level)
435
436     ! ---------------------------------------
437     ! Bathymetry at the boundaries (npt_copy)                     
438     ! ---------------------------------------
439     ! 1st step: interpolate coarse bathy on the fine grid (using partial steps or not)
440     IF( partial_steps .AND. ln_agrif_domain ) THEN                   
441        CALL Check_interp(G0,gdepw_ps_interp)
442     ELSE
443        gdepw_ps_interp = 0. * G1%gdepw_ps
444        !!CALL agrif_interp(G0%gdepw_ps,gdepw_ps_interp,'T')
445        CALL init_constant_bathy(G0%gdepw_ps,gdepw_ps_interp)
446     ENDIF
447
448     IF (.NOT.ASSOCIATED(G1%wgt))   ALLOCATE(G1%wgt(SIZE(G1%bathy_meter,1),SIZE(G1%bathy_meter,2)))
449     G1%wgt(:,:) = 0.
450     IF ((.NOT.ASSOCIATED(G0%wgt)).AND.bathy_update) THEN
451        ALLOCATE(G0%wgt(SIZE(G0%nav_lat,1),SIZE(G0%nav_lat,2)))
452        G0%wgt(:,:) = 0.
453     ENDIF
[2136]454     !
[10025]455     ! 2nd step: copy parent bathymetry at the boundaries
456     DO jj=1,nyfin   ! West and East
457        IF ( gdepw_ps_interp(nbghostcellsfine+1,jj) > 0. ) THEN
458           G1%gdepw_ps(1:boundary,jj) = gdepw_ps_interp(1:boundary,jj) 
459           G1%wgt(1:boundary,jj) = 1.
460        ELSE
461           G1%gdepw_ps(1:nbghostcellsfine+1,jj)=0. 
462        ENDIF
[2136]463        !
[10025]464        IF ( gdepw_ps_interp(nxfin-nbghostcellsfine,jj) > 0.) THEN
465           G1%gdepw_ps(nxfin-boundary+1:nxfin,jj)=gdepw_ps_interp(nxfin-boundary+1:nxfin,jj)
466           G1%wgt(nxfin-boundary+1:nxfin,jj) = 1.
467        ELSE
468           G1%gdepw_ps(nxfin-nbghostcellsfine:nxfin,jj) = 0.
469        ENDIF
470     END DO
471     !
472     DO ji=1,nxfin    ! South and North
473        IF (gdepw_ps_interp(ji,nbghostcellsfine+1)>0.) THEN
474           G1%gdepw_ps(ji,1:boundary) = gdepw_ps_interp(ji,1:boundary)
475           G1%wgt(ji,1:boundary) = 1.
476        ELSE
477           G1%gdepw_ps(ji,1:nbghostcellsfine+1)=0. 
478        ENDIF
[2136]479        !
[10025]480        IF (gdepw_ps_interp(ji,nyfin-nbghostcellsfine)>0.) THEN
481           G1%gdepw_ps(ji,nyfin-boundary+1:nyfin)=gdepw_ps_interp(ji,nyfin-boundary+1:nyfin)
482           G1%wgt(ji,nyfin-boundary+1:nyfin) = 1.
483        ELSE
484           G1%gdepw_ps(ji,nyfin-nbghostcellsfine:nyfin) = 0.
[9749]485        ENDIF
[10025]486     END DO
487     !
488     !clem: recalculate interpolation everywhere before linear connection (useless to me)
489     IF( partial_steps .AND. ln_agrif_domain ) THEN               
490        gdepw_ps_interp = 0.
491        CALL Check_interp(G0,gdepw_ps_interp)
492     ENDIF
493     !
494     ! -------------------------------------------------------
495     ! Bathymetry between boundaries and interior (npt_connect)                 
496     ! --------------------------------------------------------
497     ! Make linear connection (on npt_connect*irafx points) between the boundaries and the interior
498     IF( ln_agrif_domain ) THEN
499        boundary = (npt_copy + npt_connect)*irafx + nbghostcellsfine + 1
500     ELSE
501        boundary = (npt_copy + npt_connect)*irafx
502     ENDIF
[2136]503
[10025]504     IF( npt_connect > 0 ) THEN
[9694]505        WRITE(*,*) ' linear connection on ',npt_connect,'coarse grid points'
[2136]506
507        wghts = 1.
[10025]508        DO ji = boundary - npt_connect*irafx + 1 , boundary
509           wghts = wghts - (1. / (npt_connect*irafx + 1. ) )
510           DO jj=1,nyfin
511              IF (G1%gdepw_ps(nbghostcellsfine+1,jj) > 0.)       G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj)) 
512           END DO
513        END DO
514
[9149]515        wghts = 1.
[10025]516        DO ji = nxfin - (boundary - npt_connect*irafx), nxfin - boundary +1 , -1
517           wghts = wghts - (1. / (npt_connect*irafx + 1. ) )
518           DO jj=1,nyfin
519              IF (G1%gdepw_ps(nxfin-nbghostcellsfine,jj) > 0.)   G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj))
520           END DO
521        END DO
[2136]522
523        wghts = 1.
[10025]524        DO jj = boundary - npt_connect*irafy + 1 , boundary
525           wghts = wghts - (1. / (npt_connect*irafy + 1. ) )
526           DO ji=1,nxfin
527              IF (G1%gdepw_ps(ji,nbghostcellsfine+1) > 0.)       G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj))
528           END DO
[9149]529        END DO
530
[2136]531        wghts = 1.
[10025]532        DO jj = nyfin - (boundary - npt_connect*irafy) , nyfin - boundary +1, -1
533           wghts = wghts - (1. / (npt_connect*irafy + 1. ) )
534           DO ji=1,nxfin
535              IF (G1%gdepw_ps(ji,nyfin-nbghostcellsfine) > 0.)   G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj))
536           END DO
[9149]537        END DO
[9753]538        IF (.NOT.identical_grids) THEN
539           G1%gdepw_ps(:,:) = (1.-G1%wgt(:,:)) * G1%gdepw_ps(:,:) + gdepw_ps_interp(:,:)*G1%wgt(:,:)
540        ENDIF
[2136]541
[10025]542     ENDIF
[2136]543
[10025]544     ! replace G1%bathy_meter by G1%gdepw_ps
545     G1%bathy_meter = G1%gdepw_ps
546     !                     
547     ! --------------------
548     ! Bathymetry smoothing                 
549     ! --------------------
550     IF( smoothing .AND. (.NOT.identical_grids) ) THEN
551        ! Chanut: smoothing everywhere then discard result in connection zone
552        CALL smooth_topo(G1%gdepw_ps(:,:),nbiter)
553        WHERE (G1%wgt(:,:)==0.) G1%bathy_meter(:,:) = G1%gdepw_ps(:,:)
554     ELSE
555        WRITE(*,*) 'No smoothing process only connection is carried out'
556     ENDIF
557     !
558     ! ------------------
559     ! Remove closed seas
560     ! ------------------
561     IF (removeclosedseas) THEN
562        ALLOCATE(bathy_test(nxfin,nyfin))
563        bathy_test=0.
564        WHERE (G1%bathy_meter(1,:)    .GT.0.)   bathy_test(1,:)=1
565        WHERE (G1%bathy_meter(nxfin,:).GT.0.)   bathy_test(nxfin,:)=1
566        WHERE (G1%bathy_meter(:,1)    .GT.0.)   bathy_test(:,1)=1
567        WHERE (G1%bathy_meter(:,nyfin).GT.0.)   bathy_test(:,nyfin)=1
[2136]568
[10025]569        nbadd = 1
570        DO WHILE (nbadd.NE.0)
571           nbadd = 0
[10027]572           DO jj=2,nyfin-1
573              DO ji=2,nxfin-1
574                 IF (G1%bathy_meter(ji,jj).GT.0.) THEN
575                    IF (MAX(bathy_test(ji,jj+1),bathy_test(ji,jj-1),bathy_test(ji-1,jj),bathy_test(ji+1,jj)).EQ.1) THEN
576                       IF (bathy_test(ji,jj).NE.1.) nbadd = nbadd + 1
577                       bathy_test(ji,jj)=1.
[10025]578                    ENDIF
[2136]579
[10025]580                 ENDIF
[2136]581              ENDDO
582           ENDDO
[10025]583        ENDDO
584        WHERE (bathy_test.EQ.0.)   G1%bathy_meter = 0.
585        DEALLOCATE(bathy_test)
586     ENDIF
587     !
588     !                              ! ----------------
589     IF( partial_steps ) THEN       ! If partial steps
590        !                           ! ----------------
[2136]591        !
[10025]592        CALL get_partial_steps(G1)  ! recompute gdepw_ps for G1
[9628]593
[10025]594        IF(bathy_update)   CALL Update_Parent_Bathy( G0,G1 ) 
[2136]595        CALL set_child_name(parent_bathy_meter,child_ps)
[10025]596       
597                           status = Write_Bathy_meter(TRIM(child_ps),G1)
598        IF(bathy_update)   status = Write_Bathy_meter(TRIM(updated_parent_file),G0)
599        IF(bathy_update)   status = write_domcfg(TRIM(updated_parent_domcfg),G0)
[2136]600
601        CALL get_partial_steps(G1)
602        !
603        G1%bathy_level=NINT(G1%bathy_level)
604        !
[10025]605        ! store interpolation result in output file
606        CALL levels_to_meter(G1)   ! From levels to meters
607        IF( TRIM(parent_meshmask_file) .NE. '/NULL' )    status = Write_Bathy_level(TRIM(Childlevel_file),G1)
608        IF( TRIM(parent_meshmask_file) .NE. '/NULL' )    status = write_domcfg(TRIM(Child_domcfg),G1)
[2136]609        !
610        WRITE(*,*) '****** Bathymetry successfully created for partial cells ******'
[10025]611        STOP
[2136]612        !
[10025]613        !                           ! ----------------
614     ELSE                           ! No partial steps
615        !                           ! ----------------
616        !
617        ! convert bathymetry from meters to levels
618        CALL meter_to_levels(G1) 
619
620        IF(bathy_update)   CALL Update_Parent_Bathy( G0,G1 ) 
621        CALL set_child_name(parent_bathy_meter,child_ps)
622       
623                           status = Write_Bathy_meter(TRIM(child_ps),G1)       
624        IF(bathy_update)   status = Write_Bathy_meter(TRIM(updated_parent_file),G0)
625        IF(bathy_update)   status = write_domcfg(TRIM(updated_parent_domcfg),G0)
626        !
627        ! store interpolation result in output file
628        CALL levels_to_meter(G1)   ! From levels to meters
629        status = Write_Bathy_level(TRIM(Childlevel_file),G1)
630        status = write_domcfg(TRIM(Child_domcfg),G1)
631        !
632        WRITE(*,*) '****** Bathymetry successfully created for full cells ******'
633        STOP
634        !
[2136]635     ENDIF
636  ENDIF
637  !
638END PROGRAM create_bathy
639
640
Note: See TracBrowser for help on using the repository browser.