topography3.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Subroutine :  t o p o g r a p h y 3
!
!> @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 binary 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 binary format).
!<------------------------------------------------------------------------------
subroutine topography3(dxi, deta, z_sl, anfdatname)

use sico_types
use sico_variables

implicit none

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

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

integer(i4b) :: i, j, kc, kt, kr
integer(i4b) :: ios
character(len=256) :: ch_attr_title, ch_attr_institution, ch_attr_source, &
                      ch_attr_history, ch_attr_references

integer(i2b), dimension(0:IMAX,0:JMAX) :: maske_conv, n_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_c_conv, h_t_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
real(sp), dimension(0:IMAX,0:JMAX,0:KTMAX) :: vx_t_conv, vy_t_conv, vz_t_conv, &
                                              omega_t_conv, age_t_conv
real(sp), dimension(0:IMAX,0:JMAX,0:KRMAX) :: temp_r_conv

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

open(11, iostat=ios, &
         file=anfdatpath//'/'//trim(anfdatname), &
         status='old', form='unformatted')
if (ios /= 0) stop ' topography3: Error when opening the initial-value file!'

read(11) ch_attr_title
read(11) ch_attr_institution
read(11) ch_attr_source
read(11) ch_attr_history
read(11) ch_attr_references
read(11) time_conv
read(11) dummy_conv
read(11) z_sl_conv
read(11) xi_conv
read(11) eta_conv
read(11) sigma_level_c_conv
read(11) sigma_level_t_conv
read(11) sigma_level_r_conv
read(11) lon_conv
read(11) lat_conv
read(11) lambda_conv
read(11) phi_conv
read(11) temp_s_conv
read(11) as_perp_conv
read(11) maske_conv
read(11) n_cts_conv
read(11) zs_conv
read(11) zm_conv
read(11) zb_conv
read(11) zl_conv
read(11) h_c_conv
read(11) h_t_conv
read(11) h_conv
read(11) h_r_conv
read(11) vx_c_conv
read(11) vy_c_conv
read(11) vz_c_conv
read(11) vx_t_conv
read(11) vy_t_conv
read(11) vz_t_conv
read(11) temp_c_conv
read(11) omega_t_conv
read(11) temp_r_conv
read(11) q_bm_conv
read(11) q_tld_conv
read(11) am_perp_conv
read(11) qx_conv
read(11) qy_conv
read(11) age_c_conv
read(11) age_t_conv
read(11) dzs_dtau_conv
read(11) dzm_dtau_conv
read(11) dzb_dtau_conv
read(11) dzl_dtau_conv
read(11) dh_c_dtau_conv
read(11) dh_t_dtau_conv
read(11) dh_dtau_conv
read(11,end=110) vx_b_g_conv
read(11,end=110) vy_b_g_conv
read(11,end=110) vz_b_conv
read(11,end=110) vh_b_conv
read(11,end=110) vx_s_g_conv
read(11,end=110) vy_s_g_conv
read(11,end=110) vz_s_conv
read(11,end=110) vh_s_conv
read(11,end=110) temp_b_conv
read(11,end=110) temph_b_conv
read(11,end=110) p_b_w_conv
read(11,end=110) h_w_conv
read(11,end=110) q_gl_g_conv

go to 120

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

120 continue

close(11,status='keep')

!-------- 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)
   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)
   h_c(j,i)     = real(H_c_conv(i,j),dp)
   h_t(j,i)     = real(H_t_conv(i,j),dp)
   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
   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
   end do

end do
end do

!-------- Set topography of the relaxed bedrock --------

zl0 = 0.0_dp

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

#if (GRID==0 || GRID==1)

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

#elif GRID==2

stop ' topography3: GRID==2 not allowed for this application!'

#endif

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
!