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

Last change on this file since 12253 was 12253, checked in by clem, 10 months 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.