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

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

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