Impact of l_fixed_area flag on C grid gx3/gx1 simulations #809
Comments
|
I use the latest code with commit ID 0bf0fdc I set up the experiments by doing: ./cice.setup -c CgridFixedArea -g gx3 -m ppp6 -e intel -p 8x1 -s diag1,buildincremental As in issue #702 I use these settings: ndte=1200 note also that kstrength=0 with Pstar=27.5e3 |
|
Things are always more complicated than expected...The C grid run (dt=3600) with l_fixed_area=.true. just crashed. I get the following message: (abort_ice) error = (diagnostic_abort)ERROR: negative area (ice) |
|
@JFLemieux73, I'm sorry to hear that. I guess the documentation was right to warn that the scheme hasn't been fully tested. When I got that error, I would write out a lot of diagnostics for the cell in question. Once I understood the geometry of the situation, I could compare it to the code logic and figure out why the logic failed. Have you set bugcheck = .true. in ice_transport_remap.F90? You might get some useful info. For instance, there's some code to check whether the sum of triangle areas across the edge is equal to the prescribed edgearea. If these don't match, it wouldn't be surprising to see a negative area. One possibility is that the small midpoint displacement, which is supposed to give exactly the right area, pushes the midpoint from one side of the edge to the other. If so, we might need some logic to handle that situation differently. At the end of the day, we'll have to decide if setting l_fixed_area = .true. is a useful quick fix, or more of a bag of worms. |
|
Thanks Bill for your help. I will investigate and keep you informed. |
|
Sorry for the delay...I have some time now to investigate. It looks like it is related to an issue at the ice edge. I have set bugcheck to .true. The failure is related to a point that didn't have mass before the advection. I get the following message in the log file: The global i,j are i=49, j=12 (in southern hemisphere). Is it possible that a different condition at the ice edge (e.g. freedrift as opposed to zero velocities) would fix the problem? |
|
The 4 (C grid) velocity components for that cell are: |
|
I don't think freedrift would fix this problem. If there's no ice in the cell being advected from, even if there is ice motion at the edge of it because there is ice in the neighboring cell, then the area advected should be 0. Seems like a roundoff problem. |
|
But I think that with 3 velo components = 0 and only one v = -0.21 leads to strong divergence. I guess with freedrift velocities the divergence would be a lot smaller....this could mitigate the problem. |
|
Mitigate it in this instance, perhaps, but there will be another case where it still happens... |
|
There must be a line in the current code that makes sure that the amount of ice being advected out of a cell (or a triangle within a cell) is no greater than the amount of ice available in it, and l_fixed_area either skips that or breaks it somehow. |
|
Ok I am back working on this issue! The crash happens at istep1=48. It crashes because there is a negative mass (after the transport) at cell 50,13. Note that the initial mass (before the transport at istep1=48) is zero at cell 50,13 (ice edge) for n (cat) = 5. The negative mass is obtained in subroutine update_fields. The new mass depends on the four fluxes mflx: w1 = mflxe(i,j) - mflxe(i-1,j ) & print *, 'mass at 50,13', mm(i,j), mflxe(i,j), mflxe(i-1,j), mflxn(i,j), mflxn(i,j-1), w1 The mass is -2.07e-11 after the transport. This is caused by only one flux which is not zero mflxn(i,j-1)=mflxn(50,12)=-0.88867. |
|
So the negative mass is obtained with l_fixed_area=.true. Here I do the same 47 time steps with l_fixed_area=.true. but at istep1=48 I set l_fixed_area=.false. print *, 'mass at 50,13', mm(i,j), mflxe(i,j), mflxe(i-1,j), mflxn(i,j), mflxn(i,j-1), w1 The mass is in this case positive. The same mflxn(50,12)=-0.88867 is there but in this is case mflxe(i-1,j) is not zero = 283.01. This leads to w1<0 and mass increase in the grid cell (as opposed to a decrease with l_fixed_area=.true.). Why is mflxe(i-1,j)=0 with l_fixed_area=.true.? Is it a question of mask?...in other words mflxe(i-1,j) is not calculated as it is on the ice edge. |
|
mflxe is calc in subroutine transport_integrals. With l_fixed_area=.false. only for istep1=48 I get in transport_integrals: if (istep1 == 48 .and. nn == 5 .and. nedge == 0 .and. i == 49 .and. j == 13) then This is indeed the same flux of 283.01. Now what do I get with l_fixed_area=.true.? It does not print anything...mlfxe(49,13) is not calculated. This is why it is zero. I guess it is not in icells(ng). |
|
I think it is in subroutine locate_triangles that it is determined where the fluxes are calc. icells intent(out) ! number of cells where triarea > puny Is the problem just below this comment? ! Compute mask for edges with nonzero departure areas |
|
It looks like this is the problem. if l_fixed_area=T, the i,j where the mflx (for north and east) are calc depend on the condition edgearea(i,j) /= c0 while if l_fixed_area=F, the condition is based on the departure points being non zero. As edgearea(49,13)=0 (because uvelE(49,13)=0), mlfxe(49,13) is not calculated when l_fixed_area=T. If l_fixed_area=F and edge='east', the condition is if (dpx(i,j-1)/=c0 .or. dpy(i,j-1)/=c0 & mlfxe(49,13) will be calc because dpx(i,j-1)/=c0 (also dpy(i,j-1)/=c0). Note that the dpx and dpy are calc with the corner velocities (uvel, vvel). |
|
Ok there is no more crash with l_fixed_area=T if I use the same approach as with l_fixed_area=F to get the masks for the edges. In other words it works fine if the masks are obtained with if (trim(edge) == 'north') then |
|
Bonjour @whlipscomb, I think I found why it crashes with l_fixed_area=T. Please look at the previous comments for details. To fix the problem I would like to use the same approach for l_fixed_area=T or F to compute the masks for the edges. The masks for l_fixed_area=T would not be calculated anymore with Would it be ok? Do you have any concerns if I do this? Merci! |
|
Bonjour @JFLemieux73, Thanks for the careful and detailed analysis. It looks like you have a good solution, but I'd like to understand the dynamic/geometric configuration a little better. What is the reason that uvelE(49,13)=0? Is cell (49,13) a land cell? When the c-grid velocities are averaged to the vertices of east edge (49,13), does one vertex have a velocity of exactly zero, while the other is nonzero? And if (49,13) is a land cell, what is the source of the mass going across east edge (49,13) into cell (50,13)? Is it coming from a corner triangle located in cell (49,12) or (49,14)? What is the source of the mass that leaves cell (50,13), giving a negative mass? Since this cell initially is massless in category 5, then my guess is that the mass which is exiting the cell has entered across a different edge from a neighbor cell, and the problem is that with the current logic, CICE counts the mass leaving but not the mass entering. Finally, with the extra 'Bentsen' triangle, do the departure regions across east edge (49,13) have a net area of zero, as required? In other words, I'd like to see a picture showing the various departure triangles. It would be easier if we were sitting in a room with a piece of paper. Let me know if you'd like to try a Zoom call or something similar. |
|
Hi Bill, Cell (49,13) is not a land cell. The problem here is at the ice edge. The total concentration values are close to a_min=1e-3 (mask for momentum equation). Because aice at uvelE(49,13) (C-grid u component) is lower than a_min, then uvelE(49,13) is not calculated and is equal to 0. With the current code with l_fixed_area=T and the C grid, the edge areas (east) are given by: edgearea_e(i,j) = uvelE(i,j,iblk) * HTE(i,j,iblk) * dt The mask for the edges would not include (49,13) because uvelE(49,13)=0. This is why it does not calculate a flux in from this edge (mlfxe(49,13)=0). However, vvelN (50,12) is not zero (it is equal to -0.21) because aice > a_min. The edgearea is in this case: edgearea_n(i,j) = vvelN(i,j,iblk)*HTN(i,j,iblk) * dt. It is not equal to zero and this is why the flux is calculated for this edge (mflxn(50,12) is not equal to zero). The other two velocities (uvelE(50,13) and vvelN(50,13) are also zero so there is no flux calculated. So basically you are right: "with the current logic, CICE counts the mass leaving but not the mass entering." In this case the B grid with l_fixed_area = T would not have any problems because uvel,vvel (49,12) are not zero. The flux mlfxe(49,13) would be calculated because edgearea_e(i,j) = (uvel(i,j,iblk) + uvel(i,j-1,iblk)) & |
|
I would add also that for the B grid, given a certain velocity field, the same mask would be computed for the edges whether l_fixed_area is T or F while it is not the case with the C grid. Let's look for the east edge: For the B grid: if (l_fixed_area) then edgearea(i,j) is not equal to zero if uvel(i,j) and/or uvel(i,j-1) are not equal to zero. This is equivalent to if (dpx(i,j-1)/=c0 .or. dpy(i,j-1)/=c0 & because the departure points depend on the velocities at the corners. This is not the case right now with the C grid: the masks for the edges might be different whether l_fixed_area=T or F. |
|
So again I propose the following to solve the problem: the mask for the edges would always be calculated the same way whether l_fixed_area=T or F: icellsd = 0 |
|
This seems like a reasonable approach to me. Does it work? |
|
It seems to work. I tested the new code above in a gx3 simulation for C grid and l_fixed_area=T and it ran smoothly for 1 year. I am doing longer simulations right now as I would like to compare the ice volumes for the B and C grids. |
|
@JFLemieux73 – I think I understand now. To summarize:
If I've understood correctly, then I agree that the new logic makes sense. Thanks again for working this out. I'm curious how the Bentsen triangle logic works in this case. Since the edgearea(49,13) = 0, I would think that there has to be an extra triangle extending into (50,13). This isn't realistic, but I'd expect it doesn't make much difference near the ice edge. |
|
Thanks Bill. Good summary. I used the new code (described above) and ran 4 year gx3 simulations to see the impact on the volume with l_fixed_area = T or F. The C grid with remap and l_fixed_area = T used to crash after 2 days; now it ran without any problem for the 4 year long simulation. Note that the volume time series for the B-grid are very similar whether l_fixed_area is T or F. I therefore only show it for l_fixed_area=F (referred to as B,F in the figures). Here are the volume time series for the Arctic when using remap:
As we had seen before, remap with l_fixed_area=F leads to an overestimation (compared to the B grid) of the maximum. C grid with l_fixed_area=T (C,T) is closer to the B grid solution except in summer where there is less volume. |
|
I also looked at the impact of upwind instead of remap:
|
|
For the Antarctic and remap we get:
and with upwind:
|
|
I think it looks ok. There are bigger differences when changing the advection scheme (remap versus upwind) than when changing the grid (B versus C). |
|
If everyone agrees I will close this issue and the other ones related to the C grid noise. |
|
The code has been modified so that l_fixed_area=.true. when using the C grid. See PR #833 for details. |
|
I need to reopen this issue. As mentioned in PR #833 there was a crash after 10 months in a gx1 simulation. I think the problem is in a corner. i=27,j=40, iblk=57, iglob=26, jglob=375. Here are the 4 fluxes where it crashes: write (nu_diag,*) 'jlem flux =', mflxe(i,j),mflxe(i-1,j),mflxn(i,j),mflxn(i,j-1) |
|
Ok the models restarts at the same time level but I set l_fixed_area=F instead of T. Here are the fluxes at the same cell: Why is the first one so different?...Looks like a masking issue again. |
|
The u,v and departure points are uvel 0.000000000000000E+000 -6.380189507267686E-005 0.000000000000000E+000 0.000000000000000E+000 So only u,vel(27,39) are non-zero...makes sense because the other u,vvel touch land. |
|
As expected, uvel,vvel, dpx and dpy are the same whether l_fixed_area=T or F. |
|
if (iblk == 57) then jlem edge 11950450.3078016 -21976011.1617095 v varies smoothly but u is weird...it goes from 0.0684 to -0.0657 over one dy...this is surprising. |
|
With uvel(27,39,iblk)=-6.38e-5 and vvel(27,39,iblk)=-0.234, the departure point is in cell (28,40). The standard remapping (l_fixed_area=F) should consider a flux coming in on the eastern edge of cell (27,40). I guess with l_fixed_area=T, remapping would correct this flux because uvelE(27,40) > 0 (there is certainly a contribution for a flux going out). |
|
Another way to understand this is that in this case the edgearea_e(27,40) is not consistent with the uvel used for the departure points (i,e, uvel(27,39) and uvel(27,40) = 0). A way to fix this problem is to calculate the edgearea as done for the B grid when there is this inconsistency: edgearea_e(i,j) = uvelE(i,j,iblk) * HTE(i,j,iblk) * dt if ( abs(uvelE(i,j,iblk)) .gt. 0.05d0) then A similar approach can be used for edgearea_n(i,j). I tested it and it works. There is no more crash. This is however a bit intuitive as I dont know all the details of remapping. Moreover, the 0.05 m/s above is a bit arbitrary... I will need to understand this a bit more before pushing this in the main code and probably have a chat with Bill about this. |
|
Possible other solutions: if ( abs(uvelE(i,j,iblk)) .gt. 0.05d0) then
edgearea_e(i,j)=(edgearea_eB + edgearea_e(i,j))/2 OR edgearea_e(i,j) = alpha * edgearea_e(i,j) ! 0 < alpha < 1 endif |
|
Bill's idea: if ( abs(uvelE(i,j,iblk)) .gt. 0.05d0) then
|
|
The fix above should work but it would be better first to understand why the flux at the south edge does not match the flux at the east edge. Roundoff error? It would be good to understand better the different triangles at the south edge. |
|
I have looked at the two triangles responsible for the area flux in and out. Recall that the one on the east edge brings area in while the one on the south edge leads to area out. The other triangles do not contribute to the flux in or out (partial concentration for open water =0 in T cell 27,40) If I calculate the area only using the local coordinates (that vary between -0.5 and 0.5) I find that the area in is larger than the area out. This should not cause a negative “mass”. The problems is related to the areafact used to calculate the “true” area. Even though both triangles are pretty much at the same location (close to the U point 27,39), the area for the east edge triangle is calculated with areafact_c while areafact_r is used for the south edge triangle. Because the grid is non-uniform (and I guess distorted in this region), areafact_c is smaller than areafact_r and leads to smaller “true” area for the east edge triangle. I propose a fix. In fact I would argue that areafact should always be calculated by doing a linear interpolation using areafact_l and areafact_r. We could use the mean xp position for the interpolation. areafact = mx+b = ( areafac_r - areafac_l ) xpm + p5( areafac_r + areafac_l ) |
|
@JFLemieux73, thanks for this detailed analysis. This clearly explains the cause of the negative areas. Your analysis points to a general problem with the use of different values of areafac along different edges for what is very nearly the same triangle. I remember that when I wrote this code, I worried that it might fail on highly curved grids for unusual triangles. I'm kind of surprised it hasn't failed before. I like your proposed solution. Let's call it Solution 1. With this solution, areafact will have a value close to areafac_r from the point of view of both the east and south edges. That's what we want for a small departure triangle that lies close to the SE vertex of the cell. With Solution 1, the two areafact values (one for the east-edge triangle and one for the almost-identical south-edge triangle) will still differ slightly. The difference is partly because the midpoints are slightly different, but also because the triangles are of finite size. With the new linear formula, the areafact values would converge only in the limit that the triangle has zero area, with a midpoint lying at the vertex such that areafact = areafac_r for both edges. So we could think about a Solution 2 that might be more robust but also would be harder to code. With this solution, areafact would be determined by the location of the triangle midpoint relative to the areafac values at the four nearest vertices. In this case, areafac would be a weighted average of the values at the four vertices of cell (28,40). Then areafac would be a bilinear function of x and y, instead of a linear function of x. The areafac values for the two triangles would differ only to the extent that the midpoints aren't exactly co-located. My inclination would be to go with Solution 1 for now. If, in the future, there's an area violation resulting from an even smaller discrepancy in the area factors associated with two triangles, someone could do the extra work of coding Solution 2. But before that happens, we might already have moved on to a remapping algorithm that's better designed for a C-grid. |
|
Thanks a lot Bill. I agree that we should go with solution 1 for the moment. I quickly coded solution 1 yesterday and ran the test case that failed with the standard code. It worked with solution 1. I will code it more cleanly and test it in long simulations. |
|
These problems have been fixed in PR#849. |
C grid simulations with l_fixed_area=.false. (default) exhibit checkerboard noise (issue #792) and thicker ice (compared to B grid) in convergence zones (issue #702).
Setting l_fixed_area=.true. eliminates the checkerboard noise (issue #792).
We want to see what is the impact of l_fixed_area flag on ice thickness in gx3/gx1 simulations.
The text was updated successfully, but these errors were encountered: