program forcing

implicit none

! [Arguments]
! Constants
integer,parameter::r8 = selected_real_kind(12)
integer,parameter::regi = 256
integer,parameter::reso = 64
integer,parameter::nt = 1440
!integer,parameter::nt = 1
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) q1(nx,ny,nt,nz)
real(r8) t1(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) cape0(nx,ny,nt)
real(r8) cape00(nx,ny,nt)
real(r8) cape(nx,ny,nt)
real(r8) dcape(nx,ny,nt)
real(r8) negcape(nx,ny,nt)
real(r8) dcape0(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 readnc1("/home/wangxu/WORK4/CASPO/mb_tao_no_shallow/dcape.nc","dcape",nt,dcape0)

do i = 1,nt
    print *,"time: ",i
    write(ii,"(I4)") i
    call readnc0("/home/wangxu/WORK3/tao/m64/Qv/Qv."//trim(adjustl(ii))//adjustl(".nc"),"Qv",nz,ny,nx,q(:,:,i,:))
    call readnc0("/home/wangxu/WORK3/tao/m64/T/T."//trim(adjustl(ii))//adjustl(".nc")  ,"T" ,nz,ny,nx,t(:,:,i,:))
    call readnc0("/home/wangxu/WORK3/tao/m64/p/p."//trim(adjustl(ii))//adjustl(".nc")  ,"p" ,nz,ny,nx,p(:,:,i,:))
    call readnc0("/home/wangxu/WORK3/tao/m64/z/z."//trim(adjustl(ii))//adjustl(".nc")  ,"z" ,nz,ny,nx,z(:,:,i,:))
    call readnc0("/home/wangxu/WORK3/tao/m64/w/w."//trim(adjustl(ii))//adjustl(".nc")  ,"w" ,nz,ny,nx,w(:,:,i,:))
   
    call readnc0("/home/wangxu/WORK3/tao/m64/pf/pf."//trim(adjustl(ii))//adjustl(".nc"),"p",nz+1,ny,nx,pf(:,:,i,:))
    call readnc0("/home/wangxu/WORK3/tao/m64/zf/zf."//trim(adjustl(ii))//adjustl(".nc"),"z",nz+1,ny,nx,zf(:,:,i,:))
end do
!print *,"OK ... 1"
!print *,"OK ... 2"
do i = 1,nt
! dcape
do i1 = 1,ny
    do i2 = 1,nx
        if(i.ne.nt) then
            q1(i2,i1,i,:) = q(i2,i1,i+1,:)
            t1(i2,i1,i,:) = t(i2,i1,i+1,:)
        else
            q1(i2,i1,i,:) = q(i2,i1,i,:)
            t1(i2,i1,i,:) = t(i2,i1,i,:)
        end if
        call zm_convr(q(i2,i1,i,:),t(i2,i1,i,:),p(i2,i1,i,:),z(i2,i1,i,:), &
                      pf(i2,i1,i,:),zf(i2,i1,i,:),q1(i2,i1,i,:),t1(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),mb0(i2,i1,i),cape00(i2,i1,i),negcape(i2,i1,i),dcape0)
    end do
end do    
end do
print *,"OK ... 3"
! Write
call writenc1("data/dcape_e.nc","dcape","",0.,nt,ny,nx,dcape)
call writenc1("data/cin.nc","cin","",0.,nt,ny,nx,negcape)

print *, "Good !"
end program forcing
