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