diff --git a/src/core/MOM_checksum_packages.F90 b/src/core/MOM_checksum_packages.F90 index dc60e0888f..2eb84a1c74 100644 --- a/src/core/MOM_checksum_packages.F90 +++ b/src/core/MOM_checksum_packages.F90 @@ -75,10 +75,10 @@ subroutine MOM_state_chksum_5arg(mesg, u, v, h, uh, vh, G, GV, US, haloshift, sy scale_vel = US%L_T_to_m_s ; if (present(vel_scale)) scale_vel = vel_scale call uvchksum(mesg//" [uv]", u, v, G%HI, haloshift=hs, symmetric=sym, & - omit_corners=omit_corners, scale=scale_vel) - call hchksum(h, mesg//" h", G%HI, haloshift=hs, omit_corners=omit_corners, scale=GV%H_to_MKS) + omit_corners=omit_corners, unscale=scale_vel) + call hchksum(h, mesg//" h", G%HI, haloshift=hs, omit_corners=omit_corners, unscale=GV%H_to_MKS) call uvchksum(mesg//" [uv]h", uh, vh, G%HI, haloshift=hs, symmetric=sym, & - omit_corners=omit_corners, scale=GV%H_to_MKS*US%L_to_m**2*US%s_to_T) + omit_corners=omit_corners, unscale=GV%H_to_MKS*US%L_to_m**2*US%s_to_T) end subroutine MOM_state_chksum_5arg ! ============================================================================= @@ -110,8 +110,8 @@ subroutine MOM_state_chksum_3arg(mesg, u, v, h, G, GV, US, haloshift, symmetric, hs = 1 ; if (present(haloshift)) hs = haloshift sym = .false. ; if (present(symmetric)) sym = symmetric call uvchksum(mesg//" u", u, v, G%HI, haloshift=hs, symmetric=sym, & - omit_corners=omit_corners, scale=US%L_T_to_m_s) - call hchksum(h, mesg//" h",G%HI, haloshift=hs, omit_corners=omit_corners, scale=GV%H_to_MKS) + omit_corners=omit_corners, unscale=US%L_T_to_m_s) + call hchksum(h, mesg//" h",G%HI, haloshift=hs, omit_corners=omit_corners, unscale=GV%H_to_MKS) end subroutine MOM_state_chksum_3arg ! ============================================================================= @@ -130,22 +130,22 @@ subroutine MOM_thermo_chksum(mesg, tv, G, US, haloshift, omit_corners) hs=1 ; if (present(haloshift)) hs=haloshift if (associated(tv%T)) & - call hchksum(tv%T, mesg//" T", G%HI, haloshift=hs, omit_corners=omit_corners, scale=US%C_to_degC) + call hchksum(tv%T, mesg//" T", G%HI, haloshift=hs, omit_corners=omit_corners, unscale=US%C_to_degC) if (associated(tv%S)) & - call hchksum(tv%S, mesg//" S", G%HI, haloshift=hs, omit_corners=omit_corners, scale=US%S_to_ppt) + call hchksum(tv%S, mesg//" S", G%HI, haloshift=hs, omit_corners=omit_corners, unscale=US%S_to_ppt) if (associated(tv%frazil)) & call hchksum(tv%frazil, mesg//" frazil", G%HI, haloshift=hs, omit_corners=omit_corners, & - scale=US%Q_to_J_kg*US%R_to_kg_m3*US%Z_to_m) + unscale=US%Q_to_J_kg*US%R_to_kg_m3*US%Z_to_m) if (associated(tv%salt_deficit)) & call hchksum(tv%salt_deficit, mesg//" salt deficit", G%HI, haloshift=hs, omit_corners=omit_corners, & - scale=US%S_to_ppt*US%RZ_to_kg_m2) + unscale=US%S_to_ppt*US%RZ_to_kg_m2) if (associated(tv%varT)) & - call hchksum(tv%varT, mesg//" varT", G%HI, haloshift=hs, omit_corners=omit_corners, scale=US%C_to_degC**2) + call hchksum(tv%varT, mesg//" varT", G%HI, haloshift=hs, omit_corners=omit_corners, unscale=US%C_to_degC**2) if (associated(tv%varS)) & - call hchksum(tv%varS, mesg//" varS", G%HI, haloshift=hs, omit_corners=omit_corners, scale=US%S_to_ppt**2) + call hchksum(tv%varS, mesg//" varS", G%HI, haloshift=hs, omit_corners=omit_corners, unscale=US%S_to_ppt**2) if (associated(tv%covarTS)) & call hchksum(tv%covarTS, mesg//" covarTS", G%HI, haloshift=hs, omit_corners=omit_corners, & - scale=US%S_to_ppt*US%C_to_degC) + unscale=US%S_to_ppt*US%C_to_degC) end subroutine MOM_thermo_chksum @@ -170,18 +170,18 @@ subroutine MOM_surface_chksum(mesg, sfc_state, G, US, haloshift, symmetric) hs = 1 ; if (present(haloshift)) hs = haloshift if (allocated(sfc_state%SST)) call hchksum(sfc_state%SST, mesg//" SST", G%HI, haloshift=hs, & - scale=US%C_to_degC) + unscale=US%C_to_degC) if (allocated(sfc_state%SSS)) call hchksum(sfc_state%SSS, mesg//" SSS", G%HI, haloshift=hs, & - scale=US%S_to_ppt) + unscale=US%S_to_ppt) if (allocated(sfc_state%sea_lev)) call hchksum(sfc_state%sea_lev, mesg//" sea_lev", G%HI, & - haloshift=hs, scale=US%Z_to_m) + haloshift=hs, unscale=US%Z_to_m) if (allocated(sfc_state%Hml)) call hchksum(sfc_state%Hml, mesg//" Hml", G%HI, haloshift=hs, & - scale=US%Z_to_m) + unscale=US%Z_to_m) if (allocated(sfc_state%u) .and. allocated(sfc_state%v)) & call uvchksum(mesg//" SSU", sfc_state%u, sfc_state%v, G%HI, haloshift=hs, symmetric=sym, & - scale=US%L_T_to_m_s) + unscale=US%L_T_to_m_s) if (allocated(sfc_state%frazil)) call hchksum(sfc_state%frazil, mesg//" frazil", G%HI, & - haloshift=hs, scale=US%Q_to_J_kg*US%RZ_to_kg_m2) + haloshift=hs, unscale=US%Q_to_J_kg*US%RZ_to_kg_m2) end subroutine MOM_surface_chksum @@ -232,14 +232,14 @@ subroutine MOM_accel_chksum(mesg, CAu, CAv, PFu, PFv, diffu, diffv, G, GV, US, p ! Note that for the chksum calls to be useful for reproducing across PE ! counts, there must be no redundant points, so all variables use is..ie ! and js...je as their extent. - call uvchksum(mesg//" CA[uv]", CAu, CAv, G%HI, haloshift=0, symmetric=sym, scale=US%L_T2_to_m_s2) - call uvchksum(mesg//" PF[uv]", PFu, PFv, G%HI, haloshift=0, symmetric=sym, scale=US%L_T2_to_m_s2) - call uvchksum(mesg//" diffu", diffu, diffv, G%HI,haloshift=0, symmetric=sym, scale=US%L_T2_to_m_s2) + call uvchksum(mesg//" CA[uv]", CAu, CAv, G%HI, haloshift=0, symmetric=sym, unscale=US%L_T2_to_m_s2) + call uvchksum(mesg//" PF[uv]", PFu, PFv, G%HI, haloshift=0, symmetric=sym, unscale=US%L_T2_to_m_s2) + call uvchksum(mesg//" diffu", diffu, diffv, G%HI,haloshift=0, symmetric=sym, unscale=US%L_T2_to_m_s2) if (present(pbce)) & - call hchksum(pbce, mesg//" pbce",G%HI,haloshift=0, scale=GV%m_to_H*US%L_T_to_m_s**2) + call hchksum(pbce, mesg//" pbce",G%HI,haloshift=0, unscale=GV%m_to_H*US%L_T_to_m_s**2) if (present(u_accel_bt) .and. present(v_accel_bt)) & call uvchksum(mesg//" [uv]_accel_bt", u_accel_bt, v_accel_bt, G%HI,haloshift=0, symmetric=sym, & - scale=US%L_T2_to_m_s2) + unscale=US%L_T2_to_m_s2) end subroutine MOM_accel_chksum ! ============================================================================= diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 83b4ea17cd..4ceb14fe11 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -1262,89 +1262,89 @@ subroutine MOM_forcing_chksum(mesg, fluxes, G, US, haloshift) ! counts, there must be no redundant points, so all variables use is..ie ! and js...je as their extent. if (associated(fluxes%ustar)) & - call hchksum(fluxes%ustar, mesg//" fluxes%ustar", G%HI, haloshift=hshift, scale=US%Z_to_m*US%s_to_T) + call hchksum(fluxes%ustar, mesg//" fluxes%ustar", G%HI, haloshift=hshift, unscale=US%Z_to_m*US%s_to_T) if (associated(fluxes%tau_mag)) & - call hchksum(fluxes%tau_mag, mesg//" fluxes%tau_mag", G%HI, haloshift=hshift, scale=US%RLZ_T2_to_Pa) + call hchksum(fluxes%tau_mag, mesg//" fluxes%tau_mag", G%HI, haloshift=hshift, unscale=US%RLZ_T2_to_Pa) if (associated(fluxes%buoy)) & - call hchksum(fluxes%buoy, mesg//" fluxes%buoy ", G%HI, haloshift=hshift, scale=US%L_to_m**2*US%s_to_T**3) + call hchksum(fluxes%buoy, mesg//" fluxes%buoy ", G%HI, haloshift=hshift, unscale=US%L_to_m**2*US%s_to_T**3) if (associated(fluxes%sw)) & - call hchksum(fluxes%sw, mesg//" fluxes%sw", G%HI, haloshift=hshift, scale=US%QRZ_T_to_W_m2) + call hchksum(fluxes%sw, mesg//" fluxes%sw", G%HI, haloshift=hshift, unscale=US%QRZ_T_to_W_m2) if (associated(fluxes%sw_vis_dir)) & - call hchksum(fluxes%sw_vis_dir, mesg//" fluxes%sw_vis_dir", G%HI, haloshift=hshift, scale=US%QRZ_T_to_W_m2) + call hchksum(fluxes%sw_vis_dir, mesg//" fluxes%sw_vis_dir", G%HI, haloshift=hshift, unscale=US%QRZ_T_to_W_m2) if (associated(fluxes%sw_vis_dif)) & - call hchksum(fluxes%sw_vis_dif, mesg//" fluxes%sw_vis_dif", G%HI, haloshift=hshift, scale=US%QRZ_T_to_W_m2) + call hchksum(fluxes%sw_vis_dif, mesg//" fluxes%sw_vis_dif", G%HI, haloshift=hshift, unscale=US%QRZ_T_to_W_m2) if (associated(fluxes%sw_nir_dir)) & - call hchksum(fluxes%sw_nir_dir, mesg//" fluxes%sw_nir_dir", G%HI, haloshift=hshift, scale=US%QRZ_T_to_W_m2) + call hchksum(fluxes%sw_nir_dir, mesg//" fluxes%sw_nir_dir", G%HI, haloshift=hshift, unscale=US%QRZ_T_to_W_m2) if (associated(fluxes%sw_nir_dif)) & - call hchksum(fluxes%sw_nir_dif, mesg//" fluxes%sw_nir_dif", G%HI, haloshift=hshift, scale=US%QRZ_T_to_W_m2) + call hchksum(fluxes%sw_nir_dif, mesg//" fluxes%sw_nir_dif", G%HI, haloshift=hshift, unscale=US%QRZ_T_to_W_m2) if (associated(fluxes%lw)) & - call hchksum(fluxes%lw, mesg//" fluxes%lw", G%HI, haloshift=hshift, scale=US%QRZ_T_to_W_m2) + call hchksum(fluxes%lw, mesg//" fluxes%lw", G%HI, haloshift=hshift, unscale=US%QRZ_T_to_W_m2) if (associated(fluxes%latent)) & - call hchksum(fluxes%latent, mesg//" fluxes%latent", G%HI, haloshift=hshift, scale=US%QRZ_T_to_W_m2) + call hchksum(fluxes%latent, mesg//" fluxes%latent", G%HI, haloshift=hshift, unscale=US%QRZ_T_to_W_m2) if (associated(fluxes%latent_evap_diag)) & call hchksum(fluxes%latent_evap_diag, mesg//" fluxes%latent_evap_diag", G%HI, & - haloshift=hshift, scale=US%QRZ_T_to_W_m2) + haloshift=hshift, unscale=US%QRZ_T_to_W_m2) if (associated(fluxes%latent_fprec_diag)) & call hchksum(fluxes%latent_fprec_diag, mesg//" fluxes%latent_fprec_diag", G%HI, & - haloshift=hshift, scale=US%QRZ_T_to_W_m2) + haloshift=hshift, unscale=US%QRZ_T_to_W_m2) if (associated(fluxes%latent_frunoff_diag)) & call hchksum(fluxes%latent_frunoff_diag, mesg//" fluxes%latent_frunoff_diag", G%HI, & - haloshift=hshift, scale=US%QRZ_T_to_W_m2) + haloshift=hshift, unscale=US%QRZ_T_to_W_m2) if (associated(fluxes%sens)) & - call hchksum(fluxes%sens, mesg//" fluxes%sens", G%HI, haloshift=hshift, scale=US%QRZ_T_to_W_m2) + call hchksum(fluxes%sens, mesg//" fluxes%sens", G%HI, haloshift=hshift, unscale=US%QRZ_T_to_W_m2) if (associated(fluxes%evap)) & - call hchksum(fluxes%evap, mesg//" fluxes%evap", G%HI, haloshift=hshift, scale=US%RZ_T_to_kg_m2s) + call hchksum(fluxes%evap, mesg//" fluxes%evap", G%HI, haloshift=hshift, unscale=US%RZ_T_to_kg_m2s) if (associated(fluxes%lprec)) & - call hchksum(fluxes%lprec, mesg//" fluxes%lprec", G%HI, haloshift=hshift, scale=US%RZ_T_to_kg_m2s) + call hchksum(fluxes%lprec, mesg//" fluxes%lprec", G%HI, haloshift=hshift, unscale=US%RZ_T_to_kg_m2s) if (associated(fluxes%fprec)) & - call hchksum(fluxes%fprec, mesg//" fluxes%fprec", G%HI, haloshift=hshift, scale=US%RZ_T_to_kg_m2s) + call hchksum(fluxes%fprec, mesg//" fluxes%fprec", G%HI, haloshift=hshift, unscale=US%RZ_T_to_kg_m2s) if (associated(fluxes%vprec)) & - call hchksum(fluxes%vprec, mesg//" fluxes%vprec", G%HI, haloshift=hshift, scale=US%RZ_T_to_kg_m2s) + call hchksum(fluxes%vprec, mesg//" fluxes%vprec", G%HI, haloshift=hshift, unscale=US%RZ_T_to_kg_m2s) if (associated(fluxes%seaice_melt)) & - call hchksum(fluxes%seaice_melt, mesg//" fluxes%seaice_melt", G%HI, haloshift=hshift, scale=US%RZ_T_to_kg_m2s) + call hchksum(fluxes%seaice_melt, mesg//" fluxes%seaice_melt", G%HI, haloshift=hshift, unscale=US%RZ_T_to_kg_m2s) if (associated(fluxes%seaice_melt_heat)) & call hchksum(fluxes%seaice_melt_heat, mesg//" fluxes%seaice_melt_heat", G%HI, & - haloshift=hshift, scale=US%QRZ_T_to_W_m2) + haloshift=hshift, unscale=US%QRZ_T_to_W_m2) if (associated(fluxes%p_surf)) & - call hchksum(fluxes%p_surf, mesg//" fluxes%p_surf", G%HI, haloshift=hshift, scale=US%RL2_T2_to_Pa) + call hchksum(fluxes%p_surf, mesg//" fluxes%p_surf", G%HI, haloshift=hshift, unscale=US%RL2_T2_to_Pa) if (associated(fluxes%u10_sqr)) & - call hchksum(fluxes%u10_sqr, mesg//" fluxes%u10_sqr", G%HI, haloshift=hshift, scale=US%L_to_m**2*US%s_to_T**2) + call hchksum(fluxes%u10_sqr, mesg//" fluxes%u10_sqr", G%HI, haloshift=hshift, unscale=US%L_to_m**2*US%s_to_T**2) if (associated(fluxes%ice_fraction)) & call hchksum(fluxes%ice_fraction, mesg//" fluxes%ice_fraction", G%HI, haloshift=hshift) if (associated(fluxes%salt_flux)) & - call hchksum(fluxes%salt_flux, mesg//" fluxes%salt_flux", G%HI, haloshift=hshift, scale=US%RZ_T_to_kg_m2s) + call hchksum(fluxes%salt_flux, mesg//" fluxes%salt_flux", G%HI, haloshift=hshift, unscale=US%RZ_T_to_kg_m2s) if (associated(fluxes%TKE_tidal)) & - call hchksum(fluxes%TKE_tidal, mesg//" fluxes%TKE_tidal", G%HI, haloshift=hshift, scale=US%RZ3_T3_to_W_m2) + call hchksum(fluxes%TKE_tidal, mesg//" fluxes%TKE_tidal", G%HI, haloshift=hshift, unscale=US%RZ3_T3_to_W_m2) if (associated(fluxes%ustar_tidal)) & - call hchksum(fluxes%ustar_tidal, mesg//" fluxes%ustar_tidal", G%HI, haloshift=hshift, scale=US%Z_to_m*US%s_to_T) + call hchksum(fluxes%ustar_tidal, mesg//" fluxes%ustar_tidal", G%HI, haloshift=hshift, unscale=US%Z_to_m*US%s_to_T) if (associated(fluxes%lrunoff)) & - call hchksum(fluxes%lrunoff, mesg//" fluxes%lrunoff", G%HI, haloshift=hshift, scale=US%RZ_T_to_kg_m2s) + call hchksum(fluxes%lrunoff, mesg//" fluxes%lrunoff", G%HI, haloshift=hshift, unscale=US%RZ_T_to_kg_m2s) if (associated(fluxes%frunoff)) & - call hchksum(fluxes%frunoff, mesg//" fluxes%frunoff", G%HI, haloshift=hshift, scale=US%RZ_T_to_kg_m2s) + call hchksum(fluxes%frunoff, mesg//" fluxes%frunoff", G%HI, haloshift=hshift, unscale=US%RZ_T_to_kg_m2s) if (associated(fluxes%heat_content_lrunoff)) & call hchksum(fluxes%heat_content_lrunoff, mesg//" fluxes%heat_content_lrunoff", G%HI, & - haloshift=hshift, scale=US%QRZ_T_to_W_m2) + haloshift=hshift, unscale=US%QRZ_T_to_W_m2) if (associated(fluxes%heat_content_frunoff)) & call hchksum(fluxes%heat_content_frunoff, mesg//" fluxes%heat_content_frunoff", G%HI, & - haloshift=hshift, scale=US%QRZ_T_to_W_m2) + haloshift=hshift, unscale=US%QRZ_T_to_W_m2) if (associated(fluxes%heat_content_lprec)) & call hchksum(fluxes%heat_content_lprec, mesg//" fluxes%heat_content_lprec", G%HI, & - haloshift=hshift, scale=US%QRZ_T_to_W_m2) + haloshift=hshift, unscale=US%QRZ_T_to_W_m2) if (associated(fluxes%heat_content_fprec)) & call hchksum(fluxes%heat_content_fprec, mesg//" fluxes%heat_content_fprec", G%HI, & - haloshift=hshift, scale=US%QRZ_T_to_W_m2) + haloshift=hshift, unscale=US%QRZ_T_to_W_m2) if (associated(fluxes%heat_content_cond)) & call hchksum(fluxes%heat_content_cond, mesg//" fluxes%heat_content_cond", G%HI, & - haloshift=hshift, scale=US%QRZ_T_to_W_m2) + haloshift=hshift, unscale=US%QRZ_T_to_W_m2) if (associated(fluxes%heat_content_evap)) & call hchksum(fluxes%heat_content_evap, mesg//" fluxes%heat_content_evap", G%HI, & - haloshift=hshift, scale=US%QRZ_T_to_W_m2) + haloshift=hshift, unscale=US%QRZ_T_to_W_m2) if (associated(fluxes%heat_content_massout)) & call hchksum(fluxes%heat_content_massout, mesg//" fluxes%heat_content_massout", G%HI, & - haloshift=hshift, scale=US%QRZ_T_to_W_m2) + haloshift=hshift, unscale=US%QRZ_T_to_W_m2) if (associated(fluxes%heat_content_massin)) & call hchksum(fluxes%heat_content_massin, mesg//" fluxes%heat_content_massin", G%HI, & - haloshift=hshift, scale=US%QRZ_T_to_W_m2) + haloshift=hshift, unscale=US%QRZ_T_to_W_m2) end subroutine MOM_forcing_chksum !> Write out chksums for the driving mechanical forces. @@ -1364,17 +1364,17 @@ subroutine MOM_mech_forcing_chksum(mesg, forces, G, US, haloshift) ! and js...je as their extent. if (associated(forces%taux) .and. associated(forces%tauy)) & call uvchksum(mesg//" forces%tau[xy]", forces%taux, forces%tauy, G%HI, & - haloshift=hshift, symmetric=.true., scale=US%RLZ_T2_to_Pa) + haloshift=hshift, symmetric=.true., unscale=US%RLZ_T2_to_Pa) if (associated(forces%p_surf)) & - call hchksum(forces%p_surf, mesg//" forces%p_surf", G%HI, haloshift=hshift, scale=US%RL2_T2_to_Pa) + call hchksum(forces%p_surf, mesg//" forces%p_surf", G%HI, haloshift=hshift, unscale=US%RL2_T2_to_Pa) if (associated(forces%ustar)) & - call hchksum(forces%ustar, mesg//" forces%ustar", G%HI, haloshift=hshift, scale=US%Z_to_m*US%s_to_T) + call hchksum(forces%ustar, mesg//" forces%ustar", G%HI, haloshift=hshift, unscale=US%Z_to_m*US%s_to_T) if (associated(forces%tau_mag)) & - call hchksum(forces%tau_mag, mesg//" forces%tau_mag", G%HI, haloshift=hshift, scale=US%RLZ_T2_to_Pa) + call hchksum(forces%tau_mag, mesg//" forces%tau_mag", G%HI, haloshift=hshift, unscale=US%RLZ_T2_to_Pa) if (associated(forces%rigidity_ice_u) .and. associated(forces%rigidity_ice_v)) & call uvchksum(mesg//" forces%rigidity_ice_[uv]", forces%rigidity_ice_u, & forces%rigidity_ice_v, G%HI, haloshift=hshift, symmetric=.true., & - scale=US%L_to_m**3*US%L_to_Z*US%s_to_T, scalar_pair=.true.) + unscale=US%L_to_m**3*US%L_to_Z*US%s_to_T, scalar_pair=.true.) end subroutine MOM_mech_forcing_chksum @@ -2657,7 +2657,7 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h enddo ; enddo if (handles%id_prcme > 0) call post_data(handles%id_prcme, res, diag) if (handles%id_total_prcme > 0) then - total_transport = global_area_integral(res, G, scale=US%RZ_T_to_kg_m2s) + total_transport = global_area_integral(res, G, unscale=US%RZ_T_to_kg_m2s) call post_data(handles%id_total_prcme, total_transport, diag) endif if (handles%id_prcme_ga > 0) then @@ -2684,7 +2684,7 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h enddo ; enddo if (handles%id_net_massout > 0) call post_data(handles%id_net_massout, res, diag) if (handles%id_total_net_massout > 0) then - total_transport = global_area_integral(res, G, scale=US%RZ_T_to_kg_m2s) + total_transport = global_area_integral(res, G, unscale=US%RZ_T_to_kg_m2s) call post_data(handles%id_total_net_massout, total_transport, diag) endif endif @@ -2716,7 +2716,7 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h enddo ; enddo if (handles%id_net_massin > 0) call post_data(handles%id_net_massin, res, diag) if (handles%id_total_net_massin > 0) then - total_transport = global_area_integral(res, G, scale=US%RZ_T_to_kg_m2s) + total_transport = global_area_integral(res, G, unscale=US%RZ_T_to_kg_m2s) call post_data(handles%id_total_net_massin, total_transport, diag) endif endif @@ -2727,7 +2727,7 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h if ((handles%id_evap > 0) .and. associated(fluxes%evap)) & call post_data(handles%id_evap, fluxes%evap, diag) if ((handles%id_total_evap > 0) .and. associated(fluxes%evap)) then - total_transport = global_area_integral(fluxes%evap, G, scale=US%RZ_T_to_kg_m2s) + total_transport = global_area_integral(fluxes%evap, G, unscale=US%RZ_T_to_kg_m2s) call post_data(handles%id_total_evap, total_transport, diag) endif if ((handles%id_evap_ga > 0) .and. associated(fluxes%evap)) then @@ -2741,7 +2741,7 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h enddo ; enddo if (handles%id_precip > 0) call post_data(handles%id_precip, res, diag) if (handles%id_total_precip > 0) then - total_transport = global_area_integral(res, G, scale=US%RZ_T_to_kg_m2s) + total_transport = global_area_integral(res, G, unscale=US%RZ_T_to_kg_m2s) call post_data(handles%id_total_precip, total_transport, diag) endif if (handles%id_precip_ga > 0) then @@ -2753,7 +2753,7 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h if (associated(fluxes%lprec)) then if (handles%id_lprec > 0) call post_data(handles%id_lprec, fluxes%lprec, diag) if (handles%id_total_lprec > 0) then - total_transport = global_area_integral(fluxes%lprec, G, scale=US%RZ_T_to_kg_m2s) + total_transport = global_area_integral(fluxes%lprec, G, unscale=US%RZ_T_to_kg_m2s) call post_data(handles%id_total_lprec, total_transport, diag) endif if (handles%id_lprec_ga > 0) then @@ -2765,7 +2765,7 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h if (associated(fluxes%fprec)) then if (handles%id_fprec > 0) call post_data(handles%id_fprec, fluxes%fprec, diag) if (handles%id_total_fprec > 0) then - total_transport = global_area_integral(fluxes%fprec, G, scale=US%RZ_T_to_kg_m2s) + total_transport = global_area_integral(fluxes%fprec, G, unscale=US%RZ_T_to_kg_m2s) call post_data(handles%id_total_fprec, total_transport, diag) endif if (handles%id_fprec_ga > 0) then @@ -2777,7 +2777,7 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h if (associated(fluxes%vprec)) then if (handles%id_vprec > 0) call post_data(handles%id_vprec, fluxes%vprec, diag) if (handles%id_total_vprec > 0) then - total_transport = global_area_integral(fluxes%vprec, G, scale=US%RZ_T_to_kg_m2s) + total_transport = global_area_integral(fluxes%vprec, G, unscale=US%RZ_T_to_kg_m2s) call post_data(handles%id_total_vprec, total_transport, diag) endif if (handles%id_vprec_ga > 0) then @@ -2789,7 +2789,7 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h if (associated(fluxes%lrunoff)) then if (handles%id_lrunoff > 0) call post_data(handles%id_lrunoff, fluxes%lrunoff, diag) if (handles%id_total_lrunoff > 0) then - total_transport = global_area_integral(fluxes%lrunoff, G, scale=US%RZ_T_to_kg_m2s) + total_transport = global_area_integral(fluxes%lrunoff, G, unscale=US%RZ_T_to_kg_m2s) call post_data(handles%id_total_lrunoff, total_transport, diag) endif endif @@ -2797,7 +2797,7 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h if (associated(fluxes%frunoff)) then if (handles%id_frunoff > 0) call post_data(handles%id_frunoff, fluxes%frunoff, diag) if (handles%id_total_frunoff > 0) then - total_transport = global_area_integral(fluxes%frunoff, G, scale=US%RZ_T_to_kg_m2s) + total_transport = global_area_integral(fluxes%frunoff, G, unscale=US%RZ_T_to_kg_m2s) call post_data(handles%id_total_frunoff, total_transport, diag) endif endif @@ -2805,7 +2805,7 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h if (associated(fluxes%seaice_melt)) then if (handles%id_seaice_melt > 0) call post_data(handles%id_seaice_melt, fluxes%seaice_melt, diag) if (handles%id_total_seaice_melt > 0) then - total_transport = global_area_integral(fluxes%seaice_melt, G, scale=US%RZ_T_to_kg_m2s) + total_transport = global_area_integral(fluxes%seaice_melt, G, unscale=US%RZ_T_to_kg_m2s) call post_data(handles%id_total_seaice_melt, total_transport, diag) endif endif @@ -2815,63 +2815,63 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h if ((handles%id_heat_content_lrunoff > 0) .and. associated(fluxes%heat_content_lrunoff)) & call post_data(handles%id_heat_content_lrunoff, fluxes%heat_content_lrunoff, diag) if ((handles%id_total_heat_content_lrunoff > 0) .and. associated(fluxes%heat_content_lrunoff)) then - total_transport = global_area_integral(fluxes%heat_content_lrunoff, G, scale=US%QRZ_T_to_W_m2) + total_transport = global_area_integral(fluxes%heat_content_lrunoff, G, unscale=US%QRZ_T_to_W_m2) call post_data(handles%id_total_heat_content_lrunoff, total_transport, diag) endif if ((handles%id_heat_content_frunoff > 0) .and. associated(fluxes%heat_content_frunoff)) & call post_data(handles%id_heat_content_frunoff, fluxes%heat_content_frunoff, diag) if ((handles%id_total_heat_content_frunoff > 0) .and. associated(fluxes%heat_content_frunoff)) then - total_transport = global_area_integral(fluxes%heat_content_frunoff, G, scale=US%QRZ_T_to_W_m2) + total_transport = global_area_integral(fluxes%heat_content_frunoff, G, unscale=US%QRZ_T_to_W_m2) call post_data(handles%id_total_heat_content_frunoff, total_transport, diag) endif if ((handles%id_heat_content_lprec > 0) .and. associated(fluxes%heat_content_lprec)) & call post_data(handles%id_heat_content_lprec, fluxes%heat_content_lprec, diag) if ((handles%id_total_heat_content_lprec > 0) .and. associated(fluxes%heat_content_lprec)) then - total_transport = global_area_integral(fluxes%heat_content_lprec, G, scale=US%QRZ_T_to_W_m2) + total_transport = global_area_integral(fluxes%heat_content_lprec, G, unscale=US%QRZ_T_to_W_m2) call post_data(handles%id_total_heat_content_lprec, total_transport, diag) endif if ((handles%id_heat_content_fprec > 0) .and. associated(fluxes%heat_content_fprec)) & call post_data(handles%id_heat_content_fprec, fluxes%heat_content_fprec, diag) if ((handles%id_total_heat_content_fprec > 0) .and. associated(fluxes%heat_content_fprec)) then - total_transport = global_area_integral(fluxes%heat_content_fprec, G, scale=US%QRZ_T_to_W_m2) + total_transport = global_area_integral(fluxes%heat_content_fprec, G, unscale=US%QRZ_T_to_W_m2) call post_data(handles%id_total_heat_content_fprec, total_transport, diag) endif if ((handles%id_heat_content_vprec > 0) .and. associated(fluxes%heat_content_vprec)) & call post_data(handles%id_heat_content_vprec, fluxes%heat_content_vprec, diag) if ((handles%id_total_heat_content_vprec > 0) .and. associated(fluxes%heat_content_vprec)) then - total_transport = global_area_integral(fluxes%heat_content_vprec, G, scale=US%QRZ_T_to_W_m2) + total_transport = global_area_integral(fluxes%heat_content_vprec, G, unscale=US%QRZ_T_to_W_m2) call post_data(handles%id_total_heat_content_vprec, total_transport, diag) endif if ((handles%id_heat_content_cond > 0) .and. associated(fluxes%heat_content_cond)) & call post_data(handles%id_heat_content_cond, fluxes%heat_content_cond, diag) if ((handles%id_total_heat_content_cond > 0) .and. associated(fluxes%heat_content_cond)) then - total_transport = global_area_integral(fluxes%heat_content_cond, G, scale=US%QRZ_T_to_W_m2) + total_transport = global_area_integral(fluxes%heat_content_cond, G, unscale=US%QRZ_T_to_W_m2) call post_data(handles%id_total_heat_content_cond, total_transport, diag) endif if ((handles%id_heat_content_evap > 0) .and. associated(fluxes%heat_content_evap)) & call post_data(handles%id_heat_content_evap, fluxes%heat_content_evap, diag) if ((handles%id_total_heat_content_evap > 0) .and. associated(fluxes%heat_content_evap)) then - total_transport = global_area_integral(fluxes%heat_content_evap, G, scale=US%QRZ_T_to_W_m2) + total_transport = global_area_integral(fluxes%heat_content_evap, G, unscale=US%QRZ_T_to_W_m2) call post_data(handles%id_total_heat_content_evap, total_transport, diag) endif if ((handles%id_heat_content_massout > 0) .and. associated(fluxes%heat_content_massout)) & call post_data(handles%id_heat_content_massout, fluxes%heat_content_massout, diag) if ((handles%id_total_heat_content_massout > 0) .and. associated(fluxes%heat_content_massout)) then - total_transport = global_area_integral(fluxes%heat_content_massout, G, scale=US%QRZ_T_to_W_m2) + total_transport = global_area_integral(fluxes%heat_content_massout, G, unscale=US%QRZ_T_to_W_m2) call post_data(handles%id_total_heat_content_massout, total_transport, diag) endif if ((handles%id_heat_content_massin > 0) .and. associated(fluxes%heat_content_massin)) & call post_data(handles%id_heat_content_massin, fluxes%heat_content_massin, diag) if ((handles%id_total_heat_content_massin > 0) .and. associated(fluxes%heat_content_massin)) then - total_transport = global_area_integral(fluxes%heat_content_massin, G, scale=US%QRZ_T_to_W_m2) + total_transport = global_area_integral(fluxes%heat_content_massin, G, unscale=US%QRZ_T_to_W_m2) call post_data(handles%id_total_heat_content_massin, total_transport, diag) endif @@ -2887,7 +2887,7 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h enddo ; enddo if (handles%id_net_heat_coupler > 0) call post_data(handles%id_net_heat_coupler, res, diag) if (handles%id_total_net_heat_coupler > 0) then - total_transport = global_area_integral(res, G, scale=US%QRZ_T_to_W_m2) + total_transport = global_area_integral(res, G, unscale=US%QRZ_T_to_W_m2) call post_data(handles%id_total_net_heat_coupler, total_transport, diag) endif if (handles%id_net_heat_coupler_ga > 0) then @@ -2930,7 +2930,7 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h if (handles%id_net_heat_surface > 0) call post_data(handles%id_net_heat_surface, res, diag) if (handles%id_total_net_heat_surface > 0) then - total_transport = global_area_integral(res, G, scale=US%QRZ_T_to_W_m2) + total_transport = global_area_integral(res, G, unscale=US%QRZ_T_to_W_m2) call post_data(handles%id_total_net_heat_surface, total_transport, diag) endif if (handles%id_net_heat_surface_ga > 0) then @@ -2956,7 +2956,7 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h enddo ; enddo if (handles%id_heat_content_surfwater > 0) call post_data(handles%id_heat_content_surfwater, res, diag) if (handles%id_total_heat_content_surfwater > 0) then - total_transport = global_area_integral(res, G, scale=US%QRZ_T_to_W_m2) + total_transport = global_area_integral(res, G, unscale=US%QRZ_T_to_W_m2) call post_data(handles%id_total_heat_content_surfwater, total_transport, diag) endif endif @@ -2995,7 +2995,7 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h do j=js,je ; do i=is,ie res(i,j) = (fluxes%lw(i,j) + fluxes%latent(i,j)) + fluxes%sens(i,j) enddo ; enddo - total_transport = global_area_integral(res, G, scale=US%QRZ_T_to_W_m2) + total_transport = global_area_integral(res, G, unscale=US%QRZ_T_to_W_m2) call post_data(handles%id_total_LwLatSens, total_transport, diag) endif @@ -3020,7 +3020,7 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h call post_data(handles%id_sw_nir, fluxes%sw_nir_dir+fluxes%sw_nir_dif, diag) endif if ((handles%id_total_sw > 0) .and. associated(fluxes%sw)) then - total_transport = global_area_integral(fluxes%sw, G, scale=US%QRZ_T_to_W_m2) + total_transport = global_area_integral(fluxes%sw, G, unscale=US%QRZ_T_to_W_m2) call post_data(handles%id_total_sw, total_transport, diag) endif if ((handles%id_sw_ga > 0) .and. associated(fluxes%sw)) then @@ -3032,7 +3032,7 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h call post_data(handles%id_lw, fluxes%lw, diag) endif if ((handles%id_total_lw > 0) .and. associated(fluxes%lw)) then - total_transport = global_area_integral(fluxes%lw, G, scale=US%QRZ_T_to_W_m2) + total_transport = global_area_integral(fluxes%lw, G, unscale=US%QRZ_T_to_W_m2) call post_data(handles%id_total_lw, total_transport, diag) endif if ((handles%id_lw_ga > 0) .and. associated(fluxes%lw)) then @@ -3044,7 +3044,7 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h call post_data(handles%id_lat, fluxes%latent, diag) endif if ((handles%id_total_lat > 0) .and. associated(fluxes%latent)) then - total_transport = global_area_integral(fluxes%latent, G, scale=US%QRZ_T_to_W_m2) + total_transport = global_area_integral(fluxes%latent, G, unscale=US%QRZ_T_to_W_m2) call post_data(handles%id_total_lat, total_transport, diag) endif if ((handles%id_lat_ga > 0) .and. associated(fluxes%latent)) then @@ -3056,7 +3056,7 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h call post_data(handles%id_lat_evap, fluxes%latent_evap_diag, diag) endif if ((handles%id_total_lat_evap > 0) .and. associated(fluxes%latent_evap_diag)) then - total_transport = global_area_integral(fluxes%latent_evap_diag, G, scale=US%QRZ_T_to_W_m2) + total_transport = global_area_integral(fluxes%latent_evap_diag, G, unscale=US%QRZ_T_to_W_m2) call post_data(handles%id_total_lat_evap, total_transport, diag) endif @@ -3064,7 +3064,7 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h call post_data(handles%id_lat_fprec, fluxes%latent_fprec_diag, diag) endif if ((handles%id_total_lat_fprec > 0) .and. associated(fluxes%latent_fprec_diag)) then - total_transport = global_area_integral(fluxes%latent_fprec_diag, G, scale=US%QRZ_T_to_W_m2) + total_transport = global_area_integral(fluxes%latent_fprec_diag, G, unscale=US%QRZ_T_to_W_m2) call post_data(handles%id_total_lat_fprec, total_transport, diag) endif @@ -3072,7 +3072,7 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h call post_data(handles%id_lat_frunoff, fluxes%latent_frunoff_diag, diag) endif if (handles%id_total_lat_frunoff > 0 .and. associated(fluxes%latent_frunoff_diag)) then - total_transport = global_area_integral(fluxes%latent_frunoff_diag, G, scale=US%QRZ_T_to_W_m2) + total_transport = global_area_integral(fluxes%latent_frunoff_diag, G, unscale=US%QRZ_T_to_W_m2) call post_data(handles%id_total_lat_frunoff, total_transport, diag) endif @@ -3085,12 +3085,12 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h endif if ((handles%id_total_seaice_melt_heat > 0) .and. associated(fluxes%seaice_melt_heat)) then - total_transport = global_area_integral(fluxes%seaice_melt_heat, G, scale=US%QRZ_T_to_W_m2) + total_transport = global_area_integral(fluxes%seaice_melt_heat, G, unscale=US%QRZ_T_to_W_m2) call post_data(handles%id_total_seaice_melt_heat, total_transport, diag) endif if ((handles%id_total_sens > 0) .and. associated(fluxes%sens)) then - total_transport = global_area_integral(fluxes%sens, G, scale=US%QRZ_T_to_W_m2) + total_transport = global_area_integral(fluxes%sens, G, unscale=US%QRZ_T_to_W_m2) call post_data(handles%id_total_sens, total_transport, diag) endif if ((handles%id_sens_ga > 0) .and. associated(fluxes%sens)) then @@ -3103,7 +3103,7 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h endif if ((handles%id_total_heat_added > 0) .and. associated(fluxes%heat_added)) then - total_transport = global_area_integral(fluxes%heat_added, G, scale=US%QRZ_T_to_W_m2) + total_transport = global_area_integral(fluxes%heat_added, G, unscale=US%QRZ_T_to_W_m2) call post_data(handles%id_total_heat_added, total_transport, diag) endif @@ -3113,21 +3113,21 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h if ((handles%id_saltflux > 0) .and. associated(fluxes%salt_flux)) & call post_data(handles%id_saltflux, fluxes%salt_flux, diag) if ((handles%id_total_saltflux > 0) .and. associated(fluxes%salt_flux)) then - total_transport = ppt2mks*global_area_integral(fluxes%salt_flux, G, scale=US%RZ_T_to_kg_m2s) + total_transport = ppt2mks*global_area_integral(fluxes%salt_flux, G, unscale=US%RZ_T_to_kg_m2s) call post_data(handles%id_total_saltflux, total_transport, diag) endif if ((handles%id_saltFluxAdded > 0) .and. associated(fluxes%salt_flux_added)) & call post_data(handles%id_saltFluxAdded, fluxes%salt_flux_added, diag) if ((handles%id_total_saltFluxAdded > 0) .and. associated(fluxes%salt_flux_added)) then - total_transport = ppt2mks*global_area_integral(fluxes%salt_flux_added, G, scale=US%RZ_T_to_kg_m2s) + total_transport = ppt2mks*global_area_integral(fluxes%salt_flux_added, G, unscale=US%RZ_T_to_kg_m2s) call post_data(handles%id_total_saltFluxAdded, total_transport, diag) endif if (handles%id_saltFluxIn > 0 .and. associated(fluxes%salt_flux_in)) & call post_data(handles%id_saltFluxIn, fluxes%salt_flux_in, diag) if ((handles%id_total_saltFluxIn > 0) .and. associated(fluxes%salt_flux_in)) then - total_transport = ppt2mks*global_area_integral(fluxes%salt_flux_in, G, scale=US%RZ_T_to_kg_m2s) + total_transport = ppt2mks*global_area_integral(fluxes%salt_flux_in, G, unscale=US%RZ_T_to_kg_m2s) call post_data(handles%id_total_saltFluxIn, total_transport, diag) endif diff --git a/src/diagnostics/MOM_spatial_means.F90 b/src/diagnostics/MOM_spatial_means.F90 index 0d656edf6d..c5def4fb28 100644 --- a/src/diagnostics/MOM_spatial_means.F90 +++ b/src/diagnostics/MOM_spatial_means.F90 @@ -32,7 +32,7 @@ module MOM_spatial_means contains !> Return the global area mean of a variable. This uses reproducing sums. -function global_area_mean(var, G, scale, tmp_scale) +function global_area_mean(var, G, scale, tmp_scale, unscale) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure real, dimension(SZI_(G),SZJ_(G)), intent(in) :: var !< The variable to average in !! arbitrary, possibly rescaled units [A ~> a] @@ -41,6 +41,11 @@ function global_area_mean(var, G, scale, tmp_scale) !! units to enable the use of the reproducing sums real, optional, intent(in) :: tmp_scale !< A temporary rescaling factor for the variable !! that is reversed in the return value [a A-1 ~> 1] + real, optional, intent(in) :: unscale !< A rescaling factor for the variable [a A-1 ~> 1] + !! that converts it back to unscaled (e.g., mks) + !! units to enable the use of the reproducing sums. + !! Here scale and unscale are synonymous, but unscale + !! is preferred and takes precedence if both are present. real :: global_area_mean ! The mean of the variable in arbitrary unscaled units [a] or scaled units [A ~> a] ! Local variables @@ -53,7 +58,9 @@ function global_area_mean(var, G, scale, tmp_scale) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec temp_scale = 1.0 ; if (present(tmp_scale)) temp_scale = tmp_scale - scalefac = G%US%L_to_m**2*temp_scale ; if (present(scale)) scalefac = scalefac * scale + scalefac = G%US%L_to_m**2*temp_scale + if (present(unscale)) then ; scalefac = scalefac * unscale + elseif (present(scale)) then ; scalefac = scalefac * scale ; endif tmpForSumming(:,:) = 0. do j=js,je ; do i=is,ie @@ -148,7 +155,7 @@ end function global_area_mean_u !> Return the global area integral of a variable, by default using the masked area from the !! grid, but an alternate could be used instead. This uses reproducing sums. -function global_area_integral(var, G, scale, area, tmp_scale) +function global_area_integral(var, G, scale, area, tmp_scale, unscale) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure real, dimension(SZI_(G),SZJ_(G)), intent(in) :: var !< The variable to integrate in !! arbitrary, possibly rescaled units [A ~> a] @@ -159,6 +166,11 @@ function global_area_integral(var, G, scale, area, tmp_scale) !! any required masking [L2 ~> m2]. real, optional, intent(in) :: tmp_scale !< A temporary rescaling factor for the !! variable that is reversed in the return value [a A-1 ~> 1] + real, optional, intent(in) :: unscale !< A rescaling factor for the variable [a A-1 ~> 1] + !! that converts it back to unscaled (e.g., mks) + !! units to enable the use of the reproducing sums. + !! Here scale and unscale are synonymous, but unscale + !! is preferred and takes precedence if both are present. real :: global_area_integral !< The returned area integral, usually in the units of var times an area, !! [a m2] or [A m2 ~> a m2] depending on which optional arguments are provided @@ -172,7 +184,9 @@ function global_area_integral(var, G, scale, area, tmp_scale) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec temp_scale = 1.0 ; if (present(tmp_scale)) temp_scale = tmp_scale - scalefac = G%US%L_to_m**2*temp_scale ; if (present(scale)) scalefac = scalefac * scale + scalefac = G%US%L_to_m**2*temp_scale + if (present(unscale)) then ; scalefac = scalefac * unscale + elseif (present(scale)) then ; scalefac = scalefac * scale ; endif tmpForSumming(:,:) = 0. if (present(area)) then @@ -193,7 +207,7 @@ function global_area_integral(var, G, scale, area, tmp_scale) end function global_area_integral !> Return the layerwise global thickness-weighted mean of a variable. This uses reproducing sums. -function global_layer_mean(var, h, G, GV, scale, tmp_scale) +function global_layer_mean(var, h, G, GV, scale, tmp_scale, unscale) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: var !< The variable to average in @@ -204,6 +218,11 @@ function global_layer_mean(var, h, G, GV, scale, tmp_scale) !! units to enable the use of the reproducing sums real, optional, intent(in) :: tmp_scale !< A temporary rescaling factor for the !! variable that is reversed in the return value [a A-1 ~> 1] + real, optional, intent(in) :: unscale !< A rescaling factor for the variable [a A-1 ~> 1] + !! that converts it back to unscaled (e.g., mks) + !! units to enable the use of the reproducing sums. + !! Here scale and unscale are synonymous, but unscale + !! is preferred and takes precedence. real, dimension(SZK_(GV)) :: global_layer_mean !< The mean of the variable in the arbitrary scaled [A] !! or unscaled [a] units of var, depending on which optional !! arguments are provided @@ -224,7 +243,9 @@ function global_layer_mean(var, h, G, GV, scale, tmp_scale) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke temp_scale = 1.0 ; if (present(tmp_scale)) temp_scale = tmp_scale - scalefac = temp_scale ; if (present(scale)) scalefac = scale * temp_scale + scalefac = temp_scale + if (present(unscale)) then ; scalefac = unscale * temp_scale + elseif (present(scale)) then ; scalefac = scale * temp_scale ; endif tmpForSumming(:,:,:) = 0. ; weight(:,:,:) = 0. do k=1,nz ; do j=js,je ; do i=is,ie @@ -243,7 +264,7 @@ function global_layer_mean(var, h, G, GV, scale, tmp_scale) end function global_layer_mean !> Find the global thickness-weighted mean of a variable. This uses reproducing sums. -function global_volume_mean(var, h, G, GV, scale, tmp_scale) +function global_volume_mean(var, h, G, GV, scale, tmp_scale, unscale) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & @@ -256,6 +277,11 @@ function global_volume_mean(var, h, G, GV, scale, tmp_scale) !! units to enable the use of the reproducing sums real, optional, intent(in) :: tmp_scale !< A temporary rescaling factor for the !! variable that is reversed in the return value [a A-1 ~> 1] + real, optional, intent(in) :: unscale !< A rescaling factor for the variable [a A-1 ~> 1] + !! that converts it back to unscaled (e.g., mks) + !! units to enable the use of the reproducing sums. + !! Here scale and unscale are synonymous, but unscale + !! is preferred and takes precedence if both are present. real :: global_volume_mean !< The thickness-weighted average of var in the arbitrary scaled [A] or !! unscaled [a] units of var, depending on which optional arguments are provided @@ -271,7 +297,9 @@ function global_volume_mean(var, h, G, GV, scale, tmp_scale) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke temp_scale = 1.0 ; if (present(tmp_scale)) temp_scale = tmp_scale - scalefac = temp_scale ; if (present(scale)) scalefac = temp_scale * scale + scalefac = temp_scale + if (present(unscale)) then ; scalefac = temp_scale * unscale + elseif (present(scale)) then ; scalefac = temp_scale * scale ; endif tmpForSumming(:,:) = 0. ; sum_weight(:,:) = 0. do k=1,nz ; do j=js,je ; do i=is,ie @@ -286,7 +314,7 @@ end function global_volume_mean !> Find the global mass-weighted integral of a variable. This uses reproducing sums. -function global_mass_integral(h, G, GV, var, on_PE_only, scale, tmp_scale) +function global_mass_integral(h, G, GV, var, on_PE_only, scale, tmp_scale, unscale) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & @@ -301,6 +329,11 @@ function global_mass_integral(h, G, GV, var, on_PE_only, scale, tmp_scale) !! units to enable the use of the reproducing sums real, optional, intent(in) :: tmp_scale !< A temporary rescaling factor for the !! variable that is reversed in the return value [a A-1 ~> 1] + real, optional, intent(in) :: unscale !< A rescaling factor for the variable [a A-1 ~> 1] + !! that converts it back to unscaled (e.g., mks) + !! units to enable the use of the reproducing sums. + !! Here scale and unscale are synonymous, but unscale + !! is preferred and takes precedence if both are present. real :: global_mass_integral !< The mass-weighted integral of var (or 1) in !! kg times the arbitrary units of var [kg a] or [kg A ~> kg a] @@ -315,7 +348,9 @@ function global_mass_integral(h, G, GV, var, on_PE_only, scale, tmp_scale) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke temp_scale = 1.0 ; if (present(tmp_scale)) temp_scale = tmp_scale - scalefac = G%US%L_to_m**2*temp_scale ; if (present(scale)) scalefac = scalefac * scale + scalefac = G%US%L_to_m**2*temp_scale + if (present(unscale)) then ; scalefac = scalefac * unscale + elseif (present(scale)) then ; scalefac = scalefac * scale ; endif tmpForSumming(:,:) = 0.0 if (present(var)) then @@ -346,7 +381,7 @@ end function global_mass_integral !> Find the global mass-weighted order invariant integral of a variable in mks units, !! returning the value as an EFP_type. This uses reproducing sums. -function global_mass_int_EFP(h, G, GV, var, on_PE_only, scale) +function global_mass_int_EFP(h, G, GV, var, on_PE_only, scale, unscale) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & @@ -359,6 +394,11 @@ function global_mass_int_EFP(h, G, GV, var, on_PE_only, scale) real, optional, intent(in) :: scale !< A rescaling factor for the variable [a A-1 ~> 1] !! that converts it back to unscaled (e.g., mks) !! units to enable the use of the reproducing sums + real, optional, intent(in) :: unscale !< A rescaling factor for the variable [a A-1 ~> 1] + !! that converts it back to unscaled (e.g., mks) + !! units to enable the use of the reproducing sums. + !! Here scale and unscale are synonymous, but unscale + !! is preferred and takes precedence if both are present. type(EFP_type) :: global_mass_int_EFP !< The mass-weighted integral of var (or 1) in !! kg times the arbitrary units of var [kg a] @@ -373,7 +413,8 @@ function global_mass_int_EFP(h, G, GV, var, on_PE_only, scale) isr = is - (G%isd-1) ; ier = ie - (G%isd-1) ; jsr = js - (G%jsd-1) ; jer = je - (G%jsd-1) scalefac = GV%H_to_kg_m2 * G%US%L_to_m**2 - if (present(scale)) scalefac = scale * scalefac + if (present(unscale)) then ; scalefac = unscale * scalefac + elseif (present(scale)) then ; scalefac = scale * scalefac ; endif tmpForSum(:,:) = 0.0 if (present(var)) then @@ -395,7 +436,7 @@ end function global_mass_int_EFP !> Determine the global mean of a field along rows of constant i, returning it !! in a 1-d array using the local indexing. This uses reproducing sums. -subroutine global_i_mean(array, i_mean, G, mask, scale, tmp_scale) +subroutine global_i_mean(array, i_mean, G, mask, scale, tmp_scale, unscale) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure real, dimension(SZI_(G),SZJ_(G)), intent(in) :: array !< The variable being averaged in !! arbitrary, possibly rescaled units [A ~> a] @@ -407,14 +448,19 @@ subroutine global_i_mean(array, i_mean, G, mask, scale, tmp_scale) !! units to enable the use of the reproducing sums real, optional, intent(in) :: tmp_scale !< A rescaling factor for the internal !! calculations that is removed from the output [a A-1 ~> 1] + real, optional, intent(in) :: unscale !< A rescaling factor for the variable [a A-1 ~> 1] + !! that converts it back to unscaled (e.g., mks) + !! units to enable the use of the reproducing sums. + !! Here scale and unscale are synonymous, but unscale + !! is preferred and takes precedence if both are present. ! Local variables ! In the following comments, [A] is used to indicate the arbitrary, possibly rescaled units of the ! input array while [a] indicates the unscaled (e.g., mks) units that can be used with the reproducing sums type(EFP_type), allocatable, dimension(:) :: asum ! The masked sum of the variable in each row [a] type(EFP_type), allocatable, dimension(:) :: mask_sum ! The sum of the mask values in each row [nondim] - real :: scalefac ! A scaling factor for the variable [a A-1 ~> 1] - real :: unscale ! A factor for undoing any internal rescaling before output [A a-1 ~> 1] + real :: scalefac ! A scaling factor for the variable [a A-1 ~> 1] + real :: rescale ! A factor for redoing any internal rescaling before output [A a-1 ~> 1] real :: mask_sum_r ! The sum of the mask values in a row [nondim] integer :: is, ie, js, je, idg_off, jdg_off integer :: i, j @@ -422,10 +468,13 @@ subroutine global_i_mean(array, i_mean, G, mask, scale, tmp_scale) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec idg_off = G%idg_offset ; jdg_off = G%jdg_offset - scalefac = 1.0 ; if (present(scale)) scalefac = scale - unscale = 1.0 + scalefac = 1.0 + if (present(unscale)) then ; scalefac = unscale + elseif (present(scale)) then ; scalefac = scale ; endif + + rescale = 1.0 if (present(tmp_scale)) then ; if (tmp_scale /= 0.0) then - scalefac = scalefac * tmp_scale ; unscale = 1.0 / tmp_scale + scalefac = scalefac * tmp_scale ; rescale = 1.0 / tmp_scale endif ; endif call reset_EFP_overflow_error() @@ -479,7 +528,7 @@ subroutine global_i_mean(array, i_mean, G, mask, scale, tmp_scale) enddo endif - if (unscale /= 1.0) then ; do j=js,je ; i_mean(j) = unscale*i_mean(j) ; enddo ; endif + if (rescale /= 1.0) then ; do j=js,je ; i_mean(j) = rescale*i_mean(j) ; enddo ; endif deallocate(asum) @@ -487,7 +536,7 @@ end subroutine global_i_mean !> Determine the global mean of a field along rows of constant j, returning it !! in a 1-d array using the local indexing. This uses reproducing sums. -subroutine global_j_mean(array, j_mean, G, mask, scale, tmp_scale) +subroutine global_j_mean(array, j_mean, G, mask, scale, tmp_scale, unscale) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure real, dimension(SZI_(G),SZJ_(G)), intent(in) :: array !< The variable being averaged in !! arbitrary, possibly rescaled units [A ~> a] @@ -499,6 +548,11 @@ subroutine global_j_mean(array, j_mean, G, mask, scale, tmp_scale) !! units to enable the use of the reproducing sums real, optional, intent(in) :: tmp_scale !< A rescaling factor for the internal !! calculations that is removed from the output [a A-1 ~> 1] + real, optional, intent(in) :: unscale !< A rescaling factor for the variable [a A-1 ~> 1] + !! that converts it back to unscaled (e.g., mks) + !! units to enable the use of the reproducing sums. + !! Here scale and unscale are synonymous, but unscale + !! is preferred and takes precedence if both are present. ! Local variables ! In the following comments, [A] is used to indicate the arbitrary, possibly rescaled units of the @@ -506,18 +560,21 @@ subroutine global_j_mean(array, j_mean, G, mask, scale, tmp_scale) type(EFP_type), allocatable, dimension(:) :: asum ! The masked sum of the variable in each row [a] type(EFP_type), allocatable, dimension(:) :: mask_sum ! The sum of the mask values in each row [nondim] real :: mask_sum_r ! The sum of the mask values in a row [nondim] - real :: scalefac ! A scaling factor for the variable [a A-1 ~> 1] - real :: unscale ! A factor for undoing any internal rescaling before output [A a-1 ~> 1] + real :: scalefac ! A scaling factor for the variable [a A-1 ~> 1] + real :: rescale ! A factor for redoing any internal rescaling before output [A a-1 ~> 1] integer :: is, ie, js, je, idg_off, jdg_off integer :: i, j is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec idg_off = G%idg_offset ; jdg_off = G%jdg_offset - scalefac = 1.0 ; if (present(scale)) scalefac = scale - unscale = 1.0 + scalefac = 1.0 + if (present(unscale)) then ; scalefac = unscale + elseif (present(scale)) then ; scalefac = scale ; endif + + rescale = 1.0 if (present(tmp_scale)) then ; if (tmp_scale /= 0.0) then - scalefac = scalefac * tmp_scale ; unscale = 1.0 / tmp_scale + scalefac = scalefac * tmp_scale ; rescale = 1.0 / tmp_scale endif ; endif call reset_EFP_overflow_error() @@ -571,14 +628,14 @@ subroutine global_j_mean(array, j_mean, G, mask, scale, tmp_scale) enddo endif - if (unscale /= 1.0) then ; do i=is,ie ; j_mean(i) = unscale*j_mean(i) ; enddo ; endif + if (rescale /= 1.0) then ; do i=is,ie ; j_mean(i) = rescale*j_mean(i) ; enddo ; endif deallocate(asum) end subroutine global_j_mean !> Adjust 2d array such that area mean is zero without moving the zero contour -subroutine adjust_area_mean_to_zero(array, G, scaling, unit_scale) +subroutine adjust_area_mean_to_zero(array, G, scaling, unit_scale, unscale) type(ocean_grid_type), intent(in) :: G !< Grid structure real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: array !< 2D array to be adjusted in !! arbitrary, possibly rescaled units [A ~> a] @@ -586,6 +643,11 @@ subroutine adjust_area_mean_to_zero(array, G, scaling, unit_scale) real, optional, intent(in) :: unit_scale !< A rescaling factor for the variable [a A-1 ~> 1] !! that converts it back to unscaled (e.g., mks) !! units to enable the use of the reproducing sums + real, optional, intent(in) :: unscale !< A rescaling factor for the variable [a A-1 ~> 1] + !! that converts it back to unscaled (e.g., mks) + !! units to enable the use of the reproducing sums. + !! Here unit_scale and unscale are synonymous, but unscale + !! is preferred and takes precedence if both are present. ! Local variables ! In the following comments, [A] is used to indicate the arbitrary, possibly rescaled units of the ! input array while [a] indicates the unscaled (e.g., mks) units that can be used with the reproducing sums @@ -598,7 +660,10 @@ subroutine adjust_area_mean_to_zero(array, G, scaling, unit_scale) real :: posScale, negScale ! The scaling factor to apply to positive or negative values [nondim] integer :: i,j - scalefac = 1.0 ; if (present(unit_scale)) scalefac = unit_scale + scalefac = 1.0 + if (present(unscale)) then ; scalefac = unscale + elseif (present(unit_scale)) then ; scalefac = unit_scale ; endif + I_scalefac = 0.0 ; if (scalefac /= 0.0) I_scalefac = 1.0 / scalefac ! areaXposVals(:,:) = 0. ! This zeros out halo points. diff --git a/src/framework/MOM_checksums.F90 b/src/framework/MOM_checksums.F90 index 00e4ba4918..4e349e0fb7 100644 --- a/src/framework/MOM_checksums.F90 +++ b/src/framework/MOM_checksums.F90 @@ -102,13 +102,17 @@ module MOM_checksums contains !> Checksum a scalar field (consistent with array checksums) -subroutine chksum0(scalar, mesg, scale, logunit) +subroutine chksum0(scalar, mesg, scale, logunit, unscale) real, intent(in) :: scalar !< The array to be checksummed in !! arbitrary, possibly rescaled units [A ~> a] character(len=*), intent(in) :: mesg !< An identifying message real, optional, intent(in) :: scale !< A factor to convert this array back to unscaled units !! for checksums and output [a A-1 ~> 1] integer, optional, intent(in) :: logunit !< IO unit for checksum logging + real, optional, intent(in) :: unscale !< A factor to convert this array back to unscaled units + !! for checksums and output [a A-1 ~> 1]. + !! Here scale and unscale are synonymous, but unscale + !! takes precedence if both are present. ! Local variables ! In the following comments, [A] is used to indicate the arbitrary, possibly rescaled units @@ -122,7 +126,10 @@ subroutine chksum0(scalar, mesg, scale, logunit) if (checkForNaNs .and. is_NaN(scalar)) & call chksum_error(FATAL, 'NaN detected: '//trim(mesg)) - scaling = 1.0 ; if (present(scale)) scaling = scale + scaling = 1.0 + if (present(unscale)) then ; scaling = unscale + elseif (present(scale)) then ; scaling = scale ; endif + iounit = error_unit ; if (present(logunit)) iounit = logunit if (calculateStatistics) then @@ -141,13 +148,17 @@ end subroutine chksum0 !> Checksum a 1d array (typically a column). -subroutine zchksum(array, mesg, scale, logunit) +subroutine zchksum(array, mesg, scale, logunit, unscale) real, dimension(:), intent(in) :: array !< The array to be checksummed in !! arbitrary, possibly rescaled units [A ~> a] character(len=*), intent(in) :: mesg !< An identifying message real, optional, intent(in) :: scale !< A factor to convert this array back to unscaled units !! for checksums and output [a A-1 ~> 1] integer, optional, intent(in) :: logunit !< IO unit for checksum logging + real, optional, intent(in) :: unscale !< A factor to convert this array back to unscaled units + !! for checksums and output [a A-1 ~> 1]. + !! Here scale and unscale are synonymous, but unscale + !! takes precedence if both are present. ! Local variables ! In the following comments, [A] is used to indicate the arbitrary, possibly rescaled units @@ -165,14 +176,17 @@ subroutine zchksum(array, mesg, scale, logunit) call chksum_error(FATAL, 'NaN detected: '//trim(mesg)) endif - scaling = 1.0 ; if (present(scale)) scaling = scale + scaling = 1.0 + if (present(unscale)) then ; scaling = unscale + elseif (present(scale)) then ; scaling = scale ; endif + iounit = error_unit ; if (present(logunit)) iounit = logunit if (calculateStatistics) then - if (present(scale)) then + if (present(unscale) .or. present(scale)) then allocate(rescaled_array(LBOUND(array,1):UBOUND(array,1)), source=0.0) do k=1, size(array, 1) - rescaled_array(k) = scale * array(k) + rescaled_array(k) = scaling * array(k) enddo call subStats(rescaled_array, aMean, aMin, aMax) @@ -192,15 +206,15 @@ subroutine zchksum(array, mesg, scale, logunit) contains - integer function subchk(array, scale) + integer function subchk(array, unscale) real, dimension(:), intent(in) :: array !< The array to be checksummed in !! arbitrary, possibly rescaled units [A ~> a] - real, intent(in) :: scale !< A factor to convert this array back to unscaled units - !! for checksums and output [a A-1 ~> 1] + real, intent(in) :: unscale !< A factor to convert this array back to unscaled units + !! for checksums and output [a A-1 ~> 1] integer :: k, bc subchk = 0 do k=LBOUND(array, 1), UBOUND(array, 1) - bc = bitcount(abs(scale * array(k))) + bc = bitcount(abs(unscale * array(k))) subchk = subchk + bc enddo subchk=mod(subchk, bc_modulus) @@ -228,7 +242,7 @@ end subroutine zchksum !> Checksums on a pair of 2d arrays staggered at tracer points. subroutine chksum_pair_h_2d(mesg, arrayA, arrayB, HI, haloshift, omit_corners, & - scale, logunit, scalar_pair) + scale, logunit, scalar_pair, unscale) character(len=*), intent(in) :: mesg !< Identifying messages type(hor_index_type), target, intent(in) :: HI !< A horizontal index type real, dimension(HI%isd:,HI%jsd:), target, intent(in) :: arrayA !< The first array to be checksummed in @@ -242,6 +256,10 @@ subroutine chksum_pair_h_2d(mesg, arrayA, arrayB, HI, haloshift, omit_corners, & integer, optional, intent(in) :: logunit !< IO unit for checksum logging logical, optional, intent(in) :: scalar_pair !< If true, then the arrays describe !! a scalar, rather than vector + real, optional, intent(in) :: unscale !< A factor to convert this array back to unscaled units + !! for checksums and output [a A-1 ~> 1]. + !! Here scale and unscale are synonymous, but unscale + !! takes precedence if both are present. logical :: vector_pair integer :: turns type(hor_index_type), pointer :: HI_in @@ -271,18 +289,18 @@ subroutine chksum_pair_h_2d(mesg, arrayA, arrayB, HI, haloshift, omit_corners, & if (present(haloshift)) then call chksum_h_2d(arrayA_in, 'x '//mesg, HI_in, haloshift, omit_corners, & - scale=scale, logunit=logunit) + scale=scale, logunit=logunit, unscale=unscale) call chksum_h_2d(arrayB_in, 'y '//mesg, HI_in, haloshift, omit_corners, & - scale=scale, logunit=logunit) + scale=scale, logunit=logunit, unscale=unscale) else - call chksum_h_2d(arrayA_in, 'x '//mesg, HI_in, scale=scale, logunit=logunit) - call chksum_h_2d(arrayB_in, 'y '//mesg, HI_in, scale=scale, logunit=logunit) + call chksum_h_2d(arrayA_in, 'x '//mesg, HI_in, scale=scale, logunit=logunit, unscale=unscale) + call chksum_h_2d(arrayB_in, 'y '//mesg, HI_in, scale=scale, logunit=logunit, unscale=unscale) endif end subroutine chksum_pair_h_2d !> Checksums on a pair of 3d arrays staggered at tracer points. subroutine chksum_pair_h_3d(mesg, arrayA, arrayB, HI, haloshift, omit_corners, & - scale, logunit, scalar_pair) + scale, logunit, scalar_pair, unscale) character(len=*), intent(in) :: mesg !< Identifying messages type(hor_index_type), target, intent(in) :: HI !< A horizontal index type real, dimension(HI%isd:,HI%jsd:, :), target, intent(in) :: arrayA !< The first array to be checksummed in @@ -296,6 +314,11 @@ subroutine chksum_pair_h_3d(mesg, arrayA, arrayB, HI, haloshift, omit_corners, & integer, optional, intent(in) :: logunit !< IO unit for checksum logging logical, optional, intent(in) :: scalar_pair !< If true, then the arrays describe !! a scalar, rather than vector + real, optional, intent(in) :: unscale !< A factor to convert this array back to unscaled units + !! for checksums and output [a A-1 ~> 1]. + !! Here scale and unscale are synonymous, but unscale + !! takes precedence if both are present. + ! Local variables logical :: vector_pair integer :: turns type(hor_index_type), pointer :: HI_in @@ -325,19 +348,19 @@ subroutine chksum_pair_h_3d(mesg, arrayA, arrayB, HI, haloshift, omit_corners, & if (present(haloshift)) then call chksum_h_3d(arrayA_in, 'x '//mesg, HI_in, haloshift, omit_corners, & - scale=scale, logunit=logunit) + scale=scale, logunit=logunit, unscale=unscale) call chksum_h_3d(arrayB_in, 'y '//mesg, HI_in, haloshift, omit_corners, & - scale=scale, logunit=logunit) + scale=scale, logunit=logunit, unscale=unscale) else - call chksum_h_3d(arrayA_in, 'x '//mesg, HI_in, scale=scale, logunit=logunit) - call chksum_h_3d(arrayB_in, 'y '//mesg, HI_in, scale=scale, logunit=logunit) + call chksum_h_3d(arrayA_in, 'x '//mesg, HI_in, scale=scale, logunit=logunit, unscale=unscale) + call chksum_h_3d(arrayB_in, 'y '//mesg, HI_in, scale=scale, logunit=logunit, unscale=unscale) endif ! NOTE: automatic deallocation of array[AB]_in end subroutine chksum_pair_h_3d !> Checksums a 2d array staggered at tracer points. -subroutine chksum_h_2d(array_m, mesg, HI_m, haloshift, omit_corners, scale, logunit) +subroutine chksum_h_2d(array_m, mesg, HI_m, haloshift, omit_corners, scale, logunit, unscale) type(hor_index_type), target, intent(in) :: HI_m !< Horizontal index bounds of the model grid real, dimension(HI_m%isd:,HI_m%jsd:), target, intent(in) :: array_m !< Field array on the model grid in !! arbitrary, possibly rescaled units [A ~> a] @@ -347,6 +370,10 @@ subroutine chksum_h_2d(array_m, mesg, HI_m, haloshift, omit_corners, scale, logu real, optional, intent(in) :: scale !< A factor to convert this array back to unscaled units !! for checksums and output [a A-1 ~> 1] integer, optional, intent(in) :: logunit !< IO unit for checksum logging + real, optional, intent(in) :: unscale !< A factor to convert this array back to unscaled units + !! for checksums and output [a A-1 ~> 1]. + !! Here scale and unscale are synonymous, but unscale + !! takes precedence if both are present. ! Local variables ! In the following comments, [A] is used to indicate the arbitrary, possibly rescaled units @@ -383,15 +410,18 @@ subroutine chksum_h_2d(array_m, mesg, HI_m, haloshift, omit_corners, scale, logu ! call chksum_error(FATAL, 'NaN detected in halo: '//trim(mesg)) endif - scaling = 1.0 ; if (present(scale)) scaling = scale + scaling = 1.0 + if (present(unscale)) then ; scaling = unscale + elseif (present(scale)) then ; scaling = scale ; endif + iounit = error_unit ; if (present(logunit)) iounit = logunit if (calculateStatistics) then - if (present(scale)) then + if (present(unscale) .or. present(scale)) then allocate( rescaled_array(LBOUND(array,1):UBOUND(array,1), & LBOUND(array,2):UBOUND(array,2)), source=0.0 ) do j=HI%jsc,HI%jec ; do i=HI%isc,HI%iec - rescaled_array(i,j) = scale*array(i,j) + rescaled_array(i,j) = scaling*array(i,j) enddo ; enddo call subStats(HI, rescaled_array, aMean, aMin, aMax) deallocate(rescaled_array) @@ -445,18 +475,18 @@ subroutine chksum_h_2d(array_m, mesg, HI_m, haloshift, omit_corners, scale, logu endif contains - integer function subchk(array, HI, di, dj, scale) + integer function subchk(array, HI, di, dj, unscale) type(hor_index_type), intent(in) :: HI !< A horizontal index type real, dimension(HI%isd:,HI%jsd:), intent(in) :: array !< The array to be checksummed in !! arbitrary, possibly rescaled units [A ~> a] integer, intent(in) :: di !< i- direction array shift for this checksum integer, intent(in) :: dj !< j- direction array shift for this checksum - real, intent(in) :: scale !< A factor to convert this array back to unscaled units + real, intent(in) :: unscale !< A factor to convert this array back to unscaled units !! for checksums and output [a A-1 ~> 1] integer :: i, j, bc subchk = 0 do j=HI%jsc+dj,HI%jec+dj ; do i=HI%isc+di,HI%iec+di - bc = bitcount(abs(scale*array(i,j))) + bc = bitcount(abs(unscale*array(i,j))) subchk = subchk + bc enddo ; enddo call sum_across_PEs(subchk) @@ -491,7 +521,7 @@ end subroutine chksum_h_2d !> Checksums on a pair of 2d arrays staggered at q-points. subroutine chksum_pair_B_2d(mesg, arrayA, arrayB, HI, haloshift, symmetric, & - omit_corners, scale, logunit, scalar_pair) + omit_corners, scale, logunit, scalar_pair, unscale) character(len=*), intent(in) :: mesg !< Identifying messages type(hor_index_type), target, intent(in) :: HI !< A horizontal index type real, dimension(HI%isd:,HI%jsd:), target, intent(in) :: arrayA !< The first array to be checksummed in @@ -507,6 +537,10 @@ subroutine chksum_pair_B_2d(mesg, arrayA, arrayB, HI, haloshift, symmetric, & integer, optional, intent(in) :: logunit !< IO unit for checksum logging logical, optional, intent(in) :: scalar_pair !< If true, then the arrays describe !! a scalar, rather than vector + real, optional, intent(in) :: unscale !< A factor to convert this array back to unscaled units + !! for checksums and output [a A-1 ~> 1]. + !! Here scale and unscale are synonymous, but unscale + !! takes precedence if both are present. logical :: sym logical :: vector_pair @@ -540,21 +574,21 @@ subroutine chksum_pair_B_2d(mesg, arrayA, arrayB, HI, haloshift, symmetric, & if (present(haloshift)) then call chksum_B_2d(arrayA_in, 'x '//mesg, HI_in, haloshift, symmetric=sym, & - omit_corners=omit_corners, scale=scale, logunit=logunit) + omit_corners=omit_corners, scale=scale, logunit=logunit, unscale=unscale) call chksum_B_2d(arrayB_in, 'y '//mesg, HI_in, haloshift, symmetric=sym, & - omit_corners=omit_corners, scale=scale, logunit=logunit) + omit_corners=omit_corners, scale=scale, logunit=logunit, unscale=unscale) else - call chksum_B_2d(arrayA_in, 'x '//mesg, HI_in, symmetric=sym, scale=scale, & - logunit=logunit) - call chksum_B_2d(arrayB_in, 'y '//mesg, HI_in, symmetric=sym, scale=scale, & - logunit=logunit) + call chksum_B_2d(arrayA_in, 'x '//mesg, HI_in, symmetric=sym, & + scale=scale, logunit=logunit, unscale=unscale) + call chksum_B_2d(arrayB_in, 'y '//mesg, HI_in, symmetric=sym, & + scale=scale, logunit=logunit, unscale=unscale) endif end subroutine chksum_pair_B_2d !> Checksums on a pair of 3d arrays staggered at q-points. subroutine chksum_pair_B_3d(mesg, arrayA, arrayB, HI, haloshift, symmetric, & - omit_corners, scale, logunit, scalar_pair) + omit_corners, scale, logunit, scalar_pair, unscale) character(len=*), intent(in) :: mesg !< Identifying messages type(hor_index_type), target, intent(in) :: HI !< A horizontal index type real, dimension(HI%IsdB:,HI%JsdB:, :), target, intent(in) :: arrayA !< The first array to be checksummed in @@ -570,7 +604,11 @@ subroutine chksum_pair_B_3d(mesg, arrayA, arrayB, HI, haloshift, symmetric, & integer, optional, intent(in) :: logunit !< IO unit for checksum logging logical, optional, intent(in) :: scalar_pair !< If true, then the arrays describe !! a scalar, rather than vector - + real, optional, intent(in) :: unscale !< A factor to convert this array back to unscaled units + !! for checksums and output [a A-1 ~> 1]. + !! Here scale and unscale are synonymous, but unscale + !! takes precedence if both are present. + ! Local variables logical :: vector_pair integer :: turns type(hor_index_type), pointer :: HI_in @@ -600,20 +638,20 @@ subroutine chksum_pair_B_3d(mesg, arrayA, arrayB, HI, haloshift, symmetric, & if (present(haloshift)) then call chksum_B_3d(arrayA_in, 'x '//mesg, HI_in, haloshift, symmetric, & - omit_corners, scale=scale, logunit=logunit) + omit_corners, scale=scale, logunit=logunit, unscale=unscale) call chksum_B_3d(arrayB_in, 'y '//mesg, HI_in, haloshift, symmetric, & - omit_corners, scale=scale, logunit=logunit) + omit_corners, scale=scale, logunit=logunit, unscale=unscale) else - call chksum_B_3d(arrayA_in, 'x '//mesg, HI_in, symmetric=symmetric, scale=scale, & - logunit=logunit) - call chksum_B_3d(arrayB_in, 'y '//mesg, HI_in, symmetric=symmetric, scale=scale, & - logunit=logunit) + call chksum_B_3d(arrayA_in, 'x '//mesg, HI_in, symmetric=symmetric, & + scale=scale, logunit=logunit, unscale=unscale) + call chksum_B_3d(arrayB_in, 'y '//mesg, HI_in, symmetric=symmetric, & + scale=scale, logunit=logunit, unscale=unscale) endif end subroutine chksum_pair_B_3d !> Checksums a 2d array staggered at corner points. subroutine chksum_B_2d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners, & - scale, logunit) + scale, logunit, unscale) type(hor_index_type), target, intent(in) :: HI_m !< A horizontal index type real, dimension(HI_m%IsdB:,HI_m%JsdB:), & target, intent(in) :: array_m !< The array to be checksummed in @@ -626,6 +664,10 @@ subroutine chksum_B_2d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners, real, optional, intent(in) :: scale !< A factor to convert this array back to unscaled units !! for checksums and output [a A-1 ~> 1] integer, optional, intent(in) :: logunit !< IO unit for checksum logging + real, optional, intent(in) :: unscale !< A factor to convert this array back to unscaled units + !! for checksums and output [a A-1 ~> 1]. + !! Here scale and unscale are synonymous, but unscale + !! takes precedence if both are present. ! Local variables ! In the following comments, [A] is used to indicate the arbitrary, possibly rescaled units @@ -662,19 +704,22 @@ subroutine chksum_B_2d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners, ! call chksum_error(FATAL, 'NaN detected in halo: '//trim(mesg)) endif - scaling = 1.0 ; if (present(scale)) scaling = scale + scaling = 1.0 + if (present(unscale)) then ; scaling = unscale + elseif (present(scale)) then ; scaling = scale ; endif + iounit = error_unit ; if (present(logunit)) iounit = logunit sym_stats = .false. ; if (present(symmetric)) sym_stats = symmetric if (present(haloshift)) then ; if (haloshift > 0) sym_stats = .true. ; endif if (calculateStatistics) then - if (present(scale)) then + if (present(unscale) .or. present(scale)) then allocate( rescaled_array(LBOUND(array,1):UBOUND(array,1), & LBOUND(array,2):UBOUND(array,2)), source=0.0 ) Is = HI%isc ; if (sym_stats) Is = HI%isc-1 Js = HI%jsc ; if (sym_stats) Js = HI%jsc-1 do J=Js,HI%JecB ; do I=Is,HI%IecB - rescaled_array(I,J) = scale*array(I,J) + rescaled_array(I,J) = scaling*array(I,J) enddo ; enddo call subStats(HI, rescaled_array, sym_stats, aMean, aMin, aMax) deallocate(rescaled_array) @@ -736,19 +781,19 @@ subroutine chksum_B_2d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners, contains - integer function subchk(array, HI, di, dj, scale) + integer function subchk(array, HI, di, dj, unscale) type(hor_index_type), intent(in) :: HI !< A horizontal index type real, dimension(HI%IsdB:,HI%JsdB:), intent(in) :: array !< The array to be checksummed in !! arbitrary, possibly rescaled units [A ~> a] integer, intent(in) :: di !< i- direction array shift for this checksum integer, intent(in) :: dj !< j- direction array shift for this checksum - real, intent(in) :: scale !< A factor to convert this array back to unscaled units + real, intent(in) :: unscale !< A factor to convert this array back to unscaled units !! for checksums and output [a A-1 ~> 1] integer :: i, j, bc subchk = 0 ! This line deliberately uses the h-point computational domain. do J=HI%jsc+dj,HI%jec+dj ; do I=HI%isc+di,HI%iec+di - bc = bitcount(abs(scale*array(I,J))) + bc = bitcount(abs(unscale*array(I,J))) subchk = subchk + bc enddo ; enddo call sum_across_PEs(subchk) @@ -787,7 +832,7 @@ end subroutine chksum_B_2d !> Checksums a pair of 2d velocity arrays staggered at C-grid locations subroutine chksum_uv_2d(mesg, arrayU, arrayV, HI, haloshift, symmetric, & - omit_corners, scale, logunit, scalar_pair) + omit_corners, scale, logunit, scalar_pair, unscale) character(len=*), intent(in) :: mesg !< Identifying messages type(hor_index_type), target, intent(in) :: HI !< A horizontal index type real, dimension(HI%IsdB:,HI%jsd:), target, intent(in) :: arrayU !< The u-component array to be checksummed in @@ -803,6 +848,11 @@ subroutine chksum_uv_2d(mesg, arrayU, arrayV, HI, haloshift, symmetric, & integer, optional, intent(in) :: logunit !< IO unit for checksum logging logical, optional, intent(in) :: scalar_pair !< If true, then the arrays describe a !! a scalar, rather than vector + real, optional, intent(in) :: unscale !< A factor to convert this array back to unscaled units + !! for checksums and output [a A-1 ~> 1]. + !! Here scale and unscale are synonymous, but unscale + !! takes precedence if both are present. + ! Local variables logical :: vector_pair integer :: turns type(hor_index_type), pointer :: HI_in @@ -832,20 +882,20 @@ subroutine chksum_uv_2d(mesg, arrayU, arrayV, HI, haloshift, symmetric, & if (present(haloshift)) then call chksum_u_2d(arrayU_in, 'u '//mesg, HI_in, haloshift, symmetric, & - omit_corners, scale=scale, logunit=logunit) + omit_corners, scale=scale, logunit=logunit, unscale=unscale) call chksum_v_2d(arrayV_in, 'v '//mesg, HI_in, haloshift, symmetric, & - omit_corners, scale=scale, logunit=logunit) + omit_corners, scale=scale, logunit=logunit, unscale=unscale) else call chksum_u_2d(arrayU_in, 'u '//mesg, HI_in, symmetric=symmetric, & - scale=scale, logunit=logunit) + scale=scale, logunit=logunit, unscale=unscale) call chksum_v_2d(arrayV_in, 'v '//mesg, HI_in, symmetric=symmetric, & - scale=scale, logunit=logunit) + scale=scale, logunit=logunit, unscale=unscale) endif end subroutine chksum_uv_2d !> Checksums a pair of 3d velocity arrays staggered at C-grid locations subroutine chksum_uv_3d(mesg, arrayU, arrayV, HI, haloshift, symmetric, & - omit_corners, scale, logunit, scalar_pair) + omit_corners, scale, logunit, scalar_pair, unscale) character(len=*), intent(in) :: mesg !< Identifying messages type(hor_index_type), target, intent(in) :: HI !< A horizontal index type real, dimension(HI%IsdB:,HI%jsd:,:), target, intent(in) :: arrayU !< The u-component array to be checksummed in @@ -861,6 +911,11 @@ subroutine chksum_uv_3d(mesg, arrayU, arrayV, HI, haloshift, symmetric, & integer, optional, intent(in) :: logunit !< IO unit for checksum logging logical, optional, intent(in) :: scalar_pair !< If true, then the arrays describe a !! a scalar, rather than vector + real, optional, intent(in) :: unscale !< A factor to convert this array back to unscaled units + !! for checksums and output [a A-1 ~> 1]. + !! Here scale and unscale are synonymous, but unscale + !! takes precedence if both are present. + ! Local variables logical :: vector_pair integer :: turns type(hor_index_type), pointer :: HI_in @@ -890,20 +945,20 @@ subroutine chksum_uv_3d(mesg, arrayU, arrayV, HI, haloshift, symmetric, & if (present(haloshift)) then call chksum_u_3d(arrayU_in, 'u '//mesg, HI_in, haloshift, symmetric, & - omit_corners, scale=scale, logunit=logunit) + omit_corners, scale=scale, logunit=logunit, unscale=unscale) call chksum_v_3d(arrayV_in, 'v '//mesg, HI_in, haloshift, symmetric, & - omit_corners, scale=scale, logunit=logunit) + omit_corners, scale=scale, logunit=logunit, unscale=unscale) else call chksum_u_3d(arrayU_in, 'u '//mesg, HI_in, symmetric=symmetric, & - scale=scale, logunit=logunit) + scale=scale, logunit=logunit, unscale=unscale) call chksum_v_3d(arrayV_in, 'v '//mesg, HI_in, symmetric=symmetric, & - scale=scale, logunit=logunit) + scale=scale, logunit=logunit, unscale=unscale) endif end subroutine chksum_uv_3d !> Checksums a 2d array staggered at C-grid u points. subroutine chksum_u_2d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners, & - scale, logunit) + scale, logunit, unscale) type(hor_index_type), target, intent(in) :: HI_m !< A horizontal index type real, dimension(HI_m%IsdB:,HI_m%jsd:), target, intent(in) :: array_m !< The array to be checksummed in !! arbitrary, possibly rescaled units [A ~> a] @@ -915,6 +970,10 @@ subroutine chksum_u_2d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners, real, optional, intent(in) :: scale !< A factor to convert this array back to unscaled units !! for checksums and output [a A-1 ~> 1] integer, optional, intent(in) :: logunit !< IO unit for checksum logging + real, optional, intent(in) :: unscale !< A factor to convert this array back to unscaled units + !! for checksums and output [a A-1 ~> 1]. + !! Here scale and unscale are synonymous, but unscale + !! takes precedence if both are present. ! Local variables ! In the following comments, [A] is used to indicate the arbitrary, possibly rescaled units @@ -941,7 +1000,8 @@ subroutine chksum_u_2d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners, ! Arrays originating from v-points must be handled by vchksum allocate(array(HI%isd:HI%ied, HI%JsdB:HI%JedB)) call rotate_array(array_m, -turns, array) - call vchksum(array, mesg, HI, haloshift, symmetric, omit_corners, scale, logunit) + call vchksum(array, mesg, HI, haloshift, symmetric, omit_corners, & + scale=scale, logunit=logunit, unscale=unscale) return else allocate(array(HI%IsdB:HI%IedB, HI%jsd:HI%jed)) @@ -959,18 +1019,21 @@ subroutine chksum_u_2d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners, ! call chksum_error(FATAL, 'NaN detected in halo: '//trim(mesg)) endif - scaling = 1.0 ; if (present(scale)) scaling = scale + scaling = 1.0 + if (present(unscale)) then ; scaling = unscale + elseif (present(scale)) then ; scaling = scale ; endif + iounit = error_unit ; if (present(logunit)) iounit = logunit sym_stats = .false. ; if (present(symmetric)) sym_stats = symmetric if (present(haloshift)) then ; if (haloshift > 0) sym_stats = .true. ; endif if (calculateStatistics) then - if (present(scale)) then + if (present(unscale) .or. present(scale)) then allocate( rescaled_array(LBOUND(array,1):UBOUND(array,1), & LBOUND(array,2):UBOUND(array,2)), source=0.0 ) Is = HI%isc ; if (sym_stats) Is = HI%isc-1 do j=HI%jsc,HI%jec ; do I=Is,HI%IecB - rescaled_array(I,j) = scale*array(I,j) + rescaled_array(I,j) = scaling*array(I,j) enddo ; enddo call subStats(HI, rescaled_array, sym_stats, aMean, aMin, aMax) deallocate(rescaled_array) @@ -1039,19 +1102,19 @@ subroutine chksum_u_2d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners, contains - integer function subchk(array, HI, di, dj, scale) + integer function subchk(array, HI, di, dj, unscale) type(hor_index_type), intent(in) :: HI !< A horizontal index type real, dimension(HI%IsdB:,HI%jsd:), intent(in) :: array !< The array to be checksummed in !! arbitrary, possibly rescaled units [A ~> a] integer, intent(in) :: di !< i- direction array shift for this checksum integer, intent(in) :: dj !< j- direction array shift for this checksum - real, intent(in) :: scale !< A factor to convert this array back to unscaled units + real, intent(in) :: unscale !< A factor to convert this array back to unscaled units !! for checksums and output [a A-1 ~> 1] integer :: i, j, bc subchk = 0 ! This line deliberately uses the h-point computational domain. do j=HI%jsc+dj,HI%jec+dj ; do I=HI%isc+di,HI%iec+di - bc = bitcount(abs(scale*array(I,j))) + bc = bitcount(abs(unscale*array(I,j))) subchk = subchk + bc enddo ; enddo call sum_across_PEs(subchk) @@ -1089,7 +1152,7 @@ end subroutine chksum_u_2d !> Checksums a 2d array staggered at C-grid v points. subroutine chksum_v_2d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners, & - scale, logunit) + scale, logunit, unscale) type(hor_index_type), target, intent(in) :: HI_m !< A horizontal index type real, dimension(HI_m%isd:,HI_m%JsdB:), target, intent(in) :: array_m !< The array to be checksummed in !! arbitrary, possibly rescaled units [A ~> a] @@ -1101,6 +1164,10 @@ subroutine chksum_v_2d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners, real, optional, intent(in) :: scale !< A factor to convert this array back to unscaled units !! for checksums and output [a A-1 ~> 1] integer, optional, intent(in) :: logunit !< IO unit for checksum logging + real, optional, intent(in) :: unscale !< A factor to convert this array back to unscaled units + !! for checksums and output [a A-1 ~> 1]. + !! Here scale and unscale are synonymous, but unscale + !! takes precedence if both are present. ! Local variables ! In the following comments, [A] is used to indicate the arbitrary, possibly rescaled units @@ -1127,7 +1194,8 @@ subroutine chksum_v_2d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners, ! Arrays originating from u-points must be handled by uchksum allocate(array(HI%IsdB:HI%IedB, HI%jsd:HI%jed)) call rotate_array(array_m, -turns, array) - call uchksum(array, mesg, HI, haloshift, symmetric, omit_corners, scale, logunit) + call uchksum(array, mesg, HI, haloshift, symmetric, omit_corners, & + scale=scale, logunit=logunit, unscale=unscale) return else allocate(array(HI%isd:HI%ied, HI%JsdB:HI%JedB)) @@ -1145,18 +1213,21 @@ subroutine chksum_v_2d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners, ! call chksum_error(FATAL, 'NaN detected in halo: '//trim(mesg)) endif - scaling = 1.0 ; if (present(scale)) scaling = scale + scaling = 1.0 + if (present(unscale)) then ; scaling = unscale + elseif (present(scale)) then ; scaling = scale ; endif + iounit = error_unit ; if (present(logunit)) iounit = logunit sym_stats = .false. ; if (present(symmetric)) sym_stats = symmetric if (present(haloshift)) then ; if (haloshift > 0) sym_stats = .true. ; endif if (calculateStatistics) then - if (present(scale)) then + if (present(unscale) .or. present(scale)) then allocate( rescaled_array(LBOUND(array,1):UBOUND(array,1), & LBOUND(array,2):UBOUND(array,2)), source=0.0 ) Js = HI%jsc ; if (sym_stats) Js = HI%jsc-1 do J=Js,HI%JecB ; do i=HI%isc,HI%iec - rescaled_array(i,J) = scale*array(i,J) + rescaled_array(i,J) = scaling*array(i,J) enddo ; enddo call subStats(HI, rescaled_array, sym_stats, aMean, aMin, aMax) deallocate(rescaled_array) @@ -1225,19 +1296,19 @@ subroutine chksum_v_2d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners, contains - integer function subchk(array, HI, di, dj, scale) + integer function subchk(array, HI, di, dj, unscale) type(hor_index_type), intent(in) :: HI !< A horizontal index type real, dimension(HI%isd:,HI%JsdB:), intent(in) :: array !< The array to be checksummed in !! arbitrary, possibly rescaled units [A ~> a] integer, intent(in) :: di !< i- direction array shift for this checksum integer, intent(in) :: dj !< j- direction array shift for this checksum - real, intent(in) :: scale !< A factor to convert this array back to unscaled units + real, intent(in) :: unscale !< A factor to convert this array back to unscaled units !! for checksums and output [a A-1 ~> 1] integer :: i, j, bc subchk = 0 ! This line deliberately uses the h-point computational domain. do J=HI%jsc+dj,HI%jec+dj ; do i=HI%isc+di,HI%iec+di - bc = bitcount(abs(scale*array(i,J))) + bc = bitcount(abs(unscale*array(i,J))) subchk = subchk + bc enddo ; enddo call sum_across_PEs(subchk) @@ -1274,7 +1345,7 @@ end subroutine subStats end subroutine chksum_v_2d !> Checksums a 3d array staggered at tracer points. -subroutine chksum_h_3d(array_m, mesg, HI_m, haloshift, omit_corners, scale, logunit) +subroutine chksum_h_3d(array_m, mesg, HI_m, haloshift, omit_corners, scale, logunit, unscale) type(hor_index_type), target, intent(in) :: HI_m !< A horizontal index type real, dimension(HI_m%isd:,HI_m%jsd:,:), target, intent(in) :: array_m !< The array to be checksummed in !! arbitrary, possibly rescaled units [A ~> a] @@ -1284,6 +1355,10 @@ subroutine chksum_h_3d(array_m, mesg, HI_m, haloshift, omit_corners, scale, logu real, optional, intent(in) :: scale !< A factor to convert this array back to unscaled units !! for checksums and output [a A-1 ~> 1] integer, optional, intent(in) :: logunit !< IO unit for checksum logging + real, optional, intent(in) :: unscale !< A factor to convert this array back to unscaled units + !! for checksums and output [a A-1 ~> 1]. + !! Here scale and unscale are synonymous, but unscale + !! takes precedence if both are present. ! Local variables ! In the following comments, [A] is used to indicate the arbitrary, possibly rescaled units @@ -1320,16 +1395,19 @@ subroutine chksum_h_3d(array_m, mesg, HI_m, haloshift, omit_corners, scale, logu ! call chksum_error(FATAL, 'NaN detected in halo: '//trim(mesg)) endif - scaling = 1.0 ; if (present(scale)) scaling = scale + scaling = 1.0 + if (present(unscale)) then ; scaling = unscale + elseif (present(scale)) then ; scaling = scale ; endif + iounit = error_unit ; if (present(logunit)) iounit = logunit if (calculateStatistics) then - if (present(scale)) then + if (present(unscale) .or. present(scale)) then allocate( rescaled_array(LBOUND(array,1):UBOUND(array,1), & LBOUND(array,2):UBOUND(array,2), & LBOUND(array,3):UBOUND(array,3)), source=0.0 ) do k=1,size(array,3) ; do j=HI%jsc,HI%jec ; do i=HI%isc,HI%iec - rescaled_array(i,j,k) = scale*array(i,j,k) + rescaled_array(i,j,k) = scaling*array(i,j,k) enddo ; enddo ; enddo call subStats(HI, rescaled_array, aMean, aMin, aMax) @@ -1385,18 +1463,18 @@ subroutine chksum_h_3d(array_m, mesg, HI_m, haloshift, omit_corners, scale, logu contains - integer function subchk(array, HI, di, dj, scale) + integer function subchk(array, HI, di, dj, unscale) type(hor_index_type), intent(in) :: HI !< A horizontal index type real, dimension(HI%isd:,HI%jsd:,:), intent(in) :: array !< The array to be checksummed in !! arbitrary, possibly rescaled units [A ~> a] integer, intent(in) :: di !< i- direction array shift for this checksum integer, intent(in) :: dj !< j- direction array shift for this checksum - real, intent(in) :: scale !< A factor to convert this array back to unscaled units + real, intent(in) :: unscale !< A factor to convert this array back to unscaled units !! for checksums and output [a A-1 ~> 1] integer :: i, j, k, bc subchk = 0 do k=LBOUND(array,3),UBOUND(array,3) ; do j=HI%jsc+dj,HI%jec+dj ; do i=HI%isc+di,HI%iec+di - bc = bitcount(abs(scale*array(i,j,k))) + bc = bitcount(abs(unscale*array(i,j,k))) subchk = subchk + bc enddo ; enddo ; enddo call sum_across_PEs(subchk) @@ -1431,7 +1509,7 @@ end subroutine chksum_h_3d !> Checksums a 3d array staggered at corner points. subroutine chksum_B_3d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners, & - scale, logunit) + scale, logunit, unscale) type(hor_index_type), target, intent(in) :: HI_m !< A horizontal index type real, dimension(HI_m%IsdB:,HI_m%JsdB:,:), target, intent(in) :: array_m !< The array to be checksummed in !! arbitrary, possibly rescaled units [A ~> a] @@ -1443,6 +1521,10 @@ subroutine chksum_B_3d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners, real, optional, intent(in) :: scale !< A factor to convert this array back to unscaled units !! for checksums and output [a A-1 ~> 1] integer, optional, intent(in) :: logunit !< IO unit for checksum logging + real, optional, intent(in) :: unscale !< A factor to convert this array back to unscaled units + !! for checksums and output [a A-1 ~> 1]. + !! Here scale and unscale are synonymous, but unscale + !! takes precedence if both are present. ! Local variables ! In the following comments, [A] is used to indicate the arbitrary, possibly rescaled units @@ -1479,20 +1561,23 @@ subroutine chksum_B_3d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners, ! call chksum_error(FATAL, 'NaN detected in halo: '//trim(mesg)) endif - scaling = 1.0 ; if (present(scale)) scaling = scale + scaling = 1.0 + if (present(unscale)) then ; scaling = unscale + elseif (present(scale)) then ; scaling = scale ; endif + iounit = error_unit ; if (present(logunit)) iounit = logunit sym_stats = .false. ; if (present(symmetric)) sym_stats = symmetric if (present(haloshift)) then ; if (haloshift > 0) sym_stats = .true. ; endif if (calculateStatistics) then - if (present(scale)) then + if (present(unscale) .or. present(scale)) then allocate( rescaled_array(LBOUND(array,1):UBOUND(array,1), & LBOUND(array,2):UBOUND(array,2), & LBOUND(array,3):UBOUND(array,3)), source=0.0 ) Is = HI%isc ; if (sym_stats) Is = HI%isc-1 Js = HI%jsc ; if (sym_stats) Js = HI%jsc-1 do k=1,size(array,3) ; do J=Js,HI%JecB ; do I=Is,HI%IecB - rescaled_array(I,J,k) = scale*array(I,J,k) + rescaled_array(I,J,k) = scaling*array(I,J,k) enddo ; enddo ; enddo call subStats(HI, rescaled_array, sym_stats, aMean, aMin, aMax) deallocate(rescaled_array) @@ -1560,19 +1645,19 @@ subroutine chksum_B_3d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners, contains - integer function subchk(array, HI, di, dj, scale) + integer function subchk(array, HI, di, dj, unscale) type(hor_index_type), intent(in) :: HI !< A horizontal index type real, dimension(HI%IsdB:,HI%JsdB:,:), intent(in) :: array !< The array to be checksummed in !! arbitrary, possibly rescaled units [A ~> a] integer, intent(in) :: di !< i- direction array shift for this checksum integer, intent(in) :: dj !< j- direction array shift for this checksum - real, intent(in) :: scale !< A factor to convert this array back to unscaled units + real, intent(in) :: unscale !< A factor to convert this array back to unscaled units !! for checksums and output [a A-1 ~> 1] integer :: i, j, k, bc subchk = 0 ! This line deliberately uses the h-point computational domain. do k=LBOUND(array,3),UBOUND(array,3) ; do J=HI%jsc+dj,HI%jec+dj ; do I=HI%isc+di,HI%iec+di - bc = bitcount(abs(scale*array(I,J,k))) + bc = bitcount(abs(unscale*array(I,J,k))) subchk = subchk + bc enddo ; enddo ; enddo call sum_across_PEs(subchk) @@ -1610,7 +1695,7 @@ end subroutine chksum_B_3d !> Checksums a 3d array staggered at C-grid u points. subroutine chksum_u_3d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners, & - scale, logunit) + scale, logunit, unscale) type(hor_index_type), target, intent(in) :: HI_m !< A horizontal index type real, dimension(HI_m%isdB:,HI_m%Jsd:,:), target, intent(in) :: array_m !< The array to be checksummed in !! arbitrary, possibly rescaled units [A ~> a] @@ -1622,6 +1707,10 @@ subroutine chksum_u_3d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners, real, optional, intent(in) :: scale !< A factor to convert this array back to unscaled units !! for checksums and output [a A-1 ~> 1] integer, optional, intent(in) :: logunit !< IO unit for checksum logging + real, optional, intent(in) :: unscale !< A factor to convert this array back to unscaled units + !! for checksums and output [a A-1 ~> 1]. + !! Here scale and unscale are synonymous, but unscale + !! takes precedence if both are present. ! Local variables ! In the following comments, [A] is used to indicate the arbitrary, possibly rescaled units @@ -1648,7 +1737,8 @@ subroutine chksum_u_3d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners, ! Arrays originating from v-points must be handled by vchksum allocate(array(HI%isd:HI%ied, HI%JsdB:HI%JedB, size(array_m, 3))) call rotate_array(array_m, -turns, array) - call vchksum(array, mesg, HI, haloshift, symmetric, omit_corners, scale, logunit) + call vchksum(array, mesg, HI, haloshift, symmetric, omit_corners, & + scale=scale, logunit=logunit, unscale=unscale) return else allocate(array(HI%IsdB:HI%IedB, HI%jsd:HI%jed, size(array_m, 3))) @@ -1666,19 +1756,22 @@ subroutine chksum_u_3d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners, ! call chksum_error(FATAL, 'NaN detected in halo: '//trim(mesg)) endif - scaling = 1.0 ; if (present(scale)) scaling = scale + scaling = 1.0 + if (present(unscale)) then ; scaling = unscale + elseif (present(scale)) then ; scaling = scale ; endif + iounit = error_unit ; if (present(logunit)) iounit = logunit sym_stats = .false. ; if (present(symmetric)) sym_stats = symmetric if (present(haloshift)) then ; if (haloshift > 0) sym_stats = .true. ; endif if (calculateStatistics) then - if (present(scale)) then + if (present(unscale) .or. present(scale)) then allocate( rescaled_array(LBOUND(array,1):UBOUND(array,1), & LBOUND(array,2):UBOUND(array,2), & LBOUND(array,3):UBOUND(array,3)), source=0.0 ) Is = HI%isc ; if (sym_stats) Is = HI%isc-1 do k=1,size(array,3) ; do j=HI%jsc,HI%jec ; do I=Is,HI%IecB - rescaled_array(I,j,k) = scale*array(I,j,k) + rescaled_array(I,j,k) = scaling*array(I,j,k) enddo ; enddo ; enddo call subStats(HI, rescaled_array, sym_stats, aMean, aMin, aMax) deallocate(rescaled_array) @@ -1746,19 +1839,19 @@ subroutine chksum_u_3d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners, contains - integer function subchk(array, HI, di, dj, scale) + integer function subchk(array, HI, di, dj, unscale) type(hor_index_type), intent(in) :: HI !< A horizontal index type real, dimension(HI%IsdB:,HI%jsd:,:), intent(in) :: array !< The array to be checksummed in !! arbitrary, possibly rescaled units [A ~> a] integer, intent(in) :: di !< i- direction array shift for this checksum integer, intent(in) :: dj !< j- direction array shift for this checksum - real, intent(in) :: scale !< A factor to convert this array back to unscaled units + real, intent(in) :: unscale !< A factor to convert this array back to unscaled units !! for checksums and output [a A-1 ~> 1] integer :: i, j, k, bc subchk = 0 ! This line deliberately uses the h-point computational domain. do k=LBOUND(array,3),UBOUND(array,3) ; do j=HI%jsc+dj,HI%jec+dj ; do I=HI%isc+di,HI%iec+di - bc = bitcount(abs(scale*array(I,j,k))) + bc = bitcount(abs(unscale*array(I,j,k))) subchk = subchk + bc enddo ; enddo ; enddo call sum_across_PEs(subchk) @@ -1796,7 +1889,7 @@ end subroutine chksum_u_3d !> Checksums a 3d array staggered at C-grid v points. subroutine chksum_v_3d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners, & - scale, logunit) + scale, logunit, unscale) type(hor_index_type), target, intent(in) :: HI_m !< A horizontal index type real, dimension(HI_m%isd:,HI_m%JsdB:,:), target, intent(in) :: array_m !< The array to be checksummed in !! arbitrary, possibly rescaled units [A ~> a] @@ -1808,6 +1901,10 @@ subroutine chksum_v_3d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners, real, optional, intent(in) :: scale !< A factor to convert this array back to unscaled units !! for checksums and output [a A-1 ~> 1] integer, optional, intent(in) :: logunit !< IO unit for checksum logging + real, optional, intent(in) :: unscale !< A factor to convert this array back to unscaled units + !! for checksums and output [a A-1 ~> 1]. + !! Here scale and unscale are synonymous, but unscale + !! takes precedence if both are present. ! Local variables ! In the following comments, [A] is used to indicate the arbitrary, possibly rescaled units @@ -1834,7 +1931,8 @@ subroutine chksum_v_3d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners, ! Arrays originating from u-points must be handled by uchksum allocate(array(HI%IsdB:HI%IedB, HI%jsd:HI%jed, size(array_m, 3))) call rotate_array(array_m, -turns, array) - call uchksum(array, mesg, HI, haloshift, symmetric, omit_corners, scale, logunit) + call uchksum(array, mesg, HI, haloshift, symmetric, omit_corners, & + scale=scale, logunit=logunit, unscale=unscale) return else allocate(array(HI%isd:HI%ied, HI%JsdB:HI%JedB, size(array_m, 3))) @@ -1852,19 +1950,22 @@ subroutine chksum_v_3d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners, ! call chksum_error(FATAL, 'NaN detected in halo: '//trim(mesg)) endif - scaling = 1.0 ; if (present(scale)) scaling = scale + scaling = 1.0 + if (present(unscale)) then ; scaling = unscale + elseif (present(scale)) then ; scaling = scale ; endif + iounit = error_unit ; if (present(logunit)) iounit = logunit sym_stats = .false. ; if (present(symmetric)) sym_stats = symmetric if (present(haloshift)) then ; if (haloshift > 0) sym_stats = .true. ; endif if (calculateStatistics) then - if (present(scale)) then + if (present(unscale) .or. present(scale)) then allocate( rescaled_array(LBOUND(array,1):UBOUND(array,1), & LBOUND(array,2):UBOUND(array,2), & LBOUND(array,3):UBOUND(array,3)), source=0.0 ) Js = HI%jsc ; if (sym_stats) Js = HI%jsc-1 do k=1,size(array,3) ; do J=Js,HI%JecB ; do i=HI%isc,HI%iec - rescaled_array(i,J,k) = scale*array(i,J,k) + rescaled_array(i,J,k) = scaling*array(i,J,k) enddo ; enddo ; enddo call subStats(HI, rescaled_array, sym_stats, aMean, aMin, aMax) deallocate(rescaled_array) @@ -1932,19 +2033,19 @@ subroutine chksum_v_3d(array_m, mesg, HI_m, haloshift, symmetric, omit_corners, contains - integer function subchk(array, HI, di, dj, scale) + integer function subchk(array, HI, di, dj, unscale) type(hor_index_type), intent(in) :: HI !< A horizontal index type real, dimension(HI%isd:,HI%JsdB:,:), intent(in) :: array !< The array to be checksummed in !! arbitrary, possibly rescaled units [A ~> a] integer, intent(in) :: di !< i- direction array shift for this checksum integer, intent(in) :: dj !< j- direction array shift for this checksum - real, intent(in) :: scale !< A factor to convert this array back to unscaled units + real, intent(in) :: unscale !< A factor to convert this array back to unscaled units !! for checksums and output [a A-1 ~> 1] integer :: i, j, k, bc subchk = 0 ! This line deliberately uses the h-point computational domain. do k=LBOUND(array,3),UBOUND(array,3) ; do J=HI%jsc+dj,HI%jec+dj ; do i=HI%isc+di,HI%iec+di - bc = bitcount(abs(scale*array(i,J,k))) + bc = bitcount(abs(unscale*array(i,J,k))) subchk = subchk + bc enddo ; enddo ; enddo call sum_across_PEs(subchk)