diff --git a/.github/workflows/documentation-and-style.yml b/.github/workflows/documentation-and-style.yml index 3988db7675..c83de48159 100644 --- a/.github/workflows/documentation-and-style.yml +++ b/.github/workflows/documentation-and-style.yml @@ -18,11 +18,11 @@ jobs: continue-on-error: true - name: Install packages used when generating documentation - run: > - sudo apt-get install - python3-sphinx python3-lxml perl - texlive-binaries texlive-base bibtool tex-common texlive-bibtex-extra - graphviz + run: | + sudo apt-get update + sudo apt-get install python3-sphinx python3-lxml perl + sudo apt-get install texlive-binaries texlive-base bibtool tex-common texlive-bibtex-extra + sudo apt-get install graphviz - name: Build doxygen HTML run: | diff --git a/ac/configure.ac b/ac/configure.ac index 14249fc09a..b069fdd56f 100644 --- a/ac/configure.ac +++ b/ac/configure.ac @@ -82,12 +82,20 @@ AX_FC_CHECK_MODULE([mpi], # netCDF configuration -AC_PATH_PROG([NC_CONFIG], [nc-config]) -AS_IF([test -n "$NC_CONFIG"], - [CPPFLAGS="$CPPFLAGS -I$($NC_CONFIG --includedir)" - FCFLAGS="$FCFLAGS -I$($NC_CONFIG --includedir)" - LDFLAGS="$LDFLAGS -L$($NC_CONFIG --libdir)"], - [AC_MSG_ERROR([Could not find nc-config.])]) + +# NOTE: `nf-config --flibs` combines library paths (-L) and libraries (-l), +# even though these ought to be separated in the invocation of `ld`. +# +# We use `sed` to strip the -l and pass the -L to LDFLAGS, and rely on autoconf +# to configure the -l flags. +AC_PROG_SED + +AC_PATH_PROG([NF_CONFIG], [nf-config]) +AS_IF([test -n "$NF_CONFIG"], + [CPPFLAGS="$CPPFLAGS $($NF_CONFIG --fflags)" + FCFLAGS="$FCFLAGS $($NF_CONFIG --fflags)" + LDFLAGS="$LDFLAGS $($NF_CONFIG --flibs | $SED -e 's/-l[[^ ]]*//g')"], + [AC_MSG_ERROR([Could not find nf-config.])]) AX_FC_CHECK_MODULE([netcdf], [], [AC_MSG_ERROR([Could not find netcdf module.])]) diff --git a/ac/deps/configure.fms.ac b/ac/deps/configure.fms.ac index 714f20df6b..4ac86f6445 100644 --- a/ac/deps/configure.fms.ac +++ b/ac/deps/configure.fms.ac @@ -73,12 +73,20 @@ AC_DEFINE([use_libMPI]) # netCDF configuration -AC_PATH_PROG([NC_CONFIG], [nc-config]) -AS_IF([test -n "$NC_CONFIG"], - [CPPFLAGS="$CPPFLAGS -I$($NC_CONFIG --includedir)" - FCFLAGS="$FCFLAGS -I$($NC_CONFIG --includedir)" - LDFLAGS="$LDFLAGS -L$($NC_CONFIG --libdir)"], - [AC_MSG_ERROR([Could not find nc-config.])]) + +# NOTE: `nf-config --flibs` combines library paths (-L) and libraries (-l), +# even though these ought to be separated in the invocation of `ld`. +# +# We use `sed` to strip the -l and pass the -L to LDFLAGS, and rely on autoconf +# to configure the -l flags. +AC_PROG_SED + +AC_PATH_PROG([NF_CONFIG], [nf-config]) +AS_IF([test -n "$NF_CONFIG"], + [CPPFLAGS="$CPPFLAGS $($NF_CONFIG --fflags)" + FCFLAGS="$FCFLAGS $($NF_CONFIG --fflags)" + LDFLAGS="$LDFLAGS $($NF_CONFIG --flibs | $SED -e 's/-l[[^ ]]*//g')"], + [AC_MSG_ERROR([Could not find nf-config.])]) AX_FC_CHECK_MODULE([netcdf], [], [AC_MSG_ERROR([Could not find netcdf module.])]) diff --git a/config_src/solo_driver/MOM_driver.F90 b/config_src/solo_driver/MOM_driver.F90 index 063da0976d..9981f291b1 100644 --- a/config_src/solo_driver/MOM_driver.F90 +++ b/config_src/solo_driver/MOM_driver.F90 @@ -261,7 +261,7 @@ program MOM_main !$ call omp_set_num_threads(ocean_nthreads) !$OMP PARALLEL !$ write(6,*) "ocean_solo OMPthreading ", fms_affinity_get(), omp_get_thread_num(), omp_get_num_threads() -!$ call flush(6) +!$ flush(6) !$OMP END PARALLEL ! Read ocean_solo restart, which can override settings from the namelist. diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index d11c1ae985..99268747e1 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -1875,17 +1875,16 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, handles%id_total_saltflux = register_scalar_field('ocean_model', & 'total_salt_flux', Time, diag, & - long_name='Area integrated surface salt flux', units='kg', & + long_name='Area integrated surface salt flux', units='kg s-1', & cmor_field_name='total_sfdsi', & - cmor_units='kg s-1', & cmor_standard_name='downward_sea_ice_basal_salt_flux_area_integrated',& cmor_long_name='Downward Sea Ice Basal Salt Flux Area Integrated') handles%id_total_saltFluxIn = register_scalar_field('ocean_model', 'total_salt_Flux_In', & - Time, diag, long_name='Area integrated surface salt flux at surface from coupler', units='kg') + Time, diag, long_name='Area integrated surface salt flux at surface from coupler', units='kg s-1') handles%id_total_saltFluxAdded = register_scalar_field('ocean_model', 'total_salt_Flux_Added', & - Time, diag, long_name='Area integrated surface salt flux due to restoring or flux adjustment', units='kg') + Time, diag, long_name='Area integrated surface salt flux due to restoring or flux adjustment', units='kg s-1') end subroutine register_forcing_type_diags diff --git a/src/diagnostics/MOM_PointAccel.F90 b/src/diagnostics/MOM_PointAccel.F90 index f6326b06fa..091f88fac2 100644 --- a/src/diagnostics/MOM_PointAccel.F90 +++ b/src/diagnostics/MOM_PointAccel.F90 @@ -390,7 +390,7 @@ subroutine write_u_accel(I, j, um, hin, ADp, CDp, dt_in_T, G, GV, US, CS, vel_rp write(file,'(2/)') - call flush(file) + flush(file) endif end subroutine write_u_accel @@ -722,7 +722,7 @@ subroutine write_v_accel(i, J, vm, hin, ADp, CDp, dt_in_T, G, GV, US, CS, vel_rp write(file,'(2/)') - call flush(file) + flush(file) endif end subroutine write_v_accel diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 3936a788d0..1f55801064 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -2225,6 +2225,7 @@ subroutine set_dependent_diagnostics(MIS, ADp, CDp, G, CS) if (associated(CS%KE_dia)) then call safe_alloc_ptr(ADp%du_dt_dia,IsdB,IedB,jsd,jed,nz) call safe_alloc_ptr(ADp%dv_dt_dia,isd,ied,JsdB,JedB,nz) + call safe_alloc_ptr(CDp%diapyc_vel,isd,ied,jsd,jed,nz+1) endif if (associated(CS%uhGM_Rlay)) call safe_alloc_ptr(CDp%uhGM,IsdB,IedB,jsd,jed,nz) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 28c4c867d7..e04234c859 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -242,7 +242,7 @@ module MOM_diag_mediator !! This file is open if available_diag_doc_unit is > 0. 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 :: diag_as_chksum !< If true, log chksums in a text file instead of posting diagnostics 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 @@ -606,7 +606,7 @@ subroutine set_axes_info_dsamp(G, GV, param_file, diag_cs, id_zl_native, id_zi_n id_zl = id_zl_native ; id_zi = id_zi_native !Axes group for native downsampled diagnostics do dl=2,MAX_DSAMP_LEV - if(dl .ne. 2) call MOM_error(FATAL, "set_axes_info_dsamp: Downsample level other than 2 is not supported yet!") + if (dl /= 2) call MOM_error(FATAL, "set_axes_info_dsamp: Downsample level other than 2 is not supported yet!") if (G%symmetric) then allocate(gridLonB_dsamp(diag_cs%dsamp(dl)%isgB:diag_cs%dsamp(dl)%iegB)) allocate(gridLatB_dsamp(diag_cs%dsamp(dl)%jsgB:diag_cs%dsamp(dl)%jegB)) @@ -871,7 +871,7 @@ subroutine set_masks_for_axes_dsamp(G, diag_cs) !The downsampled mask is needed for sending out the diagnostics output via diag_manager !The non-downsampled mask is needed for downsampling the diagnostics field do dl=2,MAX_DSAMP_LEV - if(dl .ne. 2) call MOM_error(FATAL, "set_masks_for_axes_dsamp: Downsample level other than 2 is not supported!") + if (dl /= 2) call MOM_error(FATAL, "set_masks_for_axes_dsamp: Downsample level other than 2 is not supported!") do c=1, diag_cs%num_diag_coords ! Level/layer h-points in diagnostic coordinate axes => diag_cs%remap_axesTL(c) @@ -1281,7 +1281,7 @@ subroutine post_data_0d(diag_field_id, field, diag_cs, is_static) if (diag_cs%diag_as_chksum) then call chksum0(locfield, diag%debug_str, logunit=diag_cs%chksum_iounit) - else if (is_stat) then + elseif (is_stat) then used = send_data(diag%fms_diag_id, locfield) elseif (diag_cs%ave_enabled) then used = send_data(diag%fms_diag_id, locfield, diag_cs%time_end) @@ -1333,7 +1333,7 @@ subroutine post_data_1d_k(diag_field_id, field, diag_cs, is_static) if (diag_cs%diag_as_chksum) then call zchksum(locfield, diag%debug_str, logunit=diag_cs%chksum_iounit) - else if (is_stat) then + elseif (is_stat) then used = send_data(diag%fms_diag_id, locfield) elseif (diag_cs%ave_enabled) then used = send_data(diag%fms_diag_id, locfield, diag_cs%time_end, weight=diag_cs%time_int) @@ -1448,38 +1448,38 @@ subroutine post_data_2d_low(diag, field, diag_cs, is_static, mask) endif if (present(mask)) then - locmask => mask - elseif(.NOT. is_stat) then - if(associated(diag%axes%mask2d)) locmask => diag%axes%mask2d + locmask => mask + elseif (.NOT. is_stat) then + if (associated(diag%axes%mask2d)) locmask => diag%axes%mask2d endif dl=1 - if(.NOT. is_stat) dl = diag%axes%downsample_level !static field downsample i not supported yet + if (.NOT. is_stat) dl = diag%axes%downsample_level !static field downsample i not supported yet !Downsample the diag field and mask (if present) if (dl > 1) then - isv_o=isv ; jsv_o=jsv - call downsample_diag_field(locfield, locfield_dsamp, dl, diag_cs, diag,isv,iev,jsv,jev, mask) - if ((diag%conversion_factor /= 0.) .and. (diag%conversion_factor /= 1.)) deallocate( locfield ) - locfield => locfield_dsamp - if (present(mask)) then - call downsample_field_2d(locmask, locmask_dsamp, dl, MSK, locmask, diag_cs,diag,isv_o,jsv_o,isv,iev,jsv,jev) - locmask => locmask_dsamp - elseif(associated(diag%axes%dsamp(dl)%mask2d)) then - locmask => diag%axes%dsamp(dl)%mask2d - endif + isv_o = isv ; jsv_o = jsv + call downsample_diag_field(locfield, locfield_dsamp, dl, diag_cs, diag,isv,iev,jsv,jev, mask) + if ((diag%conversion_factor /= 0.) .and. (diag%conversion_factor /= 1.)) deallocate( locfield ) + locfield => locfield_dsamp + if (present(mask)) then + call downsample_field_2d(locmask, locmask_dsamp, dl, MSK, locmask, diag_cs,diag,isv_o,jsv_o,isv,iev,jsv,jev) + locmask => locmask_dsamp + elseif (associated(diag%axes%dsamp(dl)%mask2d)) then + locmask => diag%axes%dsamp(dl)%mask2d + endif endif if (diag_cs%diag_as_chksum) then if (diag%axes%is_h_point) then call hchksum(locfield, diag%debug_str, diag_cs%G%HI, & logunit=diag_cs%chksum_iounit) - else if (diag%axes%is_u_point) then + elseif (diag%axes%is_u_point) then call uchksum(locfield, diag%debug_str, diag_cs%G%HI, & logunit=diag_cs%chksum_iounit) - else if (diag%axes%is_v_point) then + elseif (diag%axes%is_v_point) then call vchksum(locfield, diag%debug_str, diag_cs%G%HI, & logunit=diag_cs%chksum_iounit) - else if (diag%axes%is_q_point) then + elseif (diag%axes%is_q_point) then call Bchksum(locfield, diag%debug_str, diag_cs%G%HI, & logunit=diag_cs%chksum_iounit) else @@ -1735,25 +1735,25 @@ subroutine post_data_3d_low(diag, field, diag_cs, is_static, mask) endif if (present(mask)) then - locmask => mask - elseif(associated(diag%axes%mask3d)) then - locmask => diag%axes%mask3d + locmask => mask + elseif (associated(diag%axes%mask3d)) then + locmask => diag%axes%mask3d endif dl=1 - if(.NOT. is_stat) dl = diag%axes%downsample_level !static field downsample i not supported yet + if (.NOT. is_stat) dl = diag%axes%downsample_level !static field downsample i not supported yet !Downsample the diag field and mask (if present) if (dl > 1) then - isv_o=isv ; jsv_o=jsv - call downsample_diag_field(locfield, locfield_dsamp, dl, diag_cs, diag,isv,iev,jsv,jev, mask) - if ((diag%conversion_factor /= 0.) .and. (diag%conversion_factor /= 1.)) deallocate( locfield ) - locfield => locfield_dsamp - if (present(mask)) then - call downsample_field_3d(locmask, locmask_dsamp, dl, MSK, locmask, diag_cs,diag,isv_o,jsv_o,isv,iev,jsv,jev) - locmask => locmask_dsamp - elseif(associated(diag%axes%dsamp(dl)%mask3d)) then - locmask => diag%axes%dsamp(dl)%mask3d - endif + isv_o = isv ; jsv_o = jsv + call downsample_diag_field(locfield, locfield_dsamp, dl, diag_cs, diag,isv,iev,jsv,jev, mask) + if ((diag%conversion_factor /= 0.) .and. (diag%conversion_factor /= 1.)) deallocate( locfield ) + locfield => locfield_dsamp + if (present(mask)) then + call downsample_field_3d(locmask, locmask_dsamp, dl, MSK, locmask, diag_cs,diag,isv_o,jsv_o,isv,iev,jsv,jev) + locmask => locmask_dsamp + elseif (associated(diag%axes%dsamp(dl)%mask3d)) then + locmask => diag%axes%dsamp(dl)%mask3d + endif endif if (diag%fms_diag_id>0) then @@ -1761,13 +1761,13 @@ subroutine post_data_3d_low(diag, field, diag_cs, is_static, mask) if (diag%axes%is_h_point) then call hchksum(locfield, diag%debug_str, diag_cs%G%HI, & logunit=diag_cs%chksum_iounit) - else if (diag%axes%is_u_point) then + elseif (diag%axes%is_u_point) then call uchksum(locfield, diag%debug_str, diag_cs%G%HI, & logunit=diag_cs%chksum_iounit) - else if (diag%axes%is_v_point) then + elseif (diag%axes%is_v_point) then call vchksum(locfield, diag%debug_str, diag_cs%G%HI, & logunit=diag_cs%chksum_iounit) - else if (diag%axes%is_q_point) then + elseif (diag%axes%is_q_point) then call Bchksum(locfield, diag%debug_str, diag_cs%G%HI, & logunit=diag_cs%chksum_iounit) else @@ -1941,12 +1941,12 @@ integer function register_diag_field(module_name, field_name, axes_in, init_time verbose, do_not_log, err_msg, interp_method, tile_count, cmor_field_name, & cmor_long_name, cmor_units, cmor_standard_name, cell_methods, & x_cell_method, y_cell_method, v_cell_method, conversion, v_extensive) - character(len=*), intent(in) :: module_name !< Name of this module, usually "ocean_model" - !! or "ice_shelf_model" - character(len=*), intent(in) :: field_name !< Name of the diagnostic field - type(axes_grp), target, intent(in) :: axes_in !< Container w/ up to 3 integer handles that - !! indicates axes for this field - type(time_type), intent(in) :: init_time !< Time at which a field is first available? + character(len=*), intent(in) :: module_name !< Name of this module, usually "ocean_model" + !! or "ice_shelf_model" + character(len=*), intent(in) :: field_name !< Name of the diagnostic field + type(axes_grp), target, intent(in) :: axes_in !< Container w/ up to 3 integer handles that + !! indicates axes for this field + type(time_type), intent(in) :: init_time !< Time at which a field is first available? character(len=*), optional, intent(in) :: long_name !< Long name of a field. character(len=*), optional, intent(in) :: units !< Units of a field. character(len=*), optional, intent(in) :: standard_name !< Standardized name associated with a field @@ -1984,7 +1984,10 @@ integer function register_diag_field(module_name, field_name, axes_in, init_time type(axes_grp), pointer :: remap_axes => null() type(axes_grp), pointer :: axes => null() integer :: dm_id, i, dl + character(len=256) :: msg, cm_string character(len=256) :: new_module_name + character(len=480) :: module_list, var_list + integer :: num_modnm, num_varnm logical :: active axes => axes_in @@ -1995,23 +1998,26 @@ integer function register_diag_field(module_name, field_name, axes_in, init_time dm_id = -1 if (axes_in%id == diag_cs%axesTL%id) then - axes => diag_cs%axesTL + axes => diag_cs%axesTL elseif (axes_in%id == diag_cs%axesBL%id) then - axes => diag_cs%axesBL + axes => diag_cs%axesBL elseif (axes_in%id == diag_cs%axesCuL%id ) then - axes => diag_cs%axesCuL + axes => diag_cs%axesCuL elseif (axes_in%id == diag_cs%axesCvL%id) then - axes => diag_cs%axesCvL + axes => diag_cs%axesCvL elseif (axes_in%id == diag_cs%axesTi%id) then - axes => diag_cs%axesTi + axes => diag_cs%axesTi elseif (axes_in%id == diag_cs%axesBi%id) then - axes => diag_cs%axesBi + axes => diag_cs%axesBi elseif (axes_in%id == diag_cs%axesCui%id ) then - axes => diag_cs%axesCui + axes => diag_cs%axesCui elseif (axes_in%id == diag_cs%axesCvi%id) then - axes => diag_cs%axesCvi + axes => diag_cs%axesCvi endif + module_list = "{"//trim(module_name) + num_modnm = 1 + ! Register the native diagnostic active = register_diag_field_expand_cmor(dm_id, module_name, field_name, axes, & init_time, long_name=long_name, units=units, missing_value=MOM_missing_value, & @@ -2023,6 +2029,21 @@ integer function register_diag_field(module_name, field_name, axes_in, init_time cell_methods=cell_methods, x_cell_method=x_cell_method, & y_cell_method=y_cell_method, v_cell_method=v_cell_method, & conversion=conversion, v_extensive=v_extensive) + if (associated(axes%xyave_axes)) then + num_varnm = 2 ; var_list = "{"//trim(field_name)//","//trim(field_name)//"_xyave" + else + num_varnm = 1 ; var_list = "{"//trim(field_name) + endif + if (present(cmor_field_name)) then + if (associated(axes%xyave_axes)) then + num_varnm = num_varnm + 2 + var_list = trim(var_list)//","//trim(cmor_field_name)//","//trim(cmor_field_name)//"_xyave" + else + num_varnm = num_varnm + 1 + var_list = trim(var_list)//","//trim(cmor_field_name) + endif + endif + var_list = trim(var_list)//"}" ! For each diagnostic coordinate register the diagnostic again under a different module name do i=1,diag_cs%num_diag_coords @@ -2032,21 +2053,21 @@ integer function register_diag_field(module_name, field_name, axes_in, init_time if (axes_in%rank == 3) then remap_axes => null() if ((axes_in%id == diag_cs%axesTL%id)) then - remap_axes => diag_cs%remap_axesTL(i) + remap_axes => diag_cs%remap_axesTL(i) elseif (axes_in%id == diag_cs%axesBL%id) then - remap_axes => diag_cs%remap_axesBL(i) + remap_axes => diag_cs%remap_axesBL(i) elseif (axes_in%id == diag_cs%axesCuL%id ) then - remap_axes => diag_cs%remap_axesCuL(i) + remap_axes => diag_cs%remap_axesCuL(i) elseif (axes_in%id == diag_cs%axesCvL%id) then - remap_axes => diag_cs%remap_axesCvL(i) + remap_axes => diag_cs%remap_axesCvL(i) elseif (axes_in%id == diag_cs%axesTi%id) then - remap_axes => diag_cs%remap_axesTi(i) + remap_axes => diag_cs%remap_axesTi(i) elseif (axes_in%id == diag_cs%axesBi%id) then - remap_axes => diag_cs%remap_axesBi(i) + remap_axes => diag_cs%remap_axesBi(i) elseif (axes_in%id == diag_cs%axesCui%id ) then - remap_axes => diag_cs%remap_axesCui(i) + remap_axes => diag_cs%remap_axesCui(i) elseif (axes_in%id == diag_cs%axesCvi%id) then - remap_axes => diag_cs%remap_axesCvi(i) + remap_axes => diag_cs%remap_axesCvi(i) endif ! When the MOM_diag_to_Z module has been obsoleted we can assume remap_axes will ! always exist but in the mean-time we have to do this check: @@ -2066,6 +2087,8 @@ integer function register_diag_field(module_name, field_name, axes_in, init_time if (active) then call diag_remap_set_active(diag_cs%diag_remap_cs(i)) endif + module_list = trim(module_list)//","//trim(new_module_name) + num_modnm = num_modnm + 1 endif ! remap_axes%needs_remapping endif ! associated(remap_axes) endif ! axes%rank == 3 @@ -2073,46 +2096,46 @@ integer function register_diag_field(module_name, field_name, axes_in, init_time !Register downsampled diagnostics do dl=2,MAX_DSAMP_LEV - ! Do not attempt to checksum the downsampled diagnostics - if (diag_cs%diag_as_chksum) cycle + ! Do not attempt to checksum the downsampled diagnostics + if (diag_cs%diag_as_chksum) cycle - new_module_name = trim(module_name)//'_d2' + new_module_name = trim(module_name)//'_d2' - if (axes_in%rank == 3 .or. axes_in%rank == 2 ) then - axes => null() - if (axes_in%id == diag_cs%axesTL%id) then - axes => diag_cs%dsamp(dl)%axesTL - elseif (axes_in%id == diag_cs%axesBL%id) then - axes => diag_cs%dsamp(dl)%axesBL - elseif (axes_in%id == diag_cs%axesCuL%id ) then - axes => diag_cs%dsamp(dl)%axesCuL - elseif (axes_in%id == diag_cs%axesCvL%id) then - axes => diag_cs%dsamp(dl)%axesCvL - elseif (axes_in%id == diag_cs%axesTi%id) then - axes => diag_cs%dsamp(dl)%axesTi - elseif (axes_in%id == diag_cs%axesBi%id) then - axes => diag_cs%dsamp(dl)%axesBi - elseif (axes_in%id == diag_cs%axesCui%id ) then - axes => diag_cs%dsamp(dl)%axesCui - elseif (axes_in%id == diag_cs%axesCvi%id) then - axes => diag_cs%dsamp(dl)%axesCvi - elseif (axes_in%id == diag_cs%axesT1%id) then - axes => diag_cs%dsamp(dl)%axesT1 - elseif (axes_in%id == diag_cs%axesB1%id) then - axes => diag_cs%dsamp(dl)%axesB1 - elseif (axes_in%id == diag_cs%axesCu1%id ) then - axes => diag_cs%dsamp(dl)%axesCu1 - elseif (axes_in%id == diag_cs%axesCv1%id) then - axes => diag_cs%dsamp(dl)%axesCv1 - else - !Niki: Should we worry about these, e.g., diag_to_Z_CS? - call MOM_error(WARNING,"register_diag_field: Could not find a proper axes for " & - //trim( new_module_name)//"-"//trim(field_name)) - endif - endif - ! Register the native diagnostic - if (associated(axes)) then - active = register_diag_field_expand_cmor(dm_id, new_module_name, field_name, axes, & + if (axes_in%rank == 3 .or. axes_in%rank == 2 ) then + axes => null() + if (axes_in%id == diag_cs%axesTL%id) then + axes => diag_cs%dsamp(dl)%axesTL + elseif (axes_in%id == diag_cs%axesBL%id) then + axes => diag_cs%dsamp(dl)%axesBL + elseif (axes_in%id == diag_cs%axesCuL%id ) then + axes => diag_cs%dsamp(dl)%axesCuL + elseif (axes_in%id == diag_cs%axesCvL%id) then + axes => diag_cs%dsamp(dl)%axesCvL + elseif (axes_in%id == diag_cs%axesTi%id) then + axes => diag_cs%dsamp(dl)%axesTi + elseif (axes_in%id == diag_cs%axesBi%id) then + axes => diag_cs%dsamp(dl)%axesBi + elseif (axes_in%id == diag_cs%axesCui%id ) then + axes => diag_cs%dsamp(dl)%axesCui + elseif (axes_in%id == diag_cs%axesCvi%id) then + axes => diag_cs%dsamp(dl)%axesCvi + elseif (axes_in%id == diag_cs%axesT1%id) then + axes => diag_cs%dsamp(dl)%axesT1 + elseif (axes_in%id == diag_cs%axesB1%id) then + axes => diag_cs%dsamp(dl)%axesB1 + elseif (axes_in%id == diag_cs%axesCu1%id ) then + axes => diag_cs%dsamp(dl)%axesCu1 + elseif (axes_in%id == diag_cs%axesCv1%id) then + axes => diag_cs%dsamp(dl)%axesCv1 + else + !Niki: Should we worry about these, e.g., diag_to_Z_CS? + call MOM_error(WARNING,"register_diag_field: Could not find a proper axes for " & + //trim(new_module_name)//"-"//trim(field_name)) + endif + endif + ! Register the native diagnostic + if (associated(axes)) then + active = register_diag_field_expand_cmor(dm_id, new_module_name, field_name, axes, & init_time, long_name=long_name, units=units, missing_value=MOM_missing_value, & range=range, mask_variant=mask_variant, standard_name=standard_name, & verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, & @@ -2122,57 +2145,75 @@ integer function register_diag_field(module_name, field_name, axes_in, init_time cell_methods=cell_methods, x_cell_method=x_cell_method, & y_cell_method=y_cell_method, v_cell_method=v_cell_method, & conversion=conversion, v_extensive=v_extensive) - endif + module_list = trim(module_list)//","//trim(new_module_name) + num_modnm = num_modnm + 1 + endif + + ! For each diagnostic coordinate register the diagnostic again under a different module name + do i=1,diag_cs%num_diag_coords + new_module_name = trim(module_name)//'_'//trim(diag_cs%diag_remap_cs(i)%diag_module_suffix)//'_d2' + + ! Register diagnostics remapped to z vertical coordinate + if (axes_in%rank == 3) then + remap_axes => null() + if ((axes_in%id == diag_cs%axesTL%id)) then + remap_axes => diag_cs%dsamp(dl)%remap_axesTL(i) + elseif (axes_in%id == diag_cs%axesBL%id) then + remap_axes => diag_cs%dsamp(dl)%remap_axesBL(i) + elseif (axes_in%id == diag_cs%axesCuL%id ) then + remap_axes => diag_cs%dsamp(dl)%remap_axesCuL(i) + elseif (axes_in%id == diag_cs%axesCvL%id) then + remap_axes => diag_cs%dsamp(dl)%remap_axesCvL(i) + elseif (axes_in%id == diag_cs%axesTi%id) then + remap_axes => diag_cs%dsamp(dl)%remap_axesTi(i) + elseif (axes_in%id == diag_cs%axesBi%id) then + remap_axes => diag_cs%dsamp(dl)%remap_axesBi(i) + elseif (axes_in%id == diag_cs%axesCui%id ) then + remap_axes => diag_cs%dsamp(dl)%remap_axesCui(i) + elseif (axes_in%id == diag_cs%axesCvi%id) then + remap_axes => diag_cs%dsamp(dl)%remap_axesCvi(i) + endif - ! For each diagnostic coordinate register the diagnostic again under a different module name - do i=1,diag_cs%num_diag_coords - new_module_name = trim(module_name)//'_'//trim(diag_cs%diag_remap_cs(i)%diag_module_suffix)//'_d2' - - ! Register diagnostics remapped to z vertical coordinate - if (axes_in%rank == 3) then - remap_axes => null() - if ((axes_in%id == diag_cs%axesTL%id)) then - remap_axes => diag_cs%dsamp(dl)%remap_axesTL(i) - elseif (axes_in%id == diag_cs%axesBL%id) then - remap_axes => diag_cs%dsamp(dl)%remap_axesBL(i) - elseif (axes_in%id == diag_cs%axesCuL%id ) then - remap_axes => diag_cs%dsamp(dl)%remap_axesCuL(i) - elseif (axes_in%id == diag_cs%axesCvL%id) then - remap_axes => diag_cs%dsamp(dl)%remap_axesCvL(i) - elseif (axes_in%id == diag_cs%axesTi%id) then - remap_axes => diag_cs%dsamp(dl)%remap_axesTi(i) - elseif (axes_in%id == diag_cs%axesBi%id) then - remap_axes => diag_cs%dsamp(dl)%remap_axesBi(i) - elseif (axes_in%id == diag_cs%axesCui%id ) then - remap_axes => diag_cs%dsamp(dl)%remap_axesCui(i) - elseif (axes_in%id == diag_cs%axesCvi%id) then - remap_axes => diag_cs%dsamp(dl)%remap_axesCvi(i) - endif - - ! When the MOM_diag_to_Z module has been obsoleted we can assume remap_axes will - ! always exist but in the mean-time we have to do this check: - ! call assert(associated(remap_axes), 'register_diag_field: remap_axes not set') - if (associated(remap_axes)) then - if (remap_axes%needs_remapping .or. remap_axes%needs_interpolating) then - active = register_diag_field_expand_cmor(dm_id, new_module_name, field_name, remap_axes, & - init_time, long_name=long_name, units=units, missing_value=MOM_missing_value, & - range=range, mask_variant=mask_variant, standard_name=standard_name, & - verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, & - interp_method=interp_method, tile_count=tile_count, & - cmor_field_name=cmor_field_name, cmor_long_name=cmor_long_name, & - cmor_units=cmor_units, cmor_standard_name=cmor_standard_name, & - cell_methods=cell_methods, x_cell_method=x_cell_method, & - y_cell_method=y_cell_method, v_cell_method=v_cell_method, & - conversion=conversion, v_extensive=v_extensive) - if (active) then - call diag_remap_set_active(diag_cs%diag_remap_cs(i)) - endif - endif ! remap_axes%needs_remapping - endif ! associated(remap_axes) - endif ! axes%rank == 3 - enddo ! i + ! When the MOM_diag_to_Z module has been obsoleted we can assume remap_axes will + ! always exist but in the mean-time we have to do this check: + ! call assert(associated(remap_axes), 'register_diag_field: remap_axes not set') + if (associated(remap_axes)) then + if (remap_axes%needs_remapping .or. remap_axes%needs_interpolating) then + active = register_diag_field_expand_cmor(dm_id, new_module_name, field_name, remap_axes, & + init_time, long_name=long_name, units=units, missing_value=MOM_missing_value, & + range=range, mask_variant=mask_variant, standard_name=standard_name, & + verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, & + interp_method=interp_method, tile_count=tile_count, & + cmor_field_name=cmor_field_name, cmor_long_name=cmor_long_name, & + cmor_units=cmor_units, cmor_standard_name=cmor_standard_name, & + cell_methods=cell_methods, x_cell_method=x_cell_method, & + y_cell_method=y_cell_method, v_cell_method=v_cell_method, & + conversion=conversion, v_extensive=v_extensive) + if (active) then + call diag_remap_set_active(diag_cs%diag_remap_cs(i)) + endif + module_list = trim(module_list)//","//trim(new_module_name) + num_modnm = num_modnm + 1 + endif ! remap_axes%needs_remapping + endif ! associated(remap_axes) + endif ! axes%rank == 3 + enddo ! i enddo + if (is_root_pe() .and. (diag_CS%available_diag_doc_unit > 0)) then + msg = '' + if (present(cmor_field_name)) msg = 'CMOR equivalent is "'//trim(cmor_field_name)//'"' + call attach_cell_methods(-1, axes, cm_string, cell_methods, & + x_cell_method, y_cell_method, v_cell_method, & + v_extensive=v_extensive) + module_list = trim(module_list)//"}" + if (num_modnm <= 1) module_list = module_name + if (num_varnm <= 1) var_list = "" + + call log_available_diag(dm_id>0, module_list, field_name, cm_string, msg, diag_CS, & + long_name, units, standard_name, variants=var_list) + endif + register_diag_field = dm_id end function register_diag_field @@ -2244,12 +2285,6 @@ logical function register_diag_field_expand_cmor(dm_id, module_name, field_name, call attach_cell_methods(fms_id, axes, cm_string, cell_methods, & x_cell_method, y_cell_method, v_cell_method, & v_extensive=v_extensive) - if (is_root_pe() .and. diag_CS%available_diag_doc_unit > 0) then - msg = '' - if (present(cmor_field_name)) msg = 'CMOR equivalent is "'//trim(cmor_field_name)//'"' - call log_available_diag(fms_id>0, module_name, field_name, cm_string, & - msg, diag_CS, long_name, units, standard_name) - endif ! Associated horizontally area-averaged diagnostic fms_xyave_id = DIAG_FIELD_NOT_FOUND if (associated(axes%xyave_axes)) then @@ -2262,12 +2297,6 @@ logical function register_diag_field_expand_cmor(dm_id, module_name, field_name, if (.not. diag_cs%diag_as_chksum) & call attach_cell_methods(fms_xyave_id, axes%xyave_axes, cm_string, & cell_methods, v_cell_method, v_extensive=v_extensive) - if (is_root_pe() .and. diag_CS%available_diag_doc_unit > 0) then - msg = '' - if (present(cmor_field_name)) msg = 'CMOR equivalent is "'//trim(cmor_field_name)//'_xyave"' - call log_available_diag(fms_xyave_id>0, module_name, trim(field_name)//'_xyave', cm_string, & - msg, diag_CS, long_name, units, standard_name) - endif endif this_diag => null() if (fms_id /= DIAG_FIELD_NOT_FOUND .or. fms_xyave_id /= DIAG_FIELD_NOT_FOUND) then @@ -2306,12 +2335,6 @@ logical function register_diag_field_expand_cmor(dm_id, module_name, field_name, call attach_cell_methods(fms_id, axes, cm_string, & cell_methods, x_cell_method, y_cell_method, v_cell_method, & v_extensive=v_extensive) - if (is_root_pe() .and. diag_CS%available_diag_doc_unit > 0) then - msg = 'native name is "'//trim(field_name)//'"' - call log_available_diag(fms_id>0, module_name, cmor_field_name, cm_string, & - msg, diag_CS, posted_cmor_long_name, posted_cmor_units, & - posted_cmor_standard_name) - endif ! Associated horizontally area-averaged diagnostic fms_xyave_id = DIAG_FIELD_NOT_FOUND if (associated(axes%xyave_axes)) then @@ -2323,12 +2346,6 @@ logical function register_diag_field_expand_cmor(dm_id, module_name, field_name, err_msg=err_msg, interp_method=interp_method, tile_count=tile_count) call attach_cell_methods(fms_xyave_id, axes%xyave_axes, cm_string, & cell_methods, v_cell_method, v_extensive=v_extensive) - if (is_root_pe() .and. diag_CS%available_diag_doc_unit > 0) then - msg = 'native name is "'//trim(field_name)//'_xyave"' - call log_available_diag(fms_xyave_id>0, module_name, trim(cmor_field_name)//'_xyave', & - cm_string, msg, diag_CS, posted_cmor_long_name, posted_cmor_units, & - posted_cmor_standard_name) - endif endif this_diag => null() if (fms_id /= DIAG_FIELD_NOT_FOUND .or. fms_xyave_id /= DIAG_FIELD_NOT_FOUND) then @@ -2381,7 +2398,7 @@ integer function register_diag_field_expand_axes(module_name, field_name, axes, if (axes%diag_cs%diag_as_chksum) then fms_id = axes%diag_cs%num_chksum_diags + 1 axes%diag_cs%num_chksum_diags = fms_id - else if (present(interp_method) .or. axes%is_h_point) then + elseif (present(interp_method) .or. axes%is_h_point) then ! If interp_method is provided we must use it if (area_id>0) then if (volume_id>0) then @@ -2504,7 +2521,7 @@ subroutine add_xyz_method(diag, axes, x_cell_method, y_cell_method, v_cell_metho if (present(v_extensive)) then if (present(v_cell_method)) call MOM_error(FATAL, "attach_cell_methods: " // & 'Vertical cell method was specified along with the vertically extensive flag.') - if(v_extensive) then + if (v_extensive) then mstr='sum' else mstr='mean' @@ -2624,7 +2641,7 @@ subroutine attach_cell_methods(id, axes, ostring, cell_methods, & ostring = trim(adjustl(ostring))//' '//trim(axis_name)//':'//trim(v_cell_method) endif elseif (present(v_extensive)) then - if(v_extensive) then + if (v_extensive) then if (axes%rank==1) then call get_diag_axis_name(axes%handles(1), axis_name) elseif (axes%rank==3) then @@ -2744,12 +2761,13 @@ function register_scalar_field(module_name, field_name, init_time, diag_cs, & ! Document diagnostics in list of available diagnostics if (is_root_pe() .and. diag_CS%available_diag_doc_unit > 0) then - call log_available_diag(associated(diag), module_name, field_name, '', '', diag_CS, & - long_name, units, standard_name) if (present(cmor_field_name)) then - call log_available_diag(associated(cmor_diag), module_name, cmor_field_name, & - '', '', diag_CS, posted_cmor_long_name, posted_cmor_units, & - posted_cmor_standard_name) + call log_available_diag(associated(diag), module_name, field_name, '', '', diag_CS, & + long_name, units, standard_name, & + variants="{"//trim(field_name)//","//trim(cmor_field_name)//"}") + else + call log_available_diag(associated(diag), module_name, field_name, '', '', diag_CS, & + long_name, units, standard_name) endif endif @@ -2891,12 +2909,13 @@ function register_static_field(module_name, field_name, axes, & ! Document diagnostics in list of available diagnostics if (is_root_pe() .and. diag_CS%available_diag_doc_unit > 0) then - call log_available_diag(associated(diag), module_name, field_name, '', '', diag_CS, & - long_name, units, standard_name) if (present(cmor_field_name)) then - call log_available_diag(associated(cmor_diag), module_name, cmor_field_name, & - '', '', diag_CS, posted_cmor_long_name, posted_cmor_units, & - posted_cmor_standard_name) + call log_available_diag(associated(diag), module_name, field_name, '', '', diag_CS, & + long_name, units, standard_name, & + variants="{"//trim(field_name)//","//trim(cmor_field_name)//"}") + else + call log_available_diag(associated(diag), module_name, field_name, '', '', diag_CS, & + long_name, units, standard_name) endif endif @@ -2910,7 +2929,7 @@ subroutine describe_option(opt_name, value, diag_CS) character(len=*), intent(in) :: value !< A character string with the setting of the option. type(diag_ctrl), intent(in) :: diag_CS !< Structure used to regulate diagnostic output - character(len=240) :: mesg + character(len=480) :: mesg integer :: len_ind len_ind = len_trim(value) ! Add error handling for long values? @@ -2944,75 +2963,46 @@ function ocean_register_diag(var_desc, G, diag_CS, day) case ("L") select case (hor_grid) - case ("q") - axes => diag_cs%axesBL - case ("h") - axes => diag_cs%axesTL - case ("u") - axes => diag_cs%axesCuL - case ("v") - axes => diag_cs%axesCvL - case ("Bu") - axes => diag_cs%axesBL - case ("T") - axes => diag_cs%axesTL - case ("Cu") - axes => diag_cs%axesCuL - case ("Cv") - axes => diag_cs%axesCvL - case ("z") - axes => diag_cs%axeszL - case default - call MOM_error(FATAL, "ocean_register_diag: " // & - "unknown hor_grid component "//trim(hor_grid)) + case ("q") ; axes => diag_cs%axesBL + case ("h") ; axes => diag_cs%axesTL + case ("u") ; axes => diag_cs%axesCuL + case ("v") ; axes => diag_cs%axesCvL + case ("Bu") ; axes => diag_cs%axesBL + case ("T") ; axes => diag_cs%axesTL + case ("Cu") ; axes => diag_cs%axesCuL + case ("Cv") ; axes => diag_cs%axesCvL + case ("z") ; axes => diag_cs%axeszL + case default ; call MOM_error(FATAL, "ocean_register_diag: " // & + "unknown hor_grid component "//trim(hor_grid)) end select case ("i") select case (hor_grid) - case ("q") - axes => diag_cs%axesBi - case ("h") - axes => diag_cs%axesTi - case ("u") - axes => diag_cs%axesCui - case ("v") - axes => diag_cs%axesCvi - case ("Bu") - axes => diag_cs%axesBi - case ("T") - axes => diag_cs%axesTi - case ("Cu") - axes => diag_cs%axesCui - case ("Cv") - axes => diag_cs%axesCvi - case ("z") - axes => diag_cs%axeszi - case default - call MOM_error(FATAL, "ocean_register_diag: " // & - "unknown hor_grid component "//trim(hor_grid)) + case ("q") ; axes => diag_cs%axesBi + case ("h") ; axes => diag_cs%axesTi + case ("u") ; axes => diag_cs%axesCui + case ("v") ; axes => diag_cs%axesCvi + case ("Bu") ; axes => diag_cs%axesBi + case ("T") ; axes => diag_cs%axesTi + case ("Cu") ; axes => diag_cs%axesCui + case ("Cv") ; axes => diag_cs%axesCvi + case ("z") ; axes => diag_cs%axeszi + case default ; call MOM_error(FATAL, "ocean_register_diag: " // & + "unknown hor_grid component "//trim(hor_grid)) end select case ("1") select case (hor_grid) - case ("q") - axes => diag_cs%axesB1 - case ("h") - axes => diag_cs%axesT1 - case ("u") - axes => diag_cs%axesCu1 - case ("v") - axes => diag_cs%axesCv1 - case ("Bu") - axes => diag_cs%axesB1 - case ("T") - axes => diag_cs%axesT1 - case ("Cu") - axes => diag_cs%axesCu1 - case ("Cv") - axes => diag_cs%axesCv1 - case default - call MOM_error(FATAL, "ocean_register_diag: " // & - "unknown hor_grid component "//trim(hor_grid)) + case ("q") ; axes => diag_cs%axesB1 + case ("h") ; axes => diag_cs%axesT1 + case ("u") ; axes => diag_cs%axesCu1 + case ("v") ; axes => diag_cs%axesCv1 + case ("Bu") ; axes => diag_cs%axesB1 + case ("T") ; axes => diag_cs%axesT1 + case ("Cu") ; axes => diag_cs%axesCu1 + case ("Cv") ; axes => diag_cs%axesCv1 + case default ; call MOM_error(FATAL, "ocean_register_diag: " // & + "unknown hor_grid component "//trim(hor_grid)) end select case default @@ -3540,7 +3530,7 @@ end subroutine alloc_diag_with_id !> Log a diagnostic to the available diagnostics file. subroutine log_available_diag(used, module_name, field_name, cell_methods_string, comment, & - diag_CS, long_name, units, standard_name) + diag_CS, long_name, units, standard_name, variants) logical, intent(in) :: used !< Whether this diagnostic was in the diag_table or not character(len=*), intent(in) :: module_name !< Name of the diagnostic module character(len=*), intent(in) :: field_name !< Name of this diagnostic field @@ -3550,26 +3540,31 @@ subroutine log_available_diag(used, module_name, field_name, cell_methods_string character(len=*), optional, intent(in) :: long_name !< CF long name of diagnostic character(len=*), optional, intent(in) :: units !< Units for diagnostic character(len=*), optional, intent(in) :: standard_name !< CF standardized name of diagnostic + character(len=*), optional, intent(in) :: variants !< Alternate modules and variable names for + !! this diagnostic and derived diagnostics ! Local variables character(len=240) :: mesg if (used) then - mesg = '"'//trim(module_name)//'", "'//trim(field_name)//'" [Used]' + mesg = '"'//trim(field_name)//'" [Used]' else - mesg = '"'//trim(module_name)//'", "'//trim(field_name)//'" [Unused]' + mesg = '"'//trim(field_name)//'" [Unused]' endif if (len(trim((comment)))>0) then write(diag_CS%available_diag_doc_unit, '(a,x,"(",a,")")') trim(mesg),trim(comment) else write(diag_CS%available_diag_doc_unit, '(a)') trim(mesg) endif + call describe_option("modules", module_name, diag_CS) if (present(long_name)) call describe_option("long_name", long_name, diag_CS) if (present(units)) call describe_option("units", units, diag_CS) if (present(standard_name)) & call describe_option("standard_name", standard_name, diag_CS) if (len(trim((cell_methods_string)))>0) & call describe_option("cell_methods", trim(cell_methods_string), diag_CS) - + if (present(variants)) then ; if (len(trim(variants)) > 0) then + call describe_option("variants", variants, diag_CS) + endif ; endif end subroutine log_available_diag !> Log the diagnostic chksum to the chksum diag file @@ -3778,8 +3773,8 @@ subroutine downsample_diag_indices_get(fo1, fo2, dl, diag_cs, isv, iev, jsv, jev !We assume that the compute domain can be subdivided to dl*dl cells, hence avoiding the need of halo updates. !We want this check to error out only if there was a downsampled diagnostics requested and about to post that is !why the check is here and not in the init routines. This check need to be done only once, hence the outer if. - if(first_check) then - if(mod(diag_cs%ie-diag_cs%is+1, dl) .ne. 0 .OR. mod(diag_cs%je-diag_cs%js+1, dl) .ne. 0) then + if (first_check) then + if (mod(diag_cs%ie-diag_cs%is+1, dl) /= 0 .OR. mod(diag_cs%je-diag_cs%js+1, dl) /= 0) then write (mesg,*) "Non-commensurate downsampled domain is not supported. "//& "Please choose a layout such that NIGLOBAL/Layout_X and NJGLOBAL/Layout_Y are both divisible by dl=",dl,& " Current domain extents: ", diag_cs%is,diag_cs%ie, diag_cs%js,diag_cs%je @@ -3849,10 +3844,10 @@ subroutine downsample_diag_field_3d(locfield, locfield_dsamp, dl, diag_cs, diag, locmask => NULL() !Get the correct indices corresponding to input field !Shape of the input diag field - f1=size(locfield,1) - f2=size(locfield,2) + f1 = size(locfield,1) + f2 = size(locfield,2) !Save the extents of the original (fine) domain - isv_o=isv;jsv_o=jsv + isv_o = isv ; jsv_o = jsv !Get the shape of the downsampled field and overwrite isv,iev,jsv,jev with them call downsample_diag_indices_get(f1,f2, dl, diag_cs,isv,iev,jsv,jev) !Set the non-downsampled mask, it must be associated and initialized @@ -3890,19 +3885,19 @@ subroutine downsample_diag_field_2d(locfield, locfield_dsamp, dl, diag_cs, diag, locmask => NULL() !Get the correct indices corresponding to input field !Shape of the input diag field - f1=size(locfield,1) - f2=size(locfield,2) + f1 = size(locfield,1) + f2 = size(locfield,2) !Save the extents of the original (fine) domain - isv_o=isv;jsv_o=jsv + isv_o = isv ; jsv_o = jsv !Get the shape of the downsampled field and overwrite isv,iev,jsv,jev with them call downsample_diag_indices_get(f1,f2, dl, diag_cs,isv,iev,jsv,jev) !Set the non-downsampled mask, it must be associated and initialized if (present(mask)) then - locmask => mask + locmask => mask elseif (associated(diag%axes%mask2d)) then - locmask => diag%axes%mask2d + locmask => diag%axes%mask2d else - call MOM_error(FATAL, "downsample_diag_field_2d: Cannot downsample without a mask!!! ") + call MOM_error(FATAL, "downsample_diag_field_2d: Cannot downsample without a mask!!! ") endif call downsample_field(locfield, locfield_dsamp, dl, diag%xyz_method, locmask, diag_cs,diag, & @@ -3991,7 +3986,7 @@ subroutine downsample_field_3d(field_in, field_out, dl, method, mask, diag_cs, d ! Fill the down sampled field on the down sampled diagnostics (almost always compuate) domain if (method == MMM) then - do k= ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d + do k=ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d i0 = isv_o+dl*(i-isv_d) j0 = jsv_o+dl*(j-jsv_d) ave = 0.0 @@ -4001,22 +3996,22 @@ subroutine downsample_field_3d(field_in, field_out, dl, method, mask, diag_cs, d weight = mask(ii,jj,k) * diag_cs%G%areaT(ii,jj) * diag_cs%h(ii,jj,k) total_weight = total_weight + weight ave = ave+field_in(ii,jj,k) * weight - enddo; enddo + enddo ; enddo field_out(i,j,k) = ave/(total_weight + eps_vol) !Avoid zero mask at all aggregating cells where ave=0.0 - enddo; enddo; enddo + enddo ; enddo ; enddo elseif (method == SSS) then !e.g., volcello - do k= ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d + do k=ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d i0 = isv_o+dl*(i-isv_d) j0 = jsv_o+dl*(j-jsv_d) ave = 0.0 do jj=j0,j0+dl-1 ; do ii=i0,i0+dl-1 weight = mask(ii,jj,k) ave = ave+field_in(ii,jj,k)*weight - enddo; enddo + enddo ; enddo field_out(i,j,k) = ave !Masked Sum (total_weight=1) - enddo; enddo; enddo - elseif(method == MMP .or. method == MMS) then !e.g., T_advection_xy - do k= ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d + enddo ; enddo ; enddo + elseif (method == MMP .or. method == MMS) then !e.g., T_advection_xy + do k=ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d i0 = isv_o+dl*(i-isv_d) j0 = jsv_o+dl*(j-jsv_d) ave = 0.0 @@ -4026,49 +4021,49 @@ subroutine downsample_field_3d(field_in, field_out, dl, method, mask, diag_cs, d weight = mask(ii,jj,k) * diag_cs%G%areaT(ii,jj) total_weight = total_weight + weight ave = ave+field_in(ii,jj,k)*weight - enddo; enddo + enddo ; enddo field_out(i,j,k) = ave / (total_weight+eps_area) !Avoid zero mask at all aggregating cells where ave=0.0 - enddo; enddo; enddo - elseif(method == PMM) then - do k= ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d + enddo ; enddo ; enddo + elseif (method == PMM) then + do k=ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d i0 = isv_o+dl*(i-isv_d) j0 = jsv_o+dl*(j-jsv_d) ave = 0.0 total_weight = 0.0 ii=i0 do jj=j0,j0+dl-1 - weight =mask(ii,jj,k) * diag_cs%G%dyCu(ii,jj) * diag_cs%h(ii,jj,k) + weight = mask(ii,jj,k) * diag_cs%G%dyCu(ii,jj) * diag_cs%h(ii,jj,k) total_weight = total_weight +weight - ave=ave+field_in(ii,jj,k)*weight + ave = ave+field_in(ii,jj,k)*weight enddo field_out(i,j,k) = ave/(total_weight+eps_face) !Avoid zero mask at all aggregating cells where ave=0.0 - enddo; enddo; enddo - elseif(method == PSS) then !e.g. umo - do k= ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d + enddo ; enddo ; enddo + elseif (method == PSS) then !e.g. umo + do k=ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d i0 = isv_o+dl*(i-isv_d) j0 = jsv_o+dl*(j-jsv_d) ave = 0.0 ii=i0 do jj=j0,j0+dl-1 - weight =mask(ii,jj,k) - ave=ave+field_in(ii,jj,k)*weight + weight = mask(ii,jj,k) + ave = ave+field_in(ii,jj,k)*weight enddo field_out(i,j,k) = ave !Masked Sum (total_weight=1) - enddo; enddo; enddo - elseif(method == SPS) then !e.g. vmo - do k= ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d + enddo ; enddo ; enddo + elseif (method == SPS) then !e.g. vmo + do k=ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d i0 = isv_o+dl*(i-isv_d) j0 = jsv_o+dl*(j-jsv_d) ave = 0.0 jj=j0 do ii=i0,i0+dl-1 - weight =mask(ii,jj,k) - ave=ave+field_in(ii,jj,k)*weight + weight = mask(ii,jj,k) + ave = ave+field_in(ii,jj,k)*weight enddo field_out(i,j,k) = ave !Masked Sum (total_weight=1) - enddo; enddo; enddo - elseif(method == MPM) then - do k= ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d + enddo ; enddo ; enddo + elseif (method == MPM) then + do k=ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d i0 = isv_o+dl*(i-isv_d) j0 = jsv_o+dl*(j-jsv_d) ave = 0.0 @@ -4077,21 +4072,21 @@ subroutine downsample_field_3d(field_in, field_out, dl, method, mask, diag_cs, d do ii=i0,i0+dl-1 weight = mask(ii,jj,k) * diag_cs%G%dxCv(ii,jj) * diag_cs%h(ii,jj,k) total_weight = total_weight + weight - ave=ave+field_in(ii,jj,k)*weight + ave = ave+field_in(ii,jj,k)*weight enddo field_out(i,j,k) = ave/(total_weight+eps_face) !Avoid zero mask at all aggregating cells where ave=0.0 - enddo; enddo; enddo - elseif(method == MSK) then !The input field is a mask, subsample + enddo ; enddo ; enddo + elseif (method == MSK) then !The input field is a mask, subsample field_out(:,:,:) = 0.0 - do k= ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d + do k=ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d i0 = isv_o+dl*(i-isv_d) j0 = jsv_o+dl*(j-jsv_d) ave = 0.0 do jj=j0,j0+dl-1 ; do ii=i0,i0+dl-1 - ave=ave+field_in(ii,jj,k) - enddo; enddo - if(ave > 0.0) field_out(i,j,k)=1.0 - enddo; enddo; enddo + ave = ave+field_in(ii,jj,k) + enddo ; enddo + if (ave > 0.0) field_out(i,j,k)=1.0 + enddo ; enddo ; enddo else write (mesg,*) " unknown sampling method: ",method call MOM_error(FATAL, "downsample_field_3d: "//trim(mesg)//" "//trim(diag%debug_str)) @@ -4154,10 +4149,10 @@ subroutine downsample_field_2d(field_in, field_out, dl, method, mask, diag_cs, d weight = mask(ii,jj)*diag_cs%G%areaT(ii,jj) total_weight = total_weight + weight ave = ave+field_in(ii,jj)*weight - enddo; enddo + enddo ; enddo field_out(i,j) = ave/(total_weight + eps_area) !Avoid zero mask at all aggregating cells where ave=0.0 - enddo; enddo - elseif(method == SSP) then ! e.g., T_dfxy_cont_tendency_2d + enddo ; enddo + elseif (method == SSP) then ! e.g., T_dfxy_cont_tendency_2d do j=jsv_d,jev_d ; do i=isv_d,iev_d i0 = isv_o+dl*(i-isv_d) j0 = jsv_o+dl*(j-jsv_d) @@ -4167,11 +4162,11 @@ subroutine downsample_field_2d(field_in, field_out, dl, method, mask, diag_cs, d ! do ii=i0,i0+dl-1 ; do jj=j0,j0+dl-1 weight = mask(ii,jj) total_weight = total_weight + weight - ave=ave+field_in(ii,jj)*weight - enddo; enddo + ave = ave+field_in(ii,jj)*weight + enddo ; enddo field_out(i,j) = ave/(total_weight+epsilon) !Avoid zero mask at all aggregating cells where ave=0.0 - enddo; enddo - elseif(method == PSP) then ! e.g., umo_2d + enddo ; enddo + elseif (method == PSP) then ! e.g., umo_2d do j=jsv_d,jev_d ; do i=isv_d,iev_d i0 = isv_o+dl*(i-isv_d) j0 = jsv_o+dl*(j-jsv_d) @@ -4181,11 +4176,11 @@ subroutine downsample_field_2d(field_in, field_out, dl, method, mask, diag_cs, d do jj=j0,j0+dl-1 weight = mask(ii,jj) total_weight = total_weight +weight - ave=ave+field_in(ii,jj)*weight + ave = ave+field_in(ii,jj)*weight enddo field_out(i,j) = ave/(total_weight+epsilon) !Avoid zero mask at all aggregating cells where ave=0.0 - enddo; enddo - elseif(method == SPP) then ! e.g., vmo_2d + enddo ; enddo + elseif (method == SPP) then ! e.g., vmo_2d do j=jsv_d,jev_d ; do i=isv_d,iev_d i0 = isv_o+dl*(i-isv_d) j0 = jsv_o+dl*(j-jsv_d) @@ -4195,11 +4190,11 @@ subroutine downsample_field_2d(field_in, field_out, dl, method, mask, diag_cs, d do ii=i0,i0+dl-1 weight = mask(ii,jj) total_weight = total_weight +weight - ave=ave+field_in(ii,jj)*weight + ave = ave+field_in(ii,jj)*weight enddo field_out(i,j) = ave/(total_weight+epsilon) !Avoid zero mask at all aggregating cells where ave=0.0 - enddo; enddo - elseif(method == PMP) then + enddo ; enddo + elseif (method == PMP) then do j=jsv_d,jev_d ; do i=isv_d,iev_d i0 = isv_o+dl*(i-isv_d) j0 = jsv_o+dl*(j-jsv_d) @@ -4209,11 +4204,11 @@ subroutine downsample_field_2d(field_in, field_out, dl, method, mask, diag_cs, d do jj=j0,j0+dl-1 weight = mask(ii,jj) * diag_cs%G%dyCu(ii,jj)!*diag_cs%h(ii,jj,1) !Niki? total_weight = total_weight +weight - ave=ave+field_in(ii,jj)*weight + ave = ave+field_in(ii,jj)*weight enddo field_out(i,j) = ave/(total_weight+eps_len) !Avoid zero mask at all aggregating cells where ave=0.0 - enddo; enddo - elseif(method == MPP) then + enddo ; enddo + elseif (method == MPP) then do j=jsv_d,jev_d ; do i=isv_d,iev_d i0 = isv_o+dl*(i-isv_d) j0 = jsv_o+dl*(j-jsv_d) @@ -4223,21 +4218,21 @@ subroutine downsample_field_2d(field_in, field_out, dl, method, mask, diag_cs, d do ii=i0,i0+dl-1 weight = mask(ii,jj)* diag_cs%G%dxCv(ii,jj)!*diag_cs%h(ii,jj,1) !Niki? total_weight = total_weight +weight - ave=ave+field_in(ii,jj)*weight + ave = ave+field_in(ii,jj)*weight enddo field_out(i,j) = ave/(total_weight+eps_len) !Avoid zero mask at all aggregating cells where ave=0.0 - enddo; enddo - elseif(method == MSK) then !The input field is a mask, subsample + enddo ; enddo + elseif (method == MSK) then !The input field is a mask, subsample field_out(:,:) = 0.0 do j=jsv_d,jev_d ; do i=isv_d,iev_d i0 = isv_o+dl*(i-isv_d) j0 = jsv_o+dl*(j-jsv_d) ave = 0.0 do jj=j0,j0+dl-1 ; do ii=i0,i0+dl-1 - ave=ave+field_in(ii,jj) - enddo; enddo - if(ave > 0.0) field_out(i,j)=1.0 - enddo; enddo + ave = ave+field_in(ii,jj) + enddo ; enddo + if (ave > 0.0) field_out(i,j)=1.0 + enddo ; enddo else write (mesg,*) " unknown sampling method: ",method call MOM_error(FATAL, "downsample_field_2d: "//trim(mesg)//" "//trim(diag%debug_str)) @@ -4276,8 +4271,8 @@ subroutine downsample_mask_2d(field_in, field_out, dl, isc_o, jsc_o, isc_d, iec_ do jj=j0,j0+dl-1 ; do ii=i0,i0+dl-1 tot_non_zero = tot_non_zero + field_in(ii,jj) enddo;enddo - if(tot_non_zero > 0.0) field_out(i,j)=1.0 - enddo; enddo + if (tot_non_zero > 0.0) field_out(i,j)=1.0 + enddo ; enddo end subroutine downsample_mask_2d !> Allocate and compute the 3d down sampled mask @@ -4305,15 +4300,15 @@ subroutine downsample_mask_3d(field_in, field_out, dl, isc_o, jsc_o, isc_d, iec_ ks = lbound(field_in,3) ; ke = ubound(field_in,3) allocate(field_out(isd_d:ied_d,jsd_d:jed_d,ks:ke)) field_out(:,:,:) = 0.0 - do k= ks,ke ; do j=jsc_d,jec_d ; do i=isc_d,iec_d + do k=ks,ke ; do j=jsc_d,jec_d ; do i=isc_d,iec_d i0 = isc_o+dl*(i-isc_d) j0 = jsc_o+dl*(j-jsc_d) tot_non_zero = 0.0 do jj=j0,j0+dl-1 ; do ii=i0,i0+dl-1 tot_non_zero = tot_non_zero + field_in(ii,jj,k) enddo;enddo - if(tot_non_zero > 0.0) field_out(i,j,k)=1.0 - enddo; enddo; enddo + if (tot_non_zero > 0.0) field_out(i,j,k)=1.0 + enddo ; enddo ; enddo end subroutine downsample_mask_3d end module MOM_diag_mediator diff --git a/src/framework/MOM_domains.F90 b/src/framework/MOM_domains.F90 index 67989f37da..46cc9c526a 100644 --- a/src/framework/MOM_domains.F90 +++ b/src/framework/MOM_domains.F90 @@ -1309,7 +1309,7 @@ subroutine MOM_domains_init(MOM_dom, param_file, symmetric, static_memory, & !$ call fms_affinity_set('OCEAN', ocean_omp_hyper_thread, ocean_nthreads) !$ call omp_set_num_threads(ocean_nthreads) !$ write(6,*) "MOM_domains_mod OMPthreading ", fms_affinity_get(), omp_get_thread_num(), omp_get_num_threads() -!$ call flush(6) +!$ flush(6) !$ endif #endif call log_param(param_file, mdl, "!SYMMETRIC_MEMORY_", MOM_dom%symmetric, & diff --git a/src/framework/MOM_write_cputime.F90 b/src/framework/MOM_write_cputime.F90 index 2c1cb3378a..4ef00707fe 100644 --- a/src/framework/MOM_write_cputime.F90 +++ b/src/framework/MOM_write_cputime.F90 @@ -110,7 +110,7 @@ subroutine MOM_write_cputime_end(CS) ! Flush and close the output files. if (is_root_pe() .and. CS%fileCPU_ascii > 0) then - call flush(CS%fileCPU_ascii) + flush(CS%fileCPU_ascii) call close_file(CS%fileCPU_ascii) endif @@ -200,7 +200,7 @@ subroutine write_cputime(day, n, CS, nmax, call_end) reday, n, (CS%cputime2 / real(CLOCKS_PER_SEC)), & d_cputime / real(CLOCKS_PER_SEC) - call flush(CS%fileCPU_ascii) + flush(CS%fileCPU_ascii) endif CS%previous_calls = CS%previous_calls + 1 diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index 6ff0184d54..16d9294308 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -2052,7 +2052,6 @@ subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS) 'm2 s-1', conversion=US%L_to_m**2*US%s_to_T, & cmor_field_name='diftrblo', & cmor_long_name='Ocean Tracer Diffusivity due to Parameterized Mesoscale Advection', & - cmor_units='m2 s-1', & cmor_standard_name='ocean_tracer_diffusivity_due_to_parameterized_mesoscale_advection') CS%id_KH_u1 = register_diag_field('ocean_model', 'KHTH_u1', diag%axesCu1, Time, & diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index edd35d61f2..7ef0877321 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -162,8 +162,6 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, param_file, CS, data_h, nz_ #include "version_variable.h" character(len=40) :: mdl = "MOM_sponge" ! This module's name. logical :: use_sponge - real, allocatable, dimension(:,:,:) :: data_hu !< thickness at u points [H ~> m or kg m-2] - real, allocatable, dimension(:,:,:) :: data_hv !< thickness at v points [H ~> m or kg m-2] real, allocatable, dimension(:,:) :: Iresttime_u !< inverse of the restoring time at u points [T-1 ~> s-1] real, allocatable, dimension(:,:) :: Iresttime_v !< inverse of the restoring time at v points [T-1 ~> s-1] logical :: bndExtrapolation = .true. ! If true, extrapolate boundaries @@ -259,43 +257,40 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, param_file, CS, data_h, nz_ "The total number of columns where sponges are applied at h points.", like_default=.true.) if (CS%sponge_uv) then - - allocate(data_hu(G%isdB:G%iedB,G%jsd:G%jed,nz_data)) ; data_hu(:,:,:) = 0.0 - allocate(data_hv(G%isd:G%ied,G%jsdB:G%jedB,nz_data)) ; data_hv(:,:,:) = 0.0 allocate(Iresttime_u(G%isdB:G%iedB,G%jsd:G%jed)) ; Iresttime_u(:,:) = 0.0 allocate(Iresttime_v(G%isd:G%ied,G%jsdB:G%jedB)) ; Iresttime_v(:,:) = 0.0 ! u points CS%num_col_u = 0 ; !CS%fldno_u = 0 do j=CS%jsc,CS%jec ; do I=CS%iscB,CS%iecB - data_hu(I,j,:) = 0.5 * (data_h(i,j,:) + data_h(i+1,j,:)) - Iresttime_u(I,j) = 0.5 * (Iresttime(i,j) + Iresttime(i+1,j)) - if ((Iresttime_u(I,j)>0.0) .and. (G%mask2dCu(I,j)>0)) & - CS%num_col_u = CS%num_col_u + 1 + Iresttime_u(I,j) = 0.5 * (Iresttime(i,j) + Iresttime(i+1,j)) + if ((Iresttime_u(I,j)>0.0) .and. (G%mask2dCu(I,j)>0)) CS%num_col_u = CS%num_col_u + 1 enddo ; enddo if (CS%num_col_u > 0) then - allocate(CS%Iresttime_col_u(CS%num_col_u)) ; CS%Iresttime_col_u(:) = 0.0 - allocate(CS%col_i_u(CS%num_col_u)) ; CS%col_i_u(:) = 0 - allocate(CS%col_j_u(CS%num_col_u)) ; CS%col_j_u(:) = 0 - - ! pass indices, restoring time to the CS structure - col = 1 - do j=CS%jsc,CS%jec ; do I=CS%iscB,CS%iecB - if ((Iresttime_u(I,j)>0.0) .and. (G%mask2dCu(I,j)>0)) then - CS%col_i_u(col) = i ; CS%col_j_u(col) = j - CS%Iresttime_col_u(col) = Iresttime_u(i,j) - col = col + 1 - endif - enddo ; enddo - - ! same for total number of arbritary layers and correspondent data - - allocate(CS%Ref_hu%p(CS%nz_data,CS%num_col_u)) - do col=1,CS%num_col_u ; do K=1,CS%nz_data - CS%Ref_hu%p(K,col) = data_hu(CS%col_i_u(col),CS%col_j_u(col),K) - enddo ; enddo + allocate(CS%Iresttime_col_u(CS%num_col_u)) ; CS%Iresttime_col_u(:) = 0.0 + allocate(CS%col_i_u(CS%num_col_u)) ; CS%col_i_u(:) = 0 + allocate(CS%col_j_u(CS%num_col_u)) ; CS%col_j_u(:) = 0 + + ! Store the column indices and restoring rates in the CS structure + col = 1 + do j=CS%jsc,CS%jec ; do I=CS%iscB,CS%iecB + if ((Iresttime_u(I,j)>0.0) .and. (G%mask2dCu(I,j)>0)) then + CS%col_i_u(col) = I ; CS%col_j_u(col) = j + CS%Iresttime_col_u(col) = Iresttime_u(I,j) + col = col + 1 + endif + enddo ; enddo + + ! same for total number of arbritary layers and correspondent data + allocate(CS%Ref_hu%p(CS%nz_data,CS%num_col_u)) + do col=1,CS%num_col_u + I = CS%col_i_u(col) ; j = CS%col_j_u(col) + do k=1,CS%nz_data + CS%Ref_hu%p(k,col) = 0.5 * (data_h(i,j,k) + data_h(i+1,j,k)) + enddo + enddo endif total_sponge_cols_u = CS%num_col_u call sum_across_PEs(total_sponge_cols_u) @@ -305,10 +300,8 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, param_file, CS, data_h, nz_ ! v points CS%num_col_v = 0 ; !CS%fldno_v = 0 do J=CS%jscB,CS%jecB; do i=CS%isc,CS%iec - data_hv(i,J,:) = 0.5 * (data_h(i,j,:) + data_h(i,j+1,:)) Iresttime_v(i,J) = 0.5 * (Iresttime(i,j) + Iresttime(i,j+1)) - if ((Iresttime_v(i,J)>0.0) .and. (G%mask2dCv(i,J)>0)) & - CS%num_col_v = CS%num_col_v + 1 + if ((Iresttime_v(i,J)>0.0) .and. (G%mask2dCv(i,J)>0)) CS%num_col_v = CS%num_col_v + 1 enddo ; enddo if (CS%num_col_v > 0) then @@ -329,9 +322,12 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, param_file, CS, data_h, nz_ ! same for total number of arbritary layers and correspondent data allocate(CS%Ref_hv%p(CS%nz_data,CS%num_col_v)) - do col=1,CS%num_col_v ; do K=1,CS%nz_data - CS%Ref_hv%p(K,col) = data_hv(CS%col_i_v(col),CS%col_j_v(col),K) - enddo ; enddo + do col=1,CS%num_col_v + i = CS%col_i_v(col) ; J = CS%col_j_v(col) + do k=1,CS%nz_data + CS%Ref_hv%p(k,col) = 0.5 * (data_h(i,j,k) + data_h(i,j+1,k)) + enddo + enddo endif total_sponge_cols_v = CS%num_col_v call sum_across_PEs(total_sponge_cols_v) @@ -473,7 +469,7 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, param_file, CS) if ((Iresttime(i,j)>0.0) .and. (G%mask2dT(i,j)>0)) then CS%col_i(col) = i ; CS%col_j(col) = j CS%Iresttime_col(col) = Iresttime(i,j) - col = col +1 + col = col + 1 endif enddo ; enddo endif @@ -505,7 +501,7 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, param_file, CS) if ((Iresttime_u(I,j)>0.0) .and. (G%mask2dCu(I,j)>0)) then CS%col_i_u(col) = i ; CS%col_j_u(col) = j CS%Iresttime_col_u(col) = Iresttime_u(i,j) - col = col +1 + col = col + 1 endif enddo ; enddo ! same for total number of arbritary layers and correspondent data @@ -531,7 +527,7 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, param_file, CS) if ((Iresttime_v(i,J)>0.0) .and. (G%mask2dCv(i,J)>0)) then CS%col_i_v(col) = i ; CS%col_j_v(col) = j CS%Iresttime_col_v(col) = Iresttime_v(i,j) - col = col +1 + col = col + 1 endif enddo ; enddo endif @@ -800,8 +796,7 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) real :: m_to_Z ! A unit conversion factor from m to Z. real, allocatable, dimension(:) :: tmp_val2 ! data values on the original grid real, dimension(SZK_(G)) :: tmp_val1 ! data values remapped to model grid - real :: hu(SZIB_(G), SZJ_(G), SZK_(G)) ! A temporary array for h at u pts - real :: hv(SZI_(G), SZJB_(G), SZK_(G)) ! A temporary array for h at v pts + real, dimension(SZK_(G)) :: h_col ! A column of thicknesses at h, u or v points [H ~> m or kg m-2] real, allocatable, dimension(:,:,:) :: sp_val ! A temporary array for fields real, allocatable, dimension(:,:,:) :: mask_z ! A temporary array for field mask at h pts real, dimension(:), allocatable :: hsrc ! Source thicknesses [Z ~> m]. @@ -830,47 +825,45 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) if (.not. present(Time)) & call MOM_error(FATAL,"apply_ALE_sponge: No time information provided") do m=1,CS%fldno - nz_data = CS%Ref_val(m)%nz_data - allocate(sp_val(G%isd:G%ied,G%jsd:G%jed,1:nz_data)) - allocate(mask_z(G%isd:G%ied,G%jsd:G%jed,1:nz_data)) - sp_val(:,:,:) = 0.0 - mask_z(:,:,:) = 0.0 - call horiz_interp_and_extrap_tracer(CS%Ref_val(m)%id, Time, 1.0, G, sp_val, mask_z, z_in, & + nz_data = CS%Ref_val(m)%nz_data + allocate(sp_val(G%isd:G%ied,G%jsd:G%jed,1:nz_data)) ; sp_val(:,:,:) = 0.0 + allocate(mask_z(G%isd:G%ied,G%jsd:G%jed,1:nz_data)) ; mask_z(:,:,:) = 0.0 + call horiz_interp_and_extrap_tracer(CS%Ref_val(m)%id, Time, 1.0, G, sp_val, mask_z, z_in, & z_edges_in, missing_value, .true., .false., .false., & spongeOnGrid=CS%SpongeDataOngrid, m_to_Z=US%m_to_Z, & answers_2018=CS%hor_regrid_answers_2018) - allocate( hsrc(nz_data) ) - allocate( tmpT1d(nz_data) ) - do c=1,CS%num_col - i = CS%col_i(c) ; j = CS%col_j(c) - CS%Ref_val(m)%p(1:nz_data,c) = sp_val(i,j,1:nz_data) - ! Build the source grid - zTopOfCell = 0. ; zBottomOfCell = 0. ; nPoints = 0; hsrc(:) = 0.0; tmpT1d(:) = -99.9 - do k=1,nz_data - if (mask_z(CS%col_i(c),CS%col_j(c),k) == 1.0) then - zBottomOfCell = -min( z_edges_in(k+1), G%bathyT(CS%col_i(c),CS%col_j(c)) ) - tmpT1d(k) = sp_val(CS%col_i(c),CS%col_j(c),k) - elseif (k>1) then - zBottomOfCell = -G%bathyT(CS%col_i(c),CS%col_j(c)) - tmpT1d(k) = tmpT1d(k-1) - else ! This next block should only ever be reached over land - tmpT1d(k) = -99.9 - endif - hsrc(k) = zTopOfCell - zBottomOfCell - if (hsrc(k)>0.) nPoints = nPoints + 1 - zTopOfCell = zBottomOfCell ! Bottom becomes top for next value of k - enddo - ! In case data is deeper than model - hsrc(nz_data) = hsrc(nz_data) + ( zTopOfCell + G%bathyT(CS%col_i(c),CS%col_j(c)) ) - CS%Ref_val(m)%h(1:nz_data,c) = GV%Z_to_H*hsrc(1:nz_data) - CS%Ref_val(m)%p(1:nz_data,c) = tmpT1d(1:nz_data) - do k=2,nz_data - ! if (mask_z(i,j,k)==0.) & - if (CS%Ref_val(m)%h(k,c) <= 0.001*GV%m_to_H) & - ! some confusion here about why the masks are not correct returning from horiz_interp - ! reverting to using a minimum thickness criteria - CS%Ref_val(m)%p(k,c) = CS%Ref_val(m)%p(k-1,c) - enddo + allocate( hsrc(nz_data) ) + allocate( tmpT1d(nz_data) ) + do c=1,CS%num_col + i = CS%col_i(c) ; j = CS%col_j(c) + do k=1,nz_data ; CS%Ref_val(m)%p(k,c) = sp_val(i,j,k) ; enddo + ! Build the source grid + zTopOfCell = 0. ; zBottomOfCell = 0. ; nPoints = 0 ; hsrc(:) = 0. ; tmpT1d(:) = -99.9 + do k=1,nz_data + if (mask_z(CS%col_i(c),CS%col_j(c),k) == 1.0) then + zBottomOfCell = -min( z_edges_in(k+1), G%bathyT(CS%col_i(c),CS%col_j(c)) ) + tmpT1d(k) = sp_val(CS%col_i(c),CS%col_j(c),k) + elseif (k>1) then + zBottomOfCell = -G%bathyT(CS%col_i(c),CS%col_j(c)) + tmpT1d(k) = tmpT1d(k-1) + else ! This next block should only ever be reached over land + tmpT1d(k) = -99.9 + endif + hsrc(k) = zTopOfCell - zBottomOfCell + if (hsrc(k)>0.) nPoints = nPoints + 1 + zTopOfCell = zBottomOfCell ! Bottom becomes top for next value of k + enddo + ! In case data is deeper than model + hsrc(nz_data) = hsrc(nz_data) + ( zTopOfCell + G%bathyT(CS%col_i(c),CS%col_j(c)) ) + CS%Ref_val(m)%h(1:nz_data,c) = GV%Z_to_H*hsrc(1:nz_data) + CS%Ref_val(m)%p(1:nz_data,c) = tmpT1d(1:nz_data) + do k=2,nz_data + ! if (mask_z(i,j,k)==0.) & + if (CS%Ref_val(m)%h(k,c) <= 0.001*GV%m_to_H) & + ! some confusion here about why the masks are not correct returning from horiz_interp + ! reverting to using a minimum thickness criteria + CS%Ref_val(m)%p(k,c) = CS%Ref_val(m)%p(k-1,c) + enddo enddo deallocate(sp_val, mask_z, hsrc, tmpT1d) enddo @@ -879,23 +872,22 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) endif allocate(tmp_val2(nz_data)) - do m=1,CS%fldno - do c=1,CS%num_col -! c is an index for the next 3 lines but a multiplier for the rest of the loop -! Therefore we use c as per C code and increment the index where necessary. - i = CS%col_i(c) ; j = CS%col_j(c) - damp = dt * CS%Iresttime_col(c) - I1pdamp = 1.0 / (1.0 + damp) + do c=1,CS%num_col + i = CS%col_i(c) ; j = CS%col_j(c) + damp = dt * CS%Iresttime_col(c) + I1pdamp = 1.0 / (1.0 + damp) + do k=1,nz ; h_col(k) = h(i,j,k) ; enddo + do m=1,CS%fldno tmp_val2(1:nz_data) = CS%Ref_val(m)%p(1:nz_data,c) if (CS%time_varying_sponges) then call remapping_core_h(CS%remap_cs, nz_data, CS%Ref_val(m)%h(1:nz_data,c), tmp_val2, & - CS%nz, h(i,j,:), tmp_val1, h_neglect, h_neglect_edge) + CS%nz, h_col, tmp_val1, h_neglect, h_neglect_edge) else - call remapping_core_h(CS%remap_cs,nz_data, CS%Ref_h%p(1:nz_data,c), tmp_val2, & - CS%nz, h(i,j,:), tmp_val1, h_neglect, h_neglect_edge) + call remapping_core_h(CS%remap_cs, nz_data, CS%Ref_h%p(1:nz_data,c), tmp_val2, & + CS%nz, h_col, tmp_val1, h_neglect, h_neglect_edge) endif - !Backward Euler method - CS%var(m)%p(i,j,1:CS%nz) = I1pdamp * (CS%var(m)%p(i,j,1:CS%nz) + tmp_val1 * damp) + ! Backward Euler method + do k=1,CS%nz ; CS%var(m)%p(i,j,k) = I1pdamp * (CS%var(m)%p(i,j,k) + tmp_val1(k) * damp) ; enddo enddo enddo @@ -907,10 +899,6 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) !enddo if (CS%sponge_uv) then - ! u points - do j=CS%jsc,CS%jec; do I=CS%iscB,CS%iecB; do k=1,nz - hu(I,j,k) = 0.5 * (h(i,j,k) + h(i+1,j,k)) - enddo ; enddo ; enddo if (CS%time_varying_sponges) then if (.not. present(Time)) & call MOM_error(FATAL,"apply_ALE_sponge: No time information provided") @@ -925,10 +913,8 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) ! call pass_var(sp_val,G%Domain) ! call pass_var(mask_z,G%Domain) do c=1,CS%num_col - ! c is an index for the next 3 lines but a multiplier for the rest of the loop - ! Therefore we use c as per C code and increment the index where necessary. i = CS%col_i(c) ; j = CS%col_j(c) - CS%Ref_val_u%p(1:nz_data,c) = sp_val(i,j,1:nz_data) + do k=1,nz_data ; CS%Ref_val_u%p(k,c) = sp_val(i,j,k) ; enddo enddo deallocate (sp_val, mask_z) @@ -945,10 +931,8 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) ! call pass_var(mask_z,G%Domain) do c=1,CS%num_col - ! c is an index for the next 3 lines but a multiplier for the rest of the loop - ! Therefore we use c as per C code and increment the index where necessary. i = CS%col_i(c) ; j = CS%col_j(c) - CS%Ref_val_v%p(1:nz_data,c) = sp_val(i,j,1:nz_data) + do k=1,nz_data ; CS%Ref_val_v%p(k,c) = sp_val(i,j,k) ; enddo enddo deallocate (sp_val, mask_z) @@ -957,42 +941,42 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) nz_data = CS%nz_data endif + ! u points do c=1,CS%num_col_u - i = CS%col_i_u(c) ; j = CS%col_j_u(c) + I = CS%col_i_u(c) ; j = CS%col_j_u(c) damp = dt * CS%Iresttime_col_u(c) I1pdamp = 1.0 / (1.0 + damp) + do k=1,nz ; h_col(k) = 0.5 * (h(i,j,k) + h(i+1,j,k)) ; enddo if (CS%time_varying_sponges) nz_data = CS%Ref_val(m)%nz_data tmp_val2(1:nz_data) = CS%Ref_val_u%p(1:nz_data,c) if (CS%time_varying_sponges) then call remapping_core_h(CS%remap_cs, nz_data, CS%Ref_val_u%h(:,c), tmp_val2, & - CS%nz, hu(i,j,:), tmp_val1, h_neglect, h_neglect_edge) + CS%nz, h_col, tmp_val1, h_neglect, h_neglect_edge) else call remapping_core_h(CS%remap_cs, nz_data, CS%Ref_hu%p(:,c), tmp_val2, & - CS%nz, hu(i,j,:), tmp_val1, h_neglect, h_neglect_edge) + CS%nz, h_col, tmp_val1, h_neglect, h_neglect_edge) endif !Backward Euler method - CS%var_u%p(i,j,:) = I1pdamp * (CS%var_u%p(i,j,:) + tmp_val1 * damp) + do k=1,CS%nz ; CS%var_u%p(I,j,k) = I1pdamp * (CS%var_u%p(I,j,k) + tmp_val1(k) * damp) ; enddo enddo ! v points - do J=CS%jscB,CS%jecB; do i=CS%isc,CS%iec; do k=1,nz - hv(i,J,k) = 0.5 * (h(i,j,k) + h(i,j+1,k)) - enddo ; enddo ; enddo do c=1,CS%num_col_v i = CS%col_i_v(c) ; j = CS%col_j_v(c) damp = dt * CS%Iresttime_col_v(c) I1pdamp = 1.0 / (1.0 + damp) + do k=1,nz ; h_col(k) = 0.5 * (h(i,j,k) + h(i,j+1,k)) ; enddo tmp_val2(1:nz_data) = CS%Ref_val_v%p(1:nz_data,c) if (CS%time_varying_sponges) then call remapping_core_h(CS%remap_cs, CS%nz_data, CS%Ref_val_v%h(:,c), tmp_val2, & - CS%nz, hv(i,j,:), tmp_val1, h_neglect, h_neglect_edge) + CS%nz, h_col, tmp_val1, h_neglect, h_neglect_edge) else call remapping_core_h(CS%remap_cs, CS%nz_data, CS%Ref_hv%p(:,c), tmp_val2, & - CS%nz, hv(i,j,:), tmp_val1, h_neglect, h_neglect_edge) + CS%nz, h_col, tmp_val1, h_neglect, h_neglect_edge) endif !Backward Euler method - CS%var_v%p(i,j,:) = I1pdamp * (CS%var_v%p(i,j,:) + tmp_val1 * damp) + do k=1,CS%nz ; CS%var_u%p(i,J,k) = I1pdamp * (CS%var_u%p(i,J,k) + tmp_val1(k) * damp) ; enddo enddo endif diff --git a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 index 079655f787..3bea0d9937 100644 --- a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 +++ b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 @@ -482,11 +482,10 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, C if (CS%ML_resort) then if (id_clock_resort>0) call cpu_clock_begin(id_clock_resort) if (CS%ML_presort_nz_conv_adj > 0) & - call convective_adjustment(h(:,1:), u, v, R0(:,1:), Rcv(:,1:), T(:,1:), & - S(:,1:), eps, d_eb, dKE_CA, cTKE, j, G, GV, US, CS, & - CS%ML_presort_nz_conv_adj) + call convective_adjustment(h, u, v, R0, Rcv, T, S, eps, d_eb, dKE_CA, cTKE, j, G, GV, & + US, CS, CS%ML_presort_nz_conv_adj) - call sort_ML(h(:,1:), R0(:,1:), eps, G, GV, CS, ksort) + call sort_ML(h, R0, eps, G, GV, CS, ksort) if (id_clock_resort>0) call cpu_clock_end(id_clock_resort) else do k=1,nz ; do i=is,ie ; ksort(i,k) = k ; enddo ; enddo @@ -495,8 +494,7 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, C ! Undergo instantaneous entrainment into the buffer layers and mixed layers ! to remove hydrostatic instabilities. Any water that is lighter than ! currently in the mixed or buffer layer is entrained. - call convective_adjustment(h(:,1:), u, v, R0(:,1:), Rcv(:,1:), T(:,1:), & - S(:,1:), eps, d_eb, dKE_CA, cTKE, j, G, GV, US, CS) + call convective_adjustment(h, u, v, R0, Rcv, T, S, eps, d_eb, dKE_CA, cTKE, j, G, GV, US, CS) do i=is,ie ; h_CA(i) = h(i,1) ; enddo if (id_clock_adjustment>0) call cpu_clock_end(id_clock_adjustment) @@ -535,18 +533,15 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, C ! Pen_SW_bnd = components to penetrative shortwave radiation call extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, & CS%H_limit_fluxes, CS%use_river_heat_content, CS%use_calving_heat_content, & - h(:,1:), T(:,1:), netMassInOut, netMassOut, Net_heat, Net_salt, Pen_SW_bnd,& + h(:,1:), T(:,1:), netMassInOut, netMassOut, Net_heat, Net_salt, Pen_SW_bnd, & tv, aggregate_FW_forcing) ! This subroutine causes the mixed layer to entrain to depth of free convection. - call mixedlayer_convection(h(:,1:), d_eb, htot, Ttot, Stot, uhtot, vhtot, & - R0_tot, Rcv_tot, u, v, T(:,1:), S(:,1:), & - R0(:,1:), Rcv(:,1:), eps, & - dR0_dT, dRcv_dT, dR0_dS, dRcv_dS, & - netMassInOut, netMassOut, Net_heat, Net_salt, & - nsw, Pen_SW_bnd, opacity_band, Conv_En, & - dKE_FC, j, ksort, G, GV, US, CS, tv, fluxes, dt, & - aggregate_FW_forcing) + call mixedlayer_convection(h, d_eb, htot, Ttot, Stot, uhtot, vhtot, R0_tot, Rcv_tot, & + u, v, T, S, R0, Rcv, eps, dR0_dT, dRcv_dT, dR0_dS, dRcv_dS, & + netMassInOut, netMassOut, Net_heat, Net_salt, & + nsw, Pen_SW_bnd, opacity_band, Conv_En, dKE_FC, & + j, ksort, G, GV, US, CS, tv, fluxes, dt, aggregate_FW_forcing) if (id_clock_conv>0) call cpu_clock_end(id_clock_conv) @@ -563,8 +558,8 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, C j, ksort, G, GV, US, CS) ! Here the mechanically driven entrainment occurs. - call mechanical_entrainment(h(:,1:), d_eb, htot, Ttot, Stot, uhtot, vhtot, & - R0_tot, Rcv_tot, u, v, T(:,1:), S(:,1:), R0(:,1:), Rcv(:,1:), eps, dR0_dT, dRcv_dT, & + call mechanical_entrainment(h, d_eb, htot, Ttot, Stot, uhtot, vhtot, & + R0_tot, Rcv_tot, u, v, T, S, R0, Rcv, eps, dR0_dT, dRcv_dT, & cMKE, Idt_diag, nsw, Pen_SW_bnd, opacity_band, TKE, & Idecay_len_TKE, j, ksort, G, GV, US, CS) @@ -803,43 +798,39 @@ end subroutine bulkmixedlayer !! is lighter than currently in the mixed- or buffer- layer is entrained. subroutine convective_adjustment(h, u, v, R0, Rcv, T, S, eps, d_eb, & dKE_CA, cTKE, j, G, GV, US, CS, nz_conv) - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. - real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]. + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]. !! The units of h are referred to as H below. - real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: u !< Zonal velocities interpolated to h + real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: u !< Zonal velocities interpolated to h !! points [L T-1 ~> m s-1]. - real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: v !< Zonal velocities interpolated to h + real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: v !< Zonal velocities interpolated to h !! points [L T-1 ~> m s-1]. - real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: T !< Layer temperatures [degC]. - real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: S !< Layer salinities [ppt]. - real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: R0 !< Potential density referenced to + real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: R0 !< Potential density referenced to !! surface pressure [R ~> kg m-3]. - real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: Rcv !< The coordinate defining potential + real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: Rcv !< The coordinate defining potential !! density [R ~> kg m-3]. - real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: d_eb !< The downward increase across a layer + real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: T !< Layer temperatures [degC]. + real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: S !< Layer salinities [ppt]. + real, dimension(SZI_(G),SZK_(GV)), intent(in) :: eps !< The negligibly small amount of water + !! that will be left in each layer [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: d_eb !< The downward increase across a layer !! in the entrainment from below [H ~> m or kg m-2]. !! Positive values go with mass gain by !! a layer. - real, dimension(SZI_(G),SZK_(GV)), intent(in) :: eps !< The negligibly small amount of water - !! that will be left in each layer [H ~> m or kg m-2]. - real, dimension(SZI_(G),SZK_(GV)), intent(out) :: dKE_CA !< The vertically integrated change in + real, dimension(SZI_(G),SZK_(GV)), intent(out) :: dKE_CA !< The vertically integrated change in !! kinetic energy due to convective !! adjustment [Z L2 T-2 ~> m3 s-2]. - real, dimension(SZI_(G),SZK_(GV)), intent(out) :: cTKE !< The buoyant turbulent kinetic energy + real, dimension(SZI_(G),SZK_(GV)), intent(out) :: cTKE !< The buoyant turbulent kinetic energy !! source due to convective adjustment !! [Z L2 T-2 ~> m3 s-2]. - integer, intent(in) :: j !< The j-index to work on. - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - type(bulkmixedlayer_CS), pointer :: CS !< The control structure for this module. - integer, optional, intent(in) :: nz_conv !< If present, the number of layers + integer, intent(in) :: j !< The j-index to work on. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(bulkmixedlayer_CS), pointer :: CS !< The control structure for this module. + integer, optional, intent(in) :: nz_conv !< If present, the number of layers !! over which to do convective adjustment !! (perhaps CS%nkml). -! This subroutine does instantaneous convective entrainment into the buffer -! layers and mixed layers to remove hydrostatic instabilities. Any water that -! is lighter than currently in the mixed- or buffer- layer is entrained. - ! Local variables real, dimension(SZI_(G)) :: & htot, & ! The total depth of the layers being considered for @@ -942,7 +933,7 @@ subroutine mixedlayer_convection(h, d_eb, htot, Ttot, Stot, uhtot, vhtot, & aggregate_FW_forcing) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. - real, dimension(SZI_(G),SZK_(GV)), & + real, dimension(SZI_(G),SZK0_(GV)), & intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]. !! The units of h are referred to as H below. real, dimension(SZI_(G),SZK_(GV)), & @@ -966,14 +957,14 @@ subroutine mixedlayer_convection(h, d_eb, htot, Ttot, Stot, uhtot, vhtot, & intent(in) :: u !< Zonal velocities interpolated to h points [L T-1 ~> m s-1]. real, dimension(SZI_(G),SZK_(GV)), & intent(in) :: v !< Zonal velocities interpolated to h points [L T-1 ~> m s-1]. - real, dimension(SZI_(G),SZK_(GV)), & + real, dimension(SZI_(G),SZK0_(GV)), & intent(in) :: T !< Layer temperatures [degC]. - real, dimension(SZI_(G),SZK_(GV)), & + real, dimension(SZI_(G),SZK0_(GV)), & intent(in) :: S !< Layer salinities [ppt]. - real, dimension(SZI_(G),SZK_(GV)), & + real, dimension(SZI_(G),SZK0_(GV)), & intent(in) :: R0 !< Potential density referenced to !! surface pressure [R ~> kg m-3]. - real, dimension(SZI_(G),SZK_(GV)), & + real, dimension(SZI_(G),SZK0_(GV)), & intent(in) :: Rcv !< The coordinate defining potential !! density [R ~> kg m-3]. real, dimension(SZI_(G),SZK_(GV)), & @@ -1503,7 +1494,7 @@ subroutine mechanical_entrainment(h, d_eb, htot, Ttot, Stot, uhtot, vhtot, & type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - real, dimension(SZI_(G),SZK_(GV)), & + real, dimension(SZI_(G),SZK0_(GV)), & intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]. real, dimension(SZI_(G),SZK_(GV)), & intent(inout) :: d_eb !< The downward increase across a layer in the @@ -1526,14 +1517,14 @@ subroutine mechanical_entrainment(h, d_eb, htot, Ttot, Stot, uhtot, vhtot, & intent(in) :: u !< Zonal velocities interpolated to h points [L T-1 ~> m s-1]. real, dimension(SZI_(G),SZK_(GV)), & intent(in) :: v !< Zonal velocities interpolated to h points [L T-1 ~> m s-1]. - real, dimension(SZI_(G),SZK_(GV)), & + real, dimension(SZI_(G),SZK0_(GV)), & intent(in) :: T !< Layer temperatures [degC]. - real, dimension(SZI_(G),SZK_(GV)), & + real, dimension(SZI_(G),SZK0_(GV)), & intent(in) :: S !< Layer salinities [ppt]. - real, dimension(SZI_(G),SZK_(GV)), & + real, dimension(SZI_(G),SZK0_(GV)), & intent(in) :: R0 !< Potential density referenced to !! surface pressure [R ~> kg m-3]. - real, dimension(SZI_(G),SZK_(GV)), & + real, dimension(SZI_(G),SZK0_(GV)), & intent(in) :: Rcv !< The coordinate defining potential !! density [R ~> kg m-3]. real, dimension(SZI_(G),SZK_(GV)), & @@ -1840,8 +1831,8 @@ end subroutine mechanical_entrainment subroutine sort_ML(h, R0, eps, G, GV, CS, ksort) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. - real, dimension(SZI_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]. - real, dimension(SZI_(G),SZK_(GV)), intent(in) :: R0 !< The potential density used to sort + real, dimension(SZI_(G),SZK0_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZK0_(GV)), intent(in) :: R0 !< The potential density used to sort !! the layers [R ~> kg m-3]. real, dimension(SZI_(G),SZK_(GV)), intent(in) :: eps !< The (small) thickness that must !! remain in each layer [H ~> m or kg m-2]. diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index fbc6bb2d16..1683e21fbe 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -30,7 +30,7 @@ module MOM_diabatic_aux #include public diabatic_aux_init, diabatic_aux_end -public make_frazil, adjust_salt, differential_diffuse_T_S, triDiagTS +public make_frazil, adjust_salt, differential_diffuse_T_S, triDiagTS, triDiagTS_Eulerian public find_uv_at_h, diagnoseMLDbyDensityDifference, applyBoundaryFluxesInOut, set_pen_shortwave public diagnoseMLDbyEnergy @@ -395,9 +395,10 @@ subroutine triDiagTS(G, GV, is, ie, js, je, hold, ea, eb, T, S) real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: S !< Layer salinities [ppt]. ! Local variables - real :: b1(SZIB_(G)), d1(SZIB_(G)) ! b1, c1, and d1 are variables used by the - real :: c1(SZIB_(G),SZK_(G)) ! tridiagonal solver. - real :: h_tr, b_denom_1 + real :: b1(SZIB_(G)) ! A variable used by the tridiagonal solver [H-1 ~> m-2 or m2 kg-1]. + real :: d1(SZIB_(G)) ! A variable used by the tridiagonal solver [nondim]. + real :: c1(SZIB_(G),SZK_(G)) ! A variable used by the tridiagonal solver [nondim]. + real :: h_tr, b_denom_1 ! Two temporary thicknesses [H ~> m or kg m-2]. integer :: i, j, k !$OMP parallel do default(shared) private(h_tr,b1,d1,c1,b_denom_1) @@ -425,9 +426,58 @@ subroutine triDiagTS(G, GV, is, ie, js, je, hold, ea, eb, T, S) enddo end subroutine triDiagTS +!> This is a simple tri-diagonal solver for T and S, with mixing across interfaces but no net +!! transfer of mass. +subroutine triDiagTS_Eulerian(G, GV, is, ie, js, je, hold, ent, T, S) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + integer, intent(in) :: is !< The start i-index to work on. + integer, intent(in) :: ie !< The end i-index to work on. + integer, intent(in) :: js !< The start j-index to work on. + integer, intent(in) :: je !< The end j-index to work on. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: hold !< The layer thicknesses before entrainment, + !! [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(in) :: ent !< The amount of fluid mixed across an interface + !! within this time step [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: T !< Layer potential temperatures [degC]. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: S !< Layer salinities [ppt]. + + ! Local variables + real :: b1(SZIB_(G)) ! A variable used by the tridiagonal solver [H-1 ~> m-2 or m2 kg-1]. + real :: d1(SZIB_(G)) ! A variable used by the tridiagonal solver [nondim]. + real :: c1(SZIB_(G),SZK_(G)) ! A variable used by the tridiagonal solver [nondim]. + real :: h_tr, b_denom_1 ! Two temporary thicknesses [H ~> m or kg m-2]. + integer :: i, j, k + + !$OMP parallel do default(shared) private(h_tr,b1,d1,c1,b_denom_1) + do j=js,je + do i=is,ie + h_tr = hold(i,j,1) + GV%H_subroundoff + b1(i) = 1.0 / (h_tr + ent(i,j,2)) + d1(i) = h_tr * b1(i) + T(i,j,1) = (b1(i)*h_tr)*T(i,j,1) + S(i,j,1) = (b1(i)*h_tr)*S(i,j,1) + enddo + do k=2,G%ke ; do i=is,ie + c1(i,k) = ent(i,j,K) * b1(i) + h_tr = hold(i,j,k) + GV%H_subroundoff + b_denom_1 = h_tr + d1(i)*ent(i,j,K) + b1(i) = 1.0 / (b_denom_1 + ent(i,j,K+1)) + d1(i) = b_denom_1 * b1(i) + T(i,j,k) = b1(i) * (h_tr*T(i,j,k) + ent(i,j,K)*T(i,j,k-1)) + S(i,j,k) = b1(i) * (h_tr*S(i,j,k) + ent(i,j,K)*S(i,j,k-1)) + enddo ; enddo + do k=G%ke-1,1,-1 ; do i=is,ie + T(i,j,k) = T(i,j,k) + c1(i,k+1)*T(i,j,k+1) + S(i,j,k) = S(i,j,k) + c1(i,k+1)*S(i,j,k+1) + enddo ; enddo + enddo +end subroutine triDiagTS_Eulerian + + !> This subroutine calculates u_h and v_h (velocities at thickness !! points), optionally using the entrainment amounts passed in as arguments. -subroutine find_uv_at_h(u, v, h, u_h, v_h, G, GV, US, ea, eb) +subroutine find_uv_at_h(u, v, h, u_h, v_h, G, GV, US, ea, eb, zero_mix) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -449,6 +499,9 @@ subroutine find_uv_at_h(u, v, h, u_h, v_h, G, GV, US, ea, eb) optional, intent(in) :: eb !< The amount of fluid entrained from the layer !! below within this time step [H ~> m or kg m-2]. !! Omitting eb is the same as setting it to 0. + logical, optional, intent(in) :: zero_mix !< If true, do the calculation of u_h and + !! v_h as though ea and eb were being supplied with + !! uniformly zero values. ! local variables real :: b_denom_1 ! The first term in the denominator of b1 [H ~> m or kg m-2]. @@ -459,7 +512,7 @@ subroutine find_uv_at_h(u, v, h, u_h, v_h, G, GV, US, ea, eb) real :: a_e(SZI_(G)), a_w(SZI_(G)) ! velocity points, ~1/2 in the open ! ocean, nondimensional. real :: sum_area, Idenom - logical :: mix_vertically + logical :: mix_vertically, zero_mixing integer :: i, j, k, is, ie, js, je, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke call cpu_clock_begin(id_clock_uv_at_h) @@ -469,9 +522,11 @@ subroutine find_uv_at_h(u, v, h, u_h, v_h, G, GV, US, ea, eb) if (present(ea) .neqv. present(eb)) call MOM_error(FATAL, & "find_uv_at_h: Either both ea and eb or neither one must be present "// & "in call to find_uv_at_h.") -!$OMP parallel do default(none) shared(is,ie,js,je,G,GV,mix_vertically,h,h_neglect, & -!$OMP eb,u_h,u,v_h,v,nz,ea) & -!$OMP private(sum_area,Idenom,a_w,a_e,a_s,a_n,b_denom_1,b1,d1,c1) + zero_mixing = .false. ; if (present(zero_mix)) zero_mixing = zero_mix + if (zero_mixing) mix_vertically = .false. + !$OMP parallel do default(none) shared(is,ie,js,je,G,GV,mix_vertically,zero_mixing,h, & + !$OMP h_neglect,ea,eb,u_h,u,v_h,v,nz) & + !$OMP private(sum_area,Idenom,a_w,a_e,a_s,a_n,b_denom_1,b1,d1,c1) do j=js,je do i=is,ie sum_area = G%areaCu(I-1,j) + G%areaCu(I,j) @@ -515,6 +570,17 @@ subroutine find_uv_at_h(u, v, h, u_h, v_h, G, GV, US, ea, eb) u_h(i,j,k) = u_h(i,j,k) + c1(i,k+1)*u_h(i,j,k+1) v_h(i,j,k) = v_h(i,j,k) + c1(i,k+1)*v_h(i,j,k+1) enddo ; enddo + elseif (zero_mixing) then + do i=is,ie + b1(i) = 1.0 / (h(i,j,1) + h_neglect) + u_h(i,j,1) = (h(i,j,1)*b1(i)) * (a_e(i)*u(I,j,1) + a_w(i)*u(I-1,j,1)) + v_h(i,j,1) = (h(i,j,1)*b1(i)) * (a_n(i)*v(i,J,1) + a_s(i)*v(i,J-1,1)) + enddo + do k=2,nz ; do i=is,ie + b1(i) = 1.0 / (h(i,j,k) + h_neglect) + u_h(i,j,k) = (h(i,j,k) * (a_e(i)*u(I,j,k) + a_w(i)*u(I-1,j,k))) * b1(i) + v_h(i,j,k) = (h(i,j,k) * (a_n(i)*v(i,J,k) + a_s(i)*v(i,J-1,k))) * b1(i) + enddo ; enddo else do k=1,nz ; do i=is,ie u_h(i,j,k) = a_e(i)*u(I,j,k) + a_w(i)*u(I-1,j,k) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index e5fe6724c5..2fb6a27542 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -13,9 +13,8 @@ module MOM_diabatic_driver use MOM_CVMix_ddiff, only : CVMix_ddiff_is_used use MOM_diabatic_aux, only : diabatic_aux_init, diabatic_aux_end, diabatic_aux_CS use MOM_diabatic_aux, only : make_frazil, adjust_salt, differential_diffuse_T_S, triDiagTS -use MOM_diabatic_aux, only : find_uv_at_h, diagnoseMLDbyDensityDifference, applyBoundaryFluxesInOut -use MOM_diabatic_aux, only : diagnoseMLDbyEnergy -use MOM_diabatic_aux, only : set_pen_shortwave +use MOM_diabatic_aux, only : triDiagTS_Eulerian, find_uv_at_h, diagnoseMLDbyDensityDifference +use MOM_diabatic_aux, only : applyBoundaryFluxesInOut, diagnoseMLDbyEnergy, set_pen_shortwave use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr use MOM_diag_mediator, only : diag_ctrl, time_type, diag_update_remap_grids use MOM_diag_mediator, only : diag_ctrl, query_averaging_enabled, enable_averages, disable_averaging @@ -62,7 +61,7 @@ module MOM_diabatic_driver use MOM_ALE_sponge, only : apply_ALE_sponge, ALE_sponge_CS use MOM_time_manager, only : time_type, real_to_time, operator(-), operator(<=) use MOM_tracer_flow_control, only : call_tracer_column_fns, tracer_flow_control_CS -use MOM_tracer_diabatic, only : tracer_vertdiff +use MOM_tracer_diabatic, only : tracer_vertdiff, tracer_vertdiff_Eulerian use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : thermo_var_ptrs, vertvisc_type, accel_diag_ptrs use MOM_variables, only : cont_diag_ptrs, MOM_thermovar_chksum, p3d @@ -174,14 +173,17 @@ module MOM_diabatic_driver !>@{ Diagnostic IDs integer :: id_cg1 = -1 ! diag handle for mode-1 speed integer, allocatable, dimension(:) :: id_cn ! diag handle for all mode speeds - integer :: id_wd = -1, id_ea = -1, id_eb = -1 ! used by layer diabatic - integer :: id_dudt_dia = -1, id_dvdt_dia = -1, id_ea_s = -1, id_eb_s = -1 - integer :: id_ea_t = -1, id_eb_t = -1 - integer :: id_Kd_heat = -1, id_Kd_salt = -1, id_Kd_interface = -1, id_Kd_ePBL = -1 - integer :: id_Tdif = -1, id_Tadv = -1, id_Sdif = -1, id_Sadv = -1 - integer :: id_MLD_003 = -1, id_MLD_0125 = -1, id_MLD_user = -1, id_mlotstsq = -1 - integer :: id_MLD_EN1 = -1, id_MLD_EN2 = -1, id_MLD_EN3= -1 - integer :: id_subMLN2 = -1 + integer :: id_ea = -1, id_eb = -1 ! used by layer diabatic + integer :: id_ea_t = -1, id_eb_t = -1, id_ea_s = -1, id_eb_s = -1 + integer :: id_Kd_heat = -1, id_Kd_salt = -1, id_Kd_int = -1, id_Kd_ePBL = -1 + integer :: id_Tdif = -1, id_Sdif = -1, id_Tadv = -1, id_Sadv = -1 + ! These are handles to diagnostics related to the mixed layer properties. + integer :: id_MLD_003 = -1, id_MLD_0125 = -1, id_MLD_user = -1, id_mlotstsq = -1 + integer :: id_MLD_EN1 = -1, id_MLD_EN2 = -1, id_MLD_EN3 = -1, id_subMLN2 = -1 + + ! These are handles to diatgnostics that are only available in non-ALE layered mode. + integer :: id_wd = -1 + integer :: id_dudt_dia = -1, id_dvdt_dia = -1 integer :: id_hf_dudt_dia_2d = -1, id_hf_dvdt_dia_2d = -1 ! diagnostic for fields prior to applying diapycnal physics @@ -502,16 +504,10 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim ! Kd_int [Z2 T-1 ~> m2 s-1]. Kd_extra_S , & ! The extra diffusivity of salinity due to double diffusion relative to ! Kd_int [Z2 T-1 ~> m2 s-1]. - zeros_h, & ! An array of zeros to handle diagnostics that should be removed. Kd_ePBL, & ! test array of diapycnal diffusivities at interfaces [Z2 T-1 ~> m2 s-1] Tdif_flx, & ! diffusive diapycnal heat flux across interfaces [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1] Sdif_flx ! diffusive diapycnal salt flux across interfaces [ppt H T-1 ~> ppt m s-1 or ppt kg m-2 s-1] - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: & - zeros_u ! An array of zeros for u-point diagnostics that should be removed. - real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: & - zeros_v ! An array of zeros for v-point diagnostics that should be removed. - logical :: in_boundary(SZI_(G)) ! True if there are no massive layers below, ! where massive is defined as sufficiently thick that ! the no-flux boundary conditions have not restricted @@ -573,10 +569,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim if (CS%use_kappa_shear .or. CS%use_CVMix_shear) then if (CS%use_geothermal) then - ! The presence of ea_s causes find_uv_at_h to use a tridiagonal solver, - ! which changes answers at the level of roundoff because ((A*B / A) /= B). - ent_s(:,:,:) = 0.0 - call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US, ent_s(:,:,1), ent_s(:,:,1)) + call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US, zero_mix=.true.) else call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US) endif @@ -902,10 +895,10 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim ent_t(i,j,K) = ent_s(i,j,K) ; ent_t(i,j,K+1) = ent_s(i,j,K+1) enddo ; enddo ; enddo if (CS%tracer_tridiag) then - call tracer_vertdiff(hold, ent_t(:,:,1:nz), ent_t(:,:,2:nz+1), dt, tv%T, G, GV) - call tracer_vertdiff(hold, ent_s(:,:,1:nz), ent_s(:,:,2:nz+1), dt, tv%S, G, GV) + call tracer_vertdiff_Eulerian(hold, ent_t, dt, tv%T, G, GV) + call tracer_vertdiff_Eulerian(hold, ent_s, dt, tv%S, G, GV) else - call triDiagTS(G, GV, is, ie, js, je, hold, ent_s(:,:,1:nz), ent_s(:,:,2:nz+1), tv%T, tv%S) + call triDiagTS_Eulerian(G, GV, is, ie, js, je, hold, ent_s, tv%T, tv%S) endif ! diagnose temperature, salinity, heat, and salt tendencies @@ -935,17 +928,17 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim call diag_update_remap_grids(CS%diag) ! Diagnose the diapycnal diffusivities and other related quantities. - if (CS%id_Kd_interface > 0) call post_data(CS%id_Kd_interface, Kd_int, CS%diag) - if (CS%id_Kd_heat > 0) call post_data(CS%id_Kd_heat, Kd_heat, CS%diag) - if (CS%id_Kd_salt > 0) call post_data(CS%id_Kd_salt, Kd_salt, CS%diag) - if (CS%id_Kd_ePBL > 0) call post_data(CS%id_Kd_ePBL, Kd_ePBL, CS%diag) - - if (CS%id_ea > 0) call post_data(CS%id_ea, ent_s(:,:,1:nz), CS%diag) - if (CS%id_eb > 0) call post_data(CS%id_eb, ent_s(:,:,2:nz+1), CS%diag) - if (CS%id_ea_t > 0) call post_data(CS%id_ea_t, ent_t(:,:,1:nz), CS%diag) - if (CS%id_eb_t > 0) call post_data(CS%id_eb_t, ent_t(:,:,2:nz+1), CS%diag) - if (CS%id_ea_s > 0) call post_data(CS%id_ea_s, ent_s(:,:,1:nz), CS%diag) - if (CS%id_eb_s > 0) call post_data(CS%id_eb_s, ent_s(:,:,2:nz+1), CS%diag) + if (CS%id_Kd_int > 0) call post_data(CS%id_Kd_int, Kd_int, CS%diag) + if (CS%id_Kd_heat > 0) call post_data(CS%id_Kd_heat, Kd_heat, CS%diag) + if (CS%id_Kd_salt > 0) call post_data(CS%id_Kd_salt, Kd_salt, CS%diag) + if (CS%id_Kd_ePBL > 0) call post_data(CS%id_Kd_ePBL, Kd_ePBL, CS%diag) + + if (CS%id_ea > 0) call post_data(CS%id_ea, ent_s(:,:,1:nz), CS%diag) + if (CS%id_eb > 0) call post_data(CS%id_eb, ent_s(:,:,2:nz+1), CS%diag) + if (CS%id_ea_t > 0) call post_data(CS%id_ea_t, ent_t(:,:,1:nz), CS%diag) + if (CS%id_eb_t > 0) call post_data(CS%id_eb_t, ent_t(:,:,2:nz+1), CS%diag) + if (CS%id_ea_s > 0) call post_data(CS%id_ea_s, ent_s(:,:,1:nz), CS%diag) + if (CS%id_eb_s > 0) call post_data(CS%id_eb_s, ent_s(:,:,2:nz+1), CS%diag) Idt = 1.0 / dt if (CS%id_Tdif > 0) then do j=js,je ; do i=is,ie @@ -1043,18 +1036,6 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim endif endif ! CS%use_sponge - !! Diagnostics that are always zero and should be removed from this routine. - if ((CS%id_wd > 0) .or. (CS%id_Tadv > 0) .or. (CS%id_Sadv > 0)) zeros_h(:,:,:) = 0.0 - if ((CS%id_dudt_dia > 0) .or. (CS%id_hf_dudt_dia_2d > 0)) zeros_u(:,:,:) = 0.0 - if ((CS%id_dvdt_dia > 0) .or. (CS%id_hf_dvdt_dia_2d > 0)) zeros_v(:,:,:) = 0.0 - if (CS%id_wd > 0) call post_data(CS%id_wd, zeros_h, CS%diag) - if (CS%id_Tadv > 0) call post_data(CS%id_Tadv, zeros_h, CS%diag) - if (CS%id_Sadv > 0) call post_data(CS%id_Sadv, zeros_h, CS%diag) - if (CS%id_dudt_dia > 0) call post_data(CS%id_dudt_dia, zeros_u, CS%diag) - if (CS%id_dvdt_dia > 0) call post_data(CS%id_dvdt_dia, zeros_v, CS%diag) - if (CS%id_hf_dudt_dia_2d > 0) call post_data(CS%id_hf_dudt_dia_2d, zeros_u(:,:,1), CS%diag) - if (CS%id_hf_dvdt_dia_2d > 0) call post_data(CS%id_hf_dvdt_dia_2d, zeros_v(:,:,1), CS%diag) - call disable_averaging(CS%diag) if (showCallTree) call callTree_leave("diabatic_ALE_legacy()") @@ -1114,13 +1095,8 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, Kd_extra_S , & ! The extra diffusivity of salinity due to double diffusion relative to ! Kd_int returned from set_diffusivity [Z2 T-1 ~> m2 s-1]. Kd_ePBL, & ! boundary layer or convective diapycnal diffusivities at interfaces [Z2 T-1 ~> m2 s-1] - zeros_h, & ! An array of zeros for h-point diagnostics that should be removed. Tdif_flx, & ! diffusive diapycnal heat flux across interfaces [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1] Sdif_flx ! diffusive diapycnal salt flux across interfaces [ppt H T-1 ~> ppt m s-1 or ppt kg m-2 s-1] - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: & - zeros_u ! An array of zeros for u-point diagnostics that should be removed. - real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: & - zeros_v ! An array of zeros for v-point diagnostics that should be removed. logical :: in_boundary(SZI_(G)) ! True if there are no massive layers below, ! where massive is defined as sufficiently thick that @@ -1183,10 +1159,7 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, if (CS%use_kappa_shear .or. CS%use_CVMix_shear) then if (CS%use_geothermal) then - ! The presence of ea_s causes find_uv_at_h to use a tridiagonal solver, - ! which changes answers at the level of roundoff because ((A*B / A) /= B). - ! Note that ent_s(:,:,:) = 0.0 - call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US, ent_s(:,:,1:nz), ent_s(:,:,1:nz)) + call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US, zero_mix=.true.) else call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US) endif @@ -1216,7 +1189,7 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, endif ! Store the diagnosed typical diffusivity at interfaces. - if (CS%id_Kd_interface > 0) call post_data(CS%id_Kd_interface, Kd_heat, CS%diag) + if (CS%id_Kd_int > 0) call post_data(CS%id_Kd_int, Kd_heat, CS%diag) ! Set diffusivities for heat and salt separately, and possibly change the meaning of Kd_heat. if (CS%double_diffuse) then @@ -1441,8 +1414,8 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, "Kd_salt (diabatic_ALE)") ! Changes T and S via the tridiagonal solver; no change to h - call tracer_vertdiff(h, ent_t(:,:,1:nz), ent_t(:,:,2:nz+1), dt, tv%T, G, GV) - call tracer_vertdiff(h, ent_s(:,:,1:nz), ent_s(:,:,2:nz+1), dt, tv%S, G, GV) + call tracer_vertdiff_Eulerian(h, ent_t, dt, tv%T, G, GV) + call tracer_vertdiff_Eulerian(h, ent_s, dt, tv%S, G, GV) ! In ALE-mode, layer thicknesses do not change. Therefore, we can use h below if (CS%diabatic_diff_tendency_diag) then @@ -1555,17 +1528,6 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, call pass_var(visc%Kv_slow, G%Domain, To_All+Omit_Corners, halo=1) call cpu_clock_end(id_clock_pass) - !! Diagnostics that are always zero and should be removed from this routine. - if ((CS%id_dudt_dia > 0) .or. (CS%id_hf_dudt_dia_2d > 0)) zeros_u(:,:,:) = 0.0 - if ((CS%id_dvdt_dia > 0) .or. (CS%id_hf_dvdt_dia_2d > 0)) zeros_v(:,:,:) = 0.0 - if ((CS%id_Tadv > 0) .or. (CS%id_Sadv > 0)) zeros_h(:,:,:) = 0.0 - if (CS%id_dudt_dia > 0) call post_data(CS%id_dudt_dia, zeros_u, CS%diag) - if (CS%id_dvdt_dia > 0) call post_data(CS%id_dvdt_dia, zeros_v, CS%diag) - if (CS%id_hf_dudt_dia_2d > 0) call post_data(CS%id_hf_dudt_dia_2d, zeros_u(:,:,1), CS%diag) - if (CS%id_hf_dvdt_dia_2d > 0) call post_data(CS%id_hf_dvdt_dia_2d, zeros_v(:,:,1), CS%diag) - if (CS%id_Tadv > 0) call post_data(CS%id_Tadv, zeros_h, CS%diag) - if (CS%id_Sadv > 0) call post_data(CS%id_Sadv, zeros_h, CS%diag) - call disable_averaging(CS%diag) if (showCallTree) call callTree_leave("diabatic_ALE()") @@ -1634,7 +1596,6 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e ! Kd_int [Z2 T-1 ~> m2 s-1]. Kd_extra_S , & ! The extra diffusivity of salinity due to double diffusion relative to ! Kd_int [Z2 T-1 ~> m2 s-1]. - zeros_h, & ! An array of zeros to handle diagnostics that should be removed. Tdif_flx, & ! diffusive diapycnal heat flux across interfaces [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1] Tadv_flx, & ! advective diapycnal heat flux across interfaces [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1] Sdif_flx, & ! diffusive diapycnal salt flux across interfaces [ppt H T-1 ~> ppt m s-1 or ppt kg m-2 s-1] @@ -2486,14 +2447,12 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e endif ! Diagnose the diapycnal diffusivities and other related quantities. - if (CS%id_Kd_interface > 0) call post_data(CS%id_Kd_interface, Kd_int, CS%diag) - if (CS%id_Kd_heat > 0) call post_data(CS%id_Kd_heat, Kd_heat, CS%diag) - if (CS%id_Kd_salt > 0) call post_data(CS%id_Kd_salt, Kd_salt, CS%diag) - if (CS%id_Kd_ePBL > 0) zeros_h(:,:,:) = 0.0 - if (CS%id_Kd_ePBL > 0) call post_data(CS%id_Kd_ePBL, zeros_h, CS%diag) + if (CS%id_Kd_int > 0) call post_data(CS%id_Kd_int, Kd_int, CS%diag) + if (CS%id_Kd_heat > 0) call post_data(CS%id_Kd_heat, Kd_heat, CS%diag) + if (CS%id_Kd_salt > 0) call post_data(CS%id_Kd_salt, Kd_salt, CS%diag) - if (CS%id_ea > 0) call post_data(CS%id_ea, ea, CS%diag) - if (CS%id_eb > 0) call post_data(CS%id_eb, eb, CS%diag) + if (CS%id_ea > 0) call post_data(CS%id_ea, ea, CS%diag) + if (CS%id_eb > 0) call post_data(CS%id_eb, eb, CS%diag) if (CS%id_dudt_dia > 0) call post_data(CS%id_dudt_dia, ADp%du_dt_dia, CS%diag) if (CS%id_dvdt_dia > 0) call post_data(CS%id_dvdt_dia, ADp%dv_dt_dia, CS%diag) @@ -3058,32 +3017,31 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di 'Layer entrainment from above per timestep', 'm', conversion=GV%H_to_m) CS%id_eb = register_diag_field('ocean_model', 'eb', diag%axesTL, Time, & 'Layer entrainment from below per timestep', 'm', conversion=GV%H_to_m) - !### if (.not.CS%useALEalgorithm) then - CS%id_wd = register_diag_field('ocean_model', 'wd', diag%axesTi, Time, & - 'Diapycnal velocity', 'm s-1', conversion=GV%H_to_m) - if (CS%id_wd > 0) call safe_alloc_ptr(CDp%diapyc_vel,isd,ied,jsd,jed,nz+1) - - CS%id_dudt_dia = register_diag_field('ocean_model', 'dudt_dia', diag%axesCuL, Time, & - 'Zonal Acceleration from Diapycnal Mixing', 'm s-2', conversion=US%L_T2_to_m_s2) - CS%id_dvdt_dia = register_diag_field('ocean_model', 'dvdt_dia', diag%axesCvL, Time, & - 'Meridional Acceleration from Diapycnal Mixing', 'm s-2', conversion=US%L_T2_to_m_s2) - - CS%id_hf_dudt_dia_2d = register_diag_field('ocean_model', 'hf_dudt_dia_2d', diag%axesCu1, Time, & - 'Depth-sum Fractional Thickness-weighted Zonal Acceleration from Diapycnal Mixing', & - 'm s-2', conversion=US%L_T2_to_m_s2) - if (CS%id_hf_dudt_dia_2d > 0) then - call safe_alloc_ptr(ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz) - call safe_alloc_ptr(ADp%du_dt_dia,IsdB,IedB,jsd,jed,nz) - endif + if (.not.CS%useALEalgorithm) then + CS%id_wd = register_diag_field('ocean_model', 'wd', diag%axesTi, Time, & + 'Diapycnal velocity', 'm s-1', conversion=GV%H_to_m) + if (CS%id_wd > 0) call safe_alloc_ptr(CDp%diapyc_vel,isd,ied,jsd,jed,nz+1) - CS%id_hf_dvdt_dia_2d = register_diag_field('ocean_model', 'hf_dvdt_dia_2d', diag%axesCv1, Time, & - 'Depth-sum Fractional Thickness-weighted Meridional Acceleration from Diapycnal Mixing', & - 'm s-2', conversion=US%L_T2_to_m_s2) - if (CS%id_hf_dvdt_dia_2d > 0) then - call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz) - call safe_alloc_ptr(ADp%dv_dt_dia,isd,ied,JsdB,JedB,nz) + CS%id_dudt_dia = register_diag_field('ocean_model', 'dudt_dia', diag%axesCuL, Time, & + 'Zonal Acceleration from Diapycnal Mixing', 'm s-2', conversion=US%L_T2_to_m_s2) + CS%id_dvdt_dia = register_diag_field('ocean_model', 'dvdt_dia', diag%axesCvL, Time, & + 'Meridional Acceleration from Diapycnal Mixing', 'm s-2', conversion=US%L_T2_to_m_s2) + + CS%id_hf_dudt_dia_2d = register_diag_field('ocean_model', 'hf_dudt_dia_2d', diag%axesCu1, Time, & + 'Depth-sum Fractional Thickness-weighted Zonal Acceleration from Diapycnal Mixing', & + 'm s-2', conversion=US%L_T2_to_m_s2) + if (CS%id_hf_dudt_dia_2d > 0) call safe_alloc_ptr(ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz) + + CS%id_hf_dvdt_dia_2d = register_diag_field('ocean_model', 'hf_dvdt_dia_2d', diag%axesCv1, Time, & + 'Depth-sum Fractional Thickness-weighted Meridional Acceleration from Diapycnal Mixing', & + 'm s-2', conversion=US%L_T2_to_m_s2) + if (CS%id_hf_dvdt_dia_2d > 0) call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz) + + if ((CS%id_dudt_dia > 0) .or. (CS%id_hf_dudt_dia_2d > 0)) & + call safe_alloc_ptr(ADp%du_dt_dia,IsdB,IedB,jsd,jed,nz) + if ((CS%id_dvdt_dia > 0) .or. (CS%id_hf_dvdt_dia_2d > 0)) & + call safe_alloc_ptr(ADp%dv_dt_dia,isd,ied,JsdB,JedB,nz) endif - !### endif if (CS%use_int_tides) then CS%id_cg1 = register_diag_field('ocean_model', 'cn1', diag%axesT1, & @@ -3102,15 +3060,19 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di CS%id_Tdif = register_diag_field('ocean_model',"Tflx_dia_diff", diag%axesTi, & Time, "Diffusive diapycnal temperature flux across interfaces", & "degC m s-1", conversion=GV%H_to_m*US%s_to_T) - CS%id_Tadv = register_diag_field('ocean_model',"Tflx_dia_adv", diag%axesTi, & - Time, "Advective diapycnal temperature flux across interfaces", & - "degC m s-1", conversion=GV%H_to_m*US%s_to_T) + if (.not.CS%useALEalgorithm) then + CS%id_Tadv = register_diag_field('ocean_model',"Tflx_dia_adv", diag%axesTi, & + Time, "Advective diapycnal temperature flux across interfaces", & + "degC m s-1", conversion=GV%H_to_m*US%s_to_T) + endif CS%id_Sdif = register_diag_field('ocean_model',"Sflx_dia_diff", diag%axesTi, & Time, "Diffusive diapycnal salnity flux across interfaces", & "psu m s-1", conversion=GV%H_to_m*US%s_to_T) - CS%id_Sadv = register_diag_field('ocean_model',"Sflx_dia_adv", diag%axesTi, & - Time, "Advective diapycnal salnity flux across interfaces", & - "psu m s-1", conversion=GV%H_to_m*US%s_to_T) + if (.not.CS%useALEalgorithm) then + CS%id_Sadv = register_diag_field('ocean_model',"Sflx_dia_adv", diag%axesTi, & + Time, "Advective diapycnal salnity flux across interfaces", & + "psu m s-1", conversion=GV%H_to_m*US%s_to_T) + endif CS%id_MLD_003 = register_diag_field('ocean_model', 'MLD_003', diag%axesT1, Time, & 'Mixed layer depth (delta rho = 0.03)', 'm', conversion=US%Z_to_m, & cmor_field_name='mlotst', cmor_long_name='Ocean Mixed Layer Thickness Defined by Sigma T', & @@ -3158,9 +3120,6 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di "stratification at the base of the mixed layer.", & units='m', default=50.0, scale=US%m_to_Z) - if (CS%id_dudt_dia > 0) call safe_alloc_ptr(ADp%du_dt_dia,IsdB,IedB,jsd,jed,nz) - if (CS%id_dvdt_dia > 0) call safe_alloc_ptr(ADp%dv_dt_dia,isd,ied,JsdB,JedB,nz) - ! diagnostics for values prior to diabatic and prior to ALE CS%id_u_predia = register_diag_field('ocean_model', 'u_predia', diag%axesCuL, Time, & 'Zonal velocity before diabatic forcing', 'm s-1', conversion=US%L_T_to_m_s) @@ -3180,7 +3139,7 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di !call set_diffusivity_init(Time, G, param_file, diag, CS%set_diff_CSp, CS%int_tide_CSp) - CS%id_Kd_interface = register_diag_field('ocean_model', 'Kd_interface', diag%axesTi, Time, & + CS%id_Kd_int = register_diag_field('ocean_model', 'Kd_interface', diag%axesTi, Time, & 'Total diapycnal diffusivity at interfaces', 'm2 s-1', conversion=US%Z2_T_to_m2_s) if (CS%use_energetic_PBL) then CS%id_Kd_ePBL = register_diag_field('ocean_model', 'Kd_ePBL', diag%axesTi, Time, & diff --git a/src/parameterizations/vertical/MOM_geothermal.F90 b/src/parameterizations/vertical/MOM_geothermal.F90 index 2ecaa4a78e..dadb5ac526 100644 --- a/src/parameterizations/vertical/MOM_geothermal.F90 +++ b/src/parameterizations/vertical/MOM_geothermal.F90 @@ -581,7 +581,7 @@ subroutine geothermal_init(Time, G, GV, US, param_file, diag, CS, useALEalgorith ! post the static geothermal heating field id = register_static_field('ocean_model', 'geo_heat', diag%axesT1, & 'Geothermal heat flux into ocean', 'W m-2', conversion=US%QRZ_T_to_W_m2, & - cmor_field_name='hfgeou', cmor_units='W m-2', & + cmor_field_name='hfgeou', & cmor_standard_name='upward_geothermal_heat_flux_at_sea_floor', & cmor_long_name='Upward geothermal heat flux at sea floor', & x_cell_method='mean', y_cell_method='mean', area_cell_method='mean') diff --git a/src/parameterizations/vertical/MOM_sponge.F90 b/src/parameterizations/vertical/MOM_sponge.F90 index 4566abcef7..dcd0ac4e02 100644 --- a/src/parameterizations/vertical/MOM_sponge.F90 +++ b/src/parameterizations/vertical/MOM_sponge.F90 @@ -477,8 +477,6 @@ subroutine apply_sponge(h, dt, G, GV, US, ea, eb, CS, Rcv_ml) endif do c=1,CS%num_col -! c is an index for the next 3 lines but a multiplier for the rest of the loop -! Therefore we use c as per C code and increment the index where necessary. i = CS%col_i(c) ; j = CS%col_j(c) damp = dt * CS%Iresttime_col(c) diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index 7786bf5b46..81f8b29a63 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -822,9 +822,7 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC) enddo if (do_any_shelf) then if (CS%harmonic_visc) then - call find_coupling_coef(a_shelf, hvel, do_i_shelf, h_harm, bbl_thick, & - kv_bbl, z_i, h_ml, dt, j, G, GV, US, CS, & - visc, forces, work_on_u=.true., OBC=OBC, shelf=.true.) + do k=1,nz ; do I=Isq,Ieq ; hvel_shelf(I,k) = hvel(I,k) ; enddo ; enddo else ! Find upwind-biased thickness near the surface. ! Perhaps this needs to be done more carefully, via find_eta. do I=Isq,Ieq ; if (do_i_shelf(I)) then @@ -850,10 +848,10 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC) endif endif ; enddo enddo - call find_coupling_coef(a_shelf, hvel_shelf, do_i_shelf, h_harm, & - bbl_thick, kv_bbl, z_i, h_ml, dt, j, G, GV, US, CS, & - visc, forces, work_on_u=.true., OBC=OBC, shelf=.true.) endif + call find_coupling_coef(a_shelf, hvel_shelf, do_i_shelf, h_harm, bbl_thick, & + kv_bbl, z_i, h_ml, dt, j, G, GV, US, CS, visc, forces, & + work_on_u=.true., OBC=OBC, shelf=.true.) do I=Isq,Ieq ; if (do_i_shelf(I)) CS%a1_shelf_u(I,j) = a_shelf(I,1) ; enddo endif endif @@ -979,7 +977,7 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC) call find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, h_ml, & dt, j, G, GV, US, CS, visc, forces, work_on_u=.false., OBC=OBC) if ( allocated(hML_v)) then - do i=is,ie ; if (do_i(i)) then ; hML_v(i,J) = h_ml(i) ; endif ; enddo + do i=is,ie ; if (do_i(i)) then ; hML_v(i,J) = h_ml(i) ; endif ; enddo endif do_any_shelf = .false. if (associated(forces%frac_shelf_v)) then @@ -990,9 +988,7 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC) enddo if (do_any_shelf) then if (CS%harmonic_visc) then - call find_coupling_coef(a_shelf, hvel, do_i_shelf, h_harm, bbl_thick, & - kv_bbl, z_i, h_ml, dt, j, G, GV, US, CS, visc, & - forces, work_on_u=.false., OBC=OBC, shelf=.true.) + do k=1,nz ; do i=is,ie ; hvel_shelf(i,k) = hvel(i,k) ; enddo ; enddo else ! Find upwind-biased thickness near the surface. ! Perhaps this needs to be done more carefully, via find_eta. do i=is,ie ; if (do_i_shelf(i)) then @@ -1018,10 +1014,10 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC) endif endif ; enddo enddo - call find_coupling_coef(a_shelf, hvel_shelf, do_i_shelf, h_harm, & - bbl_thick, kv_bbl, z_i, h_ml, dt, j, G, GV, US, CS, & - visc, forces, work_on_u=.false., OBC=OBC, shelf=.true.) endif + call find_coupling_coef(a_shelf, hvel_shelf, do_i_shelf, h_harm, bbl_thick, & + kv_bbl, z_i, h_ml, dt, j, G, GV, US, CS, visc, forces, & + work_on_u=.false., OBC=OBC, shelf=.true.) do i=is,ie ; if (do_i_shelf(i)) CS%a1_shelf_v(i,J) = a_shelf(i,1) ; enddo endif endif diff --git a/src/tracer/MOM_tracer_diabatic.F90 b/src/tracer/MOM_tracer_diabatic.F90 index ec7c025db0..6b9a12f696 100644 --- a/src/tracer/MOM_tracer_diabatic.F90 +++ b/src/tracer/MOM_tracer_diabatic.F90 @@ -13,7 +13,7 @@ module MOM_tracer_diabatic implicit none ; private #include -public tracer_vertdiff +public tracer_vertdiff, tracer_vertdiff_Eulerian public applyTracerBoundaryFluxesInOut contains @@ -41,6 +41,7 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, & real, dimension(SZI_(G),SZJ_(G)), optional,intent(in) :: btm_flux !< The (negative upward) bottom flux of the !! tracer in [CU kg m-2 T-1 ~> CU kg m-2 s-1] or !! [CU H ~> CU m or CU kg m-2] if + !! convert_flux_in is .false. real, dimension(SZI_(G),SZJ_(G)), optional,intent(inout) :: btm_reservoir !< amount of tracer in a bottom reservoir !! [CU kg m-2]; formerly [CU m] real, optional,intent(in) :: sink_rate !< rate at which the tracer sinks @@ -70,7 +71,6 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, & !! in roundoff and can be neglected [H ~> m or kg m-2]. logical :: convert_flux = .true. - integer :: i, j, k, is, ie, js, je, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke @@ -84,20 +84,17 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, & h_neglect = GV%H_subroundoff sink_dist = 0.0 if (present(sink_rate)) sink_dist = (dt*sink_rate) * GV%m_to_H -!$OMP parallel default(none) shared(is,ie,js,je,sfc_src,btm_src,sfc_flux,dt,G,GV,btm_flux, & -!$OMP sink_rate,btm_reservoir,nz,sink_dist,ea, & -!$OMP h_old,convert_flux,h_neglect,eb,tr) & -!$OMP private(sink,h_minus_dsink,b_denom_1,b1,d1,h_tr,c1) -!$OMP do + !$OMP parallel default(shared) private(sink,h_minus_dsink,b_denom_1,b1,d1,h_tr,c1) + !$OMP do do j=js,je; do i=is,ie ; sfc_src(i,j) = 0.0 ; btm_src(i,j) = 0.0 ; enddo ; enddo if (present(sfc_flux)) then if (convert_flux) then -!$OMP do + !$OMP do do j = js, je; do i = is,ie sfc_src(i,j) = (sfc_flux(i,j)*dt) * GV%kg_m2_to_H enddo ; enddo else -!$OMP do + !$OMP do do j = js, je; do i = is,ie sfc_src(i,j) = sfc_flux(i,j) enddo ; enddo @@ -105,12 +102,12 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, & endif if (present(btm_flux)) then if (convert_flux) then -!$OMP do + !$OMP do do j = js, je; do i = is,ie btm_src(i,j) = (btm_flux(i,j)*dt) * GV%kg_m2_to_H enddo ; enddo else -!$OMP do + !$OMP do do j = js, je; do i = is,ie btm_src(i,j) = btm_flux(i,j) enddo ; enddo @@ -118,7 +115,7 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, & endif if (present(sink_rate)) then -!$OMP do + !$OMP do do j=js,je ! Find the sinking rates at all interfaces, limiting them if necesary ! so that the characteristics do not cross within a timestep. @@ -186,7 +183,7 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, & endif ; enddo ; enddo enddo else -!$OMP do + !$OMP do do j=js,je do i=is,ie ; if (G%mask2dT(i,j) > -0.5) then h_tr = h_old(i,j,1) + h_neglect @@ -217,11 +214,211 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, & endif ; enddo ; enddo enddo endif - -!$OMP end parallel + !$OMP end parallel end subroutine tracer_vertdiff + +!> This subroutine solves a tridiagonal equation for the final tracer concentrations after +!! Eulerian mixing, and possibly sinking or surface and bottom sources, are applied. The sinking +!! is implemented with an fully implicit upwind advection scheme. Alternate time units can be +!! used for the timestep, surface and bottom fluxes and sink_rate provided they are all consistent. +subroutine tracer_vertdiff_Eulerian(h_old, ent, dt, tr, G, GV, & + sfc_flux, btm_flux, btm_reservoir, sink_rate, convert_flux_in) + type(ocean_grid_type), intent(in) :: G !< ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_old !< layer thickness before entrainment + !! [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: ent !< Amount of fluid mixed across interfaces + !! [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: tr !< tracer concentration in concentration units [CU] + real, intent(in) :: dt !< amount of time covered by this call [T ~> s] + real, dimension(SZI_(G),SZJ_(G)), optional,intent(in) :: sfc_flux !< surface flux of the tracer in units of + !! [CU kg m-2 T-1 ~> CU kg m-2 s-1] or + !! [CU H ~> CU m or CU kg m-2] if + !! convert_flux_in is .false. + real, dimension(SZI_(G),SZJ_(G)), optional,intent(in) :: btm_flux !< The (negative upward) bottom flux of the + !! tracer in [CU kg m-2 T-1 ~> CU kg m-2 s-1] or + !! [CU H ~> CU m or CU kg m-2] if + !! convert_flux_in is .false. + real, dimension(SZI_(G),SZJ_(G)), optional,intent(inout) :: btm_reservoir !< amount of tracer in a bottom reservoir + !! [CU kg m-2]; formerly [CU m] + real, optional,intent(in) :: sink_rate !< rate at which the tracer sinks + !! [m T-1 ~> m s-1] + logical, optional,intent(in) :: convert_flux_in !< True if the specified sfc_flux needs + !! to be integrated in time + + ! local variables + real :: sink_dist !< The distance the tracer sinks in a time step [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G)) :: & + sfc_src, & !< The time-integrated surface source of the tracer [CU H ~> CU m or CU kg m-2]. + btm_src !< The time-integrated bottom source of the tracer [CU H ~> CU m or CU kg m-2]. + real, dimension(SZI_(G)) :: & + b1, & !< b1 is used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1]. + d1 !! d1=1-c1 is used by the tridiagonal solver, nondimensional. + real :: c1(SZI_(G),SZK_(GV)) !< c1 is used by the tridiagonal solver [nondim]. + real :: h_minus_dsink(SZI_(G),SZK_(GV)) !< The layer thickness minus the + !! difference in sinking rates across the layer [H ~> m or kg m-2]. + !! By construction, 0 <= h_minus_dsink < h_work. + real :: sink(SZI_(G),SZK_(GV)+1) !< The tracer's sinking distances at the + !! interfaces, limited to prevent characteristics from + !! crossing within a single timestep [H ~> m or kg m-2]. + real :: b_denom_1 !< The first term in the denominator of b1 [H ~> m or kg m-2]. + real :: h_tr !< h_tr is h at tracer points with a h_neglect added to + !! ensure positive definiteness [H ~> m or kg m-2]. + 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]. + logical :: convert_flux = .true. + + + integer :: i, j, k, is, ie, js, je, nz + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + + if (nz == 1) then + call MOM_error(WARNING, "MOM_tracer_diabatic.F90, tracer_vertdiff called "//& + "with only one vertical level") + return + endif + + if (present(convert_flux_in)) convert_flux = convert_flux_in + h_neglect = GV%H_subroundoff + sink_dist = 0.0 + if (present(sink_rate)) sink_dist = (dt*sink_rate) * GV%m_to_H + !$OMP parallel default(shared) private(sink,h_minus_dsink,b_denom_1,b1,d1,h_tr,c1) + !$OMP do + do j=js,je; do i=is,ie ; sfc_src(i,j) = 0.0 ; btm_src(i,j) = 0.0 ; enddo ; enddo + if (present(sfc_flux)) then + if (convert_flux) then + !$OMP do + do j = js, je; do i = is,ie + sfc_src(i,j) = (sfc_flux(i,j)*dt) * GV%kg_m2_to_H + enddo ; enddo + else + !$OMP do + do j = js, je; do i = is,ie + sfc_src(i,j) = sfc_flux(i,j) + enddo ; enddo + endif + endif + if (present(btm_flux)) then + if (convert_flux) then + !$OMP do + do j = js, je; do i = is,ie + btm_src(i,j) = (btm_flux(i,j)*dt) * GV%kg_m2_to_H + enddo ; enddo + else + !$OMP do + do j = js, je; do i = is,ie + btm_src(i,j) = btm_flux(i,j) + enddo ; enddo + endif + endif + + if (present(sink_rate)) then + !$OMP do + do j=js,je + ! Find the sinking rates at all interfaces, limiting them if necesary + ! so that the characteristics do not cross within a timestep. + ! If a non-constant sinking rate were used, that would be incorprated + ! here. + if (present(btm_reservoir)) then + do i=is,ie ; sink(i,nz+1) = sink_dist ; enddo + do k=2,nz ; do i=is,ie + sink(i,K) = sink_dist ; h_minus_dsink(i,k) = h_old(i,j,k) + enddo ; enddo + else + do i=is,ie ; sink(i,nz+1) = 0.0 ; enddo + ! Find the limited sinking distance at the interfaces. + do k=nz,2,-1 ; do i=is,ie + if (sink(i,K+1) >= sink_dist) then + sink(i,K) = sink_dist + h_minus_dsink(i,k) = h_old(i,j,k) + (sink(i,K+1) - sink(i,K)) + elseif (sink(i,K+1) + h_old(i,j,k) < sink_dist) then + sink(i,K) = sink(i,K+1) + h_old(i,j,k) + h_minus_dsink(i,k) = 0.0 + else + sink(i,K) = sink_dist + h_minus_dsink(i,k) = (h_old(i,j,k) + sink(i,K+1)) - sink(i,K) + endif + enddo ; enddo + endif + do i=is,ie + sink(i,1) = 0.0 ; h_minus_dsink(i,1) = (h_old(i,j,1) + sink(i,2)) + enddo + + ! Now solve the tridiagonal equation for the tracer concentrations. + do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then + b_denom_1 = h_minus_dsink(i,1) + ent(i,j,1) + h_neglect + b1(i) = 1.0 / (b_denom_1 + ent(i,j,2)) + d1(i) = b_denom_1 * b1(i) + h_tr = h_old(i,j,1) + h_neglect + tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + sfc_src(i,j) + endif ; enddo + do k=2,nz-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then + c1(i,k) = ent(i,j,K) * b1(i) + b_denom_1 = h_minus_dsink(i,k) + d1(i) * (ent(i,j,K) + sink(i,K)) + & + h_neglect + b1(i) = 1.0 / (b_denom_1 + ent(i,j,K+1)) + d1(i) = b_denom_1 * b1(i) + h_tr = h_old(i,j,k) + h_neglect + tr(i,j,k) = b1(i) * (h_tr * tr(i,j,k) + & + (ent(i,j,K) + sink(i,K)) * tr(i,j,k-1)) + endif ; enddo ; enddo + do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then + c1(i,nz) = ent(i,j,nz) * b1(i) + b_denom_1 = h_minus_dsink(i,nz) + d1(i) * (ent(i,j,nz) + sink(i,nz)) + & + h_neglect + b1(i) = 1.0 / (b_denom_1 + ent(i,j,nz+1)) + h_tr = h_old(i,j,nz) + h_neglect + tr(i,j,nz) = b1(i) * ((h_tr * tr(i,j,nz) + btm_src(i,j)) + & + (ent(i,j,nz) + sink(i,nz)) * tr(i,j,nz-1)) + endif ; enddo + if (present(btm_reservoir)) then ; do i=is,ie ; if (G%mask2dT(i,j)>0.5) then + btm_reservoir(i,j) = btm_reservoir(i,j) + & + (sink(i,nz+1)*tr(i,j,nz)) * GV%H_to_kg_m2 + endif ; enddo ; endif + + do k=nz-1,1,-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then + tr(i,j,k) = tr(i,j,k) + c1(i,k+1)*tr(i,j,k+1) + endif ; enddo ; enddo + enddo + else + !$OMP do + do j=js,je + do i=is,ie ; if (G%mask2dT(i,j) > -0.5) then + h_tr = h_old(i,j,1) + h_neglect + b_denom_1 = h_tr + ent(i,j,1) + b1(i) = 1.0 / (b_denom_1 + ent(i,j,2)) + d1(i) = h_tr * b1(i) + tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + sfc_src(i,j) + endif + enddo + do k=2,nz-1 ; do i=is,ie ; if (G%mask2dT(i,j) > -0.5) then + c1(i,k) = ent(i,j,K) * b1(i) + h_tr = h_old(i,j,k) + h_neglect + b_denom_1 = h_tr + d1(i) * ent(i,j,K) + b1(i) = 1.0 / (b_denom_1 + ent(i,j,K+1)) + d1(i) = b_denom_1 * b1(i) + tr(i,j,k) = b1(i) * (h_tr * tr(i,j,k) + ent(i,j,K) * tr(i,j,k-1)) + endif ; enddo ; enddo + do i=is,ie ; if (G%mask2dT(i,j) > -0.5) then + c1(i,nz) = ent(i,j,nz) * b1(i) + h_tr = h_old(i,j,nz) + h_neglect + b_denom_1 = h_tr + d1(i)*ent(i,j,nz) + b1(i) = 1.0 / ( b_denom_1 + ent(i,j,nz+1)) + tr(i,j,nz) = b1(i) * (( h_tr * tr(i,j,nz) + btm_src(i,j)) + & + ent(i,j,nz) * tr(i,j,nz-1)) + endif ; enddo + do k=nz-1,1,-1 ; do i=is,ie ; if (G%mask2dT(i,j) > -0.5) then + tr(i,j,k) = tr(i,j,k) + c1(i,k+1)*tr(i,j,k+1) + endif ; enddo ; enddo + enddo + endif + !$OMP end parallel + +end subroutine tracer_vertdiff_Eulerian + + !> This routine is modeled after applyBoundaryFluxesInOut in MOM_diabatic_aux.F90 !! NOTE: Please note that in this routine sfc_flux gets set to zero to ensure that the surface !! flux of the tracer does not get applied again during a subsequent call to tracer_vertdif diff --git a/src/tracer/MOM_tracer_hor_diff.F90 b/src/tracer/MOM_tracer_hor_diff.F90 index 43ede7cff5..6506ca4927 100644 --- a/src/tracer/MOM_tracer_hor_diff.F90 +++ b/src/tracer/MOM_tracer_hor_diff.F90 @@ -11,7 +11,7 @@ module MOM_tracer_hor_diff use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type use MOM_domains, only : pass_vector use MOM_debugging, only : hchksum, uvchksum -use MOM_diabatic_driver, only : diabatic_CS +use MOM_diabatic_driver, only : diabatic_CS use MOM_EOS, only : calculate_density, EOS_type, EOS_domain use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg, is_root_pe use MOM_error_handler, only : MOM_set_verbosity, callTree_showQuery @@ -60,7 +60,7 @@ module MOM_tracer_hor_diff logical :: use_neutral_diffusion !< If true, use the neutral_diffusion module from within !! tracer_hor_diff. logical :: use_lateral_boundary_diffusion !< If true, use the lateral_boundary_diffusion module from within - !! tracer_hor_diff. + !! tracer_hor_diff. logical :: recalc_neutral_surf !< If true, recalculate the neutral surfaces if CFL has been !! exceeded type(neutral_diffusion_CS), pointer :: neutral_diffusion_CSp => NULL() !< Control structure for neutral diffusion. @@ -1513,7 +1513,7 @@ subroutine tracer_hor_diff_init(Time, G, US, param_file, diag, EOS, diabatic_CSp "USE_NEUTRAL_DIFFUSION and DIFFUSE_ML_TO_INTERIOR are mutually exclusive!") CS%use_lateral_boundary_diffusion = lateral_boundary_diffusion_init(Time, G, param_file, diag, diabatic_CSp, & CS%lateral_boundary_diffusion_CSp) - if (CS%use_neutral_diffusion .and. CS%Diffuse_ML_interior) call MOM_error(FATAL, "MOM_tracer_hor_diff: "// & + if (CS%use_lateral_boundary_diffusion .and. CS%Diffuse_ML_interior) call MOM_error(FATAL, "MOM_tracer_hor_diff: "// & "USE_LATERAL_BOUNDARY_DIFFUSION and DIFFUSE_ML_TO_INTERIOR are mutually exclusive!") call get_param(param_file, mdl, "DEBUG", CS%debug, default=.false.)