landfast ice parameterisation with lateral drag #741

mjlosch · 2023-06-27T10:31:45Z

What changes does this PR introduce?

New feature, lateral drag parameterization for landfast ice following Liu et al 2022, JGR

What is the current behaviour?

What is the new behaviour

Compiled with new flag SEAICE_ALLOW_SIDEDRAG and turned on by setting runtime parameter SEAICEsideDrag \=0, this parameterisation replaces the no-slip boundary condition for sea ice by a tunable lateral drag. If available, sub-grid-scale coastline information can be used (to be specified in tow files u/vCoastLineFile) to determine the lateral drag coefficient. Should work with all solvers, but only tested thoroughly with the default Picard/LSR solver.

Does this PR introduce a breaking change?

No

Other information:

This is a slightly modified version of @yqliu11's work, based on her seaice_lateraldrag_v3 branch.
Some details, like names can be discussed, I guess. Documentation is minimal.
No AD code yet, because the AD dynamics is not stable anyway.
Here's the paper: https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2022JC018413

Suggested addition to `tag-index`

o pkg/seaice: add new lateral drag parameterisation for landfast ice (Liu et al 2022, JGR)

- Liu et al 2021, JGR, A new parameterization of coastal drag to simulate landfast ice in deep marginal seas in the Arctic

by listing the runtime parameters and add references

antnguyen13 · 2023-06-30T09:42:20Z

how exciting!

antnguyen13 · 2024-02-29T10:15:26Z

i think it would help to provide some kind of quick script to produce an easy set of [u,v]coastlinefile following , for example, the easy form factor F1 (fig. 1b of Liu et al 2022). Otherwise, it requires quite a bit of work , even for me, to read and sort out the x,y dir and to understand the explanation of how these two files are read in, and understand what the meaning of "Coast line length "normal" to the X/Y dir" parameters NcoastX/Y inside seaice_init_fixed.F is. I've made some quick test using hfac[S,W] for the drag form factor in u,v dir, and still am a bit confused whether i mult with the correct dxG,dyG (to then be centered and div by dxC,dyC in seaice_init_fixed.F...

this is how i generate the ucoastlineFile field, does this make sense to you?

(matlab syntax):
read in hfS top layer only,
read in dxG
tmp=1-hfS (so that land is 1, ocean is 0)
tmpr=0.*tmp;
tmpx(:,2:end-1) = tmpr(:,1:end-2)+tmpr(:,2:end-1)+tmpr(:,3:end) ; %so i get values of 1 or 2, note sum is in y-dir
tmpx=tmpx.*dxg;
write this to a binary file , to be read in under the name ucoastlinefile and assigned to LcoastX in seaice_init_fixed.F

So that NcoastX is now computed by centering the field above, then div by dxC.

mjlosch · 2024-02-29T11:15:05Z

@antnguyen13 for the simple version you don't need to do anything. It is computed from the "model" coastline in seaice_init_fixed.F. If you want something more complicated, I can provide input files for the 36km arctic config.

mjlosch · 2024-02-29T11:20:53Z

In general, the "coast line length" is somewhat arbitrary. In the end it becomes some kind of roughness parameter and it really doesn't matter what it is based on. For @yqliu11 paper, we found it useful to derive this "roughness" from sub grid scale coastline as the projection onto I and J directions. There is a python script that does this for a particular coastline product. That should (have been) also be part of the publication, but it's not. I can try to find it and post it here.

antnguyen13 · 2024-02-29T21:57:15Z

Thanks Martin, if you could find the python script? I get the idea, as you said, you want the drag to act on a certain number of grid points starting from the coast line, more next to coast, less further away. I was just trying to come up with something very crude using hFac as a guide (given that all we have to do is blur the hFac toward the ocean side, then we'd get the coast points and then play with it.. the i/j are always a bit confusing to try to document.. i found from experience that it might be easier to use terminology "tangential" and "normal" ( i did this for the iceplume, after it was also called x/y, and then the user got it mixed up with XC, YC, which can be a disaster for where XC is not x-dir)

antnguyen13 · 2024-03-05T07:36:34Z

pkg/seaice/SEAICE_PARAMS.h

      COMMON /SEAICE_PARM_C/
     &   AreaFile, HsnowFile, HsaltFile, HeffFile,
-     &   uIceFile, vIceFile
+     &   uIceFile, vIceFile, uCoastLineFile, vCoastLineFile


I'd suggest putting the new fields ([u,v]CoastLIneFile, and others relevant to this PR) inside ifdef block the same way you have elsewhere

I think our convention is to have all parameters (that can be specified in a namelist) always defined (so we can catch if they are set even if the corresponding CPP-flags are not defined). So I'll keep this as it is. Similar for the related comments.

On second thought, I should probably add more code that checks if these flags (in particular uCoastLineFile and vCoastLineFile) are set consistently.

see commit 8db339f

antnguyen13 · 2024-03-05T07:36:59Z

pkg/seaice/SEAICE_PARAMS.h

@@ -508,6 +517,7 @@ C
      _RL SEAICE_drySnowAlb_south, SEAICE_wetSnowAlb_south, HO_south
      _RL SEAICE_cBasalStar, SEAICEbasalDragU0
      _RL SEAICEbasalDragK1, SEAICEbasalDragK2
+      _RL SEAICESideDrag


same comment as above, put inside ifdef block

antnguyen13 · 2024-03-05T07:37:13Z

pkg/seaice/SEAICE_PARAMS.h

@@ -566,6 +576,7 @@ C
     &    SEAICE_drySnowAlb_south, SEAICE_wetSnowAlb_south, HO_south,
     &    SEAICE_cBasalStar, SEAICEbasalDragU0,
     &    SEAICEbasalDragK1, SEAICEbasalDragK2,
+     &    SEAICESideDrag,


antnguyen13 · 2024-03-05T07:38:14Z

pkg/seaice/seaice_readparms.F

@@ -163,6 +163,7 @@ SUBROUTINE SEAICE_READPARMS( myThid )
     & SEAICE_drag_south, SEAICE_waterDrag_south,
     & SEAICE_dryIceAlb_south, SEAICE_wetIceAlb_south,
     & SEAICE_drySnowAlb_south, SEAICE_wetSnowAlb_south, HO_south,
+     & SEAICESideDrag, uCoastLineFile, vCoastLineFile,


suggest putting inside ifdef block

antnguyen13 · 2024-03-05T07:38:22Z

pkg/seaice/seaice_readparms.F

@@ -394,6 +395,7 @@ SUBROUTINE SEAICE_READPARMS( myThid )
      SEAICE_drag        = 0.001  _d 0
      OCEAN_drag         = 0.001  _d 0
      SEAICE_waterDrag   = 0.0055 _d 0
+      SEAICESideDrag     = 0.0    _d 0


antnguyen13 · 2024-03-05T07:38:29Z

pkg/seaice/seaice_readparms.F

@@ -509,6 +511,8 @@ SUBROUTINE SEAICE_READPARMS( myThid )
      HeffFile   = ' '
      uIceFile   = ' '
      vIceFile   = ' '
+      uCoastLineFile = ' '


mjlosch · 2024-03-12T10:52:34Z

@yqliu11 do you think you can share the python script that produces the u/vCoastLineFile?

yqliu11 · 2024-03-12T13:35:34Z

compute_coastline_for_grid.zip
here is the jupyer notebook file to calculate the input grid from coastline.

antnguyen13 · 2024-04-03T17:05:10Z

pkg/seaice/SEAICE.h

@@ -188,6 +188,19 @@ C     CbobC :: (linear) bottom drag coefficient for basals stress param.
      _RL CbotC      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
 # endif /* SEAICE_ALLOW_BOTTOMDRAG */

+# ifdef SEAICE_ALLOW_SIDEDRAG
+C     NcoastX/Y   :: coast line roughness (w/out units) in X and Y direction,


marking suggested swapping of X/Y to U/V (which i'll also put in the comment)

antnguyen13 · 2024-04-03T17:05:20Z

pkg/seaice/SEAICE.h

+C                    computed from the coastline length at corner points,
+C                    interpolated to U/V points, and scaled by the grid cell
+C                    width
+C     SideDragX/Y :: drag coefficients for lateral drag a parameterisation


antnguyen13 · 2024-04-03T17:07:36Z

pkg/seaice/SEAICE_PARAMS.h

@@ -324,6 +324,10 @@ C     HsaltFile         :: File containing initial sea ice salt content
 C     HeffFile          :: File containing initial sea-ice thickness
 C     uIceFile          :: File containing initial sea-ice U comp. velocity
 C     vIceFile          :: File containing initial sea-ice V comp. velocity
+C     u/vCoastLineFile  :: Files containing the length of the coastline (in m)
+C                          normal to the u/v direction per grid cell
+C                          at corner points (XG,YG) as a measure of


be careful , remove language tying to [X,Y]G here. also the fact that u/v CoastLineFile is "normal" to u/v direction, that's also confusing...

Here I find it very important that we know where the values contained in u/vCoastLineFile are defined, so we need the reference to correct coordinates which are (for corner points) (XG,YG). I see no way of changing that. One can ask, why they are defined at corner points and not at, say, centers or at U/V points. The reason is this: the tangential points for the slip conditions is a the corner points (e.g. for u you have uIce(i,j)+uIce(i,j-1) = 0 for no slip on the corner point). In the model we then average to U/V-points. I agree that this can be confusing.

The files contain a measure of roughness. This can be anything that you may want to specify. In the specific case of @yqliu11's paper, we assumed that the coastline normal to the velocity creates an obstacle (no flow) and computed the the roughness from the projection of the (sub-grid scale) coastline on the plane normal to the velocity (eq 9 and 10 in https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2022JC018413). I am not sure how to express this in the short description here. @antnguyen13 Do you have a better suggestion? We can just say "roughness in m" or "roughness length in m" (although that's maybe not so good, because "roughness length" has defined meaning in turbulent boundary layer dynamics).

antnguyen13 · 2024-04-03T18:29:26Z

pkg/seaice/seaice_init_fixed.F

+      LOGICAL readCoastLineFields
+C     Coast line length normal to the X and Y direction in meters
+C     read from file and used to compute rough parameters NcoastX/Y
+      _RL LcoastX(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)


antnguyen13 · 2024-04-03T18:29:53Z

pkg/seaice/seaice_init_fixed.F

+      ENDIF
+
+      CALL EXCH_UV_XY_RL( NcoastX, NcoastY, .FALSE., myThid )
+      CALL WRITE_FLD_XY_RL( 'Nlized_coastX', ' ', NcoastX,-1, myThid )


antnguyen13 · 2024-04-03T18:30:30Z

pkg/seaice/seaice_sidedrag_stress.F

+      _RL uIceloc     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
+      _RL vIceloc     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
+      _RL AREALoc     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
+      _RL SideDragXloc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)


antnguyen13 · 2024-04-03T18:30:42Z

pkg/seaice/seaice_sidedrag_stress.F

+      _RL AREALoc     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
+      _RL SideDragXloc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
+      _RL SideDragYloc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
+      _RL NcoastXloc  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)


same as above

antnguyen13 · 2024-04-03T18:30:58Z

pkg/seaice/seaice_sidedrag_stress.F

+         DO j=1-OLy+1,sNy+OLy-1
+          DO i=1-OLx+1,sNx+OLx-1
+           IF ( AREALoc(i,j,bi,bj) .GT. 0.01 _d 0 ) THEN
+            tmpx = 0.5 _d 0 * ( uSpd(i,j) + uSpd(i-1,j) )


mjlosch · 2024-04-04T09:30:44Z

In principle, I like the suggestion to rename the internal variables from X/Y to U/V.
There are many variables in pkg/seaice that are X/Y where they could be U/V, e.g. FORCEX/Y. Do you think that this should be changed, too?

jm-c · 2024-04-04T13:31:07Z

@mjlosch I will take a look in the coming days (for a second thought on variable names). Could you can hold on changes until then ?

mjlosch added 7 commits June 27, 2023 10:09

cleaned up version of Yuqing Liu's lateral drag parameterisation

2f601b5

- Liu et al 2021, JGR, A new parameterization of coastal drag to simulate landfast ice in deep marginal seas in the Arctic

update SEAICE_OPTIONS.h and compile SEAICE_ALLOW_SIDEDRAG code

72c0f7d

add new file required for new parameterization

6954c26

add diagnostics of lateral drag terms

c64ed62

add minimal documenation for basal and lateral drag

4973f0c

by listing the runtime parameters and add references

make lateral drag work with EVP, JFNK, and Krylov solvers

bc30cc1

fix a typo

3135970

jrscott requested review from antnguyen13 and jrscott July 11, 2023 14:58

mjlosch and others added 3 commits October 25, 2023 09:45

Merge branch upstream/master into seaice_sidedrag

871ed4f

Merge branch 'master' into seaice_sidedrag

2b6d25a

add minimal documentation for some new fields and parameters

e85160d

antnguyen13 reviewed Mar 5, 2024

View reviewed changes

check for more erroneously set parameters

8db339f

antnguyen13 reviewed Apr 3, 2024

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

landfast ice parameterisation with lateral drag #741

landfast ice parameterisation with lateral drag #741

mjlosch commented Jun 27, 2023 •

edited

antnguyen13 commented Jun 30, 2023

antnguyen13 commented Feb 29, 2024 •

edited

mjlosch commented Feb 29, 2024

mjlosch commented Feb 29, 2024

antnguyen13 commented Feb 29, 2024

antnguyen13 Mar 5, 2024

mjlosch Mar 12, 2024

mjlosch Mar 13, 2024

mjlosch Mar 14, 2024

antnguyen13 Mar 5, 2024

antnguyen13 Mar 5, 2024

antnguyen13 Mar 5, 2024

antnguyen13 Mar 5, 2024

antnguyen13 Mar 5, 2024

mjlosch commented Mar 12, 2024

yqliu11 commented Mar 12, 2024

antnguyen13 Apr 3, 2024

antnguyen13 Apr 3, 2024

antnguyen13 Apr 3, 2024

mjlosch Apr 4, 2024 •

edited

antnguyen13 Apr 3, 2024

antnguyen13 Apr 3, 2024

antnguyen13 Apr 3, 2024

antnguyen13 Apr 3, 2024

antnguyen13 Apr 3, 2024

mjlosch commented Apr 4, 2024

jm-c commented Apr 4, 2024

landfast ice parameterisation with lateral drag #741

Are you sure you want to change the base?

landfast ice parameterisation with lateral drag #741

Conversation

mjlosch commented Jun 27, 2023 • edited

What changes does this PR introduce?

What is the current behaviour?

What is the new behaviour

Does this PR introduce a breaking change?

Other information:

Suggested addition to tag-index

antnguyen13 commented Jun 30, 2023

antnguyen13 commented Feb 29, 2024 • edited

mjlosch commented Feb 29, 2024

mjlosch commented Feb 29, 2024

antnguyen13 commented Feb 29, 2024

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mjlosch commented Mar 12, 2024

yqliu11 commented Mar 12, 2024

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mjlosch Apr 4, 2024 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mjlosch commented Apr 4, 2024

jm-c commented Apr 4, 2024

mjlosch commented Jun 27, 2023 •

edited

Suggested addition to `tag-index`

antnguyen13 commented Feb 29, 2024 •

edited

mjlosch Apr 4, 2024 •

edited