! 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