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-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/>.
!<
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!-------------------------------------------------------------------------------
!> Main loop of SICOPOLIS.
!<------------------------------------------------------------------------------
subroutine sico_main_loop(delta_ts, glac_index, &
                    mean_accum, mean_accum_inv, &
                    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

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, mean_accum_inv
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, mean_accum_inv)   ! polythermal mode

#elif CALCMOD==0
   call calc_temp_cold(dxi, deta, dzeta_c, dzeta_t, dzeta_r, &
                       dtime_temp, mean_accum_inv)   ! cold-ice mode
#endif

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

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

!  ------ Flow enhancement factor

   call calc_enhance(time)

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

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

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

#if MARGIN==3
call calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, 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

!-------- Topography --------

call calc_top(time, time_init, &
              z_sl, dzsl_dtau, z_mar, &
              dtime, dxi, deta, dzeta_t, &
              mean_accum, &
              ii, jj, nn, itercount, iter_wss)

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

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