From ec9804ad137b01f19bd077e0ea77239cc13a3055 Mon Sep 17 00:00:00 2001 From: Matthew Thompson Date: Thu, 30 Nov 2023 14:43:05 -0500 Subject: [PATCH 1/2] Separate ncar_gwd into its own repo --- .../GEOSgwd_GridComp/.gitignore | 3 + .../GEOSgwd_GridComp/CMakeLists.txt | 22 +- .../GEOSgwd_GridComp/ncar_gwd/coords_1d.F90 | 148 -- .../GEOSgwd_GridComp/ncar_gwd/gw_common.F90 | 969 ------------- .../GEOSgwd_GridComp/ncar_gwd/gw_convect.F90 | 711 --------- .../GEOSgwd_GridComp/ncar_gwd/gw_drag.F90 | 283 ---- .../GEOSgwd_GridComp/ncar_gwd/gw_oro.F90 | 404 ------ .../GEOSgwd_GridComp/ncar_gwd/gw_rdg.F90 | 1291 ----------------- .../GEOSgwd_GridComp/ncar_gwd/gw_utils.F90 | 141 -- .../ncar_gwd/interpolate_data.F90 | 1094 -------------- .../ncar_gwd/linear_1d_operators.F90 | 1187 --------------- .../ncar_gwd/vdiff_lu_solver.F90 | 203 --- 12 files changed, 12 insertions(+), 6444 deletions(-) create mode 100644 GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/.gitignore delete mode 100644 GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/coords_1d.F90 delete mode 100644 GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/gw_common.F90 delete mode 100644 GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/gw_convect.F90 delete mode 100644 GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/gw_drag.F90 delete mode 100644 GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/gw_oro.F90 delete mode 100644 GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/gw_rdg.F90 delete mode 100644 GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/gw_utils.F90 delete mode 100644 GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/interpolate_data.F90 delete mode 100644 GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/linear_1d_operators.F90 delete mode 100644 GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/vdiff_lu_solver.F90 diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/.gitignore b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/.gitignore new file mode 100644 index 000000000..08c178560 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/.gitignore @@ -0,0 +1,3 @@ +@ncar_gwd/ +ncar_gwd/ +ncar_gwd@/ diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/CMakeLists.txt b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/CMakeLists.txt index c9c5140a2..b986152dc 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/CMakeLists.txt +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/CMakeLists.txt @@ -1,34 +1,30 @@ esma_set_this () +set subdirs( + ncar_gwd + ) + set (srcs gw_drag.F90 GEOS_GwdGridComp.F90 machine.F90 gwdps.f gwdc.f - ncar_gwd/linear_1d_operators.F90 - ncar_gwd/interpolate_data.F90 - ncar_gwd/vdiff_lu_solver.F90 - ncar_gwd/coords_1d.F90 - ncar_gwd/gw_utils.F90 - ncar_gwd/gw_common.F90 - ncar_gwd/gw_convect.F90 - ncar_gwd/gw_rdg.F90 - ncar_gwd/gw_oro.F90 - ncar_gwd/gw_drag.F90 ) set (resource_files GWD_GridComp.rc ) - - install( FILES ${resource_files} DESTINATION etc ) -esma_add_library (${this} SRCS ${srcs} DEPENDENCIES GEOS_Shared MAPL esmf NetCDF::NetCDF_Fortran) +esma_add_library (${this} + SRCS ${srcs} + SUBCOMPONENTS ${subdirs} + DEPENDENCIES GEOS_Shared MAPL esmf NetCDF::NetCDF_Fortran + ) # CMake has an OpenMP issue with NAG Fortran: https://gitlab.kitware.com/cmake/cmake/-/issues/21280 if (NOT CMAKE_Fortran_COMPILER_ID MATCHES "NAG") diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/coords_1d.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/coords_1d.F90 deleted file mode 100644 index 1930bd8df..000000000 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/coords_1d.F90 +++ /dev/null @@ -1,148 +0,0 @@ -module coords_1d - -! This module defines the Coords1D type, which is intended to to cache -! commonly used information derived from a collection of sets of 1-D -! coordinates. - -implicit none -private - -public :: Coords1D - -type :: Coords1D - ! Number of sets of coordinates in the object. - integer :: n = 0 - ! Number of coordinates in each set. - integer :: d = 0 - - ! All fields below will be allocated with first dimension "n". - ! The second dimension is d+1 for ifc, d for mid, del, and rdel, and - ! d-1 for dst and rdst. - - ! Cell interface coordinates. - real, allocatable :: ifc(:,:) - ! Coordinates at cell mid-points. - real, allocatable :: mid(:,:) - ! Width of cells. - real, allocatable :: del(:,:) - ! Distance between cell midpoints. - real, allocatable :: dst(:,:) - ! Reciprocals: 1/del and 1/dst. - real, allocatable :: rdel(:,:) - real, allocatable :: rdst(:,:) - contains - procedure :: section - procedure :: finalize -end type Coords1D - -interface Coords1D - module procedure new_Coords1D_from_fields - module procedure new_Coords1D_from_int -end interface - -contains - -! Constructor to create an object from existing data. -function new_Coords1D_from_fields(ifc, mid, del, dst, & - rdel, rdst) result(coords) - real, intent(in) :: ifc(:,:) - real, intent(in) :: mid(:,:) - real, intent(in) :: del(:,:) - real, intent(in) :: dst(:,:) - real, intent(in) :: rdel(:,:) - real, intent(in) :: rdst(:,:) - type(Coords1D) :: coords - - coords = allocate_coords(size(ifc, 1), size(ifc, 2) - 1) - - coords%ifc = ifc - coords%mid = mid - coords%del = del - coords%dst = dst - coords%rdel = rdel - coords%rdst = rdst - -end function new_Coords1D_from_fields - -! Constructor if you only have interface coordinates; derives all the other -! fields. -function new_Coords1D_from_int(ifc) result(coords) - real, intent(in) :: ifc(:,:) - type(Coords1D) :: coords - - coords = allocate_coords(size(ifc, 1), size(ifc, 2) - 1) - - coords%ifc = ifc - coords%mid = 0.5 * (ifc(:,:coords%d)+ifc(:,2:)) - coords%del = coords%ifc(:,2:) - coords%ifc(:,:coords%d) - coords%dst = coords%mid(:,2:) - coords%mid(:,:coords%d-1) - coords%rdel = 1./coords%del - coords%rdst = 1./coords%dst - -end function new_Coords1D_from_int - -! Create a new Coords1D object that is a subsection of some other object, -! e.g. if you want only the first m coordinates, use d_bnds=[1, m]. -! -! Originally this used pointers, but it was found to actually be cheaper -! in practice just to make a copy, especially since pointers can impede -! optimization. -function section(self, n_bnds, d_bnds) - class(Coords1D), intent(in) :: self - integer, intent(in) :: n_bnds(2), d_bnds(2) - type(Coords1D) :: section - - section = allocate_coords(n_bnds(2)-n_bnds(1)+1, d_bnds(2)-d_bnds(1)+1) - - section%ifc = self%ifc(n_bnds(1):n_bnds(2),d_bnds(1):d_bnds(2)+1) - section%mid = self%mid(n_bnds(1):n_bnds(2),d_bnds(1):d_bnds(2)) - section%del = self%del(n_bnds(1):n_bnds(2),d_bnds(1):d_bnds(2)) - section%dst = self%dst(n_bnds(1):n_bnds(2),d_bnds(1):d_bnds(2)-1) - section%rdel = self%rdel(n_bnds(1):n_bnds(2),d_bnds(1):d_bnds(2)) - section%rdst = self%rdst(n_bnds(1):n_bnds(2),d_bnds(1):d_bnds(2)-1) - -end function section - -! Quick utility to get allocate each array with the correct size. -function allocate_coords(n, d) result(coords) - integer, intent(in) :: n, d - type(Coords1D) :: coords - - coords%n = n - coords%d = d - - allocate(coords%ifc(coords%n,coords%d+1)) - allocate(coords%mid(coords%n,coords%d)) - allocate(coords%del(coords%n,coords%d)) - allocate(coords%dst(coords%n,coords%d-1)) - allocate(coords%rdel(coords%n,coords%d)) - allocate(coords%rdst(coords%n,coords%d-1)) - -end function allocate_coords - -! Deallocate and reset to initial state. -subroutine finalize(self) - class(Coords1D), intent(inout) :: self - - self%n = 0 - self%d = 0 - - call guarded_deallocate(self%ifc) - call guarded_deallocate(self%mid) - call guarded_deallocate(self%del) - call guarded_deallocate(self%dst) - call guarded_deallocate(self%rdel) - call guarded_deallocate(self%rdst) - -contains - - subroutine guarded_deallocate(array) - real, allocatable :: array(:,:) - - if (allocated(array)) deallocate(array) - - end subroutine guarded_deallocate - -end subroutine finalize - -end module coords_1d diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/gw_common.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/gw_common.F90 deleted file mode 100644 index f25d54783..000000000 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/gw_common.F90 +++ /dev/null @@ -1,969 +0,0 @@ -module gw_common - -use gw_utils, only: GW_PRC, midpoint_interp -! -! This module contains code common to different gravity wave -! parameterizations. -! -implicit none -private - -! Public interface. - -public :: GWBand - -public :: gw_common_init -public :: gw_newtonian_set -public :: gw_prof -public :: gw_drag_prof -public :: qbo_hdepth_scaling -public :: calc_taucd, momentum_flux, momentum_fixer -public :: energy_momentum_adjust, energy_change, energy_fixer -public :: hr_cf - -public :: west, east, north, south -public :: pi -public :: gravit -public :: rair -public :: cpair - -! Index the cardinal directions. -integer, parameter :: west = 1 -integer, parameter :: east = 2 -integer, parameter :: south = 3 -integer, parameter :: north = 4 - -!++jtb (03/2020) -! Some physical constants used by GW codes -!----------------------------------------- -! 3.14159... -real(GW_PRC), parameter :: pi = acos(-1._GW_PRC) -! Acceleration due to gravity. -real(GW_PRC), protected :: gravit = huge(1._GW_PRC) -! Gas constant for dry air. -real(GW_PRC), protected :: rair = huge(1._GW_PRC) -! Specific heat for dry air. -real(GW_PRC), protected :: cpair = huge(1._GW_PRC) -! rair/gravit -real(GW_PRC) :: rog = huge(1._GW_PRC) - - -! Scaling factor for generating QBO -real(GW_PRC), protected :: qbo_hdepth_scaling -! Whether or not to enforce an upper boundary condition of tau = 0. -logical :: tau_0_ubc = .false. -! Inverse Prandtl number. -real(GW_PRC) :: prndl -! Heating rate conversion factor -real(GW_PRC), protected :: hr_cf - - -! -! Private variables -! - -! Interface levels for gravity wave sources. -integer :: ktop = huge(1) - -! Background diffusivity. -real(GW_PRC), parameter :: dback = 0.05_GW_PRC - - -! Newtonian cooling coefficients. -!real, allocatable :: alpha(:) ! AOO global save/alloctable variable not thread-safe - -! -! Limits to keep values reasonable. -! - -! Minimum non-zero stress. -real(GW_PRC), parameter :: taumin = 1.e-10_GW_PRC -! Maximum allowed change in u-c (before efficiency applied). -real(GW_PRC), parameter :: umcfac = 0.5_GW_PRC -! Minimum value of (u-c)**2. -real(GW_PRC), parameter :: ubmc2mn = 0.01_GW_PRC - -! Type describing a band of wavelengths into which gravity waves can be -! emitted. -! Currently this has to have uniform spacing (i.e. adjacent elements of -! cref are exactly dc apart). -type :: GWBand - ! Dimension of the spectrum. - integer :: ngwv - ! Delta between nearest phase speeds [m/s]. - real :: dc - ! Reference speeds [m/s]. - real, allocatable :: cref(:) - ! Critical Froude number, squared (usually 1, but CAM3 used 0.5). - real :: fcrit2 - ! Horizontal wave number [1/m]. - real :: kwv - ! Effective horizontal wave number [1/m] (fcrit2*kwv). - real :: effkwv -end type GWBand - -interface GWBand - module procedure new_GWBand -end interface - -contains - -!========================================================================== - -! Constructor for a GWBand that calculates derived components. -function new_GWBand(ngwv, dc, fcrit2, wavelength) result(band) - ! Used directly to set the type's components. - integer, intent(in) :: ngwv - real, intent(in) :: dc - real, intent(in) :: fcrit2 - ! Wavelength in meters. - real, intent(in) :: wavelength - - ! Output. - type(GWBand) :: band - - ! Wavenumber index. - integer :: l - - ! Simple assignments. - band%ngwv = ngwv - band%dc = dc - band%fcrit2 = fcrit2 - - ! Uniform phase speed reference grid. - allocate(band%cref(-ngwv:ngwv)) - band%cref = [( dc * l, l = -ngwv, ngwv )] - - ! Wavenumber and effective wavenumber come from the wavelength. - band%kwv = 2._GW_PRC*pi / wavelength - band%effkwv = band%fcrit2 * band%kwv - -end function new_GWBand - -!========================================================================== - -subroutine gw_common_init( & - tau_0_ubc_in, ktop_in, gravit_in, rair_in, cpair_in, & - prndl_in, qbo_hdepth_scaling_in, hr_cf_in, errstring) - - logical, intent(in) :: tau_0_ubc_in - integer, intent(in) :: ktop_in - real, intent(in) :: gravit_in - real, intent(in) :: rair_in ! Gas constant for dry air (J kg-1 K-1) - real, intent(in) :: cpair_in ! Heat cap. for dry air (J kg-1 K-1) - real, intent(in) :: prndl_in - real, intent(in) :: qbo_hdepth_scaling_in - real, intent(in) :: hr_cf_in - ! Report any errors from this routine. - character(len=*), intent(out) :: errstring - - integer :: pver - integer :: ierr - - errstring = "" - - tau_0_ubc = tau_0_ubc_in - ktop = ktop_in - gravit = gravit_in - rair = rair_in - cpair = cpair_in - prndl = prndl_in - qbo_hdepth_scaling = qbo_hdepth_scaling_in - hr_cf = hr_cf_in - - rog = rair/gravit - -end subroutine gw_common_init - -!========================================================================== - -!!subroutine gw_prof (ncol, pver, p, t, rhoi, nm, ni) -subroutine gw_prof (ncol, pver, pint, pmid , t, rhoi, nm, ni) - !----------------------------------------------------------------------- - ! Compute profiles of background state quantities for the multiple - ! gravity wave drag parameterization. - ! - ! The parameterization is assumed to operate only where water vapor - ! concentrations are negligible in determining the density. - !----------------------------------------------------------------------- - !------------------------------Arguments-------------------------------- - ! Column and vertical dimensions. - integer, intent(in) :: ncol,pver - ! Pressure coordinates. - real, intent(in) :: pmid(ncol,pver) - real, intent(in) :: pint(ncol,pver+1) - - ! Midpoint temperatures. - real, intent(in) :: t(ncol,pver) - - ! Interface density. - real, intent(out) :: rhoi(ncol,pver+1) - ! Midpoint and interface Brunt-Vaisalla frequencies. - real, intent(out) :: nm(ncol,pver), ni(ncol,pver+1) - - !---------------------------Local Storage------------------------------- - ! Column and level indices. - integer :: i,k - - ! dt/dp - real(GW_PRC) :: dtdp - ! Brunt-Vaisalla frequency squared. - real(GW_PRC) :: n2 - - ! Interface temperature. - real(GW_PRC) :: ti(ncol,pver+1) - - ! Minimum value of Brunt-Vaisalla frequency squared. - real(GW_PRC), parameter :: n2min = 5.e-5_GW_PRC - - !------------------------------------------------------------------------ - ! Determine the interface densities and Brunt-Vaisala frequencies. - !------------------------------------------------------------------------ - - - - ! The top interface values are calculated assuming an isothermal - ! atmosphere above the top level. - k = 1 - do i = 1, ncol - ti(i,k) = t(i,k) - rhoi(i,k) = pint(i,k) / (rair*ti(i,k)) - ni(i,k) = sqrt(gravit*gravit / (cpair*ti(i,k))) - end do - - ! Interior points use centered differences. - ti(:,2:pver) = midpoint_interp(real(t,GW_PRC)) - do k = 2, pver - do i = 1, ncol - rhoi(i,k) = pint(i,k) / (rair*ti(i,k)) - dtdp = (t(i,k)-t(i,k-1)) / ( pmid(i,k)-pmid(i,k-1) ) ! * p%rdst(i,k-1) - n2 = gravit*gravit/ti(i,k) * (1._GW_PRC/cpair - rhoi(i,k)*dtdp) - ni(i,k) = sqrt(max(n2min, n2)) - end do - end do - - ! Bottom interface uses bottom level temperature, density; next interface - ! B-V frequency. - k = pver+1 - do i = 1, ncol - ti(i,k) = t(i,k-1) - rhoi(i,k) = pint(i,k) / (rair*ti(i,k)) - ni(i,k) = ni(i,k-1) - end do - - !------------------------------------------------------------------------ - ! Determine the midpoint Brunt-Vaisala frequencies. - !------------------------------------------------------------------------ - nm = midpoint_interp(ni) - -end subroutine gw_prof - -!========================================================================== -subroutine gw_drag_prof(ncol, pver, band, pint, delp, rdelp, & - src_level, tend_level, dt, t, & - piln, rhoi, nm, ni, ubm, ubi, xv, yv, & - c, kvtt, tau, utgw, vtgw, & - ttgw, gwut, alpha, ro_adjust, tau_adjust, & - kwvrdg, satfac_in, tndmax_in ) - - !----------------------------------------------------------------------- - ! Solve for the drag profile from the multiple gravity wave drag - ! parameterization. - ! 1. scan up from the wave source to determine the stress profile - ! 2. scan down the stress profile to determine the tendencies - ! => apply bounds to the tendency - ! a. from wkb solution - ! b. from computational stability constraints - ! => adjust stress on interface below to reflect actual bounded - ! tendency - !----------------------------------------------------------------------- - - use linear_1d_operators, only: TriDiagDecomp - - !------------------------------Arguments-------------------------------- - ! Column and vertical dimension. - integer, intent(in) :: ncol,pver - ! Wavelengths. - type(GWBand), intent(in) :: band - ! Pressure coordinates. - ! Interface pressures. - real, intent(in) :: pint(ncol,pver+1) - ! Delta Interface pressures. - real, intent(in) :: delp(ncol,pver) - ! Inverse of Delta Interface pressures. - real, intent(in) :: rdelp(ncol,pver) - ! Level from which gravity waves are propagated upward. - integer, intent(in) :: src_level(ncol) - ! Lowest level where wind tendencies are calculated. - integer, intent(in) :: tend_level(ncol) - ! Using tend_level > src_level allows the orographic waves to prescribe - ! wave propagation up to a certain level, but then allow wind tendencies - ! and adjustments to tau below that level. - - ! Time step. - real, intent(in) :: dt - - ! Midpoint and interface temperatures. - real, intent(in) :: t(ncol,pver) - ! Log of interface pressures. - real, intent(in) :: piln(ncol,pver+1) - ! Interface densities. - real, intent(in) :: rhoi(ncol,pver+1) - ! Midpoint and interface Brunt-Vaisalla frequencies. - real, intent(in) :: nm(ncol,pver), ni(ncol,pver+1) - ! Projection of wind at midpoints and interfaces. - real, intent(in) :: ubm(ncol,pver), ubi(ncol,pver+1) - ! Unit vectors of source wind (zonal and meridional components). - real, intent(in) :: xv(ncol), yv(ncol) - ! Wave phase speeds for each column. - real(GW_PRC), intent(in) :: c(ncol,-band%ngwv:band%ngwv) - ! Molecular thermal diffusivity. - real, intent(in) :: kvtt(ncol,pver+1) - -!++jtb -! remove q and dse for now (3/26/20) - ! Constituent array. - !real(GW_PRC), intent(in) :: q(:,:,:) - ! Dry static energy. - !real(GW_PRC), intent(in) :: dse(ncol,pver) -!--jtb - - ! Wave Reynolds stress. - real(GW_PRC), intent(inout) :: tau(ncol,-band%ngwv:band%ngwv,pver+1) - ! Zonal/meridional wind tendencies. - real, intent(out) :: utgw(ncol,pver), vtgw(ncol,pver) - ! Gravity wave heating tendency. - real, intent(out) :: ttgw(ncol,pver) -!++jtb 3/2020 - ! Gravity wave constituent tendency. - !real(GW_PRC), intent(out) :: qtgw(:,:,:) -!--jtb - - ! Gravity wave wind tendency for each wave. - real(GW_PRC), intent(out) :: gwut(ncol,pver,-band%ngwv:band%ngwv) - real, intent(in) :: alpha(:) - - ! Adjustment parameter for IGWs. - real, intent(in), optional :: & - ro_adjust(ncol,-band%ngwv:band%ngwv,pver+1) - - ! Adjustment parameter for TAU. - real, intent(in), optional :: & - tau_adjust(ncol,pver+1) - - ! Diagnosed horizontal wavenumber for ridges. - real, intent(in), optional :: & - kwvrdg(ncol) - - ! Factor for saturation calculation. Here backwards - ! compatibility. I believe it should be 1.0 (jtb). - ! Looks like it has been 2.0 for a while in CAM. - real, intent(in), optional :: & - satfac_in - - real, intent(in), optional :: tndmax_in - - !---------------------------Local storage------------------------------- - - ! Level, wavenumber, constituent and column loop indices. - integer :: k, l, m, i, kmax - - ! Lowest tendency and source levels. - integer :: kbot_tend, kbot_src - - ! "Total" and saturation diffusivity. - real(GW_PRC) :: d(ncol) - ! Imaginary part of vertical wavenumber. - real(GW_PRC) :: mi(ncol) - ! Stress after damping. - real(GW_PRC) :: taudmp(ncol) - ! Saturation stress. - real(GW_PRC) :: tausat(ncol) - ! (ub-c) and (ub-c)**2 - real(GW_PRC) :: ubmc(ncol), ubmc2(ncol) - ! Temporary ubar tendencies (overall, and at wave l). - real(GW_PRC) :: ubt(ncol,pver), ubtl(ncol) - real(GW_PRC) :: wrk(ncol) - ! Temporary effkwv - real(GW_PRC) :: effkwv(ncol) - - ! saturation factor. Defaults to 2.0 - ! unless overidden by satfac_in - real(GW_PRC) :: satfac - - real(GW_PRC) :: tndmax - - real(GW_PRC) :: near_zero = tiny(1.0_GW_PRC) - - ! LU decomposition. - type(TriDiagDecomp) :: decomp - - if (present(satfac_in)) then - satfac = satfac_in - else - satfac = 2.0 - endif - -! Maximum wind tendency from stress divergence (before efficiency applied). - if (present(tndmax_in)) then - tndmax = tndmax_in - else - tndmax = 400._GW_PRC / 86400._GW_PRC - endif - - ! Lowest levels that loops need to iterate over. - kbot_tend = maxval(tend_level) - kbot_src = maxval(src_level) - - ! Initialize gravity wave drag tendencies to zero. - - utgw = 0.0 - vtgw = 0.0 - gwut = 0.0 - ttgw = 0.0 - - ! Workaround floating point exception issues on Intel by initializing - ! everything that's first set in a where block. - mi = 0.0 - taudmp = 0.0 - tausat = 0.0 - ubmc = 0.0 - ubmc2 = 0.0 - wrk = 0.0 - - if (present(kwvrdg)) then - effkwv = kwvrdg - else - effkwv = band%effkwv - endif - - !------------------------------------------------------------------------ - ! Compute the stress profiles and diffusivities - !------------------------------------------------------------------------ - - ! Loop from bottom to top to get stress profiles. -! !$OMP parallel do default(none) shared(kbot_src,ktop,kvtt,band,ubi,c,effkwv,rhoi,ni,satfac, & -! !$OMP ro_adjust,ncol,alpha,piln,t,rog,src_level,tau_adjust,tau) & -! !$OMP private(k,d,l,i,tausat,taudmp,ubmc,ubmc2,wrk,mi) - do k = kbot_src, ktop, -1 !++ but this is in model now - - ! Determine the diffusivity for each column. - - d = dback + kvtt(:,k) - - do l = -band%ngwv, band%ngwv - - do i=1,ncol - - tausat(i) = 0.0 - taudmp(i) = 0.0 - - if (src_level(i) >= k) then - - ! Determine the absolute value of the saturation stress. - ! Define critical levels where the sign of (u-c) changes between - ! interfaces. - ubmc(i) = ubi(i,k) - c(i,l) - - ! Test to see if u-c has the same sign here as the level below. - if (ubmc(i) > 0.0 .eqv. ubi(i,k+1) > c(i,l)) then - if (ni(i,k) /= 0.0) & - tausat(i) = abs( effkwv(i) * rhoi(i,k) * ubmc(i)**3 / & - (satfac*ni(i,k)) ) - if (present(ro_adjust)) & - tausat(i) = tausat(i) * sqrt(ro_adjust(i,l,k)) - if (present(tau_adjust)) & - tausat(i) = tausat(i) * tau_adjust(i,k) - endif - - ! Compute stress for each wave. The stress at this level is the - ! min of the saturation stress and the stress at the level below - ! reduced by damping. The sign of the stress must be the same as - ! at the level below. - ubmc2(i) = max(ubmc(i)**2, ubmc2mn) - mi(i) = ni(i,k) / (2.0 * effkwv(i) * ubmc2(i)) * & ! Is this 2.0 related to satfac? - (alpha(k) + ni(i,k)**2/ubmc2(i) * d(i)) - wrk(i) = -2.0*mi(i)*rog*t(i,k)*(piln(i,k+1) - piln(i,k)) - wrk(i) = max( wrk(i), -200.0 ) * exp(wrk(i)) - taudmp(i) = tau(i,l,k+1) * exp(wrk(i)) - ! For some reason, PGI 14.1 loses bit-for-bit reproducibility if - ! we limit tau, so instead limit the arrays used to set it. - if (tausat(i) <= taumin) tausat(i) = 0.0 - if (taudmp(i) <= taumin) taudmp(i) = 0.0 - tau(i,l,k) = min(taudmp(i), tausat(i)) - - endif - end do - end do - end do - - ! Force tau at the top of the model to zero, if requested. - if (tau_0_ubc) tau(:,:,ktop) = 0.0 - - !------------------------------------------------------------------------ - ! Compute the tendencies from the stress divergence. - !------------------------------------------------------------------------ - - - ! Loop over levels from top to bottom -! !$OMP parallel do default(none) shared(kbot_tend,ktop,band,ncol,tau,delp,rdelp,c,ubm,dt,gravit,utgw,vtgw, & -! !$OMP gwut,ubt,xv,yv,tend_level,near_zero) & -! !$OMP private(k,l,i,ubtl) - do k = ktop, kbot_tend - - ! Accumulate the mean wind tendency over wavenumber. - ubt(:,k) = 0.0 - - do l = -band%ngwv, band%ngwv ! loop over wave - - do i=1,ncol - - ! Determine the wind tendency, including excess stress carried down - ! from above. - ubtl(i) = gravit * (tau(i,l,k+1)-tau(i,l,k)) * rdelp(i,k) ! p%rdel(i,k) !/1/D_pint - - - ! Apply first tendency limit to maintain numerical stability. - ! Enforce du/dt < |c-u|/dt so u-c cannot change sign - ! (u^n+1 = u^n + du/dt * dt) - ! The limiter is somewhat stricter, so that we don't come anywhere - ! near reversing c-u. - ubtl(i) = min(ubtl(i), umcfac * abs(c(i,l)-ubm(i,k)) / dt) - - if (k <= tend_level(i)) then - - ! Save tendency for each wave (for later computation of kzz). - ! sign function returns magnitude of ubtl with sign of c-ubm - ! Renders ubt/ubm check for mountain waves unecessary - gwut(i,k,l) = sign(ubtl(i), c(i,l)-ubm(i,k)) - ubt(i,k) = ubt(i,k) + gwut(i,k,l) - - end if - - end do - - end do - - do l = -band%ngwv, band%ngwv - ! Redetermine the effective stress on the interface below from the - ! wind tendency. If the wind tendency was limited above, then the - ! new stress will be smaller than the old stress, causing stress - ! divergence in the next layer down. This smoothes large stress - ! divergences downward while conserving total stress. - ! Include a protection on SMALL gwut to prevent floating point - ! issues. - !-------------------------------------------------- - do i=1,ncol - if ( abs(gwut(i,k,l)) < near_zero ) then - gwut(i,k,l) = 0.0 - end if - if (k <= tend_level(i)) then - tau(i,l,k+1) = tau(i,l,k) + & - abs(gwut(i,k,l)) * delp(i,k) / gravit - !!! abs(gwut(i,k,l)) * p%del(i,k) / gravit - end if - end do - end do - - ! Project the mean wind tendency onto the components. - do i=1,ncol - if (k <= tend_level(i)) then - utgw(i,k) = ubt(i,k) * xv(i) - vtgw(i,k) = ubt(i,k) * yv(i) - end if - end do - - ! End of level loop. - end do - - ! Evaluate second temperature tendency term: Conversion of kinetic - ! energy into thermal. -! !$OMP parallel do default(none) shared(kbot_tend,ktop,band,ttgw,ubm,c,gwut) & -! !$OMP private(k,l) - do k = ktop, kbot_tend - do l = -band%ngwv, band%ngwv - ttgw(:,k) = ttgw(:,k) - (ubm(:,k) - c(:,l)) * gwut(:,k,l) - end do - end do - ttgw = ttgw / cpair - -end subroutine gw_drag_prof - - -!========================================================================== - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -subroutine gw_newtonian_set( pver, pref, alpha ) - -use interpolate_data, only: lininterp - - integer, intent(in) :: pver - real, intent(in) :: pref( pver+1 ) -! Newtonian cooling coefficients. - real, intent(inout) :: alpha(:) ! Make alpha argument instead of global for correct behavior under threading - - ! Levels of pre-calculated Newtonian cooling (1/day). - ! The following profile is digitized from: - ! Wehrbein and Leovy (JAS, 39, 1532-1544, 1982) figure 5 - - integer :: k - integer, parameter :: nalph = 71 - - real :: alpha0(nalph) = [ & - 0.1, 0.1, 0.1, 0.1, & - 0.1, 0.1, 0.1, 0.1, & - 0.1, 0.1, 0.10133333, 0.104, & - 0.108, 0.112, 0.116, 0.12066667, & - 0.126, 0.132, 0.138, 0.144, & - 0.15133333, 0.16, 0.17, 0.18, & - 0.19, 0.19933333, 0.208, 0.216, & - 0.224, 0.232, 0.23466667, 0.232, & - 0.224, 0.216, 0.208, 0.20133333, & - 0.196, 0.192, 0.188, 0.184, & - 0.18266667, 0.184, 0.188, 0.192, & - 0.196, 0.19333333, 0.184, 0.168, & - 0.152, 0.136, 0.12133333, 0.108, & - 0.096, 0.084, 0.072, 0.061, & - 0.051, 0.042, 0.033, 0.024, & - 0.017666667, 0.014, 0.013, 0.012, & - 0.011, 0.010333333, 0.01, 0.01, & - 0.01, 0.01, 0.01 & - ] - - ! Pressure levels that were used to calculate alpha0 (hPa). - real :: palph(nalph) = [ & - 2.06115E-06, 2.74280E-06, 3.64988E-06, 4.85694E-06, & - 6.46319E-06, 8.60065E-06, 1.14450E-05, 1.52300E-05, & - 2.02667E-05, 2.69692E-05, 3.58882E-05, 4.77568E-05, & - 6.35507E-05, 8.45676E-05, 0.000112535, 0.000149752, & - 0.000199277, 0.000265180, 0.000352878, 0.000469579, & - 0.000624875, 0.000831529, 0.00110653, 0.00147247, & - 0.00195943, 0.00260744, 0.00346975, 0.00461724, & - 0.00614421, 0.00817618, 0.0108801, 0.0144783, & - 0.0192665, 0.0256382, 0.0341170, 0.0453999, & - 0.0604142, 0.0803939, 0.106981, 0.142361, & - 0.189442, 0.252093, 0.335463, 0.446404, & - 0.594036, 0.790490, 1.05192, 1.39980, & - 1.86273, 2.47875, 3.29851, 4.38936, & - 5.84098, 7.77266, 10.3432, 13.7638, & - 18.3156, 24.3728, 32.4332, 43.1593, & - 57.4326, 76.4263, 101.701, 135.335, & - 180.092, 239.651, 318.907, 424.373, & - 564.718, 751.477, 1000. & - ] - - ! pre-calculated newtonian damping: - ! * convert to 1/s - ! * ensure it is not smaller than 1e-6 - ! * convert palph from hpa to pa - do k=1,nalph - alpha0(k) = alpha0(k) / 86400.0 - alpha0(k) = max(alpha0(k), 1.e-6) - palph(k) = palph(k)*1.e2 - end do - - !allocate (alpha(pver+1)) ! Make alpha local instead of global for correct behavior under threading - - ! interpolate to current vertical grid and obtain alpha - call lininterp (alpha0, palph, nalph , alpha, pref, pver+1) - -end subroutine gw_newtonian_set - -!========================================================================== - -! Calculate Reynolds stress for waves propagating in each cardinal -! direction. - -function calc_taucd(ncol, pver, ngwv, tend_level, tau, c, xv, yv, ubi) & - result(taucd) - - ! Column and gravity wave wavenumber dimensions. - integer, intent(in) :: ncol, pver, ngwv - ! Lowest level where wind tendencies are calculated. - integer, intent(in) :: tend_level(:) - ! Wave Reynolds stress. - real(GW_PRC), intent(in) :: tau(:,-ngwv:,:) - ! Wave phase speeds for each column. - real(GW_PRC), intent(in) :: c(:,-ngwv:) - ! Unit vectors of source wind (zonal and meridional components). - real, intent(in) :: xv(:), yv(:) - ! Projection of wind at interfaces. - real, intent(in) :: ubi(:,:) - - real :: taucd(ncol,pver+1,4) - - ! Indices. - integer :: i, k, l - - ! ubi at tend_level. - real :: ubi_tend(ncol) - - ! Signed wave Reynolds stress. - real :: tausg(ncol) - - ! Reynolds stress for waves propagating behind and forward of the wind. - real :: taub(ncol) - real :: tauf(ncol) - - taucd = 0. - tausg = 0. - - ubi_tend = (/ (ubi(i,tend_level(i)+1), i = 1, ncol) /) - - do k = ktop, maxval(tend_level)+1 - - taub = 0. - tauf = 0. - - do l = -ngwv, ngwv - where (k-1 <= tend_level) - - tausg = sign(tau(:,l,k), c(:,l)-ubi(:,k)) - - where ( c(:,l) < ubi_tend ) - taub = taub + tausg - elsewhere - tauf = tauf + tausg - end where - - end where - end do - - where (k-1 <= tend_level) - where (xv > 0.) - taucd(:,k,east) = tauf * xv - taucd(:,k,west) = taub * xv - elsewhere - taucd(:,k,east) = taub * xv - taucd(:,k,west) = tauf * xv - end where - - where ( yv > 0.) - taucd(:,k,north) = tauf * yv - taucd(:,k,south) = taub * yv - elsewhere - taucd(:,k,north) = taub * yv - taucd(:,k,south) = tauf * yv - end where - end where - - end do - -end function calc_taucd - -!========================================================================== - -! Calculate the amount of momentum conveyed from below the gravity wave -! region, to the region where gravity waves are calculated. -subroutine momentum_flux(tend_level, taucd, um_flux, vm_flux) - - ! Bottom stress level. - integer, intent(in) :: tend_level(:) - ! Projected stresses. - real, intent(in) :: taucd(:,:,:) - ! Components of momentum change sourced from the bottom. - real, intent(out) :: um_flux(:), vm_flux(:) - - integer :: i - - ! Tendency for U & V below source level. - do i = 1, size(tend_level) - um_flux(i) = taucd(i,tend_level(i)+1, east) + & - taucd(i,tend_level(i)+1, west) - vm_flux(i) = taucd(i,tend_level(i)+1,north) + & - taucd(i,tend_level(i)+1,south) - end do - -end subroutine momentum_flux - -!========================================================================== - -! Subtracts a change in momentum in the gravity wave levels from wind -! tendencies in lower levels, ensuring momentum conservation. -subroutine momentum_fixer(ncol, pver, tend_level, p, um_flux, vm_flux, utgw, vtgw) - - integer, intent(in) :: ncol, pver - ! Bottom stress level. - integer, intent(in) :: tend_level(ncol) - ! Pressure coordinates. - real, intent(in) :: p(ncol,pver+1) - ! Components of momentum change sourced from the bottom. - real, intent(in) :: um_flux(ncol), vm_flux(ncol) - ! Wind tendencies. - real, intent(inout) :: utgw(ncol,pver), vtgw(ncol,pver) - - ! Indices. - integer :: i, k - ! Reciprocal of total mass. - real :: rdm(ncol) - - ! Total mass from ground to source level: rho*dz = dp/gravit - do i = 1, ncol - rdm(i) = gravit/(p(i,pver+1)-p(i,tend_level(i)+1)) - end do - - do k = minval(tend_level)+1, pver - where (k > tend_level) - utgw(:,k) = utgw(:,k) + -um_flux*rdm - vtgw(:,k) = vtgw(:,k) + -vm_flux*rdm - end where - end do - -end subroutine momentum_fixer - -!========================================================================== - -! Calculate the change in total energy from tendencies up to this point. -subroutine energy_change(ncol,pver, dt, pdel, u, v, dudt, dvdt, dtdt, de) - - integer, intent(in) :: ncol, pver - ! Time step. - real, intent(in) :: dt - ! Pressure coordinates. - real, intent(in) :: pdel(ncol,pver+1) - ! Winds at start of time step. - real, intent(in) :: u(ncol,pver), v(ncol,pver) - ! Wind tendencies. - real, intent(in) :: dudt(ncol,pver), dvdt(ncol,pver) - ! Heating tendency. - real, intent(in) :: dtdt(ncol,pver) - ! Change in energy. - real, intent(out) :: de(ncol) - - ! Level index. - integer :: k - - ! Net gain/loss of total energy in the column. - de = 0.0 - do k = 1, pver - de = de + pdel(:,k)/gravit * (cpair*dtdt(:,k) + & - dudt(:,k)*(u(:,k)+dudt(:,k)*0.5*dt) + & - dvdt(:,k)*(v(:,k)+dvdt(:,k)*0.5*dt) ) - end do - -end subroutine energy_change - -!========================================================================== - -! Subtract change in energy from the heating tendency in the levels below -! the gravity wave region. -subroutine energy_fixer(ncol, pver, tend_level, pint, de, ttgw) - - integer, intent(in) :: ncol, pver - ! Bottom stress level. - integer, intent(in) :: tend_level(ncol) - ! Pressure coordinates. - real, intent(in) :: pint(ncol,pver+1) - ! Change in energy. - real, intent(in) :: de(ncol) - ! Heating tendency. - real, intent(inout) :: ttgw(ncol,pver) - - ! Column/level indices. - integer :: i, k - ! Energy change to apply divided by all the mass it is spread across. - real :: de_dm(ncol) - - do i = 1, ncol - de_dm(i) = -de(i)*gravit/(pint(i,pver+1)-pint(i,tend_level(i)+1)) - end do - - ! Subtract net gain/loss of total energy below tend_level. - do k = minval(tend_level)+1, pver - where (k > tend_level) - ttgw(:,k) = ttgw(:,k) + de_dm/cpair - end where - end do - -end subroutine energy_fixer - -!========================================================================== - -subroutine energy_momentum_adjust(ncol, pver, band, pint, delp, u, v, dt, c, tau, & - effgw, t, ubm, ubi, xv, yv, utgw, vtgw, ttgw, & - tend_level, tndmax_in) - - integer, intent(in) :: ncol, pver - ! Wavelengths. - type(GWBand), intent(in) :: band - ! Pressure interfaces. - real, intent(in) :: pint(ncol,pver+1) - ! Pressure thickness. - real, intent(in) :: delp(ncol,pver) - ! Winds - real, intent(in) :: u(ncol,pver) ! Midpoint zonal winds. ( m s-1) - real, intent(in) :: v(ncol,pver) ! Midpoint meridional winds. ( m s-1) - ! timestep - real, intent(in) :: dt - ! Wave phase speeds for each column. - real(GW_PRC), intent(in) :: c(ncol,-band%ngwv:band%ngwv) - ! Wave Reynolds stress. - real(GW_PRC), intent(in) :: tau(ncol,-band%ngwv:band%ngwv,pver+1) - ! Efficiency - real, intent(in) :: effgw(ncol) - ! Temperature and winds - real, intent(in) :: t(ncol,pver), ubm(ncol,pver), ubi(ncol,pver) - ! projected winds. - real, intent(in) :: xv(ncol), yv(ncol) - ! tendency level index index - integer, intent(in) :: tend_level(ncol) - ! Tendency limiter - real, intent(in), optional :: tndmax_in - ! Tendencies. - real, intent(inout) :: utgw(ncol,pver) - real, intent(inout) :: vtgw(ncol,pver) - real, intent(inout) :: ttgw(ncol,pver) - - real :: taucd(ncol,pver+1,4) - real :: um_flux(ncol), vm_flux(ncol) - real :: de(ncol) - - real :: zlb,pm,rhom,cmu,fpmx,fpmy,fe,fpe,fpml,fpel,fpmt,fpet,dusrcl,dvsrcl,dtsrcl - real :: tndmax,utfac,uhtmax - integer :: ktop - integer :: kbot(ncol) - - ! Level index. - integer :: i,k,l - -! Maximum wind tendency from stress divergence (before efficiency applied). - if (present(tndmax_in)) then - tndmax = tndmax_in - else - tndmax = 400._GW_PRC / 86400._GW_PRC - endif - - ktop=1 - do i=1,ncol - !--------------------------------------------------------------------------------------- - ! Apply efficiency factor and tendency limiter to prevent unrealistically strong forcing - !--------------------------------------------------------------------------------------- - uhtmax = 0.0 - utfac = 1.0 - do k = ktop, pver - uhtmax = max(sqrt(utgw(i,k)**2 + vtgw(i,k)**2), uhtmax) - end do - if (uhtmax > tndmax) utfac = tndmax/uhtmax - do k = ktop, pver - utgw(i,k) = utgw(i,k)*effgw(i)*utfac - vtgw(i,k) = vtgw(i,k)*effgw(i)*utfac - ttgw(i,k) = ttgw(i,k)*effgw(i)*utfac - end do - end do ! i=1,ncol - -#ifdef ENERGY_ADJUST - if (band%ngwv /= 0) then - ! compute momentum and energy flux changes - taucd = calc_taucd(ncol, pver, band%ngwv, tend_level, tau, c, xv, yv, ubi) - call momentum_flux(tend_level, taucd, um_flux, vm_flux) - call energy_change(ncol,pver, dt, delp, u, v, utgw, vtgw, ttgw, de) - ! add sub-source fixers to tendencies - call momentum_fixer(ncol, pver, tend_level, pint, um_flux, vm_flux, utgw, vtgw) - call energy_fixer(ncol, pver, tend_level, pint, de, ttgw) - endif -#endif - -end subroutine energy_momentum_adjust - -!========================================================================== -end module gw_common diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/gw_convect.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/gw_convect.F90 deleted file mode 100644 index f245a8f08..000000000 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/gw_convect.F90 +++ /dev/null @@ -1,711 +0,0 @@ -module gw_convect - -! -! This module handles gravity waves from convection, and was extracted from -! gw_drag in May 2013. -! - - use gw_utils, only: GW_PRC, GW_R8, get_unit_vector, dot_2d, midpoint_interp - use gw_common, only: GWBand, qbo_hdepth_scaling, gw_drag_prof, hr_cf, & - calc_taucd, momentum_flux, momentum_fixer, & - energy_momentum_adjust, energy_change, energy_fixer - - use MAPL_ConstantsMod, only: MAPL_RGAS, MAPL_CP, MAPL_GRAV - -implicit none -private - -public :: BeresSourceDesc -public :: gw_beres_ifc -public :: gw_beres_src -public :: gw_beres_init - -real, parameter :: PI = 3.14159265358979323846 ! pi -real, parameter :: rad2deg = 180./PI - -type :: BeresSourceDesc - logical :: active - ! Whether wind speeds are shifted to be relative to storm cells. - logical :: storm_shift - ! Heating depths below this value [m] will be ignored. - real :: min_hdepth - ! Source for wave spectrum - real :: spectrum_source - ! Index for level where wind speed is used as the source speed. - real, allocatable :: k(:) - ! tendency limiter - real :: tndmax - ! Table bounds, for convenience. (Could be inferred from shape(mfcc).) - integer :: maxh - integer :: maxuh - ! Heating depths [m]. - real, allocatable :: hd(:) - ! Table of source spectra. - real, allocatable :: mfcc(:,:,:) - ! Forced background for extratropics - real, allocatable :: taubck(:,:) -end type BeresSourceDesc - - -contains - -!========================================================================== - -!------------------------------------ -subroutine gw_beres_init (file_name, band, desc, pgwv, gw_dc, fcrit2, wavelength, & - spectrum_source, min_hdepth, storm_shift, taubgnd, tndmax, active, ncol, lats) -#include - - character(len=*), intent(in) :: file_name - type(GWBand), intent(inout) :: band - - type(BeresSourceDesc), intent(inout) :: desc - - integer, intent(in) :: pgwv, ncol - real, intent(in) :: gw_dc, fcrit2, wavelength - real, intent(in) :: spectrum_source, min_hdepth, taubgnd, tndmax - logical, intent(in) :: storm_shift, active - real, intent(in) :: lats(ncol) - - ! Stuff for Beres convective gravity wave source. - real(GW_R8), allocatable :: mfcc(:,:,:), hdcc(:) - integer :: hd_mfcc , mw_mfcc, ps_mfcc, ngwv_file, ps_mfcc_mid - - ! For forced background extratropical wave speed - real :: c4, latdeg, flat_gw - real, allocatable :: c0(:), cw4(:) - integer :: i, kc - - ! Vars needed by NetCDF operators - integer :: ncid, dimid, varid, status - - status = nf_open(file_name , 0, ncid) - - status = NF_INQ_DIMID(ncid, 'PS', dimid) - IF (status .NE. NF_NOERR) CALL HANDLE_ERR(status) - status = NF_INQ_DIMLEN(ncid, dimid, ps_mfcc ) - - status = NF_INQ_DIMID(ncid, 'MW', dimid) - IF (status .NE. NF_NOERR) CALL HANDLE_ERR(status) - status = NF_INQ_DIMLEN(ncid, dimid, mw_mfcc ) - - status = NF_INQ_DIMID(ncid, 'HD', dimid) - IF (status .NE. NF_NOERR) CALL HANDLE_ERR(status) - status = NF_INQ_DIMLEN(ncid, dimid, hd_mfcc ) - - allocate( mfcc(hd_mfcc , mw_mfcc, ps_mfcc) ) - allocate( hdcc(hd_mfcc) ) - - status = NF_INQ_VARID(ncid, 'HD', varid) - IF (status .NE. NF_NOERR) CALL HANDLE_ERR(status) - status = NF_GET_VAR_DOUBLE(ncid, varid, hdcc ) - IF (status .NE. NF_NOERR) CALL HANDLE_ERR(status) - - status = NF_INQ_VARID(ncid, 'mfcc', varid) - IF (status .NE. NF_NOERR) CALL HANDLE_ERR(status) - status = NF_GET_VAR_DOUBLE(ncid, varid, mfcc ) - IF (status .NE. NF_NOERR) CALL HANDLE_ERR(status) - - status = nf_close (ncid) - - band = GWBand(pgwv, gw_dc, fcrit2, wavelength ) - - ! These dimensions; {HD,MW,PS}_MFCC, came from Beres forcing file. - - ! Get HD (heating depth) dimension. - desc%maxh = HD_MFCC - - ! Get MW (mean wind) dimension. - desc%maxuh = MW_MFCC - - ! Get PS (phase speed) dimension. - ngwv_file = PS_MFCC - - ! Number in each direction is half of total (and minus phase speed of 0). - desc%maxuh = (desc%maxuh-1)/2 - ! midpoint of spectrum in netcdf file is ps_mfcc (odd number) -1 divided by 2, plus 1 - ! E.g., ps_mfcc = 5. So, ps_mfcc_mid = 3 - ! 1 2 3 4 5 - ! -2 -1 0 +1 +2 - ps_mfcc_mid= (ngwv_file-1)/2 - - desc%active = active - if (active) then - - allocate(desc%hd(desc%maxh) , stat=status ) - - allocate(desc%mfcc(desc%maxh,-desc%maxuh:desc%maxuh,-band%ngwv:band%ngwv), stat=status ) - - desc%mfcc( : , -desc%maxuh:desc%maxuh , -band%ngwv :band%ngwv ) & - = mfcc( :, : , -band%ngwv+ps_mfcc_mid:band%ngwv+ps_mfcc_mid ) - - ! While not currently documented in the file, it uses kilometers. Convert - ! to meters. - desc%hd = hdcc * 1000.0 - - ! Source level index allocated, filled later - desc%spectrum_source = spectrum_source - allocate(desc%k(ncol)) - - desc%min_hdepth = min_hdepth - - desc%storm_shift = storm_shift - - desc%tndmax = tndmax - - ! Intialize forced background wave speeds - allocate(desc%taubck(ncol,-band%ngwv:band%ngwv)) - allocate(c0(-band%ngwv:band%ngwv)) - allocate(cw4(-band%ngwv:band%ngwv)) - desc%taubck = 0.0 - c0 = 0.0 - cw4 = 0.0 - do kc = -4,4 - c4 = 10.0*kc - cw4(kc) = exp(-(c4/30.)**2) - enddo - do kc = -band%ngwv,band%ngwv - c0(kc) = 10.0*(4.0/real(band%ngwv))*kc - desc%taubck(:,kc) = exp(-(c0(kc)/30.)**2) - enddo - do i=1,ncol - ! include forced background stress in extra tropics - ! Determine the background stress at c=0 - ! Include dependence on latitude: - latdeg = lats(i)*rad2deg - if (-15.3 < latdeg .and. latdeg < 15.3) then - flat_gw = 0.10 - else if (latdeg > -31. .and. latdeg <= -15.3) then - flat_gw = 0.10 - else if (latdeg < 31. .and. latdeg >= 15.3) then - flat_gw = 0.10 - else if (latdeg > -60. .and. latdeg <= -31.) then - flat_gw = 0.50*exp(-((abs(latdeg)-60.)/23.)**2) - else if (latdeg < 60. .and. latdeg >= 31.) then - flat_gw = 0.50*exp(-((abs(latdeg)-60.)/23.)**2) - else if (latdeg <= -60.) then - flat_gw = 0.50*exp(-((abs(latdeg)-60.)/70.)**2) - else if (latdeg >= 60.) then - flat_gw = 0.50*exp(-((abs(latdeg)-60.)/70.)**2) - end if - desc%taubck(i,:) = taubgnd*0.001*flat_gw*desc%taubck(i,:)*(sum(cw4)/sum(desc%taubck(i,:))) - enddo - deallocate( c0, cw4 ) - end if - -end subroutine gw_beres_init - -!------------------------------------ -subroutine gw_beres_src(ncol, pver, band, desc, u, v, & - netdt, zm, src_level, tend_level, tau, ubm, ubi, xv, yv, & - c, hdepth, maxq0, lats, dqcdt) -!----------------------------------------------------------------------- -! Driver for multiple gravity wave drag parameterization. -! -! The parameterization is assumed to operate only where water vapor -! concentrations are negligible in determining the density. -! -! Beres, J.H., M.J. Alexander, and J.R. Holton, 2004: "A method of -! specifying the gravity wave spectrum above convection based on latent -! heating properties and background wind". J. Atmos. Sci., Vol 61, No. 3, -! pp. 324-337. -! -!----------------------------------------------------------------------- - !!!use gw_utils, only: get_unit_vector, dot_2d, midpoint_interp - !!!use gw_common, only: GWBand, qbo_hdepth_scaling - -!------------------------------Arguments-------------------------------- - ! Column and vertical dimensions. - integer, intent(in) :: ncol, pver - - ! Wavelengths triggered by convection. - type(GWBand), intent(in) :: band - - ! Settings for convection type (e.g. deep vs shallow). - type(BeresSourceDesc), intent(in) :: desc - - ! Midpoint zonal/meridional winds. - real, intent(in) :: u(ncol,pver), v(ncol,pver) - ! Heating rate due to convection. - real, intent(in) :: netdt(:,:) - ! Midpoint altitudes. - real, intent(in) :: zm(ncol,pver) - ! latitudes. - real, intent(in) :: lats(ncol) - - ! Indices of top gravity wave source level and lowest level where wind - ! tendencies are allowed. - integer, intent(out) :: src_level(ncol) - integer, intent(out) :: tend_level(ncol) - - ! Wave Reynolds stress. - real(GW_PRC), intent(out) :: tau(ncol,-band%ngwv:band%ngwv,pver+1) - ! Projection of wind at midpoints and interfaces. - real, intent(out) :: ubm(ncol,pver) - real, intent(out) :: ubi(ncol,pver+1) - ! Unit vectors of source wind (zonal and meridional components). - real, intent(out) :: xv(ncol), yv(ncol) - ! Phase speeds. - real(GW_PRC), intent(out) :: c(ncol,-band%ngwv:band%ngwv) - - ! Heating depth [m] and maximum heating in each column. - real, intent(out) :: hdepth(ncol), maxq0(ncol) - - ! Condensate tendency due to large-scale (kg kg-1 s-1) - real, optional, intent(in) :: dqcdt(ncol,pver) ! Condensate tendency due to large-scale (kg kg-1 s-1) - -!---------------------------Local Storage------------------------------- - ! Column and level indices. - integer :: i, k - - ! Zonal/meridional wind at roughly the level where the convection occurs. - real :: uconv(ncol), vconv(ncol), ubi1d(ncol) - - ! Maximum heating rate. - real(GW_PRC) :: q0(ncol) - - ! Bottom/top heating range index. - integer :: boti(ncol), topi(ncol) - ! Index for looking up heating depth dimension in the table. - integer :: hd_idx(ncol) - ! Mean wind in heating region. - real(GW_PRC) :: uh(ncol) - ! Min/max wavenumber for critical level filtering. - integer :: Umini(ncol), Umaxi(ncol) - ! Source level tau for a column. - real(GW_PRC) :: tau0(-band%ngwv:band%ngwv) - ! Speed of convective cells relative to storm. - real(GW_PRC) :: CS(ncol) - ! Index to shift spectra relative to ground. - integer :: shift - - ! Averaging length. - real, parameter :: AL = 1.0e5 - integer :: thread - - !---------------------------------------------------------------------- - ! Initialize tau array - !---------------------------------------------------------------------- - tau = 0.0 - hdepth = 0.0 - q0 = 0.0 - tau0 = 0.0 - ubi = 0.0 - - !------------------------------------------------------------------------ - ! Determine wind and unit vectors approximately at the source level, then - ! project winds. - !------------------------------------------------------------------------ - - ! Source wind speed and direction. - do i=1,ncol - uconv(i) = u(i,desc%k(i)) - vconv(i) = v(i,desc%k(i)) - enddo - - ! Get the unit vector components and magnitude at the source level. - ubi1d = 0.0 - call get_unit_vector(uconv, vconv, xv, yv, ubi1d) - do i=1,ncol - ubi(i,desc%k(i)+1) = ubi1d(i) - enddo - - ! Project the local wind at midpoints onto the source wind. - do k = 1, pver - ubm(:,k) = dot_2d(u(:,k), v(:,k), xv, yv) - end do - - ! Compute the interface wind projection by averaging the midpoint winds. - ! Use the top level wind at the top interface. - ubi(:,1) = ubm(:,1) - - ubi(:,2:pver) = midpoint_interp(ubm) - - !----------------------------------------------------------------------- - ! Calculate heating depth. - ! - ! Heating depth is defined as the first height range from the bottom in - ! which heating rate is continuously positive. - !----------------------------------------------------------------------- - - ! First find the indices for the top and bottom of the heating range. - boti = 0 - topi = 0 - do k = pver, 1, -1 - do i = 1, ncol - if (boti(i) == 0) then - ! Detect if we are outside the maximum range (where z = 20 km). - if (zm(i,k) >= 20000.0) then - boti(i) = k - topi(i) = k - else - ! First spot where heating rate is positive. - if (netdt(i,k) > 0.0) boti(i) = k - end if - else if (topi(i) == 0) then - ! Detect if we are outside the maximum range (z = 20 km). - if (zm(i,k) >= 20000.0) then - topi(i) = k - else - ! First spot where heating rate is no longer positive. - if (netdt(i,k) <= 0.0) topi(i) = k - end if - end if - end do - ! When all done, exit. - if (all(topi /= 0)) exit - end do - - ! Heating depth in m. - hdepth = [ ( (zm(i,topi(i))-zm(i,boti(i))), i = 1, ncol ) ] - - ! J. Richter: this is an effective reduction of the GW phase speeds (needed to drive the QBO) - hdepth = hdepth*qbo_hdepth_scaling - - hd_idx = index_of_nearest(hdepth, desc%hd) - - ! hd_idx=0 signals that a heating depth is too shallow, i.e. that it is - ! either not big enough for the lowest table entry, or it is below the - ! minimum allowed for this convection type. - ! Values above the max in the table still get the highest value, though. - where (hdepth < max(desc%min_hdepth, desc%hd(1))) hd_idx = 0 - - ! Maximum heating rate. - do k = minval(topi), maxval(boti) - where (k >= topi .and. k <= boti) - q0 = max(q0, netdt(:,k)) - end where - end do - - !output max heating rate in K/day - maxq0 = q0*86400.0 - - ! Multipy by conversion factor - q0 = q0 * hr_cf - - if (desc%storm_shift) then - - ! Find the cell speed where the storm speed is > 10 m/s. - ! Storm speed is taken to be the source wind speed. - do i=1,ncol - CS(i) = sign(max(abs(ubm(i,desc%k(i)))-10.0, 0.0), ubm(i,desc%k(i))) - enddo - - ! Average wind in heating region, relative to storm cells. - uh = 0.0 - do k = minval(topi), maxval(boti) - where (k >= topi .and. k <= boti) - uh = uh + ubm(:,k)/(boti-topi+1) - end where - end do - - uh = uh - CS - - else - - ! For shallow convection, wind is relative to ground, and "heating - ! region" wind is just the source level wind. - do i=1,ncol - uh(i) = ubm(i,desc%k(i)) - enddo - - end if - - ! Limit uh to table range. - uh = min(uh, real(desc%maxuh)) - uh = max(uh, -real(desc%maxuh)) - - ! Speeds for critical level filtering. - Umini = band%ngwv - Umaxi = -band%ngwv - do k = minval(topi), maxval(boti) - where (k >= topi .and. k <= boti) - Umini = min(Umini, nint(ubm(:,k)/band%dc)) - Umaxi = max(Umaxi, nint(ubm(:,k)/band%dc)) - end where - end do - - Umini = max(Umini, -band%ngwv) - Umaxi = min(Umaxi, band%ngwv) - - !----------------------------------------------------------------------- - ! Gravity wave sources - !----------------------------------------------------------------------- - ! Start loop over all columns. - !----------------------------------------------------------------------- - do i=1,ncol - - !--------------------------------------------------------------------- - ! Look up spectrum only if the heating depth is large enough, else set - ! tau0 = 0. - !--------------------------------------------------------------------- - - if (hd_idx(i) > 0) then - - !------------------------------------------------------------------ - ! Look up the spectrum using depth and uh. - !------------------------------------------------------------------ - - tau0 = desc%mfcc(hd_idx(i),nint(uh(i)),:) - - if (desc%storm_shift) then - ! For deep convection, the wind was relative to storm cells, so - ! shift the spectrum so that it is now relative to the ground. - shift = -nint(CS(i)/band%dc) - tau0 = eoshift(tau0, shift) - end if - - ! Adjust magnitude. - tau0 = tau0*(q0(i)**2)/AL - - ! Adjust for critical level filtering. - tau0(Umini(i):Umaxi(i)) = 0.0 - - tau(i,:,topi(i)+1) = tau0 - - else - if (present(dqcdt)) then - if (dqcdt(i,desc%k(i)) > 1.e-8) then ! frontal region (large-scale forcing) - ! include forced background stress in extra tropical large-scale systems - ! Set the phase speeds and wave numbers in the direction of the source wind. - ! Set the source stress magnitude (positive only, note that the sign of the - ! stress is the same as (c-u). - tau(i,:,desc%k(i)+1) = desc%taubck(i,:) - topi(i) = desc%k(i) - endif - endif - endif - enddo - !----------------------------------------------------------------------- - ! End loop over all columns. - !----------------------------------------------------------------------- - - ! Output the source level. - src_level = topi - tend_level = topi - - ! Set phase speeds; just use reference speeds. - c = spread(band%cref, 1, ncol) - -end subroutine gw_beres_src - - - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -!!! Main Interface -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -subroutine gw_beres_ifc( band, & - ncol, pver, dt, effgw_dp, & - u, v, t, pref, pint, delp, rdelp, piln, & - zm, zi, nm, ni, rhoi, kvtt, & - netdt,desc,lats, alpha, & - utgw,vtgw,ttgw,flx_heat,dqcdt) - - type(BeresSourceDesc), intent(inout) :: desc - type(GWBand), intent(in) :: band ! I hate this variable ... it just hides information from view - integer, intent(in) :: ncol ! number of atmospheric columns - integer, intent(in) :: pver ! number of vertical layers - real, intent(in) :: dt ! Time step. - real, intent(in) :: effgw_dp - - real, intent(in) :: u(ncol,pver) ! Midpoint zonal winds. ( m s-1) - real, intent(in) :: v(ncol,pver) ! Midpoint meridional winds. ( m s-1) - real, intent(in) :: t(ncol,pver) ! Midpoint temperatures. (K) - real, intent(in) :: netdt(ncol,pver) ! Convective heating rate (K s-1) - real, intent(in) :: pref(pver+1) ! Reference pressure at interfaces (Pa !!! ) - real, intent(in) :: piln(ncol,pver+1) ! Log of interface pressures. - real, intent(in) :: pint(ncol,pver+1) ! Interface pressures. (Pa) - real, intent(in) :: delp(ncol,pver) ! Layer pressures thickness. (Pa) - real, intent(in) :: rdelp(ncol,pver) ! Inverse pressure thickness. (Pa-1) - real, intent(in) :: zm(ncol,pver) ! Midpoint altitudes above ground (m). - real, intent(in) :: zi(ncol,pver+1) ! Interface altitudes above ground (m). - real, intent(in) :: nm(ncol,pver) ! Midpoint Brunt-Vaisalla frequencies (s-1). - real, intent(in) :: ni(ncol,pver+1) ! Interface Brunt-Vaisalla frequencies (s-1). - real, intent(in) :: rhoi(ncol,pver+1) ! Interface density (kg m-3). - real, intent(in) :: kvtt(ncol,pver+1) ! Molecular thermal diffusivity. - - real, intent(in) :: lats(ncol) ! latitudes - real, intent(in) :: alpha(:) - - real, intent(out) :: utgw(ncol,pver) ! zonal wind tendency - real, intent(out) :: vtgw(ncol,pver) ! meridional wind tendency - real, intent(out) :: ttgw(ncol,pver) ! temperature tendency - real, intent(inout) :: flx_heat(ncol) ! Energy change - - real, optional, intent(in) :: dqcdt(ncol,pver) ! Condensate tendency due to large-scale (kg kg-1 s-1) - - !---------------------------Local storage------------------------------- - - integer :: k, m, nn - - real(GW_PRC), allocatable :: tau(:,:,:) ! wave Reynolds stress - ! gravity wave wind tendency for each wave - real(GW_PRC), allocatable :: gwut(:,:,:) - ! Wave phase speeds for each column - real(GW_PRC), allocatable :: c(:,:) - - ! Efficiency for a gravity wave source. - real :: effgw(ncol) - - ! Momentum fluxes used by fixer. - real :: um_flux(ncol), vm_flux(ncol) - - ! Energy change used by fixer. - real :: de(ncol) - - ! Reynolds stress for waves propagating in each cardinal direction. - real :: taucd(ncol,pver+1,4) - - ! Indices of top gravity wave source level and lowest level where wind - ! tendencies are allowed. - integer :: src_level(ncol) - integer :: tend_level(ncol) - - ! Projection of wind at midpoints and interfaces. - real :: ubm(ncol,pver) - real :: ubi(ncol,pver+1) - - ! Unit vectors of source wind (zonal and meridional components). - real :: xv(ncol) - real :: yv(ncol) - - ! Heating depth [m] and maximum heating in each column. - real :: hdepth(ncol), maxq0(ncol) - - real :: pint_adj(ncol,pver+1) - real :: zfac_layer - - character(len=1) :: cn - character(len=9) :: fname(4) - - integer :: i,j,l - - !---------------------------------------------------------------------------- - - - ! Allocate wavenumber fields. - allocate(tau(ncol,-band%ngwv:band%ngwv,pver+1)) - allocate(gwut(ncol,pver,-band%ngwv:band%ngwv)) - allocate(c(ncol,-band%ngwv:band%ngwv)) - - ! Efficiency of gravity wave momentum transfer. - ! This is really only to remove the pole points. - where (pi/2.0 - abs(lats(:ncol)) >= 1.e-4 ) !-4*epsilon(1.0)) - effgw = effgw_dp - elsewhere - effgw = 0.0 - end where - -!GEOS pressure scaling near model top - zfac_layer = 1.0e2 ! 1mb - do k=1,pver+1 - do i=1,ncol - pint_adj(i,k) = MIN(1.0,MAX(0.0,(pint(i,k)/zfac_layer)**3)) - enddo - enddo - - do k = 0, pver - ! spectrum source index - if (pref(k+1) < desc%spectrum_source) desc%k(:) = k+1 - end do - - ! Determine wave sources for Beres deep scheme - call gw_beres_src(ncol, pver, band, desc, & - u, v, netdt, zm, src_level, tend_level, tau, & - ubm, ubi, xv, yv, c, hdepth, maxq0, lats, dqcdt=dqcdt) - - ! Solve for the drag profile with orographic sources. - call gw_drag_prof(ncol, pver, band, pint, delp, rdelp, & - src_level, tend_level, dt, t, & - piln, rhoi, nm, ni, ubm, ubi, xv, yv, & - c, kvtt, tau, utgw, vtgw, & - ttgw, gwut, alpha, tau_adjust=pint_adj) - - ! Apply efficiency and limiters - call energy_momentum_adjust(ncol, pver, band, pint, delp, u, v, dt, c, tau, & - effgw, t, ubm, ubi, xv, yv, utgw, vtgw, ttgw, & - tend_level, tndmax_in=desc%tndmax) - - deallocate(tau, gwut, c) - -end subroutine gw_beres_ifc - - -!************************************************************************ -!!handle_err -!************************************************************************ -! -!!ROUTINE: handle_err -!!DESCRIPTION: error handler -!-------------------------------------------------------------------------- - -subroutine handle_err(status) - - implicit none - -#include - - integer status - - if (status .ne. nf_noerr) then - print *, nf_strerror(status) - stop 'Stopped' - endif - -end subroutine handle_err - - - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -subroutine endrun(msg) - - integer :: iulog - - character(len=*), intent(in), optional :: msg ! string to be printed - - iulog=6 - - if (present (msg)) then - write(iulog,*)'ENDRUN:', msg - else - write(iulog,*)'ENDRUN: called without a message string' - end if - - stop - -end subroutine endrun - - - - - - - - - - - - - -! Short routine to get the indices of a set of values rounded to their -! nearest points on a grid. -function index_of_nearest(x, grid) result(idx) - real, intent(in) :: x(:) - real, intent(in) :: grid(:) - - integer :: idx(size(x)) - - real :: interfaces(size(grid)-1) - integer :: i, n - - n = size(grid) - interfaces = (grid(:n-1) + grid(2:))/2.d0 - - idx = 1 - do i = 1, n-1 - where (x > interfaces(i)) idx = i + 1 - end do - -end function index_of_nearest - -end module gw_convect diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/gw_drag.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/gw_drag.F90 deleted file mode 100644 index 5d7644927..000000000 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/gw_drag.F90 +++ /dev/null @@ -1,283 +0,0 @@ - -! $Id$ -module gw_drag_ncar - -!--------------------------------------------------------------------------------- -! Purpose: -! -! Module to compute the forcing due to parameterized gravity waves. Both an -! orographic and an internal source spectrum are considered. -! -! Author: Byron Boville -! In-Sun Song -! -!--------------------------------------------------------------------------------- - - use MAPL_ConstantsMod, only: MAPL_RGAS, MAPL_GRAV - - use gw_rdg, only : gw_rdg_ifc - use gw_oro, only : gw_oro_ifc - use gw_convect, only : BeresSourceDesc, gw_beres_ifc - use gw_common, only : GWBand,gw_prof - - private ! Make default type private to the module -! -! PUBLIC: interfaces -! - public gw_intr_ncar ! interface to actual parameterization - -! -! PRIVATE: Rest of the data and interfaces are private to this module -! - real, parameter :: KWVB = 6.28e-5 ! effective horizontal wave number for background - real, parameter :: KWVBEQ = 6.28e-5/7. ! effective horizontal wave number for background - real, parameter :: KWVO = 6.28e-5 ! effective horizontal wave number for orographic - real, parameter :: FRACLDV = 0.0 ! fraction of stress deposited in low level region - - real, parameter :: MXASYM = 0.1 ! max asymmetry between tau(c) and tau(-c) - real, parameter :: MXRANGE = 0.001 ! max range of tau for all c - real, parameter :: N2MIN = 1.e-8 ! min value of bouyancy frequency - real, parameter :: FCRIT2 = 0.5 ! critical froude number - real, parameter :: OROHMIN = 10. ! min surface displacment height for orographic waves - real, parameter :: OROVMIN = 2.0 ! min wind speed for orographic waves - real, parameter :: TAUBGND = 6.4 ! background source strength (/TAUSCAL) - real, parameter :: TAUMIN = 1.e-10 ! minimum (nonzero) stress - real, parameter :: TAUSCAL = 0.001 ! scale factor for background stress source - real, parameter :: UMCFAC = 0.5 ! factor to limit tendency to prevent reversing u-c - real, parameter :: UBMC2MN = 0.01 ! min (u-c)**2 - real, parameter :: ZLDVCON = 10. ! constant for determining zldv from tau0 - - real, parameter :: ROG = MAPL_RGAS/MAPL_GRAV - real, parameter :: OROKO2 = 0.5 * KWVO ! 1/2 * horizontal wavenumber - real, parameter :: PI_GWD = 4.0*atan(1.0) ! This is *not* MAPL_PI -contains - -!=============================================================================== - - - subroutine gw_intr_ncar(pcols, pver, dt, nrdg, & - beres_dc_desc, beres_sc_desc, beres_band, oro_band, rdg_band, & - pint_dev, t_dev, u_dev, v_dev, & - ht_dc_dev, ht_sc_dev, dqcdt_dev, & - sgh_dev, mxdis_dev, hwdth_dev, clngt_dev, angll_dev, & - anixy_dev, gbxar_dev, kwvrdg_dev, effrdg_dev, pref_dev, & - pmid_dev, pdel_dev, rpdel_dev, lnpint_dev, zm_dev, rlat_dev, & - phis_dev, & - dudt_gwd_dev, dvdt_gwd_dev, dtdt_gwd_dev, & - dudt_org_dev, dvdt_org_dev, dtdt_org_dev, & - taugwdx_dev, taugwdy_dev, & - taubkgx_dev, taubkgy_dev, & - effgworo, effgwbkg, alpha, rc ) - -!----------------------------------------------------------------------- -! Interface for multiple gravity wave drag parameterization. -!----------------------------------------------------------------------- - -!------------------------------Arguments-------------------------------- - integer, intent(in ) :: pcols ! number of columns - integer, intent(in ) :: pver ! number of vertical layers -!++jtb 01/25/21 - integer, intent(in ) :: nrdg ! number of Ridges per gridbox -!--jtb - real, intent(in ) :: dt ! time step - type(GWBand), intent(inout) :: oro_band ! Band descriptor - type(GWBand), intent(inout) :: rdg_band ! Band descriptor - type(GWBand), intent(inout) :: beres_band ! Band descriptor - type(BeresSourceDesc), intent(inout) :: beres_dc_desc ! Table descriptor for DeepCu Beres scheme - type(BeresSourceDesc), intent(inout) :: beres_sc_desc ! Table descriptor for ShallowCu Beres scheme - real, intent(in ) :: effgwbkg ! tendency efficiency for background gwd (Default = 0.125) - real, intent(in ) :: effgworo ! tendency efficiency for orographic gwd (Default = 0.125) - real, intent(in ) :: pint_dev(pcols,pver+1) ! pressure at the layer edges - real, intent(in ) :: t_dev(pcols,pver) ! temperature at layers - real, intent(in ) :: u_dev(pcols,pver) ! zonal wind at layers - real, intent(in ) :: v_dev(pcols,pver) ! meridional wind at layers - real, intent(in ) :: ht_dc_dev(pcols,pver) ! DeepCu heating in layers - real, intent(in ) :: ht_sc_dev(pcols,pver) ! ShallowCu heating in layers - real, intent(in ) :: dqcdt_dev(pcols,pver) ! Condensate tendencies due to large-scale - real, intent(in ) :: sgh_dev(pcols) ! standard deviation of orography -!++jtb 01/25/21 New topo vars - real, intent(in ) :: mxdis_dev(pcols,nrdg) ! obstacle/ridge height - real, intent(in ) :: hwdth_dev(pcols,nrdg) ! obstacle width - real, intent(in ) :: clngt_dev(pcols,nrdg) ! obstacle along-crest length - real, intent(in ) :: angll_dev(pcols,nrdg) ! obstacle orientation - real, intent(in ) :: anixy_dev(pcols,nrdg) ! obstacle ansitropy param - real, intent(in ) :: gbxar_dev(pcols) ! duplicate grid box area - real, intent(in ) :: kwvrdg_dev(pcols,nrdg) ! horizontal wavenumber - real, intent(in ) :: effrdg_dev(pcols,nrdg) ! efficiency of ridge scheme -!!--jtb - real, intent(in ) :: pref_dev(pver+1) ! reference pressure at the layeredges - real, intent(in ) :: pmid_dev(pcols,pver) ! pressure at the layers - real, intent(in ) :: pdel_dev(pcols,pver) ! pressure thickness at the layers - real, intent(in ) :: rpdel_dev(pcols,pver) ! 1.0 / pdel - real, intent(in ) :: lnpint_dev(pcols,pver+1) ! log(pint) - real, intent(in ) :: zm_dev(pcols,pver) ! height above surface at layers - real, intent(in ) :: rlat_dev(pcols) ! latitude in radian - real, intent(in ) :: phis_dev(pcols) ! surface geopotential - - real, intent( out) :: dudt_gwd_dev(pcols,pver) ! zonal wind tendency at layer - real, intent( out) :: dvdt_gwd_dev(pcols,pver) ! meridional wind tendency at layer - real, intent( out) :: dtdt_gwd_dev(pcols,pver) ! temperature tendency at layer - real, intent( out) :: dudt_org_dev(pcols,pver) ! zonal wind tendency at layer due to orography GWD - real, intent( out) :: dvdt_org_dev(pcols,pver) ! meridional wind tendency at layer due to orography GWD - real, intent( out) :: dtdt_org_dev(pcols,pver) ! temperature tendency at layer due to orography GWD - real, intent( out) :: taugwdx_dev(pcols) ! zonal gravity wave surface stress - real, intent( out) :: taugwdy_dev(pcols) ! meridional gravity wave surface stress - real, intent( out) :: taubkgx_dev(pcols) ! zonal gravity wave background stress - real, intent( out) :: taubkgy_dev(pcols) ! meridional gravity wave background stress - real, intent( in) :: alpha(:) - - integer, optional, intent(out) :: RC ! return code - - -#ifndef GPU_MAXLEVS -#define GPU_MAXLEVS pver -#endif - -#ifndef MAXPGWV -#define MAXPGWV pgwv -#endif - -#define CAMGWCODE - -#ifdef CAMGWCODE -#define IFCDIMS 1:pver+1 -#else -#define IFCDIMS 0:pver -#endif - -!---------------------------Local storage------------------------------- - - integer :: i,k,kc ! loop indexes - integer :: kbotoro ! launch-level index for orographic - integer :: kbotbg ! launch-level index for background - integer :: ktopbg, ktoporo ! top interface of gwd region - integer :: kldv ! top interface of low level stress divergence region - integer :: kldvmn ! min value of kldv - integer :: ksrc ! index of top interface of source region - integer :: ksrcmn ! min value of ksrc - - real :: ttgw(pcols,pver) ! temperature tendency - real :: utgw(pcols,pver) ! zonal wind tendency - real :: vtgw(pcols,pver) ! meridional wind tendency - - real :: z(pcols) ! Surface elevation - real :: zi(pcols,pver+1) ! interface heights above ground - real :: ni(pcols,pver+1) ! interface Brunt-Vaisalla frequency - real :: nm(pcols,pver) ! midpoint Brunt-Vaisalla frequency - !!!real :: rdpldv ! 1/dp across low level divergence region - real :: rhoi(pcols,pver+1) ! interface density - real :: ti(pcols,pver+1) ! interface temperature - real :: ubi(pcols,pver+1) ! projection of wind at interfaces - real :: ubm(pcols,pver) ! projection of wind at midpoints - real :: xv(pcols) ! unit vectors of source wind (x) - real :: yv(pcols) ! unit vectors of source wind (y) - real :: kvtt(pcols,pver+1) ! Molecular thermal diffusivity. - - real :: maxq0(pcols),hdepth(pcols) - real :: flx_heat(pcols) - - real :: rdg_cd_llb - integer :: pverp, pcnst - logical :: trpd_leewv - -!----------------------------------------------------------------------------- - -! Misc dimensions needed - pverp=pver+1 - pcnst=1 - -! Initialize accumulated tendencies -! and other things ... - dtdt_gwd_dev(:,:) = 0. - dudt_gwd_dev(:,:) = 0. - dvdt_gwd_dev(:,:) = 0. - kvtt(:,:) = 0. - flx_heat(:) = 0.0 - - call gw_prof (pcols , pver, pint_dev , pmid_dev , t_dev , rhoi, nm, ni ) - - z = phis_dev/MAPL_GRAV - - zi(:,pver+1) = 0.0 - do k=2,pver - zi(:,k) = 0.5 * ( zm_dev(:,k-1) + zm_dev(:,k) ) - end do - zi(:,1) = zi(:,2) + 0.5*( zm_dev(:,1) - zm_dev(:,2) ) - - ! DeepCu BKG - if (beres_dc_desc%active .and. effgwbkg>0.0) then - call gw_beres_ifc( beres_band, & - pcols, pver, dt , effgwbkg, & - u_dev , v_dev, t_dev, & - pref_dev, pint_dev, & - pdel_dev , rpdel_dev, lnpint_dev, & - zm_dev, zi, & - nm, ni, rhoi, kvtt, & - ht_dc_dev,beres_dc_desc,rlat_dev, alpha, & - utgw, vtgw, ttgw, flx_heat) - dudt_gwd_dev = dudt_gwd_dev + utgw - dvdt_gwd_dev = dvdt_gwd_dev + vtgw - dtdt_gwd_dev = dtdt_gwd_dev + ttgw - endif - ! ShallowCu BKG - if (beres_sc_desc%active .and. effgwbkg>0.0) then - call gw_beres_ifc( beres_band, & - pcols, pver, dt , effgwbkg, & - u_dev , v_dev, t_dev, & - pref_dev, pint_dev, & - pdel_dev , rpdel_dev, lnpint_dev, & - zm_dev, zi, & - nm, ni, rhoi, kvtt, & - ht_sc_dev,beres_sc_desc,rlat_dev, alpha, & - utgw, vtgw, ttgw, flx_heat, dqcdt=dqcdt_dev) - dudt_gwd_dev = dudt_gwd_dev + utgw - dvdt_gwd_dev = dvdt_gwd_dev + vtgw - dtdt_gwd_dev = dtdt_gwd_dev + ttgw - endif - ! Orographic - if (effgworo > 0.0) then - if (nrdg > 0) then - trpd_leewv = .FALSE. - rdg_cd_llb = 1.0 - call gw_rdg_ifc( rdg_band, & - pcols, pver, pverp, pcnst, nrdg, dt, & - u_dev , v_dev, t_dev, & - pint_dev, pmid_dev, & - pdel_dev, rpdel_dev, & - lnpint_dev, zm_dev, zi, z, & - ni, nm, rhoi, & - kvtt, & - kwvrdg_dev, effrdg_dev, & - hwdth_dev, clngt_dev, gbxar_dev, & - mxdis_dev, angll_dev, anixy_dev, & - rdg_cd_llb, trpd_leewv, alpha, & - utgw, vtgw, ttgw, flx_heat) - else - call gw_oro_ifc( oro_band, & - pcols, pver, dt , effgworo, & - u_dev , v_dev, t_dev, & - pint_dev, pmid_dev, & - pdel_dev , rpdel_dev, lnpint_dev, & - zm_dev, zi, & - nm, ni, rhoi, kvtt, & - sgh_dev ,rlat_dev, alpha, & - utgw, vtgw, ttgw) - endif - dudt_org_dev = utgw - dvdt_org_dev = vtgw - dtdt_org_dev = ttgw - dudt_gwd_dev = dudt_gwd_dev + dudt_org_dev - dvdt_gwd_dev = dvdt_gwd_dev + dvdt_org_dev - dtdt_gwd_dev = dtdt_gwd_dev + dtdt_org_dev - endif - - taugwdx_dev(1:pcols) = 0.0 !zonal gravity wave surface stress - taugwdy_dev(1:pcols) = 0.0 !meridional gravity wave surface stress - taubkgx_dev(1:pcols) = 0.0 !zonal gravity wave background stress - taubkgy_dev(1:pcols) = 0.0 !meridional gravity wave background stress - - return - end subroutine gw_intr_ncar - -end module gw_drag_ncar - diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/gw_oro.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/gw_oro.F90 deleted file mode 100644 index d123d231a..000000000 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/gw_oro.F90 +++ /dev/null @@ -1,404 +0,0 @@ -module gw_oro - -! -! This module handles gravity waves from convection, and was extracted from -! gw_drag in May 2013. -! - - use gw_utils, only: GW_PRC, get_unit_vector, dot_2d, midpoint_interp - use gw_common, only: GWBand, rair, gw_drag_prof, energy_momentum_adjust - -implicit none -private - -public :: gw_oro_ifc -public :: gw_oro_src -public :: gw_oro_init - -real,parameter :: PI = 3.14159265358979323846 ! pi - -real :: gw_oro_south_fac - -contains - -!========================================================================== - -!------------------------------------ -subroutine gw_oro_init (band, gw_dc, fcrit2, wavelength, pgwv, oro_south_fac) -#include - - type(GWBand), intent(inout) :: band - real, intent(in) :: gw_dc,fcrit2,wavelength,oro_south_fac - integer, intent(in) :: pgwv - - - -! Need to call GWBand for oro waves - - band = GWBand(pgwv, gw_dc, fcrit2, wavelength ) - - gw_oro_south_fac = oro_south_fac - -end subroutine gw_oro_init - -!------------------------------------ -subroutine gw_oro_src(ncol,pver, band, & - pint, pmid, delp, & - u, v, t, sgh, zm, nm, & - src_level, tend_level, tau, ubm, ubi, xv, yv, c) - !----------------------------------------------------------------------- - ! Orographic source for multiple gravity wave drag parameterization. - ! - ! The stress is returned for a single wave with c=0, over orography. - ! For points where the orographic variance is small (including ocean), - ! the returned stress is zero. - !----------------------------------------------------------------------- - !!! use gw_utils, only: get_unit_vector, dot_2d, midpoint_interp - !!! use gw_common, only: GWBand,rair - - ! Column dimension. - integer, intent(in) :: ncol, pver - ! Band to emit orographic waves in. - ! Regardless, we will only ever emit into l = 0. - type(GWBand), intent(in) :: band - ! Pressure coordinates. - !type(Coords1D), intent(in) :: p - - ! Interface pressures. (Pa) - real, intent(in) :: pint(ncol,pver+1) - ! Midpoint pressures. (Pa) - real, intent(in) :: pmid(ncol,pver) - ! Delta Interface pressures. (Pa) - real, intent(in) :: delp(ncol,pver) - - - ! Midpoint zonal/meridional winds. - real, intent(in) :: u(ncol,pver), v(ncol,pver) - ! Midpoint temperatures. - real, intent(in) :: t(ncol,pver) - ! Standard deviation of orography. - real, intent(in) :: sgh(ncol) - ! Midpoint altitudes. - real, intent(in) :: zm(ncol,pver) - ! Midpoint Brunt-Vaisalla frequencies. - real, intent(in) :: nm(ncol,pver) - - ! Indices of top gravity wave source level and lowest level where wind - ! tendencies are allowed. - integer, intent(out) :: src_level(ncol) - integer, intent(out) :: tend_level(ncol) - - ! Wave Reynolds stress. - real(GW_PRC), intent(out) :: tau(ncol,-band%ngwv:band%ngwv,pver+1) - ! Projection of wind at midpoints and interfaces. - real, intent(out) :: ubm(ncol,pver), ubi(ncol,pver+1) - ! Unit vectors of source wind (zonal and meridional components). - real, intent(out) :: xv(ncol), yv(ncol) - ! Phase speeds. - real(GW_PRC), intent(out) :: c(ncol,-band%ngwv:band%ngwv) - - !---------------------------Local Storage------------------------------- - ! Column and level indices. - integer :: i, k - - ! Surface streamline displacement height (2*sgh). - real :: hdsp(ncol) - ! Max orographic standard deviation to use. - real :: sghmax - ! c=0 stress from orography. - real :: tauoro(ncol) - ! Averages over source region. - real :: nsrc(ncol) ! B-V frequency. - real :: rsrc(ncol) ! Density. - real :: usrc(ncol) ! Zonal wind. - real :: vsrc(ncol) ! Meridional wind. - - ! Difference in interface pressure across source region. - real :: dpsrc(ncol) - - ! Limiters (min/max values) - ! min surface displacement height for orographic waves - real, parameter :: orohmin = 10. - ! min wind speed for orographic waves - real, parameter :: orovmin = 2. - -!-------------------------------------------------------------------------- -! Average the basic state variables for the wave source over the depth of -! the orographic standard deviation. Here we assume that the appropiate -! values of wind, stability, etc. for determining the wave source are -! averages over the depth of the atmosphere penterated by the typical -! mountain. -! Reduces to the bottom midpoint values when sgh=0, such as over ocean. -!-------------------------------------------------------------------------- - - hdsp = 2.0 * sgh - - k = pver - src_level = k-1 - rsrc = pmid(:,k)/(rair*t(:,k)) * delp(:,k) - usrc = u(:,k) * delp(:,k) - vsrc = v(:,k) * delp(:,k) - nsrc = nm(:,k)* delp(:,k) - - do k = pver-1, 1, -1 - do i = 1, ncol - if (hdsp(i) > sqrt(zm(i,k)*zm(i,k+1))) then - src_level(i) = k-1 - rsrc(i) = rsrc(i) + & - pmid(i,k) / (rair*t(i,k)) * delp(i,k) - usrc(i) = usrc(i) + u(i,k) * delp(i,k) - vsrc(i) = vsrc(i) + v(i,k) * delp(i,k) - nsrc(i) = nsrc(i) + nm(i,k)* delp(i,k) - end if - end do - ! Break the loop when all source levels found. - if (all(src_level >= k)) exit - end do - - - do i = 1, ncol - dpsrc(i) = pint(i,pver+1) - pint(i,src_level(i)+1) - end do - - - rsrc = rsrc / dpsrc - usrc = usrc / dpsrc - vsrc = vsrc / dpsrc - nsrc = nsrc / dpsrc - - - ! Get the unit vector components and magnitude at the surface. - call get_unit_vector(usrc, vsrc, xv, yv, ubi(:,pver+1)) - - ! Project the local wind at midpoints onto the source wind. - do k = 1, pver - ubm(:,k) = dot_2d(u(:,k), v(:,k), xv, yv) - end do - - ! Compute the interface wind projection by averaging the midpoint winds. - ! Use the top level wind at the top interface. - ubi(:,1) = ubm(:,1) - - ubi(:,2:pver) = midpoint_interp(ubm) - - ! Determine the orographic c=0 source term following McFarlane (1987). - ! Set the source top interface index to pver, if the orographic term is - ! zero. - do i = 1, ncol - if ((ubi(i,pver+1) > orovmin) .and. (hdsp(i) > orohmin)) then - sghmax = band%fcrit2 * (ubi(i,pver+1) / nsrc(i))**2 - tauoro(i) = 0.5 * band%kwv * min(hdsp(i)**2, sghmax) * & - rsrc(i) * nsrc(i) * ubi(i,pver+1) - else - tauoro(i) = 0. - src_level(i) = pver - end if - end do - - ! Set the phase speeds and wave numbers in the direction of the source - ! wind. Set the source stress magnitude (positive only, note that the - ! sign of the stress is the same as (c-u). - tau = 0. - do k = pver, minval(src_level), -1 - where (src_level <= k) tau(:,0,k+1) = tauoro - end do - - ! Allow wind tendencies all the way to the model bottom. - tend_level = pver - - ! No spectrum; phase speed is just 0. - c = 0. - -end subroutine gw_oro_src - - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -!!! Main Interface -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -subroutine gw_oro_ifc( band, & - ncol, pver, dt, effgw_oro, & - u, v, t, pint, pmid, & - delp, rdelp, piln, & - zm, zi, nm, ni, rhoi, kvtt, & - sgh, lats, alpha, & - utgw,vtgw,ttgw) - - type(GWBand), intent(in) :: band ! I hate this variable ... it just hides information from view - integer, intent(in) :: ncol ! number of atmospheric columns - integer, intent(in) :: pver ! number of vertical layers - real, intent(in) :: dt ! Time step. - real, intent(in) :: effgw_oro - - real, intent(in) :: u(ncol,pver) ! Midpoint zonal winds. ( m s-1) - real, intent(in) :: v(ncol,pver) ! Midpoint meridional winds. ( m s-1) - real, intent(in) :: t(ncol,pver) ! Midpoint temperatures. (K) - real, intent(in) :: piln(ncol,pver+1) ! Log of interface pressures. - real, intent(in) :: pmid(ncol,pver) ! Midpoint pressures. (Pa) - real, intent(in) :: pint(ncol,pver+1) ! Interface pressures. (Pa) - real, intent(in) :: delp(ncol,pver) ! Layer pressures thickness. (Pa) - real, intent(in) :: rdelp(ncol,pver) ! Inverse pressure thickness. (Pa-1) - real, intent(in) :: zm(ncol,pver) ! Midpoint altitudes above ground (m). - real, intent(in) :: zi(ncol,pver+1) ! Interface altitudes above ground (m). - real, intent(in) :: nm(ncol,pver) ! Midpoint Brunt-Vaisalla frequencies (s-1). - real, intent(in) :: ni(ncol,pver+1) ! Interface Brunt-Vaisalla frequencies (s-1). - real, intent(in) :: rhoi(ncol,pver+1) ! Interface density (kg m-3). - real, intent(in) :: kvtt(ncol,pver+1) ! Molecular thermal diffusivity. - - real, intent(in) :: sgh(ncol) ! subgrid orographic std dev (m) - real, intent(in) :: lats(ncol) ! latitudes - real, intent(in) :: alpha(:) - - - real, intent(out) :: utgw(ncol,pver) ! zonal wind tendency - real, intent(out) :: vtgw(ncol,pver) ! meridional wind tendency - real, intent(out) :: ttgw(ncol,pver) ! temperature tendency - - !---------------------------Local storage------------------------------- - - integer :: k, m, nn - - real(GW_PRC), allocatable :: tau(:,:,:) ! wave Reynolds stress - ! gravity wave wind tendency for each wave - real(GW_PRC), allocatable :: gwut(:,:,:) - ! Wave phase speeds for each column - real(GW_PRC), allocatable :: c(:,:) - - ! Efficiency for a gravity wave source. - real :: effgw(ncol) - - ! Indices of top gravity wave source level and lowest level where wind - ! tendencies are allowed. - integer :: src_level(ncol) - integer :: tend_level(ncol) - - ! Projection of wind at midpoints and interfaces. - real :: ubm(ncol,pver) - real :: ubi(ncol,pver+1) - - ! Unit vectors of source wind (zonal and meridional components). - real :: xv(ncol) - real :: yv(ncol) - - real :: pint_adj(ncol,pver+1) - real :: zfac_layer - - character(len=1) :: cn - character(len=9) :: fname(4) - - integer :: i,j - - !---------------------------------------------------------------------------- - - - ! Allocate wavenumber fields. - allocate(tau(ncol, -band%ngwv:band%ngwv,pver+1)) - allocate(gwut(ncol,pver, -band%ngwv:band%ngwv)) - allocate(c(ncol, -band%ngwv:band%ngwv)) - -! Efficiency of gravity wave momentum transfer. - effgw(:) = effgw_oro - - -! Determine the orographic wave source - call gw_oro_src(ncol, pver, band, pint, pmid, delp, & - u, v, t, sgh, zm, nm, & - src_level, tend_level, tau, ubm, ubi, xv, yv, c) - - - do i = 1, ncol - if (lats(i) < 0.) then - tau(i,:,:) = tau(i,:,:) * gw_oro_south_fac - end if - end do - -!GEOS pressure scaling near model top - zfac_layer = 1000.0 ! 10mb - do k=1,pver+1 - do i=1,ncol - pint_adj(i,k) = MIN(1.0,MAX(0.0,(pint(i,k)/zfac_layer)**3)) - enddo - enddo - - ! Solve for the drag profile with orographic sources. - call gw_drag_prof(ncol, pver, band, pint, delp, rdelp, & - src_level, tend_level, dt, t, & - piln, rhoi, nm, ni, ubm, ubi, xv, yv, & - c, kvtt, tau, utgw, vtgw, & - ttgw, gwut, alpha, tau_adjust=pint_adj) - ! Apply efficiency and limiters - call energy_momentum_adjust(ncol, pver, band, pint, delp, u, v, dt, c, tau, & - effgw, t, ubm, ubi, xv, yv, utgw, vtgw, ttgw, tend_level) - -end subroutine gw_oro_ifc - - -!************************************************************************ -!!handle_err -!************************************************************************ -! -!!ROUTINE: handle_err -!!DESCRIPTION: error handler -!-------------------------------------------------------------------------- - -subroutine handle_err(status) - - implicit none - -#include - - integer status - - if (status .ne. nf_noerr) then - print *, nf_strerror(status) - stop 'Stopped' - endif - -end subroutine handle_err - - - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -subroutine endrun(msg) - - integer :: iulog - - character(len=*), intent(in), optional :: msg ! string to be printed - - iulog=6 - - if (present (msg)) then - write(iulog,*)'ENDRUN:', msg - else - write(iulog,*)'ENDRUN: called without a message string' - end if - - stop - -end subroutine endrun - - - - - -! Short routine to get the indices of a set of values rounded to their -! nearest points on a grid. -function index_of_nearest(x, grid) result(idx) - real, intent(in) :: x(:) - real, intent(in) :: grid(:) - - integer :: idx(size(x)) - - real :: interfaces(size(grid)-1) - integer :: i, n - - n = size(grid) - interfaces = (grid(:n-1) + grid(2:))/2. - - idx = 1 - do i = 1, n-1 - where (x > interfaces(i)) idx = i + 1 - end do - -end function index_of_nearest - -end module gw_oro diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/gw_rdg.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/gw_rdg.F90 deleted file mode 100644 index fab1e67c9..000000000 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/gw_rdg.F90 +++ /dev/null @@ -1,1291 +0,0 @@ -module gw_rdg - -! -! This module handles gravity waves from orographic sources, and was -! extracted from gw_drag in May 2013. -! -!use shr_const_mod, only: pi => shr_const_pi -use gw_common, only: gw_drag_prof, GWBand, pi,cpair,rair, & - calc_taucd, momentum_flux, momentum_fixer, & - energy_momentum_adjust, energy_change, energy_fixer -use gw_utils, only: GW_PRC, dot_2d, midpoint_interp - - -implicit none -private - -! Public interface(s) -public :: gw_rdg_init -public :: gw_rdg_ifc - -! Tunable Parameters -!-------------------- -logical :: do_divstream - -!=========================================== -! Parameters for DS2017 (do_divstream=.T.) -!=========================================== -! Amplification factor - 1.0 for -! high-drag/windstorm regime -real, protected :: C_BetaMax_DS - -! Max Ratio Fr2:Fr1 - 1.0 -real, protected :: C_GammaMax - -! Normalized limits for Fr2(Frx) function -real, protected :: Frx0 -real, protected :: Frx1 - - -!=========================================== -! Parameters for SM2000 -!=========================================== -! Amplification factor - 1.0 for -! high-drag/windstorm regime -real, protected :: C_BetaMax_SM - - - -! NOTE: Critical inverse Froude number Fr_c is -! 1./(SQRT(2.)~0.707 in SM2000 -! (should be <= 1) -real, protected :: Fr_c - - -logical :: do_smooth_regimes -logical :: do_adjust_tauoro -logical :: do_backward_compat - - -! Limiters (min/max values) -! min surface displacement height for orographic waves (m) -real, protected :: orohmin -! min wind speed for orographic waves -real, protected :: orovmin -! min stratification allowing wave behavior -real, protected :: orostratmin -! min stratification allowing wave behavior -real, protected :: orom2min -! tendency limiter for ridge scheme -real, protected :: orotndmax - -real, parameter :: PINTADJ_0 = 0.9/19.0 -real, parameter :: PINTADJ_1 = TAN(20.*pi/21.-0.5*pi) -real, parameter :: PINTADJ_2 = 0.5*pi -real, parameter :: PINTADJ_3 = 21./pi - -! Some description of GW spectrum -!type(GWBand) :: band ! I hate this variable ... it just hides information from view - - -!========================================================================== -contains -!========================================================================== - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -!! CCPP Interface routines -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - -!------------------------ -!------------------------------------ -!> \section arg_table_gw_rdg_init Argument Table -!! \htmlinclude gw_rdg_init.html -subroutine gw_rdg_init (band, gw_dc, fcrit2, wavelength, tndmax, pgwv) -#include - - type(GWBand), intent(inout) :: band ! I hate this variable ... it just hides information from view - real, intent(in) :: gw_dc,fcrit2,wavelength,tndmax - integer, intent(in) :: pgwv - - !============================================== - ! Create "Band" structure - !---------------------------------------------- - - band = GWBand(pgwv, gw_dc, fcrit2, wavelength ) - - ! Set the local variables - do_divstream = .TRUE. - C_BetaMax_DS = 0. - C_GammaMax = 2. - Frx0 = 2. - Frx1 = 3. - C_BetaMax_SM = 2. - Fr_c = fcrit2 - do_smooth_regimes = .FALSE. - do_adjust_tauoro = .TRUE. - do_backward_compat = .FALSE. - orohmin = 0.01 - orovmin = 1.0e-3 - orostratmin = 0.002 - orom2min = 0.1 - orotndmax = tndmax - -end subroutine gw_rdg_init - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - -!------------------------------------ -!> \section arg_table_gw_rdg_ifc Argument Table -!! \htmlinclude gw_rdg_ifc.html -subroutine gw_rdg_ifc( band, & - ncol, pver, pverp, pcnst, n_rdg, dt, & - u, v, t, pint, pmid, delp, rdelp, & - piln, zm, zi, z, & - ni, nm, rhoi, & - kvtt, & - kwvrdg, effrdg, & - hwdth, clngt, gbxar, & - mxdis, angll, anixy, & - rdg_cd_llb, trpd_leewv, alpha, & - utrdg, vtrdg, ttrdg, & - flx_heat ) - - !!character(len=5), intent(in) :: type ! BETA or GAMMA - type(GWBand), intent(in) :: band ! I hate this variable ... it just hides information from view - integer, intent(in) :: ncol ! number of atmospheric columns - integer, intent(in) :: pverp ! Layer Vertical dimension - integer, intent(in) :: pver ! Intfc Vertical dimension - integer, intent(in) :: pcnst ! constituent dimension - integer, intent(in) :: n_rdg - real, intent(in) :: dt ! Time step. - - real, intent(in) :: u(ncol,pver) ! Midpoint zonal winds. ( m s-1) - real, intent(in) :: v(ncol,pver) ! Midpoint meridional winds. ( m s-1) - real, intent(in) :: t(ncol,pver) ! Midpoint temperatures. (K) - real, intent(in) :: delp(ncol,pver) ! Delta(interface pressures). - real, intent(in) :: rdelp(ncol,pver) ! inverse Delta. - real, intent(in) :: pmid(ncol,pver) ! midpoint pressures. - real, intent(in) :: pint(ncol,pverp) ! interface pressures. - real, intent(in) :: piln(ncol,pverp) ! Log of interface pressures. - real, intent(in) :: zm(ncol,pver) ! Midpoint altitudes above ground (m). - real, intent(in) :: zi(ncol,pverp) ! Interface altitudes above ground (m). - real, intent(in) :: z(ncol) ! Surface elevation (m). - real, intent(in) :: kvtt(ncol,pverp) ! Molecular thermal diffusivity. - !!real, intent(in) :: q(ncol,pver,pcnst) ! Constituent array. - !!real, intent(in) :: dse(ncol,pver) ! Dry static energy. - - real, intent(in) :: nm(ncol,pver) ! Midpoint altitudes above ground (m). - real, intent(in) :: ni(ncol,pverp) ! Interface altitudes above ground (m). - real, intent(in) :: rhoi(ncol,pverp) ! Interface altitudes above ground (m). - - real, intent(in) :: kwvrdg(ncol,n_rdg) ! horiz wavenumber. - real, intent(in) :: effrdg(ncol,n_rdg) ! efficiency factor of ridge scheme. - real, intent(in) :: hwdth(ncol,n_rdg) ! width of ridges. - real, intent(in) :: clngt(ncol,n_rdg) ! length of ridges. - real, intent(in) :: gbxar(ncol) ! gridbox area - - real, intent(in) :: mxdis(ncol,n_rdg) ! Height estimate for ridge (m). - real, intent(in) :: angll(ncol,n_rdg) ! orientation of ridges. - real, intent(in) :: anixy(ncol,n_rdg) ! Anisotropy parameter. - - real, intent(in) :: rdg_cd_llb ! Drag coefficient for low-level flow - logical, intent(in) :: trpd_leewv - real, intent(in) :: alpha(:) - - - ! OUTPUTS - real, intent(out) :: utrdg(ncol,pver) ! Cum. zonal wind tendency - real, intent(out) :: vtrdg(ncol,pver) ! Cum. meridional wind tendency - real, intent(out) :: ttrdg(ncol,pver) ! Cum. temperature tendency - !!real, intent(out) :: qtrdg(ncol,pver,pcnst) ! Cum. consituent tendencies - real, intent(inout) :: flx_heat(ncol) ! Energy change - - !---------------------------Local storage------------------------------- - - integer :: i, k, m, nn, icnst - - real(GW_PRC), allocatable :: tau(:,:,:) ! wave Reynolds stress - ! gravity wave wind tendency for each wave - real(GW_PRC), allocatable :: gwut(:,:,:) - ! Wave phase speeds for each column - real(GW_PRC), allocatable :: c(:,:) - - ! Isotropic source flag [anisotropic orography]. - integer :: isoflag(ncol) - - ! Indices of top gravity wave source level and lowest level where wind - ! tendencies are allowed. - integer :: src_level(ncol) - integer :: tend_level(ncol) - integer :: bwv_level(ncol) - integer :: tlb_level(ncol) - - ! Projection of wind at midpoints and interfaces. - real :: ubm(ncol,pver) - real :: ubi(ncol,pverp) - - ! Unit vectors of source wind (zonal and meridional components). - real :: xv(ncol) - real :: yv(ncol) - - ! Averages over source region. - real :: ubmsrc(ncol) ! On-ridge wind. - real :: usrc(ncol) ! Zonal wind. - real :: vsrc(ncol) ! Meridional wind. - real :: nsrc(ncol) ! B-V frequency. - real :: rsrc(ncol) ! Density. - - ! normalized wavenumber - real :: m2src(ncol) - - ! Top of low-level flow layer. - real :: tlb(ncol) - - ! Bottom of linear wave region. - real :: bwv(ncol) - - ! Froude numbers for flow/drag regimes - real :: Fr1(ncol) - real :: Fr2(ncol) - real :: Frx(ncol) - - ! Wave Reynolds stresses at source level - real :: tauoro(ncol) - real :: taudsw(ncol) - - ! Surface streamline displacement height for linear waves. - real :: hdspwv(ncol) - - ! Surface streamline displacement height for downslope wind regime. - real :: hdspdw(ncol) - - ! Wave breaking level - real :: wbr(ncol) - - ! Momentum fluxes used by fixer. - real :: um_flux(ncol), vm_flux(ncol) - - ! Energy change used by fixer. - real :: de(ncol) - - ! Reynolds stress for waves propagating in each cardinal direction. - real :: taucd(ncol,pver+1,4) - - real :: utgw(ncol,pver) ! zonal wind tendency - real :: vtgw(ncol,pver) ! meridional wind tendency - real :: ttgw(ncol,pver) ! temperature tendency -#ifdef CAM - real :: qtgw(ncol,pver,pcnst) ! constituents tendencies -#endif - - real :: pint_adj(ncol,pver+1) - real :: zfac_layer - real :: utfac,uhtmax - - character(len=4) :: type ! BETA or GAMMA (just BETA for now) - character(len=1) :: cn - character(len=9) :: fname(4) - !---------------------------------------------------------------------------- - - ! Allocate wavenumber fields. - allocate(tau(ncol, -band%ngwv:band%ngwv , pverp)) - allocate(gwut(ncol,pver,-band%ngwv:band%ngwv )) - allocate(c(ncol,-band%ngwv:band%ngwv)) - - type='BETA' - - ! initialize accumulated momentum fluxes and tendencies - utrdg = 0. - vtrdg = 0. - ttrdg = 0. - -!GEOS pressure scaling near model top - zfac_layer = 1000.0 ! 10mb - do k=1,pver+1 - do i=1,ncol - pint_adj(i,k) = MIN(1.0,MAX(0.0,(pint(i,k)/zfac_layer)**3)) - enddo - enddo - - isoflag = 0 - - do nn = 1, n_rdg - - call gw_rdg_src(band, ncol, pver, pint, pmid, delp, & - u, v, t, mxdis(:,nn), angll(:,nn), anixy(:,nn), kwvrdg(:,nn), isoflag, zi, nm, & - src_level, tend_level, bwv_level, tlb_level, tau, ubm, ubi, xv, yv, & - ubmsrc, usrc, vsrc, nsrc, rsrc, m2src, tlb, bwv, Fr1, Fr2, Frx, c) - - call gw_rdg_belowpeak(band, ncol, pver, rdg_cd_llb, & - t, mxdis(:,nn), anixy(:,nn), kwvrdg(:,nn), & - zi, nm, ni, rhoi, & - src_level, tau, & - ubmsrc, nsrc, rsrc, m2src, tlb, bwv, Fr1, Fr2, Frx, & - tauoro, taudsw, hdspwv, hdspdw) - - call gw_rdg_break_trap(band, ncol, pver, & - zi, nm, ni, ubm, ubi, rhoi, kwvrdg(:,nn) , bwv, tlb, wbr, & - src_level, tlb_level, hdspwv, hdspdw, mxdis(:,nn), & - tauoro, taudsw, tau, & - ldo_trapped_waves=trpd_leewv) - - call gw_drag_prof(ncol, pver, band, pint, delp, rdelp, & - src_level, tend_level,dt, t, & - piln, rhoi, nm, ni, ubm, ubi, xv, yv, & - c, kvtt, tau, utgw, vtgw, & - ttgw, gwut, alpha, & - kwvrdg=kwvrdg(:,nn), satfac_in=1.0, tau_adjust=pint_adj) - - ! Apply efficiency and limiters to the totals - call energy_momentum_adjust(ncol, pver, band, pint, delp, u, v, dt, c, tau, & - effrdg(:,nn), t, ubm, ubi, xv, yv, utgw, vtgw, ttgw, & - tend_level, tndmax_in=orotndmax) - - do i=1,ncol - !------------------------------------------------------------------- - ! Apply tendency limiter to prevent unrealistically strong forcing - !------------------------------------------------------------------- - uhtmax = 0.0 - utfac = 1.0 - do k = 1, pver - ! Add the tendencies from each ridge to the totals. - utrdg(i,k) = utrdg(i,k) + utgw(i,k) - vtrdg(i,k) = vtrdg(i,k) + vtgw(i,k) - ttrdg(i,k) = ttrdg(i,k) + ttgw(i,k) - uhtmax = max(sqrt(utrdg(i,k)**2 + vtrdg(i,k)**2), uhtmax) - end do - if (uhtmax > orotndmax) utfac = orotndmax/uhtmax - do k = 1, pver - utrdg(i,k) = utrdg(i,k)*utfac - vtrdg(i,k) = vtrdg(i,k)*utfac - ttrdg(i,k) = ttrdg(i,k)*utfac - end do - end do ! i=1,ncol - -#ifdef CAM -! disable tracer mixing in GW for now. - do icnst = 1, pcnst - do k = 1, pver - qtrdg(:,k,icnst) = qtrdg(:,k,icnst) + qtgw(:,k,icnst) - end do - end do - if (nn <= 6) then - write(cn, '(i1)') nn - end if -#endif - - end do ! end of loop over multiple ridges - - if (trim(type) == 'BETA') then - fname(1) = 'TAUGWX' - fname(2) = 'TAUGWY' - fname(3) = 'UTGWORO' - fname(4) = 'VTGWORO' - else if (trim(type) == 'GAMMA') then - fname(1) = 'TAURDGGMX' - fname(2) = 'TAURDGGMY' - fname(3) = 'UTRDGGM' - fname(4) = 'VTRDGGM' - else - call endrun('gw_rdg_calc: FATAL: type must be either BETA or GAMMA'& - //' type= '//type) - end if - - deallocate(tau, gwut, c) - - end subroutine gw_rdg_ifc - - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! Non - interface subroutines -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - -!------------------------ -subroutine gw_rdg_src(band, ncol, pver , pint, pmid, delp, & - u, v, t, mxdis, angxy, anixy, kwvrdg, iso, zi, nm, & - src_level, tend_level, bwv_level ,tlb_level , tau, ubm, ubi, xv, yv, & - ubmsrc, usrc, vsrc, nsrc, rsrc, m2src, tlb, bwv, Fr1, Fr2, Frx, c) - - !----------------------------------------------------------------------- - ! Orographic source for multiple gravity wave drag parameterization. - ! - ! The stress is returned for a single wave with c=0, over orography. - ! For points where the orographic variance is small (including ocean), - ! the returned stress is zero. - !------------------------------Arguments-------------------------------- - type(GWBand), intent(in) :: band ! I hate this variable ... it just hides information from view - ! Column dimension. - integer, intent(in) :: ncol - ! Vertical dimension. - integer, intent(in) :: pver - - ! Interface pressures. (Pa) - real, intent(in) :: pint(ncol,pver+1) - ! Midpoint pressures. (Pa) - real, intent(in) :: pmid(ncol,pver) - ! Delta Interface pressures. (Pa) - real, intent(in) :: delp(ncol,pver) - - - ! Midpoint zonal/meridional winds. ( m s-1) - real, intent(in) :: u(ncol,pver), v(ncol,pver) - ! Midpoint temperatures. (K) - real, intent(in) :: t(ncol,pver) - ! Height estimate for ridge (m) [anisotropic orography]. - real, intent(in) :: mxdis(ncol) - ! Angle of ridge axis w/resp to north (degrees) [anisotropic orography]. - real, intent(in) :: angxy(ncol) - ! Anisotropy parameter [anisotropic orography]. - real, intent(in) :: anixy(ncol) - ! horiz wavenumber [anisotropic orography]. - real, intent(in) :: kwvrdg(ncol) - ! Isotropic source flag [anisotropic orography]. - integer, intent(in) :: iso(ncol) - ! Interface altitudes above ground (m). - real, intent(in) :: zi(ncol,pver+1) - ! Midpoint Brunt-Vaisalla frequencies (s-1). - real, intent(in) :: nm(ncol,pver) - - ! Indices of top gravity wave source level and lowest level where wind - ! tendencies are allowed. - integer, intent(out) :: src_level(ncol) - integer, intent(out) :: tend_level(ncol) - integer, intent(out) :: bwv_level(ncol),tlb_level(ncol) - - ! Averages over source region. - real, intent(out) :: nsrc(ncol) ! B-V frequency. - real, intent(out) :: rsrc(ncol) ! Density. - real, intent(out) :: usrc(ncol) ! Zonal wind. - real, intent(out) :: vsrc(ncol) ! Meridional wind. - real, intent(out) :: ubmsrc(ncol) ! On-ridge wind. - ! Top of low-level flow layer. - real, intent(out) :: tlb(ncol) - ! Bottom of linear wave region. - real, intent(out) :: bwv(ncol) - ! normalized wavenumber - real, intent(out) :: m2src(ncol) - - - ! Wave Reynolds stress. - real(GW_PRC), intent(out) :: tau(ncol,-band%ngwv:band%ngwv,pver+1) - ! Projection of wind at midpoints and interfaces. - real, intent(out) :: ubm(ncol,pver), ubi(ncol,pver+1) - ! Unit vectors of source wind (zonal and meridional components). - real, intent(out) :: xv(ncol), yv(ncol) - ! Phase speeds. - real(GW_PRC), intent(out) :: c(ncol,-band%ngwv:band%ngwv) - ! Froude numbers for flow/drag regimes - real, intent(out) :: Fr1(ncol), Fr2(ncol), Frx(ncol) - - !---------------------------Local Storage------------------------------- - ! Column and level indices. - integer :: i, k - - ! Surface streamline displacement height (2*sgh). - real :: hdsp(ncol) - - ! Difference in interface pressure across source region. - real :: dpsrc(ncol) - ! Thickness of downslope wind region. - real :: ddw(ncol) - ! Thickness of linear wave region. - real :: dwv(ncol) - ! Wind speed in source region. - real :: wmsrc(ncol) - - real :: ragl(ncol) - -!-------------------------------------------------------------------------- -! Check that ngwav is equal to zero, otherwise end the job -!-------------------------------------------------------------------------- - !! if (band%ngwv /= 0) call endrun(' gw_rdg_src :: ERROR - band%ngwv must be zero and it is not') - -!-------------------------------------------------------------------------- -! Average the basic state variables for the wave source over the depth of -! the orographic standard deviation. Here we assume that the appropiate -! values of wind, stability, etc. for determining the wave source are -! averages over the depth of the atmosphere penterated by the typical -! mountain. -! Reduces to the bottom midpoint values when mxdis=0, such as over ocean. -!-------------------------------------------------------------------------- - - hdsp = mxdis ! no longer multipied by 2 - - tau(:,0,:) = 0.0 - - ! Find min height to begin tendency forcing - tend_level = pver - do k = pver, 1, -1 - do i = 1, ncol - if ( 0.5*(zi(i,k+1)+zi(i,k)) <= orohmin ) then - tend_level(i) = k - end if - end do - end do - ! Find depth of "source layer" for mountain waves - ! i.e., between ground and mountain top - src_level = tend_level - do k = pver, 1, -1 - do i = 1, ncol - ! Need to have h >= z(k+1) here or code will bomb when h=0. - if ( (k <= tend_level(i)) .and. (hdsp(i) >= zi(i,k+1)) .and. (hdsp(i) < zi(i,k)) ) then - src_level(i) = k - end if - end do - end do - - rsrc = 0. - usrc = 0. - vsrc = 0. - nsrc = 0. - do i = 1, ncol - do k = pver, src_level(i), -1 - rsrc(i) = rsrc(i) + pmid(i,k) / ( rair * t(i,k))* delp(i,k) - usrc(i) = usrc(i) + u(i,k) * delp(i,k) - vsrc(i) = vsrc(i) + v(i,k) * delp(i,k) - nsrc(i) = nsrc(i) + nm(i,k)* delp(i,k) - end do - end do - - do i = 1, ncol - dpsrc(i) = pint(i,tend_level(i)) - pint(i,src_level(i)) - end do - where ( dpsrc > 0.0 ) - rsrc = rsrc / dpsrc - usrc = usrc / dpsrc - vsrc = vsrc / dpsrc - nsrc = nsrc / dpsrc - end where - - wmsrc = sqrt( usrc**2 + vsrc**2 ) - - - ! Get the unit vector components - ! Want agl=0 with U>0 to give xv=1 - - ragl = angxy * pi/180. - - ! protect from wierd "bad" angles - ! that may occur if hdsp is zero - where( hdsp <= orohmin ) - ragl = 0. - end where - - yv =-sin( ragl ) - xv = cos( ragl ) - - - ! Kluge in possible "isotropic" obstacle. - where( ( iso == 1 ) .and. (wmsrc > orovmin) ) - xv = usrc/wmsrc - yv = vsrc/wmsrc - end where - - - ! Project the local wind at midpoints into the on-ridge direction - do k = 1, pver - ubm(:,k) = dot_2d(u(:,k), v(:,k), xv, yv) - end do - ubmsrc = dot_2d(usrc , vsrc , xv, yv) - - ! Ensure on-ridge wind is positive at source level - do k = 1, pver - ubm(:,k) = sign( ubmsrc*0.+1. , ubmsrc ) * ubm(:,k) - end do - - ! Sean says just use 1. as - ! first argument - xv = sign( ubmsrc*0.+1. , ubmsrc ) * xv - yv = sign( ubmsrc*0.+1. , ubmsrc ) * yv - - ! Now make ubmsrc positive and protect - ! against zero - ubmsrc = abs(ubmsrc) - ubmsrc = max( 0.01 , ubmsrc ) - - - ! The minimum stratification allowing GW behavior - ! should really depend on horizontal scale since - ! - ! m^2 ~ (N/U)^2 - k^2 - ! - ! Should also think about parameterizing - ! trapped lee-waves. - - - ! This needs to be made constistent with later - ! treatment of nonhydrostatic effects. - m2src = ( (nsrc/(ubmsrc+0.01))**2 - kwvrdg**2 ) /((nsrc/(ubmsrc+0.01))**2) - - - !------------------------------------------------------------- - ! Calculate provisional limits (in Z [m]) for 3 regimes. This - ! will modified later if wave breaking or trapping are - ! diagnosed - ! - ! ^ - ! | *** linear propagation *** - ! (H) -------- mountain top ------------- | *** or wave breaking **** - ! | *** regimes ************* - ! (BWV)------ bottom of linear waves ---- | - ! : | - ! ******* | - ! : | - ! (TLB)--- top of flow diversion layer--- ' - ! : - ! **** flow diversion ***** - ! : - !============================================ - - !============================================ - ! For Dividing streamline para (DS2017) - !-------------------------------------------- - ! High-drag downslope wind regime exists - ! between bottom of linear waves and top of - ! flow diversion. Linear waves can only - ! attain vertical displacment of f1*U/N. So, - ! bottom of linear waves is given by - ! - ! BWV = H - Fr1*U/N - ! - ! Downslope wind layer begins at BWV and - ! extends below it until some maximum high - ! drag obstacle height Fr2*U/N is attained - ! (where Fr2 >= f1). Below downslope wind - ! there is flow diversion, so top of - ! diversion layer (TLB) is equivalent to - ! bottom of downslope wind layer and is; - ! - ! TLB = H - Fr2*U/N - ! - !----------------------------------------- - - ! Critical inverse Froude number - !----------------------------------------------- - Fr1(:) = Fr_c * 1.00 - Frx(:) = hdsp(:)*nsrc(:)/abs( ubmsrc(:) ) / Fr_c - - if ( do_divstream ) then - !------------------------------------------------ - ! Calculate Fr2(Frx) for DS2017 - !------------------------------------------------ - where(Frx <= Frx0) - Fr2(:) = Fr1(:) + Fr1(:)* C_GammaMax * anixy(:) - elsewhere((Frx > Frx0).and.(Frx <= Frx1) ) - Fr2(:) = Fr1(:) + Fr1(:)* C_GammaMax * anixy(:) & - * (Frx1 - Frx(:))/(Frx1-Frx0) - elsewhere(Frx > Frx1) - Fr2(:)=Fr1(:) - endwhere - else - !------------------------------------------ - ! Regime distinctions entirely carried by - ! amplification of taudsw (next subr) - !------------------------------------------ - Fr2(:)=Fr1(:) - end if - - - - where( m2src > orom2min ) - ddw = Fr2 * ( abs(ubmsrc) )/nsrc - elsewhere - ddw = 0. - endwhere - - - ! If TLB is less than zero then obstacle is not - ! high enough to produce an low-level diversion layer - tlb = mxdis - ddw - where( tlb < 0.) - tlb = 0. - endwhere - tlb_level = tend_level - do k = pver, pver/2, -1 - do i = 1, ncol - if ( (k <= tend_level(i)) .and. (tlb(i) > zi(i,k+1)) .and. (tlb(i) <= zi(i,k)) ) then - tlb_level(i) = k - end if - end do - end do - - - ! Find *BOTTOM* of linear wave layer (BWV) - !where ( nsrc > orostratmin ) - where( m2src > orom2min ) - dwv = Fr1 * ( abs(ubmsrc) )/nsrc - elsewhere - dwv = -9.999e9 ! if weak strat - no waves - endwhere - - bwv = mxdis - dwv - where(( bwv < 0.) .or. (dwv < 0.) ) - bwv = 0. - endwhere - bwv_level = tend_level - do k = pver,1, -1 - do i = 1, ncol - if ( (k <= tend_level(i)) .and. (bwv(i) > zi(i,k+1)) .and. (bwv(i) <= zi(i,k)) ) then - bwv_level(i) = k+1 - end if - end do - end do - - - ! Compute the interface wind projection by averaging the midpoint winds. - ! Use the top level wind at the top interface. - ubi(:,1) = ubm(:,1) - ubi(:,2:pver) = midpoint_interp(ubm) - ubi(:,pver+1) = ubm(:,pver) - - ! No spectrum; phase speed is just 0. - c = 0. - - where( m2src < orom2min ) - tlb = mxdis - tlb_level = src_level - endwhere - - -end subroutine gw_rdg_src - - -!========================================================================== - -subroutine gw_rdg_belowpeak(band, ncol, pver, rdg_cd_llb, & - t, mxdis, anixy, kwvrdg, zi, nm, ni, rhoi, & - src_level , tau, & - ubmsrc, nsrc, rsrc, m2src,tlb,bwv,Fr1,Fr2,Frx, & - tauoro,taudsw, hdspwv,hdspdw ) - - !----------------------------------------------------------------------- - ! Orographic source for multiple gravity wave drag parameterization. - ! - ! The stress is returned for a single wave with c=0, over orography. - ! For points where the orographic variance is small (including ocean), - ! the returned stress is zero. - !------------------------------Arguments-------------------------------- - type(GWBand), intent(in) :: band ! I hate this variable ... it just hides information from view - ! Column dimension. - integer, intent(in) :: ncol - ! Vertical dimension. - integer, intent(in) :: pver - ! Drag coefficient for low-level flow - real, intent(in) :: rdg_cd_llb - - - ! Midpoint temperatures. (K) - real, intent(in) :: t(ncol,pver) - ! Height estimate for ridge (m) [anisotropic orography]. - real, intent(in) :: mxdis(ncol) - ! Anisotropy parameter [0-1] [anisotropic orography]. - real, intent(in) :: anixy(ncol) - ! Inverse cross-ridge lengthscale (m-1) [anisotropic orography]. - real, intent(in) :: kwvrdg(ncol) - ! Interface altitudes above ground (m). - real, intent(in) :: zi(ncol,pver+1) - ! Midpoint Brunt-Vaisalla frequencies (s-1). - real, intent(in) :: nm(ncol,pver) - ! Interface Brunt-Vaisalla frequencies (s-1). - real, intent(in) :: ni(ncol,pver+1) - ! Interface density (kg m-3). - real, intent(in) :: rhoi(ncol,pver+1) - - ! Indices of top gravity wave source level - integer, intent(inout) :: src_level(ncol) - - ! Wave Reynolds stress. - real(GW_PRC), intent(inout) :: tau(ncol,-band%ngwv:band%ngwv,pver+1) - ! Top of low-level flow layer. - real, intent(inout) :: tlb(ncol) - ! Bottom of linear wave region. - real, intent(inout) :: bwv(ncol) - ! surface stress from linear waves. - real, intent(out) :: tauoro(ncol) - ! surface stress for downslope wind regime. - real, intent(out) :: taudsw(ncol) - - ! Surface streamline displacement height for linear waves. - real, intent(out) :: hdspwv(ncol) - ! Surface streamline displacement height for downslope wind regime. - real, intent(out) :: hdspdw(ncol) - - - - ! Froude numbers for flow/drag regimes - real, intent(in) :: Fr1(ncol), Fr2(ncol),Frx(ncol) - - ! Averages over source region. - real, intent(in) :: m2src(ncol) ! normalized non-hydro wavenumber - real, intent(in) :: nsrc(ncol) ! B-V frequency. - real, intent(in) :: rsrc(ncol) ! Density. - real, intent(in) :: ubmsrc(ncol) ! On-ridge wind. - - - !logical, intent(in), optional :: forcetlb - - !---------------------------Local Storage------------------------------- - ! Column and level indices. - integer :: i, k - - real :: Coeff_LB(ncol),tausat,ubsrcx(ncol),dswamp - real :: taulin(ncol),BetaMax - - ! ubsrcx introduced to account for situations with high shear, strong strat. - do i = 1, ncol - ubsrcx(i) = max( ubmsrc(i) , 0. ) - end do - - do i = 1, ncol - if ( m2src(i) > orom2min ) then - hdspwv(i) = min( mxdis(i) , Fr1(i) * ubsrcx(i) / nsrc(i) ) - else - hdspwv(i) = 0. - end if - end do - - if (do_divstream) then - do i = 1, ncol - if ( m2src(i) > orom2min ) then - hdspdw(i) = min( mxdis(i) , Fr2(i) * ubsrcx(i) / nsrc(i) ) - else - hdspdw(i) = 0. - end if - end do - else - do i = 1, ncol - ! Needed only to mark where a DSW occurs - if ( m2src(i) > orom2min ) then - hdspdw(i) = mxdis(i) - else - hdspdw(i) = 0. - end if - end do - end if - - ! Calculate form drag coefficient ("CD") - !-------------------------------------- - Coeff_LB = rdg_cd_llb*anixy - - ! Determine the orographic c=0 source term following McFarlane (1987). - ! Set the source top interface index to pver, if the orographic term is - ! zero. - ! - ! This formula is basically from - ! - ! tau(src) = rho * u' * w' - ! where - ! u' ~ N*h' and w' ~ U*h'/b (b="breite") - ! - ! and 1/b has been replaced with k (kwvrdg) - ! - do i = 1, ncol - if ( ( src_level(i) > 0 ) .and. ( m2src(i) > orom2min ) ) then - tauoro(i) = kwvrdg(i) * ( hdspwv(i)**2 ) * rsrc(i) * nsrc(i) & - * ubsrcx(i) - taudsw(i) = kwvrdg(i) * ( hdspdw(i)**2 ) * rsrc(i) * nsrc(i) & - * ubsrcx(i) - else - tauoro(i) = 0. - taudsw(i) = 0. - end if - end do - - if (do_divstream) then - do i = 1, ncol - taulin(i) = 0. - end do - !--------------------------------------- - ! Need linear drag when divstream is not used - !--------------------------------------- - else - do i = 1, ncol - if ( ( src_level(i) > 0 ) .and. ( m2src(i) > orom2min ) ) then - taulin(i) = kwvrdg(i) * ( mxdis(i)**2 ) * rsrc(i) * nsrc(i) & - * ubsrcx(i) - else - taulin(i) = 0. - end if - end do - end if - - if ( do_divstream ) then - ! Amplify DSW between Frx=1. and Frx=Frx1 - do i = 1,ncol - dswamp=0. - BetaMax = C_BetaMax_DS * anixy(i) - if ( (Frx(i)>1.).and.(Frx(i)<=Frx1)) then - dswamp = (Frx(i)-1.)*(Frx1-Frx(i))/(0.25*(Frx1-1.)**2) - end if - taudsw(i) = (1. + BetaMax*dswamp)*taudsw(i) - end do - else - !------------------- - ! Scinocca&McFarlane - !-------------------- - do i = 1, ncol - BetaMax = C_BetaMax_SM * anixy(i) - if ( (Frx(i) >=1.) .and. (Frx(i) < 1.5) ) then - dswamp = 2. * BetaMax * (Frx(i) -1.) - else if ( ( Frx(i) >= 1.5 ) .and. (Frx(i) < 3. ) ) then - dswamp = ( 1. + BetaMax - (0.666**2) ) * ( 0.666*(3. - Frx(i) ))**2 & - + ( 1. / Frx(i) )**2 -1. - else - dswamp = 0. - end if - if ( (Frx(i) >=1.) .and. (Frx(i) < 3.) ) then - taudsw(i) = (1. + dswamp )*taulin(i) - tauoro(i) - else - taudsw(i) = 0. - endif - ! This code defines "taudsw" as SUM of freely-propagating - ! waves +DSW enhancement. Different than in SM2000 - taudsw(i) = taudsw(i) + tauoro(i) - end do - !---------------------------------------------------- - end if - - - do i = 1, ncol - if ( m2src(i) > orom2min ) then - where ( ( zi(i,:) < mxdis(i) ) .and. ( zi(i,:) >= bwv(i) ) ) - tau(i,0,:) = tauoro(i) - else where ( ( zi(i,:) < bwv(i) ) .and. ( zi(i,:) >= tlb(i) ) ) - tau(i,0,:) = tauoro(i) +( taudsw(i)-tauoro(i) )* & - ( bwv(i) - zi(i,:) ) / & - ( bwv(i) - tlb(i) ) - endwhere - ! low-level form drag on obstacle. Quantity kwvrdg (~1/b) appears for consistency - ! with tauoro and taudsw forms. Should be weighted by L*b/A_g before applied to flow. - where ( ( zi(i,:) < tlb(i) ) .and. ( zi(i,:) >= 0. ) ) - tau(i,0,:) = taudsw(i) + & - Coeff_LB(i) * kwvrdg(i) * rsrc(i) * 0.5 * (ubsrcx(i)**2) * ( tlb(i) - zi(i,:) ) - endwhere - - if (do_smooth_regimes) then - ! This blocks accounts for case where both mxdis and tlb fall - ! between adjacent edges - do k=1,pver - if ( (zi(i,k) >= tlb(i)).and.(zi(i,k+1) < tlb(i)).and. & - (zi(i,k) >= mxdis(i)).and.(zi(i,k+1) < mxdis(i)) ) then - src_level(i) = src_level(i)-1 - tau(i,0,k) = tauoro(i) - end if - end do - end if - - else !---------------------------------------------- - ! This block allows low-level dynamics to occur - ! even if m2 is less than orom2min - where ( ( zi(i,:) < tlb(i) ) .and. ( zi(i,:) >= 0. ) ) - tau(i,0,:) = taudsw(i) + & - Coeff_LB(i) * kwvrdg(i) * rsrc(i) * 0.5 * & - (ubsrcx(i)**2) * ( tlb(i) - zi(i,:) ) - endwhere - endif - end do - - ! This may be redundant with newest version of gw_drag_prof. - ! That code reaches down to level k=src_level+1. (jtb 1/5/16) - do i = 1, ncol - k=src_level(i) - if ( ni(i,k) > orostratmin ) then - tausat = (Fr_c**2) * kwvrdg(i) * rhoi(i,k) * ubsrcx(i)**3 / & - (1.*ni(i,k)) - else - tausat = 0. - endif - tau(i,0,src_level(i)) = min( tauoro(i), tausat ) - end do - - - - ! Final clean-up. Do nothing if obstacle less than orohmin - do i = 1, ncol - if ( mxdis(i) < orohmin ) then - tau(i,0,:) = 0. - tauoro(i) = 0. - taudsw(i) = 0. - endif - end do - - ! Disable vertical propagation if Scorer param is - ! too small. - do i = 1, ncol - if ( m2src(i) <= orom2min ) then - src_level(i)=1 - endif - end do - - - -end subroutine gw_rdg_belowpeak - -!========================================================================== -subroutine gw_rdg_break_trap(band, ncol, pver, & - zi, nm, ni, ubm, ubi, rhoi, kwvrdg, bwv, tlb, wbr, & - src_level, tlb_level, & - hdspwv, hdspdw, mxdis, & - tauoro, taudsw, tau, & - ldo_trapped_waves, wdth_kwv_scale_in ) - !----------------------------------------------------------------------- - ! Parameterization of high-drag regimes and trapped lee-waves for CAM - ! - !------------------------------Arguments-------------------------------- - ! Column dimension. - type(GWBand), intent(in) :: band ! I hate this variable ... it just hides information from view - integer, intent(in) :: ncol - ! Vertical dimension. - integer, intent(in) :: pver - - ! Horz wavenumber for ridge (1/m) [anisotropic orography]. - real, intent(in) :: kwvrdg(ncol) - ! Interface altitudes above ground (m). - real, intent(in) :: zi(ncol,pver+1) - ! Midpoint Brunt-Vaisalla frequencies (s-1). - real, intent(in) :: nm(ncol,pver) - ! Interface Brunt-Vaisalla frequencies (s-1). - real, intent(in) :: ni(ncol,pver+1) - - ! Indices of gravity wave sources. - integer, intent(inout) :: src_level(ncol), tlb_level(ncol) - - ! Wave Reynolds stress. - real(GW_PRC), intent(inout) :: tau(ncol,-band%ngwv:band%ngwv,pver+1) - ! Wave Reynolds stresses at source. - real, intent(inout) :: taudsw(ncol),tauoro(ncol) - ! Projection of wind at midpoints and interfaces. - real, intent(in) :: ubm(ncol,pver) - real, intent(in) :: ubi(ncol,pver+1) - ! Interface density (kg m-3). - real, intent(in) :: rhoi(ncol,pver+1) - - ! Top of low-level flow layer. - real, intent(in) :: tlb(ncol) - ! Bottom of linear wave region. - real, intent(in) :: bwv(ncol) - - ! Surface streamline displacement height for linear waves. - real, intent(in) :: hdspwv(ncol) - ! Surface streamline displacement height for downslope wind regime. - real, intent(in) :: hdspdw(ncol) - ! Ridge height. - real, intent(in) :: mxdis(ncol) - - - ! Wave breaking level - real, intent(out) :: wbr(ncol) - - logical, intent(in), optional :: ldo_trapped_waves - real, intent(in), optional :: wdth_kwv_scale_in - - !---------------------------Local Storage------------------------------- - ! Column and level indices. - integer :: i, k, kp1, non_hydro - real:: m2(ncol,pver),delz(ncol),tausat(ncol),trn(ncol) - real:: wbrx(ncol) - real:: phswkb(ncol,pver+1) - logical :: lldo_trapped_waves - real:: wdth_kwv_scale - ! Indices of important levels. - integer :: trn_level(ncol) - - if (present(ldo_trapped_waves)) then - lldo_trapped_waves = ldo_trapped_waves - if(lldo_trapped_waves) then - non_hydro = 1 - else - non_hydro = 0 - endif - else - lldo_trapped_waves = .false. - non_hydro = 0 - endif - - if (present(wdth_kwv_scale_in)) then - wdth_kwv_scale = wdth_kwv_scale_in - else - wdth_kwv_scale = 1. - endif - - ! Calculate vertical wavenumber**2 - !--------------------------------- - m2 = (nm / (abs(ubm)+.01))**2 - do k=pver,1,-1 - m2(:,k) = m2(:,k) - non_hydro*(wdth_kwv_scale*kwvrdg)**2 - ! sweeping up, zero out m2 above first occurence - ! of m2(:,k)<=0 - kp1=min( k+1, pver ) - where( (m2(:,k) <= 0.0 ).or.(m2(:,kp1) <= 0.0 ) ) - m2(:,k) = 0. - endwhere - end do - - ! Take square root of m**2 and - ! do vertical integral to find - ! WKB phase. - !----------------------------- - m2 = SQRT( m2 ) - phswkb(:,:)=0 - do k=pver,1,-1 - where( zi(:,k) > tlb(:) ) - delz(:) = min( zi(:,k)-zi(:,k+1) , zi(:,k)-tlb(:) ) - phswkb(:,k) = phswkb(:,k+1) + m2(:,k)*delz(:) - endwhere - end do - - ! Identify top edge of layer in which phswkb reaches 3*pi/2 - ! - approximately the "breaking level" - !---------------------------------------------------------- - wbr(:)=0. - wbrx(:)=0. - if (do_smooth_regimes) then - do k=pver,1,-1 - where( (phswkb(:,k+1)<1.5*pi).and.(phswkb(:,k)>=1.5*pi) & - .and.(hdspdw(:)>hdspwv(:)) ) - wbr(:) = zi(:,k) - ! Extrapolation to make regime - ! transitions smoother - wbrx(:) = zi(:,k) - ( phswkb(:,k) - 1.5*pi ) & - / ( m2(:,k) + 1.e-6 ) - src_level(:) = k-1 - endwhere - end do - else - do k=pver,1,-1 - where( (phswkb(:,k+1)<1.5*pi).and.(phswkb(:,k)>=1.5*pi) & - .and.(hdspdw(:)>hdspwv(:)) ) - wbr(:) = zi(:,k) - src_level(:) = k - endwhere - end do - end if - - ! Adjust tauoro at new source levels if needed. - ! This is problematic if Fr_c<1.0. Not sure why. - !---------------------------------------------------------- - if (do_adjust_tauoro) then - do i = 1,ncol - if (wbr(i) > 0. ) then - tausat(i) = (Fr_c**2) * kwvrdg(i) * rhoi( i, src_level(i) ) & - * abs(ubi(i , src_level(i) ))**3 & - / ni( i , src_level(i) ) - tauoro(i) = min( tauoro(i), tausat(i) ) - end if - end do - end if - - if (do_smooth_regimes) then - do i = 1, ncol - do k=1,pver+1 - if ( ( zi(i,k) <= wbr(i) ) .and. ( zi(i,k) > tlb(i) ) ) then - tau(i,0,k) = tauoro(i) + (taudsw(i)-tauoro(i)) * & - ( wbrx(i) - zi(i,k) ) / & - ( wbrx(i) - tlb(i) ) - tau(i,0,k) = max( tau(i,0,k), tauoro(i) ) - endif - end do - end do - else - ! Following is for backwards B4B compatibility with earlier versions - ! ("N1" and "N5" -- Note: "N5" used do_backward_compat=.true.) - if (.not.do_backward_compat) then - do i = 1, ncol - do k=1,pver+1 - if ( ( zi(i,k) < wbr(i) ) .and. ( zi(i,k) >= tlb(i) ) ) then - tau(i,0,k) = tauoro(i) + (taudsw(i)-tauoro(i)) * & - ( wbr(i) - zi(i,k) ) / & - ( wbr(i) - tlb(i) ) - endif - end do - end do - else - do i = 1, ncol - do k=1,pver+1 - if ( ( zi(i,k) <= wbr(i) ) .and. ( zi(i,k) > tlb(i) ) ) then - tau(i,0,k) = tauoro(i) + (taudsw(i)-tauoro(i)) * & - ( wbr(i) - zi(i,k) ) / & - ( wbr(i) - tlb(i) ) - endif - end do - end do - end if - end if - - if (lldo_trapped_waves) then - - ! Identify top edge of layer in which Scorer param drops below 0 - ! - approximately the "turning level" - !---------------------------------------------------------- - trn(:)=1.e8 - trn_level(:) = 0 ! pver+1 - where( m2(:,pver)<= 0. ) - trn(:) = zi(:,pver) - trn_level(:) = pver - endwhere - do k=pver-1,1,-1 - where( (m2(:,k+1)> 0.).and.(m2(:,k)<= 0.) ) - trn(:) = zi(:,k) - trn_level(:) = k - endwhere - end do - - do i = 1,ncol - ! Case: Turning below mountain top - if ( (trn(i) < mxdis(i)).and.(trn_level(i)>=1) ) then - tau(i,0,:) = tau(i,0,:) - max( tauoro(i),taudsw(i) ) - tau(i,0,:) = max( tau(i,0,:) , 0. ) - tau(i,0,1:tlb_level(i))=0. - src_level(i) = 1 ! disable any more tau calculation - end if - ! Case: Turning but no breaking - if ( (wbr(i) == 0. ).and.(trn(i)>mxdis(i)).and.(trn_level(i)>=1) ) then - where ( ( zi(i,:) <= trn(i) ) .and. ( zi(i,:) >= bwv(i) ) ) - tau(i,0,:) = tauoro(i) * & - ( trn(i) - zi(i,:) ) / & - ( trn(i) - bwv(i) ) - end where - src_level(i) = 1 ! disable any more tau calculation - end if - ! Case: Turning AND breaking. Turning ABOVE breaking - if ( (wbr(i) > 0. ).and.(trn(i) >= wbr(i)).and.(trn_level(i)>=1) ) then - where ( ( zi(i,:) <= trn(i) ) .and. ( zi(i,:) >= wbr(i) ) ) - tau(i,0,:) = tauoro(i) * & - ( trn(i) - zi(i,:) ) / & - ( trn(i) - wbr(i) ) - endwhere - src_level(i) = 1 ! disable any more tau calculation - end if - ! Case: Turning AND breaking. Turning BELOW breaking - if ( (wbr(i) > 0. ).and.(trn(i) < wbr(i)).and.(trn_level(i)>=1) ) then - tauoro(i) = 0. - where ( ( zi(i,:) < wbr(i) ) .and. ( zi(i,:) >= tlb(i) ) ) - tau(i,0,:) = tauoro(i) + (taudsw(i)-tauoro(i)) * & - ( wbr(i) - zi(i,:) ) / & - ( wbr(i) - tlb(i) ) - endwhere - src_level(i) = 1 ! disable any more tau calculation - end if - end do - end if - -end subroutine gw_rdg_break_trap - - -!!!!!!!!!!!!!!!!!!!!!! -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -subroutine endrun(msg) - - integer :: iulog - - character(len=*), intent(in), optional :: msg ! string to be printed - - iulog=6 - - if (present (msg)) then - write(iulog,*)'ENDRUN:', msg - else - write(iulog,*)'ENDRUN: called without a message string' - end if - - stop - -end subroutine endrun - - -end module gw_rdg diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/gw_utils.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/gw_utils.F90 deleted file mode 100644 index 20cbb61cd..000000000 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/gw_utils.F90 +++ /dev/null @@ -1,141 +0,0 @@ -module gw_utils - -! -! This module contains utility code for the gravity wave modules. -! - -implicit none -private - -! Real kind for gravity wave parameterization. -integer,public,parameter :: GW_R4 = SELECTED_REAL_KIND(6,37) -integer,public,parameter :: GW_R8 = SELECTED_REAL_KIND(15,307) -integer,public,parameter :: GW_PRC = GW_R4 - -! Public interface - -interface get_unit_vector - module procedure get_unit_vector_r4 - module procedure get_unit_vector_r8 -end interface get_unit_vector -public :: get_unit_vector - -interface dot_2d - module procedure dot_2d_r4 - module procedure dot_2d_r8 -end interface dot_2d -public :: dot_2d - -interface midpoint_interp - module procedure midpoint_interp_r4 - module procedure midpoint_interp_r8 -end interface midpoint_interp -public :: midpoint_interp - -contains - -! Take two components of a vector, and find the unit vector components and -! total magnitude. -subroutine get_unit_vector_r4(u, v, u_n, v_n, mag) - real(GW_R4), intent(in) :: u(:) - real(GW_R4), intent(in) :: v(:) - real(GW_R4), intent(out) :: u_n(:) - real(GW_R4), intent(out) :: v_n(:) - real(GW_R4), intent(out) :: mag(:) - - integer :: i - - mag = sqrt(u*u + v*v) - - ! Has to be a loop/if instead of a where, because floating point - ! exceptions can trigger even on a masked divide-by-zero operation - ! (especially on Intel). - do i = 1, size(mag) - if (mag(i) > 0._GW_R4) then - u_n(i) = u(i)/mag(i) - v_n(i) = v(i)/mag(i) - else - u_n(i) = 0._GW_R4 - v_n(i) = 0._GW_R4 - end if - end do - -end subroutine get_unit_vector_r4 - -subroutine get_unit_vector_r8(u, v, u_n, v_n, mag) - real(GW_R8), intent(in) :: u(:) - real(GW_R8), intent(in) :: v(:) - real(GW_R8), intent(out) :: u_n(:) - real(GW_R8), intent(out) :: v_n(:) - real(GW_R8), intent(out) :: mag(:) - - integer :: i - - mag = sqrt(u*u + v*v) - - ! Has to be a loop/if instead of a where, because floating point - ! exceptions can trigger even on a masked divide-by-zero operation - ! (especially on Intel). - do i = 1, size(mag) - if (mag(i) > 0._GW_R8) then - u_n(i) = u(i)/mag(i) - v_n(i) = v(i)/mag(i) - else - u_n(i) = 0._GW_R8 - v_n(i) = 0._GW_R8 - end if - end do - -end subroutine get_unit_vector_r8 - -! Vectorized version of a 2D dot product (since the intrinsic dot_product -! is more suitable for arrays of contiguous vectors). -function dot_2d_r4(u1, v1, u2, v2) result(dot_2d) - real(GW_R4), intent(in) :: u1(:), v1(:) - real(GW_R4), intent(in) :: u2(:), v2(:) - - real(GW_R4) :: dot_2d(size(u1)) - - dot_2d = u1*u2 + v1*v2 - -end function dot_2d_r4 - -function dot_2d_r8(u1, v1, u2, v2) result(dot_2d) - real(GW_R8), intent(in) :: u1(:), v1(:) - real(GW_R8), intent(in) :: u2(:), v2(:) - - real(GW_R8) :: dot_2d(size(u1)) - - dot_2d = u1*u2 + v1*v2 - -end function dot_2d_r8 - -! Pure function that interpolates the values of the input array along -! dimension 2. This is obviously not a very generic routine, unlike, say, -! CAM's lininterp. But it's used often enough that it seems worth providing -! here. -pure function midpoint_interp_r4(arr) result(interp) - real(GW_R4), intent(in) :: arr(:,:) - real(GW_R4) :: interp(size(arr,1),size(arr,2)-1) - - integer :: i - - do i = 1, size(interp,2) - interp(:,i) = 0.5_GW_R4 * (arr(:,i)+arr(:,i+1)) - end do - -end function midpoint_interp_r4 - -pure function midpoint_interp_r8(arr) result(interp) - real(GW_R8), intent(in) :: arr(:,:) - real(GW_R8) :: interp(size(arr,1),size(arr,2)-1) - - integer :: i - - do i = 1, size(interp,2) - interp(:,i) = 0.5_GW_R8 * (arr(:,i)+arr(:,i+1)) - end do - -end function midpoint_interp_r8 - -end module gw_utils diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/interpolate_data.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/interpolate_data.F90 deleted file mode 100644 index 63bff1c86..000000000 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/interpolate_data.F90 +++ /dev/null @@ -1,1094 +0,0 @@ -module interpolate_data -! Description: -! Routines for interpolation of data in latitude, longitude and time. -! -! Author: Gathered from a number of places and put into the current format by Jim Edwards -! -! Modules Used: -! -!++jtb -!replaced with parameter statement - !!use cam_abortutils, only: endrun - !!use cam_logfile, only: iulog - - implicit none - private -! -! Public Methods: -! - - public :: interp_type, lininterp, vertinterp, bilin, get_timeinterp_factors - public :: lininterp_init, lininterp_finish - - type interp_type - real, pointer :: wgts(:) - real, pointer :: wgtn(:) - integer, pointer :: jjm(:) - integer, pointer :: jjp(:) - end type interp_type - interface lininterp - module procedure lininterp_original - module procedure lininterp_full1d - module procedure lininterp1d - module procedure lininterp2d2d - module procedure lininterp2d1d - module procedure lininterp3d2d - end interface - -integer, parameter, public :: extrap_method_zero = 0 -integer, parameter, public :: extrap_method_bndry = 1 -integer, parameter, public :: extrap_method_cycle = 2 - -integer, parameter :: iulog=6 - -contains - subroutine lininterp_full1d (arrin, yin, nin, arrout, yout, nout) - integer, intent(in) :: nin, nout - real, intent(in) :: arrin(nin), yin(nin), yout(nout) - real, intent(out) :: arrout(nout) - type (interp_type) :: interp_wgts - - call lininterp_init(yin, nin, yout, nout, extrap_method_bndry, interp_wgts) - call lininterp1d(arrin, nin, arrout, nout, interp_wgts) - call lininterp_finish(interp_wgts) - - end subroutine lininterp_full1d - - subroutine lininterp_init(yin, nin, yout, nout, extrap_method, interp_wgts, & - cyclicmin, cyclicmax) -! -! Description: -! Initialize a variable of type(interp_type) with weights for linear interpolation. -! this variable can then be used in calls to lininterp1d and lininterp2d. -! yin is a 1d array of length nin of locations to interpolate from - this array must -! be monotonic but can be increasing or decreasing -! yout is a 1d array of length nout of locations to interpolate to, this array need -! not be ordered -! extrap_method determines how to handle yout points beyond the bounds of yin -! if 0 set values outside output grid to 0 -! if 1 set to boundary value -! if 2 set to cyclic boundaries -! optional values cyclicmin and cyclicmax can be used to set the bounds of the -! cyclic mapping - these default to 0 and 360. -! - - integer, intent(in) :: nin - integer, intent(in) :: nout - real, intent(in) :: yin(:) ! input mesh - real, intent(in) :: yout(:) ! output mesh - integer, intent(in) :: extrap_method ! if 0 set values outside output grid to 0 - ! if 1 set to boundary value - ! if 2 set to cyclic boundaries - real, intent(in), optional :: cyclicmin, cyclicmax - - type (interp_type), intent(out) :: interp_wgts - real :: cmin, cmax - real :: extrap - real :: dyinwrap - real :: ratio - real :: avgdyin - integer :: i, j, icount - integer :: jj - real, pointer :: wgts(:) - real, pointer :: wgtn(:) - integer, pointer :: jjm(:) - integer, pointer :: jjp(:) - logical :: increasing - ! - ! Check validity of input coordinate arrays: must be monotonically increasing, - ! and have a total of at least 2 elements - ! - if (nin.lt.2) then - STOP ! call endrun('LININTERP: Must have at least 2 input points for interpolation') - end if - if(present(cyclicmin)) then - cmin=cyclicmin - else - cmin=0 - end if - if(present(cyclicmax)) then - cmax=cyclicmax - else - cmax=360 - end if - if(cmax<=cmin) then - STOP ! call endrun('LININTERP: cyclic min value must be < max value') - end if - increasing=.true. - icount = 0 - do j=1,nin-1 - if (yin(j).gt.yin(j+1)) icount = icount + 1 - end do - if(icount.eq.nin-1) then - increasing = .false. - icount=0 - endif - if (icount.gt.0) then - STOP ! call endrun('LININTERP: Non-monotonic input coordinate array found') - end if - allocate(interp_wgts%jjm(nout), & - interp_wgts%jjp(nout), & - interp_wgts%wgts(nout), & - interp_wgts%wgtn(nout)) - - jjm => interp_wgts%jjm - jjp => interp_wgts%jjp - wgts => interp_wgts%wgts - wgtn => interp_wgts%wgtn - - ! - ! Initialize index arrays for later checking - ! - jjm = 0 - jjp = 0 - - extrap = 0. - select case (extrap_method) - case (extrap_method_zero) - ! - ! For values which extend beyond boundaries, set weights - ! such that values will be 0. - ! - do j=1,nout - if(increasing) then - if (yout(j).lt.yin(1)) then - jjm(j) = 1 - jjp(j) = 1 - wgts(j) = 0. - wgtn(j) = 0. - extrap = extrap + 1. - else if (yout(j).gt.yin(nin)) then - jjm(j) = nin - jjp(j) = nin - wgts(j) = 0. - wgtn(j) = 0. - extrap = extrap + 1. - end if - else - if (yout(j).gt.yin(1)) then - jjm(j) = 1 - jjp(j) = 1 - wgts(j) = 0. - wgtn(j) = 0. - extrap = extrap + 1. - else if (yout(j).lt.yin(nin)) then - jjm(j) = nin - jjp(j) = nin - wgts(j) = 0. - wgtn(j) = 0. - extrap = extrap + 1. - end if - end if - end do - case (extrap_method_bndry) - ! - ! For values which extend beyond boundaries, set weights - ! such that values will just be copied. - ! - do j=1,nout - if(increasing) then - if (yout(j).le.yin(1)) then - jjm(j) = 1 - jjp(j) = 1 - wgts(j) = 1. - wgtn(j) = 0. - extrap = extrap + 1. - else if (yout(j).gt.yin(nin)) then - jjm(j) = nin - jjp(j) = nin - wgts(j) = 1. - wgtn(j) = 0. - extrap = extrap + 1. - end if - else - if (yout(j).gt.yin(1)) then - jjm(j) = 1 - jjp(j) = 1 - wgts(j) = 1. - wgtn(j) = 0. - extrap = extrap + 1. - else if (yout(j).le.yin(nin)) then - jjm(j) = nin - jjp(j) = nin - wgts(j) = 1. - wgtn(j) = 0. - extrap = extrap + 1. - end if - end if - end do - case (extrap_method_cycle) - ! - ! For values which extend beyond boundaries, set weights - ! for circular boundaries - ! - dyinwrap = yin(1) + (cmax-cmin) - yin(nin) - avgdyin = abs(yin(nin)-yin(1))/(nin-1.) - ratio = dyinwrap/avgdyin - if (ratio < 0.9 .or. ratio > 1.1) then - write(iulog,*) 'Lininterp: Bad dyinwrap value =',dyinwrap,& - ' avg=', avgdyin, yin(1),yin(nin) - STOP ! call endrun('interpolate_data') - end if - - do j=1,nout - if(increasing) then - if (yout(j) <= yin(1)) then - jjm(j) = nin - jjp(j) = 1 - wgts(j) = (yin(1)-yout(j))/dyinwrap - wgtn(j) = (yout(j)+(cmax-cmin) - yin(nin))/dyinwrap - else if (yout(j) > yin(nin)) then - jjm(j) = nin - jjp(j) = 1 - wgts(j) = (yin(1)+(cmax-cmin)-yout(j))/dyinwrap - wgtn(j) = (yout(j)-yin(nin))/dyinwrap - end if - else - if (yout(j) > yin(1)) then - jjm(j) = nin - jjp(j) = 1 - wgts(j) = (yin(1)-yout(j))/dyinwrap - wgtn(j) = (yout(j)+(cmax-cmin) - yin(nin))/dyinwrap - else if (yout(j) <= yin(nin)) then - jjm(j) = nin - jjp(j) = 1 - wgts(j) = (yin(1)+(cmax-cmin)-yout(j))/dyinwrap - wgtn(j) = (yout(j)+(cmax-cmin)-yin(nin))/dyinwrap - end if - - endif - end do - end select - - ! - ! Loop though output indices finding input indices and weights - ! - if(increasing) then - do j=1,nout - do jj=1,nin-1 - if (yout(j).gt.yin(jj) .and. yout(j).le.yin(jj+1)) then - jjm(j) = jj - jjp(j) = jj + 1 - wgts(j) = (yin(jj+1)-yout(j))/(yin(jj+1)-yin(jj)) - wgtn(j) = (yout(j)-yin(jj))/(yin(jj+1)-yin(jj)) - exit - end if - end do - end do - else - do j=1,nout - do jj=1,nin-1 - if (yout(j).le.yin(jj) .and. yout(j).gt.yin(jj+1)) then - jjm(j) = jj - jjp(j) = jj + 1 - wgts(j) = (yin(jj+1)-yout(j))/(yin(jj+1)-yin(jj)) - wgtn(j) = (yout(j)-yin(jj))/(yin(jj+1)-yin(jj)) - exit - end if - end do - end do - end if - - ! - ! Check that interp/extrap points have been found for all outputs - ! - icount = 0 - do j=1,nout - if (jjm(j).eq.0 .or. jjp(j).eq.0) icount = icount + 1 - ratio=wgts(j)+wgtn(j) - if((ratio<0.9.or.ratio>1.1).and.extrap_method.ne.0) then - write(iulog,*) j, wgts(j),wgtn(j),jjm(j),jjp(j), increasing,extrap_method - STOP ! call endrun('Bad weight computed in LININTERP_init') - end if - end do - if (icount.gt.0) then - STOP ! call endrun('LININTERP: Point found without interp indices') - end if - - end subroutine lininterp_init - - subroutine lininterp1d (arrin, n1, arrout, m1, interp_wgts) - !----------------------------------------------------------------------- - ! - ! Purpose: Do a linear interpolation from input mesh to output - ! mesh with weights as set in lininterp_init. - ! - ! - ! Author: Jim Edwards - ! - !----------------------------------------------------------------------- - !----------------------------------------------------------------------- - implicit none - !----------------------------------------------------------------------- - ! - ! Arguments - ! - integer, intent(in) :: n1 ! number of input latitudes - integer, intent(in) :: m1 ! number of output latitudes - - real, intent(in) :: arrin(n1) ! input array of values to interpolate - type(interp_type), intent(in) :: interp_wgts - real, intent(out) :: arrout(m1) ! interpolated array - - ! - ! Local workspace - ! - integer j ! latitude indices - integer, pointer :: jjm(:) - integer, pointer :: jjp(:) - - real, pointer :: wgts(:) - real, pointer :: wgtn(:) - - - jjm => interp_wgts%jjm - jjp => interp_wgts%jjp - wgts => interp_wgts%wgts - wgtn => interp_wgts%wgtn - - ! - ! Do the interpolation - ! - do j=1,m1 - arrout(j) = arrin(jjm(j))*wgts(j) + arrin(jjp(j))*wgtn(j) - end do - - return - end subroutine lininterp1d - - subroutine lininterp2d2d(arrin, n1, n2, arrout, m1, m2, wgt1, wgt2) - implicit none - !----------------------------------------------------------------------- - ! - ! Arguments - ! - integer, intent(in) :: n1, n2, m1, m2 - real, intent(in) :: arrin(n1,n2) ! input array of values to interpolate - type(interp_type), intent(in) :: wgt1, wgt2 - real, intent(out) :: arrout(m1,m2) ! interpolated array - ! - ! locals - ! - integer i,j ! indices - integer, pointer :: iim(:), jjm(:) - integer, pointer :: iip(:), jjp(:) - - real, pointer :: wgts1(:), wgts2(:) - real, pointer :: wgtn1(:), wgtn2(:) - - real :: arrtmp(n1,m2) - - - jjm => wgt2%jjm - jjp => wgt2%jjp - wgts2 => wgt2%wgts - wgtn2 => wgt2%wgtn - - iim => wgt1%jjm - iip => wgt1%jjp - wgts1 => wgt1%wgts - wgtn1 => wgt1%wgtn - - do j=1,m2 - do i=1,n1 - arrtmp(i,j) = arrin(i,jjm(j))*wgts2(j) + arrin(i,jjp(j))*wgtn2(j) - end do - end do - - do j=1,m2 - do i=1,m1 - arrout(i,j) = arrtmp(iim(i),j)*wgts1(i) + arrtmp(iip(i),j)*wgtn1(i) - end do - end do - - end subroutine lininterp2d2d - subroutine lininterp2d1d(arrin, n1, n2, arrout, m1, wgt1, wgt2, fldname) - implicit none - !----------------------------------------------------------------------- - ! - ! Arguments - ! - integer, intent(in) :: n1, n2, m1 - real, intent(in) :: arrin(n1,n2) ! input array of values to interpolate - type(interp_type), intent(in) :: wgt1, wgt2 - real, intent(out) :: arrout(m1) ! interpolated array - character(len=*), intent(in), optional :: fldname(:) - ! - ! locals - ! - integer i ! indices - integer, pointer :: iim(:), jjm(:) - integer, pointer :: iip(:), jjp(:) - - real, pointer :: wgts(:), wgte(:) - real, pointer :: wgtn(:), wgtw(:) - - jjm => wgt2%jjm - jjp => wgt2%jjp - wgts => wgt2%wgts - wgtn => wgt2%wgtn - - iim => wgt1%jjm - iip => wgt1%jjp - wgtw => wgt1%wgts - wgte => wgt1%wgtn - - do i=1,m1 - arrout(i) = arrin(iim(i),jjm(i))*wgtw(i)*wgts(i)+arrin(iip(i),jjm(i))*wgte(i)*wgts(i) + & - arrin(iim(i),jjp(i))*wgtw(i)*wgtn(i)+arrin(iip(i),jjp(i))*wgte(i)*wgtn(i) - end do - - - end subroutine lininterp2d1d - subroutine lininterp3d2d(arrin, n1, n2, n3, arrout, m1, len1, wgt1, wgt2) - implicit none - !----------------------------------------------------------------------- - ! - ! Arguments - ! - integer, intent(in) :: n1, n2, n3, m1, len1 ! m1 is to len1 as ncols is to pcols - real, intent(in) :: arrin(n1,n2,n3) ! input array of values to interpolate - type(interp_type), intent(in) :: wgt1, wgt2 - real, intent(out) :: arrout(len1, n3) ! interpolated array - - ! - ! locals - ! - integer i, k ! indices - integer, pointer :: iim(:), jjm(:) - integer, pointer :: iip(:), jjp(:) - - real, pointer :: wgts(:), wgte(:) - real, pointer :: wgtn(:), wgtw(:) - - jjm => wgt2%jjm - jjp => wgt2%jjp - wgts => wgt2%wgts - wgtn => wgt2%wgtn - - iim => wgt1%jjm - iip => wgt1%jjp - wgtw => wgt1%wgts - wgte => wgt1%wgtn - - do k=1,n3 - do i=1,m1 - arrout(i,k) = arrin(iim(i),jjm(i),k)*wgtw(i)*wgts(i)+arrin(iip(i),jjm(i),k)*wgte(i)*wgts(i) + & - arrin(iim(i),jjp(i),k)*wgtw(i)*wgtn(i)+arrin(iip(i),jjp(i),k)*wgte(i)*wgtn(i) - end do - end do - - end subroutine lininterp3d2d - - - - - subroutine lininterp_finish(interp_wgts) - type(interp_type) :: interp_wgts - - deallocate(interp_wgts%jjm, & - interp_wgts%jjp, & - interp_wgts%wgts, & - interp_wgts%wgtn) - - nullify(interp_wgts%jjm, & - interp_wgts%jjp, & - interp_wgts%wgts, & - interp_wgts%wgtn) - end subroutine lininterp_finish - - subroutine lininterp_original (arrin, yin, nlev, nlatin, arrout, & - yout, nlatout) - !----------------------------------------------------------------------- - ! - ! Purpose: Do a linear interpolation from input mesh defined by yin to output - ! mesh defined by yout. Where extrapolation is necessary, values will - ! be copied from the extreme edge of the input grid. Vectorization is over - ! the vertical (nlev) dimension. - ! - ! Method: Check validity of input, then determine weights, then do the N-S interpolation. - ! - ! Author: Jim Rosinski - ! Modified: Jim Edwards so that there is no requirement of monotonacity on the yout array - ! - !----------------------------------------------------------------------- - implicit none - !----------------------------------------------------------------------- - ! - ! Arguments - ! - integer, intent(in) :: nlev ! number of vertical levels - integer, intent(in) :: nlatin ! number of input latitudes - integer, intent(in) :: nlatout ! number of output latitudes - - real, intent(in) :: arrin(nlev,nlatin) ! input array of values to interpolate - real, intent(in) :: yin(nlatin) ! input mesh - real, intent(in) :: yout(nlatout) ! output mesh - - real, intent(out) :: arrout(nlev,nlatout) ! interpolated array - ! - ! Local workspace - ! - integer j, jj ! latitude indices - integer jjprev ! latitude indices - integer k ! level index - integer icount ! number of values - - real extrap ! percent grid non-overlap - ! - ! Dynamic - ! - integer :: jjm(nlatout) - integer :: jjp(nlatout) - - real :: wgts(nlatout) - real :: wgtn(nlatout) - ! - ! Check validity of input coordinate arrays: must be monotonically increasing, - ! and have a total of at least 2 elements - ! - if (nlatin.lt.2) then - STOP ! call endrun('LININTERP: Must have at least 2 input points for interpolation') - end if - - icount = 0 - do j=1,nlatin-1 - if (yin(j).gt.yin(j+1)) icount = icount + 1 - end do - - - if (icount.gt.0) then - STOP ! call endrun('LININTERP: Non-monotonic coordinate array(s) found') - end if - ! - ! Initialize index arrays for later checking - ! - do j=1,nlatout - jjm(j) = 0 - jjp(j) = 0 - end do - ! - ! For values which extend beyond N and S boundaries, set weights - ! such that values will just be copied. - ! - extrap = 0. - - do j=1,nlatout - if (yout(j).le.yin(1)) then - jjm(j) = 1 - jjp(j) = 1 - wgts(j) = 1. - wgtn(j) = 0. - extrap=extrap+1. - else if (yout(j).gt.yin(nlatin)) then - jjm(j) = nlatin - jjp(j) = nlatin - wgts(j) = 1. - wgtn(j) = 0. - extrap=extrap+1. - endif - end do - - ! - ! Loop though output indices finding input indices and weights - ! - do j=1,nlatout - do jj=1,nlatin-1 - if (yout(j).gt.yin(jj) .and. yout(j).le.yin(jj+1)) then - jjm(j) = jj - jjp(j) = jj + 1 - wgts(j) = (yin(jj+1)-yout(j))/(yin(jj+1)-yin(jj)) - wgtn(j) = (yout(j)-yin(jj))/(yin(jj+1)-yin(jj)) - exit - end if - end do - end do - ! - ! Check that interp/extrap points have been found for all outputs - ! - icount = 0 - do j=1,nlatout - if (jjm(j).eq.0 .or. jjp(j).eq.0) then - icount = icount + 1 - end if - end do - if (icount.gt.0) then - STOP ! call endrun('LININTERP: Point found without interp indices') - end if - ! - ! Do the interpolation - ! - do j=1,nlatout - do k=1,nlev - arrout(k,j) = arrin(k,jjm(j))*wgts(j) + arrin(k,jjp(j))*wgtn(j) - end do - end do - - return - end subroutine lininterp_original - - - subroutine bilin (arrin, xin, yin, nlondin, nlonin, & - nlevdin, nlev, nlatin, arrout, xout, & - yout, nlondout, nlonout, nlevdout, nlatout) - !----------------------------------------------------------------------- - ! - ! Purpose: - ! - ! Do a bilinear interpolation from input mesh defined by xin, yin to output - ! mesh defined by xout, yout. Circularity is assumed in the x-direction so - ! input x-grid must be in degrees east and must start from Greenwich. When - ! extrapolation is necessary in the N-S direction, values will be copied - ! from the extreme edge of the input grid. Vectorization is over the - ! longitude dimension. Input grid is assumed rectangular. Output grid - ! is assumed ragged, i.e. xout is a 2-d mesh. - ! - ! Method: Interpolate first in longitude, then in latitude. - ! - ! Author: Jim Rosinski - ! - !----------------------------------------------------------------------- - !++jtb - !!use cam_abortutils, only: endrun - !----------------------------------------------------------------------- - implicit none - !----------------------------------------------------------------------- - ! - ! Input arguments - ! - integer, intent(in) :: nlondin ! longitude dimension of input grid - integer, intent(in) :: nlonin ! number of real longitudes (input) - integer, intent(in) :: nlevdin ! vertical dimension of input grid - integer, intent(in) :: nlev ! number of vertical levels - integer, intent(in) :: nlatin ! number of input latitudes - integer, intent(in) :: nlatout ! number of output latitudes - integer, intent(in) :: nlondout ! longitude dimension of output grid - integer, intent(in) :: nlonout(nlatout) ! number of output longitudes per lat - integer, intent(in) :: nlevdout ! vertical dimension of output grid - - real, intent(in) :: arrin(nlondin,nlevdin,nlatin) ! input array of values to interpolate - real, intent(in) :: xin(nlondin) ! input x mesh - real, intent(in) :: yin(nlatin) ! input y mesh - real, intent(in) :: xout(nlondout,nlatout) ! output x mesh - real, intent(in) :: yout(nlatout) ! output y mesh - ! - ! Output arguments - ! - real, intent(out) :: arrout(nlondout,nlevdout,nlatout) ! interpolated array - ! - ! Local workspace - ! - integer :: i, ii, iw, ie, iiprev ! longitude indices - integer :: j, jj, js, jn, jjprev ! latitude indices - integer :: k ! level index - integer :: icount ! number of bad values - - real :: extrap ! percent grid non-overlap - real :: dxinwrap ! delta-x on input grid for 2-pi - real :: avgdxin ! avg input delta-x - real :: ratio ! compare dxinwrap to avgdxin - real :: sum ! sum of weights (used for testing) - ! - ! Dynamic - ! - integer :: iim(nlondout) ! interpolation index to the left - integer :: iip(nlondout) ! interpolation index to the right - integer :: jjm(nlatout) ! interpolation index to the south - integer :: jjp(nlatout) ! interpolation index to the north - - real :: wgts(nlatout) ! interpolation weight to the north - real :: wgtn(nlatout) ! interpolation weight to the north - real :: wgte(nlondout) ! interpolation weight to the north - real :: wgtw(nlondout) ! interpolation weight to the north - real :: igrid(nlonin) ! interpolation weight to the north - ! - ! Check validity of input coordinate arrays: must be monotonically increasing, - ! and have a total of at least 2 elements - ! - if (nlonin < 2 .or. nlatin < 2) then - STOP ! call endrun ('BILIN: Must have at least 2 input points for interpolation') - end if - - if (xin(1) < 0. .or. xin(nlonin) > 360.) then - STOP ! call endrun ('BILIN: Input x-grid must be between 0 and 360') - end if - - icount = 0 - do i=1,nlonin-1 - if (xin(i) >= xin(i+1)) icount = icount + 1 - end do - - do j=1,nlatin-1 - if (yin(j) >= yin(j+1)) icount = icount + 1 - end do - - do j=1,nlatout-1 - if (yout(j) >= yout(j+1)) icount = icount + 1 - end do - - do j=1,nlatout - do i=1,nlonout(j)-1 - if (xout(i,j) >= xout(i+1,j)) icount = icount + 1 - end do - end do - - if (icount > 0) then - STOP ! call endrun ('BILIN: Non-monotonic coordinate array(s) found') - end if - - if (yout(nlatout) <= yin(1) .or. yout(1) >= yin(nlatin)) then - STOP ! call endrun ('BILIN: No overlap in y-coordinates') - end if - - do j=1,nlatout - if (xout(1,j) < 0. .or. xout(nlonout(j),j) > 360.) then - STOP ! call endrun ('BILIN: Output x-grid must be between 0 and 360') - end if - - if (xout(nlonout(j),j) <= xin(1) .or. & - xout(1,j) >= xin(nlonin)) then - STOP ! call endrun ('BILIN: No overlap in x-coordinates') - end if - end do - ! - ! Initialize index arrays for later checking - ! - do j=1,nlatout - jjm(j) = 0 - jjp(j) = 0 - end do - ! - ! For values which extend beyond N and S boundaries, set weights - ! such that values will just be copied. - ! - do js=1,nlatout - if (yout(js) > yin(1)) exit - jjm(js) = 1 - jjp(js) = 1 - wgts(js) = 1. - wgtn(js) = 0. - end do - - do jn=nlatout,1,-1 - if (yout(jn) <= yin(nlatin)) exit - jjm(jn) = nlatin - jjp(jn) = nlatin - wgts(jn) = 1. - wgtn(jn) = 0. - end do - ! - ! Loop though output indices finding input indices and weights - ! - jjprev = 1 - do j=js,jn - do jj=jjprev,nlatin-1 - if (yout(j) > yin(jj) .and. yout(j) <= yin(jj+1)) then - jjm(j) = jj - jjp(j) = jj + 1 - wgts(j) = (yin(jj+1) - yout(j)) / (yin(jj+1) - yin(jj)) - wgtn(j) = (yout(j) - yin(jj)) / (yin(jj+1) - yin(jj)) - goto 30 - end if - end do - STOP ! call endrun ('BILIN: Failed to find interp values') -30 jjprev = jj - end do - - dxinwrap = xin(1) + 360. - xin(nlonin) - ! - ! Check for sane dxinwrap values. Allow to differ no more than 10% from avg - ! - avgdxin = (xin(nlonin)-xin(1))/(nlonin-1.) - ratio = dxinwrap/avgdxin - if (ratio < 0.9 .or. ratio > 1.1) then - write(iulog,*)'BILIN: Insane dxinwrap value =',dxinwrap,' avg=', avgdxin - STOP ! call endrun - end if - ! - ! Check that interp/extrap points have been found for all outputs, and that - ! they are reasonable (i.e. within 32-bit roundoff) - ! - icount = 0 - do j=1,nlatout - if (jjm(j) == 0 .or. jjp(j) == 0) icount = icount + 1 - sum = wgts(j) + wgtn(j) - if (sum < 0.99999 .or. sum > 1.00001) icount = icount + 1 - if (wgts(j) < 0. .or. wgts(j) > 1.) icount = icount + 1 - if (wgtn(j) < 0. .or. wgtn(j) > 1.) icount = icount + 1 - end do - - if (icount > 0) then - STOP ! call endrun ('BILIN: something bad in latitude indices or weights') - end if - ! - ! Do the bilinear interpolation - ! - do j=1,nlatout - ! - ! Initialize index arrays for later checking - ! - do i=1,nlondout - iim(i) = 0 - iip(i) = 0 - end do - ! - ! For values which extend beyond E and W boundaries, set weights such that - ! values will be interpolated between E and W edges of input grid. - ! - do iw=1,nlonout(j) - if (xout(iw,j) > xin(1)) exit - iim(iw) = nlonin - iip(iw) = 1 - wgtw(iw) = (xin(1) - xout(iw,j)) /dxinwrap - wgte(iw) = (xout(iw,j)+360. - xin(nlonin))/dxinwrap - end do - - do ie=nlonout(j),1,-1 - if (xout(ie,j) <= xin(nlonin)) exit - iim(ie) = nlonin - iip(ie) = 1 - wgtw(ie) = (xin(1)+360. - xout(ie,j)) /dxinwrap - wgte(ie) = (xout(ie,j) - xin(nlonin))/dxinwrap - end do - ! - ! Loop though output indices finding input indices and weights - ! - iiprev = 1 - do i=iw,ie - do ii=iiprev,nlonin-1 - if (xout(i,j) > xin(ii) .and. xout(i,j) <= xin(ii+1)) then - iim(i) = ii - iip(i) = ii + 1 - wgtw(i) = (xin(ii+1) - xout(i,j)) / (xin(ii+1) - xin(ii)) - wgte(i) = (xout(i,j) - xin(ii)) / (xin(ii+1) - xin(ii)) - goto 60 - end if - end do - STOP ! call endrun ('BILIN: Failed to find interp values') -60 iiprev = ii - end do - - icount = 0 - do i=1,nlonout(j) - if (iim(i) == 0 .or. iip(i) == 0) icount = icount + 1 - sum = wgtw(i) + wgte(i) - if (sum < 0.99999 .or. sum > 1.00001) icount = icount + 1 - if (wgtw(i) < 0. .or. wgtw(i) > 1.) icount = icount + 1 - if (wgte(i) < 0. .or. wgte(i) > 1.) icount = icount + 1 - end do - - if (icount > 0) then - write(iulog,*)'BILIN: j=',j,' Something bad in longitude indices or weights' - STOP ! call endrun - end if - ! - ! Do the interpolation, 1st in longitude then latitude - ! - do k=1,nlev - do i=1,nlonin - igrid(i) = arrin(i,k,jjm(j))*wgts(j) + arrin(i,k,jjp(j))*wgtn(j) - end do - - do i=1,nlonout(j) - arrout(i,k,j) = igrid(iim(i))*wgtw(i) + igrid(iip(i))*wgte(i) - end do - end do - end do - - - return - end subroutine bilin - - subroutine vertinterp(ncol, ncold, nlev, pmid, pout, arrin, arrout) - - !----------------------------------------------------------------------- - ! - ! Purpose: - ! Vertically interpolate input array to output pressure level - ! Copy values at boundaries. - ! - ! Method: - ! - ! Author: - ! - !----------------------------------------------------------------------- - - implicit none - - !------------------------------Arguments-------------------------------- - integer , intent(in) :: ncol ! column dimension - integer , intent(in) :: ncold ! declared column dimension - integer , intent(in) :: nlev ! vertical dimension - real, intent(in) :: pmid(ncold,nlev) ! input level pressure levels - real, intent(in) :: pout ! output pressure level - real, intent(in) :: arrin(ncold,nlev) ! input array - real, intent(out) :: arrout(ncold) ! output array (interpolated) - !-------------------------------------------------------------------------- - - !---------------------------Local variables----------------------------- - integer i,k ! indices - integer kupper(ncold) ! Level indices for interpolation - real dpu ! upper level pressure difference - real dpl ! lower level pressure difference - logical found(ncold) ! true if input levels found - logical error ! error flag - !----------------------------------------------------------------- - ! - ! Initialize index array and logical flags - ! - do i=1,ncol - found(i) = .false. - kupper(i) = 1 - end do - error = .false. - ! - ! Store level indices for interpolation. - ! If all indices for this level have been found, - ! do the interpolation - ! - do k=1,nlev-1 - do i=1,ncol - if ((.not. found(i)) .and. pmid(i,k)= pmid(i,nlev)) then - arrout(i) = arrin(i,nlev) - else if (found(i)) then - dpu = pout - pmid(i,kupper(i)) - dpl = pmid(i,kupper(i)+1) - pout - arrout(i) = (arrin(i,kupper(i) )*dpl + arrin(i,kupper(i)+1)*dpu)/(dpl + dpu) - else - error = .true. - end if - end do - ! - ! Error check - ! - if (error) then - STOP ! call endrun ('VERTINTERP: ERROR FLAG') - end if - - return - end subroutine vertinterp - - subroutine get_timeinterp_factors (cycflag, np1, cdayminus, cdayplus, cday, & - fact1, fact2, str) - !--------------------------------------------------------------------------- - ! - ! Purpose: Determine time interpolation factors (normally for a boundary dataset) - ! for linear interpolation. - ! - ! Method: Assume 365 days per year. Output variable fact1 will be the weight to - ! apply to data at calendar time "cdayminus", and fact2 the weight to apply - ! to data at time "cdayplus". Combining these values will produce a result - ! valid at time "cday". Output arguments fact1 and fact2 will be between - ! 0 and 1, and fact1 + fact2 = 1 to roundoff. - ! - ! Author: Jim Rosinski - ! - !--------------------------------------------------------------------------- - implicit none - ! - ! Arguments - ! - logical, intent(in) :: cycflag ! flag indicates whether dataset is being cycled yearly - - integer, intent(in) :: np1 ! index points to forward time slice matching cdayplus - - real, intent(in) :: cdayminus ! calendar day of rearward time slice - real, intent(in) :: cdayplus ! calendar day of forward time slice - real, intent(in) :: cday ! calenar day to be interpolated to - real, intent(out) :: fact1 ! time interpolation factor to apply to rearward time slice - real, intent(out) :: fact2 ! time interpolation factor to apply to forward time slice - - character(len=*), intent(in) :: str ! string to be added to print in case of error (normally the callers name) - ! - ! Local workspace - ! - real :: deltat ! time difference (days) between cdayminus and cdayplus - real, parameter :: daysperyear = 365. ! number of days in a year - ! - ! Initial sanity checks - ! - if (np1 == 1 .and. .not. cycflag) then - STOP ! call endrun ('GETFACTORS:'//str//' cycflag false and forward month index = Jan. not allowed') - end if - - if (np1 < 1) then - STOP ! call endrun ('GETFACTORS:'//str//' input arg np1 must be > 0') - end if - - if (cycflag) then - if ((cday < 1.) .or. (cday > (daysperyear+1.))) then - write(iulog,*) 'GETFACTORS:', str, ' bad cday=',cday - STOP ! call endrun () - end if - else - if (cday < 1.) then - write(iulog,*) 'GETFACTORS:', str, ' bad cday=',cday - STOP ! call endrun () - end if - end if - ! - ! Determine time interpolation factors. Account for December-January - ! interpolation if dataset is being cycled yearly. - ! - if (cycflag .and. np1 == 1) then ! Dec-Jan interpolation - deltat = cdayplus + daysperyear - cdayminus - if (cday > cdayplus) then ! We are in December - fact1 = (cdayplus + daysperyear - cday)/deltat - fact2 = (cday - cdayminus)/deltat - else ! We are in January - fact1 = (cdayplus - cday)/deltat - fact2 = (cday + daysperyear - cdayminus)/deltat - end if - else - deltat = cdayplus - cdayminus - fact1 = (cdayplus - cday)/deltat - fact2 = (cday - cdayminus)/deltat - end if - - if (.not. valid_timeinterp_factors (fact1, fact2)) then - write(iulog,*) 'GETFACTORS: ', str, ' bad fact1 and/or fact2=', fact1, fact2 - STOP ! call endrun () - end if - - return - end subroutine get_timeinterp_factors - - logical function valid_timeinterp_factors (fact1, fact2) - !--------------------------------------------------------------------------- - ! - ! Purpose: check sanity of time interpolation factors to within 32-bit roundoff - ! - !--------------------------------------------------------------------------- - implicit none - - real, intent(in) :: fact1, fact2 ! time interpolation factors - - valid_timeinterp_factors = .true. - - ! The fact1 .ne. fact1 and fact2 .ne. fact2 comparisons are to detect NaNs. - if (abs(fact1+fact2-1.) > 1.e-6 .or. & - fact1 > 1.000001 .or. fact1 < -1.e-6 .or. & - fact2 > 1.000001 .or. fact2 < -1.e-6 .or. & - fact1 .ne. fact1 .or. fact2 .ne. fact2) then - - valid_timeinterp_factors = .false. - end if - - return - end function valid_timeinterp_factors - -end module interpolate_data diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/linear_1d_operators.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/linear_1d_operators.F90 deleted file mode 100644 index 3ea2207c7..000000000 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/linear_1d_operators.F90 +++ /dev/null @@ -1,1187 +0,0 @@ -module linear_1d_operators -#define USE_CONTIGUOUS contiguous, -! This module provides the type "TriDiagOp" to represent operators on a 1D -! grid as tridiagonal matrices, and related types to represent boundary -! conditions. -! -! The focus is on solving diffusion equations with a finite volume method -! in one dimension, but other utility operators are provided, e.g. a second -! order approximation to the first derivative. -! -! In order to allow vectorization to occur, as well as to avoid unnecessary -! copying/reshaping of data in CAM, TriDiagOp actually represents a -! collection of independent operators that can be applied to collections of -! independent data; the innermost index is over independent systems (e.g. -! CAM columns). -! -! A simple example: -! ! First derivative operator -! op = first_derivative(coords) -! ! Convert data to its derivative (extrapolate at boundaries). -! call op%apply(data) -! -! With explicit boundary conditions: -! op = first_derivative(coords, & -! l_bndry=BoundaryFixedFlux(), & -! r_bndry=BoundaryFixedLayer(layer_distance)) -! call op%apply(data, & -! l_cond=BoundaryFlux(flux, dt, thickness), & -! r_cond=BoundaryData(boundary)) -! -! Implicit solution example: -! ! Construct diffusion matrix. -! op = diffusion_operator(coords, d) -! call op%lmult_as_diag(-dt) -! call op%add_to_diag(1.) -! ! Decompose in order to invert the operation. -! decomp = TriDiagDecomp(op) -! ! Diffuse data for one time step (fixed flux boundaries). -! call decomp%left_div(data) - -!++jtb -! replaced this set with simple call abort() -!use shr_log_mod, only: errMsg => shr_log_errMsg -!use shr_sys_mod, only: shr_sys_abort -!--jtb - -use coords_1d, only: Coords1D - -implicit none -private - -! Main type. -public :: TriDiagOp -public :: operator(+) -public :: operator(-) - -! Decomposition used for inversion (left division). -public :: TriDiagDecomp - -! Multiplies by 0. -public :: zero_operator - -! Construct identity. -public :: identity_operator - -! Produce a TriDiagOp that is simply a diagonal matrix. -public :: diagonal_operator - -! For solving the diffusion-advection equation with implicit Euler. -public :: diffusion_operator -public :: advection_operator - -! Derivatives accurate to second order on a non-uniform grid. -public :: first_derivative -public :: second_derivative - -! Boundary condition types. -public :: BoundaryType -public :: BoundaryZero -public :: BoundaryFirstOrder -public :: BoundaryExtrapolate -public :: BoundaryFixedLayer -public :: BoundaryFixedFlux - -! Boundary data types. -public :: BoundaryCond -public :: BoundaryNoData -public :: BoundaryData -public :: BoundaryFlux - -!++jtb - -! TriDiagOp represents operators that can work between nearest neighbors, -! with some extra logic at the boundaries. The implementation is a -! tridiagonal matrix plus boundary info. -type :: TriDiagOp - private - ! The number of independent systems. - integer, public :: nsys - ! The size of the matrix (number of grid cells). - integer, public :: ncel - ! Super-, sub-, and regular diagonals. - real, allocatable :: spr(:,:) - real, allocatable :: sub(:,:) - real, allocatable :: diag(:,:) - ! Buffers to hold boundary data; Details depend on the type of boundary - ! being used. - real, allocatable :: left_bound(:) - real, allocatable :: right_bound(:) - contains - ! Applies the operator to a set of data. - procedure :: apply => apply_tridiag - ! Given the off-diagonal elements, fills in the diagonal so that the - ! operator will have the constant function as an eigenvector with - ! eigenvalue 0. This is used internally as a utility for construction of - ! derivative operators. - procedure :: deriv_diag => make_tridiag_deriv_diag - ! Add/substract another tridiagonal from this one in-place (without - ! creating a temporary object). - procedure :: add => add_in_place_tridiag_ops - procedure :: subtract => subtract_in_place_tridiag_ops - ! Add input vector or scalar to the diagonal. - procedure :: scalar_add_tridiag - procedure :: diagonal_add_tridiag - generic :: add_to_diag => scalar_add_tridiag, diagonal_add_tridiag - ! Treat input vector (or scalar) as if it was the diagonal of an - ! operator, and multiply this operator on the left by that value. - procedure :: scalar_lmult_tridiag - procedure :: diagonal_lmult_tridiag - generic :: lmult_as_diag => & - scalar_lmult_tridiag, diagonal_lmult_tridiag - ! Deallocate and reset. - procedure :: finalize => tridiag_finalize -end type TriDiagOp - -interface operator(+) - module procedure add_tridiag_ops -end interface operator(+) - -interface operator(-) - module procedure subtract_tridiag_ops -end interface operator(-) - -interface TriDiagOp - module procedure new_TriDiagOp -end interface TriDiagOp - -! -! Boundary condition types for the operators. -! -! Note that BoundaryFixedLayer and BoundaryFixedFlux are the only options -! supported for backwards operation (i.e. decomp%left_div). The others are -! meant for direct application only (e.g. to find a derivative). -! -! BoundaryZero means that the operator fixes boundaries to 0. -! BoundaryFirstOrder means a one-sided approximation for the first -! derivative. -! BoundaryExtrapolate means that a second order approximation will be used, -! even at the boundaries. Boundary points do this by using their next- -! nearest neighbor to extrapolate. -! BoundaryFixedLayer means that there's an extra layer outside of the given -! grid, which must be specified when applying/inverting the operator. -! BoundaryFixedFlux is intended to provide a fixed-flux condition for -! typical advection/diffusion operators. It tweaks the edge condition -! to work on an input current rather than a value. -! -! The different types were originally implemented through polymorphism, but -! PGI required this to be done via enum instead. -integer, parameter :: zero_bndry = 0 -integer, parameter :: first_order_bndry = 1 -integer, parameter :: extrapolate_bndry = 2 -integer, parameter :: fixed_layer_bndry = 3 -integer, parameter :: fixed_flux_bndry = 4 - -type :: BoundaryType - private - integer :: bndry_type = fixed_flux_bndry - real, allocatable :: edge_width(:) - contains - procedure :: make_left - procedure :: make_right - procedure :: finalize => boundary_type_finalize -end type BoundaryType - -abstract interface - subroutine deriv_seed(del_minus, del_plus, sub, spr) - real, USE_CONTIGUOUS intent(in) :: del_minus(:) - real, USE_CONTIGUOUS intent(in) :: del_plus(:) - real, USE_CONTIGUOUS intent(out) :: sub(:) - real, USE_CONTIGUOUS intent(out) :: spr(:) - end subroutine deriv_seed -end interface - -interface BoundaryZero - module procedure new_BoundaryZero -end interface BoundaryZero - -interface BoundaryFirstOrder - module procedure new_BoundaryFirstOrder -end interface BoundaryFirstOrder - -interface BoundaryExtrapolate - module procedure new_BoundaryExtrapolate -end interface BoundaryExtrapolate - -interface BoundaryFixedLayer - module procedure new_BoundaryFixedLayer -end interface BoundaryFixedLayer - -interface BoundaryFixedFlux - module procedure new_BoundaryFixedFlux -end interface BoundaryFixedFlux - -! -! Data for boundary conditions themselves. -! -! "No data" conditions perform extrapolation, if BoundaryExtrapolate was -! the boundary type used to construct the operator. -! -! "Data" conditions contain extra data, which effectively extends the -! system with an extra cell. -! -! "Flux" conditions contain prescribed fluxes. -! -! The condition you can use depends on the boundary type from above that -! was used in the operator's construction. For BoundaryFixedLayer use -! BoundaryData. For BoundaryFixedFlux use BoundaryFlux. For everything -! else, use BoundaryNoData. - -! The switches using this enumeration used to be unnecessary due to use of -! polymorphism, but this had to be backed off due to insufficient PGI -! support for type extension. -integer, parameter :: no_data_cond = 0 -integer, parameter :: data_cond = 1 -integer, parameter :: flux_cond = 2 - -type :: BoundaryCond - private - integer :: cond_type = no_data_cond - real, allocatable :: edge_data(:) - contains - procedure :: apply_left - procedure :: apply_right - procedure :: finalize => boundary_cond_finalize -end type BoundaryCond - -! Constructors for different types of BoundaryCond. -interface BoundaryNoData - module procedure new_BoundaryNoData -end interface BoundaryNoData - -interface BoundaryData - module procedure new_BoundaryData -end interface BoundaryData - -interface BoundaryFlux - module procedure new_BoundaryFlux -end interface BoundaryFlux - -! Opaque type to hold a tridiagonal matrix decomposition. -! -! Method used is similar to Richtmyer and Morton (1967,pp 198-201), but -! the order of iteration is reversed, leading to A and C being swapped, and -! some differences in the indexing. -type :: TriDiagDecomp - private - integer :: nsys = 0 - integer :: ncel = 0 - ! These correspond to A_k, E_k, and 1 / (B_k - A_k * E_{k+1}) - real, allocatable :: ca(:,:) - real, allocatable :: ze(:,:) - real, allocatable :: dnom(:,:) -contains - procedure :: left_div => decomp_left_div - procedure :: finalize => decomp_finalize -end type TriDiagDecomp - -interface TriDiagDecomp - module procedure new_TriDiagDecomp -end interface TriDiagDecomp - -contains - -! Operator that sets to 0. -function zero_operator(nsys, ncel) result(op) - ! Sizes for operator. - integer, intent(in) :: nsys, ncel - - type(TriDiagOp) :: op - - op = TriDiagOp(nsys, ncel) - - op%spr = 0. - op%sub = 0. - op%diag = 0. - op%left_bound = 0. - op%right_bound = 0. - -end function zero_operator - -! Operator that does nothing. -function identity_operator(nsys, ncel) result(op) - ! Sizes for operator. - integer, intent(in) :: nsys, ncel - - type(TriDiagOp) :: op - - op = TriDiagOp(nsys, ncel) - - op%spr = 0. - op%sub = 0. - op%diag = 1. - op%left_bound = 0. - op%right_bound = 0. - -end function identity_operator - -! Create an operator that just does an element-wise product by some data. -function diagonal_operator(diag) result(op) - ! Data to multiply by. - real, USE_CONTIGUOUS intent(in) :: diag(:,:) - - type(TriDiagOp) :: op - - op = TriDiagOp(size(diag, 1), size(diag, 2)) - - op%spr = 0. - op%sub = 0. - op%diag = diag - op%left_bound = 0. - op%right_bound = 0. - -end function diagonal_operator - -! Diffusion matrix operator constructor. Given grid coordinates, a set of -! diffusion coefficients, and boundaries, creates a matrix corresponding -! to a finite volume representation of the operation: -! -! d/dx (d_coef * d/dx) -! -! This differs from what you would get from combining the first and second -! derivative operations, which would be more appropriate for a finite -! difference scheme that does not use grid cell averages. -function diffusion_operator(coords, d_coef, l_bndry, r_bndry) & - result(op) - ! Grid cell locations. - type(Coords1D), intent(in) :: coords - ! Diffusion coefficient defined on interfaces. - real, USE_CONTIGUOUS intent(in) :: d_coef(:,:) - ! Objects representing the kind of boundary on each side. - class(BoundaryType), target, intent(in), optional :: l_bndry, r_bndry - ! Output operator. - type(TriDiagOp) :: op - - ! Selectors to implement default boundary. - class(BoundaryType), pointer :: l_bndry_loc, r_bndry_loc - ! Fixed flux is default, no allocation/deallocation needed. - type(BoundaryType), target :: bndry_default - - ! Level index. - integer :: k - - if (present(l_bndry)) then - l_bndry_loc => l_bndry - else - l_bndry_loc => bndry_default - end if - - if (present(r_bndry)) then - r_bndry_loc => r_bndry - else - r_bndry_loc => bndry_default - end if - - ! Allocate the operator. - op = TriDiagOp(coords%n, coords%d) - - ! d_coef over the distance to the next cell gives you the matrix term for - ! flux of material between cells. Dividing by cell thickness translates - ! this to a tendency on the concentration. Hence the basic pattern is - ! d_coef*rdst*rdel. - ! - ! Boundary conditions for a fixed layer simply extend this by calculating - ! the distance to the midpoint of the extra edge layer. - - select case (l_bndry_loc%bndry_type) - case (fixed_layer_bndry) - op%left_bound = 2.*d_coef(:,1)*coords%rdel(:,1) / & - (l_bndry_loc%edge_width+coords%del(:,1)) - case default - op%left_bound = 0. - end select - - do k = 1, coords%d-1 - op%spr(:,k) = d_coef(:,k+1)*coords%rdst(:,k)*coords%rdel(:,k) - op%sub(:,k) = d_coef(:,k+1)*coords%rdst(:,k)*coords%rdel(:,k+1) - end do - - select case (r_bndry_loc%bndry_type) - case (fixed_layer_bndry) - op%right_bound = 2.*d_coef(:,coords%d+1)*coords%rdel(:,coords%d) / & - (r_bndry_loc%edge_width+coords%del(:,coords%d)) - case default - op%right_bound = 0. - end select - - ! Above, we found all off-diagonals. Now get the diagonal. - call op%deriv_diag() - -end function diffusion_operator - -! Advection matrix operator constructor. Similar to diffusion_operator, it -! constructs an operator A corresponding to: -! -! A y = d/dx (-v_coef * y) -! -! Again, this is targeted at representing this operator acting on grid-cell -! averages in a finite volume scheme, rather than a literal representation. -function advection_operator(coords, v_coef, l_bndry, r_bndry) & - result(op) - ! Grid cell locations. - type(Coords1D), intent(in) :: coords - ! Advection coefficient (effective velocity). - real, USE_CONTIGUOUS intent(in) :: v_coef(:,:) - ! Objects representing the kind of boundary on each side. - class(BoundaryType), target, intent(in), optional :: l_bndry, r_bndry - ! Output operator. - type(TriDiagOp) :: op - - ! Selectors to implement default boundary. - class(BoundaryType), pointer :: l_bndry_loc, r_bndry_loc - ! Fixed flux is default, no allocation/deallocation needed. - type(BoundaryType), target :: bndry_default - - ! Negative derivative of v. - real :: v_deriv(coords%n,coords%d) - - if (present(l_bndry)) then - l_bndry_loc => l_bndry - else - l_bndry_loc => bndry_default - end if - - if (present(r_bndry)) then - r_bndry_loc => r_bndry - else - r_bndry_loc => bndry_default - end if - - ! Allocate the operator. - op = TriDiagOp(coords%n, coords%d) - - ! Construct the operator in two stages using the product rule. First - ! create (-v * d/dx), then -dv/dx, and add the two. - ! - ! For the first part, we want to interpolate to interfaces (weighted - ! average involving del/2*dst), multiply by -v to get flux, then divide - ! by cell thickness, which gives a concentration tendency: - ! - ! (del/(2*dst))*(-v_coef)/del - ! - ! Simplifying gives -v_coef*rdst*0.5, as seen below. - - select case (l_bndry_loc%bndry_type) - case (fixed_layer_bndry) - op%left_bound = v_coef(:,1) / & - (l_bndry_loc%edge_width+coords%del(:,1)) - case default - op%left_bound = 0. - end select - - op%sub = v_coef(:,2:coords%d)*coords%rdst*0.5 - op%spr = -op%sub - - select case (r_bndry_loc%bndry_type) - case (fixed_layer_bndry) - op%right_bound = v_coef(:,coords%d+1) / & - (r_bndry_loc%edge_width+coords%del(:,coords%d)) - case default - op%right_bound = 0. - end select - - ! Above, we found all off-diagonals. Now get the diagonal. This must be - ! done at this specific point, since the other half of the operator is - ! not "derivative-like" in the sense of yielding 0 for a constant input. - call op%deriv_diag() - - ! The second half of the operator simply involves taking a first-order - ! derivative of v. Since v is on the interfaces, just use: - ! (v(k+1) - v(k))*rdel(k) - v_deriv(:,1) = v_coef(:,2)*coords%rdel(:,1) - - select case (l_bndry_loc%bndry_type) - case (fixed_layer_bndry) - v_deriv(:,1) = v_deriv(:,1) - v_coef(:,1)*coords%rdel(:,1) - end select - - v_deriv(:,2:coords%d-1) = (v_coef(:,3:coords%d) - & - v_coef(:,2:coords%d-1))*coords%rdel(:,2:coords%d-1) - - v_deriv(:,coords%d) = -v_coef(:,coords%d)*coords%rdel(:,coords%d) - - select case (r_bndry_loc%bndry_type) - case (fixed_layer_bndry) - v_deriv(:,coords%d) = v_deriv(:,coords%d) & - + v_coef(:,coords%d+1)*coords%del(:,coords%d) - end select - - ! Combine the two pieces. - op%diag = op%diag - v_deriv - -end function advection_operator - -! Second order approximation to the first and second derivatives on a non- -! uniform grid. -! -! Both operators are constructed with the same method, except for a "seed" -! function that takes local distances between points to create the -! off-diagonal terms. -function first_derivative(grid_spacing, l_bndry, r_bndry) result(op) - ! Distances between points. - real, USE_CONTIGUOUS intent(in) :: grid_spacing(:,:) - ! Boundary conditions. - class(BoundaryType), intent(in), optional :: l_bndry, r_bndry - ! Output operator. - type(TriDiagOp) :: op - - op = deriv_op_from_seed(grid_spacing, first_derivative_seed, & - l_bndry, r_bndry) - -end function first_derivative - -subroutine first_derivative_seed(del_minus, del_plus, sub, spr) - ! Distances to next and previous point. - real, USE_CONTIGUOUS intent(in) :: del_minus(:) - real, USE_CONTIGUOUS intent(in) :: del_plus(:) - ! Off-diagonal matrix terms. - real, USE_CONTIGUOUS intent(out) :: sub(:) - real, USE_CONTIGUOUS intent(out) :: spr(:) - - real :: del_sum(size(del_plus)) - - del_sum = del_plus + del_minus - - sub = - del_plus / (del_minus*del_sum) - spr = del_minus / (del_plus*del_sum) - -end subroutine first_derivative_seed - -function second_derivative(grid_spacing, l_bndry, r_bndry) result(op) - ! Distances between points. - real, USE_CONTIGUOUS intent(in) :: grid_spacing(:,:) - ! Boundary conditions. - class(BoundaryType), intent(in), optional :: l_bndry, r_bndry - ! Output operator. - type(TriDiagOp) :: op - - op = deriv_op_from_seed(grid_spacing, second_derivative_seed, & - l_bndry, r_bndry) - -end function second_derivative - -subroutine second_derivative_seed(del_minus, del_plus, sub, spr) - ! Distances to next and previous point. - real, USE_CONTIGUOUS intent(in) :: del_minus(:) - real, USE_CONTIGUOUS intent(in) :: del_plus(:) - ! Off-diagonal matrix terms. - real, USE_CONTIGUOUS intent(out) :: sub(:) - real, USE_CONTIGUOUS intent(out) :: spr(:) - - real :: del_sum(size(del_plus)) - - del_sum = del_plus + del_minus - - sub = 2. / (del_minus*del_sum) - spr = 2. / (del_plus*del_sum) - -end subroutine second_derivative_seed - -! Brains behind the first/second derivative functions. -function deriv_op_from_seed(grid_spacing, seed, l_bndry, r_bndry) result(op) - ! Distances between points. - real, USE_CONTIGUOUS intent(in) :: grid_spacing(:,:) - ! Function to locally construct matrix elements. - procedure(deriv_seed) :: seed - ! Boundary conditions. - class(BoundaryType), target, intent(in), optional :: l_bndry, r_bndry - ! Output operator. - type(TriDiagOp) :: op - - ! Selectors to implement default boundary. - class(BoundaryType), pointer :: l_bndry_loc, r_bndry_loc - ! Fixed flux is default, no allocation/deallocation needed. - type(BoundaryType), target :: bndry_default - - integer :: k - - if (present(l_bndry)) then - l_bndry_loc => l_bndry - else - l_bndry_loc => bndry_default - end if - - if (present(r_bndry)) then - r_bndry_loc => r_bndry - else - r_bndry_loc => bndry_default - end if - - ! Number of grid points is one greater than the spacing. - op = TriDiagOp(size(grid_spacing, 1), size(grid_spacing, 2) + 1) - - ! Left boundary condition. - call l_bndry_loc%make_left(grid_spacing, seed, & - op%left_bound, op%spr(:,1)) - - do k = 2, op%ncel-1 - call seed(grid_spacing(:,k-1), grid_spacing(:,k), & - op%sub(:,k-1), op%spr(:,k)) - end do - - ! Right boundary condition. - call r_bndry_loc%make_right(grid_spacing, seed, & - op%sub(:,op%ncel-1), op%right_bound) - - ! Above, we found all off-diagonals. Now get the diagonal. - call op%deriv_diag() - -end function deriv_op_from_seed - -! Boundary constructors. Most simply set an internal flag, but -! BoundaryFixedLayer accepts an argument representing the distance to the -! location where the extra layer is defined. - -function new_BoundaryZero() result(new_bndry) - type(BoundaryType) :: new_bndry - - new_bndry%bndry_type = zero_bndry - -end function new_BoundaryZero - -function new_BoundaryFirstOrder() result(new_bndry) - type(BoundaryType) :: new_bndry - - new_bndry%bndry_type = first_order_bndry - -end function new_BoundaryFirstOrder - -function new_BoundaryExtrapolate() result(new_bndry) - type(BoundaryType) :: new_bndry - - new_bndry%bndry_type = extrapolate_bndry - -end function new_BoundaryExtrapolate - -function new_BoundaryFixedLayer(width) result(new_bndry) - real, USE_CONTIGUOUS intent(in) :: width(:) - type(BoundaryType) :: new_bndry - - new_bndry%bndry_type = fixed_layer_bndry - new_bndry%edge_width = width - -end function new_BoundaryFixedLayer - -function new_BoundaryFixedFlux() result(new_bndry) - type(BoundaryType) :: new_bndry - - new_bndry%bndry_type = fixed_flux_bndry - -end function new_BoundaryFixedFlux - -! The make_left and make_right methods implement the boundary conditions -! using an input seed. - -subroutine make_left(self, grid_spacing, seed, term1, term2) - class(BoundaryType), intent(in) :: self - real, USE_CONTIGUOUS intent(in) :: grid_spacing(:,:) - procedure(deriv_seed) :: seed - real, USE_CONTIGUOUS intent(out) :: term1(:) - real, USE_CONTIGUOUS intent(out) :: term2(:) - - real :: del_plus(size(term1)), del_minus(size(term1)) - - select case (self%bndry_type) - case (zero_bndry) - term1 = 0. - term2 = 0. - case (first_order_bndry) - ! To calculate to first order, just use a really huge del_minus (i.e. - ! pretend that there's a point so far away it doesn't matter). - del_plus = grid_spacing(:,1) - del_minus = del_plus * 4. / epsilon(1.) - call seed(del_minus, del_plus, term1, term2) - case (extrapolate_bndry) - ! To extrapolate from the boundary, use distance from the nearest - ! neighbor (as usual) and the second nearest neighbor (with a negative - ! sign, since we are using two points on the same side). - del_plus = grid_spacing(:,1) - del_minus = - (grid_spacing(:,1) + grid_spacing(:,2)) - call seed(del_minus, del_plus, term1, term2) - case (fixed_layer_bndry) - ! Use edge value to extend the grid. - del_plus = grid_spacing(:,1) - del_minus = self%edge_width - call seed(del_minus, del_plus, term1, term2) - case (fixed_flux_bndry) - ! Treat grid as uniform, but then zero out the contribution from data - ! on one side (since it will be prescribed). - del_plus = grid_spacing(:,1) - del_minus = del_plus - call seed(del_minus, del_plus, term1, term2) - term1 = 0. - case default - call abort() - !call shr_sys_abort("Invalid boundary type at "// & - ! errMsg(__FILE__, __LINE__)) - end select - -end subroutine make_left - -subroutine make_right(self, grid_spacing, seed, term1, term2) - class(BoundaryType), intent(in) :: self - real, USE_CONTIGUOUS intent(in) :: grid_spacing(:,:) - procedure(deriv_seed) :: seed - real, USE_CONTIGUOUS intent(out) :: term1(:) - real, USE_CONTIGUOUS intent(out) :: term2(:) - - real :: del_plus(size(term1)), del_minus(size(term1)) - - select case (self%bndry_type) - case (zero_bndry) - term1 = 0. - term2 = 0. - case (first_order_bndry) - ! Use huge del_plus, analogous to how left boundary works. - del_minus = grid_spacing(:,size(grid_spacing, 2)) - del_plus = del_minus * 4. / epsilon(1.) - call seed(del_minus, del_plus, term1, term2) - case (extrapolate_bndry) - ! Same strategy as left boundary, but reversed. - del_plus = - (grid_spacing(:,size(grid_spacing, 2) - 1) + & - grid_spacing(:,size(grid_spacing, 2))) - del_minus = grid_spacing(:,size(grid_spacing, 2)) - call seed(del_minus, del_plus, term1, term2) - case (fixed_layer_bndry) - ! Use edge value to extend the grid. - del_plus = self%edge_width - del_minus = grid_spacing(:,size(grid_spacing, 2)) - call seed(del_minus, del_plus, term1, term2) - case (fixed_flux_bndry) - ! Uniform grid, but with edge zeroed. - del_plus = grid_spacing(:,size(grid_spacing, 2)) - del_minus = del_plus - call seed(del_minus, del_plus, term1, term2) - term2 = 0. - case default - call abort() - !call shr_sys_abort("Invalid boundary type at "// & - ! errMsg(__FILE__, __LINE__)) - end select - -end subroutine make_right - -subroutine boundary_type_finalize(self) - class(BoundaryType), intent(inout) :: self - - self%bndry_type = fixed_flux_bndry - if (allocated(self%edge_width)) deallocate(self%edge_width) - -end subroutine boundary_type_finalize - -! Constructor for TriDiagOp; this just sets the size and allocates -! arrays. -type(TriDiagOp) function new_TriDiagOp(nsys, ncel) - - integer, intent(in) :: nsys, ncel - - new_TriDiagOp%nsys = nsys - new_TriDiagOp%ncel = ncel - - allocate(new_TriDiagOp%spr(nsys,ncel-1), & - new_TriDiagOp%sub(nsys,ncel-1), & - new_TriDiagOp%diag(nsys,ncel), & - new_TriDiagOp%left_bound(nsys), & - new_TriDiagOp%right_bound(nsys)) - -end function new_TriDiagOp - -! Deallocator for TriDiagOp. -subroutine tridiag_finalize(self) - class(TriDiagOp), intent(inout) :: self - - self%nsys = 0 - self%ncel = 0 - - if (allocated(self%spr)) deallocate(self%spr) - if (allocated(self%sub)) deallocate(self%sub) - if (allocated(self%diag)) deallocate(self%diag) - if (allocated(self%left_bound)) deallocate(self%left_bound) - if (allocated(self%right_bound)) deallocate(self%right_bound) - -end subroutine tridiag_finalize - -! Boundary condition constructors. - -function new_BoundaryNoData() result(new_cond) - type(BoundaryCond) :: new_cond - - new_cond%cond_type = no_data_cond - ! No edge data, so leave it unallocated. - -end function new_BoundaryNoData - -function new_BoundaryData(data) result(new_cond) - real, USE_CONTIGUOUS intent(in) :: data(:) - type(BoundaryCond) :: new_cond - - new_cond%cond_type = data_cond - new_cond%edge_data = data - -end function new_BoundaryData - -function new_BoundaryFlux(flux, dt, spacing) result(new_cond) - real, USE_CONTIGUOUS intent(in) :: flux(:) - real, intent(in) :: dt - real, USE_CONTIGUOUS intent(in) :: spacing(:) - type(BoundaryCond) :: new_cond - - new_cond%cond_type = flux_cond - new_cond%edge_data = flux*dt/spacing - -end function new_BoundaryFlux - -! Application of input data. -! -! When no data is input, assume that any bound term is applied to the -! third element in from the edge for extrapolation. Boundary conditions -! that don't need any edge data at all can then simply set the boundary -! terms to 0. - -function apply_left(self, bound_term, array) result(delta_edge) - class(BoundaryCond), intent(in) :: self - real, USE_CONTIGUOUS intent(in) :: bound_term(:) - real, USE_CONTIGUOUS intent(in) :: array(:,:) - real :: delta_edge(size(array, 1)) - - select case (self%cond_type) - case (no_data_cond) - delta_edge = bound_term*array(:,3) - case (data_cond) - delta_edge = bound_term*self%edge_data - case (flux_cond) - delta_edge = self%edge_data - case default - call abort() - !call shr_sys_abort("Invalid boundary condition at "// & - ! errMsg(__FILE__, __LINE__)) - end select - -end function apply_left - -function apply_right(self, bound_term, array) result(delta_edge) - class(BoundaryCond), intent(in) :: self - real, USE_CONTIGUOUS intent(in) :: bound_term(:) - real, USE_CONTIGUOUS intent(in) :: array(:,:) - real :: delta_edge(size(array, 1)) - - select case (self%cond_type) - case (no_data_cond) - delta_edge = bound_term*array(:,size(array, 2)-2) - case (data_cond) - delta_edge = bound_term*self%edge_data - case (flux_cond) - delta_edge = self%edge_data - case default - call abort() - !call shr_sys_abort("Invalid boundary condition at "// & - ! errMsg(__FILE__, __LINE__)) - end select - -end function apply_right - -subroutine boundary_cond_finalize(self) - class(BoundaryCond), intent(inout) :: self - - self%cond_type = no_data_cond - if (allocated(self%edge_data)) deallocate(self%edge_data) - -end subroutine boundary_cond_finalize - -! Apply an operator and return the new data. -function apply_tridiag(self, array, l_cond, r_cond) result(output) - ! Operator to apply. - class(TriDiagOp), intent(in) :: self - ! Data to act on. - real, USE_CONTIGUOUS intent(in) :: array(:,:) - ! Objects representing boundary conditions. - class(BoundaryCond), target, intent(in), optional :: l_cond, r_cond - ! Function result. - real :: output(size(array, 1), size(array, 2)) - - ! Local objects to implement default. - class(BoundaryCond), pointer :: l_cond_loc, r_cond_loc - ! Default state is no data, no allocation/deallocation needed. - type(BoundaryCond), target :: cond_default - - ! Level index. - integer :: k - - if (present(l_cond)) then - l_cond_loc => l_cond - else - l_cond_loc => cond_default - end if - - if (present(r_cond)) then - r_cond_loc => r_cond - else - r_cond_loc => cond_default - end if - - ! Left boundary. - output(:,1) = self%diag(:,1)*array(:,1) + & - self%spr(:,1)*array(:,2) + & - l_cond_loc%apply_left(self%left_bound, array) - - do k = 2, self%ncel-1 - output(:,k) = & - self%sub(:,k-1)*array(:,k-1) + & - self%diag(:,k)*array(:,k ) + & - self%spr(:,k)*array(:,k+1) - end do - - ! Right boundary. - output(:,self%ncel) = & - self%sub(:,self%ncel-1)*array(:,self%ncel-1) + & - self%diag(:,self%ncel)*array(:,self%ncel) + & - r_cond_loc%apply_right(self%right_bound, array) - -end function apply_tridiag - -! Fill in the diagonal for a TriDiagOp for a derivative operator, where -! the off diagonal elements are already filled in. -subroutine make_tridiag_deriv_diag(self) - - class(TriDiagOp), intent(inout) :: self - - ! If a derivative operator operates on a constant function, it must - ! return 0 everywhere. To force this, make sure that all rows add to - ! zero in the matrix. - self%diag(:,:self%ncel-1) = - self%spr - self%diag(:,self%ncel) = - self%right_bound - self%diag(:,1) = self%diag(:,1) - self%left_bound - self%diag(:,2:) = self%diag(:,2:) - self%sub - -end subroutine make_tridiag_deriv_diag - -! Sum two TriDiagOp objects into a new one; this is just the addition of -! all the entries. -function add_tridiag_ops(op1, op2) result(new_op) - - type(TriDiagOp), intent(in) :: op1, op2 - type(TriDiagOp) :: new_op - - new_op = op1 - - call new_op%add(op2) - -end function add_tridiag_ops - -subroutine add_in_place_tridiag_ops(self, other) - - class(TriDiagOp), intent(inout) :: self - class(TriDiagOp), intent(in) :: other - - self%spr = self%spr + other%spr - self%sub = self%sub + other%sub - self%diag = self%diag + other%diag - - self%left_bound = self%left_bound + other%left_bound - self%right_bound = self%right_bound + other%right_bound - -end subroutine add_in_place_tridiag_ops - -! Subtract two TriDiagOp objects. -function subtract_tridiag_ops(op1, op2) result(new_op) - - type(TriDiagOp), intent(in) :: op1, op2 - type(TriDiagOp) :: new_op - - new_op = op1 - - call new_op%subtract(op2) - -end function subtract_tridiag_ops - -! Subtract two TriDiagOp objects. -subroutine subtract_in_place_tridiag_ops(self, other) - - class(TriDiagOp), intent(inout) :: self - class(TriDiagOp), intent(in) :: other - - self%spr = self%spr - other%spr - self%sub = self%sub - other%sub - self%diag = self%diag - other%diag - - self%left_bound = self%left_bound - other%left_bound - self%right_bound = self%right_bound - other%right_bound - -end subroutine subtract_in_place_tridiag_ops - -! Equivalent to adding a multiple of the identity. -subroutine scalar_add_tridiag(self, constant) - - class(TriDiagOp), intent(inout) :: self - real, intent(in) :: constant - - self%diag = self%diag + constant - -end subroutine scalar_add_tridiag - -! Equivalent to adding the diagonal operator constructed from diag_array. -subroutine diagonal_add_tridiag(self, diag_array) - - class(TriDiagOp), intent(inout) :: self - real, USE_CONTIGUOUS intent(in) :: diag_array(:,:) - - self%diag = self%diag + diag_array - -end subroutine diagonal_add_tridiag - -! Multiply a scalar by an array. -subroutine scalar_lmult_tridiag(self, constant) - - class(TriDiagOp), intent(inout) :: self - real, intent(in) :: constant - - self%spr = self%spr * constant - self%sub = self%sub * constant - self%diag = self%diag * constant - - self%left_bound = self%left_bound * constant - self%right_bound = self%right_bound * constant - -end subroutine scalar_lmult_tridiag - -! Multiply in an array as if it contained the entries of a diagonal matrix -! being multiplied from the left. -subroutine diagonal_lmult_tridiag(self, diag_array) - - class(TriDiagOp), intent(inout) :: self - real, USE_CONTIGUOUS intent(in) :: diag_array(:,:) - - self%spr = self%spr * diag_array(:,:self%ncel-1) - self%sub = self%sub * diag_array(:,2:) - self%diag = self%diag * diag_array(:,:) - - self%left_bound = self%left_bound * diag_array(:,1) - self%right_bound = self%right_bound * diag_array(:,self%ncel) - -end subroutine diagonal_lmult_tridiag - -! Decomposition constructor -! -! The equation to be solved later (with left_div) is: -! - A(k)*q(k+1) + B(k)*q(k) - C(k)*q(k-1) = D(k) -! -! The solution (effectively via LU decomposition) has the form: -! E(k) = C(k) / (B(k) - A(k)*E(k+1)) -! F(k) = (D(k) + A(k)*F(k+1)) / (B(k) - A(k)*E(k+1)) -! q(k) = E(k) * q(k-1) + F(k) -! -! Unlike Richtmyer and Morton, E and F are defined by iterating backward -! down to level 1, and then q iterates forward. -! -! E can be calculated and stored now, without knowing D. -! To calculate F later, we store A and the denominator. -function new_TriDiagDecomp(op, graft_decomp) result(decomp) - type(TriDiagOp), intent(in) :: op - type(TriDiagDecomp), intent(in), optional :: graft_decomp - - type(TriDiagDecomp) :: decomp - - integer :: k - - if (present(graft_decomp)) then - decomp%nsys = graft_decomp%nsys - decomp%ncel = graft_decomp%ncel - else - decomp%nsys = op%nsys - decomp%ncel = op%ncel - end if - - ! Simple allocation with no error checking. - allocate(decomp%ca(decomp%nsys,decomp%ncel)) - allocate(decomp%dnom(decomp%nsys,decomp%ncel)) - allocate(decomp%ze(decomp%nsys,decomp%ncel)) - - ! decomp%ca is simply the negative of the tridiagonal's superdiagonal. - decomp%ca(:,:op%ncel-1) = -op%spr - decomp%ca(:,op%ncel) = -op%right_bound - - if (present(graft_decomp)) then - ! Copy in graft_decomp beyond op%ncel. - decomp%ca(:,op%ncel+1:) = graft_decomp%ca(:,op%ncel+1:) - decomp%dnom(:,op%ncel+1:) = graft_decomp%dnom(:,op%ncel+1:) - decomp%ze(:,op%ncel+1:) = graft_decomp%ze(:,op%ncel+1:) - ! Fill in dnom edge value. - decomp%dnom(:,op%ncel) = 1. / (op%diag(:,op%ncel) - & - decomp%ca(:,op%ncel)*decomp%ze(:,op%ncel+1)) - else - ! If no grafting, the edge value of dnom comes from the diagonal. - decomp%dnom(:,op%ncel) = 1. / op%diag(:,op%ncel) - end if - - do k = op%ncel - 1, 1, -1 - decomp%ze(:,k+1) = - op%sub(:,k) * decomp%dnom(:,k+1) - decomp%dnom(:,k) = 1. / & - (op%diag(:,k) - decomp%ca(:,k)*decomp%ze(:,k+1)) - end do - - ! Don't multiply edge level by denom, because we want to leave it up to - ! the BoundaryCond object to decide what this means in left_div. - decomp%ze(:,1) = -op%left_bound - -end function new_TriDiagDecomp - -! Left-division (multiplication by inverse) using a decomposed operator. -! -! See the comment above for the constructor for a quick explanation of the -! intermediate variables. The "q" argument is "D(k)" on input and "q(k)" on -! output. -subroutine decomp_left_div(decomp, q, l_cond, r_cond) - - ! Decomposed matrix. - class(TriDiagDecomp), intent(in) :: decomp - ! Data to left-divide by the matrix. - real, USE_CONTIGUOUS intent(inout) :: q(:,:) - ! Objects representing boundary conditions. - class(BoundaryCond), intent(in), optional :: l_cond, r_cond - - ! "F" from the equation above. - real :: zf(decomp%nsys,decomp%ncel) - - ! Level index. - integer :: k - - ! Include boundary conditions. - if (present(l_cond)) then - q(:,1) = q(:,1) + l_cond%apply_left(decomp%ze(:,1), q) - end if - - if (present(r_cond)) then - q(:,decomp%ncel) = q(:,decomp%ncel) + & - r_cond%apply_right(decomp%ca(:,decomp%ncel), q) - end if - - zf(:,decomp%ncel) = q(:,decomp%ncel) * decomp%dnom(:,decomp%ncel) - - do k = decomp%ncel - 1, 1, -1 - zf(:,k) = (q(:,k) + decomp%ca(:,k)*zf(:,k+1)) * decomp%dnom(:,k) - end do - - ! Perform back substitution - - q(:,1) = zf(:,1) - - do k = 2, decomp%ncel - q(:,k) = zf(:,k) + decomp%ze(:,k)*q(:,k-1) - end do - -end subroutine decomp_left_div - -! Decomposition deallocation. -subroutine decomp_finalize(decomp) - class(TriDiagDecomp), intent(inout) :: decomp - - decomp%nsys = 0 - decomp%ncel = 0 - - if (allocated(decomp%ca)) deallocate(decomp%ca) - if (allocated(decomp%dnom)) deallocate(decomp%dnom) - if (allocated(decomp%ze)) deallocate(decomp%ze) - -end subroutine decomp_finalize - -end module linear_1d_operators diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/vdiff_lu_solver.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/vdiff_lu_solver.F90 deleted file mode 100644 index a593ba268..000000000 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/ncar_gwd/vdiff_lu_solver.F90 +++ /dev/null @@ -1,203 +0,0 @@ -module vdiff_lu_solver -#define USE_CONTIGUOUS contiguous, -! This module provides a function returning the matrix decomposition for -! an implicit finite volume solver for vertical diffusion. It accepts -! diffusion coefficients, time/grid spacing, and boundary condition -! objects, and returns a TriDiagDecomp object that can be used to diffuse -! an array for one time step with the "left_div" method. - -use coords_1d, only: Coords1D -use linear_1d_operators, only: TriDiagOp, operator(+), TriDiagDecomp - -implicit none -private - -! Public interfaces -public :: vd_lu_decomp -public :: fin_vol_lu_decomp - -contains - -! ========================================================================! - -! Designed to solve the equation: -! dq/dt = c1 q'' + c2 q' + c q - -function vd_lu_decomp(dt, dp, coef_q, coef_q_d, coef_q_d2, upper_bndry, & - lower_bndry) result(decomp) - - use linear_1d_operators, only: & - identity_operator, & - diagonal_operator, & - first_derivative, & - second_derivative, & - BoundaryType - - ! ---------------------- ! - ! Input-Output Arguments ! - ! ---------------------- ! - - ! Time step. - real, intent(in) :: dt - ! Grid spacing (deltas). - real, USE_CONTIGUOUS intent(in) :: dp(:,:) - - ! Coefficients for q, q', and q''. - real, USE_CONTIGUOUS intent(in), optional :: coef_q(:,:), & - coef_q_d(:,:), coef_q_d2(:,:) - - ! Boundary conditions (optional, default to 0 flux through boundary). - class(BoundaryType), target, intent(in), optional :: & - upper_bndry, lower_bndry - - ! Output decomposition. - type(TriDiagDecomp) :: decomp - - ! --------------- ! - ! Local Variables ! - ! --------------- ! - - ! Operator objects. - type(TriDiagOp) :: add_term - type(TriDiagOp) :: net_operator - - ! ----------------------- ! - ! Main Computation Begins ! - ! ----------------------- ! - - if (present(coef_q)) then - net_operator = diagonal_operator(1. - dt*coef_q) - else - net_operator = identity_operator(size(dp, 1), size(dp, 2) + 1) - end if - - if (present(coef_q_d)) then - add_term = first_derivative(dp, upper_bndry, lower_bndry) - call add_term%lmult_as_diag(-dt*coef_q_d) - call net_operator%add(add_term) - end if - - if (present(coef_q_d2)) then - add_term = second_derivative(dp, upper_bndry, lower_bndry) - call add_term%lmult_as_diag(-dt*coef_q_d2) - call net_operator%add(add_term) - end if - - decomp = TriDiagDecomp(net_operator) - - call net_operator%finalize() - call add_term%finalize() - -end function vd_lu_decomp - -! ========================================================================! - -! Designed to solve the equation: -! -! w * dq/dt = d/dp (D q' - v q) + c q -! -! where q is a grid-cell average, and p is the vertical coordinate -! (presumably pressure). -! -! In this function, coef_q_weight == w, coef_q_diff == D, -! coef_q_adv == v, and coef_q == c. All these are optional; omitting a -! coefficient is equivalent to setting the entire array to 0. -! -! coef_q_diff and coef_q_adv are defined at the level interfaces, while -! coef_q and coef_q_weight are grid-cell averages. - -function fin_vol_lu_decomp(dt, p, coef_q, coef_q_diff, coef_q_adv, & - coef_q_weight, upper_bndry, lower_bndry, graft_decomp) result(decomp) - - use linear_1d_operators, only: & - zero_operator, & - diagonal_operator, & - diffusion_operator, & - advection_operator, & - BoundaryType - - ! ---------------------- ! - ! Input-Output Arguments ! - ! ---------------------- ! - - ! Time step. - real, intent(in) :: dt - ! Grid spacings. - type(Coords1D), intent(in) :: p - - ! Coefficients for diffusion and advection. - ! - ! The sizes must be consistent among all the coefficients that are - ! actually present, i.e. coef_q_diff and coef_q_adv should be one level - ! bigger than coef_q and coef_q_weight, and have the same column number. - real, USE_CONTIGUOUS intent(in), optional :: coef_q(:,:), & - coef_q_diff(:,:), coef_q_adv(:,:), coef_q_weight(:,:) - - ! Boundary conditions (optional, default to 0 flux through boundary). - class(BoundaryType), target, intent(in), optional :: & - upper_bndry, lower_bndry - - ! Decomposition to graft onto. If this is provided, you can pass in - ! smaller coefficients. - type(TriDiagDecomp), intent(in), optional :: graft_decomp - - ! Output decomposition. - type(TriDiagDecomp) :: decomp - - ! --------------- ! - ! Local Variables ! - ! --------------- ! - - ! Operator objects. - type(TriDiagOp) :: add_term - type(TriDiagOp) :: net_operator - - ! ----------------------- ! - ! Main Computation Begins ! - ! ----------------------- ! - - ! A diffusion term is probably present, so start with that. Otherwise - ! start with an operator of all 0s. - - if (present(coef_q_diff)) then - net_operator = diffusion_operator(p, coef_q_diff, & - upper_bndry, lower_bndry) - else - net_operator = zero_operator(p%n, p%d) - end if - - ! Constant term (damping). - if (present(coef_q)) then - add_term = diagonal_operator(coef_q) - call net_operator%add(add_term) - end if - - ! Effective advection. - if (present(coef_q_adv)) then - add_term = advection_operator(p, coef_q_adv, & - upper_bndry, lower_bndry) - call net_operator%add(add_term) - end if - - ! We want I-dt*(w^-1)*A for a single time step, implicit method, where - ! A is the right-hand-side operator (i.e. what net_operator is now). - if (present(coef_q_weight)) then - call net_operator%lmult_as_diag(-dt/coef_q_weight) - else - call net_operator%lmult_as_diag(-dt) - end if - call net_operator%add_to_diag(1.) - - ! Decompose, grafting on an optional input decomp. The graft is a way to - ! avoid re-calculating the ending (bottom) levels when the coefficients - ! have only changed at the beginning (top), e.g. for different - ! constituents in the molecular diffusion. - decomp = TriDiagDecomp(net_operator, graft_decomp=graft_decomp) - - ! Ensure local objects are deallocated. - call net_operator%finalize() - call add_term%finalize() - -end function fin_vol_lu_decomp - -end module vdiff_lu_solver From 9f600187e1ffe8d9745269441cfa7d00f05db91d Mon Sep 17 00:00:00 2001 From: Matthew Thompson Date: Thu, 30 Nov 2023 14:57:10 -0500 Subject: [PATCH 2/2] Fix CMake --- .../GEOSphysics_GridComp/GEOSgwd_GridComp/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/CMakeLists.txt b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/CMakeLists.txt index b986152dc..841b854d4 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/CMakeLists.txt +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/CMakeLists.txt @@ -1,6 +1,6 @@ esma_set_this () -set subdirs( +set (subdirs ncar_gwd )