cpld_gridgen 1.14.0
Loading...
Searching...
No Matches
topoedits.F90
Go to the documentation of this file.
1
8
9module topoedits
10
11 use gengrid_kinds, only: real_kind,int_kind
12 use grdvars, only: ni,nj
13 use grdvars, only: wet4,dp4,minimum_depth,maximum_depth,masking_depth
14 use charstrings, only: logmsg,history
15 use netcdf
16
17 implicit none
18 private
19
20 public add_topoedits
21 public apply_topoedits
22
23contains
31
32 subroutine add_topoedits(fsrc,fdst)
33
34 character(len=*), intent(in) :: fsrc, fdst
35
36 ! local variables
37 integer :: rc,id,i,j,ii,jj,ncid,dimid,idimid,dim1(1)
38 integer :: cnt1=0, cnt2=0, icnt
39 integer(int_kind), allocatable, dimension(:) :: ieds1, jeds1, ieds2, jeds2
40 real(real_kind), allocatable, dimension(:) :: zeds1, zeds2
41
42 !---------------------------------------------------------------------
43 ! read existing topo edits
44 !---------------------------------------------------------------------
45
46 rc = nf90_open(trim(fsrc), nf90_nowrite, ncid)
47 print '(a)','using topo edits file '//trim(fsrc)//' to edit land mask '
48
49 rc = nf90_inq_dimid(ncid, 'nEdits', dimid)
50 rc = nf90_inquire_dimension(ncid, dimid, len=cnt1)
51 rc = nf90_close(ncid)
52
53 ! return the existing values
54 allocate(ieds1(cnt1)); ieds1 = 0
55 allocate(jeds1(cnt1)); jeds1 = 0
56 allocate(zeds1(cnt1)); zeds1 = 0.0
57
58 rc = nf90_open(fsrc, nf90_nowrite, ncid)
59 rc = nf90_inq_varid(ncid, 'iEdit', id)
60 rc = nf90_get_var(ncid, id, ieds1)
61 rc = nf90_inq_varid(ncid, 'jEdit', id)
62 rc = nf90_get_var(ncid, id, jeds1)
63 rc = nf90_inq_varid(ncid, 'zEdit', id)
64 rc = nf90_get_var(ncid, id, zeds1)
65 rc = nf90_close(ncid)
66
67 !---------------------------------------------------------------------
68 ! determine the number of points to be added to topo-edits file
69 ! check only j=1 now
70 !---------------------------------------------------------------------
71
72 icnt = 0
73 j = 1
74 do i = 1,ni
75 if(wet4(i,j) .eq. 1.0)icnt = icnt+1
76 end do
77 cnt2 = cnt2+icnt
78 print '(a,i4,a,i4)', 'found ',icnt,' open water points at j=1 , cnt2 = ',cnt2
79
80 cnt2 = cnt1 + cnt2
81 ! allocate space for existing+new values and copy in original values
82 allocate(ieds2(cnt2)); ieds2 = 0
83 allocate(jeds2(cnt2)); jeds2 = 0
84 allocate(zeds2(cnt2)); zeds2 = 0.0
85
86 ieds2(1:cnt1) = ieds1(1:cnt1)
87 jeds2(1:cnt1) = jeds1(1:cnt1)
88 zeds2(1:cnt1) = zeds1(1:cnt1)
89
90 !---------------------------------------------------------------------
91 ! fill in new values and write new topoedits file
92 !---------------------------------------------------------------------
93
94 icnt = cnt1
95 j = 1
96 do i = 1,ni
97 if(wet4(i,j) .eq. 1.0)then
98 icnt = icnt+1
99 ii = i-1; jj = j-1
100 ieds2(icnt) = ii
101 jeds2(icnt) = jj
102 zeds2(icnt) = 0.0
103 end if
104 end do
105
106 !do i = 1,cnt2
107 ! print '(3i5,f12.4)',i,ieds2(i),jeds2(i),zeds2(i)
108 !end do
109
110 rc = nf90_create(fdst, nf90_write, ncid)
111 print '(a)', 'writing new topo edits to '//trim(fdst)
112
113 rc = nf90_def_dim(ncid, 'nEdits', cnt2, idimid)
114
115 rc = nf90_def_var(ncid, 'ni', nf90_int, id)
116 rc = nf90_def_var(ncid, 'nj', nf90_int, id)
117
118 dim1(:) = (/idimid/)
119 rc = nf90_def_var(ncid, 'iEdit', nf90_int, dim1, id)
120 rc = nf90_def_var(ncid, 'jEdit', nf90_int, dim1, id)
121 rc = nf90_def_var(ncid, 'zEdit', nf90_float, dim1, id)
122 rc = nf90_put_att(ncid, nf90_global, 'history', trim(history))
123 rc = nf90_enddef(ncid)
124
125 rc = nf90_inq_varid(ncid, 'ni', id)
126 rc = nf90_put_var(ncid, id, ni)
127 rc = nf90_inq_varid(ncid, 'nj', id)
128 rc = nf90_put_var(ncid, id, nj)
129
130 rc = nf90_inq_varid(ncid, 'iEdit', id)
131 rc = nf90_put_var(ncid, id, ieds2)
132 rc = nf90_inq_varid(ncid, 'jEdit', id)
133 rc = nf90_put_var(ncid, id, jeds2)
134 rc = nf90_inq_varid(ncid, 'zEdit', id)
135 rc = nf90_put_var(ncid, id, zeds2)
136
137 rc = nf90_close(ncid)
138
139 !---------------------------------------------------------------------
140 ! adjust land mask by same edits used at run time
141 !---------------------------------------------------------------------
142
143 do i = 1,cnt2
144 ii = ieds2(i); jj = jeds2(i)
145 if(wet4(ii+1,jj+1) .eq. 0.0 .and. zeds2(i) .gt. 0.0) then
146 wet4(ii+1,jj+1) = 1.0
147 print '(a,2i4,a)', 'switch point ',ii+1,jj+1,' from land->ocean at runtime'
148 end if
149 if(wet4(ii+1,jj+1) .eq. 1.0 .and. zeds2(i) .eq. 0.0) then
150 wet4(ii+1,jj+1) = 0.0
151 print '(a,2i4,a)', 'switch point ',ii+1,jj+1,' from ocean->land at runtime'
152 end if
153 end do
154 deallocate(ieds1, jeds1, zeds1)
155 deallocate(ieds2, jeds2, zeds2)
156
157 end subroutine add_topoedits
158
164
165 subroutine apply_topoedits(fsrc)
166
167 character(len=*), intent(in) :: fsrc
168
169 ! local variables
170 integer :: rc,ncid,id,dimid,i,j,ii,jj,cnt1
171 integer(int_kind), allocatable, dimension(:) :: ieds1, jeds1
172 real(real_kind), allocatable, dimension(:) :: zeds1
173
174 logical :: file_exists
175 !---------------------------------------------------------------------
176 ! read and apply topo edits file, if any
177 !---------------------------------------------------------------------
178
179 inquire(file=trim(fsrc),exist=file_exists)
180 if (file_exists) then
181 rc = nf90_open(trim(fsrc), nf90_nowrite, ncid)
182 print '(a)','using topo edits file '//trim(fsrc)//' to adjust bathymetry '
183
184 rc = nf90_inq_dimid(ncid, 'nEdits', dimid)
185 rc = nf90_inquire_dimension(ncid, dimid, len=cnt1)
186 rc = nf90_close(ncid)
187
188 ! return the existing values
189 allocate(ieds1(cnt1))
190 allocate(jeds1(cnt1))
191 allocate(zeds1(cnt1))
192
193 rc = nf90_open(fsrc, nf90_nowrite, ncid)
194 rc = nf90_inq_varid(ncid, 'iEdit', id)
195 rc = nf90_get_var(ncid, id, ieds1)
196 rc = nf90_inq_varid(ncid, 'jEdit', id)
197 rc = nf90_get_var(ncid, id, jeds1)
198 rc = nf90_inq_varid(ncid, 'zEdit', id)
199 rc = nf90_get_var(ncid, id, zeds1)
200 rc = nf90_close(ncid)
201
202 ! apply topo edits from file
203 do i = 1,cnt1
204 ii = ieds1(i); jj = jeds1(i)
205 print '(a,3i5,f8.2,a,f8.2)', 'Ocean topography edit: ', i, ii+1, jj+1 , dp4(ii+1,jj+1), '->', abs(zeds1(i))
206 dp4(ii+1,jj+1) = abs(zeds1(i))
207 end do
208 deallocate(ieds1, jeds1, zeds1)
209 end if
210
211 !---------------------------------------------------------------------
212 ! limit topography
213 !---------------------------------------------------------------------
214
215 print '(a)', 'Applying topo limits to ensure that min_depth < D(x,y) < max_depth '
216 print '(a, f8.2)', 'Using min_depth = ',minimum_depth
217 print '(a, f8.2)', 'Using max_depth = ',maximum_depth
218 do j = 1,nj
219 do i = 1,ni
220 if(dp4(i,j) > min(minimum_depth,masking_depth))then
221 dp4(i,j) = min( max(dp4(i,j), minimum_depth), maximum_depth)
222 end if
223 end do
224 end do
225
226 end subroutine apply_topoedits
227end module topoedits