New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
domutl.F90 in branches/UKMO/dev_isf_remapping_UKESM_GO6package_r9314/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/UKMO/dev_isf_remapping_UKESM_GO6package_r9314/NEMOGCM/NEMO/OPA_SRC/DOM/domutl.F90 @ 9513

Last change on this file since 9513 was 9513, checked in by mathiot, 6 years ago

Add option to detect and remove subglacial lake (do not affect closed sea option)
File size: 11.1 KB
Line 
1MODULE domutl
2   !!
3   !!======================================================================
4   !!                       ***  MODULE  closea  ***
5   !! modutl : specific module to do
6   !!              - 2d flood filling (closea)
7   !!              - 3d flood filling (zgr_isf)
8   !!
9   !!======================================================================
10   !! History :   4.0  !  03-18  (P. Mathiot) Original code
11   !!----------------------------------------------------------------------
12   !!
13   !!----------------------------------------------------------------------
14   !!   fillpool    : replace all point connected to the seed by a fillvalue
15   !!----------------------------------------------------------------------
16   !!
17   USE dom_oce         ! ocean space and time domain
18   USE domngb
19   USE in_out_manager    ! I/O manager
20   USE lib_mpp
21   USE lbclnk
22
23   IMPLICIT NONE
24   PRIVATE
25
26   PUBLIC fill_pool      ! routine called by domain module
27
28   INTERFACE fill_pool
29      MODULE PROCEDURE FillPool2D, FillPool3D
30   END INTERFACE
31
32   TYPE idx
33      INTEGER :: i,j,k
34      TYPE(idx), POINTER :: next
35   END TYPE idx
36
37CONTAINS
38   SUBROUTINE FillPool3D(kiseed, kjseed, kkseed, rdta, rfill)
39      !!---------------------------------------------------------------------
40      !!                  ***  ROUTINE FillPool3D  ***
41      !!
42      !! ** Purpose :  Replace all area surrounding by mask value by kifill value
43      !!
44      !! ** Method  :  flood fill algorithm
45      !!
46      !!----------------------------------------------------------------------
47      INTEGER(KIND=4),                  INTENT(in)    :: kiseed, kjseed, kkseed
48      REAL(wp),                         INTENT(in)    :: rfill    ! filling value
49      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: rdta     ! input data
50      REAL(wp), DIMENSION(jpi,jpj,jpk)                :: rseedmap, rseedmap_b  !
51 
52      INTEGER :: ip=0                 ! size of the pile
53      INTEGER :: ii  , ij  , ik   , kii, kjj, jj, kk    ! working integer
54      INTEGER :: iip1, ijp1, ikp1
55      INTEGER :: iim1, ijm1, ikm1
56      INTEGER :: nseed
57      TYPE (idx), POINTER :: seed
58      !!----------------------------------------------------------------------
59      !
60      ! initialisation seed
61      NULLIFY(seed)
62      !
63      ! create the first seed and update nseed for all processors
64      nseed = 0
65      ii=kiseed; ij=kjseed ; ik=kkseed
66      IF ((mi0(ii) == mi1(ii)) .AND. (mj0(ij) == mj1(ij))) THEN   ! T if seed on the local domain
67         ii = mi0(ii) ; ij = mj0(ij);                             ! conversion to local index
68         IF (rdta(ii,ij,ik) > 0 ) THEN                            ! create seed if not on land
69            CALL add_head_idx(seed,ii,ij,ik)
70            rdta    (ii,ij,ik) = rfill
71            rseedmap(ii,ij,ik) = 1.
72         END IF
73      END IF
74      nseed=SUM(rseedmap); IF( lk_mpp )   CALL mpp_sum( nseed )  ! nseed =0 means on land => WARNING later on
75      !
76      ! loop until the pile size is 0 or if the pool is larger than the critical size
77      IF (nseed > 0) THEN
78         ! seed on ocean continue
79         DO WHILE ( nseed /= 0 )
80            DO WHILE ( ASSOCIATED(seed) )
81               ip=ip+1
82               ii=seed%i; ij=seed%j ; ik=seed%k ; rseedmap(ii,ij,ik)=1.
83               !
84               ! update bathy and update pile size
85               CALL del_idx(seed) 
86               !
87               ! check neighbour cells
88               iip1=ii+1 ; IF (iip1 == jpi+1 ) iip1=ii
89               iim1=ii-1 ; IF (iim1 == 0     ) iim1=ii
90               ijp1=ij+1 ; IF (ijp1 == jpj+1 ) ijp1=ij
91               ijm1=ij-1 ; IF (ijm1 == 0     ) ijm1=ij
92               ikp1=ik+1 ; IF (ikp1 == jpk+1 ) ikp1=ik
93               ikm1=ik-1 ; IF (ikm1 == 0     ) ikm1=ik
94               !
95               ! ji direction
96               IF (rdta(ii,ijp1,ik) > 0 .AND. rdta(ii,ijp1,ik) /= rfill ) THEN ; CALL add_head_idx(seed,ii,ijp1,ik) ; rdta(ii,ijp1,ik) = rfill ; rseedmap(ii,ijp1,ik)=1. ; END IF
97               IF (rdta(ii,ijm1,ik) > 0 .AND. rdta(ii,ijm1,ik) /= rfill ) THEN ; CALL add_head_idx(seed,ii,ijm1,ik) ; rdta(ii,ijm1,ik) = rfill ; rseedmap(ii,ijm1,ik)=1. ; END IF
98               !
99               ! jj direction
100               IF (rdta(iip1,ij,ik) > 0 .AND. rdta(iip1,ij,ik) /= rfill ) THEN ; CALL add_head_idx(seed,iip1,ij,ik) ; rdta(iip1,ij,ik) = rfill ; rseedmap(iip1,ij,ik)=1. ; END IF
101               IF (rdta(iim1,ij,ik) > 0 .AND. rdta(iim1,ij,ik) /= rfill ) THEN ; CALL add_head_idx(seed,iim1,ij,ik) ; rdta(iim1,ij,ik) = rfill ; rseedmap(iim1,ij,ik)=1. ; END IF
102               !
103               ! jk direction
104               IF (rdta(ii,ij,ikp1) > 0 .AND. rdta(ii,ij,ikp1) /= rfill ) THEN ; CALL add_head_idx(seed,ii,ij,ikp1) ; rdta(ii,ij,ikp1) = rfill ; rseedmap(ii,ij,ik)=1. ; END IF
105               IF (rdta(ii,ij,ikm1) > 0 .AND. rdta(ii,ij,ikm1) /= rfill ) THEN ; CALL add_head_idx(seed,ii,ij,ikm1) ; rdta(ii,ij,ikm1) = rfill ; rseedmap(ii,ij,ik)=1. ; END IF
106            END DO
107            !
108            ! exchange seed
109            nseed=SUM(rseedmap); IF( lk_mpp )   CALL mpp_sum( nseed )  ! this is the sum of all the point check this iteration
110            !
111            rseedmap_b(:,:,:)=rseedmap(:,:,:)
112            CALL lbc_lnk(rseedmap,'T',1.)
113            !
114            ! build new list of seed
115            DO ii=1,jpi
116               DO jj=1,jpj
117                  DO kk=1,jpk
118                     IF (rseedmap(ii,jj,kk) > 0.0 .AND. rseedmap(ii,jj,kk) /= rseedmap_b(ii,jj,kk) .AND. rdta(ii,jj,kk) > 0 .AND. rdta(ii,jj,kk) /= rfill) THEN
119                        CALL add_head_idx(seed,ii,jj,kk)
120                     END IF
121                  END DO
122               END DO
123            END DO
124            !
125            ! reset map of seed
126            rseedmap(:,:,:)=0.0
127            !
128         END DO
129      ELSE
130         WRITE(numout,*) 'W A R N I N G : SEED (',ii,ij,') is on land, nothing to do'
131      END IF
132 
133   END SUBROUTINE FillPool3D
134   SUBROUTINE FillPool2D(kiseed, kjseed, rdta, rfill)
135      !!---------------------------------------------------------------------
136      !!                  ***  ROUTINE FillPool3D  ***
137      !!
138      !! ** Purpose :  Replace all area surrounding by mask value by kifill value
139      !!
140      !! ** Method  :  flood fill algorithm
141      !!
142      !!----------------------------------------------------------------------
143      INTEGER(KIND=4),            INTENT(in)    :: kiseed, kjseed
144      REAL(wp),                   INTENT(in)    :: rfill    ! filling value
145      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: rdta     ! input data
146      REAL(wp), DIMENSION(jpi,jpj) :: rseedmap, rseedmap_b 
147 
148      INTEGER :: ip=0                     ! size of the pile
149      INTEGER :: ii  , ij  , jj, kii, kjj     ! working integer
150      INTEGER :: iip1, ijp1               ! working integer
151      INTEGER :: iim1, ijm1
152      INTEGER :: nseed
153      TYPE (idx), POINTER :: seed
154      !!----------------------------------------------------------------------
155      !
156      ! initialisation seed
157      NULLIFY(seed)
158      !
159      ! create the first seed and update nseed for all processors
160      nseed = 0
161      ii = kiseed ; ij = kjseed ;
162      IF ((mi0(ii) == mi1(ii)) .AND. (mj0(ij) == mj1(ij))) THEN   ! T if seed on the local domain
163         ii = mi0(ii) ; ij = mj0(ij)                              ! conversion to local index
164         IF (rdta(ii,ij) > 0 ) THEN                               ! create seed if not on land
165            CALL add_head_idx(seed,ii,ij,1)
166            rdta    (ii,ij) = rfill
167            rseedmap(ii,ij) = 1.
168         END IF
169      END IF
170      nseed=SUM(rseedmap); IF( lk_mpp )   CALL mpp_sum( nseed )  ! nseed =0 means on land => WARNING later on
171      !
172      ! loop until the pile size is 0 or if the pool is larger than the critical size
173      IF (nseed > 0) THEN
174         ! seed on ocean continue
175         DO WHILE ( nseed .NE. 0 )
176            DO WHILE ( ASSOCIATED(seed) )
177               ip=ip+1
178               ii=seed%i; ij=seed%j ; rseedmap(ii,ij)=1.
179               ! update bathy and update pile size
180               CALL del_idx(seed) 
181               !
182               ! check neighbour cells
183               iip1=ii+1 ; IF (iip1 == jpi+1 ) iip1=ii
184               iim1=ii-1 ; IF (iim1 == 0     ) iim1=ii
185               ijp1=ij+1 ; IF (ijp1 == jpj+1 ) ijp1=ij
186               ijm1=ij-1 ; IF (ijm1 == 0     ) ijm1=ij
187               !
188               ! ji direction
189               IF (rdta(ii  ,ijp1) > 0 .AND. rdta(ii  ,ijp1) /= rfill ) THEN ; CALL add_head_idx(seed,ii  ,ijp1,1) ; rdta(ii  ,ijp1) = rfill ; rseedmap(ii  ,ijp1)=1. ; END IF
190               IF (rdta(ii  ,ijm1) > 0 .AND. rdta(ii  ,ijm1) /= rfill ) THEN ; CALL add_head_idx(seed,ii  ,ijm1,1) ; rdta(ii  ,ijm1) = rfill ; rseedmap(ii  ,ijm1)=1. ; END IF
191               !
192               ! jj direction
193               IF (rdta(iip1,ij  ) > 0 .AND. rdta(iip1,ij  ) /= rfill ) THEN ; CALL add_head_idx(seed,iip1,ij  ,1) ; rdta(iip1,ij  ) = rfill ; rseedmap(iip1,ij  )=1. ; END IF
194               IF (rdta(iim1,ij  ) > 0 .AND. rdta(iim1,ij  ) /= rfill ) THEN ; CALL add_head_idx(seed,iim1,ij  ,1) ; rdta(iim1,ij  ) = rfill ; rseedmap(iim1,ij  )=1. ; END IF
195            END DO
196            !
197            ! exchange seed
198            nseed=SUM(rseedmap); IF( lk_mpp )   CALL mpp_sum( nseed )  ! this is the sum of all the point check this iteration
199            !
200            rseedmap_b(:,:)=rseedmap(:,:)
201            CALL lbc_lnk(rseedmap,'T',1.)
202            !
203            ! build new list of seed
204            DO ii=1,jpi
205               DO jj=1,jpj
206                  IF (rseedmap(ii,jj) > 0.0 .AND. rseedmap(ii,jj) /= rseedmap_b(ii,jj) .AND. rdta(ii,jj) > 0 .AND. rdta(ii,jj) /= rfill) THEN
207                     CALL add_head_idx(seed, ii, jj, 1)
208                  END IF
209               END DO
210            END DO 
211            !
212            ! reset map of seed
213            rseedmap(:,:)=0.0
214            !
215         END DO
216      ELSE
217         WRITE(numout,*) 'W A R N I N G : SEED (',ii,ij,') is on land, nothing to do'
218      END IF
219 
220   END SUBROUTINE FillPool2D
221   !
222   !
223   ! subroutine to deals with link list
224   !
225   SUBROUTINE create_idx(pt_idx, ki, kj, kk)
226      TYPE (idx), POINTER :: pt_idx
227      INTEGER, INTENT(in) :: ki, kj, kk
228      !
229      ! initialised all field to NULL()
230      NULLIFY(pt_idx)
231      !
232      ! allocate new element
233      ALLOCATE(pt_idx)
234      pt_idx%i=ki ; pt_idx%j=kj ; pt_idx%k=kk ;
235      pt_idx%next => NULL() 
236   END SUBROUTINE create_idx
237
238   SUBROUTINE add_head_idx(pt_idx, ki, kj, kk)
239      TYPE (idx), POINTER :: pt_idx
240      TYPE (idx), POINTER :: zpt_new
241      INTEGER, INTENT(in) :: ki, kj, kk
242      !
243      ! allocate new element
244      ALLOCATE(zpt_new)
245      zpt_new%i=ki ; zpt_new%j=kj ; zpt_new%k=kk ;
246      zpt_new%next => NULL()
247      !
248      ! linked new element to main list
249      zpt_new%next => pt_idx
250      pt_idx => zpt_new
251   END SUBROUTINE add_head_idx
252
253   SUBROUTINE del_idx(pt_idx)
254      TYPE (idx), POINTER :: pt_idx
255      TYPE (idx), POINTER :: zpt_tmp
256      !
257      ! delete first element
258      IF (ASSOCIATED(pt_idx%next)) THEN
259         zpt_tmp => pt_idx%next
260         DEALLOCATE(pt_idx)
261         NULLIFY(pt_idx)
262         pt_idx => zpt_tmp
263      ELSE
264         NULLIFY(pt_idx)
265      ENDIF
266   END SUBROUTINE del_idx
267   !
268END MODULE domutl
Note: See TracBrowser for help on using the repository browser.