Skip to content

Commit

Permalink
Clean up implementation following cgridDEV merge (#719)
Browse files Browse the repository at this point in the history 
* Clean up implementation following cgridDEV merge, based on comments
in #715

Refactor capping namelist, added capping_method variable string for namelist,
  options are "hibler" or "kreyscher" equivalent to 1 and 0 in old
  capping namelist.  capping variable still exists and is set
  depending on capping_method setting.  Update ice_in and other set_nml options.
Add box2001 option for ice_data_conc.  Does same things at p5, but makes
  decreases confusion in box2001 usage.  Updated set_nml options.
Remove support for ice_data_type = smallblock, bigblock, easthalf.  Remove
  a couple tests that were using these.
Capitalize last letter in variables dxt, dyt, dxe, dye, dxn, dyn, dxu, dyu ->
  dxT, dyT, dxE, dyE, dxN, dyN, dxU, dyU.  This is the general direction we
  want to go, but this was done to make "dyn" look less like "dynamics".
Modify some comment indentation in ice_init.F90
Reorder grid_average_X2Y calls in ice_dyn_evp.F90 for readability
Remove older code in ice_dyn_shared.F90 referencing dte, dte2T, tdamp2, etc.
Clean up comments in ice_history_shared.F90, ice_dyn_evp.F90, ice_transport_driver.F90,
  ice_flux.F90, ice_forcing.F90, ice_grid.F90
Remove commented out code in ice_dyn_shared.F90 around computation of iceumask
Remove commented out code in ice_forcing.F90 related to box2001_data_ocn
Update documentation.
  capping_method
  namelist table was reviewed and updated for consistency with code
  ice_ic change from default to internal
  dxt -> dxT, etc

* update andacc namelist options, comment out for now

* switch capping_method options to max and sum

* Add alt07 tests for new dynamics namelist options
  • Loading branch information
@apcraig
apcraig authored May 17, 2022
1 parent 078aab4 commit 5bc0a34
Show file tree
Hide file tree
Showing 40 changed files with 606 additions and 740 deletions.
10 changes: 5 additions & 5 deletions cicecore/cicedynB/analysis/ice_history.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2091,7 +2091,7 @@ subroutine accum_hist (dt)
use ice_blocks, only: block, get_block, nx_block, ny_block
use ice_domain, only: blocks_ice, nblocks
use ice_domain_size, only: nfsd
use ice_grid, only: tmask, lmask_n, lmask_s, dxu, dyu, grid_ice
use ice_grid, only: tmask, lmask_n, lmask_s, dxU, dyU, grid_ice
use ice_calendar, only: new_year, write_history, &
write_ic, timesecs, histfreq, nstreams, mmonth, &
new_month
Expand Down Expand Up @@ -2765,8 +2765,8 @@ subroutine accum_hist (dt)
do j = jlo, jhi
do i = ilo, ihi
if (aice(i,j,iblk) > puny) &
worka(i,j) = (rhoi*p5*(vice(i+1,j,iblk)+vice(i,j,iblk))*dyu(i,j,iblk) &
+ rhos*p5*(vsno(i+1,j,iblk)+vsno(i,j,iblk))*dyu(i,j,iblk)) &
worka(i,j) = (rhoi*p5*(vice(i+1,j,iblk)+vice(i,j,iblk))*dyU(i,j,iblk) &
+ rhos*p5*(vsno(i+1,j,iblk)+vsno(i,j,iblk))*dyU(i,j,iblk)) &
* p5*(uvel(i,j-1,iblk)+uvel(i,j,iblk))
enddo
enddo
Expand All @@ -2778,8 +2778,8 @@ subroutine accum_hist (dt)
do j = jlo, jhi
do i = ilo, ihi
if (aice(i,j,iblk) > puny) &
worka(i,j) = (rhoi*p5*(vice(i,j+1,iblk)+vice(i,j,iblk))*dxu(i,j,iblk) &
+ rhos*p5*(vsno(i,j+1,iblk)+vsno(i,j,iblk))*dxu(i,j,iblk)) &
worka(i,j) = (rhoi*p5*(vice(i,j+1,iblk)+vice(i,j,iblk))*dxU(i,j,iblk) &
+ rhos*p5*(vsno(i,j+1,iblk)+vsno(i,j,iblk))*dxU(i,j,iblk)) &
* p5*(vvel(i-1,j,iblk)+vvel(i,j,iblk))
enddo
enddo
Expand Down
2 changes: 1 addition & 1 deletion cicecore/cicedynB/analysis/ice_history_shared.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -402,7 +402,7 @@ module ice_history_shared
f_Tsfc, f_aice , &
f_uvel, f_vvel , &
f_icespd, f_icedir , &
! For now, don't allow the users to modify the CD grid quantities.
! For now, C and CD grid quantities are controlled by the generic (originally B-grid) namelist flag
! f_uvelE, f_vvelE , &
! f_icespdE, f_icedirE , &
! f_uvelN, f_vvelN , &
Expand Down
82 changes: 41 additions & 41 deletions cicecore/cicedynB/dynamics/ice_dyn_eap.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -143,7 +143,7 @@ subroutine eap (dt)
stressp_1, stressp_2, stressp_3, stressp_4, &
stressm_1, stressm_2, stressm_3, stressm_4, &
stress12_1, stress12_2, stress12_3, stress12_4
use ice_grid, only: tmask, umask, dxt, dyt, dxhy, dyhx, cxp, cyp, cxm, cym, &
use ice_grid, only: tmask, umask, dxT, dyT, dxhy, dyhx, cxp, cyp, cxm, cym, &
tarear, uarear, grid_average_X2Y, &
grid_atm_dynu, grid_atm_dynv, grid_ocn_dynu, grid_ocn_dynv
use ice_state, only: aice, vice, vsno, uvel, vvel, divu, shear, &
Expand Down Expand Up @@ -453,7 +453,7 @@ subroutine eap (dt)
indxti (:,iblk), indxtj (:,iblk), &
arlx1i, denom1, &
uvel (:,:,iblk), vvel (:,:,iblk), &
dxt (:,:,iblk), dyt (:,:,iblk), &
dxT (:,:,iblk), dyT (:,:,iblk), &
dxhy (:,:,iblk), dyhx (:,:,iblk), &
cxp (:,:,iblk), cyp (:,:,iblk), &
cxm (:,:,iblk), cym (:,:,iblk), &
Expand Down Expand Up @@ -1165,7 +1165,7 @@ subroutine stress_eap (nx_block, ny_block, &
indxti, indxtj, &
arlx1i, denom1, &
uvel, vvel, &
dxt, dyt, &
dxT, dyT, &
dxhy, dyhx, &
cxp, cyp, &
cxm, cym, &
Expand Down Expand Up @@ -1208,8 +1208,8 @@ subroutine stress_eap (nx_block, ny_block, &
strength , & ! ice strength (N/m)
uvel , & ! x-component of velocity (m/s)
vvel , & ! y-component of velocity (m/s)
dxt , & ! width of T-cell through the middle (m)
dyt , & ! height of T-cell through the middle (m)
dxT , & ! width of T-cell through the middle (m)
dyT , & ! height of T-cell through the middle (m)
dxhy , & ! 0.5*(HTE - HTW)
dyhx , & ! 0.5*(HTN - HTS)
cyp , & ! 1.5*HTE - 0.5*HTW
Expand Down Expand Up @@ -1291,34 +1291,34 @@ subroutine stress_eap (nx_block, ny_block, &
!-----------------------------------------------------------------

! 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 )
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 )
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 = 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 ) &
- 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 )
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 )

!-----------------------------------------------------------------
! Stress updated depending on strain rate and structure tensor
Expand Down Expand Up @@ -1499,17 +1499,17 @@ subroutine stress_eap (nx_block, ny_block, &
csig12se = p222*stress12_4(i,j) + ssig121 &
+ p055*stress12_2(i,j)

str12ew = p5*dxt(i,j)*(p333*ssig12e + p166*ssig12w)
str12we = p5*dxt(i,j)*(p333*ssig12w + p166*ssig12e)
str12ns = p5*dyt(i,j)*(p333*ssig12n + p166*ssig12s)
str12sn = p5*dyt(i,j)*(p333*ssig12s + p166*ssig12n)
str12ew = p5*dxT(i,j)*(p333*ssig12e + p166*ssig12w)
str12we = p5*dxT(i,j)*(p333*ssig12w + p166*ssig12e)
str12ns = p5*dyT(i,j)*(p333*ssig12n + p166*ssig12s)
str12sn = p5*dyT(i,j)*(p333*ssig12s + p166*ssig12n)

!-----------------------------------------------------------------
! for dF/dx (u momentum)
!-----------------------------------------------------------------

strp_tmp = p25*dyt(i,j)*(p333*ssigpn + p166*ssigps)
strm_tmp = p25*dyt(i,j)*(p333*ssigmn + p166*ssigms)
strp_tmp = p25*dyT(i,j)*(p333*ssigpn + p166*ssigps)
strm_tmp = p25*dyT(i,j)*(p333*ssigmn + p166*ssigms)

! northeast (i,j)
strtmp(i,j,1) = -strp_tmp - strm_tmp - str12ew &
Expand All @@ -1519,8 +1519,8 @@ subroutine stress_eap (nx_block, ny_block, &
strtmp(i,j,2) = strp_tmp + strm_tmp - str12we &
+ dxhy(i,j)*(-csigpnw + csigmnw) + dyhx(i,j)*csig12nw

strp_tmp = p25*dyt(i,j)*(p333*ssigps + p166*ssigpn)
strm_tmp = p25*dyt(i,j)*(p333*ssigms + p166*ssigmn)
strp_tmp = p25*dyT(i,j)*(p333*ssigps + p166*ssigpn)
strm_tmp = p25*dyT(i,j)*(p333*ssigms + p166*ssigmn)

! southeast (i,j+1)
strtmp(i,j,3) = -strp_tmp - strm_tmp + str12ew &
Expand All @@ -1534,8 +1534,8 @@ subroutine stress_eap (nx_block, ny_block, &
! for dF/dy (v momentum)
!-----------------------------------------------------------------

strp_tmp = p25*dxt(i,j)*(p333*ssigpe + p166*ssigpw)
strm_tmp = p25*dxt(i,j)*(p333*ssigme + p166*ssigmw)
strp_tmp = p25*dxT(i,j)*(p333*ssigpe + p166*ssigpw)
strm_tmp = p25*dxT(i,j)*(p333*ssigme + p166*ssigmw)

! northeast (i,j)
strtmp(i,j,5) = -strp_tmp + strm_tmp - str12ns &
Expand All @@ -1545,8 +1545,8 @@ subroutine stress_eap (nx_block, ny_block, &
strtmp(i,j,6) = strp_tmp - strm_tmp - str12sn &
- dyhx(i,j)*(csigpse + csigmse) + dxhy(i,j)*csig12se

strp_tmp = p25*dxt(i,j)*(p333*ssigpw + p166*ssigpe)
strm_tmp = p25*dxt(i,j)*(p333*ssigmw + p166*ssigme)
strp_tmp = p25*dxT(i,j)*(p333*ssigpw + p166*ssigpe)
strm_tmp = p25*dxT(i,j)*(p333*ssigmw + p166*ssigme)

! northwest (i+1,j)
strtmp(i,j,7) = -strp_tmp + strm_tmp + str12ns &
Expand Down
50 changes: 25 additions & 25 deletions cicecore/cicedynB/dynamics/ice_dyn_evp.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -99,7 +99,7 @@ subroutine evp (dt)
stresspT, stressmT, stress12T, &
stresspU, stressmU, stress12U
use ice_grid, only: hm, tmask, umask, umaskCD, nmask, emask, uvm, epm, npm, &
dxe, dxn, dxt, dxu, dye, dyn, dyt, dyu, &
dxE, dxN, dxT, dxU, dyE, dyN, dyT, dyU, &
ratiodxN, ratiodxNr, ratiodyE, ratiodyEr, &
dxhy, dyhx, cxp, cyp, cxm, cym, &
tarear, uarear, earear, narear, grid_average_X2Y, tarea, uarea, &
Expand Down Expand Up @@ -185,9 +185,9 @@ subroutine evp (dt)
emassdti ! mass of E-cell/dte (kg/m^2 s)

real (kind=dbl_kind), allocatable :: &
fld2(:,:,:,:) , & ! bundled fields size 2
fld3(:,:,:,:) , & ! bundled fields size 3
fld4(:,:,:,:) ! bundled fields size 4
fld2(:,:,:,:) , & ! 2 bundled fields
fld3(:,:,:,:) , & ! 3 bundled fields
fld4(:,:,:,:) ! 4 bundled fields

real (kind=dbl_kind), allocatable :: &
strengthU(:,:,:), & ! strength averaged to U points
Expand Down Expand Up @@ -326,12 +326,12 @@ subroutine evp (dt)
if (grid_ice == 'CD' .or. grid_ice == 'C') then
call grid_average_X2Y('F', tmass , 'T' , emass, 'E')
call grid_average_X2Y('F', aice_init, 'T' , aie , 'E')
call grid_average_X2Y('S', uocn , grid_ocn_dynu, uocnE, 'E')
call grid_average_X2Y('S', vocn , grid_ocn_dynv, vocnE, 'E')
call grid_average_X2Y('F', tmass , 'T' , nmass, 'N')
call grid_average_X2Y('F', aice_init, 'T' , ain , 'N')
call grid_average_X2Y('S', uocn , grid_ocn_dynu, uocnN, 'N')
call grid_average_X2Y('S', vocn , grid_ocn_dynv, vocnN, 'N')
call grid_average_X2Y('S', uocn , grid_ocn_dynu, uocnE, 'E')
call grid_average_X2Y('S', vocn , grid_ocn_dynv, vocnE, 'E')
endif
!----------------------------------------------------------------
! Set wind stress to values supplied via NEMO or other forcing
Expand Down Expand Up @@ -722,7 +722,7 @@ subroutine evp (dt)
icetmask, iceumask, &
cdn_ocn,aiu,uocnU,vocnU,forcex,forcey,Tbu, &
umassdti,fm,uarear,tarear,strintx,strinty,uvel_init,vvel_init,&
strength,uvel,vvel,dxt,dyt, &
strength,uvel,vvel,dxT,dyT, &
stressp_1 ,stressp_2, stressp_3, stressp_4, &
stressm_1 ,stressm_2, stressm_3, stressm_4, &
stress12_1,stress12_2,stress12_3,stress12_4 )
Expand Down Expand Up @@ -754,7 +754,7 @@ subroutine evp (dt)
icellt (iblk), &
indxti (:,iblk), indxtj (:,iblk), &
uvel (:,:,iblk), vvel (:,:,iblk), &
dxt (:,:,iblk), dyt (:,:,iblk), &
dxT (:,:,iblk), dyT (:,:,iblk), &
dxhy (:,:,iblk), dyhx (:,:,iblk), &
cxp (:,:,iblk), cyp (:,:,iblk), &
cxm (:,:,iblk), cym (:,:,iblk), &
Expand All @@ -776,7 +776,7 @@ subroutine evp (dt)
icellt (iblk), &
indxti (:,iblk), indxtj (:,iblk), &
uvel (:,:,iblk), vvel (:,:,iblk), &
dxt (:,:,iblk), dyt (:,:,iblk), &
dxT (:,:,iblk), dyT (:,:,iblk), &
cxp (:,:,iblk), cyp (:,:,iblk), &
cxm (:,:,iblk), cym (:,:,iblk), &
tarear (:,:,iblk), &
Expand Down Expand Up @@ -1335,7 +1335,7 @@ subroutine stress (nx_block, ny_block, &
icellt, &
indxti, indxtj, &
uvel, vvel, &
dxt, dyt, &
dxT, dyT, &
dxhy, dyhx, &
cxp, cyp, &
cxm, cym, &
Expand Down Expand Up @@ -1363,8 +1363,8 @@ subroutine stress (nx_block, ny_block, &
strength , & ! ice strength (N/m)
uvel , & ! x-component of velocity (m/s)
vvel , & ! y-component of velocity (m/s)
dxt , & ! width of T-cell through the middle (m)
dyt , & ! height of T-cell through the middle (m)
dxT , & ! width of T-cell through the middle (m)
dyT , & ! height of T-cell through the middle (m)
dxhy , & ! 0.5*(HTE - HTW)
dyhx , & ! 0.5*(HTN - HTS)
cyp , & ! 1.5*HTE - 0.5*HTW
Expand Down Expand Up @@ -1425,7 +1425,7 @@ subroutine stress (nx_block, ny_block, &
call strain_rates (nx_block, ny_block, &
i, j, &
uvel, vvel, &
dxt, dyt, &
dxT, dyT, &
cxp, cyp, &
cxm, cym, &
divune, divunw, &
Expand Down Expand Up @@ -1561,16 +1561,16 @@ subroutine stress (nx_block, ny_block, &
csig12se = p222*stress12_4(i,j) + ssig121 &
+ p055*stress12_2(i,j)

str12ew = p5*dxt(i,j)*(p333*ssig12e + p166*ssig12w)
str12we = p5*dxt(i,j)*(p333*ssig12w + p166*ssig12e)
str12ns = p5*dyt(i,j)*(p333*ssig12n + p166*ssig12s)
str12sn = p5*dyt(i,j)*(p333*ssig12s + p166*ssig12n)
str12ew = p5*dxT(i,j)*(p333*ssig12e + p166*ssig12w)
str12we = p5*dxT(i,j)*(p333*ssig12w + p166*ssig12e)
str12ns = p5*dyT(i,j)*(p333*ssig12n + p166*ssig12s)
str12sn = p5*dyT(i,j)*(p333*ssig12s + p166*ssig12n)

!-----------------------------------------------------------------
! for dF/dx (u momentum)
!-----------------------------------------------------------------
strp_tmp = p25*dyt(i,j)*(p333*ssigpn + p166*ssigps)
strm_tmp = p25*dyt(i,j)*(p333*ssigmn + p166*ssigms)
strp_tmp = p25*dyT(i,j)*(p333*ssigpn + p166*ssigps)
strm_tmp = p25*dyT(i,j)*(p333*ssigmn + p166*ssigms)

! northeast (i,j)
str(i,j,1) = -strp_tmp - strm_tmp - str12ew &
Expand All @@ -1580,8 +1580,8 @@ subroutine stress (nx_block, ny_block, &
str(i,j,2) = strp_tmp + strm_tmp - str12we &
+ dxhy(i,j)*(-csigpnw + csigmnw) + dyhx(i,j)*csig12nw

strp_tmp = p25*dyt(i,j)*(p333*ssigps + p166*ssigpn)
strm_tmp = p25*dyt(i,j)*(p333*ssigms + p166*ssigmn)
strp_tmp = p25*dyT(i,j)*(p333*ssigps + p166*ssigpn)
strm_tmp = p25*dyT(i,j)*(p333*ssigms + p166*ssigmn)

! southeast (i,j+1)
str(i,j,3) = -strp_tmp - strm_tmp + str12ew &
Expand All @@ -1594,8 +1594,8 @@ subroutine stress (nx_block, ny_block, &
!-----------------------------------------------------------------
! for dF/dy (v momentum)
!-----------------------------------------------------------------
strp_tmp = p25*dxt(i,j)*(p333*ssigpe + p166*ssigpw)
strm_tmp = p25*dxt(i,j)*(p333*ssigme + p166*ssigmw)
strp_tmp = p25*dxT(i,j)*(p333*ssigpe + p166*ssigpw)
strm_tmp = p25*dxT(i,j)*(p333*ssigme + p166*ssigmw)

! northeast (i,j)
str(i,j,5) = -strp_tmp + strm_tmp - str12ns &
Expand All @@ -1605,8 +1605,8 @@ subroutine stress (nx_block, ny_block, &
str(i,j,6) = strp_tmp - strm_tmp - str12sn &
- dyhx(i,j)*(csigpse + csigmse) + dxhy(i,j)*csig12se

strp_tmp = p25*dxt(i,j)*(p333*ssigpw + p166*ssigpe)
strm_tmp = p25*dxt(i,j)*(p333*ssigmw + p166*ssigme)
strp_tmp = p25*dxT(i,j)*(p333*ssigpw + p166*ssigpe)
strm_tmp = p25*dxT(i,j)*(p333*ssigmw + p166*ssigme)

! northwest (i+1,j)
str(i,j,7) = -strp_tmp + strm_tmp + str12ns &
Expand Down
Loading

0 comments on commit 5bc0a34

Please sign in to comment.