topography3_nc.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Subroutine :  t o p o g r a p h y 3 _ n c
!
!> @file
!!
!! Definition of the initial surface and bedrock topography
!! (including gradients) and of the horizontal grid spacings dxi, deta.
!! For initial topography from previous simulation (read in NetCDF format).
!!
!! @section Copyright
!!
!! Copyright 2009-2016 Ralf Greve
!!
!! @section License
!!
!! This file is part of SICOPOLIS.
!!
!! SICOPOLIS is free software: you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! SICOPOLIS is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with SICOPOLIS.  If not, see <http://www.gnu.org/licenses/>.
!<
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!-------------------------------------------------------------------------------
!> Definition of the initial surface and bedrock topography
!! (including gradients) and of the horizontal grid spacings dxi, deta.
!! For initial topography from previous simulation (read in NetCDF format).
!<------------------------------------------------------------------------------
subroutine topography3_nc(dxi, deta, z_sl, anfdatname)

use sico_types
use sico_variables
use sico_vars
use netcdf
use nc_check

implicit none

character(len=100), intent(in) :: anfdatname

real(dp),          intent(out) :: dxi, deta, z_sl

integer(i4b) :: i, j, kc, kt, kr, n
integer(i4b) :: ios
character :: ch_dummy
character(len=256) :: anfdatname_with_path

integer(i4b) :: ncid, ncv
!     ncid:      ID of the output file
!     ncv:       Variable ID

integer(i2b), dimension(0:IMAX,0:JMAX) :: maske_conv, n_cts_conv, kc_cts_conv
real(sp) :: time_conv, dummy_conv, z_sl_conv, &
            v_tot_conv, a_grounded_conv, a_floating_conv, &
            h_r_conv, &
            xi_conv(0:imax), eta_conv(0:jmax), &
            sigma_level_c_conv(0:kcmax), sigma_level_t_conv(0:ktmax), &
            sigma_level_r_conv(0:krmax)
real(sp), dimension(0:IMAX,0:JMAX) :: lambda_conv, phi_conv, &
            lon_conv, lat_conv, &
            temp_s_conv, as_perp_conv, &
            zs_conv, zm_conv, zb_conv, zl_conv, &
            h_cold_conv, h_temp_conv, h_conv, &
            q_bm_conv, q_tld_conv, &
            am_perp_conv, &
            qx_conv, qy_conv, &
            dzs_dtau_conv, dzm_dtau_conv, dzb_dtau_conv, dzl_dtau_conv, &
            dh_c_dtau_conv, dh_t_dtau_conv, dh_dtau_conv, &
            vx_b_g_conv, vy_b_g_conv, vz_b_conv, vh_b_conv, &
            vx_s_g_conv, vy_s_g_conv, vz_s_conv, vh_s_conv, &
            temp_b_conv, temph_b_conv, &
            p_b_w_conv, h_w_conv, q_gl_g_conv
real(sp), dimension(0:IMAX,0:JMAX,0:KCMAX) :: vx_c_conv, vy_c_conv, vz_c_conv, &
                                              temp_c_conv, age_c_conv, &
                                              enth_c_conv, omega_c_conv, &
                                              enh_c_conv
real(sp), dimension(0:IMAX,0:JMAX,0:KTMAX) :: vx_t_conv, vy_t_conv, vz_t_conv, &
                                              omega_t_conv, age_t_conv, &
                                              enth_t_conv, &
                                              enh_t_conv
real(sp), dimension(0:IMAX,0:JMAX,0:KRMAX) :: temp_r_conv

!-------- Read data from time-slice file of previous simulation --------

anfdatname_with_path = trim(anfdatpath)//'/'//trim(anfdatname)

call check( nf90_open(trim(anfdatname_with_path), nf90_nowrite, ncid) )

call check( nf90_inq_varid(ncid, 'time', ncv) )
call check( nf90_get_var(ncid, ncv, time_conv) )

if ( nf90_inq_varid(ncid, 'delta_ts', ncv) == nf90_noerr ) then
   call check( nf90_get_var(ncid, ncv, dummy_conv) )
else if ( nf90_inq_varid(ncid, 'glac_index', ncv) == nf90_noerr ) then
   call check( nf90_get_var(ncid, ncv, dummy_conv) )
else
   continue
end if

call check( nf90_inq_varid(ncid, 'z_sl', ncv) )
call check( nf90_get_var(ncid, ncv, z_sl_conv) )

call check( nf90_inq_varid(ncid, 'xi', ncv) )
call check( nf90_get_var(ncid, ncv, xi_conv) )

call check( nf90_inq_varid(ncid, 'eta', ncv) )
call check( nf90_get_var(ncid, ncv, eta_conv) )

call check( nf90_inq_varid(ncid, 'sigma_level_c', ncv) )
call check( nf90_get_var(ncid, ncv, sigma_level_c_conv) )

call check( nf90_inq_varid(ncid, 'sigma_level_t', ncv) )
call check( nf90_get_var(ncid, ncv, sigma_level_t_conv) )

call check( nf90_inq_varid(ncid, 'sigma_level_r', ncv) )
call check( nf90_get_var(ncid, ncv, sigma_level_r_conv) )

call check( nf90_inq_varid(ncid, 'lon', ncv) )
call check( nf90_get_var(ncid, ncv, lon_conv) )

call check( nf90_inq_varid(ncid, 'lat', ncv) )
call check( nf90_get_var(ncid, ncv, lat_conv) )

call check( nf90_inq_varid(ncid, 'lambda', ncv) )
call check( nf90_get_var(ncid, ncv, lambda_conv) )

call check( nf90_inq_varid(ncid, 'phi', ncv) )
call check( nf90_get_var(ncid, ncv, phi_conv) )

call check( nf90_inq_varid(ncid, 'temp_s', ncv) )
call check( nf90_get_var(ncid, ncv, temp_s_conv) )

call check( nf90_inq_varid(ncid, 'as_perp', ncv) )
call check( nf90_get_var(ncid, ncv, as_perp_conv) )

call check( nf90_inq_varid(ncid, 'maske', ncv) )
call check( nf90_get_var(ncid, ncv, maske_conv) )

call check( nf90_inq_varid(ncid, 'n_cts', ncv) )
call check( nf90_get_var(ncid, ncv, n_cts_conv) )

if ( nf90_inq_varid(ncid, 'kc_cts', ncv) == nf90_noerr ) then
   call check( nf90_get_var(ncid, ncv, kc_cts_conv) )
else
   go to 110
end if

call check( nf90_inq_varid(ncid, 'zs', ncv) )
call check( nf90_get_var(ncid, ncv, zs_conv) )

call check( nf90_inq_varid(ncid, 'zm', ncv) )
call check( nf90_get_var(ncid, ncv, zm_conv) )

call check( nf90_inq_varid(ncid, 'zb', ncv) )
call check( nf90_get_var(ncid, ncv, zb_conv) )

call check( nf90_inq_varid(ncid, 'zl', ncv) )
call check( nf90_get_var(ncid, ncv, zl_conv) )

call check( nf90_inq_varid(ncid, 'H_cold', ncv) )
call check( nf90_get_var(ncid, ncv, h_cold_conv) )

call check( nf90_inq_varid(ncid, 'H_temp', ncv) )
call check( nf90_get_var(ncid, ncv, h_temp_conv) )

call check( nf90_inq_varid(ncid, 'H', ncv) )
call check( nf90_get_var(ncid, ncv, h_conv) )

call check( nf90_inq_varid(ncid, 'H_R', ncv) )
call check( nf90_get_var(ncid, ncv, h_r_conv) )

call check( nf90_inq_varid(ncid, 'vx_c', ncv) )
call check( nf90_get_var(ncid, ncv, vx_c_conv) )

call check( nf90_inq_varid(ncid, 'vy_c', ncv) )
call check( nf90_get_var(ncid, ncv, vy_c_conv) )

call check( nf90_inq_varid(ncid, 'vz_c', ncv) )
call check( nf90_get_var(ncid, ncv, vz_c_conv) )

call check( nf90_inq_varid(ncid, 'vx_t', ncv) )
call check( nf90_get_var(ncid, ncv, vx_t_conv) )

call check( nf90_inq_varid(ncid, 'vy_t', ncv) )
call check( nf90_get_var(ncid, ncv, vy_t_conv) )

call check( nf90_inq_varid(ncid, 'vz_t', ncv) )
call check( nf90_get_var(ncid, ncv, vz_t_conv) )

call check( nf90_inq_varid(ncid, 'temp_c', ncv) )
call check( nf90_get_var(ncid, ncv, temp_c_conv) )

call check( nf90_inq_varid(ncid, 'omega_t', ncv) )
call check( nf90_get_var(ncid, ncv, omega_t_conv) )

call check( nf90_inq_varid(ncid, 'temp_r', ncv) )
call check( nf90_get_var(ncid, ncv, temp_r_conv) )

if ( nf90_inq_varid(ncid, 'enth_c', ncv) == nf90_noerr ) then
   call check( nf90_get_var(ncid, ncv, enth_c_conv) )
else
   go to 110
end if

if ( nf90_inq_varid(ncid, 'enth_t', ncv) == nf90_noerr ) then
   call check( nf90_get_var(ncid, ncv, enth_t_conv) )
else
   go to 110
end if

if ( nf90_inq_varid(ncid, 'omega_c', ncv) == nf90_noerr ) then
   call check( nf90_get_var(ncid, ncv, omega_c_conv) )
else
   go to 110
end if

if ( nf90_inq_varid(ncid, 'enh_c', ncv) == nf90_noerr ) then
   call check( nf90_get_var(ncid, ncv, enh_c_conv) )
else
   go to 110
end if

if ( nf90_inq_varid(ncid, 'enh_t', ncv) == nf90_noerr ) then
   call check( nf90_get_var(ncid, ncv, enh_t_conv) )
else
   go to 110
end if

call check( nf90_inq_varid(ncid, 'Q_bm', ncv) )
call check( nf90_get_var(ncid, ncv, q_bm_conv) )

call check( nf90_inq_varid(ncid, 'Q_tld', ncv) )
call check( nf90_get_var(ncid, ncv, q_tld_conv) )

call check( nf90_inq_varid(ncid, 'am_perp', ncv) )
call check( nf90_get_var(ncid, ncv, am_perp_conv) )

call check( nf90_inq_varid(ncid, 'qx', ncv) )
call check( nf90_get_var(ncid, ncv, qx_conv) )

call check( nf90_inq_varid(ncid, 'qy', ncv) )
call check( nf90_get_var(ncid, ncv, qy_conv) )

call check( nf90_inq_varid(ncid, 'age_c', ncv) )
call check( nf90_get_var(ncid, ncv, age_c_conv) )

call check( nf90_inq_varid(ncid, 'age_t', ncv) )
call check( nf90_get_var(ncid, ncv, age_t_conv) )

call check( nf90_inq_varid(ncid, 'dzs_dtau', ncv) )
call check( nf90_get_var(ncid, ncv, dzs_dtau_conv) )

call check( nf90_inq_varid(ncid, 'dzm_dtau', ncv) )
call check( nf90_get_var(ncid, ncv, dzm_dtau_conv) )

call check( nf90_inq_varid(ncid, 'dzb_dtau', ncv) )
call check( nf90_get_var(ncid, ncv, dzb_dtau_conv) )

call check( nf90_inq_varid(ncid, 'dzl_dtau', ncv) )
call check( nf90_get_var(ncid, ncv, dzl_dtau_conv) )

call check( nf90_inq_varid(ncid, 'dH_c_dtau', ncv) )
call check( nf90_get_var(ncid, ncv, dh_c_dtau_conv) )

call check( nf90_inq_varid(ncid, 'dH_t_dtau', ncv) )
call check( nf90_get_var(ncid, ncv, dh_t_dtau_conv) )

call check( nf90_inq_varid(ncid, 'dH_dtau', ncv) )
call check( nf90_get_var(ncid, ncv, dh_dtau_conv) )

if ( nf90_inq_varid(ncid, 'vx_b_g', ncv) == nf90_noerr ) then
   call check( nf90_get_var(ncid, ncv, vx_b_g_conv) )
else
   go to 110
end if

if ( nf90_inq_varid(ncid, 'vy_b_g', ncv) == nf90_noerr ) then
   call check( nf90_get_var(ncid, ncv, vy_b_g_conv) )
else
   go to 110
end if

if ( nf90_inq_varid(ncid, 'vz_b', ncv) == nf90_noerr ) then
   call check( nf90_get_var(ncid, ncv, vz_b_conv) )
else
   go to 110
end if

if ( nf90_inq_varid(ncid, 'vh_b', ncv) == nf90_noerr ) then
   call check( nf90_get_var(ncid, ncv, vh_b_conv) )
else
   go to 110
end if

if ( nf90_inq_varid(ncid, 'vx_s_g', ncv) == nf90_noerr ) then
   call check( nf90_get_var(ncid, ncv, vx_s_g_conv) )
else
   go to 110
end if

if ( nf90_inq_varid(ncid, 'vy_s_g', ncv) == nf90_noerr ) then
   call check( nf90_get_var(ncid, ncv, vy_s_g_conv) )
else
   go to 110
end if

if ( nf90_inq_varid(ncid, 'vz_s', ncv) == nf90_noerr ) then
   call check( nf90_get_var(ncid, ncv, vz_s_conv) )
else
   go to 110
end if

if ( nf90_inq_varid(ncid, 'vh_s', ncv) == nf90_noerr ) then
   call check( nf90_get_var(ncid, ncv, vh_s_conv) )
else
   go to 110
end if

if ( nf90_inq_varid(ncid, 'temp_b', ncv) == nf90_noerr ) then
   call check( nf90_get_var(ncid, ncv, temp_b_conv) )
else
   go to 110
end if

if ( nf90_inq_varid(ncid, 'temph_b', ncv) == nf90_noerr ) then
   call check( nf90_get_var(ncid, ncv, temph_b_conv) )
else
   go to 110
end if

if ( nf90_inq_varid(ncid, 'p_b_w', ncv) == nf90_noerr ) then
   call check( nf90_get_var(ncid, ncv, p_b_w_conv) )
else
   go to 110
end if

if ( nf90_inq_varid(ncid, 'H_w', ncv) == nf90_noerr ) then
   call check( nf90_get_var(ncid, ncv, h_w_conv) )
else
   go to 110
end if

if ( nf90_inq_varid(ncid, 'q_gl_g', ncv) == nf90_noerr ) then
   call check( nf90_get_var(ncid, ncv, q_gl_g_conv) )
else
   go to 110
end if

go to 120

110 continue
write(6,'(/1x,a)') 'topography3_nc: End-of-file condition while reading *.nc.'
write(6, '(1x,a)') '                Some variables will be undefined.'

120 continue

call check( nf90_close(ncid) )

!-------- Convert data to real*8 and years to seconds --------

z_sl = real(z_sl_conv,dp)

h_r  = real(h_r_conv,dp)

do i=0, imax
   xi(i)  = real(xi_conv(i),dp)
end do

do j=0, jmax
   eta(j) = real(eta_conv(j),dp)
end do

do i=0, imax
do j=0, jmax

   maske(j,i)   = maske_conv(i,j)
   n_cts(j,i)   = n_cts_conv(i,j)
   kc_cts(j,i)  = kc_cts_conv(i,j)
   zs(j,i)      = real(zs_conv(i,j),dp)
   zm(j,i)      = real(zm_conv(i,j),dp)
   zb(j,i)      = real(zb_conv(i,j),dp)
   zl(j,i)      = real(zl_conv(i,j),dp)
#if (CALCMOD==1)
   h_c(j,i)     = real(H_cold_conv(i,j),dp)
   h_t(j,i)     = real(H_temp_conv(i,j),dp)
#elif (CALCMOD==0 || CALCMOD==2 || CALCMOD==3 || CALCMOD==-1)
   h_c(j,i)     = real(H_conv(i,j),dp)
   h_t(j,i)     = 0.0_dp
#endif
   q_bm(j,i)    = real(Q_bm_conv(i,j),dp)/year_sec
   q_tld(j,i)   = real(Q_tld_conv(i,j),dp)/year_sec
   am_perp(j,i) = real(am_perp_conv(i,j),dp)/year_sec
   qx(j,i)      = real(qx_conv(i,j),dp)/year_sec
   qy(j,i)      = real(qy_conv(i,j),dp)/year_sec
   dzs_dtau(j,i)  = real(dzs_dtau_conv(i,j),dp)/year_sec
   dzm_dtau(j,i)  = real(dzm_dtau_conv(i,j),dp)/year_sec
   dzb_dtau(j,i)  = real(dzb_dtau_conv(i,j),dp)/year_sec
   dzl_dtau(j,i)  = real(dzl_dtau_conv(i,j),dp)/year_sec
   dh_c_dtau(j,i) = real(dH_c_dtau_conv(i,j),dp)/year_sec
   dh_t_dtau(j,i) = real(dH_t_dtau_conv(i,j),dp)/year_sec
   vx_b_g(j,i)  = real(vx_b_g_conv(i,j),dp)/year_sec
   vy_b_g(j,i)  = real(vy_b_g_conv(i,j),dp)/year_sec
   vz_b(j,i)    = real(vz_b_conv(i,j),dp)/year_sec
   vx_s_g(j,i)  = real(vx_s_g_conv(i,j),dp)/year_sec
   vy_s_g(j,i)  = real(vy_s_g_conv(i,j),dp)/year_sec
   vz_s(j,i)    = real(vz_s_conv(i,j),dp)/year_sec
   temp_b(j,i)  = real(temp_b_conv(i,j),dp)
   temph_b(j,i) = real(temph_b_conv(i,j),dp)
   p_b_w(j,i)   = real(p_b_w_conv(i,j),dp)
   h_w(j,i)     = real(H_w_conv(i,j),dp)
   q_gl_g(j,i)  = real(q_gl_g_conv(i,j),dp)/year_sec

   do kr=0, krmax
      temp_r(kr,j,i) = real(temp_r_conv(i,j,kr),dp)
   end do

   do kt=0, ktmax
      vx_t(kt,j,i)    = real(vx_t_conv(i,j,kt),dp)/year_sec
      vy_t(kt,j,i)    = real(vy_t_conv(i,j,kt),dp)/year_sec
      vz_t(kt,j,i)    = real(vz_t_conv(i,j,kt),dp)/year_sec
      omega_t(kt,j,i) = real(omega_t_conv(i,j,kt),dp)
      age_t(kt,j,i)   = real(age_t_conv(i,j,kt),dp)*year_sec
      enth_t(kt,j,i)  = real(enth_t_conv(i,j,kt),dp)
      enh_t(kt,j,i)   = real(enh_t_conv(i,j,kt),dp)
   end do

   do kc=0, kcmax
      vx_c(kc,j,i)    = real(vx_c_conv(i,j,kc),dp)/year_sec
      vy_c(kc,j,i)    = real(vy_c_conv(i,j,kc),dp)/year_sec
      vz_c(kc,j,i)    = real(vz_c_conv(i,j,kc),dp)/year_sec
      temp_c(kc,j,i)  = real(temp_c_conv(i,j,kc),dp)
      age_c(kc,j,i)   = real(age_c_conv(i,j,kc),dp)*year_sec
      enth_c(kc,j,i)  = real(enth_c_conv(i,j,kc),dp)
      omega_c(kc,j,i) = real(omega_c_conv(i,j,kc),dp)
      enh_c(kc,j,i)   = real(enh_c_conv(i,j,kc),dp)
   end do

end do
end do

!-------- Read topography of the relaxed bedrock --------

open(23, iostat=ios, &
     file=inpath//'/'//trim(ch_domain_short)//'/'//zl0_file, &
     recl=16384, status='old')

if (ios /= 0) stop ' topography3_nc: Error when opening the zl0 file!'

do n=1, 6; read(23,'(a)') ch_dummy; end do

do j=jmax, 0, -1
   read(23,*) (zl0(j,i), i=0,imax)
end do

close(23, status='keep')

!-------- Further stuff --------

dxi  = dx *1000.0_dp   ! km -> m
deta = dx *1000.0_dp   ! km -> m

do i=0, imax
do j=0, jmax
   call geo_coord(phi(j,i), lambda(j,i), xi(i), eta(j))
end do
end do

!-------- Metric tensor, gradients of the topography --------

call metric()

#if TOPOGRAD==0
call topograd_1(dxi, deta, 1)
#elif TOPOGRAD==1
call topograd_2(dxi, deta, 1)
#endif

!-------- Corresponding area of grid points --------

do i=0, imax
do j=0, jmax
   area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
end do
end do

end subroutine topography3_nc
!