From 21529b967d19098b28a80d6c66a0552cae03db80 Mon Sep 17 00:00:00 2001 From: Michael Levy Date: Fri, 12 Jul 2024 16:58:57 -0600 Subject: [PATCH] Changes following code review -- cleaned up a lot of comments and whitespace -- used source argument in more allocate statements, and deallocated more arrays -- 3D diags now have zl:mean in cell_methods attribute -- marbl_instances%domain%kmt is set once (during initialization) --- .../drivers/nuopc_cap/mom_cap_methods.F90 | 53 +++++----- .../nuopc_cap/mom_ocean_model_nuopc.F90 | 4 +- .../nuopc_cap/mom_surface_forcing_nuopc.F90 | 12 +-- .../solo_driver/MOM_surface_forcing.F90 | 30 +++--- config_src/external/MARBL/README.md | 4 +- .../external/MARBL/marbl_settings_mod.F90 | 10 -- src/core/MOM_forcing_type.F90 | 6 ++ src/tracer/MARBL_forcing_mod.F90 | 20 ++-- src/tracer/MARBL_tracers.F90 | 99 ++++++++++--------- src/tracer/MOM_tracer_registry.F90 | 7 +- 10 files changed, 120 insertions(+), 125 deletions(-) delete mode 100644 config_src/external/MARBL/marbl_settings_mod.F90 diff --git a/config_src/drivers/nuopc_cap/mom_cap_methods.F90 b/config_src/drivers/nuopc_cap/mom_cap_methods.F90 index 826969e05b..d5ec9dc259 100644 --- a/config_src/drivers/nuopc_cap/mom_cap_methods.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap_methods.F90 @@ -287,7 +287,6 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, ! nhx deposition !---- if (associated(ice_ocean_boundary%nhx_dep)) then - ice_ocean_boundary%nhx_dep(:,:) = 0._ESMF_KIND_R8 call state_getimport(importState, 'Faxa_ndep', & isc, iec, jsc, jec, ice_ocean_boundary%nhx_dep(:,:), & areacor=med2mod_areacor, esmf_ind=1, rc=rc) @@ -298,7 +297,6 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, ! noy deposition !---- if (associated(ice_ocean_boundary%noy_dep)) then - ice_ocean_boundary%noy_dep(:,:) = 0._ESMF_KIND_R8 call state_getimport(importState, 'Faxa_ndep', & isc, iec, jsc, jec, ice_ocean_boundary%noy_dep(:,:), & areacor=med2mod_areacor, esmf_ind=2, rc=rc) @@ -311,13 +309,11 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, ! in which the pointer(s) will not be associated !---- if (associated(ice_ocean_boundary%atm_co2_prog)) then - ice_ocean_boundary%atm_co2_prog(:,:) = 0._ESMF_KIND_R8 call state_getimport(importState, 'Sa_co2prog', & isc, iec, jsc, jec, ice_ocean_boundary%atm_co2_prog(:,:), rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return endif if (associated(ice_ocean_boundary%atm_co2_diag)) then - ice_ocean_boundary%atm_co2_diag(:,:) = 0._ESMF_KIND_R8 call state_getimport(importState, 'Sa_co2diag', & isc, iec, jsc, jec, ice_ocean_boundary%atm_co2_diag(:,:), rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & @@ -362,7 +358,6 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, ! dust flux from sea ice !---- if (associated(ice_ocean_boundary%seaice_dust_flux)) then - ice_ocean_boundary%seaice_dust_flux(:,:) = 0._ESMF_KIND_R8 call state_getimport(importState, 'Fioi_flxdst', & isc, iec, jsc, jec, ice_ocean_boundary%seaice_dust_flux, & areacor=med2mod_areacor, rc=rc) @@ -921,25 +916,25 @@ subroutine State_GetImport_2d(state, fldname, isc, iec, jsc, jec, output, do_sum ! determine output array and apply area correction if present n = 0 do j = jsc,jec - do i = isc,iec - n = n + 1 - if (do_sum_loc) then - if (present(areacor)) then - output(i,j) = output(i,j) + dataPtr1d(n) * areacor(n) - else - output(i,j) = output(i,j) + dataPtr1d(n) - endif + do i = isc,iec + n = n + 1 + if (do_sum_loc) then + if (present(areacor)) then + output(i,j) = output(i,j) + dataPtr1d(n) * areacor(n) else - if (present(areacor)) then - output(i,j) = dataPtr1d(n) * areacor(n) - else - output(i,j) = dataPtr1d(n) - endif + output(i,j) = output(i,j) + dataPtr1d(n) endif - enddo + else + if (present(areacor)) then + output(i,j) = dataPtr1d(n) * areacor(n) + else + output(i,j) = dataPtr1d(n) + endif + endif + enddo enddo - else if (geomtype == ESMF_GEOMTYPE_GRID) then + elseif (geomtype == ESMF_GEOMTYPE_GRID) then call state_getfldptr(state, trim(fldname), dataptr2d, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return @@ -948,15 +943,15 @@ subroutine State_GetImport_2d(state, fldname, isc, iec, jsc, jec, output, do_sum lbnd2 = lbound(dataPtr2d,2) do j = jsc, jec - j1 = j + lbnd2 - jsc - do i = isc, iec - i1 = i + lbnd1 - isc - if (do_sum_loc) then - output(i,j) = output(i,j) + dataPtr2d(i1,j1) - else - output(i,j) = dataPtr2d(i1,j1) - endif - enddo + j1 = j + lbnd2 - jsc + do i = isc, iec + i1 = i + lbnd1 - isc + if (do_sum_loc) then + output(i,j) = output(i,j) + dataPtr2d(i1,j1) + else + output(i,j) = dataPtr2d(i1,j1) + endif + enddo enddo endif diff --git a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 index df66dbadea..329f436e48 100644 --- a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 @@ -543,7 +543,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & if (OS%fluxes%fluxes_used) then - ! MNL: Necessary for averaging some MARBL forcing fields + ! enable_averages() is necessary to post forcing fields to diagnostics call enable_averages(dt_coupling, OS%Time + Ocean_coupling_time_step, OS%diag) if (do_thermo) & @@ -789,7 +789,7 @@ subroutine ocean_model_end(Ocean_sfc, Ocean_state, Time, write_restart) type(time_type), intent(in) :: Time !< The model time, used for writing restarts. logical, intent(in) :: write_restart !< true => write restart file - if (write_restart)call ocean_model_save_restart(Ocean_state, Time) + if (write_restart) call ocean_model_save_restart(Ocean_state, Time) call diag_mediator_end(Time, Ocean_state%diag, end_diag_manager=.true.) call MOM_end(Ocean_state%MOM_CSp) if (Ocean_state%use_ice_shelf) call ice_shelf_end(Ocean_state%Ice_shelf_CSp) 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 ba8fbc750c..3e8f80e265 100644 --- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 @@ -39,7 +39,7 @@ module MOM_surface_forcing_nuopc use user_revise_forcing, only : user_revise_forcing_CS use iso_fortran_env, only : int64 use MARBL_forcing_mod, only : marbl_forcing_CS, MARBL_forcing_init -use MARBL_forcing_mod, only : convert_marbl_IOB_to_forcings +use MARBL_forcing_mod, only : convert_driver_fields_to_forcings implicit none ; private @@ -585,11 +585,11 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, ! Copy MARBL-specific IOB fields into fluxes; also set some MARBL-specific forcings to other values ! (constants, values from netCDF, etc) - call convert_marbl_IOB_to_forcings(IOB%atm_fine_dust_flux, IOB%atm_coarse_dust_flux, & - IOB%seaice_dust_flux, IOB%atm_bc_flux, IOB%seaice_bc_flux, & - IOB%nhx_dep, IOB%noy_dep, IOB%atm_co2_prog, IOB%atm_co2_diag, & - IOB%afracr, IOB%swnet_afracr, IOB%ifrac_n, IOB%swpen_ifrac_n, & - Time, G, US, i0, j0, fluxes, CS%marbl_forcing_CSp) + call convert_driver_fields_to_forcings(IOB%atm_fine_dust_flux, IOB%atm_coarse_dust_flux, & + IOB%seaice_dust_flux, IOB%atm_bc_flux, IOB%seaice_bc_flux, & + IOB%nhx_dep, IOB%noy_dep, IOB%atm_co2_prog, IOB%atm_co2_diag, & + IOB%afracr, IOB%swnet_afracr, IOB%ifrac_n, IOB%swpen_ifrac_n, & + Time, G, US, i0, j0, fluxes, CS%marbl_forcing_CSp) ! wave to ocean coupling if ( associated(IOB%lamult)) then diff --git a/config_src/drivers/solo_driver/MOM_surface_forcing.F90 b/config_src/drivers/solo_driver/MOM_surface_forcing.F90 index a677896a90..3a8303e561 100644 --- a/config_src/drivers/solo_driver/MOM_surface_forcing.F90 +++ b/config_src/drivers/solo_driver/MOM_surface_forcing.F90 @@ -57,7 +57,7 @@ module MOM_surface_forcing use dumbbell_surface_forcing, only : dumbbell_surface_forcing_init, dumbbell_surface_forcing_CS use dumbbell_surface_forcing, only : dumbbell_buoyancy_forcing use MARBL_forcing_mod, only : marbl_forcing_CS, MARBL_forcing_init -use MARBL_forcing_mod, only : convert_marbl_IOB_to_forcings +use MARBL_forcing_mod, only : convert_driver_fields_to_forcings implicit none ; private @@ -1571,7 +1571,7 @@ subroutine MARBL_forcing_from_data_override(fluxes, day, G, US, CS) real, pointer, dimension(:,:) :: noy_dep =>NULL() !< Nitrogen deposition [kg/m^2/s ~> RZ/T] integer :: isc, iec, jsc, jec - ! Necessary null pointers for arguments to convert_marbl_IOB_to_forcings() + ! Necessary null pointers for arguments to convert_driver_fields_to_forcings() ! Since they are null, MARBL will not use multiple ice categories real, pointer, dimension(:,:) :: afracr =>NULL() real, pointer, dimension(:,:) :: swnet_afracr =>NULL() @@ -1595,17 +1595,9 @@ subroutine MARBL_forcing_from_data_override(fluxes, day, G, US, CS) atm_bc_flux (isc:iec,jsc:jec), & seaice_bc_flux (isc:iec,jsc:jec), & nhx_dep (isc:iec,jsc:jec), & - noy_dep (isc:iec,jsc:jec)) - - atm_co2_prog(:,:) = 0.0 - atm_co2_diag(:,:) = 0.0 - atm_fine_dust_flux(:,:) = 0.0 - atm_coarse_dust_flux(:,:) = 0.0 - atm_bc_flux(:,:) = 0.0 - seaice_dust_flux(:,:) = 0.0 - seaice_bc_flux(:,:) = 0.0 - nhx_dep(:,:) = 0.0 - noy_dep(:,:) = 0.0 + noy_dep (isc:iec,jsc:jec), & + source=0.0) + ! fluxes used directly as MARBL inputs ! (should be scaled) @@ -1614,7 +1606,7 @@ subroutine MARBL_forcing_from_data_override(fluxes, day, G, US, CS) ! fluxes used to compute MARBL inputs ! These are kept in physical units, and will be scaled appropriately in - ! convert_marbl_IOB_to_forcings() + ! convert_driver_fields_to_forcings() call data_override(G%Domain, 'atm_co2_prog', atm_co2_prog, day) call data_override(G%Domain, 'atm_co2_diag', atm_co2_diag, day) call data_override(G%Domain, 'atm_fine_dust_flux', atm_fine_dust_flux, day) @@ -1625,11 +1617,11 @@ subroutine MARBL_forcing_from_data_override(fluxes, day, G, US, CS) call data_override(G%Domain, 'nhx_dep', nhx_dep, day) call data_override(G%Domain, 'noy_dep', noy_dep, day) - call convert_marbl_IOB_to_forcings(atm_fine_dust_flux, atm_coarse_dust_flux, & - seaice_dust_flux, atm_bc_flux, seaice_bc_flux, & - nhx_dep, noy_dep, atm_co2_prog, atm_co2_diag, & - afracr, swnet_afracr, ifrac_n, swpen_ifrac_n, & - day, G, US, 0, 0, fluxes, CS%marbl_forcing_CSp) + call convert_driver_fields_to_forcings(atm_fine_dust_flux, atm_coarse_dust_flux, & + seaice_dust_flux, atm_bc_flux, seaice_bc_flux, & + nhx_dep, noy_dep, atm_co2_prog, atm_co2_diag, & + afracr, swnet_afracr, ifrac_n, swpen_ifrac_n, & + day, G, US, 0, 0, fluxes, CS%marbl_forcing_CSp) deallocate ( atm_co2_prog, & atm_co2_diag, & diff --git a/config_src/external/MARBL/README.md b/config_src/external/MARBL/README.md index 6be9e89f1b..f19f76dec8 100644 --- a/config_src/external/MARBL/README.md +++ b/config_src/external/MARBL/README.md @@ -1,5 +1,5 @@ -GFDL_ocean_BGC -============== +MARBL +===== These APIs reflect those for the MARBL library available at https://github.com/marbl-ecosys/MARBL diff --git a/config_src/external/MARBL/marbl_settings_mod.F90 b/config_src/external/MARBL/marbl_settings_mod.F90 deleted file mode 100644 index b56ba14fbe..0000000000 --- a/config_src/external/MARBL/marbl_settings_mod.F90 +++ /dev/null @@ -1,10 +0,0 @@ -!> A non-functioning template of MARBL_settings_mod -module marbl_settings_mod - - implicit none - private - - !> Molecular weight of iron - integer, public, parameter :: output_for_GCM_iopt_total_Chl_3d = 0 - -end module marbl_settings_mod \ No newline at end of file diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 6bbe14e802..6e4969142e 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -3617,6 +3617,12 @@ subroutine deallocate_forcing_type(fluxes) if (associated(fluxes%mass_berg)) deallocate(fluxes%mass_berg) if (associated(fluxes%ice_fraction)) deallocate(fluxes%ice_fraction) if (associated(fluxes%u10_sqr)) deallocate(fluxes%u10_sqr) + if (associated(fluxes%noy_dep)) deallocate(fluxes%noy_dep) + if (associated(fluxes%nhx_dep)) deallocate(fluxes%nhx_dep) + if (associated(fluxes%atm_co2)) deallocate(fluxes%atm_co2) + if (associated(fluxes%atm_alt_co2)) deallocate(fluxes%atm_alt_co2) + if (associated(fluxes%dust_flux)) deallocate(fluxes%dust_flux) + if (associated(fluxes%iron_flux)) deallocate(fluxes%iron_flux) if (associated(fluxes%fracr_cat)) deallocate(fluxes%fracr_cat) if (associated(fluxes%qsw_cat)) deallocate(fluxes%qsw_cat) diff --git a/src/tracer/MARBL_forcing_mod.F90 b/src/tracer/MARBL_forcing_mod.F90 index cbb939e8bf..ec78f54046 100644 --- a/src/tracer/MARBL_forcing_mod.F90 +++ b/src/tracer/MARBL_forcing_mod.F90 @@ -22,7 +22,7 @@ module MARBL_forcing_mod #include public :: MARBL_forcing_init -public :: convert_marbl_IOB_to_forcings +public :: convert_driver_fields_to_forcings !> Data type used to store diagnostic index returned from register_diag_field() !! For the forcing fields that can be written via post_data() @@ -161,7 +161,7 @@ subroutine MARBL_forcing_init(G, US, param_file, diag, day, inputdir, use_marbl, endif ! Register diagnostic fields for outputing forcing values - ! These fields are posted from convert_marbl_IOB_to_forcings(), and they are received + ! These fields are posted from convert_driver_fields_to_forcings(), and they are received ! in physical units so no conversion is necessary here. CS%diag_ids%atm_fine_dust = register_diag_field("ocean_model", "ATM_FINE_DUST_FLUX_CPL", & CS%diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid @@ -183,11 +183,11 @@ subroutine MARBL_forcing_init(G, US, param_file, diag, day, inputdir, use_marbl, end subroutine MARBL_forcing_init ! Note: ice fraction and u10_sqr are handled in mom_surface_forcing because of CFCs - subroutine convert_marbl_IOB_to_forcings(atm_fine_dust_flux, atm_coarse_dust_flux, & - seaice_dust_flux, atm_bc_flux, seaice_bc_flux, & - nhx_dep, noy_dep, atm_co2_prog, atm_co2_diag, & - afracr, swnet_afracr, ifrac_n, & - swpen_ifrac_n, Time, G, US, i0, j0, fluxes, CS) + subroutine convert_driver_fields_to_forcings(atm_fine_dust_flux, atm_coarse_dust_flux, & + seaice_dust_flux, atm_bc_flux, seaice_bc_flux, & + nhx_dep, noy_dep, atm_co2_prog, atm_co2_diag, & + afracr, swnet_afracr, ifrac_n, & + swpen_ifrac_n, Time, G, US, i0, j0, fluxes, CS) real, dimension(:,:), pointer, intent(in) :: atm_fine_dust_flux !< atmosphere fine dust flux from IOB !! [kg m-2 s-1] @@ -223,7 +223,7 @@ subroutine convert_marbl_IOB_to_forcings(atm_fine_dust_flux, atm_coarse_dust_flu real :: seaice_fe_bioavail_frac !< TODO: define this (local) term real :: iron_flux_conversion !< TODO: define this (local) term real :: ndep_conversion !< Combination of unit conversion factors for rescaling - !! nitrogen deposition [kg(N) m-2 s-1 ~> mol L-2 T-1] + !! nitrogen deposition [kg(N) m-2 s-1 ~> mol m-3 Z T-1] if (.not. CS%use_marbl_tracers) return @@ -373,6 +373,6 @@ subroutine convert_marbl_IOB_to_forcings(atm_fine_dust_flux, atm_coarse_dust_flu enddo; enddo endif - end subroutine convert_marbl_IOB_to_forcings + end subroutine convert_driver_fields_to_forcings -end module MARBL_forcing_mod \ No newline at end of file +end module MARBL_forcing_mod diff --git a/src/tracer/MARBL_tracers.F90 b/src/tracer/MARBL_tracers.F90 index 69dd25ff50..672d28cc07 100644 --- a/src/tracer/MARBL_tracers.F90 +++ b/src/tracer/MARBL_tracers.F90 @@ -152,11 +152,11 @@ module MARBL_tracers restoring_timescale_z_edges !< The depths of the cell interfaces in the tracer restoring timescale file [Z ~> m] real, allocatable, dimension(:) :: & restoring_timescale_dz !< The thickness of the cell layers in the tracer restoring timescale file [H ~> m] - character(len=14) :: restoring_rtau_source !< location of tracer restoring timescale data - !! valid values: file, grid_dependent + character(len=14) :: restoring_I_tau_source !< location of inverse restoring timescale data + !! valid values: file, grid_dependent character(len=200) :: restoring_file !< name of [netCDF] file containing tracer restoring data type(remapping_CS) :: restoring_remapCS !< Remapping parameters and work arrays for tracer restoring / timescale - character(len=200) :: restoring_rtau_file !< name of [netCDF] file containing tracer restoring timescale + character(len=200) :: restoring_I_tau_file !< name of [netCDF] file containing inverse restoring timescale character(len=35) :: marbl_settings_file !< name of [text] file containing MARBL settings real :: bot_flux_mix_thickness !< for bottom flux -> tendency conversion, assume uniform mixing over @@ -253,8 +253,8 @@ module MARBL_tracers type(external_field), allocatable :: id_tracer_restoring(:) !< id number for time_interp_external integer, allocatable :: tracer_restoring_ind(:) !< index of MARBL forcing field to copy !! per-tracer restoring field into - integer, allocatable :: tracer_rtau_ind(:) !< index of MARBL forcing field to copy per-tracer - !! restoring timescale into + integer, allocatable :: tracer_I_tau_ind(:) !< index of MARBL forcing field to copy per-tracer + !! inverse restoring timescale into !> Memory for storing river fluxes, tracer restoring fields, and abiotic forcing real, allocatable :: d14c(:,:) !< d14c forcing for abiotic DIC and carbon isotope tracer modules @@ -262,7 +262,7 @@ module MARBL_tracers !! (needs to be time-integrated when passed to function!) !! (dims: i, j, tracer) [conc m s-1] character(len=15), allocatable :: tracer_restoring_varname(:) !< name of variable being restored - real, allocatable :: rtau(:,:,:) !< 1 / restoring time scale for marbl tracers (dims: i, j, k) [1/s] + real, allocatable :: I_tau(:,:,:) !< inverse restoring timescale for marbl tracers (dims: i, j, k) [1/s] real, allocatable, dimension(:,:,:,:) :: restoring_in !< Restoring fields read from file !! (dims: i, j, restoring_nz, restoring_cnt) [tracer units] @@ -275,8 +275,8 @@ module MARBL_tracers ! TODO: create generic 3D forcing input type to read z coordinate + values real :: fesedflux_scale_factor !< scale factor for iron sediment flux integer :: fesedflux_nz !< number of levels in iron sediment flux file - real, allocatable, dimension(:,:,:) :: fesedflux_in !< Field to read iron sediment flux into [conc m/s] - real, allocatable, dimension(:,:,:) :: feventflux_in !< Field to read iron vent flux into [conc m/s] + real, allocatable, dimension(:,:,:) :: fesedflux_in !< Field to read iron sediment flux into [conc m s-1] + real, allocatable, dimension(:,:,:) :: feventflux_in !< Field to read iron vent flux into [conc m s-1] real, allocatable, dimension(:) :: & fesedflux_z_edges !< The depths of the cell interfaces in the input data [Z ~> m] ! TODO: this thickness does not need to be 3D, but that's a problem for future Mike @@ -303,7 +303,7 @@ subroutine configure_MARBL_tracers(GV, US, param_file, CS) character(len=256) :: log_message character(len=256) :: marbl_in_line(1) character(len=256) :: forcing_sname, field_source - integer :: m, n, nz, marbl_settings_in, read_error, rtau_count, fi + integer :: m, n, nz, marbl_settings_in, read_error, I_tau_count, fi logical :: chl_from_file, forcing_processed nz = GV%ke marbl_settings_in = 615 @@ -381,6 +381,8 @@ subroutine configure_MARBL_tracers(GV, US, param_file, CS) gcm_num_elements_surface_flux = 1, & ! FIXME: change to number of grid cells on MPI task gcm_delta_z = GV%sInterface(2:nz+1) - GV%sInterface(1:nz), gcm_zw = GV%sInterface(2:nz+1), & gcm_zt = GV%sLayer, unit_system_opt = "mks", lgcm_has_global_ops = .false.) ! FIXME: add global ops + ! Regardless of vertical grid, MOM6 will always use GV%ke levels in all columns + MARBL_instances%domain%kmt = GV%ke if (MARBL_instances%StatusLog%labort_marbl) & call MARBL_instances%StatusLog%log_error_trace("MARBL_instances%init", & "configure_MARBL_tracers") @@ -486,9 +488,9 @@ subroutine configure_MARBL_tracers(GV, US, param_file, CS) allocate(CS%tracer_restoring_varname(CS%ntr), source=' ') ! gfortran 13.2 bug? ! source = '' does not blank out strings allocate(CS%tracer_restoring_ind(CS%ntr), source=-1) - allocate(CS%tracer_rtau_ind(CS%ntr), source=-1) + allocate(CS%tracer_I_tau_ind(CS%ntr), source=-1) CS%restore_count = 0 - rtau_count = 0 + I_tau_count = 0 do m=1,size(MARBL_instances%interior_tendency_forcings) select case (trim(MARBL_instances%interior_tendency_forcings(m)%metadata%varname)) case('Dust Flux') @@ -522,8 +524,8 @@ subroutine configure_MARBL_tracers(GV, US, param_file, CS) fi = index(MARBL_instances%interior_tendency_forcings(m)%metadata%varname, & 'Restoring Inverse Timescale') if (fi > 0) then - rtau_count = rtau_count + 1 - CS%tracer_rtau_ind(rtau_count) = m + I_tau_count = I_tau_count + 1 + CS%tracer_I_tau_ind(I_tau_count) = m else write(log_message, "(A,1X,A)") & trim(MARBL_instances%interior_tendency_forcings(m)%metadata%varname), & @@ -621,8 +623,8 @@ function register_MARBL_tracers(HI, GV, US, param_file, CS, tr_Reg, restart_CS) endif ! ** Scale factor for FESEDFLUX call get_param(param_file, mdl, "MARBL_FESEDFLUX_SCALE_FACTOR", CS%fesedflux_scale_factor, & - "Conversion factor between FESEDFLUX file and MARBL units", & - units="umol/m^2/d -> mmol/m^2/s", default=0.001/86400.) + "Conversion factor between FESEDFLUX file units and MARBL units", & + units="umol m-1 d-1 -> mmol m-2 s-1", default=0.001/86400.) ! ** River fluxes call get_param(param_file, mdl, "READ_RIV_FLUXES", CS%read_riv_fluxes, & @@ -706,8 +708,8 @@ function register_MARBL_tracers(HI, GV, US, param_file, CS, tr_Reg, restart_CS) case("file") call get_param(param_file, mdl, "MARBL_TRACER_RESTORING_FILE", CS%restoring_file, & "File containing fields to restore MARBL tracers towards") - call get_param(param_file, mdl, "MARBL_TRACER_RESTORING_RTAU_SOURCE", & - CS%restoring_rtau_source, "Source of data for 1/timescale for restoring MARBL tracers") + call get_param(param_file, mdl, "MARBL_TRACER_RESTORING_I_TAU_SOURCE", & + CS%restoring_I_tau_source, "Source of data for inverse timescale for restoring MARBL tracers") ! Initialize remapping type call initialize_remapping(CS%restoring_remapCS, 'PCM', boundary_extrapolation=.false., answer_date=99991231) @@ -721,13 +723,13 @@ function register_MARBL_tracers(HI, GV, US, param_file, CS, tr_Reg, restart_CS) CS%restoring_dz(k) = (CS%restoring_z_edges(k) - CS%restoring_z_edges(kbot)) * GV%Z_to_H enddo - select case(CS%restoring_rtau_source) + select case(CS%restoring_I_tau_source) case("file") - call get_param(param_file, mdl, "MARBL_TRACER_RESTORING_RTAU_FILE", & - CS%restoring_rtau_file, & + call get_param(param_file, mdl, "MARBL_TRACER_RESTORING_I_TAU_FILE", & + CS%restoring_I_tau_file, & "File containing the inverse timescale for restoring MARBL tracers") ! Set up array for thicknesses in restoring timescale file - call read_Z_edges(CS%restoring_rtau_file, "RTAU", CS%restoring_timescale_z_edges, & + call read_Z_edges(CS%restoring_I_tau_file, "I_TAU", CS%restoring_timescale_z_edges, & CS%restoring_timescale_nz, restoring_timescale_has_edges, & restoring_timescale_use_missing, restoring_timescale_missing, scale=US%m_to_Z) allocate(CS%restoring_timescale_dz(CS%restoring_timescale_nz)) @@ -737,8 +739,8 @@ function register_MARBL_tracers(HI, GV, US, param_file, CS, tr_Reg, restart_CS) CS%restoring_timescale_z_edges(kbot)) * GV%Z_to_H enddo case DEFAULT - write(log_message, "(3A)") "'", trim(CS%restoring_rtau_source), & - "' is not a valid option for MARBL_TRACER_RESTORING_RTAU_SOURCE" + write(log_message, "(3A)") "'", trim(CS%restoring_I_tau_source), & + "' is not a valid option for MARBL_TRACER_RESTORING_I_TAU_SOURCE" call MOM_error(FATAL, log_message) end select case DEFAULT @@ -1067,7 +1069,7 @@ subroutine initialize_MARBL_tracers(restart, day, G, GV, US, h, param_file, diag endif ! Initialize external field for restoring - if (CS%restoring_rtau_source == "file") then + if (CS%restoring_I_tau_source == "file") then select case(CS%restoring_source) case("file") ! Set up array for reading in raw restoring data @@ -1077,10 +1079,10 @@ subroutine initialize_MARBL_tracers(restart, day, G, GV, US, h, param_file, diag trim(CS%tracer_restoring_varname(m)), domain=G%Domain%mpp_domain) enddo end select - select case(CS%restoring_rtau_source) + select case(CS%restoring_I_tau_source) case("file") - allocate(CS%rtau(SZI_(G), SZJ_(G), CS%restoring_timescale_nz), source=0.) - call MOM_read_data(CS%restoring_rtau_file, "RTAU", CS%rtau(:,:,:), G%Domain) + allocate(CS%I_tau(SZI_(G), SZJ_(G), CS%restoring_timescale_nz), source=0.) + call MOM_read_data(CS%restoring_I_tau_file, "RTAU", CS%I_tau(:,:,:), G%Domain) end select endif @@ -1116,13 +1118,24 @@ subroutine register_MARBL_diags(MARBL_diags, diag, day, G, id_diags) else ! 3D field ! TODO: MARBL should provide v_extensive through MARBL_diags ! (for now, FESEDFLUX is the only one that should be true) - id_diags(m)%id = register_diag_field("ocean_model", & - trim(MARBL_diags%diags(m)%short_name), & - diag%axesTL, & ! T=> tracer grid? L => layer center - day, & - trim(MARBL_diags%diags(m)%long_name), & - trim(MARBL_diags%diags(m)%units), & - v_extensive=(trim(MARBL_diags%diags(m)%short_name).eq."FESEDFLUX")) + ! Also, known issue where passing v_extensive=.false. isn't + ! treated the same as not passing v_extensive + if (trim(MARBL_diags%diags(m)%short_name).eq."FESEDFLUX") then + id_diags(m)%id = register_diag_field("ocean_model", & + trim(MARBL_diags%diags(m)%short_name), & + diag%axesTL, & ! T=> tracer grid? L => layer center + day, & + trim(MARBL_diags%diags(m)%long_name), & + trim(MARBL_diags%diags(m)%units), & + v_extensive=.true.) + else + id_diags(m)%id = register_diag_field("ocean_model", & + trim(MARBL_diags%diags(m)%short_name), & + diag%axesTL, & ! T=> tracer grid? L => layer center + day, & + trim(MARBL_diags%diags(m)%long_name), & + trim(MARBL_diags%diags(m)%units)) + endif if (id_diags(m)%id > 0) allocate(id_diags(m)%field_3d(SZI_(G),SZJ_(G), SZK_(G)), source=0.0) endif enddo @@ -1231,7 +1244,6 @@ subroutine MARBL_tracers_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke if (.not.associated(CS)) return - if (CS%ntr < 1) return ! (1) Compute surface fluxes ! FIXME: MARBL can handle computing surface fluxes for all columns simultaneously @@ -1418,8 +1430,6 @@ subroutine MARBL_tracers_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, if (G%mask2dT(i,j) == 0) cycle ! ii. Set up vertical domain and bot_flux_to_tend - ! TODO: why doesn't ke have (i,j) coordinate? - MARBL_instances%domain%kmt = GV%ke ! Calculate depth of interface by building up thicknesses from the bottom (top interface is always 0) ! MARBL wants this to be positive-down zi(GV%ke) = G%bathyT(i,j) @@ -1494,11 +1504,11 @@ subroutine MARBL_tracers_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, MARBL_instances%interior_tendency_forcings(CS%tracer_restoring_ind(m))%field_1d(1,:)) if (m==1) then call remapping_core_h(CS%restoring_remapCS, CS%restoring_timescale_nz, & - CS%restoring_timescale_dz(:), CS%rtau(i,j,:), GV%ke, dz(:), & - MARBL_instances%interior_tendency_forcings(CS%tracer_rtau_ind(m))%field_1d(1,:)) + CS%restoring_timescale_dz(:), CS%I_tau(i,j,:), GV%ke, dz(:), & + MARBL_instances%interior_tendency_forcings(CS%tracer_I_tau_ind(m))%field_1d(1,:)) else - MARBL_instances%interior_tendency_forcings(CS%tracer_rtau_ind(m))%field_1d(1,:) = & - MARBL_instances%interior_tendency_forcings(CS%tracer_rtau_ind(1))%field_1d(1,:) + MARBL_instances%interior_tendency_forcings(CS%tracer_I_tau_ind(m))%field_1d(1,:) = & + MARBL_instances%interior_tendency_forcings(CS%tracer_I_tau_ind(1))%field_1d(1,:) endif enddo @@ -1527,7 +1537,6 @@ subroutine MARBL_tracers_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, MARBL_instances%interior_tendency_forcings(CS%remin_scalef_ind)%field_1d(1,:) = 1. ! * Column Tracers - ! NOTE: POP averages previous two timesteps, should we do that too? do m=1,CS%ntr MARBL_instances%tracers(m, :) = CS%tracer_data(m)%tr(i,j,:) enddo @@ -1686,7 +1695,7 @@ subroutine MARBL_tracers_set_forcing(day_start, G, CS) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(MARBL_tracers_CS), pointer :: CS !< The control structure returned by a - ! These are Fortran parameters in POP + ! Fraction of river nutrients in refractory pools real, parameter :: DONriv_refract = 0.1 real, parameter :: DOCriv_refract = 0.2 real, parameter :: DOPriv_refract = 0.025 @@ -1961,10 +1970,10 @@ subroutine MARBL_tracers_end(CS) if (allocated(CS%RIV_FLUXES)) deallocate(CS%RIV_FLUXES) if (allocated(CS%SFO)) deallocate(CS%SFO) if (allocated(CS%tracer_restoring_ind)) deallocate(CS%tracer_restoring_ind) - if (allocated(CS%tracer_rtau_ind)) deallocate(CS%tracer_rtau_ind) + if (allocated(CS%tracer_I_tau_ind)) deallocate(CS%tracer_I_tau_ind) if (allocated(CS%fesedflux_in)) deallocate(CS%fesedflux_in) if (allocated(CS%feventflux_in)) deallocate(CS%feventflux_in) - if (allocated(CS%rtau)) deallocate(CS%rtau) + if (allocated(CS%I_tau)) deallocate(CS%I_tau) deallocate(CS) endif end subroutine MARBL_tracers_end diff --git a/src/tracer/MOM_tracer_registry.F90 b/src/tracer/MOM_tracer_registry.F90 index 8a2ad77a4c..c7d11b6030 100644 --- a/src/tracer/MOM_tracer_registry.F90 +++ b/src/tracer/MOM_tracer_registry.F90 @@ -726,7 +726,7 @@ subroutine post_tracer_transport_diagnostics(G, GV, Reg, h_diag, diag) intent(in) :: h_diag !< Layer thicknesses on which to post fields [H ~> m or kg m-2] type(diag_ctrl), intent(in) :: diag !< structure to regulate diagnostic output - integer :: i, j, k, is, ie, js, je, nz, m,khi + integer :: i, j, k, is, ie, js, je, nz, m, khi real :: frac_under_100m(SZI_(G),SZJ_(G),SZK_(GV)) real :: work2d(SZI_(G),SZJ_(G)), ztop(SZI_(G),SZJ_(G)), zbot(SZI_(G),SZJ_(G)) type(tracer_type), pointer :: Tr=>NULL() @@ -735,6 +735,9 @@ subroutine post_tracer_transport_diagnostics(G, GV, Reg, h_diag, diag) ! If any tracers are posting 100m vertical integrals, compute weights frac_under_100m(:,:,:) = 0.0 + ! khi will be the largest layer index corresponding where ztop < 100m and ztop >= 100m + ! in any column (we can reduce computation of 100m integrals by only looping through khi + ! rather than GV%ke) khi = 0 do m=1,Reg%ntr ; if (Reg%Tr(m)%registry_diags) then Tr => Reg%Tr(m) @@ -761,7 +764,7 @@ subroutine post_tracer_transport_diagnostics(G, GV, Reg, h_diag, diag) do m=1,Reg%ntr ; if (Reg%Tr(m)%registry_diags) then Tr => Reg%Tr(m) - if (Tr%id_tr_post_horzn> 0) call post_data(Tr%id_tr_post_horzn, Tr%t, diag) ! Tr%t is tracer concentration + if (Tr%id_tr_post_horzn> 0) call post_data(Tr%id_tr_post_horzn, Tr%t, diag) if (Tr%id_adx > 0) call post_data(Tr%id_adx, Tr%ad_x, diag, alt_h=h_diag) if (Tr%id_ady > 0) call post_data(Tr%id_ady, Tr%ad_y, diag, alt_h=h_diag) if (Tr%id_dfx > 0) call post_data(Tr%id_dfx, Tr%df_x, diag, alt_h=h_diag)