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-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-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
use sico_vars

#if (CALCMOD==1 || CALCMOD==0 || CALCMOD==-1)
use enth_temp_omega, only : enth_fct_temp_omega
#endif

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), dimension(0:JMAX,0:IMAX) :: h, h_cold, h_temp, dh_dtau
real(dp) :: v_tot, a_grounded, a_floating
character(len=256) :: filename
character :: ch_1000s, ch_100s, ch_10s, ch_1s

integer(i2b), dimension(0:IMAX,0:JMAX) :: maske_conv, n_cts_conv, kc_cts_conv
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)
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_cold_conv, h_temp_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, &
                                              enth_c_conv, omega_c_conv, &
                                              enh_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, &
                                              enth_t_conv, &
                                              enh_t_conv
real(sp), dimension(0:IMAX,0:JMAX,0:KRMAX) :: temp_r_conv

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(unit=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(unit=11) ch_attr_title

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

ch_attr_source = 'SICOPOLIS Version '//version
write(unit=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(unit=11) ch_attr_history

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

!-------- Ice thickness and time derivative --------

h       = h_c       + h_t
dh_dtau = dh_c_dtau + dh_t_dtau

!-------- 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

!-------- Enthalpies for the non-enthalpy methods (POLY, COLD, ISOT) --------

#if (CALCMOD==1)

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

   do kc=0, kcmax
      enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
   end do

   if ( (maske(j,i) == 0_i2b).and.(n_cts(j,i) == 1_i2b) ) then
      do kt=0, ktmax
         enth_t(kt,j,i) = enth_fct_temp_omega(temp_t_m(kt,j,i), omega_t(kt,j,i))
      end do
   else
      do kt=0, ktmax
         enth_t(kt,j,i) = enth_c(0,j,i)
      end do
   end if

end do
end do

#elif (CALCMOD==0 || CALCMOD==-1)

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

   do kc=0, kcmax
      enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
   end do

   do kt=0, ktmax
      enth_t(kt,j,i) = enth_c(0,j,i)
   end do

end do
end do

#endif

!-------- 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(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 and seconds to years --------

#if (!defined(OUT_TIMES) || OUT_TIMES==1)
time_conv = real(time/year_sec,sp)
#elif (OUT_TIMES==2)
time_conv = real((time+year_zero)/year_sec,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

   maske_conv(i,j)     = maske(j,i)
   n_cts_conv(i,j)     = n_cts(j,i)
   kc_cts_conv(i,j)    = kc_cts(j,i)

   lambda_conv(i,j)    = real(lambda(j,i),sp) 
   phi_conv(i,j)       = real(phi(j,i),sp) 
   lon_conv(i,j)       = real(lambda(j,i)*pi_180_inv,sp)   ! longitude in deg
   lon_conv(i,j)       = modulo(lon_conv(i,j)+180.0_sp, 360.0_sp)-180.0_sp
                                    ! mapping to interval [-180 deg, +180 deg)
   lat_conv(i,j)       = real(phi(j,i)   *pi_180_inv,sp)   ! latitute  in deg
   if (lat_conv(i,j) >  90.0_sp) lat_conv(i,j) =  90.0_sp
   if (lat_conv(i,j) < -90.0_sp) lat_conv(i,j) = -90.0_sp
                                 ! constraining to interval [-90 deg, +90 deg]
   temp_s_conv(i,j)    = real(temp_s(j,i),sp)
   as_perp_conv(i,j)   = real(as_perp(j,i)*year_sec,sp)
   zs_conv(i,j)        = real(zs(j,i),sp)
   zm_conv(i,j)        = real(zm(j,i),sp)
   zb_conv(i,j)        = real(zb(j,i),sp)
   zl_conv(i,j)        = real(zl(j,i),sp)
   h_cold_conv(i,j)    = real(H_cold(j,i),sp)
   h_temp_conv(i,j)    = real(H_temp(j,i),sp)
   h_conv(i,j)         = real(H(j,i),sp)
   q_bm_conv(i,j)      = real(q_bm(j,i)*year_sec,sp)
   q_tld_conv(i,j)     = real(q_tld(j,i)*year_sec,sp)
   am_perp_conv(i,j)   = real(am_perp(j,i)*year_sec,sp)
   qx_conv(i,j)        = real(qx(j,i)*year_sec,sp)
   qy_conv(i,j)        = real(qy(j,i)*year_sec,sp)
   dzs_dtau_conv(i,j)  = real(dzs_dtau(j,i)*year_sec,sp)
   dzm_dtau_conv(i,j)  = real(dzm_dtau(j,i)*year_sec,sp)
   dzb_dtau_conv(i,j)  = real(dzb_dtau(j,i)*year_sec,sp)
   dzl_dtau_conv(i,j)  = real(dzl_dtau(j,i)*year_sec,sp)
   dh_c_dtau_conv(i,j) = real(dh_c_dtau(j,i)*year_sec,sp)
   dh_t_dtau_conv(i,j) = real(dh_t_dtau(j,i)*year_sec,sp)
   dh_dtau_conv(i,j)   = real(dh_dtau(j,i)*year_sec,sp)
   vx_b_g_conv(i,j)    = real(vx_b_g(j,i)*year_sec,sp)
   vy_b_g_conv(i,j)    = real(vy_b_g(j,i)*year_sec,sp)
   vz_b_conv(i,j)      = real(vz_b(j,i)*year_sec,sp)
   vh_b_conv(i,j)      = sqrt( vx_b_g_conv(i,j)**2 + vy_b_g_conv(i,j)**2 )
   vx_s_g_conv(i,j)    = real(vx_s_g(j,i)*year_sec,sp)
   vy_s_g_conv(i,j)    = real(vy_s_g(j,i)*year_sec,sp)
   vz_s_conv(i,j)      = real(vz_s(j,i)*year_sec,sp)
   vh_s_conv(i,j)      = sqrt( vx_s_g_conv(i,j)**2 + vy_s_g_conv(i,j)**2 )
   temp_b_conv(i,j)    = real(temp_b(j,i),sp)
   temph_b_conv(i,j)   = real(temph_b(j,i),sp)
   p_b_w_conv(i,j)     = real(p_b_w(j,i),sp)
   h_w_conv(i,j)       = real(H_w(j,i),sp)
   q_gl_g_conv(i,j)    = real(q_gl_g(j,i)*year_sec,sp)

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

   do kt=0, ktmax
      vx_t_conv(i,j,kt)    = real(vx_t(kt,j,i)*year_sec,sp)
      vy_t_conv(i,j,kt)    = real(vy_t(kt,j,i)*year_sec,sp)
      vz_t_conv(i,j,kt)    = real(vz_t(kt,j,i)*year_sec,sp)
      omega_t_conv(i,j,kt) = real(omega_t(kt,j,i),sp)
      age_t_conv(i,j,kt)   = real(age_t(kt,j,i)/year_sec,sp)
      enth_t_conv(i,j,kt)  = real(enth_t(kt,j,i),sp)
      enh_t_conv(i,j,kt)   = real(enh_t(kt,j,i),sp)
   end do

   do kc=0, kcmax
      vx_c_conv(i,j,kc)    = real(vx_c(kc,j,i)*year_sec,sp)
      vy_c_conv(i,j,kc)    = real(vy_c(kc,j,i)*year_sec,sp)
      vz_c_conv(i,j,kc)    = real(vz_c(kc,j,i)*year_sec,sp)
      temp_c_conv(i,j,kc)  = real(temp_c(kc,j,i),sp)
      age_c_conv(i,j,kc)   = real(age_c(kc,j,i)/year_sec,sp)
      enth_c_conv(i,j,kc)  = real(enth_c(kc,j,i),sp)
      omega_c_conv(i,j,kc) = real(omega_c(kc,j,i),sp)
      enh_c_conv(i,j,kc)   = real(enh_c(kc,j,i),sp)
   end do

end do
end do

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

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

close(unit=11, status='keep')

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

ndat = ndat+1

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

end subroutine output1
!