output2.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Subroutine :  o u t p u t 2
!
!> @file
!!
!! Writing of time-series data on file in ASCII 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/>.
!<
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!-------------------------------------------------------------------------------
!> Writing of time-series data on file in ASCII format.
!<------------------------------------------------------------------------------
subroutine output2(time, dxi, deta, delta_ts, glac_index, z_sl)

use sico_types
use sico_variables
use sico_vars

implicit none

real(dp), intent(in) :: time, dxi, deta, delta_ts, glac_index, z_sl

integer(i4b) :: i, j
integer(i4b) :: n_base, n_tempbase
real(dp) :: time_val, &
            v_tot, v_grounded, v_floating, &
            a_tot, a_grounded, a_floating, &
            h_max, zs_max, v_temp, v_fw, v_sle, &
            a_temp, v_bm, v_tld, h_t_max, vs_max, tbh_max
real(dp) :: x_pos, y_pos
real(dp), dimension(0:JMAX,0:IMAX) :: h, h_cold, h_temp
real(dp) :: v_gr_redu, a_surf, rhosw_rho_ratio
real(dp) :: vs_help, tbh_help
real(dp) :: h_ave_sed, tbh_ave_sed, atb_sed
real(dp) :: sum_area_sed

#if (!defined(OUTSER_V_A) || OUTSER_V_A==1)

character(len=128), parameter :: &
                    fmt1 = '(1pe13.6,2(1pe13.4),/,13x,6(1pe13.4),/,'// &
                               '13x,6(1pe13.4),/)'
#elif (OUTSER_V_A==2)

character(len=128), parameter :: &
                    fmt1 = '(1pe13.6,2(1pe13.4),/,13x,6(1pe13.4),/,'// &
                           '26x,5(1pe13.4),/,26x,5(1pe13.4),/)'
#else
   stop ' sico_init: OUTSER_V_A must be either 1 or 2!'
#endif

character(len=128), parameter :: &
                    fmt2 = '(1pe13.6,2(1pe13.4),/,13x,3(1pe13.4),/)'

!-------- Ice thickness --------

h = h_c + h_t

!-------- Thickness of the cold and temperate layers --------

h_cold = 0.0_dp
h_temp = 0.0_dp

#if (CALCMOD==1)
do i=1, imax-1
do j=1, jmax-1
   h_temp(j,i) = h_t(j,i)
end do
end do
#elif (CALCMOD==0 || CALCMOD==2 || CALCMOD==3 || CALCMOD==-1)
do i=1, imax-1
do j=1, jmax-1
   h_temp(j,i) = h_c(j,i)*eaz_c_quotient(kc_cts(j,i))
end do
end do
#endif

h_cold = h - h_temp

!-------- Maximum ice elevation and thickness, ice volume,
!         sea-level equivalent, volume of the temperate ice,
!         area, area covered by temperate ice,
!         freshwater production due to melting and calving,
!         water drainage due to basal melting,
!         water drainage from the temperate layer,
!         maximum thickness of the temperate layer,
!         maximum surface velocity,
!         maximum basal temperature rel. to pmp --------

#if (defined(ANT) \
      || defined(grl) \
      || defined(scand) \
      || defined(nhem) \
      || defined(tibet) \
      || defined(asf) \
      || defined(emtp2sge) \
      || defined(xyz))   /* terrestrial ice sheet */

rhosw_rho_ratio = rho_sw/rho

#elif (defined(NMARS) || defined(SMARS))   /* Martian ice sheet */

rhosw_rho_ratio = 0.0_dp   ! dummy value

#else

stop ' output2: No valid domain specified!'

#endif

v_grounded = 0.0_dp
v_gr_redu  = 0.0_dp
v_floating = 0.0_dp
a_grounded = 0.0_dp
a_floating = 0.0_dp
h_max      = 0.0_dp
zs_max     = -1.11e+11_dp
v_temp     = 0.0_dp
v_fw       = 0.0_dp
a_temp     = 0.0_dp
v_bm       = 0.0_dp
v_tld      = 0.0_dp
h_t_max    = 0.0_dp
vs_max     = 0.0_dp
tbh_max    = -273.15_dp

do i=1, imax-1
do j=1, jmax-1

   if (maske(j,i) == 0) then   ! grounded ice

      if (zs(j,i)       > zs_max)  zs_max  = zs(j,i)
      if (h(j,i)        > h_max  ) h_max   = h(j,i)
      if (h_temp(j,i)   > h_t_max) h_t_max = h_temp(j,i)

      v_grounded = v_grounded + h(j,i)     *area(j,i)
      v_temp     = v_temp     + h_temp(j,i)*area(j,i)

#if (defined(ANT) \
      || defined(grl) \
      || defined(scand) \
      || defined(nhem) \
      || defined(tibet) \
      || defined(asf) \
      || defined(emtp2sge) \
      || defined(xyz))   /* terrestrial ice sheet */

      if (zl0(j,i) < z_sl) &
         v_gr_redu = v_gr_redu + rhosw_rho_ratio*(z_sl-zl0(j,i))*area(j,i)

#endif

      a_grounded = a_grounded + area(j,i)

      if (n_cts(j,i) /= -1) a_temp = a_temp + area(j,i)

      v_bm  = v_bm  + q_bm(j,i)  *area(j,i)
      v_tld = v_tld + q_tld(j,i) *area(j,i)

      vs_help = sqrt( &
               0.25_dp*(vx_c(kcmax,j,i)+vx_c(kcmax,j,i-1))**2 &
              +0.25_dp*(vy_c(kcmax,j,i)+vy_c(kcmax,j-1,i))**2 )
      if (vs_help > vs_max) vs_max = vs_help

      if (n_cts(j,i) >= 0) then   ! temperate base
         tbh_max = 0.0_dp
      else   ! cold base
         tbh_help = min((temp_c(0,j,i)-temp_c_m(0,j,i)), -eps)
         if (tbh_help > tbh_max) tbh_max = tbh_help
      end if

   else if (maske(j,i) == 3) then   ! floating ice
                                    ! (basal temperature assumed to be below
                                    !  the pressure melting point for pure ice)

      if (zs(j,i)    > zs_max) zs_max = zs(j,i)
      if (h(j,i)     > h_max)  h_max  = h(j,i)

      v_floating = v_floating + h(j,i)*area(j,i)
      a_floating = a_floating + area(j,i)

      v_bm       = v_bm       + q_bm(j,i) *area(j,i)

      vs_help = sqrt( &
               0.25_dp*(vx_c(kcmax,j,i)+vx_c(kcmax,j,i-1))**2 &
              +0.25_dp*(vy_c(kcmax,j,i)+vy_c(kcmax,j-1,i))**2 )
      if (vs_help > vs_max) vs_max = vs_help

      tbh_help = min((temp_c(0,j,i)-temp_c_m(0,j,i)), -eps)
      if (tbh_help > tbh_max) tbh_max = tbh_help

   end if

   if (     (maske(j,i) == 0).or.(maske_help(j,i) == 0) &
        .or.(maske(j,i) == 3).or.(maske_help(j,i) == 3) ) then
                     ! grounded or floating ice before or after the time step
      v_fw  = v_fw + accum(j,i)*area(j,i) &
                   - (dzs_dtau(j,i)-dzb_dtau(j,i))*area(j,i)
   end if

end do
end do

!-------- Conversion --------

#if (defined(ANT) \
      || defined(grl) \
      || defined(scand) \
      || defined(nhem) \
      || defined(tibet) \
      || defined(asf) \
      || defined(emtp2sge) \
      || defined(xyz))   /* terrestrial ice sheet */

a_surf = 3.61132e+14_dp              ! global ocean area, in m^2
v_sle  = (v_grounded-v_gr_redu)*(rho/rho_sw)/a_surf
                                     ! m^3 ice equiv./m^2 -> m sea water equiv.

#elif (defined(NMARS) || defined(SMARS))   /* Martian ice sheet */

a_surf = 1.4441e+14_dp            ! surface area of Mars, in m^2
v_sle  = v_grounded*(rho_i/rho_w)*(1.0_dp-frac_dust)/a_surf
                                  ! m^3 (ice+dust) equiv./m^2 -> m water equiv.

#endif

#if (!defined(OUT_TIMES) || OUT_TIMES==1)
time_val = time /year_sec   ! s -> a
#elif (OUT_TIMES==2)
time_val = (time+year_zero) /year_sec   ! s -> a
#else
stop ' output2: OUT_TIMES must be either 1 or 2!'
#endif

v_tot = v_grounded + v_floating   ! in m^3
a_tot = a_grounded + a_floating   ! in m^2

if (zs_max < -1.1e+11_dp) zs_max = 0.0_dp

vs_max = vs_max *year_sec              ! m/s -> m/a

#if (defined(ANT) \
      || defined(grl) \
      || defined(scand) \
      || defined(nhem) \
      || defined(tibet) \
      || defined(asf) \
      || defined(emtp2sge) \
      || defined(xyz))   /* terrestrial ice sheet */

v_fw  = v_fw  *year_sec *(rho/rho_w)
              ! m^3/s ice equiv. -> m^3/a water equiv.
v_bm  = v_bm  *year_sec *(rho/rho_w)
              ! m^3/s ice equiv. -> m^3/a water equiv.
v_tld = v_tld *year_sec *(rho/rho_w)
              ! m^3/s ice equiv. -> m^3/a water equiv.

#elif (defined(NMARS) || defined(SMARS))   /* Martian ice sheet */

v_fw  = v_fw  *year_sec *(rho_i/rho_w)*(1.0_dp-frac_dust)
              ! m^3/s (ice+dust) equiv. -> m^3/a water equiv.
v_bm  = v_bm  *year_sec *(rho_i/rho_w)*(1.0_dp-frac_dust)
              ! m^3/s (ice+dust) equiv. -> m^3/a water equiv.
v_tld = v_tld *year_sec *(rho_i/rho_w)*(1.0_dp-frac_dust)
              ! m^3/s (ice+dust) equiv. -> m^3/a water equiv.
#endif

!-------- Writing of data on time-series file --------

#if (defined(ANT) \
      || defined(grl) \
      || defined(scand) \
      || defined(nhem) \
      || defined(tibet) \
      || defined(asf) \
      || defined(emtp2sge) \
      || defined(xyz))   /* terrestrial ice sheet */

if (forcing_flag == 1) then

#if (!defined(OUTSER_V_A) || OUTSER_V_A==1)

   write(unit=12, fmt=trim(fmt1)) time_val, delta_ts, z_sl, &
      h_max, zs_max, v_grounded, v_temp, v_fw, v_sle, &
      a_grounded, a_temp, v_bm, v_tld, h_t_max, vs_max

#elif (OUTSER_V_A==2)

   write(unit=12, fmt=trim(fmt1)) time_val, delta_ts, z_sl, &
      v_tot, v_grounded, v_floating, a_tot, a_grounded, a_floating, &
      h_max, zs_max, v_temp, v_fw, v_sle, &
      a_temp, v_bm, v_tld, h_t_max, vs_max

#else
   stop ' sico_init: OUTSER_V_A must be either 1 or 2!'
#endif

else if (forcing_flag == 2) then

#if (!defined(OUTSER_V_A) || OUTSER_V_A==1)

   write(unit=12, fmt=trim(fmt1)) time_val, glac_index, z_sl, &
      h_max, zs_max, v_grounded, v_temp, v_fw, v_sle, &
      a_grounded, a_temp, v_bm, v_tld, h_t_max, vs_max

#elif (OUTSER_V_A==2)

   write(unit=12, fmt=trim(fmt1)) time_val, glac_index, z_sl, &
      v_tot, v_grounded, v_floating, a_tot, a_grounded, a_floating, &
      h_max, zs_max, v_temp, v_fw, v_sle, &
      a_temp, v_bm, v_tld, h_t_max, vs_max

#else
   stop ' sico_init: OUTSER_V_A must be either 1 or 2!'
#endif

else if (forcing_flag == 3) then

#if (!defined(OUTSER_V_A) || OUTSER_V_A==1)

   write(unit=12, fmt=trim(fmt1)) time_val, 1.11e+11_dp, z_sl, &
      h_max, zs_max, v_grounded, v_temp, v_fw, v_sle, &
      a_grounded, a_temp, v_bm, v_tld, h_t_max, vs_max

#elif (OUTSER_V_A==2)

   write(unit=12, fmt=trim(fmt1)) time_val, 1.11e+11_dp, z_sl, &
      v_tot, v_grounded, v_floating, a_tot, a_grounded, a_floating, &
      h_max, zs_max, v_temp, v_fw, v_sle, &
      a_temp, v_bm, v_tld, h_t_max, vs_max

#else
   stop ' sico_init: OUTSER_V_A must be either 1 or 2!'
#endif

end if

#elif (defined(NMARS) || defined(SMARS))   /* Martian ice sheet */

if (forcing_flag == 1) then

   write(unit=12, fmt=trim(fmt1)) time_val, delta_ts, z_sl, &
      h_max, zs_max, v_tot, v_temp, v_fw, tbh_max, &
      a_tot, a_temp, v_bm, v_tld, h_t_max, vs_max

else if (forcing_flag == 2) then

   write(unit=12, fmt=trim(fmt1)) time_val, glac_index, z_sl, &
      h_max, zs_max, v_tot, v_temp, v_fw, tbh_max, &
      a_tot, a_temp, v_bm, v_tld, h_t_max, vs_max

else if (forcing_flag == 3) then

   write(unit=12, fmt=trim(fmt1)) time_val, 1.11e+11_dp, z_sl, &
      h_max, zs_max, v_tot, v_temp, v_fw, tbh_max, &
      a_tot, a_temp, v_bm, v_tld, h_t_max, vs_max

end if

#endif

!-------- Special output for various domains --------

#if (defined(XYZ))

!  ------ Time-series data for the sediment region of HEINO

#if (defined(HEINO))

!    ---- Average ice thickness, average basal temperature rel. to pmp,
!         temperate basal area

h_ave_sed    = 0.0_dp
tbh_ave_sed  = 0.0_dp
atb_sed      = 0.0_dp
sum_area_sed = 0.0_dp

do i=1, imax-1
do j=1, jmax-1

   if (maske_sedi(j,i) == 7) then   ! sediment region

      sum_area_sed = sum_area_sed + area(j,i)

      h_ave_sed = h_ave_sed + area(j,i)*h(j,i)

      if (n_cts(j,i) /= -1) then   ! temperate base
         tbh_help    = 0.0_dp
      else   ! cold base
         tbh_help    = min((temp_c(0,j,i)-temp_c_m(0,j,i)), -eps)
      end if
      tbh_ave_sed = tbh_ave_sed + area(j,i)*tbh_help

      if (n_cts(j,i) /= -1) atb_sed = atb_sed + area(j,i)

   end if

end do
end do

if (sum_area_sed > eps) then
   h_ave_sed   = h_ave_sed   / sum_area_sed
   tbh_ave_sed = tbh_ave_sed / sum_area_sed
else
   stop ' output2: No sediment area found!'
end if

!    ---- Writing of data on time-series file

if (forcing_flag == 1) then
   write(unit=15, fmt=trim(fmt2)) time_val, delta_ts, z_sl, &
                                  h_ave_sed, tbh_ave_sed, atb_sed
else if (forcing_flag == 2) then
   write(unit=15, fmt=trim(fmt2)) time_val, glac_index, z_sl, &
                                  h_ave_sed, tbh_ave_sed, atb_sed
else if (forcing_flag == 3) then
   write(unit=15, fmt=trim(fmt2)) time_val, 1.11e+11_dp, z_sl, &
                                  h_ave_sed, tbh_ave_sed, atb_sed
end if

#endif

#endif

end subroutine output2
!