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

Last change on this file since 10025 was 10025, checked in by clem, 2 years ago

nesting tools are partly rewritten (mostly for create_coordinates and bathy) to get better functionality. Now you can use the nesting to either define an agrif zoom or a regional domain (for bdy purposes). Also, the nesting tools output a domain_cfg.nc that can be directly used in NEMO4 without the need of DOMAINcfg tool. The option of median average for bathymetry interpolation still does not work properly but it's not new

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