CICE-Consortium / CICE

Development repository for the CICE sea-ice model
Other
57 stars 131 forks source link

Unnecessary calculations for uvel, vvel for the C-grid #952

Closed JFLemieux73 closed 4 months ago

JFLemieux73 commented 4 months ago

In ice_dyn_evp:

call grid_average_X2Y('S', uvelE, 'E', uvel, 'U') call grid_average_X2Y('S', vvelN, 'N', vvel, 'U')

uvel(:,:,:) = uvel(:,:,:)uvm(:,:,:) vvel(:,:,:) = vvel(:,:,:)uvm(:,:,:)

Let's look at uvel = uU (same idea for vvel). In the grid_average we have:

uvel(i,j) = (uvelE(i,j) AreaE(i,j) MaskE(i,j) + uvelE(i,j+1) AreaE(i,j+1) MaskE(i,j+1)) / (AreaE(i,j) MaskE(i,j) + AreaE(i,j+1) MaskE(i,j+1))

The multiplications by the mask are not necessary. For example, if MaskE(i,j) = 0, the mask at the U point is also zero. and the next operation uvel(:,:,:) = uvel(:,:,:)*uvm(:,:,:) sets uvel to zero.

we could just have:

uvel(i,j) = (uvelE(i,j) AreaE(i,j) + uvelE(i,j+1) AreaE(i,j+1)) / (AreaE(i,j) + AreaE(i,j+1)) uvel(:,:,:) = uvel(:,:,:)*uvm(:,:,:)

apcraig commented 4 months ago

@JFLemieux73, I'd like to create a PR and get this merged ASAP. I'm happy to make the changes, do the testing, create the PR. I'm trying to figure out what needs to happen.

So, if I understand correctly, we could change the

call grid_average_X2Y('S', uvelE, 'E', uvel, 'U')

to

call grid_average_X2Y('A', uvelE, 'E', uvel, 'U')

which is an unmasked average. This is already done that way for the C grid. I see 3 places (for uvel and vvel) where this might be changed in ice_dyn_evp.F90, all in subroutine evp.

line 721+722 line 1087+1088 line 1281+1282

Is that what is needed? This should be bit-for-bit, right? Anything else?

JFLemieux73 commented 4 months ago

Hi @apcraig,

Yes you understand correctly. It should be BFB.

For the other issue I just created (#951) I think we could do it later for another release.

JFLemieux73 commented 4 months ago

FYI @eclare108213 I just verified with version 6.5 (for the paper) that changing the avg type from S to A leads to BFB results (in gx1 simulations).

apcraig commented 4 months ago

Lets do #951 too. I got the first bullet, remove the sentence in 2.1. Can you tell me what you want done in section 2.5. Or feel free to modify the file (doc/source/science_guide/sg_dynamics.rst) and then email it to me or something and I can check it with readthedocs and make it part of the PR. Thanks.

JFLemieux73 commented 4 months ago

I would like @eclare108213 to comment on this. The eqs in the doc in section 2.5 are correct for the northern hemisphere (with fm>0) for the water stress term. Hence, in the doc, waterx = uwcos(theta) - vwsin(theta). In the code, however, waterx = uwcos(theta) - vwsin(theta)*sign(c1,fm). We could do:

1) we add sign(c1,fm) in the doc for waterx, watery and ccb. But I find these eqs are already quite busy. 2) we mention that the formulation of the water stress is for the northern hemisphere and that the sin(theta) terms are multiplied by -1 in the code when the ice is in the southern hemisphere.