43 real(dp),
dimension(0:JMAX,0:IMAX),
save ::
h,
h_neu 63 stop
' >>> calc_thk_init: Replace CALCZS by CALCTHK in header file!' 64 #elif (!defined(CALCTHK)) 65 stop
' >>> calc_thk_init: Define CALCTHK in header file!' 71 #if (MARGIN==1 || MARGIN==2) /* only grounded ice */ 76 #elif (MARGIN==3) /* grounded and floating ice */ 78 where (
maske <= 1_i2b)
87 stop
' >>> calc_thk_init: MARGIN must be either 1, 2 or 3!' 100 #if (CALCTHK==1 || CALCTHK==4) /* explicit solver */ 104 #elif (CALCTHK==2 || CALCTHK==3 \ 105 || calcthk==5 || calcthk==6) /*
implicit solver */
110 stop
' >>> calc_thk_init: CALCTHK must be between 1 and 6!' 119 #elif (MB_ACCOUNT==1) 128 stop
' >>> calc_thk_init: MB_ACCOUNT must be either 0 or 1!' 140 real(dp),
intent(in) :: time, dtime, dxi, deta
141 real(dp),
intent(in) :: z_mar
144 real(dp) :: azs2, azs3
145 real(dp),
dimension(0:JMAX,0:IMAX) :: czs2, czs3
149 azs2 = dtime/(dxi*dxi)
150 azs3 = dtime/(deta*deta)
175 #if (MARGIN==1 && MB_ACCOUNT==0) 176 .and.(
maske(j,i) <= 1_i2b) &
182 + ( czs2(j,i) *(
zs(j,i+1)-
zs(j,i) ) &
183 -czs2(j,i-1)*(
zs(j,i) -
zs(j,i-1)) &
184 +czs3(j,i) *(
zs(j+1,i)-
zs(j,i) ) &
185 -czs3(j-1,i)*(
zs(j,i) -
zs(j-1,i)) ) &
207 call thk_adjust(time, dtime)
221 integer(i4b),
intent(in) :: ii((imax+1)*(jmax+1)), jj((imax+1)*(jmax+1))
222 integer(i4b),
intent(in) :: nn(0:jmax,0:imax)
223 real(dp),
intent(in) :: time, dtime, dxi, deta
224 real(dp),
intent(in) :: z_mar
225 real(dp),
intent(in) :: mean_accum
228 integer(i4b) :: k, nnz
229 real(dp) :: azs2, azs3
230 real(dp),
dimension(0:JMAX,0:IMAX) :: czs2, czs3
236 integer(i4b) :: nc, nr
237 integer(i4b),
parameter :: nmax = (imax+1)*(jmax+1)
238 integer(i4b),
parameter :: n_sprs = 10*(imax+1)*(jmax+1)
239 integer(i4b),
allocatable,
dimension(:) :: lgs_a_ptr, lgs_a_index
240 integer(i4b),
allocatable,
dimension(:) :: lgs_a_diag_index
241 integer(i4b),
allocatable,
dimension(:) :: lgs_a_index_trim
242 real(dp),
allocatable,
dimension(:) :: lgs_a_value, lgs_b_value, lgs_x_value
243 real(dp),
allocatable,
dimension(:) :: lgs_a_value_trim
250 lis_integer :: nc, nr
251 lis_integer,
parameter :: nmax = (imax+1)*(jmax+1)
252 lis_integer,
parameter :: n_sprs = 10*(imax+1)*(jmax+1)
253 lis_integer,
allocatable,
dimension(:) :: lgs_a_ptr, lgs_a_index
254 lis_integer,
allocatable,
dimension(:) :: lgs_a_diag_index
256 lis_vector :: lgs_b, lgs_x
257 lis_scalar,
allocatable,
dimension(:) :: lgs_a_value, lgs_b_value, lgs_x_value
259 character(len=256) :: ch_solver_set_option
265 azs2 = dtime/(dxi*dxi)
266 azs3 = dtime/(deta*deta)
288 #if (CALCTHK==2 || CALCTHK==3) 290 allocate(lgs_a_value(n_sprs), lgs_a_index(n_sprs), lgs_a_ptr(nmax+1))
291 allocate(lgs_a_diag_index(nmax), lgs_b_value(nmax), lgs_x_value(nmax))
310 #if (MARGIN==1 && MB_ACCOUNT==0) 311 .and.(
maske(j,i) <= 1_i2b) &
318 lgs_a_value(k) = -czs3(j-1,i)*ovi_weight &
325 lgs_a_value(k) = -czs2(j,i-1)*ovi_weight &
334 lgs_a_value(k) = 1.0_dp &
335 + (czs2(j,i)+czs2(j,i-1)+czs3(j,i)+czs3(j-1,i)) &
338 lgs_a_diag_index(nr) = k
344 lgs_a_value(k) = -czs2(j,i)*ovi_weight &
351 lgs_a_value(k) = -czs3(j,i)*ovi_weight &
355 lgs_b_value(nr) =
zs(j,i) &
357 + ( czs2(j,i)*(
zs(j,i+1)-
zs(j,i)) &
358 -czs2(j,i-1)*(
zs(j,i)-
zs(j,i-1)) &
359 +czs3(j,i)*(
zs(j+1,i)-
zs(j,i)) &
360 -czs3(j-1,i)*(
zs(j,i)-
zs(j-1,i)) ) &
361 *(1.0_dp-ovi_weight) &
369 lgs_a_value(k) = 1.0_dp
370 lgs_a_diag_index(nr) = k
372 lgs_b_value(nr) =
zb_neu(j,i)
376 lgs_x_value(nr) =
zs(j,i)
379 lgs_a_ptr(nr+1) = k+1
386 stop
' >>> calc_thk_sia_impl: CALCTHK must be either 2 or 3!' 395 allocate(lgs_a_value_trim(nnz), lgs_a_index_trim(nnz))
398 lgs_a_value_trim(k) = lgs_a_value(k)
399 lgs_a_index_trim(k) = lgs_a_index(k)
402 deallocate(lgs_a_value, lgs_a_index)
404 eps_sor = 1.0e-05_dp*mean_accum*dtime
407 lgs_a_index_trim, lgs_a_diag_index, lgs_a_ptr, &
409 nnz, nmax, omega_sor, eps_sor, lgs_x_value, ierr)
414 zs_neu(j,i) = lgs_x_value(nr)
417 deallocate(lgs_a_value_trim, lgs_a_index_trim, lgs_a_ptr)
418 deallocate(lgs_a_diag_index, lgs_b_value, lgs_x_value)
424 call lis_matrix_create(lis_comm_world, lgs_a, ierr)
425 call lis_vector_create(lis_comm_world, lgs_b, ierr)
426 call lis_vector_create(lis_comm_world, lgs_x, ierr)
428 call lis_matrix_set_size(lgs_a, 0, nmax, ierr)
429 call lis_vector_set_size(lgs_b, 0, nmax, ierr)
430 call lis_vector_set_size(lgs_x, 0, nmax, ierr)
434 do nc=lgs_a_ptr(nr), lgs_a_ptr(nr+1)-1
435 call lis_matrix_set_value(lis_ins_value, nr, lgs_a_index(nc), &
436 lgs_a_value(nc), lgs_a, ierr)
439 call lis_vector_set_value(lis_ins_value, nr, lgs_b_value(nr), lgs_b, ierr)
440 call lis_vector_set_value(lis_ins_value, nr, lgs_x_value(nr), lgs_x, ierr)
444 call lis_matrix_set_type(lgs_a, lis_matrix_csr, ierr)
445 call lis_matrix_assemble(lgs_a, ierr)
449 call lis_solver_create(solver, ierr)
451 ch_solver_set_option =
'-i bicg -p ilu '// &
452 '-maxiter 1000 -tol 1.0e-04 -initx_zeros false' 454 call lis_solver_set_option(trim(ch_solver_set_option), solver, ierr)
457 call lis_solve(lgs_a, lgs_b, lgs_x, solver, ierr)
460 call lis_solver_get_iter(solver, iter, ierr)
461 write(6,
'(10x,a,i0)')
'calc_thk_sia_impl: iter = ', iter
464 call lis_vector_gather(lgs_x, lgs_x_value, ierr)
465 call lis_matrix_destroy(lgs_a, ierr)
466 call lis_vector_destroy(lgs_b, ierr)
467 call lis_vector_destroy(lgs_x, ierr)
468 call lis_solver_destroy(solver, ierr)
473 zs_neu(j,i) = lgs_x_value(nr)
476 deallocate(lgs_a_value, lgs_a_index, lgs_a_ptr)
477 deallocate(lgs_a_diag_index, lgs_b_value, lgs_x_value)
493 call thk_adjust(time, dtime)
504 real(dp),
intent(in) :: time, dtime, dxi, deta
505 real(dp),
intent(in) :: z_mar
508 real(dp),
dimension(0:JMAX,0:IMAX) :: dt_darea
509 real(dp),
dimension(0:JMAX,0:IMAX) :: vx_m_1, vx_m_2, vy_m_1, vy_m_2
510 real(dp),
dimension(0:JMAX,0:IMAX) :: upH_x_1, upH_x_2, upH_y_1, upH_y_2
511 real(dp),
dimension(0:JMAX,0:IMAX) :: sq_g22_x_1, sq_g22_x_2
512 real(dp),
dimension(0:JMAX,0:IMAX) :: sq_g11_y_1, sq_g11_y_2
516 dt_darea = dtime/
area 523 vx_m_1(j,i) =
vx_m(j,i-1)
524 vx_m_2(j,i) =
vx_m(j,i)
525 vy_m_1(j,i) =
vy_m(j-1,i)
526 vy_m_2(j,i) =
vy_m(j,i)
528 if (vx_m_1(j,i) >= 0.0_dp)
then 529 uph_x_1(j,i) =
h(j,i-1)
531 uph_x_1(j,i) =
h(j,i)
534 if (vx_m_2(j,i) >= 0.0_dp)
then 535 uph_x_2(j,i) =
h(j,i)
537 uph_x_2(j,i) =
h(j,i+1)
540 if (vy_m_1(j,i) >= 0.0_dp)
then 541 uph_y_1(j,i) =
h(j-1,i)
543 uph_y_1(j,i) =
h(j,i)
546 if (vy_m_2(j,i) >= 0.0_dp)
then 547 uph_y_2(j,i) =
h(j,i)
549 uph_y_2(j,i) =
h(j+1,i)
564 uph_x_1(j,i) = 0.0_dp
565 uph_x_2(j,i) = 0.0_dp
566 uph_y_1(j,i) = 0.0_dp
567 uph_y_2(j,i) = 0.0_dp
569 sq_g22_x_1(j,i) = 0.0_dp
570 sq_g22_x_2(j,i) = 0.0_dp
571 sq_g11_y_1(j,i) = 0.0_dp
572 sq_g11_y_2(j,i) = 0.0_dp
582 #if (MARGIN==1 && MB_ACCOUNT==0) 583 .and.(
maske <= 1_i2b) &
589 * ( ( vx_m_2*uph_x_2*sq_g22_x_2*deta &
590 -vx_m_1*uph_x_1*sq_g22_x_1*deta ) &
591 + ( vy_m_2*uph_y_2*sq_g11_y_2*dxi &
592 -vy_m_1*uph_y_1*sq_g11_y_1*dxi ) )
606 call thk_adjust(time, dtime)
613 subroutine calc_thk_impl(time, dtime, dxi, deta, z_mar, mean_accum, &
620 integer(i4b),
intent(in) :: ii((imax+1)*(jmax+1)), jj((imax+1)*(jmax+1))
621 integer(i4b),
intent(in) :: nn(0:jmax,0:imax)
622 real(dp),
intent(in) :: time, dtime, dxi, deta
623 real(dp),
intent(in) :: z_mar
624 real(dp),
intent(in) :: mean_accum
627 integer(i4b) :: k, nnz
628 real(dp),
dimension(0:JMAX,0:IMAX) :: dt_darea
629 real(dp),
dimension(0:JMAX,0:IMAX) :: vx_m_1, vx_m_2, vy_m_1, vy_m_2
630 real(dp),
dimension(0:JMAX,0:IMAX) :: upH_x_1, upH_x_2, upH_y_1, upH_y_2
631 real(dp),
dimension(0:JMAX,0:IMAX) :: sq_g22_x_1, sq_g22_x_2
632 real(dp),
dimension(0:JMAX,0:IMAX) :: sq_g11_y_1, sq_g11_y_2
638 integer(i4b) :: nc, nr
639 integer(i4b),
parameter :: nmax = (imax+1)*(jmax+1)
640 integer(i4b),
parameter :: n_sprs = 10*(imax+1)*(jmax+1)
641 integer(i4b),
allocatable,
dimension(:) :: lgs_a_ptr, lgs_a_index
642 integer(i4b),
allocatable,
dimension(:) :: lgs_a_diag_index
643 integer(i4b),
allocatable,
dimension(:) :: lgs_a_index_trim
644 real(dp),
allocatable,
dimension(:) :: lgs_a_value, lgs_b_value, lgs_x_value
645 real(dp),
allocatable,
dimension(:) :: lgs_a_value_trim
652 lis_integer :: nc, nr
653 lis_integer,
parameter :: nmax = (imax+1)*(jmax+1)
654 lis_integer,
parameter :: n_sprs = 10*(imax+1)*(jmax+1)
655 lis_integer,
allocatable,
dimension(:) :: lgs_a_ptr, lgs_a_index
656 lis_integer,
allocatable,
dimension(:) :: lgs_a_diag_index
658 lis_vector :: lgs_b, lgs_x
659 lis_scalar,
allocatable,
dimension(:) :: lgs_a_value, lgs_b_value, lgs_x_value
661 character(len=256) :: ch_solver_set_option
667 dt_darea = dtime/
area 674 vx_m_1(j,i) =
vx_m(j,i-1)
675 vx_m_2(j,i) =
vx_m(j,i)
676 vy_m_1(j,i) =
vy_m(j-1,i)
677 vy_m_2(j,i) =
vy_m(j,i)
679 if (vx_m_1(j,i) >= 0.0_dp)
then 680 uph_x_1(j,i) =
h(j,i-1)
682 uph_x_1(j,i) =
h(j,i)
685 if (vx_m_2(j,i) >= 0.0_dp)
then 686 uph_x_2(j,i) =
h(j,i)
688 uph_x_2(j,i) =
h(j,i+1)
691 if (vy_m_1(j,i) >= 0.0_dp)
then 692 uph_y_1(j,i) =
h(j-1,i)
694 uph_y_1(j,i) =
h(j,i)
697 if (vy_m_2(j,i) >= 0.0_dp)
then 698 uph_y_2(j,i) =
h(j,i)
700 uph_y_2(j,i) =
h(j+1,i)
715 uph_x_1(j,i) = 0.0_dp
716 uph_x_2(j,i) = 0.0_dp
717 uph_y_1(j,i) = 0.0_dp
718 uph_y_2(j,i) = 0.0_dp
720 sq_g22_x_1(j,i) = 0.0_dp
721 sq_g22_x_2(j,i) = 0.0_dp
722 sq_g11_y_1(j,i) = 0.0_dp
723 sq_g11_y_2(j,i) = 0.0_dp
733 #if (CALCTHK==5 || CALCTHK==6) 735 allocate(lgs_a_value(n_sprs), lgs_a_index(n_sprs), lgs_a_ptr(nmax+1))
736 allocate(lgs_a_diag_index(nmax), lgs_b_value(nmax), lgs_x_value(nmax))
755 #if (MARGIN==1 && MB_ACCOUNT==0) 756 .and.(
maske(j,i) <= 1_i2b) &
760 k=k+1 ; nc=nn(j-1,i) ; lgs_a_index(k)=nc
761 if (vy_m_1(j,i) > 0.0_dp) &
762 lgs_a_value(k) = -dt_darea(j,i)*vy_m_1(j,i) &
763 *sq_g11_y_1(j,i)*dxi*ovi_weight
765 k=k+1 ; nc=nn(j,i-1) ; lgs_a_index(k)=nc
766 if (vx_m_1(j,i) > 0.0_dp) &
767 lgs_a_value(k) = -dt_darea(j,i)*vx_m_1(j,i) &
768 *sq_g22_x_1(j,i)*deta*ovi_weight
770 k=k+1 ; lgs_a_index(k)=nr ; lgs_a_diag_index(nr)=k
771 lgs_a_value(k) = 1.0_dp
772 if (vy_m_1(j,i) < 0.0_dp) &
773 lgs_a_value(k) = lgs_a_value(k) &
774 - dt_darea(j,i)*vy_m_1(j,i) &
775 *sq_g11_y_1(j,i)*dxi*ovi_weight
776 if (vx_m_1(j,i) < 0.0_dp) &
777 lgs_a_value(k) = lgs_a_value(k) &
778 - dt_darea(j,i)*vx_m_1(j,i) &
779 *sq_g22_x_1(j,i)*deta*ovi_weight
780 if (vx_m_2(j,i) > 0.0_dp) &
781 lgs_a_value(k) = lgs_a_value(k) &
782 + dt_darea(j,i)*vx_m_2(j,i) &
783 *sq_g22_x_2(j,i)*deta*ovi_weight
784 if (vy_m_2(j,i) > 0.0_dp) &
785 lgs_a_value(k) = lgs_a_value(k) &
786 + dt_darea(j,i)*vy_m_2(j,i) &
787 *sq_g11_y_2(j,i)*dxi*ovi_weight
789 k=k+1 ; nc=nn(j,i+1) ; lgs_a_index(k)=nc
790 if (vx_m_2(j,i) < 0.0_dp) &
791 lgs_a_value(k) = dt_darea(j,i)*vx_m_2(j,i) &
792 *sq_g22_x_2(j,i)*deta*ovi_weight
794 k=k+1 ; nc=nn(j+1,i) ; lgs_a_index(k)=nc
795 if (vy_m_2(j,i) < 0.0_dp) &
796 lgs_a_value(k) = dt_darea(j,i)*vy_m_2(j,i) &
797 *sq_g11_y_2(j,i)*dxi*ovi_weight
799 lgs_b_value(nr) =
h(j,i) &
801 -(1.0_dp-ovi_weight) &
803 * ( ( vx_m_2(j,i)*uph_x_2(j,i) &
804 *sq_g22_x_2(j,i)*deta &
805 -vx_m_1(j,i)*uph_x_1(j,i) &
806 *sq_g22_x_1(j,i)*deta ) &
807 + ( vy_m_2(j,i)*uph_y_2(j,i) &
808 *sq_g11_y_2(j,i)*dxi &
809 -vy_m_1(j,i)*uph_y_1(j,i) &
810 *sq_g11_y_1(j,i)*dxi ) )
816 lgs_a_value(k) = 1.0_dp
817 lgs_a_diag_index(nr) = k
819 lgs_b_value(nr) = 0.0_dp
823 lgs_x_value(nr) =
h(j,i)
826 lgs_a_ptr(nr+1) = k+1
833 stop
' >>> calc_thk_impl: CALCTHK must be either 5 or 6!' 842 allocate(lgs_a_value_trim(nnz), lgs_a_index_trim(nnz))
845 lgs_a_value_trim(k) = lgs_a_value(k)
846 lgs_a_index_trim(k) = lgs_a_index(k)
849 deallocate(lgs_a_value, lgs_a_index)
851 eps_sor = 1.0e-05_dp*mean_accum*dtime
854 lgs_a_index_trim, lgs_a_diag_index, lgs_a_ptr, &
856 nnz, nmax, omega_sor, eps_sor, lgs_x_value, ierr)
861 h_neu(j,i) = lgs_x_value(nr)
864 deallocate(lgs_a_value_trim, lgs_a_index_trim, lgs_a_ptr)
865 deallocate(lgs_a_diag_index, lgs_b_value, lgs_x_value)
871 call lis_matrix_create(lis_comm_world, lgs_a, ierr)
872 call lis_vector_create(lis_comm_world, lgs_b, ierr)
873 call lis_vector_create(lis_comm_world, lgs_x, ierr)
875 call lis_matrix_set_size(lgs_a, 0, nmax, ierr)
876 call lis_vector_set_size(lgs_b, 0, nmax, ierr)
877 call lis_vector_set_size(lgs_x, 0, nmax, ierr)
881 do nc=lgs_a_ptr(nr), lgs_a_ptr(nr+1)-1
882 call lis_matrix_set_value(lis_ins_value, nr, lgs_a_index(nc), &
883 lgs_a_value(nc), lgs_a, ierr)
886 call lis_vector_set_value(lis_ins_value, nr, lgs_b_value(nr), lgs_b, ierr)
887 call lis_vector_set_value(lis_ins_value, nr, lgs_x_value(nr), lgs_x, ierr)
891 call lis_matrix_set_type(lgs_a, lis_matrix_csr, ierr)
892 call lis_matrix_assemble(lgs_a, ierr)
896 call lis_solver_create(solver, ierr)
898 ch_solver_set_option =
'-i bicg -p ilu '// &
899 '-maxiter 1000 -tol 1.0e-04 -initx_zeros false' 901 call lis_solver_set_option(trim(ch_solver_set_option), solver, ierr)
904 call lis_solve(lgs_a, lgs_b, lgs_x, solver, ierr)
907 call lis_solver_get_iter(solver, iter, ierr)
908 write(6,
'(10x,a,i0)')
'calc_thk_impl: iter = ', iter
911 call lis_vector_gather(lgs_x, lgs_x_value, ierr)
912 call lis_matrix_destroy(lgs_a, ierr)
913 call lis_vector_destroy(lgs_b, ierr)
914 call lis_vector_destroy(lgs_x, ierr)
915 call lis_solver_destroy(solver, ierr)
920 h_neu(j,i) = lgs_x_value(nr)
923 deallocate(lgs_a_value, lgs_a_index, lgs_a_ptr)
924 deallocate(lgs_a_diag_index, lgs_b_value, lgs_x_value)
936 call thk_adjust(time, dtime)
949 real(dp),
intent(in) :: dtime
950 real(dp),
intent(in) :: z_mar
952 real(dp),
dimension(0:JMAX,0:IMAX) :: H_smb, H_adv
953 logical,
dimension(0:JMAX,0:IMAX) :: flag_ocean_point
969 h_smb(0,:) = 0.0_dp; h_smb(1,:) = 0.0_dp
970 h_smb(jmax-1,:) = 0.0_dp; h_smb(jmax,:) = 0.0_dp
971 h_smb(:,0) = 0.0_dp; h_smb(:,1) = 0.0_dp
972 h_smb(:,imax-1) = 0.0_dp; h_smb(:,imax) = 0.0_dp
976 h_smb(2:jmax-2,2:imax-2) &
977 = max(
h_neu(2:jmax-2,2:imax-2) + dtime*
mb_source(2:jmax-2,2:imax-2), &
986 where (
maske == 2_i2b .and. h_smb >
eps_dp) h_smb = 0.0_dp
987 #else /* margin = 2 or 3 */ 989 where (
zl0 < z_mar) h_smb = 0.0_dp
993 where (h_smb < 0.0_dp) h_smb = 0.0_dp
1010 where (
maske == 2_i2b)
1011 flag_ocean_point = .true.
1013 flag_ocean_point = .false.
1062 where (flag_ocean_point)
1080 subroutine thk_adjust(time, dtime)
1084 real(dp),
intent(in) :: time, dtime
1086 real(dp),
dimension(0:JMAX,0:IMAX) :: H_neu_tmp
1087 real(dp) :: dtime_inv
1144 stop
' >>> thk_adjust: THK_EVOL must be between 0 and 4!' 1149 dtime_inv = 1.0_dp/dtime
1153 end subroutine thk_adjust
1159 n_calc_thk_mask_update_aux)
1165 real(dp),
intent(in) :: time
1166 real(dp),
intent(in) :: dtime, dxi, deta
1167 real(dp),
intent(in) :: z_sl, z_mar
1168 integer(i2b),
intent(in) :: n_calc_thk_mask_update_aux
1170 integer(i4b) :: i, j
1171 real(dp) :: year_sec_inv
1172 real(dp),
dimension(0:JMAX,0:IMAX) :: H_neu_tmp
1173 real(dp) :: dtime_inv
1177 year_sec_inv = 1.0_dp/year_sec
1185 dtime_inv = 1.0_dp/dtime
1187 if (n_calc_thk_mask_update_aux == 1_i2b)
then 1189 else if (n_calc_thk_mask_update_aux == 2_i2b)
then 1190 call calc_thk_mask_update_aux2(time, dtime, dtime_inv, z_sl, z_mar)
1191 else if (n_calc_thk_mask_update_aux == 3_i2b)
then 1192 call calc_thk_mask_update_aux3(time, dtime, dtime_inv, z_sl, z_mar)
1194 write(6, fmt=
'(a)')
' >>> calc_thk_mask_update:' 1195 write(6, fmt=
'(a)')
' n_calc_thk_mask_update_aux has no valid value!' 1231 if ( ((
maske(j,i) == 0_i2b).or.(
maske(j,i) == 3_i2b)) &
1232 .and.(
n_cts(j,i) == 1_i2b))
then 1237 - 0.5_dp*(
vz_c(0,j,i)+
vz_t(ktmax-1,j,i))
1278 real(dp),
intent(in) :: time, dtime, dtime_inv, z_sl, z_mar
1280 integer(i4b) :: i, j, kc, kt
1281 real(dp) :: rhosw_rho_ratio, rho_rhosw_ratio
1282 logical :: flag_calving_event
1284 real(dp),
dimension(0:JMAX,0:IMAX) :: H_sea_neu, H_balance
1286 real(dp),
parameter :: H_isol_max = 1.0e+03_dp
1297 h_sea_neu = z_sl -
zl_neu 1298 h_balance = h_sea_neu*rhosw_rho_ratio
1305 if (
maske(j,i) <= 1_i2b)
then 1319 if (
h_neu(j,i) < (rhosw_rho_ratio*h_sea_neu(j,i)) )
then 1321 #if (MARINE_ICE_FORMATION==1) 1323 #elif (MARINE_ICE_FORMATION==2) 1330 #if (MARINE_ICE_CALVING==2 \ 1331 || marine_ice_calving==4 \
1332 || marine_ice_calving==6)
1333 if (
zl0(j,i) < z_mar)
maske(j,i) = 2_i2b
1334 #elif (MARINE_ICE_CALVING==3 \ 1335 || marine_ice_calving==5 \
1336 || marine_ice_calving==7)
1373 if (
maske(j,i) == 1_i2b)
then 1378 else if (
maske(j,i) == 2_i2b)
then 1383 else if (
maske(j,i) == 3_i2b)
then 1385 write(6, fmt=
'(a)')
' >>> calc_thk_mask_update_aux1:' 1386 write(6, fmt=
'(a)')
' maske(j,i)==3 not allowed for' 1387 write(6, fmt=
'(a)')
' SIA-only dynamics!' 1400 if (
maske(j,i) == 0_i2b)
then 1404 if ( ((
maske(j,i+1) == 1_i2b).or.(
maske(j,i+1) == 2_i2b)) &
1405 .and. ((
maske(j,i-1) == 1_i2b).or.(
maske(j,i-1) == 2_i2b)) &
1406 .and. ((
maske(j+1,i) == 1_i2b).or.(
maske(j+1,i) == 2_i2b)) &
1407 .and. ((
maske(j-1,i) == 1_i2b).or.(
maske(j-1,i) == 2_i2b)) &
1415 if (
n_cts(j,i) == 1)
then 1444 subroutine calc_thk_mask_update_aux2(time, dtime, dtime_inv, z_sl, z_mar)
1448 real(dp),
intent(in) :: time, dtime, dtime_inv, z_sl, z_mar
1450 integer(i4b) :: i, j, kc, kt
1451 real(dp) :: rhosw_rho_ratio, rho_rhosw_ratio
1452 logical :: flag_calving_event
1454 real(dp),
dimension(0:JMAX,0:IMAX) :: H_sea_neu, H_balance
1456 real(dp),
parameter :: H_isol_max = 1.0e+03_dp
1467 h_sea_neu = z_sl -
zl_neu 1468 h_balance = h_sea_neu*rhosw_rho_ratio
1475 if (
maske(j,i) <= 1_i2b)
then 1489 if (
h_neu(j,i) < (rhosw_rho_ratio*h_sea_neu(j,i)) )
then 1491 #if (MARINE_ICE_FORMATION==1) 1493 #elif (MARINE_ICE_FORMATION==2) 1500 #if (MARINE_ICE_CALVING==2 \ 1501 || marine_ice_calving==4 \
1502 || marine_ice_calving==6)
1503 if (
zl0(j,i) < z_mar)
maske(j,i) = 2_i2b
1504 #elif (MARINE_ICE_CALVING==3 \ 1505 || marine_ice_calving==5 \
1506 || marine_ice_calving==7)
1544 if (
maske(j,i) == 1_i2b)
then 1549 else if (
maske(j,i) == 2_i2b)
then 1554 else if (
maske(j,i) == 3_i2b)
then 1556 write(6, fmt=
'(a)')
' >>> calc_thk_mask_update_aux2:' 1557 write(6, fmt=
'(a)')
' maske(j,i)==3 not allowed for' 1558 write(6, fmt=
'(a)')
' hybrid SIA/SStA dynamics of ice sheets' 1559 write(6, fmt=
'(a)')
' without ice shelves!' 1572 if (
maske(j,i) == 0_i2b)
then 1576 if ( ((
maske(j,i+1) == 1_i2b).or.(
maske(j,i+1) == 2_i2b)) &
1577 .and. ((
maske(j,i-1) == 1_i2b).or.(
maske(j,i-1) == 2_i2b)) &
1578 .and. ((
maske(j+1,i) == 1_i2b).or.(
maske(j+1,i) == 2_i2b)) &
1579 .and. ((
maske(j-1,i) == 1_i2b).or.(
maske(j-1,i) == 2_i2b)) &
1587 if (
n_cts(j,i) == 1)
then 1610 end subroutine calc_thk_mask_update_aux2
1616 subroutine calc_thk_mask_update_aux3(time, dtime, dtime_inv, z_sl, z_mar)
1620 real(dp),
intent(in) :: time, dtime, dtime_inv, z_sl, z_mar
1622 integer(i4b) :: i, j, kc, kt
1623 real(dp) :: rhosw_rho_ratio, rho_rhosw_ratio
1624 logical :: flag_calving_event
1626 real(dp),
dimension(0:JMAX,0:IMAX) :: H_sea_neu, H_balance
1628 real(dp),
parameter :: H_isol_max = 1.0e+03_dp
1639 h_sea_neu = z_sl -
zl_neu 1640 h_balance = h_sea_neu*rhosw_rho_ratio
1649 if ( (
maske(j,i) <= 1_i2b &
1650 .and. (
maske(j,i+1)>1_i2b.or.
maske(j,i-1)>1_i2b &
1651 .or.
maske(j+1,i)>1_i2b.or.
maske(j-1,i)>1_i2b) ) &
1653 (
maske(j,i)>=2_i2b &
1654 .and. (
maske(j,i+1)==0_i2b.or.
maske(j,i-1)==0_i2b &
1655 .or.
maske(j+1,i)==0_i2b.or.
maske(j-1,i)==0_i2b) ) )
then 1659 if (
h_neu(j,i)<h_balance(j,i).and.
zl_neu(j,i)<z_sl)
then 1673 if (
zl_neu(j,i)>z_sl)
then 1687 else if (
maske(j,i) <= 1_i2b)
then 1705 if (
h_neu(j,i)<h_balance(j,i).and.
zl_neu(j,i)<z_sl)
then 1719 if (
zl_neu(j,i)>z_sl)
then 1738 #if (ICE_SHELF_CALVING==2) 1743 flag_calving_event = .false.
1748 if ( (
maske(j,i)==3_i2b) &
1750 ( (
maske(j,i+1)==2_i2b) &
1751 .or.(
maske(j,i-1)==2_i2b) &
1752 .or.(
maske(j+1,i)==2_i2b) &
1753 .or.(
maske(j-1,i)==2_i2b) ) &
1759 flag_calving_event = .true.
1766 if (.not.flag_calving_event)
exit 1770 #elif (ICE_SHELF_CALVING==9) 1772 #if (defined(MISMIPP)) 1777 if ((
maske(j,i)==3_i2b).and.(
xi(i) > lx))
then 1785 write(6, fmt=
'(a)')
' >>> calc_thk_mask_update_aux3:' 1786 write(6, fmt=
'(a)')
' Option ICE_SHELF_CALVING==9' 1787 write(6, fmt=
'(a)')
' only defined for MISMIP+!' 1812 if ( (
maske(j,i)==0_i2b) &
1814 ( (
maske(j,i+1)==3_i2b) &
1815 .or.(
maske(j,i-1)==3_i2b) &
1816 .or.(
maske(j+1,i)==3_i2b) &
1817 .or.(
maske(j-1,i)==3_i2b) ) &
1821 if ( (
maske(j,i)==3_i2b) &
1823 ( (
maske(j,i+1)==0_i2b) &
1824 .or.(
maske(j,i-1)==0_i2b) &
1825 .or.(
maske(j+1,i)==0_i2b) &
1826 .or.(
maske(j-1,i)==0_i2b) ) &
1830 if ( (
maske(j,i)==3_i2b) &
1832 ( (
maske(j,i+1)==2_i2b) &
1833 .or.(
maske(j,i-1)==2_i2b) &
1834 .or.(
maske(j+1,i)==2_i2b) &
1835 .or.(
maske(j-1,i)==2_i2b) ) &
1839 if ( (
maske(j,i)==2_i2b) &
1841 ( (
maske(j,i+1)==3_i2b) &
1842 .or.(
maske(j,i-1)==3_i2b) &
1843 .or.(
maske(j+1,i)==3_i2b) &
1844 .or.(
maske(j-1,i)==3_i2b) ) &
1854 if ( (
maske(j,i)==2_i2b) &
1855 .and. (
maske(j+1,i)==3_i2b) &
1860 if ( (
maske(j,i)==2_i2b) &
1861 .and. (
maske(j-1,i)==3_i2b) &
1870 if ( (
maske(j,i)==2_i2b) &
1871 .and. (
maske(j,i+1)==3_i2b) &
1876 if ( (
maske(j,i)==2_i2b) &
1877 .and. (
maske(j,i-1)==3_i2b) &
1889 if (
maske(j,i) == 1_i2b)
then 1894 else if (
maske(j,i) == 2_i2b)
then 1900 else if (
maske(j,i) == 3_i2b)
then 1917 if ( (
maske(j,i) == 0_i2b).or.(
maske(j,i) == 3_i2b) )
then 1922 if ( ((
maske(j,i+1) == 1_i2b).or.(
maske(j,i+1) == 2_i2b)) &
1923 .and. ((
maske(j,i-1) == 1_i2b).or.(
maske(j,i-1) == 2_i2b)) &
1924 .and. ((
maske(j+1,i) == 1_i2b).or.(
maske(j+1,i) == 2_i2b)) &
1925 .and. ((
maske(j-1,i) == 1_i2b).or.(
maske(j-1,i) == 2_i2b)) &
1933 if (
n_cts(j,i) == 1)
then 1956 end subroutine calc_thk_mask_update_aux3
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) vx_c
vx_c(kc,j,i): Velocity in x-direction in the upper (kc) ice domain (at (i+1/2,j,kc)) ...
Computation of the ice thickness.
subroutine, public calc_thk_init()
Initialisations for the ice thickness computation.
subroutine sor_sprs(lgs_a_value, lgs_a_index, lgs_a_diag_index, lgs_a_ptr, lgs_b_value, nnz, nmax, omega, eps_sor, lgs_x_value, ierr)
SOR solver for a system of linear equations lgs_a*lgs_x=lgs_b [matrix storage: compressed sparse row ...
logical, dimension(0:jmax, 0:imax) flag_inner_point
flag_inner_point(j,i): Inner-point flag. .true.: inner point, .false.: margin point ...
integer(i2b), dimension(0:jmax, 0:imax) n_cts
n_cts(j,i): Mask for thermal conditions. -1: cold ice base, 0: temperate ice base with cold ice above...
real(dp), dimension(0:jmax, 0:imax) insq_g22_sgy
insq_g22_sgy(j,i): Inverse square root of g22, at (i,j+1/2)
real(dp), dimension(0:jmax, 0:imax) vx_m
vx_m(j,i): Mean (depth-averaged) velocity in x-direction, at (i+1/2,j)
real(dp), dimension(0:jmax, 0:imax) mb_source_apl
mb_source_apl(j,i): Applied mass balance source (SMB, BMB, calving)
real(dp), parameter eps
eps: Small number
integer(i2b), dimension(0:jmax, 0:imax) maske_maxextent
maske_maxextent(j,i): Maximum ice extent mask. 0: not allowed to glaciate, 1: allowed to glaciate...
subroutine calc_thk_mask_update_aux1(time, dtime, dtime_inv, z_sl, z_mar)
Update of the ice-land-ocean mask for SIA-only dynamics of ice sheets without ice shelves...
real(dp), dimension(0:jmax, 0:imax) calv_grounded
calv_grounded(j,i): Calving rate of grounded ice
logical, dimension(0:jmax, 0:imax) flag_grounding_line_1
flag_grounding_line_1(j,i): Grounding line flag. .true.: grounding line point (grounded ice point wit...
real(dp), dimension(0:jmax, 0:imax) insq_g11_sgx
insq_g11_sgx(j,i): Inverse square root of g11, at (i+1/2,j)
Several mathematical tools used by SICOPOLIS.
real(dp), dimension(0:jmax, 0:imax) insq_g22_g
insq_g22_g(j,i): Inverse square root of g22 on grid point (i,j)
real(dp), dimension(0:jmax, 0:imax) dh_t_dtau
dH_t_dtau(j,i): Derivative of H_t by tau (time)
real(dp), dimension(0:jmax, 0:imax) dh_c_dtau
dH_c_dtau(j,i): Derivative of H_c by tau (time)
integer(i2b), dimension(0:jmax, 0:imax) maske
maske(j,i): Ice-land-ocean mask. 0: grounded ice, 1: ice-free land, 2: ocean, 3: floating ice ...
real(dp), dimension(0:jmax, 0:imax) sq_g22_sgx
sq_g22_sgx(j,i): Square root of g22, at (i+1/2,j)
subroutine, public calc_thk_impl(time, dtime, dxi, deta, z_mar, mean_accum, ii, jj, nn)
Over-implicit solver for the general ice thickness equation.
real(dp), dimension(0:jmax, 0:imax) as_perp
as_perp(j,i): Accumulation-ablation function at the ice surface (SMB)
real(dp), dimension(0:jmax, 0:imax) area
area(j,i): Area of grid cell associated with grid point (i,j)
logical, save flag_solver_explicit
real(dp), dimension(0:jmax, 0:imax) runoff
runoff(j,i): Runoff rate at the ice surface
subroutine, public calc_thk_sia_expl(time, dtime, dxi, deta, z_mar)
Explicit solver for the diffusive SIA ice surface equation.
subroutine topograd_1(dxi, deta, n_switch)
Calculation of topography gradients on the staggered grid and on the grid points (the latter by secon...
real(dp), dimension(0:jmax, 0:imax) h_t_neu
(.)_neu: New value of quantity (.) computed during an integration step
real(dp) target_topo_tau
target_topo_tau: Relaxation time for target-topography adjustment
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) vz_t
vz_t(kt,j,i): Velocity in z-direction in the lower (kt) ice domain (at (i,j,kt+1/2)) ...
subroutine, public calc_thk_expl(time, dtime, dxi, deta, z_mar)
Explicit solver for the general ice thickness equation.
real(dp), dimension(0:jmax, 0:imax) vy_m
vy_m(j,i): Mean (depth-averaged) velocity in y-direction, at (i,j+1/2)
real(dp), dimension(0:jmax, 0:imax), save h
real(dp), dimension(0:jmax, 0:imax) runoff_bot
runoff_bot(j,i): Runoff rate (accounting at ice bottom)
real(dp), dimension(0:jmax, 0:imax) zl_neu
(.)_neu: New value of quantity (.) computed during an integration step
Declarations of global variables for SICOPOLIS (for the ANT domain).
real(dp), dimension(0:jmax, 0:imax) h_diff
h_diff(j,i): Diffusivity of the SIA evolution equation of the ice surface
real(dp), dimension(0:jmax, 0:imax) runoff_top
runoff_top(j,i): Runoff rate (accounting at ice surface)
subroutine apply_mb_source(dtime, z_mar)
Ice thickness evolution due to the source term (surface mass balance, basal mass balance and calving)...
subroutine, public calc_thk_mask_update(time, dtime, dxi, deta, z_sl, z_mar, n_calc_thk_mask_update_aux)
Update of the ice-land-ocean mask etc.
Declarations of kind types for SICOPOLIS.
real(dp), dimension(0:jmax, 0:imax) zs
zs(j,i): Coordinate z of the surface topography
real(dp), dimension(0:jmax, 0:imax) h_c
H_c(j,i): Thickness of ice in the upper (kc) domain (thickness of the cold-ice layer for POLY...
subroutine topograd_2(dxi, deta, n_switch)
Calculation of topography gradients on the staggered grid and on the grid points (the latter by fourt...
real(dp), dimension(0:jmax, 0:imax), save mb_source
real(dp), dimension(0:jmax, 0:imax) dzm_dtau
dzm_dtau(j,i): Derivative of zm by tau (time)
real(dp), dimension(0:jmax, 0:imax) sq_g11_sgy
sq_g11_sgy(j,i): Square root of g11, at (i,j+1/2)
integer(i2b), dimension(0:jmax, 0:imax) maske_old
maske_old(j,i): Old value of maske (at the previous time step)
real(dp), parameter eps_dp
eps_dp: Small number to double-precision accuracy
real(dp), dimension(0:jmax, 0:imax) zs_neu
(.)_neu: New value of quantity (.) computed during an integration step
real(dp), dimension(0:jmax, 0:imax) zl0
zl0(j,i): zl for isostatically relaxed ice-free conditions
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) vz_c
vz_c(kc,j,i): Velocity in z-direction in the upper (kc) ice domain (at (i,j,kc+1/2)) ...
real(dp) time_target_topo_final
time_target_topo_final: Final time for target-topography adjustment
real(dp), dimension(0:jmax, 0:imax) h_c_neu
(.)_neu: New value of quantity (.) computed during an integration step
real(dp), dimension(0:jmax, 0:imax), save h_neu
real(dp), dimension(0:jmax, 0:imax) q_b_tot
Q_b_tot(j,i): Sum of Q_bm and Q_tld.
real(dp), dimension(0:jmax, 0:imax) disc_top
disc_top(j,i): Ice discharge rate (accounting at ice surface)
real(dp) time_target_topo_init
time_target_topo_init: Initial time for target-topography adjustment
real(dp), dimension(0:jmax, 0:imax) dzm_dxi_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
real(dp), dimension(0:jmax, 0:imax) am_perp
am_perp(j,i): Ice volume flux across the z=zm interface
Calculation of topography gradients on the staggered grid and on the grid points (including length re...
real(dp), dimension(0:jmax, 0:imax) dzm_deta_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
real(dp), dimension(0:jmax, 0:imax) as_perp_apl
as_perp_apl(j,i): Applied accumulation-ablation function (SMB)
logical, dimension(0:jmax, 0:imax) flag_calving_front_1
flag_calving_front_1(j,i): Calving front flag. .true.: calving front point (floating ice point with a...
real(dp), dimension(0:jmax, 0:imax) insq_g11_g
insq_g11_g(j,i): Inverse square root of g11 on grid point (i,j)
real(dp), dimension(0:imax) xi
xi(i): Coordinate xi (= x) of grid point i
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) vy_c
vy_c(kc,j,i): Velocity in y-direction in the upper (kc) ice domain (at (i,j+1/2,kc)) ...
real(dp) rho_sw
RHO_SW: Density of sea water.
logical, dimension(0:jmax, 0:imax) flag_calving_front_2
flag_calving_front_2(j,i): Calving front flag. .true.: calving front point (ocean point with at least...
real(dp), dimension(0:jmax, 0:imax), save mb_adjust
real(dp), dimension(0:jmax, 0:imax) zb_neu
(.)_neu: New value of quantity (.) computed during an integration step
real(dp), dimension(0:jmax, 0:imax) zm_neu
(.)_neu: New value of quantity (.) computed during an integration step
integer(i2b), dimension(0:jmax, 0:imax) mask_run
mask_run(j,i): mask indicating melt type. 2: visible (ocean, for later developments), 1: visible (grounded ice), -1: hidden on land, -2: hidden in ocean
real(dp) rho
RHO: Density of ice.
real(dp), dimension(0:jmax, 0:imax) h_target
H_target(j,i): Target topography (ice thickness)
subroutine, public calc_thk_sia_impl(time, dtime, dxi, deta, z_mar, mean_accum, ii, jj, nn)
Over-implicit solver for the diffusive SIA ice surface equation.
real(dp), dimension(0:jmax, 0:imax) zb
zb(j,i): Coordinate z of the ice base
logical, dimension(0:jmax, 0:imax) flag_grounding_line_2
flag_grounding_line_2(j,i): Grounding line flag. .true.: grounding line point (floating ice point wit...
real(dp), dimension(0:jmax, 0:imax) am_perp_st
am_perp_st(j,i): Steady-state part of am_perp (without contribution of dzm_dtau)
real(dp), dimension(0:jmax, 0:imax) dzb_dtau
dzb_dtau(j,i): Derivative of zb by tau (time)
integer(i2b), dimension(0:jmax, 0:imax) maske_target
maske_target(j,i): Target topography (ice-land-ocean mask)
real(dp), dimension(0:jmax, 0:imax) h_t
H_t(j,i): Thickness of ice in the lower (kt) domain (thickness of the temperate layer for POLY...
Declarations of global variables for SICOPOLIS.
real(dp), dimension(0:jmax, 0:imax) dzs_dtau
dzs_dtau(j,i): Derivative of zs by tau (time)
real(dp), dimension(0:jmax, 0:imax) dzl_dtau
dzl_dtau(j,i): Derivative of zl by tau (time)