Skip to content

Commit

Permalink
Bug fix for mass conservation (replace pbavg by pbsavu,pbsavv)
Browse files Browse the repository at this point in the history
  • Loading branch information
@abozec
abozec committed May 11, 2016
1 parent 4f99570 commit adfdffe
Show file tree
Hide file tree
Showing 7 changed files with 99 additions and 65 deletions.
8 changes: 4 additions & 4 deletions barotp.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -263,7 +263,7 @@ subroutine barotp(m,n)
if (iuopn(i,j).ne.0) then
!!Alex river flxloc(i,j) = ubavg(i,j,ml)*(depthu(i,j)*scuy(i,j))
flxloc(i,j) = ubavg(i,j,ml)*scuy(i,j)*
& max(0.,pbot(i,j)*onetai(i,j)+pbavg(i,j,ml))
& max(0.,pbot(i,j)*onetai(i,j)+pbsavu(i,j))
uflxba(i,j) = uflxba(i,j) +
& coeflx(lll+1,icof)*(1.0+wblpf)*flxloc(i,j)
endif !iuopn
Expand All @@ -273,7 +273,7 @@ subroutine barotp(m,n)
if (iuopn(i,j).ne.0) then
!!Alex river flxloc(i,j) = ubavg(i,j,ml)*(depthu(i,j)*scuy(i,j))
flxloc(i,j) = ubavg(i,j,ml)*scuy(i,j)*
& max(0.,pbot(i-1,j)*onetai(i-1,j)+pbavg(i-1,j,ml))
& max(0.,pbot(i-1,j)*onetai(i-1,j)+pbsavu(i,j))
uflxba(i,j) = uflxba(i,j) +
& coeflx(lll+1,icof)*(1.0+wblpf)*flxloc(i,j)
endif !iuopn
Expand All @@ -295,7 +295,7 @@ subroutine barotp(m,n)
if (ivopn(i,j).ne.0) then
!!Alex river flyloc(i,j) = vbavg(i,j,ml)*(depthv(i,j)*scvx(i,j))
flyloc(i,j) = vbavg(i,j,ml)*scvx(i,j)*
& max(0.,pbot(i,j)*onetai(i,j)+pbavg(i,j,ml))
& max(0.,pbot(i,j)*onetai(i,j)+pbsavv(i,j))
vflxba(i,j) = vflxba(i,j) +
& coeflx(lll+1,icof)*(1.0+wblpf)*flyloc(i,j)
endif !ivopn
Expand All @@ -305,7 +305,7 @@ subroutine barotp(m,n)
if (ivopn(i,j).ne.0) then
!!Alex river flyloc(i,j) = vbavg(i,j,ml)*(depthv(i,j)*scvx(i,j))
flyloc(i,j) = vbavg(i,j,ml)*scvx(i,j)*
& max(0.,pbot(i,j-1)*onetai(i,j-1)+pbavg(i,j-1,ml))
& max(0.,pbot(i,j-1)*onetai(i,j-1)+pbsavv(i,j))
vflxba(i,j) = vflxba(i,j) +
& coeflx(lll+1,icof)*(1.0+wblpf)*flyloc(i,j)
endif !ivopn
Expand Down
11 changes: 7 additions & 4 deletions cnuity.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -106,6 +106,9 @@ subroutine cnuity(m,n)
cCL/RB MPI bug correction 2011-01 !!Alex
call xctilr(oneta( 1-nbdy,1-nbdy, n),1, 1, 6,6, halo_ps)
call xctilr(oneta( 1-nbdy,1-nbdy, m),1, 1, 6,6, halo_ps)
c --- update pbavg rivers
call xctilr(pbsavu( 1-nbdy,1-nbdy ),1, 1, 6,6, halo_ps)
call xctilr(pbsavv( 1-nbdy,1-nbdy ),1, 1, 6,6, halo_ps)
c
c --- rhs: dpmixl.n
c --- lhs: util3, utotn, vtotn
Expand Down Expand Up @@ -162,7 +165,7 @@ subroutine cnuity(m,n)
#endif
else
q=min(dp(i,j,k,n),max(0.,depthu(i+1,j)-util3(i,j)))
& * (onetai(i,j) + pbavg(i,j,m) /pbot(i,j))
& * (onetai(i,j) + pbsavu(i,j) /pbot(i,j))
utotm(i,j)=ubavg(i,j,m)*scuy(i,j) ! no stock drift for rivers
endif

Expand All @@ -188,7 +191,7 @@ subroutine cnuity(m,n)
else
q=min(dp(i-1,j,k,n),max(0.,depthu(i-1,j)
& - util3(i-1,j)))
& * (onetai(i-1,j) + pbavg(i-1,j,m) /pbot(i-1,j))
& * (onetai(i-1,j) + pbsavu(i,j) /pbot(i-1,j))
utotm(i,j)=ubavg(i,j,m)*scuy(i,j) ! no stock drift for rivers
endif
uflux(i,j)=utotm(i,j)*q
Expand Down Expand Up @@ -224,7 +227,7 @@ subroutine cnuity(m,n)
else
q=min(dp(i,j ,k,n),max(0.,depthv(i,j+1)
& - util3(i,j )))
& * (onetai(i,j) + pbavg(i,j,m) /pbot(i,j))
& * (onetai(i,j) + pbsavv(i,j) /pbot(i,j))
vtotm(i,j)=vbavg(i,j,m)*scvx(i,j) ! no stock drift for rivers
endif
vflux(i,j)=vtotm(i,j)*q
Expand All @@ -249,7 +252,7 @@ subroutine cnuity(m,n)
else
q=min(dp(i,j-1,k,n),max(0.,depthv(i,j-1)
& - util3(i,j-1)))
& * (onetai(i,j-1) + pbavg(i,j-1,m) /pbot(i,j-1))
& * (onetai(i,j-1) + pbsavv(i,j) /pbot(i,j-1))
vtotm(i,j)=vbavg(i,j,m)*scvx(i,j) ! no stock drift for rivers
endif
vflux(i,j)=vtotm(i,j)*q
Expand Down
18 changes: 16 additions & 2 deletions forfun.f
View file
Original file line number Diff line number Diff line change
Expand Up @@ -234,6 +234,7 @@ subroutine forfuna
cdiag. 'wind speed (x 10 m/s)')
endif !wndflg.lt.3
c
if (flxflg.ne.3) then
call zaiopf(flnmfor(1:lgth)//'forcing.airtmp.a', 'old', 904)
if (mnproc.eq.1) then ! .b file from 1st tile only
open (unit=uoff+904,file=flnmfor(1:lgth)//'forcing.airtmp.b',
Expand All @@ -255,6 +256,7 @@ subroutine forfuna
call rdmonth(util1, 905)
cdiag call prtmsk(ip,util1,util2,idm,idm,jdm, 0.,10000.,
cdiag. 'mixing ratio (0.1 g/kg)')
endif !flxflg.ne.3
c
if (pcipf) then
call zaiopf(flnmfor(1:lgth)//'forcing.precip.a', 'old', 906)
Expand Down Expand Up @@ -680,6 +682,7 @@ subroutine forfunh(dtime)
call preambl_print(preambl)
endif !wndflg.lt.3
c
if (flxflg.ne.3) then
call zaiopf(flnmfor(1:lgth)//'forcing.airtmp.a', 'old', 904)
if (mnproc.eq.1) then ! .b file from 1st tile only
open (unit=uoff+904,file=flnmfor(1:lgth)//'forcing.airtmp.b',
Expand All @@ -695,14 +698,17 @@ subroutine forfunh(dtime)
read (uoff+905,'(a79)') preambl
endif !1st tile
call preambl_print(preambl)
endif !flxflg.ne.3
c
if (pcipf) then
call zaiopf(flnmfor(1:lgth)//'forcing.precip.a', 'old', 906)
if (mnproc.eq.1) then ! .b file from 1st tile only
open (unit=uoff+906,file=flnmfor(1:lgth)//'forcing.precip.b',
& status='old', action='read')
read (uoff+906,'(a79)') preambl
endif !1st tile
call preambl_print(preambl)
endif !pcipf
c
call zaiopf(flnmfor(1:lgth)//'forcing.radflx.a', 'old', 907)
if (mnproc.eq.1) then ! .b file from 1st tile only
Expand Down Expand Up @@ -2608,8 +2614,13 @@ subroutine rdpall(dtime0,dtime1)
& taux(1-nbdy,1-nbdy,2),
& tauy(1-nbdy,1-nbdy,2) )
endif !wndspd
call rdpall1(airtmp,dtime(904),904,mod(icall,3).eq.1)
call rdpall1(vapmix,dtime(905),905,mod(icall,3).eq.1)
if (flxflg.ne.3) then
call rdpall1(airtmp,dtime(904),904,mod(icall,3).eq.1)
call rdpall1(vapmix,dtime(905),905,mod(icall,3).eq.1)
else
dtime(904) = dtime(900)
dtime(905) = dtime(900)
endif
if (mslprf .or. flxflg.eq.6) then
call rdpall1(mslprs,dtime(899),899,mod(icall,3).eq.2)
else
Expand Down Expand Up @@ -2992,8 +3003,10 @@ subroutine rdforf(mnth,lslot)
& taux(1-nbdy,1-nbdy,lslot),
& tauy(1-nbdy,1-nbdy,lslot) )
endif !wndspd
if (flxflg.ne.3) then
call rdmonthck(airtmp(1-nbdy,1-nbdy,lslot),904,mnth)
call rdmonthck(vapmix(1-nbdy,1-nbdy,lslot),905,mnth)
endif
if (pcipf) then
call rdmonthck(precip(1-nbdy,1-nbdy,lslot),906,mnth)
endif
Expand Down Expand Up @@ -4351,3 +4364,4 @@ subroutine str2spd(wspd, tx,ty)
c> Oct 2014 - flxflg==6 requires mslprs input
c> May 2015 - smooth msl pressure
c> May 2015 - forcw0 for constant in time between inputs
c> May 2016 - flxflg==3 does not read airtmp and vapmix
70 changes: 27 additions & 43 deletions latbdt_river.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -500,9 +500,9 @@ subroutine latbdt_river(n,lll)
else
if (i.ge.i0+1-nbdy.and.i.le.i0+ii+nbdy) then
do k = max(1-nbdy,j-j0), min(jj+nbdy,j-j0+lnport(l)-1)
pline(k+j0) = depths(i-i0,k)
rspedw(k+j0,ll1) = onetai(i-i0,k)
xline(k+j0) = scuy(i-i0,k) !!Alex
pline(k+j0) = depths(i-i0,k)
rspedw(k+j0,ll1) = onetai(i-i0,k)
xline(k+j0) = scuy(i-i0,k) !!Alex
enddo
endif
endif
Expand Down Expand Up @@ -545,9 +545,9 @@ subroutine latbdt_river(n,lll)
else
if (i.ge.i0+1-nbdy.and.i.le.i0+ii+nbdy) then
do k = max(1-nbdy,j-j0), min(jj+nbdy,j-j0+lnport(l)-1)
pline(k+j0) = depths(i-i0,k)
pline(k+j0) = depths(i-i0,k)
rspede(k+j0,ll1) = onetai(i-i0,k)
xline(k+j0) = scuy(i+1-i0,k) !!Alex
xline(k+j0) = scuy(i+1-i0,k) !!Alex
!RB xline(k+j0) = scuy(i-i0,k)
enddo
endif
Expand Down Expand Up @@ -591,9 +591,9 @@ subroutine latbdt_river(n,lll)
else
if (j.ge.j0+1-nbdy.and.j.le.j0+jj+nbdy) then
do k = max(1-nbdy,i-i0), min(ii+nbdy,i-i0+lnport(l)-1)
pline(k+i0) = depths(k,j-j0)
pline(k+i0) = depths(k,j-j0)
rspedn(k+i0,ll1) = onetai(k,j-j0)
xline(k+i0) = scvx(k,j+1-j0) !!Alex
xline(k+i0) = scvx(k,j+1-j0) !!Alex
!RB xline(k+i0) = scvx(k,j-j0)
enddo
endif
Expand Down Expand Up @@ -636,9 +636,9 @@ subroutine latbdt_river(n,lll)
else
if (j.ge.j0+1-nbdy.and.j.le.j0+jj+nbdy) then
do k = max(1-nbdy,i-i0), min(ii+nbdy,i-i0+lnport(l)-1)
pline(k+i0) = depths(k,j-j0)
rspeds(k+i0,ll1) = onetai(k,j-j0)
xline(k+i0) = scvx(k,j-j0) !!Alex
pline(k+i0) = depths(k,j-j0)
rspeds(k+i0,ll1) = onetai(k,j-j0)
xline(k+i0) = scvx(k,j-j0) !!Alex
enddo
endif
endif
Expand Down Expand Up @@ -973,14 +973,14 @@ subroutine latbdt_river(n,lll)
j = jfport(l)
if (boundary_mpi) then
call xclput(pline(j), lnport(l),
& pbavg(1-nbdy,1-nbdy,n), i, j, 0,1)
& pbsavu(1-nbdy,1-nbdy ), i, j, 0,1)
call xclput(uline(j,1),lnport(l),
& ubavg(1-nbdy,1-nbdy,n), i, j, 0,1) ! normal
& ubavg (1-nbdy,1-nbdy,n), i, j, 0,1) ! normal
else
if (i.ge.i0+1-margin.and.i.le.i0+ii+margin) then
do k=max(1-margin,j-j0),min(jj+margin,j-j0+lnport(l)-1)
pbavg(i-i0,k,n) = pline(k+j0)
ubavg(i-i0,k,n) = uline(k+j0,1)
pbsavu(i-i0,k ) = pline(k+j0)
ubavg (i-i0,k,n) = uline(k+j0,1)
enddo
endif
endif ! boundary_mpi
Expand Down Expand Up @@ -1073,24 +1073,16 @@ subroutine latbdt_river(n,lll)
endif ! rarwave
j = jfport(l)
if (boundary_mpi) then
!RB call xclput(pline(j), lnport(l),
!RB & pbavg(1-nbdy,1-nbdy,n), i+1,j,0,1)
!RB call xclput(uline(j,1),lnport(l),
!RB & ubavg(1-nbdy,1-nbdy,n), i+1,j,0,1) ! normal

call xclput(pline(j), lnport(l),
& pbavg(1-nbdy,1-nbdy,n), i ,j,0,1)
call xclput(pline(j), lnport(l),
& pbsavu(1-nbdy,1-nbdy ), i+1,j,0,1)
call xclput(uline(j,1),lnport(l),
& ubavg(1-nbdy,1-nbdy,n), i+1,j,0,1) ! normal
& ubavg (1-nbdy,1-nbdy,n), i+1,j,0,1) ! normal
else
if (i.ge.i0-margin.and.i+1.le.i0+ii+margin) then
do k=max(1 -margin,j-j0),
& min(jj+margin,j-j0+lnport(l)-1)
!RB pbavg(i+1-i0,k,n) = pline(k+j0)
!RB ubavg(i+1-i0,k,n) = uline(k+j0,1)

pbavg(i -i0,k,n) = pline(k+j0)
ubavg(i+1-i0,k,n) = uline(k+j0,1)
pbsavu(i+1-i0,k ) = pline(k+j0)
ubavg (i+1-i0,k,n) = uline(k+j0,1)
enddo
endif
endif ! boundary_mpi
Expand Down Expand Up @@ -1181,24 +1173,16 @@ subroutine latbdt_river(n,lll)
endif ! rarwave
i = ifport(l)
if (boundary_mpi) then
!RB call xclput(pline(i), lnport(l),
!RB & pbavg(1-nbdy,1-nbdy,n), i, j+1,1,0)
!RB call xclput(uline(i,1),lnport(l),
!RB & vbavg(1-nbdy,1-nbdy,n), i, j+1,1,0) ! normal

call xclput(pline(i), lnport(l),
& pbavg(1-nbdy,1-nbdy,n), i, j ,1,0)
& pbsavv(1-nbdy,1-nbdy ), i, j+1,1,0)
call xclput(uline(i,1),lnport(l),
& vbavg(1-nbdy,1-nbdy,n), i, j+1,1,0) ! normal
& vbavg (1-nbdy,1-nbdy,n), i, j+1,1,0) ! normal
else
if (j.ge.j0-margin.and.j+1.le.j0+jj+margin) then
do k=max(1 -margin,i-i0),
& min(ii+margin,i-i0+lnport(l)-1)
!RB pbavg(k,j+1-j0,n) = pline(k+i0)
!RB vbavg(k,j+1-j0,n) = uline(k+i0,1)

pbavg(k,j -j0,n) = pline(k+i0)
vbavg(k,j+1-j0,n) = uline(k+i0,1)
pbsavv(k,j+1-j0 ) = pline(k+i0)
vbavg (k,j+1-j0,n) = uline(k+i0,1)
enddo
endif
endif ! boundary_mpi
Expand Down Expand Up @@ -1296,15 +1280,15 @@ subroutine latbdt_river(n,lll)
i = ifport(l)
if (boundary_mpi) then
call xclput(pline(i), lnport(l),
& pbavg(1-nbdy,1-nbdy,n), i, j, 1, 0)
& pbsavv(1-nbdy,1-nbdy), i, j, 1, 0)
call xclput(uline(i,1),lnport(l),
& vbavg(1-nbdy,1-nbdy,n), i, j, 1, 0) ! normal
& vbavg (1-nbdy,1-nbdy,n), i, j, 1, 0) ! normal
else
if (j.ge.j0+1-margin.and.j.le.j0+jj+margin) then
do k=max(1 -margin,i-i0),
& min(ii+margin,i-i0+lnport(l)-1)
pbavg(k,j-j0,n) = pline(k+i0)
vbavg(k,j-j0,n) = uline(k+i0,1)
pbsavv(k,j-j0 ) = pline(k+i0)
vbavg (k,j-j0,n) = uline(k+i0,1)
enddo
endif
endif ! boundary_mpi
Expand Down
17 changes: 17 additions & 0 deletions mod_cb_arrays.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -862,6 +862,14 @@ module mod_cb_arrays
real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,kdm) ::
#endif
& triver, sriver ! T and S river
c --- pbavg at river points
#if defined(RELO)
real,save,allocatable, dimension(:,:) ::
#else
real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) ::
#endif
& pbsavu, pbsavv
contains
Expand Down Expand Up @@ -1616,6 +1624,15 @@ subroutine cb_allocate
#endif
triver = 0.0
sriver = 0.0
c --- pbavg at river points
#if defined(RELO)
allocate(
& pbsavu(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy),
& pbsavv(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy))
call mem_stat_add( 2*(idm+2*nbdy)*(jdm+2*nbdy) )
#endif
pbsavu = 0.0
pbsavv = 0.0
end subroutine cb_allocate
Expand Down
37 changes: 27 additions & 10 deletions mod_hycom.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2798,7 +2798,11 @@ subroutine HYCOM_Run
util2(i,j)=(srfhgt(i,j)/g)**2*scp2(i,j)
endif
util3(i,j)=srfhgt(i,j)*scp2(i,j)
util4(i,j)=montg1(i,j)*scp2(i,j)
if (sshflg.ne.0) then
util4(i,j)=steric(i,j)*scp2(i,j)
else
util4(i,j)=montg1(i,j)*scp2(i,j)
endif !sshflg
sminy(j)=min(sminy(j),srfhgt(i,j))
smaxy(j)=max(smaxy(j),srfhgt(i,j))
endif !ip
Expand All @@ -2825,20 +2829,32 @@ subroutine HYCOM_Run
& '' ('',1pe8.1,'' to '',e8.1,'')'')')
& nstep,c_ydh,
& sum/(area*thref*onemm),smin/(thref*onemm),smax/(thref*onemm)
write (lp,'(i9,a,
. '' mean MontgPot (mm):'',f8.2)')
. nstep,c_ydh,
. smt/(area*thref*onemm)
call flush(lp)
write(nod,'(i9,a,
& '' mean SSH (mm):'',f8.2,
& '' ('',1pe8.1,'' to '',e8.1,'')'')')
& nstep,c_ydh,
& sum/(area*thref*onemm),smin/(thref*onemm),smax/(thref*onemm)
write(nod,'(i9,a,
. '' mean MontgPot (mm):'',f8.2)')
. nstep,c_ydh,
. smt/(area*thref*onemm)
if (sshflg.ne.0) then
c --- Non-Steric (i.e. Total - Steric) SSH
write (lp,'(i9,a,
& '' mean nsSSH (mm):'',f8.2)')
& nstep,c_ydh,
& (sum-smt)/(area*thref*onemm)
write(nod,'(i9,a,
& '' mean nsSSH (mm):'',f8.2)')
& nstep,c_ydh,
& (sum-smt)/(area*thref*onemm)
else
write (lp,'(i9,a,
& '' mean MontgPot (mm):'',f8.2)')
& nstep,c_ydh,
& smt/(area*thref*onemm)
write(nod,'(i9,a,
& '' mean MontgPot (mm):'',f8.2)')
& nstep,c_ydh,
& smt/(area*thref*onemm)
endif !sshflg
call flush(lp)
call flush(nod)
endif !1st tile
c --- NaN detection.
Expand Down Expand Up @@ -4055,3 +4071,4 @@ end module mod_hycom
c> May 2014 - use land/sea masks (e.g. ip) to skip land
c> Feb. 2015 - removed dtime0 from calculation of dtime
c> Feb. 2015 - fixed issues with bit for bit reproducability on restart
c> May 2016 - Printout mean non-steric SSH (not MontgPot) when available
Loading

0 comments on commit adfdffe

Please sign in to comment.