From 11d03b13bce3e8a2028d5c8fb80a84d0164fab44 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 21 Mar 2022 09:59:53 -0400 Subject: [PATCH 1/2] +Add unit conversion capability for restarts Added optional arguments to framework code to permit the restart files to work with dimensionally unscaled variables. All answers are bitwise identical, but there are new optional arguments for multiple public interfaces. - Added a new conversion factor element, conv, to the field_restart type to specify how a field is to be unscaled before writing, and whose reciprocal is used to scale restart fields as they are read in. - Added the new optional argument conversion to each of the register_restart_field() and register_restart_pair() routines, which is used to specify how each field will be rescaled when it is written to or read from a restart file. The default is the same as setting this to 1. - Added the new optional argument unscaled to fix_restart_scaling() and fix_restart_unit_scaling() to reset the unit conversion factors that will be save to the restart files to 1, consistent with writing out fields without scaling to the restart files. --- src/core/MOM_verticalGrid.F90 | 9 +- src/framework/MOM_restart.F90 | 178 ++++++++++++++++++----------- src/framework/MOM_unit_scaling.F90 | 13 ++- 3 files changed, 133 insertions(+), 67 deletions(-) diff --git a/src/core/MOM_verticalGrid.F90 b/src/core/MOM_verticalGrid.F90 index a340b5f80f..2df65f09aa 100644 --- a/src/core/MOM_verticalGrid.F90 +++ b/src/core/MOM_verticalGrid.F90 @@ -188,10 +188,17 @@ subroutine verticalGridInit( param_file, GV, US ) end subroutine verticalGridInit !> Set the scaling factors for restart files to the scaling factors for this run. -subroutine fix_restart_scaling(GV) +subroutine fix_restart_scaling(GV, unscaled) type(verticalGrid_type), intent(inout) :: GV !< The ocean's vertical grid structure + logical, optional, intent(in) :: unscaled !< If true, set the restart factors as though the + !! model would be unscaled, which is appropriate if the + !! scaling is undone when writing a restart file. GV%m_to_H_restart = GV%m_to_H + if (present(unscaled)) then ; if (unscaled) then + GV%m_to_H_restart = 1.0 + endif ; endif + end subroutine fix_restart_scaling !> Returns the model's thickness units, usually m or kg/m^2. diff --git a/src/framework/MOM_restart.F90 b/src/framework/MOM_restart.F90 index 1328fd676c..a5386ce7fa 100644 --- a/src/framework/MOM_restart.F90 +++ b/src/framework/MOM_restart.F90 @@ -59,6 +59,10 @@ module MOM_restart !! read from the restart file. logical :: initialized !< .true. if this field has been read from the restart file. character(len=32) :: var_name !< A name by which a variable may be queried. + real :: conv = 1.0 !< A factor by which a restart field should be multiplied before it + !! is written to a restart file, usually to convert it to MKS or + !! other standard units. When read, the restart field is multiplied + !! by the Adcroft reciprocal of this factor. end type field_restart !> A structure to store information about restart fields that are no longer used @@ -146,13 +150,15 @@ subroutine register_restart_field_as_obsolete(field_name, replacement_name, CS) end subroutine register_restart_field_as_obsolete !> Register a 3-d field for restarts, providing the metadata in a structure -subroutine register_restart_field_ptr3d(f_ptr, var_desc, mandatory, CS) +subroutine register_restart_field_ptr3d(f_ptr, var_desc, mandatory, CS, conversion) real, dimension(:,:,:), & target, intent(in) :: f_ptr !< A pointer to the field to be read or written type(vardesc), intent(in) :: var_desc !< A structure with metadata about this variable logical, intent(in) :: mandatory !< If true, the run will abort if this field is not !! successfully read from the restart file. type(MOM_restart_CS), intent(inout) :: CS !< MOM restart control struct + real, optional, intent(in) :: conversion !< A factor to multiply a restart field by + !! before it is written, 1 by default. if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // & "register_restart_field: Module must be initialized before it is used.") @@ -166,6 +172,8 @@ subroutine register_restart_field_ptr3d(f_ptr, var_desc, mandatory, CS) CS%restart_field(CS%novars)%vars = var_desc CS%restart_field(CS%novars)%mand_var = mandatory CS%restart_field(CS%novars)%initialized = .false. + CS%restart_field(CS%novars)%conv = 1.0 + if (present(conversion)) CS%restart_field(CS%novars)%conv = conversion call query_vardesc(CS%restart_field(CS%novars)%vars, & name=CS%restart_field(CS%novars)%var_name, & caller="register_restart_field_ptr3d") @@ -179,13 +187,15 @@ subroutine register_restart_field_ptr3d(f_ptr, var_desc, mandatory, CS) end subroutine register_restart_field_ptr3d !> Register a 4-d field for restarts, providing the metadata in a structure -subroutine register_restart_field_ptr4d(f_ptr, var_desc, mandatory, CS) +subroutine register_restart_field_ptr4d(f_ptr, var_desc, mandatory, CS, conversion) real, dimension(:,:,:,:), & target, intent(in) :: f_ptr !< A pointer to the field to be read or written type(vardesc), intent(in) :: var_desc !< A structure with metadata about this variable logical, intent(in) :: mandatory !< If true, the run will abort if this field is not !! successfully read from the restart file. type(MOM_restart_CS), intent(inout) :: CS !< MOM restart control struct + real, optional, intent(in) :: conversion !< A factor to multiply a restart field by + !! before it is written, 1 by default. if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // & "register_restart_field: Module must be initialized before it is used.") @@ -199,6 +209,8 @@ subroutine register_restart_field_ptr4d(f_ptr, var_desc, mandatory, CS) CS%restart_field(CS%novars)%vars = var_desc CS%restart_field(CS%novars)%mand_var = mandatory CS%restart_field(CS%novars)%initialized = .false. + CS%restart_field(CS%novars)%conv = 1.0 + if (present(conversion)) CS%restart_field(CS%novars)%conv = conversion call query_vardesc(CS%restart_field(CS%novars)%vars, & name=CS%restart_field(CS%novars)%var_name, & caller="register_restart_field_ptr4d") @@ -212,13 +224,15 @@ subroutine register_restart_field_ptr4d(f_ptr, var_desc, mandatory, CS) end subroutine register_restart_field_ptr4d !> Register a 2-d field for restarts, providing the metadata in a structure -subroutine register_restart_field_ptr2d(f_ptr, var_desc, mandatory, CS) +subroutine register_restart_field_ptr2d(f_ptr, var_desc, mandatory, CS, conversion) real, dimension(:,:), & target, intent(in) :: f_ptr !< A pointer to the field to be read or written type(vardesc), intent(in) :: var_desc !< A structure with metadata about this variable logical, intent(in) :: mandatory !< If true, the run will abort if this field is not !! successfully read from the restart file. type(MOM_restart_CS), intent(inout) :: CS !< MOM restart control struct + real, optional, intent(in) :: conversion !< A factor to multiply a restart field by + !! before it is written, 1 by default. if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // & "register_restart_field: Module must be initialized before it is used.") @@ -232,6 +246,8 @@ subroutine register_restart_field_ptr2d(f_ptr, var_desc, mandatory, CS) CS%restart_field(CS%novars)%vars = var_desc CS%restart_field(CS%novars)%mand_var = mandatory CS%restart_field(CS%novars)%initialized = .false. + CS%restart_field(CS%novars)%conv = 1.0 + if (present(conversion)) CS%restart_field(CS%novars)%conv = conversion call query_vardesc(CS%restart_field(CS%novars)%vars, & name=CS%restart_field(CS%novars)%var_name, & caller="register_restart_field_ptr2d") @@ -245,12 +261,14 @@ subroutine register_restart_field_ptr2d(f_ptr, var_desc, mandatory, CS) end subroutine register_restart_field_ptr2d !> Register a 1-d field for restarts, providing the metadata in a structure -subroutine register_restart_field_ptr1d(f_ptr, var_desc, mandatory, CS) +subroutine register_restart_field_ptr1d(f_ptr, var_desc, mandatory, CS, conversion) real, dimension(:), target, intent(in) :: f_ptr !< A pointer to the field to be read or written type(vardesc), intent(in) :: var_desc !< A structure with metadata about this variable logical, intent(in) :: mandatory !< If true, the run will abort if this field is not !! successfully read from the restart file. type(MOM_restart_CS), intent(inout) :: CS !< MOM restart control struct + real, optional, intent(in) :: conversion !< A factor to multiply a restart field by + !! before it is written, 1 by default. if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // & "register_restart_field: Module must be initialized before it is used.") @@ -264,6 +282,8 @@ subroutine register_restart_field_ptr1d(f_ptr, var_desc, mandatory, CS) CS%restart_field(CS%novars)%vars = var_desc CS%restart_field(CS%novars)%mand_var = mandatory CS%restart_field(CS%novars)%initialized = .false. + CS%restart_field(CS%novars)%conv = 1.0 + if (present(conversion)) CS%restart_field(CS%novars)%conv = conversion call query_vardesc(CS%restart_field(CS%novars)%vars, & name=CS%restart_field(CS%novars)%var_name, & caller="register_restart_field_ptr1d") @@ -277,12 +297,14 @@ subroutine register_restart_field_ptr1d(f_ptr, var_desc, mandatory, CS) end subroutine register_restart_field_ptr1d !> Register a 0-d field for restarts, providing the metadata in a structure -subroutine register_restart_field_ptr0d(f_ptr, var_desc, mandatory, CS) +subroutine register_restart_field_ptr0d(f_ptr, var_desc, mandatory, CS, conversion) real, target, intent(in) :: f_ptr !< A pointer to the field to be read or written type(vardesc), intent(in) :: var_desc !< A structure with metadata about this variable logical, intent(in) :: mandatory !< If true, the run will abort if this field is not !! successfully read from the restart file. type(MOM_restart_CS), intent(inout) :: CS !< MOM restart control struct + real, optional, intent(in) :: conversion !< A factor to multiply a restart field by + !! before it is written, 1 by default. if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // & "register_restart_field: Module must be initialized before it is used.") @@ -296,6 +318,8 @@ subroutine register_restart_field_ptr0d(f_ptr, var_desc, mandatory, CS) CS%restart_field(CS%novars)%vars = var_desc CS%restart_field(CS%novars)%mand_var = mandatory CS%restart_field(CS%novars)%initialized = .false. + CS%restart_field(CS%novars)%conv = 1.0 + if (present(conversion)) CS%restart_field(CS%novars)%conv = conversion call query_vardesc(CS%restart_field(CS%novars)%vars, & name=CS%restart_field(CS%novars)%var_name, & caller="register_restart_field_ptr0d") @@ -311,66 +335,72 @@ end subroutine register_restart_field_ptr0d !> Register a pair of rotationally equivalent 2d restart fields subroutine register_restart_pair_ptr2d(a_ptr, b_ptr, a_desc, b_desc, & - mandatory, CS) + mandatory, CS, conversion) real, dimension(:,:), target, intent(in) :: a_ptr !< First field pointer real, dimension(:,:), target, intent(in) :: b_ptr !< Second field pointer - type(vardesc), intent(in) :: a_desc !< First field descriptor - type(vardesc), intent(in) :: b_desc !< Second field descriptor - logical, intent(in) :: mandatory !< If true, abort if field is missing - type(MOM_restart_CS), intent(inout) :: CS !< MOM restart control structure + type(vardesc), intent(in) :: a_desc !< First field descriptor + type(vardesc), intent(in) :: b_desc !< Second field descriptor + logical, intent(in) :: mandatory !< If true, abort if field is missing + type(MOM_restart_CS), intent(inout) :: CS !< MOM restart control structure + real, optional, intent(in) :: conversion !< A factor to multiply a restart field by + !! before it is written, 1 by default. call lock_check(CS, a_desc) if (modulo(CS%turns, 2) /= 0) then - call register_restart_field(b_ptr, a_desc, mandatory, CS) - call register_restart_field(a_ptr, b_desc, mandatory, CS) + call register_restart_field(b_ptr, a_desc, mandatory, CS, conversion) + call register_restart_field(a_ptr, b_desc, mandatory, CS, conversion) else - call register_restart_field(a_ptr, a_desc, mandatory, CS) - call register_restart_field(b_ptr, b_desc, mandatory, CS) + call register_restart_field(a_ptr, a_desc, mandatory, CS, conversion) + call register_restart_field(b_ptr, b_desc, mandatory, CS, conversion) endif end subroutine register_restart_pair_ptr2d !> Register a pair of rotationally equivalent 3d restart fields subroutine register_restart_pair_ptr3d(a_ptr, b_ptr, a_desc, b_desc, & - mandatory, CS) + mandatory, CS, conversion) real, dimension(:,:,:), target, intent(in) :: a_ptr !< First field pointer real, dimension(:,:,:), target, intent(in) :: b_ptr !< Second field pointer - type(vardesc), intent(in) :: a_desc !< First field descriptor - type(vardesc), intent(in) :: b_desc !< Second field descriptor - logical, intent(in) :: mandatory !< If true, abort if field is missing - type(MOM_restart_CS), intent(inout) :: CS !< MOM restart control structure + type(vardesc), intent(in) :: a_desc !< First field descriptor + type(vardesc), intent(in) :: b_desc !< Second field descriptor + logical, intent(in) :: mandatory !< If true, abort if field is missing + type(MOM_restart_CS), intent(inout) :: CS !< MOM restart control structure + real, optional, intent(in) :: conversion !< A factor to multiply a restart field by + !! before it is written, 1 by default. call lock_check(CS, a_desc) if (modulo(CS%turns, 2) /= 0) then - call register_restart_field(b_ptr, a_desc, mandatory, CS) - call register_restart_field(a_ptr, b_desc, mandatory, CS) + call register_restart_field(b_ptr, a_desc, mandatory, CS, conversion) + call register_restart_field(a_ptr, b_desc, mandatory, CS, conversion) else - call register_restart_field(a_ptr, a_desc, mandatory, CS) - call register_restart_field(b_ptr, b_desc, mandatory, CS) + call register_restart_field(a_ptr, a_desc, mandatory, CS, conversion) + call register_restart_field(b_ptr, b_desc, mandatory, CS, conversion) endif end subroutine register_restart_pair_ptr3d !> Register a pair of rotationally equivalent 2d restart fields subroutine register_restart_pair_ptr4d(a_ptr, b_ptr, a_desc, b_desc, & - mandatory, CS) + mandatory, CS, conversion) real, dimension(:,:,:,:), target, intent(in) :: a_ptr !< First field pointer real, dimension(:,:,:,:), target, intent(in) :: b_ptr !< Second field pointer - type(vardesc), intent(in) :: a_desc !< First field descriptor - type(vardesc), intent(in) :: b_desc !< Second field descriptor - logical, intent(in) :: mandatory !< If true, abort if field is missing - type(MOM_restart_CS), intent(inout) :: CS !< MOM restart control structure + type(vardesc), intent(in) :: a_desc !< First field descriptor + type(vardesc), intent(in) :: b_desc !< Second field descriptor + logical, intent(in) :: mandatory !< If true, abort if field is missing + type(MOM_restart_CS), intent(inout) :: CS !< MOM restart control structure + real, optional, intent(in) :: conversion !< A factor to multiply a restart field by + !! before it is written, 1 by default. call lock_check(CS, a_desc) if (modulo(CS%turns, 2) /= 0) then - call register_restart_field(b_ptr, a_desc, mandatory, CS) - call register_restart_field(a_ptr, b_desc, mandatory, CS) + call register_restart_field(b_ptr, a_desc, mandatory, CS, conversion) + call register_restart_field(a_ptr, b_desc, mandatory, CS, conversion) else - call register_restart_field(a_ptr, a_desc, mandatory, CS) - call register_restart_field(b_ptr, b_desc, mandatory, CS) + call register_restart_field(a_ptr, a_desc, mandatory, CS, conversion) + call register_restart_field(b_ptr, b_desc, mandatory, CS, conversion) endif end subroutine register_restart_pair_ptr4d @@ -378,7 +408,7 @@ end subroutine register_restart_pair_ptr4d ! The following provide alternate interfaces to register restarts. !> Register a 4-d field for restarts, providing the metadata as individual arguments -subroutine register_restart_field_4d(f_ptr, name, mandatory, CS, longname, units, & +subroutine register_restart_field_4d(f_ptr, name, mandatory, CS, longname, units, conversion, & hor_grid, z_grid, t_grid) real, dimension(:,:,:,:), & target, intent(in) :: f_ptr !< A pointer to the field to be read or written @@ -388,6 +418,8 @@ subroutine register_restart_field_4d(f_ptr, name, mandatory, CS, longname, units type(MOM_restart_CS), intent(inout) :: CS !< MOM restart control struct character(len=*), optional, intent(in) :: longname !< variable long name character(len=*), optional, intent(in) :: units !< variable units + real, optional, intent(in) :: conversion !< A factor to multiply a restart field by + !! before it is written, 1 by default. character(len=*), optional, intent(in) :: hor_grid !< variable horizontal staggering, 'h' if absent character(len=*), optional, intent(in) :: z_grid !< variable vertical staggering, 'L' if absent character(len=*), optional, intent(in) :: t_grid !< time description: s, p, or 1, 's' if absent @@ -403,12 +435,12 @@ subroutine register_restart_field_4d(f_ptr, name, mandatory, CS, longname, units vd = var_desc(name, units=units, longname=longname, hor_grid=hor_grid, & z_grid=z_grid, t_grid=t_grid) - call register_restart_field_ptr4d(f_ptr, vd, mandatory, CS) + call register_restart_field_ptr4d(f_ptr, vd, mandatory, CS, conversion) end subroutine register_restart_field_4d !> Register a 3-d field for restarts, providing the metadata as individual arguments -subroutine register_restart_field_3d(f_ptr, name, mandatory, CS, longname, units, & +subroutine register_restart_field_3d(f_ptr, name, mandatory, CS, longname, units, conversion, & hor_grid, z_grid, t_grid) real, dimension(:,:,:), & target, intent(in) :: f_ptr !< A pointer to the field to be read or written @@ -418,6 +450,8 @@ subroutine register_restart_field_3d(f_ptr, name, mandatory, CS, longname, units type(MOM_restart_CS), intent(inout) :: CS !< MOM restart control struct character(len=*), optional, intent(in) :: longname !< variable long name character(len=*), optional, intent(in) :: units !< variable units + real, optional, intent(in) :: conversion !< A factor to multiply a restart field by + !! before it is written, 1 by default. character(len=*), optional, intent(in) :: hor_grid !< variable horizontal staggering, 'h' if absent character(len=*), optional, intent(in) :: z_grid !< variable vertical staggering, 'L' if absent character(len=*), optional, intent(in) :: t_grid !< time description: s, p, or 1, 's' if absent @@ -433,12 +467,12 @@ subroutine register_restart_field_3d(f_ptr, name, mandatory, CS, longname, units vd = var_desc(name, units=units, longname=longname, hor_grid=hor_grid, & z_grid=z_grid, t_grid=t_grid) - call register_restart_field_ptr3d(f_ptr, vd, mandatory, CS) + call register_restart_field_ptr3d(f_ptr, vd, mandatory, CS, conversion) end subroutine register_restart_field_3d !> Register a 2-d field for restarts, providing the metadata as individual arguments -subroutine register_restart_field_2d(f_ptr, name, mandatory, CS, longname, units, & +subroutine register_restart_field_2d(f_ptr, name, mandatory, CS, longname, units, conversion, & hor_grid, z_grid, t_grid) real, dimension(:,:), & target, intent(in) :: f_ptr !< A pointer to the field to be read or written @@ -448,6 +482,8 @@ subroutine register_restart_field_2d(f_ptr, name, mandatory, CS, longname, units type(MOM_restart_CS), intent(inout) :: CS !< MOM restart control struct character(len=*), optional, intent(in) :: longname !< variable long name character(len=*), optional, intent(in) :: units !< variable units + real, optional, intent(in) :: conversion !< A factor to multiply a restart field by + !! before it is written, 1 by default. character(len=*), optional, intent(in) :: hor_grid !< variable horizontal staggering, 'h' if absent character(len=*), optional, intent(in) :: z_grid !< variable vertical staggering, '1' if absent character(len=*), optional, intent(in) :: t_grid !< time description: s, p, or 1, 's' if absent @@ -466,12 +502,12 @@ subroutine register_restart_field_2d(f_ptr, name, mandatory, CS, longname, units vd = var_desc(name, units=units, longname=longname, hor_grid=hor_grid, & z_grid=zgrid, t_grid=t_grid) - call register_restart_field_ptr2d(f_ptr, vd, mandatory, CS) + call register_restart_field_ptr2d(f_ptr, vd, mandatory, CS, conversion) end subroutine register_restart_field_2d !> Register a 1-d field for restarts, providing the metadata as individual arguments -subroutine register_restart_field_1d(f_ptr, name, mandatory, CS, longname, units, & +subroutine register_restart_field_1d(f_ptr, name, mandatory, CS, longname, units, conversion, & hor_grid, z_grid, t_grid) real, dimension(:), target, intent(in) :: f_ptr !< A pointer to the field to be read or written character(len=*), intent(in) :: name !< variable name to be used in the restart file @@ -480,6 +516,8 @@ subroutine register_restart_field_1d(f_ptr, name, mandatory, CS, longname, units type(MOM_restart_CS), intent(inout) :: CS !< MOM restart control struct character(len=*), optional, intent(in) :: longname !< variable long name character(len=*), optional, intent(in) :: units !< variable units + real, optional, intent(in) :: conversion !< A factor to multiply a restart field by + !! before it is written, 1 by default. character(len=*), optional, intent(in) :: hor_grid !< variable horizontal staggering, '1' if absent character(len=*), optional, intent(in) :: z_grid !< variable vertical staggering, 'L' if absent character(len=*), optional, intent(in) :: t_grid !< time description: s, p, or 1, 's' if absent @@ -498,12 +536,12 @@ subroutine register_restart_field_1d(f_ptr, name, mandatory, CS, longname, units vd = var_desc(name, units=units, longname=longname, hor_grid=hgrid, & z_grid=z_grid, t_grid=t_grid) - call register_restart_field_ptr1d(f_ptr, vd, mandatory, CS) + call register_restart_field_ptr1d(f_ptr, vd, mandatory, CS, conversion) end subroutine register_restart_field_1d !> Register a 0-d field for restarts, providing the metadata as individual arguments -subroutine register_restart_field_0d(f_ptr, name, mandatory, CS, longname, units, & +subroutine register_restart_field_0d(f_ptr, name, mandatory, CS, longname, units, conversion, & t_grid) real, target, intent(in) :: f_ptr !< A pointer to the field to be read or written character(len=*), intent(in) :: name !< variable name to be used in the restart file @@ -512,6 +550,8 @@ subroutine register_restart_field_0d(f_ptr, name, mandatory, CS, longname, units type(MOM_restart_CS), intent(inout) :: CS !< MOM restart control struct character(len=*), optional, intent(in) :: longname !< variable long name character(len=*), optional, intent(in) :: units !< variable units + real, optional, intent(in) :: conversion !< A factor to multiply a restart field by + !! before it is written, 1 by default. character(len=*), optional, intent(in) :: t_grid !< time description: s, p, or 1, 's' if absent type(vardesc) :: vd @@ -525,7 +565,7 @@ subroutine register_restart_field_0d(f_ptr, name, mandatory, CS, longname, units vd = var_desc(name, units=units, longname=longname, hor_grid='1', & z_grid='1', t_grid=t_grid) - call register_restart_field_ptr0d(f_ptr, vd, mandatory, CS) + call register_restart_field_ptr0d(f_ptr, vd, mandatory, CS, conversion) end subroutine register_restart_field_0d @@ -910,6 +950,7 @@ subroutine save_restart(directory, time, G, CS, time_stamped, filename, GV, num_ integer :: seconds, days, year, month, hour, minute character(len=8) :: hor_grid, z_grid, t_grid ! Variable grid info. character(len=64) :: var_name ! A variable's name. + real :: conv ! Shorthand for the conversion factor real :: restart_time character(len=32) :: filename_appendix = '' ! Appendix to filename for ensemble runs integer :: length ! The length of a text string. @@ -1025,16 +1066,17 @@ subroutine save_restart(directory, time, G, CS, time_stamped, filename, GV, num_ call get_checksum_loop_ranges(G, pos, isL, ieL, jsL, jeL) endif do m=start_var,next_var-1 + conv = CS%restart_field(m)%conv if (associated(CS%var_ptr3d(m)%p)) then - check_val(m-start_var+1,1) = chksum(CS%var_ptr3d(m)%p(isL:ieL,jsL:jeL,:), turns=-turns) + check_val(m-start_var+1,1) = chksum(conv*CS%var_ptr3d(m)%p(isL:ieL,jsL:jeL,:), turns=-turns) elseif (associated(CS%var_ptr2d(m)%p)) then - check_val(m-start_var+1,1) = chksum(CS%var_ptr2d(m)%p(isL:ieL,jsL:jeL), turns=-turns) + check_val(m-start_var+1,1) = chksum(conv*CS%var_ptr2d(m)%p(isL:ieL,jsL:jeL), turns=-turns) elseif (associated(CS%var_ptr4d(m)%p)) then - check_val(m-start_var+1,1) = chksum(CS%var_ptr4d(m)%p(isL:ieL,jsL:jeL,:,:), turns=-turns) + check_val(m-start_var+1,1) = chksum(conv*CS%var_ptr4d(m)%p(isL:ieL,jsL:jeL,:,:), turns=-turns) elseif (associated(CS%var_ptr1d(m)%p)) then - check_val(m-start_var+1,1) = chksum(CS%var_ptr1d(m)%p) + check_val(m-start_var+1,1) = chksum(conv*CS%var_ptr1d(m)%p(:)) elseif (associated(CS%var_ptr0d(m)%p)) then - check_val(m-start_var+1,1) = chksum(CS%var_ptr0d(m)%p, pelist=(/PE_here()/)) + check_val(m-start_var+1,1) = chksum(conv*CS%var_ptr0d(m)%p, pelist=(/PE_here()/)) endif enddo @@ -1048,18 +1090,20 @@ subroutine save_restart(directory, time, G, CS, time_stamped, filename, GV, num_ do m=start_var,next_var-1 if (associated(CS%var_ptr3d(m)%p)) then - call MOM_write_field(IO_handle, fields(m-start_var+1), G%Domain, & - CS%var_ptr3d(m)%p, restart_time, turns=-turns) + call MOM_write_field(IO_handle, fields(m-start_var+1), G%Domain, CS%var_ptr3d(m)%p, & + restart_time, scale=CS%restart_field(m)%conv, turns=-turns) elseif (associated(CS%var_ptr2d(m)%p)) then - call MOM_write_field(IO_handle, fields(m-start_var+1), G%Domain, & - CS%var_ptr2d(m)%p, restart_time, turns=-turns) + call MOM_write_field(IO_handle, fields(m-start_var+1), G%Domain, CS%var_ptr2d(m)%p, & + restart_time, scale=CS%restart_field(m)%conv, turns=-turns) elseif (associated(CS%var_ptr4d(m)%p)) then - call MOM_write_field(IO_handle, fields(m-start_var+1), G%Domain, & - CS%var_ptr4d(m)%p, restart_time, turns=-turns) + call MOM_write_field(IO_handle, fields(m-start_var+1), G%Domain, CS%var_ptr4d(m)%p, & + restart_time, scale=CS%restart_field(m)%conv, turns=-turns) elseif (associated(CS%var_ptr1d(m)%p)) then - call MOM_write_field(IO_handle, fields(m-start_var+1), CS%var_ptr1d(m)%p, restart_time) + call MOM_write_field(IO_handle, fields(m-start_var+1), CS%var_ptr1d(m)%p, & + restart_time, scale=CS%restart_field(m)%conv) elseif (associated(CS%var_ptr0d(m)%p)) then - call MOM_write_field(IO_handle, fields(m-start_var+1), CS%var_ptr0d(m)%p, restart_time) + call MOM_write_field(IO_handle, fields(m-start_var+1), CS%var_ptr0d(m)%p, & + restart_time, scale=CS%restart_field(m)%conv) endif enddo @@ -1085,6 +1129,8 @@ subroutine restore_state(filename, directory, day, G, CS) type(MOM_restart_CS), intent(inout) :: CS !< MOM restart control struct ! Local variables + real :: scale ! A scaling factor for reading a field + real :: conv ! The output conversion factor for writing a field character(len=200) :: filepath ! The path (dir/file) to the file being opened. character(len=80) :: fname ! The name of the current file. character(len=8) :: suffix ! A suffix (like "_2") that is added to any @@ -1198,6 +1244,8 @@ subroutine restore_state(filename, directory, day, G, CS) case ('1') ; pos = 0 case default ; pos = 0 end select + conv = CS%restart_field(m)%conv + if (conv == 0.0) then ; scale = 1.0 ; else ; scale = 1.0 / conv ; endif call get_checksum_loop_ranges(G, pos, isL, ieL, jsL, jeL) do i=1, nvar @@ -1214,42 +1262,42 @@ subroutine restore_state(filename, directory, day, G, CS) if (associated(CS%var_ptr1d(m)%p)) then ! Read a 1d array, which should be invariant to domain decomposition. call MOM_read_data(unit_path(n), varname, CS%var_ptr1d(m)%p, & - timelevel=1, MOM_Domain=G%Domain) - if (is_there_a_checksum) checksum_data = chksum(CS%var_ptr1d(m)%p) + timelevel=1, scale=scale, MOM_Domain=G%Domain) + if (is_there_a_checksum) checksum_data = chksum(conv*CS%var_ptr1d(m)%p(:)) elseif (associated(CS%var_ptr0d(m)%p)) then ! Read a scalar... call MOM_read_data(unit_path(n), varname, CS%var_ptr0d(m)%p, & - timelevel=1, MOM_Domain=G%Domain) - if (is_there_a_checksum) checksum_data = chksum(CS%var_ptr0d(m)%p, pelist=(/PE_here()/)) + timelevel=1, scale=scale, MOM_Domain=G%Domain) + if (is_there_a_checksum) checksum_data = chksum(conv*CS%var_ptr0d(m)%p, pelist=(/PE_here()/)) elseif (associated(CS%var_ptr2d(m)%p)) then ! Read a 2d array. if (pos /= 0) then call MOM_read_data(unit_path(n), varname, CS%var_ptr2d(m)%p, & - G%Domain, timelevel=1, position=pos) + G%Domain, timelevel=1, position=pos, scale=scale) else ! This array is not domain-decomposed. This variant may be under-tested. call MOM_error(FATAL, & "MOM_restart does not support 2-d arrays without domain decomposition.") ! call read_data(unit_path(n), varname, CS%var_ptr2d(m)%p,no_domain=.true., timelevel=1) endif - if (is_there_a_checksum) checksum_data = chksum(CS%var_ptr2d(m)%p(isL:ieL,jsL:jeL)) + if (is_there_a_checksum) checksum_data = chksum(conv*CS%var_ptr2d(m)%p(isL:ieL,jsL:jeL)) elseif (associated(CS%var_ptr3d(m)%p)) then ! Read a 3d array. if (pos /= 0) then call MOM_read_data(unit_path(n), varname, CS%var_ptr3d(m)%p, & - G%Domain, timelevel=1, position=pos) + G%Domain, timelevel=1, position=pos, scale=scale) else ! This array is not domain-decomposed. This variant may be under-tested. call MOM_error(FATAL, & "MOM_restart does not support 3-d arrays without domain decomposition.") ! call read_data(unit_path(n), varname, CS%var_ptr3d(m)%p, no_domain=.true., timelevel=1) endif - if (is_there_a_checksum) checksum_data = chksum(CS%var_ptr3d(m)%p(isL:ieL,jsL:jeL,:)) + if (is_there_a_checksum) checksum_data = chksum(conv*CS%var_ptr3d(m)%p(isL:ieL,jsL:jeL,:)) elseif (associated(CS%var_ptr4d(m)%p)) then ! Read a 4d array. if (pos /= 0) then call MOM_read_data(unit_path(n), varname, CS%var_ptr4d(m)%p, & - G%Domain, timelevel=1, position=pos) + G%Domain, timelevel=1, position=pos, scale=scale) else ! This array is not domain-decomposed. This variant may be under-tested. call MOM_error(FATAL, & "MOM_restart does not support 4-d arrays without domain decomposition.") ! call read_data(unit_path(n), varname, CS%var_ptr4d(m)%p, no_domain=.true., timelevel=1) endif - if (is_there_a_checksum) checksum_data = chksum(CS%var_ptr4d(m)%p(isL:ieL,jsL:jeL,:,:)) + if (is_there_a_checksum) checksum_data = chksum(conv*CS%var_ptr4d(m)%p(isL:ieL,jsL:jeL,:,:)) else call MOM_error(FATAL, "MOM_restart restore_state: No pointers set for "//trim(varname)) endif diff --git a/src/framework/MOM_unit_scaling.F90 b/src/framework/MOM_unit_scaling.F90 index cd339f410c..bf8fd24b44 100644 --- a/src/framework/MOM_unit_scaling.F90 +++ b/src/framework/MOM_unit_scaling.F90 @@ -196,8 +196,11 @@ end subroutine set_unit_scaling_combos !> Set the unit scaling factors for output to restart files to the unit scaling !! factors for this run. -subroutine fix_restart_unit_scaling(US) +subroutine fix_restart_unit_scaling(US, unscaled) type(unit_scale_type), intent(inout) :: US !< A dimensional unit scaling type + logical, optional, intent(in) :: unscaled !< If true, set the restart factors as though the + !! model would be unscaled, which is appropriate if the + !! scaling is undone when writing a restart file. US%m_to_Z_restart = US%m_to_Z US%m_to_L_restart = US%m_to_L @@ -205,6 +208,14 @@ subroutine fix_restart_unit_scaling(US) US%kg_m3_to_R_restart = US%kg_m3_to_R US%J_kg_to_Q_restart = US%J_kg_to_Q + if (present(unscaled)) then ; if (unscaled) then + US%m_to_Z_restart = 1.0 + US%m_to_L_restart = 1.0 + US%s_to_T_restart = 1.0 + US%kg_m3_to_R_restart = 1.0 + US%J_kg_to_Q_restart = 1.0 + endif ; endif + end subroutine fix_restart_unit_scaling !> Deallocates a unit scaling structure. From 7c3f31bd473da8c0a9fef327bf14b0ba56aadc95 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 21 Mar 2022 10:32:36 -0400 Subject: [PATCH 2/2] +(*)Write unscaled MOM6 restart files Revised the MOM6 code so that the output restart files are always exactly the same as they would be if no dimensional rescaling is applied. The input restart files can still have rescaling, so files written by previous versions of the code still work exactly as before. This does change the output, in the sense that the restart files are unscaled and some units documents in the restart files are corrected, but the solutions themselves are bitwise identical. Also, there are new (non-optional) unit scaling type arguments to several routines. The specific changes in this commit include: - Added conversion factors to all register_restart_field() or register_restart_pair() call for variables that are subject to dimensional rescaling. - Revised the calculation of restart rescaling factors to reflect the rescaling that now occurs directly in the restore_state() call. - Added new US arguments to register_restarts_dyn_split_RK2(), MEKE_alloc_register_restart(), set_visc_register_restarts() - Added a new GV argument to mixedlayer_restrat_register_restarts() - Used the new unscaled argument in calls to fix_restart_scaling() and fix_restart_unit_scaling() - Revised several calls to register_restart_field() to avoid using the variant that uses a var_desc type, eliminating the need for some modules to reference the vardesc type or the routine var_desc(). - Added or corrected unit descriptions to the comments describing a few variables. - Noted a few probably bugs in comments with ###, but did not fix them. --- src/core/MOM.F90 | 55 +++++++--------- src/core/MOM_barotropic.F90 | 26 ++++---- src/core/MOM_dynamics_split_RK2.F90 | 66 +++++++++---------- src/core/MOM_open_boundary.F90 | 64 ++++++++++-------- src/ice_shelf/MOM_ice_shelf.F90 | 24 +++---- src/ice_shelf/MOM_ice_shelf_dynamics.F90 | 27 ++++---- src/ice_shelf/MOM_ice_shelf_initialize.F90 | 49 +++++++------- .../MOM_state_initialization.F90 | 11 ++-- src/ocean_data_assim/MOM_oda_incupd.F90 | 52 +++++++-------- src/parameterizations/lateral/MOM_MEKE.F90 | 65 +++++++++--------- .../lateral/MOM_internal_tides.F90 | 10 ++- .../lateral/MOM_mixed_layer_restrat.F90 | 28 ++++---- .../vertical/MOM_set_viscosity.F90 | 25 +++---- src/tracer/boundary_impulse_tracer.F90 | 6 +- src/user/MOM_controlled_forcing.F90 | 36 +++++----- 15 files changed, 280 insertions(+), 264 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 8f6dab2a1a..5b35c01a42 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -2422,7 +2422,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & call restart_init(param_file, restart_CSp) call set_restart_fields(GV, US, param_file, CS, restart_CSp) if (CS%split) then - call register_restarts_dyn_split_RK2(HI, GV, param_file, & + call register_restarts_dyn_split_RK2(HI, GV, US, param_file, & CS%dyn_split_RK2_CSp, restart_CSp, CS%uh, CS%vh) elseif (CS%use_RK2) then call register_restarts_dyn_unsplit_RK2(HI, GV, param_file, & @@ -2437,9 +2437,9 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & call call_tracer_register(HI, GV, US, param_file, CS%tracer_flow_CSp, & CS%tracer_Reg, restart_CSp) - call MEKE_alloc_register_restart(HI, param_file, CS%MEKE, restart_CSp) - call set_visc_register_restarts(HI, GV, param_file, CS%visc, restart_CSp) - call mixedlayer_restrat_register_restarts(HI, param_file, & + call MEKE_alloc_register_restart(HI, US, param_file, CS%MEKE, restart_CSp) + call set_visc_register_restarts(HI, GV, US, param_file, CS%visc, restart_CSp) + call mixedlayer_restrat_register_restarts(HI, GV, param_file, & CS%mixedlayer_restrat_CSp, restart_CSp) if (CS%rotate_index .and. associated(OBC_in) .and. use_temperature) then @@ -2468,7 +2468,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & ! This needs the number of tracers and to have called any code that sets whether ! reservoirs are used. - call open_boundary_register_restarts(HI, GV, CS%OBC, CS%tracer_Reg, & + call open_boundary_register_restarts(HI, GV, US, CS%OBC, CS%tracer_Reg, & param_file, restart_CSp, use_temperature) endif @@ -2865,11 +2865,9 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & if (use_frazil) then if (query_initialized(CS%tv%frazil, "frazil", restart_CSp)) then ! Test whether the dimensional rescaling has changed for heat content. - if ((US%kg_m3_to_R_restart*US%m_to_Z_restart*US%J_kg_to_Q_restart /= 0.0) .and. & - ((US%J_kg_to_Q*US%kg_m3_to_R*US%m_to_Z) /= & - (US%J_kg_to_Q_restart*US%kg_m3_to_R_restart*US%m_to_Z_restart)) ) then - QRZ_rescale = (US%J_kg_to_Q*US%kg_m3_to_R*US%m_to_Z) / & - (US%J_kg_to_Q_restart*US%kg_m3_to_R_restart*US%m_to_Z_restart) + if ((US%J_kg_to_Q_restart*US%kg_m3_to_R_restart*US%m_to_Z_restart /= 0.0) .and. & + (US%J_kg_to_Q_restart*US%kg_m3_to_R_restart*US%m_to_Z_restart /= 1.0) ) then + QRZ_rescale = 1.0 / (US%J_kg_to_Q_restart*US%kg_m3_to_R_restart*US%m_to_Z_restart) do j=js,je ; do i=is,ie CS%tv%frazil(i,j) = QRZ_rescale * CS%tv%frazil(i,j) enddo ; enddo @@ -2885,10 +2883,8 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & if (CS%p_surf_prev_set) then ! Test whether the dimensional rescaling has changed for pressure. if ((US%kg_m3_to_R_restart*US%s_to_T_restart*US%m_to_L_restart /= 0.0) .and. & - ((US%kg_m3_to_R*(US%m_to_L*US%s_to_T_restart)**2) /= & - (US%kg_m3_to_R_restart*(US%m_to_L_restart*US%s_to_T)**2)) ) then - RL2_T2_rescale = (US%kg_m3_to_R*(US%m_to_L*US%s_to_T_restart)**2) / & - (US%kg_m3_to_R_restart*(US%m_to_L_restart*US%s_to_T)**2) + (US%s_to_T_restart**2 /= US%kg_m3_to_R_restart * US%m_to_L_restart**2) ) then + RL2_T2_rescale = US%s_to_T_restart**2 / (US%kg_m3_to_R_restart*US%m_to_L_restart**2) do j=js,je ; do i=is,ie CS%p_surf_prev(i,j) = RL2_T2_rescale * CS%p_surf_prev(i,j) enddo ; enddo @@ -2901,8 +2897,8 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & if (use_ice_shelf .and. associated(CS%Hml)) then if (query_initialized(CS%Hml, "hML", restart_CSp)) then ! Test whether the dimensional rescaling has changed for depths. - if ((US%m_to_Z_restart /= 0.0) .and. (US%m_to_Z /= US%m_to_Z_restart) ) then - Z_rescale = US%m_to_Z / US%m_to_Z_restart + if ((US%m_to_Z_restart /= 0.0) .and. (US%m_to_Z_restart /= 1.0) ) then + Z_rescale = 1.0 / US%m_to_Z_restart do j=js,je ; do i=is,ie CS%Hml(i,j) = Z_rescale * CS%Hml(i,j) enddo ; enddo @@ -2911,8 +2907,8 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & endif if (query_initialized(CS%ave_ssh_ibc,"ave_ssh",restart_CSp)) then - if ((US%m_to_Z_restart /= 0.0) .and. (US%m_to_Z /= US%m_to_Z_restart) ) then - Z_rescale = US%m_to_Z / US%m_to_Z_restart + if ((US%m_to_Z_restart /= 0.0) .and. (US%m_to_Z_restart /= 1.0) ) then + Z_rescale = 1.0 / US%m_to_Z_restart do j=js,je ; do i=is,ie CS%ave_ssh_ibc(i,j) = Z_rescale * CS%ave_ssh_ibc(i,j) enddo ; enddo @@ -2967,7 +2963,6 @@ subroutine finish_MOM_initialization(Time, dirs, CS, restart_CSp) ! various unit conversion factors type(MOM_restart_CS), pointer :: restart_CSp_tmp => NULL() real, allocatable :: z_interface(:,:,:) ! Interface heights [m] - type(vardesc) :: vd call cpu_clock_begin(id_clock_init) call callTree_enter("finish_MOM_initialization()") @@ -2976,8 +2971,8 @@ subroutine finish_MOM_initialization(Time, dirs, CS, restart_CSp) G => CS%G ; GV => CS%GV ; US => CS%US !### Move to initialize_MOM? - call fix_restart_scaling(GV) - call fix_restart_unit_scaling(US) + call fix_restart_scaling(GV, unscaled=.true.) + call fix_restart_unit_scaling(US, unscaled=.true.) if (CS%use_particles) then @@ -3090,9 +3085,6 @@ subroutine set_restart_fields(GV, US, param_file, CS, restart_CSp) thickness_units = get_thickness_units(GV) flux_units = get_flux_units(GV) - u_desc = var_desc("u", "m s-1", "Zonal velocity", hor_grid='Cu') - v_desc = var_desc("v", "m s-1", "Meridional velocity", hor_grid='Cv') - if (associated(CS%tv%T)) & call register_restart_field(CS%tv%T, "Temp", .true., restart_CSp, & "Potential Temperature", "degC") @@ -3101,28 +3093,31 @@ subroutine set_restart_fields(GV, US, param_file, CS, restart_CSp) "Salinity", "PPT") call register_restart_field(CS%h, "h", .true., restart_CSp, & - "Layer Thickness", thickness_units) + "Layer Thickness", thickness_units, conversion=GV%H_to_MKS) - call register_restart_pair(CS%u, CS%v, u_desc, v_desc, .true., restart_CSp) + u_desc = var_desc("u", "m s-1", "Zonal velocity", hor_grid='Cu') + v_desc = var_desc("v", "m s-1", "Meridional velocity", hor_grid='Cv') + call register_restart_pair(CS%u, CS%v, u_desc, v_desc, .true., restart_CSp, conversion=US%L_T_to_m_s) if (associated(CS%tv%frazil)) & call register_restart_field(CS%tv%frazil, "frazil", .false., restart_CSp, & - "Frazil heat flux into ocean", "J m-2") + "Frazil heat flux into ocean", & + "J m-2", conversion=US%Q_to_J_kg*US%RZ_to_kg_m2) if (CS%interp_p_surf) then call register_restart_field(CS%p_surf_prev, "p_surf_prev", .false., restart_CSp, & - "Previous ocean surface pressure", "Pa") + "Previous ocean surface pressure", "Pa", conversion=US%RL2_T2_to_Pa) endif call register_restart_field(CS%ave_ssh_ibc, "ave_ssh", .false., restart_CSp, & - "Time average sea surface height", "meter") + "Time average sea surface height", "meter", conversion=US%Z_to_m) ! hML is needed when using the ice shelf module call get_param(param_file, '', "ICE_SHELF", use_ice_shelf, default=.false., & do_not_log=.true.) if (use_ice_shelf .and. associated(CS%Hml)) then call register_restart_field(CS%Hml, "hML", .false., restart_CSp, & - "Mixed layer thickness", "meter") + "Mixed layer thickness", "meter", conversion=US%Z_to_m) endif ! Register scalar unit conversion factors. diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 3cb1ebf399..3b5c812ba1 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -4748,8 +4748,8 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, dtbt_tmp = -1.0 if (query_initialized(CS%dtbt, "DTBT", restart_CS)) then dtbt_tmp = CS%dtbt - if ((US%s_to_T_restart /= 0.0) .and. (US%s_to_T_restart /= US%s_to_T)) & - dtbt_tmp = (US%s_to_T / US%s_to_T_restart) * CS%dtbt + if ((US%s_to_T_restart /= 0.0) .and. (US%s_to_T_restart /= 1.0)) & + dtbt_tmp = (1.0 / US%s_to_T_restart) * CS%dtbt endif ! Estimate the maximum stable barotropic time step. @@ -4909,8 +4909,8 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, CS%vbtav(i,J) = CS%vbtav(i,J) + CS%frhatv(i,J,k) * v(i,J,k) enddo ; enddo ; enddo elseif ((US%s_to_T_restart*US%m_to_L_restart /= 0.0) .and. & - (US%m_to_L*US%s_to_T_restart) /= (US%m_to_L_restart*US%s_to_T)) then - vel_rescale = (US%m_to_L*US%s_to_T_restart) / (US%m_to_L_restart*US%s_to_T) + (US%s_to_T_restart /= US%m_to_L_restart)) then + vel_rescale = US%s_to_T_restart / US%m_to_L_restart do j=js,je ; do I=is-1,ie ; CS%ubtav(I,j) = vel_rescale * CS%ubtav(I,j) ; enddo ; enddo do J=js-1,je ; do i=is,ie ; CS%vbtav(i,J) = vel_rescale * CS%vbtav(i,J) ; enddo ; enddo endif @@ -4921,8 +4921,8 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, do j=js,je ; do I=is-1,ie ; CS%ubt_IC(I,j) = CS%ubtav(I,j) ; enddo ; enddo do J=js-1,je ; do i=is,ie ; CS%vbt_IC(i,J) = CS%vbtav(i,J) ; enddo ; enddo elseif ((US%s_to_T_restart*US%m_to_L_restart /= 0.0) .and. & - (US%m_to_L*US%s_to_T_restart) /= (US%m_to_L_restart*US%s_to_T)) then - vel_rescale = (US%m_to_L*US%s_to_T_restart) / (US%m_to_L_restart*US%s_to_T) + (US%s_to_T_restart /= US%m_to_L_restart)) then + vel_rescale = US%s_to_T_restart / US%m_to_L_restart do j=js,je ; do I=is-1,ie ; CS%ubt_IC(I,j) = vel_rescale * CS%ubt_IC(I,j) ; enddo ; enddo do J=js-1,je ; do i=is,ie ; CS%vbt_IC(i,J) = vel_rescale * CS%vbt_IC(i,J) ; enddo ; enddo endif @@ -5022,11 +5022,12 @@ end subroutine barotropic_end !> This subroutine is used to register any fields from MOM_barotropic.F90 !! that should be written to or read from the restart file. -subroutine register_barotropic_restarts(HI, GV, param_file, CS, restart_CS) +subroutine register_barotropic_restarts(HI, GV, US, param_file, CS, restart_CS) type(hor_index_type), intent(in) :: HI !< A horizontal index type 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 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters. type(barotropic_CS), intent(inout) :: CS !< Barotropic control structure - type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(MOM_restart_CS), intent(inout) :: restart_CS !< MOM restart control structure ! Local variables @@ -5056,7 +5057,8 @@ subroutine register_barotropic_restarts(HI, GV, param_file, CS, restart_CS) hor_grid='u', z_grid='1') vd(3) = var_desc("vbtav","m s-1","Time mean barotropic meridional velocity",& hor_grid='v', z_grid='1') - call register_restart_pair(CS%ubtav, CS%vbtav, vd(2), vd(3), .false., restart_CS) + call register_restart_pair(CS%ubtav, CS%vbtav, vd(2), vd(3), .false., restart_CS, & + conversion=US%L_T_to_m_s) if (CS%gradual_BT_ICs) then vd(2) = var_desc("ubt_IC", "m s-1", & @@ -5065,12 +5067,12 @@ subroutine register_barotropic_restarts(HI, GV, param_file, CS, restart_CS) vd(3) = var_desc("vbt_IC", "m s-1", & longname="Next initial condition for the barotropic meridional velocity",& hor_grid='v', z_grid='1') - call register_restart_pair(CS%ubt_IC, CS%vbt_IC, vd(2), vd(3), .false., restart_CS) + call register_restart_pair(CS%ubt_IC, CS%vbt_IC, vd(2), vd(3), .false., restart_CS, & + conversion=US%L_T_to_m_s) endif - call register_restart_field(CS%dtbt, "DTBT", .false., restart_CS, & - longname="Barotropic timestep", units="seconds") + longname="Barotropic timestep", units="seconds", conversion=US%T_to_s) end subroutine register_barotropic_restarts diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index f22fb9a862..d177e63be8 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -972,9 +972,10 @@ end subroutine step_MOM_dyn_split_RK2 !> This subroutine sets up any auxiliary restart variables that are specific !! to the split-explicit time stepping scheme. All variables registered here should !! have the ability to be recreated if they are not present in a restart file. -subroutine register_restarts_dyn_split_RK2(HI, GV, param_file, CS, restart_CS, uh, vh) +subroutine register_restarts_dyn_split_RK2(HI, GV, US, param_file, CS, restart_CS, uh, vh) type(hor_index_type), intent(in) :: HI !< Horizontal index structure type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(param_file_type), intent(in) :: param_file !< parameter file type(MOM_dyn_split_RK2_CS), pointer :: CS !< module control structure type(MOM_restart_CS), intent(inout) :: restart_CS !< MOM restart control structure @@ -1016,29 +1017,32 @@ subroutine register_restarts_dyn_split_RK2(HI, GV, param_file, CS, restart_CS, u flux_units = get_flux_units(GV) if (GV%Boussinesq) then - vd(1) = var_desc("sfc",thickness_units,"Free surface Height",'h','1') + call register_restart_field(CS%eta, "sfc", .false., restart_CS, & + longname="Free surface Height", units=thickness_units, conversion=GV%H_to_mks) else - vd(1) = var_desc("p_bot",thickness_units,"Bottom Pressure",'h','1') + call register_restart_field(CS%eta, "p_bot", .false., restart_CS, & + longname="Bottom Pressure", units=thickness_units, conversion=GV%H_to_mks) endif - call register_restart_field(CS%eta, vd(1), .false., restart_CS) - vd(1) = var_desc("u2","m s-1","Auxiliary Zonal velocity",'u','L') - vd(2) = var_desc("v2","m s-1","Auxiliary Meridional velocity",'v','L') - call register_restart_pair(CS%u_av, CS%v_av, vd(1), vd(2), .false., restart_CS) + vd(1) = var_desc("u2", "m s-1", "Auxiliary Zonal velocity", 'u', 'L') + vd(2) = var_desc("v2", "m s-1", "Auxiliary Meridional velocity", 'v', 'L') + call register_restart_pair(CS%u_av, CS%v_av, vd(1), vd(2), .false., restart_CS, & + conversion=US%L_T_to_m_s) - vd(1) = var_desc("h2",thickness_units,"Auxiliary Layer Thickness",'h','L') - call register_restart_field(CS%h_av, vd(1), .false., restart_CS) + call register_restart_field(CS%h_av, "h2", .false., restart_CS, & + longname="Auxiliary Layer Thickness", units=thickness_units, conversion=GV%H_to_mks) - vd(1) = var_desc("uh",flux_units,"Zonal thickness flux",'u','L') - vd(2) = var_desc("vh",flux_units,"Meridional thickness flux",'v','L') - call register_restart_pair(uh, vh, vd(1), vd(2), .false., restart_CS) + vd(1) = var_desc("uh", flux_units, "Zonal thickness flux", 'u', 'L') + vd(2) = var_desc("vh", flux_units, "Meridional thickness flux", 'v', 'L') + call register_restart_pair(uh, vh, vd(1), vd(2), .false., restart_CS, & + conversion=GV%H_to_MKS*US%L_to_m**2*US%s_to_T) - vd(1) = var_desc("diffu","m s-2","Zonal horizontal viscous acceleration",'u','L') - vd(2) = var_desc("diffv","m s-2","Meridional horizontal viscous acceleration",'v','L') - call register_restart_pair(CS%diffu, CS%diffv, vd(1), vd(2), .false., restart_CS) + vd(1) = var_desc("diffu", "m s-2", "Zonal horizontal viscous acceleration", 'u', 'L') + vd(2) = var_desc("diffv", "m s-2", "Meridional horizontal viscous acceleration", 'v', 'L') + call register_restart_pair(CS%diffu, CS%diffv, vd(1), vd(2), .false., restart_CS, & + conversion=US%L_T2_to_m_s2) - call register_barotropic_restarts(HI, GV, param_file, CS%barotropic_CSp, & - restart_CS) + call register_barotropic_restarts(HI, GV, US, param_file, CS%barotropic_CSp, restart_CS) end subroutine register_restarts_dyn_split_RK2 @@ -1238,8 +1242,8 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param do k=1,nz ; do j=js,je ; do i=is,ie CS%eta(i,j) = CS%eta(i,j) + h(i,j,k) enddo ; enddo ; enddo - elseif ((GV%m_to_H_restart /= 0.0) .and. (GV%m_to_H_restart /= GV%m_to_H)) then - H_rescale = GV%m_to_H / GV%m_to_H_restart + elseif ((GV%m_to_H_restart /= 0.0) .and. (GV%m_to_H_restart /= 1.0)) then + H_rescale = 1.0 / GV%m_to_H_restart do j=js,je ; do i=is,ie ; CS%eta(i,j) = H_rescale * CS%eta(i,j) ; enddo ; enddo endif ! Copy eta into an output array. @@ -1251,14 +1255,12 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param if (.not. query_initialized(CS%diffu,"diffu",restart_CS) .or. & .not. query_initialized(CS%diffv,"diffv",restart_CS)) then - call horizontal_viscosity(u, v, h, CS%diffu, CS%diffv, MEKE, VarMix, & - G, GV, US, CS%hor_visc, & - OBC=CS%OBC, BT=CS%barotropic_CSp, & - TD=thickness_diffuse_CSp) + call horizontal_viscosity(u, v, h, CS%diffu, CS%diffv, MEKE, VarMix, G, GV, US, CS%hor_visc, & + OBC=CS%OBC, BT=CS%barotropic_CSp, TD=thickness_diffuse_CSp) else if ( (US%s_to_T_restart * US%m_to_L_restart /= 0.0) .and. & - (US%m_to_L * US%s_to_T_restart**2 /= US%m_to_L_restart * US%s_to_T**2) ) then - accel_rescale = (US%m_to_L * US%s_to_T_restart**2) / (US%m_to_L_restart * US%s_to_T**2) + (US%s_to_T_restart**2 /= US%m_to_L_restart) ) then + accel_rescale = US%s_to_T_restart**2 / US%m_to_L_restart do k=1,nz ; do j=js,je ; do I=G%IscB,G%IecB CS%diffu(I,j,k) = accel_rescale * CS%diffu(I,j,k) enddo ; enddo ; enddo @@ -1273,8 +1275,8 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param do k=1,nz ; do j=jsd,jed ; do I=IsdB,IedB ; CS%u_av(I,j,k) = u(I,j,k) ; enddo ; enddo ; enddo do k=1,nz ; do J=JsdB,JedB ; do i=isd,ied ; CS%v_av(i,J,k) = v(i,J,k) ; enddo ; enddo ; enddo elseif ( (US%s_to_T_restart * US%m_to_L_restart /= 0.0) .and. & - ((US%m_to_L * US%s_to_T_restart) /= (US%m_to_L_restart * US%s_to_T)) ) then - vel_rescale = (US%m_to_L * US%s_to_T_restart) / (US%m_to_L_restart * US%s_to_T) + (US%s_to_T_restart /= US%m_to_L_restart) ) then + vel_rescale = US%s_to_T_restart / US%m_to_L_restart do k=1,nz ; do j=jsd,jed ; do I=IsdB,IedB ; CS%u_av(I,j,k) = vel_rescale * CS%u_av(I,j,k) ; enddo ; enddo ; enddo do k=1,nz ; do J=JsdB,JedB ; do i=isd,ied ; CS%v_av(i,J,k) = vel_rescale * CS%v_av(i,J,k) ; enddo ; enddo ; enddo endif @@ -1291,15 +1293,13 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param else if (.not. query_initialized(CS%h_av,"h2",restart_CS)) then CS%h_av(:,:,:) = h(:,:,:) - elseif ((GV%m_to_H_restart /= 0.0) .and. (GV%m_to_H_restart /= GV%m_to_H)) then - H_rescale = GV%m_to_H / GV%m_to_H_restart + elseif ((GV%m_to_H_restart /= 0.0) .and. (GV%m_to_H_restart /= 1.0)) then + H_rescale = 1.0 / GV%m_to_H_restart do k=1,nz ; do j=js,je ; do i=is,ie ; CS%h_av(i,j,k) = H_rescale * CS%h_av(i,j,k) ; enddo ; enddo ; enddo endif if ( (GV%m_to_H_restart * US%s_to_T_restart * US%m_to_L_restart /= 0.0) .and. & - ((GV%m_to_H * US%m_to_L**2 * US%s_to_T_restart) /= & - (GV%m_to_H_restart * US%m_to_L_restart**2 * US%s_to_T)) ) then - uH_rescale = (GV%m_to_H * US%m_to_L**2 * US%s_to_T_restart) / & - (GV%m_to_H_restart * US%m_to_L_restart**2 * US%s_to_T) + (US%s_to_T_restart /= (GV%m_to_H_restart * US%m_to_L_restart**2)) ) then + uH_rescale = US%s_to_T_restart / (GV%m_to_H_restart * US%m_to_L_restart**2) do k=1,nz ; do j=js,je ; do I=G%IscB,G%IecB ; uh(I,j,k) = uH_rescale * uh(I,j,k) ; enddo ; enddo ; enddo do k=1,nz ; do J=G%JscB,G%JecB ; do i=is,ie ; vh(i,J,k) = uH_rescale * vh(i,J,k) ; enddo ; enddo ; enddo endif diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 5fc168f5a4..6ce0940ddb 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -1821,8 +1821,8 @@ subroutine open_boundary_init(G, GV, US, param_file, OBC, restart_CS) ! points per timestep, but if this were to be corrected to [L T-1 ~> m s-1] or [T-1 ~> s-1] to ! permit timesteps to change between calls to the OBC code, the following would be needed: ! if ( OBC%radiation_BCs_exist_globally .and. (US%s_to_T_restart * US%m_to_L_restart /= 0.0) .and. & -! ((US%m_to_L * US%s_to_T_restart) /= (US%m_to_L_restart * US%s_to_T)) ) then -! vel_rescale = (US%m_to_L * US%s_to_T_restart) / (US%m_to_L_restart * US%s_to_T) +! (US%s_to_T_restart /= US%m_to_L_restart) ) then +! vel_rescale = US%s_to_T_restart / US%m_to_L_restart ! if (query_initialized(OBC%rx_normal, "rx_normal", restart_CS)) then ! do k=1,nz ; do j=jsd,jed ; do I=IsdB,IedB ! OBC%rx_normal(I,j,k) = vel_rescale * OBC%rx_normal(I,j,k) @@ -1837,8 +1837,8 @@ subroutine open_boundary_init(G, GV, US, param_file, OBC, restart_CS) ! The oblique boundary condition terms have units of [L2 T-2 ~> m2 s-2] and may need to be rescaled. if ( OBC%oblique_BCs_exist_globally .and. (US%s_to_T_restart * US%m_to_L_restart /= 0.0) .and. & - ((US%m_to_L * US%s_to_T_restart) /= (US%m_to_L_restart * US%s_to_T)) ) then - vel2_rescale = (US%m_to_L * US%s_to_T_restart)**2 / (US%m_to_L_restart * US%s_to_T)**2 + (US%s_to_T_restart /= US%m_to_L_restart) ) then + vel2_rescale = US%s_to_T_restart**2 / US%m_to_L_restart**2 if (query_initialized(OBC%rx_oblique, "rx_oblique", restart_CS)) then do k=1,nz ; do j=jsd,jed ; do I=IsdB,IedB OBC%rx_oblique(I,j,k) = vel2_rescale * OBC%rx_oblique(I,j,k) @@ -4910,10 +4910,11 @@ subroutine flood_fill2(G, color, cin, cout, cland) end subroutine flood_fill2 !> Register OBC segment data for restarts -subroutine open_boundary_register_restarts(HI, GV, OBC, Reg, param_file, restart_CS, & +subroutine open_boundary_register_restarts(HI, GV, US, OBC, Reg, param_file, restart_CS, & use_temperature) type(hor_index_type), intent(in) :: HI !< Horizontal indices type(verticalGrid_type), pointer :: GV !< Container for vertical grid information + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(ocean_OBC_type), pointer :: OBC !< OBC data structure, data intent(inout) type(tracer_registry_type), pointer :: Reg !< pointer to tracer registry type(param_file_type), intent(in) :: param_file !< Parameter file handle @@ -4922,25 +4923,31 @@ subroutine open_boundary_register_restarts(HI, GV, OBC, Reg, param_file, restart ! Local variables type(vardesc) :: vd(2) integer :: m, n - character(len=100) :: mesg + character(len=100) :: mesg, var_name type(OBC_segment_type), pointer :: segment=>NULL() if (.not. associated(OBC)) & call MOM_error(FATAL, "open_boundary_register_restarts: Called with "//& "uninitialized OBC control structure") - ! *** This is a temporary work around for restarts with OBC segments. + ! ### This is a temporary work around for restarts with OBC segments. ! This implementation uses 3D arrays solely for restarts. We need ! to be able to add 2D ( x,z or y,z ) data to restarts to avoid using - ! so much memory and disk space. *** + ! so much memory and disk space. if (OBC%radiation_BCs_exist_globally) then allocate(OBC%rx_normal(HI%isdB:HI%iedB,HI%jsd:HI%jed,GV%ke), source=0.0) allocate(OBC%ry_normal(HI%isd:HI%ied,HI%jsdB:HI%jedB,GV%ke), source=0.0) - vd(1) = var_desc("rx_normal", "m s-1", "Normal Phase Speed for EW radiation OBCs", 'u', 'L') - vd(2) = var_desc("ry_normal", "m s-1", "Normal Phase Speed for NS radiation OBCs", 'v', 'L') - call register_restart_pair(OBC%rx_normal, OBC%ry_normal, vd(1), vd(2), & - .false., restart_CS) + vd(1) = var_desc("rx_normal", "gridpoint timestep-1", "Normal Phase Speed for EW radiation OBCs", 'u', 'L') + vd(2) = var_desc("ry_normal", "gridpoint timestep-1", "Normal Phase Speed for NS radiation OBCs", 'v', 'L') + call register_restart_pair(OBC%rx_normal, OBC%ry_normal, vd(1), vd(2), .false., restart_CS) + ! The rx_normal and ry_normal arrays used with radiation OBCs are currently in units of grid + ! points per timestep, but if this were to be corrected to [L T-1 ~> m s-1] or [T-1 ~> s-1] to + ! permit timesteps to change between calls to the OBC code, the following would be needed instead: + ! vd(1) = var_desc("rx_normal", "m s-1", "Normal Phase Speed for EW radiation OBCs", 'u', 'L') + ! vd(2) = var_desc("ry_normal", "m s-1", "Normal Phase Speed for NS radiation OBCs", 'v', 'L') + ! call register_restart_pair(OBC%rx_normal, OBC%ry_normal, vd(1), vd(2), .false., restart_CS, & + ! conversion=US%L_T_to_m_s) endif if (OBC%oblique_BCs_exist_globally) then @@ -4949,12 +4956,13 @@ subroutine open_boundary_register_restarts(HI, GV, OBC, Reg, param_file, restart vd(1) = var_desc("rx_oblique", "m2 s-2", "Radiation Speed Squared for EW oblique OBCs", 'u', 'L') vd(2) = var_desc("ry_oblique", "m2 s-2", "Radiation Speed Squared for NS oblique OBCs", 'v', 'L') - call register_restart_pair(OBC%rx_oblique, OBC%ry_oblique, vd(1), vd(2), & - .false., restart_CS) + call register_restart_pair(OBC%rx_oblique, OBC%ry_oblique, vd(1), vd(2), .false., & + restart_CS, conversion=US%L_T_to_m_s**2) allocate(OBC%cff_normal(HI%IsdB:HI%IedB,HI%jsdB:HI%jedB,GV%ke), source=0.0) - vd(1) = var_desc("cff_normal", "m2 s-2", "denominator for oblique OBCs", 'q', 'L') - call register_restart_field(OBC%cff_normal, vd(1), .false., restart_CS) + call register_restart_field(OBC%cff_normal, "cff_normal", .false., restart_CS, & + longname="denominator for oblique OBCs", & + units="m2 s-2", conversion=US%L_T_to_m_s**2, hor_grid="q") endif if (Reg%ntr == 0) return @@ -4978,13 +4986,13 @@ subroutine open_boundary_register_restarts(HI, GV, OBC, Reg, param_file, restart do m=1,OBC%ntr if (OBC%tracer_x_reservoirs_used(m)) then if (modulo(HI%turns, 2) /= 0) then - write(mesg,'("tres_y_",I3.3)') m - vd(1) = var_desc(mesg,"Conc", "Tracer concentration for NS OBCs",'v','L') - call register_restart_field(OBC%tres_x(:,:,:,m), vd(1), .false., restart_CS) + write(var_name,'("tres_y_",I3.3)') m + call register_restart_field(OBC%tres_x(:,:,:,m), var_name, .false., restart_CS, & + longname="Tracer concentration for NS OBCs", units="Conc", hor_grid='v') else - write(mesg,'("tres_x_",I3.3)') m - vd(1) = var_desc(mesg,"Conc", "Tracer concentration for EW OBCs",'u','L') - call register_restart_field(OBC%tres_x(:,:,:,m), vd(1), .false., restart_CS) + write(var_name,'("tres_x_",I3.3)') m + call register_restart_field(OBC%tres_x(:,:,:,m), var_name, .false., restart_CS, & + longname="Tracer concentration for EW OBCs", units="Conc", hor_grid='u') endif endif enddo @@ -4994,13 +5002,13 @@ subroutine open_boundary_register_restarts(HI, GV, OBC, Reg, param_file, restart do m=1,OBC%ntr if (OBC%tracer_y_reservoirs_used(m)) then if (modulo(HI%turns, 2) /= 0) then - write(mesg,'("tres_x_",I3.3)') m - vd(1) = var_desc(mesg,"Conc", "Tracer concentration for EW OBCs",'u','L') - call register_restart_field(OBC%tres_y(:,:,:,m), vd(1), .false., restart_CS) + write(var_name,'("tres_x_",I3.3)') m + call register_restart_field(OBC%tres_y(:,:,:,m), var_name, .false., restart_CS, & + longname="Tracer concentration for EW OBCs", units="Conc", hor_grid='u') else - write(mesg,'("tres_y_",I3.3)') m - vd(1) = var_desc(mesg,"Conc", "Tracer concentration for NS OBCs",'v','L') - call register_restart_field(OBC%tres_y(:,:,:,m), vd(1), .false., restart_CS) + write(var_name,'("tres_y_",I3.3)') m + call register_restart_field(OBC%tres_y(:,:,:,m), var_name, .false., restart_CS, & + longname="Tracer concentration for NS OBCs", units="Conc", hor_grid='v') endif endif enddo diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 13af5a936a..1e818bd298 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -1634,11 +1634,11 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, call restart_init(param_file, CS%restart_CSp, "Shelf.res") call register_restart_field(ISS%mass_shelf, "shelf_mass", .true., CS%restart_CSp, & - "Ice shelf mass", "kg m-2") + "Ice shelf mass", "kg m-2", conversion=US%RZ_to_kg_m2) call register_restart_field(ISS%area_shelf_h, "shelf_area", .true., CS%restart_CSp, & - "Ice shelf area in cell", "m2") + "Ice shelf area in cell", "m2", conversion=US%L_to_m**2) call register_restart_field(ISS%h_shelf, "h_shelf", .true., CS%restart_CSp, & - "ice sheet/shelf thickness", "m") + "ice sheet/shelf thickness", "m", conversion=US%Z_to_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", & @@ -1646,7 +1646,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, 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) + .false., CS%restart_CSp, conversion=US%RZ_T_to_kg_m2s*US%L_T_to_m_s) endif endif @@ -1663,13 +1663,13 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, if (CS%active_shelf_dynamics) then ! Allocate CS%dCS and specify additional restarts for ice shelf dynamics - call register_ice_shelf_dyn_restarts(CS%Grid_in, param_file, CS%dCS, CS%restart_CSp) + call register_ice_shelf_dyn_restarts(CS%Grid_in, US, 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 !if (.not. CS%solo_ice_sheet) then ! call register_restart_field(fluxes%ustar_shelf, "ustar_shelf", .false., CS%restart_CSp, & - ! "Friction velocity under ice shelves", "m s-1") + ! "Friction velocity under ice shelves", "m s-1", conversion=###) !endif CS%restart_output_dir = dirs%restart_output_dir @@ -1698,23 +1698,23 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, call MOM_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: Restoring ice shelf from file.") call restore_state(dirs%input_filename, dirs%restart_input_dir, Time, G, CS%restart_CSp) - if ((US%m_to_Z_restart /= 0.0) .and. (US%m_to_Z_restart /= US%m_to_Z)) then - Z_rescale = US%m_to_Z / US%m_to_Z_restart + if ((US%m_to_Z_restart /= 0.0) .and. (US%m_to_Z_restart /= 1.0)) then + Z_rescale = 1.0 / US%m_to_Z_restart do j=G%jsc,G%jec ; do i=G%isc,G%iec ISS%h_shelf(i,j) = Z_rescale * ISS%h_shelf(i,j) enddo ; enddo endif if ((US%m_to_Z_restart*US%kg_m3_to_R_restart /= 0.0) .and. & - (US%m_to_Z*US%kg_m3_to_R /= US%m_to_Z_restart * US%kg_m3_to_R_restart)) then - RZ_rescale = US%m_to_Z*US%kg_m3_to_R / (US%m_to_Z_restart * US%kg_m3_to_R_restart) + (US%m_to_Z_restart*US%kg_m3_to_R_restart /= 1.0)) then + RZ_rescale = 1.0 / (US%m_to_Z_restart * US%kg_m3_to_R_restart) do j=G%jsc,G%jec ; do i=G%isc,G%iec ISS%mass_shelf(i,j) = RZ_rescale * ISS%mass_shelf(i,j) enddo ; enddo endif - if ((US%m_to_L_restart /= 0.0) .and. (US%m_to_L_restart /= US%m_to_L)) then - L_rescale = US%m_to_L / US%m_to_L_restart + if ((US%m_to_L_restart /= 0.0) .and. (US%m_to_L_restart /= 1.0)) then + L_rescale = 1.0 / US%m_to_L_restart do j=G%jsc,G%jec ; do i=G%isc,G%iec ISS%area_shelf_h(i,j) = L_rescale**2 * ISS%area_shelf_h(i,j) enddo ; enddo diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index 8fb674e36c..606ab49d8c 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -224,8 +224,9 @@ end function quad_area !> This subroutine is used to register any fields related to the ice shelf !! dynamics that should be written to or read from the restart file. -subroutine register_ice_shelf_dyn_restarts(G, param_file, CS, restart_CS) +subroutine register_ice_shelf_dyn_restarts(G, US, param_file, CS, restart_CS) type(ocean_grid_type), intent(inout) :: G !< The grid type describing the ice shelf grid. + type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters type(ice_shelf_dyn_CS), pointer :: CS !< A pointer to the ice shelf dynamics control structure type(MOM_restart_CS), intent(inout) :: restart_CS !< MOM restart control struct @@ -275,20 +276,24 @@ subroutine register_ice_shelf_dyn_restarts(G, param_file, CS, restart_CS) allocate( CS%h_bdry_val(isd:ied,jsd:jed), source=0.0 ) ! additional restarts for ice shelf state call register_restart_field(CS%u_shelf, "u_shelf", .false., restart_CS, & - "ice sheet/shelf u-velocity", "m s-1", hor_grid='Bu') + "ice sheet/shelf u-velocity", & + units="m s-1", conversion=US%L_T_to_m_s, hor_grid='Bu') call register_restart_field(CS%v_shelf, "v_shelf", .false., restart_CS, & - "ice sheet/shelf v-velocity", "m s-1", hor_grid='Bu') + "ice sheet/shelf v-velocity", & + units="m s-1", conversion=US%L_T_to_m_s, hor_grid='Bu') call register_restart_field(CS%u_bdry_val, "u_bdry_val", .false., restart_CS, & - "ice sheet/shelf boundary u-velocity", "m s-1", hor_grid='Bu') + "ice sheet/shelf boundary u-velocity", & + units="m s-1", conversion=US%L_T_to_m_s, hor_grid='Bu') call register_restart_field(CS%v_bdry_val, "v_bdry_val", .false., restart_CS, & - "ice sheet/shelf boundary v-velocity", "m s-1", hor_grid='Bu') + "ice sheet/shelf boundary v-velocity", & + units="m s-1", conversion=US%L_T_to_m_s, hor_grid='Bu') call register_restart_field(CS%u_face_mask_bdry, "u_face_mask_bdry", .false., restart_CS, & "ice sheet/shelf boundary u-mask", "nondim", hor_grid='Bu') call register_restart_field(CS%v_face_mask_bdry, "v_face_mask_bdry", .false., restart_CS, & "ice sheet/shelf boundary v-mask", "nondim", hor_grid='Bu') call register_restart_field(CS%OD_av, "OD_av", .true., restart_CS, & - "Average open ocean depth in a cell","m") + "Average open ocean depth in a cell", "m", conversion=US%Z_to_m) call register_restart_field(CS%ground_frac, "ground_frac", .true., restart_CS, & "fractional degree of grounding", "nondim") call register_restart_field(CS%C_basal_friction, "C_basal_friction", .true., restart_CS, & @@ -296,7 +301,7 @@ subroutine register_ice_shelf_dyn_restarts(G, param_file, CS, restart_CS) call register_restart_field(CS%AGlen_visc, "AGlen_visc", .true., restart_CS, & "ice-stiffness parameter", "Pa-3 s-1") call register_restart_field(CS%h_bdry_val, "h_bdry_val", .false., restart_CS, & - "ice thickness at the boundary","m") + "ice thickness at the boundary", "m", conversion=US%Z_to_m) endif end subroutine register_ice_shelf_dyn_restarts @@ -462,16 +467,16 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ ! Take additional initialization steps, for example of dependent variables. if (active_shelf_dynamics .and. .not.new_sim) then - if ((US%m_to_Z_restart /= 0.0) .and. (US%m_to_Z_restart /= US%m_to_Z)) then - Z_rescale = US%m_to_Z / US%m_to_Z_restart + if ((US%m_to_Z_restart /= 0.0) .and. (US%m_to_Z_restart /= 1.0)) then + Z_rescale = 1.0 / US%m_to_Z_restart do j=G%jsc,G%jec ; do i=G%isc,G%iec CS%OD_av(i,j) = Z_rescale * CS%OD_av(i,j) enddo ; enddo endif if ((US%m_to_L_restart*US%s_to_T_restart /= 0.0) .and. & - (US%m_to_L_restart /= US%m_s_to_L_T*US%s_to_T_restart)) then - vel_rescale = US%m_s_to_L_T*US%s_to_T_restart / US%m_to_L_restart + (US%m_to_L_restart /= US%s_to_T_restart)) then + vel_rescale = US%s_to_T_restart / US%m_to_L_restart do J=G%jsc-1,G%jec ; do I=G%isc-1,G%iec CS%u_shelf(I,J) = vel_rescale * CS%u_shelf(I,J) CS%v_shelf(I,J) = vel_rescale * CS%v_shelf(I,J) diff --git a/src/ice_shelf/MOM_ice_shelf_initialize.F90 b/src/ice_shelf/MOM_ice_shelf_initialize.F90 index 7cc3c020a3..4e96bab8f5 100644 --- a/src/ice_shelf/MOM_ice_shelf_initialize.F90 +++ b/src/ice_shelf/MOM_ice_shelf_initialize.F90 @@ -309,7 +309,7 @@ subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_b intent(inout) :: hmask !< A mask indicating which tracer points are !! partly or fully covered by an ice-shelf real, dimension(SZDI_(G),SZDJ_(G)), & - intent(inout) :: h_shelf !< Ice-shelf thickness + intent(inout) :: h_shelf !< Ice-shelf thickness [Z ~> m] 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 @@ -330,7 +330,7 @@ subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_b call get_param(PF, mdl, "INPUT_VEL_ICE_SHELF", input_vel, & "inflow ice velocity at upstream boundary", & - units="m s-1", default=0., scale=US%m_s_to_L_T*US%m_to_Z) + units="m s-1", default=0., scale=US%m_s_to_L_T*US%m_to_Z) !### This conversion factor is wrong? call get_param(PF, mdl, "INPUT_THICK_ICE_SHELF", input_thick, & "flux thickness at upstream boundary", & units="m", default=1000., scale=US%m_to_Z) @@ -409,7 +409,7 @@ subroutine initialize_ice_flow_from_file(bed_elev,u_shelf, v_shelf,float_cond,& real, dimension(SZDI_(G),SZDJ_(G)), & intent(inout) :: float_cond !< An array indicating where the ice - !! shelf is floating: 0 if floating, 1 if not. + !! shelf is floating: 0 if floating, 1 if not. [nondim] ! real, dimension(SZDI_(G),SZDJ_(G)), & ! intent(in) :: hmask !< A mask indicating which tracer points are !! partly or fully covered by an ice-shelf @@ -461,13 +461,14 @@ subroutine initialize_ice_flow_from_file(bed_elev,u_shelf, v_shelf,float_cond,& floatfr_varname = "float_frac" - call MOM_read_data(filename, trim(ushelf_varname), u_shelf, G%Domain, position=CORNER,scale=1.0) - call MOM_read_data(filename,trim(vshelf_varname), v_shelf, G%Domain, position=CORNER,scale=1.0) -! call MOM_read_data(filename,trim(ice_visc_varname), ice_visc, G%Domain,position=CORNER,scale=1.0) - call MOM_read_data(filename,trim(floatfr_varname), float_cond, G%Domain, scale=1.) + !### I think that the following two lines should have ..., scale=US%m_s_to_L_T + call MOM_read_data(filename, trim(ushelf_varname), u_shelf, G%Domain, position=CORNER, scale=1.0) + call MOM_read_data(filename, trim(vshelf_varname), v_shelf, G%Domain, position=CORNER, scale=1.0) +! call MOM_read_data(filename, trim(ice_visc_varname), ice_visc, G%Domain,position=CORNER,scale=1.0) + call MOM_read_data(filename, trim(floatfr_varname), float_cond, G%Domain, scale=1.) - filename = trim(inputdir)//trim(bed_topo_file) - call MOM_read_data(filename,trim(bed_varname), bed_elev, G%Domain, scale=1.) + filename = trim(inputdir)//trim(bed_topo_file) + call MOM_read_data(filename,trim(bed_varname), bed_elev, G%Domain, scale=1.) ! isc = G%isc ; jsc = G%jsc ; iec = G%iec ; jec = G%jec @@ -480,18 +481,19 @@ subroutine initialize_ice_shelf_boundary_from_file(u_face_mask_bdry, v_face_mask type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure real, dimension(SZIB_(G),SZJB_(G)), & - intent(inout) :: u_face_mask_bdry !< A boundary-type mask at B-grid u faces + intent(inout) :: u_face_mask_bdry !< A boundary-type mask at B-grid u faces [nondim] real, dimension(SZIB_(G),SZJB_(G)), & - intent(inout) :: v_face_mask_bdry !< A boundary-type mask at B-grid v faces + intent(inout) :: v_face_mask_bdry !< A boundary-type mask at B-grid v faces [nondim] real, dimension(SZIB_(G),SZJB_(G)), & intent(inout) :: u_bdry_val !< The zonal ice shelf velocity at open !! boundary vertices [L T-1 ~> m s-1]. real, dimension(SZIB_(G),SZJB_(G)), & intent(inout) :: v_bdry_val !< The meridional ice shelf velocity at open + !! boundary vertices [L T-1 ~> m s-1]. real, dimension(SZDIB_(G),SZDJB_(G)), & - intent(inout) :: umask !< A mask foor ice shelf velocity + intent(inout) :: umask !< A mask for ice shelf velocity [nondim] real, dimension(SZDIB_(G),SZDJB_(G)), & - intent(inout) :: vmask !< A mask foor ice shelf velocity + intent(inout) :: vmask !< A mask for ice shelf velocity [nondim] real, dimension(SZDI_(G),SZDJ_(G)), & intent(inout) :: thickness_bdry_val !< The ice shelf thickness at open boundaries [Z ~> m] !! boundary vertices [L T-1 ~> m s-1]. @@ -499,9 +501,9 @@ subroutine initialize_ice_shelf_boundary_from_file(u_face_mask_bdry, v_face_mask intent(inout) :: h_bdry_val !< The ice shelf thickness at open boundaries [Z ~> m] real, dimension(SZDI_(G),SZDJ_(G)), & intent(inout) :: hmask !< A mask indicating which tracer points are - !! partly or fully covered by an ice-shelf + !! partly or fully covered by an ice-shelf [nondim] real, dimension(SZDI_(G),SZDJ_(G)), & - intent(in) :: h_shelf !< Ice-shelf thickness + intent(in) :: h_shelf !< Ice-shelf thickness [Z ~> m] 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 @@ -557,16 +559,17 @@ subroutine initialize_ice_shelf_boundary_from_file(u_face_mask_bdry, v_face_mask " initialize_ice_shelf_velocity_from_file: Unable to open "//trim(filename)) - call MOM_read_data(filename, trim(ufcmskbdry_varname),u_face_mask_bdry, G%Domain,position=CORNER,scale=1.0) - call MOM_read_data(filename,trim(vfcmskbdry_varname), v_face_mask_bdry, G%Domain, position=CORNER,scale=1.0) - call MOM_read_data(filename,trim(ubdryv_varname), u_bdry_val, G%Domain, position=CORNER,scale=1.0) - call MOM_read_data(filename,trim(vbdryv_varname), v_bdry_val, G%Domain, position=CORNER,scale=1.) - call MOM_read_data(filename,trim(umask_varname), umask, G%Domain, position=CORNER,scale=1.) - call MOM_read_data(filename,trim(vmask_varname), vmask, G%Domain, position=CORNER,scale=1.) - filename = trim(inputdir)//trim(icethick_file) + call MOM_read_data(filename, trim(ufcmskbdry_varname), u_face_mask_bdry, G%Domain, position=CORNER, scale=1.0) + call MOM_read_data(filename, trim(vfcmskbdry_varname), v_face_mask_bdry, G%Domain, position=CORNER, scale=1.0) + !### I think that the following two lines should have ..., scale=US%m_s_to_L_T + call MOM_read_data(filename, trim(ubdryv_varname), u_bdry_val, G%Domain, position=CORNER, scale=1.0) + call MOM_read_data(filename, trim(vbdryv_varname), v_bdry_val, G%Domain, position=CORNER, scale=1.) + call MOM_read_data(filename, trim(umask_varname), umask, G%Domain, position=CORNER, scale=1.) + call MOM_read_data(filename, trim(vmask_varname), vmask, G%Domain, position=CORNER, scale=1.) + filename = trim(inputdir)//trim(icethick_file) ! call MOM_read_data(filename, trim(h_varname), h_shelf, G%Domain, scale=US%m_to_Z) - call MOM_read_data(filename,trim(hmsk_varname), hmask, G%Domain, scale=1.) + call MOM_read_data(filename,trim(hmsk_varname), hmask, G%Domain, scale=1.) isc = G%isc ; jsc = G%jsc ; iec = G%iec ; jec = G%jec do j=jsc,jec diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 690345f988..a63a5f2e2d 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -533,13 +533,13 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & "MOM6 attempted to restart from a file from a different time than given by Time_in.") Time = Time_in endif - if ((GV%m_to_H_restart /= 0.0) .and. (GV%m_to_H_restart /= GV%m_to_H)) then - H_rescale = GV%m_to_H / GV%m_to_H_restart + if ((GV%m_to_H_restart /= 0.0) .and. (GV%m_to_H_restart /= 1.0)) then + H_rescale = 1.0 / GV%m_to_H_restart do k=1,nz ; do j=js,je ; do i=is,ie ; h(i,j,k) = H_rescale * h(i,j,k) ; enddo ; enddo ; enddo endif if ( (US%s_to_T_restart * US%m_to_L_restart /= 0.0) .and. & - ((US%m_to_L * US%s_to_T_restart) /= (US%m_to_L_restart * US%s_to_T)) ) then - vel_rescale = (US%m_to_L * US%s_to_T_restart) / (US%m_to_L_restart * US%s_to_T) + (US%s_to_T_restart /= US%m_to_L_restart) ) then + vel_rescale = US%s_to_T_restart / US%m_to_L_restart do k=1,nz ; do j=jsd,jed ; do I=IsdB,IeDB ; u(I,j,k) = vel_rescale * u(I,j,k) ; enddo ; enddo ; enddo do k=1,nz ; do J=JsdB,JedB ; do i=isd,ied ; v(i,J,k) = vel_rescale * v(i,J,k) ; enddo ; enddo ; enddo endif @@ -2778,7 +2778,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just do k=1,nz nPoints = 0 ; tempAvg = 0. ; saltAvg = 0. - do j=js,je ; do i=is,ie ; if (G%mask2dT(i,j) >= 1.0) then + do j=js,je ; do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then nPoints = nPoints + 1 tempAvg = tempAvg + tv%T(i,j,k) saltAvg = saltAvg + tv%S(i,j,k) @@ -2786,6 +2786,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just ! Horizontally homogenize data to produce perfectly "flat" initial conditions if (homogenize) then + !### These averages will not reproduce across PE layouts or grid rotation. call sum_across_PEs(nPoints) call sum_across_PEs(tempAvg) call sum_across_PEs(saltAvg) diff --git a/src/ocean_data_assim/MOM_oda_incupd.F90 b/src/ocean_data_assim/MOM_oda_incupd.F90 index af9534e943..5f38fa15d0 100644 --- a/src/ocean_data_assim/MOM_oda_incupd.F90 +++ b/src/ocean_data_assim/MOM_oda_incupd.F90 @@ -65,8 +65,8 @@ module MOM_oda_incupd !! registered by calls to set_up_oda_incupd_field type(p3d) :: Inc(MAX_FIELDS_) !< The increments to be applied to the field - type(p3d) :: Inc_u !< The increments to be applied to the u-velocities - type(p3d) :: Inc_v !< The increments to be applied to the v-velocities + type(p3d) :: Inc_u !< The increments to be applied to the u-velocities, with data in [L T-1 ~> m s-1] + type(p3d) :: Inc_v !< The increments to be applied to the v-velocities, with data in [L T-1 ~> m s-1] type(p3d) :: Ref_h !< Vertical grid on which the increments are provided @@ -99,9 +99,6 @@ subroutine initialize_oda_incupd_fixed( G, GV, US, CS, restart_CS) !! structure for this module (in/out). type(MOM_restart_CS), intent(inout) :: restart_CS !< MOM restart control struct -! This include declares and sets the variable "version". -#include "version_variable.h" - character(len=256) :: mesg if (associated(CS)) then call MOM_error(WARNING, "initialize_oda_incupd_fixed called with an associated "// & "control structure.") @@ -119,7 +116,7 @@ end subroutine initialize_oda_incupd_fixed !> This subroutine defined the number of time step for full update, stores the layer pressure !! increments and initialize remap structure. -subroutine initialize_oda_incupd( G, GV, US, param_file, CS, data_h,nz_data, restart_CS) +subroutine initialize_oda_incupd( G, GV, US, param_file, CS, data_h, nz_data, restart_CS) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -132,8 +129,8 @@ subroutine initialize_oda_incupd( G, GV, US, param_file, CS, data_h,nz_data, res !! [H ~> m or kg m-2]. type(MOM_restart_CS), intent(in) :: restart_CS !< MOM restart control struct -! This include declares and sets the variable "version". -#include "version_variable.h" + ! This include declares and sets the variable "version". +# include "version_variable.h" character(len=40) :: mdl = "MOM_oda" ! This module's name. logical :: use_oda_incupd logical :: bndExtrapolation = .true. ! If true, extrapolate boundaries @@ -231,6 +228,7 @@ subroutine initialize_oda_incupd( G, GV, US, param_file, CS, data_h,nz_data, res do j=G%jsc,G%jec; do i=G%isc,G%iec ; do k=1,CS%nz_data CS%Ref_h%p(i,j,k) = data_h(i,j,k) enddo; enddo ; enddo + !### Doing a halo update here on CS%Ref_h%p would avoid needing halo updates each timestep. ! Call the constructor for remapping control structure call initialize_remapping(CS%remap_cs, remapScheme, boundary_extrapolation=bndExtrapolation, & @@ -329,16 +327,17 @@ subroutine calc_oda_increments(h, tv, u, v, G, GV, US, CS) real, dimension(SZK_(GV)) :: tmp_val1 ! data values on the model grid real, allocatable, dimension(:) :: tmp_val2 ! data values remapped to increment grid - real, allocatable, dimension(:,:,:) :: h_obs !< h of increments - real, allocatable, dimension(:) :: tmp_h ! temporary array for corrected h_obs - real, allocatable, dimension(:) :: hu_obs,hv_obs ! A column of thicknesses at h points [H ~> m or kg m-2] - real, dimension(SZK_(GV)) :: hu, hv ! A column of thicknesses at h, u or v points [H ~> m or kg m-2] + real, allocatable, dimension(:,:,:) :: h_obs !< Layer-thicknesses of increments [H ~> m or kg m-2] + real, allocatable, dimension(:) :: tmp_h ! temporary array for corrected h_obs [H ~> m or kg m-2] + real, allocatable, dimension(:) :: hu_obs ! A column of observation-grid thicknesses at u points [H ~> m or kg m-2] + real, allocatable, dimension(:) :: hv_obs ! A column of observation-grid thicknesses at v points [H ~> m or kg m-2] + real, dimension(SZK_(GV)) :: hu, hv ! A column of thicknesses at u or v points [H ~> m or kg m-2] integer :: i, j, k, is, ie, js, je, nz, nz_data integer :: isB, ieB, jsB, jeB - real :: h_neglect, h_neglect_edge - real :: sum_h1, sum_h2 !vertical sums of h's + real :: h_neglect, h_neglect_edge ! Negligible thicknesses [H ~> m or kg m-2] + real :: sum_h1, sum_h2 ! vertical sums of h's [H ~> m or kg m-2] character(len=256) :: mesg is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke @@ -527,23 +526,24 @@ subroutine apply_oda_incupd(h, tv, u, v, dt, G, GV, US, CS) real :: m_to_Z ! A unit conversion factor from m to Z. real, allocatable, dimension(:) :: tmp_val2 ! data values on the increment grid real, dimension(SZK_(GV)) :: tmp_val1 ! data values remapped to model grid - real, dimension(SZK_(GV)) :: hu, hv ! A column of thicknesses at h, u or v points [H ~> m or kg m-2] + real, dimension(SZK_(GV)) :: hu, hv ! A column of thicknesses at u or v points [H ~> m or kg m-2] - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: tmp_t !< A temporary array for t inc. - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: tmp_s !< A temporary array for s inc. - real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: tmp_u !< A temporary array for u inc. - real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: tmp_v !< A temporary array for v inc. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: tmp_t !< A temporary array for t increments [degC] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: tmp_s !< A temporary array for s increments [ppt] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: tmp_u !< A temporary array for u increments [L T-1 ~> m s-1] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: tmp_v !< A temporary array for v increments [L T-1 ~> m s-1] real, allocatable, dimension(:,:,:) :: h_obs !< h of increments real, allocatable, dimension(:) :: tmp_h !< temporary array for corrected h_obs - real, allocatable, dimension(:) :: hu_obs,hv_obs ! A column of thicknesses at h points [H ~> m or kg m-2] + real, allocatable, dimension(:) :: hu_obs ! A column of observation-grid thicknesses at u points [H ~> m or kg m-2] + real, allocatable, dimension(:) :: hv_obs ! A column of observation-grid thicknesses at v points [H ~> m or kg m-2] integer :: i, j, k, is, ie, js, je, nz, nz_data integer :: isB, ieB, jsB, jeB ! integer :: ncount ! time step counter - real :: inc_wt ! weight of the update for this time-step - real :: h_neglect, h_neglect_edge - real :: sum_h1, sum_h2 !vertical sums of h's + real :: inc_wt ! weight of the update for this time-step [nondim] + real :: h_neglect, h_neglect_edge ! Negligible thicknesses [H ~> m or kg m-2] + real :: sum_h1, sum_h2 ! vertical sums of h's [H ~> m or kg m-2] character(len=256) :: mesg is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke @@ -735,7 +735,7 @@ subroutine apply_oda_incupd(h, tv, u, v, dt, G, GV, US, CS) deallocate(tmp_h,tmp_val2,hu_obs,hv_obs) deallocate(h_obs) - end subroutine apply_oda_incupd +end subroutine apply_oda_incupd !> Output increment if using full fields for the oda_incupd module. subroutine output_oda_incupd_inc(Time, G, GV, param_file, CS, US) @@ -771,12 +771,12 @@ subroutine output_oda_incupd_inc(Time, G, GV, param_file, CS, US) call register_restart_field(CS%Inc(2)%p, "S_inc", .true., restart_CSp_tmp, & "Salinity increment", "psu") call register_restart_field(CS%Ref_h%p, "h_obs", .true., restart_CSp_tmp, & - "Observational h", "m") + "Observational h", units=get_thickness_units(GV), conversion=GV%H_to_MKS) if (CS%uv_inc) then u_desc = var_desc("u_inc", "m s-1", "U-vel increment", hor_grid='Cu') v_desc = var_desc("v_inc", "m s-1", "V-vel increment", hor_grid='Cv') call register_restart_pair(CS%Inc_u%p, CS%Inc_v%p, u_desc, v_desc, & - .false., restart_CSp_tmp) + .false., restart_CSp_tmp, conversion=US%L_T_to_m_s) endif ! get the name of the output file diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index c786f395cf..80054ffff9 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -14,7 +14,6 @@ module MOM_MEKE use MOM_file_parser, only : read_param, get_param, log_version, param_file_type use MOM_grid, only : ocean_grid_type use MOM_hor_index, only : hor_index_type -use MOM_io, only : vardesc, var_desc use MOM_restart, only : MOM_restart_CS, register_restart_field, query_initialized use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : vertvisc_type @@ -1326,16 +1325,16 @@ logical function MEKE_init(Time, G, US, param_file, diag, CS, MEKE, restart_CS) ! Account for possible changes in dimensional scaling for variables that have been ! read from a restart file. I_T_rescale = 1.0 - if ((US%s_to_T_restart /= 0.0) .and. (US%s_to_T_restart /= US%s_to_T)) & - I_T_rescale = US%s_to_T_restart / US%s_to_T + if ((US%s_to_T_restart /= 0.0) .and. (US%s_to_T_restart /= 1.0)) & + I_T_rescale = US%s_to_T_restart L_rescale = 1.0 - if ((US%m_to_L_restart /= 0.0) .and. (US%m_to_L_restart /= US%m_to_L)) & - L_rescale = US%m_to_L / US%m_to_L_restart + if ((US%m_to_L_restart /= 0.0) .and. (US%m_to_L_restart /= 1.0)) & + L_rescale = 1.0 / US%m_to_L_restart if (L_rescale*I_T_rescale /= 1.0) then if (allocated(MEKE%MEKE)) then ; if (query_initialized(MEKE%MEKE, "MEKE_MEKE", restart_CS)) then do j=js,je ; do i=is,ie - MEKE%MEKE(i,j) = L_rescale*I_T_rescale * MEKE%MEKE(i,j) + MEKE%MEKE(i,j) = (L_rescale*I_T_rescale)**2 * MEKE%MEKE(i,j) enddo ; enddo endif ; endif endif @@ -1380,15 +1379,17 @@ logical function MEKE_init(Time, G, US, param_file, diag, CS, MEKE, restart_CS) end function MEKE_init !> Allocates memory and register restart fields for the MOM_MEKE module. -subroutine MEKE_alloc_register_restart(HI, param_file, MEKE, restart_CS) +subroutine MEKE_alloc_register_restart(HI, US, param_file, MEKE, restart_CS) ! Arguments type(hor_index_type), intent(in) :: HI !< Horizontal index structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(param_file_type), intent(in) :: param_file !< Parameter file parser structure. type(MEKE_type), intent(inout) :: MEKE !< MEKE fields type(MOM_restart_CS), intent(inout) :: restart_CS !< MOM restart control struct -! Local variables - type(vardesc) :: vd - real :: MEKE_GMcoeff, MEKE_FrCoeff, MEKE_GMECoeff, MEKE_KHCoeff, MEKE_viscCoeff_Ku, MEKE_viscCoeff_Au + + ! Local variables + real :: MEKE_GMcoeff, MEKE_FrCoeff, MEKE_GMECoeff ! Coefficients for various terms [nondim] + real :: MEKE_KHCoeff, MEKE_viscCoeff_Ku, MEKE_viscCoeff_Au ! Coefficients for various terms [nondim] logical :: Use_KH_in_MEKE logical :: useMEKE integer :: isd, ied, jsd, jed @@ -1397,13 +1398,13 @@ subroutine MEKE_alloc_register_restart(HI, param_file, MEKE, restart_CS) useMEKE = .false.; call read_param(param_file,"USE_MEKE",useMEKE) ! Read these parameters to determine what should be in the restarts - MEKE_GMcoeff =-1.; call read_param(param_file,"MEKE_GMCOEFF",MEKE_GMcoeff) - MEKE_FrCoeff =-1.; call read_param(param_file,"MEKE_FRCOEFF",MEKE_FrCoeff) - MEKE_GMEcoeff =-1.; call read_param(param_file,"MEKE_GMECOEFF",MEKE_GMEcoeff) - MEKE_KhCoeff =1.; call read_param(param_file,"MEKE_KHCOEFF",MEKE_KhCoeff) - MEKE_viscCoeff_Ku =0.; call read_param(param_file,"MEKE_VISCOSITY_COEFF_KU",MEKE_viscCoeff_Ku) - MEKE_viscCoeff_Au =0.; call read_param(param_file,"MEKE_VISCOSITY_COEFF_AU",MEKE_viscCoeff_Au) - Use_KH_in_MEKE = .false.; call read_param(param_file,"USE_KH_IN_MEKE", Use_KH_in_MEKE) + MEKE_GMcoeff = -1. ; call read_param(param_file,"MEKE_GMCOEFF",MEKE_GMcoeff) + MEKE_FrCoeff = -1. ; call read_param(param_file,"MEKE_FRCOEFF",MEKE_FrCoeff) + MEKE_GMEcoeff = -1. ; call read_param(param_file,"MEKE_GMECOEFF",MEKE_GMEcoeff) + MEKE_KhCoeff = 1. ; call read_param(param_file,"MEKE_KHCOEFF",MEKE_KhCoeff) + MEKE_viscCoeff_Ku = 0. ; call read_param(param_file,"MEKE_VISCOSITY_COEFF_KU",MEKE_viscCoeff_Ku) + MEKE_viscCoeff_Au = 0. ; call read_param(param_file,"MEKE_VISCOSITY_COEFF_AU",MEKE_viscCoeff_Au) + Use_KH_in_MEKE = .false. ; call read_param(param_file,"USE_KH_IN_MEKE", Use_KH_in_MEKE) if (.not. useMEKE) return @@ -1411,38 +1412,38 @@ subroutine MEKE_alloc_register_restart(HI, param_file, MEKE, restart_CS) call MOM_mesg("MEKE_alloc_register_restart: allocating and registering", 5) isd = HI%isd ; ied = HI%ied ; jsd = HI%jsd ; jed = HI%jed allocate(MEKE%MEKE(isd:ied,jsd:jed), source=0.0) - vd = var_desc("MEKE", "m2 s-2", hor_grid='h', z_grid='1', & - longname="Mesoscale Eddy Kinetic Energy") - call register_restart_field(MEKE%MEKE, vd, .false., restart_CS) + call register_restart_field(MEKE%MEKE, "MEKE", .false., restart_CS, & + longname="Mesoscale Eddy Kinetic Energy", units="m2 s-2", conversion=US%L_T_to_m_s**2) + if (MEKE_GMcoeff>=0.) allocate(MEKE%GM_src(isd:ied,jsd:jed), source=0.0) if (MEKE_FrCoeff>=0. .or. MEKE_GMECoeff>=0.) & allocate(MEKE%mom_src(isd:ied,jsd:jed), source=0.0) if (MEKE_GMECoeff>=0.) allocate(MEKE%GME_snk(isd:ied,jsd:jed), source=0.0) if (MEKE_KhCoeff>=0.) then allocate(MEKE%Kh(isd:ied,jsd:jed), source=0.0) - vd = var_desc("MEKE_Kh", "m2 s-1", hor_grid='h', z_grid='1', & - longname="Lateral diffusivity from Mesoscale Eddy Kinetic Energy") - call register_restart_field(MEKE%Kh, vd, .false., restart_CS) + call register_restart_field(MEKE%Kh, "MEKE_Kh", .false., restart_CS, & + longname="Lateral diffusivity from Mesoscale Eddy Kinetic Energy", & + units="m2 s-1", conversion=US%L_to_m**2*US%s_to_T) endif allocate(MEKE%Rd_dx_h(isd:ied,jsd:jed), source=0.0) if (MEKE_viscCoeff_Ku/=0.) then allocate(MEKE%Ku(isd:ied,jsd:jed), source=0.0) - vd = var_desc("MEKE_Ku", "m2 s-1", hor_grid='h', z_grid='1', & - longname="Lateral viscosity from Mesoscale Eddy Kinetic Energy") - call register_restart_field(MEKE%Ku, vd, .false., restart_CS) + call register_restart_field(MEKE%Ku, "MEKE_Ku", .false., restart_CS, & + longname="Lateral viscosity from Mesoscale Eddy Kinetic Energy", & + units="m2 s-1", conversion=US%L_to_m**2*US%s_to_T) endif if (Use_Kh_in_MEKE) then allocate(MEKE%Kh_diff(isd:ied,jsd:jed), source=0.0) - vd = var_desc("MEKE_Kh_diff", "m2 s-1",hor_grid='h',z_grid='1', & - longname="Copy of thickness diffusivity for diffusing MEKE") - call register_restart_field(MEKE%Kh_diff, vd, .false., restart_CS) + call register_restart_field(MEKE%Kh_diff, "MEKE_Kh_diff", .false., restart_CS, & + longname="Copy of thickness diffusivity for diffusing MEKE", & + units="m2 s-1", conversion=US%L_to_m**2*US%s_to_T) endif if (MEKE_viscCoeff_Au/=0.) then allocate(MEKE%Au(isd:ied,jsd:jed), source=0.0) - vd = var_desc("MEKE_Au", "m4 s-1", hor_grid='h', z_grid='1', & - longname="Lateral biharmonic viscosity from Mesoscale Eddy Kinetic Energy") - call register_restart_field(MEKE%Au, vd, .false., restart_CS) + call register_restart_field(MEKE%Au, "MEKE_Au", .false., restart_CS, & + longname="Lateral biharmonic viscosity from Mesoscale Eddy Kinetic Energy", & + units="m4 s-1", conversion=US%L_to_m**4*US%s_to_T) endif end subroutine MEKE_alloc_register_restart diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index dfbb3e0d63..7045a1a52c 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/src/parameterizations/lateral/MOM_internal_tides.F90 @@ -16,7 +16,7 @@ module MOM_internal_tides use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg, is_root_pe use MOM_file_parser, only : read_param, get_param, log_param, log_version, param_file_type use MOM_grid, only : ocean_grid_type -use MOM_io, only : slasher, vardesc, MOM_read_data, file_exists +use MOM_io, only : slasher, MOM_read_data, file_exists use MOM_restart, only : register_restart_field, MOM_restart_CS, restart_init, save_restart use MOM_spatial_means, only : global_area_mean use MOM_time_manager, only : time_type, time_type_to_real, operator(+), operator(/), operator(-) @@ -2086,7 +2086,6 @@ end subroutine PPM_limit_pos ! ! This subroutine is used to allocate and register any fields in this module ! ! that should be written to or read from the restart file. ! logical :: use_int_tides -! type(vardesc) :: vd ! integer :: num_freq, num_angle , num_mode, period_1 ! integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, a ! isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed @@ -2108,10 +2107,9 @@ end subroutine PPM_limit_pos ! call read_param(param_file, "INTERNAL_TIDE_ANGLES", num_angle) ! allocate(CS%En_restart(isd:ied, jsd:jed, num_angle), source=0.0) -! vd = vardesc("En_restart", & -! "The internal wave energy density as a function of (i,j,angle,frequency,mode)", & -! 'h','1','1',"J m-2") -! call register_restart_field(CS%En_restart, vd, .false., restart_CS) +! call register_restart_field(CS%En_restart, "En_restart", .false., restart_CS, & +! longname="The internal wave energy density as a function of (i,j,angle,frequency,mode)", & +! units="J m-2", conversion=US%RZ3_T3_to_W_m2*US%T_to_s, z_grid='1') ! end subroutine register_int_tide_restarts diff --git a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 index 6176f4aa1d..d8da752a1f 100644 --- a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 +++ b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 @@ -13,12 +13,11 @@ module MOM_mixed_layer_restrat use MOM_forcing_type, only : mech_forcing use MOM_grid, only : ocean_grid_type use MOM_hor_index, only : hor_index_type -use MOM_io, only : vardesc, var_desc use MOM_lateral_mixing_coeffs, only : VarMix_CS use MOM_restart, only : register_restart_field, query_initialized, MOM_restart_CS use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : thermo_var_ptrs -use MOM_verticalGrid, only : verticalGrid_type +use MOM_verticalGrid, only : verticalGrid_type, get_thickness_units use MOM_EOS, only : calculate_density, EOS_domain implicit none ; private @@ -922,8 +921,8 @@ logical function mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, CS, ! Rescale variables from restart files if the internal dimensional scalings have changed. if (CS%MLE_MLD_decay_time>0. .or. CS%MLE_MLD_decay_time2>0.) then if (query_initialized(CS%MLD_filtered, "MLD_MLE_filtered", restart_CS) .and. & - (GV%m_to_H_restart /= 0.0) .and. (GV%m_to_H_restart /= GV%m_to_H)) then - H_rescale = GV%m_to_H / GV%m_to_H_restart + (GV%m_to_H_restart /= 0.0) .and. (GV%m_to_H_restart /= 1.0)) then + H_rescale = 1.0 / GV%m_to_H_restart do j=G%jsc,G%jec ; do i=G%isc,G%iec CS%MLD_filtered(i,j) = H_rescale * CS%MLD_filtered(i,j) enddo ; enddo @@ -931,8 +930,8 @@ logical function mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, CS, endif if (CS%MLE_MLD_decay_time2>0.) then if (query_initialized(CS%MLD_filtered_slow, "MLD_MLE_filtered_slow", restart_CS) .and. & - (GV%m_to_H_restart /= 0.0) .and. (GV%m_to_H_restart /= GV%m_to_H)) then - H_rescale = GV%m_to_H / GV%m_to_H_restart + (GV%m_to_H_restart /= 0.0) .and. (GV%m_to_H_restart /= 1.0)) then + H_rescale = 1.0 / GV%m_to_H_restart do j=G%jsc,G%jec ; do i=G%isc,G%iec CS%MLD_filtered_slow(i,j) = H_rescale * CS%MLD_filtered_slow(i,j) enddo ; enddo @@ -945,14 +944,15 @@ logical function mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, CS, end function mixedlayer_restrat_init !> Allocate and register fields in the mixed layer restratification structure for restarts -subroutine mixedlayer_restrat_register_restarts(HI, param_file, CS, restart_CS) +subroutine mixedlayer_restrat_register_restarts(HI, GV, param_file, CS, restart_CS) ! Arguments type(hor_index_type), intent(in) :: HI !< Horizontal index structure + type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure type(param_file_type), intent(in) :: param_file !< Parameter file to parse type(mixedlayer_restrat_CS), intent(inout) :: CS !< Module control structure type(MOM_restart_CS), intent(inout) :: restart_CS !< MOM restart control struct + ! Local variables - type(vardesc) :: vd logical :: mixedlayer_restrat_init ! Check to see if this module will be used @@ -967,16 +967,16 @@ subroutine mixedlayer_restrat_register_restarts(HI, param_file, CS, restart_CS) if (CS%MLE_MLD_decay_time>0. .or. CS%MLE_MLD_decay_time2>0.) then ! CS%MLD_filtered is used to keep a running mean of the PBL's actively mixed MLD. allocate(CS%MLD_filtered(HI%isd:HI%ied,HI%jsd:HI%jed), source=0.) - vd = var_desc("MLD_MLE_filtered","m","Time-filtered MLD for use in MLE", & - hor_grid='h', z_grid='1') - call register_restart_field(CS%MLD_filtered, vd, .false., restart_CS) + call register_restart_field(CS%MLD_filtered, "MLD_MLE_filtered", .false., restart_CS, & + longname="Time-filtered MLD for use in MLE", & + units=get_thickness_units(GV), conversion=GV%H_to_MKS) endif if (CS%MLE_MLD_decay_time2>0.) then ! CS%MLD_filtered_slow is used to keep a running mean of the PBL's seasonal or winter MLD. allocate(CS%MLD_filtered_slow(HI%isd:HI%ied,HI%jsd:HI%jed), source=0.) - vd = var_desc("MLD_MLE_filtered_slow","m","c Slower time-filtered MLD for use in MLE", & - hor_grid='h', z_grid='1') - call register_restart_field(CS%MLD_filtered_slow, vd, .false., restart_CS) + call register_restart_field(CS%MLD_filtered, "MLD_MLE_filtered_slow", .false., restart_CS, & + longname="Slower time-filtered MLD for use in MLE", & + units=get_thickness_units(GV), conversion=GV%H_to_MKS) endif end subroutine mixedlayer_restrat_register_restarts diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index fb969953c4..92e7b44128 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -1807,9 +1807,10 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS) end subroutine set_viscous_ML !> Register any fields associated with the vertvisc_type. -subroutine set_visc_register_restarts(HI, GV, param_file, visc, restart_CS) +subroutine set_visc_register_restarts(HI, GV, US, param_file, visc, restart_CS) type(hor_index_type), intent(in) :: HI !< A horizontal index type 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 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time !! parameters. type(vertvisc_type), intent(inout) :: visc !< A structure containing vertical @@ -1849,20 +1850,22 @@ subroutine set_visc_register_restarts(HI, GV, param_file, visc, restart_CS) if (use_kappa_shear .or. useKPP .or. useEPBL .or. use_CVMix_shear .or. use_CVMix_conv) then call safe_alloc_ptr(visc%Kd_shear, isd, ied, jsd, jed, nz+1) call register_restart_field(visc%Kd_shear, "Kd_shear", .false., restart_CS, & - "Shear-driven turbulent diffusivity at interfaces", "m2 s-1", z_grid='i') + "Shear-driven turbulent diffusivity at interfaces", & + units="m2 s-1", conversion=US%Z2_T_to_m2_s, z_grid='i') endif if (useKPP .or. useEPBL .or. use_CVMix_shear .or. use_CVMix_conv .or. & (use_kappa_shear .and. .not.KS_at_vertex )) then call safe_alloc_ptr(visc%Kv_shear, isd, ied, jsd, jed, nz+1) call register_restart_field(visc%Kv_shear, "Kv_shear", .false., restart_CS, & - "Shear-driven turbulent viscosity at interfaces", "m2 s-1", z_grid='i') + "Shear-driven turbulent viscosity at interfaces", & + units="m2 s-1", conversion=US%Z2_T_to_m2_s, z_grid='i') endif if (use_kappa_shear .and. KS_at_vertex) then call safe_alloc_ptr(visc%TKE_turb, HI%IsdB, HI%IedB, HI%JsdB, HI%JedB, nz+1) call safe_alloc_ptr(visc%Kv_shear_Bu, HI%IsdB, HI%IedB, HI%JsdB, HI%JedB, nz+1) call register_restart_field(visc%Kv_shear_Bu, "Kv_shear_Bu", .false., restart_CS, & - "Shear-driven turbulent viscosity at vertex interfaces", "m2 s-1", & - hor_grid="Bu", z_grid='i') + "Shear-driven turbulent viscosity at vertex interfaces", & + units="m2 s-1", conversion=US%Z2_T_to_m2_s, hor_grid="Bu", z_grid='i') elseif (use_kappa_shear) then call safe_alloc_ptr(visc%TKE_turb, isd, ied, jsd, jed, nz+1) endif @@ -1882,7 +1885,7 @@ subroutine set_visc_register_restarts(HI, GV, param_file, visc, restart_CS) if (MLE_use_PBL_MLD) then call safe_alloc_ptr(visc%MLD, isd, ied, jsd, jed) call register_restart_field(visc%MLD, "MLD", .false., restart_CS, & - "Instantaneous active mixing layer depth", "m") + "Instantaneous active mixing layer depth", "m", conversion=US%Z_to_m) endif if (hfreeze >= 0.0 .and. .not.MLE_use_PBL_MLD) then @@ -2191,9 +2194,9 @@ subroutine set_visc_init(Time, G, GV, US, param_file, diag, visc, CS, restart_CS allocate(visc%nkml_visc_u(IsdB:IedB,jsd:jed), source=0.0) allocate(visc%nkml_visc_v(isd:ied,JsdB:JedB), source=0.0) CS%id_nkml_visc_u = register_diag_field('ocean_model', 'nkml_visc_u', & - diag%axesCu1, Time, 'Number of layers in viscous mixed layer at u points', 'm') + diag%axesCu1, Time, 'Number of layers in viscous mixed layer at u points', 'nondim') CS%id_nkml_visc_v = register_diag_field('ocean_model', 'nkml_visc_v', & - diag%axesCv1, Time, 'Number of layers in viscous mixed layer at v points', 'm') + diag%axesCv1, Time, 'Number of layers in viscous mixed layer at v points', 'nondim') endif call register_restart_field_as_obsolete('Kd_turb','Kd_shear', restart_CS) @@ -2202,11 +2205,9 @@ subroutine set_visc_init(Time, G, GV, US, param_file, diag, visc, CS, restart_CS ! Account for possible changes in dimensional scaling for variables that have been ! read from a restart file. Z_rescale = 1.0 - if ((US%m_to_Z_restart /= 0.0) .and. (US%m_to_Z_restart /= US%m_to_Z)) & - Z_rescale = US%m_to_Z / US%m_to_Z_restart + if (US%m_to_Z_restart /= 0.0) Z_rescale = 1.0 / US%m_to_Z_restart I_T_rescale = 1.0 - if ((US%s_to_T_restart /= 0.0) .and. (US%s_to_T_restart /= US%s_to_T)) & - I_T_rescale = US%s_to_T_restart / US%s_to_T + if (US%s_to_T_restart /= 0.0) I_T_rescale = US%s_to_T_restart Z2_T_rescale = Z_rescale**2*I_T_rescale if (Z2_T_rescale /= 1.0) then diff --git a/src/tracer/boundary_impulse_tracer.F90 b/src/tracer/boundary_impulse_tracer.F90 index 44423b5650..f85623bbd1 100644 --- a/src/tracer/boundary_impulse_tracer.F90 +++ b/src/tracer/boundary_impulse_tracer.F90 @@ -138,7 +138,7 @@ function register_boundary_impulse_tracer(HI, GV, US, param_file, CS, tr_Reg, re rem_time_ptr => CS%remaining_source_time call register_restart_field(rem_time_ptr, "bir_remain_time", & .not.CS%tracers_may_reinit, restart_CS, & - "Remaining time to apply BIR source", "s") + "Remaining time to apply BIR source", "s", conversion=US%T_to_s) CS%tr_Reg => tr_Reg CS%restart_CSp => restart_CS @@ -196,8 +196,8 @@ subroutine initialize_boundary_impulse_tracer(restart, day, G, GV, US, h, diag, endif enddo ! Tracer loop - if (restart .and. (US%s_to_T_restart /= 0.0) .and. (US%s_to_T /= US%s_to_T_restart) ) then - CS%remaining_source_time = (US%s_to_T / US%s_to_T_restart) * CS%remaining_source_time + if (restart .and. (US%s_to_T_restart /= 0.0) .and. (US%s_to_T_restart /= 1.0) ) then + CS%remaining_source_time = (1.0 / US%s_to_T_restart) * CS%remaining_source_time endif if (associated(OBC)) then diff --git a/src/user/MOM_controlled_forcing.F90 b/src/user/MOM_controlled_forcing.F90 index d136d58a19..7583485ad7 100644 --- a/src/user/MOM_controlled_forcing.F90 +++ b/src/user/MOM_controlled_forcing.F90 @@ -427,8 +427,9 @@ function periodic_real(rval, num_period) result(val_out) !> This subroutine is used to allocate and register any fields in this module !! that should be written to or read from the restart file. -subroutine register_ctrl_forcing_restarts(G, param_file, CS, restart_CS) +subroutine register_ctrl_forcing_restarts(G, US, param_file, CS, restart_CS) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(param_file_type), intent(in) :: param_file !< A structure indicating the !! open file to parse for model !! parameter values. @@ -469,9 +470,11 @@ subroutine register_ctrl_forcing_restarts(G, param_file, CS, restart_CS) allocate(CS%precip_0(isd:ied,jsd:jed), source=0.0) call register_restart_field(CS%heat_0, "Ctrl_heat", .false., restart_CS, & - longname="Control Integrative Heating", units="W m-2", z_grid='1') + longname="Control Integrative Heating", & + units="W m-2", conversion=US%QRZ_T_to_W_m2, z_grid='1') call register_restart_field(CS%precip_0, "Ctrl_precip", .false., restart_CS, & - longname="Control Integrative Precipitation", units="kg m-2 s-1", z_grid='1') + longname="Control Integrative Precipitation", & + units="kg m-2 s-1", conversion=US%RZ_T_to_kg_m2s, z_grid='1') endif if (CS%num_cycle > 0) then @@ -485,13 +488,16 @@ subroutine register_ctrl_forcing_restarts(G, param_file, CS, restart_CS) period_str = trim('p ')//trim(adjustl(period_str)) call register_restart_field(CS%heat_cyc, "Ctrl_heat_cycle", .false., restart_CS, & - longname="Cyclical Control Heating", units="W m-2", z_grid='1', t_grid=period_str) + longname="Cyclical Control Heating", & + units="W m-2", conversion=US%QRZ_T_to_W_m2, z_grid='1', t_grid=period_str) call register_restart_field(CS%precip_cyc, "Ctrl_precip_cycle", .false., restart_CS, & - longname="Cyclical Control Precipitation", units="kg m-2 s-1", z_grid='1', t_grid=period_str) + longname="Cyclical Control Precipitation", & + units="kg m-2 s-1", conversion=US%RZ_T_to_kg_m2s, z_grid='1', t_grid=period_str) call register_restart_field(CS%avg_time, "avg_time", .false., restart_CS, & - longname="Cyclical accumulated averaging time", units="sec", z_grid='1', t_grid=period_str) + longname="Cyclical accumulated averaging time", & + units="sec", conversion=US%T_to_s, z_grid='1', t_grid=period_str) call register_restart_field(CS%avg_SST_anom, "avg_SST_anom", .false., restart_CS, & - longname="Cyclical average SST Anomaly", units="deg C", z_grid='1', t_grid=period_str) + longname="Cyclical average SST Anomaly", units="degC", z_grid='1', t_grid=period_str) call register_restart_field(CS%avg_SSS_anom, "avg_SSS_anom", .false., restart_CS, & longname="Cyclical average SSS Anomaly", units="g kg-1", z_grid='1', t_grid=period_str) endif @@ -592,11 +598,9 @@ subroutine controlled_forcing_init(Time, G, US, param_file, diag, CS) ! Rescale if there are differences between the dimensional scaling of variables in ! restart files from those in use for this run. if ((US%J_kg_to_Q_restart*US%kg_m3_to_R_restart*US%m_to_Z_restart*US%s_to_T_restart /= 0.0) .and. & - ((US%J_kg_to_Q * US%kg_m3_to_R * US%m_to_Z * US%s_to_T_restart) /= & - (US%J_kg_to_Q_restart * US%kg_m3_to_R_restart * US%m_to_Z_restart * US%s_to_T)) ) then + (US%s_to_T_restart /= US%J_kg_to_Q_restart * US%kg_m3_to_R_restart * US%m_to_Z_restart) ) then ! Redo the scaling of the corrective heat fluxes to [Q R Z T-1 ~> W m-2] - QRZ_T_rescale = (US%J_kg_to_Q * US%kg_m3_to_R * US%m_to_Z * US%s_to_T_restart) / & - (US%J_kg_to_Q_restart * US%kg_m3_to_R_restart * US%m_to_Z_restart * US%s_to_T) + QRZ_T_rescale = US%s_to_T_restart / (US%J_kg_to_Q_restart * US%kg_m3_to_R_restart * US%m_to_Z_restart) if (associated(CS%heat_0)) then do j=jsc,jec ; do i=isc,iec @@ -612,11 +616,9 @@ subroutine controlled_forcing_init(Time, G, US, param_file, diag, CS) endif if ((US%kg_m3_to_R_restart * US%m_to_Z_restart * US%s_to_T_restart /= 0.0) .and. & - ((US%kg_m3_to_R * US%m_to_Z * US%s_to_T_restart) /= & - (US%kg_m3_to_R_restart * US%m_to_Z_restart * US%s_to_T)) ) then + (US%s_to_T_restart /= US%kg_m3_to_R_restart * US%m_to_Z_restart) ) then ! Redo the scaling of the corrective precipitation to [R Z T-1 ~> kg m-2 s-1] - RZ_T_rescale = (US%kg_m3_to_R * US%m_to_Z * US%s_to_T_restart) / & - (US%kg_m3_to_R_restart * US%m_to_Z_restart * US%s_to_T) + RZ_T_rescale = US%s_to_T_restart / (US%kg_m3_to_R_restart * US%m_to_Z_restart) if (associated(CS%precip_0)) then do j=jsc,jec ; do i=isc,iec @@ -632,10 +634,10 @@ subroutine controlled_forcing_init(Time, G, US, param_file, diag, CS) endif if ((CS%num_cycle > 0) .and. associated(CS%avg_time) .and. & - ((US%s_to_T_restart /= 0.0) .and. ((US%s_to_T_restart) /= US%s_to_T)) ) then + ((US%s_to_T_restart /= 0.0) .and. (US%s_to_T_restart /= 1.0)) ) then ! Redo the scaling of the accumulated times to [T ~> s] do m=1,CS%num_cycle - CS%avg_time(m) = (US%s_to_T / US%s_to_T_restart) * CS%avg_time(m) + CS%avg_time(m) = (1.0 / US%s_to_T_restart) * CS%avg_time(m) enddo endif