Commit 7725d964 authored by Sabine's avatar Sabine
Browse files

code included to interpolate hmix

parent 64810109
......@@ -116,6 +116,7 @@ subroutine advance(itime,nrelpoint,ldt,up,vp,wp, &
real :: dcwsave
real :: usigold,vsigold,wsigold,r,rs
real :: uold,vold,wold,vdepo(maxspec)
real :: h1(2)
!real uprof(nzmax),vprof(nzmax),wprof(nzmax)
!real usigprof(nzmax),vsigprof(nzmax),wsigprof(nzmax)
!real rhoprof(nzmax),rhogradprof(nzmax)
......@@ -222,6 +223,25 @@ subroutine advance(itime,nrelpoint,ldt,up,vp,wp, &
jyp=jy+1
! Determine the lower left corner and its distance to the current position
!*************************************************************************
ddx=xt-real(ix)
ddy=yt-real(jy)
rddx=1.-ddx
rddy=1.-ddy
p1=rddx*rddy
p2=ddx*rddy
p3=rddx*ddy
p4=ddx*ddy
! Calculate variables for time interpolation
!*******************************************
dt1=real(itime-memtime(1))
dt2=real(memtime(2)-itime)
dtt=1./(dt1+dt2)
! Compute maximum mixing height around particle position
!*******************************************************
......@@ -229,11 +249,16 @@ subroutine advance(itime,nrelpoint,ldt,up,vp,wp, &
if (ngrid.le.0) then
do k=1,2
mind=memind(k) ! eso: compatibility with 3-field version
do j=jy,jyp
do i=ix,ixp
if (hmix(i,j,1,mind).gt.h) h=hmix(i,j,1,mind)
end do
end do
! do j=jy,jyp
! do i=ix,ixp
! if (hmix(i,j,1,mind).gt.h) h=hmix(i,j,1,mind)
h1(k)=p1*hmix(ix ,jy ,1,mind) &
+ p2*hmix(ixp,jy ,1,mind) &
+ p3*hmix(ix ,jyp,1,mind) &
+ p4*hmix(ixp,jyp,1,mind)
! end do
! end do
end do
tropop=tropopause(nix,njy,1,1)
else
......@@ -248,6 +273,7 @@ subroutine advance(itime,nrelpoint,ldt,up,vp,wp, &
tropop=tropopausen(nix,njy,1,1,ngrid)
endif
h=(h1(1)*dt2+h1(2)*dt1)*dtt
zeta=zt/h
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment