Ignore:
Timestamp:
2019-12-16T10:59:22+01:00 (12 months ago)
Author:
clem
Message:

1) resolve some conflicts when using partial steps or not. 2) Make sure grids remain identical when they should. 3) in case of a double zoom, it is no more necessary to have the bilinear interpolation option when updating the 1st parent grid. Btw, I know these comments are unclear but it is very difficult to explain…

Location:
utils/tools/NESTING/src
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • utils/tools/NESTING/src/agrif_connect_topo.f90

    r10025 r12253  
    255255          ! 
    256256          IF (Grid%bathy_meter(i,j) .EQ. 0.0 ) THEN 
    257              Grid%bathy_level(i,j)=0. 
     257             Grid%bathy_level(i,j)=0 
    258258          ELSE 
    259259             !  
  • utils/tools/NESTING/src/agrif_create_bathy.f90

    r10381 r12253  
    8787  IF( .NOT.new_topo .AND. .NOT.partial_steps ) THEN   ! **** First option : no new topo file & no partial steps 
    8888     !                                                !                 ==> interpolate levels directly from parent file 
     89     !                                                !                     using bilinear interpolation 
    8990     !                                                ! ------------------------------------------------------------------ 
    9091     WRITE(*,*) '*** First option: no new topo file & no partial steps' 
     
    169170     DO jj=1,nyfin 
    170171        DO ji=1,nxfin 
    171            IF (G1%bathy_level(ji,jj).LT.0.) THEN 
     172           IF (G1%bathy_level(ji,jj).LT.0) THEN 
    172173              PRINT *,'error in ',ji,jj,G1%bathy_level(ji,jj) 
    173174              STOP 
     
    175176        ENDDO 
    176177     ENDDO 
    177      WHERE ((G1%bathy_level.LT.3.).AND.(G1%bathy_level.GT.0.))   G1%bathy_level=3 
     178     WHERE ((G1%bathy_level.LT.3).AND.(G1%bathy_level.GT.0))   G1%bathy_level=3 
    178179     ! 
    179180     ! remove closed seas 
     
    181182        ALLOCATE(bathy_test(nxfin,nyfin)) 
    182183 
    183         bathy_test=0. 
    184         WHERE (G1%bathy_level(1,:)     .GT.0.)   bathy_test(1,:)=1 
    185         WHERE (G1%bathy_level(nxfin,:) .GT.0.)   bathy_test(nxfin,:)=1 
    186         WHERE (G1%bathy_level(:,1)     .GT.0.)   bathy_test(:,1)=1 
    187         WHERE (G1%bathy_level(:,nyfin) .GT.0.)   bathy_test(:,nyfin)=1 
     184        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 
    188189 
    189190        nbadd = 1 
     
    192193           DO jj=2,nyfin-1 
    193194              DO ji=2,nxfin-1 
    194                  IF (G1%bathy_level(ji,jj).GT.0.) THEN 
     195                 IF (G1%bathy_level(ji,jj).GT.0) THEN 
    195196                    IF (MAX(bathy_test(ji,jj+1),bathy_test(ji,jj-1),bathy_test(ji-1,jj),bathy_test(ji+1,jj)).EQ.1) THEN 
    196                        IF (bathy_test(ji,jj).NE.1.) nbadd = nbadd + 1 
    197                        bathy_test(ji,jj)=1. 
     197                       IF (bathy_test(ji,jj).NE.1) nbadd = nbadd + 1 
     198                       bathy_test(ji,jj)=1 
    198199                    ENDIF 
    199200                 ENDIF 
     
    202203        ENDDO 
    203204 
    204         WHERE (bathy_test.EQ.0.)   G1%bathy_level = 0. 
     205        WHERE (bathy_test.EQ.0.)   G1%bathy_level = 0 
    205206 
    206207        DEALLOCATE(bathy_test) 
     
    214215     !!status = write_domcfg(TRIM(parent_domcfg_out),G0) ! do not activate it 
    215216 
    216      WRITE(*,*) '****** Bathymetry successfully created if no new topo ******' 
     217     WRITE(*,*) '****** Bathymetry successfully created if no new topo and no partial steps ******' 
    217218     STOP 
    218219     ! 
     
    277278     ! ==> It gives G1%bathy_meter from G0%bathy_meter 
    278279     ! --------------------------------------------------------------------------------- 
    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(ji-1,jj-1),G1%glamf(ji,jj-1),G1%glamf(ji,jj),G1%glamf(ji-1,jj)) 
    290               Cell_lonmax = MAX(G1%glamf(ji-1,jj-1),G1%glamf(ji,jj-1),G1%glamf(ji,jj),G1%glamf(ji-1,jj)) 
    291               Cell_latmin = MIN(G1%gphif(ji-1,jj-1),G1%gphif(ji,jj-1),G1%gphif(ji,jj),G1%gphif(ji-1,jj)) 
    292               Cell_latmax = MAX(G1%gphif(ji-1,jj-1),G1%gphif(ji,jj-1),G1%gphif(ji,jj),G1%gphif(ji-1,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(ji-1,jj-1),G1%glamf(ji,jj-1),G1%glamf(ji,jj),G1%glamf(ji-1,jj)) 
     306                 Cell_lonmax = MAX(G1%glamf(ji-1,jj-1),G1%glamf(ji,jj-1),G1%glamf(ji,jj),G1%glamf(ji-1,jj)) 
     307                 Cell_latmin = MIN(G1%gphif(ji-1,jj-1),G1%gphif(ji,jj-1),G1%gphif(ji,jj),G1%gphif(ji-1,jj)) 
     308                 Cell_latmax = MAX(G1%gphif(ji-1,jj-1),G1%gphif(ji,jj-1),G1%gphif(ji,jj),G1%gphif(ji-1,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 = iimax-1 
     335                    jjmax = jjmax-1 
     336                 ELSE 
     337                    iimax = MAX(iimin,iimax-1) 
     338                    jjmax = MAX(jjmin,jjmax-1) 
     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 
    298386              ENDDO 
    299               !                
    300               jjmin = 1 
    301               DO WHILE( G0%nav_lat(iimin,jjmin) < Cell_latmin )  
    302                  jjmin = jjmin + 1 
    303               ENDDO 
    304               !                 
    305               iimax = iimin  
    306               DO WHILE( G0%nav_lon(iimax,1) <= Cell_lonmax )  
    307                  iimax = iimax + 1 
    308                  iimax = MIN( iimax,SIZE(G0%bathy_meter,1)) 
    309               ENDDO 
    310               !                                
    311               jjmax = jjmin  
    312               DO WHILE( G0%nav_lat(iimax,jjmax) <= Cell_latmax )  
    313                  jjmax = jjmax + 1 
    314                  jjmax = MIN( jjmax,SIZE(G0%bathy_meter,2)) 
    315               ENDDO 
    316               ! 
    317               IF( ln_agrif_domain ) THEN 
    318                  iimax = iimax-1 
    319                  jjmax = jjmax-1 
    320               ELSE 
    321                  iimax = MAX(iimin,iimax-1) 
    322                  jjmax = MAX(jjmin,jjmax-1) 
    323               ENDIF 
    324               !                
    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 + 1 
    331               nyhr = jjmax - jjmin + 1                     
    332  
    333               IF( nxhr == 0 .OR. nyhr == 0 ) THEN 
    334                  ! 
    335                  trouble_points(ji,jj) = 1 
    336                  ! 
    337               ELSE 
    338                  ! 
    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                  ENDWHERE 
    345                  ! 
    346                  nxyhr = nxhr*nyhr 
    347                  IF( SUM(mask_oce) < 0.5*(nxyhr) ) THEN ! if more than half of the points are on land then bathy fine = 0 
    348                     G1%bathy_meter(ji,jj) = 0. 
    349                  ELSE 
    350                     IF( type_bathy_interp == 0 ) THEN     ! Arithmetic average 
    351                        G1%bathy_meter(ji,jj) = SUM( vardep(:,:) * mask_oce(:,:) ) / SUM( mask_oce(:,:) ) 
    352                     ELSE                                  ! Median average 
    353                        ALLOCATE(vardep1d(nxyhr)) 
    354                        vardep1d = RESHAPE(vardep,(/ nxyhr /) ) 
    355                        !!CALL ssort(vardep1d,nxyhr) 
    356                        CALL quicksort(vardep1d,1,nxyhr) 
    357                        ! 
    358                        ! Calculate median 
    359                        IF (MOD(nxyhr,2) .NE. 0) THEN 
    360                           G1%bathy_meter(ji,jj) = vardep1d( nxyhr/2 + 1 ) 
    361                        ELSE 
    362                           G1%bathy_meter(ji,jj) = 0.5 * ( vardep1d(nxyhr/2) + vardep1d(nxyhr/2+1) ) 
    363                        END IF 
    364                        DEALLOCATE(vardep1d)    
    365                     ENDIF 
    366                  ENDIF 
    367                  DEALLOCATE (mask_oce,vardep) 
    368                  ! 
    369               ENDIF 
    370387           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 
    391392           ENDIF 
    392         ENDIF 
    393  
    394         IF( .NOT. identical_grids ) THEN  
     393 
     394        ENDIF 
     395        !                                                       ! -----------------------------  
     396        IF( type_bathy_interp == 2) THEN                        ! Bilinear interpolation 
     397           !                                                    ! -----------------------------  
    395398 
    396399           ALLOCATE(masksrc(SIZE(G0%bathy_meter,1),SIZE(G0%bathy_meter,2))) 
     
    400403           !  masksrc = .false. 
    401404           !ElseWhere 
    402               masksrc = .TRUE. 
     405           masksrc = .TRUE. 
    403406           !End where                        
    404407           !             
     
    411414           DEALLOCATE(masksrc) 
    412415           DEALLOCATE(bathy_test)  
    413  
     416            
    414417        ENDIF 
    415418        !             
    416      ENDIF 
     419     ENDIF ! not identical grids 
    417420     ! --- 
    418421     ! At this stage bathymetry in meters has already been interpolated on fine grid 
  • utils/tools/NESTING/src/agrif_readwrite.f90

    r11231 r12253  
    241241    CALL write_ncdf_var('nav_lon'        ,dimnames,name,Grid%nav_lon    ,'float') 
    242242    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') 
    244244    ! 
    245245    CALL copy_ncdf_att('nav_lon'        ,TRIM(parent_bathy_level),name,MINVAL(Grid%nav_lon),MAXVAL(Grid%nav_lon)) 
     
    548548        ln_zps = 1 
    549549        ln_zco = 0 
    550         CALL zgr_zps( nx, ny, gdept_1d, gdepw_1d, e3t_1d, e3w_1d, zbathy, &        ! input 
    551         &             mbathy, e3t_0, e3u_0, e3v_0, e3f_0, e3w_0, e3uw_0, e3vw_0 )  ! output 
    552550     ELSE 
    553551        ln_zps = 0 
    554552        ln_zco = 1 
    555553     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 
    556557 
    557558     ! top/bottom levels 
     
    869870      END WHERE 
    870871 
    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 = N-1, 1, -1 
    876          zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
    877          WHERE( 0. < zbathy(:,:) .AND. zbathy(:,:) <= zdepth )   mbathy(:,:) = jk-1 
    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 = N-1, 1, -1 
     878            zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
     879            WHERE( 0. < zbathy(:,:) .AND. zbathy(:,:) <= zdepth )   mbathy(:,:) = jk-1 
     880         END DO 
     881      ELSE 
     882         DO jk = 1, N 
     883            WHERE( 0. < zbathy(:,:) .AND. zbathy(:,:) >= gdepw_1d(jk) )   mbathy(:,:) = jk-1 
     884         END DO         
     885      ENDIF 
     886       
    880887      ! Scale factors and depth at T- and W-points 
    881888      DO jk = 1, N                        ! intitialization to the reference z-coordinate 
     
    890897      !!      IF ( ln_isfcav == 1 ) CALL zgr_isf 
    891898      ! 
    892       ! Scale factors and depth at T- and W-points 
    893 !      IF ( ln_isfcav == 0 ) THEN 
     899      IF( partial_steps ) THEN 
     900         ! Scale factors and depth at T- and W-points 
     901         ! IF ( ln_isfcav == 0 ) THEN 
    894902         DO jj = 1, ny 
    895903            DO ji = 1, nx 
     
    913921                     ELSE                                        ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 
    914922                     ENDIF 
    915    !gm Bug?  check the gdepw_1d 
     923                     !gm Bug?  check the gdepw_1d 
    916924                     !       ... on ik 
    917925                     gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )   & 
     
    931939         END DO 
    932940         ! 
    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 
    943952      ! 
    944953      ! Scale factors and depth at U-, V-, UW and VW-points 
Note: See TracChangeset for help on using the changeset viewer.