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 NEMO/branches/2019/ENHANCE-03_domcfg/src – NEMO

source: NEMO/branches/2019/ENHANCE-03_domcfg/src/dombat.F90 @ 11129

Last change on this file since 11129 was 11129, checked in by mathiot, 5 years ago

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