orog_mask_tools 1.14.0
Loading...
Searching...
No Matches
inland.F90
Go to the documentation of this file.
1
4
9PROGRAM inland_mask
10 USE cs_nb
11 IMPLICIT NONE
12
13 INTEGER :: tile, i, j
14 TYPE(nb_gp_idx) :: nbs
15
16 REAL, ALLOCATABLE :: inland(:,:,:)
17 REAL, ALLOCATABLE :: land_frac(:,:,:)
18 INTEGER :: tile_beg, tile_end
19 INTEGER :: cs_res, x_res, y_res
20 CHARACTER(len=32) :: arg
21 INTEGER :: stat
22 INTEGER :: max_rd
23 REAL :: cutoff
24 CHARACTER(len=1) :: reg
25
26 LOGICAL, ALLOCATABLE :: done(:,:,:)
27
28 CALL getarg(0, arg) ! get the program name
29 IF (iargc() /= 3 .AND. iargc() /= 4) THEN
30 print*, 'Usage: '
31 print*, trim(arg), ' [cres(48,96, ...)][nonland cutoff][max recursive depth][regional(r),global(g)]'
32 stop
33 ENDIF
34
35 CALL getarg(1, arg)
36 READ(arg,*,iostat=stat) cs_res
37 CALL getarg(2, arg)
38 READ(arg,*,iostat=stat) cutoff
39 CALL getarg(3, arg)
40 READ(arg,*,iostat=stat) max_rd
41 CALL getarg(4, reg)
42
43 IF (reg /= 'r') THEN
44 tile_beg = 1; tile_end = 6
45! read in orography data (land_frac)
46 CALL read_orog(cs_res)
47
48! init inter-panel neighbor index
49 CALL idx_init(cs_res)
50
51! create a inland mask
52 CALL mark_global_inland(cs_res)
53
54! write back to the orography data files
55 CALL write_inland(cs_res)
56 ELSE
57! read in orography data (land_frac) for the SAR domain
58 CALL read_orog_reg(cs_res)
59
60! init inter-panel neighbor index
61 CALL idx_init_reg(x_res, y_res)
62
63! create a inland mask
64 CALL mark_inland_reg(cs_res)
65
66! write back to the orography data files
67 CALL write_inland_reg(cs_res)
68 ENDIF
69
70 CALL free_mem()
71
72CONTAINS
73
79SUBROUTINE mark_global_inland(cs_res)
80 INTEGER, INTENT(IN) :: cs_res
81 INTEGER :: i_seed, j_seed
82
83 ALLOCATE(done(cs_res,cs_res,6))
84 ALLOCATE(inland(cs_res,cs_res,6))
85 done = .false.
86 inland = 1.0
87
88 i_seed = cs_res/2; j_seed = cs_res/2
89 CALL mark_global_inland_rec_d(i_seed, j_seed, 2, 0)
90
91! to make sure black sea is excluded
92 i_seed = real(cs_res)/32.0*3; j_seed = i_seed
93 CALL mark_global_inland_rec_d(i_seed, j_seed, 3, 0)
94
95 DEALLOCATE(done)
96
97END SUBROUTINE mark_global_inland
98
104SUBROUTINE mark_inland_reg(cs_res)
105 INTEGER, INTENT(IN) :: cs_res
106 INTEGER :: i_seed, j_seed
107 INTEGER :: i
108
109 ALLOCATE(done(x_res,y_res,1))
110 ALLOCATE(inland(x_res,y_res,1))
111
112 done = .false.
113 inland = 1.0
114
115! set up three seeds, for the current domain
116 i_seed = 1; j_seed = y_res/2
117 CALL mark_regional_inland_rec_d(i_seed, j_seed, 1, 0)
118
119 i_seed = x_res; j_seed = y_res/2
120 CALL mark_regional_inland_rec_d(i_seed, j_seed, 1, 0)
121
122 i_seed = x_res/3; j_seed = 1
123 CALL mark_regional_inland_rec_d(i_seed, j_seed, 1, 0)
124
125 j_seed = 1
126 DO i = 1, x_res
127 CALL mark_regional_inland_rec_d(i, j_seed, 1, 0)
128 ENDDO
129
130 j_seed = y_res
131 DO i = x_res/2, x_res
132 CALL mark_regional_inland_rec_d(i, j_seed, 1, 0)
133 ENDDO
134
135! set up additional 3 seeds for ESG CONUS grid
136! i_seed = 1600; j_seed = 1040
137! i_seed = x_res - 10; j_seed = y_res
138! CALL mark_regional_inland_rec_d(i_seed, j_seed, 1, 0)
139
140! i_seed = x_res - 60; j_seed = y_res
141! CALL mark_regional_inland_rec_d(i_seed, j_seed, 1, 0)
142
143! i_seed = x_res - 275; j_seed = y_res
144! CALL mark_regional_inland_rec_d(i_seed, j_seed, 1, 0)
145
146! i_seed = 500; j_seed = 1
147! CALL mark_regional_inland_rec_d(i_seed, j_seed, 1, 0)
148
149 DEALLOCATE(done)
150
151END SUBROUTINE mark_inland_reg
152
161RECURSIVE SUBROUTINE mark_global_inland_rec_d(i, j, t, rd)
162 INTEGER, INTENT(IN) :: i, j, t, rd
163
164 TYPE(nb_gp_idx) :: nbs
165 INTEGER :: k, nrd
166
167 IF (t == 0) RETURN
168
169 IF (land_frac(i,j,t) <= 0.15) THEN
170 nrd = 1
171 ELSE
172 nrd = rd + 1
173 ENDIF
174
175 IF (nrd > max_rd) RETURN
176
177 IF (done(i,j,t)) RETURN
178 IF (land_frac(i,j,t) < cutoff) THEN
179 done(i,j,t) = .true.
180 inland(i,j,t) = 0.0
181 CALL neighbors(t, i, j, nbs)
182 ! recursively go through k neighbors
183 DO k = 1, 4
184 CALL mark_global_inland_rec_d(nbs%ijt(1,k),nbs%ijt(2,k),nbs%ijt(3,k),nrd)
185 ENDDO
186 ENDIF
187
188END SUBROUTINE mark_global_inland_rec_d
189
198RECURSIVE SUBROUTINE mark_regional_inland_rec_d(i, j, t, rd)
199 INTEGER, INTENT(IN) :: i, j, t, rd
200
201 TYPE(nb_gp_idx) :: nbs
202 INTEGER :: k, nrd
203
204 IF (t == 0) RETURN
205
206 IF (land_frac(i,j,t) <= 0.15) THEN
207 nrd = 1
208 ELSE
209 nrd = rd + 1
210 ENDIF
211
212 IF (nrd > max_rd) RETURN
213
214 IF (done(i,j,t)) RETURN
215 IF (land_frac(i,j,t) < cutoff) THEN
216 done(i,j,t) = .true.
217 inland(i,j,t) = 0.0
218 CALL neighbors_reg(i, j, nbs)
219 ! recursively go through k neighbors
220 DO k = 1, 4
221 CALL mark_regional_inland_rec_d(nbs%ijt(1,k),nbs%ijt(2,k),nbs%ijt(3,k),nrd)
222 ENDDO
223 ENDIF
224
225END SUBROUTINE mark_regional_inland_rec_d
226
232SUBROUTINE read_orog(cs_res)
233 USE netcdf
234 INTEGER, INTENT(IN) :: cs_res
235
236 INTEGER :: tile_sz, tile_num
237 INTEGER :: stat, ncid, varid
238 INTEGER :: land_frac_id, slmsk_id, geolon_id, geolat_id
239 CHARACTER(len=256) :: filename,string
240 CHARACTER(len=1) :: ich
241 CHARACTER(len=4) res_ch
242 CHARACTER(len=8) dimname
243
244 INTEGER :: i, j
245 REAL, ALLOCATABLE :: var_tmp(:,:)
246
247 ALLOCATE(var_tmp(cs_res,cs_res))
248 ALLOCATE(land_frac(cs_res,cs_res,6))
249
250 WRITE(res_ch,'(I4)') cs_res
251 DO tile_num = tile_beg, tile_end
252 WRITE(ich, '(I1)') tile_num
253 filename = "oro.C" // trim(adjustl(res_ch)) // ".tile" // ich // ".nc"
254 print *,'Read, update, and write ',trim(filename)
255 stat = nf90_open(filename, nf90_nowrite, ncid)
256 CALL nc_opchk(stat, "nf90_open oro_data.nc")
257
258! original orodata netcdf file uses (y, x) order, so we made change to match it.
259 stat = nf90_inq_varid(ncid, "land_frac", land_frac_id)
260 CALL nc_opchk(stat, "nf90_inq_varid: land_frac")
261 stat = nf90_get_var(ncid, land_frac_id, var_tmp, &
262 start = (/ 1, 1 /), count = (/ cs_res, cs_res /) )
263 CALL nc_opchk(stat, "nf90_get_var: land_frac")
264 land_frac(:,:,tile_num) = var_tmp(:,:)
265 stat = nf90_close(ncid)
266 CALL nc_opchk(stat, "nf90_close oro_data.nc")
267 ENDDO
268
269 DEALLOCATE(var_tmp)
270
271END SUBROUTINE read_orog
272
278SUBROUTINE read_orog_reg(cs_res)
279 USE netcdf
280 INTEGER, INTENT(IN) :: cs_res
281
282 INTEGER :: tile_num
283 INTEGER :: stat, ncid, x_dimid, y_dimid, varid
284 INTEGER :: land_frac_id, slmsk_id, geolon_id, geolat_id
285 CHARACTER(len=256) :: filename,string
286 CHARACTER(len=1) :: ich
287 CHARACTER(len=4) res_ch
288 CHARACTER(len=8) dimname
289
290 INTEGER :: i, j
291 REAL, ALLOCATABLE :: var_tmp(:,:)
292
293 WRITE(res_ch,'(I4)') cs_res
294 tile_num = 7
295 WRITE(ich, '(I1)') tile_num
296 filename = "oro.C" // trim(adjustl(res_ch)) // ".tile" // ich // ".nc"
297 print *,'Read, update, and write ',trim(filename)
298 stat = nf90_open(filename, nf90_nowrite, ncid)
299 CALL nc_opchk(stat, "nf90_open oro_data.nc")
300 stat = nf90_inq_dimid(ncid, "lon", x_dimid)
301 CALL nc_opchk(stat, "nf90_inq_dim: x")
302 stat = nf90_inq_dimid(ncid, "lat", y_dimid)
303 CALL nc_opchk(stat, "nf90_inq_dim: y")
304 stat = nf90_inquire_dimension(ncid,x_dimid,dimname,len=x_res)
305 CALL nc_opchk(stat,'nf90_inquire_dimension: lon')
306 stat = nf90_inquire_dimension(ncid,y_dimid,dimname,len=y_res)
307 CALL nc_opchk(stat,'nf90_inquire_dimension: lon')
308
309! original orodata netcdf file uses (y, x) order, so we made change to match it.
310 ALLOCATE(var_tmp(x_res,y_res))
311 ALLOCATE(land_frac(x_res,y_res,1))
312 stat = nf90_inq_varid(ncid, "land_frac", land_frac_id)
313 CALL nc_opchk(stat, "nf90_inq_varid: land_frac")
314 stat = nf90_get_var(ncid, land_frac_id, var_tmp, &
315 start = (/ 1, 1 /), count = (/ x_res, y_res /) )
316 CALL nc_opchk(stat, "nf90_get_var: land_frac")
317 land_frac(:,:,1) = var_tmp(:,:)
318 stat = nf90_close(ncid)
319 CALL nc_opchk(stat, "nf90_close oro_data.nc")
320
321 DEALLOCATE(var_tmp)
322
323END SUBROUTINE read_orog_reg
324
330SUBROUTINE write_inland(cs_res)
331 USE netcdf
332 INTEGER, INTENT(IN) :: cs_res
333
334 CHARACTER(len=256) :: filename
335 CHARACTER(len=1) :: ich
336 CHARACTER(len=4) res_ch
337
338 INTEGER :: tile_num
339 INTEGER :: stat, ncid, x_dimid, y_dimid, inland_id, dimids(2)
340 REAL, ALLOCATABLE :: var_tmp(:,:)
341
342 ALLOCATE(var_tmp(cs_res,cs_res))
343
344 WRITE(res_ch,'(I4)') cs_res
345 DO tile_num = tile_beg, tile_end
346 WRITE(ich, '(I1)') tile_num
347 filename = "oro.C" // trim(adjustl(res_ch)) // ".tile" // ich // ".nc"
348 print *,'write inland to ',trim(filename)
349 stat = nf90_open(filename, nf90_write, ncid)
350 CALL nc_opchk(stat, "nf90_open oro_data.nc")
351 stat = nf90_inq_dimid(ncid, "lon", x_dimid)
352 CALL nc_opchk(stat, "nf90_inq_dim: x")
353 stat = nf90_inq_dimid(ncid, "lat", y_dimid)
354 CALL nc_opchk(stat, "nf90_inq_dim: y")
355
356! original orodata netcdf file uses (y, x) order, so we made change to match it.
357 dimids = (/ x_dimid, y_dimid /)
358
359! define a new variables
360 stat = nf90_inq_varid(ncid,"inland",inland_id) !safeguard if "inland" exists
361 IF (stat /= nf90_noerr) THEN
362 stat = nf90_redef(ncid)
363 CALL nc_opchk(stat, "nf90_redef")
364 stat = nf90_def_var(ncid,"inland",nf90_float,dimids,inland_id)
365 CALL nc_opchk(stat, "nf90_def_var: inland")
366 stat = nf90_put_att(ncid, inland_id,'coordinates','geolon geolat')
367 CALL nc_opchk(stat, "nf90_put_att: inland:coordinates")
368 stat = nf90_put_att(ncid, inland_id,'description', &
369 'inland = 1 indicates grid cells away from coast')
370 CALL nc_opchk(stat, "nf90_put_att: inland:description")
371 stat = nf90_enddef(ncid)
372 CALL nc_opchk(stat, "nf90_enddef")
373 ENDIF
374
375 var_tmp(:,:) = inland(:,:,tile_num)
376 stat = nf90_put_var(ncid, inland_id, var_tmp, &
377 start = (/ 1, 1 /), count = (/ cs_res, cs_res /) )
378 CALL nc_opchk(stat, "nf90_put_var: inland")
379
380 stat = nf90_close(ncid)
381 CALL nc_opchk(stat, "nf90_close oro_data.nc")
382 ENDDO
383 DEALLOCATE(var_tmp)
384
385END SUBROUTINE write_inland
386
392SUBROUTINE write_inland_reg(cs_res)
393 USE netcdf
394 INTEGER, INTENT(IN) :: cs_res
395
396 CHARACTER(len=256) :: filename
397 CHARACTER(len=1) :: ich
398 CHARACTER(len=4) res_ch
399
400 INTEGER :: tile_num
401 INTEGER :: stat, ncid, x_dimid, y_dimid, inland_id, dimids(2)
402 REAL, ALLOCATABLE :: var_tmp(:,:)
403 CHARACTER(len=8) dimname
404
405 ALLOCATE(var_tmp(x_res,y_res))
406
407 WRITE(res_ch,'(I4)') cs_res
408 tile_num = 7
409 WRITE(ich, '(I1)') tile_num
410 filename = "oro.C" // trim(adjustl(res_ch)) // ".tile" // ich // ".nc"
411 print*,'write inland to ',trim(filename)
412 stat = nf90_open(filename, nf90_write, ncid)
413 CALL nc_opchk(stat, "nf90_open oro_data.nc")
414 stat = nf90_inq_dimid(ncid, "lon", x_dimid)
415 CALL nc_opchk(stat, "nf90_inq_dim: x")
416 stat = nf90_inq_dimid(ncid, "lat", y_dimid)
417 CALL nc_opchk(stat, "nf90_inq_dim: y")
418 stat = nf90_inquire_dimension(ncid,x_dimid,dimname,len=x_res)
419 CALL nc_opchk(stat,'nf90_inquire_dimension: lon')
420 stat = nf90_inquire_dimension(ncid,y_dimid,dimname,len=y_res)
421 CALL nc_opchk(stat,'nf90_inquire_dimension: lon')
422
423! original orodata netcdf file uses (y, x) order, so we made change to match it.
424 dimids = (/ x_dimid, y_dimid /)
425
426! define a new variables
427 stat = nf90_inq_varid(ncid,"inland",inland_id) !safeguard if "inland" exists
428 IF (stat /= nf90_noerr) THEN
429 stat = nf90_redef(ncid)
430 CALL nc_opchk(stat, "nf90_redef")
431 stat = nf90_def_var(ncid,"inland",nf90_float,dimids,inland_id)
432 CALL nc_opchk(stat, "nf90_def_var: inland")
433 stat = nf90_put_att(ncid, inland_id,'coordinates','geolon geolat')
434 CALL nc_opchk(stat, "nf90_put_att: inland:coordinates")
435 stat = nf90_put_att(ncid, inland_id,'description', &
436 'inland = 1 indicates grid cells away from coast')
437 CALL nc_opchk(stat, "nf90_put_att: inland:description")
438 stat = nf90_enddef(ncid)
439 CALL nc_opchk(stat, "nf90_enddef")
440 ENDIF
441
442 var_tmp(:,:) = inland(:,:,1)
443 stat = nf90_put_var(ncid, inland_id, var_tmp, &
444 start = (/ 1, 1 /), count = (/ x_res, y_res /) )
445 CALL nc_opchk(stat, "nf90_put_var: inland")
446
447 stat = nf90_close(ncid)
448 CALL nc_opchk(stat, "nf90_close oro_data.nc")
449
450 DEALLOCATE(var_tmp)
451
452END SUBROUTINE write_inland_reg
453
455SUBROUTINE free_mem()
456
457 DEALLOCATE(inland, land_frac)
458
459END SUBROUTINE free_mem
460
467SUBROUTINE nc_opchk(stat,opname)
468 USE netcdf
469 IMPLICIT NONE
470 INTEGER stat
471 CHARACTER(len=*) opname
472 CHARACTER(64) msg
473
474 IF (stat .NE.0) THEN
475 msg = trim(opname) // ' Error, status code and message:'
476 print*,trim(msg), stat, nf90_strerror(stat)
477 stop
478 END IF
479
480END SUBROUTINE nc_opchk
481
482END PROGRAM inland_mask
483
read_orog
subroutine read_orog(cs_res)
Read in orography (land fraction) data.
Definition inland.F90:233
mark_global_inland
subroutine mark_global_inland(cs_res)
Create inland mask for global grid.
Definition inland.F90:80
mark_inland_reg
subroutine mark_inland_reg(cs_res)
Create inland mask for regional grid.
Definition inland.F90:105
write_inland_reg
subroutine write_inland_reg(cs_res)
Write inland back to the orography data files for regional grid.
Definition inland.F90:393
nc_opchk
subroutine nc_opchk(stat, opname)
Check NetCDF return code and print error message.
Definition inland.F90:468
inland_mask
program inland_mask
This program creates the inland mask and writes it to the orography data files.
Definition inland.F90:9
mark_global_inland_rec_d
recursive subroutine mark_global_inland_rec_d(i, j, t, rd)
Recursively walk through neighbors marking inland points for global grid.
Definition inland.F90:162
mark_regional_inland_rec_d
recursive subroutine mark_regional_inland_rec_d(i, j, t, rd)
Recursively walk through neighbors marking inland points for regional grid.
Definition inland.F90:199
read_orog_reg
subroutine read_orog_reg(cs_res)
Read in orography (land fraction) data for regional grid.
Definition inland.F90:279
write_inland
subroutine write_inland(cs_res)
Write inland back to the orography data files for global grid.
Definition inland.F90:331
free_mem
subroutine free_mem()
Deallocate module arrays.
Definition inland.F90:456
cs_nb::nb_gp_idx
Neighboring cell descriptor.
Definition nb.F90:17