program calculate_sea_level_change
  !
  ! netCDF I/O subroutines
  !
  use counters_netcdf_jfl
  use interfaces_netcdf_jfl
  use definitions_netcdf_jfl
  use state_ccsm_mod
  use kinds_mod
  !
  implicit none
  !
  ! Requirements:
  ! base state file;
  ! 1 to N input files;
  ! one auxiliary grid information file.
  ! All in netCDF format.
  !
  ! Argument one is "history" or "singlevar" - REQUIRED
  !
  ! Argument two is the "base" ocn file, ie, the one from which
  !  sea level change is to be computed
  ! If arg0 is "history", then
  !  Arguments 3 to N are the ocean files for which SLC should be computed (up to 1000 files)
  ! If arg0 is "singlevar", then
  !  Argument  3 is the TEMP data
  !  Argument  4 is the SALT data
  !  Argument  5 is the SSH data
  !  Argument  6 is the RHO data
  !
  ! Example history input file:
  ! history
  ! B06.62.ocn.1000-1100Ave.nc
  ! B06.62.ocn.1000.nc
  ! B06.62.ocn.1001.nc
  ! [...]
  ! B06.62.ocn.1099.nc
  ! B06.62.ocn.1100.nc
  !
  ! Example singlevar input file:
  ! singlevar
  ! b40.1850.track1.1deg.006.pop.h.0800-ave-1299.nc
  ! b40.20th.track1.1deg.008.pop.h.RHO.192001-192912.nc
  ! b40.20th.track1.1deg.008.pop.h.SALT.192001-192912.nc
  ! b40.20th.track1.1deg.008.pop.h.SSH.192001-192912.nc
  ! b40.20th.track1.1deg.008.pop.h.TEMP.192001-192912.nc
  !
  ! The auxiliary file for PCM contains POP grid information (pcm1_grid.nc)
  !
  ! Output is a netCDF file containing sea level change
  !
  real ,parameter :: spv = -1.e34
  !
  ! Base ocean state data
  !
  integer :: ncid_base,nlon_base,nlat_base,nlev_base,ntim_base
  real, dimension(:,:) , allocatable :: base_ssh
  real, dimension(:,:,:) , allocatable :: base_temp,base_salt,base_rho
  !
  ! Ocean data for which to calculate SLC, history file format
  !
  integer :: ncid_hist,nlon_hist,nlat_hist,nlev_hist,ntim_hist
  !
  ! Ocean data for which to calculate SLC, singlevar format
  !
  integer :: ncid_temp,nlon_temp,nlat_temp,nlev_temp,ntim_temp
  integer :: ncid_salt,nlon_salt,nlat_salt,nlev_salt,ntim_salt
  integer :: ncid_ssh,nlon_ssh,nlat_ssh,nlev_ssh,ntim_ssh
  integer :: ncid_rho,nlon_rho,nlat_rho,nlev_rho,ntim_rho
  !
  ! Working arrays
  !
  real, dimension(:)      , allocatable :: ocn_lon,ocn_lat,ocn_lev
  real, dimension(:,:)    , allocatable :: ocn_ssh
  real, dimension(:,:,:)  , allocatable :: ocn_temp,ocn_salt,ocn_rho
  real (r4),dimension(:,:), allocatable :: ocn_tempd,ocn_saltd
  real (r8),dimension(:,:), allocatable :: ocn_dend,ocn_den_t,ocn_den_s
  integer, dimension(:,:), allocatable :: kmt
  logical, dimension(:,:), allocatable :: kmtm
  !
  ! Grid data
  !
  integer :: ncid_grid,nlon_grid,nlat_grid,nlev_grid
  real, dimension(:) , allocatable :: grid_dz,grid_dzw
  real, dimension(:,:) , allocatable :: grid_dxt,grid_dyt
  !
  ! Sea level change data
  !
  integer :: ncid_slc,nlon_slc,nlat_slc,nlev_slc,ntim_slc
  real, dimension(:) , allocatable :: slc_lon,slc_lat,slc_lev
  real, dimension(:,:) , allocatable :: slc,dssh,steric
  real, dimension(:,:) , allocatable :: slc_t,steric_t
  real, dimension(:,:) , allocatable :: slc_s,steric_s
  real :: gave_steric,gave_dssh,gave_slc
  real :: gave_steric_t,gave_slc_t
  real :: gave_steric_s,gave_slc_s
  !
  integer :: i,j,k,l,n,ier,length,io_status,ifile,nfile
  integer :: cnt_hist,cnt_temp,cnt_salt,cnt_ssh,cnt_rho
  integer :: cnt_base,cnt_grid,cnt_slc
  real :: rho_base1,rho_ocn1,rho_ocn_t,rho_ocn_s
  real :: ocn_area,ocn_time,time
  logical :: exists,ispcm1,isccsm3,isccsm4
  character(len=256) :: dataset_type
  character(len=256) :: base_ocean_dataset,slc_dataset
  character(len=256) :: temp_dataset,salt_dataset,rho_dataset,ssh_dataset
  character(len=256),dimension(1000) :: ocn_history
  !
  exists = .false.
  ispcm1 = .false.
  isccsm3 = .false.
  isccsm4 = .false.
  !
  read (*,*) dataset_type
  if (dataset_type == "history") then
     write(*,*) 'History file format data.'
     ishist = .true.
  elseif (dataset_type == "singlevar") then
     write(*,*) 'Single variable format data.'
     issvar = .true.
  else
     write(*,*) 'Unknown dataset format. DIE!'
     call abort()
  endif
  !
  read (*,*) base_ocean_dataset
  call check_exists(base_ocean_dataset)
  write(*,'(a,'' is base ocean state.'')') trim(base_ocean_dataset)
  !
  if (ishist) then
     io_status = 0
     ifile = 1
     do while (io_status == 0)
        read(*,*,iostat=io_status) ocn_history(ifile)
        if (io_status == 0) then
           call check_exists(ocn_history(ifile))
           ifile = ifile + 1
        endif
     enddo
     nfile = ifile - 1
  endif
  if (issvar) then
     read(*,*) rho_dataset
     read(*,*) salt_dataset
     read(*,*) ssh_dataset
     read(*,*) temp_dataset
     call check_exists(temp_dataset)
     call check_exists(salt_dataset)
     call check_exists(ssh_dataset)
     call check_exists(rho_dataset)
  endif
  !
  ! Get base ocean data
  !
  call open_cdf(ncid_base,trim(base_ocean_dataset),.true.)
  call get_dims(ncid_base)
  call get_vars(ncid_base)
  !
  ! Get base ocean dimensions
  !
  do n=1,dim_counter
     length = len_trim(dim_info(n)%name)
     if((dim_info(n)%name(:length).eq.'lat').or.(dim_info(n)%name(:length).eq.'nlat')) then
        nlat_base = dim_info(n)%length
     elseif((dim_info(n)%name(:length).eq.'lon').or.(dim_info(n)%name(:length).eq.'nlon')) then
        nlon_base = dim_info(n)%length
     elseif((dim_info(n)%name(:length).eq.'lev').or.(dim_info(n)%name(:length).eq.'z_t')) then
        nlev_base = dim_info(n)%length
     elseif(dim_info(n)%name(:length).eq.'time') then
        ntim_base = dim_info(n)%length
     endif
  enddo
  !
  if (nlev_base.eq.32) ispcm1  = .true.
  if (nlev_base.eq.40) isccsm3 = .true.
  if (nlev_base.eq.60) isccsm4 = .true.
  !
  allocate(base_temp(nlon_base,nlat_base,nlev_base))
  allocate(base_salt(nlon_base,nlat_base,nlev_base))
  allocate(base_rho (nlon_base,nlat_base,nlev_base))
  allocate(base_ssh (nlon_base,nlat_base))
  !
  allocate(ocn_tempd(nlon_base,nlat_base))
  allocate(ocn_saltd(nlon_base,nlat_base))
  allocate(ocn_dend (nlon_base,nlat_base))
  allocate(ocn_den_t(nlon_base,nlat_base))
  allocate(ocn_den_s(nlon_base,nlat_base))
  !
  if (ispcm1) then
     !
     ! Get PCM grid data from PCM grid dataset
     call check_exists('pcm1_grid.nc')
     !
     call open_cdf(ncid_grid,'pcm1_grid.nc',.true.)
     call get_dims(ncid_grid)
     call get_vars(ncid_grid)
     !
     ! Get grid ocean dimensions
     do n=1,dim_counter
        length = len_trim(dim_info(n)%name)
        if(dim_info(n)%name(:length).eq.'jc') then
           nlat_grid = dim_info(n)%length
        elseif(dim_info(n)%name(:length).eq.'ic') then
           nlon_grid = dim_info(n)%length
        elseif(dim_info(n)%name(:length).eq.'zmid') then
           nlev_grid = dim_info(n)%length
        endif
     enddo
     !
     allocate(grid_dxt(nlon_grid,nlat_grid ))
     allocate(grid_dyt(nlon_grid,nlat_grid ))
     allocate(grid_dz ( nlev_grid))
     !
     ! Get grid DXT, DYT, and DZ
     call read_var(ncid_grid,'dxt',grid_dxt)
     call read_var(ncid_grid,'dyt',grid_dyt)
     call read_var(ncid_grid,'dz' ,grid_dz )
     call close_cdf(ncid_grid)
     write(*,'(''Grid information read from pcm1_grid.nc'')')
  else
     ! CCSM3/CCSM4 data
     allocate(grid_dxt(nlon_base,nlat_base))
     allocate(grid_dyt(nlon_base,nlat_base))
     allocate(kmt     (nlon_base,nlat_base))
     allocate(kmtm    (nlon_base,nlat_base))
     allocate(grid_dz (nlev_base),grid_dzw(nlev_base))
     call read_var(ncid_base,'DXT',grid_dxt)
     call read_var(ncid_base,'DYT',grid_dyt)
     call read_var(ncid_base,'dzw',grid_dzw)
     call read_var(ncid_base,'dz' ,grid_dz)
     call read_var(ncid_base,'KMT',kmt)
  endif
  !
  ! Get base ocean data T, S, and SSH
  cnt_base = 0
  do n=1,ntim_base
     time_counter = n
     cnt_base = cnt_base + 1
     if (ispcm1) then
        call read_var(ncid_base,'temp',base_temp)
        call read_var(ncid_base,'s' ,base_salt)
        call read_var(ncid_base,'ah1' ,base_ssh)
     else
        call read_var(ncid_base,'TEMP',base_temp)
        call read_var(ncid_base,'SALT',base_salt)
        call read_var(ncid_base,'SSH' ,base_ssh)
        if (isccsm4) call read_var(ncid_base,'RHO',base_rho)
     endif
  enddo
  call close_cdf(ncid_base)
  write(*,*) 'BASE S,T: ',(base_temp(230,150,k),base_salt(230,150,k),k=1,nlev_base)
  write(*,'(a,'' base state data read.'')') trim(base_ocean_dataset)
  !
  ! Where base T/S/SSH fields are land, set to zero, correct SALT from ppt to parts
  !
  where (base_ssh  > 1.e10)
     base_ssh  = 0
  end where
  where (base_temp > 1.e10)
     base_temp = 0
  end where
  where (base_salt > 1.e10)
     base_salt = 0
  end where
  where (base_salt > 1.)
     base_salt = base_salt*1.e-3
  end where
  if (isccsm4) then
     where (base_rho > 1.e10)
        base_rho = 0
     end where
  endif
  !
  if ((isccsm3).or.(isccsm4)) then
     do k=2,nlev_base
        grid_dzw(k)=grid_dzw(k-1)+grid_dzw(k)
     end do
     call init_state('mwjf',grid_dzw)
  endif
  if (isccsm3) then
     print*,'grid_dzw=',(grid_dzw(k),k=1,40)
     print*,'grid_dz=',(grid_dz(k),k=1,40)
     do k = 1,nlev_base
        kmtm = .False.
        where (kmt >= k)
           kmtm = .true.
        end where
        ocn_tempd=base_temp(:,:,k)
        ocn_saltd=base_salt(:,:,k)   ! *0.001
        call state_ccsm(k,k,ocn_tempd,ocn_saltd,ocn_dend,kmtm,spv)
        where (kmt < k)
           ocn_dend = 100
        end where
        base_rho(:,:,k)=ocn_dend
     end do
  endif
  if (isccsm4) then
     do k = 1,nlev_base
        kmtm = .False.
        where (kmt >= k)
           kmtm = .true.
        end where
        where (kmt < k)
           base_rho(:,:,k) = 100
        end where
     end do
  endif
  !
  ! If history format, step through ocean datasets
  !
  if (ishist) then
     do ifile = 1,nfile
        slc_dataset = 'DSL.'//trim(ocn_history(ifile))
        dim_counter = 0
        var_counter = 0
        call open_cdf(ncid_hist,trim(ocn_history(ifile)),.true.)
        call get_dims(ncid_hist)
        call get_vars(ncid_hist)
        !
        ! Get ocean dataset dimensions
        do n=1,dim_counter
           length = len_trim(dim_info(n)%name)
           if((dim_info(n)%name(:length).eq.'lat').or.(dim_info(n)%name(:length).eq.'nlat')) then
              nlat_hist = dim_info(n)%length
           elseif((dim_info(n)%name(:length).eq.'lon').or.(dim_info(n)%name(:length).eq.'nlon')) then
              nlon_hist = dim_info(n)%length
           elseif((dim_info(n)%name(:length).eq.'lev').or.(dim_info(n)%name(:length).eq.'z_t')) then
              nlev_hist = dim_info(n)%length
           elseif(dim_info(n)%name(:length).eq.'time') then
              ntim_hist = dim_info(n)%length
           endif
        enddo
        !
        allocate(ocn_temp(nlon_hist,nlat_hist,nlev_hist))
        allocate(ocn_salt(nlon_hist,nlat_hist,nlev_hist))
        allocate(ocn_rho (nlon_hist,nlat_hist,nlev_hist))
        allocate(ocn_ssh (nlon_hist,nlat_hist          ))
        allocate(slc     (nlon_hist,nlat_hist          ))
        allocate(slc_t   (nlon_hist,nlat_hist          ))
        allocate(slc_s   (nlon_hist,nlat_hist          ))
        allocate(steric  (nlon_hist,nlat_hist          ))
        allocate(steric_t(nlon_hist,nlat_hist          ))
        allocate(steric_s(nlon_hist,nlat_hist          ))
        allocate(dssh    (nlon_hist,nlat_hist          ))
        allocate(ocn_lon (nlon_hist                    ))
        allocate(ocn_lat (nlat_hist                    ))
        !
        ! Get coordinate variables
        !
        if (ispcm1) then
           call read_var(ncid_hist,'lon',ocn_lon)
           call read_var(ncid_hist,'lat',ocn_lat)
        else
           do i=1,nlon_hist
              ocn_lon(i)=i
           end do
           do i=1,nlat_hist
              ocn_lat(i)=i
           end do
        endif
        !
        ! Define SLC dataset
        !
        call create_cdf(ncid_slc,slc_dataset)
        !
        call def_dim(ncid_slc,'lon',nlon_hist)
        call def_dim(ncid_slc,'lat',nlat_hist)
        call def_dim(ncid_slc,'time',0)
        !
        call def_var(ncid_slc,'lon' ,'lon' ,'float',' ','longitude')
        call def_var(ncid_slc,'lat' ,'lat' ,'float',' ','latitude')
        call def_var(ncid_slc,'time','time','float','days since 0000-01-01','time')
        call def_att(ncid_slc,'time','calendar','noleap')
        !
        call def_var(ncid_slc,'slc' ,'lon,lat,time','float','cm','total sea level change')
        call def_att(ncid_slc,'slc','missing_value',spv)
        call def_att(ncid_slc,'slc','_FillValue',spv)
        call def_var(ncid_slc,'slc_t' ,'lon,lat,time','float','cm','total sea level change thermal only')
        call def_att(ncid_slc,'slc_t','missing_value',spv)
        call def_att(ncid_slc,'slc_t','_FillValue',spv)
        call def_var(ncid_slc,'slc_s' ,'lon,lat,time','float','cm','total sea level change salinity only')
        call def_att(ncid_slc,'slc_s','missing_value',spv)
        call def_att(ncid_slc,'slc_s','_FillValue',spv)
        call def_var(ncid_slc,'delta_ssh','lon,lat,time','float','cm','sea surface ht change' )
        call def_att(ncid_slc,'delta_ssh','missing_value',spv)
        call def_att(ncid_slc,'delta_ssh','_FillValue',spv)
        call def_var(ncid_slc,'delta_steric','lon,lat,time','float','cm','steric change' )
        call def_att(ncid_slc,'delta_steric','missing_value',spv)
        call def_att(ncid_slc,'delta_steric','_FillValue',spv)
        call def_var(ncid_slc,'delta_steric_t','lon,lat,time','float','cm','steric change thermal only' )
        call def_att(ncid_slc,'delta_steric_t','missing_value',spv)
        call def_att(ncid_slc,'delta_steric_t','_FillValue',spv)
        call def_var(ncid_slc,'delta_steric_s','lon,lat,time','float','cm','steric change salinity only' )
        call def_att(ncid_slc,'delta_steric_s','missing_value',spv)
        call def_att(ncid_slc,'delta_steric_s','_FillValue',spv)
        !
        call def_var(ncid_slc,'ave_slc','time','float','cm','global ave total SLC')
        call def_att(ncid_slc,'ave_slc','missing_value',spv)
        call def_att(ncid_slc,'ave_slc','_FillValue',spv)
        call def_var(ncid_slc,'ave_slc_t','time','float','cm','global ave total SLC thermal only')
        call def_att(ncid_slc,'ave_slc_t','missing_value',spv)
        call def_att(ncid_slc,'ave_slc_t','_FillValue',spv)
        call def_var(ncid_slc,'ave_slc_s','time','float','cm','global ave total SLC salinity only')
        call def_att(ncid_slc,'ave_slc_s','missing_value',spv)
        call def_att(ncid_slc,'ave_slc_s','_FillValue',spv)
        call def_var(ncid_slc,'ave_delta_ssh','time','float','cm','global ave SSH change' )
        call def_att(ncid_slc,'ave_delta_ssh','missing_value',spv)
        call def_att(ncid_slc,'ave_delta_ssh','_FillValue',spv)
        call def_var(ncid_slc,'ave_delta_steric','time','float','cm','global ave steric change')
        call def_att(ncid_slc,'ave_delta_steric','missing_value',spv)
        call def_att(ncid_slc,'ave_delta_steric','_FillValue',spv)
        call def_var(ncid_slc,'ave_delta_steric_t','time','float','cm','global ave steric change thermal only')
        call def_att(ncid_slc,'ave_delta_steric_t','missing_value',spv)
        call def_att(ncid_slc,'ave_delta_steric_t','_FillValue',spv)
        call def_var(ncid_slc,'ave_delta_steric_s','time','float','cm','global ave steric change salinity only')
        call def_att(ncid_slc,'ave_delta_steric_s','missing_value',spv)
        call def_att(ncid_slc,'ave_delta_steric_s','_FillValue',spv)
        !
        call end_def(ncid_slc)
        !
        call write_var(ncid_slc,'lon',ocn_lon)
        call write_var(ncid_slc,'lat',ocn_lat)
        !
        ! Begin calculation
        !
        do n = 1,ntim_hist
           time_counter = n
           if (ispcm1) then
              call read_var(ncid_hist,'s'   ,ocn_salt)
              call read_var(ncid_hist,'ah1' ,ocn_ssh)
              call read_var(ncid_hist,'temp',ocn_temp)
           else
              call read_var(ncid_hist,'RHO' ,ocn_rho)
              call read_var(ncid_hist,'SSH' ,ocn_ssh)
              call read_var(ncid_hist,'SALT',ocn_salt)
              call read_var(ncid_hist,'TEMP',ocn_temp)
           endif
           call read_var(ncid_hist,'time',ocn_time)
           call write_var(ncid_slc,'time',ocn_time)
           write(*,'(a,''     data at time '',i3,'' ('',f12.3,'') read.'')') &
                trim(ocn_history(ifile)),n,ocn_time
           !
           ! Where T/S/SSH fields are land, set to zero
           !
           where (ocn_ssh  > 1.e10)
              ocn_ssh  = 0
           end where
           if (ispcm1) then
              where (ocn_temp > 1.e10)
                 ocn_temp =   0.
              end where
              where (ocn_salt > 1.e10)
                 ocn_salt =   0.
              end where
           else
              where (ocn_rho  > 1.e10)
                 ocn_rho  = 100.
              end where
              where (ocn_temp > 1.e10)
                 ocn_temp =   0.
              end where
              where (ocn_salt > 1.e10)
                 ocn_salt =   0.
              end where
              where (ocn_salt > 1.)
                 ocn_salt = ocn_salt*1.e-3
              end where
           endif
           print*,ocn_temp(30,100,1),ocn_salt(30,100,1),ocn_rho(30,100,1)
           print*,base_temp(30,100,1),base_salt(30,100,1),base_rho(30,100,1)
           !
           ! Do the state calculations
           !
           steric   = 0.
           steric_t = 0.
           steric_s = 0.
           do k = 1,nlev_hist
              if ((isccsm3).or.(isccsm4)) then
                 kmtm = .False.
                 where (kmt >= k)
                    kmtm = .true.
                 end where
                 ! thermal effect only on SLC
                 ocn_tempd=ocn_temp(:,:,k)
                 ocn_saltd=base_salt(:,:,k) !*0.001
                 call state_ccsm(k,k,ocn_tempd,ocn_saltd,ocn_dend,kmtm,spv)
                 where (kmt < k)
                    ocn_dend = 100
                 end where
                 ocn_den_t=ocn_dend
                 ! salt effect only on SLC
                 ocn_tempd=base_temp(:,:,k)
                 ocn_saltd=ocn_salt(:,:,k)
                 call state_ccsm(k,k,ocn_tempd,ocn_saltd,ocn_dend,kmtm,spv)
                 where (kmt < k)
                    ocn_dend = 100
                 end where
                 ocn_den_s=ocn_dend
                 ocn_tempd=ocn_temp(:,:,k)
                 ocn_saltd=ocn_salt(:,:,k)
                 call state_ccsm(k,k,ocn_tempd,ocn_saltd,ocn_dend,kmtm,spv)
                 where (kmt < k)
                    ocn_dend = 100
                 end where
              endif
              if (ispcm1) then
                 do j = 1,nlat_hist
                    do i = 1,nlon_hist
                       call state_pcm(k,base_temp(i,j,k),base_salt(i,j,k),rho_base1)
                       call state_pcm(k, ocn_temp(i,j,k), ocn_salt(i,j,k),rho_ocn1)
                       call state_pcm(k, ocn_temp(i,j,k),base_salt(i,j,k),rho_ocn_t)
                       call state_pcm(k,base_temp(i,j,k), ocn_salt(i,j,k),rho_ocn_s)
                       if (i.eq.230.and.j.eq.150) then
                          print*,'rho_base= ',rho_base1,' rho_hist=',rho_ocn1
                       endif
                       steric(i,j)   = steric(i,j)   + (((rho_base1/rho_ocn1 )-1.)*grid_dz(k))
                       steric_t(i,j) = steric_t(i,j) + (((rho_base1/rho_ocn_t)-1.)*grid_dz(k))
                       steric_s(i,j) = steric_s(i,j) + (((rho_base1/rho_ocn_s)-1.)*grid_dz(k))
                    enddo
                 enddo
              else
                 do j = 1,nlat_hist
                    do i = 1,nlon_hist
                       rho_base1=base_rho(i,j,k)
                       rho_ocn1=ocn_dend(i,j) !rho_ocn(i,j,k)
                       rho_ocn_t=ocn_den_t(i,j)
                       rho_ocn_s=ocn_den_s(i,j)
                       if (i.eq.230.and.j.eq.150) then
                          print*,'rho_base= ',rho_base1,' rho_hist=',rho_ocn1
                       endif
                       steric(i,j)   = steric(i,j)   - ((( rho_ocn1/rho_base1)-1.)*grid_dz(k))
                       steric_t(i,j) = steric_t(i,j) - (((rho_ocn_t/rho_base1)-1.)*grid_dz(k))
                       steric_s(i,j) = steric_s(i,j) - (((rho_ocn_s/rho_base1)-1.)*grid_dz(k))
                       ! Old AHU calc
!                       steric(i,j)   = steric(i,j)   + (((rho_base1/rho_ocn1 )-1.)*grid_dz(k))
!                       steric_t(i,j) = steric_t(i,j) + (((rho_base1/rho_ocn_t)-1.)*grid_dz(k))
!                       steric_s(i,j) = steric_s(i,j) + (((rho_base1/rho_ocn_s)-1.)*grid_dz(k))
                    enddo
                 enddo
              endif
           enddo
           write(*,'(a,'' data at time '',i3,'' ('',f12.3,'') SLC calculated.'')') trim(slc_dataset),n,ocn_time
           !
           ! Calculate sea level change, areally averaged
           !
           dssh = ocn_ssh - base_ssh
           !           slc = steric + dssh
           !           slc_t = steric_t + dssh
           !           slc_s = steric_s + dssh
           where (base_ssh == 0.)
              dssh   = spv
              !              slc    = spv ; slc_t    = spv ; slc_s    = spv
              steric = spv
              steric_t = spv
              steric_s = spv
           endwhere
           !
           ocn_area = sum(grid_dxt*grid_dyt,mask=steric/=spv)
           !
           gave_steric   = sum(steric*grid_dxt*grid_dyt,mask=steric/=spv)
           gave_steric_t = sum(steric_t*grid_dxt*grid_dyt,mask=steric/=spv)
           gave_steric_s = sum(steric_s*grid_dxt*grid_dyt,mask=steric/=spv)
           gave_dssh     = sum(dssh*grid_dxt*grid_dyt,mask=dssh/=spv)
           !
           gave_steric   = gave_steric / ocn_area
           gave_steric_t = gave_steric_t / ocn_area
           gave_steric_s = gave_steric_s / ocn_area
           gave_dssh     = gave_dssh / ocn_area
           gave_slc      = gave_steric + gave_dssh
           gave_slc_t    = gave_steric_t + gave_dssh
           gave_slc_s    = gave_steric_s + gave_dssh

           slc           = gave_steric  + dssh
           slc_t         = gave_steric_t + dssh
           slc_s         = gave_steric_s + dssh
           where (base_ssh == 0.)
              slc = spv
              slc_t = spv
              slc_s = spv
           endwhere
           print*,'ocn_area = ',ocn_area ,'gave_steric=',gave_steric
           print*,'gave_dssh=',gave_dssh,'gave_slc=',gave_slc
           !
           call write_var(ncid_slc,'slc'         ,slc)
           call write_var(ncid_slc,'slc_t'       ,slc_t)
           call write_var(ncid_slc,'slc_s'       ,slc_s)
           call write_var(ncid_slc,'delta_ssh'   ,dssh)
           call write_var(ncid_slc,'delta_steric',steric)
           call write_var(ncid_slc,'delta_steric_t',steric_t)
           call write_var(ncid_slc,'delta_steric_s',steric_s)
           !
           call write_var(ncid_slc,'ave_slc'         ,gave_slc)
           call write_var(ncid_slc,'ave_slc_t'       ,gave_slc_t)
           call write_var(ncid_slc,'ave_slc_s'       ,gave_slc_s)
           call write_var(ncid_slc,'ave_delta_ssh'   ,gave_dssh)
           call write_var(ncid_slc,'ave_delta_steric',gave_steric)
           call write_var(ncid_slc,'ave_delta_steric_t',gave_steric_t)
           call write_var(ncid_slc,'ave_delta_steric_s',gave_steric_s)
           !        call write_var(ncid_slc,'rho_base',base_rho(1,:,:))
           !
        enddo
        call close_cdf(ncid_slc)
        deallocate(ocn_temp)
        deallocate(ocn_salt)
        deallocate(ocn_ssh )
        deallocate(slc     )
        deallocate(steric  )
        deallocate(dssh    )
        deallocate(ocn_lon )
        deallocate(ocn_lat )
     enddo
  endif
  !
  ! If singlevar format, step through time
  !
  if (issvar) then
     slc_dataset = 'DSL.'//trim(temp_dataset)
!     dim_counter = 0
!     var_counter = 0
     !
     ! Get information from each of TEMP, SALT, SSH, and RHO datasets
     call open_cdf(ncid_temp,trim(temp_dataset),.true.)
     call get_dims(ncid_temp)
     call get_vars(ncid_temp)
     do n=1,dim_counter
        length = len_trim(dim_info(n)%name)
        if((dim_info(n)%name(:length).eq.'lat').or.(dim_info(n)%name(:length).eq.'nlat')) then
           nlat_temp = dim_info(n)%length
        elseif((dim_info(n)%name(:length).eq.'lon').or.(dim_info(n)%name(:length).eq.'nlon')) then
           nlon_temp = dim_info(n)%length
        elseif((dim_info(n)%name(:length).eq.'lev').or.(dim_info(n)%name(:length).eq.'z_t')) then
           nlev_temp = dim_info(n)%length
        elseif(dim_info(n)%name(:length).eq.'time') then
           ntim_temp = dim_info(n)%length
        endif
     enddo
     !
     call open_cdf(ncid_salt,trim(salt_dataset),.true.)
     call get_dims(ncid_salt)
     call get_vars(ncid_salt)
     do n=1,dim_counter
        length = len_trim(dim_info(n)%name)
        if((dim_info(n)%name(:length).eq.'lat').or.(dim_info(n)%name(:length).eq.'nlat')) then
           nlat_salt = dim_info(n)%length
        elseif((dim_info(n)%name(:length).eq.'lon').or.(dim_info(n)%name(:length).eq.'nlon')) then
           nlon_salt = dim_info(n)%length
        elseif((dim_info(n)%name(:length).eq.'lev').or.(dim_info(n)%name(:length).eq.'z_t')) then
           nlev_salt = dim_info(n)%length
        elseif(dim_info(n)%name(:length).eq.'time') then
           ntim_salt = dim_info(n)%length
        endif
     enddo
     !
     call open_cdf(ncid_rho,trim(rho_dataset),.true.)
     call get_dims(ncid_rho)
     call get_vars(ncid_rho)
     do n=1,dim_counter
        length = len_trim(dim_info(n)%name)
        if((dim_info(n)%name(:length).eq.'lat').or.(dim_info(n)%name(:length).eq.'nlat')) then
           nlat_rho = dim_info(n)%length
        elseif((dim_info(n)%name(:length).eq.'lon').or.(dim_info(n)%name(:length).eq.'nlon')) then
           nlon_rho = dim_info(n)%length
        elseif((dim_info(n)%name(:length).eq.'lev').or.(dim_info(n)%name(:length).eq.'z_t')) then
           nlev_rho = dim_info(n)%length
        elseif(dim_info(n)%name(:length).eq.'time') then
           ntim_rho = dim_info(n)%length
        endif
     enddo
     !
     call open_cdf(ncid_ssh,trim(ssh_dataset),.true.)
     call get_dims(ncid_ssh)
     call get_vars(ncid_ssh)
     do n=1,dim_counter
        length = len_trim(dim_info(n)%name)
        if((dim_info(n)%name(:length).eq.'lat').or.(dim_info(n)%name(:length).eq.'nlat')) then
           nlat_ssh = dim_info(n)%length
        elseif((dim_info(n)%name(:length).eq.'lon').or.(dim_info(n)%name(:length).eq.'nlon')) then
           nlon_ssh = dim_info(n)%length
        elseif((dim_info(n)%name(:length).eq.'lev').or.(dim_info(n)%name(:length).eq.'z_t')) then
           nlev_ssh = dim_info(n)%length
        elseif(dim_info(n)%name(:length).eq.'time') then
           ntim_ssh = dim_info(n)%length
        endif
     enddo
     !
     ! Assume all dimensions are the same (except for SSH time); put in code here to verify
     !
     allocate(ocn_temp(nlon_temp,nlat_temp,nlev_temp))
     allocate(ocn_salt(nlon_salt,nlat_salt,nlev_salt))
     allocate(ocn_rho (nlon_rho ,nlat_rho ,nlev_rho ))
     allocate(ocn_ssh (nlon_ssh ,nlat_ssh           ))
     !
     allocate(slc     (nlon_temp,nlat_temp         ))
     allocate(slc_t   (nlon_temp,nlat_temp         ))
     allocate(slc_s   (nlon_temp,nlat_temp         ))
     allocate(steric  (nlon_temp,nlat_temp         ))
     allocate(steric_t(nlon_temp,nlat_temp         ))
     allocate(steric_s(nlon_temp,nlat_temp         ))
     allocate(dssh    (nlon_temp,nlat_temp         ))
     allocate(ocn_lon (nlon_temp                  ))
     allocate(ocn_lat (nlat_temp                  ))
     !
     ! Get coordinate variables
     !
     if (ispcm1) then
        call read_var(ncid_temp,'lon',ocn_lon)
        call read_var(ncid_temp,'lat',ocn_lat)
     else
        do i=1,nlon_temp
           ocn_lon(i)=i
        end do
        do i=1,nlat_temp
           ocn_lat(i)=i
        end do
     endif
     !
     ! Define SLC dataset
     !
     call create_cdf(ncid_slc,trim(slc_dataset))
     write(*,'(''call create_cdf: '',i10,5x,a80)') ncid_slc,slc_dataset
     !
     call def_dim(ncid_slc,'lon',nlon_temp)
     call def_dim(ncid_slc,'lat',nlat_temp)
     call def_dim(ncid_slc,'time',0)
     !
     call def_var(ncid_slc,'lon' ,'lon' ,'float',' ','longitude')
     call def_var(ncid_slc,'lat' ,'lat' ,'float',' ','latitude')
     call def_var(ncid_slc,'time','time','float','days since 0000-01-01','time')
     call def_att(ncid_slc,'time','calendar','noleap')
     !
     call def_var(ncid_slc,'slc' ,'lon,lat,time','float','cm','total sea level change')
     call def_att(ncid_slc,'slc','missing_value',spv)
     call def_att(ncid_slc,'slc','_FillValue',spv)
     call def_var(ncid_slc,'slc_t' ,'lon,lat,time','float','cm','total sea level change thermal only')
     call def_att(ncid_slc,'slc_t','missing_value',spv)
     call def_att(ncid_slc,'slc_t','_FillValue',spv)
     call def_var(ncid_slc,'slc_s' ,'lon,lat,time','float','cm','total sea level change salinity only')
     call def_att(ncid_slc,'slc_s','missing_value',spv)
     call def_att(ncid_slc,'slc_s','_FillValue',spv)
     call def_var(ncid_slc,'delta_ssh','lon,lat,time','float','cm','sea surface ht change' )
     call def_att(ncid_slc,'delta_ssh','missing_value',spv)
     call def_att(ncid_slc,'delta_ssh','_FillValue',spv)
     call def_var(ncid_slc,'delta_steric','lon,lat,time','float','cm','steric change' )
     call def_att(ncid_slc,'delta_steric','missing_value',spv)
     call def_att(ncid_slc,'delta_steric','_FillValue',spv)
     call def_var(ncid_slc,'delta_steric_t','lon,lat,time','float','cm','steric change thermal only' )
     call def_att(ncid_slc,'delta_steric_t','missing_value',spv)
     call def_att(ncid_slc,'delta_steric_t','_FillValue',spv)
     call def_var(ncid_slc,'delta_steric_s','lon,lat,time','float','cm','steric change salinity only' )
     call def_att(ncid_slc,'delta_steric_s','missing_value',spv)
     call def_att(ncid_slc,'delta_steric_s','_FillValue',spv)
     !
     call def_var(ncid_slc,'ave_slc','time','float','cm','global ave total SLC')
     call def_att(ncid_slc,'ave_slc','missing_value',spv)
     call def_att(ncid_slc,'ave_slc','_FillValue',spv)
     call def_var(ncid_slc,'ave_slc_t','time','float','cm','global ave total SLC thermal only')
     call def_att(ncid_slc,'ave_slc_t','missing_value',spv)
     call def_att(ncid_slc,'ave_slc_t','_FillValue',spv)
     call def_var(ncid_slc,'ave_slc_s','time','float','cm','global ave total SLC salinity only')
     call def_att(ncid_slc,'ave_slc_s','missing_value',spv)
     call def_att(ncid_slc,'ave_slc_s','_FillValue',spv)
     call def_var(ncid_slc,'ave_delta_ssh','time','float','cm','global ave SSH change' )
     call def_att(ncid_slc,'ave_delta_ssh','missing_value',spv)
     call def_att(ncid_slc,'ave_delta_ssh','_FillValue',spv)
     call def_var(ncid_slc,'ave_delta_steric','time','float','cm','global ave steric change')
     call def_att(ncid_slc,'ave_delta_steric','missing_value',spv)
     call def_att(ncid_slc,'ave_delta_steric','_FillValue',spv)
     call def_var(ncid_slc,'ave_delta_steric_t','time','float','cm','global ave steric change thermal only')
     call def_att(ncid_slc,'ave_delta_steric_t','missing_value',spv)
     call def_att(ncid_slc,'ave_delta_steric_t','_FillValue',spv)
     call def_var(ncid_slc,'ave_delta_steric_s','time','float','cm','global ave steric change salinity only')
     call def_att(ncid_slc,'ave_delta_steric_s','missing_value',spv)
     call def_att(ncid_slc,'ave_delta_steric_s','_FillValue',spv)
     !
     call end_def(ncid_slc)
     !
     call write_var(ncid_slc,'lon',ocn_lon)
     call write_var(ncid_slc,'lat',ocn_lat)
     !
     ! Begin calculation
     !
     do n = 1,ntim_temp
        time_counter = n
        if (ispcm1) then
           call read_var(ncid_salt,'s'   ,ocn_salt)
           call read_var(ncid_ssh ,'ah1' ,ocn_ssh)
           call read_var(ncid_temp,'temp',ocn_temp)
        else
           call read_var(ncid_rho ,'RHO' ,ocn_rho)
           call read_var(ncid_ssh ,'SSH' ,ocn_ssh)
           call read_var(ncid_salt,'SALT',ocn_salt)
           call read_var(ncid_temp,'TEMP',ocn_temp)
        endif
        call read_var(ncid_temp,'time',ocn_time)
        call write_var(ncid_slc,'time' ,ocn_time)
        write(*,'(a,''     data at time '',i3,'' ('',f12.3,'') read.'')') &
             trim(temp_dataset),n,ocn_time
        !
        ! Where T/S/SSH fields are land, set to zero
        !
        where (ocn_ssh  > 1.e10)
           ocn_ssh  = 0
        end where
        if (ispcm1) then
           where (ocn_temp > 1.e10)
              ocn_temp = 0
           end where
           where (ocn_salt > 1.e10)
              ocn_salt = 0
           end where
        else
           where (ocn_rho > 1.e10)
              ocn_rho = 100
           end where
           where (ocn_temp > 1.e10)
              ocn_temp = 0
           end where
           where (ocn_salt > 1.e10)
              ocn_salt = 0
           end where
           where (ocn_salt > 1.)
              ocn_salt = ocn_salt*1.e-3
           end where
        endif
        !
        ! Do the state calculations
        !
        steric = 0.
        steric_t = 0.
        steric_s = 0.
        do k = 1,nlev_temp
           if ((isccsm3).or.(isccsm4)) then
              kmtm = .False.
              where (kmt >= k)
                 kmtm = .true.
              end where
              ! thermal effect only on SLC
              ocn_tempd=ocn_temp(:,:,k)
              ocn_saltd=base_salt(:,:,k)!*0.001
              call state_ccsm(k,k,ocn_tempd,ocn_saltd,ocn_dend,kmtm,spv)
              where (kmt < k)
                 ocn_dend = 100
              end where
              ocn_den_t=ocn_dend
              ! salt effect only on SLC
              ocn_tempd=base_temp(:,:,k)
              ocn_saltd=ocn_salt(:,:,k)
              call state_ccsm(k,k,ocn_tempd,ocn_saltd,ocn_dend,kmtm,spv)
              where (kmt < k)
                 ocn_dend = 100
              end where
              ocn_den_s=ocn_dend
              ocn_tempd=ocn_temp(:,:,k)
              ocn_saltd=ocn_salt(:,:,k)
              call state_ccsm(k,k,ocn_tempd,ocn_saltd,ocn_dend,kmtm,spv)
              where (kmt < k)
                 ocn_dend = 100
              end where
           endif
           if (ispcm1) then
              do j = 1,nlat_temp
                 do i = 1,nlon_temp
                    call state_pcm(k,base_temp(i,j,k),base_salt(i,j,k),rho_base1)
                    call state_pcm(k, ocn_temp(i,j,k), ocn_salt(i,j,k), rho_ocn1)
                    call state_pcm(k, ocn_temp(i,j,k),base_salt(i,j,k), rho_ocn_t)
                    call state_pcm(k,base_temp(i,j,k), ocn_salt(i,j,k), rho_ocn_s)
                    steric(i,j) = steric(i,j) + (((rho_base1/rho_ocn1)-1.)*grid_dz(k))
                    steric_t(i,j) = steric_t(i,j) + (((rho_base1/rho_ocn_t)-1.)*grid_dz(k))
                    steric_s(i,j) = steric_s(i,j) + (((rho_base1/rho_ocn_s)-1.)*grid_dz(k))
                 enddo
              enddo
           else
              do j = 1,nlat_temp
                 do i = 1,nlon_temp
                    rho_base1=base_rho(i,j,k)
                    rho_ocn1=ocn_dend(i,j) !rho_ocn(i,j,k)
                    !            rho_ocn1=rho_ocn(i,j,k)
                    rho_ocn_t=ocn_den_t(i,j)
                    rho_ocn_s=ocn_den_s(i,j)
                    steric(i,j) = steric(i,j) + (((rho_base1/rho_ocn1)-1.)*grid_dz(k))
                    steric_t(i,j) = steric_t(i,j) + (((rho_base1/rho_ocn_t)-1.)*grid_dz(k))
                    steric_s(i,j) = steric_s(i,j) + (((rho_base1/rho_ocn_s)-1.)*grid_dz(k))
                 enddo
              enddo
           endif
        enddo
        write(*,'(a,'' data at time '',i3,'' ('',f12.3,'') SLC calculated.'')') trim(slc_dataset),n,ocn_time
        !
        ! Calculate sea level change, areally averaged
        !
        dssh = ocn_ssh - base_ssh
        !        slc = steric + dssh
        !        slc_t = steric_t + dssh
        !        slc_s = steric_s + dssh
        where (base_ssh == 0.)
           dssh   = spv
           !              slc    = spv ; slc_t    = spv ; slc_s    = spv
           steric = spv
           steric_t = spv
           steric_s = spv
        endwhere
        !
        ocn_area = sum(grid_dxt*grid_dyt,mask=steric/=spv)
        !
        gave_steric   = sum(steric*grid_dxt*grid_dyt,mask=steric/=spv)
        gave_steric_t = sum(steric_t*grid_dxt*grid_dyt,mask=steric/=spv)
        gave_steric_s = sum(steric_s*grid_dxt*grid_dyt,mask=steric/=spv)
        gave_dssh     = sum(dssh*grid_dxt*grid_dyt,mask=dssh/=spv)
        !
        gave_steric   = gave_steric / ocn_area
        gave_steric_t = gave_steric_t / ocn_area
        gave_steric_s = gave_steric_s / ocn_area
        gave_dssh     = gave_dssh / ocn_area
        gave_slc      = gave_steric + gave_dssh
        gave_slc_t    = gave_steric_t + gave_dssh
        gave_slc_s    = gave_steric_s + gave_dssh

        slc           = gave_steric  + dssh
        slc_t         = gave_steric_t + dssh
        slc_s         = gave_steric_s + dssh
        print*,'ocn_area = ',ocn_area ,'gave_steric=',gave_steric
        print*,'gave_dssh=',gave_dssh,'gave_slc=',gave_slc
        !
        call write_var(ncid_slc,'slc'         ,slc)
        call write_var(ncid_slc,'slc_t'       ,slc_t)
        call write_var(ncid_slc,'slc_s'       ,slc_s)
        call write_var(ncid_slc,'delta_ssh'   ,dssh)
        call write_var(ncid_slc,'delta_steric',steric)
        call write_var(ncid_slc,'delta_steric_t',steric_t)
        call write_var(ncid_slc,'delta_steric_s',steric_s)
        !
        call write_var(ncid_slc,'ave_slc'         ,gave_slc)
        call write_var(ncid_slc,'ave_slc_t'       ,gave_slc_t)
        call write_var(ncid_slc,'ave_slc_s'       ,gave_slc_s)
        call write_var(ncid_slc,'ave_delta_ssh'   ,gave_dssh)
        call write_var(ncid_slc,'ave_delta_steric',gave_steric)
        call write_var(ncid_slc,'ave_delta_steric_t',gave_steric_t)
        call write_var(ncid_slc,'ave_delta_steric_s',gave_steric_s)
        !        call write_var(ncid_slc,'rho_base',base_rho(1,:,:))
        !
     enddo
     call close_cdf(ncid_slc)
     deallocate(ocn_temp)
     deallocate(ocn_salt)
     deallocate(ocn_ssh )
     deallocate(slc     )
     deallocate(steric  )
     deallocate(dssh    )
     deallocate(ocn_lon )
     deallocate(ocn_lat )
  endif
end program calculate_sea_level_change
