Skip to content

Commit

Permalink
Merge pull request NCAR#78 from SamuelTrahanNOAA/feature/per-timestep…
Browse files Browse the repository at this point in the history 
…-pgr-like-wrf

per-timestep diagnostic output: FV3 version of WRF noise parameter
  • Loading branch information
@DomHeinzeller
DomHeinzeller committed Apr 2, 2021
2 parents e2108f0 + 837f7b4 commit 6814d1a
Show file tree
Hide file tree
Showing 3 changed files with 116 additions and 4 deletions.
90 changes: 90 additions & 0 deletions atmos_model.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -226,6 +226,7 @@ module atmos_model_mod

subroutine update_atmos_radiation_physics (Atmos)
!-----------------------------------------------------------------------
implicit none
type (atmos_data_type), intent(in) :: Atmos
!--- local variables---
integer :: nb, jdat(8), rc, ierr
Expand Down Expand Up @@ -358,6 +359,10 @@ subroutine update_atmos_radiation_physics (Atmos)

endif

! Per-timestep diagnostics must be after physics but before
! flagging the first timestep.
call atmos_timestep_diagnostics(Atmos)

! Update flag for first time step of time integration
GFS_control%first_time_step = .false.

Expand All @@ -366,6 +371,91 @@ end subroutine update_atmos_radiation_physics
! </SUBROUTINE>


!#######################################################################
! <SUBROUTINE NAME="atmos_timestep_diagnostics">
!
! <OVERVIEW>
! Calculates per-timestep, domain-wide, diagnostic, information and
! prints to stdout from master rank. Must be called after physics
! update but before first_time_step flag is cleared.
! </OVERVIEW>

! <TEMPLATE>
! call atmos_timestep_diagnostics (Atmos)
! </TEMPLATE>

! <INOUT NAME="Atmos" TYPE="type(atmos_data_type)">
! Derived-type variable that contains fields needed by the flux exchange module.
! These fields describe the atmospheric grid and are needed to
! compute/exchange fluxes with other component models. All fields in this
! variable type are allocated for the global grid (without halo regions).
! </INOUT>
subroutine atmos_timestep_diagnostics(Atmos)
use mpi
implicit none
type (atmos_data_type), intent(in) :: Atmos
!--- local variables---
integer :: i, nb, count, ierror
! double precision ensures ranks and sums are not truncated
! regardless of compilation settings
double precision :: pdiff, psum, pcount, maxabs, pmaxloc(7), adiff
double precision :: sendbuf(2), recvbuf(2), global_average

if(GFS_control%print_diff_pgr) then
if(.not. GFS_control%first_time_step) then
pmaxloc = 0.0d0
recvbuf = 0.0d0
psum = 0.0d0
pcount = 0.0d0
maxabs = 0.0d0

! Put pgr stats in pmaxloc, psum, and pcount:
pmaxloc(1) = GFS_Control%tile_num
do nb = 1,ATM_block%nblks
count = size(GFS_data(nb)%Statein%pgr)
do i=1,count
pdiff = GFS_data(nb)%Statein%pgr(i)-GFS_data(nb)%Intdiag%old_pgr(i)
adiff = abs(pdiff)
psum = psum+adiff
if(adiff>=maxabs) then
maxabs=adiff
pmaxloc(2:3)=(/ ATM_block%index(nb)%ii(i), ATM_block%index(nb)%jj(i) /)
pmaxloc(4:7)=(/ pdiff, GFS_data(nb)%Statein%pgr(i), &
GFS_data(nb)%Grid%xlat(i), GFS_data(nb)%Grid%xlon(i) /)
endif
enddo
pcount = pcount+count
enddo

! Sum pgr stats from psum/pcount and convert to hPa/hour global avg:
sendbuf(1:2) = (/ psum, pcount /)
call MPI_Allreduce(sendbuf,recvbuf,2,MPI_DOUBLE_PRECISION,MPI_SUM,GFS_Control%communicator,ierror)
global_average = recvbuf(1)/recvbuf(2) * 36.0d0/GFS_control%dtp

! Get the pmaxloc for the global maximum:
sendbuf(1:2) = (/ maxabs, dble(GFS_Control%me) /)
call MPI_Allreduce(sendbuf,recvbuf,1,MPI_2DOUBLE_PRECISION,MPI_MAXLOC,GFS_Control%communicator,ierror)
call MPI_Bcast(pmaxloc,size(pmaxloc),MPI_DOUBLE_PRECISION,nint(recvbuf(2)),GFS_Control%communicator,ierror)

if(GFS_Control%me == GFS_Control%master) then
2933 format('At forecast hour ',F9.3,' mean abs pgr change is ',F16.8,' hPa/hr')
2934 format(' max abs change ',F15.10,' bar at tile=',I0,' i=',I0,' j=',I0)
2935 format(' pgr at that point',F15.10,' bar lat=',F12.6,' lon=',F12.6)
print 2933, GFS_control%fhour, global_average
print 2934, pmaxloc(4)*1d-5, nint(pmaxloc(1:3))
print 2935, pmaxloc(5)*1d-5, pmaxloc(6:7)*57.29577951308232d0 ! 180/pi
endif
endif
! old_pgr is updated every timestep, including the first one where stats aren't printed:
do nb = 1,ATM_block%nblks
GFS_data(nb)%Intdiag%old_pgr=GFS_data(nb)%Statein%pgr
enddo
endif

!-----------------------------------------------------------------------
end subroutine atmos_timestep_diagnostics
! </SUBROUTINE>

!#######################################################################
! <SUBROUTINE NAME="atmos_model_init">
!
Expand Down
17 changes: 13 additions & 4 deletions ccpp/data/GFS_typedefs.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1146,9 +1146,10 @@ module GFS_typedefs
integer :: npsdelt !< the index of surface air pressure at the previous timestep for Z-C MP in phy_f2d
integer :: ncnvwind !< the index of surface wind enhancement due to convection for MYNN SFC and RAS CNV in phy f2d

!--- debug flag
!--- debug flags
logical :: debug
logical :: pre_rad !< flag for testing purpose
logical :: print_diff_pgr !< print average change in pgr every timestep (does not need debug flag)

!--- variables modified at each time step
integer :: ipt !< index for diagnostic printout point
Expand Down Expand Up @@ -1676,6 +1677,7 @@ module GFS_typedefs
!< for black carbon, organic carbon, and sulfur dioxide ( ug/m**2/s )
real (kind=kind_phys), pointer :: aecm (:,:) => null() !< instantaneous aerosol column mass densities for
!< pm2.5, black carbon, organic carbon, sulfate, dust, sea salt ( g/m**2 )
real (kind=kind_phys), pointer :: old_pgr(:) => null() !< pgr at last timestep

! Auxiliary output arrays for debugging
real (kind=kind_phys), pointer :: aux2d(:,:) => null() !< auxiliary 2d arrays in output (for debugging)
Expand Down Expand Up @@ -3380,9 +3382,10 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
logical :: iau_filter_increments = .false. !< filter IAU increments
logical :: iau_drymassfixer = .false. !< IAU dry mass fixer

!--- debug flag
!--- debug flags
logical :: debug = .false.
logical :: pre_rad = .false. !< flag for testing purpose
logical :: print_diff_pgr = .false. !< print average change in pgr every timestep

! max and min lon and lat for critical relative humidity
integer :: max_lon=5000, max_lat=2000, min_lon=192, min_lat=94
Expand Down Expand Up @@ -3506,7 +3509,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
iau_delthrs,iaufhrs,iau_inc_files,iau_filter_increments, &
iau_drymassfixer, &
!--- debug options
debug, pre_rad, &
debug, pre_rad, print_diff_pgr, &
!--- parameter range for critical relative humidity
max_lon, max_lat, min_lon, min_lat, rhcmax, &
phys_version, &
Expand Down Expand Up @@ -4198,9 +4201,10 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
Model%iau_drymassfixer = iau_drymassfixer
if(Model%me==0) print *,' model init,iaufhrs=',Model%iaufhrs

!--- debug flag
!--- debug flags
Model%debug = debug
Model%pre_rad = pre_rad
Model%print_diff_pgr = print_diff_pgr

!--- tracer handling
Model%ntrac = size(tracer_names)
Expand Down Expand Up @@ -5610,6 +5614,11 @@ subroutine diag_create (Diag, IM, Model)
!
logical, save :: linit

if(Model%print_diff_pgr) then
allocate(Diag%old_pgr(IM))
Diag%old_pgr = clear_val
endif

!--- Radiation
allocate (Diag%fluxr (IM,Model%nfxr))
allocate (Diag%topfsw (IM))
Expand Down
13 changes: 13 additions & 0 deletions ccpp/data/GFS_typedefs.meta
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4637,6 +4637,12 @@
units = flag
dimensions = ()
type = logical
[print_diff_pgr]
standard_name = flag_to_print_pgr_differences_every_timestep
long_name = flag to print pgr differences every timestep
units = flag
dimensions = ()
type = logical
[ipt]
standard_name = index_for_diagnostic_printout
long_name = horizontal index for point used for diagnostic printout
Expand Down Expand Up @@ -7410,6 +7416,13 @@
type = real
kind = kind_phys
active = (number_of_2d_auxiliary_arrays > 0)
[old_pgr]
standard_name = surface_air_pressure_from_previous_timestep
long_name = surface air pressure from previous timestep
units = Pa
dimensions = (horizontal_loop_extent)
type = real
kind = kind_phys

########################################################################
[ccpp-table-properties]
Expand Down

0 comments on commit 6814d1a

Please sign in to comment.