Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- !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
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement