orog_mask_tools  1.13.0
lakefrac.F90
Go to the documentation of this file.
1 
4 
19 !#define DIAG_N_VERBOSE
20 #define ADD_ATT_FOR_NEW_VAR
21 PROGRAM 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
160 CONTAINS
161 
170 SUBROUTINE 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 
374 END SUBROUTINE cal_lake_frac_depth
375 
384 SUBROUTINE 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
412 END SUBROUTINE lake_cell_comp
413 
423 SUBROUTINE 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 
472 END SUBROUTINE read_cubed_sphere_grid
473 
483 SUBROUTINE 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 
541 END SUBROUTINE read_cubed_sphere_reg_grid
542 
553 SUBROUTINE 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 
599 END SUBROUTINE read_lakedata
600 
609 SUBROUTINE 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 
814 END SUBROUTINE write_lakedata_to_orodata
815 
826 SUBROUTINE 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 
1036 END SUBROUTINE write_reg_lakedata_to_orodata
1037 
1048 SUBROUTINE 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 
1104 END SUBROUTINE include_exclude_lakes
1105 
1111 SUBROUTINE 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 
1124 END SUBROUTINE nc_opchk
1125 
1126 END PROGRAM lake_frac
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
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
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_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 &#39;grid&#39; file.
Definition: lakefrac.F90:484
read_cubed_sphere_grid
subroutine read_cubed_sphere_grid(res, grid)
Read the latitude and longitude for a cubed-sphere grid from the &#39;grid&#39; files.
Definition: lakefrac.F90:424
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
lake_frac
program lake_frac
This program computes lake fraction and depth numbers for FV3 cubed sphere grid cells, from a high resolution lat/lon data set.
Definition: lakefrac.F90:21
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_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
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
nc_opchk
subroutine nc_opchk(stat, opname)
Check NetCDF return code and print error message.
Definition: inland.F90:468