!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! !!
!! GNU General Public License !!
!! !!
!! This file is part of the Flexible Modeling System (FMS). !!
!! !!
!! FMS is free software; you can redistribute it and/or modify !!
!! it and are expected to follow the terms of the GNU General Public !!
!! License as published by the Free Software Foundation. !!
!! !!
!! FMS is distributed in the hope that it will be useful, !!
!! but WITHOUT ANY WARRANTY; without even the implied warranty of !!
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !!
!! GNU General Public License for more details. !!
!! !!
!! You should have received a copy of the GNU General Public License !!
!! along with FMS; if not, write to: !!
!! Free Software Foundation, Inc. !!
!! 59 Temple Place, Suite 330 !!
!! Boston, MA 02111-1307 USA !!
!! or see: !!
!! http://www.gnu.org/licenses/gpl.txt !!
!! !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
module fv_diagnostics
! Programmer: sjl@gfdl.noaa.gov, Oct 30, 2003
use diag_manager_mod, only: diag_axis_init, register_diag_field, &
register_static_field, send_data
use constants_mod, only: GRAV, cp_air, rdgas, kappa, WTMAIR, WTMCO2
use time_manager_mod, only: time_type, get_date, get_time
use fms_mod, only: error_mesg, FATAL, stdlog, write_version_number
use mpp_domains_mod, only: mpp_define_domains, mpp_domains_init, domain2D
! use fv_pack, only: nlon, mlat, nlev, beglat, endlat, u, v, pt, delp, q, &
! omga, peln, pe, ua, va, phis, ptop, ak, bk, ks, ncnst, &
! ng_d, ng_s, ps, pk, pkz, lon, lat, lonb, latb, &
! age_tracer, age_time, fv_domain, print_freq, &
! drymadj, tracer_mass, &
! master, get_eta_level, nt_phys, do_fms_tracer_manager
use fv_pack, only: nlon, mlat, beglon, endlon, beglat, endlat, &
ptop, ak, bk, ks, ng_d, ng_s, lon, lat, lonb, latb, &
age_tracer, age_time, fv_domain, print_freq, drymadj, tracer_mass, &
master, get_eta_level, nt_phys, do_fms_tracer_manager, area
#ifdef MARS_GCM
use fv_pack, only: p_ref
use fv_phys_mod, only: mars_mass_budget
#endif MARS_GCM
! --- tracer manager ---
use tracer_manager_mod, only : get_tracer_names, get_number_tracers
use field_manager_mod, only : MODEL_ATMOS
use mpp_mod, only: mpp_sync
use fv_arrays_mod, only: fv_stack_push, fv_array_sync, fv_print_chksums
use pmaxmin_mod, only: pmaxmin, pmaxming
implicit none
integer, save :: id_ps, id_ua, id_va, id_wa, id_ta, id_pt, &
id_divg, id_vort, id_hght, id_pv, id_usus, id_tsts, &
id_vsts, id_usvs, id_h500, id_slp, id_p5km, id_aoa, &
id_dnflux, id_upflux, id_area
integer, allocatable :: id_tracer(:), id_tracer_tend(:)
! 2008/04/09 jgj: add tracer dry mass and volume mixing ratios
integer, allocatable :: id_tracer_dmmr(:)
integer, allocatable :: id_tracer_dvmr(:)
real :: missing_value = -1.e10
real, save :: ginv
logical used
type(time_type) :: fv_time
real, allocatable, save :: phalf(:)
save fv_time
!-----------------------------------------------------------------------
character(len=128) :: version = '$Id: fv_diagnostics.F90,v 17.0 2009/07/21 02:53:23 fms Exp $'
character(len=128) :: tagname = '$Name: mom4p1_pubrel_dec2009_nnz $'
contains
subroutine fv_diag_init(axes, Time)
#include "fv_arrays.h"
! FV dynamics specific diagnostics
type(time_type), intent(in) :: Time
integer :: id_lonb, id_lon, id_latb, id_lat, id_phalf, id_pfull
integer :: id_bk, id_pk, id_zsurf
character(len=8) :: mod_name = 'dynamics'
integer, intent(out) :: axes(4)
logical used
integer i, j
integer :: log_unit
real zsurf(nlon,beglat:endlat)
real pfull(nlev)
real qmin, qmax
real :: vrange(2), trange(2), arange(2), slprange(2)
! tracers
character(len=128) :: tname
character(len=256) :: tlongname, tunits
integer :: ntprog
#include "fv_point.inc"
log_unit = stdlog()
call fv_print_chksums( 'Entering fv_diag_init' )
! valid range for some fields
vrange = (/ -330., 350. /) ! winds
#ifdef MARS_GCM
trange = (/ 50., 350. /) ! temperature
#else
trange = (/ 100., 350. /) ! temperature
#endif
arange = (/ 0., 20. /) ! age-of-air
slprange = (/800., 1200./) ! sea-level-pressure
ginv = 1./GRAV
fv_time = Time
allocate ( phalf(nlev+1) )
#ifdef MARS_GCM
call get_eta_level(nlev, p_ref, pfull, phalf, 0.01)
#else
call get_eta_level(nlev, 1.E5, pfull, phalf, 0.01)
#endif
id_lonb=diag_axis_init('lonb', lonb, 'degrees_E', 'x', 'longitude edges', &
set_name=mod_name , Domain2=fv_domain )
id_lon =diag_axis_init('lon', lon, 'degrees_E', 'x', 'longitude', &
set_name=mod_name, edges=id_lonb, Domain2=fv_domain )
id_latb=diag_axis_init('latb', latb, 'degrees_N', 'y', 'latitude edges', &
set_name=mod_name, Domain2=fv_domain )
id_lat =diag_axis_init('lat', lat, 'degrees_N', 'y', 'latitude', &
set_name=mod_name, edges=id_latb, Domain2=fv_domain)
id_phalf = diag_axis_init('phalf', phalf, 'mb', 'z', &
'ref half pressure level', direction=-1, set_name=mod_name )
id_pfull = diag_axis_init('pfull', pfull, 'mb', 'z', &
'ref full pressure level', direction=-1, set_name=mod_name, &
edges=id_phalf)
axes(1) = id_lon
axes(2) = id_lat
axes(3) = id_pfull
axes(4) = id_phalf
!---- register static fields -------
id_bk = register_static_field ( mod_name, 'bk', (/id_phalf/), &
'vertical coordinate sigma value', 'none' )
id_pk = register_static_field ( mod_name, 'pk', (/id_phalf/), &
'pressure part of the hybrid coordinate', 'pascal' )
id_area = register_static_field ( mod_name, 'area', axes(1:2), &
& 'cell area elements', 'm^2')
if( id_area > 0 ) &
& used = send_data( id_area, area(beglon:endlon,beglat:endlat), Time )
do j=beglat,endlat
do i=1,nlon
zsurf(i,j) = ginv * phis(i,j)
enddo
enddo
id_zsurf = register_static_field ( mod_name, 'zsurf', axes(1:2), &
'surface height', 'm' )
!--- Send static data
if ( id_bk > 0 ) used = send_data ( id_bk, bk, Time )
if ( id_pk > 0 ) used = send_data ( id_pk, ak, Time )
if ( id_zsurf > 0 ) &
used = send_data( id_zsurf, zsurf(beglon:endlon,beglat:endlat), Time )
! if ( master ) then
! do i=1,nlev+1
! write(6,*) i, ak(i), bk(i)
! enddo
! endif
!--------------------------------------------------------------
! Register main prognostic fields: ps, (u,v), t, omega (dp/dt)
!--------------------------------------------------------------
id_ps = register_diag_field ( mod_name, 'ps', axes(1:2), Time, &
'surface pressure', 'Pa', missing_value=missing_value )
id_ua = register_diag_field ( mod_name, 'ucomp', axes(1:3), Time, &
'zonal wind', 'm/sec', missing_value=missing_value, range=vrange )
id_va = register_diag_field ( mod_name, 'vcomp', axes(1:3), Time, &
'meridional wind', 'm/sec', missing_value=missing_value, range=vrange)
id_wa = register_diag_field ( mod_name, 'omega', axes(1:3), Time, &
'omega', 'pa/sec', missing_value=missing_value )
id_dnflux = register_diag_field (mod_name, 'dnflux', axes(1:3), Time, &
'downward mass flux', 'pa/sec', missing_value=missing_value )
id_upflux = register_diag_field (mod_name, 'upflux', axes(1:3), Time, &
'upward mass flux', 'pa/sec', missing_value=missing_value )
id_ta = register_diag_field ( mod_name, 'temp', axes(1:3), Time, &
'temperature', 'deg_K', missing_value=missing_value, range=trange )
id_pt = register_diag_field ( mod_name, 'pt_dry', axes(1:3), Time, &
'dry poten temperature', 'deg_K', missing_value=missing_value )
!------------------
! Register tracers
!------------------
allocate(id_tracer(max(ncnst, nt_phys)))
id_tracer(:) = 0
! 2008/04/09 jgj: add tracer dry mass and volume mixing ratios
allocate(id_tracer_dmmr(max(ncnst, nt_phys)))
id_tracer_dmmr(:) = 0
allocate(id_tracer_dvmr(max(ncnst, nt_phys)))
id_tracer_dvmr(:) = 0
if(.NOT. do_fms_tracer_manager) then !{
! no field table
id_tracer(1) = register_diag_field ( mod_name, 'sphum', axes(1:3), &
Time, 'sphum', 'mass/air_mass', &
missing_value=missing_value )
id_tracer(2) = register_diag_field ( mod_name, 'liq_wat', axes(1:3), &
Time, 'cloud liquid water', 'mass/air_mass', &
missing_value=missing_value )
id_tracer(3) = register_diag_field ( mod_name, 'ice_wat', axes(1:3), &
Time, 'cloud ice water', 'mass/air_mass',&
missing_value=missing_value )
id_tracer(4) = register_diag_field ( mod_name, 'cld_amt', axes(1:3), &
Time, 'cloud craction', '%', &
missing_value=missing_value )
else !}{
! do_fms_traces == .TRUE.
#ifdef MARS_GCM
do i= 1, ncnst
#else
do i = 1, nt_phys
call get_tracer_names ( MODEL_ATMOS, i, tname, tlongname, tunits )
id_tracer(i) = register_diag_field ( mod_name, trim(tname), &
axes(1:3), Time, trim(tlongname), &
trim(tunits), missing_value=missing_value)
if (master) then
if (id_tracer(i) > 0) then
write(log_unit,'(a,a,a,a)') &
& 'Diagnostics available for tracer ',tname, &
' in module ', mod_name
end if
endif
enddo
! for higher tracer number we take the default missing values from FMS
do i = nt_phys+1, ncnst
#endif MARS_GCM
call get_tracer_names ( MODEL_ATMOS, i, tname, tlongname, tunits )
id_tracer(i) = register_diag_field ( mod_name, trim(tname), &
axes(1:3), Time, trim(tlongname), &
trim(tunits), missing_value=missing_value)
if (master) then
if (id_tracer(i) > 0) then
write(log_unit,'(a,a,a,a)') &
& 'Diagnostics available for tracer ',trim(tname), &
' in module ', trim(mod_name)
end if
endif
enddo
!-------
#ifdef MARS_GCM
#else
! jgj: add co2 tracer dry mass and volume mixing ratios
do i = nt_phys+1, ncnst
call get_tracer_names ( MODEL_ATMOS, i, tname, tlongname, tunits )
if (trim(tname).eq.'co2') then
id_tracer_dmmr(i) = register_diag_field ( mod_name, trim(tname)//'_dmmr', &
axes(1:3), Time, trim(tlongname)//" (dry mmr)", &
trim(tunits), missing_value=missing_value)
id_tracer_dvmr(i) = register_diag_field ( mod_name, trim(tname)//'_dvmr', &
axes(1:3), Time, trim(tlongname)//" (dry vmr)", &
'mol/mol', missing_value=missing_value)
if (master) then
if (id_tracer_dmmr(i) > 0) then
write(log_unit,'(a,a,a,a)') &
& 'Diagnostics available for co2 dry mmr ',trim(tname)//'_dmmr', &
' in module ', trim(mod_name)
end if
if (id_tracer_dvmr(i) > 0) then
write(log_unit,'(a,a,a,a)') &
& 'Diagnostics available for co2 dry vmr ',trim(tname)//'_dvmr', &
' in module ', trim(mod_name)
end if
endif
endif
enddo
#endif
!-------
endif !}
!----------------------
! Divergence on C grid
!----------------------
! Computed inside fvcore/trac2d()
id_divg = register_diag_field ( mod_name, 'divg', axes(1:3), Time, &
'divergence', '1/sec', missing_value=missing_value )
!----------------------
! Vorticity on D grid
!----------------------
id_vort = register_diag_field ( mod_name, 'vort', axes(1:3), Time, &
'vorticity', '1/sec', missing_value=missing_value )
!--------------------------
! Ertel Potential Vorticity
!--------------------------
id_pv = register_diag_field ( mod_name, 'pv', axes(1:3), Time, &
'potential vorticity', '1/sec', missing_value=missing_value )
!-------------------
! Age_of_air tracer
!-------------------
id_aoa = register_diag_field ( mod_name, 'aoa', axes(1:3), Time, &
'age_of_air tracer', 'sec', missing_value=missing_value, &
range=arange )
!---------------------------
! Height at model half level
!---------------------------
id_hght = register_diag_field ( mod_name, 'hght', &
(/id_lon,id_lat,id_phalf/), &
Time, 'height', 'm', missing_value=missing_value )
!--------------
! 500 mb Height
!--------------
id_h500 = register_diag_field ( mod_name, 'h500', axes(1:2), Time, &
'500-mb hght', 'm', missing_value=missing_value )
!-------------------
! Sea-level-pressure
!-------------------
id_slp = register_diag_field ( mod_name, 'slp', axes(1:2), Time, &
'sea-level pressure', 'mb', missing_value=missing_value, &
range=slprange )
!-------------------
! 5-km pressure
!-------------------
! This is analogous to 500-mb height. This shows the actual highs and lows
! at the 5-km physical height.
id_p5km = register_diag_field ( mod_name, 'p5km', axes(1:2), Time, &
'pressure at 5 km', 'mb', missing_value=missing_value )
!------------------
! 2nd moment stats:
!------------------
! u-wind variance
id_usus = register_diag_field ( mod_name, 'usus', axes(1:3), Time, &
'u-wind variance', '(m/sec)**2', missing_value=missing_value )
! temperature variance
id_tsts = register_diag_field ( mod_name, 'tsts', axes(1:3), Time, &
'temperature variance', 'K**2', missing_value=missing_value )
! eddy momentum transport
id_usvs = register_diag_field ( mod_name, 'usvs', axes(1:3), Time, &
'eddy momentum transport', '(m/sec)**2', missing_value=missing_value )
! eddy heat transport
id_vsts = register_diag_field ( mod_name, 'vsts', axes(1:3), Time, &
'eddy heat transport', '(m/sec)*K', missing_value=missing_value )
! make sure all send_data completed before one thread exits the routine
call fv_array_sync()
call fv_print_chksums( 'Exiting fv_diag_init' )
end subroutine fv_diag_init
!=================================================
subroutine fv_diag( Time, im, jm, km, jfirst, jlast, nq, zvir, dt_atmos, hs_phys)
use pv_module, only: vort_d, pv_entropy
#include "fv_arrays.h"
!INPUT PARAMETERS:
type(time_type), intent(in) :: Time
integer, intent(in):: im ! dimension in east-west
integer, intent(in):: jm ! dimension in North-South
integer, intent(in):: km ! number of Lagrangian layers
integer, intent(in):: jfirst ! starting latitude index for MPI
integer, intent(in):: jlast ! ending latitude index for MPI
integer, intent(in):: nq ! number of tracers
real, intent(in):: zvir
real, intent(in):: dt_atmos ! model time step in seconds
logical, intent(in):: hs_phys
integer, parameter:: kmax = 100
real height(kmax)
real log_p(kmax)
real qmax, qmin, fac
real gg
real gmean
real ztop
integer i, j, k, ic
integer yr, mon, dd, hr, mn, days, seconds
logical:: prt_minmax
real p0
! Local:
real :: a2(im,jfirst:jlast) ! work array for x-y diag
real :: ax(jfirst:jlast,km) ! work array for zonal means
! real, allocatable:: xx(:,:,:) ! work array for 2nd moments
! real, allocatable:: us(:,:,:)
! real, allocatable:: vs(:,:,:)
! real, allocatable:: wz(:,:,:)
real :: xx(im,jfirst:jlast,km)
real :: us(im,jfirst:jlast,km)
real :: vs(im,jfirst:jlast,km)
real :: wz(im,jfirst:jlast,km+1)
! 2008/04/09 jgj: add tracer dry vmr to compare to obs
real :: co2_dmmr(im,jfirst:jlast,km)
real :: co2_dvmr(im,jfirst:jlast,km)
! tracers
character(len=128) :: tname
character(len=256) :: tlongname, tunits
#ifdef use_shared_pointers
pointer( p_a2, a2 )
pointer( p_ax, ax )
pointer( p_xx, xx )
pointer( p_us, us )
pointer( p_vs, vs )
pointer( p_wz, wz )
!jgj
pointer( p_co2_dmmr, co2_dmmr )
pointer( p_co2_dvmr, co2_dvmr )
#include "fv_point.inc"
call fv_stack_push( p_a2, im*(jlast-jfirst+1) )
call fv_stack_push( p_ax, (jlast-jfirst+1)*km )
call fv_stack_push( p_xx, im*(jlast-jfirst+1)*km )
call fv_stack_push( p_us, im*(jlast-jfirst+1)*km )
call fv_stack_push( p_vs, im*(jlast-jfirst+1)*km )
call fv_stack_push( p_wz, im*(jlast-jfirst+1)*(km+1) )
!jgj
call fv_stack_push( p_co2_dmmr, im*(jlast-jfirst+1)*km )
call fv_stack_push( p_co2_dvmr, im*(jlast-jfirst+1)*km )
#endif
fv_time = Time
! fv_time used within fvcore for diagnostics purpose
if ( hs_phys ) then
call get_time (fv_time, seconds, days)
if( print_freq == 0 ) then
prt_minmax = .false.
elseif( print_freq < 0 ) then
prt_minmax = .true.
else
prt_minmax = mod(seconds, 3600*print_freq) == 0
endif
else
call get_date(fv_time, yr, mon, dd, hr, mn, seconds)
if( print_freq == 0 ) then
prt_minmax = .false.
elseif( print_freq < 0 ) then
prt_minmax = .true.
else
prt_minmax = mod(hr, print_freq) == 0 .and. mn==0 .and. seconds==0
endif
endif
if( prt_minmax ) then
if ( master ) then
if ( hs_phys ) then
write(*,*) Days, seconds
else
write(*,*) yr, mon, dd, hr, mn, seconds
endif
endif
call pmaxmin('PS', ps, qmin, qmax, im*(jlast-jfirst+1), 1, 0.01)
#ifdef MARS_GCM
call mars_mass_budget( im, jm, km, jfirst, jlast, ng_d, ps, delp, nq, q, mars_sfc_budg )
#else
! Check dry air mass:
call drymadj(im, jm, km, jfirst, jlast, ng_d, .true., kappa, &
ptop, ps, delp, pe, pk, peln, pkz, nq, &
q, .false.)
#endif MARS_GCM
call pmaxmin('U', ua, qmin, qmax, im*(jlast-jfirst+1), km, 1.)
call pmaxmin('V', va, qmin, qmax, im*(jlast-jfirst+1), km, 1.)
endif
if (id_ps > 0) used = send_data ( id_ps, ps(beglon:endlon,beglat:endlat), Time )
if ( id_ua > 0 .or. id_va > 0 ) then
if (id_ua > 0) used = send_data ( id_ua, ua(beglon:endlon,beglat:endlat,:), Time )
if (id_va > 0) used = send_data ( id_va, va(beglon:endlon,beglat:endlat,:), Time )
endif
! Send vertical velocity (pa/s)
if (id_wa > 0) used = send_data ( id_wa, omga(beglon:endlon,beglat:endlat,:), Time )
! The following sync call is required to ensure that send_data completes before
! array omga is touched by another thread.
if( id_vort>0 .or. id_pv>0 .or. id_ta>0 .or. id_pt>0 ) call fv_array_sync()
if( prt_minmax .and. km>1 ) then
call pmaxmin('W', omga, qmin, qmax, im*(jlast-jfirst+1), km, 1.)
endif
! if ( id_usus>0 .or. id_usvs>0 .or. id_vsts>0 .or. id_tsts>0 ) then
! allocate ( xx(im,jfirst:jlast,km) )
! allocate ( us(im,jfirst:jlast,km) )
! allocate ( vs(im,jfirst:jlast,km) )
! endif
!----------------------------------------------------------------
! Note: wind variances computed from A grid winds are weaker than
! native D grid due to two-grid averaging (from D to A)
!----------------------------------------------------------------
if ( id_usus > 0 .or. id_usvs > 0 ) then
!$omp parallel do default(shared) &
!$omp private(i, j, k)
do k=ksp,kep
call zsmean( ua(1,jfirst,k), im, jfirst, jlast, &
ax(jfirst,k), us(1,jfirst,k), xx(1,jfirst,k) )
enddo
call fv_array_sync()
if ( id_usus > 0 ) used = send_data ( id_usus, xx(beglon:endlon,beglat:endlat,:), Time )
! The following sync call is required to ensure that send_data completes before
! array xx is touched by another thread.
if( id_usvs>0 .or. id_tsts>0 .or. id_vsts>0 ) call fv_array_sync()
endif
if ( id_usvs > 0 .or. id_vsts > 0 ) then
!$omp parallel do default(shared) &
!$omp private(i, j, k)
do k=ksp,kep
call zsmean( va(1,jfirst,k), im, jfirst, jlast, &
ax(jfirst,k), vs(1,jfirst,k) )
if ( id_usvs > 0 ) then
do j=jfirst,jlast
do i=1,im
xx(i,j,k) = us(i,j,k) * vs(i,j,k)
enddo
enddo
endif
enddo
call fv_array_sync()
if ( id_usvs > 0 ) used = send_data ( id_usvs, xx(beglon:endlon,beglat:endlat,:), Time )
! The following sync call is required to ensure that send_data completes before
! array xx is touched by another thread.
if( id_tsts>0 .or. id_vsts>0 ) call fv_array_sync()
endif
if ( id_vort > 0 .or. id_pv > 0 ) then
call vort_d(im, jm, km, jfirst, jlast, u, v, omga, ng_s, ng_d)
if ( id_vort > 0 ) used = send_data ( id_vort, omga(beglon:endlon,beglat:endlat,:), Time )
! The following sync call is required to ensure that send_data completes before
! array omga is touched by another thread.
if( id_pv>0 .or. id_ta>0 .or. id_pt>0 ) call fv_array_sync()
endif
! Potential Vorticity
if ( id_pv > 0 ) then
! Note: this is expensive computation. don't do it too often
call pv_entropy(im, jm, km, jfirst, jlast, omga, pt, pkz, delp, ng_d, grav)
used = send_data ( id_pv, omga(beglon:endlon,beglat:endlat,:), Time )
! The following sync call is required to ensure that send_data completes before
! array omga is touched by another thread.
if( id_ta>0 .or. id_pt>0 ) call fv_array_sync()
endif
if ( id_ta > 0 .or. id_tsts > 0 .or. id_vsts > 0 ) then
!$omp parallel do default(shared) &
!$omp private(i, j, k)
do k=ksp,kep
do j=jfirst,jlast
do i=1,im
omga(i,j,k) = pt(i,j,k)
enddo
enddo
if ( id_tsts > 0 ) then
call zsmean( omga(1,jfirst,k), im, jfirst, jlast, &
ax(jfirst,k), us(1,jfirst,k), xx(1,jfirst,k) )
endif
enddo ! end parallel k-loop
call fv_array_sync()
if ( id_ta > 0) used = send_data ( id_ta, omga(beglon:endlon,beglat:endlat,:), Time )
! The following sync call is required to ensure that send_data completes before
! array omga is touched by another thread.
if( id_pt>0 ) call fv_array_sync()
if ( id_tsts > 0) used = send_data ( id_tsts, xx(beglon:endlon,beglat:endlat,:), Time )
! The following sync call is required to ensure that send_data completes before
! array xx is touched by another thread.
if( id_vsts>0 ) call fv_array_sync()
if( prt_minmax ) then
call pmaxmin('T', omga, qmin, qmax, im*(jlast-jfirst+1), km, 1.)
if(qmax > 360.) call error_mesg('fv_diagnostics:','too warm', FATAL)
#ifdef MARS_GCM
if(qmin < 50.) call error_mesg('fv_diagnostics:','too cold', FATAL)
#else
if(qmin < 100.) call error_mesg('fv_diagnostics:','too cold', FATAL)
#endif
endif
endif
if ( id_pt > 0 ) then
p0 = (1.E5) ** kappa
!$omp parallel do default(shared) &
!$omp private(i, j, k)
do k=ksp,kep
do j=jfirst,jlast
do i=1,im
omga(i,j,k) = pt(i,j,k)/pkz(i,j,k) * p0
enddo
enddo
enddo ! end parallel k-loop
call fv_array_sync()
used = send_data ( id_pt, omga(beglon:endlon,beglat:endlat,:), Time )
call fv_array_sync()
endif
if ( id_vsts > 0 ) then
!$omp parallel do default(shared) &
!$omp private(i, j, k)
do k=ksp,kep
do j=jfirst,jlast
do i=1,im
xx(i,j,k) = vs(i,j,k) * us(i,j,k)
enddo
enddo
enddo
call fv_array_sync()
used = send_data ( id_vsts, xx(beglon:endlon,beglat:endlat,:), Time )
! The following sync call is required to ensure that send_data completes before
! array xx is touched by another thread.
call fv_array_sync()
endif
if ( id_hght > 0 .or. id_h500 > 0 .or. id_slp > 0 .or. &
id_p5km > 0 ) then
#ifdef SW_DYN
gg = ginv
#else
gg = rdgas * ginv
#endif
! allocate ( wz(im,jfirst:jlast,km+1) )
call fv_array_sync() !next loops are parallel over j
!$omp parallel do default(shared) private(i,j,k)
do j=jsp,jep
! Compute height at layer edges
do i=1,im
wz(i,j,km+1) = phis(i,j) * ginv
enddo
do k=km,1,-1
do i=1,im
#ifdef SW_DYN
wz(i,j,k) = wz(i,j,k+1) + gg*(pk(i,j,k+1)-pk(i,j,k))
#else
wz(i,j,k) = wz(i,j,k+1) + gg*pt(i,j,k)*(1.+zvir*q(i,j,k,1)) &
*(peln(i,k+1,j)-peln(i,k,j))
#endif
enddo
enddo
end do
call fv_array_sync()
if(id_hght > 0) used = send_data ( id_hght, wz(beglon:endlon,beglat:endlat,:), Time )
! call pmaxmin('Ztop', wz(1,jfirst,1), qmin, qmax, im, (jlast-jfirst+1), 0.001)
! ztop = gmean(im, jm, jfirst, jlast, wz(1,jfirst,1)) * 0.001
! if( master ) write(6,*) 'Mean model top height (km) =', ztop
! Compute 5-km pressure
if(id_p5km > 0) then
height(1) = 5.E3
call get_pressure_given_height(im, jfirst, jlast, km, wz, 1, height(1), &
pt(1,beglat, nlev), a2, 0.01)
used = send_data ( id_p5km, a2(beglon:endlon,beglat:endlat), Time )
! The following sync call is required to ensure that send_data completes before
! array a2 is touched by another thread.
if( id_slp>0 .or. id_h500>0 ) call fv_array_sync()
endif
! Cumpute SLP (pressure at height=0)
if(id_slp > 0) then
height(1) = 0.
call get_pressure_given_height(im, jfirst, jlast, km, wz, 1, height(1), &
pt(1,beglat, nlev), a2, 0.01)
if( prt_minmax ) &
call pmaxmin('SLP', a2, qmin, qmax, im*(jlast-jfirst+1), 1, 1.)
used = send_data ( id_slp, a2(beglon:endlon,beglat:endlat), Time )
! The following sync call is required to ensure that send_data completes before
! array a2 is touched by another thread.
if( id_h500>0 ) call fv_array_sync()
endif
! Compute H500
if(id_h500 > 0) then
log_p(1) = log( 50000. )
call get_height_given_pressure(im, jfirst, jlast, km, wz, 1, log_p, a2)
used = send_data ( id_h500, a2(beglon:endlon,beglat:endlat), Time )
endif
! deallocate ( wz )
endif
!----------------
! Output tracers:
!----------------
if( nq /= 0 ) then
do ic=1, nq
if (id_tracer(ic) > 0) &
& used = send_data (id_tracer(ic), q(beglon:endlon,beglat:endlat,1:km,ic), Time )
if( prt_minmax ) then
call pmaxming('Q', q(1,jfirst-ng_d,1,ic), im, jm, km, &
jfirst, jlast, ng_d, ng_d, 1.0)
endif
enddo
#ifdef MARS_GCM
#else
!-------
! jgj: per SJ email (jul 17 2008): q_dry = q_moist/(1-sphum)
! mass mixing ratio: q_dry = mass_tracer/mass_dryair = mass_tracer/(mass_air - mass_water) ~ q_moist/(1-sphum)
! co2_mmr = (wco2/wair) * co2_vmr
! Note: There is a check in fv_pack.F90 that ensures that tracer number one is sphum
do ic = nt_phys+1, ncnst
call get_tracer_names ( MODEL_ATMOS, ic, tname, tlongname, tunits )
if (id_tracer_dmmr(ic) > 0 .and. trim(tname).eq.'co2') then
call fv_array_sync() !next loops are parallel over j
!$omp parallel do default(shared) private(i,j,k)
do k=1,km
do j=jsp,jep
do i=1,im
co2_dmmr(i,j,k) = q(i,j,k,ic)/(1.0-q(i,j,k,1))
enddo
enddo
enddo
call fv_array_sync()
used = send_data (id_tracer_dmmr(ic), co2_dmmr(beglon:endlon,beglat:endlat,1:km), Time )
! The following sync call is required to ensure that send_data completes before
! array co2_dmmr is touched by another thread.
call fv_array_sync()
if( prt_minmax ) then
call pmaxming('co2_dmmr', co2_dmmr(1,jfirst-ng_d,1), im, jm, km, &
jfirst, jlast, ng_d, ng_d, 1.0)
endif
endif
!! use dry vmr to compare to obs
if (id_tracer_dvmr(ic) > 0 .and. trim(tname).eq.'co2') then
call fv_array_sync() !next loops are parallel over j
!$omp parallel do default(shared) private(i,j,k)
do k=1,km
do j=jsp,jep
do i=1,im
co2_dvmr(i,j,k) = (q(i,j,k,ic)/(1.0-q(i,j,k,1))) * WTMAIR/WTMCO2 ! to compare to obs
enddo
enddo
enddo
call fv_array_sync()
used = send_data (id_tracer_dvmr(ic), co2_dvmr(beglon:endlon,beglat:endlat,1:km), Time )
! The following sync call is required to ensure that send_data completes before
! array co2_dmmr is touched by another thread.
call fv_array_sync()
if( prt_minmax ) then
call pmaxming('co2_dvmr', co2_dvmr(1,jfirst-ng_d,1), im, jm, km, &
jfirst, jlast, ng_d, ng_d, 1.0)
endif
endif
enddo
#endif
!-------
! Check tracer mass. Unit:
! if( prt_minmax ) &
! call tracer_mass(im, jm, km, jfirst, jlast, ng_d, nq, q, delp)
if ( age_tracer ) then
! Assuming age tracer is the last tracer
call age_of_air(im, jm, km, jfirst, jlast, ng_d, age_time, &
dt_atmos, peln, delp, q(1,jfirst-ng_d,1,nq), &
prt_minmax)
endif
endif
! if ( id_usus>0 .or. id_usvs>0 .or. id_vsts>0 .or. id_tsts>0 ) then
! deallocate ( xx )
! deallocate ( us )
! deallocate ( vs )
! endif
! make sure all send_data completed before one thread exits the routine
call fv_array_sync()
end subroutine fv_diag
subroutine get_pressure_given_height(im, jfirst, jlast, km, wz, kd, height, &
ts, a2, fac)
#include "fv_arrays.h"
integer, intent(in):: im, jfirst, jlast, km
integer, intent(in):: kd ! vertical dimension of the ouput height
real, intent(in):: wz(im,jfirst:jlast,km+1)
real, intent(in):: ts(im,jfirst:jlast)
real, intent(in):: height(kd) ! must be monotonically decreasing with increasing k
real, optional, intent(in):: fac
real, intent(out):: a2(im,jfirst:jlast,kd) ! pressure (pa)
! local:
integer n, i,j,k
real ptmp, tm
integer ii, mcount
#include "fv_point.inc"
call fv_array_check( LOC(wz) )
call fv_array_check( LOC(ts) )
call fv_array_check( LOC(a2) )
do n=1,kd
!$omp parallel do private(ii,i,j,k, ptmp, tm, mcount)
do j=jsp,jep
do 1000 i=1,im
if ( height(n) >= wz(i,j,km+1) ) then
!---------------------
! Search from top down
!---------------------
do k=1,km
if( height(n) < wz(i,j,k) .and. height(n) >= wz(i,j,k+1) ) then
! Found it!
ptmp = peln(i,k,j) + (peln(i,k+1,j)-peln(i,k,j)) * &
(wz(i,j,k)-height(n)) / (wz(i,j,k)-wz(i,j,k+1))
a2(i,j,n) = exp(ptmp)
go to 500
endif
enddo
else
! Method: use a mean tm based on zonal mean with height < local height
tm = 0.
mcount = 0
do ii=1,im
if (wz(ii,j,km+1) < wz(i,j,km+1) .and. wz(ii,j,km+1) > (height(n)+0.1)) then
! adding 0.1 to exclude most ocean points
! This algorithm works very well except the S pole.
mcount = mcount + 1
tm = tm + ts(ii,j)
endif
enddo
if ( mcount > im/10 + 1 ) then
tm = rdgas*ginv * tm / real(mcount)
else
! Assuming 6.5 deg/km lapse rate
tm = rdgas*ginv*(ts(i,j) + 3.25E-3*(wz(i,j,km)-height(n)))
endif
a2(i,j,n) = exp( peln(i,km+1,j) + (wz(i,j,km+1) - height(n))/tm )
endif
500 if ( present(fac) ) a2(i,j,n) = fac * a2(i,j,n)
1000 continue
enddo
enddo
call fv_array_sync()
end subroutine get_pressure_given_height
subroutine get_height_given_pressure(im, jfirst, jlast, km, wz, kd, log_p, a2)
#include "fv_arrays.h"
integer, intent(in):: im, jfirst, jlast, km
integer, intent(in):: kd ! vertical dimension of the ouput height
real, intent(in):: log_p(kd) ! must be monotonically decreasing with increasing k
! log (p)
real, intent(in):: wz(im,jfirst:jlast,km+1)
real, intent(out):: a2(im,jfirst:jlast,kd) ! height (m)
! local:
integer n, i,j,k
#include "fv_point.inc"
call fv_array_check( LOC(wz) )
call fv_array_check( LOC(a2) )
do n=1,kd
!$omp parallel do private(i,j,k)
do j=jsp,jep
do 1000 i=1,im
do k=1,km
if( log_p(n) <= peln(i,k+1,j) .and. log_p(n) >= peln(i,k,j) ) then
! Found it!
a2(i,j,n) = wz(i,j,k) + (wz(i,j,k+1) - wz(i,j,k)) * &
(log_p(n)-peln(i,k,j)) / (peln(i,k+1,j)-peln(i,k,j) )
go to 1000
endif
enddo
a2(i,j,n) = missing_value
1000 continue
enddo
enddo
call fv_array_sync()
end subroutine get_height_given_pressure
subroutine zsmean(a2, im, js, je, ux, us, usus)
integer im, js, je
real, intent(in) :: a2(im,js:je)
real, intent(out):: ux(js:je) ! zonal mean
real, intent(out):: us(im,js:je) ! departure from ux
real, optional, intent(out):: usus(im,js:je) ! us ** 2
! local
integer i, j
real rim
rim = 1. / im
do j=js,je
ux(j) = a2(1,j)
do i=2,im
ux(j) = ux(j) + a2(i,j)
enddo
ux(j) = ux(j) * rim
do i=1,im
us(i,j) = a2(i,j) - ux(j)
enddo
if ( present(usus) ) then
do i=1,im
usus(i,j) = us(i,j) ** 2
enddo
endif
enddo
end subroutine zsmean
subroutine age_of_air(im, jm, km, jfirst, jlast, ng, atime, dt, &
peln, delp, q, print_max)
use fv_arrays_mod, only: fv_array_check, isp, iep, jsp, jep, ksp, kep
integer im
integer jm
integer km
integer jfirst
integer jlast
integer ng
logical, intent(in):: print_max
! q is the age tracer
! Need to be converted to mixing ratio (mass of tracer / dry_air-mass)
! Ignore this inconsistency for now.
real, intent(in):: delp(im,jfirst:jlast,km)
real, intent(in):: peln(im,km+1,jfirst:jlast)
real, intent(inout):: atime ! accumulated time since init
real, intent(in):: dt
real, intent(inout):: q(im,jfirst-ng:jlast+ng,km)
! Local
integer i, j, k
real p_source ! source level (pa)
real ascale, ra, ryear
real tiny, pfull
real rim, qsum, qmin, qmax
real age(im,jfirst:jlast,km)
parameter ( tiny = 1.e-6 )
parameter ( p_source = 75000. )
parameter ( ascale = 5.e-6 / 60. )
parameter ( ra = 1./ascale )
parameter (ryear = 1./(365.*24.*3600) )
rim = 1./float(im)
atime = atime + dt
call fv_array_check( LOC(delp) )
call fv_array_check( LOC(peln) )
call fv_array_check( LOC(q) )
!$omp parallel do private(i, j, k, qsum, pfull)
do k=ksp,kep
do j=jfirst, jlast
do i=1,im
pfull = delp(i,j,k) / (peln(i,k+1,j)-peln(i,k,j))
if( atime < tiny ) then
q(i,j,k) = 0.
elseif( pfull >= p_source ) then
q(i,j,k) = ascale * atime
endif
enddo
if(id_aoa > 0 ) then
do i=1,im
age(i,j,k) = max(0., (atime - ra*q(i,j,k))) * ryear
enddo
endif
enddo
enddo
call fv_array_sync()
if(id_aoa > 0 ) then
used = send_data (id_aoa, age(beglon:endlon,beglat:endlat,:), fv_time)
if(print_max ) &
call pmaxmin('Age-of-Air (yr)', age, qmin, qmax, im*(jlast-jfirst+1), km, 1.)
endif
! make sure all send_data completed before one thread exits the routine
call fv_array_sync()
end subroutine age_of_air
end module fv_diagnostics