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

Last change on this file since 12253 was 12253, checked in by clem, 4 years ago

1) resolve some conflicts when using partial steps or not. 2) Make sure grids remain identical when they should. 3) in case of a double zoom, it is no more necessary to have the bilinear interpolation option when updating the 1st parent grid. Btw, I know these comments are unclear but it is very difficult to explain...

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