Skip to content

Commit

Permalink
+Rescaled the units of surface mass fluxes
Browse files Browse the repository at this point in the history 
  Rescaled the units of 7 surface mass flux elements in the forcing type,
including lprec, fprec, lrunoff, frunoff, evap, vprec, and seaice_melt.  Also
added a number of unit_scaling_type arguments to subroutines to enable this
rescaling.  All answers are bitwise identical, but the units of 7 elements of
a widely used public type have changed and there are new subroutine arguments.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Oct 4, 2019
1 parent 34193e9 commit dbad998
Show file tree
Hide file tree
Showing 23 changed files with 281 additions and 234 deletions.
37 changes: 22 additions & 15 deletions config_src/coupled_driver/MOM_surface_forcing_gfdl.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -242,6 +242,8 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, US, CS, sfc
real :: delta_sss ! temporary storage for sss diff from restoring value [ppt]
real :: delta_sst ! temporary storage for sst diff from restoring value [degC]

real :: kg_m2_s_conversion !< A combination of unit conversion factors for rescaling
!! mass fluxes [R Z s m2 kg-1 T-1 ~> 1].
real :: C_p ! heat capacity of seawater [J degC-1 kg-1]
real :: sign_for_net_FW_bug ! Should be +1. but an old bug can be recovered by using -1.

Expand All @@ -255,6 +257,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, US, CS, sfc
IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB
isr = is-isd+1 ; ier = ie-isd+1 ; jsr = js-jsd+1 ; jer = je-jsd+1

kg_m2_s_conversion = US%kg_m3_to_R*US%m_to_Z*US%T_to_s
C_p = fluxes%C_p
open_ocn_mask(:,:) = 1.0
pme_adj(:,:) = 0.0
Expand Down Expand Up @@ -372,19 +375,21 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, US, CS, sfc
delta_sss = sfc_state%SSS(i,j) - data_restore(i,j)
delta_sss = sign(1.0,delta_sss)*min(abs(delta_sss),CS%max_delta_srestore)
fluxes%vprec(i,j) = (CS%basin_mask(i,j)*open_ocn_mask(i,j)*CS%srestore_mask(i,j))* &
(US%R_to_kg_m3*CS%Rho0*CS%Flux_const) * &
(US%m_to_Z*US%T_to_s*CS%Rho0*CS%Flux_const) * &
delta_sss / (0.5*(sfc_state%SSS(i,j) + data_restore(i,j)))
endif
enddo ; enddo
if (CS%adjust_net_srestore_to_zero) then
if (CS%adjust_net_srestore_by_scaling) then
call adjust_area_mean_to_zero(fluxes%vprec, G, fluxes%vPrecGlobalScl)
call adjust_area_mean_to_zero(fluxes%vprec, G, fluxes%vPrecGlobalScl, &
unit_scale=US%R_to_kg_m3*US%Z_to_m*US%s_to_T)
fluxes%vPrecGlobalAdj = 0.
else
work_sum(is:ie,js:je) = US%L_to_m**2*G%areaT(is:ie,js:je)*fluxes%vprec(is:ie,js:je)
work_sum(is:ie,js:je) = US%L_to_m**2*G%areaT(is:ie,js:je) * &
US%R_to_kg_m3*US%Z_to_m*US%s_to_T*fluxes%vprec(is:ie,js:je)
fluxes%vPrecGlobalAdj = reproducing_sum(work_sum(:,:), isr, ier, jsr, jer) / CS%area_surf
do j=js,je ; do i=is,ie
fluxes%vprec(i,j) = ( fluxes%vprec(i,j) - fluxes%vPrecGlobalAdj ) * G%mask2dT(i,j)
fluxes%vprec(i,j) = ( fluxes%vprec(i,j) - kg_m2_s_conversion*fluxes%vPrecGlobalAdj ) * G%mask2dT(i,j)
enddo ; enddo
endif
endif
Expand All @@ -408,31 +413,31 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, US, CS, sfc
do j=js,je ; do i=is,ie

if (associated(IOB%lprec)) then
fluxes%lprec(i,j) = IOB%lprec(i-i0,j-j0) * G%mask2dT(i,j)
fluxes%lprec(i,j) = kg_m2_s_conversion * IOB%lprec(i-i0,j-j0) * G%mask2dT(i,j)
if (CS%check_no_land_fluxes) &
call check_mask_val_consistency(IOB%lprec(i-i0,j-j0), G%mask2dT(i,j), i, j, 'lprec', G)
endif

if (associated(IOB%fprec)) then
fluxes%fprec(i,j) = IOB%fprec(i-i0,j-j0) * G%mask2dT(i,j)
fluxes%fprec(i,j) = kg_m2_s_conversion * IOB%fprec(i-i0,j-j0) * G%mask2dT(i,j)
if (CS%check_no_land_fluxes) &
call check_mask_val_consistency(IOB%fprec(i-i0,j-j0), G%mask2dT(i,j), i, j, 'fprec', G)
endif

if (associated(IOB%q_flux)) then
fluxes%evap(i,j) = - IOB%q_flux(i-i0,j-j0) * G%mask2dT(i,j)
fluxes%evap(i,j) = - kg_m2_s_conversion * IOB%q_flux(i-i0,j-j0) * G%mask2dT(i,j)
if (CS%check_no_land_fluxes) &
call check_mask_val_consistency(IOB%q_flux(i-i0,j-j0), G%mask2dT(i,j), i, j, 'q_flux', G)
endif

if (associated(IOB%runoff)) then
fluxes%lrunoff(i,j) = IOB%runoff(i-i0,j-j0) * G%mask2dT(i,j)
fluxes%lrunoff(i,j) = kg_m2_s_conversion * IOB%runoff(i-i0,j-j0) * G%mask2dT(i,j)
if (CS%check_no_land_fluxes) &
call check_mask_val_consistency(IOB%runoff(i-i0,j-j0), G%mask2dT(i,j), i, j, 'runoff', G)
endif

if (associated(IOB%calving)) then
fluxes%frunoff(i,j) = IOB%calving(i-i0,j-j0) * G%mask2dT(i,j)
fluxes%frunoff(i,j) = kg_m2_s_conversion * IOB%calving(i-i0,j-j0) * G%mask2dT(i,j)
if (CS%check_no_land_fluxes) &
call check_mask_val_consistency(IOB%calving(i-i0,j-j0), G%mask2dT(i,j), i, j, 'calving', G)
endif
Expand Down Expand Up @@ -565,7 +570,8 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, US, CS, sfc
sign_for_net_FW_bug = 1.
if (CS%use_net_FW_adjustment_sign_bug) sign_for_net_FW_bug = -1.
do j=js,je ; do i=is,ie
net_FW(i,j) = (((fluxes%lprec(i,j) + fluxes%fprec(i,j)) + &
net_FW(i,j) = US%R_to_kg_m3*US%Z_to_m*US%s_to_T* &
(((fluxes%lprec(i,j) + fluxes%fprec(i,j)) + &
(fluxes%lrunoff(i,j) + fluxes%frunoff(i,j))) + &
(fluxes%evap(i,j) + fluxes%vprec(i,j)) ) * US%L_to_m**2*G%areaT(i,j)
! The following contribution appears to be calculating the volume flux of sea-ice
Expand All @@ -583,13 +589,13 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, US, CS, sfc
if (CS%adjust_net_fresh_water_by_scaling) then
call adjust_area_mean_to_zero(net_FW2, G, fluxes%netFWGlobalScl)
do j=js,je ; do i=is,ie
fluxes%vprec(i,j) = fluxes%vprec(i,j) + &
fluxes%vprec(i,j) = fluxes%vprec(i,j) + US%kg_m3_to_R*US%m_to_Z*US%T_to_s * &
(net_FW2(i,j) - net_FW(i,j)/(US%L_to_m**2*G%areaT(i,j))) * G%mask2dT(i,j)
enddo ; enddo
else
fluxes%netFWGlobalAdj = reproducing_sum(net_FW(:,:), isr, ier, jsr, jer) / CS%area_surf
do j=js,je ; do i=is,ie
fluxes%vprec(i,j) = ( fluxes%vprec(i,j) - fluxes%netFWGlobalAdj ) * G%mask2dT(i,j)
fluxes%vprec(i,j) = ( fluxes%vprec(i,j) - kg_m2_s_conversion * fluxes%netFWGlobalAdj ) * G%mask2dT(i,j)
enddo ; enddo
endif

Expand All @@ -611,7 +617,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, US, CS, sfc

if (CS%allow_flux_adjustments) then
! Apply adjustments to fluxes
call apply_flux_adjustments(G, CS, Time, fluxes)
call apply_flux_adjustments(G, US, CS, Time, fluxes)
endif

! Allow for user-written code to alter fluxes after all the above
Expand Down Expand Up @@ -1086,8 +1092,9 @@ end subroutine extract_IOB_stresses
!! - hflx_adj (Heat flux into the ocean [W m-2])
!! - sflx_adj (Salt flux into the ocean [kg salt m-2 s-1])
!! - prcme_adj (Fresh water flux into the ocean [kg m-2 s-1])
subroutine apply_flux_adjustments(G, CS, Time, fluxes)
subroutine apply_flux_adjustments(G, US, CS, Time, fluxes)
type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
type(surface_forcing_CS), pointer :: CS !< Surface forcing control structure
type(time_type), intent(in) :: Time !< Model time structure
type(forcing), intent(inout) :: fluxes !< Surface fluxes structure
Expand Down Expand Up @@ -1120,7 +1127,7 @@ subroutine apply_flux_adjustments(G, CS, Time, fluxes)
call data_override('OCN', 'prcme_adj', temp_at_h(isc:iec,jsc:jec), Time, override=overrode_h)

if (overrode_h) then ; do j=jsc,jec ; do i=isc,iec
fluxes%vprec(i,j) = fluxes%vprec(i,j) + temp_at_h(i,j)* G%mask2dT(i,j)
fluxes%vprec(i,j) = fluxes%vprec(i,j) + US%kg_m3_to_R*US%m_to_Z*US%T_to_s * temp_at_h(i,j)* G%mask2dT(i,j)
enddo ; enddo ; endif
! Not needed? ! if (overrode_h) call pass_var(fluxes%vprec, G%Domain)
end subroutine apply_flux_adjustments
Expand Down
8 changes: 4 additions & 4 deletions config_src/coupled_driver/ocean_model_MOM.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -520,7 +520,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_upda
if (OS%use_ice_shelf) &
call shelf_calc_flux(OS%sfc_state, OS%fluxes, OS%Time, dt_coupling, OS%Ice_shelf_CSp)
if (OS%icebergs_alter_ocean) &
call iceberg_fluxes(OS%grid, OS%fluxes, OS%use_ice_shelf, &
call iceberg_fluxes(OS%grid, OS%US, OS%fluxes, OS%use_ice_shelf, &
OS%sfc_state, dt_coupling, OS%marine_ice_CSp)

#ifdef _USE_GENERIC_TRACER
Expand All @@ -541,7 +541,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_upda
if (OS%use_ice_shelf) &
call shelf_calc_flux(OS%sfc_state, OS%flux_tmp, OS%Time, dt_coupling, OS%Ice_shelf_CSp)
if (OS%icebergs_alter_ocean) &
call iceberg_fluxes(OS%grid, OS%flux_tmp, OS%use_ice_shelf, &
call iceberg_fluxes(OS%grid, OS%US, OS%flux_tmp, OS%use_ice_shelf, &
OS%sfc_state, dt_coupling, OS%marine_ice_CSp)

call fluxes_accumulate(OS%flux_tmp, OS%fluxes, dt_coupling, OS%grid, weight)
Expand All @@ -554,7 +554,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_upda

! The net mass forcing is not currently used in the MOM6 dynamics solvers, so this is may be unnecessary.
if (do_dyn .and. associated(OS%forces%net_mass_src) .and. .not.OS%forces%net_mass_src_set) &
call get_net_mass_forcing(OS%fluxes, OS%grid, OS%forces%net_mass_src)
call get_net_mass_forcing(OS%fluxes, OS%grid, OS%US, OS%forces%net_mass_src)

if (OS%use_waves .and. do_thermo) then
! For now, the waves are only updated on the thermodynamics steps, because that is where
Expand Down Expand Up @@ -654,7 +654,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_upda
if (OS%fluxes%fluxes_used .and. do_thermo) then
call enable_averaging(OS%fluxes%dt_buoy_accum, OS%Time, OS%diag)
call forcing_diagnostics(OS%fluxes, OS%sfc_state, OS%fluxes%dt_buoy_accum, &
OS%grid, OS%diag, OS%forcing_CSp%handles)
OS%grid, OS%US, OS%diag, OS%forcing_CSp%handles)
call disable_averaging(OS%diag)
endif

Expand Down
23 changes: 12 additions & 11 deletions config_src/ice_solo_driver/MOM_surface_forcing.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -275,7 +275,7 @@ subroutine set_forcing(sfc_state, forcing, fluxes, day_start, day_interval, G, U

if ((CS%variable_buoyforce .or. CS%first_call_set_forcing) .and. &
(.not.CS%adiabatic)) then
call set_net_mass_forcing(fluxes, forces, G)
call set_net_mass_forcing(fluxes, forces, G, US)
endif

CS%first_call_set_forcing = .false.
Expand Down Expand Up @@ -670,7 +670,7 @@ subroutine buoyancy_forcing_from_files(sfc_state, fluxes, day, dt, G, CS)
temp(:,:), G%Domain, timelevel=time_lev)
do j=js,je ; do i=is,ie
fluxes%latent(i,j) = -hlv*temp(i,j)
fluxes%evap(i,j) = -temp(i,j)
fluxes%evap(i,j) = -US%kg_m3_to_R*US%m_to_Z*US%T_to_s * temp(i,j)
fluxes%latent_evap_diag(i,j) = fluxes%latent(i,j)

enddo ; enddo
Expand All @@ -688,20 +688,20 @@ subroutine buoyancy_forcing_from_files(sfc_state, fluxes, day, dt, G, CS)
enddo ; enddo

call MOM_read_data(trim(CS%inputdir)//trim(CS%snow_file), "snow", &
fluxes%fprec(:,:), G%Domain, timelevel=time_lev)
fluxes%fprec(:,:), G%Domain, timelevel=time_lev, scale=US%kg_m3_to_R*US%m_to_Z*US%T_to_s)
call MOM_read_data(trim(CS%inputdir)//trim(CS%precip_file), "precip", &
fluxes%lprec(:,:), G%Domain, timelevel=time_lev)
fluxes%lprec(:,:), G%Domain, timelevel=time_lev, scale=US%kg_m3_to_R*US%m_to_Z*US%T_to_s)
do j=js,je ; do i=is,ie
fluxes%lprec(i,j) = fluxes%lprec(i,j) - fluxes%fprec(i,j)
enddo ; enddo

call MOM_read_data(trim(CS%inputdir)//trim(CS%freshdischarge_file), "disch_w", &
temp(:,:), G%Domain, timelevel=time_lev_monthly)
temp(:,:), G%Domain, timelevel=time_lev_monthly, scale=US%kg_m3_to_R*US%m_to_Z*US%T_to_s)
do j=js,je ; do i=is,ie
fluxes%lrunoff(i,j) = temp(i,j)*US%m_to_L**2*G%IareaT(i,j)
enddo ; enddo
call MOM_read_data(trim(CS%inputdir)//trim(CS%freshdischarge_file), "disch_s", &
temp(:,:), G%Domain, timelevel=time_lev_monthly)
temp(:,:), G%Domain, timelevel=time_lev_monthly, scale=US%kg_m3_to_R*US%m_to_Z*US%T_to_s)
do j=js,je ; do i=is,ie
fluxes%frunoff(i,j) = temp(i,j)*US%m_to_L**2*G%IareaT(i,j)
enddo ; enddo
Expand Down Expand Up @@ -731,10 +731,11 @@ subroutine buoyancy_forcing_from_files(sfc_state, fluxes, day, dt, G, CS)
fluxes%sw(i,j) = fluxes%sw(i,j) * G%mask2dT(i,j)
fluxes%latent(i,j) = fluxes%latent(i,j) * G%mask2dT(i,j)

fluxes%heat_content_lrunoff(i,j) = fluxes%C_p*fluxes%lrunoff(i,j)*sfc_state%SST(i,j)
fluxes%heat_content_lrunoff(i,j) = fluxes%C_p*US%R_to_kg_m3*US%Z_to_m*US%s_to_T * &
fluxes%lrunoff(i,j)*sfc_state%SST(i,j)
fluxes%latent_evap_diag(i,j) = fluxes%latent_evap_diag(i,j) * G%mask2dT(i,j)
fluxes%latent_fprec_diag(i,j) = -fluxes%fprec(i,j)*hlf
fluxes%latent_frunoff_diag(i,j) = -fluxes%frunoff(i,j)*hlf
fluxes%latent_fprec_diag(i,j) = -US%R_to_kg_m3*US%Z_to_m*US%s_to_T*fluxes%fprec(i,j)*hlf
fluxes%latent_frunoff_diag(i,j) = -US%R_to_kg_m3*US%Z_to_m*US%s_to_T*fluxes%frunoff(i,j)*hlf
enddo ; enddo

endif ! time_lev /= CS%buoy_last_lev_read
Expand All @@ -745,7 +746,7 @@ subroutine buoyancy_forcing_from_files(sfc_state, fluxes, day, dt, G, CS)
if (G%mask2dT(i,j) > 0) then
fluxes%heat_restore(i,j) = G%mask2dT(i,j) * &
((CS%T_Restore(i,j) - sfc_state%SST(i,j)) * rhoXcp * CS%Flux_const)
fluxes%vprec(i,j) = - ((US%R_to_kg_m3*CS%Rho0)*CS%Flux_const) * &
fluxes%vprec(i,j) = - ((US%m_to_Z*US%T_to_s*CS%Rho0)*CS%Flux_const) * &
(CS%S_Restore(i,j) - sfc_state%SSS(i,j)) / &
(0.5*(sfc_state%SSS(i,j) + CS%S_Restore(i,j)))
else
Expand Down Expand Up @@ -878,7 +879,7 @@ subroutine buoyancy_forcing_linear(sfc_state, fluxes, day, dt, G, CS)
if (G%mask2dT(i,j) > 0) then
fluxes%heat_restore(i,j) = G%mask2dT(i,j) * &
((T_Restore - sfc_state%SST(i,j)) * (((US%R_to_kg_m3*CS%Rho0) * fluxes%C_p) * CS%Flux_const))
fluxes%vprec(i,j) = - ((US%R_to_kg_m3*CS%Rho0)*CS%Flux_const) * &
fluxes%vprec(i,j) = - ((US%m_to_Z*US%T_to_s*CS%Rho0)*CS%Flux_const) * &
(S_Restore - sfc_state%SSS(i,j)) / &
(0.5*(sfc_state%SSS(i,j) + S_Restore))
else
Expand Down
27 changes: 13 additions & 14 deletions config_src/ice_solo_driver/user_surface_forcing.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -79,11 +79,11 @@ module user_surface_forcing
! state variables.
logical :: restorebuoy ! If true, use restoring surface buoyancy forcing.
real :: Rho0 ! The density used in the Boussinesq
! approximation [kg m-3].
! approximation [R ~> kg m-3].
real :: G_Earth ! The gravitational acceleration [L2 Z-1 T-2 ~> m s-2].
real :: Flux_const ! The restoring rate at the surface [m s-1].
real :: gust_const ! A constant unresolved background gustiness
! that contributes to ustar [Pa].
! that contributes to ustar [R Z L T-1 ~> Pa].

type(diag_ctrl), pointer :: diag ! A structure that is used to regulate the
! timing of diagnostic output.
Expand Down Expand Up @@ -140,9 +140,9 @@ subroutine USER_wind_forcing(sfc_state, forces, day, G, US, CS)
! Set the surface friction velocity [Z s-1 ~> m s-1]. ustar is always positive.
if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie
! This expression can be changed if desired, but need not be.
forces%ustar(i,j) = US%m_to_Z*US%T_to_s * G%mask2dT(i,j) * sqrt(CS%gust_const/CS%Rho0 + &
US%R_to_kg_m3*US%L_T_to_m_s**2*US%Z_to_L*sqrt(0.5*(forces%taux(I-1,j)**2 + forces%taux(I,j)**2) + &
0.5*(forces%tauy(i,J-1)**2 + forces%tauy(i,J)**2))/CS%Rho0)
forces%ustar(i,j) = G%mask2dT(i,j) * sqrt(US%L_to_Z * (CS%gust_const/CS%Rho0 + &
sqrt(0.5*(forces%taux(I-1,j)**2 + forces%taux(I,j)**2) + &
0.5*(forces%tauy(i,J-1)**2 + forces%tauy(i,J)**2))/CS%Rho0))
enddo ; enddo ; endif

end subroutine USER_wind_forcing
Expand Down Expand Up @@ -174,7 +174,7 @@ subroutine USER_buoyancy_forcing(sfc_state, fluxes, day, dt, G, US, CS)
! (fprec, lrunoff and frunoff) left as arrays full of zeros.
! Evap is usually negative and precip is usually positive. All heat fluxes
! are in W m-2 and positive for heat going into the ocean. All fresh water
! fluxes are in kg m-2 s-1 and positive for water moving into the ocean.
! fluxes are in [R Z T-1 ~> kg m-2 s-1] and positive for water moving into the ocean.

real :: Temp_restore ! The temperature that is being restored toward [C].
real :: Salin_restore ! The salinity that is being restored toward [ppt]
Expand Down Expand Up @@ -250,7 +250,7 @@ subroutine USER_buoyancy_forcing(sfc_state, fluxes, day, dt, G, US, CS)
call MOM_error(FATAL, "User_buoyancy_surface_forcing: " // &
"Temperature and salinity restoring used without modification." )

rhoXcp = CS%Rho0 * fluxes%C_p
rhoXcp = US%R_to_kg_m3*CS%Rho0 * fluxes%C_p
do j=js,je ; do i=is,ie
! Set Temp_restore and Salin_restore to the temperature (in degC) and
! salinity (in ppt or PSU) that are being restored toward.
Expand All @@ -259,9 +259,8 @@ subroutine USER_buoyancy_forcing(sfc_state, fluxes, day, dt, G, US, CS)

fluxes%heat_added(i,j) = (G%mask2dT(i,j) * (rhoXcp * CS%Flux_const)) * &
(Temp_restore - sfc_state%SST(i,j))
fluxes%vprec(i,j) = - (G%mask2dT(i,j) * (CS%Rho0*CS%Flux_const)) * &
((Salin_restore - sfc_state%SSS(i,j)) / &
(0.5 * (Salin_restore + sfc_state%SSS(i,j))))
fluxes%vprec(i,j) = - (G%mask2dT(i,j) * (US%m_to_Z*US%T_to_s*CS%Rho0*CS%Flux_const)) * &
((Salin_restore - sfc_state%SSS(i,j)) / (0.5 * (Salin_restore + sfc_state%SSS(i,j))))
enddo ; enddo
else
! When modifying the code, comment out this error message. It is here
Expand All @@ -270,7 +269,7 @@ subroutine USER_buoyancy_forcing(sfc_state, fluxes, day, dt, G, US, CS)
"Buoyancy restoring used without modification." )

! The -1 is because density has the opposite sign to buoyancy.
buoy_rest_const = -1.0 * (CS%G_Earth * US%m_to_Z*US%T_to_s*CS%Flux_const) / CS%Rho0
buoy_rest_const = -1.0 * (CS%G_Earth * US%m_to_Z*US%T_to_s*CS%Flux_const) / (US%R_to_kg_m3*CS%Rho0)
do j=js,je ; do i=is,ie
! Set density_restore to an expression for the surface potential
! density [kg m-3] that is being restored toward.
Expand Down Expand Up @@ -320,10 +319,10 @@ subroutine USER_surface_forcing_init(Time, G, US, param_file, diag, CS)
"calculate accelerations and the mass for conservation "//&
"properties, or with BOUSSINSEQ false to convert some "//&
"parameters from vertical units of m to kg m-2.", &
units="kg m-3", default=1035.0)
units="kg m-3", default=1035.0, scale=US%R_to_kg_m3)
call get_param(param_file, mdl, "GUST_CONST", CS%gust_const, &
"The background gustiness in the winds.", units="Pa", &
default=0.02)
"The background gustiness in the winds.", &
units="Pa", default=0.02, scale=US%kg_m3_to_R*US%m_s_to_L_T**2*US%L_to_Z)

call get_param(param_file, mdl, "RESTOREBUOY", CS%restorebuoy, &
"If true, the buoyancy fluxes drive the model back "//&
Expand Down
Loading

0 comments on commit dbad998

Please sign in to comment.