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 @ 9753

Last change on this file since 9753 was 9753, checked in by jchanut, 6 years ago

Modify arithmetic and median bathymetry average near land so that coastline is more realistic

  • Property svn:keywords set to Id
File size: 26.2 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  :: Cell_lonmin,Cell_lonmax,Cell_latmin,Cell_latmax,wghts
56  LOGICAL :: Pacifique=.FALSE.
57  INTEGER :: boundary,xpos,ypos,iimin,iimax,jjmax,jjmin
58  INTEGER :: nbloops,nxhr,nyhr,ji,jj,nbiter,nbloopmax
59  INTEGER :: ipt,jpt,iloc,jloc
60  INTEGER, DIMENSION(2) :: i_min,i_max,j_min,j_max
61
62  TYPE(Coordinates) :: G0,G1 
63  !     
64  narg = iargc()     
65  IF (narg == 0) THEN
66     namelistname = 'namelist.input'
67  ELSE
68     CALL getarg(1,namelistname)
69  ENDIF
70  !
71  ! read input file (namelist.input)
72  !
73  CALL read_namelist(namelistname)
74  !     
75  ! define names of child grid files
76  !
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     
96     !     
97     WRITE(*,*) 'First option'
98     !read coarse grid bathymetry and coordinates file
99     !           
100     WRITE(*,*) 'No new topo file ...'
101     status = Read_Coordinates(TRIM(parent_coordinate_file),G0)   
102     status = Read_bathy_level(TRIM(parent_meshmask_file),G0)
103     !           
104     IF( imax > SIZE(G0%glamt,1) .OR. jmax > SIZE(G0%glamt,2) .OR. &
105          imax <= imin .OR. jmax <= jmin ) THEN                   
106        WRITE(*,*) 'ERROR ***** bad child grid definition ...'
107        WRITE(*,*) 'please check imin,jmin,imax,jmax,jpizoom,jpjzoom values'       
108        WRITE(*,*) ' '
109        STOP
110     ENDIF
111     !
112     !read fine grid coordinates file
113     !     
114     status = Read_Coordinates(TRIM(Child_coordinates),G1,pacifique)
115     
116     IF( SIZE(G1%nav_lon,1) .NE. nxfin .OR. SIZE(G1%nav_lon,2) .NE. nyfin ) THEN
117        !
118        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(*,*) ' '
123        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) 
131     !           
132     ALLOCATE(masksrc(numlon,numlat))
133     !
134     ! create logical array masksrc
135     !
136     WHERE(G0%bathy_level.LE.0) 
137        masksrc = .FALSE.
138     ELSEWHERE
139        masksrc = .TRUE.
140     END WHERE
141
142     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     !             
154     ALLOCATE(G1%bathy_meter(nxfin,nyfin))
155     CALL levels_to_meter(G0)
156     !             
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     !           
166     DEALLOCATE(masksrc)
167     !-----------------                                     
168     !     
169     !
170     ! compute constant bathymetry for Parent-Child bathymetry connection
171     !             
172     CALL init_constant_bathy(G0%bathy_meter,bathy_fin_constant)
173     !
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     !                 
179     G1%bathy_meter(1:boundary,:) = bathy_fin_constant(1:boundary,:)
180     G1%bathy_meter(:,1:boundary) = bathy_fin_constant(:,1:boundary)
181     G1%bathy_meter(nxfin-boundary+1:nxfin,:) = bathy_fin_constant(nxfin-boundary+1:nxfin,:)
182     G1%bathy_meter(:,nyfin-boundary+1:nyfin) = bathy_fin_constant(:,nyfin-boundary+1:nyfin)
183     !                 
184     CALL smooth_topo(G1%bathy_meter(boundary:nxfin-boundary+1,boundary:nyfin-boundary+1),nbiter)
185     !             
186     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'
205     ! read fine and coarse grids coordinates file
206     !       
207     status = Read_Coordinates(TRIM(parent_coordinate_file),G0)
208     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                   
212        WRITE(*,*) 'ERROR ***** bad child grid definition ...'
213        WRITE(*,*) 'please check imin,jmin,imax,jmax,jpizoom,jpjzoom values'       
214        WRITE(*,*) ' '
215        STOP
216     ENDIF
217     !
218
219     
220     IF( SIZE(G1%nav_lon,1) .NE. nxfin .OR. SIZE(G1%nav_lon,2) .NE. nyfin ) THEN
221        !
222        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(*,*) ' '
227        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)
235        DEALLOCATE( G0%nav_lon, G0%nav_lat )
236        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)
240        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     
250     !               
251     ALLOCATE(G1%bathy_meter(nxfin,nyfin))
252     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        !
273        ALLOCATE(trouble_points(nxfin,nyfin))
274        trouble_points = 0
275        !
276        !  POINT DETECTION
277        !
278        !                       
279        DO jj = 2,numlatfine
280           DO ji = 2,numlonfine
281              !       
282              ! FINE GRID CELLS DEFINITION               
283              !
284              Cell_lonmin = MIN(G1%glamf(ji-1,jj-1),G1%glamf(ji,jj-1),G1%glamf(ji,jj),G1%glamf(ji-1,jj))
285              Cell_lonmax = MAX(G1%glamf(ji-1,jj-1),G1%glamf(ji,jj-1),G1%glamf(ji,jj),G1%glamf(ji-1,jj))
286              Cell_latmin = MIN(G1%gphif(ji-1,jj-1),G1%gphif(ji,jj-1),G1%gphif(ji,jj),G1%gphif(ji-1,jj))
287              Cell_latmax = MAX(G1%gphif(ji-1,jj-1),G1%gphif(ji,jj-1),G1%gphif(ji,jj),G1%gphif(ji-1,jj)) 
288              !               
289              ! SEARCH FOR ALL POINTS CONTAINED IN THIS CELL
290              !
291              iimin = 1
292              DO WHILE( G0%nav_lon(iimin,1) < Cell_lonmin ) 
293                 iimin = iimin + 1
294              ENDDO
295              !               
296              jjmin = 1
297              DO WHILE( G0%nav_lat(iimin,jjmin) < Cell_latmin ) 
298                 jjmin = jjmin + 1
299              ENDDO
300              !               
301              iimax = iimin 
302              DO WHILE( G0%nav_lon(iimax,1) <= Cell_lonmax ) 
303                 iimax = iimax + 1
304              iimax = MIN( iimax,SIZE(G0%bathy_meter,1))
305              ENDDO
306              !                               
307              jjmax = jjmin 
308              DO WHILE( G0%nav_lat(iimax,jjmax) <= Cell_latmax ) 
309                 jjmax = jjmax + 1
310              jjmax = MIN( jjmax,SIZE(G0%bathy_meter,2))
311
312              ENDDO
313              !
314              iimax = iimax-1
315              jjmax = jjmax-1
316              !               
317              iimin = MAX( iimin,1 )
318              jjmin = MAX( jjmin,1 )
319              iimax = MIN( iimax,SIZE(G0%bathy_meter,1))
320              jjmax = MIN( jjmax,SIZE(G0%bathy_meter,2))
321
322              nxhr = iimax - iimin + 1
323              nyhr = jjmax - jjmin + 1                   
324
325              IF( nxhr == 0 .OR. nyhr == 0 ) THEN
326                 trouble_points(ji,jj) = 1
327              ELSE
328
329                 ALLOCATE( vardep(nxhr,nyhr) )
330                 ALLOCATE( mask_oce(nxhr,nyhr) )
331                 mask_oce = 0       
332
333                 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
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 = 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
486        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
502        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        !
526        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     
536        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
564        END DO
565        IF (.NOT.identical_grids) THEN
566           G1%gdepw_ps(:,:) = (1.-G1%wgt(:,:)) * G1%gdepw_ps(:,:) + gdepw_ps_interp(:,:)*G1%wgt(:,:)
567        ENDIF
568
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
618                    ENDIF
619                 ENDDO
620              ENDDO
621           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        !
633        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)
637
638        CALL get_partial_steps(G1)
639        !
640        G1%bathy_level=NINT(G1%bathy_level)
641        !
642        IF( TRIM(parent_meshmask_file) .NE. '/NULL' ) &
643             status = Write_Bathy_level(TRIM(Childlevel_file),G1)
644        !
645        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     !             
681  ENDIF
682  !
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
757END PROGRAM create_bathy
758
759
Note: See TracBrowser for help on using the repository browser.