! program onom_sem.f90 - file name is "o3o3_3.f90"
! J. Steppeler created
! J. Li modified 25 February 2019
!
! EH:
! program name changed to onom_sem (to be consistent)
! to compile on netcup server use these commands(see build proc makeO3o3_3):
! gfortran -std=gnu -fdefault-real-8 -cpp -fno-align-commons -g -c gkssubs2.f90
! gfortran -std=gnu -fdefault-real-8 -cpp -fno-align-commons -DNOGKS -g -o o3o3_3 o3o3_3.f90 gkssubs2.o
! When executing program o3o3_3 lot of "...dat" files are produced
! These "..dat" files are needed to produce gnuplot graphics in a second step
MODULE const
PARAMETER (ie=600,ib=6,it=7)
!PARAMETER (ie=20,ib=6,it=7)
real dx,dt,f,u00,h00,h00y
real wm2,wm1,w0,wp1,wp2
real x1(-ib:ie+ib+3)
real xp(-ib:ie+ib+3)
real a21 ,a22 ,a31 ,a32 ! coefficients for spectral to gp transform
real a2x1,a2x2,a3x1,a3x2 ! grid point values for derivatives of h2 and h3
real a2x0,a2x3,a3x0,a3x3 ! grid point values for derivatives of h2 and h3
real aglinp,aglinm ! Galerkin coefficient linear part, lower and upper interval
real a2xg0,a2xg3 ! Galerkin coefficient o2 part
real a3xg0,a3xg3 ! Galerkin coefficient o3 part
real aa2,aa3 ! coefficients for gp to spec transform
real fl1,fl2 ! coefficients for linear interpolation
real dx3 ! large grid interval
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! FOR NUMERICAL METHOD
! parameter ( imod = 1 )
parameter ( imod = 2 )
! Choice for the numerical method
! 1--control run, standard o4 differenceing
! 2--o3o3 dt = 2.5
! 3--o3o3 completely in spectral spave: not available
! 4-- not available
! 5--SEM3 dt = 1.5
! 6--SEM2
! cubic representation for h. should have no 0 space, d.h. only 2 dx wave.
! up to imod=3 in order
! imod = 5 for SEM3, but rather than Legendre polynomials linear, second and third-order componenta are used.
! There is probably an error in this. It has quite a bit of dispersion, and the CFL Condition is very large 7.8!!!.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!********************************************************
!! FOR PLOTTING
parameter ( line_number = 5 ) ! default = 5
parameter ( plot_together = 1 )
parameter ( matlab = 1 ) ! 0--close matlab files, 1--open matlab files
!!********************************************************
!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!! FOR INITIAL CONDITION
parameter ( initialconditionh = 0 ) ! defualt: 8; 0--peak solution, 4--smooth 4, 8--smooth 8, 1--constant = 4
parameter ( hconstant = 0 ) ! if hconstant = 1 then initialconditionh = 1 line_number = 1
!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
parameter ( condensation = 0 )
! CLOSE CONDENSATION
! 0--no condensation
! OPEN CONDENSATION
! 11--naive condensation method at third point
! 12--spline condensation method at third point
! 21--naive condensation method at every point
! 22--spline condensation method at every point
!!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
END MODULE const
! -------------------------------------------
MODULE field
USE const
real h (-ib:ie+ib,it) ! density
real ha2 (-ib:ie+ib,it) ! coefficient of h2
real ha3 (-ib:ie+ib,it) ! coefficient of h3
real hxx (-ib:ie+ib) ! coefficient of h2
real hta2(-ib:ie+ib) ! coefficient of h2 t
real hta3(-ib:ie+ib) ! coefficient of h2 t
real hlin(-ib:ie+ib) ! h linearly interpolated
real hho (-ib:ie+ib) ! high order part of h
real h2mid(0:ie) , htmid(0:ie), htxmid(0:ie), ht3(-ib:ie+ib)
real hmid (0:ie)
real hxmid(0:ie)
real hdiff(0:ie)
real hxxx(-ib:ie+ib) ! coefficient of h3 for SEM3
real hxxSE2 (-ib:ie+ib) ! coefficient of h2 for SEM2 ****Jinxi
END MODULE field
! ----------------------------------------------------------------------
SUBROUTINE initialize()
USE const
USE field
! USE gks_substitute ! not needed for real GKS/NCAR
!
!h2(x,dd) = (x**2-dd**2/4.)/2.;
!h3(x,dd) = (x*(x**2-dd**2/4.)/2.)/3. ! definition base
! characteristic polynomials
! initialise fields, low resolution and high resolution
dx = 1.
dx3 = 3. * dx
! weights for 4rth order differentiation
i = 6
! call mk_w(wm2,wm1,w0,wp1,wp2,x1,ie)! MOW_GC routione to compute weights for irregular grids.
! weights for regular resolution
wm2=8.33333358E-02;wm1= -0.666666687;w0=0.00000000;wp1=0.666666687;wp2= -8.33333358E-02
print*,wm2,wm1,w0,wp1,wp2,',wm1,wm2,w0,wp1,wp2,'
if (imod .eq. 1) then ! control run
dt = 2.0 ! stable stability limit control imod = 1
else if (imod .eq. 2) then
dt = 2.5 ! stable imod = 2
else if (imod .eq. 4) then
dt = 2.8 ! stable for Centred FD (imod = 4) unstable 2.9
else if (imod .eq. 5) then
dt = 1.5
else if (imod .eq. 3) then
dt = 0.5 !! smaller timesteps for scheme
else if (imod .eq. 6) then
dt = 0.5 !!
else if (imod .eq. 7) then
dt = 0.5 !!
end if
dt = 1.
print *, 2.8*.7,'2.8*.7'
u00 = 1.
h00 = 10.
x1(-ib) = -ib * dx !x1 is for u-positions
do i=-ib+1,ie+ib
x1(i) = x1(i-1) + dx
enddo
! the following is for imod=5, comment otherwise
if (imod == 5) then
do i = -ib, ie+ib, 3
xxmid = (x1(i)+x1(i+3))*.5
xx11 = -.447214*1.5 ; xx22=-xx11 ! Gauss Lobtto points
x1(i+1) = xxmid-.447214*1.5
x1(i+2) = xxmid+.447214*1.5
end do
!******************************************************************************
else if (imod == 6) then
do i = -ib, ie+ib, 2
xxmidSE2 = (x1(i)+x1(i+2))*.5
xx11SE2 = 0.0!-.447214*1.5 ; xx22=-xx11 ! Gauss Lobtto points
x1(i+1) = xxmidSE2 - xx11SE2
!x1(i+2) = xxmidSE2+1.0
end do
!******************************************************************************
end if
xp(-ib:ie+ib) = x1(-ib:ie+ib)*10./real(ie)-5.
pi = atan(1.) * 4.
h = 0.
if (initialconditionh == 0) then
h(:,1) = 0.
h(150,1) = 4.
h(149,1) = 2.*4./3.
h(151,1) = 2.*4./3.
!h(148,1) = 1.*4./3.
!h(152,1) = 1.*4./3.
!h(5,1) = 4.;
!h(4,1) = 2.;
!h(6,1) = 2.
else if (initialconditionh == 4) then
do i = -ib, ie+ib
h(i,1) = 4.*exp(-(x1(i)*dx/4.-x1(150)/4.*dx)**2)
enddo
else if (initialconditionh == 8) then
do i = -ib, ie+ib
h(i,1) = 4.*exp(-(x1(i)*dx/8.-x1(150)/8.*dx)**2)
enddo
else if (initialconditionh == 1) then
do i = -ib, ie+ib
h(i,1) = 4.0
enddo
end if
! periodic boundary conditions
h(- ib:0,1) = h(ie- ib:ie,1);
h(ie:ie+ ib,1) = h(0: ib,1)
do i = -ib, ie+ib
print*,i,h(i,1),'i,h(i,1) after bc'
enddo
! Open GKS
CALL OPNGKS
END SUBROUTINE initialize
! -------------------------------------------
SUBROUTINE finalize()
! USE gks_substitute ! not needed for real GKS/NCAR
!
! close GKS
CALL CLSGKS()
END SUBROUTINE finalize
! -------------------------------------------
PROGRAM onom_sem
use const
use field
! USE gks_substitute ! not needed for real GKS/NCAR
character(len = 80) :: filename_new_banyin1, form_banyin1
character(len = 80) :: filename_new_banyin2, form_banyin2
character(len = 80) :: filename_new_banyin3, form_banyin3
! Definitions for SEM3
eem(x) = .5-.5*x/1.5;eep(x)=.5+.5*x/1.5
eemx(x) = -.5/1.5;eepx(x)=.5/1.5
hh2(x) = (x**2-1.5**2)/2. ; hh3(x)=(x**3-x*1.5**2)/6.
hh2x(x) = (2.*x)/2. ; hh3x(x)=(3.*x**2-1.5**2)/6.
! Definitions for SEM2
eemSE2(x) = .5-.5*x/1.0;eepSE2(x)=.5+.5*x/1.0
eemSE2x(x) = -.5/1.0;eepSE2x(x)=.5/1.0
hh2SE2(x) = (x**2-1.0**2)/2.! ; hh3SE2(x)=(x**3-x*1.0**2)/6.
hh2SE2x(x) = (2.*x)/2.! ; hh3SE2x(x)=(3.*x**2-1.0**2)/6.
call initialize()
h(:,1)=0.
h(150,1) = 4.
! h(149,1) = 2.
! h(151,1) = 2.
do i=0,ie
h(i,1) = 4.*exp(-(x1(i)*dx/4.-x1(250)/4.*dx)**2)
enddo
!! Runge Kutta 4 time scheme
!nnt = 0
!ntper = 0
!nou2d = 0
!ncount = 1
! periodic boundary conditions
h(- ib:0,1) = h(ie- ib:ie,1);
h(ie:ie+ ib,1) = h(0: ib,1)
CALL SET ( 0. , 1. , 0. , 1. , -5.9, 5.9, -5. , 5. , 1 )
call gsln(3)
if (matlab == 1) then
select case (ie+ib)
case (0:9)
write(form_banyin1,'(i1)') ie+ib
case (10:99)
write(form_banyin1,'(i2)') ie+ib
case (100:999)
write(form_banyin1,'(i3)') ie+ib
case (1000:9999)
write(form_banyin1,'(i4)') ie+ib
case (10000:99999)
write(form_banyin1,'(i5)') ie+ib
end select
write(filename_new_banyin1,*) "initial",trim(form_banyin1),".txt"
open(0711,file=filename_new_banyin1)
do i = -ib, ie+ib
write(0711,*) i, h(i,1), h(i,2), h(i,3), h(i,4), h(i,5), h(i,6), h(i,7), maxval(abs(h(:,1) - h(:,6)))
end do
close(0711)
else if (matlab == 0) then
end if
! advection over a distance of 3000 dx
if (hconstant == 0) then
nnend=30000
nnend = 600 / dt ! Integration steps
nnend=150
! 1200/dt, 2/dt, 10/dt, 300/dt, 3000/dt, 30000/dt, 600/dt, 10
else if (hconstant == 1) then
nnend = 1
end if
iplot = 0
iplotrep = 0
do ntime = 1,nnend
! mass
fmass=0.
do i=0,ie-2
fmass = fmass+h(i,1) * (x1(i+1) - x1(i))
enddo
print *,ntime,'ntime',fmass,'mass',maxval(h(0:ie,1)),'max h'
h(:,2) = h(:,1)
h(:,4) = 0.
do irk = 1, 4
h(:,3) = 0.
u00 = 1.
h(-1,1) = h(1,1)
h(ie+1,1) = h(ie-1,1)
! EH: this is an empty loop - eliminate it
do i=-ib,ie
!h(i,1) = i**2!test
!h(i,1) = 0.
!print*,i,h(i,1),'i,h(i,1)'
enddo
!h(:,1) = 0.
!h(3,1) = 1.
! periodic boundary conditions
h(- ib:0,1) = h(ie- ib:ie,1);
h(ie:ie+ ib,1) = h(0: ib,1)
a21 = -1. / 2.
if (imod.eq.1)then ! control run
do i = -ib+2, ie+ib-3,3
! Mass in second order
! o2 amplitude for time derivative less accurate though more correct??
! standard differencing: non-conserving decomment for
! standard differencing (control)
! differencing from splines
h(i,3) = -u00*(wm2*h(i-2,1)+wm1*h(i-1,1)+wp1*h(i+1,1)+wp2*h(i+2,1))
h(i+1,3) = -u00*(wm2*h(i-2+1,1)+wm1*h(i-1+1,1)+wp1*h(i+1+1,1)+ &
wp2*h(i+2+1,1))
h(i+2,3) = -u00*(wm2*h(i-2+2,1)+wm1*h(i-1+2,1)+wp1*h(i+1+2,1)+ &
wp2*h(i+2+2,1))
enddo
else if (imod.eq.2) then ! o3o3 run
do i=-1,ie!test
!h(i,1)=i**4
enddo
do i = -ib+3, ie+ib-3,3
h(i,3) = -u00*(wm2*h(i-2,1)+wm1*h(i-1,1)+wp1*h(i+1,1)+wp2*h(i+2,1))
end do
do i = -ib+3, ie+ib-3,3
ha2(i,1) = -(-3.*((h(i+3,1)-h(i,1))/3.+(h(i+3,3)+h(i,3))/2.))*4./9.
!print*,I,ha2(I,1),'i,ha2',24*(i+1.5)
end do
do i = -ib+3, ie+ib-3,3
ha2(i,2)=(h(i,3)+h(i+3,3))*.5+ha2(i,1)*.5*1.5**2
enddo
do i = -ib+3, ie+ib-3,3
hhhx =(ha2(i,2)-ha2(i-3,2))/3.
hhhxp=(ha2(i+3,2)-ha2(i,2))/3.
ha3(i,1) =-(ha2(i+3,1)-ha2(i-3,1))/6.
!ha3(i,1)=0.
hlin1 = (h(i,1)*2./3. + h(i+3,1)/3.)
hlin2 = (h(i,1)*1./3. + h(i+3,1)*2./3.)
!print*,hlin1,hlin2,'hlin'
del1 = h(i+1,1)-hlin1
del2 = h(i+2,1)-hlin2
hhxx = (del1+del2)/(.5**2-1.5**2)
!print*,i,ha3(i,1),'i,ha3'
!print*,del1,del2,'del1,del2'
!print*,hhxx,'hhxx',12.*(i+1.5)**2
!ha3(i,1)= (-hhxx - (h(i+3,3)-h(i,3))/3.)!*1./1.5**2
!ha3(i,1)=-ha3(i,1)*6.
!print*,i,hhxx,'hhxx',(i+1.5)**2*i**2,(i+1.5)**2*i**2/hhxx
!print*,ha3(i,1),'ha3'
!print*,ha2(i,1),'ha2',(i+1.5)*24.
!ha3(i,1)=(hhhx+hhhxp)*6./1.5**2
end do
ha3(:,1)=ha3(:,1)!*.00
do i = -ib+3, ie+ib-3,3
! do i = -ib+9, ie+ib-9,3
!ha3(i,1)=ha3(i,1)-(h(-3,5)-2.*ha3(i,5)+ha3(i+3,5))*.25
!ha3(i,1)=(ha3(i-3,1)+ha3(i+3,1))*.5
enddo
!ha3(:,1)=0.
do i = -ib+3, ie+ib-3,3
h(i+1,3) = h(i,3)*2./3.+h(i+3,3)/3.+ha2(i,1)*.5*(.5**2-1.5**2)-ha3(i,1)*.5*(.5**2-1.5**2)*.5/3.
h(i+2,3) = h(i,3)/3.+h(i+3,3)*2./3.+ha2(i,1)*.5*(.5**2-1.5**2)+ha3(i,1)*.5*(.5**2-1.5**2)*.5/3.
!print*,i
!print*,h(i,3),h(i+1,3),h(i+2,3),-4.*i**3,-4.*(i+1)**3,-4.*(i+2)**3
end do
do i = -ib+3, ie+ib-3
!h(i,3)=h(i,3)-(h(i+1,3)-2.*h(i,3)+h(i-1,3))*.1
enddo
h(:,6)=h(:,3)
!stop
else if (imod.eq.3) then ! spectral differencing at main nodes
print*,'not available'
stop
else if (imod.eq.5) then ! SEM3
do i = -ib, ie+ib
!h(i,1) = x1(i)**3 !test
end do
! linear interpolated h
xx11 = -.447214*1.5 ; xx22=-xx11 ! Gauss lobtto points
!xx11 = -.5; xx22=.5 ! this recreates o3o3
hlin = h(:,1)
do i = -ib, ie+ib-3,3
hlin(i+1) = (hlin(i)*eem(xx11) + hlin(i+3)*eep(xx11))
hlin(i+2) = (hlin(i)*eem(xx22) + hlin(i+3)*eep(xx22))
end do
h(:,3) = 0.
h(:,6) = 0.
hxx(i) = 0.
do i = -ib+3, ie+ib-3,3
del1 = h(i+1,1)-hlin(i+1)
del2 = h(i+2,1)-hlin(i+2)
hxx(i) = (del1+del2)/(hh2(xx11)+hh2(xx22))
hxxx(i) = (del1-del2)/(hh3(xx11)-hh3(xx22))
end do
do i = -ib+3, ie+ib-3,3
hxl = (h(i+3,1)-h(i,1))/3.
hxlm = (h(i,1)-h(i-3,1))/3.
h(i,3) = -.5*(hxl-hxx(i)*1.5+hxxx(i)*hh3x(-1.5))-.5*(hxlm+hxx(i-3)*1.5+hxxx(i-3)*hh3x(-1.5))
end do
do i = -ib+3, ie+ib-3,3
!SE3
h(i+1,3) = -(h(i,1)*eemx(xx11)+h(i+3,1)*eepx(xx11)+hxx(i)*hh2x(xx11)+hxxx(i)*hh3x(xx11))
h(i+2,3) = -(h(i,1)*eemx(xx22)+h(i+3,1)*eepx(xx22)+hxx(i)*hh2x(xx22)+hxxx(i)*hh3x(xx22))
h(i+1,6) = -(h(i,1)*eemx(xx11)+h(i+3,1)*eepx(xx11)+hxx(i)*hh2x(xx11))
h(i+2,6) = -(h(i,1)*eemx(xx22)+h(i+3,1)*eepx(xx22)+hxx(i)*hh2x(xx22))
end do
!*************************************************************************************************************************
!*************************************************************************************************************************
else if (imod.eq.6) then ! SEM2
do i = -ib, ie+ib
!h(i,1) = x1(i)**3 !test
end do
! linear interpolated h **** Jinxi: changed the definition of xx11 and delete xx22
xx11SE2 = 0.0! xx11 = -.447214*1.5 ; xx22=-xx11 ! Gauss lobtto points
!xx11 = -.5; xx22=.5 ! this recreates o3o3
hlin = h(:,1)
do i = -ib, ie+ib-2,2 !**** Jinxi: changed to every 2 points and end at ie+ib-2
hlin(i+1) = (hlin(i)*eemSE2(xx11SE2) + hlin(i+2)*eepSE2(xx11SE2)) ! Jinxi: change the expressions of hlin
!hlin(i+2) = (hlin(i)*eem(xx22) + hlin(i+3)*eep(xx22))
end do
h(:,3) = 0.
h(:,6) = 0.
hxxSE2(i) = 0.
do i = -ib, ie+ib-2, 2 !**** Jinxi: changed to every 2 points and start at -ib+2 and end at ie+ib-2
del1 = h(i+1,1) - hlin(i+1)
!del2 = h(i+2,1) - hlin(i+2)
!hxx(i) = (del1 + del2) / (hh2(xx11) + hh2(xx22))
!hxxx(i) = (del1 - del2) / (hh3(xx11) - hh3(xx22))
hxxSE2(i) = - 1.0 / 1.0 * (2*h(i,1) - h(i-1,1) - h(i+1,1)) ! change the definition from Eq 13 of our paper
end do
do i = -ib, ie+ib-2,2 !**** Jinxi: changed to every 2 points and start at -ib+2 and end at ie+ib-2 and the expression of h(i,3)
!hxl = (h(i+3,1)-h(i,1))/3.
!hxlm = (h(i,1)-h(i-3,1))/3.
!************************************************************
!I copy the expression of h(i,3) from Equation (22)
!h(i,3) = (3.0/2.0)*(h(i+1,1)-h(i-1,1))/2.0 - 0.5*(h(i+2,1)-h(i-2,1))/4.0
!************************************************************
h(i,3) = -(3.0/2.0)*(h(i+1,1)-h(i-1,1))/2.0 + 0.5*(h(i+2,1)-h(i-2,1))/4.0
end do
do i = -ib, ie+ib-2,2 !**** Jinxi: changed to every 2 points and start at -ib+2 and end at ie+ib-2
h(i+1,3) = -(h(i,1)*eemSE2x(xx11SE2) + h(i+2,1)*eepSE2x(xx11SE2) + hxxSE2(i)*hh2SE2x(xx11SE2))
!h(i+1,3) = -(h(i,1)*eemx(xx11)+h(i+3,1)*eepx(xx11)+hxx(i)*hh2x(xx11)+hxxx(i)*hh3x(xx11))
!h(i+2,3) = -(h(i,1)*eemx(xx22)+h(i+3,1)*eepx(xx22)+hxx(i)*hh2x(xx22)+hxxx(i)*hh3x(xx22))
!h(i+1,6) = -(h(i,1)*eemSE2x(xx11SE2)+h(i+2,1)*eepSE2x(xx11SE2)+hxxSE2(i)*hh2SE2x(xx11SE2))
!h(i+1,6) = -(h(i,1)*eemx(xx11)+h(i+3,1)*eepx(xx11)+hxx(i)*hh2x(xx11))
!h(i+2,6) = -(h(i,1)*eemx(xx22)+h(i+3,1)*eepx(xx22)+hxx(i)*hh2x(xx22))
end do
!*******************************************************************************************************************************
!*******************************************************************************************************************************
else if (imod .eq. 7) then ! spectral differencing at main nodes, spectral time stepping
hlin = h(:,1)
!print*,h(:,1),'h(:,1)'
! linear interpolated h
do i = -ib, ie+ib-3,3
hlin(i+1) = (hlin(i)*2./3. + hlin(i+3)/3.)
hlin(i+2) = (hlin(i)*1./3. + hlin(i+3)*2./3.)
end do
h(:,3) = 0.
do i = -ib, ie+ib-3,3
hxl = (h(i+3,1)-h(i,1))/3.
h(i,3) = h(i,3)-.5*(hxl-ha2(i,1)*1.5+ha3(i,1)*3./4.)
h(i+3,3) = h(i+3,3)-.5*(hxl+ha2(i,1)*1.5+ha3(i,1)*3./4.)
end do
do i = -ib+3, ie+ib-3,3
hta2(i) = -(-3.*((h(i+3,1)-h(i,1))/3.+(h(i+3,3)+h(i,3))/2.))*4./9.
ha2(i,3) = hta2(i)
end do
do i = -ib+3, ie+ib-3,3
hta3(i) = (hta2(i+3)-hta2(i-3))/6.
ha3(i,3) = hta3(i)
end do
endif
! RK weights
if (irk.eq.1) ftim1 = dt * .5
if (irk.eq.2) ftim1 = dt * .5
if (irk.eq.3) ftim1 = dt * 1.
if (irk.eq.4) ftim1 = dt * .0
if (irk.eq.1) ftim2 = dt / 6.
if (irk.eq.2) ftim2 = dt / 3.
if (irk.eq.3) ftim2 = dt / 3.
if (irk.eq.4) ftim2 = dt / 6.
h(:,1) = h(:,3) * ftim1 + h(:,2)
h(:,4) = h(:,3) * ftim2 + h(:,4)
h(- ib:0,1) = h(ie- ib:ie,1);
h(ie:ie+ ib,1) = h(0: ib,1) ! per bc
if (imod .ne. 5) then
! ha2(:,1) = ha2(:,3) * ftim1 + h(:,2)
! ha2(:,4) = ha2(:,3) * ftim2 + h(:,4)
!! ha2(- ib:0,1) = ha2(ie- ib:ie,1);
!! ha2(ie:ie+ ib,1) = ha2(0: ib,1) ! per bc !ha2 for imod=2
! ha3(:,1) = ha3(:,3) * ftim1 + h(:,2)
! ha3(:,4) = ha3(:,3) * ftim2 + h(:,4)
! ha3(- ib:0,1) = ha3(ie- ib:ie,1);
! ha3(ie:ie+ ib,1) = ha3(0: ib,1) ! per bc !ha3 for imod=2
do i=0,ie,2
h(i,6) = h(i,1);
!******************************Jinxi
if (imod == 6) then
h(i+1,6) = h(i,1)*0.5+h(i+3,1)*0.5;
else
h(i+1,6) = h(i,1)*2./3.+h(i+3,1)/3.;
h(i+2,6) = h(i,1)/3.+h(i+3,1)*2./3.
end if
!*******************************
h(i,7) = h(i,3);
h(i+1,7) = (h(i,3)+h(i+3,3))*.5
h(i+1,7) = h(i+1,7)+hta2(i)*a21
end do
end if
end do ! end do irk=1,4
CALL CURVE ( xp(0:ie) , h(0:ie,1), ie+1 )
call frame
!CALL CURVE ( xp(0:ie) , h(0:ie,3), ie+1 )
! call frame
!CALL CURVE ( xp(0:ie) , ha3(0:ie,1), ie+1 )
! call frame
h(:,1) = h(:,2) + h(:,4)
h(:,7) = h(:,7) * 5.
ipl = 300. / dt
iplot = iplot + 1
if (plot_together == 0) then
if (iplot.eq.nnend/line_number) then
CALL SET ( 0. , 1. , 0. , 1. , -5.9 , 5.9, -5. , 5. , 1 )
call gsln(1)
if(ntime.ne.1) CALL CURVE ( xp(0:ie) , h(0:ie,1)-iplotrep*.9-0.9, ie+1 )
iplotrep=iplotrep+1
iplot=0
endif
else if (plot_together == 1) then
if (iplot.eq.nnend/line_number) then
CALL SET ( 0. , 1. , 0. , 1. , -5.9 , 5.9, -5. , 5. , 1 )
call gsln(3) ! line style
iplot=0
endif
endif
! Plotting using NCAR_GRAFICS is used, which is a freely avalable package. Any other grafics could be used as well. .
CALL SET ( 0. , 1. , 0. , 1. , -5.9 , 5.9, -5. , 5. , 1 )
! CALL HALFAX ( 10,2 , 10 , 2 , 0., 0., 1, 1 )
if(ntime.eq.1)call gsln(1)
if(ntime.eq.1) CALL CURVE ( xp(0:ie) , h(0:ie,1), ie+1 )
if(ntime.eq.1) print*,'plot'
if(ntime.eq.300)call gsln(2)
if(ntime.eq.300) CALL CURVE ( xp(0:ie) , h(0:ie,1), ie+1 )
if(ntime.eq.300) print*,'plot'
if(ntime.eq.29999)call gsln(3)
!if(ntime.eq.29999) CALL CURVE ( xp(0:ie) , h(0:ie,1), ie+1 )
if(ntime.eq.29999) print*,'plot'
if (matlab == 1) then ! output for matlab plotting
if ((ntime.eq.nnend/line_number)) then
write(filename_new_banyin2,*) "middle101.txt"
open(0712,file=filename_new_banyin2)
do i = -ib, ie+ib
write(0712,*) i, h(i,1), h(i,2), h(i,3), h(i,4), h(i,5), h(i,6), h(i,7), maxval(abs(h(:,1) - h(:,6)))
end do
close(0712)
else if ((ntime.eq.2*nnend/line_number)) then
write(filename_new_banyin2,*) "middle202.txt"
open(0712,file=filename_new_banyin2)
do i = -ib, ie+ib
write(0712,*) i, h(i,1), h(i,2), h(i,3), h(i,4), h(i,5), h(i,6), h(i,7), maxval(abs(h(:,1) - h(:,6)))
end do
close(0712)
else if ((ntime.eq.3*nnend/line_number)) then
write(filename_new_banyin2,*) "middle303.txt"
open(0712,file=filename_new_banyin2)
do i = -ib, ie+ib
write(0712,*) i, h(i,1), h(i,2), h(i,3), h(i,4), h(i,5), h(i,6), h(i,7), maxval(abs(h(:,1) - h(:,6)))
end do
close(0712)
else if ((ntime.eq.4*nnend/line_number)) then
write(filename_new_banyin2,*) "middle404.txt"
open(0712,file=filename_new_banyin2)
do i = -ib, ie+ib
write(0712,*) i, h(i,1), h(i,2), h(i,3), h(i,4), h(i,5), h(i,6), h(i,7), maxval(abs(h(:,1) - h(:,6)))
end do
close(0712)
else if ((ntime.eq.300)) then
write(filename_new_banyin2,*) "middle606.txt"
open(0712,file=filename_new_banyin2)
do i = -ib, ie+ib
write(0712,*) i, h(i,1), h(i,2), h(i,3), h(i,4), h(i,5), h(i,6), h(i,7), maxval(abs(h(:,1) - h(:,6)))
end do
close(0712)
else if ((ntime.eq.5*nnend/line_number)) then
write(0712,*) i, h(i,1), h(i,2), h(i,3), h(i,4), h(i,5), h(i,6), h(i,7), maxval(abs(h(:,1) - h(:,6)))
endif
else if (matlab == 0) then
end if
enddo ! end do ntime=...
CALL CURVE ( xp(0:ie) , h(0:ie,1), ie+1 )
call frame
call frame
if (matlab == 1) then ! output for plotting by matlab
select case (ie+ib)
case (0:9)
write(form_banyin3,'(i1)') ie+ib
case (10:99)
write(form_banyin3,'(i2)') ie+ib
case (100:999)
write(form_banyin3,'(i3)') ie+ib
case (1000:9999)
write(form_banyin3,'(i4)') ie+ib
case (10000:99999)
write(form_banyin3,'(i5)') ie+ib
end select
write(filename_new_banyin3,*) "result",trim(form_banyin3),".txt"
open(0713,file=filename_new_banyin3)
do i = -ib, ie+ib
write(0713,*) i, h(i,1), h(i,2), h(i,3), h(i,4), h(i,5), h(i,6), h(i,7), maxval(abs(h(:,1) - h(:,6)))
end do
close(0713)
else if (matlab == 0) then
end if
CALL SET ( 0. , 1. , 0. , 1. , -5.9 , 5.9, -5. , 5. , 1 )
call gsln(2)
CALL CURVE ( xp(-ib:ie+ib) , h(-ib:ie+ib,6), ie+1+2*ib )
!
! close GKS
CALL CLSGKS()
! CALL finalize()
END PROGRAM onom_sem