From 8a1e9b502bf079d1f24a2ade3023d55fa4db9534 Mon Sep 17 00:00:00 2001 From: Matthew Harrison Date: Tue, 21 Apr 2020 10:39:48 -0400 Subject: [PATCH 01/11] New option TEMP_SALT_INIT_VERTICAL_REMAP_ONLY - This commit makes available an option to initialize the model temperature and salinity from a data set which is co-located with the model horizontal grid, but for which horizontal extrapolation and vertical remapping are needed. --- src/framework/MOM_horizontal_regridding.F90 | 79 ++++++++++--------- .../MOM_state_initialization.F90 | 9 ++- 2 files changed, 49 insertions(+), 39 deletions(-) diff --git a/src/framework/MOM_horizontal_regridding.F90 b/src/framework/MOM_horizontal_regridding.F90 index 7c19d715db..ebbaf3e3b4 100644 --- a/src/framework/MOM_horizontal_regridding.F90 +++ b/src/framework/MOM_horizontal_regridding.F90 @@ -272,7 +272,7 @@ end subroutine fill_miss_2d !> Extrapolate and interpolate from a file record subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, recnum, G, tr_z, & mask_z, z_in, z_edges_in, missing_value, reentrant_x, & - tripolar_n, homogenize, m_to_Z, answers_2018) + tripolar_n, homogenize, m_to_Z, answers_2018, ongrid) character(len=*), intent(in) :: filename !< Path to file containing tracer to be !! interpolated. @@ -296,6 +296,9 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, logical, optional, intent(in) :: answers_2018 !< If true, use expressions that give the same !! answers as the code did in late 2018. Otherwise !! add parentheses for rotational symmetry. + logical, optional, intent(in) :: ongrid !< If true, then data are assumed to have been interpolated + !! to the model horizontal grid. In this case, only + !! extrapolation is performed by this routine ! Local variables real, dimension(:,:), allocatable :: tr_in, tr_inp ! A 2-d array for holding input data on @@ -314,6 +317,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, real :: roundoff ! The magnitude of roundoff, usually ~2e-16. real :: add_offset, scale_factor logical :: add_np + logical :: is_ongrid character(len=8) :: laynum type(horiz_interp_type) :: Interp integer :: is, ie, js, je ! compute domain indices @@ -336,6 +340,8 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, id_clock_read = cpu_clock_id('(Initialize tracer from Z) read', grain=CLOCK_LOOP) + is_ongrid=.false. + if (present(ongrid)) is_ongrid=ongrid if (allocated(tr_z)) deallocate(tr_z) if (allocated(mask_z)) deallocate(mask_z) @@ -418,41 +424,38 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, if (present(m_to_Z)) then ; do k=1,kd ; z_in(k) = m_to_Z * z_in(k) ; enddo ; endif ! extrapolate the input data to the north pole using the northerm-most latitude - - max_lat = maxval(lat_in) add_np=.false. - if (max_lat < 90.0) then - add_np=.true. - jdp=jd+1 - allocate(lat_inp(jdp)) - lat_inp(1:jd)=lat_in(:) - lat_inp(jd+1)=90.0 - deallocate(lat_in) - allocate(lat_in(1:jdp)) - lat_in(:)=lat_inp(:) - else - jdp=jd + jdp=jd + if (.not. is_ongrid) then + max_lat = maxval(lat_in) + if (max_lat < 90.0) then + add_np=.true. + jdp=jd+1 + allocate(lat_inp(jdp)) + lat_inp(1:jd)=lat_in(:) + lat_inp(jd+1)=90.0 + deallocate(lat_in) + allocate(lat_in(1:jdp)) + lat_in(:)=lat_inp(:) + endif endif - ! construct level cell boundaries as the mid-point between adjacent centers z_edges_in(1) = 0.0 do K=2,kd - z_edges_in(K)=0.5*(z_in(k-1)+z_in(k)) + z_edges_in(K)=0.5*(z_in(k-1)+z_in(k)) enddo z_edges_in(kd+1)=2.0*z_in(kd) - z_in(kd-1) - call horiz_interp_init() - - lon_in = lon_in*PI_180 - lat_in = lat_in*PI_180 - allocate(x_in(id,jdp),y_in(id,jdp)) - call meshgrid(lon_in,lat_in, x_in, y_in) - - lon_out(:,:) = G%geoLonT(:,:)*PI_180 - lat_out(:,:) = G%geoLatT(:,:)*PI_180 - - + if (.not. is_ongrid) then + call horiz_interp_init() + lon_in = lon_in*PI_180 + lat_in = lat_in*PI_180 + allocate(x_in(id,jdp),y_in(id,jdp)) + call meshgrid(lon_in,lat_in, x_in, y_in) + lon_out(:,:) = G%geoLonT(:,:)*PI_180 + lat_out(:,:) = G%geoLatT(:,:)*PI_180 + endif allocate(tr_in(id,jd)) ; tr_in(:,:)=0.0 allocate(tr_inp(id,jdp)) ; tr_inp(:,:)=0.0 allocate(mask_in(id,jdp)) ; mask_in(:,:)=0.0 @@ -462,7 +465,6 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, call mpp_max(max_depth) if (z_edges_in(kd+1) Date: Tue, 21 Apr 2020 11:34:37 -0400 Subject: [PATCH 02/11] Improve performance of TS ongrid initialization - local distributed read from disk in the case of ongrid initialization rather than broadcasting the entire horizontal slice --- src/framework/MOM_horizontal_regridding.F90 | 114 ++++++++++++-------- 1 file changed, 71 insertions(+), 43 deletions(-) diff --git a/src/framework/MOM_horizontal_regridding.F90 b/src/framework/MOM_horizontal_regridding.F90 index e94349ceb2..66f58b5b9d 100644 --- a/src/framework/MOM_horizontal_regridding.F90 +++ b/src/framework/MOM_horizontal_regridding.F90 @@ -448,7 +448,10 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, enddo z_edges_in(kd+1)=2.0*z_in(kd) - z_in(kd-1) - if (.not. is_ongrid) then + if (is_ongrid) then + allocate(tr_in(is:ie,js:je)) ; tr_in(:,:)=0.0 + allocate(mask_in(is:ie,js:je)) ; mask_in(:,:)=0.0 + else call horiz_interp_init() lon_in = lon_in*PI_180 lat_in = lat_in*PI_180 @@ -456,11 +459,13 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, call meshgrid(lon_in,lat_in, x_in, y_in) lon_out(:,:) = G%geoLonT(:,:)*PI_180 lat_out(:,:) = G%geoLatT(:,:)*PI_180 + allocate(tr_in(id,jd)) ; tr_in(:,:)=0.0 + allocate(tr_inp(id,jdp)) ; tr_inp(:,:)=0.0 + allocate(mask_in(id,jdp)) ; mask_in(:,:)=0.0 + allocate(last_row(id)) ; last_row(:)=0.0 endif - allocate(tr_in(id,jd)) ; tr_in(:,:)=0.0 - allocate(tr_inp(id,jdp)) ; tr_inp(:,:)=0.0 - allocate(mask_in(id,jdp)) ; mask_in(:,:)=0.0 - allocate(last_row(id)) ; last_row(:)=0.0 + + max_depth = maxval(G%bathyT) call mpp_max(max_depth) @@ -473,50 +478,72 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, ! to define the layers do k=1,kd write(laynum,'(I8)') k ; laynum = adjustl(laynum) + mask_in=0.0 + if (is_ongrid) then + start(1) = is+G%HI%idg_offset ; start(2) = js+G%HI%jdg_offset ; start(3) = k + count(1) = ie-is+1 ; count(2) = je-js+1; count(3) = 1 + rcode = NF90_GET_VAR(ncid,varid, tr_in, start, count) + if (rcode /= 0) call MOM_error(FATAL,"hinterp_and_extract_from_Fie: "//& + "error reading level "//trim(laynum)//" of variable "//& + trim(varnam)//" in file "// trim(filename)) + + do j=js,je + do i=is,ie + if (abs(tr_in(i,j)-missing_value) > abs(roundoff*missing_value)) then + mask_in(i,j) = 1.0 + tr_in(i,j) = (tr_in(i,j)*scale_factor+add_offset) * conversion + else + tr_in(i,j) = missing_value + endif + enddo + enddo - if (is_root_pe()) then - start = 1 ; start(3) = k ; count(:) = 1 ; count(1) = id ; count(2) = jd - rcode = NF90_GET_VAR(ncid,varid, tr_in, start, count) - if (rcode /= 0) call MOM_error(FATAL,"hinterp_and_extract_from_Fie: "//& - "error reading level "//trim(laynum)//" of variable "//& - trim(varnam)//" in file "// trim(filename)) + else + if (is_root_pe()) then + start = 1 ; start(3) = k ; count(:) = 1 ; count(1) = id ; count(2) = jd + rcode = NF90_GET_VAR(ncid,varid, tr_in, start, count) + if (rcode /= 0) call MOM_error(FATAL,"hinterp_and_extract_from_Fie: "//& + "error reading level "//trim(laynum)//" of variable "//& + trim(varnam)//" in file "// trim(filename)) + + if (add_np) then + last_row(:)=tr_in(:,jd); pole=0.0;npole=0.0 + do i=1,id + if (abs(tr_in(i,jd)-missing_value) > abs(roundoff*missing_value)) then + pole = pole+last_row(i) + npole = npole+1.0 + endif + enddo + if (npole > 0) then + pole=pole/npole + else + pole=missing_value + endif + tr_inp(:,1:jd) = tr_in(:,:) + tr_inp(:,jdp) = pole + else + tr_inp(:,:) = tr_in(:,:) + endif + endif - if (add_np) then - last_row(:)=tr_in(:,jd); pole=0.0;npole=0.0 + call mpp_sync() + call mpp_broadcast(tr_inp, id*jdp, root_PE()) + call mpp_sync_self() + + do j=1,jdp do i=1,id - if (abs(tr_in(i,jd)-missing_value) > abs(roundoff*missing_value)) then - pole = pole+last_row(i) - npole = npole+1.0 - endif + if (abs(tr_inp(i,j)-missing_value) > abs(roundoff*missing_value)) then + mask_in(i,j) = 1.0 + tr_inp(i,j) = (tr_inp(i,j)*scale_factor+add_offset) * conversion + else + tr_inp(i,j) = missing_value + endif enddo - if (npole > 0) then - pole=pole/npole - else - pole=missing_value - endif - tr_inp(:,1:jd) = tr_in(:,:) - tr_inp(:,jdp) = pole - else - tr_inp(:,:) = tr_in(:,:) - endif - endif + enddo - call mpp_sync() - call mpp_broadcast(tr_inp, id*jdp, root_PE()) - call mpp_sync_self() + endif - mask_in=0.0 - do j=1,jdp - do i=1,id - if (abs(tr_inp(i,j)-missing_value) > abs(roundoff*missing_value)) then - mask_in(i,j) = 1.0 - tr_inp(i,j) = (tr_inp(i,j)*scale_factor+add_offset) * conversion - else - tr_inp(i,j) = missing_value - endif - enddo - enddo ! call fms routine horiz_interp to interpolate input level data to model horizontal grid if (.not. is_ongrid) then @@ -529,9 +556,10 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, call myStats(tr_inp,missing_value, is,ie,js,je,k,'Tracer from file') endif endif + tr_out(:,:) = 0.0 if (is_ongrid) then - tr_out(is:ie,js:je)=tr_inp(is+G%HI%idg_offset:ie+G%HI%idg_offset,js+G%HI%jdg_offset:je+G%HI%jdg_offset) + tr_out(is:ie,js:je)=tr_in(is:ie,js:je) else call horiz_interp(Interp,tr_inp,tr_out(is:ie,js:je), missing_value=missing_value, new_missing_handle=.true.) endif From dc6937259d1970e8ca1d91adcff77a4a16786150 Mon Sep 17 00:00:00 2001 From: Matthew Harrison Date: Tue, 21 Apr 2020 22:27:55 -0400 Subject: [PATCH 03/11] Fix parameter description for TEMP_SALT_INIT_VERTICAL_REMAP_ONLY --- src/initialization/MOM_state_initialization.F90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 4b290a2350..beeaf6e46a 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -2117,7 +2117,9 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param "This sets the default value for the various _2018_ANSWERS parameters.", & default=.true.) call get_param(PF, mdl, "TEMP_SALT_INIT_VERTICAL_REMAP_ONLY", pre_gridded, & - "This sets the default value for the various _2018_ANSWERS parameters.", & + "If true, initial conditions are on the model horizontal grid. " //& + "Extrapolation over missing ocean values is done using an ICE-9 "//& + "procedure with vertical ALE remapping .", & default=.false.) if (useALEremapping) then call get_param(PF, mdl, "REMAPPING_2018_ANSWERS", answers_2018, & From 3b0a2a682a53b24fa9ffda6f1aaa1f3d87f6d8d2 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 30 Apr 2020 07:10:41 -0400 Subject: [PATCH 04/11] Renamed all type(surface) variables to sfc_state Renamed all type(surface) variables to sfc_state, mostly from just state, for greater clarity and standardization within the MOM6 code. All answers are bitwise identical. --- src/core/MOM_checksum_packages.F90 | 35 +++++----- src/ice_shelf/MOM_ice_shelf.F90 | 88 +++++++++++++------------- src/tracer/DOME_tracer.F90 | 6 +- src/tracer/ISOMIP_tracer.F90 | 6 +- src/tracer/MOM_OCMIP2_CFC.F90 | 16 ++--- src/tracer/MOM_generic_tracer.F90 | 12 ++-- src/tracer/MOM_tracer_flow_control.F90 | 28 ++++---- src/tracer/advection_test_tracer.F90 | 6 +- src/tracer/boundary_impulse_tracer.F90 | 6 +- src/tracer/dye_example.F90 | 6 +- src/tracer/ideal_age_example.F90 | 6 +- src/tracer/oil_tracer.F90 | 6 +- src/tracer/pseudo_salt_tracer.F90 | 4 +- src/tracer/tracer_example.F90 | 6 +- src/user/BFB_surface_forcing.F90 | 13 ++-- src/user/Idealized_Hurricane.F90 | 40 ++++++------ src/user/SCM_CVMix_tests.F90 | 8 +-- src/user/dumbbell_surface_forcing.F90 | 16 ++--- src/user/user_revise_forcing.F90 | 4 +- 19 files changed, 155 insertions(+), 157 deletions(-) diff --git a/src/core/MOM_checksum_packages.F90 b/src/core/MOM_checksum_packages.F90 index 5fb96c1c73..dfe34676a6 100644 --- a/src/core/MOM_checksum_packages.F90 +++ b/src/core/MOM_checksum_packages.F90 @@ -141,15 +141,15 @@ end subroutine MOM_thermo_chksum ! ============================================================================= !> Write out chksums for the ocean surface variables. -subroutine MOM_surface_chksum(mesg, sfc, G, haloshift, symmetric) - character(len=*), intent(in) :: mesg !< A message that appears on the chksum lines. - type(surface), intent(inout) :: sfc !< transparent ocean surface state - !! structure shared with the calling routine - !! data in this structure is intent out. - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0). - logical, optional, intent(in) :: symmetric !< If true, do checksums on the fully symmetric - !! computational domain. +subroutine MOM_surface_chksum(mesg, sfc_state, G, haloshift, symmetric) + character(len=*), intent(in) :: mesg !< A message that appears on the chksum lines. + type(surface), intent(inout) :: sfc_state !< transparent ocean surface state structure + !! shared with the calling routine data in this + !! structure is intent out. + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0). + logical, optional, intent(in) :: symmetric !< If true, do checksums on the fully symmetric + !! computational domain. integer :: hs logical :: sym @@ -157,14 +157,15 @@ subroutine MOM_surface_chksum(mesg, sfc, G, haloshift, symmetric) sym = .false. ; if (present(symmetric)) sym = symmetric hs = 1 ; if (present(haloshift)) hs = haloshift - if (allocated(sfc%SST)) call hchksum(sfc%SST, mesg//" SST",G%HI,haloshift=hs) - if (allocated(sfc%SSS)) call hchksum(sfc%SSS, mesg//" SSS",G%HI,haloshift=hs) - if (allocated(sfc%sea_lev)) call hchksum(sfc%sea_lev, mesg//" sea_lev",G%HI,haloshift=hs) - if (allocated(sfc%Hml)) call hchksum(sfc%Hml, mesg//" Hml",G%HI,haloshift=hs) - if (allocated(sfc%u) .and. allocated(sfc%v)) & - call uvchksum(mesg//" SSU", sfc%u, sfc%v, G%HI, haloshift=hs, symmetric=sym) -! if (allocated(sfc%salt_deficit)) call hchksum(sfc%salt_deficit, mesg//" salt deficit",G%HI,haloshift=hs) - if (allocated(sfc%frazil)) call hchksum(sfc%frazil, mesg//" frazil", G%HI, haloshift=hs) + if (allocated(sfc_state%SST)) call hchksum(sfc_state%SST, mesg//" SST", G%HI, haloshift=hs) + if (allocated(sfc_state%SSS)) call hchksum(sfc_state%SSS, mesg//" SSS", G%HI, haloshift=hs) + if (allocated(sfc_state%sea_lev)) call hchksum(sfc_state%sea_lev, mesg//" sea_lev", G%HI, haloshift=hs) + if (allocated(sfc_state%Hml)) call hchksum(sfc_state%Hml, mesg//" Hml", G%HI, haloshift=hs) + 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) +! if (allocated(sfc_state%salt_deficit)) & +! call hchksum(sfc_state%salt_deficit, mesg//" salt deficit", G%HI, haloshift=hs) + if (allocated(sfc_state%frazil)) call hchksum(sfc_state%frazil, mesg//" frazil", G%HI, haloshift=hs) end subroutine MOM_surface_chksum diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 0db52e57e5..819ce8bb76 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -192,8 +192,8 @@ module MOM_ice_shelf !> Calculates fluxes between the ocean and ice-shelf using the three-equations !! formulation (optional to use just two equations). !! See \ref section_ICE_SHELF_equations -subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) - type(surface), intent(inout) :: state !< A structure containing fields that +subroutine shelf_calc_flux(sfc_state, fluxes, Time, time_step, CS, forces) + type(surface), intent(inout) :: sfc_state !< A structure containing fields that !! describe the surface state of the ocean. The !! intent is only inout to allow for halo updates. type(forcing), intent(inout) :: fluxes !< structure containing pointers to any possible @@ -319,7 +319,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) ISS%salt_flux(:,:) = 0.0 ; ISS%tflux_ocn(:,:) = 0.0 ; ISS%tfreeze(:,:) = 0.0 ! define Sbdry to avoid Run-Time Check Failure, when melt is not computed. haline_driving(:,:) = 0.0 - Sbdry(:,:) = state%sss(:,:) + Sbdry(:,:) = sfc_state%sss(:,:) !update time CS%Time = Time @@ -332,16 +332,16 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) if (CS%debug) then call hchksum(fluxes%frac_shelf_h, "frac_shelf_h before apply melting", G%HI, haloshift=0) - call hchksum(state%sst, "sst before apply melting", G%HI, haloshift=0) - call hchksum(state%sss, "sss before apply melting", G%HI, haloshift=0) - call hchksum(state%u, "u_ml before apply melting", G%HI, haloshift=0) - call hchksum(state%v, "v_ml before apply melting", G%HI, haloshift=0) - call hchksum(state%ocean_mass, "ocean_mass before apply melting", G%HI, haloshift=0) + call hchksum(sfc_state%sst, "sst before apply melting", G%HI, haloshift=0) + call hchksum(sfc_state%sss, "sss before apply melting", G%HI, haloshift=0) + call hchksum(sfc_state%u, "u_ml before apply melting", G%HI, haloshift=0) + call hchksum(sfc_state%v, "v_ml before apply melting", G%HI, haloshift=0) + call hchksum(sfc_state%ocean_mass, "ocean_mass before apply melting", G%HI, haloshift=0) endif ! Calculate the friction velocity under ice shelves, using taux_shelf and tauy_shelf if possible. - if (allocated(state%taux_shelf) .and. allocated(state%tauy_shelf)) then - call pass_vector(state%taux_shelf, state%tauy_shelf, G%domain, TO_ALL, CGRID_NE) + if (allocated(sfc_state%taux_shelf) .and. allocated(sfc_state%tauy_shelf)) then + call pass_vector(sfc_state%taux_shelf, sfc_state%tauy_shelf, G%domain, TO_ALL, CGRID_NE) endif Irho0 = US%m_s_to_L_T**2*US%kg_m3_to_R / CS%Rho_ocn do j=js,je ; do i=is,ie ; if (fluxes%frac_shelf_h(i,j) > 0.0) then @@ -352,12 +352,12 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) asv2 = (ISS%area_shelf_h(i,j) + ISS%area_shelf_h(i,j+1)) I_au = 0.0 ; if (asu1 + asu2 > 0.0) I_au = 1.0 / (asu1 + asu2) I_av = 0.0 ; if (asv1 + asv2 > 0.0) I_av = 1.0 / (asv1 + asv2) - if (allocated(state%taux_shelf) .and. allocated(state%tauy_shelf)) then - taux2 = (asu1 * state%taux_shelf(I-1,j)**2 + asu2 * state%taux_shelf(I,j)**2 ) * I_au - tauy2 = (asv1 * state%tauy_shelf(i,J-1)**2 + asv2 * state%tauy_shelf(i,J)**2 ) * I_av + if (allocated(sfc_state%taux_shelf) .and. allocated(sfc_state%tauy_shelf)) then + taux2 = (asu1 * sfc_state%taux_shelf(I-1,j)**2 + asu2 * sfc_state%taux_shelf(I,j)**2 ) * I_au + tauy2 = (asv1 * sfc_state%tauy_shelf(i,J-1)**2 + asv2 * sfc_state%tauy_shelf(i,J)**2 ) * I_av endif - u2_av = US%m_s_to_L_T**2*(asu1 * state%u(I-1,j)**2 + asu2 * state%u(I,j)**2) * I_au - v2_av = US%m_s_to_L_T**2*(asv1 * state%v(i,J-1)**2 + asu2 * state%v(i,J)**2) * I_av + u2_av = US%m_s_to_L_T**2*(asu1 * sfc_state%u(I-1,j)**2 + asu2 * sfc_state%u(I,j)**2) * I_au + v2_av = US%m_s_to_L_T**2*(asv1 * sfc_state%v(i,J-1)**2 + asu2 * sfc_state%v(i,J)**2) * I_av if (taux2 + tauy2 > 0.0) then fluxes%ustar_shelf(i,j) = MAX(CS%ustar_bg, US%L_to_Z * & @@ -377,13 +377,13 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) do i=is,ie ; p_int(i) = CS%g_Earth * ISS%mass_shelf(i,j) ; enddo ! Calculate insitu densities and expansion coefficients - call calculate_density(state%sst(:,j), state%sss(:,j), p_int, Rhoml(:), & + call calculate_density(sfc_state%sst(:,j), sfc_state%sss(:,j), p_int, Rhoml(:), & CS%eqn_of_state, EOSdom) - call calculate_density_derivs(state%sst(:,j), state%sss(:,j), p_int, dR0_dT, dR0_dS, & + call calculate_density_derivs(sfc_state%sst(:,j), sfc_state%sss(:,j), p_int, dR0_dT, dR0_dS, & CS%eqn_of_state, EOSdom) do i=is,ie - if ((state%ocean_mass(i,j) > US%RZ_to_kg_m2*CS%col_mass_melt_threshold) .and. & + if ((sfc_state%ocean_mass(i,j) > US%RZ_to_kg_m2*CS%col_mass_melt_threshold) .and. & (ISS%area_shelf_h(i,j) > 0.0) .and. CS%isthermo) then if (CS%threeeq) then @@ -397,7 +397,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) ! reported ocean mixed layer thickness and the neutral Ekman depth. absf = 0.25*((abs(G%CoriolisBu(I,J)) + abs(G%CoriolisBu(I-1,J-1))) + & (abs(G%CoriolisBu(I,J-1)) + abs(G%CoriolisBu(I-1,J)))) - if (absf*US%m_to_Z*state%Hml(i,j) <= VK*ustar_h) then ; hBL_neut = US%m_to_Z*state%Hml(i,j) + if (absf*US%m_to_Z*sfc_state%Hml(i,j) <= VK*ustar_h) then ; hBL_neut = US%m_to_Z*sfc_state%Hml(i,j) else ; hBL_neut = (VK*ustar_h) / absf ; endif hBL_neut_h_molec = ZETA_N * ((hBL_neut * ustar_h) / (5.0 * CS%kv_molec)) @@ -414,9 +414,9 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) ! S_a is always < 0.0 with a realistic expression for the freezing point. S_a = CS%dTFr_dS * CS%Gamma_T_3EQ * CS%Cp - S_b = CS%Gamma_T_3EQ*CS%Cp*(CS%TFr_0_0 + CS%dTFr_dp*p_int(i) - state%sst(i,j)) - & + S_b = CS%Gamma_T_3EQ*CS%Cp*(CS%TFr_0_0 + CS%dTFr_dp*p_int(i) - sfc_state%sst(i,j)) - & CS%Lat_fusion * CS%Gamma_S_3EQ ! S_b Can take either sign, but is usually negative. - S_c = CS%Lat_fusion * CS%Gamma_S_3EQ * state%sss(i,j) ! Always >= 0 + S_c = CS%Lat_fusion * CS%Gamma_S_3EQ * sfc_state%sss(i,j) ! Always >= 0 if (S_c == 0.0) then ! The solution for fresh water. Sbdry(i,j) = 0.0 @@ -434,7 +434,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) ! Safety check if (Sbdry(i,j) < 0.) then - write(mesg,*) 'state%sss(i,j) = ',state%sss(i,j), 'S_a, S_b, S_c', S_a, S_b, S_c + write(mesg,*) 'sfc_state%sss(i,j) = ',sfc_state%sss(i,j), 'S_a, S_b, S_c', S_a, S_b, S_c call MOM_error(WARNING, mesg, .true.) write(mesg,*) 'I,J,Sbdry1,Sbdry2',i,j,Sbdry1,Sbdry2 call MOM_error(WARNING, mesg, .true.) @@ -442,7 +442,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) endif else ! Guess sss as the iteration starting point for the boundary salinity. - Sbdry(i,j) = state%sss(i,j) ; Sb_max_set = .false. + Sbdry(i,j) = sfc_state%sss(i,j) ; Sb_max_set = .false. Sb_min_set = .false. endif !find_salt_root @@ -451,8 +451,8 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) call calculate_TFreeze(Sbdry(i,j), p_int(i), ISS%tfreeze(i,j), CS%eqn_of_state, & pres_scale=US%RL2_T2_to_Pa) - dT_ustar = (ISS%tfreeze(i,j) - state%sst(i,j)) * ustar_h - dS_ustar = (Sbdry(i,j) - state%sss(i,j)) * ustar_h + dT_ustar = (ISS%tfreeze(i,j) - sfc_state%sst(i,j)) * ustar_h + dS_ustar = (Sbdry(i,j) - sfc_state%sss(i,j)) * ustar_h ! First, determine the buoyancy flux assuming no effects of stability ! on the turbulence. Following H & J '99, this limit also applies @@ -558,10 +558,10 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) else mass_exch = exch_vel_s(i,j) * CS%Rho_ocn - Sbdry_it = (state%sss(i,j) * mass_exch + CS%Salin_ice * ISS%water_flux(i,j)) / & + Sbdry_it = (sfc_state%sss(i,j) * mass_exch + CS%Salin_ice * ISS%water_flux(i,j)) / & (mass_exch + ISS%water_flux(i,j)) dS_it = Sbdry_it - Sbdry(i,j) - if (abs(dS_it) < 1e-4*(0.5*(state%sss(i,j) + Sbdry(i,j) + 1.e-10))) exit + if (abs(dS_it) < 1.0e-4*(0.5*(sfc_state%sss(i,j) + Sbdry(i,j) + 1.0e-10))) exit if (dS_it < 0.0) then ! Sbdry is now the upper bound. @@ -592,11 +592,11 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) ! is specified and large enough that the ocean salinity at the interface ! is about the same as the boundary layer salinity. - call calculate_TFreeze(state%sss(i,j), p_int(i), ISS%tfreeze(i,j), CS%eqn_of_state, & + call calculate_TFreeze(sfc_state%sss(i,j), p_int(i), ISS%tfreeze(i,j), CS%eqn_of_state, & pres_scale=US%RL2_T2_to_Pa) exch_vel_t(i,j) = CS%gamma_t - ISS%tflux_ocn(i,j) = RhoCp * exch_vel_t(i,j) * (ISS%tfreeze(i,j) - state%sst(i,j)) + ISS%tflux_ocn(i,j) = RhoCp * exch_vel_t(i,j) * (ISS%tfreeze(i,j) - sfc_state%sst(i,j)) ISS%tflux_shelf(i,j) = 0.0 ISS%water_flux(i,j) = -I_LF * ISS%tflux_ocn(i,j) Sbdry(i,j) = 0.0 @@ -607,7 +607,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) ISS%tflux_ocn(i,j) = 0.0 endif -! haline_driving(:,:) = state%sss(i,j) - Sbdry(i,j) +! haline_driving(:,:) = sfc_state%sss(i,j) - Sbdry(i,j) enddo ! i-loop enddo ! j-loop @@ -616,7 +616,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) fluxes%iceshelf_melt(:,:) = ISS%water_flux(:,:) * CS%flux_factor do j=js,je ; do i=is,ie - if ((state%ocean_mass(i,j) > US%RZ_to_kg_m2*CS%col_mass_melt_threshold) .and. & + if ((sfc_state%ocean_mass(i,j) > US%RZ_to_kg_m2*CS%col_mass_melt_threshold) .and. & (ISS%area_shelf_h(i,j) > 0.0) .and. (CS%isthermo)) then ! Set melt to zero above a cutoff pressure (CS%Rho_ocn*CS%cutoff_depth*CS%g_Earth). @@ -630,11 +630,11 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) !!!!!!!!!!!!!!!!!!!!!!!!!!!!Safety checks !!!!!!!!!!!!!!!!!!!!!!!!! !1)Check if haline_driving computed above is consistent with - ! haline_driving = state%sss - Sbdry + ! haline_driving = sfc_state%sss - Sbdry !if (fluxes%iceshelf_melt(i,j) /= 0.0) then - ! if (haline_driving(i,j) /= (state%sss(i,j) - Sbdry(i,j))) then + ! if (haline_driving(i,j) /= (sfc_state%sss(i,j) - Sbdry(i,j))) then ! write(mesg,*) 'at i,j=',i,j,' haline_driving, sss-Sbdry',haline_driving(i,j), & - ! (state%sss(i,j) - Sbdry(i,j)) + ! (sfc_state%sss(i,j) - Sbdry(i,j)) ! call MOM_error(FATAL, & ! "shelf_calc_flux: Inconsistency in melt and haline_driving"//trim(mesg)) ! endif @@ -679,7 +679,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) if (CS%debug) call MOM_forcing_chksum("Before add shelf flux", fluxes, G, CS%US, haloshift=0) - call add_shelf_flux(G, US, CS, state, fluxes) + call add_shelf_flux(G, US, CS, sfc_state, fluxes) ! now the thermodynamic data is passed on... time to update the ice dynamic quantities @@ -690,7 +690,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) ! advect the ice shelf, and advance the front. Calving will be in here somewhere as well.. ! when we decide on how to do it call update_ice_shelf(CS%dCS, ISS, G, US, US%s_to_T*time_step, Time, & - US%kg_m3_to_R*US%m_to_Z*state%ocean_mass(:,:), coupled_GL) + US%kg_m3_to_R*US%m_to_Z*sfc_state%ocean_mass(:,:), coupled_GL) endif @@ -699,12 +699,12 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) if (CS%id_area_shelf_h > 0) call post_data(CS%id_area_shelf_h, ISS%area_shelf_h, CS%diag) if (CS%id_ustar_shelf > 0) call post_data(CS%id_ustar_shelf, fluxes%ustar_shelf, CS%diag) if (CS%id_melt > 0) call post_data(CS%id_melt, fluxes%iceshelf_melt, CS%diag) - if (CS%id_thermal_driving > 0) call post_data(CS%id_thermal_driving, (state%sst-ISS%tfreeze), CS%diag) + if (CS%id_thermal_driving > 0) call post_data(CS%id_thermal_driving, (sfc_state%sst-ISS%tfreeze), CS%diag) if (CS%id_Sbdry > 0) call post_data(CS%id_Sbdry, Sbdry, CS%diag) if (CS%id_haline_driving > 0) call post_data(CS%id_haline_driving, haline_driving, CS%diag) if (CS%id_mass_flux > 0) call post_data(CS%id_mass_flux, mass_flux, CS%diag) - if (CS%id_u_ml > 0) call post_data(CS%id_u_ml, state%u, CS%diag) - if (CS%id_v_ml > 0) call post_data(CS%id_v_ml, state%v, CS%diag) + if (CS%id_u_ml > 0) call post_data(CS%id_u_ml, sfc_state%u, CS%diag) + if (CS%id_v_ml > 0) call post_data(CS%id_v_ml, sfc_state%v, CS%diag) if (CS%id_tfreeze > 0) call post_data(CS%id_tfreeze, ISS%tfreeze, CS%diag) if (CS%id_tfl_shelf > 0) call post_data(CS%id_tfl_shelf, ISS%tflux_shelf, CS%diag) if (CS%id_exch_vel_t > 0) call post_data(CS%id_exch_vel_t, exch_vel_t, CS%diag) @@ -882,11 +882,11 @@ subroutine add_shelf_pressure(G, US, CS, fluxes) end subroutine add_shelf_pressure !> Updates surface fluxes that are influenced by sub-ice-shelf melting -subroutine add_shelf_flux(G, US, CS, state, fluxes) +subroutine add_shelf_flux(G, US, CS, sfc_state, fluxes) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(ice_shelf_CS), pointer :: CS !< This module's control structure. - type(surface), intent(inout) :: state!< Surface ocean state + type(surface), intent(inout) :: sfc_state !< Surface ocean state type(forcing), intent(inout) :: fluxes !< A structure of surface fluxes that may be used/updated. ! local variables @@ -931,8 +931,8 @@ subroutine add_shelf_flux(G, US, CS, state, fluxes) ! vertical decay scale. if (CS%debug) then - if (allocated(state%taux_shelf) .and. allocated(state%tauy_shelf)) then - call uvchksum("tau[xy]_shelf", state%taux_shelf, state%tauy_shelf, & + if (allocated(sfc_state%taux_shelf) .and. allocated(sfc_state%tauy_shelf)) then + call uvchksum("tau[xy]_shelf", sfc_state%taux_shelf, sfc_state%tauy_shelf, & G%HI, haloshift=0) endif endif @@ -1023,7 +1023,7 @@ subroutine add_shelf_flux(G, US, CS, state, fluxes) ! get total ice shelf mass at (Time-dt) and (Time), in kg do j=js,je ; do i=is,ie ! Just consider the change in the mass of the floating shelf. - if ((state%ocean_mass(i,j) > US%RZ_to_kg_m2*CS%min_ocean_mass_float) .and. & + if ((sfc_state%ocean_mass(i,j) > US%RZ_to_kg_m2*CS%min_ocean_mass_float) .and. & (ISS%area_shelf_h(i,j) > 0.0)) then delta_float_mass(i,j) = ISS%mass_shelf(i,j) - last_mass_shelf(i,j) else diff --git a/src/tracer/DOME_tracer.F90 b/src/tracer/DOME_tracer.F90 index f8bc58c8d8..7396a4092a 100644 --- a/src/tracer/DOME_tracer.F90 +++ b/src/tracer/DOME_tracer.F90 @@ -338,9 +338,9 @@ end subroutine DOME_tracer_column_physics !> This subroutine extracts the surface fields from this tracer package that !! are to be shared with the atmosphere in coupled configurations. !! This particular tracer package does not report anything back to the coupler. -subroutine DOME_tracer_surface_state(state, h, G, CS) +subroutine DOME_tracer_surface_state(sfc_state, h, G, CS) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - type(surface), intent(inout) :: state !< A structure containing fields that + type(surface), intent(inout) :: sfc_state !< A structure containing fields that !! describe the surface state of the ocean. real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & intent(in) :: h !< Layer thickness [H ~> m or kg m-2]. @@ -361,7 +361,7 @@ subroutine DOME_tracer_surface_state(state, h, G, CS) ! This call loads the surface values into the appropriate array in the ! coupler-type structure. call coupler_type_set_data(CS%tr(:,:,1,m), CS%ind_tr(m), ind_csurf, & - state%tr_fields, idim=(/isd, is, ie, ied/), & + sfc_state%tr_fields, idim=(/isd, is, ie, ied/), & jdim=(/jsd, js, je, jed/) ) enddo endif diff --git a/src/tracer/ISOMIP_tracer.F90 b/src/tracer/ISOMIP_tracer.F90 index 95d451791e..c9bf98f7ff 100644 --- a/src/tracer/ISOMIP_tracer.F90 +++ b/src/tracer/ISOMIP_tracer.F90 @@ -325,9 +325,9 @@ end subroutine ISOMIP_tracer_column_physics !> This subroutine extracts the surface fields from this tracer package that !! are to be shared with the atmosphere in coupled configurations. !! This particular tracer package does not report anything back to the coupler. -subroutine ISOMIP_tracer_surface_state(state, h, G, CS) +subroutine ISOMIP_tracer_surface_state(sfc_state, h, G, CS) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - type(surface), intent(inout) :: state !< A structure containing fields that + type(surface), intent(inout) :: sfc_state !< A structure containing fields that !! describe the surface state of the ocean. real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & intent(in) :: h !< Layer thickness [H ~> m or kg m-2]. @@ -348,7 +348,7 @@ subroutine ISOMIP_tracer_surface_state(state, h, G, CS) ! This call loads the surface values into the appropriate array in the ! coupler-type structure. call coupler_type_set_data(CS%tr(:,:,1,m), CS%ind_tr(m), ind_csurf, & - state%tr_fields, idim=(/isd, is, ie, ied/), & + sfc_state%tr_fields, idim=(/isd, is, ie, ied/), & jdim=(/jsd, js, je, jed/) ) enddo endif diff --git a/src/tracer/MOM_OCMIP2_CFC.F90 b/src/tracer/MOM_OCMIP2_CFC.F90 index a95ea654f4..9aad84a6dd 100644 --- a/src/tracer/MOM_OCMIP2_CFC.F90 +++ b/src/tracer/MOM_OCMIP2_CFC.F90 @@ -542,9 +542,9 @@ end function OCMIP2_CFC_stock !> This subroutine extracts the surface CFC concentrations and other fields that !! are shared with the atmosphere to calculate CFC fluxes. -subroutine OCMIP2_CFC_surface_state(state, h, G, CS) +subroutine OCMIP2_CFC_surface_state(sfc_state, h, G, CS) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - type(surface), intent(inout) :: state !< A structure containing fields that + type(surface), intent(inout) :: sfc_state !< A structure containing fields that !! describe the surface state of the ocean. real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & intent(in) :: h !< Layer thickness [H ~> m or kg m-2]. @@ -572,8 +572,8 @@ subroutine OCMIP2_CFC_surface_state(state, h, G, CS) if (.not.associated(CS)) return do j=js,je ; do i=is,ie - ta = max(0.01, (state%SST(i,j) + 273.15) * 0.01) ! Why is this in hectoKelvin? - sal = state%SSS(i,j) ; SST = state%SST(i,j) + ta = max(0.01, (sfc_state%SST(i,j) + 273.15) * 0.01) ! Why is this in hectoKelvin? + sal = sfc_state%SSS(i,j) ; SST = sfc_state%SST(i,j) ! Calculate solubilities using Warner and Weiss (1985) DSR, vol 32. ! The final result is in mol/cm3/pptv (1 part per trillion 1e-12) ! Use Bullister and Wisegavger for CCl4. @@ -603,13 +603,13 @@ subroutine OCMIP2_CFC_surface_state(state, h, G, CS) ! These calls load these values into the appropriate arrays in the ! coupler-type structure. call coupler_type_set_data(CFC11_alpha, CS%ind_cfc_11_flux, ind_alpha, & - state%tr_fields, idim=idim, jdim=jdim) + sfc_state%tr_fields, idim=idim, jdim=jdim) call coupler_type_set_data(CFC11_Csurf, CS%ind_cfc_11_flux, ind_csurf, & - state%tr_fields, idim=idim, jdim=jdim) + sfc_state%tr_fields, idim=idim, jdim=jdim) call coupler_type_set_data(CFC12_alpha, CS%ind_cfc_12_flux, ind_alpha, & - state%tr_fields, idim=idim, jdim=jdim) + sfc_state%tr_fields, idim=idim, jdim=jdim) call coupler_type_set_data(CFC12_Csurf, CS%ind_cfc_12_flux, ind_csurf, & - state%tr_fields, idim=idim, jdim=jdim) + sfc_state%tr_fields, idim=idim, jdim=jdim) end subroutine OCMIP2_CFC_surface_state diff --git a/src/tracer/MOM_generic_tracer.F90 b/src/tracer/MOM_generic_tracer.F90 index 83c2c9a8e7..e68833c3cd 100644 --- a/src/tracer/MOM_generic_tracer.F90 +++ b/src/tracer/MOM_generic_tracer.F90 @@ -706,9 +706,9 @@ end function MOM_generic_tracer_min_max !! !! This subroutine sets up the fields that the coupler needs to calculate the !! CFC fluxes between the ocean and atmosphere. - subroutine MOM_generic_tracer_surface_state(state, h, G, CS) + subroutine MOM_generic_tracer_surface_state(sfc_state, h, G, CS) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure - type(surface), intent(inout) :: state !< A structure containing fields that + type(surface), intent(inout) :: sfc_state !< A structure containing fields that !! describe the surface state of the ocean. real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] type(MOM_generic_tracer_CS), pointer :: CS !< Pointer to the control structure for this module. @@ -727,11 +727,11 @@ subroutine MOM_generic_tracer_surface_state(state, h, G, CS) dzt(:,:,:) = CS%H_to_m * h(:,:,:) - sosga = global_area_mean(state%SSS, G) + sosga = global_area_mean(sfc_state%SSS, G) - call generic_tracer_coupler_set(state%tr_fields,& - ST=state%SST,& - SS=state%SSS,& + call generic_tracer_coupler_set(sfc_state%tr_fields,& + ST=sfc_state%SST,& + SS=sfc_state%SSS,& rho=rho0,& !nnz: required for MOM5 and previous versions. ilb=G%isd, jlb=G%jsd,& dzt=dzt,& !This is needed for the Mocsy method of carbonate system vars diff --git a/src/tracer/MOM_tracer_flow_control.F90 b/src/tracer/MOM_tracer_flow_control.F90 index 6e28477d26..86003605f7 100644 --- a/src/tracer/MOM_tracer_flow_control.F90 +++ b/src/tracer/MOM_tracer_flow_control.F90 @@ -378,9 +378,9 @@ end subroutine get_chl_from_model !> This subroutine calls the individual tracer modules' subroutines to !! specify or read quantities related to their surface forcing. -subroutine call_tracer_set_forcing(state, fluxes, day_start, day_interval, G, CS) +subroutine call_tracer_set_forcing(sfc_state, fluxes, day_start, day_interval, G, CS) - type(surface), intent(inout) :: state !< A structure containing fields that + type(surface), intent(inout) :: sfc_state !< A structure containing fields that !! describe the surface state of the !! ocean. type(forcing), intent(inout) :: fluxes !< A structure containing pointers to any @@ -396,7 +396,7 @@ subroutine call_tracer_set_forcing(state, fluxes, day_start, day_interval, G, CS if (.not. associated(CS)) call MOM_error(FATAL, "call_tracer_set_forcing"// & "Module must be initialized via call_tracer_register before it is used.") ! if (CS%use_ideal_age) & -! call ideal_age_tracer_set_forcing(state, fluxes, day_start, day_interval, & +! call ideal_age_tracer_set_forcing(sfc_state, fluxes, day_start, day_interval, & ! G, CS%ideal_age_tracer_CSp) end subroutine call_tracer_set_forcing @@ -755,8 +755,8 @@ end subroutine store_stocks !> This subroutine calls all registered tracer packages to enable them to !! add to the surface state returned to the coupler. These routines are optional. -subroutine call_tracer_surface_state(state, h, G, CS) - type(surface), intent(inout) :: state !< A structure containing fields that +subroutine call_tracer_surface_state(sfc_state, h, G, CS) + type(surface), intent(inout) :: sfc_state !< A structure containing fields that !! describe the surface state of the ocean. real, dimension(NIMEM_,NJMEM_,NKMEM_), & intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] @@ -769,24 +769,24 @@ subroutine call_tracer_surface_state(state, h, G, CS) ! Add other user-provided calls here. if (CS%use_USER_tracer_example) & - call USER_tracer_surface_state(state, h, G, CS%USER_tracer_example_CSp) + call USER_tracer_surface_state(sfc_state, h, G, CS%USER_tracer_example_CSp) if (CS%use_DOME_tracer) & - call DOME_tracer_surface_state(state, h, G, CS%DOME_tracer_CSp) + call DOME_tracer_surface_state(sfc_state, h, G, CS%DOME_tracer_CSp) if (CS%use_ISOMIP_tracer) & - call ISOMIP_tracer_surface_state(state, h, G, CS%ISOMIP_tracer_CSp) + call ISOMIP_tracer_surface_state(sfc_state, h, G, CS%ISOMIP_tracer_CSp) if (CS%use_ideal_age) & - call ideal_age_tracer_surface_state(state, h, G, CS%ideal_age_tracer_CSp) + call ideal_age_tracer_surface_state(sfc_state, h, G, CS%ideal_age_tracer_CSp) if (CS%use_regional_dyes) & - call dye_tracer_surface_state(state, h, G, CS%dye_tracer_CSp) + call dye_tracer_surface_state(sfc_state, h, G, CS%dye_tracer_CSp) if (CS%use_oil) & - call oil_tracer_surface_state(state, h, G, CS%oil_tracer_CSp) + call oil_tracer_surface_state(sfc_state, h, G, CS%oil_tracer_CSp) if (CS%use_advection_test_tracer) & - call advection_test_tracer_surface_state(state, h, G, CS%advection_test_tracer_CSp) + call advection_test_tracer_surface_state(sfc_state, h, G, CS%advection_test_tracer_CSp) if (CS%use_OCMIP2_CFC) & - call OCMIP2_CFC_surface_state(state, h, G, CS%OCMIP2_CFC_CSp) + call OCMIP2_CFC_surface_state(sfc_state, h, G, CS%OCMIP2_CFC_CSp) #ifdef _USE_GENERIC_TRACER if (CS%use_MOM_generic_tracer) & - call MOM_generic_tracer_surface_state(state, h, G, CS%MOM_generic_tracer_CSp) + call MOM_generic_tracer_surface_state(sfc_state, h, G, CS%MOM_generic_tracer_CSp) #endif end subroutine call_tracer_surface_state diff --git a/src/tracer/advection_test_tracer.F90 b/src/tracer/advection_test_tracer.F90 index 82ea38f22c..b1d657d6e2 100644 --- a/src/tracer/advection_test_tracer.F90 +++ b/src/tracer/advection_test_tracer.F90 @@ -316,9 +316,9 @@ end subroutine advection_test_tracer_column_physics !> This subroutine extracts the surface fields from this tracer package that !! are to be shared with the atmosphere in coupled configurations. !! This particular tracer package does not report anything back to the coupler. -subroutine advection_test_tracer_surface_state(state, h, G, CS) +subroutine advection_test_tracer_surface_state(sfc_state, h, G, CS) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - type(surface), intent(inout) :: state !< A structure containing fields that + type(surface), intent(inout) :: sfc_state !< A structure containing fields that !! describe the surface state of the ocean. real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & intent(in) :: h !< Layer thickness [H ~> m or kg m-2]. @@ -339,7 +339,7 @@ subroutine advection_test_tracer_surface_state(state, h, G, CS) ! This call loads the surface values into the appropriate array in the ! coupler-type structure. call coupler_type_set_data(CS%tr(:,:,1,m), CS%ind_tr(m), ind_csurf, & - state%tr_fields, idim=(/isd, is, ie, ied/), & + sfc_state%tr_fields, idim=(/isd, is, ie, ied/), & jdim=(/jsd, js, je, jed/) ) enddo endif diff --git a/src/tracer/boundary_impulse_tracer.F90 b/src/tracer/boundary_impulse_tracer.F90 index e70320a5c7..da76cb3026 100644 --- a/src/tracer/boundary_impulse_tracer.F90 +++ b/src/tracer/boundary_impulse_tracer.F90 @@ -334,9 +334,9 @@ end function boundary_impulse_stock !> This subroutine extracts the surface fields from this tracer package that !! are to be shared with the atmosphere in coupled configurations. !! This particular tracer package does not report anything back to the coupler. -subroutine boundary_impulse_tracer_surface_state(state, h, G, CS) +subroutine boundary_impulse_tracer_surface_state(sfc_state, h, G, CS) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - type(surface), intent(inout) :: state !< A structure containing fields that + type(surface), intent(inout) :: sfc_state !< A structure containing fields that !! describe the surface state of the ocean. real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & intent(in) :: h !< Layer thickness [H ~> m or kg m-2]. @@ -357,7 +357,7 @@ subroutine boundary_impulse_tracer_surface_state(state, h, G, CS) ! This call loads the surface values into the appropriate array in the ! coupler-type structure. call coupler_type_set_data(CS%tr(:,:,1,m), CS%ind_tr(m), ind_csurf, & - state%tr_fields, idim=(/isd, is, ie, ied/), & + sfc_state%tr_fields, idim=(/isd, is, ie, ied/), & jdim=(/jsd, js, je, jed/) ) enddo endif diff --git a/src/tracer/dye_example.F90 b/src/tracer/dye_example.F90 index 86a4ac7aeb..5f2f139899 100644 --- a/src/tracer/dye_example.F90 +++ b/src/tracer/dye_example.F90 @@ -372,9 +372,9 @@ end function dye_stock !> This subroutine extracts the surface fields from this tracer package that !! are to be shared with the atmosphere in coupled configurations. !! This particular tracer package does not report anything back to the coupler. -subroutine dye_tracer_surface_state(state, h, G, CS) +subroutine dye_tracer_surface_state(sfc_state, h, G, CS) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - type(surface), intent(inout) :: state !< A structure containing fields that + type(surface), intent(inout) :: sfc_state !< A structure containing fields that !! describe the surface state of the ocean. real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & intent(in) :: h !< Layer thickness [H ~> m or kg m-2]. @@ -395,7 +395,7 @@ subroutine dye_tracer_surface_state(state, h, G, CS) ! This call loads the surface values into the appropriate array in the ! coupler-type structure. call coupler_type_set_data(CS%tr(:,:,1,m), CS%ind_tr(m), ind_csurf, & - state%tr_fields, idim=(/isd, is, ie, ied/), & + sfc_state%tr_fields, idim=(/isd, is, ie, ied/), & jdim=(/jsd, js, je, jed/) ) enddo endif diff --git a/src/tracer/ideal_age_example.F90 b/src/tracer/ideal_age_example.F90 index 3ef61e1a57..8f00b0d5b9 100644 --- a/src/tracer/ideal_age_example.F90 +++ b/src/tracer/ideal_age_example.F90 @@ -420,9 +420,9 @@ end function ideal_age_stock !> This subroutine extracts the surface fields from this tracer package that !! are to be shared with the atmosphere in coupled configurations. !! This particular tracer package does not report anything back to the coupler. -subroutine ideal_age_tracer_surface_state(state, h, G, CS) +subroutine ideal_age_tracer_surface_state(sfc_state, h, G, CS) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - type(surface), intent(inout) :: state !< A structure containing fields that + type(surface), intent(inout) :: sfc_state !< A structure containing fields that !! describe the surface state of the ocean. real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & intent(in) :: h !< Layer thickness [H ~> m or kg m-2]. @@ -443,7 +443,7 @@ subroutine ideal_age_tracer_surface_state(state, h, G, CS) ! This call loads the surface values into the appropriate array in the ! coupler-type structure. call coupler_type_set_data(CS%tr(:,:,1,m), CS%ind_tr(m), ind_csurf, & - state%tr_fields, idim=(/isd, is, ie, ied/), & + sfc_state%tr_fields, idim=(/isd, is, ie, ied/), & jdim=(/jsd, js, je, jed/) ) enddo endif diff --git a/src/tracer/oil_tracer.F90 b/src/tracer/oil_tracer.F90 index 4d755497c6..c07f1c03e4 100644 --- a/src/tracer/oil_tracer.F90 +++ b/src/tracer/oil_tracer.F90 @@ -454,9 +454,9 @@ end function oil_stock !> This subroutine extracts the surface fields from this tracer package that !! are to be shared with the atmosphere in coupled configurations. !! This particular tracer package does not report anything back to the coupler. -subroutine oil_tracer_surface_state(state, h, G, CS) +subroutine oil_tracer_surface_state(sfc_state, h, G, CS) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - type(surface), intent(inout) :: state !< A structure containing fields that + type(surface), intent(inout) :: sfc_state !< A structure containing fields that !! describe the surface state of the ocean. real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & intent(in) :: h !< Layer thickness [H ~> m or kg m-2]. @@ -477,7 +477,7 @@ subroutine oil_tracer_surface_state(state, h, G, CS) ! This call loads the surface values into the appropriate array in the ! coupler-type structure. call coupler_type_set_data(CS%tr(:,:,1,m), CS%ind_tr(m), ind_csurf, & - state%tr_fields, idim=(/isd, is, ie, ied/), & + sfc_state%tr_fields, idim=(/isd, is, ie, ied/), & jdim=(/jsd, js, je, jed/) ) enddo endif diff --git a/src/tracer/pseudo_salt_tracer.F90 b/src/tracer/pseudo_salt_tracer.F90 index 5c74487c0c..95396a3b58 100644 --- a/src/tracer/pseudo_salt_tracer.F90 +++ b/src/tracer/pseudo_salt_tracer.F90 @@ -299,9 +299,9 @@ end function pseudo_salt_stock !> This subroutine extracts the surface fields from this tracer package that !! are to be shared with the atmosphere in coupled configurations. !! This particular tracer package does not report anything back to the coupler. -subroutine pseudo_salt_tracer_surface_state(state, h, G, CS) +subroutine pseudo_salt_tracer_surface_state(sfc_state, h, G, CS) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - type(surface), intent(inout) :: state !< A structure containing fields that + type(surface), intent(inout) :: sfc_state !< A structure containing fields that !! describe the surface state of the ocean. real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & intent(in) :: h !< Layer thickness [H ~> m or kg m-2]. diff --git a/src/tracer/tracer_example.F90 b/src/tracer/tracer_example.F90 index c5e8f669c6..ef16cc985d 100644 --- a/src/tracer/tracer_example.F90 +++ b/src/tracer/tracer_example.F90 @@ -405,9 +405,9 @@ end function USER_tracer_stock !> This subroutine extracts the surface fields from this tracer package that !! are to be shared with the atmosphere in coupled configurations. -subroutine USER_tracer_surface_state(state, h, G, CS) +subroutine USER_tracer_surface_state(sfc_state, h, G, CS) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure - type(surface), intent(inout) :: state !< A structure containing fields that + type(surface), intent(inout) :: sfc_state !< A structure containing fields that !! describe the surface state of the ocean. real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] @@ -428,7 +428,7 @@ subroutine USER_tracer_surface_state(state, h, G, CS) ! This call loads the surface values into the appropriate array in the ! coupler-type structure. call coupler_type_set_data(CS%tr(:,:,1,m), CS%ind_tr(m), ind_csurf, & - state%tr_fields, idim=(/isd, is, ie, ied/), & + sfc_state%tr_fields, idim=(/isd, is, ie, ied/), & jdim=(/jsd, js, je, jed/) ) enddo endif diff --git a/src/user/BFB_surface_forcing.F90 b/src/user/BFB_surface_forcing.F90 index 70d89497da..a6aae3d3f7 100644 --- a/src/user/BFB_surface_forcing.F90 +++ b/src/user/BFB_surface_forcing.F90 @@ -47,8 +47,8 @@ module BFB_surface_forcing contains !> Bouyancy forcing for the boundary-forced-basin (BFB) configuration -subroutine BFB_buoyancy_forcing(state, fluxes, day, dt, G, US, CS) - type(surface), intent(inout) :: state !< A structure containing fields that +subroutine BFB_buoyancy_forcing(sfc_state, fluxes, day, dt, G, US, CS) + type(surface), intent(inout) :: sfc_state !< A structure containing fields that !! describe the surface state of the ocean. type(forcing), intent(inout) :: fluxes !< A structure containing pointers to any !! possible forcing fields. Unused fields @@ -136,9 +136,9 @@ subroutine BFB_buoyancy_forcing(state, fluxes, day, dt, G, US, CS) Salin_restore = 0.0 fluxes%heat_added(i,j) = (G%mask2dT(i,j) * (rhoXcp * CS%Flux_const)) * & - (Temp_restore - state%SST(i,j)) + (Temp_restore - sfc_state%SST(i,j)) fluxes%vprec(i,j) = - (G%mask2dT(i,j) * (CS%Rho0*CS%Flux_const)) * & - ((Salin_restore - state%SSS(i,j)) / (0.5 * (Salin_restore + state%SSS(i,j)))) + ((Salin_restore - sfc_state%SSS(i,j)) / (0.5 * (Salin_restore + sfc_state%SSS(i,j)))) enddo ; enddo else ! When modifying the code, comment out this error message. It is here @@ -164,7 +164,7 @@ subroutine BFB_buoyancy_forcing(state, fluxes, day, dt, G, US, CS) density_restore = Temp_restore*CS%drho_dt + CS%Rho0 fluxes%buoy(i,j) = G%mask2dT(i,j) * buoy_rest_const * & - (density_restore - US%kg_m3_to_R*state%sfc_density(i,j)) + (density_restore - US%kg_m3_to_R*sfc_state%sfc_density(i,j)) enddo ; enddo endif endif ! end RESTOREBUOY @@ -195,8 +195,7 @@ subroutine BFB_surface_forcing_init(Time, G, US, param_file, diag, CS) ! Read all relevant parameters and write them to the model log. call log_version(param_file, mdl, version, "") call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", CS%use_temperature, & - "If true, Temperature and salinity are used as state "//& - "variables.", default=.true.) + "If true, Temperature and salinity are used as state variables.", default=.true.) call get_param(param_file, mdl, "G_EARTH", CS%G_Earth, & "The gravitational acceleration of the Earth.", & diff --git a/src/user/Idealized_Hurricane.F90 b/src/user/Idealized_Hurricane.F90 index 38ba0ab460..5727309926 100644 --- a/src/user/Idealized_Hurricane.F90 +++ b/src/user/Idealized_Hurricane.F90 @@ -206,8 +206,8 @@ subroutine idealized_hurricane_wind_init(Time, G, US, param_file, CS) end subroutine idealized_hurricane_wind_init !> Computes the surface wind for the idealized hurricane test cases -subroutine idealized_hurricane_wind_forcing(state, forces, day, G, US, CS) - type(surface), intent(in) :: state !< Surface state structure +subroutine idealized_hurricane_wind_forcing(sfc_state, forces, day, G, US, CS) + type(surface), intent(in) :: sfc_state !< Surface state structure type(mech_forcing), intent(inout) :: forces !< A structure with the driving mechanical forces type(time_type), intent(in) :: day !< Time in days type(ocean_grid_type), intent(inout) :: G !< Grid structure @@ -263,13 +263,13 @@ subroutine idealized_hurricane_wind_forcing(state, forces, day, G, US, CS) !> Computes taux do j=js,je do I=is-1,Ieq - Uocn = US%m_s_to_L_T * state%u(I,j)*REL_TAU_FAC + Uocn = US%m_s_to_L_T * sfc_state%u(I,j)*REL_TAU_FAC if (CS%answers_2018) then - Vocn = US%m_s_to_L_T * 0.25*(state%v(i,J)+state%v(i+1,J-1)& - +state%v(i+1,J)+state%v(i,J-1))*REL_TAU_FAC + Vocn = US%m_s_to_L_T * 0.25*(sfc_state%v(i,J)+sfc_state%v(i+1,J-1)& + +sfc_state%v(i+1,J)+sfc_state%v(i,J-1))*REL_TAU_FAC else - Vocn = US%m_s_to_L_T * 0.25*((state%v(i,J)+state%v(i+1,J-1)) +& - (state%v(i+1,J)+state%v(i,J-1))) * REL_TAU_FAC + Vocn = US%m_s_to_L_T * 0.25*((sfc_state%v(i,J)+sfc_state%v(i+1,J-1)) +& + (sfc_state%v(i+1,J)+sfc_state%v(i,J-1))) * REL_TAU_FAC endif f_local = abs(0.5*(G%CoriolisBu(I,J)+G%CoriolisBu(I,J-1)))*fbench_fac + fbench ! Calculate position as a function of time. @@ -288,13 +288,13 @@ subroutine idealized_hurricane_wind_forcing(state, forces, day, G, US, CS) do J=js-1,Jeq do i=is,ie if (CS%answers_2018) then - Uocn = US%m_s_to_L_T * 0.25*(state%u(I,j)+state%u(I-1,j+1) + & - state%u(I-1,j)+state%u(I,j+1))*REL_TAU_FAC + Uocn = US%m_s_to_L_T * 0.25*(sfc_state%u(I,j)+sfc_state%u(I-1,j+1) + & + sfc_state%u(I-1,j)+sfc_state%u(I,j+1))*REL_TAU_FAC else - Uocn = US%m_s_to_L_T * 0.25*((state%u(I,j)+state%u(I-1,j+1)) + & - (state%u(I-1,j)+state%u(I,j+1))) * REL_TAU_FAC + Uocn = US%m_s_to_L_T * 0.25*((sfc_state%u(I,j)+sfc_state%u(I-1,j+1)) + & + (sfc_state%u(I-1,j)+sfc_state%u(I,j+1))) * REL_TAU_FAC endif - Vocn = US%m_s_to_L_T * state%v(i,J)*REL_TAU_FAC + Vocn = US%m_s_to_L_T * sfc_state%v(i,J)*REL_TAU_FAC f_local = abs(0.5*(G%CoriolisBu(I-1,J)+G%CoriolisBu(I,J)))*fbench_fac + fbench ! Calculate position as a function of time. if (CS%SCM_mode) then @@ -471,8 +471,8 @@ end subroutine idealized_hurricane_wind_profile !! It is included as an additional subroutine rather than padded into the previous !! routine with flags to ease its eventual removal. Its functionality is replaced !! with the new routines and it can be deleted when answer changes are acceptable. -subroutine SCM_idealized_hurricane_wind_forcing(state, forces, day, G, US, CS) - type(surface), intent(in) :: state !< Surface state structure +subroutine SCM_idealized_hurricane_wind_forcing(sfc_state, forces, day, G, US, CS) + type(surface), intent(in) :: sfc_state !< Surface state structure type(mech_forcing), intent(inout) :: forces !< A structure with the driving mechanical forces type(time_type), intent(in) :: day !< Time in days type(ocean_grid_type), intent(inout) :: G !< Grid structure @@ -604,9 +604,9 @@ subroutine SCM_idealized_hurricane_wind_forcing(state, forces, day, G, US, CS) !/BR ! Turn off surface current for stress calculation to be ! consistent with test case. - Uocn = 0. ! state%u(I,j) - Vocn = 0. ! 0.25*( (state%v(i,J) + state%v(i+1,J-1)) & - ! +(state%v(i+1,J) + state%v(i,J-1)) ) + Uocn = 0. ! sfc_state%u(I,j) + Vocn = 0. ! 0.25*( (sfc_state%v(i,J) + sfc_state%v(i+1,J-1)) + & + ! (sfc_state%v(i+1,J) + sfc_state%v(i,J-1)) ) !/BR ! Wind vector calculated from location/direction (sin/cos flipped b/c ! cyclonic wind is 90 deg. phase shifted from position angle). @@ -633,9 +633,9 @@ subroutine SCM_idealized_hurricane_wind_forcing(state, forces, day, G, US, CS) !/BR ! See notes above do J=js-1,Jeq ; do i=is,ie - Uocn = 0. ! 0.25*( (state%u(I,j) + state%u(I-1,j+1)) & - ! +(state%u(I-1,j) + state%u(I,j+1)) ) - Vocn = 0. ! state%v(i,J) + Uocn = 0. ! 0.25*( (sfc_state%u(I,j) + sfc_state%u(I-1,j+1)) + & + ! (sfc_state%u(I-1,j) + sfc_state%u(I,j+1)) ) + Vocn = 0. ! sfc_state%v(i,J) dU = U10*sin(Adir-pie-Alph) - Uocn + U_TS dV = U10*cos(Adir-Alph) - Vocn + V_TS du10=sqrt(du**2+dv**2) diff --git a/src/user/SCM_CVMix_tests.F90 b/src/user/SCM_CVMix_tests.F90 index a63205fede..1bb1b9555e 100644 --- a/src/user/SCM_CVMix_tests.F90 +++ b/src/user/SCM_CVMix_tests.F90 @@ -200,8 +200,8 @@ subroutine SCM_CVMix_tests_surface_forcing_init(Time, G, param_file, CS) end subroutine SCM_CVMix_tests_surface_forcing_init -subroutine SCM_CVMix_tests_wind_forcing(state, forces, day, G, US, CS) - type(surface), intent(in) :: state !< Surface state structure +subroutine SCM_CVMix_tests_wind_forcing(sfc_state, forces, day, G, US, CS) + type(surface), intent(in) :: sfc_state !< Surface state structure type(mech_forcing), intent(inout) :: forces !< A structure with the driving mechanical forces type(time_type), intent(in) :: day !< Time in days type(ocean_grid_type), intent(inout) :: G !< Grid structure @@ -233,8 +233,8 @@ subroutine SCM_CVMix_tests_wind_forcing(state, forces, day, G, US, CS) end subroutine SCM_CVMix_tests_wind_forcing -subroutine SCM_CVMix_tests_buoyancy_forcing(state, fluxes, day, G, US, CS) - type(surface), intent(in) :: state !< Surface state structure +subroutine SCM_CVMix_tests_buoyancy_forcing(sfc_state, fluxes, day, G, US, CS) + type(surface), intent(in) :: sfc_state !< Surface state structure type(forcing), intent(inout) :: fluxes !< Surface fluxes structure type(time_type), intent(in) :: day !< Current model time type(ocean_grid_type), intent(inout) :: G !< Grid structure diff --git a/src/user/dumbbell_surface_forcing.F90 b/src/user/dumbbell_surface_forcing.F90 index 5be2bc9b8e..4c582dd03e 100644 --- a/src/user/dumbbell_surface_forcing.F90 +++ b/src/user/dumbbell_surface_forcing.F90 @@ -24,8 +24,7 @@ module dumbbell_surface_forcing !> Control structure for the dumbbell test case forcing type, public :: dumbbell_surface_forcing_CS ; private - logical :: use_temperature !< If true, temperature and salinity are used as - !! state variables. + logical :: use_temperature !< If true, temperature and salinity are used as state variables. logical :: restorebuoy !< If true, use restoring surface buoyancy forcing. real :: Rho0 !< The density used in the Boussinesq approximation [R ~> kg m-3]. real :: G_Earth !< The gravitational acceleration [L2 Z-1 T-2 ~> m s-2] @@ -46,8 +45,8 @@ module dumbbell_surface_forcing contains !> Surface buoyancy (heat and fresh water) fluxes for the dumbbell test case -subroutine dumbbell_buoyancy_forcing(state, fluxes, day, dt, G, US, CS) - type(surface), intent(inout) :: state !< A structure containing fields that +subroutine dumbbell_buoyancy_forcing(sfc_state, fluxes, day, dt, G, US, CS) + type(surface), intent(inout) :: sfc_state !< A structure containing fields that !! describe the surface state of the ocean. type(forcing), intent(inout) :: fluxes !< A structure containing pointers to any !! possible forcing fields. Unused fields @@ -119,7 +118,7 @@ subroutine dumbbell_buoyancy_forcing(state, fluxes, day, dt, G, US, CS) do j=js,je ; do i=is,ie if (CS%forcing_mask(i,j)>0.) then fluxes%vprec(i,j) = - (G%mask2dT(i,j) * (CS%Rho0*CS%Flux_const)) * & - ((CS%S_restore(i,j) - state%SSS(i,j)) / (0.5 * (CS%S_restore(i,j) + state%SSS(i,j)))) + ((CS%S_restore(i,j) - sfc_state%SSS(i,j)) / (0.5 * (CS%S_restore(i,j) + sfc_state%SSS(i,j)))) endif enddo ; enddo @@ -128,8 +127,8 @@ subroutine dumbbell_buoyancy_forcing(state, fluxes, day, dt, G, US, CS) end subroutine dumbbell_buoyancy_forcing !> Dynamic forcing for the dumbbell test case -subroutine dumbbell_dynamic_forcing(state, fluxes, day, dt, G, CS) - type(surface), intent(inout) :: state !< A structure containing fields that +subroutine dumbbell_dynamic_forcing(sfc_state, fluxes, day, dt, G, CS) + type(surface), intent(inout) :: sfc_state !< A structure containing fields that !! describe the surface state of the ocean. type(forcing), intent(inout) :: fluxes !< A structure containing pointers to any !! possible forcing fields. Unused fields @@ -198,8 +197,7 @@ subroutine dumbbell_surface_forcing_init(Time, G, US, param_file, diag, CS) ! Read all relevant parameters and write them to the model log. call log_version(param_file, mdl, version, "") call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", CS%use_temperature, & - "If true, Temperature and salinity are used as state "//& - "variables.", default=.true.) + "If true, Temperature and salinity are used as state variables.", default=.true.) call get_param(param_file, mdl, "G_EARTH", CS%G_Earth, & "The gravitational acceleration of the Earth.", & diff --git a/src/user/user_revise_forcing.F90 b/src/user/user_revise_forcing.F90 index d1be729734..c53451f4e8 100644 --- a/src/user/user_revise_forcing.F90 +++ b/src/user/user_revise_forcing.F90 @@ -30,8 +30,8 @@ module user_revise_forcing contains !> This subroutine sets the surface wind stresses. -subroutine user_alter_forcing(state, fluxes, day, G, CS) - type(surface), intent(in) :: state !< A structure containing fields that +subroutine user_alter_forcing(sfc_state, fluxes, day, G, CS) + type(surface), intent(in) :: sfc_state !< A structure containing fields that !! describe the surface state of the ocean. type(forcing), intent(inout) :: fluxes !< A structure containing pointers to any !! possible forcing fields. Unused fields From c6321d642e1781ad7b527d490af4749b26a6d36c Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 30 Apr 2020 08:21:35 -0400 Subject: [PATCH 05/11] +Rescaled sfc_state%u and sfc_state%v Dimensionally rescaled the surface velocity variables to [L T-1]. All answers are bitwise identical, but the units of two elements in a transparent type have changed. --- config_src/coupled_driver/ocean_model_MOM.F90 | 12 +++++------ config_src/mct_driver/mom_ocean_model_mct.F90 | 12 +++++------ .../nuopc_driver/mom_ocean_model_nuopc.F90 | 12 +++++------ src/core/MOM.F90 | 20 +++++++++---------- src/core/MOM_checksum_packages.F90 | 3 ++- src/core/MOM_variables.F90 | 4 ++-- src/diagnostics/MOM_diagnostics.F90 | 14 ++++++------- src/ice_shelf/MOM_ice_shelf.F90 | 12 +++++------ src/user/Idealized_Hurricane.F90 | 20 +++++++++---------- 9 files changed, 55 insertions(+), 54 deletions(-) diff --git a/config_src/coupled_driver/ocean_model_MOM.F90 b/config_src/coupled_driver/ocean_model_MOM.F90 index 407a11a0c3..28ac193d8d 100644 --- a/config_src/coupled_driver/ocean_model_MOM.F90 +++ b/config_src/coupled_driver/ocean_model_MOM.F90 @@ -880,22 +880,22 @@ subroutine convert_state_to_ocean_type(sfc_state, Ocean_sfc, G, US, patm, press_ if (Ocean_sfc%stagger == AGRID) then do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd - Ocean_sfc%u_surf(i,j) = G%mask2dT(i+i0,j+j0) * & + Ocean_sfc%u_surf(i,j) = G%mask2dT(i+i0,j+j0) * US%L_T_to_m_s * & 0.5*(sfc_state%u(I+i0,j+j0)+sfc_state%u(I-1+i0,j+j0)) - Ocean_sfc%v_surf(i,j) = G%mask2dT(i+i0,j+j0) * & + Ocean_sfc%v_surf(i,j) = G%mask2dT(i+i0,j+j0) * US%L_T_to_m_s * & 0.5*(sfc_state%v(i+i0,J+j0)+sfc_state%v(i+i0,J-1+j0)) enddo ; enddo elseif (Ocean_sfc%stagger == BGRID_NE) then do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd - Ocean_sfc%u_surf(i,j) = G%mask2dBu(I+i0,J+j0) * & + Ocean_sfc%u_surf(i,j) = G%mask2dBu(I+i0,J+j0) * US%L_T_to_m_s * & 0.5*(sfc_state%u(I+i0,j+j0)+sfc_state%u(I+i0,j+j0+1)) - Ocean_sfc%v_surf(i,j) = G%mask2dBu(I+i0,J+j0) * & + Ocean_sfc%v_surf(i,j) = G%mask2dBu(I+i0,J+j0) * US%L_T_to_m_s * & 0.5*(sfc_state%v(i+i0,J+j0)+sfc_state%v(i+i0+1,J+j0)) enddo ; enddo elseif (Ocean_sfc%stagger == CGRID_NE) then do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd - Ocean_sfc%u_surf(i,j) = G%mask2dCu(I+i0,j+j0)*sfc_state%u(I+i0,j+j0) - Ocean_sfc%v_surf(i,j) = G%mask2dCv(i+i0,J+j0)*sfc_state%v(i+i0,J+j0) + Ocean_sfc%u_surf(i,j) = G%mask2dCu(I+i0,j+j0)*US%L_T_to_m_s * sfc_state%u(I+i0,j+j0) + Ocean_sfc%v_surf(i,j) = G%mask2dCv(i+i0,J+j0)*US%L_T_to_m_s * sfc_state%v(i+i0,J+j0) enddo ; enddo else write(val_str, '(I8)') Ocean_sfc%stagger diff --git a/config_src/mct_driver/mom_ocean_model_mct.F90 b/config_src/mct_driver/mom_ocean_model_mct.F90 index 317a496399..0ef8a7a9d4 100644 --- a/config_src/mct_driver/mom_ocean_model_mct.F90 +++ b/config_src/mct_driver/mom_ocean_model_mct.F90 @@ -924,22 +924,22 @@ subroutine convert_state_to_ocean_type(sfc_state, Ocean_sfc, G, US, patm, press_ if (Ocean_sfc%stagger == AGRID) then do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd - Ocean_sfc%u_surf(i,j) = G%mask2dT(i+i0,j+j0) * & + Ocean_sfc%u_surf(i,j) = G%mask2dT(i+i0,j+j0) * US%L_T_to_m_s * & 0.5*(sfc_state%u(I+i0,j+j0)+sfc_state%u(I-1+i0,j+j0)) - Ocean_sfc%v_surf(i,j) = G%mask2dT(i+i0,j+j0) * & + Ocean_sfc%v_surf(i,j) = G%mask2dT(i+i0,j+j0) * US%L_T_to_m_s * & 0.5*(sfc_state%v(i+i0,J+j0)+sfc_state%v(i+i0,J-1+j0)) enddo ; enddo elseif (Ocean_sfc%stagger == BGRID_NE) then do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd - Ocean_sfc%u_surf(i,j) = G%mask2dBu(I+i0,J+j0) * & + Ocean_sfc%u_surf(i,j) = G%mask2dBu(I+i0,J+j0) * US%L_T_to_m_s * & 0.5*(sfc_state%u(I+i0,j+j0)+sfc_state%u(I+i0,j+j0+1)) - Ocean_sfc%v_surf(i,j) = G%mask2dBu(I+i0,J+j0) * & + Ocean_sfc%v_surf(i,j) = G%mask2dBu(I+i0,J+j0) * US%L_T_to_m_s * & 0.5*(sfc_state%v(i+i0,J+j0)+sfc_state%v(i+i0+1,J+j0)) enddo ; enddo elseif (Ocean_sfc%stagger == CGRID_NE) then do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd - Ocean_sfc%u_surf(i,j) = G%mask2dCu(I+i0,j+j0)*sfc_state%u(I+i0,j+j0) - Ocean_sfc%v_surf(i,j) = G%mask2dCv(i+i0,J+j0)*sfc_state%v(i+i0,J+j0) + Ocean_sfc%u_surf(i,j) = G%mask2dCu(I+i0,j+j0) * US%L_T_to_m_s * sfc_state%u(I+i0,j+j0) + Ocean_sfc%v_surf(i,j) = G%mask2dCv(i+i0,J+j0) * US%L_T_to_m_s * sfc_state%v(i+i0,J+j0) enddo ; enddo else write(val_str, '(I8)') Ocean_sfc%stagger diff --git a/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 b/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 index ffbe73881a..1a08869b77 100644 --- a/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 +++ b/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 @@ -919,22 +919,22 @@ subroutine convert_state_to_ocean_type(sfc_state, Ocean_sfc, G, US, patm, press_ if (Ocean_sfc%stagger == AGRID) then do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd - Ocean_sfc%u_surf(i,j) = G%mask2dT(i+i0,j+j0) * & + Ocean_sfc%u_surf(i,j) = G%mask2dT(i+i0,j+j0) * US%L_T_to_m_s * & 0.5*(sfc_state%u(I+i0,j+j0)+sfc_state%u(I-1+i0,j+j0)) - Ocean_sfc%v_surf(i,j) = G%mask2dT(i+i0,j+j0) * & + Ocean_sfc%v_surf(i,j) = G%mask2dT(i+i0,j+j0) * US%L_T_to_m_s * & 0.5*(sfc_state%v(i+i0,J+j0)+sfc_state%v(i+i0,J-1+j0)) enddo ; enddo elseif (Ocean_sfc%stagger == BGRID_NE) then do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd - Ocean_sfc%u_surf(i,j) = G%mask2dBu(I+i0,J+j0) * & + Ocean_sfc%u_surf(i,j) = G%mask2dBu(I+i0,J+j0) * US%L_T_to_m_s * & 0.5*(sfc_state%u(I+i0,j+j0)+sfc_state%u(I+i0,j+j0+1)) - Ocean_sfc%v_surf(i,j) = G%mask2dBu(I+i0,J+j0) * & + Ocean_sfc%v_surf(i,j) = G%mask2dBu(I+i0,J+j0) * US%L_T_to_m_s * & 0.5*(sfc_state%v(i+i0,J+j0)+sfc_state%v(i+i0+1,J+j0)) enddo ; enddo elseif (Ocean_sfc%stagger == CGRID_NE) then do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd - Ocean_sfc%u_surf(i,j) = G%mask2dCu(I+i0,j+j0)*sfc_state%u(I+i0,j+j0) - Ocean_sfc%v_surf(i,j) = G%mask2dCv(i+i0,J+j0)*sfc_state%v(i+i0,J+j0) + Ocean_sfc%u_surf(i,j) = G%mask2dCu(I+i0,j+j0) * US%L_T_to_m_s * sfc_state%u(I+i0,j+j0) + Ocean_sfc%v_surf(i,j) = G%mask2dCv(i+i0,J+j0) * US%L_T_to_m_s * sfc_state%v(i+i0,J+j0) enddo ; enddo else write(val_str, '(I8)') Ocean_sfc%stagger diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 49eed5fe1f..c47115ed68 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -3029,10 +3029,10 @@ subroutine extract_surface_state(CS, sfc_state_in) sfc_state%SSS(i,j) = CS%tv%S(i,j,1) enddo ; enddo ; endif do j=js,je ; do I=is-1,ie - sfc_state%u(I,j) = US%L_T_to_m_s * CS%u(I,j,1) + sfc_state%u(I,j) = CS%u(I,j,1) enddo ; enddo do J=js-1,je ; do i=is,ie - sfc_state%v(i,J) = US%L_T_to_m_s * CS%v(i,J,1) + sfc_state%v(i,J) = CS%v(i,J,1) enddo ; enddo else ! (CS%Hmix >= 0.0) @@ -3125,7 +3125,7 @@ subroutine extract_surface_state(CS, sfc_state_in) else dh = 0.0 endif - sfc_state%v(i,J) = sfc_state%v(i,J) + dh * US%L_T_to_m_s * CS%v(i,J,k) + sfc_state%v(i,J) = sfc_state%v(i,J) + dh * CS%v(i,J,k) depth(i) = depth(i) + dh enddo ; enddo ! Calculate the average properties of the mixed layer depth. @@ -3149,7 +3149,7 @@ subroutine extract_surface_state(CS, sfc_state_in) else dh = 0.0 endif - sfc_state%u(I,j) = sfc_state%u(I,j) + dh * US%L_T_to_m_s * CS%u(I,j,k) + sfc_state%u(I,j) = sfc_state%u(I,j) + dh * CS%u(I,j,k) depth(I) = depth(I) + dh enddo ; enddo ! Calculate the average properties of the mixed layer depth. @@ -3159,10 +3159,10 @@ subroutine extract_surface_state(CS, sfc_state_in) enddo ! end of j loop else ! Hmix_UV<=0. do j=js,je ; do I=is-1,ie - sfc_state%u(I,j) = US%L_T_to_m_s * CS%u(I,j,1) + sfc_state%u(I,j) = CS%u(I,j,1) enddo ; enddo do J=js-1,je ; do i=is,ie - sfc_state%v(i,J) = US%L_T_to_m_s * CS%v(i,J,1) + sfc_state%v(i,J) = CS%v(i,J,1) enddo ; enddo endif endif ! (CS%Hmix >= 0.0) @@ -3311,16 +3311,16 @@ subroutine extract_surface_state(CS, sfc_state_in) 'x=',G%gridLonT(ig), 'y=',G%gridLatT(jg), & 'D=',bathy_m, 'SSH=',sfc_state%sea_lev(i,j), & 'SST=',sfc_state%SST(i,j), 'SSS=',sfc_state%SSS(i,j), & - 'U-=',sfc_state%u(I-1,j), 'U+=',sfc_state%u(I,j), & - 'V-=',sfc_state%v(i,J-1), 'V+=',sfc_state%v(i,J) + 'U-=',US%L_T_to_m_s*sfc_state%u(I-1,j), 'U+=',US%L_T_to_m_s*sfc_state%u(I,j), & + 'V-=',US%L_T_to_m_s*sfc_state%v(i,J-1), 'V+=',US%L_T_to_m_s*sfc_state%v(i,J) else write(msg(1:240),'(2(a,i4,x),4(a,f8.3,x),6(a,es11.4))') & 'Extreme surface sfc_state detected: i=',ig,'j=',jg, & 'lon=',G%geoLonT(i,j), 'lat=',G%geoLatT(i,j), & 'x=',G%gridLonT(i), 'y=',G%gridLatT(j), & 'D=',bathy_m, 'SSH=',sfc_state%sea_lev(i,j), & - 'U-=',sfc_state%u(I-1,j), 'U+=',sfc_state%u(I,j), & - 'V-=',sfc_state%v(i,J-1), 'V+=',sfc_state%v(i,J) + 'U-=',US%L_T_to_m_s*sfc_state%u(I-1,j), 'U+=',US%L_T_to_m_s*sfc_state%u(I,j), & + 'V-=',US%L_T_to_m_s*sfc_state%v(i,J-1), 'V+=',US%L_T_to_m_s*sfc_state%v(i,J) endif call MOM_error(WARNING, trim(msg), all_print=.true.) elseif (numberOfErrors==9) then ! Indicate once that there are more errors diff --git a/src/core/MOM_checksum_packages.F90 b/src/core/MOM_checksum_packages.F90 index dfe34676a6..c63b15a3b4 100644 --- a/src/core/MOM_checksum_packages.F90 +++ b/src/core/MOM_checksum_packages.F90 @@ -162,7 +162,8 @@ subroutine MOM_surface_chksum(mesg, sfc_state, G, haloshift, symmetric) if (allocated(sfc_state%sea_lev)) call hchksum(sfc_state%sea_lev, mesg//" sea_lev", G%HI, haloshift=hs) if (allocated(sfc_state%Hml)) call hchksum(sfc_state%Hml, mesg//" Hml", G%HI, haloshift=hs) 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) + call uvchksum(mesg//" SSU", sfc_state%u, sfc_state%v, G%HI, haloshift=hs, symmetric=sym, & + scale=G%US%L_T_to_m_s) ! if (allocated(sfc_state%salt_deficit)) & ! call hchksum(sfc_state%salt_deficit, mesg//" salt deficit", G%HI, haloshift=hs) if (allocated(sfc_state%frazil)) call hchksum(sfc_state%frazil, mesg//" frazil", G%HI, haloshift=hs) diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index 2ac62eee5a..2b7cafdd33 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -44,8 +44,8 @@ module MOM_variables SSS, & !< The sea surface salinity [ppt ~> psu or gSalt/kg]. sfc_density, & !< The mixed layer density [kg m-3]. Hml, & !< The mixed layer depth [m]. - u, & !< The mixed layer zonal velocity [m s-1]. - v, & !< The mixed layer meridional velocity [m s-1]. + u, & !< The mixed layer zonal velocity [L T-1 ~> m s-1]. + v, & !< The mixed layer meridional velocity [L T-1 ~> m s-1]. sea_lev, & !< The sea level [m]. If a reduced surface gravity is !! used, that is compensated for in sea_lev. frazil, & !< The energy needed to heat the ocean column to the freezing point during diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 77739f3ead..c5c47b743f 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -1169,7 +1169,7 @@ subroutine post_surface_dyn_diags(IDs, G, diag, sfc_state, ssh) intent(in) :: ssh !< Time mean surface height without corrections for ice displacement [m] ! Local variables - real, dimension(SZI_(G),SZJ_(G)) :: work_2d ! A 2-d work array + real, dimension(SZI_(G),SZJ_(G)) :: speed ! The surface speed [L T-1 ~> m s-1] integer :: i, j, is, ie, js, je is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec @@ -1185,10 +1185,10 @@ subroutine post_surface_dyn_diags(IDs, G, diag, sfc_state, ssh) if (IDs%id_speed > 0) then do j=js,je ; do i=is,ie - work_2d(i,j) = sqrt(0.5*(sfc_state%u(I-1,j)**2 + sfc_state%u(I,j)**2) + & - 0.5*(sfc_state%v(i,J-1)**2 + sfc_state%v(i,J)**2)) + speed(i,j) = sqrt(0.5*(sfc_state%u(I-1,j)**2 + sfc_state%u(I,j)**2) + & + 0.5*(sfc_state%v(i,J-1)**2 + sfc_state%v(i,J)**2)) enddo ; enddo - call post_data(IDs%id_speed, work_2d, diag, mask=G%mask2dT) + call post_data(IDs%id_speed, speed, diag, mask=G%mask2dT) endif end subroutine post_surface_dyn_diags @@ -1784,11 +1784,11 @@ subroutine register_surface_diags(Time, G, US, IDs, diag, tv) long_name='Area averaged sea surface height', units='m', & standard_name='area_averaged_sea_surface_height') IDs%id_ssu = register_diag_field('ocean_model', 'SSU', diag%axesCu1, Time, & - 'Sea Surface Zonal Velocity', 'm s-1') + 'Sea Surface Zonal Velocity', 'm s-1', conversion=US%L_T_to_m_s) IDs%id_ssv = register_diag_field('ocean_model', 'SSV', diag%axesCv1, Time, & - 'Sea Surface Meridional Velocity', 'm s-1') + 'Sea Surface Meridional Velocity', 'm s-1', conversion=US%L_T_to_m_s) IDs%id_speed = register_diag_field('ocean_model', 'speed', diag%axesT1, Time, & - 'Sea Surface Speed', 'm s-1') + 'Sea Surface Speed', 'm s-1', conversion=US%L_T_to_m_s) if (associated(tv%T)) then IDs%id_sst = register_diag_field('ocean_model', 'SST', diag%axesT1, Time, & diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 819ce8bb76..3e8b1ed83b 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -334,8 +334,8 @@ subroutine shelf_calc_flux(sfc_state, fluxes, Time, time_step, CS, forces) call hchksum(fluxes%frac_shelf_h, "frac_shelf_h before apply melting", G%HI, haloshift=0) call hchksum(sfc_state%sst, "sst before apply melting", G%HI, haloshift=0) call hchksum(sfc_state%sss, "sss before apply melting", G%HI, haloshift=0) - call hchksum(sfc_state%u, "u_ml before apply melting", G%HI, haloshift=0) - call hchksum(sfc_state%v, "v_ml before apply melting", G%HI, haloshift=0) + call hchksum(sfc_state%u, "u_ml before apply melting", G%HI, haloshift=0, scale=US%L_T_to_m_s) + call hchksum(sfc_state%v, "v_ml before apply melting", G%HI, haloshift=0, scale=US%L_T_to_m_s) call hchksum(sfc_state%ocean_mass, "ocean_mass before apply melting", G%HI, haloshift=0) endif @@ -356,8 +356,8 @@ subroutine shelf_calc_flux(sfc_state, fluxes, Time, time_step, CS, forces) taux2 = (asu1 * sfc_state%taux_shelf(I-1,j)**2 + asu2 * sfc_state%taux_shelf(I,j)**2 ) * I_au tauy2 = (asv1 * sfc_state%tauy_shelf(i,J-1)**2 + asv2 * sfc_state%tauy_shelf(i,J)**2 ) * I_av endif - u2_av = US%m_s_to_L_T**2*(asu1 * sfc_state%u(I-1,j)**2 + asu2 * sfc_state%u(I,j)**2) * I_au - v2_av = US%m_s_to_L_T**2*(asv1 * sfc_state%v(i,J-1)**2 + asu2 * sfc_state%v(i,J)**2) * I_av + u2_av = (asu1 * sfc_state%u(I-1,j)**2 + asu2 * sfc_state%u(I,j)**2) * I_au + v2_av = (asv1 * sfc_state%v(i,J-1)**2 + asu2 * sfc_state%v(i,J)**2) * I_av if (taux2 + tauy2 > 0.0) then fluxes%ustar_shelf(i,j) = MAX(CS%ustar_bg, US%L_to_Z * & @@ -1617,9 +1617,9 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl CS%id_Sbdry = register_diag_field('ocean_model', 'sbdry', CS%diag%axesT1, CS%Time, & 'salinity at the ice-ocean interface.', 'psu') CS%id_u_ml = register_diag_field('ocean_model', 'u_ml', CS%diag%axesCu1, CS%Time, & - 'Eastward vel. in the boundary layer (used to compute ustar)', 'm s-1') + 'Eastward vel. in the boundary layer (used to compute ustar)', 'm s-1', conversion=US%L_T_to_m_s) CS%id_v_ml = register_diag_field('ocean_model', 'v_ml', CS%diag%axesCv1, CS%Time, & - 'Northward vel. in the boundary layer (used to compute ustar)', 'm s-1') + 'Northward vel. in the boundary layer (used to compute ustar)', 'm s-1', conversion=US%L_T_to_m_s) CS%id_exch_vel_s = register_diag_field('ocean_model', 'exch_vel_s', CS%diag%axesT1, CS%Time, & 'Sub-shelf salinity exchange velocity', 'm s-1', conversion=US%Z_to_m*US%s_to_T) CS%id_exch_vel_t = register_diag_field('ocean_model', 'exch_vel_t', CS%diag%axesT1, CS%Time, & diff --git a/src/user/Idealized_Hurricane.F90 b/src/user/Idealized_Hurricane.F90 index 5727309926..a8ec1d06ff 100644 --- a/src/user/Idealized_Hurricane.F90 +++ b/src/user/Idealized_Hurricane.F90 @@ -263,13 +263,13 @@ subroutine idealized_hurricane_wind_forcing(sfc_state, forces, day, G, US, CS) !> Computes taux do j=js,je do I=is-1,Ieq - Uocn = US%m_s_to_L_T * sfc_state%u(I,j)*REL_TAU_FAC + Uocn = sfc_state%u(I,j) * REL_TAU_FAC if (CS%answers_2018) then - Vocn = US%m_s_to_L_T * 0.25*(sfc_state%v(i,J)+sfc_state%v(i+1,J-1)& - +sfc_state%v(i+1,J)+sfc_state%v(i,J-1))*REL_TAU_FAC + Vocn = 0.25*(sfc_state%v(i,J)+sfc_state%v(i+1,J-1)& + +sfc_state%v(i+1,J)+sfc_state%v(i,J-1))*REL_TAU_FAC else - Vocn = US%m_s_to_L_T * 0.25*((sfc_state%v(i,J)+sfc_state%v(i+1,J-1)) +& - (sfc_state%v(i+1,J)+sfc_state%v(i,J-1))) * REL_TAU_FAC + Vocn =0.25*((sfc_state%v(i,J)+sfc_state%v(i+1,J-1)) +& + (sfc_state%v(i+1,J)+sfc_state%v(i,J-1))) * REL_TAU_FAC endif f_local = abs(0.5*(G%CoriolisBu(I,J)+G%CoriolisBu(I,J-1)))*fbench_fac + fbench ! Calculate position as a function of time. @@ -288,13 +288,13 @@ subroutine idealized_hurricane_wind_forcing(sfc_state, forces, day, G, US, CS) do J=js-1,Jeq do i=is,ie if (CS%answers_2018) then - Uocn = US%m_s_to_L_T * 0.25*(sfc_state%u(I,j)+sfc_state%u(I-1,j+1) + & - sfc_state%u(I-1,j)+sfc_state%u(I,j+1))*REL_TAU_FAC + Uocn = 0.25*(sfc_state%u(I,j)+sfc_state%u(I-1,j+1) + & + sfc_state%u(I-1,j)+sfc_state%u(I,j+1))*REL_TAU_FAC else - Uocn = US%m_s_to_L_T * 0.25*((sfc_state%u(I,j)+sfc_state%u(I-1,j+1)) + & - (sfc_state%u(I-1,j)+sfc_state%u(I,j+1))) * REL_TAU_FAC + Uocn = 0.25*((sfc_state%u(I,j)+sfc_state%u(I-1,j+1)) + & + (sfc_state%u(I-1,j)+sfc_state%u(I,j+1))) * REL_TAU_FAC endif - Vocn = US%m_s_to_L_T * sfc_state%v(i,J)*REL_TAU_FAC + Vocn = sfc_state%v(i,J) * REL_TAU_FAC f_local = abs(0.5*(G%CoriolisBu(I-1,J)+G%CoriolisBu(I,J)))*fbench_fac + fbench ! Calculate position as a function of time. if (CS%SCM_mode) then From 4b8fffcb73316204295005f7d00d4864ec9d4d0b Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 30 Apr 2020 11:46:03 -0400 Subject: [PATCH 06/11] +Rescaled sfc_state%sfc_density and sfc_state%Hml Dimensionally rescaled the surface density and mixed layer depth variables in the surface type to [R] and [Z], respectively. All answers are bitwise identical, but the units of two elements in a transparent type have changed. --- config_src/ice_solo_driver/MOM_surface_forcing.F90 | 8 ++++---- config_src/ice_solo_driver/user_surface_forcing.F90 | 9 ++++----- config_src/mct_driver/mom_ocean_model_mct.F90 | 2 +- config_src/nuopc_driver/mom_ocean_model_nuopc.F90 | 2 +- config_src/solo_driver/MESO_surface_forcing.F90 | 9 ++++----- config_src/solo_driver/MOM_surface_forcing.F90 | 8 ++++---- config_src/solo_driver/Neverland_surface_forcing.F90 | 8 ++++---- config_src/solo_driver/user_surface_forcing.F90 | 2 +- src/core/MOM.F90 | 6 +++--- src/core/MOM_checksum_packages.F90 | 3 ++- src/core/MOM_variables.F90 | 4 ++-- src/ice_shelf/MOM_ice_shelf.F90 | 2 +- src/user/BFB_surface_forcing.F90 | 2 +- 13 files changed, 32 insertions(+), 33 deletions(-) diff --git a/config_src/ice_solo_driver/MOM_surface_forcing.F90 b/config_src/ice_solo_driver/MOM_surface_forcing.F90 index 1e59fee863..8e218fb6c4 100644 --- a/config_src/ice_solo_driver/MOM_surface_forcing.F90 +++ b/config_src/ice_solo_driver/MOM_surface_forcing.F90 @@ -112,7 +112,7 @@ module MOM_surface_forcing real, pointer :: T_Restore(:,:) => NULL() !< temperature to damp (restore) the SST to [degC] real, pointer :: S_Restore(:,:) => NULL() !< salinity to damp (restore) the SSS [ppt] - real, pointer :: Dens_Restore(:,:) => NULL() !< density to damp (restore) surface density [kg m-3] + real, pointer :: Dens_Restore(:,:) => NULL() !< density to damp (restore) surface density [R ~> kg m-3] integer :: wind_last_lev_read = -1 !< The last time level read from the wind input files integer :: buoy_last_lev_read = -1 !< The last time level read from buoyancy input files @@ -774,7 +774,7 @@ subroutine buoyancy_forcing_from_files(sfc_state, fluxes, day, dt, G, US, CS) do j=js,je ; do i=is,ie if (G%mask2dT(i,j) > 0) then fluxes%buoy(i,j) = (CS%Dens_Restore(i,j) - sfc_state%sfc_density(i,j)) * & - (CS%G_Earth * CS%Flux_const/(US%R_to_kg_m3*CS%Rho0)) + (CS%G_Earth * CS%Flux_const / CS%Rho0) else fluxes%buoy(i,j) = 0.0 endif @@ -909,8 +909,8 @@ subroutine buoyancy_forcing_linear(sfc_state, fluxes, day, dt, G, US, CS) "RESTOREBUOY to linear not written yet.") !do j=js,je ; do i=is,ie ! if (G%mask2dT(i,j) > 0) then - ! fluxes%buoy(i,j) = US%kg_m3_to_R*(CS%Dens_Restore(i,j) - sfc_state%sfc_density(i,j)) * & - ! (CS%G_Earth * CS%Flux_const / CS%Rho0) + ! fluxes%buoy(i,j) = (CS%Dens_Restore(i,j) - sfc_state%sfc_density(i,j)) * & + ! (CS%G_Earth * CS%Flux_const / CS%Rho0) ! else ! fluxes%buoy(i,j) = 0.0 ! endif diff --git a/config_src/ice_solo_driver/user_surface_forcing.F90 b/config_src/ice_solo_driver/user_surface_forcing.F90 index 10417d4a1e..1b372bf44b 100644 --- a/config_src/ice_solo_driver/user_surface_forcing.F90 +++ b/config_src/ice_solo_driver/user_surface_forcing.F90 @@ -177,8 +177,7 @@ subroutine USER_buoyancy_forcing(sfc_state, fluxes, day, dt, G, US, CS) real :: Temp_restore ! The temperature that is being restored toward [C]. real :: Salin_restore ! The salinity that is being restored toward [ppt] - real :: density_restore ! The potential density that is being restored - ! toward [kg m-3]. + real :: density_restore ! The potential density that is being restored toward [R ~> kg m-3]. real :: rhoXcp ! The mean density times the heat capacity [Q R degC-1 ~> J m-3 degC-1]. real :: buoy_rest_const ! A constant relating density anomalies to the ! restoring buoyancy flux [L2 T-3 R-1 ~> m5 s-3 kg-1]. @@ -271,11 +270,11 @@ subroutine USER_buoyancy_forcing(sfc_state, fluxes, day, dt, G, US, CS) buoy_rest_const = -1.0 * (CS%G_Earth * CS%Flux_const) / CS%Rho0 do j=js,je ; do i=is,ie ! Set density_restore to an expression for the surface potential - ! density [kg m-3] that is being restored toward. - density_restore = 1030.0 + ! density [R ~> kg m-3] that is being restored toward. + density_restore = 1030.0*US%kg_m3_to_R fluxes%buoy(i,j) = G%mask2dT(i,j) * buoy_rest_const * & - US%kg_m3_to_R*(density_restore - sfc_state%sfc_density(i,j)) + (density_restore - sfc_state%sfc_density(i,j)) enddo ; enddo endif endif ! end RESTOREBUOY diff --git a/config_src/mct_driver/mom_ocean_model_mct.F90 b/config_src/mct_driver/mom_ocean_model_mct.F90 index 0ef8a7a9d4..1d1e9fa888 100644 --- a/config_src/mct_driver/mom_ocean_model_mct.F90 +++ b/config_src/mct_driver/mom_ocean_model_mct.F90 @@ -918,7 +918,7 @@ subroutine convert_state_to_ocean_type(sfc_state, Ocean_sfc, G, US, patm, press_ if (allocated(sfc_state%Hml)) then do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd - Ocean_sfc%OBLD(i,j) = sfc_state%Hml(i+i0,j+j0) + Ocean_sfc%OBLD(i,j) = US%Z_to_m*sfc_state%Hml(i+i0,j+j0) enddo ; enddo endif diff --git a/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 b/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 index 1a08869b77..22a4c7eaa2 100644 --- a/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 +++ b/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 @@ -913,7 +913,7 @@ subroutine convert_state_to_ocean_type(sfc_state, Ocean_sfc, G, US, patm, press_ if (allocated(sfc_state%Hml)) then do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd - Ocean_sfc%OBLD(i,j) = sfc_state%Hml(i+i0,j+j0) + Ocean_sfc%OBLD(i,j) = US%Z_to_m*sfc_state%Hml(i+i0,j+j0) enddo ; enddo endif diff --git a/config_src/solo_driver/MESO_surface_forcing.F90 b/config_src/solo_driver/MESO_surface_forcing.F90 index ebe98a3293..cc0939ac17 100644 --- a/config_src/solo_driver/MESO_surface_forcing.F90 +++ b/config_src/solo_driver/MESO_surface_forcing.F90 @@ -79,8 +79,7 @@ subroutine MESO_buoyancy_forcing(sfc_state, fluxes, day, dt, G, US, CS) real :: Temp_restore ! The temperature that is being restored toward [degC]. real :: Salin_restore ! The salinity that is being restored toward [ppt] - real :: density_restore ! The potential density that is being restored - ! toward [kg m-3]. + real :: density_restore ! The potential density that is being restored toward [R ~> kg m-3]. real :: rhoXcp ! The mean density times the heat capacity [Q R degC-1 ~> J m-3 degC-1]. real :: buoy_rest_const ! A constant relating density anomalies to the ! restoring buoyancy flux [L2 T-3 R-1 ~> m5 s-3 kg-1]. @@ -194,11 +193,11 @@ subroutine MESO_buoyancy_forcing(sfc_state, fluxes, day, dt, G, US, CS) buoy_rest_const = -1.0 * (CS%G_Earth * CS%Flux_const) / CS%Rho0 do j=js,je ; do i=is,ie ! Set density_restore to an expression for the surface potential - ! density [kg m-3] that is being restored toward. - density_restore = 1030.0 + ! density [R ~> kg m-3] that is being restored toward. + density_restore = 1030.0 * US%kg_m3_to_R fluxes%buoy(i,j) = G%mask2dT(i,j) * buoy_rest_const * & - US%kg_m3_to_R * (density_restore - sfc_state%sfc_density(i,j)) + (density_restore - sfc_state%sfc_density(i,j)) enddo ; enddo endif endif ! end RESTOREBUOY diff --git a/config_src/solo_driver/MOM_surface_forcing.F90 b/config_src/solo_driver/MOM_surface_forcing.F90 index df403712f7..173d417ff3 100644 --- a/config_src/solo_driver/MOM_surface_forcing.F90 +++ b/config_src/solo_driver/MOM_surface_forcing.F90 @@ -95,7 +95,7 @@ module MOM_surface_forcing real, pointer :: T_Restore(:,:) => NULL() !< temperature to damp (restore) the SST to [degC] real, pointer :: S_Restore(:,:) => NULL() !< salinity to damp (restore) the SSS [ppt] - real, pointer :: Dens_Restore(:,:) => NULL() !< density to damp (restore) surface density [kg m-3] + real, pointer :: Dens_Restore(:,:) => NULL() !< density to damp (restore) surface density [R ~> kg m-3] integer :: buoy_last_lev_read = -1 !< The last time level read from buoyancy input files @@ -1000,7 +1000,7 @@ subroutine buoyancy_forcing_from_files(sfc_state, fluxes, day, dt, G, US, CS) else do j=js,je ; do i=is,ie if (G%mask2dT(i,j) > 0) then - fluxes%buoy(i,j) = US%kg_m3_to_R * (CS%Dens_Restore(i,j) - sfc_state%sfc_density(i,j)) * & + fluxes%buoy(i,j) = (CS%Dens_Restore(i,j) - sfc_state%sfc_density(i,j)) * & (CS%G_Earth * CS%Flux_const / CS%Rho0) else fluxes%buoy(i,j) = 0.0 @@ -1161,7 +1161,7 @@ subroutine buoyancy_forcing_from_data_override(sfc_state, fluxes, day, dt, G, US else do j=js,je ; do i=is,ie if (G%mask2dT(i,j) > 0) then - fluxes%buoy(i,j) = US%kg_m3_to_R * (CS%Dens_Restore(i,j) - sfc_state%sfc_density(i,j)) * & + fluxes%buoy(i,j) = (CS%Dens_Restore(i,j) - sfc_state%sfc_density(i,j)) * & (CS%G_Earth * CS%Flux_const / CS%Rho0) else fluxes%buoy(i,j) = 0.0 @@ -1362,7 +1362,7 @@ subroutine buoyancy_forcing_linear(sfc_state, fluxes, day, dt, G, US, CS) "RESTOREBUOY to linear not written yet.") !do j=js,je ; do i=is,ie ! if (G%mask2dT(i,j) > 0) then - ! fluxes%buoy(i,j) = US%kg_m3_to_R * (CS%Dens_Restore(i,j) - sfc_state%sfc_density(i,j)) * & + ! fluxes%buoy(i,j) = (CS%Dens_Restore(i,j) - sfc_state%sfc_density(i,j)) * & ! (CS%G_Earth * CS%Flux_const / CS%Rho0) ! else ! fluxes%buoy(i,j) = 0.0 diff --git a/config_src/solo_driver/Neverland_surface_forcing.F90 b/config_src/solo_driver/Neverland_surface_forcing.F90 index e6b7152e86..a53eaec27e 100644 --- a/config_src/solo_driver/Neverland_surface_forcing.F90 +++ b/config_src/solo_driver/Neverland_surface_forcing.F90 @@ -148,7 +148,7 @@ subroutine Neverland_buoyancy_forcing(sfc_state, fluxes, day, dt, G, US, CS) ! Local variables real :: buoy_rest_const ! A constant relating density anomalies to the ! restoring buoyancy flux [L2 T-3 R-1 ~> m5 s-3 kg-1]. - real :: density_restore ! De + real :: density_restore ! Density being restored toward [R ~> kg m-3] integer :: i, j, is, ie, js, je integer :: isd, ied, jsd, jed @@ -199,11 +199,11 @@ subroutine Neverland_buoyancy_forcing(sfc_state, fluxes, day, dt, G, US, CS) buoy_rest_const = -1.0 * (CS%G_Earth * CS%Flux_const) / CS%Rho0 do j=js,je ; do i=is,ie ! Set density_restore to an expression for the surface potential - ! density [kg m-3] that is being restored toward. - density_restore = 1030.0 + ! density [R ~> kg m-3] that is being restored toward. + density_restore = 1030.0*US%kg_m3_to_R fluxes%buoy(i,j) = G%mask2dT(i,j) * buoy_rest_const * & - US%kg_m3_to_R*(density_restore - sfc_state%sfc_density(i,j)) + (density_restore - sfc_state%sfc_density(i,j)) enddo ; enddo endif endif ! end RESTOREBUOY diff --git a/config_src/solo_driver/user_surface_forcing.F90 b/config_src/solo_driver/user_surface_forcing.F90 index 97da89e69e..a95046fe20 100644 --- a/config_src/solo_driver/user_surface_forcing.F90 +++ b/config_src/solo_driver/user_surface_forcing.F90 @@ -226,7 +226,7 @@ subroutine USER_buoyancy_forcing(sfc_state, fluxes, day, dt, G, US, CS) density_restore = 1030.0*US%kg_m3_to_R fluxes%buoy(i,j) = G%mask2dT(i,j) * buoy_rest_const * & - (density_restore - US%kg_m3_to_R*sfc_state%sfc_density(i,j)) + (density_restore - sfc_state%sfc_density(i,j)) enddo ; enddo endif endif ! end RESTOREBUOY diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index c47115ed68..77b8bdd7d8 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -3019,7 +3019,7 @@ subroutine extract_surface_state(CS, sfc_state_in) ! copy Hml into sfc_state, so that caps can access it if (associated(CS%Hml)) then do j=js,je ; do i=is,ie - sfc_state%Hml(i,j) = CS%Hml(i,j) + sfc_state%Hml(i,j) = US%m_to_Z*CS%Hml(i,j) enddo ; enddo endif @@ -3064,7 +3064,7 @@ subroutine extract_surface_state(CS, sfc_state_in) sfc_state%SST(i,j) = sfc_state%SST(i,j) + dh * CS%tv%T(i,j,k) sfc_state%SSS(i,j) = sfc_state%SSS(i,j) + dh * CS%tv%S(i,j,k) else - sfc_state%sfc_density(i,j) = sfc_state%sfc_density(i,j) + dh * US%R_to_kg_m3*GV%Rlay(k) + sfc_state%sfc_density(i,j) = sfc_state%sfc_density(i,j) + dh * GV%Rlay(k) endif depth(i) = depth(i) + dh enddo ; enddo @@ -3088,7 +3088,7 @@ subroutine extract_surface_state(CS, sfc_state_in) sfc_state%SSS(i,j) = (sfc_state%SSS(i,j) + missing_depth*CS%tv%S(i,j,1)) * I_depth else sfc_state%sfc_density(i,j) = (sfc_state%sfc_density(i,j) + & - missing_depth*US%R_to_kg_m3*GV%Rlay(1)) * I_depth + missing_depth*GV%Rlay(1)) * I_depth endif else I_depth = 1.0 / depth(i) diff --git a/src/core/MOM_checksum_packages.F90 b/src/core/MOM_checksum_packages.F90 index c63b15a3b4..47369cf474 100644 --- a/src/core/MOM_checksum_packages.F90 +++ b/src/core/MOM_checksum_packages.F90 @@ -160,7 +160,8 @@ subroutine MOM_surface_chksum(mesg, sfc_state, G, haloshift, symmetric) if (allocated(sfc_state%SST)) call hchksum(sfc_state%SST, mesg//" SST", G%HI, haloshift=hs) if (allocated(sfc_state%SSS)) call hchksum(sfc_state%SSS, mesg//" SSS", G%HI, haloshift=hs) if (allocated(sfc_state%sea_lev)) call hchksum(sfc_state%sea_lev, mesg//" sea_lev", G%HI, haloshift=hs) - if (allocated(sfc_state%Hml)) call hchksum(sfc_state%Hml, mesg//" Hml", G%HI, haloshift=hs) + if (allocated(sfc_state%Hml)) call hchksum(sfc_state%Hml, mesg//" Hml", G%HI, haloshift=hs, & + scale=G%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=G%US%L_T_to_m_s) diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index 2b7cafdd33..b74e90ff5e 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -42,8 +42,8 @@ module MOM_variables real, allocatable, dimension(:,:) :: & SST, & !< The sea surface temperature [degC]. SSS, & !< The sea surface salinity [ppt ~> psu or gSalt/kg]. - sfc_density, & !< The mixed layer density [kg m-3]. - Hml, & !< The mixed layer depth [m]. + sfc_density, & !< The mixed layer density [R ~> kg m-3]. + Hml, & !< The mixed layer depth [Z ~> m]. u, & !< The mixed layer zonal velocity [L T-1 ~> m s-1]. v, & !< The mixed layer meridional velocity [L T-1 ~> m s-1]. sea_lev, & !< The sea level [m]. If a reduced surface gravity is diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 3e8b1ed83b..d104dfe82a 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -397,7 +397,7 @@ subroutine shelf_calc_flux(sfc_state, fluxes, Time, time_step, CS, forces) ! reported ocean mixed layer thickness and the neutral Ekman depth. absf = 0.25*((abs(G%CoriolisBu(I,J)) + abs(G%CoriolisBu(I-1,J-1))) + & (abs(G%CoriolisBu(I,J-1)) + abs(G%CoriolisBu(I-1,J)))) - if (absf*US%m_to_Z*sfc_state%Hml(i,j) <= VK*ustar_h) then ; hBL_neut = US%m_to_Z*sfc_state%Hml(i,j) + if (absf*sfc_state%Hml(i,j) <= VK*ustar_h) then ; hBL_neut = sfc_state%Hml(i,j) else ; hBL_neut = (VK*ustar_h) / absf ; endif hBL_neut_h_molec = ZETA_N * ((hBL_neut * ustar_h) / (5.0 * CS%kv_molec)) diff --git a/src/user/BFB_surface_forcing.F90 b/src/user/BFB_surface_forcing.F90 index a6aae3d3f7..88e7ae45d5 100644 --- a/src/user/BFB_surface_forcing.F90 +++ b/src/user/BFB_surface_forcing.F90 @@ -164,7 +164,7 @@ subroutine BFB_buoyancy_forcing(sfc_state, fluxes, day, dt, G, US, CS) density_restore = Temp_restore*CS%drho_dt + CS%Rho0 fluxes%buoy(i,j) = G%mask2dT(i,j) * buoy_rest_const * & - (density_restore - US%kg_m3_to_R*sfc_state%sfc_density(i,j)) + (density_restore - sfc_state%sfc_density(i,j)) enddo ; enddo endif endif ! end RESTOREBUOY From 5f7ce525b3b9b1e984cbaae4db4f37831d8696e3 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 30 Apr 2020 18:53:05 -0400 Subject: [PATCH 07/11] (*)Rescale frazil and p_surf_prev after restarts Added code to rescale frazil and p_surf_prev after restarts if the dimensional rescaling factors have changed. All answers are bitwise identical in test cases but there are some peculiar instances when incorrect behavior is avoided. --- src/core/MOM.F90 | 32 ++++++++++++++++++++++++++++++-- 1 file changed, 30 insertions(+), 2 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 77b8bdd7d8..34266b51ab 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1670,6 +1670,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & integer :: dynamics_stencil ! The computational stencil for the calculations ! in the dynamic core. real :: conv2watt, conv2salt + real :: RL2_T2_rescale, Z_rescale, QRZ_rescale ! Unit conversion factors character(len=48) :: flux_units, S_flux_units type(vardesc) :: vd_T, vd_S ! Structures describing temperature and salinity variables. @@ -2665,14 +2666,39 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & call register_obsolete_diagnostics(param_file, CS%diag) if (use_frazil) then - if (.not.query_initialized(CS%tv%frazil,"frazil",restart_CSp)) & + if (.not.query_initialized(CS%tv%frazil,"frazil",restart_CSp)) then + ! Test whether the dimensional rescaling has changed for heat content. + if ((US%kg_m3_to_R_restart*US%m_to_Z_restart*US%J_kg_to_Q_restart /= 0.0) .and. & + ((US%J_kg_to_Q*US%kg_m3_to_R*US%m_to_Z) /= & + (US%J_kg_to_Q_restart*US%kg_m3_to_R_restart*US%m_to_Z_restart)) ) then + QRZ_rescale = (US%J_kg_to_Q*US%kg_m3_to_R*US%m_to_Z) / & + (US%J_kg_to_Q_restart*US%kg_m3_to_R_restart*US%m_to_Z_restart) + do j=js,je ; do i=is,ie + CS%tv%frazil(i,j) = QRZ_rescale * CS%tv%frazil(i,j) + enddo ; enddo + endif + else CS%tv%frazil(:,:) = 0.0 + endif endif if (CS%interp_p_surf) then CS%p_surf_prev_set = query_initialized(CS%p_surf_prev,"p_surf_prev",restart_CSp) - if (CS%p_surf_prev_set) call pass_var(CS%p_surf_prev, G%domain) + if (CS%p_surf_prev_set) then + ! Test whether the dimensional rescaling has changed for pressure. + if ((US%kg_m3_to_R_restart*US%s_to_T_restart*US%m_to_L_restart /= 0.0) .and. & + ((US%kg_m3_to_R*(US%m_to_L*US%s_to_T_restart)**2) /= & + (US%kg_m3_to_R_restart*(US%m_to_L_restart*US%s_to_T)**2)) ) then + RL2_T2_rescale = (US%kg_m3_to_R*(US%m_to_L*US%s_to_T_restart)**2) / & + (US%kg_m3_to_R_restart*(US%m_to_L_restart*US%s_to_T)**2) + do j=js,je ; do i=is,ie + CS%p_surf_prev(i,j) = RL2_T2_rescale * CS%p_surf_prev(i,j) + enddo ; enddo + endif + + call pass_var(CS%p_surf_prev, G%domain) + endif endif if (.not.query_initialized(CS%ave_ssh_ibc,"ave_ssh",restart_CSp)) then @@ -2891,6 +2917,8 @@ subroutine set_restart_fields(GV, US, param_file, CS, restart_CSp) "Time unit conversion factor", "T second-1") call register_restart_field(US%kg_m3_to_R_restart, "kg_m3_to_R", .false., restart_CSp, & "Density unit conversion factor", "R m3 kg-1") + call register_restart_field(US%J_kg_to_Q_restart, "J_kg_to_Q", .false., restart_CSp, & + "Heat content unit conversion factor.", units="Q kg J-1") end subroutine set_restart_fields From a90a65ec8898a99bbef83e659b336cbf2c89fa8f Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 30 Apr 2020 18:56:45 -0400 Subject: [PATCH 08/11] +(*)Added US and m_to_BLD_units args to KPP_get_BLD Added a unit scaling type argument and an optional argument, m_to_BLD_units, to specify the units with which KPP_get_BLD returns the mixed layer depths. Also change the default units of the turbulent boundary layer depths returned from KPP_get_BLD and energetic_PBL_get_MLD to [Z]. In addition, corrected the documented units for various variables related to the recently added MOM_lateral_boundary_diffusion module. All answers are bitwise identical by default and in the MOM6-examples test cases, but some improper dimensional rescaling has been corrected for non-Boussinesq configurations using MOM_lateral_boundary_diffusion. --- .../vertical/MOM_CVMix_KPP.F90 | 16 +++- .../vertical/MOM_diabatic_driver.F90 | 18 ++-- .../vertical/MOM_energetic_PBL.F90 | 10 +-- src/tracer/MOM_lateral_boundary_diffusion.F90 | 88 ++++++++++--------- src/tracer/MOM_neutral_diffusion.F90 | 7 +- src/tracer/MOM_tracer_hor_diff.F90 | 8 +- 6 files changed, 80 insertions(+), 67 deletions(-) diff --git a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 index 6ff6046350..3b7420aa54 100644 --- a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 @@ -1365,17 +1365,25 @@ end subroutine KPP_smooth_BLD -!> Copies KPP surface boundary layer depth into BLD -subroutine KPP_get_BLD(CS, BLD, G) +!> Copies KPP surface boundary layer depth into BLD, in units of [Z ~> m] unless other units are specified. +subroutine KPP_get_BLD(CS, BLD, G, US, m_to_BLD_units) type(KPP_CS), pointer :: CS !< Control structure for !! this module type(ocean_grid_type), intent(in) :: G !< Grid structure - real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: BLD!< bnd. layer depth [m] + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: BLD !< Boundary layer depth [Z ~> m] or other units + real, optional, intent(in) :: m_to_BLD_units !< A conversion factor from meters + !! to the desired units for BLD ! Local variables + real :: scale ! A dimensional rescaling factor integer :: i,j + + scale = US%m_to_Z ; if (present(m_to_BLD_units)) scale = m_to_BLD_units + do j = G%jsc, G%jec ; do i = G%isc, G%iec - BLD(i,j) = CS%OBLdepth(i,j) + BLD(i,j) = scale * CS%OBLdepth(i,j) enddo ; enddo + end subroutine KPP_get_BLD !> Apply KPP non-local transport of surface fluxes for temperature. diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 69c800d218..3c017c0f6c 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -663,9 +663,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim Kd_salt, visc%Kv_shear, CS%KPP_NLTheat, CS%KPP_NLTscalar, Waves=Waves) if (associated(Hml)) then - !$OMP parallel default(shared) - call KPP_get_BLD(CS%KPP_CSp, Hml(:,:), G) - !$OMP end parallel + call KPP_get_BLD(CS%KPP_CSp, Hml(:,:), G, US, m_to_BLD_units=1.0) call pass_var(Hml, G%domain, halo=1) ! If visc%MLD exists, copy KPP's BLD into it if (associated(visc%MLD)) visc%MLD(:,:) = Hml(:,:) @@ -842,12 +840,12 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim CS%energetic_PBL_CSp, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves) if (associated(Hml)) then - call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, Hml(:,:), G, US) + call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, Hml(:,:), G, US, m_to_MLD_units=1.0) call pass_var(Hml, G%domain, halo=1) ! If visc%MLD exists, copy ePBL's MLD into it if (associated(visc%MLD)) visc%MLD(:,:) = Hml(:,:) elseif (associated(visc%MLD)) then - call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, visc%MLD, G, US) + call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, visc%MLD, G, US, m_to_MLD_units=1.0) call pass_var(visc%MLD, G%domain, halo=1) endif @@ -1448,9 +1446,7 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, Kd_salt, visc%Kv_shear, CS%KPP_NLTheat, CS%KPP_NLTscalar, Waves=Waves) if (associated(Hml)) then - !$OMP parallel default(shared) - call KPP_get_BLD(CS%KPP_CSp, Hml(:,:), G) - !$OMP end parallel + call KPP_get_BLD(CS%KPP_CSp, Hml(:,:), G, US, m_to_BLD_units=1.0) call pass_var(Hml, G%domain, halo=1) ! If visc%MLD exists, copy KPP's BLD into it if (associated(visc%MLD)) visc%MLD(:,:) = Hml(:,:) @@ -1573,12 +1569,12 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, CS%energetic_PBL_CSp, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves) if (associated(Hml)) then - call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, Hml(:,:), G, US) + call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, Hml(:,:), G, US, m_to_MLD_units=1.0) call pass_var(Hml, G%domain, halo=1) ! If visc%MLD exists, copy ePBL's MLD into it if (associated(visc%MLD)) visc%MLD(:,:) = Hml(:,:) elseif (associated(visc%MLD)) then - call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, visc%MLD, G, US) + call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, visc%MLD, G, US, m_to_MLD_units=1.0) call pass_var(visc%MLD, G%domain, halo=1) endif @@ -2183,7 +2179,7 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e Kd_salt, visc%Kv_shear, CS%KPP_NLTheat, CS%KPP_NLTscalar, Waves=Waves) if (associated(Hml)) then - call KPP_get_BLD(CS%KPP_CSp, Hml(:,:), G) + call KPP_get_BLD(CS%KPP_CSp, Hml(:,:), G, US, m_to_BLD_units=1.0) call pass_var(Hml, G%domain, halo=1) ! If visc%MLD exists, copy KPP's BLD into it if (associated(visc%MLD)) visc%MLD(:,:) = Hml(:,:) diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index 48b265a0e2..a9e68736e7 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -1908,19 +1908,19 @@ subroutine Mstar_Langmuir(CS, US, Abs_Coriolis, Buoyancy_Flux, UStar, BLD, Langm end subroutine Mstar_Langmuir -!> Copies the ePBL active mixed layer depth into MLD +!> Copies the ePBL active mixed layer depth into MLD, in units of [Z ~> m] unless other units are specified. subroutine energetic_PBL_get_MLD(CS, MLD, G, US, m_to_MLD_units) type(energetic_PBL_CS), pointer :: CS !< Control structure for ePBL type(ocean_grid_type), intent(in) :: G !< Grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - real, dimension(SZI_(G),SZJ_(G)), intent(out) :: MLD !< Depth of ePBL active mixing layer [m or other units] - real, optional, intent(in) :: m_to_MLD_units !< A conversion factor to the - !! desired units for MLD + real, dimension(SZI_(G),SZJ_(G)), intent(out) :: MLD !< Depth of ePBL active mixing layer [Z ~> m] or other units + real, optional, intent(in) :: m_to_MLD_units !< A conversion factor from meters + !! to the desired units for MLD ! Local variables real :: scale ! A dimensional rescaling factor integer :: i,j - scale = US%Z_to_m ; if (present(m_to_MLD_units)) scale = scale * m_to_MLD_units + scale = 1.0 ; if (present(m_to_MLD_units)) scale = US%Z_to_m * m_to_MLD_units do j=G%jsc,G%jec ; do i=G%isc,G%iec MLD(i,j) = scale*CS%ML_Depth(i,j) diff --git a/src/tracer/MOM_lateral_boundary_diffusion.F90 b/src/tracer/MOM_lateral_boundary_diffusion.F90 index d06127d0d5..8b9be533d5 100644 --- a/src/tracer/MOM_lateral_boundary_diffusion.F90 +++ b/src/tracer/MOM_lateral_boundary_diffusion.F90 @@ -134,15 +134,15 @@ subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & intent(in) :: h !< Layer thickness [H ~> m or kg m-2] - real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: Coef_x !< dt * Kh * dy / dx at u-points [m2] - real, dimension(SZI_(G),SZJB_(G)), intent(in) :: Coef_y !< dt * Kh * dx / dy at v-points [m2] + real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: Coef_x !< dt * Kh * dy / dx at u-points [L2 ~> m2] + real, dimension(SZI_(G),SZJB_(G)), intent(in) :: Coef_y !< dt * Kh * dx / dy at v-points [L2 ~> m2] real, intent(in) :: dt !< Tracer time step * I_numitts !! (I_numitts in tracer_hordiff) type(tracer_registry_type), pointer :: Reg !< Tracer registry type(lateral_boundary_diffusion_CS), intent(in) :: CS !< Control structure for this module ! Local variables - real, dimension(SZI_(G),SZJ_(G)) :: hbl !< bnd. layer depth [m] + real, dimension(SZI_(G),SZJ_(G)) :: hbl !< bnd. layer depth [H ~> m or kg m-2] real, dimension(SZI_(G),SZJ_(G),SZK_(G),CS%deg+1) :: ppoly0_coefs !< Coefficients of polynomial real, dimension(SZI_(G),SZJ_(G),SZK_(G),2) :: ppoly0_E !< Edge values from reconstructions real, dimension(SZK_(G),CS%deg+1) :: ppoly_S !< Slopes from reconstruction (placeholder) @@ -161,8 +161,9 @@ subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) Idt = 1./dt hbl(:,:) = 0. - if (ASSOCIATED(CS%KPP_CSp)) call KPP_get_BLD(CS%KPP_CSp, hbl, G) - if (ASSOCIATED(CS%energetic_PBL_CSp)) call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, hbl, G, US) + if (ASSOCIATED(CS%KPP_CSp)) call KPP_get_BLD(CS%KPP_CSp, hbl, G, US, m_to_BLD_units=GV%m_to_H) + if (ASSOCIATED(CS%energetic_PBL_CSp)) & + call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, hbl, G, US, m_to_MLD_units=GV%m_to_H) call pass_var(hbl,G%Domain) do m = 1,Reg%ntr @@ -284,7 +285,7 @@ subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) ! post tendency of tracer concentration; this step must be ! done after posting tracer content tendency, since we alter - ! the tendency array. + ! the tendency array and its units. if (tracer%id_lbdxy_conc > 0) then do k = 1, GV%ke ; do j = G%jsc,G%jec ; do i = G%isc,G%iec tendency(i,j,k) = tendency(i,j,k) / ( h(i,j,k) + GV%H_subroundoff ) @@ -302,8 +303,8 @@ real function bulk_average(boundary, nk, deg, h, hBLT, phi, ppoly0_E, ppoly0_coe integer :: boundary !< SURFACE or BOTTOM [nondim] integer :: nk !< Number of layers [nondim] integer :: deg !< Degree of polynomial [nondim] - real, dimension(nk) :: h !< Layer thicknesses [m] - real :: hBLT !< Depth of the boundary layer [m] + real, dimension(nk) :: h !< Layer thicknesses [H ~> m or kg m-2] + real :: hBLT !< Depth of the boundary layer [H ~> m or kg m-2] real, dimension(nk) :: phi !< Scalar quantity real, dimension(nk,2) :: ppoly0_E(:,:) !< Edge value of polynomial real, dimension(nk,deg+1) :: ppoly0_coefs(:,:) !< Coefficients of polynomial @@ -318,7 +319,7 @@ real function bulk_average(boundary, nk, deg, h, hBLT, phi, ppoly0_E, ppoly0_coe !! (0 if none, 1. if all). For the bottom boundary layer, this is always 1. !! because integration starts at the bottom [nondim] ! Local variables - real :: htot !< Running sum of the thicknesses (top to bottom) + real :: htot !< Running sum of the thicknesses (top to bottom) [H ~> m or kg m-2] integer :: k !< k indice @@ -364,8 +365,8 @@ end function harmonic_mean subroutine boundary_k_range(boundary, nk, h, hbl, k_top, zeta_top, k_bot, zeta_bot) integer, intent(in ) :: boundary !< SURFACE or BOTTOM [nondim] integer, intent(in ) :: nk !< Number of layers [nondim] - real, dimension(nk), intent(in ) :: h !< Layer thicknesses of the column [m] - real, intent(in ) :: hbl !< Thickness of the boundary layer [m] + real, dimension(nk), intent(in ) :: h !< Layer thicknesses of the column [H ~> m or kg m-2] + real, intent(in ) :: hbl !< Thickness of the boundary layer [H ~> m or kg m-2] !! If surface, with respect to zbl_ref = 0. !! If bottom, with respect to zbl_ref = SUM(h) integer, intent( out) :: k_top !< Index of the first layer within the boundary @@ -375,7 +376,7 @@ subroutine boundary_k_range(boundary, nk, h, hbl, k_top, zeta_top, k_bot, zeta_b real, intent( out) :: zeta_bot !< Distance of the lower layer to the boundary layer depth !! (0 at top, 1 at bottom) [nondim] ! Local variables - real :: htot + real :: htot ! Summed thickness [H ~> m or kg m-2] integer :: k ! Surface boundary layer if ( boundary == SURFACE ) then @@ -434,14 +435,14 @@ subroutine fluxes_layer_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L integer, intent(in ) :: boundary !< Which boundary layer SURFACE or BOTTOM [nondim] integer, intent(in ) :: nk !< Number of layers [nondim] integer, intent(in ) :: deg !< order of the polynomial reconstruction [nondim] - real, dimension(nk), intent(in ) :: h_L !< Layer thickness (left) [m] - real, dimension(nk), intent(in ) :: h_R !< Layer thickness (right) [m] + real, dimension(nk), intent(in ) :: h_L !< Layer thickness (left) [H ~> m or kg m-2] + real, dimension(nk), intent(in ) :: h_R !< Layer thickness (right) [H ~> m or kg m-2] real, intent(in ) :: hbl_L !< Thickness of the boundary boundary - !! layer (left) [m] + !! layer (left) [H ~> m or kg m-2] real, intent(in ) :: hbl_R !< Thickness of the boundary boundary - !! layer (right) [m] - real, intent(in ) :: area_L !< Area of the horizontal grid (left) [m^2] - real, intent(in ) :: area_R !< Area of the horizontal grid (right) [m^2] + !! layer (right) [H ~> m or kg m-2] + real, intent(in ) :: area_L !< Area of the horizontal grid (left) [L2 ~> m2] + real, intent(in ) :: area_R !< Area of the horizontal grid (right) [L2 ~> m2] real, dimension(nk), intent(in ) :: phi_L !< Tracer values (left) [conc] real, dimension(nk), intent(in ) :: phi_R !< Tracer values (right) [conc] real, dimension(nk,deg+1), intent(in ) :: ppoly0_coefs_L !< Tracer reconstruction (left) [conc] @@ -449,19 +450,22 @@ subroutine fluxes_layer_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L real, dimension(nk,2), intent(in ) :: ppoly0_E_L !< Polynomial edge values (left) [ nondim ] real, dimension(nk,2), intent(in ) :: ppoly0_E_R !< Polynomial edge values (right) [ nondim ] integer, intent(in ) :: method !< Method of polynomial integration [ nondim ] - real, intent(in ) :: khtr_u !< Horizontal diffusivities times delta t at U-point [m^2] - real, dimension(nk), intent( out) :: F_layer !< Layerwise diffusive flux at U- or V-point [m^3 conc] + real, intent(in ) :: khtr_u !< Horizontal diffusivities times delta t + !! at a velocity point [L2 ~> m2] + real, dimension(nk), intent( out) :: F_layer !< Layerwise diffusive flux at U- or V-point + !! [H L2 conc ~> m3 conc] ! Local variables - real, dimension(nk) :: h_means !< Calculate the layer-wise harmonic means [m] - real :: khtr_avg !< Thickness-weighted diffusivity at the u-point [m^2 s^-1] + real, dimension(nk) :: h_means !< Calculate the layer-wise harmonic means [H ~> m or kg m-2] + real :: khtr_avg !< Thickness-weighted diffusivity at the u-point [m^2 s^-1] !! This is just to remind developers that khtr_avg should be !! computed once khtr is 3D. - real :: heff !< Harmonic mean of layer thicknesses [m] - real :: inv_heff !< Inverse of the harmonic mean of layer thicknesses [m^[-1] + real :: heff !< Harmonic mean of layer thicknesses [H ~> m or kg m-2] + real :: inv_heff !< Inverse of the harmonic mean of layer thicknesses + !! [H-1 ~> m-1 or m2 kg-1] real :: phi_L_avg, phi_R_avg !< Bulk, thickness-weighted tracer averages (left and right column) - !! [conc m^-3 ] - real :: htot !< Total column thickness [m] + !! [conc m^-3 ] + real :: htot !< Total column thickness [H ~> m or kg m-2] integer :: k, k_bot_min, k_top_max !< k-indices, min and max for top and bottom, respectively integer :: k_top_L, k_bot_L !< k-indices left integer :: k_top_R, k_bot_R !< k-indices right @@ -547,14 +551,14 @@ subroutine fluxes_bulk_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, integer, intent(in ) :: boundary !< Which boundary layer SURFACE or BOTTOM [nondim] integer, intent(in ) :: nk !< Number of layers [nondim] integer, intent(in ) :: deg !< order of the polynomial reconstruction [nondim] - real, dimension(nk), intent(in ) :: h_L !< Layer thickness (left) [m] - real, dimension(nk), intent(in ) :: h_R !< Layer thickness (right) [m] + real, dimension(nk), intent(in ) :: h_L !< Layer thickness (left) [H ~> m or kg m-2] + real, dimension(nk), intent(in ) :: h_R !< Layer thickness (right) [H ~> m or kg m-2] real, intent(in ) :: hbl_L !< Thickness of the boundary boundary - !! layer (left) [m] + !! layer (left) [H ~> m or kg m-2] real, intent(in ) :: hbl_R !< Thickness of the boundary boundary - !! layer (left) [m] - real, intent(in ) :: area_L !< Area of the horizontal grid (left) [m^2] - real, intent(in ) :: area_R !< Area of the horizontal grid (right) [m^2] + !! layer (left) [H ~> m or kg m-2] + real, intent(in ) :: area_L !< Area of the horizontal grid (left) [L2 ~> m2] + real, intent(in ) :: area_R !< Area of the horizontal grid (right) [L2 ~> m2] real, dimension(nk), intent(in ) :: phi_L !< Tracer values (left) [conc] real, dimension(nk), intent(in ) :: phi_R !< Tracer values (right) [conc] real, dimension(nk,deg+1), intent(in ) :: ppoly0_coefs_L !< Tracer reconstruction (left) [conc] @@ -562,21 +566,25 @@ subroutine fluxes_bulk_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, real, dimension(nk,2), intent(in ) :: ppoly0_E_L !< Polynomial edge values (left) [nondim] real, dimension(nk,2), intent(in ) :: ppoly0_E_R !< Polynomial edge values (right) [nondim] integer, intent(in ) :: method !< Method of polynomial integration [nondim] - real, intent(in ) :: khtr_u !< Horizontal diffusivities times delta t at U-point [m^2] - real, intent( out) :: F_bulk !< The bulk mixed layer lateral flux [m^3 conc] - real, dimension(nk), intent( out) :: F_layer !< Layerwise diffusive flux at U-point [m^3 conc] + real, intent(in ) :: khtr_u !< Horizontal diffusivities times delta t + !! at a velocity point [L2 ~> m2] + real, intent( out) :: F_bulk !< The bulk mixed layer lateral flux + !! [H L2 conc ~> m3 conc] + real, dimension(nk), intent( out) :: F_layer !< Layerwise diffusive flux at U- or V-point + !! [H L2 conc ~> m3 conc] logical, optional, intent(in ) :: F_limit !< If True, apply a limiter ! Local variables - real, dimension(nk) :: h_means !< Calculate the layer-wise harmonic means [m] + real, dimension(nk) :: h_means !< Calculate the layer-wise harmonic means [H ~> m or kg m-2] real :: khtr_avg !< Thickness-weighted diffusivity at the u-point [m^2 s^-1] !! This is just to remind developers that khtr_avg should be !! computed once khtr is 3D. - real :: heff !< Harmonic mean of layer thicknesses [m] - real :: inv_heff !< Inverse of the harmonic mean of layer thicknesses [m^[-1] + real :: heff !< Harmonic mean of layer thicknesses [H ~> m or kg m-2] + real :: inv_heff !< Inverse of the harmonic mean of layer thicknesses + !! [H-1 ~> m-1 or m2 kg-1] real :: phi_L_avg, phi_R_avg !< Bulk, thickness-weighted tracer averages (left and right column) !! [conc m^-3 ] - real :: htot ! Total column thickness [m] + real :: htot !< Total column thickness [H ~> m or kg m-2] integer :: k, k_min, k_max !< k-indices, min and max for top and bottom, respectively integer :: k_top_L, k_bot_L !< k-indices left integer :: k_top_R, k_bot_R !< k-indices right @@ -728,7 +736,7 @@ logical function near_boundary_unit_tests( verbose ) real, dimension(nk,2) :: ppoly0_E_L, ppoly0_E_R! Polynomial edge values (left and right) [concentration] real, dimension(nk) :: h_L, h_R ! Layer thickness (left and right) [m] real :: khtr_u ! Horizontal diffusivities at U-point [m^2 s^-1] - real :: hbl_L, hbl_R ! Depth of the boundary layer (left and right) [m] + real :: hbl_L, hbl_R ! Depth of the boundary layer (left and right) [m] real :: F_bulk ! Total diffusive flux across the U point [nondim s^-1] real, dimension(nk) :: F_layer ! Diffusive flux within each layer at U-point [nondim s^-1] real :: h_u, hblt_u ! Thickness at the u-point [m] diff --git a/src/tracer/MOM_neutral_diffusion.F90 b/src/tracer/MOM_neutral_diffusion.F90 index 48678e1107..30cdec3b37 100644 --- a/src/tracer/MOM_neutral_diffusion.F90 +++ b/src/tracer/MOM_neutral_diffusion.F90 @@ -297,7 +297,7 @@ subroutine neutral_diffusion_calc_coeffs(G, GV, US, h, T, S, CS, p_surf) ! Variables used for reconstructions real, dimension(SZK_(G),2) :: ppoly_r_S ! Reconstruction slopes real, dimension(SZI_(G), SZJ_(G)) :: hEff_sum ! Summed effective face thicknesses [H ~> m or kg m-2] - real, dimension(SZI_(G),SZJ_(G)) :: hbl ! Boundary layer depth [m] + real, dimension(SZI_(G),SZJ_(G)) :: hbl ! Boundary layer depth [H ~> m or kg m-2] integer :: iMethod real, dimension(SZI_(G)) :: ref_pres ! Reference pressure used to calculate alpha/beta [R L2 T-2 ~> Pa] real, dimension(SZI_(G)) :: rho_tmp ! Routine to calculate drho_dp, returns density which is not used @@ -317,8 +317,9 @@ subroutine neutral_diffusion_calc_coeffs(G, GV, US, h, T, S, CS, p_surf) ! Check if hbl needs to be extracted if (CS%interior_only) then hbl(:,:) = 0. - if (ASSOCIATED(CS%KPP_CSp)) call KPP_get_BLD(CS%KPP_CSp, hbl, G) - if (ASSOCIATED(CS%energetic_PBL_CSp)) call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, hbl, G, US) + if (ASSOCIATED(CS%KPP_CSp)) call KPP_get_BLD(CS%KPP_CSp, hbl, G, US, m_to_BLD_units=GV%m_to_H) + if (ASSOCIATED(CS%energetic_PBL_CSp)) & + call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, hbl, G, US, m_to_MLD_units=GV%m_to_H) call pass_var(hbl, G%Domain) ! get k-indices and zeta do j=G%jsc-1, G%jec+1 ; do i=G%isc-1,G%iec+1 diff --git a/src/tracer/MOM_tracer_hor_diff.F90 b/src/tracer/MOM_tracer_hor_diff.F90 index cdbaaf28b9..02255d9424 100644 --- a/src/tracer/MOM_tracer_hor_diff.F90 +++ b/src/tracer/MOM_tracer_hor_diff.F90 @@ -141,14 +141,14 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, G, GV, US, CS, Reg, tv, do_online real, dimension(SZIB_(G),SZJ_(G)) :: & khdt_x, & ! The value of Khtr*dt times the open face width divided by ! the distance between adjacent tracer points [L2 ~> m2]. - Coef_x, & ! The coefficients relating zonal tracer differences - ! to time-integrated fluxes [H L2 ~> m3 or kg]. + Coef_x, & ! The coefficients relating zonal tracer differences to time-integrated + ! fluxes, in [L2 ~> m2] for some schemes and [H L2 ~> m3 or kg] for others. Kh_u ! Tracer mixing coefficient at u-points [L2 T-1 ~> m2 s-1]. real, dimension(SZI_(G),SZJB_(G)) :: & khdt_y, & ! The value of Khtr*dt times the open face width divided by ! the distance between adjacent tracer points [L2 ~> m2]. - Coef_y, & ! The coefficients relating meridional tracer differences - ! to time-integrated fluxes [H L2 ~> m3 or kg]. + Coef_y, & ! The coefficients relating meridional tracer differences to time-integrated + ! fluxes, in [L2 ~> m2] for some schemes and [H L2 ~> m3 or kg] for others. Kh_v ! Tracer mixing coefficient at u-points [L2 T-1 ~> m2 s-1]. real :: khdt_max ! The local limiting value of khdt_x or khdt_y [L2 ~> m2]. From d820561d7425a21fe50c2ab6fd11edc09b46b771 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 30 Apr 2020 20:14:55 -0400 Subject: [PATCH 09/11] +Rescaled visc%MLD and Hml in diabatic to [Z] Rescaled the units to [Z] of the visc%MLD element of the vertvisc_type and the Hml argument to diabatic, which in turn is an argument to bulkmixedlayer, calculate_CVMix_conv, call_tracer_column_fns, MOM_generic_tracer_column_physics and mixedlayer_restrat and stored in the MOM_control_struct as CS%Hml. In addition when visc%MLD is read from a restart file, it is rescaled if the scaling has changed between run segments. All answers are bitwise identical. --- src/core/MOM.F90 | 18 ++++++++++--- src/core/MOM_variables.F90 | 2 +- .../lateral/MOM_mixed_layer_restrat.F90 | 10 +++---- .../vertical/MOM_CVMix_conv.F90 | 18 ++++++------- .../vertical/MOM_bulk_mixed_layer.F90 | 4 +-- .../vertical/MOM_diabatic_driver.F90 | 22 ++++++++-------- .../vertical/MOM_set_viscosity.F90 | 12 ++++++++- src/tracer/MOM_generic_tracer.F90 | 26 ++++++++++--------- src/tracer/MOM_offline_main.F90 | 2 +- src/tracer/MOM_tracer_flow_control.F90 | 2 +- 10 files changed, 69 insertions(+), 47 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 34266b51ab..e3bd314273 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -184,7 +184,7 @@ module MOM !< free surface height or column mass time averaged over the last !! baroclinic dynamics time step [H ~> m or kg m-2] real, dimension(:,:), pointer :: & - Hml => NULL() !< active mixed layer depth [m] + Hml => NULL() !< active mixed layer depth [Z ~> m] real :: time_in_cycle !< The running time of the current time-stepping cycle !! in calls that step the dynamics, and also the length of !! the time integral of ssh_rint [T ~> s]. @@ -2701,6 +2701,18 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & endif endif + if (use_ice_shelf .and. associated(CS%Hml)) then + if (query_initialized(CS%Hml, "hML", restart_CSp)) then + ! Test whether the dimensional rescaling has changed for depths. + if ((US%m_to_Z_restart /= 0.0) .and. (US%m_to_Z /= US%m_to_Z_restart) ) then + Z_rescale = US%m_to_Z / US%m_to_Z_restart + do j=js,je ; do i=is,ie + CS%Hml(i,j) = Z_rescale * CS%Hml(i,j) + enddo ; enddo + endif + endif + endif + if (.not.query_initialized(CS%ave_ssh_ibc,"ave_ssh",restart_CSp)) then if (CS%split) then call find_eta(CS%h, CS%tv, G, GV, US, CS%ave_ssh_ibc, eta, eta_to_m=1.0) @@ -3047,7 +3059,7 @@ subroutine extract_surface_state(CS, sfc_state_in) ! copy Hml into sfc_state, so that caps can access it if (associated(CS%Hml)) then do j=js,je ; do i=is,ie - sfc_state%Hml(i,j) = US%m_to_Z*CS%Hml(i,j) + sfc_state%Hml(i,j) = CS%Hml(i,j) enddo ; enddo endif @@ -3205,7 +3217,7 @@ subroutine extract_surface_state(CS, sfc_state_in) enddo do k=1,nz ; do i=is,ie - depth_ml = min(CS%HFrz, CS%visc%MLD(i,j)) + depth_ml = min(CS%HFrz, US%Z_to_m*CS%visc%MLD(i,j)) if (depth(i) + h(i,j,k)*GV%H_to_m < depth_ml) then dh = h(i,j,k)*GV%H_to_m elseif (depth(i) < depth_ml) then diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index b74e90ff5e..a8a367a7fe 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -229,7 +229,7 @@ module MOM_variables real, pointer, dimension(:,:) :: nkml_visc_v => NULL() !< The number of layers in the viscous surface mixed layer at v-points [nondim]. real, pointer, dimension(:,:) :: & - MLD => NULL() !< Instantaneous active mixing layer depth in unscaled MKS units [m]. + MLD => NULL() !< Instantaneous active mixing layer depth [Z ~> m]. real, pointer, dimension(:,:,:) :: & Ray_u => NULL(), & !< The Rayleigh drag velocity to be applied to each layer at u-points [Z T-1 ~> m s-1]. Ray_v => NULL() !< The Rayleigh drag velocity to be applied to each layer at v-points [Z T-1 ~> m s-1]. diff --git a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 index c1b608b16b..3a3a25429c 100644 --- a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 +++ b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 @@ -101,7 +101,7 @@ subroutine mixedlayer_restrat(h, uhtr, vhtr, tv, forces, dt, MLD, VarMix, G, GV, type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces real, intent(in) :: dt !< Time increment [T ~> s] real, dimension(:,:), pointer :: MLD !< Mixed layer depth provided by the - !! PBL scheme [H ~> m or kg m-2] + !! PBL scheme [Z ~> m] type(VarMix_CS), pointer :: VarMix !< Container for derived fields type(mixedlayer_restrat_CS), pointer :: CS !< Module control structure @@ -131,7 +131,7 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces real, intent(in) :: dt !< Time increment [T ~> s] real, dimension(:,:), pointer :: MLD_in !< Mixed layer depth provided by the - !! PBL scheme [m] (not H) + !! PBL scheme [Z ~> m] (not H) type(VarMix_CS), pointer :: VarMix !< Container for derived fields type(mixedlayer_restrat_CS), pointer :: CS !< Module control structure ! Local variables @@ -240,7 +240,7 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var if (.not. associated(MLD_in)) call MOM_error(FATAL, "MOM_mixedlayer_restrat: "// & "Argument MLD_in was not associated!") do j = js-1, je+1 ; do i = is-1, ie+1 - MLD_fast(i,j) = (CS%MLE_MLD_stretch * GV%m_to_H) * MLD_in(i,j) + MLD_fast(i,j) = (CS%MLE_MLD_stretch * GV%Z_to_H) * MLD_in(i,j) enddo ; enddo else call MOM_error(FATAL, "MOM_mixedlayer_restrat: "// & @@ -250,8 +250,8 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var ! Apply time filter (to remove diurnal cycle) if (CS%MLE_MLD_decay_time>0.) then if (CS%debug) then - call hchksum(CS%MLD_filtered,'mixed_layer_restrat: MLD_filtered',G%HI,haloshift=1,scale=GV%H_to_m) - call hchksum(MLD_in,'mixed_layer_restrat: MLD in',G%HI,haloshift=1) + call hchksum(CS%MLD_filtered, 'mixed_layer_restrat: MLD_filtered', G%HI, haloshift=1, scale=GV%H_to_m) + call hchksum(MLD_in, 'mixed_layer_restrat: MLD in', G%HI, haloshift=1, scale=US%Z_to_m) endif aFac = CS%MLE_MLD_decay_time / ( dt + CS%MLE_MLD_decay_time ) bFac = dt / ( dt + CS%MLE_MLD_decay_time ) diff --git a/src/parameterizations/vertical/MOM_CVMix_conv.F90 b/src/parameterizations/vertical/MOM_CVMix_conv.F90 index 0b1abba577..06974095e1 100644 --- a/src/parameterizations/vertical/MOM_CVMix_conv.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_conv.F90 @@ -154,9 +154,9 @@ subroutine calculate_CVMix_conv(h, tv, G, GV, US, CS, hbl) type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]. type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure. - type(CVMix_conv_cs), pointer :: CS !< The control structure returned + type(CVMix_conv_cs), pointer :: CS !< The control structure returned !! by a previous call to CVMix_conv_init. - real, dimension(:,:), optional, pointer :: hbl!< Depth of ocean boundary layer [m] + real, dimension(SZI_(G),SZJ_(G)), intent(in) :: hbl !< Depth of ocean boundary layer [Z ~> m] ! local variables real, dimension(SZK_(G)) :: rho_lwr !< Adiabatic Water Density, this is a dummy !! variable since here convection is always @@ -172,7 +172,9 @@ subroutine calculate_CVMix_conv(h, tv, G, GV, US, CS, hbl) ! [Z s-2 R-1 ~> m4 s-2 kg-1] real :: pref ! Interface pressures [R L2 T-2 ~> Pa] real :: rhok, rhokm1 ! In situ densities of the layers above and below at the interface pressure [R ~> kg m-3] - real :: dz, dh, hcorr + real :: hbl_KPP ! The depth of the ocean boundary as used by KPP [m] + real :: dz ! A thickness [Z ~> m] + real :: dh, hcorr ! Two thicknesses [m] integer :: i, j, k g_o_rho0 = US%L_to_Z**2*US%s_to_T**2 * GV%g_Earth / GV%Rho0 @@ -180,11 +182,6 @@ subroutine calculate_CVMix_conv(h, tv, G, GV, US, CS, hbl) ! initialize dummy variables rho_lwr(:) = 0.0; rho_1d(:) = 0.0 - if (.not. associated(hbl)) then - allocate(hbl(SZI_(G), SZJ_(G))) - hbl(:,:) = 0.0 - endif - do j = G%jsc, G%jec do i = G%isc, G%iec @@ -213,7 +210,7 @@ subroutine calculate_CVMix_conv(h, tv, G, GV, US, CS, hbl) hcorr = 0.0 ! compute heights at cell center and interfaces do k=1,G%ke - dh = h(i,j,k) * GV%H_to_m ! Nominal thickness to use for increment + dh = h(i,j,k) * GV%H_to_m ! Nominal thickness to use for increment, in the units used by CVMix. dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0) hcorr = min( dh - CS%min_thickness, 0. ) ! If inflating then hcorr<0 dh = max( dh, CS%min_thickness ) ! Limit increment dh>=min_thickness @@ -222,7 +219,8 @@ subroutine calculate_CVMix_conv(h, tv, G, GV, US, CS, hbl) enddo ! gets index of the level and interface above hbl - kOBL = CVMix_kpp_compute_kOBL_depth(iFaceHeight, cellHeight,hbl(i,j)) + hbl_KPP = US%Z_to_m*hbl(i,j) ! Convert to the units used by CVMix. + kOBL = CVMix_kpp_compute_kOBL_depth(iFaceHeight, cellHeight, hbl_KPP) kv_col(:) = 0.0 ; kd_col(:) = 0.0 call CVMix_coeffs_conv(Mdiff_out=kv_col(:), & diff --git a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 index 28b1b8cc0b..358c7a7fa7 100644 --- a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 +++ b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 @@ -219,7 +219,7 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, C type(optics_type), pointer :: optics !< The structure containing the inverse of the !! vertical absorption decay scale for !! penetrating shortwave radiation [m-1]. - real, dimension(:,:), pointer :: Hml !< Active mixed layer depth [m]. + real, dimension(:,:), pointer :: Hml !< Active mixed layer depth [Z ~> m]. logical, intent(in) :: aggregate_FW_forcing !< If true, the net incoming and !! outgoing surface freshwater fluxes are !! combined before being applied, instead of @@ -592,7 +592,7 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, C CS%ML_depth(i,j) = h(i,0) * GV%H_to_m ! Rescale the diagnostic. enddo ; endif if (associated(Hml)) then ; do i=is,ie - Hml(i,j) = G%mask2dT(i,j) * (h(i,0) * GV%H_to_m) ! Rescale the diagnostic for output. + Hml(i,j) = G%mask2dT(i,j) * (h(i,0) * GV%H_to_Z) ! Rescale the diagnostic for output. enddo ; endif ! At this point, return water to the original layers, but constrained to diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 3c017c0f6c..d753afc97b 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -263,7 +263,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h !< thickness [H ~> m or kg m-2] type(thermo_var_ptrs), intent(inout) :: tv !< points to thermodynamic fields !! unused have NULL ptrs - real, dimension(:,:), pointer :: Hml !< Active mixed layer depth [m] + real, dimension(:,:), pointer :: Hml !< Active mixed layer depth [Z ~> m] type(forcing), intent(inout) :: fluxes !< points to forcing fields !! unused fields have NULL ptrs type(vertvisc_type), intent(inout) :: visc !< vertical viscosities, BBL properies, and @@ -451,7 +451,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h !< thickness [H ~> m or kg m-2] type(thermo_var_ptrs), intent(inout) :: tv !< points to thermodynamic fields !! unused have NULL ptrs - real, dimension(:,:), pointer :: Hml !< Active mixed layer depth [m] + real, dimension(:,:), pointer :: Hml !< Active mixed layer depth [Z ~> m] type(forcing), intent(inout) :: fluxes !< points to forcing fields !! unused fields have NULL ptrs type(vertvisc_type), intent(inout) :: visc !< vertical viscosities, BBL properies, and @@ -663,7 +663,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim Kd_salt, visc%Kv_shear, CS%KPP_NLTheat, CS%KPP_NLTscalar, Waves=Waves) if (associated(Hml)) then - call KPP_get_BLD(CS%KPP_CSp, Hml(:,:), G, US, m_to_BLD_units=1.0) + call KPP_get_BLD(CS%KPP_CSp, Hml(:,:), G, US) call pass_var(Hml, G%domain, halo=1) ! If visc%MLD exists, copy KPP's BLD into it if (associated(visc%MLD)) visc%MLD(:,:) = Hml(:,:) @@ -840,12 +840,12 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim CS%energetic_PBL_CSp, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves) if (associated(Hml)) then - call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, Hml(:,:), G, US, m_to_MLD_units=1.0) + call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, Hml(:,:), G, US) call pass_var(Hml, G%domain, halo=1) ! If visc%MLD exists, copy ePBL's MLD into it if (associated(visc%MLD)) visc%MLD(:,:) = Hml(:,:) elseif (associated(visc%MLD)) then - call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, visc%MLD, G, US, m_to_MLD_units=1.0) + call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, visc%MLD, G, US) call pass_var(visc%MLD, G%domain, halo=1) endif @@ -1232,7 +1232,7 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h !< thickness [H ~> m or kg m-2] type(thermo_var_ptrs), intent(inout) :: tv !< points to thermodynamic fields !! unused have NULL ptrs - real, dimension(:,:), pointer :: Hml !< Active mixed layer depth [m] + real, dimension(:,:), pointer :: Hml !< Active mixed layer depth [Z ~> m] type(forcing), intent(inout) :: fluxes !< points to forcing fields !! unused fields have NULL ptrs type(vertvisc_type), intent(inout) :: visc !< vertical viscosities, BBL properies, and @@ -1446,7 +1446,7 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, Kd_salt, visc%Kv_shear, CS%KPP_NLTheat, CS%KPP_NLTscalar, Waves=Waves) if (associated(Hml)) then - call KPP_get_BLD(CS%KPP_CSp, Hml(:,:), G, US, m_to_BLD_units=1.0) + call KPP_get_BLD(CS%KPP_CSp, Hml(:,:), G, US) call pass_var(Hml, G%domain, halo=1) ! If visc%MLD exists, copy KPP's BLD into it if (associated(visc%MLD)) visc%MLD(:,:) = Hml(:,:) @@ -1569,12 +1569,12 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, CS%energetic_PBL_CSp, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves) if (associated(Hml)) then - call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, Hml(:,:), G, US, m_to_MLD_units=1.0) + call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, Hml(:,:), G, US) call pass_var(Hml, G%domain, halo=1) ! If visc%MLD exists, copy ePBL's MLD into it if (associated(visc%MLD)) visc%MLD(:,:) = Hml(:,:) elseif (associated(visc%MLD)) then - call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, visc%MLD, G, US, m_to_MLD_units=1.0) + call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, visc%MLD, G, US) call pass_var(visc%MLD, G%domain, halo=1) endif @@ -1910,7 +1910,7 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h !< thickness [H ~> m or kg m-2] type(thermo_var_ptrs), intent(inout) :: tv !< points to thermodynamic fields !! unused have NULL ptrs - real, dimension(:,:), pointer :: Hml !< Active mixed layer depth [m] + real, dimension(:,:), pointer :: Hml !< Active mixed layer depth [Z ~> m] type(forcing), intent(inout) :: fluxes !< points to forcing fields !! unused fields have NULL ptrs type(vertvisc_type), intent(inout) :: visc !< vertical viscosities, BBL properies, and @@ -2179,7 +2179,7 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e Kd_salt, visc%Kv_shear, CS%KPP_NLTheat, CS%KPP_NLTscalar, Waves=Waves) if (associated(Hml)) then - call KPP_get_BLD(CS%KPP_CSp, Hml(:,:), G, US, m_to_BLD_units=1.0) + call KPP_get_BLD(CS%KPP_CSp, Hml(:,:), G, US) call pass_var(Hml, G%domain, halo=1) ! If visc%MLD exists, copy KPP's BLD into it if (associated(visc%MLD)) visc%MLD(:,:) = Hml(:,:) diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index 0046cd8b18..f208b9fe09 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -1844,7 +1844,7 @@ subroutine set_visc_init(Time, G, GV, US, param_file, diag, visc, CS, restart_CS integer :: i, j, k, is, ie, js, je, n integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz logical :: default_2018_answers - logical :: use_kappa_shear, adiabatic, use_omega + logical :: use_kappa_shear, adiabatic, use_omega, MLE_use_PBL_MLD logical :: use_CVMix_ddiff, differential_diffusion, use_KPP character(len=200) :: filename, tideamp_file type(OBC_segment_type), pointer :: segment => NULL() ! pointer to OBC segment type @@ -2047,6 +2047,9 @@ subroutine set_visc_init(Time, G, GV, US, param_file, diag, visc, CS, restart_CS if (CS%c_Smag < 0.0) CS%c_Smag = 0.15 endif + call get_param(param_file, mdl, "MLE_USE_PBL_MLD", MLE_use_PBL_MLD, & + default=.false., do_not_log=.true.) + if (CS%RiNo_mix .and. kappa_shear_at_vertex(param_file)) then ! This is necessary for reproduciblity across restarts in non-symmetric mode. call pass_var(visc%Kv_shear_Bu, G%Domain, position=CORNER, complete=.true.) @@ -2140,7 +2143,14 @@ subroutine set_visc_init(Time, G, GV, US, param_file, diag, visc, CS, restart_CS visc%Kv_shear_Bu(I,J,k) = Z2_T_rescale * visc%Kv_shear_Bu(I,J,k) enddo ; enddo ; enddo endif ; endif + endif + if (MLE_use_PBL_MLD .and. (Z_rescale /= 1.0)) then + if (associated(visc%MLD)) then ; if (query_initialized(visc%MLD, "MLD", restart_CS)) then + do j=js,je ; do i=is,ie + visc%MLD(i,j) = Z_rescale * visc%MLD(i,j) + enddo ; enddo + endif ; endif endif end subroutine set_visc_init diff --git a/src/tracer/MOM_generic_tracer.F90 b/src/tracer/MOM_generic_tracer.F90 index e68833c3cd..b198db3e32 100644 --- a/src/tracer/MOM_generic_tracer.F90 +++ b/src/tracer/MOM_generic_tracer.F90 @@ -388,8 +388,8 @@ end subroutine initialize_MOM_generic_tracer !! flux as a source. subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml, dt, G, GV, CS, tv, optics, & evap_CFL_limit, minimum_forcing_depth) - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure - type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + 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_(G)), & intent(in) :: h_old !< Layer thickness before entrainment [H ~> m or kg m-2]. real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & @@ -402,16 +402,16 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml, !! below during this call [H ~> m or kg m-2]. type(forcing), intent(in) :: fluxes !< A structure containing pointers to thermodynamic !! and tracer forcing fields. - real, dimension(SZI_(G),SZJ_(G)), intent(in) :: Hml !< Mixed layer depth [H ~> m or kg m-2] - real, intent(in) :: dt !< The amount of time covered by this call [s] + real, dimension(SZI_(G),SZJ_(G)), intent(in) :: Hml !< Mixed layer depth [Z ~> m] + real, intent(in) :: dt !< The amount of time covered by this call [s] type(MOM_generic_tracer_CS), pointer :: CS !< Pointer to the control structure for this module. - type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables + type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables type(optics_type), intent(in) :: optics !< The structure containing optical properties. real, optional, intent(in) :: evap_CFL_limit !< Limits how much water can be fluxed out of - !! the top layer Stored previously in diabatic CS. + !! the top layer Stored previously in diabatic CS. real, optional, intent(in) :: minimum_forcing_depth !< The smallest depth over which fluxes - !! can be applied [H ~> m or kg m-2] - ! Stored previously in diabatic CS. + !! can be applied [H ~> m or kg m-2] + ! Stored previously in diabatic CS. ! The arguments to this subroutine are redundant in that ! h_new(k) = h_old(k) + ea(k) - eb(k-1) + eb(k) - ea(k+1) @@ -423,6 +423,7 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml, real, dimension(:,:), pointer :: stf_array,trunoff_array,runoff_tracer_flux_array real :: surface_field(SZI_(G),SZJ_(G)) + real :: dz_ml(SZI_(G),SZJ_(G)) ! The mixed layer depth in the MKS units used for generic tracers [m] real :: sosga real, dimension(G%isd:G%ied,G%jsd:G%jed,G%ke) :: rho_dzt, dzt @@ -483,9 +484,10 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml, do k = 1, nk ; do j = jsc, jec ; do i = isc, iec !{ dzt(i,j,k) = GV%H_to_m * h_old(i,j,k) enddo ; enddo ; enddo !} - + dz_ml(:,:) = 0.0 do j=jsc,jec ; do i=isc,iec - surface_field(i,j) = tv%S(i,j,1) + surface_field(i,j) = tv%S(i,j,1) + dz_ml(i,j) = G%US%Z_to_m * Hml enddo ; enddo sosga = global_area_mean(surface_field, G) @@ -494,12 +496,12 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml, ! if ((G%US%L_to_m == 1.0) .and. (G%US%RZ_to_kg_m2 == 1.0) .and. (G%US%s_to_T == 1.0)) then ! Avoid unnecessary copies when no unit conversion is needed. - call generic_tracer_source(tv%T, tv%S, rho_dzt, dzt, Hml, G%isd, G%jsd, 1, dt, & + call generic_tracer_source(tv%T, tv%S, rho_dzt, dzt, dz_ml, G%isd, G%jsd, 1, dt, & G%areaT, get_diag_time_end(CS%diag), & optics%nbands, optics%max_wavelength_band, optics%sw_pen_band, optics%opacity_band, & internal_heat=tv%internal_heat, frunoff=fluxes%frunoff, sosga=sosga) else - call generic_tracer_source(tv%T, tv%S, rho_dzt, dzt, Hml, G%isd, G%jsd, 1, dt, & + call generic_tracer_source(tv%T, tv%S, rho_dzt, dzt, dz_ml, G%isd, G%jsd, 1, dt, & G%US%L_to_m**2*G%areaT(:,:), get_diag_time_end(CS%diag), & optics%nbands, optics%max_wavelength_band, optics%sw_pen_band, optics%opacity_band, & internal_heat=G%US%RZ_to_kg_m2*tv%internal_heat(:,:), & diff --git a/src/tracer/MOM_offline_main.F90 b/src/tracer/MOM_offline_main.F90 index dfdcb4c09b..65f83ecfea 100644 --- a/src/tracer/MOM_offline_main.F90 +++ b/src/tracer/MOM_offline_main.F90 @@ -176,7 +176,7 @@ module MOM_offline_main real, allocatable, dimension(:,:) :: netMassIn !< Freshwater fluxes into the ocean real, allocatable, dimension(:,:) :: netMassOut !< Freshwater fluxes out of the ocean - real, allocatable, dimension(:,:) :: mld !< Mixed layer depths at thickness points [H ~> m or kg m-2]. + real, allocatable, dimension(:,:) :: mld !< Mixed layer depths at thickness points [Z ~> m]. ! Allocatable arrays to read in entire fields during initialization real, allocatable, dimension(:,:,:,:) :: uhtr_all !< Entire field of zonal transport diff --git a/src/tracer/MOM_tracer_flow_control.F90 b/src/tracer/MOM_tracer_flow_control.F90 index 86003605f7..a9bf9a03d9 100644 --- a/src/tracer/MOM_tracer_flow_control.F90 +++ b/src/tracer/MOM_tracer_flow_control.F90 @@ -417,7 +417,7 @@ subroutine call_tracer_column_fns(h_old, h_new, ea, eb, fluxes, Hml, dt, G, GV, type(forcing), intent(in) :: fluxes !< A structure containing pointers to !! any possible forcing fields. !! Unused fields have NULL ptrs. - real, dimension(NIMEM_,NJMEM_), intent(in) :: Hml !< Mixed layer depth [H ~> m or kg m-2] + real, dimension(NIMEM_,NJMEM_), intent(in) :: Hml !< Mixed layer depth [Z ~> m] real, intent(in) :: dt !< The amount of time covered by this !! call [T ~> s] type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. From bee6be783bbc996eec72016615f4079221f5a208 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 1 May 2020 10:20:42 -0400 Subject: [PATCH 10/11] +Rescaled sfc_state%frazil and other variables Rescaled the sea_lev, frazil, melt_potential, taux_shelf, tauy_shelf, TempxPmE, salt_deficit, ocean_mass, ocean_heat, ocean_salt and internal_heat elements of the surface state type. This change includes the addition of a unit_scale_type argument to MOM_surface_chksum. Several related internal variables were rescaled as well or otherwise made less idiosyncratic, notably in extract_surface_state, MOM_forcing_chksum, forcing_diagnostics. All solutions are bitwise identical, but there are changes to the units of transparent types and a new argument to a publicly visible subroutine. --- config_src/coupled_driver/ocean_model_MOM.F90 | 10 +-- config_src/mct_driver/mom_ocean_model_mct.F90 | 14 ++-- .../nuopc_driver/mom_ocean_model_nuopc.F90 | 14 ++-- src/core/MOM.F90 | 82 +++++++++--------- src/core/MOM_PressureForce_analytic_FV.F90 | 2 +- src/core/MOM_checksum_packages.F90 | 17 ++-- src/core/MOM_forcing_type.F90 | 84 +++++++++---------- src/core/MOM_variables.F90 | 24 +++--- src/ice_shelf/MOM_ice_shelf.F90 | 36 ++++---- src/ice_shelf/MOM_marine_ice.F90 | 21 +++-- 10 files changed, 148 insertions(+), 156 deletions(-) diff --git a/config_src/coupled_driver/ocean_model_MOM.F90 b/config_src/coupled_driver/ocean_model_MOM.F90 index 28ac193d8d..05759cb7b8 100644 --- a/config_src/coupled_driver/ocean_model_MOM.F90 +++ b/config_src/coupled_driver/ocean_model_MOM.F90 @@ -862,19 +862,19 @@ subroutine convert_state_to_ocean_type(sfc_state, Ocean_sfc, G, US, patm, press_ if (present(patm)) then do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd - Ocean_sfc%sea_lev(i,j) = sfc_state%sea_lev(i+i0,j+j0) + patm(i,j) * press_to_z - Ocean_sfc%area(i,j) = US%L_to_m**2*G%areaT(i+i0,j+j0) + Ocean_sfc%sea_lev(i,j) = US%Z_to_m * sfc_state%sea_lev(i+i0,j+j0) + patm(i,j) * press_to_z + Ocean_sfc%area(i,j) = US%L_to_m**2 * G%areaT(i+i0,j+j0) enddo ; enddo else do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd - Ocean_sfc%sea_lev(i,j) = sfc_state%sea_lev(i+i0,j+j0) - Ocean_sfc%area(i,j) = US%L_to_m**2*G%areaT(i+i0,j+j0) + Ocean_sfc%sea_lev(i,j) = US%Z_to_m * sfc_state%sea_lev(i+i0,j+j0) + Ocean_sfc%area(i,j) = US%L_to_m**2 * G%areaT(i+i0,j+j0) enddo ; enddo endif if (allocated(sfc_state%frazil)) then do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd - Ocean_sfc%frazil(i,j) = sfc_state%frazil(i+i0,j+j0) + Ocean_sfc%frazil(i,j) = US%Q_to_J_kg*US%RZ_to_kg_m2 * sfc_state%frazil(i+i0,j+j0) enddo ; enddo endif diff --git a/config_src/mct_driver/mom_ocean_model_mct.F90 b/config_src/mct_driver/mom_ocean_model_mct.F90 index 1d1e9fa888..f8a4a19532 100644 --- a/config_src/mct_driver/mom_ocean_model_mct.F90 +++ b/config_src/mct_driver/mom_ocean_model_mct.F90 @@ -894,31 +894,31 @@ subroutine convert_state_to_ocean_type(sfc_state, Ocean_sfc, G, US, patm, press_ if (present(patm)) then do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd - Ocean_sfc%sea_lev(i,j) = sfc_state%sea_lev(i+i0,j+j0) + patm(i,j) * press_to_z - Ocean_sfc%area(i,j) = US%L_to_m**2*G%areaT(i+i0,j+j0) + Ocean_sfc%sea_lev(i,j) = US%Z_to_m * sfc_state%sea_lev(i+i0,j+j0) + patm(i,j) * press_to_z + Ocean_sfc%area(i,j) = US%L_to_m**2 * G%areaT(i+i0,j+j0) enddo ; enddo else do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd - Ocean_sfc%sea_lev(i,j) = sfc_state%sea_lev(i+i0,j+j0) - Ocean_sfc%area(i,j) = US%L_to_m**2*G%areaT(i+i0,j+j0) + Ocean_sfc%sea_lev(i,j) = US%Z_to_m * sfc_state%sea_lev(i+i0,j+j0) + Ocean_sfc%area(i,j) = US%L_to_m**2 * G%areaT(i+i0,j+j0) enddo ; enddo endif if (allocated(sfc_state%frazil)) then do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd - Ocean_sfc%frazil(i,j) = sfc_state%frazil(i+i0,j+j0) + Ocean_sfc%frazil(i,j) = US%Q_to_J_kg*US%RZ_to_kg_m2 * sfc_state%frazil(i+i0,j+j0) enddo ; enddo endif if (allocated(sfc_state%melt_potential)) then do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd - Ocean_sfc%melt_potential(i,j) = sfc_state%melt_potential(i+i0,j+j0) + Ocean_sfc%melt_potential(i,j) = US%Q_to_J_kg*US%RZ_to_kg_m2 * sfc_state%melt_potential(i+i0,j+j0) enddo ; enddo endif if (allocated(sfc_state%Hml)) then do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd - Ocean_sfc%OBLD(i,j) = US%Z_to_m*sfc_state%Hml(i+i0,j+j0) + Ocean_sfc%OBLD(i,j) = US%Z_to_m * sfc_state%Hml(i+i0,j+j0) enddo ; enddo endif diff --git a/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 b/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 index 22a4c7eaa2..9946aec4f9 100644 --- a/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 +++ b/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 @@ -889,31 +889,31 @@ subroutine convert_state_to_ocean_type(sfc_state, Ocean_sfc, G, US, patm, press_ if (present(patm)) then do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd - Ocean_sfc%sea_lev(i,j) = sfc_state%sea_lev(i+i0,j+j0) + patm(i,j) * press_to_z - Ocean_sfc%area(i,j) = US%L_to_m**2*G%areaT(i+i0,j+j0) + Ocean_sfc%sea_lev(i,j) = US%Z_to_m * sfc_state%sea_lev(i+i0,j+j0) + patm(i,j) * press_to_z + Ocean_sfc%area(i,j) = US%L_to_m**2 * G%areaT(i+i0,j+j0) enddo ; enddo else do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd - Ocean_sfc%sea_lev(i,j) = sfc_state%sea_lev(i+i0,j+j0) - Ocean_sfc%area(i,j) = US%L_to_m**2*G%areaT(i+i0,j+j0) + Ocean_sfc%sea_lev(i,j) = US%Z_to_m * sfc_state%sea_lev(i+i0,j+j0) + Ocean_sfc%area(i,j) = US%L_to_m**2 * G%areaT(i+i0,j+j0) enddo ; enddo endif if (allocated(sfc_state%frazil)) then do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd - Ocean_sfc%frazil(i,j) = sfc_state%frazil(i+i0,j+j0) + Ocean_sfc%frazil(i,j) = US%Q_to_J_kg*US%RZ_to_kg_m2 * sfc_state%frazil(i+i0,j+j0) enddo ; enddo endif if (allocated(sfc_state%melt_potential)) then do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd - Ocean_sfc%melt_potential(i,j) = sfc_state%melt_potential(i+i0,j+j0) + Ocean_sfc%melt_potential(i,j) = US%Q_to_J_kg*US%RZ_to_kg_m2 * sfc_state%melt_potential(i+i0,j+j0) enddo ; enddo endif if (allocated(sfc_state%Hml)) then do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd - Ocean_sfc%OBLD(i,j) = US%Z_to_m*sfc_state%Hml(i+i0,j+j0) + Ocean_sfc%OBLD(i,j) = US%Z_to_m * sfc_state%Hml(i+i0,j+j0) enddo ; enddo endif diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index e3bd314273..aff4860a21 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -296,20 +296,20 @@ module MOM !! average surface tracer properties when a bulk !! mixed layer is not used [Z ~> m], or a negative value !! if a bulk mixed layer is being used. - real :: HFrz !< If HFrz > 0, melt potential will be computed. - !! The actual depth over which melt potential is computed will - !! min(HFrz, OBLD), where OBLD is the boundary layer depth. + real :: HFrz !< If HFrz > 0, the nominal depth over which melt potential is + !! computed [Z ~> m]. The actual depth over which melt potential is + !! computed is min(HFrz, OBLD), where OBLD is the boundary layer depth. !! If HFrz <= 0 (default), melt potential will not be computed. real :: Hmix_UV !< Depth scale over which to average surface flow to !! feedback to the coupler/driver [Z ~> m] when !! bulk mixed layer is not used, or a negative value !! if a bulk mixed layer is being used. logical :: check_bad_sfc_vals !< If true, scan surface state for ridiculous values. - real :: bad_val_ssh_max !< Maximum SSH before triggering bad value message [m] + real :: bad_val_ssh_max !< Maximum SSH before triggering bad value message [Z ~> m] real :: bad_val_sst_max !< Maximum SST before triggering bad value message [degC] real :: bad_val_sst_min !< Minimum SST before triggering bad value message [degC] real :: bad_val_sss_max !< Maximum SSS before triggering bad value message [ppt] - real :: bad_val_col_thick !< Minimum column thickness before triggering bad value message [m] + real :: bad_val_col_thick !< Minimum column thickness before triggering bad value message [Z ~> m] logical :: answers_2018 !< If true, use expressions for the surface properties that recover !! the answers from the end of 2018. Otherwise, use more appropriate !! expressions that differ at roundoff for non-Boussinsq cases. @@ -1859,7 +1859,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & "If HFREEZE > 0, melt potential will be computed. The actual depth "//& "over which melt potential is computed will be min(HFREEZE, OBLD), "//& "where OBLD is the boundary layer depth. If HFREEZE <= 0 (default), "//& - "melt potential will not be computed.", units="m", default=-1.0) + "melt potential will not be computed.", units="m", default=-1.0, scale=US%m_to_Z) call get_param(param_file, "MOM", "INTERPOLATE_P_SURF", CS%interp_p_surf, & "If true, linearly interpolate the surface pressure "//& "over the coupling time step, using the specified value "//& @@ -1944,8 +1944,8 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & if (CS%check_bad_sfc_vals) then call get_param(param_file, "MOM", "BAD_VAL_SSH_MAX", CS%bad_val_ssh_max, & "The value of SSH above which a bad value message is "//& - "triggered, if CHECK_BAD_SURFACE_VALS is true.", units="m", & - default=20.0) + "triggered, if CHECK_BAD_SURFACE_VALS is true.", & + units="m", default=20.0, scale=US%m_to_Z) call get_param(param_file, "MOM", "BAD_VAL_SSS_MAX", CS%bad_val_sss_max, & "The value of SSS above which a bad value message is "//& "triggered, if CHECK_BAD_SURFACE_VALS is true.", units="PPT", & @@ -1960,8 +1960,8 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & units="deg C", default=-2.1) call get_param(param_file, "MOM", "BAD_VAL_COLUMN_THICKNESS", CS%bad_val_col_thick, & "The value of column thickness below which a bad value message is "//& - "triggered, if CHECK_BAD_SURFACE_VALS is true.", units="m", & - default=0.0) + "triggered, if CHECK_BAD_SURFACE_VALS is true.", & + units="m", default=0.0, scale=US%m_to_Z) endif call get_param(param_file, "MOM", "DEFAULT_2018_ANSWERS", default_2018_answers, & "This sets the default value for the various _2018_ANSWERS parameters.", & @@ -3000,15 +3000,14 @@ subroutine extract_surface_state(CS, sfc_state_in) real :: depth_ml !< Depth over which to average to determine mixed !! layer properties [Z ~> m] or [H ~> m or kg m-2] real :: dh !< Thickness of a layer within the mixed layer [Z ~> m] or [H ~> m or kg m-2] - real :: mass !< Mass per unit area of a layer [kg m-2] - real :: bathy_m !< The depth of bathymetry [m] (not Z), used for error checking. + real :: mass !< Mass per unit area of a layer [R Z ~> kg m-2] real :: T_freeze !< freezing temperature [degC] real :: I_depth !< The inverse of depth [Z-1 ~> m-1] or [H-1 ~> m-1 or m2 kg-1] real :: missing_depth !< The portion of depth_ml that can not be found in a column [H ~> m or kg m-2] real :: H_rescale !< A conversion factor from thickness units to the units used in the !! calculation of properties of the uppermost ocean [nondim] or [Z H-1 ~> 1 or m3 kg-1] ! After the ANSWERS_2018 flag has been obsoleted, H_rescale will be 1. - real :: delT(SZI_(CS%G)) !< Depth integral of T-T_freeze [m degC] + real :: delT(SZI_(CS%G)) !< Depth integral of T-T_freeze [Z degC ~> m degC] logical :: use_temperature !< If true, temp and saln used as state variables. integer :: i, j, k, is, ie, js, je, nz, numberOfErrors, ig, jg integer :: isd, ied, jsd, jed @@ -3049,11 +3048,11 @@ subroutine extract_surface_state(CS, sfc_state_in) sfc_state%S_is_absS = CS%tv%S_is_absS do j=js,je ; do i=is,ie - sfc_state%sea_lev(i,j) = CS%ave_ssh_ibc(i,j) + sfc_state%sea_lev(i,j) = US%m_to_Z*CS%ave_ssh_ibc(i,j) enddo ; enddo if (allocated(sfc_state%frazil) .and. associated(CS%tv%frazil)) then ; do j=js,je ; do i=is,ie - sfc_state%frazil(i,j) = US%Q_to_J_kg*US%RZ_to_kg_m2 * CS%tv%frazil(i,j) + sfc_state%frazil(i,j) = CS%tv%frazil(i,j) enddo ; enddo ; endif ! copy Hml into sfc_state, so that caps can access it @@ -3209,7 +3208,7 @@ subroutine extract_surface_state(CS, sfc_state_in) if (allocated(sfc_state%melt_potential)) then - !$OMP parallel do default(shared) private(depth_ml, dh, T_freeze, depth, delT) + !$OMP parallel do default(shared) private(depth_ml, dh, T_freeze, depth, delT) do j=js,je do i=is,ie depth(i) = 0.0 @@ -3217,9 +3216,9 @@ subroutine extract_surface_state(CS, sfc_state_in) enddo do k=1,nz ; do i=is,ie - depth_ml = min(CS%HFrz, US%Z_to_m*CS%visc%MLD(i,j)) - if (depth(i) + h(i,j,k)*GV%H_to_m < depth_ml) then - dh = h(i,j,k)*GV%H_to_m + depth_ml = min(CS%HFrz, CS%visc%MLD(i,j)) + if (depth(i) + h(i,j,k)*GV%H_to_Z < depth_ml) then + dh = h(i,j,k)*GV%H_to_Z elseif (depth(i) < depth_ml) then dh = depth_ml - depth(i) else @@ -3237,8 +3236,8 @@ subroutine extract_surface_state(CS, sfc_state_in) sfc_state%melt_potential(i,j) = 0.0 if (G%mask2dT(i,j)>0.) then - ! instantaneous melt_potential [J m-2] - sfc_state%melt_potential(i,j) = US%Q_to_J_kg*US%R_to_kg_m3 * CS%tv%C_p * GV%Rho0 * delT(i) + ! instantaneous melt_potential [Q R Z ~> J m-2] + sfc_state%melt_potential(i,j) = CS%tv%C_p * GV%Rho0 * delT(i) endif enddo enddo ! end of j loop @@ -3248,31 +3247,31 @@ subroutine extract_surface_state(CS, sfc_state_in) !$OMP parallel do default(shared) do j=js,je ; do i=is,ie ! Convert from gSalt to kgSalt - sfc_state%salt_deficit(i,j) = 0.001 * US%RZ_to_kg_m2*CS%tv%salt_deficit(i,j) + sfc_state%salt_deficit(i,j) = 0.001 * CS%tv%salt_deficit(i,j) enddo ; enddo endif if (allocated(sfc_state%TempxPmE) .and. associated(CS%tv%TempxPmE)) then !$OMP parallel do default(shared) do j=js,je ; do i=is,ie - sfc_state%TempxPmE(i,j) = US%RZ_to_kg_m2*CS%tv%TempxPmE(i,j) + sfc_state%TempxPmE(i,j) = CS%tv%TempxPmE(i,j) enddo ; enddo endif if (allocated(sfc_state%internal_heat) .and. associated(CS%tv%internal_heat)) then !$OMP parallel do default(shared) do j=js,je ; do i=is,ie - sfc_state%internal_heat(i,j) = US%RZ_to_kg_m2*CS%tv%internal_heat(i,j) + sfc_state%internal_heat(i,j) = CS%tv%internal_heat(i,j) enddo ; enddo endif if (allocated(sfc_state%taux_shelf) .and. associated(CS%visc%taux_shelf)) then !$OMP parallel do default(shared) do j=js,je ; do I=is-1,ie - sfc_state%taux_shelf(I,j) = US%RZ_T_to_kg_m2s*US%L_T_to_m_s*CS%visc%taux_shelf(I,j) + sfc_state%taux_shelf(I,j) = CS%visc%taux_shelf(I,j) enddo ; enddo endif if (allocated(sfc_state%tauy_shelf) .and. associated(CS%visc%tauy_shelf)) then !$OMP parallel do default(shared) do J=js-1,je ; do i=is,ie - sfc_state%tauy_shelf(i,J) = US%RZ_T_to_kg_m2s*US%L_T_to_m_s*CS%visc%tauy_shelf(i,J) + sfc_state%tauy_shelf(i,J) = CS%visc%tauy_shelf(i,J) enddo ; enddo endif @@ -3285,11 +3284,10 @@ subroutine extract_surface_state(CS, sfc_state_in) enddo ; enddo !$OMP parallel do default(shared) private(mass) do j=js,je ; do k=1,nz; do i=is,ie - mass = GV%H_to_kg_m2*h(i,j,k) + mass = GV%H_to_RZ*h(i,j,k) sfc_state%ocean_mass(i,j) = sfc_state%ocean_mass(i,j) + mass - sfc_state%ocean_heat(i,j) = sfc_state%ocean_heat(i,j) + mass*CS%tv%T(i,j,k) - sfc_state%ocean_salt(i,j) = sfc_state%ocean_salt(i,j) + & - mass * (1.0e-3*CS%tv%S(i,j,k)) + sfc_state%ocean_heat(i,j) = sfc_state%ocean_heat(i,j) + mass * CS%tv%T(i,j,k) + sfc_state%ocean_salt(i,j) = sfc_state%ocean_salt(i,j) + mass * (1.0e-3*CS%tv%S(i,j,k)) enddo ; enddo ; enddo else if (allocated(sfc_state%ocean_mass)) then @@ -3297,7 +3295,7 @@ subroutine extract_surface_state(CS, sfc_state_in) do j=js,je ; do i=is,ie ; sfc_state%ocean_mass(i,j) = 0.0 ; enddo ; enddo !$OMP parallel do default(shared) do j=js,je ; do k=1,nz ; do i=is,ie - sfc_state%ocean_mass(i,j) = sfc_state%ocean_mass(i,j) + GV%H_to_kg_m2*h(i,j,k) + sfc_state%ocean_mass(i,j) = sfc_state%ocean_mass(i,j) + GV%H_to_RZ*h(i,j,k) enddo ; enddo ; enddo endif if (allocated(sfc_state%ocean_heat)) then @@ -3305,7 +3303,7 @@ subroutine extract_surface_state(CS, sfc_state_in) do j=js,je ; do i=is,ie ; sfc_state%ocean_heat(i,j) = 0.0 ; enddo ; enddo !$OMP parallel do default(shared) private(mass) do j=js,je ; do k=1,nz ; do i=is,ie - mass = GV%H_to_kg_m2*h(i,j,k) + mass = GV%H_to_RZ*h(i,j,k) sfc_state%ocean_heat(i,j) = sfc_state%ocean_heat(i,j) + mass*CS%tv%T(i,j,k) enddo ; enddo ; enddo endif @@ -3314,9 +3312,8 @@ subroutine extract_surface_state(CS, sfc_state_in) do j=js,je ; do i=is,ie ; sfc_state%ocean_salt(i,j) = 0.0 ; enddo ; enddo !$OMP parallel do default(shared) private(mass) do j=js,je ; do k=1,nz ; do i=is,ie - mass = GV%H_to_kg_m2*h(i,j,k) - sfc_state%ocean_salt(i,j) = sfc_state%ocean_salt(i,j) + & - mass * (1.0e-3*CS%tv%S(i,j,k)) + mass = GV%H_to_RZ*h(i,j,k) + sfc_state%ocean_salt(i,j) = sfc_state%ocean_salt(i,j) + mass * (1.0e-3*CS%tv%S(i,j,k)) enddo ; enddo ; enddo endif endif @@ -3329,11 +3326,10 @@ subroutine extract_surface_state(CS, sfc_state_in) numberOfErrors=0 ! count number of errors do j=js,je; do i=is,ie if (G%mask2dT(i,j)>0.) then - bathy_m = CS%US%Z_to_m * G%bathyT(i,j) - localError = sfc_state%sea_lev(i,j)<=-bathy_m & - .or. sfc_state%sea_lev(i,j)>= CS%bad_val_ssh_max & - .or. sfc_state%sea_lev(i,j)<=-CS%bad_val_ssh_max & - .or. sfc_state%sea_lev(i,j) + bathy_m < CS%bad_val_col_thick + localError = sfc_state%sea_lev(i,j) <= -G%bathyT(i,j) & + .or. sfc_state%sea_lev(i,j) >= CS%bad_val_ssh_max & + .or. sfc_state%sea_lev(i,j) <= -CS%bad_val_ssh_max & + .or. sfc_state%sea_lev(i,j) + G%bathyT(i,j) < CS%bad_val_col_thick if (use_temperature) localError = localError & .or. sfc_state%SSS(i,j)<0. & .or. sfc_state%SSS(i,j)>=CS%bad_val_sss_max & @@ -3349,7 +3345,7 @@ subroutine extract_surface_state(CS, sfc_state_in) 'Extreme surface sfc_state detected: i=',ig,'j=',jg, & 'lon=',G%geoLonT(i,j), 'lat=',G%geoLatT(i,j), & 'x=',G%gridLonT(ig), 'y=',G%gridLatT(jg), & - 'D=',bathy_m, 'SSH=',sfc_state%sea_lev(i,j), & + 'D=',CS%US%Z_to_m*G%bathyT(i,j), 'SSH=',CS%US%Z_to_m*sfc_state%sea_lev(i,j), & 'SST=',sfc_state%SST(i,j), 'SSS=',sfc_state%SSS(i,j), & 'U-=',US%L_T_to_m_s*sfc_state%u(I-1,j), 'U+=',US%L_T_to_m_s*sfc_state%u(I,j), & 'V-=',US%L_T_to_m_s*sfc_state%v(i,J-1), 'V+=',US%L_T_to_m_s*sfc_state%v(i,J) @@ -3358,7 +3354,7 @@ subroutine extract_surface_state(CS, sfc_state_in) 'Extreme surface sfc_state detected: i=',ig,'j=',jg, & 'lon=',G%geoLonT(i,j), 'lat=',G%geoLatT(i,j), & 'x=',G%gridLonT(i), 'y=',G%gridLatT(j), & - 'D=',bathy_m, 'SSH=',sfc_state%sea_lev(i,j), & + 'D=',CS%US%Z_to_m*G%bathyT(i,j), 'SSH=',CS%US%Z_to_m*sfc_state%sea_lev(i,j), & 'U-=',US%L_T_to_m_s*sfc_state%u(I-1,j), 'U+=',US%L_T_to_m_s*sfc_state%u(I,j), & 'V-=',US%L_T_to_m_s*sfc_state%v(i,J-1), 'V+=',US%L_T_to_m_s*sfc_state%v(i,J) endif @@ -3377,7 +3373,7 @@ subroutine extract_surface_state(CS, sfc_state_in) endif endif - if (CS%debug) call MOM_surface_chksum("Post extract_sfc", sfc_state, G) + if (CS%debug) call MOM_surface_chksum("Post extract_sfc", sfc_state, G, US) ! Rotate sfc_state back onto the input grid, sfc_state_in if (CS%rotate_index) then diff --git a/src/core/MOM_PressureForce_analytic_FV.F90 b/src/core/MOM_PressureForce_analytic_FV.F90 index 614bf3bc8a..ecadb5c4b9 100644 --- a/src/core/MOM_PressureForce_analytic_FV.F90 +++ b/src/core/MOM_PressureForce_analytic_FV.F90 @@ -142,7 +142,7 @@ subroutine PressureForce_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p dM, & ! The barotropic adjustment to the Montgomery potential to ! account for a reduced gravity model [L2 T-2 ~> m2 s-2]. za ! The geopotential anomaly (i.e. g*e + alpha_0*pressure) at the - ! interface atop a layer [m2 s-2]. + ! interface atop a layer [L2 T-2 ~> m2 s-2]. real, dimension(SZI_(G)) :: Rho_cv_BL ! The coordinate potential density in the deepest variable ! density near-surface layer [R ~> kg m-3]. diff --git a/src/core/MOM_checksum_packages.F90 b/src/core/MOM_checksum_packages.F90 index 47369cf474..70ba32644f 100644 --- a/src/core/MOM_checksum_packages.F90 +++ b/src/core/MOM_checksum_packages.F90 @@ -132,7 +132,7 @@ subroutine MOM_thermo_chksum(mesg, tv, G, US, haloshift) if (associated(tv%T)) call hchksum(tv%T, mesg//" T", G%HI, haloshift=hs) if (associated(tv%S)) call hchksum(tv%S, mesg//" S", G%HI, haloshift=hs) if (associated(tv%frazil)) call hchksum(tv%frazil, mesg//" frazil", G%HI, haloshift=hs, & - scale=G%US%Q_to_J_kg*G%US%R_to_kg_m3*G%US%Z_to_m) + scale=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, scale=US%RZ_to_kg_m2) @@ -141,12 +141,13 @@ end subroutine MOM_thermo_chksum ! ============================================================================= !> Write out chksums for the ocean surface variables. -subroutine MOM_surface_chksum(mesg, sfc_state, G, haloshift, symmetric) +subroutine MOM_surface_chksum(mesg, sfc_state, G, US, haloshift, symmetric) character(len=*), intent(in) :: mesg !< A message that appears on the chksum lines. type(surface), intent(inout) :: sfc_state !< transparent ocean surface state structure !! shared with the calling routine data in this !! structure is intent out. type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0). logical, optional, intent(in) :: symmetric !< If true, do checksums on the fully symmetric !! computational domain. @@ -159,15 +160,17 @@ subroutine MOM_surface_chksum(mesg, sfc_state, G, haloshift, symmetric) if (allocated(sfc_state%SST)) call hchksum(sfc_state%SST, mesg//" SST", G%HI, haloshift=hs) if (allocated(sfc_state%SSS)) call hchksum(sfc_state%SSS, mesg//" SSS", G%HI, haloshift=hs) - if (allocated(sfc_state%sea_lev)) call hchksum(sfc_state%sea_lev, mesg//" sea_lev", G%HI, haloshift=hs) + if (allocated(sfc_state%sea_lev)) call hchksum(sfc_state%sea_lev, mesg//" sea_lev", G%HI, & + haloshift=hs, scale=US%Z_to_m) if (allocated(sfc_state%Hml)) call hchksum(sfc_state%Hml, mesg//" Hml", G%HI, haloshift=hs, & - scale=G%US%Z_to_m) + scale=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=G%US%L_T_to_m_s) + scale=US%L_T_to_m_s) ! if (allocated(sfc_state%salt_deficit)) & -! call hchksum(sfc_state%salt_deficit, mesg//" salt deficit", G%HI, haloshift=hs) - if (allocated(sfc_state%frazil)) call hchksum(sfc_state%frazil, mesg//" frazil", G%HI, haloshift=hs) +! call hchksum(sfc_state%salt_deficit, mesg//" salt deficit", G%HI, haloshift=hs, scale=US%RZ_to_kg_m2) + 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) end subroutine MOM_surface_chksum diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index aea62826e3..0ff9a4b287 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -97,7 +97,7 @@ module MOM_forcing_type latent_fprec_diag => NULL(), & !< latent [Q R Z T-1 ~> W m-2] from melting fprec (typically < 0) latent_frunoff_diag => NULL() !< latent [Q R Z T-1 ~> W m-2] from melting frunoff (calving) (typically < 0) - ! water mass fluxes into the ocean [kg m-2 s-1]; these fluxes impact the ocean mass + ! water mass fluxes into the ocean [R Z T-1 ~> kg m-2 s-1]; these fluxes impact the ocean mass real, pointer, dimension(:,:) :: & evap => NULL(), & !< (-1)*fresh water flux evaporated out of the ocean [R Z T-1 ~> kg m-2 s-1] lprec => NULL(), & !< precipitating liquid water into the ocean [R Z T-1 ~> kg m-2 s-1] @@ -672,7 +672,7 @@ subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, & &" exceeds total shortwave of ",1pe17.10,& &" at ",1pg11.4,"E, "1pg11.4,"N.")') & Pen_SW_tot(i), I_Cp_Hconvert*scale*dt * fluxes%sw(i,j), & - G%geoLonT(i,j),G%geoLatT(i,j) + G%geoLonT(i,j), G%geoLatT(i,j) call MOM_error(WARNING,mesg) endif endif @@ -1025,12 +1025,10 @@ subroutine MOM_forcing_chksum(mesg, fluxes, G, US, haloshift) type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type integer, optional, intent(in) :: haloshift !< shift in halo - real :: RZ_T_conversion ! A combination of scaling factors for mass fluxes [kg T m-2 s-1 R-1 Z-1 ~> 1] integer :: is, ie, js, je, nz, hshift is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke hshift = 1 ; if (present(haloshift)) hshift = haloshift - RZ_T_conversion = US%RZ_T_to_kg_m2s ! 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 @@ -1040,7 +1038,7 @@ subroutine MOM_forcing_chksum(mesg, fluxes, G, US, haloshift) 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) 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, scale=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) if (associated(fluxes%sw_vis_dif)) & @@ -1063,36 +1061,36 @@ subroutine MOM_forcing_chksum(mesg, fluxes, G, US, haloshift) call hchksum(fluxes%latent_frunoff_diag, mesg//" fluxes%latent_frunoff_diag", G%HI, & haloshift=hshift, scale=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, scale=US%QRZ_T_to_W_m2) if (associated(fluxes%evap)) & - call hchksum(fluxes%evap, mesg//" fluxes%evap",G%HI,haloshift=hshift, scale=RZ_T_conversion) + call hchksum(fluxes%evap, mesg//" fluxes%evap", G%HI, haloshift=hshift, scale=US%RZ_T_to_kg_m2s) if (associated(fluxes%lprec)) & - call hchksum(fluxes%lprec, mesg//" fluxes%lprec",G%HI,haloshift=hshift, scale=RZ_T_conversion) + call hchksum(fluxes%lprec, mesg//" fluxes%lprec", G%HI, haloshift=hshift, scale=US%RZ_T_to_kg_m2s) if (associated(fluxes%fprec)) & - call hchksum(fluxes%fprec, mesg//" fluxes%fprec",G%HI,haloshift=hshift, scale=RZ_T_conversion) + call hchksum(fluxes%fprec, mesg//" fluxes%fprec", G%HI, haloshift=hshift, scale=US%RZ_T_to_kg_m2s) if (associated(fluxes%vprec)) & - call hchksum(fluxes%vprec, mesg//" fluxes%vprec",G%HI,haloshift=hshift, scale=RZ_T_conversion) + call hchksum(fluxes%vprec, mesg//" fluxes%vprec", G%HI, haloshift=hshift, scale=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=RZ_T_conversion) + call hchksum(fluxes%seaice_melt, mesg//" fluxes%seaice_melt", G%HI, haloshift=hshift, scale=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) if (associated(fluxes%p_surf)) & call hchksum(fluxes%p_surf, mesg//" fluxes%p_surf", G%HI, haloshift=hshift , scale=US%RL2_T2_to_Pa) if (associated(fluxes%salt_flux)) & - call hchksum(fluxes%salt_flux, mesg//" fluxes%salt_flux", G%HI, haloshift=hshift, scale=RZ_T_conversion) + call hchksum(fluxes%salt_flux, mesg//" fluxes%salt_flux", G%HI, haloshift=hshift, scale=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) 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, scale=US%Z_to_m*US%s_to_T) if (associated(fluxes%lrunoff)) & - call hchksum(fluxes%lrunoff, mesg//" fluxes%lrunoff",G%HI,haloshift=hshift, scale=RZ_T_conversion) + call hchksum(fluxes%lrunoff, mesg//" fluxes%lrunoff", G%HI, haloshift=hshift, scale=US%RZ_T_to_kg_m2s) if (associated(fluxes%frunoff)) & - call hchksum(fluxes%frunoff, mesg//" fluxes%frunoff",G%HI,haloshift=hshift, scale=RZ_T_conversion) + call hchksum(fluxes%frunoff, mesg//" fluxes%frunoff", G%HI, haloshift=hshift, scale=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=RZ_T_conversion) + haloshift=hshift, scale=US%RZ_T_to_kg_m2s) 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) @@ -2065,7 +2063,6 @@ subroutine copy_common_forcing_fields(forces, fluxes, G, skip_pres) type(ocean_grid_type), intent(in) :: G !< grid type logical, optional, intent(in) :: skip_pres !< If present and true, do not copy pressure fields. - real :: taux2, tauy2 ! Squared wind stress components [Pa2]. logical :: do_pres integer :: i, j, is, ie, js, je is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec @@ -2203,7 +2200,6 @@ subroutine copy_back_forcing_fields(fluxes, forces, G) type(mech_forcing), intent(inout) :: forces !< A structure with the driving mechanical forces type(ocean_grid_type), intent(in) :: G !< grid type - real :: taux2, tauy2 ! Squared wind stress components [Pa2]. integer :: i, j, is, ie, js, je is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec @@ -2290,9 +2286,8 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h real, dimension(SZI_(diag%G),SZJ_(diag%G)) :: res real :: total_transport ! for diagnosing integrated boundary transport real :: ave_flux ! for diagnosing averaged boundary flux - real :: C_p ! seawater heat capacity [J degC-1 kg-1] real :: RZ_T_conversion ! A combination of scaling factors for mass fluxes [kg T m-2 s-1 R-1 Z-1 ~> 1] - real :: I_dt ! inverse time step [s-1] + real :: I_dt ! inverse time step [T-1 ~> s-1] real :: ppt2mks ! conversion between ppt and mks integer :: turns ! Number of index quarter turns integer :: i,j,is,ie,js,je @@ -2312,9 +2307,8 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h fluxes => fluxes_in endif - C_p = US%Q_to_J_kg*fluxes%C_p RZ_T_conversion = US%RZ_T_to_kg_m2s - I_dt = 1.0 / (US%T_to_s*fluxes%dt_buoy_accum) + I_dt = 1.0 / fluxes%dt_buoy_accum ppt2mks = 1e-3 is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec @@ -2341,7 +2335,7 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h call post_data(handles%id_total_prcme, total_transport, diag) endif if (handles%id_prcme_ga > 0) then - ave_flux = global_area_mean(res,G) + ave_flux = global_area_mean(res, G) call post_data(handles%id_prcme_ga, ave_flux, diag) endif endif @@ -2365,7 +2359,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) + total_transport = global_area_integral(res, G) call post_data(handles%id_total_net_massout, total_transport, diag) endif endif @@ -2401,7 +2395,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) + total_transport = global_area_integral(res, G) call post_data(handles%id_total_net_massin, total_transport, diag) endif endif @@ -2412,11 +2406,11 @@ 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=RZ_T_conversion) + total_transport = global_area_integral(fluxes%evap, G, scale=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 - ave_flux = global_area_mean(fluxes%evap, G, scale=RZ_T_conversion) + ave_flux = global_area_mean(fluxes%evap, G, scale=US%RZ_T_to_kg_m2s) call post_data(handles%id_evap_ga, ave_flux, diag) endif @@ -2426,11 +2420,11 @@ 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) + total_transport = global_area_integral(res, G) call post_data(handles%id_total_precip, total_transport, diag) endif if (handles%id_precip_ga > 0) then - ave_flux = global_area_mean(res,G) + ave_flux = global_area_mean(res, G) call post_data(handles%id_precip_ga, ave_flux, diag) endif endif @@ -2438,11 +2432,11 @@ 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=RZ_T_conversion) + total_transport = global_area_integral(fluxes%lprec, G, scale=US%RZ_T_to_kg_m2s) call post_data(handles%id_total_lprec, total_transport, diag) endif if (handles%id_lprec_ga > 0) then - ave_flux = global_area_mean(fluxes%lprec, G, scale=RZ_T_conversion) + ave_flux = global_area_mean(fluxes%lprec, G, scale=US%RZ_T_to_kg_m2s) call post_data(handles%id_lprec_ga, ave_flux, diag) endif endif @@ -2450,11 +2444,11 @@ 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=RZ_T_conversion) + total_transport = global_area_integral(fluxes%fprec, G, scale=US%RZ_T_to_kg_m2s) call post_data(handles%id_total_fprec, total_transport, diag) endif if (handles%id_fprec_ga > 0) then - ave_flux = global_area_mean(fluxes%fprec, G, scale=RZ_T_conversion) + ave_flux = global_area_mean(fluxes%fprec, G, scale=US%RZ_T_to_kg_m2s) call post_data(handles%id_fprec_ga, ave_flux, diag) endif endif @@ -2462,11 +2456,11 @@ 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=RZ_T_conversion) + total_transport = global_area_integral(fluxes%vprec, G, scale=US%RZ_T_to_kg_m2s) call post_data(handles%id_total_vprec, total_transport, diag) endif if (handles%id_vprec_ga > 0) then - ave_flux = global_area_mean(fluxes%vprec, G, scale=RZ_T_conversion) + ave_flux = global_area_mean(fluxes%vprec, G, scale=US%RZ_T_to_kg_m2s) call post_data(handles%id_vprec_ga, ave_flux, diag) endif endif @@ -2474,7 +2468,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=RZ_T_conversion) + total_transport = global_area_integral(fluxes%lrunoff, G, scale=US%RZ_T_to_kg_m2s) call post_data(handles%id_total_lrunoff, total_transport, diag) endif endif @@ -2482,7 +2476,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=RZ_T_conversion) + total_transport = global_area_integral(fluxes%frunoff, G, scale=US%RZ_T_to_kg_m2s) call post_data(handles%id_total_frunoff, total_transport, diag) endif endif @@ -2490,7 +2484,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=RZ_T_conversion) + total_transport = global_area_integral(fluxes%seaice_melt, G, scale=US%RZ_T_to_kg_m2s) call post_data(handles%id_total_seaice_melt, total_transport, diag) endif endif @@ -2549,7 +2543,7 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h 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, scale=US%QRZ_T_to_W_m2) call post_data(handles%id_total_heat_content_massout, total_transport, diag) endif @@ -2590,9 +2584,9 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h if (associated(fluxes%sens)) res(i,j) = res(i,j) + fluxes%sens(i,j) if (associated(fluxes%SW)) res(i,j) = res(i,j) + fluxes%sw(i,j) if (associated(fluxes%seaice_melt_heat)) res(i,j) = res(i,j) + fluxes%seaice_melt_heat(i,j) - if (allocated(sfc_state%frazil)) res(i,j) = res(i,j) + US%W_m2_to_QRZ_T*sfc_state%frazil(i,j) * I_dt + if (allocated(sfc_state%frazil)) res(i,j) = res(i,j) + sfc_state%frazil(i,j) * I_dt !if (associated(sfc_state%TempXpme)) then - ! res(i,j) = res(i,j) + sfc_state%TempXpme(i,j) * US%Q_to_J_kg*fluxes%C_p * I_dt + ! res(i,j) = res(i,j) + sfc_state%TempXpme(i,j) * fluxes%C_p * I_dt !else if (associated(fluxes%heat_content_lrunoff)) & res(i,j) = res(i,j) + fluxes%heat_content_lrunoff(i,j) @@ -2629,7 +2623,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) = 0.0 ! if (associated(sfc_state%TempXpme)) then - ! res(i,j) = res(i,j) + sfc_state%TempXpme(i,j) * US%Q_to_J_kg*fluxes%C_p * I_dt + ! res(i,j) = res(i,j) + sfc_state%TempXpme(i,j) * fluxes%C_p * I_dt ! else if (associated(fluxes%heat_content_lrunoff)) res(i,j) = res(i,j) + fluxes%heat_content_lrunoff(i,j) if (associated(fluxes%heat_content_frunoff)) res(i,j) = res(i,j) + fluxes%heat_content_frunoff(i,j) @@ -2800,21 +2794,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=RZ_T_conversion) + total_transport = ppt2mks*global_area_integral(fluxes%salt_flux, G, scale=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=RZ_T_conversion) + total_transport = ppt2mks*global_area_integral(fluxes%salt_flux_added, G, scale=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=RZ_T_conversion) + total_transport = ppt2mks*global_area_integral(fluxes%salt_flux_in, G, scale=US%RZ_T_to_kg_m2s) call post_data(handles%id_total_saltFluxIn, total_transport, diag) endif diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index a8a367a7fe..97e5b36db5 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -46,23 +46,23 @@ module MOM_variables Hml, & !< The mixed layer depth [Z ~> m]. u, & !< The mixed layer zonal velocity [L T-1 ~> m s-1]. v, & !< The mixed layer meridional velocity [L T-1 ~> m s-1]. - sea_lev, & !< The sea level [m]. If a reduced surface gravity is + sea_lev, & !< The sea level [Z ~> m]. If a reduced surface gravity is !! used, that is compensated for in sea_lev. frazil, & !< The energy needed to heat the ocean column to the freezing point during - !! the call to step_MOM [J m-2]. - melt_potential, & !< Instantaneous amount of heat that can be used to melt sea ice [J m-2]. + !! the call to step_MOM [Q R Z ~> J m-2]. + melt_potential, & !< Instantaneous amount of heat that can be used to melt sea ice [Q R Z ~> J m-2]. !! This is computed w.r.t. surface freezing temperature. - ocean_mass, & !< The total mass of the ocean [kg m-2]. - ocean_heat, & !< The total heat content of the ocean in [degC kg m-2]. - ocean_salt, & !< The total salt content of the ocean in [kgSalt m-2]. - taux_shelf, & !< The zonal stresses on the ocean under shelves [Pa]. - tauy_shelf, & !< The meridional stresses on the ocean under shelves [Pa]. + ocean_mass, & !< The total mass of the ocean [R Z ~> kg m-2]. + ocean_heat, & !< The total heat content of the ocean in [degC R Z ~> degC kg m-2]. + ocean_salt, & !< The total salt content of the ocean in [kgSalt kg-1 R Z ~> kgSalt m-2]. + taux_shelf, & !< The zonal stresses on the ocean under shelves [R L Z T-2 ~> Pa]. + tauy_shelf, & !< The meridional stresses on the ocean under shelves [R L Z T-2 ~> Pa]. TempxPmE, & !< The net inflow of water into the ocean times the temperature at which this - !! inflow occurs during the call to step_MOM [degC kg m-2]. - salt_deficit, & !< The salt needed to maintain the ocean column at a minimum - !! salinity of 0.01 PSU over the call to step_MOM [kgSalt m-2]. + !! inflow occurs during the call to step_MOM [degC R Z ~> degC kg m-2]. + salt_deficit, & !< The salt needed to maintain the ocean column above a minimum + !! salinity over the call to step_MOM [kgSalt kg-1 R Z ~> kgSalt m-2]. internal_heat !< Any internal or geothermal heat sources that are applied to the ocean - !! integrated over the call to step_MOM [degC kg m-2]. + !! integrated over the call to step_MOM [degC R Z ~> degC kg m-2]. logical :: T_is_conT = .false. !< If true, the temperature variable SST is actually the !! conservative temperature in [degC]. logical :: S_is_absS = .false. !< If true, the salinity variable SSS is actually the diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index d104dfe82a..891d6b3ea7 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -273,13 +273,12 @@ subroutine shelf_calc_flux(sfc_state, fluxes, Time, time_step, CS, forces) real :: wB_flux_new, dDwB_dwB_in real :: I_Gam_T, I_Gam_S real :: dG_dwB ! The derivative of Gam_turb with wB [T3 Z-2 ~> s3 m-2] - real :: taux2, tauy2 ! The squared surface stresses [Pa]. + real :: taux2, tauy2 ! The squared surface stresses [R2 L2 Z2 T-4 ~> Pa2]. real :: u2_av, v2_av ! The ice-area weighted average squared ocean velocities [L2 T-2 ~> m2 s-2] real :: asu1, asu2 ! Ocean areas covered by ice shelves at neighboring u- real :: asv1, asv2 ! and v-points [L2 ~> m2]. real :: I_au, I_av ! The Adcroft reciprocals of the ice shelf areas at adjacent points [L-2 ~> m-2] - real :: Irho0 ! The inverse of the mean density times unit conversion factors that - ! arise because state uses MKS units [L2 m s2 kg-1 T-2 ~> m3 kg-1]. + real :: Irho0 ! The inverse of the mean density times a unit conversion factor [R-1 L Z-1 ~> m3 kg-1] logical :: Sb_min_set, Sb_max_set logical :: update_ice_vel ! If true, it is time to update the ice shelf velocities. logical :: coupled_GL ! If true, the grouding line position is determined based on @@ -336,14 +335,15 @@ subroutine shelf_calc_flux(sfc_state, fluxes, Time, time_step, CS, forces) call hchksum(sfc_state%sss, "sss before apply melting", G%HI, haloshift=0) call hchksum(sfc_state%u, "u_ml before apply melting", G%HI, haloshift=0, scale=US%L_T_to_m_s) call hchksum(sfc_state%v, "v_ml before apply melting", G%HI, haloshift=0, scale=US%L_T_to_m_s) - call hchksum(sfc_state%ocean_mass, "ocean_mass before apply melting", G%HI, haloshift=0) + call hchksum(sfc_state%ocean_mass, "ocean_mass before apply melting", G%HI, haloshift=0, & + scale=US%RZ_to_kg_m2) endif ! Calculate the friction velocity under ice shelves, using taux_shelf and tauy_shelf if possible. if (allocated(sfc_state%taux_shelf) .and. allocated(sfc_state%tauy_shelf)) then call pass_vector(sfc_state%taux_shelf, sfc_state%tauy_shelf, G%domain, TO_ALL, CGRID_NE) endif - Irho0 = US%m_s_to_L_T**2*US%kg_m3_to_R / CS%Rho_ocn + Irho0 = US%Z_to_L / CS%Rho_ocn do j=js,je ; do i=is,ie ; if (fluxes%frac_shelf_h(i,j) > 0.0) then taux2 = 0.0 ; tauy2 = 0.0 ; u2_av = 0.0 ; v2_av = 0.0 asu1 = (ISS%area_shelf_h(i-1,j) + ISS%area_shelf_h(i,j)) @@ -383,7 +383,7 @@ subroutine shelf_calc_flux(sfc_state, fluxes, Time, time_step, CS, forces) CS%eqn_of_state, EOSdom) do i=is,ie - if ((sfc_state%ocean_mass(i,j) > US%RZ_to_kg_m2*CS%col_mass_melt_threshold) .and. & + if ((sfc_state%ocean_mass(i,j) > CS%col_mass_melt_threshold) .and. & (ISS%area_shelf_h(i,j) > 0.0) .and. CS%isthermo) then if (CS%threeeq) then @@ -616,7 +616,7 @@ subroutine shelf_calc_flux(sfc_state, fluxes, Time, time_step, CS, forces) fluxes%iceshelf_melt(:,:) = ISS%water_flux(:,:) * CS%flux_factor do j=js,je ; do i=is,ie - if ((sfc_state%ocean_mass(i,j) > US%RZ_to_kg_m2*CS%col_mass_melt_threshold) .and. & + if ((sfc_state%ocean_mass(i,j) > CS%col_mass_melt_threshold) .and. & (ISS%area_shelf_h(i,j) > 0.0) .and. (CS%isthermo)) then ! Set melt to zero above a cutoff pressure (CS%Rho_ocn*CS%cutoff_depth*CS%g_Earth). @@ -655,7 +655,7 @@ subroutine shelf_calc_flux(sfc_state, fluxes, Time, time_step, CS, forces) endif ! area_shelf_h enddo ; enddo ! i- and j-loops - ! mass flux [kg s-1], part of ISOMIP diags. + ! mass flux [R Z L2 T-1 ~> kg s-1], part of ISOMIP diags. mass_flux(:,:) = 0.0 mass_flux(:,:) = ISS%water_flux(:,:) * ISS%area_shelf_h(:,:) @@ -690,7 +690,7 @@ subroutine shelf_calc_flux(sfc_state, fluxes, Time, time_step, CS, forces) ! advect the ice shelf, and advance the front. Calving will be in here somewhere as well.. ! when we decide on how to do it call update_ice_shelf(CS%dCS, ISS, G, US, US%s_to_T*time_step, Time, & - US%kg_m3_to_R*US%m_to_Z*sfc_state%ocean_mass(:,:), coupled_GL) + sfc_state%ocean_mass(:,:), coupled_GL) endif @@ -890,10 +890,10 @@ subroutine add_shelf_flux(G, US, CS, sfc_state, fluxes) type(forcing), intent(inout) :: fluxes !< A structure of surface fluxes that may be used/updated. ! local variables - real :: frac_shelf !< The fractional area covered by the ice shelf [nondim]. - real :: frac_open !< The fractional area of the ocean that is not covered by the ice shelf [nondim]. - real :: delta_mass_shelf!< Change in ice shelf mass over one time step [kg s-1] - real :: balancing_flux !< The fresh water flux that balances the integrated melt flux [R Z T-1 ~> kg m-2 s-1] + real :: frac_shelf !< The fractional area covered by the ice shelf [nondim]. + real :: frac_open !< The fractional area of the ocean that is not covered by the ice shelf [nondim]. + real :: delta_mass_shelf !< Change in ice shelf mass over one time step [R Z m2 T-1 ~> kg s-1] + real :: balancing_flux !< The fresh water flux that balances the integrated melt flux [R Z T-1 ~> kg m-2 s-1] real :: balancing_area !< total area where the balancing flux is applied [m2] type(time_type) :: dTime !< The time step as a time_type type(time_type) :: Time0 !< The previous time (Time-dt) @@ -905,7 +905,7 @@ subroutine add_shelf_flux(G, US, CS, sfc_state, fluxes) !! the two timesteps at (Time) and (Time-dt) [R Z ~> kg m-2]. real, dimension(SZDI_(G),SZDJ_(G)) :: last_h_shelf !< Ice shelf thickness [Z ~> m] !! at at previous time (Time-dt) - real, dimension(SZDI_(G),SZDJ_(G)) :: last_hmask !< Ice shelf mask + real, dimension(SZDI_(G),SZDJ_(G)) :: last_hmask !< Ice shelf mask [nondim] !! at at previous time (Time-dt) real, dimension(SZDI_(G),SZDJ_(G)) :: last_area_shelf_h !< Ice shelf area [L2 ~> m2] !! at at previous time (Time-dt) @@ -933,7 +933,7 @@ subroutine add_shelf_flux(G, US, CS, sfc_state, fluxes) if (CS%debug) then if (allocated(sfc_state%taux_shelf) .and. allocated(sfc_state%tauy_shelf)) then call uvchksum("tau[xy]_shelf", sfc_state%taux_shelf, sfc_state%tauy_shelf, & - G%HI, haloshift=0) + G%HI, haloshift=0, scale=US%RZ_T_to_kg_m2s*US%L_T_to_m_s) endif endif @@ -1023,7 +1023,7 @@ subroutine add_shelf_flux(G, US, CS, sfc_state, fluxes) ! get total ice shelf mass at (Time-dt) and (Time), in kg do j=js,je ; do i=is,ie ! Just consider the change in the mass of the floating shelf. - if ((sfc_state%ocean_mass(i,j) > US%RZ_to_kg_m2*CS%min_ocean_mass_float) .and. & + if ((sfc_state%ocean_mass(i,j) > CS%min_ocean_mass_float) .and. & (ISS%area_shelf_h(i,j) > 0.0)) then delta_float_mass(i,j) = ISS%mass_shelf(i,j) - last_mass_shelf(i,j) else @@ -1051,7 +1051,7 @@ subroutine add_shelf_flux(G, US, CS, sfc_state, fluxes) enddo ; enddo balancing_area = global_area_integral(bal_frac, G) - if (balancing_area > 0.0) then + if (balancing_area > 0.0) then !### Examine whether the rescaling should be inside the parenthesis. balancing_flux = US%kg_m2s_to_RZ_T*(global_area_integral(ISS%water_flux, G, scale=US%RZ_T_to_kg_m2s, & area=ISS%area_shelf_h) + & delta_mass_shelf ) / balancing_area @@ -1123,7 +1123,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl character(len=240) :: Tideamp_file real :: utide ! A tidal velocity [L T-1 ~> m s-1] real :: col_thick_melt_thresh ! An ocean column thickness below which iceshelf melting - ! does not occur [m] + ! does not occur [Z ~> m] if (associated(CS)) then call MOM_error(FATAL, "MOM_ice_shelf.F90, initialize_ice_shelf: "// & "called with an associated control structure.") diff --git a/src/ice_shelf/MOM_marine_ice.F90 b/src/ice_shelf/MOM_marine_ice.F90 index 30121d0c8e..64d4dbfdab 100644 --- a/src/ice_shelf/MOM_marine_ice.F90 +++ b/src/ice_shelf/MOM_marine_ice.F90 @@ -48,8 +48,8 @@ subroutine iceberg_forces(G, forces, use_ice_shelf, sfc_state, time_step, CS) type(surface), intent(inout) :: sfc_state !< A structure containing fields that !! describe the surface state of the ocean. logical, intent(in) :: use_ice_shelf !< If true, this configuration uses ice shelves. - real, intent(in) :: time_step !< The coupling time step [s]. - type(marine_ice_CS), pointer :: CS !< Pointer to the control structure for MOM_marine_ice + real, intent(in) :: time_step !< The coupling time step [s]. + type(marine_ice_CS), pointer :: CS !< Pointer to the control structure for MOM_marine_ice real :: kv_rho_ice ! The viscosity of ice divided by its density [L4 Z-2 T-1 R-1 ~> m5 kg-1 s-1]. integer :: i, j, is, ie, js, je @@ -107,18 +107,17 @@ subroutine iceberg_fluxes(G, US, fluxes, use_ice_shelf, sfc_state, time_step, CS !! describe the surface state of the ocean. logical, intent(in) :: use_ice_shelf !< If true, this configuration uses ice shelves. real, intent(in) :: time_step !< The coupling time step [s]. - type(marine_ice_CS), pointer :: CS !< Pointer to the control structure for MOM_marine_ice + type(marine_ice_CS), pointer :: CS !< Pointer to the control structure for MOM_marine_ice real :: fraz ! refreezing rate [R Z T-1 ~> kg m-2 s-1] - real :: I_dt_LHF ! The inverse of the timestep times the latent heat of fusion times unit conversion - ! factors because sfc_state is in MKS units [R Z m2 J-1 T-1 ~> kg J-1 s-1]. + real :: I_dt_LHF ! The inverse of the timestep times the latent heat of fusion times [Q-1 T-1 ~> kg J-1 s-1]. integer :: i, j, is, ie, js, je, isd, ied, jsd, jed is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec isd = G%isd ; jsd = G%jsd ; ied = G%ied ; jed = G%jed - !This routine adds iceberg data to the ice shelf data (if ice shelf is used) - !which can then be used to change the top of ocean boundary condition used in - !the ocean model. This routine is taken from the add_shelf_flux subroutine - !within the ice shelf model. + ! This routine adds iceberg data to the ice shelf data (if ice shelf is used) + ! which can then be used to change the top of ocean boundary condition used in + ! the ocean model. This routine is taken from the add_shelf_flux subroutine + ! within the ice shelf model. if (.not.associated(CS)) return if (.not.(associated(fluxes%area_berg) .and. associated(fluxes%ustar_berg) .and. & @@ -139,7 +138,7 @@ subroutine iceberg_fluxes(G, US, fluxes, use_ice_shelf, sfc_state, time_step, CS !Zero'ing out other fluxes under the tabular icebergs if (CS%berg_area_threshold >= 0.) then - I_dt_LHF = US%W_m2_to_QRZ_T / (time_step * CS%latent_heat_fusion) + I_dt_LHF = 1.0 / (US%s_to_T*time_step * CS%latent_heat_fusion) do j=jsd,jed ; do i=isd,ied if (fluxes%frac_shelf_h(i,j) > CS%berg_area_threshold) then ! Only applying for ice shelf covering most of cell. @@ -149,7 +148,7 @@ subroutine iceberg_fluxes(G, US, fluxes, use_ice_shelf, sfc_state, time_step, CS if (associated(fluxes%latent)) fluxes%latent(i,j) = 0.0 if (associated(fluxes%evap)) fluxes%evap(i,j) = 0.0 - ! Add frazil formation diagnosed by the ocean model [J m-2] in the + ! Add frazil formation diagnosed by the ocean model [Q R Z ~> J m-2] in the ! form of surface layer evaporation [R Z T-1 ~> kg m-2 s-1]. Update lprec in the ! control structure for diagnostic purposes. From 3a04ca549a982f3dcb504a82fd41ce5744c5725a Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Mon, 4 May 2020 16:38:56 -0400 Subject: [PATCH 11/11] Scalar diag conversion; Q and R dimension scaling This patch adds support for enthalpy (Q) and density (R) dimensional scaling into the test suite. It also resolves an issue with conversion scaling in scalars (post_data_0d) which were not being applied, and make it impossible to verify the dimensions of scalar diagnostics. --- .testing/Makefile | 14 +++++++++----- src/framework/MOM_diag_mediator.F90 | 12 +++++++++--- 2 files changed, 18 insertions(+), 8 deletions(-) diff --git a/.testing/Makefile b/.testing/Makefile index 0d73979204..99672268c3 100644 --- a/.testing/Makefile +++ b/.testing/Makefile @@ -1,4 +1,5 @@ SHELL = bash +.SUFFIXES: # User-defined configuration -include config.mk @@ -30,9 +31,9 @@ MKMF_CPP = "-Duse_libMPI -Duse_netCDF -DSPMD" # Environment # TODO: This info ought to be determined by CMake, automake, etc. -#MKMF_TEMPLATE ?= .testing/linux-ubuntu-xenial-gnu.mk -MKMF_TEMPLATE ?= build/mkmf/templates/ncrc-gnu.mk -#MKMF_TEMPLATE ?= build/mkmf/templates/ncrc-intel.mk +#MKMF_TEMPLATE ?= linux-ubuntu-xenial-gnu.mk +MKMF_TEMPLATE ?= deps/mkmf/templates/ncrc-gnu.mk +#MKMF_TEMPLATE ?= deps/mkmf/templates/ncrc-intel.mk #--- # Test configuration @@ -41,6 +42,7 @@ MKMF_TEMPLATE ?= build/mkmf/templates/ncrc-gnu.mk BUILDS = symmetric asymmetric repro openmp CONFIGS := $(wildcard tc*) TESTS = grids layouts restarts nans dims openmps rotations +DIMS = t l h z q r # REPRO tests enable reproducibility with optimization, and often do not match # the DEBUG results in older GCCs and vendor compilers, so we can optionally @@ -199,7 +201,7 @@ test.restarts: $(foreach c,$(CONFIGS),$(c).restart) test.repros: $(foreach c,$(CONFIGS),$(c).repro $(c).repro.diag) test.openmps: $(foreach c,$(CONFIGS),$(c).openmp $(c).openmp.diag) test.nans: $(foreach c,$(CONFIGS),$(c).nan $(c).nan.diag) -test.dims: $(foreach c,$(CONFIGS),$(foreach d,t l h z,$(c).dim.$(d) $(c).dim.$(d).diag)) +test.dims: $(foreach c,$(CONFIGS),$(foreach d,$(DIMS),$(c).dim.$(d) $(c).dim.$(d).diag)) test.regressions: $(foreach c,$(CONFIGS),$(c).regression $(c).regression.diag) ! ls -1 results/*/*.reg @@ -220,7 +222,7 @@ $(eval $(call CMP_RULE,rotate,symmetric rotate)) $(eval $(call CMP_RULE,repro,symmetric repro)) $(eval $(call CMP_RULE,openmp,symmetric openmp)) $(eval $(call CMP_RULE,nan,symmetric nan)) -$(foreach d,t l h z,$(eval $(call CMP_RULE,dim.$(d),symmetric dim.$(d)))) +$(foreach d,$(DIMS),$(eval $(call CMP_RULE,dim.$(d),symmetric dim.$(d)))) # Custom comparison rules @@ -295,6 +297,8 @@ $(eval $(call STAT_RULE,dim.t,symmetric,,T_RESCALE_POWER=11,,1)) $(eval $(call STAT_RULE,dim.l,symmetric,,L_RESCALE_POWER=11,,1)) $(eval $(call STAT_RULE,dim.h,symmetric,,H_RESCALE_POWER=11,,1)) $(eval $(call STAT_RULE,dim.z,symmetric,,Z_RESCALE_POWER=11,,1)) +$(eval $(call STAT_RULE,dim.q,symmetric,,Q_RESCALE_POWER=11,,1)) +$(eval $(call STAT_RULE,dim.r,symmetric,,R_RESCALE_POWER=11,,1)) # Restart tests require significant preprocessing, and are handled separately. results/%/ocean.stats.restart: build/symmetric/MOM6 diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 0fcee7624a..03de6405fe 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -1212,6 +1212,7 @@ subroutine post_data_0d(diag_field_id, field, diag_cs, is_static) logical, optional, intent(in) :: is_static !< If true, this is a static field that is always offered. ! Local variables + real :: locfield logical :: used, is_stat type(diag_type), pointer :: diag => null() @@ -1223,13 +1224,18 @@ subroutine post_data_0d(diag_field_id, field, diag_cs, is_static) call assert(diag_field_id < diag_cs%next_free_diag_id, & 'post_data_0d: Unregistered diagnostic id') diag => diag_cs%diags(diag_field_id) + do while (associated(diag)) + locfield = field + if (diag%conversion_factor /= 0.) & + locfield = locfield * diag%conversion_factor + if (diag_cs%diag_as_chksum) then - call chksum0(field, diag%debug_str, logunit=diag_cs%chksum_iounit) + call chksum0(locfield, diag%debug_str, logunit=diag_cs%chksum_iounit) else if (is_stat) then - used = send_data(diag%fms_diag_id, field) + used = send_data(diag%fms_diag_id, locfield) elseif (diag_cs%ave_enabled) then - used = send_data(diag%fms_diag_id, field, diag_cs%time_end) + used = send_data(diag%fms_diag_id, locfield, diag_cs%time_end) endif diag => diag%next enddo