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 @ 11201

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

ENHANCE-03_domcfg : add management of closed seas in domain cfg by flood filling and lat/lon seed instead of i/j box definition (ticket #2143)
File size: 11.2 KB
RevLine 
[11201]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
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      rseedmap = 0.
66      ii=kiseed; ij=kjseed ; ik=kkseed
67      IF ((mi0(ii) == mi1(ii)) .AND. (mj0(ij) == mj1(ij))) THEN   ! T if seed on the local domain
68         ii = mi0(ii) ; ij = mj0(ij);                             ! conversion to local index
69         IF (rdta(ii,ij,ik) > 0 ) THEN                            ! create seed if not on land
70            CALL add_head_idx(seed,ii,ij,ik)
71            rdta    (ii,ij,ik) = rfill
72            rseedmap(ii,ij,ik) = 1.
73         END IF
74      END IF
75      nseed=SUM(rseedmap); IF( lk_mpp )   CALL mpp_sum('domutil', nseed )  ! nseed =0 means on land => WARNING later on
76      !
77      ! loop until the pile size is 0 or if the pool is larger than the critical size
78      IF (nseed > 0) THEN
79         ! seed on ocean continue
80         DO WHILE ( nseed /= 0 )
81            DO WHILE ( ASSOCIATED(seed) )
82               ip=ip+1
83               ii=seed%i; ij=seed%j ; ik=seed%k ; rseedmap(ii,ij,ik)=1.
84               !
85               ! update bathy and update pile size
86               CALL del_head_idx(seed) 
87               !
88               ! check neighbour cells
89               iip1=ii+1 ; IF (iip1 == jpi+1 ) iip1=ii
90               iim1=ii-1 ; IF (iim1 == 0     ) iim1=ii
91               ijp1=ij+1 ; IF (ijp1 == jpj+1 ) ijp1=ij
92               ijm1=ij-1 ; IF (ijm1 == 0     ) ijm1=ij
93               ikp1=ik+1 ; IF (ikp1 == jpk+1 ) ikp1=ik
94               ikm1=ik-1 ; IF (ikm1 == 0     ) ikm1=ik
95               !
96               ! ji direction
97               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
98               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
99               !
100               ! jj direction
101               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
102               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
103               !
104               ! jk direction
105               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
106               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
107            END DO
108            !
109            ! exchange seed
110            nseed=SUM(rseedmap); IF( lk_mpp )   CALL mpp_sum('domutil', nseed )  ! this is the sum of all the point check this iteration
111            !
112            rseedmap_b(:,:,:)=rseedmap(:,:,:)
113            CALL lbc_lnk('domutil', rseedmap, 'T', 1.)
114            !
115            ! build new list of seed
116            DO ii=1,jpi
117               DO jj=1,jpj
118                  DO kk=1,jpk
119                     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
120                        CALL add_head_idx(seed,ii,jj,kk)
121                     END IF
122                  END DO
123               END DO
124            END DO
125            !
126            ! reset map of seed
127            rseedmap(:,:,:)=0.0
128            !
129         END DO
130      ELSE
131         WRITE(numout,*) 'W A R N I N G : SEED (',ii,ij,') is on land, nothing to do'
132      END IF
133 
134   END SUBROUTINE FillPool3D
135
136   SUBROUTINE FillPool2D(kiseed, kjseed, rdta, rfill)
137      !!---------------------------------------------------------------------
138      !!                  ***  ROUTINE FillPool3D  ***
139      !!
140      !! ** Purpose :  Replace all area surrounding by mask value by kifill value
141      !!
142      !! ** Method  :  flood fill algorithm
143      !!
144      !!----------------------------------------------------------------------
145      INTEGER,                    INTENT(in)    :: kiseed, kjseed
146      REAL(wp),                   INTENT(in)    :: rfill    ! filling value
147      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: rdta     ! input data
148      REAL(wp), DIMENSION(jpi,jpj) :: rseedmap, rseedmap_b 
149 
150      INTEGER :: ip=0                     ! size of the pile
151      INTEGER :: ii  , ij  , jj, kii, kjj     ! working integer
152      INTEGER :: iip1, ijp1               ! working integer
153      INTEGER :: iim1, ijm1
154      INTEGER :: nseed
155      TYPE (idx), POINTER :: seed
156      !!----------------------------------------------------------------------
157      !
158      ! initialisation seed
159      NULLIFY(seed)
160      !
161      ! create the first seed and update nseed for all processors
162      nseed = 0
163      rseedmap = 0.
164      ii = kiseed ; ij = kjseed ;
165      IF ((mi0(ii) == mi1(ii)) .AND. (mj0(ij) == mj1(ij))) THEN   ! T if seed on the local domain
166         ii = mi0(ii) ; ij = mj0(ij)                              ! conversion to local index
167         IF (rdta(ii,ij) > 0 ) THEN                               ! create seed if not on land
168            CALL add_head_idx(seed,ii,ij,1)
169            rdta    (ii,ij) = rfill
170            rseedmap(ii,ij) = 1.
171         END IF
172      END IF
173      nseed=SUM(rseedmap); IF( lk_mpp )   CALL mpp_sum('domutil', nseed )  ! nseed =0 means on land => WARNING later on
174      !
175      ! loop until the pile size is 0 or if the pool is larger than the critical size
176      IF (nseed > 0) THEN
177         ! seed on ocean continue
178         DO WHILE ( nseed .NE. 0 )
179            DO WHILE ( ASSOCIATED(seed) )
180               ip=ip+1
181               ii=seed%i; ij=seed%j ; rseedmap(ii,ij)=1.
182               ! update pile size
183               CALL del_head_idx(seed) 
184               !
185               ! check neighbour cells
186               iip1=ii+1 ; IF (iip1 == jpi+1 ) iip1=ii
187               iim1=ii-1 ; IF (iim1 == 0     ) iim1=ii
188               ijp1=ij+1 ; IF (ijp1 == jpj+1 ) ijp1=ij
189               ijm1=ij-1 ; IF (ijm1 == 0     ) ijm1=ij
190               !
191               ! ji direction
192               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
193               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
194               !
195               ! jj direction
196               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
197               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
198            END DO
199            !
200            ! exchange seed
201            nseed=SUM(rseedmap); IF( lk_mpp )   CALL mpp_sum('domutil', nseed )  ! this is the sum of all the point check this iteration
202            !
203            rseedmap_b(:,:)=rseedmap(:,:)
204            CALL lbc_lnk('domutil', rseedmap, 'T', 1.)
205            !
206            ! build new list of seed
207            DO ii=1,jpi
208               DO jj=1,jpj
209                  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
210                     CALL add_head_idx(seed, ii, jj, 1)
211                  END IF
212               END DO
213            END DO 
214            !
215            ! reset map of seed
216            rseedmap(:,:)=0.0
217            !
218         END DO
219      ELSE
220         WRITE(numout,*) 'W A R N I N G : SEED (',ii,ij,') is on land, nothing to do', narea
221      END IF
222 
223   END SUBROUTINE FillPool2D
224   !
225   !
226   ! subroutine to deals with link list
227   !
228   SUBROUTINE create_idx(pt_idx, ki, kj, kk)
229      TYPE (idx), POINTER :: pt_idx
230      INTEGER, INTENT(in) :: ki, kj, kk
231      !
232      ! initialised all field to NULL()
233      NULLIFY(pt_idx)
234      !
235      ! allocate new element
236      ALLOCATE(pt_idx)
237      pt_idx%i=ki ; pt_idx%j=kj ; pt_idx%k=kk ;
238      pt_idx%next => NULL() 
239   END SUBROUTINE create_idx
240
241   SUBROUTINE add_head_idx(pt_idx, ki, kj, kk)
242      TYPE (idx), POINTER :: pt_idx
243      TYPE (idx), POINTER :: zpt_new
244      INTEGER, INTENT(in) :: ki, kj, kk
245      !
246      ! allocate new element
247      ALLOCATE(zpt_new)
248      zpt_new%i=ki ; zpt_new%j=kj ; zpt_new%k=kk ;
249      zpt_new%next => NULL()
250      !
251      ! linked new element to main list
252      zpt_new%next => pt_idx
253      pt_idx => zpt_new
254   END SUBROUTINE add_head_idx
255
256   SUBROUTINE del_head_idx(pt_idx)
257      TYPE (idx), POINTER :: pt_idx
258      TYPE (idx), POINTER :: zpt_tmp
259      !
260      ! delete first element
261      IF (ASSOCIATED(pt_idx%next)) THEN
262         zpt_tmp => pt_idx%next
263         DEALLOCATE(pt_idx)
264         NULLIFY(pt_idx)
265         pt_idx => zpt_tmp
266      ELSE
267         NULLIFY(pt_idx)
268      ENDIF
269   END SUBROUTINE del_head_idx
270   !
271END MODULE domutil
Note: See TracBrowser for help on using the repository browser.