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-2013 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 netcdf

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
real(sp) :: time_conv, dummy_conv, z_sl_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), &
            lambda_conv(0:jmax,0:imax), phi_conv(0:jmax,0:imax), &
            temp_s_conv(0:jmax,0:imax), as_perp_conv(0:jmax,0:imax), &
            zs_conv(0:jmax,0:imax), zm_conv(0:jmax,0:imax), &
            zb_conv(0:jmax,0:imax), zl_conv(0:jmax,0:imax), &
            h_c_conv(0:jmax,0:imax), h_t_conv(0:jmax,0:imax), &
            h_conv(0:jmax,0:imax), &
            q_bm_conv(0:jmax,0:imax), q_tld_conv(0:jmax,0:imax), &
            am_perp_conv(0:jmax,0:imax), &
            qx_conv(0:jmax,0:imax), qy_conv(0:jmax,0:imax), &
            temp_r_conv(0:krmax,0:jmax,0:imax), &
            vx_c_conv(0:kcmax,0:jmax,0:imax), &
            vy_c_conv(0:kcmax,0:jmax,0:imax), &
            vz_c_conv(0:kcmax,0:jmax,0:imax), &
            vx_t_conv(0:ktmax,0:jmax,0:imax), &
            vy_t_conv(0:ktmax,0:jmax,0:imax), &
            vz_t_conv(0:ktmax,0:jmax,0:imax), &
            temp_c_conv(0:kcmax,0:jmax,0:imax), &
            age_c_conv(0:kcmax,0:jmax,0:imax), &
            omega_t_conv(0:ktmax,0:jmax,0:imax), &
            age_t_conv(0:ktmax,0:jmax,0:imax), &
            dzs_dtau_conv(0:jmax,0:imax), dzm_dtau_conv(0:jmax,0:imax), &
            dzb_dtau_conv(0:jmax,0:imax), dzl_dtau_conv(0:jmax,0:imax), &
            dh_c_dtau_conv(0:jmax,0:imax), dh_t_dtau_conv(0:jmax,0:imax), &
            dh_dtau_conv(0:jmax,0:imax), &
            vx_b_g_conv(0:jmax,0:imax), vy_b_g_conv(0:jmax,0:imax), &
            vz_b_conv(0:jmax,0:imax), &
            vx_s_g_conv(0:jmax,0:imax), vy_s_g_conv(0:jmax,0:imax), &
            vz_s_conv(0:jmax,0:imax), &
            temp_b_conv(0:jmax,0:imax), temph_b_conv(0:jmax,0:imax), &
            p_b_w_conv(0:jmax,0:imax), h_w_conv(0:jmax,0:imax), &
            q_gl_g_conv(0:jmax,0:imax)
character :: ch_dummy

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

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

call check( nf90_open(anfdatpath//'/'//trim(anfdatname), 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, '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) )

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

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_c', ncv) )
call check( nf90_get_var(ncid, ncv, h_c_conv) )

call check( nf90_inq_varid(ncid, 'H_t', ncv) )
call check( nf90_get_var(ncid, ncv, h_t_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) )

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, '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, '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 double precision --------

z_sl = real(z_sl_conv,dp)

h_r  = real(h_r_conv,dp)

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

   xi(i)        = real(xi_conv(i),dp)
   eta(j)       = real(eta_conv(j),dp)
   zs(j,i)      = real(zs_conv(j,i),dp)
   zm(j,i)      = real(zm_conv(j,i),dp)
   zb(j,i)      = real(zb_conv(j,i),dp)
   zl(j,i)      = real(zl_conv(j,i),dp)
   h_c(j,i)     = real(H_c_conv(j,i),dp)
   h_t(j,i)     = real(H_t_conv(j,i),dp)
   q_bm(j,i)    = real(Q_bm_conv(j,i),dp)
   q_tld(j,i)   = real(Q_tld_conv(j,i),dp)
   am_perp(j,i) = real(am_perp_conv(j,i),dp)
   qx(j,i)      = real(qx_conv(j,i),dp)
   qy(j,i)      = real(qy_conv(j,i),dp)
   dzs_dtau(j,i)  = real(dzs_dtau_conv(j,i),dp)
   dzm_dtau(j,i)  = real(dzm_dtau_conv(j,i),dp)
   dzb_dtau(j,i)  = real(dzb_dtau_conv(j,i),dp)
   dzl_dtau(j,i)  = real(dzl_dtau_conv(j,i),dp)
   dh_c_dtau(j,i) = real(dH_c_dtau_conv(j,i),dp)
   dh_t_dtau(j,i) = real(dH_t_dtau_conv(j,i),dp)
   vx_b_g(j,i)  = real(vx_b_g_conv(j,i),dp)
   vy_b_g(j,i)  = real(vy_b_g_conv(j,i),dp)
   vz_b(j,i)    = real(vz_b_conv(j,i),dp)
   vx_s_g(j,i)  = real(vx_s_g_conv(j,i),dp)
   vy_s_g(j,i)  = real(vy_s_g_conv(j,i),dp)
   vz_s(j,i)    = real(vz_s_conv(j,i),dp)
   temp_b(j,i)  = real(temp_b_conv(j,i),dp)
   temph_b(j,i) = real(temph_b_conv(j,i),dp)
   p_b_w(j,i)   = real(p_b_w_conv(j,i),dp)
   h_w(j,i)     = real(H_w_conv(j,i),dp)
   q_gl_g(j,i)  = real(q_gl_g_conv(j,i),dp)

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

   do kt=0, ktmax
      vx_t(kt,j,i)    = real(vx_t_conv(kt,j,i),dp)
      vy_t(kt,j,i)    = real(vy_t_conv(kt,j,i),dp)
      vz_t(kt,j,i)    = real(vz_t_conv(kt,j,i),dp)
      omega_t(kt,j,i) = real(omega_t_conv(kt,j,i),dp)
      age_t(kt,j,i)   = real(age_t_conv(kt,j,i),dp)
   end do

   do kc=0, kcmax
      vx_c(kc,j,i)    = real(vx_c_conv(kc,j,i),dp)
      vy_c(kc,j,i)    = real(vy_c_conv(kc,j,i),dp)
      vz_c(kc,j,i)    = real(vz_c_conv(kc,j,i),dp)
      temp_c(kc,j,i)  = real(temp_c_conv(kc,j,i),dp)
      age_c(kc,j,i)   = real(age_c_conv(kc,j,i),dp)
   end do

end do
end do

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

open(23, iostat=ios, &
     file=inpath//'/tibet/'//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  = dlambda *pi_180   ! deg -> rad
deta = dphi    *pi_180   ! deg -> rad

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

contains

!-------------------------------------------------------------------------------
!> NetCDF error capturing.
!<------------------------------------------------------------------------------
   subroutine check(status)
      integer(i4b), intent (in) :: status
      if(status /= nf90_noerr) then 
         write(6,'(a)') trim(nf90_strerror(status))
         stop ' topography3_nc: Stopped due to netcdf error!'
      end if
   end subroutine check  

end subroutine topography3_nc
!