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

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

ticket #2129 : cleaning domcfg
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#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, lagrif=.FALSE. )
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            i_min(1)=max(i_min(1),1)
164            !     
165            IF(i_min(1)-2 >= 1 .AND. i_max(1)+3 <= SIZE(lon_new1D,1) ) THEN
166               i_min(1) = i_min(1)-2
167               i_max(1) = i_max(1)+3
168            ENDIF
169            tabdim1 = i_max(1) - i_min(1) + 1
170            !
171            IF(j_min(1)-2 >= 1 .AND. j_max(1)+3 <= SIZE(lat_new1D,1) ) THEN
172               j_min(1) = j_min(1)-2
173               j_max(1) = j_max(1)+3
174            ENDIF
175            tabdim2 = j_max(1) - j_min(1) + 1
176            !
177            ALLOCATE(coarselon  (tabdim1,tabdim2), STAT=status(1)) 
178            ALLOCATE(coarselat  (tabdim1,tabdim2), STAT=status(2)) 
179            ALLOCATE(coarsebathy(tabdim1,tabdim2), STAT=status(3)) 
180
181            IF( SUM(status) /= 0 ) CALL ctl_stop( 'STOP', 'dom_bat : unable to allocate arrays' ) 
182            !
183            DO jj = 1,tabdim2
184               coarselon(:,jj) = lon_new1D(i_min(1):i_max(1))
185            END DO
186            !     
187            DO ji = 1,tabdim1
188              coarselat(ji,:) = lat_new1D(j_min(1):j_max(1)) 
189            END DO
190            !
191            CALL iom_get(inum,jpdom_unknown,bathyname,coarsebathy, &
192            &                  kstart=(/i_min(1),j_min(1)/),kcount=(/tabdim1,tabdim2/))
193            !
194         ENDIF   ! > 180
195   
196         DEALLOCATE(lon_new1D) ; DEALLOCATE(lat_new1D)
197         CALL iom_close (inum)
198         
199         coarsebathy = coarsebathy *rn_scale
200        ! reset land to 0
201         WHERE (coarsebathy < 0.)       
202          coarsebathy=0.
203         ENDWHERE 
204 
205      ENDIF   ! external
206
207      IF(lwp) THEN
208         WRITE(numout,*) 'Interpolation of high resolution bathymetry on child grid'
209         IF( nn_interp == 0 ) THEN
210            WRITE(numout,*) 'Arithmetic average ...'
211         ELSE IF( nn_interp == 1 ) THEN
212            WRITE(numout,*) 'Median average ...'
213         ELSE IF( nn_interp == 2 ) THEN     
214            WRITE(numout,*) 'Bilinear interpolation ...'
215         ELSE     
216            WRITE(*,*) 'bad value for nn_interp variable ( must be 0,1 or 2 )'
217            STOP
218         ENDIF
219      ENDIF 
220      bathy(:,:) = 0.
221      !
222      !------------------------------------
223      !MEDIAN AVERAGE or ARITHMETIC AVERAGE
224      !------------------------------------
225      !
226      IF( nn_interp == 0 .OR. nn_interp == 1 ) THEN 
227         !
228         ALLOCATE(trouble_points(jpi,jpj))
229         trouble_points = 0
230         !
231         !  POINT DETECTION
232         !
233         !                       
234         DO jj = 1,jpj
235            DO ji = 1,jpi
236               !     
237               ! FINE GRID CELLS DEFINITION 
238               ji1=max(ji-1,1)
239               jj1=max(jj-1,1)             
240             
241               !
242               Cell_lonmin(ji,jj) = MIN(zglamf(ji1,jj1),zglamf(ji,jj1),zglamf(ji,jj),zglamf(ji1,jj))
243               Cell_lonmax(ji,jj) = MAX(zglamf(ji1,jj1),zglamf(ji,jj1),zglamf(ji,jj),zglamf(ji1,jj))
244               Cell_latmin(ji,jj) = MIN(gphif(ji1,jj1),gphif(ji,jj1),gphif(ji,jj),gphif(ji1,jj))
245               Cell_latmax(ji,jj) = MAX(gphif(ji1,jj1),gphif(ji,jj1),gphif(ji,jj),gphif(ji1,jj))
246               IF( ABS(Cell_lonmax(ji,jj) - Cell_lonmin(ji,jj) ) > 180 ) THEN
247                    zdel = Cell_lonmax(ji,jj)
248                    Cell_lonmax(ji,jj) = Cell_lonmin(ji,jj)
249                    Cell_lonmin(ji,jj) = zdel-360
250               ENDIF
251               !               
252               ! SEARCH FOR ALL POINTS CONTAINED IN THIS CELL
253               !   
254     !      ENDDO
255     !    ENDDO   
256      !   CALL lbc_lnk( 'dom_bat', Cell_lonmin, 'T', 1. )
257      !   CALL lbc_lnk( 'dom_bat', Cell_lonmax, 'T', 1. )
258      !   CALL lbc_lnk( 'dom_bat', Cell_latmin, 'T', 1. )
259      !   CALL lbc_lnk( 'dom_bat', Cell_latmax, 'T', 1. )
260
261
262      !   DO jj = 2,jpj
263      !      DO ji = 2,jpi
264               iimin = 1
265               DO WHILE( coarselon(iimin,1) < Cell_lonmin(ji,jj) ) 
266                  iimin = iimin + 1
267             !     IF ( iimin .LE. 1 ) THEN
268             !     iimin = 1
269             !     EXIT
270             !     ENDIF
271               ENDDO
272               !               
273               jjmin = 1
274               DO WHILE( coarselat(iimin,jjmin) < Cell_latmin(ji,jj) ) 
275                  jjmin = jjmin + 1
276             !     IF ( iimin .LE. 1 ) THEN
277             !     iimin = 1
278             !     EXIT
279             !     ENDIF
280               ENDDO
281               jjmin=max(1,jjmin)
282               !               
283               iimax = iimin 
284               DO WHILE( coarselon(iimax,1)<= Cell_lonmax(ji,jj) ) 
285                  iimax = iimax + 1
286                  IF ( iimax .GE. SIZE(coarsebathy,1) ) THEN
287                  iimax = MIN( iimax,SIZE(coarsebathy,1))
288                  EXIT
289                ENDIF
290               ENDDO
291               !                               
292               jjmax = jjmin 
293               DO WHILE( coarselat(iimax,jjmax) <= Cell_latmax(ji,jj) ) 
294                  jjmax = jjmax + 1
295                  IF ( jjmax .GE. SIZE(coarsebathy,2) ) THEN
296                  jjmax = MIN( jjmax,SIZE(coarsebathy,2))
297                  EXIT
298                ENDIF
299               ENDDO
300               !
301            !   iimax = iimax-1
302            !   jjmax = jjmax-1
303               !               
304               iimin = MAX( iimin,1 )
305               jjmin = MAX( jjmin,1 )
306               iimax = MIN( iimax,SIZE(coarsebathy,1))
307               jjmax = MIN( jjmax,SIZE(coarsebathy,2))
308
309               nxhr = iimax - iimin + 1
310               nyhr = jjmax - jjmin + 1   
311
312
313 
314               IF( nxhr == 0 .OR. nyhr == 0 ) THEN
315                  trouble_points(ji,jj) = 1
316               ELSE
317                  ALLOCATE( vardep(nxhr,nyhr) )
318                  ALLOCATE( mask_oce(nxhr,nyhr) )
319                  mask_oce = 0       
320                  vardep(:,:) = coarsebathy(iimin:iimax,jjmin:jjmax)
321                  WHERE( vardep(:,:) .GT. 0. ) 
322                     mask_oce = 1
323                  ENDWHERE
324                  nxyhr = nxhr*nyhr
325!                 IF( SUM(mask_oce) == 0 ) THEN
326                   IF( SUM(mask_oce) < 0.5*(nxhr*nyhr) ) THEN
327                     bathy(ji,jj) = 0.
328                   ELSE
329                     IF( nn_interp == 0 ) THEN
330                        ! Arithmetic average                   
331                        bathy(ji,jj) = SUM (vardep(:,:)*mask_oce(:,:))/SUM(mask_oce)
332                     ELSE
333                        ! Median average       
334                        !
335                        ALLOCATE(vardep1d(nxhr*nyhr))
336                        vardep1d = RESHAPE(vardep,(/ nxhr*nyhr /) )
337                       ! CALL ssort(vardep1d,nxhr*nyhr)
338                        CALL quicksort(vardep1d,1,nxhr*nyhr)
339                        !
340                        ! Calculate median
341                        !
342                        IF (MOD(SUM(mask_oce),2) .NE. 0) THEN
343                           bathy(ji,jj) = vardep1d( nxyhr/2 + 1 )
344                        ELSE
345                           bathy(ji,jj) =0.5 * ( vardep1d(nxyhr/2) + vardep1d(nxyhr/2+1) )
346                        END IF
347                        DEALLOCATE(vardep1d)   
348                     ENDIF
349                  ENDIF
350                  DEALLOCATE (mask_oce,vardep)
351               ENDIF
352            ENDDO
353         ENDDO
354
355         IF( SUM( trouble_points ) > 0 ) THEN
356            CALL ctl_warn ('too much empty cells, proceed to bilinear interpolation')
357            nn_interp = 2
358         ENDIF
359      ENDIF
360
361     !
362     ! create logical array masksrc
363     !
364      IF( nn_interp == 2) THEN 
365         !
366         identical_grids = .FALSE.
367
368# ifdef key_agrif
369         IF( Agrif_Parent(jpi) == jpi  .AND.   &
370             Agrif_Parent(jpj) == jpj  ) THEN
371            IF( MAXVAL( ABS(coarselat(:,:)- gphit(:,:)) ) < 0.0001 .AND.   &
372                 MAXVAL( ABS(coarselon(:,:)- gphit(:,:)) ) < 0.0001 ) THEN
373               bathy(:,:)  = coarsebathy(:,:) 
374                IF(lwp) WRITE(numout,*) 'same grid between ', bathyname,' and child domain'                   
375               identical_grids = .TRUE.                         
376            ENDIF
377         ENDIF
378# endif
379         IF( .NOT. identical_grids ) THEN 
380            ALLOCATE(masksrc(SIZE(coarsebathy,1),SIZE(coarsebathy,2)))
381            ALLOCATE(bathy_test(jpi,jpj))
382            !
383            !                    Where(G0%bathy_meter.le.0.00001)
384            !                  masksrc = .false.
385            !              ElseWhere
386            !
387            masksrc = .TRUE.
388            !
389            !              End where                       
390            !           
391            ! compute remapping matrix thanks to SCRIP package           
392            !                                 
393            CALL get_remap_matrix(coarselat,gphit,   &
394                 coarselon,zglamt,   &
395                 masksrc,matrix,src_add,dst_add)
396            CALL make_remap(coarsebathy,bathy_test,jpi,jpj, &
397                 matrix,src_add,dst_add) 
398            !                                 
399            bathy= bathy_test               
400            !           
401            DEALLOCATE(masksrc)
402            DEALLOCATE(bathy_test) 
403         ENDIF
404        !           
405      ENDIF
406      CALL lbc_lnk( 'dom_bat', bathy, 'T', 1. )
407
408       ! Correct South and North
409#if defined key_agrif
410      IF( ln_bry_south ) THEN 
411         IF( (nbondj == -1).OR.(nbondj == 2) ) THEN
412           bathy(:,1)=bathy(:,2)
413         ENDIF
414      ELSE
415            bathy(:,1) = 0.
416      ENDIF
417#else
418      IF ((nbondj == -1).OR.(nbondj == 2)) THEN
419            bathy(:,1)=bathy(:,2)
420      ENDIF
421#endif
422
423      IF ((nbondj == 1).OR.(nbondj == 2)) THEN
424         bathy(:,jpj)=bathy(:,jpj-1)
425      ENDIF
426
427      ! Correct West and East
428      IF (jperio /=1) THEN
429         IF ((nbondi == -1).OR.(nbondi == 2)) THEN
430            bathy(1,:)=bathy(2,:)
431         ENDIF
432         IF ((nbondi == 1).OR.(nbondi == 2)) THEN
433         bathy(jpi,:)=bathy(jpi-1,:)
434         ENDIF
435      ENDIF
436
437
438   END SUBROUTINE dom_bat
439
440END MODULE dombat
Note: See TracBrowser for help on using the repository browser.