forked from kthyng/tracpy
-
Notifications
You must be signed in to change notification settings - Fork 1
/
pos.f95
executable file
·157 lines (141 loc) · 6.25 KB
/
pos.f95
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
subroutine pos_orgn(ijk,ia,ja,ka,r0,r1,ds,rr,uflux,vflux,wflux,ff,imt,jmt,km,do3d,doturb,upr)
!====================================================================
! computes new position (r0 --> r1) of trajectory after time ds
! the new coordinate is still on one of the faces of box at ia,ja,ka
!
! Input:
!
! ijk : considered direction (1=zonal,2=meridional,3=vertical)
! ia,ja,ka : original position in integers
! r0 : original non-dimensional position in the ijk-direction
! of particle (fractions of a grid box side in the
! corresponding direction)
! ds : crossing time to reach the grid box wall (units=s/m3)
! rr : time interpolation constant between 0 and 1. Controls how much
! : of earlier time step is used in interpolation.
! uflux : u velocity (zonal) flux field, two time steps [ixjxkxt]
! vflux : v velocity (meridional) flux field, two time steps [ixjxkxt]
! wflux : w velocity (vertical) flux field, two time steps [kxt]
! ff : time direction. ff=1 forward, ff=-1 backward
! imt,jmt,km : grid index sizing constants in (x,y,z), are for
! horizontal and vertical rho grid [scalar]
! do3d : Flag to set whether to use 3d velocities or not
! doturb : Flag to set whether or not to use turb/diff and which kind if so
! upr : parameterized turbulent velocities u', v', and w'
! optional because only used if using turb flag for diffusion
! size [6,2]. The 2nd dimension is for two time steps.
! The 1st dimension is: [u'_ia,u'_ia-1,v'_ja,v'_ja-1,w'_ka,w'_ka-1]
!
! Output:
!
! r1 : the new position (coordinate)
!
! Other parameters used in function:
! rg : rg=1-rr for time interpolation between time steps. Controls how much
! : of later time step is used in interpolation.
! uu : time-interpolated flux at ia/ja/ka (depending on ijk)
! um : time-interpolated flux at ia-1/ja-1/ka-1 (depending on ijk)
! ii : generic index for grid index for whichever direction, ijk
! im : generic index for grid index -1 for whichever direction, ijk.
! Is only used in the i direction for whatever reason.
! nsm=1,nsp=2 : Time index. nsm picks out the earlier bounding time step and
! nsp picks out the later bounding time step for interpolation.
!====================================================================
implicit none
integer, intent(in) :: ijk,ia,ja,ka,ff,imt,jmt,km
integer, intent(in) :: do3d, doturb
real(kind=8), intent(in) :: r0,rr,ds
real(kind=8), intent(in), dimension(imt-1,jmt,km,2) :: uflux
real(kind=8), intent(in), dimension(imt,jmt-1,km,2) :: vflux
real(kind=8), intent(in), dimension(0:km,2) :: wflux
real*8, optional, intent(in), dimension(6,2) :: upr
real(kind=8), intent(out) :: r1
real*8 :: rg,uu,um
integer :: ii,im,nsm=1,nsp=2
rg=1.d0-rr
if(do3d==0) then
if(ijk.eq.3) then
r1=r0
return
endif
endif
if(ijk.eq.1) then
ii=ia
im=ia-1
! KMT: This appears to be a periodic boundary condition. I think this condition
! would arise only when the drifter is at the very west x edge, and I am not sure
! what the boundary condition here should be. Why isn't there one in the v/y direction?
! if(im.eq.0) im=imt
uu=(rg*uflux(ia,ja,ka,nsp)+rr*uflux(ia,ja,ka,nsm))*ff
um=(rg*uflux(im,ja,ka,nsp)+rr*uflux(im,ja,ka,nsm))*ff
if(doturb==1) then
if(r0.ne.dble(ii)) then
uu=uu+upr(1,2)
else
uu=uu+upr(1,1)
! add u' from previous iterative time step if on box wall
endif
if(r0.ne.dble(im)) then
um=um+upr(2,2)
else
um=um+upr(2,1)
! add u' from previous iterative time step if on box wall
endif
endif
else if(ijk.eq.2) then
ii=ja
uu=(rg*vflux(ia,ja ,ka,nsp)+rr*vflux(ia,ja ,ka,nsm))*ff
um=(rg*vflux(ia,ja-1,ka,nsp)+rr*vflux(ia,ja-1,ka,nsm))*ff
if(doturb==1) then
if(r0.ne.dble(ja )) then
uu=uu+upr(3,2)
else
uu=uu+upr(3,1)
! add u' from previous iterative time step if on box wall
endif
if(r0.ne.dble(ja-1)) then
um=um+upr(4,2)
else
um=um+upr(4,1)
! add u' from previous iterative time step if on box wall
endif
endif
else if(ijk.eq.3) then
ii=ka
! #ifdef full_wflux
! uu=wflux(ia ,ja ,ka ,nsm)
! um=wflux(ia ,ja ,ka-1 ,nsm)
! #else
uu=rg*wflux(ka ,nsp)+rr*wflux(ka ,nsm)
um=rg*wflux(ka-1,nsp)+rr*wflux(ka-1,nsm)
! #endif
if(do3d==1 .and. doturb==1) then
if(r0.ne.dble(ka )) then
uu=uu+upr(5,2)
else
uu=uu+upr(5,1)
! add u' from previous iterative time step if on box wall
endif
if(r0.ne.dble(ka-1)) then
uu=uu+upr(6,2)
else
uu=uu+upr(6,1)
! add u' from previous iterative time step if on box wall
endif
endif
endif
!
! note: consider in future to improve the code below for accuracy
! in case of um-uu = small; also see subroutine cross
if(um.ne.uu) then
r1= (r0+(-dble(ii-1) + um/(uu-um))) * dexp( (uu-um)*ds ) + dble(ii-1) - um/(uu-um)
else
r1=r0+uu*ds
endif
! print '(a,i1,a,f6.2,a,f6.2,a,e6.1,a,f4.2)','ijk=',ijk,' r0=',r0,' r1=',r1,' ds=',ds,' rr=',rr
! print '(a,f8.1,a,f8.1)','uu=',uu,' um=',um
! print '(a,f8.1,a,f8.1,a)','vflux(ia,ja,ka,1)=',vflux(ia,ja,ka,1),' vflux(ia ,ja,ka,2)=',vflux(ia ,ja,ka,2)
! print '(a,f7.1,a,f6.1,a,f5.2,a,f5.2)', 'tt=',tt,' dt=',dt,' ts=',ts,' tss=',tss
!if(abs(um/(uu-um)).gt.1.d10) print *,'possible precision problem?',um/(uu-um),uu,um,ijk,ia,ja,ka,r0,r1,ds,rr
return
end subroutine pos_orgn