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.

Sign up for GitHub

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

Jump to bottom

Add CAM diagnostic schemes #105

Merged
merged 13 commits into from
Oct 17, 2024
Merged

Add CAM diagnostic schemes #105

Show file tree
Hide file tree
Changes from 6 commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
cc42526
add cam diagnostics scheme(s)
Jul 11, 2024
ca7794b
add scheme diagnostics template
Jul 24, 2024
78b8bfb
use parsable temp name
Jul 24, 2024
6363a53
consistent error message length
Jul 29, 2024
61bd9d3
Merge remote-tracking branch 'ESCOMP/development' into diagnostics
peverwhee Aug 23, 2024
8d44484
address review comments
Sep 16, 2024
9b4a79a
fix standard names/ic file names
Sep 25, 2024
06dee83
split diagnostics
Sep 27, 2024
e89e912
rename cam_diagnostics to sima_diagnostics
Oct 2, 2024
ff29507
merge to head of development
Oct 2, 2024
a3ee11b
fix cam7 sdf
Oct 10, 2024
9f8fc41
merge to head of develop; move sima_diagnostics to schemes directory
Oct 10, 2024
17d4465
update NamesNotInDictionary to reflect new directory structure
Oct 10, 2024
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
211 changes: 211 additions & 0 deletions cam_diagnostics/cam_diagnostics.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,211 @@
module cam_diagnostics

use ccpp_kinds, only: kind_phys
use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t

implicit none
private
save

public :: cam_state_diagnostics_init ! init routine
public :: cam_state_diagnostics_run ! main routine
public :: cam_tend_diagnostics_init ! init routine
public :: cam_tend_diagnostics_run ! main routine
cacraigucar marked this conversation as resolved.
Show resolved Hide resolved

character(len=65) :: const_std_names(4) = &
(/'water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water ', &
'cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water', &
'rain_mixing_ratio_wrt_moist_air_and_condensed_water ', &
'cloud_ice_mixing_ratio_wrt_moist_air_and_condensed_water '/)

character(len=6) :: const_diag_names(4) = (/'Q ', &
'CLDLIQ', &
'CLDICE', &
'RAINQM'/)
Copy link
Collaborator

Choose a reason for hiding this comment

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

Shouldn't the contents of these two variables be dictated by whatever is populating the constituent object and loop over what is in that object?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

we're looping over the constituents object, yes. But this is a mapping of the standard name to the diagnostic name (which is not contained within the object). Not all of these will necessarily be "used". Does that make any sense?

Copy link
Collaborator

Choose a reason for hiding this comment

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

This still feels "fragile" to me, as the diagnostic name index may not match the constituent index. Right now in ZM, I am assuming that the tendency array maps one-to-one to the constituent array. I would naively think that the diagnostic array would also use the same indexing. Let me know if we should have a discussion on this.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Hi @cacraigucar I'm not quite following. Is the problem that someone could modify one of the lists (const_std_names or const_diag_names) and then they'd be off? To clarify, when I'm looping over the constituents, I don't expect them to be in this exact order. This is just a mapping between standard name and diagnostic name (which is not contained in the properties object). Let me know if I'm still missing something!

Copy link
Member

Choose a reason for hiding this comment

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

Hi @peverwhee, @cacraigucar. If I understand this fragility correctly it might be possible that user error could cause the lists to become out-of-sync. I wonder of ways to mitigate this, e.g., checking the lists are of same length at initialization (to prevent the obvious typo of adding one but not the other), but all bets are off if an entry is added to the wrong place in one of the lists.

My immediate thought was to assign these into a Nx2 array instead of two N-sized arrays but a quick Google didn't reveal if Fortran had syntax like this, which would prevent typing errors as much as possible:

const_std_and_diag_names(4,2) = (/ /'water_vapor_mixing_ratio..', 'Q'/, ...
                                   /'cloud_liquid_water...', 'CLDLIQ'/, ...
                                )

Copy link
Contributor

Choose a reason for hiding this comment

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

When I do this, I always use a parameter for each dimension size.
The only example of this I found in CAM code is from the very ancient (well before my time) src/physics/camrt/radae.F90

Newer code, such as src/dynamics/se/dycore/fvm_consistent_se_cslam.F90 does have a parameter for the changeable dimension which I think addresses your main safety concern.

Copy link
Member

Choose a reason for hiding this comment

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

Thanks! I like this in fvm_consistent_se_cslam.F90, which I guess is what you're referring to: it's as close to safe and 'automatic' as we can get before the clean-up issue #137 is addressed:

    integer, parameter :: num_area=5
    integer, dimension(num_area*4), parameter :: idy_shift_tmp = (/-1, 0, 0, 0,-1,&  !iside=1
                                                                   -1,-1, 0, 1, 1,&  !iside=2
                                                                    1, 0, 0, 0, 1,&  !iside=3
                                                                    1, 1, 0,-1,-1/)  !iside=4

    integer, dimension(num_area,4), parameter :: idy_shift = RESHAPE(idy_shift_tmp,(/num_area,4/))

Copy link
Collaborator

Choose a reason for hiding this comment

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

While this conversation has been good, my concern was not actually about there being two arrays, but rather about the necessity of having them at all. I had a google chat with Courtney which resulted in #137 and will eventually eliminate the need for these arrays.

Copy link
Member

Choose a reason for hiding this comment

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

Thanks @cacraigucar, my apologies for misunderstanding the original issue. Until #137 becomes a more permanent fix the current code looks good!

Copy link
Collaborator

Choose a reason for hiding this comment

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

No apology needed. Nobody understood what my concern was until I met with Courtney and explained it in person!


CONTAINS

!> \section arg_table_cam_state_diagnostics_init Argument Table
!! \htmlinclude cam_state_diagnostics_init.html
subroutine cam_state_diagnostics_init(const_props, errmsg, errflg)
use cam_history, only: history_add_field
use cam_history_support, only: horiz_only

type(ccpp_constituent_prop_ptr_t), intent(in) :: const_props(:)
character(len=512), intent(out) :: errmsg
integer, intent(out) :: errflg

! Local variables:

integer :: const_idx, name_idx
integer :: const_num_found
character(len=512) :: standard_name

errmsg = ''
errflg = 0

! Add state fields
cacraigucar marked this conversation as resolved.
Show resolved Hide resolved
call history_add_field('PS', 'surface_pressure', horiz_only, 'avg', 'Pa')
call history_add_field('PSDRY', 'surface_pressure_of_dry_air', horiz_only, 'avg', 'Pa')
call history_add_field('PHIS', 'surface_geopotential', horiz_only, 'lst', 'Pa')
call history_add_field('T', 'air_temperature', 'lev', 'avg', 'K')
call history_add_field('U', 'eastward_wind', 'lev', 'avg', 'm s-1')
call history_add_field('V', 'northward_wind', 'lev', 'avg', 'm s-1')
call history_add_field('DSE', 'dry_static_energy', 'lev', 'avg', 'm s-1')
call history_add_field('OMEGA', 'lagrangian_tendency_of_air_pressure', 'lev', 'avg', 'Pa s-1')
call history_add_field('PMID', 'air_pressure', 'lev', 'avg', 'Pa')
call history_add_field('PMIDDRY', 'air_pressure_of_dry_air', 'lev', 'avg', 'Pa')
call history_add_field('PDEL', 'air_pressure_thickness', 'lev', 'avg', 'Pa')
call history_add_field('PDELDRY', 'air_pressure_thickness_of_dry_air', 'lev', 'avg', 'Pa')
call history_add_field('RPDEL', 'reciprocal_of_air_pressure_thickness', 'lev', 'avg', 'Pa-1')
call history_add_field('RPDELDRY', 'reciprocal_of_air_pressure_thickness_of_dry_air', 'lev', 'avg', 'Pa-1')
call history_add_field('LNPMID', 'ln_air_pressure', 'lev', 'avg', '1')
call history_add_field('LNPMIDDRY', 'ln_air_pressure_of_dry_air', 'lev', 'avg', '1')
call history_add_field('EXNER', 'reciprocal_of_dimensionless_exner_function_wrt_surface_air_pressure','lev', 'avg', '1')
call history_add_field('ZM', 'geopotential_height_wrt_surface', 'lev', 'avg', 'm')
call history_add_field('PINT', 'air_pressure_at_interfaces', 'ilev', 'avg', 'Pa')
call history_add_field('PINTDRY', 'air_pressure_of_dry_air_at_interfaces', 'ilev', 'avg', 'Pa')
call history_add_field('LNPINT', 'ln_air_pressure_at_interfaces', 'ilev', 'avg', '1')
call history_add_field('LNPINTDRY', 'ln_air_pressure_of_dry_air_at_interfaces', 'ilev', 'avg', '1')
call history_add_field('ZI', 'geopotential_height_wrt_surface_at_interfaces', 'ilev', 'avg', 'm')
! Add constituent fields
const_num_found = 0
do const_idx = 1, size(const_props)
call const_props(const_idx)%standard_name(standard_name, errflg, errmsg)
do name_idx = 1, size(const_std_names)
if (trim(standard_name) == trim(const_std_names(name_idx))) then
call history_add_field(trim(const_diag_names(name_idx)), trim(const_std_names(name_idx)), 'lev', 'avg', 'kg kg-1', mixing_ratio='wet')
const_num_found = const_num_found + 1
end if
end do
if (const_num_found == size(const_std_names)) then
exit
end if
end do

end subroutine cam_state_diagnostics_init

!> \section arg_table_cam_state_diagnostics_run Argument Table
!! \htmlinclude cam_state_diagnostics_run.html
subroutine cam_state_diagnostics_run(ps, psdry, phis, T, u, v, dse, omega, &
pmid, pmiddry, pdel, pdeldry, rpdel, rpdeldry, lnpmid, lnpmiddry, &
inv_exner, zm, pint, pintdry, lnpint, lnpintdry, zi, const_array, &
const_props, errmsg, errflg)

use cam_history, only: history_out_field
!------------------------------------------------
! Input / output parameters
!------------------------------------------------
! State variables
real(kind_phys), intent(in) :: ps(:) ! surface pressure
real(kind_phys), intent(in) :: psdry(:) ! surface pressure of dry air
real(kind_phys), intent(in) :: phis(:) ! surface geopotential
real(kind_phys), intent(in) :: T(:,:) ! air temperature
real(kind_phys), intent(in) :: u(:,:) ! eastward wind (x wind)
real(kind_phys), intent(in) :: v(:,:) ! northward wind (y wind)
real(kind_phys), intent(in) :: dse(:,:) ! dry static energy
real(kind_phys), intent(in) :: omega(:,:) ! lagrangian tendency of air pressure
real(kind_phys), intent(in) :: pmid(:,:) ! air pressure
real(kind_phys), intent(in) :: pmiddry(:,:) ! air pressure of dry air
real(kind_phys), intent(in) :: pdel(:,:) ! air pressure thickness
real(kind_phys), intent(in) :: pdeldry(:,:) ! air pressure thickness of dry air
real(kind_phys), intent(in) :: rpdel(:,:) ! reciprocal of air pressure thickness
real(kind_phys), intent(in) :: rpdeldry(:,:) ! reciprocal of air pressure thickness of dry air
real(kind_phys), intent(in) :: lnpmid(:,:) ! ln air pressure
real(kind_phys), intent(in) :: lnpmiddry(:,:) ! ln air pressure of dry air
real(kind_phys), intent(in) :: inv_exner(:,:) ! inverse exner function wrt surface pressure
real(kind_phys), intent(in) :: zm(:,:) ! geopotential height wrt surface
real(kind_phys), intent(in) :: pint(:,:) ! air pressure at interfaces
real(kind_phys), intent(in) :: pintdry(:,:) ! air pressure of dry air at interfaces
real(kind_phys), intent(in) :: lnpint(:,:) ! ln air pressure at interfaces
real(kind_phys), intent(in) :: lnpintdry(:,:) ! ln air pressure of dry air at interfaces
real(kind_phys), intent(in) :: zi(:,:) ! geopotential height wrt surface at interfaces
! Constituent variables
real(kind_phys), intent(in) :: const_array(:,:,:)
type(ccpp_constituent_prop_ptr_t), intent(in) :: const_props(:)
! CCPP error handling variables
character(len=512), intent(out) :: errmsg
integer, intent(out) :: errflg

character(len=512) :: standard_name
integer :: const_idx, name_idx
integer :: const_num_found

errmsg = ''
errflg = 0

! Capture state fields
call history_out_field('PS' , ps)
call history_out_field('PSDRY' , psdry)
call history_out_field('PHIS' , phis)
call history_out_field('T' , T)
call history_out_field('U' , u)
call history_out_field('V' , v)
call history_out_field('DSE' , dse)
call history_out_field('OMEGA' , omega)
call history_out_field('PMID' , pmid)
call history_out_field('PMIDDRY' , pmiddry)
call history_out_field('PDEL' , pdel)
call history_out_field('PDELDRY' , pdeldry)
call history_out_field('RPDEL' , rpdel)
call history_out_field('RPDELDRY' , rpdeldry)
call history_out_field('LNPMID' , lnpmid)
call history_out_field('LNPMIDDRY', lnpmiddry)
call history_out_field('EXNER' , inv_exner)
call history_out_field('ZM' , zm)
call history_out_field('PINT' , pint)
call history_out_field('PINTDRY' , pintdry)
call history_out_field('LNPINT' , lnpint)
call history_out_field('LNPINTDRY', lnpintdry)
call history_out_field('ZI' , zi)

! Capture constituent fields
const_num_found = 0
do const_idx = 1, size(const_props)
call const_props(const_idx)%standard_name(standard_name, errflg, errmsg)
do name_idx = 1, size(const_std_names)
if (trim(standard_name) == trim(const_std_names(name_idx))) then
call history_out_field(trim(const_diag_names(name_idx)), const_array(:,:,const_idx))
const_num_found = const_num_found + 1
end if
end do
if (const_num_found == size(const_std_names)) then
exit
end if
end do

end subroutine cam_state_diagnostics_run

!> \section arg_table_cam_tend_diagnostics_init Argument Table
!! \htmlinclude cam_tend_diagnostics_init.html
subroutine cam_tend_diagnostics_init(errmsg, errflg)
use cam_history, only: history_add_field
character(len=512), intent(out) :: errmsg
integer, intent(out) :: errflg

! Add tendency fields
call history_add_field('TTEND', 'tendency_of_air_temperature_due_to_model_physics', 'lev', 'avg', 'K s-1')
call history_add_field('UTEND', 'tendency_of_eastward_wind_due_to_model_physics', 'lev', 'avg', 'm s-2')
call history_add_field('VTEND', 'tendency_of_northward_wind_due_to_model_physics', 'lev', 'avg', 'm s-2')
nusbaume marked this conversation as resolved.
Show resolved Hide resolved

end subroutine cam_tend_diagnostics_init

!> \section arg_table_cam_tend_diagnostics_run Argument Table
!! \htmlinclude cam_tend_diagnostics_run.html
subroutine cam_tend_diagnostics_run(dTdt_total, dudt_total, dvdt_total, errmsg, errflg)
use cam_history, only: history_out_field
! Tendency variables
real(kind_phys), intent(in) :: dTdt_total(:,:) ! tendency of air temperature due to model physics
real(kind_phys), intent(in) :: dudt_total(:,:) ! tendency of eastward wind due to model physics
real(kind_phys), intent(in) :: dvdt_total(:,:) ! tendency of northward wind due to model physics
cacraigucar marked this conversation as resolved.
Show resolved Hide resolved
character(len=512), intent(out) :: errmsg
integer, intent(out) :: errflg

! Capture tendency fields
call history_out_field('TTEND', dTdt_total)
call history_out_field('UTEND', dudt_total)
call history_out_field('VTEND', dvdt_total)

end subroutine cam_tend_diagnostics_run
!=======================================================================
end module cam_diagnostics
Loading