From 8f7df084dd8262aebe126db6981abb0ad9c042f1 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 24 Jan 2024 06:40:34 -0500 Subject: [PATCH 01/15] (*)non-Boussinesq Z unit update_OBC_segment_data Revised the update_OBC_segment_data code to keep the z-space input data in height units ([Z ~> m]) rather than rescaling them to thickness units. This change means that the non-Boussinesq open boundary condition calculations avoid using the Boussinesq reference density. Associated with this change, 5 internal variables in update_OBC_segment_data were renamed to reflect that they are layer vertical extents instead of thicknesses. Also added or amended comments describing the purpose and units of 25 real variables in this module and corrected a handful of lines in this module that were not adhering to the guidance in the MOM6 style guide. Answers are bitwise identical in any Boussinesq cases. However, answers will change in any non-Boussinesq cases that use open boundary conditions that use time_interp_extrnal to read in Z-space data. --- src/core/MOM_open_boundary.F90 | 250 +++++++++++++++++---------------- 1 file changed, 129 insertions(+), 121 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 76ac477906..94320a30c7 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -91,28 +91,36 @@ module MOM_open_boundary character(len=32) :: name !< a name identifier for the segment data character(len=8) :: genre !< an identifier for the segment data real :: scale !< A scaling factor for converting input data to - !! the internal units of this field - real, allocatable :: buffer_src(:,:,:) !< buffer for segment data located at cell faces - !! and on the original vertical grid. The values for tracers should - !! have the same units as the field they are being applied to? + !! the internal units of this field. For salinity this would + !! be in units of [S ppt-1 ~> 1] + real, allocatable :: buffer_src(:,:,:) !< buffer for segment data located at cell faces and on + !! the original vertical grid in the internally scaled + !! units for the field in question, such as [L T-1 ~> m s-1] + !! for a velocity or [S ~> ppt] for salinity. integer :: nk_src !< Number of vertical levels in the source data real, allocatable :: dz_src(:,:,:) !< vertical grid cell spacing of the incoming segment - !! data, set in [Z ~> m] then scaled to [H ~> m or kg m-2] - real, allocatable :: buffer_dst(:,:,:) !< buffer src data remapped to the target vertical grid. - !! The values for tracers should have the same units as the field - !! they are being applied to? - real :: value !< constant value if not read from file - real :: resrv_lfac_in = 1. !< reservoir inverse length scale factor for IN direction per field - !< the general 1/Lscale_IN is multiplied by this factor for each tracer - real :: resrv_lfac_out= 1. !< reservoir inverse length scale factor for OUT direction per field - !< the general 1/Lscale_OUT is multiplied by this factor for each tracer + !! data in [Z ~> m]. + real, allocatable :: buffer_dst(:,:,:) !< buffer src data remapped to the target vertical grid + !! in the internally scaled units for the field in + !! question, such as [L T-1 ~> m s-1] for a velocity or + !! [S ~> ppt] for salinity. + real :: value !< A constant value for the inflow concentration if not read + !! from file, in the internal units of a field, such as [S ~> ppt] + !! for salinity. + real :: resrv_lfac_in = 1. !< The reservoir inverse length scale factor for the inward + !! direction per field [nondim]. The general 1/Lscale_in is + !! multiplied by this factor for a specific tracer. + real :: resrv_lfac_out= 1. !< The reservoir inverse length scale factor for the outward + !! direction per field [nondim]. The general 1/Lscale_out is + !! multiplied by this factor for a specific tracer. end type OBC_segment_data_type !> Tracer on OBC segment data structure, for putting into a segment tracer registry. type, public :: OBC_segment_tracer_type real, allocatable :: t(:,:,:) !< tracer concentration array in rescaled units, !! like [S ~> ppt] for salinity. - real :: OBC_inflow_conc = 0.0 !< tracer concentration for generic inflows + real :: OBC_inflow_conc = 0.0 !< tracer concentration for generic inflows in rescaled units, + !! like [S ~> ppt] for salinity. character(len=32) :: name !< tracer name used for error messages type(tracer_type), pointer :: Tr => NULL() !< metadata describing the tracer real, allocatable :: tres(:,:,:) !< tracer reservoir array in rescaled units, @@ -380,7 +388,7 @@ module MOM_open_boundary !> Control structure for open boundaries that read from files. !! Probably lots to update here. type, public :: file_OBC_CS ; private - real :: tide_flow = 3.0e6 !< Placeholder for now... + real :: tide_flow = 3.0e6 !< Placeholder for now..., perhaps in [m3 s-1]? end type file_OBC_CS !> Type to carry something (what??) for the OBC registry. @@ -402,8 +410,8 @@ module MOM_open_boundary character(len=128) :: tracer_name !< tracer name character(len=128) :: tracer_src_file !< tracer source file for BC character(len=128) :: tracer_src_field !< name of the field in source file to extract BC - real :: lfac_in !< multiplicative factor for inbound tracer reservoir length scale - real :: lfac_out !< multiplicative factor for outbound tracer reservoir length scale + real :: lfac_in !< multiplicative factor for inbound tracer reservoir length scale [nondim] + real :: lfac_out !< multiplicative factor for outbound tracer reservoir length scale [nondim] end type external_tracers_segments_props integer :: id_clock_pass !< A CPU time clock @@ -811,7 +819,7 @@ subroutine initialize_segment_data(G, GV, US, OBC, PF) obgc_segments_props_list => OBC%obgc_segments_props !pointer to the head node do m=1,segment%num_fields - if (m .le. num_fields) then + if (m <= num_fields) then !These are tracers with segments specified in MOM6 style override files call parse_segment_data_str(trim(segstr), m, trim(fields(m)), value, filename, fieldname) else @@ -824,8 +832,8 @@ subroutine initialize_segment_data(G, GV, US, OBC, PF) segment%field(m)%resrv_lfac_in,segment%field(m)%resrv_lfac_out) !Make sure the obgc tracer is not specified in the MOM6 param file too. do mm=1,num_fields - if(trim(fields(m)) == trim(fields(mm))) then - if(is_root_pe()) & + if (trim(fields(m)) == trim(fields(mm))) then + if (is_root_pe()) & call MOM_error(FATAL,"MOM_open_boundary:initialize_segment_data(): obgc tracer " //trim(fields(m))// & " appears in OBC_SEGMENT_XXX_DATA string in MOM6 param file. This is not supported!") endif @@ -858,24 +866,24 @@ subroutine initialize_segment_data(G, GV, US, OBC, PF) call MOM_error(WARNING, mesg // " " // trim(filename) // " " // trim(fieldname)) call MOM_error(FATAL,'segment data are not on the supergrid') endif - siz2(1)=1 + siz2(1) = 1 if (siz(1)>1) then if (OBC%brushcutter_mode) then - siz2(1)=(siz(1)-1)/2 + siz2(1) = (siz(1)-1)/2 else - siz2(1)=siz(1) + siz2(1) = siz(1) endif endif - siz2(2)=1 + siz2(2) = 1 if (siz(2)>1) then if (OBC%brushcutter_mode) then - siz2(2)=(siz(2)-1)/2 + siz2(2) = (siz(2)-1)/2 else - siz2(2)=siz(2) + siz2(2) = siz(2) endif endif - siz2(3)=siz(3) + siz2(3) = siz(3) if (segment%is_E_or_W) then if (segment%field(m)%name == 'V') then @@ -986,7 +994,7 @@ subroutine initialize_segment_data(G, GV, US, OBC, PF) allocate(segment%field(m)%dz_src(isd:ied,JsdB:JedB,siz(3))) endif endif - segment%field(m)%dz_src(:,:,:)=0.0 + segment%field(m)%dz_src(:,:,:) = 0.0 segment%field(m)%nk_src=siz(3) segment%field(m)%dz_handle = init_external_field(trim(filename), trim(fieldname), & ignore_axis_atts=.true., threading=SINGLE_FILE) @@ -1746,7 +1754,8 @@ subroutine parse_segment_data_str(segment_str, idx, var, value, filename, fieldn character(len=*), intent(out) :: filename !< The name of the input file if using "file" method character(len=*), intent(out) :: fieldname !< The name of the variable in the input file if using !! "file" method - real, optional, intent(out) :: value !< A constant value if using the "value" method + real, optional, intent(out) :: value !< A constant value if using the "value" method in various + !! units but without the internal rescaling [various units] ! Local variables character(len=128) :: word1, word2, word3, method @@ -1814,8 +1823,7 @@ subroutine parse_for_tracer_reservoirs(OBC, PF, use_temperature) ! At this point, just search for TEMP and SALT as tracers 1 and 2. do m=1,num_fields - call parse_segment_data_str(trim(segstr), m, trim(fields(m)), & - value, filename, fieldname) + call parse_segment_data_str(trim(segstr), m, trim(fields(m)), value, filename, fieldname) if (trim(filename) /= 'none') then if (fields(m) == 'TEMP') then if (segment%is_E_or_W_2) then @@ -1877,8 +1885,6 @@ subroutine open_boundary_init(G, GV, US, param_file, OBC, restart_CS) type(MOM_restart_CS), intent(in) :: restart_CS !< Restart structure, data intent(inout) ! Local variables - real :: vel2_rescale ! A rescaling factor for squared velocities from the representation in - ! a restart file to the internal representation in this run. integer :: i, j, k, isd, ied, jsd, jed, nz, m integer :: IsdB, IedB, JsdB, JedB isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nz = GV%ke @@ -1972,9 +1978,9 @@ end subroutine open_boundary_end !> Sets the slope of bathymetry normal to an open boundary to zero. subroutine open_boundary_impose_normal_slope(OBC, G, depth) - type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure - type(dyn_horgrid_type), intent(in) :: G !< Ocean grid structure - real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: depth !< Bathymetry at h-points + type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure + type(dyn_horgrid_type), intent(in) :: G !< Ocean grid structure + real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: depth !< Bathymetry at h-points, in [Z ~> m] or other units ! Local variables integer :: i, j, n type(OBC_segment_type), pointer :: segment => NULL() @@ -3367,11 +3373,11 @@ end subroutine open_boundary_apply_normal_flow !> Applies zero values to 3d u,v fields on OBC segments subroutine open_boundary_zero_normal_flow(OBC, G, GV, u, v) ! Arguments - type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure - type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure - type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure - real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: u !< u field to update on open boundaries - real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: v !< v field to update on open boundaries + type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure + type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: u !< u field to update on open boundaries [arbitrary] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: v !< v field to update on open boundaries [arbitrary] ! Local variables integer :: i, j, k, n type(OBC_segment_type), pointer :: segment => NULL() @@ -3528,7 +3534,7 @@ subroutine set_tracer_data(OBC, tv, h, G, GV, PF) type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure type(ocean_OBC_type), target, intent(in) :: OBC !< Open boundary structure type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamics structure - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Thickness + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2] type(param_file_type), intent(in) :: PF !< Parameter file handle type(OBC_segment_type), pointer :: segment => NULL() ! pointer to segment type list @@ -3791,7 +3797,7 @@ subroutine open_boundary_test_extern_h(G, GV, OBC, h) if (.not. associated(OBC)) return - silly_h = GV%Z_to_H * OBC%silly_h + silly_h = GV%Z_to_H * OBC%silly_h ! This rescaling is here because GV was initialized after OBC. do n = 1, OBC%number_of_segments do k = 1, GV%ke @@ -3844,18 +3850,18 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) integer :: ishift, jshift ! offsets for staggered locations real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: dz ! Distance between the interfaces around a layer [Z ~> m] real, dimension(:,:,:), allocatable, target :: tmp_buffer ! A buffer for input data [various units] - real, dimension(:), allocatable :: h_stack ! Thicknesses at corner points [H ~> m or kg m-2] + real, dimension(:), allocatable :: dz_stack ! Distance between the interfaces at corner points [Z ~> m] integer :: is_obc2, js_obc2 integer :: i_seg_offset, j_seg_offset - real :: net_H_src ! Total thickness of the incoming flow in the source field [H ~> m or kg m-2] - real :: net_H_int ! Total thickness of the incoming flow in the model [H ~> m or kg m-2] + real :: net_dz_src ! Total vertical extent of the incoming flow in the source field [Z ~> m] + real :: net_dz_int ! Total vertical extent of the incoming flow in the model [Z ~> m] real :: scl_fac ! A scaling factor to compensate for differences in total thicknesses [nondim] real :: tidal_vel ! Interpolated tidal velocity at the OBC points [L T-1 ~> m s-1] real :: tidal_elev ! Interpolated tidal elevation at the OBC points [Z ~> m] real, allocatable :: normal_trans_bt(:,:) ! barotropic transport [H L2 T-1 ~> m3 s-1] integer :: turns ! Number of index quarter turns real :: time_delta ! Time since tidal reference date [T ~> s] - real :: h_neglect, h_neglect_edge ! Small thicknesses [H ~> m or kg m-2] + real :: dz_neglect, dz_neglect_edge ! Small thicknesses [Z ~> m] is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed @@ -3868,12 +3874,12 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) if (OBC%add_tide_constituents) time_delta = US%s_to_T * time_type_to_real(Time - OBC%time_ref) - if (OBC%remap_answer_date >= 20190101) then - h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff - elseif (GV%Boussinesq) then - h_neglect = GV%m_to_H * 1.0e-30 ; h_neglect_edge = GV%m_to_H * 1.0e-10 + if (GV%Boussinesq .and. (OBC%remap_answer_date < 20190101)) then + dz_neglect = US%m_to_Z * 1.0e-30 ; dz_neglect_edge = US%m_to_Z * 1.0e-10 + elseif (GV%semi_Boussinesq .and. (OBC%remap_answer_date < 20190101)) then + dz_neglect = GV%kg_m2_to_H*GV%H_to_Z * 1.0e-30 ; dz_neglect_edge = GV%kg_m2_to_H*GV%H_to_Z * 1.0e-10 else - h_neglect = GV%kg_m2_to_H * 1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H * 1.0e-10 + dz_neglect = GV%dZ_subroundoff ; dz_neglect_edge = GV%dZ_subroundoff endif if (OBC%number_of_segments >= 1) then @@ -3937,16 +3943,16 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) enddo endif - allocate(h_stack(GV%ke), source=0.0) + allocate(dz_stack(GV%ke), source=0.0) do m = 1,segment%num_fields !This field may not require a high frequency OBC segment update and might be allowed !a less frequent update as set by the parameter update_OBC_period_max in MOM.F90. !Cycle if it is not the time to update OBC segment data for this field. if (trim(segment%field(m)%genre) == 'obgc' .and. (.not. OBC%update_OBC_seg_data)) cycle if (segment%field(m)%use_IO) then - siz(1)=size(segment%field(m)%buffer_src,1) - siz(2)=size(segment%field(m)%buffer_src,2) - siz(3)=size(segment%field(m)%buffer_src,3) + siz(1) = size(segment%field(m)%buffer_src,1) + siz(2) = size(segment%field(m)%buffer_src,2) + siz(3) = size(segment%field(m)%buffer_src,3) if (.not.allocated(segment%field(m)%buffer_dst)) then if (siz(3) /= segment%field(m)%nk_src) call MOM_error(FATAL,'nk_src inconsistency') if (segment%field(m)%nk_src > 1) then @@ -3990,7 +3996,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) endif endif endif - segment%field(m)%buffer_dst(:,:,:)=0.0 + segment%field(m)%buffer_dst(:,:,:) = 0.0 endif ! read source data interpolated to the current model time ! NOTE: buffer is sized for vertex points, but may be used for faces @@ -4149,6 +4155,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) endif endif + ! The units of ...%dz_src are no longer changed from [Z ~> m] to [H ~> m or kg m-2] here. call adjustSegmentEtaToFitBathymetry(G,GV,US,segment,m) if (segment%is_E_or_W) then @@ -4160,44 +4167,44 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) do J=max(js_obc,jsd),min(je_obc,jed-1) ! Using the h remapping approach ! Pretty sure we need to check for source/target grid consistency here - segment%field(m)%buffer_dst(I,J,:)=0.0 ! initialize remap destination buffer + segment%field(m)%buffer_dst(I,J,:) = 0.0 ! initialize remap destination buffer if (G%mask2dCu(I,j)>0. .and. G%mask2dCu(I,j+1)>0.) then - h_stack(:) = 0.5*(h(i+ishift,j,:) + h(i+ishift,j+1,:)) + dz_stack(:) = 0.5*(dz(i+ishift,j,:) + dz(i+ishift,j+1,:)) call remapping_core_h(OBC%remap_CS, & - segment%field(m)%nk_src,segment%field(m)%dz_src(I,J,:), & + segment%field(m)%nk_src, segment%field(m)%dz_src(I,J,:), & segment%field(m)%buffer_src(I,J,:), & - GV%ke, h_stack, segment%field(m)%buffer_dst(I,J,:), & - h_neglect, h_neglect_edge) + GV%ke, dz_stack, segment%field(m)%buffer_dst(I,J,:), & + dz_neglect, dz_neglect_edge) elseif (G%mask2dCu(I,j)>0.) then - h_stack(:) = h(i+ishift,j,:) + dz_stack(:) = dz(i+ishift,j,:) call remapping_core_h(OBC%remap_CS, & - segment%field(m)%nk_src,segment%field(m)%dz_src(I,J,:), & + segment%field(m)%nk_src, segment%field(m)%dz_src(I,J,:), & segment%field(m)%buffer_src(I,J,:), & - GV%ke, h_stack, segment%field(m)%buffer_dst(I,J,:), & - h_neglect, h_neglect_edge) + GV%ke, dz_stack, segment%field(m)%buffer_dst(I,J,:), & + dz_neglect, dz_neglect_edge) elseif (G%mask2dCu(I,j+1)>0.) then - h_stack(:) = h(i+ishift,j+1,:) + dz_stack(:) = dz(i+ishift,j+1,:) call remapping_core_h(OBC%remap_CS, & - segment%field(m)%nk_src,segment%field(m)%dz_src(I,j,:), & + segment%field(m)%nk_src, segment%field(m)%dz_src(I,j,:), & segment%field(m)%buffer_src(I,J,:), & - GV%ke, h_stack, segment%field(m)%buffer_dst(I,J,:), & - h_neglect, h_neglect_edge) + GV%ke, dz_stack, segment%field(m)%buffer_dst(I,J,:), & + dz_neglect, dz_neglect_edge) endif enddo else do j=js_obc+1,je_obc ! Using the h remapping approach ! Pretty sure we need to check for source/target grid consistency here - segment%field(m)%buffer_dst(I,j,:)=0.0 ! initialize remap destination buffer + segment%field(m)%buffer_dst(I,j,:) = 0.0 ! initialize remap destination buffer if (G%mask2dCu(I,j)>0.) then - net_H_src = sum( segment%field(m)%dz_src(I,j,:) ) - net_H_int = sum( h(i+ishift,j,:) ) - scl_fac = net_H_int / net_H_src + net_dz_src = sum( segment%field(m)%dz_src(I,j,:) ) + net_dz_int = sum( dz(i+ishift,j,:) ) + scl_fac = net_dz_int / net_dz_src call remapping_core_h(OBC%remap_CS, & - segment%field(m)%nk_src, scl_fac*segment%field(m)%dz_src(I,j,:), & + segment%field(m)%nk_src, scl_fac*segment%field(m)%dz_src(I,j,:), & segment%field(m)%buffer_src(I,j,:), & - GV%ke, h(i+ishift,j,:), segment%field(m)%buffer_dst(I,j,:), & - h_neglect, h_neglect_edge) + GV%ke, dz(i+ishift,j,:), segment%field(m)%buffer_dst(I,j,:), & + dz_neglect, dz_neglect_edge) endif enddo endif @@ -4208,46 +4215,46 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) if (segment%field(m)%name == 'U' .or. segment%field(m)%name == 'DUDY') then ! Do q points for the whole segment do I=max(is_obc,isd),min(ie_obc,ied-1) - segment%field(m)%buffer_dst(I,J,:)=0.0 ! initialize remap destination buffer + segment%field(m)%buffer_dst(I,J,:) = 0.0 ! initialize remap destination buffer if (G%mask2dCv(i,J)>0. .and. G%mask2dCv(i+1,J)>0.) then ! Using the h remapping approach ! Pretty sure we need to check for source/target grid consistency here - h_stack(:) = 0.5*(h(i,j+jshift,:) + h(i+1,j+jshift,:)) + dz_stack(:) = 0.5*(dz(i,j+jshift,:) + dz(i+1,j+jshift,:)) call remapping_core_h(OBC%remap_CS, & - segment%field(m)%nk_src,segment%field(m)%dz_src(I,J,:), & + segment%field(m)%nk_src, segment%field(m)%dz_src(I,J,:), & segment%field(m)%buffer_src(I,J,:), & - GV%ke, h_stack, segment%field(m)%buffer_dst(I,J,:), & - h_neglect, h_neglect_edge) + GV%ke, dz_stack, segment%field(m)%buffer_dst(I,J,:), & + dz_neglect, dz_neglect_edge) elseif (G%mask2dCv(i,J)>0.) then - h_stack(:) = h(i,j+jshift,:) + dz_stack(:) = dz(i,j+jshift,:) call remapping_core_h(OBC%remap_CS, & - segment%field(m)%nk_src,segment%field(m)%dz_src(I,J,:), & + segment%field(m)%nk_src, segment%field(m)%dz_src(I,J,:), & segment%field(m)%buffer_src(I,J,:), & - GV%ke, h_stack, segment%field(m)%buffer_dst(I,J,:), & - h_neglect, h_neglect_edge) + GV%ke, dz_stack, segment%field(m)%buffer_dst(I,J,:), & + dz_neglect, dz_neglect_edge) elseif (G%mask2dCv(i+1,J)>0.) then - h_stack(:) = h(i+1,j+jshift,:) + dz_stack(:) = dz(i+1,j+jshift,:) call remapping_core_h(OBC%remap_CS, & - segment%field(m)%nk_src,segment%field(m)%dz_src(I,J,:), & + segment%field(m)%nk_src, segment%field(m)%dz_src(I,J,:), & segment%field(m)%buffer_src(I,J,:), & - GV%ke, h_stack, segment%field(m)%buffer_dst(I,J,:), & - h_neglect, h_neglect_edge) + GV%ke, dz_stack, segment%field(m)%buffer_dst(I,J,:), & + dz_neglect, dz_neglect_edge) endif enddo else do i=is_obc+1,ie_obc ! Using the h remapping approach ! Pretty sure we need to check for source/target grid consistency here - segment%field(m)%buffer_dst(i,J,:)=0.0 ! initialize remap destination buffer + segment%field(m)%buffer_dst(i,J,:) = 0.0 ! initialize remap destination buffer if (G%mask2dCv(i,J)>0.) then - net_H_src = sum( segment%field(m)%dz_src(i,J,:) ) - net_H_int = sum( h(i,j+jshift,:) ) - scl_fac = net_H_int / net_H_src + net_dz_src = sum( segment%field(m)%dz_src(i,J,:) ) + net_dz_int = sum( dz(i,j+jshift,:) ) + scl_fac = net_dz_int / net_dz_src call remapping_core_h(OBC%remap_CS, & - segment%field(m)%nk_src, scl_fac*segment%field(m)%dz_src(i,J,:), & + segment%field(m)%nk_src, scl_fac* segment%field(m)%dz_src(i,J,:), & segment%field(m)%buffer_src(i,J,:), & - GV%ke, h(i,j+jshift,:), segment%field(m)%buffer_dst(i,J,:), & - h_neglect, h_neglect_edge) + GV%ke, dz(i,j+jshift,:), segment%field(m)%buffer_dst(i,J,:), & + dz_neglect, dz_neglect_edge) endif enddo endif @@ -4496,7 +4503,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) endif elseif (trim(segment%field(m)%genre) == 'obgc') then nt=get_tracer_index(segment,trim(segment%field(m)%name)) - if(nt .lt. 0) then + if (nt < 0) then call MOM_error(FATAL,"update_OBC_segment_data: Did not find tracer "//trim(segment%field(m)%name)) endif if (allocated(segment%field(m)%buffer_dst)) then @@ -4516,7 +4523,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) endif enddo ! end field loop - deallocate(h_stack) + deallocate(dz_stack) deallocate(normal_trans_bt) enddo ! end segment loop @@ -4698,16 +4705,19 @@ subroutine register_segment_tracer(tr_ptr, ntr_index, param_file, GV, segment, & type(OBC_segment_type), intent(inout) :: segment !< current segment data structure real, optional, intent(in) :: OBC_scalar !< If present, use scalar value for segment tracer !! inflow concentration, including any rescaling to - !! put the tracer concentration into its internal units. + !! put the tracer concentration into its internal units, + !! like [S ~> ppt] for salinity. logical, optional, intent(in) :: OBC_array !< If true, use array values for segment tracer !! inflow concentration. real, optional, intent(in) :: scale !< A scaling factor that should be used with any - !! data that is read in, to convert it to the internal - !! units of this tracer. + !! data that is read in to convert it to the internal + !! units of this tracer, in units like [S ppt-1 ~> 1] + !! for salinity. integer, optional, intent(in) :: fd_index !< index of segment tracer in the input field ! Local variables - real :: rescale ! A multiplicative correction to the scaling factor. + real :: rescale ! A multiplicatively corrected scaling factor, in units like [S ppt-1 ~> 1] for + ! salinity, or other various units depending on what rescaling has occurred previously. integer :: ntseg, m, isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB character(len=256) :: mesg ! Message for error messages. @@ -4824,8 +4834,8 @@ subroutine set_obgc_segments_props(OBC,tr_name,obc_src_file_name,obc_src_field_n character(len=*), intent(in) :: tr_name !< Tracer name character(len=*), intent(in) :: obc_src_file_name !< OBC source file name character(len=*), intent(in) :: obc_src_field_name !< name of the field in the source file - real, intent(in) :: lfac_in !< factors for tracer reservoir length scales - real, intent(in) :: lfac_out !< factors for tracer reservoir length scales + real, intent(in) :: lfac_in !< factors for tracer reservoir inbound length scales [nondim] + real, intent(in) :: lfac_out !< factors for tracer reservoir outbound length scales [nondim] type(external_tracers_segments_props),pointer :: node_ptr => NULL() !pointer to type that keeps ! the tracer segment properties @@ -4848,8 +4858,8 @@ subroutine get_obgc_segments_props(node, tr_name,obc_src_file_name,obc_src_field character(len=*), intent(out) :: tr_name !< Tracer name character(len=*), intent(out) :: obc_src_file_name !< OBC source file name character(len=*), intent(out) :: obc_src_field_name !< name of the field in the source file - real, intent(out) :: lfac_in !< multiplicative factor for inbound reservoir length scale - real, intent(out) :: lfac_out !< multiplicative factor for outbound reservoir length scale + real, intent(out) :: lfac_in !< multiplicative factor for inbound reservoir length scale [nondim] + real, intent(out) :: lfac_out !< multiplicative factor for outbound reservoir length scale [nondim] tr_name = trim(node%tracer_name) obc_src_file_name = trim(node%tracer_src_file) obc_src_field_name = trim(node%tracer_src_field) @@ -4890,13 +4900,14 @@ subroutine fill_obgc_segments(G, GV, OBC, tr_ptr, tr_name) type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure type(ocean_OBC_type), pointer :: OBC !< Open boundary structure - real, dimension(:,:,:), pointer :: tr_ptr !< Pointer to tracer field - character(len=*), intent(in) :: tr_name!< Tracer name + real, dimension(:,:,:), pointer :: tr_ptr !< Pointer to tracer field in scaled concentration + !! units, like [S ~> ppt] for salinity. + character(len=*), intent(in) :: tr_name !< Tracer name ! Local variables integer :: isd, ied, IsdB, IedB, jsd, jed, JsdB, JedB, n, nz, nt integer :: i, j, k type(OBC_segment_type), pointer :: segment => NULL() ! pointer to segment type list - real :: I_scale + real :: I_scale ! A factor that unscales the internal units of a tracer, like [ppt S-1 ~> 1] for salinity if (.not. associated(OBC)) return call pass_var(tr_ptr, G%Domain) @@ -4905,7 +4916,7 @@ subroutine fill_obgc_segments(G, GV, OBC, tr_ptr, tr_name) segment => OBC%segment(n) if (.not. segment%on_pe) cycle nt=get_tracer_index(segment,tr_name) - if(nt .lt. 0) then + if (nt < 0) then call MOM_error(FATAL,"fill_obgc_segments: Did not find tracer "// tr_name) endif isd = segment%HI%isd ; ied = segment%HI%ied @@ -5414,7 +5425,7 @@ subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg) do m=1,segment%tr_Reg%ntseg ntr_id = segment%tr_reg%Tr(m)%ntr_index fd_id = segment%tr_reg%Tr(m)%fd_index - if(fd_id == -1) then + if (fd_id == -1) then resrv_lfac_out = 1.0 resrv_lfac_in = 1.0 else @@ -5458,7 +5469,7 @@ subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg) do m=1,segment%tr_Reg%ntseg ntr_id = segment%tr_reg%Tr(m)%ntr_index fd_id = segment%tr_reg%Tr(m)%fd_index - if(fd_id == -1) then + if (fd_id == -1) then resrv_lfac_out = 1.0 resrv_lfac_in = 1.0 else @@ -5500,7 +5511,8 @@ subroutine remap_OBC_fields(G, GV, h_old, h_new, OBC, PCM_cell) ! Local variables type(OBC_segment_type), pointer :: segment => NULL() ! A pointer to the various segments, used just for shorthand. - real :: tr_column(GV%ke) ! A column of updated tracer concentrations [CU ~> Conc] + real :: tr_column(GV%ke) ! A column of updated tracer concentrations in internally scaled units. + ! For salinity the units would be [S ~> ppt]. real :: r_norm_col(GV%ke) ! A column of updated radiation rates, in grid points per timestep [nondim] real :: rxy_col(GV%ke) ! A column of updated radiation rates for oblique OBCs [L2 T-2 ~> m2 s-2] real :: h1(GV%ke) ! A column of source grid layer thicknesses [H ~> m or kg m-2] @@ -5706,14 +5718,14 @@ subroutine adjustSegmentEtaToFitBathymetry(G, GV, US, segment,fld) allocate(eta(is:ie,js:je,nz+1)) contractions=0; dilations=0 do j=js,je ; do i=is,ie - eta(i,j,1)=0.0 ! segment data are assumed to be located on a static grid + eta(i,j,1) = 0.0 ! segment data are assumed to be located on a static grid ! For remapping calls, the entire column will be dilated ! by a factor equal to the ratio of the sum of the geopotential referenced ! source data thicknesses, and the current model thicknesses. This could be ! an issue to be addressed, for instance if we are placing open boundaries ! under ice shelf cavities. do k=2,nz+1 - eta(i,j,k)=eta(i,j,k-1)-segment%field(fld)%dz_src(i,j,k-1) + eta(i,j,k) = eta(i,j,k-1) - segment%field(fld)%dz_src(i,j,k-1) enddo ! The normal slope at the boundary is zero by a ! previous call to open_boundary_impose_normal_slope @@ -5741,7 +5753,7 @@ subroutine adjustSegmentEtaToFitBathymetry(G, GV, US, segment,fld) dilations = dilations + 1 ! expand bottom-most cell only eta(i,j,nz+1) = -segment%dZtot(i,j) - segment%field(fld)%dz_src(i,j,nz)= eta(i,j,nz)-eta(i,j,nz+1) + segment%field(fld)%dz_src(i,j,nz) = eta(i,j,nz) - eta(i,j,nz+1) ! if (eta(i,j,1) <= eta(i,j,nz+1)) then ! do k=1,nz ; segment%field(fld)%dz_src(i,j,k) = (eta(i,j,1) + G%bathyT(i,j)) / real(nz) ; enddo ! else @@ -5750,10 +5762,6 @@ subroutine adjustSegmentEtaToFitBathymetry(G, GV, US, segment,fld) ! endif !do k=nz,2,-1 ; eta(i,j,K) = eta(i,j,K+1) + segment%field(fld)%dz_src(i,j,k) ; enddo endif - ! Now convert thicknesses to units of H. - do k=1,nz - segment%field(fld)%dz_src(i,j,k) = segment%field(fld)%dz_src(i,j,k)*GV%Z_to_H - enddo enddo ; enddo ! can not do communication call here since only PEs on the current segment are here From d0e9c25a08db5ddf1e13f0b875f26259095e3d31 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 24 Jan 2024 11:28:21 -0500 Subject: [PATCH 02/15] *Corrected CFC diagnostics and non-Bous CFC fluxes Removed the duplicative multiplication by the Boussinesq density and related unit conversion factors in the calculation of the diagnosed mole concentration of the "cfc11" and "cfc12" diagnostics. This unit rescaling is duplicative of the conversion factors in the register_diag_field calls for these diagnostics (at about line 263). This will change these two diagnostics by a factor of RHO_0. This bug arose from separate and independent changes coming from GFDL and NCAR to add the same missing rescaling factor, and it might have gone undetected because of the vast range of valid values for CFC concentrations or because all of our regression testing uses calendar dates that precede the invention and first release of CFCs in the 1930s. Although this commit will correct the code that is on the dev/gfdl branch, care will have to be taken when merging in modified versions of this file from other branches that might have also corrected this bug to avoid reintroducing it (perhaps in reverse). Also replaced the expression for the CFC flux rescaling factor for use with KPP_NonLocalTransport with GV%RZ_to_H, which is equivalent to the previous expression in Boussinesq mode but avoids a multiplication and division by the Boussinesq reference density in non-Boussinesq mode. All Boussinesq answers are identical, but some non-Boussinesq CFCs could be altered in the last bits. This PR will change diagnostics of CFC concentrations by about 3 orders of magnitude. The internal representation of the CFCs is bitwise identical in all Boussinesq cases, but they may change at roundoff in non-Boussinesq cases that use KPP. --- src/tracer/MOM_CFC_cap.F90 | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/tracer/MOM_CFC_cap.F90 b/src/tracer/MOM_CFC_cap.F90 index 16506b41c3..489948a63c 100644 --- a/src/tracer/MOM_CFC_cap.F90 +++ b/src/tracer/MOM_CFC_cap.F90 @@ -361,7 +361,6 @@ subroutine CFC_cap_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US, C ! Local variables real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_work ! Used so that h can be modified [H ~> m or kg m-2] - real :: flux_scale ! A dimensional rescaling factor for fluxes [H R-1 Z-1 ~> m3 kg-1 or nondim] integer :: i, j, k, is, ie, js, je, nz, m is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke @@ -371,13 +370,11 @@ subroutine CFC_cap_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US, C ! Compute KPP nonlocal term if necessary if (present(KPP_CSp)) then if (associated(KPP_CSp) .and. present(nonLocalTrans)) then - flux_scale = GV%Z_to_H / GV%rho0 - do m=1,NTR call KPP_NonLocalTransport(KPP_CSp, G, GV, h_old, nonLocalTrans, & CS%CFC_data(m)%sfc_flux(:,:), dt, CS%diag, & CS%CFC_data(m)%tr_ptr, CS%CFC_data(m)%conc(:,:,:), & - flux_scale=flux_scale) + flux_scale=GV%RZ_to_H) enddo endif endif From fd5696ba07d81c9841c7710de4dd7553ca22ee84 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Thu, 29 Feb 2024 12:05:59 -0500 Subject: [PATCH 03/15] Pinning package requirements for readthedocs Minor version updates to multiple packages used in the generation of the documentation introduce dependencies needing Sphinx >= 5.0, which breaks the sphinx extensions we use in documenting the MOM6 APIs. I have added versions for all the packages needed to keep things working with Sphinx4 for now, but we really do need to find a way to work with the newer versions. More pinningaof of requirements for readthedocs More pinning of requirements for readthedocs --- docs/requirements.txt | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/docs/requirements.txt b/docs/requirements.txt index ff627c61c7..b38dbc34b7 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -10,3 +10,9 @@ six future # Old Sphinx requires an old Jinja2 jinja2<3.1 +sphinxcontrib_applehelp<1.0.8 +sphinxcontrib_devhelp<1.0.6 +sphinxcontrib_htmlhelp<2.0.5 +sphinxcontrib_qthelp<1.0.7 +sphinxcontrib_serializinghtml<1.0.7 +alabaster<0.7.14 From f9da6733155d60b152ae9df35b028a093b341ad2 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Tue, 16 Jan 2024 16:21:40 -0500 Subject: [PATCH 04/15] Create BUILD and WORK directory macros This patch introduces two new macros, BUILD and WORK, to permit relocation of the build/ and work/ directories. It also makes the following smaller changes: * deps/ is now defined by the DEPS macro. If unset, deps/ is placed in the BUILD directory. * results/ is moved into WORK. * Compiler flags which track directories now use $(abspath ...) to allow for arbitrary paths. * GitHub CI paths were adjusted to support these new settings. * DO_* flags are now used as on/off with ifdef testing, rather than checking for `true` values. * mkmf macros have been removed from the coupled test config. * The default FMS infra has been changed to FMS2 in all components, including the configure.ac outside of .testing. This work will enable testing of multiple FMS libraries in our CI. --- .github/actions/testing-setup/action.yml | 2 +- .testing/Makefile | 456 ++++++++++++----------- ac/configure.ac | 4 +- 3 files changed, 244 insertions(+), 218 deletions(-) diff --git a/.github/actions/testing-setup/action.yml b/.github/actions/testing-setup/action.yml index a15dd6d0a2..60499d4be1 100644 --- a/.github/actions/testing-setup/action.yml +++ b/.github/actions/testing-setup/action.yml @@ -28,7 +28,7 @@ runs: run: | echo "::group::Compile FMS library" cd .testing - REPORT_ERROR_LOGS=true make deps/lib/libFMS.a -s -j + REPORT_ERROR_LOGS=true make build/deps/lib/libFMS.a -s -j echo "::endgroup::" - name: Compile MOM6 in symmetric memory mode diff --git a/.testing/Makefile b/.testing/Makefile index f7101a3463..a1461fe7ab 100644 --- a/.testing/Makefile +++ b/.testing/Makefile @@ -46,10 +46,10 @@ # LDFLAGS_USER User-defined linker flags (used for all MOM/FMS builds) # # Experiment Configuration: -# BUILDS Executables to be built by `make` or `make all` +# EXECS Executables to be built by `make` or `make all` # CONFIGS Model configurations to test (default: `tc*`) -# TESTS Tests to run # DIMS Dimensional scaling tests +# TESTS Tests to run # # Regression repository ("target") configuration: # MOM_TARGET_SLUG URL slug (minus domain) of the target repo @@ -57,19 +57,19 @@ # MOM_TARGET_LOCAL_BRANCH Target branch name # (NOTE: These would typically be configured by a CI.) # -# Paths for stages: -# WORKSPACE Location to place work/ and results/ directories (i.e. where to run the model) -# -#---- +# Output paths: +# BUILD Compiled executables and libraries +# DEPS Compiled dependencies +# WORK Test model output # TODO: POSIX shell compatibility SHELL = bash -# No implicit rules -.SUFFIXES: +# No implicit rules, suffixes, or variables +MAKEFLAGS += -rR -# No implicit variables -MAKEFLAGS += -R +# Determine the MOM6 autoconf srcdir +AC_SRCDIR := $(dir $(abspath $(lastword $(MAKEFILE_LIST))))../ac # User-defined configuration -include config.mk @@ -104,10 +104,7 @@ FCFLAGS_FMS ?= $(FCFLAGS_DEBUG) LDFLAGS_COVERAGE ?= --coverage LDFLAGS_USER ?= -# Set to `true` to require identical results from DEBUG and REPRO builds -# NOTE: Many compilers (Intel, GCC on ARM64) do not produce identical results -# across DEBUG and REPRO builds (as defined below), so we disable on -# default. +# Set to verify identical DEBUG and REPRO results DO_REPRO_TESTS ?= # Enable profiling @@ -116,9 +113,12 @@ DO_PROFILE ?= # Enable code coverage runs DO_COVERAGE ?= -# Enable code coverage runs +# Enable unit tests DO_UNIT_TESTS ?= +# Check for regressions with target branch +DO_REGRESSION_TESTS ?= + # Report failure if coverage report is not uploaded REQUIRE_COVERAGE_UPLOAD ?= @@ -128,52 +128,69 @@ REPORT_ERROR_LOGS ?= # Time measurement (configurable by the CI) TIME ?= time +# Legacy external work directory +#WORKSPACE ?= +WORKSPACE ?= . + +# Set directories for build/ and work/ +#BUILD ?= $(WORKSPACE)build +#DEPS ?= $(BUILD)/deps +#WORK ?= $(WORKSPACE)work +BUILD ?= $(WORKSPACE)/build +DEPS ?= $(BUILD)/deps +WORK ?= $(WORKSPACE)/work # Experiment configuration -BUILDS ?= symmetric/MOM6 asymmetric/MOM6 openmp/MOM6 +EXECS ?= symmetric/MOM6 asymmetric/MOM6 openmp/MOM6 CONFIGS ?= $(wildcard tc*) -TESTS ?= grid layout rotate restart openmp nan $(foreach d,$(DIMS),dim.$(d)) DIMS ?= t l h z q r +TESTS ?= grid layout rotate restart openmp nan $(foreach d,$(DIMS),dim.$(d)) + +# Unit test executables +UNIT_EXECS ?= \ + $(basename $(notdir $(wildcard ../config_src/drivers/unit_tests/*.F90))) + +# Timing test executables +TIMING_EXECS ?= \ + $(basename $(notdir $(wildcard ../config_src/drivers/timing_tests/*.F90))) -# Default is to place work/ and results/ in current directory -WORKSPACE ?= . #--- # Test configuration -# REPRO tests enable reproducibility with optimization, and often do not match -# the DEBUG results in older GCCs and vendor compilers, so we can optionally -# disable them. -ifeq ($(DO_REPRO_TESTS), true) - BUILDS += repro/MOM6 +# Set if either DO_COVERAGE or DO_UNIT_TESTS is set +run_unit_tests = + +# REPRO and DEBUG equivalence +ifdef DO_REPRO_TESTS + EXECS += repro/MOM6 TESTS += repro endif # Profiling -ifeq ($(DO_PROFILE), true) - BUILDS += opt/MOM6 opt_target/MOM6 +ifdef DO_PROFILE + EXECS += opt/MOM6 opt_target/MOM6 endif # Coverage -ifeq ($(DO_COVERAGE), true) - BUILDS += cov/MOM6 +ifdef DO_COVERAGE + EXECS += cov/MOM6 + run_unit_execs = yes endif -# Unit testing (or coverage) -UNIT_EXECS ?= $(basename $(notdir $(wildcard ../config_src/drivers/unit_tests/*.F90) ) ) -TIMING_EXECS ?= $(basename $(notdir $(wildcard ../config_src/drivers/timing_tests/*.F90) ) ) -ifneq (X$(DO_COVERAGE)$(DO_UNIT_TESTS)X, XX) - BUILDS += $(foreach e, $(UNIT_EXECS), unit/$(e)) +# Unit test executables +ifdef DO_UNIT_TESTS + run_unit_tests = yes endif -ifeq ($(DO_PROFILE), false) - BUILDS += opt/MOM6 opt_target/MOM6 +# If either coverage or unit tests are enabled, build the unit test execs +ifdef run_unit_tests + EXECS += $(foreach e, $(UNIT_EXECS), unit/$(e)) endif - -DO_REGRESSION_TESTS ?= -ifeq ($(DO_REGRESSION_TESTS), true) - BUILDS += target/MOM6 +# Regression testing +ifdef DO_REGRESSION_TESTS + EXECS += target/MOM6 TESTS += regression MOM_TARGET_SLUG ?= NOAA-GFDL/MOM6 @@ -182,13 +199,14 @@ ifeq ($(DO_REGRESSION_TESTS), true) MOM_TARGET_LOCAL_BRANCH ?= dev/gfdl MOM_TARGET_BRANCH := origin/$(MOM_TARGET_LOCAL_BRANCH) - TARGET_CODEBASE = build/target_codebase + TARGET_CODEBASE = $(BUILD)/target_codebase else MOM_TARGET_URL = MOM_TARGET_BRANCH = TARGET_CODEBASE = endif + # List of source files to link this Makefile's dependencies to model Makefiles # Assumes a depth of two, and the following extensions: F90 inc c h # (1): Root directory @@ -202,31 +220,30 @@ MOM_SOURCE = \ $(wildcard ../config_src/ext*/*/*.F90) TARGET_SOURCE = \ - $(call SOURCE,build/target_codebase/src) \ - $(wildcard build/target_codebase/config_src/drivers/solo_driver/*.F90) \ - $(wildcard build/target_codebase/config_src/ext*/*.F90) + $(call SOURCE,$(BUILD)/target_codebase/src) \ + $(wildcard $(BUILD)/target_codebase/config_src/drivers/solo_driver/*.F90) \ + $(wildcard $(BUILD)target_codebase/config_src/ext*/*.F90) -# NOTE: Current default framework is FMS1, but this could change. -ifeq ($(FRAMEWORK), fms2) - MOM_SOURCE +=$(wildcard ../config_src/infra/FMS2/*.F90) - TARGET_SOURCE += $(wildcard build/target_codebase/config_src/infra/FMS2/*.F90) -else +ifeq ($(FRAMEWORK), fms1) MOM_SOURCE += $(wildcard ../config_src/infra/FMS1/*.F90) - TARGET_SOURCE += $(wildcard build/target_codebase/config_src/infra/FMS1/*.F90) + TARGET_SOURCE += $(wildcard $(BUILD)/target_codebase/config_src/infra/FMS1/*.F90) +else + MOM_SOURCE +=$(wildcard ../config_src/infra/FMS2/*.F90) + TARGET_SOURCE += $(wildcard $(BUILD)/target_codebase/config_src/infra/FMS2/*.F90) endif -FMS_SOURCE = $(call SOURCE,deps/fms/src) +FMS_SOURCE = $(call SOURCE,$(DEPS)/fms/src) -#--- -# Rules + +## Rules .PHONY: all build.regressions build.prof -all: $(foreach b,$(BUILDS),build/$(b)) -build.regressions: $(foreach b,symmetric target,build/$(b)/MOM6) -build.prof: $(foreach b,opt opt_target,build/$(b)/MOM6) +all: $(foreach b,$(EXECS),$(BUILD)/$(b)) +build.regressions: $(foreach b,symmetric target,$(BUILD)/$(b)/MOM6) +build.prof: $(foreach b,opt opt_target,$(BUILD)/$(b)/MOM6) # Executable -.PRECIOUS: $(foreach b,$(BUILDS),build/$(b)) +.PRECIOUS: $(foreach b,$(EXECS),$(BUILD)/$(b)) # Compiler flags @@ -234,9 +251,9 @@ build.prof: $(foreach b,opt opt_target,build/$(b)/MOM6) # .testing dependencies # TODO: We should probably build TARGET with the FMS that it was configured # to use. But for now we use the same FMS over all builds. -FCFLAGS_DEPS = -I../../deps/include -LDFLAGS_DEPS = -L../../deps/lib -PATH_DEPS = PATH="${PATH}:../../deps/bin" +FCFLAGS_DEPS = -I$(abspath $(DEPS)/include) +LDFLAGS_DEPS = -L$(abspath $(DEPS)/lib) +PATH_DEPS = PATH="${PATH}:$(abspath $(DEPS)/bin)" # Define the build targets in terms of the traditional DEBUG/REPRO/etc labels @@ -254,82 +271,96 @@ COV_LDFLAGS := LDFLAGS="$(LDFLAGS_COVERAGE) $(LDFLAGS_DEPS) $(LDFLAGS_USER)" # Environment variable configuration MOM_ENV := $(PATH_FMS) -build/symmetric/Makefile: MOM_ENV += $(SYMMETRIC_FCFLAGS) $(MOM_LDFLAGS) -build/asymmetric/Makefile: MOM_ENV += $(ASYMMETRIC_FCFLAGS) $(MOM_LDFLAGS) \ - MOM_MEMORY=../../../config_src/memory/dynamic_nonsymmetric/MOM_memory.h -build/repro/Makefile: MOM_ENV += $(REPRO_FCFLAGS) $(MOM_LDFLAGS) -build/openmp/Makefile: MOM_ENV += $(OPENMP_FCFLAGS) $(MOM_LDFLAGS) -build/target/Makefile: MOM_ENV += $(TARGET_FCFLAGS) $(MOM_LDFLAGS) -build/opt/Makefile: MOM_ENV += $(OPT_FCFLAGS) $(MOM_LDFLAGS) -build/opt_target/Makefile: MOM_ENV += $(OPT_FCFLAGS) $(MOM_LDFLAGS) -build/coupled/Makefile: MOM_ENV += $(SYMMETRIC_FCFLAGS) $(MOM_LDFLAGS) -build/nuopc/Makefile: MOM_ENV += $(SYMMETRIC_FCFLAGS) $(MOM_LDFLAGS) -build/cov/Makefile: MOM_ENV += $(COV_FCFLAGS) $(COV_LDFLAGS) -build/unit/Makefile: MOM_ENV += $(COV_FCFLAGS) $(COV_LDFLAGS) -build/timing/Makefile: MOM_ENV += $(OPT_FCFLAGS) $(MOM_LDFLAGS) +$(BUILD)/symmetric/Makefile: MOM_ENV += $(SYMMETRIC_FCFLAGS) $(MOM_LDFLAGS) +$(BUILD)/asymmetric/Makefile: MOM_ENV += $(ASYMMETRIC_FCFLAGS) $(MOM_LDFLAGS) \ + MOM_MEMORY=$(AC_SRCDIR)/../config_src/memory/dynamic_nonsymmetric/MOM_memory.h +$(BUILD)/repro/Makefile: MOM_ENV += $(REPRO_FCFLAGS) $(MOM_LDFLAGS) +$(BUILD)/openmp/Makefile: MOM_ENV += $(OPENMP_FCFLAGS) $(MOM_LDFLAGS) +$(BUILD)/target/Makefile: MOM_ENV += $(TARGET_FCFLAGS) $(MOM_LDFLAGS) +$(BUILD)/opt/Makefile: MOM_ENV += $(OPT_FCFLAGS) $(MOM_LDFLAGS) +$(BUILD)/opt_target/Makefile: MOM_ENV += $(OPT_FCFLAGS) $(MOM_LDFLAGS) +$(BUILD)/coupled/Makefile: MOM_ENV += $(SYMMETRIC_FCFLAGS) $(MOM_LDFLAGS) +$(BUILD)/nuopc/Makefile: MOM_ENV += $(SYMMETRIC_FCFLAGS) $(MOM_LDFLAGS) +$(BUILD)/cov/Makefile: MOM_ENV += $(COV_FCFLAGS) $(COV_LDFLAGS) +$(BUILD)/unit/Makefile: MOM_ENV += $(COV_FCFLAGS) $(COV_LDFLAGS) +$(BUILD)/timing/Makefile: MOM_ENV += $(OPT_FCFLAGS) $(MOM_LDFLAGS) # Configure script flags MOM_ACFLAGS := --with-framework=$(FRAMEWORK) -build/openmp/Makefile: MOM_ACFLAGS += --enable-openmp -build/coupled/Makefile: MOM_ACFLAGS += --with-driver=FMS_cap -build/nuopc/Makefile: MOM_ACFLAGS += --with-driver=nuopc_cap -build/unit/Makefile: MOM_ACFLAGS += --with-driver=unit_tests -build/timing/Makefile: MOM_ACFLAGS += --with-driver=timing_tests - -# Fetch regression target source code -build/target/Makefile: | $(TARGET_CODEBASE) -build/opt_target/Makefile: | $(TARGET_CODEBASE) - - -# Define source code dependencies -build/target_codebase/configure: $(TARGET_SOURCE) +$(BUILD)/openmp/Makefile: MOM_ACFLAGS += --enable-openmp +$(BUILD)/coupled/Makefile: MOM_ACFLAGS += --with-driver=FMS_cap +$(BUILD)/nuopc/Makefile: MOM_ACFLAGS += --with-driver=nuopc_cap +$(BUILD)/unit/Makefile: MOM_ACFLAGS += --with-driver=unit_tests +$(BUILD)/timing/Makefile: MOM_ACFLAGS += --with-driver=timing_tests # Build executables -build/unit/test_%: build/unit/Makefile FORCE +$(BUILD)/unit/test_%: $(BUILD)/unit/Makefile FORCE cd $(@D) && $(TIME) $(MAKE) $(@F) -j -build/unit/Makefile: $(foreach e,$(UNIT_EXECS),../config_src/drivers/unit_tests/$(e).F90) -build/timing/time_%: build/timing/Makefile FORCE +$(BUILD)/unit/Makefile: $(foreach e,$(UNIT_EXECS),../config_src/drivers/unit_tests/$(e).F90) +$(BUILD)/timing/time_%: $(BUILD)/timing/Makefile FORCE cd $(@D) && $(TIME) $(MAKE) $(@F) -j -build/timing/Makefile: $(foreach e,$(TIMING_EXECS),../config_src/drivers/timing_tests/$(e).F90) -build/%/MOM6: build/%/Makefile FORCE +$(BUILD)/timing/Makefile: $(foreach e,$(TIMING_EXECS),../config_src/drivers/timing_tests/$(e).F90) +$(BUILD)/%/MOM6: $(BUILD)/%/Makefile FORCE cd $(@D) && $(TIME) $(MAKE) $(@F) -j -FORCE: ; +FORCE: -# Use autoconf to construct the Makefile for each target -.PRECIOUS: build/%/Makefile -build/%/Makefile: ../ac/configure ../ac/Makefile.in deps/lib/libFMS.a - mkdir -p $(@D) - cd $(@D) \ - && $(MOM_ENV) ../../../ac/configure $(MOM_ACFLAGS) \ - || (cat config.log && false) +## Use autoconf to construct the Makefile for each target +# TODO: This could all be moved to a top-level MOM6 Makefile +.PRECIOUS: $(BUILD)/%/Makefile +.PRECIOUS: $(BUILD)/%/Makefile.in +.PRECIOUS: $(BUILD)/%/configure +.PRECIOUS: $(BUILD)/%/config.status +.PRECIOUS: $(BUILD)/%/configure.ac +.PRECIOUS: $(BUILD)/%/m4/ +$(BUILD)/%/Makefile: $(BUILD)/%/Makefile.in $(BUILD)/%/config.status + cd $(@D) && ./config.status -../ac/configure: ../ac/configure.ac ../ac/m4 - autoreconf -i $< +$(BUILD)/%/config.status: $(BUILD)/%/configure $(DEPS)/lib/libFMS.a + cd $(@D) && $(MOM_ENV) ./configure -n --srcdir=$(AC_SRCDIR) $(MOM_ACFLAGS) \ + || (cat config.log && false) +$(BUILD)/%/Makefile.in: ../ac/Makefile.in | $(BUILD)/%/ + cp ../ac/Makefile.in $(@D) + +$(BUILD)/%/configure: $(BUILD)/%/configure.ac $(BUILD)/%/m4/ + autoreconf -if $(@D) + +$(BUILD)/%/configure.ac: ../ac/configure.ac | $(BUILD)/%/ + cp ../ac/configure.ac $(@D) + +$(BUILD)/%/m4/: ../ac/m4/ | $(BUILD)/%/ + cp -r ../ac/m4 $(@D) + +ALL_EXECS = symmetric asymmetric repro openmp target opt opt_target coupled \ + nuopc cov unit timing +$(foreach b,$(ALL_EXECS),$(BUILD)/$(b)/): + mkdir -p $@ # Fetch the regression target codebase -build/target/Makefile build/opt_target/Makefile: \ - $(TARGET_CODEBASE)/ac/configure deps/lib/libFMS.a - mkdir -p $(@D) - cd $(@D) \ - && $(MOM_ENV) ../../$(TARGET_CODEBASE)/ac/configure $(MOM_ACFLAGS) \ + +$(BUILD)/target/config.status: $(BUILD)/target/configure $(DEPS)/lib/libFMS.a + cd $(@D) && $(MOM_ENV) ./configure -n \ + --srcdir=$(abspath $(BUILD))/target_codebase/ac $(MOM_ACFLAGS) \ || (cat config.log && false) +$(BUILD)/target/Makefile.in: | $(TARGET_CODEBASE) $(BUILD)/target/ + cp $(TARGET_CODEBASE)/ac/Makefile.in $(@D) -$(TARGET_CODEBASE)/ac/configure: $(TARGET_CODEBASE) - autoreconf -i $/dev/null)" || rm -rf $(WORKSPACE)/results/$(1) +.PRECIOUS: $(foreach b,$(3),$(WORK)/$(1)/$(b)/ocean.stats) +$(1).$(2): $(foreach b,$(3),$(WORK)/$(1)/$(b)/ocean.stats) + @test "$$(shell ls -A $(WORK)/results/$(1) 2>/dev/null)" || rm -rf $(WORK)/results/$(1) @cmp $$^ || !( \ - mkdir -p $(WORKSPACE)/results/$(1); \ - (diff $$^ | tee $(WORKSPACE)/results/$(1)/ocean.stats.$(2).diff | head -n 20) ; \ + mkdir -p $(WORK)/results/$(1); \ + (diff $$^ | tee $(WORK)/results/$(1)/ocean.stats.$(2).diff | head -n 20) ; \ echo -e "$(FAIL): Solutions $(1).$(2) have changed." \ ) @echo -e "$(PASS): Solutions $(1).$(2) agree." -.PRECIOUS: $(foreach b,$(3),$(WORKSPACE)/work/$(1)/$(b)/chksum_diag) -$(1).$(2).diag: $(foreach b,$(3),$(WORKSPACE)/work/$(1)/$(b)/chksum_diag) +.PRECIOUS: $(foreach b,$(3),$(WORK)/$(1)/$(b)/chksum_diag) +$(1).$(2).diag: $(foreach b,$(3),$(WORK)/$(1)/$(b)/chksum_diag) @cmp $$^ || !( \ - mkdir -p $(WORKSPACE)/results/$(1); \ - (diff $$^ | tee $(WORKSPACE)/results/$(1)/chksum_diag.$(2).diff | head -n 20) ; \ + mkdir -p $(WORK)/results/$(1); \ + (diff $$^ | tee $(WORK)/results/$(1)/chksum_diag.$(2).diff | head -n 20) ; \ echo -e "$(FAIL): Diagnostics $(1).$(2).diag have changed." \ ) @echo -e "$(PASS): Diagnostics $(1).$(2).diag agree." @@ -473,17 +497,17 @@ $(foreach d,$(DIMS),$(eval $(call CMP_RULE,$(1),dim.$(d),symmetric dim.$(d)))) endef $(foreach c,$(CONFIGS),$(eval $(call CONFIG_DIM_RULE,$(c)))) -# Custom comparison rules +# Custom comparison rules # Restart tests only compare the final stat record -.PRECIOUS: $(foreach b,symmetric restart target,$(WORKSPACE)/work/%/$(b)/ocean.stats) -%.restart: $(foreach b,symmetric restart,$(WORKSPACE)/work/%/$(b)/ocean.stats) - @test "$(shell ls -A $(WORKSPACE)/results/$* 2>/dev/null)" || rm -rf $(WORKSPACE)/results/$* +.PRECIOUS: $(foreach b,symmetric restart target,$(WORK)/%/$(b)/ocean.stats) +%.restart: $(foreach b,symmetric restart,$(WORK)/%/$(b)/ocean.stats) + @test "$(shell ls -A $(WORK)/results/$* 2>/dev/null)" || rm -rf $(WORK)/results/$* @cmp $(foreach f,$^,<(tr -s ' ' < $(f) | cut -d ' ' -f3- | tail -n 1)) \ || !( \ - mkdir -p $(WORKSPACE)/results/$*; \ - (diff $^ | tee $(WORKSPACE)/results/$*/chksum_diag.restart.diff | head -n 20) ; \ + mkdir -p $(WORK)/results/$*; \ + (diff $^ | tee $(WORK)/results/$*/chksum_diag.restart.diff | head -n 20) ; \ echo -e "$(FAIL): Solutions $*.restart have changed." \ ) @echo -e "$(PASS): Solutions $*.restart agree." @@ -491,21 +515,22 @@ $(foreach c,$(CONFIGS),$(eval $(call CONFIG_DIM_RULE,$(c)))) # TODO: chksum_diag parsing of restart files # stats rule is unchanged, but we cannot use CMP_RULE to generate it. -%.regression: $(foreach b,symmetric target,$(WORKSPACE)/work/%/$(b)/ocean.stats) - @test "$(shell ls -A $(WORKSPACE)/results/$* 2>/dev/null)" || rm -rf $(WORKSPACE)/results/$* +%.regression: $(foreach b,symmetric target,$(WORK)/%/$(b)/ocean.stats) + @test "$(shell ls -A $(WORK)/results/$* 2>/dev/null)" || rm -rf $(WORK)/results/$* @cmp $^ || !( \ - mkdir -p $(WORKSPACE)/results/$*; \ - (diff $^ | tee $(WORKSPACE)/results/$*/ocean.stats.regression.diff | head -n 20) ; \ + mkdir -p $(WORK)/results/$*; \ + (diff $^ | tee $(WORK)/results/$*/ocean.stats.regression.diff | head -n 20) ; \ echo -e "$(FAIL): Solutions $*.regression have changed." \ ) @echo -e "$(PASS): Solutions $*.regression agree." # Regression testing only checks for changes in existing diagnostics -%.regression.diag: $(foreach b,symmetric target,$(WORKSPACE)/work/%/$(b)/chksum_diag) +.PRECIOUS: $(WORK)/%/target/chksum_diag +%.regression.diag: $(foreach b,symmetric target,$(WORK)/%/$(b)/chksum_diag) @! diff $^ | grep "^[<>]" | grep "^>" > /dev/null \ || ! (\ - mkdir -p $(WORKSPACE)/results/$*; \ - (diff $^ | tee $(WORKSPACE)/results/$*/chksum_diag.regression.diff | head -n 20) ; \ + mkdir -p $(WORK)/results/$*; \ + (diff $^ | tee $(WORK)/results/$*/chksum_diag.regression.diff | head -n 20) ; \ echo -e "$(FAIL): Diagnostics $*.regression.diag have changed." \ ) @cmp $^ || ( \ @@ -534,7 +559,7 @@ tc4/configure: tc4/configure.ac #--- # Test run output files -# Rule to build $(WORKSPACE)/work//{ocean.stats,chksum_diag}. +# Rule to build $(WORK)//{ocean.stats,chksum_diag}. # $(1): Test configuration name # $(2): Executable type # $(3): Enable coverage flag @@ -543,15 +568,14 @@ tc4/configure: tc4/configure.ac # $(6): Number of MPI ranks define STAT_RULE -$(WORKSPACE)/work/%/$(1)/ocean.stats $(WORKSPACE)/work/%/$(1)/chksum_diag: build/$(2)/MOM6 | preproc +$(WORK)/%/$(1)/ocean.stats $(WORK)/%/$(1)/chksum_diag: $(BUILD)/$(2)/MOM6 | preproc @echo "Running test $$*.$(1)..." mkdir -p $$(@D) cp -RL $$*/* $$(@D) - mkdir -p $$(@D)/RESTART echo -e "$(4)" > $$(@D)/MOM_override - rm -f $(WORKSPACE)/results/$$*/std.$(1).{out,err} + rm -f $(WORK)/results/$$*/std.$(1).{out,err} cd $$(@D) \ - && $(TIME) $(5) $(MPIRUN) -n $(6) $(abspath $$<) 2> std.err > std.out \ + && $(TIME) $(5) $(MPIRUN) -n $(6) $$(abspath $$<) 2> std.err > std.out \ || !( \ mkdir -p ../../../results/$$*/ ; \ cat std.out | tee ../../../results/$$*/std.$(1).out | tail -n 40 ; \ @@ -561,8 +585,8 @@ $(WORKSPACE)/work/%/$(1)/ocean.stats $(WORKSPACE)/work/%/$(1)/chksum_diag: build ) @echo -e "$(DONE): $$*.$(1); no runtime errors." if [ $(3) ]; then \ - mkdir -p $(WORKSPACE)/results/$$* ; \ - cd build/$(2) ; \ + mkdir -p $(WORK)/results/$$* ; \ + cd $(BUILD)/$(2) ; \ gcov -b *.gcda > gcov.$$*.$(1).out ; \ find -name "*.gcov" -exec sed -i -r 's/^( *[0-9]*)\*:/ \1:/g' {} \; ; \ fi @@ -585,12 +609,12 @@ codecov: .PHONY: report.cov report.cov: run.cov codecov - ./codecov $(CODECOV_TOKEN_ARG) -R build/cov -Z -f "*.gcov" \ - > build/cov/codecov.out \ - 2> build/cov/codecov.err \ + ./codecov $(CODECOV_TOKEN_ARG) -R $(BUILD)/cov -Z -f "*.gcov" \ + > $(BUILD)/cov/codecov.out \ + 2> $(BUILD)/cov/codecov.err \ && echo -e "${MAGENTA}Report uploaded to codecov.${RESET}" \ || { \ - cat build/cov/codecov.err ; \ + cat $(BUILD)/cov/codecov.err ; \ echo -e "${RED}Failed to upload report.${RESET}" ; \ if [ "$(REQUIRE_COVERAGE_UPLOAD)" = true ] ; then false ; fi ; \ } @@ -620,7 +644,7 @@ $(eval $(call STAT_RULE,cov,cov,true,,,1)) # 2. Convert DAYMAX from TIMEUNIT to seconds # 3. Apply seconds to `ocean_solo_nml` inside input.nml. # NOTE: Assumes that runtime set by DAYMAX, will fail if set by input.nml -$(WORKSPACE)/work/%/restart/ocean.stats: build/symmetric/MOM6 | preproc +$(WORK)/%/restart/ocean.stats: $(BUILD)/symmetric/MOM6 | preproc rm -rf $(@D) mkdir -p $(@D) cp -RL $*/* $(@D) @@ -634,7 +658,7 @@ $(WORKSPACE)/work/%/restart/ocean.stats: build/symmetric/MOM6 | preproc && halfperiod=$$(awk -v t=$${daymax} -v dt=$${timeunit} 'BEGIN {printf "%.f", 0.5*t*dt}') \ && printf "\n&ocean_solo_nml\n seconds = $${halfperiod}\n/\n" >> input.nml # Remove any previous archived output - rm -f $(WORKSPACE)/results/$*/std.restart{1,2}.{out,err} + rm -f $(WORK)/results/$*/std.restart{1,2}.{out,err} # Run the first half-period cd $(@D) && $(TIME) $(MPIRUN) -n 1 $(abspath $<) 2> std1.err > std1.out \ || !( \ @@ -660,7 +684,7 @@ $(WORKSPACE)/work/%/restart/ocean.stats: build/symmetric/MOM6 | preproc # Not a true rule; only call this after `make test` to summarize test results. .PHONY: test.summary test.summary: - @./tools/report_test_results.sh $(WORKSPACE)/results + @./tools/report_test_results.sh $(WORK)/results #--- @@ -668,48 +692,50 @@ test.summary: # NOTE: Using file parser gcov report as a proxy for test completion .PHONY: run.cov.unit -run.cov.unit: build/unit/MOM_file_parser_tests.F90.gcov +run.cov.unit: $(BUILD)/unit/MOM_file_parser_tests.F90.gcov .PHONY: build.unit -build.unit: $(foreach f, $(UNIT_EXECS), build/unit/$(f)) +build.unit: $(foreach f, $(UNIT_EXECS), $(BUILD)/unit/$(f)) .PHONY: run.unit run.unit: $(foreach f, $(UNIT_EXECS), work/unit/$(f).out) .PHONY: build.timing -build.timing: $(foreach f, $(TIMING_EXECS), build/timing/$(f)) +build.timing: $(foreach f, $(TIMING_EXECS), $(BUILD)/timing/$(f)) .PHONY: run.timing run.timing: $(foreach f, $(TIMING_EXECS), work/timing/$(f).out) .PHONY: show.timing show.timing: $(foreach f, $(TIMING_EXECS), work/timing/$(f).show) -$(WORKSPACE)/work/timing/%.show: +$(WORK)/timing/%.show: ./tools/disp_timing.py $(@:.show=.out) + # Invoke the above unit/timing rules for a "target" code # Invoke with appropriate macros defines, i.e. -# make build.timing_target MOM_TARGET_URL=... MOM_TARGET_BRANCH=... TARGET_CODEBASE=build/target_codebase -# make run.timing_target TARGET_CODEBASE=build/target_codebase +# make build.timing_target MOM_TARGET_URL=... MOM_TARGET_BRANCH=... TARGET_CODEBASE=$(BUILD)/target_codebase +# make run.timing_target TARGET_CODEBASE=$(BUILD)/target_codebase TIMING_TARGET_EXECS ?= $(basename $(notdir $(wildcard $(TARGET_CODEBASE)/config_src/drivers/timing_tests/*.F90) ) ) .PHONY: build.timing_target -build.timing_target: $(foreach f, $(TIMING_TARGET_EXECS), $(TARGET_CODEBASE)/.testing/build/timing/$(f)) +build.timing_target: $(foreach f, $(TIMING_TARGET_EXECS), $(TARGET_CODEBASE)/.testing/$(BUILD)/timing/$(f)) .PHONY: run.timing_target run.timing_target: $(foreach f, $(TIMING_TARGET_EXECS), $(TARGET_CODEBASE)/.testing/work/timing/$(f).out) .PHONY: compare.timing compare.timing: $(foreach f, $(filter $(TIMING_EXECS),$(TIMING_TARGET_EXECS)), work/timing/$(f).compare) -$(WORKSPACE)/work/timing/%.compare: $(TARGET_CODEBASE) +$(WORK)/timing/%.compare: $(TARGET_CODEBASE) ./tools/disp_timing.py -r $(TARGET_CODEBASE)/.testing/$(@:.compare=.out) $(@:.compare=.out) $(TARGET_CODEBASE)/.testing/%: | $(TARGET_CODEBASE) cd $(TARGET_CODEBASE)/.testing && make $* + # General rule to run a unit test executable -# Pattern is to run build/unit/executable and direct output to executable.out -$(WORKSPACE)/work/unit/%.out: build/unit/% +# Pattern is to run $(BUILD)/unit/executable and direct output to executable.out +$(WORK)/unit/%.out: $(BUILD)/unit/% @mkdir -p $(@D) cd $(@D) ; $(TIME) $(MPIRUN) -n 1 $(abspath $<) 2> >(tee $*.err) > $*.out -$(WORKSPACE)/work/unit/test_MOM_file_parser.out: build/unit/test_MOM_file_parser +$(WORK)/unit/test_MOM_file_parser.out: $(BUILD)/unit/test_MOM_file_parser if [ $(REPORT_COVERAGE) ]; then \ - find build/unit -name *.gcda -exec rm -f '{}' \; ; \ + find $(BUILD)/unit -name *.gcda -exec rm -f '{}' \; ; \ fi mkdir -p $(@D) cd $(@D) \ @@ -727,31 +753,31 @@ $(WORKSPACE)/work/unit/test_MOM_file_parser.out: build/unit/test_MOM_file_parser ) # NOTE: .gcov actually depends on .gcda, but .gcda is produced with std.out -# TODO: Replace $(WORKSPACE)/work/unit/std.out with *.gcda? -build/unit/MOM_file_parser_tests.F90.gcov: $(WORKSPACE)/work/unit/test_MOM_file_parser.out +# TODO: Replace $(WORK)/unit/std.out with *.gcda? +$(BUILD)/unit/MOM_file_parser_tests.F90.gcov: $(WORK)/unit/test_MOM_file_parser.out cd $(@D) \ && gcov -b *.gcda > gcov.unit.out find $(@D) -name "*.gcov" -exec sed -i -r 's/^( *[0-9]*)\*:/ \1:/g' {} \; .PHONY: report.cov.unit -report.cov.unit: build/unit/MOM_file_parser_tests.F90.gcov codecov - ./codecov $(CODECOV_TOKEN_ARG) -R build/unit -f "*.gcov" -Z -n "Unit tests" \ - > build/unit/codecov.out \ - 2> build/unit/codecov.err \ +report.cov.unit: $(BUILD)/unit/MOM_file_parser_tests.F90.gcov codecov + ./codecov $(CODECOV_TOKEN_ARG) -R $(BUILD)/unit -f "*.gcov" -Z -n "Unit tests" \ + > $(BUILD)/unit/codecov.out \ + 2> $(BUILD)/unit/codecov.err \ && echo -e "${MAGENTA}Report uploaded to codecov.${RESET}" \ || { \ - cat build/unit/codecov.err ; \ + cat $(BUILD)/unit/codecov.err ; \ echo -e "${RED}Failed to upload report.${RESET}" ; \ if [ "$(REQUIRE_COVERAGE_UPLOAD)" = true ] ; then false ; fi ; \ } -$(WORKSPACE)/work/timing/%.out: build/timing/% FORCE +$(WORK)/timing/%.out: $(BUILD)/timing/% FORCE @mkdir -p $(@D) @echo Running $< in $(@D) @cd $(@D) ; $(TIME) $(MPIRUN) -n 1 $(abspath $<) 2> $*.err > $*.out -#--- -# Profiling based on FMS clocks + +## Profiling based on FMS clocks PCONFIGS = p0 @@ -759,17 +785,17 @@ PCONFIGS = p0 profile: $(foreach p,$(PCONFIGS), prof.$(p)) .PHONY: prof.p0 -prof.p0: $(WORKSPACE)/work/p0/opt/clocks.json $(WORKSPACE)/work/p0/opt_target/clocks.json +prof.p0: $(WORK)/p0/opt/clocks.json $(WORK)/p0/opt_target/clocks.json python tools/compare_clocks.py $^ -$(WORKSPACE)/work/p0/%/clocks.json: $(WORKSPACE)/work/p0/%/std.out +$(WORK)/p0/%/clocks.json: $(WORK)/p0/%/std.out python tools/parse_fms_clocks.py -d $(@D) $^ > $@ \ || !( rm $@ ) -$(WORKSPACE)/work/p0/opt/std.out: build/opt/MOM6 -$(WORKSPACE)/work/p0/opt_target/std.out: build/opt_target/MOM6 +$(WORK)/p0/opt/std.out: $(BUILD)/opt/MOM6 +$(WORK)/p0/opt_target/std.out: $(BUILD)/opt_target/MOM6 -$(WORKSPACE)/work/p0/%/std.out: +$(WORK)/p0/%/std.out: mkdir -p $(@D) cp -RL p0/* $(@D) mkdir -p $(@D)/RESTART @@ -778,8 +804,7 @@ $(WORKSPACE)/work/p0/%/std.out: && $(MPIRUN) -n 1 $(abspath $<) 2> std.err > std.out -#--- -# Profiling based on perf output +## Profiling based on perf output # TODO: This expects the -e flag, can I handle it in the command? PERF_EVENTS ?= @@ -788,16 +813,16 @@ PERF_EVENTS ?= perf: $(foreach p,$(PCONFIGS), perf.$(p)) .PHONY: prof.p0 -perf.p0: $(WORKSPACE)/work/p0/opt/profile.json $(WORKSPACE)/work/p0/opt_target/profile.json +perf.p0: $(WORK)/p0/opt/profile.json $(WORK)/p0/opt_target/profile.json python tools/compare_perf.py $^ -$(WORKSPACE)/work/p0/%/profile.json: $(WORKSPACE)/work/p0/%/perf.data +$(WORK)/p0/%/profile.json: $(WORK)/p0/%/perf.data python tools/parse_perf.py -f $< > $@ -$(WORKSPACE)/work/p0/opt/perf.data: build/opt/MOM6 -$(WORKSPACE)/work/p0/opt_target/perf.data: build/opt_target/MOM6 +$(WORK)/p0/opt/perf.data: $(BUILD)/opt/MOM6 +$(WORK)/p0/opt_target/perf.data: $(BUILD)/opt_target/MOM6 -$(WORKSPACE)/work/p0/%/perf.data: +$(WORK)/p0/%/perf.data: mkdir -p $(@D) cp -RL p0/* $(@D) mkdir -p $(@D)/RESTART @@ -810,25 +835,26 @@ $(WORKSPACE)/work/p0/%/perf.data: || cat std.perf.err -#---- +## Cleanup # NOTE: These tests assert that we are in the .testing directory. .PHONY: clean clean: clean.build clean.stats - @[ $$(basename $$(pwd)) = .testing ] - rm -rf deps + rm -rf $(BUILD) .PHONY: clean.build clean.build: @[ $$(basename $$(pwd)) = .testing ] - rm -rf build + for b in $(ALL_EXECS); do \ + rm -rf $(BUILD)/$${b}; \ + done .PHONY: clean.stats clean.stats: @[ $$(basename $$(pwd)) = .testing ] - rm -rf $(WORKSPACE)/work $(WORKSPACE)/results + rm -rf $(WORK) .PHONY: clean.preproc diff --git a/ac/configure.ac b/ac/configure.ac index 9d87240506..c774c39129 100644 --- a/ac/configure.ac +++ b/ac/configure.ac @@ -82,13 +82,13 @@ AS_IF([test "x$with_driver" != "x"], # Select the model framework (default: FMS1) # NOTE: We can phase this out after the FMS1 I/O has been removed from FMS and # replace with a detection test. For now, it is a user-defined switch. -MODEL_FRAMEWORK=${srcdir}/config_src/infra/FMS1 +MODEL_FRAMEWORK=${srcdir}/config_src/infra/FMS2 AC_ARG_WITH([framework], AS_HELP_STRING([--with-framework=fms1|fms2], [Select the model framework])) AS_CASE(["$with_framework"], [fms1], [MODEL_FRAMEWORK=${srcdir}/config_src/infra/FMS1], [fms2], [MODEL_FRAMEWORK=${srcdir}/config_src/infra/FMS2], - [MODEL_FRAMEWORK=${srcdir}/config_src/infra/FMS1] + [MODEL_FRAMEWORK=${srcdir}/config_src/infra/FMS2] ) From 2ae638504b5eb9c1d492a18fa8c4ab471f17153a Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 22 Dec 2023 09:34:36 -0500 Subject: [PATCH 05/15] Refactor set_viscous_BBL channel drag calculations Refactored set_viscous_BBL to separate out the routines setting the open interface lengths used for the channel drag, shortening a 1070 line long routine to 915 lines and reducing the scope of a number of temporary variables. A number of logical branch points have been moved outside of the innermost do loops. This refactoring will also make it easier to provide alternatives to some of the solvers that do not use the trigonometric functions to solve for the roots of a cubic expression and avoiding the issues noted at NOAA-GFDL/MOM6/issues/483. All answers are bitwise identical and public interfaces are unchanged. --- .../vertical/MOM_set_viscosity.F90 | 511 +++++++++++------- 1 file changed, 322 insertions(+), 189 deletions(-) diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index b207b1ff1c..bf0db51a44 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -20,6 +20,7 @@ module MOM_set_visc use MOM_grid, only : ocean_grid_type use MOM_hor_index, only : hor_index_type use MOM_interface_heights, only : thickness_to_dz +! use MOM_intrinsic_functions, only : cuberoot use MOM_io, only : slasher, MOM_read_data, vardesc, var_desc use MOM_kappa_shear, only : kappa_shear_is_used, kappa_shear_at_vertex use MOM_open_boundary, only : ocean_OBC_type, OBC_segment_type, OBC_NONE, OBC_DIRECTION_E @@ -258,40 +259,15 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) real :: crv ! crv is the curvature of the bottom depth across a ! cell, times the cell width squared [Z ~> m]. real :: crv_3 ! crv/3 [Z ~> m]. - real :: C24_crv ! 24/crv [Z-1 ~> m-1]. real :: slope ! The absolute value of the bottom depth slope across ! a cell times the cell width [Z ~> m]. - real :: apb_4a, ax2_3apb ! Various nondimensional ratios of crv and slope [nondim]. - real :: a2x48_apb3, Iapb, Ibma_2 ! Combinations of crv (a) and slope (b) [Z-1 ~> m-1] - ! All of the following "volumes" have units of vertical heights because they are normalized - ! by the full horizontal area of a velocity cell. real :: Vol_bbl_chan ! The volume of the bottom boundary layer as used in the channel ! drag parameterization, normalized by the full horizontal area ! of the velocity cell [Z ~> m]. - real :: Vol_open ! The cell volume above which it is open [Z ~> m]. - real :: Vol_direct ! With less than Vol_direct [Z ~> m], there is a direct - ! solution of a cubic equation for L. - real :: Vol_2_reg ! The cell volume above which there are two separate - ! open areas that must be integrated [Z ~> m]. - real :: vol ! The volume below the interface whose normalized - ! width is being sought [Z ~> m]. - real :: vol_below ! The volume below the interface below the one that - ! is currently under consideration [Z ~> m]. - real :: Vol_err ! The error in the volume with the latest estimate of - ! L, or the error for the interface below [Z ~> m]. - real :: Vol_quit ! The volume error below which to quit iterating [Z ~> m]. - real :: Vol_tol ! A volume error tolerance [Z ~> m]. + real :: vol_below(SZK_(GV)+1) ! The volume below each interface, normalized by the full + ! horizontal area of a velocity cell [Z ~> m]. real :: L(SZK_(GV)+1) ! The fraction of the full cell width that is open at ! the depth of each interface [nondim]. - real :: L_direct ! The value of L above volume Vol_direct [nondim]. - real :: L_max, L_min ! Upper and lower bounds on the correct value for L [nondim]. - real :: Vol_err_max ! The volume error for the upper bound on the correct value for L [Z ~> m] - real :: Vol_err_min ! The volume error for the lower bound on the correct value for L [Z ~> m] - real :: Vol_0 ! A deeper volume with known width L0 [Z ~> m]. - real :: L0 ! The value of L above volume Vol_0 [nondim]. - real :: dVol ! vol - Vol_0 [Z ~> m]. - real :: dV_dL2 ! The partial derivative of volume with L squared - ! evaluated at L=L0 [Z ~> m]. real :: h_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected [H ~> m or kg m-2]. real :: dz_neglect ! A vertical distance that is so small it is usually lost @@ -315,14 +291,9 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) real, parameter :: C1_3 = 1.0/3.0, C1_6 = 1.0/6.0, C1_12 = 1.0/12.0 ! Rational constants [nondim] real :: C2pi_3 ! An irrational constant, 2/3 pi. [nondim] real :: tmp ! A temporary variable, sometimes in [Z ~> m] - real :: tmp_val_m1_to_p1 ! A temporary variable [nondim] - real :: curv_tol ! Numerator of curvature cubed, used to estimate - ! accuracy of a single L(:) Newton iteration [Z5 ~> m5] - logical :: use_L0, do_one_L_iter ! Control flags for L(:) Newton iteration logical :: use_BBL_EOS, do_i(SZIB_(G)) integer, dimension(2) :: EOSdom ! The computational domain for the equation of state integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, m, n, K2, nkmb, nkml - integer :: itt, maxitt=20 type(ocean_OBC_type), pointer :: OBC => NULL() is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke @@ -332,7 +303,6 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) dz_neglect = GV%dZ_subroundoff Rho0x400_G = 400.0*(GV%H_to_RZ / (US%L_to_Z**2 * GV%g_Earth)) - Vol_quit = (0.9*GV%Angstrom_Z + dz_neglect) C2pi_3 = 8.0*atan(1.0)/3.0 if (.not.CS%initialized) call MOM_error(FATAL,"MOM_set_viscosity(BBL): "//& @@ -442,8 +412,7 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) !$OMP Isq,Ieq,Jsq,Jeq,h_neglect,dz_neglect,Rho0x400_G,C2pi_3, & !$OMP U_bg_sq,cdrag_sqrt,cdrag_sqrt_H,cdrag_sqrt_H_RL, & !$OMP cdrag_L_to_H,cdrag_RL_to_H,use_BBL_EOS,BBL_thick_max, & - !$OMP OBC,maxitt,D_u,D_v,mask_u,mask_v,pbv) & - !$OMP firstprivate(Vol_quit) + !$OMP OBC,D_u,D_v,mask_u,mask_v,pbv) do j=Jsq,Jeq ; do m=1,2 if (m==1) then @@ -865,12 +834,11 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) if (CS%body_force_drag) bbl_thick = dz_bbl_drag(i) if (CS%Channel_drag) then - ! The drag within the bottommost Vol_bbl_chan is applied as a part of - ! an enhanced bottom viscosity, while above this the drag is applied - ! directly to the layers in question as a Rayleigh drag term. - ! Restrict the volume over which the channel drag is applied. - if (CS%Chan_drag_max_vol >= 0.0) Vol_bbl_chan = min(Vol_bbl_chan, CS%Chan_drag_max_vol) + vol_below(nz+1) = 0.0 + do K=nz,1,-1 + vol_below(K) = vol_below(K+1) + dz_vel(i,k) + enddo !### The harmonic mean edge depths here are not invariant to offsets! if (m==1) then @@ -887,162 +855,33 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) Dm = 2.0 * D_vel * tmp / (D_vel + tmp) endif if (Dm > Dp) then ; tmp = Dp ; Dp = Dm ; Dm = tmp ; endif - crv_3 = (Dp + Dm - 2.0*D_vel) ; crv = 3.0*crv_3 slope = Dp - Dm + ! If the curvature is small enough, there is no reason not to assume ! a uniformly sloping or flat bottom. if (abs(crv) < 1e-2*(slope + CS%BBL_thick_min)) crv = 0.0 - ! Each cell extends from x=-1/2 to 1/2, and has a topography - ! given by D(x) = crv*x^2 + slope*x + D - crv/12. - - ! Calculate the volume above which the entire cell is open and the - ! other volumes at which the equation that is solved for L changes. - if (crv > 0.0) then - if (slope >= crv) then - Vol_open = D_vel - Dm ; Vol_2_reg = Vol_open - else - tmp = slope/crv - Vol_open = 0.25*slope*tmp + C1_12*crv - Vol_2_reg = 0.5*tmp**2 * (crv - C1_3*slope) - endif - ! Define some combinations of crv & slope for later use. - C24_crv = 24.0/crv ; Iapb = 1.0/(crv+slope) - apb_4a = (slope+crv)/(4.0*crv) ; a2x48_apb3 = (48.0*(crv*crv))*(Iapb**3) - ax2_3apb = 2.0*C1_3*crv*Iapb - elseif (crv == 0.0) then - Vol_open = 0.5*slope - if (slope > 0) Iapb = 1.0/slope - else ! crv < 0.0 - Vol_open = D_vel - Dm - if (slope >= -crv) then - Iapb = 1.0e30*US%Z_to_m ; if (slope+crv /= 0.0) Iapb = 1.0/(crv+slope) - Vol_direct = 0.0 ; L_direct = 0.0 ; C24_crv = 0.0 - else - C24_crv = 24.0/crv ; Iapb = 1.0/(crv+slope) - L_direct = 1.0 + slope/crv ! L_direct < 1 because crv < 0 - Vol_direct = -C1_6*crv*L_direct**3 - endif - Ibma_2 = 2.0 / (slope - crv) - endif - L(nz+1) = 0.0 ; vol = 0.0 ; Vol_err = 0.0 ; BBL_visc_frac = 0.0 - ! Determine the normalized open length at each interface. - do K=nz,1,-1 - vol_below = vol - - vol = vol + dz_vel(i,k) - h_vel_pos = h_vel(i,k) + h_neglect - - if (vol >= Vol_open) then ; L(K) = 1.0 - elseif (crv == 0) then ! The bottom has no curvature. - L(K) = sqrt(2.0*vol*Iapb) - elseif (crv > 0) then - ! There may be a minimum depth, and there are - ! analytic expressions for L for all cases. - if (vol < Vol_2_reg) then - ! In this case, there is a contiguous open region and - ! vol = 0.5*L^2*(slope + crv/3*(3-4L)). - if (a2x48_apb3*vol < 1e-8) then ! Could be 1e-7? - ! There is a very good approximation here for massless layers. - L0 = sqrt(2.0*vol*Iapb) ; L(K) = L0*(1.0 + ax2_3apb*L0) - else - L(K) = apb_4a * (1.0 - & - 2.0 * cos(C1_3*acos(a2x48_apb3*vol - 1.0) - C2pi_3)) - endif - ! To check the answers. - ! Vol_err = 0.5*(L(K)*L(K))*(slope + crv_3*(3.0-4.0*L(K))) - vol - else ! There are two separate open regions. - ! vol = slope^2/4crv + crv/12 - (crv/12)*(1-L)^2*(1+2L) - ! At the deepest volume, L = slope/crv, at the top L = 1. - !L(K) = 0.5 - cos(C1_3*acos(1.0 - C24_crv*(Vol_open - vol)) - C2pi_3) - tmp_val_m1_to_p1 = 1.0 - C24_crv*(Vol_open - vol) - tmp_val_m1_to_p1 = max(-1., min(1., tmp_val_m1_to_p1)) - L(K) = 0.5 - cos(C1_3*acos(tmp_val_m1_to_p1) - C2pi_3) - ! To check the answers. - ! Vol_err = Vol_open - 0.25*crv_3*(1.0+2.0*L(K)) * (1.0-L(K))**2 - vol - endif - else ! a < 0. - if (vol <= Vol_direct) then - ! Both edges of the cell are bounded by walls. - L(K) = (-0.25*C24_crv*vol)**C1_3 - else - ! x_R is at 1/2 but x_L is in the interior & L is found by solving - ! vol = 0.5*L^2*(slope + crv/3*(3-4L)) - - ! Vol_err = 0.5*(L(K+1)*L(K+1))*(slope + crv_3*(3.0-4.0*L(K+1))) - vol_below - ! Change to ... - ! if (min(vol_below + Vol_err, vol) <= Vol_direct) then ? - if (vol_below + Vol_err <= Vol_direct) then - L0 = L_direct ; Vol_0 = Vol_direct - else - L0 = L(K+1) ; Vol_0 = vol_below + Vol_err - ! Change to Vol_0 = min(vol_below + Vol_err, vol) ? - endif - - ! Try a relatively simple solution that usually works well - ! for massless layers. - dV_dL2 = 0.5*(slope+crv) - crv*L0 ; dVol = (vol-Vol_0) - ! dV_dL2 = 0.5*(slope+crv) - crv*L0 ; dVol = max(vol-Vol_0, 0.0) - - use_L0 = .false. - do_one_L_iter = .false. - if (CS%answer_date < 20190101) then - curv_tol = GV%Angstrom_Z*dV_dL2**2 & - * (0.25 * dV_dL2 * GV%Angstrom_Z - crv * L0 * dVol) - do_one_L_iter = (crv * crv * dVol**3) < curv_tol - else - ! The following code is more robust when GV%Angstrom_H=0, but - ! it changes answers. - use_L0 = (dVol <= 0.) - - Vol_tol = max(0.5 * GV%Angstrom_Z + dz_neglect, 1e-14 * vol) - Vol_quit = max(0.9 * GV%Angstrom_Z + dz_neglect, 1e-14 * vol) + ! Determine the normalized open length (L) at each interface. + if (crv == 0.0) then + call find_L_open_uniform_slope(vol_below, Dp, Dm, L, GV) + elseif (crv > 0.0) then + call find_L_open_concave_analytic(vol_below, D_vel, Dp, Dm, C2pi_3, L, GV) + else ! crv < 0.0 + call find_L_open_convex(vol_below, D_vel, Dp, Dm, L, GV, US, CS) + endif ! end of crv<0 cases. - curv_tol = Vol_tol * dV_dL2**2 & - * (dV_dL2 * Vol_tol - 2.0 * crv * L0 * dVol) - do_one_L_iter = (crv * crv * dVol**3) < curv_tol - endif + ! Determine the Rayliegh drag contributions. - if (use_L0) then - L(K) = L0 - Vol_err = 0.5*(L(K)*L(K))*(slope + crv_3*(3.0-4.0*L(K))) - vol - elseif (do_one_L_iter) then - ! One iteration of Newton's method should give an estimate - ! that is accurate to within Vol_tol. - L(K) = sqrt(L0*L0 + dVol / dV_dL2) - Vol_err = 0.5*(L(K)*L(K))*(slope + crv_3*(3.0-4.0*L(K))) - vol - else - if (dV_dL2*(1.0-L0*L0) < dVol + & - dV_dL2 * (Vol_open - Vol)*Ibma_2) then - L_max = sqrt(1.0 - (Vol_open - Vol)*Ibma_2) - else - L_max = sqrt(L0*L0 + dVol / dV_dL2) - endif - L_min = sqrt(L0*L0 + dVol / (0.5*(slope+crv) - crv*L_max)) + ! The drag within the bottommost Vol_bbl_chan is applied as a part of an enhanced bottom + ! viscosity, while above this the drag is applied directly to the layers in question as a + ! Rayleigh drag term. - Vol_err_min = 0.5*(L_min**2)*(slope + crv_3*(3.0-4.0*L_min)) - vol - Vol_err_max = 0.5*(L_max**2)*(slope + crv_3*(3.0-4.0*L_max)) - vol - ! if ((abs(Vol_err_min) <= Vol_quit) .or. (Vol_err_min >= Vol_err_max)) then - if (abs(Vol_err_min) <= Vol_quit) then - L(K) = L_min ; Vol_err = Vol_err_min - else - L(K) = sqrt((L_min**2*Vol_err_max - L_max**2*Vol_err_min) / & - (Vol_err_max - Vol_err_min)) - do itt=1,maxitt - Vol_err = 0.5*(L(K)*L(K))*(slope + crv_3*(3.0-4.0*L(K))) - vol - if (abs(Vol_err) <= Vol_quit) exit - ! Take a Newton's method iteration. This equation has proven - ! robust enough not to need bracketing. - L(K) = L(K) - Vol_err / (L(K)* (slope + crv - 2.0*crv*L(K))) - ! This would be a Newton's method iteration for L^2: - ! L(K) = sqrt(L(K)*L(K) - Vol_err / (0.5*(slope+crv) - crv*L(K))) - enddo - endif ! end of iterative solver - endif ! end of 1-boundary alternatives. - endif ! end of a<0 cases. - endif + ! Restrict the volume over which the channel drag is applied from the previously determined valu.e + if (CS%Chan_drag_max_vol >= 0.0) Vol_bbl_chan = min(Vol_bbl_chan, CS%Chan_drag_max_vol) + BBL_visc_frac = 0.0 + do K=nz,1,-1 !modify L(K) for porous barrier parameterization if (m==1) then ; L(K) = L(K)*pbv%por_layer_widthU(I,j,K) else ; L(K) = L(K)*pbv%por_layer_widthV(i,J,K); endif @@ -1050,8 +889,8 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) ! Determine the drag contributing to the bottom boundary layer ! and the Rayleigh drag that acts on each layer. if (L(K) > L(K+1)) then - if (vol_below < Vol_bbl_chan) then - BBL_frac = (1.0-vol_below/Vol_bbl_chan)**2 + if (vol_below(K+1) < Vol_bbl_chan) then + BBL_frac = (1.0-vol_below(K+1)/Vol_bbl_chan)**2 BBL_visc_frac = BBL_visc_frac + BBL_frac*(L(K) - L(K+1)) else BBL_frac = 0.0 @@ -1063,6 +902,7 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) cdrag_conv = cdrag_L_to_H endif + h_vel_pos = h_vel(i,k) + h_neglect if (m==1) then ; Cell_width = G%dy_Cu(I,j)*pbv%por_face_areaU(I,j,k) else ; Cell_width = G%dx_Cv(i,J)*pbv%por_face_areaV(i,J,k) ; endif gam = 1.0 - L(K+1)/L(K) @@ -1085,7 +925,7 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) else ; visc%Ray_v(i,J,k) = 0.0 ; endif endif - enddo ! k loop to determine L(K). + enddo ! k loop to determine visc%Ray_[uv]. ! Set the near-bottom viscosity to a value which will give ! the correct stress when the shear occurs over bbl_thick. @@ -1202,6 +1042,299 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) end subroutine set_viscous_BBL +!> Determine the normalized open length of each interface, given the edge depths and normalized +!! volumes below each interface. +subroutine find_L_open_uniform_slope(vol_below, Dp, Dm, L, GV) + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + real, dimension(SZK_(GV)+1), intent(in) :: vol_below !< The volume below each interface, normalized by + !! the full horizontal area of a velocity cell [Z ~> m] + real, intent(in) :: Dp !< The larger of the two depths at the edge + !! of a velocity cell [Z ~> m] + real, intent(in) :: Dm !< The smaller of the two depths at the edge + !! of a velocity cell [Z ~> m] + real, dimension(SZK_(GV)+1), intent(out) :: L !< The fraction of the full cell width that is open at + !! the depth of each interface [nondim] + + ! Local variables + real :: slope ! The absolute value of the bottom depth slope across a cell times the cell width [Z ~> m]. + real :: I_slope ! The inverse of the normalized slope [Z-1 ~> m-1] + real :: Vol_open ! The cell volume above which it is open [Z ~> m]. + integer :: K, nz + + nz = GV%ke + + slope = abs(Dp - Dm) + if (slope == 0.0) then + L(1:nz) = 1.0 ; L(nz+1) = 0.0 + else + Vol_open = 0.5*slope + I_slope = 1.0 / slope + + L(nz+1) = 0.0 + do K=nz,1,-1 + if (vol_below(K) >= Vol_open) then ; L(K) = 1.0 + else + ! With a uniformly sloping bottom, the calculation of L(K) is the solution of a simple quadratic equation. + L(K) = sqrt(2.0*vol_below(K)*I_slope) + endif + enddo + endif + +end subroutine find_L_open_uniform_slope + +!> Determine the normalized open length of each interface for concave bathymetry (from the ocean perspective) using +!! analytic expressions. In this case there can be two separate open regions. +subroutine find_L_open_concave_analytic(vol_below, D_vel, Dp, Dm, C2pi_3, L, GV) + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + real, dimension(SZK_(GV)+1), intent(in) :: vol_below !< The volume below each interface, normalized by + !! the full horizontal area of a velocity cell [Z ~> m] + real, intent(in) :: D_vel !< The average bottom depth at a velocity point [Z ~> m] + real, intent(in) :: Dp !< The larger of the two depths at the edge + !! of a velocity cell [Z ~> m] + real, intent(in) :: Dm !< The smaller of the two depths at the edge + !! of a velocity cell [Z ~> m] + real, intent(in) :: C2pi_3 !< An irrational constant, 2/3 pi [nondim] + real, dimension(SZK_(GV)+1), intent(out) :: L !< The fraction of the full cell width that is open at + !! the depth of each interface [nondim] + + ! Local variables + real :: crv ! crv is the curvature of the bottom depth across a + ! cell, times the cell width squared [Z ~> m]. + real :: crv_3 ! crv/3 [Z ~> m]. + real :: slope ! The absolute value of the bottom depth slope across + ! a cell times the cell width [Z ~> m]. + ! The following "volumes" have units of vertical heights because they are normalized + ! by the full horizontal area of a velocity cell. + real :: Vol_open ! The cell volume above which the face is fully is open [Z ~> m]. + real :: Vol_2_reg ! The cell volume above which there are two separate + ! open areas that must be integrated [Z ~> m]. + real :: C24_crv ! 24/crv [Z-1 ~> m-1]. + real :: apb_4a, ax2_3apb ! Various nondimensional ratios of crv and slope [nondim]. + real :: a2x48_apb3, Iapb ! Combinations of crv (a) and slope (b) [Z-1 ~> m-1] + real :: L0 ! A linear estimate of L approprate for tiny volumes [nondim]. + real :: slope_crv ! The slope divided by the curvature [nondim] + real :: tmp_val_m1_to_p1 ! A temporary variable [nondim] + real, parameter :: C1_3 = 1.0/3.0, C1_12 = 1.0/12.0 ! Rational constants [nondim] + integer :: K, nz + + nz = GV%ke + + ! Each cell extends from x=-1/2 to 1/2, and has a topography + ! given by D(x) = crv*x^2 + slope*x + D_vel - crv/12. + crv_3 = (Dp + Dm - 2.0*D_vel) ; crv = 3.0*crv_3 + if (crv < 0.0) call MOM_error(FATAL, "find_L_open_concave should only be called with a positive curvature.") + slope = Dp - Dm + + ! Calculate the volume above which the entire cell is open and the volume at which the + ! equation that is solved for L changes because there are two separate open regions. + if (slope >= crv) then + Vol_open = D_vel - Dm ; Vol_2_reg = Vol_open + else + slope_crv = slope / crv + Vol_open = 0.25*slope*slope_crv + C1_12*crv + Vol_2_reg = 0.5*slope_crv**2 * (crv - C1_3*slope) + endif + ! Define some combinations of crv & slope for later use. + C24_crv = 24.0/crv ; Iapb = 1.0/(crv+slope) + apb_4a = (slope+crv)/(4.0*crv) ; a2x48_apb3 = (48.0*(crv*crv))*(Iapb**3) + ax2_3apb = 2.0*C1_3*crv*Iapb + + L(nz+1) = 0.0 + ! Determine the normalized open length (L) at each interface. + do K=nz,1,-1 + if (vol_below(K) >= Vol_open) then ! The whole cell is open. + L(K) = 1.0 + elseif (vol_below(K) < Vol_2_reg) then + ! In this case, there is a contiguous open region and + ! vol_below(K) = 0.5*L^2*(slope + crv/3*(3-4L)). + if (a2x48_apb3*vol_below(K) < 1e-8) then ! Could be 1e-7? + ! There is a very good approximation here for massless layers. + L0 = sqrt(2.0*vol_below(K)*Iapb) ; L(K) = L0*(1.0 + ax2_3apb*L0) + else + L(K) = apb_4a * (1.0 - & + 2.0 * cos(C1_3*acos(a2x48_apb3*vol_below(K) - 1.0) - C2pi_3)) + endif + ! To check the answers. + ! Vol_err = 0.5*(L(K)*L(K))*(slope + crv_3*(3.0-4.0*L(K))) - vol_below(K) + else ! There are two separate open regions. + ! vol_below(K) = slope^2/4crv + crv/12 - (crv/12)*(1-L)^2*(1+2L) + ! At the deepest volume, L = slope/crv, at the top L = 1. + ! L(K) = 0.5 - cos(C1_3*acos(1.0 - C24_crv*(Vol_open - vol_below(K))) - C2pi_3) + tmp_val_m1_to_p1 = 1.0 - C24_crv*(Vol_open - vol_below(K)) + tmp_val_m1_to_p1 = max(-1., min(1., tmp_val_m1_to_p1)) + L(K) = 0.5 - cos(C1_3*acos(tmp_val_m1_to_p1) - C2pi_3) + ! To check the answers. + ! Vol_err = Vol_open - 0.25*crv_3*(1.0+2.0*L(K)) * (1.0-L(K))**2 - vol_below(K) + endif + enddo ! k loop to determine L(K) in the concave case + +end subroutine find_L_open_concave_analytic + + +!> Determine the normalized open length of each interface for convex bathymetry (from the ocean perspective) using +!! Newton's method iterations. In this case there is a single open region with the minimum depth +!! at one edge of the cell. +subroutine find_L_open_convex(vol_below, D_vel, Dp, Dm, L, GV, US, CS) + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + real, dimension(SZK_(GV)+1), intent(in) :: vol_below !< The volume below each interface, normalized by + !! the full horizontal area of a velocity cell [Z ~> m] + real, intent(in) :: D_vel !< The average bottom depth at a velocity point [Z ~> m] + real, intent(in) :: Dp !< The larger of the two depths at the edge + !! of a velocity cell [Z ~> m] + real, intent(in) :: Dm !< The smaller of the two depths at the edge + !! of a velocity cell [Z ~> m] + real, dimension(SZK_(GV)+1), intent(out) :: L !< The fraction of the full cell width that is open at + !! the depth of each interface [nondim] + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(set_visc_CS), intent(in) :: CS !< The control structure returned by a previous + !! call to set_visc_init. + + ! Local variables + real :: crv ! crv is the curvature of the bottom depth across a + ! cell, times the cell width squared [Z ~> m]. + real :: crv_3 ! crv/3 [Z ~> m]. + real :: slope ! The absolute value of the bottom depth slope across + ! a cell times the cell width [Z ~> m]. + ! All of the following "volumes" have units of vertical heights because they are normalized + ! by the full horizontal area of a velocity cell. + real :: Vol_err ! The error in the volume with the latest estimate of + ! L, or the error for the interface below [Z ~> m]. + real :: Vol_quit ! The volume error below which to quit iterating [Z ~> m]. + real :: Vol_tol ! A volume error tolerance [Z ~> m]. + real :: Vol_open ! The cell volume above which the face is fully open [Z ~> m]. + real :: Vol_direct ! With less than Vol_direct [Z ~> m], there is a direct + ! solution of a cubic equation for L. + real :: Vol_err_max ! The volume error for the upper bound on the correct value for L [Z ~> m] + real :: Vol_err_min ! The volume error for the lower bound on the correct value for L [Z ~> m] + real :: Vol_0 ! A deeper volume with known width L0 [Z ~> m]. + real :: dVol ! vol - Vol_0 [Z ~> m]. + real :: dV_dL2 ! The partial derivative of volume with L squared + ! evaluated at L=L0 [Z ~> m]. + real :: L_direct ! The value of L above volume Vol_direct [nondim]. + real :: L_max, L_min ! Upper and lower bounds on the correct value for L [nondim]. + real :: L0 ! The value of L above volume Vol_0 [nondim]. + real :: Iapb, Ibma_2 ! Combinations of crv (a) and slope (b) [Z-1 ~> m-1] + real :: C24_crv ! 24/crv [Z-1 ~> m-1]. + real :: curv_tol ! Numerator of curvature cubed, used to estimate + ! accuracy of a single L(:) Newton iteration [Z5 ~> m5] + real, parameter :: C1_3 = 1.0/3.0, C1_6 = 1.0/6.0 ! Rational constants [nondim] + logical :: use_L0, do_one_L_iter ! Control flags for L(:) Newton iteration + integer :: K, nz, itt, maxitt=20 + + nz = GV%ke + + ! Each cell extends from x=-1/2 to 1/2, and has a topography + ! given by D(x) = crv*x^2 + slope*x + D_vel - crv/12. + crv_3 = (Dp + Dm - 2.0*D_vel) ; crv = 3.0*crv_3 + if (crv > 0.0) call MOM_error(FATAL, "find_L_open_convex should only be called with a negative curvature.") + slope = Dp - Dm + + ! Calculate the volume above which the entire cell is open and the volume at which the + ! equation that is solved for L changes because there is a direct solution. + Vol_open = D_vel - Dm + if (slope >= -crv) then + Iapb = 1.0e30*US%Z_to_m ; if (slope+crv /= 0.0) Iapb = 1.0/(crv+slope) + Vol_direct = 0.0 ; L_direct = 0.0 ; C24_crv = 0.0 + else + C24_crv = 24.0/crv ; Iapb = 1.0/(crv+slope) + L_direct = 1.0 + slope/crv ! L_direct < 1 because crv < 0 + Vol_direct = -C1_6*crv*L_direct**3 + endif + Ibma_2 = 2.0 / (slope - crv) + + if (CS%answer_date < 20190101) Vol_quit = (0.9*GV%Angstrom_Z + GV%dZ_subroundoff) + + L(nz+1) = 0.0 ; Vol_err = 0.0 + ! Determine the normalized open length (L) at each interface. + do K=nz,1,-1 + if (vol_below(K) >= Vol_open) then + L(K) = 1.0 + elseif (vol_below(K) <= Vol_direct) then + ! Both edges of the cell are bounded by walls. + ! if (CS%answer_date < 20240101)) then + L(K) = (-0.25*C24_crv*vol_below(K))**C1_3 + ! else + ! L(K) = cuberoot(-0.25*C24_crv*vol_below(K)) + ! endif + else + ! x_R is at 1/2 but x_L is in the interior & L is found by iteratively solving + ! vol_below(K) = 0.5*L^2*(slope + crv/3*(3-4L)) + + ! Vol_err = 0.5*(L(K+1)*L(K+1))*(slope + crv_3*(3.0-4.0*L(K+1))) - vol_below(K+1) + ! Change to ... + ! if (min(vol_below(K+1) + Vol_err, vol_below(K)) <= Vol_direct) then ? + if (vol_below(K+1) + Vol_err <= Vol_direct) then + L0 = L_direct ; Vol_0 = Vol_direct + else + L0 = L(K+1) ; Vol_0 = vol_below(K+1) + Vol_err + ! Change to Vol_0 = min(vol_below(K+1) + Vol_err, vol_below(K)) ? + endif + + ! Try a relatively simple solution that usually works well + ! for massless layers. + dV_dL2 = 0.5*(slope+crv) - crv*L0 ; dVol = (vol_below(K)-Vol_0) + ! dV_dL2 = 0.5*(slope+crv) - crv*L0 ; dVol = max(vol_below(K)-Vol_0, 0.0) + + use_L0 = .false. + do_one_L_iter = .false. + if (CS%answer_date < 20190101) then + curv_tol = GV%Angstrom_Z*dV_dL2**2 & + * (0.25 * dV_dL2 * GV%Angstrom_Z - crv * L0 * dVol) + do_one_L_iter = (crv * crv * dVol**3) < curv_tol + else + ! The following code is more robust when GV%Angstrom_H=0, but + ! it changes answers. + use_L0 = (dVol <= 0.) + + Vol_tol = max(0.5 * GV%Angstrom_Z + GV%dZ_subroundoff, 1e-14 * vol_below(K)) + Vol_quit = max(0.9 * GV%Angstrom_Z + GV%dZ_subroundoff, 1e-14 * vol_below(K)) + + curv_tol = Vol_tol * dV_dL2**2 & + * (dV_dL2 * Vol_tol - 2.0 * crv * L0 * dVol) + do_one_L_iter = (crv * crv * dVol**3) < curv_tol + endif + + if (use_L0) then + L(K) = L0 + Vol_err = 0.5*(L(K)*L(K))*(slope + crv_3*(3.0-4.0*L(K))) - vol_below(K) + elseif (do_one_L_iter) then + ! One iteration of Newton's method should give an estimate + ! that is accurate to within Vol_tol. + L(K) = sqrt(L0*L0 + dVol / dV_dL2) + Vol_err = 0.5*(L(K)*L(K))*(slope + crv_3*(3.0-4.0*L(K))) - vol_below(K) + else + if (dV_dL2*(1.0-L0*L0) < dVol + & + dV_dL2 * (Vol_open - vol_below(K))*Ibma_2) then + L_max = sqrt(1.0 - (Vol_open - vol_below(K))*Ibma_2) + else + L_max = sqrt(L0*L0 + dVol / dV_dL2) + endif + L_min = sqrt(L0*L0 + dVol / (0.5*(slope+crv) - crv*L_max)) + + Vol_err_min = 0.5*(L_min**2)*(slope + crv_3*(3.0-4.0*L_min)) - vol_below(K) + Vol_err_max = 0.5*(L_max**2)*(slope + crv_3*(3.0-4.0*L_max)) - vol_below(K) + ! if ((abs(Vol_err_min) <= Vol_quit) .or. (Vol_err_min >= Vol_err_max)) then + if (abs(Vol_err_min) <= Vol_quit) then + L(K) = L_min ; Vol_err = Vol_err_min + else + L(K) = sqrt((L_min**2*Vol_err_max - L_max**2*Vol_err_min) / & + (Vol_err_max - Vol_err_min)) + do itt=1,maxitt + Vol_err = 0.5*(L(K)*L(K))*(slope + crv_3*(3.0-4.0*L(K))) - vol_below(K) + if (abs(Vol_err) <= Vol_quit) exit + ! Take a Newton's method iteration. This equation has proven + ! robust enough not to need bracketing. + L(K) = L(K) - Vol_err / (L(K)* (slope + crv - 2.0*crv*L(K))) + ! This would be a Newton's method iteration for L^2: + ! L(K) = sqrt(L(K)*L(K) - Vol_err / (0.5*(slope+crv) - crv*L(K))) + enddo + endif ! end of iterative solver + endif ! end of 1-boundary alternatives. + endif ! end of 0, 1- and 2- boundary cases. + enddo ! k loop to determine L(K) in the convex case + +end subroutine find_L_open_convex + !> This subroutine finds a thickness-weighted value of v at the u-points. function set_v_at_u(v, h, G, GV, i, j, k, mask2dCv, OBC) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure From da02ffc601ff62a8ab800eb1f759a5105780d071 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 29 Dec 2023 23:23:56 -0500 Subject: [PATCH 06/15] +Add find_L_open_concave_iterative Added the new routine find_L_open_concave_iterative to use iterative Newton's method approaches with appropriate limits to solve the cubic equation for the fractional open face lengths at interfaces that are used by the CHANNEL_DRAG code. These solutions are analogous to those given by the previous expressions that are now in find_L_open_concave_trigonometric, and the two differ at close to roundoff, but the new method is completely independent of the transcendental function library, thereby addressing dev/gfdl MOM6 issue #483. This new routine is called when the new runtime parameter TRIG_CHANNEL_DRAG_WIDTHS is set to false, but by default the previous answers are recovered. By default all answers are bitwise identical, but there is a new runtime parameter in some MOM_parameter_doc files. --- .../vertical/MOM_set_viscosity.F90 | 389 +++++++++++++++++- 1 file changed, 370 insertions(+), 19 deletions(-) diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index bf0db51a44..b25d660957 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -20,7 +20,7 @@ module MOM_set_visc use MOM_grid, only : ocean_grid_type use MOM_hor_index, only : hor_index_type use MOM_interface_heights, only : thickness_to_dz -! use MOM_intrinsic_functions, only : cuberoot +use MOM_intrinsic_functions, only : cuberoot use MOM_io, only : slasher, MOM_read_data, vardesc, var_desc use MOM_kappa_shear, only : kappa_shear_is_used, kappa_shear_at_vertex use MOM_open_boundary, only : ocean_OBC_type, OBC_segment_type, OBC_NONE, OBC_DIRECTION_E @@ -101,6 +101,8 @@ module MOM_set_visc real :: omega_frac !< When setting the decay scale for turbulence, use this !! fraction of the absolute rotation rate blended with the local !! value of f, as sqrt((1-of)*f^2 + of*4*omega^2) [nondim] + logical :: concave_trigonometric_L !< If true, use trigonometric expressions to determine the + !! fractional open interface lengths for concave topography. integer :: answer_date !< The vintage of the order of arithmetic and expressions in the set !! viscosity calculations. Values below 20190101 recover the answers !! from the end of 2018, while higher values use updated and more robust @@ -208,11 +210,11 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) ! from lateral lengths to layer thicknesses [H L-1 ~> nondim or kg m-3]. real :: cdrag_sqrt_H_RL ! Square root of the drag coefficient, times a unit conversion factor from ! density times lateral lengths to layer thicknesses [H L-1 R-1 ~> m3 kg-1 or nondim] - real :: cdrag_L_to_H ! The drag coeffient times conversion factors from lateral + real :: cdrag_L_to_H ! The drag coefficient times conversion factors from lateral ! distance to thickness units [H L-1 ~> nondim or kg m-3] - real :: cdrag_RL_to_H ! The drag coeffient times conversion factors from density times lateral + real :: cdrag_RL_to_H ! The drag coefficient times conversion factors from density times lateral ! distance to thickness units [H L-1 R-1 ~> m3 kg-1 or nondim] - real :: cdrag_conv ! The drag coeffient times a combination of static conversion factors and in + real :: cdrag_conv ! The drag coefficient times a combination of static conversion factors and in ! situ density or Boussinesq reference density [H L-1 ~> nondim or kg m-3] real :: oldfn ! The integrated energy required to ! entrain up to the bottom of the layer, @@ -268,6 +270,11 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) ! horizontal area of a velocity cell [Z ~> m]. real :: L(SZK_(GV)+1) ! The fraction of the full cell width that is open at ! the depth of each interface [nondim]. + real :: L_trig(SZK_(GV)+1) ! The fraction of the full cell width that is open at + ! the depth of each interface from trigonometric expressions [nondim]. + real :: dL_trig_itt(SZK_(GV)+1) ! The difference between estimates of the fraction of the full cell + ! width that is open at the depth of each interface [nondim]. + real :: max_dL_trig_itt ! The largest difference between L and L_trig, for debugging [nondim] real :: h_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected [H ~> m or kg m-2]. real :: dz_neglect ! A vertical distance that is so small it is usually lost @@ -866,18 +873,32 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) if (crv == 0.0) then call find_L_open_uniform_slope(vol_below, Dp, Dm, L, GV) elseif (crv > 0.0) then - call find_L_open_concave_analytic(vol_below, D_vel, Dp, Dm, C2pi_3, L, GV) + if (CS%concave_trigonometric_L) then + call find_L_open_concave_trigonometric(vol_below, D_vel, Dp, Dm, C2pi_3, L, GV) + else + call find_L_open_concave_iterative(vol_below, D_vel, Dp, Dm, L, GV) + if (CS%debug) then + call find_L_open_concave_trigonometric(vol_below, D_vel, Dp, Dm, C2pi_3, L_trig, GV) + max_dL_trig_itt = 0.0 + do K=1,nz+1 + dL_trig_itt(K) = L_trig(K) - L(K) + if (abs(dL_trig_itt(K)) > abs(max_dL_trig_itt)) max_dL_trig_itt = dL_trig_itt(K) + enddo + if (abs(max_dL_trig_itt) > 1.0e-12) & + K = nz+1 ! This is here to use with a debugger only. + endif + endif else ! crv < 0.0 call find_L_open_convex(vol_below, D_vel, Dp, Dm, L, GV, US, CS) endif ! end of crv<0 cases. - ! Determine the Rayliegh drag contributions. + ! Determine the Rayleigh drag contributions. ! The drag within the bottommost Vol_bbl_chan is applied as a part of an enhanced bottom ! viscosity, while above this the drag is applied directly to the layers in question as a ! Rayleigh drag term. - ! Restrict the volume over which the channel drag is applied from the previously determined valu.e + ! Restrict the volume over which the channel drag is applied from the previously determined value. if (CS%Chan_drag_max_vol >= 0.0) Vol_bbl_chan = min(Vol_bbl_chan, CS%Chan_drag_max_vol) BBL_visc_frac = 0.0 @@ -1082,9 +1103,9 @@ subroutine find_L_open_uniform_slope(vol_below, Dp, Dm, L, GV) end subroutine find_L_open_uniform_slope -!> Determine the normalized open length of each interface for concave bathymetry (from the ocean perspective) using -!! analytic expressions. In this case there can be two separate open regions. -subroutine find_L_open_concave_analytic(vol_below, D_vel, Dp, Dm, C2pi_3, L, GV) +!> Determine the normalized open length of each interface for concave bathymetry (from the ocean perspective) +!! using trigonometric expressions. In this case there can be two separate open regions. +subroutine find_L_open_concave_trigonometric(vol_below, D_vel, Dp, Dm, C2pi_3, L, GV) type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. real, dimension(SZK_(GV)+1), intent(in) :: vol_below !< The volume below each interface, normalized by !! the full horizontal area of a velocity cell [Z ~> m] @@ -1111,7 +1132,7 @@ subroutine find_L_open_concave_analytic(vol_below, D_vel, Dp, Dm, C2pi_3, L, GV) real :: C24_crv ! 24/crv [Z-1 ~> m-1]. real :: apb_4a, ax2_3apb ! Various nondimensional ratios of crv and slope [nondim]. real :: a2x48_apb3, Iapb ! Combinations of crv (a) and slope (b) [Z-1 ~> m-1] - real :: L0 ! A linear estimate of L approprate for tiny volumes [nondim]. + real :: L0 ! A linear estimate of L appropriate for tiny volumes [nondim]. real :: slope_crv ! The slope divided by the curvature [nondim] real :: tmp_val_m1_to_p1 ! A temporary variable [nondim] real, parameter :: C1_3 = 1.0/3.0, C1_12 = 1.0/12.0 ! Rational constants [nondim] @@ -1122,7 +1143,7 @@ subroutine find_L_open_concave_analytic(vol_below, D_vel, Dp, Dm, C2pi_3, L, GV) ! Each cell extends from x=-1/2 to 1/2, and has a topography ! given by D(x) = crv*x^2 + slope*x + D_vel - crv/12. crv_3 = (Dp + Dm - 2.0*D_vel) ; crv = 3.0*crv_3 - if (crv < 0.0) call MOM_error(FATAL, "find_L_open_concave should only be called with a positive curvature.") + if (crv <= 0.0) call MOM_error(FATAL, "find_L_open_concave should only be called with a positive curvature.") slope = Dp - Dm ! Calculate the volume above which the entire cell is open and the volume at which the @@ -1168,12 +1189,337 @@ subroutine find_L_open_concave_analytic(vol_below, D_vel, Dp, Dm, C2pi_3, L, GV) endif enddo ! k loop to determine L(K) in the concave case -end subroutine find_L_open_concave_analytic +end subroutine find_L_open_concave_trigonometric + + + +!> Determine the normalized open length of each interface for concave bathymetry (from the ocean perspective) using +!! iterative methods to solve the relevant cubic equations. In this case there can be two separate open regions. +subroutine find_L_open_concave_iterative(vol_below, D_vel, Dp, Dm, L, GV) + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + real, dimension(SZK_(GV)+1), intent(in) :: vol_below !< The volume below each interface, normalized by + !! the full horizontal area of a velocity cell [Z ~> m] + real, intent(in) :: D_vel !< The average bottom depth at a velocity point [Z ~> m] + real, intent(in) :: Dp !< The larger of the two depths at the edge + !! of a velocity cell [Z ~> m] + real, intent(in) :: Dm !< The smaller of the two depths at the edge + !! of a velocity cell [Z ~> m] + real, dimension(SZK_(GV)+1), intent(out) :: L !< The fraction of the full cell width that is open at + !! the depth of each interface [nondim] + + ! Local variables + real :: crv ! crv is the curvature of the bottom depth across a + ! cell, times the cell width squared [Z ~> m]. + real :: crv_3 ! crv/3 [Z ~> m]. + real :: slope ! The absolute value of the bottom depth slope across + ! a cell times the cell width [Z ~> m]. + + ! The following "volumes" have units of vertical heights because they are normalized + ! by the full horizontal area of a velocity cell. + real :: Vol_open ! The cell volume above which the face is fully is open [Z ~> m]. + real :: Vol_2_reg ! The cell volume above which there are two separate + ! open areas that must be integrated [Z ~> m]. + real :: L_2_reg ! The value of L when vol_below is Vol_2_reg [nondim] + real :: vol_inflect_1 ! The volume at which there is an inflection point in the expression + ! relating L to vol_err when there is a single open region [Z ~> m] + real :: vol_inflect_2 ! The volume at which there is an inflection point in the expression + ! relating L to vol_err when there are two open regions [Z ~> m] + + real :: L_inflect_1 ! The value of L that sits at an inflection point in the expression + ! relating L to vol_err when there is a single open region [nondim] + real :: L_inflect_2 ! The value of L that sits at an inflection point in the expression + ! relating L to vol_err when there is are two open regions [nondim] + real :: L_max, L_min ! Maximum and minimum bounds on the solution for L for an interface [nondim] + real :: vol_err ! The difference between the volume below an interface for a given value + ! of L and the target value [Z ~> m] + real :: dVol_dL ! The partial derivative of the volume below with L [Z ~> m] + real :: vol_err_max ! The value of vol_err when L is L_max [Z ~> m] + + ! The following combinations of slope and crv are reused across layers, and hence are pre-calculated + ! for efficiency. All are non-negative. + real :: Icrvpslope ! The inverse of the sum of crv and slope [Z-1 ~> m-1] + real :: slope_crv ! The slope divided by the curvature [nondim] + ! These are only used if the slope exceeds or matches the curvature. + real :: smc ! The slope minus the curvature [Z ~> m] + real :: C3c_m_s ! 3 times the curvature minus the slope [Z ~> m] + real :: I_3c_m_s ! The inverse of 3 times the curvature minus the slope [Z-1 ~> m-1] + ! These are only used if the curvature exceeds the slope. + real :: C4_crv ! The inverse of a quarter of the curvature [Z-1 ~> m-1] + real :: sxcms_c ! The slope times the difference between the curvature and slope + ! divided by the curvature [Z ~> m] + real :: slope2_4crv ! A quarter of the slope squared divided by the curvature [Z ~> m] + real :: I_3s_m_c ! The inverse of 3 times the slope minus the curvature [Z-1 ~> m-1] + real :: C3s_m_c ! 3 times the slope minus the curvature [Z ~> m] + + real, parameter :: C1_3 = 1.0 / 3.0, C1_12 = 1.0 / 12.0 ! Rational constants [nondim] + integer :: K, nz, itt + integer, parameter :: max_itt = 10 + + nz = GV%ke + + ! Each cell extends from x=-1/2 to 1/2, and has a topography + ! given by D(x) = crv*x^2 + slope*x + D_vel - crv/12. + + crv_3 = (Dp + Dm - 2.0*D_vel) ; crv = 3.0*crv_3 + if (crv <= 0.0) call MOM_error(FATAL, "find_L_open_concave should only be called with a positive curvature.") + slope = Dp - Dm + + ! Calculate the volume above which the entire cell is open and the volume at which the + ! equation that is solved for L changes because there are two separate open regions. + if (slope >= crv) then + Vol_open = D_vel - Dm ; Vol_2_reg = Vol_open + L_2_reg = 1.0 + if (crv + slope >= 4.0*crv) then + L_inflect_1 = 1.0 ; Vol_inflect_1 = Vol_open + else + slope_crv = slope / crv + L_inflect_1 = 0.25 + 0.25*slope_crv + vol_inflect_1 = 0.25*C1_12 * ((slope_crv + 1.0)**2 * (slope + crv)) + endif + ! Precalculate some combinations of crv & slope for later use. + smc = slope - crv + C3c_m_s = 3.0*crv - slope + if (C3c_m_s > 2.0*smc) I_3c_m_s = 1.0 / C3c_m_s + else + slope_crv = slope / crv + Vol_open = 0.25*slope*slope_crv + C1_12*crv + Vol_2_reg = 0.5*slope_crv**2 * (crv - C1_3*slope) + L_2_reg = slope_crv + + ! The inflection point is useful to know because below the inflection point + ! Newton's method converges monotonically from above and conversely above it. + ! These are the inflection point values of L and vol_below with a single open segment. + vol_inflect_1 = 0.25*C1_12 * ((slope_crv + 1.0)**2 * (slope + crv)) + L_inflect_1 = 0.25 + 0.25*slope_crv + ! These are the inflection point values of L and vol_below when there are two open segments. + ! Vol_inflect_2 = Vol_open - 0.125 * crv_3, which is equivalent to: + vol_inflect_2 = 0.25*slope*slope_crv + 0.125*crv_3 + L_inflect_2 = 0.5 + ! Precalculate some combinations of crv & slope for later use. + C4_crv = 4.0 / crv + slope2_4crv = 0.25 * slope * slope_crv + sxcms_c = slope_crv*(crv - slope) + C3s_m_c = 3.0*slope - crv + if (C3s_m_c > 2.0*sxcms_c) I_3s_m_c = 1.0 / C3s_m_c + endif + ! Define some combinations of crv & slope for later use. + Icrvpslope = 1.0 / (crv+slope) + + L(nz+1) = 0.0 + ! Determine the normalized open length (L) at each interface. + do K=nz,1,-1 + if (vol_below(K) >= Vol_open) then ! The whole cell is open. + L(K) = 1.0 + elseif (vol_below(K) < Vol_2_reg) then + ! In this case, there is a single contiguous open region from x=1/2-L to 1/2. + ! Changing the horizontal variable in the expression from D(x) to D(L) gives: + ! x(L) = 1/2 - L + ! D(L) = crv*(0.5 - L)^2 + slope*(0.5 - L) + D_vel - crv/12 + ! D(L) = crv*L^2 - crv*L + crv/4 + slope*(1/2 - L) + D_vel - crv/12 + ! D(L) = crv*L^2 - (slope+crv)*L + slope/2 + D_vel + crv/6 + ! D(0) = slope/2 + D_vel + crv/6 = (Dp - Dm)/2 + D_vel + (Dp + Dm - 2*D_vel)/2 = Dp + ! D(1) = crv - slope - crv + slope/2 + Dvel + crv/6 = D_vel - slope/2 + crv/6 = Dm + ! + ! vol_below = integral(y = 0 to L) D(y) dy - L * D(L) + ! = crv/3*L^3 - (slope+crv)/2*L^2 + (slope/2 + D_vel + crv/6)*L - + ! (crv*L^2 - (slope+crv)*L + slope/2 + D_vel + crv/6) * L + ! = -2/3 * crv * L^3 + 1/2 * (slope+crv) * L^2 + ! vol_below(K) = 0.5*L(K)**2*(slope + crv_3*(3-4*L(K))) + ! L(K) is between L(K+1) and slope_crv. + L_max = min(L_2_reg, 1.0) + if (vol_below(K) <= vol_inflect_1) L_max = min(L_max, L_inflect_1) + + L_min = L(K+1) + if (vol_below(K) >= vol_inflect_1) L_min = max(L_min, L_inflect_1) + + ! Ignoring the cubic term gives an under-estimate but is very accurate for near bottom + ! layers, so use this as a potential floor. + if (2.0*vol_below(K)*Icrvpslope > L_min**2) L_min = sqrt(2.0*vol_below(K)*Icrvpslope) + + ! Start with L_min in most cases. + L(k) = L_min + + if (vol_below(K) <= vol_inflect_1) then + ! Starting with L_min below L_inflect_1, only the first overshooting iteration of Newton's + ! method needs bounding. + L(k) = L_min + vol_err = 0.5*L(K)**2 * (slope + crv*(1.0 - 4.0*C1_3*L(K))) - vol_below(K) + ! If vol_err is 0 or positive (perhaps due to roundoff in L(K+1)), L_min is already the best solution. + if (vol_err < 0.0) then + dVol_dL = L(K) * (slope + crv*(1.0 - 2.0*L(k))) + if (L(K)*dVol_dL > vol_err + L_max*dVol_dL) then + L(K) = L_max + else + L(K) = L(K) - (vol_err / dVol_dL) + endif + + ! Subsequent iterations of Newton's method do not need bounds. + do itt=1,max_itt + vol_err = 0.5*L(K)**2 * (slope + crv*(1.0 - 4.0*C1_3*L(K))) - vol_below(K) + dVol_dL = L(K) * (slope + crv*(1.0 - 2.0*L(k))) + if (abs(vol_err) < max(1.0e-15*L(K), 1.0e-25)*dVol_dL) exit + L(K) = L(K) - (vol_err / dVol_dL) + enddo + endif + else ! (vol_below(K) > vol_inflect_1) + ! Iteration from below converges monotonically, but we need to deal with the case where we are + ! close to the peak of the topography and Newton's method mimics the convergence of bisection. + + ! Evaluate the error when L(K) = L_min as a possible first guess. + L(k) = L_min + vol_err = 0.5*L(K)**2 * (slope + crv*(1.0 - 4.0*C1_3*L(K))) - vol_below(K) + ! If vol_err is 0 or positive (perhaps due to roundoff in L(K+1)), L_min is already the best solution. + if (vol_err < 0.0) then + + ! These two upper estimates deal with the possibility that this point may be near + ! the upper extrema, where the error term might be approximately parabolic and + ! Newton's method would converge slowly like simple bisection. + if (slope < crv) then + ! if ((L_2_reg - L_min)*(3.0*slope - crv) > 2.0*slope_crv*(crv-slope)) then + if ((L_2_reg - L_min)*C3s_m_c > 2.0*sxcms_c) then + ! There is a decent upper estimate of L from the approximate quadratic equation found + ! by examining the error expressions at L ~= L_2_reg and ignoring the cubic term. + L_max = (slope_crv*(2.0*slope) - sqrt(sxcms_c**2 + & + 2.0*C3s_m_c*(Vol_2_reg - vol_below(K))) ) * I_3s_m_c + ! The line above is equivalent to: + ! L_max = (slope_crv*(2.0*slope) - sqrt(slope_crv**2*(crv-slope)**2 + & + ! 2.0*(3.0*slope - crv)*(Vol_2_reg - vol_below(K))) ) / & + ! (3.0*slope - crv) + else + L_max = slope_crv + endif + else ! (slope >= crv) + if ((1.0 - L_min)*C3c_m_s > 2.0*smc) then + ! There is a decent upper estimate of L from the approximate quadratic equation found + ! by examining the error expressions at L ~= 1 and ignoring the cubic term. + L_max = ( 2.0*crv - sqrt(smc**2 + 2.0*C3c_m_s * (Vol_open - vol_below(K))) ) * I_3c_m_s + ! The line above is equivalent to: + ! L_max = ( 2.0*crv - sqrt((slope - crv)**2 + 2.0*(3.0*crv - slope) * (Vol_open - vol_below(K))) ) / & + ! (3.0*crv - slope) + else + L_max = 1.0 + endif + endif + Vol_err_max = 0.5*L_max**2 * (slope + crv*(1.0 - 4.0*C1_3*L_max)) - vol_below(K) + if (Vol_err_max < 0.0) call MOM_error(FATAL, & + "Vol_err_max should never be negative in find_L_open_concave_iterative.") + if ((Vol_err_max < abs(Vol_err)) .and. (L_max < 1.0)) then + ! Start with 1 bounded Newton's method step from L_max + dVol_dL = L_max * (slope + crv*(1.0 - 2.0*L_max)) + L(K) = max(L_min, L_max - (vol_err_max / dVol_dL) ) + ! else ! Could use the fact that Vol_err is known to take an iteration? + endif + ! Subsequent iterations of Newton's method do not need bounds. + do itt=1,max_itt + vol_err = 0.5*L(K)**2 * (slope + crv*(1.0 - 4.0*C1_3*L(K))) - vol_below(K) + dVol_dL = L(K) * (slope + crv*(1.0 - 2.0*L(k))) + if (abs(vol_err) < max(1.0e-15*L(K), 1.0e-25)*dVol_dL) exit + L(K) = L(K) - (vol_err / dVol_dL) + enddo + endif -!> Determine the normalized open length of each interface for convex bathymetry (from the ocean perspective) using -!! Newton's method iterations. In this case there is a single open region with the minimum depth -!! at one edge of the cell. + endif + + ! To check the answers. + ! Vol_err = 0.5*(L(K)*L(K))*(slope + crv_3*(3.0-4.0*L(K))) - vol_below(K) + else ! There are two separate open regions. + ! vol_below(K) = slope^2/(4*crv) + crv/12 - (crv/12)*(1-L)^2*(1+2L) + ! At the deepest volume, L = slope/crv, at the top L = 1. + + ! To check the answers. + ! Vol_err = Vol_open - 0.25*crv_3*(1.0+2.0*L(K)) * (1.0-L(K))**2 - vol_below(K) + ! or equivalently: + ! Vol_err = Vol_open - 0.25*crv_3*(3.0-2.0*(1.0-L(K))) * (1.0-L(K))**2 - vol_below(K) + ! ! Note that: Vol_open = 0.25*slope*slope_crv + C1_12*crv + ! Vol_err = 0.25*slope*slope_crv + 0.25*crv_3*( 1.0 - (1.0 + 2.0*L(K)) * (1.0-L(K))**2 ) - vol_below(K) + ! Vol_err = 0.25*crv_3*L(K)**2*( 3.0 - 2.0*L(K) ) + 0.25*slope*slope_crv - vol_below(K) + + ! Derivation of the L_max limit below: + ! Vol_open - vol_below(K) = 0.25*crv_3*(3.0-2.0*(1.0-L(K))) * (1.0-L(K))**2 + ! (3.0-2.0*(1.0-L(K))) * (1.0-L(K))**2 = (Vol_open - vol_below(K)) / (0.25*crv_3) + ! When 1-L(K) << 1: + ! 3.0 * (1.0-L_max)**2 = (Vol_open - vol_below(K)) / (0.25*crv_3) + ! (1.0-L_max)**2 = (Vol_open - vol_below(K)) / (0.25*crv) + + ! Derivation of the L_min limit below: + ! Vol_err = 0.25*crv_3*L(K)**2*( 3.0 - 2.0*L(K) ) + 0.25*slope*slope_crv - vol_below(K) + ! crv*L(K)**2*( 1.0 - 2.0*C1_3*L(K) ) = 4.0*vol_below(K) - slope*slope_crv + ! When L(K) << 1: + ! crv*L_min**2 = 4.0*vol_below(K) - slope*slope_crv + ! L_min = sqrt((4.0*vol_below(K) - slope*slope_crv)/crv) + ! Noting that L(K) >= slope_crv, when L(K)-slope_crv << 1: + ! (crv + 2.0*C1_3*slope)*L_min**2 = 4.0*vol_below(K) - slope*slope_crv + ! L_min = sqrt((4.0*vol_below(K) - slope*slope_crv)/(crv + 2.0*C1_3*slope)) + + if (vol_below(K) <= Vol_inflect_2) then + ! Newton's Method would converge monotonically from above, but overshoot from below. + L_min = max(L(K+1), L_2_reg) ! L_2_reg = slope_crv + ! This under-estimate of L(K) is accurate for L ~= slope_crv: + if ((4.0*vol_below(K) - slope*slope_crv) > (crv + 2.0*C1_3*slope)*L_min**2) & + L_min = max(L_min, sqrt((4.0*vol_below(K) - slope*slope_crv) / (crv + 2.0*C1_3*slope))) + L_max = 0.5 ! = L_inflect_2 + + ! Starting with L_min below L_inflect_2, only the first overshooting iteration of Newton's + ! method needs bounding. + L(k) = L_min + Vol_err = crv_3*L(K)**2*( 0.75 - 0.5*L(K) ) + (slope2_4crv - vol_below(K)) + + ! If vol_err is 0 or positive (perhaps due to roundoff in L(K+1)), L_min is already the best solution. + if (vol_err < 0.0) then + dVol_dL = 0.5*crv * (L(K) * (1.0 - L(K))) + if (L(K)*dVol_dL >= vol_err + L_max*dVol_dL) then + L(K) = L_max + else + L(K) = L(K) - (vol_err / dVol_dL) + endif + ! Subsequent iterations of Newton's method do not need bounds. + do itt=1,max_itt + Vol_err = crv_3 * (L(K)**2 * (0.75 - 0.5*L(K))) + (slope2_4crv - vol_below(K)) + dVol_dL = 0.5*crv * (L(K)*(1.0 - L(K))) + if (abs(vol_err) < max(1.0e-15*L(K), 1.0e-25)*dVol_dL) exit + L(K) = L(K) - (vol_err / dVol_dL) + enddo + endif + else ! (vol_below(K) > Vol_inflect_2) + ! Newton's Method would converge monotonically from below, but overshoots from above, and + ! we may need to deal with the case where we are close to the peak of the topography. + L_min = max(L(K+1), 0.5) + L(k) = L_min + + Vol_err = crv_3 * (L(K)**2 * ( 0.75 - 0.5*L(K))) + (slope2_4crv - vol_below(K)) + ! If vol_err is 0 or positive (perhaps due to roundoff in L(K+1)), L(k) is already the best solution. + if (Vol_err < 0.0) then + ! This over-estimate of L(K) is accurate for L ~= 1: + L_max = 1.0 - sqrt( (Vol_open - vol_below(K)) * C4_crv ) + Vol_err_max = crv_3 * (L_max**2 * ( 0.75 - 0.5*L_max)) + (slope2_4crv - vol_below(K)) + if (Vol_err_max < 0.0) call MOM_error(FATAL, & + "Vol_err_max should never be negative in find_L_open_concave_iterative.") + if ((Vol_err_max < abs(Vol_err)) .and. (L_max < 1.0)) then + ! Start with 1 bounded Newton's method step from L_max + dVol_dL = 0.5*crv * (L_max * (1.0 - L_max)) + L(K) = max(L_min, L_max - (vol_err_max / dVol_dL) ) + ! else ! Could use the fact that Vol_err is known to take an iteration? + endif + + ! Subsequent iterations of Newton's method do not need bounds. + do itt=1,max_itt + Vol_err = crv_3 * (L(K)**2 * ( 0.75 - 0.5*L(K))) + (slope2_4crv - vol_below(K)) + dVol_dL = 0.5*crv * (L(K) * (1.0 - L(K))) + if (abs(vol_err) < max(1.0e-15*L(K), 1.0e-25)*dVol_dL) exit + L(K) = L(K) - (vol_err / dVol_dL) + enddo + endif + endif + + endif + enddo ! k loop to determine L(K) in the concave case + +end subroutine find_L_open_concave_iterative + +!> Determine the normalized open length of each interface for convex bathymetry (from the ocean +!! perspective) using Newton's method iterations. In this case there is a single open region +!! with the minimum depth at one edge of the cell. subroutine find_L_open_convex(vol_below, D_vel, Dp, Dm, L, GV, US, CS) type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. real, dimension(SZK_(GV)+1), intent(in) :: vol_below !< The volume below each interface, normalized by @@ -1226,7 +1572,7 @@ subroutine find_L_open_convex(vol_below, D_vel, Dp, Dm, L, GV, US, CS) ! Each cell extends from x=-1/2 to 1/2, and has a topography ! given by D(x) = crv*x^2 + slope*x + D_vel - crv/12. crv_3 = (Dp + Dm - 2.0*D_vel) ; crv = 3.0*crv_3 - if (crv > 0.0) call MOM_error(FATAL, "find_L_open_convex should only be called with a negative curvature.") + if (crv >= 0.0) call MOM_error(FATAL, "find_L_open_convex should only be called with a negative curvature.") slope = Dp - Dm ! Calculate the volume above which the entire cell is open and the volume at which the @@ -1252,7 +1598,7 @@ subroutine find_L_open_convex(vol_below, D_vel, Dp, Dm, L, GV, US, CS) elseif (vol_below(K) <= Vol_direct) then ! Both edges of the cell are bounded by walls. ! if (CS%answer_date < 20240101)) then - L(K) = (-0.25*C24_crv*vol_below(K))**C1_3 + L(K) = (-0.25*C24_crv*vol_below(K))**C1_3 ! else ! L(K) = cuberoot(-0.25*C24_crv*vol_below(K)) ! endif @@ -2592,8 +2938,13 @@ subroutine set_visc_init(Time, G, GV, US, param_file, diag, visc, CS, restart_CS "default is to use the same value as SMAG_LAP_CONST if "//& "it is defined, or 0.15 if it is not. The value used is "//& "also 0.15 if the specified value is negative.", & - units="nondim", default=cSmag_chan_dflt) + units="nondim", default=cSmag_chan_dflt, do_not_log=.not.CS%Channel_drag) if (CS%c_Smag < 0.0) CS%c_Smag = 0.15 + + call get_param(param_file, mdl, "TRIG_CHANNEL_DRAG_WIDTHS", CS%concave_trigonometric_L, & + "If true, use trigonometric expressions to determine the fractional open "//& + "interface lengths for concave topography.", & + default=.true., do_not_log=.not.CS%Channel_drag) endif Chan_max_thick_dflt = -1.0*US%m_to_Z From e1ea75847e97b4dc2a2dd375d4826afe493b70c0 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 15 Jan 2024 16:02:41 -0500 Subject: [PATCH 07/15] Add the debugging routine test_L_open_concave Added the new debugging or testing subroutine test_L_open_concave along with extra calls when DEBUG = True that can be used to demonstrate that the iterative solver in find_L_open_concave_iterative is substantially more accurate but mathematically equivalent to the solver in find_L_open_concave_trigonometric. This extra code is only called in debugging mode, and it probably should be deleted in a separate commit after find_L_open_concave_iterative has been accepted onto the dev/gfdl branch of MOM6. All answers are bitwise identical and no output or input is changed. --- .../vertical/MOM_set_viscosity.F90 | 112 +++++++++++++++++- 1 file changed, 109 insertions(+), 3 deletions(-) diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index b25d660957..1778b8d870 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -270,11 +270,19 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) ! horizontal area of a velocity cell [Z ~> m]. real :: L(SZK_(GV)+1) ! The fraction of the full cell width that is open at ! the depth of each interface [nondim]. + ! The next 9 variables are only used for debugging. real :: L_trig(SZK_(GV)+1) ! The fraction of the full cell width that is open at ! the depth of each interface from trigonometric expressions [nondim]. + real :: vol_err_trig(SZK_(GV)+1) ! The error in the volume below based on L_trig [Z ~> m] + real :: vol_err_iter(SZK_(GV)+1) ! The error in the volume below based on L_iter [Z ~> m] + real :: norm_err_trig(SZK_(GV)+1) ! vol_err_trig normalized by vol_below [nondim] + real :: norm_err_iter(SZK_(GV)+1) ! vol_err_iter normalized by vol_below [nondim] real :: dL_trig_itt(SZK_(GV)+1) ! The difference between estimates of the fraction of the full cell ! width that is open at the depth of each interface [nondim]. real :: max_dL_trig_itt ! The largest difference between L and L_trig, for debugging [nondim] + real :: max_norm_err_trig ! The largest magnitude value of norm_err_trig in a column [nondim] + real :: max_norm_err_iter ! The largest magnitude value of norm_err_iter in a column [nondim] + real :: h_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected [H ~> m or kg m-2]. real :: dz_neglect ! A vertical distance that is so small it is usually lost @@ -878,14 +886,29 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) else call find_L_open_concave_iterative(vol_below, D_vel, Dp, Dm, L, GV) if (CS%debug) then + ! The tests in this block reveal that the iterative and trigonometric solutions are + ! mathematically equiavalent, but in some cases the iterative solution is consistent + ! at roundoff, but that the trigonmetric solutions have errors that can be several + ! orders of magnitude larger in some cases. call find_L_open_concave_trigonometric(vol_below, D_vel, Dp, Dm, C2pi_3, L_trig, GV) - max_dL_trig_itt = 0.0 + call test_L_open_concave(vol_below, D_vel, Dp, Dm, L_trig, vol_err_trig, GV) + call test_L_open_concave(vol_below, D_vel, Dp, Dm, L, vol_err_iter, GV) + max_dL_trig_itt = 0.0 ; max_norm_err_trig = 0.0 ; max_norm_err_iter = 0.0 + norm_err_trig(:) = 0.0 ; norm_err_iter(:) = 0.0 do K=1,nz+1 dL_trig_itt(K) = L_trig(K) - L(K) if (abs(dL_trig_itt(K)) > abs(max_dL_trig_itt)) max_dL_trig_itt = dL_trig_itt(K) + norm_err_trig(K) = vol_err_trig(K) / (vol_below(K) + dz_neglect) + norm_err_iter(K) = vol_err_iter(K) / (vol_below(K) + dz_neglect) + if (abs(norm_err_trig(K)) > abs(max_norm_err_trig)) max_norm_err_trig = norm_err_trig(K) + if (abs(norm_err_iter(K)) > abs(max_norm_err_iter)) max_norm_err_iter = norm_err_iter(K) enddo - if (abs(max_dL_trig_itt) > 1.0e-12) & - K = nz+1 ! This is here to use with a debugger only. + if (abs(max_dL_trig_itt) > 1.0e-13) & + K = nz+1 ! This is here only to use as a break point for a debugger. + if (abs(max_norm_err_trig) > 1.0e-13) & + K = nz+1 ! This is here only to use as a break point for a debugger. + if (abs(max_norm_err_iter) > 1.0e-13) & + K = nz+1 ! This is here only to use as a break point for a debugger. endif endif else ! crv < 0.0 @@ -1517,6 +1540,89 @@ subroutine find_L_open_concave_iterative(vol_below, D_vel, Dp, Dm, L, GV) end subroutine find_L_open_concave_iterative + + +!> Test the validity the normalized open lengths of each interface for concave bathymetry (from the ocean perspective) +!! by evaluating and returing the relevant cubic equations. +subroutine test_L_open_concave(vol_below, D_vel, Dp, Dm, L, vol_err, GV) + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + real, dimension(SZK_(GV)+1), intent(in) :: vol_below !< The volume below each interface, normalized by + !! the full horizontal area of a velocity cell [Z ~> m] + real, intent(in) :: D_vel !< The average bottom depth at a velocity point [Z ~> m] + real, intent(in) :: Dp !< The larger of the two depths at the edge + !! of a velocity cell [Z ~> m] + real, intent(in) :: Dm !< The smaller of the two depths at the edge + !! of a velocity cell [Z ~> m] + real, dimension(SZK_(GV)+1), intent(in) :: L !< The fraction of the full cell width that is open at + !! the depth of each interface [nondim] + real, dimension(SZK_(GV)+1), intent(out) :: vol_err !< The difference between vol_below and the + !! value obtained from using L in the cubic equation [Z ~> m] + + ! Local variables + real :: crv ! crv is the curvature of the bottom depth across a + ! cell, times the cell width squared [Z ~> m]. + real :: crv_3 ! crv/3 [Z ~> m]. + real :: slope ! The absolute value of the bottom depth slope across + ! a cell times the cell width [Z ~> m]. + + ! The following "volumes" have units of vertical heights because they are normalized + ! by the full horizontal area of a velocity cell. + real :: Vol_open ! The cell volume above which the face is fully is open [Z ~> m]. + real :: Vol_2_reg ! The cell volume above which there are two separate + ! open areas that must be integrated [Z ~> m]. + real :: L_2_reg ! The value of L when vol_below is Vol_2_reg [nondim] + + ! The following combinations of slope and crv are reused across layers, and hence are pre-calculated + ! for efficiency. All are non-negative. + real :: slope_crv ! The slope divided by the curvature [nondim] + ! These are only used if the curvature exceeds the slope. + real :: slope2_4crv ! A quarter of the slope squared divided by the curvature [Z ~> m] + + real, parameter :: C1_3 = 1.0 / 3.0, C1_12 = 1.0 / 12.0 ! Rational constants [nondim] + integer :: K, nz + + nz = GV%ke + + ! Each cell extends from x=-1/2 to 1/2, and has a topography + ! given by D(x) = crv*x^2 + slope*x + D_vel - crv/12. + + crv_3 = (Dp + Dm - 2.0*D_vel) ; crv = 3.0*crv_3 + if (crv <= 0.0) call MOM_error(FATAL, "test_L_open_concave should only be called with a positive curvature.") + slope = Dp - Dm + + ! Calculate the volume above which the entire cell is open and the volume at which the + ! equation that is solved for L changes because there are two separate open regions. + if (slope >= crv) then + Vol_open = D_vel - Dm ; Vol_2_reg = Vol_open + L_2_reg = 1.0 + if (crv + slope >= 4.0*crv) then + slope_crv = 1.0 + else + slope_crv = slope / crv + endif + else + slope_crv = slope / crv + Vol_open = 0.25*slope*slope_crv + C1_12*crv + Vol_2_reg = 0.5*slope_crv**2 * (crv - C1_3*slope) + L_2_reg = slope_crv + endif + slope2_4crv = 0.25 * slope * slope_crv + + ! Determine the volume error based on the normalized open length (L) at each interface. + Vol_err(nz+1) = 0.0 + do K=nz,1,-1 + if (L(K) >= 1.0) then + Vol_err(K) = max(Vol_open - vol_below(K), 0.0) + elseif (L(K) <= L_2_reg) then + vol_err(K) = 0.5*L(K)**2 * (slope + crv*(1.0 - 4.0*C1_3*L(K))) - vol_below(K) + else ! There are two separate open regions. + Vol_err(K) = crv_3 * (L(K)**2 * ( 0.75 - 0.5*L(K))) + (slope2_4crv - vol_below(K)) + endif + enddo ! k loop to determine L(K) in the concave case + +end subroutine test_L_open_concave + + !> Determine the normalized open length of each interface for convex bathymetry (from the ocean !! perspective) using Newton's method iterations. In this case there is a single open region !! with the minimum depth at one edge of the cell. From 714d2da17ed29d34f73e88513192acb9439e1061 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Tue, 5 Mar 2024 07:10:00 -0500 Subject: [PATCH 08/15] Minor refactoring in set_viscous_BBL Carried out minor refactoring in set_viscous_BBL as suggested by the reviews of this PR, including the elimination of some unnecessary error handling and the replacement of C2pi_3 as an argument to find_L_open_concave_trigonometric with an internal parameter. All answers are bitwise identical and there are no changes to publicly visible interfaces. --- .../vertical/MOM_set_viscosity.F90 | 29 +++++++------------ 1 file changed, 11 insertions(+), 18 deletions(-) diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index 1778b8d870..60298f9226 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -260,7 +260,6 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) real :: Dp, Dm ! The depths at the edges of a velocity cell [Z ~> m]. real :: crv ! crv is the curvature of the bottom depth across a ! cell, times the cell width squared [Z ~> m]. - real :: crv_3 ! crv/3 [Z ~> m]. real :: slope ! The absolute value of the bottom depth slope across ! a cell times the cell width [Z ~> m]. real :: Vol_bbl_chan ! The volume of the bottom boundary layer as used in the channel @@ -304,7 +303,6 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) real :: h_sum ! The sum of the thicknesses of the layers below the one being ! worked on [H ~> m or kg m-2]. real, parameter :: C1_3 = 1.0/3.0, C1_6 = 1.0/6.0, C1_12 = 1.0/12.0 ! Rational constants [nondim] - real :: C2pi_3 ! An irrational constant, 2/3 pi. [nondim] real :: tmp ! A temporary variable, sometimes in [Z ~> m] logical :: use_BBL_EOS, do_i(SZIB_(G)) integer, dimension(2) :: EOSdom ! The computational domain for the equation of state @@ -318,7 +316,6 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) dz_neglect = GV%dZ_subroundoff Rho0x400_G = 400.0*(GV%H_to_RZ / (US%L_to_Z**2 * GV%g_Earth)) - C2pi_3 = 8.0*atan(1.0)/3.0 if (.not.CS%initialized) call MOM_error(FATAL,"MOM_set_viscosity(BBL): "//& "Module must be initialized before it is used.") @@ -424,7 +421,7 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) if (allocated(visc%Ray_v)) visc%Ray_v(:,:,:) = 0.0 !$OMP parallel do default(private) shared(u,v,h,dz,tv,visc,G,GV,US,CS,Rml,nz,nkmb,nkml,K2, & - !$OMP Isq,Ieq,Jsq,Jeq,h_neglect,dz_neglect,Rho0x400_G,C2pi_3, & + !$OMP Isq,Ieq,Jsq,Jeq,h_neglect,dz_neglect,Rho0x400_G, & !$OMP U_bg_sq,cdrag_sqrt,cdrag_sqrt_H,cdrag_sqrt_H_RL, & !$OMP cdrag_L_to_H,cdrag_RL_to_H,use_BBL_EOS,BBL_thick_max, & !$OMP OBC,D_u,D_v,mask_u,mask_v,pbv) @@ -870,7 +867,7 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) Dm = 2.0 * D_vel * tmp / (D_vel + tmp) endif if (Dm > Dp) then ; tmp = Dp ; Dp = Dm ; Dm = tmp ; endif - crv_3 = (Dp + Dm - 2.0*D_vel) ; crv = 3.0*crv_3 + crv = 3.0*(Dp + Dm - 2.0*D_vel) slope = Dp - Dm ! If the curvature is small enough, there is no reason not to assume @@ -882,15 +879,15 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) call find_L_open_uniform_slope(vol_below, Dp, Dm, L, GV) elseif (crv > 0.0) then if (CS%concave_trigonometric_L) then - call find_L_open_concave_trigonometric(vol_below, D_vel, Dp, Dm, C2pi_3, L, GV) + call find_L_open_concave_trigonometric(vol_below, D_vel, Dp, Dm, L, GV) else call find_L_open_concave_iterative(vol_below, D_vel, Dp, Dm, L, GV) if (CS%debug) then ! The tests in this block reveal that the iterative and trigonometric solutions are - ! mathematically equiavalent, but in some cases the iterative solution is consistent + ! mathematically equivalent, but in some cases the iterative solution is consistent ! at roundoff, but that the trigonmetric solutions have errors that can be several ! orders of magnitude larger in some cases. - call find_L_open_concave_trigonometric(vol_below, D_vel, Dp, Dm, C2pi_3, L_trig, GV) + call find_L_open_concave_trigonometric(vol_below, D_vel, Dp, Dm, L_trig, GV) call test_L_open_concave(vol_below, D_vel, Dp, Dm, L_trig, vol_err_trig, GV) call test_L_open_concave(vol_below, D_vel, Dp, Dm, L, vol_err_iter, GV) max_dL_trig_itt = 0.0 ; max_norm_err_trig = 0.0 ; max_norm_err_iter = 0.0 @@ -1128,7 +1125,7 @@ end subroutine find_L_open_uniform_slope !> Determine the normalized open length of each interface for concave bathymetry (from the ocean perspective) !! using trigonometric expressions. In this case there can be two separate open regions. -subroutine find_L_open_concave_trigonometric(vol_below, D_vel, Dp, Dm, C2pi_3, L, GV) +subroutine find_L_open_concave_trigonometric(vol_below, D_vel, Dp, Dm, L, GV) type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. real, dimension(SZK_(GV)+1), intent(in) :: vol_below !< The volume below each interface, normalized by !! the full horizontal area of a velocity cell [Z ~> m] @@ -1137,7 +1134,6 @@ subroutine find_L_open_concave_trigonometric(vol_below, D_vel, Dp, Dm, C2pi_3, L !! of a velocity cell [Z ~> m] real, intent(in) :: Dm !< The smaller of the two depths at the edge !! of a velocity cell [Z ~> m] - real, intent(in) :: C2pi_3 !< An irrational constant, 2/3 pi [nondim] real, dimension(SZK_(GV)+1), intent(out) :: L !< The fraction of the full cell width that is open at !! the depth of each interface [nondim] @@ -1159,6 +1155,7 @@ subroutine find_L_open_concave_trigonometric(vol_below, D_vel, Dp, Dm, C2pi_3, L real :: slope_crv ! The slope divided by the curvature [nondim] real :: tmp_val_m1_to_p1 ! A temporary variable [nondim] real, parameter :: C1_3 = 1.0/3.0, C1_12 = 1.0/12.0 ! Rational constants [nondim] + real, parameter :: C2pi_3 = 8.0*atan(1.0)/3.0 ! An irrational constant, 2/3 pi. [nondim] integer :: K, nz nz = GV%ke @@ -1166,7 +1163,6 @@ subroutine find_L_open_concave_trigonometric(vol_below, D_vel, Dp, Dm, C2pi_3, L ! Each cell extends from x=-1/2 to 1/2, and has a topography ! given by D(x) = crv*x^2 + slope*x + D_vel - crv/12. crv_3 = (Dp + Dm - 2.0*D_vel) ; crv = 3.0*crv_3 - if (crv <= 0.0) call MOM_error(FATAL, "find_L_open_concave should only be called with a positive curvature.") slope = Dp - Dm ! Calculate the volume above which the entire cell is open and the volume at which the @@ -1284,7 +1280,6 @@ subroutine find_L_open_concave_iterative(vol_below, D_vel, Dp, Dm, L, GV) ! given by D(x) = crv*x^2 + slope*x + D_vel - crv/12. crv_3 = (Dp + Dm - 2.0*D_vel) ; crv = 3.0*crv_3 - if (crv <= 0.0) call MOM_error(FATAL, "find_L_open_concave should only be called with a positive curvature.") slope = Dp - Dm ! Calculate the volume above which the entire cell is open and the volume at which the @@ -1424,8 +1419,8 @@ subroutine find_L_open_concave_iterative(vol_below, D_vel, Dp, Dm, L, GV) endif endif Vol_err_max = 0.5*L_max**2 * (slope + crv*(1.0 - 4.0*C1_3*L_max)) - vol_below(K) - if (Vol_err_max < 0.0) call MOM_error(FATAL, & - "Vol_err_max should never be negative in find_L_open_concave_iterative.") + ! if (Vol_err_max < 0.0) call MOM_error(FATAL, & + ! "Vol_err_max should never be negative in find_L_open_concave_iterative.") if ((Vol_err_max < abs(Vol_err)) .and. (L_max < 1.0)) then ! Start with 1 bounded Newton's method step from L_max dVol_dL = L_max * (slope + crv*(1.0 - 2.0*L_max)) @@ -1516,8 +1511,8 @@ subroutine find_L_open_concave_iterative(vol_below, D_vel, Dp, Dm, L, GV) ! This over-estimate of L(K) is accurate for L ~= 1: L_max = 1.0 - sqrt( (Vol_open - vol_below(K)) * C4_crv ) Vol_err_max = crv_3 * (L_max**2 * ( 0.75 - 0.5*L_max)) + (slope2_4crv - vol_below(K)) - if (Vol_err_max < 0.0) call MOM_error(FATAL, & - "Vol_err_max should never be negative in find_L_open_concave_iterative.") + ! if (Vol_err_max < 0.0) call MOM_error(FATAL, & + ! "Vol_err_max should never be negative in find_L_open_concave_iterative.") if ((Vol_err_max < abs(Vol_err)) .and. (L_max < 1.0)) then ! Start with 1 bounded Newton's method step from L_max dVol_dL = 0.5*crv * (L_max * (1.0 - L_max)) @@ -1587,7 +1582,6 @@ subroutine test_L_open_concave(vol_below, D_vel, Dp, Dm, L, vol_err, GV) ! given by D(x) = crv*x^2 + slope*x + D_vel - crv/12. crv_3 = (Dp + Dm - 2.0*D_vel) ; crv = 3.0*crv_3 - if (crv <= 0.0) call MOM_error(FATAL, "test_L_open_concave should only be called with a positive curvature.") slope = Dp - Dm ! Calculate the volume above which the entire cell is open and the volume at which the @@ -1678,7 +1672,6 @@ subroutine find_L_open_convex(vol_below, D_vel, Dp, Dm, L, GV, US, CS) ! Each cell extends from x=-1/2 to 1/2, and has a topography ! given by D(x) = crv*x^2 + slope*x + D_vel - crv/12. crv_3 = (Dp + Dm - 2.0*D_vel) ; crv = 3.0*crv_3 - if (crv >= 0.0) call MOM_error(FATAL, "find_L_open_convex should only be called with a negative curvature.") slope = Dp - Dm ! Calculate the volume above which the entire cell is open and the volume at which the From 6153d97ad5d6ed0c81e95382710cac313ef64973 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Sat, 9 Mar 2024 10:50:07 -0500 Subject: [PATCH 09/15] Compute dz in the halo for OBC segments It was blowing up with "forrtl: error (65): floating invalid" when accessing dz in the halo at the boundary, but just sometimes. My default layout is trouble while my testing layout of 48 cores is not. --- src/core/MOM.F90 | 5 ++++- src/tracer/MOM_tracer_advect.F90 | 3 ++- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 9d64b9bf17..53f20bb10b 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1436,8 +1436,11 @@ subroutine step_MOM_tracer_dyn(CS, G, GV, US, h, Time_local) CS%tracer_diff_CSp, CS%tracer_Reg, CS%tv) if (CS%debug) call MOM_tracer_chksum("Post-diffuse ", CS%tracer_Reg, G) if (showCallTree) call callTree_waypoint("finished tracer advection/diffusion (step_MOM)") - call update_segment_tracer_reservoirs(G, GV, CS%uhtr, CS%vhtr, h, CS%OBC, & + if (associated(CS%OBC)) then + call pass_vector(CS%uhtr, CS%vhtr, G%Domain) + call update_segment_tracer_reservoirs(G, GV, CS%uhtr, CS%vhtr, h, CS%OBC, & CS%t_dyn_rel_adv, CS%tracer_Reg) + endif call cpu_clock_end(id_clock_tracer) ; call cpu_clock_end(id_clock_thermo) call cpu_clock_begin(id_clock_other) ; call cpu_clock_begin(id_clock_diagnostics) diff --git a/src/tracer/MOM_tracer_advect.F90 b/src/tracer/MOM_tracer_advect.F90 index efe6397de0..ef2c3125cd 100644 --- a/src/tracer/MOM_tracer_advect.F90 +++ b/src/tracer/MOM_tracer_advect.F90 @@ -123,7 +123,8 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, x_first x_first = (MOD(G%first_direction,2) == 0) ! increase stencil size for Colella & Woodward PPM - if (CS%usePPM .and. .not. CS%useHuynh) stencil = 3 +! if (CS%usePPM .and. .not. CS%useHuynh) stencil = 3 + if (CS%usePPM) stencil = 3 ntr = Reg%ntr Idt = 1.0 / dt From 0ff03aea4badc0634d0f024e96d1385255724e29 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Tue, 6 Feb 2024 16:32:58 -0500 Subject: [PATCH 10/15] +Fix myStats global indexing segmentation faults Revised the interfaces to the myStats routine in the horizontal_regridding module to avoid segmentation faults due to inconsistent horizontal indices and array extents in global indexing mode. Rather than passing in absolute array extents to work on, an ocean grid type is now passed as an argument to myStats, with the new optional full_halo argument used to capture the case where the tracer statistics are being taken over the full data domain. The most frequently encountered problems occurred when the hard-coded debug variable in the horiz_interp_and_extrap_tracer routines are changed from false to true. When global indexing is not used, this revised work exactly as before, but when it is used with global indexing, it avoids segmentation faults that were preventing the model from running in some cases with all debugging enabled. --- src/framework/MOM_horizontal_regridding.F90 | 44 +++++++++++-------- .../MOM_tracer_initialization_from_Z.F90 | 2 +- 2 files changed, 26 insertions(+), 20 deletions(-) diff --git a/src/framework/MOM_horizontal_regridding.F90 b/src/framework/MOM_horizontal_regridding.F90 index 34d0b73cb9..205dd6d7be 100644 --- a/src/framework/MOM_horizontal_regridding.F90 +++ b/src/framework/MOM_horizontal_regridding.F90 @@ -44,27 +44,32 @@ module MOM_horizontal_regridding contains !> Write to the terminal some basic statistics about the k-th level of an array -subroutine myStats(array, missing, is, ie, js, je, k, mesg, scale) - real, dimension(:,:), intent(in) :: array !< input array in arbitrary units [A ~> a] - real, intent(in) :: missing !< missing value in arbitrary units [A ~> a] - integer, intent(in) :: is !< Start index in i - integer, intent(in) :: ie !< End index in i - integer, intent(in) :: js !< Start index in j - integer, intent(in) :: je !< End index in j - integer, intent(in) :: k !< Level to calculate statistics for - character(len=*), intent(in) :: mesg !< Label to use in message - real, optional, intent(in) :: scale !< A scaling factor for output [a A-1 ~> 1] +subroutine myStats(array, missing, G, k, mesg, scale, full_halo) + type(ocean_grid_type), intent(in) :: G !< Ocean grid type + real, dimension(SZI_(G),SZJ_(G)), & + intent(in) :: array !< input array in arbitrary units [A ~> a] + real, intent(in) :: missing !< missing value in arbitrary units [A ~> a] + integer, intent(in) :: k !< Level to calculate statistics for + character(len=*), intent(in) :: mesg !< Label to use in message + real, optional, intent(in) :: scale !< A scaling factor for output [a A-1 ~> 1] + logical, optional, intent(in) :: full_halo !< If present and true, test values on the whole + !! array rather than just the computational domain. ! Local variables real :: minA ! Minimum value in the array in the arbitrary units of the input array [A ~> a] real :: maxA ! Maximum value in the array in the arbitrary units of the input array [A ~> a] real :: scl ! A factor for undoing any scaling of the array statistics for output [a A-1 ~> 1] - integer :: i,j + integer :: i, j, is, ie, js, je logical :: found character(len=120) :: lMesg scl = 1.0 ; if (present(scale)) scl = scale minA = 9.E24 / scl ; maxA = -9.E24 / scl ; found = .false. + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + if (present(full_halo)) then ; if (full_halo) then + is = G%isd ; ie = G%ied ; js = G%jsd ; je = G%jed + endif ; endif + do j=js,je ; do i=is,ie if (array(i,j) /= array(i,j)) stop 'Nan!' if (abs(array(i,j)-missing) > 1.e-6*abs(missing)) then @@ -309,7 +314,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, recnum, G, tr real, dimension(:,:), allocatable :: tr_in !< A 2-d array for holding input data on its !! native horizontal grid, with units that change !! as the input data is interpreted [a] then [A ~> a] - real, dimension(:,:,:), allocatable :: tr_in_full !< A 3-d array for holding input data on the + real, dimension(:,:,:), allocatable :: tr_in_full !< A 3-d array for holding input data on the !! model horizontal grid, with units that change !! as the input data is interpreted [a] then [A ~> a] real, dimension(:,:), allocatable :: tr_inp !< Native horizontal grid data extended to the poles @@ -332,7 +337,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, recnum, G, tr real :: npole ! The number of points contributing to the pole value [nondim] real :: missing_val_in ! The missing value in the input field [a] real :: roundoff ! The magnitude of roundoff, usually ~2e-16 [nondim] - real :: add_offset, scale_factor ! File-specific conversion factors. + real :: add_offset, scale_factor ! File-specific conversion factors [a] or [nondim] integer :: ans_date ! The vintage of the expressions and order of arithmetic to use logical :: found_attr logical :: add_np @@ -545,7 +550,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, recnum, G, tr endif if (debug) then - call myStats(tr_inp, missing_value, 1, id, 1, jd, k, 'Tracer from file', scale=I_scale) + call myStats(tr_inp, missing_value, G, k, 'Tracer from file', scale=I_scale, full_halo=.true.) endif call run_horiz_interp(Interp, tr_inp, tr_out(is:ie,js:je), missing_value=missing_value) @@ -568,11 +573,12 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, recnum, G, tr (mask_out(i,j) < 1.0)) & fill(i,j) = 1.0 enddo ; enddo + call pass_var(fill, G%Domain) call pass_var(good, G%Domain) if (debug) then - call myStats(tr_out, missing_value, is, ie, js, je, k, 'variable from horiz_interp()', scale=I_scale) + call myStats(tr_out, missing_value, G, k, 'variable from horiz_interp()', scale=I_scale) endif ! Horizontally homogenize data to produce perfectly "flat" initial conditions @@ -589,7 +595,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, recnum, G, tr call fill_miss_2d(tr_outf, good2, fill2, tr_prev, G, dtr_iter_stop, answer_date=ans_date) if (debug) then - call myStats(tr_outf, missing_value, is, ie, js, je, k, 'field from fill_miss_2d()', scale=I_scale) + call myStats(tr_outf, missing_value, G, k, 'field from fill_miss_2d()', scale=I_scale) endif tr_z(:,:,k) = tr_outf(:,:) * G%mask2dT(:,:) @@ -850,7 +856,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(field, Time, G, tr_z, mask_z, & endif if (debug) then - call myStats(tr_inp, missing_value, 1, id, 1, jd, k, 'Tracer from file', scale=I_scale) + call myStats(tr_inp, missing_value, G, k, 'Tracer from file', scale=I_scale, full_halo=.true.) endif tr_out(:,:) = 0.0 @@ -878,7 +884,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(field, Time, G, tr_z, mask_z, & call pass_var(good, G%Domain) if (debug) then - call myStats(tr_out, missing_value, is, ie, js, je, k, 'variable from horiz_interp()', scale=I_scale) + call myStats(tr_out, missing_value, G, k, 'variable from horiz_interp()', scale=I_scale) endif ! Horizontally homogenize data to produce perfectly "flat" initial conditions @@ -897,7 +903,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(field, Time, G, tr_z, mask_z, & ! if (debug) then ! call hchksum(tr_outf, 'field from fill_miss_2d ', G%HI, scale=I_scale) -! call myStats(tr_outf, missing_value, is, ie, js, je, k, 'field from fill_miss_2d()', scale=I_scale) +! call myStats(tr_outf, missing_value, G, k, 'field from fill_miss_2d()', scale=I_scale) ! endif tr_z(:,:,k) = tr_outf(:,:) * G%mask2dT(:,:) diff --git a/src/initialization/MOM_tracer_initialization_from_Z.F90 b/src/initialization/MOM_tracer_initialization_from_Z.F90 index 5a172b5d97..d28a925c03 100644 --- a/src/initialization/MOM_tracer_initialization_from_Z.F90 +++ b/src/initialization/MOM_tracer_initialization_from_Z.F90 @@ -217,7 +217,7 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_ deallocate( h1 ) do k=1,nz - call myStats(tr(:,:,k), missing_value, is, ie, js, je, k, 'Tracer from ALE()') + call myStats(tr(:,:,k), missing_value, G, k, 'Tracer from ALE()') enddo call cpu_clock_end(id_clock_ALE) endif ! useALEremapping From 0bd4c160b2d42564effe11918134f597dfb05634 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 7 Feb 2024 20:46:41 -0500 Subject: [PATCH 11/15] (*)Fix apply_topography_edits_from_file indexing Corrected bugs in the horizontal indexing in apply_topography_edits_from_file that led to differing answers depending on the value of GLOBAL_INDEXING. This change gives identical results when GLOBAL_INDEXING is used, as can be seen by noting that G%idg_offset = G%isd_global + G%isd and that without global indexing G%isd = 1, but with global indexing the new expressions give the same answers as without them. Because global indexing is not typically used, answers are not changed for any cases in the MOM6-examples test suite. --- src/initialization/MOM_shared_initialization.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/initialization/MOM_shared_initialization.F90 b/src/initialization/MOM_shared_initialization.F90 index 46d0448699..821232b80d 100644 --- a/src/initialization/MOM_shared_initialization.F90 +++ b/src/initialization/MOM_shared_initialization.F90 @@ -256,8 +256,8 @@ subroutine apply_topography_edits_from_file(D, G, param_file, US) call close_file_to_read(ncid, topo_edits_file) do n = 1, n_edits - i = ig(n) - G%isd_global + 2 ! +1 for python indexing and +1 for ig-isd_global+1 - j = jg(n) - G%jsd_global + 2 + i = ig(n) - G%idg_offset + 1 ! +1 for python indexing + j = jg(n) - G%jdg_offset + 1 if (i>=G%isc .and. i<=G%iec .and. j>=G%jsc .and. j<=G%jec) then if (new_depth(n) /= mask_depth) then write(stdout,'(a,3i5,f8.2,a,f8.2,2i4)') & From 9e78615cb396250120bde2a6e6bb4426cd58bbbd Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 4 Feb 2024 09:08:35 -0500 Subject: [PATCH 12/15] Earlier initialization thickness halo updates Moved the post-initialization halo updates for thicknesses and temperatures and salinities in initialize_MOM to occur immediately after the last point where they are modified and before the remapped diagnostic grids are set up. This does not change any existing answers but it will enable the future use of the thermo_var_ptr type in calls to thickness_to_dz when setting up remapped diagnostics, and perhaps elsewhere in the initialization of other auxiliary variables. All answers are bitwise identical. --- src/core/MOM.F90 | 43 ++++++++++++++++++++++++------------------- 1 file changed, 24 insertions(+), 19 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 53f20bb10b..2bf2b11b1a 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -3117,6 +3117,29 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & endif if ( CS%use_ALE_algorithm ) call ALE_updateVerticalGridType( CS%ALE_CSp, GV ) + ! The basic state variables have now been fully initialized, so update their halos and + ! calculate any derived thermodynmics quantities. + + !--- set up group pass for u,v,T,S and h. pass_uv_T_S_h also is used in step_MOM + call cpu_clock_begin(id_clock_pass_init) + dynamics_stencil = min(3, G%Domain%nihalo, G%Domain%njhalo) + call create_group_pass(pass_uv_T_S_h, CS%u, CS%v, G%Domain, halo=dynamics_stencil) + if (use_temperature) then + call create_group_pass(pass_uv_T_S_h, CS%tv%T, G%Domain, halo=dynamics_stencil) + call create_group_pass(pass_uv_T_S_h, CS%tv%S, G%Domain, halo=dynamics_stencil) + endif + call create_group_pass(pass_uv_T_S_h, CS%h, G%Domain, halo=dynamics_stencil) + + call do_group_pass(pass_uv_T_S_h, G%Domain) + if (associated(CS%tv%p_surf)) call pass_var(CS%tv%p_surf, G%Domain, halo=dynamics_stencil) + call cpu_clock_end(id_clock_pass_init) + + ! Update derived thermodynamic quantities. + if (allocated(CS%tv%SpV_avg)) then + call calc_derived_thermo(CS%tv, CS%h, G, GV, US, halo=dynamics_stencil, debug=CS%debug) + endif + + diag => CS%diag ! Initialize the diag mediator. call diag_mediator_init(G, GV, US, GV%ke, param_file, diag, doc_file_dir=dirs%output_directory) @@ -3137,7 +3160,6 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & ! Whenever thickness/T/S changes let the diag manager know, target grids ! for vertical remapping may need to be regenerated. - ! FIXME: are h, T, S updated at the same time? Review these for T, S updates. call diag_update_remap_grids(diag) ! Setup the diagnostic grid storage types @@ -3287,30 +3309,13 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & call ALE_register_diags(Time, G, GV, US, diag, CS%ALE_CSp) endif - !--- set up group pass for u,v,T,S and h. pass_uv_T_S_h also is used in step_MOM + ! Do any necessary halo updates on any auxiliary variables that have been initialized. call cpu_clock_begin(id_clock_pass_init) - dynamics_stencil = min(3, G%Domain%nihalo, G%Domain%njhalo) - call create_group_pass(pass_uv_T_S_h, CS%u, CS%v, G%Domain, halo=dynamics_stencil) - if (use_temperature) then - call create_group_pass(pass_uv_T_S_h, CS%tv%T, G%Domain, halo=dynamics_stencil) - call create_group_pass(pass_uv_T_S_h, CS%tv%S, G%Domain, halo=dynamics_stencil) - endif - call create_group_pass(pass_uv_T_S_h, CS%h, G%Domain, halo=dynamics_stencil) - - call do_group_pass(pass_uv_T_S_h, G%Domain) - - ! Update derived thermodynamic quantities. - if (associated(CS%tv%p_surf)) call pass_var(CS%tv%p_surf, G%Domain, halo=dynamics_stencil) - if (allocated(CS%tv%SpV_avg)) then - call calc_derived_thermo(CS%tv, CS%h, G, GV, US, halo=dynamics_stencil, debug=CS%debug) - endif - if (associated(CS%visc%Kv_shear)) & call pass_var(CS%visc%Kv_shear, G%Domain, To_All+Omit_Corners, halo=1) if (associated(CS%visc%Kv_slow)) & call pass_var(CS%visc%Kv_slow, G%Domain, To_All+Omit_Corners, halo=1) - call cpu_clock_end(id_clock_pass_init) ! This subroutine initializes any tracer packages. From 1ae2a3a35ecf9d26a0b4550850d310a9440d6c38 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 5 Feb 2024 07:57:17 -0500 Subject: [PATCH 13/15] +(*)Refactor MOM_diag_remap for global indices Refactored MOM_diag_remap to work with global indices and to move logical branches outside of do loops. This was done by adding internal routines that set the loop indices consistently with the NE c-grid convention used throughout the MOM6 code and converting the optionally associated mask pointer into an optional argument. It also simplifies the logic of many of the expressions within the remapping code. There is also a new element, Z_based_coord, in the diag_remap_ctrl type, to indicate whether the remapping is working in thickness or height units, but for now it is always set to false. The function set_h_neglect or set_dz_neglect is used to set the negligible thicknesses used for remapping in diag_remap_update and diag_remap_do_remap, depending on whether the remapping is being done in thickness or vertical height coordinates. Diag_remap_init has a new vertical_grid_type argument, and diag_remap_do_remap has a new unit_scale_type argument. For REMAPPING_ANSWER_DATES later than 20240201, diag_remap_updated does an explicit sum to determine the total water column thickness, rather than using sum function, which is indeterminate of the order of the sums. For some compilers, this could change the vertical grids used for remapping diagnostics at roundoff, but no such change was detected for any of the compilers used with the MOM6 regression test suite. All answers and diagnostics in cases that worked before are bitwise identical, but there are new arguments to two publicly visible interfaces. --- src/framework/MOM_diag_mediator.F90 | 4 +- src/framework/MOM_diag_remap.F90 | 720 +++++++++++++++++----------- 2 files changed, 439 insertions(+), 285 deletions(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 2c71a93e42..61bfa93046 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -1583,7 +1583,7 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask, alt_h) if (id_clock_diag_remap>0) call cpu_clock_begin(id_clock_diag_remap) allocate(remapped_field(size(field,1), size(field,2), diag%axes%nz)) call diag_remap_do_remap(diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number), & - diag_cs%G, diag_cs%GV, h_diag, staggered_in_x, staggered_in_y, & + diag_cs%G, diag_cs%GV, diag_cs%US, h_diag, staggered_in_x, staggered_in_y, & diag%axes%mask3d, field, remapped_field) if (id_clock_diag_remap>0) call cpu_clock_end(id_clock_diag_remap) if (associated(diag%axes%mask3d)) then @@ -3207,7 +3207,7 @@ subroutine diag_mediator_init(G, GV, US, nz, param_file, diag_cs, doc_file_dir) allocate(diag_cs%diag_remap_cs(diag_cs%num_diag_coords)) ! Initialize each diagnostic vertical coordinate do i=1, diag_cs%num_diag_coords - call diag_remap_init(diag_cs%diag_remap_cs(i), diag_coords(i), answer_date=remap_answer_date) + call diag_remap_init(diag_cs%diag_remap_cs(i), diag_coords(i), answer_date=remap_answer_date, GV=GV) enddo deallocate(diag_coords) endif diff --git a/src/framework/MOM_diag_remap.F90 b/src/framework/MOM_diag_remap.F90 index 0e2e594403..cd9e22aae2 100644 --- a/src/framework/MOM_diag_remap.F90 +++ b/src/framework/MOM_diag_remap.F90 @@ -15,42 +15,14 @@ !! the diagnostic is written out. -! NOTE: In the following functions, the fields are passed using 1-based -! indexing, which requires special handling within the grid index loops. +! NOTE: In the following functions, the fields are initially passed using 1-based +! indexing, which are then passed to separate private internal routines that shift +! the indexing to use the same indexing conventions used elsewhere in the MOM6 code. ! -! * diag_remap_do_remap -! * vertically_reintegrate_diag_field -! * vertically_interpolate_diag_field -! * horizontally_average_diag_field -! -! Symmetric grids add an additional row of western and southern points to u- -! and v-grids. Non-symmetric grids are 1-based and symmetric grids are -! zero-based, allowing the same expressions to be used when accessing the -! fields. But if u- or v-points become 1-indexed, as in these functions, then -! the stencils must be re-assessed. -! -! For interpolation between h and u grids, we use the following relations: -! -! h->u: f_u(ig) = 0.5 * (f_h( ig ) + f_h(ig+1)) -! f_u(i1) = 0.5 * (f_h(i1-1) + f_h( i1 )) -! -! u->h: f_h(ig) = 0.5 * (f_u(ig-1) + f_u( ig )) -! f_h(i1) = 0.5 * (f_u( i1 ) + f_u(i1+1)) -! -! where ig is the grid index and i1 is the 1-based index. That is, a 1-based -! u-point is ahead of its matching h-point in non-symmetric mode, but behind -! its matching h-point in non-symmetric mode. -! -! We can combine these expressions by applying to ig a -1 shift on u-grids and -! a +1 shift on h-grids in symmetric mode. -! -! We do not adjust the h-point indices, since they are assumed to be 1-based. -! This is only correct when global indexing is disabled. If global indexing is -! enabled, then all indices will need to be defined relative to the data -! domain. -! -! Finally, note that the mask input fields are pointers to arrays which are -! zero-indexed, and do not need any corrections over grid index loops. +! * diag_remap_do_remap, which calls do_remap +! * vertically_reintegrate_diag_field, which calls vertically_reintegrate_field +! * vertically_interpolate_diag_field, which calls vertically_interpolate_field +! * horizontally_average_diag_field, which calls horizontally_average_field module MOM_diag_remap @@ -70,10 +42,9 @@ module MOM_diag_remap use MOM_EOS, only : EOS_type use MOM_remapping, only : remapping_CS, initialize_remapping, remapping_core_h use MOM_remapping, only : interpolate_column, reintegrate_column -use MOM_regridding, only : regridding_CS, initialize_regridding -use MOM_regridding, only : end_regridding +use MOM_regridding, only : regridding_CS, initialize_regridding, end_regridding use MOM_regridding, only : set_regrid_params, get_regrid_size -use MOM_regridding, only : getCoordinateInterfaces +use MOM_regridding, only : getCoordinateInterfaces, set_h_neglect, set_dz_neglect use MOM_regridding, only : get_zlike_CS, get_sigma_CS, get_rho_CS use regrid_consts, only : coordinateMode use coord_zlike, only : build_zstar_column @@ -83,6 +54,8 @@ module MOM_diag_remap implicit none ; private +#include "MOM_memory.h" + public diag_remap_ctrl public diag_remap_init, diag_remap_end, diag_remap_update, diag_remap_do_remap public diag_remap_configure_axes, diag_remap_axes_configured @@ -104,14 +77,19 @@ module MOM_diag_remap logical :: used = .false. !< Whether this coordinate actually gets used. integer :: vertical_coord = 0 !< The vertical coordinate that we remap to character(len=10) :: vertical_coord_name ='' !< The coordinate name as understood by ALE + logical :: Z_based_coord = .false. !< If true, this coordinate is based on remapping of + !! geometric distances across layers (in [Z ~> m]) rather + !! than layer thicknesses (in [H ~> m or kg m-2]). This + !! distinction only matters in non-Boussinesq mode. character(len=16) :: diag_coord_name = '' !< A name for the purpose of run-time parameters character(len=8) :: diag_module_suffix = '' !< The suffix for the module to appear in diag_table type(remapping_CS) :: remap_cs !< Remapping control structure use for this axes type(regridding_CS) :: regrid_cs !< Regridding control structure that defines the coordinates for this axes integer :: nz = 0 !< Number of vertical levels used for remapping - real, dimension(:,:,:), allocatable :: h !< Remap grid thicknesses [H ~> m or kg m-2] - real, dimension(:,:,:), allocatable :: h_extensive !< Remap grid thicknesses for extensive - !! variables [H ~> m or kg m-2] + real, dimension(:,:,:), allocatable :: h !< Remap grid thicknesses in [H ~> m or kg m-2] or + !! vertical extents in [Z ~> m], depending on the setting of Z_based_coord. + real, dimension(:,:,:), allocatable :: h_extensive !< Remap grid thicknesses in [H ~> m or kg m-2] or + !! vertical extents in [Z ~> m] for remapping extensive variables integer :: interface_axes_id = 0 !< Vertical axes id for remapping at interfaces integer :: layer_axes_id = 0 !< Vertical axes id for remapping on layers integer :: answer_date !< The vintage of the order of arithmetic and expressions @@ -124,7 +102,7 @@ module MOM_diag_remap contains !> Initialize a diagnostic remapping type with the given vertical coordinate. -subroutine diag_remap_init(remap_cs, coord_tuple, answer_date) +subroutine diag_remap_init(remap_cs, coord_tuple, answer_date, GV) type(diag_remap_ctrl), intent(inout) :: remap_cs !< Diag remapping control structure character(len=*), intent(in) :: coord_tuple !< A string in form of !! MODULE_SUFFIX PARAMETER_SUFFIX COORDINATE_NAME @@ -132,11 +110,19 @@ subroutine diag_remap_init(remap_cs, coord_tuple, answer_date) !! to use for remapping. Values below 20190101 recover !! the answers from 2018, while higher values use more !! robust forms of the same remapping expressions. + type(verticalGrid_type), intent(in) :: GV !< The ocean vertical grid structure, used here to evaluate + !! whether the model is in non-Boussinesq mode. remap_cs%diag_module_suffix = trim(extractWord(coord_tuple, 1)) remap_cs%diag_coord_name = trim(extractWord(coord_tuple, 2)) remap_cs%vertical_coord_name = trim(extractWord(coord_tuple, 3)) remap_cs%vertical_coord = coordinateMode(remap_cs%vertical_coord_name) + remap_cs%Z_based_coord = .false. + ! if ( & ! (.not.(GV%Boussinesq .or. GV%semi_Boussinesq)) .and. & + ! ((remap_cs%vertical_coord == coordinateMode('ZSTAR')) .or. & + ! (remap_cs%vertical_coord == coordinateMode('SIGMA'))) ) & + ! remap_cs%Z_based_coord = .true. + remap_cs%configured = .false. remap_cs%initialized = .false. remap_cs%used = .false. @@ -274,31 +260,46 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state, h_targe type(ocean_grid_type), pointer :: G !< The ocean's grid type type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - real, dimension(:,:,:), intent(in) :: h !< New thickness [H ~> m or kg m-2] - real, dimension(:,:,:), intent(in) :: T !< New temperatures [C ~> degC] - real, dimension(:,:,:), intent(in) :: S !< New salinities [S ~> ppt] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: h !< New thickness in [H ~> m or kg m-2] or [Z ~> m], depending + !! on the value of remap_cs%Z_based_coord + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: T !< New temperatures [C ~> degC] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: S !< New salinities [S ~> ppt] type(EOS_type), intent(in) :: eqn_of_state !< A pointer to the equation of state - real, dimension(:,:,:), intent(inout) :: h_target !< The new diagnostic thicknesses [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),remap_cs%nz), & + intent(inout) :: h_target !< The new diagnostic thicknesses in [H ~> m or kg m-2] + !! or [Z ~> m], depending on the value of remap_cs%Z_based_coord ! Local variables - real, dimension(remap_cs%nz + 1) :: zInterfaces ! Interface positions [H ~> m or kg m-2] - real :: h_neglect, h_neglect_edge ! Negligible thicknesses [H ~> m or kg m-2] - integer :: i, j, k, nz - - ! Note that coordinateMode('LAYER') is never 'configured' so will - ! always return here. - if (.not. remap_cs%configured) then - return - endif - - if (remap_cs%answer_date >= 20190101) then - h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff - elseif (GV%Boussinesq) then - h_neglect = GV%m_to_H*1.0e-30 ; h_neglect_edge = GV%m_to_H*1.0e-10 + real, dimension(remap_cs%nz + 1) :: zInterfaces ! Interface positions [H ~> m or kg m-2] or [Z ~> m] + real :: h_neglect, h_neglect_edge ! Negligible thicknesses [H ~> m or kg m-2] or [Z ~> m] + real :: bottom_depth(SZI_(G),SZJ_(G)) ! The depth of the bathymetry in [H ~> m or kg m-2] or [Z ~> m] + real :: h_tot(SZI_(G),SZJ_(G)) ! The total thickness of the water column [H ~> m or kg m-2] or [Z ~> m] + real :: Z_unit_scale ! A conversion factor from Z-units the internal work units in this routine, + ! in units of [H Z-1 ~> 1 or kg m-3] or [nondim], depending on remap_cs%Z_based_coord. + integer :: i, j, k, is, ie, js, je, nz + + ! Note that coordinateMode('LAYER') is never 'configured' so will always return here. + if (.not. remap_cs%configured) return + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + + ! Set the bottom depth and negligible thicknesses used in the coordinate remapping in the right units. + if (remap_cs%Z_based_coord) then + h_neglect = set_dz_neglect(GV, US, remap_cs%answer_date, h_neglect_edge) + Z_unit_scale = 1.0 + do j=js-1,je+1 ; do i=is-1,ie+1 + bottom_depth(i,j) = G%bathyT(i,j) + G%Z_ref + enddo ; enddo else - h_neglect = GV%kg_m2_to_H*1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H*1.0e-10 + h_neglect = set_h_neglect(GV, remap_cs%answer_date, h_neglect_edge) + Z_unit_scale = GV%Z_to_H + do j=js-1,je+1 ; do i=is-1,ie+1 + bottom_depth(i,j) = GV%Z_to_H * (G%bathyT(i,j) + G%Z_ref) + enddo ; enddo endif - nz = remap_cs%nz if (.not. remap_cs%initialized) then ! Initialize remapping and regridding on the first call @@ -307,145 +308,203 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state, h_targe remap_cs%initialized = .true. endif + ! Calculate the total thickness of the water column, if it is needed, + if ((remap_cs%vertical_coord == coordinateMode('ZSTAR')) .or. & + (remap_cs%vertical_coord == coordinateMode('SIGMA'))) then + if (remap_CS%answer_date >= 20240201) then + ! Avoid using sum to have a specific order for the vertical sums. + ! For some compilers, the explicit expression gives the same answers as the sum function. + h_tot(:,:) = 0.0 + do k=1,GV%ke ; do j=js-1,je+1 ; do i=is-1,ie+1 + h_tot(i,j) = h_tot(i,j) + h(i,j,k) + enddo ; enddo ; enddo + else + do j=js-1,je+1 ; do i=is-1,ie+1 + h_tot(i,j) = sum(h(i,j,:)) + enddo ; enddo + endif + endif + ! Calculate remapping thicknesses for different target grids based on ! nominal/target interface locations. This happens for every call on the ! assumption that h, T, S has changed. - do j=G%jsc-1, G%jec+1 ; do i=G%isc-1, G%iec+1 - if (G%mask2dT(i,j)==0.) then - h_target(i,j,:) = 0. - cycle - endif + h_target(:,:,:) = 0.0 - if (remap_cs%vertical_coord == coordinateMode('ZSTAR')) then + nz = remap_cs%nz + if (remap_cs%vertical_coord == coordinateMode('ZSTAR')) then + do j=js-1,je+1 ; do i=is-1,ie+1 ; if (G%mask2dT(i,j) > 0.0) then + ! This function call can work with the last 4 arguments all in units of [Z ~> m] or [H ~> kg m-2]. call build_zstar_column(get_zlike_CS(remap_cs%regrid_cs), & - GV%Z_to_H*(G%bathyT(i,j)+G%Z_ref), sum(h(i,j,:)), & - zInterfaces, zScale=GV%Z_to_H) - elseif (remap_cs%vertical_coord == coordinateMode('SIGMA')) then + bottom_depth(i,j), h_tot(i,j), zInterfaces, zScale=Z_unit_scale) + do k=1,nz ; h_target(i,j,k) = zInterfaces(K) - zInterfaces(K+1) ; enddo + endif ; enddo ; enddo + elseif (remap_cs%vertical_coord == coordinateMode('SIGMA')) then + do j=js-1, je+1 ; do i=is-1,ie+1 ; if (G%mask2dT(i,j) > 0.0) then + ! This function call can work with the last 3 arguments all in units of [Z ~> m] or [H ~> kg m-2]. call build_sigma_column(get_sigma_CS(remap_cs%regrid_cs), & - GV%Z_to_H*(G%bathyT(i,j)+G%Z_ref), sum(h(i,j,:)), zInterfaces) - elseif (remap_cs%vertical_coord == coordinateMode('RHO')) then + bottom_depth(i,j), h_tot(i,j), zInterfaces) + do k=1,nz ; h_target(i,j,k) = zInterfaces(K) - zInterfaces(K+1) ; enddo + endif ; enddo ; enddo + elseif (remap_cs%vertical_coord == coordinateMode('RHO')) then + do j=js-1,je+1 ; do i=is-1,ie+1 ; if (G%mask2dT(i,j) > 0.0) then + ! This function call can work with 5 arguments in units of [Z ~> m] or [H ~> kg m-2]. call build_rho_column(get_rho_CS(remap_cs%regrid_cs), GV%ke, & - GV%Z_to_H*(G%bathyT(i,j)+G%Z_ref), h(i,j,:), T(i,j,:), S(i,j,:), & + bottom_depth(i,j), h(i,j,:), T(i,j,:), S(i,j,:), & eqn_of_state, zInterfaces, h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) - elseif (remap_cs%vertical_coord == coordinateMode('HYCOM1')) then -! call build_hycom1_column(remap_cs%regrid_cs, nz, & -! GV%Z_to_H*(G%bathyT(i,j)+G%Z_ref), sum(h(i,j,:)), zInterfaces) - call MOM_error(FATAL,"diag_remap_update: HYCOM1 coordinate not coded for diagnostics yet!") - endif - do k = 1,nz - h_target(i,j,k) = zInterfaces(k) - zInterfaces(k+1) - enddo - enddo ; enddo + do k=1,nz ; h_target(i,j,k) = zInterfaces(K) - zInterfaces(K+1) ; enddo + endif ; enddo ; enddo + elseif (remap_cs%vertical_coord == coordinateMode('HYCOM1')) then + call MOM_error(FATAL,"diag_remap_update: HYCOM1 coordinate not coded for diagnostics yet!") +! do j=js-1,je+1 ; do i=is-1,ie+1 ; if (G%mask2dT(i,j) > 0.0) then +! call build_hycom1_column(remap_cs%regrid_cs, nz, & +! bottom_depth(i,j), h_tot(i,j), zInterfaces) +! do k=1,nz ; h_target(i,j,k) = zInterfaces(K) - zInterfaces(K+1) ; enddo +! endif ; enddo ; enddo + endif end subroutine diag_remap_update !> Remap diagnostic field to alternative vertical grid. -subroutine diag_remap_do_remap(remap_cs, G, GV, h, staggered_in_x, staggered_in_y, & +subroutine diag_remap_do_remap(remap_cs, G, GV, US, h, staggered_in_x, staggered_in_y, & mask, field, remapped_field) - type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coodinate control structure - type(ocean_grid_type), intent(in) :: G !< Ocean grid structure - type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure - real, dimension(:,:,:), intent(in) :: h !< The current thicknesses [H ~> m or kg m-2] - logical, intent(in) :: staggered_in_x !< True is the x-axis location is at u or q points - logical, intent(in) :: staggered_in_y !< True is the y-axis location is at v or q points - real, dimension(:,:,:), pointer :: mask !< A mask for the field [nondim] - real, dimension(:,:,:), intent(in) :: field(:,:,:) !< The diagnostic field to be remapped [A] - real, dimension(:,:,:), intent(inout) :: remapped_field !< Field remapped to new coordinate [A] + type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coordinate control structure + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, dimension(:,:,:), intent(in) :: h !< The current thicknesses [H ~> m or kg m-2] or [Z ~> m], + !! depending on the value of remap_CS%Z_based_coord + logical, intent(in) :: staggered_in_x !< True is the x-axis location is at u or q points + logical, intent(in) :: staggered_in_y !< True is the y-axis location is at v or q points + real, dimension(:,:,:), pointer :: mask !< A mask for the field [nondim]. + real, dimension(:,:,:), intent(in) :: field(:,:,:) !< The diagnostic field to be remapped [A] + real, dimension(:,:,:), intent(out) :: remapped_field !< Field remapped to new coordinate [A] + ! Local variables - real, dimension(remap_cs%nz) :: h_dest ! Destination thicknesses [H ~> m or kg m-2] - real, dimension(size(h,3)) :: h_src ! A column of source thicknesses [H ~> m or kg m-2] - real :: h_neglect, h_neglect_edge ! Negligible thicknesses [H ~> m or kg m-2] - integer :: nz_src, nz_dest - integer :: i, j !< Grid index - integer :: i1, j1 !< 1-based index - integer :: i_lo, i_hi, j_lo, j_hi !< (uv->h) interpolation indices - integer :: shift !< Symmetric offset for 1-based indexing + integer :: isdf, jsdf !< The starting i- and j-indices in memory for field call assert(remap_cs%initialized, 'diag_remap_do_remap: remap_cs not initialized.') call assert(size(field, 3) == size(h, 3), & 'diag_remap_do_remap: Remap field and thickness z-axes do not match.') - if (remap_cs%answer_date >= 20190101) then - h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff - elseif (GV%Boussinesq) then - h_neglect = GV%m_to_H*1.0e-30 ; h_neglect_edge = GV%m_to_H*1.0e-10 + isdf = G%isd ; if (staggered_in_x) Isdf = G%IsdB + jsdf = G%jsd ; if (staggered_in_y) Jsdf = G%JsdB + + if (associated(mask)) then + call do_remap(remap_cs, G, GV, US, isdf, jsdf, h, staggered_in_x, staggered_in_y, & + field, remapped_field, mask(:,:,1)) else - h_neglect = GV%kg_m2_to_H*1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H*1.0e-10 + call do_remap(remap_cs, G, GV, US, isdf, jsdf, h, staggered_in_x, staggered_in_y, & + field, remapped_field) + endif + +end subroutine diag_remap_do_remap + +!> The internal routine to remap a diagnostic field to an alternative vertical grid. +subroutine do_remap(remap_cs, G, GV, US, isdf, jsdf, h, staggered_in_x, staggered_in_y, & + field, remapped_field, mask) + type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coordinate control structure + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + integer, intent(in) :: isdf !< The starting i-index in memory for field + integer, intent(in) :: jsdf !< The starting j-index in memory for field + real, dimension(G%isd:,G%jsd:,:), & + intent(in) :: h !< The current thicknesses [H ~> m or kg m-2] or [Z ~> m], + !! depending on the value of remap_CS%Z_based_coord + logical, intent(in) :: staggered_in_x !< True is the x-axis location is at u or q points + logical, intent(in) :: staggered_in_y !< True is the y-axis location is at v or q points + real, dimension(isdf:,jsdf:,:), & + intent(in) :: field !< The diagnostic field to be remapped [A] + real, dimension(isdf:,jsdf:,:), & + intent(out) :: remapped_field !< Field remapped to new coordinate [A] + real, dimension(isdf:,jsdf:), & + optional, intent(in) :: mask !< A mask for the field [nondim] + + ! Local variables + real, dimension(remap_cs%nz) :: h_dest ! Destination thicknesses [H ~> m or kg m-2] or [Z ~> m] + real, dimension(size(h,3)) :: h_src ! A column of source thicknesses [H ~> m or kg m-2] or [Z ~> m] + real :: h_neglect, h_neglect_edge ! Negligible thicknesses [H ~> m or kg m-2] or [Z ~> m] + integer :: nz_src, nz_dest ! The number of layers on the native and remapped grids + integer :: i, j ! Grid index + + if (remap_cs%Z_based_coord) then + h_neglect = set_dz_neglect(GV, US, remap_cs%answer_date, h_neglect_edge) + else + h_neglect = set_h_neglect(GV, remap_cs%answer_date, h_neglect_edge) endif nz_src = size(field,3) nz_dest = remap_cs%nz remapped_field(:,:,:) = 0. - ! Symmetric grid offset under 1-based indexing; see header for details. - shift = 0 ; if (G%symmetric) shift = 1 - if (staggered_in_x .and. .not. staggered_in_y) then ! U-points - do j=G%jsc, G%jec - do I=G%iscB, G%iecB - I1 = I - G%isdB + 1 - i_lo = I1 - shift; i_hi = i_lo + 1 - if (associated(mask)) then - if (mask(I,j,1) == 0.) cycle - endif - h_src(:) = 0.5 * (h(i_lo,j,:) + h(i_hi,j,:)) - h_dest(:) = 0.5 * (remap_cs%h(i_lo,j,:) + remap_cs%h(i_hi,j,:)) - call remapping_core_h(remap_cs%remap_cs, & - nz_src, h_src(:), field(I1,j,:), & - nz_dest, h_dest(:), remapped_field(I1,j,:), & - h_neglect, h_neglect_edge) - enddo - enddo + if (present(mask)) then + do j=G%jsc,G%jec ; do I=G%IscB,G%IecB ; if (mask(I,j) > 0.) then + h_src(:) = 0.5 * (h(i,j,:) + h(i+1,j,:)) + h_dest(:) = 0.5 * (remap_cs%h(i,j,:) + remap_cs%h(i+1,j,:)) + call remapping_core_h(remap_cs%remap_cs, nz_src, h_src(:), field(I,j,:), & + nz_dest, h_dest(:), remapped_field(I,j,:), h_neglect, h_neglect_edge) + endif ; enddo ; enddo + else + do j=G%jsc,G%jec ; do I=G%IscB,G%IecB + h_src(:) = 0.5 * (h(i,j,:) + h(i+1,j,:)) + h_dest(:) = 0.5 * (remap_cs%h(i,j,:) + remap_cs%h(i+1,j,:)) + call remapping_core_h(remap_cs%remap_cs, nz_src, h_src(:), field(I,j,:), & + nz_dest, h_dest(:), remapped_field(I,j,:), h_neglect, h_neglect_edge) + enddo ; enddo + endif elseif (staggered_in_y .and. .not. staggered_in_x) then ! V-points - do J=G%jscB, G%jecB - J1 = J - G%jsdB + 1 - j_lo = J1 - shift; j_hi = j_lo + 1 - do i=G%isc, G%iec - if (associated(mask)) then - if (mask(i,J,1) == 0.) cycle - endif - h_src(:) = 0.5 * (h(i,j_lo,:) + h(i,j_hi,:)) - h_dest(:) = 0.5 * (remap_cs%h(i,j_lo,:) + remap_cs%h(i,j_hi,:)) - call remapping_core_h(remap_cs%remap_cs, & - nz_src, h_src(:), field(i,J1,:), & - nz_dest, h_dest(:), remapped_field(i,J1,:), & - h_neglect, h_neglect_edge) - enddo - enddo + if (present(mask)) then + do J=G%jscB,G%jecB ; do i=G%isc,G%iec ; if (mask(i,j) > 0.) then + h_src(:) = 0.5 * (h(i,j,:) + h(i,j+1,:)) + h_dest(:) = 0.5 * (remap_cs%h(i,j,:) + remap_cs%h(i,j+1,:)) + call remapping_core_h(remap_cs%remap_cs, nz_src, h_src(:), field(i,J,:), & + nz_dest, h_dest(:), remapped_field(i,J,:), h_neglect, h_neglect_edge) + endif ; enddo ; enddo + else + do J=G%jscB,G%jecB ; do i=G%isc,G%iec + h_src(:) = 0.5 * (h(i,j,:) + h(i,j+1,:)) + h_dest(:) = 0.5 * (remap_cs%h(i,j,:) + remap_cs%h(i,j+1,:)) + call remapping_core_h(remap_cs%remap_cs, nz_src, h_src(:), field(i,J,:), & + nz_dest, h_dest(:), remapped_field(i,J,:), h_neglect, h_neglect_edge) + enddo ; enddo + endif elseif ((.not. staggered_in_x) .and. (.not. staggered_in_y)) then ! H-points - do j=G%jsc, G%jec - do i=G%isc, G%iec - if (associated(mask)) then - if (mask(i,j,1) == 0.) cycle - endif - h_src(:) = h(i,j,:) - h_dest(:) = remap_cs%h(i,j,:) - call remapping_core_h(remap_cs%remap_cs, & - nz_src, h_src(:), field(i,j,:), & - nz_dest, h_dest(:), remapped_field(i,j,:), & + if (present(mask)) then + do j=G%jsc,G%jec ; do i=G%isc,G%iec ; if (mask(i,j) > 0.) then + call remapping_core_h(remap_cs%remap_cs, nz_src, h(i,j,:), field(i,j,:), & + nz_dest, remap_cs%h(i,j,:), remapped_field(i,j,:), & h_neglect, h_neglect_edge) - enddo - enddo + endif ; enddo ; enddo + else + do j=G%jsc,G%jec ; do i=G%isc,G%iec + call remapping_core_h(remap_cs%remap_cs, nz_src, h(i,j,:), field(i,j,:), & + nz_dest, remap_cs%h(i,j,:), remapped_field(i,j,:), & + h_neglect, h_neglect_edge) + enddo ; enddo + endif else call assert(.false., 'diag_remap_do_remap: Unsupported axis combination') endif -end subroutine diag_remap_do_remap +end subroutine do_remap !> Calculate masks for target grid subroutine diag_remap_calc_hmask(remap_cs, G, mask) - type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coodinate control structure + type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coordinate control structure type(ocean_grid_type), intent(in) :: G !< Ocean grid structure - real, dimension(:,:,:), intent(out) :: mask !< h-point mask for target grid [nondim] + real, dimension(G%isd:,G%jsd:,:), & + intent(out) :: mask !< h-point mask for target grid [nondim] + ! Local variables - real, dimension(remap_cs%nz) :: h_dest ! Destination thicknesses [H ~> m or kg m-2] + real, dimension(remap_cs%nz) :: h_dest ! Destination thicknesses [H ~> m or kg m-2] or [Z ~> m] integer :: i, j, k logical :: mask_vanished_layers - real :: h_tot ! Sum of all thicknesses [H ~> m or kg m-2] - real :: h_err ! An estimate of a negligible thickness [H ~> m or kg m-2] + real :: h_tot ! Sum of all thicknesses [H ~> m or kg m-2] or [Z ~> m] + real :: h_err ! An estimate of a negligible thickness [H ~> m or kg m-2] or [Z ~> m] call assert(remap_cs%initialized, 'diag_remap_calc_hmask: remap_cs not initialized.') @@ -453,7 +512,7 @@ subroutine diag_remap_calc_hmask(remap_cs, G, mask) mask_vanished_layers = (remap_cs%vertical_coord == coordinateMode('ZSTAR')) mask(:,:,:) = 0. - do j=G%jsc-1, G%jec+1 ; do i=G%isc-1, G%iec+1 + do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1 if (G%mask2dT(i,j)>0.) then if (mask_vanished_layers) then h_dest(:) = remap_cs%h(i,j,:) @@ -482,166 +541,239 @@ end subroutine diag_remap_calc_hmask !> Vertically re-grid an already vertically-integrated diagnostic field to alternative vertical grid. subroutine vertically_reintegrate_diag_field(remap_cs, G, h, h_target, staggered_in_x, staggered_in_y, & mask, field, reintegrated_field) - type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coodinate control structure - type(ocean_grid_type), intent(in) :: G !< Ocean grid structure - real, dimension(:,:,:), intent(in) :: h !< The thicknesses of the source grid [H ~> m or kg m-2] - real, dimension(:,:,:), intent(in) :: h_target !< The thicknesses of the target grid [H ~> m or kg m-2] - logical, intent(in) :: staggered_in_x !< True is the x-axis location is at u or q points - logical, intent(in) :: staggered_in_y !< True is the y-axis location is at v or q points - real, dimension(:,:,:), pointer :: mask !< A mask for the field [nondim] - real, dimension(:,:,:), intent(in) :: field !< The diagnostic field to be remapped [A] - real, dimension(:,:,:), intent(inout) :: reintegrated_field !< Field argument remapped to alternative coordinate [A] + type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coordinate control structure + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + real, dimension(:,:,:), intent(in) :: h !< The thicknesses of the source grid [H ~> m or kg m-2] or [Z ~> m] + real, dimension(:,:,:), intent(in) :: h_target !< The thicknesses of the target grid [H ~> m or kg m-2] or [Z ~> m] + logical, intent(in) :: staggered_in_x !< True is the x-axis location is at u or q points + logical, intent(in) :: staggered_in_y !< True is the y-axis location is at v or q points + real, dimension(:,:,:), pointer :: mask !< A mask for the field [nondim]. Note that because this + !! is a pointer it retains its declared indexing conventions. + real, dimension(:,:,:), intent(in) :: field !< The diagnostic field to be remapped [A] + real, dimension(:,:,:), intent(out) :: reintegrated_field !< Field argument remapped to alternative coordinate [A] + ! Local variables - real, dimension(remap_cs%nz) :: h_dest ! Destination thicknesses [H ~> m or kg m-2] - real, dimension(size(h,3)) :: h_src ! A column of source thicknesses [H ~> m or kg m-2] - integer :: nz_src, nz_dest - integer :: i, j !< Grid index - integer :: i1, j1 !< 1-based index - integer :: i_lo, i_hi, j_lo, j_hi !< (uv->h) interpolation indices - integer :: shift !< Symmetric offset for 1-based indexing + integer :: isdf, jsdf !< The starting i- and j-indices in memory for field call assert(remap_cs%initialized, 'vertically_reintegrate_diag_field: remap_cs not initialized.') call assert(size(field, 3) == size(h, 3), & 'vertically_reintegrate_diag_field: Remap field and thickness z-axes do not match.') + isdf = G%isd ; if (staggered_in_x) Isdf = G%IsdB + jsdf = G%jsd ; if (staggered_in_y) Jsdf = G%JsdB + + if (associated(mask)) then + call vertically_reintegrate_field(remap_cs, G, isdf, jsdf, h, h_target, staggered_in_x, staggered_in_y, & + field, reintegrated_field, mask(:,:,1)) + else + call vertically_reintegrate_field(remap_cs, G, isdf, jsdf, h, h_target, staggered_in_x, staggered_in_y, & + field, reintegrated_field) + endif + +end subroutine vertically_reintegrate_diag_field + +!> The internal routine to vertically re-grid an already vertically-integrated diagnostic field to +!! an alternative vertical grid. +subroutine vertically_reintegrate_field(remap_cs, G, isdf, jsdf, h, h_target, staggered_in_x, staggered_in_y, & + field, reintegrated_field, mask) + type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coordinate control structure + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + integer, intent(in) :: isdf !< The starting i-index in memory for field + integer, intent(in) :: jsdf !< The starting j-index in memory for field + real, dimension(G%isd:,G%jsd:,:), & + intent(in) :: h !< The thicknesses of the source grid [H ~> m or kg m-2] or [Z ~> m] + real, dimension(G%isd:,G%jsd:,:), & + intent(in) :: h_target !< The thicknesses of the target grid [H ~> m or kg m-2] or [Z ~> m] + logical, intent(in) :: staggered_in_x !< True is the x-axis location is at u or q points + logical, intent(in) :: staggered_in_y !< True is the y-axis location is at v or q points + real, dimension(isdf:,jsdf:,:), & + intent(in) :: field !< The diagnostic field to be remapped [A] + real, dimension(isdf:,jsdf:,:), & + intent(out) :: reintegrated_field !< Field argument remapped to alternative coordinate [A] + real, dimension(isdf:,jsdf:), & + optional, intent(in) :: mask !< A mask for the field [nondim] + + ! Local variables + real, dimension(remap_cs%nz) :: h_dest ! Destination thicknesses [H ~> m or kg m-2] or [Z ~> m] + real, dimension(size(h,3)) :: h_src ! A column of source thicknesses [H ~> m or kg m-2] or [Z ~> m] + integer :: nz_src, nz_dest ! The number of layers on the native and remapped grids + integer :: i, j ! Grid index + nz_src = size(field,3) nz_dest = remap_cs%nz reintegrated_field(:,:,:) = 0. - ! Symmetric grid offset under 1-based indexing; see header for details. - shift = 0 ; if (G%symmetric) shift = 1 - if (staggered_in_x .and. .not. staggered_in_y) then ! U-points - do j=G%jsc, G%jec - do I=G%iscB, G%iecB - I1 = I - G%isdB + 1 - i_lo = I1 - shift; i_hi = i_lo + 1 - if (associated(mask)) then - if (mask(I,j,1) == 0.) cycle - endif - h_src(:) = 0.5 * (h(i_lo,j,:) + h(i_hi,j,:)) - h_dest(:) = 0.5 * (h_target(i_lo,j,:) + h_target(i_hi,j,:)) - call reintegrate_column(nz_src, h_src, field(I1,j,:), & - nz_dest, h_dest, reintegrated_field(I1,j,:)) - enddo - enddo + if (present(mask)) then + do j=G%jsc,G%jec ; do I=G%IscB,G%IecB ; if (mask(I,j) > 0.0) then + h_src(:) = 0.5 * (h(i,j,:) + h(i+1,j,:)) + h_dest(:) = 0.5 * (h_target(i,j,:) + h_target(i+1,j,:)) + call reintegrate_column(nz_src, h_src, field(I,j,:), & + nz_dest, h_dest, reintegrated_field(I,j,:)) + endif ; enddo ; enddo + else + do j=G%jsc,G%jec ; do I=G%IscB,G%IecB + h_src(:) = 0.5 * (h(i,j,:) + h(i+1,j,:)) + h_dest(:) = 0.5 * (h_target(i,j,:) + h_target(i+1,j,:)) + call reintegrate_column(nz_src, h_src, field(I,j,:), & + nz_dest, h_dest, reintegrated_field(I,j,:)) + enddo ; enddo + endif elseif (staggered_in_y .and. .not. staggered_in_x) then ! V-points - do J=G%jscB, G%jecB - J1 = J - G%jsdB + 1 - j_lo = J1 - shift; j_hi = j_lo + 1 - do i=G%isc, G%iec - if (associated(mask)) then - if (mask(i,J,1) == 0.) cycle - endif - h_src(:) = 0.5 * (h(i,j_lo,:) + h(i,j_hi,:)) - h_dest(:) = 0.5 * (h_target(i,j_lo,:) + h_target(i,j_hi,:)) - call reintegrate_column(nz_src, h_src, field(i,J1,:), & - nz_dest, h_dest, reintegrated_field(i,J1,:)) - enddo - enddo + if (present(mask)) then + do J=G%jscB,G%jecB ; do i=G%isc,G%iec ; if (mask(i,J) > 0.0) then + h_src(:) = 0.5 * (h(i,j,:) + h(i,j+1,:)) + h_dest(:) = 0.5 * (h_target(i,j,:) + h_target(i,j+1,:)) + call reintegrate_column(nz_src, h_src, field(i,J,:), & + nz_dest, h_dest, reintegrated_field(i,J,:)) + endif ; enddo ; enddo + else + do J=G%jscB,G%jecB ; do i=G%isc,G%iec + h_src(:) = 0.5 * (h(i,j,:) + h(i,j+1,:)) + h_dest(:) = 0.5 * (h_target(i,j,:) + h_target(i,j+1,:)) + call reintegrate_column(nz_src, h_src, field(i,J,:), & + nz_dest, h_dest, reintegrated_field(i,J,:)) + enddo ; enddo + endif elseif ((.not. staggered_in_x) .and. (.not. staggered_in_y)) then ! H-points - do j=G%jsc, G%jec - do i=G%isc, G%iec - if (associated(mask)) then - if (mask(i,j,1) == 0.) cycle - endif - h_src(:) = h(i,j,:) - h_dest(:) = h_target(i,j,:) - call reintegrate_column(nz_src, h_src, field(i,j,:), & - nz_dest, h_dest, reintegrated_field(i,j,:)) - enddo - enddo + if (present(mask)) then + do j=G%jsc,G%jec ; do i=G%isc,G%iec ; if (mask(i,J) > 0.0) then + call reintegrate_column(nz_src, h(i,j,:), field(i,j,:), & + nz_dest, h_target(i,j,:), reintegrated_field(i,j,:)) + endif ; enddo ; enddo + else + do j=G%jsc,G%jec ; do i=G%isc,G%iec + call reintegrate_column(nz_src, h(i,j,:), field(i,j,:), & + nz_dest, h_target(i,j,:), reintegrated_field(i,j,:)) + enddo ; enddo + endif else call assert(.false., 'vertically_reintegrate_diag_field: Q point remapping is not coded yet.') endif -end subroutine vertically_reintegrate_diag_field +end subroutine vertically_reintegrate_field !> Vertically interpolate diagnostic field to alternative vertical grid. subroutine vertically_interpolate_diag_field(remap_cs, G, h, staggered_in_x, staggered_in_y, & mask, field, interpolated_field) - type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coodinate control structure + type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coordinate control structure type(ocean_grid_type), intent(in) :: G !< Ocean grid structure - real, dimension(:,:,:), intent(in) :: h !< The current thicknesses [H ~> m or kg m-2] + real, dimension(:,:,:), intent(in) :: h !< The current thicknesses [H ~> m or kg m-2] or [Z ~> m], + !! depending on the value of remap_cs%Z_based_coord logical, intent(in) :: staggered_in_x !< True is the x-axis location is at u or q points logical, intent(in) :: staggered_in_y !< True is the y-axis location is at v or q points - real, dimension(:,:,:), pointer :: mask !< A mask for the field [nondim] + real, dimension(:,:,:), pointer :: mask !< A mask for the field [nondim]. Note that because this + !! is a pointer it retains its declared indexing conventions. real, dimension(:,:,:), intent(in) :: field !< The diagnostic field to be remapped [A] real, dimension(:,:,:), intent(inout) :: interpolated_field !< Field argument remapped to alternative coordinate [A] + ! Local variables - real, dimension(remap_cs%nz) :: h_dest ! Destination thicknesses [H ~> m or kg m-2] - real, dimension(size(h,3)) :: h_src ! A column of source thicknesses [H ~> m or kg m-2] - integer :: nz_src, nz_dest - integer :: i, j !< Grid index - integer :: i1, j1 !< 1-based index - integer :: i_lo, i_hi, j_lo, j_hi !< (uv->h) interpolation indices - integer :: shift !< Symmetric offset for 1-based indexing + integer :: isdf, jsdf !< The starting i- and j-indices in memory for field call assert(remap_cs%initialized, 'vertically_interpolate_diag_field: remap_cs not initialized.') call assert(size(field, 3) == size(h, 3)+1, & 'vertically_interpolate_diag_field: Remap field and thickness z-axes do not match.') + isdf = G%isd ; if (staggered_in_x) Isdf = G%IsdB + jsdf = G%jsd ; if (staggered_in_y) Jsdf = G%JsdB + + if (associated(mask)) then + call vertically_interpolate_field(remap_cs, G, isdf, jsdf, h, staggered_in_x, staggered_in_y, & + field, interpolated_field, mask(:,:,1)) + else + call vertically_interpolate_field(remap_cs, G, isdf, jsdf, h, staggered_in_x, staggered_in_y, & + field, interpolated_field) + endif + +end subroutine vertically_interpolate_diag_field + +!> Internal routine to vertically interpolate a diagnostic field to an alternative vertical grid. +subroutine vertically_interpolate_field(remap_cs, G, isdf, jsdf, h, staggered_in_x, staggered_in_y, & + field, interpolated_field, mask) + type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coordinate control structure + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + integer, intent(in) :: isdf !< The starting i-index in memory for field + integer, intent(in) :: jsdf !< The starting j-index in memory for field + real, dimension(G%isd:,G%jsd:,:), & + intent(in) :: h !< The current thicknesses [H ~> m or kg m-2] or [Z ~> m], + !! depending on the value of remap_cs%Z_based_coord + logical, intent(in) :: staggered_in_x !< True is the x-axis location is at u or q points + logical, intent(in) :: staggered_in_y !< True is the y-axis location is at v or q points + real, dimension(isdf:,jsdf:,:), & + intent(in) :: field !< The diagnostic field to be remapped [A] + real, dimension(isdf:,jsdf:,:), & + intent(out) :: interpolated_field !< Field argument remapped to alternative coordinate [A] + real, dimension(isdf:,jsdf:), & + optional, intent(in) :: mask !< A mask for the field [nondim] + + ! Local variables + real, dimension(remap_cs%nz) :: h_dest ! Destination thicknesses [H ~> m or kg m-2] or [Z ~> m] + real, dimension(size(h,3)) :: h_src ! A column of source thicknesses [H ~> m or kg m-2] or [Z ~> m] + integer :: nz_src, nz_dest ! The number of layers on the native and remapped grids + integer :: i, j !< Grid index + interpolated_field(:,:,:) = 0. nz_src = size(h,3) nz_dest = remap_cs%nz - ! Symmetric grid offset under 1-based indexing; see header for details. - shift = 0 ; if (G%symmetric) shift = 1 - if (staggered_in_x .and. .not. staggered_in_y) then ! U-points - do j=G%jsc, G%jec - do I=G%iscB, G%iecB - I1 = I - G%isdB + 1 - i_lo = I1 - shift; i_hi = i_lo + 1 - if (associated(mask)) then - if (mask(I,j,1) == 0.) cycle - endif - h_src(:) = 0.5 * (h(i_lo,j,:) + h(i_hi,j,:)) - h_dest(:) = 0.5 * (remap_cs%h(i_lo,j,:) + remap_cs%h(i_hi,j,:)) - call interpolate_column(nz_src, h_src, field(I1,j,:), & - nz_dest, h_dest, interpolated_field(I1,j,:), .true.) - enddo - enddo + if (present(mask)) then + do j=G%jsc,G%jec ; do I=G%IscB,G%IecB ; if (mask(I,j) > 0.0) then + h_src(:) = 0.5 * (h(i,j,:) + h(i+1,j,:)) + h_dest(:) = 0.5 * (remap_cs%h(i,j,:) + remap_cs%h(i+1,j,:)) + call interpolate_column(nz_src, h_src, field(I,j,:), & + nz_dest, h_dest, interpolated_field(I,j,:), .true.) + endif ; enddo ; enddo + else + do j=G%jsc,G%jec ; do I=G%IscB,G%IecB + h_src(:) = 0.5 * (h(i,j,:) + h(i+1,j,:)) + h_dest(:) = 0.5 * (remap_cs%h(i,j,:) + remap_cs%h(i+1,j,:)) + call interpolate_column(nz_src, h_src, field(I,j,:), & + nz_dest, h_dest, interpolated_field(I,j,:), .true.) + enddo ; enddo + endif elseif (staggered_in_y .and. .not. staggered_in_x) then ! V-points - do J=G%jscB, G%jecB - J1 = J - G%jsdB + 1 - j_lo = J1 - shift; j_hi = j_lo + 1 - do i=G%isc, G%iec - if (associated(mask)) then - if (mask(i,J,1) == 0.) cycle - endif - h_src(:) = 0.5 * (h(i,j_lo,:) + h(i,j_hi,:)) - h_dest(:) = 0.5 * (remap_cs%h(i,j_lo,:) + remap_cs%h(i,j_hi,:)) - call interpolate_column(nz_src, h_src, field(i,J1,:), & - nz_dest, h_dest, interpolated_field(i,J1,:), .true.) - enddo - enddo + if (present(mask)) then + do J=G%jscB,G%jecB ; do i=G%isc,G%iec ; if (mask(I,j) > 0.0) then + h_src(:) = 0.5 * (h(i,j,:) + h(i,j+1,:)) + h_dest(:) = 0.5 * (remap_cs%h(i,j,:) + remap_cs%h(i,j+1,:)) + call interpolate_column(nz_src, h_src, field(i,J,:), & + nz_dest, h_dest, interpolated_field(i,J,:), .true.) + endif ; enddo ; enddo + else + do J=G%jscB,G%jecB ; do i=G%isc,G%iec + h_src(:) = 0.5 * (h(i,j,:) + h(i,j+1,:)) + h_dest(:) = 0.5 * (remap_cs%h(i,j,:) + remap_cs%h(i,j+1,:)) + call interpolate_column(nz_src, h_src, field(i,J,:), & + nz_dest, h_dest, interpolated_field(i,J,:), .true.) + enddo ; enddo + endif elseif ((.not. staggered_in_x) .and. (.not. staggered_in_y)) then ! H-points - do j=G%jsc, G%jec - do i=G%isc, G%iec - if (associated(mask)) then - if (mask(i,j,1) == 0.) cycle - endif - h_src(:) = h(i,j,:) - h_dest(:) = remap_cs%h(i,j,:) - call interpolate_column(nz_src, h_src, field(i,j,:), & - nz_dest, h_dest, interpolated_field(i,j,:), .true.) - enddo - enddo + if (present(mask)) then + do j=G%jsc,G%jec ; do i=G%isc,G%iec ; if (mask(i,j) > 0.0) then + call interpolate_column(nz_src, h(i,j,:), field(i,j,:), & + nz_dest, remap_cs%h(i,j,:), interpolated_field(i,j,:), .true.) + endif ; enddo ; enddo + else + do j=G%jsc,G%jec ; do i=G%isc,G%iec + call interpolate_column(nz_src, h(i,j,:), field(i,j,:), & + nz_dest, remap_cs%h(i,j,:), interpolated_field(i,j,:), .true.) + enddo ; enddo + endif else call assert(.false., 'vertically_interpolate_diag_field: Q point remapping is not coded yet.') endif -end subroutine vertically_interpolate_diag_field +end subroutine vertically_interpolate_field -!> Horizontally average field +!> Horizontally average a diagnostic field subroutine horizontally_average_diag_field(G, GV, h, staggered_in_x, staggered_in_y, & is_layer, is_extensive, & field, averaged_field, & @@ -654,8 +786,37 @@ subroutine horizontally_average_diag_field(G, GV, h, staggered_in_x, staggered_i logical, intent(in) :: is_layer !< True if the z-axis location is at h points logical, intent(in) :: is_extensive !< True if the z-direction is spatially integrated (over layers) real, dimension(:,:,:), intent(in) :: field !< The diagnostic field to be remapped [A] - real, dimension(:), intent(inout) :: averaged_field !< Field argument horizontally averaged [A] - logical, dimension(:), intent(inout) :: averaged_mask !< Mask for horizontally averaged field [nondim] + real, dimension(:), intent(out) :: averaged_field !< Field argument horizontally averaged [A] + logical, dimension(:), intent(out) :: averaged_mask !< Mask for horizontally averaged field [nondim] + + ! Local variables + integer :: isdf, jsdf !< The starting i- and j-indices in memory for field + + isdf = G%isd ; if (staggered_in_x) Isdf = G%IsdB + jsdf = G%jsd ; if (staggered_in_y) Jsdf = G%JsdB + + call horizontally_average_field(G, GV, isdf, jsdf, h, staggered_in_x, staggered_in_y, & + is_layer, is_extensive, field, averaged_field, averaged_mask) + +end subroutine horizontally_average_diag_field + +!> Horizontally average a diagnostic field +subroutine horizontally_average_field(G, GV, isdf, jsdf, h, staggered_in_x, staggered_in_y, & + is_layer, is_extensive, field, averaged_field, averaged_mask) + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean vertical grid structure + integer, intent(in) :: isdf !< The starting i-index in memory for field + integer, intent(in) :: jsdf !< The starting j-index in memory for field + real, dimension(G%isd:,G%jsd:,:), & + intent(in) :: h !< The current thicknesses [H ~> m or kg m-2] + logical, intent(in) :: staggered_in_x !< True if the x-axis location is at u or q points + logical, intent(in) :: staggered_in_y !< True if the y-axis location is at v or q points + logical, intent(in) :: is_layer !< True if the z-axis location is at h points + logical, intent(in) :: is_extensive !< True if the z-direction is spatially integrated (over layers) + real, dimension(isdf:,jsdf:,:), & + intent(in) :: field !< The diagnostic field to be remapped [A] + real, dimension(:), intent(out) :: averaged_field !< Field argument horizontally averaged [A] + logical, dimension(:), intent(out) :: averaged_mask !< Mask for horizontally averaged field [nondim] ! Local variables real :: volume(G%isc:G%iec, G%jsc:G%jec, size(field,3)) ! The area [m2], volume [m3] or mass [kg] of each cell. @@ -670,7 +831,6 @@ subroutine horizontally_average_diag_field(G, GV, h, staggered_in_x, staggered_i type(EFP_type), dimension(2*size(field,3)) :: sums_EFP ! Sums of volume or stuff by layer real :: height ! An average thickness attributed to an velocity point [H ~> m or kg m-2] integer :: i, j, k, nz - integer :: i1, j1 !< 1-based index nz = size(field, 3) @@ -686,26 +846,23 @@ subroutine horizontally_average_diag_field(G, GV, h, staggered_in_x, staggered_i stuff_sum(k) = 0. if (is_extensive) then do j=G%jsc, G%jec ; do I=G%isc, G%iec - I1 = I - G%isdB + 1 volume(I,j,k) = (G%US%L_to_m**2 * G%areaCu(I,j)) * G%mask2dCu(I,j) - stuff(I,j,k) = volume(I,j,k) * field(I1,j,k) + stuff(I,j,k) = volume(I,j,k) * field(I,j,k) enddo ; enddo else ! Intensive do j=G%jsc, G%jec ; do I=G%isc, G%iec - I1 = i - G%isdB + 1 height = 0.5 * (h(i,j,k) + h(i+1,j,k)) volume(I,j,k) = (G%US%L_to_m**2 * G%areaCu(I,j)) & * (GV%H_to_MKS * height) * G%mask2dCu(I,j) - stuff(I,j,k) = volume(I,j,k) * field(I1,j,k) + stuff(I,j,k) = volume(I,j,k) * field(I,j,k) enddo ; enddo endif enddo else ! Interface do k=1,nz do j=G%jsc, G%jec ; do I=G%isc, G%iec - I1 = I - G%isdB + 1 volume(I,j,k) = (G%US%L_to_m**2 * G%areaCu(I,j)) * G%mask2dCu(I,j) - stuff(I,j,k) = volume(I,j,k) * field(I1,j,k) + stuff(I,j,k) = volume(I,j,k) * field(I,j,k) enddo ; enddo enddo endif @@ -715,26 +872,23 @@ subroutine horizontally_average_diag_field(G, GV, h, staggered_in_x, staggered_i do k=1,nz if (is_extensive) then do J=G%jsc, G%jec ; do i=G%isc, G%iec - J1 = J - G%jsdB + 1 volume(i,J,k) = (G%US%L_to_m**2 * G%areaCv(i,J)) * G%mask2dCv(i,J) - stuff(i,J,k) = volume(i,J,k) * field(i,J1,k) + stuff(i,J,k) = volume(i,J,k) * field(i,J,k) enddo ; enddo else ! Intensive do J=G%jsc, G%jec ; do i=G%isc, G%iec - J1 = J - G%jsdB + 1 height = 0.5 * (h(i,j,k) + h(i,j+1,k)) volume(i,J,k) = (G%US%L_to_m**2 * G%areaCv(i,J)) & * (GV%H_to_MKS * height) * G%mask2dCv(i,J) - stuff(i,J,k) = volume(i,J,k) * field(i,J1,k) + stuff(i,J,k) = volume(i,J,k) * field(i,J,k) enddo ; enddo endif enddo else ! Interface do k=1,nz do J=G%jsc, G%jec ; do i=G%isc, G%iec - J1 = J - G%jsdB + 1 volume(i,J,k) = (G%US%L_to_m**2 * G%areaCv(i,J)) * G%mask2dCv(i,J) - stuff(i,J,k) = volume(i,J,k) * field(i,J1,k) + stuff(i,J,k) = volume(i,J,k) * field(i,J,k) enddo ; enddo enddo endif @@ -794,6 +948,6 @@ subroutine horizontally_average_diag_field(G, GV, h, staggered_in_x, staggered_i endif enddo -end subroutine horizontally_average_diag_field +end subroutine horizontally_average_field end module MOM_diag_remap From d311b1da1b4d7ff562bac80272604366de270ab5 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 5 Feb 2024 09:56:35 -0500 Subject: [PATCH 14/15] +(*)Use z-space remapping for some diagnostics Do remapping in Z-space for some remapped diagnostics, depending on which coordinate is used. The subroutine thickness_to_dz is used with the thermo_vars_type to do the rescaling properly in non-Boussinesq mode. A new thermo_var_ptrs argument was added to diag_mediator_init, replacing three other arguments that had the same information, and its call from MOM.F90 was modified accordingly. The various calls to the diag_remap routines from post_data_3d and diag_update_remap_grids were modified depending on whether a z-unit or h-unit vertical grid is being remapped to. All answers and diagnostics are identical in Boussinesq mode, but some remapped diagnostics are changed (by not using expressions that depend on the Boussinesq reference density) in the non-Boussinesq mode. There are altered or augmented public arguments to two publicly visible routines. --- src/core/MOM.F90 | 5 +- src/framework/MOM_diag_mediator.F90 | 158 +++++++++++++++++++++------- src/framework/MOM_diag_remap.F90 | 13 +-- 3 files changed, 131 insertions(+), 45 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 2bf2b11b1a..094e534312 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -3152,14 +3152,15 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & ! Set up pointers within diag mediator control structure, ! this needs to occur _after_ CS%h etc. have been allocated. - call diag_set_state_ptrs(CS%h, CS%T, CS%S, CS%tv%eqn_of_state, diag) + call diag_set_state_ptrs(CS%h, CS%tv, diag) ! This call sets up the diagnostic axes. These are needed, ! e.g. to generate the target grids below. call set_axes_info(G, GV, US, param_file, diag) ! Whenever thickness/T/S changes let the diag manager know, target grids - ! for vertical remapping may need to be regenerated. + ! for vertical remapping may need to be regenerated. In non-Boussinesq mode, + ! calc_derived_thermo needs to be called before diag_update_remap_grids. call diag_update_remap_grids(diag) ! Setup the diagnostic grid storage types diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 61bfa93046..6c90e26b98 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -21,15 +21,18 @@ module MOM_diag_mediator use MOM_diag_remap, only : diag_remap_configure_axes, diag_remap_axes_configured use MOM_diag_remap, only : diag_remap_diag_registration_closed, diag_remap_set_active use MOM_EOS, only : EOS_type -use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe, assert +use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe, assert, callTree_showQuery +use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_grid, only : ocean_grid_type +use MOM_interface_heights, only : thickness_to_dz use MOM_io, only : slasher, vardesc, query_vardesc, MOM_read_data use MOM_io, only : get_filename_appendix use MOM_safe_alloc, only : safe_alloc_ptr, safe_alloc_alloc use MOM_string_functions, only : lowercase use MOM_time_manager, only : time_type use MOM_unit_scaling, only : unit_scale_type +use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type implicit none ; private @@ -134,7 +137,7 @@ module MOM_diag_mediator !> Contains an array to store a diagnostic target grid type, private :: diag_grids_type - real, dimension(:,:,:), allocatable :: h !< Target grid for remapped coordinate + real, dimension(:,:,:), allocatable :: h !< Target grid for remapped coordinate [H ~> m or kg m-2] or [Z ~> m] end type diag_grids_type !> Stores all the remapping grids and the model's native space thicknesses @@ -238,6 +241,7 @@ module MOM_diag_mediator integer :: chksum_iounit = -1 !< The unit number of a diagnostic documentation file. !! This file is open if available_diag_doc_unit is > 0. logical :: diag_as_chksum !< If true, log chksums in a text file instead of posting diagnostics + logical :: show_call_tree !< Display the call tree while running. Set by VERBOSITY level. logical :: grid_space_axes !< If true, diagnostic horizontal coordinates axes are in grid space. ! The following fields are used for the output of the data. integer :: is !< The start i-index of cell centers within the computational domain @@ -311,7 +315,9 @@ module MOM_diag_mediator real, dimension(:,:,:), pointer :: h => null() !< The thicknesses needed for remapping [H ~> m or kg m-2] real, dimension(:,:,:), pointer :: T => null() !< The temperatures needed for remapping [C ~> degC] real, dimension(:,:,:), pointer :: S => null() !< The salinities needed for remapping [S ~> ppt] - type(EOS_type), pointer :: eqn_of_state => null() !< The equation of state type + type(EOS_type), pointer :: eqn_of_state => null() !< The equation of state type + type(thermo_var_ptrs), pointer :: tv => null() !< A sturcture with thermodynamic variables that are + !! are used to convert thicknesses to vertical extents type(ocean_grid_type), pointer :: G => null() !< The ocean grid type type(verticalGrid_type), pointer :: GV => null() !< The model's vertical ocean grid type(unit_scale_type), pointer :: US => null() !< A dimensional unit scaling type @@ -329,7 +335,7 @@ module MOM_diag_mediator integer :: num_chksum_diags real, dimension(:,:,:), allocatable :: h_begin !< Layer thicknesses at the beginning of the timestep used - !! for remapping of extensive variables + !! for remapping of extensive variables [H ~> m or kg m-2] end type diag_ctrl @@ -1525,13 +1531,15 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask, alt_h) ! Local variables type(diag_type), pointer :: diag => null() real, dimension(:,:,:), allocatable :: remapped_field - logical :: staggered_in_x, staggered_in_y + logical :: staggered_in_x, staggered_in_y, dz_diag_needed, dz_begin_needed real, dimension(:,:,:), pointer :: h_diag => NULL() + real, dimension(diag_cs%G%isd:diag_cS%G%ied, diag_cs%G%jsd:diag_cS%G%jed, diag_cs%GV%ke) :: & + dz_diag, & ! Layer vertical extents for remapping [Z ~> m] + dz_begin ! Layer vertical extents for remapping extensive quantities [Z ~> m] if (id_clock_diag_mediator>0) call cpu_clock_begin(id_clock_diag_mediator) - ! For intensive variables only, we can choose to use a different diagnostic grid - ! to map to + ! For intensive variables only, we can choose to use a different diagnostic grid to map to if (present(alt_h)) then h_diag => alt_h else @@ -1542,6 +1550,32 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask, alt_h) ! grids, and post each. call assert(diag_field_id < diag_cs%next_free_diag_id, & 'post_data_3d: Unregistered diagnostic id') + + if (diag_cs%show_call_tree) & + call callTree_enter("post_data_3d("//trim(diag_cs%diags(diag_field_id)%debug_str)//")") + + ! Find out whether there are any z-based diagnostics + diag => diag_cs%diags(diag_field_id) + dz_diag_needed = .false. ; dz_begin_needed = .false. + do while (associated(diag)) + if (diag%v_extensive .and. .not.diag%axes%is_native) then + if (diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number)%Z_based_coord) & + dz_begin_needed = .true. + elseif (diag%axes%needs_remapping .or. diag%axes%needs_interpolating) then + if (diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number)%Z_based_coord) & + dz_diag_needed = .true. + endif + diag => diag%next + enddo + + ! Determine the diagnostic grid spacing in height units, if it is needed. + if (dz_diag_needed) then + call thickness_to_dz(h_diag, diag_cs%tv, dz_diag, diag_cs%G, diag_cs%GV, diag_cs%US, halo_size=1) + endif + if (dz_begin_needed) then + call thickness_to_dz(diag_cs%h_begin, diag_cs%tv, dz_begin, diag_cs%G, diag_cs%GV, diag_cs%US, halo_size=1) + endif + diag => diag_cs%diags(diag_field_id) do while (associated(diag)) call assert(associated(diag%axes), 'post_data_3d: axes is not associated') @@ -1557,11 +1591,17 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask, alt_h) if (id_clock_diag_remap>0) call cpu_clock_begin(id_clock_diag_remap) allocate(remapped_field(size(field,1), size(field,2), diag%axes%nz)) - call vertically_reintegrate_diag_field( & - diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number), diag_cs%G, & - diag_cs%h_begin, & - diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number)%h_extensive, & - staggered_in_x, staggered_in_y, diag%axes%mask3d, field, remapped_field) + if (diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number)%Z_based_coord) then + call vertically_reintegrate_diag_field( & + diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number), diag_cs%G, & + dz_begin, diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number)%h_extensive, & + staggered_in_x, staggered_in_y, diag%axes%mask3d, field, remapped_field) + else + call vertically_reintegrate_diag_field( & + diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number), diag_cs%G, & + diag_cs%h_begin, diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number)%h_extensive, & + staggered_in_x, staggered_in_y, diag%axes%mask3d, field, remapped_field) + endif if (id_clock_diag_remap>0) call cpu_clock_end(id_clock_diag_remap) if (associated(diag%axes%mask3d)) then ! Since 3d masks do not vary in the vertical, just use as much as is @@ -1582,9 +1622,15 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask, alt_h) if (id_clock_diag_remap>0) call cpu_clock_begin(id_clock_diag_remap) allocate(remapped_field(size(field,1), size(field,2), diag%axes%nz)) - call diag_remap_do_remap(diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number), & - diag_cs%G, diag_cs%GV, diag_cs%US, h_diag, staggered_in_x, staggered_in_y, & - diag%axes%mask3d, field, remapped_field) + if (diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number)%Z_based_coord) then + call diag_remap_do_remap(diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number), & + diag_cs%G, diag_cs%GV, diag_cs%US, dz_diag, staggered_in_x, staggered_in_y, & + diag%axes%mask3d, field, remapped_field) + else + call diag_remap_do_remap(diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number), & + diag_cs%G, diag_cs%GV, diag_cs%US, h_diag, staggered_in_x, staggered_in_y, & + diag%axes%mask3d, field, remapped_field) + endif if (id_clock_diag_remap>0) call cpu_clock_end(id_clock_diag_remap) if (associated(diag%axes%mask3d)) then ! Since 3d masks do not vary in the vertical, just use as much as is @@ -1605,10 +1651,15 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask, alt_h) if (id_clock_diag_remap>0) call cpu_clock_begin(id_clock_diag_remap) allocate(remapped_field(size(field,1), size(field,2), diag%axes%nz+1)) - call vertically_interpolate_diag_field(diag_cs%diag_remap_cs( & - diag%axes%vertical_coordinate_number), & - diag_cs%G, h_diag, staggered_in_x, staggered_in_y, & - diag%axes%mask3d, field, remapped_field) + if (diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number)%Z_based_coord) then + call vertically_interpolate_diag_field(diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number), & + diag_cs%G, dz_diag, staggered_in_x, staggered_in_y, & + diag%axes%mask3d, field, remapped_field) + else + call vertically_interpolate_diag_field(diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number), & + diag_cs%G, h_diag, staggered_in_x, staggered_in_y, & + diag%axes%mask3d, field, remapped_field) + endif if (id_clock_diag_remap>0) call cpu_clock_end(id_clock_diag_remap) if (associated(diag%axes%mask3d)) then ! Since 3d masks do not vary in the vertical, just use as much as is @@ -1628,6 +1679,9 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask, alt_h) enddo if (id_clock_diag_mediator>0) call cpu_clock_end(id_clock_diag_mediator) + if (diag_cs%show_call_tree) & + call callTree_leave("post_data_3d("//trim(diag_cs%diags(diag_field_id)%debug_str)//")") + end subroutine post_data_3d !> Make a real 3-d array diagnostic available for averaging or output @@ -1926,8 +1980,7 @@ subroutine post_xy_average(diag_cs, diag, field) call horizontally_average_diag_field(diag_cs%G, diag_cs%GV, diag_cs%h, & staggered_in_x, staggered_in_y, & diag%axes%is_layer, diag%v_extensive, & - field, & - averaged_field, averaged_mask) + field, averaged_field, averaged_mask) else nz = size(field, 3) coord = diag%axes%vertical_coordinate_number @@ -3168,6 +3221,8 @@ subroutine diag_mediator_init(G, GV, US, nz, param_file, diag_cs, doc_file_dir) call initialize_diag_type(diag_cs%diags(i)) enddo + diag_cs%show_call_tree = callTree_showQuery() + ! Read all relevant parameters and write them to the model log. call log_version(param_file, mdl, version, "") @@ -3223,7 +3278,7 @@ subroutine diag_mediator_init(G, GV, US, nz, param_file, diag_cs, doc_file_dir) if (diag_cs%diag_as_chksum) & diag_cs%num_chksum_diags = 0 - ! Keep pointers grid, h, T, S needed diagnostic remapping + ! Keep pointers to the grid, h, T, S needed for diagnostic remapping diag_cs%G => G diag_cs%GV => GV diag_cs%US => US @@ -3231,6 +3286,7 @@ subroutine diag_mediator_init(G, GV, US, nz, param_file, diag_cs, doc_file_dir) diag_cs%T => null() diag_cs%S => null() diag_cs%eqn_of_state => null() + diag_cs%tv => null() allocate(diag_cs%h_begin(G%isd:G%ied,G%jsd:G%jed,nz)) #if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__) @@ -3343,18 +3399,18 @@ subroutine diag_mediator_init(G, GV, US, nz, param_file, diag_cs, doc_file_dir) end subroutine diag_mediator_init !> Set pointers to the default state fields used to remap diagnostics. -subroutine diag_set_state_ptrs(h, T, S, eqn_of_state, diag_cs) +subroutine diag_set_state_ptrs(h, tv, diag_cs) real, dimension(:,:,:), target, intent(in ) :: h !< the model thickness array [H ~> m or kg m-2] - real, dimension(:,:,:), target, intent(in ) :: T !< the model temperature array [C ~> degC] - real, dimension(:,:,:), target, intent(in ) :: S !< the model salinity array [S ~> ppt] - type(EOS_type), target, intent(in ) :: eqn_of_state !< Equation of state structure + type(thermo_var_ptrs), target, intent(in ) :: tv !< A sturcture with thermodynamic variables that are + !! are used to convert thicknesses to vertical extents type(diag_ctrl), intent(inout) :: diag_cs !< diag mediator control structure ! Keep pointers to h, T, S needed for the diagnostic remapping diag_cs%h => h - diag_cs%T => T - diag_cs%S => S - diag_cs%eqn_of_state => eqn_of_state + diag_cs%T => tv%T + diag_cs%S => tv%S + diag_cs%eqn_of_state => tv%eqn_of_state + diag_cs%tv => tv end subroutine @@ -3374,11 +3430,15 @@ subroutine diag_update_remap_grids(diag_cs, alt_h, alt_T, alt_S, update_intensiv logical, optional, intent(in ) :: update_extensive !< If true (not default), update the grids used for !! intensive diagnostics ! Local variables - integer :: i + integer :: m real, dimension(:,:,:), pointer :: h_diag => NULL() ! The layer thickneses for diagnostics [H ~> m or kg m-2] real, dimension(:,:,:), pointer :: T_diag => NULL() ! The layer temperatures for diagnostics [C ~> degC] real, dimension(:,:,:), pointer :: S_diag => NULL() ! The layer salinities for diagnostics [S ~> ppt] - logical :: update_intensive_local, update_extensive_local + real, dimension(diag_cs%G%isd:diag_cS%G%ied, diag_cs%G%jsd:diag_cS%G%jed, diag_cs%GV%ke) :: & + dz_diag ! Layer vertical extents for remapping [Z ~> m] + logical :: update_intensive_local, update_extensive_local, dz_diag_needed + + if (diag_cs%show_call_tree) call callTree_enter("diag_update_remap_grids()") ! Set values based on optional input arguments if (present(alt_h)) then @@ -3415,17 +3475,38 @@ subroutine diag_update_remap_grids(diag_cs, alt_h, alt_T, alt_S, update_intensiv "diagnostic structure have been overridden") endif + ! Determine the diagnostic grid spacing in height units, if it is needed. + dz_diag_needed = .false. + if (update_intensive_local .or. update_extensive_local) then + do m=1, diag_cs%num_diag_coords + if (diag_cs%diag_remap_cs(m)%Z_based_coord) dz_diag_needed = .true. + enddo + endif + if (dz_diag_needed) then + call thickness_to_dz(h_diag, diag_cs%tv, dz_diag, diag_cs%G, diag_cs%GV, diag_cs%US, halo_size=1) + endif + if (update_intensive_local) then - do i=1, diag_cs%num_diag_coords - call diag_remap_update(diag_cs%diag_remap_cs(i), diag_cs%G, diag_cs%GV, diag_cs%US, h_diag, T_diag, S_diag, & - diag_cs%eqn_of_state, diag_cs%diag_remap_cs(i)%h) + do m=1, diag_cs%num_diag_coords + if (diag_cs%diag_remap_cs(m)%Z_based_coord) then + call diag_remap_update(diag_cs%diag_remap_cs(m), diag_cs%G, diag_cs%GV, diag_cs%US, dz_diag, T_diag, S_diag, & + diag_cs%eqn_of_state, diag_cs%diag_remap_cs(m)%h) + else + call diag_remap_update(diag_cs%diag_remap_cs(m), diag_cs%G, diag_cs%GV, diag_cs%US, h_diag, T_diag, S_diag, & + diag_cs%eqn_of_state, diag_cs%diag_remap_cs(m)%h) + endif enddo endif if (update_extensive_local) then diag_cs%h_begin(:,:,:) = diag_cs%h(:,:,:) - do i=1, diag_cs%num_diag_coords - call diag_remap_update(diag_cs%diag_remap_cs(i), diag_cs%G, diag_cs%GV, diag_cs%US, h_diag, T_diag, S_diag, & - diag_cs%eqn_of_state, diag_cs%diag_remap_cs(i)%h_extensive) + do m=1, diag_cs%num_diag_coords + if (diag_cs%diag_remap_cs(m)%Z_based_coord) then + call diag_remap_update(diag_cs%diag_remap_cs(m), diag_cs%G, diag_cs%GV, diag_cs%US, dz_diag, T_diag, S_diag, & + diag_cs%eqn_of_state, diag_cs%diag_remap_cs(m)%h_extensive) + else + call diag_remap_update(diag_cs%diag_remap_cs(m), diag_cs%G, diag_cs%GV, diag_cs%US, h_diag, T_diag, S_diag, & + diag_cs%eqn_of_state, diag_cs%diag_remap_cs(m)%h_extensive) + endif enddo endif @@ -3437,6 +3518,8 @@ subroutine diag_update_remap_grids(diag_cs, alt_h, alt_T, alt_S, update_intensiv if (id_clock_diag_grid_updates>0) call cpu_clock_end(id_clock_diag_grid_updates) + if (diag_cs%show_call_tree) call callTree_leave("diag_update_remap_grids()") + end subroutine diag_update_remap_grids !> Sets up the 2d and 3d masks for native diagnostics @@ -4182,6 +4265,7 @@ subroutine downsample_field_3d(field_in, field_out, dl, method, mask, diag_cs, d allocate(field_out(1:f1,1:f2,ks:ke)) ! Fill the down sampled field on the down sampled diagnostics (almost always compuate) domain + !### The averaging used here is not rotationally invariant. if (method == MMM) then do k=ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d i0 = isv_o+dl*(i-isv_d) diff --git a/src/framework/MOM_diag_remap.F90 b/src/framework/MOM_diag_remap.F90 index cd9e22aae2..ace50242a5 100644 --- a/src/framework/MOM_diag_remap.F90 +++ b/src/framework/MOM_diag_remap.F90 @@ -103,7 +103,7 @@ module MOM_diag_remap !> Initialize a diagnostic remapping type with the given vertical coordinate. subroutine diag_remap_init(remap_cs, coord_tuple, answer_date, GV) - type(diag_remap_ctrl), intent(inout) :: remap_cs !< Diag remapping control structure + type(diag_remap_ctrl), intent(inout) :: remap_cs !< Diag remapping control structure character(len=*), intent(in) :: coord_tuple !< A string in form of !! MODULE_SUFFIX PARAMETER_SUFFIX COORDINATE_NAME integer, intent(in) :: answer_date !< The vintage of the order of arithmetic and expressions @@ -118,10 +118,11 @@ subroutine diag_remap_init(remap_cs, coord_tuple, answer_date, GV) remap_cs%vertical_coord_name = trim(extractWord(coord_tuple, 3)) remap_cs%vertical_coord = coordinateMode(remap_cs%vertical_coord_name) remap_cs%Z_based_coord = .false. - ! if ( & ! (.not.(GV%Boussinesq .or. GV%semi_Boussinesq)) .and. & - ! ((remap_cs%vertical_coord == coordinateMode('ZSTAR')) .or. & - ! (remap_cs%vertical_coord == coordinateMode('SIGMA'))) ) & - ! remap_cs%Z_based_coord = .true. + if (.not.(GV%Boussinesq .or. GV%semi_Boussinesq) .and. & + ((remap_cs%vertical_coord == coordinateMode('ZSTAR')) .or. & + (remap_cs%vertical_coord == coordinateMode('SIGMA')) .or. & + (remap_cs%vertical_coord == coordinateMode('RHO'))) ) & + remap_cs%Z_based_coord = .true. remap_cs%configured = .false. remap_cs%initialized = .false. @@ -295,7 +296,7 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state, h_targe enddo ; enddo else h_neglect = set_h_neglect(GV, remap_cs%answer_date, h_neglect_edge) - Z_unit_scale = GV%Z_to_H + Z_unit_scale = GV%Z_to_H ! This branch is not used in fully non-Boussinesq mode. do j=js-1,je+1 ; do i=is-1,ie+1 bottom_depth(i,j) = GV%Z_to_H * (G%bathyT(i,j) + G%Z_ref) enddo ; enddo From 611f575fbee7a27686f44832ea8d1197e8765386 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Tue, 6 Feb 2024 05:36:36 -0500 Subject: [PATCH 15/15] +(*)Fix downsample_mask indexing bugs Corrected indexing problems in downsample_mask that would cause masked reduced-resolution diagnostics to work improperly when the model is in symmetric memory mode or when global indexing is used. This involves passing the starting index in memory of the native-grid field being downsampled in all calls to downsample_mask. In global-indexing mode, this change avoids a series of segmentation faults that stopped the model runs when compiled for debugging. All solutions are bitwise identical in all cases but some down-scaled diagnostics will change in some memory modes, hopefully becoming consistent across all memory modes (although this has yet to be tested in all modes). --- src/framework/MOM_diag_mediator.F90 | 80 +++++++++++++++++------------ 1 file changed, 46 insertions(+), 34 deletions(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 6c90e26b98..eeb239859d 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -872,43 +872,51 @@ subroutine set_masks_for_axes_dsamp(G, diag_cs) do c=1, diag_cs%num_diag_coords ! Level/layer h-points in diagnostic coordinate axes => diag_cs%remap_axesTL(c) - call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesTL(c)%dsamp(dl)%mask3d, dl,G%isc, G%jsc, & + call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesTL(c)%dsamp(dl)%mask3d, & + dl, G%isc, G%jsc, G%isd, G%jsd, & G%HId2%isc, G%HId2%iec, G%HId2%jsc, G%HId2%jec, G%HId2%isd, G%HId2%ied, G%HId2%jsd, G%HId2%jed) diag_cs%dsamp(dl)%remap_axesTL(c)%mask3d => axes%mask3d !set non-downsampled mask ! Level/layer u-points in diagnostic coordinate axes => diag_cs%remap_axesCuL(c) - call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesCuL(c)%dsamp(dl)%mask3d, dl,G%IscB,G%JscB, & - G%HId2%IscB,G%HId2%IecB,G%HId2%jsc, G%HId2%jec,G%HId2%IsdB,G%HId2%IedB,G%HId2%jsd, G%HId2%jed) + call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesCuL(c)%dsamp(dl)%mask3d, & + dl, G%IscB, G%jsc, G%IsdB, G%jsd, & + G%HId2%IscB, G%HId2%IecB, G%HId2%jsc, G%HId2%jec, G%HId2%IsdB, G%HId2%IedB, G%HId2%jsd, G%HId2%jed) diag_cs%dsamp(dl)%remap_axesCul(c)%mask3d => axes%mask3d !set non-downsampled mask ! Level/layer v-points in diagnostic coordinate axes => diag_cs%remap_axesCvL(c) - call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesCvL(c)%dsamp(dl)%mask3d, dl,G%isc ,G%JscB, & - G%HId2%isc ,G%HId2%iec, G%HId2%JscB,G%HId2%JecB,G%HId2%isd ,G%HId2%ied, G%HId2%JsdB,G%HId2%JedB) + call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesCvL(c)%dsamp(dl)%mask3d, & + dl, G%isc, G%JscB, G%isd, G%JsdB, & + G%HId2%isc, G%HId2%iec, G%HId2%JscB, G%HId2%JecB, G%HId2%isd, G%HId2%ied, G%HId2%JsdB, G%HId2%JedB) diag_cs%dsamp(dl)%remap_axesCvL(c)%mask3d => axes%mask3d !set non-downsampled mask ! Level/layer q-points in diagnostic coordinate axes => diag_cs%remap_axesBL(c) - call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesBL(c)%dsamp(dl)%mask3d, dl,G%IscB,G%JscB, & - G%HId2%IscB,G%HId2%IecB,G%HId2%JscB,G%HId2%JecB,G%HId2%IsdB,G%HId2%IedB,G%HId2%JsdB,G%HId2%JedB) + call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesBL(c)%dsamp(dl)%mask3d, & + dl, G%IscB, G%JscB, G%IsdB, G%JsdB, & + G%HId2%IscB, G%HId2%IecB, G%HId2%JscB, G%HId2%JecB, G%HId2%IsdB, G%HId2%IedB, G%HId2%JsdB, G%HId2%JedB) diag_cs%dsamp(dl)%remap_axesBL(c)%mask3d => axes%mask3d !set non-downsampled mask ! Interface h-points in diagnostic coordinate (w-point) axes => diag_cs%remap_axesTi(c) - call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesTi(c)%dsamp(dl)%mask3d, dl,G%isc, G%jsc, & + call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesTi(c)%dsamp(dl)%mask3d, & + dl, G%isc, G%jsc, G%isd, G%jsd, & G%HId2%isc, G%HId2%iec, G%HId2%jsc, G%HId2%jec, G%HId2%isd, G%HId2%ied, G%HId2%jsd, G%HId2%jed) diag_cs%dsamp(dl)%remap_axesTi(c)%mask3d => axes%mask3d !set non-downsampled mask ! Interface u-points in diagnostic coordinate axes => diag_cs%remap_axesCui(c) - call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesCui(c)%dsamp(dl)%mask3d, dl,G%IscB,G%JscB, & - G%HId2%IscB,G%HId2%IecB,G%HId2%jsc, G%HId2%jec,G%HId2%IsdB,G%HId2%IedB,G%HId2%jsd, G%HId2%jed) + call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesCui(c)%dsamp(dl)%mask3d, & + dl, G%IscB, G%jsc, G%IsdB, G%jsd, & + G%HId2%IscB, G%HId2%IecB, G%HId2%jsc, G%HId2%jec, G%HId2%IsdB, G%HId2%IedB, G%HId2%jsd, G%HId2%jed) diag_cs%dsamp(dl)%remap_axesCui(c)%mask3d => axes%mask3d !set non-downsampled mask ! Interface v-points in diagnostic coordinate axes => diag_cs%remap_axesCvi(c) - call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesCvi(c)%dsamp(dl)%mask3d, dl,G%isc ,G%JscB, & - G%HId2%isc ,G%HId2%iec, G%HId2%JscB,G%HId2%JecB,G%HId2%isd ,G%HId2%ied, G%HId2%JsdB,G%HId2%JedB) + call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesCvi(c)%dsamp(dl)%mask3d, & + dl, G%isc, G%JscB, G%isd, G%JsdB, & + G%HId2%isc, G%HId2%iec, G%HId2%JscB, G%HId2%JecB, G%HId2%isd, G%HId2%ied, G%HId2%JsdB, G%HId2%JedB) diag_cs%dsamp(dl)%remap_axesCvi(c)%mask3d => axes%mask3d !set non-downsampled mask ! Interface q-points in diagnostic coordinate axes => diag_cs%remap_axesBi(c) - call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesBi(c)%dsamp(dl)%mask3d, dl,G%IscB,G%JscB, & - G%HId2%IscB,G%HId2%IecB,G%HId2%JscB,G%HId2%JecB,G%HId2%IsdB,G%HId2%IedB,G%HId2%JsdB,G%HId2%JedB) + call downsample_mask(axes%mask3d, diag_cs%dsamp(dl)%remap_axesBi(c)%dsamp(dl)%mask3d, & + dl, G%IscB, G%JscB, G%IsdB, G%JsdB, & + G%HId2%IscB, G%HId2%IecB, G%HId2%JscB, G%HId2%JecB, G%HId2%IsdB, G%HId2%IedB, G%HId2%JsdB, G%HId2%JedB) diag_cs%dsamp(dl)%remap_axesBi(c)%mask3d => axes%mask3d !set non-downsampled mask enddo enddo @@ -3998,13 +4006,13 @@ subroutine downsample_diag_masks_set(G, nz, diag_cs) do dl=2,MAX_DSAMP_LEV ! 2d mask - call downsample_mask(G%mask2dT, diag_cs%dsamp(dl)%mask2dT, dl,G%isc, G%jsc, & + call downsample_mask(G%mask2dT, diag_cs%dsamp(dl)%mask2dT, dl, G%isc, G%jsc, G%isd, G%jsd, & G%HId2%isc, G%HId2%iec, G%HId2%jsc, G%HId2%jec, G%HId2%isd, G%HId2%ied, G%HId2%jsd, G%HId2%jed) - call downsample_mask(G%mask2dBu,diag_cs%dsamp(dl)%mask2dBu, dl,G%IscB,G%JscB, & - G%HId2%IscB,G%HId2%IecB,G%HId2%JscB,G%HId2%JecB,G%HId2%IsdB,G%HId2%IedB,G%HId2%JsdB,G%HId2%JedB) - call downsample_mask(G%mask2dCu,diag_cs%dsamp(dl)%mask2dCu, dl,G%IscB,G%JscB, & - G%HId2%IscB,G%HId2%IecB,G%HId2%jsc, G%HId2%jec,G%HId2%IsdB,G%HId2%IedB,G%HId2%jsd, G%HId2%jed) - call downsample_mask(G%mask2dCv,diag_cs%dsamp(dl)%mask2dCv, dl,G%isc ,G%JscB, & + call downsample_mask(G%mask2dBu, diag_cs%dsamp(dl)%mask2dBu, dl,G%IscB, G%JscB, G%IsdB, G%JsdB, & + G%HId2%IscB,G%HId2%IecB, G%HId2%JscB,G%HId2%JecB,G%HId2%IsdB,G%HId2%IedB,G%HId2%JsdB,G%HId2%JedB) + call downsample_mask(G%mask2dCu, diag_cs%dsamp(dl)%mask2dCu, dl, G%IscB, G%jsc, G%IsdB, G%jsd, & + G%HId2%IscB,G%HId2%IecB, G%HId2%jsc, G%HId2%jec,G%HId2%IsdB,G%HId2%IedB,G%HId2%jsd, G%HId2%jed) + call downsample_mask(G%mask2dCv, diag_cs%dsamp(dl)%mask2dCv, dl,G %isc ,G%JscB, G%isd, G%JsdB, & G%HId2%isc ,G%HId2%iec, G%HId2%JscB,G%HId2%JecB,G%HId2%isd ,G%HId2%ied, G%HId2%JsdB,G%HId2%JedB) ! 3d native masks are needed by diag_manager but the native variables ! can only be masked 2d - for ocean points, all layers exists. @@ -4517,10 +4525,12 @@ end subroutine downsample_field_2d !> Allocate and compute the 2d down sampled mask !! The masks are down sampled based on a minority rule, i.e., a coarse cell is open (1) !! if at least one of the sub-cells are open, otherwise it's closed (0) -subroutine downsample_mask_2d(field_in, field_out, dl, isc_o, jsc_o, isc_d, iec_d, jsc_d, jec_d, & - isd_d, ied_d, jsd_d, jed_d) - real, dimension(:,:), intent(in) :: field_in !< Original field to be down sampled - real, dimension(:,:), pointer :: field_out !< Down sampled field +subroutine downsample_mask_2d(field_in, field_out, dl, isc_o, jsc_o, isd_o, jsd_o, & + isc_d, iec_d, jsc_d, jec_d, isd_d, ied_d, jsd_d, jed_d) + integer, intent(in) :: isd_o !< Original data domain i-start index + integer, intent(in) :: jsd_o !< Original data domain j-start index + real, dimension(isd_o:,jsd_o:), intent(in) :: field_in !< Original field to be down sampled in arbitrary units [A] + real, dimension(:,:), pointer :: field_out !< Down sampled field mask [nondim] integer, intent(in) :: dl !< Level of down sampling integer, intent(in) :: isc_o !< Original i-start index integer, intent(in) :: jsc_o !< Original j-start index @@ -4528,13 +4538,13 @@ subroutine downsample_mask_2d(field_in, field_out, dl, isc_o, jsc_o, isc_d, iec_ integer, intent(in) :: iec_d !< Computational i-end index of down sampled data integer, intent(in) :: jsc_d !< Computational j-start index of down sampled data integer, intent(in) :: jec_d !< Computational j-end index of down sampled data - integer, intent(in) :: isd_d !< Computational i-start index of down sampled data - integer, intent(in) :: ied_d !< Computational i-end index of down sampled data - integer, intent(in) :: jsd_d !< Computational j-start index of down sampled data - integer, intent(in) :: jed_d !< Computational j-end index of down sampled data + integer, intent(in) :: isd_d !< Data domain i-start index of down sampled data + integer, intent(in) :: ied_d !< Data domain i-end index of down sampled data + integer, intent(in) :: jsd_d !< Data domain j-start index of down sampled data + integer, intent(in) :: jed_d !< Data domain j-end index of down sampled data ! Locals integer :: i,j,ii,jj,i0,j0 - real :: tot_non_zero + real :: tot_non_zero ! The sum of values in the down-scaled cell [A] ! down sampled mask = 0 unless the mask value of one of the down sampling cells is 1 allocate(field_out(isd_d:ied_d,jsd_d:jed_d)) field_out(:,:) = 0.0 @@ -4552,10 +4562,12 @@ end subroutine downsample_mask_2d !> Allocate and compute the 3d down sampled mask !! The masks are down sampled based on a minority rule, i.e., a coarse cell is open (1) !! if at least one of the sub-cells are open, otherwise it's closed (0) -subroutine downsample_mask_3d(field_in, field_out, dl, isc_o, jsc_o, isc_d, iec_d, jsc_d, jec_d, & - isd_d, ied_d, jsd_d, jed_d) - real, dimension(:,:,:), intent(in) :: field_in !< Original field to be down sampled - real, dimension(:,:,:), pointer :: field_out !< down sampled field +subroutine downsample_mask_3d(field_in, field_out, dl, isc_o, jsc_o, isd_o, jsd_o, & + isc_d, iec_d, jsc_d, jec_d, isd_d, ied_d, jsd_d, jed_d) + integer, intent(in) :: isd_o !< Original data domain i-start index + integer, intent(in) :: jsd_o !< Original data domain j-start index + real, dimension(isd_o:,jsd_o:,:), intent(in) :: field_in !< Original field to be down sampled in arbitrary units [A] + real, dimension(:,:,:), pointer :: field_out !< down sampled field mask [nondim] integer, intent(in) :: dl !< Level of down sampling integer, intent(in) :: isc_o !< Original i-start index integer, intent(in) :: jsc_o !< Original j-start index @@ -4569,7 +4581,7 @@ subroutine downsample_mask_3d(field_in, field_out, dl, isc_o, jsc_o, isc_d, iec_ integer, intent(in) :: jed_d !< Computational j-end index of down sampled data ! Locals integer :: i,j,ii,jj,i0,j0,k,ks,ke - real :: tot_non_zero + real :: tot_non_zero ! The sum of values in the down-scaled cell [A] ! down sampled mask = 0 unless the mask value of one of the down sampling cells is 1 ks = lbound(field_in,3) ; ke = ubound(field_in,3) allocate(field_out(isd_d:ied_d,jsd_d:jed_d,ks:ke))