FMS  2025.04
Flexible Modeling System
topography.F90
1 !***********************************************************************
2 !* Apache License 2.0
3 !*
4 !* This file is part of the GFDL Flexible Modeling System (FMS).
5 !*
6 !* Licensed under the Apache License, Version 2.0 (the "License");
7 !* you may not use this file except in compliance with the License.
8 !* You may obtain a copy of the License at
9 !*
10 !* http://www.apache.org/licenses/LICENSE-2.0
11 !*
12 !* FMS is distributed in the hope that it will be useful, but WITHOUT
13 !* WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied;
14 !* without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
15 !* PARTICULAR PURPOSE. See the License for the specific language
16 !* governing permissions and limitations under the License.
17 !***********************************************************************
18 !> @defgroup topography_mod topography_mod
19 !> @ingroup topography
20 !> @brief Routines for creating land surface topography fields and land-water masks
21 !! for latitude-longitude grids.
22 !> @author Bruce Wyman
23 !!
24 !! This module generates realistic mountains and land-water masks
25 !! on a specified latitude-longitude grid by interpolating from the
26 !! 1/6 degree Navy mean topography and percent water data sets.
27 !! The fields that can be generated are mean and standard deviation
28 !! of topography within the specified grid boxes; and land-ocean (or
29 !! water) mask and land-ocean (or water) fractional area.
30 !!
31 !! The interpolation scheme conserves the area-weighted average
32 !! of the input data by using module horiz_interp.
33 !!
34 !! The interfaces get_gaussian_topog and gaussian_topog_init are documented
35 !! in @ref gaussian_topog_mod
36 
37 module topography_mod
38 
39 use gaussian_topog_mod, only: gaussian_topog_init, get_gaussian_topog
40 use horiz_interp_mod, only: horiz_interp_type, horiz_interp_new, &
42 
43 use fms_mod, only: check_nml_error, stdlog, &
44  mpp_pe, mpp_root_pe, write_version_number, &
45  error_mesg, fatal, note, &
46  mpp_error
47 ! required for fms2_io
48 use fms2_io_mod, only: read_data, fmsnetcdffile_t, file_exists, open_file
49 !-----------------------------------------------------------------------
50 
51 use constants_mod, only: pi
52 use mpp_mod, only: input_nml_file
53 use platform_mod, only: r4_kind, r8_kind, fms_path_len
54 
55 implicit none
56 private
57 
58 public :: topography_init, &
63 
64 !> @brief Returns a "realistic" mean surface height field.
65 !!
66 !> Returns realistic mountains on a latitude-longtude grid.
67 !! The returned field is the mean topography for the given grid boxes.
68 !! Computed using a conserving area-weighted interpolation.
69 !! The current input data set is the 1/6 degree Navy mean topography.
70 !!
71 !! @param blon The longitude (in radians) at grid box boundaries.
72 !! @param blat The latitude (in radians) at grid box boundaries.
73 !! @param zmean The mean surface height (meters). The size of this
74 !! field must be size(blon)-1 by size(blat)-1.
75 !! @return A logical value of TRUE is returned if the surface height field
76 !! was successfully created. A value of FALSE may be returned if the
77 !! input topography data set was not readable.
78 !!
79 !! <br>Example usage:
80 !! @code{.F90} flag = get_topog_mean ( blon, blat, zmean )@endcode
81 !> @ingroup topography_mod
82 
83 interface get_topog_mean
84  module procedure get_topog_mean_1d_r4, get_topog_mean_1d_r8
85  module procedure get_topog_mean_2d_r4, get_topog_mean_2d_r8
86 end interface get_topog_mean
87 
88 !> Returns a standard deviation of higher resolution topography with
89 !! the given model grid boxes.
90 !!
91 !> Returns the standard deviation of the "finer" input topography data set,
92 !! currently the Navy 1/6 degree mean topography data, within the
93 !! boundaries of the given input grid.
94 !!
95 !! @param blon The longitude (in radians) at grid box boundaries.
96 !! @param blat The latitude (in radians) at grid box boundaries.
97 !! @param [out] stdev The standard deviation of surface height (in meters) within
98 !! given input model grid boxes.
99 !! The size of this field must be size(blon)-1 by size(blat)-1.
100 !!
101 !! @return A logical value of TRUE is returned if the output field was
102 !! successfully created. A value of FALSE may be returned if the
103 !! input topography data set was not readable.
104 !!
105 !! Example usage:
106 !! @code{.F90} flag = get_topog_stdev( blon, blat, stdev ) @code
107 !> @ingroup topography_mod
109  module procedure get_topog_stdev_1d_r4, get_topog_stdev_1d_r8
110  module procedure get_topog_stdev_2d_r4, get_topog_stdev_2d_r8
111 end interface get_topog_stdev
112 
113 !> @brief Returns fractional area covered by ocean in a grid box.
114 !! Returns fractional area covered by ocean in the given model grid boxes.
115 !!
116 !! @param blon The longitude (in radians) at grid box boundaries.
117 !! @param blat The latitude (in radians) at grid box boundaries.
118 !! @param ocean_frac The fractional amount (0 to 1) of ocean in a grid box.
119 !! The size of this field must be size(blon)-1 by size(blat)-1.
120 !!
121 !! @return A logical value of TRUE is returned if the output field
122 !! was successfully created. A value of FALSE may be returned
123 !! if the Navy 1/6 degree percent water data set was not readable.
124 !!
125 !! Example usage:
126 !! @code{.F90} flag = get_ocean_frac ( blon, blat, ocean_frac ) @endcode
127 !> @ingroup topography_mod
128 interface get_ocean_frac
129  module procedure get_ocean_frac_1d_r4, get_ocean_frac_1d_r8
130  module procedure get_ocean_frac_2d_r4, get_ocean_frac_2d_r8
131 end interface get_ocean_frac
132 
133 !> @brief Returns a land-ocean mask in a grid box.
134 !!
135 !> Returns a land-ocean mask in the given model grid boxes.
136 !!
137 !! @param blon The longitude (in radians) at grid box boundaries.
138 !! @param blat The latitude (in radians) at grid box boundaries.
139 !! @param ocean_frac The fractional amount (0 to 1) of ocean in a grid box.
140 !! The size of this field must be size(blon)-1 by size(blat)-1.
141 !!
142 !! @return A logical value of TRUE is returned if the output field
143 !! was successfully created. A value of FALSE may be returned
144 !! if the Navy 1/6 degree percent water data set was not readable.
145 !!
146 !! Example code:
147 !! @code{.F90} flag = get_ocean_mask( blon, blat, ocean_mask ) @endcode
148 !> @ingroup topography_mod
149 interface get_ocean_mask
150  module procedure get_ocean_mask_1d_r4, get_ocean_mask_1d_r8
151  module procedure get_ocean_mask_2d_r4, get_ocean_mask_2d_r8
152 end interface get_ocean_mask
153 
154 !> @brief Returns fractional area covered by water.
155 !!
156 !> Returns the percent of water in a grid box.
157 !!
158 !! @param blon The longitude (in radians) at grid box boundaries.
159 !! @param blat The latitude (in radians) at grid box boundaries.
160 !! @param [out] water_frac The fractional amount (0 to 1) of water in a grid box.
161 !! The size of this field must be size(blon)-1 by size(blat)-1.
162 !!
163 !! @return A logical value of TRUE is returned if the output field
164 !! was successfully created. A value of FALSE may be returned
165 !! if the Navy 1/6 degree percent water data set was not readable.
166 !!
167 !! <br>Example usage:<br> @code{.F90} flag = get_water_frac ( blon, blat, water_frac ) @endcode
168 !> @ingroup topography_mod
169 interface get_water_frac
170  module procedure get_water_frac_1d_r4, get_water_frac_1d_r8
171  module procedure get_water_frac_2d_r4, get_water_frac_2d_r8
172 end interface get_water_frac
173 
174 !> @brief Returns a land-water mask in a grid box.
175 !!
176 !> Returns a land-water mask in the given model grid boxes.
177 !!
178 !! @param blon The longitude (in radians) at grid box boundaries.
179 !! @param blat The latitude (in radians) at grid box boundaries.
180 !! @param water_mask A binary mask for water (true) or land (false).
181 !! The size of this field must be size(blon)-1 by size(blat)-1.
182 !!
183 !! @return A logical value of TRUE is returned if the output field
184 !! was successfully created. A value of FALSE may be returned
185 !! if the Navy 1/6 degree percent water data set was not readable.
186 !!
187 !! Example usage: @code{.F90}flag = get_water_mask( blon, blat, water_mask ) @endcode
188 !> @ingroup topography_mod
189 interface get_water_mask
190  module procedure get_water_mask_1d_r4, get_water_mask_1d_r8
191  module procedure get_water_mask_2d_r4, get_water_mask_2d_r8
192 end interface get_water_mask
193 
194 interface interp_topog
195  module procedure interp_topog_1d_r4, interp_topog_1d_r8
196  module procedure interp_topog_2d_r4, interp_topog_2d_r8
197 end interface interp_topog
198 
199 interface find_indices
200  module procedure find_indices_r4, find_indices_r8
201 end interface find_indices
202 
203 interface input_data
204  module procedure input_data_r4, input_data_r8
205 end interface input_data
206 
207 interface interp_water
208  module procedure interp_water_1d_r4, interp_water_1d_r8
209  module procedure interp_water_2d_r4, interp_water_2d_r8
210 end interface interp_water
211 
213  module procedure determine_ocean_points_r4, determine_ocean_points_r8
214 end interface determine_ocean_points
215 
216 !> @addtogroup topography_mod
217 !> @{
218 
219 character(len=FMS_PATH_LEN) :: topog_file = 'DATA/navy_topography.data', &
220  water_file = 'DATA/navy_pctwater.data'
221 namelist /topography_nml/ topog_file, water_file
222 
223 integer, parameter :: TOPOG_INDEX = 1
224 integer, parameter :: WATER_INDEX = 2
225 logical :: file_is_opened(2) = .false.
226 type(fmsnetcdffile_t) :: fileobj(2) !< needed for fms2_io
227 
228 !-----------------------------------------------------------------------
229 ! --- resolution of the topography data set ---
230 ! <DATASET NAME="">
231 ! This module uses the 1/6 degree U.S. Navy mean topography
232 ! and percent water data sets.
233 !
234 ! These data sets have been re-formatted to separate 32-bit IEEE files.
235 ! The names of these files is specified by the <LINK SRC="#NAMELIST">namelist</LINK> input.
236 !
237 !The format for both files is as follows:
238 ! <PRE>
239 ! record = 1 nlon, nlat
240 ! record = 2 blon, blat
241 ! record = 3 data
242 ! </PRE>
243 !where:
244 ! <PRE>
245 ! nlon, nlat = The number of longitude and latitude points
246 ! in the horizontal grid. For the 1/6 degree
247 ! data sets this is 2160 x 1080. [integer]
248 ! blon, blat = The longitude and latitude grid box boundaries in degrees.
249 ! [real :: blon(nlon+1), blat(nlat+1)]
250 !
251 ! data = The topography or percent water data.
252 ! [real :: data(nlon,nlat)]
253 ! </PRE>
254 ! </DATASET>
255 integer :: ipts, jpts
256 integer, parameter :: compute_stdev = 123 ! use this flag to
257  ! compute st dev
258 
259 !-----------------------------------------------------------------------
260 
261 ! Include variable "version" to be written to log file.
262 #include<file_version.h>
263 
264 logical :: module_is_initialized = .false.
265 
266 !-----------------------------------------------------------------------
267 
268 contains
269 
270 !#######################################################################
271 
272 subroutine topography_init ()
273  if ( module_is_initialized ) return
274  call write_version_number("TOPOGRAPHY_MOD", version)
275  call read_namelist
276  module_is_initialized = .true.
277 end subroutine topography_init
278 
279 !#######################################################################
280 !#######################################################################
281 !################## private interfaces below here ##################
282 !#######################################################################
283 
284 function open_topog_file ( )
285 logical :: open_topog_file
286 real(kind=r4_kind) :: r_ipts, r_jpts
287 integer :: namelen
288 
289 namelen = len(trim(topog_file))
290  if ( file_exists(topog_file) .AND. topog_file(namelen-2:namelen) == '.nc') then
291  if (mpp_pe() == mpp_root_pe()) call mpp_error ('topography_mod', &
292  'Reading NetCDF formatted input data file: '//trim(topog_file), note)
293  if(.not. file_is_opened(topog_index) ) then
294  if(.not. open_file(fileobj(topog_index), topog_file, 'read' )) then
295  call mpp_error(fatal, 'topography_mod: Error in opening file '//trim(topog_file))
296  endif
297  endif
298 
299  call read_data(fileobj(topog_index), 'ipts', r_ipts)
300  call read_data(fileobj(topog_index), 'jpts', r_jpts)
301  ipts = nint(r_ipts)
302  jpts = nint(r_jpts)
303  open_topog_file = .true.
304  file_is_opened(topog_index) = .true.
305  else
306  open_topog_file = .false.
307  endif
308 
309 end function open_topog_file
310 
311 function open_water_file ( )
312 logical :: open_water_file
313 real(kind=r4_kind) :: r_ipts, r_jpts
314 integer :: namelen
315 
316 namelen = len(trim(water_file))
317 if ( file_exists(water_file) .AND. water_file(namelen-2:namelen) == '.nc') then
318  if (mpp_pe() == mpp_root_pe()) call mpp_error ('topography_mod', &
319  'Reading NetCDF formatted input data file: '//trim(water_file), note)
320  if(.not. file_is_opened(water_index) ) then
321  if(.not. open_file(fileobj(water_index), water_file, 'read' )) then
322  call mpp_error(fatal, 'topography_mod: Error in opening file '//trim(water_file))
323  endif
324  endif
325 
326  call read_data(fileobj(water_index), 'ipts', r_ipts)
327  call read_data(fileobj(water_index), 'jpts', r_jpts)
328  ipts = nint(r_ipts)
329  jpts = nint(r_jpts)
330  open_water_file = .true.
331  file_is_opened(water_index) = .true.
332  else
333  open_water_file = .false.
334  endif
335 
336 end function open_water_file
337 
338 
339 !#######################################################################
340 
341 !> @brief Reads the namelist file, write namelist to log file,
342 !! and initializes constants
343 subroutine read_namelist
344 
345 integer :: iunit, ierr, io
346 
347 ! read namelist
348 
349 read (input_nml_file, topography_nml, iostat=io)
350 ierr = check_nml_error(io,'topography_nml')
351 
352 ! write version and namelist to log file
353 
354 if (mpp_pe() == mpp_root_pe()) then
355  iunit = stdlog()
356  write (iunit, nml=topography_nml)
357 endif
358 
359 end subroutine read_namelist
360 
361 #include "topography_r4.fh"
362 #include "topography_r8.fh"
363 
364 end module topography_mod
365 
366 ! <INFO>
367 
368 ! <TESTPROGRAM NAME="">
369 !
370 ! To run this program you will need the topography and percent water
371 ! data sets and use the following namelist (in the input nml file).
372 !
373 ! &amp;gaussian_topog_nml
374 ! height = 5000., 3000., 3000., 3000.,
375 ! olon = 90., 255., 285., 0.,
376 ! olat = 45., 45., -15., -90.,
377 ! wlon = 15., 10., 5., 180.,
378 ! wlat = 15., 25., 25., 20., /
379 !
380 ! program test
381 !
382 ! ! test program for topography and gaussian_topog modules
383 ! <PRE>
384 ! use topography_mod
385 ! implicit none
386 !
387 ! integer, parameter :: nlon=24, nlat=18
388 ! real :: x(nlon), y(nlat), xb(nlon+1), yb(nlat+1), z(nlon,nlat)
389 ! real :: hpi, rtd
390 ! integer :: i,j
391 ! logical :: a
392 !
393 ! ! gaussian mountain parameters
394 ! real, parameter :: ht=4000.
395 ! real, parameter :: x0=90., y0=45. ! origin in degrees
396 ! real, parameter :: xw=15., yw=15. ! half-width in degees
397 ! real, parameter :: xr=30., yr= 0. ! ridge-width in degrees
398 !
399 ! ! create lat/lon grid in radians
400 ! hpi = acos(0.0)
401 ! rtd = 90./hpi ! rad to deg
402 ! do i=1,nlon
403 ! xb(i) = 4.*hpi*real(i-1)/real(nlon)
404 ! enddo
405 ! xb(nlon+1) = xb(1)+4.*hpi
406 ! yb(1) = -hpi
407 ! do j=2,nlat
408 ! yb(j) = yb(j-1) + 2.*hpi/real(nlat)
409 ! enddo
410 ! yb(nlat+1) = hpi
411 ! ! mid-point of grid boxes
412 ! x(1:nlon) = 0.5*(xb(1:nlon)+xb(2:nlon+1))
413 ! y(1:nlat) = 0.5*(yb(1:nlat)+yb(2:nlat+1))
414 ! ! test topography_mod routines
415 ! a = get_topog_mean(xb,yb,z)
416 ! call printz ('get_topog_mean')
417 !
418 ! a = get_water_frac(xb,yb,z)
419 ! z = z*100. ! in percent
420 ! call printz ('get_water_frac')
421 !
422 ! a = get_ocean_frac(xb,yb,z)
423 ! z = z*100. ! in percent
424 ! call printz ('get_ocean_frac')
425 !
426 ! ! test gaussian_topog_mod routines
427 ! a = .true.
428 ! z = get_gaussian_topog(x,y,ht,x0,y0,xw,yw,xr,yr)
429 ! call printz ('get_gaussian_topog')
430 !
431 ! call gaussian_topog_init (x,y,z)
432 ! call printz ('gaussian_topog_init')
433 !
434 ! contains
435 !
436 ! ! simple printout of topog/water array
437 ! subroutine printz (lab)
438 ! character(len=*), intent(in) :: lab
439 ! if (a) then
440 ! print '(/a)', trim(lab)
441 ! else
442 ! print '(/a)', 'no data available: '//trim(lab)
443 ! return
444 ! endif
445 ! ! print full grid
446 ! print '(3x,25i5)', (nint(x(i)*rtd),i=1,nlon)
447 ! do j=nlat,1,-1
448 ! print '(i3,25i5)', nint(y(j)*rtd), (nint(z(i,j)),i=1,nlon)
449 ! enddo
450 ! end subroutine printz
451 !
452 ! end program test
453 ! </PRE>
454 ! </TESTPROGRAM>
455 
456 ! <BUG>
457 ! Water mask produces some possible erroneous water points along
458 ! the coast of Antarctic (at about 90W).
459 ! </BUG>
460 
461 ! <FUTURE>Use of netcdf data sets. </FUTURE>
462 ! <FUTURE>Incorporate other topography and ocean data sets. </FUTURE>
463 !
464 ! </INFO>
465 !> @}
466 ! close documentation grouping
Opens a given netcdf or domain file.
Definition: fms2_io.F90:121
Read data from a defined field in a file.
Definition: fms2_io.F90:291
integer function, public check_nml_error(IOSTAT, NML_NAME)
Checks the iostat argument that is returned after reading a namelist and determines if the error code...
Definition: fms.F90:523
subroutine, public write_version_number(version, tag, unit)
Prints to the log file (or a specified unit) the version id string and tag name.
Definition: fms.F90:701
subroutine, public error_mesg(routine, message, level)
Print notes, warnings and error messages; terminates program for warning and error messages....
Definition: fms.F90:441
subroutine, public horiz_interp_del(Interp)
Deallocates memory used by "horiz_interp_type" variables. Must be called before reinitializing with h...
Subroutine for performing the horizontal interpolation between two grids.
integer function stdlog()
This function returns the current standard fortran unit numbers for log messages. Log messages,...
Definition: mpp_util.inc:58
integer function mpp_pe()
Returns processor ID.
Definition: mpp_util.inc:406
Error handler.
Definition: mpp.F90:381
subroutine read_namelist
Reads the namelist file, write namelist to log file, and initializes constants.
Definition: topography.F90:344
type(fmsnetcdffile_t), dimension(2) fileobj
needed for fms2_io
Definition: topography.F90:226
Returns fractional area covered by ocean in a grid box. Returns fractional area covered by ocean in t...
Definition: topography.F90:128
Returns a land-ocean mask in a grid box.
Definition: topography.F90:149
Returns a "realistic" mean surface height field.
Definition: topography.F90:83
Returns fractional area covered by water.
Definition: topography.F90:169
Returns a land-water mask in a grid box.
Definition: topography.F90:189
Allocates space and initializes a derived-type variable that contains pre-computed interpolation indi...
Returns a standard deviation of higher resolution topography with the given model grid boxes.
Definition: topography.F90:108