/home/cholod/wd/peps/ROMS/romsagrif_r215/Roms_tools/Roms_Agrif/set_weights.F

Go to the documentation of this file.

#include "cppdefs.h"
#ifdef SOLVE3D
c--# define TEST_WEIGHTS

 
      subroutine set_weights
      implicit none
# include "param.h"
# include "scalars.h"
      integer i,j, iter
# ifdef M2FILTER_POWER
      real p,q,r, scale
# else
      real arg
# endif
# ifdef TEST_WEIGHTS
      real zeta(4), DUon, Zt_avg1, DU_avg2, cff1,cff2,cff3
# endif
      real*QUAD sum, shft, cff

# ifdef AGRIF
#  include "zoom.h"
      real t1,t2,t3,t4,t11,t12
      real cff10,cff1m
# endif     
      
      nfast=0
      do i=1,2*ndtfast                  ! Reset both sets of
        weight(1,i)=0.                  ! weights to all-zeros.
        weight(2,i)=0.
      enddo
 
# ifdef M2FILTER_POWER
!                                        ! Possible settings of
! Power-function shaped filters          ! parameters to yield the
!                                        ! second-order accuracy:
!      f(xi)=xi^p*(1-xi^q)-r*xi          !
!                                        ! p   q   r,I3=3I2-2  I2=1
! where xi=scale*i/ndtfast; scale, p,    ! -------------------------
! q, r, and normalization are chosen to  !  2.  1.    0.1181   0.169
! yield correct zeroth- (normalization), !  2.  2.    0.1576   0.234
! first- (consistency), and second-order !  2.  3.    0.1772   0.266
! moments, resulting in second-order     !  2.  4.    0.1892   0.284
! temporal accuracy for time-averaged    !  2.  5.    0.1976   0.296
! barotropic motions resolved by         !  2.  6.    0.2039   0.304
! baroclinic time step.                  !  2.  8.    0.2129   0.314
!
      p=2.0
      q=4.0
      r=0.25
      
      scale=(p+1.)*(p+q+1.)/((p+2.)*(p+q+2.)*float(ndtfast))
      do iter=1,16
        nfast=0
        do i=1,2*ndtfast
          cff=scale*float(i)
          weight(1,i)=cff**p - cff**(p+q) - r*cff
          if (i.gt.ndtfast .and. weight(1,i).lt.0.) then
            weight(1,i)=0.
          else
            nfast=i                  ! Find center of gravity of
          endif                      ! the primary weighting shape
        enddo                        ! function and iteratively 
        sum=0.                       ! adjust "scale" to place the 
        shft=0.                      ! centroid exactly at
        do i=1,nfast                 ! "ndtfast".
          sum=sum+weight(1,i)
          shft=shft+weight(1,i)*float(i)
        enddo
        scale=scale * shft/(sum*float(ndtfast))
c**     write(*,*) shft/sum, ndtfast
      enddo
# elif defined M2FILTER_FLAT
      cff=1./float(ndtfast)               ! FLAT filter
      do i=1,2*ndtfast                    ! with width of 2
        arg=cff*float(i-ndtfast)
        if (2.*abs(arg) .lt. 1.) then
          weight(1,i)=1.
          nfast=i
        endif
      enddo
            
# elif defined M2FILTER_COSINE
      cff=pi/float(ndtfast)               ! cos**2 shaped
      do i=1,2*ndtfast                    ! filter 
        arg=cff*float(i-ndtfast)
        if (2.*abs(arg) .lt. pi) then
          weight(1,i)=(cos(arg))**2
          nfast=i
        endif
      enddo
# else
      cff=pi/float(ndtfast)
      do i=1,2*ndtfast
        arg=cff*float(i-ndtfast)
        if (2.*abs(arg) .lt. pi) then
!          weight(1,i)=???
          weight(1,i)=(cos(arg))**2       ! cos**2 shaped
          nfast=i                         ! filter (DEFAULT).
        endif
      enddo
# endif
!
! Post-processing of primary weights: Although it is assumed that
! the initial settings of the primary weights has its center of
! gravity "reasonably close" to ndtfast, it may be not so according
! to the discrete rules of integration. The following procedure is
! designed to put the senter of gravity exactly to ndtfast by
! computing mismatch "ndtfast-shft" and applying basically an
! upstream advection of weights to eliminate the mismatch
! iteratively. Once this procedure is complete primary weights are
! normalized.
!
      do iter=1,ndtfast                    ! Find center of gravity
        sum=0.                             ! of the primary weights
        shft=0.                            ! and subsequently
        do i=1,nfast                       ! calculate the mismatch
          sum=sum+weight(1,i)              ! to be compensated
          shft=shft+float(i)*weight(1,i)
        enddo
        shft=shft/sum
        cff=float(ndtfast)-shft
# ifdef TEST_WEIGHTS
        MPI_master_only write(*,'(A,2I4,2F22.18)')
     &  'centering weights:', iter, ndtfast, shft, cff
# endif
        if (cff .gt. 1.) then              ! Apply advection step
          nfast=nfast+1                    ! using either whole, or
          do i=nfast,2,-1                  ! fractional shifts.
            weight(1,i)=weight(1,i-1)
          enddo                            ! Note that none of
          weight(1,1)=0.                   ! the four loops here
        elseif (cff .gt. 0.) then          ! is reversible.
          sum=1.-cff
          do i=nfast,2,-1
            weight(1,i)=sum*weight(1,i)+cff*weight(1,i-1)
          enddo
          weight(1,1)=sum*weight(1,1)
        elseif (cff .lt. -1.) then
          nfast=nfast-1
          do i=1,nfast,+1
            weight(1,i)=weight(1,i+1)
          enddo
          weight(1,nfast+1)=0.
        elseif (cff .lt. 0.) then
          sum=1.+cff
          do i=1,nfast-1,+1
            weight(1,i)=sum*weight(1,i)-cff*weight(1,i+1)
          enddo
          weight(1,nfast)=sum*weight(1,nfast)
        endif
      enddo
                                           ! Set SECONDARY weights
      do j=1,nfast                         ! assuming that backward
        cff=weight(1,j)                    ! Euler time step is used
        do i=1,j                           ! for free surface. NOTE
          weight(2,i)=weight(2,i)+cff      ! that array weight(2,i)
        enddo                              ! is assumed to have all-
      enddo                                ! zero status at entry in
      sum=0.                               ! this segment of code.
      cff=0.
      do i=1,nfast                         ! Normalize both sets of
        sum=sum+weight(1,i)                ! weights.
        cff=cff+weight(2,i)
      enddo
      sum=1./sum
      cff=1./cff
      do i=1,nfast
        weight(1,i)=sum*weight(1,i)
        weight(2,i)=cff*weight(2,i)
      enddo
            
# ifdef AGRIF
       
      weight2 = 0.
       
      weight2(nfast,:) = weight(1,:)
      
            
      do j=nfast-1,1,-1
      
      cff10 = weight2(j+1,0)
      cff1m = weight2(j+1,j+1)
      
      t1 = (1+cff10*(-1.+(j+1)))*ndtfast
     &        +(-1+cff10)*nfast
      t2 = (cff10-cff1m)*(j+1)*ndtfast-
     & (-1.+cff10)*(-1.+cff1m+cff1m*(j+1))*nfast
     
      t3 = t1/t2
      
      t4 = (1.-t3*(1-cff1m))/(1.-cff10)
           
      do i=0,j
        weight2(j,i) = t3*weight2(j+1,i)+t4*weight2(j+1,i+1)
      enddo
                      
      weight2(j,j+1:) = 0.
 
      enddo
                   
# endif
       
 
# ifdef TEST_WEIGHTS
      MPI_master_only write(*,'(/1x,A,I3,4x,A,I3/4x,A,2(12x,A))')
     &               'Time Splitting weights: ndtfast =', ndtfast,
     &               'nfast =',  nfast,   'primary',  'secondary',
     &               'accumulated-to-current-step'
      cff=0.
      cff1=0.
      cff2=0.
      cff3=0.
      sum=0.
      shft=0.
      do i=1,nfast
        cff=cff   + weight(1,i)
        cff1=cff1 + weight(1,i)*float(i)
        cff2=cff2 + weight(1,i)*float(i*i)
        cff3=cff3 + weight(1,i)*float(i*i*i)
        sum=sum+weight(2,i)
        shft=shft + float(i)*weight(2,i)
        MPI_master_only write(*,'(I3,4F19.16)')
     &    i, weight(1,i), weight(2,i), cff, sum
      enddo
      cff1=cff1/float(ndtfast)
      cff2=cff2/float(ndtfast*ndtfast)
      cff3=cff3/float(ndtfast*ndtfast*ndtfast)
      shft=(shft-0.5)/float(ndtfast)

      MPI_master_only write(*,'(3(/A,2F14.10),F14.10/A,2I4,F8.4)')
     &            '  Checking integrals (must be 1, 1)', cff,  sum,
     &            '        centroid  (must be 1, ~0.5)', cff1, shft,
     &            ' second and third moments (>~1, ~1)', cff2, cff3,
     &                        cff3-(3*cff2-2.),
     &            '     ndtfast, nfast, nfast/ndtfast ',
     &                  ndtfast, nfast, float(nfast)/float(ndtfast)

      if (cff2 .lt. 1.001) then
        MPI_master_only write(*,'(/1x,A/)')
     &  'WARNING: unstable weights, reduce parameter "r".'
      endif
 
      MPI_master_only write(*,'(/A/A,1x,A,15x,A,7x,A/)')
     &  'Checking consistency of primary and secondary weights...',
     &                     'step',  'zeta_avr1', 'dt*divVBAR_avr2',
     &                     'zeta_avr1-dt*divVBAR_avr2'
 
      do j=1,nfast              !--> Loop over realizations
        next_kstp=1
        zeta(1)=0.              ! set initial conditions
        zeta(2)=0.              ! for each realization.
        do iif=1,nfast

 
          if (iif.eq.j) then
            DUon=1.
          else                  ! Set up delta function over
            DUon=0.             ! over ensamble of individual
          endif                 ! realizations j.
 
#  undef FORW_BAK
#  ifdef FORW_BAK

          kstp=knew
          knew=kstp+1
          if (knew.gt.4) knew=1
          call step2d_FB_imitator (zeta, DUon, Zt_avg1, DU_avg2)

#  else
          kstp=next_kstp
          knew=3
          call step2d_imitator (zeta, DUon, Zt_avg1, DU_avg2)
 
          knew=3-kstp
          next_kstp=knew
          call step2d_imitator (zeta, DUon, Zt_avg1, DU_avg2)
#  endif          
        enddo
        write(*,'(I4,2F24.18,F24.20)') j, Zt_avg1, dt*DU_avg2,
     &                                     Zt_avg1-dt*DU_avg2
      enddo
      stop
                                       ! Return values
      iif=1                            ! of time indices
      kstp=1                           ! back to their
      knew=1                           ! default initial
      return                           ! values, as set
      end                              ! by "init_scalars".
 
 
      subroutine step2d_FB_imitator (zeta, DUon, Zt_avg1, DU_avg2)
      implicit none
      real zeta(4), DUon, Zt_avg1, DU_avg2, zeta_new, cff1,cff2
#  include "param.h"
#  include "scalars.h"

      zeta(knew)=zeta(kstp) + dtfast*DUon

      cff1=weight(1,iif)
      cff2=weight(2,iif)
      if (FIRST_2D_STEP) then
        Zt_avg1=cff1*zeta(knew)
        DU_avg2=cff2*DUon
      else
        Zt_avg1=Zt_avg1+cff1*zeta(knew)
        DU_avg2=DU_avg2+cff2*DUon
      endif

      return
      end

 
      subroutine step2d_imitator (zeta, DUon, Zt_avg1, DU_avg2)
      implicit none
      real zeta(3), DUon, Zt_avg1, DU_avg2, zeta_new
#  include "param.h"
#  include "scalars.h"
 
      if (knew.eq.3) then
        krhs=kstp
        PREDICTOR_2D_STEP=.true.
      else
        krhs=3
        PREDICTOR_2D_STEP=.false.
      endif
!
! Imitate fast time averaging procedure. This module MUST BE
! consistent (to every minor detail, including startup process)
! with the fast-time-averaging procedure in step2d_tile (see
! file "step2D_LF_AM3.F").
!
 
      if (PREDICTOR_2D_STEP) then
        if (FIRST_2D_STEP) then
          Zt_avg1=0.
          DU_avg2=0.
        else
          Zt_avg1=Zt_avg1 + weight(1,iif-1)*zeta(krhs)
        endif
      else
        DU_avg2=DU_avg2 + weight(2,iif)*DUon
      endif
!
! Imitate the actual fast-time stepping for free-surface elevation.
! This module must be consistent with the time stepping in step2d,
! see file "step2d_LF_AM3.F".
!
      if (PREDICTOR_2D_STEP) then
 
      else
        zeta(knew)=zeta(kstp) + dtfast*DUon
      endif
!
! Finalize computation of barotropic mode averages.
!
      if (iif.eq.nfast .and. .not.PREDICTOR_2D_STEP) then
        Zt_avg1=Zt_avg1 + weight(1,iif)*zeta(knew)
      endif
 
# else
      MPI_master_only write(*,'(/1x,A,I3,4x,A,I3/)')
     &   'Time splitting: ndtfast =', ndtfast, 'nfast =', nfast
 
# endif   /* TEST_WEIGHTS */
      return
      end
 
#else
      subroutine set_weights_empty
      end
#endif     /* SOLVE3D */
 

---