orog_mask_tools  1.4.0
 All Data Structures Files Functions Variables
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
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 = 119.99444445
43  REAL*8, PARAMETER :: delta = 1.0 / 119.99444445
44 
45  CHARACTER(len=32) :: arg
46  CHARACTER(len=256) :: lakedata_path
47  INTEGER :: stat
48 
49  CALL getarg(0, arg) ! get the program name
50  IF (iargc() /= 3 .AND. iargc() /= 4) THEN
51  print*, 'Usage: ', trim(arg), &
52  ' [tile_num (0:all tiles, 7:regional)] [resolution (48,96, ...)] [path to lake data file]'
53  print*, 'Or: ', trim(arg), &
54  ' [tile_num (0:all tiles, 7:regional)] [resolution (48,96, ...)] [path to lake data file] [lake_cutoff]'
55  stop
56  ENDIF
57  CALL getarg(1, arg)
58  READ(arg,*,iostat=stat) tile_req
59  CALL getarg(2, arg)
60  READ(arg,*,iostat=stat) cs_res
61  CALL getarg(3, lakedata_path)
62 
63  IF (iargc() == 3) THEN
64  lake_cutoff = 0.20
65  ELSE
66  CALL getarg(4, arg)
67  READ(arg,*,iostat=stat) lake_cutoff
68  ENDIF
69 
70  IF (tile_req == 0) THEN
71  tile_beg = 1; tile_end = 6
72  print*, 'Process tile 1 - 6 at resolution C',cs_res
73  ELSE IF (tile_req /= 7) THEN
74  tile_beg = tile_req; tile_end = tile_req
75  print*, 'Process tile',tile_req, 'at resolution C',cs_res
76  ELSE
77  tile_beg = 1; tile_end = 1
78  print*, 'Process regional tile (tile', tile_req, ') at resolution C',cs_res
79  ENDIF
80 
81  ! read in grid spec data for each tile and concatenate them together
82 
83  ncsmp = (2*cs_res+1)*(2*cs_res+1)*6
84  IF (tile_req /= 7) THEN
85  print*, 'Read in cubed sphere grid information ... ',ncsmp,'pairs of lat/lons'
86  ENDIF
87 
88  IF (tile_req /= 7) THEN
89  ALLOCATE(cs_grid(ncsmp, 2))
90  CALL read_cubed_sphere_grid(cs_res, cs_grid)
91  ELSE
92  CALL read_cubed_sphere_reg_grid(cs_res, cs_grid, 3, res_x, res_y)
93  ENDIF
94  ! compute source grid
95  ALLOCATE(src_grid_lon(nlon), src_grid_lat(nlat))
96  DO i = 1, nlon
97  src_grid_lon(i) = -180.0 + delta * (i-1)
98  ENDDO
99  DO i = 1, nlat
100  src_grid_lat(i) = 90.0 - delta * (i-1)
101  ENDDO
102 
103  ! read in lake data file
104 ! sfcdata_path = '/scratch1/NCEPDEV/global/glopara/fix/fix_orog/'
105  lakedata_path = trim(lakedata_path) // "/"
106  ALLOCATE(lakestatus(nlon*nlat),lakedepth(nlon*nlat))
107  print*, 'Read in lake data file ...'
108  CALL read_lakedata(lakedata_path,lakestatus,lakedepth,nlat,nlon)
109 
110  ! calculate fraction numbers for all cs-cells
111  ncscp = cs_res*cs_res*6
112  ALLOCATE(cs_lakestatus(ncscp))
113  ALLOCATE(cs_lakedepth(ncscp))
114 
115  print*, 'Calculate lake fraction and average depth for cubed-sphere cells ...'
116  CALL cal_lake_frac_depth(lakestatus,cs_lakestatus,lakedepth,cs_lakedepth)
117 
118  ! write lake status (in terms of fraction) and lake depth(average)
119  ! to a netcdf file
120  IF (tile_req /= 7) THEN
121  print*, 'Write lake fraction/depth on cubed sphere grid cells to netCDF files ...'
122  CALL write_lakedata_to_orodata(cs_res, cs_lakestatus, cs_lakedepth)
123  ELSE
124  print*, 'Write lake fraction/depth on regional FV3 grid cells to a netCDF file ...'
125  CALL write_reg_lakedata_to_orodata(cs_res, res_x, res_y, cs_lakestatus, cs_lakedepth)
126  ENDIF
127 
128  DEALLOCATE(cs_lakestatus,cs_lakedepth)
129  DEALLOCATE(cs_grid)
130  DEALLOCATE(lakestatus,lakedepth)
131  DEALLOCATE(src_grid_lat, src_grid_lon)
132 
133  stop
134 CONTAINS
135 
144 SUBROUTINE cal_lake_frac_depth(lakestat,cs_lakestat,lakedpth,cs_lakedpth)
145  INTEGER*1, INTENT(IN) :: lakestat(:)
146  INTEGER*2, INTENT(IN) :: lakedpth(:)
147  REAL, INTENT(OUT) :: cs_lakestat(:), cs_lakedpth(:)
148 
149  REAL*8 lolf(2), lort(2), uplf(2), uprt(2), sd_ltmn(4), sd_ltmx(4)
150  REAL*8 :: v(2,4), p(2)
151  REAL :: latmin1, latmax1
152  REAL :: latmin, latmax, lonmin, lonmax, lontmp, lat_sz_max, lon_sz_max
153  INTEGER :: tile_num, i, j, gp, row, col, cs_grid_idx, cs_data_idx
154  INTEGER :: sidex_res, sidey_res, sidex_sz, sidey_sz
155  INTEGER :: stride_lat, stride_lon
156  INTEGER :: src_grid_lat_beg,src_grid_lat_end,src_grid_lon_beg,src_grid_lon_end
157  INTEGER :: src_grid_lon_beg1,src_grid_lon_end1,src_grid_lon_beg2,src_grid_lon_end2
158  INTEGER :: grid_ct, lake_ct, co_gc, tmp
159 
160  REAL*8 :: lake_dpth_sum
161  LOGICAL :: two_section, enclosure_cnvx
162 
163  IF (tile_req /= 7) THEN
164  sidex_res = cs_res; sidey_res = cs_res
165  ELSE
166  sidex_res = res_x; sidey_res = res_y
167  ENDIF
168 
169  sidex_sz = 2*sidex_res+1; sidey_sz = 2*sidey_res+1
170 
171  stride_lat = 1
172 
173  lat_sz_max = 0.0
174  lon_sz_max = 0.0
175 
176  cs_lakestat = 0.0
177 
178  DO tile_num = tile_beg, tile_end
179  row = 2 + sidex_sz*(tile_num-1); col = 2
180  DO gp = 1, sidex_res*sidey_res
181  two_section = .false.
182  cs_grid_idx = (row-1)*sidex_sz+col
183  cs_data_idx = (tile_num-1)*sidex_res*sidey_res+gp
184  IF (abs(cs_grid(cs_grid_idx,1)) > 80.0 ) THEN !ignore lakes in very high latitude
185  cs_lakestat(cs_data_idx) = 0.0
186  cs_lakedpth(cs_data_idx) = 0.0
187  ! move to next cs cell
188  col = col + 2
189  IF (col > sidex_sz) THEN
190  col = 2
191  row = row + 2
192  ENDIF
193  cycle
194  ENDIF
195  ! get the four corners of the cs cell
196  lolf(1) = cs_grid(cs_grid_idx-sidex_sz-1, 1)
197  lolf(2) = cs_grid(cs_grid_idx-sidex_sz-1, 2)
198  IF (lolf(2) > 180.0) lolf(2) = lolf(2) - 360.0
199  lort(1) = cs_grid(cs_grid_idx-sidex_sz+1, 1)
200  lort(2) = cs_grid(cs_grid_idx-sidex_sz+1, 2)
201  IF (lort(2) > 180.0) lort(2) = lort(2) - 360.0
202  uplf(1) = cs_grid(cs_grid_idx+sidex_sz-1,1)
203  uplf(2) = cs_grid(cs_grid_idx+sidex_sz-1,2)
204  IF (uplf(2) > 180.0) uplf(2) = uplf(2) - 360.0
205  uprt(1) = cs_grid(cs_grid_idx+sidex_sz+1,1)
206  uprt(2) = cs_grid(cs_grid_idx+sidex_sz+1,2)
207 
208  v(1,1) = lolf(1); v(2,1) = lolf(2)
209  v(1,2) = lort(1); v(2,2) = lort(2)
210  v(1,3) = uprt(1); v(2,3) = uprt(2)
211  v(1,4) = uplf(1); v(2,4) = uplf(2)
212  v(:,:) = v(:,:) * d2r
213 
214  IF (uprt(2) > 180.0) uprt(2) = uprt(2) - 360.0
215  ! gather the candidate indices in lakestat
216 #ifdef LIMIT_CAL
217  CALL find_limit(lolf, lort, sd_ltmn(1), sd_ltmx(1))
218  CALL find_limit(lort, uprt, sd_ltmn(2), sd_ltmx(2))
219  CALL find_limit(uprt, uplf, sd_ltmn(3), sd_ltmx(3))
220  CALL find_limit(uplf, lolf, sd_ltmn(4), sd_ltmx(4))
221  latmin = min(sd_ltmn(1),min(sd_ltmn(2),min(sd_ltmn(3),sd_ltmn(4))))
222  latmax = max(sd_ltmx(1),max(sd_ltmx(2),max(sd_ltmx(3),sd_ltmx(4))))
223 #endif
224  latmin = min(lolf(1),min(lort(1),min(uplf(1),uprt(1))))
225  latmax = max(lolf(1),max(lort(1),max(uplf(1),uprt(1))))
226  lonmin = min(lolf(2),min(lort(2),min(uplf(2),uprt(2))))
227  lonmax = max(lolf(2),max(lort(2),max(uplf(2),uprt(2))))
228 ! lat_sz_max = max(lat_sz_max, (latmax-latmin))
229 ! lon_sz_max = max(lon_sz_max, (lonmax-lonmin))
230 
231  src_grid_lat_beg = nint((90.0-latmax)*gppdeg+0.5)
232  src_grid_lat_end = nint((90.0-latmin)*gppdeg+0.5)
233  src_grid_lon_beg = nint((180.0+lonmin)*gppdeg+0.5)
234  src_grid_lon_end = nint((180.0+lonmax)*gppdeg+0.5)
235 
236  IF (src_grid_lat_beg > src_grid_lat_end) THEN
237  print*,'switch!!!'
238  tmp = src_grid_lat_beg
239  src_grid_lat_beg = src_grid_lat_end
240  src_grid_lat_end = tmp
241  ENDIF
242  IF (src_grid_lon_beg > src_grid_lon_end) THEN
243  print*,'switch!!!'
244  tmp = src_grid_lon_beg
245  src_grid_lon_beg = src_grid_lon_end
246  src_grid_lon_end = tmp
247  ENDIF
248  IF ((src_grid_lon_end - src_grid_lon_beg) > nlon*0.75) THEN
249  two_section = .true.
250  src_grid_lon_beg1 = src_grid_lon_end
251  src_grid_lon_end1 = nlon
252  src_grid_lon_beg2 = 1
253  src_grid_lon_end2 = src_grid_lon_beg
254  ENDIF
255 
256 #ifdef DIAG_N_VERBOSE
257  print*, 'cell centre lat/lon =', &
258  gp, cs_res*cs_res, cs_grid(cs_grid_idx,1),cs_grid(cs_grid_idx,2)
259  print*, 'lat index range and stride', &
260  src_grid_lat_beg,src_grid_lat_end,stride_lat
261  print*, 'lat range ', &
262  src_grid_lat(src_grid_lat_beg),src_grid_lat(src_grid_lat_end)
263 #endif
264  lake_ct = 0; grid_ct = 0
265  lake_dpth_sum = 0.0
266  DO j = src_grid_lat_beg, src_grid_lat_end, stride_lat
267  stride_lon = int(1.0/cos(src_grid_lat(j)*d2r)*REAL(stride_lat))
268 #ifdef DIAG_N_VERBOSE
269  IF (j == src_grid_lat_beg) THEN
270  print*, 'lon index range and stride', &
271  src_grid_lon_beg,src_grid_lon_end,stride_lon
272  print*, 'lon range ', &
273  src_grid_lon(src_grid_lon_beg),src_grid_lon(src_grid_lon_end)
274  IF (two_section == .true.) THEN
275  print*, 'section1 index lon range and stride', &
276  src_grid_lon_beg1,src_grid_lon_end1,stride_lon
277  print*, 'section1 lon range ', &
278  src_grid_lon(src_grid_lon_beg1),src_grid_lon(src_grid_lon_end1)
279  print*, 'section2 index lon range and stride', &
280  src_grid_lon_beg2,src_grid_lon_end2,stride_lon
281  print*, 'section2 lon range ', &
282  src_grid_lon(src_grid_lon_beg2),src_grid_lon(src_grid_lon_end2)
283  ENDIF
284  ENDIF
285 #endif
286  IF (two_section .EQV. .false.) THEN
287  DO i = src_grid_lon_beg, src_grid_lon_end, stride_lon
288  p(1) = src_grid_lat(j); p(2) = src_grid_lon(i)
289  p(:) = p(:)*d2r
290  IF(enclosure_cnvx(v, 4, p, co_gc) .EQV. .true.) THEN
291  grid_ct = grid_ct+1
292  IF (lakestat((j-1)*nlon+i) /= 0) THEN
293  lake_ct = lake_ct+1
294  IF (lakedpth((j-1)*nlon+i) < 0) THEN
295  IF (lakestat((j-1)*nlon+i) == 4) THEN
296  lake_dpth_sum = lake_dpth_sum+30.0
297  ELSE
298  lake_dpth_sum = lake_dpth_sum+100.0
299  ENDIF
300  ELSE
301  lake_dpth_sum = lake_dpth_sum+REAL(lakedpth((j-1)*nlon+i))
302  ENDIF
303  ENDIF
304  ENDIF
305  ENDDO
306  ELSE
307  DO i = src_grid_lon_beg1, src_grid_lon_end1, stride_lon
308  p(1) = src_grid_lat(j); p(2) = src_grid_lon(i)
309  p(:) = p(:)*d2r
310  IF(enclosure_cnvx(v, 4, p, co_gc) .EQV. .true.) THEN
311  grid_ct = grid_ct+1
312  IF (lakestat((j-1)*nlon+i) /= 0) THEN
313  lake_ct = lake_ct+1
314  IF (lakedpth((j-1)*nlon+i) < 0) THEN
315  IF (lakestat((j-1)*nlon+i) == 4) THEN
316  lake_dpth_sum = lake_dpth_sum+30.0
317  ELSE
318  lake_dpth_sum = lake_dpth_sum+100.0
319  ENDIF
320  ELSE
321  lake_dpth_sum = lake_dpth_sum+REAL(lakedpth((j-1)*nlon+i))
322  ENDIF
323  ENDIF
324  ENDIF
325  ENDDO
326  DO i = src_grid_lon_beg2, src_grid_lon_end2, stride_lon
327  p(1) = src_grid_lat(j); p(2) = src_grid_lon(i)
328  p(:) = p(:)*d2r
329  IF(enclosure_cnvx(v, 4, p, co_gc) .EQV. .true.) THEN
330  grid_ct = grid_ct+1
331  IF (lakestat((j-1)*nlon+i) /= 0) THEN
332  lake_ct = lake_ct+1
333  IF (lakedpth((j-1)*nlon+i) < 0) THEN
334  IF (lakestat((j-1)*nlon+i) == 4) THEN
335  lake_dpth_sum = lake_dpth_sum+30.0
336  ELSE
337  lake_dpth_sum = lake_dpth_sum+100.0
338  ENDIF
339  ELSE
340  lake_dpth_sum = lake_dpth_sum+REAL(lakedpth((j-1)*nlon+i))
341  ENDIF
342  ENDIF
343  ENDIF
344  ENDDO
345  ENDIF
346  ENDDO
347  cs_lakestat(cs_data_idx)=REAL(lake_ct)/REAL(grid_ct)
348  IF (lake_ct /= 0) THEN
349  cs_lakedpth(cs_data_idx)=lake_dpth_sum/REAL(lake_ct)/10.0 !convert to meter
350  ELSE
351  cs_lakedpth(cs_data_idx)=0.0
352  ENDIF
353 #ifdef DIAG_N_VERBOSE
354  print*, 'tile_num, row, col:', tile_num, row, col
355  print*, 'grid_ct, lake_ct = ', grid_ct, lake_ct
356  print*, 'lake_frac= ', cs_lakestat(cs_data_idx)
357  print*, 'lake_depth (avg) = ', cs_lakedpth(cs_data_idx)
358 #endif
359 
360  ! move to the next control volume
361  col = col + 2
362  IF (col > sidex_sz) THEN
363  col = 2
364  row = row + 2
365  ENDIF
366  ENDDO
367  print "('*'$)" ! progress '*'
368  ENDDO
369  print*, ''
370 
371 END SUBROUTINE cal_lake_frac_depth
372 
382 SUBROUTINE read_cubed_sphere_grid(res, grid)
383  INTEGER, INTENT(IN) :: res
384  REAL, INTENT(OUT) :: grid(:,:)
385 
386  REAL*8, ALLOCATABLE :: lat(:), lon(:)
387  INTEGER :: side_sz, tile_sz, ncid, varid
388  INTEGER :: i, tile_beg, tile_end, stat
389  CHARACTER(len=256) :: gridfile_path,gridfile
390  CHARACTER(len=1) ich
391  CHARACTER(len=4) res_ch
392 
393  side_sz = 2*res+1
394  tile_sz = side_sz*side_sz
395  ALLOCATE(lat(tile_sz), lon(tile_sz))
396 
397  IF (tile_req == 0) THEN
398  tile_beg = 1; tile_end = 6
399  ELSE
400  tile_beg = tile_req; tile_end = tile_req
401  ENDIF
402  WRITE(res_ch,'(I4)') res
403 ! gridfile_path = "/scratch1/NCEPDEV/global/glopara/fix/fix_fv3/C"//trim(adjustl(res_ch))//"/"
404  gridfile_path = "./"
405  DO i = tile_beg, tile_end
406  WRITE(ich, '(I1)') i
407  gridfile = trim(gridfile_path)//"C"//trim(adjustl(res_ch))//"_grid.tile"//ich//".nc"
408 
409  print*, 'Open cubed sphere grid spec file ', trim(gridfile)
410 
411  stat = nf90_open(trim(gridfile), nf90_nowrite, ncid)
412  CALL nc_opchk(stat,'nf90_open')
413 
414  stat = nf90_inq_varid(ncid,'x',varid)
415  CALL nc_opchk(stat,'nf90_inq_lon')
416  stat = nf90_get_var(ncid,varid,lon,start=(/1,1/),count=(/side_sz,side_sz/))
417  CALL nc_opchk(stat,'nf90_get_var_lon')
418  stat = nf90_inq_varid(ncid,'y',varid)
419  CALL nc_opchk(stat,'nf90_inq_lat')
420  stat = nf90_get_var(ncid,varid,lat,start=(/1,1/),count=(/side_sz,side_sz/))
421  CALL nc_opchk(stat,'nf90_get_var_lat')
422  stat = nf90_close(ncid)
423  CALL nc_opchk(stat,'nf90_close')
424 
425  tile_beg = (i-1)*tile_sz+1
426  tile_end = i*tile_sz
427  grid(tile_beg:tile_end,1) = lat(1:tile_sz)
428  grid(tile_beg:tile_end,2) = lon(1:tile_sz)
429  END DO
430  DEALLOCATE(lat,lon)
431 
432 END SUBROUTINE read_cubed_sphere_grid
433 
443 SUBROUTINE read_cubed_sphere_reg_grid(res, grid, halo_depth, res_x, res_y)
444  INTEGER, INTENT(IN) :: res, halo_depth
445  INTEGER, INTENT(OUT) :: res_x, res_y
446  REAL, ALLOCATABLE, INTENT(OUT) :: grid(:,:)
447 
448  REAL*8, ALLOCATABLE :: lat(:), lon(:)
449  INTEGER :: sidex_sz, sidey_sz, tile_sz, ncid, varid, dimid
450  INTEGER :: x_start, y_start
451  INTEGER :: nxp, nyp, stat
452  CHARACTER(len=256) :: gridfile_path,gridfile
453  CHARACTER(len=1) ich
454  CHARACTER(len=4) res_ch
455  CHARACTER(len=8) dimname
456 
457  WRITE(res_ch,'(I4)') res
458  gridfile_path = './'
459  gridfile = trim(gridfile_path)//"C"//trim(adjustl(res_ch))//"_grid.tile7.nc"
460 
461  print*, 'Open cubed sphere grid spec file ', trim(gridfile)
462 
463  stat = nf90_open(trim(gridfile), nf90_nowrite, ncid)
464  CALL nc_opchk(stat,'nf90_open')
465 
466  stat = nf90_inq_dimid(ncid,'nxp',dimid)
467  CALL nc_opchk(stat,'nf90_inq_dimid: nxp')
468  stat = nf90_inquire_dimension(ncid,dimid,dimname,len=nxp)
469  CALL nc_opchk(stat,'nf90_inquire_dimension: nxp')
470 
471  stat = nf90_inq_dimid(ncid,'nyp',dimid)
472  CALL nc_opchk(stat,'nf90_inq_dimid: nyp')
473  stat = nf90_inquire_dimension(ncid,dimid,dimname,len=nyp)
474  CALL nc_opchk(stat,'nf90_inquire_dimension: nyp')
475 
476  sidex_sz = nxp
477  sidey_sz = nyp
478  tile_sz = sidex_sz*sidey_sz
479  ALLOCATE(lat(tile_sz), lon(tile_sz))
480 
481 ! x_start = halo_depth+1; y_start = halo_depth+1
482  x_start = 1; y_start = 1
483  stat = nf90_inq_varid(ncid,'x',varid)
484  CALL nc_opchk(stat,'nf90_inq_lon')
485  stat = nf90_get_var(ncid,varid,lon,start=(/x_start,y_start/),count=(/sidex_sz,sidey_sz/))
486  CALL nc_opchk(stat,'nf90_get_var_lon')
487  stat = nf90_inq_varid(ncid,'y',varid)
488  CALL nc_opchk(stat,'nf90_inq_lat')
489  stat = nf90_get_var(ncid,varid,lat,start=(/x_start,y_start/),count=(/sidex_sz,sidey_sz/))
490  CALL nc_opchk(stat,'nf90_get_var_lat')
491  stat = nf90_close(ncid)
492  CALL nc_opchk(stat,'nf90_close')
493 
494  ALLOCATE(grid(tile_sz,2))
495  grid(1:tile_sz,1) = lat(1:tile_sz)
496  grid(1:tile_sz,2) = lon(1:tile_sz)
497 
498  res_x = sidex_sz/2; res_y = sidey_sz/2
499  DEALLOCATE(lat,lon)
500 
501 END SUBROUTINE read_cubed_sphere_reg_grid
502 
513 SUBROUTINE read_lakedata(lakedata_path,lake_stat,lake_dpth,nlat,nlon)
514  CHARACTER(len=256), INTENT(IN) :: lakedata_path
515  INTEGER*1, INTENT(OUT) :: lake_stat(:)
516  INTEGER*2, INTENT(OUT) :: lake_dpth(:)
517  INTEGER, INTENT(IN) :: nlat, nlon
518 
519  CHARACTER(len=256) lakefile
520  INTEGER :: data_sz, i
521 
522  data_sz = nlon*nlat
523 
524  ! read in 30 arc seconds lake data
525  lakefile = trim(lakedata_path) // "GlobalLakeStatus.dat"
526  OPEN(10, file=lakefile, form='unformatted', access='direct', recl=data_sz*1)
527  READ(10,rec=1) lake_stat
528  CLOSE(10)
529 
530  lakefile = trim(lakedata_path) // "GlobalLakeDepth.dat"
531  OPEN(10, file=lakefile, form='unformatted', access='direct', recl=data_sz*2)
532  READ(10,rec=1) lake_dpth
533  CLOSE(10)
534 
535 END SUBROUTINE read_lakedata
536 
545 SUBROUTINE write_lakedata_to_orodata(cs_res, cs_lakestat, cs_lakedpth)
546  USE netcdf
547  INTEGER, INTENT(IN) :: cs_res
548  REAL, INTENT(IN) :: cs_lakestat(:)
549  REAL, INTENT(IN) :: cs_lakedpth(:)
550 
551  INTEGER :: tile_sz, tile_num
552  INTEGER :: stat, ncid, x_dimid, y_dimid, varid, dimids(2)
553  INTEGER :: lake_frac_id, lake_depth_id
554  INTEGER :: land_frac_id, slmsk_id, inland_id, geolon_id, geolat_id
555  CHARACTER(len=256) :: filename,string
556  CHARACTER(len=1) :: ich
557  CHARACTER(len=4) res_ch
558  REAL :: lake_frac(cs_res*cs_res),lake_depth(cs_res*cs_res)
559  REAL :: geolon(cs_res*cs_res),geolat(cs_res*cs_res)
560  REAL :: land_frac(cs_res*cs_res),slmsk(cs_res*cs_res),inland(cs_res*cs_res)
561  real, parameter :: epsil=1.e-6 ! numerical min for lake_frac/land_frac
562  real :: land_cutoff=1.e-4 ! land_frac=0 if it is < land_cutoff
563 
564  INTEGER :: i, j
565 
566 ! include "netcdf.inc"
567  tile_sz = cs_res*cs_res
568 
569  WRITE(res_ch,'(I4)') cs_res
570  DO tile_num = tile_beg, tile_end
571  WRITE(ich, '(I1)') tile_num
572 ! filename = "C" // trim(adjustl(res_ch)) // "_oro_data.tile" // ich // ".nc"
573 ! filename = "oro_data.tile" // ich // ".nc"
574  filename = "oro.C" // trim(adjustl(res_ch)) // ".tile" // ich // ".nc"
575  print *,'Read, update, and write ',trim(filename)
576  stat = nf90_open(filename, nf90_write, ncid)
577  CALL nc_opchk(stat, "nf90_open oro_data.nc")
578  stat = nf90_inq_dimid(ncid, "lon", x_dimid)
579  CALL nc_opchk(stat, "nf90_inq_dim: x")
580  stat = nf90_inq_dimid(ncid, "lat", y_dimid)
581  CALL nc_opchk(stat, "nf90_inq_dim: y")
582 ! dimids = (/ y_dimid, x_dimid /)
583 ! original orodata netcdf file uses (y, x) order, so we made change to match it.
584  dimids = (/ x_dimid, y_dimid /)
585  stat = nf90_inq_varid(ncid, "land_frac", land_frac_id)
586  CALL nc_opchk(stat, "nf90_inq_varid: land_frac")
587  stat = nf90_inq_varid(ncid, "slmsk", slmsk_id)
588  CALL nc_opchk(stat, "nf90_inq_varid: slmsk")
589 ! define 2 new variables
590  stat = nf90_redef(ncid)
591  CALL nc_opchk(stat, "nf90_redef")
592  stat = nf90_def_var(ncid,"lake_frac",nf90_float,dimids,lake_frac_id)
593  CALL nc_opchk(stat, "nf90_def_var: lake_frac")
594 #ifdef ADD_ATT_FOR_NEW_VAR
595  stat = nf90_put_att(ncid, lake_frac_id,'coordinates','geolon geolat')
596  CALL nc_opchk(stat, "nf90_put_att: lake_frac:coordinates")
597  stat = nf90_put_att(ncid, lake_frac_id,'long_name','lake fraction')
598  CALL nc_opchk(stat, "nf90_put_att: lake_frac:long_name")
599  stat = nf90_put_att(ncid, lake_frac_id,'unit','fraction')
600  CALL nc_opchk(stat, "nf90_put_att: lake_frac:unit")
601  write(string,'(a,f5.2)') 'based on GLDBv2 (Choulga et al. 2014); missing lakes &
602  added based on land_frac in this dataset; lake_frac cutoff is',lake_cutoff
603  stat = nf90_put_att(ncid, lake_frac_id,'description',trim(string))
604  CALL nc_opchk(stat, "nf90_put_att: lake_frac:description")
605 #endif
606  stat = nf90_def_var(ncid,"lake_depth",nf90_float,dimids,lake_depth_id)
607  CALL nc_opchk(stat, "nf90_def_var: lake_depth")
608 #ifdef ADD_ATT_FOR_NEW_VAR
609  stat = nf90_put_att(ncid, lake_depth_id,'coordinates','geolon geolat')
610  CALL nc_opchk(stat, "nf90_put_att: lake_depth:coordinates")
611  stat = nf90_put_att(ncid, lake_depth_id,'long_name','lake depth')
612  CALL nc_opchk(stat, "nf90_put_att: lake_depth:long_name")
613  stat = nf90_put_att(ncid, lake_depth_id,'unit','meter')
614  CALL nc_opchk(stat, "nf90_put_att: lake_depth:long_name")
615  stat = nf90_put_att(ncid, lake_depth_id,'description', &
616  'based on GLDBv2 (Choulga et al. 2014); missing depth set to 10m &
617  (except to 211m in Caspian Sea); spurious large pos. depths are left unchanged.')
618  CALL nc_opchk(stat, "nf90_put_att: lake_depth:description")
619 #endif
620 
621  write(string,'(a,es8.1)') 'land_frac and lake_frac are adjusted such that their sum '// &
622  'is 1 at points where inland=1; land_frac cutoff is',land_cutoff
623  stat = nf90_put_att(ncid, land_frac_id,'description',trim(string))
624  CALL nc_opchk(stat, "nf90_put_att: land_frac:description")
625 
626  write(string,'(a)') 'slmsk = nint(land_frac)'
627  stat = nf90_put_att(ncid, slmsk_id,'description',trim(string))
628  CALL nc_opchk(stat, "nf90_put_att: slmsk:description")
629 
630  stat = nf90_enddef(ncid)
631  CALL nc_opchk(stat, "nf90_enddef")
632 
633 ! read in geolon and geolat and 2 variables from orog data file
634  stat = nf90_inq_varid(ncid, "geolon", geolon_id)
635  CALL nc_opchk(stat, "nf90_inq_varid: geolon")
636  stat = nf90_inq_varid(ncid, "geolat", geolat_id)
637  CALL nc_opchk(stat, "nf90_inq_varid: geolat")
638  stat = nf90_inq_varid(ncid, "land_frac", land_frac_id)
639  CALL nc_opchk(stat, "nf90_inq_varid: land_frac")
640  stat = nf90_inq_varid(ncid, "slmsk", slmsk_id)
641  CALL nc_opchk(stat, "nf90_inq_varid: slmsk")
642  stat = nf90_inq_varid(ncid, "inland", inland_id)
643  CALL nc_opchk(stat, "nf90_inq_varid: inland")
644  stat = nf90_get_var(ncid, geolon_id, geolon, &
645  start = (/ 1, 1 /), count = (/ cs_res, cs_res /) )
646  CALL nc_opchk(stat, "nf90_get_var: geolon")
647  stat = nf90_get_var(ncid, geolat_id, geolat, &
648  start = (/ 1, 1 /), count = (/ cs_res, cs_res /) )
649  CALL nc_opchk(stat, "nf90_get_var: geolat")
650  stat = nf90_get_var(ncid, land_frac_id, land_frac, &
651  start = (/ 1, 1 /), count = (/ cs_res, cs_res /) )
652  CALL nc_opchk(stat, "nf90_get_var: land_frac")
653  stat = nf90_get_var(ncid, slmsk_id, slmsk, &
654  start = (/ 1, 1 /), count = (/ cs_res, cs_res /) )
655  CALL nc_opchk(stat, "nf90_get_var: slmsk")
656  stat = nf90_get_var(ncid, inland_id, inland, &
657  start = (/ 1, 1 /), count = (/ cs_res, cs_res /) )
658  CALL nc_opchk(stat, "nf90_get_var: inland")
659 
660  lake_frac(:) = cs_lakestat((tile_num-1)*tile_sz+1:tile_num*tile_sz)
661  lake_depth(:) = cs_lakedepth((tile_num-1)*tile_sz+1:tile_num*tile_sz)
662 
663 ! add Caspian Sea and Aral Sea to lake_frac and lake_depth
664  IF (tile_num == 2 .or. tile_num == 3) THEN
665  DO i = 1, tile_sz
666  IF (land_frac(i) < 0.9 .AND. lake_frac(i) < 0.1) THEN
667  IF (geolat(i) > 35.0 .AND. geolat(i) <= 50.0 .AND. &
668  geolon(i) > 45.0 .AND. geolon(i) <= 55.0) THEN
669  lake_frac(i) = 1.-land_frac(i)
670  lake_depth(i) = 211.0
671  ENDIF
672  IF (geolat(i) > 35.0 .AND. geolat(i) <= 50.0 .AND. &
673  geolon(i) > 57.0 .AND. geolon(i) <= 63.0) THEN
674  lake_frac(i) = 1.-land_frac(i)
675  lake_depth(i) = 10.0
676  ENDIF
677  ENDIF
678  ENDDO
679  ENDIF
680 
681 ! adjust land_frac and lake_frac, and make sure land_frac+lake_frac=1 at inland points
682  DO i = 1, tile_sz
683  if (lake_frac(i) >= lake_cutoff) then ! non-zero lake_frac dominates over land_frac
684  land_frac(i) = max(0., min(1., 1.-lake_frac(i)))
685  elseif (inland(i) == 1.) then ! land_frac dominates over lake_frac at inland points
686  lake_frac(i) = max(0., min(1., 1.-land_frac(i)))
687  end if
688 
689 ! epsil is "numerical" nonzero min for lake_frac/land_frac
690 ! lake_cutoff/land_cutoff is practical min for lake_frac/land_frac
691  if (min(lake_cutoff,land_cutoff) < epsil) then
692  print *,'lake_cutoff/land_cutoff cannot be smaller than epsil, reset...'
693  lake_cutoff=max(epsil,lake_cutoff)
694  land_cutoff=max(epsil,land_cutoff)
695  end if
696 
697  if (lake_frac(i)< lake_cutoff) lake_frac(i)=0.
698  if (lake_frac(i)>1.-lake_cutoff) lake_frac(i)=1.
699  if (inland(i) == 1.) land_frac(i) = 1.-lake_frac(i)
700  if (land_frac(i)< land_cutoff) land_frac(i)=0.
701  if (land_frac(i)>1.-land_cutoff) land_frac(i)=1.
702 
703  if (lake_frac(i) < lake_cutoff) then
704  lake_depth(i)=0.
705  elseif (lake_frac(i) >= lake_cutoff .and. lake_depth(i)==0.) then
706  lake_depth(i)=10.
707  end if
708  slmsk(i) = nint(land_frac(i))
709  ENDDO
710 ! write 2 new variables
711  stat = nf90_put_var(ncid, lake_frac_id, lake_frac, &
712  start = (/ 1, 1 /), count = (/ cs_res, cs_res /) )
713  CALL nc_opchk(stat, "nf90_put_var: lake_frac")
714 
715  stat = nf90_put_var(ncid, lake_depth_id, lake_depth, &
716  start = (/ 1, 1 /), count = (/ cs_res, cs_res /) )
717  CALL nc_opchk(stat, "nf90_put_var: lake_depth")
718 
719 ! write back 2 modified variables: land_frac and slmsk
720  stat = nf90_put_var(ncid, land_frac_id, land_frac, &
721  start = (/ 1, 1 /), count = (/ cs_res, cs_res /) )
722  CALL nc_opchk(stat, "nf90_put_var: land_frac")
723 
724  stat = nf90_put_var(ncid, slmsk_id, slmsk, &
725  start = (/ 1, 1 /), count = (/ cs_res, cs_res /) )
726  CALL nc_opchk(stat, "nf90_put_var: slmsk")
727 
728  stat = nf90_close(ncid)
729  CALL nc_opchk(stat, "nf90_close")
730  ENDDO
731 
732 END SUBROUTINE write_lakedata_to_orodata
733 
744 SUBROUTINE write_reg_lakedata_to_orodata(cs_res, tile_x_dim, tile_y_dim, cs_lakestat, cs_lakedpth)
745  USE netcdf
746  INTEGER, INTENT(IN) :: cs_res, tile_x_dim, tile_y_dim
747  REAL, INTENT(IN) :: cs_lakestat(:)
748  REAL, INTENT(IN) :: cs_lakedpth(:)
749 
750  INTEGER :: tile_sz, tile_num
751  INTEGER :: stat, ncid, x_dimid, y_dimid, varid, dimids(2)
752  INTEGER :: lake_frac_id, lake_depth_id
753  INTEGER :: land_frac_id, slmsk_id, geolon_id, geolat_id, inland_id
754  CHARACTER(len=256) :: filename,string
755  CHARACTER(len=1) :: ich
756  CHARACTER(len=4) res_ch
757 
758  REAL, ALLOCATABLE :: lake_frac(:), lake_depth(:), geolon(:), geolat(:)
759  REAL, ALLOCATABLE :: land_frac(:), slmsk(:), inland(:)
760 
761  real, parameter :: epsil=1.e-6 ! numerical min for lake_frac/land_frac
762  real :: land_cutoff=1.e-6 ! land_frac=0 if it is < land_cutoff
763 
764  INTEGER :: i, j, var_id
765 
766 ! include "netcdf.inc"
767  tile_sz = tile_x_dim*tile_y_dim
768 
769  ALLOCATE(lake_frac(tile_sz), lake_depth(tile_sz))
770  ALLOCATE(geolon(tile_sz), geolat(tile_sz))
771  ALLOCATE(land_frac(tile_sz), slmsk(tile_sz), inland(tile_sz))
772 
773  WRITE(res_ch,'(I4)') cs_res
774  tile_num = 7
775  WRITE(ich, '(I1)') tile_num
776 ! filename = "C" // trim(adjustl(res_ch)) // "_oro_data.tile" // ich // ".halo4.nc"
777  filename = "oro.C" // trim(adjustl(res_ch)) // ".tile" // ich // ".nc"
778  print*, 'Open and write regional orography data netCDF file ', trim(filename)
779  stat = nf90_open(filename, nf90_write, ncid)
780  CALL nc_opchk(stat, "nf90_open oro_data.nc")
781  stat = nf90_inq_dimid(ncid, "lon", x_dimid)
782  CALL nc_opchk(stat, "nf90_inq_dim: x")
783  stat = nf90_inq_dimid(ncid, "lat", y_dimid)
784  CALL nc_opchk(stat, "nf90_inq_dim: y")
785  dimids = (/ x_dimid, y_dimid /)
786 
787  stat = nf90_inq_varid(ncid, "land_frac", land_frac_id)
788  CALL nc_opchk(stat, "nf90_inq_varid: land_frac")
789  stat = nf90_inq_varid(ncid, "slmsk", slmsk_id)
790  CALL nc_opchk(stat, "nf90_inq_varid: slmsk")
791 
792 ! define 2 new variables, lake_frac and lake_depth
793  stat = nf90_redef(ncid)
794  CALL nc_opchk(stat, "nf90_redef")
795 
796  IF (nf90_inq_varid(ncid, "lake_frac",lake_frac_id) /= 0) THEN
797  stat = nf90_def_var(ncid,"lake_frac",nf90_float,dimids,lake_frac_id)
798  CALL nc_opchk(stat, "nf90_def_var: lake_frac")
799 #ifdef ADD_ATT_FOR_NEW_VAR
800  stat = nf90_put_att(ncid, lake_frac_id,'coordinates','geolon geolat')
801  CALL nc_opchk(stat, "nf90_put_att: lake_frac:coordinates")
802  stat = nf90_put_att(ncid, lake_frac_id,'long_name','lake fraction')
803  CALL nc_opchk(stat, "nf90_put_att: lake_frac:long_name")
804  stat = nf90_put_att(ncid, lake_frac_id,'unit','fraction')
805  CALL nc_opchk(stat, "nf90_put_att: lake_frac:unit")
806  stat = nf90_put_att(ncid, lake_frac_id,'description', &
807  'based on GLDBv2 (Choulga et al. 2014); missing lakes &
808  added based on land_frac in this dataset.')
809  CALL nc_opchk(stat, "nf90_put_att: lake_frac:description")
810 #endif
811  ENDIF
812  IF (nf90_inq_varid(ncid, "lake_depth",lake_depth_id) /= 0) THEN
813  stat = nf90_def_var(ncid,"lake_depth",nf90_float,dimids,lake_depth_id)
814  CALL nc_opchk(stat, "nf90_def_var: lake_depth")
815 #ifdef ADD_ATT_FOR_NEW_VAR
816  stat = nf90_put_att(ncid, lake_depth_id,'coordinates','geolon geolat')
817  CALL nc_opchk(stat, "nf90_put_att: lake_depth:coordinates")
818  stat = nf90_put_att(ncid, lake_depth_id,'long_name','lake depth')
819  CALL nc_opchk(stat, "nf90_put_att: lake_depth:long_name")
820  stat = nf90_put_att(ncid, lake_depth_id,'unit','meter')
821  CALL nc_opchk(stat, "nf90_put_att: lake_depth:long_name")
822  stat = nf90_put_att(ncid, lake_depth_id,'description', &
823  'based on GLDBv2 (Choulga et al. 2014); missing depth set to 10m &
824  (except to 211m in Caspian Sea); spurious large pos. depths are left unchanged.')
825  CALL nc_opchk(stat, "nf90_put_att: lake_depth:description")
826 #endif
827  ENDIF
828  write(string,'(a,es8.1)') 'land_frac and lake_frac are adjusted '// &
829  'such that their sum is 1 at points where inland=1; land_frac '// &
830  'cutoff is ',land_cutoff
831  stat = nf90_put_att(ncid, land_frac_id,'description',trim(string))
832  CALL nc_opchk(stat, "nf90_put_att: land_frac:description")
833 
834  write(string,'(a)') 'slmsk = nint(land_frac)'
835  stat = nf90_put_att(ncid, slmsk_id,'description',trim(string))
836  CALL nc_opchk(stat, "nf90_put_att: slmsk:description")
837 
838  stat = nf90_enddef(ncid)
839  CALL nc_opchk(stat, "nf90_enddef")
840 
841 ! read in geolon and geolat and 2 variables from orog data file
842  stat = nf90_inq_varid(ncid, "geolon", geolon_id)
843  CALL nc_opchk(stat, "nf90_inq_varid: geolon")
844  stat = nf90_inq_varid(ncid, "geolat", geolat_id)
845  CALL nc_opchk(stat, "nf90_inq_varid: geolat")
846  stat = nf90_inq_varid(ncid, "land_frac", land_frac_id)
847  CALL nc_opchk(stat, "nf90_inq_varid: land_frac")
848  stat = nf90_inq_varid(ncid, "slmsk", slmsk_id)
849  CALL nc_opchk(stat, "nf90_inq_varid: slmsk")
850  stat = nf90_inq_varid(ncid, "inland", inland_id)
851  CALL nc_opchk(stat, "nf90_inq_varid: inland")
852 
853  stat = nf90_get_var(ncid, geolon_id, geolon, &
854  start = (/ 1, 1 /), count = (/ tile_x_dim, tile_y_dim /) )
855  CALL nc_opchk(stat, "nf90_get_var: geolon")
856  stat = nf90_get_var(ncid, geolat_id, geolat, &
857  start = (/ 1, 1 /), count = (/ tile_x_dim, tile_y_dim /) )
858  CALL nc_opchk(stat, "nf90_get_var: geolat")
859  stat = nf90_get_var(ncid, land_frac_id, land_frac, &
860  start = (/ 1, 1 /), count = (/ tile_x_dim, tile_y_dim /) )
861  CALL nc_opchk(stat, "nf90_get_var: land_frac")
862  stat = nf90_get_var(ncid, slmsk_id, slmsk, &
863  start = (/ 1, 1 /), count = (/ tile_x_dim, tile_y_dim /) )
864  CALL nc_opchk(stat, "nf90_get_var: slmsk")
865  stat = nf90_get_var(ncid, inland_id, inland, &
866  start = (/ 1, 1 /), count = (/ tile_x_dim, tile_y_dim /) )
867  CALL nc_opchk(stat, "nf90_get_var: inland")
868 
869  tile_num = 1
870  lake_frac(:) = cs_lakestat((tile_num-1)*tile_sz+1:tile_num*tile_sz)
871  lake_depth(:) = cs_lakedepth((tile_num-1)*tile_sz+1:tile_num*tile_sz)
872 
873 ! add Caspian Sea and Aral Sea to lake_frac and lake_depth
874  DO i = 1, tile_sz
875  IF (land_frac(i) < 0.9 .AND. lake_frac(i) < 0.1) THEN
876  IF (geolat(i) > 35.0 .AND. geolat(i) <= 50.0 .AND. &
877  geolon(i) > 45.0 .AND. geolon(i) <= 55.0) THEN
878  lake_frac(i) = 1.-land_frac(i)
879  lake_depth(i) = 211.0
880  ENDIF
881  IF (geolat(i) > 35.0 .AND. geolat(i) <= 50.0 .AND. &
882  geolon(i) > 57.0 .AND. geolon(i) <= 63.0) THEN
883  lake_frac(i) = 1.-land_frac(i)
884  lake_depth(i) = 10.0
885  ENDIF
886  ENDIF
887  ENDDO
888 
889 ! adjust land_frac and lake_frac, and make sure land_frac+lake_frac=1 at inland points
890  DO i = 1, tile_sz
891  if (lake_frac(i) >= lake_cutoff) then ! non-zero lake_frac dominates over land_frac
892  land_frac(i) = max(0., min(1., 1.-lake_frac(i)))
893  elseif (inland(i) == 1.) then ! land_frac dominates over lake_frac at inland points
894  lake_frac(i) = max(0., min(1., 1.-land_frac(i)))
895  end if
896 
897 ! epsil is "numerical" nonzero min for lake_frac/land_frac
898 ! lake_cutoff/land_cutoff is practical min for lake_frac/land_frac
899  if (min(lake_cutoff,land_cutoff) < epsil) then
900  print *,'lake_cutoff/land_cutoff cannot be smaller than epsil, reset...'
901  lake_cutoff=max(epsil,lake_cutoff)
902  land_cutoff=max(epsil,land_cutoff)
903  end if
904 
905  if (lake_frac(i)< lake_cutoff) lake_frac(i)=0.
906  if (lake_frac(i)>1.-lake_cutoff) lake_frac(i)=1.
907  if (inland(i) == 1.) land_frac(i) = 1.-lake_frac(i)
908  if (land_frac(i)< land_cutoff) land_frac(i)=0.
909  if (land_frac(i)>1.-land_cutoff) land_frac(i)=1.
910 
911  if (lake_frac(i) < lake_cutoff) then
912  lake_depth(i)=0.
913  elseif (lake_frac(i) >= lake_cutoff .and. lake_depth(i)==0.) then
914  lake_depth(i)=10.
915  end if
916  slmsk(i) = nint(land_frac(i))
917  ENDDO
918 
919 ! write 2 new variables
920  stat = nf90_put_var(ncid, lake_frac_id, lake_frac, &
921  start = (/ 1, 1 /), count = (/ tile_x_dim, tile_y_dim /) )
922  CALL nc_opchk(stat, "nf90_put_var: lake_frac")
923 
924  stat = nf90_put_var(ncid, lake_depth_id, lake_depth, &
925  start = (/ 1, 1 /), count = (/ tile_x_dim, tile_y_dim /) )
926  CALL nc_opchk(stat, "nf90_put_var: lake_depth")
927 
928 ! write back 2 modified variables: land_frac and slmsk
929  stat = nf90_put_var(ncid, land_frac_id, land_frac, &
930  start = (/ 1, 1 /), count = (/ tile_x_dim, tile_y_dim /) )
931  CALL nc_opchk(stat, "nf90_put_var: land_frac")
932 
933  stat = nf90_put_var(ncid, slmsk_id, slmsk, &
934  start = (/ 1, 1 /), count = (/ tile_x_dim, tile_y_dim /) )
935  CALL nc_opchk(stat, "nf90_put_var: slmsk")
936 
937  stat = nf90_close(ncid)
938  CALL nc_opchk(stat, "nf90_close")
939 
940  DEALLOCATE(lake_frac, lake_depth)
941  DEALLOCATE(geolon, geolat)
942  DEALLOCATE(land_frac, slmsk)
943 
944 END SUBROUTINE write_reg_lakedata_to_orodata
945 
951 SUBROUTINE nc_opchk(stat,opname)
952  USE netcdf
953  IMPLICIT NONE
954  INTEGER stat
955  CHARACTER(len=*) opname
956  CHARACTER(64) msg
957 
958  IF (stat .NE.0) THEN
959  msg = trim(opname) // ' Error, status code and message:'
960  print*,trim(msg), stat, nf90_strerror(stat)
961  stop
962  END IF
963 
964 END SUBROUTINE nc_opchk
965 
966 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:513
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:545
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:443
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:382
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:15
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
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:744
enclosure_cnvx
LOGICAL function enclosure_cnvx(v, n, p, co_gc)
Test whether a given point &#39;p&#39; is inside a convex spherical polygon defined with a series of &#39;n&#39; vert...
Definition: enclosure_cnvx.F90:49
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:144
nc_opchk
subroutine nc_opchk(stat, opname)
Check NetCDF return code and print error message.
Definition: inland.F90:437