-
Notifications
You must be signed in to change notification settings - Fork 15
/
Copy pathprogram_setup.F90
319 lines (273 loc) · 13.7 KB
/
program_setup.F90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
!> @file
!! @brief Set up program execution.
!! @author Larissa Reames CIWRO/NOAA/NSSL/FRDD
!> This module contains code to read the setup namelist file.
!!
!! @author Larissa Reames CIWRO/NOAA/NSSL/FRDD
module program_setup
use esmf
use ESMF_LogPublicMod
use constants_module
use utils_mod
use misc_definitions_module
implicit none
private
type(ESMF_LogKind_Flag), public :: LogType
! Namelist variables
character(len=500), public :: grid_file_input_grid = "NULL" !< Full path of MPAS file containing grid information
character(len=500), public :: diag_file_input_grid = "NULL" !< Full path of input diagnostic MPAS data
character(len=500), public :: hist_file_input_grid = "NULL" !< Full path of input history MPAS data
character(len=500), public :: file_target_grid = "NULL" !<Full path of file containing target
!<grid information for target_grid_type='file'
character(len=500), public :: output_file = "NULL" !< Full path of output file
logical, public :: interp_diag = .false. !< Read data from diag file?
logical, public :: interp_hist = .false. !< Read data from hist file?
logical, public :: wrf_mod_vars = .false.!< Whether to modify variable values/dimensions
!< to conform to WRF format. Set to true for
!< UPP-compatible output
character(len=500), public :: target_grid_type !< Grid type to interpolate data to
!< Valid options: 'file', 'lambert',
!< 'mercator','polar',lat-lon'
character(len=500), public :: block_decomp_file = "NULL" !< Full path to MPAS grid-specific block decomposition file
real(kind=4), public :: missing_value = -9999. !< Value to be used for mapped data outside of input data
!! These entries are only valid when target_grid_type is not 'file'
logical, public :: is_regional = .true. !< Is the output grid regional or global?
!< Default: True
integer, public :: i_target !< # staggered east-west grid points in target grid
integer, public :: j_target !< # staggered north-south grid points in target grid
real, public :: truelat1 = NAN !< First true latitude (all projections)
real, public :: truelat2 = NAN !< Second true latitude (LCC only)
real, public :: stand_lon = NAN !< Longitude parallel to y-axis (-180->180E)
real, public :: dx = NAN !< Grid cell east-west dimension(meters or deg for
!< target_grid_type='lat-lon')
real, public :: dy = NAN !< Grid cell north-south dimension(meters or deg for
!< target_grid_type='lat-lon')
real, public :: ref_lat !< Latitude of reference point
real, public :: ref_lon !< Longitude of reference point
real, public :: ref_x !< Grid-relative e-w index of reference point
!< Defaults to grid center (nx/2)
real, public :: ref_y !< Grid-relative n-s index of reference point
!< Defaults to grid center (ny/2)
!These aren't namelist variables but they're created directly from them
real, public :: dxkm !< grid-cell east-west dimension (meters)
real, public :: dykm !< grid-cell north-south dimension (meters)
real, public :: dlondeg !< grid-cell east-west dimension (deg)
real, public :: dlatdeg !< grid-cell north-south dimension (deg)
real, public :: known_lat !< Latitude of reference point
real, public :: known_lon !< Longitude of reference point
real, public :: known_x !< Grid-relative e-w index of reference point
real, public :: known_y !< Grid-relative n-s index of reference point
real, public :: pole_lat = -90.0 !< Latitude of pole for target grid projection
real, public :: pole_lon = 0.0 !< Longitude of pole for target grid projection
real, public :: cone !< Projection cone (valid only for proj_code=PROJ_LC)
real, public :: hemi !< Projection hemisphere (valid only for proj_code=PROJ_LC or
!< PROJ_CASSINI)
integer, public :: proj_code !< Integer code corresponding to the requested
!< target grid projection type
character(len=500), public :: map_proj_char !< Map projection name
logical, public :: interp_as_bundle = .true. !< If true, use ESMF_FieldBundleRegridStore and interpolate the fields as a bundle
!< If false, use ESMFRegridStore and interpolate individual fields.
!< .false. seems faster and less memory intensive
!< Currently, only applies to conservative regridding
public :: read_setup_namelist
contains
!> Reads program configuration namelist.
!!
!! @oaran unum unit number to use for namelist file (default 41)
!! @param filename the name of the configuration file (defaults to
!! ./fort.41).
!! @author Larissa Reames CIWRO/NOAA/NSSL/FRDD
subroutine read_setup_namelist(unum, filename)
implicit none
character(len=*), intent(in), optional :: filename
integer, intent(in), optional :: unum
character(:), allocatable :: filename_to_use
character(len=500) :: map_proj
integer :: unit_to_use
logical :: esmf_log,decomp_exists
integer :: is, ie, ierr
!Namelist variables that are used to create global variables
integer :: nx,ny
namelist /config/ grid_file_input_grid, diag_file_input_grid, hist_file_input_grid, &
file_target_grid, output_file, interp_diag, interp_hist, &
wrf_mod_vars, esmf_log,target_grid_type,nx,ny,dx,dy,ref_lat,ref_lon,ref_x,ref_y,&
truelat1,truelat2,stand_lon,is_regional,pole_lat,pole_lon, interp_as_bundle,block_decomp_file, &
missing_value
ref_x = NAN
ref_y = NAN
ref_lat = NAN
ref_lon = NAN
dx = NAN
dy = NAN
! pole_lat = 90.0
! pole_lon = 0.0
nx = 0
ny = 0
!print*,"- READ SETUP NAMELIST"
if (present(filename)) then
filename_to_use = filename
else
filename_to_use = "./fort.41"
endif
if (present(unum)) then
unit_to_use = unum
else
unit_to_use = 41
endif
open(unit_to_use, file=filename_to_use, iostat=ierr)
if (ierr /= 0) call error_handler("OPENING SETUP NAMELIST.", ierr)
read(unit_to_use, nml=config, iostat=ierr)
if (ierr /= 0) call error_handler("READING SETUP NAMELIST.", ierr)
close (unit_to_use)
if (esmf_log) then
LogType = ESMF_LOGKIND_MULTI_ON_ERROR
else
LogType = ESMF_LOGKIND_NONE
endif
!if (block_decomp_file=='NULL') then
!call error_handler("block_decomp_file IS REQUIRED BUT IS MISSING IN NAMELIST.", -1)
!else
if(.not. block_decomp_file=='NULL') then
inquire(file=block_decomp_file, exist=decomp_exists)
if (.not. decomp_exists) then
call error_handler("block_decomp_file DOES NOT EXIST.",-1)
endif
endif
if (trim(target_grid_type) .ne. 'file') then
dxkm = dx
dykm = dy
known_lat = ref_lat
known_lon = ref_lon
known_x = ref_x
known_y = ref_y
i_target = nx-1
j_target = ny-1
map_proj = to_upper(target_grid_type)
!print*, map_proj
! Assign parameters to module variables
if ((index(map_proj, 'LAMBERT') /= 0) .and. &
(len_trim(map_proj) == len('LAMBERT'))) then
proj_code = PROJ_LC
map_proj_char = 'Lambert Conformal'
IF (truelat1 .LT. 0.) THEN
hemi = -1.0
ELSE
hemi = 1.0
ENDIF
else if ((index(map_proj, 'MERCATOR') /= 0) .and. &
(len_trim(map_proj) == len('MERCATOR'))) then
proj_code = PROJ_MERC
map_proj_char = 'Mercator'
else if ((index(map_proj, 'POLAR') /= 0) .and. &
(len_trim(map_proj) == len('POLAR'))) then
proj_code = PROJ_PS
map_proj_char = "Polar Stereographic"
else if ((index(map_proj, 'LAT-LON') /= 0) .and. &
(len_trim(map_proj) == len('LAT-LON'))) then
proj_code = PROJ_LATLON
map_proj_char = 'Lat/Lon'
else if ((index(map_proj, 'ROTATED LAT-LON') /= 0) .and. &
(len_trim(map_proj) == len('ROTATED LAT-LON'))) then
proj_code = PROJ_CASSINI
map_proj_char = 'Cylindrical Equidistant'
if (known_lat > 0.) then
pole_lat = 90. - known_lat
pole_lon = 180.
stand_lon = -known_lon
endif
dlondeg = dx
dlatdeg = dy
IF (truelat1 .LT. 0.) THEN
hemi = -1.0
ELSE
hemi = 1.0
ENDIF
else
call error_handler('In namelist, invalid target_grid_type specified. Valid '// &
'projections are "lambert", "mercator", "polar", '// &
'"lat-lon", and "rotated lat-lon".',ERROR_CODE)
end if
if (proj_code == PROJ_LATLON) then
! If no dx,dy specified, assume global grid
if (dx == NAN .and. dy == NAN) then
if (is_regional) then
call error_handler("For lat-lon projection, if dx/dy are not specified "// &
"a global grid is assumed. Please set dx/dy if a regional grid is desired, "//&
"or change is_regional to .false. if a global grid is desired.", ERROR_CODE)
endif
dlondeg = 360. / (i_target) ! Here, we really do not want e_we-s_we+1
dlatdeg = 180. / (j_target) ! Here, we really do not want e_we-s_we+1
known_x = 1.
known_y = 1.
known_lon = stand_lon + dlondeg/2.
known_lat = -90. + dlatdeg/2.
dxkm = EARTH_RADIUS_M * PI * 2.0 / (i_target)
dykm = EARTH_RADIUS_M * PI / (j_target)
! If dx,dy specified, however, assume regional grid
else
if (.not. is_regional) then
call error_handler("For lat-lon projection, if dx/dy are specified "// &
"a regional grid is assumed. Please unset dx/dy if a global grid is desired, "//&
"or change is_regional to .true. if a regional grid is desired.", ERROR_CODE)
endif
dlatdeg = dy
dlondeg = dx
dxkm = dlondeg * EARTH_RADIUS_M * PI * 2.0 / 360.0
dykm = dlatdeg * EARTH_RADIUS_M * PI * 2.0 / 360.0
!print*, "dxkm, dykm = ", dxkm, dykm
if (known_lat == NAN .or. known_lon == NAN) then
call error_handler('For lat-lon projection, if dx/dy are specified, '// &
'a regional domain is assumed, and a ref_lat,ref_lon must also be specified',ERROR_CODE)
end if
end if
end if
! Manually set truelat2 = truelat1 if truelat2 not specified for Lambert
if (proj_code==PROJ_LC) then
if (truelat2 == NAN) then
if (truelat1 == NAN) call error_handler("No TRUELAT1 specified for Lambert conformal projection.",ERROR_CODE)
truelat2 = truelat1
endif
call lc_cone(truelat1, truelat2, cone)
endif
! If the user hasn't supplied a known_x and known_y, assume the center of domain 1
if (known_x == NAN .and. known_y == NAN) then
known_x = real(i_target+1) / 2.
known_y = real(j_target+1) / 2.
! print*, known_x, known_y
else if (known_x == NAN .or. known_y == NAN) then
call error_handler('In namelist, neither or both of ref_x, ref_y must be specified.',ERROR_CODE)
end if
endif
return
end subroutine read_setup_namelist
!> Compute the cone factor of a Lambert Conformal projection
!!
!! @param truelat1 Projection true latitude 1
!! @param truelat2 Projection true latitude 2
!! @param cone LCC Projection cone (output)
!! @author Larissa Reames CIWRO/NOAA/NSSL/FRDD
SUBROUTINE lc_cone(truelat1, truelat2, cone)
IMPLICIT NONE
! Input Args
REAL, INTENT(IN) :: truelat1 ! (-90 -> 90 degrees N)
REAL, INTENT(IN) :: truelat2 ! " " " " "
! Output Args
REAL, INTENT(OUT) :: cone
! Locals
! BEGIN CODE
! First, see if this is a secant or tangent projection. For tangent
! projections, truelat1 = truelat2 and the cone is tangent to the
! Earth's surface at this latitude. For secant projections, the cone
! intersects the Earth's surface at each of the distinctly different
! latitudes
IF (ABS(truelat1-truelat2) .GT. 0.1) THEN
cone = ALOG10(COS(truelat1*rad_per_deg)) - &
ALOG10(COS(truelat2*rad_per_deg))
cone = cone /(ALOG10(TAN((45.0 - ABS(truelat1)/2.0) * rad_per_deg)) - &
ALOG10(TAN((45.0 - ABS(truelat2)/2.0) * rad_per_deg)))
ELSE
cone = SIN(ABS(truelat1)*rad_per_deg )
ENDIF
RETURN
END SUBROUTINE lc_cone
end module program_setup