diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 4bdd7bf17d..a13d4121ab 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1424,29 +1424,6 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) call find_obsolete_params(param_file) -#ifdef SYMMETRIC_MEMORY_ - symmetric = .true. -#else - symmetric = .false. -#endif -#ifdef STATIC_MEMORY_ - call MOM_domains_init(G%domain, param_file, symmetric=symmetric, & - static_memory=.true., NIHALO=NIHALO_, NJHALO=NJHALO_, & - NIGLOBAL=NIGLOBAL_, NJGLOBAL=NJGLOBAL_, NIPROC=NIPROC_, & - NJPROC=NJPROC_) -#else - call MOM_domains_init(G%domain, param_file, symmetric=symmetric) -#endif - call callTree_waypoint("domains initialized (initialize_MOM)") - - call MOM_checksums_init(param_file) - - call diag_mediator_infrastructure_init() - call MOM_io_init(param_file) - call MOM_grid_init(G, param_file) - call verticalGridInit( param_file, CS%GV ) - GV => CS%GV - ! Read relevant parameters and write them to the model log. call log_version(param_file, "MOM", version, "") call get_param(param_file, "MOM", "VERBOSITY", verbosity, & @@ -1656,9 +1633,41 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) if (CS%adiabatic .and. CS%bulkmixedlayer) call MOM_error(FATAL, & "MOM: ADIABATIC and BULKMIXEDLAYER can not both be defined.") + ! Set up the model domain and grids. +#ifdef SYMMETRIC_MEMORY_ + symmetric = .true. +#else + symmetric = .false. +#endif +#ifdef STATIC_MEMORY_ + call MOM_domains_init(G%domain, param_file, symmetric=symmetric, & + static_memory=.true., NIHALO=NIHALO_, NJHALO=NJHALO_, & + NIGLOBAL=NIGLOBAL_, NJGLOBAL=NJGLOBAL_, NIPROC=NIPROC_, & + NJPROC=NJPROC_) +#else + call MOM_domains_init(G%domain, param_file, symmetric=symmetric) +#endif + call callTree_waypoint("domains initialized (initialize_MOM)") + + call MOM_checksums_init(param_file) + + call diag_mediator_infrastructure_init() + call MOM_io_init(param_file) + call MOM_grid_init(G, param_file) + + call create_dyn_horgrid(dG, G%HI) + dG%first_direction = G%first_direction + dG%bathymetry_at_vel = G%bathymetry_at_vel + call clone_MOM_domain(G%Domain, dG%Domain) + + call verticalGridInit( param_file, CS%GV ) + GV => CS%GV + dG%g_Earth = GV%g_Earth + + ! Allocate the auxiliary non-symmetric domain for debugging or I/O purposes. - if (CS%debug .or. G%symmetric) & - call clone_MOM_domain(G%Domain, G%Domain_aux, symmetric=.false.) + if (CS%debug .or. dG%symmetric) & + call clone_MOM_domain(dG%Domain, dG%Domain_aux, symmetric=.false.) call MOM_timing_init(CS) @@ -1666,11 +1675,11 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) ! Copy a common variable from the vertical grid to the horizontal grid. ! Consider removing this later? - G%ke = GV%ke +! G%ke = GV%ke - is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke - isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed - IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB + is = dG%isc ; ie = dG%iec ; js = dG%jsc ; je = dG%jec ; nz = GV%ke + isd = dG%isd ; ied = dG%ied ; jsd = dG%jsd ; jed = dG%jed + IsdB = dG%IsdB ; IedB = dG%IedB ; JsdB = dG%JsdB ; JedB = dG%JedB ! Allocate and initialize space for primary MOM variables. ALLOC_(CS%u(IsdB:IedB,jsd:jed,nz)) ; CS%u(:,:,:) = 0.0 @@ -1785,17 +1794,28 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) call callTree_waypoint("restart registration complete (initialize_MOM)") call cpu_clock_begin(id_clock_MOM_init) - call MOM_initialize_fixed(G, CS%OBC, param_file, write_geom_files, dirs%output_directory) + call MOM_initialize_fixed(dG, CS%OBC, param_file, write_geom_files, dirs%output_directory) call callTree_waypoint("returned from MOM_initialize_fixed() (initialize_MOM)") call MOM_initialize_coord(GV, param_file, write_geom_files, & - dirs%output_directory, CS%tv, G%max_depth) + dirs%output_directory, CS%tv, dG%max_depth) call callTree_waypoint("returned from MOM_initialize_coord() (initialize_MOM)") if (CS%use_ALE_algorithm) then - call ALE_init(param_file, GV, G%max_depth, CS%ALE_CSp) + call ALE_init(param_file, GV, dG%max_depth, CS%ALE_CSp) call callTree_waypoint("returned from ALE_init() (initialize_MOM)") endif + ! Shift from using the temporary dynamic grid type to using the final (potentially + ! static) ocean grid type. + ! call clone_MOM_domain(dG%Domain, CS%G%Domain) + ! call MOM_grid_init(CS%G, param_file) + + call copy_dyngrid_to_MOM_grid(dg, G) + ! Allocate the auxiliary non-symmetric domain for debugging or I/O purposes. + if (CS%debug .or. G%symmetric) & + call clone_MOM_domain(G%Domain, G%Domain_aux, symmetric=.false.) + G%ke = GV%ke + call MOM_initialize_state(CS%u, CS%v, CS%h, CS%tv, Time, G, GV, param_file, & dirs, CS%restart_CSp, CS%ALE_CSp, CS%tracer_Reg, & CS%sponge_CSp, CS%ALE_sponge_CSp, CS%OBC, Time_in) @@ -1805,20 +1825,6 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) ! From this point, there may be pointers being set, so the final grid type ! that will persist through the run has to be used. - ! Shift from using the temporary dynamic grid type to using the final (potentially - ! static) ocean grid type. - ! call clone_MOM_domain(dG%Domain, CS%G%Domain) - ! call MOM_grid_init(CS%G, param_file) - ! call copy_dyngrid_to_MOM_grid(dg, CS%G) - ! Allocate the auxiliary non-symmetric domain for debugging or I/O purposes. - ! if (CS%debug .or. CS%G%symmetric) & - ! call clone_MOM_domain(CS%G%Domain, CS%G%Domain_aux, symmetric=.false.) - - ! ! Copy a common variable from the vertical grid to the horizontal grid. - ! ! Consider removing this later? - ! CS%G%ke = GV%ke - ! G => CS%G - if (test_grid_copy) then ! Copy the data from the temporary grid to the dyn_hor_grid to CS%G. call create_dyn_horgrid(dG, G%HI) @@ -1885,12 +1891,12 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) diag => CS%diag ! Initialize the diag mediator. - call diag_mediator_init(G, param_file, diag, doc_file_dir=dirs%output_directory) + call diag_mediator_init(G, GV%ke, param_file, diag, doc_file_dir=dirs%output_directory) ! Initialize the diagnostics mask arrays. ! This step has to be done after call to MOM_initialize_state ! and before MOM_diagnostics_init - call diag_masks_set(G, CS%missing, diag) + call diag_masks_set(G, GV%ke, CS%missing, diag) ! Set up a pointers h within diag mediator control structure, ! this needs to occur _after_ CS%h has been allocated. diff --git a/src/core/MOM_CoriolisAdv.F90 b/src/core/MOM_CoriolisAdv.F90 index be81ed8a8e..5f229bfebf 100644 --- a/src/core/MOM_CoriolisAdv.F90 +++ b/src/core/MOM_CoriolisAdv.F90 @@ -875,7 +875,7 @@ subroutine CoriolisAdv_init(Time, G, param_file, diag, AD, CS) CS%diag => diag ; CS%Time => Time ! Read all relevant parameters and write them to the model log. - call log_version(param_file, mod, version) + call log_version(param_file, mod, version, "") call get_param(param_file, mod, "NOSLIP", CS%no_slip, & "If true, no slip boundary conditions are used; otherwise \n"//& "free slip boundary conditions are assumed. The \n"//& diff --git a/src/core/MOM_PressureForce.F90 b/src/core/MOM_PressureForce.F90 index d1a502e8ca..5edfdd6606 100644 --- a/src/core/MOM_PressureForce.F90 +++ b/src/core/MOM_PressureForce.F90 @@ -144,7 +144,7 @@ subroutine PressureForce_init(Time, G, GV, param_file, diag, CS, tides_CSp) else ; allocate(CS) ; endif ! Read all relevant parameters and write them to the model log. - call log_version(param_file, mod, version) + call log_version(param_file, mod, version, "") call get_param(param_file, mod, "ANALYTIC_FV_PGF", CS%Analytic_FV_PGF, & "If true the pressure gradient forces are calculated \n"//& "with a finite volume form that analytically integrates \n"//& diff --git a/src/core/MOM_PressureForce_Montgomery.F90 b/src/core/MOM_PressureForce_Montgomery.F90 index 10aeeea1c0..eea083f9ab 100644 --- a/src/core/MOM_PressureForce_Montgomery.F90 +++ b/src/core/MOM_PressureForce_Montgomery.F90 @@ -967,7 +967,7 @@ subroutine PressureForce_Mont_init(Time, G, GV, param_file, diag, CS, tides_CSp) endif mod = "MOM_PressureForce_Mont" - call log_version(param_file, mod, version) + call log_version(param_file, mod, version, "") call get_param(param_file, mod, "RHO_0", CS%Rho0, & "The mean ocean density used with BOUSSINESQ true to \n"//& "calculate accelerations and the mass for conservation \n"//& diff --git a/src/core/MOM_PressureForce_analytic_FV.F90 b/src/core/MOM_PressureForce_analytic_FV.F90 index 4e547154bc..2d92fa6a3f 100644 --- a/src/core/MOM_PressureForce_analytic_FV.F90 +++ b/src/core/MOM_PressureForce_analytic_FV.F90 @@ -871,7 +871,7 @@ subroutine PressureForce_AFV_init(Time, G, GV, param_file, diag, CS, tides_CSp) endif mod = "MOM_PressureForce_AFV" - call log_version(param_file, mod, version) + call log_version(param_file, mod, version, "") call get_param(param_file, mod, "RHO_0", CS%Rho0, & "The mean ocean density used with BOUSSINESQ true to \n"//& "calculate accelerations and the mass for conservation \n"//& diff --git a/src/core/MOM_continuity.F90 b/src/core/MOM_continuity.F90 index cefb7ee938..912dde87a7 100644 --- a/src/core/MOM_continuity.F90 +++ b/src/core/MOM_continuity.F90 @@ -190,7 +190,7 @@ subroutine continuity_init(Time, G, GV, param_file, diag, CS) allocate(CS) ! Read all relevant parameters and write them to the model log. - call log_version(param_file, mod, version) + call log_version(param_file, mod, version, "") call get_param(param_file, mod, "CONTINUITY_SCHEME", tmpstr, & "CONTINUITY_SCHEME selects the discretization for the \n"//& "continuity solver. The only valid value currently is: \n"//& diff --git a/src/core/MOM_continuity_PPM.F90 b/src/core/MOM_continuity_PPM.F90 index 84f24e4ced..86d05a15b4 100644 --- a/src/core/MOM_continuity_PPM.F90 +++ b/src/core/MOM_continuity_PPM.F90 @@ -192,7 +192,7 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, CS, uhbt, vhbt, OBC, apply_OBC_u_flather_east = .false. ; apply_OBC_u_flather_west = .false. apply_OBC_v_flather_north = .false. ; apply_OBC_v_flather_south = .false. - if (present(OBC)) then ; if (associated(OBC)) then + if (present(OBC)) then ; if (associated(OBC)) then ; if (OBC%this_pe) then apply_OBC_u_flather_east = OBC%apply_OBC_u_flather_east apply_OBC_u_flather_west = OBC%apply_OBC_u_flather_west apply_OBC_v_flather_north = OBC%apply_OBC_v_flather_north @@ -202,7 +202,7 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, CS, uhbt, vhbt, OBC, if (apply_OBC_u_flather_east .or. apply_OBC_u_flather_west .or. & apply_OBC_v_flather_north .or. apply_OBC_v_flather_south) & h_input(:,:,:) = hin(:,:,:) - endif ; endif + endif ; endif ; endif if (present(visc_rem_u) .neqv. present(visc_rem_v)) call MOM_error(FATAL, & "MOM_continuity_PPM: Either both visc_rem_u and visc_rem_v or neither"// & @@ -424,9 +424,9 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, CS, LB, uhbt, OBC, & use_visc_rem = present(visc_rem_u) apply_OBC_u = .false. ; set_BT_cont = .false. if (present(BT_cont)) set_BT_cont = (associated(BT_cont)) - if (present(OBC)) then ; if (associated(OBC)) then + if (present(OBC)) then ; if (associated(OBC)) then ; if (OBC%this_pe) then apply_OBC_u = OBC%apply_OBC_u - endif ; endif + endif ; endif ; endif ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh ; nz = G%ke CFL_dt = CS%CFL_limit_adjust / dt @@ -1181,9 +1181,9 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, CS, LB, vhbt, OBC, & use_visc_rem = present(visc_rem_v) apply_OBC_v = .false. ; set_BT_cont = .false. if (present(BT_cont)) set_BT_cont = (associated(BT_cont)) - if (present(OBC)) then ; if (associated(OBC)) then + if (present(OBC)) then ; if (associated(OBC)) then ; if (OBC%this_pe) then apply_OBC_v = OBC%apply_OBC_v - endif ; endif + endif ; endif ; endif ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh ; nz = G%ke CFL_dt = CS%CFL_limit_adjust / dt @@ -2163,7 +2163,7 @@ subroutine continuity_PPM_init(Time, G, GV, param_file, diag, CS) allocate(CS) ! Read all relevant parameters and write them to the model log. - call log_version(param_file, mod, version) + call log_version(param_file, mod, version, "") call get_param(param_file, mod, "MONOTONIC_CONTINUITY", CS%monotonic, & "If true, CONTINUITY_PPM uses the Colella and Woodward \n"//& "monotonic limiter. The default (false) is to use a \n"//& diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 5b58b10009..357ef2cc52 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -1,5 +1,5 @@ ! This file is part of MOM6. See LICENSE.md for the license. -!> Controls where open boundary conditions are applied +!> Controls where open boundary conditions are applied module MOM_open_boundary ! This file is part of MOM6. See LICENSE.md for the license. @@ -11,6 +11,7 @@ module MOM_open_boundary use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, WARNING use MOM_file_parser, only : get_param, log_version, param_file_type, log_param use MOM_grid, only : ocean_grid_type +use MOM_dyn_horgrid, only : dyn_horgrid_type use MOM_io, only : EAST_FACE, NORTH_FACE use MOM_io, only : slasher, read_data use MOM_tracer_registry, only : add_tracer_OBC_values, tracer_registry_type @@ -100,6 +101,7 @@ module MOM_open_boundary real :: rx_max !< The maximum magnitude of the baroclinic radiation !! velocity (or speed of characteristics), in m s-1. The !! default value is 10 m s-1. + logical :: this_pe !< Is there an open boundary on this tile? end type ocean_OBC_type integer :: id_clock_pass @@ -112,11 +114,9 @@ module MOM_open_boundary !> Enables OBC module and reads configuration parameters subroutine open_boundary_config(G, param_file, OBC) - type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(dyn_horgrid_type), intent(in) :: G !< Ocean grid structure type(param_file_type), intent(in) :: param_file !< Parameter file handle type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure - ! Local variables - logical :: flather_east, flather_west, flather_north, flather_south allocate(OBC) @@ -164,9 +164,9 @@ end subroutine open_boundary_config !> Initialize open boundary control structure subroutine open_boundary_init(G, param_file, OBC) - type(ocean_grid_type), intent(in) :: G !< Ocean grid structure - type(param_file_type), intent(in) :: param_file !< Parameter file handle - type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(param_file_type), intent(in) :: param_file !< Parameter file handle + type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure ! Local variables if (.not.associated(OBC)) return @@ -244,7 +244,7 @@ end subroutine open_boundary_end !> Sets the slope of bathymetry normal to an open bounndary to zero. subroutine open_boundary_impose_normal_slope(OBC, G, depth) type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure - type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(dyn_horgrid_type), intent(in) :: G !< Ocean grid structure real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: depth !< Bathymetry at h-points ! Local variables integer :: i, j @@ -270,7 +270,7 @@ end subroutine open_boundary_impose_normal_slope !> Reconcile masks and open boundaries, deallocate OBC on PEs where it is not needed subroutine open_boundary_impose_land_mask(OBC, G) type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure - type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(dyn_horgrid_type), intent(in) :: G !< Ocean grid structure ! Local variables integer :: i, j logical :: any_U, any_V @@ -306,9 +306,9 @@ subroutine open_boundary_impose_land_mask(OBC, G) ! bathymetry inside the boundary was do shallow and flagged as land. if (OBC%OBC_mask_u(I,j)) any_U = .true. enddo ; enddo - if (.not. any_U) then - deallocate(OBC%OBC_mask_u) - endif +! if (.not. any_U) then +! deallocate(OBC%OBC_mask_u) +! endif endif any_V = .false. @@ -316,12 +316,14 @@ subroutine open_boundary_impose_land_mask(OBC, G) do J=G%JsdB,G%JedB ; do i=G%isd,G%ied if (OBC%OBC_mask_v(i,J)) any_V = .true. enddo ; enddo - if (.not. any_V) then - deallocate(OBC%OBC_mask_v) - endif +! if (.not. any_V) then +! deallocate(OBC%OBC_mask_v) +! endif endif - if (.not.(any_U .or. any_V)) call open_boundary_dealloc(OBC) +! if (.not.(any_U .or. any_V)) call open_boundary_dealloc(OBC) + OBC%this_pe = .true. + if (.not.(any_U .or. any_V)) OBC%this_pe = .false. end subroutine open_boundary_impose_land_mask @@ -442,7 +444,7 @@ end subroutine Radiation_Open_Bdry_Conds !> Sets the domain boundaries as Flather open boundaries using the original !! Flather run-time logicals subroutine set_Flather_positions(G, OBC) - type(ocean_grid_type), intent(inout) :: G + type(dyn_horgrid_type), intent(inout) :: G type(ocean_OBC_type), pointer :: OBC ! Local variables integer :: east_boundary, west_boundary, north_boundary, south_boundary @@ -566,8 +568,6 @@ subroutine set_Flather_positions(G, OBC) enddo ; enddo endif - ! If there are no OBC points on this PE, there is no reason to keep the OBC - ! type, and it could be deallocated. end subroutine set_Flather_positions !> Sets the initial definitions of the characteristic open boundary conditions. @@ -595,8 +595,8 @@ subroutine set_Flather_data(OBC, tv, h, G, PF, tracer_Reg) real, pointer, dimension(:,:,:) :: & OBC_T_u => NULL(), & ! These arrays should be allocated and set to OBC_T_v => NULL(), & ! specify the values of T and S that should come - OBC_S_u => NULL(), & - OBC_S_v => NULL() + OBC_S_u => NULL(), & + OBC_S_v => NULL() is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed @@ -639,7 +639,7 @@ subroutine set_Flather_data(OBC, tv, h, G, PF, tracer_Reg) if (.not.associated(OBC%eta_outer_v)) then allocate(OBC%eta_outer_v(isd:ied,JsdB:JedB)) ; OBC%eta_outer_v(:,:) = 0.0 endif - + if (read_OBC_uv) then call read_data(filename, 'ubt', OBC%ubt_outer, & domain=G%Domain%mpp_domain, position=EAST_FACE) diff --git a/src/core/MOM_transcribe_grid.F90 b/src/core/MOM_transcribe_grid.F90 index a3d6e75375..def4c1f606 100644 --- a/src/core/MOM_transcribe_grid.F90 +++ b/src/core/MOM_transcribe_grid.F90 @@ -139,6 +139,7 @@ subroutine copy_dyngrid_to_MOM_grid(dG, oG) oG%south_lat = dG%south_lat ; oG%west_lon = dG%west_lon oG%len_lat = dG%len_lat ; oG%len_lon = dG%len_lon oG%Rad_Earth = dG%Rad_Earth ; oG%max_depth = dG%max_depth + oG%g_Earth = dG%g_Earth ! Update the halos in case the dynamic grid has smaller halos than the ocean grid. call pass_var(oG%areaT, oG%Domain) @@ -288,7 +289,8 @@ subroutine copy_MOM_grid_to_dyngrid(oG, dG) dG%areaT_global = oG%areaT_global ; dG%IareaT_global = oG%IareaT_global dG%south_lat = oG%south_lat ; dG%west_lon = oG%west_lon dG%len_lat = oG%len_lat ; dG%len_lon = oG%len_lon - dG%Rad_Earth = oG%Rad_Earth ; dG%max_depth = oG%max_depth + dG%Rad_Earth = oG%Rad_Earth ; dG%max_depth = oG%max_depth + dG%g_Earth = oG%g_Earth ! Update the halos in case the dynamic grid has smaller halos than the ocean grid. call pass_var(dG%areaT, dG%Domain) diff --git a/src/diagnostics/MOM_diag_to_Z.F90 b/src/diagnostics/MOM_diag_to_Z.F90 index 8eb34d9604..a1ddc3adb2 100644 --- a/src/diagnostics/MOM_diag_to_Z.F90 +++ b/src/diagnostics/MOM_diag_to_Z.F90 @@ -988,7 +988,7 @@ subroutine MOM_diag_to_Z_init(Time, G, GV, param_file, diag, CS) CS%diag => diag ! Read parameters and write them to the model log. - call log_version(param_file, mod, version) + call log_version(param_file, mod, version, "") ! Read in z-space info from a NetCDF file. call get_param(param_file, mod, "Z_OUTPUT_GRID_FILE", zgrid_file, & "The file that specifies the vertical grid for \n"//& diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 5dda3dce55..86e6812969 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -326,7 +326,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, fluxes, & ! volume mean potential temperature if (CS%id_thetaoga>0) then - thetaoga = global_volume_mean(tv%T, h, G) + thetaoga = global_volume_mean(tv%T, h, G, GV) call post_data(CS%id_thetaoga, thetaoga, CS%diag) endif @@ -341,7 +341,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, fluxes, & ! volume mean salinity if (CS%id_soga>0) then - soga = global_volume_mean(tv%S, h, G) + soga = global_volume_mean(tv%S, h, G, GV) call post_data(CS%id_soga, soga, CS%diag) endif @@ -356,13 +356,13 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, fluxes, & ! layer mean potential temperature if (CS%id_temp_layer_ave>0) then - temp_layer_ave = global_layer_mean(tv%T, h, G) + temp_layer_ave = global_layer_mean(tv%T, h, G, GV) call post_data_1d_k(CS%id_temp_layer_ave, temp_layer_ave, CS%diag) endif ! layer mean salinity if (CS%id_salt_layer_ave>0) then - salt_layer_ave = global_layer_mean(tv%S, h, G) + salt_layer_ave = global_layer_mean(tv%S, h, G, GV) call post_data_1d_k(CS%id_salt_layer_ave, salt_layer_ave, CS%diag) endif diff --git a/src/diagnostics/MOM_sum_output.F90 b/src/diagnostics/MOM_sum_output.F90 index d8101ca9c3..d0aa3cd65c 100644 --- a/src/diagnostics/MOM_sum_output.F90 +++ b/src/diagnostics/MOM_sum_output.F90 @@ -592,7 +592,8 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, CS, tracer_CSp) energypath_nc = trim(CS%energyfile) // ".nc" if (day > CS%Start_time) then call reopen_file(CS%fileenergy_nc, trim(energypath_nc), vars, & - num_nc_fields, G, CS%fields, SINGLE_FILE, CS%timeunit, GV=GV) + num_nc_fields, CS%fields, SINGLE_FILE, CS%timeunit, & + G=G, GV=GV) else call create_file(CS%fileenergy_nc, trim(energypath_nc), vars, & num_nc_fields, CS%fields, SINGLE_FILE, CS%timeunit, & diff --git a/src/equation_of_state/MOM_EOS.F90 b/src/equation_of_state/MOM_EOS.F90 index e05d868c88..980074db69 100644 --- a/src/equation_of_state/MOM_EOS.F90 +++ b/src/equation_of_state/MOM_EOS.F90 @@ -311,34 +311,34 @@ end subroutine calculate_2_densities subroutine int_specific_vol_dp(T, S, p_t, p_b, alpha_ref, HI, EOS, & dza, intp_dza, intx_dza, inty_dza, halo_size) !> The horizontal index structure - type(hor_index_type), intent(in) :: HI + type(hor_index_type), intent(in) :: HI !> Potential temperature referenced to the surface (degC) - real, dimension(SZDI_(HI),SZDJ_(HI)), intent(in) :: T + real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), intent(in) :: T !> Salinity (PSU) - real, dimension(SZDI_(HI),SZDJ_(HI)), intent(in) :: S + real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), intent(in) :: S !> Pressure at the top of the layer in Pa. - real, dimension(SZDI_(HI),SZDJ_(HI)), intent(in) :: p_t + real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), intent(in) :: p_t !> Pressure at the bottom of the layer in Pa. - real, dimension(SZDI_(HI),SZDJ_(HI)), intent(in) :: p_b + real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), intent(in) :: p_b !> A mean specific volume that is subtracted out to reduce the magnitude of !! each of the integrals, m3 kg-1. The calculation is mathematically identical !! with different values of alpha_ref, but this reduces the effects of roundoff. - real, intent(in) :: alpha_ref + real, intent(in) :: alpha_ref !> Equation of state structure - type(EOS_type), pointer :: EOS + type(EOS_type), pointer :: EOS !> The change in the geopotential anomaly across the layer, in m2 s-2. - real, dimension(SZDI_(HI),SZDJ_(HI)), intent(out) :: dza + real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), intent(out) :: dza !> The integral in pressure through the layer of the geopotential anomaly !! relative to the anomaly at the bottom of the layer, in Pa m2 s-2. - real, dimension(SZDI_(HI),SZDJ_(HI)), optional, intent(out) :: intp_dza + real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), optional, intent(out) :: intp_dza !> The integral in x of the difference between the geopotential anomaly at the !! top and bottom of the layer divided by the x grid spacing, in m2 s-2. - real, dimension(SZDIB_(HI),SZDJ_(HI)), optional, intent(out) :: intx_dza + real, dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), optional, intent(out) :: intx_dza !> The integral in y of the difference between the geopotential anomaly at the !! top and bottom of the layer divided by the y grid spacing, in m2 s-2. - real, dimension(SZDI_(HI),SZDJB_(HI)), optional, intent(out) :: inty_dza + real, dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), optional, intent(out) :: inty_dza !> The width of halo points on which to calculate dza. - integer, optional, intent(in) :: halo_size + integer, optional, intent(in) :: halo_size if (.not.associated(EOS)) call MOM_error(FATAL, & "int_specific_vol_dp called with an unassociated EOS_type EOS.") @@ -374,13 +374,13 @@ subroutine int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, EOS, & !> Ocean horizontal index structures for the output arrays type(hor_index_type), intent(in) :: HIO !> Potential temperature referenced to the surface (degC) - real, dimension(SZDI_(HII),SZDJ_(HII)), intent(in) :: T + real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), intent(in) :: T !> Salinity (PSU) - real, dimension(SZDI_(HII),SZDJ_(HII)), intent(in) :: S + real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), intent(in) :: S !> Height at the top of the layer in m. - real, dimension(SZDI_(HII),SZDJ_(HII)), intent(in) :: z_t + real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), intent(in) :: z_t !> Height at the bottom of the layer in m. - real, dimension(SZDI_(HII),SZDJ_(HII)), intent(in) :: z_b + real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), intent(in) :: z_b !> A mean density, in kg m-3, that is subtracted out to reduce the magnitude !! of each of the integrals. (The pressure is calculated as p~=-z*rho_0*G_e.) real, intent(in) :: rho_ref @@ -392,16 +392,16 @@ subroutine int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, EOS, & !> Equation of state structure type(EOS_type), pointer :: EOS !> The change in the pressure anomaly across the layer, in Pa. - real, dimension(SZDI_(HIO),SZDJ_(HIO)), intent(out) :: dpa + real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), intent(out) :: dpa !> The integral through the thickness of the layer of the pressure anomaly !! relative to the anomaly at the top of the layer, in Pa m. - real, dimension(SZDI_(HIO),SZDJ_(HIO)), optional, intent(out) :: intz_dpa + real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), optional, intent(out) :: intz_dpa !> The integral in x of the difference between the pressure anomaly at the !! top and bottom of the layer divided by the x grid spacing, in Pa. - real, dimension(SZDIB_(HIO),SZDJ_(HIO)), optional, intent(out) :: intx_dpa + real, dimension(HIO%IsdB:HIO%IedB,HIO%jsd:HIO%jed), optional, intent(out) :: intx_dpa !> The integral in y of the difference between the pressure anomaly at the !! top and bottom of the layer divided by the y grid spacing, in Pa. - real, dimension(SZDI_(HIO),SZDJB_(HIO)), optional, intent(out) :: inty_dpa + real, dimension(HIO%isd:HIO%ied,HIO%JsdB:HIO%JedB), optional, intent(out) :: inty_dpa if (.not.associated(EOS)) call MOM_error(FATAL, & "int_density_dz called with an unassociated EOS_type EOS.") @@ -446,7 +446,7 @@ subroutine EOS_init(param_file, EOS) if (.not.associated(EOS)) call EOS_allocate(EOS) ! Read all relevant parameters and write them to the model log. - call log_version(param_file, mod, version) + call log_version(param_file, mod, version, "") call get_param(param_file, mod, "EQN_OF_STATE", tmpstr, & "EQN_OF_STATE determines which ocean equation of state \n"//& @@ -561,14 +561,19 @@ end subroutine EOS_use_linear subroutine int_density_dz_generic(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, & EOS, dpa, intz_dpa, intx_dpa, inty_dpa) - type(hor_index_type), intent(in) :: HII, HIO - real, dimension(SZDI_(HII),SZDJ_(HII)), intent(in) :: T, S, z_t, z_b - real, intent(in) :: rho_ref, rho_0, G_e - type(EOS_type), pointer :: EOS !< Equation of state structure - real, dimension(SZDI_(HIO),SZDJ_(HIO)), intent(out) :: dpa - real, dimension(SZDI_(HIO),SZDJ_(HIO)), optional, intent(out) :: intz_dpa - real, dimension(SZDIB_(HIO),SZDJ_(HIO)), optional, intent(out) :: intx_dpa - real, dimension(SZDI_(HIO),SZDJB_(HIO)), optional, intent(out) :: inty_dpa + type(hor_index_type), intent(in) :: HII, HIO + real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), & + intent(in) :: T, S, z_t, z_b + real, intent(in) :: rho_ref, rho_0, G_e + type(EOS_type), pointer :: EOS !< Equation of state structure + real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), & + intent(out) :: dpa + real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), & + optional, intent(out) :: intz_dpa + real, dimension(HIO%IsdB:HIO%IedB,HIO%jsd:HIO%jed), & + optional, intent(out) :: intx_dpa + real, dimension(HIO%isd:HIO%ied,HIO%JsdB:HIO%JedB), & + optional, intent(out) :: inty_dpa ! This subroutine calculates (by numerical quadrature) integrals of ! pressure anomalies across layers, which are required for calculating the ! finite-volume form pressure accelerations in a Boussinesq model. The one @@ -897,16 +902,21 @@ subroutine int_density_dz_generic_plm (T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, & rho_0, G_e, H_subroundoff, bathyT, HII, HIO, EOS, dpa, & intz_dpa, intx_dpa, inty_dpa, & useMassWghtInterp) - type(hor_index_type), intent(in) :: HII, HIO - real, dimension(SZDI_(HII),SZDJ_(HII)), intent(in) :: T_t, T_b, S_t, S_b, z_t, z_b - real, intent(in) :: rho_ref, rho_0, G_e, H_subroundoff - real, dimension(SZDI_(HII),SZDJ_(HII)), intent(in) :: bathyT - type(EOS_type), pointer :: EOS !< Equation of state structure - real, dimension(SZDI_(HIO),SZDJ_(HIO)), intent(out) :: dpa - real, dimension(SZDI_(HIO),SZDJ_(HIO)), optional, intent(out) :: intz_dpa - real, dimension(SZDIB_(HIO),SZDJ_(HIO)), optional, intent(out) :: intx_dpa - real, dimension(SZDI_(HIO),SZDJB_(HIO)), optional, intent(out) :: inty_dpa - logical, optional, intent(in) :: useMassWghtInterp + type(hor_index_type), intent(in) :: HII, HIO + real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), & + intent(in) :: T_t, T_b, S_t, S_b, z_t, z_b + real, intent(in) :: rho_ref, rho_0, G_e, H_subroundoff + real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), intent(in) :: bathyT + type(EOS_type), pointer :: EOS !< Equation of state structure + real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), & + intent(out) :: dpa + real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), & + optional, intent(out) :: intz_dpa + real, dimension(HIO%IsdB:HIO%IedB,HIO%jsd:HIO%jed), & + optional, intent(out) :: intx_dpa + real, dimension(HIO%isd:HIO%ied,HIO%JsdB:HIO%JedB), & + optional, intent(out) :: inty_dpa + logical, optional, intent(in) :: useMassWghtInterp ! This subroutine calculates (by numerical quadrature) integrals of ! pressure anomalies across layers, which are required for calculating the ! finite-volume form pressure accelerations in a Boussinesq model. The one @@ -1294,16 +1304,19 @@ subroutine int_density_dz_generic_ppm (T, T_t, T_b, S, S_t, S_b, & z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, & EOS, dpa, intz_dpa, intx_dpa, inty_dpa) - type(hor_index_type), intent(in) :: HII, HIO - real, dimension(SZDI_(HII),SZDJ_(HII)), intent(in) :: T, T_t, T_b, S, S_t, S_b, & - z_t, z_b + type(hor_index_type), intent(in) :: HII, HIO + real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), & + intent(in) :: T, T_t, T_b, S, S_t, S_b, z_t, z_b real, intent(in) :: rho_ref, rho_0, G_e - type(EOS_type), pointer :: EOS !< Equation of state structure - real, dimension(SZDI_(HIO),SZDJ_(HIO)), intent(out) :: dpa - real, dimension(SZDI_(HIO),SZDJ_(HIO)), optional, intent(out) :: intz_dpa - - real, dimension(SZDIB_(HIO),SZDJB_(HIO)), optional, intent(out) :: intx_dpa - real, dimension(SZDI_(HIO),SZDJB_(HIO)), optional, intent(out) :: inty_dpa + type(EOS_type), pointer :: EOS !< Equation of state structure + real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), & + intent(out) :: dpa + real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), & + optional, intent(out) :: intz_dpa + real, dimension(HIO%IsdB:HIO%IedB,HIO%JsdB:HIO%JedB), & + optional, intent(out) :: intx_dpa + real, dimension(HIO%isd:HIO%ied,HIO%JsdB:HIO%JedB), & + optional, intent(out) :: inty_dpa ! This subroutine calculates (by numerical quadrature) integrals of ! pressure anomalies across layers, which are required for calculating the ! finite-volume form pressure accelerations in a Boussinesq model. The one @@ -1543,14 +1556,19 @@ end subroutine int_density_dz_generic_ppm subroutine int_density_dz_generic_plm_analytic (T_t, T_b, S_t, S_b, z_t, & z_b, rho_ref, rho_0, G_e, HI, EOS, dpa, intz_dpa, intx_dpa, inty_dpa) - type(hor_index_type), intent(in) :: HI - real, dimension(SZDI_(HI),SZDJ_(HI)), intent(in) :: T_t, T_b, S_t, S_b, z_t, z_b - real, intent(in) :: rho_ref, rho_0, G_e - type(EOS_type), pointer :: EOS !< Equation of state structure - real, dimension(SZDI_(HI),SZDJ_(HI)), intent(out) :: dpa - real, dimension(SZDI_(HI),SZDJ_(HI)), optional, intent(out) :: intz_dpa - real, dimension(SZDIB_(HI),SZDJ_(HI)), optional, intent(out) :: intx_dpa - real, dimension(SZDI_(HI),SZDJB_(HI)), optional, intent(out) :: inty_dpa + type(hor_index_type), intent(in) :: HI + real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & + intent(in) :: T_t, T_b, S_t, S_b, z_t, z_b + real, intent(in) :: rho_ref, rho_0, G_e + type(EOS_type), pointer :: EOS !< Equation of state structure + real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & + intent(out) :: dpa + real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & + optional, intent(out) :: intz_dpa + real, dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), & + optional, intent(out) :: intx_dpa + real, dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), & + optional, intent(out) :: inty_dpa ! This subroutine calculates (by numerical quadrature) integrals of ! pressure anomalies across layers, which are required for calculating the ! finite-volume form pressure accelerations in a Boussinesq model. The one @@ -1949,15 +1967,20 @@ end subroutine evaluate_shape_quadratic subroutine int_spec_vol_dp_generic(T, S, p_t, p_b, alpha_ref, HI, EOS, & dza, intp_dza, intx_dza, inty_dza, halo_size) - type(hor_index_type), intent(in) :: HI - real, dimension(SZDI_(HI),SZDJ_(HI)), intent(in) :: T, S, p_t, p_b - real, intent(in) :: alpha_ref - type(EOS_type), pointer :: EOS !< Equation of state structure - real, dimension(SZDI_(HI),SZDJ_(HI)), intent(out) :: dza - real, dimension(SZDI_(HI),SZDJ_(HI)), optional, intent(out) :: intp_dza - real, dimension(SZDIB_(HI),SZDJ_(HI)), optional, intent(out) :: intx_dza - real, dimension(SZDI_(HI),SZDJB_(HI)), optional, intent(out) :: inty_dza - integer, optional, intent(in) :: halo_size + type(hor_index_type), intent(in) :: HI + real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & + intent(in) :: T, S, p_t, p_b + real, intent(in) :: alpha_ref + type(EOS_type), pointer :: EOS !< Equation of state structure + real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & + intent(out) :: dza + real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & + optional, intent(out) :: intp_dza + real, dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), & + optional, intent(out) :: intx_dza + real, dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), & + optional, intent(out) :: inty_dza + integer, optional, intent(in) :: halo_size ! This subroutine calculates analytical and nearly-analytical integrals in ! pressure across layers of geopotential anomalies, which are required for ! calculating the finite-volume form pressure accelerations in a non-Boussinesq diff --git a/src/equation_of_state/MOM_EOS_Wright.F90 b/src/equation_of_state/MOM_EOS_Wright.F90 index 97e0b90e8b..778dff90ae 100644 --- a/src/equation_of_state/MOM_EOS_Wright.F90 +++ b/src/equation_of_state/MOM_EOS_Wright.F90 @@ -259,13 +259,18 @@ end subroutine calculate_2_densities_wright subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, & dpa, intz_dpa, intx_dpa, inty_dpa) - type(hor_index_type), intent(in) :: HII, HIO - real, dimension(SZDI_(HII),SZDJ_(HII)), intent(in) :: T, S, z_t, z_b - real, intent(in) :: rho_ref, rho_0, G_e - real, dimension(SZDI_(HIO),SZDJ_(HIO)), intent(out) :: dpa - real, dimension(SZDI_(HIO),SZDJ_(HIO)), optional, intent(out) :: intz_dpa - real, dimension(SZDIB_(HIO),SZDJ_(HIO)), optional, intent(out) :: intx_dpa - real, dimension(SZDI_(HIO),SZDJB_(HIO)), optional, intent(out) :: inty_dpa + type(hor_index_type), intent(in) :: HII, HIO + real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), & + intent(in) :: T, S, z_t, z_b + real, intent(in) :: rho_ref, rho_0, G_e + real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), & + intent(out) :: dpa + real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), & + optional, intent(out) :: intz_dpa + real, dimension(HIO%IsdB:HIO%IedB,HIO%jsd:HIO%jed), & + optional, intent(out) :: intx_dpa + real, dimension(HIO%isd:HIO%ied,HIO%JsdB:HIO%JedB), & + optional, intent(out) :: inty_dpa ! This subroutine calculates analytical and nearly-analytical integrals of ! pressure anomalies across layers, which are required for calculating the ! finite-volume form pressure accelerations in a Boussinesq model. @@ -293,7 +298,7 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, ! pressure anomaly at the top and bottom of the layer ! divided by the y grid spacing, in Pa. - real, dimension(SZDI_(HII),SZDJ_(HII)) :: al0_2d, p0_2d, lambda_2d + real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed) :: al0_2d, p0_2d, lambda_2d real :: al0, p0, lambda real :: eps, eps2, rho_anom, rem real :: w_left, w_right, intz(5) @@ -393,14 +398,19 @@ end subroutine int_density_dz_wright subroutine int_spec_vol_dp_wright(T, S, p_t, p_b, alpha_ref, HI, dza, & intp_dza, intx_dza, inty_dza, halo_size) - type(hor_index_type), intent(in) :: HI - real, dimension(SZDI_(HI),SZDJ_(HI)), intent(in) :: T, S, p_t, p_b - real, intent(in) :: alpha_ref - real, dimension(SZDI_(HI),SZDJ_(HI)), intent(out) :: dza - real, dimension(SZDI_(HI),SZDJ_(HI)), optional, intent(out) :: intp_dza - real, dimension(SZDIB_(HI),SZDJ_(HI)), optional, intent(out) :: intx_dza - real, dimension(SZDI_(HI),SZDJB_(HI)), optional, intent(out) :: inty_dza - integer, optional, intent(in) :: halo_size + type(hor_index_type), intent(in) :: HI + real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & + intent(in) :: T, S, p_t, p_b + real, intent(in) :: alpha_ref + real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & + intent(out) :: dza + real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & + optional, intent(out) :: intp_dza + real, dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), & + optional, intent(out) :: intx_dza + real, dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), & + optional, intent(out) :: inty_dza + integer, optional, intent(in) :: halo_size ! This subroutine calculates analytical and nearly-analytical integrals in ! pressure across layers of geopotential anomalies, which are required for ! calculating the finite-volume form pressure accelerations in a non-Boussinesq @@ -431,7 +441,7 @@ subroutine int_spec_vol_dp_wright(T, S, p_t, p_b, alpha_ref, HI, dza, & ! divided by the y grid spacing, in m2 s-2. ! (in,opt) halo_size - The width of halo points on which to calculate dza. - real, dimension(SZDI_(HI),SZDJ_(HI)) :: al0_2d, p0_2d, lambda_2d + real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed) :: al0_2d, p0_2d, lambda_2d real :: al0, p0, lambda real :: alpha_anom, dp, p_ave real :: rem, eps, eps2 diff --git a/src/equation_of_state/MOM_EOS_linear.F90 b/src/equation_of_state/MOM_EOS_linear.F90 index 090e12f5ac..96e2468ae5 100644 --- a/src/equation_of_state/MOM_EOS_linear.F90 +++ b/src/equation_of_state/MOM_EOS_linear.F90 @@ -213,14 +213,19 @@ end subroutine calculate_2_densities_linear subroutine int_density_dz_linear(T, S, z_t, z_b, rho_ref, rho_0_pres, G_e, HII, HIO, & Rho_T0_S0, dRho_dT, dRho_dS, dpa, intz_dpa, intx_dpa, inty_dpa) - type(hor_index_type), intent(in) :: HII, HIO - real, dimension(SZDI_(HII),SZDJ_(HII)), intent(in) :: T, S, z_t, z_b - real, intent(in) :: rho_ref, rho_0_pres, G_e - real, intent(in) :: Rho_T0_S0, dRho_dT, dRho_dS - real, dimension(SZDI_(HIO),SZDJ_(HIO)), intent(out) :: dpa - real, dimension(SZDI_(HIO),SZDJ_(HIO)), optional, intent(out) :: intz_dpa - real, dimension(SZDIB_(HIO),SZDJ_(HIO)), optional, intent(out) :: intx_dpa - real, dimension(SZDI_(HIO),SZDJB_(HIO)), optional, intent(out) :: inty_dpa + type(hor_index_type), intent(in) :: HII, HIO + real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), & + intent(in) :: T, S, z_t, z_b + real, intent(in) :: rho_ref, rho_0_pres, G_e + real, intent(in) :: Rho_T0_S0, dRho_dT, dRho_dS + real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), & + intent(out) :: dpa + real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), & + optional, intent(out) :: intz_dpa + real, dimension(HIO%IsdB:HIO%IedB,HIO%jsd:HIO%jed), & + optional, intent(out) :: intx_dpa + real, dimension(HIO%isd:HIO%ied,HIO%JsdB:HIO%JedB), & + optional, intent(out) :: inty_dpa ! This subroutine calculates analytical and nearly-analytical integrals of ! pressure anomalies across layers, which are required for calculating the ! finite-volume form pressure accelerations in a Boussinesq model. @@ -292,14 +297,19 @@ end subroutine int_density_dz_linear subroutine int_spec_vol_dp_linear(T, S, p_t, p_b, alpha_ref, HI, Rho_T0_S0, & dRho_dT, dRho_dS, dza, intp_dza, intx_dza, inty_dza, halo_size) - type(hor_index_type), intent(in) :: HI - real, dimension(SZDI_(HI),SZDJ_(HI)), intent(in) :: T, S, p_t, p_b - real, intent(in) :: alpha_ref - real, intent(in) :: Rho_T0_S0, dRho_dT, dRho_dS - real, dimension(SZDI_(HI),SZDJ_(HI)), intent(out) :: dza - real, dimension(SZDI_(HI),SZDJ_(HI)), optional, intent(out) :: intp_dza - real, dimension(SZDIB_(HI),SZDJ_(HI)), optional, intent(out) :: intx_dza - real, dimension(SZDI_(HI),SZDJB_(HI)), optional, intent(out) :: inty_dza + type(hor_index_type), intent(in) :: HI + real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & + intent(in) :: T, S, p_t, p_b + real, intent(in) :: alpha_ref + real, intent(in) :: Rho_T0_S0, dRho_dT, dRho_dS + real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & + intent(out) :: dza + real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & + optional, intent(out) :: intp_dza + real, dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), & + optional, intent(out) :: intx_dza + real, dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), & + optional, intent(out) :: inty_dza integer, optional, intent(in) :: halo_size ! This subroutine calculates analytical and nearly-analytical integrals in ! pressure across layers of geopotential anomalies, which are required for diff --git a/src/framework/MOM_checksums.F90 b/src/framework/MOM_checksums.F90 index c30f4a38d2..30a0c1398d 100644 --- a/src/framework/MOM_checksums.F90 +++ b/src/framework/MOM_checksums.F90 @@ -1156,10 +1156,11 @@ function totalStuff(G, hThick, stuff) real, dimension(G%isd:,G%jsd:,:), intent(in) :: hThick !< The array of thicknesses to use as weights real, dimension(G%isd:,G%jsd:,:), intent(in) :: stuff !< The array of stuff to be summed real :: totalStuff - integer :: i, j, k + integer :: i, j, k, nz + nz = size(hThick,3) totalStuff = 0. - do k = 1, G%ke ; do j = G%jsc, G%jec ; do i = G%isc, G%iec + do k = 1, nz ; do j = G%jsc, G%jec ; do i = G%isc, G%iec totalStuff = totalStuff + hThick(i,j,k) * stuff(i,j,k) * G%areaT(i,j) enddo ; enddo ; enddo call sum_across_PEs(totalStuff) @@ -1185,10 +1186,11 @@ subroutine totalTandS(G, hThick, temperature, salinity, mesg) logical, save :: firstCall = .true. real :: thisH, thisT, thisS, delH, delT, delS - integer :: i, j, k + integer :: i, j, k, nz + nz = size(hThick,3) thisH = 0. - do k = 1, G%ke ; do j = G%jsc, G%jec ; do i = G%isc, G%iec + do k = 1, nz ; do j = G%jsc, G%jec ; do i = G%isc, G%iec thisH = thisH + hThick(i,j,k) * G%areaT(i,j) enddo ; enddo ; enddo call sum_across_PEs(thisH) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 233877f152..8cc208a232 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -204,7 +204,7 @@ subroutine set_axes_info(G, GV, param_file, diag_cs, set_vertical) integer :: id_xq, id_yq, id_zl, id_zi, id_xh, id_yh, id_zzl, id_zzi integer :: k, nz integer :: nzi(4) - real :: zlev(G%ke), zinter(G%ke+1) + real :: zlev(GV%ke), zinter(GV%ke+1) logical :: set_vert character(len=200) :: inputdir, string, filename, varname, dimname @@ -215,7 +215,7 @@ subroutine set_axes_info(G, GV, param_file, diag_cs, set_vertical) set_vert = .true. ; if (present(set_vertical)) set_vert = set_vertical ! Read all relevant parameters and write them to the model log. - call log_version(param_file, mod, version) + call log_version(param_file, mod, version, "") if(G%symmetric) then id_xq = diag_axis_init('xq', G%gridLonB(G%isgB:G%iegB), G%x_axis_units, 'x', & @@ -234,7 +234,7 @@ subroutine set_axes_info(G, GV, param_file, diag_cs, set_vertical) 'h point nominal latitude', Domain2=G%Domain%mpp_domain) if (set_vert) then - nz = G%ke + nz = GV%ke zinter(1:nz+1) = GV%sInterface(1:nz+1) zlev(1:nz) = GV%sLayer(1:nz) id_zl = diag_axis_init('zl', zlev, trim(GV%zAxisUnits), 'z', & @@ -262,7 +262,7 @@ subroutine set_axes_info(G, GV, param_file, diag_cs, set_vertical) if (len_trim(string) > 0) then if (trim(string) == 'UNIFORM') then ! initialise a uniform coordinate with depth - nzi(1) = G%ke + 1 + nzi(1) = GV%ke + 1 allocate(diag_cs%zi_remap(nzi(1))) allocate(diag_cs%zl_remap(nzi(1) - 1)) @@ -828,14 +828,14 @@ end subroutine remap_diag_to_z !! height changes. subroutine diag_update_target_grids(diag_cs) type(diag_ctrl), intent(inout) :: diag_cs !< Diagnostics control structure + ! Local variables - real, dimension(size(diag_cs%h, 3)) :: h_src type(ocean_grid_type), pointer :: G real :: depth integer :: nz_src, nz_dest integer :: i, j, k logical :: force, h_changed - real, dimension(diag_cs%G%ke) :: h_src_1d + real, dimension(size(diag_cs%h, 3)) :: h_src_1d real, dimension(diag_cs%nz_remap) :: h_zout_1d real, dimension(diag_cs%nz_remap+1) :: z_out_1d @@ -1654,11 +1654,16 @@ subroutine diag_mediator_infrastructure_init(err_msg) call diag_manager_init(err_msg=err_msg) end subroutine diag_mediator_infrastructure_init -subroutine diag_mediator_init(G, param_file, diag_cs, doc_file_dir) - type(ocean_grid_type), target, intent(inout) :: G - type(param_file_type), intent(in) :: param_file - type(diag_ctrl), intent(inout) :: diag_cs - character(len=*), optional, intent(in) :: doc_file_dir +!> diag_mediator_init initializes the MOM diag_mediator and opens the available +!! diagnostics file, if appropriate. +subroutine diag_mediator_init(G, nz, param_file, diag_cs, doc_file_dir) + type(ocean_grid_type), target, intent(inout) :: G !< The ocean grid type. + integer, intent(in) :: nz !< The number of layers in the model's native grid. + type(param_file_type), intent(in) :: param_file !< Parameter file structure + type(diag_ctrl), intent(inout) :: diag_cs !< A pointer to a type with many variables + !! used for diagnostics + character(len=*), optional, intent(in) :: doc_file_dir !< A directory in which to create the + !! file ! This subroutine initializes the diag_mediator and the diag_manager. ! The grid type should have its dimensions set by this point, but it @@ -1693,7 +1698,7 @@ subroutine diag_mediator_init(G, param_file, diag_cs, doc_file_dir) diag_cs%do_z_remapping_on_T = .false. diag_cs%remapping_initialized = .false. #if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__) - allocate(diag_cs%h_old(G%isd:G%ied,G%jsd:G%jed,G%ke)) + allocate(diag_cs%h_old(G%isd:G%ied,G%jsd:G%jed,nz)) diag_cs%h_old(:,:,:) = 0.0 #endif @@ -1755,11 +1760,13 @@ subroutine diag_set_thickness_ptr(h, diag_cs) end subroutine -subroutine diag_masks_set(G, missing_value, diag_cs) -! Setup the 2d masks for diagnostics - type(ocean_grid_type), target, intent(in) :: G - real, intent(in) :: missing_value - type(diag_ctrl), pointer :: diag_cs +!> diag_masks_set sets up the 2d and 3d masks for diagnostics +subroutine diag_masks_set(G, nz, missing_value, diag_cs) + type(ocean_grid_type), target, intent(in) :: G !< The ocean grid type. + integer, intent(in) :: nz !< The number of layers in the model's native grid. + real, intent(in) :: missing_value !< A value to use for masked points. + type(diag_ctrl), pointer :: diag_cs !< A pointer to a type with many variables + !! used for diagnostics ! Local variables integer :: k @@ -1767,21 +1774,21 @@ subroutine diag_masks_set(G, missing_value, diag_cs) diag_cs%mask2dBu=> G%mask2dBu diag_cs%mask2dCu=> G%mask2dCu diag_cs%mask2dCv=> G%mask2dCv - allocate(diag_cs%mask3dTL(G%isd:G%ied,G%jsd:G%jed,1:G%ke)) - allocate(diag_cs%mask3dBuL(G%IsdB:G%IedB,G%JsdB:G%JedB,1:G%ke)) - allocate(diag_cs%mask3dCuL(G%IsdB:G%IedB,G%jsd:G%jed,1:G%ke)) - allocate(diag_cs%mask3dCvL(G%isd:G%ied,G%JsdB:G%JedB,1:G%ke)) - do k = 1,G%ke + allocate(diag_cs%mask3dTL(G%isd:G%ied,G%jsd:G%jed,1:nz)) + allocate(diag_cs%mask3dBuL(G%IsdB:G%IedB,G%JsdB:G%JedB,1:nz)) + allocate(diag_cs%mask3dCuL(G%IsdB:G%IedB,G%jsd:G%jed,1:nz)) + allocate(diag_cs%mask3dCvL(G%isd:G%ied,G%JsdB:G%JedB,1:nz)) + do k=1,nz diag_cs%mask3dTL(:,:,k) = diag_cs%mask2dT (:,:) diag_cs%mask3dBuL(:,:,k) = diag_cs%mask2dBu(:,:) diag_cs%mask3dCuL(:,:,k) = diag_cs%mask2dCu(:,:) diag_cs%mask3dCvL(:,:,k) = diag_cs%mask2dCv(:,:) enddo - allocate(diag_cs%mask3dTi(G%isd:G%ied,G%jsd:G%jed,1:G%ke+1)) - allocate(diag_cs%mask3dBui(G%IsdB:G%IedB,G%JsdB:G%JedB,1:G%ke+1)) - allocate(diag_cs%mask3dCui(G%IsdB:G%IedB,G%jsd:G%jed,1:G%ke+1)) - allocate(diag_cs%mask3dCvi(G%isd:G%ied,G%JsdB:G%JedB,1:G%ke+1)) - do k = 1,G%ke+1 + allocate(diag_cs%mask3dTi(G%isd:G%ied,G%jsd:G%jed,1:nz+1)) + allocate(diag_cs%mask3dBui(G%IsdB:G%IedB,G%JsdB:G%JedB,1:nz+1)) + allocate(diag_cs%mask3dCui(G%IsdB:G%IedB,G%jsd:G%jed,1:nz+1)) + allocate(diag_cs%mask3dCvi(G%isd:G%ied,G%JsdB:G%JedB,1:nz+1)) + do k=1,nz+1 diag_cs%mask3dTi(:,:,k) = diag_cs%mask2dT (:,:) diag_cs%mask3dBui(:,:,k) = diag_cs%mask2dBu(:,:) diag_cs%mask3dCui(:,:,k) = diag_cs%mask2dCu(:,:) diff --git a/src/framework/MOM_domains.F90 b/src/framework/MOM_domains.F90 index cc44aadf11..3689f7511d 100644 --- a/src/framework/MOM_domains.F90 +++ b/src/framework/MOM_domains.F90 @@ -941,7 +941,7 @@ subroutine MOM_domains_init(MOM_dom, param_file, symmetric, static_memory, & endif ! Read all relevant parameters and write them to the model log. - call log_version(param_file, mod, version) + call log_version(param_file, mod, version, "") call get_param(param_file, mod, "REENTRANT_X", reentrant_x, & "If true, the domain is zonally reentrant.", default=.true.) call get_param(param_file, mod, "REENTRANT_Y", reentrant_y, & diff --git a/src/framework/MOM_file_parser.F90 b/src/framework/MOM_file_parser.F90 index 397084cbc5..aceace13ae 100644 --- a/src/framework/MOM_file_parser.F90 +++ b/src/framework/MOM_file_parser.F90 @@ -284,7 +284,7 @@ subroutine close_param_file(CS, quiet_close, component) ! Log the parameters for the parser. mod = "MOM_file_parser" - call log_version(CS, mod, version) + call log_version(CS, mod, version, "") call log_param(CS, mod, "SEND_LOG_TO_STDOUT", & CS%log_to_stdout, & "If true, all log messages are also sent to stdout.", & diff --git a/src/framework/MOM_io.F90 b/src/framework/MOM_io.F90 index 7a1741d7b3..91cf9137bf 100644 --- a/src/framework/MOM_io.F90 +++ b/src/framework/MOM_io.F90 @@ -7,8 +7,9 @@ module MOM_io use MOM_error_handler, only : MOM_error, NOTE, FATAL, WARNING use MOM_domains, only : MOM_domain_type use MOM_file_parser, only : log_version, param_file_type -use MOM_string_functions, only : lowercase use MOM_grid, only : ocean_grid_type +use MOM_dyn_horgrid, only : dyn_horgrid_type +use MOM_string_functions, only : lowercase use MOM_verticalGrid, only : verticalGrid_type use ensemble_manager_mod, only : get_ensemble_id @@ -76,7 +77,7 @@ module MOM_io !> Routine creates a new NetCDF file. It also sets up !! structures that describe this file and variables that will !! later be written to this file. Type for describing a variable, typically a tracer -subroutine create_file(unit, filename, vars, novars, fields, threading, timeunit, G, GV) +subroutine create_file(unit, filename, vars, novars, fields, threading, timeunit, G, dG, GV) integer, intent(out) :: unit !< unit id of an open file or -1 on a !! nonwriting PE with single file output character(len=*), intent(in) :: filename !< full path to the file to create @@ -86,8 +87,11 @@ subroutine create_file(unit, filename, vars, novars, fields, threading, timeunit integer, optional, intent(in) :: threading !< SINGLE_FILE or MULTIPLE real, optional, intent(in) :: timeunit !< length, in seconds, of the units for time. The !! default value is 86400.0, for 1 day. - type(ocean_grid_type), optional, intent(in) :: G !< ocean grid structure structure, which is - !! required if the new file uses any + type(ocean_grid_type), optional, intent(in) :: G !< ocean horizontal grid structure; G or dG + !! is required if the new file uses any + !! horizontal grid axes. + type(dyn_horgrid_type), optional, intent(in) :: dG !< dynamic horizontal grid structure; G or dG + !! is required if the new file uses any !! horizontal grid axes. type(verticalGrid_type), optional, intent(in) :: GV !< ocean vertical grid structure, which is !! required if the new file uses any @@ -95,15 +99,21 @@ subroutine create_file(unit, filename, vars, novars, fields, threading, timeunit logical :: use_lath, use_lonh, use_latq, use_lonq, use_time logical :: use_layer, use_int, use_periodic - logical :: one_file + logical :: one_file, domain_set type(axistype) :: axis_lath, axis_latq, axis_lonh, axis_lonq type(axistype) :: axis_layer, axis_int, axis_time, axis_periodic type(axistype) :: axes(4) + type(MOM_domain_type), pointer :: Domain => NULL() type(domain1d) :: x_domain, y_domain integer :: numaxes, pack, thread, k + integer :: isg, ieg, jsg, jeg, IsgB, IegB, JsgB, JegB integer :: var_periods, num_periods=0 real, dimension(:), allocatable :: period_val - character(len=40) :: time_units + real, pointer, dimension(:) :: & + gridLatT => NULL(), & ! The latitude or longitude of T or B points for + gridLatB => NULL(), & ! the purpose of labeling the output axes. + gridLonT => NULL(), gridLonB => NULL() + character(len=40) :: time_units, x_axis_units, y_axis_units character(len=8) :: t_grid, t_grid_read use_lath = .false. ; use_lonh = .false. @@ -114,15 +124,32 @@ subroutine create_file(unit, filename, vars, novars, fields, threading, timeunit thread = SINGLE_FILE if (PRESENT(threading)) thread = threading - one_file = .true. + domain_set = .false. if (present(G)) then - one_file = ((thread == SINGLE_FILE) .or. .not.G%Domain%use_io_layout) + domain_set = .true. ; Domain => G%Domain + gridLatT => G%gridLatT ; gridLatB => G%gridLatB + gridLonT => G%gridLonT ; gridLonB => G%gridLonB + x_axis_units = G%x_axis_units ; y_axis_units = G%y_axis_units + isg = G%isg ; ieg = G%ieg ; jsg = G%jsg ; jeg = G%jeg + IsgB = G%IsgB ; IegB = G%IegB ; JsgB = G%JsgB ; JegB = G%JegB + elseif (present(dG)) then + domain_set = .true. ; Domain => dG%Domain + gridLatT => dG%gridLatT ; gridLatB => dG%gridLatB + gridLonT => dG%gridLonT ; gridLonB => dG%gridLonB + x_axis_units = dG%x_axis_units ; y_axis_units = dG%y_axis_units + isg = dG%isg ; ieg = dG%ieg ; jsg = dG%jsg ; jeg = dG%jeg + IsgB = dG%IsgB ; IegB = dG%IegB ; JsgB = dG%JsgB ; JegB = dG%JegB + endif + + one_file = .true. + if (domain_set) then + one_file = ((thread == SINGLE_FILE) .or. .not.Domain%use_io_layout) endif if (one_file) then call open_file(unit, filename, MPP_OVERWR, MPP_NETCDF, threading=thread) else - call open_file(unit, filename, MPP_OVERWR, MPP_NETCDF, domain=G%Domain%mpp_domain) + call open_file(unit, filename, MPP_OVERWR, MPP_NETCDF, domain=Domain%mpp_domain) endif ! Define the coordinates. @@ -183,10 +210,10 @@ subroutine create_file(unit, filename, vars, novars, fields, threading, timeunit end do if ((use_lath .or. use_lonh .or. use_latq .or. use_lonq)) then - if (.not.present(G)) call MOM_error(FATAL, "create_file: "//& - "An ocean_grid_type is required to create a file with a vertical coordinate.") + if (.not.domain_set) call MOM_error(FATAL, "create_file: "//& + "An ocean_grid_type or dyn_horgrid_type is required to create a file with a horizontal coordinate.") - call mpp_get_domain_components(G%Domain%mpp_domain, x_domain, y_domain) + call mpp_get_domain_components(Domain%mpp_domain, x_domain, y_domain) endif if ((use_layer .or. use_int) .and. .not.present(GV)) call MOM_error(FATAL, & "create_file: A vertical grid type is required to create a file with a vertical coordinate.") @@ -194,20 +221,20 @@ subroutine create_file(unit, filename, vars, novars, fields, threading, timeunit ! Specify all optional arguments to mpp_write_meta: name, units, longname, cartesian, calendar, sense, domain, data, min) ! Otherwise if optional arguments are added to mpp_write_meta the compiler may (and in case of GNU is) get confused and crash. if (use_lath) & - call mpp_write_meta(unit, axis_lath, name="lath", units=G%y_axis_units, longname="Latitude", & - cartesian='Y', domain = y_domain, data=G%gridLatT(G%jsg:G%jeg)) + call mpp_write_meta(unit, axis_lath, name="lath", units=y_axis_units, longname="Latitude", & + cartesian='Y', domain = y_domain, data=gridLatT(jsg:jeg)) if (use_lonh) & - call mpp_write_meta(unit, axis_lonh, name="lonh", units=G%x_axis_units, longname="Longitude", & - cartesian='X', domain = x_domain, data=G%gridLonT(G%isg:G%ieg)) + call mpp_write_meta(unit, axis_lonh, name="lonh", units=x_axis_units, longname="Longitude", & + cartesian='X', domain = x_domain, data=gridLonT(isg:ieg)) if (use_latq) & - call mpp_write_meta(unit, axis_latq, name="latq", units=G%y_axis_units, longname="Latitude", & - cartesian='Y', domain = y_domain, data=G%gridLatB(G%JsgB:G%JegB)) + call mpp_write_meta(unit, axis_latq, name="latq", units=y_axis_units, longname="Latitude", & + cartesian='Y', domain = y_domain, data=gridLatB(JsgB:JegB)) if (use_lonq) & - call mpp_write_meta(unit, axis_lonq, name="lonq", units=G%x_axis_units, longname="Longitude", & - cartesian='X', domain = x_domain, data=G%gridLonB(G%IsgB:G%IegB)) + call mpp_write_meta(unit, axis_lonq, name="lonq", units=x_axis_units, longname="Longitude", & + cartesian='X', domain = x_domain, data=gridLonB(IsgB:IegB)) if (use_layer) & call mpp_write_meta(unit, axis_layer, name="Layer", units=trim(GV%zAxisUnits), & @@ -305,24 +332,30 @@ end subroutine create_file !! does not find the file, a new file is created. It also sets up !! structures that describe this file and the variables that will !! later be written to this file. -subroutine reopen_file(unit, filename, vars, novars, G, fields, threading, timeunit, GV) +subroutine reopen_file(unit, filename, vars, novars, fields, threading, timeunit, G, dG, GV) integer, intent(out) :: unit !< unit id of an open file or -1 on a !! nonwriting PE with single file output character(len=*), intent(in) :: filename !< full path to the file to create type(vardesc), intent(in) :: vars(:) !< structures describing fields written to filename integer, intent(in) :: novars !< number of fields written to filename - type(ocean_grid_type), intent(in) :: G !< ocean grid structure type(fieldtype), intent(inout) :: fields(:) !< array of fieldtypes for each variable integer, optional, intent(in) :: threading !< SINGLE_FILE or MULTIPLE real, optional, intent(in) :: timeunit !< length, in seconds, of the units for time. The !! default value is 86400.0, for 1 day. + type(ocean_grid_type), optional, intent(in) :: G !< ocean horizontal grid structure; G or dG + !! is required if a new file uses any + !! horizontal grid axes. + type(dyn_horgrid_type), optional, intent(in) :: dG !< dynamic horizontal grid structure; G or dG + !! is required if a new file uses any + !! horizontal grid axes. type(verticalGrid_type), optional, intent(in) :: GV !< ocean vertical grid structure, which is - !! required if the new file uses any + !! required if a new file uses any !! vertical grid axes. + type(MOM_domain_type), pointer :: Domain => NULL() character(len=200) :: check_name, mesg integer :: length, ndim, nvar, natt, ntime, thread - logical :: exists + logical :: exists, one_file, domain_set thread = SINGLE_FILE if (PRESENT(threading)) thread = threading @@ -335,12 +368,26 @@ subroutine reopen_file(unit, filename, vars, novars, G, fields, threading, timeu inquire(file=check_name,EXIST=exists) if (.not.exists) then - call create_file(unit, filename, vars, novars, fields, threading, timeunit, G=G, GV=GV) + call create_file(unit, filename, vars, novars, fields, threading, timeunit, & + G=G, dG=dG, GV=GV) else - if ((thread == SINGLE_FILE) .or. .not.G%Domain%use_io_layout) then + + domain_set = .false. + if (present(G)) then + domain_set = .true. ; Domain => G%Domain + elseif (present(dG)) then + domain_set = .true. ; Domain => dG%Domain + endif + + one_file = .true. + if (domain_set) then + one_file = ((thread == SINGLE_FILE) .or. .not.Domain%use_io_layout) + endif + + if (one_file) then call open_file(unit, filename, MPP_APPEND, MPP_NETCDF, threading=thread) else - call open_file(unit, filename, MPP_APPEND, MPP_NETCDF, domain=G%Domain%mpp_domain) + call open_file(unit, filename, MPP_APPEND, MPP_NETCDF, domain=Domain%mpp_domain) endif if (unit < 0) return diff --git a/src/framework/MOM_restart.F90 b/src/framework/MOM_restart.F90 index 57a8be9451..217aeaa0f5 100644 --- a/src/framework/MOM_restart.F90 +++ b/src/framework/MOM_restart.F90 @@ -759,7 +759,7 @@ subroutine save_restart(directory, time, G, CS, time_stamped, filename, GV) num_files = 0 next_var = 0 - nz = G%ke + nz = 1 ; if (present(GV)) nz = GV%ke call get_time(time,seconds,days) restart_time = real(days) + real(seconds)/86400.0 @@ -1237,7 +1237,7 @@ subroutine restart_init(param_file, CS, restart_root) allocate(CS) ! Read all relevant parameters and write them to the model log. - call log_version(param_file, mod, version) + call log_version(param_file, mod, version, "") call get_param(param_file, mod, "PARALLEL_RESTARTFILES", & CS%parallel_restartfiles, & "If true, each processor writes its own restart file, \n"//& diff --git a/src/framework/MOM_spatial_means.F90 b/src/framework/MOM_spatial_means.F90 index 0bce719151..54f58b7269 100644 --- a/src/framework/MOM_spatial_means.F90 +++ b/src/framework/MOM_spatial_means.F90 @@ -27,6 +27,7 @@ module MOM_spatial_means use MOM_error_handler, only : MOM_error, NOTE, WARNING, FATAL, is_root_pe use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_grid, only : ocean_grid_type +use MOM_verticalGrid, only : verticalGrid_type implicit none ; private @@ -74,21 +75,22 @@ function global_area_integral(var,G) end function global_area_integral -function global_layer_mean(var,h,G) - type(ocean_grid_type), intent(in) :: G - real, dimension(SZI_(G), SZJ_(G), SZK_(G)), intent(in) :: var - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h - real, dimension(SZK_(G)) :: global_layer_mean +function global_layer_mean(var, h, G, GV) + type(ocean_grid_type), intent(in) :: G + type(verticalGrid_type), intent(in) :: GV + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: var + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h + real, dimension(SZK_(GV)) :: global_layer_mean - real, dimension(SZI_(G), SZJ_(G), SZK_(G)) :: tmpForSumming, weight - real, dimension(SZK_(G)) :: scalarij, weightij - real, dimension(SZK_(G)) :: global_temp_scalar, global_weight_scalar + real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: tmpForSumming, weight + real, dimension(SZK_(GV)) :: scalarij, weightij + real, dimension(SZK_(GV)) :: global_temp_scalar, global_weight_scalar integer :: i, j, k, is, ie, js, je, nz - is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke tmpForSumming(:,:,:) = 0. ; weight(:,:,:) = 0. - do k=1, nz ; do j=js,je ; do i=is, ie + do k=1,nz ; do j=js,je ; do i=is,ie weight(i,j,k) = h(i,j,k) * (G%areaT(i,j) * G%mask2dT(i,j)) tmpForSumming(i,j,k) = var(i,j,k) * weight(i,j,k) enddo ; enddo ; enddo @@ -102,23 +104,27 @@ function global_layer_mean(var,h,G) end function global_layer_mean -function global_volume_mean(var,h,G) - type(ocean_grid_type), intent(in) :: G - real, dimension(SZI_(G), SZJ_(G), SZK_(G)), intent(in) :: var - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h +function global_volume_mean(var, h, G, GV) + type(ocean_grid_type), intent(in) :: G + type(verticalGrid_type), intent(in) :: GV + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: var + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h real :: global_volume_mean - real, dimension(SZI_(G), SZJ_(G), SZK_(G)) :: tmpForSumming, weight + real :: weight_here + real, dimension(SZI_(G), SZJ_(G)) :: tmpForSumming, sum_weight integer :: i, j, k, is, ie, js, je, nz - is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke - tmpForSumming(:,:,:) = 0. ; weight(:,:,:) = 0. + tmpForSumming(:,:) = 0. ; sum_weight(:,:) = 0. - do k=1, nz ; do j=js,je ; do i=is, ie - weight(i,j,k) = h(i,j,k) * (G%areaT(i,j) * G%mask2dT(i,j)) - tmpForSumming(i,j,k) = var(i,j,k) * weight(i,j,k) + do k=1,nz ; do j=js,je ; do i=is,ie + weight_here = h(i,j,k) * (G%areaT(i,j) * G%mask2dT(i,j)) + tmpForSumming(i,j) = tmpForSumming(i,j) + var(i,j,k) * weight_here + sum_weight(i,j) = sum_weight(i,j) + weight_here enddo ; enddo ; enddo - global_volume_mean = (reproducing_sum(sum(tmpForSumming,3))) / (reproducing_sum(sum(weight,3))) + global_volume_mean = (reproducing_sum(tmpForSumming)) / & + (reproducing_sum(sum_weight)) end function global_volume_mean @@ -137,7 +143,7 @@ subroutine global_i_mean(array, i_mean, G, mask) ! (in) G - The ocean's grid structure. ! (in) mask - An array used for weighting the i-mean. - type(EFP_type), allocatable, dimension(:) :: sum, mask_sum + type(EFP_type), allocatable, dimension(:) :: asum, mask_sum real :: mask_sum_r integer :: is, ie, js, je, idg_off, jdg_off integer :: i, j @@ -147,23 +153,23 @@ subroutine global_i_mean(array, i_mean, G, mask) call reset_EFP_overflow_error() - allocate(sum(G%jsg:G%jeg)) + allocate(asum(G%jsg:G%jeg)) if (present(mask)) then allocate(mask_sum(G%jsg:G%jeg)) do j=G%jsg,G%jeg - sum(j) = real_to_EFP(0.0) ; mask_sum(j) = real_to_EFP(0.0) + asum(j) = real_to_EFP(0.0) ; mask_sum(j) = real_to_EFP(0.0) enddo do i=is,ie ; do j=js,je - sum(j+jdg_off) = sum(j+jdg_off) + real_to_EFP(array(i,j)*mask(i,j)) + asum(j+jdg_off) = asum(j+jdg_off) + real_to_EFP(array(i,j)*mask(i,j)) mask_sum(j+jdg_off) = mask_sum(j+jdg_off) + real_to_EFP(mask(i,j)) enddo ; enddo if (query_EFP_overflow_error()) call MOM_error(FATAL, & "global_i_mean overflow error occurred before sums across PEs.") - call EFP_list_sum_across_PEs(sum(G%jsg:G%jeg), G%jeg-G%jsg+1) + call EFP_list_sum_across_PEs(asum(G%jsg:G%jeg), G%jeg-G%jsg+1) call EFP_list_sum_across_PEs(mask_sum(G%jsg:G%jeg), G%jeg-G%jsg+1) if (query_EFP_overflow_error()) call MOM_error(FATAL, & @@ -172,32 +178,32 @@ subroutine global_i_mean(array, i_mean, G, mask) do j=js,je mask_sum_r = EFP_to_real(mask_sum(j+jdg_off)) if (mask_sum_r == 0.0 ) then ; i_mean(j) = 0.0 ; else - i_mean(j) = EFP_to_real(sum(j+jdg_off)) / mask_sum_r + i_mean(j) = EFP_to_real(asum(j+jdg_off)) / mask_sum_r endif enddo deallocate(mask_sum) else - do j=G%jsg,G%jeg ; sum(j) = real_to_EFP(0.0) ; enddo + do j=G%jsg,G%jeg ; asum(j) = real_to_EFP(0.0) ; enddo do i=is,ie ; do j=js,je - sum(j+jdg_off) = sum(j+jdg_off) + real_to_EFP(array(i,j)) + asum(j+jdg_off) = asum(j+jdg_off) + real_to_EFP(array(i,j)) enddo ; enddo if (query_EFP_overflow_error()) call MOM_error(FATAL, & "global_i_mean overflow error occurred before sum across PEs.") - call EFP_list_sum_across_PEs(sum(G%jsg:G%jeg), G%jeg-G%jsg+1) + call EFP_list_sum_across_PEs(asum(G%jsg:G%jeg), G%jeg-G%jsg+1) if (query_EFP_overflow_error()) call MOM_error(FATAL, & "global_i_mean overflow error occurred during sum across PEs.") do j=js,je - i_mean(j) = EFP_to_real(sum(j+jdg_off)) / real(G%ieg-G%isg+1) + i_mean(j) = EFP_to_real(asum(j+jdg_off)) / real(G%ieg-G%isg+1) enddo endif - deallocate(sum) + deallocate(asum) end subroutine global_i_mean @@ -215,7 +221,7 @@ subroutine global_j_mean(array, j_mean, G, mask) ! (in) G - The ocean's grid structure. ! (in) mask - An array used for weighting the j-mean. - type(EFP_type), allocatable, dimension(:) :: sum, mask_sum + type(EFP_type), allocatable, dimension(:) :: asum, mask_sum real :: mask_sum_r integer :: is, ie, js, je, idg_off, jdg_off integer :: i, j @@ -225,23 +231,23 @@ subroutine global_j_mean(array, j_mean, G, mask) call reset_EFP_overflow_error() - allocate(sum(G%isg:G%ieg)) + allocate(asum(G%isg:G%ieg)) if (present(mask)) then allocate (mask_sum(G%isg:G%ieg)) do i=G%isg,G%ieg - sum(i) = real_to_EFP(0.0) ; mask_sum(i) = real_to_EFP(0.0) + asum(i) = real_to_EFP(0.0) ; mask_sum(i) = real_to_EFP(0.0) enddo do i=is,ie ; do j=js,je - sum(i+idg_off) = sum(i+idg_off) + real_to_EFP(array(i,j)*mask(i,j)) + asum(i+idg_off) = asum(i+idg_off) + real_to_EFP(array(i,j)*mask(i,j)) mask_sum(i+idg_off) = mask_sum(i+idg_off) + real_to_EFP(mask(i,j)) enddo ; enddo if (query_EFP_overflow_error()) call MOM_error(FATAL, & "global_j_mean overflow error occurred before sums across PEs.") - call EFP_list_sum_across_PEs(sum(G%isg:G%ieg), G%ieg-G%isg+1) + call EFP_list_sum_across_PEs(asum(G%isg:G%ieg), G%ieg-G%isg+1) call EFP_list_sum_across_PEs(mask_sum(G%isg:G%ieg), G%ieg-G%isg+1) if (query_EFP_overflow_error()) call MOM_error(FATAL, & @@ -250,32 +256,32 @@ subroutine global_j_mean(array, j_mean, G, mask) do i=is,ie mask_sum_r = EFP_to_real(mask_sum(i+idg_off)) if (mask_sum_r == 0.0 ) then ; j_mean(i) = 0.0 ; else - j_mean(i) = EFP_to_real(sum(i+idg_off)) / mask_sum_r + j_mean(i) = EFP_to_real(asum(i+idg_off)) / mask_sum_r endif enddo deallocate(mask_sum) else - do i=G%isg,G%ieg ; sum(i) = real_to_EFP(0.0) ; enddo + do i=G%isg,G%ieg ; asum(i) = real_to_EFP(0.0) ; enddo do i=is,ie ; do j=js,je - sum(i+idg_off) = sum(i+idg_off) + real_to_EFP(array(i,j)) + asum(i+idg_off) = asum(i+idg_off) + real_to_EFP(array(i,j)) enddo ; enddo if (query_EFP_overflow_error()) call MOM_error(FATAL, & "global_j_mean overflow error occurred before sum across PEs.") - call EFP_list_sum_across_PEs(sum(G%isg:G%ieg), G%ieg-G%isg+1) + call EFP_list_sum_across_PEs(asum(G%isg:G%ieg), G%ieg-G%isg+1) if (query_EFP_overflow_error()) call MOM_error(FATAL, & "global_j_mean overflow error occurred during sum across PEs.") do i=is,ie - j_mean(i) = EFP_to_real(sum(i+idg_off)) / real(G%jeg-G%jsg+1) + j_mean(i) = EFP_to_real(asum(i+idg_off)) / real(G%jeg-G%jsg+1) enddo endif - deallocate(sum) + deallocate(asum) end subroutine global_j_mean diff --git a/src/framework/MOM_write_cputime.F90 b/src/framework/MOM_write_cputime.F90 index 570a97ef7d..43fb3519bf 100644 --- a/src/framework/MOM_write_cputime.F90 +++ b/src/framework/MOM_write_cputime.F90 @@ -100,7 +100,7 @@ subroutine MOM_write_cputime_init(param_file, directory, Input_start_time, CS) endif ! Read all relevant parameters and write them to the model log. - call log_version(param_file, mod, version) + call log_version(param_file, mod, version, "") call get_param(param_file, mod, "MAXCPU", CS%maxcpu, & "The maximum amount of cpu time per processor for which \n"//& "MOM should run before saving a restart file and \n"//& diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index ef2b752cd9..14cde0c674 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -99,7 +99,9 @@ module MOM_ice_shelf use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr use MOM_diag_mediator, only : diag_mediator_init, set_diag_mediator_grid use MOM_diag_mediator, only : diag_ctrl, time_type, enable_averaging, disable_averaging -use MOM_domains, only : MOM_domains_init, pass_var, pass_vector, TO_ALL, CGRID_NE, BGRID_NE +use MOM_domains, only : MOM_domains_init, clone_MOM_domain +use MOM_domains, only : pass_var, pass_vector, TO_ALL, CGRID_NE, BGRID_NE +use MOM_dyn_horgrid, only : dyn_horgrid_type, create_dyn_horgrid, destroy_dyn_horgrid use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING, is_root_pe use MOM_file_parser, only : read_param, get_param, log_param, log_version, param_file_type use MOM_grid, only : MOM_grid_init, ocean_grid_type @@ -113,6 +115,7 @@ module MOM_ice_shelf use MOM_restart, only : register_restart_field, query_initialized, save_restart use MOM_restart, only : restart_init, restore_state, MOM_restart_CS use MOM_time_manager, only : time_type, set_time, time_type_to_real +use MOM_transcribe_grid, only : copy_dyngrid_to_MOM_grid, copy_MOM_grid_to_dyngrid !use MOM_variables, only : forcing, surface use MOM_variables, only : surface use MOM_forcing_type, only : forcing, allocate_forcing_type @@ -152,6 +155,7 @@ module MOM_ice_shelf type, public :: ice_shelf_CS ; private type(MOM_restart_CS), pointer :: restart_CSp => NULL() type(ocean_grid_type) :: grid !< Grid for the ice-shelf model +! type(dyn_horgrid_type), pointer :: dG !< Dynamic grid for the ice-shelf model type(ocean_grid_type), pointer :: ocn_grid => NULL() !< A pointer to the ocean model grid ! The rest is private real :: flux_factor = 1.0 @@ -685,7 +689,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS) fluxes%iceshelf_melt = CS%lprec * (86400.0*365.0/CS%density_ice) if (CS%DEBUG) then - call hchksum (CS%h_shelf, "melt rate", G, haloshift=0) + call hchksum (CS%h_shelf, "melt rate", G%HI, haloshift=0) endif if (CS%shelf_mass_is_dynamic) then @@ -830,10 +834,10 @@ subroutine add_shelf_flux(G, CS, state, fluxes) if (CS%debug) then if (associated(state%taux_shelf)) then - call uchksum(state%taux_shelf, "taux_shelf", G, haloshift=0) + call uchksum(state%taux_shelf, "taux_shelf", G%HI, haloshift=0) endif if (associated(state%tauy_shelf)) then - call vchksum(state%tauy_shelf, "tauy_shelf", G, haloshift=0) + call vchksum(state%tauy_shelf, "tauy_shelf", G%HI, haloshift=0) endif endif @@ -896,7 +900,7 @@ subroutine add_shelf_flux(G, CS, state, fluxes) endif enddo ; enddo if (CS%debug) then - call hchksum(fluxes%ustar_shelf, "ustar_shelf", G, haloshift=0) + call hchksum(fluxes%ustar_shelf, "ustar_shelf", G%HI, haloshift=0) endif ! If the shelf mass is changing, the fluxes%rigidity_ice_[uv] needs to be @@ -973,10 +977,10 @@ end subroutine add_shelf_flux ! if (CS%debug) then ! if (associated(state%taux_shelf)) then -! call uchksum(state%taux_shelf, "taux_shelf", G, haloshift=0) +! call uchksum(state%taux_shelf, "taux_shelf", G%HI, haloshift=0) ! endif ! if (associated(state%tauy_shelf)) then -! call vchksum(state%tauy_shelf, "tauy_shelf", G, haloshift=0) +! call vchksum(state%tauy_shelf, "tauy_shelf", G%HI, haloshift=0) ! endif ! endif @@ -1021,7 +1025,7 @@ end subroutine add_shelf_flux ! enddo ; enddo ! if (CS%debug) then -! call hchksum(fluxes%ustar_shelf, "ustar_shelf", G, haloshift=0) +! call hchksum(fluxes%ustar_shelf, "ustar_shelf", G%HI, haloshift=0) ! endif ! ! If the shelf mass is changing, the fluxes%rigidity_ice_[uv] needs to be @@ -1057,6 +1061,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, fluxes, Ti type(ocean_grid_type), pointer :: G, OG ! Convenience pointers type(directories) :: dirs type(vardesc) :: vd + type(dyn_horgrid_type), pointer :: dG => NULL() real :: cdrag, drag_bg_vel logical :: new_sim, save_IC, var_force ! This include declares and sets the variable "version". @@ -1089,7 +1094,11 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, fluxes, Ti call MOM_domains_init(CS%grid%domain, param_file, min_halo=wd_halos, symmetric=GRID_SYM_) ! call diag_mediator_init(CS%grid,param_file,CS%diag) ! this needs to be fixed - will probably break when not using coupled driver 0 call MOM_grid_init(CS%grid, param_file) - call set_grid_metrics(CS%grid, param_file) + + call create_dyn_horgrid(dG, CS%grid%HI) + call clone_MOM_domain(CS%grid%Domain, dG%Domain) + + call set_grid_metrics(dG, param_file) ! call set_diag_mediator_grid(CS%grid, CS%diag) ! The ocean grid is possibly different @@ -1396,9 +1405,12 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, fluxes, Ti endif ! Set up the bottom depth, G%D either analytically or from file - call MOM_initialize_topography(G%bathyT, G%max_depth, G, param_file) + call MOM_initialize_topography(G%bathyT, G%max_depth, dG, param_file) ! Set up the Coriolis parameter, G%f, usually analytically. - call MOM_initialize_rotation(G%CoriolisBu, G, param_file) + call MOM_initialize_rotation(G%CoriolisBu, dG, param_file) + call copy_dyngrid_to_MOM_grid(dG, CS%grid) + + call destroy_dyn_horgrid(dG) ! Set up the restarts. call restart_init(param_file, CS%restart_CSp, "Shelf.res") @@ -2031,8 +2043,8 @@ subroutine ice_shelf_advect(CS, time_step, melt_rate, Time) call pass_var(CS%hmask, G%domain) if (CS%DEBUG) then - call hchksum (CS%h_shelf, "h after front", G, haloshift=3) - call hchksum (CS%h_shelf, "shelf area after front", G, haloshift=3) + call hchksum (CS%h_shelf, "h after front", G%HI, haloshift=3) + call hchksum (CS%h_shelf, "shelf area after front", G%HI, haloshift=3) endif @@ -2257,8 +2269,8 @@ subroutine ice_shelf_solve_outer (CS, u, v, FE, iters, time) if (CS%DEBUG) then - call qchksum (u, "u shelf", G, haloshift=2) - call qchksum (v, "v shelf", G, haloshift=2) + call qchksum (u, "u shelf", G%HI, haloshift=2) + call qchksum (v, "v shelf", G%HI, haloshift=2) endif if (is_root_pe()) print *,"linear solve done",iters," iterations" @@ -5920,7 +5932,7 @@ subroutine ice_shelf_temp(CS, time_step, melt_rate, Time) call pass_var(CS%tmask, G%domain) if (CS%DEBUG) then - call hchksum (CS%t_shelf, "temp after front", G, haloshift=3) + call hchksum (CS%t_shelf, "temp after front", G%HI, haloshift=3) endif end subroutine ice_shelf_temp diff --git a/src/initialization/MOM_coord_initialization.F90 b/src/initialization/MOM_coord_initialization.F90 index c200829b8e..3f341a25b7 100644 --- a/src/initialization/MOM_coord_initialization.F90 +++ b/src/initialization/MOM_coord_initialization.F90 @@ -53,7 +53,7 @@ subroutine MOM_initialize_coord(GV, PF, write_geom, output_dir, tv, max_depth) nz = GV%ke call callTree_enter("MOM_initialize_coord(), MOM_coord_initialization.F90") -! call log_version(PF, mod, version) + call log_version(PF, mod, version, "") call get_param(PF, mod, "DEBUG", debug, default=.false.) ! Set-up the layer densities, GV%Rlay, and reduced gravities, GV%g_prime. diff --git a/src/initialization/MOM_fixed_initialization.F90 b/src/initialization/MOM_fixed_initialization.F90 index ec8f85d278..c7ffc4bd53 100644 --- a/src/initialization/MOM_fixed_initialization.F90 +++ b/src/initialization/MOM_fixed_initialization.F90 @@ -8,11 +8,12 @@ module MOM_fixed_initialization use MOM_coms, only : max_across_PEs use MOM_domains, only : pass_var, pass_vector, sum_across_PEs, broadcast use MOM_domains, only : root_PE, To_All, SCALAR_PAIR, CGRID_NE, AGRID +use MOM_dyn_horgrid, only : dyn_horgrid_type use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, WARNING, is_root_pe use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint use MOM_file_parser, only : get_param, read_param, log_param, param_file_type use MOM_file_parser, only : log_version -use MOM_grid, only : ocean_grid_type, isPointInCell +! use MOM_grid, only : ocean_grid_type use MOM_io, only : close_file, create_file, fieldtype, file_exists use MOM_io, only : open_file, read_data, read_axis_data, SINGLE_FILE, MULTIPLE use MOM_io, only : slasher, vardesc, write_field, var_desc @@ -48,7 +49,7 @@ module MOM_fixed_initialization !> MOM_initialize_fixed sets up time-invariant quantities related to MOM6's !! horizontal grid, bathymetry, and the Coriolis parameter. subroutine MOM_initialize_fixed(G, OBC, PF, write_geom, output_dir) - type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. + type(dyn_horgrid_type), intent(inout) :: G !< The ocean's grid structure. type(ocean_OBC_type), pointer :: OBC !< Open boundary structure. type(param_file_type), intent(in) :: PF !< A structure indicating the open file !! to parse for model parameter values. @@ -63,7 +64,7 @@ subroutine MOM_initialize_fixed(G, OBC, PF, write_geom, output_dir) #include "version_variable.h" call callTree_enter("MOM_initialize_fixed(), MOM_fixed_initialization.F90") - call log_version(PF, mod, version) + call log_version(PF, mod, version, "") call get_param(PF, mod, "DEBUG", debug, default=.false.) call get_param(PF, mod, "INPUTDIR", inputdir, & @@ -183,15 +184,11 @@ subroutine MOM_initialize_fixed(G, OBC, PF, write_geom, output_dir) end subroutine MOM_initialize_fixed ! ----------------------------------------------------------------------------- - +!> MOM_initialize_rotation makes the appropriate call to set up the Coriolis parameter. subroutine MOM_initialize_rotation(f, G, PF) - type(ocean_grid_type), intent(in) :: G - real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), intent(out) :: f - type(param_file_type), intent(in) :: PF -! Arguments: f - the Coriolis parameter in s-1. Intent out. -! (in) G - The ocean's grid structure. -! (in) PF - A structure indicating the open file to parse for -! model parameter values. + type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type + real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), intent(out) :: f !< The Coriolis parameter in s-1 + type(param_file_type), intent(in) :: PF !< Parameter file structure ! This subroutine makes the appropriate call to set up the Coriolis parameter. ! This is a separate subroutine so that it can be made public and shared with @@ -221,9 +218,11 @@ end subroutine MOM_initialize_rotation !> Calculates the components of grad f (Coriolis parameter) subroutine MOM_calculate_grad_Coriolis(dF_dx, dF_dy, G) - type(ocean_grid_type), intent(inout) :: G !< Grid type - real, dimension(SZI_(G),SZJ_(G)), intent(out) :: dF_dx !< x-component of grad f - real, dimension(SZI_(G),SZJ_(G)), intent(out) :: dF_dy !< y-component of grad f + type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type + real, dimension(G%isd:G%ied,G%jsd:G%jed), & + intent(out) :: dF_dx !< x-component of grad f + real, dimension(G%isd:G%ied,G%jsd:G%jed), & + intent(out) :: dF_dy !< y-component of grad f ! Local variables integer :: i,j real :: f1, f2 @@ -238,15 +237,13 @@ subroutine MOM_calculate_grad_Coriolis(dF_dx, dF_dy, G) call pass_vector(dF_dx, dF_dy, G%Domain, stagger=AGRID) end subroutine MOM_calculate_grad_Coriolis +!> MOM_initialize_topography makes the appropriate call to set up the bathymetry. subroutine MOM_initialize_topography(D, max_depth, G, PF) - type(ocean_grid_type), intent(in) :: G - real, dimension(SZI_(G),SZJ_(G)), intent(out) :: D - real, intent(out) :: max_depth - type(param_file_type), intent(in) :: PF -! Arguments: D - the bottom depth in m. Intent out. -! (in) G - The ocean's grid structure. -! (in) PF - A structure indicating the open file to parse for -! model parameter values. + type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type + real, dimension(G%isd:G%ied,G%jsd:G%jed), & + intent(out) :: D !< Ocean bottom depth in m + type(param_file_type), intent(in) :: PF !< Parameter file structure + real, intent(out) :: max_depth !< Maximum depth of model in m ! This subroutine makes the appropriate call to set up the bottom depth. ! This is a separate subroutine so that it can be made public and shared with @@ -285,13 +282,13 @@ subroutine MOM_initialize_topography(D, max_depth, G, PF) case ("bowl"); call initialize_topography_named(D, G, PF, config, max_depth) case ("halfpipe"); call initialize_topography_named(D, G, PF, config, max_depth) case ("DOME"); call DOME_initialize_topography(D, G, PF, max_depth) - case ("ISOMIP"); call ISOMIP_initialize_topography(D, G, PF, max_depth) + case ("ISOMIP"); call ISOMIP_initialize_topography(D, G, PF, max_depth) case ("benchmark"); call benchmark_initialize_topography(D, G, PF, max_depth) case ("DOME2D"); call DOME2d_initialize_topography(D, G, PF, max_depth) case ("sloshing"); call sloshing_initialize_topography(D, G, PF, max_depth) case ("seamount"); call seamount_initialize_topography(D, G, PF, max_depth) - case ("Phillips"); call Phillips_initialize_topography(D, G, PF) - case ("USER"); call user_initialize_topography(D, G, PF) + case ("Phillips"); call Phillips_initialize_topography(D, G, PF, max_depth) + case ("USER"); call user_initialize_topography(D, G, PF, max_depth) case default ; call MOM_error(FATAL,"MOM_initialize_topography: "// & "Unrecognized topography setup '"//trim(config)//"'") end select @@ -311,7 +308,7 @@ end subroutine MOM_initialize_topography ! ----------------------------------------------------------------------------- function diagnoseMaximumDepth(D,G) - type(ocean_grid_type), intent(in) :: G + type(dyn_horgrid_type), intent(in) :: G real, dimension(SZI_(G),SZJ_(G)), intent(in) :: D real :: diagnoseMaximumDepth ! Local variables @@ -328,8 +325,9 @@ end function diagnoseMaximumDepth !> Read gridded depths from file subroutine initialize_topography_from_file(D, G, param_file) - type(ocean_grid_type), intent(in) :: G !< Grid structure - real, dimension(SZI_(G),SZJ_(G)), intent(out) :: D !< Bottom depth (positive, in m) + type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type + real, dimension(G%isd:G%ied,G%jsd:G%jed), & + intent(out) :: D !< Ocean bottom depth in m type(param_file_type), intent(in) :: param_file !< Parameter file structure ! Local variables character(len=200) :: filename, topo_file, inputdir ! Strings for file/path @@ -369,9 +367,11 @@ end subroutine initialize_topography_from_file !> Applies a list of topography overrides read from a netcdf file subroutine apply_topography_edits_from_file(D, G, param_file) - type(ocean_grid_type), intent(in) :: G !< Grid structure - real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: D !< Bottom depth (positive, in m) + type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type + real, dimension(G%isd:G%ied,G%jsd:G%jed), & + intent(inout) :: D !< Ocean bottom depth in m type(param_file_type), intent(in) :: param_file !< Parameter file structure + ! Local variables character(len=200) :: topo_edits_file, inputdir ! Strings for file/path character(len=40) :: mod = "apply_topography_edits_from_file" ! This subroutine's name. @@ -477,12 +477,16 @@ subroutine apply_topography_edits_from_file(D, G, param_file) end subroutine apply_topography_edits_from_file ! ----------------------------------------------------------------------------- +!> initialized the bathymetry based on one of several named idealized configurations subroutine initialize_topography_named(D, G, param_file, topog_config, max_depth) - type(ocean_grid_type), intent(in) :: G - real, dimension(SZI_(G),SZJ_(G)), intent(out) :: D - type(param_file_type), intent(in) :: param_file - character(len=*), intent(in) :: topog_config - real, intent(in) :: max_depth + type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type + real, dimension(G%isd:G%ied,G%jsd:G%jed), & + intent(out) :: D !< Ocean bottom depth in m + type(param_file_type), intent(in) :: param_file !< Parameter file structure + character(len=*), intent(in) :: topog_config !< The name of an idealized + !! topographic configuration + real, intent(in) :: max_depth !< Maximum depth of model in m + ! Arguments: D - the bottom depth in m. Intent out. ! (in) G - The ocean's grid structure. ! (in) param_file - A structure indicating the open file to parse for @@ -589,11 +593,13 @@ end subroutine initialize_topography_named ! ----------------------------------------------------------------------------- ! ----------------------------------------------------------------------------- +!> limit_topography ensures that min_depth < D(x,y) < max_depth subroutine limit_topography(D, G, param_file, max_depth) - type(ocean_grid_type), intent(in) :: G - real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: D - type(param_file_type), intent(in) :: param_file - real, intent(in) :: max_depth + type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type + real, dimension(G%isd:G%ied,G%jsd:G%jed), & + intent(inout) :: D !< Ocean bottom depth in m + type(param_file_type), intent(in) :: param_file !< Parameter file structure + real, intent(in) :: max_depth !< Maximum depth of model in m ! Arguments: D - the bottom depth in m. Intent in/out. ! (in) G - The ocean's grid structure. ! (in) param_file - A structure indicating the open file to parse for @@ -638,7 +644,7 @@ end subroutine limit_topography ! ----------------------------------------------------------------------------- subroutine set_rotation_planetary(f, G, param_file) - type(ocean_grid_type), intent(in) :: G + type(dyn_horgrid_type), intent(in) :: G real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), intent(out) :: f type(param_file_type), intent(in) :: param_file ! Arguments: f - Coriolis parameter (vertical component) in s^-1 @@ -667,7 +673,7 @@ end subroutine set_rotation_planetary ! ----------------------------------------------------------------------------- subroutine set_rotation_beta_plane(f, G, param_file) - type(ocean_grid_type), intent(in) :: G + type(dyn_horgrid_type), intent(in) :: G real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), intent(out) :: f type(param_file_type), intent(in) :: param_file ! Arguments: f - Coriolis parameter (vertical component) in s^-1 @@ -713,7 +719,7 @@ end subroutine set_rotation_beta_plane ! ----------------------------------------------------------------------------- subroutine reset_face_lengths_named(G, param_file, name) - type(ocean_grid_type), intent(inout) :: G + type(dyn_horgrid_type), intent(inout) :: G type(param_file_type), intent(in) :: param_file character(len=*), intent(in) :: name ! This subroutine sets the open face lengths at selected points to restrict @@ -838,7 +844,7 @@ end subroutine reset_face_lengths_named ! ----------------------------------------------------------------------------- subroutine reset_face_lengths_file(G, param_file) - type(ocean_grid_type), intent(inout) :: G + type(dyn_horgrid_type), intent(inout) :: G type(param_file_type), intent(in) :: param_file ! This subroutine sets the open face lengths at selected points to restrict ! passages to their observed widths. @@ -906,7 +912,7 @@ end subroutine reset_face_lengths_file ! ----------------------------------------------------------------------------- subroutine reset_face_lengths_list(G, param_file) - type(ocean_grid_type), intent(inout) :: G + type(dyn_horgrid_type), intent(inout) :: G type(param_file_type), intent(in) :: param_file ! This subroutine sets the open face lengths at selected points to restrict ! passages to their observed widths. @@ -1165,7 +1171,7 @@ end subroutine read_face_length_list ! ----------------------------------------------------------------------------- subroutine set_velocity_depth_max(G) - type(ocean_grid_type), intent(inout) :: G + type(dyn_horgrid_type), intent(inout) :: G ! This subroutine sets the 4 bottom depths at velocity points to be the ! maximum of the adjacent depths. integer :: i, j @@ -1183,7 +1189,7 @@ end subroutine set_velocity_depth_max ! ----------------------------------------------------------------------------- subroutine compute_global_grid_integrals(G) - type(ocean_grid_type), intent(inout) :: G + type(dyn_horgrid_type), intent(inout) :: G ! Subroutine to pre-compute global integrals of grid quantities for ! later use in reporting diagnostics integer :: i,j @@ -1204,7 +1210,7 @@ end subroutine compute_global_grid_integrals ! ----------------------------------------------------------------------------- subroutine set_velocity_depth_min(G) - type(ocean_grid_type), intent(inout) :: G + type(dyn_horgrid_type), intent(inout) :: G ! This subroutine sets the 4 bottom depths at velocity points to be the ! minimum of the adjacent depths. integer :: i, j @@ -1222,9 +1228,9 @@ end subroutine set_velocity_depth_min ! ----------------------------------------------------------------------------- subroutine write_ocean_geometry_file(G, param_file, directory) - type(ocean_grid_type), intent(inout) :: G - type(param_file_type), intent(in) :: param_file - character(len=*), intent(in) :: directory + type(dyn_horgrid_type), intent(inout) :: G + type(param_file_type), intent(in) :: param_file + character(len=*), intent(in) :: directory ! This subroutine writes out a file containing all of the ocean geometry ! and grid data uses by the MOM ocean model. ! Arguments: G - The ocean's grid structure. Effectively intent in. @@ -1305,7 +1311,7 @@ subroutine write_ocean_geometry_file(G, param_file, directory) if (multiple_files) file_threading = MULTIPLE call create_file(unit, trim(filepath), vars, nFlds_used, fields, & - file_threading, G=G) + file_threading, dG=G) do J=Jsq,Jeq; do I=Isq,Ieq; out_q(I,J) = G%geoLatBu(I,J); enddo; enddo call write_field(unit, fields(1), G%Domain%mpp_domain, out_q) diff --git a/src/initialization/MOM_grid_initialize.F90 b/src/initialization/MOM_grid_initialize.F90 index 6bb268cd3a..540c2abd7a 100644 --- a/src/initialization/MOM_grid_initialize.F90 +++ b/src/initialization/MOM_grid_initialize.F90 @@ -72,10 +72,11 @@ module MOM_grid_initialize use MOM_domains, only : To_North, To_South, To_East, To_West use MOM_domains, only : MOM_define_domain, MOM_define_IO_domain use MOM_domains, only : MOM_domain_type +use MOM_dyn_horgrid, only : dyn_horgrid_type, set_derived_dyn_horgrid use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, is_root_pe use MOM_error_handler, only : callTree_enter, callTree_leave use MOM_file_parser, only : get_param, log_param, log_version, param_file_type -use MOM_grid, only : ocean_grid_type, set_derived_metrics +! use MOM_grid, only : ocean_grid_type, set_derived_metrics use MOM_io, only : read_data, slasher, file_exists use MOM_io, only : CORNER, NORTH_FACE, EAST_FACE @@ -107,7 +108,7 @@ module MOM_grid_initialize !! grid. The bathymetry, land-sea mask and any restricted channel widths are !! not know yet, so these are set later. subroutine set_grid_metrics(G, param_file) - type(ocean_grid_type), intent(inout) :: G !< The horizontal grid structure + type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type type(param_file_type), intent(in) :: param_file !< Parameter file structure ! Arguments: ! (inout) G - The ocean's grid structure. @@ -155,7 +156,8 @@ subroutine set_grid_metrics(G, param_file) ! Calculate derived metrics (i.e. reciprocals and products) call callTree_enter("set_derived_metrics(), MOM_grid_initialize.F90") - call set_derived_metrics(G) +! call set_derived_metrics(G) + call set_derived_dyn_horgrid(G) call callTree_leave("set_derived_metrics()") if (debug) call grid_metrics_chksum('MOM_grid_init/set_grid_metrics',G) @@ -169,7 +171,7 @@ end subroutine set_grid_metrics !! debugging. subroutine grid_metrics_chksum(parent, G) character(len=*), intent(in) :: parent !< A string identifying the caller - type(ocean_grid_type), intent(in) :: G !< The horizontal grid structure + type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type real, dimension(G%isd :G%ied ,G%jsd :G%jed ) :: tempH real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB) :: tempQ @@ -271,9 +273,10 @@ end subroutine grid_metrics_chksum ! ------------------------------------------------------------------------------ -subroutine set_grid_metrics_from_mosaic(G,param_file) - type(ocean_grid_type), intent(inout) :: G - type(param_file_type), intent(in) :: param_file +!> et_grid_metrics_from_mosaic sets the grid metrics from a mosaic file. +subroutine set_grid_metrics_from_mosaic(G, param_file) + type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type + type(param_file_type), intent(in) :: param_file !< Parameter file structure ! This subroutine sets the grid metrics from a mosaic file. ! ! Arguments: @@ -511,8 +514,9 @@ end subroutine set_grid_metrics_from_mosaic ! ------------------------------------------------------------------------------ subroutine set_grid_metrics_cartesian(G, param_file) - type(ocean_grid_type), intent(inout) :: G - type(param_file_type), intent(in) :: param_file + type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type + type(param_file_type), intent(in) :: param_file !< Parameter file structure + ! Arguments: ! (inout) G - The ocean's grid structure. ! (in) param_file - A structure indicating the open file to parse for @@ -646,8 +650,9 @@ end subroutine set_grid_metrics_cartesian ! ------------------------------------------------------------------------------ subroutine set_grid_metrics_spherical(G, param_file) - type(ocean_grid_type), intent(inout) :: G - type(param_file_type), intent(in) :: param_file + type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type + type(param_file_type), intent(in) :: param_file !< Parameter file structure + ! Arguments: ! (inout) G - The ocean's grid structure. ! (in) param_file - A structure indicating the open file to parse for @@ -785,8 +790,9 @@ end subroutine set_grid_metrics_spherical ! ------------------------------------------------------------------------------ subroutine set_grid_metrics_mercator(G, param_file) - type(ocean_grid_type), intent(inout) :: G - type(param_file_type), intent(in) :: param_file + type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type + type(param_file_type), intent(in) :: param_file !< Parameter file structure + ! Arguments: ! (inout) G - The ocean's grid structure. ! (in) param_file - A structure indicating the open file to parse for @@ -1327,8 +1333,9 @@ end function Adcroft_reciprocal !> initialize_masks initializes the grid masks and any metrics that come !! with masks already applied. subroutine initialize_masks(G, PF) - type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure - type(param_file_type), intent(in) :: PF !< The param_file handle + type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type + type(param_file_type), intent(in) :: PF !< Parameter file structure + ! Arguments: ! (inout) G - The ocean's grid structure. ! (in) PF - A structure indicating the open file to parse for diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index cf239d6685..7b227c8cc3 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -157,7 +157,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, PF, dirs, & IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB call callTree_enter("MOM_initialize_state(), MOM_state_initialization.F90") - call log_version(PF, mod, version) + call log_version(PF, mod, version, "") call get_param(PF, mod, "DEBUG", debug, default=.false.) new_sim = .false. @@ -1747,7 +1747,6 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, PF, dirs) ! model parameter values. ! (in) dirs - A structure containing several relevant directory paths. - type(ocean_grid_type), intent(inout) :: G real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(out) :: h type(thermo_var_ptrs), intent(inout) :: tv @@ -1764,12 +1763,8 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, PF, dirs) ! This include declares and sets the variable "version". #include "version_variable.h" - - character(len=40) :: mod = "MOM_initialize_layers_from_Z" ! This module's name. - - integer :: is, ie, js, je, nz ! compute domain indices integer :: isc,iec,jsc,jec ! global compute domain indices integer :: isg, ieg, jsg, jeg ! global extent @@ -1832,7 +1827,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, PF, dirs) PI_180=atan(1.0)/45. call callTree_enter(trim(mod)//"(), MOM_state_initialization.F90") - call log_version(PF, mod, version) + call log_version(PF, mod, version, "") new_sim = .false. if ((dirs%input_filename(1:1) == 'n') .and. & diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 897810c527..78939b1f75 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -298,10 +298,10 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, h_neglect = GV%H_subroundoff h_neglect3 = h_neglect**3 - if (present(OBC)) then ; if (associated(OBC)) then + if (present(OBC)) then ; if (associated(OBC)) then ; if (OBC%this_pe) then apply_OBC = OBC%apply_OBC_u_flather_east .or. OBC%apply_OBC_u_flather_west .or. & OBC%apply_OBC_v_flather_north .or. OBC%apply_OBC_v_flather_south - endif ; endif + endif ; endif ; endif if (.not.associated(CS)) call MOM_error(FATAL, & "MOM_hor_visc: Module must be initialized before it is used.") diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index ce5f0778be..f6007e3d1f 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -93,7 +93,7 @@ subroutine initialize_ALE_sponge(Iresttime, data_h, nz_data, G, param_file, CS) endif ! Set default, read and log parameters - call log_version(param_file, mod, version) + call log_version(param_file, mod, version, "") call get_param(param_file, mod, "SPONGE", use_sponge, & "If true, sponges may be applied anywhere in the domain. \n"//& "The exact location and properties of those sponges are \n"//& diff --git a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 index 9a813d95f6..815aebec9a 100644 --- a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 +++ b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 @@ -118,9 +118,9 @@ module MOM_bulk_mixed_layer integer :: ML_presort_nz_conv_adj ! If ML_resort is true, do convective ! adjustment on this many layers (starting from the ! top) before sorting the remaining layers. - logical :: use_omega ! If true, use the absolute rotation rate instead - ! of the vertical component of rotation when - ! setting the decay scale for turbulence. + real :: omega_frac ! When setting the decay scale for turbulence, use + ! this fraction of the absolute rotation rate blended + ! with the local value of f, as sqrt((1-of)*f^2 + of*4*omega^2). logical :: correct_absorption ! If true, the depth at which penetrating ! shortwave radiation is absorbed is corrected by ! moving some of the heating upward in the water @@ -1368,7 +1368,7 @@ subroutine find_starting_TKE(htot, h_CA, fluxes, Conv_En, cTKE, dKE_FC, dKE_CA, is = G%isc ; ie = G%iec ; nz = G%ke diag_wt = dt * Idt_diag - if (CS%use_omega) absf = 2.0*CS%omega + if (CS%omega_frac >= 1.0) absf = 2.0*CS%omega do i=is,ie U_Star = fluxes%ustar(i,j) if (associated(fluxes%ustar_shelf) .and. associated(fluxes%frac_shelf_h)) then @@ -1378,9 +1378,12 @@ subroutine find_starting_TKE(htot, h_CA, fluxes, Conv_En, cTKE, dKE_FC, dKE_CA, endif if (U_Star < CS%ustar_min) U_Star = CS%ustar_min - if (.not.CS%use_omega) & + if (CS%omega_frac < 1.0) then absf = 0.25*((abs(G%CoriolisBu(I,J)) + abs(G%CoriolisBu(I-1,J-1))) + & (abs(G%CoriolisBu(I,J-1)) + abs(G%CoriolisBu(I-1,J)))) + if (CS%omega_frac > 0.0) & + absf = sqrt(CS%omega_frac*4.0*CS%omega**2 + (1.0-CS%omega_frac)*absf**2) + endif absf_Ustar = absf / U_Star Idecay_len_TKE(i) = (absf_Ustar * CS%TKE_decay) * GV%H_to_m @@ -3355,8 +3358,9 @@ subroutine bulkmixedlayer_init(Time, G, GV, param_file, diag, CS) ! This include declares and sets the variable "version". #include "version_variable.h" character(len=40) :: mod = "MOM_mixed_layer" ! This module's name. + real :: omega_frac_dflt integer :: isd, ied, jsd, jed - logical :: use_temperature + logical :: use_temperature, use_omega isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed if (associated(CS)) then @@ -3446,10 +3450,20 @@ subroutine bulkmixedlayer_init(Time, G, GV, param_file, diag, CS) call get_param(param_file, mod, "OMEGA",CS%omega, & "The rotation rate of the earth.", units="s-1", & default=7.2921e-5) - call get_param(param_file, mod, "ML_USE_OMEGA", CS%use_omega, & + call get_param(param_file, mod, "ML_USE_OMEGA", use_omega, & "If true, use the absolute rotation rate instead of the \n"//& "vertical component of rotation when setting the decay \n"//& - "scale for turbulence.", default=.false.) + "scale for turbulence.", default=.false., do_not_log=.true.) + omega_frac_dflt = 0.0 + if (use_omega) then + call MOM_error(WARNING, "ML_USE_OMEGA is depricated; use ML_OMEGA_FRAC=1.0 instead.") + omega_frac_dflt = 1.0 + endif + call get_param(param_file, mod, "ML_OMEGA_FRAC", CS%omega_frac, & + "When setting the decay scale for turbulence, use this \n"//& + "fraction of the absolute rotation rate blended with the \n"//& + "local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).", & + units="nondim", default=omega_frac_dflt) call get_param(param_file, mod, "ML_RESORT", CS%ML_resort, & "If true, resort the topmost layers by potential density \n"//& "before the mixed layer calculations.", default=.false.) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 391fc1f316..ca1a325fb0 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -84,6 +84,9 @@ module MOM_diabatic_driver logical :: use_geothermal !< If true, apply geothermal heating. logical :: use_int_tides !< If true, use the code that advances a separate set !! of equations for the internal tide energy density. + logical :: ePBL_is_additive !< If true, the diffusivity from ePBL is added to all + !! other diffusivities. Otherwise, the larger of kappa- + !! shear and ePBL diffusivities are used. integer :: nMode = 1 !< Number of baroclinic modes to consider logical :: int_tide_source_test !< If true, apply an arbitrary generation site !! for internal tide testing (BDM) @@ -324,6 +327,7 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, GV, CS) integer :: ig, jg ! global indices for testing testing itide point source (BDM) logical :: avg_enabled ! for testing internal tides (BDM) + real :: Kd_add_here ! An added diffusivity in m2/s is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB @@ -741,12 +745,18 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, GV, CS) ! Augment the diffusivities due to those diagnosed in energetic_PBL. do K=2,nz ; do j=js,je ; do i=is,ie - Ent_int = Kd_ePBL(i,j,K) * (GV%m_to_H**2 * dt) / & - (0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect) + if (CS%ePBL_is_additive) then + Kd_add_here = Kd_ePBL(i,j,K) + visc%Kv_turb(i,j,K) = visc%Kv_turb(i,j,K) + Kd_ePBL(i,j,K) + else + Kd_add_here = max(Kd_ePBL(i,j,K) - visc%Kd_turb(i,j,K), 0.0) + visc%Kv_turb(i,j,K) = max(visc%Kv_turb(i,j,K), Kd_ePBL(i,j,K)) + endif + Ent_int = Kd_add_here * (GV%m_to_H**2 * dt) / & + (0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect) eb(i,j,k-1) = eb(i,j,k-1) + Ent_int ea(i,j,k) = ea(i,j,k) + Ent_int - visc%Kv_turb(i,j,K) = visc%Kv_turb(i,j,K) + Kd_ePBL(i,j,K) - Kd_int(i,j,K) = Kd_int(i,j,K) + Kd_ePBL(i,j,K) + Kd_int(i,j,K) = Kd_int(i,j,K) + Kd_add_here ! for diagnostics Kd_heat(i,j,K) = Kd_heat(i,j,K) + Kd_int(i,j,K) @@ -1751,6 +1761,10 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, "If true, use an implied energetics planetary boundary \n"//& "layer scheme to determine the diffusivity and viscosity \n"//& "in the surface boundary layer.", default=.false.) + call get_param(param_file, mod, "EPBL_IS_ADDITIVE", CS%ePBL_is_additive, & + "If true, the diffusivity from ePBL is added to all\n"//& + "other diffusivities. Otherwise, the larger of kappa-\n"//& + "shear and ePBL diffusivities are used.", default=.true.) call get_param(param_file, mod, "DOUBLE_DIFFUSION", differentialDiffusion, & "If true, apply parameterization of double-diffusion.", & default=.false. ) diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index eb1d15fced..f56aca2e39 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -99,9 +99,9 @@ module MOM_energetic_PBL ! problems, in m s-1. If the value is small enough, ! this should not affect the solution. real :: omega ! The Earth's rotation rate, in s-1. - logical :: use_omega ! If true, use the absolute rotation rate instead - ! of the vertical component of rotation when - ! setting the decay scale for turbulence. + real :: omega_frac ! When setting the decay scale for turbulence, use + ! this fraction of the absolute rotation rate blended + ! with the local value of f, as sqrt((1-of)*f^2 + of*4*omega^2). real :: wstar_ustar_coef ! A ratio relating the efficiency with which ! convectively released energy is converted to a ! turbulent velocity, relative to mechanically @@ -443,10 +443,12 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & endif if (U_Star < CS%ustar_min) U_Star = CS%ustar_min - if (CS%use_omega) then ; absf(i) = 2.0*CS%omega + if (CS%omega_frac >= 1.0) then ; absf(i) = 2.0*CS%omega else absf(i) = 0.25*((abs(G%CoriolisBu(I,J)) + abs(G%CoriolisBu(I-1,J-1))) + & (abs(G%CoriolisBu(I,J-1)) + abs(G%CoriolisBu(I-1,J)))) + if (CS%omega_frac > 0.0) & + absf = sqrt(CS%omega_frac*4.0*CS%omega**2 + (1.0-CS%omega_frac)*absf**2) endif mech_TKE(i) = (dt*CS%mstar*GV%Rho0)*((U_Star**3)) @@ -1117,8 +1119,9 @@ subroutine energetic_PBL_init(Time, G, GV, param_file, diag, CS) ! This include declares and sets the variable "version". #include "version_variable.h" character(len=40) :: mod = "MOM_energetic_PBL" ! This module's name. + real :: omega_frac_dflt integer :: isd, ied, jsd, jed - logical :: use_temperature + logical :: use_temperature, use_omega isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed if (associated(CS)) then @@ -1157,10 +1160,20 @@ subroutine energetic_PBL_init(Time, G, GV, param_file, diag, CS) call get_param(param_file, mod, "OMEGA",CS%omega, & "The rotation rate of the earth.", units="s-1", & default=7.2921e-5) - call get_param(param_file, mod, "ML_USE_OMEGA", CS%use_omega, & + call get_param(param_file, mod, "ML_USE_OMEGA", use_omega, & "If true, use the absolute rotation rate instead of the \n"//& "vertical component of rotation when setting the decay \n"//& - "scale for turbulence.", default=.false.) + "scale for turbulence.", default=.false., do_not_log=.true.) + omega_frac_dflt = 0.0 + if (use_omega) then + call MOM_error(WARNING, "ML_USE_OMEGA is depricated; use ML_OMEGA_FRAC=1.0 instead.") + omega_frac_dflt = 1.0 + endif + call get_param(param_file, mod, "ML_OMEGA_FRAC", CS%omega_frac, & + "When setting the decay scale for turbulence, use this \n"//& + "fraction of the absolute rotation rate blended with the \n"//& + "local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).", & + units="nondim", default=omega_frac_dflt) call get_param(param_file, mod, "WSTAR_USTAR_COEF", CS%wstar_ustar_coef, & "A ratio relating the efficiency with which convectively\n"//& "released energy is converted to a turbulent velocity,\n"//& diff --git a/src/parameterizations/vertical/MOM_internal_tide_input.F90 b/src/parameterizations/vertical/MOM_internal_tide_input.F90 index e637a48d68..ba1f90fc06 100644 --- a/src/parameterizations/vertical/MOM_internal_tide_input.F90 +++ b/src/parameterizations/vertical/MOM_internal_tide_input.F90 @@ -325,7 +325,7 @@ subroutine int_tide_input_init(Time, G, GV, param_file, diag, CS, itide) CS%diag => diag ! Read all relevant parameters and write them to the model log. - call log_version(param_file, mod, version) + call log_version(param_file, mod, version, "") call get_param(param_file, mod, "INPUTDIR", CS%inputdir, default=".") CS%inputdir = slasher(CS%inputdir) diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index cc8059355a..1cab09531f 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -255,6 +255,9 @@ module MOM_set_diffusivity logical :: ML_use_omega ! If true, use absolute rotation rate instead ! of the vertical component of rotation when ! setting the decay scale for mixed layer turbulence. + real :: ML_omega_frac ! When setting the decay scale for turbulence, use + ! this fraction of the absolute rotation rate blended + ! with the local value of f, as f^2 ~= (1-of)*f^2 + of*4*omega^2. logical :: user_change_diff ! If true, call user-defined code to change diffusivity. logical :: useKappaShear ! If true, use the kappa_shear module to find the ! shear-driven diapycnal diffusivity. @@ -1791,11 +1794,13 @@ subroutine add_MLrad_diffusivity(h, fluxes, j, G, GV, CS, Kd, TKE_to_Kd, Kd_int) do k=1,kml ; do i=is,ie ; h_ml(i) = h_ml(i) + GV%H_to_m*h(i,j,k) ; enddo ; enddo do i=is,ie ; if (do_i(i)) then - if (CS%ML_use_omega) then + if (CS%ML_omega_frac >= 1.0) then f_sq = 4.0*Omega2 else f_sq = 0.25*((G%CoriolisBu(I,J)**2 + G%CoriolisBu(I-1,J-1)**2) + & (G%CoriolisBu(I,J-1)**2 + G%CoriolisBu(I-1,J)**2)) + if (CS%ML_omega_frac > 0.0) & + f_sq = CS%ML_omega_frac*4.0*Omega2 + (1.0-CS%ML_omega_frac)*f_sq endif ustar_sq = max(fluxes%ustar(i,j), CS%ustar_min)**2 @@ -2509,13 +2514,14 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp real :: decay_length, utide, zbot, hamp type(vardesc) :: vd - logical :: read_tideamp + logical :: read_tideamp, ML_use_omega ! This include declares and sets the variable "version". #include "version_variable.h" character(len=40) :: mod = "MOM_set_diffusivity" ! This module's name. character(len=20) :: tmpstr character(len=200) :: filename, tideamp_file, h2_file, Niku_TKE_input_file real :: Niku_scale ! local variable for scaling the Nikurashin TKE flux data + real :: omega_frac_dflt integer :: i, j, is, ie, js, je if (associated(CS)) then @@ -2538,7 +2544,7 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp ! Read all relevant parameters and write them to the model log. - call log_version(param_file, mod, version) + call log_version(param_file, mod, version, "") call get_param(param_file, mod, "INPUTDIR", CS%inputdir, default=".") CS%inputdir = slasher(CS%inputdir) @@ -2587,11 +2593,20 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp call get_param(param_file, mod, "TKE_DECAY", CS%TKE_decay, & "The ratio of the natural Ekman depth to the TKE decay scale.", & units="nondim", default=2.5) - call get_param(param_file, mod, "ML_USE_OMEGA", CS%ML_use_omega, & + call get_param(param_file, mod, "ML_USE_OMEGA", ML_use_omega, & "If true, use the absolute rotation rate instead of the \n"//& "vertical component of rotation when setting the decay \n"//& - "scale for turbulence in the mixed layer. \n"//& - "This is only used if ML_RADIATION is true.", default=.false.) + "scale for turbulence.", default=.false., do_not_log=.true.) + omega_frac_dflt = 0.0 + if (ML_use_omega) then + call MOM_error(WARNING, "ML_USE_OMEGA is depricated; use ML_OMEGA_FRAC=1.0 instead.") + omega_frac_dflt = 1.0 + endif + call get_param(param_file, mod, "ML_OMEGA_FRAC", CS%ML_omega_frac, & + "When setting the decay scale for turbulence, use this \n"//& + "fraction of the absolute rotation rate blended with the \n"//& + "local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).", & + units="nondim", default=omega_frac_dflt) endif call get_param(param_file, mod, "BOTTOMDRAGLAW", CS%bottomdraglaw, & diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index 77434af853..6edfdb6565 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -116,9 +116,9 @@ module MOM_set_visc ! this should not affect the solution. real :: TKE_decay ! The ratio of the natural Ekman depth to the TKE ! decay scale, nondimensional. - logical :: use_omega ! If true, use the absolute rotation rate instead - ! of the vertical component of rotation when - ! setting the decay scale for turbulence. + real :: omega_frac ! When setting the decay scale for turbulence, use + ! this fraction of the absolute rotation rate blended + ! with the local value of f, as sqrt((1-of)*f^2 + of*4*omega^2). logical :: debug ! If true, write verbose checksums for debugging purposes. type(diag_ctrl), pointer :: diag ! A structure that is used to regulate the ! timing of diagnostic output. @@ -1064,8 +1064,11 @@ subroutine set_viscous_ML(u, v, h, tv, fluxes, visc, dt, G, GV, CS) vhtot(I) = 0.25 * dt_Rho0 * ((fluxes%tauy(i,J) + fluxes%tauy(i+1,J-1)) + & (fluxes%tauy(i,J-1) + fluxes%tauy(i+1,J))) - if (CS%use_omega) then ; absf = 2.0*CS%omega - else ; absf = 0.5*(abs(G%CoriolisBu(I,J)) + abs(G%CoriolisBu(I,J-1))) ; endif + if (CS%omega_frac >= 1.0) then ; absf = 2.0*CS%omega ; else + absf = 0.5*(abs(G%CoriolisBu(I,J)) + abs(G%CoriolisBu(I,J-1))) + if (CS%omega_frac > 0.0) & + absf = sqrt(CS%omega_frac*4.0*CS%omega**2 + (1.0-CS%omega_frac)*absf**2) + endif U_Star = max(CS%ustar_min, 0.5 * (fluxes%ustar(i,j) + fluxes%ustar(i+1,j))) Idecay_len_TKE(I) = ((absf / U_Star) * CS%TKE_decay) * H_to_m endif @@ -1304,8 +1307,12 @@ subroutine set_viscous_ML(u, v, h, tv, fluxes, visc, dt, G, GV, CS) uhtot(i) = 0.25 * dt_Rho0 * ((fluxes%taux(I,j) + fluxes%tauy(I-1,j+1)) + & (fluxes%taux(I-1,j) + fluxes%tauy(I,j+1))) - if (CS%use_omega) then ; absf = 2.0*CS%omega - else ; absf = 0.5*(abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J))) ; endif + if (CS%omega_frac >= 1.0) then ; absf = 2.0*CS%omega ; else + absf = 0.5*(abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J))) + if (CS%omega_frac > 0.0) & + absf = sqrt(CS%omega_frac*4.0*CS%omega**2 + (1.0-CS%omega_frac)*absf**2) + endif + U_Star = max(CS%ustar_min, 0.5 * (fluxes%ustar(i,j) + fluxes%ustar(i,j+1))) Idecay_len_TKE(i) = ((absf / U_Star) * CS%TKE_decay) * H_to_m @@ -1617,8 +1624,9 @@ subroutine set_visc_init(Time, G, GV, param_file, diag, visc, CS) ! for this module real :: Csmag_chan_dflt, smag_const1, TKE_decay_dflt, bulk_Ri_ML_dflt real :: Kv_background + real :: omega_frac_dflt integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz - logical :: use_kappa_shear, adiabatic, differential_diffusion + logical :: use_kappa_shear, adiabatic, differential_diffusion, use_omega ! This include declares and sets the variable "version". #include "version_variable.h" character(len=40) :: mod = "MOM_set_visc" ! This module's name. @@ -1695,10 +1703,20 @@ subroutine set_visc_init(Time, G, GV, param_file, diag, visc, CS) "mixed layer viscosity. By default, \n"//& "TKE_DECAY_VISC = TKE_DECAY or 0.", units="nondim", & default=TKE_decay_dflt) - call get_param(param_file, mod, "ML_USE_OMEGA", CS%use_omega, & + call get_param(param_file, mod, "ML_USE_OMEGA", use_omega, & "If true, use the absolute rotation rate instead of the \n"//& "vertical component of rotation when setting the decay \n"//& - "scale for turbulence.", default=.false.) + "scale for turbulence.", default=.false., do_not_log=.true.) + omega_frac_dflt = 0.0 + if (use_omega) then + call MOM_error(WARNING, "ML_USE_OMEGA is depricated; use ML_OMEGA_FRAC=1.0 instead.") + omega_frac_dflt = 1.0 + endif + call get_param(param_file, mod, "ML_OMEGA_FRAC", CS%omega_frac, & + "When setting the decay scale for turbulence, use this \n"//& + "fraction of the absolute rotation rate blended with the \n"//& + "local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).", & + units="nondim", default=omega_frac_dflt) call get_param(param_file, mod, "OMEGA", CS%omega, & "The rotation rate of the earth.", units="s-1", & default=7.2921e-5) diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index 7f32001ace..cca1f7807a 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -388,7 +388,7 @@ subroutine vertvisc(u, v, h, fluxes, visc, dt, OBC, ADp, CDp, G, GV, CS, & call vertvisc_limit_vel(u, v, h, ADp, CDp, fluxes, visc, dt, G, GV, CS) ! Here the velocities associated with open boundary conditions are applied. - if (associated(OBC)) then + if (associated(OBC)) then ; if (OBC%this_pe) then if (OBC%apply_OBC_u) then do k=1,nz ; do j=G%jsc,G%jec ; do I=Isq,Ieq if (OBC%OBC_mask_u(I,j) .and. (OBC%OBC_kind_u(I,j) == OBC_SIMPLE)) & @@ -401,7 +401,7 @@ subroutine vertvisc(u, v, h, fluxes, visc, dt, OBC, ADp, CDp, G, GV, CS, & v(i,J,k) = OBC%v(i,J,k) enddo ; enddo ; enddo endif - endif + endif ; endif ! Offer diagnostic fields for averaging. if (CS%id_du_dt_visc > 0) & call post_data(CS%id_du_dt_visc, ADp%du_dt_visc, CS%diag) diff --git a/src/tracer/MOM_tracer_advect.F90 b/src/tracer/MOM_tracer_advect.F90 index 674a5f1e41..9686dc7eb2 100644 --- a/src/tracer/MOM_tracer_advect.F90 +++ b/src/tracer/MOM_tracer_advect.F90 @@ -470,7 +470,7 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & enddo ; enddo endif ! usePPM - if (associated(OBC)) then ; if (OBC%apply_OBC_u) then + if (associated(OBC)) then ; if (OBC%this_pe) then ; if (OBC%apply_OBC_u) then do_any_i = .false. do I=is-1,ie do_i(I) = .false. @@ -491,7 +491,7 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & flux_x(I,m) = uhh(I)*Tr(m)%OBC_in_u(I,j,k) else ; flux_x(I,m) = uhh(I)*Tr(m)%OBC_inflow_conc ; endif endif ; enddo ; enddo ; endif - endif ; endif + endif ; endif ; endif ! Calculate new tracer concentration in each cell after accounting ! for the i-direction fluxes. @@ -730,7 +730,7 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & enddo ; enddo endif ! usePPM - if (associated(OBC)) then ; if (OBC%apply_OBC_v) then + if (associated(OBC)) then ; if (OBC%this_pe) then ; if (OBC%apply_OBC_v) then do_any_i = .false. do i=is,ie do_i(i) = .false. @@ -751,7 +751,7 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & flux_y(i,m,J) = vhh(i,J)*Tr(m)%OBC_in_v(i,J,k) else ; flux_y(i,m,J) = vhh(i,J)*Tr(m)%OBC_inflow_conc ; endif endif ; enddo ; enddo ; endif - endif ; endif + endif ; endif ; endif else ! not domore_v. do i=is,ie ; vhh(i,J) = 0.0 ; enddo do m=1,ntr ; do i=is,ie ; flux_y(i,m,J) = 0.0 ; enddo ; enddo diff --git a/src/user/DOME2d_initialization.F90 b/src/user/DOME2d_initialization.F90 index fcf8367bcb..fab4e2f2dc 100644 --- a/src/user/DOME2d_initialization.F90 +++ b/src/user/DOME2d_initialization.F90 @@ -1,6 +1,7 @@ module DOME2d_initialization use MOM_ALE_sponge, only : ALE_sponge_CS, set_up_ALE_sponge_field, initialize_ALE_sponge +use MOM_dyn_horgrid, only : dyn_horgrid_type use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_get_input, only : directories @@ -33,10 +34,11 @@ module DOME2d_initialization !> Initialize topography with a shelf and slope in a 2D domain subroutine DOME2d_initialize_topography ( D, G, param_file, max_depth ) ! Arguments - type(ocean_grid_type), intent(in) :: G !< Ocean grid structure - real, dimension(SZI_(G),SZJ_(G)), intent(out) :: D !< Ocean bottom depth - type(param_file_type), intent(in) :: param_file !< Parameter file structure - real, intent(in) :: max_depth !< Maximum depth of model + type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type + real, dimension(G%isd:G%ied,G%jsd:G%jed), & + intent(out) :: D !< Ocean bottom depth in m + type(param_file_type), intent(in) :: param_file !< Parameter file structure + real, intent(in) :: max_depth !< Maximum depth of model in m ! Local variables integer :: i, j real :: x, bay_depth, l1, l2 @@ -82,8 +84,8 @@ end subroutine DOME2d_initialize_topography !> Initialize thicknesses according to coordinate mode subroutine DOME2d_initialize_thickness ( h, G, GV, param_file ) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure - real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: h !< Layer thicknesses type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: h !< Layer thicknesses type(param_file_type), intent(in) :: param_file !< Parameter file structure ! Local variables real :: e0(SZK_(G)) ! The resting interface heights, in m, usually ! diff --git a/src/user/DOME_initialization.F90 b/src/user/DOME_initialization.F90 index 44d0ec94d0..334fbd8208 100644 --- a/src/user/DOME_initialization.F90 +++ b/src/user/DOME_initialization.F90 @@ -20,6 +20,7 @@ module DOME_initialization !*********************************************************************** use MOM_sponge, only : sponge_CS, set_up_sponge_field, initialize_sponge +use MOM_dyn_horgrid, only : dyn_horgrid_type use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, is_root_pe use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_get_input, only : directories @@ -45,11 +46,11 @@ module DOME_initialization ! ----------------------------------------------------------------------------- !> This subroutine sets up the DOME topography subroutine DOME_initialize_topography(D, G, param_file, max_depth) - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - real, intent(out), dimension(SZI_(G),SZJ_(G)) :: D !< The bottom depth in m. - type(param_file_type), intent(in) :: param_file !< A structure indicating the open file - !! to parse for model parameter values. - real, intent(in) :: max_depth !< Maximum depth. + type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type + real, dimension(G%isd:G%ied,G%jsd:G%jed), & + intent(out) :: D !< Ocean bottom depth in m + type(param_file_type), intent(in) :: param_file !< Parameter file structure + real, intent(in) :: max_depth !< Maximum depth of model in m real :: min_depth ! The minimum and maximum depths in m. ! This include declares and sets the variable "version". @@ -92,7 +93,7 @@ end subroutine DOME_initialize_topography subroutine DOME_initialize_thickness(h, G, GV, param_file) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. - real, intent(out), dimension(SZI_(G),SZJ_(G), SZK_(G)) :: h !< The thickness that is being initialized. + real, intent(out), dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h !< The thickness that is being initialized. type(param_file_type), intent(in) :: param_file !< A structure indicating the open file !! to parse for model parameter values. @@ -234,7 +235,7 @@ end subroutine DOME_initialize_sponges !> Set the positions of the open boundary needed for the DOME experiment. subroutine DOME_set_OBC_positions(G, param_file, OBC) - type(ocean_grid_type), intent(in) :: G !< Grid structure. + type(dyn_horgrid_type), intent(in) :: G !< Grid structure. type(param_file_type), intent(in) :: param_file !< Parameter file handle. type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure. ! Local variables diff --git a/src/user/ISOMIP_initialization.F90 b/src/user/ISOMIP_initialization.F90 index f4f8e4d326..4a21b31324 100644 --- a/src/user/ISOMIP_initialization.F90 +++ b/src/user/ISOMIP_initialization.F90 @@ -20,6 +20,7 @@ module ISOMIP_initialization !*********************************************************************** use MOM_ALE_sponge, only : ALE_sponge_CS, set_up_ALE_sponge_field, initialize_ALE_sponge +use MOM_dyn_horgrid, only : dyn_horgrid_type use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, is_root_pe, WARNING use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_get_input, only : directories @@ -58,12 +59,11 @@ module ISOMIP_initialization !> Initialization of topography subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth) - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - real, intent(out), dimension(SZI_(G),SZJ_(G)) :: D !< The bottom depth in m. - type(param_file_type), intent(in) :: param_file !< A structure indicating the - !! open file to parse for model - !! parameter values. - real, intent(in) :: max_depth !< Maximum depth. + type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type + real, dimension(G%isd:G%ied,G%jsd:G%jed), & + intent(out) :: D !< Ocean bottom depth in m + type(param_file_type), intent(in) :: param_file !< Parameter file structure + real, intent(in) :: max_depth !< Maximum depth of model in m ! This subroutine sets up the ISOMIP topography real :: min_depth ! The minimum and maximum depths in m. @@ -146,7 +146,7 @@ end subroutine ISOMIP_initialize_topography subroutine ISOMIP_initialize_thickness ( h, G, GV, param_file, tv ) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. - real, intent(out), dimension(SZI_(G),SZJ_(G), SZK_(G)) :: h !< The thickness that is being + real, intent(out), dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h !< The thickness that is being !! initialized. type(param_file_type), intent(in) :: param_file !< A structure indicating the !! open file to parse for model diff --git a/src/user/Phillips_initialization.F90 b/src/user/Phillips_initialization.F90 index ad176ae4d3..3f85c2ee06 100644 --- a/src/user/Phillips_initialization.F90 +++ b/src/user/Phillips_initialization.F90 @@ -20,6 +20,7 @@ module Phillips_initialization !*********************************************************************** use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, is_root_pe +use MOM_dyn_horgrid, only : dyn_horgrid_type use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_get_input, only : directories use MOM_grid, only : ocean_grid_type @@ -49,9 +50,9 @@ module Phillips_initialization !> Initialize thickness field. subroutine Phillips_initialize_thickness(h, G, GV, param_file) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - real, intent(out), dimension(SZI_(G),SZJ_(G), SZK_(G)) :: h !< The thickness that is - !! being initialized. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(out) :: h !< The thickness that is being initialized. type(param_file_type), intent(in) :: param_file !< A structure indicating the !! open file to parse for model !! parameter values. @@ -273,12 +274,12 @@ function sech(x) end function sech !> Initialize topography. -subroutine Phillips_initialize_topography(D, G, param_file) - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - real, intent(out), dimension(SZI_(G),SZJ_(G)) :: D !< The bottom depth in m. - type(param_file_type), intent(in) :: param_file !< A structure indicating the - !! open file to parse for model - !! parameter values. +subroutine Phillips_initialize_topography(D, G, param_file, max_depth) + type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type + real, dimension(G%isd:G%ied,G%jsd:G%jed), & + intent(out) :: D !< Ocean bottom depth in m + type(param_file_type), intent(in) :: param_file !< Parameter file structure + real, intent(in) :: max_depth !< Maximum depth of model in m real :: PI, Htop, Wtop, Ltop, offset, dist, & x1, x2, x3, x4, y1, y2 @@ -314,7 +315,7 @@ subroutine Phillips_initialize_topography(D, G, param_file) D(i,j) = 2.0/3.0*Htop*sin(PI*(G%geoLonT(i,j)-x3)/(x4-x3))**2 & *sin(PI*(G%geoLatT(i,j)-y1)/(y2-y1))**2 end if - D(i,j)=G%max_depth-D(i,j) + D(i,j)=max_depth-D(i,j) enddo; enddo end subroutine Phillips_initialize_topography diff --git a/src/user/benchmark_initialization.F90 b/src/user/benchmark_initialization.F90 index c39c9db0f4..6415db7a13 100644 --- a/src/user/benchmark_initialization.F90 +++ b/src/user/benchmark_initialization.F90 @@ -20,6 +20,7 @@ module benchmark_initialization !*********************************************************************** use MOM_sponge, only : sponge_CS, set_up_sponge_field, initialize_sponge +use MOM_dyn_horgrid, only : dyn_horgrid_type use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, is_root_pe use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_get_input, only : directories @@ -42,11 +43,11 @@ module benchmark_initialization ! ----------------------------------------------------------------------------- !> This subroutine sets up the benchmark test case topography. subroutine benchmark_initialize_topography(D, G, param_file, max_depth) - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - real, dimension(SZI_(G),SZJ_(G)), intent(out) :: D !< the bottom depth in m. - type(param_file_type), intent(in) :: param_file !< A structure indicating the open - !! file to parse for model parameter values. - real, intent(in) :: max_depth !< The Maximum depth. + type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type + real, dimension(G%isd:G%ied,G%jsd:G%jed), & + intent(out) :: D !< Ocean bottom depth in m + type(param_file_type), intent(in) :: param_file !< Parameter file structure + real, intent(in) :: max_depth !< Maximum depth of model in m ! This subroutine sets up the benchmark test case topography real :: min_depth ! The minimum and maximum depths in m. @@ -93,7 +94,7 @@ end subroutine benchmark_initialize_topography subroutine benchmark_initialize_thickness(h, G, GV, param_file, eqn_of_state, P_ref) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. - real, intent(out), dimension(SZI_(G),SZJ_(G), SZK_(G)) :: h !< The thickness that is being + real, intent(out), dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h !< The thickness that is being !! initialized. type(param_file_type), intent(in) :: param_file !< A structure indicating the open !! file to parse for model diff --git a/src/user/seamount_initialization.F90 b/src/user/seamount_initialization.F90 index e25d05c5c1..3fa5cb91c7 100644 --- a/src/user/seamount_initialization.F90 +++ b/src/user/seamount_initialization.F90 @@ -20,6 +20,7 @@ module seamount_initialization !*********************************************************************** use MOM_domains, only : sum_across_PEs +use MOM_dyn_horgrid, only : dyn_horgrid_type use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, is_root_pe use MOM_file_parser, only : get_param, param_file_type use MOM_get_input, only : directories @@ -57,12 +58,11 @@ module seamount_initialization !> Initialization of topography. subroutine seamount_initialize_topography ( D, G, param_file, max_depth ) ! Arguments - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - real, intent(out), dimension(SZI_(G),SZJ_(G)) :: D !< The bottom depth in m. - type(param_file_type), intent(in) :: param_file !< A structure indicating the - !! open file to parse for model - !! parameter values. - real, intent(in) :: max_depth !< Maximum depth. + type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type + real, dimension(G%isd:G%ied,G%jsd:G%jed), & + intent(out) :: D !< Ocean bottom depth in m + type(param_file_type), intent(in) :: param_file !< Parameter file structure + real, intent(in) :: max_depth !< Maximum depth of model in m ! Local variables integer :: i, j @@ -99,9 +99,9 @@ end subroutine seamount_initialize_topography !! This subroutine initializes the layer thicknesses to be uniform. subroutine seamount_initialize_thickness ( h, G, GV, param_file ) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - real, intent(out), dimension(SZI_(G),SZJ_(G), SZK_(G)) :: h !< The thicknesses being - !! initialized. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + real, intent(out), dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h !< The thicknesses being + !! initialized. type(param_file_type), intent(in) :: param_file !< A structure indicating the !! open file to parse for model !! parameter values. diff --git a/src/user/sloshing_initialization.F90 b/src/user/sloshing_initialization.F90 index e0a48fa135..34bc7edc46 100644 --- a/src/user/sloshing_initialization.F90 +++ b/src/user/sloshing_initialization.F90 @@ -20,6 +20,7 @@ module sloshing_initialization !*********************************************************************** use MOM_domains, only : sum_across_PEs +use MOM_dyn_horgrid, only : dyn_horgrid_type use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, is_root_pe use MOM_file_parser, only : get_param, param_file_type use MOM_get_input, only : directories @@ -53,12 +54,11 @@ module sloshing_initialization !> Initialization of topography. subroutine sloshing_initialize_topography ( D, G, param_file, max_depth ) - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - real, intent(out), dimension(SZI_(G),SZJ_(G)) :: D !< The bottom depth in m. - type(param_file_type), intent(in) :: param_file !< A structure indicating the - !! open file to parse for model - !! parameter values. - real, intent(in) :: max_depth !< Maximum depth. + type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type + real, dimension(G%isd:G%ied,G%jsd:G%jed), & + intent(out) :: D !< Ocean bottom depth in m + type(param_file_type), intent(in) :: param_file !< Parameter file structure + real, intent(in) :: max_depth !< Maximum depth of model in m ! Local variables integer :: i, j @@ -84,9 +84,9 @@ end subroutine sloshing_initialize_topography !! left and minimum value on the right. This sets off a regular sloshing motion. subroutine sloshing_initialize_thickness ( h, G, GV, param_file ) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - real, intent(out), dimension(SZI_(G),SZJ_(G), SZK_(G)) :: h !< The thicknesses being - !! initialized. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + real, intent(out), dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h !< The thicknesses being + !! initialized. type(param_file_type), intent(in) :: param_file !< A structure indicating the !! open file to parse for model !! parameter values. diff --git a/src/user/user_change_diffusivity.F90 b/src/user/user_change_diffusivity.F90 index ccab939e87..b58b8bce71 100644 --- a/src/user/user_change_diffusivity.F90 +++ b/src/user/user_change_diffusivity.F90 @@ -226,7 +226,7 @@ subroutine user_change_diff_init(Time, G, param_file, diag, CS) CS%diag => diag ! Read all relevant parameters and write them to the model log. - call log_version(param_file, mod, version) + call log_version(param_file, mod, version, "") call get_param(param_file, mod, "USER_KD_ADD", CS%Kd_add, & "A user-specified additional diffusivity over a range of \n"//& "latitude and density.", units="m2 s-1", default=0.0) diff --git a/src/user/user_initialization.F90 b/src/user/user_initialization.F90 index 1481c8845c..c0868ef543 100644 --- a/src/user/user_initialization.F90 +++ b/src/user/user_initialization.F90 @@ -20,6 +20,7 @@ module user_initialization !*********************************************************************** use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, is_root_pe +use MOM_dyn_horgrid, only : dyn_horgrid_type use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_get_input, only : directories use MOM_grid, only : ocean_grid_type @@ -70,12 +71,13 @@ subroutine USER_set_coord(Rlay, g_prime, GV, param_file, eqn_of_state) end subroutine USER_set_coord !> Initialize topography. -subroutine USER_initialize_topography(D, G, param_file) - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - real, intent(out), dimension(SZI_(G),SZJ_(G)) :: D !< The bottom depth in m. - type(param_file_type), intent(in) :: param_file !< A structure indicating the - !! open file to parse for model - !! parameter values. +subroutine USER_initialize_topography(D, G, param_file, max_depth) + type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type + real, dimension(G%isd:G%ied,G%jsd:G%jed), & + intent(out) :: D !< Ocean bottom depth in m + type(param_file_type), intent(in) :: param_file !< Parameter file structure + real, intent(in) :: max_depth !< Maximum depth of model in m + call MOM_error(FATAL, & "USER_initialization.F90, USER_initialize_topography: " // & "Unmodified user routine called - you must edit the routine to use it") @@ -89,7 +91,7 @@ end subroutine USER_initialize_topography !> initialize thicknesses. subroutine USER_initialize_thickness(h, G, param_file, T) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - real, intent(out), dimension(SZI_(G),SZJ_(G), SZK_(G)) :: h !< The thicknesses being + real, intent(out), dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h !< The thicknesses being !! initialized. type(param_file_type), intent(in) :: param_file !< A structure indicating the !! open file to parse for model @@ -199,7 +201,7 @@ end subroutine USER_initialize_sponges !> This subroutine sets the location of open boundaries. subroutine USER_set_OBC_positions(G, param_file, OBC) - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(dyn_horgrid_type), intent(in) :: G !< The ocean's grid structure. type(param_file_type), intent(in) :: param_file !< A structure indicating the !! open file to parse for model !! parameter values.