sico_main_loop.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Subroutine :  s i c o _ m a i n _ l o o p
!
!> @file
!!
!! Main loop of SICOPOLIS.
!!
!! @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/>.
!<
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!-------------------------------------------------------------------------------
!> Main loop of SICOPOLIS.
!<------------------------------------------------------------------------------
subroutine sico_main_loop(delta_ts, glac_index, &
                    mean_accum, &
                    dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, &
                    time, time_init, time_end, time_output, &
                    dxi, deta, dzeta_c, dzeta_t, dzeta_r, &
                    z_sl, dzsl_dtau, z_mar, &
                    ii, jj, nn, &
                    ndat2d, ndat3d, n_output, &
                    runname)

use sico_types
use sico_variables
use sico_vars

implicit none

integer(i4b),       intent(in)    :: ii((imax+1)*(jmax+1)), &
                                     jj((imax+1)*(jmax+1)), &
                                     nn(0:jmax,0:imax)
integer(i4b),       intent(in)    :: n_output
real(dp),           intent(in)    :: mean_accum
real(dp),           intent(in)    :: dtime, dtime_temp, dtime_wss, &
                                     dtime_out, dtime_ser
real(dp),           intent(in)    :: time_init, time_end, time_output(100)
real(dp),           intent(in)    :: dxi, deta, dzeta_c, dzeta_t, dzeta_r
character(len=100), intent(in)    :: runname

integer(i4b),       intent(inout) :: ndat2d, ndat3d
real(dp),           intent(inout) :: delta_ts, glac_index
real(dp),           intent(inout) :: time
real(dp),           intent(inout) :: z_sl, dzsl_dtau, z_mar

integer(i4b) :: i, j, kc, kt, kr, n
integer(i4b) :: itercount, itercount_max
integer(i4b) :: iter_temp, iter_wss, iter_ser, iter_out, iter_output(100)
real(dp)     :: dtime_temp_inv
logical      :: flag_3d_output

!-------- Begin of main loop --------

write(unit=6, fmt='(a)') ' '

itercount_max = nint((time_end-time_init)/dtime)

iter_temp = nint(dtime_temp/dtime)
#if (REBOUND==2)
iter_wss = nint(dtime_wss/dtime)
#endif
iter_ser  = nint(dtime_ser/dtime)
#if (OUTPUT==1 || OUTPUT==3)
iter_out  = nint(dtime_out/dtime)
#endif
#if (OUTPUT==2 || OUTPUT==3)
do n=1, n_output
   iter_output(n) = nint((time_output(n)-time_init)/dtime)
end do
#endif

main_loop : do itercount=1, itercount_max

write(unit=6, fmt='(i7)') itercount

!-------- Update of time --------

time      = time_init + real(itercount,dp)*dtime

!-------- Save old mask --------

maske_help = maske

!-------- Boundary conditions --------

call boundary(time, dtime, dxi, deta, delta_ts, glac_index, &
              z_sl, dzsl_dtau, z_mar)

!-------- Temperature, water content, age, flow enhancement factor --------

if ( mod(itercount, iter_temp) == 0 ) then

   write(unit=6, fmt='(a)') ' -------- Computation of T'

!  ------ Temperature, water content, age

#if (CALCMOD==1)
   call calc_temp_poly(dxi, deta, dzeta_c, dzeta_t, dzeta_r, dtime_temp)
                  ! polythermal scheme (POLY)

#elif (CALCMOD==0)
   call calc_temp_cold(dxi, deta, dzeta_c, dzeta_t, dzeta_r, dtime_temp)
                  ! cold-ice scheme (COLD)

#elif (CALCMOD==2 || CALCMOD==3)
   call calc_temp_enth(dxi, deta, dzeta_c, dzeta_t, dzeta_r, dtime_temp)
                  ! enthalpy scheme (ENTC or ENTM)

#elif (CALCMOD==-1)
   call calc_temp_const()   ! isothermal scheme (ISOT)

#else
   stop ' sico_main_loop: Parameter CALCMOD must be either -1, 0, 1, 2 or 3!'
#endif

!  ------ Time derivative of H_t (further time derivatives are
!         computed in subroutine calc_top_xxx)

   dtime_temp_inv = 1.0_dp/dtime_temp

   dh_t_dtau      = (h_t_neu - h_t)*dtime_temp_inv

!  ------ New values --> old values

   n_cts   = n_cts_neu
   kc_cts  = kc_cts_neu
   zm      = zm_neu
   h_c     = h_c_neu
   h_t     = h_t_neu
   temp_c  = temp_c_neu
   age_c   = age_c_neu
   omega_t = omega_t_neu
   age_t   = age_t_neu
   temp_r  = temp_r_neu

#if (CALCMOD==2 || CALCMOD==3)
   enth_c  = enth_c_neu
   enth_t  = enth_t_neu
   omega_c = omega_c_neu
#endif

!  ------ Flow enhancement factor

   call calc_enhance(time)

end if
!     End of computation of temperature, water content, age and
!     enhancement factor

!-------- Velocity --------

#if (DYNAMICS==1 || DYNAMICS==2)

call calc_vxy_b_sia(time, z_sl)
call calc_vxy_sia(dzeta_c, dzeta_t)

#if (MARGIN==3 || DYNAMICS==2)
call calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)
#endif

call calc_vz_sia(dxi, deta, dzeta_c, dzeta_t)

#if (MARGIN==3)
call calc_vz_ssa(z_sl, dxi, deta, dzeta_c)
#endif

#elif (DYNAMICS==0)

call calc_vxy_static()
call calc_vz_static()

#else
   stop ' sico_main_loop: DYNAMICS must be either 0, 1 or 2!'
#endif

call calc_dxyz(dxi, deta, dzeta_c, dzeta_t)

!-------- Glacial isostatic adjustment and ice topography --------

call calc_gia(time, z_sl, dtime, dxi, deta, itercount, iter_wss)

#if (MARGIN==3)       /* coupled SIA/SSA or SIA/SStA/SSA dynamics */
call calc_top_3(time, z_sl, z_mar, dtime, dxi, deta, mean_accum, ii, jj, nn)
#elif (DYNAMICS==2)   /* hybrid SIA/SStA dynamics */
call calc_top_2(time, z_sl, z_mar, dtime, dxi, deta, mean_accum, ii, jj, nn)
#else                 /* SIA-only dynamics */
call calc_top_1(time, z_sl, z_mar, dtime, dxi, deta, mean_accum, ii, jj, nn)
#endif

!  ------ New values --> old values

zs = zs_neu
zm = zm_neu
zb = zb_neu
zl = zl_neu
h_c= h_c_neu
h_t= h_t_neu

!-------- Melting temperature --------

call calc_temp_melt()

!-------- Basal temperature --------

call calc_temp_bas()

!-------- Basal melting rate --------

call calc_qbm(time, z_sl, dzeta_c, dzeta_r)

!-------- Effective thickness of subglacial water  --------

call calc_water_bas(z_sl)

!-------- Data output --------

#if (OUTPUT==1)

if ( mod(itercount, iter_out) == 0 ) then

#if (ERGDAT==0)
   flag_3d_output = .false.
#elif (ERGDAT==1)
   flag_3d_output = .true.
#endif

#if (NETCDF==1)
   call output1(runname, time, delta_ts, glac_index, z_sl, &
                flag_3d_output, ndat2d, ndat3d)
#elif (NETCDF==2)
   call output_nc(runname, time, delta_ts, glac_index, z_sl, &
                  flag_3d_output, ndat2d, ndat3d)
#endif

end if

#elif (OUTPUT==2)

do n=1, n_output

   if (itercount == iter_output(n)) then

#if (ERGDAT==0)
      flag_3d_output = .false.
#elif (ERGDAT==1)
      flag_3d_output = .true.
#endif

#if (NETCDF==1)
      call output1(runname, time, delta_ts, glac_index, z_sl, &
                   flag_3d_output, ndat2d, ndat3d)
#elif (NETCDF==2)
      call output_nc(runname, time, delta_ts, glac_index, z_sl, &
                     flag_3d_output, ndat2d, ndat3d)
#endif

   end if

end do

#elif (OUTPUT==3)

if ( mod(itercount, iter_out) == 0 ) then

   flag_3d_output = .false.

#if (NETCDF==1)
   call output1(runname, time, delta_ts, glac_index, z_sl, &
                flag_3d_output, ndat2d, ndat3d)
#elif (NETCDF==2)
   call output_nc(runname, time, delta_ts, glac_index, z_sl, &
                  flag_3d_output, ndat2d, ndat3d)
#endif

end if

do n=1, n_output

   if (itercount == iter_output(n)) then

      flag_3d_output = .true.

#if (NETCDF==1)
      call output1(runname, time, delta_ts, glac_index, z_sl, &
                   flag_3d_output, ndat2d, ndat3d)
#elif (NETCDF==2)
      call output_nc(runname, time, delta_ts, glac_index, z_sl, &
                     flag_3d_output, ndat2d, ndat3d)
#endif

   end if

end do

#endif

if ( mod(itercount, iter_ser) == 0 ) then
   call output2(time, dxi, deta, delta_ts, glac_index, z_sl)
#if (OUTSER==2)
   call output3(time, delta_ts, glac_index, z_sl)
#endif
#if (OUTSER==3)
   call output4(time, dxi, deta, delta_ts, glac_index, z_sl)
#endif
#if (defined(ASF) && OUTSER==4)
   call output5(time, dxi, deta, delta_ts, glac_index, z_sl)
#endif

end if

end do main_loop   ! End of main loop

end subroutine sico_main_loop
!