FMS 2025.01-dev
Flexible Modeling System
Loading...
Searching...
No Matches
topography.inc
1!***********************************************************************
2!* GNU Lesser General Public License
3!*
4!* This file is part of the GFDL Flexible Modeling System (FMS).
5!*
6!* FMS is free software: you can redistribute it and/or modify it under
7!* the terms of the GNU Lesser General Public License as published by
8!* the Free Software Foundation, either version 3 of the License, or (at
9!* your option) any later version.
10!*
11!* FMS is distributed in the hope that it will be useful, but WITHOUT
12!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13!* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
14!* for more details.
15!*
16!* You should have received a copy of the GNU Lesser General Public
17!* License along with FMS. If not, see <http://www.gnu.org/licenses/>.
18!***********************************************************************
19!> @defgroup topography_mod topography_mod
20!> @ingroup topography
21!> @brief Routines for creating land surface topography fields and land-water masks
22!! for latitude-longitude grids.
23!> @author Bruce Wyman
24!!
25!! This module generates realistic mountains and land-water masks
26!! on a specified latitude-longitude grid by interpolating from the
27!! 1/6 degree Navy mean topography and percent water data sets.
28!! The fields that can be generated are mean and standard deviation
29!! of topography within the specified grid boxes; and land-ocean (or
30!! water) mask and land-ocean (or water) fractional area.
31!!
32!! The interpolation scheme conserves the area-weighted average
33!! of the input data by using module horiz_interp.
34!!
35!! The interfaces get_gaussian_topog and gaussian_topog_init are documented
36!! in @ref gaussian_topog_mod
37
38!#######################################################################
39
40 !> @brief Returns a "realistic" mean surface height field.
41 !!
42 !> Returns realistic mountains on a latitude-longtude grid.
43 !! The returned field is the mean topography for the given grid boxes.
44 !! Computed using a conserving area-weighted interpolation.
45 !! The current input data set is the 1/6 degree Navy mean topography.
46 !!
47 !! @returns A logical value of true is returned if the surface height field was successfully
48 !! created. A value of false may be returned if the input topography data set was not readable.
49 !!
50 !! @throws FATAL, shape(zmean) is not equal to (/size(blon)-1,size(blat)-1/)
51 !! Check the input grid size and output field size.
52
53 function get_topog_mean_1d_(blon, blat, zmean)
54
55 real(kind=fms_top_kind_), intent(in), dimension(:) :: blon !< Longitude (radians) at grid box boundaries
56 real(kind=fms_top_kind_), intent(in), dimension(:) :: blat !< Latitude (radians) at grid box boundaries
57 real(kind=fms_top_kind_), intent(out), dimension(:,:) :: zmean !< Mean surface height(meters). Size must be
58 !! size(blon)-1 by size(blat)-1
59 logical :: GET_TOPOG_MEAN_1D_
60 integer, parameter :: lkind = fms_top_kind_
61
62!-----------------------------------------------------------------------
63 if (.not. module_is_initialized) call topography_init()
64
65 if ( any(shape(zmean(:,:)) /= (/size(blon(:))-1,size(blat(:))-1/)) ) &
66 call error_mesg('get_topog_mean_1d','shape(zmean) is not&
67 & equal to (/size(blon)-1,size(blat)-1/))', fatal)
68
69 get_topog_mean_1d_ = open_topog_file()
70
71 if (get_topog_mean_1d_) call interp_topog(blon, blat, zmean)
72
73!-----------------------------------------------------------------------
74
75 end function get_topog_mean_1d_
76
77!############################################################
78
79 function get_topog_mean_2d_(blon, blat, zmean)
80
81 real(kind=fms_top_kind_), intent(in), dimension(:,:) :: blon, blat
82 real(kind=fms_top_kind_), intent(out), dimension(:,:) :: zmean
83 logical :: GET_TOPOG_MEAN_2D_
84 integer, parameter :: lkind = fms_top_kind_
85!-----------------------------------------------------------------------
86 if (.not. module_is_initialized) call topography_init()
87
88 if ( any(shape(zmean(:,:)) /= (/size(blon,1)-1,size(blon,2)-1/)) .or. &
89 any(shape(zmean(:,:)) /= (/size(blat,1)-1,size(blat,2)-1/)) ) &
90 call error_mesg('get_topog_mean_2d','shape(zmean) is not&
91 & equal to (/size(blon,1)-1,size(blon,2)-1/))', fatal)
92
93 get_topog_mean_2d_ = open_topog_file()
94
95 if (get_topog_mean_2d_) call interp_topog(blon, blat, zmean)
96!-----------------------------------------------------------------------
97
98 end function get_topog_mean_2d_
99
100!#######################################################################
101
102 !> @brief Returns a standard deviation of higher resolution topography with
103 !! the given model grid boxes.
104 !!
105 !> Returns the standard deviation of the "finer" input topography data set,
106 !! currently the Navy 1/6 degree mean topography data, within the
107 !! boundaries of the given input grid.
108 !!
109 !! @returns A logical value of true if the output field was successfully created and false
110 !! if the input topography data set was not readable.
111 function get_topog_stdev_1d_(blon, blat, stdev)
112
113 real(kind=fms_top_kind_), intent(in), dimension(:) :: blon !< Longitude (radians) at grid box boundaries
114 real(kind=fms_top_kind_), intent(in), dimension(:) :: blat !< Latitude (radians) at grid box boundaries
115 real(kind=fms_top_kind_), intent(out), dimension(:,:) :: stdev !< The standard deviation of surface height (in
116 !! meters) within given input model grid boxes. Size must be
117 !! size(blon)-1 by size(blat)-1
118 logical :: GET_TOPOG_STDEV_1D_
119 integer, parameter :: lkind = fms_top_kind_
120!-----------------------------------------------------------------------
121 if (.not. module_is_initialized) call topography_init()
122
123 if ( any(shape(stdev(:,:)) /= (/size(blon(:))-1,size(blat(:))-1/)) ) &
124 call error_mesg('get_topog_stdev','shape(stdev) is not&
125 & equal to (/size(blon)-1,size(blat)-1/))', fatal)
126
127 get_topog_stdev_1d_ = open_topog_file()
128
129 if (get_topog_stdev_1d_) call interp_topog(blon, blat, &
130 stdev, flag=compute_stdev)
131
132!-----------------------------------------------------------------------
133
134 end function get_topog_stdev_1d_
135
136!#######################################################################
137
138 function get_topog_stdev_2d_(blon, blat, stdev)
139
140 real(FMS_TOP_KIND_), intent(in), dimension(:,:) :: blon, blat
141 real(FMS_TOP_KIND_), intent(out), dimension(:,:) :: stdev
142 logical :: GET_TOPOG_STDEV_2D_
143 integer, parameter :: lkind = fms_top_kind_
144!-----------------------------------------------------------------------
145 if (.not. module_is_initialized) call topography_init()
146
147 if ( any(shape(stdev(:,:)) /= (/size(blon,1)-1,size(blon,2)-1/)) .or. &
148 any(shape(stdev(:,:)) /= (/size(blat,1)-1,size(blat,2)-1/)) ) &
149 call error_mesg('get_topog_stdev_2d','shape(stdev) is not&
150 & equal to (/size(blon,1)-1,size(blon,2)-1/))', fatal)
151
152 get_topog_stdev_2d_ = open_topog_file()
153
154 if (get_topog_stdev_2d_) call interp_topog(blon, blat, &
155 stdev, flag=compute_stdev)
156!-----------------------------------------------------------------------
157
158 end function get_topog_stdev_2d_
159
160!#######################################################################
161
162 !> @brief Returns fractional area covered by ocean in a grid box.
163 !> @returns A logical value of true if the output field was successfully created. A value of false
164 !! may be returned if the Navy 1/6 degree percent water data set was not readable.
165 function get_ocean_frac_1d_(blon, blat, ocean_frac)
166
167 real(FMS_TOP_KIND_), intent(in), dimension(:) :: blon !< Longitude (radians) at grid box boundaries
168 real(FMS_TOP_KIND_), intent(in), dimension(:) :: blat !< Latitude (radians) at grid box boundaries
169 real(FMS_TOP_KIND_), intent(out), dimension(:,:) :: ocean_frac !< The fractional amount (0-1) of ocean in a grid
170 !! box. The size must be size(blon)-1 by size(blat)-1
171 logical :: GET_OCEAN_FRAC_1D_
172 integer, parameter :: lkind = fms_top_kind_
173!-----------------------------------------------------------------------
174 if (.not. module_is_initialized) call topography_init()
175
176 if ( any(shape(ocean_frac(:,:)) /= (/size(blon(:))-1,size(blat(:))-1/)) ) &
177 call error_mesg('get_ocean_frac','shape(ocean_frac) is not&
178 & equal to (/size(blon)-1,size(blat)-1/))', fatal)
179
180 get_ocean_frac_1d_ = open_water_file()
181 if(get_ocean_frac_1d_) call interp_water( blon, blat, &
182 ocean_frac, do_ocean=.true. )
183
184!-----------------------------------------------------------------------
185
186 end function get_ocean_frac_1d_
187
188!#######################################################################
189
190 function get_ocean_frac_2d_(blon, blat, ocean_frac)
191
192 real(kind=fms_top_kind_), intent(in), dimension(:,:) :: blon, blat
193 real(kind=fms_top_kind_), intent(out), dimension(:,:) :: ocean_frac
194 logical :: GET_OCEAN_FRAC_2D_
195 integer, parameter :: lkind = fms_top_kind_
196!-----------------------------------------------------------------------
197 if (.not. module_is_initialized) call topography_init()
198
199 if ( any(shape(ocean_frac(:,:)) /= (/size(blon,1)-1,size(blon,2)-1/)) .or. &
200 any(shape(ocean_frac(:,:)) /= (/size(blat,1)-1,size(blat,2)-1/)) ) &
201 call error_mesg('get_ocean_frac_2d','shape(ocean_frac) is not&
202 & equal to (/size(blon,1)-1,size(blon,2)-1/))', fatal)
203
204 get_ocean_frac_2d_ = open_water_file()
205 if(get_ocean_frac_2d_) call interp_water( blon, blat, &
206 ocean_frac, do_ocean=.true. )
207
208!-----------------------------------------------------------------------
209
210 end function get_ocean_frac_2d_
211
212!#######################################################################
213
214 !> @brief Returns a land-ocean mask in a grid box.
215 !> @returns A logical value of true if the output field was successfully created. A value of false
216 !! may be returned if the Navy 1/6 degree percent water data set was not readable.
217 function get_ocean_mask_1d_(blon, blat, ocean_mask)
218
219 real(kind=fms_top_kind_), intent(in), dimension(:) :: blon !< Longitude (radians) at grid box boundaries
220 real(kind=fms_top_kind_), intent(in), dimension(:) :: blat !< Latitude (radians) at grid box boundaries
221 logical, intent(out), dimension(:,:) :: ocean_mask !< Mask for ocean in a grid box.
222 !! The size must be size(blon)-1 by size(blat)-1
223 logical :: GET_OCEAN_MASK_1D_
224 real(kind=fms_top_kind_), dimension(size(ocean_mask,1),size(ocean_mask,2)) :: ocean_frac
225 integer, parameter :: lkind = fms_top_kind_
226!-----------------------------------------------------------------------
227 if (.not. module_is_initialized) call topography_init()
228
229 if (get_ocean_frac(blon, blat, ocean_frac) ) then
230 where (ocean_frac > 0.50_lkind)
231 ocean_mask = .true.
232 elsewhere
233 ocean_mask = .false.
234 end where
235 get_ocean_mask_1d_ = .true.
236 else
237 get_ocean_mask_1d_ = .false.
238 endif
239!-----------------------------------------------------------------------
240
241 end function get_ocean_mask_1d_
242
243!#######################################################################
244
245 function get_ocean_mask_2d_(blon, blat, ocean_mask)
246
247 real(kind=fms_top_kind_), intent(in), dimension(:,:) :: blon, blat
248 logical, intent(out), dimension(:,:) :: ocean_mask
249 logical :: GET_OCEAN_MASK_2D_
250 real(kind=fms_top_kind_), dimension(size(ocean_mask,1),size(ocean_mask,2)) :: ocean_frac
251 integer, parameter :: lkind = fms_top_kind_
252!-----------------------------------------------------------------------
253 if (.not. module_is_initialized) call topography_init()
254
255 if (get_ocean_frac(blon, blat, ocean_frac) ) then
256 where (ocean_frac > 0.50_lkind)
257 ocean_mask = .true.
258 elsewhere
259 ocean_mask = .false.
260 end where
261 get_ocean_mask_2d_ = .true.
262 else
263 get_ocean_mask_2d_ = .false.
264 endif
265
266!-----------------------------------------------------------------------
267
268 end function get_ocean_mask_2d_
269
270 !> @brief Returns the percent of water in a grid box.
271 !> @returns A logical value of true if the output field was successfully created. A value of false
272 !! may be returned if the Navy 1/6 degree percent water data set was not readable.
273 !!
274 !! @throws FATAL, shape(water_frac) is not equal to (/size(blon)-1,size(blat)-1/)
275 !! Check the input grid size and output field size.
276 function get_water_frac_1d_(blon, blat, water_frac)
277 real(kind=fms_top_kind_), intent(in), dimension(:) :: blon !< The longitude (in radians) at grid box boundaries.
278 real(kind=fms_top_kind_), intent(in), dimension(:) :: blat !< The latitude (in radians) at grid box boundaries.
279 real(kind=fms_top_kind_), intent(out), dimension(:,:) :: water_frac !< The fractional amount (0 to 1) of water in a
280 !! grid box. The size of this field must be size(blon)-1 by size(blat)-1.
281 logical :: GET_WATER_FRAC_1D_
282 integer, parameter :: lkind = fms_top_kind_
283
284!-----------------------------------------------------------------------
285 if (.not. module_is_initialized) call topography_init()
286
287 if ( any(shape(water_frac(:,:)) /= (/size(blon(:))-1,size(blat(:))-1/)) ) &
288 call error_mesg('get_water_frac_1d','shape(water_frac) is not&
289 & equal to (/size(blon)-1,size(blat)-1/))', fatal)
290
291 get_water_frac_1d_ = open_water_file()
292 if(get_water_frac_1d_) call interp_water( blon, blat, water_frac )
293
294!-----------------------------------------------------------------------
295
296 end function get_water_frac_1d_
297
298!#######################################################################
299
300 function get_water_frac_2d_(blon, blat, water_frac)
301
302 real(kind=fms_top_kind_), intent(in), dimension(:,:) :: blon !< The longitude (in radians) at grid box boundaries.
303 real(kind=fms_top_kind_), intent(in), dimension(:,:) :: blat !< The latitude (in radians) at grid box boundaries.
304 real(kind=fms_top_kind_), intent(out), dimension(:,:) :: water_frac !< The fractional amount (0 to 1) of water in a
305 !! grid box. The size of this field must be size(blon)-1 by size(blat)-1.
306 logical :: GET_WATER_FRAC_2D_
307 integer, parameter :: lkind = fms_top_kind_
308
309!-----------------------------------------------------------------------
310 if (.not. module_is_initialized) call topography_init()
311
312 if ( any(shape(water_frac(:,:)) /= (/size(blon,1)-1,size(blon,2)-1/)) .or. &
313 any(shape(water_frac(:,:)) /= (/size(blat,1)-1,size(blat,2)-1/)) ) &
314 call error_mesg('get_water_frac_2d','shape(water_frac) is not&
315 & equal to (/size(blon,1)-1,size(blon,2)-1/))', fatal)
316
317 get_water_frac_2d_ = open_water_file()
318 if(get_water_frac_2d_) call interp_water( blon, blat, water_frac )
319
320!-----------------------------------------------------------------------
321
322 end function get_water_frac_2d_
323
324!#######################################################################
325
326 !> @brief Returns a land-water mask in the given model grid boxes.
327 !> @return A logical value of TRUE is returned if the output field
328 !! was successfully created. A value of FALSE may be returned
329 !! if the Navy 1/6 degree percent water data set was not readable.
330 function get_water_mask_1d_(blon, blat, water_mask)
331
332 real(kind=fms_top_kind_), intent(in), dimension(:) :: blon !< The longitude (in radians) at grid box boundaries.
333 real(kind=fms_top_kind_), intent(in), dimension(:) :: blat !< The latitude (in radians) at grid box boundaries.
334 logical, intent(out), dimension(:,:) :: water_mask !< A binary mask for water (true) or land (false).
335 !! The size of this field must be size(blon)-1 by size(blat)-1.
336 logical :: GET_WATER_MASK_1D_
337
338 real(kind=fms_top_kind_), dimension(size(water_mask,1),size(water_mask,2)) :: water_frac
339 integer, parameter :: lkind = fms_top_kind_
340!-----------------------------------------------------------------------
341 if (.not. module_is_initialized) call topography_init()
342
343 if (get_water_frac(blon, blat, water_frac) ) then
344 where (water_frac > 0.50_lkind)
345 water_mask = .true.
346 elsewhere
347 water_mask = .false.
348 end where
349 get_water_mask_1d_ = .true.
350 else
351 get_water_mask_1d_ = .false.
352 endif
353!-----------------------------------------------------------------------
354
355 end function get_water_mask_1d_
356
357!#######################################################################
358
359 function get_water_mask_2d_(blon, blat, water_mask)
360
361 real(kind=fms_top_kind_), intent(in), dimension(:,:) :: blon !< The longitude (in radians) at grid box boundaries.
362 real(kind=fms_top_kind_), intent(in), dimension(:,:) :: blat !< The latitude (in radians) at grid box boundaries.
363 logical, intent(out), dimension(:,:) :: water_mask !< A binary mask for water (true) or land (false).
364 !! The size of this field must be size(blon)-1 by size(blat)-1.
365 logical :: GET_WATER_MASK_2D_
366 real(kind=fms_top_kind_), dimension(size(water_mask,1),size(water_mask,2)) :: water_frac
367 integer, parameter :: lkind = fms_top_kind_
368!-----------------------------------------------------------------------
369 if (.not. module_is_initialized) call topography_init()
370
371 if (get_water_frac(blon, blat, water_frac) ) then
372 where (water_frac > 0.50_lkind)
373 water_mask = .true.
374 elsewhere
375 water_mask = .false.
376 end where
377 get_water_mask_2d_ = .true.
378 else
379 get_water_mask_2d_ = .false.
380 endif
381
382!-----------------------------------------------------------------------
383
384 end function get_water_mask_2d_
385
386!#######################################################################
387
388subroutine interp_topog_1d_( blon, blat, zout, flag)
389real(kind=fms_top_kind_), intent(in) :: blon(:), blat(:)
390real(kind=fms_top_kind_), intent(out) :: zout(:,:)
391integer, intent(in), optional :: flag
392
393real(kind=fms_top_kind_) :: xdat(ipts+1), ydat(jpts+1)
394real(kind=fms_top_kind_) :: zdat(ipts,jpts)
395real(kind=fms_top_kind_) :: zout2(size(zout,1),size(zout,2))
396integer, parameter :: lkind = fms_top_kind_
397
398 call input_data(topog_index, xdat, ydat, zdat)
399
400 call horiz_interp( zdat, xdat, ydat, blon, blat, zout )
401
402! compute standard deviation if necessary.
403 if (present(flag)) then
404 if (flag == compute_stdev) then
405 zdat = zdat*zdat
406 call horiz_interp ( zdat, xdat, ydat, blon, blat, zout2 )
407 zout = zout2 - zout*zout
408 where (zout > 0.0_lkind)
409 zout = sqrt( zout )
410 elsewhere
411 zout = 0.0_lkind
412 endwhere
413 endif
414 endif
415
416end subroutine interp_topog_1d_
417
418!#######################################################################
419
420subroutine interp_topog_2d_( blon, blat, zout, flag )
421real(kind=fms_top_kind_), intent(in) :: blon(:,:), blat(:,:)
422real(kind=fms_top_kind_), intent(out) :: zout(:,:)
423integer, intent(in), optional :: flag
424
425real(kind=fms_top_kind_) :: xdat(ipts+1), ydat(jpts+1)
426real(kind=fms_top_kind_) :: zdat(ipts,jpts)
427real(kind=fms_top_kind_) :: zout2(size(zout,1),size(zout,2))
428integer, parameter :: lkind = fms_top_kind_
429integer :: js, je
430type (horiz_interp_type) :: Interp
431
432 call input_data(topog_index, xdat, ydat, zdat)
433
434 call find_indices( minval(blat), maxval(blat), ydat, js, je )
435
436 call horiz_interp_new ( interp, xdat, ydat(js:je+1), blon, blat )
437 call horiz_interp ( interp, zdat(:,js:je), zout )
438
439! compute standard deviation if necessary
440 if (present(flag)) then
441 if (flag == compute_stdev) then
442 zdat = zdat*zdat
443 call horiz_interp ( interp, zdat(:,js:je), zout2 )
444 zout = zout2 - zout*zout
445 where (zout > 0.0_lkind)
446 zout = sqrt( zout )
447 elsewhere
448 zout = 0.0_lkind
449 endwhere
450 endif
451 endif
452
453 call horiz_interp_del( interp )
454
455end subroutine interp_topog_2d_
456
457!#######################################################################
458
459subroutine find_indices_( ybeg, yend, ydat, js, je )
460real(kind=fms_top_kind_), intent(in) :: ybeg, yend, ydat(:)
461integer, intent(out) :: js, je
462integer :: j
463integer, parameter :: lkind = fms_top_kind_
464
465 js = 1
466 do j = 1, size(ydat(:))-1
467 if (ybeg >= ydat(j) .and. ybeg <= ydat(j+1)) then
468 js = j
469 exit
470 endif
471 enddo
472
473 je = size(ydat(:))-1
474 do j = js, size(ydat(:))-1
475 if (yend >= ydat(j) .and. yend <= ydat(j+1)) then
476 je = j
477 exit
478 endif
479 enddo
480
481!print '(a,i2,2(a,f10.5),2(a,i4))', "PE=",mpp_pe()," phs=",ybeg," phn=",yend," js=",js," je=",je
482
483end subroutine find_indices_
484
485!#######################################################################
486subroutine input_data_( indx, xdat, ydat, zdat )
487integer, intent(in) :: indx
488real(kind=fms_top_kind_), intent(out) :: xdat(ipts+1), ydat(jpts+1), zdat(ipts,jpts)
489
490if( file_is_opened(indx) ) then
491 call read_data(fileobj(indx), 'xdat', xdat)
492 call read_data(fileobj(indx), 'ydat', ydat)
493 call read_data(fileobj(indx), 'zdat', zdat)
494endif
495
496end subroutine input_data_
497
498!#######################################################################
499
500subroutine interp_water_1d_( blon, blat, zout, do_ocean )
501real(kind=fms_top_kind_), intent(in) :: blon(:), blat(:)
502real(kind=fms_top_kind_), intent(out) :: zout(:,:)
503logical, intent(in), optional :: do_ocean
504real(kind=fms_top_kind_) :: xdat(ipts+1), ydat(jpts+1), zdat(ipts,jpts)
505 call input_data ( water_index, xdat, ydat, zdat )
506
507! only use designated ocean points
508if (present(do_ocean)) then
509 if (do_ocean) call determine_ocean_points (zdat)
510endif
511
512! interpolate onto output grid
513call horiz_interp ( zdat, xdat, ydat, blon, blat, zout )
514
515end subroutine interp_water_1d_
516
517!#######################################################################
518
519subroutine interp_water_2d_( blon, blat, zout, do_ocean )
520real(kind=fms_top_kind_), intent(in) :: blon(:,:), blat(:,:)
521real(kind=fms_top_kind_), intent(out) :: zout(:,:)
522logical, intent(in), optional :: do_ocean
523real(kind=fms_top_kind_) :: xdat(ipts+1), ydat(jpts+1), zdat(ipts,jpts)
524
525call input_data ( water_index, xdat, ydat, zdat )
526
527! only use designated ocean points
528if (present(do_ocean)) then
529 if (do_ocean) call determine_ocean_points (zdat)
530endif
531
532! interpolate onto output grid
533call horiz_interp ( zdat, xdat, ydat, blon, blat, zout )
534
535end subroutine interp_water_2d_
536
537!#######################################################################
538
539subroutine determine_ocean_points_(pctwater)
540real(kind=fms_top_kind_), intent(inout) :: pctwater(:,:)
541logical :: ocean(size(pctwater,1),size(pctwater,2))
542integer :: i, j, m, n, im, ip, jm, jp, new
543integer, parameter :: lkind = fms_top_kind_
544real(kind=fms_top_kind_) :: ocean_pct_crit = 0.500_lkind
545
546
547! resolution of the grid
548m = size(pctwater,1)
549n = size(pctwater,2)
550
551! the 1/6 degree navy percent water data set
552! designates ocean grid boxes as 100 percent water
553! all other grid boxes have <= 99 percent water
554
555! set a mask for ocean grid boxes
556ocean = (pctwater > 0.999_lkind)
557new = count(ocean)
558
559! set land grid boxes that have sufficient amount of water
560! to ocean grid boxes when they are adjacent to ocean points
561! iterate until there are no new ocean points
562do
563 if (new == 0) exit
564 new = 0
565
566 do j = 1, n
567 do i = 1, m
568 if (.not.ocean(i,j) .and. pctwater(i,j) > ocean_pct_crit) then
569 im = i-1; ip = i+1; jm = j-1; jp = j+1
570 if (im == 0) im = m
571 if (ip == m+1) ip = 1
572 if (jm == 0) jm = 1
573 if (jp == n+1) jp = n
574 ! check the 8 grid boxes that surround this grid box
575 if (ocean(im,j ) .or. ocean(ip,j ) .or. ocean(i ,jm) .or. ocean(i ,jp) .or. &
576 ocean(im,jm) .or. ocean(ip,jm) .or. ocean(ip,jp) .or. ocean(im,jp)) then
577 ocean(i,j) = .true.
578 new = new + 1
579 endif
580 endif
581 enddo
582 enddo
583 !print *, 'new=',new
584
585enddo
586
587! final step is to elimate water percentage if land
588where (.not.ocean) pctwater = 0.0_lkind
589
590end subroutine determine_ocean_points_
591
592!> @}
593! close documentation grouping