1module oceanMask_mod
2 ! MODULE oceanMask_mod (prefix='ocm' category='7. Low-level data objects')
3 !
4 !:Purpose: The horizontal mask indicating land (.false.) and water (.true.) grid points.
5 ! This mask is either:
6 ! * In the case of variables on ocean depth levels, it varies with vertical level.
7 ! * In other cases it is a single 2D field used for all variables.
8 !
9 use midasMpi_mod
10 use kdTree2_mod
11 use horizontalCoord_mod
12 use verticalCoord_mod
13 use utilities_mod
14 implicit none
15 save
16 private
17
18 ! public structure definition
19 public :: struct_ocm
20
21 ! public subroutines and functions
22 public :: ocm_readMaskFromFile, ocm_deallocate
23 public :: ocm_copyMask, ocm_communicateMask
24 public :: ocm_farFromLand
25 public :: ocm_copyToInt, ocm_copyFromInt
26
27 type struct_ocm
28 ! This is the derived type of the ocean mask object
29 integer :: nLev
30 logical, pointer :: mask(:,:,:) => null()
31 logical :: maskPresent = .false.
32 type(struct_hco), pointer :: hco
33 end type struct_ocm
34
35 integer, external :: get_max_rss
36
37 contains
38
39 !--------------------------------------------------------------------------
40 ! ocm_readMaskFromFile
41 !--------------------------------------------------------------------------
42 subroutine ocm_readMaskFromFile(oceanMask, hco, vco, filename)
43 !
44 !:Purpose: Check if any mask fields exist for surface or ocean depth levels.
45 !
46 !
47 implicit none
48
49 ! Arguments:
50 type(struct_ocm), intent(inout) :: oceanMask
51 type(struct_hco), pointer, intent(in) :: hco
52 type(struct_vco), intent(in) :: vco
53 character(len=*), intent(in) :: fileName
54
55 ! Locals:
56 integer :: nulfile, ierr, ip1, ni_file, nj_file, nk_file
57 integer :: ikey, levIndex
58 integer :: fnom, fstouv, fclos, fstfrm, fstluk, fstinf, fstsui
59 integer, allocatable :: mask(:,:)
60 integer :: maxkeys
61
62 ! Check if any mask is present in file, return if not
63 if ( .not. utl_varNamePresentInFile(' ',fileName_opt=trim(fileName),typvar_opt='@@') ) then
64 return
65 end if
66
67 !- Open input file
68 nulfile = 0
69 write(*,*) 'ocm_readMaskFromFile: file name = ',trim(fileName)
70 ierr = fnom(nulfile,trim(fileName),'RND+OLD+R/O',0)
71
72 if ( ierr >= 0 ) then
73 maxkeys = fstouv(nulfile,'RND+OLD')
74 else
75 call utl_abort('ocm_readMaskFromFile: problem opening input file')
76 end if
77
78 if (nulfile == 0 ) then
79 call utl_abort('ocm_readMaskFromFile: unit number for input file not valid')
80 end if
81
82 ! Read mask for all fields
83 if ( vco%nLev_depth > 0 ) then
84 lev_loop: do levIndex = 1, vco%nLev_depth
85
86 ip1 = vco%ip1_depth(levIndex)
87
88 ! Make sure that the mask for this variable has the same grid size as hco
89 ikey = fstinf(nulfile, ni_file, nj_file, nk_file, &
90 -1, ' ', ip1, -1, -1, '@@', ' ')
91
92 if (ikey < 0) then
93 write(*,*) 'ocm_readMaskFromFile: searched for mask with ip1 = ', ip1
94 call utl_abort('ocm_readMaskFromFile: cannot find mask for this ip1 in file ' // trim(fileName))
95 end if
96
97 do while (ni_file /= hco%ni .or. nj_file /= hco%nj)
98 ikey = fstsui(nulfile, ni_file, nj_file, nk_file)
99 end do
100
101 if (ni_file == hco%ni .and. nj_file == hco%nj) then
102 write(*,*) 'ocm_readMaskFromFile: read mask for ip1 = ', ip1
103 if (.not. associated(oceanMask%mask)) then
104 call ocm_allocate(oceanMask,hco,vco%nLev_depth)
105 end if
106 if (.not. allocated(mask)) allocate(mask(hco%ni,hco%nj))
107 ierr = fstluk(mask(:,:), ikey, ni_file, nj_file, nk_file)
108 call ocm_copyFromInt(oceanMask,mask,levIndex)
109 if (ierr < 0) then
110 call utl_abort('ocm_readMaskFromFile: error when reading mask record')
111 end if
112 else
113 ! Special cases for variables that are on a different horizontal grid in LAM (e.g. TG)
114 write(*,*)
115 write(*,*) 'ocm_readMaskFromFile: mask is on a different horizontal grid'
116 write(*,*) 'ni = ', ni_file, hco%ni, ', nj = ', nj_file, hco%nj
117 call utl_abort('ocm_readMaskFromFile: This is not allowed at the moment')
118 end if
119 end do lev_loop
120
121 else
122 ! ***No depth levels, so just look for any mask field***
123
124 ! Make sure that the mask for this variable has the same grid size as hco
125 ikey = fstinf(nulfile, ni_file, nj_file, nk_file, &
126 -1, ' ', -1, -1, -1, '@@', ' ')
127
128 if (ikey < 0) then
129 call utl_abort('ocm_readMaskFromFile: cannot find any mask in file ' // trim(fileName))
130 end if
131
132 do while (ni_file /= hco%ni .or. nj_file /= hco%nj)
133 ikey = fstsui(nulfile, ni_file, nj_file, nk_file)
134 end do
135
136 if (ni_file == hco%ni .and. nj_file == hco%nj) then
137 write(*,*) 'ocm_readMaskFromFile: reading mask'
138 if (.not. associated(oceanMask%mask)) then
139 call ocm_allocate(oceanMask,hco,1)
140 if (.not. allocated(mask)) allocate(mask(hco%ni,hco%nj))
141 ierr = fstluk(mask(:,:), ikey, ni_file, nj_file, nk_file)
142 call ocm_copyFromInt(oceanMask,mask,1)
143 if (ierr < 0) then
144 call utl_abort('ocm_readMaskFromFile: error when reading mask record')
145 end if
146 end if
147 else
148 ! Special cases for variables that are on a different horizontal grid in LAM (e.g. TG)
149 write(*,*)
150 write(*,*) 'ocm_readMaskFromFile: mask is on a different horizontal grid'
151 write(*,*) 'ni = ', ni_file, hco%ni, ', nj = ', nj_file, hco%nj
152 call utl_abort('ocm_readMaskFromFile: This is not allowed at the moment')
153 end if
154
155 end if
156
157 ierr = fstfrm(nulfile)
158 ierr = fclos(nulfile)
159
160 if (allocated(mask)) deallocate(mask)
161
162 end subroutine ocm_readMaskFromFile
163
164 !--------------------------------------------------------------------------
165 ! ocm_farFromLand
166 !--------------------------------------------------------------------------
167 function ocm_farFromLand(oceanMask, levIndex, lon, lat, distanceToLand) result(farFromLand)
168 !
169 !:Purpose: Determine if the supplied lat/lat location is far from the
170 ! nearest land based on the supplied 'distanceToLand'.
171 !
172 implicit none
173
174 ! Arguments:
175 type(struct_ocm), intent(in) :: oceanMask
176 integer, intent(in) :: levIndex
177 real(8), intent(in) :: lon
178 real(8), intent(in) :: lat
179 real(8), intent(in) :: distanceToLand
180 ! Result:
181 logical :: farFromLand
182
183 ! Locals:
184 integer, parameter :: maxNumLocalGridPointsSearch = 200000
185 type(kdtree2), save, pointer :: tree => null()
186 integer :: ni, nj, xIndex, yIndex, gridIndex
187 integer :: numTotalLandPoints, numLocalGridPointsFound
188 real(kdkind), allocatable :: positionArray(:,:)
189 type(kdtree2_result) :: searchResults(maxNumLocalGridPointsSearch)
190 real(kdkind) :: searchRadiusSquared
191 real(kdkind) :: refPosition(3)
192
193 ! do some basic checks
194 if (.not.oceanMask%maskPresent .or. .not.associated(oceanMask%mask)) then
195 call utl_abort('ocm_farFromLand: mask is not allocated')
196 end if
197 if (levIndex < 0 .or. levIndex > oceanMask%nLev) then
198 call utl_abort('ocm_farFromLand: specified levIndex not valid')
199 end if
200
201 ni = oceanMask%hco%ni
202 nj = oceanMask%hco%nj
203 searchRadiusSquared = (1.1D0 * distanceToLand)**2
204
205 ! create the kdtree on the first call
206 if (.not. associated(tree)) then
207 write(*,*) 'ocm_farFromLand: start creating kdtree'
208 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
209
210 numTotalLandPoints = count(.not. oceanMask%mask(:,:,levIndex))
211 allocate(positionArray(3,numTotalLandPoints))
212
213 gridIndex = 0
214 do xIndex = 1, ni
215 do yIndex = 1, nj
216 if (.not. oceanMask%mask(xIndex,yIndex,levIndex)) then
217 gridIndex = gridIndex + 1
218 positionArray(:,gridIndex) = &
219 kdtree2_3dPosition(real(oceanMask%hco%lon2d_4(xIndex,yIndex),8), &
220 real(oceanMask%hco%lat2d_4(xIndex,yIndex),8))
221 end if
222 end do
223 end do
224 tree => kdtree2_create(positionArray, sort=.true., rearrange=.true.)
225 write(*,*) 'ocm_farFromLand: done creating kdtree'
226 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
227
228 end if
229
230 ! do the search
231 refPosition(:) = kdtree2_3dPosition(lon, lat)
232
233 call kdtree2_r_nearest(tp=tree, qv=refPosition, r2=searchRadiusSquared, &
234 nfound=numLocalGridPointsFound, &
235 nalloc=maxNumLocalGridPointsSearch, results=searchResults)
236 if (numLocalGridPointsFound > maxNumLocalGridPointsSearch) then
237 call utl_abort('ocm_farFromLand: the parameter maxNumLocalGridPointsSearch must be increased')
238 end if
239 if (numLocalGridPointsFound == 0) then
240 farFromLand = .true.
241 else
242 farFromLand = ( sqrt(searchResults(1)%dis) > distanceToLand )
243 end if
244
245 end function ocm_farFromLand
246
247 !--------------------------------------------------------------------------
248 ! ocm_copyMask
249 !--------------------------------------------------------------------------
250 subroutine ocm_copyMask(oceanMask_in, oceanMask_out, beSilent_opt)
251 !
252 !:Purpose: Copy the mask data from one instance of oceanMask to
253 ! another. If the destination instance is not already
254 ! allocated, then this will also be done.
255 !
256 implicit none
257
258 ! Arguments:
259 type(struct_ocm), intent(in) :: oceanMask_in
260 logical, optional, intent(in) :: beSilent_opt
261 type(struct_ocm), intent(inout) :: oceanMask_out
262
263 ! Locals:
264 logical :: beSilent
265
266 if ( present(beSilent_opt) ) then
267 beSilent = beSilent_opt
268 else
269 beSilent = .false.
270 end if
271
272 if (.not.oceanMask_in%maskPresent .or. .not.associated(oceanMask_in%mask)) then
273 if ( .not. beSilent ) write(*,*) 'ocm_copyMask: no input mask, do nothing'
274 return
275 end if
276
277 if (.not.associated(oceanMask_out%mask)) then
278 call ocm_allocate(oceanMask_out, oceanMask_in%hco, oceanMask_in%nLev)
279 end if
280
281 if ( .not. beSilent ) write(*,*) 'ocm_copyMask: copying over the horizontal mask'
282 oceanMask_out%mask(:,:,:) = oceanMask_in%mask(:,:,:)
283 oceanMask_out%maskPresent = .true.
284
285 end subroutine ocm_copyMask
286
287 !--------------------------------------------------------------------------
288 ! ocm_communicateMask
289 !--------------------------------------------------------------------------
290 subroutine ocm_communicateMask(oceanMask)
291 !
292 !:Purpose: Copy mask fields from task 0 to all others
293 !
294 implicit none
295
296 ! Arguments:
297 type(struct_ocm), intent(inout) :: oceanMask
298
299 ! Locals:
300 integer :: ierr
301 type(struct_hco), pointer :: hco_temp
302
303 write(*,*) 'ocm_communicateMask: starting'
304
305 call rpn_comm_bcast(oceanMask%maskPresent, 1, &
306 'MPI_LOGICAL', 0, 'GRID', ierr)
307 if (.not.oceanMask%maskPresent) then
308 write(*,*) 'ocm_communicateMask: mask not present, return'
309 return
310 end if
311
312 call rpn_comm_bcast(oceanMask%nLev, 1, &
313 'MPI_INTEGER', 0, 'GRID', ierr)
314
315 ! special treatment of hco object since EZscintID not properly communicated
316 nullify(hco_temp)
317 if (mmpi_myid > 0 .and. associated(oceanMask%hco)) then
318 hco_temp => oceanMask%hco
319 nullify(oceanMask%hco)
320 end if
321 call hco_mpiBcast(oceanMask%hco)
322
323 if (associated(hco_temp)) then
324 call hco_deallocate(oceanMask%hco)
325 oceanMask%hco => hco_temp
326 end if
327
328 if (.not.associated(oceanMask%mask)) then
329 call ocm_allocate(oceanMask,oceanMask%hco,oceanMask%nLev)
330 end if
331 call rpn_comm_bcast(oceanMask%mask, oceanMask%hco%ni*oceanMask%hco%nj*1, &
332 'MPI_LOGICAL', 0, 'GRID', ierr)
333
334 write(*,*) 'ocm_communicateMask: finished'
335
336 end subroutine ocm_communicateMask
337
338 !--------------------------------------------------------------------------
339 ! ocm_allocate
340 !--------------------------------------------------------------------------
341 subroutine ocm_allocate(oceanMask,hco,nLev)
342 !
343 !:Purpose: Allocate the object, if it isn't already.
344 !
345 implicit none
346
347 ! Arguments:
348 type(struct_ocm), intent(inout) :: oceanMask
349 type(struct_hco), pointer, intent(in) :: hco
350 integer, intent(in) :: nLev
351
352 if (.not.associated(oceanMask%mask)) then
353 allocate(oceanMask%mask(hco%ni,hco%nj,nLev))
354 oceanMask%maskPresent = .true.
355 oceanMask%hco => hco
356 oceanMask%nLev = nLev
357 end if
358
359 end subroutine ocm_allocate
360
361 !--------------------------------------------------------------------------
362 ! ocm_deallocate
363 !--------------------------------------------------------------------------
364 subroutine ocm_deallocate(oceanMask)
365 !
366 !:Purpose: Deallocate object.
367 !
368 implicit none
369
370 ! Arguments:
371 type(struct_ocm), intent(inout) :: oceanMask
372
373 if (associated(oceanMask%mask)) then
374 deallocate(oceanMask%mask)
375 nullify(oceanMask%mask)
376 nullify(oceanMask%hco)
377 oceanMask%maskPresent = .false.
378 oceanMask%nLev = 0
379 end if
380
381 end subroutine ocm_deallocate
382
383 !--------------------------------------------------------------------------
384 ! ocm_logicalToInt
385 !--------------------------------------------------------------------------
386 subroutine ocm_copyToInt(oceanMask, intArray, maskLev)
387 !
388 !:Purpose: Convert the selected level of the logical oceanMask
389 ! object into integer values where true is 1 and false is 0.
390 !
391 implicit none
392
393 ! Arguments:
394 type(struct_ocm), intent(inout) :: oceanMask
395 integer, intent(out) :: intArray(:,:)
396 integer, intent(in) :: maskLev
397
398 intArray(:,:) = 0
399 where(oceanMask%mask(:,:,maskLev)) intArray(:,:) = 1
400
401 end subroutine ocm_copyToInt
402
403 !--------------------------------------------------------------------------
404 ! ocm_logicalToInt
405 !--------------------------------------------------------------------------
406 subroutine ocm_copyFromInt(oceanMask, intArray, maskLev)
407 !
408 !:Purpose: Convert a 2D integer array into logical values
409 ! where true is 1 and false is 0 and copy into
410 ! the selected level of the oceanMask object.
411 !
412 implicit none
413
414 ! Arguments:
415 type(struct_ocm), intent(inout) :: oceanMask
416 integer, intent(in) :: intArray(:,:)
417 integer, intent(in) :: maskLev
418
419 oceanMask%mask(:,:,maskLev) = .false.
420 where(intArray(:,:) == 1) oceanMask%mask(:,:,maskLev) = .true.
421
422 end subroutine ocm_copyFromInt
423
424end module oceanMask_mod
