      program proj
      implicit none
      integer    i,j
      integer nx,ny,loop,nt
C     nx, ny is the grid size dictated by the terrain database
      parameter (nx=899,ny=599,nt=1000000)
              
C     decide how often to output the solution and when to end 
C     the simulation
      integer  runtype, outputevery, finishtime
      parameter    (runtype=1, outputevery=500,finishtime=10000)
C     input parameters
      real    gammaf, rivspeed,coastspeed
      parameter    (gammaf=0.02,rivspeed=4.,coastspeed=10.)
      real    seatrav
      parameter    (seatrav=40.)
      real    v0,K0
      parameter    (v0=12.5,K0=3.5)
C     step length
      real    dt
      parameter    (dt=0.003)
C     dimensionless paramter R (cf. Reynolds number)
      real    R
      parameter    (R=gammaf/v0)
C     dimensionless parameters Mgamma (cf Mach number)
      real    Mgammar
      parameter    (Mgammar=rivspeed/gammaf)
      real    Mgammac
      parameter    (Mgammac=coastspeed/gammaf)
C     set altitude and latitude dependence
      real    hv,hn0,av,an0,latgrad
      parameter    (hv=1200.,hn0=1100.,av=100.,an0=50.,latgrad=0.)
      real    latgradriv,latgradE
      parameter    (latgradriv=0.6,latgradE=0.)
      double precision omp_get_wtime,t1,t2
                                                                                                                            
C      define variables
       real    unew(0:nx,0:ny),uold(0:nx,0:ny),a(0:nx,0:ny)
       real    v(0:nx,0:ny),n0(0:nx,0:ny),gamma(0:nx,0:ny)
       real    lat(0:nx,0:ny),lat2(0:nx,0:ny),lat3(0:nx,0:ny) 
       real    advelx(0:nx,0:ny),advely(0:nx,0:ny)
       real    veg(0:nx,0:ny),seas(0:nx,0:ny)
       real    cadvelx(0:nx,0:ny),cadvely(0:nx,0:ny)
       real    utime(0:nx,0:ny),utimecont(0:nx,0:ny)
       real    t,tm,dx,dy,pi,pr,counter,rad,multi
       real    dtheta,dphi,phi(0:nx),theta(0:ny),denom(0:ny)
       real    gradnx(0:nx,0:ny),gradny(0:nx,0:ny)
       real    nhatx(0:nx,0:ny),nhaty(0:nx,0:ny)
       real    vx(0:nx,0:ny),vy(0:nx,0:ny),cvx(0:nx,0:ny),cvy(0:nx,0:ny)
       real    dotprod(0:nx,0:ny),dotprodc(0:nx,0:ny),sq(0:nx,0:ny)

C      open files for outputting
       open(unit=90,file='ModelOutput.dat',form='unformatted') 

C      initialise variables
       multi=1.0
       do i=0,nx
       do j=0,ny
              utime(i,j)=0.0
              utimecont(i,j)=0.0
              veg(i,j)=100.0
       enddo
       enddo
       pi=4.0*atan(1.0)   
       rad=6378.1
       pr=0.0
       counter=1.0
       t=0.0
       tm=2.0
                                      
C      read in altitude and seas data and set up grid
       call geta (nx,ny,runtype,a,theta,phi)
       dtheta=(theta(0)-theta(ny))/ny
       dphi=(phi(nx)-phi(0))/nx
       call getseas (nx,ny,seas)
C      calculate latitude dependences
       call getlat (nx,ny,latgrad,latgradE,lat)
       call getlat2 (nx,ny,theta,lat2)
       call getlat3 (nx,ny,latgradriv,lat3)

C      calculate the spatial dependence of parameters
       call getv (nx,ny,a,lat,seas,hv,av,seatrav,v) 
       call getn0 (nx,ny,a,veg,lat,hn0,an0,n0)
       call getgamma (nx,ny,a,gamma)  

C      read in river and coasts vectors
       call getrvecs(nx,ny,advelx,advely)
       call getcvecs(nx,ny,cadvelx,cadvely)
       call getpolaric(nx,ny,runtype,dtheta,dphi,a,theta,phi,uold)
C      Quench Advection according to geography
       do j=0,ny 
       do i=0,nx
          cadvelx(i,j)=lat2(i,j)*cadvelx(i,j)*(0.5+
     +                        (-0.5*tanh((a(i,j)-1000)/100)))
          cadvely(i,j)=lat2(i,j)*cadvely(i,j)*(0.5+
     +                        (-0.5*tanh((a(i,j)-1000)/100)))
          advelx(i,j)=lat3(i,j)*advelx(i,j)
          advely(i,j)=lat3(i,j)*advely(i,j)
       enddo 
       enddo 
       

C     initialise the time since first arrival arrays
      do j=0,ny
        do i=0,nx
          if(uold(i,j).gt.0.01)then
            utime(i,j)=-1.
            utimecont(i,j)=-1
          endif
        enddo 
      enddo 
      
      do j=0,ny      
            denom(j)=rad*rad*sin(theta(j))*sin(theta(j))
      enddo 

C     Start timeing the loop
!$    t1 = omp_get_wtime()
      
C     begin time stepping loop
      do while (t.lt.(finishtime*gammaf))
!$OMP PARALLEL PRIVATE(i,j)
!$OMP do
      do j=0,ny
         do i=0,nx
C           calculate boundary conditions
            if(i.eq.  0)then
            unew(i,j)=(4*uold(1,j)-uold(2,j))/3
            elseif(i.eq.  nx)then
            unew(i,j)=(4*uold(nx-1,j)-uold(nx-2,j))/3
            elseif(j.eq.  0)then
            unew(i,j)=(4*uold(i,1)-uold(i,2))/3
            elseif(j.eq.  ny)then
            unew(i,j)=(4*uold(i,ny-1)-uold(i,ny-2))/3
          else
C          Calculate the direction of the advection
           gradnx(i,j)=(uold(i+1,j)-uold(i-1,j))/(2.0*dphi)
           gradny(i,j)=(uold(i,j+1)-uold(i,j-1))/(2.0*dtheta)
           sq(i,j)=sqrt(gradnx(i,j)**2+gradny(i,j)**2)
           if(sq(i,j).le.0.00000001)then
             vx(i,j)=0.0
             vy(i,j)=0.0
             cvx(i,j)=0.0
             cvy(i,j)=0.0     
           else
             nhatx(i,j)=-gradnx(i,j)/sq(i,j)
             nhaty(i,j)=-gradny(i,j)/sq(i,j)
           dotprod(i,j)=nhatx(i,j)*advelx(i,j)+nhaty(i,j)*advely(i,j)
           dotprodc(i,j)=nhatx(i,j)*cadvelx(i,j)+nhaty(i,j)*cadvely(i,j)
             vx(i,j)=sign(multi,dotprod(i,j))*advelx(i,j)
             vy(i,j)=sign(multi,dotprod(i,j))*advely(i,j)
             cvx(i,j)=sign(multi,dotprodc(i,j))*cadvelx(i,j)
             cvy(i,j)=sign(multi,dotprodc(i,j))*cadvely(i,j)
           endif

C      update the population density via the reaction-diffusion-
C      advection equation

       unew(i,j)=uold(i,j)+dt*((1.0/R)*(v(i,j)*
     + (uold(i,j+1)-2*uold(i,j)+uold(i,j-1))/(dtheta*dtheta)/(rad*rad)+
     + v(i,j)*(uold(i+1,j)-2*uold(i,j)+uold(i-1,j))/(dphi*dphi)/denom(j)
     + +v(i,j)*(uold(i,j+1)-uold(i,j-1))/(2*dtheta)*
     + cos(theta(j))/(rad*rad)/sin(theta(j))+       
     + (((v(i,j+1)-v(i,j-1))*(uold(i,j+1)-uold(i,j-1)))/
     + (4*dtheta*dtheta))/rad**2       
     + +(((v(i+1,j)-v(i-1,j))*(uold(i+1,j)-uold(i-1,j)))/
     + (4*dphi*dphi))/(denom(j)))+ 
     + gamma(i,j)*uold(i,j)*(1-(uold(i,j)/n0(i,j)))-
     + (Mgammar*(vx(i,j)*(uold(i+1,j)-uold(i-1,j))/
     + (2*dphi*sin(theta(j)))       
     + +vy(i,j)*(uold(i,j+1)-uold(i,j-1))/(2*dtheta))/rad
     + +Mgammac*(cvx(i,j)*(uold(i+1,j)-uold(i-1,j))/
     + (2*dphi*sin(theta(j)))
     + +cvy(i,j)*(uold(i,j+1)-uold(i,j-1))/
     + (2*dtheta))/rad))
                                                                              

          endif
        enddo 
      enddo 
!$OMP enddo
!$OMP END PARALLEL


C      update time
       t=t+dt
!$OMP PARALLEL PRIVATE(i,j)
!$OMP do
C       Do updates for use in the next timestep
        do j=0,ny
        do i=0,nx
        if (unew(i,j).lt.0.0)then
        uold(i,j)=0.0
        else
        uold(i,j)=unew(i,j)
        endif
        enddo
        enddo
!$OMP enddo
!$OMP END PARALLEL
C     calculate the time since first arrival (continuous)
      call getutime(nx,ny,t,uold,utimecont)
      
C     calculate the time since first arrival (discrete)
      if(t.ge.counter*gammaf*outputevery)then
      call getutime(nx,ny,t,uold,utime)
      counter=counter+1
      endif

C     end the time stepping loop
      enddo         
C     Stop timeing the loop
!$     t2 = omp_get_wtime()
!$    write(*,*) 'time=',t2-t1
      
C     finalise the time of arrival arrays
      do j=0,ny
      do i=0,nx
      if(utimecont(i,j).lt.0.0)then
      utimecont(i,j)=0.0
      endif
      if(utime(i,j).lt.0.0)then
      utime(i,j)=0.0
      endif
      enddo 
      enddo 
      
C     write out the solution
      write(90)utimecont
      
C     tidy up
      close(unit=90)
      stop                             
      end


      subroutine getutime(nx,ny,t,uold,utime)
C     Calculate the when the population reaches 
C     each grid cell
      implicit none
      integer    i,j,nx,ny
      real   t,uold(0:nx,0:ny),utime(0:nx,0:ny)
      do j=0,ny
      do i=0,nx
      if(uold(i,j).gt.0.01.and.utime(i,j).eq.  0.0)then
      utime(i,j)=t
      endif
      enddo 
      enddo 
      end 
            
      subroutine getpolaric(nx,ny,runtype,dtheta,dphi,a,theta,phi,uold)
C     initial conditions set up
      implicit none
      integer    nx,ny,i,j,runtype
      real    uold(0:nx,0:ny),dtheta,dphi
      real    theta(0:ny),phi(0:nx),a(0:nx,0:ny)
      
      do i=0,nx
         do j=0,ny
           phi(i)=phi(0)+i*dphi
           theta(j)=theta(0)-j*dtheta
        uold(i,j)=1.0*exp(-(1.0/0.00075)*((theta(j)-0.925)**2+
     + (phi(i)-0.7)**2))

           if(a(i,j).eq.  0.0)then
             uold(i,j)=0.0
            elseif(uold(i,j).le.0.01)then
             uold(i,j)=0.0
           endif                                        
         enddo 
      enddo 
      end 
                                       

                                                                                      

      Subroutine getv (nx,ny,a,lat,seas,hv,av,seatrav,v)
C     Calculation of diffusivity dependent on geography
      implicit none
      Integer    nx,ny,i,j
      Real    v(0:nx,0:ny),a(0:nx,0:ny),lat(0:nx,0:ny)
      Real    seas(0:nx,0:ny),av,hv,seatrav
      do j=0,ny
         do i=0,nx                             
            if (a(i,j).le.0.0)then
               v(i,j)=lat(i,j)*exp(-seas(i,j)/seatrav)                         
           else
               v(i,j)=1.0*lat(i,j)*(0.5+(-0.5*tanh((a(i,j)-hv)/av)))      
            endif
         enddo 
      enddo    
      return
      end 



   
      Subroutine getn0 (nx,ny,a,veg,lat,hn0,an0,n0)
C     Calculation of carrying capacity dependent on geography
      implicit none
      Integer    nx,ny,i,j
      Real    n0(0:nx,0:ny),a(0:nx,0:ny),lat(0:nx,0:ny)
      Real    veg(0:nx,0:ny),hn0,an0
      do j=0,ny
         do i=0,nx
          if (a(i,j).le.0.0)then
             n0(i,j)=0.0001
          else
             n0(i,j)=lat(i,j)*(0.5+(-0.5*tanh((a(i,j)-hn0)/an0)))
          endif
          if(n0(i,j).lt.0.0001)then
          n0(i,j)=0.0001
          endif
        enddo 
      enddo 
      return
      end 

      Subroutine getgamma (nx,ny,a,gamma)
C     Calculation of growth rate dependent on geography
      Integer    nx,ny
      Real    gamma(0:nx,0:ny),a(0:nx,0:ny)
      do j=0,ny
         do i=0,nx 
            if(a(i,j).le.0)then
            gamma(i,j)=0.0
            else
            gamma(i,j)=1.0
            endif 
         enddo 
      enddo    


      return
      end 

      subroutine getlat(nx,ny,latgrad,latgradE,lat)
C     Calculation of the latitude dependence
      implicit none
      integer    i,j,nx,ny
      real    lat(0:nx,0:ny),latgrad, EW(0:nx,0:ny),latgradE
      open(unit=1,file='EW.out',form='formatted')
      read(1,*) EW
      close(1)
       do j=0,ny
        do i=0,nx
        if(EW(i,j).eq.0)then
             lat(i,j)=max(0.0,1.-latgrad*real(j)/real(ny))
        else
             lat(i,j)=max(0.0,1.-latgradE*real(j)/real(ny))
        endif
        enddo 
      enddo 
      return
      end 

      subroutine getlat3(nx,ny,latgradriv,lat3)
C     Calculation of the latitude dependence for rivers
      implicit none
      integer    i,j,nx,ny
      real    lat3(0:nx,0:ny),latgradriv
      do j=0,ny
        do i=0,nx
             lat3(i,j)=max(0.0,1.-latgradriv*real(j)/real(ny))
        enddo 
      enddo 
      return
      end 

      subroutine getlat2(nx,ny,theta,lat)
C     Calculation of the latitude dependence for coasts
      implicit none
      integer    i,j,nx,ny
      real    lat(0:nx,0:ny),theta(0:ny)
      do j=0,ny
        do i=0,nx
             lat(i,j)=0.5+0.5*tanh((theta(j)-0.7)/0.05)
        enddo 
      enddo 
      return
      end


      Subroutine getrvecs (nx,ny,advelx,advely)
C     read in the vector of Rivers
      implicit none
      Integer    nx,ny
      Real    advelx(0:nx,0:ny),advely(0:nx,0:ny)
      open(unit=1,file='rvecsbigdirected2.out',form='formatted')
      read(1,*) advelx
      read(1,*) advely
      close(1)
      return
      end


      Subroutine getcvecs (nx,ny,cadvelx,cadvely)
C     read in the vector of coasts
      Implicit none
      Integer    nx,ny
      Real    cadvelx(0:nx,0:ny),cadvely(0:nx,0:ny)
      open(unit=1,file='cvecsnarrowbigdirected2.out',form='formatted')
      read(1,*) cadvelx
      read(1,*) cadvely
      close(1)
      return
      end

      Subroutine getseas (nx,ny,seas)
C     read in the shortest distance from land
      implicit none
      Integer    nx,ny
      Real    seas(0:nx,0:ny)
      open(unit=1,file='seasbig.dat',form='formatted')
      read(1,*) seas
      close(1)
      return
      end 
    

      Subroutine geta (nx,ny,runtype,a,theta,phi)
C     read in the altitude, latitude and longitude arrays
      implicit none
      Integer    nx,ny,i,j,runtype
      Real    a(0:nx,0:ny),theta(0:ny),phi(0:nx)
      real  pi
      pi=4.0*atan(1.0)
      if(runtype.eq.1)then
      open(unit=1,file='transformtemp.out',form='formatted')
      elseif (runtype.eq.2) then
      open(unit=1,file='transformextrasmooth.out',form='formatted')
      endif
      read(1,*) a
      read(1,*) phi
      read(1,*) theta
      close(1) 

C        change to radians from degrees
C        change to co-latitude from latitude
         do i=0,ny
            theta(i)=(Pi/2.0)-(theta(i)/180)*pi
         enddo 
C        change to radians from degrees
         do j=0,nx
            phi(j)=(phi(j)/180)*pi
         enddo      
      
      return
      end 
      
         
