sico_main_loop_m.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Module :  s i c o _ m a i n _ l o o p _ m
!
!> @file
!!
!! Main loop of SICOPOLIS.
!!
!! @section Copyright
!!
!! Copyright 2009-2017 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.
!<------------------------------------------------------------------------------
module sico_main_loop_m

  use sico_types_m
  use sico_variables_m
  use sico_vars_m

  implicit none

  public

contains

!-------------------------------------------------------------------------------
!> Main routine of sico_main_loop_m: 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 boundary_m

#if (CALCMOD==0 || CALCMOD==1 || CALCMOD==-1)
  use calc_temp_m
#elif (CALCMOD==2 || CALCMOD==3)
  use calc_temp_enth_m
#endif

  use calc_enhance_m
  use calc_vxy_m
  use calc_vz_m
  use calc_dxyz_m
  use calc_gia_m
  use calc_thk_m
  use calc_temp_melt_bas_m
  use calc_bas_melt_m
  use calc_thk_water_bas_m
  use output_m

  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/)') ' -------- sico_main_loop --------'

  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='(2x,i0)') itercount

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

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

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

  maske_old = 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='(10x,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 between -1 and 3!'
#endif

!  ------ Time derivative of H_t (further time derivatives are
!         computed in subroutine calc_thk_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

#if (ENHMOD==1)
     call calc_enhance_1()
#elif (ENHMOD==2)
     call calc_enhance_2()
#elif (ENHMOD==3)
     call calc_enhance_3(time)
#elif (ENHMOD==4)
     call calc_enhance_4()
#elif (ENHMOD==5)
     call calc_enhance_5()
#else
     stop ' >>> sico_main_loop: Parameter ENHMOD must be between 1 and 5!'
#endif

  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_grounded(dxi, deta, dzeta_c, dzeta_t)

#if (MARGIN==3)
  call calc_vz_floating(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)

  call calc_thk_init()

#if ((MARGIN==3 || DYNAMICS==2) && (CALCTHK==1 || CALCTHK==2 || CALCTHK==3))
    write(6, fmt='(a)') ' >>> sico_main_loop:'
    write(6, fmt='(a)') '          Non-SIA dynamics combined with'
    write(6, fmt='(a)') '          SIA ice thickness solver!'
    stop
#endif

#if (CALCTHK==1)
    call calc_thk_sia_expl(time, dtime, dxi, deta, z_mar)
#elif (CALCTHK==2 || CALCTHK==3)
    call calc_thk_sia_impl(time, dtime, dxi, deta, z_mar, mean_accum, ii, jj, nn)
#elif (CALCTHK==4)
    call calc_thk_expl(time, dtime, dxi, deta, z_mar)
#elif (CALCTHK==5 || CALCTHK==6)
    call calc_thk_impl(time, dtime, dxi, deta, z_mar, mean_accum, ii, jj, nn)
#else
    stop ' >>> sico_main_loop: Parameter CALCTHK must be between 1 and 6!'
#endif

#if (MARGIN==3)       /* coupled SIA/SSA or SIA/SStA/SSA dynamics */
    call calc_thk_mask_update(time, dtime, dxi, deta, z_sl, z_mar, 3_i2b)
#elif (DYNAMICS==2)   /* hybrid SIA/SStA dynamics */
    call calc_thk_mask_update(time, dtime, dxi, deta, z_sl, z_mar, 2_i2b)
#else                 /* SIA-only dynamics */
#if (CALCTHK==1 || CALCTHK==2 || CALCTHK==3)
    call calc_thk_mask_update(time, dtime, dxi, deta, z_sl, z_mar, 1_i2b)
#elif (CALCTHK==4 || CALCTHK==5 || CALCTHK==6)
    call calc_thk_mask_update(time, dtime, dxi, deta, z_sl, z_mar, 2_i2b)
#endif
#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_thk_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

     call output1(runname, time, delta_ts, glac_index, z_sl, &
                  flag_3d_output, ndat2d, ndat3d)

  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

        call output1(runname, time, delta_ts, glac_index, z_sl, &
                     flag_3d_output, ndat2d, ndat3d)

     end if

  end do

#elif (OUTPUT==3)

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

     flag_3d_output = .false.

     call output1(runname, time, delta_ts, glac_index, z_sl, &
                  flag_3d_output, ndat2d, ndat3d)

  end if

  do n=1, n_output

     if (itercount == iter_output(n)) then

        flag_3d_output = .true.

        call output1(runname, time, delta_ts, glac_index, z_sl, &
                     flag_3d_output, ndat2d, ndat3d)

     end if

  end do

#endif

  if ( mod(itercount, iter_ser) == 0 ) then
     call output2(time, dxi, deta, delta_ts, glac_index, z_sl)
     call output4(time, dxi, deta, delta_ts, glac_index, z_sl)
#if (defined(ASF) && WRITE_SER_FILE_STAKES>0)
     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

!-------------------------------------------------------------------------------

end module sico_main_loop_m
!