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

source: utils/tools_AGRIF_CMEMS_2020/DOMAINcfg/src/dombat.F90 @ 10727

Last change on this file since 10727 was 10727, checked in by rblod, 5 years ago

new nesting tools (attempt) and brutal cleaning of DOMAINcfg, see ticket #2129
File size: 9.4 KB
Line 
1MODULE dombat
2
3   USE oce               ! ocean variables
4   USE dom_oce           ! ocean domain
5!   USE closea            ! closed seas
6   !
7   USE in_out_manager    ! I/O manager
8   USE iom               ! I/O library
9   USE lbclnk            ! ocean lateral boundary conditions (or mpp link)
10   USE lib_mpp           ! distributed memory computing library
11   USE wrk_nemo          ! Memory allocation
12   USE timing            ! Timing
13   USE agrif_modutil
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, isize, jsize, id, ji, jj
26      INTEGER :: tabdim1, tabdim2, nxhr, nyhr, nxyhr
27      INTEGER, DIMENSION(2) :: ddims
28      INTEGER, DIMENSION(3) :: status
29      INTEGER, DIMENSION(:,:), ALLOCATABLE :: trouble_points ,  vardep,mask_oce
30      INTEGER :: iimin,iimax,jjmin,jjmax
31      INTEGER, DIMENSION(1) :: i_min,i_max
32      INTEGER, DIMENSION(1) ::j_min,j_max
33      REAL(wp), DIMENSION(jpi,jpj) :: myglamf
34      INTEGER,DIMENSION(:)  ,POINTER     ::   src_add,dst_add => NULL()
35      REAL(wp), DIMENSION(:)  ,ALLOCATABLE ::   vardep1d, lon_new1D,lat_new1D 
36      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bathy_new, lat_new, lon_new, bathy_test
37      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: coarselon, coarselat, coarsebathy
38      REAL(wp) :: Cell_lonmin, Cell_lonmax, Cell_latmin, Cell_latmax
39      LOGICAL :: identical_grids
40      LOGICAL, DIMENSION(:,:), ALLOCATABLE     ::   masksrc
41      REAL*8, DIMENSION(:,:),POINTER     ::   matrix,interpdata => NULL() 
42      LOGICAL :: lonlat_2D
43
44      CHARACTER(32) :: bathyfile, bathyname, lonname,latname       
45
46     lonlat_2D=.false.
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, lagrif=.FALSE. )
54      id = iom_varid( inum, bathyname, ddims )
55     
56      status=-1
57      ALLOCATE(lon_new  (ddims(1),ddims(2)), STAT=status(1)) 
58      ALLOCATE(lat_new  (ddims(1),ddims(2)), STAT=status(2)) 
59      ALLOCATE(bathy_new(ddims(1),ddims(2)), STAT=status(3)) 
60      IF( sum(status) /= 0 )   CALL ctl_stop( 'STOP', 'dom_bat : unable to allocate arrays' )
61
62      IF (lonlat_2D) THEN
63        CALL iom_get  ( inum, jpdom_unknown, lonname, lon_new )
64        CALL iom_get  ( inum, jpdom_unknown, latname, lat_new )
65      ELSE
66        ALLOCATE(lon_new1D(ddims(1)), lat_new1D(ddims(2)))
67        CALL iom_get  ( inum, jpdom_unknown, lonname, lon_new1D )
68        CALL iom_get  ( inum, jpdom_unknown, latname, lat_new1D )
69        DO ji=1, ddims(1)
70           lon_new(ji,:)=lon_new1D(ji)
71        ENDDO             
72        DO ji=1, ddims(2)
73           lat_new(:,ji)=lat_new1D(ji)
74        ENDDO             
75      ENDIF 
76      CALL iom_get  ( inum, jpdom_unknown, bathyname, bathy_new )
77      CALL iom_close (inum)
78      WHERE (bathy_new > 0.)       
79        bathy_new=0.
80      ENDWHERE 
81      bathy_new=-bathy_new 
82
83       ! Eventually add here a pre-selection of the area (coarselon/lat)
84       i_min=10  ; j_min=10
85       i_max= ddims(1)-10 ; j_max=ddims(2)-10 
86
87       tabdim1 = i_max(1) -  i_min(1) + 1 
88       tabdim2 = j_max(1) - j_min(1) + 1
89       !
90
91       ALLOCATE(coarselon(tabdim1,tabdim2))
92       ALLOCATE(coarselat(tabdim1,tabdim2))
93       ALLOCATE(coarsebathy(tabdim1,tabdim2))         
94       !
95       WHERE( lon_new < 0. )   lon_new = lon_new + 360.
96       myglamf=glamf
97       WHERE( myglamf < 0. )   myglamf = myglamf + 360.
98
99       coarselat(:,:)   =  lat_new  (i_min(1):i_max(1), j_min(1):j_max(1))
100       coarselon  (:,:) =  lon_new  (i_min(1):i_max(1), j_min(1):j_max(1))
101       coarsebathy(:,:) =  bathy_new(i_min(1):i_max(1), j_min(1):j_max(1))
102
103       IF( nn_interp == 0 .OR. nn_interp == 1 ) THEN   ! arithmetic or median averages
104        !                                                            ! -----------------------------
105        ALLOCATE(trouble_points(jpi,jpj))
106        trouble_points(:,:) = 0
107        !                       
108        DO jj = 2, jpj
109           DO ji = 2, jpi
110              !       
111              ! fine grid cell extension               
112              Cell_lonmin = MIN(myglamf(ji-1,jj-1),myglamf(ji,jj-1),myglamf(ji,jj),myglamf(ji-1,jj))
113              Cell_lonmax = MAX(myglamf(ji-1,jj-1),myglamf(ji,jj-1),myglamf(ji,jj),myglamf(ji-1,jj))
114              Cell_latmin = MIN(gphif(ji-1,jj-1),gphif(ji,jj-1),gphif(ji,jj),gphif(ji-1,jj))
115              Cell_latmax = MAX(gphif(ji-1,jj-1),gphif(ji,jj-1),gphif(ji,jj),gphif(ji-1,jj)) 
116              !               
117              ! look for points in G0 (bathy dataset) contained in the fine grid cells 
118              iimin = 1
119              DO WHILE( coarselon(iimin,1) < Cell_lonmin ) 
120                 iimin = iimin + 1
121              ENDDO
122              !               
123              jjmin = 1
124              DO WHILE( coarselat(iimin,jjmin) < Cell_latmin ) 
125                 jjmin = jjmin + 1
126              ENDDO
127             !               
128              iimax = iimin 
129              DO WHILE( coarselon(iimax,1) <= Cell_lonmax ) 
130                 iimax = iimax + 1
131                 iimax = MIN( iimax,SIZE(coarsebathy,1))
132              ENDDO
133              !                               
134              jjmax = jjmin 
135              DO WHILE( coarselat(iimax,jjmax) <= Cell_latmax ) 
136                 jjmax = jjmax + 1
137                 jjmax = MIN( jjmax,SIZE(coarsebathy,2))
138              ENDDO
139              !
140              IF( .NOT. Agrif_Root() ) THEN
141                 iimax = iimax-1
142                 jjmax = jjmax-1
143              ELSE
144                 iimax = MAX(iimin,iimax-1)
145                 jjmax = MAX(jjmin,jjmax-1)
146              ENDIF
147              !               
148              iimin = MAX( iimin,1 )
149              jjmin = MAX( jjmin,1 )
150              iimax = MIN( iimax,SIZE(coarsebathy,1))
151              jjmax = MIN( jjmax,SIZE(coarsebathy,2))
152
153              nxhr = iimax - iimin + 1
154              nyhr = jjmax - jjmin + 1                   
155
156               
157              IF( nxhr == 0 .OR. nyhr == 0 ) THEN
158                 !
159                 trouble_points(ji,jj) = 1
160                 !
161              ELSE
162                 !
163                 ALLOCATE( vardep(nxhr,nyhr), mask_oce(nxhr,nyhr) )
164                 vardep(:,:) = coarsebathy(iimin:iimax,jjmin:jjmax)
165                 !
166                 WHERE( vardep(:,:) .GT. 0. )   ;   mask_oce = 1 ;
167                 ELSEWHERE                      ;   mask_oce = 0 ;
168                 ENDWHERE
169                 !
170                 nxyhr = nxhr*nyhr
171                 IF( SUM(mask_oce) < 0.5*(nxyhr) ) THEN ! if more than half of the points are on land then bathy fine = 0
172                    bathy(ji,jj) = 0.
173                 ELSE
174                    IF( nn_interp == 0 ) THEN     ! Arithmetic average
175                      bathy(ji,jj) = SUM( vardep(:,:) * mask_oce(:,:) ) / SUM( mask_oce(:,:) )
176                    ELSE                                  ! Median average
177                       ALLOCATE(vardep1d(nxyhr))
178                       vardep1d = RESHAPE(vardep,(/ nxyhr /) )
179                       !!CALL ssort(vardep1d,nxyhr)
180                       CALL quicksort(vardep1d,1,nxyhr)
181                       !
182                       ! Calculate median
183                       IF (MOD(nxyhr,2) .NE. 0) THEN
184                          bathy(ji,jj) = vardep1d( nxyhr/2 + 1 )
185                       ELSE
186                          bathy(ji,jj) = 0.5 * ( vardep1d(nxyhr/2) + vardep1d(nxyhr/2+1) )
187                       END IF
188                       DEALLOCATE(vardep1d)   
189                    ENDIF
190                 ENDIF
191                 DEALLOCATE (mask_oce,vardep)
192                 !
193              ENDIF
194           ENDDO
195        ENDDO
196
197        IF( SUM( trouble_points ) > 0 ) THEN
198           PRINT*,'too much empty cells, proceed to bilinear interpolation'
199           nn_interp = 2
200           stop
201        ENDIF           
202     ENDIF
203
204#undef MYTEST
205#ifdef MYTEST
206         IF( nn_interp == 2) THEN                        ! Bilinear interpolation
207        !                                                    ! -----------------------------
208        identical_grids = .FALSE.
209
210        IF( SIZE(coarselat,1) == jpi .AND. SIZE(coarselat,2) == jpj  .AND.   &
211            SIZE(coarselon,1) == jpj .AND. SIZE(coarselon,2) == jpj ) THEN
212           IF( MAXVAL( ABS(coarselat(:,:)- gphit(:,:)) ) < 0.0001 .AND.   &
213               MAXVAL( ABS(coarselon(:,:)- glamt(:,:)) ) < 0.0001 ) THEN
214              PRINT*,'same grid between ', cn_topo,' and child domain'   
215              bathy = bathy_new
216              identical_grids = .TRUE.                         
217           ENDIF
218        ENDIF
219
220        IF( .NOT. identical_grids ) THEN
221
222           ALLOCATE(masksrc(SIZE(bathy_new,1),SIZE(bathy_new,2)))
223           ALLOCATE(bathy_test(jpi,jpj))
224           !
225           !Where(G0%bathy_meter.le.0.00001)
226           !  masksrc = .false.
227           !ElseWhere
228              masksrc = .TRUE.
229           !End where                       
230           !           
231           ! compute remapping matrix thanks to SCRIP package           
232           CALL get_remap_matrix(coarselat,gphit,coarselon,glamt,masksrc,matrix,src_add,dst_add)
233           CALL make_remap(bathy_new,bathy_test,jpi,jpj,matrix,src_add,dst_add) 
234           !                                 
235           bathy = bathy_test               
236           !           
237           DEALLOCATE(masksrc)
238           DEALLOCATE(bathy_test) 
239
240        ENDIF
241        !           
242     ENDIF
243#endif
244   END SUBROUTINE dom_bat
245
246END MODULE dombat
Note: See TracBrowser for help on using the repository browser.