diff --git a/.testing/Makefile b/.testing/Makefile index 21da6cfde4..02f6557c09 100644 --- a/.testing/Makefile +++ b/.testing/Makefile @@ -84,6 +84,9 @@ FCFLAGS_COVERAGE ?= # - FMS cannot be built with the same aggressive initialization flags as MOM6, # so FCFLAGS_INIT is used to provide additional MOM6 configuration. +# User-defined LDFLAGS (applied to all builds and FMS) +LDFLAGS_USER ?= + # Set to `true` to require identical results from DEBUG and REPRO builds # NOTE: Many compilers (Intel, GCC on ARM64) do not yet produce identical # results across DEBUG and REPRO builds (as defined below), so we disable on @@ -217,8 +220,8 @@ REPRO_FCFLAGS := FCFLAGS="$(FCFLAGS_REPRO) $(FCFLAGS_FMS)" OPENMP_FCFLAGS := FCFLAGS="$(FCFLAGS_DEBUG) $(FCFLAGS_INIT) $(FCFLAGS_FMS)" TARGET_FCFLAGS := FCFLAGS="$(FCFLAGS_DEBUG) $(FCFLAGS_INIT) $(FCFLAGS_FMS)" -MOM_LDFLAGS := LDFLAGS="$(LDFLAGS_FMS)" -SYMMETRIC_LDFLAGS := LDFLAGS="$(COVERAGE) $(LDFLAGS_FMS)" +MOM_LDFLAGS := LDFLAGS="$(LDFLAGS_FMS) $(LDFLAGS_USER)" +SYMMETRIC_LDFLAGS := LDFLAGS="$(COVERAGE) $(LDFLAGS_FMS) $(LDFLAGS_USER)" # Environment variable configuration @@ -286,7 +289,7 @@ $(TARGET_CODEBASE)/ac/configure: $(TARGET_CODEBASE) $(TARGET_CODEBASE): git clone --recursive $(MOM_TARGET_URL) $@ - cd $@ && git checkout $(MOM_TARGET_BRANCH) + cd $@ && git checkout --recurse-submodules $(MOM_TARGET_BRANCH) #--- diff --git a/config_src/drivers/ice_solo_driver/ice_shelf_driver.F90 b/config_src/drivers/ice_solo_driver/ice_shelf_driver.F90 index bd64050a6f..959e4676d0 100644 --- a/config_src/drivers/ice_solo_driver/ice_shelf_driver.F90 +++ b/config_src/drivers/ice_solo_driver/ice_shelf_driver.F90 @@ -185,6 +185,8 @@ program Shelf_main endif endif + call Get_MOM_Input(param_file, dirs) + ! Read ocean_solo restart, which can override settings from the namelist. if (file_exists(trim(dirs%restart_input_dir)//'ice_solo.res')) then call open_ASCII_file(unit, trim(dirs%restart_input_dir)//'ice_solo.res', action=READONLY_FILE) @@ -215,7 +217,6 @@ program Shelf_main Start_time = real_to_time(0.0) endif - call Get_MOM_Input(param_file, dirs) ! Determining the internal unit scaling factors for this run. call unit_scaling_init(param_file, US) diff --git a/config_src/drivers/mct_cap/mom_surface_forcing_mct.F90 b/config_src/drivers/mct_cap/mom_surface_forcing_mct.F90 index e5e1309a91..e675170575 100644 --- a/config_src/drivers/mct_cap/mom_surface_forcing_mct.F90 +++ b/config_src/drivers/mct_cap/mom_surface_forcing_mct.F90 @@ -496,8 +496,8 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, if (associated(fluxes%frunoff)) then fluxes%latent(i,j) = fluxes%latent(i,j) - & IOB%rofi_flux(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion - fluxes%latent_frunoff_diag(i,j) = - G%mask2dT(i,j) * & - IOB%rofi_flux(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion + fluxes%latent_frunoff_diag(i,j) = -G%mask2dT(i,j) & + * IOB%rofi_flux(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion endif ! contribution from evaporation if (associated(IOB%q_flux)) then diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90 index 2f6dfa86d3..4776183c90 100644 --- a/config_src/drivers/nuopc_cap/mom_cap.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap.F90 @@ -53,7 +53,7 @@ module MOM_cap_mod use ESMF, only: ESMF_Grid, ESMF_GridCreate, ESMF_GridAddCoord use ESMF, only: ESMF_GridGetCoord, ESMF_GridAddItem, ESMF_GridGetItem use ESMF, only: ESMF_GridComp, ESMF_GridCompSetEntryPoint, ESMF_GridCompGet -use ESMF, only: ESMF_LogFoundError, ESMF_LogWrite, ESMF_LogSetError +use ESMF, only: ESMF_LogWrite, ESMF_LogSetError use ESMF, only: ESMF_LOGERR_PASSTHRU, ESMF_KIND_R8, ESMF_RC_VAL_WRONG use ESMF, only: ESMF_GEOMTYPE_MESH, ESMF_GEOMTYPE_GRID, ESMF_SUCCESS use ESMF, only: ESMF_METHOD_INITIALIZE, ESMF_MethodRemove, ESMF_State @@ -79,6 +79,7 @@ module MOM_cap_mod use ESMF, only: ESMF_AlarmCreate, ESMF_ClockGetAlarmList, ESMF_AlarmList_Flag use ESMF, only: ESMF_AlarmGet, ESMF_AlarmIsCreated, ESMF_ALARMLIST_ALL, ESMF_AlarmIsEnabled use ESMF, only: ESMF_STATEITEM_NOTFOUND, ESMF_FieldWrite +use ESMF, only: ESMF_END_ABORT, ESMF_Finalize use ESMF, only: operator(==), operator(/=), operator(+), operator(-) ! TODO ESMF_GridCompGetInternalState does not have an explicit Fortran interface. @@ -140,6 +141,7 @@ module MOM_cap_mod logical :: profile_memory = .true. logical :: grid_attach_area = .false. logical :: use_coldstart = .true. +logical :: use_mommesh = .true. character(len=128) :: scalar_field_name = '' integer :: scalar_field_count = 0 integer :: scalar_field_idx_grid_nx = 0 @@ -153,7 +155,7 @@ module MOM_cap_mod type(ESMF_GeomType_Flag) :: geomtype = ESMF_GEOMTYPE_MESH #else logical :: cesm_coupled = .false. -type(ESMF_GeomType_Flag) :: geomtype = ESMF_GEOMTYPE_GRID +type(ESMF_GeomType_Flag) :: geomtype #endif character(len=8) :: restart_mode = 'alarms' @@ -353,6 +355,25 @@ subroutine InitializeP0(gcomp, importState, exportState, clock, rc) write(logmsg,*) use_coldstart call ESMF_LogWrite('MOM_cap:use_coldstart = '//trim(logmsg), ESMF_LOGMSG_INFO) + use_mommesh = .true. + call NUOPC_CompAttributeGet(gcomp, name="use_mommesh", value=value, & + isPresent=isPresent, isSet=isSet, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (isPresent .and. isSet) use_mommesh=(trim(value)=="true") + write(logmsg,*) use_mommesh + call ESMF_LogWrite('MOM_cap:use_mommesh = '//trim(logmsg), ESMF_LOGMSG_INFO) + + if(use_mommesh)then + geomtype = ESMF_GEOMTYPE_MESH + call NUOPC_CompAttributeGet(gcomp, name='mesh_ocn', isPresent=isPresent, isSet=isSet, rc=rc) + if (.not. isPresent .and. .not. isSet) then + call ESMF_LogWrite('geomtype set to mesh but mesh_ocn is not specified', ESMF_LOGMSG_INFO) + call ESMF_Finalize(endflag=ESMF_END_ABORT) + endif + else + geomtype = ESMF_GEOMTYPE_GRID + endif + end subroutine !> Called by NUOPC to advertise import and export fields. "Advertise" @@ -443,18 +464,11 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) !--------------------------------- call ESMF_VMGet(vm, pet=localPet, peCount=localPeCount, rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & - line=__LINE__, & - file=__FILE__)) & - return ! bail out - + if (ChkErr(rc,__LINE__,u_FILE_u)) return if(localPeCount == 1) then call NUOPC_CompAttributeGet(gcomp, "nthreads", value=cvalue, rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & - line=__LINE__, & - file=__FILE__)) & - return ! bail out + if (ChkErr(rc,__LINE__,u_FILE_u)) return read(cvalue,*) nthrds else nthrds = localPeCount @@ -879,18 +893,11 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) !--------------------------------- call ESMF_VMGet(vm, pet=localPet, peCount=localPeCount, rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & - line=__LINE__, & - file=__FILE__)) & - return ! bail out - + if (ChkErr(rc,__LINE__,u_FILE_u)) return if(localPeCount == 1) then call NUOPC_CompAttributeGet(gcomp, "nthreads", value=cvalue, rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & - line=__LINE__, & - file=__FILE__)) & - return ! bail out + if (ChkErr(rc,__LINE__,u_FILE_u)) return read(cvalue,*) nthrds else nthrds = localPeCount @@ -1054,7 +1061,7 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) ! Determine mesh areas for regridding call ESMF_MeshGet(Emesh, numOwnedElements=numOwnedElements, spatialDim=spatialDim, rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + if (ChkErr(rc,__LINE__,u_FILE_u)) return allocate (mod2med_areacor(numOwnedElements)) allocate (med2mod_areacor(numOwnedElements)) @@ -1064,11 +1071,11 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) #ifdef CESMCOUPLED ! Determine model areas and flux correction factors (module variables in mom_) call ESMF_StateGet(exportState, itemName=trim(fldsFrOcn(2)%stdname), field=lfield, rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + if (ChkErr(rc,__LINE__,u_FILE_u)) return call ESMF_FieldRegridGetArea(lfield, rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + if (ChkErr(rc,__LINE__,u_FILE_u)) return call ESMF_FieldGet(lfield, farrayPtr=dataPtr_mesh_areas, rc=rc) - if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return + if (ChkErr(rc,__LINE__,u_FILE_u)) return allocate(mesh_areas(numOwnedElements)) allocate(model_areas(numOwnedElements)) diff --git a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 index d65f6be4e1..7168823fbc 100644 --- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 @@ -796,7 +796,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS) endif forces%ustar(i,j) = sqrt(gustiness*Irho0 + Irho0*tau_mag) enddo ; enddo - + call pass_vector(forces%taux, forces%tauy, G%Domain, halo=1) elseif (wind_stagger == AGRID) then call pass_vector(taux_at_h, tauy_at_h, G%Domain, To_All+Omit_Corners, stagger=AGRID, halo=1) @@ -822,7 +822,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS) forces%ustar(i,j) = sqrt(gustiness*Irho0 + Irho0 * G%mask2dT(i,j) * & sqrt(taux_at_h(i,j)**2 + tauy_at_h(i,j)**2)) enddo ; enddo - + call pass_vector(forces%taux, forces%tauy, G%Domain, halo=1) else ! C-grid wind stresses. if (G%symmetric) & call fill_symmetric_edges(forces%taux, forces%tauy, G%Domain) diff --git a/config_src/drivers/solo_driver/MOM_driver.F90 b/config_src/drivers/solo_driver/MOM_driver.F90 index 9c222bb0bb..8edad7fa05 100644 --- a/config_src/drivers/solo_driver/MOM_driver.F90 +++ b/config_src/drivers/solo_driver/MOM_driver.F90 @@ -250,6 +250,10 @@ program MOM_main ! This call sets the number and affinity of threads with openMP. !$ call set_MOM_thread_affinity(ocean_nthreads, use_hyper_thread) + ! This call is required to initiate dirs%restart_input_dir for ocean_solo.res + ! The contents of dirs will be reread in initialize_MOM. + call get_MOM_input(dirs=dirs) + ! Read ocean_solo restart, which can override settings from the namelist. if (file_exists(trim(dirs%restart_input_dir)//'ocean_solo.res')) then call open_ASCII_file(unit, trim(dirs%restart_input_dir)//'ocean_solo.res', action=READONLY_FILE) diff --git a/config_src/infra/FMS1/MOM_diag_manager_infra.F90 b/config_src/infra/FMS1/MOM_diag_manager_infra.F90 index 702c464814..18c80cf24c 100644 --- a/config_src/infra/FMS1/MOM_diag_manager_infra.F90 +++ b/config_src/infra/FMS1/MOM_diag_manager_infra.F90 @@ -20,9 +20,9 @@ module MOM_diag_manager_infra use diag_manager_mod, only : register_diag_field_fms => register_diag_field use diag_manager_mod, only : register_static_field_fms => register_static_field use diag_manager_mod, only : get_diag_field_id_fms => get_diag_field_id -use time_manager_mod, only : time_type +use MOM_time_manager, only : time_type use MOM_domain_infra, only : MOM_domain_type -use MOM_error_handler, only : MOM_error, FATAL, WARNING +use MOM_error_infra, only : MOM_error => MOM_err, FATAL, WARNING implicit none ; private diff --git a/config_src/infra/FMS1/MOM_domain_infra.F90 b/config_src/infra/FMS1/MOM_domain_infra.F90 index 86e85e60a6..fc39777a2f 100644 --- a/config_src/infra/FMS1/MOM_domain_infra.F90 +++ b/config_src/infra/FMS1/MOM_domain_infra.F90 @@ -4,7 +4,7 @@ module MOM_domain_infra ! This file is part of MOM6. See LICENSE.md for the license. use MOM_coms_infra, only : PE_here, root_PE, num_PEs -use MOM_cpu_clock, only : cpu_clock_begin, cpu_clock_end +use MOM_cpu_clock_infra, only : cpu_clock_begin, cpu_clock_end use MOM_error_infra, only : MOM_error=>MOM_err, NOTE, WARNING, FATAL use mpp_domains_mod, only : domain2D, domain1D @@ -1689,6 +1689,8 @@ subroutine clone_MD_to_d2D(MD_in, mpp_domain, min_halo, halo_size, symmetric, & if ((MD_in%io_layout(1) + MD_in%io_layout(2) > 0) .and. & (MD_in%layout(1)*MD_in%layout(2) > 1)) then call mpp_define_io_domain(mpp_domain, MD_in%io_layout) + else + call mpp_define_io_domain(mpp_domain, (/ 1, 1 /) ) endif end subroutine clone_MD_to_d2D diff --git a/pkg/CVMix-src b/pkg/CVMix-src index 534fc38e75..9423197f89 160000 --- a/pkg/CVMix-src +++ b/pkg/CVMix-src @@ -1 +1 @@ -Subproject commit 534fc38e759fcb8a2563fa0dc4c0bbf81f758db9 +Subproject commit 9423197f894112edfcb1502245f7d7b873d551f9 diff --git a/src/core/MOM_CoriolisAdv.F90 b/src/core/MOM_CoriolisAdv.F90 index 8dab711d32..7cbc1eb262 100644 --- a/src/core/MOM_CoriolisAdv.F90 +++ b/src/core/MOM_CoriolisAdv.F90 @@ -76,8 +76,10 @@ module MOM_CoriolisAdv integer :: id_rvxu = -1, id_rvxv = -1 ! integer :: id_hf_gKEu = -1, id_hf_gKEv = -1 integer :: id_hf_gKEu_2d = -1, id_hf_gKEv_2d = -1 + integer :: id_intz_gKEu_2d = -1, id_intz_gKEv_2d = -1 ! integer :: id_hf_rvxu = -1, id_hf_rvxv = -1 integer :: id_hf_rvxu_2d = -1, id_hf_rvxv_2d = -1 + integer :: id_intz_rvxu_2d = -1, id_intz_rvxv_2d = -1 !>@} end type CoriolisAdv_CS @@ -219,12 +221,17 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS) real, allocatable, dimension(:,:) :: & hf_gKEu_2d, hf_gKEv_2d, & ! Depth sum of hf_gKEu, hf_gKEv [L T-2 ~> m s-2]. hf_rvxu_2d, hf_rvxv_2d ! Depth sum of hf_rvxu, hf_rvxv [L T-2 ~> m s-2]. + !real, allocatable, dimension(:,:,:) :: & ! hf_gKEu, hf_gKEv, & ! accel. due to KE gradient x fract. thickness [L T-2 ~> m s-2]. ! hf_rvxu, hf_rvxv ! accel. due to RV x fract. thickness [L T-2 ~> m s-2]. ! 3D diagnostics hf_gKEu etc. are commented because there is no clarity on proper remapping grid option. ! The code is retained for degugging purposes in the future. +! Diagnostics for depth-integrated momentum budget terms + real, dimension(SZIB_(G),SZJ_(G)) :: intz_gKEu_2d, intz_rvxv_2d ! [L2 T-2 ~> m2 s-2]. + real, dimension(SZI_(G),SZJB_(G)) :: intz_gKEv_2d, intz_rvxu_2d ! [L2 T-2 ~> m2 s-2]. + ! To work, the following fields must be set outside of the usual ! is to ie range before this subroutine is called: ! v(is-1:ie+2,js-1:je+1), u(is-1:ie+1,js-1:je+2), h(is-1:ie+2,js-1:je+2), @@ -883,6 +890,22 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS) deallocate(hf_gKEv_2d) endif + if (CS%id_intz_gKEu_2d > 0) then + intz_gKEu_2d(:,:) = 0.0 + do k=1,nz ; do j=js,je ; do I=Isq,Ieq + intz_gKEu_2d(I,j) = intz_gKEu_2d(I,j) + AD%gradKEu(I,j,k) * AD%diag_hu(I,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_intz_gKEu_2d, intz_gKEu_2d, CS%diag) + endif + + if (CS%id_intz_gKEv_2d > 0) then + intz_gKEv_2d(:,:) = 0.0 + do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + intz_gKEv_2d(i,J) = intz_gKEv_2d(i,J) + AD%gradKEv(i,J,k) * AD%diag_hv(i,J,k) + enddo ; enddo ; enddo + call post_data(CS%id_intz_gKEv_2d, intz_gKEv_2d, CS%diag) + endif + !if (CS%id_hf_rvxv > 0) then ! allocate(hf_rvxv(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke)) ! do k=1,nz ; do j=js,je ; do I=Isq,Ieq @@ -918,6 +941,22 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS) call post_data(CS%id_hf_rvxu_2d, hf_rvxu_2d, CS%diag) deallocate(hf_rvxu_2d) endif + + if (CS%id_intz_rvxv_2d > 0) then + intz_rvxv_2d(:,:) = 0.0 + do k=1,nz ; do j=js,je ; do I=Isq,Ieq + intz_rvxv_2d(I,j) = intz_rvxv_2d(I,j) + AD%rv_x_v(I,j,k) * AD%diag_hu(I,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_intz_rvxv_2d, intz_rvxv_2d, CS%diag) + endif + + if (CS%id_intz_rvxu_2d > 0) then + intz_rvxu_2d(:,:) = 0.0 + do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + intz_rvxu_2d(i,J) = intz_rvxu_2d(i,J) + AD%rv_x_u(i,J,k) * AD%diag_hv(i,J,k) + enddo ; enddo ; enddo + call post_data(CS%id_intz_rvxu_2d, intz_rvxu_2d, CS%diag) + endif endif end subroutine CorAdCalc @@ -1163,24 +1202,24 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS) 'Potential Vorticity', 'm-1 s-1', conversion=GV%m_to_H*US%s_to_T) CS%id_gKEu = register_diag_field('ocean_model', 'gKEu', diag%axesCuL, Time, & - 'Zonal Acceleration from Grad. Kinetic Energy', 'm-1 s-2', conversion=US%L_T2_to_m_s2) + 'Zonal Acceleration from Grad. Kinetic Energy', 'm s-2', conversion=US%L_T2_to_m_s2) if (CS%id_gKEu > 0) call safe_alloc_ptr(AD%gradKEu,IsdB,IedB,jsd,jed,nz) CS%id_gKEv = register_diag_field('ocean_model', 'gKEv', diag%axesCvL, Time, & - 'Meridional Acceleration from Grad. Kinetic Energy', 'm-1 s-2', conversion=US%L_T2_to_m_s2) + 'Meridional Acceleration from Grad. Kinetic Energy', 'm s-2', conversion=US%L_T2_to_m_s2) if (CS%id_gKEv > 0) call safe_alloc_ptr(AD%gradKEv,isd,ied,JsdB,JedB,nz) CS%id_rvxu = register_diag_field('ocean_model', 'rvxu', diag%axesCvL, Time, & - 'Meridional Acceleration from Relative Vorticity', 'm-1 s-2', conversion=US%L_T2_to_m_s2) + 'Meridional Acceleration from Relative Vorticity', 'm s-2', conversion=US%L_T2_to_m_s2) if (CS%id_rvxu > 0) call safe_alloc_ptr(AD%rv_x_u,isd,ied,JsdB,JedB,nz) CS%id_rvxv = register_diag_field('ocean_model', 'rvxv', diag%axesCuL, Time, & - 'Zonal Acceleration from Relative Vorticity', 'm-1 s-2', conversion=US%L_T2_to_m_s2) + 'Zonal Acceleration from Relative Vorticity', 'm s-2', conversion=US%L_T2_to_m_s2) if (CS%id_rvxv > 0) call safe_alloc_ptr(AD%rv_x_v,IsdB,IedB,jsd,jed,nz) !CS%id_hf_gKEu = register_diag_field('ocean_model', 'hf_gKEu', diag%axesCuL, Time, & ! 'Fractional Thickness-weighted Zonal Acceleration from Grad. Kinetic Energy', & - ! 'm-1 s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2) + ! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2) !if (CS%id_hf_gKEu > 0) then ! call safe_alloc_ptr(AD%gradKEu,IsdB,IedB,jsd,jed,nz) ! call safe_alloc_ptr(AD%diag_hfrac_u,IsdB,IedB,jsd,jed,nz) @@ -1188,15 +1227,15 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS) !CS%id_hf_gKEv = register_diag_field('ocean_model', 'hf_gKEv', diag%axesCvL, Time, & ! 'Fractional Thickness-weighted Meridional Acceleration from Grad. Kinetic Energy', & - ! 'm-1 s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2) + ! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2) !if (CS%id_hf_gKEv > 0) then ! call safe_alloc_ptr(AD%gradKEv,isd,ied,JsdB,JedB,nz) - ! call safe_alloc_ptr(AD%diag_hfrac_v,isd,ied,Jsd,JedB,nz) + ! call safe_alloc_ptr(AD%diag_hfrac_v,isd,ied,JsdB,JedB,nz) !endif CS%id_hf_gKEu_2d = register_diag_field('ocean_model', 'hf_gKEu_2d', diag%axesCu1, Time, & 'Depth-sum Fractional Thickness-weighted Zonal Acceleration from Grad. Kinetic Energy', & - 'm-1 s-2', conversion=US%L_T2_to_m_s2) + 'm s-2', conversion=US%L_T2_to_m_s2) if (CS%id_hf_gKEu_2d > 0) then call safe_alloc_ptr(AD%gradKEu,IsdB,IedB,jsd,jed,nz) call safe_alloc_ptr(AD%diag_hfrac_u,IsdB,IedB,jsd,jed,nz) @@ -1204,10 +1243,26 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS) CS%id_hf_gKEv_2d = register_diag_field('ocean_model', 'hf_gKEv_2d', diag%axesCv1, Time, & 'Depth-sum Fractional Thickness-weighted Meridional Acceleration from Grad. Kinetic Energy', & - 'm-1 s-2', conversion=US%L_T2_to_m_s2) + 'm s-2', conversion=US%L_T2_to_m_s2) if (CS%id_hf_gKEv_2d > 0) then call safe_alloc_ptr(AD%gradKEv,isd,ied,JsdB,JedB,nz) - call safe_alloc_ptr(AD%diag_hfrac_v,isd,ied,Jsd,JedB,nz) + call safe_alloc_ptr(AD%diag_hfrac_v,isd,ied,JsdB,JedB,nz) + endif + + CS%id_intz_gKEu_2d = register_diag_field('ocean_model', 'intz_gKEu_2d', diag%axesCu1, Time, & + 'Depth-integral of Zonal Acceleration from Grad. Kinetic Energy', & + 'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2) + if (CS%id_intz_gKEu_2d > 0) then + call safe_alloc_ptr(AD%gradKEu,IsdB,IedB,jsd,jed,nz) + call safe_alloc_ptr(AD%diag_hu,IsdB,IedB,jsd,jed,nz) + endif + + CS%id_intz_gKEv_2d = register_diag_field('ocean_model', 'intz_gKEv_2d', diag%axesCv1, Time, & + 'Depth-integral of Meridional Acceleration from Grad. Kinetic Energy', & + 'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2) + if (CS%id_intz_gKEv_2d > 0) then + call safe_alloc_ptr(AD%gradKEv,isd,ied,JsdB,JedB,nz) + call safe_alloc_ptr(AD%diag_hv,isd,ied,JsdB,JedB,nz) endif !CS%id_hf_rvxu = register_diag_field('ocean_model', 'hf_rvxu', diag%axesCvL, Time, & @@ -1215,7 +1270,7 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS) ! 'm-1 s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2) !if (CS%id_hf_rvxu > 0) then ! call safe_alloc_ptr(AD%rv_x_u,isd,ied,JsdB,JedB,nz) - ! call safe_alloc_ptr(AD%diag_hfrac_v,isd,ied,Jsd,JedB,nz) + ! call safe_alloc_ptr(AD%diag_hfrac_v,isd,ied,JsdB,JedB,nz) !endif !CS%id_hf_rvxv = register_diag_field('ocean_model', 'hf_rvxv', diag%axesCuL, Time, & @@ -1227,21 +1282,37 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS) !endif CS%id_hf_rvxu_2d = register_diag_field('ocean_model', 'hf_rvxu_2d', diag%axesCv1, Time, & - 'Fractional Thickness-weighted Meridional Acceleration from Relative Vorticity', & - 'm-1 s-2', conversion=US%L_T2_to_m_s2) + 'Depth-sum Fractional Thickness-weighted Meridional Acceleration from Relative Vorticity', & + 'm s-2', conversion=US%L_T2_to_m_s2) if (CS%id_hf_rvxu_2d > 0) then call safe_alloc_ptr(AD%rv_x_u,isd,ied,JsdB,JedB,nz) - call safe_alloc_ptr(AD%diag_hfrac_v,isd,ied,Jsd,JedB,nz) + call safe_alloc_ptr(AD%diag_hfrac_v,isd,ied,JsdB,JedB,nz) endif CS%id_hf_rvxv_2d = register_diag_field('ocean_model', 'hf_rvxv_2d', diag%axesCu1, Time, & 'Depth-sum Fractional Thickness-weighted Zonal Acceleration from Relative Vorticity', & - 'm-1 s-2', conversion=US%L_T2_to_m_s2) + 'm s-2', conversion=US%L_T2_to_m_s2) if (CS%id_hf_rvxv_2d > 0) then call safe_alloc_ptr(AD%rv_x_v,IsdB,IedB,jsd,jed,nz) call safe_alloc_ptr(AD%diag_hfrac_u,IsdB,IedB,jsd,jed,nz) endif + CS%id_intz_rvxu_2d = register_diag_field('ocean_model', 'intz_rvxu_2d', diag%axesCv1, Time, & + 'Depth-integral of Meridional Acceleration from Relative Vorticity', & + 'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2) + if (CS%id_intz_rvxu_2d > 0) then + call safe_alloc_ptr(AD%rv_x_u,isd,ied,JsdB,JedB,nz) + call safe_alloc_ptr(AD%diag_hv,isd,ied,JsdB,JedB,nz) + endif + + CS%id_intz_rvxv_2d = register_diag_field('ocean_model', 'intz_rvxv_2d', diag%axesCu1, Time, & + 'Depth-integral of Fractional Thickness-weighted Zonal Acceleration from Relative Vorticity', & + 'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2) + if (CS%id_intz_rvxv_2d > 0) then + call safe_alloc_ptr(AD%rv_x_v,IsdB,IedB,jsd,jed,nz) + call safe_alloc_ptr(AD%diag_hu,IsdB,IedB,jsd,jed,nz) + endif + end subroutine CoriolisAdv_init !> Destructor for coriolisadv_cs diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 628f8e1c39..17e7ebee40 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -955,6 +955,17 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, Datv(i,J) = 0.0 ; bt_rem_v(i,J) = 0.0 ; vhbt0(i,J) = 0.0 enddo ; enddo + if (CS%linear_wave_drag) then + !$OMP parallel do default(shared) + do j=CS%jsdw,CS%jedw ; do I=CS%isdw-1,CS%iedw + Rayleigh_u(I,j) = 0.0 + enddo ; enddo + !$OMP parallel do default(shared) + do J=CS%jsdw-1,CS%jedw ; do i=CS%isdw,CS%iedw + Rayleigh_v(i,J) = 0.0 + enddo ; enddo + endif + ! Copy input arrays into their wide-halo counterparts. if (interp_eta_PF) then !$OMP parallel do default(shared) @@ -1900,20 +1911,28 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, endif ; enddo ; enddo !$OMP end do nowait endif + !$OMP do schedule(static) do J=jsv-1,jev ; do i=isv-1,iev+1 vel_prev = vbt(i,J) vbt(i,J) = bt_rem_v(i,J) * (vbt(i,J) + & - dtbt * ((BT_force_v(i,J) + Cor_v(i,J)) + PFv(i,J))) + dtbt * ((BT_force_v(i,J) + Cor_v(i,J)) + PFv(i,J))) + if (abs(vbt(i,J)) < CS%vel_underflow) vbt(i,J) = 0.0 vbt_trans(i,J) = trans_wt1*vbt(i,J) + trans_wt2*vel_prev + enddo ; enddo - if (CS%linear_wave_drag) then + if (CS%linear_wave_drag) then + !$OMP do schedule(static) + do J=jsv-1,jev ; do i=isv-1,iev+1 v_accel_bt(i,J) = v_accel_bt(i,J) + wt_accel(n) * & ((Cor_v(i,J) + PFv(i,J)) - vbt(i,J)*Rayleigh_v(i,J)) - else + enddo ; enddo + else + !$OMP do schedule(static) + do J=jsv-1,jev ; do i=isv-1,iev+1 v_accel_bt(i,J) = v_accel_bt(i,J) + wt_accel(n) * (Cor_v(i,J) + PFv(i,J)) - endif - enddo ; enddo + enddo ; enddo + endif if (integral_BT_cont) then !$OMP do schedule(static) @@ -1969,6 +1988,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, endif ; enddo ; enddo !$OMP end do nowait endif + !$OMP do schedule(static) do j=jsv,jev ; do I=isv-1,iev vel_prev = ubt(I,j) @@ -1976,15 +1996,23 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, dtbt * ((BT_force_u(I,j) + Cor_u(I,j)) + PFu(I,j))) if (abs(ubt(I,j)) < CS%vel_underflow) ubt(I,j) = 0.0 ubt_trans(I,j) = trans_wt1*ubt(I,j) + trans_wt2*vel_prev + enddo ; enddo + !$OMP end do nowait - if (CS%linear_wave_drag) then + if (CS%linear_wave_drag) then + !$OMP do schedule(static) + do j=jsv,jev ; do I=isv-1,iev u_accel_bt(I,j) = u_accel_bt(I,j) + wt_accel(n) * & - ((Cor_u(I,j) + PFu(I,j)) - ubt(I,j)*Rayleigh_u(I,j)) - else + ((Cor_u(I,j) + PFu(I,j)) - ubt(I,j)*Rayleigh_u(I,j)) + enddo ; enddo + !$OMP end do nowait + else + !$OMP do schedule(static) + do j=jsv,jev ; do I=isv-1,iev u_accel_bt(I,j) = u_accel_bt(I,j) + wt_accel(n) * (Cor_u(I,j) + PFu(I,j)) - endif - enddo ; enddo - !$OMP end do nowait + enddo ; enddo + !$OMP end do nowait + endif if (integral_BT_cont) then !$OMP do schedule(static) @@ -2046,14 +2074,20 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, dtbt * ((BT_force_u(I,j) + Cor_u(I,j)) + PFu(I,j))) if (abs(ubt(I,j)) < CS%vel_underflow) ubt(I,j) = 0.0 ubt_trans(I,j) = trans_wt1*ubt(I,j) + trans_wt2*vel_prev + enddo ; enddo - if (CS%linear_wave_drag) then + if (CS%linear_wave_drag) then + !$OMP do schedule(static) + do j=jsv-1,jev+1 ; do I=isv-1,iev u_accel_bt(I,j) = u_accel_bt(I,j) + wt_accel(n) * & ((Cor_u(I,j) + PFu(I,j)) - ubt(I,j)*Rayleigh_u(I,j)) - else + enddo ; enddo + else + !$OMP do schedule(static) + do j=jsv-1,jev+1 ; do I=isv-1,iev u_accel_bt(I,j) = u_accel_bt(I,j) + wt_accel(n) * (Cor_u(I,j) + PFu(I,j)) - endif - enddo ; enddo + enddo ; enddo + endif if (integral_BT_cont) then !$OMP do schedule(static) @@ -2128,15 +2162,24 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, dtbt * ((BT_force_v(i,J) + Cor_v(i,J)) + PFv(i,J))) if (abs(vbt(i,J)) < CS%vel_underflow) vbt(i,J) = 0.0 vbt_trans(i,J) = trans_wt1*vbt(i,J) + trans_wt2*vel_prev + enddo ; enddo + !$OMP end do nowait - if (CS%linear_wave_drag) then + if (CS%linear_wave_drag) then + !$OMP do schedule(static) + do J=jsv-1,jev ; do i=isv,iev v_accel_bt(i,J) = v_accel_bt(i,J) + wt_accel(n) * & - ((Cor_v(i,J) + PFv(i,J)) - vbt(i,J)*Rayleigh_v(i,J)) - else + ((Cor_v(i,J) + PFv(i,J)) - vbt(i,J)*Rayleigh_v(i,J)) + enddo ; enddo + !$OMP end do nowait + else + !$OMP do schedule(static) + do J=jsv-1,jev ; do i=isv,iev v_accel_bt(i,J) = v_accel_bt(i,J) + wt_accel(n) * (Cor_v(i,J) + PFv(i,J)) - endif - enddo ; enddo - !$OMP end do nowait + enddo ; enddo + !$OMP end do nowait + endif + if (integral_BT_cont) then !$OMP do schedule(static) do J=jsv-1,jev ; do i=isv,iev @@ -2632,6 +2675,17 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, enddo ; enddo ; enddo endif + if ((present(ADp)) .and. (present(BT_cont)) .and. (associated(ADp%diag_hu))) then + do k=1,nz ; do j=js,je ; do I=is-1,ie + ADp%diag_hu(I,j,k) = BT_cont%h_u(I,j,k) + enddo ; enddo ; enddo + endif + if ((present(ADp)) .and. (present(BT_cont)) .and. (associated(ADp%diag_hv))) then + do k=1,nz ; do J=js-1,je ; do i=is,ie + ADp%diag_hv(i,J,k) = BT_cont%h_v(i,J,k) + enddo ; enddo ; enddo + endif + if (G%nonblocking_updates) then if (find_etaav) call complete_group_pass(CS%pass_etaav, G%Domain) call complete_group_pass(CS%pass_ubta_uhbta, G%Domain) diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index a89514cde4..ef7da5c291 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -161,14 +161,17 @@ module MOM_dynamics_split_RK2 integer :: id_CAu = -1, id_CAv = -1 ! integer :: id_hf_PFu = -1, id_hf_PFv = -1 integer :: id_hf_PFu_2d = -1, id_hf_PFv_2d = -1 + integer :: id_intz_PFu_2d = -1, id_intz_PFv_2d = -1 ! integer :: id_hf_CAu = -1, id_hf_CAv = -1 integer :: id_hf_CAu_2d = -1, id_hf_CAv_2d = -1 + integer :: id_intz_CAu_2d = -1, id_intz_CAv_2d = -1 ! Split scheme only. integer :: id_uav = -1, id_vav = -1 integer :: id_u_BT_accel = -1, id_v_BT_accel = -1 ! integer :: id_hf_u_BT_accel = -1, id_hf_v_BT_accel = -1 integer :: id_hf_u_BT_accel_2d = -1, id_hf_v_BT_accel_2d = -1 + integer :: id_intz_u_BT_accel_2d = -1, id_intz_v_BT_accel_2d = -1 !>@} type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the @@ -336,6 +339,12 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s hf_CAu_2d, hf_CAv_2d, & ! Depth integeral of hf_CAu, hf_CAv [L T-2 ~> m s-2]. hf_u_BT_accel_2d, hf_v_BT_accel_2d ! Depth integeral of hf_u_BT_accel, hf_v_BT_accel + real, dimension(SZIB_(G),SZJ_(G)) :: & + intz_PFu_2d, intz_CAu_2d, intz_u_BT_accel_2d ! [L2 T-2 ~> m2 s-2]. + + real, dimension(SZI_(G),SZJB_(G)) :: & + intz_PFv_2d, intz_CAv_2d, intz_v_BT_accel_2d ! [L2 T-2 ~> m2 s-2]. + real :: dt_pred ! The time step for the predictor part of the baroclinic time stepping [T ~> s]. logical :: dyn_p_surf @@ -897,6 +906,21 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s ! enddo ; enddo ; enddo ! call post_data(CS%id_hf_PFv, hf_PFv, CS%diag) !endif + if (CS%id_intz_PFu_2d > 0) then + intz_PFu_2d(:,:) = 0.0 + do k=1,nz ; do j=js,je ; do I=Isq,Ieq + intz_PFu_2d(I,j) = intz_PFu_2d(I,j) + CS%PFu(I,j,k) * CS%ADp%diag_hu(I,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_intz_PFu_2d, intz_PFu_2d, CS%diag) + endif + if (CS%id_intz_PFv_2d > 0) then + intz_PFv_2d(:,:) = 0.0 + do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + intz_PFv_2d(i,J) = intz_PFv_2d(i,J) + CS%PFv(i,J,k) * CS%ADp%diag_hv(i,J,k) + enddo ; enddo ; enddo + call post_data(CS%id_intz_PFv_2d, intz_PFv_2d, CS%diag) + endif + if (CS%id_hf_PFu_2d > 0) then allocate(hf_PFu_2d(G%IsdB:G%IedB,G%jsd:G%jed)) hf_PFu_2d(:,:) = 0.0 @@ -930,6 +954,21 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s ! enddo ; enddo ; enddo ! call post_data(CS%id_hf_CAv, hf_CAv, CS%diag) !endif + if (CS%id_intz_CAu_2d > 0) then + intz_CAu_2d(:,:) = 0.0 + do k=1,nz ; do j=js,je ; do I=Isq,Ieq + intz_CAu_2d(I,j) = intz_CAu_2d(I,j) + CS%CAu(I,j,k) * CS%ADp%diag_hu(I,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_intz_CAu_2d, intz_CAu_2d, CS%diag) + endif + if (CS%id_intz_CAv_2d > 0) then + intz_CAv_2d(:,:) = 0.0 + do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + intz_CAv_2d(i,J) = intz_CAv_2d(i,J) + CS%CAv(i,J,k) * CS%ADp%diag_hv(i,J,k) + enddo ; enddo ; enddo + call post_data(CS%id_intz_CAv_2d, intz_CAv_2d, CS%diag) + endif + if (CS%id_hf_CAu_2d > 0) then allocate(hf_CAu_2d(G%IsdB:G%IedB,G%jsd:G%jed)) hf_CAu_2d(:,:) = 0.0 @@ -963,6 +1002,21 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s ! enddo ; enddo ; enddo ! call post_data(CS%id_hf_v_BT_accel, hf_v_BT_accel, CS%diag) !endif + if (CS%id_intz_u_BT_accel_2d > 0) then + intz_u_BT_accel_2d(:,:) = 0.0 + do k=1,nz ; do j=js,je ; do I=Isq,Ieq + intz_u_BT_accel_2d(I,j) = intz_u_BT_accel_2d(I,j) + CS%u_accel_bt(I,j,k) * CS%ADp%diag_hu(I,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_intz_u_BT_accel_2d, intz_u_BT_accel_2d, CS%diag) + endif + if (CS%id_intz_v_BT_accel_2d > 0) then + intz_v_BT_accel_2d(:,:) = 0.0 + do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + intz_v_BT_accel_2d(i,J) = intz_v_BT_accel_2d(i,J) + CS%v_accel_bt(i,J,k) * CS%ADp%diag_hv(i,J,k) + enddo ; enddo ; enddo + call post_data(CS%id_intz_v_BT_accel_2d, intz_v_BT_accel_2d, CS%diag) + endif + if (CS%id_hf_u_BT_accel_2d > 0) then allocate(hf_u_BT_accel_2d(G%IsdB:G%IedB,G%jsd:G%jed)) hf_u_BT_accel_2d(:,:) = 0.0 @@ -1366,7 +1420,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param !CS%id_hf_PFv = register_diag_field('ocean_model', 'hf_PFv', diag%axesCvL, Time, & ! 'Fractional Thickness-weighted Meridional Pressure Force Acceleration', 'm s-2', & ! v_extensive=.true., conversion=US%L_T2_to_m_s2) - !if(CS%id_hf_PFv > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz) + !if(CS%id_hf_PFv > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz) !CS%id_hf_CAu = register_diag_field('ocean_model', 'hf_CAu', diag%axesCuL, Time, & ! 'Fractional Thickness-weighted Zonal Coriolis and Advective Acceleration', 'm s-2', & @@ -1376,7 +1430,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param !CS%id_hf_CAv = register_diag_field('ocean_model', 'hf_CAv', diag%axesCvL, Time, & ! 'Fractional Thickness-weighted Meridional Coriolis and Advective Acceleration', 'm s-2', & ! v_extensive=.true., conversion=US%L_T2_to_m_s2) - !if(CS%id_hf_CAv > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz) + !if(CS%id_hf_CAv > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz) CS%id_hf_PFu_2d = register_diag_field('ocean_model', 'hf_PFu_2d', diag%axesCu1, Time, & 'Depth-sum Fractional Thickness-weighted Zonal Pressure Force Acceleration', 'm s-2', & @@ -1386,7 +1440,17 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param CS%id_hf_PFv_2d = register_diag_field('ocean_model', 'hf_PFv_2d', diag%axesCv1, Time, & 'Depth-sum Fractional Thickness-weighted Meridional Pressure Force Acceleration', 'm s-2', & conversion=US%L_T2_to_m_s2) - if(CS%id_hf_PFv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz) + if(CS%id_hf_PFv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz) + + CS%id_intz_PFu_2d = register_diag_field('ocean_model', 'intz_PFu_2d', diag%axesCu1, Time, & + 'Depth-integral of Zonal Pressure Force Acceleration', 'm2 s-2', & + conversion=GV%H_to_m*US%L_T2_to_m_s2) + if(CS%id_intz_PFu_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hu,IsdB,IedB,jsd,jed,nz) + + CS%id_intz_PFv_2d = register_diag_field('ocean_model', 'intz_PFv_2d', diag%axesCv1, Time, & + 'Depth-integral of Meridional Pressure Force Acceleration', 'm2 s-2', & + conversion=GV%H_to_m*US%L_T2_to_m_s2) + if(CS%id_intz_PFv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,JsdB,JedB,nz) CS%id_hf_CAu_2d = register_diag_field('ocean_model', 'hf_CAu_2d', diag%axesCu1, Time, & 'Depth-sum Fractional Thickness-weighted Zonal Coriolis and Advective Acceleration', 'm s-2', & @@ -1396,7 +1460,17 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param CS%id_hf_CAv_2d = register_diag_field('ocean_model', 'hf_CAv_2d', diag%axesCv1, Time, & 'Depth-sum Fractional Thickness-weighted Meridional Coriolis and Advective Acceleration', 'm s-2', & conversion=US%L_T2_to_m_s2) - if(CS%id_hf_CAv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz) + if(CS%id_hf_CAv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz) + + CS%id_intz_CAu_2d = register_diag_field('ocean_model', 'intz_CAu_2d', diag%axesCu1, Time, & + 'Depth-integral of Zonal Coriolis and Advective Acceleration', 'm2 s-2', & + conversion=GV%H_to_m*US%L_T2_to_m_s2) + if(CS%id_intz_CAu_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hu,IsdB,IedB,jsd,jed,nz) + + CS%id_intz_CAv_2d = register_diag_field('ocean_model', 'intz_CAv_2d', diag%axesCv1, Time, & + 'Depth-integral of Meridional Coriolis and Advective Acceleration', 'm2 s-2', & + conversion=GV%H_to_m*US%L_T2_to_m_s2) + if(CS%id_intz_CAv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,JsdB,JedB,nz) CS%id_uav = register_diag_field('ocean_model', 'uav', diag%axesCuL, Time, & 'Barotropic-step Averaged Zonal Velocity', 'm s-1', conversion=US%L_T_to_m_s) @@ -1416,7 +1490,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param !CS%id_hf_v_BT_accel = register_diag_field('ocean_model', 'hf_v_BT_accel', diag%axesCvL, Time, & ! 'Fractional Thickness-weighted Barotropic Anomaly Meridional Acceleration', 'm s-2', & ! v_extensive=.true., conversion=US%L_T2_to_m_s2) - !if(CS%id_hf_v_BT_accel > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz) + !if(CS%id_hf_v_BT_accel > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz) CS%id_hf_u_BT_accel_2d = register_diag_field('ocean_model', 'hf_u_BT_accel_2d', diag%axesCu1, Time, & 'Depth-sum Fractional Thickness-weighted Barotropic Anomaly Zonal Acceleration', 'm s-2', & @@ -1426,7 +1500,17 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param CS%id_hf_v_BT_accel_2d = register_diag_field('ocean_model', 'hf_v_BT_accel_2d', diag%axesCv1, Time, & 'Depth-sum Fractional Thickness-weighted Barotropic Anomaly Meridional Acceleration', 'm s-2', & conversion=US%L_T2_to_m_s2) - if(CS%id_hf_v_BT_accel_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz) + if(CS%id_hf_v_BT_accel_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz) + + CS%id_intz_u_BT_accel_2d = register_diag_field('ocean_model', 'intz_u_BT_accel_2d', diag%axesCu1, Time, & + 'Depth-integral of Barotropic Anomaly Zonal Acceleration', 'm2 s-2', & + conversion=GV%H_to_m*US%L_T2_to_m_s2) + if(CS%id_intz_u_BT_accel_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hu,IsdB,IedB,jsd,jed,nz) + + CS%id_intz_v_BT_accel_2d = register_diag_field('ocean_model', 'intz_v_BT_accel_2d', diag%axesCv1, Time, & + 'Depth-integral of Barotropic Anomaly Meridional Acceleration', 'm2 s-2', & + conversion=GV%H_to_m*US%L_T2_to_m_s2) + if(CS%id_intz_v_BT_accel_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,JsdB,JedB,nz) id_clock_Cor = cpu_clock_id('(Ocean Coriolis & mom advection)', grain=CLOCK_MODULE) id_clock_continuity = cpu_clock_id('(Ocean continuity equation)', grain=CLOCK_MODULE) diff --git a/src/core/MOM_isopycnal_slopes.F90 b/src/core/MOM_isopycnal_slopes.F90 index e1f573f6ea..98b5b10998 100644 --- a/src/core/MOM_isopycnal_slopes.F90 +++ b/src/core/MOM_isopycnal_slopes.F90 @@ -37,14 +37,14 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, & !! thermodynamic variables real, intent(in) :: dt_kappa_smooth !< A smoothing vertical diffusivity !! times a smoothing timescale [Z2 ~> m2]. - real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: slope_x !< Isopycnal slope in i-direction [nondim] - real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(inout) :: slope_y !< Isopycnal slope in j-direction [nondim] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: slope_x !< Isopycnal slope in i-dir [Z L-1 ~> nondim] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(inout) :: slope_y !< Isopycnal slope in j-dir [Z L-1 ~> nondim] real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), & optional, intent(inout) :: N2_u !< Brunt-Vaisala frequency squared at - !! interfaces between u-points [T-2 ~> s-2] + !! interfaces between u-points [L2 Z-2 T-2 ~> s-2] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), & optional, intent(inout) :: N2_v !< Brunt-Vaisala frequency squared at - !! interfaces between u-points [T-2 ~> s-2] + !! interfaces between v-points [L2 Z-2 T-2 ~> s-2] integer, optional, intent(in) :: halo !< Halo width over which to compute type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure. @@ -86,7 +86,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, & real :: drdz ! Vertical density gradient [R Z-1 ~> kg m-4]. real :: Slope ! The slope of density surfaces, calculated in a way ! that is always between -1 and 1. - real :: mag_grad2 ! The squared magnitude of the 3-d density gradient [R2 L-2 ~> kg2 m-8]. + real :: mag_grad2 ! The squared magnitude of the 3-d density gradient [R2 Z-2 ~> kg2 m-8]. real :: slope2_Ratio ! The ratio of the slope squared to slope_max squared. real :: h_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected [H ~> m or kg m-2]. @@ -94,7 +94,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, & real :: dz_neglect ! A change in interface heighs that is so small it is usually lost ! in roundoff and can be neglected [Z ~> m]. logical :: use_EOS ! If true, density is calculated from T & S using an equation of state. - real :: G_Rho0 ! The gravitational acceleration divided by density [Z2 T-2 R-1 ~> m5 kg-2 s-2] + real :: G_Rho0 ! The gravitational acceleration divided by density [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1] real :: Z_to_L ! A conversion factor between from units for e to the ! units for lateral distances. real :: L_to_Z ! A conversion factor between from units for lateral distances @@ -134,7 +134,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, & present_N2_u = PRESENT(N2_u) present_N2_v = PRESENT(N2_v) - G_Rho0 = (US%L_to_Z*L_to_Z*GV%g_Earth) / GV%Rho0 + G_Rho0 = GV%g_Earth / GV%Rho0 if (present_N2_u) then do j=js,je ; do I=is-1,ie N2_u(I,j,1) = 0. @@ -248,17 +248,17 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, & ! This estimate of slope is accurate for small slopes, but bounded ! to be between -1 and 1. - mag_grad2 = drdx**2 + (L_to_Z*drdz)**2 + mag_grad2 = (Z_to_L*drdx)**2 + drdz**2 if (mag_grad2 > 0.0) then slope_x(I,j,K) = drdx / sqrt(mag_grad2) else ! Just in case mag_grad2 = 0 ever. slope_x(I,j,K) = 0.0 endif - if (present_N2_u) N2_u(I,j,k) = G_Rho0 * drdz * G%mask2dCu(I,j) ! Square of buoyancy frequency [T-2 ~> s-2] + if (present_N2_u) N2_u(I,j,k) = G_Rho0 * drdz * G%mask2dCu(I,j) ! Square of buoyancy freq. [L2 Z-2 T-2 ~> s-2] else ! With .not.use_EOS, the layers are constant density. - slope_x(I,j,K) = (Z_to_L*(e(i,j,K)-e(i+1,j,K))) * G%IdxCu(I,j) + slope_x(I,j,K) = (e(i,j,K)-e(i+1,j,K)) * G%IdxCu(I,j) endif if (local_open_u_BC) then l_seg = OBC%segnum_u(I,j) @@ -351,17 +351,17 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, & ! This estimate of slope is accurate for small slopes, but bounded ! to be between -1 and 1. - mag_grad2 = drdy**2 + (L_to_Z*drdz)**2 + mag_grad2 = (Z_to_L*drdy)**2 + drdz**2 if (mag_grad2 > 0.0) then slope_y(i,J,K) = drdy / sqrt(mag_grad2) else ! Just in case mag_grad2 = 0 ever. slope_y(i,J,K) = 0.0 endif - if (present_N2_v) N2_v(i,J,k) = G_Rho0 * drdz * G%mask2dCv(i,J) ! Square of buoyancy frequency [T-2 ~> s-2] + if (present_N2_v) N2_v(i,J,k) = G_Rho0 * drdz * G%mask2dCv(i,J) ! Square of buoyancy freq. [L2 Z-2 T-2 ~> s-2] else ! With .not.use_EOS, the layers are constant density. - slope_y(i,J,K) = (Z_to_L*(e(i,j,K)-e(i,j+1,K))) * G%IdyCv(i,J) + slope_y(i,J,K) = (e(i,j,K)-e(i,j+1,K)) * G%IdyCv(i,J) endif if (local_open_v_BC) then l_seg = OBC%segnum_v(i,J) diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index d81cf28e17..886ee77510 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -188,6 +188,9 @@ module MOM_variables real, pointer :: diag_hfrac_u(:,:,:) => NULL() !< Fractional layer thickness at u points real, pointer :: diag_hfrac_v(:,:,:) => NULL() !< Fractional layer thickness at v points + real, pointer :: diag_hu(:,:,:) => NULL() !< layer thickness at u points + real, pointer :: diag_hv(:,:,:) => NULL() !< layer thickness at v points + end type accel_diag_ptrs !> Pointers to arrays with transports, which can later be used for derived diagnostics, like energy balances. diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 47d322dfa0..bb5243892d 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -1767,7 +1767,7 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag ! call safe_alloc_ptr(CS%dv_dt,isd,ied,JsdB,JedB,nz) ! call register_time_deriv(lbound(MIS%v), MIS%v, CS%dv_dt, CS) ! endif - ! call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz) + ! call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz) !endif CS%id_hf_du_dt_2d = register_diag_field('ocean_model', 'hf_dudt_2d', diag%axesCu1, Time, & @@ -1787,7 +1787,7 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag call safe_alloc_ptr(CS%dv_dt,isd,ied,JsdB,JedB,nz) call register_time_deriv(lbound(MIS%v), MIS%v, CS%dv_dt, CS) endif - call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz) + call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz) endif ! layer thickness variables diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index d6390c9453..87ee49d449 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -722,6 +722,17 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step, CS) endif endif + ! Melting has been computed, now is time to update thickness and mass with dynamic ice shelf + if (CS%active_shelf_dynamics) then + call change_thickness_using_melt(ISS, G, US, US%s_to_T*time_step, fluxes, CS%density_ice, CS%debug) + + if (CS%debug) then + call hchksum(ISS%h_shelf, "h_shelf after change thickness using melt", G%HI, haloshift=0, scale=US%Z_to_m) + call hchksum(ISS%mass_shelf, "mass_shelf after change thickness using melt", G%HI, haloshift=0, & + scale=US%RZ_to_kg_m2) + endif + endif + if (CS%debug) call MOM_forcing_chksum("Before add shelf flux", fluxes, G, CS%US, haloshift=0) call add_shelf_flux(G, US, CS, sfc_state, fluxes) diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index a70ae45137..cf6845599b 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -24,6 +24,7 @@ module MOM_ice_shelf_dynamics use MOM_ice_shelf_state, only : ice_shelf_state use MOM_coms, only : reproducing_sum, sum_across_PEs, max_across_PEs, min_across_PEs use MOM_checksums, only : hchksum, qchksum +use MOM_ice_shelf_initialize, only : initialize_ice_shelf_boundary_channel,initialize_ice_flow_from_file implicit none ; private @@ -44,7 +45,10 @@ module MOM_ice_shelf_dynamics !! on q-points (B grid) [L T-1 ~> m s-1] real, pointer, dimension(:,:) :: v_shelf => NULL() !< the meridional velocity of the ice shelf/sheet !! on q-points (B grid) [L T-1 ~> m s-1] - + real, pointer, dimension(:,:) :: taudx_shelf => NULL() !< the driving stress of the ice shelf/sheet + !! on q-points (C grid) [Pa ~> Pa] + real, pointer, dimension(:,:) :: taudy_shelf => NULL() !< the meridional stress of the ice shelf/sheet + !! on q-points (C grid) [Pa ~> Pa] real, pointer, dimension(:,:) :: u_face_mask => NULL() !< mask for velocity boundary conditions on the C-grid !! u-face - this is because the FEM cares about FACES THAT GET INTEGRATED OVER, !! not vertices. Will represent boundary conditions on computational boundary @@ -152,12 +156,14 @@ module MOM_ice_shelf_dynamics !>@{ Diagnostic handles integer :: id_u_shelf = -1, id_v_shelf = -1, id_t_shelf = -1, & + id_taudx_shelf = -1, id_taudy_shelf = -1, & id_ground_frac = -1, id_col_thick = -1, id_OD_av = -1, & id_u_mask = -1, id_v_mask = -1, id_t_mask = -1 !>@} ! ids for outputting intermediate thickness in advection subroutine (debugging) - !integer :: id_h_after_uflux = -1, id_h_after_vflux = -1, id_h_after_adv = -1 - + !>@{ Diagnostic handles for debugging + integer :: id_h_after_uflux = -1, id_h_after_vflux = -1, id_h_after_adv = -1, id_visc_shelf = -1 + !>@} type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to control diagnostic output. end type ice_shelf_dyn_CS @@ -250,7 +256,8 @@ subroutine register_ice_shelf_dyn_restarts(G, param_file, CS, restart_CS) allocate( CS%basal_traction(isd:ied,jsd:jed) ) ; CS%basal_traction(:,:) = 0.0 allocate( CS%OD_av(isd:ied,jsd:jed) ) ; CS%OD_av(:,:) = 0.0 allocate( CS%ground_frac(isd:ied,jsd:jed) ) ; CS%ground_frac(:,:) = 0.0 - + allocate( CS%taudx_shelf(Isd:Ied,Jsd:Jed) ) ; CS%taudx_shelf(:,:) = 0.0 + allocate( CS%taudy_shelf(Isd:Ied,Jsd:Jed) ) ; CS%taudy_shelf(:,:) = 0.0 ! additional restarts for ice shelf state call register_restart_field(CS%u_shelf, "u_shelf", .false., restart_CS, & "ice sheet/shelf u-velocity", "m s-1", hor_grid='Bu') @@ -258,6 +265,10 @@ subroutine register_ice_shelf_dyn_restarts(G, param_file, CS, restart_CS) "ice sheet/shelf v-velocity", "m s-1", hor_grid='Bu') call register_restart_field(CS%t_shelf, "t_shelf", .true., restart_CS, & "ice sheet/shelf vertically averaged temperature", "deg C") + call register_restart_field(CS%taudx_shelf, "taudx_shelf", .true., restart_CS, & + "ice sheet/shelf taudx-driving stress", "kPa") + call register_restart_field(CS%taudy_shelf, "taudy_shelf", .true., restart_CS, & + "ice sheet/shelf taudy-driving stress", "kPa") call register_restart_field(CS%OD_av, "OD_av", .true., restart_CS, & "Average open ocean depth in a cell","m") call register_restart_field(CS%ground_frac, "ground_frac", .true., restart_CS, & @@ -367,20 +378,20 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ call get_param(param_file, mdl, "A_GLEN_ISOTHERM", CS%A_glen_isothermal, & "Ice viscosity parameter in Glen's Law", & - units="Pa-3 yr-1", default=9.461e-18, scale=1.0/(365.0*86400.0)) + units="Pa-3 s-1", default=2.2261e-25, scale=1.0) ! This default is equivalent to 3.0001e-25 Pa-3 s-1, appropriate at about -10 C. call get_param(param_file, mdl, "GLEN_EXPONENT", CS%n_glen, & "nonlinearity exponent in Glen's Law", & units="none", default=3.) call get_param(param_file, mdl, "MIN_STRAIN_RATE_GLEN", CS%eps_glen_min, & "min. strain rate to avoid infinite Glen's law viscosity", & - units="a-1", default=1.e-12, scale=US%T_to_s/(365.0*86400.0)) + units="s-1", default=1.e-19, scale=US%T_to_s) call get_param(param_file, mdl, "BASAL_FRICTION_EXP", CS%n_basal_fric, & "Exponent in sliding law \tau_b = C u^(n_basal_fric)", & units="none", fail_if_missing=.true.) call get_param(param_file, mdl, "BASAL_FRICTION_COEFF", CS%C_basal_friction, & "Coefficient in sliding law \tau_b = C u^(n_basal_fric)", & - units="Pa (m yr-1)-(n_basal_fric)", scale=US%kg_m2s_to_RZ_T*((365.0*86400.0)**CS%n_basal_fric), & + units="Pa (m s-1)^(n_basal_fric)", scale=US%kg_m2s_to_RZ_T**CS%n_basal_fric, & fail_if_missing=.true.) call get_param(param_file, mdl, "DENSITY_ICE", CS%density_ice, & "A typical density of ice.", units="kg m-3", default=917.0, scale=US%kg_m3_to_R) @@ -400,10 +411,11 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ call get_param(param_file, mdl, "SHELF_MOVING_FRONT", CS%moving_shelf_front, & "Specify whether to advance shelf front (and calve).", & - default=.true.) + default=.false.) call get_param(param_file, mdl, "CALVE_TO_MASK", CS%calve_to_mask, & "If true, do not allow an ice shelf where prohibited by a mask.", & default=.false.) + endif call get_param(param_file, mdl, "MIN_THICKNESS_SIMPLE_CALVE", CS%min_thickness_simple_calve, & "Min thickness rule for the VERY simple calving law",& @@ -411,10 +423,8 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ ! Allocate memory in the ice shelf dynamics control structure that was not ! previously allocated for registration for restarts. - ! OVS vertically integrated Temperature if (active_shelf_dynamics) then - ! DNG allocate( CS%u_bdry_val(Isdq:Iedq,Jsdq:Jedq) ) ; CS%u_bdry_val(:,:) = 0.0 allocate( CS%v_bdry_val(Isdq:Iedq,Jsdq:Jedq) ) ; CS%v_bdry_val(:,:) = 0.0 allocate( CS%t_bdry_val(isd:ied,jsd:jed) ) ; CS%t_bdry_val(:,:) = -15.0 @@ -516,46 +526,63 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ call pass_var(CS%calve_mask,G%domain) endif -! call init_boundary_values(CS, G, time, ISS%hmask, CS%input_flux, CS%input_thickness, new_sim) - - if (new_sim) then - call MOM_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: initialize ice velocity.") - call update_OD_ffrac_uncoupled(CS, G, ISS%h_shelf(:,:)) - call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf, iters, Time) - - if (CS%id_u_shelf > 0) call post_data(CS%id_u_shelf, CS%u_shelf, CS%diag) - if (CS%id_v_shelf > 0) call post_data(CS%id_v_shelf, CS%v_shelf,CS%diag) - endif + call initialize_ice_shelf_boundary_channel(CS%u_face_mask_bdry, CS%v_face_mask_bdry, & + CS%u_flux_bdry_val, CS%v_flux_bdry_val, CS%u_bdry_val, CS%v_bdry_val, CS%u_shelf, CS%v_shelf,& + CS%h_bdry_val, & + CS%thickness_bdry_val, ISS%hmask, ISS%h_shelf, G, & + US, param_file ) + call pass_var(ISS%hmask, G%domain) + call pass_var(CS%h_bdry_val, G%domain) + call pass_var(CS%thickness_bdry_val, G%domain) + call pass_var(CS%u_bdry_val, G%domain) + call pass_var(CS%v_bdry_val, G%domain) + call pass_var(CS%u_face_mask_bdry, G%domain) + call pass_var(CS%v_face_mask_bdry, G%domain) + !call init_boundary_values(CS, G, time, ISS%hmask, CS%input_flux, CS%input_thickness, new_sim) + call update_velocity_masks(CS, G, ISS%hmask, CS%umask, CS%vmask, CS%u_face_mask, CS%v_face_mask) ! Register diagnostics. - CS%id_u_shelf = register_diag_field('ocean_model','u_shelf',CS%diag%axesCu1, Time, & + CS%id_u_shelf = register_diag_field('ice_shelf_model','u_shelf',CS%diag%axesCu1, Time, & 'x-velocity of ice', 'm yr-1', conversion=365.0*86400.0*US%L_T_to_m_s) - CS%id_v_shelf = register_diag_field('ocean_model','v_shelf',CS%diag%axesCv1, Time, & + CS%id_v_shelf = register_diag_field('ice_shelf_model','v_shelf',CS%diag%axesCv1, Time, & 'y-velocity of ice', 'm yr-1', conversion=365.0*86400.0*US%L_T_to_m_s) - CS%id_u_mask = register_diag_field('ocean_model','u_mask',CS%diag%axesCu1, Time, & + CS%id_taudx_shelf = register_diag_field('ice_shelf_model','taudx_shelf',CS%diag%axesT1, Time, & + 'x-driving stress of ice', 'kPa', conversion=1.e-9*US%L_T_to_m_s) + CS%id_taudy_shelf = register_diag_field('ice_shelf_model','taudy_shelf',CS%diag%axesT1, Time, & + 'x-driving stress of ice', 'kPa', conversion=1.e-9*US%L_T_to_m_s) + CS%id_u_mask = register_diag_field('ice_shelf_model','u_mask',CS%diag%axesCu1, Time, & 'mask for u-nodes', 'none') - CS%id_v_mask = register_diag_field('ocean_model','v_mask',CS%diag%axesCv1, Time, & + CS%id_v_mask = register_diag_field('ice_shelf_model','v_mask',CS%diag%axesCv1, Time, & 'mask for v-nodes', 'none') -! CS%id_surf_elev = register_diag_field('ocean_model','ice_surf',CS%diag%axesT1, Time, & -! 'ice surf elev', 'm') - CS%id_ground_frac = register_diag_field('ocean_model','ice_ground_frac',CS%diag%axesT1, Time, & + CS%id_ground_frac = register_diag_field('ice_shelf_model','ice_ground_frac',CS%diag%axesT1, Time, & 'fraction of cell that is grounded', 'none') - CS%id_col_thick = register_diag_field('ocean_model','col_thick',CS%diag%axesT1, Time, & + + CS%id_col_thick = register_diag_field('ice_shelf_model','col_thick',CS%diag%axesT1, Time, & 'ocean column thickness passed to ice model', 'm', conversion=US%Z_to_m) - CS%id_OD_av = register_diag_field('ocean_model','OD_av',CS%diag%axesT1, Time, & + CS%id_visc_shelf = register_diag_field('ice_shelf_model','ice_visc',CS%diag%axesT1, Time, & + 'viscosity', 'm', conversion=1e-6*US%Z_to_m) + CS%id_OD_av = register_diag_field('ice_shelf_model','OD_av',CS%diag%axesT1, Time, & 'intermediate ocean column thickness passed to ice model', 'm', conversion=US%Z_to_m) - !CS%id_h_after_uflux = register_diag_field('ocean_model','h_after_uflux',CS%diag%axesh1, Time, & - ! 'thickness after u flux ', 'none') - !CS%id_h_after_vflux = register_diag_field('ocean_model','h_after_vflux',CS%diag%axesh1, Time, & - ! 'thickness after v flux ', 'none') - !CS%id_h_after_adv = register_diag_field('ocean_model','h_after_adv',CS%diag%axesh1, Time, & - ! 'thickness after front adv ', 'none') - -!!! OVS vertically integrated temperature - CS%id_t_shelf = register_diag_field('ocean_model','t_shelf',CS%diag%axesT1, Time, & - 'T of ice', 'oC') - CS%id_t_mask = register_diag_field('ocean_model','tmask',CS%diag%axesT1, Time, & - 'mask for T-nodes', 'none') + CS%id_h_after_uflux = register_diag_field('ice_shelf_model','h_after_uflux',CS%diag%axesT1, Time, & + 'thickness after u flux ', 'none') + CS%id_h_after_vflux = register_diag_field('ice_shelf_model','h_after_vflux',CS%diag%axesT1, Time, & + 'thickness after v flux ', 'none') + CS%id_h_after_adv = register_diag_field('ice_shelf_model','h_after_adv',CS%diag%axesT1, Time, & + 'thickness after front adv ', 'none') + if (new_sim) then + call MOM_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: initialize ice velocity.") + call update_OD_ffrac_uncoupled(CS, G, ISS%h_shelf(:,:)) + call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf,CS%taudx_shelf,CS%taudy_shelf, iters, Time) + if (CS%id_u_shelf > 0) call post_data(CS%id_u_shelf, CS%u_shelf, CS%diag) + if (CS%id_v_shelf > 0) call post_data(CS%id_v_shelf, CS%v_shelf,CS%diag) + if (CS%id_taudx_shelf > 0) call post_data(CS%id_taudx_shelf, CS%taudx_shelf,CS%diag) + if (CS%id_taudy_shelf > 0) call post_data(CS%id_taudy_shelf, CS%taudy_shelf,CS%diag) + if (CS%id_visc_shelf > 0) call post_data(CS%id_visc_shelf, CS%ice_visc,CS%diag) +! CS%id_t_shelf = register_diag_field('ocean_model','t_shelf',CS%diag%axesT1, Time, & +! 'T of ice', 'oC') +! CS%id_t_mask = register_diag_field('ocean_model','tmask',CS%diag%axesT1, Time, & +! 'mask for T-nodes', 'none') + endif endif end subroutine initialize_ice_shelf_dyn @@ -592,8 +619,7 @@ subroutine initialize_diagnostic_fields(CS, ISS, G, US, Time) enddo enddo - call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf, iters, dummy_time) - + call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf,CS%taudx_shelf,CS%taudy_shelf, iters, Time) end subroutine initialize_diagnostic_fields !> This function returns the global maximum advective timestep that can be taken based on the current @@ -652,7 +678,7 @@ subroutine update_ice_shelf(CS, ISS, G, US, time_step, Time, ocean_mass, coupled coupled_GL = .false. if (present(ocean_mass) .and. present(coupled_grounding)) coupled_GL = coupled_grounding - +! call ice_shelf_advect(CS, ISS, G, time_step, Time) CS%elapsed_velocity_time = CS%elapsed_velocity_time + time_step if (CS%elapsed_velocity_time >= CS%velocity_update_time_step) update_ice_vel = .true. @@ -663,8 +689,9 @@ subroutine update_ice_shelf(CS, ISS, G, US, time_step, Time, ocean_mass, coupled call update_OD_ffrac_uncoupled(CS, G, ISS%h_shelf(:,:)) endif + if (update_ice_vel) then - call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf, iters, Time) + call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf,CS%taudx_shelf,CS%taudy_shelf, iters, Time) endif call ice_shelf_temp(CS, ISS, G, US, time_step, ISS%water_flux, Time) @@ -675,8 +702,11 @@ subroutine update_ice_shelf(CS, ISS, G, US, time_step, Time, ocean_mass, coupled if (CS%id_u_shelf > 0) call post_data(CS%id_u_shelf, CS%u_shelf, CS%diag) if (CS%id_v_shelf > 0) call post_data(CS%id_v_shelf, CS%v_shelf, CS%diag) if (CS%id_t_shelf > 0) call post_data(CS%id_t_shelf,CS%t_shelf,CS%diag) + if (CS%id_taudx_shelf > 0) call post_data(CS%id_taudx_shelf, CS%taudx_shelf, CS%diag) + if (CS%id_taudy_shelf > 0) call post_data(CS%id_taudy_shelf, CS%taudy_shelf, CS%diag) if (CS%id_ground_frac > 0) call post_data(CS%id_ground_frac, CS%ground_frac,CS%diag) if (CS%id_OD_av >0) call post_data(CS%id_OD_av, CS%OD_av,CS%diag) + if (CS%id_visc_shelf > 0) call post_data(CS%id_visc_shelf, CS%ice_visc,CS%diag) if (CS%id_u_mask > 0) call post_data(CS%id_u_mask,CS%umask,CS%diag) if (CS%id_v_mask > 0) call post_data(CS%id_v_mask,CS%vmask,CS%diag) @@ -738,16 +768,16 @@ subroutine ice_shelf_advect(CS, ISS, G, time_step, Time) call ice_shelf_advect_thickness_x(CS, G, LB, time_step, ISS%hmask, ISS%h_shelf, h_after_uflux, uh_ice) ! call enable_averages(time_step, Time, CS%diag) -! call pass_var(h_after_uflux, G%domain) -! if (CS%id_h_after_uflux > 0) call post_data(CS%id_h_after_uflux, h_after_uflux, CS%diag) + call pass_var(h_after_uflux, G%domain) + if (CS%id_h_after_uflux > 0) call post_data(CS%id_h_after_uflux, h_after_uflux, CS%diag) ! call disable_averaging(CS%diag) LB%ish = G%isc ; LB%ieh = G%iec ; LB%jsh = G%jsc ; LB%jeh = G%jec call ice_shelf_advect_thickness_y(CS, G, LB, time_step, ISS%hmask, h_after_uflux, h_after_vflux, vh_ice) ! call enable_averages(time_step, Time, CS%diag) -! call pass_var(h_after_vflux, G%domain) -! if (CS%id_h_after_vflux > 0) call post_data(CS%id_h_after_vflux, h_after_vflux, CS%diag) + call pass_var(h_after_vflux, G%domain) + if (CS%id_h_after_vflux > 0) call post_data(CS%id_h_after_vflux, h_after_vflux, CS%diag) ! call disable_averaging(CS%diag) do j=jsd,jed @@ -777,7 +807,9 @@ subroutine ice_shelf_advect(CS, ISS, G, time_step, Time) end subroutine ice_shelf_advect -subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, iters, time) +!>This subroutine computes u- and v-velocities of the ice shelf iterating on non-linear ice viscosity +!subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, iters, time) + subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf,taudx,taudy, iters, Time) type(ice_shelf_dyn_CS), intent(inout) :: CS !< The ice shelf dynamics control structure type(ice_shelf_state), intent(in) :: ISS !< A structure with elements that describe !! the ice-shelf state @@ -790,7 +822,10 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, iters, time) integer, intent(out) :: iters !< The number of iterations used in the solver. type(time_type), intent(in) :: Time !< The current model time - real, dimension(SZDIB_(G),SZDJB_(G)) :: taudx, taudy ! Driving stresses at q-points [R L3 Z T-2 ~> kg m s-2] + real, dimension(SZDIB_(G),SZDJB_(G)), & + intent(out) :: taudx !< Driving x-stress at q-points [R L3 Z T-2 ~> kg m s-2] + real, dimension(SZDIB_(G),SZDJB_(G)), & + intent(out) :: taudy !< Driving y-stress at q-points [R L3 Z T-2 ~> kg m s-2] real, dimension(SZDIB_(G),SZDJB_(G)) :: u_bdry_cont ! Boundary u-stress contribution [R L3 Z T-2 ~> kg m s-2] real, dimension(SZDIB_(G),SZDJB_(G)) :: v_bdry_cont ! Boundary v-stress contribution [R L3 Z T-2 ~> kg m s-2] real, dimension(SZDIB_(G),SZDJB_(G)) :: Au, Av ! The retarding lateral stress contributions [R L3 Z T-2 ~> kg m s-2] @@ -824,10 +859,20 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, iters, time) ! need to make these conditional on GL interpolation float_cond(:,:) = 0.0 ; H_node(:,:) = 0.0 + CS%ground_frac(:,:) = 0.0 allocate(Phisub(nsub,nsub,2,2,2,2)) ; Phisub(:,:,:,:,:,:) = 0.0 - call calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, CS%OD_av) + do j=G%jsc,G%jec + do i=G%isc,G%iec + if (rhoi_rhow * ISS%h_shelf(i,j) - G%bathyT(i,j) > 0) then + float_cond(i,j) = 1.0 + CS%ground_frac(i,j) = 1.0 + endif + enddo + enddo + call calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, CS%OD_av) + call pass_vector(taudx, taudy, G%domain, TO_ALL, BGRID_NE) ! This is to determine which cells contain the grounding line, the criterion being that the cell ! is ice-covered, with some nodes floating and some grounded flotation condition is estimated by ! assuming topography is cellwise constant and H is bilinear in a cell; floating where @@ -868,8 +913,8 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, iters, time) enddo ; enddo call calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) - call pass_var(CS%ice_visc, G%domain) + call calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf) call pass_var(CS%basal_traction, G%domain) ! This makes sure basal stress is only applied when it is supposed to be @@ -885,7 +930,7 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, iters, time) call CG_action(Au, Av, u_shlf, v_shlf, Phi, Phisub, CS%umask, CS%vmask, ISS%hmask, H_node, & CS%ice_visc, float_cond, G%bathyT, CS%basal_traction, & G, US, G%isc-1, G%iec+1, G%jsc-1, G%jec+1, rhoi_rhow) - + call pass_vector(Au,Av,G%domain) if (CS%nonlin_solve_err_mode == 1) then err_init = 0 ; err_tempu = 0 ; err_tempv = 0 do J=G%IscB,G%JecB ; do I=G%IscB,G%IecB @@ -921,6 +966,7 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, iters, time) call calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) call pass_var(CS%ice_visc, G%domain) + call calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf) call pass_var(CS%basal_traction, G%domain) ! makes sure basal stress is only applied when it is supposed to be @@ -987,8 +1033,11 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, iters, time) call MOM_mesg(mesg, 5) if (err_max <= CS%nonlinear_tolerance * err_init) then + write(mesg,*) "ice_shelf_solve_outer: nonlinear fractional residual = ", err_max/err_init + call MOM_mesg(mesg) write(mesg,*) "ice_shelf_solve_outer: exiting nonlinear solve after ",iter," iterations" - call MOM_mesg(mesg, 5) +! call MOM_mesg(mesg, 5) + call MOM_mesg(mesg) exit endif @@ -1074,7 +1123,7 @@ subroutine ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H rhoi_rhow = CS%density_ice / CS%density_ocean_avg Zu(:,:) = 0 ; Zv(:,:) = 0 ; DIAGu(:,:) = 0 ; DIAGv(:,:) = 0 - Ru(:,:) = 0 ; Rv(:,:) = 0 ; Au(:,:) = 0 ; Av(:,:) = 0 + Ru(:,:) = 0 ; Rv(:,:) = 0 ; Au(:,:) = 0 ; Av(:,:) = 0 ; RHSu(:,:) = 0 ; RHSv(:,:) = 0 Du(:,:) = 0 ; Dv(:,:) = 0 ; ubd(:,:) = 0 ; vbd(:,:) = 0 dot_p1 = 0 @@ -1126,8 +1175,8 @@ subroutine ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H do j=jsdq,jedq do i=isdq,iedq - if (CS%umask(I,J) == 1) Zu(I,J) = Ru(I,J) / DIAGu(I,J) - if (CS%vmask(I,J) == 1) Zv(I,J) = Rv(I,J) / DIAGv(I,J) + if (CS%umask(I,J) == 1 .AND.(DIAGu(I,J)/=0)) Zu(I,J) = Ru(I,J) / DIAGu(I,J) + if (CS%vmask(I,J) == 1 .AND.(DIAGv(I,J)/=0)) Zv(I,J) = Rv(I,J) / DIAGv(I,J) enddo enddo @@ -1162,7 +1211,7 @@ subroutine ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H ! Au, Av valid region moves in by 1 - + call pass_vector(Au,Av,G%domain, TO_ALL, BGRID_NE) sum_vec(:,:) = 0.0 ; sum_vec_2(:,:) = 0.0 do j=jscq,jecq ; do i=iscq,iecq @@ -1206,10 +1255,10 @@ subroutine ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H do j=jsdq,jedq do i=isdq,iedq - if (CS%umask(I,J) == 1) then + if (CS%umask(I,J) == 1 .AND.(DIAGu(I,J)/=0)) then Zu(I,J) = Ru(I,J) / DIAGu(I,J) endif - if (CS%vmask(I,J) == 1) then + if (CS%vmask(I,J) == 1 .AND.(DIAGv(I,J)/=0)) then Zv(I,J) = Rv(I,J) / DIAGv(I,J) endif enddo @@ -1264,9 +1313,11 @@ subroutine ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H cg_halo = cg_halo - 1 if (cg_halo == 0) then - ! pass vectors + ! pass vectors call pass_vector(Du, Dv, G%domain, TO_ALL, BGRID_NE) call pass_vector(u_shlf, v_shlf, G%domain, TO_ALL, BGRID_NE) + call pass_var(u_shlf, G%domain) + call pass_var(v_shlf, G%domain) call pass_vector(Ru, Rv, G%domain, TO_ALL, BGRID_NE) cg_halo = 3 endif @@ -1733,7 +1784,7 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) BASE ! basal elevation of shelf/stream [Z ~> m]. - real :: rho, rhow ! Ice and ocean densities [R ~> kg m-3] + real :: rho, rhow, rhoi_rhow ! Ice and ocean densities [R ~> kg m-3] real :: sx, sy ! Ice shelf top slopes [Z L-1 ~> m s-1] real :: neumann_val ! [R Z L2 T-2 ~> kg s-2] real :: dxh, dyh ! Local grid spacing [L ~> m] @@ -1747,21 +1798,35 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) iscq = G%iscB ; iecq = G%iecB ; jscq = G%jscB ; jecq = G%jecB isd = G%isd ; jsd = G%jsd iegq = G%iegB ; jegq = G%jegB - gisc = G%domain%nihalo+1 ; gjsc = G%domain%njhalo+1 - giec = G%domain%niglobal+G%domain%nihalo ; gjec = G%domain%njglobal+G%domain%njhalo +! gisc = G%domain%nihalo+1 ; gjsc = G%domain%njhalo+1 + gisc = 1 ; gjsc = 1 +! giec = G%domain%niglobal+G%domain%nihalo ; gjec = G%domain%njglobal+G%domain%njhalo + giec = G%domain%niglobal ; gjec = G%domain%njglobal is = iscq - 1; js = jscq - 1 i_off = G%idg_offset ; j_off = G%jdg_offset rho = CS%density_ice rhow = CS%density_ocean_avg grav = CS%g_Earth - + rhoi_rhow = rho/rhow ! prelim - go through and calculate S ! or is this faster? BASE(:,:) = -G%bathyT(:,:) + OD(:,:) S(:,:) = BASE(:,:) + ISS%h_shelf(:,:) + ! check whether the ice is floating or grounded + do j=jsc-G%domain%njhalo,jec+G%domain%njhalo + do i=isc-G%domain%nihalo,iec+G%domain%nihalo + +! if (ISS%h_shelf(i,j) < rhow/rho * G%bathyT(i,j)) then + if (rhoi_rhow * ISS%h_shelf(i,j) - G%bathyT(i,j) <= 0) then + S(i,j)=(1 - rhoi_rhow)*ISS%h_shelf(i,j) + endif + + + enddo + enddo do j=jsc-1,jec+1 do i=isc-1,iec+1 cnt = 0 @@ -1774,7 +1839,7 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) ! calculate sx if ((i+i_off) == gisc) then ! at left computational bdry - if (ISS%hmask(i+1,j) == 1) then + if (ISS%hmask(i+1,j) == 1) then sx = (S(i+1,j)-S(i,j))/dxh else sx = 0 @@ -1841,23 +1906,28 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) endif ! SW vertex - taudx(I-1,J-1) = taudx(I-1,J-1) - .25 * rho * grav * ISS%h_shelf(i,j) * sx * G%areaT(i,j) - taudy(I-1,J-1) = taudy(I-1,J-1) - .25 * rho * grav * ISS%h_shelf(i,j) * sy * G%areaT(i,j) - + if (ISS%hmask(I-1,J-1) == 1) then + taudx(I-1,J-1) = taudx(I-1,J-1) - .25 * rho * grav * ISS%h_shelf(i,j) * sx * G%areaT(i,j) + taudy(I-1,J-1) = taudy(I-1,J-1) - .25 * rho * grav * ISS%h_shelf(i,j) * sy * G%areaT(i,j) + endif ! SE vertex - taudx(I,J-1) = taudx(I,J-1) - .25 * rho * grav * ISS%h_shelf(i,j) * sx * G%areaT(i,j) - taudy(I,J-1) = taudy(I,J-1) - .25 * rho * grav * ISS%h_shelf(i,j) * sy * G%areaT(i,j) - + if (ISS%hmask(I,J-1) == 1) then + taudx(I,J-1) = taudx(I,J-1) - .25 * rho * grav * ISS%h_shelf(i,j) * sx * G%areaT(i,j) + taudy(I,J-1) = taudy(I,J-1) - .25 * rho * grav * ISS%h_shelf(i,j) * sy * G%areaT(i,j) + endif ! NW vertex - taudx(I-1,J) = taudx(I-1,J) - .25 * rho * grav * ISS%h_shelf(i,j) * sx * G%areaT(i,j) - taudy(I-1,J) = taudy(I-1,J) - .25 * rho * grav * ISS%h_shelf(i,j) * sy * G%areaT(i,j) - + if (ISS%hmask(I-1,J) == 1) then + taudx(I-1,J) = taudx(I-1,J) - .25 * rho * grav * ISS%h_shelf(i,j) * sx * G%areaT(i,j) + taudy(I-1,J) = taudy(I-1,J) - .25 * rho * grav * ISS%h_shelf(i,j) * sy * G%areaT(i,j) + endif ! NE vertex - taudx(I,J) = taudx(I,J) - .25 * rho * grav * ISS%h_shelf(i,j) * sx * G%areaT(i,j) - taudy(I,J) = taudy(I,J) - .25 * rho * grav * ISS%h_shelf(i,j) * sy * G%areaT(i,j) - + if (ISS%hmask(I,J) == 1) then + taudx(I,J) = taudx(I,J) - .25 * rho * grav * ISS%h_shelf(i,j) * sx * G%areaT(i,j) + taudy(I,J) = taudy(I,J) - .25 * rho * grav * ISS%h_shelf(i,j) * sy * G%areaT(i,j) + endif if (CS%ground_frac(i,j) == 1) then - neumann_val = .5 * grav * (rho * ISS%h_shelf(i,j)**2 - rhow * G%bathyT(i,j)**2) +! neumann_val = .5 * grav * (rho * ISS%h_shelf(i,j)**2 - rhow * G%bathyT(i,j)**2) + neumann_val = .5 * grav * (1-rho/rhow) * rho * ISS%h_shelf(i,j)**2 else neumann_val = .5 * grav * (1-rho/rhow) * rho * ISS%h_shelf(i,j)**2 endif @@ -1977,7 +2047,7 @@ subroutine CG_action(uret, vret, u_shlf, v_shlf, Phi, Phisub, umask, vmask, hmas intent(inout) :: uret !< The retarding stresses working at u-points [R L3 Z T-2 ~> kg m s-2]. real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), & intent(inout) :: vret !< The retarding stresses working at v-points [R L3 Z T-2 ~> kg m s-2]. - real, dimension(SZDI_(G),SZDJ_(G),8,4), & + real, dimension(8,4,SZDI_(G),SZDJ_(G)), & intent(in) :: Phi !< The gradients of bilinear basis elements at Gaussian !! quadrature points surrounding the cell vertices [L-1 ~> m-1]. real, dimension(:,:,:,:,:,:), & @@ -2081,7 +2151,7 @@ subroutine CG_action(uret, vret, u_shlf, v_shlf, Phi, Phisub, umask, vmask, hmas v_shlf(I-1,J) * Phi(6,2*(jq-1)+iq,i,j) + & v_shlf(I,J) * Phi(8,2*(jq-1)+iq,i,j) - do iphi=1,2 ; do jphi=1,2 ; Itgt = I-2+iphi ; Jtgt = J-2-jphi + do iphi=1,2 ; do jphi=1,2 ; Itgt = I-2+iphi ; ;Jtgt = J-2+jphi if (umask(Itgt,Jtgt) == 1) uret(Itgt,Jtgt) = uret(Itgt,Jtgt) + 0.25 * ice_visc(i,j) * & ((4*ux+2*vy) * Phi(2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq,i,j) + & (uy+vx) * Phi(2*(2*(jphi-1)+iphi),2*(jq-1)+iq,i,j)) @@ -2215,7 +2285,7 @@ subroutine matrix_diagonal(CS, G, US, float_cond, H_node, ice_visc, basal_trac, ! Phi(2*i-1,j) gives d(Phi_i)/dx at quadrature point j ! Phi(2*i,j) gives d(Phi_i)/dy at quadrature point j - do iq=1,2 ; do jq=1,2 ; do iphi=1,2 ; do jphi=1,2 ; Itgt = I-2+iphi ; Jtgt = J-2-jphi + do iq=1,2 ; do jq=1,2 ; do iphi=1,2 ; do jphi=1,2 ; Itgt = I-2+iphi ; Jtgt = J-2+jphi ilq = 1 ; if (iq == iphi) ilq = 2 jlq = 1 ; if (jq == jphi) jlq = 2 @@ -2259,7 +2329,7 @@ subroutine matrix_diagonal(CS, G, US, float_cond, H_node, ice_visc, basal_trac, if (float_cond(i,j) == 1) then Hcell(:,:) = H_node(i-1:i,j-1:j) call CG_diagonal_subgrid_basal(Phisub, Hcell, G%bathyT(i,j), dens_ratio, sub_ground) - do iphi=1,2 ; do jphi=1,2 ; Itgt = I-2+iphi ; Jtgt = J-2-jphi + do iphi=1,2 ; do jphi=1,2 ; Itgt = I-2+iphi ; Jtgt = J-2+jphi if (CS%umask(Itgt,Jtgt) == 1) then u_diagonal(Itgt,Jtgt) = u_diagonal(Itgt,Jtgt) + sub_ground(iphi,jphi) * basal_trac(i,j) v_diagonal(Itgt,Jtgt) = v_diagonal(Itgt,Jtgt) + sub_ground(iphi,jphi) * basal_trac(i,j) @@ -2400,7 +2470,7 @@ subroutine apply_boundary_values(CS, ISS, G, US, time, Phisub, H_node, ice_visc, CS%v_bdry_val(I-1,J) * Phi(6,2*(jq-1)+iq) + & CS%v_bdry_val(I,J) * Phi(8,2*(jq-1)+iq) - do iphi=1,2 ; do jphi=1,2 ; Itgt = I-2+iphi ; Jtgt = J-2-jphi + do iphi=1,2 ; do jphi=1,2 ; Itgt = I-2+iphi ; Jtgt = J-2+jphi ilq = 1 ; if (iq == iphi) ilq = 2 jlq = 1 ; if (jq == jphi) jlq = 2 @@ -2447,6 +2517,8 @@ subroutine apply_boundary_values(CS, ISS, G, US, time, Phisub, H_node, ice_visc, endif endif ; enddo ; enddo + call pass_vector(u_bdry_contr, v_bdry_contr, G%domain, TO_ALL, BGRID_NE) + end subroutine apply_boundary_values !> Update depth integrated viscosity, based on horizontal strain rates, and also update the @@ -2468,12 +2540,15 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) ! also this subroutine updates the nonlinear part of the basal traction ! this may be subject to change later... to make it "hybrid" - - integer :: i, j, iscq, iecq, jscq, jecq, isd, jsd, ied, jed, iegq, jegq - integer :: giec, gjec, gisc, gjsc, cnt, isc, jsc, iec, jec, is, js +! real, dimension(SZDIB_(G),SZDJB_(G)) :: eII, ux, uy, vx, vy + integer :: i, j, iscq, iecq, jscq, jecq, isd, jsd, ied, jed, iegq, jegq, iq, jq + integer :: giec, gjec, gisc, gjsc, cnt, isc, jsc, iec, jec, is, js, i_off, j_off real :: Visc_coef, n_g - real :: ux, uy, vx, vy, eps_min ! Velocity shears [T-1 ~> s-1] - real :: umid, vmid, unorm ! Velocities [L T-1 ~> m s-1] + real :: ux, uy, vx, vy + real :: eps_min, dxh, dyh ! Velocity shears [T-1 ~> s-1] + real, dimension(8,4) :: Phi + real, dimension(2) :: xquad +! real :: umid, vmid, unorm ! Velocities [L T-1 ~> m s-1] isc = G%isc ; jsc = G%jsc ; iec = G%iec ; jec = G%jec iscq = G%iscB ; iecq = G%iecB ; jscq = G%jscB ; jecq = G%jecB @@ -2482,22 +2557,65 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) gisc = G%domain%nihalo+1 ; gjsc = G%domain%njhalo+1 giec = G%domain%niglobal+gisc ; gjec = G%domain%njglobal+gjsc is = iscq - 1; js = jscq - 1 + i_off = G%idg_offset ; j_off = G%jdg_offset n_g = CS%n_glen; eps_min = CS%eps_glen_min - Visc_coef = US%kg_m2s_to_RZ_T*US%m_to_L*US%Z_to_L*(CS%A_glen_isothermal)**(1./CS%n_glen) - - do j=jsd+1,jed-1 - do i=isd+1,ied-1 - - if (ISS%hmask(i,j) == 1) then - ux = ((u_shlf(I,J) + u_shlf(I,J-1)) - (u_shlf(I-1,J) + u_shlf(I-1,J-1))) / (2*G%dxT(i,j)) - vx = ((v_shlf(I,J) + v_shlf(I,J-1)) - (v_shlf(I-1,J) + v_shlf(I-1,J-1))) / (2*G%dxT(i,j)) - uy = ((u_shlf(I,J) + u_shlf(I-1,J)) - (u_shlf(I,J-1) + u_shlf(I-1,J-1))) / (2*G%dyT(i,j)) - vy = ((v_shlf(I,J) + v_shlf(I-1,J)) - (v_shlf(I,J-1) + v_shlf(I-1,J-1))) / (2*G%dyT(i,j)) + Visc_coef = US%kg_m2s_to_RZ_T*US%m_to_L*US%Z_to_L*(CS%A_glen_isothermal)**(-1./CS%n_glen) + do j=jsc,jec + do i=isc,iec + + if ((ISS%hmask(i,j) == 1) .OR. (ISS%hmask(i,j) == 3)) then + ux = ((u_shlf(I,J) + (u_shlf(I,J-1) + u_shlf(I,J+1))) - & + (u_shlf(I-1,J) + (u_shlf(I-1,J-1) + u_shlf(I-1,J+1)))) / (3*G%dxT(i,j)) + vx = ((v_shlf(I,J) + v_shlf(I,J-1) + v_shlf(I,J+1)) - & + (v_shlf(I-1,J) + (v_shlf(I-1,J-1) + v_shlf(I-1,J+1)))) / (3*G%dxT(i,j)) + uy = ((u_shlf(I,J) + (u_shlf(I-1,J) + u_shlf(I+1,J))) - & + (u_shlf(I,J-1) + (u_shlf(I-1,J-1) + u_shlf(I+1,J-1)))) / (3*G%dyT(i,j)) + vy = ((v_shlf(I,J) + (v_shlf(I-1,J)+ v_shlf(I+1,J))) - & + (v_shlf(I,J-1) + (v_shlf(I-1,J-1)+ v_shlf(I+1,J-1)))) / (3*G%dyT(i,j)) CS%ice_visc(i,j) = 0.5 * Visc_coef * (G%areaT(i,j) * ISS%h_shelf(i,j)) * & (US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2))**((1.-n_g)/(2.*n_g)) + endif + enddo + enddo +end subroutine calc_shelf_visc + +subroutine calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf) + type(ice_shelf_dyn_CS), intent(inout) :: CS !< A pointer to the ice shelf control structure + type(ice_shelf_state), intent(in) :: ISS !< A structure with elements that describe + !! the ice-shelf state + type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf. + type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors + real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), & + intent(inout) :: u_shlf !< The zonal ice shelf velocity [L T-1 ~> m s-1]. + real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), & + intent(inout) :: v_shlf !< The meridional ice shelf velocity [L T-1 ~> m s-1]. + +! also this subroutine updates the nonlinear part of the basal traction + +! this may be subject to change later... to make it "hybrid" + + integer :: i, j, iscq, iecq, jscq, jecq, isd, jsd, ied, jed, iegq, jegq + integer :: giec, gjec, gisc, gjsc, cnt, isc, jsc, iec, jec, is, js + real :: umid, vmid, unorm, eps_min ! Velocities [L T-1 ~> m s-1] + + isc = G%isc ; jsc = G%jsc ; iec = G%iec ; jec = G%jec + iscq = G%iscB ; iecq = G%iecB ; jscq = G%jscB ; jecq = G%jecB + isd = G%isd ; jsd = G%jsd ; ied = G%ied ; jed = G%jed + iegq = G%iegB ; jegq = G%jegB + gisc = G%domain%nihalo+1 ; gjsc = G%domain%njhalo+1 + giec = G%domain%niglobal+gisc ; gjec = G%domain%njglobal+gjsc + is = iscq - 1; js = jscq - 1 + + eps_min = CS%eps_glen_min + + + do j=jsd+1,jed + do i=isd+1,ied + + if ((ISS%hmask(i,j) == 1) .OR. (ISS%hmask(i,j) == 3)) then umid = ((u_shlf(I,J) + u_shlf(I-1,J-1)) + (u_shlf(I,J-1) + u_shlf(I-1,J))) * 0.25 vmid = ((v_shlf(I,J) + v_shlf(I-1,J-1)) + (v_shlf(I,J-1) + v_shlf(I-1,J))) * 0.25 unorm = sqrt(umid**2 + vmid**2 + eps_min**2*(G%dxT(i,j)**2 + G%dyT(i,j)**2)) @@ -2506,7 +2624,7 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) enddo enddo -end subroutine calc_shelf_visc +end subroutine calc_shelf_taub subroutine update_OD_ffrac(CS, G, US, ocean_mass, find_avg) type(ice_shelf_dyn_CS), intent(inout) :: CS !< A pointer to the ice shelf control structure @@ -2674,8 +2792,18 @@ subroutine bilinear_shape_fn_grid(G, i, j, Phi) xquad(2:4:2) = .5 * (1+sqrt(1./3)) ; yquad(3:4) = .5 * (1+sqrt(1./3)) do qpoint=1,4 - a = G%dxCv(i,J-1) * (1-yquad(qpoint)) + G%dxCv(i,J) * yquad(qpoint) ! d(x)/d(x*) - d = G%dyCu(I-1,j) * (1-xquad(qpoint)) + G%dyCu(I,j) * xquad(qpoint) ! d(y)/d(y*) + if (J>1) then + a = G%dxCv(i,J-1) * (1-yquad(qpoint)) + G%dxCv(i,J) * yquad(qpoint) ! d(x)/d(x*) + else + a= G%dxCv(i,J) !* yquad(qpoint) ! d(x)/d(x*) + endif + if (I>1) then + d = G%dyCu(I-1,j) * (1-xquad(qpoint)) + G%dyCu(I,j) * xquad(qpoint) ! d(y)/d(y*) + else + d = G%dyCu(I,j) !* xquad(qpoint) + endif +! a = G%dxCv(i,J-1) * (1-yquad(qpoint)) + G%dxCv(i,J) * yquad(qpoint) ! d(x)/d(x*) +! d = G%dyCu(I-1,j) * (1-xquad(qpoint)) + G%dyCu(I,j) * xquad(qpoint) ! d(y)/d(y*) do node=1,4 xnode = 2-mod(node,2) ; ynode = ceiling(REAL(node)/2) @@ -2799,16 +2927,18 @@ subroutine update_velocity_masks(CS, G, hmask, umask, vmask, u_face_mask, v_face if (hmask(i,j) == 1) then - umask(I-1:I,j-1:j) = 1. - vmask(I-1:I,j-1:j) = 1. + umask(I,j) = 1. + vmask(I,j) = 1. do k=0,1 select case (int(CS%u_face_mask_bdry(I-1+k,j))) case (3) - umask(I-1+k,J-1:J)=3. - vmask(I-1+k,J-1:J)=0. + ! vmask(I-1+k,J-1)=0. u_face_mask(I-1+k,j)=3. + umask(I-1+k,J)=3. + !vmask(I-1+k,J)=0. + vmask(I-1+k,J)=3. case (2) u_face_mask(I-1+k,j)=2. case (4) @@ -2829,8 +2959,10 @@ subroutine update_velocity_masks(CS, G, hmask, umask, vmask, u_face_mask, v_face select case (int(CS%v_face_mask_bdry(i,J-1+k))) case (3) - vmask(I-1:I,J-1+k)=3. - umask(I-1:I,J-1+k)=0. + vmask(I-1,J-1+k)=3. + umask(I-1,J-1+k)=0. + vmask(I,J-1+k)=3. + umask(I,J-1+k)=0. v_face_mask(i,J-1+k)=3. case (2) v_face_mask(i,J-1+k)=2. @@ -2848,21 +2980,6 @@ subroutine update_velocity_masks(CS, G, hmask, umask, vmask, u_face_mask, v_face end select enddo - !if (CS%u_face_mask_bdry(I-1,j) >= 0) then ! Western boundary - ! u_face_mask(I-1,j) = CS%u_face_mask_bdry(I-1,j) - ! umask(I-1,J-1:J) = 3. - ! vmask(I-1,J-1:J) = 0. - !endif - - !if (j_off+j == gjsc+1) then ! SoutherN boundary - ! v_face_mask(i,J-1) = 0. - ! umask(I-1:I,J-1) = 0. - ! vmask(I-1:I,J-1) = 0. - !elseif (j_off+j == gjec) then ! Northern boundary - ! v_face_mask(i,J) = 0. - ! umask(I-1:I,J) = 0. - ! vmask(I-1:I,J) = 0. - !endif if (i < G%ied) then if ((hmask(i+1,j) == 0) .OR. (hmask(i+1,j) == 2)) then @@ -2984,7 +3101,6 @@ subroutine ice_shelf_temp(CS, ISS, G, US, time_step, melt_rate, Time) intent(in) :: melt_rate !< basal melt rate [R Z T-1 ~> kg m-2 s-1] type(time_type), intent(in) :: Time !< The current model time -! 5/23/12 OVS ! This subroutine takes the velocity (on the Bgrid) and timesteps ! (HT)_t = - div (uHT) + (adot Tsurf -bdot Tbot) once and then calculates T=HT/H ! @@ -3020,12 +3136,6 @@ subroutine ice_shelf_temp(CS, ISS, G, US, time_step, melt_rate, Time) TH(i,j) = CS%t_shelf(i,j)*ISS%h_shelf(i,j) enddo ; enddo -! call enable_averages(time_step, Time, CS%diag) -! call pass_var(h_after_uflux, G%domain) -! call pass_var(h_after_vflux, G%domain) -! if (CS%id_h_after_uflux > 0) call post_data(CS%id_h_after_uflux, h_after_uflux, CS%diag) -! if (CS%id_h_after_vflux > 0) call post_data(CS%id_h_after_vflux, h_after_vflux, CS%diag) -! call disable_averaging(CS%diag) call ice_shelf_advect_temp_x(CS, G, time_step, ISS%hmask, TH, th_after_uflux) call ice_shelf_advect_temp_y(CS, G, time_step, ISS%hmask, th_after_uflux, th_after_vflux) diff --git a/src/ice_shelf/MOM_ice_shelf_initialize.F90 b/src/ice_shelf/MOM_ice_shelf_initialize.F90 index 9fe8028ac6..7a59d1586d 100644 --- a/src/ice_shelf/MOM_ice_shelf_initialize.F90 +++ b/src/ice_shelf/MOM_ice_shelf_initialize.F90 @@ -16,8 +16,9 @@ module MOM_ice_shelf_initialize #include -!MJHpublic initialize_ice_shelf_boundary, initialize_ice_thickness public initialize_ice_thickness +public initialize_ice_shelf_boundary_channel +public initialize_ice_flow_from_file ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with @@ -128,10 +129,6 @@ subroutine initialize_ice_thickness_from_file(h_shelf, area_shelf_h, hmask, G, U call MOM_read_data(filename, trim(thickness_varname), h_shelf, G%Domain, scale=US%m_to_Z) call MOM_read_data(filename,trim(area_varname), area_shelf_h, G%Domain, scale=US%m_to_L**2) -! call get_param(PF, mdl, "ICE_BOUNDARY_CONFIG", config, & -! "This specifies how the ice domain boundary is specified", & -! fail_if_missing=.true.) - isc = G%isc ; jsc = G%jsc ; iec = G%iec ; jec = G%jec do j=jsc,jec @@ -224,7 +221,6 @@ subroutine initialize_ice_thickness_channel(h_shelf, area_shelf_h, hmask, G, US, if (G%geoLonCu(i-1,j) >= edge_pos) then ! Everything past the edge is open ocean. -! mass_shelf(i,j) = 0.0 area_shelf_h(i,j) = 0.0 hmask (i,j) = 0.0 h_shelf (i,j) = 0.0 @@ -240,11 +236,7 @@ subroutine initialize_ice_thickness_channel(h_shelf, area_shelf_h, hmask, G, US, if (G%geoLonT(i,j) > slope_pos) then h_shelf(i,j) = min_draft -! mass_shelf(i,j) = Rho_ocean * min_draft else -! mass_shelf(i,j) = Rho_ocean * (min_draft + & -! (CS%max_draft - CS%min_draft) * & -! min(1.0, (c1*(slope_pos - G%geoLonT(i,j)))**2) ) h_shelf(i,j) = (min_draft + & (max_draft - min_draft) * & min(1.0, (c1*(slope_pos - G%geoLonT(i,j)))**2) ) @@ -262,165 +254,198 @@ subroutine initialize_ice_thickness_channel(h_shelf, area_shelf_h, hmask, G, US, end subroutine initialize_ice_thickness_channel -!BEGIN MJH -! subroutine initialize_ice_shelf_boundary(u_face_mask_bdry, v_face_mask_bdry, & -! u_flux_bdry_val, v_flux_bdry_val, u_bdry_val, v_bdry_val, h_bdry_val, & -! hmask, G, US, PF ) - -! type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure -! real, dimension(SZIB_(G),SZJ_(G)), & -! intent(inout) :: u_face_mask_bdry !< A boundary-type mask at C-grid u faces -! real, dimension(SZIB_(G),SZJ_(G)), & -! intent(inout) :: u_flux_bdry_val !< The boundary thickness flux through -! !! C-grid u faces [L Z T-1 ~> m2 s-1]. -! real, dimension(SZI_(G),SZJB_(G)), & -! intent(inout) :: v_face_mask_bdry !< A boundary-type mask at C-grid v faces -! real, dimension(SZI_(G),SZJB_(G)), & -! intent(inout) :: v_flux_bdry_val !< The boundary thickness flux through -! !! C-grid v faces [L Z T-1 ~> m2 s-1]. -! real, dimension(SZIB_(G),SZJB_(G)), & -! intent(inout) :: u_bdry_val !< The zonal ice shelf velocity at open -! !! boundary vertices [L T-1 ~> m s-1]. -! real, dimension(SZIB_(G),SZJB_(G)), & -! intent(inout) :: v_bdry_val !< The meridional ice shelf velocity at open -! !! boundary vertices [L T-1 ~> m s-1]. -! real, dimension(SZDI_(G),SZDJ_(G)), & -! intent(inout) :: h_bdry_val !< The ice shelf thickness at open boundaries [Z ~> m] -! real, dimension(SZDI_(G),SZDJ_(G)), & -! intent(inout) :: hmask !< A mask indicating which tracer points are -! !! partly or fully covered by an ice-shelf -! type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors -! type(param_file_type), intent(in) :: PF !< A structure to parse for run-time parameters - -! character(len=40) :: mdl = "initialize_ice_shelf_boundary" ! This subroutine's name. -! character(len=200) :: config -! logical flux_bdry - -! call get_param(PF, mdl, "ICE_BOUNDARY_CONFIG", config, & -! "This specifies how the ice domain boundary is specified. "//& -! "valid values include CHANNEL, FILE and USER.", & -! fail_if_missing=.true.) -! call get_param(PF, mdl, "ICE_BOUNDARY_FLUX_CONDITION", flux_bdry, & -! "This specifies whether mass input is a dirichlet or "//& -! "flux condition", default=.true.) - -! select case ( trim(config) ) -! case ("CHANNEL") -! call initialize_ice_shelf_boundary_channel(u_face_mask_bdry, & -! v_face_mask_bdry, u_flux_bdry_val, v_flux_bdry_val, & -! u_bdry_val, v_bdry_val, h_bdry_val, hmask, G, & -! flux_bdry, PF) -! case ("FILE"); call MOM_error(FATAL,"MOM_initialize: "// & -! "Unrecognized topography setup "//trim(config)) -! case ("USER"); call MOM_error(FATAL,"MOM_initialize: "// & -! "Unrecognized topography setup "//trim(config)) -! case default ; call MOM_error(FATAL,"MOM_initialize: "// & -! "Unrecognized topography setup "//trim(config)) -! end select - -! end subroutine initialize_ice_shelf_boundary - -! subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_bdry, & -! u_flux_bdry_val, v_flux_bdry_val, u_bdry_val, v_bdry_val, h_bdry_val, & -! hmask, G, flux_bdry, US, PF ) - -! type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure -! real, dimension(SZIB_(G),SZJ_(G)), & -! intent(inout) :: u_face_mask_bdry !< A boundary-type mask at C-grid u faces -! real, dimension(SZIB_(G),SZJ_(G)), & -! intent(inout) :: u_flux_bdry_val !< The boundary thickness flux through -! !! C-grid u faces [L Z T-1 ~> m2 s-1]. -! real, dimension(SZI_(G),SZJB_(G)), & -! intent(inout) :: v_face_mask_bdry !< A boundary-type mask at C-grid v faces -! real, dimension(SZI_(G),SZJB_(G)), & -! intent(inout) :: v_flux_bdry_val !< The boundary thickness flux through -! !! C-grid v faces [L Z T-1 ~> m2 s-1]. -! real, dimension(SZIB_(G),SZJB_(G)), & -! intent(inout) :: u_bdry_val !< The zonal ice shelf velocity at open - !! boundary vertices [L T-1 ~> m s-1]. -! real, dimension(SZIB_(G),SZJB_(G)), & -! intent(inout) :: v_bdry_val !< The meridional ice shelf velocity at open +!> Initialize ice shelf boundary conditions for a channel configuration +subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_bdry, & + u_flux_bdry_val, v_flux_bdry_val, u_bdry_val, v_bdry_val, u_shelf, v_shelf, h_bdry_val, & + thickness_bdry_val, hmask, h_shelf, G,& + US, PF ) + + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + real, dimension(SZIB_(G),SZJ_(G)), & + intent(inout) :: u_face_mask_bdry !< A boundary-type mask at C-grid u faces + real, dimension(SZIB_(G),SZJ_(G)), & + intent(inout) :: u_flux_bdry_val !< The boundary thickness flux through + !! C-grid u faces [L Z T-1 ~> m2 s-1]. + real, dimension(SZI_(G),SZJB_(G)), & + intent(inout) :: v_face_mask_bdry !< A boundary-type mask at C-grid v faces + real, dimension(SZI_(G),SZJB_(G)), & + intent(inout) :: v_flux_bdry_val !< The boundary thickness flux through + !! C-grid v faces [L Z T-1 ~> m2 s-1]. + real, dimension(SZIB_(G),SZJB_(G)), & + intent(inout) :: u_bdry_val !< The zonal ice shelf velocity at open !! boundary vertices [L T-1 ~> m s-1]. -! real, dimension(SZDI_(G),SZDJ_(G)), & -! intent(inout) :: h_bdry_val !< The ice shelf thickness at open boundaries [Z ~> m] -! real, dimension(SZDI_(G),SZDJ_(G)), & -! intent(inout) :: hmask !< A mask indicating which tracer points are -! !! partly or fully covered by an ice-shelf -! logical, intent(in) :: flux_bdry !< If true, use mass fluxes as the boundary value. -! type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors -! type(param_file_type), intent(in) :: PF !< A structure to parse for run-time parameters - -! character(len=40) :: mdl = "initialize_ice_shelf_boundary_channel" ! This subroutine's name. -! integer :: i, j, isd, jsd, is, js, iegq, jegq, giec, gjec, gisc, gjsc, isc, jsc, iec, jec, ied, jed -! real :: input_thick ! The input ice shelf thickness [Z ~> m] -! real :: input_flux ! The input ice flux per unit length [L Z T-1 ~> m2 s-1] -! real :: lenlat, len_stress - -! call get_param(PF, mdl, "LENLAT", lenlat, fail_if_missing=.true.) - -! call get_param(PF, mdl, "INPUT_FLUX_ICE_SHELF", input_flux, & -! "volume flux at upstream boundary", & -! units="m2 s-1", default=0., scale=US%m_s_to_L_T*US%m_to_Z) -! call get_param(PF, mdl, "INPUT_THICK_ICE_SHELF", input_thick, & -! "flux thickness at upstream boundary", & -! units="m", default=1000., scale=US%m_to_Z) -! call get_param(PF, mdl, "LEN_SIDE_STRESS", len_stress, & -! "maximum position of no-flow condition in along-flow direction", & -! units="km", default=0.) - -! call MOM_mesg(mdl//": setting boundary") - -! isd = G%isd ; ied = G%ied -! jsd = G%jsd ; jed = G%jed -! isc = G%isc ; jsc = G%jsc ; iec = G%iec ; jec = G%jec -! gisc = G%Domain%nihalo ; gjsc = G%Domain%njhalo -! giec = G%Domain%niglobal+gisc ; gjec = G%Domain%njglobal+gjsc - -! do j=jsd,jed -! do i=isd,ied - -! ! upstream boundary - set either dirichlet or flux condition - -! if ((i+G%idg_offset) == G%domain%nihalo+1) then -! if (flux_bdry) then -! u_face_mask_bdry(i-1,j) = 4.0 -! u_flux_bdry_val(i-1,j) = input_flux -! else -! hmask(i-1,j) = 3.0 -! h_bdry_val(i-1,j) = input_thick -! u_face_mask_bdry(i-1,j) = 3.0 -! u_bdry_val(i-1,j-1) = (1 - ((G%geoLatBu(i-1,j-1) - 0.5*lenlat)*2./lenlat)**2) * & -! 1.5 * input_flux / input_thick -! u_bdry_val(i-1,j) = (1 - ((G%geoLatBu(i-1,j) - 0.5*lenlat)*2./lenlat)**2) * & -! 1.5 * input_flux / input_thick -! endif -! endif - -! ! side boundaries: no flow - -! if (G%jdg_offset+j == gjsc+1) then !bot boundary -! if (len_stress == 0. .OR. G%geoLonCv(i,j-1) <= len_stress) then -! v_face_mask_bdry(i,j-1) = 0. -! else -! v_face_mask_bdry(i,j-1) = 1. -! endif -! elseif (G%jdg_offset+j == gjec) then !top boundary -! if (len_stress == 0. .OR. G%geoLonCv(i,j-1) <= len_stress) then -! v_face_mask_bdry(i,j) = 0. -! else -! v_face_mask_bdry(i,j) = 1. -! endif -! endif - -! ! downstream boundary - CFBC - -! if (i+G%idg_offset == giec) then -! u_face_mask_bdry(i,j) = 2.0 -! endif - -! enddo -! enddo - -!END MJH end subroutine initialize_ice_shelf_boundary_channel + real, dimension(SZIB_(G),SZJB_(G)), & + intent(inout) :: v_bdry_val !< The meridional ice shelf velocity at open + real, dimension(SZIB_(G),SZJB_(G)), & + intent(inout) :: u_shelf !< The zonal ice shelf velocity [L T-1 ~> m s-1]. + real, dimension(SZIB_(G),SZJB_(G)), & + intent(inout) :: v_shelf !< The meridional ice shelf velocity [L T-1 ~> m s-1]. + real, dimension(SZDI_(G),SZDJ_(G)), & + intent(inout) :: thickness_bdry_val !< The ice shelf thickness at open boundaries [Z ~> m] + !! boundary vertices [L T-1 ~> m s-1]. + real, dimension(SZDI_(G),SZDJ_(G)), & + intent(inout) :: h_bdry_val !< The ice shelf thickness at open boundaries [Z ~> m] + real, dimension(SZDI_(G),SZDJ_(G)), & + intent(inout) :: hmask !< A mask indicating which tracer points are + !! partly or fully covered by an ice-shelf + real, dimension(SZDI_(G),SZDJ_(G)), & + intent(inout) :: h_shelf !< Ice-shelf thickness + type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors + type(param_file_type), intent(in) :: PF !< A structure to parse for run-time parameters + + character(len=40) :: mdl = "initialize_ice_shelf_boundary_channel" ! This subroutine's name. + integer :: i, j, isd, jsd, is, js, iegq, jegq, giec, gjec, gisc, gjsc,gisd,gjsd, isc, jsc, iec, jec, ied, jed + real :: input_thick ! The input ice shelf thickness [Z ~> m] + real :: input_vel ! The input ice velocity per [L Z T-1 ~> m s-1] + real :: lenlat, len_stress, westlon, lenlon, southlat ! The input positions of the channel boundarises + + + call get_param(PF, mdl, "LENLAT", lenlat, fail_if_missing=.true.) + + call get_param(PF, mdl, "LENLON", lenlon, fail_if_missing=.true.) + + call get_param(PF, mdl, "WESTLON", westlon, fail_if_missing=.true.) + + call get_param(PF, mdl, "SOUTHLAT", southlat, fail_if_missing=.true.) + + call get_param(PF, mdl, "INPUT_VEL_ICE_SHELF", input_vel, & + "inflow ice velocity at upstream boundary", & + units="m s-1", default=0., scale=US%m_s_to_L_T*US%m_to_Z) + call get_param(PF, mdl, "INPUT_THICK_ICE_SHELF", input_thick, & + "flux thickness at upstream boundary", & + units="m", default=1000., scale=US%m_to_Z) + call get_param(PF, mdl, "LEN_SIDE_STRESS", len_stress, & + "maximum position of no-flow condition in along-flow direction", & + units="km", default=0.) + + call MOM_mesg(mdl//": setting boundary") + + isd = G%isd ; ied = G%ied + jsd = G%jsd ; jed = G%jed + isc = G%isc ; jsc = G%jsc ; iec = G%iec ; jec = G%jec + gjsd = G%Domain%njglobal ; gisd = G%Domain%niglobal + gisc = G%Domain%nihalo ; gjsc = G%Domain%njhalo + giec = G%Domain%niglobal+gisc ; gjec = G%Domain%njglobal+gjsc + + !---------b.c.s based on geopositions ----------------- + do j=jsc,jec+1 + do i=isc-1,iec+1 + ! upstream boundary - set either dirichlet or flux condition + + if (G%geoLonBu(i,j) == westlon) then + hmask(i+1,j) = 3.0 + h_bdry_val(i+1,j) = h_shelf(i+1,j) + thickness_bdry_val(i+1,j) = h_bdry_val(i+0*1,j) + u_face_mask_bdry(i+1,j) = 3.0 + u_bdry_val(i+1,j) = input_vel*(1-16.0*((G%geoLatBu(i-1,j)/lenlat-0.5))**4) !velocity distribution + endif + + + ! side boundaries: no flow + if (G%geoLatBu(i,j-1) == southlat) then !bot boundary + if (len_stress == 0. .OR. G%geoLonCv(i,j) <= len_stress) then + v_face_mask_bdry(i,j+1) = 0. + u_face_mask_bdry(i,j) = 3. + u_bdry_val(i,j) = 0. + v_bdry_val(i,j) = 0. + else + v_face_mask_bdry(i,j+1) = 1. + u_face_mask_bdry(i,j) = 3. + u_bdry_val(i,j) = 0. + v_bdry_val(i,j) = 0. + endif + elseif (G%geoLatBu(i,j-1) == southlat+lenlat) then !top boundary + if (len_stress == 0. .OR. G%geoLonCv(i,j) <= len_stress) then + v_face_mask_bdry(i,j-1) = 0. + u_face_mask_bdry(i,j-1) = 3. + else + v_face_mask_bdry(i,j-1) = 3. + u_face_mask_bdry(i,j-1) = 3. + endif + endif + + ! downstream boundary - CFBC + if (G%geoLonBu(i,j) == westlon+lenlon) then + u_face_mask_bdry(i-1,j) = 2.0 + endif + + enddo + enddo +end subroutine initialize_ice_shelf_boundary_channel + + +!> Initialize ice shelf flow from file +subroutine initialize_ice_flow_from_file(u_shelf, v_shelf,ice_visc,float_cond,& + hmask,h_shelf, G, US, PF) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + real, dimension(SZDI_(G),SZDJ_(G)), & + intent(inout) :: u_shelf !< The ice shelf u velocity [Z ~> m T ~>s]. + real, dimension(SZDI_(G),SZDJ_(G)), & + intent(inout) :: v_shelf !< The ice shelf v velocity [Z ~> m T ~> s]. + real, dimension(SZDI_(G),SZDJ_(G)), & + intent(inout) :: ice_visc !< The ice shelf viscosity [Pa ~> m T ~> s]. + real, dimension(SZDI_(G),SZDJ_(G)), & + intent(inout) :: float_cond !< An array indicating where the ice + !! shelf is floating: 0 if floating, 1 if not. + real, dimension(SZDI_(G),SZDJ_(G)), & + intent(in) :: hmask !< A mask indicating which tracer points are + !! partly or fully covered by an ice-shelf + real, dimension(SZDI_(G),SZDJ_(G)), & + intent(in) :: h_shelf !< A mask indicating which tracer points are + !! partly or fully covered by an ice-shelf + type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors + type(param_file_type), intent(in) :: PF !< A structure to parse for run-time parameters + + ! This subroutine reads ice thickness and area from a file and puts it into + ! h_shelf [Z ~> m] and area_shelf_h [L2 ~> m2] (and dimensionless) and updates hmask + character(len=200) :: filename,vel_file,inputdir ! Strings for file/path + character(len=200) :: ushelf_varname, vshelf_varname, ice_visc_varname, floatfr_varname ! Variable name in file + character(len=40) :: mdl = "initialize_ice_velocity_from_file" ! This subroutine's name. + integer :: i, j, isc, jsc, iec, jec + real :: len_sidestress, mask, udh + + call MOM_mesg(" MOM_ice_shelf_init_profile.F90, initialize_velocity_from_file: reading velocity") + + call get_param(PF, mdl, "INPUTDIR", inputdir, default=".") + inputdir = slasher(inputdir) + call get_param(PF, mdl, "ICE_VELOCITY_FILE", vel_file, & + "The file from which the velocity is read.", & + default="ice_shelf_vel.nc") + call get_param(PF, mdl, "LEN_SIDE_STRESS", len_sidestress, & + "position past which shelf sides are stress free.", & + default=0.0, units="axis_units") + + filename = trim(inputdir)//trim(vel_file) + call log_param(PF, mdl, "INPUTDIR/THICKNESS_FILE", filename) + call get_param(PF, mdl, "ICE_U_VEL_VARNAME", ushelf_varname, & + "The name of the thickness variable in ICE_VELOCITY_FILE.", & + default="u_shelf") + call get_param(PF, mdl, "ICE_V_VEL_VARNAME", vshelf_varname, & + "The name of the thickness variable in ICE_VELOCITY_FILE.", & + default="v_shelf") + call get_param(PF, mdl, "ICE_VISC_VARNAME", ice_visc_varname, & + "The name of the thickness variable in ICE_VELOCITY_FILE.", & + default="viscosity") + + if (.not.file_exists(filename, G%Domain)) call MOM_error(FATAL, & + " initialize_ice_shelf_velocity_from_file: Unable to open "//trim(filename)) + + floatfr_varname = "float_frac" + + call MOM_read_data(filename, trim(ushelf_varname), u_shelf, G%Domain, scale=1.0) + call MOM_read_data(filename,trim(vshelf_varname), v_shelf, G%Domain, scale=1.0) + call MOM_read_data(filename,trim(ice_visc_varname), ice_visc, G%Domain, scale=1.0) + call MOM_read_data(filename,trim(floatfr_varname), float_cond, G%Domain, scale=1.) + + isc = G%isc ; jsc = G%jsc ; iec = G%iec ; jec = G%jec + + do j=jsc,jec + do i=isc,iec + if (hmask(i,j) == 1.) then + ice_visc(i,j) = ice_visc(i,j) * (G%areaT(i,j) * h_shelf(i,j)) + endif + enddo + enddo +end subroutine initialize_ice_flow_from_file end module MOM_ice_shelf_initialize diff --git a/src/initialization/MOM_shared_initialization.F90 b/src/initialization/MOM_shared_initialization.F90 index 70ef0768d5..ee80bbdace 100644 --- a/src/initialization/MOM_shared_initialization.F90 +++ b/src/initialization/MOM_shared_initialization.F90 @@ -196,6 +196,7 @@ subroutine apply_topography_edits_from_file(D, G, param_file, US) character(len=40) :: mdl = "apply_topography_edits_from_file" ! This subroutine's name. integer :: i, j, n, ncid, n_edits, i_file, j_file, ndims, sizes(8) logical :: found + logical :: topo_edits_change_mask call callTree_enter(trim(mdl)//"(), MOM_shared_initialization.F90") @@ -206,6 +207,9 @@ subroutine apply_topography_edits_from_file(D, G, param_file, US) call get_param(param_file, mdl, "TOPO_EDITS_FILE", topo_edits_file, & "The file from which to read a list of i,j,z topography overrides.", & default="") + call get_param(param_file, mdl, "ALLOW_LANDMASK_CHANGES", topo_edits_change_mask, & + "If true, allow topography overrides to change land mask.", & + default=.false.) if (len_trim(topo_edits_file)==0) return @@ -250,8 +254,14 @@ subroutine apply_topography_edits_from_file(D, G, param_file, US) '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 - call MOM_error(FATAL, trim(mdl)//': A zero depth edit would change the land mask and '//& - "is not allowed in"//trim(topo_edits_file)) + if (topo_edits_change_mask) then + 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 + call MOM_error(FATAL, ' apply_topography_edits_from_file: '//& + "A zero depth edit would change the land mask and is not allowed in"//trim(topo_edits_file)) + endif endif endif enddo diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index fc6d97040a..5bbc495e93 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -188,6 +188,7 @@ module MOM_hor_visc integer :: id_diffu = -1, id_diffv = -1 ! integer :: id_hf_diffu = -1, id_hf_diffv = -1 integer :: id_hf_diffu_2d = -1, id_hf_diffv_2d = -1 + integer :: id_intz_diffu_2d = -1, id_intz_diffv_2d = -1 integer :: id_Ah_h = -1, id_Ah_q = -1 integer :: id_Kh_h = -1, id_Kh_q = -1 integer :: id_GME_coeff_h = -1, id_GME_coeff_q = -1 @@ -276,8 +277,9 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, grad_d2vel_mag_h, & ! Magnitude of the Laplacian of the velocity vector, squared [L-2 T-2 ~> m-2 s-2] boundary_mask_h ! A mask that zeroes out cells with at least one land edge [nondim] - real, allocatable, dimension(:,:) :: hf_diffu_2d ! Depth sum of hf_diffu [L T-2 ~> m s-2] - real, allocatable, dimension(:,:) :: hf_diffv_2d ! Depth sum of hf_diffv [L T-2 ~> m s-2] + real, allocatable, dimension(:,:) :: hf_diffu_2d, hf_diffv_2d ! Depth sum of hf_diffu, hf_diffv [L T-2 ~> m s-2] + real, dimension(SZIB_(G),SZJ_(G)) :: intz_diffu_2d ! Depth-integral of diffu [L2 T-2 ~> m2 s-2] + real, dimension(SZI_(G),SZJB_(G)) :: intz_diffv_2d ! Depth-integral of diffv [L2 T-2 ~> m2 s-2] real, dimension(SZIB_(G),SZJB_(G)) :: & dvdx, dudy, & ! components in the shearing strain [T-1 ~> s-1] @@ -1662,6 +1664,21 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, deallocate(hf_diffv_2d) endif + if (present(ADp) .and. (CS%id_intz_diffu_2d > 0)) then + intz_diffu_2d(:,:) = 0.0 + do k=1,nz ; do j=js,je ; do I=Isq,Ieq + intz_diffu_2d(I,j) = intz_diffu_2d(I,j) + diffu(I,j,k) * ADp%diag_hu(I,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_intz_diffu_2d, intz_diffu_2d, CS%diag) + endif + if (present(ADp) .and. (CS%id_intz_diffv_2d > 0)) then + intz_diffv_2d(:,:) = 0.0 + do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + intz_diffv_2d(i,J) = intz_diffv_2d(i,J) + diffv(i,J,k) * ADp%diag_hv(i,J,k) + enddo ; enddo ; enddo + call post_data(CS%id_intz_diffv_2d, intz_diffv_2d, CS%diag) + endif + end subroutine horizontal_viscosity !> Allocates space for and calculates static variables used by horizontal_viscosity(). @@ -2378,6 +2395,20 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, MEKE, ADp) call safe_alloc_ptr(ADp%diag_hfrac_v,G%isd,G%ied,G%JsdB,G%JedB,GV%ke) endif + CS%id_intz_diffu_2d = register_diag_field('ocean_model', 'intz_diffu_2d', diag%axesCu1, Time, & + 'Depth-integral of Zonal Acceleration from Horizontal Viscosity', 'm2 s-2', & + conversion=GV%H_to_m*US%L_T2_to_m_s2) + if ((CS%id_intz_diffu_2d > 0) .and. (present(ADp))) then + call safe_alloc_ptr(ADp%diag_hu,G%IsdB,G%IedB,G%jsd,G%jed,GV%ke) + endif + + CS%id_intz_diffv_2d = register_diag_field('ocean_model', 'intz_diffv_2d', diag%axesCv1, Time, & + 'Depth-integral of Meridional Acceleration from Horizontal Viscosity', 'm2 s-2', & + conversion=GV%H_to_m*US%L_T2_to_m_s2) + if ((CS%id_intz_diffv_2d > 0) .and. (present(ADp))) then + call safe_alloc_ptr(ADp%diag_hv,G%isd,G%ied,G%JsdB,G%JedB,GV%ke) + endif + if (CS%biharmonic) then CS%id_Ah_h = register_diag_field('ocean_model', 'Ahh', diag%axesTL, Time, & 'Biharmonic Horizontal Viscosity at h Points', 'm4 s-1', conversion=US%L_to_m**4*US%s_to_T, & diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index e3a6f1599e..7d95c43b98 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -1125,16 +1125,18 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) if (CS%calculate_Eady_growth_rate .and. CS%use_stored_slopes) then CS%id_N2_u = register_diag_field('ocean_model', 'N2_u', diag%axesCui, Time, & 'Square of Brunt-Vaisala frequency, N^2, at u-points, as used in Visbeck et al.', & - 's-2', conversion=US%s_to_T**2) + 's-2', conversion=(US%L_to_Z*US%s_to_T)**2) CS%id_N2_v = register_diag_field('ocean_model', 'N2_v', diag%axesCvi, Time, & 'Square of Brunt-Vaisala frequency, N^2, at v-points, as used in Visbeck et al.', & - 's-2', conversion=US%s_to_T**2) + 's-2', conversion=(US%L_to_Z*US%s_to_T)**2) endif if (CS%use_stored_slopes) then CS%id_S2_u = register_diag_field('ocean_model', 'S2_u', diag%axesCu1, Time, & - 'Depth average square of slope magnitude, S^2, at u-points, as used in Visbeck et al.', 'nondim') + 'Depth average square of slope magnitude, S^2, at u-points, as used in Visbeck et al.', & + 'nondim', conversion=US%Z_to_L**2) CS%id_S2_v = register_diag_field('ocean_model', 'S2_v', diag%axesCv1, Time, & - 'Depth average square of slope magnitude, S^2, at v-points, as used in Visbeck et al.', 'nondim') + 'Depth average square of slope magnitude, S^2, at v-points, as used in Visbeck et al.', & + 'nondim', conversion=US%Z_to_L**2) endif oneOrTwo = 1.0 diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index 8c6a90ba9c..02a49a2a1a 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -40,7 +40,7 @@ module MOM_thickness_diffuse real :: max_Khth_CFL !< Maximum value of the diffusive CFL for thickness diffusion real :: Khth_Min !< Minimum value of Khth [L2 T-1 ~> m2 s-1] real :: Khth_Max !< Maximum value of Khth [L2 T-1 ~> m2 s-1], or 0 for no max - real :: slope_max !< Slopes steeper than slope_max are limited in some way [nondim]. + real :: slope_max !< Slopes steeper than slope_max are limited in some way [Z L-1 ~> nondim]. real :: kappa_smooth !< Vertical diffusivity used to interpolate more !! sensible values of T & S into thin layers [Z2 T-1 ~> m2 s-1]. logical :: thickness_diffuse !< If true, interfaces heights are diffused. @@ -83,8 +83,8 @@ module MOM_thickness_diffuse type(diag_ctrl), pointer :: diag => NULL() !< structure used to regulate timing of diagnostics real, pointer :: GMwork(:,:) => NULL() !< Work by thickness diffusivity [R Z L2 T-3 ~> W m-2] - real, pointer :: diagSlopeX(:,:,:) => NULL() !< Diagnostic: zonal neutral slope [nondim] - real, pointer :: diagSlopeY(:,:,:) => NULL() !< Diagnostic: zonal neutral slope [nondim] + real, pointer :: diagSlopeX(:,:,:) => NULL() !< Diagnostic: zonal neutral slope [Z L-1 ~> nondim] + real, pointer :: diagSlopeY(:,:,:) => NULL() !< Diagnostic: zonal neutral slope [Z L-1 ~> nondim] real, dimension(:,:,:), pointer :: & KH_u_GME => NULL(), & !< interface height diffusivities in u-columns [L2 T-1 ~> m2 s-1] @@ -578,8 +578,8 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV !! the isopycnal slopes are taken directly from !! the interface slopes without consideration of !! density gradients [nondim]. - real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), optional, intent(in) :: slope_x !< Isopycnal slope at u-points - real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), optional, intent(in) :: slope_y !< Isopycnal slope at v-points + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), optional, intent(in) :: slope_x !< Isopyc. slope at u [Z L-1 ~> nondim] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), optional, intent(in) :: slope_y !< Isopyc. slope at v [Z L-1 ~> nondim] ! Local variables real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: & T, & ! The temperature (or density) [degC], with the values in @@ -660,7 +660,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV real :: Slope ! The slope of density surfaces, calculated in a way ! that is always between -1 and 1, nondimensional. real :: mag_grad2 ! The squared magnitude of the 3-d density gradient [R2 L-2 ~> kg2 m-8]. - real :: I_slope_max2 ! The inverse of slope_max squared, nondimensional. + real :: I_slope_max2 ! The inverse of slope_max squared [L2 Z-2 ~> nondim]. real :: h_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected [H ~> m or kg m-2]. real :: h_neglect2 ! h_neglect^2 [H2 ~> m2 or kg2 m-4]. @@ -919,7 +919,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV ! This estimate of slope is accurate for small slopes, but bounded ! to be between -1 and 1. - mag_grad2 = drdx**2 + (US%L_to_Z*drdz)**2 + mag_grad2 = (US%Z_to_L*drdx)**2 + drdz**2 if (mag_grad2 > 0.0) then Slope = drdx / sqrt(mag_grad2) slope2_Ratio_u(I,K) = Slope**2 * I_slope_max2 @@ -933,7 +933,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV ! that ignore density gradients along layers. if (present_int_slope_u) then Slope = (1.0 - int_slope_u(I,j,K)) * Slope + & - int_slope_u(I,j,K) * US%Z_to_L*((e(i+1,j,K)-e(i,j,K)) * G%IdxCu(I,j)) + int_slope_u(I,j,K) * ((e(i+1,j,K)-e(i,j,K)) * G%IdxCu(I,j)) slope2_Ratio_u(I,K) = (1.0 - int_slope_u(I,j,K)) * slope2_Ratio_u(I,K) endif @@ -942,7 +942,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV if (CS%id_slope_x > 0) CS%diagSlopeX(I,j,k) = Slope ! Estimate the streamfunction at each interface [Z L2 T-1 ~> m3 s-1]. - Sfn_unlim_u(I,K) = -((KH_u(I,j,K)*G%dy_Cu(I,j))*US%L_to_Z*Slope) + Sfn_unlim_u(I,K) = -((KH_u(I,j,K)*G%dy_Cu(I,j))*Slope) ! Avoid moving dense water upslope from below the level of ! the bottom on the receiving side. @@ -968,10 +968,10 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV if (present_slope_x) then Slope = slope_x(I,j,k) else - Slope = US%Z_to_L*((e(i,j,K)-e(i+1,j,K))*G%IdxCu(I,j)) * G%mask2dCu(I,j) + Slope = ((e(i,j,K)-e(i+1,j,K))*G%IdxCu(I,j)) * G%mask2dCu(I,j) endif if (CS%id_slope_x > 0) CS%diagSlopeX(I,j,k) = Slope - Sfn_unlim_u(I,K) = ((KH_u(I,j,K)*G%dy_Cu(I,j))*US%L_to_Z*Slope) + Sfn_unlim_u(I,K) = ((KH_u(I,j,K)*G%dy_Cu(I,j))*Slope) hN2_u(I,K) = GV%g_prime(K) endif ! if (use_EOS) else ! if (k > nk_linear) @@ -993,10 +993,14 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV ! Solve an elliptic equation for the streamfunction following Ferrari et al., 2010. do I=is-1,ie if (G%mask2dCu(I,j)>0.) then - Sfn_unlim_u(I,:) = ( 1. + CS%FGNV_scale ) * Sfn_unlim_u(I,:) + do K=2,nz + Sfn_unlim_u(I,K) = (1. + CS%FGNV_scale) * Sfn_unlim_u(I,K) + enddo call streamfn_solver(nz, c2_h_u(I,:), hN2_u(I,:), Sfn_unlim_u(I,:)) else - Sfn_unlim_u(I,:) = 0. + do K=2,nz + Sfn_unlim_u(I,K) = 0. + enddo endif enddo endif @@ -1185,7 +1189,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV ! This estimate of slope is accurate for small slopes, but bounded ! to be between -1 and 1. - mag_grad2 = drdy**2 + (US%L_to_Z*drdz)**2 + mag_grad2 = (US%Z_to_L*drdy)**2 + drdz**2 if (mag_grad2 > 0.0) then Slope = drdy / sqrt(mag_grad2) slope2_Ratio_v(i,K) = Slope**2 * I_slope_max2 @@ -1199,7 +1203,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV ! that ignore density gradients along layers. if (present_int_slope_v) then Slope = (1.0 - int_slope_v(i,J,K)) * Slope + & - int_slope_v(i,J,K) * US%Z_to_L*((e(i,j+1,K)-e(i,j,K)) * G%IdyCv(i,J)) + int_slope_v(i,J,K) * ((e(i,j+1,K)-e(i,j,K)) * G%IdyCv(i,J)) slope2_Ratio_v(i,K) = (1.0 - int_slope_v(i,J,K)) * slope2_Ratio_v(i,K) endif @@ -1208,7 +1212,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV if (CS%id_slope_y > 0) CS%diagSlopeY(I,j,k) = Slope ! Estimate the streamfunction at each interface [Z L2 T-1 ~> m3 s-1]. - Sfn_unlim_v(i,K) = -((KH_v(i,J,K)*G%dx_Cv(i,J))*US%L_to_Z*Slope) + Sfn_unlim_v(i,K) = -((KH_v(i,J,K)*G%dx_Cv(i,J))*Slope) ! Avoid moving dense water upslope from below the level of ! the bottom on the receiving side. @@ -1234,10 +1238,10 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV if (present_slope_y) then Slope = slope_y(i,J,k) else - Slope = US%Z_to_L*((e(i,j,K)-e(i,j+1,K))*G%IdyCv(i,J)) * G%mask2dCv(i,J) + Slope = ((e(i,j,K)-e(i,j+1,K))*G%IdyCv(i,J)) * G%mask2dCv(i,J) endif if (CS%id_slope_y > 0) CS%diagSlopeY(I,j,k) = Slope - Sfn_unlim_v(i,K) = ((KH_v(i,J,K)*G%dx_Cv(i,J))*US%L_to_Z*Slope) + Sfn_unlim_v(i,K) = ((KH_v(i,J,K)*G%dx_Cv(i,J))*Slope) hN2_v(i,K) = GV%g_prime(K) endif ! if (use_EOS) else ! if (k > nk_linear) @@ -1259,10 +1263,14 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV ! Solve an elliptic equation for the streamfunction following Ferrari et al., 2010. do i=is,ie if (G%mask2dCv(i,J)>0.) then - Sfn_unlim_v(i,:) = ( 1. + CS%FGNV_scale ) * Sfn_unlim_v(i,:) + do K=2,nz + Sfn_unlim_v(i,K) = (1. + CS%FGNV_scale) * Sfn_unlim_v(i,K) + enddo call streamfn_solver(nz, c2_h_v(i,:), hN2_v(i,:), Sfn_unlim_v(i,:)) else - Sfn_unlim_v(i,:) = 0. + do K=2,nz + Sfn_unlim_v(i,K) = 0. + enddo endif enddo endif @@ -1947,7 +1955,7 @@ subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS) "longer than DT, or 0 to use DT.", units="s", default=0.0, scale=US%s_to_T) call get_param(param_file, mdl, "KHTH_SLOPE_MAX", CS%slope_max, & "A slope beyond which the calculated isopycnal slope is "//& - "not reliable and is scaled away.", units="nondim", default=0.01) + "not reliable and is scaled away.", units="nondim", default=0.01, scale=US%L_to_Z) call get_param(param_file, mdl, "KD_SMOOTH", CS%kappa_smooth, & "A diapycnal diffusivity that is used to interpolate "//& "more sensible values of T & S into thin layers.", & @@ -2065,10 +2073,10 @@ subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS) 'm2 s-1', conversion=US%L_to_m**2*US%s_to_T) CS%id_slope_x = register_diag_field('ocean_model', 'neutral_slope_x', diag%axesCui, Time, & - 'Zonal slope of neutral surface', 'nondim') + 'Zonal slope of neutral surface', 'nondim', conversion=US%Z_to_L) if (CS%id_slope_x > 0) call safe_alloc_ptr(CS%diagSlopeX,G%IsdB,G%IedB,G%jsd,G%jed,GV%ke+1) CS%id_slope_y = register_diag_field('ocean_model', 'neutral_slope_y', diag%axesCvi, Time, & - 'Meridional slope of neutral surface', 'nondim') + 'Meridional slope of neutral surface', 'nondim', conversion=US%Z_to_L) if (CS%id_slope_y > 0) call safe_alloc_ptr(CS%diagSlopeY,G%isd,G%ied,G%JsdB,G%JedB,GV%ke+1) CS%id_sfn_x = register_diag_field('ocean_model', 'GM_sfn_x', diag%axesCui, Time, & 'Parameterized Zonal Overturning Streamfunction', & diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 649fc725de..c96edb785c 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -2982,7 +2982,7 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di CS%id_hf_dvdt_dia_2d = register_diag_field('ocean_model', 'hf_dvdt_dia_2d', diag%axesCv1, Time, & 'Depth-sum Fractional Thickness-weighted Meridional Acceleration from Diapycnal Mixing', & 'm s-2', conversion=US%L_T2_to_m_s2) - if (CS%id_hf_dvdt_dia_2d > 0) call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz) + if (CS%id_hf_dvdt_dia_2d > 0) call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz) if ((CS%id_dudt_dia > 0) .or. (CS%id_hf_dudt_dia_2d > 0)) & call safe_alloc_ptr(ADp%du_dt_dia,IsdB,IedB,jsd,jed,nz) diff --git a/src/parameterizations/vertical/MOM_tidal_mixing.F90 b/src/parameterizations/vertical/MOM_tidal_mixing.F90 index b870dff1af..512179445b 100644 --- a/src/parameterizations/vertical/MOM_tidal_mixing.F90 +++ b/src/parameterizations/vertical/MOM_tidal_mixing.F90 @@ -882,7 +882,6 @@ subroutine calculate_CVMix_tidal(h, j, G, GV, US, CS, N2_int, Kd_lay, Kd_int, Kv ! summing q_i*TidalConstituent_i over the number of constituents. call CVMix_compute_SchmittnerCoeff( nlev = GV%ke, & energy_flux = tidal_qe_md(:), & - rho = rho_fw, & SchmittnerCoeff = Schmittner_coeff, & exp_hab_zetar = exp_hab_zetar, & CVmix_tidal_params_user = CS%CVMix_tidal_params) @@ -896,7 +895,6 @@ subroutine calculate_CVMix_tidal(h, j, G, GV, US, CS, N2_int, Kd_lay, Kd_int, Kv Tdiff_out = Kd_tidal, & Nsqr = N2_int_i, & OceanDepth = -iFaceHeight(GV%ke+1), & - vert_dep = vert_dep, & nlev = GV%ke, & max_nlev = GV%ke, & SchmittnerCoeff = Schmittner_coeff, & diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index 7e2c2c6926..c7c9355f97 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -1817,7 +1817,7 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & !if (CS%id_hf_dv_dt_visc > 0) then ! call safe_alloc_ptr(CS%hf_dv_dt_visc,isd,ied,JsdB,JedB,nz) ! call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz) - ! call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz) + ! call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz) !endif CS%id_hf_du_dt_visc_2d = register_diag_field('ocean_model', 'hf_du_dt_visc_2d', diag%axesCu1, Time, & @@ -1833,7 +1833,7 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & conversion=US%L_T2_to_m_s2) if (CS%id_hf_dv_dt_visc_2d > 0) then call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz) - call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz) + call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz) endif if ((len_trim(CS%u_trunc_file) > 0) .or. (len_trim(CS%v_trunc_file) > 0)) &