Ignore:
Timestamp:
2019-12-18T14:57:53+01:00 (10 months ago)
Author:
clem
Message:

debug: make sure bathymetry at the boundaries of the child domain is the same as the parent's one when using the same bathymetry for both parent and child grids

File:
1 edited

Legend:

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

    r12264 r12273  
    348348  ENDIF 
    349349  ! 
    350   IF( new_topo ) THEN 
    351      ! 2nd step: copy parent bathymetry at the boundaries 
    352      DO jj=1,nyfin   ! West and East 
    353         IF ( gdepw_ps_interp(nbghostcellsfine+1,jj) > 0. ) THEN 
    354            G1%gdepw_ps(1:boundary,jj) = gdepw_ps_interp(1:boundary,jj)  
    355            G1%wgt(1:boundary,jj) = 1. 
    356         ELSE 
    357            G1%gdepw_ps(1:nbghostcellsfine+1,jj)=0.  
    358         ENDIF 
    359         ! 
    360         IF ( gdepw_ps_interp(nxfin-nbghostcellsfine,jj) > 0.) THEN 
    361            G1%gdepw_ps(nxfin-boundary+1:nxfin,jj)=gdepw_ps_interp(nxfin-boundary+1:nxfin,jj) 
    362            G1%wgt(nxfin-boundary+1:nxfin,jj) = 1. 
    363         ELSE 
    364            G1%gdepw_ps(nxfin-nbghostcellsfine:nxfin,jj) = 0. 
    365         ENDIF 
     350!!$  IF( new_topo ) THEN ! clem: no, do it even when there is no new topo 
     351  ! 2nd step: copy parent bathymetry at the boundaries 
     352  DO jj=1,nyfin   ! West and East 
     353     IF ( gdepw_ps_interp(nbghostcellsfine+1,jj) > 0. ) THEN 
     354        G1%gdepw_ps(1:boundary,jj) = gdepw_ps_interp(1:boundary,jj)  
     355        G1%wgt(1:boundary,jj) = 1. 
     356     ELSE 
     357        G1%gdepw_ps(1:nbghostcellsfine+1,jj)=0.  
     358     ENDIF 
     359     ! 
     360     IF ( gdepw_ps_interp(nxfin-nbghostcellsfine,jj) > 0.) THEN 
     361        G1%gdepw_ps(nxfin-boundary+1:nxfin,jj)=gdepw_ps_interp(nxfin-boundary+1:nxfin,jj) 
     362        G1%wgt(nxfin-boundary+1:nxfin,jj) = 1. 
     363     ELSE 
     364        G1%gdepw_ps(nxfin-nbghostcellsfine:nxfin,jj) = 0. 
     365     ENDIF 
     366  END DO 
     367  ! 
     368  DO ji=1,nxfin    ! South and North 
     369     IF (gdepw_ps_interp(ji,nbghostcellsfine+1)>0.) THEN 
     370        G1%gdepw_ps(ji,1:boundary) = gdepw_ps_interp(ji,1:boundary) 
     371        G1%wgt(ji,1:boundary) = 1. 
     372     ELSE 
     373        G1%gdepw_ps(ji,1:nbghostcellsfine+1)=0.  
     374     ENDIF 
     375     ! 
     376     IF (gdepw_ps_interp(ji,nyfin-nbghostcellsfine)>0.) THEN 
     377        G1%gdepw_ps(ji,nyfin-boundary+1:nyfin)=gdepw_ps_interp(ji,nyfin-boundary+1:nyfin) 
     378        G1%wgt(ji,nyfin-boundary+1:nyfin) = 1. 
     379     ELSE 
     380        G1%gdepw_ps(ji,nyfin-nbghostcellsfine:nyfin) = 0. 
     381     ENDIF 
     382  END DO 
     383  ! 
     384  !clem: recalculate interpolation everywhere before linear connection (useless to me??) 
     385  IF( ln_agrif_domain ) THEN                 
     386     gdepw_ps_interp = 0. 
     387     CALL Check_interp(G0,gdepw_ps_interp) 
     388  ENDIF 
     389  ! 
     390  ! ------------------------------------------------------- 
     391  ! Bathymetry between boundaries and interior (npt_connect)                  
     392  ! -------------------------------------------------------- 
     393  ! Make linear connection (on npt_connect*irafx points) between the boundaries and the interior 
     394  IF( ln_agrif_domain ) THEN 
     395     boundary = (npt_copy + npt_connect)*irafx + nbghostcellsfine + 1 
     396  ELSE 
     397     boundary = (npt_copy + npt_connect)*irafx 
     398  ENDIF 
     399 
     400  IF( npt_connect > 0 ) THEN 
     401     WRITE(*,*) ' linear connection on ',npt_connect,'coarse grid points' 
     402 
     403     wghts = 1. 
     404     DO ji = boundary - npt_connect*irafx + 1 , boundary 
     405        wghts = wghts - (1. / (npt_connect*irafx + 1. ) ) 
     406        DO jj=1,nyfin 
     407           IF (G1%gdepw_ps(nbghostcellsfine+1,jj) > 0.)       G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj))   
     408        END DO 
    366409     END DO 
    367      ! 
    368      DO ji=1,nxfin    ! South and North 
    369         IF (gdepw_ps_interp(ji,nbghostcellsfine+1)>0.) THEN 
    370            G1%gdepw_ps(ji,1:boundary) = gdepw_ps_interp(ji,1:boundary) 
    371            G1%wgt(ji,1:boundary) = 1. 
    372         ELSE 
    373            G1%gdepw_ps(ji,1:nbghostcellsfine+1)=0.  
    374         ENDIF 
    375         ! 
    376         IF (gdepw_ps_interp(ji,nyfin-nbghostcellsfine)>0.) THEN 
    377            G1%gdepw_ps(ji,nyfin-boundary+1:nyfin)=gdepw_ps_interp(ji,nyfin-boundary+1:nyfin) 
    378            G1%wgt(ji,nyfin-boundary+1:nyfin) = 1. 
    379         ELSE 
    380            G1%gdepw_ps(ji,nyfin-nbghostcellsfine:nyfin) = 0. 
    381         ENDIF 
     410 
     411     wghts = 1. 
     412     DO ji = nxfin - (boundary - npt_connect*irafx), nxfin - boundary +1 , -1 
     413        wghts = wghts - (1. / (npt_connect*irafx + 1. ) ) 
     414        DO jj=1,nyfin 
     415           IF (G1%gdepw_ps(nxfin-nbghostcellsfine,jj) > 0.)   G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj)) 
     416        END DO 
    382417     END DO 
    383      ! 
    384      !clem: recalculate interpolation everywhere before linear connection (useless to me??) 
    385      IF( ln_agrif_domain ) THEN                 
    386         gdepw_ps_interp = 0. 
    387         CALL Check_interp(G0,gdepw_ps_interp) 
    388      ENDIF 
    389      ! 
    390      ! ------------------------------------------------------- 
    391      ! Bathymetry between boundaries and interior (npt_connect)                  
    392      ! -------------------------------------------------------- 
    393      ! Make linear connection (on npt_connect*irafx points) between the boundaries and the interior 
    394      IF( ln_agrif_domain ) THEN 
    395         boundary = (npt_copy + npt_connect)*irafx + nbghostcellsfine + 1 
    396      ELSE 
    397         boundary = (npt_copy + npt_connect)*irafx 
    398      ENDIF 
    399  
    400      IF( npt_connect > 0 ) THEN 
    401         WRITE(*,*) ' linear connection on ',npt_connect,'coarse grid points' 
    402  
    403         wghts = 1. 
    404         DO ji = boundary - npt_connect*irafx + 1 , boundary 
    405            wghts = wghts - (1. / (npt_connect*irafx + 1. ) ) 
    406            DO jj=1,nyfin 
    407               IF (G1%gdepw_ps(nbghostcellsfine+1,jj) > 0.)       G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj))   
    408            END DO 
     418 
     419     wghts = 1. 
     420     DO jj = boundary - npt_connect*irafy + 1 , boundary 
     421        wghts = wghts - (1. / (npt_connect*irafy + 1. ) ) 
     422        DO ji=1,nxfin 
     423           IF (G1%gdepw_ps(ji,nbghostcellsfine+1) > 0.)       G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj)) 
    409424        END DO 
    410  
    411         wghts = 1. 
    412         DO ji = nxfin - (boundary - npt_connect*irafx), nxfin - boundary +1 , -1 
    413            wghts = wghts - (1. / (npt_connect*irafx + 1. ) ) 
    414            DO jj=1,nyfin 
    415               IF (G1%gdepw_ps(nxfin-nbghostcellsfine,jj) > 0.)   G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj)) 
    416            END DO 
     425     END DO 
     426 
     427     wghts = 1. 
     428     DO jj = nyfin - (boundary - npt_connect*irafy) , nyfin - boundary +1, -1 
     429        wghts = wghts - (1. / (npt_connect*irafy + 1. ) ) 
     430        DO ji=1,nxfin 
     431           IF (G1%gdepw_ps(ji,nyfin-nbghostcellsfine) > 0.)   G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj)) 
    417432        END DO 
    418  
    419         wghts = 1. 
    420         DO jj = boundary - npt_connect*irafy + 1 , boundary 
    421            wghts = wghts - (1. / (npt_connect*irafy + 1. ) ) 
    422            DO ji=1,nxfin 
    423               IF (G1%gdepw_ps(ji,nbghostcellsfine+1) > 0.)       G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj)) 
    424            END DO 
    425         END DO 
    426  
    427         wghts = 1. 
    428         DO jj = nyfin - (boundary - npt_connect*irafy) , nyfin - boundary +1, -1 
    429            wghts = wghts - (1. / (npt_connect*irafy + 1. ) ) 
    430            DO ji=1,nxfin 
    431               IF (G1%gdepw_ps(ji,nyfin-nbghostcellsfine) > 0.)   G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj)) 
    432            END DO 
    433         END DO 
    434         IF (.NOT.identical_grids) THEN 
    435            G1%gdepw_ps(:,:) = (1.-G1%wgt(:,:)) * G1%gdepw_ps(:,:) + gdepw_ps_interp(:,:)*G1%wgt(:,:) 
    436         ENDIF 
    437  
    438      ENDIF 
    439   ENDIF 
     433     END DO 
     434     IF (.NOT.identical_grids) THEN 
     435        G1%gdepw_ps(:,:) = (1.-G1%wgt(:,:)) * G1%gdepw_ps(:,:) + gdepw_ps_interp(:,:)*G1%wgt(:,:) 
     436     ENDIF 
     437 
     438  ENDIF 
     439!!$  ENDIF 
    440440 
    441441  ! replace G1%bathy_meter by G1%gdepw_ps 
Note: See TracChangeset for help on using the changeset viewer.