Changeset 9628 for utils/tools

Timestamp:

2018-05-24T11:16:05+02:00 (6 years ago)

Author:

jchanut

Message:

Remove bathymetry connection if open boundary is closed, #2089

Location:

utils/tools/NESTING/src

Files:

• agrif_create_bathy.f90 (modified) (8 diffs)
• agrif_readwrite.f90 (modified) (1 diff)
• agrif_types.f90 (modified) (1 diff)

utils/tools/NESTING/src/agrif_create_bathy.f90

@@ -53 +53 @@
   REAL*8, DIMENSION(:,:),POINTER :: gdepw_ps_interp => NULL() 
   REAL*8, DIMENSION(:,:),POINTER :: save_gdepw,rx,ry,maskedtopo
-  REAL*8, DIMENSION(:,:),POINTER :: wtab  => NULL()
   REAL*8  :: Cell_lonmin,Cell_lonmax,Cell_latmin,Cell_latmax,wghts
   LOGICAL :: Pacifique=.FALSE.
@@ -460 +459 @@
         !
         boundary = connectionsize*irafx + nbghostcellsfine + 1                     
-        G1%gdepw_ps(1:boundary,:) = gdepw_ps_interp(1:boundary,:)
-        G1%gdepw_ps(:,1:boundary) = gdepw_ps_interp(:,1:boundary)
-        G1%gdepw_ps(nxfin-boundary+1:nxfin,:) = gdepw_ps_interp(nxfin-boundary+1:nxfin,:)
-        G1%gdepw_ps(:,nyfin-boundary+1:nyfin) = gdepw_ps_interp(:,nyfin-boundary+1:nyfin)
-
-
+! J chanut: exclude matching if no open boundaries
+        IF (.NOT.ASSOCIATED(G1%wgt)) &
+             ALLOCATE(G1%wgt(SIZE(G1%bathy_meter,1),SIZE(G1%bathy_meter,2)))
+        G1%wgt(:,:) = 0.
+
+        DO jj=1,nyfin
+          ! West
+          IF (gdepw_ps_interp(nbghostcellsfine+1,jj)>0.) THEN
+             G1%gdepw_ps(1:boundary,jj) = gdepw_ps_interp(1:boundary,jj) 
+             G1%wgt(1:boundary,jj) = 1.
+          ELSE
+             G1%gdepw_ps(1:nbghostcellsfine+1,jj)=0. 
+          ENDIF 
+          ! East
+          IF (gdepw_ps_interp(nxfin-nbghostcellsfine,jj)>0.) THEN
+             G1%gdepw_ps(nxfin-boundary+1:nxfin,jj)=gdepw_ps_interp(nxfin-boundary+1:nxfin,jj)
+             G1%wgt(nxfin-boundary+1:nxfin,jj) = 1.
+          ELSE
+             G1%gdepw_ps(nxfin-nbghostcellsfine:nxfin,jj) = 0.
+          ENDIF
+        END DO
+        DO ji=1,nxfin
+          ! South 
+          IF (gdepw_ps_interp(ji,nbghostcellsfine+1)>0.) THEN
+             G1%gdepw_ps(ji,1:boundary) = gdepw_ps_interp(ji,1:boundary)
+             G1%wgt(ji,1:boundary) = 1.
+          ELSE
+             G1%gdepw_ps(ji,1:nbghostcellsfine+1)=0. 
+          ENDIF
+          ! North
+          IF (gdepw_ps_interp(ji,nyfin-nbghostcellsfine)>0.) THEN
+             G1%gdepw_ps(ji,nyfin-boundary+1:nyfin)=gdepw_ps_interp(ji,nyfin-boundary+1:nyfin)
+             G1%wgt(ji,nyfin-boundary+1:nyfin) = 1.
+          ELSE
+             G1%gdepw_ps(ji,nyfin-nbghostcellsfine:nyfin) = 0.
+          ENDIF
+        END DO
+        
+!        G1%gdepw_ps(1:boundary,:) = gdepw_ps_interp(1:boundary,:)
+!        G1%gdepw_ps(:,1:boundary) = gdepw_ps_interp(:,1:boundary)
+!        G1%gdepw_ps(nxfin-boundary+1:nxfin,:) = gdepw_ps_interp(nxfin-boundary+1:nxfin,:)
+!        G1%gdepw_ps(:,nyfin-boundary+1:nyfin) = gdepw_ps_interp(:,nyfin-boundary+1:nyfin)
         !                   
 
@@ -485 +520 @@
         !
         !
-        IF (.NOT.ASSOCIATED(wtab)) &
-             ALLOCATE(wtab(SIZE(G1%bathy_meter,1),SIZE(G1%bathy_meter,2)))
-        wtab(:,:) = 0.
         wghts = 1.
         DO ji = boundary + 1 , boundary + nb_connection_pts * irafx
             wghts = wghts - (1. / (nb_connection_pts*irafx - 1. ) )
-            DO jj=boundary+1,nyfin-boundary
-               wtab(ji,jj) = MAX(wghts, wtab(ji,jj))  
+            DO jj=1,nyfin
+               IF (G1%gdepw_ps(nbghostcellsfine+1,jj) > 0.) THEN
+                  G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj))  
+               ENDIF
             END DO
         END DO 
@@ -499 +533 @@
         DO ji = nxfin - boundary , nxfin - boundary - nb_connection_pts * irafx +1 , -1
             wghts = wghts - (1. / (nb_connection_pts*irafx - 1. ) )
-            DO jj=boundary+1,nyfin-boundary
-               wtab(ji,jj) = MAX(wghts, wtab(ji,jj))
+            DO jj=1,nyfin
+               IF (G1%gdepw_ps(nxfin-nbghostcellsfine,jj) > 0.) THEN
+                  G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj))
+               ENDIF
             END DO
         END DO  
@@ -507 +543 @@
         DO jj = boundary + 1 , boundary + nb_connection_pts * irafy
             wghts = wghts - (1. / (nb_connection_pts*irafy - 1. ) )
-            DO ji=boundary+1,nxfin-boundary
-               wtab(ji,jj) = MAX(wghts, wtab(ji,jj))
+            DO ji=1,nxfin
+               IF (G1%gdepw_ps(ji,nbghostcellsfine+1) > 0.) THEN
+                  G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj))
+               ENDIF
             END DO
         END DO
@@ -515 +553 @@
         DO jj = nyfin - boundary , nyfin - boundary - nb_connection_pts * irafy +1, -1
             wghts = wghts - (1. / (nb_connection_pts*irafy - 1. ) )
-            DO ji=boundary+1,nxfin-boundary
-               wtab(ji,jj) = MAX(wghts, wtab(ji,jj))
+            DO ji=1,nxfin
+               IF (G1%gdepw_ps(ji,nyfin-nbghostcellsfine) > 0.) THEN
+                  G1%wgt(ji,jj) = MAX(wghts, G1%wgt(ji,jj))
+               ENDIF
             END DO
         END DO
 
-        G1%gdepw_ps(:,:) = (1.-wtab(:,:)) * G1%gdepw_ps(:,:) + gdepw_ps_interp(:,:)*wtab(:,:)
+        G1%gdepw_ps(:,:) = (1.-G1%wgt(:,:)) * G1%gdepw_ps(:,:) + gdepw_ps_interp(:,:)*G1%wgt(:,:)
 
         G1%bathy_meter = G1%gdepw_ps
@@ -526 +566 @@
         connectionsize = 2 
         ! 
-        IF(smoothing) THEN 
+! Chanut: remove smoothing if child grid bathymetry is found to already exist
+        IF(smoothing.AND.(.NOT.identical_grids)) THEN 
 
            !
            ! Smoothing to connect the connection zone (3 + nb_connection_pts coarse grid cells) and the interior domain
            !
-           boundary = (connectionsize+nb_connection_pts)*irafx + nbghostcellsfine + 1 
-           CALL smooth_topo(G1%gdepw_ps(boundary:nxfin-boundary+1,boundary:nyfin-boundary+1),nbiter)
-           G1%bathy_meter = G1%gdepw_ps                         
+! Chanut: smoothing everywhere then discard result in connection zone
+           CALL smooth_topo(G1%gdepw_ps(1:nxfin,1:nyfin),nbiter)
+           WHERE (G1%wgt(:,:)==0.) G1%bathy_meter(:,:) = G1%gdepw_ps(:,:)
+!           boundary = (connectionsize+nb_connection_pts)*irafx + nbghostcellsfine + 1 
+!           CALL smooth_topo(G1%gdepw_ps(boundary:nxfin-boundary+1,boundary:nyfin-boundary+1),nbiter)
+!           G1%bathy_meter = G1%gdepw_ps                         
         ENDIF
 
@@ -578 +622 @@
         ENDIF
         !
+        ! Chanut: Compute partial step bathymetry once more
+        CALL get_partial_steps(G1)                 ! compute gdepw_ps for G1
+
         IF(bathy_update) CALL Update_Parent_Bathy( G0,G1 )                  
         !

utils/tools/NESTING/src/agrif_readwrite.f90

@@ -460 +460 @@
     CALL Write_Ncdf_var('nav_lat'         ,dimnames,name,Grid%nav_lat    ,'float')
     CALL Write_Ncdf_var(parent_batmet_name,dimnames,name,Grid%bathy_meter,'float')
+    CALL Write_Ncdf_var('weight'          ,dimnames,name,Grid%wgt        ,'float')
     !
     CALL Copy_Ncdf_att('nav_lon'         ,TRIM(parent_bathy_meter),name,MINVAL(Grid%nav_lon),MAXVAL(Grid%nav_lon))

utils/tools/NESTING/src/agrif_types.f90

@@ -22 +22 @@
      REAL*8, DIMENSION(:,:), POINTER :: bathy_level => NULL()
      REAL*8, DIMENSION(:,:), POINTER :: bathy_meter => NULL()
+     REAL*8, DIMENSION(:,:), POINTER :: wgt => NULL()
      REAL*8, DIMENSION(:,:,:),POINTER :: fmask,umask,vmask,tmask => NULL()
      REAL*8, DIMENSION(:,:,:),POINTER :: e3t_ps,e3w_ps,gdept_ps,gdepwps => NULL()