ocean_merge  1.12.0
 All Files Functions
merge_lake_ocnmsk.f90
Go to the documentation of this file.
1 
7 
18  use netcdf
19 
20  implicit none
21 
22  character(len=120) :: pth1
23  character(len=120) :: pth2,pth3
24  character(len=10) :: atmres,ocnres
25  real, parameter :: min_land=1.e-4, def_lakedp=10.
26  ! this variable is now renamed as binary_lake and is passed in from the name
27  ! list
28  ! logical, parameter :: int_lake=.true.
29  ! all instances of int_lake was changed to binary_lake
30  integer :: binary_lake
31 
32  character(len=250) :: flnm
33  integer :: ncid,ndims,nvars,natts,lat,lon,v1id,v2id,v3id,v4id,start(2),count(2),i,j,latid,lonid,ncid4, dims(2),tile,nodp_pt
34  integer :: lake_pt,vlat
35  real, allocatable :: lake_frac(:,:),lake_depth(:,:),land_frac(:,:),ocn_frac(:,:),slmsk(:,:),lat2d(:,:)
36 
37 
38  call read_nml(pth1, pth2, atmres, ocnres, pth3,binary_lake)
39 
40  print *, pth1
41  nodp_pt=0
42  lake_pt=0
43  do tile=1,6
44  write(flnm,'(5a,i1,a)') trim(pth1),trim(atmres),'.',trim(ocnres),'.tile',tile,'.nc'
45  call handle_err(nf90_open(flnm, nf90_nowrite, ncid))
46  call handle_err(nf90_inquire(ncid, ndimensions=ndims, nvariables=nvars, nattributes=natts))
47  write(6,*) 'flnm_ocn=',flnm,' ncid=',ncid, ' ndims=',ndims, 'nvars=',nvars,' natts=',natts
48  call handle_err(nf90_inq_dimid(ncid, 'grid_xt', latid)) ! RM: lon is no longer in this file; try grid_xt
49  call handle_err(nf90_inq_dimid(ncid, 'grid_yt', lonid)) ! RM: lat is no longer in this file; try grid_tyt
50  call handle_err(nf90_inquire_dimension(ncid, latid, len=lat))
51  call handle_err(nf90_inquire_dimension(ncid, lonid, len=lon))
52  if (tile==1) then
53  write(6,*) 'lat=',lat,'lon=',lon
54  allocate (lake_frac(lon,lat),lake_depth(lon,lat),land_frac(lon,lat),ocn_frac(lon,lat),slmsk(lon,lat),lat2d(lon,lat))
55  start(1:2) = (/1,1/)
56  count(1:2) = (/lon,lat/)
57  end if
58  call handle_err(nf90_inq_varid(ncid, 'land_frac', v1id))
59  call handle_err(nf90_get_var(ncid, v1id, ocn_frac, start=start, count=count))
60 
61  !write(flnm,'(3a,i1,a)') trim(pth2),trim(atmres),'_orog_data.tile',tile,'.nc'
62  write(flnm,'(4a,i1,a)') trim(pth2),'oro.',trim(atmres),'.tile',tile,'.nc'
63  print *,' flnm2=',trim(flnm)
64  call handle_err(nf90_open(flnm, nf90_nowrite, ncid))
65  call handle_err(nf90_inquire(ncid, ndimensions=ndims, nvariables=nvars, nattributes=natts))
66  write(6,*) 'flnm_lake=',flnm,' ncid=',ncid, ' ndims=',ndims, 'nvars=',nvars,' natts=',natts
67  !if (tile==1) then
68  ! call handle_err (nf90_inq_dimid (ncid, 'lat', latid))
69  ! call handle_err (nf90_inq_dimid (ncid, 'lon', lonid))
70  ! call handle_err (nf90_inquire_dimension (ncid, latid, len=lat))
71  ! call handle_err (nf90_inquire_dimension (ncid, lonid, len=lon))
72  ! write(6,*) 'lat=',lat,'lon=',lon
73  ! start(1:2) = (/1,1/)
74  ! count(1:2) = (/lon,lat/)
75  !end if
76  call handle_err(nf90_inq_varid(ncid, 'lake_frac', v2id))
77  call handle_err(nf90_inq_varid(ncid, 'lake_depth',v3id))
78  call handle_err(nf90_inq_varid(ncid, 'geolat' ,vlat))
79  call handle_err(nf90_get_var(ncid, v2id, lake_frac, start=start, count=count))
80  call handle_err(nf90_get_var(ncid, v3id, lake_depth,start=start, count=count))
81  call handle_err(nf90_get_var(ncid, vlat, lat2d, start=start, count=count))
82 
83  do i=1,lon
84  do j=1,lat
85  if (binary_lake.eq.1) lake_frac(i,j)=nint(lake_frac(i,j)) ! using integer lake_frac
86  if (lat2d(i,j).le.-60.) lake_frac(i,j)=0. ! ignore lakes on Antarctica
87  land_frac(i,j)=1.-ocn_frac(i,j)
88  if (land_frac(i,j) < min_land) land_frac(i,j)=0. ! ignore land < min_land
89  if (land_frac(i,j) > 1.-min_land) land_frac(i,j)=1. ! ignore water < min_land
90  if (1.-land_frac(i,j) > 0.) lake_frac(i,j)=0. ! ocn dominates
91 
92  if (lake_frac(i,j) > 0.) then
93  lake_pt=lake_pt+1 ! calculating total lake points
94  if (binary_lake.eq.1) then
95  land_frac(i,j)=0.
96  else
97  land_frac(i,j)=1.-lake_frac(i,j)
98  end if
99  if (lake_depth(i,j) <= 0.) then
100  lake_depth(i,j)=def_lakedp ! set missing lake depth to default value
101  nodp_pt=nodp_pt+1 ! calculating total lake points without depth
102  end if
103  else
104  lake_depth(i,j)=0.
105  end if
106 ! slmsk(i,j) = ceiling(land_frac(i,j))! "ceiling is used for orog smoothing"
107  slmsk(i,j) = nint(land_frac(i,j)) ! nint got the land pts correct
108  end do
109  end do
110 
111  write(flnm,'(4a,i1,a)') trim(atmres),'.',trim(ocnres),'.tile',tile,'.nc'
112  print *,'output=',trim(flnm)
113  call handle_err(nf90_create(path=trim(pth3)//trim(flnm), &
114  cmode=or(nf90_clobber, nf90_64bit_offset), ncid=ncid4)) ! netcdf3
115 
116  call handle_err(nf90_def_dim(ncid4,'lon', lon, dims(1)))
117  call handle_err(nf90_def_dim(ncid4,'lat', lat, dims(2)))
118  call handle_err(nf90_def_var(ncid4,'land_frac', nf90_float, dims(1:2), v1id))
119  call handle_err(nf90_def_var(ncid4,'lake_frac', nf90_float, dims(1:2), v2id))
120  call handle_err(nf90_def_var(ncid4,'lake_depth',nf90_float, dims(1:2), v3id))
121  call handle_err(nf90_def_var(ncid4,'slmsk', nf90_float, dims(1:2), v4id))
122 
123  call handle_err(nf90_enddef(ncid4))
124 
125  call handle_err(nf90_put_var(ncid4, v1id,land_frac))
126  call handle_err(nf90_put_var(ncid4, v2id,lake_frac))
127  call handle_err(nf90_put_var(ncid4, v3id,lake_depth))
128  call handle_err(nf90_put_var(ncid4, v4id,slmsk))
129  call handle_err(nf90_close(ncid4))
130 
131 ! do i=1,48
132 ! do j=1,48
133 ! write(*,'(2i4,4f6.1)') i,j,land_frac(i,j),lake_frac(i,j),lake_depth(i,j),slmsk(i,j)
134 ! end do
135 ! end do
136 
137  end do ! tile
138  write(*,'(a,i8,a,i8,a)') 'total lake point ',lake_pt,' where ',nodp_pt,' has no depth'
139 
140 end program merge_lake_ocnmsk
141 
146 subroutine handle_err (ret)
147  use netcdf
148  integer, intent(in) :: ret
149 
150  if (ret /= nf90_noerr) then
151  write(6,*) nf90_strerror(ret)
152  stop 999
153  end if
154 end subroutine handle_err
155 
166 subroutine read_nml(ocean_mask_dir, lake_mask_dir, atmres,ocnres,out_dir,binary_lake)
167 
168  integer :: unit=7, io_status
169 
170  character(len=120), intent(out) :: ocean_mask_dir
171  character(len=120), intent(out) :: lake_mask_dir
172  character(len=120), intent(out) :: out_dir
173  character(len=10), intent(out) :: atmres,ocnres
174  integer, intent(out):: binary_lake
175 
176  namelist/mask_nml/ocean_mask_dir, lake_mask_dir, atmres, ocnres,out_dir,binary_lake
177  open(unit=unit, file='input.nml', iostat=io_status )
178  read(unit,mask_nml, iostat=io_status )
179  close(unit)
180  if (io_status > 0) then
181  print *,'Error reading input.nml'
182  call handle_err(-1)
183  end if
184 end subroutine read_nml
subroutine handle_err(ret)
Handle netCDF errors.
program merge_lake_ocnmsk
Determine the water mask by merging the lake mask with the mapped ocean mask from MOM6...
subroutine read_nml(ocean_mask_dir, lake_mask_dir, atmres, ocnres, out_dir, binary_lake)
Read program namelist.