diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 64ee27ab7d..fe44a68e0b 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -691,7 +691,7 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) if (CS%id_T_preale > 0) call post_data(CS%id_T_preale, CS%tv%T, CS%diag) if (CS%id_S_preale > 0) call post_data(CS%id_S_preale, CS%tv%S, CS%diag) if (CS%id_e_preale > 0) then - call find_eta(h, CS%tv, G%g_Earth, G, GV, eta_preale) + call find_eta(h, CS%tv, GV%g_Earth, G, GV, eta_preale) call post_data(CS%id_e_preale, eta_preale, CS%diag) endif @@ -985,7 +985,7 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) if (CS%id_T_predia > 0) call post_data(CS%id_T_predia, CS%tv%T, CS%diag) if (CS%id_S_predia > 0) call post_data(CS%id_S_predia, CS%tv%S, CS%diag) if (CS%id_e_predia > 0) then - call find_eta(h, CS%tv, G%g_Earth, G, GV, eta_predia) + call find_eta(h, CS%tv, GV%g_Earth, G, GV, eta_predia) call post_data(CS%id_e_predia, eta_predia, CS%diag) endif @@ -1185,7 +1185,7 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) ! diagnosed ssh for non-Bouss; call "find_eta" for this ! purpose. tot_wt_ssh = tot_wt_ssh + dt - call find_eta(h, CS%tv, G%g_Earth, G, GV, ssh, eta_av) + call find_eta(h, CS%tv, GV%g_Earth, G, GV, ssh, eta_av) do j=js,je ; do i=is,ie CS%ave_ssh(i,j) = CS%ave_ssh(i,j) + dt*ssh(i,j) enddo ; enddo @@ -1206,13 +1206,13 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) CS%ave_ssh(i,j) = CS%ave_ssh(i,j)*Itot_wt_ssh enddo ; enddo + call enable_averaging(dt*n_max,Time_local, CS%diag) ! area mean SSH if (CS%id_ssh_ga > 0) then ssh_ga = global_area_mean(CS%ave_ssh, G) call post_data(CS%id_ssh_ga, ssh_ga, CS%diag) endif - call enable_averaging(dt*n_max,Time_local, CS%diag) I_time_int = 1.0/(dt*n_max) if (CS%id_ssh > 0) & call post_data(CS%id_ssh, CS%ave_ssh, CS%diag, mask=G%mask2dT) @@ -1228,7 +1228,7 @@ subroutine step_MOM(fluxes, state, Time_start, time_interval, CS) enddo ; enddo if (ASSOCIATED(fluxes%p_surf)) then do j=js,je ; do i=is,ie - zos(i,j) = zos(i,j)+G%mask2dT(i,j)*fluxes%p_surf(i,j)/(GV%Rho0 * G%g_Earth) + zos(i,j) = zos(i,j)+G%mask2dT(i,j)*fluxes%p_surf(i,j)/(GV%Rho0 * GV%g_Earth) enddo ; enddo endif zos_area_mean = global_area_mean(zos, G) @@ -1445,6 +1445,7 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) default=.false.) CS%legacy_split = .false. endif + call get_param(param_file, "MOM", "ENABLE_THERMODYNAMICS", CS%use_temperature, & "If true, Temperature and salinity are used as state \n"//& "variables.", default=.true.) @@ -1462,27 +1463,22 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) "true. This assumes that KD = KDML = 0.0 and that \n"//& "there is no buoyancy forcing, but makes the model \n"//& "faster by eliminating subroutine calls.", default=.false.) + call get_param(param_file, "MOM", "DO_DYNAMICS", CS%do_dynamics, & + "If False, skips the dynamics calls that update u & v, as well as\n"//& + "the gravity wave adjustment to h. This is a fragile feature and\n"//& + "thus undocumented.", default=.true., do_not_log=.true. ) + + call get_param(param_file, "MOM", "USE_REGRIDDING", CS%use_ALE_algorithm , & + "If True, use the ALE algorithm (regridding/remapping).\n"//& + "If False, use the layered isopycnal algorithm.", default=.false. ) call get_param(param_file, "MOM", "BULKMIXEDLAYER", CS%bulkmixedlayer, & "If true, use a Kraus-Turner-like bulk mixed layer \n"//& "with transitional buffer layers. Layers 1 through \n"//& "NKML+NKBL have variable densities. There must be at \n"//& "least NKML+NKBL+1 layers if BULKMIXEDLAYER is true. \n"//& - "The default is the same setting as ENABLE_THERMODYNAMICS.", & - default=CS%use_temperature) - call get_param(param_file, "MOM", "USE_REGRIDDING", & - CS%use_ALE_algorithm , & - "If True, use the ALE algorithm (regridding/remapping).\n"//& - "If False, use the layered isopycnal algorithm.", default=.false. ) - if (CS%use_ALE_algorithm) then - if (CS%bulkmixedlayer) call MOM_error(FATAL, & - "MOM: BULKMIXEDLAYER can not currently be used with the ALE algorithm.") - if (.not. CS%use_temperature) call MOM_error(FATAL, & - "MOM: At this time, USE_EOS should be True when using the ALE algorithm.") - endif - call get_param(param_file, "MOM", "DO_DYNAMICS", CS%do_dynamics, & - "If False, skips the dynamics calls that update u & v, as well as\n"//& - "the gravity wave adjustment to h. This is a fragile feature and\n"//& - "thus undocumented.", default=.true., do_not_log=.true. ) + "BULKMIXEDLAYER can not be used with USE_REGRIDDING. \n"//& + "The default is influenced by ENABLE_THERMODYNAMICS.", & + default=CS%use_temperature .and. .not.CS%use_ALE_algorithm) call get_param(param_file, "MOM", "THICKNESSDIFFUSE", CS%thickness_diffuse, & "If true, interface heights are diffused with a \n"//& "coefficient of KHTH.", default=.false.) @@ -1625,7 +1621,11 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) write_geom_files = ((write_geom==2) .or. ((write_geom==1) .and. & ((dirs%input_filename(1:1)=='n') .and. (LEN_TRIM(dirs%input_filename)==1)))) - ! Check for inconsistent settings. + ! Check for inconsistent parameter settings. + if (CS%use_ALE_algorithm .and. CS%bulkmixedlayer) call MOM_error(FATAL, & + "MOM: BULKMIXEDLAYER can not currently be used with the ALE algorithm.") + if (CS%use_ALE_algorithm .and. .not.CS%use_temperature) call MOM_error(FATAL, & + "MOM: At this time, USE_EOS should be True when using the ALE algorithm.") if (CS%adiabatic .and. CS%use_temperature) call MOM_error(WARNING, & "MOM: ADIABATIC and ENABLE_THERMODYNAMICS both defined is usually unwise.") if (use_EOS .and. .not.CS%use_temperature) call MOM_error(FATAL, & @@ -1633,6 +1633,8 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) if (CS%adiabatic .and. CS%bulkmixedlayer) call MOM_error(FATAL, & "MOM: ADIABATIC and BULKMIXEDLAYER can not both be defined.") + call callTree_waypoint("MOM parameters read (initialize_MOM)") + ! Set up the model domain and grids. #ifdef SYMMETRIC_MEMORY_ symmetric = .true. @@ -1662,13 +1664,16 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) call verticalGridInit( param_file, CS%GV ) GV => CS%GV - dG%g_Earth = GV%g_Earth + dG%g_Earth = GV%g_Earth ; G%g_Earth = GV%g_Earth ! Allocate the auxiliary non-symmetric domain for debugging or I/O purposes. if (CS%debug .or. dG%symmetric) & call clone_MOM_domain(dG%Domain, dG%Domain_aux, symmetric=.false.) + call callTree_waypoint("grids initialized (initialize_MOM)") + + call MOM_timing_init(CS) call tracer_registry_init(param_file, CS%tracer_Reg) @@ -2052,9 +2057,9 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) if (.not.query_initialized(CS%ave_ssh,"ave_ssh",CS%restart_CSp)) then if (CS%split) then - call find_eta(CS%h, CS%tv, G%g_Earth, G, GV, CS%ave_ssh, eta) + call find_eta(CS%h, CS%tv, GV%g_Earth, G, GV, CS%ave_ssh, eta) else - call find_eta(CS%h, CS%tv, G%g_Earth, G, GV, CS%ave_ssh) + call find_eta(CS%h, CS%tv, GV%g_Earth, G, GV, CS%ave_ssh) endif endif if (CS%split) deallocate(eta) @@ -2099,7 +2104,7 @@ subroutine finish_MOM_initialization(Time, dirs, CS, fluxes) allocate(restart_CSp_tmp) restart_CSp_tmp = CS%restart_CSp allocate(z_interface(SZI_(G),SZJ_(G),SZK_(G)+1)) - call find_eta(CS%h, CS%tv, G%g_Earth, G, GV, z_interface) + call find_eta(CS%h, CS%tv, GV%g_Earth, G, GV, z_interface) vd = var_desc("eta","meter","Interface heights",z_grid='i') call register_restart_field(z_interface, vd, .true., restart_CSp_tmp) @@ -2844,7 +2849,7 @@ subroutine calculate_surface_state(state, u, v, h, ssh, G, GV, CS, p_atm) state%sea_lev => ssh if (present(p_atm)) then ; if (ASSOCIATED(p_atm)) then - IgR0 = 1.0 / (GV%Rho0 * G%g_Earth) + IgR0 = 1.0 / (GV%Rho0 * GV%g_Earth) do j=js,je ; do i=is,ie ssh(i,j) = ssh(i,j) + p_atm(i,j) * IgR0 enddo ; enddo diff --git a/src/core/MOM_PressureForce_Montgomery.F90 b/src/core/MOM_PressureForce_Montgomery.F90 index eea083f9ab..c0418aeb0a 100644 --- a/src/core/MOM_PressureForce_Montgomery.F90 +++ b/src/core/MOM_PressureForce_Montgomery.F90 @@ -195,7 +195,7 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, CS, p_atm, pbce, "can no longer be used with a compressible EOS. Use #define ANALYTIC_FV_PGF.") endif - I_gEarth = 1.0 / G%g_Earth + I_gEarth = 1.0 / GV%g_Earth dp_neglect = GV%H_to_Pa * GV%H_subroundoff !$OMP parallel default(none) shared(nz,alpha_Lay,GV,dalpha_int) !$OMP do @@ -263,12 +263,12 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, CS, p_atm, pbce, call calc_tidal_forcing(CS%Time, SSH, e_tidal, G, CS%tides_CSp) !$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,geopot_bot,G,e_tidal) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - geopot_bot(i,j) = -G%g_Earth*(e_tidal(i,j) + G%bathyT(i,j)) + geopot_bot(i,j) = -GV%g_Earth*(e_tidal(i,j) + G%bathyT(i,j)) enddo ; enddo else !$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,geopot_bot,G) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - geopot_bot(i,j) = -G%g_Earth*G%bathyT(i,j) + geopot_bot(i,j) = -GV%g_Earth*G%bathyT(i,j) enddo ; enddo endif @@ -366,7 +366,7 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, CS, p_atm, pbce, ! Note that ddM/dPb = alpha_star(i,j,1) if (present(pbce)) then - call Set_pbce_nonBouss(p, tv_tmp, G, GV, G%g_Earth, CS%GFS_scale, pbce, & + call Set_pbce_nonBouss(p, tv_tmp, G, GV, GV%g_Earth, CS%GFS_scale, pbce, & alpha_star) endif @@ -510,7 +510,7 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, CS, p_atm, pbce, eta h_neglect = GV%H_subroundoff * GV%H_to_m I_Rho0 = 1.0/CS%Rho0 - G_Rho0 = G%g_Earth/GV%Rho0 + G_Rho0 = GV%g_Earth/GV%Rho0 if (CS%tides) then ! Determine the surface height anomaly for calculating self attraction @@ -619,7 +619,7 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, CS, p_atm, pbce, eta endif ! use_EOS if (present(pbce)) then - call Set_pbce_Bouss(e, tv_tmp, G, GV, G%g_Earth, CS%Rho0, CS%GFS_scale, pbce, & + call Set_pbce_Bouss(e, tv_tmp, G, GV, GV%g_Earth, CS%Rho0, CS%GFS_scale, pbce, & rho_star) endif @@ -1000,7 +1000,7 @@ subroutine PressureForce_Mont_init(Time, G, GV, param_file, diag, CS, tides_CSp) endif CS%GFS_scale = 1.0 - if (GV%g_prime(1) /= G%g_Earth) CS%GFS_scale = GV%g_prime(1) / G%g_Earth + if (GV%g_prime(1) /= GV%g_Earth) CS%GFS_scale = GV%g_prime(1) / GV%g_Earth call log_param(param_file, mod, "GFS / G_EARTH", CS%GFS_scale) diff --git a/src/core/MOM_PressureForce_analytic_FV.F90 b/src/core/MOM_PressureForce_analytic_FV.F90 index 2d92fa6a3f..21aabfabbc 100644 --- a/src/core/MOM_PressureForce_analytic_FV.F90 +++ b/src/core/MOM_PressureForce_analytic_FV.F90 @@ -252,7 +252,7 @@ subroutine PressureForce_AFV_nonBouss(h, tv, PFu, PFv, G, GV, CS, p_atm, pbce, e enddo ; enddo ; enddo !$OMP end parallel - I_gEarth = 1.0 / G%g_Earth + I_gEarth = 1.0 / GV%g_Earth if (use_EOS) then ! With a bulk mixed layer, replace the T & S of any layers that are @@ -323,7 +323,7 @@ subroutine PressureForce_AFV_nonBouss(h, tv, PFu, PFv, G, GV, CS, p_atm, pbce, e !$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,za,alpha_ref,p,G,dza) do j=Jsq,Jeq+1 do i=Isq,Ieq+1 - za(i,j) = alpha_ref*p(i,j,nz+1) - G%g_Earth*G%bathyT(i,j) + za(i,j) = alpha_ref*p(i,j,nz+1) - GV%g_Earth*G%bathyT(i,j) enddo do k=nz,1,-1 ; do i=Isq,Ieq+1 za(i,j) = za(i,j) + dza(i,j,k) @@ -339,7 +339,7 @@ subroutine PressureForce_AFV_nonBouss(h, tv, PFu, PFv, G, GV, CS, p_atm, pbce, e call calc_tidal_forcing(CS%Time, SSH, e_tidal, G, CS%tides_CSp) !$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,za,G,e_tidal) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - za(i,j) = za(i,j) - G%g_Earth*e_tidal(i,j) + za(i,j) = za(i,j) - GV%g_Earth*e_tidal(i,j) enddo ; enddo endif @@ -439,7 +439,7 @@ subroutine PressureForce_AFV_nonBouss(h, tv, PFu, PFv, G, GV, CS, p_atm, pbce, e enddo if (present(pbce)) then - call set_pbce_nonBouss(p, tv_tmp, G, GV, G%g_Earth, CS%GFS_scale, pbce) + call set_pbce_nonBouss(p, tv_tmp, G, GV, GV%g_Earth, CS%GFS_scale, pbce) endif if (present(eta)) then @@ -574,7 +574,7 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, GV, CS, ALE_CSp, p_atm, p PRScheme = pressureReconstructionScheme(ALE_CSp) h_neglect = GV%H_subroundoff I_Rho0 = 1.0/GV%Rho0 - G_Rho0 = G%g_Earth/GV%Rho0 + G_Rho0 = GV%g_Earth/GV%Rho0 rho_ref = CS%Rho0 if (CS%tides) then @@ -711,12 +711,12 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, GV, CS, ALE_CSp, p_atm, p if (use_p_atm) then do jb=Jsq_bk,Jeq_bk+1 ; do ib=Isq_bk,Ieq_bk+1 i = ib+ioff_bk ; j = jb+joff_bk - pa_bk(ib,jb) = (rho_ref*G%g_Earth)*e(i,j,1) + p_atm(i,j) + pa_bk(ib,jb) = (rho_ref*GV%g_Earth)*e(i,j,1) + p_atm(i,j) enddo ; enddo else do jb=Jsq_bk,Jeq_bk+1 ; do ib=Isq_bk,Ieq_bk+1 i = ib+ioff_bk ; j = jb+joff_bk - pa_bk(ib,jb) = (rho_ref*G%g_Earth)*e(i,j,1) + pa_bk(ib,jb) = (rho_ref*GV%g_Earth)*e(i,j,1) enddo ; enddo endif do jb=js_bk,je_bk ; do Ib=Isq_bk,Ieq_bk @@ -740,28 +740,28 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, GV, CS, ALE_CSp, p_atm, p if ( PRScheme == PRESSURE_RECONSTRUCTION_PLM ) then call int_density_dz_generic_plm ( T_t(:,:,k), T_b(:,:,k), & S_t(:,:,k), S_b(:,:,k), e(:,:,K), e(:,:,K+1), & - rho_ref, CS%Rho0, G%g_Earth, & + rho_ref, CS%Rho0, GV%g_Earth, & GV%H_subroundoff, G%bathyT, G%HI, G%Block(n), & tv%eqn_of_state, dpa_bk, intz_dpa_bk, intx_dpa_bk, inty_dpa_bk, & useMassWghtInterp = CS%useMassWghtInterp) elseif ( PRScheme == PRESSURE_RECONSTRUCTION_PPM ) then call int_density_dz_generic_ppm ( tv%T(:,:,k), T_t(:,:,k), T_b(:,:,k), & tv%S(:,:,k), S_t(:,:,k), S_b(:,:,k), e(:,:,K), e(:,:,K+1), & - rho_ref, CS%Rho0, G%G_Earth, & + rho_ref, CS%Rho0, GV%g_Earth, & G%HI, G%Block(n), tv%eqn_of_state, dpa_bk, intz_dpa_bk, & intx_dpa_bk, inty_dpa_bk) endif else call int_density_dz(tv_tmp%T(:,:,k), tv_tmp%S(:,:,k), & e(:,:,K), e(:,:,K+1), & - rho_ref, CS%Rho0, G%g_Earth, G%HI, G%Block(n), tv%eqn_of_state, & + rho_ref, CS%Rho0, GV%g_Earth, G%HI, G%Block(n), tv%eqn_of_state, & dpa_bk, intz_dpa_bk, intx_dpa_bk, inty_dpa_bk ) endif intz_dpa_bk(:,:) = intz_dpa_bk(:,:)*GV%m_to_H else do jb=Jsq_bk,Jeq_bk+1 ; do ib=Isq_bk,Ieq_bk+1 i = ib+ioff_bk ; j = jb+joff_bk - dz_bk(ib,jb) = G%g_Earth*GV%H_to_m*h(i,j,k) + dz_bk(ib,jb) = GV%g_Earth*GV%H_to_m*h(i,j,k) dpa_bk(ib,jb) = (GV%Rlay(k) - rho_ref)*dz_bk(ib,jb) intz_dpa_bk(ib,jb) = 0.5*(GV%Rlay(k) - rho_ref)*dz_bk(ib,jb)*h(i,j,k) enddo ; enddo @@ -813,7 +813,7 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, GV, CS, ALE_CSp, p_atm, p enddo if (present(pbce)) then - call set_pbce_Bouss(e, tv_tmp, G, GV, G%g_Earth, CS%Rho0, CS%GFS_scale, pbce) + call set_pbce_Bouss(e, tv_tmp, G, GV, GV%g_Earth, CS%Rho0, CS%GFS_scale, pbce) endif if (present(eta)) then @@ -890,7 +890,7 @@ subroutine PressureForce_AFV_init(Time, G, GV, param_file, diag, CS, tides_CSp) endif CS%GFS_scale = 1.0 - if (GV%g_prime(1) /= G%g_Earth) CS%GFS_scale = GV%g_prime(1) / G%g_Earth + if (GV%g_prime(1) /= GV%g_Earth) CS%GFS_scale = GV%g_prime(1) / GV%g_Earth call log_param(param_file, mod, "GFS / G_EARTH", CS%GFS_scale) diff --git a/src/core/MOM_dynamics_legacy_split.F90 b/src/core/MOM_dynamics_legacy_split.F90 index 3ad5c3363c..e35f197db5 100644 --- a/src/core/MOM_dynamics_legacy_split.F90 +++ b/src/core/MOM_dynamics_legacy_split.F90 @@ -476,7 +476,7 @@ subroutine step_MOM_dyn_legacy_split(u, v, h, tv, visc, & CS%ALE_CSp, p_surf, CS%pbce, CS%eta_PF) if (dyn_p_surf) then if (GV%Boussinesq) then - Pa_to_eta = 1.0 / (GV%Rho0*G%g_Earth) + Pa_to_eta = 1.0 / (GV%Rho0*GV%g_Earth) else Pa_to_eta = 1.0 / GV%H_to_Pa endif diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index e2f9998b6a..f657381c51 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -696,7 +696,7 @@ subroutine calculateBuoyancyFlux1d(G, GV, fluxes, optics, h, Temp, Salt, tv, j, depthBeforeScalingFluxes = max( GV%Angstrom, 1.e-30*GV%m_to_H ) pressure(:) = 0. ! Ignore atmospheric pressure - GoRho = G%g_Earth / GV%Rho0 + GoRho = GV%g_Earth / GV%Rho0 start = 1 + G%isc - G%isd npts = 1 + G%iec - G%isc @@ -977,13 +977,13 @@ subroutine register_forcing_type_diags(Time, diag, use_temperature, handles) cmor_long_name='Water Evaporation Flux Where Ice Free Ocean over Sea') ! smg: seaice_melt field requires updates to the sea ice model - handles%id_seaice_melt = register_diag_field('ocean_model', 'seaice_melt', & - diag%axesT1, Time, 'water flux to ocean from sea ice melt(> 0) or form(< 0)', & - 'kilogram/(meter^2 * second)', & - standard_name='water_flux_into_sea_water_due_to_sea_ice_thermodynamics', & - cmor_field_name='fsitherm', cmor_units='kg m-2 s-1', & - cmor_standard_name='water_flux_into_sea_water_due_to_sea_ice_thermodynamics',& - cmor_long_name='water flux to ocean from sea ice melt(> 0) or form(< 0)') + !handles%id_seaice_melt = register_diag_field('ocean_model', 'seaice_melt', & + ! diag%axesT1, Time, 'water flux to ocean from sea ice melt(> 0) or form(< 0)', & + ! 'kilogram/(meter^2 * second)', & + ! standard_name='water_flux_into_sea_water_due_to_sea_ice_thermodynamics', & + ! cmor_field_name='fsitherm', cmor_units='kg m-2 s-1', & + ! cmor_standard_name='water_flux_into_sea_water_due_to_sea_ice_thermodynamics',& + ! cmor_long_name='water flux to ocean from sea ice melt(> 0) or form(< 0)') handles%id_precip = register_diag_field('ocean_model', 'precip', diag%axesT1, Time, & 'Liquid + frozen precipitation into ocean', 'kilogram/(meter^2 * second)') @@ -1039,12 +1039,13 @@ subroutine register_forcing_type_diags(Time, diag, use_temperature, handles) cmor_standard_name='water_evaporation_flux_area_integrated', & cmor_long_name='Evaporation Where Ice Free Ocean over Sea Area Integrated') - handles%id_total_seaice_melt = register_scalar_field('ocean_model', 'total_seaice_melt', Time, diag, & - long_name='Area integrated sea ice melt (>0) or form (<0)', units='kg/s', & - standard_name='water_flux_into_sea_water_due_to_sea_ice_thermodynamics_area_integrated', & - cmor_field_name='total_fsitherm', cmor_units='kg s-1', & - cmor_standard_name='water_flux_into_sea_water_due_to_sea_ice_thermodynamics_area_integrated', & - cmor_long_name='Water Melt/Form from Sea Ice Area Integrated') + ! seaice_melt field requires updates to the sea ice model + !handles%id_total_seaice_melt = register_scalar_field('ocean_model', 'total_seaice_melt', Time, diag, & + ! long_name='Area integrated sea ice melt (>0) or form (<0)', units='kg/s', & + ! standard_name='water_flux_into_sea_water_due_to_sea_ice_thermodynamics_area_integrated', & + ! cmor_field_name='total_fsitherm', cmor_units='kg s-1', & + ! cmor_standard_name='water_flux_into_sea_water_due_to_sea_ice_thermodynamics_area_integrated', & + ! cmor_long_name='Water Melt/Form from Sea Ice Area Integrated') handles%id_total_precip = register_scalar_field('ocean_model', 'total_precip', Time, diag, & long_name='Area integrated liquid+frozen precip into ocean', units='kg/s') diff --git a/src/core/MOM_grid.F90 b/src/core/MOM_grid.F90 index fd8f252eb4..de44bd8f7f 100644 --- a/src/core/MOM_grid.F90 +++ b/src/core/MOM_grid.F90 @@ -151,9 +151,10 @@ module MOM_grid !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! !> MOM_grid_init initializes the ocean grid array sizes and grid memory. -subroutine MOM_grid_init(G, param_file) +subroutine MOM_grid_init(G, param_file, HI) type(ocean_grid_type), intent(inout) :: G !< The horizontal grid type type(param_file_type), intent(in) :: param_file !< Parameter file handle + type(hor_index_type), optional, intent(in) :: HI !< A hor_index_type for array extents ! Arguments: G - The ocean's grid structure. ! (in) param_file - A structure indicating the open file to parse for ! model parameter values. @@ -173,9 +174,9 @@ subroutine MOM_grid_init(G, param_file) ! Read all relevant parameters and write them to the model log. call log_version(param_file, mod_nm, version, & "Parameters providing information about the lateral grid.") - call get_param(param_file, "MOM", "G_EARTH", G%g_Earth, & - "The gravitational acceleration of the Earth.", & - units="m s-2", default = 9.80) +! call get_param(param_file, "MOM", "G_EARTH", G%g_Earth, & +! "The gravitational acceleration of the Earth.", & +! units="m s-2", default = 9.80) call get_param(param_file, mod_nm, "GLOBAL_INDEXING", global_indexing, & "If true, use a global lateral indexing convention, so \n"//& "that corresponding points on different processors have \n"//& @@ -204,17 +205,33 @@ subroutine MOM_grid_init(G, param_file) "in the y-direction on each processor (for openmp).", default=1, & layoutParam=.true.) - call hor_index_init(G%Domain, G%HI, param_file, & - local_indexing=.not.global_indexing) - - ! get_domain_extent ensures that domains start at 1 for compatibility between - ! static and dynamically allocated arrays, unless global_indexing is true. - call get_domain_extent(G%Domain, G%isc, G%iec, G%jsc, G%jec, & - G%isd, G%ied, G%jsd, G%jed, & - G%isg, G%ieg, G%jsg, G%jeg, & - G%idg_offset, G%jdg_offset, G%symmetric, & - local_indexing=.not.global_indexing) - G%isd_global = G%isd+G%idg_offset ; G%jsd_global = G%jsd+G%jdg_offset + if (present(HI)) then + G%HI = HI + + G%isc = HI%isc ; G%iec = HI%iec ; G%jsc = HI%jsc ; G%jec = HI%jec + G%isd = HI%isd ; G%ied = HI%ied ; G%jsd = HI%jsd ; G%jed = HI%jed + G%isg = HI%isg ; G%ieg = HI%ieg ; G%jsg = HI%jsg ; G%jeg = HI%jeg + + G%IscB = HI%IscB ; G%IecB = HI%IecB ; G%JscB = HI%JscB ; G%JecB = HI%JecB + G%IsdB = HI%IsdB ; G%IedB = HI%IedB ; G%JsdB = HI%JsdB ; G%JedB = HI%JedB + G%IsgB = HI%IsgB ; G%IegB = HI%IegB ; G%JsgB = HI%JsgB ; G%JegB = HI%JegB + + G%idg_offset = HI%idg_offset ; G%jdg_offset = HI%jdg_offset + G%isd_global = G%isd + HI%idg_offset ; G%jsd_global = G%jsd + HI%jdg_offset + G%symmetric = HI%symmetric + else + call hor_index_init(G%Domain, G%HI, param_file, & + local_indexing=.not.global_indexing) + + ! get_domain_extent ensures that domains start at 1 for compatibility between + ! static and dynamically allocated arrays, unless global_indexing is true. + call get_domain_extent(G%Domain, G%isc, G%iec, G%jsc, G%jec, & + G%isd, G%ied, G%jsd, G%jed, & + G%isg, G%ieg, G%jsg, G%jeg, & + G%idg_offset, G%jdg_offset, G%symmetric, & + local_indexing=.not.global_indexing) + G%isd_global = G%isd+G%idg_offset ; G%jsd_global = G%jsd+G%jdg_offset + endif G%nonblocking_updates = G%Domain%nonblocking_updates diff --git a/src/core/MOM_isopycnal_slopes.F90 b/src/core/MOM_isopycnal_slopes.F90 index c5432934b6..6e12fa1bd4 100644 --- a/src/core/MOM_isopycnal_slopes.F90 +++ b/src/core/MOM_isopycnal_slopes.F90 @@ -94,7 +94,7 @@ subroutine calc_isoneutral_slopes(G, GV, h, e, tv, dt_kappa_smooth, & present_N2_u = PRESENT(N2_u) present_N2_v = PRESENT(N2_v) - G_Rho0 = G%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. diff --git a/src/diagnostics/MOM_diag_to_Z.F90 b/src/diagnostics/MOM_diag_to_Z.F90 index a1ddc3adb2..4c019bd9d8 100644 --- a/src/diagnostics/MOM_diag_to_Z.F90 +++ b/src/diagnostics/MOM_diag_to_Z.F90 @@ -153,7 +153,11 @@ function global_z_mean(var,G,CS,tracer) global_weight_scalar = reproducing_sum(weight,sums=weightij) do k=1, nz - global_z_mean(k) = scalarij(k) / weightij(k) + if (scalarij(k) == 0) then + global_z_mean(k) = 0.0 + else + global_z_mean(k) = scalarij(k) / weightij(k) + endif enddo end function global_z_mean diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 86e6812969..9667ba03cf 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -243,7 +243,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, fluxes, & call calculate_derivs(dt, G, CS) if (ASSOCIATED(CS%e)) then - call find_eta(h, tv, G%g_Earth, G, GV, CS%e, eta_bt) + call find_eta(h, tv, GV%g_Earth, G, GV, CS%e, eta_bt) if (CS%id_e > 0) call post_data(CS%id_e, CS%e, CS%diag) endif @@ -253,7 +253,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, fluxes, & CS%e_D(i,j,k) = CS%e(i,j,k) + G%bathyT(i,j) enddo ; enddo ; enddo else - call find_eta(h, tv, G%g_Earth, G, GV, CS%e_D, eta_bt) + call find_eta(h, tv, GV%g_Earth, G, GV, CS%e_D, eta_bt) do k=1,nz+1 ; do j=js,je ; do i=is,ie CS%e_D(i,j,k) = CS%e_D(i,j,k) + G%bathyT(i,j) enddo ; enddo ; enddo @@ -304,7 +304,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, fluxes, & do k=1,nz ! pressure for EOS at the layer center (Pa) do i=is,ie - pressure_1d(i) = pressure_1d(i) + 0.5*(G%G_Earth*GV%H_to_kg_m2)*h(i,j,k) + pressure_1d(i) = pressure_1d(i) + 0.5*(GV%g_Earth*GV%H_to_kg_m2)*h(i,j,k) enddo ! store in-situ density (kg/m3) in diag_tmp3d call calculate_density(tv%T(:,j,k),tv%S(:,j,k), pressure_1d, & @@ -315,7 +315,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, fluxes, & enddo ! pressure for EOS at the bottom interface (Pa) do i=is,ie - pressure_1d(i) = pressure_1d(i) + 0.5*(G%G_Earth*GV%H_to_kg_m2)*h(i,j,k) + pressure_1d(i) = pressure_1d(i) + 0.5*(GV%g_Earth*GV%H_to_kg_m2)*h(i,j,k) enddo enddo ! k @@ -707,7 +707,7 @@ subroutine calculate_vertical_integrals(h, tv, fluxes, G, GV, CS) endif if (CS%id_col_ht > 0) then - call find_eta(h, tv, G%g_Earth, G, GV, z_top) + call find_eta(h, tv, GV%g_Earth, G, GV, z_top) do j=js,je ; do i=is,ie z_bot(i,j) = z_top(i,j) + G%bathyT(i,j) enddo ; enddo @@ -718,7 +718,7 @@ subroutine calculate_vertical_integrals(h, tv, fluxes, G, GV, CS) do j=js,je ; do i=is,ie ; mass(i,j) = 0.0 ; enddo ; enddo if (GV%Boussinesq) then if (associated(tv%eqn_of_state)) then - IG_Earth = 1.0 / G%g_Earth + IG_Earth = 1.0 / GV%g_Earth ! do j=js,je ; do i=is,ie ; z_bot(i,j) = -P_SURF(i,j)/GV%H_to_Pa ; enddo ; enddo do j=js,je ; do i=is,ie ; z_bot(i,j) = 0.0 ; enddo ; enddo do k=1,nz @@ -727,7 +727,7 @@ subroutine calculate_vertical_integrals(h, tv, fluxes, G, GV, CS) z_bot(i,j) = z_top(i,j) - GV%H_to_m*h(i,j,k) enddo ; enddo call int_density_dz(tv%T(:,:,k), tv%S(:,:,k), & - z_top, z_bot, 0.0, GV%H_to_kg_m2, G%g_Earth, & + z_top, z_bot, 0.0, GV%H_to_kg_m2, GV%g_Earth, & G%HI, G%HI, tv%eqn_of_state, dpress) do j=js,je ; do i=is,ie mass(i,j) = mass(i,j) + dpress(i,j) * IG_Earth @@ -753,7 +753,7 @@ subroutine calculate_vertical_integrals(h, tv, fluxes, G, GV, CS) ! where pso is the sea water pressure at sea water surface ! note that pso is equivalent to fluxes%p_surf do j=js,je ; do i=is,ie - btm_pres(i,j) = mass(i,j) * G%g_Earth + btm_pres(i,j) = mass(i,j) * GV%g_Earth if (ASSOCIATED(fluxes%p_surf)) then btm_pres(i,j) = btm_pres(i,j) + fluxes%p_surf(i,j) endif diff --git a/src/diagnostics/MOM_obsolete_diagnostics.F90 b/src/diagnostics/MOM_obsolete_diagnostics.F90 index 9100367ab5..4b498292d7 100644 --- a/src/diagnostics/MOM_obsolete_diagnostics.F90 +++ b/src/diagnostics/MOM_obsolete_diagnostics.F90 @@ -6,7 +6,8 @@ module MOM_obsolete_diagnostics use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe use MOM_file_parser, only : param_file_type, log_version, get_param -use MOM_diag_mediator, only : diag_ctrl, register_static_field +use MOM_diag_mediator, only : diag_ctrl +use diag_manager_mod, only : register_static_field_fms=>register_static_field implicit none ; private @@ -61,7 +62,7 @@ subroutine register_obsolete_diagnostics(param_file, diag) end subroutine register_obsolete_diagnostics !> Fakes a register of a diagnostic to find out if an obsolete -!! parameter appears in the diag_table. +!! parameter appears in the diag_table. logical function found_in_diagtable(diag, varName, newVarName) type(diag_ctrl), intent(in) :: diag character(len=*), intent(in) :: varName @@ -69,14 +70,11 @@ logical function found_in_diagtable(diag, varName, newVarName) ! Local integer :: handle ! Integer handle returned from diag_manager - ! Note that using register_static_field() means the diagnostic - ! does not appear in the available diagnostics list. If we change - ! register_static_field to record available diagnostics also then - ! we will have to call the FMS version of register_field instead. - ! IOW, this is a leveraging a "feature" that might go away. -AJA - handle = register_static_field('ocean_model', varName, diag%axesT1, & - 'Obsolete parameter', 'N/A') - + ! We use register_static_field_fms() instead of register_static_field() so + ! that the diagnostic does not appear in the available diagnostics list. + handle = register_static_field_fms('ocean_model', varName, & + diag%axesT1%handles, 'Obsolete parameter', 'N/A') + found_in_diagtable = (handle>0) if (handle>0 .and. is_root_pe()) then diff --git a/src/diagnostics/MOM_sum_output.F90 b/src/diagnostics/MOM_sum_output.F90 index d0aa3cd65c..da6c0f5b63 100644 --- a/src/diagnostics/MOM_sum_output.F90 +++ b/src/diagnostics/MOM_sum_output.F90 @@ -474,7 +474,7 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, CS, tracer_CSp) enddo ; enddo ; enddo mass_tot = reproducing_sum(tmp1, sums=mass_lay, EFP_sum=mass_EFP) - call find_eta(h, tv, G%g_Earth, G, GV, eta) + call find_eta(h, tv, GV%g_Earth, G, GV, eta) do k=1,nz ; do j=js,je ; do i=is,ie tmp1(i,j,k) = (eta(i,j,K)-eta(i,j,K+1)) * areaTm(i,j) enddo ; enddo ; enddo @@ -494,7 +494,7 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, CS, tracer_CSp) enddo ; enddo enddo if (CS%do_APE_calc) then - call find_eta(h, tv, G%g_Earth, G, GV, eta) + call find_eta(h, tv, GV%g_Earth, G, GV, eta) do k=1,nz vol_lay(k) = 0.0 diff --git a/src/diagnostics/MOM_wave_speed.F90 b/src/diagnostics/MOM_wave_speed.F90 index 9079124e6a..b84ec40432 100644 --- a/src/diagnostics/MOM_wave_speed.F90 +++ b/src/diagnostics/MOM_wave_speed.F90 @@ -160,10 +160,10 @@ subroutine wave_speed(h, tv, G, GV, cg1, CS, full_halos) endif ; endif S => tv%S ; T => tv%T - g_Rho0 = G%g_Earth/GV%Rho0 + g_Rho0 = GV%g_Earth/GV%Rho0 use_EOS = associated(tv%eqn_of_state) - H_to_pres = G%g_Earth * GV%Rho0 + H_to_pres = GV%g_Earth * GV%Rho0 H_to_m = GV%H_to_m rescale = 1024.0**4 ; I_rescale = 1.0/rescale @@ -515,10 +515,10 @@ subroutine wave_speeds(h, tv, G, GV, nmodes, cn, CS, full_halos) endif ; endif S => tv%S ; T => tv%T - g_Rho0 = G%g_Earth/GV%Rho0 + g_Rho0 = GV%g_Earth/GV%Rho0 use_EOS = associated(tv%eqn_of_state) - H_to_pres = G%g_Earth * GV%Rho0 + H_to_pres = GV%g_Earth * GV%Rho0 H_to_m = GV%H_to_m min_h_frac = tol1 / real(nz) diff --git a/src/diagnostics/MOM_wave_structure.F90 b/src/diagnostics/MOM_wave_structure.F90 index f0767255ed..7c47d88971 100644 --- a/src/diagnostics/MOM_wave_structure.F90 +++ b/src/diagnostics/MOM_wave_structure.F90 @@ -213,10 +213,10 @@ subroutine wave_structure(h, tv, G, GV, cn, ModeNum, freq, CS, En, full_halos) Pi = (4.0*atan(1.0)) S => tv%S ; T => tv%T - g_Rho0 = G%g_Earth/GV%Rho0 + g_Rho0 = GV%g_Earth/GV%Rho0 use_EOS = associated(tv%eqn_of_state) - H_to_pres = G%g_Earth * GV%Rho0 + H_to_pres = GV%g_Earth * GV%Rho0 H_to_m = GV%H_to_m rescale = 1024.0**4 ; I_rescale = 1.0/rescale diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 14cde0c674..b4533b3de5 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -866,7 +866,6 @@ subroutine add_shelf_flux(G, CS, state, fluxes) tauy2 = (asv1 * state%tauy_shelf(i,j-1)**2 + & asv2 * state%tauy_shelf(i,j)**2 ) / (asv1 + asv2) fluxes%ustar(i,j) = MAX(CS%ustar_bg, sqrt(Irho0 * sqrt(taux2 + tauy2))) - if (associated(fluxes%sw)) fluxes%sw(i,j) = 0.0 if (associated(fluxes%lw)) fluxes%lw(i,j) = 0.0 if (associated(fluxes%latent)) fluxes%latent(i,j) = 0.0 @@ -1389,16 +1388,9 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, fluxes, Ti ! Allocate the arrays for passing ice-shelf data through the forcing type. if (.not. solo_ice_sheet) then if (is_root_pe()) print *,"initialize_ice_shelf: allocating fluxes" - allocate( fluxes%frac_shelf_h(isd:ied,jsd:jed) ) ; fluxes%frac_shelf_h(:,:) = 0.0 - allocate( fluxes%frac_shelf_u(Isdq:Iedq,jsd:jed) ) ; fluxes%frac_shelf_u(:,:) = 0.0 - allocate( fluxes%frac_shelf_v(isd:ied,Jsdq:Jedq) ) ; fluxes%frac_shelf_v(:,:) = 0.0 - allocate( fluxes%ustar_shelf(isd:ied,jsd:jed) ) ; fluxes%ustar_shelf(:,:) = 0.0 - allocate( fluxes%iceshelf_melt(isd:ied,jsd:jed) ) ; fluxes%iceshelf_melt(:,:) = 0.0 - if (.not.associated(fluxes%p_surf)) then - allocate( fluxes%p_surf(isd:ied,jsd:jed) ) ; fluxes%p_surf(:,:) = 0.0 - endif - allocate( fluxes%rigidity_ice_u(Isdq:Iedq,jsd:jed) ) ; fluxes%rigidity_ice_u(:,:) = 0.0 - allocate( fluxes%rigidity_ice_v(isd:ied,Jsdq:Jedq) ) ; fluxes%rigidity_ice_v(:,:) = 0.0 + !GM, the following is necessary to make ENERGETICS_SFC_PBL work + call allocate_forcing_type(G, fluxes, ustar=.true., shelf=.true., & + press=.true., water=.true., heat=.true.) else if (is_root_pe()) print *,"allocating fluxes in solo mode" call allocate_forcing_type(G, fluxes, ustar=.true., shelf=.true., press=.true.) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 274b3e0bf8..4f8ff9fbd3 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -713,7 +713,7 @@ subroutine convert_thickness(h, G, GV, param_file, tv) Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB max_itt = 10 Boussinesq = GV%Boussinesq - I_gEarth = 1.0 / G%g_Earth + I_gEarth = 1.0 / GV%g_Earth if (Boussinesq) then call MOM_error(FATAL,"Not yet converting thickness with Boussinesq approx.") @@ -728,7 +728,7 @@ subroutine convert_thickness(h, G, GV, param_file, tv) call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_top(:,j), rho, & is, ie-is+1, tv%eqn_of_state) do i=is,ie - p_bot(i,j) = p_top(i,j) + G%g_Earth * h(i,j,k) * rho(i) + p_bot(i,j) = p_top(i,j) + GV%g_Earth * h(i,j,k) * rho(i) enddo enddo @@ -742,7 +742,7 @@ subroutine convert_thickness(h, G, GV, param_file, tv) ! The hydrostatic equation is linear to such a ! high degree that no bounds-checking is needed. do i=is,ie - p_bot(i,j) = p_bot(i,j) + rho(i) * (G%g_Earth*h(i,j,k) - dz_geo(i,j)) + p_bot(i,j) = p_bot(i,j) + rho(i) * (GV%g_Earth*h(i,j,k) - dz_geo(i,j)) enddo enddo ; endif enddo @@ -808,7 +808,7 @@ subroutine depress_surface(h, G, GV, param_file, tv) enddo ; enddo ; endif ! Convert thicknesses to interface heights. - call find_eta(h, tv, G%g_Earth, G, GV, eta) + call find_eta(h, tv, GV%g_Earth, G, GV, eta) do j=js,je ; do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then ! if (eta_sfc(i,j) < eta(i,j,nz+1)) then @@ -902,7 +902,7 @@ subroutine trim_for_ice(PF, G, GV, ALE_CSp, tv, h) endif do j=G%jsc,G%jec ; do i=G%isc,G%iec - call cut_off_column_top(GV%ke, tv, GV%Rho0, G%G_earth, G%bathyT(i,j), min_thickness, & + call cut_off_column_top(GV%ke, tv, GV%Rho0, GV%g_Earth, G%bathyT(i,j), min_thickness, & tv%T(i,j,:), T_t(i,j,:), T_b(i,j,:), tv%S(i,j,:), S_t(i,j,:), S_b(i,j,:), & p_surf(i,j), h(i,j,:), remap_CS) enddo ; enddo @@ -2530,14 +2530,14 @@ subroutine MOM_state_init_tests(G, GV, tv) S_t(k) = 35.-(0./500.)*e(k) S(k) = 35.+(0./500.)*z(k) S_b(k) = 35.-(0./500.)*e(k+1) - call calculate_density(0.5*(T_t(k)+T_b(k)), 0.5*(S_t(k)+S_b(k)), -GV%Rho0*G%G_earth*z(k), rho(k), tv%eqn_of_state) - P_tot = P_tot + G%G_earth * rho(k) * h(k) + call calculate_density(0.5*(T_t(k)+T_b(k)), 0.5*(S_t(k)+S_b(k)), -GV%Rho0*GV%g_Earth*z(k), rho(k), tv%eqn_of_state) + P_tot = P_tot + GV%g_Earth * rho(k) * h(k) enddo P_t = 0. do k = 1, nk call find_depth_of_pressure_in_cell(T_t(k), T_b(k), S_t(k), S_b(k), e(K), e(K+1), & - P_t, 0.5*P_tot, GV%Rho0, G%G_earth, tv%eqn_of_state, P_b, z_out) + P_t, 0.5*P_tot, GV%Rho0, GV%g_Earth, tv%eqn_of_state, P_b, z_out) write(0,*) k,P_t,P_b,0.5*P_tot,e(K),e(K+1),z_out P_t = P_b enddo @@ -2547,7 +2547,7 @@ subroutine MOM_state_init_tests(G, GV, tv) write(0,*) ' ==================================================================== ' write(0,*) '' write(0,*) h - call cut_off_column_top(nk, tv, GV%Rho0, G%G_earth, -e(nk+1), GV%Angstrom, & + call cut_off_column_top(nk, tv, GV%Rho0, GV%g_Earth, -e(nk+1), GV%Angstrom, & T, T_t, T_b, S, S_t, S_b, 0.5*P_tot, h, remap_CS) write(0,*) h diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index daa27c6b8c..56b972e9e9 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -371,7 +371,7 @@ subroutine calc_slope_functions(h, tv, dt, G, GV, CS) if (.not. ASSOCIATED(CS)) call MOM_error(FATAL, "MOM_lateral_mixing_coeffs.F90, calc_slope_functions:"//& "Module must be initialized before it is used.") - call find_eta(h, tv, G%g_Earth, G, GV, e, halo_size=2) + call find_eta(h, tv, GV%g_Earth, G, GV, e, halo_size=2) if (CS%use_variable_mixing) then if (CS%use_stored_slopes) then call calc_isoneutral_slopes(G, GV, h, e, tv, dt*CS%kappa_smooth, & diff --git a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 index ab939f2cdd..3f35afa4cc 100644 --- a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 +++ b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 @@ -223,7 +223,7 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, fluxes, dt, MLD, G, GV, uDml(:) = 0.0 ; vDml(:) = 0.0 I4dt = 0.25 / dt - g_Rho0 = G%g_Earth/GV%Rho0 + g_Rho0 = GV%g_Earth/GV%Rho0 h_neglect = GV%H_subroundoff dz_neglect = GV%H_subroundoff*GV%H_to_m @@ -465,7 +465,7 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, fluxes, dt, G, GV, CS) uDml(:) = 0.0 ; vDml(:) = 0.0 I4dt = 0.25 / dt - g_Rho0 = G%g_Earth/GV%Rho0 + g_Rho0 = GV%g_Earth/GV%Rho0 use_EOS = associated(tv%eqn_of_state) h_neglect = GV%H_subroundoff dz_neglect = GV%H_subroundoff*GV%H_to_m diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index 03569e6870..1290e6432a 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -202,7 +202,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, MEKE, VarMix, CDp, CS (dt*(G%IdxCv(i,J)*G%IdxCv(i,J) + G%IdyCv(i,J)*G%IdyCv(i,J))) enddo ; enddo - call find_eta(h, tv, G%g_Earth, G, GV, e, halo_size=1) + call find_eta(h, tv, GV%g_Earth, G, GV, e, halo_size=1) ! Set the diffusivities. !$OMP parallel default(none) shared(is,ie,js,je,Khth_Loc_u,CS,use_VarMix,VarMix, & @@ -546,7 +546,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, dt, G, GV, MEK H_to_m = GV%H_to_m ; m_to_H = GV%m_to_H I4dt = 0.25 / dt I_slope_max2 = 1.0 / (CS%slope_max**2) - G_scale = G%g_Earth * H_to_m + G_scale = GV%g_Earth * H_to_m h_neglect = GV%H_subroundoff ; h_neglect2 = h_neglect**2 dz_neglect = GV%H_subroundoff*H_to_m diff --git a/src/parameterizations/vertical/MOM_KPP.F90 b/src/parameterizations/vertical/MOM_KPP.F90 index bebdad862a..7bc9b74187 100644 --- a/src/parameterizations/vertical/MOM_KPP.F90 +++ b/src/parameterizations/vertical/MOM_KPP.F90 @@ -448,7 +448,7 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & #endif ! some constants - GoRho = G%g_Earth / GV%Rho0 + GoRho = GV%g_Earth / GV%Rho0 nonLocalTrans(:,:) = 0.0 if (CS%id_Kd_in > 0) call post_data(CS%id_Kd_in, Kt, CS%diag) diff --git a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 index 815aebec9a..7ba5896c99 100644 --- a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 +++ b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 @@ -563,7 +563,7 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, CS, & ! rivermix_depth = The prescribed depth over which to mix river inflow ! drho_ds = The gradient of density wrt salt at the ambient surface salinity. ! Sriver = 0 (i.e. rivers are assumed to be pure freshwater) - RmixConst = 0.5*CS%rivermix_depth*G%g_Earth*Irho0**2 + RmixConst = 0.5*CS%rivermix_depth*GV%g_Earth*Irho0**2 do i=is,ie TKE_river(i) = max(0.0, RmixConst*dR0_dS(i)* & (fluxes%lrunoff(i,j) + fluxes%frunoff(i,j)) * S(i,1)) @@ -908,7 +908,7 @@ subroutine convective_adjustment(h, u, v, R0, Rcv, T, S, eps, d_eb, & integer :: is, ie, nz, i, k, k1, nzc, nkmb is = G%isc ; ie = G%iec ; nz = G%ke - g_H2_2Rho0 = (G%g_Earth * GV%H_to_m**2) / (2.0 * GV%Rho0) + g_H2_2Rho0 = (GV%g_Earth * GV%H_to_m**2) / (2.0 * GV%Rho0) nzc = nz ; if (present(nz_conv)) nzc = nz_conv nkmb = CS%nkml+CS%nkbl @@ -1067,7 +1067,7 @@ subroutine mixedlayer_convection(h, d_eb, htot, Ttot, Stot, uhtot, vhtot, & Angstrom = GV%Angstrom C1_3 = 1.0/3.0 ; C1_6 = 1.0/6.0 - g_H2_2Rho0 = (G%g_Earth * GV%H_to_m**2) / (2.0 * GV%Rho0) + g_H2_2Rho0 = (GV%g_Earth * GV%H_to_m**2) / (2.0 * GV%Rho0) Idt = 1.0/dt is = G%isc ; ie = G%iec ; nz = G%ke @@ -1582,7 +1582,7 @@ subroutine mechanical_entrainment(h, d_eb, htot, Ttot, Stot, uhtot, vhtot, & C1_3 = 1.0/3.0 ; C1_6 = 1.0/6.0 ; C1_24 = 1.0/24.0 H_to_m = GV%H_to_m - g_H_2Rho0 = (G%g_Earth * H_to_m) / (2.0 * GV%Rho0) + g_H_2Rho0 = (GV%g_Earth * H_to_m) / (2.0 * GV%Rho0) Hmix_min = CS%Hmix_min * GV%m_to_H h_neglect = GV%H_subroundoff is = G%isc ; ie = G%iec ; nz = G%ke @@ -2323,8 +2323,8 @@ subroutine mixedlayer_detrain_2(h, T, S, R0, Rcv, RcvTgt, dt, dt_diag, d_ea, j, kb1 = CS%nkml+1; kb2 = CS%nkml+2 nkmb = CS%nkml+CS%nkbl h_neglect = GV%H_subroundoff - G_2 = 0.5*G%g_Earth - Rho0xG = GV%Rho0 * G%g_Earth + G_2 = 0.5*GV%g_Earth + Rho0xG = GV%Rho0 * GV%g_Earth Idt_H2 = GV%H_to_m**2 / dt_diag I2Rho0 = 0.5 / GV%Rho0 Angstrom = GV%Angstrom @@ -3122,8 +3122,8 @@ subroutine mixedlayer_detrain_1(h, T, S, R0, Rcv, RcvTgt, dt, dt_diag, d_ea, d_e "CS%nkbl must be 1 in mixedlayer_detrain_1.") Idt = 1.0/dt dt_Time = dt/Timescale - g_H2_2Rho0dt = (G%g_Earth * GV%H_to_m**2) / (2.0 * GV%Rho0 * dt_diag) - g_H2_2dt = (G%g_Earth * GV%H_to_m**2) / (2.0 * dt_diag) + g_H2_2Rho0dt = (GV%g_Earth * GV%H_to_m**2) / (2.0 * GV%Rho0 * dt_diag) + g_H2_2dt = (GV%g_Earth * GV%H_to_m**2) / (2.0 * dt_diag) ! Move detrained water into the buffer layer. do k=1,CS%nkml diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index b11cd75421..587e55f6ce 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -746,7 +746,7 @@ subroutine diagnoseMLDbyDensityDifference(id_MLD, h, tv, densityDiff, G, GV, dia subMLN2(:,j) = 0. rho1(:) = 0. d1(:) = 0. - pRef_N2(:) = G%g_Earth * GV%Rho0 * h(:,j,1) * GV%H_to_m ! Boussinesq approximation!!!! ????? + pRef_N2(:) = GV%g_Earth * GV%Rho0 * h(:,j,1) * GV%H_to_m ! Boussinesq approximation!!!! ????? endif do k = 2, nz dKm1(:) = dK(:) ! Depth of center of layer K-1 @@ -754,7 +754,7 @@ subroutine diagnoseMLDbyDensityDifference(id_MLD, h, tv, densityDiff, G, GV, dia ! Stratification, N2, immediately below the mixed layer, averaged over at least 50 m. if (id_N2>0) then - pRef_N2(:) = pRef_N2(:) + G%g_Earth * GV%Rho0 * h(:,j,k) * GV%H_to_m ! Boussinesq approximation!!!! ????? + pRef_N2(:) = pRef_N2(:) + GV%g_Earth * GV%Rho0 * h(:,j,k) * GV%H_to_m ! Boussinesq approximation!!!! ????? call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pRef_N2, rhoAtK, is, ie-is+1, tv%eqn_of_state) do i = is, ie if (MLD(i,j)>0. .and. subMLN2(i,j)==0.) then ! This block is below the mixed layer @@ -762,10 +762,10 @@ subroutine diagnoseMLDbyDensityDifference(id_MLD, h, tv, densityDiff, G, GV, dia rho1(i) = rhoAtK(i) d1(i) = dK(i) ! Use pressure at the bottom of the upper layer used in calculating d/dz rho - pRef_N2(i) = pRef_N2(i) + G%g_Earth * GV%Rho0 * h(i,j,k) * GV%H_to_m ! Boussinesq approximation!!!! ????? + pRef_N2(i) = pRef_N2(i) + GV%g_Earth * GV%Rho0 * h(i,j,k) * GV%H_to_m ! Boussinesq approximation!!!! ????? endif if (d1(i)>0. .and. dK(i)-d1(i)>=dz_subML) then - subMLN2(i,j) = G%g_Earth/ GV%Rho0 * (rho1(i)-rhoAtK(i)) / (d1(i) - dK(i)) + subMLN2(i,j) = GV%g_Earth/ GV%Rho0 * (rho1(i)-rhoAtK(i)) / (d1(i) - dK(i)) endif endif enddo ! i-loop @@ -789,7 +789,7 @@ subroutine diagnoseMLDbyDensityDifference(id_MLD, h, tv, densityDiff, G, GV, dia if ((MLD(i,j)==0.) .and. (deltaRhoAtK(i)0 .and. subMLN2(i,j)==0. .and. d1(i)>0. .and. dK(i)-d1(i)>0.) then ! ! Use what ever stratification we can, measured over what ever distance is available - ! subMLN2(i,j) = G%g_Earth/ GV%Rho0 * (rho1(i)-rhoAtK(i)) / (d1(i) - dK(i)) + ! subMLN2(i,j) = GV%g_Earth/ GV%Rho0 * (rho1(i)-rhoAtK(i)) / (d1(i) - dK(i)) ! endif enddo enddo ! j-loop @@ -862,8 +862,8 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, dt, fluxes, optics, ea, h, tv, & Idt = 1.0/dt calculate_energetics = (present(cTKE) .and. present(dSV_dT) .and. present(dSV_dS)) - I_G_Earth = 1.0 / G%G_earth - g_Hconv2 = G%G_earth * GV%H_to_kg_m2**2 + I_G_Earth = 1.0 / GV%g_Earth + g_Hconv2 = GV%g_Earth * GV%H_to_kg_m2**2 if (present(cTKE)) cTKE(:,:,:) = 0.0 @@ -908,7 +908,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, dt, fluxes, optics, ea, h, tv, & do i=is,ie ; pres(i) = 0.0 ; enddo ! Add surface pressure? do k=1,nz do i=is,ie - d_pres(i) = G%g_Earth * GV%H_to_kg_m2 * h2d(i,k) + d_pres(i) = GV%g_Earth * GV%H_to_kg_m2 * h2d(i,k) p_lay(i) = pres(i) + 0.5*d_pres(i) pres(i) = pres(i) + d_pres(i) enddo @@ -1216,13 +1216,14 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, dt, fluxes, optics, ea, h, tv, & end subroutine applyBoundaryFluxesInOut -subroutine diabatic_aux_init(Time, G, GV, param_file, diag, CS, use_ePBL) +subroutine diabatic_aux_init(Time, G, GV, param_file, diag, CS, useALEalgorithm, use_ePBL) type(time_type), intent(in) :: Time type(ocean_grid_type), intent(in) :: G type(verticalGrid_type), intent(in) :: GV type(param_file_type), intent(in) :: param_file type(diag_ctrl), target, intent(inout) :: diag type(diabatic_aux_CS), pointer :: CS + logical, intent(in) :: useALEalgorithm logical, intent(in) :: use_ePBL ! Arguments: @@ -1306,40 +1307,40 @@ subroutine diabatic_aux_init(Time, G, GV, param_file, diag, CS, use_ePBL) CS%use_calving_heat_content = .false. endif - CS%id_createdH = register_diag_field('ocean_model',"created_H",diag%axesT1, & - Time, "The volume flux added to stop the ocean from drying out and becoming negative in depth", & - "meter second-1") - if (CS%id_createdH>0) allocate(CS%createdH(isd:ied,jsd:jed)) - - - ! diagnostic for heating of a grid cell from convergence of SW heat into the cell - CS%id_penSW_diag = register_diag_field('ocean_model', 'rsdoabsorb', & - diag%axesTL, Time, 'Convergence of Penetrative Shortwave Flux in Sea Water Layer',& - 'Watt meter-2', standard_name='net_rate_of_absorption_of_shortwave_energy_in_ocean_layer') - - ! diagnostic for penetrative SW heat flux at top interface of tracer cell (nz+1 interfaces) - ! k=1 gives penetrative SW at surface; SW(k=nz+1)=0 (no penetration through rock). - CS%id_penSWflux_diag = register_diag_field('ocean_model', 'rsdo', & - diag%axesTi, Time, 'Downwelling Shortwave Flux in Sea Water at Grid Cell Upper Interface',& - 'Watt meter-2', standard_name='downwelling_shortwave_flux_in_sea_water') - - ! need both arrays for the SW diagnostics (one for flux, one for convergence) - if (CS%id_penSW_diag>0 .or. CS%id_penSWflux_diag>0) then - allocate(CS%penSW_diag(isd:ied,jsd:jed,nz)) - CS%penSW_diag(:,:,:) = 0.0 - allocate(CS%penSWflux_diag(isd:ied,jsd:jed,nz+1)) - CS%penSWflux_diag(:,:,:) = 0.0 - endif - + if (useALEalgorithm) then + CS%id_createdH = register_diag_field('ocean_model',"created_H",diag%axesT1, & + Time, "The volume flux added to stop the ocean from drying out and becoming negative in depth", & + "meter second-1") + if (CS%id_createdH>0) allocate(CS%createdH(isd:ied,jsd:jed)) + + ! diagnostic for heating of a grid cell from convergence of SW heat into the cell + CS%id_penSW_diag = register_diag_field('ocean_model', 'rsdoabsorb', & + diag%axesTL, Time, 'Convergence of Penetrative Shortwave Flux in Sea Water Layer',& + 'Watt meter-2', standard_name='net_rate_of_absorption_of_shortwave_energy_in_ocean_layer') + + ! diagnostic for penetrative SW heat flux at top interface of tracer cell (nz+1 interfaces) + ! k=1 gives penetrative SW at surface; SW(k=nz+1)=0 (no penetration through rock). + CS%id_penSWflux_diag = register_diag_field('ocean_model', 'rsdo', & + diag%axesTi, Time, 'Downwelling Shortwave Flux in Sea Water at Grid Cell Upper Interface',& + 'Watt meter-2', standard_name='downwelling_shortwave_flux_in_sea_water') + + ! need both arrays for the SW diagnostics (one for flux, one for convergence) + if (CS%id_penSW_diag>0 .or. CS%id_penSWflux_diag>0) then + allocate(CS%penSW_diag(isd:ied,jsd:jed,nz)) + CS%penSW_diag(:,:,:) = 0.0 + allocate(CS%penSWflux_diag(isd:ied,jsd:jed,nz+1)) + CS%penSWflux_diag(:,:,:) = 0.0 + endif - ! diagnostic for non-downwelling SW radiation (i.e., SW absorbed at ocean surface) - CS%id_nonpenSW_diag = register_diag_field('ocean_model', 'nonpenSW', & - diag%axesT1, Time, & - 'Non-downwelling SW radiation (i.e., SW absorbed in ocean surface with LW,SENS,LAT)',& - 'Watt meter-2', standard_name='nondownwelling_shortwave_flux_in_sea_water') - if (CS%id_nonpenSW_diag > 0) then - allocate(CS%nonpenSW_diag(isd:ied,jsd:jed)) - CS%nonpenSW_diag(:,:) = 0.0 + ! diagnostic for non-downwelling SW radiation (i.e., SW absorbed at ocean surface) + CS%id_nonpenSW_diag = register_diag_field('ocean_model', 'nonpenSW', & + diag%axesT1, Time, & + 'Non-downwelling SW radiation (i.e., SW absorbed in ocean surface with LW,SENS,LAT)',& + 'Watt meter-2', standard_name='nondownwelling_shortwave_flux_in_sea_water') + if (CS%id_nonpenSW_diag > 0) then + allocate(CS%nonpenSW_diag(isd:ied,jsd:jed)) + CS%nonpenSW_diag(:,:) = 0.0 + endif endif id_clock_uv_at_h = cpu_clock_id('(Ocean find_uv_at_h)', grain=CLOCK_ROUTINE) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index ca1a325fb0..57176c2249 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -1944,8 +1944,11 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, !call set_diffusivity_init(Time, G, param_file, diag, CS%set_diff_CSp, diag_to_Z_CSp, CS%int_tide_CSp) CS%id_Kd_interface = register_diag_field('ocean_model', 'Kd_interface', diag%axesTi, Time, & 'Total diapycnal diffusivity at interfaces', 'meter2 second-1') - CS%id_Kd_ePBL = register_diag_field('ocean_model', 'Kd_ePBL', diag%axesTi, Time, & - 'ePBL diapycnal diffusivity at interfaces', 'meter2 second-1') + if (CS%use_energetic_PBL) then + CS%id_Kd_ePBL = register_diag_field('ocean_model', 'Kd_ePBL', diag%axesTi, Time, & + 'ePBL diapycnal diffusivity at interfaces', 'meter2 second-1') + endif + CS%id_Kd_heat = register_diag_field('ocean_model', 'Kd_heat', diag%axesTi, Time, & 'Total diapycnal diffusivity for heat at interfaces', 'meter2 second-1', & cmor_field_name='difvho', cmor_units='m2 s-1', & @@ -1981,120 +1984,118 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, ! diagnostics for tendencies of temp and saln due to diabatic processes; ! available only for ALE algorithm. + if (CS%useALEalgorithm) then + CS%id_diabatic_diff_temp_tend = register_diag_field('ocean_model', & + 'diabatic_diff_temp_tendency', diag%axesTL, Time, & + 'Diabatic diffusion temperature tendency', 'Degree C per second') + if(CS%id_diabatic_diff_temp_tend > 0) then + CS%diabatic_diff_tendency_diag = .true. + endif - CS%id_diabatic_diff_temp_tend = register_diag_field('ocean_model', & - 'diabatic_diff_temp_tendency', diag%axesTL, Time, & - 'Diabatic diffusion temperature tendency', 'Degree C per second') - if(CS%id_diabatic_diff_temp_tend > 0) then - CS%diabatic_diff_tendency_diag = .true. - endif - - CS%id_diabatic_diff_saln_tend = register_diag_field('ocean_model',& - 'diabatic_diff_saln_tendency', diag%axesTL, Time, & - 'Diabatic diffusion salinity tendency', 'PPT per second') - if(CS%id_diabatic_diff_saln_tend > 0) then - CS%diabatic_diff_tendency_diag = .true. - endif - - CS%id_diabatic_diff_heat_tend = register_diag_field('ocean_model', & - 'diabatic_heat_tendency', diag%axesTL, Time, & - 'Diabatic diffusion heat tendency', & - 'Watts/m2',cmor_field_name='opottempdiff', cmor_units='W m-2', & - cmor_standard_name= & - 'tendency_of_sea_water_potential_temperature_expressed_as_heat_content_due_to_parameterized_dianeutral_mixing',& - cmor_long_name = & - 'Tendency of sea water potential temperature expressed as heat content due to parameterized dianeutral mixing') - if(CS%id_diabatic_diff_heat_tend > 0) then - CS%diabatic_diff_tendency_diag = .true. - endif - - CS%id_diabatic_diff_salt_tend = register_diag_field('ocean_model', & - 'diabatic_salt_tendency', diag%axesTL, Time, & - 'Diabatic diffusion of salt tendency', & - 'kg/(m2 * s)',cmor_field_name='osaltdiff', cmor_units='kg m-2 s-1', & - cmor_standard_name= & - 'tendency_of_sea_water_salinity_expressed_as_salt_content_due_to_parameterized_dianeutral_mixing', & - cmor_long_name = & - 'Tendency of sea water salinity expressed as salt content due to parameterized dianeutral mixing') - if(CS%id_diabatic_diff_salt_tend > 0) then - CS%diabatic_diff_tendency_diag = .true. - endif + CS%id_diabatic_diff_saln_tend = register_diag_field('ocean_model',& + 'diabatic_diff_saln_tendency', diag%axesTL, Time, & + 'Diabatic diffusion salinity tendency', 'PPT per second') + if(CS%id_diabatic_diff_saln_tend > 0) then + CS%diabatic_diff_tendency_diag = .true. + endif - ! This diagnostic should equal to roundoff if all is working well. - CS%id_diabatic_diff_heat_tend_2d = register_diag_field('ocean_model', & - 'diabatic_heat_tendency_2d', diag%axesT1, Time, & - 'Depth integrated diabatic diffusion heat tendency', & - 'Watts/m2',cmor_field_name='opottempdiff_2d', cmor_units='W m-2', & - cmor_standard_name= & - 'tendency_of_sea_water_potential_temperature_expressed_as_heat_content_due_to_parameterized_dianeutral_mixing_depth_integrated',& - cmor_long_name = & - 'Tendency of sea water potential temperature expressed as heat content due to parameterized dianeutral mixing depth integrated') - if(CS%id_diabatic_diff_heat_tend_2d > 0) then - CS%diabatic_diff_tendency_diag = .true. - endif + CS%id_diabatic_diff_heat_tend = register_diag_field('ocean_model', & + 'diabatic_heat_tendency', diag%axesTL, Time, & + 'Diabatic diffusion heat tendency', & + 'Watts/m2',cmor_field_name='opottempdiff', cmor_units='W m-2', & + cmor_standard_name= & + 'tendency_of_sea_water_potential_temperature_expressed_as_heat_content_due_to_parameterized_dianeutral_mixing',& + cmor_long_name = & + 'Tendency of sea water potential temperature expressed as heat content due to parameterized dianeutral mixing') + if(CS%id_diabatic_diff_heat_tend > 0) then + CS%diabatic_diff_tendency_diag = .true. + endif - ! This diagnostic should equal to roundoff if all is working well. - CS%id_diabatic_diff_salt_tend_2d = register_diag_field('ocean_model', & - 'diabatic_salt_tendency_2d', diag%axesT1, Time, & - 'Depth integrated diabatic diffusion salt tendency', & - 'kg/(m2 * s)',cmor_field_name='osaltdiff_2d', cmor_units='kg m-2 s-1', & - cmor_standard_name= & - 'tendency_of_sea_water_salinity_expressed_as_salt_content_due_to_parameterized_dianeutral_mixing_depth_integrated',& - cmor_long_name = & - 'Tendency of sea water salinity expressed as salt content due to parameterized dianeutral mixing depth integrated') - if(CS%id_diabatic_diff_salt_tend_2d > 0) then - CS%diabatic_diff_tendency_diag = .true. - endif + CS%id_diabatic_diff_salt_tend = register_diag_field('ocean_model', & + 'diabatic_salt_tendency', diag%axesTL, Time, & + 'Diabatic diffusion of salt tendency', & + 'kg/(m2 * s)',cmor_field_name='osaltdiff', cmor_units='kg m-2 s-1', & + cmor_standard_name= & + 'tendency_of_sea_water_salinity_expressed_as_salt_content_due_to_parameterized_dianeutral_mixing', & + cmor_long_name = & + 'Tendency of sea water salinity expressed as salt content due to parameterized dianeutral mixing') + if(CS%id_diabatic_diff_salt_tend > 0) then + CS%diabatic_diff_tendency_diag = .true. + endif + ! This diagnostic should equal to roundoff if all is working well. + CS%id_diabatic_diff_heat_tend_2d = register_diag_field('ocean_model', & + 'diabatic_heat_tendency_2d', diag%axesT1, Time, & + 'Depth integrated diabatic diffusion heat tendency', & + 'Watts/m2',cmor_field_name='opottempdiff_2d', cmor_units='W m-2', & + cmor_standard_name= & + 'tendency_of_sea_water_potential_temperature_expressed_as_heat_content_due_to_parameterized_dianeutral_mixing_depth_integrated',& + cmor_long_name = & + 'Tendency of sea water potential temperature expressed as heat content due to parameterized dianeutral mixing depth integrated') + if(CS%id_diabatic_diff_heat_tend_2d > 0) then + CS%diabatic_diff_tendency_diag = .true. + endif - ! diagnostics for tendencies of temp and saln due to boundary forcing; - ! available only for ALE algorithm. + ! This diagnostic should equal to roundoff if all is working well. + CS%id_diabatic_diff_salt_tend_2d = register_diag_field('ocean_model', & + 'diabatic_salt_tendency_2d', diag%axesT1, Time, & + 'Depth integrated diabatic diffusion salt tendency', & + 'kg/(m2 * s)',cmor_field_name='osaltdiff_2d', cmor_units='kg m-2 s-1', & + cmor_standard_name= & + 'tendency_of_sea_water_salinity_expressed_as_salt_content_due_to_parameterized_dianeutral_mixing_depth_integrated',& + cmor_long_name = & + 'Tendency of sea water salinity expressed as salt content due to parameterized dianeutral mixing depth integrated') + if(CS%id_diabatic_diff_salt_tend_2d > 0) then + CS%diabatic_diff_tendency_diag = .true. + endif - CS%id_boundary_forcing_temp_tend = register_diag_field('ocean_model',& - 'boundary_forcing_temp_tendency', diag%axesTL, Time, & - 'Boundary forcing temperature tendency', 'Degree C per second') - if(CS%id_boundary_forcing_temp_tend > 0) then - CS%boundary_forcing_tendency_diag = .true. - endif + ! diagnostics for tendencies of temp and saln due to boundary forcing; + ! available only for ALE algorithm. + CS%id_boundary_forcing_temp_tend = register_diag_field('ocean_model',& + 'boundary_forcing_temp_tendency', diag%axesTL, Time, & + 'Boundary forcing temperature tendency', 'Degree C per second') + if(CS%id_boundary_forcing_temp_tend > 0) then + CS%boundary_forcing_tendency_diag = .true. + endif - CS%id_boundary_forcing_saln_tend = register_diag_field('ocean_model',& - 'boundary_forcing_saln_tendency', diag%axesTL, Time, & - 'Boundary forcing saln tendency', 'PPT per second') - if(CS%id_boundary_forcing_saln_tend > 0) then - CS%boundary_forcing_tendency_diag = .true. - endif + CS%id_boundary_forcing_saln_tend = register_diag_field('ocean_model',& + 'boundary_forcing_saln_tendency', diag%axesTL, Time, & + 'Boundary forcing saln tendency', 'PPT per second') + if(CS%id_boundary_forcing_saln_tend > 0) then + CS%boundary_forcing_tendency_diag = .true. + endif - CS%id_boundary_forcing_heat_tend = register_diag_field('ocean_model',& - 'boundary_forcing_heat_tendency', diag%axesTL, Time, & - 'Boundary forcing heat tendency','Watts/m2') - if(CS%id_boundary_forcing_heat_tend > 0) then - CS%boundary_forcing_tendency_diag = .true. - endif + CS%id_boundary_forcing_heat_tend = register_diag_field('ocean_model',& + 'boundary_forcing_heat_tendency', diag%axesTL, Time, & + 'Boundary forcing heat tendency','Watts/m2') + if(CS%id_boundary_forcing_heat_tend > 0) then + CS%boundary_forcing_tendency_diag = .true. + endif - CS%id_boundary_forcing_salt_tend = register_diag_field('ocean_model',& - 'boundary_forcing_salt_tendency', diag%axesTL, Time, & - 'Boundary forcing salt tendency','kg m-2 s-1') - if(CS%id_boundary_forcing_salt_tend > 0) then - CS%boundary_forcing_tendency_diag = .true. - endif + CS%id_boundary_forcing_salt_tend = register_diag_field('ocean_model',& + 'boundary_forcing_salt_tendency', diag%axesTL, Time, & + 'Boundary forcing salt tendency','kg m-2 s-1') + if(CS%id_boundary_forcing_salt_tend > 0) then + CS%boundary_forcing_tendency_diag = .true. + endif - ! This diagnostic should equal to surface heat flux if all is working well. - CS%id_boundary_forcing_heat_tend_2d = register_diag_field('ocean_model',& - 'boundary_forcing_heat_tendency_2d', diag%axesT1, Time, & - 'Depth integrated boundary forcing of ocean heat','Watts/m2') - if(CS%id_boundary_forcing_heat_tend_2d > 0) then - CS%boundary_forcing_tendency_diag = .true. - endif + ! This diagnostic should equal to surface heat flux if all is working well. + CS%id_boundary_forcing_heat_tend_2d = register_diag_field('ocean_model',& + 'boundary_forcing_heat_tendency_2d', diag%axesT1, Time, & + 'Depth integrated boundary forcing of ocean heat','Watts/m2') + if(CS%id_boundary_forcing_heat_tend_2d > 0) then + CS%boundary_forcing_tendency_diag = .true. + endif - ! This diagnostic should equal to surface salt flux if all is working well. - CS%id_boundary_forcing_salt_tend_2d = register_diag_field('ocean_model',& - 'boundary_forcing_salt_tendency_2d', diag%axesT1, Time, & - 'Depth integrated boundary forcing of ocean salt','kg m-2 s-1') - if(CS%id_boundary_forcing_salt_tend_2d > 0) then - CS%boundary_forcing_tendency_diag = .true. + ! This diagnostic should equal to surface salt flux if all is working well. + CS%id_boundary_forcing_salt_tend_2d = register_diag_field('ocean_model',& + 'boundary_forcing_salt_tendency_2d', diag%axesT1, Time, & + 'Depth integrated boundary forcing of ocean salt','kg m-2 s-1') + if(CS%id_boundary_forcing_salt_tend_2d > 0) then + CS%boundary_forcing_tendency_diag = .true. + endif endif - ! diagnostics for tendencies of temp and heat due to frazil ! diagnostic for tendency of temp due to frazil @@ -2167,7 +2168,7 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, ! initialize the auxiliary diabatic driver module call diabatic_aux_init(Time, G, GV, param_file, diag, CS%diabatic_aux_CSp, & - CS%use_energetic_PBL) + CS%useALEalgorithm, CS%use_energetic_PBL) ! initialize the boundary layer modules if (CS%bulkmixedlayer) & diff --git a/src/parameterizations/vertical/MOM_diapyc_energy_req.F90 b/src/parameterizations/vertical/MOM_diapyc_energy_req.F90 index f1f614f516..2713a99a96 100644 --- a/src/parameterizations/vertical/MOM_diapyc_energy_req.F90 +++ b/src/parameterizations/vertical/MOM_diapyc_energy_req.F90 @@ -184,7 +184,7 @@ subroutine diapyc_energy_req_calc(h_in, T_in, S_in, Kd, energy_Kd, dt, tv, G, GV nz = G%ke h_neglect = GV%H_subroundoff - I_G_Earth = 1.0 / G%G_earth + I_G_Earth = 1.0 / GV%g_Earth surface_BL = .true. ; bottom_BL = .true. ; debug = .true. dPEa_dKd(:) = 0.0 ; dPEa_dKd_est(:) = 0.0 ; dPEa_dKd_err(:) = 0.0 ; dPEa_dKd_err_norm(:) = 0.0 ; dPEa_dKd_trunc(:) = 0.0 @@ -195,7 +195,7 @@ subroutine diapyc_energy_req_calc(h_in, T_in, S_in, Kd, energy_Kd, dt, tv, G, GV T0(k) = T_in(k) ; S0(k) = S_in(k) h_tr(k) = h_in(k) htot = htot + h_tr(k) - pres(K+1) = pres(K) + G%g_Earth * GV%H_to_kg_m2 * h_tr(k) + pres(K+1) = pres(K) + GV%g_Earth * GV%H_to_kg_m2 * h_tr(k) p_lay(k) = 0.5*(pres(K) + pres(K+1)) enddo do k=1,nz diff --git a/src/parameterizations/vertical/MOM_diffConvection.F90 b/src/parameterizations/vertical/MOM_diffConvection.F90 index acafdba6ee..e7ba6e668e 100644 --- a/src/parameterizations/vertical/MOM_diffConvection.F90 +++ b/src/parameterizations/vertical/MOM_diffConvection.F90 @@ -115,7 +115,7 @@ subroutine diffConvection_calculate(CS, G, GV, h, Temp, Salt, EOS, Kd_int) real, dimension( G%ke+1 ) :: Kd_1d ! Vertical diffusivity at interfaces (m2/s) real :: GoRho, pRef, rhoK, rhoKm1 - GoRho = G%g_Earth / GV%Rho0 + GoRho = GV%g_Earth / GV%Rho0 N2_1d( 1 ) = 0. N2_1d( G%ke+1 ) = 0. @@ -127,7 +127,7 @@ subroutine diffConvection_calculate(CS, G, GV, h, Temp, Salt, EOS, Kd_int) pRef = 0. ! Ignore atmospheric pressure do K = 2, G%ke ! Pressure at interface K is incremented by mass of level above - pRef = pRef + G%g_Earth * GV%Rho0 * h(i,j,k-1) * GV%H_to_m ! Boussinesq approximation!!!! ????? + pRef = pRef + GV%g_Earth * GV%Rho0 * h(i,j,k-1) * GV%H_to_m ! Boussinesq approximation!!!! ????? ! Compute Brunt-Vaisala frequency (static stability) on interfaces call calculate_density(Temp(i,j,k), Salt(i,j,k), pRef, rhoK, EOS) call calculate_density(Temp(i,j,k-1), Salt(i,j,k-1), pRef, rhoKm1, EOS) diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index f56aca2e39..b14f2389cd 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -499,7 +499,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & pres(i,1) = 0.0 do k=1,nz dMass = GV%H_to_kg_m2 * h(i,k) - dPres = G%G_Earth * dMass + dPres = GV%g_Earth * dMass dT_to_dPE(i,k) = (dMass * (pres(i,K) + 0.5*dPres)) * dSV_dT(i,j,k) dS_to_dPE(i,k) = (dMass * (pres(i,K) + 0.5*dPres)) * dSV_dS(i,j,k) dT_to_dColHt(i,k) = dMass * dSV_dT(i,j,k) diff --git a/src/parameterizations/vertical/MOM_entrain_diffusive.F90 b/src/parameterizations/vertical/MOM_entrain_diffusive.F90 index 2234be15b9..7bfc0c1070 100644 --- a/src/parameterizations/vertical/MOM_entrain_diffusive.F90 +++ b/src/parameterizations/vertical/MOM_entrain_diffusive.F90 @@ -287,7 +287,7 @@ subroutine entrainment_diffusive(u, v, h, tv, fluxes, dt, G, GV, CS, ea, eb, & H_to_m = GV%H_to_m ; m_to_H = GV%m_to_H tolerance = m_to_H * CS%Tolerance_Ent - g_2dt = 0.5 * G%g_Earth / dt + g_2dt = 0.5 * GV%g_Earth / dt kmb = GV%nk_rho_varies K2 = max(kmb+1,2) ; kb_min = K2 if (.not. CS%bulkmixedlayer) then diff --git a/src/parameterizations/vertical/MOM_internal_tide_input.F90 b/src/parameterizations/vertical/MOM_internal_tide_input.F90 index ba1f90fc06..8201e19962 100644 --- a/src/parameterizations/vertical/MOM_internal_tide_input.F90 +++ b/src/parameterizations/vertical/MOM_internal_tide_input.F90 @@ -194,7 +194,7 @@ subroutine find_N2_bottom(h, tv, T_f, S_f, h2, fluxes, G, GV, N2_bot) logical :: do_i(SZI_(G)), do_any integer :: i, j, k, is, ie, js, je, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke - G_Rho0 = G%g_Earth / GV%Rho0 + G_Rho0 = GV%g_Earth / GV%Rho0 ! Find the (limited) density jump across each interface. do i=is,ie diff --git a/src/parameterizations/vertical/MOM_kappa_shear.F90 b/src/parameterizations/vertical/MOM_kappa_shear.F90 index c5d0527fe7..a67e5aed02 100644 --- a/src/parameterizations/vertical/MOM_kappa_shear.F90 +++ b/src/parameterizations/vertical/MOM_kappa_shear.F90 @@ -319,7 +319,7 @@ subroutine Calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, & new_kappa = .true. ; if (present(initialize_all)) new_kappa = initialize_all Ri_crit = CS%Rino_crit - gR0 = GV%Rho0*G%g_Earth ; g_R0 = G%g_Earth/GV%Rho0 + gR0 = GV%Rho0*GV%g_Earth ; g_R0 = GV%g_Earth/GV%Rho0 k0dt = dt*CS%kappa_0 dz_massless = 0.1*sqrt(k0dt) diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index 1cab09531f..ffe6164ecc 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -846,6 +846,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & if (CS%id_Kd_itidal > 0) call post_data(CS%id_Kd_itidal, dd%Kd_itidal, CS%diag) if (CS%id_Kd_Niku > 0) call post_data(CS%id_Kd_Niku, dd%Kd_Niku, CS%diag) if (CS%id_Kd_lowmode> 0) call post_data(CS%id_Kd_lowmode, dd%Kd_lowmode, CS%diag) + if (CS%id_Fl_lowmode> 0) call post_data(CS%id_Fl_lowmode, dd%Fl_lowmode, CS%diag) if (CS%id_Kd_user > 0) call post_data(CS%id_Kd_user, dd%Kd_user, CS%diag) if (CS%id_Kd_Work > 0) call post_data(CS%id_Kd_Work, dd%Kd_Work, CS%diag) if (CS%id_Kd_Itidal_Work > 0) & @@ -996,7 +997,7 @@ subroutine find_TKE_to_Kd(h, tv, dRho_int, N2_lay, j, dt, G, GV, CS, & I_dt = 1.0/dt Omega2 = CS%Omega**2 - G_Rho0 = G%g_Earth / GV%Rho0 + G_Rho0 = GV%g_Earth / GV%Rho0 H_neglect = GV%H_subroundoff I_Rho0 = 1.0/GV%Rho0 @@ -1121,11 +1122,11 @@ subroutine find_TKE_to_Kd(h, tv, dRho_int, N2_lay, j, dt, G, GV, CS, & ! ### This should be 1 / G_Earth * (delta rho_InSitu) ! kappa_max = I_dt * dRho_int(i,K+1) * maxEnt(i,k) * & ! (GV%H_to_m*h(i,j,k) + dh_max) / dRho_lay - ! maxTKE(i,k) = G%g_Earth * dRho_lay * kappa_max + ! maxTKE(i,k) = GV%g_Earth * dRho_lay * kappa_max ! dRho_int should already be non-negative, so the max is redundant? dh_max = maxEnt(i,k) * (1.0 + dsp1_ds(i,k)) dRho_lay = 0.5 * max(dRho_int(i,K) + dRho_int(i,K+1), 0.0) - maxTKE(i,k) = I_dt * ((G%g_Earth * I_Rho0) * & + maxTKE(i,k) = I_dt * ((GV%g_Earth * I_Rho0) * & (0.5*max(dRho_int(i,K+1) + dsp1_ds(i,k)*dRho_int(i,K),0.0))) * & ((GV%H_to_m*h(i,j,k) + dh_max) * maxEnt(i,k)) TKE_to_Kd(i,k) = 1.0 / (G_Rho0 * dRho_lay + & @@ -1172,7 +1173,7 @@ subroutine find_N2(h, tv, T_f, S_f, fluxes, j, G, GV, CS, dRho_int, & integer :: i, k, is, ie, nz is = G%isc ; ie = G%iec ; nz = G%ke - G_Rho0 = G%g_Earth / GV%Rho0 + G_Rho0 = GV%g_Earth / GV%Rho0 H_neglect = GV%H_subroundoff ! Find the (limited) density jump across each interface. @@ -1444,7 +1445,7 @@ subroutine add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, & if (associated(visc%Ray_u) .and. associated(visc%Ray_v)) Rayleigh_drag = .true. I_Rho0 = 1.0/GV%Rho0 - R0_g = GV%Rho0/G%g_Earth + R0_g = GV%Rho0/GV%g_Earth do K=2,nz ; Rint(K) = 0.5*(GV%Rlay(k-1)+GV%Rlay(k)) ; enddo @@ -2433,7 +2434,7 @@ subroutine set_density_ratios(h, tv, kb, G, GV, CS, j, ds_dsp1, rho_0) enddo if (CS%bulkmixedlayer) then - g_R0 = G%g_Earth/GV%Rho0 + g_R0 = GV%g_Earth/GV%Rho0 kmb = GV%nk_rho_varies eps = 0.1 do i=is,ie ; p_ref(i) = tv%P_Ref ; enddo @@ -3006,81 +3007,85 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp CS%Decay_scale_factor_lee = -9.e99 ! This should never be used if CS%Lee_wave_dissipation = False endif - CS%id_TKE_itidal = register_diag_field('ocean_model','TKE_itidal',diag%axesT1,Time, & - 'Internal Tide Driven Turbulent Kinetic Energy', 'Watt meter-2') - CS%id_maxTKE = register_diag_field('ocean_model','maxTKE',diag%axesTL,Time, & - 'Maximum layer TKE', 'meter3 second-3') - CS%id_TKE_to_Kd = register_diag_field('ocean_model','TKE_to_Kd',diag%axesTL,Time, & - 'Convert TKE to Kd', 'second2 meter') + if (CS%Int_tide_dissipation .or. CS%Lee_wave_dissipation .or. & + CS%Lowmode_itidal_dissipation) then - CS%id_Nb = register_diag_field('ocean_model','Nb',diag%axesT1,Time, & - 'Bottom Buoyancy Frequency', 'sec-1') + CS%id_TKE_itidal = register_diag_field('ocean_model','TKE_itidal',diag%axesT1,Time, & + 'Internal Tide Driven Turbulent Kinetic Energy', 'Watt meter-2') + CS%id_maxTKE = register_diag_field('ocean_model','maxTKE',diag%axesTL,Time, & + 'Maximum layer TKE', 'meter3 second-3') + CS%id_TKE_to_Kd = register_diag_field('ocean_model','TKE_to_Kd',diag%axesTL,Time, & + 'Convert TKE to Kd', 'second2 meter') - CS%id_Kd_itidal = register_diag_field('ocean_model','Kd_itides',diag%axesTi,Time, & - 'Internal Tide Driven Diffusivity', 'meter2 sec-1') + CS%id_Nb = register_diag_field('ocean_model','Nb',diag%axesT1,Time, & + 'Bottom Buoyancy Frequency', 'sec-1') - CS%id_Kd_lowmode = register_diag_field('ocean_model','Kd_lowmode',diag%axesTi,Time, & - 'Internal Tide Driven Diffusivity (from propagating low modes)', 'meter2 sec-1') + CS%id_Kd_itidal = register_diag_field('ocean_model','Kd_itides',diag%axesTi,Time, & + 'Internal Tide Driven Diffusivity', 'meter2 sec-1') - CS%id_Fl_itidal = register_diag_field('ocean_model','Fl_itides',diag%axesTi,Time, & - 'Vertical flux of tidal turbulent dissipation', 'meter3 sec-3') + CS%id_Kd_lowmode = register_diag_field('ocean_model','Kd_lowmode',diag%axesTi,Time, & + 'Internal Tide Driven Diffusivity (from propagating low modes)', 'meter2 sec-1') - CS%id_Fl_lowmode = register_diag_field('ocean_model','Fl_lowmode',diag%axesTi,Time, & - 'Vertical flux of tidal turbulent dissipation (from propagating low modes)', 'meter3 sec-3') + CS%id_Fl_itidal = register_diag_field('ocean_model','Fl_itides',diag%axesTi,Time, & + 'Vertical flux of tidal turbulent dissipation', 'meter3 sec-3') - CS%id_Polzin_decay_scale = register_diag_field('ocean_model','Polzin_decay_scale',diag%axesT1,Time, & - 'Vertical decay scale for the tidal turbulent dissipation with Polzin scheme', 'meter') + CS%id_Fl_lowmode = register_diag_field('ocean_model','Fl_lowmode',diag%axesTi,Time, & + 'Vertical flux of tidal turbulent dissipation (from propagating low modes)', 'meter3 sec-3') - CS%id_Polzin_decay_scale_scaled = register_diag_field('ocean_model','Polzin_decay_scale_scaled',diag%axesT1,Time, & - 'Vertical decay scale for the tidal turbulent dissipation with Polzin scheme, scaled by N2_bot/N2_meanz', 'meter') + CS%id_Polzin_decay_scale = register_diag_field('ocean_model','Polzin_decay_scale',diag%axesT1,Time, & + 'Vertical decay scale for the tidal turbulent dissipation with Polzin scheme', 'meter') - CS%id_N2_bot = register_diag_field('ocean_model','N2_b',diag%axesT1,Time, & - 'Bottom Buoyancy frequency squared', 's-2') + CS%id_Polzin_decay_scale_scaled = register_diag_field('ocean_model','Polzin_decay_scale_scaled',diag%axesT1,Time, & + 'Vertical decay scale for the tidal turbulent dissipation with Polzin scheme, scaled by N2_bot/N2_meanz', 'meter') - CS%id_N2_meanz = register_diag_field('ocean_model','N2_meanz',diag%axesT1,Time, & - 'Buoyancy frequency squared averaged over the water column', 's-2') + CS%id_N2_bot = register_diag_field('ocean_model','N2_b',diag%axesT1,Time, & + 'Bottom Buoyancy frequency squared', 's-2') - CS%id_Kd_Work = register_diag_field('ocean_model','Kd_Work',diag%axesTL,Time, & - 'Work done by Diapycnal Mixing', 'Watts m-2') + CS%id_N2_meanz = register_diag_field('ocean_model','N2_meanz',diag%axesT1,Time, & + 'Buoyancy frequency squared averaged over the water column', 's-2') - CS%id_Kd_Itidal_Work = register_diag_field('ocean_model','Kd_Itidal_Work',diag%axesTL,Time, & - 'Work done by Internal Tide Diapycnal Mixing', 'Watts m-2') + CS%id_Kd_Work = register_diag_field('ocean_model','Kd_Work',diag%axesTL,Time, & + 'Work done by Diapycnal Mixing', 'Watts m-2') - CS%id_Kd_Niku_Work = register_diag_field('ocean_model','Kd_Nikurashin_Work',diag%axesTL,Time, & - 'Work done by Nikurashin Lee Wave Drag Scheme', 'Watts m-2') + CS%id_Kd_Itidal_Work = register_diag_field('ocean_model','Kd_Itidal_Work',diag%axesTL,Time, & + 'Work done by Internal Tide Diapycnal Mixing', 'Watts m-2') - CS%id_Kd_Lowmode_Work = register_diag_field('ocean_model','Kd_Lowmode_Work',diag%axesTL,Time, & - 'Work done by Internal Tide Diapycnal Mixing (low modes)', 'Watts m-2') + CS%id_Kd_Niku_Work = register_diag_field('ocean_model','Kd_Nikurashin_Work',diag%axesTL,Time, & + 'Work done by Nikurashin Lee Wave Drag Scheme', 'Watts m-2') - CS%id_N2 = register_diag_field('ocean_model','N2',diag%axesTi,Time, & - 'Buoyancy frequency squared', 'sec-2', cmor_field_name='obvfsq', & - cmor_units='s-2', cmor_long_name='Square of seawater buoyancy frequency',& - cmor_standard_name='square_of_brunt_vaisala_frequency_in_sea_water') + CS%id_Kd_Lowmode_Work = register_diag_field('ocean_model','Kd_Lowmode_Work',diag%axesTL,Time, & + 'Work done by Internal Tide Diapycnal Mixing (low modes)', 'Watts m-2') + + CS%id_N2 = register_diag_field('ocean_model','N2',diag%axesTi,Time, & + 'Buoyancy frequency squared', 'sec-2', cmor_field_name='obvfsq', & + cmor_units='s-2', cmor_long_name='Square of seawater buoyancy frequency',& + cmor_standard_name='square_of_brunt_vaisala_frequency_in_sea_water') - if (CS%user_change_diff) & - CS%id_Kd_user = register_diag_field('ocean_model','Kd_user',diag%axesTi,Time, & - 'User-specified Extra Diffusivity', 'meter2 sec-1') - - if (associated(diag_to_Z_CSp)) then - vd = var_desc("N2", "second-2",& - "Buoyancy frequency, interpolated to z", z_grid='z') - CS%id_N2_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) - vd = var_desc("Kd_itides","meter2 second-1", & - "Internal Tide Driven Diffusivity, interpolated to z", z_grid='z') - CS%id_Kd_itidal_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) - if (CS%Lee_wave_dissipation) then - vd = var_desc("Kd_Nikurashin", "meter2 second-1", & - "Lee Wave Driven Diffusivity, interpolated to z", z_grid='z') - CS%id_Kd_Niku_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) - endif - if (CS%Lowmode_itidal_dissipation) then - vd = var_desc("Kd_lowmode","meter2 second-1", & - "Internal Tide Driven Diffusivity (from low modes), interpolated to z",& - z_grid='z') - CS%id_Kd_lowmode_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) - endif if (CS%user_change_diff) & - CS%id_Kd_user_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) + CS%id_Kd_user = register_diag_field('ocean_model','Kd_user',diag%axesTi,Time, & + 'User-specified Extra Diffusivity', 'meter2 sec-1') + + if (associated(diag_to_Z_CSp)) then + vd = var_desc("N2", "second-2",& + "Buoyancy frequency, interpolated to z", z_grid='z') + CS%id_N2_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) + vd = var_desc("Kd_itides","meter2 second-1", & + "Internal Tide Driven Diffusivity, interpolated to z", z_grid='z') + CS%id_Kd_itidal_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) + if (CS%Lee_wave_dissipation) then + vd = var_desc("Kd_Nikurashin", "meter2 second-1", & + "Lee Wave Driven Diffusivity, interpolated to z", z_grid='z') + CS%id_Kd_Niku_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) + endif + if (CS%Lowmode_itidal_dissipation) then + vd = var_desc("Kd_lowmode","meter2 second-1", & + "Internal Tide Driven Diffusivity (from low modes), interpolated to z",& + z_grid='z') + CS%id_Kd_lowmode_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) + endif + if (CS%user_change_diff) & + CS%id_Kd_user_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) + endif endif call get_param(param_file, mod, "DOUBLE_DIFFUSION", CS%double_diffusion, & diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index 6edfdb6565..92586a922f 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -287,7 +287,7 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, CS) Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB nkmb = GV%nk_rho_varies ; nkml = GV%nkml h_neglect = GV%H_subroundoff - Rho0x400_G = 400.0*(GV%Rho0/G%g_Earth)*GV%m_to_H + Rho0x400_G = 400.0*(GV%Rho0/GV%g_Earth)*GV%m_to_H Vol_quit = 0.9*GV%Angstrom + h_neglect H_to_m = GV%H_to_m ; m_to_H = GV%m_to_H C2pi_3 = 8.0*atan(1.0)/3.0 @@ -990,7 +990,7 @@ subroutine set_viscous_ML(u, v, h, tv, fluxes, visc, dt, G, GV, CS) "Module must be initialized before it is used.") if (.not.(CS%dynamic_viscous_ML .or. associated(fluxes%frac_shelf_h))) return - Rho0x400_G = 400.0*(GV%Rho0/G%g_Earth)*GV%m_to_H + Rho0x400_G = 400.0*(GV%Rho0/GV%g_Earth)*GV%m_to_H U_bg_sq = CS%drag_bg_vel * CS%drag_bg_vel cdrag_sqrt=sqrt(CS%cdrag) @@ -998,7 +998,7 @@ subroutine set_viscous_ML(u, v, h, tv, fluxes, visc, dt, G, GV, CS) dt_Rho0 = dt/GV%H_to_kg_m2 h_neglect = GV%H_subroundoff h_tiny = 2.0*GV%Angstrom + h_neglect - g_H_Rho0 = (G%g_Earth * GV%H_to_m) / GV%Rho0 + g_H_Rho0 = (GV%g_Earth * GV%H_to_m) / GV%Rho0 H_to_m = GV%H_to_m ; m_to_H = GV%m_to_H if (associated(fluxes%frac_shelf_h)) then diff --git a/src/parameterizations/vertical/MOM_shortwave_abs.F90 b/src/parameterizations/vertical/MOM_shortwave_abs.F90 index 54e1c53d36..e384efba55 100644 --- a/src/parameterizations/vertical/MOM_shortwave_abs.F90 +++ b/src/parameterizations/vertical/MOM_shortwave_abs.F90 @@ -174,7 +174,7 @@ subroutine absorbRemainingSW(G, GV, h, opacity_band, nsw, j, dt, H_limit_fluxes, C1_6 = 1.0 / 6.0 ; C1_60 = 1.0 / 60.0 TKE_calc = (present(TKE) .and. present(dSV_dT)) - g_Hconv2 = G%G_earth * GV%H_to_kg_m2**2 + g_Hconv2 = GV%g_Earth * GV%H_to_kg_m2**2 h_heat(:) = 0.0 if (present(htot)) then ; do i=is,ie ; h_heat(i) = htot(i) ; enddo ; endif diff --git a/src/tracer/ISOMIP_tracer.F90 b/src/tracer/ISOMIP_tracer.F90 index 44d985ae05..6b826f40d5 100644 --- a/src/tracer/ISOMIP_tracer.F90 +++ b/src/tracer/ISOMIP_tracer.F90 @@ -167,13 +167,15 @@ end function register_ISOMIP_tracer !> Initializes the NTR tracer fields in tr(:,:,:,:) ! and it sets up the tracer output. -subroutine initialize_ISOMIP_tracer(restart, day, G, GV, h, OBC, CS, ALE_sponge_CSp, & - diag_to_Z_CSp) +subroutine initialize_ISOMIP_tracer(restart, day, G, GV, h, diag, OBC, CS, & + ALE_sponge_CSp, diag_to_Z_CSp) + type(ocean_grid_type), intent(in) :: G !< Grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. logical, intent(in) :: restart !< .true. if the fields have already been read from a restart file. type(time_type), target, intent(in) :: day !< Time of the start of the run. real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness, in m or kg m-2. + type(diag_ctrl), target, intent(in) :: diag type(ocean_OBC_type), pointer :: OBC !< This open boundary condition type specifies whether, where, and what open boundary conditions are used. This is not being used for now. type(ISOMIP_tracer_CS), pointer :: CS !< The control structure returned by a previous call to ISOMIP_register_tracer. type(ALE_sponge_CS), pointer :: ALE_sponge_CSp !< A pointer to the control structure for the sponges, if they are in use. Otherwise this may be unassociated. @@ -207,6 +209,7 @@ subroutine initialize_ISOMIP_tracer(restart, day, G, GV, h, OBC, CS, ALE_sponge_ h_neglect = GV%H_subroundoff CS%Time => day + CS%diag => diag if (.not.restart) then if (len_trim(CS%tracer_IC_file) >= 1) then @@ -341,7 +344,7 @@ subroutine ISOMIP_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, G if (fluxes%iceshelf_melt(i,j) > 0.0) then ! CS%tr(i,j,1,m) = (dt*fluxes%iceshelf_melt(i,j)/(365.0 * 86400.0) & ! + h_old(i,j,1)*CS%tr(i,j,1,m))/h_new(i,j,1) - CS%tr(i,j,1,m) = 1.0 + CS%tr(i,j,1:3,m) = 1.0 endif enddo ; enddo enddo diff --git a/src/tracer/MOM_tracer_flow_control.F90 b/src/tracer/MOM_tracer_flow_control.F90 index bc0127faf5..3d43765d9d 100644 --- a/src/tracer/MOM_tracer_flow_control.F90 +++ b/src/tracer/MOM_tracer_flow_control.F90 @@ -257,7 +257,7 @@ subroutine tracer_flow_control_init(restart, day, G, GV, h, param_file, diag, OB call initialize_DOME_tracer(restart, day, G, GV, h, diag, OBC, CS%DOME_tracer_CSp, & sponge_CSp, diag_to_Z_CSp) if (CS%use_ISOMIP_tracer) & - call initialize_ISOMIP_tracer(restart, day, G, GV, h, OBC, CS%ISOMIP_tracer_CSp, & + call initialize_ISOMIP_tracer(restart, day, G, GV, h, diag, OBC, CS%ISOMIP_tracer_CSp, & ALE_sponge_CSp, diag_to_Z_CSp) if (CS%use_ideal_age) & call initialize_ideal_age_tracer(restart, day, G, GV, h, diag, OBC, CS%ideal_age_tracer_CSp, & diff --git a/src/tracer/MOM_tracer_hor_diff.F90 b/src/tracer/MOM_tracer_hor_diff.F90 index 73b581d299..33dd6c55fb 100644 --- a/src/tracer/MOM_tracer_hor_diff.F90 +++ b/src/tracer/MOM_tracer_hor_diff.F90 @@ -1397,8 +1397,10 @@ subroutine tracer_hor_diff_init(Time, G, param_file, diag, CS, CSnd) cmor_field_name='diftrelo', cmor_units='m2 sec-1', & cmor_standard_name= 'ocean_tracer_epineutral_laplacian_diffusivity', & cmor_long_name = 'Ocean Tracer Epineutral Laplacian Diffusivity') - CS%id_CFL = register_diag_field('ocean_model', 'CFL_lateral_diff', diag%axesT1, Time,& - 'Grid CFL number for lateral/neutral tracer diffusion', 'dimensionless') + if (CS%check_diffusive_CFL) then + CS%id_CFL = register_diag_field('ocean_model', 'CFL_lateral_diff', diag%axesT1, Time,& + 'Grid CFL number for lateral/neutral tracer diffusion', 'dimensionless') + endif end subroutine tracer_hor_diff_init diff --git a/src/user/DOME_initialization.F90 b/src/user/DOME_initialization.F90 index 3e6e1a205c..c1487cc3f9 100644 --- a/src/user/DOME_initialization.F90 +++ b/src/user/DOME_initialization.F90 @@ -359,7 +359,7 @@ subroutine DOME_set_Open_Bdry_Conds(OBC, tv, G, GV, param_file, tr_Reg) endif if (apply_OBC_v) then - g_prime_tot = (G%g_Earth/GV%Rho0)*2.0 + g_prime_tot = (GV%g_Earth/GV%Rho0)*2.0 Def_Rad = sqrt(D_edge*g_prime_tot) / (1.0e-4*1000.0) tr_0 = (-D_edge*sqrt(D_edge*g_prime_tot)*0.5e3*Def_Rad) * GV%m_to_H diff --git a/src/user/Rossby_front_2d_initialization.F90 b/src/user/Rossby_front_2d_initialization.F90 index 149ff3db25..29b3e7adf5 100644 --- a/src/user/Rossby_front_2d_initialization.F90 +++ b/src/user/Rossby_front_2d_initialization.F90 @@ -166,7 +166,7 @@ subroutine Rossby_front_initialize_velocity(u, v, h, G, GV, param_file) do j = G%jsc,G%jec ; do I = G%isc-1,G%iec+1 f = 0.5*( G%CoriolisBu(I,j) + G%CoriolisBu(I,j-1) ) dUdT = 0.0 ; if (abs(f) > 0.0) & - dUdT = ( G%g_Earth * dRho_dT ) / ( f * GV%Rho0 ) + dUdT = ( GV%g_Earth * dRho_dT ) / ( f * GV%Rho0 ) Dml = Hml( G, G%geoLatT(i,j) ) Ty = dTdy( G, T_range, G%geoLatT(i,j) ) zi = 0.