Created
March 28, 2017 23:51
-
-
Save m1t9/69628ef68f6c39685aa09d766f200b05 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
!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 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment