Skip to content
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

Update grid averaging for tmass, aice, uvelT, vvelT #762

Merged
merged 2 commits into from
Sep 29, 2022
Merged

Update grid averaging for tmass, aice, uvelT, vvelT #762

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
8d52c4e
Update grid averaging for tmass, aice, uvelT, vvelT
apcraig Sep 22, 2022
d2f04bd
Refactor uvelT, vvelT implementation
apcraig Sep 23, 2022
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
64 changes: 16 additions & 48 deletions cicecore/cicedynB/analysis/ice_history_shared.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -156,57 +156,25 @@ module ice_history_shared
igrdz(nvar_grdz) ! true if category/vertical grid field is written

character (len=25), public, parameter :: &
tcstr = 'area: tarea' , & ! vcellmeas for T cell quantities
ucstr = 'area: uarea' , & ! vcellmeas for U cell quantities
ncstr = 'area: narea' , & ! vcellmeas for N cell quantities
ecstr = 'area: earea' , & ! vcellmeas for E cell quantities
tstr2D = 'TLON TLAT time' , & ! vcoord for T cell quantities, 2D
ustr2D = 'ULON ULAT time' , & ! vcoord for U cell quantities, 2D
nstr2D = 'NLON NLAT time' , & ! vcoord for N cell quantities, 2D
estr2D = 'ELON ELAT time' , & ! vcoord for E cell quantities, 2D
tstr3Dz = 'TLON TLAT VGRDi time',& ! vcoord for T cell quantities, 3D
ustr3Dz = 'ULON ULAT VGRDi time',& ! vcoord for U cell quantities, 3D
nstr3Dz = 'NLON NLAT VGRDi time',& ! vcoord for N cell quantities, 3D
estr3Dz = 'ELON ELAT VGRDi time',& ! vcoord for E cell quantities, 3D
tstr3Dc = 'TLON TLAT NCAT time',& ! vcoord for T cell quantities, 3D
ustr3Dc = 'ULON ULAT NCAT time',& ! vcoord for U cell quantities, 3D
nstr3Dc = 'NLON NLAT NCAT time',& ! vcoord for N cell quantities, 3D
estr3Dc = 'ELON ELAT NCAT time',& ! vcoord for E cell quantities, 3D
tstr3Db = 'TLON TLAT VGRDb time',& ! vcoord for T cell quantities, 3D
ustr3Db = 'ULON ULAT VGRDb time',& ! vcoord for U cell quantities, 3D
nstr3Db = 'NLON NLAT VGRDb time',& ! vcoord for N cell quantities, 3D
estr3Db = 'ELON ELAT VGRDb time',& ! vcoord for E cell quantities, 3D
tstr3Da = 'TLON TLAT VGRDa time',& ! vcoord for T cell quantities, 3D
ustr3Da = 'ULON ULAT VGRDa time',& ! vcoord for U cell quantities, 3D
nstr3Da = 'NLON NLAT VGRDa time',& ! vcoord for N cell quantities, 3D
estr3Da = 'ELON ELAT VGRDa time',& ! vcoord for E cell quantities, 3D
tstr3Df = 'TLON TLAT NFSD time',& ! vcoord for T cell quantities, 3D
ustr3Df = 'ULON ULAT NFSD time',& ! vcoord for U cell quantities, 3D
nstr3Df = 'NLON NLAT NFSD time',& ! vcoord for N cell quantities, 3D
estr3Df = 'ELON ELAT NFSD time',& ! vcoord for E cell quantities, 3D

!ferret
! T grids
tcstr = 'area: tarea' , & ! vcellmeas for T cell quantities
tstr2D = 'TLON TLAT time' , & ! vcoord for T cell, 2D
tstr3Dc = 'TLON TLAT NCAT time', & ! vcoord for T cell, 3D, ncat
tstr3Da = 'TLON TLAT VGRDa time', & ! vcoord for T cell, 3D, ice-snow-bio
tstr3Db = 'TLON TLAT VGRDb time', & ! vcoord for T cell, 3D, ice-bio
tstr3Df = 'TLON TLAT NFSD time', & ! vcoord for T cell, 3D, fsd
tstr4Di = 'TLON TLAT VGRDi NCAT', & ! vcoord for T cell, 4D, ice
ustr4Di = 'ULON ULAT VGRDi NCAT', & ! vcoord for U cell, 4D, ice
nstr4Di = 'NLON NLAT VGRDi NCAT', & ! vcoord for N cell, 4D, ice
estr4Di = 'ELON ELAT VGRDi NCAT', & ! vcoord for E cell, 4D, ice
tstr4Ds = 'TLON TLAT VGRDs NCAT', & ! vcoord for T cell, 4D, snow
ustr4Ds = 'ULON ULAT VGRDs NCAT', & ! vcoord for U cell, 4D, snow
nstr4Ds = 'NLON NLAT VGRDs NCAT', & ! vcoord for N cell, 4D, snow
estr4Ds = 'ELON ELAT VGRDs NCAT', & ! vcoord for E cell, 4D, snow
tstr4Df = 'TLON TLAT NFSD NCAT', & ! vcoord for T cell, 4D, fsd
ustr4Df = 'ULON ULAT NFSD NCAT', & ! vcoord for U cell, 4D, fsd
nstr4Df = 'NLON NLAT NFSD NCAT', & ! vcoord for N cell, 4D, fsd
estr4Df = 'ELON ELAT NFSD NCAT' ! vcoord for E cell, 4D, fsd
!ferret
! tstr4Di = 'TLON TLAT VGRDi NCAT time', & ! ferret can not handle time
! ustr4Di = 'ULON ULAT VGRDi NCAT time', & ! index on 4D variables.
! tstr4Ds = 'TLON TLAT VGRDs NCAT time', & ! Use 'ferret' lines instead
! ustr4Ds = 'ULON ULAT VGRDs NCAT time', & ! (below also)
! tstr4Db = 'TLON TLAT VGRDb NCAT time', &
! ustr4Db = 'ULON ULAT VGRDb NCAT time', &
! tstr4Df = 'TLON TLAT NFSD NCAT time', &
! ustr4Df = 'ULON ULAT NFSD NCAT time', &
! U grids
ucstr = 'area: uarea' , & ! vcellmeas for U cell quantities
ustr2D = 'ULON ULAT time' , & ! vcoord for U cell, 2D
! N grids
ncstr = 'area: narea' , & ! vcellmeas for N cell quantities
nstr2D = 'NLON NLAT time' , & ! vcoord for N cell, 2D
! E grids
ecstr = 'area: earea' , & ! vcellmeas for E cell quantities
estr2D = 'ELON ELAT time' ! vcoord for E cell, 2D

!---------------------------------------------------------------
! flags: write to output file if true or histfreq value
Expand Down
4 changes: 2 additions & 2 deletions cicecore/cicedynB/dynamics/ice_dyn_eap.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -270,8 +270,8 @@ subroutine eap (dt)
! convert fields from T to U grid
!-----------------------------------------------------------------

call grid_average_X2Y('F', tmass , 'T' , umass, 'U')
call grid_average_X2Y('F', aice_init, 'T' , aiU , 'U')
call grid_average_X2Y('S', tmass , 'T' , umass, 'U')
call grid_average_X2Y('S', aice_init, 'T' , aiU , 'U')
call grid_average_X2Y('S', uocn , grid_ocn_dynu, uocnU , 'U')
call grid_average_X2Y('S', vocn , grid_ocn_dynv, vocnU , 'U')
call grid_average_X2Y('S', ss_tltx , grid_ocn_dynu, ss_tltxU, 'U')
Expand Down
12 changes: 6 additions & 6 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 @@ -317,22 +317,22 @@ subroutine evp (dt)
! convert fields from T to U grid
!-----------------------------------------------------------------

call grid_average_X2Y('F', tmass , 'T' , umass , 'U')
call grid_average_X2Y('F', aice_init, 'T' , aiU , 'U')
call grid_average_X2Y('S', tmass , 'T' , umass , 'U')
call grid_average_X2Y('S', aice_init, 'T' , aiU , 'U')
Copy link
Contributor

@eclare108213 eclare108213 Sep 23, 2022 •
edited
Loading

Choose a reason for hiding this comment

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

The change here for tmass is right: to get the correct mass at the U point, we need to compute
sum(rhoi*vice*aice*T cell area)/sum(aice*T cell area)
over the T cells [ignoring snow for now], since vice has units of volume per unit area.

But I do not think this is correct for aice_init. Since ice concentration is assumed to be uniformly distributed over each grid cell, the concentration at U should be
sum(aice_init*T cell area)/sum(T cell area)
but I think the calculation here is
sum(aice_init*aice*T cell area)/sum(aice*T cell area).

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Just to clarify, The "S" mapping is masked area weighted while the "A" mapping is just area weighted. When you use "aice" above, is that the ice fraction? There are no ice fraction weighted averages built into CICE. If we want to do that, we can use user defined averaging weights. Again,

"S":  sum(var*Tmask*Tarea)/sum(Tmask*Tarea)
"A":  sum(var*Tarea)/sum(Tarea)
"F":  sum(var*Tarea_var)/Tarea_dst

Where Tarea and Tmask are the static grid cell areas and masks computed at initialization. If we do want to have fraction weighted mapping, we can do that, but it's not the standard "S", "A", or "F". To do that we'd do

mask_s(:,:,;) = 1
wght_s(:,:,:) = Tarea(:,:,:) * aice(:,:,:)
call grid_average_X2Y('A', tmass, 'T', wght_s, mask_s, umass, 'U')

and that would do

umass = sum(tmass*Tarea*aice)/sum(Tarea*aice)

At this point, we NEVER weigh the grid averaging by the ice fraction anywhere. Is that something we should be doing?

Copy link
Contributor

Choose a reason for hiding this comment

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

Sorry, I misunderstood "area weighted" in the doc to mean ice-area-weighted, but it clearly says grid cell area. My mistake. Let me think about this some more - what I wrote above is not quite correct.

Copy link
Contributor

Choose a reason for hiding this comment

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

Correction to my comment above:

To get the correct mass at the U point, we need to compute
sum(rhoi*vice*Tarea)/sum(Tarea)
over the T cells [ignoring snow for now], since vice has units of volume per unit area and already incorporates aice (i.e. vice = hi * aice). So the tmass line is correct.

Since ice concentration is assumed to be uniformly distributed over each grid cell, the initial concentration at U should be
sum(aice_init*Tarea)/sum(Tarea)
and that is what's being done in the S mapping.

So I agree with both of these modifications. Sorry for the confusion!

Copy link
Contributor Author

Choose a reason for hiding this comment

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

The "S" mapping includes the mask, the "A" mapping does not. I think we need to change "S" to "A" for umass and aiU if we want an unmasked weighted average as written above. This only matters for gridcells near land. So should the average be an average of only ocean cells ("S" mapping) or an average of all ocean/land cells ("A" mapping)?

One other thing that I just want reiterate. Maybe there is a place for fraction weighted mapping in CICE. We do that in the coupler when we map ice coupling fields to other grids. We fraction area weight the mapping so a gridcell with fraction=0.9 has much greater weight than a gridcell with fraction=0.2 when interpolating from ice to atmosphere for example. You can also think of it as a fraction weighted merge. How we do that depends whether it's a flux (which has to be strictly conserved and has a particular mapping requirement) or a state (which could be done many ways). Maybe we need to do the same in CICE in some places for some fields?? I guess it partly depends whether we can define an ice fraction on the grid we're mapping from. We'd have to think about it a bit.

Copy link
Contributor

Choose a reason for hiding this comment

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

I think S is correct in this case, with the mask included, although it would be good for @JFLemieux73 @dabail10 to confirm that it's consistent with other assumptions in the C-grid approach. On a B-grid, the velocity being computed would be zero if any of the T cells is land (i.e. the point is moot).

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Great, thanks for the clarification.

Copy link
Contributor

Choose a reason for hiding this comment

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

I also wondered whether ice area weighting could/should be included as part of the averaging infrastructure. For coupling, the actual averaging happens in the coupler with other components' data, so we can't control the entire average anyhow. Once the model starts being coupled with different grids, we should check whether we are calculating the coupling fields appropriately, but I think you've already thought through all of that. In the ice model, we "aggregate" across ice thickness categories, but that's only within a particular cell, not for shifting things around horizontally (there are vertical grid remappings). The tracer loading/unloading also involves a lot of multiplying/dividing by "weights" (including ice area and other tracers). Neither of these cases involve grid shifts. So: I don't think we should worry about it unless we come across a particular case in which we need it.

call grid_average_X2Y('S', uocn , grid_ocn_dynu, uocnU , 'U')
call grid_average_X2Y('S', vocn , grid_ocn_dynv, vocnU , 'U')
call grid_average_X2Y('S', ss_tltx , grid_ocn_dynu, ss_tltxU, 'U')
call grid_average_X2Y('S', ss_tlty , grid_ocn_dynv, ss_tltyU, 'U')

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', tmass , 'T' , emass , 'E')
call grid_average_X2Y('S', 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('S', ss_tltx , grid_ocn_dynu, ss_tltxE, 'E')
call grid_average_X2Y('S', ss_tlty , grid_ocn_dynv, ss_tltyE, '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', tmass , 'T' , nmass , 'N')
call grid_average_X2Y('S', 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', ss_tltx , grid_ocn_dynu, ss_tltxN, 'N')
Expand Down
4 changes: 2 additions & 2 deletions cicecore/cicedynB/dynamics/ice_dyn_vp.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -304,8 +304,8 @@ subroutine implicit_solver (dt)
! convert fields from T to U grid
!-----------------------------------------------------------------

call grid_average_X2Y('F',tmass , 'T', umass, 'U')
call grid_average_X2Y('F',aice_init, 'T', aiU , 'U')
call grid_average_X2Y('S',tmass , 'T', umass, 'U')
call grid_average_X2Y('S',aice_init, 'T', aiU , 'U')
call grid_average_X2Y('S',uocn , grid_ocn_dynu, uocnU , 'U')
call grid_average_X2Y('S',vocn , grid_ocn_dynv, vocnU , 'U')
call grid_average_X2Y('S',ss_tltx, grid_ocn_dynu, ss_tltxU, 'U')
Expand Down
6 changes: 3 additions & 3 deletions cicecore/cicedynB/general/ice_forcing.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4129,8 +4129,8 @@ subroutine ocn_data_ncar_init_3D

work1(:,:,:) = ocn_frc_m(:,:,:,n ,m)
work2(:,:,:) = ocn_frc_m(:,:,:,n+1,m)
call grid_average_X2Y('F',work1,'T',ocn_frc_m(:,:,:,n ,m),'U')
call grid_average_X2Y('F',work2,'T',ocn_frc_m(:,:,:,n+1,m),'U')
call grid_average_X2Y('A',work1,'T',ocn_frc_m(:,:,:,n ,m),'U')
call grid_average_X2Y('A',work2,'T',ocn_frc_m(:,:,:,n+1,m),'U')

enddo ! month loop
enddo ! field loop
Expand Down Expand Up @@ -4373,7 +4373,7 @@ subroutine ocn_data_hadgem(dt)
use ice_domain, only: nblocks
use ice_domain_size, only: max_blocks
use ice_flux, only: sst, uocn, vocn
use ice_grid, only: grid_average_X2Y, ANGLET
use ice_grid, only: ANGLET

real (kind=dbl_kind), intent(in) :: &
dt ! time step
Expand Down
4 changes: 4 additions & 0 deletions cicecore/cicedynB/general/ice_state.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -109,6 +109,8 @@ module ice_state
public :: &
uvel , & ! x-component of velocity on U grid (m/s)
vvel , & ! y-component of velocity on U grid (m/s)
uvelT , & ! x-component of velocity on T grid (m/s)
vvelT , & ! y-component of velocity on T grid (m/s)
uvelE , & ! x-component of velocity on E grid (m/s)
vvelE , & ! y-component of velocity on E grid (m/s)
uvelN , & ! x-component of velocity on N grid (m/s)
Expand Down Expand Up @@ -155,6 +157,8 @@ subroutine alloc_state
aice0 (nx_block,ny_block,max_blocks) , & ! concentration of open water
uvel (nx_block,ny_block,max_blocks) , & ! x-component of velocity on U grid (m/s)
vvel (nx_block,ny_block,max_blocks) , & ! y-component of velocity on U grid (m/s)
uvelT (nx_block,ny_block,max_blocks) , & ! x-component of velocity on T grid (m/s)
vvelT (nx_block,ny_block,max_blocks) , & ! y-component of velocity on T grid (m/s)
uvelE (nx_block,ny_block,max_blocks) , & ! x-component of velocity on E grid (m/s)
vvelE (nx_block,ny_block,max_blocks) , & ! y-component of velocity on E grid (m/s)
uvelN (nx_block,ny_block,max_blocks) , & ! x-component of velocity on N grid (m/s)
Expand Down
30 changes: 19 additions & 11 deletions cicecore/cicedynB/general/ice_step_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -209,7 +209,7 @@ subroutine step_therm1 (dt, iblk)
Qa_iso, Qref_iso, fiso_evap, HDO_ocn, H2_16O_ocn, H2_18O_ocn
use ice_grid, only: lmask_n, lmask_s, tmask
use ice_state, only: aice, aicen, aicen_init, vicen_init, &
vice, vicen, vsno, vsnon, trcrn, uvel, vvel, vsnon_init
vice, vicen, vsno, vsnon, trcrn, uvelT, vvelT, vsnon_init
#ifdef CICE_IN_NEMO
use ice_state, only: aice_init
#endif
Expand Down Expand Up @@ -251,8 +251,8 @@ subroutine step_therm1 (dt, iblk)
tr_pond_lvl, tr_pond_topo, calc_Tsfc, highfreq, tr_snow

real (kind=dbl_kind) :: &
uvel_center, & ! cell-centered velocity, x component (m/s)
vvel_center, & ! cell-centered velocity, y component (m/s)
uvelTij, & ! cell-centered velocity, x component (m/s)
vvelTij, & ! cell-centered velocity, y component (m/s)
puny ! a very small number

real (kind=dbl_kind), dimension(n_aero,2,ncat) :: &
Expand Down Expand Up @@ -337,13 +337,11 @@ subroutine step_therm1 (dt, iblk)
do i = ilo, ihi

if (highfreq) then ! include ice velocity in calculation of wind stress
uvel_center = p25*(uvel(i,j ,iblk) + uvel(i-1,j ,iblk) & ! cell-centered velocity
+ uvel(i,j-1,iblk) + uvel(i-1,j-1,iblk)) ! assumes wind components
vvel_center = p25*(vvel(i,j ,iblk) + vvel(i-1,j ,iblk) & ! are also cell-centered
+ vvel(i,j-1,iblk) + vvel(i-1,j-1,iblk))
uvelTij = uvelT(i,j,iblk)
vvelTij = vvelT(i,j,iblk)
else
uvel_center = c0 ! not used
vvel_center = c0
uvelTij = c0
vvelTij = c0
endif ! highfreq

if (tr_snow) then
Expand Down Expand Up @@ -391,8 +389,8 @@ subroutine step_therm1 (dt, iblk)
vicen = vicen (i,j,:,iblk), &
vsno = vsno (i,j, iblk), &
vsnon = vsnon (i,j,:,iblk), &
uvel = uvel_center , &
vvel = vvel_center , &
uvel = uvelTij , &
vvel = vvelTij , &
Tsfc = trcrn (i,j,nt_Tsfc,:,iblk), &
zqsn = trcrn (i,j,nt_qsno:nt_qsno+nslyr-1,:,iblk), &
zqin = trcrn (i,j,nt_qice:nt_qice+nilyr-1,:,iblk), &
Expand Down Expand Up @@ -944,6 +942,8 @@ subroutine step_dyn_horiz (dt)
use ice_dyn_vp, only: implicit_solver
use ice_dyn_shared, only: kdyn
use ice_flux, only: init_history_dyn
use ice_grid, only: grid_average_X2Y
use ice_state, only: uvel, vvel, uvelT, vvelT
use ice_transport_driver, only: advection, transport_upwind, transport_remap

real (kind=dbl_kind), intent(in) :: &
Expand All @@ -961,6 +961,14 @@ subroutine step_dyn_horiz (dt)
if (kdyn == 2) call eap (dt)
if (kdyn == 3) call implicit_solver (dt)

!-----------------------------------------------------------------
! Compute uvelT, vvelT
! only needed for highfreq, but compute anyway
!-----------------------------------------------------------------

call grid_average_X2Y('A', uvel, 'U', uvelT, 'T')
call grid_average_X2Y('A', vvel, 'U', vvelT, 'T')

!-----------------------------------------------------------------
! Horizontal ice transport
!-----------------------------------------------------------------
Expand Down
7 changes: 5 additions & 2 deletions cicecore/cicedynB/infrastructure/ice_restart_driver.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -285,10 +285,10 @@ subroutine restartfile (ice_ic)
stresspU, stressmU, stress12U
use ice_flux, only: coszen
use ice_grid, only: tmask, grid_type, grid_ice, &
iceumask, iceemask, icenmask
iceumask, iceemask, icenmask, grid_average_X2Y
use ice_state, only: trcr_depend, aice, vice, vsno, trcr, &
aice0, aicen, vicen, vsnon, trcrn, aice_init, uvel, vvel, &
uvelE, vvelE, uvelN, vvelN, &
uvelE, vvelE, uvelN, vvelN, uvelT, vvelT, &
trcr_base, nt_strata, n_trcr_strata

character (*), optional :: ice_ic
Expand Down Expand Up @@ -402,6 +402,9 @@ subroutine restartfile (ice_ic)
'vvelN',1,diag,field_loc_Nface, field_type_vector)
endif

call grid_average_X2Y('A', uvel, 'U', uvelT, 'T')
call grid_average_X2Y('A', vvel, 'U', vvelT, 'T')

!-----------------------------------------------------------------
! radiation fields
!-----------------------------------------------------------------
Expand Down