!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! !! !! 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 bgrid_diagnostics_mod !----------------------------------------------------------------------- use bgrid_horiz_mod, only: horiz_grid_type use bgrid_vert_mod, only: vert_grid_type, compute_pres_full, & compute_pres_half, compute_pres_depth use bgrid_masks_mod, only: grid_mask_type use bgrid_prog_var_mod, only: prog_var_type use bgrid_change_grid_mod, only: change_grid, TEMP_GRID, WIND_GRID use diag_manager_mod, only: diag_axis_init, register_diag_field, & register_static_field, send_data use time_manager_mod, only: time_type use fms_mod, only: file_exist, open_namelist_file, & error_mesg, NOTE, check_nml_error, & mpp_pe, mpp_root_pe, stdlog, & close_file, write_version_number use constants_mod, only: GRAV, KAPPA, RDGAS use field_manager_mod, only: MODEL_ATMOS use tracer_manager_mod, only: get_tracer_names, get_number_tracers implicit none private public :: bgrid_diagnostics_init, & bgrid_diagnostics, & bgrid_diagnostics_tend !----------------------------------------------------------------------- !------------------------- axis names ---------------------------------- character(len=8) :: axiset = 'dynamics' character(len=8) :: mod_name = 'dynamics' !----------------------------------------------------------------------- integer, parameter :: MXTR = 10 real, parameter :: GINV = 1./GRAV !----------------------------------------------------------------------- ! axis and field identifiers for the diag manager integer :: id_hlonb, id_hlon , id_hlatb, id_hlat , & id_vlonb, id_vlon , id_vlatb, id_vlat , & id_phalf, id_pfull, id_hlat_wgt, id_vlat_wgt integer :: id_bk , id_pk , id_zsurf, id_res , id_wspd, & id_ps , id_ucomp, id_vcomp, id_temp , id_pres_full, & id_omega, id_div , id_vor , id_pgfx , id_pgfy, & id_udt , id_vdt , id_tdt , id_pres_half, & id_alm , id_aph , id_theta, id_mfew, id_mfns, id_slp integer, allocatable :: id_tracer(:), id_tracer_tend(:) integer :: id_ucomp_sq, id_vcomp_sq, id_temp_sq, id_omega_sq, & id_ucomp_vcomp, id_omega_temp !----------------------------------------------------------------------- ! need to save surface geopotential height argument to initialization call ! the surface height will be needed for computing sea level pressure real, pointer, dimension(:,:) :: zsurfg !----------------------------------------------------------------------- character(len=128) :: version = '$Id: bgrid_diagnostics.F90,v 11.0 2004/09/28 19:07:18 fms Exp $' character(len=128) :: tag = '$Name: mom4p1_pubrel_dec2009_nnz $' !----------------------------------------------------------------------- contains !####################################################################### subroutine bgrid_diagnostics_init ( Time, Hgrid, Vgrid, Var, & fis, res, & mass_axes, vel_axes ) !----------------------------------------------------------------------- ! setup/write netcdf metadata and static fields !----------------------------------------------------------------------- ! Time = current/initial time ! Hgrid = horizontal grid constants ! Vgrid = vertical grid constants ! fis = geopotential height of the surface ! res = reciprocal of eta at the surface ! mass_axes = axis identifiers for the temperature (mass) grid ! vel_axes = axis identifiers for the velocity grid ! NOTE: The axes identifiers are for the lon, lat, pfull, and ! phalf axes, respectively. They are returned by the diag_manager. !----------------------------------------------------------------------- type(time_type), intent(in) :: Time type(horiz_grid_type), intent(in) :: Hgrid type (vert_grid_type), intent(in) :: Vgrid type (prog_var_type), intent(in) :: Var real, intent(in), dimension(Hgrid%ilb:,Hgrid%jlb:), target :: fis, res integer, dimension(4), intent(out) :: mass_axes, vel_axes !----------------------------------------------------------------------- real, dimension(Hgrid%Tmp%isg:Hgrid%Tmp%ieg+1) :: hlonb real, dimension(Hgrid%Vel%isg:Hgrid%Vel%ieg+1) :: vlonb real, dimension(Hgrid%Tmp%jsg:Hgrid%Tmp%jeg+1) :: hlatb real, dimension(Hgrid%Vel%jsg:Hgrid%Vel%jeg+1) :: vlatb real, dimension(Hgrid%Tmp%isg:Hgrid%Tmp%ieg) :: hlon real, dimension(Hgrid%Vel%isg:Hgrid%Vel%ieg) :: vlon real, dimension(Hgrid%Tmp%jsg:Hgrid%Tmp%jeg) :: hlat real, dimension(Hgrid%Vel%jsg:Hgrid%Vel%jeg) :: vlat real, dimension(Hgrid%Tmp%js:Hgrid%Tmp%je) :: hlat_wgt real, dimension(Hgrid%Vel%js:Hgrid%Vel%je) :: vlat_wgt real, dimension(1,1) :: psurf real, dimension(1,1,Vgrid%nlev) :: pfull real, dimension(1,1,Vgrid%nlev+1) :: phalf real :: vrange(2), trange(2), prange(2) real :: rad2deg integer :: i, j, n, unit, io, ierr, ntprog integer :: isg, ieg, hsg, heg, vsg, veg integer :: is, ie, hs, he, vs, ve integer :: uflx_axes(4), vflx_axes(4) logical :: used character(len=128) :: tname character(len=256) :: longname, units !--------------------------- set up axes ------------------------------- ! compute grid indices is = Hgrid % Tmp % is; ie = Hgrid % Tmp % ie hs = Hgrid % Tmp % js; he = Hgrid % Tmp % je vs = Hgrid % Vel % js; ve = Hgrid % Vel % je ! global grid indices isg = Hgrid % Tmp % isg; ieg = Hgrid % Tmp % ieg hsg = Hgrid % Tmp % jsg; heg = Hgrid % Tmp % jeg vsg = Hgrid % Vel % jsg; veg = Hgrid % Vel % jeg ! grid box boundaries in degrees rad2deg = 90./acos(0.0) hlonb(isg:ieg+1) = Hgrid % Tmp % blong(isg:ieg+1) * rad2deg hlatb(hsg:heg+1) = Hgrid % Tmp % blatg(hsg:heg+1) * rad2deg vlonb(isg:ieg+1) = Hgrid % Vel % blong(isg:ieg+1) * rad2deg vlatb(vsg:veg+1) = Hgrid % Vel % blatg(vsg:veg+1) * rad2deg ! grid box centers in degrees do i = isg, ieg hlon(i) = 0.5*(hlonb(i)+hlonb(i+1)) vlon(i) = 0.5*(vlonb(i)+vlonb(i+1)) enddo do j = hsg, heg hlat(j) = 0.5*(hlatb(j)+hlatb(j+1)) enddo do j = vsg, veg vlat(j) = 0.5*(vlatb(j)+vlatb(j+1)) enddo ! compute a reference profile of pressure based on psurf = 1000 hPa psurf = reshape ( (/ 100000. /), (/ 1, 1 /) ) call compute_pres_full (Vgrid, psurf, pfull) call compute_pres_half (Vgrid, psurf, phalf) ! in units of hPa pfull = pfull*0.01 phalf = phalf*0.01 !----- initialize mass axes ------ id_hlonb = diag_axis_init ( 'lonb', hlonb, 'degrees_E', 'x', & 'longitude edges', set_name='atmos', & Domain2=Hgrid%Tmp%Domain_nohalo ) id_hlon = diag_axis_init ( 'lon', hlon, 'degrees_E', 'x', & 'longitude', set_name='atmos', & edges=id_hlonb, & Domain2=Hgrid%Tmp%Domain_nohalo ) id_hlatb = diag_axis_init ( 'latb', hlatb, 'degrees_N', 'y', & 'latitude edges', set_name='atmos', & Domain2=Hgrid%Tmp%Domain_nohalo ) id_hlat = diag_axis_init ( 'lat', hlat, 'degrees_N', 'y', & 'latitude', set_name='atmos', & edges=id_hlatb, & Domain2=Hgrid%Tmp%Domain_nohalo ) !----- initialize velocity axes ------ id_vlonb = diag_axis_init ( 'vlonb', vlonb, 'degrees_E', 'x', & 'longitude edges', set_name='atmos', & Domain2=Hgrid%Vel%Domain_nohalo ) id_vlon = diag_axis_init ( 'vlon', vlon, 'degrees_E', 'x', & 'longitude', set_name='atmos', & edges=id_vlonb, & Domain2=Hgrid%Vel%Domain_nohalo ) id_vlatb = diag_axis_init ( 'vlatb', vlatb, 'degrees_N', 'y', & 'latitude edges', set_name='atmos', & Domain2=Hgrid%Vel%Domain_nohalo ) id_vlat = diag_axis_init ( 'vlat', vlat, 'degrees_N', 'y', & 'latitude', set_name='atmos', & edges=id_vlatb, & Domain2=Hgrid%Vel%Domain_nohalo ) !----- initialize vertical axes ----- id_phalf = diag_axis_init ( 'phalf', phalf(1,1,:), 'hPa', 'z', & 'approx half pressure level', & direction=-1, set_name='atmos' ) id_pfull = diag_axis_init ( 'pfull', pfull(1,1,:), 'hPa', 'z', & 'approx full pressure level', & direction=-1, edges=id_phalf, & set_name='atmos' ) !----------------------------------------------------------------------- !-------- initialize and output variables with no time axis ------------ mass_axes = (/ id_hlon, id_hlat, id_pfull, id_phalf /) vel_axes = (/ id_vlon, id_vlat, id_pfull, id_phalf /) uflx_axes = (/ id_vlon, id_hlat, id_pfull, id_phalf /) vflx_axes = (/ id_hlon, id_vlat, id_pfull, id_phalf /) ! valid range for some fields vrange = (/ -400., +400. /) ! momentum trange = (/ 100., 400. /) ! temperature prange = (/ -1., 107500. /) ! pressure !----------------------------------------------------------------------- !---- 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/), & 'vertical coordinate reference pressure value (ak*pref)', 'pascals' ) id_zsurf = register_static_field ( mod_name, 'zsurf', mass_axes(1:2),& 'surface height', 'm' ) id_res = register_static_field ( mod_name, 'res', mass_axes(1:2), & 'reciprocal of sigma/eta at the surface', 'none' ) id_alm = register_static_field ( mod_name, 'alm', mass_axes(1:2), & 'actual longitudes for temperature grid', 'degrees_E' ) id_aph = register_static_field ( mod_name, 'aph', mass_axes(1:2), & 'actual latitudes for temperature grid', 'degrees_N' ) ! these changes cannot be implemented until changes to diag_manager ! initialize fields useful for computing offline global averages ! ! id_hlat_wgt = register_static_field ( mod_name, 'lat_wgt', & ! (/id_hlat/), 'latitude weight for mass grid', 'none' ) ! ! id_vlat_wgt = register_static_field ( mod_name, 'vlat_wgt', & ! (/id_vlat/), 'latitude weight for momentum grid', 'none' ) if ( id_bk > 0 ) & used = send_data ( id_bk, Vgrid%eta, Time ) if ( id_pk > 0 ) & used = send_data ( id_pk, Vgrid%peta, Time ) if ( id_zsurf > 0 ) & used = send_data ( id_zsurf, fis(is:ie,hs:he)*GINV, Time ) if ( id_res > 0 ) & used = send_data ( id_res, res(is:ie,hs:he), Time ) if ( id_alm > 0 ) & used = send_data ( id_alm, Hgrid%Tmp%alm(is:ie,hs:he)*rad2deg, Time ) if ( id_aph > 0 ) & used = send_data ( id_aph, Hgrid%Tmp%aph(is:ie,hs:he)*rad2deg, Time ) ! if ( id_hlat_wgt > 0 ) then ! hlat_wgt = sin(Hgrid%Tmp%blatg(hs+1:he+1))-sin(Hgrid%Tmp%blatg(hs:he)) ! used = send_data ( id_hlat_wgt, hlat_wgt, Time ) ! endif ! ! if ( id_vlat_wgt > 0 ) then ! vlat_wgt = sin(Hgrid%Vel%blatg(vs+1:ve+1))-sin(Hgrid%Vel%blatg(vs:ve)) ! used = send_data ( id_vlat_wgt, vlat_wgt, Time ) ! endif !---- register non-static fields ------- id_ps = register_diag_field ( mod_name, 'ps', mass_axes(1:2), & Time, 'surface pressure', 'pascals' ) id_slp = register_diag_field ( mod_name, 'slp', mass_axes(1:2), & Time, 'sea level pressure', 'pascals' ) id_ucomp = register_diag_field ( mod_name, 'ucomp', vel_axes(1:3), & Time, 'zonal wind component', 'm/sec', & missing_value=vrange(1), range=vrange ) id_vcomp = register_diag_field ( mod_name, 'vcomp', vel_axes(1:3), & Time, 'meridional wind component', 'm/sec', & missing_value=vrange(1), range=vrange ) id_temp = register_diag_field ( mod_name, 'temp', mass_axes(1:3), & Time, 'temperature', 'deg_k', & missing_value=trange(1), range=trange ) id_pres_full = register_diag_field ( mod_name, 'pres_full', mass_axes(1:3), & Time, 'pressure at full model levels', 'pascals', & missing_value=prange(1), range=prange ) ! pressure at half levels id_pres_half = register_diag_field ( mod_name, 'pres_half', & (/ id_hlon, id_hlat, id_phalf /), Time, & 'pressure at half model levels', 'pascals', & missing_value=prange(1), range=prange ) id_omega = register_diag_field ( mod_name, 'omega', mass_axes(1:3),& Time, 'omega vertical velocity', & 'pascals/sec', & missing_value=-999. ) id_theta = register_diag_field ( mod_name, 'theta', mass_axes(1:3), & Time, 'potential temperature', 'deg_k', & missing_value=-999. ) id_mfew = register_diag_field ( mod_name, 'mfew', uflx_axes(1:3), & Time, 'Zonal mass flux', 'Pa-m2/s', & missing_value=-1.e30 ) id_mfns = register_diag_field ( mod_name, 'mfns', vflx_axes(1:3), & Time, 'Meridional mass flux', 'Pa-m2/s', & missing_value=-1.e30 ) ! write version (to log file) call write_version_number (version,tag) ! register diagnostics for all tracers allocate (id_tracer(Var%ntrace)) if (mpp_pe() == mpp_root_pe()) write(stdlog(),100) trim(mod_name) do n = 1, Var%ntrace call get_tracer_names ( MODEL_ATMOS, n, tname, longname, units ) if (mpp_pe() == mpp_root_pe()) write(stdlog(),110) trim(tname),trim(longname),trim(units) id_tracer(n) = register_diag_field ( mod_name, trim(tname), & mass_axes(1:3), Time, trim(longname), & trim(units), missing_value=-999. ) enddo 100 format ('Diagnostics for the following tracer fields are available for module name = ',a) 110 format (3x,a,' (',a,'; ',a,')') !-------- register second-moment quantities ------- ! (for now we are only saving fields on the same grids) id_ucomp_sq = register_diag_field ( mod_name, 'ucomp_sq', vel_axes(1:3), & Time, 'zonal wind component squared', 'm2/s2', & missing_value=-1., range=(/0.,vrange(2)**2/) ) id_vcomp_sq = register_diag_field ( mod_name, 'vcomp_sq', vel_axes(1:3), & Time, 'meridional wind component squared', 'm2/s2', & missing_value=-1., range=(/0.,vrange(2)**2/) ) id_temp_sq = register_diag_field ( mod_name, 'temp_sq', mass_axes(1:3), & Time, 'temperature squared', 'deg_K**2', & missing_value=-1., range=(/0.,trange(2)**2/) ) id_omega_sq = register_diag_field ( mod_name, 'omega_sq', mass_axes(1:3),& Time, 'omega vertical velocity squared', & 'Pa**2/s**2', missing_value=-999. ) id_ucomp_vcomp = register_diag_field ( mod_name, 'ucomp_vcomp', vel_axes(1:3),& Time, 'zonal times meridional wind components', 'm2/s2', & missing_value=-1. ) id_omega_temp = register_diag_field ( mod_name, 'omega_temp', mass_axes(1:3),& Time, 'omega vertical velocity time temperature',& 'Pascals*deg_K/sec', missing_value=-999. ) !-------- wind speed, divergence, vorticity ---------------------------- id_wspd = register_diag_field ( mod_name, 'wspd', vel_axes(1:3), & Time, 'wind speed', 'm/s', missing_value=-999.,& range=(/0.,vrange(2)/) ) id_div = register_diag_field ( mod_name, 'div', mass_axes(1:3), & Time, 'divergence', '1/s', missing_value=-999. ) id_vor = register_diag_field ( mod_name, 'vor', mass_axes(1:3), & Time, 'relative vorticity', '1/s', missing_value=-999. ) !--------- pressure gradient components (NOT USED) --------------------- ! id_pgfx = register_diag_field ( mod_name, 'pgfx', vel_axes(1:3), & ! Time, 'zonal pressure gradient force', & ! 'm/s2', missing_value=-999. ) ! id_pgfy = register_diag_field ( mod_name, 'pgfy', vel_axes(1:3), & ! Time, 'meridional pressure gradient force', & ! 'm/s2', missing_value=-999. ) !----------------------------------------------------------------------- ! -------- tendencies --------- id_udt = register_diag_field ( mod_name, 'udt_dyn', vel_axes(1:3), & Time, 'zonal wind tendency for dynamics', & 'm/s2', missing_value=-999. ) id_vdt = register_diag_field ( mod_name, 'vdt_dyn', vel_axes(1:3), & Time, 'meridional wind tendency for dynamics', & 'm/s2', missing_value=-999. ) id_tdt = register_diag_field ( mod_name, 'tdt_dyn', mass_axes(1:3), & Time, 'temperature tendency for dynamics', & 'deg_k/sec', missing_value=-999. ) ! tendencies for prognostic tracers only call get_number_tracers ( MODEL_ATMOS, num_prog=ntprog ) allocate (id_tracer_tend(ntprog)) do n = 1, ntprog call get_tracer_names ( MODEL_ATMOS, n, tname, longname, units ) tname = trim(tname) //'_dt_dyn' longname = trim(longname)//' tendency for dynamics' units = trim(units) //'/s' if (units == 'none') units = '1/sec' if (mpp_pe() == mpp_root_pe()) write(stdlog(),110) trim(tname),trim(longname),trim(units) id_tracer_tend(n) = register_diag_field ( mod_name, trim(tname), & mass_axes(1:3), Time, trim(longname), & trim(units), missing_value=-999. ) enddo ! save surface geopotential height for computing sea level pressure zsurfg => fis !----------------------------------------------------------------------- end subroutine bgrid_diagnostics_init !####################################################################### subroutine bgrid_diagnostics ( Hgrid, Vgrid, Var, Masks, Time, & omega, div, mfew, mfns ) !----------------------------------------------------------------------- ! write netcdf fields !----------------------------------------------------------------------- ! Hgrid = horizontal grid constants ! Vgrid = vertical grid constants ! Var = prognostic variables at diagnostics Time ! Masks = grid box masks for step-mountain topography ! Time = diagnostics time ! omega = omega (vertical velocity) diagnostic !----------------------------------------------------------------------- type(horiz_grid_type), intent(in) :: Hgrid type (vert_grid_type), intent(in) :: Vgrid type (prog_var_type), intent(in) :: Var type(grid_mask_type), intent(in) :: Masks type(time_type), intent(in) :: Time real, intent(in), dimension(Hgrid%ilb:,Hgrid%jlb:,:) :: & omega, div, mfew, mfns ! real, intent(in), dimension(Hgrid%ilb:,Hgrid%jlb:,:), optional ::& ! div, pgfx, pgfy real, dimension(Hgrid%ilb:Hgrid%iub,Hgrid%jlb:Hgrid%jub) :: slp real, dimension(Hgrid%ilb:Hgrid%iub,Hgrid%jlb:Hgrid%jub, & Vgrid%nlev) :: wspd, vor, dp, udp, vdp real, dimension(Hgrid%ilb:Hgrid%iub,Hgrid%jlb:Hgrid%jub, & Vgrid%nlev+1) :: ph logical, dimension(Hgrid%ilb:Hgrid%iub,Hgrid%jlb:Hgrid%jub, & Vgrid%nlev+1) :: lmask !----------------------------------------------------------------------- integer :: is, ie, hs, he, vs, ve, n, j, k logical :: used !----------------------------------------------------------------------- is = Hgrid % Tmp % is; ie = Hgrid % Tmp % ie hs = Hgrid % Tmp % js; he = Hgrid % Tmp % je vs = Hgrid % Vel % js; ve = Hgrid % Vel % je !----------------------------------------------------------------------- !---------------- surface fields --------------------------------------- if ( id_ps > 0 ) & used = send_data ( id_ps , Var%ps(is:ie,hs:he), Time ) !---------------- 3d momentum fields (u & v) --------------------------- if ( id_ucomp > 0 ) & used = send_data ( id_ucomp, Var%u(is:ie,vs:ve,:), Time, & mask=Masks%Vel%mask(is:ie,vs:ve,:) > 0.5 ) if ( id_vcomp > 0 ) & used = send_data ( id_vcomp, Var%v(is:ie,vs:ve,:), Time, & mask=Masks%Vel%mask(is:ie,vs:ve,:) > 0.5 ) if ( id_temp > 0 ) & used = send_data ( id_temp, Var%t(is:ie,hs:he,:), Time, & mask=Masks%Tmp%mask(is:ie,hs:he,:) > 0.5 ) do n = 1, Var%ntrace if ( id_tracer(n) > 0 ) & used = send_data ( id_tracer(n), Var%r(is:ie,hs:he,:,n), Time, & mask=Masks%Tmp%mask(is:ie,hs:he,:) > 0.5 ) enddo if ( id_omega > 0 ) & used = send_data ( id_omega, omega(is:ie,hs:he,:), Time, & mask=Masks%Tmp%mask(is:ie,hs:he,:) > 0.5 ) ! pressure at full levels ! Note: not computational efficient to recompute pfull if ( id_pres_full > 0 ) then call compute_pres_full (Vgrid, Var%pssl, dp) used = send_data ( id_pres_full, dp(is:ie,hs:he,:), Time, & mask=Masks%Tmp%mask(is:ie,hs:he,:) > 0.5 ) endif ! pressure at half (interface) levels ! Note: not computational efficient to recompute phalf if ( id_pres_half > 0 ) then call compute_pres_half (Vgrid, Var%pssl, ph) lmask(is:ie,hs:he,1) = .true. lmask(is:ie,hs:he,2:Vgrid%nlev+1) = Masks%Tmp%mask(is:ie,hs:he,:) > 0.5 used = send_data ( id_pres_half, ph(is:ie,hs:he,:), Time, & mask=lmask(is:ie,hs:he,:) ) endif !--- sea level pressure --- if ( id_slp > 0 ) then if ( id_pres_full <= 0 ) call compute_pres_full (Vgrid, Var%pssl, dp) call sea_level_pressure ( Var%ps, zsurfg, dp, Var%t, slp ) used = send_data ( id_slp, slp(is:ie,hs:he), Time ) endif ! potential temperature (compute pfull if necessary) if ( id_theta > 0 ) then if ( id_pres_full <= 0 .and. id_slp <= 0 ) call compute_pres_full (Vgrid, Var%pssl, dp) dp = Var%t * (1000.e2/dp)**KAPPA used = send_data ( id_theta, dp(is:ie,hs:he,:), Time, & mask=Masks%Tmp%mask(is:ie,hs:he,:) > 0.5 ) endif !--------- second moment quantities ---------- if ( id_ucomp_sq > 0 ) & used = send_data ( id_ucomp_sq, Var%u(is:ie,vs:ve,:)**2, Time, & mask=Masks%Vel%mask(is:ie,vs:ve,:) > 0.5 ) if ( id_vcomp_sq > 0 ) & used = send_data ( id_vcomp_sq, Var%v(is:ie,vs:ve,:)**2, Time, & mask=Masks%Vel%mask(is:ie,vs:ve,:) > 0.5 ) if ( id_temp_sq > 0 ) & used = send_data ( id_temp_sq, Var%t(is:ie,hs:he,:)**2, Time, & mask=Masks%Tmp%mask(is:ie,hs:he,:) > 0.5 ) if ( id_omega_sq > 0 ) & used = send_data ( id_omega_sq, omega(is:ie,hs:he,:)**2, Time, & mask=Masks%Tmp%mask(is:ie,hs:he,:) > 0.5 ) if ( id_ucomp_vcomp > 0 ) used = send_data ( id_ucomp_vcomp, & Var%u(is:ie,vs:ve,:)*Var%v(is:ie,vs:ve,:), Time, & mask=Masks%Vel%mask(is:ie,vs:ve,:) > 0.5 ) if ( id_omega_temp > 0 ) used = send_data ( id_omega_temp, & omega(is:ie,hs:he,:)*Var%t(is:ie,hs:he,:), Time, & mask=Masks%Tmp%mask(is:ie,hs:he,:) > 0.5 ) !------------ wind speed, divergence, vorticity ------------------------ if ( id_wspd > 0 ) then wspd(is:ie,vs:ve,:) = sqrt & ( Var%u(is:ie,vs:ve,:)*Var%u(is:ie,vs:ve,:) + & Var%v(is:ie,vs:ve,:)*Var%v(is:ie,vs:ve,:) ) used = send_data ( id_wspd, wspd(is:ie,vs:ve,:), Time, & mask=Masks%Vel%mask(is:ie,vs:ve,:) > 0.5 ) endif if ( id_div > 0 ) then used = send_data ( id_div, div(is:ie,hs:he,:), Time, & mask=Masks%Tmp%mask(is:ie,hs:he,:) > 0.5 ) endif if ( id_vor > 0 ) then !if ( id_vor > 0 .or. id_div > 0 ) then !--- precompute quantities common to both vor and div --- call compute_pres_depth (Vgrid, Var%pssl, dp) call change_grid (Hgrid, TEMP_GRID, WIND_GRID, dp, udp) vdp = Var%v * udp ! note: using udp to store dp at vel pts udp = Var%u * udp !if ( id_vor > 0 ) then call compute_vorticity (Hgrid, dp, udp, vdp, vor ) used = send_data ( id_vor, vor(is:ie,hs:he,:), Time, & mask=Masks%Tmp%mask(is:ie,hs:he,:) > 0.5 ) !endif !if ( id_div > 0 ) then ! call compute_divergence (Hgrid, dp, udp, vdp, div ) ! used = send_data ( id_div, div(is:ie,hs:he,:), Time, & ! mask=Masks%Tmp%mask(is:ie,hs:he,:) > 0.5 ) !endif endif !------- mass fluxes (without topography masks) ----------- if ( id_mfew > 0 ) then used = send_data ( id_mfew, mfew(is:ie,hs:he,:), Time ) endif if ( id_mfns > 0 ) then used = send_data ( id_mfns, mfns(is:ie,vs:ve,:), Time ) endif !--------------- pressure gradient components -------------------------- !------------------------ NOT USED ------------------------------------- ! if ( id_pgfx > 0 .and. present(pgfx) ) & ! used = send_data ( id_pgfx, pgfx(is:ie,vs:ve,:), Time, & ! mask=Masks%Vel%mask(is:ie,vs:ve,:) > 0.5 ) ! if ( id_pgfy > 0 .and. present(pgfy) ) & ! used = send_data ( id_pgfy, pgfy(is:ie,vs:ve,:), Time, & ! mask=Masks%Vel%mask(is:ie,vs:ve,:) > 0.5 ) !----------------------------------------------------------------------- end subroutine bgrid_diagnostics !####################################################################### subroutine bgrid_diagnostics_tend ( Hgrid, Var_dt, Masks, Time ) !----------------------------------------------------------------------- ! Hgrid = horizontal grid constants ! Var_dt = prognostic variables tendencies FROM ONLY THE DYNAMICS ! Masks = grid box masks for step-mountain topography ! Time = diagnostics time !----------------------------------------------------------------------- type(horiz_grid_type), intent(in) :: Hgrid type (prog_var_type), intent(in) :: Var_dt type(grid_mask_type), intent(in) :: Masks type(time_type), intent(in) :: Time !----------------------------------------------------------------------- integer :: is, ie, hs, he, vs, ve, n logical :: used !----------------------------------------------------------------------- ! compute domain indices is = Hgrid % Tmp % is; ie = Hgrid % Tmp % ie hs = Hgrid % Tmp % js; he = Hgrid % Tmp % je vs = Hgrid % Vel % js; ve = Hgrid % Vel % je !----------------------------------------------------------------------- !---------------- 3d prognostic fields --------------------------- if ( id_udt > 0 ) & used = send_data ( id_udt, Var_dt%u(is:ie,vs:ve,:), Time, & mask=Masks%Vel%mask(is:ie,vs:ve,:) > 0.5 ) if ( id_vdt > 0 ) & used = send_data ( id_vdt, Var_dt%v(is:ie,vs:ve,:), Time, & mask=Masks%Vel%mask(is:ie,vs:ve,:) > 0.5 ) if ( id_tdt > 0 ) & used = send_data ( id_tdt, Var_dt%t(is:ie,hs:he,:), Time, & mask=Masks%Tmp%mask(is:ie,hs:he,:) > 0.5 ) do n = 1, Var_dt%ntrace if ( id_tracer_tend(n) > 0 ) & used = send_data ( id_tracer_tend(n), Var_dt%r(is:ie,hs:he,:,n), & Time, mask=Masks%Tmp%mask(is:ie,hs:he,:) > 0.5 ) enddo !----------------------------------------------------------------------- end subroutine bgrid_diagnostics_tend !####################################################################### subroutine compute_vorticity ( Hgrid, dp, udp, vdp, vor ) !----------------------------------------------------------------------- ! Computes relative vorticity on B-grid ! Hgrid = horizontal grid constants ! dp = pressure thickness of model layers at temperature points ! udp = zonal wind, u * dp, at velocity points ! vdp = meridional wind, v * dp, at velocity points ! vor = relative vorticity (1/s) at temperature points !----------------------------------------------------------------------- type(horiz_grid_type), intent(in) :: Hgrid real, intent(in), dimension(Hgrid%ilb:,Hgrid%jlb:,:) :: dp, udp, vdp real, intent(out), dimension(Hgrid%ilb:,Hgrid%jlb:,:) :: vor real,dimension(Hgrid%ilb:Hgrid%iub,Hgrid%jlb:Hgrid%jub) :: & vdy, udx, few, fns integer :: i, j, k, is, ie, js, je is = Hgrid % Tmp % is; ie = Hgrid % Tmp % ie js = Hgrid % Tmp % js; je = Hgrid % Tmp % je do k = 1, size(dp,3) do j = js-1, je vdy(:,j) = vdp(:,j,k)*Hgrid%Vel%dy udx(:,j) = udp(:,j,k)*Hgrid%Vel%dx(j) enddo do j = js, je do i = is-1, ie fns(i,j) = (vdy(i,j-1)+vdy(i,j))*0.5 enddo enddo do j = js-1, je do i = is, ie few(i,j) = (udx(i-1,j)+udx(i,j))*0.5 enddo enddo ! ------ vorticity ------ do j = js, je do i = is, ie vor(i,j,k)=((fns(i,j)-fns(i-1,j))-(few(i,j)-few(i,j-1))) & /(dp(i,j,k)*Hgrid%Tmp%area(j)) enddo enddo enddo end subroutine compute_vorticity !####################################################################### subroutine compute_divergence ( Hgrid, dp, udp, vdp, div ) !----------------------------------------------------------------------- ! Computes divergence on B-grid ! Hgrid = horizontal grid constants ! dp = pressure thickness of model layers at temperature points ! udp = zonal wind, u * dp, at velocity points ! vdp = meridional wind, v * dp, at velocity points ! div = divergence (1/s) at temperature points !----------------------------------------------------------------------- type(horiz_grid_type), intent(in) :: Hgrid real, intent(in), dimension(Hgrid%ilb:,Hgrid%jlb:,:) :: dp, udp, vdp real, intent(out), dimension(Hgrid%ilb:,Hgrid%jlb:,:) :: div real,dimension(Hgrid%ilb:Hgrid%iub,Hgrid%jlb:Hgrid%jub) :: & udy, vdx, few, fns integer :: i, j, k, is, ie, js, je is = Hgrid % Tmp % is; ie = Hgrid % Tmp % ie js = Hgrid % Tmp % js; je = Hgrid % Tmp % je do k = 1, size(dp,3) do j = js-1, je udy(:,j) = udp(:,j,k)*Hgrid%Vel%dy vdx(:,j) = vdp(:,j,k)*Hgrid%Vel%dx(j) enddo do j = js, je do i = is-1, ie few(i,j) = (udy(i,j-1)+udy(i,j))*0.5 enddo enddo do j = js-1, je do i = is, ie fns(i,j) = (vdx(i-1,j)+vdx(i,j))*0.5 enddo enddo ! ------ divergence ------ do j = js, je do i = is, ie div(i,j,k)=((few(i,j)+fns(i,j))-(few(i-1,j)+fns(i,j-1))) & /(dp(i,j,k)*Hgrid%Tmp%area(j)) enddo enddo enddo end subroutine compute_divergence !####################################################################### subroutine sea_level_pressure ( psurf, zsurf, pfull, tfull, slp ) real, intent(in), dimension(:,:) :: psurf, zsurf real, intent(in), dimension(:,:,:) :: pfull, tfull real, intent(out), dimension(:,:) :: slp ! psurf = surface pressure ! zsurf = surface geopotential height in meters^2/sec^2 ! pfull = pressure at full model levels ! tfull = temperature at full model levels ! slp = sea level pressure in pascals integer :: i, j, k, kr real :: sig, tbot real, parameter :: TLAPSE = 6.5e-3 real, parameter :: GORG = GRAV/(RDGAS*TLAPSE) real, parameter :: MRGOG = -1./GORG do j = 1, size(psurf,2) do i = 1, size(psurf,1) if ( abs(zsurf(i,j)) > 0.0001 ) then !---- get ref level for temp ---- do k = 1, size(tfull,3) sig = pfull(i,j,k)/psurf(i,j) if ( sig > 0.8 ) then kr = k exit endif enddo tbot = tfull(i,j,kr) * sig ** MRGOG slp(i,j) = psurf(i,j) * ( 1.0 + TLAPSE * zsurf(i,j) / (tbot*GRAV) ) ** GORG else slp(i,j) = psurf(i,j) endif enddo enddo end subroutine sea_level_pressure !####################################################################### end module bgrid_diagnostics_mod