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 8135 – NEMO

Changeset 8135


Ignore:
Timestamp:
2017-06-05T11:01:20+02:00 (7 years ago)
Author:
timgraham
Message:

Added changes to code in NST_SRC from dev_r5803_UKMO_AGRIF_Vert_interp

Location:
branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_oce.F90

    r5656 r8135  
    1717 
    1818   PUBLIC agrif_oce_alloc ! routine called by nemo_init in nemogcm.F90 
     19   PUBLIC reconstructandremap 
    1920 
    2021   !                                              !!* Namelist namagrif: AGRIF parameters 
     
    104105   END FUNCTION agrif_oce_alloc 
    105106 
     107      subroutine reconstructandremap(tabin,hin,tabout,hout,N,Nout)       
     108      implicit none 
     109      integer N, Nout 
     110      real tabin(N), tabout(Nout) 
     111      real hin(N), hout(Nout) 
     112      real coeffremap(N,3),zwork(N,3) 
     113      real zwork2(N+1,3) 
     114      integer k 
     115      double precision, parameter :: dsmll=1.0d-8   
     116      real q,q01,q02,q001,q002,q0 
     117      real z_win(1:N+1), z_wout(1:Nout+1) 
     118      real,parameter :: dpthin = 1.D-3 
     119      integer :: k1, kbox, ktop, ka, kbot 
     120      real :: tsum, qbot, rpsum, zbox, ztop, zthk, zbot, offset, qtop 
     121 
     122      z_win(1)=0.; z_wout(1)= 0. 
     123      do k=1,N 
     124       z_win(k+1)=z_win(k)+hin(k) 
     125      enddo  
     126       
     127      do k=1,Nout 
     128       z_wout(k+1)=z_wout(k)+hout(k)        
     129      enddo        
     130 
     131        do k=2,N 
     132          zwork(k,1)=1./(hin(k-1)+hin(k)) 
     133        enddo 
     134         
     135        do k=2,N-1 
     136          q0 = 1./(hin(k-1)+hin(k)+hin(k+1)) 
     137          zwork(k,2)=hin(k-1)+2.*hin(k)+hin(k+1) 
     138          zwork(k,3)=q0 
     139        enddo        
     140       
     141        do k= 2,N 
     142        zwork2(k,1)=zwork(k,1)*(tabin(k)-tabin(k-1)) 
     143        enddo 
     144         
     145        coeffremap(:,1) = tabin(:) 
     146  
     147         do k=2,N-1 
     148        q001 = hin(k)*zwork2(k+1,1) 
     149        q002 = hin(k)*zwork2(k,1)         
     150        if (q001*q002 < 0) then 
     151          q001 = 0. 
     152          q002 = 0. 
     153        endif 
     154        q=zwork(k,2) 
     155        q01=q*zwork2(k+1,1) 
     156        q02=q*zwork2(k,1) 
     157        if (abs(q001) > abs(q02)) q001 = q02 
     158        if (abs(q002) > abs(q01)) q002 = q01 
     159         
     160        q=(q001-q002)*zwork(k,3) 
     161        q001=q001-q*hin(k+1) 
     162        q002=q002+q*hin(k-1) 
     163         
     164        coeffremap(k,3)=coeffremap(k,1)+q001 
     165        coeffremap(k,2)=coeffremap(k,1)-q002 
     166         
     167        zwork2(k,1)=(2.*q001-q002)**2 
     168        zwork2(k,2)=(2.*q002-q001)**2 
     169        enddo 
     170         
     171        do k=1,N 
     172        if     (k.eq.1 .or. k.eq.N .or.   hin(k).le.dpthin) then 
     173        coeffremap(k,3) = coeffremap(k,1) 
     174        coeffremap(k,2) = coeffremap(k,1) 
     175        zwork2(k,1) = 0. 
     176        zwork2(k,2) = 0. 
     177        endif 
     178        enddo 
     179         
     180        do k=2,N 
     181        q002=max(zwork2(k-1,2),dsmll) 
     182        q001=max(zwork2(k,1),dsmll) 
     183        zwork2(k,3)=(q001*coeffremap(k-1,3)+q002*coeffremap(k,2))/(q001+q002) 
     184        enddo 
     185         
     186        zwork2(1,3) = 2*coeffremap(1,1)-zwork2(2,3) 
     187        zwork2(N+1,3)=2*coeffremap(N,1)-zwork2(N,3) 
     188  
     189        do k=1,N 
     190        q01=zwork2(k+1,3)-coeffremap(k,1) 
     191        q02=coeffremap(k,1)-zwork2(k,3) 
     192        q001=2.*q01 
     193        q002=2.*q02 
     194        if (q01*q02<0) then 
     195          q01=0. 
     196          q02=0. 
     197        elseif (abs(q01)>abs(q002)) then 
     198          q01=q002 
     199        elseif (abs(q02)>abs(q001)) then 
     200          q02=q001 
     201        endif 
     202        coeffremap(k,2)=coeffremap(k,1)-q02 
     203        coeffremap(k,3)=coeffremap(k,1)+q01 
     204        enddo 
     205 
     206      zbot=0.0 
     207      kbot=1 
     208      do k=1,Nout 
     209        ztop=zbot  !top is bottom of previous layer 
     210        ktop=kbot 
     211        if     (ztop.ge.z_win(ktop+1)) then 
     212          ktop=ktop+1 
     213        endif 
     214         
     215        zbot=z_wout(k+1) 
     216        zthk=zbot-ztop 
     217 
     218        if     (zthk.gt.dpthin .and. ztop.lt.z_wout(Nout+1)) then 
     219 
     220          kbot=ktop 
     221          do while (z_win(kbot+1).lt.zbot.and.kbot.lt.N) 
     222            kbot=kbot+1 
     223          enddo 
     224          zbox=zbot 
     225          do k1= k+1,Nout 
     226            if     (z_wout(k1+1)-z_wout(k1).gt.dpthin) then 
     227              exit !thick layer 
     228            else 
     229              zbox=z_wout(k1+1)  !include thin adjacent layers 
     230              if     (zbox.eq.z_wout(Nout+1)) then 
     231                exit !at bottom 
     232              endif 
     233            endif 
     234          enddo 
     235          zthk=zbox-ztop 
     236 
     237          kbox=ktop 
     238          do while (z_win(kbox+1).lt.zbox.and.kbox.lt.N) 
     239            kbox=kbox+1 
     240          enddo 
     241           
     242          if     (ktop.eq.kbox) then 
     243 
     244 
     245            if     (z_wout(k)  .ne.z_win(kbox)   .or.z_wout(k+1).ne.z_win(kbox+1)     ) then 
     246 
     247              if     (hin(kbox).gt.dpthin) then 
     248                q001 = (zbox-z_win(kbox))/hin(kbox) 
     249                q002 = (ztop-z_win(kbox))/hin(kbox) 
     250                q01=q001**2+q002**2+q001*q002+1.-2.*(q001+q002) 
     251                q02=q01-1.+(q001+q002) 
     252                q0=1.-q01-q02 
     253              else 
     254                q0 = 1.0 
     255                q01 = 0. 
     256                q02 = 0. 
     257              endif 
     258          tabout(k)=q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3) 
     259               
     260            else 
     261            tabout(k) = tabin(kbox) 
     262               
     263            endif  
     264 
     265          else 
     266 
     267            if     (ktop.le.k .and. kbox.ge.k) then 
     268              ka = k 
     269            elseif (kbox-ktop.ge.3) then 
     270              ka = (kbox+ktop)/2 
     271            elseif (hin(ktop).ge.hin(kbox)) then 
     272              ka = ktop 
     273            else 
     274              ka = kbox 
     275            endif !choose ka 
     276 
     277            offset=coeffremap(ka,1) 
     278 
     279            qtop = z_win(ktop+1)-ztop !partial layer thickness 
     280            if     (hin(ktop).gt.dpthin) then 
     281              q=(ztop-z_win(ktop))/hin(ktop) 
     282              q01=q*(q-1.) 
     283              q02=q01+q 
     284              q0=1-q01-q02             
     285            else 
     286              q0 = 1. 
     287              q01 = 0. 
     288              q02 = 0. 
     289            endif 
     290             
     291            tsum =((q0*coeffremap(ktop,1)+q01*coeffremap(ktop,2)+q02*coeffremap(ktop,3))-offset)*qtop 
     292 
     293            do k1= ktop+1,kbox-1 
     294              tsum =tsum +(coeffremap(k1,1)-offset)*hin(k1) 
     295            enddo !k1 
     296 
     297             
     298            qbot = zbox-z_win(kbox) !partial layer thickness 
     299            if     (hin(kbox).gt.dpthin) then 
     300              q=qbot/hin(kbox) 
     301              q01=(q-1.)**2 
     302              q02=q01-1.+q 
     303              q0=1-q01-q02                             
     304            else 
     305              q0 = 1.0 
     306              q01 = 0. 
     307              q02 = 0. 
     308            endif 
     309            
     310            tsum = tsum +((q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3))-offset)*qbot 
     311             
     312            rpsum=1.0d0/zthk 
     313              tabout(k)=offset+tsum*rpsum 
     314               
     315          endif !single or multiple layers 
     316        else 
     317        if (k==1) then 
     318        write(*,'(a7,i4,i4,3f12.5)')'problem = ',N,Nout,zthk,z_wout(k+1),hout(1) 
     319        endif 
     320         tabout(k) = tabout(k-1) 
     321           
     322        endif !normal:thin layer 
     323      enddo !k 
     324             
     325      return 
     326      end subroutine reconstructandremap 
     327       
    106328#endif 
    107329   !!====================================================================== 
  • branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90

    r7646 r8135  
    576576      !!----------------------------------------------------------------------   
    577577      ! 
     578      return 
     579       
    578580      zalpha = REAL( Agrif_NbStepint() + Agrif_IRhot() - 1, wp ) / REAL( Agrif_IRhot(), wp ) 
    579581      IF( zalpha > 1. )   zalpha = 1. 
     
    603605      REAL(wp) ::   zrhox , zalpha1, zalpha2, zalpha3 
    604606      REAL(wp) ::   zalpha4, zalpha5, zalpha6, zalpha7 
    605       LOGICAL  ::   western_side, eastern_side,northern_side,southern_side 
    606       !!---------------------------------------------------------------------- 
    607       ! 
     607      LOGICAL :: western_side, eastern_side,northern_side,southern_side 
     608! VERTICAL REFINEMENT BEGIN 
     609      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: ptab_child 
     610      REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin 
     611      REAL(wp) :: h_in(k1:k2) 
     612      REAL(wp) :: h_out(1:jpk) 
     613      INTEGER :: N_in, N_out 
     614      REAL(wp) :: h_diff 
     615      REAL(wp) :: zrhoxy 
     616! VERTICAL REFINEMENT END 
     617 
     618      zrhoxy = Agrif_rhox()*Agrif_rhoy() 
    608619      IF (before) THEN          
    609          ptab(i1:i2,j1:j2,k1:k2,n1:n2) = tsn(i1:i2,j1:j2,k1:k2,n1:n2) 
    610       ELSE 
     620         DO jn = n1,n2-1 
     621            DO jk=k1,k2 
     622               DO jj=j1,j2 
     623                  DO ji=i1,i2 
     624                     ptab(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) 
     625                  END DO 
     626               END DO 
     627            END DO 
     628         END DO 
     629         DO jk=k1,k2 
     630            DO jj=j1,j2 
     631               DO ji=i1,i2 
     632                  ptab(ji,jj,jk,n2) = tmask(ji,jj,jk) * e1e2t(ji,jj) * e3t_n(ji,jj,jk)  
     633               END DO 
     634            END DO 
     635         END DO 
     636 
     637      ELSE  
     638! VERTICAL REFINEMENT BEGIN 
     639 
     640         ptab_child(:,:,:,:) = 0. 
     641         do jj=j1,j2 
     642         do ji=i1,i2 
     643           N_in = 0 
     644           DO jk=k1,k2 !k2 = jpk of parent grid 
     645             IF (ptab(ji,jj,jk,n2) == 0) EXIT 
     646             N_in = N_in + 1 
     647             tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1)/ptab(ji,jj,jk,n2) 
     648             h_in(N_in) = ptab(ji,jj,jk,n2)/(e1e2t(ji,jj)*zrhoxy) 
     649           END DO 
     650           N_out = 0 
     651           DO jk=1,jpk ! jpk of child grid 
     652             IF (tmask(ji,jj,jk) == 0) EXIT ! TODO: Will not work with ISF. !This doesn't seem to work at the moment in GYRE but is OK in overflow model 
     653             N_out = N_out + 1 
     654             h_out(jk) = e3t_n(ji,jj,jk) !Child grid scale factors. Could multiply by e1e2t here instead of division above 
     655           ENDDO 
     656           IF (N_in > 0) THEN 
     657             h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
     658!             if (h_diff > 0) then 
     659!               h_in(N_in+1) = h_diff 
     660!               N_in = N_in + 1 
     661!             else 
     662!               h_out(N_out+1) = -h_diff 
     663!               N_out = N_out + 1 
     664!             endif  
     665           ptab(ji,jj,k2,:) = ptab(ji,jj,k2-1,:) !what is this line for????? 
     666           do jn=1,jpts 
     667             call reconstructandremap(tabin(1:N_in,jn),h_in,ptab_child(ji,jj,1:N_out,jn),h_out,N_in,N_out) 
     668           enddo 
     669!           if (abs(h_diff) > 1000.) then 
     670!           do jn=1,jpts 
     671!             do jk=1,N_out 
     672!             print *,'AVANT APRES = ',ji,jj,jk,N_out,ptab(ji,jj,jk,jn),ptab_child(ji,jj,jk,jn) 
     673!             enddo 
     674!           enddo 
     675!         endif 
     676           ENDIF 
     677         enddo 
     678         enddo 
     679    
     680 
     681! VERTICAL REFINEMENT END 
     682 
    611683         ! 
    612684         western_side  = (nb == 1).AND.(ndir == 1) 
     
    640712         IF( eastern_side ) THEN 
    641713            DO jn = 1, jpts 
    642                tsa(nlci,j1:j2,k1:k2,jn) = zalpha1 * ptab(nlci,j1:j2,k1:k2,jn) + zalpha2 * ptab(nlci-1,j1:j2,k1:k2,jn) 
     714               tsa(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) 
    643715               DO jk = 1, jpkm1 
    644716                  DO jj = jmin,jmax 
     
    660732         IF( northern_side ) THEN             
    661733            DO jn = 1, jpts 
    662                tsa(i1:i2,nlcj,k1:k2,jn) = zalpha1 * ptab(i1:i2,nlcj,k1:k2,jn) + zalpha2 * ptab(i1:i2,nlcj-1,k1:k2,jn) 
     734               tsa(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) 
    663735               DO jk = 1, jpkm1 
    664736                  DO ji = imin,imax 
     
    680752         IF( western_side ) THEN             
    681753            DO jn = 1, jpts 
    682                tsa(1,j1:j2,k1:k2,jn) = zalpha1 * ptab(1,j1:j2,k1:k2,jn) + zalpha2 * ptab(2,j1:j2,k1:k2,jn) 
     754               tsa(1,j1:j2,1:jpk,jn) = zalpha1 * ptab_child(1,j1:j2,1:jpk,jn) + zalpha2 * ptab_child(2,j1:j2,1:jpk,jn) 
    683755               DO jk = 1, jpkm1 
    684756                  DO jj = jmin,jmax 
     
    699771         IF( southern_side ) THEN            
    700772            DO jn = 1, jpts 
    701                tsa(i1:i2,1,k1:k2,jn) = zalpha1 * ptab(i1:i2,1,k1:k2,jn) + zalpha2 * ptab(i1:i2,2,k1:k2,jn) 
     773               tsa(i1:i2,1,1:jpk,jn) = zalpha1 * ptab_child(i1:i2,1,1:jpk,jn) + zalpha2 * ptab_child(i1:i2,2,1:jpk,jn) 
    702774               DO jk = 1, jpk       
    703775                  DO ji=imin,imax 
     
    720792         ! East south 
    721793         IF ((eastern_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN 
    722             tsa(nlci-1,2,:,:) = ptab(nlci-1,2,:,:) 
     794            tsa(nlci-1,2,:,:) = ptab_child(nlci-1,2,:,1:jpts) 
    723795         ENDIF 
    724796         ! East north 
    725797         IF ((eastern_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN 
    726             tsa(nlci-1,nlcj-1,:,:) = ptab(nlci-1,nlcj-1,:,:) 
     798            tsa(nlci-1,nlcj-1,:,:) = ptab_child(nlci-1,nlcj-1,:,1:jpts) 
    727799         ENDIF 
    728800         ! West south 
    729801         IF ((western_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN 
    730             tsa(2,2,:,:) = ptab(2,2,:,:) 
     802            tsa(2,2,:,:) = ptab_child(2,2,:,1:jpts) 
    731803         ENDIF 
    732804         ! West north 
    733805         IF ((western_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN 
    734             tsa(2,nlcj-1,:,:) = ptab(2,nlcj-1,:,:) 
     806            tsa(2,nlcj-1,:,:) = ptab_child(2,nlcj-1,:,1:jpts) 
    735807         ENDIF 
    736808         ! 
     
    771843      !!---------------------------------------------------------------------- 
    772844      !!   *** ROUTINE interpun *** 
    773       !!---------------------------------------------------------------------- 
    774       INTEGER                               , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2 
    775       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab 
    776       LOGICAL                               , INTENT(in   ) ::   before 
    777       ! 
    778       INTEGER  ::   ji, jj, jk 
    779       REAL(wp) ::   zrhoy   
    780       !!---------------------------------------------------------------------- 
    781       ! 
    782       IF( before ) THEN  
    783          DO jk = k1, jpk 
    784             ptab(i1:i2,j1:j2,jk) = e2u(i1:i2,j1:j2) * e3u_n(i1:i2,j1:j2,jk) * un(i1:i2,j1:j2,jk) 
     845      !!---------------------------------------------     
     846      !! 
     847      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2 
     848      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab 
     849      LOGICAL, INTENT(in) :: before 
     850      INTEGER, INTENT(in) :: nb , ndir 
     851      !! 
     852      INTEGER :: ji,jj,jk 
     853      REAL(wp) :: zrhoy 
     854! VERTICAL REFINEMENT BEGIN 
     855      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ptab_child 
     856      REAL(wp), DIMENSION(k1:k2) :: tabin 
     857      REAL(wp) :: h_in(k1:k2) 
     858      REAL(wp) :: h_out(1:jpk) 
     859      INTEGER :: N_in, N_out 
     860      REAL(wp) :: h_diff 
     861      LOGICAL :: western_side, eastern_side 
     862      INTEGER :: iref 
     863 
     864! VERTICAL REFINEMENT END 
     865      !!---------------------------------------------     
     866      ! 
     867      zrhoy = Agrif_rhoy() 
     868      IF (before) THEN  
     869         !We can't use zero as the special value because we need to include zeros 
     870         !when interpolating the scale factors 
     871         IF(Agrif_UseSpecialValue) THEN  
     872             Agrif_SpecialValue = -999._wp 
     873         ELSE 
     874             Agrif_SpecialValue = 0._wp 
     875         ENDIF 
     876         DO jk=1,jpk 
     877            DO jj=j1,j2 
     878               DO ji=i1,i2 
     879                  ptab(ji,jj,jk,1) = (e2u(ji,jj) * e3u_n(ji,jj,jk) * un(ji,jj,jk)*umask(ji,jj,jk)) - & 
     880                                   & ((umask(ji,jj,jk)-1) * Agrif_SpecialValue) 
     881                  ptab(ji,jj,jk,2) = (umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk)) 
     882               END DO 
     883            END DO 
    785884         END DO 
    786885      ELSE 
    787          zrhoy = Agrif_Rhoy() 
     886! VERTICAL REFINEMENT BEGIN 
     887         western_side  = (nb == 1).AND.(ndir == 1) 
     888         eastern_side  = (nb == 1).AND.(ndir == 2) 
     889 
     890         Agrif_SpecialValue = 0._wp ! reset specialvalue to zero now interpolation completed 
     891 
     892         ptab_child(:,:,:) = 0. 
     893         DO jj=j1,j2 
     894            DO ji=i1,i2 
     895               iref = ji 
     896               IF (western_side) iref = 2 
     897               IF (eastern_side) iref = nlci-2 
     898 
     899               N_in = 0 
     900               DO jk=k1,k2 
     901                  IF (ptab(ji,jj,jk,2) == 0) EXIT 
     902                  N_in = N_in + 1 
     903                  tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2) 
     904                  h_in(N_in) = ptab(ji,jj,jk,2)/(e2u(ji,jj) * zrhoy)  
     905              ENDDO 
     906          
     907              IF (N_in == 0) THEN 
     908                 ptab_child(ji,jj,:) = 0. 
     909                 CYCLE 
     910              ENDIF 
     911          
     912              N_out = 0 
     913              DO jk=1,jpk 
     914                 if (umask(ji,jj,jk) == 0) EXIT 
     915                 N_out = N_out + 1 
     916                 h_out(N_out) = e3u_n(ji,jj,jk) 
     917              ENDDO 
     918          
     919              IF (N_out == 0) THEN 
     920                 ptab_child(ji,jj,:) = 0. 
     921                 CYCLE 
     922              ENDIF 
     923          
     924!              IF (N_in * N_out > 0) THEN 
     925!                 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
     926! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly 
     927!                 if (h_diff < 0.) then 
     928!                    print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in)) 
     929!                    stop 
     930!                 endif 
     931!              ENDIF 
     932              call reconstructandremap(tabin(1:N_in),h_in(1:N_in),ptab_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out) 
     933          
     934            ENDDO 
     935         ENDDO 
     936 
     937! in the following 
     938! remove division of ua by fs e3u (already done) and also zrhoy and e2u 
     939! VERTICAL REFINEMENT END 
     940 
    788941         DO jk = 1, jpkm1 
    789942            DO jj=j1,j2 
    790                ua(i1:i2,jj,jk) = ptab(i1:i2,jj,jk) / ( zrhoy * e2u(i1:i2,jj) * e3u_n(i1:i2,jj,jk) ) 
     943               ua(i1:i2,jj,jk) = ptab_child(i1:i2,jj,jk) 
     944!/(zrhoy*e2u(i1:i2,jj))) 
    791945            END DO 
    792946         END DO 
     
    800954      !!   *** ROUTINE interpvn *** 
    801955      !!---------------------------------------------------------------------- 
    802       INTEGER                               , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2 
    803       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab 
    804       LOGICAL                               , INTENT(in   ) ::   before 
    805       ! 
    806       INTEGER  ::   ji, jj, jk 
    807       REAL(wp) ::   zrhox   
    808       !!---------------------------------------------------------------------- 
     956      ! 
     957      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2 
     958      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab 
     959      LOGICAL, INTENT(in) :: before 
     960      INTEGER, INTENT(in) :: nb , ndir 
     961      ! 
     962      INTEGER :: ji,jj,jk 
     963      REAL(wp) :: zrhox 
     964! VERTICAL REFINEMENT BEGIN 
     965      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ptab_child 
     966      REAL(wp), DIMENSION(k1:k2) :: tabin 
     967      REAL(wp) :: h_in(k1:k2) 
     968      REAL(wp) :: h_out(1:jpk) 
     969      INTEGER :: N_in, N_out 
     970      REAL(wp) :: h_diff 
     971      LOGICAL :: northern_side,southern_side 
     972      INTEGER :: jref 
     973 
     974! VERTICAL REFINEMENT END 
     975      !!---------------------------------------------     
    809976      !       
    810       IF( before ) THEN       !interpv entre 1 et k2 et interpv2d en jpkp1 
    811          DO jk = k1, jpk 
    812             ptab(i1:i2,j1:j2,jk) = e1v(i1:i2,j1:j2) * e3v_n(i1:i2,j1:j2,jk) * vn(i1:i2,j1:j2,jk) 
    813          END DO 
    814       ELSE           
    815          zrhox= Agrif_Rhox() 
    816          DO jk = 1, jpkm1 
    817             va(i1:i2,j1:j2,jk) = ptab(i1:i2,j1:j2,jk) / ( zrhox * e1v(i1:i2,j1:j2) * e3v_n(i1:i2,j1:j2,jk) ) 
     977      zrhox = Agrif_rhox() 
     978      IF (before) THEN           
     979         IF(Agrif_UseSpecialValue) THEN  
     980             Agrif_SpecialValue = -999._wp 
     981         ELSE 
     982             Agrif_SpecialValue = 0._wp 
     983         ENDIF 
     984         DO jk=k1,k2 
     985            DO jj=j1,j2 
     986               DO ji=i1,i2 
     987                  ptab(ji,jj,jk,1) = e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk) 
     988                  ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v_n(ji,jj,jk) 
     989               END DO 
     990            END DO 
     991         END DO 
     992      ELSE         
     993! VERTICAL REFINEMENT BEGIN 
     994         ptab_child(:,:,:) = 0. 
     995         southern_side = (nb == 2).AND.(ndir == 1) 
     996         northern_side = (nb == 2).AND.(ndir == 2) 
     997 
     998         Agrif_SpecialValue = 0._wp !Reset special value to zero now interpolation is done 
     999 
     1000         do jj=j1,j2 
     1001            jref = jj 
     1002            IF (southern_side) jref = 2 
     1003            IF (northern_side) jref = nlcj-2 
     1004            do ji=i1,i2 
     1005 
     1006               N_in = 0 
     1007               do jk=k1,k2 
     1008                  if (ptab(ji,jj,jk,2) == 0) EXIT 
     1009                  N_in = N_in + 1 
     1010                  tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2) 
     1011                  h_in(N_in) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox) 
     1012               enddo 
     1013               IF (N_in == 0) THEN 
     1014                  ptab_child(ji,jj,:) = 0. 
     1015                  CYCLE 
     1016               ENDIF 
     1017          
     1018               N_out = 0 
     1019               do jk=1,jpk 
     1020                  if (vmask(ji,jref,jk) == 0) EXIT 
     1021                  N_out = N_out + 1 
     1022                  h_out(N_out) = e3v_n(ji,jj,jk) 
     1023               enddo 
     1024               IF (N_out == 0) THEN 
     1025                 ptab_child(ji,jj,:) = 0. 
     1026                 CYCLE 
     1027               ENDIF 
     1028 
     1029!               IF (N_in * N_out > 0) THEN 
     1030!                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
     1031! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly 
     1032!                  if (h_diff < 0.) then 
     1033!                     print *,'CHECK YOUR BATHY interpvn...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in)) 
     1034!                     stop 
     1035!                  endif 
     1036!               ENDIF 
     1037               call reconstructandremap(tabin(1:N_in),h_in(1:N_in),ptab_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out) 
     1038          
     1039            enddo 
     1040         enddo 
     1041! in the following 
     1042! remove division of va by fs e3v, zrhox and e1v (already done) 
     1043! VERTICAL REFINEMENT END 
     1044         DO jk=1,jpkm1 
     1045            DO jj=j1,j2 
     1046               va(i1:i2,jj,jk) = ptab_child(i1:i2,jj,jk) 
     1047            END DO 
    8181048         END DO 
    8191049      ENDIF 
  • branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90

    r7646 r8135  
    158158      INTEGER, INTENT(in) :: kt 
    159159      !        
     160      return 
     161       
    160162      IF( ( Agrif_NbStepint() .NE. 0 ) .AND. (kt /= 0) ) RETURN 
    161163#  if defined TWO_WAY 
     
    185187      !! 
    186188      INTEGER :: ji,jj,jk,jn 
    187       !!--------------------------------------------- 
    188       ! 
    189       IF (before) THEN 
    190          DO jn = n1,n2 
     189! VERTICAL REFINEMENT BEGIN 
     190      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: tabres_child 
     191      REAL(wp) :: h_in(k1:k2) 
     192      REAL(wp) :: h_out(1:jpk) 
     193      INTEGER  :: N_in, N_out 
     194      REAL(wp) :: h_diff 
     195      REAL(wp) :: zrho_xy 
     196      REAL(wp) :: tabin(k1:k2,n1:n2) 
     197! VERTICAL REFINEMENT END 
     198      !!--------------------------------------------- 
     199      ! 
     200      IF (before) THEN 
     201         zrho_xy = Agrif_rhox() * Agrif_rhoy()  
     202         DO jn = n1,n2-1 
    191203            DO jk=k1,k2 
    192204               DO jj=j1,j2 
    193205                  DO ji=i1,i2 
    194                      tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) 
     206                     tabres(ji,jj,jk,jn) = zrho_xy * tsn(ji,jj,jk,jn) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) 
    195207                  END DO 
    196208               END DO 
    197209            END DO 
    198210         END DO 
    199       ELSE 
     211         DO jk=k1,k2 
     212            DO jj=j1,j2 
     213               DO ji=i1,i2 
     214                  tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * zrho_xy * e1e2t(ji,jj) * e3t_n(ji,jj,jk)  
     215               END DO 
     216            END DO 
     217         END DO 
     218          
     219      ELSE 
     220! VERTICAL REFINEMENT BEGIN 
     221         tabres_child(:,:,:,:) = 0. 
     222 
     223         DO jj=j1,j2 
     224         DO ji=i1,i2 
     225           N_in = 0 
     226           DO jk=k1,k2 !k2 = jpk of child grid 
     227             IF (tabres(ji,jj,jk,n2) == 0) EXIT 
     228             N_in = N_in + 1 
     229             tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2) 
     230             h_in(N_in) = tabres(ji,jj,jk,n2)/e1e2t(ji,jj) 
     231           ENDDO 
     232           N_out = 0 
     233           DO jk=1,jpk ! jpk of parent grid 
     234             IF (tmask(ji,jj,jk) == 0) EXIT ! TODO: Will not work with ISF 
     235             N_out = N_out + 1 
     236             h_out(N_out) = e3t_n(ji,jj,jk) !Parent grid scale factors. Could multiply by e1e2t here instead of division above 
     237           ENDDO 
     238           IF (N_in > 0) THEN 
     239             h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
     240! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly 
     241!             IF abs(h_diff > 1.e-8) THEN 
     242!               N_in = N_in + 1 
     243!               h_in(N_in) = h_diff 
     244!               tabin(N_in,:) = tabin(N_in-1,:) 
     245             IF (h_diff < 0) THEN 
     246             print *,'CHECK YOUR bathy T points ...',ji,jj,h_diff,sum(h_in(1:N_in)),sum(h_out(1:N_out)) 
     247             print *,'Nval = ',N_out,mbathy(ji,jj) 
     248             print *,'BATHY = ',gdepw_0(ji,jj,mbathy(ji,jj)+1),sum(e3t_0(ji,jj,1:mbathy(ji,jj))) 
     249             STOP 
     250!               N_out = N_out + 1 
     251!               h_out(N_out) = - h_diff 
     252             ENDIF 
     253             DO jn=n1,n2-1 
     254               CALL reconstructandremap(tabin(1:N_in,jn),h_in(1:N_in),tabres_child(ji,jj,1:N_out,jn),h_out(1:N_out),N_in,N_out) 
     255             ENDDO 
     256          ENDIF 
     257         ENDDO 
     258         ENDDO 
     259          
     260! WARNING : 
     261! tabres replaced by tabres_child in the following 
     262! k1 k2 replaced by 1 jpk 
     263! VERTICAL REFINEMENT END 
     264 
    200265         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 
    201266            ! Add asselin part 
    202             DO jn = n1,n2 
    203                DO jk=k1,k2 
     267            DO jn = n1,n2-1 
     268               DO jk=1,jpk 
    204269                  DO jj=j1,j2 
    205270                     DO ji=i1,i2 
    206                         IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN 
     271                        IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN 
    207272                           tsb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) &  
    208                                  & + atfp * ( tabres(ji,jj,jk,jn) & 
     273                                 & + atfp * ( tabres_child(ji,jj,jk,jn) & 
    209274                                 &             - tsn(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 
    210275                        ENDIF 
     
    214279            ENDDO 
    215280         ENDIF 
    216          DO jn = n1,n2 
    217             DO jk=k1,k2 
     281         DO jn = n1,n2-1 
     282            DO jk=1,jpk 
    218283               DO jj=j1,j2 
    219284                  DO ji=i1,i2 
    220                      IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN  
    221                         tsn(ji,jj,jk,jn) = tabres(ji,jj,jk,jn) * tmask(ji,jj,jk) 
     285                     IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN  
     286                        tsn(ji,jj,jk,jn) = tabres_child(ji,jj,jk,jn) * tmask(ji,jj,jk) 
    222287                     END IF 
    223288                  END DO 
     
    230295 
    231296 
    232    SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, before ) 
     297   SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
    233298      !!--------------------------------------------- 
    234299      !!           *** ROUTINE updateu *** 
    235300      !!--------------------------------------------- 
    236       INTEGER                               , INTENT(in   ) :: i1, i2, j1, j2, k1, k2 
    237       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres 
    238       LOGICAL                               , INTENT(in   ) :: before 
     301      INTEGER                                 , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2 
     302      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,2), INTENT(inout) :: tabres 
     303      LOGICAL                                 , INTENT(in   ) :: before 
    239304      ! 
    240305      INTEGER  ::   ji, jj, jk 
    241306      REAL(wp) ::   zrhoy 
     307! VERTICAL REFINEMENT BEGIN 
     308      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child 
     309      REAL(wp) :: h_in(k1:k2) 
     310      REAL(wp) :: h_out(1:jpk) 
     311      INTEGER :: N_in, N_out 
     312      REAL(wp) :: h_diff 
     313      REAL(wp) :: tabin(k1:k2) 
     314! VERTICAL REFINEMENT END 
    242315      !!--------------------------------------------- 
    243316      !  
    244317      IF( before ) THEN 
    245318         zrhoy = Agrif_Rhoy() 
    246          DO jk = k1, k2 
    247             tabres(i1:i2,j1:j2,jk) = zrhoy * e2u(i1:i2,j1:j2) * e3u_n(i1:i2,j1:j2,jk) * un(i1:i2,j1:j2,jk) 
    248          END DO 
    249       ELSE 
    250319         DO jk=k1,k2 
    251320            DO jj=j1,j2 
    252321               DO ji=i1,i2 
    253                   tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e2u(ji,jj) / e3u_n(ji,jj,jk) 
     322                  tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) * un(ji,jj,jk) 
     323                  tabres(ji,jj,jk,2) = umask(ji,jj,jk) * zrhoy * e2u(ji,jj) * e3u_n(ji,jj,jk)  
     324               END DO 
     325            END DO 
     326         END DO 
     327      ELSE 
     328! VERTICAL REFINEMENT BEGIN 
     329         tabres_child(:,:,:) = 0. 
     330          
     331         DO jj=j1,j2 
     332         DO ji=i1,i2 
     333           N_in = 0 
     334           DO jk=k1,k2 !k2=jpk of child grid 
     335             IF (tabres(ji,jj,jk,2) == 0) EXIT 
     336             N_in = N_in + 1 
     337             tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 
     338             h_in(N_in) = tabres(ji,jj,jk,2)/e2u(ji,jj) 
     339           ENDDO 
     340           N_out = 0 
     341           DO jk=1,jpk 
     342             IF (umask(ji,jj,jk) == 0) EXIT 
     343             N_out = N_out + 1 
     344             h_out(N_out) = e3u_n(ji,jj,jk) 
     345           ENDDO 
     346           IF (N_in * N_out > 0) THEN 
     347             h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
     348! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly 
     349             if (h_diff < 0.) then 
     350             print *,'CHECK YOUR BATHY ...' 
     351             stop 
     352!             else ! Extends with 0 
     353!             N_in = N_in + 1 
     354!             tabin(N_in) = 0. 
     355!             h_in(N_in) = h_diff 
     356             endif 
     357             CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out) 
     358          ENDIF 
     359         ENDDO 
     360         ENDDO 
     361          
     362! WARNING : 
     363! tabres replaced by tabres_child in the following 
     364! k1 k2 replaced by 1 jpk 
     365! remove division by fs e3u (already done) 
     366! VERTICAL REFINEMENT END 
     367 
     368         DO jk=1,jpk 
     369            DO jj=j1,j2 
     370               DO ji=i1,i2 
     371!Following line now replaced by division higher up in vertical refinement case 
     372!                  tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e2u(ji,jj) / e3u_n(ji,jj,jk) 
    254373                  ! 
    255374                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 
    256375                     ub(ji,jj,jk) = ub(ji,jj,jk) &  
    257                            & + atfp * ( tabres(ji,jj,jk) - un(ji,jj,jk) ) * umask(ji,jj,jk) 
     376                           & + atfp * ( tabres_child(ji,jj,jk) - un(ji,jj,jk) ) * umask(ji,jj,jk) 
    258377                  ENDIF 
    259378                  ! 
    260                   un(ji,jj,jk) = tabres(ji,jj,jk) * umask(ji,jj,jk) 
     379                  un(ji,jj,jk) = tabres_child(ji,jj,jk) * umask(ji,jj,jk) 
    261380               END DO 
    262381            END DO 
     
    271390      !!           *** ROUTINE updatev *** 
    272391      !!--------------------------------------------- 
    273       INTEGER :: i1,i2,j1,j2,k1,k2 
     392      INTEGER :: i1,i2,j1,j2,k1,k2,n1,n2 
    274393      INTEGER :: ji,jj,jk 
    275       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: tabres 
     394      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,2) :: tabres 
    276395      LOGICAL :: before 
    277396      !! 
    278397      REAL(wp) :: zrhox 
     398! VERTICAL REFINEMENT BEGIN 
     399      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child 
     400      REAL(wp) :: h_in(k1:k2) 
     401      REAL(wp) :: h_out(1:jpk) 
     402      INTEGER :: N_in, N_out 
     403      REAL(wp) :: h_diff 
     404      REAL(wp) :: tabin(k1:k2) 
     405! VERTICAL REFINEMENT END 
    279406      !!---------------------------------------------       
    280407      ! 
     
    284411            DO jj=j1,j2 
    285412               DO ji=i1,i2 
    286                   tabres(ji,jj,jk) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk) 
    287                END DO 
    288             END DO 
    289          END DO 
    290       ELSE 
    291          DO jk=k1,k2 
     413                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) * vn(ji,jj,jk) 
     414                  tabres(ji,jj,jk,2) = vmask(ji,jj,jk) * zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk)  
     415               END DO 
     416            END DO 
     417         END DO 
     418      ELSE 
     419! VERTICAL REFINEMENT BEGIN 
     420         tabres_child(:,:,:) = 0. 
     421          
     422         DO jj=j1,j2 
     423         DO ji=i1,i2 
     424           N_in = 0 
     425           DO jk=k1,k2 
     426             IF (tabres(ji,jj,jk,2) == 0) EXIT 
     427             N_in = N_in + 1 
     428             tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 
     429             h_in(N_in) = tabres(ji,jj,jk,2)/e1v(ji,jj) 
     430           ENDDO 
     431           N_out = 0 
     432           DO jk=1,jpk 
     433             IF (vmask(ji,jj,jk) == 0) EXIT 
     434             N_out = N_out + 1 
     435             h_out(N_out) = e3v_n(ji,jj,jk) 
     436           ENDDO 
     437           IF (N_in * N_out > 0) THEN 
     438             h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
     439! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly 
     440             if (h_diff < 0.) then 
     441             print *,'CHECK YOUR BATHY ...' 
     442             stop 
     443!             else ! Extends with 0 
     444!             N_in = N_in + 1 
     445!             tabin(N_in) = 0. 
     446!             h_in(N_in) = h_diff 
     447             endif 
     448             CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out) 
     449          ENDIF 
     450         ENDDO 
     451         ENDDO 
     452          
     453! WARNING : 
     454! tabres replaced by tabres_child in the following 
     455! k1 k2 replaced by 1 jpk 
     456! remove division by fs e3v (already done) 
     457! VERTICAL REFINEMENT END 
     458 
     459         DO jk=k1,jpk 
    292460            DO jj=j1,j2 
    293461               DO ji=i1,i2 
    294                   tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e1v(ji,jj) / e3v_n(ji,jj,jk) 
     462!Following line now replaced by division higher up in vertical refinement case 
     463!                  tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e1v(ji,jj) / e3v_n(ji,jj,jk) 
    295464                  ! 
    296465                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 
    297466                     vb(ji,jj,jk) = vb(ji,jj,jk) &  
    298                            & + atfp * ( tabres(ji,jj,jk) - vn(ji,jj,jk) ) * vmask(ji,jj,jk) 
     467                           & + atfp * ( tabres_child(ji,jj,jk) - vn(ji,jj,jk) ) * vmask(ji,jj,jk) 
    299468                  ENDIF 
    300469                  ! 
    301                   vn(ji,jj,jk) = tabres(ji,jj,jk) * vmask(ji,jj,jk) 
     470                  vn(ji,jj,jk) = tabres_child(ji,jj,jk) * vmask(ji,jj,jk) 
    302471               END DO 
    303472            END DO 
  • branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_user.F90

    r7761 r8135  
    346346   ! 1. Declaration of the type of variable which have to be interpolated 
    347347   !--------------------------------------------------------------------- 
    348    CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts/),tsn_id) 
     348   CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts+1/),tsn_id) 
    349349   CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts/),tsn_sponge_id) 
    350350 
    351    CALL agrif_declare_variable((/1,2,0/),(/2,3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),un_interp_id) 
    352    CALL agrif_declare_variable((/2,1,0/),(/3,2,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),vn_interp_id) 
    353    CALL agrif_declare_variable((/1,2,0/),(/2,3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),un_update_id) 
    354    CALL agrif_declare_variable((/2,1,0/),(/3,2,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),vn_update_id) 
     351   CALL agrif_declare_variable((/1,2,0,0/),(/2,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),un_interp_id) 
     352   CALL agrif_declare_variable((/2,1,0,0/),(/3,2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),vn_interp_id) 
     353   CALL agrif_declare_variable((/1,2,0,0/),(/2,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),un_update_id) 
     354   CALL agrif_declare_variable((/2,1,0,0/),(/3,2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),vn_update_id) 
    355355   CALL agrif_declare_variable((/1,2,0/),(/2,3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),un_sponge_id) 
    356356   CALL agrif_declare_variable((/2,1,0/),(/3,2,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),vn_sponge_id) 
Note: See TracChangeset for help on using the changeset viewer.