Skip to content

umaskCD and more #33

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
apcraig merged 15 commits into from
Nov 24, 2021
Merged

umaskCD and more #33

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
cf65c39
Changes in ice_dyn* are interpolation of uvelE/vvelN to B grid. ice_t…
TillRasmussen Nov 17, 2021
89e401a
changed according to comments. changed average from F to S and. comme…
TillRasmussen Nov 17, 2021
24d98db
Changes in ice_dyn* are interpolation of uvelE/vvelN to B grid. ice_t…
TillRasmussen Nov 17, 2021
8849523
changed according to comments. changed average from F to S and. comme…
TillRasmussen Nov 17, 2021
ae40964
Create local uvmCD that is projected into umaskCD and icellu.
TillRasmussen Nov 19, 2021
3855b8c
Added call to routine
TillRasmussen Nov 19, 2021
004aa9c
bux fixes
TillRasmussen Nov 19, 2021
fcd291b
last bugfixes. Code compiles and is bit for bit with code before.
TillRasmussen Nov 20, 2021
710967f
First move towards merge with merge of commit 39
TillRasmussen Nov 21, 2021
9b1ec14
update ice_dyn_evp
TillRasmussen Nov 21, 2021
f5d7d89
It is now bit for bit. Number of mask is reduced from 2 to one.
TillRasmussen Nov 21, 2021
15c16f5
Merge branch 'cgridDEV' into cgrid-something
TillRasmussen Nov 21, 2021
65d46c0
moved trim
TillRasmussen Nov 23, 2021
82c90b2
Based on Quicksuite: B grid bit for bit. CD grid passes 13 of 16. 2 f…
TillRasmussen Nov 23, 2021
de8f6e9
Merge branch 'cgrid-something' of github.com:TillRasmussen/CICE into …
TillRasmussen Nov 24, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
205 changes: 122 additions & 83 deletions cicecore/cicedynB/dynamics/ice_dyn_evp.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -98,11 +98,11 @@ subroutine evp (dt)
stress12_1, stress12_2, stress12_3, stress12_4, &
stresspT, stressmT, stress12T, &
stresspU, stressmU, stress12U
use ice_grid, only: tmask, umask, nmask, emask, uvm, epm, npm, &
use ice_grid, only: hm, tmask, umask, umaskCD, nmask, emask, uvm, epm, npm, &
dxe, dxn, dxt, dxu, dye, dyn, dyt, dyu, &
ratiodxN, ratiodxNr, ratiodyE, ratiodyEr, &
dxhy, dyhx, cxp, cyp, cxm, cym, &
tarear, uarear, earea, narea, tinyarea, grid_average_X2Y, &
tarear, uarear, earea, narea, tinyarea, grid_average_X2Y, tarea, &
grid_type, grid_system
use ice_state, only: aice, vice, vsno, uvel, vvel, uvelN, vvelN, &
uvelE, vvelE, divu, shear, &
Expand Down Expand Up @@ -341,35 +341,67 @@ subroutine evp (dt)
ihi = this_block%ihi
jlo = this_block%jlo
jhi = this_block%jhi

call dyn_prep2 (nx_block, ny_block, &
ilo, ihi, jlo, jhi, &
icellt(iblk), icellu(iblk), &
indxti (:,iblk), indxtj (:,iblk), &
indxui (:,iblk), indxuj (:,iblk), &
aiu (:,:,iblk), umass (:,:,iblk), &
umassdti (:,:,iblk), fcor_blk (:,:,iblk), &
umask (:,:,iblk), &
uocn (:,:,iblk), vocn (:,:,iblk), &
strairx (:,:,iblk), strairy (:,:,iblk), &
ss_tltx (:,:,iblk), ss_tlty (:,:,iblk), &
icetmask (:,:,iblk), iceumask (:,:,iblk), &
fm (:,:,iblk), dt, &
strtltx (:,:,iblk), strtlty (:,:,iblk), &
strocnx (:,:,iblk), strocny (:,:,iblk), &
strintx (:,:,iblk), strinty (:,:,iblk), &
taubx (:,:,iblk), tauby (:,:,iblk), &
waterx (:,:,iblk), watery (:,:,iblk), &
forcex (:,:,iblk), forcey (:,:,iblk), &
stressp_1 (:,:,iblk), stressp_2 (:,:,iblk), &
stressp_3 (:,:,iblk), stressp_4 (:,:,iblk), &
stressm_1 (:,:,iblk), stressm_2 (:,:,iblk), &
stressm_3 (:,:,iblk), stressm_4 (:,:,iblk), &
stress12_1(:,:,iblk), stress12_2(:,:,iblk), &
stress12_3(:,:,iblk), stress12_4(:,:,iblk), &
uvel_init (:,:,iblk), vvel_init (:,:,iblk), &
uvel (:,:,iblk), vvel (:,:,iblk), &
Tbu (:,:,iblk))

if (trim(grid_system) == 'B') then
call dyn_prep2 (nx_block, ny_block, &
ilo, ihi, jlo, jhi, &
icellt(iblk), icellu(iblk), &
indxti (:,iblk), indxtj (:,iblk), &
indxui (:,iblk), indxuj (:,iblk), &
aiu (:,:,iblk), umass (:,:,iblk), &
umassdti (:,:,iblk), fcor_blk (:,:,iblk), &
umask (:,:,iblk), &
uocn (:,:,iblk), vocn (:,:,iblk), &
strairx (:,:,iblk), strairy (:,:,iblk), &
ss_tltx (:,:,iblk), ss_tlty (:,:,iblk), &
icetmask (:,:,iblk), iceumask (:,:,iblk), &
fm (:,:,iblk), dt, &
strtltx (:,:,iblk), strtlty (:,:,iblk), &
strocnx (:,:,iblk), strocny (:,:,iblk), &
strintx (:,:,iblk), strinty (:,:,iblk), &
taubx (:,:,iblk), tauby (:,:,iblk), &
waterx (:,:,iblk), watery (:,:,iblk), &
forcex (:,:,iblk), forcey (:,:,iblk), &
stressp_1 (:,:,iblk), stressp_2 (:,:,iblk), &
stressp_3 (:,:,iblk), stressp_4 (:,:,iblk), &
stressm_1 (:,:,iblk), stressm_2 (:,:,iblk), &
stressm_3 (:,:,iblk), stressm_4 (:,:,iblk), &
stress12_1(:,:,iblk), stress12_2(:,:,iblk), &
stress12_3(:,:,iblk), stress12_4(:,:,iblk), &
uvel_init (:,:,iblk), vvel_init (:,:,iblk), &
uvel (:,:,iblk), vvel (:,:,iblk), &
Tbu (:,:,iblk))

elseif (trim(grid_system) == 'CD') then
call dyn_prep2 (nx_block, ny_block, &
ilo, ihi, jlo, jhi, &
icellt(iblk), icellu(iblk), &
indxti (:,iblk), indxtj (:,iblk), &
indxui (:,iblk), indxuj (:,iblk), &
aiu (:,:,iblk), umass (:,:,iblk), &
umassdti (:,:,iblk), fcor_blk (:,:,iblk), &
umaskCD (:,:,iblk), &
Copy link
Copy Markdown
Collaborator

@phil-blain phil-blain Nov 23, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's unclear to me why we need yet another call to dyn_prep2 here. I thought that all CD grid quantities were being dealt with in the two calls to dyn_prep2 below (for 'E' and 'N' quantities, line 440 and 473 in this version). Here your are adding a new call, identical ( as far as I can tell) to the B-grid case, apart from this new umaskCD array instead of umask... But on the CD grid, we are not using most of these arrays that are sent as arguments here...

Copy link
Copy Markdown
Author

@TillRasmussen TillRasmussen Nov 23, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

umaskCD is defined different, icellu will be different (stress_u). It can be that calling dyn_prep2 is too much and a call to a simpler function could be used.

uocn (:,:,iblk), vocn (:,:,iblk), &
strairx (:,:,iblk), strairy (:,:,iblk), &
ss_tltx (:,:,iblk), ss_tlty (:,:,iblk), &
icetmask (:,:,iblk), iceumask (:,:,iblk), &
fm (:,:,iblk), dt, &
strtltx (:,:,iblk), strtlty (:,:,iblk), &
strocnx (:,:,iblk), strocny (:,:,iblk), &
strintx (:,:,iblk), strinty (:,:,iblk), &
taubx (:,:,iblk), tauby (:,:,iblk), &
waterx (:,:,iblk), watery (:,:,iblk), &
forcex (:,:,iblk), forcey (:,:,iblk), &
stressp_1 (:,:,iblk), stressp_2 (:,:,iblk), &
stressp_3 (:,:,iblk), stressp_4 (:,:,iblk), &
stressm_1 (:,:,iblk), stressm_2 (:,:,iblk), &
stressm_3 (:,:,iblk), stressm_4 (:,:,iblk), &
stress12_1(:,:,iblk), stress12_2(:,:,iblk), &
stress12_3(:,:,iblk), stress12_4(:,:,iblk), &
uvel_init (:,:,iblk), vvel_init (:,:,iblk), &
uvel (:,:,iblk), vvel (:,:,iblk), &
Tbu (:,:,iblk))
endif

!-----------------------------------------------------------------
! ice strength
Expand Down Expand Up @@ -412,8 +444,8 @@ subroutine evp (dt)
indxti (:,iblk), indxtj (:,iblk), &
indxni (:,iblk), indxnj (:,iblk), &
aiN (:,:,iblk), nmass (:,:,iblk), &
nmassdti (:,:,iblk), fcorN_blk (:,:,iblk), &
nmask (:,:,iblk), &
nmassdti (:,:,iblk), fcorN_blk (:,:,iblk),&
nmask (:,:,iblk), &
uocnN (:,:,iblk), vocnN (:,:,iblk), &
strairxN (:,:,iblk), strairyN (:,:,iblk), &
ss_tltxN (:,:,iblk), ss_tltyN (:,:,iblk), &
Expand Down Expand Up @@ -445,7 +477,7 @@ subroutine evp (dt)
indxti (:,iblk), indxtj (:,iblk), &
indxei (:,iblk), indxej (:,iblk), &
aiE (:,:,iblk), emass (:,:,iblk), &
emassdti (:,:,iblk), fcorE_blk (:,:,iblk), &
emassdti (:,:,iblk), fcorE_blk (:,:,iblk),&
emask (:,:,iblk), &
uocnE (:,:,iblk), vocnE (:,:,iblk), &
strairxE (:,:,iblk), strairyE (:,:,iblk), &
Expand Down Expand Up @@ -643,10 +675,9 @@ subroutine evp (dt)
rdg_conv (:,:,iblk), rdg_shear (:,:,iblk), &
strtmp (:,:,:) )

!-----------------------------------------------------------------
! momentum equation
!-----------------------------------------------------------------

!-----------------------------------------------------------------
! momentum equation
!-----------------------------------------------------------------
call stepu (nx_block, ny_block, &
icellu (iblk), Cdn_ocn (:,:,iblk), &
indxui (:,iblk), indxuj (:,iblk), &
Expand Down Expand Up @@ -688,13 +719,14 @@ subroutine evp (dt)
uvel (:,:,iblk), vvel (:,:,iblk), &
dxE (:,:,iblk), dyN (:,:,iblk), &
dxU (:,:,iblk), dyU (:,:,iblk), &
tarea (:,:,iblk), &
ratiodxN (:,:,iblk), ratiodxNr (:,:,iblk), &
ratiodyE (:,:,iblk), ratiodyEr (:,:,iblk), &
epm (:,:,iblk), npm (:,:,iblk), &
uvm (:,:,iblk), &
hm (:,:,iblk), uvm (:,:,iblk), &
zetax2T (:,:,iblk), etax2T (:,:,iblk), &
stresspU (:,:,iblk), stressmU (:,:,iblk), &
stress12U (:,:,iblk))
stress12U (:,:,iblk))

call div_stress (nx_block, ny_block, & ! E point
ksub, icelle(iblk), &
Expand Down Expand Up @@ -977,7 +1009,7 @@ subroutine stress (nx_block, ny_block, &
rdg_conv, rdg_shear, &
str )

use ice_dyn_shared, only: strain_rates, deformations, viscous_coeffs_and_rep_pressure
use ice_dyn_shared, only: strain_rates, deformations, viscous_coeffs_and_rep_pressure_T

integer (kind=int_kind), intent(in) :: &
nx_block, ny_block, & ! block dimensions
Expand Down Expand Up @@ -1041,16 +1073,15 @@ subroutine stress (nx_block, ny_block, &
str12ew, str12we, str12ns, str12sn , &
strp_tmp, strm_tmp, tmp

logical :: capping ! of the viscous coef
real(kind=dbl_kind),parameter :: capping = c1 ! of the viscous coef

character(len=*), parameter :: subname = '(stress)'

!-----------------------------------------------------------------
! Initialize
!-----------------------------------------------------------------

str(:,:,:) = c0
capping = .true. ! could be later included in ice_in

do ij = 1, icellt
i = indxti(ij)
Expand Down Expand Up @@ -1079,17 +1110,27 @@ subroutine stress (nx_block, ny_block, &
!-----------------------------------------------------------------
! viscous coefficients and replacement pressure
!-----------------------------------------------------------------

call viscous_coeffs_and_rep_pressure (strength(i,j), tinyarea(i,j),&
Deltane, Deltanw, &
Deltasw, Deltase, &
zetax2ne, zetax2nw, &
zetax2sw, zetax2se, &
etax2ne, etax2nw, &
etax2sw, etax2se, &
rep_prsne, rep_prsnw, &
rep_prssw, rep_prsse, &
capping)

call viscous_coeffs_and_rep_pressure_T (strength(i,j), tinyarea(i,j),&
Deltane, zetax2ne, &
etax2ne, rep_prsne, &
capping)

call viscous_coeffs_and_rep_pressure_T (strength(i,j), tinyarea(i,j),&
Deltanw, zetax2nw, &
etax2nw, rep_prsnw, &
capping)

call viscous_coeffs_and_rep_pressure_T (strength(i,j), tinyarea(i,j),&
Deltasw, zetax2sw, &
etax2sw, rep_prssw, &
capping)

call viscous_coeffs_and_rep_pressure_T (strength(i,j), tinyarea(i,j),&
Deltase, zetax2se, &
etax2se, rep_prsse, &
capping)

Comment on lines -1082 to +1133
Copy link
Copy Markdown
Collaborator

@phil-blain phil-blain Nov 23, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That is not strictly necessary, right ? it's just cleaning up the B-grid implementation. So maybe it could be done in a later step ? If it's bit for bit, then all for the better and I guess it does not matter much...


!-----------------------------------------------------------------
! the stresses ! kg/s^2
Expand Down Expand Up @@ -1345,15 +1386,14 @@ subroutine stress_T (nx_block, ny_block, &
divT, tensionT, shearT, DeltaT, & ! strain rates at T point
rep_prsT ! replacement pressure at T point

logical :: capping ! of the viscous coef
real(kind=dbl_kind), parameter :: capping = c1 ! of the viscous coef

character(len=*), parameter :: subname = '(stress_T)'

!-----------------------------------------------------------------
! Initialize
!-----------------------------------------------------------------

capping = .true. ! could be later included in ice_in

do ij = 1, icellt
i = indxti(ij)
Expand Down Expand Up @@ -1435,15 +1475,16 @@ subroutine stress_U (nx_block, ny_block, &
uvelU, vvelU, &
dxE, dyN, &
dxU, dyU, &
tarea, &
ratiodxN, ratiodxNr, &
ratiodyE, ratiodyEr, &
epm, npm, uvm, &
epm, npm, hm, uvm, &
Comment on lines -1440 to +1481
Copy link
Copy Markdown
Collaborator

@phil-blain phil-blain Nov 23, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it would be preferable to restrict ourselves to 2 arguments per line, for better readability.

zetax2T, etax2T, &
stresspU, stressmU, &
stress12U )
stress12U )

use ice_dyn_shared, only: strain_rates_U!, &
! viscous_coeffs_and_rep_pressure_U
use ice_dyn_shared, only: strain_rates_U, &
viscous_coeffs_and_rep_pressure_T2U

integer (kind=int_kind), intent(in) :: &
nx_block, ny_block, & ! block dimensions
Expand All @@ -1456,21 +1497,23 @@ subroutine stress_U (nx_block, ny_block, &

real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
uvelE , & ! x-component of velocity (m/s) at the E point
vvelE , & ! y-component of velocity (m/s) at the N point
uvelN , & ! x-component of velocity (m/s) at the E point
vvelE , & ! y-component of velocity (m/s) at the E point
uvelN , & ! x-component of velocity (m/s) at the N point
vvelN , & ! y-component of velocity (m/s) at the N point
uvelU , & ! x-component of velocity (m/s) at the U point
vvelU , & ! y-component of velocity (m/s) at the U point
dxE , & ! width of E-cell through the middle (m)
dxE , & ! width of E-cell through the middle (m)
dyN , & ! height of N-cell through the middle (m)
dxU , & ! width of U-cell through the middle (m)
dxU , & ! width of U-cell through the middle (m)
dyU , & ! height of U-cell through the middle (m)
ratiodxN , & ! -dxN(i+1,j)/dxN(i,j) for BCs
ratiodxNr, & ! -dxN(i,j)/dxN(i+1,j) for BCs
ratiodyE , & ! -dyE(i,j+1)/dyE(i,j) for BCs
ratiodyEr, & ! -dyE(i,j)/dyE(i,j+1) for BCs
tarea , & ! area of T-cell (m^2)
ratiodxN , & ! -dxN(i+1,j)/dxN(i,j) factor for BCs across coastline
ratiodxNr, & ! -dxN(i,j)/dxN(i+1,j) factor for BCs across coastline
ratiodyE , & ! -dyE(i,j+1)/dyE(i,j) factor for BCs across coastline
ratiodyEr, & ! -dyE(i,j)/dyE(i,j+1) factor for BCs across coastline
epm , & ! E-cell mask
npm , & ! E-cell mask
hm , & ! T-cell mask
uvm , & ! U-cell mask
zetax2T , & ! 2*zeta at the T point
etax2T ! 2*eta at the T point
Expand Down Expand Up @@ -1517,19 +1560,15 @@ subroutine stress_U (nx_block, ny_block, &
! viscous coefficients and replacement pressure at T point
!-----------------------------------------------------------------

! COMING SOON!!!
zetax2U = c0
etax2U = c0
rep_prsU = c0
Copy link
Copy Markdown
Collaborator

@eclare108213 eclare108213 Nov 22, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It looks like these initializations were lost in the new code modifications. Or were they moved somewhere else?

Copy link
Copy Markdown
Author

@TillRasmussen TillRasmussen Nov 22, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

they dont need to be initialized. They are calculated in the function after as intent out.

Copy link
Copy Markdown
Owner

@apcraig apcraig Nov 22, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These lines were added as a temporary fix when viscous_coeffs_and_rep_pressure_U was not called. Now it is. I think this change is correct in the big picture.

Copy link
Copy Markdown
Collaborator

@eclare108213 eclare108213 Nov 23, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ok

! call viscous_coeffs_and_rep_pressure_U (zetax2T(i,j), zetax2T(i,j+1), &
! zetax2T(i+1,j+1),zetax2T(i+1,j), &
! etax2T(i,j), etax2T(i,j+1), &
! etax2T(i+1,j+1), etax2T(i+1,j), &
! hm(i,j), hm(i,j+1), &
! hm(i+1,j+1), hm(i+1,j), &
! tarea(i,j), tarea(i,j+1), &
! tarea(i+1,j+1), tarea(i+1,j), &
! DeltaU )
call viscous_coeffs_and_rep_pressure_T2U (zetax2T(i ,j ), zetax2T(i ,j+1), &
zetax2T(i+1,j+1), zetax2T(i+1,j ), &
etax2T (i ,j ), etax2T (i ,j+1), &
etax2T (i+1,j+1), etax2T (i+1,j ), &
hm (i ,j ), hm (i ,j+1), &
hm (i+1,j+1), hm (i+1,j ), &
tarea (i ,j ), tarea (i ,j+1), &
tarea (i+1,j+1), tarea (i+1,j ), &
DeltaU,zetax2U, etax2U, rep_prsU)

!-----------------------------------------------------------------
! the stresses ! kg/s^2
Expand Down
Loading