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

Last change on this file since 10381 was 10381, checked in by clem, 23 months ago

attempt to correct several bugs in the NESTING tools. With this version, domcfg.nc file should be written correctly and the partial steps should be taken into account.

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