Skip to content

Commit

Permalink
Simplified the vertical pressure conversion interfaces so that only t…
Browse files Browse the repository at this point in the history 
…he column interface

does actual computation. This version includes preliminary code for adding support
for dry mass vertical coordinate but this is incomplete.
  • Loading branch information
@jlaucar
jlaucar committed Jun 29, 2020
1 parent c40bf6a commit 1982c75
Show file tree
Hide file tree
Showing 2 changed files with 107 additions and 102 deletions.
94 changes: 93 additions & 1 deletion models/CAM-SE-SHARED/model_mod.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -290,6 +290,12 @@ module model_mod
logical :: print_details = .true.


!SENote TEMPORARY LOGICAL THAT SHOULD GO IN MODEL NAMELIST FOR NOW?
logical :: DRY_MASS_VERTICAL_COORDINATE = .false.





contains

Expand Down Expand Up @@ -2100,7 +2106,15 @@ subroutine find_se_vertical_levels(ens_handle, ens_size, column, vert_val, &
return
endif

call build_cam_pressure_columns(ens_size, surf_pressure, ref_nlevels, pressure_array)
if(DRY_MASS_VERTICAL_COORDINATE) then
call build_dry_mass_pressure_columns(ens_handle, ens_size, ref_nlevels, column, surf_pressure, &
pressure_array, status1)
else
call build_cam_pressure_columns(ens_size, surf_pressure, ref_nlevels, pressure_array)
endif




do imember = 1, ens_size
call pressure_to_level(ref_nlevels, pressure_array(:, imember), vert_val, &
Expand Down Expand Up @@ -2190,6 +2204,84 @@ subroutine find_se_vertical_levels(ens_handle, ens_size, column, vert_val, &
end subroutine find_se_vertical_levels


!-----------------------------------------------------------------------
!> Compute pressure column for the dry mass vertical coordinate option
!>
!> this version does all ensemble members at once.


subroutine build_dry_mass_pressure_columns(ens_handle, ens_size, nlevels, column, surf_pressure, pressure, status)
type(ensemble_type), intent(in) :: ens_handle
integer, intent(in) :: ens_size
integer, intent(in) :: nlevels
integer, intent(in) :: column
real(r8), intent(in) :: surf_pressure(ens_size)
real(r8), intent(out) :: pressure(nlevels, ens_size)
integer, intent(out) :: status

real(r8), dimension(ens_size) :: specific_humidity, cldliq, cldice, sum_specific_water_ratios
real(r8) :: sum_dry_mix_ratio(nlevels, ens_size)
integer :: k, istatus

! Building pressure columns for dry mass vertical coordinate; Add in references to Lauritzen and pointer
! to the document on the algorithm

! Begin by getting column values for specific mixing ratio for water vapor (specific humidity), cloud liquid
! and cloud ice. For more accuracy could also include rain, snow, and any other tracers.
! Note that other tracers in the dry mass cam/se have a dry mixing ratio, not specific mixing ratio (although
! Peter plans to change this for consistency).

!SENote: The following should all be 0, but A's at levels 2, 3, 4, and 5 are not.
! This is confirmed in the caminput.nc file and is a PROBLEM.
! Some tests on the A and B coefficients
!do k = 1, nlevels
!write(*, *) k, (grid_data%hyai%vals(k) + grid_data%hyai%vals(k+1))/2.0_r8 - grid_data%hyam%vals(k)
!write(*, *) k, (grid_data%hybi%vals(k) + grid_data%hybi%vals(k+1))/2.0_r8 - grid_data%hybm%vals(k)
!enddo
!stop
! Another test
write(*, *) sum(grid_data%hyam%vals(:)), grid_data%hyai%vals(33) - grid_data%hyai%vals(1)
stop



pressure = MISSING_R8
status = 1

! Need the water tracer specific mixing ratios for ever level in the column to compute their mass sum
do k = 1, nlevels
! Specific Humidity
call get_se_values_from_single_level(ens_handle, ens_size, QTY_SPECIFIC_HUMIDITY, column, k, &
specific_humidity, istatus)

if (istatus < 0) return

! Cloud liquid
call get_se_values_from_single_level(ens_handle, ens_size, QTY_CLOUD_LIQUID_WATER, column, k, &
cldliq, istatus)

if (istatus < 0) return


! Cloud ice
call get_se_values_from_single_level(ens_handle, ens_size, QTY_CLOUD_ICE, column, k, &
cldice, istatus)

if (istatus < 0) return

! Compute the sum of the dry mixing ratio of dry air plus all the water tracers (ref. to notes)
sum_specific_water_ratios = specific_humidity + cldliq + cldice
sum_dry_mix_ratio(k, :) = 1.0_r8 + sum_specific_water_ratios / (1.0_r8 - sum_specific_water_ratios)

enddo

! Compute the dry mass at the bottom of the column



end subroutine build_dry_mass_pressure_columns


!-----------------------------------------------------------------------
!> Compute the heights at pressure midpoints
!>
Expand Down
115 changes: 14 additions & 101 deletions models/cam-common-code/cam_common_code_mod.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -115,14 +115,6 @@ module cam_common_code_mod
real(r8), allocatable :: std_atm_pres_col(:)


interface single_pressure_value
module procedure single_pressure_value_int
module procedure single_pressure_value_real
end interface




contains


Expand Down Expand Up @@ -693,35 +685,6 @@ end function generic_height_to_pressure


!--------------------------------------------------------------------
!SENote: This does not appear to be used at this point.
!SENote: This means that we don't need to use P0 at all?

!> using the cam eta arrays, convert a pressure directly to model level
!> use P0 as surface, ignore elevation.
!> This is not currently being used but might be useful.

function generic_cam_pressure_to_cam_level(pressure, status)
real(r8), intent(in) :: pressure
integer, intent(out) :: status
real(r8) :: generic_cam_pressure_to_cam_level

integer :: lev1, lev2
real(r8) :: fract
real(r8) :: pressure_array(ref_nlevels)

generic_cam_pressure_to_cam_level = MISSING_R8

call single_pressure_column(ref_surface_pressure, ref_nlevels, pressure_array)

call pressure_to_level(ref_nlevels, pressure_array, pressure, &
lev1, lev2, fract, status)
if (status /= 0) return

generic_cam_pressure_to_cam_level = lev1 + fract

end function generic_cam_pressure_to_cam_level

!-----------------------------------------------------------------------
!> in cam level 1 is at the model top, level N is the lowest level
!> our convention in this code is: between levels a fraction of 0
!> is 100% level 1, and fraction of 1 is 100% level 2.
Expand Down Expand Up @@ -836,12 +799,9 @@ subroutine single_pressure_column(surface_pressure, n_levels, pressure_array)

integer :: k

! Set midpoint pressures. This array mirrors the order of the
! cam model levels: 1 is the model top, N is the bottom.

do k=1, n_levels
pressure_array(k) = single_pressure_value_int(surface_pressure, k)
enddo
! Set midpoint pressures.
pressure_array(1:n_levels) = ref_surface_pressure * grid_data%hyam%vals(1:n_levels) + &
surface_pressure * grid_data%hybm%vals(1:n_levels)

end subroutine single_pressure_column

Expand All @@ -863,6 +823,7 @@ subroutine init_discard_high_obs(no_obs_assim_above_level)
integer :: table_type
character(len=16) :: out_fmt, out_fmt1, pres_fmt
real(r8) :: no_assim_above_scaleh
real(r8) :: temp_p_col(ref_nlevels)

! pick the better table:
! one is more accurate for the lower atmosphere,
Expand All @@ -882,7 +843,9 @@ subroutine init_discard_high_obs(no_obs_assim_above_level)
'Discarding observations higher than model level ', no_obs_assim_above_level
call error_handler(E_MSG, 'init_discard_high_obs', string1, source, revision, revdate)

no_assim_above_pressure = single_pressure_value(ref_surface_pressure, no_obs_assim_above_level)
!SENote: Getting rid of single pressure call, converting to only column
call single_pressure_column(ref_surface_pressure, ref_nlevels, temp_p_col)
no_assim_above_pressure = temp_p_col(no_obs_assim_above_level)
write(string1, pres_fmt) &
' ... which is equivalent to pressure level ', no_assim_above_pressure, ' Pascals'
call error_handler(E_MSG, 'init_discard_high_obs', string1, source, revision, revdate)
Expand Down Expand Up @@ -939,7 +902,7 @@ subroutine init_damping_ramp_info(model_damping_ends_at_level)
vertical_localization_type == VERTISPRESSURE) out_fmt = '(A,E12.5,A)'

! convert to vertical localization units
call convert_vertical_level_generic(real(model_damping_ends_at_level, r8), &
call convert_vertical_level_generic(model_damping_ends_at_level, &
vertical_localization_type, ramp_end, string3, no_norm=.false.)

! check for conversion errors
Expand All @@ -950,7 +913,7 @@ subroutine init_damping_ramp_info(model_damping_ends_at_level)
endif

! this value only used for print statement, unused otherwise
call convert_vertical_level_generic(1.0_r8, vertical_localization_type, &
call convert_vertical_level_generic(1, vertical_localization_type, &
model_top, string3, no_norm=.false.)

! check for conversion errors
Expand Down Expand Up @@ -1011,7 +974,7 @@ end subroutine init_sign_of_vert_units
!> it uses generic values to do a vertical conversion.

subroutine convert_vertical_level_generic(level_value, want_vert_type, out_value, out_label, no_norm)
real(r8), intent(in) :: level_value
integer, intent(in) :: level_value
integer, intent(in) :: want_vert_type
real(r8), intent(out) :: out_value
character(len=*), intent(out), optional :: out_label
Expand All @@ -1020,7 +983,7 @@ subroutine convert_vertical_level_generic(level_value, want_vert_type, out_value
character(len=*), parameter :: routine = 'convert_vertical_level_generic'

integer :: status
real(r8) :: tmp_val
real(r8) :: tmp_val, temp_p_col(ref_nlevels)
logical :: no_norm_flag

if (present(no_norm)) then
Expand All @@ -1030,11 +993,12 @@ subroutine convert_vertical_level_generic(level_value, want_vert_type, out_value
endif

if (want_vert_type == VERTISLEVEL) then
out_value = level_value
out_value = real(level_value, r8)
if (present(out_label)) out_label = 'levels'
else
! convert to the requested units. start by going to pressure
tmp_val = single_pressure_value(ref_surface_pressure, level_value)
call single_pressure_column(ref_surface_pressure, ref_nlevels, temp_p_col)
tmp_val = temp_p_col(level_value)

select case (want_vert_type)
case (VERTISPRESSURE)
Expand All @@ -1059,57 +1023,6 @@ subroutine convert_vertical_level_generic(level_value, want_vert_type, out_value

end subroutine convert_vertical_level_generic

!-----------------------------------------------------------------------
!> Compute pressure at one level given the surface pressure
!> cam model levels: 1 is the model top, N is the bottom.
!> in this version of the routine level is integer/whole value

function single_pressure_value_int(surface_pressure, level)

real(r8), intent(in) :: surface_pressure ! in pascals
integer, intent(in) :: level
real(r8) :: single_pressure_value_int

! cam model levels: 1 is the model top, N is the bottom.

single_pressure_value_int = ref_surface_pressure * grid_data%hyam%vals(level) + &
surface_pressure * grid_data%hybm%vals(level)

end function single_pressure_value_int

!-----------------------------------------------------------------------
!> Compute pressure at one level given the surface pressure
!> cam model levels: 1 is the model top, N is the bottom.
!> fraction = 0 is full level 1, fraction = 1 is full level 2
!> level is real/fractional value


function single_pressure_value_real(surface_pressure, level)

real(r8), intent(in) :: surface_pressure ! in pascals
real(r8), intent(in) :: level
real(r8) :: single_pressure_value_real

integer :: k
real(r8) :: fract, pres1, pres2

k = int(level)
fract = level - int(level)

if (k /= ref_nlevels) then
pres1 = single_pressure_value_int(surface_pressure, k)
pres2 = single_pressure_value_int(surface_pressure, k+1)
else
pres1 = single_pressure_value_int(surface_pressure, k-1)
pres2 = single_pressure_value_int(surface_pressure, k)
fract = 1.0_r8
endif

single_pressure_value_real = (pres1 * (1.0_r8 - fract)) + &
pres2 * (fract)

end function single_pressure_value_real

!-----------------------------------------------------------------------
!> return the level indices and fraction across the level.
!> level 1 is model top, level N is model bottom.
Expand Down

0 comments on commit 1982c75

Please sign in to comment.