Skip to content

GitLab

Commit ad0a46eb authored by Sabine's avatar Sabine
Browse files

Merge branch 'watercycle' into testmerge_h2o

parents 2bfec12f fd6eec30
Hide whitespace changes
Inline Side-by-side
src/boundcond_domainfill.f90
View file @ ad0a46eb
......@@ -79,14 +79,16 @@ subroutine boundcond_domainfill(itime,loutend)
!********************************************************************
do i=1,numpart
if (itra1(i).eq.itime) then
if ((ytra1(i).gt.real(ny_sn(2))).or. &
(ytra1(i).lt.real(ny_sn(1)))) itra1(i)=-999999999
if (((.not.xglobal).or.(nx_we(2).ne.(nx-2))).and. &
((xtra1(i).lt.real(nx_we(1))).or. &
(xtra1(i).gt.real(nx_we(2))))) itra1(i)=-999999999
endif
if (itra1(i).ne.-999999999) numactiveparticles= &
if (WATERCYCLE.eqv..FALSE.) then
if (itra1(i).eq.itime) then
if ((ytra1(i).gt.real(ny_sn(2))).or. &
(ytra1(i).lt.real(ny_sn(1)))) itra1(i)=-999999999
if (((.not.xglobal).or.(nx_we(2).ne.(nx-2))).and. &
((xtra1(i).lt.real(nx_we(1))).or. &
(xtra1(i).gt.real(nx_we(2))))) itra1(i)=-999999999
endif
endif
if (itra1(i).ne.-999999999) numactiveparticles= &
numactiveparticles+1
end do
......@@ -308,6 +310,11 @@ subroutine boundcond_domainfill(itime,loutend)
itramem(ipart)=itra1(ipart)
itrasplit(ipart)=itra1(ipart)+ldirect*itsplit
xmass1(ipart,1)=xmassperparticle
if (WATERCYCLE) then
status_q(ipart)=-1
! write(*,*) 'calculate watercycle in domainfilling',numactiveparticles,ipart
call calculate_watercycle(ipart,itime)
endif
if (mdomainfill.eq.2) xmass1(ipart,1)= &
xmass1(ipart,1)*pvpart*48./29.*ozonescale/10.**9
else
......@@ -540,6 +547,11 @@ subroutine boundcond_domainfill(itime,loutend)
itramem(ipart)=itra1(ipart)
itrasplit(ipart)=itra1(ipart)+ldirect*itsplit
xmass1(ipart,1)=xmassperparticle
if (WATERCYCLE) then
status_q(ipart)=-1
! write(*,*) 'calculate watercycle in domainfilling 1',numactiveparticles,ipart
call calculate_watercycle(ipart,itime)
endif
if (mdomainfill.eq.2) xmass1(ipart,1)= &
xmass1(ipart,1)*pvpart*48./29.*ozonescale/10.**9
else
......
src/calcpar.f90
View file @ ad0a46eb
......@@ -139,6 +139,9 @@ subroutine calcpar(n,uuh,vvh,pvh,metdata_format)
tth(ix,jy,llev,n),ustar(ix,jy,1,n),sshf(ix,jy,1,n),akm,bkm,akz(llev),metdata_format)
else
llev=0
! ol=obukhov(ps(ix,jy,1,n),tt2(ix,jy,1,n),td2(ix,jy,1,n), &
! tth(ix,jy,2,n),ustar(ix,jy,1,n),sshf(ix,jy,1,n),akm,bkm,akz(llev),metdata_format)
! ESO: akz(0) is out of bounds
ol=obukhov(ps(ix,jy,1,n),tt2(ix,jy,1,n),td2(ix,jy,1,n), &
tth(ix,jy,2,n),ustar(ix,jy,1,n),sshf(ix,jy,1,n),akm,bkm,akzdummy,metdata_format)
end if
......
src/calculate_watercycle.f90 0 → 100644
View file @ ad0a46eb
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
! *
! This file is part of FLEXPART. *
! *
! FLEXPART is free software: you can redistribute it and/or modify *
! it under the terms of the GNU General Public License as published by*
! the Free Software Foundation, either version 3 of the License, or *
! (at your option) any later version. *
! *
! FLEXPART is distributed in the hope that it will be useful, *
! but WITHOUT ANY WARRANTY; without even the implied warranty of *
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
! GNU General Public License for more details. *
! *
! You should have received a copy of the GNU General Public License *
! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
!**********************************************************************
subroutine calculate_watercycle(partnumber,itime)
!*****************************************************************************
! *
! calculate difference in q and dump it to the dry/wetdepofield *
! VAR: real, allocatable, dimension(:) :: xtra1_q, ytra1_q, val_q *
! *
!*****************************************************************************
use par_mod
use com_mod
implicit none
real(kind=dp) :: jul
integer :: itime,i,j,k,jjjjmmdd,ihmmss,numshortout,numshortall
integer :: ix,jy,ixp,jyp,ind,indz,indzp,il,m,indexh,itage,nage
integer :: interp_time, ngrid
integer :: partnumber,forparticle,totalparticle ! -1 inititalize q, partnumber real do deposition
real :: dt1,dt2,dtt,ddx,ddy,rddx,rddy,p1,p2,p3,p4,xtn,ytn
real centerposx_real, centerposy_real, centerposx, centerposy, diff, xl(2), yl(2)
real qv1(2),qvprof(2),qvi,dz,dz1,dz2
real e_minus_p1(2), e_minus_pi
real precip_reference(360,180),eminusp_saved
integer ixr,jyr,idumx,idumy,istep, watersynctime
! watersynctime=3600*3 ! was loutstep before, but for monthly run not possible
watersynctime=3600
! watersynctime=loutstep
if ((partnumber.eq.1).and.(itime.eq.0)) write(*,*) 'Watersynctime ',watersynctime,lsynctime
! Some variables needed for temporal interpolation
!*************************************************
dt1=real(itime-memtime(1))
dt2=real(memtime(2)-itime)
dtt=1./(dt1+dt2)
! if called from timemanager a loop over all particle is necessary
! if particle is released, the initialization has just to be done
! for the single particle
!******************************************************************
if (partnumber.eq.-1) then
forparticle=numpart
totalparticle=0
istep=int(itime/watersynctime)
write(*,*) 'Model step:',istep
! open(43,file='precip_ref_input.dat',status='old')
! write(*,*) 'Reading eminusp from file'
! do 75 jyr=1,90
! do 75 ixr=1,360
!75 read(43,'(2i4,e15.3)') idumx,idumy,precip_reference(ixr,jyr)
! close(43)
! open(43,file='precip_ref_plus_eminusp.dat')
! write(*,*) 'Writing eminusp from file'
! do 76 jyr=1,90
! do 76 ixr=1,360
!6 write(43,'(2i4,3e15.3)') ixr,jyr,precip_reference(ixr,jyr), &
! e_minus_p(ixr,jyr,1), e_minus_p(ixr,jyr+90,1)
! close(43)
! write(*,*) 'Written eminusp from file'
else
forparticle=1
endif
! Loop about all particles
!*************************
do j=1,forparticle
if (partnumber.eq.-1) then
i=j ! loop over all particle
else
i=partnumber ! just for one particle
endif
! write(*,*) 'particle: ',i,xtra1_q(i),xtra1(i),itime,itra1(i),itramem(i),status_q(i)
! Take only valid particles, itime invalid when particle released!
!*****************************************************************
if (itra1(i).eq.itime) then ! trajectory itra1 is not -999999
! Determine age class of the particle - nage is used for the kernel
!******************************************************************
itage=abs(itra1(i)-itramem(i))
do nage=1,nageclass
if (itage.lt.lage(nage)) goto 33
end do
33 continue
! Determine which nesting level to be used
!*****************************************
ngrid=0
do k=numbnests,1,-1
if ((xtra1(i).gt.xln(k)).and.(xtra1(i).lt.xrn(k)).and. &
(ytra1(i).gt.yln(k)).and.(ytra1(i).lt.yrn(k))) then
ngrid=k
goto 23
endif
end do
23 continue
! Determine nested grid coordinates
!**********************************
if (ngrid.gt.0) then
xtn=(xtra1(i)-xln(ngrid))*xresoln(ngrid)
ytn=(ytra1(i)-yln(ngrid))*yresoln(ngrid)
ix=int(xtn)
jy=int(ytn)
ddy=ytn-real(jy)
ddx=xtn-real(ix)
else
ix=int(xtra1(i))
jy=int(ytra1(i))
ddx=xtra1(i)-real(ix)
ddy=ytra1(i)-real(jy)
endif
ixp=ix+1
jyp=jy+1
if (jy.eq.numygrid) jyp=numygrid
rddx=1.-ddx
rddy=1.-ddy
p1=rddx*rddy
p2=ddx*rddy
p3=rddx*ddy
p4=ddx*ddy
! Interpolate specific humidity - code from partoutput.f90
!*********************************************************
do il=2,nz
if (height(il).gt.ztra1(i)) then
indz=il-1
indzp=il
goto 6
endif
end do
6 continue
dz1=ztra1(i)-height(indz)
dz2=height(indzp)-ztra1(i)
dz=1./(dz1+dz2)
do ind=indz,indzp
do m=1,2
indexh=memind(m)
if (ngrid.eq.0) then
qv1(m)=p1*qv(ix ,jy ,ind,indexh) &
+p2*qv(ixp,jy ,ind,indexh) &
+p3*qv(ix ,jyp,ind,indexh) &
+p4*qv(ixp,jyp,ind,indexh)
else
qv1(m)=p1*qvn(ix ,jy ,ind,indexh,ngrid) &
+p2*qvn(ixp,jy ,ind,indexh,ngrid) &
+p3*qvn(ix ,jyp,ind,indexh,ngrid) &
+p4*qvn(ixp,jyp,ind,indexh,ngrid)
endif
end do
qvprof(ind-indz+1)=(qv1(1)*dt2+qv1(2)*dt1)*dtt
end do
qvi=(dz1*qvprof(2)+dz2*qvprof(1))*dz
if ( (mod(itage,watersynctime).eq.0).and.(itage.gt.0) ) then !trajectory is loutstep after simulation start, not a new particle
do m=1,2
indexh=memind(m)
if (ngrid.eq.0) then
e_minus_p1(m)=p1*e_minus_p(ix ,jy ,indexh) &
+p2*e_minus_p(ixp,jy ,indexh) &
+p3*e_minus_p(ix ,jyp,indexh) &
+p4*e_minus_p(ixp,jyp,indexh)
else
e_minus_p1(m)=p1*e_minus_p_nest(ix ,jy ,indexh,ngrid) &
+p2*e_minus_p_nest(ixp,jy ,indexh,ngrid) &
+p3*e_minus_p_nest(ix ,jyp,indexh,ngrid) &
+p4*e_minus_p_nest(ixp,jyp,indexh,ngrid)
endif
end do
e_minus_pi=(e_minus_p1(1)*dt2+dt1*e_minus_p1(2))*dtt
eminusp_saved=e_minus_pi
! Calculate humidity differences
!*******************************
diff=(qvi-val_q(i))*xmass1(i,1)*ldirect
! Determine center of mass position on earth and average height
!**************************************************************
xl(1)=xlon0+xtra1_q(i)*dx
xl(2)=xlon0+xtra1(i)*dx
yl(1)=ylat0+ytra1_q(i)*dy
yl(2)=ylat0+ytra1(i)*dy
! Do this only for the 2nd timestep
! there has to be a deltaq which is negative and the column has to have precipitation
! e_minus_p units is mm/day .. convert to mm/h
! e_minus_pi=precip_reference(ix+1,jy-90+2)*(-24.)
! First timestep and criteria fullfilled, keep it or throw it?
if ( ((diff.ge.0).or.((e_minus_pi/24).gt.-.5)) .and. (status_q(i).eq.1) ) then
! if ( ((diff.ge.0).or.((e_minus_pi/24).gt.-2.0)) .and. (status_q(i).eq.1) ) then
! if ( ((diff.ge.0)) .and. (status_q(i).eq.1) ) then
! if (status_q(i).eq.-100) then ! never! all trajectories.
! write(*,*) 'terminated: ',i,diff,e_minus_pi,itage,itramem(i),status_q(i),qvi,val_q(i),xmass1(i,1)
status_q(i)=-9 ! this particle is invalid
itra1(i)=-999999999
else
totalparticle=totalparticle+1
call centerofmass(xl,yl,2,centerposx_real,centerposy_real)
centerposx=(centerposx_real-xlon0)/dx
centerposy=(centerposy_real-ylat0)/dy
! write(*,*) 'accounted: ',i,itime,itra1(i),itage,diff,itramem(i),status_q(i)
if (diff.gt.0) then
call drydepokernel(1,diff,centerposx,centerposy,nage,1)
if (nested_output.eq.1) &
call drydepokernel_nest(1,diff,centerposx,centerposy,nage,1)
else
call wetdepokernel(1,abs(diff),centerposx,centerposy,nage,1)
if (nested_output.eq.1) &
call wetdepokernel_nest(1,abs(diff),centerposx,centerposy,nage,1)
endif
! if (istep.lt.100) then
! waterfieldp(istep,ix+1,jy+1)=waterfieldp(istep,ix+1,jy+1)+xmass1(i,1)
! endif
! waterfielde(ix+1,jy+1)=waterfielde(ix+1,jy+1)+xmass1(i,1)
xtra1_q(i)=xtra1(i)
ytra1_q(i)=ytra1(i)
val_q(i)=qvi
status_q(i)=2
endif
endif ! at louttimestep interval
if (status_q(i).eq.-1) then ! first time, a value is saved
! write(*,*) 'saved: ',i,xtra1(i),ytra1(i),qvi,itime,itage,itramem(i),status_q(i)
xtra1_q(i)=xtra1(i)
ytra1_q(i)=ytra1(i)
val_q(i)=qvi
status_q(i)=1
endif
endif ! for valid trajectories - itra1.eq.itime
end do
if ((partnumber.eq.-1).and.(totalparticle.gt.0)) then
write (*,*) 'calculated watercycle, ',totalparticle,numpart,partnumber,forparticle,itime
endif
end subroutine calculate_watercycle
src/com_mod.f90
View file @ ad0a46eb
......@@ -341,6 +341,9 @@ module com_mod
real :: lsm(0:nxmax-1,0:nymax-1)
real :: xlanduse(0:nxmax-1,0:nymax-1,numclass)
real :: waterfieldp(100,360,180),waterfielde(360,180)
! oro [m] orography of the ECMWF model
! excessoro excess orography mother domain
! lsm land sea mask of the ECMWF model
......@@ -495,6 +498,8 @@ module com_mod
real,allocatable,dimension(:,:,:,:,:) :: uun, vvn, wwn, ttn, qvn, pvn,&
& rhon, drhodzn, tthn, qvhn, clwcn, ciwcn, clwn, clwchn, ciwchn
real,allocatable,dimension(:,:,:,:) :: ctwcn
real,allocatable, dimension (:,:,:,:) :: uul,vvl
real,allocatable, dimension (:,:,:,:,:) :: uuln,vvln
integer,allocatable,dimension(:,:,:,:) :: cloudshn
integer(kind=1),allocatable,dimension(:,:,:,:,:) :: cloudsn
......@@ -583,6 +588,7 @@ module com_mod
logical :: DEP,DRYDEP,DRYDEPSPEC(maxspec),WETDEP,WETDEPSPEC(maxspec),&
& OHREA,ASSSPEC
logical :: DRYBKDEP,WETBKDEP
logical :: WATERCYCLE
! numxgrid,numygrid number of grid points in x,y-direction
! numxgridn,numygridn number of grid points in x,y-direction for nested output grid
......@@ -604,6 +610,7 @@ module com_mod
! ASSSPEC .true., if there are two species asscoiated
! DRYBKDEP,WETBKDEP .true., for bkwd runs, where mass deposited and source regions is calculated - either for dry or for wet deposition
! (i.e. transfer of mass between these two occurs
! WATERCYCLE .true. bkwd domainfilling run, which gives qdiff fields in the dry and wetdeposition fields
......@@ -675,6 +682,11 @@ module com_mod
real, allocatable, dimension(:,:) :: xmass1
real, allocatable, dimension(:,:) :: xscav_frac1
real, allocatable, dimension(:,:,:) :: e_minus_p ! E-P field for the watercycle
real, allocatable, dimension(:,:,:,:) :: e_minus_p_nest ! E-P field for the watercycle
real, allocatable, dimension(:) :: xtra1_q, ytra1_q, val_q
integer, allocatable, dimension(:) :: status_q
! eso: Moved from timemanager
real, allocatable, dimension(:) :: uap,ucp,uzp,us,vs,ws
integer(kind=2), allocatable, dimension(:) :: cbt
......
src/concoutput.f90
View file @ ad0a46eb
......@@ -281,7 +281,7 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
do ix=0,numxgrid-1
! WET DEPOSITION
if ((WETDEP).and.(ldirect.gt.0)) then
if ((WETDEP).and.(ldirect.gt.0).or.(WATERCYCLE)) then
do l=1,nclassunc
auxgrid(l)=wetgridunc(ix,jy,ks,kp,l,nage)
end do
......@@ -300,7 +300,7 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
endif
! DRY DEPOSITION
if ((DRYDEP).and.(ldirect.gt.0)) then
if ((DRYDEP).and.(ldirect.gt.0).or.(WATERCYCLE)) then
do l=1,nclassunc
auxgrid(l)=drygridunc(ix,jy,ks,kp,l,nage)
end do
......@@ -359,11 +359,12 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
sp_count_r=0
sp_fact=-1.
sp_zer=.true.
if ((ldirect.eq.1).and.(WETDEP)) then
if (((ldirect.eq.1).and.(WETDEP)).or.(WATERCYCLE)) then
do jy=0,numygrid-1
do ix=0,numxgrid-1
!oncentraion greater zero
if (wetgrid(ix,jy).gt.smallnum) then
! write(*,*) ix,jy,wetgrid(ix,jy)
if (sp_zer.eqv..true.) then ! first non zero value
sp_count_i=sp_count_i+1
sparse_dump_i(sp_count_i)=ix+jy*numxgrid
......@@ -371,8 +372,13 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
sp_fact=sp_fact*(-1.)
endif
sp_count_r=sp_count_r+1
if (WATERCYCLE) then
sparse_dump_r(sp_count_r)= &
sp_fact*wetgrid(ix,jy)/area(ix,jy)
else
sparse_dump_r(sp_count_r)= &
sp_fact*1.e12*wetgrid(ix,jy)/area(ix,jy)
endif
! sparse_dump_u(sp_count_r)=
!+ 1.e12*wetgridsigma(ix,jy,ks,kp,nage)/area(ix,jy)
else ! concentration is zero
......@@ -396,7 +402,7 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
sp_count_r=0
sp_fact=-1.
sp_zer=.true.
if ((ldirect.eq.1).and.(DRYDEP)) then
if (((ldirect.eq.1).and.(DRYDEP)).or.(WATERCYCLE)) then
do jy=0,numygrid-1
do ix=0,numxgrid-1
if (drygrid(ix,jy).gt.smallnum) then
......@@ -407,9 +413,14 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
sp_fact=sp_fact*(-1.)
endif
sp_count_r=sp_count_r+1
if (WATERCYCLE) then
sparse_dump_r(sp_count_r)= &
sp_fact*drygrid(ix,jy)/area(ix,jy)
else
sparse_dump_r(sp_count_r)= &
sp_fact* &
1.e12*drygrid(ix,jy)/area(ix,jy)
endif
! sparse_dump_u(sp_count_r)=
!+ 1.e12*drygridsigma(ix,jy,ks,kp,nage)/area(ix,jy)
else ! concentration is zero
......@@ -438,7 +449,7 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
do kz=1,numzgrid
do jy=0,numygrid-1
do ix=0,numxgrid-1
if (grid(ix,jy,kz).gt.smallnum) then
if ((grid(ix,jy,kz).gt.smallnum).and.(.not.WATERCYCLE)) then
if (sp_zer.eqv..true.) then ! first non zero value
sp_count_i=sp_count_i+1
sparse_dump_i(sp_count_i)= &
......@@ -644,26 +655,10 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
! Reinitialization of grid
!*************************
! do ks=1,nspec
! do kp=1,maxpointspec_act
! do i=1,numreceptor
! creceptor(i,ks)=0.
! end do
! do jy=0,numygrid-1
! do ix=0,numxgrid-1
! do l=1,nclassunc
! do nage=1,nageclass
! do kz=1,numzgrid
! gridunc(ix,jy,kz,ks,kp,l,nage)=0.
! end do
! end do
! end do
! end do
! end do
! end do
! end do
creceptor(:,:)=0.
gridunc(:,:,:,:,:,:,:)=0.
if (WATERCYCLE) then
wetgridunc(:,:,:,:,:,:)=0
drygridunc(:,:,:,:,:,:)=0
endif
end subroutine concoutput
src/drydepokernel.f90
View file @ ad0a46eb
......@@ -84,7 +84,7 @@ subroutine drydepokernel(nunc,deposit,x,y,nage,kp)
if (.not.lusekerneloutput) then
do ks=1,nspec
if ((abs(deposit(ks)).gt.0).and.DRYDEPSPEC(ks)) then
if ((abs(deposit(ks)).gt.0).and.DRYDEPSPEC(ks).or.WATERCYCLE) then
if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
(jy.le.numygrid-1)) then
drygridunc(ix,jy,ks,kp,nunc,nage)= &
......@@ -99,7 +99,7 @@ subroutine drydepokernel(nunc,deposit,x,y,nage,kp)
!**********************************************
do ks=1,nspec
if ((abs(deposit(ks)).gt.0).and.DRYDEPSPEC(ks)) then
if ((abs(deposit(ks)).gt.0).and.DRYDEPSPEC(ks).or.WATERCYCLE) then
if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
(jy.le.numygrid-1)) then
......
src/euler_rain.f90 0 → 100644
View file @ ad0a46eb
subroutine euler_rain(itime)
!*******************************************************************************
! Compute E-P between 2 time steps from moisture budget in an atmospheric column
! 2 Terms: 1. Moisture column tendency between the 2 time steps, *
! 2. Divergence of time-mean wind between the 2 time steps *
!*******************************************************************************
! Variables: *
! *
! uul, vvl horiz. wind components on original eta levels *
! *
! *
!*******************************************************************************
use par_mod
use com_mod
integer itime
integer ix,jy,kz,k,n,ixp,jyp,ixm,jym
real pih,ylat,ylatp,ylatm
parameter(pih=pi/180.)
real xmasscolumn(2),pint(nzmax,2)
real pdiff(0:nxmax-1,0:nymax-1,nzmax,2),gradx(2),grady(2),gradterm
real divx,divy,div(nzmax,3),divwater(nzmax,3)
real divxwater,divywater,gradtermwater
real watercolumn(2),watermass(nzmax,2)
real dwatercolumndt(0:nxmax-1,0:nymax-1)
real divwatercolumn(0:nxmax-1,0:nymax-1,3)
! write(89,*) xlon0+dx,ylat0+dy,nx-2,ny-2,dx,dy,1,10000.,3
! Loop over the whole grid
!*************************
do 10 jy=0,ny-1
do 10 ix=0,nx-1
!*****************************************************************
! 1. Compute the precipitable water tendency in atmospheric column
! and the fluxes needed for the divergence term (step 2)
!*****************************************************************
do 20 k=1,2
n=memind(k)
watercolumn(n)=0.
! Compute mass of total atmospheric column (xmasscolumn)
! and water vapor mass in atmospheric column (watercolumn)
!*********************************************************
xmasscolumn(n)=ps(ix,jy,1,n)/ga
do 21 kz=1,nwz
21 pint(kz,n)=akm(kz)+bkm(kz)*ps(ix,jy,1,n)
do 20 kz=2,nwz
pdiff(ix,jy,kz,n)=pint(kz-1,n)-pint(kz,n)
watermass(kz,n)=qvh(ix,jy,kz,n)*pdiff(ix,jy,kz,n)/ga
20 watercolumn(n)=watercolumn(n)+watermass(kz,n)
! Compute the tendency
!*********************
dwatercolumndt(ix,jy)=(watercolumn(memind(2))- &
watercolumn(memind(1)))/float(memtime(2)-memtime(1))
10 continue
!**************************************************
! 2. Compute the precipitable water divergence term
!**************************************************