Skip to content

Commit

Permalink
Removed blocks of commented code. Added parentheses in calc_shelf_visc
Browse files Browse the repository at this point in the history
  • Loading branch information
@OlgaSergienko
OlgaSergienko committed Mar 9, 2021
1 parent 43dadc1 commit b47e493
Show file tree
Hide file tree
Showing 2 changed files with 18 additions and 222 deletions.
56 changes: 17 additions & 39 deletions src/ice_shelf/MOM_ice_shelf_dynamics.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -579,12 +579,11 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_
if (CS%id_taudx_shelf > 0) call post_data(CS%id_taudx_shelf, CS%taudx_shelf,CS%diag)
if (CS%id_taudy_shelf > 0) call post_data(CS%id_taudy_shelf, CS%taudy_shelf,CS%diag)
if (CS%id_visc_shelf > 0) call post_data(CS%id_visc_shelf, CS%ice_visc,CS%diag)
endif
!!! OVS vertically integrated temperature
! CS%id_t_shelf = register_diag_field('ocean_model','t_shelf',CS%diag%axesT1, Time, &
! 'T of ice', 'oC')
! CS%id_t_mask = register_diag_field('ocean_model','tmask',CS%diag%axesT1, Time, &
! 'mask for T-nodes', 'none')
endif
endif

end subroutine initialize_ice_shelf_dyn
Expand Down Expand Up @@ -875,9 +874,7 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf,taudx,taudy, ite
enddo

call calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, CS%OD_av)
call pass_vector(taudx, taudy, G%domain, TO_ALL, BGRID_NE) !OVS 02/01/21
! call pass_var(taudx, G%Domain) !OVS 01/21/21
! call pass_var(taudy, G%Domain) !OVS 01/21/21
call pass_vector(taudx, taudy, G%domain, TO_ALL, BGRID_NE)
! This is to determine which cells contain the grounding line, the criterion being that the cell
! is ice-covered, with some nodes floating and some grounded flotation condition is estimated by
! assuming topography is cellwise constant and H is bilinear in a cell; floating where
Expand Down Expand Up @@ -917,12 +914,10 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf,taudx,taudy, ite
call bilinear_shape_fn_grid(G, i, j, Phi(:,:,i,j))
enddo ; enddo

call calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) !OVS 02/24/21
call calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf)
call pass_var(CS%ice_visc, G%domain)
! call pass_vector(CS%ice_visc, G%domain, TO_ALL, BGRID_NE) !OVS 02/11/21
call calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf)
call pass_var(CS%basal_traction, G%domain)
! call pass_vector(CS%ice_visc,CS%basal_traction, G%domain, TO_ALL, BGRID_NE) !OVS 02/11/21

! This makes sure basal stress is only applied when it is supposed to be
do j=G%jsd,G%jed ; do i=G%isd,G%ied
Expand All @@ -937,7 +932,7 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf,taudx,taudy, ite
call CG_action(Au, Av, u_shlf, v_shlf, Phi, Phisub, CS%umask, CS%vmask, ISS%hmask, H_node, &
CS%ice_visc, float_cond, G%bathyT, CS%basal_traction, &
G, US, G%isc-1, G%iec+1, G%jsc-1, G%jec+1, rhoi_rhow)
call pass_vector(Au,Av,G%domain) !OVS pass Au and Av
call pass_vector(Au,Av,G%domain)
if (CS%nonlin_solve_err_mode == 1) then
err_init = 0 ; err_tempu = 0 ; err_tempv = 0
do J=G%IscB,G%JecB ; do I=G%IscB,G%IecB
Expand Down Expand Up @@ -971,7 +966,7 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf,taudx,taudy, ite
write(mesg,*) "ice_shelf_solve_outer: linear solve done in ",iters," iterations"
call MOM_mesg(mesg, 5)

call calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) !OVS 02/24/21
call calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf)
call pass_var(CS%ice_visc, G%domain)
call calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf)
call pass_var(CS%basal_traction, G%domain)
Expand Down Expand Up @@ -1823,10 +1818,8 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD)
S(:,:) = BASE(:,:) + ISS%h_shelf(:,:)

! check whether the ice is floating or grounded
! do j=jsc-1,jec+1 !OVS 02/02/21
! do i=isc-1,iec+1 !OVS 02/02/21
do j=jsc-G%domain%njhalo,jec+G%domain%njhalo !OVS 02/02/21
do i=isc-G%domain%nihalo,iec+G%domain%nihalo !OVS 02/02/21
do j=jsc-G%domain%njhalo,jec+G%domain%njhalo
do i=isc-G%domain%nihalo,iec+G%domain%nihalo

! if (ISS%h_shelf(i,j) < rhow/rho * G%bathyT(i,j)) then
if (rhoi_rhow * ISS%h_shelf(i,j) - G%bathyT(i,j) <= 0) then
Expand Down Expand Up @@ -2160,7 +2153,7 @@ subroutine CG_action(uret, vret, u_shlf, v_shlf, Phi, Phisub, umask, vmask, hmas
v_shlf(I-1,J) * Phi(6,2*(jq-1)+iq,i,j) + &
v_shlf(I,J) * Phi(8,2*(jq-1)+iq,i,j)

do iphi=1,2 ; do jphi=1,2 ; Itgt = I-2+iphi ; ;Jtgt = J-2+jphi !Jtgt = J-2-jphi !OVS fix index
do iphi=1,2 ; do jphi=1,2 ; Itgt = I-2+iphi ; ;Jtgt = J-2+jphi
if (umask(Itgt,Jtgt) == 1) uret(Itgt,Jtgt) = uret(Itgt,Jtgt) + 0.25 * ice_visc(i,j) * &
((4*ux+2*vy) * Phi(2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq,i,j) + &
(uy+vx) * Phi(2*(2*(jphi-1)+iphi),2*(jq-1)+iq,i,j))
Expand Down Expand Up @@ -2338,7 +2331,7 @@ subroutine matrix_diagonal(CS, G, US, float_cond, H_node, ice_visc, basal_trac,
if (float_cond(i,j) == 1) then
Hcell(:,:) = H_node(i-1:i,j-1:j)
call CG_diagonal_subgrid_basal(Phisub, Hcell, G%bathyT(i,j), dens_ratio, sub_ground)
do iphi=1,2 ; do jphi=1,2 ; Itgt = I-2+iphi ; Jtgt = J-2+jphi !Jtgt = J-2-jphi !OVS fix index
do iphi=1,2 ; do jphi=1,2 ; Itgt = I-2+iphi ; Jtgt = J-2+jphi
if (CS%umask(Itgt,Jtgt) == 1) then
u_diagonal(Itgt,Jtgt) = u_diagonal(Itgt,Jtgt) + sub_ground(iphi,jphi) * basal_trac(i,j)
v_diagonal(Itgt,Jtgt) = v_diagonal(Itgt,Jtgt) + sub_ground(iphi,jphi) * basal_trac(i,j)
Expand Down Expand Up @@ -2479,7 +2472,7 @@ subroutine apply_boundary_values(CS, ISS, G, US, time, Phisub, H_node, ice_visc,
CS%v_bdry_val(I-1,J) * Phi(6,2*(jq-1)+iq) + &
CS%v_bdry_val(I,J) * Phi(8,2*(jq-1)+iq)

do iphi=1,2 ; do jphi=1,2 ; Itgt = I-2+iphi ; Jtgt = J-2+jphi !Jtgt = J-2-jphi !OVS fix index
do iphi=1,2 ; do jphi=1,2 ; Itgt = I-2+iphi ; Jtgt = J-2+jphi
ilq = 1 ; if (iq == iphi) ilq = 2
jlq = 1 ; if (jq == jphi) jlq = 2

Expand Down Expand Up @@ -2590,14 +2583,14 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf)
! enddo
! call pass_vector(ux, uy, G%domain, TO_ALL, BGRID_NE)
! call pass_vector(vx, vy, G%domain, TO_ALL, BGRID_NE)
ux = ((u_shlf(I,J) + u_shlf(I,J-1) + u_shlf(I,J+1)) - &
(u_shlf(I-1,J) + u_shlf(I-1,J-1) + u_shlf(I-1,J+1))) / (3*G%dxT(i,j))
ux = ((u_shlf(I,J) + (u_shlf(I,J-1) + u_shlf(I,J+1))) - &
(u_shlf(I-1,J) + (u_shlf(I-1,J-1) + u_shlf(I-1,J+1)))) / (3*G%dxT(i,j))
vx = ((v_shlf(I,J) + v_shlf(I,J-1) + v_shlf(I,J+1)) - &
(v_shlf(I-1,J) + v_shlf(I-1,J-1) + v_shlf(I-1,J+1))) / (3*G%dxT(i,j))
uy = ((u_shlf(I,J) + u_shlf(I-1,J) + u_shlf(I+1,J)) - &
(u_shlf(I,J-1) + u_shlf(I-1,J-1) + u_shlf(I+1,J-1))) / (3*G%dyT(i,j))
vy = ((v_shlf(I,J) + v_shlf(I-1,J)+ v_shlf(I+1,J)) - &
(v_shlf(I,J-1) + v_shlf(I-1,J-1)+ v_shlf(I+1,J-1))) / (3*G%dyT(i,j))
(v_shlf(I-1,J) + (v_shlf(I-1,J-1) + v_shlf(I-1,J+1)))) / (3*G%dxT(i,j))
uy = ((u_shlf(I,J) + (u_shlf(I-1,J) + u_shlf(I+1,J))) - &
(u_shlf(I,J-1) + (u_shlf(I-1,J-1) + u_shlf(I+1,J-1)))) / (3*G%dyT(i,j))
vy = ((v_shlf(I,J) + (v_shlf(I-1,J)+ v_shlf(I+1,J))) - &
(v_shlf(I,J-1) + (v_shlf(I-1,J-1)+ v_shlf(I+1,J-1)))) / (3*G%dyT(i,j))
! ux = ((u_shlf(I,J) + u_shlf(I,J-1)) - (u_shlf(I-1,J) + u_shlf(I-1,J-1))) / (2*G%dxT(i,j))
! vx = ((v_shlf(I,J) + v_shlf(I,J-1)) - (v_shlf(I-1,J) + v_shlf(I-1,J-1))) / (2*G%dxT(i,j))
! uy = ((u_shlf(I,J) + u_shlf(I-1,J)) - (u_shlf(I,J-1) + u_shlf(I-1,J-1))) / (2*G%dyT(i,j))
Expand Down Expand Up @@ -3020,21 +3013,6 @@ subroutine update_velocity_masks(CS, G, hmask, umask, vmask, u_face_mask, v_face
end select
enddo

!if (CS%u_face_mask_bdry(I-1,j) >= 0) then ! Western boundary
! u_face_mask(I-1,j) = CS%u_face_mask_bdry(I-1,j)
! umask(I-1,J-1:J) = 3.
! vmask(I-1,J-1:J) = 0.
!endif

!if (j_off+j == gjsc+1) then ! SoutherN boundary
! v_face_mask(i,J-1) = 0.
! umask(I-1:I,J-1) = 0.
! vmask(I-1:I,J-1) = 0.
!elseif (j_off+j == gjec) then ! Northern boundary
! v_face_mask(i,J) = 0.
! umask(I-1:I,J) = 0.
! vmask(I-1:I,J) = 0.
!endif

if (i < G%ied) then
if ((hmask(i+1,j) == 0) .OR. (hmask(i+1,j) == 2)) then
Expand Down
184 changes: 1 addition & 183 deletions src/ice_shelf/MOM_ice_shelf_initialize.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -397,190 +397,8 @@ subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_b

enddo
enddo
! if (.not. G%symmetric) then
!! do j=G%jsd,G%jed
!! do i=G%isd,G%ied
! do j=jsc-1,jec+1
! do i=isc-1,iec+1
!! if (((i+G%idg_offset) == (G%domain%nihalo+1)).and.(u_face_mask_bdry(I-1,j) == 3)) then
! if (u_face_mask_bdry(I-1,j) == 3) then
! u_shelf(I-1,J-1) = u_bdry_val(I-1,J-1)
! u_shelf(I-1,J) = u_bdry_val(I-1,J)
! v_shelf(I-1,J-1) = v_bdry_val(I-1,J-1)
! v_shelf(I-1,J) = v_bdry_val(I-1,J)
! endif
!! if (((j+G%jdg_offset) == (G%domain%njhalo+1)).and.(v_face_mask_bdry(i,J-1) == 3)) then
! if (v_face_mask_bdry(I,j-1) == 3) then
! u_shelf(I-1,J-1) = u_bdry_val(I-1,J-1)
! u_shelf(I,J-1) = u_bdry_val(I,J-1)
! v_shelf(I-1,J-1) = v_bdry_val(I-1,J-1)
! v_shelf(I,J-1) = v_bdry_val(I,J-1)
! endif
! enddo
! enddo
! endif
end subroutine initialize_ice_shelf_boundary_channel

!BEGIN MJH
! subroutine initialize_ice_shelf_boundary(u_face_mask_bdry, v_face_mask_bdry, &
! u_flux_bdry_val, v_flux_bdry_val, u_bdry_val, v_bdry_val, h_bdry_val, &
! hmask, G, US, PF )

! type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
! real, dimension(SZIB_(G),SZJ_(G)), &
! intent(inout) :: u_face_mask_bdry !< A boundary-type mask at C-grid u faces
! real, dimension(SZIB_(G),SZJ_(G)), &
! intent(inout) :: u_flux_bdry_val !< The boundary thickness flux through
! !! C-grid u faces [L Z T-1 ~> m2 s-1].
! real, dimension(SZI_(G),SZJB_(G)), &
! intent(inout) :: v_face_mask_bdry !< A boundary-type mask at C-grid v faces
! real, dimension(SZI_(G),SZJB_(G)), &
! intent(inout) :: v_flux_bdry_val !< The boundary thickness flux through
! !! C-grid v faces [L Z T-1 ~> m2 s-1].
! real, dimension(SZIB_(G),SZJB_(G)), &
! intent(inout) :: u_bdry_val !< The zonal ice shelf velocity at open
! !! boundary vertices [L T-1 ~> m s-1].
! real, dimension(SZIB_(G),SZJB_(G)), &
! intent(inout) :: v_bdry_val !< The meridional ice shelf velocity at open
! !! boundary vertices [L T-1 ~> m s-1].
! real, dimension(SZDI_(G),SZDJ_(G)), &
! intent(inout) :: h_bdry_val !< The ice shelf thickness at open boundaries [Z ~> m]
! real, dimension(SZDI_(G),SZDJ_(G)), &
! intent(inout) :: hmask !< A mask indicating which tracer points are
! !! partly or fully covered by an ice-shelf
! type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors
! type(param_file_type), intent(in) :: PF !< A structure to parse for run-time parameters

! character(len=40) :: mdl = "initialize_ice_shelf_boundary" ! This subroutine's name.
! character(len=200) :: config
! logical flux_bdry

! call get_param(PF, mdl, "ICE_BOUNDARY_CONFIG", config, &
! "This specifies how the ice domain boundary is specified. "//&
! "valid values include CHANNEL, FILE and USER.", &
! fail_if_missing=.true.)
! call get_param(PF, mdl, "ICE_BOUNDARY_FLUX_CONDITION", flux_bdry, &
! "This specifies whether mass input is a dirichlet or "//&
! "flux condition", default=.true.)

! select case ( trim(config) )
! case ("CHANNEL")
! call initialize_ice_shelf_boundary_channel(u_face_mask_bdry, &
! v_face_mask_bdry, u_flux_bdry_val, v_flux_bdry_val, &
! u_bdry_val, v_bdry_val, h_bdry_val, hmask, G, &
! flux_bdry, PF)
! case ("FILE"); call MOM_error(FATAL,"MOM_initialize: "// &
! "Unrecognized topography setup "//trim(config))
! case ("USER"); call MOM_error(FATAL,"MOM_initialize: "// &
! "Unrecognized topography setup "//trim(config))
! case default ; call MOM_error(FATAL,"MOM_initialize: "// &
! "Unrecognized topography setup "//trim(config))
! end select

! end subroutine initialize_ice_shelf_boundary

! subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_bdry, &
! u_flux_bdry_val, v_flux_bdry_val, u_bdry_val, v_bdry_val, h_bdry_val, &
! hmask, G, flux_bdry, US, PF )

! type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
! real, dimension(SZIB_(G),SZJ_(G)), &
! intent(inout) :: u_face_mask_bdry !< A boundary-type mask at C-grid u faces
! real, dimension(SZIB_(G),SZJ_(G)), &
! intent(inout) :: u_flux_bdry_val !< The boundary thickness flux through
! !! C-grid u faces [L Z T-1 ~> m2 s-1].
! real, dimension(SZI_(G),SZJB_(G)), &
! intent(inout) :: v_face_mask_bdry !< A boundary-type mask at C-grid v faces
! real, dimension(SZI_(G),SZJB_(G)), &
! intent(inout) :: v_flux_bdry_val !< The boundary thickness flux through
! !! C-grid v faces [L Z T-1 ~> m2 s-1].
! real, dimension(SZIB_(G),SZJB_(G)), &
! intent(inout) :: u_bdry_val !< The zonal ice shelf velocity at open
!! boundary vertices [L T-1 ~> m s-1].
! real, dimension(SZIB_(G),SZJB_(G)), &
! intent(inout) :: v_bdry_val !< The meridional ice shelf velocity at open
!! boundary vertices [L T-1 ~> m s-1].
! real, dimension(SZDI_(G),SZDJ_(G)), &
! intent(inout) :: h_bdry_val !< The ice shelf thickness at open boundaries [Z ~> m]
! real, dimension(SZDI_(G),SZDJ_(G)), &
! intent(inout) :: hmask !< A mask indicating which tracer points are
! !! partly or fully covered by an ice-shelf
! logical, intent(in) :: flux_bdry !< If true, use mass fluxes as the boundary value.
! type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors
! type(param_file_type), intent(in) :: PF !< A structure to parse for run-time parameters

! character(len=40) :: mdl = "initialize_ice_shelf_boundary_channel" ! This subroutine's name.
! integer :: i, j, isd, jsd, is, js, iegq, jegq, giec, gjec, gisc, gjsc, isc, jsc, iec, jec, ied, jed
! real :: input_thick ! The input ice shelf thickness [Z ~> m]
! real :: input_flux ! The input ice flux per unit length [L Z T-1 ~> m2 s-1]
! real :: lenlat, len_stress

! call get_param(PF, mdl, "LENLAT", lenlat, fail_if_missing=.true.)

! call get_param(PF, mdl, "INPUT_FLUX_ICE_SHELF", input_flux, &
! "volume flux at upstream boundary", &
! units="m2 s-1", default=0., scale=US%m_s_to_L_T*US%m_to_Z)
! call get_param(PF, mdl, "INPUT_THICK_ICE_SHELF", input_thick, &
! "flux thickness at upstream boundary", &
! units="m", default=1000., scale=US%m_to_Z)
! call get_param(PF, mdl, "LEN_SIDE_STRESS", len_stress, &
! "maximum position of no-flow condition in along-flow direction", &
! units="km", default=0.)

! call MOM_mesg(mdl//": setting boundary")

! isd = G%isd ; ied = G%ied
! jsd = G%jsd ; jed = G%jed
! isc = G%isc ; jsc = G%jsc ; iec = G%iec ; jec = G%jec
! gisc = G%Domain%nihalo ; gjsc = G%Domain%njhalo
! giec = G%Domain%niglobal+gisc ; gjec = G%Domain%njglobal+gjsc

! do j=jsd,jed
! do i=isd,ied

! ! upstream boundary - set either dirichlet or flux condition

! if ((i+G%idg_offset) == G%domain%nihalo+1) then
! if (flux_bdry) then
! u_face_mask_bdry(i-1,j) = 4.0
! u_flux_bdry_val(i-1,j) = input_flux
! else
! hmask(i-1,j) = 3.0
! h_bdry_val(i-1,j) = input_thick
! u_face_mask_bdry(i-1,j) = 3.0
! u_bdry_val(i-1,j-1) = (1 - ((G%geoLatBu(i-1,j-1) - 0.5*lenlat)*2./lenlat)**2) * &
! 1.5 * input_flux / input_thick
! u_bdry_val(i-1,j) = (1 - ((G%geoLatBu(i-1,j) - 0.5*lenlat)*2./lenlat)**2) * &
! 1.5 * input_flux / input_thick
! endif
! endif

! ! side boundaries: no flow

! if (G%jdg_offset+j == gjsc+1) then !bot boundary
! if (len_stress == 0. .OR. G%geoLonCv(i,j-1) <= len_stress) then
! v_face_mask_bdry(i,j-1) = 0.
! else
! v_face_mask_bdry(i,j-1) = 1.
! endif
! elseif (G%jdg_offset+j == gjec) then !top boundary
! if (len_stress == 0. .OR. G%geoLonCv(i,j-1) <= len_stress) then
! v_face_mask_bdry(i,j) = 0.
! else
! v_face_mask_bdry(i,j) = 1.
! endif
! endif

! ! downstream boundary - CFBC

! if (i+G%idg_offset == giec) then
! u_face_mask_bdry(i,j) = 2.0
! endif

! enddo
! enddo

!END MJH end subroutine initialize_ice_shelf_boundary_channel

!> Initialize ice shelf flow from file
subroutine initialize_ice_flow_from_file(u_shelf, v_shelf,ice_visc,float_cond,&
Expand Down Expand Up @@ -642,7 +460,7 @@ subroutine initialize_ice_flow_from_file(u_shelf, v_shelf,ice_visc,float_cond,&

! call MOM_read_data(filename, trim(ushelf_varname), u_shelf, G%Domain, scale=1.0) !/(365.0*86400.0))
! call MOM_read_data(filename,trim(vshelf_varname), v_shelf, G%Domain, scale=1.0) !/(365.0*86400.0))
call MOM_read_data(filename,trim(ice_visc_varname), ice_visc, G%Domain, scale=1.0) !*(365.0*86400.0))
call MOM_read_data(filename,trim(ice_visc_varname), ice_visc, G%Domain, scale=1.0)
call MOM_read_data(filename,trim(floatfr_varname), float_cond, G%Domain, scale=1.)
! call get_param(PF, mdl, "ICE_BOUNDARY_CONFIG", config, &
! "This specifies how the ice domain boundary is specified", &
Expand Down

0 comments on commit b47e493

Please sign in to comment.