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.
Changeset 12253 for utils/tools/NESTING/src/agrif_create_bathy.f90 – NEMO

Ignore:
Timestamp:
2019-12-16T10:59:22+01:00 (4 years 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...

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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 
Note: See TracChangeset for help on using the changeset viewer.