oceanmodeling / WW3

WAVEWATCH III
Other
0 stars 0 forks source link

Improve WW3 mesh cap to allow using vortex formulation in ocean side #11

Open janahaddad opened 4 months ago

janahaddad commented 4 months ago

Background

[based on @uturuncoglu email] UFS Coastal is working on coupling WW3 and SCHISM. We're able to run 2D configurations of SCHISM + WW3 with mesh cap (or NCAR cap) without any issues, but we are preparing the way for running in 3D mode (https://github.com/oceanmodeling/ufs-coastal/issues/32). This requires additional fields from WW3.

The aim is to add missing fields to the WW3 cap (I think they are already available through OASIS complier) and complete the wiring with the SCHSIM (the interaction with wave with new formulation is already implemented) to enable 3d coupling

Info needed to proceed

@deniseworthen we will follow up for a description of what we're searching for here !

cc @josephzhang8 @uturuncoglu @yunfangsun @saeed-moghimi-noaa @anntsay

DeniseWorthen commented 4 months ago

The cpl_scalars are used when CMEPs write history or other files and allows the fields on the mesh to be written as 2-D fields. This simply removes the cpl_scalars from the export fields list, since in the unstructured case it doesn't make sense.

janahaddad commented 4 months ago

from @DeniseWorthen email:

In wav_import_export there are several "calc" routines which we moved over from the old esmf cap to the new mesh cap. I think some of the fields I see listed (like the stokes-drift spectrum, calcstokes3d) are there in the old cap, but were not added to the new cap. I think it is straightforward to add them to the new cap.

Some of the other fields---well, I have no idea what internal WW3 fields correspond to the needed export. But, if it is calculated as a "history" field in WW3, I think you could make them available by adding additional "calc" routines in wav_import_export. For example, in wav_grdout, I see

varatts( "BHD ", "BHD ", "Bernoulli head (J term) ", "m2 s-2 ", " ", .false.)

and if I look in w3iogomd, I see how that term is calculated on LN 1653. I think you can find the bottom momentum fluxes also.

As for units, I'm not 100% sure I got them right when I wrote the wav_grdout. Ana has found a couple of errors (I still have on my list to fix these!). I think you'll need to rely on wave modelers to be entirely sure of the units.

janahaddad commented 4 months ago

@uturuncoglu @DeniseWorthen transferring conversation on this issue from email to GH! thx

DeniseWorthen commented 4 months ago

For point of clarify, you do want any export fields to be calculated in the cap, and not just use them from the internal fields that WW3 might calculate as an output variable. The reason is that the output variables are only calculated at some set frequency, and you want the export fields to be updated at each modelAdvance. See https://github.com/NOAA-EMC/WW3/issues/843

uturuncoglu commented 4 months ago

@DeniseWorthen Yes, it would be nice to have them in the cap and updated with coupling frequency. It would be also nice to have them in the output but we could still output import and export states for debugging purpose.

janahaddad commented 4 months ago

from @DeniseWorthen: w3ounfmetamd.F90 is script with BHD, etc. units as above...

saeed-moghimi-noaa commented 4 months ago

Denise Worthen - NOAA Affiliate 1:51 PM All the mesh-cap files are named as wav_*F90

saeed-moghimi-noaa commented 4 months ago

Denise Worthen - NOAA Affiliate 1:53 PM wav_import_export.F90

saeed-moghimi-noaa commented 4 months ago

where to add field to w3: image

saeed-moghimi-noaa commented 4 months ago

image

saeed-moghimi-noaa commented 4 months ago

See in this file (in ufs-weather there is one) in application level: something like this: Denise Worthen - NOAA Affiliate 1:56 PM fd_ufs.yaml

saeed-moghimi-noaa commented 4 months ago

This file is not reserved for ufs-weather (coastal):

Panagiotis Velissariou - NOAA Affiliate 2:00 PM ./tests/parm/fd_ufs.yaml

This one is the one we suppose to use the one under: tests/parm/fd_ufs.yaml

image

saeed-moghimi-noaa commented 4 months ago
saeed-moghimi-noaa commented 4 months ago

Just to document chats from the meeting:

Jana Haddad - NOAA Affiliate 1:33 PM https://github.com/oceanmodeling/WW3/issues/11 Denise Worthen - NOAA Affiliate 1:44 PM w3ounfmetamd.F90 Ali Abdolali 1:46 PM Are we doing 2D or 3D coupling or both? And could you pass the list of variabls you need? for 2D, we need SXX, XYY and SXY, but for 3D, we need more Ali Salimi - NOAA Affiliate 1:47 PM https://docs.google.com/document/d/140_xAMOb30iljGL09C_S7tNAKGDY2j8GEjPb7Pg3Xxg/edit Jana Haddad - NOAA Affiliate 1:47 PM https://docs.google.com/document/d/140_xAMOb30iljGL09C_S7tNAKGDY2j8GEjPb7Pg3Xxg/edit Ali Salimi - NOAA Affiliate 1:47 PM Ali I think they are here Ali Abdolali 1:47 PM thanks Ali Abdolali 1:48 PM SO, as Denise said, for ufs-weather-model, we already have a few variables for coupling. Do you want Ufuk to add the rest to the same location? I would recommend do not go back and forth between wmesmf, OASIS, cstorm, ... because as Denise said, they do not show up magically. Ali Abdolali 1:51 PM the same should apply to import fields to WW3 mesh cap Denise Worthen - NOAA Affiliate 1:51 PM All the mesh-cap files are named as wav_*F90 Denise Worthen - NOAA Affiliate 1:53 PM wav_import_export.F90 Denise Worthen - NOAA Affiliate 1:56 PM fd_ufs.yaml Panagiotis Velissariou - NOAA Affiliate 2:00 PM ./tests/parm/fd_ufs.yaml Denise Worthen - NOAA Affiliate 2:00 PM can I ask---your coastal oceanmodel is now a fork of ufs-wm? Panagiotis Velissariou - NOAA Affiliate 2:00 PM ./CMEPS-interface/CMEPS/mediator/fd_cesm.yaml Denise Worthen - NOAA Affiliate 2:02 PM tests/parm/fd_ufs.yaml Ali Abdolali 2:06 PM where is the caluclation routine? in wave_import? Ali Abdolali 2:16 PM FYI, Denise does not consider herself a WW3 developer :) Ali Salimi - NOAA Affiliate 2:19 PM We are not seeing Denise Worthen - NOAA Affiliate 2:22 PM w3outg is the SR where the output variables are calculated

josephzhang8 commented 4 months ago

Thank you all! Here is the list of arrays we need from WW3:

    wave_hs(1:np)=WW3__OHS !Sig wave height [m]
    wave_dir(1:np)=WW3__DIR !mean wave dir [deg]
    wave_tm1(1:np)=WW3_T0M1 !mean wave period [s]
    wave_wnm(1:np)=WW3__WNM !mean wave number [1/m]
    wave_pres(1:np)=WW3__BHD !wave-induced Bernoulli head pressure [N/m or Pa?]
    wave_stokes_x(1:np)=WW3_USSX !Stokes drift [m/s]
    wave_stokes_y(1:np)=WW3_USSY
    wave_ocean_flux_x(1:np)=WW3_TWOX !wave-ocean mom flux [m2/s2]
    wave_ocean_flux_y(1:np)=WW3_TWOY
    wave_flux_friction_x(1:np)=WW3_TBBX !Momentum flux due to bottom friction [m2/s2]
    wave_flux_friction_y(1:np)=WW3_TBBY
    wave_orbu(1:np)=WW3_UBRX !near bed orbital vel [m/s]
    wave_orbv(1:np)=WW3_UBRY
DeniseWorthen commented 4 months ago

Take the bottom momentum fields as an example: wave_flux_friction_x,y, the "momentum flux due to bottom friction" in m2 s-2. I can see this field is exported via the oasis coupler (in w3gocmmd.F90).

It appears to be calculated in w3sbt4md.F90, at LN 556-570:

! 5. Fills output arrays and estimates the energy and momentum loss
!
DO IK=1, NK
  CONST2=DDEN(IK)/CG(IK) &         !Jacobian to get energy in band
       *GRAV/(SIG(IK)/WN(IK))    ! coefficient to get momentum
  DSUM(IK)=-FW*UORB*CBETA(IK)      !*WSUB(ISUB)
  DO ITH=1,NTH
    IS=ITH+(IK-1)*NTH
    D(IS)=DSUM(IK)
    TEMP2=CONST2*D(IS)*A(IS)
    TAUBBL(1) = TAUBBL(1) - TEMP2*ECOS(IS)
    TAUBBL(2) = TAUBBL(2) - TEMP2*ESIN(IS)
    S(IS)=D(IS)*A(IS)
  END DO
END DO

This is the code (and any preceding code on which it depends) which would need to be added in wav_import_export.F90 as a new "calc" routine. So, the steps would be (I've invented a field name here):

1) add the fields to the list of exported fields

call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_taubot_x')
call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_taubot_y')

2) add in export_fields (a SR in wav_import_export) code similar to the other export variables

a) define pointers to the esmf fields:

    real(r8), pointer :: tauxbot(:)
    real(r8), pointer :: tauybot(:)

b) check that the field is required in the export state using and get pointers to the fields

   if ( state_fldchk(exportState, 'Sw_taubot_x') .and. &
         state_fldchk(exportState, 'Sw_taubot_y') )then
      if (ChkErr(rc,__LINE__,u_FILE_u)) return
      call state_getfldptr(exportState, 'Sw_taubot_x', tauxbot, rc=rc)
      if (ChkErr(rc,__LINE__,u_FILE_u)) return
      call state_getfldptr(exportState, 'Sw_taubot_y', tauybot, rc=rc)
      if (ChkErr(rc,__LINE__,u_FILE_u)) return

c) Create a "calc" routine to calculate the fields, replicating whatever is required from the w3sbt4md.F90 routine.

call CalcBottomStress(input argument list, tauxbot, tauybot)

At this point, the bottom stress should be ready to be exported by WW3 at every model advance. There is more to be done however to get them sent to another component:

1) the new field names need to be defined in the field-dictionary used by all components (in UFS-WM this is in tests/parm/fd_ufs.yaml) 2) the fields need to be added to the FieldExchange module in CMEP. 3) the fields need to be requested and realized by whatever model is expecting these fields as imports.

janahaddad commented 4 months ago

Summary from meeting 5/28:

Next steps:

need to eventually merge import_export work that's already done for the 2d coupling:

in the meantime:

  1. find where Joseph's list of fields above are calc'd as output fields
  2. note the variable name being used & find that in subroutine that calc's the output: w3outg is the SR where the output variables are calculated
  3. copy how it's being calc'd in that output routine and paste it into wave_import_export.f90. Then go into the cap and expose/export the variables to WW3

attendees: @DeniseWorthen @saeed-moghimi-noaa @josephzhang8 @pvelissariou1 @aliabdolali @AliS-Noaa Saeideh Banihashemi @yunfangsun

janahaddad commented 4 months ago

Here is the history showing Ufuk's addition of radiation stress components to wav_import_export.f90, on the UFS Coastal side:

He added to advertise_fields SR:

      call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_wavsuu')
      call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_wavsuv')
      call fldlist_add(fldsFrWav_num, fldsFrWav, 'Sw_wavsvv')

and corrected those names in export_fields SR.

Looks like those changes have not been merged to UFS-WM's WW3 branch. See here.

DeniseWorthen commented 4 months ago

I need to amend what I wrote previously, because it appears that taubbl might be calculated prognostically at each model advance in the w3srcemd routine. In that case, I believe it would not be required to add a new calc routine, but you could simply 'use' the taubbl values as calculated by the model (they're in w3adatmd). Then in export_fields, you'd have instead of calling a new calc routine in step 2c above,

       do jsea=1, nseal_cpl
        call init_get_isea(isea, jsea)
        ix  = mapsf(isea,1)
        iy  = mapsf(isea,2)
        if (mapsta(iy,ix) == 1) then
           tauxbot(jsea) = taubbl(1)
           tauybot(jsea) = taubbl(2)
        end if
      end do
josephzhang8 commented 4 months ago

Thx @DeniseWorthen for your instructions! The calculation of tau*bot you showed above seems to suggest it's from the multi-grid (not UG) version of WW3?

DeniseWorthen commented 4 months ago

@josephzhang8 Thanks, I should have looked more carefully. I think it should actually just be

        do jsea=1, nseal_cpl
           tauxbot(jsea) = taubbl(jsea,1)
           tauybot(jsea) = taubbl(jsea,2)
        end do

since taubbl allocated as (nsealm,2)

josephzhang8 commented 4 months ago

That makes sense. Thx @DeniseWorthen

josephzhang8 commented 2 months ago

The list of WW3 variables SCHISM needs (from OASIS ): wave_hs(1:np)=WW3OHS !Sig wave height [m] wave_dir(1:np)=WW3__DIR !mean wave dir [deg] wave_tm1(1:np)=WW3_T0M1 !mean wave period [s] wave_wnm(1:np)=WW3WNM !mean wave number [1/m] wave_pres(1:np)=WW3__BHD !wave-induced Bernoulli head pressure [N/m or Pa?] wave_stokes_x(1:np)=WW3_USSX !Stokes drift [m/s] wave_stokes_y(1:np)=WW3_USSY wave_ocean_flux_x(1:np)=WW3_TWOX !wave-ocean mom flux [m2/s2] wave_ocean_flux_y(1:np)=WW3_TWOY wave_flux_friction_x(1:np)=WW3_TBBX !Momentum flux due to bottom friction [m2/s2] wave_flux_friction_y(1:np)=WW3_TBBY wave_orbu(1:np)=WW3_UBRX !near bed orbital vel [m/s] wave_orbv(1:np)=WW3_UBRY

josephzhang8 commented 2 months ago

@aliabdolali many thanks for taking time help us set up the test!

You can find the setup for SCHISM-WWM for the analytical case at: https://columbia.vims.edu/schism/schism_verification_tests/Test_WWM_Analytical/

The web may warn you about potential risk; just accept and you'll see the dir (our IT guy is still looking to upgrade). Let me know if you have any questions. Thx!

uturuncoglu commented 2 months ago

@DeniseWorthen @josephzhang8 I think we could use this issue for WW3 implementation. I'll open another one under SCHSIM-ESMF to track issues specifically for ocean component.

DeniseWorthen commented 2 months ago

@uturuncoglu I used the coastal_ike_shinnecock_atm2sch2ww3_gnu to create a sandbox run directory. I created some feature branches for my work off of your hotfix branch. Then I tried to compile the app in debug mode

./compile.sh hercules '-DAPP=CSTLSW -DUSE_ATMOS=ON -DUSE_WW3=ON -DNO_PARMETIS=OFF -DOLDIO=ON -DPDLIB=ON -DBUILD_UTILS=ON -DDEBUG=ON' xflds.gnu.db gnu NO NO 2>&1 | tee xflds.gnu.db.compile.log

I also compiled w/o DEBUG.

Using the same sandbox run directory, with the debug compile I get an message

54:  B_JGS_BLOCK_GAUSS_SEIDEL is used but the Jacobi solver is not choosen 
Please set JGS_USE_JACOBI .eqv. .true.

But the non-debug compile runs fine. This seems very odd to me since I'm using the same exact same mod_def.ww3 file, which is where this namelist setting would be encoded.

Then I tried intel debug, and that gave me a segfault at this line in SCHISM

 8:  3 0x0000000003b29c9d schism_step_()  /work/noaa/nems/dworthen/ufs-coastal/SCHISM-interface/SCHISM/src/Hydro/schism_step.F90:1924

I always try to start development in debug mode. Have you ever tried either GNU or Intel debug?

uturuncoglu commented 2 months ago

@DeniseWorthen We have GNU test under rt_coastal.conf. but not debug. There could be issue with compiling -DDEBUG. Let me check.

uturuncoglu commented 2 months ago

@DeniseWorthen Is this error coming from WW3 tor SCHSIM. Do you have any idea? @josephzhang8 I wonder if you testing SCHISM with debug flags?

DeniseWorthen commented 2 months ago

@uturuncoglu The error about the Jacobi setting is from WW3. It looks like you're running the implicit time-stepping. But both debug and non-debug should be getting this information from the mod_def.ww3. That is what is so odd. In UFS, we definitely test in GNU debug + PDLIB, but not with implicit.

The intel debug error is definitely SCHISM.

uturuncoglu commented 2 months ago

@DeniseWorthen I have no information about those test cases. I just ported them from old version of UFS coastal and make it available through the RTs. You are right it is using implicit and if I remember there was a reason for it but I could remember. Maybe there is a bug in WW3 side that prevents it to compile with debug option when implicit is used. I could try to run DATM+SCHSIM case to see I could reproduce the issue. If I could, I'll try to solve in SCHSIM side.

aliabdolali commented 2 months ago

could you share your ww3_grid.nml and namelist.nml with me, the type of solver is defined there, and when you execute ww3_grid, it should generate mod_def.ww3 with such selection. EMC did not used implicit Jacobi solver due to their b4b limit, but for coastal application, we should use it.

DeniseWorthen commented 2 months ago

I think this is where the mod_def is being created. The one there is identical to what I have in my sandbox.

/work2/noaa/nems/tufuk/RT/NEMSfv3gfs/input-data-20240501/WW3/Ike/atm2sch2wav
josephzhang8 commented 2 months ago

@DeniseWorthen : _step error seems to be triggered by following code:

if defined USE_WWM || defined USE_WW3

    if(RADFLAG.eq.'VOR') then
      tmpx = 0.d0; tmpy = 0.d0;
      do k=kbs(i),nvrt-1
        tmpx = tmpx + (zs(k+1,i)-zs(k,i))*(stokes_hvel_side(1,k,i)+stokes_hvel_side(1,k+1,i))/2.d0
        tmpy = tmpy + (zs(k+1,i)-zs(k,i))*(stokes_hvel_side(2,k,i)+stokes_hvel_side(2,k+1,i))/2.d0
      enddo !k

! n1 = isidenode(1,i); n2 = isidenode(2,i) htot = (eta2(n1)+eta2(n2))/2 + dps(i) uth(:,i) = uth(:,i)-tmpx/max(0.01d0,htot) vth(:,i) = vth(:,i)-tmpy/max(0.01d0,htot) endif !RADFLAG

endif

I init'ed RADFLAG to be 'VOR' (new 3D coupling which is not ready yet) so that's the reason. I'll set it to 'LON' now and let you know when the code is ready. Sorry for the confusion.

josephzhang8 commented 2 months ago

I've pushed a new master; can you plz pull and test. Thanks.

DeniseWorthen commented 2 months ago

The intel debug job is now running, thanks. It took me a bit to find the right repo.

Still really curious why the GNU debug seems not to be recognizing settings in the mod_def.ww3 file. But I can stick w/ intel debug for now.

DeniseWorthen commented 2 months ago

@uturuncoglu Did you check that the radstr fields in the export state are the same as in the internal WW3 history file? In my test they do not match.

uturuncoglu commented 2 months ago

@DeniseWorthen No I did not check them specifically. Do you have any plot that shows the difference? Are they round-off type difference. Does spatial pattern match? Does WW3 doing some kind of time averaging or unit conversion when it is writing to history? If you remember you fixed something in the WW3 output in my PR. Probably, that does not available in your WW3 version. Maybe it is but not sure.

josephzhang8 commented 2 months ago

@DeniseWorthen @uturuncoglu The choices of 2D/3D coupling is controlled by the flag RADFLAG. Since the cap will call some routines inside SCHISM, we need to somehow let the cap know this flag. I assume we can add it to the cap? The flag can be accessed via use schism_glbl. Let me know if this is doable. Thx.

DeniseWorthen commented 2 months ago

@uturuncoglu I put in the fix for the history write, along w/ another one to fix the time axis in the history files.

The differences are not round-off, one set of fields is ~1/2 the magnitude of the others. I think the mediator fields were bigger than the history fields. Let me dig in some more. I know I have plotting scripts for the unstructured I can dig up.

@josephzhang8 At this stage, I don't need to see or fix what SCHISM is doing w/ the radstr fields. I'm just trying to match the wave export fields with what is being calculated inside of WW3 itself. But in general, we can add a SCHISM config attribute to set optional use cases for SCHISM.

josephzhang8 commented 2 months ago

Thx @DeniseWorthen. I think we'll need the flag inside the cap, b/c it needs to call different routines of SCHISM. We could potential add a top level wrapper routine before triaging inside SCHISM.

DeniseWorthen commented 2 months ago

@uturuncoglu If I look at just the output in the regression test (coastal_ike_shinnecock_atm2sch2ww3_intel), it does not have the fix for the WW3 history files, so only the XY field is correct in the history file. Comparing the maximum values of WAVIMP_SW_WAVSUv with sxy, I get

Screenshot 2024-08-02 at 12 34 27 PM
uturuncoglu commented 2 months ago

@DeniseWorthen Yes. It does not have that fix. At this point, I am not pointing any feature branch in the model side. I was plaining to update the model once the PR is merged. It will also change the baseline in our side too.

DeniseWorthen commented 2 months ago

I've gotten the history file 2d radiation stresses to agree w/ the export stresses now. I've also got the significant wave height exporting OK.

@josephzhang8 I remember someone mentioning an issue w/ the units of the BHD field. Can you remind me of the correct unit?

josephzhang8 commented 2 months ago

That's great, thx @DeniseWorthen!

The unit for BHD should be m*m/s/s

DeniseWorthen commented 2 months ago

I've run into an issue w/ the fields

wave_flux_friction_x(1:np)=WW3_TBBX !Momentum flux due to bottom friction [m2/s2]
wave_flux_friction_y(1:np)=WW3_TBBY

It doesn't appear that the BT4 switch is being used, which is where the taubbl is calculated. Right now, it is showing up as zero in the history files.

https://github.com/oceanmodeling/WW3/blob/80a5ad544d8997931dde7fa7c8ebad56612ed299/model/src/w3srcemd.F90#L1296-L1299

DeniseWorthen commented 2 months ago

@josephzhang8 Do you know who created the mod_def file for the coastal_ike_shinnecock_atm2sch2ww3 configuration? I believe we will need a new mod_def generated with the BT4 switch present in order to have those fields available for export.

uturuncoglu commented 2 months ago

@DeniseWorthen I created it. I am not sure BT4 will be available through the switch file or it could be just added by running ww3_grid tool.

aliabdolali commented 2 months ago

@DeniseWorthen I created it. I am not sure BT4 will be available through the switch file or it could be just added by running ww3_grid tool.

BT4 should be in the switch and model should be recompiled, and ww3_grid should be re-executed to activate BT4 physics.

DeniseWorthen commented 2 months ago

@aliabdolali I see in /work2/noaa/nems/tufuk/RT/NEMSfv3gfs/input-data-20240501/WW3/Ike/atm2sch2wav/ww3_grid.inp a setting for SBT1. Do you know what we should use for BT4?

$ Bottom friction  - - - - - - - - - - - - - - - - - - - - - - - - - -
$   JONSWAP             : Namelist SBT1
$                           GAMMA   : As it says.
$  &SBT1 GAMMA = 0.15 /
$
DeniseWorthen commented 2 months ago

@uturuncoglu The current switch_meshcap_pdlib is currently using BT1

https://github.com/oceanmodeling/WW3/blob/80a5ad544d8997931dde7fa7c8ebad56612ed299/model/bin/switch_meshcap_pdlib#L1

If Coastal wants to use BT4, then you'll probably need a new switch file (maybe switch_meshcap_pdlib_bt4).