Skip to content

Commit

Permalink
+Use GV%g_Earth in place of G%g_Earth
Browse files Browse the repository at this point in the history 
  Replaced G%g_Earth with GV%g_Earth whereever possible, and read G_EARTH only
from the vertical grid module.  Also reordered the reading of BULKMIXEDLAYER and
USE_REGRIDDING, and changed the default value of BULKMIXEDLAYER to false when
USE_REGRIDDING is true.  (Default settings should not lead to fatal errors!)
All answers are bitwise identical in configurations that would actually work,
although one default value sometimes changes and most MOM_parameter_doc files
are changed.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Jun 29, 2016
1 parent 671aa5f commit 3e80bcb
Show file tree
Hide file tree
Showing 29 changed files with 153 additions and 131 deletions.
57 changes: 31 additions & 26 deletions src/core/MOM.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -690,7 +690,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

Expand Down Expand Up @@ -984,7 +984,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

Expand Down Expand Up @@ -1184,7 +1184,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
Expand Down Expand Up @@ -1227,7 +1227,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)
Expand Down Expand Up @@ -1444,6 +1444,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.)
Expand All @@ -1461,27 +1462,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.)
Expand Down Expand Up @@ -1624,14 +1620,20 @@ 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, &
"MOM: ENABLE_THERMODYNAMICS must be defined to use USE_EOS.")
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.
Expand Down Expand Up @@ -1661,13 +1663,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)
Expand Down Expand Up @@ -2051,9 +2056,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)
Expand Down Expand Up @@ -2098,7 +2103,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)

Expand Down Expand Up @@ -2843,7 +2848,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
Expand Down
14 changes: 7 additions & 7 deletions src/core/MOM_PressureForce_Montgomery.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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)

Expand Down
26 changes: 13 additions & 13 deletions src/core/MOM_PressureForce_analytic_FV.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand All @@ -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

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)

Expand Down
2 changes: 1 addition & 1 deletion src/core/MOM_dynamics_legacy_split.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -475,7 +475,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
Expand Down
2 changes: 1 addition & 1 deletion src/core/MOM_forcing_type.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
Loading

0 comments on commit 3e80bcb

Please sign in to comment.