Skip to content

Commit

Permalink
Fix the code to run timeanalytical stepping
Browse files Browse the repository at this point in the history
  • Loading branch information
@brorfred
brorfred committed Nov 24, 2015
1 parent 6592ff2 commit 91d8e45
Show file tree
Hide file tree
Showing 3 changed files with 151 additions and 61 deletions.
154 changes: 114 additions & 40 deletions src/interp.f95
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,56 +65,130 @@ subroutine interp(ib,jb,kb,x1,y1,z1,temp,salt,dens,ns)
az=(dble(kp)-z1)

! temperature, salinity, density calculation
tppp=tem(ip,jp,kp,ns)
if(tppp.eq.0.) tppp=tem(ib,jb,kb,ns)
tppm=tem(ip,jp,kn,ns)
if(tppm.eq.0.) tppm=tem(ib,jb,kb,ns)
tpmp=tem(ip,jm,kp,ns)
if(tpmp.eq.0.) tpmp=tem(ib,jb,kb,ns)
tpmm=tem(ip,jm,kn,ns)
if(tpmm.eq.0.) tpmm=tem(ib,jb,kb,ns)
tmpp=tem(im,jp,kp,ns)
if(tmpp.eq.0.) tmpp=tem(ib,jb,kb,ns)
tmpm=tem(im,jp,kn,ns)
if(tmpm.eq.0.) tmpm=tem(ib,jb,kb,ns)
tmmp=tem(im,jm,kp,ns)
if(tmmp.eq.0.) tmmp=tem(ib,jb,kb,ns)
tmmm=tem(im,jm,kn,ns)
if(tmmm.eq.0.) tmmm=tem(ib,jb,kb,ns)
! tppp=tem(ip,jp,kp,ns)
! if(tppp.eq.0.) tppp=tem(ib,jb,kb,ns)
! tppm=tem(ip,jp,kn,ns)
! if(tppm.eq.0.) tppm=tem(ib,jb,kb,ns)
! tpmp=tem(ip,jm,kp,ns)
! if(tpmp.eq.0.) tpmp=tem(ib,jb,kb,ns)
! tpmm=tem(ip,jm,kn,ns)
! if(tpmm.eq.0.) tpmm=tem(ib,jb,kb,ns)
! tmpp=tem(im,jp,kp,ns)
! if(tmpp.eq.0.) tmpp=tem(ib,jb,kb,ns)
! tmpm=tem(im,jp,kn,ns)
! if(tmpm.eq.0.) tmpm=tem(ib,jb,kb,ns)
! tmmp=tem(im,jm,kp,ns)
! if(tmmp.eq.0.) tmmp=tem(ib,jb,kb,ns)
! tmmm=tem(im,jm,kn,ns)
! if(tmmm.eq.0.) tmmm=tem(ib,jb,kb,ns)
!
! sppp=sal(ip,jp,kp,ns)
! if(sppp.eq.0.) sppp=sal(ib,jb,kb,ns)
! sppm=sal(ip,jp,kn,ns)
! if(sppm.eq.0.) sppm=sal(ib,jb,kb,ns)
! spmp=sal(ip,jm,kp,ns)
! if(spmp.eq.0.) spmp=sal(ib,jb,kb,ns)
! spmm=sal(ip,jm,kn,ns)
! if(spmm.eq.0.) spmm=sal(ib,jb,kb,ns)
! smpp=sal(im,jp,kp,ns)
! if(smpp.eq.0.) smpp=sal(ib,jb,kb,ns)
! smpm=sal(im,jp,kn,ns)
! if(smpm.eq.0.) smpm=sal(ib,jb,kb,ns)
! smmp=sal(im,jm,kp,ns)
! if(smmp.eq.0.) smmp=sal(ib,jb,kb,ns)
! smmm=sal(im,jm,kn,ns)
! if(smmm.eq.0.) smmm=sal(ib,jb,kb,ns)
!
! rppp=rho(ip,jp,kp,ns)
! if(rppp.eq.0.) rppp=rho(ib,jb,kb,ns)
! rppm=rho(ip,jp,kn,ns)
! if(rppm.eq.0.) rppm=rho(ib,jb,kb,ns)
! rpmp=rho(ip,jm,kp,ns)
! if(rpmp.eq.0.) rpmp=rho(ib,jb,kb,ns)
! rpmm=rho(ip,jm,kn,ns)
! if(rpmm.eq.0.) rpmm=rho(ib,jb,kb,ns)
! rmpp=rho(im,jp,kp,ns)
! if(rmpp.eq.0.) rmpp=rho(ib,jb,kb,ns)
! rmpm=rho(im,jp,kn,ns)
! if(rmpm.eq.0.) rmpm=rho(ib,jb,kb,ns)
! rmmp=rho(im,jm,kp,ns)
! if(rmmp.eq.0.) rmmp=rho(ib,jb,kb,ns)
! rmmm=rho(im,jm,kn,ns)
! if(rmmm.eq.0.) rmmm=rho(ib,jb,kb,ns)

tppp=tem(ip,jp,kp,ns)
sppp=sal(ip,jp,kp,ns)
if(sppp.eq.0.) sppp=sal(ib,jb,kb,ns)
sppm=sal(ip,jp,kn,ns)
if(sppm.eq.0.) sppm=sal(ib,jb,kb,ns)
spmp=sal(ip,jm,kp,ns)
if(spmp.eq.0.) spmp=sal(ib,jb,kb,ns)
spmm=sal(ip,jm,kn,ns)
if(spmm.eq.0.) spmm=sal(ib,jb,kb,ns)
smpp=sal(im,jp,kp,ns)
if(smpp.eq.0.) smpp=sal(ib,jb,kb,ns)
smpm=sal(im,jp,kn,ns)
if(smpm.eq.0.) smpm=sal(ib,jb,kb,ns)
smmp=sal(im,jm,kp,ns)
if(smmp.eq.0.) smmp=sal(ib,jb,kb,ns)
smmm=sal(im,jm,kn,ns)
if(smmm.eq.0.) smmm=sal(ib,jb,kb,ns)

rppp=rho(ip,jp,kp,ns)
if(rppp.eq.0.) rppp=rho(ib,jb,kb,ns)
if(tppp==0. .and. sppp==0.) then
tppp=tem(ip,jp,kn,ns)
sppp=sal(ip,jp,kn,ns)
rppp=rho(ip,jp,kn,ns)
endif

tppm=tem(ip,jp,kn,ns)
sppm=sal(ip,jp,kn,ns)
rppm=rho(ip,jp,kn,ns)
if(rppm.eq.0.) rppm=rho(ib,jb,kb,ns)
if(tppm==0. .and. sppm==0.) then
tppm=tem(ip,jp,kn,ns)
sppm=sal(ip,jp,kn,ns)
rppm=rho(ip,jp,kn,ns)
endif

tpmp=tem(ip,jm,kp,ns)
spmp=sal(ip,jm,kp,ns)
rpmp=rho(ip,jm,kp,ns)
if(rpmp.eq.0.) rpmp=rho(ib,jb,kb,ns)
if(tpmp==0. .and. spmp==0.) then
tpmp=tem(ip,jp,kn,ns)
spmp=sal(ip,jp,kn,ns)
rpmp=rho(ip,jp,kn,ns)
endif

tpmm=tem(ip,jm,kn,ns)
spmm=sal(ip,jm,kn,ns)
rpmm=rho(ip,jm,kn,ns)
if(rpmm.eq.0.) rpmm=rho(ib,jb,kb,ns)
if(tpmm==0. .and. spmm==0.) then
tpmm=tem(ip,jp,kn,ns)
spmm=sal(ip,jp,kn,ns)
rpmm=rho(ip,jp,kn,ns)
endif

tmpp=tem(im,jp,kp,ns)
smpp=sal(im,jp,kp,ns)
rmpp=rho(im,jp,kp,ns)
if(rmpp.eq.0.) rmpp=rho(ib,jb,kb,ns)
if(tmpp==0. .and. smpp==0.) then
tmpp=tem(ip,jp,kn,ns)
smpp=sal(ip,jp,kn,ns)
rmpp=rho(ip,jp,kn,ns)
endif

tmpm=tem(im,jp,kn,ns)
smpm=sal(im,jp,kn,ns)
rmpm=rho(im,jp,kn,ns)
if(rmpm.eq.0.) rmpm=rho(ib,jb,kb,ns)
if(tmpm==0. .and. smpm==0.) then
tmpm=tem(ip,jp,kn,ns)
smpm=sal(ip,jp,kn,ns)
rmpm=rho(ip,jp,kn,ns)
endif

tmmp=tem(im,jm,kp,ns)
smmp=sal(im,jm,kp,ns)
rmmp=rho(im,jm,kp,ns)
if(rmmp.eq.0.) rmmp=rho(ib,jb,kb,ns)
if(tmmp==0. .and. smmp==0.) then
tmmp=tem(ip,jp,kn,ns)
smmp=sal(ip,jp,kn,ns)
rmmp=rho(ip,jp,kn,ns)
endif

tmmm=tem(im,jm,kn,ns)
smmm=sal(im,jm,kn,ns)
rmmm=rho(im,jm,kn,ns)
if(rmmm.eq.0.) rmmm=rho(ib,jb,kb,ns)
if(tmmm==0. .and. smmm ==0.) then
tmmm=tem(ip,jp,kn,ns)
smmm=sal(ip,jp,kn,ns)
rmmm=rho(ip,jp,kn,ns)
endif



temp=tppp*(1.-ax)*(1.-ay)*(1.-az) &
+ tmpp* ax *(1.-ay)*(1.-az) &
Expand Down
21 changes: 20 additions & 1 deletion src/loop_pos.f95
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,9 @@ module mod_pos
USE mod_time, only: intrpb, intrpbg
USE mod_streamfunctions
USE mod_psi

USE mod_tempsalt
USE mod_traj, only: ntrac

IMPLICIT none

contains
Expand Down Expand Up @@ -328,6 +330,23 @@ subroutine pos(ia,iam,ja,ka,ib,jb,kb,x0,y0,z0,x1,y1,z1)
call pos_time(1,ia,ja,ka,x0,x1)
call pos_time(2,ia,ja,ka,y0,y1)
call pos_time(3,ia,ja,ka,z0,z1)
if(dble(int(x0))/=x0 .and. abs(int(x0)-int(x1))>=1.) then
print *,'pangx=',x0,x1,y0,y1,z0,z1
! print *,'ds',dse,dsw,dsn,dss,dsu,dsd,dsc,dsmin
x1=x0
! stop 3854
endif
if(dble(int(y0))/=y0 .and. abs(int(y0)-int(y1))>=1.) then
print *,'pangy=',y0,y1
! stop 3854
y1=y0
endif
if(dble(int(z0))/=z0 .and. abs(int(z0)-int(z1))>=1.) then
print *,'pangz=',ntrac,z0,z1
z1=z0
stop 3854
endif
! print *,'x0x1',x0,x1,y0,y1,z0,z1
#else
! If there is no spatial solution, i.e. a convergence zone
if(dse==UNDEF .and. dsw==UNDEF .and. dsn==UNDEF .and. &
Expand Down
37 changes: 17 additions & 20 deletions src/vertvel.f95
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,25 +24,29 @@ subroutine vertvel(ia,iam,ja,ka)
#if defined twodim || explicit_w
return
#else

n1=min(nsm,nsp)
n2=max(nsm,nsp)

kloop: do k=1,ka
uu = intrpg * uflux(ia ,ja ,k,nsp) + intrpr * uflux(ia ,ja ,k,nsm)
um = intrpg * uflux(iam,ja ,k,nsp) + intrpr * uflux(iam,ja ,k,nsm)
vv = intrpg * vflux(ia ,ja ,k,nsp) + intrpr * vflux(ia ,ja ,k,nsm)
vm = intrpg * vflux(ia ,ja-1,k,nsp) + intrpr * vflux(ia ,ja-1,k,nsm)
#if defined zgrid3Dt
do n=n1,n2
! time change of the mass the in grid box
wflux(k,n) = wflux(k-1,n) - ff * &
( uflux(ia,ja,k,n) - uflux(iam, ja, k, n) + &
vflux(ia,ja,k,n) - vflux(ia, ja-1, k, n) + &
(dzt(ia,ja,k,nsp)-dzt(ia,ja,k,nsm))*dxdy(ia,ja)/tseas )
( uflux(ia,ja,k,n) - uflux(iam, ja, k, n) &
+ vflux(ia,ja,k,n) - vflux(ia, ja-1, k, n) &
#if defined ifs
+ (dzt(ia,ja,k,n2)-dzt(ia,ja,k,n1))*dxdy(ia,ja)/tseas )
#else
- (dzt(ia,ja,k,n2)-dzt(ia,ja,k,n1))*dxdy(ia,ja)/tseas )
#endif
! ska det vara plus för ifs och minus för orca???
enddo
#else
#if defined full_wflux
uu = intrpg * uflux(ia ,ja ,k,nsp) + intrpr * uflux(ia ,ja ,k,nsm)
um = intrpg * uflux(iam,ja ,k,nsp) + intrpr * uflux(iam,ja ,k,nsm)
vv = intrpg * vflux(ia ,ja ,k,nsp) + intrpr * vflux(ia ,ja ,k,nsm)
vm = intrpg * vflux(ia ,ja-1,k,nsp) + intrpr * vflux(ia ,ja-1,k,nsm)
wflux(ia,ja,k,nsm)=wflux(ia,ja,k-1,nsm) - ff * ( uu - um + vv - vm )
#else
do n=n1,n2
Expand All @@ -55,19 +59,12 @@ subroutine vertvel(ia,iam,ja,ka)
end do kloop


#ifdef ifs
! Make sure the vertical velocity is always zero at the bottom and below
#if defined ifs
wflux(0,:) = 0.d0
wflux(km,:) = 0.d0
#endif

! Make sure the vertical velocity is always zero below and at the bottom
#ifdef orca12
#elif defined orca1 || orca12
do k=0,KM-kmt(ia,ja)
do n=n1,n2
if(wflux(k,n)/=0.d0) then
wflux(k,n)=0.d0
endif
enddo
wflux(k,:)=0.d0
enddo
#endif

Expand Down

0 comments on commit 91d8e45

Please sign in to comment.