Skip to content

Commit

Permalink
Merge pull request #76 from jiandewang/feature/update-to-GFDL-20210914
Browse files Browse the repository at this point in the history 
update to gfdl 20210914
  • Loading branch information
@jiandewang
jiandewang authored Sep 21, 2021
2 parents f8a8e4c + 29016c2 commit 14ca4a1
Show file tree
Hide file tree
Showing 20 changed files with 1,284 additions and 233 deletions.
3 changes: 2 additions & 1 deletion config_src/drivers/mct_cap/mom_surface_forcing_mct.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1245,7 +1245,8 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt,
"If true, use a 2-dimensional gustiness supplied from "//&
"an input file", default=.false.)
call get_param(param_file, mdl, "GUST_CONST", CS%gust_const, &
"The background gustiness in the winds.", units="Pa", default=0.0)
"The background gustiness in the winds.", units="Pa", default=0.0, &
scale=US%kg_m3_to_R*US%m_s_to_L_T**2*US%L_to_Z)
if (CS%read_gust_2d) then
call get_param(param_file, mdl, "GUST_2D_FILE", gust_file, &
"The file in which the wind gustiness is found in "//&
Expand Down
74 changes: 41 additions & 33 deletions config_src/drivers/nuopc_cap/mom_cap.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,7 @@ module MOM_cap_mod
use MOM_grid, only: ocean_grid_type, get_global_grid_size
use MOM_ocean_model_nuopc, only: ice_ocean_boundary_type
use MOM_ocean_model_nuopc, only: ocean_model_restart, ocean_public_type, ocean_state_type
use MOM_ocean_model_nuopc, only: ocean_model_init_sfc
use MOM_ocean_model_nuopc, only: ocean_model_init_sfc, ocean_model_flux_init
use MOM_ocean_model_nuopc, only: ocean_model_init, update_ocean_model, ocean_model_end
use MOM_ocean_model_nuopc, only: get_ocean_grid, get_eps_omesh, query_ocean_state
use MOM_cap_time, only: AlarmInit
Expand Down Expand Up @@ -422,6 +422,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
logical :: isPresent, isPresentDiro, isPresentLogfile, isSet
logical :: existflag
logical :: use_waves ! If true, the wave modules are active.
character(len=40) :: wave_method ! Wave coupling method.
integer :: userRc
integer :: localPet
integer :: localPeCount
Expand Down Expand Up @@ -650,6 +651,10 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
ocean_public%is_ocean_pe = .true.
call ocean_model_init(ocean_public, ocean_state, time0, time_start, input_restart_file=trim(restartfiles))

! GMM, this call is not needed for NCAR. Check with EMC.
! If this can be deleted, perhaps we should also delete ocean_model_flux_init
call ocean_model_flux_init(ocean_state)

call ocean_model_init_sfc(ocean_state, ocean_public)

call mpp_get_compute_domain(ocean_public%domain, isc, iec, jsc, jec)
Expand All @@ -669,6 +674,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
Ice_ocean_boundary% seaice_melt_heat (isc:iec,jsc:jec),&
Ice_ocean_boundary% seaice_melt (isc:iec,jsc:jec), &
Ice_ocean_boundary% mi (isc:iec,jsc:jec), &
Ice_ocean_boundary% ice_fraction (isc:iec,jsc:jec), &
Ice_ocean_boundary% u10_sqr (isc:iec,jsc:jec), &
Ice_ocean_boundary% p (isc:iec,jsc:jec), &
Ice_ocean_boundary% lrunoff_hflx (isc:iec,jsc:jec), &
Ice_ocean_boundary% frunoff_hflx (isc:iec,jsc:jec), &
Expand All @@ -690,25 +697,32 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
Ice_ocean_boundary%seaice_melt = 0.0
Ice_ocean_boundary%seaice_melt_heat= 0.0
Ice_ocean_boundary%mi = 0.0
Ice_ocean_boundary%ice_fraction = 0.0
Ice_ocean_boundary%u10_sqr = 0.0
Ice_ocean_boundary%p = 0.0
Ice_ocean_boundary%lrunoff_hflx = 0.0
Ice_ocean_boundary%frunoff_hflx = 0.0
Ice_ocean_boundary%lrunoff = 0.0
Ice_ocean_boundary%frunoff = 0.0

call query_ocean_state(ocean_state, use_waves=use_waves)
call query_ocean_state(ocean_state, use_waves=use_waves, wave_method=wave_method)
if (use_waves) then
call query_ocean_state(ocean_state, NumWaveBands=Ice_ocean_boundary%num_stk_bands)
allocate ( Ice_ocean_boundary% ustk0 (isc:iec,jsc:jec), &
Ice_ocean_boundary% vstk0 (isc:iec,jsc:jec), &
Ice_ocean_boundary% ustkb (isc:iec,jsc:jec,Ice_ocean_boundary%num_stk_bands), &
Ice_ocean_boundary% vstkb (isc:iec,jsc:jec,Ice_ocean_boundary%num_stk_bands), &
Ice_ocean_boundary%stk_wavenumbers (Ice_ocean_boundary%num_stk_bands))
Ice_ocean_boundary%ustk0 = 0.0
Ice_ocean_boundary%vstk0 = 0.0
call query_ocean_state(ocean_state, WaveNumbers=Ice_ocean_boundary%stk_wavenumbers, unscale=.true.)
Ice_ocean_boundary%ustkb = 0.0
Ice_ocean_boundary%vstkb = 0.0
if (wave_method == "EFACTOR") then
allocate( Ice_ocean_boundary%lamult(isc:iec,jsc:jec) )
Ice_ocean_boundary%lamult = 0.0
else
allocate ( Ice_ocean_boundary% ustk0 (isc:iec,jsc:jec), &
Ice_ocean_boundary% vstk0 (isc:iec,jsc:jec), &
Ice_ocean_boundary% ustkb (isc:iec,jsc:jec,Ice_ocean_boundary%num_stk_bands), &
Ice_ocean_boundary% vstkb (isc:iec,jsc:jec,Ice_ocean_boundary%num_stk_bands), &
Ice_ocean_boundary%stk_wavenumbers (Ice_ocean_boundary%num_stk_bands))
Ice_ocean_boundary%ustk0 = 0.0
Ice_ocean_boundary%vstk0 = 0.0
call query_ocean_state(ocean_state, WaveNumbers=Ice_ocean_boundary%stk_wavenumbers, unscale=.true.)
Ice_ocean_boundary%ustkb = 0.0
Ice_ocean_boundary%vstkb = 0.0
endif
endif
! Consider adding this:
! if (.not.use_waves) Ice_ocean_boundary%num_stk_bands = 0
Expand All @@ -722,18 +736,6 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
call fld_list_add(fldsFrOcn_num, fldsFrOcn, trim(scalar_field_name), "will_provide")
end if

if (cesm_coupled) then
!call fld_list_add(fldsToOcn_num, fldsToOcn, "Sw_lamult" , "will provide")
!call fld_list_add(fldsToOcn_num, fldsToOcn, "Sw_ustokes" , "will provide")
!call fld_list_add(fldsToOcn_num, fldsToOcn, "Sw_vstokes" , "will provide")
!call fld_list_add(fldsToOcn_num, fldsToOcn, "Sw_hstokes" , "will provide")
!call fld_list_add(fldsToOcn_num, fldsToOcn, "Fioi_melth" , "will provide")
!call fld_list_add(fldsToOcn_num, fldsToOcn, "Fioi_meltw" , "will provide")
!call fld_list_add(fldsFrOcn_num, fldsFrOcn, "So_fswpen" , "will provide")
else
!call fld_list_add(fldsToOcn_num, fldsToOcn, "mass_of_overlying_sea_ice" , "will provide")
!call fld_list_add(fldsFrOcn_num, fldsFrOcn, "sea_lev" , "will provide")
endif

!--------- import fields -------------
call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_salt_rate" , "will provide") ! from ice
Expand All @@ -751,21 +753,27 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
call fld_list_add(fldsToOcn_num, fldsToOcn, "inst_pres_height_surface" , "will provide")
call fld_list_add(fldsToOcn_num, fldsToOcn, "Foxx_rofl" , "will provide") !-> liquid runoff
call fld_list_add(fldsToOcn_num, fldsToOcn, "Foxx_rofi" , "will provide") !-> ice runoff
call fld_list_add(fldsToOcn_num, fldsToOcn, "Si_ifrac" , "will provide") !-> ice fraction
call fld_list_add(fldsToOcn_num, fldsToOcn, "So_duu10n" , "will provide") !-> wind^2 at 10m
call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_fresh_water_to_ocean_rate", "will provide")
call fld_list_add(fldsToOcn_num, fldsToOcn, "net_heat_flx_to_ocn" , "will provide")
!These are not currently used and changing requires a nuopc dictionary change
!call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_runoff_heat_flx" , "will provide")
!call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_calving_heat_flx" , "will provide")
if (use_waves) then
if (Ice_ocean_boundary%num_stk_bands > 3) then
call MOM_error(FATAL, "Number of Stokes Bands > 3, NUOPC cap not set up for this")
if (wave_method == "EFACTOR") then
call fld_list_add(fldsToOcn_num, fldsToOcn, "Sw_lamult" , "will provide")
else
if (Ice_ocean_boundary%num_stk_bands > 3) then
call MOM_error(FATAL, "Number of Stokes Bands > 3, NUOPC cap not set up for this")
endif
call fld_list_add(fldsToOcn_num, fldsToOcn, "eastward_partitioned_stokes_drift_1" , "will provide")
call fld_list_add(fldsToOcn_num, fldsToOcn, "northward_partitioned_stokes_drift_1", "will provide")
call fld_list_add(fldsToOcn_num, fldsToOcn, "eastward_partitioned_stokes_drift_2" , "will provide")
call fld_list_add(fldsToOcn_num, fldsToOcn, "northward_partitioned_stokes_drift_2", "will provide")
call fld_list_add(fldsToOcn_num, fldsToOcn, "eastward_partitioned_stokes_drift_3" , "will provide")
call fld_list_add(fldsToOcn_num, fldsToOcn, "northward_partitioned_stokes_drift_3", "will provide")
endif
call fld_list_add(fldsToOcn_num, fldsToOcn, "eastward_partitioned_stokes_drift_1" , "will provide")
call fld_list_add(fldsToOcn_num, fldsToOcn, "northward_partitioned_stokes_drift_1", "will provide")
call fld_list_add(fldsToOcn_num, fldsToOcn, "eastward_partitioned_stokes_drift_2" , "will provide")
call fld_list_add(fldsToOcn_num, fldsToOcn, "northward_partitioned_stokes_drift_2", "will provide")
call fld_list_add(fldsToOcn_num, fldsToOcn, "eastward_partitioned_stokes_drift_3" , "will provide")
call fld_list_add(fldsToOcn_num, fldsToOcn, "northward_partitioned_stokes_drift_3", "will provide")
endif

!--------- export fields -------------
Expand Down Expand Up @@ -1103,7 +1111,7 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc)
k = k + 1 ! Increment position within gindex
if (mask(k) /= 0) then
mesh_areas(k) = dataPtr_mesh_areas(k)
model_areas(k) = ocean_grid%AreaT(i,j) / ocean_grid%Rad_Earth**2
model_areas(k) = ocean_grid%US%L_to_m**2 * ocean_grid%AreaT(i,j) / ocean_grid%Rad_Earth**2
mod2med_areacor(k) = model_areas(k) / mesh_areas(k)
med2mod_areacor(k) = mesh_areas(k) / model_areas(k)
end if
Expand Down
27 changes: 26 additions & 1 deletion config_src/drivers/nuopc_cap/mom_cap_methods.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -250,12 +250,37 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary,
!----
! Note - preset values to 0, if field does not exist in importState, then will simply return
! and preset value will be used

ice_ocean_boundary%mi(:,:) = 0._ESMF_KIND_R8
call state_getimport(importState, 'mass_of_overlying_ice', &
isc, iec, jsc, jec, ice_ocean_boundary%mi,rc=rc)
if (ChkErr(rc,__LINE__,u_FILE_u)) return

!----
! sea-ice fraction
!----
ice_ocean_boundary%ice_fraction(:,:) = 0._ESMF_KIND_R8
call state_getimport(importState, 'Si_ifrac', &
isc, iec, jsc, jec, ice_ocean_boundary%ice_fraction, rc=rc)
if (ChkErr(rc,__LINE__,u_FILE_u)) return

!----
! 10m wind squared
!----
ice_ocean_boundary%u10_sqr(:,:) = 0._ESMF_KIND_R8
call state_getimport(importState, 'So_duu10n', &
isc, iec, jsc, jec, ice_ocean_boundary%u10_sqr, rc=rc)
if (ChkErr(rc,__LINE__,u_FILE_u)) return

!----
! Langmuir enhancement factor
!----
if ( associated(ice_ocean_boundary%lamult) ) then
ice_ocean_boundary%lamult (:,:) = 0._ESMF_KIND_R8
call state_getimport(importState, 'Sw_lamult', &
isc, iec, jsc, jec, ice_ocean_boundary%lamult, rc=rc)
if (ChkErr(rc,__LINE__,u_FILE_u)) return
endif

!----
! Partitioned Stokes Drift Components
!----
Expand Down
31 changes: 24 additions & 7 deletions config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,6 @@ module MOM_ocean_model_nuopc
use MOM_time_manager, only : operator(/=), operator(<=), operator(>=)
use MOM_time_manager, only : operator(<), real_to_time_type, time_type_to_real
use time_interp_external_mod,only : time_interp_external_init
use MOM_tracer_flow_control, only : call_tracer_register, tracer_flow_control_init
use MOM_tracer_flow_control, only : call_tracer_flux_init
use MOM_unit_scaling, only : unit_scale_type
use MOM_variables, only : surface
Expand Down Expand Up @@ -147,6 +146,7 @@ module MOM_ocean_model_nuopc
integer :: nstep = 0 !< The number of calls to update_ocean.
logical :: use_ice_shelf !< If true, the ice shelf model is enabled.
logical,public :: use_waves !< If true use wave coupling.
character(len=40) :: wave_method !< Wave coupling method.

logical :: icebergs_alter_ocean !< If true, the icebergs can change ocean the
!! ocean dynamics and forcing fluxes.
Expand Down Expand Up @@ -242,14 +242,17 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i
!! tracer fluxes, and can be used to spawn related
!! internal variables in the ice model.
character(len=*), optional, intent(in) :: input_restart_file !< If present, name of restart file to read

! Local variables
real :: Rho0 ! The Boussinesq ocean density, in kg m-3.
real :: G_Earth ! The gravitational acceleration in m s-2.
real :: HFrz !< If HFrz > 0 (m), melt potential will be computed.
!! The actual depth over which melt potential is computed will
!! min(HFrz, OBLD), where OBLD is the boundary layer depth.
!! If HFrz <= 0 (default), melt potential will not be computed.
logical :: use_melt_pot!< If true, allocate melt_potential array
logical :: use_melt_pot !< If true, allocate melt_potential array
logical :: use_CFC !< If true, allocated arrays for surface CFCs.


! This include declares and sets the variable "version".
#include "version_variable.h"
Expand Down Expand Up @@ -368,13 +371,17 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i
use_melt_pot=.false.
endif

call get_param(param_file, mdl, "USE_CFC_CAP", use_CFC, &
default=.false., do_not_log=.true.)

! Consider using a run-time flag to determine whether to do the diagnostic
! vertical integrals, since the related 3-d sums are not negligible in cost.
call allocate_surface_state(OS%sfc_state, OS%grid, use_temperature, &
do_integrals=.true., gas_fields_ocn=gas_fields_ocn, use_meltpot=use_melt_pot)
do_integrals=.true., gas_fields_ocn=gas_fields_ocn, &
use_meltpot=use_melt_pot, use_cfcs=use_CFC)

call surface_forcing_init(Time_in, OS%grid, OS%US, param_file, OS%diag, &
OS%forcing_CSp, OS%restore_salinity, OS%restore_temp)
OS%forcing_CSp, OS%restore_salinity, OS%restore_temp, OS%use_waves)

if (OS%use_ice_shelf) then
call initialize_ice_shelf(param_file, OS%grid, OS%Time, OS%ice_shelf_CSp, &
Expand All @@ -386,8 +393,12 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i
call allocate_forcing_type(OS%grid, OS%fluxes, shelf=.true.)
endif

call allocate_forcing_type(OS%grid, OS%fluxes, waves=.true.)
call get_param(param_file, mdl, "USE_WAVES", OS%Use_Waves, &
"If true, enables surface wave modules.", default=.false.)
if (OS%Use_Waves) then
call get_param(param_file, mdl, "WAVE_METHOD", OS%wave_method, default="EMPTY", do_not_log=.true.)
endif
! MOM_wave_interface_init is called regardless of the value of USE_WAVES because
! it also initializes statistical waves.
call MOM_wave_interface_init(OS%Time, OS%grid, OS%GV, OS%US, param_file, OS%Waves, OS%diag)
Expand All @@ -413,6 +424,8 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i

endif

call extract_surface_state(OS%MOM_CSp, OS%sfc_state)

call close_param_file(param_file)
call diag_mediator_close_registration(OS%diag)

Expand Down Expand Up @@ -572,7 +585,9 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, &
call set_net_mass_forcing(OS%fluxes, OS%forces, OS%grid, OS%US)

if (OS%use_waves) then
call Update_Surface_Waves(OS%grid, OS%GV, OS%US, OS%time, ocean_coupling_time_step, OS%waves, OS%forces)
if (OS%wave_method /= "EFACTOR") then
call Update_Surface_Waves(OS%grid, OS%GV, OS%US, OS%time, ocean_coupling_time_step, OS%waves, OS%forces)
endif
endif

if (OS%nstep==0) then
Expand Down Expand Up @@ -1002,15 +1017,16 @@ end subroutine ocean_model_flux_init

!> This interface allows certain properties that are stored in the ocean_state_type to be
!! obtained.
subroutine query_ocean_state(OS, use_waves, NumWaveBands, Wavenumbers, unscale)
subroutine query_ocean_state(OS, use_waves, NumWaveBands, Wavenumbers, unscale, wave_method)
type(ocean_state_type), intent(in) :: OS !< The structure with the complete ocean state
logical, optional, intent(out) :: use_waves !< Indicates whether surface waves are in use
integer, optional, intent(out) :: NumWaveBands !< If present, this gives the number of
!! wavenumber partitions in the wave discretization
real, dimension(:), optional, intent(out) :: Wavenumbers !< If present, this gives the characteristic
!! wavenumbers of the wave discretization [m-1 or Z-1 ~> m-1]
logical, optional, intent(in) :: unscale !< If present and true, undo any dimensional
logical, optional, intent(in) :: unscale !< If present and true, undo any dimensional
!! rescaling and return dimensional values in MKS units
character(len=40), optional, intent(out) :: wave_method !< Wave coupling method.

logical :: undo_scaling
undo_scaling = .false. ; if (present(unscale)) undo_scaling = unscale
Expand All @@ -1022,6 +1038,7 @@ subroutine query_ocean_state(OS, use_waves, NumWaveBands, Wavenumbers, unscale)
elseif (present(Wavenumbers)) then
call query_wave_properties(OS%Waves, WaveNumbers=WaveNumbers)
endif
if (present(wave_method)) wave_method = OS%wave_method

end subroutine query_ocean_state

Expand Down
Loading

0 comments on commit 14ca4a1

Please sign in to comment.