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_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src – NEMO

source: utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/dombat.F90 @ 13024

Last change on this file since 13024 was 13024, checked in by rblod, 4 years ago

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