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/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. diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index c1e5392d5e..bcb17589dc 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -1659,11 +1659,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", & @@ -1671,7 +1671,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 @@ -1688,13 +1688,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 @@ -1723,23 +1723,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 5eccf92976..608a672a46 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 751a6953e6..d3202a9f2a 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) @@ -407,7 +407,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] 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 @@ -453,13 +453,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 @@ -472,18 +473,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]. @@ -491,9 +493,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 @@ -546,15 +548,16 @@ 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(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 c03b3e7453..d337e039e2 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 @@ -2791,7 +2791,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) @@ -2799,6 +2799,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 6c276fe8d0..f6d8d05dcf 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(-) @@ -2105,7 +2105,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 @@ -2127,10 +2126,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