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.
dombat.F90 in utils/tools/DOMAINcfg/src – NEMO

source: utils/tools/DOMAINcfg/src/dombat.F90 @ 15331

Last change on this file since 15331 was 15331, checked in by jchanut, 3 years ago

#2638, add closed seas filling algorithm (inside AGRIF zooms only)
File size: 18.8 KB
Line 
1MODULE dombat
2
3   USE dom_oce           ! ocean domain
4   !
5   USE in_out_manager    ! I/O manager
6   USE iom               ! I/O library
7   USE lbclnk            ! ocean lateral boundary conditions (or mpp link)
8   USE lib_mpp           ! distributed memory computing library
9#if defined key_agrif
10   USE agrif_modutil
11   USE agrif_parameters
12#endif   
13   USE bilinear_interp
14   USE dom_oce
15
16   IMPLICIT NONE
17   PRIVATE
18
19   PUBLIC   dom_bat        ! called by dom_zgr.F90
20
21CONTAINS
22
23   SUBROUTINE dom_bat
24
25      INTEGER :: inum, id, ji, jj,ji1,jj1
26      INTEGER :: iimin,iimax,jjmin,jjmax
27      INTEGER :: tabdim1, tabdim2, nxhr, nyhr, nxyhr
28      INTEGER :: nbadd, istart, iend, jstart, jend
29      INTEGER, DIMENSION(2) :: ddims
30      INTEGER, DIMENSION(3) :: status
31      INTEGER, DIMENSION(1) :: i_min,i_max
32      INTEGER, DIMENSION(1) :: j_min,j_max
33      INTEGER, DIMENSION(:)  ,ALLOCATABLE     ::   src_add,dst_add 
34      INTEGER, DIMENSION(:,:), ALLOCATABLE :: trouble_points , vardep,mask_oce
35      REAL(wp) ,DIMENSION(jpi,jpj):: Cell_lonmin, Cell_lonmax, Cell_latmin, Cell_latmax
36      REAL(wp) ::zdel
37      REAL(wp), DIMENSION(:)  , ALLOCATABLE ::   lon_new1D , lat_new1D, vardep1d
38      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   coarselon, coarselat, coarsebathy, bathy_test
39      REAL(wp), DIMENSION(:,:),ALLOCATABLE     ::   matrix,interpdata 
40      LOGICAL :: lonlat_2D, ln_pacifiq
41      LOGICAL :: identical_grids
42      LOGICAL, DIMENSION(:,:), ALLOCATABLE     ::   masksrc
43      REAL(wp), DIMENSION(jpi,jpj) :: zglamt, zgphi, zglamu, zglamv, zglamf
44      REAL(wp) :: zshift
45 
46      CHARACTER(32) :: bathyfile, bathyname, lonname,latname       
47
48      bathyfile=TRIM(cn_topo)
49      bathyname=TRIM(cn_bath)
50      lonname=TRIM(cn_lon)
51      latname=TRIM(cn_lat)
52   
53      CALL iom_open( bathyfile, inum, ldiof=.TRUE. )
54     
55      ! check if lon/lat are 2D arrays
56      id = iom_varid( inum, lonname, ddims )
57      IF (ddims(2)==0) THEN
58         lonlat_2D = .FALSE.
59      ELSE
60         lonlat_2D = .TRUE.
61      ENDIF   
62     
63      id = iom_varid( inum, bathyname, ddims )
64      ln_pacifiq = .FALSE.
65      zglamt(:,:) = glamt(:,:)
66      zglamu(:,:) = glamu(:,:)
67      zglamv(:,:) = glamv(:,:)
68      zglamf(:,:) = glamf(:,:)
69
70      IF( glamt(1,1) .GT. glamt(jpi,jpj) ) ln_pacifiq =.false. 
71
72      zshift = 0.
73      IF( ln_pacifiq  ) THEN
74         zshift = 0.!Abs(minval(glamt)) +0.1
75         WHERE ( glamt < 0 )
76            zglamt = zglamt + zshift + 360.
77         END WHERE
78         WHERE ( glamu < 0 )
79            zglamu = zglamu + zshift +360.
80         END WHERE
81         WHERE ( glamv < 0 )
82            zglamv = zglamv + zshift +360.
83         END WHERE
84         WHERE ( glamf < 0 )
85            zglamf = zglamf + zshift +360.
86         END WHERE
87      ENDIF
88
89      status=-1
90
91      IF (lonlat_2D) THEN
92         ! here implicitly it's old topo database (orca format)
93         ALLOCATE(coarselon  (ddims(1),ddims(2)), STAT=status(1)) 
94         ALLOCATE(coarselat  (ddims(1),ddims(2)), STAT=status(2)) 
95         ALLOCATE(coarsebathy(ddims(1),ddims(2)), STAT=status(3)) 
96         IF( sum(status) /= 0 )   CALL ctl_stop( 'STOP', 'dom_bat : unable to allocate arrays' )
97         CALL iom_get  ( inum, jpdom_unknown, lonname, coarselon )
98         CALL iom_get  ( inum, jpdom_unknown, latname, coarselat )
99         CALL iom_get  ( inum, jpdom_unknown, bathyname, coarsebathy )
100         CALL iom_close (inum)
101         IF( ln_pacifiq ) THEN
102       !     WHERE(coarselon < 0.00001)
103                coarselon = Coarselon + zshift
104       !      END WHERE
105         ENDIF     
106         ! equivalent to new database
107      ELSE
108         ALLOCATE(lon_new1D(ddims(1)), lat_new1D(ddims(2)))
109         CALL iom_get  ( inum, jpdom_unknown, lonname, lon_new1D )
110         CALL iom_get  ( inum, jpdom_unknown, latname, lat_new1D )
111         IF( ln_pacifiq ) THEN
112            WHERE(lon_new1D < 0.00001) 
113                lon_new1D = lon_new1D +360.!zshift
114             END WHERE
115         ENDIF               
116         zdel =  0.00   
117         IF( MAXVAL(zglamf) > 180 + zshift ) THEN 
118            !         
119         !   WHERE( lon_new1D < 0 )
120         !       lon_new1D = lon_new1D + 360.
121         !   END WHERE
122            !     
123            i_min = MAXLOC(lon_new1D,mask = lon_new1D < MINVAL(zglamf(1:jpi-1,1:jpj-1)) )
124            i_max = MINLOC(lon_new1D,mask = lon_new1D > MAXVAL(zglamf(1:jpi-1,1:jpj-1)) )                   
125            j_min = MAXLOC(lat_new1D,mask = lat_new1D < MINVAL( gphif(1:jpi-1,1:jpj-1)) )
126            j_max = MINLOC(lat_new1D,mask = lat_new1D > MAXVAL( gphif(1:jpi-1,1:jpj-1)) )
127            !
128            tabdim1 = ( SIZE(lon_new1D) - i_min(1) + 1 ) + i_max(1)                   
129            !         
130            IF(j_min(1)-2 >= 1 .AND. j_max(1)+3 <= SIZE(lat_new1D,1) ) THEN
131               j_min(1) = j_min(1) - 2
132               j_max(1) = j_max(1)+ 3
133            ENDIF
134            tabdim2 = j_max(1) - j_min(1) + 1
135            !
136            ALLOCATE(coarselon  (tabdim1,tabdim2), STAT=status(1))
137            ALLOCATE(coarselat  (tabdim1,tabdim2), STAT=status(2))
138            ALLOCATE(Coarsebathy(tabdim1,tabdim2), STAT=status(3)) 
139
140            IF( SUM(status) /= 0 ) CALL ctl_stop( 'STOP', 'dom_bat : unable to allocate arrays' )         
141            !
142            DO ji = 1,tabdim1
143               coarselat(ji,:) = lat_new1D(j_min(1):j_max(1))
144            END DO
145           
146            !
147            DO jj = 1, tabdim2                                 
148               coarselon(1:SIZE(lon_new1D)-i_min(1)+1      ,jj) = lon_new1D(i_min(1):SIZE(lon_new1D))
149               coarselon(2+SIZE(lon_new1D)-i_min(1):tabdim1,jj) = lon_new1D(1:i_max(1))   
150            END DO
151            !
152            CALL iom_get(inum, jpdom_unknown, bathyname,coarsebathy(1:SIZE(lon_new1D)-i_min(1)+1,:), &
153            kstart=(/i_min(1),j_min(1)/), kcount=(/SIZE(lon_new1D)-i_min(1)+1,tabdim2/))   ! +1?
154            !
155            CALL iom_get(inum, jpdom_unknown, bathyname,coarsebathy(2+SIZE(lon_new1D)-i_min(1):tabdim1,:), &
156            kstart=(/1,j_min(1)/),kcount=(/i_max(1),tabdim2/))               
157            !
158         ELSE
159          !  WHERE( lon_new1D > (180. + zshift) )   lon_new1D = lon_new1D - 360.
160            i_min = MAXLOC(lon_new1D,mask = lon_new1D < MINVAL(zglamf)-zdel)
161            i_max = MINLOC(lon_new1D,mask = lon_new1D > MAXVAL(zglamf)+zdel)
162            j_min = MAXLOC(lat_new1D,mask = lat_new1D < MINVAL(gphif)-zdel)
163            j_max = MINLOC(lat_new1D,mask = lat_new1D > MAXVAL(gphif)+zdel)
164
165            ! jc: to prevent from issues with grids that exactly pass through
166            ! longitude= +-180:
167            IF (i_max(1)==0) i_max(1)=ddims(1)
168         
169            i_min(1)=max(i_min(1),1)
170            !     
171            IF(i_min(1)-2 >= 1 .AND. i_max(1)+3 <= SIZE(lon_new1D,1) ) THEN
172               i_min(1) = i_min(1)-2
173               i_max(1) = i_max(1)+3
174            ENDIF
175            tabdim1 = i_max(1) - i_min(1) + 1
176            !
177            IF(j_min(1)-2 >= 1 .AND. j_max(1)+3 <= SIZE(lat_new1D,1) ) THEN
178               j_min(1) = j_min(1)-2
179               j_max(1) = j_max(1)+3
180            ENDIF
181            tabdim2 = j_max(1) - j_min(1) + 1
182            !
183            ALLOCATE(coarselon  (tabdim1,tabdim2), STAT=status(1)) 
184            ALLOCATE(coarselat  (tabdim1,tabdim2), STAT=status(2)) 
185            ALLOCATE(coarsebathy(tabdim1,tabdim2), STAT=status(3)) 
186
187            IF( SUM(status) /= 0 ) CALL ctl_stop( 'STOP', 'dom_bat : unable to allocate arrays' ) 
188            !
189            DO jj = 1,tabdim2
190               coarselon(:,jj) = lon_new1D(i_min(1):i_max(1))
191            END DO
192            !     
193            DO ji = 1,tabdim1
194              coarselat(ji,:) = lat_new1D(j_min(1):j_max(1)) 
195            END DO
196            !
197            CALL iom_get(inum,jpdom_unknown,bathyname,coarsebathy, &
198            &                  kstart=(/i_min(1),j_min(1)/),kcount=(/tabdim1,tabdim2/))
199            !
200         ENDIF   ! > 180
201   
202         DEALLOCATE(lon_new1D) ; DEALLOCATE(lat_new1D)
203         CALL iom_close (inum)
204         
205         coarsebathy = coarsebathy *rn_scale
206        ! reset land to 0
207         WHERE (coarsebathy < 0.)       
208          coarsebathy=0.
209         ENDWHERE 
210 
211      ENDIF   ! external
212
213      IF(lwp) THEN
214         WRITE(numout,*) 'Interpolation of high resolution bathymetry on child grid'
215         IF( nn_interp == 0 ) THEN
216            WRITE(numout,*) 'Arithmetic average ...'
217         ELSE IF( nn_interp == 1 ) THEN
218            WRITE(numout,*) 'Median average ...'
219         ELSE IF( nn_interp == 2 ) THEN     
220            WRITE(numout,*) 'Bilinear interpolation ...'
221         ELSE     
222            WRITE(*,*) 'bad value for nn_interp variable ( must be 0,1 or 2 )'
223            STOP
224         ENDIF
225      ENDIF 
226      bathy(:,:) = 0.
227      !
228      !------------------------------------
229      !MEDIAN AVERAGE or ARITHMETIC AVERAGE
230      !------------------------------------
231      !
232      IF( nn_interp == 0 .OR. nn_interp == 1 ) THEN 
233         !
234         ALLOCATE(trouble_points(jpi,jpj))
235         trouble_points = 0
236         !
237         !  POINT DETECTION
238         !
239         !                       
240         DO jj = 1,jpj
241            DO ji = 1,jpi
242               !     
243               ! FINE GRID CELLS DEFINITION 
244               ji1=max(ji-1,1)
245               jj1=max(jj-1,1)             
246             
247               !
248               Cell_lonmin(ji,jj) = MIN(zglamf(ji1,jj1),zglamf(ji,jj1),zglamf(ji,jj),zglamf(ji1,jj))
249               Cell_lonmax(ji,jj) = MAX(zglamf(ji1,jj1),zglamf(ji,jj1),zglamf(ji,jj),zglamf(ji1,jj))
250               Cell_latmin(ji,jj) = MIN(gphif(ji1,jj1),gphif(ji,jj1),gphif(ji,jj),gphif(ji1,jj))
251               Cell_latmax(ji,jj) = MAX(gphif(ji1,jj1),gphif(ji,jj1),gphif(ji,jj),gphif(ji1,jj))
252               IF( ABS(Cell_lonmax(ji,jj) - Cell_lonmin(ji,jj) ) > 180 ) THEN
253                    zdel = Cell_lonmax(ji,jj)
254                    Cell_lonmax(ji,jj) = Cell_lonmin(ji,jj)
255                    Cell_lonmin(ji,jj) = zdel-360
256               ENDIF
257               !               
258               ! SEARCH FOR ALL POINTS CONTAINED IN THIS CELL
259               !   
260     !      ENDDO
261     !    ENDDO   
262      !   CALL lbc_lnk( 'dom_bat', Cell_lonmin, 'T', 1. )
263      !   CALL lbc_lnk( 'dom_bat', Cell_lonmax, 'T', 1. )
264      !   CALL lbc_lnk( 'dom_bat', Cell_latmin, 'T', 1. )
265      !   CALL lbc_lnk( 'dom_bat', Cell_latmax, 'T', 1. )
266
267
268      !   DO jj = 2,jpj
269      !      DO ji = 2,jpi
270               iimin = 1
271               DO WHILE( coarselon(iimin,1) < Cell_lonmin(ji,jj) ) 
272                  iimin = iimin + 1
273             !     IF ( iimin .LE. 1 ) THEN
274             !     iimin = 1
275             !     EXIT
276             !     ENDIF
277               ENDDO
278               !               
279               jjmin = 1
280               DO WHILE( coarselat(iimin,jjmin) < Cell_latmin(ji,jj) ) 
281                  jjmin = jjmin + 1
282             !     IF ( iimin .LE. 1 ) THEN
283             !     iimin = 1
284             !     EXIT
285             !     ENDIF
286               ENDDO
287               jjmin=max(1,jjmin)
288               !               
289               iimax = iimin 
290               DO WHILE( coarselon(iimax,1)<= Cell_lonmax(ji,jj) ) 
291                  iimax = iimax + 1
292                  IF ( iimax .GE. SIZE(coarsebathy,1) ) THEN
293                  iimax = MIN( iimax,SIZE(coarsebathy,1))
294                  EXIT
295                ENDIF
296               ENDDO
297               !                               
298               jjmax = jjmin 
299               DO WHILE( coarselat(iimax,jjmax) <= Cell_latmax(ji,jj) ) 
300                  jjmax = jjmax + 1
301                  IF ( jjmax .GE. SIZE(coarsebathy,2) ) THEN
302                  jjmax = MIN( jjmax,SIZE(coarsebathy,2))
303                  EXIT
304                ENDIF
305               ENDDO
306               !
307            !   iimax = iimax-1
308            !   jjmax = jjmax-1
309               !               
310               iimin = MAX( iimin,1 )
311               jjmin = MAX( jjmin,1 )
312               iimax = MIN( iimax,SIZE(coarsebathy,1))
313               jjmax = MIN( jjmax,SIZE(coarsebathy,2))
314
315               nxhr = iimax - iimin + 1
316               nyhr = jjmax - jjmin + 1   
317
318
319 
320               IF( nxhr == 0 .OR. nyhr == 0 ) THEN
321                  trouble_points(ji,jj) = 1
322               ELSE
323                  ALLOCATE( vardep(nxhr,nyhr) )
324                  ALLOCATE( mask_oce(nxhr,nyhr) )
325                  mask_oce = 0       
326                  vardep(:,:) = coarsebathy(iimin:iimax,jjmin:jjmax)
327                  WHERE( vardep(:,:) .GT. 0. ) 
328                     mask_oce = 1
329                  ENDWHERE
330                  nxyhr = nxhr*nyhr
331!                 IF( SUM(mask_oce) == 0 ) THEN
332                   IF( SUM(mask_oce) < 0.5*(nxhr*nyhr) ) THEN
333                     bathy(ji,jj) = 0.
334                   ELSE
335                     IF( nn_interp == 0 ) THEN
336                        ! Arithmetic average                   
337                        bathy(ji,jj) = SUM (vardep(:,:)*mask_oce(:,:))/SUM(mask_oce)
338                     ELSE
339                        ! Median average       
340                        !
341                        ALLOCATE(vardep1d(nxhr*nyhr))
342                        vardep1d = RESHAPE(vardep,(/ nxhr*nyhr /) )
343                       ! CALL ssort(vardep1d,nxhr*nyhr)
344                        CALL quicksort(vardep1d,1,nxhr*nyhr)
345                        !
346                        ! Calculate median
347                        !
348                        IF (MOD(SUM(mask_oce),2) .NE. 0) THEN
349                           bathy(ji,jj) = vardep1d( nxyhr/2 + 1 )
350                        ELSE
351                           bathy(ji,jj) =0.5 * ( vardep1d(nxyhr/2) + vardep1d(nxyhr/2+1) )
352                        END IF
353                        DEALLOCATE(vardep1d)   
354                     ENDIF
355                  ENDIF
356                  DEALLOCATE (mask_oce,vardep)
357               ENDIF
358            ENDDO
359         ENDDO
360
361         IF( SUM( trouble_points ) > 0 ) THEN
362            CALL ctl_warn ('too much empty cells, proceed to bilinear interpolation')
363            nn_interp = 2
364         ENDIF
365      ENDIF
366
367     !
368     ! create logical array masksrc
369     !
370      IF( nn_interp == 2) THEN 
371         !
372         identical_grids = .FALSE.
373
374# ifdef key_agrif
375         IF( Agrif_Parent(jpi) == jpi  .AND.   &
376             Agrif_Parent(jpj) == jpj  ) THEN
377            IF( MAXVAL( ABS(coarselat(:,:)- gphit(:,:)) ) < 0.0001 .AND.   &
378                 MAXVAL( ABS(coarselon(:,:)- gphit(:,:)) ) < 0.0001 ) THEN
379               bathy(:,:)  = coarsebathy(:,:) 
380                IF(lwp) WRITE(numout,*) 'same grid between ', bathyname,' and child domain'                   
381               identical_grids = .TRUE.                         
382            ENDIF
383         ENDIF
384# endif
385         IF( .NOT. identical_grids ) THEN 
386            ALLOCATE(masksrc(SIZE(coarsebathy,1),SIZE(coarsebathy,2)))
387            ALLOCATE(bathy_test(jpi,jpj))
388            !
389            !                    Where(G0%bathy_meter.le.0.00001)
390            !                  masksrc = .false.
391            !              ElseWhere
392            !
393            masksrc = .TRUE.
394            !
395            !              End where                       
396            !           
397            ! compute remapping matrix thanks to SCRIP package           
398            !                                 
399            CALL get_remap_matrix(coarselat,gphit,   &
400                 coarselon,zglamt,   &
401                 masksrc,matrix,src_add,dst_add)
402            CALL make_remap(coarsebathy,bathy_test,jpi,jpj, &
403                 matrix,src_add,dst_add) 
404            !                                 
405            bathy= bathy_test               
406            !           
407            DEALLOCATE(masksrc)
408            DEALLOCATE(bathy_test) 
409         ENDIF
410        !           
411      ENDIF
412      CALL lbc_lnk( 'dom_bat', bathy, 'T', 1.,kfillmode = jpfillcopy)
413
414#if defined key_agrif
415      IF (ln_remove_closedseas.AND.(.NOT.Agrif_Root())) THEN
416         ALLOCATE(bathy_test(jpi,jpj))
417         bathy_test(:,:) = 0._wp 
418         !
419         ! --- West --- !
420         IF(lk_west) THEN
421            istart = nn_hls + 2
422            iend   = nn_hls + nbghostcells 
423            DO ji = mi0(istart), mi1(iend)
424               DO jj = 1, jpj
425                  IF ( bathy (ji,jj)/=0._wp )   bathy_test(ji,jj) = 1._wp
426               END DO
427            END DO
428         ENDIF
429         !
430         ! --- East --- !
431         IF(lk_east) THEN
432            istart = jpiglo - ( nn_hls + nbghostcells -1 ) 
433            iend   = jpiglo - ( nn_hls + 1 )
434            DO ji = mi0(istart), mi1(iend)
435               DO jj = 1, jpj
436                  IF ( bathy (ji,jj)/=0._wp )   bathy_test(ji,jj) = 1._wp
437               END DO
438            END DO
439         ENDIF
440         !
441         ! --- South --- !
442         IF(lk_south) THEN
443            jstart = nn_hls + 2
444            jend   = nn_hls + nbghostcells
445            DO jj = mj0(jstart), mj1(jend)
446               DO ji = 1, jpi
447                  IF ( bathy (ji,jj)/=0._wp )   bathy_test(ji,jj) = 1._wp
448               END DO
449            END DO
450         ENDIF
451         !
452         ! --- North --- !
453         IF(lk_north) THEN
454            jstart = jpjglo - ( nn_hls + nbghostcells -1 )
455            jend   = jpjglo - ( nn_hls + 1 )
456            DO jj = mj0(jstart), mj1(jend)
457               DO ji = 1, jpi
458                  IF ( bathy (ji,jj)/=0._wp )  bathy_test(ji,jj) = 1._wp
459               END DO
460            END DO
461         ENDIF
462         
463         nbadd = 1
464         DO WHILE ( nbadd/=0 )
465            nbadd = 0
466            DO ji = 1+nn_hls, jpi-nn_hls
467               DO jj = 1+nn_hls, jpj-nn_hls           
468                  IF (bathy(ji,jj) > 0._wp) THEN
469                     IF (MAX(bathy_test(ji,jj+1),bathy_test(ji,jj-1), & 
470                     &       bathy_test(ji-1,jj),bathy_test(ji+1,jj))==1._wp)  THEN
471                        IF (bathy_test(ji,jj)/=1._wp) nbadd = nbadd + 1
472                        bathy_test(ji,jj)=1._wp
473                     ENDIF
474                  ENDIF
475               END DO
476            END DO
477            IF( lk_mpp )   CALL mpp_sum('dom_bat', nbadd )
478            CALL lbc_lnk( 'dom_bat', bathy_test, 'T', 1.,kfillmode = jpfillcopy)
479
480         END DO
481
482         WHERE(bathy_test==0._wp) bathy = 0._wp
483         DEALLOCATE(bathy_test)
484      ENDIF
485#endif
486
487
488       ! Correct South and North
489! #if defined key_agrif
490!       IF( lk_south ) THEN 
491!          IF( (nbondj == -1).OR.(nbondj == 2) ) THEN
492!            bathy(:,1)=bathy(:,2)
493!          ENDIF
494!       ELSE
495!             bathy(:,1) = 0.
496!       ENDIF
497! #else
498!       IF ((nbondj == -1).OR.(nbondj == 2)) THEN
499!             bathy(:,1)=bathy(:,2)
500!       ENDIF
501! #endif
502
503!       IF ((nbondj == 1).OR.(nbondj == 2)) THEN
504!          bathy(:,jpj)=bathy(:,jpj-1)
505!       ENDIF
506
507!       ! Correct West and East
508!       IF (jperio /=1) THEN
509!          IF ((nbondi == -1).OR.(nbondi == 2)) THEN
510!             bathy(1,:)=bathy(2,:)
511!          ENDIF
512!          IF ((nbondi == 1).OR.(nbondi == 2)) THEN
513!          bathy(jpi,:)=bathy(jpi-1,:)
514!          ENDIF
515!       ENDIF
516
517
518   END SUBROUTINE dom_bat
519
520END MODULE dombat
Note: See TracBrowser for help on using the repository browser.