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