Changeset 9000

Timestamp:

2017-12-13T10:09:47+01:00 (6 years ago)

Author:

timgraham

Message:

Added vertical refinement to interptrn

File:

• branches/2017/dev_METO_MERCATOR_2017_agrif/NEMOGCM/NEMO/NST_SRC/agrif_top_interp.F90 (modified) (7 diffs)

branches/2017/dev_METO_MERCATOR_2017_agrif/NEMOGCM/NEMO/NST_SRC/agrif_top_interp.F90

@@ -50 +50 @@
       !
       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
-      INTEGER :: imin, imax, jmin, jmax
+      INTEGER  ::   imin, imax, jmin, jmax, N_in, N_out
       REAL(wp) ::   zrhox , zalpha1, zalpha2, zalpha3
       REAL(wp) ::   zalpha4, zalpha5, zalpha6, zalpha7
       LOGICAL :: western_side, eastern_side,northern_side,southern_side
-
+      ! vertical interpolation:
+      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: ptab_child
+      REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin
+      REAL(wp), DIMENSION(k1:k2) :: h_in
+      REAL(wp), DIMENSION(1:jpk) :: h_out(1:jpk)
+      REAL(wp) :: h_diff, zrhoxy
+
+      zrhoxy = Agrif_rhox()*Agrif_rhoy()
       IF (before) THEN         
-         ptab(i1:i2,j1:j2,k1:k2,n1:n2) = trn(i1:i2,j1:j2,k1:k2,n1:n2)
+         DO jn = 1,jpts
+            DO jk=k1,k2
+               DO jj=j1,j2
+                 DO ji=i1,i2
+                       ptab(ji,jj,jk,jn) = trn(ji,jj,jk,jn)
+                 END DO
+               END DO
+            END DO
+         END DO
+# if defined key_vertical
+        DO jk=k1,k2
+           DO jj=j1,j2
+              DO ji=i1,i2
+                 ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) 
+              END DO
+           END DO
+        END DO
+# endif
+
       ELSE
          !
@@ -63 +88 @@
          southern_side = (nb == 2).AND.(ndir == 1)
          northern_side = (nb == 2).AND.(ndir == 2)
+
+# if defined key_vertical              
+         DO jj=j1,j2
+            DO ji=i1,i2
+               iref = ji
+               jref = jj
+               if(western_side) iref=MAX(2,ji)
+               if(eastern_side) iref=MIN(nlci-1,ji)
+               if(southern_side) jref=MAX(2,jj)
+               if(northern_side) jref=MIN(nlcj-1,jj)
+               N_in = 0
+               DO jk=k1,k2 !k2 = jpk of parent grid
+                  IF (ptab(ji,jj,jk,n2) == 0) EXIT
+                  N_in = N_in + 1
+                  tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1)
+                  h_in(N_in) = ptab(ji,jj,jk,n2)
+               END DO
+               N_out = 0
+               DO jk=1,jpk ! jpk of child grid
+                  IF (tmask(iref,jref,jk) == 0) EXIT 
+                  N_out = N_out + 1
+                  h_out(jk) = e3t_n(iref,jref,jk)
+               ENDDO
+               IF (N_in > 0) THEN
+                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
+                  DO jn=1,jptra
+                     call reconstructandremap(tabin(1:N_in,jn),h_in,ptab_child(ji,jj,1:N_out,jn),h_out,N_in,N_out)
+                  ENDDO
+               ENDIF
+            ENDDO
+         ENDDO
+# else
+         ptab_child(i1:i2,j1:j2,1:jpk,1:jptra) = ptab(i1:i2,j1:j2,1:jpk,1:jptra)
+# endif
+
          !
          zrhox = Agrif_Rhox()
@@ -89 +149 @@
          IF( eastern_side) THEN
             DO jn = 1, jptra
-               tra(nlci,j1:j2,k1:k2,jn) = zalpha1 * ptab(nlci,j1:j2,k1:k2,jn) + zalpha2 * ptab(nlci-1,j1:j2,k1:k2,jn)
+               tra(nlci,j1:j2,1:jpk,jn) = zalpha1 * ptab_child(nlci,j1:j2,1:jpk,jn) + zalpha2 * ptab_child(nlci-1,j1:j2,1:jpk,jn)
                DO jk = 1, jpkm1
                   DO jj = jmin,jmax
@@ -108 +168 @@
          IF( northern_side ) THEN            
             DO jn = 1, jptra
-               tra(i1:i2,nlcj,k1:k2,jn) = zalpha1 * ptab(i1:i2,nlcj,k1:k2,jn) + zalpha2 * ptab(i1:i2,nlcj-1,k1:k2,jn)
+               tra(i1:i2,nlcj,1:jpk,jn) = zalpha1 * ptab_child(i1:i2,nlcj,1:jpk,jn) + zalpha2 * ptab_child(i1:i2,nlcj-1,1:jpk,jn)
                DO jk = 1, jpkm1
                   DO ji = imin,imax
@@ -127 +187 @@
          IF( western_side) THEN            
             DO jn = 1, jptra
-               tra(1,j1:j2,k1:k2,jn) = zalpha1 * ptab(1,j1:j2,k1:k2,jn) + zalpha2 * ptab(2,j1:j2,k1:k2,jn)
+               tra(1,j1:j2,1:jpk,jn) = zalpha1 * ptab_child(1,j1:j2,1:jpk,jn) + zalpha2 * ptab_child(2,j1:j2,1:jpk,jn)
                DO jk = 1, jpkm1
                   DO jj = jmin,jmax
@@ -145 +205 @@
          IF( southern_side ) THEN           
             DO jn = 1, jptra
-               tra(i1:i2,1,k1:k2,jn) = zalpha1 * ptab(i1:i2,1,k1:k2,jn) + zalpha2 * ptab(i1:i2,2,k1:k2,jn)
-               DO jk=1,jpk      
+               tra(i1:i2,1,1:jpk,jn) = zalpha1 * ptab_child(i1:i2,1,1:jpk,jn) + zalpha2 * ptab_child(i1:i2,2,1:jpk,jn)
+               DO jk=1,jpkm1
                   DO ji=imin,imax
                      IF( vmask(ji,2,jk) == 0.e0 ) THEN
@@ -165 +225 @@
          ! East south
          IF ((eastern_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN
-            tra(nlci-1,2,:,:) = ptab(nlci-1,2,:,:)
+            tra(nlci-1,2,:,:) = ptab_child(nlci-1,2,:,:)
          ENDIF
          ! East north
          IF ((eastern_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN
-            tra(nlci-1,nlcj-1,:,:) = ptab(nlci-1,nlcj-1,:,:)
+            tra(nlci-1,nlcj-1,:,:) = ptab_child(nlci-1,nlcj-1,:,:)
          ENDIF
          ! West south
          IF ((western_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN
-            tra(2,2,:,:) = ptab(2,2,:,:)
+            tra(2,2,:,:) = ptab_child(2,2,:,:)
          ENDIF
          ! West north
          IF ((western_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN
-            tra(2,nlcj-1,:,:) = ptab(2,nlcj-1,:,:)
+            tra(2,nlcj-1,:,:) = ptab_child(2,nlcj-1,:,:)
          ENDIF
          !