diff --git a/config_src/mct_driver/mom_ocean_model_mct.F90 b/config_src/mct_driver/mom_ocean_model_mct.F90 index 2f94c9b7f9..5a04739971 100644 --- a/config_src/mct_driver/mom_ocean_model_mct.F90 +++ b/config_src/mct_driver/mom_ocean_model_mct.F90 @@ -56,7 +56,7 @@ module MOM_ocean_model_mct use coupler_types_mod, only : coupler_type_set_diags, coupler_type_send_data use mpp_domains_mod, only : domain2d, mpp_get_layout, mpp_get_global_domain use mpp_domains_mod, only : mpp_define_domains, mpp_get_compute_domain, mpp_get_data_domain -use fms_mod, only : stdout +use MOM_io, only : stdout use mpp_mod, only : mpp_chksum use MOM_EOS, only : gsw_sp_from_sr, gsw_pt_from_ct use MOM_wave_interface, only : wave_parameters_CS, MOM_wave_interface_init @@ -409,10 +409,6 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i call close_param_file(param_file) call diag_mediator_close_registration(OS%diag) - - if (is_root_pe()) & - write(*,'(/12x,a/)') '======== COMPLETED MOM INITIALIZATION ========' - call callTree_leave("ocean_model_init(") end subroutine ocean_model_init @@ -1053,20 +1049,18 @@ subroutine ocean_public_type_chksum(id, timestep, ocn) integer, intent(in) :: timestep !< The number of elapsed timesteps type(ocean_public_type), intent(in) :: ocn !< A structure containing various publicly !! visible ocean surface fields. - integer :: n, m, outunit - - outunit = stdout() - - write(outunit,*) "BEGIN CHECKSUM(ocean_type):: ", id, timestep - write(outunit,100) 'ocean%t_surf ',mpp_chksum(ocn%t_surf ) - write(outunit,100) 'ocean%s_surf ',mpp_chksum(ocn%s_surf ) - write(outunit,100) 'ocean%u_surf ',mpp_chksum(ocn%u_surf ) - write(outunit,100) 'ocean%v_surf ',mpp_chksum(ocn%v_surf ) - write(outunit,100) 'ocean%sea_lev ',mpp_chksum(ocn%sea_lev) - write(outunit,100) 'ocean%frazil ',mpp_chksum(ocn%frazil ) - write(outunit,100) 'ocean%melt_potential ',mpp_chksum(ocn%melt_potential) - - call coupler_type_write_chksums(ocn%fields, outunit, 'ocean%') + integer :: n, m + + write(stdout,*) "BEGIN CHECKSUM(ocean_type):: ", id, timestep + write(stdout,100) 'ocean%t_surf ',mpp_chksum(ocn%t_surf ) + write(stdout,100) 'ocean%s_surf ',mpp_chksum(ocn%s_surf ) + write(stdout,100) 'ocean%u_surf ',mpp_chksum(ocn%u_surf ) + write(stdout,100) 'ocean%v_surf ',mpp_chksum(ocn%v_surf ) + write(stdout,100) 'ocean%sea_lev ',mpp_chksum(ocn%sea_lev) + write(stdout,100) 'ocean%frazil ',mpp_chksum(ocn%frazil ) + write(stdout,100) 'ocean%melt_potential ',mpp_chksum(ocn%melt_potential) + + call coupler_type_write_chksums(ocn%fields, stdout, 'ocean%') 100 FORMAT(" CHECKSUM::",A20," = ",Z20) end subroutine ocean_public_type_chksum diff --git a/config_src/mct_driver/mom_surface_forcing_mct.F90 b/config_src/mct_driver/mom_surface_forcing_mct.F90 index 92b5d148bb..82105e040e 100644 --- a/config_src/mct_driver/mom_surface_forcing_mct.F90 +++ b/config_src/mct_driver/mom_surface_forcing_mct.F90 @@ -34,10 +34,10 @@ module MOM_surface_forcing_mct use coupler_types_mod, only : coupler_type_initialized, coupler_type_spawn use coupler_types_mod, only : coupler_type_copy_data use data_override_mod, only : data_override_init, data_override -use fms_mod, only : stdout use mpp_mod, only : mpp_chksum use time_interp_external_mod, only : init_external_field, time_interp_external use time_interp_external_mod, only : time_interp_external_init +use MOM_io, only: stdout implicit none ; private @@ -1361,37 +1361,35 @@ subroutine ice_ocn_bnd_type_chksum(id, timestep, iobt) !! ocean in a coupled model whose checksums are reported ! local variables - integer :: n,m, outunit - - outunit = stdout() - - write(outunit,*) "BEGIN CHECKSUM(ice_ocean_boundary_type):: ", id, timestep - write(outunit,100) 'iobt%u_flux ' , mpp_chksum( iobt%u_flux ) - write(outunit,100) 'iobt%v_flux ' , mpp_chksum( iobt%v_flux ) - write(outunit,100) 'iobt%t_flux ' , mpp_chksum( iobt%t_flux ) - write(outunit,100) 'iobt%q_flux ' , mpp_chksum( iobt%q_flux ) - write(outunit,100) 'iobt%salt_flux ' , mpp_chksum( iobt%salt_flux ) - write(outunit,100) 'iobt%seaice_melt_heat' , mpp_chksum( iobt%seaice_melt_heat) - write(outunit,100) 'iobt%seaice_melt ' , mpp_chksum( iobt%seaice_melt ) - write(outunit,100) 'iobt%lw_flux ' , mpp_chksum( iobt%lw_flux ) - write(outunit,100) 'iobt%sw_flux_vis_dir' , mpp_chksum( iobt%sw_flux_vis_dir) - write(outunit,100) 'iobt%sw_flux_vis_dif' , mpp_chksum( iobt%sw_flux_vis_dif) - write(outunit,100) 'iobt%sw_flux_nir_dir' , mpp_chksum( iobt%sw_flux_nir_dir) - write(outunit,100) 'iobt%sw_flux_nir_dif' , mpp_chksum( iobt%sw_flux_nir_dif) - write(outunit,100) 'iobt%lprec ' , mpp_chksum( iobt%lprec ) - write(outunit,100) 'iobt%fprec ' , mpp_chksum( iobt%fprec ) - write(outunit,100) 'iobt%runoff ' , mpp_chksum( iobt%runoff ) - write(outunit,100) 'iobt%calving ' , mpp_chksum( iobt%calving ) - write(outunit,100) 'iobt%p ' , mpp_chksum( iobt%p ) + integer :: n,m + + write(stdout,*) "BEGIN CHECKSUM(ice_ocean_boundary_type):: ", id, timestep + write(stdout,100) 'iobt%u_flux ' , mpp_chksum( iobt%u_flux ) + write(stdout,100) 'iobt%v_flux ' , mpp_chksum( iobt%v_flux ) + write(stdout,100) 'iobt%t_flux ' , mpp_chksum( iobt%t_flux ) + write(stdout,100) 'iobt%q_flux ' , mpp_chksum( iobt%q_flux ) + write(stdout,100) 'iobt%salt_flux ' , mpp_chksum( iobt%salt_flux ) + write(stdout,100) 'iobt%seaice_melt_heat' , mpp_chksum( iobt%seaice_melt_heat) + write(stdout,100) 'iobt%seaice_melt ' , mpp_chksum( iobt%seaice_melt ) + write(stdout,100) 'iobt%lw_flux ' , mpp_chksum( iobt%lw_flux ) + write(stdout,100) 'iobt%sw_flux_vis_dir' , mpp_chksum( iobt%sw_flux_vis_dir) + write(stdout,100) 'iobt%sw_flux_vis_dif' , mpp_chksum( iobt%sw_flux_vis_dif) + write(stdout,100) 'iobt%sw_flux_nir_dir' , mpp_chksum( iobt%sw_flux_nir_dir) + write(stdout,100) 'iobt%sw_flux_nir_dif' , mpp_chksum( iobt%sw_flux_nir_dif) + write(stdout,100) 'iobt%lprec ' , mpp_chksum( iobt%lprec ) + write(stdout,100) 'iobt%fprec ' , mpp_chksum( iobt%fprec ) + write(stdout,100) 'iobt%runoff ' , mpp_chksum( iobt%runoff ) + write(stdout,100) 'iobt%calving ' , mpp_chksum( iobt%calving ) + write(stdout,100) 'iobt%p ' , mpp_chksum( iobt%p ) if (associated(iobt%ustar_berg)) & - write(outunit,100) 'iobt%ustar_berg ' , mpp_chksum( iobt%ustar_berg ) + write(stdout,100) 'iobt%ustar_berg ' , mpp_chksum( iobt%ustar_berg ) if (associated(iobt%area_berg)) & - write(outunit,100) 'iobt%area_berg ' , mpp_chksum( iobt%area_berg ) + write(stdout,100) 'iobt%area_berg ' , mpp_chksum( iobt%area_berg ) if (associated(iobt%mass_berg)) & - write(outunit,100) 'iobt%mass_berg ' , mpp_chksum( iobt%mass_berg ) + write(stdout,100) 'iobt%mass_berg ' , mpp_chksum( iobt%mass_berg ) 100 FORMAT(" CHECKSUM::",A20," = ",Z20) - call coupler_type_write_chksums(iobt%fluxes, outunit, 'iobt%') + call coupler_type_write_chksums(iobt%fluxes, stdout, 'iobt%') end subroutine ice_ocn_bnd_type_chksum diff --git a/config_src/mct_driver/ocn_comp_mct.F90 b/config_src/mct_driver/ocn_comp_mct.F90 index 741ce832e8..1872fff335 100644 --- a/config_src/mct_driver/ocn_comp_mct.F90 +++ b/config_src/mct_driver/ocn_comp_mct.F90 @@ -42,6 +42,7 @@ module ocn_comp_mct use MOM_constants, only: CELSIUS_KELVIN_OFFSET use MOM_domains, only: AGRID, BGRID_NE, CGRID_NE, pass_vector use mpp_domains_mod, only: mpp_get_compute_domain +use MOM_io, only: stdout ! Previously inlined - now in separate modules use MOM_ocean_model_mct, only: ocean_public_type, ocean_state_type @@ -88,7 +89,6 @@ module ocn_comp_mct type(cpl_indices_type) :: ind !< Variable IDs logical :: sw_decomp !< Controls whether shortwave is decomposed into 4 components real :: c1, c2, c3, c4 !< Coeffs. used in the shortwave decomposition i/o - integer :: stdout !< standard output unit. (by default, points to ocn.log.* ) character(len=384) :: pointer_filename !< Name of the ascii file that contains the path !! and filename of the latest restart file. end type MCT_MOM_Data @@ -194,14 +194,14 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) call shr_file_getLogUnit (shrlogunit) call shr_file_getLogLevel(shrloglev) - glb%stdout = shr_file_getUnit() ! get an unused unit number + stdout = shr_file_getUnit() ! get an unused unit number ! open the ocn_modelio.nml file and then open a log file associated with stdout ocn_modelio_name = 'ocn_modelio.nml' // trim(inst_suffix) - call shr_file_setIO(ocn_modelio_name,glb%stdout) + call shr_file_setIO(ocn_modelio_name,stdout) ! set the shr log io unit number - call shr_file_setLogUnit(glb%stdout) + call shr_file_setLogUnit(stdout) end if call set_calendar_type(NOLEAP) !TODO: confirm this @@ -218,23 +218,23 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) ! Debugging clocks if (debug .and. is_root_pe()) then - write(glb%stdout,*) 'ocn_init_mct, current time: y,m,d-',year,month,day,'h,m,s=',hour,minute,seconds + write(stdout,*) 'ocn_init_mct, current time: y,m,d-',year,month,day,'h,m,s=',hour,minute,seconds call ESMF_ClockGet(EClock, StartTime=time_var, rc=rc) call ESMF_TimeGet(time_var, yy=year, mm=month, dd=day, h=hour, m=minute, s=seconds, rc=rc) - write(glb%stdout,*) 'ocn_init_mct, start time: y,m,d-',year,month,day,'h,m,s=',hour,minute,seconds + write(stdout,*) 'ocn_init_mct, start time: y,m,d-',year,month,day,'h,m,s=',hour,minute,seconds call ESMF_ClockGet(EClock, StopTime=time_var, rc=rc) call ESMF_TimeGet(time_var, yy=year, mm=month, dd=day, h=hour, m=minute, s=seconds, rc=rc) - write(glb%stdout,*) 'ocn_init_mct, stop time: y,m,d-',year,month,day,'h,m,s=',hour,minute,seconds + write(stdout,*) 'ocn_init_mct, stop time: y,m,d-',year,month,day,'h,m,s=',hour,minute,seconds call ESMF_ClockGet(EClock, PrevTime=time_var, rc=rc) call ESMF_TimeGet(time_var, yy=year, mm=month, dd=day, h=hour, m=minute, s=seconds, rc=rc) - write(glb%stdout,*) 'ocn_init_mct, previous time: y,m,d-',year,month,day,'h,m,s=',hour,minute,seconds + write(stdout,*) 'ocn_init_mct, previous time: y,m,d-',year,month,day,'h,m,s=',hour,minute,seconds call ESMF_ClockGet(EClock, TimeStep=ocn_cpl_interval, rc=rc) call ESMF_TimeIntervalGet(ocn_cpl_interval, yy=year, mm=month, d=day, s=seconds, sn=seconds_n, sd=seconds_d, rc=rc) - write(glb%stdout,*) 'ocn_init_mct, time step: y,m,d-',year,month,day,'s,sn,sd=',seconds,seconds_n,seconds_d + write(stdout,*) 'ocn_init_mct, time step: y,m,d-',year,month,day,'s,sn,sd=',seconds,seconds_n,seconds_d endif npes = num_pes() @@ -298,7 +298,7 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) ! read name of restart file in the pointer file nu = shr_file_getUnit() restart_pointer_file = trim(glb%pointer_filename) - if (is_root_pe()) write(glb%stdout,*) 'Reading ocn pointer file: ',restart_pointer_file + if (is_root_pe()) write(stdout,*) 'Reading ocn pointer file: ',restart_pointer_file restartfile = ""; restartfiles = ""; open(nu, file=restart_pointer_file, form='formatted', status='unknown') do @@ -316,13 +316,13 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) enddo close(nu) if (is_root_pe()) then - write(glb%stdout,*) 'Reading restart file(s): ',trim(restartfiles) + write(stdout,*) 'Reading restart file(s): ',trim(restartfiles) end if call shr_file_freeUnit(nu) call ocean_model_init(glb%ocn_public, glb%ocn_state, time0, time_start, input_restart_file=trim(restartfiles)) endif if (is_root_pe()) then - write(glb%stdout,'(/12x,a/)') '======== COMPLETED MOM INITIALIZATION ========' + write(stdout,'(/12x,a/)') '======== COMPLETED MOM INITIALIZATION ========' end if ! Initialize ocn_state%sfc_state out of sight @@ -383,7 +383,7 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) ncouple_per_day = seconds_in_day / ocn_cpl_dt mom_cpl_dt = seconds_in_day / ncouple_per_day if (mom_cpl_dt /= ocn_cpl_dt) then - write(glb%stdout,*) 'ERROR mom_cpl_dt and ocn_cpl_dt must be identical' + write(stdout,*) 'ERROR mom_cpl_dt and ocn_cpl_dt must be identical' call exit(0) end if @@ -457,7 +457,7 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o) if (is_root_pe()) then call shr_file_getLogUnit(shrlogunit) call shr_file_getLogLevel(shrloglev) - call shr_file_setLogUnit(glb%stdout) + call shr_file_setLogUnit(stdout) endif ! Query the beginning time of the current coupling interval @@ -484,7 +484,7 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o) if (runtype /= "continue" .and. runtype /= "branch") then if (debug .and. is_root_pe()) then - write(glb%stdout,*) 'doubling first interval duration!' + write(stdout,*) 'doubling first interval duration!' endif ! shift back the start time by one coupling interval (to align the start time with other components) @@ -500,19 +500,19 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o) if (debug .and. is_root_pe()) then call ESMF_ClockGet(EClock, CurrTime=time_var, rc=rc) call ESMF_TimeGet(time_var, yy=year, mm=month, dd=day, h=hour, m=minute, s=seconds, rc=rc) - write(glb%stdout,*) 'ocn_run_mct, current time: y,m,d-',year,month,day,'h,m,s=',hour,minute,seconds + write(stdout,*) 'ocn_run_mct, current time: y,m,d-',year,month,day,'h,m,s=',hour,minute,seconds call ESMF_ClockGet(EClock, StartTime=time_var, rc=rc) call ESMF_TimeGet(time_var, yy=year, mm=month, dd=day, h=hour, m=minute, s=seconds, rc=rc) - write(glb%stdout,*) 'ocn_run_mct, start time: y,m,d-',year,month,day,'h,m,s=',hour,minute,seconds + write(stdout,*) 'ocn_run_mct, start time: y,m,d-',year,month,day,'h,m,s=',hour,minute,seconds call ESMF_ClockGet(EClock, StopTime=time_var, rc=rc) call ESMF_TimeGet(time_var, yy=year, mm=month, dd=day, h=hour, m=minute, s=seconds, rc=rc) - write(glb%stdout,*) 'ocn_run_mct, stop time: y,m,d-',year,month,day,'h,m,s=',hour,minute,seconds + write(stdout,*) 'ocn_run_mct, stop time: y,m,d-',year,month,day,'h,m,s=',hour,minute,seconds call ESMF_ClockGet(EClock, PrevTime=time_var, rc=rc) call ESMF_TimeGet(time_var, yy=year, mm=month, dd=day, h=hour, m=minute, s=seconds, rc=rc) - write(glb%stdout,*) 'ocn_run_mct, previous time: y,m,d-',year,month,day,'h,m,s=',hour,minute,seconds + write(stdout,*) 'ocn_run_mct, previous time: y,m,d-',year,month,day,'h,m,s=',hour,minute,seconds call ESMF_ClockGet(EClock, TimeStep=ocn_cpl_interval, rc=rc) call ESMF_TimeIntervalGet(ocn_cpl_interval, yy=year, mm=month, d=day, s=seconds, sn=seconds_n, sd=seconds_d, rc=rc) - write(glb%stdout,*) 'ocn_init_mct, time step: y,m,d-',year,month,day,'s,sn,sd=',seconds,seconds_n,seconds_d + write(stdout,*) 'ocn_init_mct, time step: y,m,d-',year,month,day,'s,sn,sd=',seconds,seconds_n,seconds_d endif ! set the cdata pointers: @@ -525,10 +525,10 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o) !glb%sw_decomp = .false. !END TODO: if (glb%sw_decomp) then - call ocn_import(x2o_o%rattr, glb%ind, glb%grid, Ice_ocean_boundary, glb%ocn_public, glb%stdout, Eclock, & + call ocn_import(x2o_o%rattr, glb%ind, glb%grid, Ice_ocean_boundary, glb%ocn_public, stdout, Eclock, & c1=glb%c1, c2=glb%c2, c3=glb%c3, c4=glb%c4) else - call ocn_import(x2o_o%rattr, glb%ind, glb%grid, Ice_ocean_boundary, glb%ocn_public, glb%stdout, Eclock ) + call ocn_import(x2o_o%rattr, glb%ind, glb%grid, Ice_ocean_boundary, glb%ocn_public, stdout, Eclock ) end if ! Update internal ocean @@ -540,7 +540,7 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o) !--- write out intermediate restart file when needed. ! Check alarms for flag to write restart at end of day write_restart_at_eod = seq_timemgr_RestartAlarmIsOn(EClock) - if (debug .and. is_root_pe()) write(glb%stdout,*) 'ocn_run_mct, write_restart_at_eod=', write_restart_at_eod + if (debug .and. is_root_pe()) write(stdout,*) 'ocn_run_mct, write_restart_at_eod=', write_restart_at_eod if (write_restart_at_eod) then ! case name @@ -575,7 +575,7 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o) endif close(nu) - write(glb%stdout,*) 'ocn restart pointer file written: ',trim(restartname) + write(stdout,*) 'ocn restart pointer file written: ',trim(restartname) endif call shr_file_freeUnit(nu) @@ -761,7 +761,7 @@ end subroutine ocn_domain_mct else if (trim(starttype) == trim(seq_infodata_start_type_brnch)) then get_runtype = "branch" else - write(glb%stdout,*) 'ocn_comp_mct ERROR: unknown starttype' + write(stdout,*) 'ocn_comp_mct ERROR: unknown starttype' call exit(0) end if return diff --git a/config_src/nuopc_driver/mom_cap.F90 b/config_src/nuopc_driver/mom_cap.F90 index c2a2e98838..fc6bb5035e 100644 --- a/config_src/nuopc_driver/mom_cap.F90 +++ b/config_src/nuopc_driver/mom_cap.F90 @@ -353,7 +353,7 @@ subroutine InitializeP0(gcomp, importState, exportState, clock, rc) file=__FILE__)) & return if (isPresent .and. isSet) then - read(value, '(i)', iostat=iostat) scalar_field_count + read(value, *, iostat=iostat) scalar_field_count if (iostat /= 0) then call ESMF_LogSetError(ESMF_RC_ARG_BAD, & msg=subname//": ScalarFieldCount not an integer: "//trim(value), & @@ -376,7 +376,7 @@ subroutine InitializeP0(gcomp, importState, exportState, clock, rc) file=__FILE__)) & return if (isPresent .and. isSet) then - read(value, '(i)', iostat=iostat) scalar_field_idx_grid_nx + read(value, *, iostat=iostat) scalar_field_idx_grid_nx if (iostat /= 0) then call ESMF_LogSetError(ESMF_RC_ARG_BAD, & msg=subname//": ScalarFieldIdxGridNX not an integer: "//trim(value), & @@ -399,7 +399,7 @@ subroutine InitializeP0(gcomp, importState, exportState, clock, rc) file=__FILE__)) & return if (isPresent .and. isSet) then - read(value, '(i)', iostat=iostat) scalar_field_idx_grid_ny + read(value, *, iostat=iostat) scalar_field_idx_grid_ny if (iostat /= 0) then call ESMF_LogSetError(ESMF_RC_ARG_BAD, & msg=subname//": ScalarFieldIdxGridNY not an integer: "//trim(value), & @@ -1434,14 +1434,14 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) !--------------------------------- if (len_trim(scalar_field_name) > 0) then - call State_SetScalar(dble(nxg),scalar_field_idx_grid_nx, exportState, localPet, & + call State_SetScalar(real(nxg,ESMF_KIND_R8),scalar_field_idx_grid_nx, exportState, localPet, & scalar_field_name, scalar_field_count, rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & return - call State_SetScalar(dble(nyg),scalar_field_idx_grid_ny, exportState, localPet, & + call State_SetScalar(real(nyg,ESMF_KIND_R8),scalar_field_idx_grid_ny, exportState, localPet, & scalar_field_name, scalar_field_count, rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & diff --git a/src/ALE/MOM_remapping.F90 b/src/ALE/MOM_remapping.F90 index 71ba83f3ba..1b3c5884de 100644 --- a/src/ALE/MOM_remapping.F90 +++ b/src/ALE/MOM_remapping.F90 @@ -13,8 +13,7 @@ module MOM_remapping use PLM_functions, only : PLM_reconstruction, PLM_boundary_extrapolation use PPM_functions, only : PPM_reconstruction, PPM_boundary_extrapolation use PQM_functions, only : PQM_reconstruction, PQM_boundary_extrapolation_v1 - -use iso_fortran_env, only : stdout=>output_unit, stderr=>error_unit +use MOM_io, only : stdout, stderr implicit none ; private @@ -1636,7 +1635,7 @@ logical function remapping_unit_tests(verbose) h_neglect = hNeglect_dflt h_neglect_edge = hNeglect_dflt ; if (answers_2018) h_neglect_edge = 1.0e-10 - write(*,*) '==== MOM_remapping: remapping_unit_tests =================' + write(stdout,*) '==== MOM_remapping: remapping_unit_tests =================' remapping_unit_tests = .false. ! Normally return false thisTest = .false. @@ -1645,19 +1644,19 @@ logical function remapping_unit_tests(verbose) err=x0(i)-0.75*real(i-1) if (abs(err)>real(i-1)*epsilon(err)) thisTest = .true. enddo - if (thisTest) write(*,*) 'remapping_unit_tests: Failed buildGridFromH() 1' + if (thisTest) write(stdout,*) 'remapping_unit_tests: Failed buildGridFromH() 1' remapping_unit_tests = remapping_unit_tests .or. thisTest call buildGridFromH(n1, h1, x1) do i=1,n1+1 err=x1(i)-real(i-1) if (abs(err)>real(i-1)*epsilon(err)) thisTest = .true. enddo - if (thisTest) write(*,*) 'remapping_unit_tests: Failed buildGridFromH() 2' + if (thisTest) write(stdout,*) 'remapping_unit_tests: Failed buildGridFromH() 2' remapping_unit_tests = remapping_unit_tests .or. thisTest thisTest = .false. call initialize_remapping(CS, 'PPM_H4', answers_2018=answers_2018) - if (verbose) write(*,*) 'h0 (test data)' + if (verbose) write(stdout,*) 'h0 (test data)' if (verbose) call dumpGrid(n0,h0,x0,u0) call dzFromH1H2( n0, h0, n1, h1, dx1 ) @@ -1666,9 +1665,9 @@ logical function remapping_unit_tests(verbose) err=u1(i)-8.*(0.5*real(1+n1)-real(i)) if (abs(err)>real(n1-1)*epsilon(err)) thisTest = .true. enddo - if (verbose) write(*,*) 'h1 (by projection)' + if (verbose) write(stdout,*) 'h1 (by projection)' if (verbose) call dumpGrid(n1,h1,x1,u1) - if (thisTest) write(*,*) 'remapping_unit_tests: Failed remapping_core_w()' + if (thisTest) write(stdout,*) 'remapping_unit_tests: Failed remapping_core_w()' remapping_unit_tests = remapping_unit_tests .or. thisTest thisTest = .false. @@ -1690,7 +1689,7 @@ logical function remapping_unit_tests(verbose) err=u1(i)-8.*(0.5*real(1+n1)-real(i)) if (abs(err)>2.*epsilon(err)) thisTest = .true. enddo - if (thisTest) write(*,*) 'remapping_unit_tests: Failed remapByProjection()' + if (thisTest) write(stdout,*) 'remapping_unit_tests: Failed remapByProjection()' remapping_unit_tests = remapping_unit_tests .or. thisTest thisTest = .false. @@ -1698,14 +1697,14 @@ logical function remapping_unit_tests(verbose) call remapByDeltaZ( n0, h0, u0, ppoly0_E, ppoly0_coefs, & n1, x1-x0(1:n1+1), & INTEGRATION_PPM, u1, hn1, h_neglect ) - if (verbose) write(*,*) 'h1 (by delta)' + if (verbose) write(stdout,*) 'h1 (by delta)' if (verbose) call dumpGrid(n1,h1,x1,u1) hn1=hn1-h1 do i=1,n1 err=u1(i)-8.*(0.5*real(1+n1)-real(i)) if (abs(err)>2.*epsilon(err)) thisTest = .true. enddo - if (thisTest) write(*,*) 'remapping_unit_tests: Failed remapByDeltaZ() 1' + if (thisTest) write(stdout,*) 'remapping_unit_tests: Failed remapByDeltaZ() 1' remapping_unit_tests = remapping_unit_tests .or. thisTest thisTest = .false. @@ -1715,19 +1714,19 @@ logical function remapping_unit_tests(verbose) call remapByDeltaZ( n0, h0, u0, ppoly0_E, ppoly0_coefs, & n2, dx2, & INTEGRATION_PPM, u2, hn2, h_neglect ) - if (verbose) write(*,*) 'h2' + if (verbose) write(stdout,*) 'h2' if (verbose) call dumpGrid(n2,h2,x2,u2) - if (verbose) write(*,*) 'hn2' + if (verbose) write(stdout,*) 'hn2' if (verbose) call dumpGrid(n2,hn2,x2,u2) do i=1,n2 err=u2(i)-8./2.*(0.5*real(1+n2)-real(i)) if (abs(err)>2.*epsilon(err)) thisTest = .true. enddo - if (thisTest) write(*,*) 'remapping_unit_tests: Failed remapByDeltaZ() 2' + if (thisTest) write(stdout,*) 'remapping_unit_tests: Failed remapByDeltaZ() 2' remapping_unit_tests = remapping_unit_tests .or. thisTest - if (verbose) write(*,*) 'Via sub-cells' + if (verbose) write(stdout,*) 'Via sub-cells' thisTest = .false. call remap_via_sub_cells( n0, h0, u0, ppoly0_E, ppoly0_coefs, & n2, h2, INTEGRATION_PPM, .false., u2, err ) @@ -1737,7 +1736,7 @@ logical function remapping_unit_tests(verbose) err=u2(i)-8./2.*(0.5*real(1+n2)-real(i)) if (abs(err)>2.*epsilon(err)) thisTest = .true. enddo - if (thisTest) write(*,*) 'remapping_unit_tests: Failed remap_via_sub_cells() 2' + if (thisTest) write(stdout,*) 'remapping_unit_tests: Failed remap_via_sub_cells() 2' remapping_unit_tests = remapping_unit_tests .or. thisTest call remap_via_sub_cells( n0, h0, u0, ppoly0_E, ppoly0_coefs, & @@ -1748,9 +1747,9 @@ logical function remapping_unit_tests(verbose) 3, (/2.25,1.5,1./), INTEGRATION_PPM, .false., u2, err ) if (verbose) call dumpGrid(3,h2,x2,u2) - if (.not. remapping_unit_tests) write(*,*) 'Pass' + if (.not. remapping_unit_tests) write(stdout,*) 'Pass' - write(*,*) '===== MOM_remapping: new remapping_unit_tests ==================' + write(stdout,*) '===== MOM_remapping: new remapping_unit_tests ==================' deallocate(ppoly0_E, ppoly0_S, ppoly0_coefs) allocate(ppoly0_coefs(5,6)) @@ -1879,7 +1878,7 @@ logical function remapping_unit_tests(verbose) deallocate(ppoly0_E, ppoly0_S, ppoly0_coefs) - if (.not. remapping_unit_tests) write(*,*) 'Pass' + if (.not. remapping_unit_tests) write(stdout,*) 'Pass' end function remapping_unit_tests diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index e7a3e54c4e..6930b2d4cb 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -2627,7 +2627,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & endif call tracer_advect_init(Time, G, US, param_file, diag, CS%tracer_adv_CSp) - call tracer_hor_diff_init(Time, G, US, param_file, diag, CS%tv%eqn_of_state, CS%diabatic_CSp, & + call tracer_hor_diff_init(Time, G, GV, US, param_file, diag, CS%tv%eqn_of_state, CS%diabatic_CSp, & CS%tracer_diff_CSp) call lock_tracer_registry(CS%tracer_Reg) diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 3936a788d0..18c9fc96c4 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -110,7 +110,8 @@ module MOM_diagnostics KE_dia => NULL() !< KE source from diapycnal diffusion [H L2 T-3 ~> m3 s-3] !>@{ Diagnostic IDs - integer :: id_u = -1, id_v = -1, id_h = -1 + integer :: id_u = -1, id_v = -1, id_h = -1 + integer :: id_usq = -1, id_vsq = -1, id_uv = -1 integer :: id_e = -1, id_e_D = -1 integer :: id_du_dt = -1, id_dv_dt = -1 ! integer :: id_hf_du_dt = -1, id_hf_dv_dt = -1 @@ -230,6 +231,10 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & !! calculating interface heights [H ~> m or kg m-2]. ! Local variables + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: usq ! squared eastward velocity [L2 T-2 ~> m2 s-2] + real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: vsq ! squared northward velocity [L2 T-2 ~> m2 s-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: uv ! u x v at h-points [L2 T-2 ~> m2 s-2] + integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state integer i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb @@ -337,6 +342,28 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & if (CS%id_h > 0) call post_data(CS%id_h, h, CS%diag) + if (CS%id_usq > 0) then + do k=1,nz ; do j=js,je ; do I=Isq,Ieq + usq(I,j,k) = u(I,j,k) * u(I,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_usq, usq, CS%diag) + endif + + if (CS%id_vsq > 0) then + do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + vsq(i,J,k) = v(i,J,k) * v(i,J,k) + enddo ; enddo ; enddo + call post_data(CS%id_vsq, vsq, CS%diag) + endif + + if (CS%id_uv > 0) then + do k=1,nz ; do j=js,je ; do i=is,ie + uv(i,j,k) = (0.5*(u(I-1,j,k) + u(I,j,k))) * & + (0.5*(v(i,J-1,k) + v(i,J,k))) + enddo ; enddo ; enddo + call post_data(CS%id_uv, uv, CS%diag) + endif + if (associated(CS%e)) then call find_eta(h, tv, G, GV, US, CS%e, eta_bt) if (CS%id_e > 0) call post_data(CS%id_e, CS%e, CS%diag) @@ -1641,6 +1668,13 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag CS%id_v = register_diag_field('ocean_model', 'v', diag%axesCvL, Time, & 'Meridional velocity', 'm s-1', conversion=US%L_T_to_m_s, cmor_field_name='vo', & cmor_standard_name='sea_water_y_velocity', cmor_long_name='Sea Water Y Velocity') + CS%id_usq = register_diag_field('ocean_model', 'usq', diag%axesCuL, Time, & + 'Zonal velocity squared', 'm2 s-2', conversion=US%L_T_to_m_s**2) + CS%id_vsq = register_diag_field('ocean_model', 'vsq', diag%axesCvL, Time, & + 'Meridional velocity squared', 'm2 s-2', conversion=US%L_T_to_m_s**2) + CS%id_uv = register_diag_field('ocean_model', 'uv', diag%axesTL, Time, & + 'Product between zonal and meridional velocities at h-points', 'm2 s-2', & + conversion=US%L_T_to_m_s**2) CS%id_h = register_diag_field('ocean_model', 'h', diag%axesTL, Time, & 'Layer Thickness', thickness_units, v_extensive=.true., conversion=convert_H) diff --git a/src/diagnostics/MOM_sum_output.F90 b/src/diagnostics/MOM_sum_output.F90 index 2d4fb7e06f..e1b1b8efaf 100644 --- a/src/diagnostics/MOM_sum_output.F90 +++ b/src/diagnostics/MOM_sum_output.F90 @@ -12,7 +12,7 @@ module MOM_sum_output use MOM_forcing_type, only : forcing use MOM_grid, only : ocean_grid_type use MOM_interface_heights, only : find_eta -use MOM_io, only : create_file, fieldtype, flush_file, open_file, reopen_file +use MOM_io, only : create_file, fieldtype, flush_file, open_file, reopen_file, stdout use MOM_io, only : file_exists, slasher, vardesc, var_desc, write_field, get_filename_appendix use MOM_io, only : APPEND_FILE, ASCII_FILE, SINGLE_FILE, WRITEONLY_FILE use MOM_open_boundary, only : ocean_OBC_type, OBC_segment_type @@ -827,11 +827,11 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, OBC, dt_ if (is_root_pe()) then if (CS%use_temperature) then - write(*,'(A," ",A,": En ",ES12.6, ", MaxCFL ", F8.5, ", Mass ", & + write(stdout,'(A," ",A,": En ",ES12.6, ", MaxCFL ", F8.5, ", Mass ", & & ES18.12, ", Salt ", F15.11,", Temp ", F15.11)') & trim(date_str), trim(n_str), En_mass, max_CFL(1), mass_tot, salin, temp else - write(*,'(A," ",A,": En ",ES12.6, ", MaxCFL ", F8.5, ", Mass ", & + write(stdout,'(A," ",A,": En ",ES12.6, ", MaxCFL ", F8.5, ", Mass ", & & ES18.12)') & trim(date_str), trim(n_str), En_mass, max_CFL(1), mass_tot endif @@ -853,39 +853,39 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, OBC, dt_ endif if (CS%ntrunc > 0) then - write(*,'(A," Energy/Mass:",ES12.5," Truncations ",I0)') & + write(stdout,'(A," Energy/Mass:",ES12.5," Truncations ",I0)') & trim(date_str), En_mass, CS%ntrunc endif if (CS%write_stocks) then - write(*,'(" Total Energy: ",Z16.16,ES24.16)') toten, toten - write(*,'(" Total Mass: ",ES24.16,", Change: ",ES24.16," Error: ",ES12.5," (",ES8.1,")")') & + write(stdout,'(" Total Energy: ",Z16.16,ES24.16)') toten, toten + write(stdout,'(" Total Mass: ",ES24.16,", Change: ",ES24.16," Error: ",ES12.5," (",ES8.1,")")') & mass_tot, mass_chg, mass_anom, mass_anom/mass_tot if (CS%use_temperature) then if (Salt == 0.) then - write(*,'(" Total Salt: ",ES24.16,", Change: ",ES24.16," Error: ",ES12.5)') & + write(stdout,'(" Total Salt: ",ES24.16,", Change: ",ES24.16," Error: ",ES12.5)') & Salt*0.001, Salt_chg*0.001, Salt_anom*0.001 else - write(*,'(" Total Salt: ",ES24.16,", Change: ",ES24.16," Error: ",ES12.5," (",ES8.1,")")') & + write(stdout,'(" Total Salt: ",ES24.16,", Change: ",ES24.16," Error: ",ES12.5," (",ES8.1,")")') & Salt*0.001, Salt_chg*0.001, Salt_anom*0.001, Salt_anom/Salt endif if (Heat == 0.) then - write(*,'(" Total Heat: ",ES24.16,", Change: ",ES24.16," Error: ",ES12.5)') & + write(stdout,'(" Total Heat: ",ES24.16,", Change: ",ES24.16," Error: ",ES12.5)') & Heat, Heat_chg, Heat_anom else - write(*,'(" Total Heat: ",ES24.16,", Change: ",ES24.16," Error: ",ES12.5," (",ES8.1,")")') & + write(stdout,'(" Total Heat: ",ES24.16,", Change: ",ES24.16," Error: ",ES12.5," (",ES8.1,")")') & Heat, Heat_chg, Heat_anom, Heat_anom/Heat endif endif do m=1,nTr_stocks - write(*,'(" Total ",a,": ",ES24.16,X,a)') & + write(stdout,'(" Total ",a,": ",ES24.16,X,a)') & trim(Tr_names(m)), Tr_stocks(m), trim(Tr_units(m)) if (Tr_minmax_got(m)) then - write(*,'(64X,"Global Min:",ES24.16,X,"at: (", f7.2,","f7.2,","f8.2,")" )') & + write(stdout,'(64X,"Global Min:",ES24.16,X,"at: (", f7.2,","f7.2,","f8.2,")" )') & Tr_min(m),Tr_min_x(m),Tr_min_y(m),Tr_min_z(m) - write(*,'(64X,"Global Max:",ES24.16,X,"at: (", f7.2,","f7.2,","f8.2,")" )') & + write(stdout,'(64X,"Global Max:",ES24.16,X,"at: (", f7.2,","f7.2,","f8.2,")" )') & Tr_max(m),Tr_max_x(m),Tr_max_y(m),Tr_max_z(m) endif diff --git a/src/framework/MOM_diag_vkernels.F90 b/src/framework/MOM_diag_vkernels.F90 index b7c1130521..3d6e3e3f65 100644 --- a/src/framework/MOM_diag_vkernels.F90 +++ b/src/framework/MOM_diag_vkernels.F90 @@ -4,7 +4,7 @@ module MOM_diag_vkernels ! This file is part of MOM6. See LICENSE.md for the license. -use iso_fortran_env, only : stdout=>output_unit, stderr=>error_unit +use MOM_io, only : stdout, stderr implicit none ; private diff --git a/src/framework/MOM_document.F90 b/src/framework/MOM_document.F90 index 15d0839ee9..ff0934ac55 100644 --- a/src/framework/MOM_document.F90 +++ b/src/framework/MOM_document.F90 @@ -661,7 +661,7 @@ end function real_string !> Returns a character string of a comma-separated, compact formatted, reals !> e.g. "1., 2., 5*3., 5.E2", that give the list of values. function real_array_string(vals, sep) - character(len=1320) :: real_array_string !< The output string listing vals + character(len=:) ,allocatable :: real_array_string !< The output string listing vals real, intent(in) :: vals(:) !< The array of values to record character(len=*), & optional, intent(in) :: sep !< The separator between successive values, @@ -669,10 +669,10 @@ function real_array_string(vals, sep) ! Returns a character string of a comma-separated, compact formatted, reals ! e.g. "1., 2., 5*3., 5.E2" ! Local variables - integer :: j, n, b, ns + integer :: j, n, ns logical :: doWrite character(len=10) :: separator - n=1 ; doWrite=.true. ; real_array_string='' ; b=1 + n=1 ; doWrite=.true. ; real_array_string='' if (present(sep)) then separator=sep ; ns=len(sep) else @@ -687,16 +687,15 @@ function real_array_string(vals, sep) endif endif if (doWrite) then - if (b>1) then ! Write separator if a number has already been written - write(real_array_string(b:),'(A)') separator - b=b+ns + if(len(real_array_string)>0) then ! Write separator if a number has already been written + real_array_string = real_array_string // separator(1:ns) endif if (n>1) then - write(real_array_string(b:),'(A,"*",A)') trim(int_string(n)),trim(real_string(vals(j))) + real_array_string = real_array_string // trim(int_string(n)) // "*" // trim(real_string(vals(j))) else - write(real_array_string(b:),'(A)') trim(real_string(vals(j))) + real_array_string = real_array_string // trim(real_string(vals(j))) endif - n=1 ; b=len_trim(real_array_string)+1 + n=1 endif enddo end function real_array_string diff --git a/src/framework/MOM_file_parser.F90 b/src/framework/MOM_file_parser.F90 index 2e7a14dbe4..522b0958c1 100644 --- a/src/framework/MOM_file_parser.F90 +++ b/src/framework/MOM_file_parser.F90 @@ -1406,14 +1406,13 @@ subroutine log_param_real_array(CS, modulename, varname, value, desc, & logical, optional, intent(in) :: like_default !< If present and true, log this parameter as !! though it has the default value, even if there is no default. - character(len=1320) :: mesg + character(len=:), allocatable :: mesg character(len=240) :: myunits !write(mesg, '(" ",a," ",a,": ",ES19.12,99(",",ES19.12))') & !write(mesg, '(" ",a," ",a,": ",G,99(",",G))') & ! trim(modulename), trim(varname), value - write(mesg, '(" ",a," ",a,": ",a)') & - trim(modulename), trim(varname), trim(left_reals(value)) + mesg = " " // trim(modulename) // " " // trim(varname) // ": " // trim(left_reals(value)) if (is_root_pe()) then if (CS%log_open) write(CS%stdlog,'(a)') trim(mesg) if (CS%log_to_stdout) write(CS%stdout,'(a)') trim(mesg) diff --git a/src/framework/MOM_io.F90 b/src/framework/MOM_io.F90 index c516c96e86..d13dddc3c7 100644 --- a/src/framework/MOM_io.F90 +++ b/src/framework/MOM_io.F90 @@ -33,6 +33,7 @@ module MOM_io use mpp_io_mod, only : get_file_fields=>mpp_get_fields, get_file_times=>mpp_get_times use mpp_io_mod, only : io_infra_init=>mpp_io_init +use iso_fortran_env, only : stdout_iso=>output_unit, stderr_iso=>error_unit use netcdf implicit none ; private @@ -84,6 +85,9 @@ module MOM_io module procedure MOM_read_vector_2d end interface +integer, public :: stdout = stdout_iso !< standard output unit +integer, public :: stderr = stderr_iso !< standard output unit + contains !> Routine creates a new NetCDF file. It also sets up diff --git a/src/framework/MOM_random.F90 b/src/framework/MOM_random.F90 index 14800df9aa..161236572c 100644 --- a/src/framework/MOM_random.F90 +++ b/src/framework/MOM_random.F90 @@ -11,7 +11,7 @@ module MOM_random use MersenneTwister_mod, only : getRandomReal ! Generates a random number use MersenneTwister_mod, only : getRandomPositiveInt ! Generates a random positive integer -use iso_fortran_env, only : stdout=>output_unit, stderr=>error_unit +use MOM_io, only : stdout, stderr implicit none ; private diff --git a/src/framework/MOM_string_functions.F90 b/src/framework/MOM_string_functions.F90 index 1293499930..5c04a77b7d 100644 --- a/src/framework/MOM_string_functions.F90 +++ b/src/framework/MOM_string_functions.F90 @@ -142,13 +142,13 @@ function left_reals(r,sep) real, intent(in) :: r(:) !< The array of real variables to convert to a string character(len=*), optional, intent(in) :: sep !< The separator between !! successive values, by default it is ', '. - character(len=1320) :: left_reals !< The output string + character(len=:), allocatable :: left_reals !< The output string - integer :: j, n, b, ns + integer :: j, n, ns logical :: doWrite character(len=10) :: separator - n=1 ; doWrite=.true. ; left_reals='' ; b=1 + n=1 ; doWrite=.true. ; left_reals='' if (present(sep)) then separator=sep ; ns=len(sep) else @@ -163,16 +163,15 @@ function left_reals(r,sep) endif endif if (doWrite) then - if (b>1) then ! Write separator if a number has already been written - write(left_reals(b:),'(A)') separator - b=b+ns + if (len(left_reals)>0) then ! Write separator if a number has already been written + left_reals = left_reals // separator(1:ns) endif if (n>1) then - write(left_reals(b:),'(A,"*",A)') trim(left_int(n)),trim(left_real(r(j))) + left_reals = left_reals // trim(left_int(n)) // "*" // trim(left_real(r(j))) else - write(left_reals(b:),'(A)') trim(left_real(r(j))) + left_reals = left_reals // trim(left_real(r(j))) endif - n=1 ; b=len_trim(left_reals)+1 + n=1 endif enddo end function left_reals diff --git a/src/initialization/MOM_grid_initialize.F90 b/src/initialization/MOM_grid_initialize.F90 index 88130857c7..2c9445ae3e 100644 --- a/src/initialization/MOM_grid_initialize.F90 +++ b/src/initialization/MOM_grid_initialize.F90 @@ -14,7 +14,7 @@ module MOM_grid_initialize use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, is_root_pe use MOM_error_handler, only : callTree_enter, callTree_leave use MOM_file_parser, only : get_param, log_param, log_version, param_file_type -use MOM_io, only : MOM_read_data, read_data, slasher, file_exists +use MOM_io, only : MOM_read_data, read_data, slasher, file_exists, stdout use MOM_io, only : CORNER, NORTH_FACE, EAST_FACE use MOM_unit_scaling, only : unit_scale_type @@ -806,14 +806,14 @@ subroutine set_grid_metrics_mercator(G, param_file, US) y_q = find_root(Int_dj_dy, dy_dj, GP, jd, 0.0, -1.0*PI_2, PI_2, itt2) G%gridLatB(J) = y_q*180.0/PI ! if (is_root_pe()) & - ! write(*, '("J, y_q = ",I4,ES14.4," itts = ",I4)') j, y_q, itt2 + ! write(stdout, '("J, y_q = ",I4,ES14.4," itts = ",I4)') j, y_q, itt2 enddo do j=G%jsg,G%jeg jd = fnRef + (j - jRef) - 0.5 y_h = find_root(Int_dj_dy, dy_dj, GP, jd, 0.0, -1.0*PI_2, PI_2, itt1) G%gridLatT(j) = y_h*180.0/PI ! if (is_root_pe()) & - ! write(*, '("j, y_h = ",I4,ES14.4," itts = ",I4)') j, y_h, itt1 + ! write(stdout, '("j, y_h = ",I4,ES14.4," itts = ",I4)') j, y_h, itt1 enddo do J=JsdB+J_off,JedB+J_off jd = fnRef + (J - jRef) diff --git a/src/initialization/MOM_shared_initialization.F90 b/src/initialization/MOM_shared_initialization.F90 index 51676fb54d..7528f3f33e 100644 --- a/src/initialization/MOM_shared_initialization.F90 +++ b/src/initialization/MOM_shared_initialization.F90 @@ -11,7 +11,7 @@ module MOM_shared_initialization use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, WARNING, is_root_pe use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint use MOM_file_parser, only : get_param, log_param, param_file_type, log_version -use MOM_io, only : close_file, create_file, fieldtype, file_exists +use MOM_io, only : close_file, create_file, fieldtype, file_exists, stdout use MOM_io, only : MOM_read_data, MOM_read_vector, SINGLE_FILE, MULTIPLE use MOM_io, only : slasher, vardesc, write_field, var_desc use MOM_string_functions, only : uppercase @@ -282,7 +282,7 @@ subroutine apply_topography_edits_from_file(D, G, param_file, US) j = jg(n) - G%jsd_global + 2 if (i>=G%isc .and. i<=G%iec .and. j>=G%jsc .and. j<=G%jec) then if (new_depth(n)/=0.) then - write(*,'(a,3i5,f8.2,a,f8.2,2i4)') & + write(stdout,'(a,3i5,f8.2,a,f8.2,2i4)') & 'Ocean topography edit: ',n,ig(n),jg(n),D(i,j)/m_to_Z,'->',abs(new_depth(n)),i,j D(i,j) = abs(m_to_Z*new_depth(n)) ! Allows for height-file edits (i.e. converts negatives) else @@ -995,10 +995,10 @@ subroutine reset_face_lengths_list(G, param_file, US) G%dy_Cu(I,j) = G%mask2dCu(I,j) * m_to_L*min(L_to_m*G%dyCu(I,j), max(u_width(npt), 0.0)) if (j>=G%jsc .and. j<=G%jec .and. I>=G%isc .and. I<=G%iec) then ! Limit messages/checking to compute domain if ( G%mask2dCu(I,j) == 0.0 ) then - write(*,'(A,2F8.2,A,4F8.2,A)') "read_face_lengths_list : G%mask2dCu=0 at ",lat,lon," (",& + write(stdout,'(A,2F8.2,A,4F8.2,A)') "read_face_lengths_list : G%mask2dCu=0 at ",lat,lon," (",& u_lat(1,npt), u_lat(2,npt), u_lon(1,npt), u_lon(2,npt),") so grid metric is unmodified." else - write(*,'(A,2F8.2,A,4F8.2,A5,F9.2,A1)') & + write(stdout,'(A,2F8.2,A,4F8.2,A5,F9.2,A1)') & "read_face_lengths_list : Modifying dy_Cu gridpoint at ",lat,lon," (",& u_lat(1,npt), u_lat(2,npt), u_lon(1,npt), u_lon(2,npt),") to ",L_to_m*G%dy_Cu(I,j),"m" endif @@ -1024,10 +1024,10 @@ subroutine reset_face_lengths_list(G, param_file, US) G%dx_Cv(i,J) = G%mask2dCv(i,J) * m_to_L*min(L_to_m*G%dxCv(i,J), max(v_width(npt), 0.0)) if (i>=G%isc .and. i<=G%iec .and. J>=G%jsc .and. J<=G%jec) then ! Limit messages/checking to compute domain if ( G%mask2dCv(i,J) == 0.0 ) then - write(*,'(A,2F8.2,A,4F8.2,A)') "read_face_lengths_list : G%mask2dCv=0 at ",lat,lon," (",& + write(stdout,'(A,2F8.2,A,4F8.2,A)') "read_face_lengths_list : G%mask2dCv=0 at ",lat,lon," (",& v_lat(1,npt), v_lat(2,npt), v_lon(1,npt), v_lon(2,npt),") so grid metric is unmodified." else - write(*,'(A,2F8.2,A,4F8.2,A5,F9.2,A1)') & + write(stdout,'(A,2F8.2,A,4F8.2,A5,F9.2,A1)') & "read_face_lengths_list : Modifying dx_Cv gridpoint at ",lat,lon," (",& v_lat(1,npt), v_lat(2,npt), v_lon(1,npt), v_lon(2,npt),") to ",L_to_m*G%dx_Cv(I,j),"m" endif diff --git a/src/tracer/MOM_lateral_boundary_diffusion.F90 b/src/tracer/MOM_lateral_boundary_diffusion.F90 index 465174f676..2b7a5646cc 100644 --- a/src/tracer/MOM_lateral_boundary_diffusion.F90 +++ b/src/tracer/MOM_lateral_boundary_diffusion.F90 @@ -6,24 +6,25 @@ module MOM_lateral_boundary_diffusion ! This file is part of MOM6. See LICENSE.md for the license. use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end -use MOM_cpu_clock, only : CLOCK_MODULE, CLOCK_ROUTINE -use MOM_domains, only : pass_var +use MOM_cpu_clock, only : CLOCK_MODULE +use MOM_checksums, only : hchksum +use MOM_domains, only : pass_var, sum_across_PEs use MOM_diag_mediator, only : diag_ctrl, time_type use MOM_diag_mediator, only : post_data, register_diag_field -use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg, is_root_pe +use MOM_diag_vkernels, only : reintegrate_column +use MOM_error_handler, only : MOM_error, FATAL, is_root_pe use MOM_file_parser, only : get_param, log_version, param_file_type -use MOM_file_parser, only : openParameterBlock, closeParameterBlock use MOM_grid, only : ocean_grid_type use MOM_remapping, only : remapping_CS, initialize_remapping -use MOM_remapping, only : extract_member_remapping_CS, build_reconstructions_1d -use MOM_remapping, only : average_value_ppoly, remappingSchemesDoc, remappingDefaultScheme +use MOM_remapping, only : extract_member_remapping_CS, remapping_core_h +use MOM_remapping, only : remappingSchemesDoc, remappingDefaultScheme use MOM_tracer_registry, only : tracer_registry_type, tracer_type use MOM_unit_scaling, only : unit_scale_type use MOM_verticalGrid, only : verticalGrid_type use MOM_CVMix_KPP, only : KPP_get_BLD, KPP_CS use MOM_energetic_PBL, only : energetic_PBL_get_MLD, energetic_PBL_CS use MOM_diabatic_driver, only : diabatic_CS, extract_diabatic_member -use iso_fortran_env, only : stdout=>output_unit, stderr=>error_unit +use MOM_io, only : stdout, stderr implicit none ; private @@ -36,45 +37,49 @@ module MOM_lateral_boundary_diffusion #include !> Sets parameters for lateral boundary mixing module. -type, public :: lateral_boundary_diffusion_CS ; private - integer :: method !< Determine which of the three methods calculate - !! and apply near boundary layer fluxes - !! 1. Along layer - !! 2. Bulk-layer approach (not recommended) - integer :: deg !< Degree of polynomial reconstruction - integer :: surface_boundary_scheme !< Which boundary layer scheme to use - !! 1. ePBL; 2. KPP - logical :: limiter !< Controls wether a flux limiter is applied. - !! Only valid when method = 2. - logical :: linear !< If True, apply a linear transition at the base/top of the boundary. - !! The flux will be fully applied at k=k_min and zero at k=k_max. - - type(remapping_CS) :: remap_CS !< Control structure to hold remapping configuration - type(KPP_CS), pointer :: KPP_CSp => NULL() !< KPP control structure needed to get BLD - type(energetic_PBL_CS), pointer :: energetic_PBL_CSp => NULL() !< ePBL control structure needed to get BLD - type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to - !! regulate the timing of diagnostic output. -end type lateral_boundary_diffusion_CS +type, public :: lbd_CS ; private + logical :: debug !< If true, write verbose checksums for debugging. + integer :: deg !< Degree of polynomial reconstruction. + integer :: surface_boundary_scheme !< Which boundary layer scheme to use + !! 1. ePBL; 2. KPP + logical :: limiter !< Controls whether a flux limiter is applied in the + !! native grid (default is true). + logical :: limiter_remap !< Controls whether a flux limiter is applied in the + !! remapped grid (default is false). + logical :: linear !< If True, apply a linear transition at the base/top of the boundary. + !! The flux will be fully applied at k=k_min and zero at k=k_max. + real :: H_subroundoff !< A thickness that is so small that it can be added to a thickness of + !! Angstrom or larger without changing it at the bit level [H ~> m or kg m-2]. + !! If Angstrom is 0 or exceedingly small, this is negligible compared to 1e-17 m. + type(remapping_CS) :: remap_CS !< Control structure to hold remapping configuration. + type(KPP_CS), pointer :: KPP_CSp => NULL() !< KPP control structure needed to get BLD. + type(energetic_PBL_CS), pointer :: energetic_PBL_CSp => NULL() !< ePBL control structure needed to get BLD. + type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to + !! regulate the timing of diagnostic output. +end type lbd_CS ! This include declares and sets the variable "version". #include "version_variable.h" -character(len=40) :: mdl = "MOM_lateral_boundary_diffusion" !< Name of this module +character(len=40) :: mdl = "MOM_lateral_boundary_diffusion" !< Name of this module +integer :: id_clock_lbd !< CPU clock for lbd contains !> Initialization routine that reads runtime parameters and sets up pointers to other control structures that might be !! needed for lateral boundary diffusion. -logical function lateral_boundary_diffusion_init(Time, G, param_file, diag, diabatic_CSp, CS) +logical function lateral_boundary_diffusion_init(Time, G, GV, param_file, diag, diabatic_CSp, CS) type(time_type), target, intent(in) :: Time !< Time structure type(ocean_grid_type), intent(in) :: G !< Grid structure + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure type(param_file_type), intent(in) :: param_file !< Parameter file structure type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure type(diabatic_CS), pointer :: diabatic_CSp !< KPP control structure needed to get BLD - type(lateral_boundary_diffusion_CS), pointer :: CS !< Lateral boundary mixing control structure + type(lbd_CS), pointer :: CS !< Lateral boundary mixing control structure ! local variables - character(len=80) :: string ! Temporary strings - logical :: boundary_extrap + character(len=80) :: string ! Temporary strings + integer :: ke, nk ! Number of levels in the LBD and native grids, respectively + logical :: boundary_extrap ! controls if boundary extrapolation is used in the LBD code if (ASSOCIATED(CS)) then call MOM_error(FATAL, "lateral_boundary_diffusion_init called with associated control structure.") @@ -90,11 +95,11 @@ logical function lateral_boundary_diffusion_init(Time, G, param_file, diag, diab call get_param(param_file, mdl, "USE_LATERAL_BOUNDARY_DIFFUSION", lateral_boundary_diffusion_init, & "If true, enables the lateral boundary tracer's diffusion module.", & default=.false.) - if (.not. lateral_boundary_diffusion_init) return allocate(CS) CS%diag => diag + CS%H_subroundoff = GV%H_subroundoff call extract_diabatic_member(diabatic_CSp, KPP_CSp=CS%KPP_CSp) call extract_diabatic_member(diabatic_CSp, energetic_PBL_CSp=CS%energetic_PBL_CSp) @@ -104,18 +109,13 @@ logical function lateral_boundary_diffusion_init(Time, G, param_file, diag, diab endif ! Read all relevant parameters and write them to the model log. - call get_param(param_file, mdl, "LATERAL_BOUNDARY_METHOD", CS%method, & - "Determine how to apply boundary lateral diffusion of tracers: \n"//& - "1. Along layer approach \n"//& - "2. Bulk layer approach (this option is not recommended)", default=1) - if (CS%method == 2) then - call get_param(param_file, mdl, "APPLY_LIMITER", CS%limiter, & - "If True, apply a flux limiter in the LBD. This is only available \n"//& - "when LATERAL_BOUNDARY_METHOD=2.", default=.false.) - endif call get_param(param_file, mdl, "LBD_LINEAR_TRANSITION", CS%linear, & "If True, apply a linear transition at the base/top of the boundary. \n"//& "The flux will be fully applied at k=k_min and zero at k=k_max.", default=.false.) + call get_param(param_file, mdl, "APPLY_LIMITER", CS%limiter, & + "If True, apply a flux limiter in the native grid.", default=.true.) + call get_param(param_file, mdl, "APPLY_LIMITER_REMAP", CS%limiter_remap, & + "If True, apply a flux limiter in the remapped grid.", default=.false.) call get_param(param_file, mdl, "LBD_BOUNDARY_EXTRAP", boundary_extrap, & "Use boundary extrapolation in LBD code", & default=.false.) @@ -124,20 +124,27 @@ logical function lateral_boundary_diffusion_init(Time, G, param_file, diag, diab "for vertical remapping for all variables. "//& "It can be one of the following schemes: "//& trim(remappingSchemesDoc), default=remappingDefaultScheme) - call initialize_remapping( CS%remap_CS, string, boundary_extrapolation = boundary_extrap ) + call initialize_remapping( CS%remap_CS, string, boundary_extrapolation = boundary_extrap ,& + check_reconstruction = .false., check_remapping = .false., answers_2018 = .false.) call extract_member_remapping_CS(CS%remap_CS, degree=CS%deg) + call get_param(param_file, mdl, "LBD_DEBUG", CS%debug, & + "If true, write out verbose debugging data in the LBD module.", & + default=.false.) + + id_clock_lbd = cpu_clock_id('(Ocean LBD)', grain=CLOCK_MODULE) end function lateral_boundary_diffusion_init !> Driver routine for calculating lateral diffusive fluxes near the top and bottom boundaries. -!! Two different methods are available: -!! Method 1: lower order representation, calculate fluxes from bulk layer integrated quantities. -!! Method 2: more straight forward, diffusion is applied layer by layer using only information -!! from neighboring cells. +!! Diffusion is applied using only information from neighboring cells, as follows: +!! 1) remap tracer to a z* grid (LBD grid) +!! 2) calculate diffusive tracer fluxes (F) in the LBD grid using a layer by layer approach +!! 3) remap fluxes to the native grid +!! 4) update tracer by adding the divergence of F subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) - type(ocean_grid_type), intent(inout) :: G !< Grid type - type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(ocean_grid_type), intent(inout) :: G !< Grid type + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & intent(in) :: h !< Layer thickness [H ~> m or kg m-2] real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: Coef_x !< dt * Kh * dy / dx at u-points [L2 ~> m2] @@ -145,118 +152,102 @@ subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) real, intent(in) :: dt !< Tracer time step * I_numitts !! (I_numitts in tracer_hordiff) type(tracer_registry_type), pointer :: Reg !< Tracer registry - type(lateral_boundary_diffusion_CS), intent(in) :: CS !< Control structure for this module + type(lbd_CS), pointer :: CS !< Control structure for this module ! Local variables - real, dimension(SZI_(G),SZJ_(G)) :: hbl !< bnd. layer depth [H ~> m or kg m-2] - real, dimension(SZI_(G),SZJ_(G),SZK_(G),CS%deg+1) :: ppoly0_coefs !< Coefficients of polynomial - real, dimension(SZI_(G),SZJ_(G),SZK_(G),2) :: ppoly0_E !< Edge values from reconstructions - real, dimension(SZK_(G),CS%deg+1) :: ppoly_S !< Slopes from reconstruction (placeholder) + real, dimension(SZI_(G),SZJ_(G)) :: hbl !< bnd. layer depth [H ~> m or kg m-2] real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: uFlx !< Zonal flux of tracer [conc m^3] - real, dimension(SZIB_(G),SZJ_(G)) :: uFLx_bulk !< Total calculated bulk-layer u-flux for the tracer real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: vFlx !< Meridional flux of tracer [conc m^3] - real, dimension(SZI_(G),SZJB_(G)) :: vFlx_bulk !< Total calculated bulk-layer v-flux for the tracer real, dimension(SZIB_(G),SZJ_(G)) :: uwork_2d !< Layer summed u-flux transport real, dimension(SZI_(G),SZJB_(G)) :: vwork_2d !< Layer summed v-flux transport - real, dimension(SZI_(G),SZJ_(G),G%ke) :: tendency !< tendency array for diagn + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: tendency !< tendency array for diagnostic real, dimension(SZI_(G),SZJ_(G)) :: tendency_2d !< depth integrated content tendency for diagn - type(tracer_type), pointer :: Tracer => NULL() !< Pointer to the current tracer - integer :: remap_method !< Reconstruction method - integer :: i,j,k,m !< indices to loop over + type(tracer_type), pointer :: tracer => NULL() !< Pointer to the current tracer + real, dimension(SZK_(GV)) :: tracer_1d !< 1d-array used to remap tracer change to native grid + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: tracer_old !< local copy of the initial tracer concentration, + !! only used to compute tendencies. + real, dimension(SZI_(G),SZJ_(G)) :: tracer_int, tracer_end + !< integrated tracer in the native grid, before and after + ! LBD is applied. + integer :: i, j, k, m !< indices to loop over real :: Idt !< inverse of the time step [s-1] + real :: tmp1, tmp2 !< temporary variables + call cpu_clock_begin(id_clock_lbd) Idt = 1./dt - hbl(:,:) = 0. if (ASSOCIATED(CS%KPP_CSp)) call KPP_get_BLD(CS%KPP_CSp, hbl, G, US, m_to_BLD_units=GV%m_to_H) - if (ASSOCIATED(CS%energetic_PBL_CSp)) & - call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, hbl, G, US, m_to_MLD_units=GV%m_to_H) - + if (ASSOCIATED(CS%energetic_PBL_CSp)) call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, hbl, G, US, & + m_to_MLD_units=GV%m_to_H) call pass_var(hbl,G%Domain) do m = 1,Reg%ntr + ! current tracer tracer => Reg%tr(m) + if (CS%debug) then + call hchksum(tracer%t, "before LBD "//tracer%name,G%HI) + endif + ! for diagnostics - if (tracer%id_lbdxy_conc > 0 .or. tracer%id_lbdxy_cont > 0 .or. tracer%id_lbdxy_cont_2d > 0) then + if (tracer%id_lbdxy_conc > 0 .or. tracer%id_lbdxy_cont > 0 .or. tracer%id_lbdxy_cont_2d > 0 .or. CS%debug) then tendency(:,:,:) = 0.0 + tracer_old(:,:,:) = tracer%t(:,:,:) endif - ! Interpolate state to interface - do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1 - call build_reconstructions_1d( CS%remap_CS, G%ke, h(i,j,:), tracer%t(i,j,:), ppoly0_coefs(i,j,:,:), & - ppoly0_E(i,j,:,:), ppoly_S, remap_method, GV%H_subroundoff, GV%H_subroundoff) - enddo ; enddo - - ! Diffusive fluxes in the i-direction + ! Diffusive fluxes in the i- and j-direction uFlx(:,:,:) = 0. vFlx(:,:,:) = 0. - uFlx_bulk(:,:) = 0. - vFlx_bulk(:,:) = 0. - - ! Method #1 (layer by layer) - if (CS%method == 1) then - do j=G%jsc,G%jec - do i=G%isc-1,G%iec - if (G%mask2dCu(I,j)>0.) then - call fluxes_layer_method(SURFACE, GV%ke, CS%deg, h(I,j,:), h(I+1,j,:), hbl(I,j), hbl(I+1,j), & - G%areaT(I,j), G%areaT(I+1,j), tracer%t(I,j,:), tracer%t(I+1,j,:), ppoly0_coefs(I,j,:,:), & - ppoly0_coefs(I+1,j,:,:), ppoly0_E(I,j,:,:), ppoly0_E(I+1,j,:,:), remap_method, Coef_x(I,j), & - uFlx(I,j,:), CS%linear) - endif - enddo - enddo - do J=G%jsc-1,G%jec - do i=G%isc,G%iec - if (G%mask2dCv(i,J)>0.) then - call fluxes_layer_method(SURFACE, GV%ke, CS%deg, h(i,J,:), h(i,J+1,:), hbl(i,J), hbl(i,J+1), & - G%areaT(i,J), G%areaT(i,J+1), tracer%t(i,J,:), tracer%t(i,J+1,:), ppoly0_coefs(i,J,:,:), & - ppoly0_coefs(i,J+1,:,:), ppoly0_E(i,J,:,:), ppoly0_E(i,J+1,:,:), remap_method, Coef_y(i,J), & - vFlx(i,J,:), CS%linear) - endif - enddo - enddo - ! Method #2 (bulk approach) - elseif (CS%method == 2) then - do j=G%jsc,G%jec - do i=G%isc-1,G%iec - if (G%mask2dCu(I,j)>0.) then - call fluxes_bulk_method(SURFACE, GV%ke, CS%deg, h(I,j,:), h(I+1,j,:), hbl(I,j), hbl(I+1,j), & - G%areaT(I,j), G%areaT(I+1,j), tracer%t(I,j,:), tracer%t(I+1,j,:), & - ppoly0_coefs(I,j,:,:), ppoly0_coefs(I+1,j,:,:), ppoly0_E(I,j,:,:), & - ppoly0_E(I+1,j,:,:), remap_method, Coef_x(I,j), uFlx_bulk(I,j), uFlx(I,j,:), CS%limiter, & - CS%linear) - endif - enddo + ! LBD layer by layer + do j=G%jsc,G%jec + do i=G%isc-1,G%iec + if (G%mask2dCu(I,j)>0.) then + call fluxes_layer_method(SURFACE, G%ke, hbl(I,j), hbl(I+1,j), & + h(I,j,:), h(I+1,j,:), tracer%t(I,j,:), tracer%t(I+1,j,:), & + Coef_x(I,j), uFlx(I,j,:), G%areaT(I,j), G%areaT(I+1,j), CS) + endif enddo - do J=G%jsc-1,G%jec - do i=G%isc,G%iec - if (G%mask2dCv(i,J)>0.) then - call fluxes_bulk_method(SURFACE, GV%ke, CS%deg, h(i,J,:), h(i,J+1,:), hbl(i,J), hbl(i,J+1), & - G%areaT(i,J), G%areaT(i,J+1), tracer%t(i,J,:), tracer%t(i,J+1,:), & - ppoly0_coefs(i,J,:,:), ppoly0_coefs(i,J+1,:,:), ppoly0_E(i,J,:,:), & - ppoly0_E(i,J+1,:,:), remap_method, Coef_y(i,J), vFlx_bulk(i,J), vFlx(i,J,:), CS%limiter, & - CS%linear) - endif - enddo + enddo + do J=G%jsc-1,G%jec + do i=G%isc,G%iec + if (G%mask2dCv(i,J)>0.) then + call fluxes_layer_method(SURFACE, GV%ke, hbl(i,J), hbl(i,J+1), & + h(i,J,:), h(i,J+1,:), tracer%t(i,J,:), tracer%t(i,J+1,:), & + Coef_y(i,J), vFlx(i,J,:), G%areaT(i,J), G%areaT(i,J+1), CS) + endif enddo - ! Post tracer bulk diags - if (tracer%id_lbd_bulk_dfx>0) call post_data(tracer%id_lbd_bulk_dfx, uFlx_bulk*Idt, CS%diag) - if (tracer%id_lbd_bulk_dfy>0) call post_data(tracer%id_lbd_bulk_dfy, vFlx_bulk*Idt, CS%diag) - endif + enddo ! Update the tracer fluxes do k=1,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec if (G%mask2dT(i,j)>0.) then tracer%t(i,j,k) = tracer%t(i,j,k) + (( (uFlx(I-1,j,k)-uFlx(I,j,k)) ) + ( (vFlx(i,J-1,k)-vFlx(i,J,k) ) ))* & - (G%IareaT(i,j)/( h(i,j,k) + GV%H_subroundoff)) - + G%IareaT(i,j) / ( h(i,j,k) + GV%H_subroundoff ) if (tracer%id_lbdxy_conc > 0 .or. tracer%id_lbdxy_cont > 0 .or. tracer%id_lbdxy_cont_2d > 0 ) then - tendency(i,j,k) = (( (uFlx(I-1,j,k)-uFlx(I,j,k)) ) + ( (vFlx(i,J-1,k)-vFlx(i,J,k) ) )) * G%IareaT(i,j) * Idt + tendency(i,j,k) = (tracer%t(i,j,k)-tracer_old(i,j,k)) * Idt endif - endif enddo ; enddo ; enddo + if (CS%debug) then + call hchksum(tracer%t, "after LBD "//tracer%name,G%HI) + tracer_int(:,:) = 0.0; tracer_end(:,:) = 0.0 + ! tracer (native grid) before and after LBD + do j=G%jsc,G%jec ; do i=G%isc,G%iec + do k=1,GV%ke + tracer_int(i,j) = tracer_int(i,j) + tracer_old(i,j,k) * & + (h(i,j,k)*(G%mask2dT(i,j)*G%areaT(i,j))) + tracer_end(i,j) = tracer_end(i,j) + tracer%t(i,j,k) * & + (h(i,j,k)*(G%mask2dT(i,j)*G%areaT(i,j))) + enddo + enddo; enddo + + tmp1 = SUM(tracer_int) + tmp2 = SUM(tracer_end) + call sum_across_PEs(tmp1) + call sum_across_PEs(tmp2) + if (is_root_pe()) write(*,*)'Total '//tracer%name//' before/after LBD:', tmp1, tmp2 + endif + ! Post the tracer diagnostics if (tracer%id_lbd_dfx>0) call post_data(tracer%id_lbd_dfx, uFlx*Idt, CS%diag) if (tracer%id_lbd_dfy>0) call post_data(tracer%id_lbd_dfy, vFlx*Idt, CS%diag) @@ -297,66 +288,16 @@ subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) ! the tendency array and its units. if (tracer%id_lbdxy_conc > 0) then do k=1,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec - tendency(i,j,k) = tendency(i,j,k) / ( h(i,j,k) + GV%H_subroundoff ) + tendency(i,j,k) = tendency(i,j,k) / ( h(i,j,k) + CS%H_subroundoff ) enddo ; enddo ; enddo call post_data(tracer%id_lbdxy_conc, tendency, CS%diag) endif enddo -end subroutine lateral_boundary_diffusion - -!< Calculate bulk layer value of a scalar quantity as the thickness weighted average -real function bulk_average(boundary, nk, deg, h, hBLT, phi, ppoly0_E, ppoly0_coefs, method, k_top, zeta_top, k_bot, & - zeta_bot) - integer :: boundary !< SURFACE or BOTTOM [nondim] - integer :: nk !< Number of layers [nondim] - integer :: deg !< Degree of polynomial [nondim] - real, dimension(nk) :: h !< Layer thicknesses [H ~> m or kg m-2] - real :: hBLT !< Depth of the boundary layer [H ~> m or kg m-2] - real, dimension(nk) :: phi !< Scalar quantity - real, dimension(nk,2) :: ppoly0_E !< Edge value of polynomial - real, dimension(nk,deg+1) :: ppoly0_coefs!< Coefficients of polynomial - integer :: method !< Remapping scheme to use - - integer :: k_top !< Index of the first layer within the boundary - real :: zeta_top !< Fraction of the layer encompassed by the bottom boundary layer - !! (0 if none, 1. if all). For the surface, this is always 0. because - !! integration starts at the surface [nondim] - integer :: k_bot !< Index of the last layer within the boundary - real :: zeta_bot !< Fraction of the layer encompassed by the surface boundary layer - !! (0 if none, 1. if all). For the bottom boundary layer, this is always 1. - !! because integration starts at the bottom [nondim] - ! Local variables - real :: htot !< Running sum of the thicknesses (top to bottom) - integer :: k !< k indice - - - htot = 0. - bulk_average = 0. - if (hblt == 0.) return - if (boundary == SURFACE) then - htot = (h(k_bot) * zeta_bot) - bulk_average = average_value_ppoly( nk, phi, ppoly0_E, ppoly0_coefs, method, k_bot, 0., zeta_bot) * htot - do k = k_bot-1,1,-1 - bulk_average = bulk_average + phi(k)*h(k) - htot = htot + h(k) - enddo - elseif (boundary == BOTTOM) then - htot = (h(k_top) * zeta_top) - ! (note 1-zeta_top because zeta_top is the fraction of the layer) - bulk_average = average_value_ppoly( nk, phi, ppoly0_E, ppoly0_coefs, method, k_top, (1.-zeta_top), 1.) * htot - do k = k_top+1,nk - bulk_average = bulk_average + phi(k)*h(k) - htot = htot + h(k) - enddo - else - call MOM_error(FATAL, "bulk_average: a valid boundary type must be provided.") - endif + call cpu_clock_end(id_clock_lbd) - bulk_average = bulk_average / hBLT - -end function bulk_average +end subroutine lateral_boundary_diffusion !> Calculate the harmonic mean of two quantities !! See \ref section_harmonic_mean. @@ -370,6 +311,185 @@ real function harmonic_mean(h1,h2) endif end function harmonic_mean +!> Returns the location of the minimum value in a 1D array +!! between indices s and e. +integer function find_minimum(x, s, e) + integer, intent(in) :: s, e !< start and end indices + real, dimension(e), intent(in) :: x !< 1D array to be checked + + ! local variables + real :: minimum + integer :: location + integer :: i + + minimum = x(s) ! assume the first is the min + location = s ! record its position + do i = s+1, e ! start with next elements + if (x(i) < minimum) then ! if x(i) less than the min? + minimum = x(i) ! Yes, a new minimum found + location = i ! record its position + end if + enddo + find_minimum = location ! return the position +end function find_minimum + +!> Swaps the values of its two formal arguments. +subroutine swap(a, b) + real, intent(inout) :: a, b !< values to be swaped + + ! local variables + real :: tmp + tmp = a + a = b + b = tmp +end subroutine swap + +!> Receives a 1D array x and sorts it into ascending order. +subroutine sort(x, n) + real, dimension(n), intent(inout) :: x !< 1D array to be sorted + integer, intent(in ) :: n !< # of pts in the array + + ! local variables + integer :: i, location + + do i = 1, n-1 + location = find_minimum(x, i, n) ! find min from this to last + call swap(x(i), x(location)) ! swap this and the minimum + enddo +end subroutine sort + +!> Returns the unique values in a 1D array. +subroutine unique(val, n, val_unique, val_max) + integer, intent(in ) :: n !< # of pts in the array. + real, dimension(n), intent(in ) :: val !< 1D array to be checked. + real, dimension(:), allocatable, intent(inout) :: val_unique !< Returned 1D array with unique values. + real, optional, intent(in ) :: val_max !< sets the maximum value in val_unique to + !! this value. + ! local variables + real, dimension(n) :: tmp + integer :: i, j, ii + real :: min_val, max_val + logical :: limit + + limit = .false. + if (present(val_max)) then + limit = .true. + if (val_max > MAXVAL(val)) then + if (is_root_pe()) write(*,*)'val_max, MAXVAL(val)',val_max, MAXVAL(val) + call MOM_error(FATAL,"Houston, we've had a problem in unique (val_max cannot be > MAXVAL(val))") + endif + endif + + tmp(:) = 0. + min_val = MINVAL(val)-1 + max_val = MAXVAL(val) + i = 0 + do while (min_valmin_val) + tmp(i) = min_val + enddo + ii = i + if (limit) then + do j=1,ii + if (tmp(j) <= val_max) i = j + enddo + endif + allocate(val_unique(i), source=tmp(1:i)) +end subroutine unique + + +!> Given layer thicknesses (and corresponding interfaces) and BLDs in two adjacent columns, +!! return a set of 1-d layer thicknesses whose interfaces cover all interfaces in the left +!! and right columns plus the two BLDs. This can be used to accurately remap tracer tendencies +!! in both columns. +subroutine merge_interfaces(nk, h_L, h_R, hbl_L, hbl_R, H_subroundoff, h) + integer, intent(in ) :: nk !< Number of layers [nondim] + real, dimension(nk), intent(in ) :: h_L !< Layer thicknesses in the left column [H ~> m or kg m-2] + real, dimension(nk), intent(in ) :: h_R !< Layer thicknesses in the right column [H ~> m or kg m-2] + real, intent(in ) :: hbl_L !< Thickness of the boundary layer in the left column + !! [H ~> m or kg m-2] + real, intent(in ) :: hbl_R !< Thickness of the boundary layer in the right column + !! [H ~> m or kg m-2] + real, intent(in ) :: H_subroundoff !< GV%H_subroundoff [H ~> m or kg m-2] + real, dimension(:), allocatable, intent(inout) :: h !< Combined thicknesses [H ~> m or kg m-2] + + ! Local variables + integer :: n !< Number of layers in eta_all + real, dimension(nk+1) :: eta_L, eta_R!< Interfaces in the left and right coloumns + real, dimension(:), allocatable :: eta_all !< Combined interfaces in the left/right columns + hbl_L and hbl_R + real, dimension(:), allocatable :: eta_unique !< Combined interfaces (eta_L, eta_R), possibly hbl_L and hbl_R + real :: min_depth !< Minimum depth + real :: max_depth !< Maximum depth + real :: max_bld !< Deepest BLD + integer :: k, kk, nk1 !< loop indices (k and kk) and array size (nk1) + + n = (2*nk)+3 + allocate(eta_all(n)) + ! compute and merge interfaces + eta_L(:) = 0.0; eta_R(:) = 0.0; eta_all(:) = 0.0 + kk = 0 + do k=2,nk+1 + eta_L(k) = eta_L(k-1) + h_L(k-1) + eta_R(k) = eta_R(k-1) + h_R(k-1) + kk = kk + 2 + eta_all(kk) = eta_L(k) + eta_all(kk+1) = eta_R(k) + enddo + + ! add hbl_L and hbl_R into eta_all + eta_all(kk+2) = hbl_L + eta_all(kk+3) = hbl_R + + ! find maximum depth + min_depth = MIN(MAXVAL(eta_L), MAXVAL(eta_R)) + max_bld = MAX(hbl_L, hbl_R) + max_depth = MIN(min_depth, max_bld) + + ! sort eta_all + call sort(eta_all, n) + ! remove duplicates from eta_all and sets maximum depth + call unique(eta_all, n, eta_unique, max_depth) + + nk1 = SIZE(eta_unique) + allocate(h(nk1-1)) + do k=1,nk1-1 + h(k) = (eta_unique(k+1) - eta_unique(k)) + H_subroundoff + enddo +end subroutine merge_interfaces + +!> Calculates the maximum flux that can leave a cell and uses that to apply a +!! limiter to F_layer. +subroutine flux_limiter(F_layer, area_L, area_R, phi_L, phi_R, h_L, h_R) + real, intent(inout) :: F_layer !< Tracer flux to be checked [H L2 conc ~> m3 conc] + real, intent(in ) :: area_L, area_R !< Area of left and right cells [L2 ~> m2] + real, intent(in ) :: h_L, h_R !< Thickness of left and right cells [H ~> m or kg m-2] + real, intent(in ) :: phi_L, phi_R !< Tracer concentration in the left and right cells + !! [conc] + + ! local variables + real :: F_max !< maximum flux allowed + ! limit the flux to 0.2 of the tracer *gradient* + ! Why 0.2? + ! t=0 t=inf + ! 0 .2 + ! 0 1 0 .2.2.2 + ! 0 .2 + ! + F_max = -0.2 * ((area_R*(phi_R*h_R))-(area_L*(phi_L*h_L))) + + if ( SIGN(1.,F_layer) == SIGN(1., F_max)) then + ! Apply flux limiter calculated above + if (F_max >= 0.) then + F_layer = MIN(F_layer,F_max) + else + F_layer = MAX(F_layer,F_max) + endif + else + F_layer = 0.0 + endif +end subroutine flux_limiter + !> Find the k-index range corresponding to the layers that are within the boundary-layer region subroutine boundary_k_range(boundary, nk, h, hbl, k_top, zeta_top, k_bot, zeta_bot) integer, intent(in ) :: boundary !< SURFACE or BOTTOM [nondim] @@ -435,372 +555,174 @@ subroutine boundary_k_range(boundary, nk, h, hbl, k_top, zeta_top, k_bot, zeta_b end subroutine boundary_k_range - !> Calculate the lateral boundary diffusive fluxes using the layer by layer method. -!! See \ref section_method1 -subroutine fluxes_layer_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, & - ppoly0_coefs_L, ppoly0_coefs_R, ppoly0_E_L, ppoly0_E_R, method, khtr_u, & - F_layer, linear_decay) - - integer, intent(in ) :: boundary !< Which boundary layer SURFACE or BOTTOM [nondim] - integer, intent(in ) :: nk !< Number of layers [nondim] - integer, intent(in ) :: deg !< order of the polynomial reconstruction [nondim] - real, dimension(nk), intent(in ) :: h_L !< Layer thickness (left) [H ~> m or kg m-2] - real, dimension(nk), intent(in ) :: h_R !< Layer thickness (right) [H ~> m or kg m-2] - real, intent(in ) :: hbl_L !< Thickness of the boundary boundary - !! layer (left) [H ~> m or kg m-2] - real, intent(in ) :: hbl_R !< Thickness of the boundary boundary - !! layer (right) [H ~> m or kg m-2] - real, intent(in ) :: area_L !< Area of the horizontal grid (left) [L2 ~> m2] - real, intent(in ) :: area_R !< Area of the horizontal grid (right) [L2 ~> m2] - real, dimension(nk), intent(in ) :: phi_L !< Tracer values (left) [conc] - real, dimension(nk), intent(in ) :: phi_R !< Tracer values (right) [conc] - real, dimension(nk,deg+1), intent(in ) :: ppoly0_coefs_L !< Tracer reconstruction (left) [conc] - real, dimension(nk,deg+1), intent(in ) :: ppoly0_coefs_R !< Tracer reconstruction (right) [conc] - real, dimension(nk,2), intent(in ) :: ppoly0_E_L !< Polynomial edge values (left) [nondim] - real, dimension(nk,2), intent(in ) :: ppoly0_E_R !< Polynomial edge values (right) [nondim] - integer, intent(in ) :: method !< Method of polynomial integration [nondim] - real, intent(in ) :: khtr_u !< Horizontal diffusivities times delta t - !! at a velocity point [L2 ~> m2] - real, dimension(nk), intent( out) :: F_layer !< Layerwise diffusive flux at U- or V-point - !! [H L2 conc ~> m3 conc] - logical, optional, intent(in ) :: linear_decay !< If True, apply a linear transition at the base of - !! the boundary layer +!! See \ref section_method +subroutine fluxes_layer_method(boundary, ke, hbl_L, hbl_R, h_L, h_R, phi_L, phi_R, & + khtr_u, F_layer, area_L, area_R, CS) + + integer, intent(in ) :: boundary !< Which boundary layer SURFACE or BOTTOM [nondim] + integer, intent(in ) :: ke !< Number of layers in the native grid [nondim] + real, intent(in ) :: hbl_L !< Thickness of the boundary boundary + !! layer (left) [H ~> m or kg m-2] + real, intent(in ) :: hbl_R !< Thickness of the boundary boundary + !! layer (right) [H ~> m or kg m-2] + real, dimension(ke), intent(in ) :: h_L !< Thicknesses in the native grid (left) [H ~> m or kg m-2] + real, dimension(ke), intent(in ) :: h_R !< Thicknesses in the native grid (right) [H ~> m or kg m-2] + real, dimension(ke), intent(in ) :: phi_L !< Tracer values in the native grid (left) [conc] + real, dimension(ke), intent(in ) :: phi_R !< Tracer values in the native grid (right) [conc] + real, intent(in ) :: khtr_u !< Horizontal diffusivities times delta t + !! at a velocity point [L2 ~> m2] + real, dimension(ke), intent( out) :: F_layer !< Layerwise diffusive flux at U- or V-point in the native + !! grid [H L2 conc ~> m3 conc] + real, intent(in ) :: area_L !< Area of the horizontal grid (left) [L2 ~> m2] + real, intent(in ) :: area_R !< Area of the horizontal grid (right) [L2 ~> m2] + type(lbd_CS), pointer :: CS !< Lateral diffusion control structure + !! the boundary layer ! Local variables - real, dimension(nk) :: h_means !< Calculate the layer-wise harmonic means [H ~> m or kg m-2] - real :: khtr_avg !< Thickness-weighted diffusivity at the u-point [m^2 s^-1] - !! This is just to remind developers that khtr_avg should be - !! computed once khtr is 3D. - real :: heff !< Harmonic mean of layer thicknesses [H ~> m or kg m-2] - real :: inv_heff !< Inverse of the harmonic mean of layer thicknesses - !! [H-1 ~> m-1 or m2 kg-1] - real :: phi_L_avg, phi_R_avg !< Bulk, thickness-weighted tracer averages (left and right column) - !! [conc m^-3 ] - real :: htot !< Total column thickness [H ~> m or kg m-2] - real :: heff_tot !< Total effective column thickness in the transition layer [m] - integer :: k, k_bot_min, k_top_max !< k-indices, min and max for bottom and top, respectively - integer :: k_bot_max, k_top_min !< k-indices, max and min for bottom and top, respectively - integer :: k_bot_diff, k_top_diff !< different between left and right k-indices for bottom and top, respectively - integer :: k_top_L, k_bot_L !< k-indices left - integer :: k_top_R, k_bot_R !< k-indices right - real :: zeta_top_L, zeta_top_R !< distance from the top of a layer to the boundary - !! layer depth [nondim] - real :: zeta_bot_L, zeta_bot_R !< distance from the bottom of a layer to the boundary - !!layer depth [nondim] - real :: h_work_L, h_work_R !< dummy variables - real :: hbl_min !< minimum BLD (left and right) [m] - real :: wgt !< weight to be used in the linear transition to the interior [nondim] - real :: a !< coefficient to be used in the linear transition to the interior [nondim] - logical :: linear !< True if apply a linear transition + real, dimension(:), allocatable :: dz_top !< The LBD z grid to be created [L ~ m] + real, dimension(:), allocatable :: phi_L_z !< Tracer values in the ztop grid (left) [conc] + real, dimension(:), allocatable :: phi_R_z !< Tracer values in the ztop grid (right) [conc] + real, dimension(:), allocatable :: F_layer_z !< Diffusive flux at U- or V-point in the ztop grid [H L2 conc ~> m3 conc] + real, dimension(ke) :: h_vel !< Thicknesses at u- and v-points in the native grid + !! The harmonic mean is used to avoid zero values [H ~> m or kg m-2] + real :: khtr_avg !< Thickness-weighted diffusivity at the u-point [m^2 s^-1] + !! This is just to remind developers that khtr_avg should be + !! computed once khtr is 3D. + real :: htot !< Total column thickness [H ~> m or kg m-2] + integer :: k, k_bot_min, k_top_max !< k-indices, min and max for bottom and top, respectively + integer :: k_bot_max, k_top_min !< k-indices, max and min for bottom and top, respectively + integer :: k_bot_diff, k_top_diff !< different between left and right k-indices for bottom and top, respectively + integer :: k_top_L, k_bot_L !< k-indices left native grid + integer :: k_top_R, k_bot_R !< k-indices right native grid + real :: zeta_top_L, zeta_top_R !< distance from the top of a layer to the boundary + !! layer depth in the native grid [nondim] + real :: zeta_bot_L, zeta_bot_R !< distance from the bottom of a layer to the boundary + !!layer depth in the native grid [nondim] + real :: hbl_min !< minimum BLD (left and right) [m] + real :: wgt !< weight to be used in the linear transition to the interior [nondim] + real :: a !< coefficient to be used in the linear transition to the interior [nondim] + real :: tmp1, tmp2 !< dummy variables + real :: htot_max !< depth below which no fluxes should be applied + integer :: nk !< number of layers in the LBD grid F_layer(:) = 0.0 if (hbl_L == 0. .or. hbl_R == 0.) then return endif - linear = .false. - if (PRESENT(linear_decay)) then - linear = linear_decay - endif + ! Define vertical grid, dz_top + call merge_interfaces(ke, h_L(:), h_R(:), hbl_L, hbl_R, CS%H_subroundoff, dz_top) + nk = SIZE(dz_top) - ! Calculate vertical indices containing the boundary layer - call boundary_k_range(boundary, nk, h_L, hbl_L, k_top_L, zeta_top_L, k_bot_L, zeta_bot_L) - call boundary_k_range(boundary, nk, h_R, hbl_R, k_top_R, zeta_top_R, k_bot_R, zeta_bot_R) + ! allocate arrays + allocate(phi_L_z(nk)); phi_L_z(:) = 0.0 + allocate(phi_R_z(nk)); phi_R_z(:) = 0.0 + allocate(F_layer_z(nk)); F_layer_z(:) = 0.0 + + ! remap tracer to dz_top + call remapping_core_h(CS%remap_cs, ke, h_L(:), phi_L(:), nk, dz_top(:), phi_L_z(:)) + call remapping_core_h(CS%remap_cs, ke, h_R(:), phi_R(:), nk, dz_top(:), phi_R_z(:)) + + ! Calculate vertical indices containing the boundary layer in dz_top + call boundary_k_range(boundary, nk, dz_top, hbl_L, k_top_L, zeta_top_L, k_bot_L, zeta_bot_L) + call boundary_k_range(boundary, nk, dz_top, hbl_R, k_top_R, zeta_top_R, k_bot_R, zeta_bot_R) if (boundary == SURFACE) then k_bot_min = MIN(k_bot_L, k_bot_R) k_bot_max = MAX(k_bot_L, k_bot_R) k_bot_diff = (k_bot_max - k_bot_min) - ! make sure left and right k indices span same range - if (k_bot_min .ne. k_bot_L) then - k_bot_L = k_bot_min - zeta_bot_L = 1.0 - endif - if (k_bot_min .ne. k_bot_R) then - k_bot_R= k_bot_min - zeta_bot_R = 1.0 - endif - - h_work_L = (h_L(k_bot_L) * zeta_bot_L) - h_work_R = (h_R(k_bot_R) * zeta_bot_R) - - phi_L_avg = average_value_ppoly( nk, phi_L, ppoly0_E_L, ppoly0_coefs_L, method, k_bot_L, 0., zeta_bot_L) - phi_R_avg = average_value_ppoly( nk, phi_R, ppoly0_E_R, ppoly0_coefs_R, method, k_bot_R, 0., zeta_bot_R) - heff = harmonic_mean(h_work_L, h_work_R) ! tracer flux where the minimum BLD intersets layer - ! GMM, khtr_avg should be computed once khtr is 3D - if ((linear) .and. (k_bot_diff .gt. 1)) then + if ((CS%linear) .and. (k_bot_diff .gt. 1)) then ! apply linear decay at the base of hbl do k = k_bot_min,1,-1 - heff = harmonic_mean(h_L(k), h_R(k)) - F_layer(k) = -(heff * khtr_u) * (phi_R(k) - phi_L(k)) + F_layer_z(k) = -(dz_top(k) * khtr_u) * (phi_R_z(k) - phi_L_z(k)) + if (CS%limiter_remap) call flux_limiter(F_layer_z(k), area_L, area_R, phi_L_z(k), & + phi_R_z(k), dz_top(k), dz_top(k)) enddo - ! heff_total - heff_tot = 0.0 + htot = 0.0 do k = k_bot_min+1,k_bot_max, 1 - heff_tot = heff_tot + harmonic_mean(h_L(k), h_R(k)) + htot = htot + dz_top(k) enddo - a = -1.0/heff_tot - heff_tot = 0.0 + a = -1.0/htot + htot = 0. do k = k_bot_min+1,k_bot_max, 1 - heff = harmonic_mean(h_L(k), h_R(k)) - wgt = (a*(heff_tot + (heff * 0.5))) + 1.0 - F_layer(k) = -(heff * khtr_u) * (phi_R(k) - phi_L(k)) * wgt - heff_tot = heff_tot + heff + wgt = (a*(htot + (dz_top(k) * 0.5))) + 1.0 + F_layer_z(k) = -(dz_top(k) * khtr_u) * (phi_R_z(k) - phi_L_z(k)) * wgt + htot = htot + dz_top(k) + if (CS%limiter_remap) call flux_limiter(F_layer_z(k), area_L, area_R, phi_L_z(k), & + phi_R_z(k), dz_top(k), dz_top(k)) enddo else - F_layer(k_bot_min) = -(heff * khtr_u) * (phi_R_avg - phi_L_avg) - do k = k_bot_min-1,1,-1 - heff = harmonic_mean(h_L(k), h_R(k)) - F_layer(k) = -(heff * khtr_u) * (phi_R(k) - phi_L(k)) + do k = k_bot_min,1,-1 + F_layer_z(k) = -(dz_top(k) * khtr_u) * (phi_R_z(k) - phi_L_z(k)) + if (CS%limiter_remap) call flux_limiter(F_layer_z(k), area_L, area_R, phi_L_z(k), & + phi_R_z(k), dz_top(k), dz_top(k)) enddo endif endif - if (boundary == BOTTOM) then - ! TODO: GMM add option to apply linear decay - k_top_max = MAX(k_top_L, k_top_R) - ! make sure left and right k indices span same range - if (k_top_max .ne. k_top_L) then - k_top_L = k_top_max - zeta_top_L = 1.0 - endif - if (k_top_max .ne. k_top_R) then - k_top_R= k_top_max - zeta_top_R = 1.0 - endif - - h_work_L = (h_L(k_top_L) * zeta_top_L) - h_work_R = (h_R(k_top_R) * zeta_top_R) - - phi_L_avg = average_value_ppoly( nk, phi_L, ppoly0_E_L, ppoly0_coefs_L, method, k_top_L, 1.0-zeta_top_L, 1.0) - phi_R_avg = average_value_ppoly( nk, phi_R, ppoly0_E_R, ppoly0_coefs_R, method, k_top_R, 1.0-zeta_top_R, 1.0) - heff = harmonic_mean(h_work_L, h_work_R) - - ! tracer flux where the minimum BLD intersets layer - F_layer(k_top_max) = (-heff * khtr_u) * (phi_R_avg - phi_L_avg) - - do k = k_top_max+1,nk - heff = harmonic_mean(h_L(k), h_R(k)) - F_layer(k) = -(heff * khtr_u) * (phi_R(k) - phi_L(k)) - enddo - endif - -end subroutine fluxes_layer_method - -!> Apply the lateral boundary diffusive fluxes calculated from a 'bulk model' -!! See \ref section_method2 -subroutine fluxes_bulk_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, ppoly0_coefs_L, & - ppoly0_coefs_R, ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer, F_limit, & - linear_decay) - - integer, intent(in ) :: boundary !< Which boundary layer SURFACE or BOTTOM [nondim] - integer, intent(in ) :: nk !< Number of layers [nondim] - integer, intent(in ) :: deg !< order of the polynomial reconstruction [nondim] - real, dimension(nk), intent(in ) :: h_L !< Layer thickness (left) [H ~> m or kg m-2] - real, dimension(nk), intent(in ) :: h_R !< Layer thickness (right) [H ~> m or kg m-2] - real, intent(in ) :: hbl_L !< Thickness of the boundary boundary - !! layer (left) [H ~> m or kg m-2] - real, intent(in ) :: hbl_R !< Thickness of the boundary boundary - !! layer (left) [H ~> m or kg m-2] - real, intent(in ) :: area_L !< Area of the horizontal grid (left) [L2 ~> m2] - real, intent(in ) :: area_R !< Area of the horizontal grid (right) [L2 ~> m2] - real, dimension(nk), intent(in ) :: phi_L !< Tracer values (left) [conc] - real, dimension(nk), intent(in ) :: phi_R !< Tracer values (right) [conc] - real, dimension(nk,deg+1), intent(in ) :: ppoly0_coefs_L !< Tracer reconstruction (left) [conc] - real, dimension(nk,deg+1), intent(in ) :: ppoly0_coefs_R !< Tracer reconstruction (right) [conc] - real, dimension(nk,2), intent(in ) :: ppoly0_E_L !< Polynomial edge values (left) [nondim] - real, dimension(nk,2), intent(in ) :: ppoly0_E_R !< Polynomial edge values (right) [nondim] - integer, intent(in ) :: method !< Method of polynomial integration [nondim] - real, intent(in ) :: khtr_u !< Horizontal diffusivities times delta t - !! at a velocity point [L2 ~> m2] - real, intent( out) :: F_bulk !< The bulk mixed layer lateral flux - !! [H L2 conc ~> m3 conc] - real, dimension(nk), intent( out) :: F_layer !< Layerwise diffusive flux at U- or V-point - !! [H L2 conc ~> m3 conc] - logical, optional, intent(in ) :: F_limit !< If True, apply a limiter - logical, optional, intent(in ) :: linear_decay !< If True, apply a linear transition at the base of - !! the boundary layer +! TODO, boundary == BOTTOM +! if (boundary == BOTTOM) then +! ! TODO: GMM add option to apply linear decay +! k_top_max = MAX(k_top_L, k_top_R) +! ! make sure left and right k indices span same range +! if (k_top_max .ne. k_top_L) then +! k_top_L = k_top_max +! zeta_top_L = 1.0 +! endif +! if (k_top_max .ne. k_top_R) then +! k_top_R= k_top_max +! zeta_top_R = 1.0 +! endif +! +! ! tracer flux where the minimum BLD intersets layer +! F_layer(k_top_max) = (-heff * khtr_u) * (phi_R_avg - phi_L_avg) +! +! do k = k_top_max+1,nk +! F_layer_z(k) = -(heff * khtr_u) * (phi_R_z(k) - phi_L_z(k)) +! enddo +! endif + + ! thicknesses at velocity points + do k = 1,ke + h_vel(k) = harmonic_mean(h_L(k), h_R(k)) + enddo - ! Local variables - real, dimension(nk) :: h_means !< Calculate the layer-wise harmonic means [H ~> m or kg m-2] - real :: khtr_avg !< Thickness-weighted diffusivity at the u-point [m^2 s^-1] - !! This is just to remind developers that khtr_avg should be - !! computed once khtr is 3D. - real :: heff !< Harmonic mean of layer thicknesses [H ~> m or kg m-2] - real :: heff_tot !< Total effective column thickness in the transition layer [m] - real :: inv_heff !< Inverse of the harmonic mean of layer thicknesses - !! [H-1 ~> m-1 or m2 kg-1] - real :: phi_L_avg, phi_R_avg !< Bulk, thickness-weighted tracer averages (left and right column) - !! [conc m^-3 ] - real :: htot ! Total column thickness [H ~> m or kg m-2] - integer :: k, k_min, k_max !< k-indices, min and max for top and bottom, respectively - integer :: k_diff !< difference between k_max and k_min - integer :: k_top_L, k_bot_L !< k-indices left - integer :: k_top_R, k_bot_R !< k-indices right - real :: zeta_top_L, zeta_top_R !< distance from the top of a layer to the - !! boundary layer [nondim] - real :: zeta_bot_L, zeta_bot_R !< distance from the bottom of a layer to the - !! boundary layer [nondim] - real :: h_work_L, h_work_R !< dummy variables - real :: F_max !< The maximum amount of flux that can leave a - !! cell [m^3 conc] - logical :: limiter !< True if flux limiter should be applied - logical :: linear !< True if apply a linear transition - real :: hfrac !< Layer fraction wrt sum of all layers [nondim] - real :: dphi !< tracer gradient [conc m^-3] - real :: wgt !< weight to be used in the linear transition to the - !! interior [nondim] - real :: a !< coefficient to be used in the linear transition to the - !! interior [nondim] - - F_bulk = 0. - F_layer(:) = 0. - if (hbl_L == 0. .or. hbl_R == 0.) then - return - endif + ! remap flux to h_vel (native grid) + call reintegrate_column(nk, dz_top(:), F_layer_z(:), ke, h_vel(:), 0.0, F_layer(:)) - limiter = .false. - if (PRESENT(F_limit)) then - limiter = F_limit - endif - linear = .false. - if (PRESENT(linear_decay)) then - linear = linear_decay + ! used to avoid fluxes below hbl + if (CS%linear) then + htot_max = MAX(hbl_L, hbl_R) + else + htot_max = MIN(hbl_L, hbl_R) endif - ! Calculate vertical indices containing the boundary layer - call boundary_k_range(boundary, nk, h_L, hbl_L, k_top_L, zeta_top_L, k_bot_L, zeta_bot_L) - call boundary_k_range(boundary, nk, h_R, hbl_R, k_top_R, zeta_top_R, k_bot_R, zeta_bot_R) - - ! Calculate bulk averages of various quantities - phi_L_avg = bulk_average(boundary, nk, deg, h_L, hbl_L, phi_L, ppoly0_E_L, ppoly0_coefs_L, method, k_top_L, & - zeta_top_L, k_bot_L, zeta_bot_L) - phi_R_avg = bulk_average(boundary, nk, deg, h_R, hbl_R, phi_R, ppoly0_E_R, ppoly0_coefs_R, method, k_top_R, & - zeta_top_R, k_bot_R, zeta_bot_R) - ! Calculate the 'bulk' diffusive flux from the bulk averaged quantities - ! GMM, khtr_avg should be computed once khtr is 3D - heff = harmonic_mean(hbl_L, hbl_R) - F_bulk = -(khtr_u * heff) * (phi_R_avg - phi_L_avg) - ! Calculate the layerwise sum of the vertical effective thickness. This is different than the heff calculated - ! above, but is used as a way to decompose the fluxes onto the individual layers - h_means(:) = 0. - if (boundary == SURFACE) then - k_min = MIN(k_bot_L, k_bot_R) - k_max = MAX(k_bot_L, k_bot_R) - k_diff = (k_max - k_min) - if ((linear) .and. (k_diff .gt. 1)) then - do k=1,k_min - h_means(k) = harmonic_mean(h_L(k),h_R(k)) - enddo - ! heff_total - heff_tot = 0.0 - do k = k_min+1,k_max, 1 - heff_tot = heff_tot + harmonic_mean(h_L(k), h_R(k)) - enddo - - a = -1.0/heff_tot - heff_tot = 0.0 - ! fluxes will decay linearly at base of hbl - do k = k_min+1,k_max, 1 - heff = harmonic_mean(h_L(k), h_R(k)) - wgt = (a*(heff_tot + (heff * 0.5))) + 1.0 - h_means(k) = harmonic_mean(h_L(k), h_R(k)) * wgt - heff_tot = heff_tot + heff - enddo - else - ! left hand side - if (k_bot_L == k_min) then - h_work_L = h_L(k_min) * zeta_bot_L - else - h_work_L = h_L(k_min) - endif - - ! right hand side - if (k_bot_R == k_min) then - h_work_R = h_R(k_min) * zeta_bot_R - else - h_work_R = h_R(k_min) - endif - - h_means(k_min) = harmonic_mean(h_work_L,h_work_R) - - do k=1,k_min-1 - h_means(k) = harmonic_mean(h_L(k),h_R(k)) - enddo - endif - - elseif (boundary == BOTTOM) then - !TODO, GMM linear decay is not implemented here - k_max = MAX(k_top_L, k_top_R) - ! left hand side - if (k_top_L == k_max) then - h_work_L = h_L(k_max) * zeta_top_L - else - h_work_L = h_L(k_max) + tmp1 = 0.0; tmp2 = 0.0 + do k = 1,ke + ! apply flux_limiter + if (CS%limiter .and. F_layer(k) /= 0.) then + call flux_limiter(F_layer(k), area_L, area_R, phi_L(k), phi_R(k), h_L(k), h_R(k)) endif - ! right hand side - if (k_top_R == k_max) then - h_work_R = h_R(k_max) * zeta_top_R - else - h_work_R = h_R(k_max) + ! if tracer point is below htot_max, set flux to zero + if (MAX(tmp1+(h_L(k)*0.5), tmp2+(h_R(k)*0.5)) > htot_max) then + F_layer(k) = 0. endif - h_means(k_max) = harmonic_mean(h_work_L,h_work_R) + tmp1 = tmp1 + h_L(k) + tmp2 = tmp2 + h_R(k) + enddo - do k=nk,k_max+1,-1 - h_means(k) = harmonic_mean(h_L(k),h_R(k)) - enddo - endif + ! deallocated arrays + deallocate(dz_top) + deallocate(phi_L_z) + deallocate(phi_R_z) + deallocate(F_layer_z) - if ( SUM(h_means) == 0. .or. F_bulk == 0.) then - return - ! Decompose the bulk flux onto the individual layers - else - ! Initialize remaining thickness - inv_heff = 1./SUM(h_means) - do k=1,nk - if ((h_means(k) > 0.) .and. (phi_L(k) /= phi_R(k))) then - hfrac = h_means(k)*inv_heff - F_layer(k) = F_bulk * hfrac - - if (limiter) then - ! limit the flux to 0.2 of the tracer *gradient* - ! Why 0.2? - ! t=0 t=inf - ! 0 .2 - ! 0 1 0 .2.2.2 - ! 0 .2 - ! - F_max = -0.2 * ((area_R*(phi_R(k)*h_R(k)))-(area_L*(phi_L(k)*h_R(k)))) - - ! check if bulk flux (or F_layer) and F_max have same direction - if ( SIGN(1.,F_bulk) == SIGN(1., F_max)) then - ! Apply flux limiter calculated above - if (F_max >= 0.) then - F_layer(k) = MIN(F_layer(k),F_max) - else - F_layer(k) = MAX(F_layer(k),F_max) - endif - else - ! do not apply a flux on this layer - F_layer(k) = 0. - endif - else - dphi = -(phi_R(k) - phi_L(k)) - if (.not. SIGN(1.,F_bulk) == SIGN(1., dphi)) then - ! upgradient, do not apply a flux on this layer - F_layer(k) = 0. - endif - endif ! limited - endif - enddo - endif - -end subroutine fluxes_bulk_method +end subroutine fluxes_layer_method !> Unit tests for near-boundary horizontal mixing logical function near_boundary_unit_tests( verbose ) @@ -808,32 +730,33 @@ logical function near_boundary_unit_tests( verbose ) ! Local variables integer, parameter :: nk = 2 ! Number of layers - integer, parameter :: deg = 1 ! Degree of reconstruction (linear here) - integer, parameter :: method = 1 ! Method used for integrating polynomials + real, dimension(nk+1) :: eta1 ! Updated interfaces with one extra value [m] + real, dimension(:), allocatable :: h1 ! Upates layer thicknesses [m] real, dimension(nk) :: phi_L, phi_R ! Tracer values (left and right column) [ nondim m^-3 ] - real, dimension(nk) :: phi_L_avg, phi_R_avg ! Bulk, thickness-weighted tracer averages (left and right column) - real, dimension(nk,deg+1) :: phi_pp_L, phi_pp_R ! Coefficients for the linear pseudo-reconstructions - ! [ nondim m^-3 ] - - real, dimension(nk,2) :: ppoly0_E_L, ppoly0_E_R! Polynomial edge values (left and right) [concentration] real, dimension(nk) :: h_L, h_R ! Layer thickness (left and right) [m] real :: khtr_u ! Horizontal diffusivities at U-point [m^2 s^-1] real :: hbl_L, hbl_R ! Depth of the boundary layer (left and right) [m] - real :: F_bulk ! Total diffusive flux across the U point [nondim s^-1] real, dimension(nk) :: F_layer ! Diffusive flux within each layer at U-point [nondim s^-1] - real :: h_u, hblt_u ! Thickness at the u-point [m] - real :: khtr_avg ! Thickness-weighted diffusivity at the u-point [m^2 s^-1] - real :: heff ! Harmonic mean of layer thicknesses [m] - real :: inv_heff ! Inverse of the harmonic mean of layer thicknesses [m^[-1] character(len=120) :: test_name ! Title of the unit test integer :: k_top ! Index of cell containing top of boundary real :: zeta_top ! Nondimension position integer :: k_bot ! Index of cell containing bottom of boundary real :: zeta_bot ! Nondimension position - real :: area_L,area_R ! Area of grid cell [m^2] - area_L = 1.; area_R = 1. ! Set to unity for all unit tests + type(lbd_CS), pointer :: CS + + allocate(CS) + ! fill required fields in CS + CS%linear=.false. + call initialize_remapping( CS%remap_CS, 'PLM', boundary_extrapolation = .true. ,& + check_reconstruction = .true., check_remapping = .true.) + call extract_member_remapping_CS(CS%remap_CS, degree=CS%deg) + CS%H_subroundoff = 1.0E-20 + CS%debug=.false. + CS%limiter=.false. + CS%limiter_remap=.false. near_boundary_unit_tests = .false. + write(stdout,*) '==== MOM_lateral_boundary_diffusion =======================' ! Unit tests for boundary_k_range test_name = 'Surface boundary spans the entire top cell' @@ -890,215 +813,153 @@ logical function near_boundary_unit_tests( verbose ) near_boundary_unit_tests = near_boundary_unit_tests .or. & test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 2, 0.25, 2, 0., test_name, verbose) - ! All cases in this section have hbl which are equal to the column thicknesses - test_name = 'Equal hbl and same layer thicknesses (gradient from right to left)' - hbl_L = 10; hbl_R = 10 - h_L = (/5.,5./) ; h_R = (/5.,5./) - phi_L = (/0.,0./) ; phi_R = (/1.,1./) - phi_pp_L(1,1) = 0.; phi_pp_L(1,2) = 0. - phi_pp_L(2,1) = 0.; phi_pp_L(2,2) = 0. - phi_pp_R(1,1) = 1.; phi_pp_R(1,2) = 0. - phi_pp_R(2,1) = 1.; phi_pp_R(2,2) = 0. - ppoly0_E_L(1,1) = 0.; ppoly0_E_L(1,2) = 0. - ppoly0_E_L(2,1) = 0.; ppoly0_E_L(2,2) = 0. - ppoly0_E_R(1,1) = 1.; ppoly0_E_R(1,2) = 1. - ppoly0_E_R(2,1) = 1.; ppoly0_E_R(2,2) = 1. - khtr_u = 1. - ! Without limiter - call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, phi_pp_R, & - ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer) + if (.not. near_boundary_unit_tests) write(stdout,*) 'Passed boundary_k_range' + + ! unit tests for sorting array and finding unique values + test_name = 'Sorting array' + eta1 = (/1., 0., 0.1/) + call sort(eta1, nk+1) near_boundary_unit_tests = near_boundary_unit_tests .or. & - test_layer_fluxes( verbose, nk, test_name, F_layer, (/-5.0,-5.0/) ) + test_layer_fluxes( verbose, nk+1, test_name, eta1, (/0., 0.1, 1./) ) - ! same as above, but with limiter - call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, phi_pp_R, & - ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer, .true.) + test_name = 'Unique values' + call unique((/0., 1., 1., 2./), nk+2, h1) near_boundary_unit_tests = near_boundary_unit_tests .or. & - test_layer_fluxes( verbose, nk, test_name, F_layer, (/-1.0,-1.0/) ) + test_layer_fluxes( verbose, nk+1, test_name, h1, (/0., 1., 2./) ) + deallocate(h1) - test_name = 'Equal hbl and same layer thicknesses (gradient from left to right)' - hbl_L = 10.; hbl_R = 10. - h_L = (/5.,5./) ; h_R = (/5.,5./) - phi_L = (/1.,1./) ; phi_R = (/0.,0./) - phi_pp_L(1,1) = 1.; phi_pp_L(1,2) = 0. - phi_pp_L(2,1) = 1.; phi_pp_L(2,2) = 0. - phi_pp_R(1,1) = 0.; phi_pp_R(1,2) = 0. - phi_pp_R(2,1) = 0.; phi_pp_R(2,2) = 0. - ppoly0_E_L(1,1) = 1.; ppoly0_E_L(1,2) = 1. - ppoly0_E_L(2,1) = 1.; ppoly0_E_L(2,2) = 1. - ppoly0_E_R(1,1) = 0.; ppoly0_E_R(1,2) = 0. - ppoly0_E_R(2,1) = 0.; ppoly0_E_R(2,2) = 0. - khtr_u = 1. - call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, phi_pp_R,& - ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer) + test_name = 'Unique values with maximum depth' + call unique((/0., 1., 1., 2., 3./), nk+3, h1, 2.) near_boundary_unit_tests = near_boundary_unit_tests .or. & - test_layer_fluxes( verbose, nk, test_name, F_layer, (/5.0,5.0/) ) + test_layer_fluxes( verbose, nk+1, test_name, h1, (/0., 1., 2./) ) + deallocate(h1) - test_name = 'Equal hbl and same layer thicknesses (no gradient)' - hbl_L = 10; hbl_R = 10 - h_L = (/5.,5./) ; h_R = (/5.,5./) - phi_L = (/1.,1./) ; phi_R = (/1.,1./) - phi_pp_L(1,1) = 1.; phi_pp_L(1,2) = 0. - phi_pp_L(2,1) = 1.; phi_pp_L(2,2) = 0. - phi_pp_R(1,1) = 1.; phi_pp_R(1,2) = 0. - phi_pp_R(2,1) = 1.; phi_pp_R(2,2) = 0. - ppoly0_E_L(1,1) = 1.; ppoly0_E_L(1,2) = 1. - ppoly0_E_L(2,1) = 1.; ppoly0_E_L(2,2) = 1. - ppoly0_E_R(1,1) = 1.; ppoly0_E_R(1,2) = 1. - ppoly0_E_R(2,1) = 1.; ppoly0_E_R(2,2) = 1. - khtr_u = 1. - call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, phi_pp_R,& - ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer) + if (.not. near_boundary_unit_tests) write(stdout,*) 'Passed sort and unique' + + ! unit tests for merge_interfaces + test_name = 'h_L = h_R and BLD_L = BLD_R' + call merge_interfaces(nk, (/1., 2./), (/1., 2./), 1.5, 1.5, CS%H_subroundoff, h1) near_boundary_unit_tests = near_boundary_unit_tests .or. & - test_layer_fluxes( verbose, nk, test_name, F_layer, (/0.0,0.0/) ) + test_layer_fluxes( verbose, nk, test_name, h1, (/1., 0.5/) ) + deallocate(h1) - test_name = 'Equal hbl and different layer thicknesses (gradient right to left)' - hbl_L = 16.; hbl_R = 16. - h_L = (/10.,6./) ; h_R = (/6.,10./) - phi_L = (/0.,0./) ; phi_R = (/1.,1./) - phi_pp_L(1,1) = 0.; phi_pp_L(1,2) = 0. - phi_pp_L(2,1) = 0.; phi_pp_L(2,2) = 0. - phi_pp_R(1,1) = 1.; phi_pp_R(1,2) = 0. - phi_pp_R(2,1) = 1.; phi_pp_R(2,2) = 0. - ppoly0_E_L(1,1) = 0.; ppoly0_E_L(1,2) = 0. - ppoly0_E_L(2,1) = 0.; ppoly0_E_L(2,2) = 0. - ppoly0_E_R(1,1) = 1.; ppoly0_E_R(1,2) = 1. - ppoly0_E_R(2,1) = 1.; ppoly0_E_R(2,2) = 1. - khtr_u = 1. - call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, phi_pp_R,& - ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer) + test_name = 'h_L = h_R and BLD_L /= BLD_R' + call merge_interfaces(nk, (/1., 2./), (/1., 2./), 0.5, 1.5, CS%H_subroundoff, h1) near_boundary_unit_tests = near_boundary_unit_tests .or. & - test_layer_fluxes( verbose, nk, test_name, F_layer, (/-8.0,-8.0/) ) - - test_name = 'Equal hbl and same layer thicknesses (diagonal tracer values)' - hbl_L = 10.; hbl_R = 10. - h_L = (/5.,5./) ; h_R = (/5.,5./) - phi_L = (/1.,0./) ; phi_R = (/0.,1./) - phi_pp_L(1,1) = 1.; phi_pp_L(1,2) = 0. - phi_pp_L(2,1) = 0.; phi_pp_L(2,2) = 0. - phi_pp_R(1,1) = 0.; phi_pp_R(1,2) = 0. - phi_pp_R(2,1) = 1.; phi_pp_R(2,2) = 0. - ppoly0_E_L(1,1) = 1.; ppoly0_E_L(1,2) = 1. - ppoly0_E_L(2,1) = 0.; ppoly0_E_L(2,2) = 0. - ppoly0_E_R(1,1) = 0.; ppoly0_E_R(1,2) = 0. - ppoly0_E_R(2,1) = 1.; ppoly0_E_R(2,2) = 1. - khtr_u = 1. - call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, phi_pp_R,& - ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer) + test_layer_fluxes( verbose, nk+1, test_name, h1, (/0.5, 0.5, 0.5/) ) + deallocate(h1) + + test_name = 'h_L /= h_R and BLD_L = BLD_R' + call merge_interfaces(nk, (/1., 3./), (/2., 2./), 1.5, 1.5, CS%H_subroundoff, h1) near_boundary_unit_tests = near_boundary_unit_tests .or. & - test_layer_fluxes( verbose, nk, test_name, F_layer, (/0.0,0.0/) ) + test_layer_fluxes( verbose, nk, test_name, h1, (/1., 0.5/) ) + deallocate(h1) - test_name = 'Different hbl and different layer thicknesses (gradient from right to left)' - hbl_L = 12; hbl_R = 20 - h_L = (/6.,6./) ; h_R = (/10.,10./) - phi_L = (/0.,0./) ; phi_R = (/1.,1./) - phi_pp_L(1,1) = 0.; phi_pp_L(1,2) = 0. - phi_pp_L(2,1) = 0.; phi_pp_L(2,2) = 0. - phi_pp_R(1,1) = 1.; phi_pp_R(1,2) = 0. - phi_pp_R(2,1) = 1.; phi_pp_R(2,2) = 0. - ppoly0_E_L(1,1) = 0.; ppoly0_E_L(1,2) = 0. - ppoly0_E_L(2,1) = 0.; ppoly0_E_L(2,2) = 0. - ppoly0_E_R(1,1) = 1.; ppoly0_E_R(1,2) = 1. - ppoly0_E_R(2,1) = 1.; ppoly0_E_R(2,2) = 1. - khtr_u = 1. - call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, phi_pp_R,& - ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer) + test_name = 'h_L /= h_R and BLD_L /= BLD_R' + call merge_interfaces(nk, (/1., 3./), (/2., 2./), 0.5, 1.5, CS%H_subroundoff, h1) near_boundary_unit_tests = near_boundary_unit_tests .or. & - test_layer_fluxes( verbose, nk, test_name, F_layer, (/-7.5,-7.5/) ) + test_layer_fluxes( verbose, nk+1, test_name, h1, (/0.5, 0.5, 0.5/) ) + deallocate(h1) - ! Cases where hbl < column thickness (polynomial coefficients specified for pseudo-linear reconstruction) + test_name = 'Left deeper than right, h_L /= h_R and BLD_L /= BLD_R' + call merge_interfaces(nk, (/2., 3./), (/2., 2./), 1.0, 2.0, CS%H_subroundoff, h1) + near_boundary_unit_tests = near_boundary_unit_tests .or. & + test_layer_fluxes( verbose, nk, test_name, h1, (/1., 1./) ) + deallocate(h1) - test_name = 'hbl < column thickness, hbl same, constant concentration each column' - hbl_L = 2; hbl_R = 2 - h_L = (/1.,2./) ; h_R = (/1.,2./) + test_name = 'Left has zero thickness, h_L /= h_R and BLD_L = BLD_R' + call merge_interfaces(nk, (/4., 0./), (/2., 2./), 2.0, 2.0, CS%H_subroundoff, h1) + near_boundary_unit_tests = near_boundary_unit_tests .or. & + test_layer_fluxes( verbose, nk-1, test_name, h1, (/2./) ) + deallocate(h1) + + test_name = 'Left has zero thickness, h_L /= h_R and BLD_L /= BLD_R' + call merge_interfaces(nk, (/4., 0./), (/2., 2./), 1.0, 2.0, CS%H_subroundoff, h1) + near_boundary_unit_tests = near_boundary_unit_tests .or. & + test_layer_fluxes( verbose, nk, test_name, h1, (/1., 1./) ) + deallocate(h1) + + test_name = 'Right has zero thickness, h_L /= h_R and BLD_L = BLD_R' + call merge_interfaces(nk, (/2., 2./), (/0., 4./), 2.0, 2.0, CS%H_subroundoff, h1) + near_boundary_unit_tests = near_boundary_unit_tests .or. & + test_layer_fluxes( verbose, nk-1, test_name, h1, (/2./) ) + deallocate(h1) + + test_name = 'Right has zero thickness, h_L /= h_R and BLD_L /= BLD_R' + call merge_interfaces(nk, (/2., 2./), (/0., 4./), 1.0, 2.0, CS%H_subroundoff, h1) + near_boundary_unit_tests = near_boundary_unit_tests .or. & + test_layer_fluxes( verbose, nk, test_name, h1, (/1., 1./) ) + deallocate(h1) + + test_name = 'Right deeper than left, h_L /= h_R and BLD_L = BLD_R' + call merge_interfaces(nk+1, (/2., 2., 0./), (/2., 2., 1./), 4., 4., CS%H_subroundoff, h1) + near_boundary_unit_tests = near_boundary_unit_tests .or. & + test_layer_fluxes( verbose, nk, test_name, h1, (/2., 2./) ) + deallocate(h1) + + test_name = 'Right and left small values at bottom, h_L /= h_R and BLD_L = BLD_R' + call merge_interfaces(nk+2, (/2., 2., 1., 1./), (/1., 1., .5, .5/), 3., 3., CS%H_subroundoff, h1) + near_boundary_unit_tests = near_boundary_unit_tests .or. & + test_layer_fluxes( verbose, nk+2, test_name, h1, (/1., 1., .5, .5/) ) + deallocate(h1) + + if (.not. near_boundary_unit_tests) write(stdout,*) 'Passed merge interfaces' + + ! All cases in this section have hbl which are equal to the column thicknesses + test_name = 'Equal hbl and same layer thicknesses (gradient from right to left)' + hbl_L = 2.; hbl_R = 2. + h_L = (/2.,2./) ; h_R = (/2.,2./) phi_L = (/0.,0./) ; phi_R = (/1.,1./) - phi_pp_L(1,1) = 0.; phi_pp_L(1,2) = 0. - phi_pp_L(2,1) = 0.; phi_pp_L(2,2) = 0. - phi_pp_R(1,1) = 1.; phi_pp_R(1,2) = 0. - phi_pp_R(2,1) = 1.; phi_pp_R(2,2) = 0. - ppoly0_E_L(1,1) = 0.; ppoly0_E_L(1,2) = 0. - ppoly0_E_L(2,1) = 0.; ppoly0_E_L(2,2) = 0. - ppoly0_E_R(1,1) = 1.; ppoly0_E_R(1,2) = 1. - ppoly0_E_R(2,1) = 1.; ppoly0_E_R(2,2) = 1. khtr_u = 1. - call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, phi_pp_R,& - ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer) + call fluxes_layer_method(SURFACE, nk, hbl_L, hbl_R, h_L, h_R, phi_L, phi_R, & + khtr_u, F_layer, 1., 1., CS) near_boundary_unit_tests = near_boundary_unit_tests .or. & - test_layer_fluxes( verbose, nk, test_name, F_layer, (/-1.,-1./) ) + test_layer_fluxes( verbose, nk, test_name, F_layer, (/-2.0,0.0/) ) - test_name = 'hbl < column thickness, hbl same, linear profile right' - hbl_L = 2; hbl_R = 2 - h_L = (/1.,2./) ; h_R = (/1.,2./) - phi_L = (/0.,0./) ; phi_R = (/0.5,2./) - phi_pp_L(1,1) = 0.; phi_pp_L(1,2) = 0. - phi_pp_L(2,1) = 0.; phi_pp_L(2,2) = 0. - phi_pp_R(1,1) = 0.; phi_pp_R(1,2) = 1. - phi_pp_R(2,1) = 1.; phi_pp_R(2,2) = 2. - khtr_u = 1. - ppoly0_E_L(1,1) = 0.; ppoly0_E_L(1,2) = 0. - ppoly0_E_L(2,1) = 0.; ppoly0_E_L(2,2) = 0. - ppoly0_E_R(1,1) = 0.; ppoly0_E_R(1,2) = 1. - ppoly0_E_R(2,1) = 1.; ppoly0_E_R(2,2) = 3. - call fluxes_bulk_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, phi_pp_R,& - ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer) + test_name = 'Equal hbl and same layer thicknesses (gradient from left to right)' + hbl_L = 2.; hbl_R = 2. + h_L = (/2.,2./) ; h_R = (/2.,2./) + phi_L = (/2.,1./) ; phi_R = (/1.,1./) + khtr_u = 0.5 + call fluxes_layer_method(SURFACE, nk, hbl_L, hbl_R, h_L, h_R, phi_L, phi_R, & + khtr_u, F_layer, 1., 1., CS) near_boundary_unit_tests = near_boundary_unit_tests .or. & - test_layer_fluxes( verbose, nk, test_name, F_layer, (/-1.,-1./) ) + test_layer_fluxes( verbose, nk, test_name, F_layer, (/1.0,0.0/) ) test_name = 'hbl < column thickness, hbl same, linear profile right, khtr=2' hbl_L = 2; hbl_R = 2 h_L = (/1.,2./) ; h_R = (/1.,2./) phi_L = (/0.,0./) ; phi_R = (/0.5,2./) - phi_pp_L(1,1) = 0.; phi_pp_L(1,2) = 0. - phi_pp_L(2,1) = 0.; phi_pp_L(2,2) = 0. - phi_pp_R(1,1) = 0.; phi_pp_R(1,2) = 1. - phi_pp_R(2,1) = 1.; phi_pp_R(2,2) = 2. khtr_u = 2. - ppoly0_E_L(1,1) = 0.; ppoly0_E_L(1,2) = 0. - ppoly0_E_L(2,1) = 0.; ppoly0_E_L(2,2) = 0. - ppoly0_E_R(1,1) = 0.; ppoly0_E_R(1,2) = 1. - ppoly0_E_R(2,1) = 1.; ppoly0_E_R(2,2) = 3. - call fluxes_layer_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, & - phi_pp_R, ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_layer) + call fluxes_layer_method(SURFACE, nk, hbl_L, hbl_R, h_L, h_R, phi_L, phi_R, & + khtr_u, F_layer, 1., 1., CS) near_boundary_unit_tests = near_boundary_unit_tests .or. & - test_layer_fluxes( verbose, nk, test_name, F_layer, (/-1.,-3./) ) + test_layer_fluxes( verbose, nk, test_name, F_layer, (/-1.0,-4.0/) ) - ! unit tests for layer by layer method - test_name = 'Different hbl and different column thicknesses (gradient from right to left)' + test_name = 'Different hbl and different column thicknesses (zero gradient)' hbl_L = 12; hbl_R = 20 h_L = (/6.,6./) ; h_R = (/10.,10./) - phi_L = (/0.,0./) ; phi_R = (/1.,1./) - phi_pp_L(1,1) = 0.; phi_pp_L(1,2) = 0. - phi_pp_L(2,1) = 0.; phi_pp_L(2,2) = 0. - phi_pp_R(1,1) = 1.; phi_pp_R(1,2) = 0. - phi_pp_R(2,1) = 1.; phi_pp_R(2,2) = 0. - ppoly0_E_L(1,1) = 0.; ppoly0_E_L(1,2) = 0. - ppoly0_E_L(2,1) = 0.; ppoly0_E_L(2,2) = 0. - ppoly0_E_R(1,1) = 1.; ppoly0_E_R(1,2) = 1. - ppoly0_E_R(2,1) = 1.; ppoly0_E_R(2,2) = 1. + phi_L = (/1.,1./) ; phi_R = (/1.,1./) khtr_u = 1. - call fluxes_layer_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, & - phi_pp_R, ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_layer) + call fluxes_layer_method(SURFACE, nk, hbl_L, hbl_R, h_L, h_R, phi_L, phi_R, & + khtr_u, F_layer, 1., 1., CS) near_boundary_unit_tests = near_boundary_unit_tests .or. & - test_layer_fluxes( verbose, nk, test_name, F_layer, (/-7.5,-7.5/) ) - - test_name = 'Different hbl and different column thicknesses (linear profile right)' - - hbl_L = 15; hbl_R = 6 - h_L = (/10.,10./) ; h_R = (/12.,10./) - phi_L = (/0.,0./) ; phi_R = (/1.,3./) - phi_pp_L(1,1) = 0.; phi_pp_L(1,2) = 0. - phi_pp_L(2,1) = 0.; phi_pp_L(2,2) = 0. - phi_pp_R(1,1) = 0.; phi_pp_R(1,2) = 2. - phi_pp_R(2,1) = 2.; phi_pp_R(2,2) = 2. - ppoly0_E_L(1,1) = 0.; ppoly0_E_L(1,2) = 0. - ppoly0_E_L(2,1) = 0.; ppoly0_E_L(2,2) = 0. - ppoly0_E_R(1,1) = 0.; ppoly0_E_R(1,2) = 2. - ppoly0_E_R(2,1) = 2.; ppoly0_E_R(2,2) = 4. + test_layer_fluxes( verbose, nk, test_name, F_layer, (/0.,0./) ) + + test_name = 'Different hbl and different column thicknesses (gradient from left to right)' + + hbl_L = 15; hbl_R = 10. + h_L = (/10.,5./) ; h_R = (/10.,0./) + phi_L = (/1.,1./) ; phi_R = (/0.,0./) khtr_u = 1. - call fluxes_layer_method(SURFACE, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, phi_pp_L, & - phi_pp_R, ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_layer) + call fluxes_layer_method(SURFACE, nk, hbl_L, hbl_R, h_L, h_R, phi_L, phi_R, & + khtr_u, F_layer, 1., 1., CS) + near_boundary_unit_tests = near_boundary_unit_tests .or. & - test_layer_fluxes( verbose, nk, test_name, F_layer, (/-3.75,0.0/) ) + test_layer_fluxes( verbose, nk, test_name, F_layer, (/10.,0.0/) ) + + if (.not. near_boundary_unit_tests) write(stdout,*) 'Passed fluxes_layer_method' + end function near_boundary_unit_tests !> Returns true if output of near-boundary unit tests does not match correct computed values @@ -1111,16 +972,15 @@ logical function test_layer_fluxes(verbose, nk, test_name, F_calc, F_ans) real, dimension(nk), intent(in) :: F_ans !< Fluxes of the unitless tracer calculated by hand [s^-1] ! Local variables integer :: k - integer, parameter :: stdunit = stdout test_layer_fluxes = .false. do k=1,nk if ( F_calc(k) /= F_ans(k) ) then test_layer_fluxes = .true. - write(stdunit,*) "MOM_lateral_boundary_diffusion, UNIT TEST FAILED: ", test_name - write(stdunit,10) k, F_calc(k), F_ans(k) + write(stdout,*) "MOM_lateral_boundary_diffusion, UNIT TEST FAILED: ", test_name + write(stdout,10) k, F_calc(k), F_ans(k) elseif (verbose) then - write(stdunit,10) k, F_calc(k), F_ans(k) + write(stdout,10) k, F_calc(k), F_ans(k) endif enddo @@ -1141,19 +1001,17 @@ logical function test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, k_top_a character(len=80) :: test_name !< Name of the unit test logical :: verbose !< If true always print output - integer, parameter :: stdunit = stdout - test_boundary_k_range = k_top .ne. k_top_ans test_boundary_k_range = test_boundary_k_range .or. (zeta_top .ne. zeta_top_ans) test_boundary_k_range = test_boundary_k_range .or. (k_bot .ne. k_bot_ans) test_boundary_k_range = test_boundary_k_range .or. (zeta_bot .ne. zeta_bot_ans) - if (test_boundary_k_range) write(stdunit,*) "UNIT TEST FAILED: ", test_name + if (test_boundary_k_range) write(stdout,*) "UNIT TEST FAILED: ", test_name if (test_boundary_k_range .or. verbose) then - write(stdunit,20) "k_top", k_top, "k_top_ans", k_top_ans - write(stdunit,20) "k_bot", k_bot, "k_bot_ans", k_bot_ans - write(stdunit,30) "zeta_top", zeta_top, "zeta_top_ans", zeta_top_ans - write(stdunit,30) "zeta_bot", zeta_bot, "zeta_bot_ans", zeta_bot_ans + write(stdout,20) "k_top", k_top, "k_top_ans", k_top_ans + write(stdout,20) "k_bot", k_bot, "k_bot_ans", k_bot_ans + write(stdout,30) "zeta_top", zeta_top, "zeta_top_ans", zeta_top_ans + write(stdout,30) "zeta_bot", zeta_bot, "zeta_bot_ans", zeta_bot_ans endif 20 format(A,"=",i3,X,A,"=",i3) @@ -1169,24 +1027,28 @@ end function test_boundary_k_range !! The LBD framework accounts for the effects of diabatic mesoscale fluxes !! within surface and bottom boundary layers. Unlike the equivalent adiabatic !! fluxes, which is applied along neutral density surfaces, LBD is purely -!! horizontal. +!! horizontal. To assure that diffusive fluxes are strictly horizontal +!! regardless of the vertical coordinate system, this method relies on +!! regridding/remapping techniques. !! -!! The bottom boundary layer fluxes remain to be implemented, although most +!! The bottom boundary layer fluxes remain to be implemented, although some !! of the steps needed to do so have already been added and tested. !! -!! Boundary lateral diffusion can be applied using one of the three methods: +!! Boundary lateral diffusion is applied as follows: !! -!! * [Method #1: Along layer](@ref section_method2) (default); -!! * [Method #2: Bulk layer](@ref section_method1); +!! 1) remap tracer to a z* grid (LBD grid) +!! 2) calculate diffusive tracer fluxes (F) in the LBD grid using a layer by layer approach (@ref section_method) +!! 3) remap fluxes to the native grid +!! 4) update tracer by adding the divergence of F !! -!! A brief summary of these methods is provided below. +!! \subsection section_method Along layer approach !! -!! \subsection section_method1 Along layer approach (Method #1) +!! Here diffusion is applied layer by layer using only information from neighboring cells. !! -!! This is the recommended and more straight forward method where diffusion is -!! applied layer by layer using only information from neighboring cells. +!! Step #1: define vertical grid using interfaces and surface boundary layers from left and right +!! columns (see merge_interfaces). !! -!! Step #1: compute vertical indices containing boundary layer (boundary_k_range). +!! Step #2: compute vertical indices containing boundary layer (boundary_k_range). !! For the TOP boundary layer, these are: !! !! k_top, k_bot, zeta_top, zeta_bot @@ -1195,9 +1057,7 @@ end function test_boundary_k_range !! !! \f[ F_{k} = -KHTR \times h_{eff}(k) \times (\phi_R(k) - \phi_L(k)), \f] !! where h_eff is the [harmonic mean](@ref section_harmonic_mean) of the layer thickness -!! in the left and right columns. This method does not require a limiter since KHTR -!! is already limted based on a diffusive CFL condition prior to the call of this -!! module. +!! in the left and right columns. !! !! Step #3: option to linearly decay the flux from k_bot_min to k_bot_max: !! @@ -1206,44 +1066,9 @@ end function test_boundary_k_range !! layer depth (k_bot_min) and the lower interface of the layer containing the !! maximum layer depth (k_bot_max). !! -!! \subsection section_method2 Bulk layer approach (Method #2) -!! -!! Apply the lateral boundary diffusive fluxes calculated from a 'bulk model'.This -!! is a lower order representation (Kraus-Turner like approach) which assumes that -!! eddies are acting along well mixed layers (i.e., eddies do not know care about -!! vertical tracer gradients within the boundary layer). -!! -!! Step #1: compute vertical indices containing boundary layer (boundary_k_range). -!! For the TOP boundary layer, these are: -!! -!! k_top, k_bot, zeta_top, zeta_bot -!! -!! Step #2: compute bulk averages (thickness weighted) tracer averages (phi_L and phi_R), -!! then calculate the bulk diffusive flux (F_{bulk}): -!! -!! \f[ F_{bulk} = -KHTR \times h_{eff} \times (\phi_R - \phi_L), \f] -!! where h_eff is the [harmonic mean](@ref section_harmonic_mean) of the boundary layer depth -!! in the left and right columns (\f[ HBL_L \f] and \f[ HBL_R \f], respectively). -!! -!! Step #3: decompose F_bulk onto individual layers: -!! -!! \f[ F_{layer}(k) = F_{bulk} \times h_{frac}(k) , \f] -!! -!! where h_{frac} is -!! -!! \f[ h_{frac}(k) = h_u(k) \times \frac{1}{\sum(h_u)}. \f] -!! -!! h_u is the [harmonic mean](@ref section_harmonic_mean) of thicknesses at each layer. -!! Special care (layer reconstruction) must be taken at k_min = min(k_botL, k_bot_R). -!! -!! Step #4: option to linearly decay the flux from k_bot_min to k_bot_max: -!! -!! If LBD_LINEAR_TRANSITION = True and k_bot_diff > 1, the diffusive flux will decay -!! linearly between the top interface of the layer containing the minimum boundary -!! layer depth (k_bot_min) and the lower interface of the layer containing the -!! maximum layer depth (k_bot_max). -!! -!! Step #5: limit the tracer flux so that 1) only down-gradient fluxes are applied, +!! Step #4: remap the fluxes back to the native grid. This is done at velocity points, whose vertical grid +!! is determined using [harmonic mean](@ref section_harmonic_mean). To assure monotonicity, +!! tracer fluxes are limited so that 1) only down-gradient fluxes are applied, !! and 2) the flux cannot be larger than F_max, which is defined using the tracer !! gradient: !! diff --git a/src/tracer/MOM_neutral_diffusion.F90 b/src/tracer/MOM_neutral_diffusion.F90 index 086caf390f..e995ca1972 100644 --- a/src/tracer/MOM_neutral_diffusion.F90 +++ b/src/tracer/MOM_neutral_diffusion.F90 @@ -28,8 +28,7 @@ module MOM_neutral_diffusion use MOM_energetic_PBL, only : energetic_PBL_get_MLD, energetic_PBL_CS use MOM_diabatic_driver, only : diabatic_CS, extract_diabatic_member use MOM_lateral_boundary_diffusion, only : boundary_k_range, SURFACE, BOTTOM - -use iso_fortran_env, only : stdout=>output_unit, stderr=>error_unit +use MOM_io, only : stdout, stderr implicit none ; private @@ -228,9 +227,6 @@ logical function neutral_diffusion_init(Time, G, US, param_file, diag, EOS, diab default = .true.) endif - ! Store a rescaling factor for use in diagnostic messages. - CS%R_to_kg_m3 = US%R_to_kg_m3 - if (CS%interior_only) then call extract_diabatic_member(diabatic_CSp, KPP_CSp=CS%KPP_CSp) call extract_diabatic_member(diabatic_CSp, energetic_PBL_CSp=CS%energetic_PBL_CSp) @@ -238,6 +234,8 @@ logical function neutral_diffusion_init(Time, G, US, param_file, diag, EOS, diab call MOM_error(FATAL,"NDIFF_INTERIOR_ONLY is true, but no valid boundary layer scheme was found") endif endif + ! Store a rescaling factor for use in diagnostic messages. + CS%R_to_kg_m3 = US%R_to_kg_m3 ! call get_param(param_file, mdl, "KHTR", CS%KhTr, & ! "The background along-isopycnal tracer diffusivity.", & @@ -316,11 +314,10 @@ subroutine neutral_diffusion_calc_coeffs(G, GV, US, h, T, S, CS, p_surf) ! Check if hbl needs to be extracted if (CS%interior_only) then - hbl(:,:) = 0. if (ASSOCIATED(CS%KPP_CSp)) call KPP_get_BLD(CS%KPP_CSp, hbl, G, US, m_to_BLD_units=GV%m_to_H) - if (ASSOCIATED(CS%energetic_PBL_CSp)) & - call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, hbl, G, US, m_to_MLD_units=GV%m_to_H) - call pass_var(hbl, G%Domain) + if (ASSOCIATED(CS%energetic_PBL_CSp)) call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, hbl, G, US, & + m_to_MLD_units=GV%m_to_H) + call pass_var(hbl,G%Domain) ! get k-indices and zeta do j=G%jsc-1, G%jec+1 ; do i=G%isc-1,G%iec+1 call boundary_k_range(SURFACE, G%ke, h(i,j,:), hbl(i,j), k_top(i,j), zeta_top(i,j), k_bot(i,j), zeta_bot(i,j)) @@ -435,7 +432,7 @@ subroutine neutral_diffusion_calc_coeffs(G, GV, US, h, T, S, CS, p_surf) if (CS%interior_only) then if (.not. CS%stable_cell(i,j,k_bot(i,j))) zeta_bot(i,j) = -1. ! set values in the surface and bottom boundary layer to false. - do k = 1, k_bot(i,j)-1 + do k = 1, k_bot(i,j) CS%stable_cell(i,j,k) = .false. enddo endif @@ -481,7 +478,7 @@ subroutine neutral_diffusion_calc_coeffs(G, GV, US, h, T, S, CS, p_surf) call find_neutral_surface_positions_continuous(G%ke, & CS%Pint(i,j,:), CS%Tint(i,j,:), CS%Sint(i,j,:), CS%dRdT(i,j,:), CS%dRdS(i,j,:), & CS%Pint(i,j+1,:), CS%Tint(i,j+1,:), CS%Sint(i,j+1,:), CS%dRdT(i,j+1,:), CS%dRdS(i,j+1,:), & - CS%vPoL(i,J,:), CS%vPoR(i,J,:), CS%vKoL(i,J,:), CS%vKoR(i,J,:), CS%vhEff(i,J,:), & + CS%vPoL(i,J,:), CS%vPoR(i,J,:), CS%vKoL(i,J,:), CS%vKoR(i,J,:), CS%vhEff(i,J,:), & k_bot(i,J), k_bot(i,J+1), zeta_bot(i,J), zeta_bot(i,J+1)) else call find_neutral_surface_positions_discontinuous(CS, G%ke, & diff --git a/src/tracer/MOM_tracer_hor_diff.F90 b/src/tracer/MOM_tracer_hor_diff.F90 index 43ede7cff5..d3afed2f75 100644 --- a/src/tracer/MOM_tracer_hor_diff.F90 +++ b/src/tracer/MOM_tracer_hor_diff.F90 @@ -23,7 +23,7 @@ module MOM_tracer_hor_diff use MOM_neutral_diffusion, only : neutral_diffusion_init, neutral_diffusion_end use MOM_neutral_diffusion, only : neutral_diffusion_CS use MOM_neutral_diffusion, only : neutral_diffusion_calc_coeffs, neutral_diffusion -use MOM_lateral_boundary_diffusion, only : lateral_boundary_diffusion_CS, lateral_boundary_diffusion_init +use MOM_lateral_boundary_diffusion, only : lbd_CS, lateral_boundary_diffusion_init use MOM_lateral_boundary_diffusion, only : lateral_boundary_diffusion use MOM_tracer_registry, only : tracer_registry_type, tracer_type, MOM_tracer_chksum use MOM_unit_scaling, only : unit_scale_type @@ -64,7 +64,7 @@ module MOM_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. - type(lateral_boundary_diffusion_CS), pointer :: lateral_boundary_diffusion_CSp => NULL() !< Control structure for + type(lbd_CS), pointer :: lateral_boundary_diffusion_CSp => NULL() !< Control structure for !! lateral boundary mixing. type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to !! regulate the timing of diagnostic output. @@ -1430,9 +1430,10 @@ end subroutine tracer_epipycnal_ML_diff !> Initialize lateral tracer diffusion module -subroutine tracer_hor_diff_init(Time, G, US, param_file, diag, EOS, diabatic_CSp, CS) +subroutine tracer_hor_diff_init(Time, G, GV, US, param_file, diag, EOS, diabatic_CSp, CS) type(time_type), target, intent(in) :: Time !< current model time type(ocean_grid_type), intent(in) :: G !< ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(diag_ctrl), target, intent(inout) :: diag !< diagnostic control type(EOS_type), target, intent(in) :: EOS !< Equation of state CS @@ -1511,7 +1512,7 @@ subroutine tracer_hor_diff_init(Time, G, US, param_file, diag, EOS, diabatic_CSp diabatic_CSp, CS%neutral_diffusion_CSp ) if (CS%use_neutral_diffusion .and. CS%Diffuse_ML_interior) call MOM_error(FATAL, "MOM_tracer_hor_diff: "// & "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%use_lateral_boundary_diffusion = lateral_boundary_diffusion_init(Time, G, GV, 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: "// & "USE_LATERAL_BOUNDARY_DIFFUSION and DIFFUSE_ML_TO_INTERIOR are mutually exclusive!")