program forcing

implicit none

! [Arguments]
! Constants
integer,parameter::r8 = selected_real_kind(12)
integer,parameter::regi = 512
integer,parameter::reso = 64
integer,parameter::nt = 160
integer,parameter::nz = 40
integer,parameter::ny = regi/reso
integer,parameter::nx = regi/reso

real(r8),parameter::cp = 1004.64_r8
real(r8),parameter::g  = 9.80616_r8
! Variables
real(r8) q(nx,ny,nt,nz)
real(r8) t(nx,ny,nt,nz)
real(r8) p(nx,ny,nt,nz)
real(r8) z(nx,ny,nt,nz)
real(r8) w(nx,ny,nt,nz)

real(r8) pf(nx,ny,nt,nz+1)
real(r8) zf(nx,ny,nt,nz+1)

real(r8) utx(nx,ny,nt,nz)
real(r8) vty(nx,ny,nt,nz)
real(r8) wtz(nx,ny,nt,nz)
real(r8) uqx(nx,ny,nt,nz)
real(r8) vqy(nx,ny,nt,nz)
real(r8) wqz(nx,ny,nt,nz)
real(r8) ptendt(nx,ny,nt,nz)
real(r8) ptendq(nx,ny,nt,nz)

real(r8) cape0(nx,ny,nt)
real(r8) cape00(nx,ny,nt)
real(r8) cape(nx,ny,nt)
real(r8) dcape(nx,ny,nt)
real(r8) mb(nx,ny,nt)
real(r8) mb0(nx,ny,nt)
real(r8) dadt(nx,ny,nt)

integer jb(nx,ny,nt)
integer jt0(nx,ny,nt)
integer jt(nx,ny,nt)
integer jt3(nx,ny,nt)

integer i
integer i1
integer i2
integer k

character(10)::ii
! [Routine]
! Forcing
call readnc0("/home/wangxu/WORK3/jiwen/m64/qv.nc","qv",nt,nz,ny,nx,q)
call readnc0("/home/wangxu/WORK3/jiwen/m64/tk.nc","tk",nt,nz,ny,nx,t)
call readnc0("/home/wangxu/WORK3/jiwen/m64/p.nc" ,"p" ,nt,nz,ny,nx,p)
call readnc0("/home/wangxu/WORK3/jiwen/m64/z.nc" ,"z" ,nt,nz,ny,nx,z)
call readnc0("/home/wangxu/WORK3/jiwen/m64/w.nc" ,"w" ,nt,nz,ny,nx,w)

call readnc0("/home/wangxu/WORK3/jiwen/m64/pf.nc","p" ,nt,nz+1,ny,nx,pf)
call readnc0("/home/wangxu/WORK3/jiwen/m64/zf.nc","z" ,nt,nz+1,ny,nx,zf)

call readnc0("/home/wangxu/WORK3/jiwen/m64/utx.nc","utx",nt,nz,ny,nx,utx)
call readnc0("/home/wangxu/WORK3/jiwen/m64/vty.nc","vty",nt,nz,ny,nx,vty)
call readnc0("/home/wangxu/WORK3/jiwen/m64/wtz.nc","wtz",nt,nz,ny,nx,wtz)

call readnc0("/home/wangxu/WORK3/jiwen/m64/uqx.nc","uqx",nt,nz,ny,nx,uqx)
call readnc0("/home/wangxu/WORK3/jiwen/m64/vqy.nc","vqy",nt,nz,ny,nx,vqy)
call readnc0("/home/wangxu/WORK3/jiwen/m64/wqz.nc","wqz",nt,nz,ny,nx,wqz)
print *,"OK ... 1"
! ptendt, ptendq
ptendt = -1.*(utx+vty+wtz+w*g/cp)
ptendq = -1.*(uqx+vqy+wqz)
print *,"OK ... 2"
! dcape
do i = 1,nt
    print*,"time:",i
    do i1 = 1,ny
        do i2 = 1,nx
            call zm_convr(q(i2,i1,i,:),t(i2,i1,i,:),p(i2,i1,i,:)/100.,z(i2,i1,i,:), &
                          pf(i2,i1,i,:)/100.,zf(i2,i1,i,:),ptendq(i2,i1,i,:),ptendt(i2,i1,i,:), &
                          cape0(i2,i1,i),cape(i2,i1,i),dcape(i2,i1,i), &
                          jb(i2,i1,i),jt0(i2,i1,i),jt(i2,i1,i),mb(i2,i1,i),dadt(i2,i1,i),jt3(i2,i1,i),mb0(i2,i1,i),cape00(i2,i1,i))
        end do    
    end do
end do    
print *,"OK ... 3"
! Write
call writenc1("data/cape0.nc","cape","",0.,nt,ny,nx,cape0)
call writenc1("data/cape.nc","cape","",0.,nt,ny,nx,cape00)
call writenc1("data/dcape.nc","dcape","",0.,nt,ny,nx,dcape)
call writenc1("data/mb.nc","mb","",0.,nt,ny,nx,mb)
call writenc1("data/mb0.nc","mb","",0.,nt,ny,nx,mb0)
call writenc1("data/dadt.nc","dadt","",0.,nt,ny,nx,dadt)

call writenc2("data/jb.nc","jb","",0,nt,ny,nx,jb)
call writenc2("data/jt.nc","jt","",0,nt,ny,nx,jt)
call writenc2("data/jt3.nc","jt","",0,nt,ny,nx,jt3)

print *, "Good !"
end program forcing
