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 branches/2017/dev_merge_2017/NEMOGCM/TOOLS/NESTING/src – NEMO

source: branches/2017/dev_merge_2017/NEMOGCM/TOOLS/NESTING/src/agrif_create_bathy.f90 @ 9149

Last change on this file since 9149 was 9149, checked in by jchanut, 7 years ago

NESTING TOOLS
Add new e3t definition, #1981
Correct transition weights, #1998
Correct arithmetic and median average interpolations, #1999
Changed number of coarse grid bathymetry points inside child grid domain (2 now, instead of 3)
  • Property svn:keywords set to Id
File size: 24.1 KB
Line 
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
15  !
16  IMPLICIT NONE
17  !
18  !************************************************************************
19  !                           *
20  ! PROGRAM  CREATE_BATHY                 *
21  !                           *
22  ! program to implement bathymetry interpolation to generate     *
23  ! child grid bathymetry file                  *
24  !                           *
25  ! various options :                     *
26  !                           *
27  ! 1- Interpolation directly from parent bathymetry file (z-coord)  *
28  ! 2- Use new topo file in meters (for example etopo)      *
29  !                           *
30  ! vertical coordinates permitted : z-coord and partial steps    *
31  ! sigma coordinates is not yet implemented          *
32  !                           *
33  !Interpolation is carried out using bilinear interpolation      *
34  !routine from SCRIP package or median average          *     
35  !                           *
36  !http://climate.lanl.gov/Software/SCRIP/            *
37  !************************************************************************
38  !
39  ! variables declaration
40  !     
41  CHARACTER(len=80) :: namelistname
42  CHARACTER*100 :: Childmeter_file,Childlevel_file,Child_coordinates,child_ps     
43  LOGICAL,DIMENSION(:,:),POINTER :: masksrc => NULL() 
44  LOGICAL :: identical_grids     
45  INTEGER,DIMENSION(:,:),ALLOCATABLE ::mask_oce,trouble_points
46  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
55  REAL*8, DIMENSION(:,:),POINTER :: wtab  => NULL()
56  REAL*8  :: Cell_lonmin,Cell_lonmax,Cell_latmin,Cell_latmax,wghts
57  LOGICAL :: Pacifique=.FALSE.
58  INTEGER :: boundary,xpos,ypos,iimin,iimax,jjmax,jjmin
59  INTEGER :: nbloops,nxhr,nyhr,ji,jj,nbiter,nbloopmax
60  INTEGER :: ipt,jpt,iloc,jloc
61  INTEGER, DIMENSION(2) :: i_min,i_max,j_min,j_max
62
63  TYPE(Coordinates) :: G0,G1 
64  !     
65  narg = iargc()     
66  IF (narg == 0) THEN
67     namelistname = 'namelist.input'
68  ELSE
69     CALL getarg(1,namelistname)
70  ENDIF
71  !
72  ! read input file (namelist.input)
73  !
74  CALL read_namelist(namelistname)
75  !     
76  ! define names of child grid files
77  !
78  CALL set_child_name(parent_coordinate_file,Child_coordinates) 
79  IF( TRIM(parent_meshmask_file) .NE. '/NULL' ) &
80       CALL set_child_name(parent_meshmask_file,Childlevel_file)           
81  !
82  !
83  !
84  !
85  !
86  !------------------------------------------------------------------
87  ! ****First option : no new topo file / no partial steps
88  ! interpolate levels directly from parent file
89  !------------------------------------------------------------------
90  !
91  !
92  !
93  !
94  !
95  !
96  IF(.NOT.new_topo .AND. .NOT.partial_steps) THEN     
97     !     
98     WRITE(*,*) 'First option'
99     !read coarse grid bathymetry and coordinates file
100     !           
101     WRITE(*,*) 'No new topo file ...'
102     status = Read_Coordinates(TRIM(parent_coordinate_file),G0)   
103     status = Read_bathy_level(TRIM(parent_meshmask_file),G0)
104     !           
105     IF( imax > SIZE(G0%glamt,1) .OR. jmax > SIZE(G0%glamt,2) .OR. &
106          imax <= imin .OR. jmax <= jmin ) THEN                   
107        WRITE(*,*) 'ERROR ***** bad child grid definition ...'
108        WRITE(*,*) 'please check imin,jmin,imax,jmax,jpizoom,jpjzoom values'       
109        WRITE(*,*) ' '
110        STOP
111     ENDIF
112     !
113     !read fine grid coordinates file
114     !     
115     status = Read_Coordinates(TRIM(Child_coordinates),G1,pacifique)
116     ! 
117     IF( SIZE(G1%nav_lon,1) .NE. nxfin .OR. SIZE(G1%nav_lon,2) .NE. nyfin ) THEN
118        !
119        WRITE(*,*) 'ERROR ***** bad child coordinates file ...'
120        WRITE(*,*) ' '
121        WRITE(*,*) 'please check that child coordinates file '
122        WRITE(*,*) 'has been created with the same namelist '
123        WRITE(*,*) ' '
124        STOP
125        !
126     ENDIF
127     !
128     numlat =  SIZE(G0%nav_lat,2)
129     numlon =  SIZE(G0%nav_lat,1)   
130     numlatfine =  SIZE(G1%nav_lat,2)
131     numlonfine =  SIZE(G1%nav_lat,1) 
132     !           
133     ALLOCATE(masksrc(numlon,numlat))
134     !
135     ! create logical array masksrc
136     !
137     WHERE(G0%bathy_level.LE.0) 
138        masksrc = .FALSE.
139     ELSEWHERE
140        masksrc = .TRUE.
141     END WHERE
142
143     IF ( Pacifique ) THEN
144        WHERE(G0%nav_lon < 0.001) 
145           G0%nav_lon = G0%nav_lon + 360.
146        END WHERE
147     ENDIF
148
149
150     !-----------------         
151     ! compute remapping matrix thanks to SCRIP package
152     !
153     ! remapping process
154     !             
155     ALLOCATE(G1%bathy_meter(nxfin,nyfin))
156     CALL levels_to_meter(G0)
157     !             
158     !             Call levels_to_meter(G1)
159     !             
160     CALL get_remap_matrix(G0%nav_lat,G1%nav_lat,   &
161          G0%nav_lon,G1%nav_lon,   &
162          masksrc,matrix,src_add,dst_add)
163     CALL make_remap(G0%bathy_meter,G1%bathy_meter,nxfin,nyfin, &
164          matrix,src_add,dst_add) 
165     !                                 
166     !           
167     DEALLOCATE(masksrc)
168     !-----------------                                     
169     !     
170     !
171     ! compute constant bathymetry for Parent-Child bathymetry connection
172     !             
173     CALL init_constant_bathy(G0%bathy_meter,bathy_fin_constant)
174     !
175     boundary = connectionsize*irafx + nbghostcellsfine 
176     !
177     ! connection carried out by copying parent grid values for the fine points
178     ! corresponding to 3 coarse grid cells at the boundaries
179     !                 
180     G1%bathy_meter(1:boundary,:) = bathy_fin_constant(1:boundary,:)
181     G1%bathy_meter(:,1:boundary) = bathy_fin_constant(:,1:boundary)
182     G1%bathy_meter(nxfin-boundary+1:nxfin,:) = bathy_fin_constant(nxfin-boundary+1:nxfin,:)
183     G1%bathy_meter(:,nyfin-boundary+1:nyfin) = bathy_fin_constant(:,nyfin-boundary+1:nyfin)
184     !                 
185     CALL smooth_topo(G1%bathy_meter(boundary:nxfin-boundary+1,boundary:nyfin-boundary+1),nbiter)
186     !             
187     CALL meter_to_levels(G1)
188     !             
189     DEALLOCATE(bathy_fin_constant)
190     !
191     !
192     !------------------------------------------------------------------
193     ! ****Second option : new topo file or/and partial steps     
194     !------------------------------------------------------------------
195     !
196     !
197     !
198     !
199     !
200     !
201     !
202     !
203  ELSE
204     !
205     WRITE(*,*) 'Second option : partial steps'
206     ! read fine and coarse grids coordinates file
207     !       
208     status = Read_Coordinates(TRIM(parent_coordinate_file),G0)
209     status = Read_Coordinates(TRIM(Child_coordinates),G1,Pacifique)
210     !                       
211     IF( imax > SIZE(G0%nav_lon,1) .OR. jmax > SIZE(G0%nav_lon,2) .OR. &
212          imax <= imin .OR. jmax <= jmin ) THEN                   
213        WRITE(*,*) 'ERROR ***** bad child grid definition ...'
214        WRITE(*,*) 'please check imin,jmin,imax,jmax,jpizoom,jpjzoom values'       
215        WRITE(*,*) ' '
216        STOP
217     ENDIF
218     !
219
220     ! 
221     IF( SIZE(G1%nav_lon,1) .NE. nxfin .OR. SIZE(G1%nav_lon,2) .NE. nyfin ) THEN
222        !
223        WRITE(*,*) 'ERROR ***** bad child coordinates file ...'
224        WRITE(*,*) ' '
225        WRITE(*,*) 'please check that child coordinates file '
226        WRITE(*,*) 'has been created with the same namelist '
227        WRITE(*,*) ' '
228        STOP
229        !
230     ENDIF
231     !     
232     ! read coarse grid bathymetry file
233     !---
234     IF( new_topo ) THEN
235        WRITE(*,*) 'use new topo file : ',TRIM(elevation_database)
236        DEALLOCATE( G0%nav_lon, G0%nav_lat )
237        status = Read_bathy_meter(TRIM(elevation_database),G0,G1,Pacifique)
238     ELSE
239        WRITE(*,*) 'no new topo file'
240        status = Read_Bathymeter(TRIM(parent_bathy_meter),G0)
241        IF(Pacifique) THEN
242           WHERE(G0%nav_lon < 0.0001) 
243              G0%nav_lon = G0%nav_lon + 360.
244           END WHERE
245        ENDIF
246     ENDIF
247     !---           
248     numlatfine =  SIZE(G1%nav_lat,2)
249     numlonfine =  SIZE(G1%nav_lat,1) 
250     ! 
251     !               
252     ALLOCATE(G1%bathy_meter(nxfin,nyfin))
253     G1%bathy_meter(:,:)=0.                       
254
255     WRITE(*,*) 'Interpolation of high resolution bathymetry on child grid'
256
257     IF( type_bathy_interp == 0 ) THEN
258        WRITE(*,*) 'Arithmetic average ...'
259     ELSE IF( type_bathy_interp == 1 ) THEN
260        WRITE(*,*) 'Median average ...'
261     ELSE IF( type_bathy_interp == 2 ) THEN     
262        WRITE(*,*) 'Bilinear interpolation ...'
263     ELSE     
264        WRITE(*,*) 'bad value for type_bathy_interp variable ( must be 0,1 or 2 )'
265        STOP
266     ENDIF
267     !
268     !************************************
269     !MEDIAN AVERAGE or ARITHMETIC AVERAGE
270     !************************************
271     !
272     IF( type_bathy_interp == 0 .OR. type_bathy_interp == 1 ) THEN 
273        !
274        ALLOCATE(trouble_points(nxfin,nyfin))
275        trouble_points = 0
276        !
277        !  POINT DETECTION
278        !
279        !                       
280        DO jj = 2,numlatfine
281           DO ji = 2,numlonfine
282              !       
283              ! FINE GRID CELLS DEFINITION               
284              !
285              Cell_lonmin = MIN(G1%glamf(ji-1,jj-1),G1%glamf(ji,jj-1),G1%glamf(ji,jj),G1%glamf(ji-1,jj))
286              Cell_lonmax = MAX(G1%glamf(ji-1,jj-1),G1%glamf(ji,jj-1),G1%glamf(ji,jj),G1%glamf(ji-1,jj))
287              Cell_latmin = MIN(G1%gphif(ji-1,jj-1),G1%gphif(ji,jj-1),G1%gphif(ji,jj),G1%gphif(ji-1,jj))
288              Cell_latmax = MAX(G1%gphif(ji-1,jj-1),G1%gphif(ji,jj-1),G1%gphif(ji,jj),G1%gphif(ji-1,jj))                   
289              !               
290              ! SEARCH FOR ALL POINTS CONTAINED IN THIS CELL
291              !
292              iimin = 1
293              DO WHILE( G0%nav_lon(iimin,1) < Cell_lonmin ) 
294                 iimin = iimin + 1
295              ENDDO
296              !               
297              jjmin = 1
298              DO WHILE( G0%nav_lat(iimin,jjmin) < Cell_latmin ) 
299                 jjmin = jjmin + 1
300              ENDDO
301              !               
302              iimax = iimin 
303              DO WHILE( G0%nav_lon(iimax,1) <= Cell_lonmax ) 
304                 iimax = iimax + 1
305              iimax = MIN( iimax,SIZE(G0%bathy_meter,1))
306              ENDDO
307              !                               
308              jjmax = jjmin 
309              DO WHILE( G0%nav_lat(iimax,jjmax) <= Cell_latmax ) 
310                 jjmax = jjmax + 1
311              jjmax = MIN( jjmax,SIZE(G0%bathy_meter,2))
312
313              ENDDO
314              !
315              iimax = iimax-1
316              jjmax = jjmax-1
317              !               
318              iimin = MAX( iimin,1 )
319              jjmin = MAX( jjmin,1 )
320              iimax = MIN( iimax,SIZE(G0%bathy_meter,1))
321              jjmax = MIN( jjmax,SIZE(G0%bathy_meter,2))
322
323              nxhr = iimax - iimin + 1
324              nyhr = jjmax - jjmin + 1                   
325
326              IF( nxhr == 0 .OR. nyhr == 0 ) THEN
327                 trouble_points(ji,jj) = 1
328              ELSE
329
330                 ALLOCATE( vardep(nxhr,nyhr) )
331                 ALLOCATE( mask_oce(nxhr,nyhr) )
332                 mask_oce = 0       
333
334                 vardep(:,:) = G0%bathy_meter(iimin:iimax,jjmin:jjmax)
335
336                 WHERE( vardep(:,:) .GT. 0. )  mask_oce = 1
337
338                 IF( SUM(mask_oce) == 0 ) THEN
339                    G1%bathy_meter(ji,jj) = 0.
340                 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(:,:))
348                       ALLOCATE(vardep1d(nxhr*nyhr))
349                       vardep1d = RESHAPE(vardep,(/ nxhr*nyhr /) )
350                       CALL ssort(vardep1d,nxhr*nyhr)
351                       !
352                       ! Calculate median
353                       !
354                       IF (MOD(SUM(mask_oce),2) .NE. 0) THEN
355                          G1%bathy_meter(ji,jj) = vardep1d( SUM(mask_oce)/2 + 1)
356                       ELSE
357                          G1%bathy_meter(ji,jj) = ( vardep1d(SUM(mask_oce)/2) + vardep1d(SUM(mask_oce)/2+1) )/2.0
358                       END IF
359                       DEALLOCATE(vardep1d)       
360                    ENDIF
361                 ENDIF
362                 DEALLOCATE (mask_oce,vardep)
363
364              ENDIF
365           ENDDO
366        ENDDO
367
368        IF( SUM( trouble_points ) > 0 ) THEN
369           PRINT*,'too much empty cells, proceed to bilinear interpolation !!'
370           type_bathy_interp = 2
371        ENDIF
372
373     ENDIF
374
375     !
376     ! create logical array masksrc
377     !
378     IF( type_bathy_interp == 2) THEN 
379        !
380
381        !           
382        identical_grids = .FALSE.
383
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
388           IF( MAXVAL( ABS(G0%nav_lat(:,:)- G1%nav_lat(:,:)) ) < 0.0001 .AND.   &
389                MAXVAL( ABS(G0%nav_lon(:,:)- G1%nav_lon(:,:)) ) < 0.0001 ) THEN
390              G1%bathy_meter = G0%bathy_meter 
391              PRINT*,'same grid between ',elevation_database,' and child domain'   
392              identical_grids = .TRUE.                         
393           ENDIF
394        ENDIF
395
396
397        IF( .NOT. identical_grids ) THEN
398
399           ALLOCATE(masksrc(SIZE(G0%bathy_meter,1),SIZE(G0%bathy_meter,2)))
400           ALLOCATE(bathy_test(nxfin,nyfin))
401           !
402           !                    Where(G0%bathy_meter.le.0.00001)
403           !                   masksrc = .false.
404           !               ElseWhere
405           !
406           masksrc = .TRUE.
407           !
408           !               End where                       
409           !           
410           ! 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) 
417           !                                 
418           G1%bathy_meter = bathy_test               
419           !           
420           DEALLOCATE(masksrc)
421           DEALLOCATE(bathy_test) 
422
423        ENDIF
424        !           
425     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        !                 
441        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 = connectionsize*irafx + nbghostcellsfine                     
462        G1%gdepw_ps(1:boundary,:) = gdepw_ps_interp(1:boundary,:)
463        G1%gdepw_ps(:,1:boundary) = gdepw_ps_interp(:,1:boundary)
464        G1%gdepw_ps(nxfin-boundary+1:nxfin,:) = gdepw_ps_interp(nxfin-boundary+1:nxfin,:)
465        G1%gdepw_ps(:,nyfin-boundary+1:nyfin) = gdepw_ps_interp(:,nyfin-boundary+1:nyfin)
466
467
468        !                   
469
470        IF(.NOT. smoothing) WRITE(*,*) 'No smoothing process only connection is carried out'
471        WRITE(*,*) ' linear connection on ',nb_connection_pts,'coarse grid points'
472
473        connectionsize = 2 + nb_connection_pts 
474        !           
475        gdepw_ps_interp = 0.
476        CALL Check_interp(G0,gdepw_ps_interp)      ! interpolation in connection zone (3 coarse grid cells)
477        !
478        !
479        !
480        !
481        !                    LINEAR CONNECTION
482        !
483        !
484        !
485        !
486        !
487        IF (.NOT.ASSOCIATED(wtab)) &
488             ALLOCATE(wtab(SIZE(G1%bathy_meter,1),SIZE(G1%bathy_meter,2)))
489        wtab(:,:) = 0.
490        wghts = 1.
491        DO ji = boundary + 1 , boundary + nb_connection_pts * irafx
492            wghts = wghts - (1. / (nb_connection_pts*irafx - 1. ) )
493            DO jj=boundary+1,nyfin-boundary
494               wtab(ji,jj) = MAX(wghts, wtab(ji,jj)) 
495            END DO
496        END DO
497     
498        wghts = 1.
499        DO ji = nxfin - boundary , nxfin - boundary - nb_connection_pts * irafx +1 , -1
500            wghts = wghts - (1. / (nb_connection_pts*irafx - 1. ) )
501            DO jj=boundary+1,nyfin-boundary
502               wtab(ji,jj) = MAX(wghts, wtab(ji,jj))
503            END DO
504        END DO 
505
506        wghts = 1.
507        DO jj = boundary + 1 , boundary + nb_connection_pts * irafy
508            wghts = wghts - (1. / (nb_connection_pts*irafy - 1. ) )
509            DO ji=boundary+1,nxfin-boundary
510               wtab(ji,jj) = MAX(wghts, wtab(ji,jj))
511            END DO
512        END DO
513
514        wghts = 1.
515        DO jj = nyfin - boundary , nyfin - boundary - nb_connection_pts * irafy +1, -1
516            wghts = wghts - (1. / (nb_connection_pts*irafy - 1. ) )
517            DO ji=boundary+1,nxfin-boundary
518               wtab(ji,jj) = MAX(wghts, wtab(ji,jj))
519            END DO
520        END DO
521
522        G1%gdepw_ps(:,:) = (1.-wtab(:,:)) * G1%gdepw_ps(:,:) + gdepw_ps_interp(:,:)*wtab(:,:)
523
524        G1%bathy_meter = G1%gdepw_ps
525        !                     
526        connectionsize = 2 
527        !
528        IF(smoothing) THEN 
529
530           !
531           ! Smoothing to connect the connection zone (3 + nb_connection_pts coarse grid cells) and the interior domain
532           !
533           boundary = (connectionsize+nb_connection_pts)*irafx + nbghostcellsfine 
534           CALL smooth_topo(G1%gdepw_ps(boundary:nxfin-boundary+1,boundary:nyfin-boundary+1),nbiter)
535           G1%bathy_meter = G1%gdepw_ps                         
536        ENDIF
537
538
539        !
540        ! 
541        ! Remove closed seas
542        !                           
543        IF (removeclosedseas) THEN
544           ALLOCATE(bathy_test(nxfin,nyfin))
545           bathy_test=0.
546           WHERE (G1%bathy_meter(1,:).GT.0.)
547              bathy_test(1,:)=1
548           END WHERE
549           WHERE (G1%bathy_meter(nxfin,:).GT.0.)
550              bathy_test(nxfin,:)=1
551           END WHERE
552           WHERE (G1%bathy_meter(:,1).GT.0.)
553              bathy_test(:,1)=1
554           END WHERE
555           WHERE (G1%bathy_meter(:,nyfin).GT.0.)
556              bathy_test(:,nyfin)=1
557           END WHERE
558           nbadd = 1
559           DO WHILE (nbadd.NE.0)
560              nbadd = 0
561              DO j=2,nyfin-1
562                 DO i=2,nxfin-1
563                    IF (G1%bathy_meter(i,j).GT.0.) THEN
564                       IF (MAX(bathy_test(i,j+1),bathy_test(i,j-1), &
565                            bathy_test(i-1,j),bathy_test(i+1,j)).EQ.1) THEN
566                          IF (bathy_test(i,j).NE.1.) nbadd = nbadd + 1
567                          bathy_test(i,j)=1.
568                       ENDIF
569
570                    ENDIF
571                 ENDDO
572              ENDDO
573           ENDDO
574           WHERE (bathy_test.EQ.0.)
575              G1%bathy_meter = 0.
576           END WHERE
577           DEALLOCATE(bathy_test)
578        ENDIF
579        !
580        IF(bathy_update) CALL Update_Parent_Bathy( G0,G1 )                 
581        !
582        CALL set_child_name(parent_bathy_meter,child_ps)
583        status = Write_Bathy_meter(TRIM(child_ps),G1)       
584
585        IF(bathy_update) status = Write_Bathy_meter(TRIM(updated_parent_file),G0)
586
587        CALL get_partial_steps(G1)
588        !
589        G1%bathy_level=NINT(G1%bathy_level)
590        !
591        IF( TRIM(parent_meshmask_file) .NE. '/NULL' ) &
592             status = Write_Bathy_level(TRIM(Childlevel_file),G1)
593        !
594        WRITE(*,*) '****** Bathymetry successfully created for partial cells ******'
595        WRITE(*,*) ' '
596        !
597        STOP         
598     ENDIF
599     !           
600     !--------------------------------end if partial steps------------------------------------
601     !
602     !
603     status = Read_bathy_level(TRIM(parent_meshmask_file),G0)
604     !           
605     CALL levels_to_meter(G0)
606     !
607     ! compute constant bathymetry for Parent-Child bathymetry connection
608     !             
609     WHERE( G0%bathy_meter < 0.000001 ) G0%bathy_meter = 0.
610
611     CALL init_constant_bathy(G0%bathy_meter,bathy_fin_constant)
612     !
613     boundary = connectionsize*irafx + nbghostcellsfine   
614     !             
615     G1%bathy_meter(1:boundary,:) = bathy_fin_constant(1:boundary,:)
616     G1%bathy_meter(:,1:boundary) = bathy_fin_constant(:,1:boundary)
617     G1%bathy_meter(nxfin-boundary+1:nxfin,:) = bathy_fin_constant(nxfin-boundary+1:nxfin,:)
618     G1%bathy_meter(:,nyfin-boundary+1:nyfin) = bathy_fin_constant(:,nyfin-boundary+1:nyfin)
619     !
620     ! bathymetry smoothing
621     !                 
622     CALL smooth_topo(G1%bathy_meter(boundary:nxfin-boundary+1,boundary:nyfin-boundary+1),nbiter)
623     !
624     ! convert bathymetry from meters to levels
625     !
626     CALL meter_to_levels(G1) 
627     !           
628     DEALLOCATE(G1%bathy_meter)           
629     !             
630  ENDIF
631  !
632  !
633  ! make connection thanks to constant and interpolated bathymetry
634  !
635  !     
636  G1%bathy_level=NINT(G1%bathy_level)
637  !       
638  DO j=1,nyfin
639     DO i=1,nxfin
640        IF (g1%bathy_level(i,j).LT.0.) THEN
641           PRINT *,'error in ',i,j,g1%bathy_level(i,j)
642        ENDIF
643     ENDDO
644  ENDDO
645  !       
646  WHERE ((G1%bathy_level.LT.3.).AND.(G1%bathy_level.GT.0.))
647     G1%bathy_level=3
648  END WHERE
649  !
650  ! possibility to remove closed seas
651  !     
652  IF (removeclosedseas) THEN
653     ALLOCATE(bathy_test(nxfin,nyfin))
654
655     bathy_test=0.
656     WHERE (G1%bathy_level(1,:).GT.0.)
657        bathy_test(1,:)=1
658     END WHERE
659
660     WHERE (G1%bathy_level(nxfin,:).GT.0.)
661        bathy_test(nxfin,:)=1
662     END WHERE
663
664     WHERE (G1%bathy_level(:,1).GT.0.)
665        bathy_test(:,1)=1
666     END WHERE
667
668     WHERE (G1%bathy_level(:,nyfin).GT.0.)
669        bathy_test(:,nyfin)=1
670     END WHERE
671
672     nbadd = 1
673
674     DO WHILE (nbadd.NE.0)
675        nbadd = 0
676        DO j=2,nyfin-1
677           DO i=2,nxfin-1
678              IF (G1%bathy_level(i,j).GT.0.) THEN
679                 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
680                    IF (bathy_test(i,j).NE.1.) nbadd = nbadd + 1
681                    bathy_test(i,j)=1.
682                 ENDIF
683
684              ENDIF
685           ENDDO
686        ENDDO
687
688     ENDDO
689
690     WHERE (bathy_test.EQ.0.)
691        G1%bathy_level = 0.
692     END WHERE
693     DEALLOCATE(bathy_test)           
694  ENDIF
695
696
697  !
698  ! store interpolation result in output file
699  !
700  status = Write_Bathy_level(TRIM(Childlevel_file),G1)
701
702  WRITE(*,*) '****** Bathymetry successfully created for full cells ******'
703  WRITE(*,*) ' '
704
705  STOP
706END PROGRAM create_bathy
707
708
Note: See TracBrowser for help on using the repository browser.