#include <misc.h>
#include <params.h>
subroutine physpkg(phys_state, gw, ztodt, phys_tend, pbuf) 1,50
!-----------------------------------------------------------------------
!
! Purpose:
! Loop over time, calling driving routines for physics
!
! Method:
! COUP_CSM and must be checked in order to invoke the proper calling
! sequence for running the CSM model
!
! Author:
! Original version: CCM3
!-----------------------------------------------------------------------
use shr_kind_mod, only: r8 => shr_kind_r8
use pmgrid, only: plon, plat, masterproc
use tracers, only: set_state_pdry
use ppgrid, only: pcols, pver
use buffer, only: pblht, tpert, qpert, qrs, qrl
use check_energy, only: check_energy_gmean
use comsrf
use comsrfdiag
#ifdef COUP_CSM
use ccsm_msg, only: ccsmave, dorecv, dosend, ccsmsnd, ccsmrcv
#else
use atm_lndMod, only: atmlnd_drv
#endif
#ifdef SPMD
use mpishorthand
#endif
use phys_buffer, only: pbuf_size_max, pbuf_fld, pbuf_allocate, pbuf_deallocate, &
pbuf_update_tim_idx
use phys_grid, only: get_ncols_p, get_lat_all_p, get_lon_all_p
use physics_types, only: physics_state, physics_tend
use diagnostics, only: diag_surf
use time_manager, only: get_nstep, is_first_step, is_first_restart_step, &
is_end_curr_month, get_curr_date, is_end_curr_day
use physconst, only: stebol
use dycore, only: dycore_is
#if (!defined COUP_CSM)
use ice_constants, only: TfrezK
#endif
use history, only: outfld
#if (!defined COUP_CSM)
#if (!defined COUP_SOM)
use sst_data, only: sstint
use ice_data, only: iceint
#endif
#endif
!-----------------------------------------------------------------------
implicit none
!-----------------------------------------------------------------------
#include <comctl.h>
!-----------------------------------------------------------------------
#include <comsol.h>
!-----------------------------------------------------------------------
!
! Input arguments
!
real(r8), intent(in) :: gw(plat) ! Gaussian weights
real(r8), intent(in) :: ztodt ! physics time step unless nstep=0
!
! Input/Output arguments
!
type(physics_state), intent(inout), dimension(begchunk:endchunk) :: phys_state
type(physics_tend ), intent(inout), dimension(begchunk:endchunk) :: phys_tend
type(pbuf_fld), intent(inout), dimension(pbuf_size_max) :: pbuf
!
!---------------------------Local workspace-----------------------------
!
integer :: i,m,lat,c,lchnk ! indices
integer :: lats(pcols) ! array of latitude indices
integer :: lons(pcols) ! array of longitude indices
integer :: ncol ! number of columns
integer :: nstep ! current timestep number
integer :: ncdate ! current date in integer format [yyyymmdd]
integer :: ncsec ! current time of day [seconds]
integer :: yr, mon, day ! year, month, and day components of a date
real(r8) fsds(pcols,begchunk:endchunk) ! Surface solar down flux
real(r8) :: tmp(pcols,begchunk:endchunk)
!-----------------------------------------------------------------------
call t_startf ('physpkg_st')
nstep = get_nstep()
call pbuf_allocate('physpkg')
! Compute total energy of input state and previous output state
call t_startf ('chk_en_gmean')
call check_energy_gmean(phys_state, pbuf, ztodt, nstep)
call t_stopf ('chk_en_gmean')
!-----------------------------------------------------------------------
! Advance time information
!-----------------------------------------------------------------------
call advnce()
call t_stopf ('physpkg_st')
!
! set the state vector dry quantities
!
if ( .not. dycore_is('LR') )then
! for LR, this is done in d_p_coupling since dynamics is called first
!$OMP PARALLEL DO PRIVATE (C,NCOL)
do c=begchunk, endchunk
call set_state_pdry( phys_state(c))
end do
endif ! .not. dycore_is('LR')
#ifdef TRACER_CHECK
call gavglook ('before tphysbc DRY', phys_state, gw)
#endif
!-----------------------------------------------------------------------
! Tendency physics before flux coupler invokation
!-----------------------------------------------------------------------
!
#if (defined BFB_CAM_SCAM_IOP )
do c=begchunk, endchunk
call outfld('Tg',srfflx_state2d(c)%ts,pcols ,c )
end do
#endif
call t_startf ('bc_physics')
!$OMP PARALLEL DO PRIVATE (C,NCOL)
do c=begchunk, endchunk
call t_startf ('tphysbc')
call tphysbc (ztodt, pblht(1,c), tpert(1,c), &
srfflx_state2d(c)%ts, srfflx_state2d(c)%sst, &
qpert(1,1,c), surface_state2d(c)%precl, &
surface_state2d(c)%precc, surface_state2d(c)%precsl,&
surface_state2d(c)%precsc, &
srfflx_state2d(c)%asdir, srfflx_state2d(c)%asdif, &
srfflx_state2d(c)%aldir, srfflx_state2d(c)%aldif, &
snowhland(1,c), &
qrs(1,1,c), qrl(1,1,c), surface_state2d(c)%flwds, &
fsns(1,c), fsnt(1,c), &
flns(1,c), flnt(1,c), srfflx_state2d(c)%lwup, &
surface_state2d(c)%srfrad, surface_state2d(c)%sols, &
surface_state2d(c)%soll, surface_state2d(c)%solsd, &
surface_state2d(c)%solld, &
phys_state(c), phys_tend(c), &
pbuf, prcsnw(1,c), fsds(1,c), landm(1,c), landfrac(1,c), &
ocnfrac(1,c),icefrac(1,c))
call t_stopf ('tphysbc')
if (dosw .or. dolw) then
call output_flns_fsns_fluxes(surface_state2d(c),c)
end if
#if ( ! defined COUP_CSM )
!
! zero surface fluxes at beginning of each time step. Land Ocean and Ice
! processes will will write into process specific flux variables
! at the end of the time step these separate fluxes will be combined over the
! entire grid
!
call srfflx_state_reset (srfflx_state2d(c))
#endif
end do
call t_stopf ('bc_physics')
#ifdef TRACER_CHECK
call gavglook ('between DRY', phys_state, gw)
#endif
#if ( ! defined COUP_CSM )
!
!-----------------------------------------------------------------------
! Determine surface quantities - no flux coupler
!-----------------------------------------------------------------------
!
if (.not. aqua_planet) then
!
! Call land model driving routine
!
#ifdef TIMING_BARRIERS
call t_startf ('sync_tphysbc_lnd')
call mpibarrier (mpicom)
call t_stopf ('sync_tphysbc_lnd')
#endif
call t_startf ('atmlnd_drv')
#if ( defined SCAM )
if (landfrac(1,begchunk).gt.0) &
#endif
call atmlnd_drv(nstep, iradsw, eccen, obliqr, lambm0,&
mvelpp,surface_state2d,srfflx_parm2d)
call t_stopf ('atmlnd_drv')
#ifdef TIMING_BARRIERS
call t_startf ('sync_after_lnd')
call mpibarrier (mpicom)
call t_stopf ('sync_after_lnd')
#endif
!
! save off albedos and longwave for som offline vars
!
!$OMP PARALLEL DO PRIVATE (C,NCOL,I)
do c=begchunk,endchunk
ncol = get_ncols_p(c)
do i=1,ncol
if (landfrac(i,c) > 0.) then
asdirlnd(i,c) = srfflx_parm2d(c)%asdir(i)
asdiflnd(i,c) = srfflx_parm2d(c)%asdif(i)
aldirlnd(i,c) = srfflx_parm2d(c)%aldir(i)
aldiflnd(i,c) = srfflx_parm2d(c)%aldif(i)
lwuplnd(i,c) = srfflx_parm2d(c)%lwup(i)
else
asdirlnd(i,c) = 0.
asdiflnd(i,c) = 0.
aldirlnd(i,c) = 0.
aldiflnd(i,c) = 0.
lwuplnd(i,c) = 0.
end if
end do
!
!output shf/lhf fluxes for land model
!
call output_shf_lhf_fluxes(srfflx_parm2d(c), c, ncol, landfrac(1,c), 'LND')
call update_srf_fluxes (srfflx_state2d(c), srfflx_parm2d(c), landfrac(1,c), ncol)
end do
end if ! end of not aqua_planet if block
#if (defined COUP_SOM)
!
! Set ocean surface quantities - ocn model internal to atm
!
if (is_end_curr_day ()) then
call print_coverage ('icefrac', ' million km^2', icefrac, 1.d-12)
do c=begchunk,endchunk
ncol = get_ncols_p(c)
do i=1,ncol
tmp(i,c) = icefrac(i,c)*sicthk(i,c)
end do
end do
call print_coverage ('icevol ', ' 10^13m^3', tmp, 1.d-13)
do c=begchunk,endchunk
ncol = get_ncols_p(c)
do i=1,ncol
tmp(i,c) = icefrac(i,c)*snowhice(i,c)
end do
end do
call print_coverage ('snowvol', ' 10^13m^3', tmp, 1.d-13)
end if
call t_startf ('somint')
call somint ()
call t_stopf ('somint')
call t_startf ('somoce')
call somoce (surface_state2d, srfflx_parm2d_ocn)
call t_stopf ('somoce')
#else
call t_startf ('sstint')
call sstint ()
call t_stopf ('sstint')
!
! iceint may change ocean fraction, so call it before camoce
!
call t_startf ('iceint')
call iceint ()
call t_stopf ('iceint')
call t_startf ('camoce')
call camoce (surface_state2d, srfflx_parm2d_ocn)
call t_stopf ('camoce')
#endif
!
! Set ice surface quantities - icn model internal to atm
!
call t_startf('camice')
call camice (surface_state2d, srfflx_parm2d)
call t_stopf('camice')
!
! output shf/lhf fluxes for ice/ocn/som_offline
!
!$OMP PARALLEL DO PRIVATE (C, NCOL, I)
do c=begchunk,endchunk
ncol = get_ncols_p(c)
do i=1,ncol
if(icefrac(i,c) > 0.) then
tsice_rad(i,c) = sqrt(sqrt(srfflx_parm2d(c)%lwup(i)/stebol))
else
tsice_rad(i,c) = TfrezK
endif
end do
call output_shf_lhf_fluxes (srfflx_parm2d(c), c, ncol, icefrac(1,c), 'ICE')
call output_shf_lhf_fluxes (srfflx_parm2d_ocn(c), c, ncol, ocnfrac(1,c), 'OCN')
call output_shfoi_lhfoi_fluxes (srfflx_parm2d_ocn(c), srfflx_parm2d(c), c)
!JR SOM case: Have to wait to call update routine till after both ocean and ice have
!JR operated, since the fractions can change internal to the parameterization
do i = 1, ncol
srfflx_state2d(c)%sst(i) = srfflx_parm2d_ocn(c)%ts(i)
enddo
call update_srf_fluxes (srfflx_state2d(c), srfflx_parm2d_ocn(c), ocnfrac(1,c), ncol)
call update_srf_fluxes (srfflx_state2d(c), srfflx_parm2d(c), icefrac(1,c), ncol)
end do
#endif
#if ( defined COUP_CSM )
!
!-----------------------------------------------------------------------
! Determine surface quantities using csm flux coupler
!-----------------------------------------------------------------------
!
! If send data to flux coupler only on radiation time steps:
!
if (flxave) then
!
! Average the precipitation input to lsm between radiation calls.
!
call ccsmave(iradsw, nstep, dosw)
!
! Use solar radiation flag to determine data exchange steps
! with flux coupler. This processes are not independent since
! instantaneous radiative fluxes are passed, valid over the
! interval to the next radiation calculation. The same
! considerations apply to the long and shortwave fluxes, so
! the intervals must be the same. Data is received from the
! coupler one step after it is sent.
!
if (nstep == 0) then
dorecv = .true.
dosend = .true.
else if (nstep == 1) then
dorecv = .false.
dosend = .false.
else if ( (nstep == 2) .and. (iradsw == 1) ) then
dorecv = .true.
dosend = dosw
else
dorecv = dosend
dosend = dosw
end if
endif
!
! If send data to flux coupler on every time step
!
if (.not. flxave) then
if (nstep /= 1) then
dorecv = .true.
dosend = .true.
else
dorecv = .false.
dosend = .false.
endif
endif
!
! Send/recv data to/from the csm flux coupler.
!
if (dosend) call ccsmsnd ( )
if (dorecv) call ccsmrcv ( )
#endif
!
!-----------------------------------------------------------------------
! Tendency physics after coupler
! Not necessary at terminal timestep.
!-----------------------------------------------------------------------
!
call t_startf ('ac_physics')
!$OMP PARALLEL DO PRIVATE (C, NCOL)
do c=begchunk,endchunk
ncol = get_ncols_p(c)
!
! surface diagnostics for history files
!
call diag_surf (srfflx_state2d(c), surface_state2d(c), icefrac(1,c), ocnfrac(1,c), landfrac(1,c), &
sicthk(1,c), snowhland(1,c), snowhice(1,c), tsice(1,c), trefmxav(1,c), &
trefmnav(1,c) )
call t_startf ('tphysac')
call tphysac (ztodt, pblht(1,c), qpert(1,1,c), tpert(1,c), srfflx_state2d(c)%shf, &
srfflx_state2d(c)%wsx,srfflx_state2d(c)%wsy, srfflx_state2d(c)%cflx, sgh(1,c), srfflx_state2d(c)%lhf, &
landfrac(1,c), snowhland(1,c),srfflx_state2d(c)%tref, surface_state2d(c)%precc, surface_state2d(c)%precl, &
surface_state2d(c)%precsc, surface_state2d(c)%precsl, phys_state(c), phys_tend(c), pbuf, &
ocnfrac(1,c), fsds(1,c), icefrac(1,c), fv(1,c), ram1(1,c))
call t_stopf ('tphysac')
end do ! Chunk loop
call t_stopf('ac_physics')
#ifdef TRACER_CHECK
call gavglook ('after tphysac FV:WET)', phys_state, gw )
#endif
call pbuf_deallocate('physpkg')
call pbuf_update_tim_idx()
end subroutine physpkg
subroutine gavglook (title, state, gw) 3,17
!
! process info from state vectors. Only useful when data in all chunks are in sync
! e.g. before and after tphysac and tphysbc
!
use physics_types, only: physics_state, physics_tend
! use comsrf, only: srfflx_state
use phys_grid, only: get_ncols_p, get_lat_all_p, get_lon_all_p
use ppgrid, only: begchunk, endchunk
use ppgrid, only: pcols, pver
use constituents, only: ppcnst, cnst_name
use shr_kind_mod, only: r8 => shr_kind_r8
use pmgrid, only: plon, plat, masterproc
use physconst, only: gravit
use time_manager, only: dtime
#ifdef SPMD
use mpishorthand
#endif
implicit none
! arguments
character(len=*), intent(in) :: title
type(physics_state), intent(in), dimension(begchunk:endchunk) :: state
real(r8), intent(in) :: gw(plat) ! Gaussian weights
! local
integer i, lat, c, lon, k
integer :: lats(pcols) ! array of latitude indices
integer :: lons(pcols) ! array of longitude indices
integer m
integer :: ncol ! number of columns
real(r8) twodfld(plon,plat,ppcnst) ! summed at each grid point
real(r8) twodfle(plon,plat,ppcnst) ! summed at each grid point
real(r8) twodflx(plon,plat,ppcnst) ! summed at each grid point
real(r8) twodfly(plon,plat,ppcnst) ! summed at each grid point
#ifdef SPMD
real(r8) :: twodfld_glob(plon,plat,ppcnst) ! global summed at each grid point
real(r8) :: twodfle_glob(plon,plat,ppcnst) ! global summed at each grid point
real(r8) :: twodflx_glob(plon,plat,ppcnst) ! global summed at each grid point
real(r8) :: twodfly_glob(plon,plat,ppcnst) ! global summed at each grid point
#endif
real(r8) :: zonal(plat), zonalw(plat) ! summed along each latitude
real(r8) gavg, gavgw
real(r8) col, wmin, wmax, colw
!!--------------------------------------------------------
! operations on each processor
twodfld(:,:,:) = 0.
twodfle(:,:,:) = 0.
twodflx(:,:,:) = 0.
twodfly(:,:,:) = 0.
do c=begchunk, endchunk
ncol = get_ncols_p(c)
call get_lat_all_p(c, ncol, lats)
call get_lon_all_p(c, ncol, lons)
do m = 1,ppcnst
do i=1,ncol
lat = lats(i)
lon = lons(i)
col = 0.
colw = 0.
! fluxcol = 0.
wmax = -1.e36
wmin = 1.e36
do k = 1,pver
col = col + state(c)%pdeldry(i,k)*state(c)%q(i,k,m)*gw(lats(i))
colw = colw + state(c)%pdel(i,k) *state(c)%q(i,k,m)*gw(lats(i))
wmax = max(wmax,state(c)%q(i,k,m))
wmin = min(wmin,state(c)%q(i,k,m))
end do ! k
col = col/gravit
colw = colw/gravit
twodfld(lons(i),lats(i),m) = twodfld(lons(i),lats(i),m) + col
twodfle(lons(i),lats(i),m) = twodfle(lons(i),lats(i),m) + colw
twodflx(lons(i),lats(i),m) = twodflx(lons(i),lats(i),m) + wmin
twodfly(lons(i),lats(i),m) = twodfly(lons(i),lats(i),m) + wmax
enddo ! i
enddo ! m
end do ! c
! move data to masterproc
#ifdef SPMD
#ifdef TIMING_BARRIERS
call mpibarrier (mpicom)
#endif
call mpisum(twodfld, twodfld_glob, plon*plat*ppcnst, mpir8, 0, mpicom)
call mpisum(twodfle, twodfle_glob, plon*plat*ppcnst, mpir8, 0, mpicom)
call mpisum(twodflx, twodflx_glob, plon*plat*ppcnst, mpir8, 0, mpicom)
call mpisum(twodfly, twodfly_glob, plon*plat*ppcnst, mpir8, 0, mpicom)
if (masterproc) then
twodfld(:,:,:) = twodfld_glob(:,:,:)
twodfle(:,:,:) = twodfle_glob(:,:,:)
twodflx(:,:,:) = twodflx_glob(:,:,:)
twodfly(:,:,:) = twodfly_glob(:,:,:)
endif
#endif
! process the data
if (masterproc) then
do m = 1,ppcnst
wmax = -1.e36
wmin = 1.e36
do lat=1,plat
zonal(lat) = 0.
zonalw(lat) = 0.
do i=1,plon
zonal(lat) = zonal(lat) + twodfld(i,lat,m)
zonalw(lat) = zonalw(lat) + twodfle(i,lat,m)
wmax = max(wmax,twodfly(i,lat,m))
wmin = min(wmin,twodflx(i,lat,m))
end do
end do
gavg = 0.
gavgw = 0.
do lat=1,plat
gavg = gavg + zonal(lat)
gavgw = gavgw + zonalw(lat)
end do
gavg = gavg/(2.*plon)
gavgw = gavgw/(2.*plon)
write (6,66) trim(title)//' m=',m,'name='//trim(cnst_name(m))//' gavg dry, wet, min, max ' &
, gavg, gavgw,wmin,wmax
66 format (a24,i2,a36,1p,4g25.14)
end do
endif
end subroutine gavglook