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.

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

NCAR to main candidate branch (2022-06-03) #1571

Conversation

gustavo-marques
Copy link
Collaborator

In addition to bug fixes and general code improvements, this PR includes the following features:

  • Option to use the wave-averaged Boussinesq equations;
  • Option to pass enthalpy fluxes via NUOPC cap;
  • Improvements in the correction to the computation of cell-averaged density in various parts of the model;
  • Add KPP nonlocal term to passive tracers.

I have verified that this PR does not change answers in MOM6-examples using intel/19.1.1.


New parameters

  • ENTHALPY_FROM_COUPLER (only via NUOPC cap)
  • CFC_APPLY_NONLOCAL_TRANSPORT
  • STANLEY parameterization
    • USE_STANLEY_ISO
    • STANLEY_COEFF
    • USE_STANLEY_GM
    • STOCH_EOS
    • USE_STANLEY_ISO
    • USE_STANLEY_ML
    • USE_STANLEY_PGF

Obsolete parameters

  • STANLEY_PRM_DET_COEFF
  • PGF_STANLEY_T2_DET_COEFF

New diagnostics

  • NUOPC cap with ENTHALPY_FROM_COUPLER=True
    • heat_content_evap
    • total_heat_content_evap
  • Wave-averaged Boussinesq equations
    • CAu_Stokes
    • CAv_Stokes
    • dvdt_Stokes
    • dudt_Stokes
    • 3d_stokes_y_from_ddt
    • 3d_stokes_x_from_ddt
    • PFv_Stokes
    • PFu_Stokes

Obsolete diagnostics

  • heat_content_icemelt
  • total_heat_content_icemelt

The enthalpy associated with the mass from sea ice formation/melting is already accounted for in seaice_melt_heat. A note explaining this has been added to the code;


Features


Bug fixes


Testing, and cleanup


Contributors:

breichl and others added 30 commits May 21, 2021 12:29
- Add Craik-Leibovich terms, vorticity and dynamic pressure components.
- Add exploratory Stokes ddt version in MOM_wave_interface and MOM.F90
- Add tendency diagnostics for Stokes terms w/ control and option for diagnostic only Stokes force computation in MOM_CoriolisAdv.
- This only makes sense to add on the dynamics timestep.
Fix tidal mixing coeff diags when they are unavailable.
With these changes, arbitrary number of of partitioned stokes drift components
may be imported from ww3 when coupled within cesm. After the planned unification
of ww3 nuopc caps of CESM and EMC, the old (fixed) approach may be removed.
profile to restart file when surfbands wave coupling mode is on.
add the optional waves_csp arg to initialize_MOM
calls so as to ensure that the waves_csp is allocated.
This commit also has several misc doxygen fixes.
in MOM_CoriolisAdv. Also remove the Passive_Stokes_VF
block since it's never set to true.
…ments incorrectly override dvdx and dudy when surfbands wave coupling is on.
- Adding logical to prevent using terms if not allocated
- Adding optional diagnostic only output of term
- Moving time increment to within dynamics step instead of in main step MOOM loop.
- Correcting time tendency term to reflect time between Stokes drift updates and extent of dynamics time loop.
@alperaltuntas
Copy link
Collaborator

Regarding the answer changes observed by EMC, I am inviting @breichl to chime on why it may be the case. Brandon, this PR brings in the refactorings and bug fixes we have worked on several months ago. I believe EMC folks have the following wave set up, and they are seeing answer changes:

WAVE_METHOD = "SURFACE_BANDS"
SURFBAND_SOURCE = "COUPLER"
STK_BAND_COUPLER = 3
SURFBAND_WAVENUMBERS = 0.04, 0.11, 0.3305

@jiandewang
Copy link
Collaborator

Regarding the answer changes observed by EMC, I am inviting @breichl to chime on why it may be the case. Brandon, this PR brings in the refactorings and bug fixes we have worked on several months ago. I believe EMC folks have the following wave set up, and they are seeing answer changes:

WAVE_METHOD = "SURFACE_BANDS"
SURFBAND_SOURCE = "COUPLER"
STK_BAND_COUPLER = 3
SURFBAND_WAVENUMBERS = 0.04, 0.11, 0.3305

yes we are using the above settings, see
https://github.com/ufs-community/ufs-weather-model/blob/develop/tests/parm/MOM_input_template_025#L757-L776

Address comment from reviewer by adding units to covTS and varS.
* Add ``implicit none ; private`` to this module;
* Put module variables into the control structure for this module;
* Add the description of the units for all real variables;
* Add a consistent two-point indent throughout the module .

TODO:
Without further modifications, adding ``private`` to the control
structure of this module will break the model. Currently, MOM.F90
needs access to ``use_stoch_eos``, ``stanley_coeff``, and some of
the diagnostic ids.
@jiandewang
Copy link
Collaborator

Regarding the answer changes observed by EMC, I am inviting @breichl to chime on why it may be the case. Brandon, this PR brings in the refactorings and bug fixes we have worked on several months ago. I believe EMC folks have the following wave set up, and they are seeing answer changes:

WAVE_METHOD = "SURFACE_BANDS"
SURFBAND_SOURCE = "COUPLER"
STK_BAND_COUPLER = 3
SURFBAND_WAVENUMBERS = 0.04, 0.11, 0.3305

I further tested UFS by turning wave off and got same answer as using main branch, so the issue is in the updated wave code

@alperaltuntas
Copy link
Collaborator

@jiandewang I have set up an experiment that is similar to your configuration as close as possible, and pinpointed all of the bugfixes that cause the answer changes. Applying these changes in the main branch results in bitwise identical ocean.stats files when compared to the results of this candidate branch.

(1) In mom_surface_forcing_nuopc.F90:

diff --git a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90
index da4c15e52..8b33e91e5 100644
--- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90
+++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90
@@ -902,7 +902,8 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
         forces%ustkb(i,j,istk) = IOB%ustkb(i-I0,j-J0,istk)
         forces%vstkb(i,j,istk) = IOB%vstkb(i-I0,j-J0,istk)
       enddo; enddo
-      call pass_vector(forces%ustkb(:,:,istk),forces%vstkb(:,:,istk), G%domain )
+      call pass_var(forces%ustkb(:,:,istk), G%domain )
+      call pass_var(forces%vstkb(:,:,istk), G%domain )
     enddo
   endif

(2) In MOM_wave_interface.F90:

diff --git a/src/user/MOM_wave_interface.F90 b/src/user/MOM_wave_interface.F90
index ab8a693ba..d887ecfa7 100644
--- a/src/user/MOM_wave_interface.F90
+++ b/src/user/MOM_wave_interface.F90
@@ -333,7 +333,7 @@ subroutine MOM_wave_interface_init(time, G, GV, US, param_file, CS, diag )
                  "or the model will fail.", default=1)
       allocate( CS%WaveNum_Cen(CS%NumBands), source=0.0 )
       allocate( CS%STKx0(G%isdB:G%iedB,G%jsd:G%jed,CS%NumBands), source=0.0 )
-      allocate( CS%STKy0(G%isdB:G%iedB,G%jsd:G%jed,CS%NumBands), source=0.0 )
+      allocate( CS%STKy0(G%isd:G%ied,G%jsdB:G%jedB,CS%NumBands), source=0.0 )
       CS%PartitionMode = 0
       call get_param(param_file, mdl, "SURFBAND_WAVENUMBERS", CS%WaveNum_Cen, &
            "Central wavenumbers for surface Stokes drift bands.", &
@@ -449,26 +449,24 @@ subroutine query_wave_properties(CS, NumBands, WaveNumbers, US)
 end subroutine query_wave_properties
 
 !> Subroutine that handles updating of surface wave/Stokes drift related properties
-subroutine Update_Surface_Waves(G, GV, US, Day, dt, CS, forces)
+subroutine Update_Surface_Waves(G, GV, US, Time_present, dt, CS, forces)
   type(wave_parameters_CS), pointer    :: CS  !< Wave parameter Control structure
   type(ocean_grid_type), intent(inout) :: G   !< Grid structure
   type(verticalGrid_type), intent(in)  :: GV  !< Vertical grid structure
   type(unit_scale_type),   intent(in)  :: US  !< A dimensional unit scaling type
-  type(time_type),         intent(in)  :: Day !< Current model time
+  type(time_type),         intent(in)  :: Time_present !< Current model time
   type(time_type),         intent(in)  :: dt  !< Timestep as a time-type
   type(mech_forcing),      intent(in), optional  :: forces !< MOM_forcing_type
   ! Local variables
+  type(time_type) :: Stokes_Time
   integer :: ii, jj, b
-  type(time_type) :: Day_Center
-
-  ! Computing central time of time step
-  Day_Center = Day + DT/2
 
   if (CS%WaveMethod == TESTPROF) then
     ! Do nothing
   elseif (CS%WaveMethod == SURFBANDS) then
     if (CS%DataSource == DATAOVR) then
-      call Surface_Bands_by_data_override(day_center, G, GV, US, CS)
+      Stokes_Time = Time_present + dt/2
+      call Surface_Bands_by_data_override(Stokes_Time, G, GV, US, CS)
     elseif (CS%DataSource == COUPLER) then
       if (.not.present(FORCES)) then
         call MOM_error(FATAL,"The option SURFBAND = COUPLER can not be used with "//&

@jiandewang
Copy link
Collaborator

@alperaltuntas thanks, I will try that and get back to you. Just for confirmation: STOKES_VF, STOKES_PGF and STOKES_DDT need to be True in my test, right ?

@alperaltuntas
Copy link
Collaborator

Just for confirmation: STOKES_VF, STOKES_PGF and STOKES_DDT need to be True in my test, right ?

No, please set those to False as in your original configuration.

@jiandewang
Copy link
Collaborator

jiandewang commented Jun 15, 2022

@alperaltuntas
for changes here
remove call pass_vector(forces%ustkb(:,:,istk),forces%vstkb(:,:,istk), G%domain )
add call pass_var(forces%ustkb(:,:,istk), G%domain )
add call pass_var(forces%vstkb(:,:,istk), G%domain )
they are already in the candidate branch
https://github.com/gustavo-marques/MOM6/blob/dev-ncar-main-candidate-2022-06-03/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90#L918-L919

I am bit confused here. I guess your git diff with main branch is based on your own branch which is not based on candidate branch

@alperaltuntas
Copy link
Collaborator

alperaltuntas commented Jun 15, 2022

I applied those diffs (bugfixes) to the main branch so as to make the main branch results the same as the results of this candidate. If you'd like to do the opposite, i.e., make the results of this candidate branch the same as the results of the main branch, then I guess you can revert those changes in the candidate branch instead. Essentially, what I did was to fix the bugs in main branch by applying those diffs. You can instead do the opposite and reintroduce those bugs to this main candidate to make the results the same.

@jiandewang
Copy link
Collaborator

I applied those diffs (bugfixes) to the main branch so as to make the main branch results the same as the results of this candidate. If you'd like to do the opposite, i.e., make the results of this candidate branch the same as the results of the main branch, then I guess you can revert those changes in the candidate branch instead. Essentially, what I did was to fix the bugs in main branch by applying those diffs. You can instead do the opposite and reintroduce those bugs to this main candidate to make the results the same.

@alperaltuntas: still not clear here:
(1) in your test, you modified main branch code (bugfixes) and got the same result as candidate, but they differ from main branch, in other words this candidate changed NCAR results, right ?
(2) can you generate a branch based on candidate but revert the changes for me to test ?
(3) is there a way to have a flag like "keep-wave-bug" in the code so that we can have options ?

@alperaltuntas
Copy link
Collaborator

alperaltuntas commented Jun 15, 2022

@jiandewang, please apply the following change in this candidate branch and see if you can get bitwise identical results when compared to the main.

diff --git a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90
index c45a59c22..9e14d678f 100644
--- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90
+++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90
@@ -915,8 +915,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
         forces%ustkb(i,j,istk) = IOB%ustkb(i-I0,j-J0,istk)
         forces%vstkb(i,j,istk) = IOB%vstkb(i-I0,j-J0,istk)
       enddo; enddo
-      call pass_var(forces%ustkb(:,:,istk), G%domain )
-      call pass_var(forces%vstkb(:,:,istk), G%domain )
+      call pass_vector(forces%ustkb(:,:,istk),forces%vstkb(:,:,istk), G%domain )
     enddo
   endif

Let me know if you still encounter answer changes.

@jiandewang
Copy link
Collaborator

@jiandewang, please apply the following change in this candidate branch and see if you can get bitwise identical results when compared to the main.

diff --git a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90
index c45a59c22..9e14d678f 100644
--- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90
+++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90
@@ -915,8 +915,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
         forces%ustkb(i,j,istk) = IOB%ustkb(i-I0,j-J0,istk)
         forces%vstkb(i,j,istk) = IOB%vstkb(i-I0,j-J0,istk)
       enddo; enddo
-      call pass_var(forces%ustkb(:,:,istk), G%domain )
-      call pass_var(forces%vstkb(:,:,istk), G%domain )
+      call pass_vector(forces%ustkb(:,:,istk),forces%vstkb(:,:,istk), G%domain )
     enddo
   endif

Let me know if you still encounter answer changes.

now I got the same answer as from main branch. Why ustkb and vstkb need to use pass_var instead of pass as vector ? I need document this in UFS as we will need to create a new baseline.

@alperaltuntas
Copy link
Collaborator

The rationale is that both ustkb and vstkb are received from the wave component at T points, so during halo updates, they need to be treated as scalars and not vectors.

@jiandewang
Copy link
Collaborator

jiandewang commented Jun 16, 2022

The rationale is that both ustkb and vstkb are received from the wave component at T points, so during halo updates, they need to be treated as scalars and not vectors.

now the final question is do we need to have a wavebug flag here ? It will be fine either way for UFS.

@alperaltuntas
Copy link
Collaborator

@jiandewang If that's okay, my preference would be to not a have a wavebug flag, for the sake of keeping things more concise.

@jiandewang
Copy link
Collaborator

@jiandewang If that's okay, my preference would be to not a have a wavebug flag, for the sake of keeping things more concise.

@alperaltuntas sounds good, one of my team member is doing test in his system, will approve this PR once got confirmation.

@Hallberg-NOAA
Copy link
Collaborator

@alperaltuntas and @jiandewang, in the past we have opted not to introduce bug-flag in cases where the code was demonstrably not-self-consistent (e.g., changing layouts or changing between symmetric and non-symmetric memory leading to answer changes), because in that case it is not clear what previous solution needs to be recovered. In this case I could easily envision that there are internal self-inconsistencies with the previous code that fall into this category.

Copy link
Collaborator

@sanAkel sanAkel left a comment

Choose a reason for hiding this comment

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

This code in this PR gives the same answers within GEOS.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.