Skip to content

Commit

Permalink
Update CICE to consortium master (NOAA-EMC#23)
Browse files Browse the repository at this point in the history 
updates include:

* deprecate upwind advection (CICE-Consortium#508)
* add implicit VP solver (CICE-Consortium#491)
  • Loading branch information
@DeniseWorthen
DeniseWorthen authored Nov 10, 2020
1 parent ac617cd commit f773ef3
Show file tree
Hide file tree
Showing 45 changed files with 4,622 additions and 269 deletions.
44 changes: 9 additions & 35 deletions cicecore/cicedynB/dynamics/ice_dyn_eap.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -122,7 +122,8 @@ subroutine eap (dt)
use ice_dyn_shared, only: fcor_blk, ndte, dtei, &
denom1, uvel_init, vvel_init, arlx1i, &
dyn_prep1, dyn_prep2, stepu, dyn_finish, &
basal_stress_coeff, basalstress
basal_stress_coeff, basalstress, &
stack_velocity_field, unstack_velocity_field
use ice_flux, only: rdg_conv, strairxT, strairyT, &
strairx, strairy, uocn, vocn, ss_tltx, ss_tlty, iceumask, fm, &
strtltx, strtlty, strocnx, strocny, strintx, strinty, taubx, tauby, &
Expand Down Expand Up @@ -354,11 +355,6 @@ subroutine eap (dt)
vicen = vicen (i,j,:,iblk), &
strength = strength(i,j, iblk) )
enddo ! ij

! load velocity into array for boundary updates
fld2(:,:,1,iblk) = uvel(:,:,iblk)
fld2(:,:,2,iblk) = vvel(:,:,iblk)

enddo ! iblk
!$TCXOMP END PARALLEL DO

Expand All @@ -370,18 +366,12 @@ subroutine eap (dt)
call ice_HaloUpdate (strength, halo_info, &
field_loc_center, field_type_scalar)
! velocities may have changed in dyn_prep2
call stack_velocity_field(uvel, vvel, fld2)
call ice_HaloUpdate (fld2, halo_info, &
field_loc_NEcorner, field_type_vector)
call unstack_velocity_field(fld2, uvel, vvel)
call ice_timer_stop(timer_bound)

! unload
!$OMP PARALLEL DO PRIVATE(iblk)
do iblk = 1, nblocks
uvel(:,:,iblk) = fld2(:,:,1,iblk)
vvel(:,:,iblk) = fld2(:,:,2,iblk)
enddo
!$OMP END PARALLEL DO

if (maskhalo_dyn) then
call ice_timer_start(timer_bound)
halomask = 0
Expand Down Expand Up @@ -472,10 +462,6 @@ subroutine eap (dt)
uvel (:,:,iblk), vvel (:,:,iblk), &
Tbu (:,:,iblk))

! load velocity into array for boundary updates
fld2(:,:,1,iblk) = uvel(:,:,iblk)
fld2(:,:,2,iblk) = vvel(:,:,iblk)

!-----------------------------------------------------------------
! evolution of structure tensor A
!-----------------------------------------------------------------
Expand All @@ -501,6 +487,7 @@ subroutine eap (dt)
enddo
!$TCXOMP END PARALLEL DO

call stack_velocity_field(uvel, vvel, fld2)
call ice_timer_start(timer_bound)
if (maskhalo_dyn) then
call ice_HaloUpdate (fld2, halo_info_mask, &
Expand All @@ -510,14 +497,7 @@ subroutine eap (dt)
field_loc_NEcorner, field_type_vector)
endif
call ice_timer_stop(timer_bound)

! unload
!$OMP PARALLEL DO PRIVATE(iblk)
do iblk = 1, nblocks
uvel(:,:,iblk) = fld2(:,:,1,iblk)
vvel(:,:,iblk) = fld2(:,:,2,iblk)
enddo
!$OMP END PARALLEL DO
call unstack_velocity_field(fld2, uvel, vvel)

enddo ! subcycling

Expand Down Expand Up @@ -556,16 +536,12 @@ end subroutine eap
!=======================================================================

! Initialize parameters and variables needed for the eap dynamics
! (based on init_evp)
! (based on init_dyn)

subroutine init_eap (dt)
subroutine init_eap

use ice_blocks, only: nx_block, ny_block
use ice_domain, only: nblocks
use ice_dyn_shared, only: init_evp

real (kind=dbl_kind), intent(in) :: &
dt ! time step

! local variables

Expand Down Expand Up @@ -595,8 +571,6 @@ subroutine init_eap (dt)
file=__FILE__, line=__LINE__)
phi = pi/c12 ! diamond shaped floe smaller angle (default phi = 30 deg)

call init_evp (dt)

!$OMP PARALLEL DO PRIVATE(iblk,i,j)
do iblk = 1, nblocks
do j = 1, ny_block
Expand Down Expand Up @@ -1321,7 +1295,7 @@ subroutine stress_eap (nx_block, ny_block, &
tensionse = -cym(i,j)*uvel(i ,j-1) - dyt(i,j)*uvel(i-1,j-1) &
+ cxp(i,j)*vvel(i ,j-1) - dxt(i,j)*vvel(i ,j )

! shearing strain rate = e_12
! shearing strain rate = 2*e_12
shearne = -cym(i,j)*vvel(i ,j ) - dyt(i,j)*vvel(i-1,j ) &
- cxm(i,j)*uvel(i ,j ) - dxt(i,j)*uvel(i ,j-1)
shearnw = -cyp(i,j)*vvel(i-1,j ) + dyt(i,j)*vvel(i ,j ) &
Expand Down
115 changes: 38 additions & 77 deletions cicecore/cicedynB/dynamics/ice_dyn_evp.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -94,7 +94,7 @@ subroutine evp (dt)
ice_timer_start, ice_timer_stop, timer_evp_1d, timer_evp_2d
use ice_dyn_evp_1d, only: ice_dyn_evp_1d_copyin, ice_dyn_evp_1d_kernel, &
ice_dyn_evp_1d_copyout
use ice_dyn_shared, only: kevp_kernel
use ice_dyn_shared, only: kevp_kernel, stack_velocity_field, unstack_velocity_field

real (kind=dbl_kind), intent(in) :: &
dt ! time step
Expand Down Expand Up @@ -297,10 +297,6 @@ subroutine evp (dt)
strength = strength(i,j, iblk) )
enddo ! ij

! load velocity into array for boundary updates
fld2(:,:,1,iblk) = uvel(:,:,iblk)
fld2(:,:,2,iblk) = vvel(:,:,iblk)

enddo ! iblk
!$TCXOMP END PARALLEL DO

Expand All @@ -312,18 +308,12 @@ subroutine evp (dt)
call ice_HaloUpdate (strength, halo_info, &
field_loc_center, field_type_scalar)
! velocities may have changed in dyn_prep2
call stack_velocity_field(uvel, vvel, fld2)
call ice_HaloUpdate (fld2, halo_info, &
field_loc_NEcorner, field_type_vector)
call unstack_velocity_field(fld2, uvel, vvel)
call ice_timer_stop(timer_bound)

! unload
!$OMP PARALLEL DO PRIVATE(iblk)
do iblk = 1, nblocks
uvel(:,:,iblk) = fld2(:,:,1,iblk)
vvel(:,:,iblk) = fld2(:,:,2,iblk)
enddo
!$OMP END PARALLEL DO

if (maskhalo_dyn) then
call ice_timer_start(timer_bound)
halomask = 0
Expand Down Expand Up @@ -442,13 +432,10 @@ subroutine evp (dt)
uvel_init(:,:,iblk), vvel_init(:,:,iblk),&
uvel (:,:,iblk), vvel (:,:,iblk), &
Tbu (:,:,iblk))

! load velocity into array for boundary updates
fld2(:,:,1,iblk) = uvel(:,:,iblk)
fld2(:,:,2,iblk) = vvel(:,:,iblk)
enddo
!$TCXOMP END PARALLEL DO

call stack_velocity_field(uvel, vvel, fld2)
call ice_timer_start(timer_bound)
if (maskhalo_dyn) then
call ice_HaloUpdate (fld2, halo_info_mask, &
Expand All @@ -458,14 +445,7 @@ subroutine evp (dt)
field_loc_NEcorner, field_type_vector)
endif
call ice_timer_stop(timer_bound)

! unload
!$OMP PARALLEL DO PRIVATE(iblk)
do iblk = 1, nblocks
uvel(:,:,iblk) = fld2(:,:,1,iblk)
vvel(:,:,iblk) = fld2(:,:,2,iblk)
enddo
!$OMP END PARALLEL DO
call unstack_velocity_field(fld2, uvel, vvel)

enddo ! subcycling
endif ! kevp_kernel
Expand Down Expand Up @@ -599,6 +579,8 @@ subroutine stress (nx_block, ny_block, &
rdg_conv, rdg_shear, &
str )

use ice_dyn_shared, only: strain_rates, deformations

integer (kind=int_kind), intent(in) :: &
nx_block, ny_block, & ! block dimensions
ksub , & ! subcycling step
Expand Down Expand Up @@ -676,58 +658,20 @@ subroutine stress (nx_block, ny_block, &
! strain rates
! NOTE these are actually strain rates * area (m^2/s)
!-----------------------------------------------------------------
! divergence = e_11 + e_22
divune = cyp(i,j)*uvel(i ,j ) - dyt(i,j)*uvel(i-1,j ) &
+ cxp(i,j)*vvel(i ,j ) - dxt(i,j)*vvel(i ,j-1)
divunw = cym(i,j)*uvel(i-1,j ) + dyt(i,j)*uvel(i ,j ) &
+ cxp(i,j)*vvel(i-1,j ) - dxt(i,j)*vvel(i-1,j-1)
divusw = cym(i,j)*uvel(i-1,j-1) + dyt(i,j)*uvel(i ,j-1) &
+ cxm(i,j)*vvel(i-1,j-1) + dxt(i,j)*vvel(i-1,j )
divuse = cyp(i,j)*uvel(i ,j-1) - dyt(i,j)*uvel(i-1,j-1) &
+ cxm(i,j)*vvel(i ,j-1) + dxt(i,j)*vvel(i ,j )

! tension strain rate = e_11 - e_22
tensionne = -cym(i,j)*uvel(i ,j ) - dyt(i,j)*uvel(i-1,j ) &
+ cxm(i,j)*vvel(i ,j ) + dxt(i,j)*vvel(i ,j-1)
tensionnw = -cyp(i,j)*uvel(i-1,j ) + dyt(i,j)*uvel(i ,j ) &
+ cxm(i,j)*vvel(i-1,j ) + dxt(i,j)*vvel(i-1,j-1)
tensionsw = -cyp(i,j)*uvel(i-1,j-1) + dyt(i,j)*uvel(i ,j-1) &
+ cxp(i,j)*vvel(i-1,j-1) - dxt(i,j)*vvel(i-1,j )
tensionse = -cym(i,j)*uvel(i ,j-1) - dyt(i,j)*uvel(i-1,j-1) &
+ cxp(i,j)*vvel(i ,j-1) - dxt(i,j)*vvel(i ,j )

! shearing strain rate = e_12
shearne = -cym(i,j)*vvel(i ,j ) - dyt(i,j)*vvel(i-1,j ) &
- cxm(i,j)*uvel(i ,j ) - dxt(i,j)*uvel(i ,j-1)
shearnw = -cyp(i,j)*vvel(i-1,j ) + dyt(i,j)*vvel(i ,j ) &
- cxm(i,j)*uvel(i-1,j ) - dxt(i,j)*uvel(i-1,j-1)
shearsw = -cyp(i,j)*vvel(i-1,j-1) + dyt(i,j)*vvel(i ,j-1) &
- cxp(i,j)*uvel(i-1,j-1) + dxt(i,j)*uvel(i-1,j )
shearse = -cym(i,j)*vvel(i ,j-1) - dyt(i,j)*vvel(i-1,j-1) &
- cxp(i,j)*uvel(i ,j-1) + dxt(i,j)*uvel(i ,j )

! Delta (in the denominator of zeta, eta)
Deltane = sqrt(divune**2 + ecci*(tensionne**2 + shearne**2))
Deltanw = sqrt(divunw**2 + ecci*(tensionnw**2 + shearnw**2))
Deltasw = sqrt(divusw**2 + ecci*(tensionsw**2 + shearsw**2))
Deltase = sqrt(divuse**2 + ecci*(tensionse**2 + shearse**2))

!-----------------------------------------------------------------
! on last subcycle, save quantities for mechanical redistribution
!-----------------------------------------------------------------
if (ksub == ndte) then
divu(i,j) = p25*(divune + divunw + divuse + divusw) * tarear(i,j)
tmp = p25*(Deltane + Deltanw + Deltase + Deltasw) * tarear(i,j)
rdg_conv(i,j) = -min(divu(i,j),c0)
rdg_shear(i,j) = p5*(tmp-abs(divu(i,j)))

! diagnostic only
! shear = sqrt(tension**2 + shearing**2)
shear(i,j) = p25*tarear(i,j)*sqrt( &
(tensionne + tensionnw + tensionse + tensionsw)**2 &
+ (shearne + shearnw + shearse + shearsw)**2)

endif
call strain_rates (nx_block, ny_block, &
i, j, &
uvel, vvel, &
dxt, dyt, &
cxp, cyp, &
cxm, cym, &
divune, divunw, &
divuse, divusw, &
tensionne, tensionnw, &
tensionse, tensionsw, &
shearne, shearnw, &
shearse, shearsw, &
Deltane, Deltanw, &
Deltase, Deltasw )

!-----------------------------------------------------------------
! strength/Delta ! kg/s
Expand Down Expand Up @@ -902,6 +846,23 @@ subroutine stress (nx_block, ny_block, &

enddo ! ij

!-----------------------------------------------------------------
! on last subcycle, save quantities for mechanical redistribution
!-----------------------------------------------------------------
if (ksub == ndte) then
call deformations (nx_block , ny_block , &
icellt , &
indxti , indxtj , &
uvel , vvel , &
dxt , dyt , &
cxp , cyp , &
cxm , cym , &
tarear , &
shear , divu , &
rdg_conv , rdg_shear )

endif

end subroutine stress

!=======================================================================
Expand Down
4 changes: 2 additions & 2 deletions cicecore/cicedynB/dynamics/ice_dyn_evp_1d.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -326,7 +326,7 @@ subroutine stress_i(NA_len, &
! tension strain rate = e_11 - e_22
tensionne = -cym*uvel(iw) - dyt(iw)*tmp_uvel_ee &
+ cxm*vvel(iw) + dxt(iw)*tmp_vvel_se
! shearing strain rate = e_12
! shearing strain rate = 2*e_12
shearne = -cym*vvel(iw) - dyt(iw)*tmp_vvel_ee &
- cxm*uvel(iw) - dxt(iw)*tmp_uvel_se
! Delta (in the denominator of zeta, eta)
Expand Down Expand Up @@ -614,7 +614,7 @@ subroutine stress_l(NA_len, tarear, &
tensionse = -cym*tmp_uvel_se - dyt(iw)*tmp_uvel_ne &
+ cxp*tmp_vvel_se - dxt(iw)*vvel(iw)

! shearing strain rate = e_12
! shearing strain rate = 2*e_12
shearne = -cym*vvel(iw) - dyt(iw)*tmp_vvel_ee &
- cxm*uvel(iw) - dxt(iw)*tmp_uvel_se
shearnw = -cyp*tmp_vvel_ee + dyt(iw)*vvel(iw) &
Expand Down
Loading

0 comments on commit f773ef3

Please sign in to comment.