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