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 NEMO/branches/2019/ENHANCE-03_domcfg/src – NEMO

source: NEMO/branches/2019/ENHANCE-03_domcfg/src/domutl.F90 @ 11604

Last change on this file since 11604 was 11604, checked in by mathiot, 5 years ago

ENHANCE-03_domcfg: remove useless variable in domclo and domutil + cosmetics changes (ticket #2143)
File size: 11.8 KB
Line 
1MODULE domutil
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,                          INTENT(in)    :: kiseed, kjseed, kkseed ! seed
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   ! map of seed (use for processor communication)
51 
52      INTEGER :: ii  , ij  , ik   , kii, kjj, jj, kk    ! working integer
53      INTEGER :: iip1, ijp1, ikp1                       ! working integer
54      INTEGER :: iim1, ijm1, ikm1                       ! working integer
55      INTEGER :: nseed                                  ! size of the stack
56      TYPE (idx), POINTER :: seed
57      !!----------------------------------------------------------------------
58      !
59      ! initialisation seed
60      NULLIFY(seed)
61      !
62      ! create the first seed and update nseed for all processors
63      nseed = 0
64      rseedmap = 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('domutil', nseed )  ! nseed =0 means on land => WARNING later on
75      !
76      ! loop until the stack 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               ii=seed%i; ij=seed%j ; ik=seed%k ; rseedmap(ii,ij,ik)=1.
82               !
83               ! update bathy and update stack size
84               CALL del_head_idx(seed) 
85               !
86               ! check neighbour cells
87               iip1=ii+1 ; IF (iip1 == jpi+1 ) iip1=ii
88               iim1=ii-1 ; IF (iim1 == 0     ) iim1=ii
89               ijp1=ij+1 ; IF (ijp1 == jpj+1 ) ijp1=ij
90               ijm1=ij-1 ; IF (ijm1 == 0     ) ijm1=ij
91               ikp1=ik+1 ; IF (ikp1 == jpk+1 ) ikp1=ik
92               ikm1=ik-1 ; IF (ikm1 == 0     ) ikm1=ik
93               !
94               ! ji direction
95               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
96               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
97               !
98               ! jj direction
99               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
100               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
101               !
102               ! jk direction
103               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
104               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
105            END DO
106            !
107            ! exchange seed
108            nseed=SUM(rseedmap); IF( lk_mpp )   CALL mpp_sum('domutil', nseed )  ! this is the sum of all the point check this iteration
109            !
110            rseedmap_b(:,:,:)=rseedmap(:,:,:)
111            CALL lbc_lnk('domutil', rseedmap, 'T', 1.)
112            !
113            ! build new list of seed
114            DO ii=1,jpi
115               DO jj=1,jpj
116                  DO kk=1,jpk
117                     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
118                        CALL add_head_idx(seed,ii,jj,kk)
119                     END IF
120                  END DO
121               END DO
122            END DO
123            !
124            ! reset map of seed
125            rseedmap(:,:,:)=0.0
126            !
127         END DO
128      ELSE
129         WRITE(numout,*) 'W A R N I N G : SEED (',ii,ij,') is on land, nothing to do'
130      END IF
131 
132   END SUBROUTINE FillPool3D
133
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,                      INTENT(in)    :: kiseed, kjseed ! seed
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          ! location of new seed (used for processor exchange)
147 
148      INTEGER :: ii  , ij  , jj, kii, kjj     ! working integer
149      INTEGER :: iip1, ijp1                   ! working integer
150      INTEGER :: iim1, ijm1                   ! working integer
151      INTEGER :: nseed                        ! size of the stack
152      TYPE (idx), POINTER :: seed
153      !!----------------------------------------------------------------------
154      !
155      ! initialisation seed
156      NULLIFY(seed)
157      !
158      ! create the first seed and update nseed for all processors
159      nseed = 0
160      rseedmap = 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('domutil', nseed )  ! nseed =0 means on land => WARNING later on
171      !
172      ! loop until the stack 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               ii=seed%i; ij=seed%j ; rseedmap(ii,ij)=1.
178               ! update stack size
179               CALL del_head_idx(seed) 
180               !
181               ! check neighbour cells
182               iip1=ii+1 ; IF (iip1 == jpi+1 ) iip1=ii
183               iim1=ii-1 ; IF (iim1 == 0     ) iim1=ii
184               ijp1=ij+1 ; IF (ijp1 == jpj+1 ) ijp1=ij
185               ijm1=ij-1 ; IF (ijm1 == 0     ) ijm1=ij
186               !
187               ! ji direction
188               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
189               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
190               !
191               ! jj direction
192               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
193               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
194            END DO
195            !
196            ! exchange seed
197            nseed=SUM(rseedmap); IF( lk_mpp )   CALL mpp_sum('domutil', nseed )  ! this is the sum of all the point check this iteration
198            !
199            rseedmap_b(:,:)=rseedmap(:,:)
200            CALL lbc_lnk('domutil', rseedmap, 'T', 1.)
201            !
202            ! build new list of seed
203            DO ii=1,jpi
204               DO jj=1,jpj
205                  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
206                     CALL add_head_idx(seed, ii, jj, 1)
207                  END IF
208               END DO
209            END DO 
210            !
211            ! reset map of seed
212            rseedmap(:,:)=0.0
213            !
214         END DO
215      ELSE
216         WRITE(numout,*) 'W A R N I N G : SEED (',ii,ij,') is on land, nothing to do', narea
217      END IF
218 
219   END SUBROUTINE FillPool2D
220   !
221   !
222   ! subroutine to deals with link list
223   !
224   SUBROUTINE add_head_idx(pt_idx, ki, kj, kk)
225      !!---------------------------------------------------------------------
226      !!                  ***  ROUTINE add_head_idx  ***
227      !!
228      !! ** Purpose : add one element in the linked list
229      !!
230      !! ** Method  :  allocate one element, then point %next to the linked list
231      !!----------------------------------------------------------------------
232      TYPE (idx), POINTER :: pt_idx
233      TYPE (idx), POINTER :: zpt_new
234      INTEGER, INTENT(in) :: ki, kj, kk
235      !
236      ! allocate new element
237      ALLOCATE(zpt_new)
238      zpt_new%i=ki ; zpt_new%j=kj ; zpt_new%k=kk ;
239      zpt_new%next => NULL()
240      !
241      ! linked new element to main list
242      zpt_new%next => pt_idx
243      pt_idx => zpt_new
244   END SUBROUTINE add_head_idx
245
246   SUBROUTINE del_head_idx(pt_idx)
247      !!---------------------------------------------------------------------
248      !!                  ***  ROUTINE del_head_idx  ***
249      !!
250      !! ** Purpose : delete one element in the linked list
251      !!
252      !! ** Method  : move the pointer to the next node
253      !!----------------------------------------------------------------------
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_head_idx
267   !
268END MODULE domutil
Note: See TracBrowser for help on using the repository browser.