diff --git a/.testing/tc4/MOM_input b/.testing/tc4/MOM_input index 04598a9dc9..a9fde1a6c7 100644 --- a/.testing/tc4/MOM_input +++ b/.testing/tc4/MOM_input @@ -15,34 +15,31 @@ DT_THERM = 3600.0 ! [s] default = 300.0 ! The thermodynamic and tracer advection time step. Ideally DT_THERM should be ! an integer multiple of DT and less than the forcing or coupling time-step, ! unless THERMO_SPANS_COUPLING is true, in which case DT_THERM can be an integer - ! multiple of the coupling timestep. By default DT_THERM is set to DT. + ! multiple of the coupling timestep. By default DT_THERM is set to DT. C_P = 3925.0 ! [J kg-1 K-1] default = 3991.86795711963 ! The heat capacity of sea water, approximated as a constant. This is only used - ! if ENABLE_THERMODYNAMICS is true. The default value is from the TEOS-10 - ! definition of conservative temperature. + ! if ENABLE_THERMODYNAMICS is true. The default value is from the TEOS-10 +INTERPOLATE_P_SURF = False ! [Boolean] default = False + ! If true, linearly interpolate the surface pressure over the coupling time + ! step, using the specified value at the end of the step. SAVE_INITIAL_CONDS = False ! [Boolean] default = False ! If true, write the initial conditions to a file given by IC_OUTPUT_FILE. - +ICE_THICKNESS_FILE = "shelf.nc" ! + ! The file from which the ice bathymetry and area are read. +ICE_AREA_VARNAME = "area" ! + ! The name of the area variable in ICE_THICKNESS_FILE. ! === module MOM_domains === REENTRANT_X = False ! [Boolean] default = True ! If true, the domain is zonally reentrant. +REENTRANT_Y = False ! [Boolean] default = False + ! If true, the domain is meridionally reentrant. NIGLOBAL = 14 ! ! The total number of thickness grid points in the x-direction in the physical ! domain. With STATIC_MEMORY_ this is set in MOM_memory.h at compile time. -NJGLOBAL = 10 ! +NJGLOBAL = 10 ! ! The total number of thickness grid points in the y-direction in the physical ! domain. With STATIC_MEMORY_ this is set in MOM_memory.h at compile time. -! === module MOM_hor_index === -! Sets the horizontal array index types. - -! === module MOM_verticalGrid === -! Parameters providing information about the vertical grid. -NK = 2 ! [nondim] - ! The number of model layers. - -! === module MOM_fixed_initialization === - ! === module MOM_grid_init === GRID_CONFIG = "mosaic" ! ! A character string that determines the method for defining the horizontal @@ -65,8 +62,9 @@ TOPO_CONFIG = "file" ! ! wall at the southern face. ! halfpipe - a zonally uniform channel with a half-sine ! profile in the meridional direction. + ! bbuilder - build topography from list of functions. ! benchmark - use the benchmark test case topography. - ! Neverland - use the Neverland test case topography. + ! Neverworld - use the Neverworld test case topography. ! DOME - use a slope and channel configuration for the ! DOME sill-overflow test case. ! ISOMIP - use a slope and channel configuration for the @@ -80,38 +78,59 @@ TOPO_CONFIG = "file" ! ! Phillips - ACC-like idealized topography used in the Phillips config. ! dense - Denmark Strait-like dense water formation and overflow. ! USER - call a user modified routine. -!MAXIMUM_DEPTH = 100.0 ! [m] - ! The (diagnosed) maximum depth of the ocean. - -! === module MOM_open_boundary === -! Controls where open boundaries are located, what kind of boundary condition to impose, and what data to apply, -! if any. ROTATION = "betaplane" ! default = "2omegasinlat" ! This specifies how the Coriolis parameter is specified: ! 2omegasinlat - Use twice the planetary rotation rate ! times the sine of latitude. ! betaplane - Use a beta-plane or f-plane. ! USER - call a user modified routine. -F_0 = 1.0E-04 ! [s-1] default = 0.0 - ! The reference value of the Coriolis parameter with the betaplane option. -! === module MOM_tracer_registry === +! === module MOM_open_boundary === +! Controls where open boundaries are located, what kind of boundary condition to impose, and what data to apply, if any. +OBC_NUMBER_OF_SEGMENTS = 0 ! default = 0 + ! The number of open boundary segments. +OBC_FREESLIP_VORTICITY = True ! [Boolean] default = False + ! If true, sets the normal gradient of tangential velocity to + ! zero in the relative vorticity on open boundaries. This cannot + ! be true if OBC_ZERO_VORTICITY is True. +OBC_FREESLIP_STRAIN = True ! [Boolean] default = False + ! If true, sets the normal gradient of tangential velocity to + ! zero in the strain use in the stress tensor on open boundaries. This cannot + ! be true if OBC_ZERO_STRAIN is True. +OBC_ZERO_BIHARMONIC = True ! [Boolean] default = False + ! If true, zeros the Laplacian of flow on open boundaries in the biharmonic + ! viscosity term. +OBC_SEGMENT_001 = "I=N,J=N:0,FLATHER,ORLANSKI" ! +OBC_SEGMENT_001_DATA = "U=value:0.0,V=value:0.0,TEMP=value:0.0,SALT=value:35.0,SSH=value:0.0" +MASK_OUTSIDE_OBCS = False ! [Boolean] default = False + ! If true, set the areas outside open boundaries to be land. + +! === module MOM_verticalGrid === +! Parameters providing information about the vertical grid. +!ANGSTROM = 1.0E-15 ! [m] default = 1.0E-10 + ! The minimum layer thickness, usually one-Angstrom. +NK = 2 ! [nondim] + ! The number of model layers. ! === module MOM_EOS === EQN_OF_STATE = "LINEAR" ! default = "WRIGHT" ! EQN_OF_STATE determines which ocean equation of state should be used. ! Currently, the valid choices are "LINEAR", "UNESCO", "WRIGHT", "NEMO" and ! "TEOS10". This is only used if USE_EOS is true. +RHO_T0_S0 = 1000.0 ! [kg m-3] default = 1000.0 + ! When EQN_OF_STATE=LINEAR, this is the density at T=0, S=0. +DRHO_DT = -0.2 ! [kg m-3 K-1] default = -0.2 + ! When EQN_OF_STATE=LINEAR, this is the partial derivative of density with + ! temperature. DRHO_DS = 0.0 ! [kg m-3 PSU-1] default = 0.8 ! When EQN_OF_STATE=LINEAR, this is the partial derivative of density with ! salinity. - -! === module MOM_restart === - +F_0 = 1.0E-4 ! [s-1] default = 0.0 + ! The reference value of the Coriolis parameter with the betaplane option. ! === module MOM_tracer_flow_control === ! === module MOM_coord_initialization === -COORD_CONFIG = "linear" ! +COORD_CONFIG = "linear" ! default = "none" ! This specifies how layers are to be defined: ! ALE or none - used to avoid defining layers in ALE mode ! file - read coordinate information from the file @@ -139,13 +158,19 @@ REGRIDDING_COORDINATE_MODE = "Z*" ! default = "LAYER" ! HYCOM1 - HyCOM-like hybrid coordinate ! SLIGHT - stretched coordinates above continuous isopycnal ! ADAPTIVE - optimize for smooth neutral density surfaces -!ALE_RESOLUTION = 2*50.0 ! [m] +REGRIDDING_COORDINATE_UNITS = "m" ! default = "m" + ! Units of the regridding coordinate. +!ALE_RESOLUTION = 2*1.0 ! [m] ! The distribution of vertical resolution for the target ! grid used for Eulerian-like coordinates. For example, ! in z-coordinate mode, the parameter is a list of level ! thicknesses (in m). In sigma-coordinate mode, the list ! is of non-dimensional fractions of the water column. -REMAPPING_SCHEME = "PPM_IH4" ! default = "PLM" +!TARGET_DENSITIES = 1000.0, 1001.0, 1002.0 ! [kg m^-3] + ! RHO target densities for interfaces +!MIN_THICKNESS = 1.0E-12 ! [m] default = 0.001 + ! When regridding, this is the minimum layer thickness allowed. +REMAPPING_SCHEME = "PPM_IH4" ! default = "PLM" ! This sets the reconstruction scheme used for vertical remapping for all ! variables. It can be one of the following schemes: PCM (1st-order ! accurate) @@ -154,6 +179,10 @@ REMAPPING_SCHEME = "PPM_IH4" ! default = "PLM" ! PPM_IH4 (3rd-order accurate) ! PQM_IH4IH3 (4th-order accurate) ! PQM_IH6IH5 (5th-order accurate) +REMAP_AFTER_INITIALIZATION = True ! [Boolean] default = True + ! If true, applies regridding and remapping immediately after initialization so + ! that the state is ALE consistent. This is a legacy step and should not be + ! needed if the initialization is consistent with the coordinate mode. ! === module MOM_grid === ! Parameters providing information about the lateral grid. @@ -168,11 +197,54 @@ TEMP_SALT_Z_INIT_FILE = "temp_salt_ic.nc" ! default = "temp_salt_z.nc" ! The name of the z-space input file used to initialize temperatures (T) and ! salinities (S). If T and S are not in the same file, TEMP_Z_INIT_FILE and ! SALT_Z_INIT_FILE must be set. +TEMP_Z_INIT_FILE = "temp_salt_ic.nc" ! default = "temp_salt_ic.nc" + ! The name of the z-space input file used to initialize temperatures, only. +SALT_Z_INIT_FILE = "temp_salt_ic.nc" ! default = "temp_salt_ic.nc" + ! The name of the z-space input file used to initialize temperatures, only. +Z_INIT_FILE_PTEMP_VAR = "ptemp" ! default = "ptemp" + ! The name of the potential temperature variable in TEMP_Z_INIT_FILE. +Z_INIT_FILE_SALT_VAR = "salt" ! default = "salt" + ! The name of the salinity variable in SALT_Z_INIT_FILE. +Z_INIT_HOMOGENIZE = False ! [Boolean] default = False + ! If True, then horizontally homogenize the interpolated initial conditions. Z_INIT_ALE_REMAPPING = True ! [Boolean] default = False ! If True, then remap straight to model coordinate from file. +Z_INIT_REMAPPING_SCHEME = "PPM_IH4" ! default = "PPM_IH4" + ! The remapping scheme to use if using Z_INIT_ALE_REMAPPING is True. +Z_INIT_REMAP_GENERAL = False ! [Boolean] default = False + ! If false, only initializes to z* coordinates. If true, allows initialization + ! directly to general coordinates. +TEMP_SALT_INIT_VERTICAL_REMAP_ONLY = False ! [Boolean] default = False + ! If true, initial conditions are on the model horizontal grid. Extrapolation + ! over missing ocean values is done using an ICE-9 procedure with vertical ALE + ! remapping . +TRIM_IC_FOR_P_SURF = False ! [Boolean] default = False + ! If true, cuts way the top of the column for initial conditions at the depth + ! where the hydrostatic pressure matches the imposed surface pressure which is + ! read from file. +SURFACE_PRESSURE_FILE = "shelf.nc" ! + ! The initial condition file for the surface pressure exerted by ice. +SURFACE_PRESSURE_VAR = "thick" ! [Pa] default = "" + ! The initial condition variable for the surface pressure exerted by ice. +SURFACE_PRESSURE_SCALE = 9800. ! [file dependent] default = 1.0 + ! A scaling factor to convert SURFACE_PRESSURE_VAR from file + ! SURFACE_PRESSURE_FILE into a surface pressure. +TRIMMING_USES_REMAPPING = False ! [Boolean] default = False + ! When trimming the column, also remap T and S. SPONGE = True ! [Boolean] default = False ! If true, sponges may be applied anywhere in the domain. The exact location and ! properties of those sponges are specified via SPONGE_CONFIG. +SPONGE_CONFIG = "file" ! default = "file" + ! A string that sets how the sponges are configured: + ! file - read sponge properties from the file + ! specified by (SPONGE_FILE). + ! ISOMIP - apply ale sponge in the ISOMIP case + ! RGC - apply sponge in the rotating_gravity_current case + ! DOME - use a slope and channel configuration for the + ! DOME sill-overflow test case. + ! BFB - Sponge at the southern boundary of the domain + ! for buoyancy-forced basin case. + ! USER - call a user modified routine. SPONGE_DAMPING_FILE = "sponge.nc" ! ! The name of the file with the sponge damping rates. SPONGE_STATE_FILE = "temp_salt_ic.nc" ! default = "sponge.nc" @@ -181,6 +253,10 @@ SPONGE_PTEMP_VAR = "ptemp" ! default = "PTEMP" ! The name of the potential temperature variable in SPONGE_STATE_FILE. SPONGE_SALT_VAR = "salt" ! default = "SALT" ! The name of the salinity variable in SPONGE_STATE_FILE. +SPONGE_ETA_VAR = "ETA" ! default = "ETA" + ! The name of the interface height variable in SPONGE_STATE_FILE. +SPONGE_IDAMP_VAR = "IDAMP" ! default = "IDAMP" + ! The name of the inverse damping rate variable in SPONGE_DAMPING_FILE. NEW_SPONGES = True ! [of sponge restoring data.] default = False ! Set True if using the newer sponging code which performs on-the-fly regridding ! in lat-lon-time. @@ -188,8 +264,6 @@ NEW_SPONGES = True ! [of sponge restoring data.] default = False ! === module MOM_sponge === SPONGE_DATA_ONGRID = True ! [Boolean] default = False ! When defined, the incoming sponge data are assumed to be on the model grid -!Total sponge columns at h points = 100 ! - ! The total number of columns where sponges are applied at h points. ! === module MOM_diag_mediator === @@ -256,14 +330,16 @@ BOUND_CORIOLIS = True ! [Boolean] default = False ! option is always effectively false with CORIOLIS_EN_DIS defined and ! CORIOLIS_SCHEME set to SADOURNY75_ENERGY. +DYNAMIC_SURFACE_PRESSURE = False ! [Boolean] default = False + ! If true, add a dynamic pressure due to a viscous ice shelf, for instance. + ! === module MOM_PressureForce === -! === module MOM_PressureForce_AFV === -RECONSTRUCT_FOR_PRESSURE = False ! [Boolean] default = True - ! If True, use vertical reconstruction of T & S within the integrals of the FV - ! pressure gradient calculation. If False, use the constant-by-layer algorithm. - ! The default is set by USE_REGRIDDING. + RECONSTRUCT_FOR_PRESSURE = False ! [Boolean] default = True + ! If True, use vertical reconstruction of T & S within the integrals of the FV + ! pressure gradient calculation. If False, use the constant-by-layer algorithm. + ! The default is set by USE_REGRIDDING. ! === module MOM_hor_visc === SMAGORINSKY_AH = True ! [Boolean] default = False ! If true, use a biharmonic Smagorinsky nonlinear eddy viscosity. @@ -274,7 +350,7 @@ SMAG_BI_CONST = 0.03 ! [nondim] default = 0.0 DIRECT_STRESS = True ! [Boolean] default = False ! If true, the wind stress is distributed over the topmost HMIX_STRESS of fluid ! (like in HYCOM), and KVML may be set to a very small value. -HMIX_FIXED = 20.0 ! [m] +HMIX_FIXED = 20. ! [m] ! The prescribed depth over which the near-surface viscosity and diffusivity are ! elevated when the bulk mixed layer is not used. KVML = 0.01 ! [m2 s-1] default = 1.0E-04 @@ -308,18 +384,12 @@ DTBT = 10.0 ! [s or nondim] default = -0.98 ! === module MOM_diabatic_driver === ! The following parameters are used for diabatic processes. - -! === module MOM_CVMix_KPP === -! This is the MOM wrapper to CVMix:KPP -! See http://cvmix.github.io/ - -! === module MOM_tidal_mixing === -! Vertical Tidal Mixing Parameterization - -! === module MOM_CVMix_conv === -! Parameterization of enhanced mixing due to convection via CVMix - -! === module MOM_entrain_diffusive === +USE_LEGACY_DIABATIC_DRIVER = True ! [Boolean] default = True + ! If true, use a legacy version of the diabatic subroutine. This is temporary + ! and is needed to avoid change in answers. +ENERGETICS_SFC_PBL = False ! [Boolean] default = False + ! If true, use an implied energetics planetary boundary layer scheme to + ! determine the diffusivity and viscosity in the surface boundary layer. ! === module MOM_set_diffusivity === BBL_EFFIC = 0.0 ! [nondim] default = 0.2 @@ -356,12 +426,19 @@ KD = 0.0 ! [m2 s-1] ! This module implements neutral diffusion of tracers ! === module MOM_sum_output === +CALCULATE_APE = False ! [Boolean] default = True + ! If true, calculate the available potential energy of the interfaces. Setting + ! this to false reduces the memory footprint of high-PE-count models + ! dramatically. MAXTRUNC = 5000 ! [truncations save_interval-1] default = 0 ! The run will be stopped, and the day set to a very large value if the velocity ! is truncated more than MAXTRUNC times between energy saves. Set MAXTRUNC to 0 ! to stop if there is any truncation of velocities. DATE_STAMPED_STDOUT = False ! [Boolean] default = True ! If true, use dates (not times) in messages to stdout +ENERGYSAVEDAYS = 0.125 ! [days] default = 1.0 + ! The interval in units of TIMEUNIT between saves of the energies of the run and + ! other globally summed diagnostics. ! === module MOM_surface_forcing === VARIABLE_WINDS = False ! [Boolean] default = True @@ -377,18 +454,107 @@ WIND_CONFIG = "zero" ! ! options include (file), (2gyre), (1gyre), (gyres), (zero), and (USER). ! === module MOM_restart === +ICE_SHELF = False ! [Boolean] default = False + ! If true, enables the ice shelf model. + +! === module MOM_domains min_halo === + +! === module MOM_grid === +! Parameters providing information about the lateral grid. + +! === module MOM_hor_index === +! Sets the horizontal array index types. +! === module MOM_grid_init === + +! === module MOM_ice_shelf === +DEBUG_IS = True ! [Boolean] default = False + ! If true, write verbose debugging messages for the ice shelf. +DYNAMIC_SHELF_MASS = False ! [Boolean] default = False + ! If true, the ice sheet mass can evolve with time. +SHELF_THERMO = True ! [Boolean] default = False + ! If true, use a thermodynamically interactive ice shelf. +SHELF_THREE_EQN = True ! [Boolean] default = True + ! If true, use the three equation expression of consistency to calculate the + ! fluxes at the ice-ocean interface. +SHELF_INSULATOR = False ! [Boolean] default = False + ! If true, the ice shelf is a perfect insulatior (no conduction). +MELTING_CUTOFF_DEPTH = 0.0 ! [m] default = 0.0 + ! Depth above which the melt is set to zero (it must be >= 0) Default value + ! won't affect the solution. +CONST_SEA_LEVEL = False ! [Boolean] default = False + ! If true, apply evaporative, heat and salt fluxes in the sponge region. This + ! will avoid a large increase in sea level. This option is needed for some of + ! the ISOMIP+ experiments (Ocean3 and Ocean4). IMPORTANT: it is not currently + ! possible to do prefect restarts using this flag. +SHELF_3EQ_GAMMA = True ! [Boolean] default = False + ! If true, user specifies a constant nondimensional heat-transfer coefficient + ! (GAMMA_T_3EQ), from which the default salt-transfer coefficient is set as + ! GAMMA_T_3EQ/35. This is used with SHELF_THREE_EQN. +SHELF_2EQ_GAMMA_T = 0.1 ! [m s-1] + ! If SHELF_THREE_EQN is false, this the fixed turbulent exchange velocity at the + ! ice-ocean interface. +SHELF_3EQ_GAMMA_T = 0.22 ! [nondim] default = 0.022 + ! Nondimensional heat-transfer coefficient. +SHELF_3EQ_GAMMA_S = 0.0062857142857143 ! [nondim] default = 0.0062857142857143 + ! Nondimensional salt-transfer coefficient. +ICE_SHELF_MASS_FROM_FILE = False ! [Boolean] default = False + ! Read the mass of the ice shelf (every time step) from a file. +C_P_ICE = 2100.0 ! [J kg-1 K-1] default = 2100.0 + ! The heat capacity of ice. +ICE_SHELF_FLUX_FACTOR = 1.0 ! [none] default = 1.0 + ! Non-dimensional factor applied to shelf thermodynamic fluxes. +KV_ICE = 1.0E+10 ! [m2 s-1] default = 1.0E+10 + ! The viscosity of the ice. +KV_MOLECULAR = 1.95E-06 ! [m2 s-1] default = 1.95E-06 + ! The molecular kinimatic viscosity of sea water at the freezing temperature. +ICE_SHELF_SALINITY = 0.0 ! [psu] default = 0.0 + ! The salinity of the ice inside the ice shelf. +ICE_SHELF_TEMPERATURE = -15.0 ! [degC] default = -15.0 + ! The temperature at the center of the ice shelf. +KD_SALT_MOLECULAR = 8.02E-10 ! [m2 s-1] default = 8.02E-10 + ! The molecular diffusivity of salt in sea water at the freezing point. +KD_TEMP_MOLECULAR = 1.41E-07 ! [m2 s-1] default = 1.41E-07 + ! The molecular diffusivity of heat in sea water at the freezing point. +COL_THICK_MELT_THRESHOLD = 0.0 ! [m] default = 0.0 + ! The minimum ocean column thickness where melting is allowed. +READ_TIDEAMP = False ! [Boolean] default = False + ! If true, read a file (given by TIDEAMP_FILE) containing the tidal amplitude + ! with INT_TIDE_DISSIPATION. +UTIDE = 0.0 ! [m s-1] default = 0.0 + ! The constant tidal amplitude used with INT_TIDE_DISSIPATION. + +! === module MOM_EOS === +DENSITY_ICE = 1000.0 ! [kg m-3] default = 900.0 + ! A typical density of ice. +MIN_THICKNESS_SIMPLE_CALVE = 0.0 ! [m] default = 0.0 + ! Min thickness rule for the very simple calving law +USTAR_SHELF_BG = 0.0 ! [m s-1] default = 0.0 + ! The minimum value of ustar under ice shelves. +CDRAG_SHELF = 0.003 ! [nondim] default = 0.003 + ! CDRAG is the drag coefficient relating the magnitude of the velocity field to + ! the surface stress. +DRAG_BG_VEL_SHELF = 0.01 ! [m s-1] default = 0.0 + ! DRAG_BG_VEL is either the assumed bottom velocity (with LINEAR_DRAG) or an + ! unresolved velocity that is combined with the resolved velocity to estimate + ! the velocity magnitude. + +! === module MOM_restart === +ICE_PROFILE_CONFIG = "FILE" ! + ! This specifies how the initial ice profile is specified. Valid values are: + ! CHANNEL, FILE, and USER. +LEN_SIDE_STRESS = 0.0 ! [axis_units] default = 0.0 + ! position past which shelf sides are stress free. +ICE_THICKNESS_VARNAME = "thick" ! default = "h_shelf" + ! The name of the thickness variable in ICE_THICKNESS_FILE. +SHELF_IC_OUTPUT_FILE = "MOM_Shelf_IC" ! default = "MOM_Shelf_IC" + ! The name-root of the output file for the ice shelf initial conditions. ! === module MOM_main (MOM_driver) === -DAYMAX = 0.25 ! [days] +DAYMAX = 0.25 ! [days] ! The final time of the whole simulation, in units of TIMEUNIT seconds. This ! also sets the potential end time of the present run segment if the end time is ! not set via ocean_solo_nml in input.nml. - -ENERGYSAVEDAYS = 0.125 ! [days] default = 1.44E+04 - ! The interval in units of TIMEUNIT between saves of the - ! energies of the run and other globally summed diagnostics. - -RESTART_CONTROL = 3 ! default = 1 +RESTART_CONTROL = 0 ! default = 1 ! An integer whose bits encode which restart files are written. Add 2 (bit 1) ! for a time-stamped file, and odd (bit 0) for a non-time-stamped file. A ! non-time-stamped restart file is saved at the end of the run segment for any @@ -402,25 +568,5 @@ MAXCPU = 2.88E+04 ! [wall-clock seconds] default = -1.0 ! not desired, use a negative value for MAXCPU. MAXCPU has units of wall-clock ! seconds, so the actual CPU time used is larger by a factor of the number of ! processors used. - -! === module MOM_file_parser === - DIAG_AS_CHKSUM = True DEBUG = True - -USE_PSURF_IN_EOS = False ! [Boolean] default = False -DEFAULT_2018_ANSWERS = True ! [Boolean] default = True -GRID_ROTATION_ANGLE_BUGS = True ! [Boolean] default = True -INTERPOLATE_RES_FN = True ! [Boolean] default = True -GILL_EQUATORIAL_LD = False ! [Boolean] default = False -USE_GM_WORK_BUG = True ! [Boolean] default = True -FIX_UNSPLIT_DT_VISC_BUG = False ! [Boolean] default = False -REMAP_UV_USING_OLD_ALG = True ! [Boolean] default = True -USE_LAND_MASK_FOR_HVISC = False ! [Boolean] default = False -KAPPA_SHEAR_ITER_BUG = True ! [Boolean] default = True -KAPPA_SHEAR_ALL_LAYER_TKE_BUG = True ! [Boolean] default = True -USE_MLD_ITERATION = False ! [Boolean] default = False -PEN_SW_ABSORB_MINTHICK = 0.001 ! [m] default = 0.001 -GUST_CONST = 0.02 ! [Pa] default = 0.02 -FIX_USTAR_GUSTLESS_BUG = False ! [Boolean] default = False - diff --git a/.testing/tc4/build_grid.py b/.testing/tc4/build_grid.py index 7f1be74efd..1223f5e635 100644 --- a/.testing/tc4/build_grid.py +++ b/.testing/tc4/build_grid.py @@ -9,7 +9,6 @@ topo_ = np.zeros((ny, nx)) + depth0 f_topo = nc.Dataset('topog.nc', 'w', format='NETCDF3_CLASSIC') -ny, nx = topo_.shape f_topo.createDimension('ny', ny) f_topo.createDimension('nx', nx) f_topo.createDimension('ntiles', 1) @@ -28,6 +27,7 @@ rad_deg = np.pi / 180. dx[:] = (rad_deg * Re * (x[:, 1:] - x[:, 0:-1]) * np.cos(0.5*rad_deg*(y[:, 0:-1] + y[:, 1:]))) +#dx[:] = (rad_deg * Re * (x[:, 1:] - x[:, 0:-1])) dy[:] = rad_deg * Re * (y[1:, :] - y[0:-1, :]) f_sg = nc.Dataset('ocean_hgrid.nc', 'w', format='NETCDF3_CLASSIC') @@ -74,3 +74,17 @@ f_sg.variables['tile'][:] = str_ f_sg.sync() f_sg.close() + +thick_ = np.zeros((ny, nx)) +area_ = area[::2,::2]+area[::2,1::2]+area[1::2,::2]+area[1::2,1::2] +thick_[:,:int(nx/2)]=0.5*depth0 +area_[thick_==0.]=0. +f_thick = nc.Dataset('shelf.nc', 'w', format='NETCDF3_CLASSIC') +f_thick.createDimension('ny', ny) +f_thick.createDimension('nx', nx) +f_thick.createVariable('thick', 'f8', ('ny', 'nx')) +f_thick.createVariable('area', 'f8', ('ny', 'nx')) +f_thick.variables['thick'][:] = thick_ +f_thick.variables['area'][:] = area_ +f_thick.sync() +f_thick.close() diff --git a/config_src/coupled_driver/ocean_model_MOM.F90 b/config_src/coupled_driver/ocean_model_MOM.F90 index 082099158c..ad83ceaf18 100644 --- a/config_src/coupled_driver/ocean_model_MOM.F90 +++ b/config_src/coupled_driver/ocean_model_MOM.F90 @@ -48,6 +48,7 @@ module ocean_model_mod use MOM_verticalGrid, only : verticalGrid_type use MOM_ice_shelf, only : initialize_ice_shelf, shelf_calc_flux, ice_shelf_CS use MOM_ice_shelf, only : add_shelf_forces, ice_shelf_end, ice_shelf_save_restart +use MOM_IS_diag_mediator, only : diag_IS_ctrl => diag_ctrl, diag_mediator_IS_end=>diag_mediator_end use coupler_types_mod, only : coupler_1d_bc_type, coupler_2d_bc_type use coupler_types_mod, only : coupler_type_spawn, coupler_type_write_chksums use coupler_types_mod, only : coupler_type_initialized, coupler_type_copy_data @@ -179,13 +180,13 @@ module ocean_model_mod !! processes before time stepping the dynamics. type(directories) :: dirs !< A structure containing several relevant directory paths. - type(mech_forcing) :: forces !< A structure with the driving mechanical surface forces - type(forcing) :: fluxes !< A structure containing pointers to - !! the thermodynamic ocean forcing fields. - type(forcing) :: flux_tmp !< A secondary structure containing pointers to the + type(mech_forcing), pointer :: forces => NULL() !< A structure with the driving mechanical surface forces + type(forcing), pointer :: fluxes => NULL() !< A structure containing pointers to + !! the thermodynamic ocean forcing fields. + type(forcing), pointer :: flux_tmp => NULL() !< A secondary structure containing pointers to the !! ocean forcing fields for when multiple coupled !! timesteps are taken per thermodynamic step. - type(surface) :: sfc_state !< A structure containing pointers to + type(surface), pointer :: sfc_state => NULL() !< A structure containing pointers to !! the ocean surface state fields. type(ocean_grid_type), pointer :: & grid => NULL() !< A pointer to a grid structure containing metrics @@ -213,7 +214,10 @@ module ocean_model_mod restart_CSp => NULL() !< A pointer set to the restart control structure !! that will be used for MOM restart files. type(diag_ctrl), pointer :: & - diag => NULL() !< A pointer to the diagnostic regulatory structure + diag => NULL() !< A pointer to the diagnostic regulatory structure + type(diag_IS_ctrl), pointer :: & + diag_IS => NULL() !< A pointer to the diagnostic regulatory structure + !! for the ice shelf module. end type ocean_state_type contains @@ -267,6 +271,9 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, wind_stagger, gas endif allocate(OS) + allocate(OS%fluxes) + allocate(OS%flux_tmp) + OS%is_ocean_pe = Ocean_sfc%is_ocean_pe if (.not.OS%is_ocean_pe) return @@ -355,6 +362,7 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, wind_stagger, gas use_melt_pot=.false. endif + allocate(OS%sfc_state) call allocate_surface_state(OS%sfc_state, OS%grid, use_temperature, do_integrals=.true., & gas_fields_ocn=gas_fields_ocn, use_meltpot=use_melt_pot) @@ -366,9 +374,11 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, wind_stagger, gas OS%forcing_CSp) endif + allocate(OS%forces) + if (OS%use_ice_shelf) then call initialize_ice_shelf(param_file, OS%grid, OS%Time, OS%ice_shelf_CSp, & - OS%diag, OS%forces, OS%fluxes) + OS%diag_IS, OS%forces, OS%fluxes) endif if (OS%icebergs_alter_ocean) then call marine_ice_init(OS%Time, OS%grid, param_file, OS%diag, OS%marine_ice_CSp) @@ -717,6 +727,8 @@ subroutine ocean_model_end(Ocean_sfc, Ocean_state, Time) call ocean_model_save_restart(Ocean_state, Time) call diag_mediator_end(Time, Ocean_state%diag) + if (Ocean_state%use_ice_shelf) & + call diag_mediator_IS_end(Time, Ocean_state%diag_IS) call MOM_end(Ocean_state%MOM_CSp) if (Ocean_state%use_ice_shelf) call ice_shelf_end(Ocean_state%Ice_shelf_CSp) end subroutine ocean_model_end diff --git a/config_src/solo_driver/MOM_driver.F90 b/config_src/solo_driver/MOM_driver.F90 index ba52d9c02a..0641dc248f 100644 --- a/config_src/solo_driver/MOM_driver.F90 +++ b/config_src/solo_driver/MOM_driver.F90 @@ -28,6 +28,7 @@ program MOM_main use MOM_cpu_clock, only : CLOCK_COMPONENT use MOM_diag_mediator, only : enable_averaging, disable_averaging, diag_mediator_end use MOM_diag_mediator, only : diag_ctrl, diag_mediator_close_registration + use MOM_IS_diag_mediator, only : diag_IS_ctrl=>diag_ctrl, diag_mediator_IS_end=>diag_mediator_end use MOM, only : initialize_MOM, step_MOM, MOM_control_struct, MOM_end use MOM, only : extract_surface_state, finish_MOM_initialization use MOM, only : get_MOM_state_elements, MOM_state_is_synchronized @@ -61,7 +62,7 @@ program MOM_main use MOM_verticalGrid, only : verticalGrid_type use MOM_write_cputime, only : write_cputime, MOM_write_cputime_init use MOM_write_cputime, only : write_cputime_start_clock, write_cputime_CS - + use MOM_get_input, only : get_MOM_input use ensemble_manager_mod, only : ensemble_manager_init, get_ensemble_size use ensemble_manager_mod, only : ensemble_pelist_setup use mpp_mod, only : set_current_pelist => mpp_set_current_pelist @@ -70,7 +71,6 @@ program MOM_main use MOM_ice_shelf, only : initialize_ice_shelf, ice_shelf_end, ice_shelf_CS use MOM_ice_shelf, only : shelf_calc_flux, add_shelf_forces, ice_shelf_save_restart -! , add_shelf_flux_forcing, add_shelf_flux_IOB use MOM_wave_interface, only: wave_parameters_CS, MOM_wave_interface_init use MOM_wave_interface, only: MOM_wave_interface_init_lite, Update_Surface_Waves @@ -80,22 +80,22 @@ program MOM_main #include ! A structure with the driving mechanical surface forces - type(mech_forcing) :: forces + type(mech_forcing), pointer :: forces => NULL() ! A structure containing pointers to the thermodynamic forcing fields ! at the ocean surface. - type(forcing) :: fluxes + type(forcing), pointer :: fluxes => NULL() ! A structure containing pointers to the ocean surface state fields. - type(surface) :: sfc_state + type(surface), pointer :: sfc_state => NULL() ! A pointer to a structure containing metrics and related information. - type(ocean_grid_type), pointer :: grid - type(verticalGrid_type), pointer :: GV + type(ocean_grid_type), pointer :: grid => NULL() + type(verticalGrid_type), pointer :: GV => NULL() ! A pointer to a structure containing dimensional unit scaling factors. - type(unit_scale_type), pointer :: US + type(unit_scale_type), pointer :: US => NULL() ! If .true., use the ice shelf model for part of the domain. - logical :: use_ice_shelf + logical :: use_ice_shelf = .false. ! If .true., use surface wave coupling logical :: use_waves = .false. @@ -198,8 +198,10 @@ program MOM_main type(MOM_restart_CS), pointer :: & restart_CSp => NULL() !< A pointer to the restart control structure !! that will be used for MOM restart files. - type(diag_ctrl), pointer :: & - diag => NULL() !< A pointer to the diagnostic regulatory structure + type(diag_ctrl), pointer :: & + diag => NULL() !< A pointer to the diagnostic regulatory structure + type(diag_IS_ctrl), pointer :: & + diag_IS => NULL() !< A pointer to the diagnostic regulatory structure !----------------------------------------------------------------------- character(len=4), parameter :: vers_num = 'v2.0' @@ -219,6 +221,8 @@ program MOM_main call MOM_infra_init() ; call io_infra_init() + allocate(forces,fluxes,sfc_state) + ! Initialize the ensemble manager. If there are no settings for ensemble_size ! in input.nml(ensemble.nml), these should not do anything. In coupled ! configurations, this all occurs in the external driver. @@ -297,16 +301,34 @@ program MOM_main ! In this case, the segment starts at a time fixed by ocean_solo.res segment_start_time = set_date(date(1),date(2),date(3),date(4),date(5),date(6)) Time = segment_start_time - call initialize_MOM(Time, Start_time, param_file, dirs, MOM_CSp, restart_CSp, & - segment_start_time, offline_tracer_mode=offline_tracer_mode, & - diag_ptr=diag, tracer_flow_CSp=tracer_flow_CSp) else ! In this case, the segment starts at a time read from the MOM restart file ! or left as Start_time by MOM_initialize. Time = Start_time + endif + + ! Read paths and filenames from namelist and store in "dirs". + ! Also open the parsed input parameter file(s) and setup param_file. + call get_MOM_input(param_file, dirs) + + call get_param(param_file, mod_name, "ICE_SHELF", use_ice_shelf, & + "If true, enables the ice shelf model.", default=.false.) + if (use_ice_shelf) then + ! These arrays are not initialized in most solo cases, but are needed + ! when using an ice shelf + call initialize_ice_shelf(param_file, grid, Time, ice_shelf_CSp, & + diag_IS, forces, fluxes, sfc_state) + endif + call close_param_file(param_file) + + if (sum(date) >= 0) then + call initialize_MOM(Time, Start_time, param_file, dirs, MOM_CSp, restart_CSp, & + segment_start_time, offline_tracer_mode=offline_tracer_mode, & + diag_ptr=diag, tracer_flow_CSp=tracer_flow_CSp,ice_shelf_CSp=ice_shelf_CSp) + else call initialize_MOM(Time, Start_time, param_file, dirs, MOM_CSp, restart_CSp, & offline_tracer_mode=offline_tracer_mode, diag_ptr=diag, & - tracer_flow_CSp=tracer_flow_CSp) + tracer_flow_CSp=tracer_flow_CSp,ice_shelf_CSp=ice_shelf_CSp) endif call get_MOM_state_elements(MOM_CSp, G=grid, GV=GV, US=US, C_p_scaled=fluxes%C_p) @@ -320,14 +342,6 @@ program MOM_main surface_forcing_CSp, tracer_flow_CSp) call callTree_waypoint("done surface_forcing_init") - call get_param(param_file, mod_name, "ICE_SHELF", use_ice_shelf, & - "If true, enables the ice shelf model.", default=.false.) - if (use_ice_shelf) then - ! These arrays are not initialized in most solo cases, but are needed - ! when using an ice shelf - call initialize_ice_shelf(param_file, grid, Time, ice_shelf_CSp, & - diag, forces, fluxes) - endif call get_param(param_file,mod_name,"USE_WAVES",Use_Waves,& "If true, enables surface wave modules.",default=.false.) @@ -434,6 +448,7 @@ program MOM_main call close_param_file(param_file) call diag_mediator_close_registration(diag) + ! Write out a time stamp file. if (calendar_type /= NO_CALENDAR) then call open_file(unit, 'time_stamp.out', form=ASCII_FILE, action=APPEND_FILE, & @@ -482,7 +497,7 @@ program MOM_main if (use_ice_shelf) then call shelf_calc_flux(sfc_state, fluxes, Time, dt_forcing, ice_shelf_CSp) - call add_shelf_forces(grid, US, Ice_shelf_CSp, forces) + call add_shelf_forces(grid, US, Ice_shelf_CSp, forces, external_call=.true.) endif fluxes%fluxes_used = .false. fluxes%dt_buoy_accum = US%s_to_T*dt_forcing @@ -652,6 +667,9 @@ program MOM_main call callTree_waypoint("End MOM_main") call diag_mediator_end(Time, diag, end_diag_manager=.true.) + if (use_ice_shelf) then + call diag_mediator_IS_end(Time, diag_IS) + endif if (cpu_steps > 0) call write_cputime(Time, ns-1, write_CPU_CSp, call_end=.true.) call cpu_clock_end(termClock) diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index f130c2977a..156a30b2d7 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -331,7 +331,7 @@ subroutine ALE_main( G, GV, US, h, u, v, tv, Reg, CS, OBC, dt, frac_shelf_h) type(ALE_CS), pointer :: CS !< Regridding parameters and options type(ocean_OBC_type), pointer :: OBC !< Open boundary structure real, optional, intent(in) :: dt !< Time step between calls to ALE_main [T ~> s] - real, dimension(:,:), optional, pointer :: frac_shelf_h !< Fractional ice shelf coverage + real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: frac_shelf_h !< Fractional ice shelf coverage [nondim] ! Local variables real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: dzRegrid ! The change in grid interface positions real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: eta_preale @@ -341,10 +341,7 @@ subroutine ALE_main( G, GV, US, h, u, v, tv, Reg, CS, OBC, dt, frac_shelf_h) nk = GV%ke; isc = G%isc; iec = G%iec; jsc = G%jsc; jec = G%jec - ice_shelf = .false. - if (present(frac_shelf_h)) then - if (associated(frac_shelf_h)) ice_shelf = .true. - endif + ice_shelf = present(frac_shelf_h) if (CS%show_call_tree) call callTree_enter("ALE_main(), MOM_ALE.F90") @@ -621,7 +618,7 @@ subroutine ALE_build_grid( G, GV, regridCS, remapCS, h, tv, debug, frac_shelf_h real, dimension(SZI_(G),SZJ_(G), SZK_(GV)), intent(inout) :: h !< Current 3D grid obtained after the !! last time step [H ~> m or kg-2] logical, optional, intent(in) :: debug !< If true, show the call tree - real, dimension(:,:), optional, pointer :: frac_shelf_h !< Fractional ice shelf coverage + real, dimension(SZI_(G),SZJ_(G)), optional, intent(in):: frac_shelf_h !< Fractional ice shelf coverage [nondim] ! Local variables integer :: nk, i, j, k real, dimension(SZI_(G), SZJ_(G), SZK_(GV)+1) :: dzRegrid ! The change in grid interface positions @@ -631,10 +628,7 @@ subroutine ALE_build_grid( G, GV, regridCS, remapCS, h, tv, debug, frac_shelf_h show_call_tree = .false. if (present(debug)) show_call_tree = debug if (show_call_tree) call callTree_enter("ALE_build_grid(), MOM_ALE.F90") - use_ice_shelf = .false. - if (present(frac_shelf_h)) then - if (associated(frac_shelf_h)) use_ice_shelf = .true. - endif + use_ice_shelf = present(frac_shelf_h) ! Build new grid. The new grid is stored in h_new. The old grid is h. ! Both are needed for the subsequent remapping of variables. diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 2a77cb06fe..c4bb4926c4 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -837,7 +837,7 @@ subroutine regridding_main( remapCS, CS, G, GV, h, tv, h_new, dzInterface, frac_ type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamical variables (T, S, ...) real, dimension(SZI_(G),SZJ_(G), CS%nk), intent(inout) :: h_new !< New 3D grid consistent with target coordinate real, dimension(SZI_(G),SZJ_(G), CS%nk+1), intent(inout) :: dzInterface !< The change in position of each interface - real, dimension(:,:), optional, pointer :: frac_shelf_h !< Fractional ice shelf coverage + real, dimension(SZI_(G),SZJ_(G)), optional, intent(in ) :: frac_shelf_h !< Fractional ice shelf coverage logical, optional, intent(in ) :: conv_adjust !< If true, do convective adjustment ! Local variables real :: trickGnuCompiler @@ -847,45 +847,31 @@ subroutine regridding_main( remapCS, CS, G, GV, h, tv, h_new, dzInterface, frac_ do_convective_adjustment = .true. if (present(conv_adjust)) do_convective_adjustment = conv_adjust - use_ice_shelf = .false. - if (present(frac_shelf_h)) then - if (associated(frac_shelf_h)) use_ice_shelf = .true. - endif + use_ice_shelf = present(frac_shelf_h) select case ( CS%regridding_scheme ) case ( REGRIDDING_ZSTAR ) - if (use_ice_shelf) then - call build_zstar_grid( CS, G, GV, h, dzInterface, frac_shelf_h ) - else - call build_zstar_grid( CS, G, GV, h, dzInterface ) - endif + call build_zstar_grid( CS, G, GV, h, dzInterface, frac_shelf_h ) call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new) - case ( REGRIDDING_SIGMA_SHELF_ZSTAR) - call build_zstar_grid( CS, G, GV, h, dzInterface ) + call build_zstar_grid( CS, G, GV, h, dzInterface, frac_shelf_h ) call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new) - case ( REGRIDDING_SIGMA ) call build_sigma_grid( CS, G, GV, h, dzInterface ) call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new) - case ( REGRIDDING_RHO ) if (do_convective_adjustment) call convective_adjustment(G, GV, h, tv) - call build_rho_grid( G, GV, G%US, h, tv, dzInterface, remapCS, CS ) + call build_rho_grid( G, GV, G%US, h, tv, dzInterface, remapCS, CS, frac_shelf_h ) call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new) - case ( REGRIDDING_ARBITRARY ) call build_grid_arbitrary( G, GV, h, dzInterface, trickGnuCompiler, CS ) call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new) - case ( REGRIDDING_HYCOM1 ) - call build_grid_HyCOM1( G, GV, G%US, h, tv, h_new, dzInterface, CS ) - + call build_grid_HyCOM1( G, GV, G%US, h, tv, h_new, dzInterface, CS, frac_shelf_h ) case ( REGRIDDING_SLIGHT ) call build_grid_SLight( G, GV, G%US, h, tv, dzInterface, CS ) call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new) - case ( REGRIDDING_ADAPTIVE ) call build_grid_adaptive(G, GV, G%US, h, tv, dzInterface, remapCS, CS) call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new) @@ -1164,21 +1150,24 @@ subroutine build_zstar_grid( CS, G, GV, h, dzInterface, frac_shelf_h) real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] real, dimension(SZI_(G),SZJ_(G), CS%nk+1), intent(inout) :: dzInterface !< The change in interface depth !! [H ~> m or kg m-2]. - real, dimension(:,:), optional, pointer :: frac_shelf_h !< Fractional ice shelf coverage [nondim]. + real, dimension(SZI_(G),SZJ_(G)), optional,intent(in) :: frac_shelf_h !< Fractional + !! ice shelf coverage [nondim]. ! Local variables - real :: nominalDepth, totalThickness, dh ! Depths and thicknesses [H ~> m or kg m-2] - real, dimension(SZK_(GV)+1) :: zOld, zNew ! Coordinate interface heights [H ~> m or kg m-2] - integer :: i, j, k, nz + integer :: i, j, k + integer :: nz + real :: nominalDepth, totalThickness, dh + real, dimension(SZK_(GV)+1) :: zOld, zNew + real :: minThickness logical :: ice_shelf nz = GV%ke - ice_shelf = .false. - if (present(frac_shelf_h)) then - if (associated(frac_shelf_h)) ice_shelf = .true. - endif + minThickness = CS%min_thickness + ice_shelf = present(frac_shelf_h) -!$OMP parallel do default(none) shared(G,GV,dzInterface,CS,nz,h,frac_shelf_h,ice_shelf) & -!$OMP private(nominalDepth,totalThickness,zNew,dh,zOld) +!$OMP parallel do default(none) shared(G,GV,dzInterface,CS,nz,h,frac_shelf_h, & +!$OMP ice_shelf,minThickness) & +!$OMP private(nominalDepth,totalThickness, & +!$OMP zNew,dh,zOld) do j = G%jsc-1,G%jec+1 do i = G%isc-1,G%iec+1 @@ -1324,7 +1313,7 @@ end subroutine build_sigma_grid ! Build grid based on target interface densities !------------------------------------------------------------------------------ !> This routine builds a new grid based on a given set of target interface densities. -subroutine build_rho_grid( G, GV, US, h, tv, dzInterface, remapCS, CS ) +subroutine build_rho_grid( G, GV, US, h, tv, dzInterface, remapCS, CS, frac_shelf_h ) !------------------------------------------------------------------------------ ! This routine builds a new grid based on a given set of target interface ! densities (these target densities are computed by taking the mean value @@ -1343,24 +1332,26 @@ subroutine build_rho_grid( G, GV, US, h, tv, dzInterface, remapCS, CS ) ! Arguments type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure real, dimension(SZI_(G),SZJ_(G), SZK_(GV)+1), intent(inout) :: dzInterface !< The change in interface depth !! [H ~> m or kg m-2] type(remapping_CS), intent(in) :: remapCS !< The remapping control structure type(regridding_CS), intent(in) :: CS !< Regridding control structure - + real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: frac_shelf_h !< Fractional + !! ice shelf coverage [nodim] ! Local variables integer :: nz integer :: i, j, k real :: nominalDepth ! Depth of the bottom of the ocean, positive downward [H ~> m or kg m-2] real, dimension(SZK_(GV)+1) :: zOld, zNew ! Old and new interface heights [H ~> m or kg m-2] real :: h_neglect, h_neglect_edge ! Negligible thicknesses [H ~> m or kg m-2] + real :: totalThickness ! Total thicknesses [H ~> m or kg m-2] #ifdef __DO_SAFETY_CHECKS__ - real :: totalThickness real :: dh #endif + logical :: ice_shelf if (.not.CS%remap_answers_2018) then h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff @@ -1371,6 +1362,7 @@ subroutine build_rho_grid( G, GV, US, h, tv, dzInterface, remapCS, CS ) endif nz = GV%ke + ice_shelf = present(frac_shelf_h) if (.not.CS%target_density_set) call MOM_error(FATAL, "build_rho_grid: "//& "Target densities must be set before build_rho_grid is called.") @@ -1388,9 +1380,27 @@ subroutine build_rho_grid( G, GV, US, h, tv, dzInterface, remapCS, CS ) ! Local depth (G%bathyT is positive) nominalDepth = G%bathyT(i,j)*GV%Z_to_H - call build_rho_column(CS%rho_CS, nz, nominalDepth, h(i, j, :), & + ! Determine total water column thickness + totalThickness = 0.0 + do k=1,nz + totalThickness = totalThickness + h(i,j,k) + enddo + ! Determine absolute interface positions + zOld(nz+1) = - nominalDepth + do k = nz,1,-1 + zOld(k) = zOld(k+1) + h(i,j,k) + enddo + + if (ice_shelf) then + call build_rho_column(CS%rho_CS, nz, nominalDepth, h(i, j, :), & + tv%T(i, j, :), tv%S(i, j, :), tv%eqn_of_state, zNew, & + z_rigid_top = totalThickness - nominalDepth, eta_orig = zOld(1), & + h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) + else + call build_rho_column(CS%rho_CS, nz, nominalDepth, h(i, j, :), & tv%T(i, j, :), tv%S(i, j, :), tv%eqn_of_state, zNew, & h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) + endif if (CS%integrate_downward_for_e) then zOld(1) = 0. @@ -1457,7 +1467,7 @@ end subroutine build_rho_grid !! \remark { Based on Bleck, 2002: An oceanice general circulation model framed in !! hybrid isopycnic-Cartesian coordinates, Ocean Modelling 37, 55-88. !! http://dx.doi.org/10.1016/S1463-5003(01)00012-9 } -subroutine build_grid_HyCOM1( G, GV, US, h, tv, h_new, dzInterface, CS ) +subroutine build_grid_HyCOM1( G, GV, US, h, tv, h_new, dzInterface, CS, frac_shelf_h ) type(ocean_grid_type), intent(in) :: G !< Grid structure type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -1466,16 +1476,19 @@ subroutine build_grid_HyCOM1( G, GV, US, h, tv, h_new, dzInterface, CS ) type(regridding_CS), intent(in) :: CS !< Regridding control structure real, dimension(SZI_(G),SZJ_(G),CS%nk), intent(inout) :: h_new !< New layer thicknesses [H ~> m or kg m-2] real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< Changes in interface position + real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: frac_shelf_h !< Fractional + !! ice shelf coverage [nodim] ! Local variables real, dimension(SZK_(GV)+1) :: z_col ! Source interface positions relative to the surface [H ~> m or kg m-2] real, dimension(CS%nk+1) :: z_col_new ! New interface positions relative to the surface [H ~> m or kg m-2] real, dimension(SZK_(GV)+1) :: dz_col ! The realized change in z_col [H ~> m or kg m-2] - real, dimension(SZK_(GV)) :: p_col ! Layer center pressure [R L2 T-2 ~> Pa] - real :: ref_pres ! The reference pressure [R L2 T-2 ~> Pa] + real, dimension(SZK_(GV)) :: p_col ! Layer center pressure [Pa] integer :: i, j, k, nki - real :: depth + real :: depth, nominalDepth real :: h_neglect, h_neglect_edge + real :: z_top_col, totalThickness + logical :: ice_shelf if (.not.CS%remap_answers_2018) then h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff @@ -1489,24 +1502,35 @@ subroutine build_grid_HyCOM1( G, GV, US, h, tv, h_new, dzInterface, CS ) "Target densities must be set before build_grid_HyCOM1 is called.") nki = min(GV%ke, CS%nk) + ice_shelf = present(frac_shelf_h) ! Build grid based on target interface densities do j = G%jsc-1,G%jec+1 ; do i = G%isc-1,G%iec+1 if (G%mask2dT(i,j)>0.) then - depth = G%bathyT(i,j) * GV%Z_to_H + nominalDepth = G%bathyT(i,j) * GV%Z_to_H - z_col(1) = 0. ! Work downward rather than bottom up + if (ice_shelf) then + totalThickness = 0.0 + do k=1,GV%ke + totalThickness = totalThickness + h(i,j,k) * GV%Z_to_H + enddo + z_top_col = max(nominalDepth-totalThickness,0.0) + else + z_top_col = 0.0 + endif + + z_col(1) = z_top_col ! Work downward rather than bottom up do K = 1, GV%ke z_col(K+1) = z_col(K) + h(i,j,k) p_col(k) = tv%P_Ref + CS%compressibility_fraction * & ( 0.5 * ( z_col(K) + z_col(K+1) ) * (GV%H_to_RZ*GV%g_Earth) - tv%P_Ref ) enddo - call build_hycom1_column(CS%hycom_CS, tv%eqn_of_state, GV%ke, depth, & - h(i,j,:), tv%T(i,j,:), tv%S(i,j,:), p_col, & - z_col, z_col_new, zScale=GV%Z_to_H, & - h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) + call build_hycom1_column(CS%hycom_CS, tv%eqn_of_state, GV%ke, nominalDepth, & + h(i,j,:), tv%T(i,j,:), tv%S(i,j,:), p_col, & + z_col, z_col_new, zScale=GV%Z_to_H, & + h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) ! Calculate the final change in grid position after blending new and old grids call filtered_grid_motion( CS, GV%ke, z_col, z_col_new, dz_col ) diff --git a/src/ALE/coord_rho.F90 b/src/ALE/coord_rho.F90 index dce802ff3c..c1e35ac314 100644 --- a/src/ALE/coord_rho.F90 +++ b/src/ALE/coord_rho.F90 @@ -88,7 +88,7 @@ end subroutine set_rho_params !! 1. Density profiles are calculated on the source grid. !! 2. Positions of target densities (for interfaces) are found by interpolation. subroutine build_rho_column(CS, nz, depth, h, T, S, eqn_of_state, z_interface, & - h_neglect, h_neglect_edge) + z_rigid_top, eta_orig, h_neglect, h_neglect_edge) type(rho_CS), intent(in) :: CS !< coord_rho control structure integer, intent(in) :: nz !< Number of levels on source grid (i.e. length of h, T, S) real, intent(in) :: depth !< Depth of ocean bottom (positive downward) [H ~> m or kg m-2] @@ -97,7 +97,11 @@ subroutine build_rho_column(CS, nz, depth, h, T, S, eqn_of_state, z_interface, & real, dimension(nz), intent(in) :: S !< Salinity for source column [ppt] type(EOS_type), pointer :: eqn_of_state !< Equation of state structure real, dimension(CS%nk+1), & - intent(inout) :: z_interface !< Absolute positions of interfaces + intent(inout) :: z_interface !< Absolute positions of interfaces + real, optional, intent(in) :: z_rigid_top !< The height of a rigid top (positive upward in the same + !! units as depth) [Z ~> m] or [H ~> m or kg m-2] + real, optional, intent(in) :: eta_orig !< The actual original height of the top in the same + !! units as depth) [Z ~> m] or [H ~> m or kg m-2] real, optional, intent(in) :: h_neglect !< A negligibly small width for the purpose !! of cell reconstructions [H ~> m or kg m-2] real, optional, intent(in) :: h_neglect_edge !< A negligibly small width for the purpose @@ -112,10 +116,22 @@ subroutine build_rho_column(CS, nz, depth, h, T, S, eqn_of_state, z_interface, & real, dimension(nz+1) :: xTmp ! Temporary positions [H ~> m or kg m-2] real, dimension(CS%nk) :: h_new ! New thicknesses [H ~> m or kg m-2] real, dimension(CS%nk+1) :: x1 ! Interface heights [H ~> m or kg m-2] + real :: z0_top, eta ! Thicknesses or heights [Z ~> m] or [H ~> m or kg m-2] ! Construct source column with vanished layers removed (stored in h_nv) call copy_finite_thicknesses(nz, h, CS%min_thickness, count_nonzero_layers, h_nv, mapping) + z0_top = 0. + eta=0.0 + if (present(z_rigid_top)) then + z0_top = z_rigid_top + eta=z0_top + if (present(eta_orig)) then + eta=eta_orig + endif + endif + + if (count_nonzero_layers > 1) then xTmp(1) = 0.0 do k = 1,count_nonzero_layers diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 0e736e7312..c6e016a65a 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -141,6 +141,7 @@ module MOM use MOM_offline_main, only : offline_fw_fluxes_into_ocean, offline_fw_fluxes_out_ocean use MOM_offline_main, only : offline_advection_layer, offline_transport_end use MOM_ALE, only : ale_offline_tracer_final, ALE_main_offline +use MOM_ice_shelf, only : ice_shelf_CS, ice_shelf_query implicit none ; private @@ -259,8 +260,8 @@ module MOM !! calculated, and if it is 0, dtbt is calculated every step. type(time_type) :: dtbt_reset_interval !< A time_time representation of dtbt_reset_period. type(time_type) :: dtbt_reset_time !< The next time DTBT should be calculated. - - + real, dimension(:,:), pointer :: frac_shelf_h => NULL() !< fraction of total area occupied + !! by ice shelf [nondim] real, dimension(:,:,:), pointer :: & h_pre_dyn => NULL(), & !< The thickness before the transports [H ~> m or kg m-2]. T_pre_dyn => NULL(), & !< Temperature before the transports [degC]. @@ -574,6 +575,8 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS !---------- Initiate group halo pass of the forcing fields call cpu_clock_begin(id_clock_pass) + if (.not.associated(forces%taux) .or. .not.associated(forces%tauy)) & + call MOM_error(FATAL,'step_MOM:forces%taux,tauy not associated') call create_group_pass(pass_tau_ustar_psurf, forces%taux, forces%tauy, G%Domain) if (associated(forces%ustar)) & call create_group_pass(pass_tau_ustar_psurf, forces%ustar, G%Domain) @@ -927,10 +930,9 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS ! De-rotate fluxes and copy back to the input, since they can be changed. if (CS%rotate_index) then call rotate_forcing(fluxes, fluxes_in, -turns) - + call rotate_mech_forcing(forces, -turns, forces_in) call deallocate_mech_forcing(forces) deallocate(forces) - call deallocate_forcing_type(fluxes) deallocate(fluxes) endif @@ -1242,7 +1244,7 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, & if (showCallTree) call callTree_enter("step_MOM_thermo(), MOM.F90") use_ice_shelf = .false. - if (associated(fluxes%frac_shelf_h)) use_ice_shelf = .true. + if (associated(CS%frac_shelf_h)) use_ice_shelf = .true. call enable_averages(dtdia, Time_end_thermo, CS%diag) @@ -1318,7 +1320,7 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, & call cpu_clock_begin(id_clock_ALE) if (use_ice_shelf) then call ALE_main(G, GV, US, h, u, v, tv, CS%tracer_Reg, CS%ALE_CSp, CS%OBC, & - dtdia, fluxes%frac_shelf_h) + dtdia, CS%frac_shelf_h) else call ALE_main(G, GV, US, h, u, v, tv, CS%tracer_Reg, CS%ALE_CSp, CS%OBC, dtdia) endif @@ -1600,7 +1602,7 @@ end subroutine step_offline !! initializing the ocean state variables, and initializing subsidiary modules subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & Time_in, offline_tracer_mode, input_restart_file, diag_ptr, & - count_calls, tracer_flow_CSp) + count_calls, tracer_flow_CSp, ice_shelf_CSp) type(time_type), target, intent(inout) :: Time !< model time, set in this routine type(time_type), intent(in) :: Time_init !< The start time for the coupled model's calendar type(param_file_type), intent(out) :: param_file !< structure indicating parameter file to parse @@ -1621,6 +1623,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & logical, optional, intent(in) :: count_calls !< If true, nstep_tot counts the number of !! calls to step_MOM instead of the number of !! dynamics timesteps. + type(ice_shelf_CS), optional, pointer :: ice_shelf_CSp !< A pointer to an ice shelf control structure ! local variables type(ocean_grid_type), pointer :: G => NULL() ! A pointer to the metric grid use for the run type(ocean_grid_type), pointer :: G_in => NULL() ! Pointer to the input grid @@ -1636,6 +1639,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & ! Initial state on the input index map real, allocatable, dimension(:,:,:) :: u_in, v_in, h_in + real, allocatable, dimension(:,:), target :: frac_shelf_in real, allocatable, dimension(:,:,:), target :: T_in, S_in type(ocean_OBC_type), pointer :: OBC_in => NULL() type(sponge_CS), pointer :: sponge_in_CSp => NULL() @@ -1649,9 +1653,8 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & real :: dtbt ! The barotropic timestep [s] real, allocatable, dimension(:,:) :: eta ! free surface height or column mass [H ~> m or kg m-2] - real, allocatable, dimension(:,:) :: area_shelf_h ! area occupied by ice shelf [L2 ~> m2] - real, dimension(:,:), allocatable, target :: frac_shelf_h ! fraction of total area occupied by ice shelf [nondim] - real, dimension(:,:), pointer :: shelf_area => NULL() + real, allocatable, dimension(:,:) :: area_shelf_in ! area occupied by ice shelf [L2 ~> m2] +! real, dimension(:,:), pointer :: shelf_area => NULL() type(MOM_restart_CS), pointer :: restart_CSp_tmp => NULL() type(group_pass_type) :: tmp_pass_uv_T_S_h, pass_uv_T_S_h @@ -2032,19 +2035,11 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & "initialize_MOM: A bulk mixed layer can only be used with T & S as "//& "state variables. Add USE_EOS = True to MOM_input.") - call get_param(param_file, 'MOM', "ICE_SHELF", use_ice_shelf, default=.false., do_not_log=.true.) - if (use_ice_shelf) then - inputdir = "." ; call get_param(param_file, 'MOM', "INPUTDIR", inputdir) - inputdir = slasher(inputdir) - call get_param(param_file, 'MOM', "ICE_THICKNESS_FILE", ice_shelf_file, & - "The file from which the ice bathymetry and area are read.", & - fail_if_missing=.true.) - call get_param(param_file, 'MOM', "ICE_AREA_VARNAME", area_varname, & - "The name of the area variable in ICE_THICKNESS_FILE.", & - fail_if_missing=.true.) + use_ice_shelf=.false. + if (present(ice_shelf_CSp)) then + if (associated(ice_shelf_CSp)) use_ice_shelf=.true. endif - CS%ensemble_ocean=.false. call get_param(param_file, "MOM", "ENSEMBLE_OCEAN", CS%ensemble_ocean, & "If False, The model is being run in serial mode as a single realization. "//& @@ -2348,10 +2343,14 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & call copy_dyngrid_to_MOM_grid(dG, G, US) call destroy_dyn_horgrid(dG) endif + call MOM_grid_init(G_in, param_file, US, HI_in, bathymetry_at_vel=bathy_at_vel) call copy_dyngrid_to_MOM_grid(dG_in, G_in, US) call destroy_dyn_horgrid(dG_in) + if (.not. CS%rotate_index) & + G => G_in + ! Set a few remaining fields that are specific to the ocean grid type. call set_first_direction(G, first_direction) ! Allocate the auxiliary non-symmetric domain for debugging or I/O purposes. @@ -2368,6 +2367,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & allocate(u_in(G_in%IsdB:G_in%IedB, G_in%jsd:G_in%jed, nz)) allocate(v_in(G_in%isd:G_in%ied, G_in%JsdB:G_in%JedB, nz)) allocate(h_in(G_in%isd:G_in%ied, G_in%jsd:G_in%jed, nz)) + u_in(:,:,:) = 0.0 v_in(:,:,:) = 0.0 h_in(:,:,:) = GV%Angstrom_H @@ -2382,9 +2382,28 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & CS%tv%S => S_in endif - call MOM_initialize_state(u_in, v_in, h_in, CS%tv, Time, G_in, GV, US, & - param_file, dirs, restart_CSp, CS%ALE_CSp, CS%tracer_Reg, & - sponge_in_CSp, ALE_sponge_in_CSp, OBC_in, Time_in) + if (use_ice_shelf) then + allocate(frac_shelf_in(G_in%isd:G_in%ied, G_in%jsd:G_in%jed)) + frac_shelf_in(:,:) = 0.0 + allocate(CS%frac_shelf_h(isd:ied, jsd:jed)) + CS%frac_shelf_h(:,:) = 0.0 + call ice_shelf_query(ice_shelf_CSp,G,CS%frac_shelf_h) + ! MOM_initialize_state is using the unrotated metric + call rotate_array(CS%frac_shelf_h, -turns, frac_shelf_in) + ! TODO: Verify that pass_var is needed here? + call pass_var(frac_shelf_in, G_in%Domain) + endif + + if (use_ice_shelf) then + call MOM_initialize_state(u_in, v_in, h_in, CS%tv, Time, G_in, GV, US, & + param_file, dirs, restart_CSp, CS%ALE_CSp, CS%tracer_Reg, & + sponge_in_CSp, ALE_sponge_in_CSp, OBC_in, Time_in, & + frac_shelf_h=frac_shelf_in) + else + call MOM_initialize_state(u_in, v_in, h_in, CS%tv, Time, G_in, GV, US, & + param_file, dirs, restart_CSp, CS%ALE_CSp, CS%tracer_Reg, & + sponge_in_CSp, ALE_sponge_in_CSp, OBC_in, Time_in) + endif if (use_temperature) then CS%tv%T => CS%T @@ -2415,12 +2434,29 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & deallocate(T_in) deallocate(S_in) endif + if (use_ice_shelf) & + deallocate(frac_shelf_in) else - call MOM_initialize_state(CS%u, CS%v, CS%h, CS%tv, Time, G, GV, US, & - param_file, dirs, restart_CSp, CS%ALE_CSp, CS%tracer_Reg, & - CS%sponge_CSp, CS%ALE_sponge_CSp, CS%OBC, Time_in) + if (use_ice_shelf) then + allocate(CS%frac_shelf_h(isd:ied, jsd:jed)) + CS%frac_shelf_h(:,:) = 0.0 + call ice_shelf_query(ice_shelf_CSp,G,CS%frac_shelf_h) + + call MOM_initialize_state(CS%u, CS%v, CS%h, CS%tv, Time, G, GV, US, & + param_file, dirs, restart_CSp, CS%ALE_CSp, CS%tracer_Reg, & + CS%sponge_CSp, CS%ALE_sponge_CSp, CS%OBC, Time_in, & + frac_shelf_h=CS%frac_shelf_h) + else + call MOM_initialize_state(CS%u, CS%v, CS%h, CS%tv, Time, G, GV, US, & + param_file, dirs, restart_CSp, CS%ALE_CSp, CS%tracer_Reg, & + CS%sponge_CSp, CS%ALE_sponge_CSp, CS%OBC, Time_in) + endif endif + if (use_ice_shelf) & + call hchksum(CS%frac_shelf_h, "MOM:frac_shelf_h", G%HI, & + haloshift=0) + call cpu_clock_end(id_clock_MOM_init) call callTree_waypoint("returned from MOM_initialize_state() (initialize_MOM)") @@ -2472,24 +2508,8 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & call adjustGridForIntegrity(CS%ALE_CSp, G, GV, CS%h ) call callTree_waypoint("Calling ALE_main() to remap initial conditions (initialize_MOM)") if (use_ice_shelf) then - filename = trim(inputdir)//trim(ice_shelf_file) - if (.not.file_exists(filename, G%Domain)) call MOM_error(FATAL, & - "MOM: Unable to open "//trim(filename)) - - allocate(area_shelf_h(isd:ied,jsd:jed)) - allocate(frac_shelf_h(isd:ied,jsd:jed)) - call MOM_read_data(filename, trim(area_varname), area_shelf_h, G%Domain, scale=US%m_to_L**2) - ! initialize frac_shelf_h with zeros (open water everywhere) - frac_shelf_h(:,:) = 0.0 - ! compute fractional ice shelf coverage of h - do j=jsd,jed ; do i=isd,ied - if (G%areaT(i,j) > 0.0) & - frac_shelf_h(i,j) = area_shelf_h(i,j) / G%areaT(i,j) - enddo ; enddo - ! pass to the pointer - shelf_area => frac_shelf_h call ALE_main(G, GV, US, CS%h, CS%u, CS%v, CS%tv, CS%tracer_Reg, CS%ALE_CSp, & - CS%OBC, frac_shelf_h=shelf_area) + CS%OBC, frac_shelf_h=CS%frac_shelf_h) else call ALE_main( G, GV, US, CS%h, CS%u, CS%v, CS%tv, CS%tracer_Reg, CS%ALE_CSp, CS%OBC) endif @@ -3018,8 +3038,7 @@ subroutine extract_surface_state(CS, sfc_state_in) type(verticalGrid_type), pointer :: GV => NULL() !< structure containing vertical grid info type(unit_scale_type), pointer :: US => NULL() !< structure containing various unit conversion factors type(surface), pointer :: sfc_state => NULL() ! surface state on the model grid - real, dimension(:,:,:), pointer :: & - h => NULL() !< h : layer thickness [H ~> m or kg m-2] + real, dimension(:,:,:), pointer :: h => NULL() !< h : layer thickness [H ~> m or kg m-2] real :: depth(SZI_(CS%G)) !< Distance from the surface in depth units [Z ~> m] or [H ~> m or kg m-2] real :: depth_ml !< Depth over which to average to determine mixed !! layer properties [Z ~> m] or [H ~> m or kg m-2] @@ -3037,6 +3056,7 @@ subroutine extract_surface_state(CS, sfc_state_in) integer :: isd, ied, jsd, jed integer :: iscB, iecB, jscB, jecB, isdB, iedB, jsdB, jedB logical :: localError + logical :: use_iceshelves character(240) :: msg integer :: turns ! Number of quarter turns @@ -3050,20 +3070,26 @@ subroutine extract_surface_state(CS, sfc_state_in) use_temperature = associated(CS%tv%T) + use_iceshelves=.false. + if (associated(CS%frac_shelf_h)) use_iceshelves = .true. + turns = 0 if (CS%rotate_index) & turns = G%HI%turns - if (.not.sfc_state_in%arrays_allocated) & + if (.not.sfc_state_in%arrays_allocated) then ! Consider using a run-time flag to determine whether to do the vertical ! integrals, since the 3-d sums are not negligible in cost. call allocate_surface_state(sfc_state_in, G_in, use_temperature, & - do_integrals=.true., omit_frazil=.not.associated(CS%tv%frazil)) + do_integrals=.true., omit_frazil=.not.associated(CS%tv%frazil),& + use_iceshelves=use_iceshelves) + endif if (CS%rotate_index) then allocate(sfc_state) call allocate_surface_state(sfc_state, G, use_temperature, & - do_integrals=.true., omit_frazil=.not.associated(CS%tv%frazil)) + do_integrals=.true., omit_frazil=.not.associated(CS%tv%frazil),& + use_iceshelves=use_iceshelves) else sfc_state => sfc_state_in endif diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 5e42a9575f..47ac680ae9 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -2173,7 +2173,8 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, scale=US%L_T_to_m_s*US%s_to_T) call uvchksum(trim(mesg)//" BT_force_[uv]", BT_force_u, BT_force_v, CS%debug_BT_HI, haloshift=iev-ie, & scale=US%L_T_to_m_s*US%s_to_T) - call uvchksum(trim(mesg)//" BT_rem_[uv]", BT_rem_u, BT_rem_v, CS%debug_BT_HI, haloshift=iev-ie) + call uvchksum(trim(mesg)//" BT_rem_[uv]", BT_rem_u, BT_rem_v, CS%debug_BT_HI, & + haloshift=iev-ie, scalar_pair=.true.) call uvchksum(trim(mesg)//" [uv]bt", ubt, vbt, CS%debug_BT_HI, haloshift=iev-ie, & scale=US%L_T_to_m_s) call uvchksum(trim(mesg)//" [uv]bt_trans", ubt_trans, vbt_trans, CS%debug_BT_HI, haloshift=iev-ie, & diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index ed4b8d1ba2..d11c1ae985 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -1135,8 +1135,9 @@ subroutine MOM_mech_forcing_chksum(mesg, forces, G, US, haloshift) if (associated(forces%ustar)) & call hchksum(forces%ustar, mesg//" forces%ustar", G%HI, haloshift=hshift, scale=US%Z_to_m*US%s_to_T) if (associated(forces%rigidity_ice_u) .and. associated(forces%rigidity_ice_v)) & - call uvchksum(mesg//" forces%rigidity_ice_[uv]", forces%rigidity_ice_u, forces%rigidity_ice_v, & - G%HI, haloshift=hshift, symmetric=.true., scale=US%L_to_m**3*US%L_to_Z*US%s_to_T) + call uvchksum(mesg//" forces%rigidity_ice_[uv]", forces%rigidity_ice_u, & + forces%rigidity_ice_v, G%HI, haloshift=hshift, symmetric=.true., & + scale=US%L_to_m**3*US%L_to_Z*US%s_to_T, scalar_pair=.true.) end subroutine MOM_mech_forcing_chksum diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index 0b225f0bf7..7267edcd9d 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -457,7 +457,7 @@ subroutine rotate_surface_state(sfc_state_in, G_in, sfc_state, G, turns) if (use_temperature) then call rotate_array(sfc_state_in%ocean_heat, turns, sfc_state%ocean_heat) call rotate_array(sfc_state_in%ocean_salt, turns, sfc_state%ocean_salt) - call rotate_array(sfc_state_in%SSS, turns, sfc_state%TempxPmE) + call rotate_array(sfc_state_in%SSS, turns, sfc_state%SSS) call rotate_array(sfc_state_in%salt_deficit, turns, sfc_state%salt_deficit) call rotate_array(sfc_state_in%internal_heat, turns, sfc_state%internal_heat) endif @@ -468,9 +468,9 @@ subroutine rotate_surface_state(sfc_state_in, G_in, sfc_state, G, turns) sfc_state%taux_shelf, sfc_state%tauy_shelf) endif - if (use_temperature .and. allocated(sfc_state_in%frazil)) & + if (use_temperature .and. allocated(sfc_state_in%frazil)) then call rotate_array(sfc_state_in%frazil, turns, sfc_state%frazil) - + endif ! Scalar transfers sfc_state%T_is_conT = sfc_state_in%T_is_conT sfc_state%S_is_absS = sfc_state_in%S_is_absS diff --git a/src/framework/MOM_horizontal_regridding.F90 b/src/framework/MOM_horizontal_regridding.F90 index 66f58b5b9d..5ff50f7aed 100644 --- a/src/framework/MOM_horizontal_regridding.F90 +++ b/src/framework/MOM_horizontal_regridding.F90 @@ -732,6 +732,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t allocate(z_in(kd),z_edges_in(kd+1)) allocate(tr_z(isd:ied,jsd:jed,kd), mask_z(isd:ied,jsd:jed,kd)) + tr_z(:,:,:)=0.0;mask_z(:,:,:)=0.0 call mpp_get_axis_data(axes_data(3), z_in) @@ -770,7 +771,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t allocate(mask_in(id,jdp)) ; mask_in(:,:)=0.0 allocate(last_row(id)) ; last_row(:)=0.0 else - allocate(data_in(isd:ied,jsd:jed,kd)) + allocate(data_in(isd:ied,jsd:jed,kd)); data_in(:,:,:)=0.0 endif ! construct level cell boundaries as the mid-point between adjacent centers z_edges_in(1) = 0.0 @@ -910,15 +911,16 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t enddo ! kd else - call time_interp_external(fms_id, Time, data_in, verbose=.true., turns=turns) - do k=1,kd - do j=js,je - do i=is,ie - tr_z(i,j,k)=data_in(i,j,k) - if (abs(tr_z(i,j,k)-missing_value) < abs(roundoff*missing_value)) mask_z(i,j,k) = 0. - enddo + call time_interp_external(fms_id, Time, data_in, verbose=.true., turns=turns) + do k=1,kd + do j=js,je + do i=is,ie + tr_z(i,j,k)=data_in(i,j,k) + mask_z(i,j,k)=1.0 + if (abs(tr_z(i,j,k)-missing_value) < abs(roundoff*missing_value)) mask_z(i,j,k) = 0. enddo enddo + enddo endif end subroutine horiz_interp_and_extrap_tracer_fms_id diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 66fd873f67..7a223efd9e 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -4,13 +4,16 @@ module MOM_ice_shelf ! This file is part of MOM6. See LICENSE.md for the license. - +use MOM_array_transform, only : rotate_array use MOM_constants, only : hlf use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end use MOM_cpu_clock, only : CLOCK_COMPONENT, CLOCK_ROUTINE -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, diag_ctrl, time_type -use MOM_diag_mediator, only : enable_averages, enable_averaging, disable_averaging +use MOM_coms, only : num_PEs +use MOM_IS_diag_mediator, only : post_data, register_diag_field=>register_MOM_IS_diag_field, safe_alloc_ptr +use MOM_IS_diag_mediator, only : set_axes_info +use MOM_IS_diag_mediator, only : diag_mediator_init, set_diag_mediator_grid, diag_ctrl, time_type +use MOM_IS_diag_mediator, only : enable_averages, enable_averaging, disable_averaging +use MOM_IS_diag_mediator, only : diag_mediator_infrastructure_init, diag_mediator_close_registration use MOM_domains, only : MOM_domains_init, clone_MOM_domain use MOM_domains, only : pass_var, pass_vector, TO_ALL, CGRID_NE, BGRID_NE, CORNER use MOM_dyn_horgrid, only : dyn_horgrid_type, create_dyn_horgrid, destroy_dyn_horgrid @@ -19,21 +22,25 @@ module MOM_ice_shelf 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 use MOM_grid_initialize, only : set_grid_metrics +use MOM_hor_index, only : hor_index_type, hor_index_init +use MOM_hor_index, only : rotate_hor_index use MOM_fixed_initialization, only : MOM_initialize_topography use MOM_fixed_initialization, only : MOM_initialize_rotation use user_initialization, only : user_initialize_topography use MOM_io, only : field_exists, file_exists, MOM_read_data, write_version_number -use MOM_io, only : slasher, fieldtype +use MOM_io, only : slasher, fieldtype, vardesc, var_desc use MOM_io, only : write_field, close_file, SINGLE_FILE, MULTIPLE use MOM_restart, only : register_restart_field, query_initialized, save_restart -use MOM_restart, only : restart_init, restore_state, MOM_restart_CS +use MOM_restart, only : restart_init, restore_state, MOM_restart_CS, register_restart_pair use MOM_time_manager, only : time_type, time_type_to_real, real_to_time, operator(>), operator(-) use MOM_transcribe_grid, only : copy_dyngrid_to_MOM_grid, copy_MOM_grid_to_dyngrid +use MOM_transcribe_grid, only : rotate_dyngrid use MOM_unit_scaling, only : unit_scale_type, unit_scaling_init, fix_restart_unit_scaling -use MOM_variables, only : surface +use MOM_variables, only : surface, allocate_surface_state +use MOM_variables, only : rotate_surface_state use MOM_forcing_type, only : forcing, allocate_forcing_type, MOM_forcing_chksum use MOM_forcing_type, only : mech_forcing, allocate_mech_forcing, MOM_mech_forcing_chksum -use MOM_forcing_type, only : copy_common_forcing_fields +use MOM_forcing_type, only : copy_common_forcing_fields, rotate_forcing, rotate_mech_forcing use MOM_get_input, only : directories, Get_MOM_input use MOM_EOS, only : calculate_density, calculate_density_derivs, calculate_TFreeze, EOS_domain use MOM_EOS, only : EOS_type, EOS_init @@ -54,13 +61,13 @@ module MOM_ice_shelf implicit none ; private #include -#ifdef SYMMETRIC_LAND_ICE +#ifdef SYMMETRIC_MEMORY_ # define GRID_SYM_ .true. #else # define GRID_SYM_ .false. #endif -public shelf_calc_flux, add_shelf_flux, initialize_ice_shelf, ice_shelf_end +public shelf_calc_flux, initialize_ice_shelf, ice_shelf_end, ice_shelf_query public ice_shelf_save_restart, solo_step_ice_shelf, add_shelf_forces ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional @@ -72,11 +79,17 @@ module MOM_ice_shelf type, public :: ice_shelf_CS ; private ! Parameters type(MOM_restart_CS), pointer :: restart_CSp => NULL() !< A pointer to the restart control - !! structure for the ice shelves - type(ocean_grid_type) :: grid !< Grid for the ice-shelf model + !! structure for the ice shelves + type(ocean_grid_type), pointer :: Grid_in => NULL() !< un-rotated input grid metric + type(hor_index_type), pointer :: HI_in => NULL() !< Pointer to a horizontal indexing structure for + !! incoming data which has not been rotated. + type(hor_index_type), pointer :: HI => NULL() !< Pointer to a horizontal indexing structure for + !! incoming data which has not been rotated. + logical :: rotate_index = .false. !< True if index map is rotated + integer :: turns !< The number of quarter turns for rotation testing. + type(ocean_grid_type), pointer :: Grid => NULL() !< Grid for the ice-shelf model type(unit_scale_type), pointer :: & US => NULL() !< A structure containing various unit conversion factors - !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 !< A factor that can be used to turn off ice shelf @@ -181,28 +194,28 @@ module MOM_ice_shelf logical :: debug !< If true, write verbose checksums for debugging purposes !! and use reproducible sums + integer :: id_clock_shelf=-1 !< CPU Clock for the ice shelf code + integer :: id_clock_pass=-1 !< CPU Clock for group pass calls end type ice_shelf_CS -integer :: id_clock_shelf !< CPU Clock for the ice shelf code -integer :: id_clock_pass !< CPU Clock for group pass calls + contains !> Calculates fluxes between the ocean and ice-shelf using the three-equations !! formulation (optional to use just two equations). !! See \ref section_ICE_SHELF_equations -subroutine shelf_calc_flux(sfc_state, fluxes, Time, time_step, CS, forces) - type(surface), intent(inout) :: sfc_state !< A structure containing fields that +subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step, CS) + type(surface), target, intent(inout) :: sfc_state_in !< A structure containing fields that !! describe the surface state of the ocean. The !! intent is only inout to allow for halo updates. - type(forcing), intent(inout) :: fluxes !< structure containing pointers to any possible + type(forcing), pointer :: fluxes_in !< structure containing pointers to any possible !! thermodynamic or mass-flux forcing fields. type(time_type), intent(in) :: Time !< Start time of the fluxes. real, intent(in) :: time_step !< Length of time over which these fluxes !! will be applied [s]. type(ice_shelf_CS), pointer :: CS !< A pointer to the control structure returned !! by a previous call to initialize_ice_shelf. - type(mech_forcing), optional, intent(inout) :: forces !< A structure with the driving mechanical forces ! Local variables type(ocean_grid_type), pointer :: G => NULL() !< The grid structure used by the ice shelf. @@ -211,6 +224,9 @@ subroutine shelf_calc_flux(sfc_state, fluxes, Time, time_step, CS, forces) type(ice_shelf_state), pointer :: ISS => NULL() !< A structure with elements that describe !! the ice-shelf state + type(surface), pointer :: sfc_state => NULL() + type(forcing), pointer :: fluxes => NULL() + real, dimension(SZI_(CS%grid)) :: & Rhoml, & !< Ocean mixed layer density [R ~> kg m-3]. dR0_dT, & !< Partial derivative of the mixed layer density @@ -282,6 +298,7 @@ subroutine shelf_calc_flux(sfc_state, fluxes, Time, time_step, CS, forces) logical :: update_ice_vel ! If true, it is time to update the ice shelf velocities. logical :: coupled_GL ! If true, the grouding line position is determined based on ! coupled ice-ocean dynamics. + logical :: use_temperature = .true. ! real, parameter :: c2_3 = 2.0/3.0 character(len=160) :: mesg ! The text of an error message @@ -290,11 +307,21 @@ subroutine shelf_calc_flux(sfc_state, fluxes, Time, time_step, CS, forces) if (.not. associated(CS)) call MOM_error(FATAL, "shelf_calc_flux: "// & "initialize_ice_shelf must be called before shelf_calc_flux.") - call cpu_clock_begin(id_clock_shelf) + call cpu_clock_begin(CS%id_clock_shelf) G => CS%grid ; US => CS%US ISS => CS%ISS + if (CS%rotate_index) then + allocate(sfc_state) + call rotate_surface_state(sfc_state_in,CS%Grid_in, sfc_state,CS%Grid,CS%turns) + allocate(fluxes) + call allocate_forcing_type(fluxes_in,G,fluxes) + call rotate_forcing(fluxes_in,fluxes,CS%turns) + else + sfc_state=>sfc_state_in + fluxes=>fluxes_in + endif ! useful parameters is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; ied = G%ied ; jed = G%jed I_ZETA_N = 1.0 / ZETA_N @@ -329,12 +356,12 @@ subroutine shelf_calc_flux(sfc_state, fluxes, Time, time_step, CS, forces) endif if (CS%debug) then - call hchksum(fluxes%frac_shelf_h, "frac_shelf_h before apply melting", G%HI, haloshift=0) - call hchksum(sfc_state%sst, "sst before apply melting", G%HI, haloshift=0) - call hchksum(sfc_state%sss, "sss before apply melting", G%HI, haloshift=0) - call hchksum(sfc_state%u, "u_ml before apply melting", G%HI, haloshift=0, scale=US%L_T_to_m_s) - call hchksum(sfc_state%v, "v_ml before apply melting", G%HI, haloshift=0, scale=US%L_T_to_m_s) - call hchksum(sfc_state%ocean_mass, "ocean_mass before apply melting", G%HI, haloshift=0, & + call hchksum(fluxes_in%frac_shelf_h, "frac_shelf_h before apply melting", CS%Grid_in%HI, haloshift=0) + call hchksum(sfc_state_in%sst, "sst before apply melting", CS%Grid_in%HI, haloshift=0) + call hchksum(sfc_state_in%sss, "sss before apply melting", CS%Grid_in%HI, haloshift=0) + call uchksum(sfc_state_in%u, "u_ml before apply melting", CS%Grid_in%HI, haloshift=0, scale=US%L_T_to_m_s) + call vchksum(sfc_state_in%v, "v_ml before apply melting", CS%Grid_in%HI, haloshift=0, scale=US%L_T_to_m_s) + call hchksum(sfc_state_in%ocean_mass, "ocean_mass before apply melting", CS%Grid_in%HI, haloshift=0, & scale=US%RZ_to_kg_m2) endif @@ -564,13 +591,17 @@ subroutine shelf_calc_flux(sfc_state, fluxes, Time, time_step, CS, forces) if (dS_it < 0.0) then ! Sbdry is now the upper bound. - if (Sb_max_set .and. (Sbdry(i,j) > Sb_max)) & - call MOM_error(FATAL,"shelf_calc_flux: Irregular iteration for Sbdry (max).") + if (Sb_max_set) then + if (Sbdry(i,j) > Sb_max) & + call MOM_error(FATAL,"shelf_calc_flux: Irregular iteration for Sbdry (max).") + endif Sb_max = Sbdry(i,j) ; dS_max = dS_it ; Sb_max_set = .true. else ! Sbdry is now the lower bound. - if (Sb_min_set .and. (Sbdry(i,j) < Sb_min)) & - call MOM_error(FATAL, "shelf_calc_flux: Irregular iteration for Sbdry (min).") - Sb_min = Sbdry(i,j) ; dS_min = dS_it ; Sb_min_set = .true. + if (Sb_min_set) then + if (Sbdry(i,j) < Sb_min) & + call MOM_error(FATAL, "shelf_calc_flux: Irregular iteration for Sbdry (min).") + endif + Sb_min = Sbdry(i,j) ; dS_min = dS_it ; Sb_min_set = .true. endif ! dS_it < 0.0 if (Sb_min_set .and. Sb_max_set) then @@ -626,7 +657,7 @@ subroutine shelf_calc_flux(sfc_state, fluxes, Time, time_step, CS, forces) fluxes%iceshelf_melt(i,j) = 0.0 endif ! Compute haline driving, which is one of the diags. used in ISOMIP - haline_driving(i,j) = (ISS%water_flux(i,j) * Sbdry(i,j)) / (CS%Rho_ocn * exch_vel_s(i,j)) + if (exch_vel_s(i,j)>0.) haline_driving(i,j) = (ISS%water_flux(i,j) * Sbdry(i,j)) / (CS%Rho_ocn * exch_vel_s(i,j)) !!!!!!!!!!!!!!!!!!!!!!!!!!!!Safety checks !!!!!!!!!!!!!!!!!!!!!!!!! !1)Check if haline_driving computed above is consistent with @@ -659,10 +690,10 @@ subroutine shelf_calc_flux(sfc_state, fluxes, Time, time_step, CS, forces) enddo ; enddo ! i- and j-loops if (CS%active_shelf_dynamics .or. CS%override_shelf_movement) then - call cpu_clock_begin(id_clock_pass) + call cpu_clock_begin(CS%id_clock_pass) call pass_var(ISS%area_shelf_h, G%domain, complete=.false.) call pass_var(ISS%mass_shelf, G%domain) - call cpu_clock_end(id_clock_pass) + call cpu_clock_end(CS%id_clock_pass) endif ! Melting has been computed, now is time to update thickness and mass @@ -712,12 +743,14 @@ subroutine shelf_calc_flux(sfc_state, fluxes, Time, time_step, CS, forces) if (CS%id_h_mask > 0) call post_data(CS%id_h_mask,ISS%hmask,CS%diag) call disable_averaging(CS%diag) - if (present(forces)) then - call add_shelf_forces(G, US, CS, forces, do_shelf_area=(CS%active_shelf_dynamics .or. & - CS%override_shelf_movement)) + + call cpu_clock_end(CS%id_clock_shelf) + + if (CS%rotate_index) then +! call rotate_surface_state(sfc_state,CS%Grid, sfc_state_in,CS%Grid_in,-CS%turns) + call rotate_forcing(fluxes,fluxes_in,-CS%turns) endif - call cpu_clock_end(id_clock_shelf) if (CS%debug) call MOM_forcing_chksum("End of shelf calc flux", fluxes, G, CS%US, haloshift=0) @@ -772,39 +805,65 @@ end subroutine change_thickness_using_melt !> This subroutine adds the mechanical forcing fields and perhaps shelf areas, based on !! the ice state in ice_shelf_CS. -subroutine add_shelf_forces(G, US, CS, forces, do_shelf_area) - type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. +subroutine add_shelf_forces(Ocn_grid, US, CS, forces_in, do_shelf_area, external_call) + type(ocean_grid_type), intent(in) :: Ocn_grid !< The ocean's grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(ice_shelf_CS), pointer :: CS !< This module's control structure. - type(mech_forcing), intent(inout) :: forces !< A structure with the driving mechanical forces + type(mech_forcing), pointer :: forces_in !< A structure with the driving mechanical forces logical, optional, intent(in) :: do_shelf_area !< If true find the shelf-covered areas. - + logical, optional, intent(in) :: external_call !< If true the incoming forcing type + !! is using the input grid metric and needs + !! to be rotated. + type(ocean_grid_type), pointer :: G => NULL() !< A pointer to the ocean grid metric. + type(mech_forcing), pointer :: forces => NULL() !< A structure with the driving mechanical forces real :: kv_rho_ice ! The viscosity of ice divided by its density [L4 T-1 R-1 Z-2 ~> m5 kg-1 s-1]. real :: press_ice ! The pressure of the ice shelf per unit area of ocean (not ice) [R L2 T-2 ~> Pa]. logical :: find_area ! If true find the shelf areas at u & v points. + logical :: rotate = .false. type(ice_shelf_state), pointer :: ISS => NULL() ! A structure with elements that describe ! the ice-shelf state integer :: i, j, is, ie, js, je, isd, ied, jsd, jed + + if (present(external_call)) rotate=external_call + + if (CS%rotate_index .and. rotate) then + if ((Ocn_grid%isc /= CS%Grid_in%isc) .or. (Ocn_grid%iec /= CS%Grid_in%iec) .or. & + (Ocn_grid%jsc /= CS%Grid_in%jsc) .or. (Ocn_grid%jec /= CS%Grid_in%jec)) & + call MOM_error(FATAL,"add_shelf_forces: Incompatible Ocean and Ice shelf grids.") + + allocate(forces) + call allocate_mech_forcing(forces_in, CS%Grid, forces) + call rotate_mech_forcing(forces_in, CS%turns, forces) + else + if ((Ocn_grid%isc /= CS%Grid%isc) .or. (Ocn_grid%iec /= CS%Grid%iec) .or. & + (Ocn_grid%jsc /= CS%Grid%jsc) .or. (Ocn_grid%jec /= CS%Grid%jec)) & + call MOM_error(FATAL,"add_shelf_forces: Incompatible Ocean and Ice shelf grids.") + + forces=>forces_in + endif + + G=>CS%Grid + + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec isd = G%isd ; jsd = G%jsd ; ied = G%ied ; jed = G%jed - if ((CS%grid%isc /= G%isc) .or. (CS%grid%iec /= G%iec) .or. & - (CS%grid%jsc /= G%jsc) .or. (CS%grid%jec /= G%jec)) & - call MOM_error(FATAL,"add_shelf_forces: Incompatible ocean and ice shelf grids.") - ISS => CS%ISS find_area = .true. ; if (present(do_shelf_area)) find_area = do_shelf_area if (find_area) then ! The frac_shelf is set over the widest possible area. Could it be smaller? +! do j=js,je ; do I=is-1,ie do j=jsd,jed ; do I=isd,ied-1 forces%frac_shelf_u(I,j) = 0.0 if ((G%areaT(i,j) + G%areaT(i+1,j) > 0.0)) & ! .and. (G%areaCu(I,j) > 0.0)) & forces%frac_shelf_u(I,j) = (ISS%area_shelf_h(i,j) + ISS%area_shelf_h(i+1,j)) / & (G%areaT(i,j) + G%areaT(i+1,j)) enddo ; enddo +! do J=js-1,je ; do i=is,ie do J=jsd,jed-1 ; do i=isd,ied forces%frac_shelf_v(i,J) = 0.0 if ((G%areaT(i,j) + G%areaT(i,j+1) > 0.0)) & ! .and. (G%areaCv(i,J) > 0.0)) & @@ -843,29 +902,41 @@ subroutine add_shelf_forces(G, US, CS, forces, do_shelf_area) enddo ; enddo if (CS%debug) then - call uvchksum("rigidity_ice_[uv]", forces%rigidity_ice_u, forces%rigidity_ice_v, & - G%HI, symmetric=.true., scale=US%L_to_m**3*US%L_to_Z*US%s_to_T) - call uvchksum("frac_shelf_[uv]", forces%frac_shelf_u, forces%frac_shelf_v, & - G%HI, symmetric=.true.) + call uvchksum("rigidity_ice_[uv]", forces%rigidity_ice_u, & + forces%rigidity_ice_v, CS%Grid%HI, symmetric=.true., & + scale=US%L_to_m**3*US%L_to_Z*US%s_to_T, scalar_pair=.true.) + call uvchksum("frac_shelf_[uv]", forces%frac_shelf_u, & + forces%frac_shelf_v, CS%Grid%HI, symmetric=.true., & + scalar_pair=.true.) + endif + + if (CS%rotate_index .and. rotate) then + call rotate_mech_forcing(forces, -CS%turns, forces_in) + ! TODO: deallocate mech forcing? endif end subroutine add_shelf_forces !> This subroutine adds the ice shelf pressure to the fluxes type. -subroutine add_shelf_pressure(G, US, CS, fluxes) - type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - type(ice_shelf_CS), intent(in) :: CS !< This module's control structure. - type(forcing), intent(inout) :: fluxes !< A structure of surface fluxes that may be updated. +subroutine add_shelf_pressure(Ocn_grid, US, CS, fluxes) + type(ocean_grid_type), intent(in) :: Ocn_grid !< The ocean's grid structure. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(ice_shelf_CS), intent(in) :: CS !< This module's control structure. + type(forcing), pointer :: fluxes !< A structure of surface fluxes that may be updated. + type(ocean_grid_type), pointer :: G => NULL() ! A pointer to ocean's grid structure. real :: press_ice !< The pressure of the ice shelf per unit area of ocean (not ice) [R L2 T-2 ~> Pa]. integer :: i, j, is, ie, js, je, isd, ied, jsd, jed + + + G=>CS%Grid is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec if ((CS%grid%isc /= G%isc) .or. (CS%grid%iec /= G%iec) .or. & (CS%grid%jsc /= G%jsc) .or. (CS%grid%jec /= G%jec)) & call MOM_error(FATAL,"add_shelf_pressure: Incompatible ocean and ice shelf grids.") + do j=js,je ; do i=is,ie press_ice = (CS%ISS%area_shelf_h(i,j) * G%IareaT(i,j)) * (CS%g_Earth * CS%ISS%mass_shelf(i,j)) if (associated(fluxes%p_surf)) then @@ -886,7 +957,7 @@ subroutine add_shelf_flux(G, US, CS, sfc_state, fluxes) type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(ice_shelf_CS), pointer :: CS !< This module's control structure. type(surface), intent(inout) :: sfc_state !< Surface ocean state - type(forcing), intent(inout) :: fluxes !< A structure of surface fluxes that may be used/updated. + type(forcing), pointer :: fluxes !< A structure of surface fluxes that may be used/updated. ! local variables real :: frac_shelf !< The fractional area covered by the ice shelf [nondim]. @@ -1080,15 +1151,19 @@ end subroutine add_shelf_flux !> Initializes shelf model data, parameters and diagnostics -subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fluxes, Time_in, solo_ice_sheet_in) +subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, & + fluxes_in, sfc_state_in, Time_in, solo_ice_sheet_in) type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters type(ocean_grid_type), pointer :: ocn_grid !< The calling ocean model's horizontal grid structure type(time_type), intent(inout) :: Time !< The clock that that will indicate the model time type(ice_shelf_CS), pointer :: CS !< A pointer to the ice shelf control structure - type(diag_ctrl), target, intent(in) :: diag !< A structure that is used to regulate the diagnostic output. - type(forcing), optional, intent(inout) :: fluxes !< A structure containing pointers to any possible - !! thermodynamic or mass-flux forcing fields. - type(mech_forcing), optional, intent(inout) :: forces !< A structure with the driving mechanical forces + type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the diagnostic output. + type(mech_forcing), optional, pointer :: forces_in !< A structure with the driving mechanical forces + type(forcing), optional, pointer :: fluxes_in !< A structure containing pointers to any possible + !! thermodynamic or mass-flux forcing fields. + type(surface), target, optional, intent(inout) :: sfc_state_in !< A structure containing fields that + !! describe the surface state of the ocean. The + !! intent is only inout to allow for halo updates. type(time_type), optional, intent(in) :: Time_in !< The time at initialization. logical, optional, intent(in) :: solo_ice_sheet_in !< If present, this indicates whether !! a solo ice-sheet driver. @@ -1100,6 +1175,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl !! the ice-shelf state type(directories) :: dirs type(dyn_horgrid_type), pointer :: dG => NULL() + type(dyn_horgrid_type), pointer :: dG_in => NULL() real :: Z_rescale ! A rescaling factor for heights from the representation in ! a restart file to the internal representation in this run. real :: RZ_rescale ! A rescaling factor for mass loads from the representation in @@ -1119,10 +1195,18 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl integer :: i, j, is, ie, js, je, isd, ied, jsd, jed, Isdq, Iedq, Jsdq, Jedq integer :: wd_halos(2) logical :: read_TideAmp, shelf_mass_is_dynamic, debug + logical :: global_indexing character(len=240) :: Tideamp_file real :: utide ! A tidal velocity [L T-1 ~> m s-1] real :: col_thick_melt_thresh ! An ocean column thickness below which iceshelf melting - ! does not occur [Z ~> m] + ! does not occur [Z ~> m] + real, allocatable, dimension(:,:) :: tmp2d ! Temporary array for storing ice shelf input data + + type(mech_forcing), pointer :: forces => NULL() + type(forcing), pointer :: fluxes => NULL() + type(surface), pointer :: sfc_state => NULL() + type(vardesc) :: u_desc, v_desc + if (associated(CS)) then call MOM_error(FATAL, "MOM_ice_shelf.F90, initialize_ice_shelf: "// & "called with an associated control structure.") @@ -1135,37 +1219,86 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl ! MOM's grid and infrastructure. call Get_MOM_Input(dirs=dirs) + call diag_mediator_infrastructure_init() + ! Determining the internal unit scaling factors for this run. call unit_scaling_init(param_file, CS%US) + call get_param(param_file, mdl, "ROTATE_INDEX", CS%rotate_index, & + "Enable rotation of the horizontal indices.", default=.false., & + debuggingParam=.true.) + ! Set up the ice-shelf domain and grid wd_halos(:)=0 - 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, CS%US) + allocate(CS%Grid_in) + call MOM_domains_init(CS%Grid_in%domain, param_file, min_halo=wd_halos, symmetric=GRID_SYM_,& + domain_name='MOM_Ice_Shelf_in') + allocate(CS%HI_in) + call hor_index_init(CS%Grid_in%Domain, CS%HI_in, param_file, & + local_indexing=.not.global_indexing) + call MOM_grid_init(CS%Grid_in, param_file, CS%US, CS%HI_in) + + if (CS%rotate_index) then + ! TODO: Index rotation currently only works when index rotation does not + ! change the MPI rank of each domain. Resolving this will require a + ! modification to FMS PE assignment. + ! For now, we only permit single-core runs. + + if (num_PEs() /= 1) & + call MOM_error(FATAL, "Index rotation is only supported on one PE.") + + call get_param(param_file, mdl, "INDEX_TURNS", CS%turns, & + "Number of counterclockwise quarter-turn index rotations.", & + default=1, debuggingParam=.true.) + ! NOTE: If indices are rotated, then CS%Grid and CS%Grid_in must both be initialized. + ! If not rotated, then CS%Grid_in and CS%Ggrid are the same grid. + allocate(CS%Grid) + allocate(CS%HI) + call clone_MOM_domain(CS%Grid_in%Domain, CS%Grid%Domain,turns=CS%turns) + call rotate_hor_index(CS%HI_in, CS%turns, CS%HI) + call MOM_grid_init(CS%Grid, param_file, CS%US, CS%HI) + call create_dyn_horgrid(dG, CS%HI) + call create_dyn_horgrid(dG_in, CS%HI_in) + call clone_MOM_domain(CS%Grid_in%Domain, dG_in%Domain) + ! Set up the bottom depth, G%D either analytically or from file + call set_grid_metrics(dG_in,param_file,CS%US) + call MOM_initialize_topography(dG_in%bathyT, CS%Grid_in%max_depth, dG_in, param_file) + call rescale_dyn_horgrid_bathymetry(dG_in, CS%US%Z_to_m) + call rotate_dyngrid(dG_in, dG, CS%US, CS%turns) + call copy_dyngrid_to_MOM_grid(dG,CS%Grid,CS%US) + else + ! this needs to be fixed - will probably break when not using coupled driver 0 + CS%Grid=>CS%Grid_in + CS%HI=>CS%HI_in + call create_dyn_horgrid(dG, CS%Grid%HI) + call clone_MOM_domain(CS%Grid%Domain,dG%Domain) + call set_grid_metrics(dG,param_file,CS%US) + ! Set up the bottom depth, G%D either analytically or from file + call MOM_initialize_topography(dG%bathyT, CS%Grid%max_depth, dG, param_file) + call rescale_dyn_horgrid_bathymetry(dG, CS%US%Z_to_m) + call copy_dyngrid_to_MOM_grid(dG,CS%Grid,CS%US) + endif + G=>CS%Grid - call create_dyn_horgrid(dG, CS%grid%HI) - call clone_MOM_domain(CS%grid%Domain, dG%Domain) + allocate(diag) + call diag_mediator_init(G, param_file,diag,component='MOM_IceShelf') + ! This call sets up the diagnostic axes. These are needed, + ! e.g. to generate the target grids below. + call set_axes_info(G, param_file, diag) - call set_grid_metrics(dG, param_file, CS%US) - ! call set_diag_mediator_grid(CS%grid, CS%diag) + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + isd = G%isd ; jsd = G%jsd ; ied = G%ied ; jed = G%jed + Isdq = G%IsdB ; Iedq = G%IedB ; Jsdq = G%JsdB ; Jedq = G%JedB ! The ocean grid possibly uses different symmetry. if (associated(ocn_grid)) then ; CS%ocn_grid => ocn_grid else ; CS%ocn_grid => CS%grid ; endif ! Convenience pointers - G => CS%grid OG => CS%ocn_grid US => CS%US - - if (is_root_pe()) then - write(0,*) 'OG: ', OG%isd, OG%isc, OG%iec, OG%ied, OG%jsd, OG%jsc, OG%jsd, OG%jed - write(0,*) 'IG: ', G%isd, G%isc, G%iec, G%ied, G%jsd, G%jsc, G%jsd, G%jed - endif - - CS%diag => diag + CS%diag=>diag ! Are we being called from the solo ice-sheet driver? When called by the ocean ! model solo_ice_sheet_in is not preset. @@ -1174,9 +1307,6 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl if (present(Time_in)) Time = Time_in - is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec - isd = G%isd ; jsd = G%jsd ; ied = G%ied ; jed = G%jed - Isdq = G%IsdB ; Iedq = G%IedB ; Jsdq = G%JsdB ; Jedq = G%JedB CS%override_shelf_movement = .false. ; CS%active_shelf_dynamics = .false. @@ -1335,6 +1465,25 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl call get_param(param_file, mdl, "READ_TIDEAMP", read_TIDEAMP, & "If true, read a file (given by TIDEAMP_FILE) containing "//& "the tidal amplitude with INT_TIDE_DISSIPATION.", default=.false.) + call get_param(param_file, "MOM", "GLOBAL_INDEXING", global_indexing, & + "If true, use a global lateral indexing convention, so "//& + "that corresponding points on different processors have "//& + "the same index. This does not work with static memory.", & + default=.false., layoutParam=.true.) + + + if (PRESENT(sfc_state_in)) then + allocate(sfc_state) + ! assuming frazil is enabled in ocean. This could break some configurations? + call allocate_surface_state(sfc_state_in, CS%Grid_in, use_temperature=.true.,& + do_integrals=.true.,omit_frazil=.false.,use_iceshelves=.true.) + if (CS%rotate_index) then + call rotate_surface_state(sfc_state_in,CS%Grid_in, sfc_state,CS%Grid,CS%turns) + else + sfc_state=>sfc_state_in + endif + endif + call safe_alloc_ptr(CS%utide,isd,ied,jsd,jed) ; CS%utide(:,:) = 0.0 @@ -1346,7 +1495,17 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".") inputdir = slasher(inputdir) TideAmp_file = trim(inputdir) // trim(TideAmp_file) - call MOM_read_data(TideAmp_file, 'tideamp', CS%utide, G%domain, timelevel=1, scale=US%m_s_to_L_T) + if (CS%rotate_index) then + allocate(tmp2d(CS%HI_in%isd:CS%HI_in%ied,CS%HI_in%jsd:CS%HI_in%jed));tmp2d(:,:)=0.0 + allocate(tmp2d(CS%Grid_in%isd:CS%Grid_in%ied,CS%Grid_in%jsd:CS%Grid_in%jed));tmp2d(:,:)=0.0 + + + call MOM_read_data(TideAmp_file, 'tideamp', tmp2d, CS%Grid_in%domain, timelevel=1, scale=US%m_s_to_L_T) + call rotate_array(tmp2d,CS%turns, CS%utide) + deallocate(tmp2d) + else + call MOM_read_data(TideAmp_file, 'tideamp', CS%utide, CS%Grid%domain, timelevel=1, scale=US%m_s_to_L_T) + endif else call get_param(param_file, mdl, "UTIDE", utide, & "The constant tidal amplitude used with INT_TIDE_DISSIPATION.", & @@ -1404,30 +1563,47 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl ! GMM: the following assures that water/heat fluxes are just allocated ! when SHELF_THERMO = True. These fluxes are necessary if one wants to ! use either ENERGETICS_SFC_PBL (ALE mode) or BULKMIXEDLAYER (layer mode). - if (present(fluxes)) & - call allocate_forcing_type(CS%ocn_grid, fluxes, ustar=.true., shelf=.true., & - press=.true., water=CS%isthermo, heat=CS%isthermo) - if (present(forces)) & - call allocate_mech_forcing(CS%ocn_grid, forces, ustar=.true., shelf=.true., press=.true.) + if (present(fluxes_in)) then + call allocate_forcing_type(CS%Grid_in, fluxes_in, ustar=.true., shelf=.true., & + press=.true., water=CS%isthermo, heat=CS%isthermo) + if (CS%rotate_index) then + allocate(fluxes) + call allocate_forcing_type(fluxes_in, CS%Grid, fluxes) + call rotate_forcing(fluxes_in, fluxes, CS%turns) + else + fluxes=>fluxes_in + endif + endif + if (present(forces_in)) then + call allocate_mech_forcing(CS%Grid_in, forces_in, ustar=.true., shelf=.true., press=.true.) + if (CS%rotate_index) then + allocate(forces) + call allocate_mech_forcing(forces_in, CS%Grid, forces) + call rotate_mech_forcing(forces_in, CS%turns, forces) + else + forces=>forces_in + endif + endif else call MOM_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: allocating fluxes in solo mode.") - if (present(fluxes)) & - call allocate_forcing_type(G, fluxes, ustar=.true., shelf=.true., press=.true.) - if (present(forces)) & - call allocate_mech_forcing(G, forces, ustar=.true., shelf=.true., press=.true.) + if (present(fluxes_in)) then + call allocate_forcing_type(CS%Grid_in, fluxes_in, ustar=.true., shelf=.true., press=.true.) + if (CS%rotate_index) then + allocate(fluxes) + call allocate_forcing_type(fluxes_in, CS%Grid, fluxes) + call rotate_forcing(fluxes_in, fluxes, CS%turns) + endif + endif + if (present(forces_in)) then + call allocate_mech_forcing(CS%Grid_in, forces_in, ustar=.true., shelf=.true., press=.true.) + if (CS%rotate_index) then + allocate(forces) + call allocate_mech_forcing(forces_in, CS%Grid, forces) + call rotate_mech_forcing(forces_in, CS%turns, forces) + endif + endif endif - ! Set up the bottom depth, G%D either analytically or from file - call MOM_initialize_topography(dG%bathyT, G%max_depth, dG, param_file) - call rescale_dyn_horgrid_bathymetry(dG, US%Z_to_m) - - ! Set up the Coriolis parameter, G%f, usually analytically. - call MOM_initialize_rotation(dG%CoriolisBu, dG, param_file, US) - ! This copies grid elements, including bathyT and CoriolisBu from dG to CS%grid. - call copy_dyngrid_to_MOM_grid(dG, CS%grid, US) - - call destroy_dyn_horgrid(dG) - ! Set up the restarts. call restart_init(param_file, CS%restart_CSp, "Shelf.res") call register_restart_field(ISS%mass_shelf, "shelf_mass", .true., CS%restart_CSp, & @@ -1436,6 +1612,19 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl "Ice shelf area in cell", "m2") call register_restart_field(ISS%h_shelf, "h_shelf", .true., CS%restart_CSp, & "ice sheet/shelf thickness", "m") + if (PRESENT(sfc_state_in)) then + if (allocated(sfc_state%taux_shelf) .and. allocated(sfc_state%tauy_shelf)) then + u_desc = var_desc("taux_shelf", "Pa", "the zonal stress on the ocean under ice shelves", & + hor_grid='Cu',z_grid='1') + v_desc = var_desc("tauy_shelf", "Pa", "the meridional stress on the ocean under ice shelves", & + hor_grid='Cv',z_grid='1') + call register_restart_pair(sfc_state%taux_shelf, sfc_state%tauy_shelf, u_desc,v_desc, & + .false., CS%restart_CSp) + endif + endif + + call register_restart_field(ISS%h_shelf, "_shelf", .true., CS%restart_CSp, & + "ice sheet/shelf thickness", "m") call register_restart_field(US%m_to_Z_restart, "m_to_Z", .false., CS%restart_CSp, & "Height unit conversion factor", "Z meter-1") call register_restart_field(US%m_to_L_restart, "m_to_L", .false., CS%restart_CSp, & @@ -1449,7 +1638,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl if (CS%active_shelf_dynamics) then ! Allocate CS%dCS and specify additional restarts for ice shelf dynamics - call register_ice_shelf_dyn_restarts(G, param_file, CS%dCS, CS%restart_CSp) + call register_ice_shelf_dyn_restarts(CS%Grid_in, param_file, CS%dCS, CS%restart_CSp) endif !GMM - I think we do not need to save ustar_shelf and iceshelf_melt in the restart file @@ -1464,6 +1653,11 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl if ((dirs%input_filename(1:1) == 'n') .and. & (LEN_TRIM(dirs%input_filename) == 1)) new_sim = .true. + ISS%area_shelf_h(:,:)=0.0 + ISS%h_shelf(:,:)=0.0 + ISS%hmask(:,:)=0.0 + ISS%mass_shelf(:,:)=0.0 + if (CS%override_shelf_movement .and. CS%mass_from_file) then ! initialize the ids for reading shelf mass from a netCDF @@ -1471,7 +1665,8 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl if (new_sim) then ! new simulation, initialize ice thickness as in the static case - call initialize_ice_thickness(ISS%h_shelf, ISS%area_shelf_h, ISS%hmask, G, US, param_file) + call initialize_ice_thickness(ISS%h_shelf, ISS%area_shelf_h, ISS%hmask, CS%Grid, CS%Grid_in, US, param_file, & + CS%rotate_index, CS%turns) ! next make sure mass is consistent with thickness do j=G%jsd,G%jed ; do i=G%isd,G%ied @@ -1497,16 +1692,20 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl endif if (new_sim .and. (.not. (CS%override_shelf_movement .and. CS%mass_from_file))) then - - ! This model is initialized internally or from a file. - call initialize_ice_thickness(ISS%h_shelf, ISS%area_shelf_h, ISS%hmask, G, US, param_file) - + ! This model is initialized internally or from a file. + call initialize_ice_thickness(ISS%h_shelf, ISS%area_shelf_h, ISS%hmask, CS%Grid, CS%Grid_in, US, param_file,& + CS%rotate_index, CS%turns) ! next make sure mass is consistent with thickness do j=G%jsd,G%jed ; do i=G%isd,G%ied if ((ISS%hmask(i,j) == 1) .or. (ISS%hmask(i,j) == 2)) then ISS%mass_shelf(i,j) = ISS%h_shelf(i,j)*CS%density_ice endif enddo ; enddo + if (CS%debug) then + call hchksum(ISS%mass_shelf, "IS init: mass_shelf", G%HI, haloshift=0, scale=US%RZ_to_kg_m2) + call hchksum(ISS%area_shelf_h, "IS init: area_shelf", G%HI, haloshift=0, scale=US%L_to_m*US%L_to_m) + call hchksum(ISS%hmask, "IS init: hmask", G%HI, haloshift=0) + endif ! else ! Previous block for new_sim=.T., this block restores the state. elseif (.not.new_sim) then @@ -1539,15 +1738,23 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl endif ! .not. new_sim +! do j=G%jsc,G%jec ; do i=G%isc,G%iec +! ISS%area_shelf_h(i,j) = ISS%area_shelf_h(i,j)*G%mask2dT(i,j) +! enddo; enddo + CS%Time = Time - call cpu_clock_begin(id_clock_pass) + CS%id_clock_shelf = cpu_clock_id('Ice shelf', grain=CLOCK_COMPONENT) + CS%id_clock_pass = cpu_clock_id(' Ice shelf halo updates', grain=CLOCK_ROUTINE) + + call cpu_clock_begin(CS%id_clock_pass) call pass_var(ISS%area_shelf_h, G%domain) call pass_var(ISS%h_shelf, G%domain) call pass_var(ISS%mass_shelf, G%domain) call pass_var(ISS%hmask, G%domain) call pass_var(G%bathyT, G%domain) - call cpu_clock_end(id_clock_pass) + call cpu_clock_end(CS%id_clock_pass) + do j=jsd,jed ; do i=isd,ied if (ISS%area_shelf_h(i,j) > G%areaT(i,j)) then @@ -1555,18 +1762,19 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl ISS%area_shelf_h(i,j) = G%areaT(i,j) endif enddo ; enddo - if (present(fluxes)) then ; do j=jsd,jed ; do i=isd,ied - if (G%areaT(i,j) > 0.0) fluxes%frac_shelf_h(i,j) = ISS%area_shelf_h(i,j) / G%areaT(i,j) + if (present(fluxes_in)) then ; do j=jsd,jed ; do i=isd,ied + if (G%areaT(i,j)>0.) fluxes%frac_shelf_h(i,j) = ISS%area_shelf_h(i,j) / G%areaT(i,j) enddo ; enddo ; endif if (CS%debug) then call hchksum(fluxes%frac_shelf_h, "IS init: frac_shelf_h", G%HI, haloshift=0) + call hchksum(ISS%area_shelf_h, "IS init: area_shelf_h", G%HI, haloshift=0, scale=US%L_to_m*US%L_to_m) endif - if (present(forces)) & + if (present(forces_in)) & call add_shelf_forces(G, US, CS, forces, do_shelf_area=.not.CS%solo_ice_sheet) - if (present(fluxes)) call add_shelf_pressure(G, US, CS, fluxes) + if (present(fluxes_in)) call add_shelf_pressure(ocn_grid, US, CS, fluxes) if (CS%active_shelf_dynamics .and. .not.CS%isthermo) then ISS%water_flux(:,:) = 0.0 @@ -1586,18 +1794,18 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl if (save_IC .and. .not.((dirs%input_filename(1:1) == 'r') .and. & (LEN_TRIM(dirs%input_filename) == 1))) then - call save_restart(dirs%output_directory, CS%Time, G, & + call save_restart(dirs%output_directory, CS%Time, CS%Grid_in, & CS%restart_CSp, filename=IC_file) endif - CS%id_area_shelf_h = register_diag_field('ocean_model', 'area_shelf_h', CS%diag%axesT1, CS%Time, & + CS%id_area_shelf_h = register_diag_field('ice_shelf_model', 'area_shelf_h', CS%diag%axesT1, CS%Time, & 'Ice Shelf Area in cell', 'meter-2', conversion=US%L_to_m**2) - CS%id_shelf_mass = register_diag_field('ocean_model', 'shelf_mass', CS%diag%axesT1, CS%Time, & + CS%id_shelf_mass = register_diag_field('ice_shelf_model', 'shelf_mass', CS%diag%axesT1, CS%Time, & 'mass of shelf', 'kg/m^2', conversion=US%RZ_to_kg_m2) - CS%id_h_shelf = register_diag_field('ocean_model', 'h_shelf', CS%diag%axesT1, CS%Time, & + CS%id_h_shelf = register_diag_field('ice_shelf_model', 'h_shelf', CS%diag%axesT1, CS%Time, & 'ice shelf thickness', 'm', conversion=US%Z_to_m) - CS%id_mass_flux = register_diag_field('ocean_model', 'mass_flux', CS%diag%axesT1,& + CS%id_mass_flux = register_diag_field('ice_shelf_model', 'mass_flux', CS%diag%axesT1,& CS%Time, 'Total mass flux of freshwater across the ice-ocean interface.', & 'kg/s', conversion=US%RZ_T_to_kg_m2s*US%L_to_m**2) @@ -1606,35 +1814,39 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl else ! use original eq. meltrate_conversion = 86400.0*365.0*US%Z_to_m*US%s_to_T / CS%density_ice endif - CS%id_melt = register_diag_field('ocean_model', 'melt', CS%diag%axesT1, CS%Time, & + CS%id_melt = register_diag_field('ice_shelf_model', 'melt', CS%diag%axesT1, CS%Time, & 'Ice Shelf Melt Rate', 'm yr-1', conversion= meltrate_conversion) - CS%id_thermal_driving = register_diag_field('ocean_model', 'thermal_driving', CS%diag%axesT1, CS%Time, & + CS%id_thermal_driving = register_diag_field('ice_shelf_model', 'thermal_driving', CS%diag%axesT1, CS%Time, & 'pot. temp. in the boundary layer minus freezing pot. temp. at the ice-ocean interface.', 'Celsius') - CS%id_haline_driving = register_diag_field('ocean_model', 'haline_driving', CS%diag%axesT1, CS%Time, & + CS%id_haline_driving = register_diag_field('ice_shelf_model', 'haline_driving', CS%diag%axesT1, CS%Time, & 'salinity in the boundary layer minus salinity at the ice-ocean interface.', 'psu') - CS%id_Sbdry = register_diag_field('ocean_model', 'sbdry', CS%diag%axesT1, CS%Time, & + CS%id_Sbdry = register_diag_field('ice_shelf_model', 'sbdry', CS%diag%axesT1, CS%Time, & 'salinity at the ice-ocean interface.', 'psu') - CS%id_u_ml = register_diag_field('ocean_model', 'u_ml', CS%diag%axesCu1, CS%Time, & + CS%id_u_ml = register_diag_field('ice_shelf_model', 'u_ml', CS%diag%axesCu1, CS%Time, & 'Eastward vel. in the boundary layer (used to compute ustar)', 'm s-1', conversion=US%L_T_to_m_s) - CS%id_v_ml = register_diag_field('ocean_model', 'v_ml', CS%diag%axesCv1, CS%Time, & + CS%id_v_ml = register_diag_field('ice_shelf_model', 'v_ml', CS%diag%axesCv1, CS%Time, & 'Northward vel. in the boundary layer (used to compute ustar)', 'm s-1', conversion=US%L_T_to_m_s) - CS%id_exch_vel_s = register_diag_field('ocean_model', 'exch_vel_s', CS%diag%axesT1, CS%Time, & + CS%id_exch_vel_s = register_diag_field('ice_shelf_model', 'exch_vel_s', CS%diag%axesT1, CS%Time, & 'Sub-shelf salinity exchange velocity', 'm s-1', conversion=US%Z_to_m*US%s_to_T) - CS%id_exch_vel_t = register_diag_field('ocean_model', 'exch_vel_t', CS%diag%axesT1, CS%Time, & + CS%id_exch_vel_t = register_diag_field('ice_shelf_model', 'exch_vel_t', CS%diag%axesT1, CS%Time, & 'Sub-shelf thermal exchange velocity', 'm s-1' , conversion=US%Z_to_m*US%s_to_T) - CS%id_tfreeze = register_diag_field('ocean_model', 'tfreeze', CS%diag%axesT1, CS%Time, & + CS%id_tfreeze = register_diag_field('ice_shelf_model', 'tfreeze', CS%diag%axesT1, CS%Time, & 'In Situ Freezing point at ice shelf interface', 'degC') - CS%id_tfl_shelf = register_diag_field('ocean_model', 'tflux_shelf', CS%diag%axesT1, CS%Time, & + CS%id_tfl_shelf = register_diag_field('ice_shelf_model', 'tflux_shelf', CS%diag%axesT1, CS%Time, & 'Heat conduction into ice shelf', 'W m-2', conversion=-US%QRZ_T_to_W_m2) - CS%id_ustar_shelf = register_diag_field('ocean_model', 'ustar_shelf', CS%diag%axesT1, CS%Time, & + CS%id_ustar_shelf = register_diag_field('ice_shelf_model', 'ustar_shelf', CS%diag%axesT1, CS%Time, & 'Fric vel under shelf', 'm/s', conversion=US%Z_to_m*US%s_to_T) if (CS%active_shelf_dynamics) then - CS%id_h_mask = register_diag_field('ocean_model', 'h_mask', CS%diag%axesT1, CS%Time, & + CS%id_h_mask = register_diag_field('ice_shelf_model', 'h_mask', CS%diag%axesT1, CS%Time, & 'ice shelf thickness mask', 'none') endif + call diag_mediator_close_registration(CS%diag) + - id_clock_shelf = cpu_clock_id('Ice shelf', grain=CLOCK_COMPONENT) - id_clock_pass = cpu_clock_id(' Ice shelf halo updates', grain=CLOCK_ROUTINE) + if (present(fluxes_in) .and. CS%rotate_index) & + call rotate_forcing(fluxes, fluxes_in, -CS%turns) + if (present(forces_in) .and. CS%rotate_index) & + call rotate_mech_forcing(forces, -CS%turns, forces_in) end subroutine initialize_ice_shelf @@ -1689,7 +1901,7 @@ subroutine initialize_shelf_mass(G, param_file, CS, ISS, new_sim) call log_param(param_file, mdl, "INPUTDIR/SHELF_FILE", filename) CS%id_read_mass = init_external_field(filename, shelf_mass_var, & - domain=G%Domain%mpp_domain, verbose=CS%debug) + domain=CS%Grid_in%Domain%mpp_domain, verbose=CS%debug) if (read_shelf_area) then call get_param(param_file, mdl, "SHELF_AREA_VAR", shelf_area_var, & @@ -1697,10 +1909,10 @@ subroutine initialize_shelf_mass(G, param_file, CS, ISS, new_sim) default="shelf_area") CS%id_read_area = init_external_field(filename,shelf_area_var, & - domain=G%Domain%mpp_domain) + domain=CS%Grid_in%Domain%mpp_domain) endif - if (.not.file_exists(filename, G%Domain)) call MOM_error(FATAL, & + if (.not.file_exists(filename, CS%Grid_in%Domain)) call MOM_error(FATAL, & " initialize_shelf_mass: Unable to open "//trim(filename)) case ("zero") @@ -1729,9 +1941,23 @@ subroutine update_shelf_mass(G, US, CS, ISS, Time) ! local variables integer :: i, j, is, ie, js, je + real, allocatable, dimension(:,:) :: tmp2d ! Temporary array for storing ice shelf input data + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec - call time_interp_external(CS%id_read_mass, Time, ISS%mass_shelf) + + if (CS%rotate_index) then + allocate(tmp2d(CS%Grid_in%isc:CS%Grid_in%iec,CS%Grid_in%jsc:CS%Grid_in%jec));tmp2d(:,:)=0.0 + else + allocate(tmp2d(is:ie,js:je)) + endif + + + + call time_interp_external(CS%id_read_mass, Time, tmp2d) + call rotate_array(tmp2d,CS%turns, ISS%mass_shelf) + deallocate(tmp2d) + ! This should only be done if time_interp_external did an update. do j=js,je ; do i=is,ie ISS%mass_shelf(i,j) = US%kg_m3_to_R*US%m_to_Z * ISS%mass_shelf(i,j) ! Rescale after time_interp @@ -1762,6 +1988,28 @@ subroutine update_shelf_mass(G, US, CS, ISS, Time) end subroutine update_shelf_mass +!> Save the ice shelf restart file +subroutine ice_shelf_query(CS, G, frac_shelf_h) + type(ice_shelf_CS), pointer :: CS !< ice shelf control structure + type(ocean_grid_type), intent(in) :: G !< A pointer to an ocean grid control structure. + real, optional, dimension(SZI_(G),SZJ_(G)) :: frac_shelf_h !< + !< Ice shelf area fraction [nodim]. + + logical :: do_frac=.false. + integer :: i,j + + if (present(frac_shelf_h)) do_frac=.true. + + if (do_frac) then + do j=G%jsd,G%jed + do i=G%isd,G%ied + frac_shelf_h(i,j)=0.0 + if (G%areaT(i,j)>0.) frac_shelf_h(i,j) = CS%ISS%area_shelf_h(i,j) / G%areaT(i,j) + enddo + enddo + endif + +end subroutine ice_shelf_query !> Save the ice shelf restart file subroutine ice_shelf_save_restart(CS, Time, directory, time_stamped, filename_suffix) type(ice_shelf_CS), pointer :: CS !< ice shelf control structure @@ -1781,7 +2029,7 @@ subroutine ice_shelf_save_restart(CS, Time, directory, time_stamped, filename_su if (present(directory)) then ; restart_dir = directory else ; restart_dir = CS%restart_output_dir ; endif - call save_restart(restart_dir, Time, CS%grid, CS%restart_CSp, time_stamped) + call save_restart(restart_dir, Time, CS%grid_in, CS%restart_CSp, time_stamped) end subroutine ice_shelf_save_restart diff --git a/src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 b/src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 new file mode 100644 index 0000000000..ab4245bd83 --- /dev/null +++ b/src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 @@ -0,0 +1,758 @@ +!> Convenient wrappers to the FMS diag_manager interfaces with additional diagnostic capabilies. +module MOM_IS_diag_mediator + +! This file is a part of SIS2. See LICENSE.md for the license. + +use MOM_grid, only : ocean_grid_type + +use MOM_coms, only : PE_here +use MOM_error_handler, only : MOM_error, FATAL, is_root_pe, assert +use MOM_file_parser, only : get_param, log_param, log_version, param_file_type +use MOM_safe_alloc, only : safe_alloc_ptr, safe_alloc_alloc +use MOM_string_functions, only : lowercase, uppercase, slasher +use MOM_time_manager, only : time_type + +use diag_manager_mod, only : diag_manager_init +use diag_manager_mod, only : send_data, diag_axis_init,EAST,NORTH +use diag_manager_mod, only : register_diag_field_fms=>register_diag_field +use diag_manager_mod, only : register_static_field_fms=>register_static_field + +implicit none ; private + +public diag_mediator_infrastructure_init +public set_axes_info, post_data, register_MOM_IS_diag_field, time_type +public safe_alloc_ptr, safe_alloc_alloc +public enable_averaging, disable_averaging, query_averaging_enabled +public enable_averages +public diag_mediator_init, diag_mediator_end, set_diag_mediator_grid +public diag_mediator_close_registration, get_diag_time_end +public diag_axis_init, register_static_field + +!> Make a diagnostic available for averaging or output. +!interface post_data +! module procedure post_data_2d +!end interface post_data + +!> 2D/3D axes type to contain 1D axes handles and pointers to masks +type, public :: axesType + character(len=15) :: id !< The id string for this particular combination of handles. + integer :: rank !< Number of dimensions in the list of axes. + integer, dimension(:), allocatable :: handles !< Handles to 1D axes. + type(diag_ctrl), pointer :: diag_cs => null() !< A structure that is used to regulate diagnostic output +end type axesType + +!> This type is used to represent a diagnostic at the diag_mediator level. +type, private :: diag_type + logical :: in_use !< This diagnostic is in use + integer :: fms_diag_id !< underlying FMS diag id + character(len=24) :: name !< The diagnostic name + real :: conversion_factor = 0. !< A factor to multiply data by before posting to FMS, if non-zero. + real, pointer, dimension(:,:) :: mask2d => null() !< A 2-d mask on the data domain for this diagnostic + real, pointer, dimension(:,:) :: mask2d_comp => null() !< A 2-d mask on the computational domain for this diagnostic +end type diag_type + +!> The SIS_diag_ctrl data type contains times to regulate diagnostics along with masks and +!! axes to use with diagnostics, and a list of structures with data about each diagnostic. +type, public :: diag_ctrl + integer :: doc_unit = -1 !< The unit number of a diagnostic documentation file. + !! This file is open if doc_unit is > 0. + + ! The following fields are used for the output of the data. + ! These give the computational-domain sizes, and are relative to a start value + ! of 1 in memory for the tracer-point arrays. + integer :: is !< The start i-index of cell centers within the computational domain + integer :: ie !< The end i-index of cell centers within the computational domain + integer :: js !< The start j-index of cell centers within the computational domain + integer :: je !< The end j-index of cell centers within the computational domain + ! These give the memory-domain sizes, and can be start at any value on each PE. + integer :: isd !< The start i-index of cell centers within the data domain + integer :: ied !< The end i-index of cell centers within the data domain + integer :: jsd !< The start j-index of cell centers within the data domain + integer :: jed !< The end j-index of cell centers within the data domain + real :: time_int !< The time interval in s for any fields that are offered for averaging. + type(time_type) :: time_end !< The end time of the valid interval for any offered field. + logical :: ave_enabled = .false. !< .true. if averaging is enabled. + + !>@{ The following are 3D and 2D axis groups defined for output. The names indicate + !! the horizontal locations (B, T, Cu, or Cv), vertical locations (L, i, or 1) and + !! thickness categories (c, c0, or 1). + type(axesType) :: axesBL, axesTL, axesCuL, axesCvL + type(axesType) :: axesBi, axesTi, axesCui, axesCvi + type(axesType) :: axesBc, axesTc, axesCuc, axesCvc + type(axesType) :: axesBc0, axesTc0, axesCuc0, axesCvc0 + type(axesType) :: axesB1, axesT1, axesCu1, axesCv1 + !!@} + + ! Mask arrays for diagnostics + real, dimension(:,:), pointer :: mask2dT => null() !< 2D mask array for cell-center points + real, dimension(:,:), pointer :: mask2dBu => null() !< 2D mask array for cell-corners + real, dimension(:,:), pointer :: mask2dCu => null() !< 2D mask array for east-faces + real, dimension(:,:), pointer :: mask2dCv => null() !< 2D mask array for north-faces + !> Computational domain mask arrays for diagnostics. + real, dimension(:,:), pointer :: mask2dT_comp => null() + +#define DIAG_ALLOC_CHUNK_SIZE 15 + type(diag_type), dimension(:), allocatable :: diags !< The array of diagnostics + integer :: next_free_diag_id !< The next unused diagnostic ID + !> default missing value to be sent to ALL diagnostics registerations + real :: missing_value = -1.0e34 + +end type diag_ctrl + +contains + +!> Set up the grid and axis information for use by the ice shelf model. +subroutine set_axes_info(G, param_file, diag_cs, axes_set_name) + type(ocean_grid_type), intent(inout) :: G !< The horizontal grid type + type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + type(diag_ctrl), intent(inout) :: diag_cs !< A structure that is used to regulate diagnostic output + character(len=*), optional, intent(in) :: axes_set_name !< A name to use for this set of axes. + !! The default is "ice". +! This subroutine sets up the grid and axis information for use by the ice shelf model. + + ! Local variables + integer :: id_xq, id_yq, id_zl, id_zi, id_xh, id_yh, id_ct, id_ct0 + integer :: k + logical :: Cartesian_grid + character(len=80) :: grid_config, units_temp, set_name +! This include declares and sets the variable "version". +#include "version_variable.h" + character(len=40) :: mdl = "MOM_IS_diag_mediator" ! This module's name. + + set_name = "ice_shelf" ; if (present(axes_set_name)) set_name = trim(axes_set_name) + + ! Read all relevant parameters and write them to the model log. + call log_version(param_file, mdl, version) + call get_param(param_file, mdl, "GRID_CONFIG", grid_config, & + "The method for defining the horizontal grid. Valid "//& + "entries include:\n"//& + "\t file - read the grid from GRID_FILE \n"//& + "\t mosaic - read the grid from a mosaic grid file \n"//& + "\t cartesian - a Cartesian grid \n"//& + "\t spherical - a spherical grid \n"//& + "\t mercator - a Mercator grid", fail_if_missing=.true.) + + G%x_axis_units = "degrees_E" + G%y_axis_units = "degrees_N" + if (index(lowercase(trim(grid_config)),"cartesian") > 0) then + ! This is a cartesian grid, and may have different axis units. + Cartesian_grid = .true. + call get_param(param_file, mdl, "AXIS_UNITS", units_temp, & + "The units for the x- and y- axis labels. AXIS_UNITS "//& + "should be defined as 'k' for km, 'm' for m, or 'd' "//& + "for degrees of latitude and longitude (the default). "//& + "Except on a Cartesian grid, only degrees are currently "//& + "implemented.", default='degrees') + if (units_temp(1:1) == 'k') then + G%x_axis_units = "kilometers" ; G%y_axis_units = "kilometers" + elseif (units_temp(1:1) == 'm') then + G%x_axis_units = "meters" ; G%y_axis_units = "meters" + endif + call log_param(param_file, mdl, "explicit AXIS_UNITS", G%x_axis_units) + else + Cartesian_grid = .false. + endif + + if (G%symmetric) then + id_xq = diag_axis_init('xB', G%gridLonB(G%isgB:G%iegB), G%x_axis_units, 'x', & + 'Boundary point nominal longitude',set_name=set_name, & + Domain2=G%Domain%mpp_domain, domain_position=EAST) + id_yq = diag_axis_init('yB', G%gridLatB(G%jsgB:G%jegB), G%y_axis_units, 'y', & + 'Boundary point nominal latitude', set_name=set_name, & + Domain2=G%Domain%mpp_domain, domain_position=NORTH) + + else + id_xq = diag_axis_init('xB', G%gridLonB(G%isg:G%ieg), G%x_axis_units, 'x', & + 'Boundary point nominal longitude',set_name=set_name, & + Domain2=G%Domain%mpp_domain, domain_position=EAST) + id_yq = diag_axis_init('yB', G%gridLatB(G%jsg:G%jeg), G%y_axis_units, 'y', & + 'Boundary point nominal latitude', set_name=set_name, & + Domain2=G%Domain%mpp_domain, domain_position=NORTH) + + endif + id_xh = diag_axis_init('xT', G%gridLonT(G%isg:G%ieg), G%x_axis_units, 'x', & + 'T point nominal longitude', set_name=set_name, & + Domain2=G%Domain%mpp_domain) + id_yh = diag_axis_init('yT', G%gridLatT(G%jsg:G%jeg), G%y_axis_units, 'y', & + 'T point nominal latitude', set_name=set_name, & + Domain2=G%Domain%mpp_domain) + + ! Axis groupings for 2-D arrays. + call defineAxes(diag_cs, [id_xh, id_yh], diag_cs%axesT1) + call defineAxes(diag_cs, [id_xq, id_yq], diag_cs%axesB1) + call defineAxes(diag_cs, [id_xq, id_yh], diag_cs%axesCu1) + call defineAxes(diag_cs, [id_xh, id_yq], diag_cs%axesCv1) + +end subroutine set_axes_info + +!> Define an a group of axes from a list of handles +subroutine defineAxes(diag_cs, handles, axes) + ! Defines "axes" from list of handle and associates mask + type(diag_ctrl), target, intent(in) :: diag_cs !< A structure that is used to regulate diagnostic output + integer, dimension(:), intent(in) :: handles !< A set of axis handles that define the axis group + type(axesType), intent(out) :: axes !< A group of axes that is set up here + + ! Local variables + integer :: n + n = size(handles) + if (n<1 .or. n>3) call MOM_error(FATAL,"defineAxes: wrong size for list of handles!") + allocate( axes%handles(n) ) + axes%id = i2s(handles, n) ! Identifying string + axes%rank = n + axes%handles(:) = handles(:) + axes%diag_cs => diag_cs ! A (circular) link back to the MOM_IS_diag_ctrl structure +end subroutine defineAxes + +!> Set up the current grid for the diag mediator +subroutine set_diag_mediator_grid(G, diag_cs) + type(ocean_grid_type), intent(inout) :: G !< The horizontal grid type + type(diag_ctrl), intent(inout) :: diag_cs !< A structure that is used to regulate diagnostic output + + diag_cs%is = G%isc - (G%isd-1) ; diag_cs%ie = G%iec - (G%isd-1) + diag_cs%js = G%jsc - (G%jsd-1) ; diag_cs%je = G%jec - (G%jsd-1) + diag_cs%isd = G%isd ; diag_cs%ied = G%ied ; diag_cs%jsd = G%jsd ; diag_cs%jed = G%jed +end subroutine set_diag_mediator_grid + +!> Offer a 2d diagnostic field for output or averaging +subroutine post_data(diag_field_id, field, diag_cs, is_static, mask) + integer, intent(in) :: diag_field_id !< the id for an output variable returned by a + !! previous call to register_diag_field. + real, target, intent(in) :: field(:,:) !< The 2-d array being offered for output or averaging. + type(diag_ctrl), target, & + intent(in) :: diag_cs !< A structure that is used to regulate diagnostic output + logical, optional, intent(in) :: is_static !< If true, this is a static field that is always offered. + logical, optional, intent(in) :: mask(:,:) !< If present, use this logical array as the data mask. + + ! Local variables + real, dimension(:,:), pointer :: locfield + logical :: used, is_stat + logical :: i_data, j_data + integer :: isv, iev, jsv, jev, i, j + integer :: fms_diag_id + type(diag_type), pointer :: diag => NULL() + + locfield => NULL() + is_stat = .false. ; if (present(is_static)) is_stat = is_static + + ! Get a pointer to the diag type for this id, and the FMS-level diag id. + call assert(diag_field_id < diag_cs%next_free_diag_id, & + 'post_data: Unregistered diagnostic id') + diag => diag_cs%diags(diag_field_id) + fms_diag_id = diag%fms_diag_id + + ! Determine the proper array indices, noting that because of the (:,:) + ! declaration of field, symmetric arrays are using a SW-grid indexing, + ! but non-symmetric arrays are using a NE-grid indexing. Send_data + ! actually only uses the difference between ie and is to determine + ! the output data size and assumes that halos are symmetric. + isv = diag_cs%is ; iev = diag_cs%ie ; jsv = diag_cs%js ; jev = diag_cs%je + + if ( size(field,1) == diag_cs%ied-diag_cs%isd +1 ) then + isv = diag_cs%is ; iev = diag_cs%ie ; i_data = .true. ! Data domain + elseif ( size(field,1) == diag_cs%ied-diag_cs%isd +2 ) then + isv = diag_cs%is ; iev = diag_cs%ie+1 ; i_data = .true. ! Symmetric data domain + elseif ( size(field,1) == diag_cs%ie-diag_cs%is +1 ) then + isv = 1 ; iev = diag_cs%ie + 1-diag_cs%is ; i_data = .false. ! Computational domain + elseif ( size(field,1) == diag_cs%ie-diag_cs%is +2 ) then + isv = 1 ; iev = diag_cs%ie + 2-diag_cs%is ; i_data = .false. ! Symmetric computational domain + else + call MOM_error(FATAL,"post_MOM_IS_data_2d: peculiar size in i-direction of "//trim(diag%name)) + endif + if ( size(field,2) == diag_cs%jed-diag_cs%jsd +1 ) then + jsv = diag_cs%js ; jev = diag_cs%je ; j_data = .true. ! Data domain + elseif ( size(field,2) == diag_cs%jed-diag_cs%jsd +2 ) then + jsv = diag_cs%js ; jev = diag_cs%je+1 ; j_data = .true. ! Symmetric data domain + elseif ( size(field,2) == diag_cs%je-diag_cs%js +1 ) then + jsv = 1 ; jev = diag_cs%je + 1-diag_cs%js ; j_data = .false. ! Computational domain + elseif ( size(field,1) == diag_cs%je-diag_cs%js +2 ) then + jsv = 1 ; jev = diag_cs%je + 2-diag_cs%js ; j_data = .false. ! Symmetric computational domain + else + call MOM_error(FATAL,"post_MOM_IS_data_2d: peculiar size in j-direction "//trim(diag%name)) + endif + + if ((diag%conversion_factor /= 0.) .and. (diag%conversion_factor /= 1.)) then + allocate( locfield( lbound(field,1):ubound(field,1), lbound(field,2):ubound(field,2) ) ) + do j=jsv,jev ; do i=isv,iev + if (field(i,j) == diag_cs%missing_value) then + locfield(i,j) = diag_cs%missing_value + else + locfield(i,j) = field(i,j) * diag%conversion_factor + endif + enddo ; enddo + locfield(isv:iev,jsv:jev) = field(isv:iev,jsv:jev) * diag%conversion_factor + else + locfield => field + endif + + ! Handle cases where the data and computational domain are the same size. + if (diag_cs%ied-diag_cs%isd == diag_cs%ie-diag_cs%is) i_data = j_data + if (diag_cs%jed-diag_cs%jsd == diag_cs%je-diag_cs%js) j_data = i_data + + if (present(mask)) then + if ((size(field,1) /= size(mask,1)) .or. & + (size(field,2) /= size(mask,2))) then + call MOM_error(FATAL, "post_MOM_IS_data_2d: post_MOM_IS_data called with a mask "//& + "that does not match the size of field "//trim(diag%name)) + endif + elseif ( i_data .NEQV. j_data ) then + call MOM_error(FATAL, "post_MOM_IS_data_2d: post_MOM_IS_data called for "//& + trim(diag%name)//" with mixed computational and data domain array sizes.") + endif + + if (is_stat) then + if (present(mask)) then + used = send_data(fms_diag_id, locfield, & + is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, mask=mask) + elseif(i_data .and. associated(diag%mask2d)) then + used = send_data(fms_diag_id, locfield, & + is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, rmask=diag%mask2d) + elseif((.not.i_data) .and. associated(diag%mask2d_comp)) then + used = send_data(fms_diag_id, locfield, & + is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, rmask=diag%mask2d_comp) + else + used = send_data(fms_diag_id, locfield, & + is_in=isv, js_in=jsv, ie_in=iev, je_in=jev) + endif + elseif (diag_cs%ave_enabled) then + if (present(mask)) then + used = send_data(fms_diag_id, locfield, diag_cs%time_end, & + is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, & + weight=diag_cs%time_int, mask=mask) + elseif(i_data .and. associated(diag%mask2d)) then + used = send_data(fms_diag_id, locfield, diag_cs%time_end, & + is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, & + weight=diag_cs%time_int, rmask=diag%mask2d) + elseif((.not.i_data) .and. associated(diag%mask2d_comp)) then + used = send_data(fms_diag_id, locfield, diag_cs%time_end, & + is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, & + weight=diag_cs%time_int, rmask=diag%mask2d_comp) + else + used = send_data(fms_diag_id, locfield, diag_cs%time_end, & + is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, & + weight=diag_cs%time_int) + endif + endif + + if ((diag%conversion_factor /= 0.) .and. (diag%conversion_factor /= 1.) ) deallocate( locfield ) + +end subroutine post_data + + +!> Enable the accumulation of time averages over the specified time interval. +subroutine enable_averaging(time_int_in, time_end_in, diag_cs) + real, intent(in) :: time_int_in !< The time interval over which any values +! !! that are offered are valid [s]. + type(time_type), intent(in) :: time_end_in !< The end time of the valid interval. + type(diag_ctrl), intent(inout) :: diag_cs !< A structure that is used to regulate diagnostic output +! This subroutine enables the accumulation of time averages over the +! specified time interval. + +! if (num_file==0) return + diag_cs%time_int = time_int_in + diag_cs%time_end = time_end_in + diag_cs%ave_enabled = .true. +end subroutine enable_averaging + +! Put a block on averaging any offered fields. +subroutine disable_averaging(diag_cs) + type(diag_ctrl), intent(inout) :: diag_cs !< A structure that is used to regulate diagnostic output + + diag_cs%time_int = 0.0 + diag_cs%ave_enabled = .false. + +end subroutine disable_averaging + +!> Enable the accumulation of time averages over the specified time interval in time units. +subroutine enable_averages(time_int, time_end, diag_CS, T_to_s) + real, intent(in) :: time_int !< The time interval over which any values + !! that are offered are valid [T ~> s]. + type(time_type), intent(in) :: time_end !< The end time of the valid interval. + type(diag_ctrl), intent(inout) :: diag_CS !< A structure that is used to regulate diagnostic output + real, optional, intent(in) :: T_to_s !< A conversion factor for time_int to [s]. +! This subroutine enables the accumulation of time averages over the specified time interval. + + if (present(T_to_s)) then + diag_cs%time_int = time_int*T_to_s +! elseif (associated(diag_CS%US)) then +! diag_cs%time_int = time_int*diag_CS%US%T_to_s + else + diag_cs%time_int = time_int + endif + diag_cs%time_end = time_end + diag_cs%ave_enabled = .true. +end subroutine enable_averages + +!> Indicate whether averaging diagnostics is currently enabled +logical function query_averaging_enabled(diag_cs, time_int, time_end) + type(diag_ctrl), intent(in) :: diag_cs !< A structure that is used to regulate diagnostic output + real, optional, intent(out) :: time_int !< The current setting of diag_cs%time_int [s]. + type(time_type), optional, intent(out) :: time_end !< The current setting of diag_cs%time_end. + + if (present(time_int)) time_int = diag_cs%time_int + if (present(time_end)) time_end = diag_cs%time_end + query_averaging_enabled = diag_cs%ave_enabled +end function query_averaging_enabled + +subroutine diag_mediator_infrastructure_init(err_msg) + ! This subroutine initializes the FMS diag_manager. + character(len=*), optional, intent(out) :: err_msg !< An error message + + call diag_manager_init(err_msg=err_msg) +end subroutine diag_mediator_infrastructure_init + +!> diag_mediator_init initializes the MOM diag_mediator and opens the available + +!> Return the currently specified valid end time for diagnostics +function get_diag_time_end(diag_cs) + type(diag_ctrl), intent(in) :: diag_cs !< A structure that is used to regulate diagnostic output + type(time_type) :: get_diag_time_end + +! This function returns the valid end time for diagnostics that are handled +! outside of the MOM6 infrastructure, such as via the generic tracer code. + + get_diag_time_end = diag_cs%time_end +end function get_diag_time_end + +!> Returns the "MOM_IS_diag_mediator" handle for a group of diagnostics derived from one field. +function register_MOM_IS_diag_field(module_name, field_name, axes, init_time, & + long_name, units, missing_value, range, mask_variant, standard_name, & + verbose, do_not_log, err_msg, interp_method, tile_count, conversion) result (register_diag_field) + integer :: register_diag_field !< The returned diagnostic handle + character(len=*), intent(in) :: module_name !< Name of this module, usually "ice_model" + character(len=*), intent(in) :: field_name !< Name of the diagnostic field + type(axesType), intent(in) :: axes !< The axis group for this field + type(time_type), intent(in) :: init_time !< Time at which a field is first available? + character(len=*), optional, intent(in) :: long_name !< Long name of a field. + character(len=*), optional, intent(in) :: units !< Units of a field. + character(len=*), optional, intent(in) :: standard_name !< Standardized name associated with a field + real, optional, intent(in) :: missing_value !< A value that indicates missing values. + real, optional, intent(in) :: range(2) !< Valid range of a variable (not used in MOM?) + logical, optional, intent(in) :: mask_variant !< If true a logical mask must be provided with + !! post_data calls (not used in MOM?) + logical, optional, intent(in) :: verbose !< If true, FMS is verbose (not used in MOM?) + logical, optional, intent(in) :: do_not_log !< If true, do not log something (not used in MOM?) + character(len=*), optional, intent(out):: err_msg !< String into which an error message might be + !! placed (not used in MOM?) + character(len=*), optional, intent(in) :: interp_method !< If 'none' indicates the field should not + !! be interpolated as a scalar + integer, optional, intent(in) :: tile_count !< no clue (not used in MOM_IS?) + real, optional, intent(in) :: conversion !< A value to multiply data by before writing to file + + ! Local variables + character(len=240) :: mesg + real :: MOM_missing_value + integer :: primary_id, fms_id + type(diag_ctrl), pointer :: diag_cs => NULL() ! A structure that is used + ! to regulate diagnostic output + type(diag_type), pointer :: diag => NULL() + + MOM_missing_value = axes%diag_cs%missing_value + if(present(missing_value)) MOM_missing_value = missing_value + + diag_cs => axes%diag_cs + primary_id = -1 + + fms_id = register_diag_field_fms(module_name, field_name, axes%handles, & + init_time, long_name=long_name, units=units, missing_value=MOM_missing_value, & + range=range, mask_variant=mask_variant, standard_name=standard_name, & + verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, & + interp_method=interp_method, tile_count=tile_count) + if (fms_id > 0) then + primary_id = get_new_diag_id(diag_cs) + diag => diag_cs%diags(primary_id) + diag%fms_diag_id = fms_id + if (len(field_name) > len(diag%name)) then + diag%name = field_name(1:len(diag%name)) + else ; diag%name = field_name ; endif + + if (present(conversion)) diag%conversion_factor = conversion + endif + + if (is_root_pe() .and. diag_CS%doc_unit > 0) then + if (primary_id > 0) then + mesg = '"'//trim(module_name)//'", "'//trim(field_name)//'" [Used]' + else + mesg = '"'//trim(module_name)//'", "'//trim(field_name)//'" [Unused]' + endif + write(diag_CS%doc_unit, '(a)') trim(mesg) + if (present(long_name)) call describe_option("long_name", long_name, diag_CS) + if (present(units)) call describe_option("units", units, diag_CS) + if (present(standard_name)) & + call describe_option("standard_name", standard_name, diag_CS) + endif + + !Decide what mask to use based on the axes info + if (primary_id > 0) then + !3d masks + !2d masks + if (axes%rank == 2) then + diag%mask2d => null() ; diag%mask2d_comp => null() + if (axes%id == diag_cs%axesT1%id) then + diag%mask2d => diag_cs%mask2dT + diag%mask2d_comp => diag_cs%mask2dT_comp + elseif (axes%id == diag_cs%axesB1%id) then + diag%mask2d => diag_cs%mask2dBu + elseif (axes%id == diag_cs%axesCu1%id) then + diag%mask2d => diag_cs%mask2dCu + elseif (axes%id == diag_cs%axesCv1%id) then + diag%mask2d => diag_cs%mask2dCv + ! else + ! call SIS_error(FATAL, "SIS_diag_mediator:register_diag_field: " // & + ! "unknown axes for diagnostic variable "//trim(field_name)) + endif + else + call MOM_error(FATAL, "MOM_IS_diag_mediator:register_diag_field: " // & + "unknown axes for diagnostic variable "//trim(field_name)) + endif + endif ! if (primary_id>-1) + + register_diag_field = primary_id + +end function register_MOM_IS_diag_field + +!> Registers a static diagnostic, returning an integer handle +function register_static_field(module_name, field_name, axes, & + long_name, units, missing_value, range, mask_variant, standard_name, & + do_not_log, interp_method, tile_count) + integer :: register_static_field !< The returned diagnostic handle + character(len=*), intent(in) :: module_name !< Name of this module, usually "ice_model" + character(len=*), intent(in) :: field_name !< Name of the diagnostic field + type(axesType), intent(in) :: axes !< The axis group for this field + character(len=*), optional, intent(in) :: long_name !< Long name of a field. + character(len=*), optional, intent(in) :: units !< Units of a field. + character(len=*), optional, intent(in) :: standard_name !< Standardized name associated with a field + real, optional, intent(in) :: missing_value !< A value that indicates missing values. + real, optional, intent(in) :: range(2) !< Valid range of a variable (not used in MOM?) + logical, optional, intent(in) :: mask_variant !< If true a logical mask must be provided with + !! post_data calls (not used in MOM?) + logical, optional, intent(in) :: do_not_log !< If true, do not log something (not used in MOM?) + character(len=*), optional, intent(in) :: interp_method !< If 'none' indicates the field should not + !! be interpolated as a scalar + integer, optional, intent(in) :: tile_count !< no clue (not used in MOM_IS?) + + ! Local variables + character(len=240) :: mesg + real :: MOM_missing_value + integer :: primary_id, fms_id + type(diag_ctrl), pointer :: diag_cs !< A structure that is used to regulate diagnostic output + + MOM_missing_value = axes%diag_cs%missing_value + if(present(missing_value)) MOM_missing_value = missing_value + + diag_cs => axes%diag_cs + primary_id = -1 + + fms_id = register_static_field_fms(module_name, field_name, axes%handles, & + long_name=long_name, units=units, missing_value=MOM_missing_value, & + range=range, mask_variant=mask_variant, standard_name=standard_name, & + do_not_log=do_not_log, & + interp_method=interp_method, tile_count=tile_count) + if (fms_id > 0) then + primary_id = get_new_diag_id(diag_cs) + diag_cs%diags(primary_id)%fms_diag_id = fms_id + endif + + register_static_field = primary_id + +end function register_static_field + +!> Add a description of an option to the documentation file +subroutine describe_option(opt_name, value, diag_CS) + character(len=*), intent(in) :: opt_name !< The name of the option + character(len=*), intent(in) :: value !< The value of the option + type(diag_ctrl), intent(in) :: diag_CS !< Diagnostic being documented + + ! Local variables + character(len=240) :: mesg + integer :: start_ind = 1, end_ind, len_ind + + len_ind = len_trim(value) + + mesg = " ! "//trim(opt_name)//": "//trim(value) + write(diag_CS%doc_unit, '(a)') trim(mesg) +end subroutine describe_option + +!> Convert the first n elements (up to 3) of an integer array to an underscore delimited string. +function i2s(a, n_in) + integer, dimension(:), intent(in) :: a !< The array of integers to translate + integer, optional , intent(in) :: n_in !< The number of elements to translate, by default all + character(len=15) :: i2s !< The returned string + + ! Local variables + character(len=15) :: i2s_temp + integer :: i,n + + n=size(a) + if(present(n_in)) n = n_in + + i2s = '' + do i=1,n + write (i2s_temp, '(I4.4)') a(i) + i2s = trim(i2s) //'_'// trim(i2s_temp) + enddo + i2s = adjustl(i2s) +end function i2s + +!> Initialize the MOM_IS diag_mediator and opens the available diagnostics file. +subroutine diag_mediator_init(G, param_file, diag_cs, component, err_msg, & + doc_file_dir) + type(ocean_grid_type), intent(inout) :: G !< The horizontal grid type + type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + type(diag_ctrl), intent(inout) :: diag_cs !< A structure that is used to regulate diagnostic output + character(len=*), optional, intent(in) :: component !< An opitonal component name + character(len=*), optional, intent(out) :: err_msg !< A string for a returned error message + 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 + ! is not necessary that the metrics and axis labels be set up yet. + + ! Local variables + integer :: ios, new_unit + logical :: opened, new_file + character(len=8) :: this_pe + character(len=240) :: doc_file, doc_file_dflt, doc_path + character(len=40) :: doc_file_param + character(len=40) :: mdl = "MOM_IS_diag_mediator" ! This module's name. + + call diag_manager_init(err_msg=err_msg) + + ! Allocate list of all diagnostics + allocate(diag_cs%diags(DIAG_ALLOC_CHUNK_SIZE)) + diag_cs%next_free_diag_id = 1 + diag_cs%diags(:)%in_use = .false. + + diag_cs%is = G%isc - (G%isd-1) ; diag_cs%ie = G%iec - (G%isd-1) + diag_cs%js = G%jsc - (G%jsd-1) ; diag_cs%je = G%jec - (G%jsd-1) + diag_cs%isd = G%isd ; diag_cs%ied = G%ied ; diag_cs%jsd = G%jsd ; diag_cs%jed = G%jed + + if (is_root_pe() .and. (diag_CS%doc_unit < 0)) then + if (present(component)) then + doc_file_dflt = trim(component)//".available_diags" + doc_file_param = trim(uppercase(component))//"_AVAILABLE_DIAGS_FILE" + else + write(this_pe,'(i6.6)') PE_here() + doc_file_dflt = "MOM_IS.available_diags."//this_pe + doc_file_param = "AVAILABLE_MOM_IS_DIAGS_FILE" + endif + call get_param(param_file, mdl, trim(doc_file_param), doc_file, & + "A file into which to write a list of all available "//& + "ice shelf diagnostics that can be included in a diag_table.", & + default=doc_file_dflt) + if (len_trim(doc_file) > 0) then + new_file = .true. ; if (diag_CS%doc_unit /= -1) new_file = .false. + ! Find an unused unit number. + do new_unit=512,42,-1 + inquire( new_unit, opened=opened) + if (.not.opened) exit + enddo + + if (opened) call MOM_error(FATAL, & + "diag_mediator_init failed to find an unused unit number.") + + doc_path = doc_file + if (present(doc_file_dir)) then ; if (len_trim(doc_file_dir) > 0) then + doc_path = trim(slasher(doc_file_dir))//trim(doc_file) + endif ; endif + + diag_CS%doc_unit = new_unit + + if (new_file) then + open(diag_CS%doc_unit, file=trim(doc_path), access='SEQUENTIAL', form='FORMATTED', & + action='WRITE', status='REPLACE', iostat=ios) + else ! This file is being reopened, and should be appended. + open(diag_CS%doc_unit, file=trim(doc_path), access='SEQUENTIAL', form='FORMATTED', & + action='WRITE', status='OLD', position='APPEND', iostat=ios) + endif + inquire(diag_CS%doc_unit, opened=opened) + if ((.not.opened) .or. (ios /= 0)) then + call MOM_error(FATAL, "Failed to open available diags file "//trim(doc_path)//".") + endif + endif + endif + + call diag_masks_set(G, -1.0e34, diag_cs) + +end subroutine diag_mediator_init + +subroutine diag_masks_set(G, missing_value, diag_cs) +! Setup the 2d masks for diagnostics + type(ocean_grid_type), target, intent(in) :: G !< The horizontal grid type + real, intent(in) :: missing_value !< A fill value for missing points + type(diag_ctrl), intent(inout) :: diag_cs !< A structure that is used to regulate diagnostic output + + ! Local variables + integer :: i, j, k, NkIce, CatIce + + + diag_cs%mask2dT => G%mask2dT + diag_cs%mask2dBu => G%mask2dBu + diag_cs%mask2dCu => G%mask2dCu + diag_cs%mask2dCv => G%mask2dCv + + allocate(diag_cs%mask2dT_comp(G%isc:G%iec,G%jsc:G%jec)) + do j=G%jsc,G%jec ; do i=G%isc,G%iec + diag_cs%mask2dT_comp(i,j) = diag_cs%mask2dT(i,j) + enddo ; enddo + + + diag_cs%missing_value = missing_value + +end subroutine diag_masks_set + +!> Prevent the registration of additional diagnostics, so that the creation of files can occur +subroutine diag_mediator_close_registration(diag_CS) + type(diag_ctrl), intent(inout) :: diag_CS !< A structure that is used to regulate diagnostic output + + if (diag_CS%doc_unit > -1) then + close(diag_CS%doc_unit) ; diag_CS%doc_unit = -2 + endif + +end subroutine diag_mediator_close_registration + +!> Deallocate memory associated with the MOM_IS diag mediator +subroutine diag_mediator_end(time, diag_CS) + type(time_type), intent(in) :: time !< The current model time + type(diag_ctrl), intent(inout) :: diag_CS !< A structure that is used to regulate diagnostic output + + if (diag_CS%doc_unit > -1) then + close(diag_CS%doc_unit) ; diag_CS%doc_unit = -3 + endif + +end subroutine diag_mediator_end + +!> Allocate a new diagnostic id, noting that it may be necessary to expand the diagnostics array. +function get_new_diag_id(diag_cs) + + integer :: get_new_diag_id !< The returned ID for the new diagnostic + type(diag_ctrl), intent(inout) :: diag_cs !< A structure that is used to regulate diagnostic output + + ! Local variables + type(diag_type), dimension(:), allocatable :: tmp + integer :: i + + if (diag_cs%next_free_diag_id > size(diag_cs%diags)) then + call assert(diag_cs%next_free_diag_id - size(diag_cs%diags) == 1, & + 'get_new_diag_id: inconsistent diag id') + + ! Increase the size of diag_cs%diags and copy data over. + ! Do not use move_alloc() because it is not supported by Fortran 90 + allocate(tmp(size(diag_cs%diags))) + tmp(:) = diag_cs%diags(:) + deallocate(diag_cs%diags) + allocate(diag_cs%diags(size(tmp) + DIAG_ALLOC_CHUNK_SIZE)) + diag_cs%diags(1:size(tmp)) = tmp(:) + deallocate(tmp) + + ! Initialise new part of the diag array. + do i=diag_cs%next_free_diag_id, size(diag_cs%diags) + diag_cs%diags(i)%in_use = .false. + enddo + endif + + get_new_diag_id = diag_cs%next_free_diag_id + diag_cs%next_free_diag_id = diag_cs%next_free_diag_id + 1 + +end function get_new_diag_id + +end module MOM_IS_diag_mediator diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index 0c9fe4e77e..f038190753 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -6,9 +6,9 @@ module MOM_ice_shelf_dynamics use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end use MOM_cpu_clock, only : CLOCK_COMPONENT, CLOCK_ROUTINE -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_averages, disable_averaging +use MOM_IS_diag_mediator, only : post_data, register_diag_field=>register_MOM_IS_diag_field, safe_alloc_ptr +use MOM_IS_diag_mediator, only : diag_mediator_init, set_diag_mediator_grid +use MOM_IS_diag_mediator, only : diag_ctrl, time_type, enable_averages, disable_averaging use MOM_domains, only : MOM_domains_init, clone_MOM_domain use MOM_domains, only : pass_var, pass_vector, TO_ALL, CGRID_NE, BGRID_NE, CORNER use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING, is_root_pe diff --git a/src/ice_shelf/MOM_ice_shelf_initialize.F90 b/src/ice_shelf/MOM_ice_shelf_initialize.F90 index 20479531a8..13eee57539 100644 --- a/src/ice_shelf/MOM_ice_shelf_initialize.F90 +++ b/src/ice_shelf/MOM_ice_shelf_initialize.F90 @@ -4,6 +4,8 @@ module MOM_ice_shelf_initialize ! This file is part of MOM6. See LICENSE.md for the license. use MOM_grid, only : ocean_grid_type +use MOM_array_transform, only : rotate_array +use MOM_hor_index, only : hor_index_type use MOM_file_parser, only : get_param, read_param, log_param, param_file_type use MOM_io, only: MOM_read_data, file_exists, slasher use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING, is_root_pe @@ -25,8 +27,9 @@ module MOM_ice_shelf_initialize contains !> Initialize ice shelf thickness -subroutine initialize_ice_thickness(h_shelf, area_shelf_h, hmask, G, US, PF) +subroutine initialize_ice_thickness(h_shelf, area_shelf_h, hmask, G, G_in, US, PF, rotate_index, turns) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(ocean_grid_type), intent(in) :: G_in !< The ocean's unrotated grid structure real, dimension(SZDI_(G),SZDJ_(G)), & intent(inout) :: h_shelf !< The ice shelf thickness [Z ~> m]. real, dimension(SZDI_(G),SZDJ_(G)), & @@ -36,23 +39,47 @@ subroutine initialize_ice_thickness(h_shelf, area_shelf_h, hmask, G, US, PF) !! partly or fully covered by an ice-shelf type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors type(param_file_type), intent(in) :: PF !< A structure to parse for run-time parameters + logical, intent(in), optional :: rotate_index !< If true, this is a rotation test + integer, intent(in), optional :: turns !< Number of turns for rotation test integer :: i, j character(len=40) :: mdl = "initialize_ice_thickness" ! This subroutine's name. character(len=200) :: config + logical :: rotate = .false. + real, allocatable, dimension(:,:) :: tmp1_2d ! Temporary array for storing ice shelf input data + real, allocatable, dimension(:,:) :: tmp2_2d ! Temporary array for storing ice shelf input data + real, allocatable, dimension(:,:) :: tmp3_2d ! Temporary array for storing ice shelf input data call get_param(PF, mdl, "ICE_PROFILE_CONFIG", config, & "This specifies how the initial ice profile is specified. "//& "Valid values are: CHANNEL, FILE, and USER.", & fail_if_missing=.true.) - select case ( trim(config) ) - case ("CHANNEL"); call initialize_ice_thickness_channel (h_shelf, area_shelf_h, hmask, G, US, PF) - case ("FILE"); call initialize_ice_thickness_from_file (h_shelf, area_shelf_h, hmask, G, US, PF) - case ("USER"); call USER_init_ice_thickness (h_shelf, area_shelf_h, hmask, G, US, PF) - case default ; call MOM_error(FATAL,"MOM_initialize: "// & - "Unrecognized ice profile setup "//trim(config)) - end select + if (PRESENT(rotate_index)) rotate=rotate_index + + if (rotate) then + allocate(tmp1_2d(G_in%isd:G_in%ied,G_in%jsd:G_in%jed));tmp1_2d(:,:)=0.0 + allocate(tmp2_2d(G_in%isd:G_in%ied,G_in%jsd:G_in%jed));tmp2_2d(:,:)=0.0 + allocate(tmp3_2d(G_in%isd:G_in%ied,G_in%jsd:G_in%jed));tmp3_2d(:,:)=0.0 + select case ( trim(config) ) + case ("CHANNEL"); call initialize_ice_thickness_channel (tmp1_2d, tmp2_2d, tmp3_2d, G_in, US, PF) + case ("FILE"); call initialize_ice_thickness_from_file (tmp1_2d, tmp2_2d, tmp3_2d, G_in, US, PF) + case ("USER"); call USER_init_ice_thickness (tmp1_2d, tmp2_2d, tmp3_2d, G_in, US, PF) + case default ; call MOM_error(FATAL,"MOM_initialize: "// & + "Unrecognized ice profile setup "//trim(config)) + end select + call rotate_array(tmp1_2d,turns, h_shelf) + call rotate_array(tmp2_2d,turns, area_shelf_h) + call rotate_array(tmp3_2d,turns, hmask) + else + select case ( trim(config) ) + case ("CHANNEL"); call initialize_ice_thickness_channel (h_shelf, area_shelf_h, hmask, G, US, PF) + case ("FILE"); call initialize_ice_thickness_from_file (h_shelf, area_shelf_h, hmask, G, US, PF) + case ("USER"); call USER_init_ice_thickness (h_shelf, area_shelf_h, hmask, G, US, PF) + case default ; call MOM_error(FATAL,"MOM_initialize: "// & + "Unrecognized ice profile setup "//trim(config)) + end select + endif end subroutine initialize_ice_thickness @@ -77,7 +104,7 @@ subroutine initialize_ice_thickness_from_file(h_shelf, area_shelf_h, hmask, G, U integer :: i, j, isc, jsc, iec, jec real :: len_sidestress, mask, udh - call MOM_mesg(" MOM_ice_shelf_init_profile.F90, initialize_thickness_from_file: reading thickness") + call MOM_mesg("Initialize_ice_thickness_from_file: reading thickness") call get_param(PF, mdl, "INPUTDIR", inputdir, default=".") inputdir = slasher(inputdir) @@ -99,7 +126,6 @@ subroutine initialize_ice_thickness_from_file(h_shelf, area_shelf_h, hmask, G, U if (.not.file_exists(filename, G%Domain)) call MOM_error(FATAL, & " initialize_topography_from_file: Unable to open "//trim(filename)) - call MOM_read_data(filename, trim(thickness_varname), h_shelf, G%Domain, scale=US%m_to_Z) call MOM_read_data(filename,trim(area_varname), area_shelf_h, G%Domain, scale=US%m_to_L**2) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 7972b51fe4..e29687c1c9 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -113,7 +113,7 @@ module MOM_state_initialization !! conditions or by reading them from a restart (or saves) file. subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & restart_CS, ALE_CSp, tracer_Reg, sponge_CSp, & - ALE_sponge_CSp, OBC, Time_in) + ALE_sponge_CSp, OBC, Time_in, frac_shelf_h) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -140,7 +140,9 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & type(ALE_sponge_CS), pointer :: ALE_sponge_CSp !< The ALE sponge control structure. type(ocean_OBC_type), pointer :: OBC !< The open boundary condition control structure. type(time_type), optional, intent(in) :: Time_in !< Time at the start of the run segment. - !! Time_in overrides any value set for Time. + real, dimension(SZI_(G),SZJ_(G)), & + optional, intent(in) :: frac_shelf_h !< The fraction of the grid cell covered + !! by a floating ice shelf [nondim]. ! Local variables character(len=200) :: filename ! The name of an input file. character(len=200) :: filename2 ! The name of an input files. @@ -171,6 +173,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & logical :: debug ! If true, write debugging output. logical :: debug_obc ! If true, do debugging calls related to OBCs. logical :: debug_layers = .false. + logical :: use_ice_shelf character(len=80) :: mesg ! This include declares and sets the variable "version". #include "version_variable.h" @@ -200,6 +203,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & use_EOS = associated(tv%eqn_of_state) use_OBC = associated(OBC) if (use_EOS) eos => tv%eqn_of_state + use_ice_shelf=PRESENT(frac_shelf_h) !==================================================================== ! Initialize temporally evolving fields, either as initial @@ -234,8 +238,8 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & if (.NOT.use_temperature) call MOM_error(FATAL,"MOM_initialize_state : "//& "use_temperature must be true if INIT_LAYERS_FROM_Z_FILE is true") - call MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_params=just_read) - + call MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_params=just_read,& + frac_shelf_h=frac_shelf_h) else ! Initialize thickness, h. call get_param(PF, mdl, "THICKNESS_CONFIG", config, & @@ -1967,7 +1971,7 @@ end subroutine set_velocity_depth_min !> This subroutine determines the isopycnal or other coordinate interfaces and !! layer potential temperatures and salinities directly from a z-space file on !! a latitude-longitude grid. -subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_params) +subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_params, frac_shelf_h) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & intent(out) :: h !< Layer thicknesses being initialized [H ~> m or kg m-2] @@ -1979,6 +1983,9 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param !! to parse for model parameter values. logical, optional, intent(in) :: just_read_params !< If present and true, this call will !! only read parameters without changing h. + real, dimension(SZI_(G),SZJ_(G)), & + optional, intent(in) :: frac_shelf_h !< The fraction of the grid cell covered + !! by a floating ice shelf [nondim]. ! Local variables character(len=200) :: filename !< The name of an input file containing temperature @@ -2012,8 +2019,6 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param real :: eps_Z ! A negligibly thin layer thickness [Z ~> m]. real :: eps_rho ! A negligibly small density difference [R ~> kg m-3]. real :: PI_180 ! for conversion from degrees to radians - - real, dimension(:,:), pointer :: shelf_area => NULL() real :: Hmix_default ! The default initial mixed layer depth [m]. real :: Hmix_depth ! The mixed layer depth in the initial condition [Z ~> m]. real :: dilate ! A dilation factor to match topography [nondim] @@ -2041,8 +2046,6 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param ! Local variables for ALE remapping real, dimension(:), allocatable :: hTarget ! Target thicknesses [Z ~> m]. - real, dimension(:,:), allocatable :: area_shelf_h ! Shelf-covered area per grid cell [L2 ~> m2] - real, dimension(:,:), allocatable, target :: frac_shelf_h ! Fractional shelf area per grid cell [nondim] real, dimension(:,:,:), allocatable, target :: tmpT1dIn, tmpS1dIn real, dimension(:,:,:), allocatable :: tmp_mask_in real, dimension(:,:,:), allocatable :: h1 ! Thicknesses [H ~> m or kg m-2]. @@ -2085,6 +2088,8 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param reentrant_x = .false. ; call get_param(PF, mdl, "REENTRANT_X", reentrant_x, default=.true.) tripolar_n = .false. ; call get_param(PF, mdl, "TRIPOLAR_N", tripolar_n, default=.false.) + use_ice_shelf = present(frac_shelf_h) + call get_param(PF, mdl, "TEMP_SALT_Z_INIT_FILE",filename, & "The name of the z-space input file used to initialize "//& "temperatures (T) and salinities (S). If T and S are not "//& @@ -2144,17 +2149,6 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param "If true, use the order of arithmetic for horizonal regridding that recovers "//& "the answers from the end of 2018. Otherwise, use rotationally symmetric "//& "forms of the same expressions.", default=default_2018_answers) - call get_param(PF, mdl, "ICE_SHELF", use_ice_shelf, default=.false.) - if (use_ice_shelf) then - call get_param(PF, mdl, "ICE_THICKNESS_FILE", ice_shelf_file, & - "The file from which the ice bathymetry and area are read.", & - fail_if_missing=.not.just_read, do_not_log=just_read) - shelf_file = trim(inputdir)//trim(ice_shelf_file) - if (.not.just_read) call log_param(PF, mdl, "INPUTDIR/THICKNESS_FILE", shelf_file) - call get_param(PF, mdl, "ICE_AREA_VARNAME", area_varname, & - "The name of the area variable in ICE_THICKNESS_FILE.", & - fail_if_missing=.not.just_read, do_not_log=just_read) - endif if (.not.useALEremapping) then call get_param(PF, mdl, "ADJUST_THICKNESS", correct_thickness, & "If true, all mass below the bottom removed if the "//& @@ -2217,8 +2211,6 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param do k=1,size(Z_edges_in,1) ; Z_edges_in(k) = -Z_edges_in(k) ; enddo allocate(rho_z(isd:ied,jsd:jed,kd)) - allocate(area_shelf_h(isd:ied,jsd:jed)) - allocate(frac_shelf_h(isd:ied,jsd:jed)) ! Convert T&S to Absolute Salinity and Conservative Temperature if using TEOS10 or NEMO call convert_temp_salt_for_TEOS10(temp_z, salt_z, G%HI, kd, mask_z, eos) @@ -2234,25 +2226,6 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param call pass_var(mask_z,G%Domain) call pass_var(rho_z,G%Domain) - ! This is needed for building an ALE grid under ice shelves - if (use_ice_shelf) then - if (.not.file_exists(shelf_file, G%Domain)) call MOM_error(FATAL, & - "MOM_temp_salt_initialize_from_Z: Unable to open shelf file "//trim(shelf_file)) - - call MOM_read_data(shelf_file, trim(area_varname), area_shelf_h, G%Domain, scale=US%m_to_L**2) - - ! Initialize frac_shelf_h with zeros (open water everywhere) - frac_shelf_h(:,:) = 0.0 - ! Compute fractional ice shelf coverage of h - do j=jsd,jed ; do i=isd,ied - if (G%areaT(i,j) > 0.0) & - frac_shelf_h(i,j) = area_shelf_h(i,j) / G%areaT(i,j) - enddo ; enddo - ! Pass to the pointer for use as an argument to regridding_main - shelf_area => frac_shelf_h - - endif - ! Done with horizontal interpolation. ! Now remap to model coordinates if (useALEremapping) then @@ -2330,7 +2303,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param GV_loc%ke = nkd allocate( dz_interface(isd:ied,jsd:jed,nkd+1) ) ! Need for argument to regridding_main() but is not used if (use_ice_shelf) then - call regridding_main( remapCS, regridCS, G, GV_loc, h1, tv_loc, h, dz_interface, shelf_area ) + call regridding_main( remapCS, regridCS, G, GV_loc, h1, tv_loc, h, dz_interface, frac_shelf_h ) else call regridding_main( remapCS, regridCS, G, GV_loc, h1, tv_loc, h, dz_interface ) endif @@ -2436,7 +2409,8 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param endif deallocate(z_in, z_edges_in, temp_z, salt_z, mask_z) - deallocate(rho_z, area_shelf_h, frac_shelf_h) + deallocate(rho_z) + call pass_var(h, G%Domain) call pass_var(tv%T, G%Domain) diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index 27aa43274b..f66bf95e1b 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -847,11 +847,11 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) ! Build the source grid zTopOfCell = 0. ; zBottomOfCell = 0. ; nPoints = 0; hsrc(:) = 0.0; tmpT1d(:) = -99.9 do k=1,nz_data - if (mask_z(CS%col_i(c),CS%col_j(c),k) == 1.0) then - zBottomOfCell = -min( z_edges_in(k+1), G%bathyT(CS%col_i(c),CS%col_j(c)) ) - tmpT1d(k) = sp_val(CS%col_i(c),CS%col_j(c),k) + if (mask_z(i,j,k) == 1.0) then + zBottomOfCell = -min( z_edges_in(k+1), G%bathyT(i,j) ) + tmpT1d(k) = sp_val(i,j,k) elseif (k>1) then - zBottomOfCell = -G%bathyT(CS%col_i(c),CS%col_j(c)) + zBottomOfCell = -G%bathyT(i,j) tmpT1d(k) = tmpT1d(k-1) else ! This next block should only ever be reached over land tmpT1d(k) = -99.9 @@ -861,9 +861,9 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) zTopOfCell = zBottomOfCell ! Bottom becomes top for next value of k enddo ! In case data is deeper than model - hsrc(nz_data) = hsrc(nz_data) + ( zTopOfCell + G%bathyT(CS%col_i(c),CS%col_j(c)) ) - CS%Ref_val(CS%fldno)%h(1:nz_data,c) = GV%Z_to_H*hsrc(1:nz_data) - CS%Ref_val(CS%fldno)%p(1:nz_data,c) = tmpT1d(1:nz_data) + hsrc(nz_data) = hsrc(nz_data) + ( zTopOfCell + G%bathyT(i,j) ) + CS%Ref_val(m)%h(1:nz_data,c) = GV%Z_to_H*hsrc(1:nz_data) + CS%Ref_val(m)%p(1:nz_data,c) = tmpT1d(1:nz_data) do k=2,nz_data ! if (mask_z(i,j,k)==0.) & if (CS%Ref_val(m)%h(k,c) <= 0.001*GV%m_to_H) & @@ -895,7 +895,7 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) CS%nz, h(i,j,:), tmp_val1, h_neglect, h_neglect_edge) endif !Backward Euler method - CS%var(m)%p(i,j,1:CS%nz) = I1pdamp * (CS%var(m)%p(i,j,1:CS%nz) + tmp_val1 * damp) + CS%var(m)%p(i,j,1:CS%nz) = I1pdamp * (CS%var(m)%p(i,j,1:CS%nz) + tmp_val1(1:CS%nz) * damp) enddo enddo @@ -1185,7 +1185,7 @@ subroutine ALE_sponge_end(CS) if (associated(CS%Iresttime_col_v)) deallocate(CS%Iresttime_col_v) do m=1,CS%fldno - if (associated(CS%Ref_val(CS%fldno)%p)) deallocate(CS%Ref_val(CS%fldno)%p) + if (associated(CS%Ref_val(m)%p)) deallocate(CS%Ref_val(m)%p) enddo deallocate(CS) diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index 7786bf5b46..9def9ae0a1 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -673,6 +673,7 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC) real :: topfn ! A function which goes from 1 at the top to 0 much more ! than Htbl into the interior. real :: z2 ! The distance from the bottom, normalized by Hbbl, nondim. + real :: z2_sq ! z2 squared, used for reproducible evaluation of z2**6. real :: z2_wt ! A nondimensional (0-1) weight used when calculating z2. real :: z_clear ! The clearance of an interface above the surrounding topography [H ~> m or kg m-2]. real :: a_cpl_max ! The maximum drag doefficient across interfaces, set so that it will be @@ -689,6 +690,7 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC) zi_dir ! A trinary logical array indicating which thicknesses to use for ! finding z_clear. integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB ; nz = G%ke @@ -844,7 +846,8 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC) z2_wt = 1.0 ; if (zh(I) * I_HTbl(I) < 2.0*CS%harm_BL_val) & z2_wt = max(0.0, min(1.0, zh(I) * I_HTbl(I) * I_valBL - 1.0)) z2 = z2_wt * (max(zh(I), Ztop_min(I) - min(zcol(i),zcol(i+1))) * I_HTbl(I)) - topfn = 1.0 / (1.0 + 0.09*z2**6) + z2_sq = z2**2 + topfn = 1.0 / (1.0 + 0.09 * (z2_sq * z2_sq * z2_sq)) hvel_shelf(I,k) = min(hvel(I,k), (1.0-topfn)*h_arith(I,k) + topfn*h_harm(I,k)) endif endif @@ -1012,7 +1015,8 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC) z2_wt = 1.0 ; if (zh(i) * I_HTbl(i) < 2.0*CS%harm_BL_val) & z2_wt = max(0.0, min(1.0, zh(i) * I_HTbl(i) * I_valBL - 1.0)) z2 = z2_wt * (max(zh(i), Ztop_min(i) - min(zcol1(i),zcol2(i))) * I_HTbl(i)) - topfn = 1.0 / (1.0 + 0.09*z2**6) + z2_sq = z2**2 + topfn = 1.0 / (1.0 + 0.09 * (z2_sq * z2_sq * z2_sq)) hvel_shelf(i,k) = min(hvel(i,k), (1.0-topfn)*h_arith(i,k) + topfn*h_harm(i,k)) endif endif @@ -1056,7 +1060,6 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC) endif enddo ! end of v-point j loop - if (CS%debug) then call uvchksum("vertvisc_coef h_[uv]", CS%h_u, CS%h_v, G%HI, haloshift=0, & scale=GV%H_to_m, scalar_pair=.true.) @@ -1146,10 +1149,14 @@ subroutine find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, real :: botfn ! A function that is 1 at the bottom and small far from it [nondim] real :: topfn ! A function that is 1 at the top and small far from it [nondim] real :: kv_top ! A viscosity associated with the top boundary layer [Z2 T-1 ~> m2 s-1] + real :: zt_sq ! Square of elements of z_t, used for reproducible evaluation + ! of z_i**6. [H2 ~> m2 or kg2 m-4] or [nondim] logical :: do_shelf, do_OBCs integer :: i, k, is, ie, max_nk integer :: nz + ! testing + a_cpl(:,:) = 0.0 Kv_tot(:,:) = 0.0 @@ -1294,7 +1301,8 @@ subroutine find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, do K=2,nz ; do i=is,ie ; if (do_i(i)) then z_t(i) = z_t(i) + hvel(i,k-1) / tbl_thick(i) - topfn = 1.0 / (1.0 + 0.09 * z_t(i)**6) + zt_sq = z_t(i)**2 + topfn = 1.0 / (1.0 + 0.09 * (zt_sq * zt_sq * zt_sq)) r = 0.5*(hvel(i,k)+hvel(i,k-1)) if (r > tbl_thick(i)) then