program poptrans
  ! Calculate POP zonal mean temperature, from
  ! PCM1/PCTM/CCSM2 POP history tapes
  !
  implicit none
  include 'netcdf.inc'
  !
  ! f90 poptrans.f90 -L/contrib/lib -I/contrib/include -lnetcdf
  ! modified by tc 3/1999 - psi computed from the w method
  ! modified by gs 6/1999 - f90 and output is now netCDF
  ! modified by gs 1/2000 - do each level of U,V,T one at a time to fit
  !                         into memory of T3E
  ! modified by gs 4/2001 - get grid info from pcm1_grid.nc file
  ! Number of single-field format input files per analysis run, i.e.:
  ! 1: b40.20th.track1.1deg.005.pop.h.PD.185001-185912.nc
  ! 2: b40.20th.track1.1deg.005.pop.h.QFLUX.185001-185912.nc
  ! 3: b40.20th.track1.1deg.005.pop.h.SALT.185001-185912.nc
  ! 4: b40.20th.track1.1deg.005.pop.h.SFWF.185001-185912.nc
  ! 5: b40.20th.track1.1deg.005.pop.h.SHF.185001-185912.nc
  ! 6: b40.20th.track1.1deg.005.pop.h.SSH.185001-185912.nc
  ! 7: b40.20th.track1.1deg.005.pop.h.TEMP.185001-185912.nc
  !
  integer,parameter::nphi=360
  integer,parameter::nd =  3,nv = 28,ndv=nd+nv

  real   ,parameter::phi_start=-89.75,phi_end=89.75
  real   ,parameter::spval=1.e30
  real   ,parameter::cal2joule = 4.1868
  !
  integer::varxtype,varndims,varnatts,status,itime,imt,jmt,km,ocnncid
  integer::ocnncid_pd,ocnncid_qflux,ocnncid_salt,ocnncid_sfwf,ocnncid_shf,ocnncid_ssh,ocnncid_temp
  integer::nvars_pd,nvars_qflux,nvars_salt,nvars_sfwf,nvars_shf,nvars_ssh,nvars_temp
  integer::ndims_pd,ndims_qflux,ndims_salt,ndims_sfwf,ndims_shf,ndims_ssh,ndims_temp
  integer::timelen_pd,timelen_qflux,timelen_salt,timelen_sfwf,timelen_shf,timelen_ssh,timelen_temp
  integer::tranncid,ngatts,unlimdimid
  integer::latd,levd,timd,vdim,varid,gridid,nvars,ndims
  integer::i,j,k,l,m,n,jc,ip1,im1,jp1,jm1,iv
  integer,dimension(ndv)::varlev,vrdims
  integer,dimension(4)::istart,icount,vardimids
  integer,dimension(3)::ostart,ocount,vdims
  integer,dimension(20)::nc_dimlen
  integer::nnn
  !
  real::pi,r2d,dphi
  real::wtkb,wtk,fuew,fuww,fvnw,fvsw,work
  real::fuetw,fuwtw,fvntw,fvstw
  integer,dimension(:,:),allocatable::amaskt,imaskt,pmaskt,smaskt,kmu,kmt,region_mask
  real,dimension(:,:),allocatable::ulat,htn,hte,dxu,dyu,tlon,tlat,tarea
  real,dimension(:,:),allocatable::saln,dens,t,shf,sfwf,ssh,qflx
  real,dimension(:,:),allocatable::gt1,at1,it1,pt1,gs1,as1,is1,ps1
  real,dimension(:,:),allocatable::gt2,at2,it2,pt2,gd1,ad1,id1,pd1
  real,dimension(:),allocatable::ght,aht,iht,pht,gsf,asf,isf,psf
  real,dimension(:),allocatable::gsh,ash,ish,psh
  real,dimension(:),allocatable::areag,areaa,areai,areap
  real,dimension(:,:,:),allocatable::w,htw,fue,fvn,fuet,fvnt
  real,dimension(:),allocatable::psia,mhta,psii,mhti,psip,mhtp
  real,dimension(:),allocatable::dz,depth
  real,dimension(nphi)::latr,htrans
  real,dimension(5000)::time_pd,time_qflux,time_salt,time_sfwf,time_shf,time_ssh,time_temp,time
  !
  character(len=132)::varname,title,ncfile,nc_dimname,gridfile
  character(len=132)::infile_pd,infile_qflux,infile_salt,infile_sfwf,infile_shf,infile_ssh,infile_temp
  character(len=132),dimension(ndv)::varnam,vrunit,vrlong,vrtype,vrgrid
  character(len=132),dimension(ndv)::varstr
  !
  logical::does_exist
  !
  varstr = (/&
       !VARNAM     Lv VRLONG                             T  D  Gr VRUNITS
  'lat         1 Latitude                           F  1  XX degrees_north        ',&
       'lev         1 Level                              F  1  XX m                    ',&
       'time        1 Time                               F  1  XX days since 0000-01-01',&
       'gtemp       1 Global zonal mean temperature      F  3  TS degree C             ',&
       'atemp       1 Atlantic zonal mean temperature    F  3  TS degree C             ',&
       'itemp       1 Indian zonal mean temperature      F  3  TS degree C             ',&
       'ptemp       1 Pacific zonal mean temperature     F  3  TS degree C             ',&
       'gsalt       1 Global zonal mean salinity         F  3  TS g/kg                 ',&
       'asalt       1 Atlantic zonal mean salinity       F  3  TS g/kg                 ',&
       'isalt       1 Indian zonal mean salinity         F  3  TS g/kg                 ',&
       'psalt       1 Pacific zonal mean salinity        F  3  TS g/kg                 ',&
       'gdens       1 Global zonal mean density          F  3  TS kg/m^3               ',&
       'adens       1 Atlantic zonal mean density        F  3  TS kg/m^3               ',&
       'idens       1 Indian zonal mean density          F  3  TS kg/m^3               ',&
       'pdens       1 Pacific zonal mean density         F  3  TS kg/m^3               ',&
       'gsht        1 Global zonal mean surf. ht flux    F  2  TS w/m^2                ',&
       'asht        1 Atlantic zonal mean surf. ht flux  F  2  TS w/m^2                ',&
       'psht        1 Pacific zonal mean surf. ht flux   F  2  TS w/m^2                ',&
       'isht        1 Indian zonal mean surf. ht flux    F  2  TS w/m^2                ',&
       'gsfwf       1 Global zonal mean sf frshwt flx    F  2  TS mm/day               ',&
       'asfwf       1 Atlantic zonal mean sf frshwt flx  F  2  TS mm/day               ',&
       'psfwf       1 Pacific zonal mean sf frshwt flx   F  2  TS mm/day               ',&
       'isfwf       1 Indian zonal mean sf frshwt flx    F  2  TS mm/day               ',&
       'gssh        1 Global zonal mean sea surf ht      F  2  TS cm                   ',&
       'assh        1 Atlantic zonal mean sea surf ht    F  2  TS cm                   ',&
       'pssh        1 Pacific zonal mean sea surf ht     F  2  TS cm                   ',&
       'issh        1 Indian zonal mean sea surf ht      F  2  TS cm                   ',&
       'areag       1 global ocean area by latitude      F  2  TS cm^2                 ',&
       'areaa       1 atlantic ocean area by latitude    F  2  TS cm^2                 ',&
       'areai       1 Indian ocean area by latitude      F  2  TS cm^2                 ',&
       'areap       1 Pacific ocean area by latitude     F  2  TS cm^2                 '/)
  !
  write(*,'(''Input PD file: '',$)')
  read (*,'(a)') infile_pd
  write(*,*)
  call check_exists(infile_pd)
  !
  write(*,'(''Input QFLUX file: '',$)')
  read (*,'(a)') infile_qflux
  write(*,*)
  call check_exists(infile_qflux)
  !
  write(*,'(''Input SALT file: '',$)')
  read (*,'(a)') infile_salt
  write(*,*)
  call check_exists(infile_salt)
  !
  write(*,'(''Input SFWF file: '',$)')
  read (*,'(a)') infile_sfwf
  write(*,*)
  call check_exists(infile_sfwf)
  !
  write(*,'(''Input SHF file: '',$)')
  read (*,'(a)') infile_shf
  write(*,*)
  call check_exists(infile_shf)
  !
  write(*,'(''Input SSH file: '',$)')
  read (*,'(a)') infile_ssh
  write(*,*)
  call check_exists(infile_ssh)
  !
  write(*,'(''Input TEMP file: '',$)')
  read (*,'(a)') infile_temp
  write(*,*)
  call check_exists(infile_temp)
  write(*,'(7(a,/))') trim(infile_pd),trim(infile_qflux),trim(infile_salt),trim(infile_sfwf),trim(infile_shf),trim(infile_ssh),trim(infile_temp)
  !
  status = nf_open(infile_pd,nf_nowrite,ocnncid_pd)                        ; if (status /= nf_noerr) write(*,*) nf_strerror(status)
  status = nf_inq(ocnncid_pd,ndims_pd,nvars_pd,ngatts,unlimdimid)          ; if (status /= nf_noerr) write(*,*) nf_strerror(status)
  status = nf_open(infile_qflux,nf_nowrite,ocnncid_qflux)                  ; if (status /= nf_noerr) write(*,*) nf_strerror(status)
  status = nf_inq(ocnncid_qflux,ndims_qflux,nvars_qflux,ngatts,unlimdimid) ; if (status /= nf_noerr) write(*,*) nf_strerror(status)
  status = nf_open(infile_salt,nf_nowrite,ocnncid_salt)                    ; if (status /= nf_noerr) write(*,*) nf_strerror(status)
  status = nf_inq(ocnncid_salt,ndims_salt,nvars_salt,ngatts,unlimdimid)    ; if (status /= nf_noerr) write(*,*) nf_strerror(status)
  status = nf_open(infile_sfwf,nf_nowrite,ocnncid_sfwf)                    ; if (status /= nf_noerr) write(*,*) nf_strerror(status)
  status = nf_inq(ocnncid_sfwf,ndims_sfwf,nvars_sfwf,ngatts,unlimdimid)    ; if (status /= nf_noerr) write(*,*) nf_strerror(status)
  status = nf_open(infile_shf,nf_nowrite,ocnncid_shf)                      ; if (status /= nf_noerr) write(*,*) nf_strerror(status)
  status = nf_inq(ocnncid_shf,ndims_shf,nvars_shf,ngatts,unlimdimid)       ; if (status /= nf_noerr) write(*,*) nf_strerror(status)
  status = nf_open(infile_ssh,nf_nowrite,ocnncid_ssh)                      ; if (status /= nf_noerr) write(*,*) nf_strerror(status)
  status = nf_inq(ocnncid_ssh,ndims_ssh,nvars_ssh,ngatts,unlimdimid)       ; if (status /= nf_noerr) write(*,*) nf_strerror(status)
  status = nf_open(infile_temp,nf_nowrite,ocnncid_temp)                    ; if (status /= nf_noerr) write(*,*) nf_strerror(status)
  status = nf_inq(ocnncid_temp,ndims_temp,nvars_temp,ngatts,unlimdimid)    ; if (status /= nf_noerr) write(*,*) nf_strerror(status)

  ! Find out I and J and K dimensions
  do i = 1,ndims_pd
     status = nf_inq_dim(ocnncid_pd,i,nc_dimname,nc_dimlen(i))
     if (status /= nf_noerr) write(*,*) nf_strerror(status)
     if (trim(nc_dimname) == 'lon')  imt = nc_dimlen(i)
     if (trim(nc_dimname) == 'lat')  jmt = nc_dimlen(i)
     if (trim(nc_dimname) == 'lev')   km = nc_dimlen(i)
     if (trim(nc_dimname) == 'nlon') imt = nc_dimlen(i)
     if (trim(nc_dimname) == 'nlat') jmt = nc_dimlen(i)
     if (trim(nc_dimname) == 'z_t')   km = nc_dimlen(i)
  enddo
  !
  allocate(amaskt(imt,jmt),imaskt(imt,jmt),pmaskt(imt,jmt),smaskt(imt,jmt))
  allocate(kmu(imt,jmt),kmt(imt,jmt),region_mask(imt,jmt))
  allocate(ulat(imt,jmt),htn(imt,jmt),hte(imt,jmt))
  allocate(dxu(imt,jmt),dyu(imt,jmt),tlon(imt,jmt),tlat(imt,jmt),tarea(imt,jmt))
  allocate(saln(imt,jmt),dens(imt,jmt),t(imt,jmt))
  allocate(shf(imt,jmt),sfwf(imt,jmt),ssh(imt,jmt),qflx(imt,jmt))
  allocate(gt1(km,nphi),at1(km,nphi),gt2(km,nphi),at2(km,nphi))
  allocate(it1(km,nphi),pt1(km,nphi),it2(km,nphi),pt2(km,nphi))
  allocate(gd1(km,nphi),ad1(km,nphi),gs1(km,nphi),as1(km,nphi))
  allocate(id1(km,nphi),pd1(km,nphi),is1(km,nphi),ps1(km,nphi))
  allocate(ght(nphi),aht(nphi),pht(nphi),iht(nphi))
  allocate(gsf(nphi),asf(nphi),psf(nphi),isf(nphi))
  allocate(gsh(nphi),ash(nphi),psh(nphi),ish(nphi))
  allocate(areag(nphi),areaa(nphi),areai(nphi),areap(nphi))
  allocate(w(imt,jmt,km),htw(imt,jmt,km),fue(imt,jmt,km),fvn(imt,jmt,km))
  allocate(fuet(imt,jmt,km),fvnt(imt,jmt,km))
  allocate(psia(km),mhta(km),psii(km),mhti(km),psip(km),mhtp(km))
  allocate(dz(km),depth(km))
  gt1 = 0.;at1 = 0.;it1 = 0.;pt1 = 0.
  gt2 = 0.;at2 = 0.;it2 = 0.;pt2 = 0.
  gs1 = 0.;as1 = 0.;is1 = 0.;ps1 = 0.
  gd1 = 0.;ad1 = 0.;id1 = 0.;pd1 = 0.
  ght = 0.;aht = 0.;iht = 0.;pht = 0.
  gsf = 0.;asf = 0.;isf = 0.;psf = 0.
  gsh = 0.;ash = 0.;ish = 0.;psh = 0.

  ! Find out how many time samples there are in each file
  do i = 1,nvars_pd
     status = nf_inq_var(ocnncid_pd,i,varname,varxtype,varndims,vardimids,varnatts)
     if (status /= nf_noerr) write(*,*) nf_strerror(status)
     if ((trim(varname) == 'time').or.(trim(varname) ==  'TIME')) then
        status = nf_inq_dim(ocnncid_pd,unlimdimid,varname,timelen_pd) ; if (status /= nf_noerr) write(*,*) nf_strerror(status)
        status = nf_get_var_real(ocnncid_pd,i,time_pd)                ; if (status /= nf_noerr) write(*,*) nf_strerror(status)
     endif
  enddo
  do i = 1,nvars_qflux
     status = nf_inq_var(ocnncid_qflux,i,varname,varxtype,varndims,vardimids,varnatts)
     if (status /= nf_noerr) write(*,*) nf_strerror(status)
     if ((trim(varname) == 'time').or.(trim(varname) ==  'TIME')) then
        status = nf_inq_dim(ocnncid_qflux,unlimdimid,varname,timelen_qflux) ; if (status /= nf_noerr) write(*,*) nf_strerror(status)
        status = nf_get_var_real(ocnncid_qflux,i,time_qflux)                ; if (status /= nf_noerr) write(*,*) nf_strerror(status)
     endif
  enddo
  do i = 1,nvars_salt
     status = nf_inq_var(ocnncid_salt,i,varname,varxtype,varndims,vardimids,varnatts)
     if (status /= nf_noerr) write(*,*) nf_strerror(status)
     if ((trim(varname) == 'time').or.(trim(varname) ==  'TIME')) then
        status = nf_inq_dim(ocnncid_salt,unlimdimid,varname,timelen_salt) ; if (status /= nf_noerr) write(*,*) nf_strerror(status)
        status = nf_get_var_real(ocnncid_salt,i,time_salt)                ; if (status /= nf_noerr) write(*,*) nf_strerror(status)
     endif
  enddo
  do i = 1,nvars_sfwf
     status = nf_inq_var(ocnncid_sfwf,i,varname,varxtype,varndims,vardimids,varnatts)
     if (status /= nf_noerr) write(*,*) nf_strerror(status)
     if ((trim(varname) == 'time').or.(trim(varname) ==  'TIME')) then
        status = nf_inq_dim(ocnncid_sfwf,unlimdimid,varname,timelen_sfwf) ; if (status /= nf_noerr) write(*,*) nf_strerror(status)
        status = nf_get_var_real(ocnncid_sfwf,i,time_sfwf)                ; if (status /= nf_noerr) write(*,*) nf_strerror(status)
     endif
  enddo
  do i = 1,nvars_shf
     status = nf_inq_var(ocnncid_shf,i,varname,varxtype,varndims,vardimids,varnatts)
     if (status /= nf_noerr) write(*,*) nf_strerror(status)
     if ((trim(varname) == 'time').or.(trim(varname) ==  'TIME')) then
        status = nf_inq_dim(ocnncid_shf,unlimdimid,varname,timelen_shf) ; if (status /= nf_noerr) write(*,*) nf_strerror(status)
        status = nf_get_var_real(ocnncid_shf,i,time_shf)                ; if (status /= nf_noerr) write(*,*) nf_strerror(status)
     endif
  enddo
  do i = 1,nvars_ssh
     status = nf_inq_var(ocnncid_ssh,i,varname,varxtype,varndims,vardimids,varnatts)
     if (status /= nf_noerr) write(*,*) nf_strerror(status)
     if ((trim(varname) == 'time').or.(trim(varname) ==  'TIME')) then
        status = nf_inq_dim(ocnncid_ssh,unlimdimid,varname,timelen_ssh) ; if (status /= nf_noerr) write(*,*) nf_strerror(status)
        status = nf_get_var_real(ocnncid_ssh,i,time_ssh)                ; if (status /= nf_noerr) write(*,*) nf_strerror(status)
     endif
  enddo
  do i = 1,nvars_temp
     status = nf_inq_var(ocnncid_temp,i,varname,varxtype,varndims,vardimids,varnatts)
     if (status /= nf_noerr) write(*,*) nf_strerror(status)
     if ((trim(varname) == 'time').or.(trim(varname) ==  'TIME')) then
        status = nf_inq_dim(ocnncid_temp,unlimdimid,varname,timelen_temp) ; if (status /= nf_noerr) write(*,*) nf_strerror(status)
        status = nf_get_var_real(ocnncid_temp,i,time_temp)                ; if (status /= nf_noerr) write(*,*) nf_strerror(status)
     endif
  enddo
  write(*,'(7i6)') timelen_pd,timelen_qflux,timelen_salt,timelen_sfwf,timelen_shf,timelen_ssh,timelen_temp
  !
  istart =(/   1,  1,  1,  1/)
  icount =(/ imt,jmt,  1,  1/)
  ostart =(/   1,  1,  1/)
  ocount =(/nphi, km,  1/)
  !
  pi=atan(1.)*4.
  r2d=180./pi
  !
  if (jmt == 288) gridfile = 'pcm1_grid.nc'
  if (jmt == 320) gridfile = 'pcm2_grid.nc'
  if (jmt == 384) gridfile = 'ccsm2_ocn.nc'
  if (km.eq.60) then
     gridfile = infile_pd
  else
     inquire(file=trim(gridfile),exist=does_exist)
     if (.not.(does_exist)) then
        write(*,*) 'Need to link in ',trim(gridfile),' as follows: '
        write(*,*) 'ln -s /fs/cgd/data0/strandwg/ferret/data/',trim(gridfile),' .'
        write(*,*) 'Stopping.'
        stop
     endif
  endif
  !
  ! Get grid info
  !
  status = nf_open(trim(gridfile),nf_nowrite,gridid)    ; if (status /= nf_noerr) write(*,*) nf_strerror(status)
  status = nf_inq(gridid,ndims,nvars,ngatts,unlimdimid) ; if (status /= nf_noerr) write(*,*) nf_strerror(status)
  do i = 1,nvars
     status = nf_inq_var(gridid,i,varname,varxtype,varndims,vardimids,varnatts)
     if (status /= nf_noerr) write(*,*) nf_strerror(status)
     select case (trim(varname))
     case ('kmt','KMT')
        status = nf_get_var_int(gridid,i,kmt)
        if (status /= nf_noerr) write(*,*) nf_strerror(status)
     case ('amaskm')
        status = nf_get_var_int(gridid,i,amaskt)
        if (status /= nf_noerr) write(*,*) nf_strerror(status)
     case ('imaskt')
        status = nf_get_var_int(gridid,i,imaskt)
        if (status /= nf_noerr) write(*,*) nf_strerror(status)
     case ('pmaskt')
        status = nf_get_var_int(gridid,i,pmaskt)
        if (status /= nf_noerr) write(*,*) nf_strerror(status)
     case ('smaskt')
        status = nf_get_var_int(gridid,i,smaskt)
        if (status /= nf_noerr) write(*,*) nf_strerror(status)
     case ('REGION_MASK')
        status = nf_get_var_int(gridid,i,region_mask)
        if (status /= nf_noerr) write(*,*) nf_strerror(status)
        amaskt = 0 ; imaskt = 0 ; pmaskt = 0 ; smaskt = 0
        where (region_mask == 2) ; pmaskt = 1 ; endwhere
        where (region_mask == 3) ; imaskt = 1 ; endwhere
        where (region_mask == 6) ; amaskt = 1 ; endwhere
        where (region_mask == 8) ; amaskt = 1 ; endwhere
	where (region_mask == 9) ; amaskt = 1 ; endwhere
        where (region_mask == 1) ; smaskt = 1 ; endwhere
     case ('lat','ULAT')
        status = nf_get_var_real(gridid,i,ulat)
        if (status /= nf_noerr) write(*,*) nf_strerror(status)
        ulat = ulat / r2d
     case ('lon','TLON','tlong','TLONG')
        status = nf_get_var_real(gridid,i,tlon)
        if (status /= nf_noerr) write(*,*) nf_strerror(status)
        !        tlon = tlon / r2d
     case ('htn','HTN')
        status = nf_get_var_real(gridid,i,htn)
        if (status /= nf_noerr) write(*,*) nf_strerror(status)
     case ('hte','HTE')
        status = nf_get_var_real(gridid,i,hte)
        if (status /= nf_noerr) write(*,*) nf_strerror(status)
     case ('dz')
        status = nf_get_var_real(gridid,i,dz)
        if (status /= nf_noerr) write(*,*) nf_strerror(status)
     case default
     end select
  enddo
  status = nf_close(gridid)
  if (status /= nf_noerr) write(*,*) nf_strerror(status)
  !
  do j=1,jmt
     do i=1,imt
        im1=mod(i-2+imt,imt)+1
        ip1=mod(i,imt)+1
        jm1=mod(j-2+jmt,jmt)+1
        jp1=mod(j,jmt)+1
        kmu(i,j)=min(kmt(i,j),kmt(ip1,j),kmt(i,jp1),kmt(ip1,jp1))
        tlat(i,j)=0.25*(ulat(i,j)+ulat(im1,j)+ulat(i,jm1)+ulat(im1,jm1))*r2d
        tarea(i,j)=0.25*(htn(i,j)+htn(i,jm1))*(hte(i,j)+hte(im1,j))
        dxu(i,j)=0.5*(htn(i,j)+htn(ip1,j))
        dyu(i,j)=0.5*(hte(i,j)+hte(i,jp1))
     enddo
  enddo

  depth(1)=dz(1)
  do k=2,km
     depth(k)=(depth(k-1)+dz(k))
  enddo
  depth = depth/100 ! cm to m

  dphi = (phi_end - phi_start)/float(nphi - 1)
  do n = 1,nphi
     latr(n) = (phi_start+float(n-1)*dphi)
  enddo
  write(*,*) 'Grid information read and computed successfully.'

  ! Create new netCDF file for transports
  ncfile = 'TSDsurf_'//trim(infile_pd)
  do i = 1,ndv
     read(varstr(i),'(a10,i3,1x,a34,1x,a1,2x,i1,2x,a2,1x,a)') &
          varnam(i),varlev(i),vrlong(i),vrtype(i),vrdims(i),vrgrid(i),vrunit(i)
  enddo
  status = nf_create(trim(ncfile),nf_clobber,tranncid)
  !
  status = nf_def_dim(tranncid,varnam(1),nphi,latd)
  status = nf_def_dim(tranncid,varnam(2),km,levd)
  status = nf_def_dim(tranncid,varnam(3),NF_UNLIMITED,timd)
  !
  do i = 1,nd
     if (i == 1) vdim = latd
     if (i == 2) vdim = levd
     if (i == 3) vdim = timd
     status = nf_def_var(tranncid,varnam(i),ncfloat,vrdims(i),vdim,varid)
     status = nf_put_att_text(tranncid,varid,'units',len_trim(vrunit(i)),vrunit(i))
     status = nf_put_att_text(tranncid,varid,'long_name',len_trim(vrlong(i)),vrlong(i))
     if (i == 2) status = nf_put_att_text(tranncid,varid,'positive',4,'down')
  enddo
  !
  do i = 1,nv
     j = (i+nd)
     vdims(1) = latd
     if (vrdims(j) == 3) then
        vdims(2) = levd
        vdims(3) = timd
     else
        vdims(2) = timd
     endif
     status = nf_def_var(tranncid,varnam(j),ncfloat,vrdims(j),vdims,varid)
     status = nf_put_att_text(tranncid,varid,'units',len_trim(vrunit(j)),vrunit(j))
     status = nf_put_att_text(tranncid,varid,'long_name',len_trim(vrlong(j)),vrlong(j))
     status = nf_put_att_real(tranncid,varid,'missing_value',ncfloat,1,spval)
     status = nf_put_att_real(tranncid,varid,'_FillValue',ncfloat,1,spval)
  enddo
  !
  title  = 'zonal mean fields from '//trim(infile_pd)
  status = nf_put_att_text(tranncid,nf_global,'title',len_trim(title),title)
  status = nf_enddef(tranncid)
  !
  status = nf_inq_varid(tranncid,trim(varnam(1)),varid)
  status = nf_put_var_real(tranncid,varid,latr)
  if (status /= nf_noerr) write(*,*) nf_strerror(status)
  status = nf_inq_varid(tranncid,trim(varnam(2)),varid)
  status = nf_put_var_real(tranncid,varid,depth)
  if (status /= nf_noerr) write(*,*) nf_strerror(status)
  !
  ! Now do the calculation
  !
  do itime = 1,timelen_pd
     do k=1,km
        ! get T,S
        saln = spval ; dens = spval ; t = spval
        ssh  = spval ; sfwf = spval ; shf = spval ; qflx = spval
        istart(3) = k ; istart(4) = itime
        !
        do iv = 1,nvars_pd
           status = nf_inq_var(ocnncid_pd,iv,varname,varxtype,varndims,vardimids,varnatts) ; if (status /= nf_noerr) write(*,*) nf_strerror(status)
           select case(trim(varname))
           case ('pd','PD')
              status = nf_get_vara_real(ocnncid_pd,iv,istart,icount,dens) ; if (status /= nf_noerr) write(*,*) nf_strerror(status)
           case default
           end select
        enddo
        do iv = 1,nvars_temp
           status = nf_inq_var(ocnncid_temp,iv,varname,varxtype,varndims,vardimids,varnatts) ; if (status /= nf_noerr) write(*,*) nf_strerror(status)
           select case(trim(varname))
           case ('temp','TEMP')
              status = nf_get_vara_real(ocnncid_temp,iv,istart,icount,t)
              if (status /= nf_noerr) write(*,*) nf_strerror(status)
           case default
           end select
        enddo
        do iv = 1,nvars_salt
           status = nf_inq_var(ocnncid_salt,iv,varname,varxtype,varndims,vardimids,varnatts) ; if (status /= nf_noerr) write(*,*) nf_strerror(status)
           select case(trim(varname))
           case ('salt','SALT','S','s')
              status = nf_get_vara_real(ocnncid_salt,iv,istart,icount,saln)
              if (status /= nf_noerr) write(*,*) nf_strerror(status)
           case default
           end select
        enddo
        if (k==1) then
           istart(3) = itime
           do iv = 1,nvars_shf
              status = nf_inq_var(ocnncid_shf,iv,varname,varxtype,varndims,vardimids,varnatts) ; if (status /= nf_noerr) write(*,*) nf_strerror(status)
              select case(trim(varname))
	      case ('shf','SHF')
	         status = nf_get_vara_real(ocnncid_shf,iv,istart,icount,shf)
                 if (status /= nf_noerr) write(*,*) nf_strerror(status)
              case default
              end select
           enddo
           do iv = 1,nvars_ssh
              status = nf_inq_var(ocnncid_ssh,iv,varname,varxtype,varndims,vardimids,varnatts) ; if (status /= nf_noerr) write(*,*) nf_strerror(status)
              select case(trim(varname))
	      case ('ah1','SSH')
		 status = nf_get_vara_real(ocnncid_ssh,iv,istart,icount,ssh)
                 if (status /= nf_noerr) write(*,*) nf_strerror(status)
              case default
              end select
           enddo
           do iv = 1,nvars_sfwf
              status = nf_inq_var(ocnncid_sfwf,iv,varname,varxtype,varndims,vardimids,varnatts) ; if (status /= nf_noerr) write(*,*) nf_strerror(status)
              select case(trim(varname))
	      case ('ssf','SFWF')
		 status = nf_get_vara_real(ocnncid_sfwf,iv,istart,icount,sfwf)
                 if (status /= nf_noerr) write(*,*) nf_strerror(status)
              case default
              end select
           enddo
           do iv = 1,nvars_qflux
              status = nf_inq_var(ocnncid_qflux,iv,varname,varxtype,varndims,vardimids,varnatts) ; if (status /= nf_noerr) write(*,*) nf_strerror(status)
              select case(trim(varname))
	      case ('QFLUX')
		 status = nf_get_vara_real(ocnncid_qflux,iv,istart,icount,qflx)
                 if (status /= nf_noerr) write(*,*) nf_strerror(status)
              case default
              end select
           enddo
        end if
        if (jmt == 288) then
           saln=saln*1000.
           !	print*,'before density',k,t(25,150),saln(25,150),dens(25,150)
           call density(imt,jmt,t,saln,dens,kmt,k,0.)
        end if
        !        if (jmt == 320 .and. k == 1 ) then
        !       qflx=qflx/3.34e5
        !       sfwf=(sfwf-qflx)*8.64*1.e4   !convert to mm/day
        !       end if

        !        if (k.eq.1) print*,'sfwf(200,200)=',sfwf(200,200),'shf(200,200)=',shf(200,200)

        !        print*,'after density',k,t(25,150),saln(25,150),dens(25,150)
        ! Clean them up a wee bit
        where (abs(dens)>1.e10) ; dens = 0. ; end where
        where (abs(saln)>1.e10) ; saln = 0. ; end where
        where (abs(t)>1.e10) ; t = 0. ; end where
        if (k == 1) then
           where (abs(ssh)>1.e10) ; ssh = 0. ; end where
           where (abs(shf)>1.e10) ; shf = 0. ; end where
           where (abs(qflx)>1.e10) ; shf = 0. ; end where
           where (abs(sfwf)>1.e10) ; sfwf = 0. ; end where
           where (kmt<k) ; ssh(:,:) = 0.; shf(:,:) = 0.  ; end where
           where (kmt<k) ; sfwf(:,:) = 0.; end where
        end if

        where (kmt<k) ; dens(:,:) = 0. ; saln(:,:) = 0. ; end where
        where (kmt<k) ; t(:,:) = 0. ; end where

        if (jmt == 320 .and. k == 1 ) then
           qflx=qflx/3.34e5
           sfwf=(sfwf-qflx)*8.64*1.e4   !convert to mm/day
        end if
        !        if (k.eq.1) print*,'sfwf(200,200)=',sfwf(200,200),'shf(200,200)=',shf(200,200)
        ! precompute flux information
        ! compute dxu,dyu,fue,fvn like wcalc in pop
        ! fue, fvn is the flux through the east and north edge, u*dy or v*dx
        do j=1,jmt
           do i=1,imt
              im1=mod(i-2+imt,imt)+1
              ip1=mod(i,imt)+1
              jm1=mod(j-2+jmt,jmt)+1
              jp1=mod(j,jmt)+1
              !
              fuet(i,j,k)=t(i,j)*tarea(i,j)    !temperature
              fue(i,j,k)=saln(i,j)*tarea(i,j)     !salinity
     	      fvn(i,j,k)=dens(i,j)*tarea(i,j)     !density
           enddo
        enddo
        if (k == 1) then
           do j=1,jmt
              do i=1,imt
                 fvnt(i,j,1)=shf(i,j)*tarea(i,j)
                 fvnt(i,j,2)=sfwf(i,j)*tarea(i,j)
                 fvnt(i,j,3)=ssh(i,j)*tarea(i,j)
              enddo
           enddo
        endif
     enddo ! do k = 1,km
     print*,"test case"
     !
     gt1=0.;gt2=0.;gs1=0.;gd1=0.
     at1=0.;at2=0.;as1=0.;ad1=0.
     it1=0.;it2=0.;is1=0.;id1=0.
     pt1=0.;pt2=0.;ps1=0.;pd1=0.

     !
     do n = 2,nphi
        !
  	nnn=0
        do j=1,jmt
           do i=1,imt
              if (tlat(i,j)<latr(n).and.tlat(i,j)>=latr(n-1)) then
                 nnn=nnn+1
                 do k=1,km
                    if (k<=kmt(i,j)) then
                       gt1(k,n)=gt1(k,n)+fuet(i,j,k)
                       gt2(k,n)=gt2(k,n)+tarea(i,j)
                       gs1(k,n)=gs1(k,n)+fue(i,j,k)
                       gd1(k,n)=gd1(k,n)+fvn(i,j,k)
		       if (k == 1) then
                          ght(n)=ght(n)+fvnt(i,j,1)
                          gsf(n)=gsf(n)+fvnt(i,j,2)
                          gsh(n)=gsh(n)+fvnt(i,j,3)
	               end if
                       if (amaskt(i,j) == 1 .or. (smaskt(i,j)==1.and.(tlon(i,j)>290. .or. tlon(i,j) <= 20.))) then
                          at1(k,n)=at1(k,n)+fuet(i,j,k)
                          at2(k,n)=at2(k,n)+tarea(i,j)
                          as1(k,n)=as1(k,n)+fue(i,j,k)
                          ad1(k,n)=ad1(k,n)+fvn(i,j,k)
                          if (k == 1) then
                             aht(n)=aht(n)+fvnt(i,j,1)
                             asf(n)=asf(n)+fvnt(i,j,2)
                             ash(n)=ash(n)+fvnt(i,j,3)
                          end if
                       endif
                       if (imaskt(i,j) == 1.or.(smaskt(i,j)==1.and.tlon(i,j)> 20. .and.tlon(i,j) <= 145.)) then
                          it1(k,n)=it1(k,n)+fuet(i,j,k)
                          it2(k,n)=it2(k,n)+tarea(i,j)
                          is1(k,n)=is1(k,n)+fue(i,j,k)
                          id1(k,n)=id1(k,n)+fvn(i,j,k)
                          if (k == 1) then
                             iht(n)=iht(n)+fvnt(i,j,1)
                             isf(n)=isf(n)+fvnt(i,j,2)
                             ish(n)=ish(n)+fvnt(i,j,3)
                          end if
                       endif
                       if (pmaskt(i,j) == 1 .or. (smaskt(i,j)==1.and.(tlon(i,j) >145. .and. tlon(i,j) <= 290.))) then
                          pt1(k,n)=pt1(k,n)+fuet(i,j,k)
                          pt2(k,n)=pt2(k,n)+tarea(i,j)
                          ps1(k,n)=ps1(k,n)+fue(i,j,k)
                          pd1(k,n)=pd1(k,n)+fvn(i,j,k)
                          if (k == 1) then
                             pht(n)=pht(n)+fvnt(i,j,1)
                             psf(n)=psf(n)+fvnt(i,j,2)
                             psh(n)=psh(n)+fvnt(i,j,3)
                          end if
                       endif

                    endif
                 enddo ! do k=1,km
              endif

           enddo ! do i=1,imt
        enddo ! do j=1,jmt

        if (n.eq.200) print*,'asf(200)=',asf(200)
        if (gt2(1,n).gt.0.) then
           ght(n)=ght(n)/gt2(1,n)
           gsf(n)=gsf(n)/gt2(1,n)
           gsh(n)=gsh(n)/gt2(1,n)
        else
           ght(n)=1.e30
           gsf(n)=1.e30
           gsh(n)=1.e30
        end if
        if (at2(1,n).gt.0.) then
           aht(n)=aht(n)/at2(1,n)
           asf(n)=asf(n)/at2(1,n)
           ash(n)=ash(n)/at2(1,n)
        else
           aht(n)=1.e30
           asf(n)=1.e30
           ash(n)=1.e30
        end if
        if (it2(1,n).gt.0.) then
           iht(n)=iht(n)/it2(1,n)
           isf(n)=isf(n)/it2(1,n)
           ish(n)=ish(n)/it2(1,n)
        else
           iht(n)=1.e30
           isf(n)=1.e30
           ish(n)=1.e30
        end if
        if (pt2(1,n).gt.0.) then
           pht(n)=pht(n)/pt2(1,n)
           psf(n)=psf(n)/pt2(1,n)
           psh(n)=psh(n)/pt2(1,n)
        else
           pht(n)=1.e30
           psf(n)=1.e30
           psh(n)=1.e30
        end if
        areag=gt2(1,:)
        areaa=at2(1,:)
        areai=it2(1,:)
        areap=pt2(1,:)

        do k=1,km
           if (gt2(k,n).gt.0.) then
              gt1(k,n)=gt1(k,n)/gt2(k,n)
              gs1(k,n)=gs1(k,n)/gt2(k,n)
              gd1(k,n)=gd1(k,n)/gt2(k,n)
           else
              gt1(k,n)=1.e30
              gs1(k,n)=1.e30
              gd1(k,n)=1.e30
           end if
           if (at2(k,n).gt.0.) then
              at1(k,n)=at1(k,n)/at2(k,n)
              as1(k,n)=as1(k,n)/at2(k,n)
              ad1(k,n)=ad1(k,n)/at2(k,n)
           else
              at1(k,n)=1.e30
              as1(k,n)=1.e30
              ad1(k,n)=1.e30
           end if
           if (it2(k,n).gt.0.) then
              it1(k,n)=it1(k,n)/it2(k,n)
              is1(k,n)=is1(k,n)/it2(k,n)
              id1(k,n)=id1(k,n)/it2(k,n)
           else
              it1(k,n)=1.e30
              is1(k,n)=1.e30
              id1(k,n)=1.e30
           end if
           if (pt2(k,n).gt.0.) then
              pt1(k,n)=pt1(k,n)/pt2(k,n)
              ps1(k,n)=ps1(k,n)/pt2(k,n)
              pd1(k,n)=pd1(k,n)/pt2(k,n)
           else
              pt1(k,n)=1.e30
              ps1(k,n)=1.e30
              pd1(k,n)=1.e30
           end if

	end do
        if (nnn.le.10 .and.latr(n).gt.-75.5) then
           gt1(:,n)=gt1(:,n-1)
           gs1(:,n)=gs1(:,n-1)
           gd1(:,n)=gd1(:,n-1)
           at1(:,n)=at1(:,n-1)
           as1(:,n)=as1(:,n-1)
           ad1(:,n)=ad1(:,n-1)
           it1(:,n)=it1(:,n-1)
           is1(:,n)=is1(:,n-1)
           id1(:,n)=id1(:,n-1)
           pt1(:,n)=pt1(:,n-1)
           ps1(:,n)=ps1(:,n-1)
           pd1(:,n)=pd1(:,n-1)
           ght(n)=ght(n-1)
           aht(n)=aht(n-1)
           iht(n)=iht(n-1)
           pht(n)=pht(n-1)
           gsh(n)=gsh(n-1)
           ash(n)=ash(n-1)
           ish(n)=ish(n-1)
           psh(n)=psh(n-1)
           gsf(n)=gsf(n-1)
           asf(n)=asf(n-1)
           isf(n)=isf(n-1)
           psf(n)=psf(n-1)
           areag(n)=areag(n-1)
           areaa(n)=areaa(n-1)
           areai(n)=areai(n-1)
           areap(n)=areap(n-1)
        end if
        if (latr(n).gt.-17..and.latr(n).lt.-15.5) then
           gt1(:,n)=gt1(:,n-1)
           gs1(:,n)=gs1(:,n-1)
           gd1(:,n)=gd1(:,n-1)
           at1(:,n)=at1(:,n-1)
           as1(:,n)=as1(:,n-1)
           ad1(:,n)=ad1(:,n-1)
           it1(:,n)=it1(:,n-1)
           is1(:,n)=is1(:,n-1)
           id1(:,n)=id1(:,n-1)
           pt1(:,n)=pt1(:,n-1)
           ps1(:,n)=ps1(:,n-1)
           pd1(:,n)=pd1(:,n-1)
           ght(n)=ght(n-1)
           aht(n)=aht(n-1)
           iht(n)=iht(n-1)
           pht(n)=pht(n-1)
           gsh(n)=gsh(n-1)
           ash(n)=ash(n-1)
           ish(n)=ish(n-1)
           psh(n)=psh(n-1)
           gsf(n)=gsf(n-1)
           asf(n)=asf(n-1)
           isf(n)=isf(n-1)
           psf(n)=psf(n-1)
           areag(n)=areag(n-1)
           areaa(n)=areaa(n-1)
           areai(n)=areai(n-1)
           areap(n)=areap(n-1)
	end if
        if (jmt.eq.288.and.it1(1,n).gt.1.e10.and.latr(n).gt.2.) then
           it1(:,n)=it1(:,n-1)
           is1(:,n)=is1(:,n-1)
           id1(:,n)=id1(:,n-1)
           iht(n)=iht(n-1)
           ish(n)=ish(n-1)
           isf(n)=isf(n-1)
           areai(n)=areai(n-1)
        end if
     enddo ! do n=1,nphi

     print*,'asf(200)=',asf(200)
     ! Add variable data to netCDF file
     status = nf_inq_varid(tranncid,trim(varnam(3)),varid)
     if (status /= nf_noerr) write(*,*) nf_strerror(status)
     status = nf_put_vara_real(tranncid,varid,itime,1,time_pd(itime))
     if (status /= nf_noerr) write(*,*) nf_strerror(status)
     do l = 1,nv
        m = l+nd
        ostart(1) = 1 ; ocount(1) = nphi
        status = nf_inq_varid(tranncid,varnam(m),varid)
        if (vrdims(m) == 3) then
           ostart(2) = 1     ; ocount(2) =   km
           ostart(3) = itime ; ocount(3) =    1
        else
           ostart(2) = itime ; ocount(2) =    1
        endif
        !        if (vrdims(m) == 3) then
        !           ostart(1) = 1 ; ocount(1) = nphi
        !           ostart(2) = 1     ; ocount(2) =   km
        !           ostart(3) = itime ; ocount(3) =    1
        !        elseif (vrdims(m) == 2) then
        !           ostart(1) = 1 ; ocount(1) = nphi
        !           ostart(2) = itime ; ocount(2) =    1
        !        else
        !           ostart(1) = itime ; ocount(1) = 1
        !        endif

        select case (trim(varnam(m)))
        case ('gtemp')
           status = nf_put_vara_real(tranncid,varid,ostart,ocount,transpose(gt1))
           if (status /= nf_noerr) write(*,*) nf_strerror(status)
        case ('atemp')
           status = nf_put_vara_real(tranncid,varid,ostart,ocount,transpose(at1))
           if (status /= nf_noerr) write(*,*) nf_strerror(status)
        case ('itemp')
           status = nf_put_vara_real(tranncid,varid,ostart,ocount,transpose(it1))
           if (status /= nf_noerr) write(*,*) nf_strerror(status)
        case ('ptemp')
           status = nf_put_vara_real(tranncid,varid,ostart,ocount,transpose(pt1))
           if (status /= nf_noerr) write(*,*) nf_strerror(status)
        case ('gsalt')
           status = nf_put_vara_real(tranncid,varid,ostart,ocount,transpose(gs1))
           if (status /= nf_noerr) write(*,*) nf_strerror(status)
        case ('asalt')
           status = nf_put_vara_real(tranncid,varid,ostart,ocount,transpose(as1))
           if (status /= nf_noerr) write(*,*) nf_strerror(status)
        case ('isalt')
           status = nf_put_vara_real(tranncid,varid,ostart,ocount,transpose(is1))
           if (status /= nf_noerr) write(*,*) nf_strerror(status)
        case ('psalt')
           status = nf_put_vara_real(tranncid,varid,ostart,ocount,transpose(ps1))
           if (status /= nf_noerr) write(*,*) nf_strerror(status)
        case ('gdens')
           status = nf_put_vara_real(tranncid,varid,ostart,ocount,transpose(gd1))
           if (status /= nf_noerr) write(*,*) nf_strerror(status)
        case ('adens')
           status = nf_put_vara_real(tranncid,varid,ostart,ocount,transpose(ad1))
           if (status /= nf_noerr) write(*,*) nf_strerror(status)
        case ('idens')
           status = nf_put_vara_real(tranncid,varid,ostart,ocount,transpose(id1))
           if (status /= nf_noerr) write(*,*) nf_strerror(status)
        case ('pdens')
           status = nf_put_vara_real(tranncid,varid,ostart,ocount,transpose(pd1))
           if (status /= nf_noerr) write(*,*) nf_strerror(status)
        case ('gsht')
           status = nf_put_vara_real(tranncid,varid,ostart,ocount,ght)
           if (status /= nf_noerr) write(*,*) nf_strerror(status)
        case ('asht')
           status = nf_put_vara_real(tranncid,varid,ostart,ocount,aht)
           if (status /= nf_noerr) write(*,*) nf_strerror(status)
        case ('isht')
           status = nf_put_vara_real(tranncid,varid,ostart,ocount,iht)
           if (status /= nf_noerr) write(*,*) nf_strerror(status)
        case ('psht')
           status = nf_put_vara_real(tranncid,varid,ostart,ocount,pht)
           if (status /= nf_noerr) write(*,*) nf_strerror(status)
        case ('gsfwf')
           status = nf_put_vara_real(tranncid,varid,ostart,ocount,gsf)
           if (status /= nf_noerr) write(*,*) nf_strerror(status)
        case ('asfwf')
           status = nf_put_vara_real(tranncid,varid,ostart,ocount,asf)
           if (status /= nf_noerr) write(*,*) nf_strerror(status)
        case ('isfwf')
           status = nf_put_vara_real(tranncid,varid,ostart,ocount,isf)
           if (status /= nf_noerr) write(*,*) nf_strerror(status)
        case ('psfwf')
           status = nf_put_vara_real(tranncid,varid,ostart,ocount,psf)
           if (status /= nf_noerr) write(*,*) nf_strerror(status)
        case ('gssh')
           status = nf_put_vara_real(tranncid,varid,ostart,ocount,gsh)
           if (status /= nf_noerr) write(*,*) nf_strerror(status)
        case ('assh')
           status = nf_put_vara_real(tranncid,varid,ostart,ocount,ash)
           if (status /= nf_noerr) write(*,*) nf_strerror(status)
        case ('issh')
           status = nf_put_vara_real(tranncid,varid,ostart,ocount,ish)
           if (status /= nf_noerr) write(*,*) nf_strerror(status)
        case ('pssh')
           status = nf_put_vara_real(tranncid,varid,ostart,ocount,psh)
           if (status /= nf_noerr) write(*,*) nf_strerror(status)
        case ('areag')
           status = nf_put_vara_real(tranncid,varid,ostart,ocount,areag)
           if (status /= nf_noerr) write(*,*) nf_strerror(status)
        case ('areaa')
           status = nf_put_vara_real(tranncid,varid,ostart,ocount,areaa)
           if (status /= nf_noerr) write(*,*) nf_strerror(status)
        case ('areai')
           status = nf_put_vara_real(tranncid,varid,ostart,ocount,areai)
           if (status /= nf_noerr) write(*,*) nf_strerror(status)
        case ('areap')
           status = nf_put_vara_real(tranncid,varid,ostart,ocount,areap)
           if (status /= nf_noerr) write(*,*) nf_strerror(status)
        case default
           write(*,*) 'Cannot happen.'
        end select
     enddo
  enddo
  status = nf_close(tranncid)
end program poptrans
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine density(idm,jdm,tempp,salnp,den,kmt,klev,refp)
  integer idm,jdm,klev
  real tempp(idm,jdm),salnp(idm,jdm),ref_p
  integer kmt(idm,jdm)
  real den(idm,jdm)
  !===================================================================
  ! This program is used to calculate the density of each ocean grid
  ! point. Now the default is to calculate the potential density of
  ! the sea water -- sigma-0.
  ! The inputs are the dimensions of i, j, k direction, temperature
  ! and salinity, ocean grid mask (kmt). The outputs are density
  ! (sigma), temperature and salionity.
  !
  ! Coded by Aixue Hu, May 9, 2002.
  !
  ! input:
  !	idm -- i dimension of the input array
  !	jdm -- j dimension of the input array
  !	kmax -- k dimension of the input array
  !	tempp -- input potential temperature
  !	salnp -- input salinity (g/kg, or psu)
  !	kmt -- ocean depth file
  !	sigma -- output potential density
  !	ref_p -- reference pressure (decibars)
  !===================================================================

  integer i,j,k,lev
  real sigma(idm,jdm)


  ! --- calculate the potential density

  do j=1,jdm
     do i=1,idm
        if(kmt(i,j).ge.1.and.klev.le.kmt(i,j)) then
           call svannew(salnp(i,j),tempp(i,j),ref_p,sig0)
           den(i,j) = sig0
        else
           den(i,j) = 1.e30
        end if
     enddo
  enddo

  return
end subroutine density

!     subroutine dens(s,t,p0,sigma)

!     data c1,c2,c3,c4,c5,c6,c7/
!    . -1.36471E-01, 4.68181E-02, 8.07004E-01,-7.45353E-03,-2.94418E-03,
!    .  3.43570E-05, 3.48658E-05/
!
! --- sigma-theta as a function of temp (deg c) and salinity (mil)
! --- (friedrich-levitus 3rd degree polynomial fit)
!
!     sigma=(c1+c3*s+t*(c2+c5*s+t*(c4+c7*s+c6*t)))*1.e-3
!     return
!     end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine svannew(s,t,p0,sigma)
  !*
  !* Specific Volume Anomaly (steric anomaly) based on 1980 equation
  !* of state for seawater and 1978 practical salinity scale.
  !* References:
  !*        millero, et al (1980) deep-sea res.,27a,255-264
  !*        millero and poisson 1981,deep-sea res.,28a pp 625-629.
  !*   (both  references are also found in unesco report 38 -- 1981)
  !* ................................................................
  !*  The type of density anomaly (sigma) computed depends on the type of
  !*  temperature and pressure supplied:
  !*
  !*   in situ density:  t is in situ temperature
  !*                     p0 is in situ pressure
  !*
  !*   potential density: p0 is reference pressure
  !*                      t is temperature referenced already to p0
  !*
  !*   sigma-t: : p0 = 0
  !*              t is in situ temperature
  !* ................................................................
  !* units:
  !*       p0  :    pressure  [or ref pressure] (decibars)
  !*        t  :    temperature  [or pot temp]  (deg C)
  !*        s  :    salinity (ipss-78)
  !*     svan  :    spec. vol. anom.  (m**3/kg *1.0e-8)
  !*     sigma :    density anomaly   (kg/m**3)
  !* ................................................................
  !* check value: svan=981.3021 e-8 m**3/kg.  for s = 40 (ipss-78) ,
  !* t = 40 deg c, p0= 10000 decibars.
  !* check value: sigma = 59.82037  kg/m**3 for s = 40 (ipss-78) ,
  !* t = 40 deg c, p0= 10000 decibars.
  !* ................................................................
  !      implicit real*8 (a-h,o-z)

  real p,t,s,sig,sr,r1,r2,r3,r4
  real a,b,c,d,e,a1,b1,aw,bw,k,k0,kw,k35
  !*
  equivalence (e,d,b1), (bw,b,r3), (c,a1,r2)
  equivalence (aw,a,r1),(kw,k0,k)
  !*
  !      data r3500 /1028.1063/
  !      data r4 /4.8314e-4/    ! referred to as C in millero and poisson 1981
  !      data dr350 /28.106331/
  r3500=1028.1063
  r4   =4.8314e-4
  dr350=28.106331
  !*
  !*
  !* convert pressure to bars and take square root of salinity.
  !*
  p = p0 / 10.
  sr = sqrt(abs(s))
  !*
  !* pure water density at atmospheric pressure
  !*   bigg p.h.,(1967) br. j. applied physics 8 pp 521-537.
  !*
  r1 = ((((6.536332e-9*t-1.120083e-6)*t+1.001685e-4)*t &
       &      -9.095290e-3)*t+6.793952e-2)*t-28.263737
  !*
  !* seawater density at atmospheric press.
  !*    coefficients involving salinity:
  !*      r2 = A   in notation of millero and poisson 1981
  !*      r3 = B
  !*
  r2 = (((5.3875e-9*t-8.2467e-7)*t+7.6438e-5)*t-4.0899e-3)*t&
       &     +8.24493e-1
  !*
  r3 = (-1.6546e-6*t+1.0227e-4)*t-5.72466e-3
  !*
  !*  international one-atmosphere equation of state of seawater
  !*
  sig = (r4*s + r3*sr + r2)*s + r1
  !*
  !*  specific volume at atmospheric pressure
  !*
  v350p = 1.0 / r3500
  sva = -sig * v350p / (r3500+sig)
  sigma = sig + dr350
  !*
  !*  scale specific vol. anomaly to normally reported units
  !*
  svan = sva * 1.0e+8
  !*      svan = sva * 1.0e+7 !+8 crash out put
  if (p.eq.0.) return
  !*
  !* ............................................................
  !      high pressure equation of state for seawater
  !        millero, et al , 1980 dsr 27a, pp 255-264
  !               constant notation follows article
  !* ............................................................
  !* compute compression terms ...
  !*
  e = (9.1697e-10*t+2.0816e-8)*t-9.9348e-7
  bw = (5.2787e-8*t-6.12293e-6)*t+3.47718e-5
  b = bw + e * s
  !*
  d = 1.91075e-4
  c = (-1.6078e-6*t-1.0981e-5)*t+2.2838e-3
  aw = ((-5.77905e-7*t+1.16092e-4)*t+1.43713e-3)*t &
       &      -0.1194975
  a = (d * sr + c) * s + aw
  !*
  b1 = (-5.3009e-4*t+1.6483e-2)*t+7.944e-2
  a1 = ((-6.1670e-5*t+1.09987e-2)*t-0.603459)*t+54.6746
  kw = (((-5.155288e-5*t+1.360477e-2)*t-2.327105)*t &
       &      +148.4206)*t-1930.06
  k0 = (b1*sr + a1)*s + kw
  !*
  !* evaluate pressure polynomial
  !* ............................................................
  !c   k equals the secant bulk modulus of seawater
  !   dk = k(s,t,p) - k(35,0,p)
  !   k35 = k(35,0,p)
  !* ............................................................
  dk = (b * p + a) * p + k0
  k35 = (5.03217e-5*p+3.359406)*p+21582.27
  gam = p / k35
  pk = 1.0 - gam
  sva = sva*pk + (v350p+sva)*p*dk/(k35*(k35+dk))
  !*
  !*  scale specific vol. anamoly to normally reported units...
  !*
  svan = sva*1.0e+8
  !*      svan = sva*1.0e+7 !if +8 the out put crash because I/O format
  v350p = v350p * pk
  !* ....................................................................
  !* compute density anamoly with respect to 1000.0 kg/m**3
  !*  1) dr350: density anomaly at 35 (ipss-78), 0 deg. c and 0 decibars
  !*  2) dr35p: density anomaly 35 (ipss-78), 0 deg. c ,  pres. variation
  !*  3) dvan : density anomaly variations involving specfic vol. anamoly
  !* ....................................................................
  dr35p = gam / v350p
  dvan = sva / (v350p * (v350p+sva))
  sigma = dr350 + dr35p - dvan
  !*
  return
end subroutine svannew
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine check_exists(file)
  !
  ! Check and make sure 'file' exists
  !
  character(len=*),intent(in)::file
  !
  logical::exists
  !
  inquire(file=trim(file),exist=exists)
  if (.not.(exists)) then
     write(*,'(a,'' not found. STOP.'')') trim(file)
     call abort()
  endif
end subroutine check_exists
