!Model of migration of krill in Scotia sea (date 2000)
!The first approximation
!Pleshanov D. A. 05.2015
program Migrationofkrillmodel

use ifport
implicit none

integer nDay,ix,jy,kz,i1,n1,nt,i,j,iday,k,iitime,zz
integer ix1,ix2,jy1,jy2,kz1,kz2
integer n2, j2, jp, i2, n3, c
integer p,pt,pp,jj

integer*2 s1,s2
real*4 rr

!for trilinear interpolation
real U11,U12,U1,V11,V12,V1
real U21,U22,U2,V21,V22,V2
real U31,U32,U3,V31,V32,V3
real U41,U42,U4,V41,V42,V4
real U51,U52,U5,V51,V52,V5
real U61,U62,U6,V61,V62,V6
real U71,U72,U7,V71,V72,V7
real U81,U82,U8,V81,V82,V8
real iix1,iix2,jjy1,jjy2,kkz1,kkz2
real U,V,Ui,Vi
real fi,sigma1,sigma2,a1,a2,rd
real*4 rann

real iitime2,kzz,zm
real dt,begX,begY,dX,dY,dZ,R,pi,A,b,Cr,nt1,crz
real pointX,pointY,pointZ,dpx,dpy,dpz
real pointXm,pointYm,dpointXm,dpointYm
real pointXrand,pointYrand,pointZrand

character*30 infileU, infileV, inkrill, outkrill,crout, crout2
character*1 inU,inV
character*2 Uname2, Vname2,cr3
character*3 Uname1,Vname1,cr2
character*4 frmt, cr1

include 'krill.c'   !file with input variables
dimension U(ix,jy,kz),V(ix,jy,kz),Cr(n2,i2),crz(n2,i2)

a1=.001   !for random U and V
a2=.001   !for random Z
pi=3.141593
!R=6371.032 !km
R=6371032.0 !m
!input and output file names
data inU/'u'/
data inV/'v'/
data frmt/'.dat'/
data Uname1/'U00'/
data Vname1/'V00'/
data Uname2/'U0'/
data Vname2/'V0'/
data cr1/'Cr00'/
data cr2/'Cr0'/
data cr3/'Cr'/

print*, 'Model of migration of krill in Scotia sea (date 2000)'
print*, 'The first approximation'
print*, ' '
print*, 'Read input variables'

!open and read file with another variables
open(1,file='Krill.ini',status='old')
read(1,*) nDay
read(1,*) dt
read(1,*) begX
read(1,*) begY
read(1,*) dX
read(1,*) dY
read(1,*) dz
read(1,*) inKrill
close(1)

print*, 'Read input file with krill'
open(1,file=inKrill,status='old')   !read input krill data file
open(2,file='krill00.dat')

do j=1,n1,1
    read(1,*) (Cr(j,i),i=1,i1)
enddo
do j=1,n1
    Cr(j,4)=(pi/2.)*(1.-(2.*(Cr(j,3)-1.)/zm))
enddo

!write input file with krill on 0-250m
print*, 'Write input file'
jj=2
do j2=2,200,4
    do j=1,n1,1
        jp=(jj-1)*n1+j
        Cr(jp,1)=Cr(j,1)
        Cr(jp,2)=Cr(j,2)
        Cr(jp,3)=j2*1.
        Cr(jp,4)=(pi/2.)*(1.-(2.*Cr(jp,3)/zm))
    enddo
    jj=jj+1
enddo

do j=1,n2
    write(2,*) (Cr(j,i),i=1,4)
enddo

close(1)
close(2)

print*, 'Calculation began'
print*, 'Write output files:'
do iday=1,nDay-1    !an external step at a time
    !write names of input files with U and V
    !write names of output files
    if (iday.ge.1.and.iday.le.9) then
        write(inFileU,'(a3,i1,a4)') Uname1,iday,frmt
        write(inFileV,'(a3,i1,a4)') Vname1,iday,frmt
        write(Crout,'(a4,i1,a4)') cr1,iday,frmt
    else
        if (iday.ge.10.and.iday.le.99) then
            write(inFileU,'(a2,i2,a4)') Uname2,iday,frmt
            write(inFileV,'(a2,i2,a4)') Vname2,iday,frmt
            write(Crout,'(a3,i2,a4)') cr2,iday,frmt
        else
            write(inFileU,'(a1,i3,a4)') inU,iday,frmt
            write(inFileV,'(a1,i3,a4)') inV,iday,frmt
            write(Crout,'(a2,i3,a4)') cr3,iday,frmt
        endif
    endif
    open(1,file=inFileU,status='old')
    open(2,file=inFileV,status='old')

    do k=1,kz,1     !read input massive with U and V
        do j=1,jy,1
            read(1,*) (U(i,j,k),i=1,ix,1)
            read(2,*) (V(i,j,k),i=1,ix,1)
        enddo
    enddo
    close(1)
    close(2)

    !internal step at a time
    do iitime=0,nt,dt
        do i=1,n2   !for markers
            pointX=Cr(i,1)
            pointY=Cr(i,2)
            pointZ=Cr(i,3)
            fi=Cr(i,4)

            !write index
            dpX=(pointX-begX)/dx
            dpY=(pointY-begY)/dy
            dpZ=pointZ/dz

            ix1=int(dpX)+1
            if (ix1.ge.475) goto 404
            ix2=ix1+1
            iix1=real(ix1)-1.
            iix2=real(ix2)-1.

            jy1=int(dpY)+1
            if (jy1.ge.114) goto 404
            jy2=jy1+1
            jjy1=real(jy1)-1.
            jjy2=real(jy2)-1.

            kz1=int(dpZ)+1
            kz2=kz1+1
            kkz1=real(kz1)-1.
            kkz2=real(kz2)-1.

            !trilinear interpolation in point
            !/((iix2-iix1)*(jjy2-jjy1)*(kkz2-kkz1))
            U11=U(ix1,jy1,kz1)
            U12=(iix2-dpX)*(jjy2-dpY)*(kkz2-dpz)
            U1=U11*U12
            U21=U(ix1,jy1,kz2)
            U22=(iix2-dpX)*(jjy2-dpY)*(dpz-kkz1)
            U2=U21*U22
            U31=U(ix1,jy2,kz1)
            U32=(iix2-dpX)*(dpY-jjy1)*(kkz2-dpZ)
            U3=U31*U32
            U41=U(ix1,jy2,kz2)
            U42=(iix2-dpX)*(dpY-jjy1)*(dpZ-kkz1)
            U4=U41*U42
            U51=U(ix2,jy1,kz1)
            U52=(dpX-iix1)*(jjy2-dpY)*(kkz2-dpZ)
            U5=U51*U52
            U61=U(ix2,jy1,kz2)
            U62=(dpX-iix1)*(jjy2-dpY)*(dpZ-kkz1)
            U6=U61*U62
            U71=U(ix2,jy2,kz1)
            U72=(dpX-iix1)*(dpY-jjy1)*(kkz2-dpZ)
            U7=U71*U72
            U81=U(ix2,jy2,kz2)
            U82=(dpX-iix1)*(dpy-jjy1)*(dpz-kkz1)
            U8=U81*U82
            Ui=(U1+U2+U3+U4+U5+U6+U7+U8)  !rezalt U in point

            V11=V(ix1,jy1,kz1)
            V12=(iix2-dpX)*(jjy2-dpY)*(kkz2-dpz)
            V1=V11*V12
            V21=U(ix1,jy1,kz2)
            V22=(iix2-dpX)*(jjy2-dpY)*(dpz-kkz1)
            V2=V21*V22
            V31=V(ix1,jy2,kz1)
            V32=(iix2-dpX)*(dpY-jjy1)*(kkz2-dpZ)
            V3=V31*V32
            V41=V(ix1,jy2,kz2)
            V42=(iix2-dpX)*(dpY-jjy1)*(dpZ-kkz1)
            V4=V41*V42
            V51=V(ix2,jy1,kz1)
            V52=(dpX-iix1)*(jjy2-dpY)*(kkz2-dpZ)
            V5=V51*V52
            V61=V(ix2,jy1,kz2)
            V62=(dpX-iix1)*(jjy2-dpY)*(dpZ-kkz1)
            V6=V61*V62
            V71=V(ix2,jy2,kz1)
            V72=(dpX-iix1)*(dpY-jjy1)*(kkz2-dpZ)
            V7=V71*V72
            V81=V(ix2,jy2,kz2)
            V82=(dpX-iix1)*(dpy-jjy1)*(dpz-kkz1)
            V8=V81*V82
            Vi=(V1+V2+V3+V4+V5+V6+V7+V8)  !rezalt V in point
            if (Ui.ge.1700) then
                print*, begX, begY
                print*, U1, U3, U5, U7
                print*, U2, U4, U6, U8
                print*, pointX, pointY
                print*, Ui, dpX, dpY, dpZ
                pause 1
            endif
            if (Vi.ge.1700) then
                print*, begX, begY
                print*, V1, V3, V5, V7
                print*, V2, V4, V6, V8
                print*, pointX, pointY
                print*, Vi, dpX, dpY, kkz1, dpZ, kkz2
                pause 2
            endif

            !convert geographic coordinates to meters
            pointXm=(pointX-begX)*pi/180.*R*Cos(pointY*pi/180.)
            pointYm=(pointY-begY)*pi/180.*R

            !Random part
            c=1313*(iday+iitime)
            call seed(c)
            sigma1=sqrt(Ui*Ui+Vi*Vi)*dt*pi
            sigma2=5.*pi*dt
                call randu(s1,s2,rr)
                rd=2.*rr-1.
                pointXrand=sigma1*rd*a1
                call randu(s1,s2,rr)
                rd=2.*rr-1.
                pointYrand=sigma1*rd*a1
                call randu(s1,s2,rr)
                rd=2.*rr-1.
                pointZrand=sigma2*rd*a2

            !calculation U and V
            dpointXm=Ui*dt/1000.
            dpointYm=Vi*dt/1000.
            pointXm=pointXm+dpointXm+pointXrand
            pointYm=pointYm+dpointYm+pointYrand

            !calculation Z
            kzz=real(kz)
            iitime2=real(iitime)
            pointZ=125.-125.*sin((iitime2/nt1)*2.*pi+fi)+pointZrand

            !convert meters to geographic coordinates
            pointY=(pointYm/pi)*(180./R)+begY
            pointX=(pointXm/pi)*(180./R)/Cos(pointY*pi/180.)+begX

            if (pointY.ge.(-52.96)) pointY=-52.96
            if (pointX.ge.(-26.0)) pointX=-26.0
            if (pointY.le.(-62.0)) pointY=-62.0
            if (pointX.le.(-64.0)) pointX=-64.0
            if (pointZ.le.0.001) pointZ=0.001
            if (pointZ.ge.250.0) pointZ=249.9

            !output massive with new dat
            Cr(i,1)=pointX
            Cr(i,2)=pointY
            Cr(i,3)=pointZ
404         continue

        enddo
    enddo


    !write output data
    !open(1,file=Crout,status='new')
    !do j=1,n2,1
    !    write(1,*) (Cr(j,i),i=1,3,1)
    !enddo
    !close(1)

    do p=0,80,10
        pp=1
        if (iday.ge.1.and.iday.le.9) then
            do j=1,n2,1
                if (Cr(j,3).gt.p.and.Cr(j,3).le.(p+10)) then
                    pt=(p+10)/10
                    pt=int(pt)
                    write(Crout2,'(a4,i1,a1,i1,a4)') cr1,iday,'_',pt,frmt
                    Crz(pp,1)=Cr(j,1)
                    Crz(pp,2)=Cr(j,2)
                    Crz(pp,3)=Cr(j,3)
                    open(1,file=crout2)
                    write(1,*) (Crz(pp,i),i=1,3)
                    pp=pp+1
                endif
            enddo
        else
            if (iday.ge.10.and.iday.le.99) then
                do j=1,n2,1
                    if (Cr(j,3).gt.p.and.Cr(j,3).le.(p+10)) then
                        pt=(p+10)/10
                        pt=int(pt)
                        write(Crout2,'(a3,i2,a1,i1,a4)') cr2,iday,'_',pt,frmt
                        Crz(pp,1)=Cr(j,1)
                        Crz(pp,2)=Cr(j,2)
                        Crz(pp,3)=Cr(j,3)
                        open(1,file=crout2)
                        write(1,*) (Crz(pp,i),i=1,3)
                        pp=pp+1
                    endif
                enddo
            else
                if (iday.ge.100.and.iday.le.999) then
                    do j=1,n2,1
                        if (Cr(j,3).ge.p.and.Cr(j,3).le.(p+10)) then
                            pt=(p+10)/10
                            pt=int(pt)
                            write(Crout2,'(a2,i3,a1,i1,a4)') cr3,iday,'_',pt,frmt
                            Crz(pp,1)=Cr(j,1)
                            Crz(pp,2)=Cr(j,2)
                            Crz(pp,3)=Cr(j,3)
                            open(1,file=crout2)
                            write(1,*) (Crz(pp,i),i=1,3)
                            pp=pp+1
                        endif
                    enddo
                endif
            endif
        endif
    enddo

    do p=90,190,10
        pp=1
        if (iday.ge.1.and.iday.le.9) then
            do j=1,n2,1
                if (Cr(j,3).gt.p.and.Cr(j,3).le.(p+10)) then
                    pt=(p+10)/10
                    pt=int(pt)
                    write(Crout2,'(a4,i1,a1,i2,a4)') cr1,iday,'_',pt,frmt
                    Crz(pp,1)=Cr(j,1)
                    Crz(pp,2)=Cr(j,2)
                    Crz(pp,3)=Cr(j,3)
                    open(1,file=crout2)
                    write(1,*) (Crz(pp,i),i=1,3)
                    pp=pp+1
                endif
            enddo
        else
            if (iday.ge.10.and.iday.le.99) then
                do j=1,n2,1
                    if (Cr(j,3).gt.p.and.Cr(j,3).le.(p+10)) then
                        pt=(p+10)/10
                        pt=int(pt)
                        write(Crout2,'(a3,i2,a1,i2,a4)') cr2,iday,'_',pt,frmt
                        Crz(pp,1)=Cr(j,1)
                        Crz(pp,2)=Cr(j,2)
                        Crz(pp,3)=Cr(j,3)
                        open(1,file=crout2)
                        write(1,*) (Crz(pp,i),i=1,3)
                        pp=pp+1
                    endif
                enddo
            else
                if (iday.ge.100.and.iday.le.999) then
                    do j=1,n2,1
                        if (Cr(j,3).ge.p.and.Cr(j,3).le.(p+10)) then
                            pt=(p+10)/10
                            pt=int(pt)
                            write(Crout2,'(a2,i3,a1,i2,a4)') cr3,iday,'_',pt,frmt
                            Crz(pp,1)=Cr(j,1)
                            Crz(pp,2)=Cr(j,2)
                            Crz(pp,3)=Cr(j,3)
                            open(1,file=crout2)
                            write(1,*) (Crz(pp,i),i=1,3)
                            pp=pp+1
                        endif
                    enddo
                endif
            endif
        endif
    enddo
    print*, 'Day', iday, 'complete'
enddo

print*, 'The process is completed'
print*, 'Press enter to exit'
read(*,*)
pause
endprogram