module met 
    use const
    use comm
    use var
implicit real*8 (a-h,o-z)
#   include "pnetcdf.inc"

contains

subroutine metgen
    dimension fg(j2t),tanphig(j2t),		  &
     yg(j0t),yvg(j0t),yvdegg(j0t),ydegg(j0t),     &
     csg(j0t),csvg(j1t),ocsg(j0t),ocsvg(j1t),     &
     dxg(j0t),dxvg(j1t),odxg(j0t),odxvg(j1t),     &
     dyg(j0t),dyvg(j1t),odyg(j0t),odyvg(j1t)      


    real*8 dd,ddv,flag  !! kilmer

    ! tlz= depth of bottom z-level
    ! ztop= depth of interface between top two layers
    ! dydx= ratio of dy to dx at each latitude (constant)
    ! dxmnut= longitude resolution (minutes)
    ! r0= earth's radius (cm)
    ! dxdeg= longitudinal resolution (degrees)
    data tlz_depth,d,ztop,dydx,dxmnut,y0deg,r0/ &
      5800.E2,.04, 20.E2,  1.,   60.,-86.7,6.4E8/  ! kilmer

    dxmnut=360.*60./dble(i2t)
    dxdeg=dxmnut/60.
    !dydx=dydx/dxdeg


    westdeg=360./dble(ngx0)
    westdeg=westdeg*dble(mylon)
    if (westdeg .ge. 360.) westdeg=westdeg-360.

!    open(9,file='data/westdeg'//grd)
!    rewind 9
!    write(9,'(f9.4)') westdeg
!    close(9)
    call wr8nc_westdeg !dog
!    open(18,file='data/latitudes'//grd,form='unformatted')
!    open(10,file='data/yzgrid'//grd,form='unformatted')
!    rewind 10
!    rewind 18

    tmp=.01*tlz_depth
    dzz=tmp/(2.*dble(k1))
    coeff_c=.001
    coeff_d=d
    coeff_a=tmp*(1.-coeff_d)/(1.-exp(coeff_c*tmp))
    coeff_b=-coeff_a
    zz=0.
    do k=1,k0+k1
	z(k)=100.*(coeff_a+coeff_b*exp(coeff_c*zz)+coeff_d*zz)
	zz=zz+dzz
    enddo	     
    if(myid.eq.0)write(*,'(10f8.0)') z

    ! calculate vertical metrics
    do k=1,k1
	odz(k)=1./(z(2*k+1)-z(2*k-1))
    enddo	     
    ! odzw(1,i,j) and odzw(k0,i,j) are not used
    do k=2,k1
	odzw(k)=1./(z(2*k)-z(2*k-2))
    enddo	     


    it=mylon*i2
    jt=mylat*j2

    ! degrees to radians
    pi_180=3.141592654/180.
    yvdegg(1)=y0deg
    ydegg(1)=y0deg
    yvg(1)=0.
    dyg(1)=0.
    odyg(1)=0.
    ! equatorial dx
    dx0=r0*dxdeg*pi_180
    ! kilmer
    !dd=21.0/100.0
    !ddv=dd/100.0
    !flag=1.0
    ! kilmer end
    do j=2,j0t
	! initial guess
	!      yvdegg(j)=yvdegg(j-1)+dxdeg*dydx*cos(yvdegg(j-1)*pi_180)
	yvdegg(j)=yvdegg(j-1)+1.476-dxdeg*dydx*cos(yvdegg(j-1)*pi_180) ! kilmer
	!      print*,j,yvdegg(j)
	dxg(j)=dx0*cos(.5*(yvdegg(j-1)+yvdegg(j))*pi_180)

	tmp=yvdegg(j)
	!      yvdegg(j)=yvdegg(j-1)+dydx*dxg(j)/(r0*pi_180)
	yvdegg(j)=yvdegg(j-1)+1.476-dydx*dxg(j)/(r0*pi_180) ! kilmer
	ydegg(j)=.5*(yvdegg(j-1)+yvdegg(j))
	csg(j)=cos(ydegg(j)*pi_180)
	dxrad=dxdeg*pi_180*csg(j)
	dxg(j)=dxrad*r0
	dyrad=(yvdegg(j)-yvdegg(j-1))*pi_180
	temp=dyrad/dxrad
	dydeg=yvdegg(j)-yvdegg(j-1)

	dyg(j)=dydeg*r0*pi_180
	dyvg(j-1)=(ydegg(j)-ydegg(j-1))*r0*pi_180
	odyg(j)=1./dyg(j)
	odyvg(j-1)=1./dyvg(j-1)
	yvg(j)=yvg(j-1)+dyg(j)
	yg(j)=.5*(yvg(j)+yvg(j-1))
	ocsg(j)=1./csg(j)
	odxg(j)=1./dxg(j)
	csvg(j-1)=cos(pi_180*yvdegg(j-1))
	ocsvg(j-1)=1./csvg(j-1)
	dxvg(j-1)=dx0*csvg(j-1)
	odxvg(j-1)=1./dxvg(j-1)
	! kilmer
	!odxvg(j-1)=1./dxvg(j-1)
	!if (j>j0t-100) then
	!    dd=dd+ddv
	!    flag=0.0
	!endif
	!if (flag==1.0) dd=dd-ddv
	!if (dd<0.0) dd=0.0
    enddo
    ! kilmer end
    y1deg=yvdegg(j1)


    jt=mylat*j2
    do j=1,j1
	csv(j)=csvg(jt+j)
	ocsv(j)=ocsvg(jt+j)
	dxv(j)=dxvg(jt+j)
	odxv(j)=odxvg(jt+j)
	dyv(j)=dyvg(jt+j)
	odyv(j)=odyvg(jt+j)
    enddo
    do j=1,j0
	y(j)=yg(jt+j)
	yv(j)=yvg(jt+j)
	yvdeg(j)=yvdegg(jt+j)
	ydeg(j)=ydegg(jt+j)
	cs(j)=csg(jt+j)
	ocs(j)=ocsg(jt+j)
	dx(j)=dxg(jt+j)
	odx(j)=odxg(jt+j)
	dy(j)=dyg(jt+j)
	ody(j)=odyg(jt+j)
    enddo
!    open(16,file='data/cell_centers'//grd)
     if (myid ==0) write(*,16) (j,ydeg(j),j=2,j1)
    16   format(' j ','cell center latitude'/(i3,f11.4))
!    write(18) ydeg

    ! ----------------------------
    ! set coriolis parameter array
    ! ----------------------------
    omega2=3.141592654/21600.
    do j=2,j1
	! spherical curvature parameter
	tanphi(j-1)=tan(pi_180*ydeg(j))/r0
	f(j-1)=omega2*sin(pi_180*ydeg(j))
    enddo
    do i=1,i0
	xdeg(i)=westdeg+(i-1.5)*(dxmnut/60.)
	!w         xdeg(i)=(i-1.5)*(dxmnut/60.)
	if(xdeg(i).ge.360.)then
	    xdeg(i)=xdeg(i)-360.
	    elseif(xdeg(i).lt.0.)then
	    xdeg(i)=xdeg(i)+360.
	endif
    enddo
    write(*,'(10f9.4)') (xdeg(i),i=2,i1)

!    write(10)  j0,k0
    write(*,*) j0,k0

!    write(10) tlz_depth,z,odz,odzw,dxmnut,xdeg,y,yv,yvdeg,ydeg,cs,csv,ocs,ocsv,dx,dxv,odx,odxv,dy,dyv,ody,odyv
    call wr8nc_yzgrid(tlz_depth,dxmnut) !dog
!    open(11,file='data/ydeg'//grd,form='unformatted')
!    write(11) ydeg
    call wr8nc_ydeg !dog
    ! kilmer
    !if (myid==0) then
    !	open(3333,file='data/global_ydegg.log',form='formatted')
    !	write(3333,322) (ydegg(j),j=1,j0t)
    !	close(3333)
    !endif
    !322  format(/'global ydegg'/(10f10.4))
    ! kilmer end
    if(mylon.eq.0)then
	write(*,20) (yvdeg(j),j=1,j1)
	20   format(/'model latitude grid lines'/(10f8.4))
	write(*,21) f
	21   format('coriolis parameter'/(1p,10e9.2))
	write(*,99) (j,yvdeg(j),j=1,j1)
	99   format('j,yvdeg(j)=',i3,f9.4)
    endif
end subroutine

subroutine wr8nc_westdeg
use const, only:ngx0,ngy0
use comm, only: mylon,myid,m_cart
 implicit none
integer :: rank,ncid,ierr,nout
 INTEGER(KIND=MPI_OFFSET_KIND) :: start(1),count(1),bufcount=1,dimm
 integer :: x_dim
 integer :: dims=1,dimids(1),varid
 real*8 :: mydata(1)

dimm = ngx0*ngy0
    mydata(1)=360./real(ngx0)
    mydata(1)=mydata(1)*real(mylon)
    if (mydata(1) > 360.) mydata(1)=mydata(1)-360.

 ierr=nfmpi_create(m_cart,'data/westdeg.nc',nf_64bit_offset,mpi_info_null,ncid)
 nout = nfmpi_def_dim(ncid,"myid+1",dimm,x_dim)
 dimids = (/x_dim/)
 nout = nfmpi_def_var(ncid,"westdeg",NF_DOUBLE,dims,dimids,varid)
 nout = nfmpi_enddef(ncid)
 start = (/myid+1/)
 count = (/1/)
 nout = nfmpi_put_vara_all(ncid,varid,start,count,mydata,bufcount,MPI_DOUBLE)
 nout = nfmpi_close(ncid)

end subroutine wr8nc_westdeg

subroutine wr8nc_ydeg
use const, only: ngy0,j0
use comm, only: myid,m_cart
implicit none
character :: fn*12
integer :: ncid,ierr,ns
integer :: j0id
integer :: iydeg
integer(kind=mpi_offset_kind) :: dim_nn
integer(kind=mpi_offset_kind) :: ns_i,nc_i
      !output file name
      fn = 'data/ydeg.nc'
      !create ncfile and id
      ierr=nfmpi_create(m_cart,fn,nf_64bit_offset,mpi_info_null,ncid)
      !define dimension and assign dimID
      dim_nn = j0*ngy0
      ierr=nfmpi_def_dim(ncid,"j0_global",dim_nn,j0id)

      !define variable
      ns = j0id
      ierr=nfmpi_def_var(ncid,"ydeg",nf_double,1,ns,iydeg)
      ierr=nfmpi_enddef(ncid)
      !put variable
      ns_i = mylat*j0+1
      nc_i = j0
      ierr=nfmpi_put_vara_double_all(ncid,iydeg,ns_i,nc_i,ydeg)

      ierr=nfmpi_close(ncid)

end subroutine wr8nc_ydeg

subroutine wr8nc_yzgrid(tlz_depth,dxmnut)
use const, only: i0,j0,j1,k0,k1,ngx0,ngy0
use comm, only: m_cart,mylat,mylon
implicit none

!declare pnetcdf var
character :: fn*14
integer :: ierr,ncid,ns
integer :: siid,k0id,k1id,kkid,i0id,j0id,j1id
integer :: itlz,idxm,iz,iodz,iodzw,ixdeg,iij,iik
integer :: iy,iyv,iyvd,iydeg,ics,iocs,icsv,iocsv,idx,iodx
integer :: idxv,iodxv,idy,idyv,iody,iodyv
integer(kind=mpi_offset_kind) :: dim_nn(7)
integer(kind=mpi_offset_kind) :: ns_i,nc_i
real*8 :: tlz_depth,dxmnut
integer :: tmp(1)
      !output file name
      fn = 'data/yzgrid.nc'
      !create ncfile and id
      ierr=nfmpi_create(m_cart,fn,nf_64bit_offset,mpi_info_null,ncid)

      !define dimension and assign dimID
      dim_nn(1) = 1
      dim_nn(2) = k0
      dim_nn(3) = k1
      dim_nn(4) = k0+k1
      dim_nn(5) = i0*ngx0
      dim_nn(6) = j0*ngy0
      dim_nn(7) = j1*ngy0
      ierr=nfmpi_def_dim(ncid,"pt_val",dim_nn(1),siid)
      ierr=nfmpi_def_dim(ncid,"k0",dim_nn(2),k0id)
      ierr=nfmpi_def_dim(ncid,"k1",dim_nn(3),k1id)
      ierr=nfmpi_def_dim(ncid,"k0+k1",dim_nn(4),kkid)
      ierr=nfmpi_def_dim(ncid,"i0_global",dim_nn(5),i0id)
      ierr=nfmpi_def_dim(ncid,"j0_global",dim_nn(6),j0id)
      ierr=nfmpi_def_dim(ncid,"j1_global",dim_nn(7),j1id)

      !define variables and assign varID
      ns = siid
      ierr=nfmpi_def_var(ncid,"j0",nf_int,1,ns,iij)
      ierr=nfmpi_def_var(ncid,"k0",nf_int,1,ns,iik)
      ierr=nfmpi_def_var(ncid,"tlz_depth",nf_double,1,ns,itlz)
      ierr=nfmpi_def_var(ncid,"dxmnut",nf_double,1,ns,idxm)
      ns = kkid
      ierr=nfmpi_def_var(ncid,"z",nf_double,1,ns,iz)
      ns = k1id
      ierr=nfmpi_def_var(ncid,"odz",nf_double,1,ns,iodz)
      ns = k0id
      ierr=nfmpi_def_var(ncid,"odzw",nf_double,1,ns,iodzw)
      ns = i0id
      ierr=nfmpi_def_var(ncid,"xdeg",nf_double,1,ns,ixdeg)
      ns = j0id
      ierr=nfmpi_def_var(ncid,"y",nf_double,1,ns,iy)
      ierr=nfmpi_def_var(ncid,"yv",nf_double,1,ns,iyv)
      ierr=nfmpi_def_var(ncid,"yvdeg",nf_double,1,ns,iyvd)
      ierr=nfmpi_def_var(ncid,"ydeg",nf_double,1,ns,iydeg)
      ierr=nfmpi_def_var(ncid,"cs",nf_double,1,ns,ics)
      ierr=nfmpi_def_var(ncid,"ocs",nf_double,1,ns,iocs)
      ierr=nfmpi_def_var(ncid,"dx",nf_double,1,ns,idx)
      ierr=nfmpi_def_var(ncid,"odx",nf_double,1,ns,iodx)
      ierr=nfmpi_def_var(ncid,"dy",nf_double,1,ns,idy)
      ierr=nfmpi_def_var(ncid,"ody",nf_double,1,ns,iody)
      ns = j1id
      ierr=nfmpi_def_var(ncid,"csv",nf_double,1,ns,icsv)
      ierr=nfmpi_def_var(ncid,"ocsv",nf_double,1,ns,iocsv)
      ierr=nfmpi_def_var(ncid,"dxv",nf_double,1,ns,idxv)
      ierr=nfmpi_def_var(ncid,"odxv",nf_double,1,ns,iodxv)
      ierr=nfmpi_def_var(ncid,"dyv",nf_double,1,ns,idyv)
      ierr=nfmpi_def_var(ncid,"odyv",nf_double,1,ns,iodyv)
      ierr=nfmpi_enddef(ncid)

      !put variables
      ns_i = 1
      nc_i = 1
      ierr=nfmpi_begin_indep_data(ncid)
      tmp(1) = j0
      ierr=nfmpi_put_var_int(ncid,iij,tmp)
      tmp(1) = k0
      ierr=nfmpi_put_var_int(ncid,iik,tmp)
      ierr=nfmpi_put_var_double(ncid,itlz,tlz_depth)
      ierr=nfmpi_put_var_double(ncid,idxm,dxmnut)
      nc_i = k0+k1
      ierr=nfmpi_put_var_double(ncid,iz,z)
      nc_i = k1
      ierr=nfmpi_put_var_double(ncid,iodz,odz)
      nc_i = k0
      ierr=nfmpi_put_var_double(ncid,iodzw,odzw)
      ierr=nfmpi_end_indep_data(ncid)

      ns_i = mylon*i0+1
      nc_i = i0
      ierr=nfmpi_put_vara_double_all(ncid,ixdeg,ns_i,nc_i,xdeg)
      ns_i = mylat*j0+1
      nc_i = j0
      ierr=nfmpi_put_vara_double_all(ncid,iy,ns_i,nc_i,y)
      ierr=nfmpi_put_vara_double_all(ncid,iyv,ns_i,nc_i,yv)
      ierr=nfmpi_put_vara_double_all(ncid,iyvd,ns_i,nc_i,yvdeg)
      ierr=nfmpi_put_vara_double_all(ncid,iydeg,ns_i,nc_i,ydeg)
      ierr=nfmpi_put_vara_double_all(ncid,ics,ns_i,nc_i,cs)
      ierr=nfmpi_put_vara_double_all(ncid,iocs,ns_i,nc_i,ocs)
      ierr=nfmpi_put_vara_double_all(ncid,idx,ns_i,nc_i,dx)
      ierr=nfmpi_put_vara_double_all(ncid,iodx,ns_i,nc_i,odx)
      ierr=nfmpi_put_vara_double_all(ncid,idy,ns_i,nc_i,dy)
      ierr=nfmpi_put_vara_double_all(ncid,iody,ns_i,nc_i,ody)
      ns_i = mylat*j1+1
      nc_i = j1
      ierr=nfmpi_put_vara_double_all(ncid,icsv,ns_i,nc_i,csv)
      ierr=nfmpi_put_vara_double_all(ncid,iocsv,ns_i,nc_i,ocsv)
      ierr=nfmpi_put_vara_double_all(ncid,idxv,ns_i,nc_i,dxv)
      ierr=nfmpi_put_vara_double_all(ncid,iodxv,ns_i,nc_i,odxv)
      ierr=nfmpi_put_vara_double_all(ncid,idyv,ns_i,nc_i,dyv)
      ierr=nfmpi_put_vara_double_all(ncid,iodyv,ns_i,nc_i,odyv)

      ierr=nfmpi_close(ncid)

end subroutine wr8nc_yzgrid



end module
