Changeset 12253 for utils/tools/NESTING/src
 Timestamp:
 20191216T10:59:22+01:00 (12 months ago)
 Location:
 utils/tools/NESTING/src
 Files:

 3 edited
Legend:
 Unmodified
 Added
 Removed

utils/tools/NESTING/src/agrif_connect_topo.f90
r10025 r12253 255 255 ! 256 256 IF (Grid%bathy_meter(i,j) .EQ. 0.0 ) THEN 257 Grid%bathy_level(i,j)=0 .257 Grid%bathy_level(i,j)=0 258 258 ELSE 259 259 ! 
utils/tools/NESTING/src/agrif_create_bathy.f90
r10381 r12253 87 87 IF( .NOT.new_topo .AND. .NOT.partial_steps ) THEN ! **** First option : no new topo file & no partial steps 88 88 ! ! ==> interpolate levels directly from parent file 89 ! ! using bilinear interpolation 89 90 ! !  90 91 WRITE(*,*) '*** First option: no new topo file & no partial steps' … … 169 170 DO jj=1,nyfin 170 171 DO ji=1,nxfin 171 IF (G1%bathy_level(ji,jj).LT.0 .) THEN172 IF (G1%bathy_level(ji,jj).LT.0) THEN 172 173 PRINT *,'error in ',ji,jj,G1%bathy_level(ji,jj) 173 174 STOP … … 175 176 ENDDO 176 177 ENDDO 177 WHERE ((G1%bathy_level.LT.3 .).AND.(G1%bathy_level.GT.0.)) G1%bathy_level=3178 WHERE ((G1%bathy_level.LT.3).AND.(G1%bathy_level.GT.0)) G1%bathy_level=3 178 179 ! 179 180 ! remove closed seas … … 181 182 ALLOCATE(bathy_test(nxfin,nyfin)) 182 183 183 bathy_test=0 .184 WHERE (G1%bathy_level(1,:) .GT.0 .) bathy_test(1,:)=1185 WHERE (G1%bathy_level(nxfin,:) .GT.0 .) bathy_test(nxfin,:)=1186 WHERE (G1%bathy_level(:,1) .GT.0 .) bathy_test(:,1)=1187 WHERE (G1%bathy_level(:,nyfin) .GT.0 .) bathy_test(:,nyfin)=1184 bathy_test=0 185 WHERE (G1%bathy_level(1,:) .GT.0) bathy_test(1,:)=1 186 WHERE (G1%bathy_level(nxfin,:) .GT.0) bathy_test(nxfin,:)=1 187 WHERE (G1%bathy_level(:,1) .GT.0) bathy_test(:,1)=1 188 WHERE (G1%bathy_level(:,nyfin) .GT.0) bathy_test(:,nyfin)=1 188 189 189 190 nbadd = 1 … … 192 193 DO jj=2,nyfin1 193 194 DO ji=2,nxfin1 194 IF (G1%bathy_level(ji,jj).GT.0 .) THEN195 IF (G1%bathy_level(ji,jj).GT.0) THEN 195 196 IF (MAX(bathy_test(ji,jj+1),bathy_test(ji,jj1),bathy_test(ji1,jj),bathy_test(ji+1,jj)).EQ.1) THEN 196 IF (bathy_test(ji,jj).NE.1 .) nbadd = nbadd + 1197 bathy_test(ji,jj)=1 .197 IF (bathy_test(ji,jj).NE.1) nbadd = nbadd + 1 198 bathy_test(ji,jj)=1 198 199 ENDIF 199 200 ENDIF … … 202 203 ENDDO 203 204 204 WHERE (bathy_test.EQ.0.) G1%bathy_level = 0 .205 WHERE (bathy_test.EQ.0.) G1%bathy_level = 0 205 206 206 207 DEALLOCATE(bathy_test) … … 214 215 !!status = write_domcfg(TRIM(parent_domcfg_out),G0) ! do not activate it 215 216 216 WRITE(*,*) '****** Bathymetry successfully created if no new topo ******'217 WRITE(*,*) '****** Bathymetry successfully created if no new topo and no partial steps ******' 217 218 STOP 218 219 ! … … 277 278 ! ==> It gives G1%bathy_meter from G0%bathy_meter 278 279 !  279 ! !  280 IF( type_bathy_interp == 0 .OR. type_bathy_interp == 1 ) THEN ! arithmetic or median averages 281 ! !  282 ALLOCATE(trouble_points(nxfin,nyfin)) 283 trouble_points(:,:) = 0 284 ! 285 DO jj = 2, nyfin 286 DO ji = 2, nxfin 287 ! 288 ! fine grid cell extension 289 Cell_lonmin = MIN(G1%glamf(ji1,jj1),G1%glamf(ji,jj1),G1%glamf(ji,jj),G1%glamf(ji1,jj)) 290 Cell_lonmax = MAX(G1%glamf(ji1,jj1),G1%glamf(ji,jj1),G1%glamf(ji,jj),G1%glamf(ji1,jj)) 291 Cell_latmin = MIN(G1%gphif(ji1,jj1),G1%gphif(ji,jj1),G1%gphif(ji,jj),G1%gphif(ji1,jj)) 292 Cell_latmax = MAX(G1%gphif(ji1,jj1),G1%gphif(ji,jj1),G1%gphif(ji,jj),G1%gphif(ji1,jj)) 293 ! 294 ! look for points in G0 (bathy dataset) contained in the fine grid cells 295 iimin = 1 296 DO WHILE( G0%nav_lon(iimin,1) < Cell_lonmin ) 297 iimin = iimin + 1 280 identical_grids = .FALSE. 281 282 IF( SIZE(G0%nav_lat,1) == SIZE(G1%nav_lat,1) .AND. SIZE(G0%nav_lat,2) == SIZE(G1%nav_lat,2) .AND. & 283 & SIZE(G0%nav_lon,1) == SIZE(G1%nav_lon,1) .AND. SIZE(G0%nav_lon,2) == SIZE(G1%nav_lon,2) ) THEN 284 IF( MAXVAL( ABS(G0%nav_lat(:,:) G1%nav_lat(:,:)) ) < 0.0001 .AND. & 285 & MAXVAL( ABS(G0%nav_lon(:,:) G1%nav_lon(:,:)) ) < 0.0001 ) THEN 286 WRITE(*,*) '' 287 WRITE(*,*) 'same grid between parent and child domains => NO INTERPOLATION' 288 WRITE(*,*) '' 289 G1%bathy_meter = G0%bathy_meter 290 identical_grids = .TRUE. 291 ENDIF 292 ENDIF 293 294 IF( .NOT. identical_grids ) THEN 295 ! !  296 IF( type_bathy_interp == 0 .OR. type_bathy_interp == 1 ) THEN ! arithmetic or median averages 297 ! !  298 ALLOCATE(trouble_points(nxfin,nyfin)) 299 trouble_points(:,:) = 0 300 ! 301 DO jj = 2, nyfin 302 DO ji = 2, nxfin 303 ! 304 ! fine grid cell extension 305 Cell_lonmin = MIN(G1%glamf(ji1,jj1),G1%glamf(ji,jj1),G1%glamf(ji,jj),G1%glamf(ji1,jj)) 306 Cell_lonmax = MAX(G1%glamf(ji1,jj1),G1%glamf(ji,jj1),G1%glamf(ji,jj),G1%glamf(ji1,jj)) 307 Cell_latmin = MIN(G1%gphif(ji1,jj1),G1%gphif(ji,jj1),G1%gphif(ji,jj),G1%gphif(ji1,jj)) 308 Cell_latmax = MAX(G1%gphif(ji1,jj1),G1%gphif(ji,jj1),G1%gphif(ji,jj),G1%gphif(ji1,jj)) 309 ! 310 ! look for points in G0 (bathy dataset) contained in the fine grid cells 311 iimin = 1 312 DO WHILE( G0%nav_lon(iimin,1) < Cell_lonmin ) 313 iimin = iimin + 1 314 ENDDO 315 ! 316 jjmin = 1 317 DO WHILE( G0%nav_lat(iimin,jjmin) < Cell_latmin ) 318 jjmin = jjmin + 1 319 ENDDO 320 ! 321 iimax = iimin 322 DO WHILE( G0%nav_lon(iimax,1) <= Cell_lonmax ) 323 iimax = iimax + 1 324 iimax = MIN( iimax,SIZE(G0%bathy_meter,1)) 325 ENDDO 326 ! 327 jjmax = jjmin 328 DO WHILE( G0%nav_lat(iimax,jjmax) <= Cell_latmax ) 329 jjmax = jjmax + 1 330 jjmax = MIN( jjmax,SIZE(G0%bathy_meter,2)) 331 ENDDO 332 ! 333 IF( ln_agrif_domain ) THEN 334 iimax = iimax1 335 jjmax = jjmax1 336 ELSE 337 iimax = MAX(iimin,iimax1) 338 jjmax = MAX(jjmin,jjmax1) 339 ENDIF 340 ! 341 iimin = MAX( iimin,1 ) 342 jjmin = MAX( jjmin,1 ) 343 iimax = MIN( iimax,SIZE(G0%bathy_meter,1)) 344 jjmax = MIN( jjmax,SIZE(G0%bathy_meter,2)) 345 346 nxhr = iimax  iimin + 1 347 nyhr = jjmax  jjmin + 1 348 349 IF( nxhr == 0 .OR. nyhr == 0 ) THEN 350 ! 351 trouble_points(ji,jj) = 1 352 ! 353 ELSE 354 ! 355 ALLOCATE( vardep(nxhr,nyhr), mask_oce(nxhr,nyhr) ) 356 vardep(:,:) = G0%bathy_meter(iimin:iimax,jjmin:jjmax) 357 ! 358 WHERE( vardep(:,:) .GT. 0. ) ; mask_oce = 1 ; 359 ELSEWHERE ; mask_oce = 0 ; 360 ENDWHERE 361 ! 362 nxyhr = nxhr*nyhr 363 IF( SUM(mask_oce) < 0.5*(nxyhr) ) THEN ! if more than half of the points are on land then bathy fine = 0 364 G1%bathy_meter(ji,jj) = 0. 365 ELSE 366 IF( type_bathy_interp == 0 ) THEN ! Arithmetic average 367 G1%bathy_meter(ji,jj) = SUM( vardep(:,:) * mask_oce(:,:) ) / SUM( mask_oce(:,:) ) 368 ELSE ! Median average 369 ALLOCATE(vardep1d(nxyhr)) 370 vardep1d = RESHAPE(vardep,(/ nxyhr /) ) 371 !!CALL ssort(vardep1d,nxyhr) 372 CALL quicksort(vardep1d,1,nxyhr) 373 ! 374 ! Calculate median 375 IF (MOD(nxyhr,2) .NE. 0) THEN 376 G1%bathy_meter(ji,jj) = vardep1d( nxyhr/2 + 1 ) 377 ELSE 378 G1%bathy_meter(ji,jj) = 0.5 * ( vardep1d(nxyhr/2) + vardep1d(nxyhr/2+1) ) 379 END IF 380 DEALLOCATE(vardep1d) 381 ENDIF 382 ENDIF 383 DEALLOCATE (mask_oce,vardep) 384 ! 385 ENDIF 298 386 ENDDO 299 !300 jjmin = 1301 DO WHILE( G0%nav_lat(iimin,jjmin) < Cell_latmin )302 jjmin = jjmin + 1303 ENDDO304 !305 iimax = iimin306 DO WHILE( G0%nav_lon(iimax,1) <= Cell_lonmax )307 iimax = iimax + 1308 iimax = MIN( iimax,SIZE(G0%bathy_meter,1))309 ENDDO310 !311 jjmax = jjmin312 DO WHILE( G0%nav_lat(iimax,jjmax) <= Cell_latmax )313 jjmax = jjmax + 1314 jjmax = MIN( jjmax,SIZE(G0%bathy_meter,2))315 ENDDO316 !317 IF( ln_agrif_domain ) THEN318 iimax = iimax1319 jjmax = jjmax1320 ELSE321 iimax = MAX(iimin,iimax1)322 jjmax = MAX(jjmin,jjmax1)323 ENDIF324 !325 iimin = MAX( iimin,1 )326 jjmin = MAX( jjmin,1 )327 iimax = MIN( iimax,SIZE(G0%bathy_meter,1))328 jjmax = MIN( jjmax,SIZE(G0%bathy_meter,2))329 330 nxhr = iimax  iimin + 1331 nyhr = jjmax  jjmin + 1332 333 IF( nxhr == 0 .OR. nyhr == 0 ) THEN334 !335 trouble_points(ji,jj) = 1336 !337 ELSE338 !339 ALLOCATE( vardep(nxhr,nyhr), mask_oce(nxhr,nyhr) )340 vardep(:,:) = G0%bathy_meter(iimin:iimax,jjmin:jjmax)341 !342 WHERE( vardep(:,:) .GT. 0. ) ; mask_oce = 1 ;343 ELSEWHERE ; mask_oce = 0 ;344 ENDWHERE345 !346 nxyhr = nxhr*nyhr347 IF( SUM(mask_oce) < 0.5*(nxyhr) ) THEN ! if more than half of the points are on land then bathy fine = 0348 G1%bathy_meter(ji,jj) = 0.349 ELSE350 IF( type_bathy_interp == 0 ) THEN ! Arithmetic average351 G1%bathy_meter(ji,jj) = SUM( vardep(:,:) * mask_oce(:,:) ) / SUM( mask_oce(:,:) )352 ELSE ! Median average353 ALLOCATE(vardep1d(nxyhr))354 vardep1d = RESHAPE(vardep,(/ nxyhr /) )355 !!CALL ssort(vardep1d,nxyhr)356 CALL quicksort(vardep1d,1,nxyhr)357 !358 ! Calculate median359 IF (MOD(nxyhr,2) .NE. 0) THEN360 G1%bathy_meter(ji,jj) = vardep1d( nxyhr/2 + 1 )361 ELSE362 G1%bathy_meter(ji,jj) = 0.5 * ( vardep1d(nxyhr/2) + vardep1d(nxyhr/2+1) )363 END IF364 DEALLOCATE(vardep1d)365 ENDIF366 ENDIF367 DEALLOCATE (mask_oce,vardep)368 !369 ENDIF370 387 ENDDO 371 ENDDO 372 373 IF( SUM( trouble_points ) > 0 ) THEN 374 PRINT*,'too much empty cells, proceed to bilinear interpolation' 375 type_bathy_interp = 2 376 ENDIF 377 378 ENDIF 379 ! !  380 IF( type_bathy_interp == 2) THEN ! Bilinear interpolation 381 ! !  382 identical_grids = .FALSE. 383 384 IF( SIZE(G0%nav_lat,1) == SIZE(G1%nav_lat,1) .AND. SIZE(G0%nav_lat,2) == SIZE(G1%nav_lat,2) .AND. & 385 SIZE(G0%nav_lon,1) == SIZE(G1%nav_lon,1) .AND. SIZE(G0%nav_lon,2) == SIZE(G1%nav_lon,2) ) THEN 386 IF( MAXVAL( ABS(G0%nav_lat(:,:) G1%nav_lat(:,:)) ) < 0.0001 .AND. & 387 MAXVAL( ABS(G0%nav_lon(:,:) G1%nav_lon(:,:)) ) < 0.0001 ) THEN 388 PRINT*,'same grid between ',elevation_database,' and child domain' 389 G1%bathy_meter = G0%bathy_meter 390 identical_grids = .TRUE. 388 389 IF( SUM( trouble_points ) > 0 ) THEN 390 PRINT*,'too much empty cells, proceed to bilinear interpolation' 391 type_bathy_interp = 2 391 392 ENDIF 392 ENDIF 393 394 IF( .NOT. identical_grids ) THEN 393 394 ENDIF 395 ! !  396 IF( type_bathy_interp == 2) THEN ! Bilinear interpolation 397 ! !  395 398 396 399 ALLOCATE(masksrc(SIZE(G0%bathy_meter,1),SIZE(G0%bathy_meter,2))) … … 400 403 ! masksrc = .false. 401 404 !ElseWhere 402 405 masksrc = .TRUE. 403 406 !End where 404 407 ! … … 411 414 DEALLOCATE(masksrc) 412 415 DEALLOCATE(bathy_test) 413 416 414 417 ENDIF 415 418 ! 416 ENDIF 419 ENDIF ! not identical grids 417 420 !  418 421 ! At this stage bathymetry in meters has already been interpolated on fine grid 
utils/tools/NESTING/src/agrif_readwrite.f90
r11231 r12253 241 241 CALL write_ncdf_var('nav_lon' ,dimnames,name,Grid%nav_lon ,'float') 242 242 CALL write_ncdf_var('nav_lat' ,dimnames,name,Grid%nav_lat ,'float') 243 CALL write_ncdf_var(parent_level_name,dimnames,name, Grid%bathy_level,'float')243 CALL write_ncdf_var(parent_level_name,dimnames,name,NINT(Grid%bathy_level),'integer') 244 244 ! 245 245 CALL copy_ncdf_att('nav_lon' ,TRIM(parent_bathy_level),name,MINVAL(Grid%nav_lon),MAXVAL(Grid%nav_lon)) … … 548 548 ln_zps = 1 549 549 ln_zco = 0 550 CALL zgr_zps( nx, ny, gdept_1d, gdepw_1d, e3t_1d, e3w_1d, zbathy, & ! input551 & mbathy, e3t_0, e3u_0, e3v_0, e3f_0, e3w_0, e3uw_0, e3vw_0 ) ! output552 550 ELSE 553 551 ln_zps = 0 554 552 ln_zco = 1 555 553 ENDIF 554 555 CALL zgr_zps( nx, ny, gdept_1d, gdepw_1d, e3t_1d, e3w_1d, zbathy, & ! input 556 & mbathy, e3t_0, e3u_0, e3v_0, e3f_0, e3w_0, e3uw_0, e3vw_0 ) ! output 556 557 557 558 ! top/bottom levels … … 869 870 END WHERE 870 871 871 ! Compute mbathy for ocean points (i.e. the number of ocean levels) 872 ! find the number of ocean levels such that the last level thickness 873 ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where 874 ! e3t_1d is the reference level thickness 875 DO jk = N1, 1, 1 876 zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 877 WHERE( 0. < zbathy(:,:) .AND. zbathy(:,:) <= zdepth ) mbathy(:,:) = jk1 878 END DO 879 872 IF( partial_steps ) THEN 873 ! Compute mbathy for ocean points (i.e. the number of ocean levels) 874 ! find the number of ocean levels such that the last level thickness 875 ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where 876 ! e3t_1d is the reference level thickness 877 DO jk = N1, 1, 1 878 zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 879 WHERE( 0. < zbathy(:,:) .AND. zbathy(:,:) <= zdepth ) mbathy(:,:) = jk1 880 END DO 881 ELSE 882 DO jk = 1, N 883 WHERE( 0. < zbathy(:,:) .AND. zbathy(:,:) >= gdepw_1d(jk) ) mbathy(:,:) = jk1 884 END DO 885 ENDIF 886 880 887 ! Scale factors and depth at T and Wpoints 881 888 DO jk = 1, N ! intitialization to the reference zcoordinate … … 890 897 !! IF ( ln_isfcav == 1 ) CALL zgr_isf 891 898 ! 892 ! Scale factors and depth at T and Wpoints 893 ! IF ( ln_isfcav == 0 ) THEN 899 IF( partial_steps ) THEN 900 ! Scale factors and depth at T and Wpoints 901 ! IF ( ln_isfcav == 0 ) THEN 894 902 DO jj = 1, ny 895 903 DO ji = 1, nx … … 913 921 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 914 922 ENDIF 915 !gm Bug? check the gdepw_1d923 !gm Bug? check the gdepw_1d 916 924 ! ... on ik 917 925 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1)  gdepw_1d(ik) ) & … … 931 939 END DO 932 940 ! 933 ! DO jj = 1, ny 934 ! DO ji = 1, nx 935 ! ik = mbathy(ji,jj) 936 ! IF( ik > 0 ) THEN ! ocean point only 937 ! e3tp (ji,jj) = e3t_0(ji,jj,ik) 938 ! e3wp (ji,jj) = e3w_0(ji,jj,ik) 939 ! ENDIF 940 ! END DO 941 ! END DO 942 ! END IF 941 ! DO jj = 1, ny 942 ! DO ji = 1, nx 943 ! ik = mbathy(ji,jj) 944 ! IF( ik > 0 ) THEN ! ocean point only 945 ! e3tp (ji,jj) = e3t_0(ji,jj,ik) 946 ! e3wp (ji,jj) = e3w_0(ji,jj,ik) 947 ! ENDIF 948 ! END DO 949 ! END DO 950 ! END IF 951 ENDIF 943 952 ! 944 953 ! Scale factors and depth at U, V, UW and VWpoints
Note: See TracChangeset
for help on using the changeset viewer.