From 2a5e98a538cd78df54b0e4daa43a4f0ffa3c22d3 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 11 Oct 2019 18:28:41 -0400 Subject: [PATCH] +Rescaled lengths in SIS_hor_grid_type Rescaled all of the lengths, inverse lengths, areas and inverse areas in the SIS_hor_grid_type, and added unit_scale_type elements to the fast and slow ice control structures and as a arguments to multiple routines. Answers in the Baltic test case are bitwise identical, but there are numerous public interface with added mandatory arguments. --- src/SIS_continuity.F90 | 161 ++++++++++++++++--------------- src/SIS_ctrl_types.F90 | 8 +- src/SIS_dyn_bgrid.F90 | 46 ++++----- src/SIS_dyn_cgrid.F90 | 102 +++++++++++--------- src/SIS_dyn_trans.F90 | 132 +++++++++++++------------ src/SIS_fixed_initialization.F90 | 16 +-- src/SIS_hor_grid.F90 | 59 +++++------ src/SIS_slow_thermo.F90 | 10 +- src/SIS_sum_output.F90 | 24 ++--- src/SIS_tracer_advect.F90 | 99 ++++++++++--------- src/SIS_transport.F90 | 59 ++++++----- src/SIS_utils.F90 | 42 ++++---- src/ice_age_tracer.F90 | 3 +- src/ice_model.F90 | 50 ++++++---- src/ice_type.F90 | 13 ++- src/slab_ice.F90 | 14 +-- src/specified_ice.F90 | 8 +- 17 files changed, 462 insertions(+), 384 deletions(-) diff --git a/src/SIS_continuity.F90 b/src/SIS_continuity.F90 index 6a0f8e73..bafdde0c 100644 --- a/src/SIS_continuity.F90 +++ b/src/SIS_continuity.F90 @@ -18,6 +18,7 @@ module SIS_continuity use SIS_diag_mediator, only : time_type, SIS_diag_ctrl use MOM_error_handler, only : SIS_error=>MOM_error, FATAL, WARNING, is_root_pe use MOM_file_parser, only : get_param, log_version, param_file_type +use MOM_unit_scaling, only : unit_scale_type use SIS_hor_grid, only : SIS_hor_grid_type use ice_grid, only : ice_grid_type ! use MOM_variables, only : ocean_OBC_type, OBC_SIMPLE @@ -63,7 +64,7 @@ module SIS_continuity !> ice_continuity time steps the category thickness changes due to advection, !! using a monotonically limited, directionally split PPM scheme. -subroutine ice_continuity(u, v, hin, h, uh, vh, dt, G, IG, CS) +subroutine ice_continuity(u, v, hin, h, uh, vh, dt, G, US, IG, CS) type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type type(ice_grid_type), intent(inout) :: IG !< The sea-ice specific grid type real, dimension(SZIB_(G),SZJ_(G)), & @@ -81,6 +82,7 @@ subroutine ice_continuity(u, v, hin, h, uh, vh, dt, G, IG, CS) intent(out) :: vh !< Volume flux through meridional faces = v*h*dx !! [H m2 s-1 ~> kg s-1]. real, intent(in) :: dt !< Time increment [s] + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors type(SIS_continuity_CS), pointer :: CS !< The control structure returned by a !! previous call to SIS_continuity_init. ! This subroutine time steps the category thicknesses, using a monotonically @@ -113,17 +115,17 @@ subroutine ice_continuity(u, v, hin, h, uh, vh, dt, G, IG, CS) do j=js,je ; do k=1,nCat ; do I=is-1,ie if (u(I,j) >= 0.0) then ; h_up = hin(i,j,k) else ; h_up = hin(i+1,j,k) ; endif - uh(I,j,k) = G%dy_Cu(I,j) * u(I,j) * h_up + uh(I,j,k) = US%L_to_m*G%dy_Cu(I,j) * u(I,j) * h_up enddo ; enddo ; enddo !$OMP do do J=js-1,je ; do k=1,nCat ; do i=is,ie if (v(i,J) >= 0.0) then ; h_up = hin(i,j,k) else ; h_up = hin(i,j+1,k) ; endif - vh(i,J,k) = G%dx_Cv(i,J) * v(i,J) * h_up + vh(i,J,k) = US%L_to_m*G%dx_Cv(i,J) * v(i,J) * h_up enddo ; enddo ; enddo !$OMP do do j=js,je ; do k=1,nCat ; do i=is,ie - h(i,j,k) = hin(i,j,k) - dt* G%IareaT(i,j) * & + h(i,j,k) = hin(i,j,k) - dt* US%m_to_L**2*G%IareaT(i,j) * & ((uh(I,j,k) - uh(I-1,j,k)) + (vh(i,J,k) - vh(i,J-1,k))) if (h(i,j,k) < 0.0) then @@ -135,12 +137,12 @@ subroutine ice_continuity(u, v, hin, h, uh, vh, dt, G, IG, CS) ! First, advect zonally. LB%ish = G%isc ; LB%ieh = G%iec LB%jsh = G%jsc-stensil ; LB%jeh = G%jec+stensil - call zonal_mass_flux(u, dt, G, IG, CS, LB, hin, uh) + call zonal_mass_flux(u, dt, G, US, IG, CS, LB, hin, uh) call cpu_clock_begin(id_clock_update) !$OMP parallel do default(shared) do j=LB%jsh,LB%jeh ; do k=1,nCat ; do i=LB%ish,LB%ieh - h(i,j,k) = hin(i,j,k) - dt* G%IareaT(i,j) * (uh(I,j,k) - uh(I-1,j,k)) + h(i,j,k) = hin(i,j,k) - dt* US%m_to_L**2*G%IareaT(i,j) * (uh(I,j,k) - uh(I-1,j,k)) if (h(i,j,k) < 0.0) then call SIS_error(FATAL, & 'Negative thickness encountered in u-pass of ice_continuity().') @@ -152,12 +154,12 @@ subroutine ice_continuity(u, v, hin, h, uh, vh, dt, G, IG, CS) ! Now advect meridionally, using the updated thicknesses to determine ! the fluxes. - call meridional_mass_flux(v, dt, G, IG, CS, LB, h, vh) + call meridional_mass_flux(v, dt, G, US, IG, CS, LB, h, vh) call cpu_clock_begin(id_clock_update) !$OMP parallel do default(shared) do j=LB%jsh,LB%jeh ; do k=1,nCat ; do i=LB%ish,LB%ieh - h(i,j,k) = h(i,j,k) - dt*G%IareaT(i,j) * (vh(i,J,k) - vh(i,J-1,k)) + h(i,j,k) = h(i,j,k) - dt*US%m_to_L**2*G%IareaT(i,j) * (vh(i,J,k) - vh(i,J-1,k)) if (h(i,j,k) < 0.0) then call SIS_error(FATAL, & 'Negative thickness encountered in v-pass of ice_continuity().') @@ -170,12 +172,12 @@ subroutine ice_continuity(u, v, hin, h, uh, vh, dt, G, IG, CS) LB%ish = G%isc-stensil ; LB%ieh = G%iec+stensil LB%jsh = G%jsc ; LB%jeh = G%jec - call meridional_mass_flux(v, dt, G, IG, CS, LB, hin, vh) + call meridional_mass_flux(v, dt, G, US, IG, CS, LB, hin, vh) call cpu_clock_begin(id_clock_update) !$OMP parallel do default(shared) do j=LB%jsh,LB%jeh ; do k=1,nCat ; do i=LB%ish,LB%ieh - h(i,j,k) = hin(i,j,k) - dt*G%IareaT(i,j) * (vh(i,J,k) - vh(i,J-1,k)) + h(i,j,k) = hin(i,j,k) - dt*US%m_to_L**2*G%IareaT(i,j) * (vh(i,J,k) - vh(i,J-1,k)) if (h(i,j,k) < 0.0) then call SIS_error(FATAL, & 'Negative thickness encountered in v-pass of ice_continuity().') @@ -186,12 +188,12 @@ subroutine ice_continuity(u, v, hin, h, uh, vh, dt, G, IG, CS) ! Now advect zonally, using the updated thicknesses to determine ! the fluxes. LB%ish = G%isc ; LB%ieh = G%iec ; LB%jsh = G%jsc ; LB%jeh = G%jec - call zonal_mass_flux(u, dt, G, IG, CS, LB, h, uh) + call zonal_mass_flux(u, dt, G, US, IG, CS, LB, h, uh) call cpu_clock_begin(id_clock_update) !$OMP parallel do default(shared) do j=LB%jsh,LB%jeh ; do k=1,nCat ; do i=LB%ish,LB%ieh - h(i,j,k) = h(i,j,k) - dt* G%IareaT(i,j) * (uh(I,j,k) - uh(I-1,j,k)) + h(i,j,k) = h(i,j,k) - dt* US%m_to_L**2*G%IareaT(i,j) * (uh(I,j,k) - uh(I-1,j,k)) if (h(i,j,k) < 0.0) then call SIS_error(FATAL, & 'Negative thickness encountered in u-pass of ice_continuity().') @@ -205,13 +207,14 @@ end subroutine ice_continuity !> ice_cover_transport advects the total fractional ice cover and limits them not to exceed 1. -subroutine ice_cover_transport(u, v, cvr, dt, G, IG, CS) +subroutine ice_cover_transport(u, v, cvr, dt, G, US, IG, CS) type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type type(ice_grid_type), intent(inout) :: IG !< The sea-ice specific grid type real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: u !< Zonal ice velocity [m s-1]. real, dimension(SZI_(G),SZJB_(G)), intent(in) :: v !< Meridional ice velocity [m s-1]. real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: cvr !< Fractional ice cover [nondim]. real, intent(in) :: dt !< Time increment [s] + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors type(SIS_continuity_CS), pointer :: CS !< The control structure returned by a !! previous call to SIS_continuity_init. @@ -243,17 +246,17 @@ subroutine ice_cover_transport(u, v, cvr, dt, G, IG, CS) do j=js,je ; do I=is-1,ie if (u(I,j) >= 0.0) then ; cvr_up = cvr(i,j) else ; cvr_up = cvr(i+1,j) ; endif - ucvr(I,j) = G%dy_Cu(I,j) * u(I,j) * cvr_up + ucvr(I,j) = US%L_to_m*G%dy_Cu(I,j) * u(I,j) * cvr_up enddo ; enddo !$OMP do do J=js-1,je ; do i=is,ie if (v(i,J) >= 0.0) then ; cvr_up = cvr(i,j) else ; cvr_up = cvr(i,j+1) ; endif - vcvr(i,J) = G%dx_Cv(i,J) * v(i,J) * cvr_up + vcvr(i,J) = US%L_to_m*G%dx_Cv(i,J) * v(i,J) * cvr_up enddo ; enddo !$OMP do do j=js,je ; do i=is,ie - cvr(i,j) = cvr(i,j) - dt* G%IareaT(i,j) * & + cvr(i,j) = cvr(i,j) - dt* US%m_to_L**2*G%IareaT(i,j) * & ((ucvr(I,j) - ucvr(I-1,j)) + (vcvr(i,J) - vcvr(i,J-1))) if (cvr(i,j) < 0.0) call SIS_error(FATAL, & 'Negative ice cover encountered in ice_cover_transport().') @@ -262,12 +265,12 @@ subroutine ice_cover_transport(u, v, cvr, dt, G, IG, CS) elseif (x_first) then ! First, advect zonally. LB%ish = G%isc ; LB%ieh = G%iec ; LB%jsh = G%jsc-stensil ; LB%jeh = G%jec+stensil - call zonal_mass_flux(u, dt, G, IG, CS, LB, htot_in=cvr, uh_tot=ucvr) + call zonal_mass_flux(u, dt, G, US, IG, CS, LB, htot_in=cvr, uh_tot=ucvr) call cpu_clock_begin(id_clock_update) !$OMP parallel do default(shared) do j=LB%jsh,LB%jeh ; do i=LB%ish,LB%ieh - cvr(i,j) = cvr(i,j) - G%IareaT(i,j) * (dt*(ucvr(I,j) - ucvr(I-1,j))) + cvr(i,j) = cvr(i,j) - US%m_to_L**2*G%IareaT(i,j) * (dt*(ucvr(I,j) - ucvr(I-1,j))) if (cvr(i,j) < 0.0) call SIS_error(FATAL, & 'Negative ice cover encountered in u-pass of ice_cover_transport().') enddo ; enddo @@ -275,12 +278,12 @@ subroutine ice_cover_transport(u, v, cvr, dt, G, IG, CS) LB%ish = G%isc ; LB%ieh = G%iec ; LB%jsh = G%jsc ; LB%jeh = G%jec ! Now advect meridionally, using the updated ice covers to determine the fluxes. - call meridional_mass_flux(v, dt, G, IG, CS, LB, htot_in=cvr, vh_tot=vcvr) + call meridional_mass_flux(v, dt, G, US, IG, CS, LB, htot_in=cvr, vh_tot=vcvr) call cpu_clock_begin(id_clock_update) !$OMP parallel do default(shared) do j=LB%jsh,LB%jeh ; do i=LB%ish,LB%ieh - cvr(i,j) = max(1.0, cvr(i,j) - dt*G%IareaT(i,j) * (vcvr(i,J) - vcvr(i,J-1))) + cvr(i,j) = max(1.0, cvr(i,j) - dt*US%m_to_L**2*G%IareaT(i,j) * (vcvr(i,J) - vcvr(i,J-1))) if (cvr(i,j) < 0.0) call SIS_error(FATAL, & 'Negative ice cover encountered in v-pass of ice_cover_transport().') enddo ; enddo @@ -289,12 +292,12 @@ subroutine ice_cover_transport(u, v, cvr, dt, G, IG, CS) else ! .not. x_first ! First, advect meridionally, so set the loop bounds accordingly. LB%ish = G%isc-stensil ; LB%ieh = G%iec+stensil ; LB%jsh = G%jsc ; LB%jeh = G%jec - call meridional_mass_flux(v, dt, G, IG, CS, LB, htot_in=cvr, vh_tot=vcvr) + call meridional_mass_flux(v, dt, G, US, IG, CS, LB, htot_in=cvr, vh_tot=vcvr) call cpu_clock_begin(id_clock_update) !$OMP parallel do default(shared) do j=LB%jsh,LB%jeh ; do i=LB%ish,LB%ieh - cvr(i,j) = cvr(i,j) - dt*G%IareaT(i,j) * (vcvr(i,J) - vcvr(i,J-1)) + cvr(i,j) = cvr(i,j) - dt*US%m_to_L**2*G%IareaT(i,j) * (vcvr(i,J) - vcvr(i,J-1)) if (cvr(i,j) < 0.0) call SIS_error(FATAL, & 'Negative ice cover encountered in v-pass of ice_cover_transport().') enddo ; enddo @@ -302,12 +305,12 @@ subroutine ice_cover_transport(u, v, cvr, dt, G, IG, CS) ! Now advect zonally, using the updated ice covers to determine the fluxes. LB%ish = G%isc ; LB%ieh = G%iec ; LB%jsh = G%jsc ; LB%jeh = G%jec - call zonal_mass_flux(u, dt, G, IG, CS, LB, htot_in=cvr, uh_tot=ucvr) + call zonal_mass_flux(u, dt, G, US, IG, CS, LB, htot_in=cvr, uh_tot=ucvr) call cpu_clock_begin(id_clock_update) !$OMP parallel do default(shared) do j=LB%jsh,LB%jeh ; do i=LB%ish,LB%ieh - cvr(i,j) = max(1.0, cvr(i,j) - dt* G%IareaT(i,j) * (ucvr(I,j) - ucvr(I-1,j))) + cvr(i,j) = max(1.0, cvr(i,j) - dt* US%m_to_L**2*G%IareaT(i,j) * (ucvr(I,j) - ucvr(I-1,j))) if (cvr(i,j) < 0.0) call SIS_error(FATAL, & 'Negative ice cover encountered in u-pass of ice_cover_transport().') enddo ; enddo @@ -322,7 +325,7 @@ end subroutine ice_cover_transport !! thickness categories due to advection, using a monotonically limited, directionally split PPM !! scheme or simple upwind 2-d scheme. It may also update the ice thickness, using fluxes that are !! proportional to the total fluxes times the ice mass divided by the total mass in the upwind cell. -subroutine summed_continuity(u, v, h_in, h, uh, vh, dt, G, IG, CS, h_ice) +subroutine summed_continuity(u, v, h_in, h, uh, vh, dt, G, US, IG, CS, h_ice) type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type type(ice_grid_type), intent(inout) :: IG !< The sea-ice specific grid type real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: u !< Zonal ice velocity [m s-1]. @@ -336,6 +339,7 @@ subroutine summed_continuity(u, v, h_in, h, uh, vh, dt, G, IG, CS, h_ice) real, dimension(SZI_(G),SZJB_(G)), intent(out) :: vh !< Total mass flux through meridional faces !! = v*h*dx [H m2 s-1 ~> kg s-1]. real, intent(in) :: dt !< Time increment [s] + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors type(SIS_continuity_CS), pointer :: CS !< The control structure returned by a !! previous call to SIS_continuity_init. real, dimension(SZI_(G),SZJ_(G)), optional, intent(inout) :: h_ice !< Total ice mass per unit cell @@ -376,13 +380,13 @@ subroutine summed_continuity(u, v, h_in, h, uh, vh, dt, G, IG, CS, h_ice) do j=js,je ; do I=is-1,ie if (u(I,j) >= 0.0) then ; h_up = h_in(i,j) else ; h_up = h_in(i+1,j) ; endif - uh(I,j) = G%dy_Cu(I,j) * u(I,j) * h_up + uh(I,j) = US%L_to_m*G%dy_Cu(I,j) * u(I,j) * h_up enddo ; enddo !$OMP do do J=js-1,je ; do i=is,ie if (v(i,J) >= 0.0) then ; h_up = h_in(i,j) else ; h_up = h_in(i,j+1) ; endif - vh(i,J) = G%dx_Cv(i,J) * v(i,J) * h_up + vh(i,J) = US%L_to_m*G%dx_Cv(i,J) * v(i,J) * h_up enddo ; enddo if (present(h_ice)) then !$OMP do @@ -399,13 +403,13 @@ subroutine summed_continuity(u, v, h_in, h, uh, vh, dt, G, IG, CS, h_ice) enddo ; enddo !$OMP do do j=js,je ; do i=is,ie - h_ice(i,j) = h_ice(i,j) - (dt * G%IareaT(i,j)) * & + h_ice(i,j) = h_ice(i,j) - (dt * US%m_to_L**2*G%IareaT(i,j)) * & ((uh_ice(I,j) - uh_ice(I-1,j)) + (vh_ice(i,J) - vh_ice(i,J-1))) enddo ; enddo endif !$OMP do do j=js,je ; do i=is,ie - h(i,j) = h_in(i,j) - (dt * G%IareaT(i,j)) * & + h(i,j) = h_in(i,j) - (dt * US%m_to_L**2*G%IareaT(i,j)) * & ((uh(I,j) - uh(I-1,j)) + (vh(i,J) - vh(i,J-1))) ! if (h(i,j) < 0.0) call SIS_error(FATAL, & ! 'Negative thickness encountered in ice_total_continuity().') @@ -417,7 +421,7 @@ subroutine summed_continuity(u, v, h_in, h, uh, vh, dt, G, IG, CS, h_ice) elseif (x_first) then ! First, advect zonally. LB%ish = G%isc ; LB%ieh = G%iec ; LB%jsh = G%jsc-stensil ; LB%jeh = G%jec+stensil - call zonal_mass_flux(u, dt, G, IG, CS, LB, htot_in=h_in, uh_tot=uh) + call zonal_mass_flux(u, dt, G, US, IG, CS, LB, htot_in=h_in, uh_tot=uh) call cpu_clock_begin(id_clock_update) @@ -430,14 +434,14 @@ subroutine summed_continuity(u, v, h_in, h, uh, vh, dt, G, IG, CS, h_ice) else ; uh_ice(I,j) = 0.0 ; endif enddo do i=LB%ish,LB%ieh - h_ice(i,j) = h_ice(i,j) - (dt * G%IareaT(i,j)) * (uh_ice(I,j) - uh_ice(I-1,j)) + h_ice(i,j) = h_ice(i,j) - (dt * US%m_to_L**2*G%IareaT(i,j)) * (uh_ice(I,j) - uh_ice(I-1,j)) enddo enddo endif !$OMP parallel do default(shared) do j=LB%jsh,LB%jeh ; do i=LB%ish,LB%ieh - h(i,j) = h_in(i,j) - (dt * G%IareaT(i,j)) * (uh(I,j) - uh(I-1,j)) + h(i,j) = h_in(i,j) - (dt * US%m_to_L**2*G%IareaT(i,j)) * (uh(I,j) - uh(I-1,j)) ! if (h(i,j) < 0.0) call SIS_error(FATAL, & ! 'Negative thickness encountered in u-pass of ice_total_continuity().') ! if (present(h_ice)) then ; if (h_ice(i,j) > h(i,j)) then @@ -448,7 +452,7 @@ subroutine summed_continuity(u, v, h_in, h, uh, vh, dt, G, IG, CS, h_ice) LB%ish = G%isc ; LB%ieh = G%iec ; LB%jsh = G%jsc ; LB%jeh = G%jec ! Now advect meridionally, using the updated thicknesses to determine the fluxes. - call meridional_mass_flux(v, dt, G, IG, CS, LB, htot_in=h, vh_tot=vh) + call meridional_mass_flux(v, dt, G, US, IG, CS, LB, htot_in=h, vh_tot=vh) call cpu_clock_begin(id_clock_update) if (present(h_ice)) then @@ -460,13 +464,13 @@ subroutine summed_continuity(u, v, h_in, h, uh, vh, dt, G, IG, CS, h_ice) enddo ; enddo !$OMP parallel do default(shared) do j=LB%jsh,LB%jeh ; do i=LB%ish,LB%ieh - h_ice(i,j) = h_ice(i,j) - (dt * G%IareaT(i,j)) * (vh_ice(i,J) - vh_ice(i,J-1)) + h_ice(i,j) = h_ice(i,j) - (dt * US%m_to_L**2*G%IareaT(i,j)) * (vh_ice(i,J) - vh_ice(i,J-1)) enddo ; enddo endif !$OMP parallel do default(shared) do j=LB%jsh,LB%jeh ; do i=LB%ish,LB%ieh - h(i,j) = h(i,j) - (dt * G%IareaT(i,j)) * (vh(i,J) - vh(i,J-1)) + h(i,j) = h(i,j) - (dt * US%m_to_L**2*G%IareaT(i,j)) * (vh(i,J) - vh(i,J-1)) if (h(i,j) < 0.0) call SIS_error(FATAL, & 'Negative thickness encountered in v-pass of ice_total_continuity().') ! if (present(h_ice)) then ; if (h_ice(i,j) > h(i,j)) then @@ -478,7 +482,7 @@ subroutine summed_continuity(u, v, h_in, h, uh, vh, dt, G, IG, CS, h_ice) else ! .not. x_first ! First, advect meridionally, so set the loop bounds accordingly. LB%ish = G%isc-stensil ; LB%ieh = G%iec+stensil ; LB%jsh = G%jsc ; LB%jeh = G%jec - call meridional_mass_flux(v, dt, G, IG, CS, LB, htot_in=h_in, vh_tot=vh) + call meridional_mass_flux(v, dt, G, US, IG, CS, LB, htot_in=h_in, vh_tot=vh) call cpu_clock_begin(id_clock_update) if (present(h_ice)) then @@ -490,13 +494,13 @@ subroutine summed_continuity(u, v, h_in, h, uh, vh, dt, G, IG, CS, h_ice) enddo ; enddo !$OMP parallel do default(shared) do j=LB%jsh,LB%jeh ; do i=LB%ish,LB%ieh - h_ice(i,j) = h_ice(i,j) - (dt * G%IareaT(i,j)) * (vh_ice(i,J) - vh_ice(i,J-1)) + h_ice(i,j) = h_ice(i,j) - (dt * US%m_to_L**2*G%IareaT(i,j)) * (vh_ice(i,J) - vh_ice(i,J-1)) enddo ; enddo endif !$OMP parallel do default(shared) do j=LB%jsh,LB%jeh ; do i=LB%ish,LB%ieh - h(i,j) = h_in(i,j) - (dt * G%IareaT(i,j)) * (vh(i,J) - vh(i,J-1)) + h(i,j) = h_in(i,j) - (dt * US%m_to_L**2*G%IareaT(i,j)) * (vh(i,J) - vh(i,J-1)) ! if (h(i,j) < 0.0) call SIS_error(FATAL, & ! 'Negative thickness encountered in v-pass of ice_total_continuity().') ! if (present(h_ice)) then ; if (h_ice(i,j) > h(i,j)) then @@ -507,7 +511,7 @@ subroutine summed_continuity(u, v, h_in, h, uh, vh, dt, G, IG, CS, h_ice) ! Now advect zonally, using the updated thicknesses to determine the fluxes. LB%ish = G%isc ; LB%ieh = G%iec ; LB%jsh = G%jsc ; LB%jeh = G%jec - call zonal_mass_flux(u, dt, G, IG, CS, LB, htot_in=h, uh_tot=uh) + call zonal_mass_flux(u, dt, G, US, IG, CS, LB, htot_in=h, uh_tot=uh) call cpu_clock_begin(id_clock_update) @@ -520,14 +524,14 @@ subroutine summed_continuity(u, v, h_in, h, uh, vh, dt, G, IG, CS, h_ice) else ; uh_ice(I,j) = 0.0 ; endif enddo do i=LB%ish,LB%ieh - h_ice(i,j) = h_ice(i,j) - (dt * G%IareaT(i,j)) * (uh_ice(I,j) - uh_ice(I-1,j)) + h_ice(i,j) = h_ice(i,j) - (dt * US%m_to_L**2*G%IareaT(i,j)) * (uh_ice(I,j) - uh_ice(I-1,j)) enddo enddo endif !$OMP parallel do default(shared) do j=LB%jsh,LB%jeh ; do i=LB%ish,LB%ieh - h(i,j) = h(i,j) - (dt * G%IareaT(i,j)) * (uh(I,j) - uh(I-1,j)) + h(i,j) = h(i,j) - (dt * US%m_to_L**2*G%IareaT(i,j)) * (uh(I,j) - uh(I-1,j)) if (h(i,j) < 0.0) call SIS_error(FATAL, & 'Negative thickness encountered in u-pass of ice_continuity().') ! if (present(h_ice)) then ; if (h_ice(i,j) > h(i,j)) then @@ -543,7 +547,7 @@ end subroutine summed_continuity !> proportionate_continuity time steps the category thickness changes due to advection, !! using input total mass fluxes with the fluxes proprotionate to the relative upwind !! thicknesses. -subroutine proportionate_continuity(h_tot_in, uh_tot, vh_tot, dt, G, IG, CS, & +subroutine proportionate_continuity(h_tot_in, uh_tot, vh_tot, dt, G, US, IG, CS, & h1, uh1, vh1, h2, uh2, vh2, h3, uh3, vh3) type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type type(ice_grid_type), intent(inout) :: IG !< The sea-ice specific grid type @@ -554,6 +558,7 @@ subroutine proportionate_continuity(h_tot_in, uh_tot, vh_tot, dt, G, IG, CS, & real, dimension(SZI_(G),SZJB_(G)), intent(in) :: vh_tot !< Total mass flux through meridional faces !! [H m2 s-1 ~> kg s-1]. real, intent(in) :: dt !< Time increment [s] + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors type(SIS_continuity_CS), pointer :: CS !< The control structure returned by a !! previous call to SIS_continuity_init. real, dimension(SZI_(G),SZJ_(G),SZCAT_(IG)), & @@ -617,7 +622,7 @@ subroutine proportionate_continuity(h_tot_in, uh_tot, vh_tot, dt, G, IG, CS, & call merid_proportionate_fluxes(vh_tot, I_htot, h1, vh1, G, IG, LB) !$OMP parallel do default(shared) do j=js,je ; do k=1,nCat ; do i=is,ie - h1(i,j,k) = h1(i,j,k) - G%IareaT(i,j) * (dt * & + h1(i,j,k) = h1(i,j,k) - US%m_to_L**2*G%IareaT(i,j) * (dt * & ((uh1(I,j,k) - uh1(I-1,j,k)) + (vh1(i,J,k) - vh1(i,J-1,k)))) enddo ; enddo ; enddo endif @@ -626,7 +631,7 @@ subroutine proportionate_continuity(h_tot_in, uh_tot, vh_tot, dt, G, IG, CS, & call merid_proportionate_fluxes(vh_tot, I_htot, h2, vh2, G, IG, LB) !$OMP parallel do default(shared) do j=js,je ; do k=1,nCat ; do i=is,ie - h2(i,j,k) = h2(i,j,k) - G%IareaT(i,j) * (dt * & + h2(i,j,k) = h2(i,j,k) - US%m_to_L**2*G%IareaT(i,j) * (dt * & ((uh2(I,j,k) - uh2(I-1,j,k)) + (vh2(i,J,k) - vh2(i,J-1,k)))) enddo ; enddo ; enddo endif @@ -635,7 +640,7 @@ subroutine proportionate_continuity(h_tot_in, uh_tot, vh_tot, dt, G, IG, CS, & call merid_proportionate_fluxes(vh_tot, I_htot, h3, vh3, G, IG, LB) !$OMP parallel do default(shared) do j=js,je ; do k=1,nCat ; do i=is,ie - h3(i,j,k) = h3(i,j,k) - G%IareaT(i,j) * (dt * & + h3(i,j,k) = h3(i,j,k) - US%m_to_L**2*G%IareaT(i,j) * (dt * & ((uh3(I,j,k) - uh3(I-1,j,k)) + (vh3(i,J,k) - vh3(i,J-1,k)))) enddo ; enddo ; enddo endif @@ -647,27 +652,27 @@ subroutine proportionate_continuity(h_tot_in, uh_tot, vh_tot, dt, G, IG, CS, & call zonal_proportionate_fluxes(uh_tot, I_htot, h1, uh1, G, IG, LB) !$OMP parallel do default(shared) do j=LB%jsh,LB%jeh ; do k=1,nCat ; do i=LB%ish,LB%ieh - h1(i,j,k) = h1(i,j,k) - G%IareaT(i,j) * (dt * (uh1(I,j,k) - uh1(I-1,j,k))) + h1(i,j,k) = h1(i,j,k) - US%m_to_L**2*G%IareaT(i,j) * (dt * (uh1(I,j,k) - uh1(I-1,j,k))) enddo ; enddo ; enddo endif if (present(h2)) then call zonal_proportionate_fluxes(uh_tot, I_htot, h2, uh2, G, IG, LB) !$OMP parallel do default(shared) do j=LB%jsh,LB%jeh ; do k=1,nCat ; do i=LB%ish,LB%ieh - h2(i,j,k) = h2(i,j,k) - G%IareaT(i,j) * (dt * (uh2(I,j,k) - uh2(I-1,j,k))) + h2(i,j,k) = h2(i,j,k) - US%m_to_L**2*G%IareaT(i,j) * (dt * (uh2(I,j,k) - uh2(I-1,j,k))) enddo ; enddo ; enddo endif if (present(h3)) then call zonal_proportionate_fluxes(uh_tot, I_htot, h3, uh3, G, IG, LB) !$OMP parallel do default(shared) do j=LB%jsh,LB%jeh ; do k=1,nCat ; do i=LB%ish,LB%ieh - h3(i,j,k) = h3(i,j,k) - G%IareaT(i,j) * (dt * (uh3(I,j,k) - uh3(I-1,j,k))) + h3(i,j,k) = h3(i,j,k) - US%m_to_L**2*G%IareaT(i,j) * (dt * (uh3(I,j,k) - uh3(I-1,j,k))) enddo ; enddo ; enddo endif !$OMP parallel do default(shared) do j=LB%jsh,LB%jeh ; do i=LB%ish,LB%ieh - h_tot(i,j) = h_tot_in(i,j) - dt* G%IareaT(i,j) * (uh_tot(I,j) - uh_tot(I-1,j)) + h_tot(i,j) = h_tot_in(i,j) - dt* US%m_to_L**2*G%IareaT(i,j) * (uh_tot(I,j) - uh_tot(I-1,j)) if (h_tot(i,j) < 0.0) call SIS_error(FATAL, & 'Negative thickness encountered in u-pass of proportionate_continuity().') I_htot(i,j) = 0.0 ; if (h_tot(i,j) > 0.0) I_htot(i,j) = 1.0 / h_tot(i,j) @@ -679,27 +684,27 @@ subroutine proportionate_continuity(h_tot_in, uh_tot, vh_tot, dt, G, IG, CS, & call merid_proportionate_fluxes(vh_tot, I_htot, h1, vh1, G, IG, LB) !$OMP parallel do default(shared) do j=js,je ; do k=1,nCat ; do i=is,ie - h1(i,j,k) = h1(i,j,k) - G%IareaT(i,j) * (dt * (vh1(i,J,k) - vh1(i,J-1,k)) ) + h1(i,j,k) = h1(i,j,k) - US%m_to_L**2*G%IareaT(i,j) * (dt * (vh1(i,J,k) - vh1(i,J-1,k)) ) enddo ; enddo ; enddo endif if (present(h2)) then call merid_proportionate_fluxes(vh_tot, I_htot, h2, vh2, G, IG, LB) !$OMP parallel do default(shared) do j=js,je ; do k=1,nCat ; do i=is,ie - h2(i,j,k) = h2(i,j,k) - G%IareaT(i,j) * (dt * (vh2(i,J,k) - vh2(i,J-1,k)) ) + h2(i,j,k) = h2(i,j,k) - US%m_to_L**2*G%IareaT(i,j) * (dt * (vh2(i,J,k) - vh2(i,J-1,k)) ) enddo ; enddo ; enddo endif if (present(h3)) then call merid_proportionate_fluxes(vh_tot, I_htot, h3, vh3, G, IG, LB) !$OMP parallel do default(shared) do j=js,je ; do k=1,nCat ; do i=is,ie - h3(i,j,k) = h3(i,j,k) - G%IareaT(i,j) * (dt * (vh3(i,J,k) - vh3(i,J-1,k)) ) + h3(i,j,k) = h3(i,j,k) - US%m_to_L**2*G%IareaT(i,j) * (dt * (vh3(i,J,k) - vh3(i,J-1,k)) ) enddo ; enddo ; enddo endif !$OMP parallel do default(shared) do j=LB%jsh,LB%jeh ; do i=LB%ish,LB%ieh - h_tot(i,j) = h_tot(i,j) - dt* G%IareaT(i,j) * (vh_tot(i,J) - vh_tot(i,J-1)) + h_tot(i,j) = h_tot(i,j) - dt* US%m_to_L**2*G%IareaT(i,j) * (vh_tot(i,J) - vh_tot(i,J-1)) if (h_tot(i,j) < 0.0) call SIS_error(FATAL, & 'Negative thickness encountered in v-pass of proportionate_continuity().') ! I_htot(i,j) = 0.0 ; if (h_tot(i,j) > 0.0) I_htot(i,j) = 1.0 / h_tot(i,j) @@ -713,27 +718,27 @@ subroutine proportionate_continuity(h_tot_in, uh_tot, vh_tot, dt, G, IG, CS, & call merid_proportionate_fluxes(vh_tot, I_htot, h1, vh1, G, IG, LB) !$OMP parallel do default(shared) do j=js,je ; do k=1,nCat ; do i=is,ie - h1(i,j,k) = h1(i,j,k) - G%IareaT(i,j) * (dt * (vh1(i,J,k) - vh1(i,J-1,k)) ) + h1(i,j,k) = h1(i,j,k) - US%m_to_L**2*G%IareaT(i,j) * (dt * (vh1(i,J,k) - vh1(i,J-1,k)) ) enddo ; enddo ; enddo endif if (present(h2)) then call merid_proportionate_fluxes(vh_tot, I_htot, h2, vh2, G, IG, LB) !$OMP parallel do default(shared) do j=js,je ; do k=1,nCat ; do i=is,ie - h2(i,j,k) = h2(i,j,k) - G%IareaT(i,j) * (dt * (vh2(i,J,k) - vh2(i,J-1,k)) ) + h2(i,j,k) = h2(i,j,k) - US%m_to_L**2*G%IareaT(i,j) * (dt * (vh2(i,J,k) - vh2(i,J-1,k)) ) enddo ; enddo ; enddo endif if (present(h3)) then call merid_proportionate_fluxes(vh_tot, I_htot, h3, vh3, G, IG, LB) !$OMP parallel do default(shared) do j=js,je ; do k=1,nCat ; do i=is,ie - h3(i,j,k) = h3(i,j,k) - G%IareaT(i,j) * (dt * (vh3(i,J,k) - vh3(i,J-1,k)) ) + h3(i,j,k) = h3(i,j,k) - US%m_to_L**2*G%IareaT(i,j) * (dt * (vh3(i,J,k) - vh3(i,J-1,k)) ) enddo ; enddo ; enddo endif !$OMP parallel do default(shared) do j=LB%jsh,LB%jeh ; do i=LB%ish,LB%ieh - h_tot(i,j) = h_tot(i,j) - dt* G%IareaT(i,j) * (vh_tot(i,J) - vh_tot(i,J-1)) + h_tot(i,j) = h_tot(i,j) - dt* US%m_to_L**2*G%IareaT(i,j) * (vh_tot(i,J) - vh_tot(i,J-1)) if (h_tot(i,j) < 0.0) call SIS_error(FATAL, & 'Negative thickness encountered in v-pass of proportionate_continuity().') I_htot(i,j) = 0.0 ; if (h_tot(i,j) > 0.0) I_htot(i,j) = 1.0 / h_tot(i,j) @@ -746,27 +751,27 @@ subroutine proportionate_continuity(h_tot_in, uh_tot, vh_tot, dt, G, IG, CS, & call zonal_proportionate_fluxes(uh_tot, I_htot, h1, uh1, G, IG, LB) !$OMP parallel do default(shared) do j=LB%jsh,LB%jeh ; do k=1,nCat ; do i=LB%ish,LB%ieh - h1(i,j,k) = h1(i,j,k) - G%IareaT(i,j) * (dt * (uh1(I,j,k) - uh1(I-1,j,k))) + h1(i,j,k) = h1(i,j,k) - US%m_to_L**2*G%IareaT(i,j) * (dt * (uh1(I,j,k) - uh1(I-1,j,k))) enddo ; enddo ; enddo endif if (present(h2)) then call zonal_proportionate_fluxes(uh_tot, I_htot, h2, uh2, G, IG, LB) !$OMP parallel do default(shared) do j=LB%jsh,LB%jeh ; do k=1,nCat ; do i=LB%ish,LB%ieh - h2(i,j,k) = h2(i,j,k) - G%IareaT(i,j) * (dt * (uh2(I,j,k) - uh2(I-1,j,k))) + h2(i,j,k) = h2(i,j,k) - US%m_to_L**2*G%IareaT(i,j) * (dt * (uh2(I,j,k) - uh2(I-1,j,k))) enddo ; enddo ; enddo endif if (present(h3)) then call zonal_proportionate_fluxes(uh_tot, I_htot, h3, uh3, G, IG, LB) !$OMP parallel do default(shared) do j=LB%jsh,LB%jeh ; do k=1,nCat ; do i=LB%ish,LB%ieh - h3(i,j,k) = h3(i,j,k) - G%IareaT(i,j) * (dt * (uh3(I,j,k) - uh3(I-1,j,k))) + h3(i,j,k) = h3(i,j,k) - US%m_to_L**2*G%IareaT(i,j) * (dt * (uh3(I,j,k) - uh3(I-1,j,k))) enddo ; enddo ; enddo endif !$OMP parallel do default(shared) do j=LB%jsh,LB%jeh ; do i=LB%ish,LB%ieh - h_tot(i,j) = h_tot_in(i,j) - dt* G%IareaT(i,j) * (uh_tot(I,j) - uh_tot(I-1,j)) + h_tot(i,j) = h_tot_in(i,j) - dt* US%m_to_L**2*G%IareaT(i,j) * (uh_tot(I,j) - uh_tot(I-1,j)) if (h_tot(i,j) < 0.0) call SIS_error(FATAL, & 'Negative thickness encountered in u-pass of proportionate_continuity().') ! I_htot(i,j) = 0.0 ; if (h_tot(i,j) > 0.0) I_htot(i,j) = 1.0 / h_tot(i,j) @@ -834,12 +839,13 @@ end subroutine merid_proportionate_fluxes !> Calculates the mass or volume fluxes through the zonal !! faces, and other related quantities. -subroutine zonal_mass_flux(u, dt, G, IG, CS, LB, h_in, uh, htot_in, uh_tot) +subroutine zonal_mass_flux(u, dt, G, US, IG, CS, LB, h_in, uh, htot_in, uh_tot) type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type type(ice_grid_type), intent(inout) :: IG !< The sea-ice specific grid type real, dimension(SZIB_(G),SZJ_(G)), & intent(in) :: u !< Zonal ice velocity [m s-1]. real, intent(in) :: dt !< Time increment [s] + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors type(SIS_continuity_CS), pointer :: CS !< The control structure returned by a !! previous call to SIS_continuity_init. type(loop_bounds_type), intent(in) :: LB !< A structure with the active loop bounds. @@ -911,24 +917,24 @@ subroutine zonal_mass_flux(u, dt, G, IG, CS, LB, h_in, uh, htot_in, uh_tot) do I=ish-1,ieh ! Set new values of uh and duhdu. if (u(I,j) > 0.0) then - if (CS%vol_CFL) then ; CFL = (u(I,j) * dt) * (G%dy_Cu(I,j) * G%IareaT(i,j)) - else ; CFL = u(I,j) * dt * G%IdxT(i,j) ; endif + if (CS%vol_CFL) then ; CFL = (u(I,j) * dt) * (US%L_to_m*G%dy_Cu(I,j) * US%m_to_L**2*G%IareaT(i,j)) + else ; CFL = u(I,j) * dt * US%m_to_L*G%IdxT(i,j) ; endif curv_3 = hL(i,j) + hR(i,j) - 2.0*htot(i,j) - uhtot(I) = G%dy_Cu(I,j) * u(I,j) * & + uhtot(I) = US%L_to_m*G%dy_Cu(I,j) * u(I,j) * & (hR(i,j) + CFL * (0.5*(hL(i,j) - hR(i,j)) + curv_3*(CFL - 1.5))) ! h_marg = hR(i,j) + CFL * ((hL(i,j) - hR(i,j)) + 3.0*curv_3*(CFL - 1.0)) elseif (u(I,j) < 0.0) then - if (CS%vol_CFL) then ; CFL = (-u(I,j) * dt) * (G%dy_Cu(I,j) * G%IareaT(i+1,j)) - else ; CFL = -u(I,j) * dt * G%IdxT(i+1,j) ; endif + if (CS%vol_CFL) then ; CFL = (-u(I,j) * dt) * (US%L_to_m*G%dy_Cu(I,j) * US%m_to_L**2*G%IareaT(i+1,j)) + else ; CFL = -u(I,j) * dt * US%m_to_L*G%IdxT(i+1,j) ; endif curv_3 = hL(i+1,j) + hR(i+1,j) - 2.0*htot(i+1,j) - uhtot(I) = G%dy_Cu(I,j) * u(I,j) * & + uhtot(I) = US%L_to_m*G%dy_Cu(I,j) * u(I,j) * & (hL(i+1,j) + CFL * (0.5*(hR(i+1,j)-hL(i+1,j)) + curv_3*(CFL - 1.5))) ! h_marg = hL(i+1) + CFL * ((hR(i+1,j)-hL(i+1,j)) + 3.0*curv_3*(CFL - 1.0)) else uhtot(I) = 0.0 ! h_marg = 0.5 * (hl(i+1,j) + hr(i,j)) endif -! duhdu(I,j) = G%dy_Cu(I,j) * h_marg ! * visc_rem(I) +! duhdu(I,j) = US%L_to_m*G%dy_Cu(I,j) * h_marg ! * visc_rem(I) enddo ! Partition the transports by category in proportion to their relative masses. @@ -956,12 +962,13 @@ end subroutine zonal_mass_flux !> Calculates the mass or volume fluxes through the meridional !! faces, and other related quantities. -subroutine meridional_mass_flux(v, dt, G, IG, CS, LB, h_in, vh, htot_in, vh_tot) +subroutine meridional_mass_flux(v, dt, G, US, IG, CS, LB, h_in, vh, htot_in, vh_tot) type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type type(ice_grid_type), intent(inout) :: IG !< The sea-ice specific grid type real, dimension(SZI_(G),SZJB_(G)), & intent(in) :: v !< Meridional ice velocity [m s-1]. real, intent(in) :: dt !< Time increment [s] + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors type(SIS_continuity_CS), pointer :: CS !< The control structure returned by a !! previous call to SIS_continuity_init. type(loop_bounds_type), intent(in) :: LB !< A structure with the active loop bounds. @@ -1035,24 +1042,24 @@ subroutine meridional_mass_flux(v, dt, G, IG, CS, LB, h_in, vh, htot_in, vh_tot) ! This sets vh and dvhdv. do i=ish,ieh if (v(i,J) > 0.0) then - if (CS%vol_CFL) then ; CFL = (v(i,J) * dt) * (G%dx_Cv(i,J) * G%IareaT(i,j)) - else ; CFL = v(i,J) * dt * G%IdyT(i,j) ; endif + if (CS%vol_CFL) then ; CFL = US%m_to_L*(v(i,J) * dt) * (G%dx_Cv(i,J) * G%IareaT(i,j)) + else ; CFL = v(i,J) * dt * US%m_to_L*G%IdyT(i,j) ; endif curv_3 = hL(i,j) + hR(i,j) - 2.0*htot(i,j) - vhtot(i) = G%dx_Cv(i,J) * v(i,J) * ( hR(i,j) + CFL * & + vhtot(i) = US%L_to_m*G%dx_Cv(i,J) * v(i,J) * ( hR(i,j) + CFL * & (0.5*(hL(i,j) - hR(i,j)) + curv_3*(CFL - 1.5)) ) ! h_marg = hR(i,j) + CFL * ((hL(i,j) - hR(i,j)) + 3.0*curv_3*(CFL - 1.0)) elseif (v(i,J) < 0.0) then - if (CS%vol_CFL) then ; CFL = (-v(i,J) * dt) * (G%dx_Cv(i,J) * G%IareaT(i,j+1)) - else ; CFL = -v(i,J) * dt * G%IdyT(i,j+1) ; endif + if (CS%vol_CFL) then ; CFL = US%m_to_L*(-v(i,J) * dt) * (G%dx_Cv(i,J) * G%IareaT(i,j+1)) + else ; CFL = -v(i,J) * dt * US%m_to_L*G%IdyT(i,j+1) ; endif curv_3 = hL(i,j+1) + hR(i,j+1) - 2.0*htot(i,j+1) - vhtot(i) = G%dx_Cv(i,J) * v(i,J) * ( hL(i,j+1) + CFL * & + vhtot(i) = US%L_to_m*G%dx_Cv(i,J) * v(i,J) * ( hL(i,j+1) + CFL * & (0.5*(hR(i,j+1)-hL(i,j+1)) + curv_3*(CFL - 1.5)) ) ! h_marg = hL(i,j+1) + CFL * ((hR(i,j+1)-hL(i,j+1)) + 3.0*curv_3*(CFL - 1.0)) else vhtot(i) = 0.0 ! h_marg = 0.5 * (hl(i,j+1) + hr(i,j)) endif - ! dvhdv(i) = G%dx_Cv(i,J) * h_marg ! * visc_rem(i) + ! dvhdv(i) = US%L_to_m*G%dx_Cv(i,J) * h_marg ! * visc_rem(i) enddo ! Partition the transports by category in proportion to their relative masses. diff --git a/src/SIS_ctrl_types.F90 b/src/SIS_ctrl_types.F90 index 974efa48..a7772b1c 100644 --- a/src/SIS_ctrl_types.F90 +++ b/src/SIS_ctrl_types.F90 @@ -26,6 +26,7 @@ module SIS_ctrl_types use MOM_file_parser, only : param_file_type use MOM_hor_index, only : hor_index_type use MOM_time_manager, only : time_type, time_type_to_real +use MOM_unit_scaling, only : unit_scale_type use SIS_diag_mediator, only : SIS_diag_ctrl, post_data=>post_SIS_data use SIS_diag_mediator, only : register_SIS_diag_field, register_static_field use SIS_sum_output, only : SIS_sum_out_CS @@ -74,6 +75,7 @@ module SIS_ctrl_types !! shortwave radiation. type(SIS_hor_grid_type), pointer :: G => NULL() !< A structure containing metrics and grid info. + type(unit_scale_type), pointer :: US => NULL() !< A structure containing various unit conversion factors. type(ice_grid_type), pointer :: IG => NULL() !< A structure containing sea-ice specific grid info. type(simple_OSS_type), pointer :: sOSS => NULL() !< A structure containing the arrays !! that describe the ocean's surface state, as it is revealed @@ -148,6 +150,7 @@ module SIS_ctrl_types type(SIS_diag_ctrl) :: diag !< A structure that regulates diagnostics. type(SIS_hor_grid_type), pointer :: G => NULL() !< A structure containing metrics and grid info. + type(unit_scale_type), pointer :: US => NULL() !< A structure containing various unit conversion factors. type(ice_grid_type), pointer :: IG => NULL() !< A structure containing sea-ice specific grid info. type(ocean_sfc_state_type), pointer :: OSS => NULL() !< A structure containing the arrays !! that describe the ocean's surface state, as it is revealed @@ -175,7 +178,7 @@ module SIS_ctrl_types !> ice_diagnostics_init does the registration for a variety of sea-ice model !! diagnostics and saves several static diagnotic fields. -subroutine ice_diagnostics_init(IOF, OSS, FIA, G, IG, diag, Time, Cgrid) +subroutine ice_diagnostics_init(IOF, OSS, FIA, G, US, IG, diag, Time, Cgrid) type(ice_ocean_flux_type), intent(inout) :: IOF !< A structure containing fluxes from the ice to !! the ocean that are calculated by the ice model. type(ocean_sfc_state_type), intent(inout) :: OSS !< A structure containing the arrays that describe @@ -183,6 +186,7 @@ subroutine ice_diagnostics_init(IOF, OSS, FIA, G, IG, diag, Time, Cgrid) type(fast_ice_avg_type), intent(inout) :: FIA !< A type containing averages of fields !! (mostly fluxes) over the fast updates type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors type(ice_grid_type), intent(in) :: IG !< The sea-ice specific grid type type(SIS_diag_ctrl), intent(in) :: diag !< A structure that is used to regulate diagnostic output type(time_type), intent(inout) :: Time !< The sea-ice model's clock, @@ -366,7 +370,7 @@ subroutine ice_diagnostics_init(IOF, OSS, FIA, G, IG, diag, Time, Cgrid) I_area_Earth = 1.0 / (16.0*atan(1.0)*G%Rad_Earth**2) !$OMP parallel do default(none) shared(isc,iec,jsc,jec,G,I_area_Earth,tmp_diag) do j=jsc,jec ; do i=isc,iec - tmp_diag(i,j) = (G%areaT(i,j) * G%mask2dT(i,j)) * I_area_Earth + tmp_diag(i,j) = (US%L_to_m**2*G%areaT(i,j) * G%mask2dT(i,j)) * I_area_Earth enddo ; enddo call post_data(id_cell_area, tmp_diag, diag, is_static=.true.) endif diff --git a/src/SIS_dyn_bgrid.F90 b/src/SIS_dyn_bgrid.F90 index 153a88ad..90b22e3d 100644 --- a/src/SIS_dyn_bgrid.F90 +++ b/src/SIS_dyn_bgrid.F90 @@ -18,6 +18,7 @@ module SIS_dyn_bgrid use MOM_error_handler, only : SIS_error=>MOM_error, FATAL, WARNING, SIS_mesg=>MOM_mesg use MOM_file_parser, only : get_param, log_param, read_param, log_version, param_file_type use MOM_hor_index, only : hor_index_type +use MOM_unit_scaling, only : unit_scale_type use SIS_hor_grid, only : SIS_hor_grid_type use fms_io_mod, only : register_restart_field, restart_file_type use mpp_domains_mod, only : domain2D @@ -245,7 +246,7 @@ end subroutine find_ice_strength !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! !> SIS_B_dynamics takes a single dynamics timestep with EVP subcycles subroutine SIS_B_dynamics(ci, misp, mice, ui, vi, uo, vo, & - fxat, fyat, sea_lev, fxoc, fyoc, do_ridging, rdg_rate, dt_slow, G, CS) + fxat, fyat, sea_lev, fxoc, fyoc, do_ridging, rdg_rate, dt_slow, G, US, CS) type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type real, dimension(SZI_(G),SZJ_(G)), intent(in ) :: ci !< Sea ice concentration [nondim] @@ -267,6 +268,7 @@ subroutine SIS_B_dynamics(ci, misp, mice, ui, vi, uo, vo, & real, dimension(SZIB_(G),SZJB_(G)), intent( out) :: rdg_rate !< ridging rate from drift state in UNITS? real, intent(in ) :: dt_slow !< The amount of time over which the ice !! dynamics are to be advanced [s]. + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors type(SIS_B_dyn_CS), pointer :: CS !< The control structure for this module ! Local variables @@ -340,15 +342,15 @@ subroutine SIS_B_dynamics(ci, misp, mice, ui, vi, uo, vo, & dt_Rheo = dt_slow/EVP_steps do J=jsc-1,jec ; do I=isc-1,iec - dydx(I,J) = 0.5*((G%dyT(i+1,j+1) - G%dyT(i,j+1)) + (G%dyT(i+1,j) - G%dyT(i,j))) - dxdy(I,J) = 0.5*((G%dxT(i+1,j+1) - G%dxT(i+1,j)) + (G%dxT(i,j+1) - G%dxT(i,j))) + dydx(I,J) = 0.5*US%L_to_m*((G%dyT(i+1,j+1) - G%dyT(i,j+1)) + (G%dyT(i+1,j) - G%dyT(i,j))) + dxdy(I,J) = 0.5*US%L_to_m*((G%dxT(i+1,j+1) - G%dxT(i+1,j)) + (G%dxT(i,j+1) - G%dxT(i,j))) enddo ; enddo do j=jsc,jec ; do i=isc,iec - grid_fac1(i,j) = (G%dxCv(i,J)-G%dxCv(i,J-1))*G%IdyT(i,j) - grid_fac2(i,j) = (G%dyCu(I,j)-G%dyCu(I-1,j))*G%IdxT(i,j) - grid_fac3(i,j) = 0.5*G%dyT(i,j) * G%IdxT(i,j) - grid_fac4(i,j) = 0.5*G%dxT(i,j) * G%IdyT(i,j) + grid_fac1(i,j) = US%L_to_m*(G%dxCv(i,J)-G%dxCv(i,J-1))*US%m_to_L*G%IdyT(i,j) + grid_fac2(i,j) = US%L_to_m*(G%dyCu(I,j)-G%dyCu(I-1,j))*US%m_to_L*G%IdxT(i,j) + grid_fac3(i,j) = 0.5*US%L_to_m*G%dyT(i,j) * US%m_to_L*G%IdxT(i,j) + grid_fac4(i,j) = 0.5*US%L_to_m*G%dxT(i,j) * US%m_to_L*G%IdyT(i,j) enddo ; enddo !TOM> check where ice is present @@ -359,9 +361,9 @@ subroutine SIS_B_dynamics(ci, misp, mice, ui, vi, uo, vo, & ! sea level slope force do J=jsc-1,jec ; do I=isc-1,iec sldx(I,J) = -dt_Rheo*G%g_Earth*(0.5*((sea_lev(i+1,j+1)-sea_lev(i,j+1)) & - + (sea_lev(i+1,j)-sea_lev(i,j)))) * G%IdxBu(i,J) + + (sea_lev(i+1,j)-sea_lev(i,j)))) * US%m_to_L*G%IdxBu(i,J) sldy(I,J) = -dt_Rheo*G%g_Earth*(0.5*((sea_lev(i+1,j+1)-sea_lev(i+1,j)) & - + (sea_lev(i,j+1)-sea_lev(i,j)))) * G%IdyBu(I,J) + + (sea_lev(i,j+1)-sea_lev(i,j)))) * US%m_to_L*G%IdyBu(I,J) enddo ; enddo ! put ice/snow mass and concentration on v-grid, first finding mass on t-grid. @@ -386,10 +388,10 @@ subroutine SIS_B_dynamics(ci, misp, mice, ui, vi, uo, vo, & ! This is H&D97, Eq. 44, with their E_0 = 0.25. I_2dt_Rheo = 1.0 / (2.0*dt_Rheo) do j=jsc,jec ; do i=isc,iec - if (G%dxT(i,j) < G%dyT(i,j) ) then - edt(i,j) = I_2dt_Rheo * (G%dxT(i,j)**2 * mice(i,j)) + if (US%L_to_m*G%dxT(i,j) < US%L_to_m*G%dyT(i,j) ) then + edt(i,j) = I_2dt_Rheo * (US%L_to_m**2*G%dxT(i,j)**2 * mice(i,j)) else - edt(i,j) = I_2dt_Rheo * (G%dyT(i,j)**2 * mice(i,j)) + edt(i,j) = I_2dt_Rheo * (US%L_to_m**2*G%dyT(i,j)**2 * mice(i,j)) endif enddo ; enddo endif @@ -425,14 +427,14 @@ subroutine SIS_B_dynamics(ci, misp, mice, ui, vi, uo, vo, & do j=jsc,jec ; do i=isc,iec strn11(i,j) = (0.5*((ui(I,J)-ui(I-1,J)) + (ui(I,J-1)-ui(I-1,J-1))) + & 0.25*((vi(I,J)+vi(I-1,J-1)) + (vi(I,J-1)+vi(I-1,J))) * & - grid_fac1(i,j)) * G%IdxT(i,j) + grid_fac1(i,j)) * US%m_to_L*G%IdxT(i,j) strn22(i,j) = (0.5*((vi(I,J)-vi(I,J-1)) + (vi(I-1,J)-vi(I-1,J-1))) + & 0.25*((ui(I,J)+ui(I,J-1)) + (ui(I-1,J)+ui(I-1,J-1))) * & - grid_fac2(i,j)) * G%IdyT(i,j) - strn12(i,j) = 0.5*grid_fac3(i,j) * & + grid_fac2(i,j)) * US%m_to_L*G%IdyT(i,j) + strn12(i,j) = 0.5*grid_fac3(i,j) * US%m_to_L*& ( (vi(I,J)*G%IdyBu(I,J) - vi(I-1,J)*G%IdyBu(I-1,J)) + & (vi(I,J-1)*G%IdyBu(I,J-1) - vi(I-1,J-1)*G%IdyBu(I-1,J-1)) ) + & - 0.5*grid_fac4(i,j) * & + 0.5*grid_fac4(i,j) * US%m_to_L*& ( (ui(I,J)*G%IdxBu(I,J) - ui(I,J-1)*G%IdxBu(I,J-1)) + & (ui(I-1,J)*G%IdxBu(I-1,J) - ui(I-1,J-1)*G%IdxBu(I-1,J-1)) ) enddo ; enddo @@ -513,20 +515,20 @@ subroutine SIS_B_dynamics(ci, misp, mice, ui, vi, uo, vo, & ! ! first, timestep explicit parts (ice, wind & ocean part of water stress) ! - tmp1 = 0.5*((CS%sig12(i+1,j+1)*G%dxT(i+1,j+1) - CS%sig12(i+1,j)*G%dxT(i+1,j)) + & + tmp1 = 0.5*US%L_to_m*((CS%sig12(i+1,j+1)*G%dxT(i+1,j+1) - CS%sig12(i+1,j)*G%dxT(i+1,j)) + & (CS%sig12(i,j+1)*G%dxT(i,j+1) - CS%sig12(i,j)*G%dxT(i,j)) ) - tmp2 = 0.5*((CS%sig11(i+1,j+1)*G%dyT(i+1,j+1) - CS%sig11(i,j+1)*G%dyT(i,j+1)) + & + tmp2 = 0.5*US%L_to_m*((CS%sig11(i+1,j+1)*G%dyT(i+1,j+1) - CS%sig11(i,j+1)*G%dyT(i,j+1)) + & (CS%sig11(i+1,j)*G%dyT(i+1,j) - CS%sig11(i,j)*G%dyT(i,j)) ) - tmp6 = 0.5*((CS%sig12(i+1,j+1)*G%dyT(i+1,j+1) - CS%sig12(i,j+1)*G%dyT(i,j+1)) + & + tmp6 = 0.5*US%L_to_m*((CS%sig12(i+1,j+1)*G%dyT(i+1,j+1) - CS%sig12(i,j+1)*G%dyT(i,j+1)) + & (CS%sig12(i+1,j)*G%dyT(i+1,j) - CS%sig12(i,j)*G%dyT(i,j)) ) - tmp7 = 0.5*((CS%sig22(i+1,j+1)*G%dxT(i+1,j+1) - CS%sig22(i+1,j)*G%dxT(i+1,j)) + & + tmp7 = 0.5*US%L_to_m*((CS%sig22(i+1,j+1)*G%dxT(i+1,j+1) - CS%sig22(i+1,j)*G%dxT(i+1,j)) + & (CS%sig22(i,j+1)*G%dxT(i,j+1) - CS%sig22(i,j)*G%dxT(i,j))) tmp3 = 0.25*((CS%sig12(i+1,j+1)+CS%sig12(i,j)) + (CS%sig12(i+1,j)+CS%sig12(i,j+1)) ) tmp4 = 0.25*((CS%sig22(i+1,j+1)+CS%sig22(i,j)) + (CS%sig22(i+1,j)+CS%sig22(i,j+1)) ) tmp5 = 0.25*((CS%sig11(i+1,j+1)+CS%sig11(i,j)) + (CS%sig11(i+1,j)+CS%sig11(i,j+1)) ) - fxic_now = ( (tmp1 + tmp2) + (tmp3*dxdy(I,J) - tmp4*dydx(I,J)) ) * G%IareaBu(I,J) - fyic_now = ( (tmp6 + tmp7) + (tmp3*dydx(I,J) - tmp5*dxdy(I,J)) ) * G%IareaBu(I,J) + fxic_now = ( (tmp1 + tmp2) + (tmp3*dxdy(I,J) - tmp4*dydx(I,J)) ) * US%m_to_L**2*G%IareaBu(I,J) + fyic_now = ( (tmp6 + tmp7) + (tmp3*dydx(I,J) - tmp5*dxdy(I,J)) ) * US%m_to_L**2*G%IareaBu(I,J) !### REWRITE TO AVOID COMPLEX EXPRESSIONS. ui(I,J) = ui(I,J) + (fxic_now + civ(I,J)*fxat(I,J) + & diff --git a/src/SIS_dyn_cgrid.F90 b/src/SIS_dyn_cgrid.F90 index e1f96d1e..b93715ad 100644 --- a/src/SIS_dyn_cgrid.F90 +++ b/src/SIS_dyn_cgrid.F90 @@ -27,6 +27,7 @@ module SIS_dyn_cgrid use MOM_io, only : open_file, APPEND_FILE, ASCII_FILE, MULTIPLE, SINGLE_FILE use MOM_time_manager, only : time_type, real_to_time, operator(+), operator(-) use MOM_time_manager, only : set_date, get_time, get_date +use MOM_unit_scaling, only : unit_scale_type use SIS_hor_grid, only : SIS_hor_grid_type use fms_io_mod, only : register_restart_field, restart_file_type use fms_io_mod, only : restore_state, query_initialized @@ -433,7 +434,7 @@ end subroutine find_ice_strength !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! !> SIS_C_dynamics takes a single dynamics timestep with EVP subcycles subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, & - fxat, fyat, sea_lev, fxoc, fyoc, dt_slow, G, CS) + fxat, fyat, sea_lev, fxoc, fyoc, dt_slow, G, US, CS) type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type real, dimension(SZI_(G),SZJ_(G)), intent(in ) :: ci !< Sea ice concentration [nondim] @@ -453,6 +454,7 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, & real, dimension(SZI_(G),SZJB_(G)), intent( out) :: fyoc !< Meridional ice stress on ocean [Pa] real, intent(in ) :: dt_slow !< The amount of time over which the ice !! dynamics are to be advanced [s]. + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors type(SIS_C_dyn_CS), pointer :: CS !< The control structure for this module ! Local variables @@ -514,7 +516,7 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, & ! divided by the sum of the ocean areas around a point [m-2]. q, & ! A potential-vorticity-like field for the ice, the Coriolis parameter ! divided by a spatially averaged mass per unit area [s-1 m2 kg-1]. - dx2B, dy2B, & ! dx^2 or dy^2 at B points [m2]. + dx2B, dy2B, & ! dx^2 or dy^2 at B points [L2 ~> m2]. dx_dyB, dy_dxB ! dx/dy or dy_dx at B points [nondim]. real, dimension(SZIB_(G),SZJ_(G)) :: & azon, bzon, & ! _zon & _mer are the values of the Coriolis force which @@ -528,7 +530,7 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, & real :: fxic_now, fyic_now ! ice internal stress convergence [kg m-1 s-2]. real :: drag_u, drag_v ! Drag rates with the ocean at u & v points [kg m-2 s-1]. real :: drag_max ! A maximum drag rate allowed in the ocean [kg m-2 s-1]. - real :: tot_area ! The sum of the area of the four neighboring cells [m2]. + real :: tot_area ! The sum of the area of the four neighboring cells [L2 ~> m2]. real :: dxharm ! The harmonic mean of the x- and y- grid spacings [m]. real :: muq2, mvq2 ! The product of the u- and v-face masses per unit cell ! area surrounding a vorticity point [kg2 m-4]. @@ -568,7 +570,7 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, & real :: m_neglect ! A tiny mass per unit area [kg m-2]. real :: m_neglect2 ! A tiny mass per unit area squared [kg2 m-4]. real :: m_neglect4 ! A tiny mass per unit area to the 4th power [kg4 m-8]. - real :: sum_area ! The sum of ocean areas around a vorticity point [m2]. + real :: sum_area ! The sum of ocean areas around a vorticity point [L2 ~> m2]. type(time_type) :: & time_it_start, & ! The starting time of the iteratve steps. @@ -661,8 +663,8 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, & !$OMP do do j=jsc,jec do I=isc-1,iec ; if (G%dy_Cu(I,j) > 0.0) then - ui_min_trunc(I,j) = (-CS%CFL_trunc) * G%areaT(i+1,j) / (dt_slow*G%dy_Cu(I,j)) - ui_max_trunc(I,j) = CS%CFL_trunc * G%areaT(i,j) / (dt_slow*G%dy_Cu(I,j)) + ui_min_trunc(I,j) = (-CS%CFL_trunc) * US%L_to_m*G%areaT(i+1,j) / (dt_slow*G%dy_Cu(I,j)) + ui_max_trunc(I,j) = CS%CFL_trunc * US%L_to_m*G%areaT(i,j) / (dt_slow*G%dy_Cu(I,j)) endif ; enddo do I=isc-1,iec ; u_IC(I,j) = ui(I,j) ; enddo enddo @@ -670,8 +672,8 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, & !$OMP do do J=jsc-1,jec do i=isc,iec ; if (G%dx_Cv(i,J) > 0.0) then - vi_min_trunc(i,J) = (-CS%CFL_trunc) * G%areaT(i,j+1) / (dt_slow*G%dx_Cv(i,J)) - vi_max_trunc(i,J) = CS%CFL_trunc * G%areaT(i,j) / (dt_slow*G%dx_Cv(i,J)) + vi_min_trunc(i,J) = (-CS%CFL_trunc) * US%L_to_m*G%areaT(i,j+1) / (dt_slow*G%dx_Cv(i,J)) + vi_max_trunc(i,J) = CS%CFL_trunc * US%L_to_m*G%areaT(i,j) / (dt_slow*G%dx_Cv(i,J)) endif ; enddo do i=isc,iec ; v_IC(i,J) = vi(i,j) ; enddo enddo @@ -684,7 +686,7 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, & ! Precompute pres_mice and the minimum value of del_sh for stability. pres_mice(i,j) = CS%p0_rho*exp(-CS%c0*max(1.0-ci(i,j),0.0)) - dxharm = 2.0*G%dxT(i,j)*G%dyT(i,j) / (G%dxT(i,j) + G%dyT(i,j)) + dxharm = 2.0*US%L_to_m*G%dxT(i,j)*G%dyT(i,j) / (G%dxT(i,j) + G%dyT(i,j)) ! Setting a minimum value of del_sh is sufficient to guarantee numerical ! stability of the overall time-stepping for the velocities and stresses. ! Setting a minimum value of the shear magnitudes is equivalent to setting @@ -700,7 +702,7 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, & ! Ensure that the input stresses are not larger than could be justified by ! the ice pressure now, as the ice might have melted or been advected away ! during the thermodynamic and transport phases. - call limit_stresses(pres_mice, mice, CS%str_d, CS%str_t, CS%str_s, G, CS) + call limit_stresses(pres_mice, mice, CS%str_d, CS%str_t, CS%str_s, G, US, CS) ! Zero out ice velocities with no mass. !$OMP do @@ -732,10 +734,10 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, & !$OMP do do J=jsc-1,jec ; do I=isc-1,iec if (CS%weak_coast_stress) then - sum_area = (G%areaT(i,j) + G%areaT(i+1,j+1)) + (G%areaT(i,j+1) + G%areaT(i+1,j)) + sum_area = US%L_to_m**2*((G%areaT(i,j) + G%areaT(i+1,j+1)) + (G%areaT(i,j+1) + G%areaT(i+1,j))) else - sum_area = (G%mask2dT(i,j)*G%areaT(i,j) + G%mask2dT(i+1,j+1)*G%areaT(i+1,j+1)) + & - (G%mask2dT(i,j+1)*G%areaT(i,j+1) + G%mask2dT(i+1,j)*G%areaT(i+1,j)) + sum_area = US%L_to_m**2*((G%mask2dT(i,j)*G%areaT(i,j) + G%mask2dT(i+1,j+1)*G%areaT(i+1,j+1)) + & + (G%mask2dT(i,j+1)*G%areaT(i,j+1) + G%mask2dT(i+1,j)*G%areaT(i+1,j))) endif if (sum_area <= 0.0) then ! This is a land point. @@ -783,7 +785,7 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, & !$OMP end do nowait !$OMP do do J=jsc-1,jec ; do I=isc-1,iec - tot_area = (G%areaT(i,j) + G%areaT(i+1,j+1)) + (G%areaT(i+1,j) + G%areaT(i,j+1)) + tot_area = ((G%areaT(i,j) + G%areaT(i+1,j+1)) + (G%areaT(i+1,j) + G%areaT(i,j+1))) q(I,J) = G%CoriolisBu(I,J) * tot_area / & (((G%areaT(i,j) * mis(i,j) + G%areaT(i+1,j+1) * mis(i+1,j+1)) + & (G%areaT(i+1,j) * mis(i+1,j) + G%areaT(i,j+1) * mis(i,j+1))) + tot_area * m_neglect) @@ -801,7 +803,7 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, & I1_f2dt2_u(I,j) = 1.0 / ( 1.0 + dt * f2dt_u(I,j) ) ! Calculate the zonal acceleration due to the sea level slope. - PFu(I,j) = -G%g_Earth*(sea_lev(i+1,j)-sea_lev(i,j)) * G%IdxCu(I,j) + PFu(I,j) = -G%g_Earth*(sea_lev(i+1,j)-sea_lev(i,j)) * US%m_to_L*G%IdxCu(I,j) enddo ; enddo !$OMP end do nowait !$OMP do @@ -817,7 +819,7 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, & I1_f2dt2_v(i,J) = 1.0 / ( 1.0 + dt * f2dt_v(i,J) ) ! Calculate the meridional acceleration due to the sea level slope. - PFv(i,J) = -G%g_Earth*(sea_lev(i,j+1)-sea_lev(i,j)) * G%IdyCv(i,J) + PFv(i,J) = -G%g_Earth*(sea_lev(i,j+1)-sea_lev(i,j)) * US%m_to_L*G%IdyCv(i,J) enddo ; enddo !$OMP end parallel @@ -852,21 +854,21 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, & !$OMP dx_dyB,dy_dxB,ui,vi) do J=jsc-halo_sh_Ds,jec+halo_sh_Ds-1 ; do I=isc-halo_sh_Ds,iec+halo_sh_Ds-1 ! This uses a no-slip boundary condition. - sh_Ds(I,J) = (2.0-G%mask2dBu(I,J)) * & + sh_Ds(I,J) = (2.0-G%mask2dBu(I,J)) * US%m_to_L * & (dx_dyB(I,J)*(ui(I,j+1)*G%IdxCu(I,j+1) - ui(I,j)*G%IdxCu(I,j)) + & dy_dxB(I,J)*(vi(i+1,J)*G%IdyCv(i+1,J) - vi(i,J)*G%IdyCv(i,J))) enddo ; enddo if (halo_sh_Ds < 2) call pass_var(sh_Ds, G%Domain, position=CORNER) !$OMP parallel do default(none) shared(isc,iec,jsc,jec,sh_Dt,sh_Dd,dy_dxT,dx_dyT,G,ui,vi) do j=jsc-1,jec+1 ; do i=isc-1,iec+1 - sh_Dt(i,j) = (dy_dxT(i,j)*(G%IdyCu(I,j) * ui(I,j) - & - G%IdyCu(I-1,j)*ui(I-1,j)) - & + sh_Dt(i,j) = US%m_to_L*(dy_dxT(i,j)*(G%IdyCu(I,j) * ui(I,j) - & + G%IdyCu(I-1,j)*ui(I-1,j)) - & dx_dyT(i,j)*(G%IdxCv(i,J) * vi(i,J) - & G%IdxCv(i,J-1)*vi(i,J-1))) - sh_Dd(i,j) = (G%IareaT(i,j)*(G%dyCu(I,j) * ui(I,j) - & - G%dyCu(I-1,j)*ui(I-1,j)) + & - G%IareaT(i,j)*(G%dxCv(i,J) * vi(i,J) - & - G%dxCv(i,J-1)*vi(i,J-1))) + sh_Dd(i,j) = US%m_to_L*(G%IareaT(i,j)*(G%dyCu(I,j) * ui(I,j) - & + G%dyCu(I-1,j)*ui(I-1,j)) + & + G%IareaT(i,j)*(G%dxCv(i,J) * vi(i,J) - & + G%dxCv(i,J-1)*vi(i,J-1))) enddo ; enddo if (CS%project_ci) then @@ -931,8 +933,8 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, & do J=jsc-1,jec ; do I=isc-1,iec ! zeta is already set to 0 over land. CS%str_s(I,J) = I_1pdt_T * ( CS%str_s(I,J) + (I_EC2 * dt_2Tdamp) * & - ( ((G%areaT(i,j)*zeta(i,j) + G%areaT(i+1,j+1)*zeta(i+1,j+1)) + & - (G%areaT(i+1,j)*zeta(i+1,j) + G%areaT(i,j+1)*zeta(i,j+1))) * & + ( US%L_to_m**2*((G%areaT(i,j)*zeta(i,j) + G%areaT(i+1,j+1)*zeta(i+1,j+1)) + & + (G%areaT(i+1,j)*zeta(i+1,j) + G%areaT(i,j+1)*zeta(i,j+1))) * & mi_ratio_A_q(I,J) * sh_Ds(I,J) ) ) enddo ; enddo @@ -958,7 +960,7 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, & ! Evaluate 1/m x.Div(m strain). This expressions include all metric terms ! for an orthogonal grid. The str_d term integrates out to no curl, while ! str_s & str_t terms impose no divergence and do not act on solid body rotation. - fxic_now = G%IdxCu(I,j) * (CS%str_d(i+1,j) - CS%str_d(i,j)) + & + fxic_now = US%m_to_L*G%IdxCu(I,j) * (CS%str_d(i+1,j) - CS%str_d(i,j)) + US%m_to_L * & (G%IdyCu(I,j)*(dy2T(i+1,j)*CS%str_t(i+1,j) - & dy2T(i,j) *CS%str_t(i,j)) + & G%IdxCu(I,j)*(dx2B(I,J) *CS%str_s(I,J) - & @@ -1010,12 +1012,12 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, & ! sum accelerations to take averages. fxic(I,j) = fxic(I,j) + fxic_now - if (CS%id_fix_d>0) fxic_d(I,j) = fxic_d(I,j) + G%mask2dCu(I,j) * & + if (CS%id_fix_d>0) fxic_d(I,j) = fxic_d(I,j) + G%mask2dCu(I,j) * US%m_to_L * & G%IdxCu(I,j) * (CS%str_d(i+1,j) - CS%str_d(i,j)) - if (CS%id_fix_t>0) fxic_t(I,j) = fxic_t(I,j) + G%mask2dCu(I,j) * & + if (CS%id_fix_t>0) fxic_t(I,j) = fxic_t(I,j) + G%mask2dCu(I,j) * US%m_to_L * & G%IdyCu(I,j)*(dy2T(i+1,j)* CS%str_t(i+1,j) - & dy2T(i,j) * CS%str_t(i,j) ) * G%IareaCu(I,j) - if (CS%id_fix_s>0) fxic_s(I,j) = fxic_s(I,j) + G%mask2dCu(I,j) * & + if (CS%id_fix_s>0) fxic_s(I,j) = fxic_s(I,j) + G%mask2dCu(I,j) * US%m_to_L * & G%IdxCu(I,j)*(dx2B(I,J) *CS%str_s(I,J) - & dx2B(I,J-1)*CS%str_s(I,J-1)) * G%IareaCu(I,j) @@ -1040,7 +1042,7 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, & ! Evaluate 1/m y.Div(m strain). This expressions include all metric terms ! for an orthogonal grid. The str_d term integrates out to no curl, while ! str_s & str_t terms impose no divergence and do not act on solid body rotation. - fyic_now = G%IdyCv(i,J) * (CS%str_d(i,j+1)-CS%str_d(i,j)) + & + fyic_now = US%m_to_L*G%IdyCv(i,J) * (CS%str_d(i,j+1)-CS%str_d(i,j)) + US%m_to_L * & (-G%IdxCv(i,J)*(dx2T(i,j+1)*CS%str_t(i,j+1) - & dx2T(i,j) *CS%str_t(i,j)) + & G%IdyCv(i,J)*(dy2B(I,J) *CS%str_s(I,J) - & @@ -1094,11 +1096,11 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, & fyic(i,J) = fyic(i,J) + fyic_now if (CS%id_fiy_d>0) fyic_d(i,J) = fyic_d(i,J) + G%mask2dCv(i,J) * & - G%IdyCv(i,J) * (CS%str_d(i,j+1)-CS%str_d(i,j)) - if (CS%id_fiy_t>0) fyic_t(i,J) = fyic_t(i,J) + G%mask2dCv(i,J) * & + US%m_to_L*G%IdyCv(i,J) * (CS%str_d(i,j+1)-CS%str_d(i,j)) + if (CS%id_fiy_t>0) fyic_t(i,J) = fyic_t(i,J) + G%mask2dCv(i,J) * US%m_to_L * & (G%IdxCv(i,J)*(dx2T(i,j+1)*(-CS%str_t(i,j+1)) - & dx2T(i,j) *(-CS%str_t(i,j))) ) * G%IareaCv(i,J) - if (CS%id_fiy_s>0) fyic_s(i,J) = fyic_s(i,J) + G%mask2dCv(i,J) * & + if (CS%id_fiy_s>0) fyic_s(i,J) = fyic_s(i,J) + G%mask2dCv(i,J) * US%m_to_L * & (G%IdyCv(i,J)*(dy2B(I,J) *CS%str_s(I,J) - & dy2B(I-1,J)*CS%str_s(I-1,J)) ) * G%IareaCv(i,J) @@ -1127,7 +1129,7 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, & call post_SIS_data(CS%id_sigi_hifreq, diag_val, CS%diag) endif if (CS%id_sigii_hifreq>0) then - call find_sigII(mice, ci_proj, CS%str_t, CS%str_s, diag_val, G, CS) + call find_sigII(mice, ci_proj, CS%str_t, CS%str_s, diag_val, G, US, CS) call post_SIS_data(CS%id_sigii_hifreq, diag_val, CS%diag) endif if (CS%id_ci_hifreq>0) call post_SIS_data(CS%id_ci_hifreq, ci_proj, CS%diag) @@ -1197,7 +1199,7 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, & if (mi_u(I,j) > m_neglect) then CS%ntrunc = CS%ntrunc + 1 call write_u_trunc(I, j, ui, u_IC, uo, mis, fxoc, fxic, Cor_u, & - PFu, fxat, dt_slow, G, CS) + PFu, fxat, dt_slow, G, US, CS) endif if (ui(I,j) < ui_min_trunc(I,j)) then ui(I,j) = 0.95 * ui_min_trunc(I,j) @@ -1223,7 +1225,7 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, & if (mi_v(i,J) > m_neglect) then CS%ntrunc = CS%ntrunc + 1 call write_v_trunc(i, J, vi, v_IC, vo, mis, fyoc, fyic, Cor_v, & - PFv, fyat, dt_slow, G, CS) + PFv, fyat, dt_slow, G, US, CS) endif if (vi(i,J) < vi_min_trunc(i,J)) then vi(i,J) = 0.95 * vi_min_trunc(i,J) @@ -1286,7 +1288,7 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, & call post_SIS_data(CS%id_sigi, diag_val, CS%diag) endif if (CS%id_sigii>0) then - call find_sigII(mice, ci, CS%str_t, CS%str_s, diag_val, G, CS) + call find_sigII(mice, ci, CS%str_t, CS%str_s, diag_val, G, US, CS) call post_SIS_data(CS%id_sigii, diag_val, CS%diag) endif if (CS%id_stren>0) then @@ -1348,7 +1350,7 @@ end subroutine SIS_C_dynamics !> limit_stresses ensures that the input stresses are not larger than could be justified by the ice !! pressure now, as the ice might have melted or been advected away during the thermodynamic and !! transport phases, or the ice flow convergence or divergence may have altered the ice concentration. -subroutine limit_stresses(pres_mice, mice, str_d, str_t, str_s, G, CS, limit) +subroutine limit_stresses(pres_mice, mice, str_d, str_t, str_s, G, US, CS, limit) type(SIS_hor_grid_type), intent(in) :: G !< The horizontal grid type real, dimension(SZI_(G),SZJ_(G)), intent(in) :: pres_mice !< The ice internal pressure per !! unit column mass [N m kg-1]. @@ -1357,6 +1359,7 @@ subroutine limit_stresses(pres_mice, mice, str_d, str_t, str_s, G, CS, limit) real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: str_d !< The divergence stress tensor component [Pa m]. real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: str_t !< The tension stress tensor component [Pa m]. real, dimension(SZIB_(G),SZJB_(G)), intent(inout) :: str_s !< The shearing stress tensor component [Pa m]. + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors type(SIS_C_dyn_CS), pointer :: CS !< The control structure for this module real, optional, intent(in) :: limit !< A factor by which the strength limits are changed. @@ -1368,7 +1371,7 @@ subroutine limit_stresses(pres_mice, mice, str_d, str_t, str_s, G, CS, limit) ! Local variables real :: pressure ! The internal ice pressure at a point [Pa]. real :: pres_avg ! The average of the internal ice pressures around a point [Pa]. - real :: sum_area ! The sum of ocean areas around a vorticity point [m2]. + real :: sum_area ! The sum of ocean areas around a vorticity point [L2 ~> m2]. real :: I_2EC ! 1/(2*EC), where EC is the yield curve axis ratio. real :: lim ! A local copy of the factor by which the limits are changed. real :: lim_2 ! The limit divided by 2. @@ -1493,21 +1496,22 @@ end subroutine find_sigI !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! !> find_sigII finds the second stress invariant -subroutine find_sigII(mi, ci, str_t, str_s, sigII, G, CS) +subroutine find_sigII(mi, ci, str_t, str_s, sigII, G, US, CS) type(SIS_hor_grid_type), intent(in) :: G !< The horizontal grid type real, dimension(SZI_(G),SZJ_(G)), intent(in) :: mi !< Mass per unit ocean area of sea ice [kg m-2] real, dimension(SZI_(G),SZJ_(G)), intent(in) :: ci !< Sea ice concentration [nondim] real, dimension(SZI_(G),SZJ_(G)), intent(in) :: str_t !< The tension stress tensor component, [Pa m] real, dimension(SZIB_(G),SZJB_(G)), intent(in) :: str_s !< The shearing stress tensor component [Pa m]. real, dimension(SZI_(G),SZJ_(G)), intent(out) :: sigII !< The second stress invariant [nondim]. + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors type(SIS_C_dyn_CS), pointer :: CS !< The control structure for this module real, dimension(SZI_(G),SZJ_(G)) :: & strength ! The ice strength [Pa m]. real, dimension(SZIB_(G),SZJB_(G)) :: & str_s_ss ! Str_s divided by the sum of the neighboring ice strengths. - real :: strength_sum ! The sum of the 4 neighboring strengths [Pa m]. - real :: sum_area ! The sum of ocean areas around a vorticity point [m2]. + real :: strength_sum ! The sum of the 4 neighboring strengths [L2 Pa m-1 ~> Pa m]. + real :: sum_area ! The sum of ocean areas around a vorticity point [L2 ~> m2]. integer :: i, j, isc, iec, jsc, jec isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec @@ -1647,7 +1651,7 @@ end subroutine SIS_C_dyn_read_alt_restarts !> write_u_trunc is used to record the location of any pseudo-zonal velocity !! truncations and related fields. subroutine write_u_trunc(I, j, ui, u_IC, uo, mis, fxoc, fxic, Cor_u, PFu, fxat, & - dt_slow, G, CS) + dt_slow, G, US, CS) integer, intent(in) :: I !< The i-index of the column to report on integer, intent(in) :: j !< The j-index of the column to report on type(SIS_hor_grid_type), intent(in) :: G !< The horizontal grid type @@ -1661,6 +1665,7 @@ subroutine write_u_trunc(I, j, ui, u_IC, uo, mis, fxoc, fxic, Cor_u, PFu, fxat, real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: PFu !< The zonal Pressure force accleration [m s-2]. real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: fxat !< The zonal wind stress [Pa]. real, intent(in) :: dt_slow !< The slow ice dynamics timestep [s]. + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors type(SIS_C_dyn_CS), pointer :: CS !< The control structure for this module real :: dt_mi, CFL @@ -1685,9 +1690,9 @@ subroutine write_u_trunc(I, j, ui, u_IC, uo, mis, fxoc, fxic, Cor_u, PFu, fxat, file = CS%u_file if (ui(I,j) > 0.0) then - CFL = (ui(I,j) * (dt_slow*G%dy_Cu(I,j))) / G%areaT(i,j) + CFL = (ui(I,j) * (dt_slow*US%m_to_L*G%dy_Cu(I,j))) / G%areaT(i,j) else - CFL = (ui(I,j) * (dt_slow*G%dy_Cu(I,j))) / G%areaT(i+1,j) + CFL = (ui(I,j) * (dt_slow*US%m_to_L*G%dy_Cu(I,j))) / G%areaT(i+1,j) endif @@ -1717,7 +1722,7 @@ end subroutine write_u_trunc !> write_v_trunc is used to record the location of any pseudo-meridional velocity !! truncations and related fields. subroutine write_v_trunc(i, J, vi, v_IC, vo, mis, fyoc, fyic, Cor_v, PFv, fyat, & - dt_slow, G, CS) + dt_slow, G, US, CS) integer, intent(in) :: i !< The i-index of the column to report on integer, intent(in) :: J !< The j-index of the column to report on type(SIS_hor_grid_type), intent(in) :: G !< The horizontal grid type @@ -1731,6 +1736,7 @@ subroutine write_v_trunc(i, J, vi, v_IC, vo, mis, fyoc, fyic, Cor_v, PFv, fyat, real, dimension(SZI_(G),SZJB_(G)), intent(in) :: PFv !< The meridional pressure force accleration [m s-2]. real, dimension(SZI_(G),SZJB_(G)), intent(in) :: fyat !< The meridional wind stress [Pa]. real, intent(in) :: dt_slow !< The slow ice dynamics timestep [s]. + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors type(SIS_C_dyn_CS), pointer :: CS !< The control structure for this module real :: dt_mi, CFL @@ -1756,9 +1762,9 @@ subroutine write_v_trunc(i, J, vi, v_IC, vo, mis, fyoc, fyic, Cor_v, PFv, fyat, file = CS%v_file if (vi(i,J) > 0.0) then - CFL = (vi(i,J) * (dt_slow*G%dx_Cv(i,J))) / G%areaT(i,j) + CFL = (vi(i,J) * (dt_slow*US%m_to_L*G%dx_Cv(i,J))) / G%areaT(i,j) else - CFL = (vi(i,J) * (dt_slow*G%dx_Cv(i,J))) / G%areaT(i,j+1) + CFL = (vi(i,J) * (dt_slow*US%m_to_L*G%dx_Cv(i,J))) / G%areaT(i,j+1) endif call get_date(CS%Time, yr, mo, day, hr, minute, sec) diff --git a/src/SIS_dyn_trans.F90 b/src/SIS_dyn_trans.F90 index a86092d1..ab7a7177 100644 --- a/src/SIS_dyn_trans.F90 +++ b/src/SIS_dyn_trans.F90 @@ -24,6 +24,7 @@ module SIS_dyn_trans use MOM_time_manager, only : time_type, time_type_to_real, real_to_time use MOM_time_manager, only : operator(+), operator(-) use MOM_time_manager, only : operator(>), operator(*), operator(/), operator(/=) +use MOM_unit_scaling, only : unit_scale_type use MOM_EOS, only : EOS_type, calculate_density_derivs use coupler_types_mod, only: coupler_type_initialized, coupler_type_send_data @@ -279,7 +280,7 @@ end subroutine update_icebergs !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! !> SIS_dynamics_trans makes the calls to do ice dynamics and mass and tracer transport -subroutine SIS_dynamics_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, IG, tracer_CSp) +subroutine SIS_dynamics_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, US, IG, tracer_CSp) type(ice_state_type), intent(inout) :: IST !< A type describing the state of the sea ice type(ocean_sfc_state_type), intent(in) :: OSS !< A structure containing the arrays that describe !! the ocean's surface state for the ice model. @@ -289,6 +290,7 @@ subroutine SIS_dynamics_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, I !! the ocean that are calculated by the ice model. real, intent(in) :: dt_slow !< The slow ice dynamics timestep [s]. type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors type(ice_grid_type), intent(inout) :: IG !< The sea-ice specific grid type type(dyn_trans_CS), pointer :: CS !< The control structure for the SIS_dyn_trans module type(icebergs), pointer :: icebergs_CS !< A control structure for the iceberg model. @@ -366,13 +368,13 @@ subroutine SIS_dynamics_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, I if (CS%merged_cont) then ! Convert the category-resolved ice state into the simplified 2-d ice state. ! This should be called after a thermodynamic step or if ice_transport was called. - call convert_IST_to_simple_state(IST, CS%DS2d, CS%CAS, G, IG, CS) + call convert_IST_to_simple_state(IST, CS%DS2d, CS%CAS, G, US, IG, CS) ! Update the category-merged dynamics and use the merged continuity equation. - call SIS_merged_dyn_cont(OSS, FIA, IOF, CS%DS2d, dt_adv_cycle, Time_cycle_start, G, IG, CS) + call SIS_merged_dyn_cont(OSS, FIA, IOF, CS%DS2d, dt_adv_cycle, Time_cycle_start, G, US, IG, CS) ! Complete the category-resolved mass and tracer transport and update the ice state type. - call complete_IST_transport(CS%DS2d, CS%CAS, IST, dt_adv_cycle, G, IG, CS) + call complete_IST_transport(CS%DS2d, CS%CAS, IST, dt_adv_cycle, G, US, IG, CS) else ! (.not.CS%merged_cont) @@ -398,7 +400,7 @@ subroutine SIS_dynamics_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, I enddo ; enddo ! Determine the whole-cell averaged mass of snow and ice. - call ice_state_to_cell_ave_state(IST, G, IG, CS%SIS_transport_CSp, CS%CAS) + call ice_state_to_cell_ave_state(IST, G, US, IG, CS%SIS_transport_CSp, CS%CAS) endif if (.not.CS%Warsaw_sum_order) then do j=jsd,jed ; do i=isd,ied ; ice_free(i,j) = max(1.0 - ice_cover(i,j), 0.0) ; enddo ; enddo @@ -425,7 +427,7 @@ subroutine SIS_dynamics_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, I ! This block of code must be executed if ice_cover and ice_free or the various wind ! stresses were updated. call set_wind_stresses_C(FIA, ice_cover, ice_free, WindStr_x_Cu, WindStr_y_Cv, & - WindStr_x_ocn_Cu, WindStr_y_ocn_Cv, G, CS%complete_ice_cover) + WindStr_x_ocn_Cu, WindStr_y_ocn_Cv, G, US, CS%complete_ice_cover) if (CS%debug) then @@ -446,11 +448,11 @@ subroutine SIS_dynamics_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, I if (CS%Warsaw_sum_order) then call SIS_C_dynamics(1.0-ice_free(:,:), misp_sum, mi_sum, IST%u_ice_C, IST%v_ice_C, & OSS%u_ocn_C, OSS%v_ocn_C, WindStr_x_Cu, WindStr_y_Cv, OSS%sea_lev, & - str_x_ice_ocn_Cu, str_y_ice_ocn_Cv, dt_slow_dyn, G, CS%SIS_C_dyn_CSp) + str_x_ice_ocn_Cu, str_y_ice_ocn_Cv, dt_slow_dyn, G, US, CS%SIS_C_dyn_CSp) else call SIS_C_dynamics(ice_cover, misp_sum, mi_sum, IST%u_ice_C, IST%v_ice_C, & OSS%u_ocn_C, OSS%v_ocn_C, WindStr_x_Cu, WindStr_y_Cv, OSS%sea_lev, & - str_x_ice_ocn_Cu, str_y_ice_ocn_Cv, dt_slow_dyn, G, CS%SIS_C_dyn_CSp) + str_x_ice_ocn_Cu, str_y_ice_ocn_Cv, dt_slow_dyn, G, US, CS%SIS_C_dyn_CSp) endif call mpp_clock_end(iceClocka) @@ -486,7 +488,7 @@ subroutine SIS_dynamics_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, I ! This block of code must be executed if ice_cover and ice_free or the various wind ! stresses were updated. call set_wind_stresses_B(FIA, ice_cover, ice_free, WindStr_x_B, WindStr_y_B, & - WindStr_x_ocn_B, WindStr_y_ocn_B, G, CS%complete_ice_cover) + WindStr_x_ocn_B, WindStr_y_ocn_B, G, US, CS%complete_ice_cover) if (CS%debug) then call Bchksum_pair("[uv]_ice_B before dynamics", IST%u_ice_B, IST%v_ice_B, G) @@ -504,12 +506,12 @@ subroutine SIS_dynamics_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, I call SIS_B_dynamics(1.0-ice_free(:,:), misp_sum, mi_sum, IST%u_ice_B, IST%v_ice_B, & OSS%u_ocn_B, OSS%v_ocn_B, WindStr_x_B, WindStr_y_B, OSS%sea_lev, & str_x_ice_ocn_B, str_y_ice_ocn_B, CS%do_ridging, & - rdg_rate(isc:iec,jsc:jec), dt_slow_dyn, G, CS%SIS_B_dyn_CSp) + rdg_rate(isc:iec,jsc:jec), dt_slow_dyn, G, US, CS%SIS_B_dyn_CSp) else call SIS_B_dynamics(ice_cover, misp_sum, mi_sum, IST%u_ice_B, IST%v_ice_B, & OSS%u_ocn_B, OSS%v_ocn_B, WindStr_x_B, WindStr_y_B, OSS%sea_lev, & str_x_ice_ocn_B, str_y_ice_ocn_B, CS%do_ridging, & - rdg_rate(isc:iec,jsc:jec), dt_slow_dyn, G, CS%SIS_B_dyn_CSp) + rdg_rate(isc:iec,jsc:jec), dt_slow_dyn, G, US, CS%SIS_B_dyn_CSp) endif call mpp_clock_end(iceClocka) @@ -573,16 +575,16 @@ subroutine SIS_dynamics_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, I if (CS%debug) call uvchksum("Before ice_transport [uv]_ice_C", IST%u_ice_C, IST%v_ice_C, G) call enable_SIS_averaging(dt_slow_dyn, Time_cycle_start + real_to_time(nds*dt_slow_dyn), CS%diag) - call ice_cat_transport(CS%CAS, IST%TrReg, dt_slow_dyn, CS%adv_substeps, G, IG, CS%SIS_transport_CSp, & + call ice_cat_transport(CS%CAS, IST%TrReg, dt_slow_dyn, CS%adv_substeps, G, US, IG, CS%SIS_transport_CSp, & uc=IST%u_ice_C, vc=IST%v_ice_C) if (DS2d%nts==0) then if (CS%do_ridging) then - call finish_ice_transport(CS%CAS, IST, IST%TrReg, G, IG, CS%SIS_transport_CSp, & + call finish_ice_transport(CS%CAS, IST, IST%TrReg, G, US, IG, CS%SIS_transport_CSp, & rdg_rate=DS2d%avg_ridge_rate) DS2d%ridge_rate_count = 0. ; DS2d%avg_ridge_rate(:,:) = 0.0 else - call finish_ice_transport(CS%CAS, IST, IST%TrReg, G, IG, CS%SIS_transport_CSp) + call finish_ice_transport(CS%CAS, IST, IST%TrReg, G, US, IG, CS%SIS_transport_CSp) endif endif call mpp_clock_end(iceClock8) @@ -590,7 +592,7 @@ subroutine SIS_dynamics_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, I endif ! (.not.CS%merged_cont) if (CS%column_check .and. (DS2d%nts==0)) & - call write_ice_statistics(IST, CS%Time, CS%n_calls, G, IG, CS%sum_output_CSp, & + call write_ice_statistics(IST, CS%Time, CS%n_calls, G, US, IG, CS%sum_output_CSp, & message=" Post_transport")! , check_column=.true.) enddo ! nac = 1,nadv_cycle @@ -599,7 +601,7 @@ subroutine SIS_dynamics_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, I call finish_ocean_top_stresses(IOF, G) ! Do diagnostics and update some information for the atmosphere. - call ice_state_cleanup(IST, OSS, IOF, dt_slow, G, IG, CS, tracer_CSp) + call ice_state_cleanup(IST, OSS, IOF, dt_slow, G, US, IG, CS, tracer_CSp) end subroutine SIS_dynamics_trans @@ -607,7 +609,7 @@ end subroutine SIS_dynamics_trans !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! !> SIS_multi_dyn_trans makes the calls to do ice dynamics and mass and tracer transport as !! appropriate for a dynamic and advective update cycle with multiple calls. -subroutine SIS_multi_dyn_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, IG, tracer_CSp, & +subroutine SIS_multi_dyn_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, US, IG, tracer_CSp, & start_cycle, end_cycle, cycle_length) type(ice_state_type), intent(inout) :: IST !< A type describing the state of the sea ice type(ocean_sfc_state_type), intent(in) :: OSS !< A structure containing the arrays that describe @@ -618,6 +620,7 @@ subroutine SIS_multi_dyn_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, !! the ocean that are calculated by the ice model. real, intent(in) :: dt_slow !< The slow ice dynamics timestep [s]. type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors type(ice_grid_type), intent(inout) :: IG !< The sea-ice specific grid type type(dyn_trans_CS), pointer :: CS !< The control structure for the SIS_dyn_trans module type(icebergs), pointer :: icebergs_CS !< A control structure for the iceberg model. @@ -658,22 +661,22 @@ subroutine SIS_multi_dyn_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, ! Convert the category-resolved ice state into the simplified 2-d ice state. ! This should be called after a thermodynamic step or if ice_transport was called. if ((nac > 1) .or. cycle_start) & - call convert_IST_to_simple_state(IST, CS%DS2d, CS%CAS, G, IG, CS) + call convert_IST_to_simple_state(IST, CS%DS2d, CS%CAS, G, US, IG, CS) ! Update the category-merged dynamics and use the merged continuity equation. ! This could be called as many times as necessary. Time_cycle_start = CS%Time - real_to_time((nadv_cycle-(nac-1))*dt_adv_cycle) end_of_cycle = (nac < nadv_cycle) .or. cycle_end - call SIS_merged_dyn_cont(OSS, FIA, IOF, CS%DS2d, dt_adv_cycle, Time_cycle_start, G, IG, CS, & + call SIS_merged_dyn_cont(OSS, FIA, IOF, CS%DS2d, dt_adv_cycle, Time_cycle_start, G, US, IG, CS, & end_call=end_of_cycle) ! Complete the category-resolved mass and tracer transport and update the ice state type. ! This must be done before the next thermodynamic step. if (end_of_cycle) & - call complete_IST_transport(CS%DS2d, CS%CAS, IST, dt_adv_cycle, G, IG, CS) + call complete_IST_transport(CS%DS2d, CS%CAS, IST, dt_adv_cycle, G, US, IG, CS) if (CS%column_check .and. IST%valid_IST) & ! This is just here from early debugging exercises, - call write_ice_statistics(IST, CS%Time, CS%n_calls, G, IG, CS%sum_output_CSp, & + call write_ice_statistics(IST, CS%Time, CS%n_calls, G, US, IG, CS%sum_output_CSp, & message=" Post_transport")! , check_column=.true.) enddo ! nac=0,nadv_cycle-1 @@ -683,19 +686,20 @@ subroutine SIS_multi_dyn_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, ! This must be done before returning control to the atmosphere and before writing any diagnostics. if (cycle_end) & - call ice_state_cleanup(IST, OSS, IOF, dt_diags, G, IG, CS, tracer_CSp) + call ice_state_cleanup(IST, OSS, IOF, dt_diags, G, US, IG, CS, tracer_CSp) end subroutine SIS_multi_dyn_trans !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! !> Complete the category-resolved mass and tracer transport and update the ice state type. -subroutine complete_IST_transport(DS2d, CAS, IST, dt_adv_cycle, G, IG, CS) +subroutine complete_IST_transport(DS2d, CAS, IST, dt_adv_cycle, G, US, IG, CS) type(ice_state_type), intent(inout) :: IST !< A type describing the state of the sea ice type(dyn_state_2d), intent(inout) :: DS2d !< A simplified 2-d description of the ice state !! integrated across thickness categories and layers. type(cell_average_state_type), intent(inout) :: CAS !< A structure with ocean-cell averaged masses. real, intent(in) :: dt_adv_cycle !< The time since the last IST transport [s]. type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors type(ice_grid_type), intent(inout) :: IG !< The sea-ice specific grid type type(dyn_trans_CS), pointer :: CS !< The control structure for the SIS_dyn_trans module @@ -708,17 +712,17 @@ subroutine complete_IST_transport(DS2d, CAS, IST, dt_adv_cycle, G, IG, CS) call mpp_clock_begin(iceClock8) ! Do the transport of mass and tracers by category and vertical layer. - call ice_cat_transport(CS%CAS, IST%TrReg, dt_adv_cycle, DS2d%nts, G, IG, & + call ice_cat_transport(CS%CAS, IST%TrReg, dt_adv_cycle, DS2d%nts, G, US, IG, & CS%SIS_transport_CSp, mca_tot=DS2d%mca_step(:,:,0:DS2d%nts), & uh_tot=DS2d%uh_step(:,:,1:DS2d%nts), vh_tot=DS2d%vh_step(:,:,1:DS2d%nts)) ! Convert the cell-averaged state back to the ice-state type, adjusting the ! category mass distributions, doing ridging, and updating the partition sizes. if (CS%do_ridging) then - call finish_ice_transport(CS%CAS, IST, IST%TrReg, G, IG, CS%SIS_transport_CSp, & + call finish_ice_transport(CS%CAS, IST, IST%TrReg, G, US, IG, CS%SIS_transport_CSp, & rdg_rate=DS2d%avg_ridge_rate) DS2d%ridge_rate_count = 0. ; DS2d%avg_ridge_rate(:,:) = 0.0 else - call finish_ice_transport(CS%CAS, IST, IST%TrReg, G, IG, CS%SIS_transport_CSp) + call finish_ice_transport(CS%CAS, IST, IST%TrReg, G, US, IG, CS%SIS_transport_CSp) endif DS2d%nts = 0 ! There is no outstanding transport to be done and IST is up-to-date. @@ -739,7 +743,7 @@ end subroutine complete_IST_transport !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! !> Do final checks to set a consistent ice state and write diagnostics as appropriate. -subroutine ice_state_cleanup(IST, OSS, IOF, dt_slow, G, IG, CS, tracer_CSp) +subroutine ice_state_cleanup(IST, OSS, IOF, dt_slow, G, US, IG, CS, tracer_CSp) type(ice_state_type), intent(inout) :: IST !< A type describing the state of the sea ice type(ocean_sfc_state_type), intent(in) :: OSS !< A structure containing the arrays that describe !! the ocean's surface state for the ice model. @@ -747,6 +751,7 @@ subroutine ice_state_cleanup(IST, OSS, IOF, dt_slow, G, IG, CS, tracer_CSp) !! the ocean that are calculated by the ice model. real, intent(in) :: dt_slow !< The slow ice dynamics timestep [s]. type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors type(ice_grid_type), intent(inout) :: IG !< The sea-ice specific grid type type(dyn_trans_CS), pointer :: CS !< The control structure for the SIS_dyn_trans module type(SIS_tracer_flow_control_CS), optional, pointer :: tracer_CSp !< The structure for controlling @@ -777,11 +782,11 @@ subroutine ice_state_cleanup(IST, OSS, IOF, dt_slow, G, IG, CS, tracer_CSp) if (CS%bounds_check) call IST_bounds_check(IST, G, IG, "End of ice_state_cleanup", OSS=OSS) if (CS%Time + real_to_time(0.5*dt_slow) > CS%write_ice_stats_time) then - call write_ice_statistics(IST, CS%Time, CS%n_calls, G, IG, CS%sum_output_CSp, & + call write_ice_statistics(IST, CS%Time, CS%n_calls, G, US, IG, CS%sum_output_CSp, & tracer_CSp=tracer_CSp) CS%write_ice_stats_time = CS%write_ice_stats_time + CS%ice_stats_interval elseif (CS%column_check) then - call write_ice_statistics(IST, CS%Time, CS%n_calls, G, IG, CS%sum_output_CSp) + call write_ice_statistics(IST, CS%Time, CS%n_calls, G, US, IG, CS%sum_output_CSp) endif call mpp_clock_end(iceClock9) @@ -790,12 +795,13 @@ end subroutine ice_state_cleanup !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! !> Convert the category-resolved ice state into the simplified 2-d ice state and a cell averaged state. -subroutine convert_IST_to_simple_state(IST, DS2d, CAS, G, IG, CS) +subroutine convert_IST_to_simple_state(IST, DS2d, CAS, G, US, IG, CS) type(ice_state_type), intent(inout) :: IST !< A type describing the state of the sea ice type(dyn_state_2d), intent(inout) :: DS2d !< A simplified 2-d description of the ice state !! integrated across thickness categories and layers. type(cell_average_state_type), intent(inout) :: CAS !< A structure with ocean-cell averaged masses. type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors type(ice_grid_type), intent(inout) :: IG !< The sea-ice specific grid type type(dyn_trans_CS), pointer :: CS !< The control structure for the SIS_dyn_trans module @@ -838,7 +844,7 @@ subroutine convert_IST_to_simple_state(IST, DS2d, CAS, G, IG, CS) endif ! Determine the whole-cell averaged mass of snow and ice. - call ice_state_to_cell_ave_state(IST, G, IG, CS%SIS_transport_CSp, CAS) + call ice_state_to_cell_ave_state(IST, G, US, IG, CS%SIS_transport_CSp, CAS) IST%valid_IST = .false. @@ -849,7 +855,7 @@ end subroutine convert_IST_to_simple_state !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! !> Update the category-merged ice state and call the merged continuity update. -subroutine SIS_merged_dyn_cont(OSS, FIA, IOF, DS2d, dt_cycle, Time_start, G, IG, CS, end_call) +subroutine SIS_merged_dyn_cont(OSS, FIA, IOF, DS2d, dt_cycle, Time_start, G, US, IG, CS, end_call) type(ocean_sfc_state_type), intent(in) :: OSS !< A structure containing the arrays that describe !! the ocean's surface state for the ice model. type(fast_ice_avg_type), intent(in) :: FIA !< A type containing averages of fields @@ -861,6 +867,7 @@ subroutine SIS_merged_dyn_cont(OSS, FIA, IOF, DS2d, dt_cycle, Time_start, G, IG, real, intent(in) :: dt_cycle !< The slow ice dynamics timestep [s]. type(time_type), intent(in) :: TIme_start !< The starting time for this update cycle. type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors type(ice_grid_type), intent(inout) :: IG !< The sea-ice specific grid type type(dyn_trans_CS), pointer :: CS !< The control structure for the SIS_dyn_trans module logical, optional, intent(in) :: end_call !< If present and false, this call is @@ -932,7 +939,7 @@ subroutine SIS_merged_dyn_cont(OSS, FIA, IOF, DS2d, dt_cycle, Time_start, G, IG, ! This block of code must be executed if ice_cover and ice_free or the various wind ! stresses were updated. call set_wind_stresses_C(FIA, DS2d%ice_cover, ice_free, WindStr_x_Cu, WindStr_y_Cv, & - WindStr_x_ocn_Cu, WindStr_y_ocn_Cv, G, CS%complete_ice_cover) + WindStr_x_ocn_Cu, WindStr_y_ocn_Cv, G, US, CS%complete_ice_cover) if (CS%debug) then call uvchksum("Before SIS_C_dynamics [uv]_ice_C", DS2d%u_ice_C, DS2d%v_ice_C, G) @@ -951,7 +958,7 @@ subroutine SIS_merged_dyn_cont(OSS, FIA, IOF, DS2d, dt_cycle, Time_start, G, IG, if (CS%do_ridging) rdg_rate(:,:) = 0.0 call SIS_C_dynamics(DS2d%ice_cover, DS2d%mca_step(:,:,DS2d%nts), DS2d%mi_sum, DS2d%u_ice_C, DS2d%v_ice_C, & OSS%u_ocn_C, OSS%v_ocn_C, WindStr_x_Cu, WindStr_y_Cv, OSS%sea_lev, & - str_x_ice_ocn_Cu, str_y_ice_ocn_Cv, dt_slow_dyn, G, CS%SIS_C_dyn_CSp) + str_x_ice_ocn_Cu, str_y_ice_ocn_Cv, dt_slow_dyn, G, US, CS%SIS_C_dyn_CSp) call mpp_clock_end(iceClocka) if (CS%debug) call uvchksum("After ice_dynamics [uv]_ice_C", DS2d%u_ice_C, DS2d%v_ice_C, G) @@ -980,7 +987,7 @@ subroutine SIS_merged_dyn_cont(OSS, FIA, IOF, DS2d, dt_cycle, Time_start, G, IG, ! stresses were updated. call set_wind_stresses_B(FIA, DS2d%ice_cover, ice_free, WindStr_x_B, WindStr_y_B, & - WindStr_x_ocn_B, WindStr_y_ocn_B, G, CS%complete_ice_cover) + WindStr_x_ocn_B, WindStr_y_ocn_B, G, US, CS%complete_ice_cover) if (CS%debug) then call Bchksum_pair("[uv]_ice_B before dynamics", DS2d%u_ice_B, DS2d%v_ice_B, G) @@ -997,7 +1004,7 @@ subroutine SIS_merged_dyn_cont(OSS, FIA, IOF, DS2d, dt_cycle, Time_start, G, IG, call SIS_B_dynamics(DS2d%ice_cover, DS2d%mca_step(:,:,DS2d%nts), DS2d%mi_sum, DS2d%u_ice_B, DS2d%v_ice_B, & OSS%u_ocn_B, OSS%v_ocn_B, WindStr_x_B, WindStr_y_B, OSS%sea_lev, & str_x_ice_ocn_B, str_y_ice_ocn_B, CS%do_ridging, & - rdg_rate(isc:iec,jsc:jec), dt_slow_dyn, G, CS%SIS_B_dyn_CSp) + rdg_rate(isc:iec,jsc:jec), dt_slow_dyn, G, US, CS%SIS_B_dyn_CSp) call mpp_clock_end(iceClocka) if (CS%debug) call Bchksum_pair("After dynamics [uv]_ice_B", DS2d%u_ice_B, DS2d%v_ice_B, G) @@ -1056,15 +1063,15 @@ subroutine SIS_merged_dyn_cont(OSS, FIA, IOF, DS2d, dt_cycle, Time_start, G, IG, if ((n < ndyn_steps*CS%adv_substeps) .or. continuing_call) then ! Some of the work is not needed for the last step before cat_ice_transport. call summed_continuity(DS2d%u_ice_C, DS2d%v_ice_C, DS2d%mca_step(:,:,n-1), DS2d%mca_step(:,:,n), & - DS2d%uh_step(:,:,n), DS2d%vh_step(:,:,n), dt_adv, G, IG, CS%continuity_CSp, & + DS2d%uh_step(:,:,n), DS2d%vh_step(:,:,n), dt_adv, G, US, IG, CS%continuity_CSp, & h_ice=DS2d%mi_sum) - call ice_cover_transport(DS2d%u_ice_C, DS2d%v_ice_C, DS2d%ice_cover, dt_adv, G, IG, CS%cover_trans_CSp) + call ice_cover_transport(DS2d%u_ice_C, DS2d%v_ice_C, DS2d%ice_cover, dt_adv, G, US, IG, CS%cover_trans_CSp) call pass_var(DS2d%mi_sum, G%Domain, complete=.false.) call pass_var(DS2d%ice_cover, G%Domain, complete=.false.) call pass_var(DS2d%mca_step(:,:,n), G%Domain, complete=.true.) else call summed_continuity(DS2d%u_ice_C, DS2d%v_ice_C, DS2d%mca_step(:,:,n-1), DS2d%mca_step(:,:,n), & - DS2d%uh_step(:,:,n), DS2d%vh_step(:,:,n), dt_adv, G, IG, CS%continuity_CSp) + DS2d%uh_step(:,:,n), DS2d%vh_step(:,:,n), dt_adv, G, US, IG, CS%continuity_CSp) endif enddo DS2d%nts = DS2d%nts + CS%adv_substeps @@ -1076,7 +1083,7 @@ end subroutine SIS_merged_dyn_cont !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! !> slab_ice_dynamics_trans makes the calls to do the slab ice version of dynamics and mass and tracer transport -subroutine slab_ice_dyn_trans(IST, OSS, FIA, IOF, dt_slow, CS, G, IG, tracer_CSp) +subroutine slab_ice_dyn_trans(IST, OSS, FIA, IOF, dt_slow, CS, G, US, IG, tracer_CSp) type(ice_state_type), intent(inout) :: IST !< A type describing the state of the sea ice type(ocean_sfc_state_type), intent(in) :: OSS !< A structure containing the arrays that describe @@ -1087,6 +1094,7 @@ subroutine slab_ice_dyn_trans(IST, OSS, FIA, IOF, dt_slow, CS, G, IG, tracer_CSp !! the ocean that are calculated by the ice model. real, intent(in) :: dt_slow !< The slow ice dynamics timestep [s]. type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors type(ice_grid_type), intent(inout) :: IG !< The sea-ice specific grid type type(dyn_trans_CS), pointer :: CS !< The control structure for the SIS_dyn_trans module type(SIS_tracer_flow_control_CS), pointer :: tracer_CSp !< The structure for controlling calls to @@ -1160,7 +1168,7 @@ subroutine slab_ice_dyn_trans(IST, OSS, FIA, IOF, dt_slow, CS, G, IG, tracer_CSp ! This block of code must be executed if ice_cover and ice_free or the various wind ! stresses were updated. call set_wind_stresses_C(FIA, IST%part_size(:,:,1), IST%part_size(:,:,0), WindStr_x_Cu, WindStr_y_Cv, & - WindStr_x_ocn_Cu, WindStr_y_ocn_Cv, G, CS%complete_ice_cover) + WindStr_x_ocn_Cu, WindStr_y_ocn_Cv, G, US, CS%complete_ice_cover) if (CS%debug) then call uvchksum("Before SIS_C_dynamics [uv]_ice_C", IST%u_ice_C, IST%v_ice_C, G) @@ -1209,7 +1217,7 @@ subroutine slab_ice_dyn_trans(IST, OSS, FIA, IOF, dt_slow, CS, G, IG, tracer_CSp ! stresses were updated. call set_wind_stresses_B(FIA, IST%part_size(:,:,1), IST%part_size(:,:,0), WindStr_x_B, WindStr_y_B, & - WindStr_x_ocn_B, WindStr_y_ocn_B, G, CS%complete_ice_cover) + WindStr_x_ocn_B, WindStr_y_ocn_B, G, US, CS%complete_ice_cover) if (CS%debug) then call Bchksum_pair("[uv]_ice_B before dynamics", IST%u_ice_B, IST%v_ice_B, G) @@ -1273,17 +1281,17 @@ subroutine slab_ice_dyn_trans(IST, OSS, FIA, IOF, dt_slow, CS, G, IG, tracer_CSp call enable_SIS_averaging(dt_slow_dyn, CS%Time - real_to_time((ndyn_steps-nds)*dt_slow_dyn), CS%diag) call slab_ice_advect(IST%u_ice_C, IST%v_ice_C, IST%mH_ice(:,:,1), 4.0*IG%kg_m2_to_H, & - dt_slow_dyn, G, IST%part_size(:,:,1), nsteps=CS%adv_substeps) + dt_slow_dyn, G, US, IST%part_size(:,:,1), nsteps=CS%adv_substeps) call mpp_clock_end(iceClock8) if (CS%column_check) & - call write_ice_statistics(IST, CS%Time, CS%n_calls, G, IG, CS%sum_output_CSp, & + call write_ice_statistics(IST, CS%Time, CS%n_calls, G, US, IG, CS%sum_output_CSp, & message=" Post_transport")! , check_column=.true.) enddo ! nds=1,ndyn_steps call finish_ocean_top_stresses(IOF, G) - call ice_state_cleanup(IST, OSS, IOF, dt_slow, G, IG, CS, tracer_CSp) + call ice_state_cleanup(IST, OSS, IOF, dt_slow, G, US, IG, CS, tracer_CSp) end subroutine slab_ice_dyn_trans @@ -1804,7 +1812,7 @@ end subroutine set_ocean_top_stress_C2 !> set_wind_stresses_C determines the wind stresses on the ice and open ocean with !! a C-grid staggering of the points. subroutine set_wind_stresses_C(FIA, ice_cover, ice_free, WindStr_x_Cu, WindStr_y_Cv, & - WindStr_x_ocn_Cu, WindStr_y_ocn_Cv, G, max_ice_cover) + WindStr_x_ocn_Cu, WindStr_y_ocn_Cv, G, US, max_ice_cover) type(fast_ice_avg_type), intent(in) :: FIA !< A type containing averages of fields !! (mostly fluxes) over the fast updates type(SIS_hor_grid_type), intent(in) :: G !< The horizontal grid type @@ -1818,6 +1826,7 @@ subroutine set_wind_stresses_C(FIA, ice_cover, ice_free, WindStr_x_Cu, WindStr_y real, dimension(SZI_(G),SZJB_(G)), intent(out) :: & WindStr_y_Cv, & !< Meridional wind stress averaged over the ice categores on C-grid v-points [Pa]. WindStr_y_ocn_Cv !< Meridional wind stress on the ice-free ocean on C-grid v-points [Pa]. + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors real, intent(in) :: max_ice_cover !< The fractional ice coverage !! that is close enough to 1 to be complete for the purpose of calculating !! wind stresses [nondim]. @@ -1872,18 +1881,18 @@ subroutine set_wind_stresses_C(FIA, ice_cover, ice_free, WindStr_x_Cu, WindStr_y !$OMP parallel default(shared) private(weights,I_wts) !$OMP do do j=jsc-1,jec+1 ; do I=isc-1,iec - weights = (G%areaT(i,j)*ice_cover(i,j) + G%areaT(i+1,j)*ice_cover(i+1,j)) + weights = US%L_to_m**2*(G%areaT(i,j)*ice_cover(i,j) + G%areaT(i+1,j)*ice_cover(i+1,j)) if (G%mask2dCu(I,j) * weights > 0.0) then ; I_wts = 1.0 / weights - WindStr_x_Cu(I,j) = G%mask2dCu(I,j) * & + WindStr_x_Cu(I,j) = G%mask2dCu(I,j) *US%L_to_m**2* & (G%areaT(i,j) * ice_cover(i,j) * WindStr_x_A(i,j) + & G%areaT(i+1,j)*ice_cover(i+1,j)*WindStr_x_A(i+1,j)) * I_wts else WindStr_x_Cu(I,j) = 0.0 endif - weights = (G%areaT(i,j)*ice_free(i,j) + G%areaT(i+1,j)*ice_free(i+1,j)) + weights = US%L_to_m**2*(G%areaT(i,j)*ice_free(i,j) + G%areaT(i+1,j)*ice_free(i+1,j)) if (G%mask2dCu(I,j) * weights > 0.0) then ; I_wts = 1.0 / weights - WindStr_x_ocn_Cu(I,j) = G%mask2dCu(I,j) * & + WindStr_x_ocn_Cu(I,j) = G%mask2dCu(I,j) *US%L_to_m**2* & (G%areaT(i,j) * ice_free(i,j) * WindStr_x_ocn_A(i,j) + & G%areaT(i+1,j)*ice_free(i+1,j)*WindStr_x_ocn_A(i+1,j)) * I_wts else @@ -1893,18 +1902,18 @@ subroutine set_wind_stresses_C(FIA, ice_cover, ice_free, WindStr_x_Cu, WindStr_y !$OMP end do nowait !$OMP do do J=jsc-1,jec ; do i=isc-1,iec+1 - weights = (G%areaT(i,j)*ice_cover(i,j) + G%areaT(i,j+1)*ice_cover(i,j+1)) + weights = US%L_to_m**2*(G%areaT(i,j)*ice_cover(i,j) + G%areaT(i,j+1)*ice_cover(i,j+1)) if (G%mask2dCv(i,J) * weights > 0.0) then ; I_wts = 1.0 / weights - WindStr_y_Cv(i,J) = G%mask2dCv(i,J) * & + WindStr_y_Cv(i,J) = G%mask2dCv(i,J) *US%L_to_m**2* & (G%areaT(i,j) * ice_cover(i,j) * WindStr_y_A(i,j) + & G%areaT(i,j+1)*ice_cover(i,j+1)*WindStr_y_A(i,j+1)) * I_wts else WindStr_y_Cv(i,J) = 0.0 endif - weights = (G%areaT(i,j)*ice_free(i,j) + G%areaT(i,j+1)*ice_free(i,j+1)) + weights = US%L_to_m**2*(G%areaT(i,j)*ice_free(i,j) + G%areaT(i,j+1)*ice_free(i,j+1)) if (weights > 0.0) then ; I_wts = 1.0 / weights - WindStr_y_ocn_Cv(i,J) = G%mask2dCv(i,J) * & + WindStr_y_ocn_Cv(i,J) = G%mask2dCv(i,J) *US%L_to_m**2* & (G%areaT(i,j) * ice_free(i,j) * WindStr_y_ocn_A(i,j) + & G%areaT(i,j+1)*ice_free(i,j+1)*WindStr_y_ocn_A(i,j+1)) * I_wts else @@ -1920,7 +1929,7 @@ end subroutine set_wind_stresses_C !> set_wind_stresses_B determines the wind stresses on the ice and open ocean with !! a B-grid staggering of the points. subroutine set_wind_stresses_B(FIA, ice_cover, ice_free, WindStr_x_B, WindStr_y_B, & - WindStr_x_ocn_B, WindStr_y_ocn_B, G, max_ice_cover) + WindStr_x_ocn_B, WindStr_y_ocn_B, G, US, max_ice_cover) type(fast_ice_avg_type), intent(in) :: FIA !< A type containing averages of fields !! (mostly fluxes) over the fast updates type(SIS_hor_grid_type), intent(in) :: G !< The horizontal grid type @@ -1933,6 +1942,7 @@ subroutine set_wind_stresses_B(FIA, ice_cover, ice_free, WindStr_x_B, WindStr_y_ WindStr_y_B, & !< averaged over the ice categories on a B-grid [Pa]. WindStr_x_ocn_B, & !< Zonal wind stress on the ice-free ocean on a B-grid [Pa]. WindStr_y_ocn_B !< Meridional wind stress on the ice-free ocean on a B-grid [Pa]. + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors real, intent(in) :: max_ice_cover !< The fractional ice coverage !! that is close enough to 1 to be complete for the purpose of calculating !! wind stresses [nondim]. @@ -1984,30 +1994,30 @@ subroutine set_wind_stresses_B(FIA, ice_cover, ice_free, WindStr_x_B, WindStr_y_ !$OMP parallel do default(shared) private(weights,I_wts) do J=jsc-1,jec ; do I=isc-1,iec ; if (G%mask2dBu(I,J) > 0.0) then - weights = ((G%areaT(i+1,j+1)*ice_cover(i+1,j+1) + G%areaT(i,j)*ice_cover(i,j)) + & + weights = US%L_to_m**2*((G%areaT(i+1,j+1)*ice_cover(i+1,j+1) + G%areaT(i,j)*ice_cover(i,j)) + & (G%areaT(i+1,j)*ice_cover(i+1,j) + G%areaT(i,j+1)*ice_cover(i,j+1)) ) I_wts = 0.0 ; if (weights > 0.0) I_wts = 1.0 / weights - WindStr_x_B(I,J) = G%mask2dBu(I,J) * & + WindStr_x_B(I,J) = G%mask2dBu(I,J) *US%L_to_m**2* & ((G%areaT(i+1,j+1)*ice_cover(i+1,j+1)*WindStr_x_A(i+1,j+1) + & G%areaT(i,j) * ice_cover(i,j) * WindStr_x_A(i,j)) + & (G%areaT(i+1,j) * ice_cover(i+1,j) * WindStr_x_A(i+1,j) + & G%areaT(i,j+1) * ice_cover(i,j+1) * WindStr_x_A(i,j+1)) ) * I_wts - WindStr_y_B(I,J) = G%mask2dBu(I,J) * & + WindStr_y_B(I,J) = G%mask2dBu(I,J) *US%L_to_m**2* & ((G%areaT(i+1,j+1)*ice_cover(i+1,j+1)*WindStr_y_A(i+1,j+1) + & G%areaT(i,j) * ice_cover(i,j) * WindStr_y_A(i,j)) + & (G%areaT(i+1,j) * ice_cover(i+1,j) * WindStr_y_A(i+1,j) + & G%areaT(i,j+1) * ice_cover(i,j+1) * WindStr_y_A(i,j+1)) ) * I_wts - weights = ((G%areaT(i+1,j+1)*ice_free(i+1,j+1) + G%areaT(i,j)*ice_free(i,j)) + & + weights = US%L_to_m**2*((G%areaT(i+1,j+1)*ice_free(i+1,j+1) + G%areaT(i,j)*ice_free(i,j)) + & (G%areaT(i+1,j)*ice_free(i+1,j) + G%areaT(i,j+1)*ice_free(i,j+1)) ) I_wts = 0.0 ; if (weights > 0.0) I_wts = 1.0 / weights - WindStr_x_ocn_B(I,J) = G%mask2dBu(I,J) * & + WindStr_x_ocn_B(I,J) = G%mask2dBu(I,J) * US%L_to_m**2*& ((G%areaT(i+1,j+1)*ice_free(i+1,j+1)*WindStr_x_ocn_A(i+1,j+1) + & G%areaT(i,j) * ice_free(i,j) * WindStr_x_ocn_A(i,j)) + & (G%areaT(i+1,j) * ice_free(i+1,j) * WindStr_x_ocn_A(i+1,j) + & G%areaT(i,j+1) * ice_free(i,j+1) * WindStr_x_ocn_A(i,j+1)) ) * I_wts - WindStr_y_ocn_B(I,J) = G%mask2dBu(I,J) * & + WindStr_y_ocn_B(I,J) = G%mask2dBu(I,J) *US%L_to_m**2* & ((G%areaT(i+1,j+1)*ice_free(i+1,j+1)*WindStr_y_ocn_A(i+1,j+1) + & G%areaT(i,j) * ice_free(i,j) * WindStr_y_ocn_A(i,j)) + & (G%areaT(i+1,j) * ice_free(i+1,j) * WindStr_y_ocn_A(i+1,j) + & diff --git a/src/SIS_fixed_initialization.F90 b/src/SIS_fixed_initialization.F90 index d377a609..3df873a4 100644 --- a/src/SIS_fixed_initialization.F90 +++ b/src/SIS_fixed_initialization.F90 @@ -20,6 +20,7 @@ module SIS_fixed_initialization use MOM_shared_initialization, only : reset_face_lengths_named, reset_face_lengths_file, reset_face_lengths_list use MOM_shared_initialization, only : read_face_length_list, set_velocity_depth_max, set_velocity_depth_min use MOM_shared_initialization, only : compute_global_grid_integrals, write_ocean_geometry_file +use MOM_unit_scaling, only : unit_scale_type use netcdf @@ -32,8 +33,9 @@ module SIS_fixed_initialization !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! !> SIS_initialize_fixed sets up time-invariant quantities related to SIS's !! horizontal grid, bathymetry, restricted channel widths and the Coriolis parameter. -subroutine SIS_initialize_fixed(G, PF, write_geom, output_dir) +subroutine SIS_initialize_fixed(G, US, PF, write_geom, output_dir) type(dyn_horgrid_type), intent(inout) :: G !< The ocean's grid structure. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(param_file_type), intent(in) :: PF !< A structure indicating the open file !! to parse for model parameter values. logical, intent(in) :: write_geom !< If true, write grid geometry files. @@ -60,7 +62,7 @@ subroutine SIS_initialize_fixed(G, PF, write_geom, output_dir) inputdir = slasher(inputdir) ! Set up the parameters of the physical domain (i.e. the grid), G - call set_grid_metrics(G, PF) + call set_grid_metrics(G, PF, US) ! Set up the bottom depth, G%bathyT, either analytically or from a file call SIS_initialize_topography(G%bathyT, G%max_depth, G, PF) @@ -95,9 +97,9 @@ subroutine SIS_initialize_fixed(G, PF, write_geom, output_dir) default="none") select case ( trim(config) ) case ("none") - case ("list") ; call reset_face_lengths_list(G, PF) - case ("file") ; call reset_face_lengths_file(G, PF) - case ("global_1deg") ; call reset_face_lengths_named(G, PF, trim(config)) + case ("list") ; call reset_face_lengths_list(G, PF, US) + case ("file") ; call reset_face_lengths_file(G, PF, US) + case ("global_1deg") ; call reset_face_lengths_named(G, PF, trim(config), US) case default ; call MOM_error(FATAL, "SIS_initialize_fixed: "// & "Unrecognized channel configuration "//trim(config)) end select @@ -110,8 +112,8 @@ subroutine SIS_initialize_fixed(G, PF, write_geom, output_dir) call MOM_calculate_grad_Coriolis(G%dF_dx, G%dF_dy, G) if (debug) then call Bchksum(G%CoriolisBu, "SIS_initialize_fixed: f ", G%HI) - call hchksum(G%dF_dx, "SIS_initialize_fixed: dF_dx ", G%HI) - call hchksum(G%dF_dy, "SIS_initialize_fixed: dF_dy ", G%HI) + call hchksum(G%dF_dx, "SIS_initialize_fixed: dF_dx ", G%HI, scale=US%m_to_L) + call hchksum(G%dF_dy, "SIS_initialize_fixed: dF_dy ", G%HI, scale=US%m_to_L) endif call initialize_grid_rotation_angle(G, PF) diff --git a/src/SIS_hor_grid.F90 b/src/SIS_hor_grid.F90 index fa296aeb..14240cd4 100644 --- a/src/SIS_hor_grid.F90 +++ b/src/SIS_hor_grid.F90 @@ -16,6 +16,7 @@ module SIS_hor_grid use MOM_domains, only : MOM_domains_init, clone_MOM_domain use MOM_error_handler, only : SIS_error=>MOM_error, FATAL, WARNING, SIS_mesg=>MOM_mesg use MOM_file_parser, only : get_param, log_param, log_version, param_file_type +use MOM_unit_scaling, only : unit_scale_type implicit none ; private @@ -75,12 +76,12 @@ module SIS_hor_grid mask2dT, & !< 0 for land points and 1 for ocean points on the h-grid [nondim]. geoLatT, & !< The geographic latitude at q points in degrees of latitude or m. geoLonT, & !< The geographic longitude at q points in degrees of longitude or m. - dxT, & !< dxT is delta x at h points [m]. - IdxT, & !< 1/dxT [m-1]. - dyT, & !< dyT is delta y at h points [m], and IdyT is 1/dyT [m-1]. - IdyT, & !< dyT is delta y at h points [m], and IdyT is 1/dyT [m-1]. - areaT, & !< The area of an h-cell [m2]. - IareaT !< 1/areaT [m-2]. + dxT, & !< dxT is delta x at h points [L ~> m]. + IdxT, & !< 1/dxT [L-1 ~> m-1]. + dyT, & !< dyT is delta y at h points [L ~> m]. + IdyT, & !< IdyT is 1/dyT [L-1 ~> m-1]. + areaT, & !< The area of an h-cell [L2 ~> m2]. + IareaT !< 1/areaT [L-2 ~> m-2]. real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: sin_rot !< The sine of the angular rotation between the local model grid northward !! and the true northward directions. @@ -92,36 +93,36 @@ module SIS_hor_grid mask2dCu, & !< 0 for boundary points and 1 for ocean points on the u grid [nondim]. geoLatCu, & !< The geographic latitude at u points in degrees of latitude or m. geoLonCu, & !< The geographic longitude at u points in degrees of longitude or m. - dxCu, & !< dxCu is delta x at u points [m]. - IdxCu, & !< 1/dxCu [m-1]. - dyCu, & !< dyCu is delta y at u points [m]. - IdyCu, & !< 1/dyCu [m-1]. - dy_Cu, & !< The unblocked lengths of the u-faces of the h-cell in m. - IareaCu, & !< The masked inverse areas of u-grid cells [m2]. - areaCu !< The areas of the u-grid cells [m2]. + dxCu, & !< dxCu is delta x at u points [L ~> m]. + IdxCu, & !< 1/dxCu [L-1 ~> m-1]. + dyCu, & !< dyCu is delta y at u points [L ~> m]. + IdyCu, & !< 1/dyCu [L-1 ~> m-1]. + dy_Cu, & !< The unblocked lengths of the u-faces of the h-cell [L ~> m]. + IareaCu, & !< The masked inverse areas of u-grid cells [L-2 ~> m-2]. + areaCu !< The areas of the u-grid cells [L2 ~> m2]. real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_) :: & mask2dCv, & !< 0 for boundary points and 1 for ocean points on the v grid [nondim]. geoLatCv, & !< The geographic latitude at v points in degrees of latitude or m. geoLonCv, & !< The geographic longitude at v points in degrees of longitude or m. - dxCv, & !< dxCv is delta x at v points [m]. - IdxCv, & !< 1/dxCv [m-1]. - dyCv, & !< dyCv is delta y at v points [m]. - IdyCv, & !< 1/dyCv [m-1]. - dx_Cv, & !< The unblocked lengths of the v-faces of the h-cell in m. - IareaCv, & !< The masked inverse areas of v-grid cells [m2]. - areaCv !< The areas of the v-grid cells [m2]. + dxCv, & !< dxCv is delta x at v points [L ~> m]. + IdxCv, & !< 1/dxCv [L-1 ~> m-1]. + dyCv, & !< dyCv is delta y at v points [L ~> m]. + IdyCv, & !< 1/dyCv [L-1 ~> m-1]. + dx_Cv, & !< The unblocked lengths of the v-faces of the h-cell [L ~> m]. + IareaCv, & !< The masked inverse areas of v-grid cells [L-2 ~> m-2]. + areaCv !< The areas of the v-grid cells [L2 ~> m2]. real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: & mask2dBu, & !< 0 for boundary points and 1 for ocean points on the q grid [nondim]. geoLatBu, & !< The geographic latitude at q points in degrees of latitude or m. geoLonBu, & !< The geographic longitude at q points in degrees of longitude or m. - dxBu, & !< dxBu is delta x at q points [m]. - IdxBu, & !< 1/dxBu [m-1]. - dyBu, & !< dyBu is delta y at q points [m]. - IdyBu, & !< 1/dyBu [m-1]. - areaBu, & !< areaBu is the area of a q-cell [m2] - IareaBu !< IareaBu = 1/areaBu [m-2]. + dxBu, & !< dxBu is delta x at q points [L ~> m]. + IdxBu, & !< 1/dxBu [L-1 ~> m-1]. + dyBu, & !< dyBu is delta y at q points [L ~> m]. + IdyBu, & !< 1/dyBu [L-1 ~> m-1]. + areaBu, & !< areaBu is the area of a q-cell [L2 ~> m2] + IareaBu !< IareaBu = 1/areaBu [L-2 ~> m-2]. real, pointer, dimension(:) :: gridLatT => NULL() !< The latitude of T points for the purpose of labeling the output axes. @@ -145,10 +146,12 @@ module SIS_hor_grid real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: & CoriolisBu !< The Coriolis parameter at corner points [s-1]. real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: & - df_dx, & !< Derivative d/dx f (Coriolis parameter) at h-points [s-1 m-1]. - df_dy !< Derivative d/dy f (Coriolis parameter) at h-points [s-1 m-1]. + df_dx, & !< Derivative d/dx f (Coriolis parameter) at h-points [s-1 L-1 ~> s-1 m-1]. + df_dy !< Derivative d/dy f (Coriolis parameter) at h-points [s-1 L-1 ~> s-1 m-1]. real :: g_Earth !< The gravitational acceleration [m s-2]. + type(unit_scale_type), pointer :: US => NULL() !< A dimensional unit scaling type + ! These variables are for block structures. integer :: nblocks !< The number of sub-PE blocks on this PE type(hor_index_type), pointer :: Block(:) => NULL() !< Index ranges for each block diff --git a/src/SIS_slow_thermo.F90 b/src/SIS_slow_thermo.F90 index ac29f68b..c9ecfa06 100644 --- a/src/SIS_slow_thermo.F90 +++ b/src/SIS_slow_thermo.F90 @@ -33,6 +33,7 @@ module SIS_slow_thermo use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint use MOM_file_parser, only : get_param, log_param, log_version, param_file_type use MOM_hor_index, only : hor_index_type +use MOM_unit_scaling, only : unit_scale_type use MOM_EOS, only : EOS_type, calculate_density_derivs use coupler_types_mod, only : coupler_type_spawn, coupler_type_initialized @@ -303,7 +304,7 @@ end subroutine post_flux_diagnostics !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! !> slow_thermodynamics takes care of slow ice thermodynamics and mass changes -subroutine slow_thermodynamics(IST, dt_slow, CS, OSS, FIA, XSF, IOF, G, IG) +subroutine slow_thermodynamics(IST, dt_slow, CS, OSS, FIA, XSF, IOF, G, US, IG) type(ice_state_type), intent(inout) :: IST !< A type describing the state of the sea ice real, intent(in) :: dt_slow !< The thermodynamic step [s]. @@ -318,6 +319,7 @@ subroutine slow_thermodynamics(IST, dt_slow, CS, OSS, FIA, XSF, IOF, G, IG) type(ice_ocean_flux_type), intent(inout) :: IOF !< A structure containing fluxes from the ice to !! the ocean that are calculated by the ice model. type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors type(ice_grid_type), intent(inout) :: IG !< The sea-ice specific grid type ! Local variables @@ -363,7 +365,7 @@ subroutine slow_thermodynamics(IST, dt_slow, CS, OSS, FIA, XSF, IOF, G, IG) call mpp_clock_begin(iceClock7) call accumulate_input_1(IST, FIA, OSS, dt_slow, G, IG, CS%sum_output_CSp) if (CS%column_check) & - call write_ice_statistics(IST, CS%Time, CS%n_calls, G, IG, CS%sum_output_CSp, & + call write_ice_statistics(IST, CS%Time, CS%n_calls, G, US, IG, CS%sum_output_CSp, & message=" SIS_slow_thermo", check_column=.true.) call mpp_clock_end(iceClock7) @@ -467,13 +469,13 @@ subroutine slow_thermodynamics(IST, dt_slow, CS, OSS, FIA, XSF, IOF, G, IG) if (associated(XSF)) call add_excess_fluxes(IOF, XSF, G) if (CS%column_check) & - call write_ice_statistics(IST, CS%Time, CS%n_calls, G, IG, CS%sum_output_CSp, & + call write_ice_statistics(IST, CS%Time, CS%n_calls, G, US, IG, CS%sum_output_CSp, & message=" Post_thermo A", check_column=.true.) call adjust_ice_categories(IST%mH_ice, IST%mH_snow, IST%mH_pond, IST%part_size, & IST%TrReg, G, IG, CS%SIS_transport_CSp) !Niki: add ridging? if (CS%column_check) & - call write_ice_statistics(IST, CS%Time, CS%n_calls, G, IG, CS%sum_output_CSp, & + call write_ice_statistics(IST, CS%Time, CS%n_calls, G, US, IG, CS%sum_output_CSp, & message=" Post_thermo B ", check_column=.true.) end subroutine slow_thermodynamics diff --git a/src/SIS_sum_output.F90 b/src/SIS_sum_output.F90 index b8e6b6b2..d6d0cb1e 100644 --- a/src/SIS_sum_output.F90 +++ b/src/SIS_sum_output.F90 @@ -24,6 +24,7 @@ module SIS_sum_output use MOM_string_functions, only : slasher use MOM_time_manager, only : time_type, get_time, operator(>), operator(-) use MOM_time_manager, only : get_date, get_calendar_type, NO_CALENDAR +use MOM_unit_scaling, only : unit_scale_type use SIS_types, only : ice_state_type, ice_ocean_flux_type, fast_ice_avg_type use SIS_types, only : ocean_sfc_state_type use SIS_hor_grid, only : SIS_hor_grid_type @@ -205,11 +206,12 @@ end subroutine SIS_sum_output_end !> Write out the sea ice statistics of the total sea-ice mass, heat and salt by !! hemisphere and other globally integrated quantities. -subroutine write_ice_statistics(IST, day, n, G, IG, CS, message, check_column, tracer_CSp) +subroutine write_ice_statistics(IST, day, n, G, US, IG, CS, message, check_column, tracer_CSp) type(ice_state_type), intent(inout) :: IST !< A type describing the state of the sea ice type(time_type), intent(inout) :: day !< The current model time. integer, intent(in) :: n !< The time step number of the current execution type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors type(ice_grid_type), intent(inout) :: IG !< The sea-ice specific grid type type(SIS_sum_out_CS), pointer :: CS !< The control structure returned by a previous !! call to SIS_sum_output_init @@ -428,7 +430,7 @@ subroutine write_ice_statistics(IST, day, n, G, IG, CS, message, check_column, t do j=js,je ; do i=is,ie hem = 1 ; if (G%geolatT(i,j) < 0.0) hem = 2 do k=1,ncat ; if (G%mask2dT(i,j) * IST%part_size(i,j,k) > 0.0) then - area_pt = G%areaT(i,j) * G%mask2dT(i,j) * IST%part_size(i,j,k) + area_pt = US%L_to_m**2*G%areaT(i,j) * G%mask2dT(i,j) * IST%part_size(i,j,k) ice_area(i,j,hem) = ice_area(i,j,hem) + area_pt col_mass(i,j,hem) = col_mass(i,j,hem) + area_pt * IG%H_to_kg_m2 * & @@ -446,11 +448,11 @@ subroutine write_ice_statistics(IST, day, n, G, IG, CS, message, check_column, t enddo endif ; enddo if (allocated(IST%snow_to_ocn)) then ; if (IST%snow_to_ocn(i,j) > 0.0) then - area_pt = G%areaT(i,j) * G%mask2dT(i,j) + area_pt = US%L_to_m**2*G%areaT(i,j) * G%mask2dT(i,j) col_mass(i,j,hem) = col_mass(i,j,hem) + area_pt * IST%snow_to_ocn(i,j) col_heat(i,j,hem) = col_heat(i,j,hem) + area_pt * (IST%snow_to_ocn(i,j) * IST%enth_snow_to_ocn(i,j)) endif ; endif - if (ice_area(i,j,hem) > 0.1*G%AreaT(i,j)) ice_extent(i,j,hem) = G%AreaT(i,j) + if (ice_area(i,j,hem) > 0.1*US%L_to_m**2*G%areaT(i,j)) ice_extent(i,j,hem) = US%L_to_m**2*G%areaT(i,j) enddo ; enddo Area = reproducing_sum(ice_area, sums=Area_NS) @@ -471,24 +473,24 @@ subroutine write_ice_statistics(IST, day, n, G, IG, CS, message, check_column, t if (IST%Cgrid_dyn) then if (allocated(IST%u_ice_C)) then ; do j=js,je ; do I=is-1,ie if (IST%u_ice_C(I,j) < 0.0) then - CFL_trans = (-IST%u_ice_C(I,j) * dt_CFL) * (G%dy_Cu(I,j) * G%IareaT(i+1,j)) + CFL_trans = (-IST%u_ice_C(I,j) * dt_CFL) * (G%dy_Cu(I,j) * US%m_to_L*G%IareaT(i+1,j)) else - CFL_trans = (IST%u_ice_C(I,j) * dt_CFL) * (G%dy_Cu(I,j) * G%IareaT(i,j)) + CFL_trans = (IST%u_ice_C(I,j) * dt_CFL) * (G%dy_Cu(I,j) * US%m_to_L*G%IareaT(i,j)) endif max_CFL = max(max_CFL, CFL_trans) enddo ; enddo ; endif if (allocated(IST%v_ice_C)) then ; do J=js-1,je ; do i=is,ie if (IST%v_ice_C(i,J) < 0.0) then - CFL_trans = (-IST%v_ice_C(i,J) * dt_CFL) * (G%dx_Cv(i,J) * G%IareaT(i,j+1)) + CFL_trans = (-IST%v_ice_C(i,J) * dt_CFL) * (G%dx_Cv(i,J) * US%m_to_L*G%IareaT(i,j+1)) else - CFL_trans = (IST%v_ice_C(i,J) * dt_CFL) * (G%dx_Cv(i,J) * G%IareaT(i,j)) + CFL_trans = (IST%v_ice_C(i,J) * dt_CFL) * (G%dx_Cv(i,J) * US%m_to_L*G%IareaT(i,j)) endif max_CFL = max(max_CFL, CFL_trans) enddo ; enddo ; endif elseif (allocated(IST%u_ice_B) .and. allocated(IST%v_ice_B)) then do J=js-1,je ; do I=is-1,ie - CFL_u = abs(IST%u_ice_B(I,J)) * dt_CFL * G%IdxBu(I,J) - CFL_v = abs(IST%v_ice_B(I,J)) * dt_CFL * G%IdyBu(I,J) + CFL_u = abs(IST%u_ice_B(I,J)) * dt_CFL * US%m_to_L * G%IdxBu(I,J) + CFL_v = abs(IST%v_ice_B(I,J)) * dt_CFL * US%m_to_L * G%IdyBu(I,J) max_CFL = max(max_CFL, CFL_u, CFL_v) enddo ; enddo endif @@ -507,7 +509,7 @@ subroutine write_ice_statistics(IST, day, n, G, IG, CS, message, check_column, t CS%heat_prev_EFP = heat_EFP ; CS%net_heat_in_EFP = real_to_EFP(0.0) else do j=js,je ; do i=is,ie - area_h = G%areaT(i,j) * G%mask2dT(i,j) + area_h = US%L_to_m**2*G%areaT(i,j) * G%mask2dT(i,j) CS%water_in_col(i,j) = area_h * CS%water_in_col(i,j) CS%heat_in_col(i,j) = area_h * CS%heat_in_col(i,j) CS%salt_in_col(i,j) = area_h * CS%salt_in_col(i,j) diff --git a/src/SIS_tracer_advect.F90 b/src/SIS_tracer_advect.F90 index fd530512..1a7afa66 100644 --- a/src/SIS_tracer_advect.F90 +++ b/src/SIS_tracer_advect.F90 @@ -10,6 +10,7 @@ module SIS_tracer_advect use MOM_domains, only : pass_var, pass_vector, sum_across_PEs, max_across_PEs use MOM_error_handler, only : SIS_error=>MOM_error, FATAL, WARNING, SIS_mesg=>MOM_mesg use MOM_file_parser, only : get_param, log_version, param_file_type +use MOM_unit_scaling, only : unit_scale_type use SIS_hor_grid, only : SIS_hor_grid_type use ice_grid, only : ice_grid_type use SIS_tracer_registry, only : SIS_tracer_registry_type, SIS_tracer_type, SIS_tracer_chksum @@ -44,7 +45,7 @@ module SIS_tracer_advect contains !> advect_SIS_tracers manages the advection of either the snow or ice tracers -subroutine advect_SIS_tracers(h_prev, h_end, uhtr, vhtr, dt, G, IG, CS, TrReg, snow_tr ) ! (, OBC) +subroutine advect_SIS_tracers(h_prev, h_end, uhtr, vhtr, dt, G, US, IG, CS, TrReg, snow_tr ) ! (, OBC) type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type type(ice_grid_type), intent(in) :: IG !< The sea-ice specific grid type real, dimension(SZI_(G),SZJ_(G),SZCAT_(IG)), & @@ -60,6 +61,7 @@ subroutine advect_SIS_tracers(h_prev, h_end, uhtr, vhtr, dt, G, IG, CS, TrReg, s intent(in) :: vhtr !< Accumulated volume or mass fluxes through !! meridional faces [H m2 s-1 ~> kg s-1]. real, intent(in) :: dt !< Time increment [s]. + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors type(SIS_tracer_advect_CS), pointer :: CS !< The control structure returned by a previous !! call to SIS_tracer_advect_init. type(SIS_tracer_registry_type), pointer :: TrReg !< A pointer to the SIS tracer registry. @@ -80,15 +82,15 @@ subroutine advect_SIS_tracers(h_prev, h_end, uhtr, vhtr, dt, G, IG, CS, TrReg, s call cpu_clock_begin(id_clock_advect) if (snow_tr) then if (CS%use_upwind2d) then - call advect_upwind_2d(TrReg%Tr_snow, h_prev, h_end, uhtr, vhtr, ntr, dt, G, IG) + call advect_upwind_2d(TrReg%Tr_snow, h_prev, h_end, uhtr, vhtr, ntr, dt, G, US, IG) else - call advect_tracer(TrReg%Tr_snow, h_prev, h_end, uhtr, vhtr, ntr, dt, G, IG, CS) + call advect_tracer(TrReg%Tr_snow, h_prev, h_end, uhtr, vhtr, ntr, dt, G, US, IG, CS) endif else if (CS%use_upwind2d) then - call advect_upwind_2d(TrReg%Tr_ice, h_prev, h_end, uhtr, vhtr, ntr, dt, G, IG) + call advect_upwind_2d(TrReg%Tr_ice, h_prev, h_end, uhtr, vhtr, ntr, dt, G, US, IG) else - call advect_tracer(TrReg%Tr_ice, h_prev, h_end, uhtr, vhtr, ntr, dt, G, IG, CS) + call advect_tracer(TrReg%Tr_ice, h_prev, h_end, uhtr, vhtr, ntr, dt, G, US, IG, CS) endif endif call cpu_clock_end(id_clock_advect) @@ -97,7 +99,7 @@ end subroutine advect_SIS_tracers !> This subroutine time steps the tracer concentrations using a monotonic, conservative, !! weakly diffusive scheme. -subroutine advect_tracer(Tr, h_prev, h_end, uhtr, vhtr, ntr, dt, G, IG, CS) ! (, OBC) +subroutine advect_tracer(Tr, h_prev, h_end, uhtr, vhtr, ntr, dt, G, US, IG, CS) ! (, OBC) type(SIS_tracer_type), dimension(ntr), & intent(inout) :: Tr !< The tracer concentrations being advected type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type @@ -116,6 +118,7 @@ subroutine advect_tracer(Tr, h_prev, h_end, uhtr, vhtr, ntr, dt, G, IG, CS) ! (, !! meridional faces [H m2 s-1 ~> kg s-1]. real, intent(in) :: dt !< Time increment [s]. integer, intent(in) :: ntr !< The number of tracers to advect + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors type(SIS_tracer_advect_CS), pointer :: CS !< The control structure returned by a previous !! call to SIS_tracer_advect_init. ! type(ocean_OBC_type), pointer :: OBC ! < This open boundary condition type specifies @@ -192,10 +195,10 @@ subroutine advect_tracer(Tr, h_prev, h_end, uhtr, vhtr, ntr, dt, G, IG, CS) ! (, ! bit of extra mass to avoid nonsensical tracer concentrations. This will ! lead rarely to a very slight non-conservation of tracers, but not mass. do j=js,je; do i=is,ie - hprev(i,j,k) = G%areaT(i,j) * (h_prev(i,j,k) + & + hprev(i,j,k) = US%L_to_m**2*G%areaT(i,j) * (h_prev(i,j,k) + & max(0.0, 1.0e-13*h_prev(i,j,k) - h_end(i,j,k))) if (h_end(i,j,k) - h_prev(i,j,k) + ((uhr(I,j,k) - uhr(I-1,j,k)) + & - (vhr(i,J,k) - vhr(i,J-1,k))) * G%IareaT(i,j) > & + (vhr(i,J,k) - vhr(i,J-1,k))) * US%m_to_L**2*G%IareaT(i,j) > & 1e-10*(h_end(i,j,k) + h_prev(i,j,k))) then !$OMP critical call SIS_error(WARNING, "Apparently inconsistent h_prev, h_end, uhr and vhr in advect_tracer.") @@ -206,12 +209,12 @@ subroutine advect_tracer(Tr, h_prev, h_end, uhtr, vhtr, ntr, dt, G, IG, CS) ! (, !$OMP end do nowait !$OMP do do j=jsd,jed ; do I=isd,ied-1 - uh_neglect(I,j) = h_neglect*MIN(G%areaT(i,j),G%areaT(i+1,j)) + uh_neglect(I,j) = US%L_to_m**2*h_neglect*MIN(G%areaT(i,j),G%areaT(i+1,j)) enddo ; enddo !$OMP end do nowait !$OMP do do J=jsd,jed-1 ; do i=isd,ied - vh_neglect(i,J) = h_neglect*MIN(G%areaT(i,j),G%areaT(i,j+1)) + vh_neglect(i,J) = US%L_to_m**2*h_neglect*MIN(G%areaT(i,j),G%areaT(i,j+1)) enddo ; enddo !$OMP end do nowait !$OMP do @@ -304,12 +307,12 @@ subroutine advect_tracer(Tr, h_prev, h_end, uhtr, vhtr, ntr, dt, G, IG, CS) ! (, if (x_first) then ! First, advect zonally. call advect_x(Tr, hprev, uhr, uh_neglect, domore_u, ntr, nL_max, Idt, & - isv, iev, jsv-stensil, jev+stensil, k, G, IG, & + isv, iev, jsv-stensil, jev+stensil, k, G, US, IG, & CS%usePPM, CS%usePCM) !(, OBC) ! Next, advect meridionally. call advect_y(Tr, hprev, vhr, vh_neglect, domore_v, ntr, nL_max, Idt, & - isv, iev, jsv, jev, k, G, IG, CS%usePPM, CS%usePCM) !(, OBC) + isv, iev, jsv, jev, k, G, US, IG, CS%usePPM, CS%usePCM) !(, OBC) domore_k(k) = 0 do j=jsv-stensil,jev+stensil ; if (domore_u(j,k)) domore_k(k) = 1 ; enddo @@ -317,12 +320,12 @@ subroutine advect_tracer(Tr, h_prev, h_end, uhtr, vhtr, ntr, dt, G, IG, CS) ! (, else ! First, advect meridionally. call advect_y(Tr, hprev, vhr, vh_neglect, domore_v, ntr, nL_max, Idt, & - isv-stensil, iev+stensil, jsv, jev, k, G, IG, & + isv-stensil, iev+stensil, jsv, jev, k, G, US, IG, & CS%usePPM, CS%usePCM) !(, OBC) ! Next, advect zonally. call advect_x(Tr, hprev, uhr, uh_neglect, domore_u, ntr, nL_max, Idt, & - isv, iev, jsv, jev, k, G, IG, CS%usePPM, CS%usePCM) !(, OBC) + isv, iev, jsv, jev, k, G, US, IG, CS%usePPM, CS%usePCM) !(, OBC) domore_k(k) = 0 do j=jsv,jev ; if (domore_u(j,k)) domore_k(k) = 1 ; enddo @@ -349,7 +352,7 @@ subroutine advect_tracer(Tr, h_prev, h_end, uhtr, vhtr, ntr, dt, G, IG, CS) ! (, end subroutine advect_tracer !> advect_scalar does advection of a single scalar tracer field. -subroutine advect_scalar(scalar, h_prev, h_end, uhtr, vhtr, dt, G, IG, CS) ! (, OBC) +subroutine advect_scalar(scalar, h_prev, h_end, uhtr, vhtr, dt, G, US, IG, CS) ! (, OBC) type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type type(ice_grid_type), intent(in) :: IG !< The sea-ice specific grid type real, dimension(SZI_(G),SZJ_(G),SZCAT_(IG)), & @@ -367,6 +370,7 @@ subroutine advect_scalar(scalar, h_prev, h_end, uhtr, vhtr, dt, G, IG, CS) ! (, intent(in) :: vhtr !< Accumulated volume or mass fluxes through !! meridional faces [H m2 s-1 ~> kg s-1]. real, intent(in) :: dt !< Time increment [s]. + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors type(SIS_tracer_advect_CS), pointer :: CS !< The control structure returned by a previous !! call to SIS_tracer_advect_init. @@ -428,9 +432,9 @@ subroutine advect_scalar(scalar, h_prev, h_end, uhtr, vhtr, dt, G, IG, CS) ! (, enddo ; enddo do j=js,je ; do i=is,ie - vol_end = (G%areaT(i,j) * h_end(i,j,k)) + vol_end = (US%L_to_m**2*G%areaT(i,j) * h_end(i,j,k)) Ivol_end = 0.0 ; if (vol_end > 0.0) Ivol_end = 1.0 / vol_end - scalar(i,j,k) = ( (G%areaT(i,j)*h_prev(i,j,k))*scalar(i,j,k) - & + scalar(i,j,k) = ( (US%L_to_m**2*G%areaT(i,j)*h_prev(i,j,k))*scalar(i,j,k) - & ((flux_U2d_x(I,j) - flux_U2d_x(I-1,j)) + & (flux_U2d_y(i,J) - flux_U2d_y(i,J-1))) ) * Ivol_end enddo ; enddo @@ -472,10 +476,10 @@ subroutine advect_scalar(scalar, h_prev, h_end, uhtr, vhtr, dt, G, IG, CS) ! (, ! bit of extra mass to avoid nonsensical tracer concentrations. This will ! lead rarely to a very slight non-conservation of tracers, but not mass. do i=is,ie ; do j=js,je - hprev(i,j,k) = G%areaT(i,j) * (h_prev(i,j,k) + & + hprev(i,j,k) = US%L_to_m**2*G%areaT(i,j) * (h_prev(i,j,k) + & max(0.0, 1.0e-13*h_prev(i,j,k) - h_end(i,j,k))) if (h_end(i,j,k) - h_prev(i,j,k) + ((uhr(I,j,k) - uhr(I-1,j,k)) + & - (vhr(i,J,k) - vhr(i,J-1,k))) * G%IareaT(i,j) > & + (vhr(i,J,k) - vhr(i,J-1,k))) * US%m_to_L**2*G%IareaT(i,j) > & 1e-10*(h_end(i,j,k) + h_prev(i,j,k))) then !$OMP critical call SIS_error(WARNING, "Apparently inconsistent h_prev, h_end, uhr and vhr in advect_tracer.") @@ -486,12 +490,12 @@ subroutine advect_scalar(scalar, h_prev, h_end, uhtr, vhtr, dt, G, IG, CS) ! (, !$OMP end do nowait !$OMP do do j=jsd,jed ; do I=isd,ied-1 - uh_neglect(I,j) = h_neglect*MIN(G%areaT(i,j),G%areaT(i+1,j)) + uh_neglect(I,j) = US%L_to_m**2*h_neglect*MIN(G%areaT(i,j),G%areaT(i+1,j)) enddo ; enddo !$OMP end do nowait !$OMP do do J=jsd,jed-1 ; do i=isd,ied - vh_neglect(i,J) = h_neglect*MIN(G%areaT(i,j),G%areaT(i,j+1)) + vh_neglect(i,J) = US%L_to_m**2*h_neglect*MIN(G%areaT(i,j),G%areaT(i,j+1)) enddo ; enddo !$OMP end parallel @@ -552,11 +556,11 @@ subroutine advect_scalar(scalar, h_prev, h_end, uhtr, vhtr, dt, G, IG, CS) ! (, if (x_first) then ! First, advect zonally. call advect_scalar_x(scalar, hprev, uhr, uh_neglect, domore_u, Idt, & - isv, iev, jsv-stensil, jev+stensil, k, G, IG, CS%usePPM, CS%usePCM) !(, OBC) + isv, iev, jsv-stensil, jev+stensil, k, G, US, IG, CS%usePPM, CS%usePCM) !(, OBC) ! Next, advect meridionally. call advect_scalar_y(scalar, hprev, vhr, vh_neglect, domore_v, Idt, & - isv, iev, jsv, jev, k, G, IG, CS%usePPM, CS%usePCM) !(, OBC) + isv, iev, jsv, jev, k, G, US, IG, CS%usePPM, CS%usePCM) !(, OBC) domore_k(k) = 0 do j=jsv-stensil,jev+stensil ; if (domore_u(j,k)) domore_k(k) = 1 ; enddo @@ -564,11 +568,11 @@ subroutine advect_scalar(scalar, h_prev, h_end, uhtr, vhtr, dt, G, IG, CS) ! (, else ! First, advect meridionally. call advect_scalar_y(scalar, hprev, vhr, vh_neglect, domore_v, Idt, & - isv-stensil, iev+stensil, jsv, jev, k, G, IG, CS%usePPM, CS%usePCM) !(, OBC) + isv-stensil, iev+stensil, jsv, jev, k, G, US, IG, CS%usePPM, CS%usePCM) !(, OBC) ! Next, advect zonally. call advect_scalar_x(scalar, hprev, uhr, uh_neglect, domore_u, Idt, & - isv, iev, jsv, jev, k, G, IG, CS%usePPM, CS%usePCM) !(, OBC) + isv, iev, jsv, jev, k, G, US, IG, CS%usePPM, CS%usePCM) !(, OBC) domore_k(k) = 0 do j=jsv,jev ; if (domore_u(j,k)) domore_k(k) = 1 ; enddo @@ -598,7 +602,7 @@ end subroutine advect_scalar !> advect_scalar_x does 1-d flux-form advection in the x-direction !! using a monotonic piecewise constant, linear, or parabolic scheme. subroutine advect_scalar_x(scalar, hprev, uhr, uh_neglect, domore_u, Idt, & - is, ie, js, je, k, G, IG, usePPM, usePCM) ! (, OBC) + is, ie, js, je, k, G, US, IG, usePPM, usePCM) ! (, OBC) type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type type(ice_grid_type), intent(in) :: IG !< The sea-ice specific grid type real, dimension(SZI_(G),SZJ_(G),SZCAT_(IG)), & @@ -622,6 +626,7 @@ subroutine advect_scalar_x(scalar, hprev, uhr, uh_neglect, domore_u, Idt, & integer, intent(in) :: js !< The starting tracer j-index to work on integer, intent(in) :: je !< The ending tracer j-index to work on integer, intent(in) :: k !< The thickness category to work on + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors logical, intent(in) :: usePPM !< If true, use PPM tracer advection instead of PLM. logical, intent(in) :: usePCM !< If true, use PCM tracer advection instead of PLM. ! This subroutine does 1-d flux-form advection in the zonal direction using @@ -710,17 +715,17 @@ subroutine advect_scalar_x(scalar, hprev, uhr, uh_neglect, domore_u, Idt, & haddW(i) = 0.0 ; haddE(i) = 0.0 if (hnew <= 0.0) then hnew = 0.0 ; do_i(i) = .false. - elseif (hnew < h_neglect*G%areaT(i,j)) then + elseif (hnew < h_neglect*US%L_to_m**2*G%areaT(i,j)) then ! Add a bit of thickness with tracer concentrations that are ! proportional to the mass associated with fluxes and the previous ! mass in the cell. - h_add = h_neglect*G%areaT(i,j) - hnew + h_add = h_neglect*US%L_to_m**2*G%areaT(i,j) - hnew I_htot = 1.0 / (hlst(i) + (abs(uhh(I)) + abs(uhh(I-1)))) hlst(i) = hlst(i) + h_add*(hlst(i)*I_htot) haddW(i) = h_add * (abs(uhh(I-1))*I_htot) haddE(i) = h_add * (abs(uhh(I))*I_htot) - Ihnew(i) = 1.0 / (h_neglect*G%areaT(i,j)) + Ihnew(i) = 1.0 / (h_neglect*US%L_to_m**2*G%areaT(i,j)) else Ihnew(i) = 1.0 / hnew endif @@ -743,7 +748,7 @@ end subroutine advect_scalar_x !> advect_x does 1-d flux-form advection of multiple tracers in the x-direction !! using a monotonic piecewise constant, linear, or parabolic scheme. subroutine advect_x(Tr, hprev, uhr, uh_neglect, domore_u, ntr, nL_max, Idt, & - is, ie, js, je, k, G, IG, usePPM, usePCM) ! (, OBC) + is, ie, js, je, k, G, US, IG, usePPM, usePCM) ! (, OBC) type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type type(ice_grid_type), intent(in) :: IG !< The sea-ice specific grid type type(SIS_tracer_type), dimension(ntr), & @@ -769,6 +774,7 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, domore_u, ntr, nL_max, Idt, & integer, intent(in) :: js !< The starting tracer j-index to work on integer, intent(in) :: je !< The ending tracer j-index to work on integer, intent(in) :: k !< The thickness category to work on + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors logical, intent(in) :: usePPM !< If true, use PPM tracer advection instead of PLM. logical, intent(in) :: usePCM !< If true, use PCM tracer advection instead of PLM. ! This subroutine does 1-d flux-form advection in the zonal direction using @@ -865,17 +871,17 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, domore_u, ntr, nL_max, Idt, & haddW(i) = 0.0 ; haddE(i) = 0.0 if (hnew <= 0.0) then hnew = 0.0 ; do_i(i) = .false. - elseif (hnew < h_neglect*G%areaT(i,j)) then + elseif (hnew < h_neglect*US%L_to_m**2*G%areaT(i,j)) then ! Add a bit of thickness with tracer concentrations that are ! proportional to the mass associated with fluxes and the previous ! mass in the cell. - h_add = h_neglect*G%areaT(i,j) - hnew + h_add = h_neglect*US%L_to_m**2*G%areaT(i,j) - hnew I_htot = 1.0 / (hlst(i) + (abs(uhh(I)) + abs(uhh(I-1)))) hlst(i) = hlst(i) + h_add*(hlst(i)*I_htot) haddW(i) = h_add * (abs(uhh(I-1))*I_htot) haddE(i) = h_add * (abs(uhh(I))*I_htot) - Ihnew(i) = 1.0 / (h_neglect*G%areaT(i,j)) + Ihnew(i) = 1.0 / (h_neglect*US%L_to_m**2*G%areaT(i,j)) else Ihnew(i) = 1.0 / hnew endif @@ -1071,7 +1077,7 @@ end subroutine kernel_PPMH3_Tr_x !> advect_scalar_y does 1-d flux-form advection in the y-direction using a !! monotonic piecewise constant, linear, or parabolic scheme. subroutine advect_scalar_y(scalar, hprev, vhr, vh_neglect, domore_v, Idt, & - is, ie, js, je, k, G, IG, usePPM, usePCM) ! (, OBC) + is, ie, js, je, k, G, US, IG, usePPM, usePCM) ! (, OBC) type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type type(ice_grid_type), intent(in) :: IG !< The sea-ice specific grid type real, dimension(SZI_(G),SZJ_(G),SZCAT_(IG)), & @@ -1095,6 +1101,7 @@ subroutine advect_scalar_y(scalar, hprev, vhr, vh_neglect, domore_v, Idt, & integer, intent(in) :: js !< The starting tracer j-index to work on integer, intent(in) :: je !< The ending tracer j-index to work on integer, intent(in) :: k !< The thickness category to work on + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors logical, intent(in) :: usePPM !< If true, use PPM tracer advection instead of PLM. logical, intent(in) :: usePCM !< If true, use PCM tracer advection instead of PLM. ! This subroutine does 1-d flux-form advection using a monotonic piecewise linear scheme. @@ -1192,17 +1199,17 @@ subroutine advect_scalar_y(scalar, hprev, vhr, vh_neglect, domore_v, Idt, & haddS(i) = 0.0 ; haddN(i) = 0.0 if (hnew <= 0.0) then hnew = 0.0 ; do_i(i) = .false. - elseif (hnew < h_neglect*G%areaT(i,j)) then + elseif (hnew < h_neglect*US%L_to_m**2*G%areaT(i,j)) then ! Add a tiny bit of thickness with tracer concentrations that are ! proportional to the mass associated with fluxes and the previous ! mass in the cell. - h_add = h_neglect*G%areaT(i,j) - hnew + h_add = h_neglect*US%L_to_m**2*G%areaT(i,j) - hnew I_htot = 1.0 / (hlst(i) + (abs(vhh(i,J)) + abs(vhh(i,J-1)))) hlst(i) = hlst(i) + h_add*(hlst(i)*I_htot) haddS(i) = h_add * (abs(vhh(i,J-1))*I_htot) haddN(i) = h_add * (abs(vhh(i,J))*I_htot) - Ihnew(i) = 1.0 / (h_neglect*G%areaT(i,j)) + Ihnew(i) = 1.0 / (h_neglect*US%L_to_m**2*G%areaT(i,j)) else Ihnew(i) = 1.0 / hnew endif @@ -1224,7 +1231,7 @@ end subroutine advect_scalar_y !> advect_y does 1-d flux-form advection of multiple tracers in the y-direction !! using a monotonic piecewise constant, linear, or parabolic scheme. subroutine advect_y(Tr, hprev, vhr, vh_neglect, domore_v, ntr, nL_max, Idt, & - is, ie, js, je, k, G, IG, usePPM, usePCM) ! (, OBC) + is, ie, js, je, k, G, US, IG, usePPM, usePCM) ! (, OBC) type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type type(ice_grid_type), intent(in) :: IG !< The sea-ice specific grid type type(SIS_tracer_type), dimension(ntr), & @@ -1250,6 +1257,7 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, domore_v, ntr, nL_max, Idt, & integer, intent(in) :: js !< The starting tracer j-index to work on integer, intent(in) :: je !< The ending tracer j-index to work on integer, intent(in) :: k !< The thickness category to work on + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors logical, intent(in) :: usePPM !< If true, use PPM tracer advection instead of PLM. logical, intent(in) :: usePCM !< If true, use PCM tracer advection instead of PLM. ! This subroutine does 1-d flux-form advection using a monotonic piecewise @@ -1350,17 +1358,17 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, domore_v, ntr, nL_max, Idt, & haddS(i) = 0.0 ; haddN(i) = 0.0 if (hnew <= 0.0) then hnew = 0.0 ; do_i(i) = .false. - elseif (hnew < h_neglect*G%areaT(i,j)) then + elseif (hnew < h_neglect*US%L_to_m**2*G%areaT(i,j)) then ! Add a tiny bit of thickness with tracer concentrations that are ! proportional to the mass associated with fluxes and the previous ! mass in the cell. - h_add = h_neglect*G%areaT(i,j) - hnew + h_add = h_neglect*US%L_to_m**2*G%areaT(i,j) - hnew I_htot = 1.0 / (hlst(i) + (abs(vhh(i,J)) + abs(vhh(i,J-1)))) hlst(i) = hlst(i) + h_add*(hlst(i)*I_htot) haddS(i) = h_add * (abs(vhh(i,J-1))*I_htot) haddN(i) = h_add * (abs(vhh(i,J))*I_htot) - Ihnew(i) = 1.0 / (h_neglect*G%areaT(i,j)) + Ihnew(i) = 1.0 / (h_neglect*US%L_to_m**2*G%areaT(i,j)) else Ihnew(i) = 1.0 / hnew endif @@ -1572,7 +1580,7 @@ end subroutine kernel_PPMH3_Tr_y !> Advect tracers laterally within their categories using 2-d upwind advection. -subroutine advect_upwind_2d(Tr, h_prev, h_end, uhtr, vhtr, ntr, dt, G, IG) +subroutine advect_upwind_2d(Tr, h_prev, h_end, uhtr, vhtr, ntr, dt, G, US, IG) type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type type(ice_grid_type), intent(in) :: IG !< The sea-ice specific grid type type(SIS_tracer_type), dimension(ntr), & @@ -1591,6 +1599,7 @@ subroutine advect_upwind_2d(Tr, h_prev, h_end, uhtr, vhtr, ntr, dt, G, IG) !! meridional faces [H m2 s-1 ~> kg s-1]. real, intent(in) :: dt !< Time increment [s]. integer, intent(in) :: ntr !< The number of tracers to advect + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors real, dimension(SZIB_(G),SZJ_(G)) :: flux_x ! x-direction tracer fluxes [Conc kg] real, dimension(SZI_(G),SZJB_(G)) :: flux_y ! y-direction tracer fluxes [Conc kg] @@ -1604,7 +1613,7 @@ subroutine advect_upwind_2d(Tr, h_prev, h_end, uhtr, vhtr, ntr, dt, G, IG) ! Reconstruct the old value of h ??? ! if (h_prev(i,j,k) > 0.0) then - ! h_last(i,j,k) = h_end(i,j,k) + dt * G%IareaT(i,j) * & + ! h_last(i,j,k) = h_end(i,j,k) + dt * US%m_to_L**2*G%IareaT(i,j) * & ! ((uh(I,j,k) - uh(I-1,j,k)) + (vh(i,J,k) - vh(i,J-1,k))) ! For now this is just non-directionally split upwind advection. @@ -1622,9 +1631,9 @@ subroutine advect_upwind_2d(Tr, h_prev, h_end, uhtr, vhtr, ntr, dt, G, IG) enddo ; enddo do j=js,je ; do i=is,ie - vol_end = (G%areaT(i,j) * h_end(i,j,k)) + vol_end = (US%L_to_m**2*G%areaT(i,j) * h_end(i,j,k)) Ivol_end = 0.0 ; if (vol_end > 0.0) Ivol_end = 1.0 / vol_end - Tr(m)%t(i,j,k,l) = ( (G%areaT(i,j)*h_prev(i,j,k))*Tr(m)%t(i,j,k,l) - & + Tr(m)%t(i,j,k,l) = ( (US%L_to_m**2*G%areaT(i,j)*h_prev(i,j,k))*Tr(m)%t(i,j,k,l) - & ((flux_x(I,j) - flux_x(I-1,j)) + & (flux_y(i,J) - flux_y(i,J-1))) ) * Ivol_end enddo ; enddo diff --git a/src/SIS_transport.F90 b/src/SIS_transport.F90 index a94b5850..3c69add0 100644 --- a/src/SIS_transport.F90 +++ b/src/SIS_transport.F90 @@ -10,6 +10,7 @@ module SIS_transport use MOM_file_parser, only : get_param, log_param, read_param, log_version, param_file_type use MOM_hor_index, only : hor_index_type use MOM_obsolete_params, only : obsolete_logical, obsolete_real +use MOM_unit_scaling, only : unit_scale_type use SIS_continuity, only : SIS_continuity_init, SIS_continuity_end use SIS_continuity, only : continuity=>ice_continuity, SIS_continuity_CS use SIS_continuity, only : summed_continuity, proportionate_continuity @@ -106,7 +107,7 @@ module SIS_transport !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! !> ice_cat_transport does ice transport of mass and tracers by thickness category -subroutine ice_cat_transport(CAS, TrReg, dt_slow, nsteps, G, IG, CS, uc, vc, mca_tot, uh_tot, vh_tot) +subroutine ice_cat_transport(CAS, TrReg, dt_slow, nsteps, G, US, IG, CS, uc, vc, mca_tot, uh_tot, vh_tot) type(cell_average_state_type), intent(inout) :: CAS !< A structure with ocean-cell averaged masses. type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type type(ice_grid_type), intent(inout) :: IG !< The sea-ice specific grid type @@ -115,6 +116,7 @@ subroutine ice_cat_transport(CAS, TrReg, dt_slow, nsteps, G, IG, CS, uc, vc, mca !! ice dynamics are to be advanced [s]. integer, intent(in) :: nsteps !< The number of advective iterations !! to use within this time step. + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors type(SIS_transport_CS), pointer :: CS !< A pointer to the control structure for this module real, dimension(SZIB_(G),SZJ_(G)), optional, intent(in) :: uc !< The zonal ice velocity [m s-1]. real, dimension(SZI_(G),SZJB_(G)), optional, intent(in) :: vc !< The meridional ice velocity [m s-1]. @@ -179,20 +181,20 @@ subroutine ice_cat_transport(CAS, TrReg, dt_slow, nsteps, G, IG, CS, uc, vc, mca if (merged_cont) then call proportionate_continuity(mca_tot(:,:,n-1), uh_tot(:,:,n), vh_tot(:,:,n), & - dt_adv, G, IG, CS%continuity_CSp, & + dt_adv, G, US, IG, CS%continuity_CSp, & h1=CAS%m_ice, uh1=uh_ice, vh1=vh_ice, & h2=CAS%m_snow, uh2=uh_snow, vh2=vh_snow, & h3=CAS%m_pond, uh3=uh_pond, vh3=vh_pond) else - call continuity(uc, vc, mca0_ice, CAS%m_ice, uh_ice, vh_ice, dt_adv, G, IG, CS%continuity_CSp) - call continuity(uc, vc, mca0_snow, CAS%m_snow, uh_snow, vh_snow, dt_adv, G, IG, CS%continuity_CSp) - call continuity(uc, vc, mca0_pond, CAS%m_pond, uh_pond, vh_pond, dt_adv, G, IG, CS%continuity_CSp) + call continuity(uc, vc, mca0_ice, CAS%m_ice, uh_ice, vh_ice, dt_adv, G, US, IG, CS%continuity_CSp) + call continuity(uc, vc, mca0_snow, CAS%m_snow, uh_snow, vh_snow, dt_adv, G, US, IG, CS%continuity_CSp) + call continuity(uc, vc, mca0_pond, CAS%m_pond, uh_pond, vh_pond, dt_adv, G, US, IG, CS%continuity_CSp) endif - call advect_scalar(CAS%mH_ice, mca0_ice, CAS%m_ice, uh_ice, vh_ice, dt_adv, G, IG, CS%SIS_thick_adv_CSp) - call advect_SIS_tracers(mca0_ice, CAS%m_ice, uh_ice, vh_ice, dt_adv, G, IG, & + call advect_scalar(CAS%mH_ice, mca0_ice, CAS%m_ice, uh_ice, vh_ice, dt_adv, G, US, IG, CS%SIS_thick_adv_CSp) + call advect_SIS_tracers(mca0_ice, CAS%m_ice, uh_ice, vh_ice, dt_adv, G, US, IG, & CS%SIS_tr_adv_CSp, TrReg, snow_tr=.false.) - call advect_SIS_tracers(mca0_snow, CAS%m_snow, uh_snow, vh_snow, dt_adv, G, IG, & + call advect_SIS_tracers(mca0_snow, CAS%m_snow, uh_snow, vh_snow, dt_adv, G, US, IG, & CS%SIS_tr_adv_CSp, TrReg, snow_tr=.true.) ! Accumulated diagnostics @@ -214,12 +216,13 @@ end subroutine ice_cat_transport !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! !> finish_ice_transport completes the ice transport and thickness class redistribution -subroutine finish_ice_transport(CAS, IST, TrReg, G, IG, CS, rdg_rate) +subroutine finish_ice_transport(CAS, IST, TrReg, G, US, IG, CS, rdg_rate) type(cell_average_state_type), intent(inout) :: CAS !< A structure with ocean-cell averaged masses. type(ice_state_type), intent(inout) :: IST !< A type describing the state of the sea ice type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type type(ice_grid_type), intent(inout) :: IG !< The sea-ice specific grid type type(SIS_tracer_registry_type), pointer :: TrReg !< The registry of SIS ice and snow tracers. + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors type(SIS_transport_CS), pointer :: CS !< A pointer to the control structure for this module real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: rdg_rate !< The ice ridging rate [s-1]. @@ -249,7 +252,7 @@ subroutine finish_ice_transport(CAS, IST, TrReg, G, IG, CS, rdg_rate) isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; nCat = IG%CatIce ! Convert the ocean-cell averaged properties back into the ice_state_type. - call cell_ave_state_to_ice_state(CAS, G, IG, CS, IST, TrReg) + call cell_ave_state_to_ice_state(CAS, G, US, IG, CS, IST, TrReg) ! Compress the ice where the fractional coverage exceeds 1, starting with the ! thinnest category, in what amounts to a minimalist version of a sea-ice @@ -321,8 +324,8 @@ subroutine finish_ice_transport(CAS, IST, TrReg, G, IG, CS, rdg_rate) call pass_var(IST%mH_ice, G%Domain, complete=.true.) if (CS%check_conservation) then - call get_total_mass(IST, G, IG, tot_ice, tot_snow, scale=IG%H_to_kg_m2) - call get_total_enthalpy(IST, G, IG, enth_ice, enth_snow, scale=IG%H_to_kg_m2) + call get_total_mass(IST, G, US, IG, tot_ice, tot_snow, scale=IG%H_to_kg_m2) + call get_total_enthalpy(IST, G, US, IG, enth_ice, enth_snow, scale=IG%H_to_kg_m2) if (is_root_pe()) then I_tot_ice = abs(EFP_to_real(CAS%tot_ice)) @@ -373,7 +376,7 @@ subroutine finish_ice_transport(CAS, IST, TrReg, G, IG, CS, rdg_rate) ! if (CS%id_rdgo>0) call post_SIS_data(CS%id_rdgo, rdg_open, diag) ! if (CS%id_rdgv>0) then ! do j=jsc,jec ; do i=isc,iec -! tmp2d(i,j) = rdg_vosh(i,j) * G%areaT(i,j) * G%mask2dT(i,j) +! tmp2d(i,j) = rdg_vosh(i,j) * US%L_to_m**2*G%areaT(i,j) * G%mask2dT(i,j) ! enddo ; enddo ! call post_SIS_data(CS%id_rdgv, tmp2d, diag) ! endif @@ -386,9 +389,10 @@ end subroutine finish_ice_transport !> Determine the whole-cell averaged mass of snow and ice by thickness category based !! on the information in the ice state type. -subroutine ice_state_to_cell_ave_state(IST, G, IG, CS, CAS) +subroutine ice_state_to_cell_ave_state(IST, G, US, IG, CS, CAS) type(ice_state_type), intent(in) :: IST !< A type describing the state of the sea ice type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors type(ice_grid_type), intent(in) :: IG !< The sea-ice specific grid type type(SIS_transport_CS), pointer :: CS !< A pointer to the control structure for this module type(cell_average_state_type), intent(inout) :: CAS !< A structure with ocean-cell averaged masses. @@ -447,16 +451,17 @@ subroutine ice_state_to_cell_ave_state(IST, G, IG, CS, CAS) if (allocated(CAS%vh_sum)) CAS%vh_sum(:,:) = 0.0 if (CS%check_conservation) then ! mw/new - need to update this for pond ? - call get_total_mass(IST, G, IG, CAS%tot_ice, CAS%tot_snow, scale=IG%H_to_kg_m2) - call get_total_enthalpy(IST, G, IG, CAS%enth_ice, CAS%enth_snow, scale=IG%H_to_kg_m2) + call get_total_mass(IST, G, US, IG, CAS%tot_ice, CAS%tot_snow, scale=IG%H_to_kg_m2) + call get_total_enthalpy(IST, G, US, IG, CAS%enth_ice, CAS%enth_snow, scale=IG%H_to_kg_m2) endif end subroutine ice_state_to_cell_ave_state !> Convert the ocean-cell averaged properties back into the ice_state_type. -subroutine cell_ave_state_to_ice_state(CAS, G, IG, CS, IST, TrReg) +subroutine cell_ave_state_to_ice_state(CAS, G, US, IG, CS, IST, TrReg) type(cell_average_state_type), intent(inout) :: CAS !< A structure with ocean-cell averaged masses. type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors type(ice_grid_type), intent(in) :: IG !< The sea-ice specific grid type type(SIS_transport_CS), pointer :: CS !< A pointer to the control structure for this module type(ice_state_type), intent(inout) :: IST !< A type describing the state of the sea ice @@ -482,7 +487,7 @@ subroutine cell_ave_state_to_ice_state(CAS, G, IG, CS, IST, TrReg) do j=jsc,jec ; do k=1,nCat ; do i=isc,iec if (CAS%m_ice(i,j,k) > 0.0) then if (CS%roll_factor * (CAS%mH_ice(i,j,k)*IG%H_to_kg_m2/CS%Rho_Ice)**3 > & - (CAS%m_ice(i,j,k)*IG%H_to_kg_m2/CS%Rho_Ice)*G%areaT(i,j)) then + (CAS%m_ice(i,j,k)*IG%H_to_kg_m2/CS%Rho_Ice)*US%L_to_m**2*G%areaT(i,j)) then ! This ice is thicker than it is wide even if all the ice in a grid ! cell is collected into a single cube, so it will roll. Any snow on ! top will simply be redistributed into a thinner layer, although it @@ -490,7 +495,7 @@ subroutine cell_ave_state_to_ice_state(CAS, G, IG, CS, IST, TrReg) ! thinner so that it melts faster, but it should never be made thinner ! than IG%mH_cat_bound(1). CAS%mH_ice(i,j,k) = max((CS%Rho_ice*IG%kg_m2_to_H) * & - sqrt((CAS%m_ice(i,j,k)*G%areaT(i,j)) / & + sqrt((CAS%m_ice(i,j,k)*US%L_to_m**2*G%areaT(i,j)) / & (CS%roll_factor * CAS%mH_ice(i,j,k)) ), IG%mH_cat_bound(1)) endif @@ -967,9 +972,10 @@ subroutine compress_ice(part_sz, mH_ice, mH_snow, mH_pond, TrReg, G, IG, CS, CAS end subroutine compress_ice !> get_total_mass determines the globally integrated mass of snow and ice -subroutine get_total_mass(IST, G, IG, tot_ice, tot_snow, tot_pond, scale) +subroutine get_total_mass(IST, G, US, IG, tot_ice, tot_snow, tot_pond, scale) type(ice_state_type), intent(in) :: IST !< A type describing the state of the sea ice type(SIS_hor_grid_type), intent(in) :: G !< The horizontal grid type + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors type(ice_grid_type), intent(in) :: IG !< The sea-ice specific grid type type(EFP_type), intent(out) :: tot_ice !< The globally integrated total ice [kg]. type(EFP_type), intent(out) :: tot_snow !< The globally integrated total snow [kg]. @@ -987,12 +993,12 @@ subroutine get_total_mass(IST, G, IG, tot_ice, tot_snow, tot_pond, scale) sum_ice(:,:) = 0.0 sum_snow(:,:) = 0.0 do k=1,IG%CatIce ; do j=jsc,jec ; do i=isc,iec - sum_ice(i,j) = sum_ice(i,j) + G%areaT(i,j) * & + sum_ice(i,j) = sum_ice(i,j) + US%L_to_m**2*G%areaT(i,j) * & (IST%part_size(i,j,k) * (H_to_units*IST%mH_ice(i,j,k))) - sum_snow(i,j) = sum_snow(i,j) + G%areaT(i,j) * & + sum_snow(i,j) = sum_snow(i,j) + US%L_to_m**2*G%areaT(i,j) * & (IST%part_size(i,j,k) * (H_to_units*IST%mH_snow(i,j,k))) if (present(tot_pond)) & - sum_pond(i,j) = sum_pond(i,j) + G%areaT(i,j) * & + sum_pond(i,j) = sum_pond(i,j) + US%L_to_m**2*G%areaT(i,j) * & (IST%part_size(i,j,k) * (H_to_units*IST%mH_pond(i,j,k))) enddo ; enddo ; enddo @@ -1047,9 +1053,10 @@ subroutine cell_mass_from_CAS(CAS, G, IG, mca, scale) end subroutine cell_mass_from_CAS !> get_total_enthalpy determines the globally integrated enthalpy of snow and ice -subroutine get_total_enthalpy(IST, G, IG, enth_ice, enth_snow, scale) +subroutine get_total_enthalpy(IST, G, US, IG, enth_ice, enth_snow, scale) type(ice_state_type), intent(in) :: IST !< A type describing the state of the sea ice type(SIS_hor_grid_type), intent(in) :: G !< The horizontal grid type + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors type(ice_grid_type), intent(in) :: IG !< The sea-ice specific grid type type(EFP_type), intent(out) :: enth_ice !< The globally integrated total ice enthalpy [J]. type(EFP_type), intent(out) :: enth_snow !< The globally integrated total snow enthalpy [J]. @@ -1078,11 +1085,11 @@ subroutine get_total_enthalpy(IST, G, IG, enth_ice, enth_snow, scale) I_Nk = 1.0 / IG%NkIce do m=1,IG%NkIce ; do k=1,IG%CatIce ; do j=jsc,jec ; do i=isc,iec - sum_enth_ice(i,j) = sum_enth_ice(i,j) + (G%areaT(i,j) * & + sum_enth_ice(i,j) = sum_enth_ice(i,j) + (US%L_to_m**2*G%areaT(i,j) * & (((H_to_units*IST%mH_ice(i,j,k))*IST%part_size(i,j,k))*I_Nk)) * heat_ice(i,j,k,m) enddo ; enddo ; enddo ; enddo do k=1,IG%CatIce ; do j=jsc,jec ; do i=isc,iec - sum_enth_snow(i,j) = sum_enth_snow(i,j) + (G%areaT(i,j) * & + sum_enth_snow(i,j) = sum_enth_snow(i,j) + (US%L_to_m**2*G%areaT(i,j) * & ((H_to_units*IST%mH_snow(i,j,k))*IST%part_size(i,j,k))) * heat_snow(i,j,k,1) enddo ; enddo ; enddo !### What about sum_enth_pond? diff --git a/src/SIS_utils.F90 b/src/SIS_utils.F90 index 67ff9751..ed37edde 100644 --- a/src/SIS_utils.F90 +++ b/src/SIS_utils.F90 @@ -8,6 +8,7 @@ module SIS_utils use MOM_error_handler, only : SIS_error=>MOM_error, FATAL, WARNING, SIS_mesg=>MOM_mesg use MOM_error_handler, only : is_root_pe use MOM_time_manager, only : time_type, get_date, get_time, set_date, operator(-) +use MOM_unit_scaling, only : unit_scale_type use SIS_diag_mediator, only : post_SIS_data, SIS_diag_ctrl use SIS_debugging, only : hchksum, Bchksum, uvchksum, hchksum_pair, Bchksum_pair use SIS_debugging, only : check_redundant_B @@ -104,12 +105,12 @@ subroutine ice_line(Time, cn_ocn, sst, G) do j=jsc,jec ; do i=isc,iec x(i,j) = 0.0 if (cn_ocn(i,j)<0.85 .and. n*G%geoLatT(i,j)>0.0) & - x(i,j) = G%mask2dT(i,j)*G%areaT(i,j) + x(i,j) = G%mask2dT(i,j)*G%US%L_to_m**2*G%areaT(i,j) enddo ; enddo gx((n+3)/2) = g_sum(x(isc:iec,jsc:jec))/1e12 enddo gx(3) = g_sum(sst(isc:iec,jsc:jec)*G%mask2dT(isc:iec,jsc:jec)*G%areaT(isc:iec,jsc:jec)) / & - (g_sum(G%mask2dT(isc:iec,jsc:jec)*G%areaT(isc:iec,jsc:jec)) + 1e-10) + (g_sum(G%mask2dT(isc:iec,jsc:jec)*G%areaT(isc:iec,jsc:jec)) + G%US%m_to_L**2*1e-10) ! ! print info every 5 days ! @@ -288,9 +289,10 @@ subroutine post_avg_4d(id, val, part, diag, G, mask, scale, offset, wtd) end subroutine post_avg_4d !> Write checksums of the elements of the sea-ice grid -subroutine ice_grid_chksum(G, haloshift) - type(SIS_hor_grid_type), optional, intent(inout) :: G !< The horizontal grid type - integer, optional, intent(in) :: haloshift !< The size of the halo to check +subroutine ice_grid_chksum(G, US, haloshift) + type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors + integer, optional, intent(in) :: haloshift !< The size of the halo to check integer :: isc, iec, jsc, jec, hs isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec @@ -301,10 +303,10 @@ subroutine ice_grid_chksum(G, haloshift) call hchksum(G%geoLatT, "G%geoLatT", G%HI, haloshift=hs) call hchksum(G%geoLonT, "G%geoLonT", G%HI, haloshift=hs) - call hchksum_pair("G%d[xy]T", G%dxT, G%dyT, G, halos=hs) - call hchksum_pair("G%Id[xy]T", G%IdxT, G%IdyT, G, halos=hs) - call hchksum(G%areaT, "G%areaT", G%HI, haloshift=hs) - call hchksum(G%IareaT, "G%IareaT", G%HI, haloshift=hs) + call hchksum_pair("G%d[xy]T", G%dxT, G%dyT, G, halos=hs, scale=US%L_to_m) + call hchksum_pair("G%Id[xy]T", G%IdxT, G%IdyT, G, halos=hs, scale=US%m_to_L) + call hchksum(G%areaT, "G%areaT", G%HI, haloshift=hs, scale=US%L_to_m**2) + call hchksum(G%IareaT, "G%IareaT", G%HI, haloshift=hs, scale=US%m_to_L**2) call hchksum(G%mask2dT, "G%mask2dT", G%HI, haloshift=hs) call hchksum(G%cos_rot, "G%cos_rot", G%HI) call hchksum(G%sin_rot, "G%sin_rot", G%HI) @@ -314,11 +316,11 @@ subroutine ice_grid_chksum(G, haloshift) call Bchksum(G%geoLatBu, "G%geoLatBu", G%HI, haloshift=hs) call Bchksum(G%geoLonBu, "G%geoLonBu", G%HI, haloshift=hs) - call Bchksum_pair("G%d[xy]Bu", G%dxBu, G%dyBu, G, halos=hs, scalars=.true.) - call Bchksum_pair("G%Id[xy]Bu", G%IdxBu, G%IdyBu, G, halos=hs, scalars=.true.) + call Bchksum_pair("G%d[xy]Bu", G%dxBu, G%dyBu, G, halos=hs, scalars=.true., scale=US%L_to_m) + call Bchksum_pair("G%Id[xy]Bu", G%IdxBu, G%IdyBu, G, halos=hs, scalars=.true., scale=US%m_to_L) - call Bchksum(G%areaBu, "G%areaBu", G%HI, haloshift=hs) - call Bchksum(G%IareaBu, "G%IareaBu", G%HI, haloshift=hs) + call Bchksum(G%areaBu, "G%areaBu", G%HI, haloshift=hs, scale=US%L_to_m**2) + call Bchksum(G%IareaBu, "G%IareaBu", G%HI, haloshift=hs, scale=US%m_to_L**2) call check_redundant_B("G%areaBu", G%areaBu, G, isc-1, iec+1, jsc-1, jec+1) call check_redundant_B("G%IareaBu", G%IareaBu, G, isc-1, iec+1, jsc-1, jec+1) @@ -328,17 +330,17 @@ subroutine ice_grid_chksum(G, haloshift) call uvchksum("G%geoLatC[uv]", G%geoLatCu, G%geoLatCv, G, halos=hs) call uvchksum("G%geolonC[uv]", G%geoLonCu, G%geoLonCv, G, halos=hs) - call uvchksum("G%d[xy]C[uv]", G%dxCu, G%dyCv, G, halos=hs, scalars=.true.) - call uvchksum("G%d[yx]C[uv]", G%dyCu, G%dxCv, G, halos=hs, scalars=.true.) - call uvchksum("G%Id[xy]C[uv]", G%IdxCu, G%IdyCv, G, halos=hs, scalars=.true.) - call uvchksum("G%Id[yx]C[uv]", G%IdyCu, G%IdxCv, G, halos=hs, scalars=.true.) + call uvchksum("G%d[xy]C[uv]", G%dxCu, G%dyCv, G, halos=hs, scalars=.true., scale=US%L_to_m) + call uvchksum("G%d[yx]C[uv]", G%dyCu, G%dxCv, G, halos=hs, scalars=.true., scale=US%L_to_m) + call uvchksum("G%Id[xy]C[uv]", G%IdxCu, G%IdyCv, G, halos=hs, scalars=.true., scale=US%m_to_L) + call uvchksum("G%Id[yx]C[uv]", G%IdyCu, G%IdxCv, G, halos=hs, scalars=.true., scale=US%m_to_L) - call uvchksum("G%areaC[uv]", G%areaCu, G%areaCv, G, halos=hs) - call uvchksum("G%IareaC[uv]", G%IareaCu, G%IareaCv, G, halos=hs) + call uvchksum("G%areaC[uv]", G%areaCu, G%areaCv, G, halos=hs, scale=US%L_to_m**2) + call uvchksum("G%IareaC[uv]", G%IareaCu, G%IareaCv, G, halos=hs, scale=US%m_to_L**2) call hchksum(G%bathyT, "G%bathyT", G%HI, haloshift=hs) call Bchksum(G%CoriolisBu, "G%CoriolisBu", G%HI, haloshift=hs) - call hchksum_pair("G%dF_d[xy]", G%dF_dx, G%dF_dy, G, halos=hs) + call hchksum_pair("G%dF_d[xy]", G%dF_dx, G%dF_dy, G, halos=hs, scale=US%m_to_L) end subroutine ice_grid_chksum diff --git a/src/ice_age_tracer.F90 b/src/ice_age_tracer.F90 index df256100..6a74f6d8 100644 --- a/src/ice_age_tracer.F90 +++ b/src/ice_age_tracer.F90 @@ -21,6 +21,7 @@ module ice_age_tracer use MOM_error_handler, only : SIS_error=>MOM_error, FATAL, WARNING use MOM_error_handler, only : SIS_mesg=>MOM_mesg use MOM_string_functions, only : slasher +use MOM_unit_scaling, only : unit_scale_type use fms_mod, only : read_data use fms_io_mod, only : register_restart_field, restore_state @@ -417,7 +418,7 @@ function ice_age_stock(mi, stocks, G, IG, CS, names, units) avg_tr = avg_tr/IG%NkIce stocks(tr) = stocks(tr) + avg_tr * & - (G%mask2dT(i,j) * G%areaT(i,j) * mi(i,j,k)) + (G%mask2dT(i,j) * G%US%m_to_L**2*G%areaT(i,j) * mi(i,j,k)) enddo ; enddo ; enddo enddo diff --git a/src/ice_model.F90 b/src/ice_model.F90 index c4b65be1..32a61179 100644 --- a/src/ice_model.F90 +++ b/src/ice_model.F90 @@ -46,6 +46,8 @@ module ice_model_mod use MOM_time_manager, only : time_type, time_type_to_real, real_to_time use MOM_time_manager, only : operator(+), operator(-) use MOM_time_manager, only : operator(>), operator(*), operator(/), operator(/=) +use MOM_unit_scaling, only : unit_scale_type, unit_scaling_init +use MOM_unit_scaling, only : unit_scaling_end, fix_restart_unit_scaling use coupler_types_mod, only : coupler_1d_bc_type, coupler_2d_bc_type, coupler_3d_bc_type use coupler_types_mod, only : coupler_type_spawn, coupler_type_initialized @@ -170,6 +172,7 @@ subroutine update_ice_slow_thermo(Ice) type(ice_state_type), pointer :: sIST => NULL() type(fast_ice_avg_type), pointer :: FIA => NULL() type(ice_rad_type), pointer :: Rad => NULL() + type(unit_scale_type), pointer :: US => NULL() real :: dt_slow ! The time step over which to advance the model. integer :: i, j, i2, j2, i_off, j_off @@ -177,7 +180,7 @@ subroutine update_ice_slow_thermo(Ice) "The pointer to Ice%sCS must be associated in update_ice_slow_thermo.") sIST => Ice%sCS%IST ; sIG => Ice%sCS%IG ; sG => Ice%sCS%G ; FIA => Ice%sCS%FIA - Rad => Ice%sCS%Rad + Rad => Ice%sCS%Rad ; US => Ice%sCS%US call mpp_clock_begin(iceClock) ; call mpp_clock_begin(ice_clock_slow) ! Advance the slow PE clock to give the end time of the slow timestep. There @@ -239,7 +242,7 @@ subroutine update_ice_slow_thermo(Ice) endif call slow_thermodynamics(sIST, dt_slow, Ice%sCS%slow_thermo_CSp, Ice%sCS%OSS, FIA, & - Ice%sCS%XSF, Ice%sCS%IOF, sG, sIG) + Ice%sCS%XSF, Ice%sCS%IOF, sG, US, sIG) if (Ice%sCS%debug) then call Ice_public_type_chksum("Before set_ocean_top_fluxes", Ice, check_slow=.true.) call IOF_chksum("Before set_ocean_top_fluxes", Ice%sCS%IOF, sG) @@ -270,13 +273,14 @@ subroutine update_ice_dynamics_trans(Ice, time_step, start_cycle, end_cycle, cyc type(SIS_hor_grid_type), pointer :: sG => NULL() type(ice_state_type), pointer :: sIST => NULL() type(fast_ice_avg_type), pointer :: FIA => NULL() + type(unit_scale_type), pointer :: US => NULL() real :: dt_slow ! The time step over which to advance the model. logical :: do_multi_trans, cycle_start if (.not.associated(Ice%sCS)) call SIS_error(FATAL, & "The pointer to Ice%sCS must be associated in update_ice_dynamics_trans.") - sIST => Ice%sCS%IST ; sIG => Ice%sCS%IG ; sG => Ice%sCS%G ; FIA => Ice%sCS%FIA + sIST => Ice%sCS%IST ; sIG => Ice%sCS%IG ; sG => Ice%sCS%G ; FIA => Ice%sCS%FIA ; US => Ice%sCS%US dt_slow = time_type_to_real(Ice%sCS%Time_step_slow) if (present(time_step)) dt_slow = time_type_to_real(time_step) cycle_start = .true. ; if (present(start_cycle)) cycle_start = start_cycle @@ -308,17 +312,17 @@ subroutine update_ice_dynamics_trans(Ice, time_step, start_cycle, end_cycle, cyc if (Ice%sCS%specified_ice) then ! There is no ice dynamics or transport. call specified_ice_dynamics(sIST, Ice%sCS%OSS, FIA, Ice%sCS%IOF, dt_slow, & - Ice%sCS%specified_ice_CSp, sG, sIG) + Ice%sCS%specified_ice_CSp, sG, US, sIG) elseif (do_multi_trans) then call SIS_multi_dyn_trans(sIST, Ice%sCS%OSS, FIA, Ice%sCS%IOF, dt_slow, Ice%sCS%dyn_trans_CSp, & - Ice%icebergs, sG, sIG, Ice%sCS%SIS_tracer_flow_CSp, & + Ice%icebergs, sG, US, sIG, Ice%sCS%SIS_tracer_flow_CSp, & start_cycle, end_cycle, cycle_length) elseif (Ice%sCS%slab_ice) then ! Use a very old slab ice model. call slab_ice_dyn_trans(sIST, Ice%sCS%OSS, FIA, Ice%sCS%IOF, dt_slow, Ice%sCS%dyn_trans_CSp, & - sG, sIG, Ice%sCS%SIS_tracer_flow_CSp) + sG, US, sIG, Ice%sCS%SIS_tracer_flow_CSp) else ! This is the typical branch used by SIS2. call SIS_dynamics_trans(sIST, Ice%sCS%OSS, FIA, Ice%sCS%IOF, dt_slow, Ice%sCS%dyn_trans_CSp, & - Ice%icebergs, sG, sIG, Ice%sCS%SIS_tracer_flow_CSp) + Ice%icebergs, sG, US, sIG, Ice%sCS%SIS_tracer_flow_CSp) endif ! Set up the stresses and surface pressure in the externally visible structure Ice. @@ -1612,8 +1616,8 @@ subroutine ice_model_init(Ice, Time_Init, Time, Time_step_fast, Time_step_slow, !! tracer fluxes, and can be used to spawn related !! internal variables in the ice model. -! This include declares and sets the variable "version". -#include "version_variable.h" + ! This include declares and sets the variable "version". +# include "version_variable.h" real :: enth_spec_snow, enth_spec_ice real, allocatable :: S_col(:) real :: pi ! pi = 3.1415926... calculated as 4*atan(1) @@ -1631,6 +1635,7 @@ subroutine ice_model_init(Ice, Time_Init, Time, Time_step_fast, Time_step_slow, type(hor_index_type) :: sHI ! A hor_index_type for array extents on slow_ice_PEs. type(dyn_horgrid_type), pointer :: dG => NULL() + type(unit_scale_type), pointer :: US => NULL() ! These pointers are used only for coding convenience on slow PEs. type(SIS_hor_grid_type), pointer :: sG => NULL() type(MOM_domain_type), pointer :: sGD => NULL() @@ -1770,6 +1775,10 @@ subroutine ice_model_init(Ice, Time_Init, Time, Time_step_fast, Time_step_slow, ! Read all relevant parameters and write them to the model log. call log_version(param_file, mdl, version, "") + + ! Determining the internal unit scaling factors for this run. + call unit_scaling_init(param_file, US) + call get_param(param_file, mdl, "SPECIFIED_ICE", specified_ice, & "If true, the ice is specified and there is no dynamics.", & default=.false.) @@ -2007,6 +2016,7 @@ subroutine ice_model_init(Ice, Time_Init, Time, Time_step_fast, Time_step_slow, if (.not.associated(Ice%sCS)) allocate(Ice%sCS) if (.not.associated(Ice%sCS%IG)) allocate(Ice%sCS%IG) if (.not.associated(Ice%sCS%IST)) allocate(Ice%sCS%IST) + Ice%sCS%US => US Ice%sCS%Time = Time ! Set some pointers for convenience. @@ -2054,11 +2064,12 @@ subroutine ice_model_init(Ice, Time_Init, Time, Time_step_fast, Time_step_slow, call clone_MOM_domain(sGD, dG%Domain) ! Set the bathymetry, Coriolis parameter, open channel widths and masks. - call SIS_initialize_fixed(dG, param_file, write_geom_files, dirs%output_directory) + call SIS_initialize_fixed(dG, US, param_file, write_geom_files, dirs%output_directory) call set_hor_grid(sG, param_file, global_indexing=global_indexing) call copy_dyngrid_to_SIS_horgrid(dG, sG) call destroy_dyn_horgrid(dG) + Ice%sCS%G%US => US ! Allocate and register fields for restarts. @@ -2132,7 +2143,7 @@ subroutine ice_model_init(Ice, Time_Init, Time, Time_step_fast, Time_step_slow, isc = sHI%isc ; iec = sHI%iec ; jsc = sHI%jsc ; jec = sHI%jec i_off = LBOUND(Ice%area,1) - sHI%isc ; j_off = LBOUND(Ice%area,2) - sHI%jsc do j=jsc,jec ; do i=isc,iec ; i2 = i+i_off ; j2 = j+j_off - Ice%area(i2,j2) = sG%areaT(i,j) * sG%mask2dT(i,j) + Ice%area(i2,j2) = US%L_to_m**2 * sG%areaT(i,j) * sG%mask2dT(i,j) enddo ; enddo endif ! slow_ice_PE @@ -2154,6 +2165,8 @@ subroutine ice_model_init(Ice, Time_Init, Time, Time_step_fast, Time_step_slow, if (.not.associated(Ice%fCS%IST)) allocate(Ice%fCS%IST) endif + Ice%fCS%US => US + if (single_IST) then Ice%fCS%G => Ice%sCS%G fG => Ice%fCS%G @@ -2181,11 +2194,12 @@ subroutine ice_model_init(Ice, Time_Init, Time, Time_step_fast, Time_step_slow, call clone_MOM_domain(fGD, dG%Domain) ! Set the bathymetry, Coriolis parameter, open channel widths and masks. - call SIS_initialize_fixed(dG, param_file, .false., dirs%output_directory) + call SIS_initialize_fixed(dG, US, param_file, .false., dirs%output_directory) call set_hor_grid(Ice%fCS%G, param_file, global_indexing=global_indexing) call copy_dyngrid_to_SIS_horgrid(dG, Ice%fCS%G) call destroy_dyn_horgrid(dG) + Ice%fCS%G%US => US endif Ice%fCS%bounds_check = bounds_check @@ -2550,7 +2564,7 @@ subroutine ice_model_init(Ice, Time_Init, Time, Time_step_fast, Time_step_slow, ! in the restart files have been initialized. Now call the initialization ! routines for any dependent sub-modules. - call ice_diagnostics_init(Ice%sCS%IOF, Ice%sCS%OSS, Ice%sCS%FIA, sG, sIG, & + call ice_diagnostics_init(Ice%sCS%IOF, Ice%sCS%OSS, Ice%sCS%FIA, sG, US, sIG, & Ice%sCS%diag, Ice%sCS%Time, Cgrid=sIST%Cgrid_dyn) Ice%axes(1:3) = Ice%sCS%diag%axesTc0%handles(1:3) @@ -2599,7 +2613,7 @@ subroutine ice_model_init(Ice, Time_Init, Time, Time_step_fast, Time_step_slow, sGD%X_flags, sGD%Y_flags, time_type_to_real(Time_step_slow), & Time, sG%geoLonBu(isc:iec,jsc:jec), sG%geoLatBu(isc:iec,jsc:jec), & sG%mask2dT(isc-1:iec+1,jsc-1:jec+1), & - sG%dxCv(isc-1:iec+1,jsc-1:jec+1), sG%dyCu(isc-1:iec+1,jsc-1:jec+1), & + US%L_to_m*sG%dxCv(isc-1:iec+1,jsc-1:jec+1), US%L_to_m*sG%dyCu(isc-1:iec+1,jsc-1:jec+1), & Ice%area, sG%cos_rot(isc-1:iec+1,jsc-1:jec+1), & sG%sin_rot(isc-1:iec+1,jsc-1:jec+1), maskmap=sGD%maskmap ) else @@ -2608,20 +2622,20 @@ subroutine ice_model_init(Ice, Time_Init, Time, Time_step_fast, Time_step_slow, sGD%X_flags, sGD%Y_flags, time_type_to_real(Time_step_slow), & Time, sG%geoLonBu(isc:iec,jsc:jec), sG%geoLatBu(isc:iec,jsc:jec), & sG%mask2dT(isc-1:iec+1,jsc-1:jec+1), & - sG%dxCv(isc-1:iec+1,jsc-1:jec+1), sG%dyCu(isc-1:iec+1,jsc-1:jec+1), & + US%L_to_m*sG%dxCv(isc-1:iec+1,jsc-1:jec+1), US%L_to_m*sG%dyCu(isc-1:iec+1,jsc-1:jec+1), & Ice%area, sG%cos_rot(isc-1:iec+1,jsc-1:jec+1), & sG%sin_rot(isc-1:iec+1,jsc-1:jec+1) ) endif endif ! Do any error checking here. - if (Ice%sCS%debug) call ice_grid_chksum(sG, haloshift=1) + if (Ice%sCS%debug) call ice_grid_chksum(sG, US, haloshift=1) if (specified_ice) then - call write_ice_statistics(sIST, Ice%sCS%Time, 0, sG, sIG, & + call write_ice_statistics(sIST, Ice%sCS%Time, 0, sG, US, sIG, & specified_ice_sum_output_CS(Ice%sCS%specified_ice_CSp)) else - call write_ice_statistics(sIST, Ice%sCS%Time, 0, sG, sIG, & + call write_ice_statistics(sIST, Ice%sCS%Time, 0, sG, US, sIG, & SIS_dyn_trans_sum_output_CS(Ice%sCS%dyn_trans_CSp)) endif endif ! slow_ice_PE diff --git a/src/ice_type.F90 b/src/ice_type.F90 index 7dd5bd2a..1b5a6de9 100644 --- a/src/ice_type.F90 +++ b/src/ice_type.F90 @@ -21,6 +21,7 @@ module ice_type_mod use MOM_file_parser, only : param_file_type use MOM_hor_index, only : hor_index_type use MOM_time_manager, only : time_type, time_type_to_real +use MOM_unit_scaling, only : unit_scale_type use SIS_debugging, only : chksum use SIS_diag_mediator, only : SIS_diag_ctrl, post_data=>post_SIS_data use SIS_diag_mediator, only : register_SIS_diag_field @@ -151,7 +152,9 @@ module ice_type_mod !< A pointer to the SIS fast ice update control structure type(SIS_slow_CS), pointer :: sCS => NULL() !< A pointer to the SIS slow ice update control structure - type(restart_file_type), pointer :: Ice_restart => NULL() + type(unit_scale_type), pointer :: US => NULL() + !< structure containing various unit conversion factors + type(restart_file_type), pointer :: Ice_restart => NULL() !< A pointer to the slow ice restart control structure type(restart_file_type), pointer :: Ice_fast_restart => NULL() !< A pointer to the fast ice restart control structure @@ -571,7 +574,7 @@ subroutine ice_stock_pe(Ice, index, value) value = 0.0 do k=1,ncat ; do j=jsc,jec ; do i=isc,iec value = value + kg_H * (IST%mH_ice(i,j,k) + (IST%mH_snow(i,j,k) + IST%mH_pond(i,j,k))) * & - IST%part_size(i,j,k) * (G%areaT(i,j)*G%mask2dT(i,j)) + IST%part_size(i,j,k) * (G%US%L_to_m**2*G%areaT(i,j)*G%mask2dT(i,j)) enddo ; enddo ; enddo case (ISTOCK_HEAT) @@ -579,13 +582,13 @@ subroutine ice_stock_pe(Ice, index, value) if (slab_ice) then do k=1,ncat ; do j=jsc,jec ; do i=isc,iec if (IST%part_size(i,j,k)*IST%mH_ice(i,j,k) > 0.0) then - value = value - (G%areaT(i,j)*G%mask2dT(i,j)) * IST%part_size(i,j,k) * & + value = value - (G%US%L_to_m**2*G%areaT(i,j)*G%mask2dT(i,j)) * IST%part_size(i,j,k) * & (kg_H * IST%mH_ice(i,j,k)) * LI endif enddo ; enddo ; enddo else !### Should this be changed to raise the temperature to 0 degC? do k=1,ncat ; do j=jsc,jec ; do i=isc,iec - part_wt = (G%areaT(i,j)*G%mask2dT(i,j)) * IST%part_size(i,j,k) + part_wt = (G%US%L_to_m**2*G%areaT(i,j)*G%mask2dT(i,j)) * IST%part_size(i,j,k) if (part_wt*IST%mH_ice(i,j,k) > 0.0) then value = value - (part_wt * (kg_H * IST%mH_snow(i,j,k))) * & Energy_melt_enthS(IST%enth_snow(i,j,k,1), 0.0, IST%ITV) @@ -601,7 +604,7 @@ subroutine ice_stock_pe(Ice, index, value) !There is no salt in the snow. value = 0.0 do m=1,NkIce ; do k=1,ncat ; do j=jsc,jec ; do i=isc,iec - value = value + (IST%part_size(i,j,k) * (G%areaT(i,j)*G%mask2dT(i,j))) * & + value = value + (IST%part_size(i,j,k) * (G%US%L_to_m**2*G%areaT(i,j)*G%mask2dT(i,j))) * & (0.001*(kg_H_Nk*IST%mH_ice(i,j,k))) * IST%sal_ice(i,j,k,m) enddo ; enddo ; enddo ; enddo diff --git a/src/slab_ice.F90 b/src/slab_ice.F90 index e30dd958..2d4b80ba 100644 --- a/src/slab_ice.F90 +++ b/src/slab_ice.F90 @@ -10,6 +10,7 @@ module slab_ice ! use MOM_file_parser, only : get_param, log_param, read_param, log_version, param_file_type use MOM_hor_index, only : hor_index_type use MOM_obsolete_params, only : obsolete_logical, obsolete_real +use MOM_unit_scaling, only : unit_scale_type ! use SIS_diag_mediator, only : post_SIS_data, query_SIS_averaging_enabled, SIS_diag_ctrl ! use SIS_diag_mediator, only : register_diag_field=>register_SIS_diag_field, time_type ! use SIS_diag_mediator, only : safe_alloc_alloc @@ -27,7 +28,7 @@ module slab_ice !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! !> Advect an ice tracer or the thickness using a very old slab-ice algorithm !! dating back to the Manabe model. -subroutine slab_ice_advect(uc, vc, trc, stop_lim, dt_slow, G, part_sz, nsteps) +subroutine slab_ice_advect(uc, vc, trc, stop_lim, dt_slow, G, US, part_sz, nsteps) type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type real, dimension(SZIB_(G),SZJ_(G)), intent(in ) :: uc !< x-face advecting velocity [m s-1] real, dimension(SZI_(G),SZJB_(G)), intent(in ) :: vc !< y-face advecting velocity [m s-1] @@ -37,6 +38,7 @@ subroutine slab_ice_advect(uc, vc, trc, stop_lim, dt_slow, G, part_sz, nsteps) real, intent(in ) :: stop_lim !< A tracer amount below which to !! stop advection, in the same units as tr [Conc] real, intent(in ) :: dt_slow !< The time covered by this call [s]. + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: part_sz !< A part size that is set based on !! whether trc (which may be mass) exceeds 0. integer, optional, intent(in ) :: nsteps !< The number of advective substeps. @@ -60,9 +62,9 @@ subroutine slab_ice_advect(uc, vc, trc, stop_lim, dt_slow, G, part_sz, nsteps) if ( avg > stop_lim .and. uc(I,j) * dif > 0.0) then uflx(I,j) = 0.0 elseif ( uc(i,j) > 0.0 ) then - uflx(I,j) = uc(I,j) * trc(i,j) * G%dy_Cu(I,j) + uflx(I,j) = uc(I,j) * trc(i,j) * US%L_to_m*G%dy_Cu(I,j) else - uflx(I,j) = uc(I,j) * trc(i+1,j) * G%dy_Cu(I,j) + uflx(I,j) = uc(I,j) * trc(i+1,j) * US%L_to_m*G%dy_Cu(I,j) endif enddo ; enddo @@ -72,15 +74,15 @@ subroutine slab_ice_advect(uc, vc, trc, stop_lim, dt_slow, G, part_sz, nsteps) if (avg > stop_lim .and. vc(i,J) * dif > 0.0) then vflx(i,J) = 0.0 elseif ( vc(i,J) > 0.0 ) then - vflx(i,J) = vc(i,J) * trc(i,j) * G%dx_Cv(i,J) + vflx(i,J) = vc(i,J) * trc(i,j) * US%L_to_m*G%dx_Cv(i,J) else - vflx(i,J) = vc(i,J) * trc(i,j+1) * G%dx_Cv(i,J) + vflx(i,J) = vc(i,J) * trc(i,j+1) * US%L_to_m*G%dx_Cv(i,J) endif enddo ; enddo do j=jsc,jec ; do i=isc,iec trc(i,j) = trc(i,j) + dt_adv * ((uflx(I-1,j) - uflx(I,j)) + & - (vflx(i,J-1) - vflx(i,J)) ) * G%IareaT(i,j) + (vflx(i,J-1) - vflx(i,J)) ) * US%m_to_L**2*G%IareaT(i,j) enddo ; enddo call pass_var(trc, G%Domain) diff --git a/src/specified_ice.F90 b/src/specified_ice.F90 index a8eeb96a..20ace05d 100644 --- a/src/specified_ice.F90 +++ b/src/specified_ice.F90 @@ -10,10 +10,11 @@ module specified_ice use MOM_domains, only : AGRID, BGRID_NE, CGRID_NE use MOM_error_handler, only : SIS_error=>MOM_error, FATAL, WARNING, SIS_mesg=>MOM_mesg use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint -use MOM_file_parser, only : get_param, read_param, log_param, log_version, param_file_type +use MOM_file_parser, only : get_param, read_param, log_param, log_version, param_file_type use MOM_time_manager, only : time_type, time_type_to_real, real_to_time use MOM_time_manager, only : operator(+), operator(-) use MOM_time_manager, only : operator(>), operator(*), operator(/), operator(/=) +use MOM_unit_scaling, only : unit_scale_type use SIS_diag_mediator, only : enable_SIS_averaging, disable_SIS_averaging use SIS_diag_mediator, only : query_SIS_averaging_enabled, SIS_diag_ctrl @@ -57,7 +58,7 @@ module specified_ice !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! !> specified_ice_dynamics does an update of ice dynamic quantities with specified ice. -subroutine specified_ice_dynamics(IST, OSS, FIA, IOF, dt_slow, CS, G, IG) +subroutine specified_ice_dynamics(IST, OSS, FIA, IOF, dt_slow, CS, G, US, IG) type(ice_state_type), intent(inout) :: IST !< A type describing the state of the sea ice type(ocean_sfc_state_type), intent(in) :: OSS !< A structure containing the arrays that describe @@ -68,6 +69,7 @@ subroutine specified_ice_dynamics(IST, OSS, FIA, IOF, dt_slow, CS, G, IG) !! the ocean that are calculated by the ice model. real, intent(in) :: dt_slow !< The slow ice dynamics timestep [s]. type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type + type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors type(ice_grid_type), intent(inout) :: IG !< The sea-ice specific grid type type(specified_ice_CS), pointer :: CS !< The control structure for the specified_ice module @@ -99,7 +101,7 @@ subroutine specified_ice_dynamics(IST, OSS, FIA, IOF, dt_slow, CS, G, IG) if (CS%bounds_check) call IST_bounds_check(IST, G, IG, "End of specified_ice_dynamics", OSS=OSS) if (CS%Time + real_to_time(0.5*dt_slow) > CS%write_ice_stats_time) then - call write_ice_statistics(IST, CS%Time, CS%n_calls, G, IG, CS%sum_output_CSp) + call write_ice_statistics(IST, CS%Time, CS%n_calls, G, US, IG, CS%sum_output_CSp) CS%write_ice_stats_time = CS%write_ice_stats_time + CS%ice_stats_interval endif