output1.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Subroutine :  o u t p u t 1
!
!> @file
!!
!! Writing of time-slice files 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/>.
!<
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!-------------------------------------------------------------------------------
!> Writing of time-slice files in binary format.
!<------------------------------------------------------------------------------
subroutine output1(runname, time, delta_ts, glac_index, z_sl, &
                   flag_3d_output, ndat2d, ndat3d)

use sico_types
use sico_variables

implicit none

real(dp), intent (in) :: time, delta_ts, glac_index, z_sl
character (len=100), intent (in) :: runname
logical, intent (in) :: flag_3d_output

integer(i4b), intent (inout) :: ndat2d, ndat3d

integer(i4b) :: i, j, kc, kt, kr
integer(i4b) :: ndat, ndat_help, ndat_1000s, ndat_100s, ndat_10s, ndat_1s
real(dp) :: v_tot, a_grounded, a_floating
real(sp) :: time_conv, delta_ts_conv, glac_index_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), &
            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) :: filename
character :: ch_1000s, ch_100s, ch_10s, ch_1s

integer(i4b) :: ios
character (len= 16) :: ch_date, ch_time, ch_zone
character (len=256) :: ch_attr_title, ch_attr_institution, ch_attr_source, &
                       ch_attr_history, ch_attr_references

!-------- Create consecutively numbered file names --------

if (flag_3d_output) then
   ndat = ndat3d
else
   ndat = ndat2d
end if

if (ndat > 9999) stop ' output1: Too many time-slice files!'

ndat_help  = ndat
ndat_1000s = ndat_help/1000
ndat_help  = ndat_help-ndat_1000s*1000
ndat_100s  = ndat_help/100
ndat_help  = ndat_help-ndat_100s*100
ndat_10s   = ndat_help/10
ndat_help  = ndat_help-ndat_10s*10
ndat_1s    = ndat_help

ch_1000s = char(ndat_1000s+ichar('0'))
ch_100s  = char(ndat_100s +ichar('0'))
ch_10s   = char(ndat_10s  +ichar('0'))
ch_1s    = char(ndat_1s   +ichar('0'))

if (flag_3d_output) then
   filename = trim(runname)//ch_1000s//ch_100s//ch_10s//ch_1s//'.erg'
else
   filename = trim(runname)//'_2d_'//ch_1000s//ch_100s//ch_10s//ch_1s//'.erg'
end if

!-------- Open time-slice file runname.erg --------

open(11, iostat=ios, &
         file=outpath//'/'//trim(filename), &
         status='new', form='unformatted')
if (ios /= 0) stop ' output1: Error when opening an erg file!'

!-------- Write attributes on file --------

ch_attr_title = 'Time-slice output no. '//ch_1000s//ch_100s//ch_10s//ch_1s// &
                'of simulation '//trim(runname)
write(11) ch_attr_title

ch_attr_institution = 'Institute of Low Temperature Science, '// &
                      'Hokkaido University, Sapporo, Japan'
write(11) ch_attr_institution

ch_attr_source = 'SICOPOLIS Version '//version
write(11) ch_attr_source

call date_and_time(ch_date, ch_time, ch_zone)
ch_attr_history = ch_date(1:4)//'-'//ch_date(5:6)//'-'//ch_date(7:8)//' '// &
                  ch_time(1:2)//':'//ch_time(3:4)//':'//ch_time(5:6)//' '// &
                  ch_zone(1:3)//':'//ch_zone(4:5)//' - Data produced'
write(11) ch_attr_history

ch_attr_references = 'http://sicopolis.greveweb.net/'
write(11) ch_attr_references

!-------- Computation of the ice volume,
!         the area covered by grounded ice
!         and the area covered by floating ice --------

v_tot      = 0.0_dp
a_grounded = 0.0_dp
a_floating = 0.0_dp

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

   if ( (maske(j,i)==0).or.(maske(j,i)==3) ) &   ! grounded or floating ice
      v_tot  = v_tot + (h_c(j,i)+h_t(j,i))*area(j,i)

   if (maske(j,i)==0) &   ! grounded ice
      a_grounded = a_grounded + area(j,i)

   if (maske(j,i)==3) &   ! floating ice
      a_floating = a_floating + area(j,i)

end do
end do

!-------- Convert data to real*4 --------

#if ( !defined(OUT_TIMES) || OUT_TIMES==1 )
time_conv = real(time,sp)
#elif OUT_TIMES == 2
time_conv = real((time+year_zero),sp)
#else
stop ' output1: OUT_TIMES must be either 1 or 2!'
#endif

delta_ts_conv   = real(delta_ts,sp)
glac_index_conv = real(glac_index,sp)
z_sl_conv       = real(z_sl,sp)
v_tot_conv      = real(v_tot,sp)
a_grounded_conv = real(a_grounded,sp)
a_floating_conv = real(a_floating,sp)
h_r_conv        = real(h_r,sp)

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

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

do kc=0, kcmax
   sigma_level_c_conv(kc) = real(eaz_c_quotient(kc),sp)
end do

do kt=0, ktmax
   sigma_level_t_conv(kt) = real(zeta_t(kt),sp)
end do

do kr=0, krmax
   sigma_level_r_conv(kr) = real(kr,sp)/real(krmax,sp)
end do

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

   lambda_conv(j,i)    = real(lambda(j,i),sp)
   phi_conv(j,i)       = real(phi(j,i),sp)
   temp_s_conv(j,i)    = real(temp_s(j,i),sp)
   as_perp_conv(j,i)   = real(as_perp(j,i),sp)
   zs_conv(j,i)        = real(zs(j,i),sp)
   zm_conv(j,i)        = real(zm(j,i),sp)
   zb_conv(j,i)        = real(zb(j,i),sp)
   zl_conv(j,i)        = real(zl(j,i),sp)
   h_c_conv(j,i)       = real(H_c(j,i),sp)
   h_t_conv(j,i)       = real(H_t(j,i),sp)
   h_conv(j,i)         = real(H_c(j,i)+H_t(j,i),sp)
   q_bm_conv(j,i)      = real(Q_bm(j,i),sp)
   q_tld_conv(j,i)     = real(Q_tld(j,i),sp)
   am_perp_conv(j,i)   = real(am_perp(j,i),sp)
   qx_conv(j,i)        = real(qx(j,i),sp)
   qy_conv(j,i)        = real(qy(j,i),sp)
   dzs_dtau_conv(j,i)  = real(dzs_dtau(j,i),sp)
   dzm_dtau_conv(j,i)  = real(dzm_dtau(j,i),sp)
   dzb_dtau_conv(j,i)  = real(dzb_dtau(j,i),sp)
   dzl_dtau_conv(j,i)  = real(dzl_dtau(j,i),sp)
   dh_c_dtau_conv(j,i) = real(dH_c_dtau(j,i),sp)
   dh_t_dtau_conv(j,i) = real(dH_t_dtau(j,i),sp)
   dh_dtau_conv(j,i)   = real(dH_c_dtau(j,i)+dH_t_dtau(j,i),sp)
   vx_b_g_conv(j,i)    = real(vx_b_g(j,i),sp)
   vy_b_g_conv(j,i)    = real(vy_b_g(j,i),sp)
   vz_b_conv(j,i)      = real(vz_b(j,i),sp)
   vx_s_g_conv(j,i)    = real(vx_s_g(j,i),sp)
   vy_s_g_conv(j,i)    = real(vy_s_g(j,i),sp)
   vz_s_conv(j,i)      = real(vz_s(j,i),sp)
   temp_b_conv(j,i)    = real(temp_b(j,i),sp)
   temph_b_conv(j,i)   = real(temph_b(j,i),sp)
   p_b_w_conv(j,i)     = real(p_b_w(j,i),sp)
   h_w_conv(j,i)       = real(H_w(j,i),sp)
   q_gl_g_conv(j,i)    = real(q_gl_g(j,i),sp)

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

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

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

end do
end do

!-------- Write data on file --------

write(11) time_conv
if (forcing_flag == 1) then
   write(11) delta_ts_conv
else if (forcing_flag == 2) then
   write(11) glac_index_conv
else if (forcing_flag == 3) then
   glac_index_conv = 1.11e+11   ! dummy value
   write(11) glac_index_conv
end if
write(11) z_sl_conv
write(11) xi_conv
write(11) eta_conv
write(11) sigma_level_c_conv
write(11) sigma_level_t_conv
write(11) sigma_level_r_conv
write(11) lambda_conv
write(11) phi_conv
write(11) temp_s_conv
write(11) as_perp_conv
write(11) maske
write(11) n_cts
write(11) zs_conv
write(11) zm_conv
write(11) zb_conv
write(11) zl_conv
write(11) h_c_conv
write(11) h_t_conv
write(11) h_conv
write(11) h_r_conv
if (flag_3d_output) then
   write(11) vx_c_conv
   write(11) vy_c_conv
   write(11) vz_c_conv
   write(11) vx_t_conv
   write(11) vy_t_conv
   write(11) vz_t_conv
   write(11) temp_c_conv
   write(11) omega_t_conv
   write(11) temp_r_conv
end if
write(11) q_bm_conv
write(11) q_tld_conv
write(11) am_perp_conv
write(11) qx_conv
write(11) qy_conv
if (flag_3d_output) then
   write(11) age_c_conv
   write(11) age_t_conv
end if
write(11) dzs_dtau_conv
write(11) dzm_dtau_conv
write(11) dzb_dtau_conv
write(11) dzl_dtau_conv
write(11) dh_c_dtau_conv
write(11) dh_t_dtau_conv
write(11) dh_dtau_conv
write(11) vx_b_g_conv
write(11) vy_b_g_conv
write(11) vz_b_conv
write(11) vx_s_g_conv
write(11) vy_s_g_conv
write(11) vz_s_conv
write(11) temp_b_conv
write(11) temph_b_conv
write(11) p_b_w_conv
write(11) h_w_conv
write(11) q_gl_g_conv
write(11) v_tot_conv
write(11) a_grounded_conv
write(11) a_floating_conv

close(11, status='keep')

!-------- Increase file counter --------

ndat = ndat+1

if (flag_3d_output) then
   ndat3d = ndat
else
   ndat2d = ndat
end if

end subroutine output1
!