Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Added diagnostics for partial derivatives of density #1142

Merged
merged 6 commits into from
Jul 3, 2020
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
23 changes: 22 additions & 1 deletion src/diagnostics/MOM_diagnostics.F90
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ module MOM_diagnostics
use MOM_diag_mediator, only : diag_save_grids, diag_restore_grids, diag_copy_storage_to_diag
use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type
use MOM_domains, only : To_North, To_East
use MOM_EOS, only : calculate_density, int_density_dz, EOS_domain
use MOM_EOS, only : calculate_density, calculate_density_derivs, int_density_dz, EOS_domain
use MOM_EOS, only : gsw_sp_from_sr, gsw_pt_from_ct
use MOM_error_handler, only : MOM_error, FATAL, WARNING
use MOM_file_parser, only : get_param, log_version, param_file_type
Expand Down Expand Up @@ -134,6 +134,7 @@ module MOM_diagnostics
integer :: id_pbo = -1
integer :: id_thkcello = -1, id_rhoinsitu = -1
integer :: id_rhopot0 = -1, id_rhopot2 = -1
integer :: id_drho_dT = -1, id_drho_dS = -1
integer :: id_h_pre_sync = -1
!>@}
!> The control structure for calculating wave speed.
Expand Down Expand Up @@ -619,6 +620,22 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
enddo
if (CS%id_rhoinsitu > 0) call post_data(CS%id_rhoinsitu, Rcv, CS%diag)
endif

if (CS%id_drho_dT > 0 .or. CS%id_drho_dS > 0) then
!$OMP parallel do default(shared) private(pressure_1d)
do j=js,je
pressure_1d(:) = 0. ! Start at p=0 Pa at surface
do k=1,nz
pressure_1d(:) = pressure_1d(:) + 0.5 * h(:,j,k) * GV%H_to_Pa ! Pressure in middle of layer k
! To avoid storing more arrays, put drho_dT into Rcv, and drho_dS into work3d
call calculate_density_derivs(tv%T(:,j,k),tv%S(:,j,k),pressure_1d, &
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You need to import this function in module specification (L18-19).

(Was this able to build for you?)

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Aha, this was my error. I forgot to push my latest commit after making that change 🤦

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

All good 🎉
BTW do you think you are able to sort out the OpenMP error as well? If not, then I can look into it.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think I have the knowledge to sort that. I'd be happy if you would talk me through it though.

Is it failing because of something in the changes that I've made? Or have I done something wrong in not merging my branch with dev/gfdl correctly in the first place?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The OpenMP directive preceding the block (line starting with !$OMP) has an error because you have default(none) and haven't specified if work_3d is shared across threads or not.

The directive could be fixed by adding work_3d to the OpenMP directive. But if you did not mean to add it, then it might be simpler to just remove it.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not knowing how OpenMP functions, I don't have a sense for whether it is necessary to have work_3d shared across threads. I was just copying the convention of the other density calculations in including this line.

I think, in fact, that the line as present in my PR was copied from a previous version of MOM_diagnostic.F90, as I see that the OpenMP directives for the other density variables (e.g. rhoinsitu) have changed. I wonder if this line should instead just be consistent with that for rhoinsitu? i.e. !$OMP parallel do default(shared) private(pressure_1d)

I will implement that just now and you can tell me if I'm missing something.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That looks good, and also looks like it passed. (Just waiting on the ARM tests...)
The regression fail is expected since this is a new diagnostic. I'll push this to Gaea and if it passes then I will merge it.

Thanks for fixing that up!

Rcv(:,j,k),work_3d(:,j,k),is,ie-is+1, tv%eqn_of_state)
pressure_1d(:) = pressure_1d(:) + 0.5 * h(:,j,k) * GV%H_to_Pa ! Pressure at bottom of layer k
enddo
enddo
if (CS%id_drho_dT > 0) call post_data(CS%id_drho_dT, Rcv, CS%diag)
if (CS%id_drho_dS > 0) call post_data(CS%id_drho_dS, work_3d, CS%diag)
endif
endif

if ((CS%id_cg1>0) .or. (CS%id_Rd1>0) .or. (CS%id_cfl_cg1>0) .or. &
Expand Down Expand Up @@ -1600,6 +1617,10 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag
'Potential density referenced to 2000 dbar', 'kg m-3', conversion=US%R_to_kg_m3)
CS%id_rhoinsitu = register_diag_field('ocean_model', 'rhoinsitu', diag%axesTL, Time, &
'In situ density', 'kg m-3', conversion=US%R_to_kg_m3)
CS%id_drho_dT = register_diag_field('ocean_model', 'drho_dT', diag%axesTL, Time, &
'Partial derivative of rhoinsitu with respect to temperature (alpha)', 'kg m-3 degC-1')
CS%id_drho_dS = register_diag_field('ocean_model', 'drho_dS', diag%axesTL, Time, &
'Partial derivative of rhoinsitu with respect to salinity (beta)', 'kg^2 g-1 m-3')

CS%id_du_dt = register_diag_field('ocean_model', 'dudt', diag%axesCuL, Time, &
'Zonal Acceleration', 'm s-2', conversion=US%L_T2_to_m_s2)
Expand Down