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, 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 (len=256) :: ch_attr_title, ch_attr_institution, ch_attr_source, &
                       ch_attr_history, ch_attr_references
character :: ch_dummy

!-------- 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) lambda_conv
read(11) phi_conv
read(11) temp_s_conv
read(11) as_perp_conv
read(11) maske
read(11) n_cts
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) vx_s_g_conv
read(11,end=110) vy_s_g_conv
read(11,end=110) vz_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,'(/1x,a)') 'topography3: End-of-file condition while reading *.erg.'
write(6, '(1x,a)') '             Some variables will be undefined.'

120 continue

close(11,status='keep')

!-------- 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//'/smars/'//zl0_file, &
     recl=8192, status='old')

if (ios /= 0) stop ' topography3: 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
   zl0(j,i) = zl0(j,i)*1000.0_dp   ! km --> m
   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
!