orog_mask_tools 1.14.0
Loading...
Searching...
No Matches
lakefrac.F90
Go to the documentation of this file.
1
4
19!#define DIAG_N_VERBOSE
20#define ADD_ATT_FOR_NEW_VAR
21PROGRAM lake_frac
22 USE netcdf
23 IMPLICIT NONE
24
25 CHARACTER(len=256) :: sfcdata_path
26 INTEGER :: cs_res, ncsmp, ncscp, i
27 INTEGER :: res_x, res_y
28
29 INTEGER*1, ALLOCATABLE :: lakestatus(:)
30 INTEGER*2, ALLOCATABLE :: lakedepth(:)
31 REAL, ALLOCATABLE :: cs_grid(:,:)
32 REAL, ALLOCATABLE :: cs_lakestatus(:), cs_lakedepth(:)
33 REAL, ALLOCATABLE :: src_grid_lon(:), src_grid_lat(:)
34
35 INTEGER :: tile_req, tile_beg, tile_end, binary_lake
36 REAL :: lake_cutoff
37
38 INTEGER, PARAMETER :: nlat = 21600, nlon = 43200
39 REAL, PARAMETER :: d2r = acos(-1.0) / 180.0
40 REAL, PARAMETER :: r2d = 180.0 /acos(-1.0)
41 REAL, PARAMETER :: pi = acos(-1.0)
42 real*8, PARAMETER :: gppdeg = 120.0
43 real*8, PARAMETER :: delta = 1.0 / 120.0
44
45 CHARACTER(len=32) :: arg, lakestatus_srce, lakedepth_srce
46 CHARACTER(len=256) :: lakedata_path
47 INTEGER :: stat
48
49 CALL getarg(0, arg) ! get the program name
50 IF (iargc() /= 5 .AND. iargc() /= 7) THEN
51 print*, 'Usage: ', trim(arg), &
52 ' [tile_num (0:all tiles, 7:regional)] [resolution (48,96, ...)] &
53 [lake data path] [lake status source] [lake depth source]'
54 print*, 'Or: ', trim(arg), &
55 ' [tile_num (0:all tiles, 7:regional)] [resolution (48,96, ...)] &
56 [lake data path] [lake status source] [lake depth source] [lake_cutoff] [binary_lake]'
57 stop -1
58 ENDIF
59 CALL getarg(1, arg)
60 READ(arg,*,iostat=stat) tile_req
61 CALL getarg(2, arg)
62 READ(arg,*,iostat=stat) cs_res
63 CALL getarg(3, lakedata_path)
64 CALL getarg(4, lakestatus_srce)
65 CALL getarg(5, lakedepth_srce)
66
67 IF (iargc() == 5) THEN
68 lake_cutoff = 0.50
69 binary_lake = 1
70 ELSE
71 CALL getarg(6, arg)
72 READ(arg,*,iostat=stat) lake_cutoff
73 CALL getarg(7, arg)
74 READ(arg,*,iostat=stat) binary_lake
75 ENDIF
76
77 print*, 'lake status source:', trim(lakestatus_srce)
78 print*, 'lake depth source:', trim(lakedepth_srce)
79 print*, 'lake cutoff:', lake_cutoff
80 print*, 'binary lake:', binary_lake
81
82 IF (tile_req == 0) THEN
83 tile_beg = 1; tile_end = 6
84 print*, 'Process tile 1 - 6 at resolution C',cs_res
85 ELSE IF (tile_req /= 7) THEN
86 tile_beg = tile_req; tile_end = tile_req
87 print*, 'Process tile',tile_req, 'at resolution C',cs_res
88 ELSE
89 tile_beg = 1; tile_end = 1
90 print*, 'Process regional tile (tile', tile_req, ') at resolution C',cs_res
91 ENDIF
92
93 ! read in grid spec data for each tile and concatenate them together
94
95 ncsmp = (2*cs_res+1)*(2*cs_res+1)*6
96 IF (tile_req /= 7) THEN
97 print*, 'Read in cubed sphere grid information ... ',ncsmp,'pairs of lat/lons'
98 ENDIF
99
100 IF (tile_req /= 7) THEN
101 ALLOCATE(cs_grid(ncsmp, 2))
102 CALL read_cubed_sphere_grid(cs_res, cs_grid)
103 ELSE
104 CALL read_cubed_sphere_reg_grid(cs_res, cs_grid, 3, res_x, res_y)
105 ENDIF
106
107 ! allocate and compute source grid
108 ALLOCATE(src_grid_lon(nlon), src_grid_lat(nlat))
109
110 IF (lakestatus_srce == "GLDBV2" .OR. lakestatus_srce == "GLDBV3") THEN
111 ! GLDB data points are at the lower right corners of the grid cells
112 DO i = 1, nlon
113 src_grid_lon(i) = -180.0 + delta * i
114 ENDDO
115 DO i = 1, nlat
116 src_grid_lat(i) = 90.0 - delta * i
117 ENDDO
118 ENDIF
119
120 IF (lakestatus_srce == "MODISP" .OR. lakestatus_srce == "VIIRS") THEN
121 ! GLDB data points are at the upprt left corners of the grid cells
122 DO i = 1, nlon
123 src_grid_lon(i) = -180.0 + delta * (i-1)
124 ENDDO
125 DO i = 1, nlat
126 src_grid_lat(i) = 90.0 - delta * (i-1)
127 ENDDO
128 ENDIF
129
130 ! read in lake data file
131 lakedata_path = trim(lakedata_path) // "/"
132 ALLOCATE(lakestatus(nlon*nlat),lakedepth(nlon*nlat))
133 print*, 'Read in lake data file ...'
134 CALL read_lakedata(lakedata_path,lakestatus,lakedepth,nlat,nlon)
135
136 ! calculate fraction numbers for all cs-cells
137 ncscp = cs_res*cs_res*6
138 ALLOCATE(cs_lakestatus(ncscp))
139 ALLOCATE(cs_lakedepth(ncscp))
140
141 print*, 'Calculate lake fraction and average depth for cubed-sphere cells ...'
142 CALL cal_lake_frac_depth(lakestatus,cs_lakestatus,lakedepth,cs_lakedepth)
143
144 ! write lake status (in terms of fraction) and lake depth(average, in meters)
145 ! to a netcdf file
146 IF (tile_req /= 7) THEN
147 print*, 'Write lake fraction/depth on cubed sphere grid cells to netCDF files ...'
148 CALL write_lakedata_to_orodata(cs_res, cs_lakestatus, cs_lakedepth)
149 ELSE
150 print*, 'Write lake fraction/depth on regional FV3 grid cells to a netCDF file ...'
151 CALL write_reg_lakedata_to_orodata(cs_res, res_x, res_y, cs_lakestatus, cs_lakedepth)
152 ENDIF
153
154 DEALLOCATE(cs_lakestatus,cs_lakedepth)
155 DEALLOCATE(cs_grid)
156 DEALLOCATE(lakestatus,lakedepth)
157 DEALLOCATE(src_grid_lat, src_grid_lon)
158
159 stop 0
160CONTAINS
161
170SUBROUTINE cal_lake_frac_depth(lakestat,cs_lakestat,lakedpth,cs_lakedpth)
171 INTEGER*1, INTENT(IN) :: lakestat(:)
172 INTEGER*2, INTENT(IN) :: lakedpth(:)
173 REAL, INTENT(OUT) :: cs_lakestat(:), cs_lakedpth(:)
174
175 real*8 lolf(2), lort(2), uplf(2), uprt(2), sd_ltmn(4), sd_ltmx(4)
176 real*8 :: v(2,4), p(2)
177 REAL :: latmin1, latmax1
178 REAL :: latmin, latmax, lonmin, lonmax, lontmp, lat_sz_max, lon_sz_max
179 INTEGER :: tile_num, i, j, gp, row, col, cs_grid_idx, cs_data_idx
180 INTEGER :: sidex_res, sidey_res, sidex_sz, sidey_sz
181 INTEGER :: stride_lat, stride_lon
182 INTEGER :: src_grid_lat_beg,src_grid_lat_end,src_grid_lon_beg,src_grid_lon_end
183 INTEGER :: src_grid_lon_beg1,src_grid_lon_end1,src_grid_lon_beg2,src_grid_lon_end2
184 INTEGER :: grid_ct, lake_ct, co_gc, tmp
185
186 INTEGER*1 :: lkst
187 INTEGER*2 :: lkdp
188 real*8 :: lake_dpth_sum, lake_avg_frac
189 LOGICAL :: two_section, enclosure_cnvx
190
191 IF (tile_req /= 7) THEN
192 sidex_res = cs_res; sidey_res = cs_res
193 ELSE
194 sidex_res = res_x; sidey_res = res_y
195 ENDIF
196
197 sidex_sz = 2*sidex_res+1; sidey_sz = 2*sidey_res+1
198
199 stride_lat = 1
200
201 lat_sz_max = 0.0
202 lon_sz_max = 0.0
203
204 cs_lakestat = 0.0
205
206 DO tile_num = tile_beg, tile_end
207 row = 2 + sidex_sz*(tile_num-1); col = 2
208 DO gp = 1, sidex_res*sidey_res
209 two_section = .false.
210 cs_grid_idx = (row-1)*sidex_sz+col
211 cs_data_idx = (tile_num-1)*sidex_res*sidey_res+gp
212 IF (abs(cs_grid(cs_grid_idx,1)) > 80.0 ) THEN !ignore lakes in very high latitude
213 cs_lakestat(cs_data_idx) = 0.0
214 cs_lakedpth(cs_data_idx) = 0.0
215 ! move to next cs cell
216 col = col + 2
217 IF (col > sidex_sz) THEN
218 col = 2
219 row = row + 2
220 ENDIF
221 cycle
222 ENDIF
223 ! get the four corners of the cs cell
224 lolf(1) = cs_grid(cs_grid_idx-sidex_sz-1, 1)
225 lolf(2) = cs_grid(cs_grid_idx-sidex_sz-1, 2)
226 IF (lolf(2) > 180.0) lolf(2) = lolf(2) - 360.0
227 lort(1) = cs_grid(cs_grid_idx-sidex_sz+1, 1)
228 lort(2) = cs_grid(cs_grid_idx-sidex_sz+1, 2)
229 IF (lort(2) > 180.0) lort(2) = lort(2) - 360.0
230 uplf(1) = cs_grid(cs_grid_idx+sidex_sz-1,1)
231 uplf(2) = cs_grid(cs_grid_idx+sidex_sz-1,2)
232 IF (uplf(2) > 180.0) uplf(2) = uplf(2) - 360.0
233 uprt(1) = cs_grid(cs_grid_idx+sidex_sz+1,1)
234 uprt(2) = cs_grid(cs_grid_idx+sidex_sz+1,2)
235
236 v(1,1) = lolf(1); v(2,1) = lolf(2)
237 v(1,2) = lort(1); v(2,2) = lort(2)
238 v(1,3) = uprt(1); v(2,3) = uprt(2)
239 v(1,4) = uplf(1); v(2,4) = uplf(2)
240 v(:,:) = v(:,:) * d2r
241
242 IF (uprt(2) > 180.0) uprt(2) = uprt(2) - 360.0
243 ! gather the candidate indices in lakestat
244#ifdef LIMIT_CAL
245 CALL find_limit (lolf, lort, sd_ltmn(1), sd_ltmx(1))
246 CALL find_limit (lort, uprt, sd_ltmn(2), sd_ltmx(2))
247 CALL find_limit (uprt, uplf, sd_ltmn(3), sd_ltmx(3))
248 CALL find_limit (uplf, lolf, sd_ltmn(4), sd_ltmx(4))
249 latmin = min(sd_ltmn(1),min(sd_ltmn(2),min(sd_ltmn(3),sd_ltmn(4))))
250 latmax = max(sd_ltmx(1),max(sd_ltmx(2),max(sd_ltmx(3),sd_ltmx(4))))
251#endif
252 latmin = min(lolf(1),min(lort(1),min(uplf(1),uprt(1))))
253 latmax = max(lolf(1),max(lort(1),max(uplf(1),uprt(1))))
254 lonmin = min(lolf(2),min(lort(2),min(uplf(2),uprt(2))))
255 lonmax = max(lolf(2),max(lort(2),max(uplf(2),uprt(2))))
256! lat_sz_max = max(lat_sz_max, (latmax-latmin))
257! lon_sz_max = max(lon_sz_max, (lonmax-lonmin))
258
259 src_grid_lat_beg = nint((90.0-latmax)*gppdeg+0.5)
260 src_grid_lat_end = nint((90.0-latmin)*gppdeg+0.5)
261 src_grid_lon_beg = nint((180.0+lonmin)*gppdeg+0.5)
262 src_grid_lon_end = nint((180.0+lonmax)*gppdeg+0.5)
263
264 IF (src_grid_lat_beg > src_grid_lat_end) THEN
265 tmp = src_grid_lat_beg
266 src_grid_lat_beg = src_grid_lat_end
267 src_grid_lat_end = tmp
268 ENDIF
269 IF (src_grid_lon_beg > src_grid_lon_end) THEN
270 tmp = src_grid_lon_beg
271 src_grid_lon_beg = src_grid_lon_end
272 src_grid_lon_end = tmp
273 ENDIF
274 IF ((src_grid_lon_end - src_grid_lon_beg) > nlon*0.75) THEN
275 two_section = .true.
276 src_grid_lon_beg1 = src_grid_lon_end
277 src_grid_lon_end1 = nlon
278 src_grid_lon_beg2 = 1
279 src_grid_lon_end2 = src_grid_lon_beg
280 ENDIF
281
282#ifdef DIAG_N_VERBOSE
283 print*, 'cell centre lat/lon =', &
284 gp, cs_res*cs_res, cs_grid(cs_grid_idx,1),cs_grid(cs_grid_idx,2)
285 print*, 'lat index range and stride', &
286 src_grid_lat_beg,src_grid_lat_end,stride_lat
287 print*, 'lat range ', &
288 src_grid_lat(src_grid_lat_beg),src_grid_lat(src_grid_lat_end)
289#endif
290 lake_ct = 0; grid_ct = 0
291 lake_dpth_sum = 0.0
292 lake_avg_frac = 0.0
293 DO j = src_grid_lat_beg, src_grid_lat_end, stride_lat
294 stride_lon = int(1.0/cos(src_grid_lat(j)*d2r)*real(stride_lat))
295#ifdef DIAG_N_VERBOSE
296 IF (j == src_grid_lat_beg) THEN
297 print*, 'lon index range and stride', &
298 src_grid_lon_beg,src_grid_lon_end,stride_lon
299 print*, 'lon range ', &
300 src_grid_lon(src_grid_lon_beg),src_grid_lon(src_grid_lon_end)
301 IF (two_section .eqv. .true.) THEN
302 print*, 'section1 index lon range and stride', &
303 src_grid_lon_beg1,src_grid_lon_end1,stride_lon
304 print*, 'section1 lon range ', &
305 src_grid_lon(src_grid_lon_beg1),src_grid_lon(src_grid_lon_end1)
306 print*, 'section2 index lon range and stride', &
307 src_grid_lon_beg2,src_grid_lon_end2,stride_lon
308 print*, 'section2 lon range ', &
309 src_grid_lon(src_grid_lon_beg2),src_grid_lon(src_grid_lon_end2)
310 ENDIF
311 ENDIF
312#endif
313 IF (two_section .eqv. .false.) THEN
314 DO i = src_grid_lon_beg, src_grid_lon_end, stride_lon
315 p(1) = src_grid_lat(j); p(2) = src_grid_lon(i)
316 p(:) = p(:)*d2r
317 IF(enclosure_cnvx(v, 4, p, co_gc) .eqv. .true.) THEN
318 grid_ct = grid_ct+1
319 lkst = lakestat((j-1)*nlon+i); lkdp = lakedpth((j-1)*nlon+i)
320 CALL lake_cell_comp(lkst, lkdp, lake_ct, lake_avg_frac, lake_dpth_sum)
321 ENDIF
322 ENDDO
323 ELSE
324 DO i = src_grid_lon_beg1, src_grid_lon_end1, stride_lon
325 p(1) = src_grid_lat(j); p(2) = src_grid_lon(i)
326 p(:) = p(:)*d2r
327 IF(enclosure_cnvx(v, 4, p, co_gc) .eqv. .true.) THEN
328 grid_ct = grid_ct+1
329 lkst = lakestat((j-1)*nlon+i); lkdp = lakedpth((j-1)*nlon+i)
330 CALL lake_cell_comp(lkst, lkdp, lake_ct, lake_avg_frac, lake_dpth_sum)
331 ENDIF
332 ENDDO
333 DO i = src_grid_lon_beg2, src_grid_lon_end2, stride_lon
334 p(1) = src_grid_lat(j); p(2) = src_grid_lon(i)
335 p(:) = p(:)*d2r
336 IF(enclosure_cnvx(v, 4, p, co_gc) .eqv. .true.) THEN
337 grid_ct = grid_ct+1
338 lkst = lakestat((j-1)*nlon+i); lkdp = lakedpth((j-1)*nlon+i)
339 CALL lake_cell_comp(lkst, lkdp, lake_ct, lake_avg_frac, lake_dpth_sum)
340 ENDIF
341 ENDDO
342 ENDIF
343 ENDDO
344 IF (lakestatus_srce == "GLDBV3" .OR. lakestatus_srce == "GLDBV2" .OR. &
345 lakestatus_srce == "VIIRS" ) THEN
346 cs_lakestat(cs_data_idx)=real(lake_ct)/real(grid_ct)
347 ENDIF
348 IF (lakestatus_srce == "MODISP") THEN
349 cs_lakestat(cs_data_idx)=lake_avg_frac/real(grid_ct)
350 ENDIF
351 IF (lake_ct /= 0) THEN
352 cs_lakedpth(cs_data_idx)=lake_dpth_sum/real(lake_ct)/10.0 !convert to meter
353 ELSE
354 cs_lakedpth(cs_data_idx)=0.0
355 ENDIF
356#ifdef DIAG_N_VERBOSE
357 print*, 'tile_num, row, col:', tile_num, row, col
358 print*, 'grid_ct, lake_ct = ', grid_ct, lake_ct
359 print*, 'lake_frac= ', cs_lakestat(cs_data_idx)
360 print*, 'lake_depth (avg) = ', cs_lakedpth(cs_data_idx)
361#endif
362
363 ! move to the next control volume
364 col = col + 2
365 IF (col > sidex_sz) THEN
366 col = 2
367 row = row + 2
368 ENDIF
369 ENDDO
370 print "('*'$)" ! progress '*'
371 ENDDO
372 print*, ''
373
374END SUBROUTINE cal_lake_frac_depth
375
384SUBROUTINE lake_cell_comp(lkst, lkdp, lake_ct, lake_avg_frac, lake_dpth_sum)
385 INTEGER*1, INTENT(IN) :: lkst
386 INTEGER*2, INTENT(IN) :: lkdp
387 INTEGER, INTENT(OUT) :: lake_ct
388 real*8, INTENT(OUT) :: lake_avg_frac, lake_dpth_sum
389
390 IF (lkst /= 0) THEN ! lake point
391 lake_ct = lake_ct+1
392 IF (lakestatus_srce == "GLDBV3" .OR. lakestatus_srce == "GLDBV2") THEN
393 IF (lkdp <= 0) THEN
394 IF (lkst == 4) THEN
395 lake_dpth_sum = lake_dpth_sum+30.0
396 ELSE
397 lake_dpth_sum = lake_dpth_sum+100.0
398 ENDIF
399 ELSE
400 lake_dpth_sum = lake_dpth_sum+real(lkdp)
401 ENDIF
402 ENDIF
403 IF (lakestatus_srce == "MODISP" .OR. lakestatus_srce == "VIIRS") THEN
404 lake_avg_frac = lake_avg_frac + real(lkst) / 100.0
405 IF (lkdp <= 0) THEN
406 lake_dpth_sum = lake_dpth_sum+100.0
407 ELSE
408 lake_dpth_sum = lake_dpth_sum+real(lkdp)
409 ENDIF
410 ENDIF
411 ENDIF
412END SUBROUTINE lake_cell_comp
413
423SUBROUTINE read_cubed_sphere_grid(res, grid)
424 INTEGER, INTENT(IN) :: res
425 REAL, INTENT(OUT) :: grid(:,:)
426
427 real*8, ALLOCATABLE :: lat(:), lon(:)
428 INTEGER :: side_sz, tile_sz, ncid, varid
429 INTEGER :: i, tile_beg, tile_end, stat
430 CHARACTER(len=256) :: gridfile_path,gridfile
431 CHARACTER(len=1) ich
432 CHARACTER(len=4) res_ch
433
434 side_sz = 2*res+1
435 tile_sz = side_sz*side_sz
436 ALLOCATE(lat(tile_sz), lon(tile_sz))
437
438 IF (tile_req == 0) THEN
439 tile_beg = 1; tile_end = 6
440 ELSE
441 tile_beg = tile_req; tile_end = tile_req
442 ENDIF
443 WRITE(res_ch,'(I4)') res
444 gridfile_path = "./"
445 DO i = tile_beg, tile_end
446 WRITE(ich, '(I1)') i
447 gridfile = trim(gridfile_path)//"C"//trim(adjustl(res_ch))//"_grid.tile"//ich//".nc"
448
449 print*, 'Open cubed sphere grid spec file ', trim(gridfile)
450
451 stat = nf90_open(trim(gridfile), nf90_nowrite, ncid)
452 CALL nc_opchk(stat,'nf90_open')
453
454 stat = nf90_inq_varid(ncid,'x',varid)
455 CALL nc_opchk(stat,'nf90_inq_lon')
456 stat = nf90_get_var(ncid,varid,lon,start=(/1,1/),count=(/side_sz,side_sz/))
457 CALL nc_opchk(stat,'nf90_get_var_lon')
458 stat = nf90_inq_varid(ncid,'y',varid)
459 CALL nc_opchk(stat,'nf90_inq_lat')
460 stat = nf90_get_var(ncid,varid,lat,start=(/1,1/),count=(/side_sz,side_sz/))
461 CALL nc_opchk(stat,'nf90_get_var_lat')
462 stat = nf90_close(ncid)
463 CALL nc_opchk(stat,'nf90_close')
464
465 tile_beg = (i-1)*tile_sz+1
466 tile_end = i*tile_sz
467 grid(tile_beg:tile_end,1) = lat(1:tile_sz)
468 grid(tile_beg:tile_end,2) = lon(1:tile_sz)
469 END DO
470 DEALLOCATE(lat,lon)
471
472END SUBROUTINE read_cubed_sphere_grid
473
483SUBROUTINE read_cubed_sphere_reg_grid(res, grid, halo_depth, res_x, res_y)
484 INTEGER, INTENT(IN) :: res, halo_depth
485 INTEGER, INTENT(OUT) :: res_x, res_y
486 REAL, ALLOCATABLE, INTENT(OUT) :: grid(:,:)
487
488 real*8, ALLOCATABLE :: lat(:), lon(:)
489 INTEGER :: sidex_sz, sidey_sz, tile_sz, ncid, varid, dimid
490 INTEGER :: x_start, y_start
491 INTEGER :: nxp, nyp, stat
492 CHARACTER(len=256) :: gridfile_path,gridfile
493 CHARACTER(len=1) ich
494 CHARACTER(len=4) res_ch
495 CHARACTER(len=8) dimname
496
497 WRITE(res_ch,'(I4)') res
498 gridfile_path = './'
499 gridfile = trim(gridfile_path)//"C"//trim(adjustl(res_ch))//"_grid.tile7.nc"
500
501 print*, 'Open cubed sphere grid spec file ', trim(gridfile)
502
503 stat = nf90_open(trim(gridfile), nf90_nowrite, ncid)
504 CALL nc_opchk(stat,'nf90_open')
505
506 stat = nf90_inq_dimid(ncid,'nxp',dimid)
507 CALL nc_opchk(stat,'nf90_inq_dimid: nxp')
508 stat = nf90_inquire_dimension(ncid,dimid,dimname,len=nxp)
509 CALL nc_opchk(stat,'nf90_inquire_dimension: nxp')
510
511 stat = nf90_inq_dimid(ncid,'nyp',dimid)
512 CALL nc_opchk(stat,'nf90_inq_dimid: nyp')
513 stat = nf90_inquire_dimension(ncid,dimid,dimname,len=nyp)
514 CALL nc_opchk(stat,'nf90_inquire_dimension: nyp')
515
516 sidex_sz = nxp
517 sidey_sz = nyp
518 tile_sz = sidex_sz*sidey_sz
519 ALLOCATE(lat(tile_sz), lon(tile_sz))
520
521! x_start = halo_depth+1; y_start = halo_depth+1
522 x_start = 1; y_start = 1
523 stat = nf90_inq_varid(ncid,'x',varid)
524 CALL nc_opchk(stat,'nf90_inq_lon')
525 stat = nf90_get_var(ncid,varid,lon,start=(/x_start,y_start/),count=(/sidex_sz,sidey_sz/))
526 CALL nc_opchk(stat,'nf90_get_var_lon')
527 stat = nf90_inq_varid(ncid,'y',varid)
528 CALL nc_opchk(stat,'nf90_inq_lat')
529 stat = nf90_get_var(ncid,varid,lat,start=(/x_start,y_start/),count=(/sidex_sz,sidey_sz/))
530 CALL nc_opchk(stat,'nf90_get_var_lat')
531 stat = nf90_close(ncid)
532 CALL nc_opchk(stat,'nf90_close')
533
534 ALLOCATE(grid(tile_sz,2))
535 grid(1:tile_sz,1) = lat(1:tile_sz)
536 grid(1:tile_sz,2) = lon(1:tile_sz)
537
538 res_x = sidex_sz/2; res_y = sidey_sz/2
539 DEALLOCATE(lat,lon)
540
541END SUBROUTINE read_cubed_sphere_reg_grid
542
553SUBROUTINE read_lakedata(lakedata_path,lake_stat,lake_dpth,nlat,nlon)
554 CHARACTER(len=256), INTENT(IN) :: lakedata_path
555 INTEGER*1, INTENT(OUT) :: lake_stat(:)
556 INTEGER*2, INTENT(OUT) :: lake_dpth(:)
557 INTEGER, INTENT(IN) :: nlat, nlon
558
559 CHARACTER(len=256) lakefile
560 INTEGER :: data_sz, i
561
562 data_sz = nlon*nlat
563
564! read in 30 arc seconds lake data
565 ! lake fraction data
566 lakefile = trim(lakedata_path) // "GlobalLakeStatus_GLDBv3release.dat"
567 IF (lakestatus_srce == "GLDBV2") THEN
568 lakefile = trim(lakedata_path) // "GlobalLakeStatus.dat"
569 ENDIF
570 IF (lakestatus_srce == "GLDBV3") THEN
571 lakefile = trim(lakedata_path) // "GlobalLakeStatus_GLDBv3release.dat"
572 ENDIF
573 IF (lakestatus_srce == "MODISP") THEN
574! lakefile = trim(lakedata_path) // "GlobalLakeStatus_MODIS15s.dat"
575 lakefile = trim(lakedata_path) // "GlobalLakeStatus_MODISp.dat"
576 ENDIF
577 IF (lakestatus_srce == "VIIRS") THEN
578 lakefile = trim(lakedata_path) // "GlobalLakeStatus_VIIRS.dat"
579 ENDIF
580 OPEN(10, file=lakefile, form='unformatted', access='direct', status='old', recl=data_sz*1)
581 READ(10,rec=1) lake_stat
582 CLOSE(10)
583
584 ! lake depth data
585 lakefile = trim(lakedata_path) // "GlobalLakeDepth_GLDBv3release.dat"
586 IF (lakedepth_srce == "GLDBV2") THEN
587 lakefile = trim(lakedata_path) // "GlobalLakeDepth.dat"
588 ENDIF
589 IF (lakedepth_srce == "GLDBV3") THEN
590 lakefile = trim(lakedata_path) // "GlobalLakeDepth_GLDBv3release.dat"
591 ENDIF
592 IF (lakedepth_srce == "GLOBATHY") THEN
593 lakefile = trim(lakedata_path) // "GlobalLakeDepth_GLOBathy.dat"
594 ENDIF
595 OPEN(10, file=lakefile, form='unformatted', access='direct', status='old', recl=data_sz*2)
596 READ(10,rec=1) lake_dpth
597 CLOSE(10)
598
599END SUBROUTINE read_lakedata
600
609SUBROUTINE write_lakedata_to_orodata(cs_res, cs_lakestat, cs_lakedpth)
610 USE netcdf
611 INTEGER, INTENT(IN) :: cs_res
612 REAL, INTENT(IN) :: cs_lakestat(:)
613 REAL, INTENT(IN) :: cs_lakedpth(:)
614
615 INTEGER :: tile_sz, tile_num
616 INTEGER :: stat, ncid, x_dimid, y_dimid, varid, dimids(2)
617 INTEGER :: lake_frac_id, lake_depth_id
618 INTEGER :: land_frac_id, slmsk_id, inland_id, geolon_id, geolat_id
619 CHARACTER(len=256) :: filename,string,lakeinfo
620 CHARACTER(len=1) :: ich
621 CHARACTER(len=4) res_ch
622 REAL :: lake_frac(cs_res*cs_res),lake_depth(cs_res*cs_res)
623 REAL :: geolon(cs_res*cs_res),geolat(cs_res*cs_res)
624 REAL :: land_frac(cs_res*cs_res),slmsk(cs_res*cs_res),inland(cs_res*cs_res)
625 real, parameter :: epsil=1.e-6 ! numerical min for lake_frac/land_frac
626 real :: land_cutoff=1.e-4 ! land_frac=0 if it is < land_cutoff
627
628 INTEGER :: i, j
629
630 tile_sz = cs_res*cs_res
631
632 WRITE(res_ch,'(I4)') cs_res
633 DO tile_num = tile_beg, tile_end
634 WRITE(ich, '(I1)') tile_num
635! filename = "C" // trim(adjustl(res_ch)) // "_oro_data.tile" // ich // ".nc"
636! filename = "oro_data.tile" // ich // ".nc"
637 filename = "oro.C" // trim(adjustl(res_ch)) // ".tile" // ich // ".nc"
638 print *,'Read, update, and write ',trim(filename)
639 stat = nf90_open(filename, nf90_write, ncid)
640 CALL nc_opchk(stat, "nf90_open oro_data.nc")
641 stat = nf90_inq_dimid(ncid, "lon", x_dimid)
642 CALL nc_opchk(stat, "nf90_inq_dim: x")
643 stat = nf90_inq_dimid(ncid, "lat", y_dimid)
644 CALL nc_opchk(stat, "nf90_inq_dim: y")
645! dimids = (/ y_dimid, x_dimid /)
646! original orodata netcdf file uses (y, x) order, so we made change to match it.
647 dimids = (/ x_dimid, y_dimid /)
648 stat = nf90_inq_varid(ncid, "land_frac", land_frac_id)
649 CALL nc_opchk(stat, "nf90_inq_varid: land_frac")
650 stat = nf90_inq_varid(ncid, "slmsk", slmsk_id)
651 CALL nc_opchk(stat, "nf90_inq_varid: slmsk")
652! define 2 new variables
653 stat = nf90_redef(ncid)
654 CALL nc_opchk(stat, "nf90_redef")
655 stat = nf90_def_var(ncid,"lake_frac",nf90_float,dimids,lake_frac_id)
656 CALL nc_opchk(stat, "nf90_def_var: lake_frac")
657#ifdef ADD_ATT_FOR_NEW_VAR
658 stat = nf90_put_att(ncid, lake_frac_id,'coordinates','geolon geolat')
659 CALL nc_opchk(stat, "nf90_put_att: lake_frac:coordinates")
660 stat = nf90_put_att(ncid, lake_frac_id,'long_name','lake fraction')
661 CALL nc_opchk(stat, "nf90_put_att: lake_frac:long_name")
662 stat = nf90_put_att(ncid, lake_frac_id,'unit','fraction')
663 CALL nc_opchk(stat, "nf90_put_att: lake_frac:unit")
664 write(lakeinfo,'(a,f4.2,a,i1)') ' lake_frac cutoff=',lake_cutoff,'; binary_lake=',binary_lake
665 IF (lakestatus_srce == "GLDBV3") THEN
666 write(string,'(2a)') 'based on GLDBv3 (Choulga et al. 2019); missing lakes & added based on land_frac in this dataset;', &
667 trim(lakeinfo)
668 ELSE IF (lakestatus_srce == "GLDBV2") THEN
669 write(string,'(2a)') 'based on GLDBv2 (Choulga et al. 2014); missing lakes & added based on land_frac in this dataset;', &
670 trim(lakeinfo)
671 ELSE IF (lakestatus_srce == "MODISP") THEN
672 write(string,'(2a)') 'based on MODIS (2011-2015) product updated with two &
673 Landsat products: the JRC water product (2016-2020) and the GLC-FCS30 (2020); &
674 the source data set was created by Chengquan Huang of UMD;',trim(lakeinfo)
675 ELSE IF (lakestatus_srce == "VIIRS") THEN
676 write(string,'(2a)') 'based on multi-year VIIRS global surface type &
677 classification map (2012-2019); the source data set was created by &
678 Chengquan Huang of UMD and Michael Barlage of NOAA;',trim(lakeinfo)
679 ENDIF
680 stat = nf90_put_att(ncid, lake_frac_id,'description',trim(string))
681 CALL nc_opchk(stat, "nf90_put_att: lake_frac:description")
682#endif
683 stat = nf90_def_var(ncid,"lake_depth",nf90_float,dimids,lake_depth_id)
684 CALL nc_opchk(stat, "nf90_def_var: lake_depth")
685#ifdef ADD_ATT_FOR_NEW_VAR
686 stat = nf90_put_att(ncid, lake_depth_id,'coordinates','geolon geolat')
687 CALL nc_opchk(stat, "nf90_put_att: lake_depth:coordinates")
688 stat = nf90_put_att(ncid, lake_depth_id,'long_name','lake depth')
689 CALL nc_opchk(stat, "nf90_put_att: lake_depth:long_name")
690 stat = nf90_put_att(ncid, lake_depth_id,'unit','meter')
691 CALL nc_opchk(stat, "nf90_put_att: lake_depth:long_name")
692 IF (lakedepth_srce == "GLDBV3") THEN
693 stat = nf90_put_att(ncid, lake_depth_id,'description', &
694 'based on GLDBv3 (Choulga et al. 2019); missing depth set to 10m &
695 (except to 211m in Caspian Sea); spurious large pos. depths are left unchanged.')
696 ELSE IF (lakedepth_srce == "GLDBV2") THEN
697 stat = nf90_put_att(ncid, lake_depth_id,'description', &
698 'based on GLDBv2 (Choulga et al. 2014); missing depth set to 10m &
699 (except to 211m in Caspian Sea); spurious large pos. depths are left unchanged.')
700 ELSE IF (lakedepth_srce == "GLOBATHY") THEN
701 stat = nf90_put_att(ncid, lake_depth_id,'description', &
702 'based on GLOBathy data resampled and projected to the MODIS domain.')
703 ENDIF
704 CALL nc_opchk(stat, "nf90_put_att: lake_depth:description")
705#endif
706
707 write(string,'(a,es8.1)') 'land_frac and lake_frac are adjusted such that &
708 their sum is 1 at points where inland=1; land_frac cutoff is',land_cutoff
709 stat = nf90_put_att(ncid, land_frac_id,'description',trim(string))
710 CALL nc_opchk(stat, "nf90_put_att: land_frac:description")
711
712 write(string,'(a)') 'slmsk = nint(land_frac)'
713 stat = nf90_put_att(ncid, slmsk_id,'description',trim(string))
714 CALL nc_opchk(stat, "nf90_put_att: slmsk:description")
715
716 stat = nf90_enddef(ncid)
717 CALL nc_opchk(stat, "nf90_enddef")
718
719! read in geolon and geolat and 2 variables from orog data file
720 stat = nf90_inq_varid(ncid, "geolon", geolon_id)
721 CALL nc_opchk(stat, "nf90_inq_varid: geolon")
722 stat = nf90_inq_varid(ncid, "geolat", geolat_id)
723 CALL nc_opchk(stat, "nf90_inq_varid: geolat")
724 stat = nf90_inq_varid(ncid, "land_frac", land_frac_id)
725 CALL nc_opchk(stat, "nf90_inq_varid: land_frac")
726 stat = nf90_inq_varid(ncid, "slmsk", slmsk_id)
727 CALL nc_opchk(stat, "nf90_inq_varid: slmsk")
728 stat = nf90_inq_varid(ncid, "inland", inland_id)
729 CALL nc_opchk(stat, "nf90_inq_varid: inland")
730 stat = nf90_get_var(ncid, geolon_id, geolon, &
731 start = (/ 1, 1 /), count = (/ cs_res, cs_res /) )
732 CALL nc_opchk(stat, "nf90_get_var: geolon")
733 stat = nf90_get_var(ncid, geolat_id, geolat, &
734 start = (/ 1, 1 /), count = (/ cs_res, cs_res /) )
735 CALL nc_opchk(stat, "nf90_get_var: geolat")
736 stat = nf90_get_var(ncid, land_frac_id, land_frac, &
737 start = (/ 1, 1 /), count = (/ cs_res, cs_res /) )
738 CALL nc_opchk(stat, "nf90_get_var: land_frac")
739 stat = nf90_get_var(ncid, slmsk_id, slmsk, &
740 start = (/ 1, 1 /), count = (/ cs_res, cs_res /) )
741 CALL nc_opchk(stat, "nf90_get_var: slmsk")
742 stat = nf90_get_var(ncid, inland_id, inland, &
743 start = (/ 1, 1 /), count = (/ cs_res, cs_res /) )
744 CALL nc_opchk(stat, "nf90_get_var: inland")
745
746 lake_frac(:) = cs_lakestat((tile_num-1)*tile_sz+1:tile_num*tile_sz)
747 lake_depth(:) = cs_lakedepth((tile_num-1)*tile_sz+1:tile_num*tile_sz)
748
749! include Caspian Sea and Aral Sea if GLDB data set is used, and
750! exclude lakes in the coastal areas of Antarctica if MODIS data set is used
751 CALL include_exclude_lakes(lake_frac,land_frac,lake_depth,geolat,geolon,tile_num)
752
753! epsil is "numerical" nonzero min for lake_frac/land_frac
754! lake_cutoff/land_cutoff is practical min for lake_frac/land_frac
755 IF (min(lake_cutoff,land_cutoff) < epsil) then
756 print *,'lake_cutoff/land_cutoff cannot be smaller than epsil, reset...'
757 lake_cutoff=max(epsil,lake_cutoff)
758 land_cutoff=max(epsil,land_cutoff)
759 ENDIF
760
761! adjust land_frac and lake_frac, and make sure land_frac+lake_frac=1 at inland points
762 DO i = 1, tile_sz
763 if (land_frac(i)< land_cutoff) land_frac(i)=0.
764 if (land_frac(i)>1.-land_cutoff) land_frac(i)=1.
765
766 if (inland(i) /= 1.) then
767 lake_frac(i) = 0.
768 endif
769
770 if (lake_frac(i) < lake_cutoff) then
771 lake_frac(i)=0.
772 elseif (binary_lake == 1) then
773 lake_frac(i)=1.
774 end if
775 if (lake_frac(i) > 1.-epsil) lake_frac(i)=1.
776 ENDDO
777
778! finalize land_frac/slmsk based on modified lake_frac
779 DO i = 1, tile_sz
780 if (inland(i) == 1.) then
781 land_frac(i) = 1. - lake_frac(i)
782 end if
783
784 if (lake_frac(i) < lake_cutoff) then
785 lake_depth(i)=0.
786 elseif (lake_frac(i) >= lake_cutoff .and. lake_depth(i)==0.) then
787 lake_depth(i)=10.
788 end if
789 slmsk(i) = nint(land_frac(i))
790 ENDDO
791
792! write 2 new variables
793 stat = nf90_put_var(ncid, lake_frac_id, lake_frac, &
794 start = (/ 1, 1 /), count = (/ cs_res, cs_res /) )
795 CALL nc_opchk(stat, "nf90_put_var: lake_frac")
796
797 stat = nf90_put_var(ncid, lake_depth_id, lake_depth, &
798 start = (/ 1, 1 /), count = (/ cs_res, cs_res /) )
799 CALL nc_opchk(stat, "nf90_put_var: lake_depth")
800
801! write back 2 modified variables: land_frac and slmsk
802 stat = nf90_put_var(ncid, land_frac_id, land_frac, &
803 start = (/ 1, 1 /), count = (/ cs_res, cs_res /) )
804 CALL nc_opchk(stat, "nf90_put_var: land_frac")
805
806 stat = nf90_put_var(ncid, slmsk_id, slmsk, &
807 start = (/ 1, 1 /), count = (/ cs_res, cs_res /) )
808 CALL nc_opchk(stat, "nf90_put_var: slmsk")
809
810 stat = nf90_close(ncid)
811 CALL nc_opchk(stat, "nf90_close")
812 ENDDO
813
814END SUBROUTINE write_lakedata_to_orodata
815
826SUBROUTINE write_reg_lakedata_to_orodata(cs_res, tile_x_dim, tile_y_dim, cs_lakestat, cs_lakedpth)
827 USE netcdf
828 INTEGER, INTENT(IN) :: cs_res, tile_x_dim, tile_y_dim
829 REAL, INTENT(IN) :: cs_lakestat(:)
830 REAL, INTENT(IN) :: cs_lakedpth(:)
831
832 INTEGER :: tile_sz, tile_num
833 INTEGER :: stat, ncid, x_dimid, y_dimid, varid, dimids(2)
834 INTEGER :: lake_frac_id, lake_depth_id
835 INTEGER :: land_frac_id, slmsk_id, geolon_id, geolat_id, inland_id
836 CHARACTER(len=256) :: filename,string
837 CHARACTER(len=1) :: ich
838 CHARACTER(len=4) res_ch
839
840 REAL, ALLOCATABLE :: lake_frac(:), lake_depth(:), geolon(:), geolat(:)
841 REAL, ALLOCATABLE :: land_frac(:), slmsk(:), inland(:)
842
843 real, parameter :: epsil=1.e-6 ! numerical min for lake_frac/land_frac
844 real :: land_cutoff=1.e-6 ! land_frac=0 if it is < land_cutoff
845
846 INTEGER :: i, j, var_id
847
848! include "netcdf.inc"
849 tile_sz = tile_x_dim*tile_y_dim
850
851 ALLOCATE(lake_frac(tile_sz), lake_depth(tile_sz))
852 ALLOCATE(geolon(tile_sz), geolat(tile_sz))
853 ALLOCATE(land_frac(tile_sz), slmsk(tile_sz), inland(tile_sz))
854
855 WRITE(res_ch,'(I4)') cs_res
856 tile_num = 7
857 WRITE(ich, '(I1)') tile_num
858! filename = "C" // trim(adjustl(res_ch)) // "_oro_data.tile" // ich // ".halo4.nc"
859 filename = "oro.C" // trim(adjustl(res_ch)) // ".tile" // ich // ".nc"
860 print*, 'Open and write regional orography data netCDF file ', trim(filename)
861 stat = nf90_open(filename, nf90_write, ncid)
862 CALL nc_opchk(stat, "nf90_open oro_data.nc")
863 stat = nf90_inq_dimid(ncid, "lon", x_dimid)
864 CALL nc_opchk(stat, "nf90_inq_dim: x")
865 stat = nf90_inq_dimid(ncid, "lat", y_dimid)
866 CALL nc_opchk(stat, "nf90_inq_dim: y")
867 dimids = (/ x_dimid, y_dimid /)
868
869 stat = nf90_inq_varid(ncid, "land_frac", land_frac_id)
870 CALL nc_opchk(stat, "nf90_inq_varid: land_frac")
871 stat = nf90_inq_varid(ncid, "slmsk", slmsk_id)
872 CALL nc_opchk(stat, "nf90_inq_varid: slmsk")
873
874! define 2 new variables, lake_frac and lake_depth
875 stat = nf90_redef(ncid)
876 CALL nc_opchk(stat, "nf90_redef")
877
878 IF (nf90_inq_varid(ncid, "lake_frac",lake_frac_id) /= 0) THEN
879 stat = nf90_def_var(ncid,"lake_frac",nf90_float,dimids,lake_frac_id)
880 CALL nc_opchk(stat, "nf90_def_var: lake_frac")
881#ifdef ADD_ATT_FOR_NEW_VAR
882 stat = nf90_put_att(ncid, lake_frac_id,'coordinates','geolon geolat')
883 CALL nc_opchk(stat, "nf90_put_att: lake_frac:coordinates")
884 stat = nf90_put_att(ncid, lake_frac_id,'long_name','lake fraction')
885 CALL nc_opchk(stat, "nf90_put_att: lake_frac:long_name")
886 stat = nf90_put_att(ncid, lake_frac_id,'unit','fraction')
887 CALL nc_opchk(stat, "nf90_put_att: lake_frac:unit")
888 IF (lakestatus_srce == "GLDBV3") THEN
889 write(string,'(a,es8.1)') 'based on GLDBv3 (Choulga et al. 2019); missing lakes &
890 added based on land_frac in this dataset; lake_frac cutoff:',lake_cutoff
891 ELSE IF (lakestatus_srce == "GLDBV2") THEN
892 write(string,'(a,es8.1)') 'based on GLDBv2 (Choulga et al. 2014); missing lakes &
893 added based on land_frac in this dataset; lake_frac cutoff:',lake_cutoff
894 ELSE IF (lakestatus_srce == "MODISP") THEN
895 write(string,'(a,es8.1)') 'based on MODIS (2011-2015) product updated with two &
896 Landsat products: the JRC water product (2016-2020) and the GLC-FCS30 (2020); &
897 the source data set was created by Chengquan Huang of UMD; &
898 lake_frac cutoff:',lake_cutoff
899 ELSE IF (lakestatus_srce == "VIIRS") THEN
900 write(string,'(a,es8.1)') 'based on multi-year VIIRS global surface type &
901 classification map (2012-2019); the source data set was created by &
902 Chengquan Huang of UMD and Michael Barlage of NOAA; &
903 lake_frac cutoff:',lake_cutoff
904 ENDIF
905 stat = nf90_put_att(ncid, lake_frac_id,'description',trim(string))
906 CALL nc_opchk(stat, "nf90_put_att: lake_frac:description")
907#endif
908 ENDIF
909 IF (nf90_inq_varid(ncid, "lake_depth",lake_depth_id) /= 0) THEN
910 stat = nf90_def_var(ncid,"lake_depth",nf90_float,dimids,lake_depth_id)
911 CALL nc_opchk(stat, "nf90_def_var: lake_depth")
912#ifdef ADD_ATT_FOR_NEW_VAR
913 stat = nf90_put_att(ncid, lake_depth_id,'coordinates','geolon geolat')
914 CALL nc_opchk(stat, "nf90_put_att: lake_depth:coordinates")
915 stat = nf90_put_att(ncid, lake_depth_id,'long_name','lake depth')
916 CALL nc_opchk(stat, "nf90_put_att: lake_depth:long_name")
917 stat = nf90_put_att(ncid, lake_depth_id,'unit','meter')
918 CALL nc_opchk(stat, "nf90_put_att: lake_depth:long_name")
919 IF (lakedepth_srce == "GLDBV3") THEN
920 stat = nf90_put_att(ncid, lake_depth_id,'description', &
921 'based on GLDBv3 (Choulga et al. 2019); missing depth set to 10m &
922 (except to 211m in Caspian Sea); spurious large pos. depths are left unchanged.')
923 ELSE IF (lakedepth_srce == "GLDBV2") THEN
924 stat = nf90_put_att(ncid, lake_depth_id,'description', &
925 'based on GLDBv2 (Choulga et al. 2014); missing depth set to 10m &
926 (except to 211m in Caspian Sea); spurious large pos. depths are left unchanged.')
927 ELSE IF (lakedepth_srce == "GLOBATHY") THEN
928 stat = nf90_put_att(ncid, lake_depth_id,'description', &
929 'based on GLOBathy data resampled and projected to the MODIS domain.')
930 ENDIF
931 CALL nc_opchk(stat, "nf90_put_att: lake_depth:description")
932 ENDIF
933#endif
934 write(string,'(a,es8.1)') 'land_frac and lake_frac are adjusted such that &
935 their sum is 1 at points where inland=1; land_frac cutoff is',land_cutoff
936 stat = nf90_put_att(ncid, land_frac_id,'description',trim(string))
937 CALL nc_opchk(stat, "nf90_put_att: land_frac:description")
938 write(string,'(a)') 'slmsk = nint(land_frac)'
939 stat = nf90_put_att(ncid, slmsk_id,'description',trim(string))
940 CALL nc_opchk(stat, "nf90_put_att: slmsk:description")
941
942 stat = nf90_enddef(ncid)
943 CALL nc_opchk(stat, "nf90_enddef")
944
945! read in geolon and geolat and 2 variables from orog data file
946 stat = nf90_inq_varid(ncid, "geolon", geolon_id)
947 CALL nc_opchk(stat, "nf90_inq_varid: geolon")
948 stat = nf90_inq_varid(ncid, "geolat", geolat_id)
949 CALL nc_opchk(stat, "nf90_inq_varid: geolat")
950 stat = nf90_inq_varid(ncid, "land_frac", land_frac_id)
951 CALL nc_opchk(stat, "nf90_inq_varid: land_frac")
952 stat = nf90_inq_varid(ncid, "slmsk", slmsk_id)
953 CALL nc_opchk(stat, "nf90_inq_varid: slmsk")
954 stat = nf90_inq_varid(ncid, "inland", inland_id)
955 CALL nc_opchk(stat, "nf90_inq_varid: inland")
956
957 stat = nf90_get_var(ncid, geolon_id, geolon, &
958 start = (/ 1, 1 /), count = (/ tile_x_dim, tile_y_dim /) )
959 CALL nc_opchk(stat, "nf90_get_var: geolon")
960 stat = nf90_get_var(ncid, geolat_id, geolat, &
961 start = (/ 1, 1 /), count = (/ tile_x_dim, tile_y_dim /) )
962 CALL nc_opchk(stat, "nf90_get_var: geolat")
963 stat = nf90_get_var(ncid, land_frac_id, land_frac, &
964 start = (/ 1, 1 /), count = (/ tile_x_dim, tile_y_dim /) )
965 CALL nc_opchk(stat, "nf90_get_var: land_frac")
966 stat = nf90_get_var(ncid, slmsk_id, slmsk, &
967 start = (/ 1, 1 /), count = (/ tile_x_dim, tile_y_dim /) )
968 CALL nc_opchk(stat, "nf90_get_var: slmsk")
969 stat = nf90_get_var(ncid, inland_id, inland, &
970 start = (/ 1, 1 /), count = (/ tile_x_dim, tile_y_dim /) )
971 CALL nc_opchk(stat, "nf90_get_var: inland")
972
973 tile_num = 1
974 lake_frac(:) = cs_lakestat((tile_num-1)*tile_sz+1:tile_num*tile_sz)
975 lake_depth(:) = cs_lakedepth((tile_num-1)*tile_sz+1:tile_num*tile_sz)
976
977! include Caspian Sea and Aral Sea if GLDB data set is used, and
978! exclude lakes in the coastal areas of Antarctica if MODIS data set is used
979 CALL include_exclude_lakes(lake_frac,land_frac,lake_depth,geolat,geolon,tile_num)
980
981! epsil is "numerical" nonzero min for lake_frac/land_frac
982! lake_cutoff/land_cutoff is practical min for lake_frac/land_frac
983 IF (min(lake_cutoff,land_cutoff) < epsil) then
984 print *,'lake_cutoff/land_cutoff cannot be smaller than epsil, reset...'
985 lake_cutoff=max(epsil,lake_cutoff)
986 land_cutoff=max(epsil,land_cutoff)
987 ENDIF
988
989! adjust land_frac and lake_frac, and make sure land_frac+lake_frac=1 at inland points
990 DO i = 1, tile_sz
991 if (land_frac(i)< land_cutoff) land_frac(i)=0.
992 if (land_frac(i)>1.-land_cutoff) land_frac(i)=1.
993
994 if (inland(i) /= 1.) then
995 lake_frac(i) = 0.
996 endif
997
998 if (lake_frac(i) < lake_cutoff) lake_frac(i)=0.
999 if (lake_frac(i) > 1.-epsil) lake_frac(i)=1.
1000 ENDDO
1001
1002 DO i = 1, tile_sz
1003 if (inland(i) == 1.) then
1004 land_frac(i) = 1. - lake_frac(i)
1005 end if
1006
1007 if (lake_frac(i) < lake_cutoff) then
1008 lake_depth(i)=0.
1009 elseif (lake_frac(i) >= lake_cutoff .and. lake_depth(i)==0.) then
1010 lake_depth(i)=10.
1011 end if
1012 slmsk(i) = nint(land_frac(i))
1013 ENDDO
1014
1015! write 2 new variables
1016 stat = nf90_put_var(ncid, lake_frac_id, lake_frac, &
1017 start = (/ 1, 1 /), count = (/ tile_x_dim, tile_y_dim /) )
1018 CALL nc_opchk(stat, "nf90_put_var: lake_frac")
1019
1020 stat = nf90_put_var(ncid, lake_depth_id, lake_depth, &
1021 start = (/ 1, 1 /), count = (/ tile_x_dim, tile_y_dim /) )
1022 CALL nc_opchk(stat, "nf90_put_var: lake_depth")
1023
1024! write back 2 modified variables: land_frac and slmsk
1025 stat = nf90_put_var(ncid, land_frac_id, land_frac, &
1026 start = (/ 1, 1 /), count = (/ tile_x_dim, tile_y_dim /) )
1027 CALL nc_opchk(stat, "nf90_put_var: land_frac")
1028
1029 stat = nf90_put_var(ncid, slmsk_id, slmsk, &
1030 start = (/ 1, 1 /), count = (/ tile_x_dim, tile_y_dim /) )
1031 CALL nc_opchk(stat, "nf90_put_var: slmsk")
1032
1033 stat = nf90_close(ncid)
1034 CALL nc_opchk(stat, "nf90_close")
1035
1036END SUBROUTINE write_reg_lakedata_to_orodata
1037
1048SUBROUTINE include_exclude_lakes(lake_frac,land_frac,lake_depth,geolat,geolon,tile_num)
1049 REAL, INTENT(INOUT) :: lake_frac(cs_res*cs_res), lake_depth(cs_res*cs_res)
1050 REAL, INTENT(IN) :: land_frac(cs_res*cs_res)
1051 REAL, INTENT(IN) :: geolat(cs_res*cs_res), geolon(cs_res*cs_res)
1052 INTEGER, INTENT(IN) :: tile_num
1053
1054 INTEGER :: tile_sz
1055
1056 tile_sz = cs_res*cs_res
1057! add Caspian Sea and Aral Sea
1058 IF (tile_num == 2 .OR. tile_num == 3 .OR. tile_num == 7) THEN
1059 IF (lakedepth_srce == "GLDBV3" .OR. lakedepth_srce == "GLDBV2") THEN
1060 IF (lakestatus_srce == "GLDBV3" .OR. lakestatus_srce == "GLDBV2") THEN
1061 DO i = 1, tile_sz
1062 IF (land_frac(i) < 0.9 .AND. lake_frac(i) < 0.1) THEN
1063 IF (geolat(i) > 35.0 .AND. geolat(i) <= 50.0 .AND. &
1064 geolon(i) > 45.0 .AND. geolon(i) <= 55.0) THEN
1065 lake_frac(i) = 1.-land_frac(i)
1066 lake_depth(i) = 211.0
1067 ENDIF
1068 IF (geolat(i) > 35.0 .AND. geolat(i) <= 50.0 .AND. &
1069 geolon(i) > 57.0 .AND. geolon(i) <= 63.0) THEN
1070 lake_frac(i) = 1.-land_frac(i)
1071 lake_depth(i) = 10.0
1072 ENDIF
1073 ENDIF
1074 ENDDO
1075 ENDIF
1076 IF (lakestatus_srce == "MODISP" .OR. lakestatus_srce == "VIIRS") THEN
1077 DO i = 1, tile_sz
1078 IF (land_frac(i) < 0.9) THEN
1079 IF (geolat(i) > 35.0 .AND. geolat(i) <= 50.0 .AND. &
1080 geolon(i) > 45.0 .AND. geolon(i) <= 55.0) THEN
1081 lake_depth(i) = 211.0
1082 ENDIF
1083 IF (geolat(i) > 35.0 .AND. geolat(i) <= 50.0 .AND. &
1084 geolon(i) > 57.0 .AND. geolon(i) <= 63.0) THEN
1085 lake_depth(i) = 10.0
1086 ENDIF
1087 ENDIF
1088 ENDDO
1089 ENDIF
1090 ENDIF
1091 ENDIF
1092! remove lakes below 60 deg south
1093 IF (tile_num == 6) THEN
1094 IF (lakestatus_srce == "MODISP" .OR. lakestatus_srce == "VIIRS") THEN
1095 DO i = 1, tile_sz
1096 IF (geolat(i) < -60.0) THEN
1097 lake_frac(i) = 0.0
1098 lake_depth(i) = 0.0
1099 ENDIF
1100 ENDDO
1101 ENDIF
1102 ENDIF
1103
1104END SUBROUTINE include_exclude_lakes
1105
1111SUBROUTINE nc_opchk(stat,opname)
1112 USE netcdf
1113 IMPLICIT NONE
1114 INTEGER stat
1115 CHARACTER(len=*) opname
1116 CHARACTER(64) msg
1117
1118 IF (stat .NE.0) THEN
1119 msg = trim(opname) // ' Error, status code and message:'
1120 print*,trim(msg), stat, nf90_strerror(stat)
1121 stop -1
1122 END IF
1123
1124END SUBROUTINE nc_opchk
1125
1126END PROGRAM lake_frac
find_limit
subroutine find_limit(p1_in, p2_in, latmin, latmax)
Given two points on a cubed-sphere grid, compute the maximum and minimum latitudinal extent of the re...
Definition find_limit.F90:16
nc_opchk
subroutine nc_opchk(stat, opname)
Check NetCDF return code and print error message.
Definition inland.F90:468
cal_lake_frac_depth
subroutine cal_lake_frac_depth(lakestat, cs_lakestat, lakedpth, cs_lakedpth)
Calculate lake fraction and depth on the model grid from high-resolution data.
Definition lakefrac.F90:171
include_exclude_lakes
subroutine include_exclude_lakes(lake_frac, land_frac, lake_depth, geolat, geolon, tile_num)
Include Caspian Sea and Aral Sea if GLDB dataset is used, and exclude lakes in the coastal areas of A...
Definition lakefrac.F90:1049
write_lakedata_to_orodata
subroutine write_lakedata_to_orodata(cs_res, cs_lakestat, cs_lakedpth)
Write lake depth and fraction to an existing model orography file.
Definition lakefrac.F90:610
read_lakedata
subroutine read_lakedata(lakedata_path, lake_stat, lake_dpth, nlat, nlon)
Read a high-resolution lake depth dataset, and a corresponding lake status dataset which provides a s...
Definition lakefrac.F90:554
read_cubed_sphere_reg_grid
subroutine read_cubed_sphere_reg_grid(res, grid, halo_depth, res_x, res_y)
Read the latitude and longitude for a regional grid from the 'grid' file.
Definition lakefrac.F90:484
lake_cell_comp
subroutine lake_cell_comp(lkst, lkdp, lake_ct, lake_avg_frac, lake_dpth_sum)
Compute cumulatively the lake fraction and lake depth for a cell.
Definition lakefrac.F90:385
read_cubed_sphere_grid
subroutine read_cubed_sphere_grid(res, grid)
Read the latitude and longitude for a cubed-sphere grid from the 'grid' files.
Definition lakefrac.F90:424
write_reg_lakedata_to_orodata
subroutine write_reg_lakedata_to_orodata(cs_res, tile_x_dim, tile_y_dim, cs_lakestat, cs_lakedpth)
Write lake depth and fraction to an existing model orography file.
Definition lakefrac.F90:827