diff --git a/.circleci/config.yml b/.circleci/config.yml index 3c95aa11..e4568429 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -1,92 +1,64 @@ version: 2.1 -executors: - gcc-build-env: - docker: - - image: gmao/ubuntu20-geos-env-mkl:v6.2.8-openmpi_4.0.6-gcc_11.2.0 - auth: - username: $DOCKERHUB_USER - password: $DOCKERHUB_AUTH_TOKEN - environment: - OMPI_ALLOW_RUN_AS_ROOT: 1 - OMPI_ALLOW_RUN_AS_ROOT_CONFIRM: 1 - OMPI_MCA_btl_vader_single_copy_mechanism: none - resource_class: large - #MEDIUM# resource_class: medium +orbs: + circleci-tools: geos-esm/circleci-tools@0.13.0 workflows: - version: 2.1 build-and-test: jobs: - build-GOCART2G: + name: build-GOCART2G-on-<< matrix.compiler >> + matrix: + parameters: + compiler: [gfortran, ifort] context: - docker-hub-creds - build-GEOSgcm: + name: build-GEOSgcm-on-<< matrix.compiler >> + matrix: + parameters: + compiler: [gfortran, ifort] context: - docker-hub-creds jobs: build-GOCART2G: - executor: gcc-build-env + parameters: + compiler: + type: string + executor: circleci-tools/<< parameters.compiler >> working_directory: /root/project steps: - checkout: path: GOCART - - run: - name: "Versions etc" - command: mpirun --version && gfortran --version && echo $BASEDIR && pwd && ls - - run: - name: "Mepo clone external repos" - command: | - cd ${CIRCLE_WORKING_DIRECTORY}/GOCART - mepo clone - mepo status - - run: - name: "CMake" - command: | - cd ${CIRCLE_WORKING_DIRECTORY}/GOCART - mkdir build - cd build - cmake .. -DBASEDIR=$BASEDIR/Linux -DCMAKE_Fortran_COMPILER=gfortran -DCMAKE_BUILD_TYPE=Debug -DUSE_F2PY=OFF -DMPIEXEC_PREFLAGS='--oversubscribe' - - run: - name: "Build GOCART2G_GridComp" - command: | - cd ${CIRCLE_WORKING_DIRECTORY}/GOCART/build - make -j"$(nproc)" GOCART2G_GridComp - #MEDIUM# make -j4 GOCART2G_GridComp + - circleci-tools/versions: + compiler: << parameters.compiler >> + - circleci-tools/mepoclone: + repo: GOCART + - circleci-tools/cmake: + repo: GOCART + compiler: << parameters.compiler >> + - circleci-tools/buildtarget: + repo: GOCART + target: GOCART2G_GridComp + - circleci-tools/compress_artifacts + - store_artifacts: + path: /logfiles build-GEOSgcm: - executor: gcc-build-env + parameters: + compiler: + type: string + executor: circleci-tools/<< parameters.compiler >> working_directory: /root/project steps: - - run: - name: "Checkout GEOSgcm fixture" - command: | - cd ${CIRCLE_WORKING_DIRECTORY} - git clone https://github.com/GEOS-ESM/GEOSgcm.git - - run: - name: "Mepo clone external repos" - command: | - cd ${CIRCLE_WORKING_DIRECTORY}/GEOSgcm - mepo clone - mepo develop GEOSgcm_GridComp GEOSgcm_App - mepo status - - run: - name: "Mepo checkout GOCART branch(es)" - command: | - cd ${CIRCLE_WORKING_DIRECTORY}/GEOSgcm - mepo checkout-if-exists ${CIRCLE_BRANCH} - mepo status - - run: - name: "CMake" - command: | - cd ${CIRCLE_WORKING_DIRECTORY}/GEOSgcm - mkdir build - cd build - cmake .. -DBASEDIR=$BASEDIR/Linux -DCMAKE_Fortran_COMPILER=gfortran -DCMAKE_BUILD_TYPE=Debug -DCMAKE_INSTALL_PREFIX=${CIRCLE_WORKING_DIRECTORY}/workspace/install-GEOSgcm -DUSE_F2PY=OFF - - run: - name: "Build and install" - command: | - cd ${CIRCLE_WORKING_DIRECTORY}/GEOSgcm/build - make -j"$(nproc)" install - #MEDIUM# make -j4 install + - circleci-tools/checkout_fixture + - circleci-tools/mepoclone + - circleci-tools/mepodevelop + - circleci-tools/checkout_if_exists + - circleci-tools/cmake: + compiler: << parameters.compiler >> + - circleci-tools/buildinstall + - circleci-tools/compress_artifacts + - store_artifacts: + path: /logfiles diff --git a/.editorconfig b/.editorconfig new file mode 100644 index 00000000..9f1cd658 --- /dev/null +++ b/.editorconfig @@ -0,0 +1,32 @@ +# Global Editor Config for MAPL +# +# This is an ini style configuration. See http://editorconfig.org/ for more information on this file. +# +# Top level editor config. +root = true + +# Always use Unix style new lines with new line ending on every file and trim whitespace +[*] +end_of_line = lf +insert_final_newline = true +trim_trailing_whitespace = true + +# Python: PEP8 defines 4 spaces for indentation +[*.py] +indent_style = space +indent_size = 4 + +# YAML format, 2 spaces +[{*.yaml,*.yml}] +indent_style = space +indent_size = 2 + +# CMake (from KitWare: https://github.com/Kitware/CMake/blob/master/.editorconfig) +[{CMakeLists.txt,*.cmake,*.rst}] +indent_style = space +indent_size = 2 + +# Markdown +[*.md] +trim_trailing_whitespace = true +indent_style = space diff --git a/.github/CODEOWNERS b/.github/CODEOWNERS index 47d7ab45..67d9db28 100644 --- a/.github/CODEOWNERS +++ b/.github/CODEOWNERS @@ -7,28 +7,12 @@ * @GEOS-ESM/aerosol-team # NOAA Extras are Tom, Weiyuan for now -/CCPP/ @tclune @weiyuan-jiang -/ESMF/Aerosol_GridComp/ @tclune @weiyuan-jiang -/ESMF/NUOPC/ @tclune @weiyuan-jiang -/ESMF/Shared/ @tclune @weiyuan-jiang -/ESMF/UFS/ @tclune @weiyuan-jiang - -/ESMF/GOCART_GridComp/ @GEOS-ESM/aerosol-team - -# GOCART2G is Elliot for now -/ESMF/GOCART2G_GridComp/ @GEOS-ESM/aerosol-team +/CCPP/ @GEOS-ESM/mapl-team @amdasilva +/ESMF/Aerosol_GridComp/ @GEOS-ESM/mapl-team @amdasilva +/ESMF/NUOPC/ @GEOS-ESM/mapl-team @amdasilva +/ESMF/UFS/ @GEOS-ESM/mapl-team @amdasilva # The CMake Team should know/approve these -/.github/ @GEOS-ESM/cmake-team -/.circleci/ @GEOS-ESM/cmake-team -/.codebuild/ @GEOS-ESM/cmake-team - -# The GEOS CMake Team should be notified about and approve config changes -/config/ @GEOS-ESM/cmake-team - -# The Python Transition Team will own Python files -# until the Python 3 transition is completed -*.py @GEOS-ESM/python-transition-team - -# The GEOS CMake Team is the CODEOWNER for the CMakeLists.txt files in this repository -CMakeLists.txt @GEOS-ESM/cmake-team +/.github/ @GEOS-ESM/cmake-team @amdasilva +/.circleci/ @GEOS-ESM/cmake-team @amdasilva +/.codebuild/ @GEOS-ESM/cmake-team @amdasilva diff --git a/CHANGELOG.md b/CHANGELOG.md index 614933c5..7189343c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,29 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +### Added + +- Add `.editorconfig` file + - This matches the styles currently used in MAPL (2 space indents in CMake and yaml, 4 spaces for Python) + +### Changed + +- Update `CODEOWNERS` file to make approvals less restrictive +- Updated the CircleCI to use circleci-tools 0.13.0 orb + - Moves CI to use Baselibs 6.2.13 needed by MAPL development +- Update `components.yaml` to be in line with GEOSgcm v10.22.1 +- Updates to support Spack +- Changed the handling of state variable names in multiple instances of component (see Issue #93) +- Major refactoring of Mie table class. (see Issue #96) + - Renamed Chem_MieTableMod.F90 --> GOCART2G_Mie2GMod.F90 + - Renamed module Chem_MieTableMod2G --> GOCART2G_Mie2GMod + - Introduced object oriented design with type-bound methods + - renamed some components/arguments for clarity + - eliminated extraneous container data type that is not needed under new GOCART design. +- Cleaned up optional keyword arguments in call to Mie calculator for aerosol + radiative forcing calculation; zero diff change +- Simplified loading of radiation MieTables. + ## [2.0.6] - 2021-04-28 ### Fixed @@ -40,8 +63,6 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Changed -- Cleaned up optional keyword arguments in call to Mie calculator for aerosol - radiative forcing calculation; zero diff change - Updated FENGSHA dust flux according to Webb et al., Aeolian Res. 42 (2020) 100560 ## [2.0.2] - 2021-01-07 diff --git a/ESMF/Apps/GOCART2G_AopMod.F90 b/ESMF/Apps/GOCART2G_AopMod.F90 new file mode 100644 index 00000000..efba6aff --- /dev/null +++ b/ESMF/Apps/GOCART2G_AopMod.F90 @@ -0,0 +1,452 @@ + +#include "MAPL_Generic.h" + +!------------------------------------------------------------------------- +! NASA/GSFC, Global Modeling & Assimilation Office, Code 900.3 ! +!------------------------------------------------------------------------- +!BOP +! +! !MODULE: Chem_AodMod --- Aerosol Optical Depth Calculator +! +! !INTERFACE: +! + + module Chem_AodMod + +! !USES: + + + use ESMF + use MAPL + + use Chem_MieTableMod + use Chem_MieMod + use Chem_RegistryMod + + implicit none + +! +! !PUBLIC MEMBER FUNCTIONS: +! + PUBLIC Chem_AodCalculator ! Calculates AOD given mixing ratios + PUBLIC Chem_rEffCalculator ! Calculates Aerosol Effective Radius + +! +! !DESCRIPTION: +! +! This module implements an Aerosol Optical Depth calculator. +! +! !REVISION HISTORY: +! +! 29Nov2010 da Silva Revamped to use Simple Bundle construct. +! 12Jan2011 da Silva Reff capability. +! +!EOP +!------------------------------------------------------------------------- + + interface Chem_AodCalculator + Module Procedure Chem_AodCalculator2d + end interface Chem_AodCalculator + +CONTAINS + +!----------------------------------------------------------------------------- +!BOP + +! !IROUTINE: Chem_AodCalculator ---- Mie Calculator given Aerosol Mixing Ratio +! +! !INTERFACE: +! + subroutine Chem_AodCalculator2d ( y, x, Mie, verbose, rc ) + + +! !ARGUMENTS: + + type(MAPL_SimpleBundle), intent(inout) :: y ! Already created Bundle + type(MAPL_SimpleBundle), intent(in) :: x ! Aerosol concentrations + type(Chem_Mie), intent(inout) :: Mie ! Mie tables, etc + + logical, OPTIONAL, intent(in) :: verbose + integer, OPTIONAL, intent(out) :: rc + +! !DESCRIPTION: Compute AOD given aerosol mixing ratio. +! +!EOP +!----------------------------------------------------------------------------- + + integer iq, idxTable, im, jm, km, iRH, i, j, k, n + real(ESMF_KIND_R4) :: tau, delm, rh + logical :: verbose_ + + __Iam__("Chem_AodCalculator") + + if ( present(verbose) ) then + verbose_ = verbose .and. MAPL_AM_I_ROOT() + else + verbose_ = .FALSE. + end if + + im = size(x%coords%lons,1) + jm = size(x%coords%lons,2) + km = size(x%coords%levs) + + iRH = MAPL_SimpleBundleGetIndex(x,'RH',rank=3,rc=STATUS,quiet=.TRUE.) + if ( STATUS /= 0 ) then + iRH = MAPL_SimpleBundleGetIndex(x,'RH2',rank=3,__RC__) + end if + +! Consistency check +! ----------------- + if ( Mie%nch /= size(y%coords%levs) ) then + __raise__(MAPL_RC_ERROR,"mieTable/y has inconsistent number of channels") + end if + if ( .not. associated(x%coords%lcv%delp) ) then + __raise__(MAPL_RC_ERROR,"x%coords%lcv must have delp set") + end if + +! Initialize output arrays to zero +! -------------------------------- + y%r3(1)%q = 0.0 + +! Loop over aerosol species +! ------------------------- + do iq = 1, x%n3d + + idxTable = Chem_MieQueryIdx(Mie,x%r3(iq)%name,rc) + + if (idxTable == -1) then + if (verbose_) & + print *, '[-] Skipping '//trim(x%r3(iq)%name)//' contribution' + + cycle + end if + + if ( rc/=0 ) then + __raise__(MAPL_RC_ERROR,"cannot get Mie index for "//trim(x%r3(iq)%name)) + end if + + if (verbose_) & + print *, '[+] Adding '//trim(x%r3(iq)%name)//' contribution' + +! Loop over channel, x, y, z +! -------------------------- + do n = 1, Mie%nch + do k = 1, km + do j = 1, jm + do i = 1, im + + delm = x%r3(iq)%q(i,j,k) * x%coords%lcv%delp(i,j,k) / MAPL_GRAV + rh = x%r3(iRH)%q(i,j,k) + + call Chem_MieQuery(Mie, idxTable, float(n), delm, rh, tau=tau) + + y%r3(1)%q(i,j,n) = y%r3(1)%q(i,j,n) + tau + + end do ! longitudes + end do ! latitudes + end do ! levels + end do ! channels + + end do ! aerosol tracers + + + if (verbose_) & + print *, '[x] All done!' + + rc = 0 + + end subroutine Chem_AodCalculator2d + +!----------------------------------------------------------------------------- +!BOP + +! !IROUTINE: Chem_rEffCalculator ---- Compute Aerosol Effective Radius +! +! !INTERFACE: +! + subroutine Chem_rEffCalculator ( y, x, Mie, verbose, rc ) + + +! !ARGUMENTS: + + type(MAPL_SimpleBundle), intent(inout) :: y ! Already created Bundle + type(MAPL_SimpleBundle), intent(in) :: x ! Aerosol concentrations + type(Chem_Mie), intent(inout) :: Mie ! Mie tables, etc + + logical, OPTIONAL, intent(in) :: verbose + integer, OPTIONAL, intent(out) :: rc + +! !DESCRIPTION: Return Aerosol Effective Radius. +! +!EOP +!----------------------------------------------------------------------------- + + integer iq, idxTable, im, jm, km, iRH, i, j, k + real(ESMF_KIND_R4) :: rEff, delm, rh + logical :: verbose_ + + __Iam__("Chem_rEffCalculator") + + if ( present(verbose) ) then + verbose_ = verbose .and. MAPL_AM_I_ROOT() + else + verbose_ = .FALSE. + end if + + im = size(x%coords%lons,1) + jm = size(x%coords%lons,2) + km = size(x%coords%levs) + + iRH = MAPL_SimpleBundleGetIndex(x,'rh',rank=3,quiet=.true.,rc=STATUS) + if ( STATUS /= 0 ) then + iRH = MAPL_SimpleBundleGetIndex(x,'RH',rank=3,__RC__) + end if + +! Consistency check +! ----------------- + if ( .not. associated(x%coords%lcv%delp) ) then + __raise__(MAPL_RC_ERROR,"x%coords%lcv must have delp set") + end if + + delm = 0.0 + +! Loop over aerosol species +! ------------------------- + do iq = 1, x%n3d + + idxTable = Chem_MieQueryIdx(Mie,x%r3(iq)%name,rc) + if(idxTable == -1) cycle + if ( rc/=0 ) then + __raise__(MAPL_RC_ERROR,"cannot get Mie index for "//trim(x%r3(iq)%name)) + end if + + if (verbose_) & + print *, '[+] Adding '//trim(x%r3(iq)%name)//' contribution' + +! Loop over x, y, z +! ----------------- + do k = 1, km + do j = 1, jm + do i = 1, im + + rh = x%r3(iRH)%q(i,j,k) + + call Chem_MieQuery(Mie, idxTable, 1.0, delm, rh, rEff=rEff) + + y%r3(iq)%q(i,j,k) = rEff + + + end do ! longitudes + end do ! latitudes + end do ! levels + + end do ! aerosol tracers + + + if (verbose_) & + print *, '[x] All done!' + + rc = 0 + + end subroutine Chem_rEffCalculator + +!----------------------------------------------------------------------------- +!BOP + +! !IROUTINE: Chem_ConcCalculator ---- Compute Aerosol Effective Radius +! +! !INTERFACE: +! + subroutine Chem_ConcCalculator ( y, x, Mie, verbose, rc ) + + +! !ARGUMENTS: + + type(MAPL_SimpleBundle), intent(inout) :: y ! Already created Bundle + type(MAPL_SimpleBundle), intent(in) :: x ! Aerosol concentrations + type(Chem_Mie), intent(inout) :: Mie ! Mie tables, etc + + logical, OPTIONAL, intent(in) :: verbose + integer, OPTIONAL, intent(out) :: rc + +! !DESCRIPTION: Return Aerosol Concentration bundle. +! +!EOP +!----------------------------------------------------------------------------- + + integer iq, im, jm, km, i, j, k, idxTable + real(ESMF_KIND_R4) :: delm + logical :: verbose_ + + _Iam_("Chem_ConcCalculator") + + if ( present(verbose) ) then + verbose_ = verbose .and. MAPL_AM_I_ROOT() + else + verbose_ = .FALSE. + end if + + im = size(x%coords%lons,1) + jm = size(x%coords%lons,2) + km = size(x%coords%levs) + +! Consistency check +! ----------------- + if ( .not. associated(x%coords%lcv%delp) ) then + __raise__(MAPL_RC_ERROR,"x%coords%lcv must have delp set") + end if + +! Loop over aerosol species +! ------------------------- + do iq = 1, x%n3d + + idxTable = Chem_MieQueryIdx(Mie,x%r3(iq)%name,rc) + if(idxTable == -1) cycle + if ( rc/=0 ) then + __raise__(MAPL_RC_ERROR,"cannot get Mie index for "//trim(x%r3(iq)%name)) + end if + + if (verbose_) & + print *, '[+] Adding '//trim(x%r3(iq)%name)//' contribution' + +! Loop over x, y, z +! ----------------- + do k = 1, km + do j = 1, jm + do i = 1, im + + delm = x%r3(iq)%q(i,j,k) * x%coords%lcv%delp(i,j,k) / MAPL_GRAV + + y%r3(iq)%q(i,j,k) = delm + + end do ! longitudes + end do ! latitudes + end do ! levels + + end do ! aerosol tracers + + + if (verbose_) & + print *, '[x] All done!' + + rc = 0 + + end subroutine Chem_ConcCalculator + +!----------------------------------------------------------------------------- +!BOP + +! !IROUTINE: Chem_ExtCalculator ---- Compute 3D Aerosol Extinction +! +! !INTERFACE: +! + subroutine Chem_ExtCalculator ( y, x, Mie, verbose, rc ) + + +! !ARGUMENTS: + + type(MAPL_SimpleBundle), intent(inout) :: y ! Aerosol extinction + type(MAPL_SimpleBundle), intent(in) :: x ! Aerosol mixing ratio + type(Chem_Mie), intent(inout) :: Mie ! Mie tables, etc + + logical, OPTIONAL, intent(in) :: verbose + integer, OPTIONAL, intent(out) :: rc + +! !DESCRIPTION: Return Aerosol Extinction bundle. +! +!EOP +!----------------------------------------------------------------------------- + + integer :: iq, im, jm, km, i, j, k, n, idxTable + integer :: iTau, iExt, iBck, iRH, iRHO + real(ESMF_KIND_R4) :: delc, delm, delz, rh, tau, bck + logical :: verbose_ + + __Iam__("Chem_ExtCalculator") + + if ( present(verbose) ) then + verbose_ = verbose .and. MAPL_AM_I_ROOT() + else + verbose_ = .FALSE. + end if + + im = size(x%coords%lons,1) + jm = size(x%coords%lons,2) + km = size(x%coords%levs) + + iRH = MAPL_SimpleBundleGetIndex(x,'rh',rank=3,quiet=.true.,rc=STATUS) + if ( STATUS /= 0 ) then + iRH = MAPL_SimpleBundleGetIndex(x,'RH',rank=3,__RC__) + end if + + iRHO = MAPL_SimpleBundleGetIndex(x,'airdens',rank=3,quiet=.true.,rc=STATUS) + if ( STATUS /= 0 ) then + iRHO = MAPL_SimpleBundleGetIndex(x,'AIRDENS',rank=3,__RC__) + end if + + iTau = MAPL_SimpleBundleGetIndex(y,'tau', rank=3,quiet=.true.,rc=STATUS) + iExt = MAPL_SimpleBundleGetIndex(y,'extinction',rank=3,quiet=.true.,rc=STATUS) + iBck = MAPL_SimpleBundleGetIndex(y,'backscat', rank=3,quiet=.true.,rc=STATUS) + +! Consistency check +! ----------------- + if ( .not. associated(x%coords%lcv%delp) ) then + __raise__(MAPL_RC_ERROR,"x%coords%lcv must have delp set") + end if + if ( Mie%nch /= 1 ) then + __raise__(MAPL_RC_ERROR,"for now, Mie must have a single channel") + end if + n = 1 + +! Loop over aerosol species +! ------------------------- + if ( iTau>0 ) y%r3(iTau)%q = 0.0 + if ( iExt>0 ) y%r3(iExt)%q = 0.0 + if ( iBck>0 ) y%r3(iBck)%q = 0.0 + + do iq = 1, x%n3d + + idxTable = Chem_MieQueryIdx(Mie,x%r3(iq)%name,rc) + if(idxTable == -1) cycle + if ( rc/=0 ) then + __raise__(MAPL_RC_ERROR,"cannot get Mie index for "//trim(x%r3(iq)%name)) + end if + + if (verbose_) & + print *, '[+] Adding '//trim(x%r3(iq)%name)//' contribution' + +! Loop over x, y, z +! ----------------- + do k = 1, km + do j = 1, jm + do i = 1, im + + delm = x%r3(iq)%q(i,j,k) * x%coords%lcv%delp(i,j,k) / MAPL_GRAV + delc = x%r3(iq)%q(i,j,k) * x%r3(iRHO)%q(i,j,k) ! concentration + delz = x%coords%lcv%delp(i,j,k) / (MAPL_GRAV*x%r3(iRHO)%q(i,j,k)) + rh = x%r3(iRH)%q(i,j,k) + + call Chem_MieQuery(Mie, idxTable, float(n), delm, rh, tau=tau, bbck=bck ) + + if ( iTau>0 ) y%r3(iTau)%q(i,j,k) = y%r3(iTau)%q(i,j,k) + tau + + if ( iExt>0 ) y%r3(iExt)%q(i,j,k) = y%r3(iExt)%q(i,j,k) + tau / delz + + if ( iBck>0 ) y%r3(iBck)%q(i,j,k) = y%r3(iBck)%q(i,j,k) + bck * delc + + end do ! longitudes + end do ! latitudes + end do ! levels + + end do ! aerosol tracers + + + if (verbose_) & + print *, '[x] All done!' + + rc = 0 + + end subroutine Chem_ExtCalculator + + end module Chem_AodMod + diff --git a/ESMF/Apps/aop_calculator.F90 b/ESMF/Apps/aop_calculator.F90 new file mode 100644 index 00000000..34e715bb --- /dev/null +++ b/ESMF/Apps/aop_calculator.F90 @@ -0,0 +1,188 @@ +# include "MAPL_Generic.h" + + program ext_calculator + +!----------------------------------------------------------------------- +! NASA/GSFC, Global Modeling and Assimilation Office, Code 610.1 ! +!----------------------------------------------------------------------- +!BOP +! +! !ROUTINE: ext_calculator --- Extinction calculator +! +! !INTERFACE: +! +! Usage: ext_calculator.x +! +! !USES: +! + use ESMF + use MAPL_Mod + use Chem_SimpleBundleMod + use Chem_RegistryMod + use Chem_MieMod + use Chem_AodMod + + use m_die + + implicit NONE + +! !DESCRIPTION: This is a parallel version of the 3D AOD Calculator. +! +! !REVISION HISTORY: +! +! 18Jun2011 da Silva Derived from mpana_aod.x +! +!EOP +!----------------------------------------------------------------------- + + character(len=*), parameter :: myname = 'ext_calculator' + +! Local variables +! --------------- + integer :: rc + integer :: nymd=0, nhms=0 + integer :: yy, mm, dd, h, m, s + + logical :: verbose = .TRUE. + + character(len=ESMF_MAXSTR) :: aer_registry, ext_registry + +! Control variables and obervations +! --------------------------------- + type (MAPL_SimpleBundle) :: q_f ! aerosol mixing ratio + type (MAPL_SimpleBundle) :: y_f ! extinction parameters + + type (Chem_Mie) :: Mie ! Mie Tables, etc + type (Chem_Registry) :: aerReg ! Registry with many species + type (Chem_Registry) :: extReg ! Registry with XX tracers + +! Basic ESMF objects +! ------------------ + type(ESMF_Config) :: CF ! Resource file + type(ESMF_Grid) :: etaGrid ! Eta Grid (lon, lat, eta) + type(ESMF_Time) :: Time + type(ESMF_VM) :: VM + + integer :: Nx, Ny ! Layout + integer :: IM_World, JM_World, LM_WORLD ! Global Grid dimensions + + call Main() + +CONTAINS + +!............................................................................................... + Subroutine Main() + + __Iam__('ext_calculator') + +! Initialize the ESMF. For performance reasons, it is important +! to turn OFF ESMF's automatic logging feature +! ------------------------------------------------------------- + call ESMF_Initialize (LogKindFlag=ESMF_LOGKIND_NONE, VM=VM, __RC__) + call ESMF_CalendarSetDefault ( ESMF_CALKIND_GREGORIAN, __RC__ ) + + if ( MAPL_am_I_root() ) then + print * + print *, ' --------------------------------------' + print *, ' 3D Extinction Calculator' + print *, ' --------------------------------------' + print * + end if + +! ------------------- +! ESMF Grid, Etc +! ------------------- + +! Load resources +! -------------- + CF = ESMF_ConfigCreate(__RC__) + call ESMF_ConfigLoadFile(CF, fileName='ext.rc', __RC__) + +! World grid dimensions and layout +! -------------------------------- + call ESMF_ConfigGetAttribute(CF, IM_World, Label='IM_World:', __RC__ ) + call ESMF_ConfigGetAttribute(CF, JM_World, Label='JM_World:', __RC__ ) + call ESMF_ConfigGetAttribute(CF, LM_World, Label='LM_World:', __RC__ ) + call ESMF_ConfigGetAttribute(CF, Nx, Label='Layout_Nx:', __RC__ ) + call ESMF_ConfigGetAttribute(CF, Ny, Label='Layout_Ny:', __RC__ ) + call ESMF_ConfigGetAttribute(CF, verbose, Label='verbose:', __RC__ ) + +! Create global grids: +! ------------------- + etaGrid = MAPL_LatLonGridCreate (Name='etaGrid', & + Nx = Nx, Ny = Ny, & + IM_World = IM_World, & + JM_World = JM_World, & + LM_World = LM_World, & + __RC__ ) + +! Validate grid +! ------------- + call ESMF_GridValidate(etaGrid,__RC__) + +! Get date/time from CF +! --------------------- + call ESMF_ConfigGetAttribute(CF, nymd, Label='nymd:', __RC__ ) + call ESMF_ConfigGetAttribute(CF, nhms, Label='nhms:', __RC__ ) + +! Create ESMF Time +! ---------------- + yy = nymd/10000; mm = (nymd-yy*10000) / 100; dd = nymd - (10000*yy + mm*100) + h = nhms/10000; m = (nhms - h*10000) / 100; s = nhms - (10000*h + m*100) + call ESMF_TimeSet(Time, yy=yy, mm=mm, dd=dd, h=h, m=m, s=s) + +! ------------------- +! Gridded Background +! ------------------- + +! Registries +! ---------- + call ESMF_ConfigGetAttribute(CF, aer_registry, Label='aer_registry:', __RC__ ) + aerReg = Chem_RegistryCreate ( rc, aer_registry ) + if ( rc == 0 ) then + if ( MAPL_AM_I_ROOT() ) then + call Chem_RegistryPrint(aerReg) + end if + else + call die(myname,'could not read Chem Registry for INPUT species') + end if + + call ESMF_ConfigGetAttribute(CF, ext_registry, Label='ext_registry:', __RC__ ) + extReg = Chem_RegistryCreate ( rc, ext_registry ) + if ( rc == 0 ) then + if ( MAPL_AM_I_ROOT() ) then + call Chem_RegistryPrint(extReg) + end if + else + call die(myname,'could not read Chem Registry for OUTPUT extinction parameters') + end if + +! Read aerosol mixing ratio +! ------------------------- + q_f = Chem_SimpleBundleRead (CF, 'aer_filename', etaGrid, & + time=Time, verbose=verbose, __RC__ ) + +! Perform Mie calculation +! ----------------------- + y_f = Chem_SimpleBundleCreate ('ext', extReg, etaGrid, __RC__ ) + Mie = Chem_MieCreate(CF, chemReg=aerReg, __RC__) + + call Chem_ExtCalculator (y_f, q_f, Mie, verbose, __RC__) + + call MAPL_SimpleBundlePrint(q_f) + call MAPL_SimpleBundlePrint(y_f) + +! Write file with AOD/Extinction output +! ------------------------------------ + call Chem_SimpleBundleWrite (y_f, CF, 'ext_filename', Time, __RC__ ) + +! All done +! -------- + call ESMF_Finalize ( __RC__ ) + + end subroutine Main + + end program ext_calculator + + + diff --git a/ESMF/GOCART2G_GridComp/CA2G_GridComp/CA2G_GridCompMod.F90 b/ESMF/GOCART2G_GridComp/CA2G_GridComp/CA2G_GridCompMod.F90 index 42b18d9a..460123ad 100644 --- a/ESMF/GOCART2G_GridComp/CA2G_GridComp/CA2G_GridCompMod.F90 +++ b/ESMF/GOCART2G_GridComp/CA2G_GridComp/CA2G_GridCompMod.F90 @@ -11,7 +11,7 @@ module CA2G_GridCompMod ! !USES: use ESMF use MAPL - use Chem_MieTableMod2G + use GOCART2G_MieMod use Chem_AeroGeneric use iso_c_binding, only: c_loc, c_f_pointer, c_ptr @@ -365,8 +365,10 @@ subroutine Initialize (GC, import, export, clock, RC) integer :: HDT ! model timestep (secs) real, pointer, dimension(:,:,:) :: ple logical :: data_driven - integer :: NUM_BANDS logical :: bands_are_present + integer, allocatable, dimension(:) :: channels_ + integer :: nmom_ + character(len=ESMF_MAXSTR) :: file_ __Iam__('Initialize') @@ -516,48 +518,21 @@ subroutine Initialize (GC, import, export, clock, RC) ! Create Radiation Mie Table ! -------------------------- - call MAPL_GetResource (MAPL, NUM_BANDS, 'NUM_BANDS:', __RC__) - -! Get file names for the optical tables - call ESMF_ConfigGetAttribute (cfg, self%rad_MieTable(instance)%optics_file, & - label="aerosol_radBands_optics_file:", __RC__ ) - - allocate (self%rad_MieTable(instance)%channels(NUM_BANDS), __STAT__ ) - self%rad_MieTable(instance)%nch = NUM_BANDS - - call ESMF_ConfigFindLabel(cfg, label="BANDS:", isPresent=bands_are_present, __RC__) - - if (bands_are_present) then - call ESMF_ConfigGetAttribute (cfg, self%rad_MieTable(instance)%channels, label= "BANDS:", & - count=self%rad_MieTable(instance)%nch, __RC__) - else - do i = 1, NUM_BANDS - self%rad_MieTable(instance)%channels(i) = i - end do - endif - - allocate (self%rad_MieTable(instance)%mie_aerosol, __STAT__) - self%rad_MieTable(instance)%mie_aerosol = Chem_MieTableCreate (self%rad_MieTable(instance)%optics_file, __RC__) - call Chem_MieTableRead (self%rad_MieTable(instance)%mie_aerosol, NUM_BANDS, self%rad_MieTable(instance)%channels, __RC__) + call ESMF_ConfigGetAttribute (cfg, file_, label="aerosol_radBands_optics_file:", __RC__ ) + self%rad_Mie = GOCART2G_Mie(trim(file_), __RC__) ! Create Diagnostics Mie Table ! ----------------------------- ! Get file names for the optical tables - call ESMF_ConfigGetAttribute (cfg, self%diag_MieTable(instance)%optics_file, & + call ESMF_ConfigGetAttribute (cfg, file_, & label="aerosol_monochromatic_optics_file:", __RC__ ) - call ESMF_ConfigGetAttribute (cfg, self%diag_MieTable(instance)%nmom, label="n_moments:", default=0, __RC__) - + call ESMF_ConfigGetAttribute (cfg, nmom_, label="n_moments:", default=0, __RC__) i = ESMF_ConfigGetLen (universal_cfg, label='aerosol_monochromatic_optics_wavelength_in_nm_from_LUT:', __RC__) - self%diag_MieTable(instance)%nch = i - allocate (self%diag_MieTable(instance)%channels(self%diag_MieTable(instance)%nch), __STAT__ ) - call ESMF_ConfigGetAttribute (universal_cfg, self%diag_MieTable(instance)%channels, & + allocate (channels_(i), __STAT__ ) + call ESMF_ConfigGetAttribute (universal_cfg, channels_, & label= "aerosol_monochromatic_optics_wavelength_in_nm_from_LUT:", __RC__) - - allocate (self%diag_MieTable(instance)%mie_aerosol, __STAT__) - self%diag_MieTable(instance)%mie_aerosol = Chem_MieTableCreate (self%diag_MieTable(instance)%optics_file, __RC__ ) - call Chem_MieTableRead (self%diag_MieTable(instance)%mie_aerosol, self%diag_MieTable(instance)%nch, & - self%diag_MieTable(instance)%channels*1.e-9, rc=status, nmom=self%diag_MieTable(instance)%nmom) - VERIFY_(status) + self%diag_Mie = GOCART2G_Mie(trim(file_), channels_*1.e-9, nmom=nmom_, __RC__) + deallocate(channels_) ! Finish creating AERO state ! -------------------------- @@ -836,14 +811,14 @@ subroutine Run1 (GC, import, export, clock, RC) nhms, self%cdt) end if - call CAEmission (self%diag_MieTable(self%instance), self%km, self%nbins, self%cdt, & + call CAEmission (self%diag_Mie, self%km, self%nbins, self%cdt, & MAPL_GRAV, GCsuffix, self%ratPOM, & self%eAircraftfuel, aircraft_fuel_src, & aviation_lto_src, aviation_cds_src, & aviation_crs_src, self%fHydrophobic, zpbl, t, airdens, rh2, & intPtr_philic, intPtr_phobic, delp, self%aviation_layers, biomass_src, & biogvoc_src, eocant1_src, eocant2_src, oc_ship_src, biofuel_src, & - CAEM, CAEMAN, CAEMBB, CAEMBF, CAEMBG, __RC__ ) + EM, EMAN, EMBB, EMBF, EMBG, __RC__ ) ! Read any pointwise emissions, if requested ! ------------------------------------------ @@ -974,8 +949,8 @@ subroutine Run2 (GC, import, export, clock, RC) where (1.01 * pSOA_VOC > MAPL_UNDEF) pSOA_VOC = 0.0 intPtr_philic = intPtr_philic + self%cdt * pSOA_VOC/airdens - if (associated(CAPSOA)) & - CAPSOA = sum(pSOA_VOC*delp/airdens/MAPL_GRAV, 3) + if (associated(PSOA)) & + PSOA = sum(pSOA_VOC*delp/airdens/MAPL_GRAV, 3) end if if (trim(comp_name) == 'CA.br') then @@ -983,14 +958,14 @@ subroutine Run2 (GC, import, export, clock, RC) where (1.01 * pSOA_VOC > MAPL_UNDEF) pSOA_VOC = 0.0 intPtr_philic = intPtr_philic + self%cdt * pSOA_VOC/airdens - if (associated(CAPSOA)) & - CAPSOA = sum(pSOA_VOC*delp/airdens/MAPL_GRAV, 3) + if (associated(PSOA)) & + PSOA = sum(pSOA_VOC*delp/airdens/MAPL_GRAV, 3) end if ! Ad Hoc transfer of hydrophobic to hydrophilic aerosols ! Following Chin's parameterization, the rate constant is ! k = 4.63e-6 s-1 (.4 day-1; e-folding time = 2.5 days) - call phobicTophilic (intPtr_phobic, intPtr_philic, CAHYPHIL, self%km, self%cdt, MAPL_GRAV, delp, __RC__) + call phobicTophilic (intPtr_phobic, intPtr_philic, HYPHIL, self%km, self%cdt, MAPL_GRAV, delp, __RC__) ! CA Settling ! ----------- @@ -1000,7 +975,7 @@ subroutine Run2 (GC, import, export, clock, RC) call Chem_Settling (self%km, self%klid, n, self%rhFlag, self%cdt, MAPL_GRAV, & self%radius(n)*1.e-6, self%rhop(n), int_ptr, t, airdens, & - rh2, zle, delp, CASD, __RC__) + rh2, zle, delp, SD, __RC__) end do ! CA Deposition @@ -1018,21 +993,21 @@ subroutine Run2 (GC, import, export, clock, RC) dqa = 0. dqa = max(0.0, int_ptr(:,:,self%km)*(1.-exp(-drydepositionfrequency*self%cdt))) int_ptr(:,:,self%km) = int_ptr(:,:,self%km) - dqa - if (associated(CADP)) then - CADP(:,:,n) = dqa * delp(:,:,self%km) / MAPL_GRAV / self%cdt + if (associated(DP)) then + DP(:,:,n) = dqa * delp(:,:,self%km) / MAPL_GRAV / self%cdt end if end do ! Large-scale Wet Removal ! ------------------------------- ! Hydrophobic mode (first tracer) is not removed - if (associated(CAWT)) CAWT(:,:,1)=0.0 + if (associated(WT)) WT(:,:,1)=0.0 KIN = .true. ! Hydrophilic mode (second tracer) is removed fwet = 1. call WetRemovalGOCART2G (self%km, self%klid, self%nbins, self%nbins, 2, self%cdt, GCsuffix, & - KIN, MAPL_GRAV, fwet, CAphilic, ple, t, airdens, & - pfl_lsan, pfi_lsan, cn_prcp, ncn_prcp, CAWT, __RC__) + KIN, MAPL_GRAV, fwet, philic, ple, t, airdens, & + pfl_lsan, pfi_lsan, cn_prcp, ncn_prcp, WT, __RC__) ! Compute diagnostics ! ------------------- @@ -1040,15 +1015,14 @@ subroutine Run2 (GC, import, export, clock, RC) int_arr(:,:,:,1) = intPtr_phobic int_arr(:,:,:,2) = intPtr_philic - call Aero_Compute_Diags (mie_table=self%diag_MieTable(self%instance), km=self%km, klid=self%klid, nbegin=1, nbins=2, & - channels=self%diag_MieTable(self%instance)%channels*1.0e-9, & + call Aero_Compute_Diags (mie=self%diag_Mie, km=self%km, klid=self%klid, nbegin=1, nbins=2, & wavelengths_profile=self%wavelengths_profile*1.0e-9, & wavelengths_vertint=self%wavelengths_vertint*1.0e-9, aerosol=int_arr, grav=MAPL_GRAV, & tmpu=t, rhoa=airdens, rh=rh2, u=u, v=v, delp=delp, ple=ple, tropp=tropp, & - sfcmass=CASMASS, colmass=CACMASS, mass=CAMASS,& - exttau=CAEXTTAU,stexttau=CASTEXTTAU, scatau=CASCATAU, stscatau=CASTSCATAU,& - fluxu=CAFLUXU, fluxv=CAFLUXV, & - conc=CACONC, extcoef=CAEXTCOEF, scacoef=CASCACOEF, angstrom=CAANGSTR, aerindx=CAAERIDX,& + sfcmass=SMASS, colmass=CMASS, mass=MASS,& + exttau=EXTTAU,stexttau=STEXTTAU, scatau=SCATAU, stscatau=STSCATAU,& + fluxu=FLUXU, fluxv=FLUXV, & + conc=CONC, extcoef=EXTCOEF, scacoef=SCACOEF, angstrom=ANGSTR, aerindx=AERIDX,& NO3nFlag=.false., __RC__) @@ -1153,8 +1127,7 @@ subroutine aerosol_optics(state, rc) integer :: instance integer :: n, nbins integer :: i1, j1, i2, j2, km - integer :: band, offset - integer, parameter :: n_bands = 1 + integer :: band integer :: i, j, k @@ -1176,7 +1149,6 @@ subroutine aerosol_optics(state, rc) ! -------------- band = 0 call ESMF_AttributeGet(state, name='band_for_aerosol_optics', value=band, __RC__) - offset = band - n_bands ! Pressure at layer edges ! ------------------------ @@ -1223,7 +1195,7 @@ subroutine aerosol_optics(state, rc) address = transfer(opaque_self, address) call c_f_pointer(address, self) - call mie_ (self%rad_MieTable(instance), nbins, n_bands, offset, q_4d, rh, ext_s, ssa_s, asy_s, __RC__) + call mie_ (self%rad_Mie, nbins, band, q_4d, rh, ext_s, ssa_s, asy_s, __RC__) call ESMF_AttributeGet(state, name='extinction_in_air_due_to_ambient_aerosol', value=fld_name, __RC__) if (fld_name /= '') then @@ -1250,14 +1222,13 @@ subroutine aerosol_optics(state, rc) contains - subroutine mie_(mie_table, nbins, nb, offset, q, rh, bext_s, bssa_s, basym_s, rc) + subroutine mie_(mie, nbins, band, q, rh, bext_s, bssa_s, basym_s, rc) implicit none - type(Chem_Mie), intent(inout) :: mie_table ! mie table + type(GOCART2G_Mie), intent(inout) :: mie ! mie table integer, intent(in ) :: nbins ! number of bins - integer, intent(in ) :: nb ! number of bands - integer, intent(in ) :: offset ! bands offset + integer, intent(in ) :: band ! channel real, intent(in ) :: q(:,:,:,:) ! aerosol mass mixing ratio, kg kg-1 real, intent(in ) :: rh(:,:,:) ! relative humidity real(kind=8), intent( out) :: bext_s (size(ext_s,1),size(ext_s,2),size(ext_s,3)) @@ -1277,7 +1248,8 @@ subroutine mie_(mie_table, nbins, nb, offset, q, rh, bext_s, bssa_s, basym_s, rc basym_s = 0.0d0 do l = 1, nbins - call Chem_MieQuery(mie_table, l, real(offset+1.), q(:,:,:,l), rh, tau=bext, gasym=gasym, ssa=bssa) + !tau is converted to bext + call mie%Query( band, l, q(:,:,:,l), rh, tau=bext, gasym=gasym, ssa=bssa, __RC__) bext_s = bext_s + bext ! extinction bssa_s = bssa_s + (bssa*bext) ! scattering extinction @@ -1318,7 +1290,7 @@ subroutine monochromatic_aerosol_optics(state, rc) integer :: instance integer :: n, nbins integer :: i1, j1, i2, j2, km - real :: wavelength, mieTable_index + real :: wavelength integer :: i, j, k __Iam__('CA2G::monochromatic_aerosol_optics') @@ -1339,21 +1311,6 @@ subroutine monochromatic_aerosol_optics(state, rc) ! -------------------- call ESMF_AttributeGet(state, name='wavelength_for_aerosol_optics', value=wavelength, __RC__) -! Get wavelength index for Mie Table -! ---------------------------------- -! Channel values are 4.7e-7 5.5e-7 6.7e-7 8.7e-7 [meter]. Their indices are 1,2,3,4 respectively. - if ((wavelength .ge. 4.69e-7) .and. (wavelength .le. 4.71e-7)) then - mieTable_index = 1. - else if ((wavelength .ge. 5.49e-7) .and. (wavelength .le. 5.51e-7)) then - mieTable_index = 2. - else if ((wavelength .ge. 6.69e-7) .and. (wavelength .le. 6.71e-7)) then - mieTable_index = 3. - else if ((wavelength .ge. 8.68e-7) .and. (wavelength .le. 8.71e-7)) then - mieTable_index = 4. - else - print*,trim(Iam),' : wavelengths of ',wavelength,' is an invalid value.' - return - end if ! Pressure at layer edges ! ------------------------ @@ -1400,14 +1357,8 @@ subroutine monochromatic_aerosol_optics(state, rc) call c_f_pointer(address, self) do n = 1, nbins - do i = 1, i2 - do j = 1, j2 - do k = 1, km - call Chem_MieQuery(self%diag_MieTable(instance), n, mieTable_index, q_4d(i,j,k,n), rh(i,j,k), tau=tau(i,j,k), __RC__) - tau_s(i,j,k) = tau_s(i,j,k) + tau(i,j,k) - end do - end do - end do + call self%diag_Mie%Query(wavelength, n, q_4d(:,:,:,n), rh, tau=tau, __RC__) + tau_s = tau_s + tau end do call ESMF_AttributeGet(state, name='monochromatic_extinction_in_air_due_to_ambient_aerosol', value=fld_name, __RC__) diff --git a/ESMF/GOCART2G_GridComp/CA2G_GridComp/CA2G_StateSpecs.rc b/ESMF/GOCART2G_GridComp/CA2G_GridComp/CA2G_StateSpecs.rc index 00af8374..bb3a6cd0 100644 --- a/ESMF/GOCART2G_GridComp/CA2G_GridComp/CA2G_StateSpecs.rc +++ b/ESMF/GOCART2G_GridComp/CA2G_GridComp/CA2G_StateSpecs.rc @@ -74,32 +74,32 @@ category: EXPORT #---------------------------------------------------------------------------------------- NAME | UNITS | DIMS| VLOC| UNGRIDDED | LONG NAME #---------------------------------------------------------------------------------------- - CAMASS* | kg kg-1 | xyz | C | | Carbonaceous Aerosol Mass Mixing Ratio - CACONC* | kg m-3 | xyz | C | | Carbonaceous Aerosol Mass Concentration - CAEXTCOEF* | m-1 | xyz | C | size(self%wavelengths_profile) | Carbonaceous Aerosol Extinction Coefficient - CASCACOEF* | m-1 | xyz | C | size(self%wavelengths_profile) | Carbonaceous Aerosol Scattering Coefficient + *MASS | kg kg-1 | xyz | C | | Carbonaceous Aerosol Mass Mixing Ratio + *CONC | kg m-3 | xyz | C | | Carbonaceous Aerosol Mass Concentration + *EXTCOEF | m-1 | xyz | C | size(self%wavelengths_profile) | Carbonaceous Aerosol Extinction Coefficient + *SCACOEF | m-1 | xyz | C | size(self%wavelengths_profile) | Carbonaceous Aerosol Scattering Coefficient #...........|............|.....|.....|.......|............................................ - CAEM* | kg m-2 s-1 | xy | N | nbins | Carbonaceous Aerosol Emission (Bin %d) - CASD* | kg m-2 s-1 | xy | N | nbins | Carbonaceous Aerosol Sedimentation (Bin %d) - CADP* | kg m-2 s-1 | xy | N | nbins | Carbonaceous Aerosol Dry Deposition (Bin %d) - CAWT* | kg m-2 s-1 | xy | N | nbins | Carbonaceous Aerosol Wet Deposition (Bin %d) - CASV* | kg m-2 s-1 | xy | N | nbins | Carbonaceous Aerosol Convective Scavenging (Bin %d) - CAEMAN* | kg m-2 s-1 | xy | N | | Carbonaceous Aerosol Anthropogenic Emissions - CAEMBB* | kg m-2 s-1 | xy | N | | Carbonaceous Aerosol Biomass Burning Emissions - CAEMBF* | kg m-2 s-1 | xy | N | | Carbonaceous Aerosol Biofuel Emissions - CAEMBG* | kg m-2 s-1 | xy | N | | Carbonaceous Aerosol Biogenic Emissions - CAHYPHIL* | kg m-2 s-1 | xy | N | | Carbonaceous Aerosol Hydrophobic to Hydrophilic - CAPSOA* | kg m-2 s-1 | xy | N | | Carbonaceous Aerosol SOA Production - CASMASS* | kg m-3 | xy | N | | Carbonaceous Aerosol Surface Mass Concentration - CACMASS* | kg m-2 | xy | N | | Carbonaceous Aerosol Column Mass Density - CAEXTTAU* | 1 | xy | N | size(self%wavelengths_vertint) | Carbonaceous Aerosol Extinction AOT - CASTEXTTAU*| 1 | xy | N | size(self%wavelengths_vertint) | Carbonaceous Aerosol Extinction AOT Stratosphere - CASCATAU* | 1 | xy | N | size(self%wavelengths_vertint) | Carbonaceous Aerosol Scattering AOT - CASTSCATAU*| 1 | xy | N | size(self%wavelengths_vertint) | Carbonaceous Aerosol Scattering AOT Stratosphere - CAANGSTR* | 1 | xy | N | | Carbonaceous Aerosol Angstrom parameter [470-870 nm] - CAFLUXU* | kg m-1 s-1 | xy | N | | Carbonaceous Aerosol column u-wind mass flux - CAFLUXV* | kg m-1 s-1 | xy | N | | Carbonaceous Aerosol column v-wind mass flux - CAAERIDX* | 1 | xy | N | | Carbonaceous Aerosol TOMS UV Aerosol Index + *EM | kg m-2 s-1 | xy | N | nbins | Carbonaceous Aerosol Emission (Bin %d) + *SD | kg m-2 s-1 | xy | N | nbins | Carbonaceous Aerosol Sedimentation (Bin %d) + *DP | kg m-2 s-1 | xy | N | nbins | Carbonaceous Aerosol Dry Deposition (Bin %d) + *WT | kg m-2 s-1 | xy | N | nbins | Carbonaceous Aerosol Wet Deposition (Bin %d) + *SV | kg m-2 s-1 | xy | N | nbins | Carbonaceous Aerosol Convective Scavenging (Bin %d) + *EMAN | kg m-2 s-1 | xy | N | | Carbonaceous Aerosol Anthropogenic Emissions + *EMBB | kg m-2 s-1 | xy | N | | Carbonaceous Aerosol Biomass Burning Emissions + *EMBF | kg m-2 s-1 | xy | N | | Carbonaceous Aerosol Biofuel Emissions + *EMBG | kg m-2 s-1 | xy | N | | Carbonaceous Aerosol Biogenic Emissions + *HYPHIL | kg m-2 s-1 | xy | N | | Carbonaceous Aerosol Hydrophobic to Hydrophilic + *PSOA | kg m-2 s-1 | xy | N | | Carbonaceous Aerosol SOA Production + *SMASS | kg m-3 | xy | N | | Carbonaceous Aerosol Surface Mass Concentration + *CMASS | kg m-2 | xy | N | | Carbonaceous Aerosol Column Mass Density + *EXTTAU | 1 | xy | N | size(self%wavelengths_vertint) | Carbonaceous Aerosol Extinction AOT + *STEXTTAU| 1 | xy | N | size(self%wavelengths_vertint) | Carbonaceous Aerosol Extinction AOT Stratosphere + *SCATAU | 1 | xy | N | size(self%wavelengths_vertint) | Carbonaceous Aerosol Scattering AOT + *STSCATAU| 1 | xy | N | size(self%wavelengths_vertint) | Carbonaceous Aerosol Scattering AOT Stratosphere + *ANGSTR | 1 | xy | N | | Carbonaceous Aerosol Angstrom parameter [470-870 nm] + *FLUXU | kg m-1 s-1 | xy | N | | Carbonaceous Aerosol column u-wind mass flux + *FLUXV | kg m-1 s-1 | xy | N | | Carbonaceous Aerosol column v-wind mass flux + *AERIDX | 1 | xy | N | | Carbonaceous Aerosol TOMS UV Aerosol Index @@ -109,8 +109,8 @@ category: INTERNAL #---------------------------------------------------------------------------------------- NAME | UNITS | DIMS | VLOC | RESTART | ADD2EXPORT | FRIENDLYTO | LONG NAME #---------------------------------------------------------------------------------------- - CAphobic* |kg kg-1| xyz | C | MAPL_RestartOptional | T | DYNAMICS:TURBULENCE:MOIST | Carbonaceous Aerosol Mixing Ratio - CAphilic* |kg kg-1| xyz | C | MAPL_RestartOptional | T | DYNAMICS:TURBULENCE:MOIST | Carbonaceous Aerosol Mixing Ratio + *phobic |kg kg-1| xyz | C | MAPL_RestartOptional | T | DYNAMICS:TURBULENCE:MOIST | Carbonaceous Aerosol Mixing Ratio + *philic |kg kg-1| xyz | C | MAPL_RestartOptional | T | DYNAMICS:TURBULENCE:MOIST | Carbonaceous Aerosol Mixing Ratio #******************************************************** # diff --git a/ESMF/GOCART2G_GridComp/DU2G_GridComp/DU2G_GridCompMod.F90 b/ESMF/GOCART2G_GridComp/DU2G_GridComp/DU2G_GridCompMod.F90 index 83d17cff..150ff746 100644 --- a/ESMF/GOCART2G_GridComp/DU2G_GridComp/DU2G_GridCompMod.F90 +++ b/ESMF/GOCART2G_GridComp/DU2G_GridComp/DU2G_GridCompMod.F90 @@ -11,7 +11,7 @@ module DU2G_GridCompMod ! !USES: use ESMF use MAPL - use Chem_MieTableMod2G + use GOCART2G_MieMod use Chem_AeroGeneric use iso_c_binding, only: c_loc, c_f_pointer, c_ptr @@ -373,7 +373,9 @@ subroutine Initialize (GC, import, export, clock, RC) logical :: data_driven integer :: NUM_BANDS logical :: bands_are_present - + integer, allocatable, dimension(:) :: channels_ + integer :: nmom_ + character(len=ESMF_MAXSTR) :: file_ __Iam__('Initialize') !**************************************************************************** @@ -516,49 +518,21 @@ subroutine Initialize (GC, import, export, clock, RC) ! Create Radiation Mie Table ! -------------------------- - call MAPL_GetResource (MAPL, NUM_BANDS, 'NUM_BANDS:', __RC__) - -! Get file names for the optical tables - call ESMF_ConfigGetAttribute (cfg, self%rad_MieTable(instance)%optics_file, & - label="aerosol_radBands_optics_file:", __RC__ ) - - allocate (self%rad_MieTable(instance)%channels(NUM_BANDS), __STAT__ ) - self%rad_MieTable(instance)%nch = NUM_BANDS - - call ESMF_ConfigFindLabel(cfg, label="BANDS:", isPresent=bands_are_present, __RC__) - - if (bands_are_present) then - call ESMF_ConfigGetAttribute (cfg, self%rad_MieTable(instance)%channels, label= "BANDS:", & - count=self%rad_MieTable(instance)%nch, __RC__) - else - do i = 1, NUM_BANDS - self%rad_MieTable(instance)%channels(i) = i - end do - endif - - allocate (self%rad_MieTable(instance)%mie_aerosol, __STAT__) - self%rad_MieTable(instance)%mie_aerosol = Chem_MieTableCreate (self%rad_MieTable(instance)%optics_file, __RC__) - call Chem_MieTableRead (self%rad_MieTable(instance)%mie_aerosol, NUM_BANDS, self%rad_MieTable(instance)%channels, __RC__) + call ESMF_ConfigGetAttribute (cfg, file_, label="aerosol_radBands_optics_file:", __RC__ ) + self%rad_Mie = GOCART2G_Mie(trim(file_), __RC__) ! Create Diagnostics Mie Table ! ----------------------------- ! Get file names for the optical tables - call ESMF_ConfigGetAttribute (cfg, self%diag_MieTable(instance)%optics_file, & + call ESMF_ConfigGetAttribute (cfg, file_, & label="aerosol_monochromatic_optics_file:", __RC__ ) - call ESMF_ConfigGetAttribute (cfg, self%diag_MieTable(instance)%nmom, label="n_moments:", default=0, __RC__) - + call ESMF_ConfigGetAttribute (cfg, nmom_, label="n_moments:", default=0, __RC__) i = ESMF_ConfigGetLen (universal_cfg, label='aerosol_monochromatic_optics_wavelength_in_nm_from_LUT:', __RC__) - self%diag_MieTable(instance)%nch = i - allocate (self%diag_MieTable(instance)%channels(self%diag_MieTable(instance)%nch), __STAT__ ) - call ESMF_ConfigGetAttribute (universal_cfg, self%diag_MieTable(instance)%channels, & - label= "aerosol_monochromatic_optics_wavelength_in_nm_from_LUT:", __RC__) - - allocate (self%diag_MieTable(instance)%mie_aerosol, __STAT__) - self%diag_MieTable(instance)%mie_aerosol = Chem_MieTableCreate (self%diag_MieTable(instance)%optics_file, __RC__ ) - call Chem_MieTableRead (self%diag_MieTable(instance)%mie_aerosol, self%diag_MieTable(instance)%nch, & - self%diag_MieTable(instance)%channels*1.e-9, rc=status, nmom=self%diag_MieTable(instance)%nmom) - VERIFY_(status) - + allocate (channels_(i), __STAT__ ) + call ESMF_ConfigGetAttribute (universal_cfg, channels_, & + label= "aerosol_monochromatic_optics_wavelength_in_nm_from_LUT:", __RC__) + self%diag_Mie = GOCART2G_Mie(trim(file_), channels_*1.e-9, nmom=nmom_, __RC__) + deallocate(channels_) ! Mie Table instance/index call ESMF_AttributeSet (aero, name='mie_table_instance', value=instance, __RC__) @@ -978,8 +952,8 @@ subroutine Run2 (GC, import, export, clock, RC) ! Compute diagnostics ! ------------------- ! Certain variables are multiplied by 1.0e-9 to convert from nanometers to meters - call Aero_Compute_Diags (self%diag_MieTable(self%instance), self%km, self%klid, 1, self%nbins, self%rlow, & - self%rup, self%diag_MieTable(self%instance)%channels*1.0e-9, self%wavelengths_profile*1.0e-9, & + call Aero_Compute_Diags (self%diag_Mie, self%km, self%klid, 1, self%nbins, self%rlow, & + self%rup, self%wavelengths_profile*1.0e-9, & self%wavelengths_vertint*1.0e-9, DU, MAPL_GRAV, t, airdens, & rh2, u, v, delp, ple,tropp, & DUSMASS, DUCMASS, DUMASS, DUEXTTAU, DUSTEXTTAU, DUSCATAU,DUSTSCATAU, & @@ -1076,8 +1050,7 @@ subroutine aerosol_optics(state, rc) integer :: instance integer :: n, nbins integer :: i1, j1, i2, j2, km - integer :: band, offset - integer, parameter :: n_bands = 1 + integer :: band integer :: k @@ -1093,7 +1066,6 @@ subroutine aerosol_optics(state, rc) ! -------------- band = 0 call ESMF_AttributeGet (state, name='band_for_aerosol_optics', value=band, __RC__) - offset = band - n_bands ! Pressure at layer edges ! ------------------------ @@ -1139,7 +1111,7 @@ subroutine aerosol_optics(state, rc) address = transfer(opaque_self, address) call c_f_pointer(address, self) - call mie_ (self%rad_MieTable(instance), nbins, n_bands, offset, q_4d, rh, ext_s, ssa_s, asy_s, __RC__) + call mie_ (self%rad_Mie, nbins, band, q_4d, rh, ext_s, ssa_s, asy_s, __RC__) call ESMF_AttributeGet (state, name='extinction_in_air_due_to_ambient_aerosol', value=fld_name, __RC__) if (fld_name /= '') then @@ -1166,14 +1138,13 @@ subroutine aerosol_optics(state, rc) contains - subroutine mie_(mie_table, nbins, nb, offset, q, rh, bext_s, bssa_s, basym_s, rc) + subroutine mie_(mie, nbins, band, q, rh, bext_s, bssa_s, basym_s, rc) implicit none - - type(Chem_Mie), intent(inout) :: mie_table ! mie table + + type(GOCART2G_Mie), intent(inout) :: mie ! mie table integer, intent(in ) :: nbins ! number of bins - integer, intent(in ) :: nb ! number of bands - integer, intent(in ) :: offset ! bands offset + integer, intent(in ) :: band ! channel real, intent(in ) :: q(:,:,:,:) ! aerosol mass mixing ratio, kg kg-1 real, intent(in ) :: rh(:,:,:) ! relative humidity real(kind=DP), intent( out) :: bext_s (size(ext_s,1),size(ext_s,2),size(ext_s,3)) @@ -1194,8 +1165,8 @@ subroutine mie_(mie_table, nbins, nb, offset, q, rh, bext_s, bssa_s, basym_s, rc basym_s = 0.0d0 do l = 1, nbins - call Chem_MieQuery(mie_table, l, real(offset+1.), q(:,:,:,l), rh, tau=bext, gasym=gasym, ssa=bssa) - + ! tau is converted to bext + call mie%Query(band, l, q(:,:,:,l), rh, tau=bext, gasym=gasym, ssa=bssa, __RC__) bext_s = bext_s + bext ! extinction bssa_s = bssa_s + (bssa*bext) ! scattering extinction basym_s = basym_s + gasym*(bssa*bext) ! asymetry parameter multiplied by scatering extiction @@ -1231,7 +1202,7 @@ subroutine monochromatic_aerosol_optics(state, rc) integer :: instance integer :: n, nbins, k integer :: i1, j1, i2, j2, km, i, j - real :: wavelength, mieTable_index + real :: wavelength __Iam__('DU2G::monochromatic_aerosol_optics') @@ -1246,21 +1217,6 @@ subroutine monochromatic_aerosol_optics(state, rc) wavelength = 0. call ESMF_AttributeGet (state, name='wavelength_for_aerosol_optics', value=wavelength, __RC__) -! Get wavelength index for Mie Table -! Channel values are 4.7e-7 5.5e-7 6.7e-7 8.7e-7 [meter]. Their indices are 1,2,3,4 respectively. - if ((wavelength .ge. 4.69e-7) .and. (wavelength .le. 4.71e-7)) then - mieTable_index = 1. - else if ((wavelength .ge. 5.49e-7) .and. (wavelength .le. 5.51e-7)) then - mieTable_index = 2. - else if ((wavelength .ge. 6.69e-7) .and. (wavelength .le. 6.71e-7)) then - mieTable_index = 3. - else if ((wavelength .ge. 8.68e-7) .and. (wavelength .le. 8.71e-7)) then - mieTable_index = 4. - else - print*,trim(Iam),' : wavelength of ',wavelength,' is an invalid value.' - return - end if - ! Pressure at layer edges ! ------------------------ call ESMF_AttributeGet (state, name='air_pressure_for_aerosol_optics', value=fld_name, __RC__) @@ -1308,14 +1264,8 @@ subroutine monochromatic_aerosol_optics(state, rc) call c_f_pointer(address, self) do n = 1, nbins - do i = 1, i2 - do j = 1, j2 - do k = 1, km - call Chem_MieQuery(self%diag_MieTable(instance), n, mieTable_index, q_4d(i,j,k,n), rh(i,j,k), tau=tau(i,j,k), __RC__) - tau_s(i,j,k) = tau_s(i,j,k) + tau(i,j,k) - end do - end do - end do + call self%diag_Mie%Query(wavelength, n, q_4d(:,:,:,n), rh, tau=tau, __RC__) + tau_s = tau_s + tau end do call ESMF_AttributeGet (state, name='monochromatic_extinction_in_air_due_to_ambient_aerosol', value=fld_name, __RC__) diff --git a/ESMF/GOCART2G_GridComp/GA_Environment/GA_EnvironmentMod.F90 b/ESMF/GOCART2G_GridComp/GA_Environment/GA_EnvironmentMod.F90 index 487f493f..bfac0c34 100644 --- a/ESMF/GOCART2G_GridComp/GA_Environment/GA_EnvironmentMod.F90 +++ b/ESMF/GOCART2G_GridComp/GA_Environment/GA_EnvironmentMod.F90 @@ -4,7 +4,7 @@ module GA_EnvironmentMod use ESMF use MAPL - use Chem_MieTableMod2G + use GOCART2G_MieMod implicit none private @@ -12,7 +12,7 @@ module GA_EnvironmentMod public :: GA_Environment type :: GA_Environment - type(Chem_Mie), dimension(2) :: rad_MieTable, diag_MieTable + type(GOCART2G_Mie) :: rad_Mie, diag_Mie real, allocatable :: radius(:) ! particle effective radius [um] real, allocatable :: rhop(:) ! soil class density [kg m-3] real, allocatable :: fscav(:) ! scavenging efficiency diff --git a/ESMF/GOCART2G_GridComp/NI2G_GridComp/NI2G_GridCompMod.F90 b/ESMF/GOCART2G_GridComp/NI2G_GridComp/NI2G_GridCompMod.F90 index d9901b66..82fa5c14 100644 --- a/ESMF/GOCART2G_GridComp/NI2G_GridComp/NI2G_GridCompMod.F90 +++ b/ESMF/GOCART2G_GridComp/NI2G_GridComp/NI2G_GridCompMod.F90 @@ -11,7 +11,7 @@ module NI2G_GridCompMod ! !USES: use ESMF use MAPL - use Chem_MieTableMod2G + use GOCART2G_MieMod use Chem_AeroGeneric use iso_c_binding, only: c_loc, c_f_pointer, c_ptr @@ -310,7 +310,9 @@ subroutine Initialize (GC, IMPORT, EXPORT, CLOCK, RC) real, dimension(4) :: Vect_Hcts ! real, allocatable, dimension(:) :: rmedDU, rmedSS, fnumDU, fnumSS integer :: itemCount - + integer, allocatable, dimension(:) :: channels_ + integer :: nmom_ + character(len=ESMF_MAXSTR) :: file_ __Iam__('Initialize') !**************************************************************************** @@ -475,51 +477,24 @@ subroutine Initialize (GC, IMPORT, EXPORT, CLOCK, RC) end if self%instance = instance - + ! Create Radiation Mie Table ! -------------------------- - call MAPL_GetResource (MAPL, NUM_BANDS, 'NUM_BANDS:', __RC__) - -! Get file names for the optical tables - call ESMF_ConfigGetAttribute (cfg, self%rad_MieTable(instance)%optics_file, & - label="aerosol_radBands_optics_file:", __RC__ ) - - allocate (self%rad_MieTable(instance)%channels(NUM_BANDS), __STAT__ ) - self%rad_MieTable(instance)%nch = NUM_BANDS - - call ESMF_ConfigFindLabel(cfg, label="BANDS:", isPresent=bands_are_present, __RC__) - - if (bands_are_present) then - call ESMF_ConfigGetAttribute (cfg, self%rad_MieTable(instance)%channels, label= "BANDS:", & - count=self%rad_MieTable(instance)%nch, __RC__) - else - do i = 1, NUM_BANDS - self%rad_MieTable(instance)%channels(i) = i - end do - endif - - allocate (self%rad_MieTable(instance)%mie_aerosol, __STAT__) - self%rad_MieTable(instance)%mie_aerosol = Chem_MieTableCreate (self%rad_MieTable(instance)%optics_file, __RC__) - call Chem_MieTableRead (self%rad_MieTable(instance)%mie_aerosol, NUM_BANDS, self%rad_MieTable(instance)%channels, __RC__) + call ESMF_ConfigGetAttribute (cfg, file_, label="aerosol_radBands_optics_file:", __RC__ ) + self%rad_Mie = GOCART2G_Mie(trim(file_), __RC__) ! Create Diagnostics Mie Table ! ----------------------------- ! Get file names for the optical tables - call ESMF_ConfigGetAttribute (cfg, self%diag_MieTable(instance)%optics_file, & + call ESMF_ConfigGetAttribute (cfg, file_, & label="aerosol_monochromatic_optics_file:", __RC__ ) - call ESMF_ConfigGetAttribute (cfg, self%diag_MieTable(instance)%nmom, label="n_moments:", default=0, __RC__) - + call ESMF_ConfigGetAttribute (cfg, nmom_, label="n_moments:", default=0, __RC__) i = ESMF_ConfigGetLen (universal_cfg, label='aerosol_monochromatic_optics_wavelength_in_nm_from_LUT:', __RC__) - self%diag_MieTable(instance)%nch = i - allocate (self%diag_MieTable(instance)%channels(self%diag_MieTable(instance)%nch), __STAT__ ) - call ESMF_ConfigGetAttribute (universal_cfg, self%diag_MieTable(instance)%channels, & + allocate (channels_(i), __STAT__ ) + call ESMF_ConfigGetAttribute (universal_cfg, channels_, & label= "aerosol_monochromatic_optics_wavelength_in_nm_from_LUT:", __RC__) - - allocate (self%diag_MieTable(instance)%mie_aerosol, __STAT__) - self%diag_MieTable(instance)%mie_aerosol = Chem_MieTableCreate (self%diag_MieTable(instance)%optics_file, __RC__ ) - call Chem_MieTableRead (self%diag_MieTable(instance)%mie_aerosol, self%diag_MieTable(instance)%nch, & - self%diag_MieTable(instance)%channels*1.e-9, rc=status, nmom=self%diag_MieTable(instance)%nmom) - VERIFY_(status) + self%diag_Mie = GOCART2G_Mie(trim(file_), channels_*1.e-9, nmom=nmom_, __RC__) + deallocate(channels_) ! Mie Table instance/index call ESMF_AttributeSet(aero, name='mie_table_instance', value=instance, __RC__) @@ -955,8 +930,8 @@ subroutine Run2 (GC, import, export, clock, RC) allocate(aerosol(ubound(NH4a,1), ubound(NH4a,2), ubound(NH4a,3), 3), __STAT__) aerosol(:,:,:,:) = 0.0 aerosol(:,:,:,1) = NH4a - call Aero_Compute_Diags (mie_table=self%diag_MieTable(self%instance), km=self%km, klid=self%klid, nbegin=1, & - nbins=1, channels=self%diag_MieTable(self%instance)%channels*1.0e-9, & + call Aero_Compute_Diags (mie=self%diag_Mie, km=self%km, klid=self%klid, nbegin=1, & + nbins=1, & wavelengths_profile=self%wavelengths_profile*1.0e-9, & wavelengths_vertint=self%wavelengths_vertint*1.0e-9, & aerosol=aerosol, grav=MAPL_GRAV, tmpu=t, rhoa=airdens, rh=rh2, u=u, v=v, & @@ -964,8 +939,8 @@ subroutine Run2 (GC, import, export, clock, RC) sfcmass=NH4SMASS, colmass=NH4CMASS, mass=NH4MASS, conc=NH4CONC, __RC__) aerosol(:,:,:,1) = NH3 - call Aero_Compute_Diags (mie_table=self%diag_MieTable(self%instance), km=self%km, klid=self%klid, nbegin=1, & - nbins=1, channels=self%diag_MieTable(self%instance)%channels*1.0e-9, & + call Aero_Compute_Diags (mie=self%diag_Mie, km=self%km, klid=self%klid, nbegin=1, & + nbins=1, & wavelengths_profile=self%wavelengths_profile*1.0e-9, & wavelengths_vertint=self%wavelengths_vertint*1.0e-9, & aerosol=aerosol, grav=MAPL_GRAV, tmpu=t, rhoa=airdens, rh=rh2, u=u, v=v, & @@ -973,8 +948,8 @@ subroutine Run2 (GC, import, export, clock, RC) sfcmass=NH3SMASS, colmass=NH3CMASS, mass=NH3MASS, conc=NH3CONC, __RC__) aerosol(:,:,:,1) = NO3an1 - call Aero_Compute_Diags (mie_table=self%diag_MieTable(self%instance), km=self%km, klid=self%klid, nbegin=1, & - nbins=1, channels=self%diag_MieTable(self%instance)%channels*1.0e-9, & + call Aero_Compute_Diags (mie=self%diag_Mie, km=self%km, klid=self%klid, nbegin=1, & + nbins=1, & wavelengths_profile=self%wavelengths_profile*1.0e-9, & wavelengths_vertint=self%wavelengths_vertint*1.0e-9, & aerosol=aerosol, grav=MAPL_GRAV, tmpu=t, rhoa=airdens, rh=rh2, u=u, v=v, & @@ -986,8 +961,8 @@ subroutine Run2 (GC, import, export, clock, RC) aerosol(:,:,:,1) = NO3an1 aerosol(:,:,:,2) = NO3an2 aerosol(:,:,:,3) = NO3an3 - call Aero_Compute_Diags (mie_table=self%diag_MieTable(self%instance), km=self%km, klid=self%klid, nbegin=1, & - nbins=3, channels=self%diag_MieTable(self%instance)%channels*1.0e-9, & + call Aero_Compute_Diags (mie=self%diag_Mie, km=self%km, klid=self%klid, nbegin=1, & + nbins=3, & wavelengths_profile=self%wavelengths_profile*1.0e-9, & wavelengths_vertint=self%wavelengths_vertint*1.0e-9, & aerosol=aerosol, grav=MAPL_GRAV, tmpu=t, rhoa=airdens, rh=rh2, u=u, v=v, & @@ -1089,8 +1064,7 @@ subroutine aerosol_optics(state, rc) integer :: instance integer :: n, nbins integer :: i1, j1, i2, j2, km - integer :: band, offset - integer, parameter :: n_bands = 1 + integer :: band integer :: i, j, k @@ -1112,7 +1086,6 @@ subroutine aerosol_optics(state, rc) ! -------------- band = 0 call ESMF_AttributeGet(state, name='band_for_aerosol_optics', value=band, __RC__) - offset = band - n_bands ! Pressure at layer edges ! ------------------------ @@ -1159,7 +1132,7 @@ subroutine aerosol_optics(state, rc) address = transfer(opaque_self, address) call c_f_pointer(address, self) - call mie_ (self%rad_MieTable(instance), nbins, n_bands, offset, q_4d, rh, ext_s, ssa_s, asy_s, __RC__) + call mie_ (self%rad_Mie, nbins, band, q_4d, rh, ext_s, ssa_s, asy_s, __RC__) call ESMF_AttributeGet(state, name='extinction_in_air_due_to_ambient_aerosol', value=fld_name, __RC__) if (fld_name /= '') then @@ -1187,34 +1160,34 @@ subroutine aerosol_optics(state, rc) contains ! subroutine mie_(mie_table, aerosol_names, nb, offset, q, rh, bext_s, bssa_s, basym_s, rc) - subroutine mie_(mie_table, nbins, nb, offset, q, rh, bext_s, bssa_s, basym_s, rc) + subroutine mie_(mie, nbins, band, q, rh, bext_s, bssa_s, basym_s, rc) implicit none - type(Chem_Mie), intent(inout) :: mie_table ! mie table + type(GOCART2G_Mie) , intent(inout) :: mie ! mie table integer, intent(in ) :: nbins ! number of bins - integer, intent(in ) :: nb ! number of bands - integer, intent(in ) :: offset ! bands offset + integer, intent(in ) :: band ! channel real, intent(in ) :: q(:,:,:,:) ! aerosol mass mixing ratio, kg kg-1 real, intent(in ) :: rh(:,:,:) ! relative humidity real(kind=8), intent( out) :: bext_s (size(ext_s,1),size(ext_s,2),size(ext_s,3)) real(kind=8), intent( out) :: bssa_s (size(ext_s,1),size(ext_s,2),size(ext_s,3)) real(kind=8), intent( out) :: basym_s(size(ext_s,1),size(ext_s,2),size(ext_s,3)) - integer, intent(out) :: rc + integer, intent(out) :: rc ! local - integer :: l - real :: bext (size(ext_s,1),size(ext_s,2),size(ext_s,3)) ! extinction - real :: bssa (size(ext_s,1),size(ext_s,2),size(ext_s,3)) ! SSA - real :: gasym(size(ext_s,1),size(ext_s,2),size(ext_s,3)) ! asymmetry parameter + integer :: l + real :: bext (size(ext_s,1),size(ext_s,2),size(ext_s,3)) ! extinction + real :: bssa (size(ext_s,1),size(ext_s,2),size(ext_s,3)) ! SSA + real :: gasym(size(ext_s,1),size(ext_s,2),size(ext_s,3)) ! asymmetry parameter __Iam__('NI2G::aerosol_optics::mie_') - bext_s = 0.0d0 - bssa_s = 0.0d0 - basym_s = 0.0d0 + bext_s = 0.0d0 + bssa_s = 0.0d0 + basym_s = 0.0d0 do l = 1, nbins - call Chem_MieQuery(mie_table, l, real(offset+1.), q(:,:,:,l), rh, tau=bext, gasym=gasym, ssa=bssa) + !tau is converted to bext + call mie%Query(band, l, q(:,:,:,l), rh, tau=bext, gasym=gasym, ssa=bssa, __RC__) bext_s = bext_s + bext ! extinction bssa_s = bssa_s + (bssa*bext) ! scattering extinction @@ -1256,7 +1229,7 @@ subroutine monochromatic_aerosol_optics(state, rc) integer :: instance integer :: n, nbins integer :: i1, j1, i2, j2, km - real :: wavelength, mieTable_index + real :: wavelength integer :: i, j, k __Iam__('NI2G:: monochromatic_aerosol_optics') @@ -1277,21 +1250,6 @@ subroutine monochromatic_aerosol_optics(state, rc) ! -------------- call ESMF_AttributeGet(state, name='wavelength_for_aerosol_optics', value=wavelength, __RC__) -! Get wavelength index for Mie Table -! ---------------------------------- -! Channel values are 4.7e-7 5.5e-7 6.7e-7 8.7e-7 [meter]. Their indices are 1,2,3,4 respectively. - if ((wavelength .ge. 4.69e-7) .and. (wavelength .le. 4.71e-7)) then - mieTable_index = 1. - else if ((wavelength .ge. 5.49e-7) .and. (wavelength .le. 5.51e-7)) then - mieTable_index = 2. - else if ((wavelength .ge. 6.69e-7) .and. (wavelength .le. 6.71e-7)) then - mieTable_index = 3. - else if ((wavelength .ge. 8.68e-7) .and. (wavelength .le. 8.71e-7)) then - mieTable_index = 4. - else - print*,trim(Iam),' : wavelength of ',wavelength,' is an invalid value.' - return - end if ! Pressure at layer edges ! ------------------------ @@ -1338,14 +1296,8 @@ subroutine monochromatic_aerosol_optics(state, rc) call c_f_pointer(address, self) do n = 1, nbins - do i = 1, i2 - do j = 1, j2 - do k = 1, km - call Chem_MieQuery(self%diag_MieTable(instance), n, mieTable_index, q_4d(i,j,k,n), rh(i,j,k), tau=tau(i,j,k), __RC__) - tau_s(i,j,k) = tau_s(i,j,k) + tau(i,j,k) - end do - end do - end do + call self%diag_Mie%Query(wavelength, n, q_4d(:,:,:,n), rh, tau=tau, __RC__) + tau_s = tau_s + tau end do call ESMF_AttributeGet(state, name='monochromatic_extinction_in_air_due_to_ambient_aerosol', value=fld_name, __RC__) diff --git a/ESMF/GOCART2G_GridComp/SS2G_GridComp/SS2G_GridCompMod.F90 b/ESMF/GOCART2G_GridComp/SS2G_GridComp/SS2G_GridCompMod.F90 index badbff97..81e63e64 100644 --- a/ESMF/GOCART2G_GridComp/SS2G_GridComp/SS2G_GridCompMod.F90 +++ b/ESMF/GOCART2G_GridComp/SS2G_GridComp/SS2G_GridCompMod.F90 @@ -11,7 +11,7 @@ module SS2G_GridCompMod ! !USES: use ESMF use MAPL - use Chem_MieTableMod2G + use GOCART2G_MieMod use Chem_AeroGeneric use iso_c_binding, only: c_loc, c_f_pointer, c_ptr @@ -328,7 +328,9 @@ subroutine Initialize (GC, IMPORT, EXPORT, CLOCK, RC) integer :: NUM_BANDS logical :: bands_are_present real, pointer, dimension(:,:,:) :: ple - + integer, allocatable, dimension(:) :: channels_ + integer :: nmom_ + character(len=ESMF_MAXSTR) :: file_ __Iam__('Initialize') !**************************************************************************** @@ -464,48 +466,21 @@ subroutine Initialize (GC, IMPORT, EXPORT, CLOCK, RC) ! Create Radiation Mie Table ! -------------------------- - call MAPL_GetResource (MAPL, NUM_BANDS, 'NUM_BANDS:', __RC__) - -! Get file names for the optical tables - call ESMF_ConfigGetAttribute (cfg, self%rad_MieTable(instance)%optics_file, & - label="aerosol_radBands_optics_file:", __RC__ ) - - allocate (self%rad_MieTable(instance)%channels(NUM_BANDS), __STAT__ ) - self%rad_MieTable(instance)%nch = NUM_BANDS - - call ESMF_ConfigFindLabel(cfg, label="BANDS:", isPresent=bands_are_present, __RC__) - - if (bands_are_present) then - call ESMF_ConfigGetAttribute (cfg, self%rad_MieTable(instance)%channels, label= "BANDS:", & - count=self%rad_MieTable(instance)%nch, __RC__) - else - do i = 1, NUM_BANDS - self%rad_MieTable(instance)%channels(i) = i - end do - endif - - allocate (self%rad_MieTable(instance)%mie_aerosol, __STAT__) - self%rad_MieTable(instance)%mie_aerosol = Chem_MieTableCreate (self%rad_MieTable(instance)%optics_file, __RC__) - call Chem_MieTableRead (self%rad_MieTable(instance)%mie_aerosol, NUM_BANDS, self%rad_MieTable(instance)%channels, __RC__) + call ESMF_ConfigGetAttribute (cfg, file_, label="aerosol_radBands_optics_file:", __RC__ ) + self%rad_Mie = GOCART2G_Mie(trim(file_), __RC__) ! Create Diagnostics Mie Table ! ----------------------------- ! Get file names for the optical tables - call ESMF_ConfigGetAttribute (cfg, self%diag_MieTable(instance)%optics_file, & + call ESMF_ConfigGetAttribute (cfg, file_, & label="aerosol_monochromatic_optics_file:", __RC__ ) - call ESMF_ConfigGetAttribute (cfg, self%diag_MieTable(instance)%nmom, label="n_moments:", default=0, __RC__) + call ESMF_ConfigGetAttribute (cfg, nmom_, label="n_moments:", default=0, __RC__) i = ESMF_ConfigGetLen (universal_cfg, label='aerosol_monochromatic_optics_wavelength_in_nm_from_LUT:', __RC__) - self%diag_MieTable(instance)%nch = i - allocate (self%diag_MieTable(instance)%channels(self%diag_MieTable(instance)%nch), __STAT__ ) - call ESMF_ConfigGetAttribute (universal_cfg, self%diag_MieTable(instance)%channels, & + allocate (channels_(i), __STAT__ ) + call ESMF_ConfigGetAttribute (universal_cfg, channels_, & label= "aerosol_monochromatic_optics_wavelength_in_nm_from_LUT:", __RC__) - - allocate (self%diag_MieTable(instance)%mie_aerosol, __STAT__) - self%diag_MieTable(instance)%mie_aerosol = Chem_MieTableCreate (self%diag_MieTable(instance)%optics_file, __RC__ ) - call Chem_MieTableRead (self%diag_MieTable(instance)%mie_aerosol, self%diag_MieTable(instance)%nch, & - self%diag_MieTable(instance)%channels*1.e-9, rc=status, nmom=self%diag_MieTable(instance)%nmom) - VERIFY_(status) - + self%diag_Mie = GOCART2G_Mie(trim(file_), channels_*1.e-9, nmom=nmom_, __RC__) + deallocate(channels_) ! Mie Table instance/index call ESMF_AttributeSet(aero, name='mie_table_instance', value=instance, __RC__) @@ -831,8 +806,8 @@ subroutine Run2 (GC, import, export, clock, RC) ! Compute diagnostics ! ------------------- ! Certain variables are multiplied by 1.0e-9 to convert from nanometers to meters - call Aero_Compute_Diags (self%diag_MieTable(self%instance), self%km, self%klid, 1, self%nbins, self%rlow, & - self%rup, self%diag_MieTable(self%instance)%channels*1.0e-9, self%wavelengths_profile*1.0e-9, & + call Aero_Compute_Diags (self%diag_Mie, self%km, self%klid, 1, self%nbins, self%rlow, & + self%rup, self%wavelengths_profile*1.0e-9, & self%wavelengths_vertint*1.0e-9, SS, MAPL_GRAV, t, airdens,rh2, u, v, & delp, ple, tropp,SSSMASS, SSCMASS, SSMASS, SSEXTTAU,SSSTEXTTAU, SSSCATAU,SSSTSCATAU, & SSSMASS25, SSCMASS25, SSMASS25, SSEXTT25, SSSCAT25, & @@ -932,8 +907,7 @@ subroutine aerosol_optics(state, rc) integer :: instance integer :: n, nbins integer :: i1, j1, i2, j2, km - integer :: band, offset - integer, parameter :: n_bands = 1 + integer :: band integer :: k __Iam__('SS2G::aerosol_optics') @@ -948,7 +922,6 @@ subroutine aerosol_optics(state, rc) ! -------------- band = 0 call ESMF_AttributeGet(state, name='band_for_aerosol_optics', value=band, __RC__) - offset = band - n_bands ! Pressure at layer edges ! ------------------------ @@ -995,7 +968,7 @@ subroutine aerosol_optics(state, rc) address = transfer(opaque_self, address) call c_f_pointer(address, self) - call mie_ (self%rad_MieTable(instance), nbins, n_bands, offset, q_4d, rh, ext_s, ssa_s, asy_s, __RC__) + call mie_ (self%rad_Mie, nbins, band, q_4d, rh, ext_s, ssa_s, asy_s, __RC__) call ESMF_AttributeGet(state, name='extinction_in_air_due_to_ambient_aerosol', value=fld_name, __RC__) if (fld_name /= '') then @@ -1022,14 +995,13 @@ subroutine aerosol_optics(state, rc) contains - subroutine mie_(mie_table, nbins, nb, offset, q, rh, bext_s, bssa_s, basym_s, rc) + subroutine mie_(mie, nbins, band, q, rh, bext_s, bssa_s, basym_s, rc) implicit none - type(Chem_Mie), intent(inout) :: mie_table ! mie table + type(GOCART2G_Mie), intent(inout) :: mie ! mie table integer, intent(in ) :: nbins ! number of bins - integer, intent(in ) :: nb ! number of bands - integer, intent(in ) :: offset ! bands offset + integer, intent(in ) :: band ! channel real, intent(in ) :: q(:,:,:,:) ! aerosol mass mixing ratio, kg kg-1 real, intent(in ) :: rh(:,:,:) ! relative humidity real(kind=8), intent( out) :: bext_s (size(ext_s,1),size(ext_s,2),size(ext_s,3)) @@ -1045,12 +1017,13 @@ subroutine mie_(mie_table, nbins, nb, offset, q, rh, bext_s, bssa_s, basym_s, rc __Iam__('SS2G::aerosol_optics::mie_') - bext_s = 0.0d0 - bssa_s = 0.0d0 - basym_s = 0.0d0 + bext_s = 0.0d0 + bssa_s = 0.0d0 + basym_s = 0.0d0 do l = 1, nbins - call Chem_MieQuery(mie_table, l, real(offset+1.), q(:,:,:,l), rh, tau=bext, gasym=gasym, ssa=bssa) + ! tau is converted to bext + call mie%Query(band, l, q(:,:,:,l), rh, tau=bext, gasym=gasym, ssa=bssa, __RC__) bext_s = bext_s + bext ! extinction bssa_s = bssa_s + (bssa*bext) ! scattering extinction @@ -1087,7 +1060,7 @@ subroutine monochromatic_aerosol_optics(state, rc) integer :: instance integer :: n, nbins, k integer :: i1, j1, i2, j2, km, i, j - real :: wavelength, mieTable_index + real :: wavelength __Iam__('SS2G::monochromatic_aerosol_optics') @@ -1102,22 +1075,6 @@ subroutine monochromatic_aerosol_optics(state, rc) wavelength = 0. call ESMF_AttributeGet (state, name='wavelength_for_aerosol_optics', value=wavelength, __RC__) -! Get wavelength index for Mie Table -! ---------------------------------- -! Channel values are 4.7e-7 5.5e-7 6.7e-7 8.7e-7 [meter]. Their indices are 1,2,3,4 respectively. - if ((wavelength .ge. 4.69e-7) .and. (wavelength .le. 4.71e-7)) then - mieTable_index = 1. - else if ((wavelength .ge. 5.49e-7) .and. (wavelength .le. 5.51e-7)) then - mieTable_index = 2. - else if ((wavelength .ge. 6.69e-7) .and. (wavelength .le. 6.71e-7)) then - mieTable_index = 3. - else if ((wavelength .ge. 8.68e-7) .and. (wavelength .le. 8.71e-7)) then - mieTable_index = 4. - else - print*,trim(Iam),' : wavelength of ',wavelength,' is an invalid value.' - return - end if - ! Pressure at layer edges ! ------------------------ call ESMF_AttributeGet (state, name='air_pressure_for_aerosol_optics', value=fld_name, __RC__) @@ -1165,15 +1122,8 @@ subroutine monochromatic_aerosol_optics(state, rc) call c_f_pointer(address, self) do n = 1, nbins - do i = 1, i2 - do j = 1, j2 - do k = 1, km - call Chem_MieQuery(self%diag_MieTable(instance), n, mieTable_index, q_4d(i,j,k,n), rh(i,j,k), tau=tau(i,j,k), __RC__) - tau_s = tau_s + tau - tau = 0. - end do - end do - end do + call self%diag_Mie%Query(wavelength, n, q_4d(:,:,:,n), rh, tau=tau, __RC__) + tau_s = tau_s + tau end do call ESMF_AttributeGet (state, name='monochromatic_extinction_in_air_due_to_ambient_aerosol', value=fld_name, __RC__) diff --git a/ESMF/GOCART2G_GridComp/SU2G_GridComp/SU2G_GridCompMod.F90 b/ESMF/GOCART2G_GridComp/SU2G_GridComp/SU2G_GridCompMod.F90 index 546fde66..09a21a97 100644 --- a/ESMF/GOCART2G_GridComp/SU2G_GridComp/SU2G_GridCompMod.F90 +++ b/ESMF/GOCART2G_GridComp/SU2G_GridComp/SU2G_GridCompMod.F90 @@ -11,7 +11,7 @@ module SU2G_GridCompMod ! !USES: use ESMF use MAPL - use Chem_MieTableMod2G + use GOCART2G_MieMod use Chem_AeroGeneric use iso_c_binding, only: c_loc, c_f_pointer, c_ptr @@ -401,7 +401,9 @@ subroutine Initialize (GC, IMPORT, EXPORT, CLOCK, RC) integer :: year, month, day, hh, mm, ss real, dimension(4) :: Vect_Hcts - + integer, allocatable, dimension(:) :: channels_ + integer :: nmom_ + character(len=ESMF_MAXSTR) :: file_ __Iam__('Initialize') !**************************************************************************** @@ -576,48 +578,21 @@ subroutine Initialize (GC, IMPORT, EXPORT, CLOCK, RC) ! Create Radiation Mie Table ! -------------------------- - call MAPL_GetResource (MAPL, NUM_BANDS, 'NUM_BANDS:', __RC__) - -! Get file names for the optical tables - call ESMF_ConfigGetAttribute (cfg, self%rad_MieTable(instance)%optics_file, & - label="aerosol_radBands_optics_file:", __RC__ ) - - allocate (self%rad_MieTable(instance)%channels(NUM_BANDS), __STAT__ ) - self%rad_MieTable(instance)%nch = NUM_BANDS - - call ESMF_ConfigFindLabel(cfg, label="BANDS:", isPresent=bands_are_present, __RC__) - - if (bands_are_present) then - call ESMF_ConfigGetAttribute (cfg, self%rad_MieTable(instance)%channels, label= "BANDS:", & - count=self%rad_MieTable(instance)%nch, __RC__) - else - do i = 1, NUM_BANDS - self%rad_MieTable(instance)%channels(i) = i - end do - endif - - allocate (self%rad_MieTable(instance)%mie_aerosol, __STAT__) - self%rad_MieTable(instance)%mie_aerosol = Chem_MieTableCreate (self%rad_MieTable(instance)%optics_file, __RC__) - call Chem_MieTableRead (self%rad_MieTable(instance)%mie_aerosol, NUM_BANDS, self%rad_MieTable(instance)%channels, __RC__) + call ESMF_ConfigGetAttribute (cfg, file_, label="aerosol_radBands_optics_file:", __RC__ ) + self%rad_Mie = GOCART2G_Mie(trim(file_), __RC__) ! Create Diagnostics Mie Table ! ----------------------------- ! Get file names for the optical tables - call ESMF_ConfigGetAttribute (cfg, self%diag_MieTable(instance)%optics_file, & + call ESMF_ConfigGetAttribute (cfg, file_, & label="aerosol_monochromatic_optics_file:", __RC__ ) - call ESMF_ConfigGetAttribute (cfg, self%diag_MieTable(instance)%nmom, label="n_moments:", default=0, __RC__) - + call ESMF_ConfigGetAttribute (cfg, nmom_, label="n_moments:", default=0, __RC__) i = ESMF_ConfigGetLen (universal_cfg, label='aerosol_monochromatic_optics_wavelength_in_nm_from_LUT:', __RC__) - self%diag_MieTable(instance)%nch = i - allocate (self%diag_MieTable(instance)%channels(self%diag_MieTable(instance)%nch), __STAT__ ) - call ESMF_ConfigGetAttribute (universal_cfg, self%diag_MieTable(instance)%channels, & + allocate (channels_(i), __STAT__ ) + call ESMF_ConfigGetAttribute (universal_cfg, channels_, & label= "aerosol_monochromatic_optics_wavelength_in_nm_from_LUT:", __RC__) - - allocate (self%diag_MieTable(instance)%mie_aerosol, __STAT__) - self%diag_MieTable(instance)%mie_aerosol = Chem_MieTableCreate (self%diag_MieTable(instance)%optics_file, __RC__ ) - call Chem_MieTableRead (self%diag_MieTable(instance)%mie_aerosol, self%diag_MieTable(instance)%nch, & - self%diag_MieTable(instance)%channels*1.e-9, rc=status, nmom=self%diag_MieTable(instance)%nmom) - VERIFY_(status) + self%diag_Mie = GOCART2G_Mie(trim(file_), channels_*1.e-9, nmom=nmom_, __RC__) + deallocate(channels_) ! Mie Table instance/index call ESMF_AttributeSet(aero, name='mie_table_instance', value=instance, __RC__) @@ -1108,8 +1083,7 @@ subroutine Run2 (GC, import, export, clock, RC) ! Certain variables are multiplied by 1.0e-9 to convert from nanometers to meters call SU_Compute_Diags ( self%km, self%klid, self%radius(nSO4), self%sigma(nSO4), self%rhop(nSO4), & - MAPL_GRAV, MAPL_PI, nSO4, self%diag_MieTable(self%instance), & - self%diag_MieTable(self%instance)%channels*1.0e-9, & + MAPL_GRAV, MAPL_PI, nSO4, self%diag_Mie, & self%wavelengths_profile*1.0e-9, self%wavelengths_vertint*1.0e-9, & t, airdens, delp, ple,tropp, rh2, u, v, DMS, SO2, SO4, dummyMSA, & DMSSMASS, DMSCMASS, & @@ -1206,8 +1180,7 @@ subroutine aerosol_optics(state, rc) integer :: instance integer :: n, nbins integer :: i1, j1, i2, j2, km - integer :: band, offset - integer, parameter :: n_bands = 1 + integer :: band integer :: i, j, k @@ -1229,7 +1202,6 @@ subroutine aerosol_optics(state, rc) ! -------------- band = 0 call ESMF_AttributeGet(state, name='band_for_aerosol_optics', value=band, __RC__) - offset = band - n_bands ! Pressure at layer edges ! ------------------------ @@ -1276,7 +1248,7 @@ subroutine aerosol_optics(state, rc) address = transfer(opaque_self, address) call c_f_pointer(address, self) - call mie_ (self%rad_MieTable(instance), nbins, n_bands, offset, q_4d, rh, ext_s, ssa_s, asy_s, __RC__) + call mie_ (self%rad_Mie, nbins, band, q_4d, rh, ext_s, ssa_s, asy_s, __RC__) call ESMF_AttributeGet(state, name='extinction_in_air_due_to_ambient_aerosol', value=fld_name, __RC__) if (fld_name /= '') then @@ -1304,14 +1276,13 @@ subroutine aerosol_optics(state, rc) contains ! subroutine mie_(mie_table, aerosol_names, nb, offset, q, rh, bext_s, bssa_s, basym_s, rc) - subroutine mie_(mie_table, nbins, nb, offset, q, rh, bext_s, bssa_s, basym_s, rc) + subroutine mie_(mie, nbins, band, q, rh, bext_s, bssa_s, basym_s, rc) implicit none - type(Chem_Mie), intent(inout) :: mie_table ! mie table + type(GOCART2G_Mie), intent(inout) :: mie ! mie table integer, intent(in ) :: nbins ! number of bins - integer, intent(in ) :: nb ! number of bands - integer, intent(in ) :: offset ! bands offset + integer, intent(in ) :: band ! channel real, intent(in ) :: q(:,:,:,:) ! aerosol mass mixing ratio, kg kg-1 real, intent(in ) :: rh(:,:,:) ! relative humidity real(kind=8), intent( out) :: bext_s (size(ext_s,1),size(ext_s,2),size(ext_s,3)) @@ -1324,15 +1295,15 @@ subroutine mie_(mie_table, nbins, nb, offset, q, rh, bext_s, bssa_s, basym_s, rc real :: bssa (size(ext_s,1),size(ext_s,2),size(ext_s,3)) ! SSA real :: gasym(size(ext_s,1),size(ext_s,2),size(ext_s,3)) ! asymmetry parameter - __Iam__('SU2G::aerosol_optics::mie_') - bext_s = 0.0d0 - bssa_s = 0.0d0 - basym_s = 0.0d0 + bext_s = 0.0d0 + bssa_s = 0.0d0 + basym_s = 0.0d0 do l = 1, nbins - call Chem_MieQuery(mie_table, l, real(offset+1.), q(:,:,:,l), rh, tau=bext, gasym=gasym, ssa=bssa) + ! tau is converted to bext + call mie%Query(band, l, q(:,:,:,l), rh, tau=bext, gasym=gasym, ssa=bssa, __RC__) bext_s = bext_s + bext ! extinction bssa_s = bssa_s + (bssa*bext) ! scattering extinction @@ -1373,7 +1344,7 @@ subroutine monochromatic_aerosol_optics(state, rc) integer :: instance integer :: n, nbins integer :: i1, j1, i2, j2, km - real :: wavelength, mieTable_index + real :: wavelength integer :: i, j, k __Iam__('SU2G::monochromatic_aerosol_optics') @@ -1394,21 +1365,6 @@ subroutine monochromatic_aerosol_optics(state, rc) ! -------------- call ESMF_AttributeGet(state, name='wavelength_for_aerosol_optics', value=wavelength, __RC__) -! Get wavelength index for Mie Table -! ---------------------------------- -! Channel values are 4.7e-7 5.5e-7 6.7e-7 8.7e-7 [meter]. Their indices are 1,2,3,4 respectively. - if ((wavelength .ge. 4.69e-7) .and. (wavelength .le. 4.71e-7)) then - mieTable_index = 1. - else if ((wavelength .ge. 5.49e-7) .and. (wavelength .le. 5.51e-7)) then - mieTable_index = 2. - else if ((wavelength .ge. 6.69e-7) .and. (wavelength .le. 6.71e-7)) then - mieTable_index = 3. - else if ((wavelength .ge. 8.68e-7) .and. (wavelength .le. 8.71e-7)) then - mieTable_index = 4. - else - print*,trim(Iam),' : wavelengths of ',wavelength,' is an invalid value.' - return - end if ! Pressure at layer edges ! ------------------------ @@ -1429,7 +1385,7 @@ subroutine monochromatic_aerosol_optics(state, rc) allocate(tau_s(i1:i2, j1:j2, km), & tau(i1:i2, j1:j2, km), __STAT__) tau_s = 0.0 - tau = 0.0 + tau = 0.0 allocate(q_4d(i1:i2, j1:j2, km, nbins), __STAT__) @@ -1455,14 +1411,8 @@ subroutine monochromatic_aerosol_optics(state, rc) call c_f_pointer(address, self) do n = 1, nbins - do i = 1, i2 - do j = 1, j2 - do k = 1, km - call Chem_MieQuery(self%diag_MieTable(instance), n, mieTable_index, q_4d(i,j,k,n), rh(i,j,k), tau=tau(i,j,k), __RC__) - tau_s(i,j,k) = tau_s(i,j,k) + tau(i,j,k) - end do - end do - end do + call self%diag_Mie%Query(wavelength, n, q_4d(:,:,:,n), rh, tau=tau, __RC__) + tau_s = tau_s + tau end do call ESMF_AttributeGet(state, name='monochromatic_extinction_in_air_due_to_ambient_aerosol', value=fld_name, __RC__) diff --git a/ESMF/GOCART_GridComp/CFC_GridComp/CMakeLists.txt b/ESMF/GOCART_GridComp/CFC_GridComp/CMakeLists.txt index 7d5828fc..7284fa66 100644 --- a/ESMF/GOCART_GridComp/CFC_GridComp/CMakeLists.txt +++ b/ESMF/GOCART_GridComp/CFC_GridComp/CMakeLists.txt @@ -2,8 +2,7 @@ esma_set_this () esma_add_library (${this} SRCS ${this}Mod.F90 - DEPENDENCIES Chem_Shared Chem_Base GMAO_mpeu MAPL_cfio_r4) -target_include_directories (${this} PUBLIC ${INC_ESMF} ${INC_NETCDF}) + DEPENDENCIES Chem_Shared Chem_Base GMAO_mpeu MAPL_cfio_r4 esmf NetCDF::NetCDF_Fortran) esma_generate_gocart_code (${this} "-B\;-C\;-N\;GOCART") diff --git a/ESMF/GOCART_GridComp/CH4_GridComp/CMakeLists.txt b/ESMF/GOCART_GridComp/CH4_GridComp/CMakeLists.txt index 716b598c..2d73b255 100644 --- a/ESMF/GOCART_GridComp/CH4_GridComp/CMakeLists.txt +++ b/ESMF/GOCART_GridComp/CH4_GridComp/CMakeLists.txt @@ -2,9 +2,8 @@ esma_set_this () esma_add_library (${this} SRCS ${this}Mod.F90 - DEPENDENCIES Chem_Shared Chem_Base GMAO_mpeu MAPL_cfio_r4 + DEPENDENCIES Chem_Shared Chem_Base GMAO_mpeu MAPL_cfio_r4 NetCDF::NetCDF_Fortran ) -target_include_directories (${this} PUBLIC ${INC_NETCDF}) esma_generate_gocart_code (${this} "-B\;-E\;-F\;-N\;GOCART") file (GLOB_RECURSE rc_files CONFIGURE_DEPENDS RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.rc) diff --git a/ESMF/GOCART_GridComp/CMakeLists.txt b/ESMF/GOCART_GridComp/CMakeLists.txt index ffd39452..cce3ec13 100644 --- a/ESMF/GOCART_GridComp/CMakeLists.txt +++ b/ESMF/GOCART_GridComp/CMakeLists.txt @@ -16,16 +16,15 @@ set (srcs set (resource_files GOCARTdata_AerRegistry.rc - ) + ) -install( FILES ${resource_files} +install( FILES ${resource_files} DESTINATION etc ) -set (dependencies Chem_Base Chem_Shared MAPL GMAO_mpeu) +set (dependencies Chem_Base Chem_Shared MAPL GMAO_mpeu esmf) esma_add_library (${this} SRCS ${srcs} SUBCOMPONENTS ${alldirs} - DEPENDENCIES ${dependencies} - INCLUDES ${INC_ESMF}) + DEPENDENCIES ${dependencies}) diff --git a/ESMF/GOCART_GridComp/CO2_GridComp/CMakeLists.txt b/ESMF/GOCART_GridComp/CO2_GridComp/CMakeLists.txt index 7d5828fc..7284fa66 100644 --- a/ESMF/GOCART_GridComp/CO2_GridComp/CMakeLists.txt +++ b/ESMF/GOCART_GridComp/CO2_GridComp/CMakeLists.txt @@ -2,8 +2,7 @@ esma_set_this () esma_add_library (${this} SRCS ${this}Mod.F90 - DEPENDENCIES Chem_Shared Chem_Base GMAO_mpeu MAPL_cfio_r4) -target_include_directories (${this} PUBLIC ${INC_ESMF} ${INC_NETCDF}) + DEPENDENCIES Chem_Shared Chem_Base GMAO_mpeu MAPL_cfio_r4 esmf NetCDF::NetCDF_Fortran) esma_generate_gocart_code (${this} "-B\;-C\;-N\;GOCART") diff --git a/ESMF/GOCART_GridComp/CO_GridComp/CMakeLists.txt b/ESMF/GOCART_GridComp/CO_GridComp/CMakeLists.txt index c19b9c6d..2f6da177 100644 --- a/ESMF/GOCART_GridComp/CO_GridComp/CMakeLists.txt +++ b/ESMF/GOCART_GridComp/CO_GridComp/CMakeLists.txt @@ -2,8 +2,7 @@ esma_set_this () esma_add_library (${this} SRCS ${this}Mod.F90 - DEPENDENCIES Chem_Shared Chem_Base GMAO_mpeu MAPL_cfio_r4) -target_include_directories (${this} PUBLIC ${INC_ESMF} ${INC_NETCDF}) + DEPENDENCIES Chem_Shared Chem_Base GMAO_mpeu MAPL_cfio_r4 esmf NetCDF::NetCDF_Fortran) esma_generate_gocart_code (${this} "-B\;-E\;-F\;-N\;GOCART") diff --git a/ESMF/GOCART_GridComp/O3_GridComp/CMakeLists.txt b/ESMF/GOCART_GridComp/O3_GridComp/CMakeLists.txt index 0debe058..fd882c84 100644 --- a/ESMF/GOCART_GridComp/O3_GridComp/CMakeLists.txt +++ b/ESMF/GOCART_GridComp/O3_GridComp/CMakeLists.txt @@ -2,8 +2,7 @@ esma_set_this() esma_add_library (${this} SRCS ${this}Mod.F90 - DEPENDENCIES Chem_Shared Chem_Base GMAO_mpeu MAPL_cfio_r4) -target_include_directories (${this} PUBLIC ${INC_NETCDF}) + DEPENDENCIES Chem_Shared Chem_Base GMAO_mpeu MAPL_cfio_r4 NetCDF::NetCDF_Fortran) esma_generate_gocart_code (${this} "-B\;-C\;-N\;GOCART") diff --git a/ESMF/GOCART_GridComp/Rn_GridComp/CMakeLists.txt b/ESMF/GOCART_GridComp/Rn_GridComp/CMakeLists.txt index c19b9c6d..2f6da177 100644 --- a/ESMF/GOCART_GridComp/Rn_GridComp/CMakeLists.txt +++ b/ESMF/GOCART_GridComp/Rn_GridComp/CMakeLists.txt @@ -2,8 +2,7 @@ esma_set_this () esma_add_library (${this} SRCS ${this}Mod.F90 - DEPENDENCIES Chem_Shared Chem_Base GMAO_mpeu MAPL_cfio_r4) -target_include_directories (${this} PUBLIC ${INC_ESMF} ${INC_NETCDF}) + DEPENDENCIES Chem_Shared Chem_Base GMAO_mpeu MAPL_cfio_r4 esmf NetCDF::NetCDF_Fortran) esma_generate_gocart_code (${this} "-B\;-E\;-F\;-N\;GOCART") diff --git a/Process_Library/CMakeLists.txt b/Process_Library/CMakeLists.txt index be48bcee..4ca4da3c 100644 --- a/Process_Library/CMakeLists.txt +++ b/Process_Library/CMakeLists.txt @@ -1,7 +1,7 @@ esma_set_this () set (srcs - Chem_MieTableMod2G.F90 + GOCART2G_MieMod.F90 GOCART2G_Process.F90 ) diff --git a/Process_Library/Chem_MieTableMod2G.F90 b/Process_Library/Chem_MieTableMod2G.F90 deleted file mode 100644 index 2fc9514e..00000000 --- a/Process_Library/Chem_MieTableMod2G.F90 +++ /dev/null @@ -1,866 +0,0 @@ -!BOP -! -! !MODULE: Chem_MieTableMod --- Reader for aerosol mie tables -! -! !INTERFACE: -! - - module Chem_MieTableMod2G - -! !USES: - implicit none - include "netcdf.inc" ! Required for Mie tables stored as NCDF files - -! !PUBLIC TYPES: -! - private - public Chem_MieTable ! Holds Mie Lookup Tables - public Chem_Mie ! Holds Chem_MieTable and other information - -! -! !PUBLIC MEMBER FUNCTIONS: -! - public Chem_MieTableCreate ! Constructor - public Chem_MieTableRead ! Read the mie table from the file - public Chem_MieQuery -! -! !DESCRIPTION: -! -! This module read the mie aerosol tables. -! -! !REVISION HISTORY: -! -! 23Mar2005 Colarco - Initial code. -! 31Mar2005 Todling - Declared netcdf nf_ routines as external (OSF1) -! Removed # from include netcdf.inc -!EOP -!------------------------------------------------------------------------- -! -! Mie LUT table -! Will be reduced from input files to the desired channels -! -------- - - type Chem_MieTable - character(len=255) :: mietablename - integer :: nlambda ! number of wavelengths in table - integer :: nrh ! number of RH values in table - integer :: nbin ! number of size bins in table - integer :: nMom ! number of moments of phase function - integer :: nPol ! number of elements of scattering phase matrix - real, pointer :: lambda(:) => null() ! wavelengths [m] - real, pointer :: rh(:) => null() ! RH values [fraction] - real, pointer :: reff(:,:) => null() ! effective radius [m] - real, pointer :: bext(:,:,:) => null() ! bext values [m2 kg-1] - real, pointer :: bsca(:,:,:) => null() ! bsca values [m2 kg-1] - real, pointer :: bbck(:,:,:) => null() ! bbck values [m2 kg-1] - real, pointer :: g(:,:,:) => null() ! asymmetry parameter - real, pointer :: pback(:,:,:,:) => null() ! Backscatter phase function - real, pointer :: pmom(:,:,:,:,:) => null( ) ! moments of phase function - real, pointer :: gf(:,:) => null() ! hygroscopic growth factor - real, pointer :: rhop(:,:) => null() ! wet particle density [kg m-3] - real, pointer :: rhod(:,:) => null() ! wet particle density [kg m-3] - real, pointer :: vol(:,:) => null() ! wet particle volume [m3 kg-1] - real, pointer :: area(:,:) => null() ! wet particle cross section [m2 kg-1] - real, pointer :: refr(:,:,:) => null() ! real part of refractive index - real, pointer :: refi(:,:,:) => null() ! imaginary part of refractive index - - integer :: rhi(991) ! pointer to rh map - real :: rha(991) ! slope on rh map - end type Chem_MieTable - - type Chem_Mie - integer :: nch ! number of channels - integer :: nMom=0 ! number of moments (phase function) - integer :: nPol=0 ! number of moments (phase function) - real, pointer :: channels(:) ! wavelengths - - character(len=255) :: rcfile - character(len=255) :: optics_file - ! mie tables -- dim(nch,nrh,nbin) - type(Chem_MieTable), pointer :: mie_aerosol => null() - - integer :: nq ! number of tracers - character(len=255), pointer :: vname(:) => null() ! possibly remove lines 67-71 vname,vindex,vtable - integer, pointer :: vindex(:) => null() - type(Chem_MieTable), pointer :: vtable(:) => null() - ! mapping of vtable for given idx - type(Chem_MieTable), pointer :: vtableUse => null() - end type Chem_Mie - -# ifndef HAS_NETCDF3 - external nf_open, nf_inq_dimid, nf_inq_dimlen, nf_inq_varid, & - nf_get_var_double, nf_close -#endif - - - interface Chem_MieQuery - module procedure Chem_MieQueryByInt - module procedure Chem_MieQueryByIntWithpmom !possibly remove - end interface - - -CONTAINS - -!------------------------------------------------------------------------- -! NASA/GSFC, Global Modeling and Assimilation Office, Code 900.3 ! -!------------------------------------------------------------------------- -!BOP -! -! !IROUTINE: Chem_MieTableCreate --- Construct Chemistry Registry -! -! !INTERFACE: -! - - Function Chem_MieTableCreate ( rcfile, rc ) - - implicit none - type(Chem_MieTable) Chem_MieTableCreate - -! !INPUT PARAMETERS: - - character(len=*) :: rcfile ! Mie table file name - -! !OUTPUT PARAMETERS: - - integer, intent(out) :: rc ! Error return code: - ! 0 - all is well - ! 1 - - -! !DESCRIPTION: -! -! -! !REVISION HISTORY: -! -! 09Mar2005 da Silva API, prologues. -! -!EOP -!------------------------------------------------------------------------- - character(len=*), parameter :: myname = 'Chem_MieTableCreate' - - type(Chem_MieTable) :: this - -! _Iam_("Chem_MieTableCreate") - - rc = 0 - - this%mietablename = rcfile - -! Note: The actual allocation is done when reading because dimensions are -! read from file - -! All done -! -------- - Chem_MieTableCreate = this - - return - - end Function Chem_MieTableCreate - - -!------------------------------------------------------------------------- -! NASA/GSFC, Global Modeling and Assimilation Office, Code 900.3 ! -!------------------------------------------------------------------------- -!BOP -! -! !IROUTINE: Chem_MieTableRead --- Read and fill in the Mie table, interpolated -! to the requested channels -! -! !INTERFACE: -! - SUBROUTINE Chem_MieTableRead ( this, nch, channels, rc, nmom ) - -! !INPUT PARAMETERS: - - IMPLICIT none - TYPE(Chem_MieTable), intent(inout) :: this - integer, intent(in) :: nch ! number of channels to interpolate table to - real, intent(in) :: channels(:) ! channels to interpolate table to - integer, OPTIONAL, intent(in) :: nmom ! number of moments to keep (default=0) - integer, intent(out) :: rc ! return code - - -! !DESCRIPTION: -! -! Fills in the Mie table -! -! !REVISION HISTORY: -! -! 23Mar2005 Colarco -! -!EOP -!------------------------------------------------------------------------- - - character(len=*), parameter :: myname = 'Chem_MieTableRead' - - integer :: ncid, idimid, ivarid, n, i, j, ip1 - integer :: nch_table, nrh_table, nbin_table, nmom_table, nPol_table -! Tables are hard-wired as single precision - real*8, pointer :: channels_table(:), rh_table(:), reff_table(:,:), & - bext_table(:,:,:), bsca_table(:,:,:), & - bbck_table(:,:,:), g_table(:,:,:), & - pmom_table(:,:,:,:,:), pback_table(:,:,:,:), & - gf_table(:,:), rhop_table(:,:), rhod_table(:,:), & - vol_table(:,:), area_table(:,:), & - refr_table(:,:,:), refi_table(:,:,:) - - real :: yerr - integer :: nmom_, imom, ipol - integer :: status - -#define NF_VERIFY_(expr) rc = expr; if (rc /= 0) return -#define __NF_STAT__ stat=status); NF_VERIFY_(status - - rc = 0 - -! Whether or not doing phase function -! ----------------------------------- - if ( present(nmom) ) then - nmom_ = nmom - else - nmom_ = 0 - end if - -! Set up nPol_table for reading backscatter phase function -! This will get overwritten if pmoments is requested -! -------------------------------------------------------- - nPol_table = 6 - -! Open the table and get the dimensions -! ------------------------------------- - rc = nf_open(this%mietablename, NF_NOWRITE, ncid) - IF ( rc /= 0 ) THEN - print *, 'nf_open '//this%mietablename//' RETURN CODE=', rc - NF_VERIFY_(rc) - END IF - -! RH -! -- - NF_VERIFY_(nf_inq_dimid(ncid,'rh',idimid)) - NF_VERIFY_(nf_inq_dimlen(ncid,idimid,nrh_table)) - -! Channels -! -------- - NF_VERIFY_(nf_inq_dimid(ncid,'lambda',idimid)) - NF_VERIFY_(nf_inq_dimlen(ncid,idimid,nch_table)) - -! Dry Effective radius -! -------------------- - NF_VERIFY_(nf_inq_dimid(ncid,'radius',idimid)) - NF_VERIFY_(nf_inq_dimlen(ncid,idimid,nbin_table)) - -! Moments of phase function -! ------------------------- - if ( nmom_ > 0 ) then - NF_VERIFY_(nf_inq_dimid(ncid,'nMom',idimid)) - NF_VERIFY_(nf_inq_dimlen(ncid,idimid,nmom_table)) - if ( nmom_ > nmom_table ) then -! rc = 99 - print*,'Error: nmom_ > nmom_table, see:'//myname - NF_VERIFY_(1) - end if - NF_VERIFY_(nf_inq_dimid(ncid,'nPol',idimid)) - NF_VERIFY_(nf_inq_dimlen(ncid,idimid,nPol_table)) - endif - -! Get the table contents -! ------------------------------------- -! allocate ( channels_table(nch_table), rh_table(nrh_table), & -! bext_table(nch_table,nrh_table,nbin_table), & -! bsca_table(nch_table,nrh_table,nbin_table), & -! bbck_table(nch_table,nrh_table,nbin_table), & -! g_table(nch_table,nrh_table,nbin_table), stat = rc ) - - allocate(channels_table(nch_table), __NF_STAT__) - allocate(rh_table(nrh_table), __NF_STAT__) - allocate(reff_table(nrh_table,nbin_table), __NF_STAT__) - allocate(bext_table(nch_table,nrh_table,nbin_table), __NF_STAT__) - allocate(bsca_table(nch_table,nrh_table,nbin_table), __NF_STAT__) - allocate(bbck_table(nch_table,nrh_table,nbin_table), __NF_STAT__) - allocate(g_table(nch_table,nrh_table,nbin_table), stat = rc ) - allocate(pback_table(nch_table,nrh_table,nbin_table,nPol_table), __NF_STAT__) - allocate(gf_table(nrh_table,nbin_table), __NF_STAT__) - allocate(rhop_table(nrh_table,nbin_table), __NF_STAT__) - allocate(rhod_table(nrh_table,nbin_table), __NF_STAT__) - allocate(vol_table(nrh_table,nbin_table), __NF_STAT__) - allocate(area_table(nrh_table,nbin_table), __NF_STAT__) - allocate(refr_table(nch_table,nrh_table,nbin_table), __NF_STAT__) - allocate(refi_table(nch_table,nrh_table,nbin_table), __NF_STAT__) - - if ( nmom_ > 0 ) then - allocate(pmom_table(nch_table,nrh_table,nbin_table,nmom_table,nPol_table), __NF_STAT__) - end if - NF_VERIFY_(nf_inq_varid(ncid,'lambda',ivarid)) - NF_VERIFY_(nf_get_var_double(ncid,ivarid,channels_table)) - NF_VERIFY_(nf_inq_varid(ncid,'rEff',ivarid)) - NF_VERIFY_(nf_get_var_double(ncid,ivarid,reff_table)) - NF_VERIFY_(nf_inq_varid(ncid,'bext',ivarid)) - NF_VERIFY_(nf_get_var_double(ncid,ivarid,bext_table)) - NF_VERIFY_(nf_inq_varid(ncid,'bsca',ivarid)) - NF_VERIFY_(nf_get_var_double(ncid,ivarid,bsca_table)) - NF_VERIFY_(nf_inq_varid(ncid,'bbck',ivarid)) - NF_VERIFY_(nf_get_var_double(ncid,ivarid,bbck_table)) - NF_VERIFY_(nf_inq_varid(ncid,'g',ivarid)) - NF_VERIFY_(nf_get_var_double(ncid,ivarid,g_table)) - NF_VERIFY_(nf_inq_varid(ncid,'rh',ivarid)) - NF_VERIFY_(nf_get_var_double(ncid,ivarid,rh_table)) - - ! TODO: we need to look at these NF_NOERR checks -! Get the backscatter phase function values - rc = nf_inq_varid(ncid,'pback',ivarid) - if(rc .ne. NF_NOERR) then ! pback not in table, fill in dummy variable - pback_table = 1. - else - NF_VERIFY_(nf_get_var_double(ncid,ivarid,pback_table)) - endif - - if ( nmom_ > 0 ) then - NF_VERIFY_(nf_inq_varid(ncid,'pmom',ivarid)) - NF_VERIFY_(nf_get_var_double(ncid,ivarid,pmom_table)) - end if - -! Aerosol optical properties not necessarily stored in all versions of the tables -! ---------------------- -! Particle growth factor - rc = nf_inq_varid(ncid,'growth_factor',ivarid) - if(rc .ne. NF_NOERR) then ! not in table, fill in dummy variable - gf_table = -999. - else - NF_VERIFY_(nf_get_var_double(ncid,ivarid,gf_table)) - endif - -! Wet particle density - rc = nf_inq_varid(ncid,'rhop',ivarid) - if(rc .ne. NF_NOERR) then ! not in table, fill in dummy variable - rhop_table = -999. - else - NF_VERIFY_(nf_get_var_double(ncid,ivarid,rhop_table)) - endif - -! Dry particle density (will be pulled from wet particle radius) - rc = nf_inq_varid(ncid,'rhop',ivarid) - if(rc .ne. NF_NOERR) then ! not in table, fill in dummy variable - rhod_table = -999. - else - rc = nf_get_var_double(ncid,ivarid,rhod_table) - do i = 1, nrh_table - rhod_table(i,:) = rhod_table(1,:) - enddo - if (rc /=0) return - endif - -! Wet particle real part of refractive index - rc = nf_inq_varid(ncid,'refreal',ivarid) - if(rc .ne. NF_NOERR) then ! not in table, fill in dummy variable - refr_table = -999. - else - NF_VERIFY_(nf_get_var_double(ncid,ivarid,refr_table)) - endif - -! Wet particle imaginary part of refractive index (ensure positive) - rc = nf_inq_varid(ncid,'refimag',ivarid) - if(rc .ne. NF_NOERR) then ! not in table, fill in dummy variable - refi_table = -999. - else - NF_VERIFY_(nf_get_var_double(ncid,ivarid,refi_table)) - refi_table = abs(refi_table) - endif - -! Wet particle volume [m3 kg-1] -! Ratio of wet to dry volume is gf^3, hence the following - vol_table = gf_table**3 / rhod_table - -! Wet particle cross sectional area [m2 kg-1] -! Assume area is volume divided by (4./3.*reff) - area_table = vol_table / (4./3.*reff_table) - -! Close the table file -! ------------------------------------- - NF_VERIFY_(nf_close(ncid)) - - -! Setup the table to be returned -! ------------------------------------- - this%nlambda = nch - this%nrh = nrh_table - this%nbin = nbin_table - this%nMom = nmom_ -! if ( nmom_ > 0 ) this%nPol = nPol_table - this%nPol = nPol_table - -! allocate ( this%lambda(this%nLambda), this%rh(this%nrh), & -! this%bext(this%nLambda,this%nrh,this%nbin), & -! this%bsca(this%nLambda,this%nrh,this%nbin), & -! this%bbck(this%nLambda,this%nrh,this%nbin), & -! this%g(this%nLambda,this%nrh,this%nbin), & -! stat = rc ) - - allocate (this%lambda(this%nLambda), __NF_STAT__) - allocate (this%rh(this%nrh), __NF_STAT__) - allocate (this%reff(this%nrh,this%nbin), __NF_STAT__) - allocate (this%bext(this%nrh,this%nLambda,this%nbin), __NF_STAT__) - allocate (this%bsca(this%nrh,this%nLambda,this%nbin), __NF_STAT__) - allocate (this%bbck(this%nrh,this%nLambda,this%nbin), __NF_STAT__) - allocate (this%g(this%nrh,this%nLambda,this%nbin), __NF_STAT__) - allocate (this%pback(this%nrh,this%nLambda,this%nbin,this%nPol), __NF_STAT__) - if ( nmom_ > 0 ) then - allocate (this%pmom(this%nrh,this%nLambda,this%nbin,this%nMom,this%nPol), __NF_STAT__) - end if - allocate (this%gf(this%nrh,this%nbin), __NF_STAT__) - allocate (this%rhop(this%nrh,this%nbin), __NF_STAT__) - allocate (this%rhod(this%nrh,this%nbin), __NF_STAT__) - allocate (this%vol(this%nrh,this%nbin), __NF_STAT__) - allocate (this%area(this%nrh,this%nbin), __NF_STAT__) - allocate (this%refr(this%nrh,this%nLambda,this%nbin), __NF_STAT__) - allocate (this%refi(this%nrh,this%nLambda,this%nbin), __NF_STAT__) - -! Preserve the full RH structure of the input table - this%rh(:) = rh_table(:) - -! Insert the requested channels in the output table - this%lambda(:) = channels(:) - -! Insert rEff (moist effective radius) - this%reff(:,:) = reff_table(:,:) - -! Now we linearly interpolate the input table to the output table grid -! of requested channels - do j = 1, this%nbin - do i = 1, this%nrh - do n = 1, this%nlambda - call polint(channels_table,bext_table(:,i,j),nch_table, & - this%lambda(n),this%bext(i,n,j),yerr) - call polint(channels_table,bsca_table(:,i,j),nch_table, & - this%lambda(n),this%bsca(i,n,j),yerr) - call polint(channels_table,bbck_table(:,i,j),nch_table, & - this%lambda(n),this%bbck(i,n,j),yerr) - call polint(channels_table,g_table(:,i,j),nch_table, & - this%lambda(n),this%g(i,n,j),yerr) - call polint(channels_table,refr_table(:,i,j),nch_table, & - this%lambda(n),this%refr(i,n,j),yerr) - call polint(channels_table,refi_table(:,i,j),nch_table, & - this%lambda(n),this%refi(i,n,j),yerr) - do ipol = 1, this%nPol - call polint(channels_table,pback_table(:,i,j,ipol),nch_table, & - this%lambda(n),this%pback(i,n,j,ipol),yerr) - end do - if ( nmom_ > 0 ) then - do imom = 1, this%nMom - do ipol = 1, this%nPol - call polint(channels_table,pmom_table(:,i,j,imom,ipol),nch_table, & - this%lambda(n),this%pmom(i,n,j,imom,ipol),yerr) - end do - end do - end if - enddo - enddo - enddo - -! Insert growth factor - this%gf(:,:) = gf_table(:,:) - -! Wet particle density [kg m-3] - this%rhop(:,:) = rhop_table(:,:) - -! Dry particle density [kg m-3] - this%rhod(:,:) = rhod_table(:,:) - -! Volume [m3 kg-1] - this%vol(:,:) = vol_table(:,:) - -! Area [m2 kg-1] - this%area(:,:) = area_table(:,:) - -! Now we do a mapping of the RH from the input table to some high -! resolution representation. This is to spare us the need to -! do a full-up interpolation later on. -! RH input from the table is scaled 0 - 0.99 -! We resolve the map to 0 - 0.990 in steps of 0.001 (991 total steps) - do j = 1, 991 - do i = this%nrh, 1, -1 - if( (j-1) .ge. int(this%rh(i)*1000)) then - ip1 = i + 1 - this%rhi(j) = i - if(ip1 .gt. this%nrh) then - this%rha(j) = 0. - else - this%rha(j) = ( (j-1)/1000. - this%rh(i)) & - / ( this%rh(ip1)- this%rh(i)) - endif - exit - endif - enddo -! print *, j, this%rhi(j), this%rha(j), this%rh(this%rhi(j)) - enddo - -! deallocate (channels_table, rh_table, bext_table, bsca_table, & -! bbck_table, g_table, stat = rc ) - - deallocate (channels_table, __NF_STAT__) - deallocate (rh_table, __NF_STAT__) - deallocate (reff_table, __NF_STAT__) - deallocate (bext_table, __NF_STAT__) - deallocate (bsca_table, __NF_STAT__) - deallocate (bbck_table, __NF_STAT__) - deallocate (g_table, __NF_STAT__) - deallocate (pback_table, __NF_STAT__) - if ( nmom_ > 0 ) then - deallocate (pmom_table, __NF_STAT__) - endif - deallocate (gf_table, __NF_STAT__) - deallocate (rhop_table, __NF_STAT__) - deallocate (rhod_table, __NF_STAT__) - deallocate (vol_table, __NF_STAT__) - deallocate (area_table, __NF_STAT__) - deallocate (refr_table, __NF_STAT__) - deallocate (refi_table, __NF_STAT__) - -return - -contains - - subroutine polint(x,y,n,xWant,yWant,yErr) - integer :: n -! recall, table hard-wired single precision - real*8 :: x(n),y(n) - real :: xWant, yWant, yErr - -! given array x(n) of independent variables and array y(n) of dependent -! variables, compute the linear interpolated result yWant at xWant and return -! with a dummy error estimate yErr. Hacked up from Numerical Recipes Chapter 3 - - integer :: i, j - real :: dx, slope - character(len=255) :: msg - -! on out of bounds, set i to lower or upper limit - i = 0 - if(xWant .lt. x(1)) then -! write(msg,*) "in polint, wanted: ", xWant, ", got lower bound: ", x(1) -! call warn(myname,msg) -! if (mapl_am_i_root()) print *,'in polint(), wnted: ', xWant, ', got lower bound: ', x(1) - i = 1 - endif - if(xWant .gt. x(n)) then -! write(msg,*) "in polint, wanted: ", xWant, ", got upper bound: ", x(n) -! call warn(myname,msg) -! if (mapl_am_i_root()) print *,'in polint(), wnted: ', xWant, ', got upper bound: ', x(n) - - i = n - endif - -! if i is still zero find i less than xWant - if(i .eq. 0) then - do j = 1, n - if(xWant .ge. x(j)) i = j - enddo - endif - -! slope - if(i .eq. n) then - slope = 0. - else - slope = (y(i+1)-y(i)) / (x(i+1)-x(i)) - endif - dx = xWant - x(i) - yWant = y(i) + slope*dx - - yErr = 0. - - return - end subroutine polint - -END SUBROUTINE Chem_MieTableRead -!---------------------------------------------------------------------------- - -!BOP -! -! !IROUTINE: Chem_MieQueryByInt --- Return Tau, SSA, etc -! -! -! !INTERFACE: -! - impure elemental subroutine Chem_MieQueryByInt ( this, idx, channel, q_mass, rh, & - tau, ssa, gasym, bext, bsca, bbck, & - reff, p11, p22, gf, rhop, rhod, & - vol, area, refr, refi, rc ) - -! !INPUT PARAMETERS: - - type(Chem_Mie), target, intent(in ) :: this - integer, intent(in ) :: idx ! variable index on Chem_Mie - real, intent(in ) :: channel ! channel number - real, intent(in ) :: q_mass ! aerosol mass [kg/m2], - real, intent(in ) :: rh ! relative himidity - -! !OUTPUT PARAMETERS: - - real, optional, intent(out) :: tau ! aerol extinction optical depth - real, optional, intent(out) :: ssa ! single scattering albedo - real, optional, intent(out) :: gasym ! asymmetry parameter - real, optional, intent(out) :: bext ! mass extinction efficiency [m2 (kg dry mass)-1] - real, optional, intent(out) :: bsca ! mass scattering efficiency [m2 (kg dry mass)-1] - real, optional, intent(out) :: bbck ! mass backscatter efficiency [m2 (kg dry mass)-1] - real, optional, intent(out) :: reff ! effective radius (micron) - real, optional, intent(out) :: p11 ! P11 phase function at backscatter - real, optional, intent(out) :: p22 ! P22 phase function at backscatter - real, optional, intent(out) :: gf ! Growth factor (ratio of wet to dry radius) - real, optional, intent(out) :: rhop ! Wet particle density [kg m-3] - real, optional, intent(out) :: rhod ! Dry particle density [kg m-3] - real, optional, intent(out) :: vol ! Wet particle volume [m3 kg-1] - real, optional, intent(out) :: area ! Wet particle cross section [m2 kg-1] - real, optional, intent(out) :: refr ! Wet particle real part of ref. index - real, optional, intent(out) :: refi ! Wet particle imag. part of ref. index - integer, optional, intent(out) :: rc ! error code - -! !DESCRIPTION: -! -! Returns requested parameters from the Mie tables, as a function -! of species, relative humidity, and channel -! -! Notes: Needs some checking, and I still force an interpolation step - -! -! !REVISION HISTORY: -! -! 23Mar2005 Colarco -! 11Jul2005 da Silva Standardization. -! -!EOP -!------------------------------------------------------------------------- - - - integer :: ICHANNEL, TYPE - integer :: irh, irhp1, isnap - real :: rhUse, arh - real :: bextIn, bscaIn, bbckIn, gasymIn, p11In, p22In, & - gfIn, rhopIn, rhodIn, volIn, areaIn, & - refrIn, refiIn - type(Chem_MieTable), pointer :: TABLE - - character(len=*), parameter :: Iam = 'Chem_MieQuery' - - if ( present(rc) ) rc = 0 - - ICHANNEL = nint(CHANNEL) -! TABLE => this%vtableUse - Table => this%mie_aerosol - TYPE = idx - -! ASSERT_(TYPE>0) -! ASSERT_(ICHANNEL>=LBOUND(TABLE%bext,1)) -! ASSERT_(ICHANNEL<=UBOUND(TABLE%bext,1)) - -! Now map the input RH to the high resolution hash table for RH - rhUse = max(rh,0.) - rhUse = min(rh,0.99) - isnap = int((rhUse+0.001)*1000.) - if(isnap .lt. 1) isnap = 1 - arh = TABLE%rha( isnap ) - irh = TABLE%rhi( isnap ) - irhp1 = irh+1 - if(irhp1 .gt. TABLE%nrh) irhp1 = TABLE%nrh - -! Now linearly interpolate the input table for the requested aerosol and -! channel; rh is the relative humidity. - - if(present(bext) .or. present(tau) .or. present(ssa) ) then - bextIn = TABLE%bext(irh ,ichannel,TYPE) * (1.-arh) & - + TABLE%bext(irhp1,ichannel,TYPE) * arh - endif - - if(present(bsca) .or. present(ssa) ) then - bscaIn = TABLE%bsca(irh ,ichannel,TYPE) * (1.-arh) & - + TABLE%bsca(irhp1,ichannel,TYPE) * arh - endif - - if(present(bbck)) then - bbckIn = TABLE%bbck(irh ,ichannel,TYPE) * (1.-arh) & - + TABLE%bbck(irhp1,ichannel,TYPE) * arh - endif - - if(present(gasym)) then - gasymIn = TABLE%g(irh ,ichannel,TYPE) * (1.-arh) & - + TABLE%g(irhp1,ichannel,TYPE) * arh - endif - - if(present(rEff) ) then - rEff = TABLE%rEff(irh ,TYPE) * (1.-arh) & - + TABLE%rEff(irhp1,TYPE) * arh - rEff = 1.E6 * rEff ! convert to microns - endif - -! if(present(pmom)) then -! pmom(:,:) = TABLE%pmom(irh ,ichannel,TYPE,:,:) * (1.-arh) & -! + TABLE%pmom(irhp1,ichannel,TYPE,:,:) * arh -! endif - -! if(present(pmom)) then -! call Chem_MieQueryByIntWithpmom(this, idx, channel, q_mass, rh, pmom) -! endif - - - if(present(p11) ) then - p11In = TABLE%pback(irh ,ichannel,TYPE,1) * (1.-arh) & - + TABLE%pback(irhp1,ichannel,TYPE,1) * arh - endif - - if(present(p22) ) then - p22In = TABLE%pback(irh ,ichannel,TYPE,5) * (1.-arh) & - + TABLE%pback(irhp1,ichannel,TYPE,5) * arh - endif - - if(present(gf) ) then - gfIn = TABLE%gf(irh ,TYPE) * (1.-arh) & - + TABLE%gf(irhp1,TYPE) * arh - endif - - if(present(rhod) ) then - rhodIn = TABLE%rhod(1 ,TYPE) - endif - - if(present(vol) ) then - volIn = TABLE%vol(irh ,TYPE) * (1.-arh) & - + TABLE%vol(irhp1,TYPE) * arh - endif - - if(present(area) ) then - areaIn = TABLE%area(irh ,TYPE) * (1.-arh) & - + TABLE%area(irhp1,TYPE) * arh - endif - - if(present(refr) .or. present(tau) .or. present(ssa) ) then - refrIn = TABLE%refr(irh ,ichannel,TYPE) * (1.-arh) & - + TABLE%refr(irhp1,ichannel,TYPE) * arh - endif - - if(present(refi) .or. present(tau) .or. present(ssa) ) then - refiIn = TABLE%refi(irh ,ichannel,TYPE) * (1.-arh) & - + TABLE%refi(irhp1,ichannel,TYPE) * arh - endif - -! Fill the requested outputs - if(present(tau )) tau = bextIn * q_mass - if(present(ssa )) ssa = bscaIn/bextIn - if(present(bext )) bext = bextIn - if(present(bsca )) bsca = bscaIn - if(present(bbck )) bbck = bbckIn - if(present(gasym)) gasym = gasymIn - if(present(p11 )) p11 = p11In - if(present(p22 )) p22 = p22In - if(present(gf )) gf = gfIn - if(present(rhop )) rhop = rhopIn - if(present(rhod )) rhod = rhodIn - if(present(vol )) vol = volIn - if(present(area )) area = areaIn - if(present(refr )) refr = refrIn - if(present(refi )) refi = refiIn - -! All Done -!---------- - end subroutine Chem_MieQueryByInt - - -!------------------------------------------------------------------------- -!BOP -! -! !IROUTINE: Chem_MieQueryByIntWithpmom --- Return Tau, SSA, etc -! -! -! !INTERFACE: -! - subroutine Chem_MieQueryByIntWithpmom ( this, idx, channel, q_mass, rh, & - tau, ssa, gasym, bext, bsca, bbck, & - reff, pmom, p11, p22, gf, rhop, rhod, & - vol, area, refr, refi, rc ) - -! !INPUT PARAMETERS: - - type(Chem_Mie), target, intent(in ) :: this - integer, intent(in ) :: idx ! variable index on Chem_Mie - real, intent(in ) :: channel ! channel number - real, intent(in ) :: q_mass ! aerosol mass [kg/m2], - real, intent(in ) :: rh ! relative himidity - -! !OUTPUT PARAMETERS: - - real, optional, intent(out) :: tau ! aerol extinction optical depth - real, optional, intent(out) :: ssa ! single scattering albedo - real, optional, intent(out) :: gasym ! asymmetry parameter - real, optional, intent(out) :: bext ! mass extinction efficiency [m2 (kg dry mass)-1] - real, optional, intent(out) :: bsca ! mass scattering efficiency [m2 (kg dry mass)-1] - real, optional, intent(out) :: bbck ! mass backscatter efficiency [m2 (kg dry mass)-1] - real, optional, intent(out) :: reff ! effective radius (micron) - real, intent(out) :: pmom(:,:) - real, optional, intent(out) :: p11 ! P11 phase function at backscatter - real, optional, intent(out) :: p22 ! P22 phase function at backscatter - real, optional, intent(out) :: gf ! Growth factor (ratio of wet to dry radius) - real, optional, intent(out) :: rhop ! Wet particle density [kg m-3] - real, optional, intent(out) :: rhod ! Dry particle density [kg m-3] - real, optional, intent(out) :: vol ! Wet particle volume [m3 kg-1] - real, optional, intent(out) :: area ! Wet particle cross section [m2 kg-1] - real, optional, intent(out) :: refr ! Wet particle real part of ref. index - real, optional, intent(out) :: refi ! Wet particle imag. part of ref. index - integer, optional, intent(out) :: rc ! error code - -! !DESCRIPTION: -! -! Returns requested parameters from the Mie tables, as a function -! of species, relative humidity, and channel -! -! Notes: Needs some checking, and I still force an interpolation step - -! -! !REVISION HISTORY: -! -! 23Mar2005 Colarco -! 11Jul2005 da Silva Standardization. -! -!EOP -!------------------------------------------------------------------------- - - - integer :: ICHANNEL, TYPE - integer :: irh, irhp1, isnap - real :: rhUse, arh - type(Chem_MieTable), pointer :: TABLE - - character(len=*), parameter :: Iam = 'Chem_MieQueryByIntWithpmom' - - if ( present(rc) ) rc = 0 - - ICHANNEL = nint(CHANNEL) - TABLE => this%vtableUse - TYPE = idx - -! Now map the input RH to the high resolution hash table for RH - rhUse = max(rh,0.) - rhUse = min(rh,0.99) - isnap = int((rhUse+0.001)*1000.) - if(isnap .lt. 1) isnap = 1 - arh = TABLE%rha( isnap ) - irh = TABLE%rhi( isnap ) - irhp1 = irh+1 - if(irhp1 .gt. TABLE%nrh) irhp1 = TABLE%nrh - -! Now linearly interpolate the input table for the requested aerosol and -! channel; rh is the relative humidity. - call Chem_MieQuery ( this, idx, channel, q_mass, rh, & - tau, ssa, gasym, bext, bsca, bbck, & - reff, p11, p22, gf, rhop, rhod, & - vol, area, refr, refi, rc ) - NF_VERIFY_(rc) - - pmom(:,:) = TABLE%pmom(irh ,ichannel,TYPE,:,:) * (1.-arh) & - + TABLE%pmom(irhp1,ichannel,TYPE,:,:) * arh - - - -! All Done -!---------- - end subroutine Chem_MieQueryByIntWithpmom - - - - - end module Chem_MieTableMod2G - diff --git a/Process_Library/GOCART2G_MieMod.F90 b/Process_Library/GOCART2G_MieMod.F90 new file mode 100644 index 00000000..4f00dfca --- /dev/null +++ b/Process_Library/GOCART2G_MieMod.F90 @@ -0,0 +1,571 @@ +!BOP +! +! !MODULE: GOCART2G_MieMod --- Reader for aerosol mie tables +! +! !INTERFACE: +! + +module GOCART2G_MieMod +! !USES: + use netcdf + implicit none + +! !PUBLIC TYPES: +! + private + public GOCART2G_Mie ! Holds Chem_MieTable and other information + +! +! !PUBLIC MEMBER FUNCTIONS: +! +! +! !DESCRIPTION: +! +! This module read the mie aerosol tables. +! +! !REVISION HISTORY: +! +! 23Mar2005 Colarco - Initial code. +! 31Mar2005 Todling - Declared netcdf nf_ routines as external (OSF1) +! Removed # from include netcdf.inc +! 16Feb2022 da Silva Refactoring from old Chem_MieTable and Chem_MieMod; +! these functionalities have been merged. +! +!EOP +!------------------------------------------------------------------------- +! +! Mie LUT table +! Will be reduced from input files to the desired channels +! -------- + integer, parameter :: NRH_BINS = 991 + + type GOCART2G_Mie + + private + character(len=:), allocatable :: table_name + integer :: nch ! number of channels in table (replacement of nlamfda) + integer :: nrh ! number of RH values in table + integer :: nbin ! number of size bins in table + integer :: nMom ! number of moments of phase function + integer :: nPol ! number of elements of scattering phase matrix + + ! c=channel, r=rh, b=bin, m=moments, p=nPol + real, allocatable :: wavelengths(:) ! (c) wavelengths [m] + real, allocatable :: rh(:) ! (r) RH values [fraction] + real, allocatable :: reff(:,:) ! (r,b) effective radius [m] + real, allocatable :: bext(:,:,:) ! (r,c,b) bext values [m2 kg-1] + real, allocatable :: bsca(:,:,:) ! (r,c,b) bsca values [m2 kg-1] + real, allocatable :: bbck(:,:,:) ! (r,c,b) bbck values [m2 kg-1] + real, allocatable :: g(:,:,:) ! (r,c,b) asymmetry parameter +!ams real, allocatable :: pback(:,:,:,:) ! (r,c,b,p) Backscatter phase function + real, allocatable :: p11(:,:,:) ! (r,c,b) Backscatter phase function, index 1 + real, allocatable :: p22(:,:,:) ! (r,c,b) Backscatter phase function, index 5 + real, allocatable :: pmom(:,:,:,:,:) ! (r,c,b,m,p) moments of phase function + real, allocatable :: gf(:,:) ! (r,b) hygroscopic growth factor + real, allocatable :: rhop(:,:) ! (r,b) wet particle density [kg m-3] + real, allocatable :: rhod(:,:) ! (r,b) wet particle density [kg m-3] + real, allocatable :: vol(:,:) ! (r,b) wet particle volume [m3 kg-1] + real, allocatable :: area(:,:) ! (r,b) wet particle cross section [m2 kg-1] + real, allocatable :: refr(:,:,:) ! (r,c,b) real part of refractive index + real, allocatable :: refi(:,:,:) ! (r,c,b) imaginary part of refractive index + + integer :: rhi(NRH_BINS) ! pointer to rh LUT + real :: rha(NRH_BINS) ! slope on rh LUT + + CONTAINS + + procedure :: QueryByWavelength_1d + procedure :: QueryByWavelength_2d + procedure :: QueryByWavelength_3d + procedure :: QueryByChannel_1d + procedure :: QueryByChannel_2d + procedure :: QueryByChannel_3d + generic :: Query => QueryByWavelength_1d, & + QueryByWavelength_2d, & + QueryByWavelength_3d, & + QueryByChannel_1d, & + QueryByChannel_2d, & + QueryByChannel_3d + procedure :: getChannel + procedure :: getWavelength + + end type GOCART2G_Mie + + interface GOCART2G_Mie + module procedure GOCART2G_MieCreate + end interface GOCART2G_Mie + +CONTAINS + +!------------------------------------------------------------------------- +! NASA/GSFC, Global Modeling and Assimilation Office, Code 900.3 ! +!------------------------------------------------------------------------- +!BOP +! +! !IROUTINE: GOCART2G_MieCreate --- Construct Chemistry Registry +! +! +! !DESCRIPTION: +! +! +! !REVISION HISTORY: +! +! 09Mar2005 da Silva API, prologues. +! + + type(GOCART2G_Mie) function GOCART2G_MieCreate ( rcfile, wavelengths, nmom, rc ) result (this) + +! !INPUT PARAMETERS: + + character(len=*), intent(in) :: rcfile ! Mie table file name + real, optional, intent(in) :: wavelengths(:) + integer, optional,intent(in) :: nmom + +! !OUTPUT PARAMETERS: + + integer, intent(out) :: rc ! Error return code: + ! 0 - all is well +! !DESCRIPTION: +! +! Fills in the Mie table +! +! !REVISION HISTORY: +! +! 23Mar2005 Colarco +! +!EOP +!------------------------------------------------------------------------- + character(len=*), parameter :: myname = 'GOCART2G_MieCreate' + integer :: nch + integer :: ncid, idimid, ivarid, n, i, j, ip1 + integer :: nch_table, nrh_table, nbin_table, nmom_table, nPol_table +! Tables are hard-wired as single precision + real, allocatable :: channels_table(:), rh_table(:), reff_table(:,:), & + bext_table(:,:,:), bsca_table(:,:,:), & + bbck_table(:,:,:), g_table(:,:,:), & + pmom_table(:,:,:,:,:),pback_table(:,:,:,:), & + gf_table(:,:), rhop_table(:,:), rhod_table(:,:),& + vol_table(:,:), area_table(:,:), & + refr_table(:,:,:), refi_table(:,:,:) + + real, allocatable :: pback(:,:,:,:) ! (r,c,b,p) Backscatter phase function + + real :: yerr + integer :: nmom_, imom, ipol + integer :: status + +#define NF_VERIFY_(expr) rc = expr; if (rc /= 0) return +#define __NF_STAT__ stat=status); NF_VERIFY_(status + + rc = 0 + this%table_name = rcfile + + +! Whether or not doing phase function +! ----------------------------------- + if ( present(nmom) ) then + nmom_ = nmom + else + nmom_ = 0 + end if + +! Set up nPol_table for reading backscatter phase function +! This will get overwritten if pmoments is requested +! -------------------------------------------------------- + nPol_table = 6 + +! Open the table and get the dimensions +! ------------------------------------- + rc = nf90_open(this%table_name, NF90_NOWRITE, ncid) + IF ( rc /= 0 ) THEN + print *, 'nf90_open '//this%table_name//' RETURN CODE=', rc + NF_VERIFY_(rc) + END IF + +! RH +! -- + NF_VERIFY_(nf90_inq_dimid(ncid,'rh',idimid)) + NF_VERIFY_(nf90_inquire_dimension(ncid,idimid,len=nrh_table)) + +! Channels +! -------- + NF_VERIFY_(nf90_inq_dimid(ncid,'lambda',idimid)) + NF_VERIFY_(nf90_inquire_dimension(ncid,idimid,len=nch_table)) + + if (present(wavelengths) ) then + nch = size(wavelengths) + else + nch = nch_table + end if + +! Dry Effective radius +! -------------------- + NF_VERIFY_(nf90_inq_dimid(ncid,'radius',idimid)) + NF_VERIFY_(nf90_inquire_dimension(ncid,idimid,len=nbin_table)) + +! Moments of phase function +! ------------------------- + if ( nmom_ > 0 ) then + NF_VERIFY_(nf90_inq_dimid(ncid,'nMom',idimid)) + NF_VERIFY_(nf90_inquire_dimension(ncid,idimid,len=nmom_table)) + if ( nmom_ > nmom_table ) then +! rc = 99 + print*,'Error: nmom_ > nmom_table, see:'//myname + NF_VERIFY_(1) + end if + NF_VERIFY_(nf90_inq_dimid(ncid,'nPol',idimid)) + NF_VERIFY_(nf90_inquire_dimension(ncid,idimid,len=nPol_table)) + endif + +! Get the table contents +! ------------------------------------- + + allocate(channels_table(nch_table), __NF_STAT__) + allocate(rh_table(nrh_table), __NF_STAT__) + allocate(reff_table(nrh_table,nbin_table), __NF_STAT__) + allocate(bext_table(nch_table,nrh_table,nbin_table), __NF_STAT__) + allocate(bsca_table(nch_table,nrh_table,nbin_table), __NF_STAT__) + allocate(bbck_table(nch_table,nrh_table,nbin_table), __NF_STAT__) + allocate(g_table(nch_table,nrh_table,nbin_table), stat = rc ) + allocate(pback_table(nch_table,nrh_table,nbin_table,nPol_table), __NF_STAT__) + allocate(gf_table(nrh_table,nbin_table), __NF_STAT__) + allocate(rhop_table(nrh_table,nbin_table), __NF_STAT__) + allocate(rhod_table(nrh_table,nbin_table), __NF_STAT__) + allocate(vol_table(nrh_table,nbin_table), __NF_STAT__) + allocate(area_table(nrh_table,nbin_table), __NF_STAT__) + allocate(refr_table(nch_table,nrh_table,nbin_table), __NF_STAT__) + allocate(refi_table(nch_table,nrh_table,nbin_table), __NF_STAT__) + + if ( nmom_ > 0 ) then + allocate(pmom_table(nch_table,nrh_table,nbin_table,nmom_table,nPol_table), __NF_STAT__) + end if + NF_VERIFY_(nf90_inq_varid(ncid,'lambda',ivarid)) + NF_VERIFY_(nf90_get_var(ncid,ivarid,channels_table)) + NF_VERIFY_(nf90_inq_varid(ncid,'rEff',ivarid)) + NF_VERIFY_(nf90_get_var(ncid,ivarid,reff_table)) + NF_VERIFY_(nf90_inq_varid(ncid,'bext',ivarid)) + NF_VERIFY_(nf90_get_var(ncid,ivarid,bext_table)) + NF_VERIFY_(nf90_inq_varid(ncid,'bsca',ivarid)) + NF_VERIFY_(nf90_get_var(ncid,ivarid,bsca_table)) + NF_VERIFY_(nf90_inq_varid(ncid,'bbck',ivarid)) + NF_VERIFY_(nf90_get_var(ncid,ivarid,bbck_table)) + NF_VERIFY_(nf90_inq_varid(ncid,'g',ivarid)) + NF_VERIFY_(nf90_get_var(ncid,ivarid,g_table)) + NF_VERIFY_(nf90_inq_varid(ncid,'rh',ivarid)) + NF_VERIFY_(nf90_get_var(ncid,ivarid,rh_table)) + + ! TODO: we need to look at these NF90_NOERR checks +! Get the backscatter phase function values + rc = nf90_inq_varid(ncid,'pback',ivarid) + if(rc .ne. NF90_NOERR) then ! pback not in table, fill in dummy variable + pback_table = 1. + else + NF_VERIFY_(nf90_get_var(ncid,ivarid,pback_table)) + endif + + if ( nmom_ > 0 ) then + NF_VERIFY_(nf90_inq_varid(ncid,'pmom',ivarid)) + NF_VERIFY_(nf90_get_var(ncid,ivarid,pmom_table)) + end if + +! Aerosol optical properties not necessarily stored in all versions of the tables +! ---------------------- +! Particle growth factor + rc = nf90_inq_varid(ncid,'growth_factor',ivarid) + if(rc .ne. NF90_NOERR) then ! not in table, fill in dummy variable + gf_table = -999. + else + NF_VERIFY_(nf90_get_var(ncid,ivarid,gf_table)) + endif + +! Wet particle density + rc = nf90_inq_varid(ncid,'rhop',ivarid) + if(rc .ne. NF90_NOERR) then ! not in table, fill in dummy variable + rhop_table = -999. + else + NF_VERIFY_(nf90_get_var(ncid,ivarid,rhop_table)) + endif + +! Dry particle density (will be pulled from wet particle radius) + rc = nf90_inq_varid(ncid,'rhop',ivarid) + if(rc .ne. NF90_NOERR) then ! not in table, fill in dummy variable + rhod_table = -999. + else + rc = nf90_get_var(ncid,ivarid,rhod_table) + do i = 1, nrh_table + rhod_table(i,:) = rhod_table(1,:) + enddo + if (rc /=0) return + endif + +! Wet particle real part of refractive index + rc = nf90_inq_varid(ncid,'refreal',ivarid) + if(rc .ne. NF90_NOERR) then ! not in table, fill in dummy variable + refr_table = -999. + else + NF_VERIFY_(nf90_get_var(ncid,ivarid,refr_table)) + endif + +! Wet particle imaginary part of refractive index (ensure positive) + rc = nf90_inq_varid(ncid,'refimag',ivarid) + if(rc .ne. NF90_NOERR) then ! not in table, fill in dummy variable + refi_table = -999. + else + NF_VERIFY_(nf90_get_var(ncid,ivarid,refi_table)) + refi_table = abs(refi_table) + endif + +! Wet particle volume [m3 kg-1] +! Ratio of wet to dry volume is gf^3, hence the following + vol_table = gf_table**3 / rhod_table + +! Wet particle cross sectional area [m2 kg-1] +! Assume area is volume divided by (4./3.*reff) + area_table = vol_table / (4./3.*reff_table) + +! Close the table file +! ------------------------------------- + NF_VERIFY_(nf90_close(ncid)) + + +! Setup the table to be returned +! ------------------------------------- + this%nch = nch + this%nrh = nrh_table + this%nbin = nbin_table + this%nMom = nmom_ + this%nPol = nPol_table + + allocate (this%bext(this%nrh,this%nch,this%nbin), __NF_STAT__) + allocate (this%bsca(this%nrh,this%nch,this%nbin), __NF_STAT__) + allocate (this%bbck(this%nrh,this%nch,this%nbin), __NF_STAT__) + allocate (this%g(this%nrh,this%nch,this%nbin), __NF_STAT__) + allocate (pback(this%nrh,this%nch,this%nbin,this%nPol), __NF_STAT__) + if ( nmom_ > 0 ) then + allocate (this%pmom(this%nrh,this%nch,this%nbin,this%nMom,this%nPol), __NF_STAT__) + end if + allocate (this%refr(this%nrh,this%nch,this%nbin), __NF_STAT__) + allocate (this%refi(this%nrh,this%nch,this%nbin), __NF_STAT__) + +! Preserve the full RH structure of the input table + this%rh = rh_table ! assignment does allocation + +! Insert rEff (moist effective radius) + this%reff = reff_table + +! Insert growth factor + this%gf = gf_table + +! Wet particle density [kg m-3] + this%rhop = rhop_table + +! Dry particle density [kg m-3] + this%rhod = rhod_table + +! Volume [m3 kg-1] + this%vol = vol_table + +! Area [m2 kg-1] + this%area = area_table + +! Insert the requested channels in the output table + if ( present(wavelengths) ) then + this%wavelengths = wavelengths + else + this%wavelengths = channels_table + endif + +! Now we linearly interpolate the input table to the output table grid +! of requested channels + if ( present(wavelengths) ) then + do j = 1, this%nbin + do i = 1, this%nrh + do n = 1, this%nch + call polint(channels_table,bext_table(:,i,j),nch_table, & + this%wavelengths(n),this%bext(i,n,j),yerr) + call polint(channels_table,bsca_table(:,i,j),nch_table, & + this%wavelengths(n),this%bsca(i,n,j),yerr) + call polint(channels_table,bbck_table(:,i,j),nch_table, & + this%wavelengths(n),this%bbck(i,n,j),yerr) + call polint(channels_table,g_table(:,i,j),nch_table, & + this%wavelengths(n),this%g(i,n,j),yerr) + call polint(channels_table,refr_table(:,i,j),nch_table, & + this%wavelengths(n),this%refr(i,n,j),yerr) + call polint(channels_table,refi_table(:,i,j),nch_table, & + this%wavelengths(n),this%refi(i,n,j),yerr) + do ipol = 1, this%nPol + call polint(channels_table,pback_table(:,i,j,ipol),nch_table, & + this%wavelengths(n),pback(i,n,j,ipol),yerr) + end do + if ( nmom_ > 0 ) then + do imom = 1, this%nMom + do ipol = 1, this%nPol + call polint(channels_table,pmom_table(:,i,j,imom,ipol),nch_table, & + this%wavelengths(n),this%pmom(i,n,j,imom,ipol),yerr) + enddo + enddo + endif + enddo + enddo + enddo + else !(no wavelength) + !swap the order + this%bext = reshape(bext_table, [nrh_table, nch, nbin_table],order =[2,1,3]) + this%bsca = reshape(bsca_table, [nrh_table, nch, nbin_table],order =[2,1,3]) + this%bbck = reshape(bbck_table, [nrh_table, nch, nbin_table],order =[2,1,3]) + this%g = reshape( g_table, [nrh_table, nch, nbin_table],order =[2,1,3]) + this%refr = reshape(refr_table, [nrh_table, nch, nbin_table],order =[2,1,3]) + this%refi = reshape(refi_table, [nrh_table, nch, nbin_table],order =[2,1,3]) + pback = reshape(pback_table,[nrh_table, nch, nbin_table, npol_table],order =[2,1,3,4]) + if ( nmom_ > 0 ) then + this%pmom = reshape(pmom_table,[nrh_table,nch, nbin_table, nmom_, npol_table], order = [2,1,3,4,5]) + endif + endif + +! Pick p11, p12 + this%p11 = pback(:,:,:,1) + this%p22 = pback(:,:,:,5) + +! Now we do a mapping of the RH from the input table to some high +! resolution representation. This is to spare us the need to +! do a full-up interpolation later on. +! RH input from the table is scaled 0 - 0.99 +! We resolve the map to 0 - 0.990 in steps of 0.001 (991 total steps) + do j = 1, NRH_BINS + do i = this%nrh, 1, -1 + if ((j-1) .ge. int(this%rh(i)*1000)) then + ip1 = i + 1 + this%rhi(j) = i + if (ip1 .gt. this%nrh) then + this%rha(j) = 0. + else + this%rha(j) = ( (j-1)/1000. - this%rh(i)) & + / ( this%rh(ip1)- this%rh(i)) + endif + exit + endif + enddo +! print *, j, this%rhi(j), this%rha(j), this%rh(this%rhi(j)) + enddo + + return + + contains + + subroutine polint(x,y,n,xWant,yWant,yErr) + integer :: n + ! recall, table hard-wired single precision + real :: x(n),y(n) + real :: xWant, yWant, yErr + + ! given array x(n) of independent variables and array y(n) of dependent + ! variables, compute the linear interpolated result yWant at xWant and return + ! with a dummy error estimate yErr. Hacked up from Numerical Recipes Chapter 3 + + integer :: i, j + real :: dx, slope + character(len=255) :: msg + + ! on out of bounds, set i to lower or upper limit + i = 0 + if(xWant .lt. x(1)) then + ! write(msg,*) "in polint, wanted: ", xWant, ", got lower bound: ", x(1) + ! call warn(myname,msg) + ! if (mapl_am_i_root()) print *,'in polint(), wnted: ', xWant, ', got lower bound: ', x(1) + i = 1 + endif + if(xWant .gt. x(n)) then + ! write(msg,*) "in polint, wanted: ", xWant, ", got upper bound: ", x(n) + ! call warn(myname,msg) + ! if (mapl_am_i_root()) print *,'in polint(), wnted: ', xWant, ', got upper bound: ', x(n) + + i = n + endif + + ! if i is still zero find i less than xWant + if (i .eq. 0) then + do j = 1, n + if(xWant .ge. x(j)) i = j + enddo + endif + + ! slope + if (i .eq. n) then + slope = 0. + else + slope = (y(i+1)-y(i)) / (x(i+1)-x(i)) + endif + dx = xWant - x(i) + yWant = y(i) + slope*dx + + yErr = 0. + + end subroutine polint + + end function GOCART2G_MieCreate + +! +! Query subroutines +! + +#define RANK_ 1 +#include "MieQuery.H" +#undef RANK_ + +#define RANK_ 2 +#include "MieQuery.H" +#undef RANK_ + +#define RANK_ 3 +#include "MieQuery.H" +#undef RANK_ + + integer function getChannel(this, wavelength, rc) result (ch) + class (GOCART2G_Mie), intent(in) :: this + real, intent(in) :: wavelength + integer, optional, intent(out) :: rc + real, parameter :: w_tol = 1.e-9 + integer :: i + + ch = -1 + do i = 1, this%nch + if (abs(this%wavelengths(i)-wavelength) <= w_tol) then + ch = i + exit + endif + enddo + + if (present(rc)) rc = 0 + + if (ch < 0) then + !$omp critical + print*, "wavelength of ",wavelength, " is an invalid value." + !$omp end critical + if (present(rc)) rc = -1 + endif + + end function getChannel + + real function getWavelength(this, ith_channel, rc) result (wavelength) + class (GOCART2G_Mie), intent(in) :: this + integer, intent(in) :: ith_channel + integer, optional, intent(out) :: rc + real, parameter :: w_tol = 1.e-9 + integer :: i + + if (present(rc)) rc = 0 + + if (ith_channel <=0 .or. ith_channel > this%nch ) then + !$omp critical + print*, "The channel of ",ith_channel, " is an invalid channel number." + !$omp end critical + if (present(rc)) rc = -1 + wavelength = -1. ! meanlingless nagative + return + endif + + wavelength = this%wavelengths(ith_channel) + + end function getWavelength + +end module GOCART2G_MieMod diff --git a/Process_Library/GOCART2G_Process.F90 b/Process_Library/GOCART2G_Process.F90 index 71dc8f51..c5ccfcc0 100644 --- a/Process_Library/GOCART2G_Process.F90 +++ b/Process_Library/GOCART2G_Process.F90 @@ -17,7 +17,7 @@ module GOCART2G_Process ! !USES: ! Only instrinsic fortran types and functions are allowed. - use Chem_MieTableMod2G + use GOCART2G_MieMod use, intrinsic :: iso_fortran_env, only: IOSTAT_END implicit none @@ -3249,7 +3249,7 @@ end subroutine UpdateAerosolState ! !IROUTINE: Aero_Compute_Diags - Calculate aerosol diagnostics ! ! !INTERFACE: - subroutine Aero_Compute_Diags (mie_table, km, klid, nbegin, nbins, rlow, rup, channels, & + subroutine Aero_Compute_Diags (mie, km, klid, nbegin, nbins, rlow, rup, & wavelengths_profile, wavelengths_vertint, aerosol, & grav, tmpu, rhoa, rh, u, v, delp, ple,tropp, & sfcmass, colmass, mass, exttau, stexttau, scatau, stscatau,& @@ -3262,12 +3262,11 @@ subroutine Aero_Compute_Diags (mie_table, km, klid, nbegin, nbins, rlow, rup, ch implicit NONE ! !INPUT PARAMETERS: - type(Chem_Mie), intent(in) :: mie_table ! mie table + type(GOCART2G_Mie), intent(in) :: mie ! mie table integer, intent(in) :: km, nbegin, nbins integer, intent(in) :: klid ! index for pressure lid real, optional, dimension(:), intent(in) :: rlow ! bin radii - low bounds real, optional, dimension(:), intent(in) :: rup ! bin radii - upper bounds - real, dimension(:), intent(in) :: channels real, dimension(:), intent(in) :: wavelengths_profile real, dimension(:), intent(in) :: wavelengths_vertint real, dimension(:,:,:,:), intent(in) :: aerosol ! @@ -3320,11 +3319,10 @@ subroutine Aero_Compute_Diags (mie_table, km, klid, nbegin, nbins, rlow, rup, ch ! !Local Variables character(len=*), parameter :: myname = 'Aero_Compute_Diags' - integer :: i, j, k, n, w, ios, nch + integer :: i, j, k, n, w, ios, status integer :: i1 =1, i2, j1=1, j2 - real :: ilam550, ilam470, ilam870 - real, allocatable, dimension(:) :: wavelengths_index_profile, wavelengths_index_vertint - real :: tau, ssa + integer :: ilam470, ilam870 + real, allocatable, dimension(:,:,:) :: tau, ssa ! real :: fPMfm(nbins) ! fraction of bin with particles diameter < 1.0 um ! real :: fPM25(nbins) ! fraction of bin with particles diameter < 2.5 um real, dimension(:), allocatable :: fPMfm ! fraction of bin with particles diameter < 1.0 um @@ -3345,7 +3343,6 @@ subroutine Aero_Compute_Diags (mie_table, km, klid, nbegin, nbins, rlow, rup, ch ! Initialize local variables ! -------------------------- - nch = size(channels) i2 = size(rhoa,1) j2 = size(rhoa,2) allocate(fPMfm(nbins),source=0.0) @@ -3353,64 +3350,20 @@ subroutine Aero_Compute_Diags (mie_table, km, klid, nbegin, nbins, rlow, rup, ch ! Get the wavelength indices ! -------------------------- -! Must provide ilam550 for AOT calculation - allocate(wavelengths_index_profile(size(wavelengths_profile))) - allocate(wavelengths_index_vertint(size(wavelengths_vertint))) - wavelengths_index_profile = 0. - wavelengths_index_vertint = 0. - ilam550 = 1. - ilam470 = 0. - ilam870 = 0. - if(nch .gt. 1) then - do i = 1, nch - if ( channels(i) .ge. 5.49e-7 .and. & - channels(i) .le. 5.51e-7) ilam550 = i - if ( channels(i) .ge. 4.69e-7 .and. & - channels(i) .le. 4.71e-7) ilam470 = i - if ( channels(i) .ge. 8.69e-7 .and. & - channels(i) .le. 8.71e-7) ilam870 = i - enddo - endif - ! Channel values are 4.7e-7 5.5e-7 6.7e-7 8.7e-7 [meter]. Their indices are 1,2,3,4 respectively. - do i = 1, size(wavelengths_profile) - if ((wavelengths_profile(i) .ge. 5.49e-7) .and. (wavelengths_profile(i) .le. 5.51e-7)) then - wavelengths_index_profile(i) = 2. - else if ((wavelengths_profile(i) .ge. 4.69e-7) .and. (wavelengths_profile(i) .le. 4.71e-7)) then - wavelengths_index_profile(i) = 1. - else if ((wavelengths_profile(i) .ge. 6.69e-7) .and. (wavelengths_profile(i) .le. 6.71e-7)) then - wavelengths_index_profile(i) = 3. - else if ((wavelengths_profile(i) .ge. 8.68e-7) .and. (wavelengths_profile(i) .le. 8.71e-7)) then - wavelengths_index_profile(i) = 4. - else - print*,'wavelengths_profile of ',wavelengths_profile(i),' is an invalid value.' - return - end if - end do + ilam470 = mie%getChannel(4.70e-7) + if(ilam470 <= 0) ilam470 = 0 - ! Channel values are 4.7e-7 5.5e-7 6.7e-7 8.7e-7 [meter]. Their indices are 1,2,3,4 respectively. - do i = 1, size(wavelengths_vertint) - if ((wavelengths_vertint(i) .ge. 5.49e-7) .and. (wavelengths_vertint(i) .le. 5.51e-7)) then - wavelengths_index_vertint(i) = 2. - else if ((wavelengths_vertint(i) .ge. 4.69e-7) .and. (wavelengths_vertint(i) .le. 4.71e-7)) then - wavelengths_index_vertint(i) = 1. - else if ((wavelengths_vertint(i) .ge. 6.69e-7) .and. (wavelengths_vertint(i) .le. 6.71e-7)) then - wavelengths_index_vertint(i) = 3. - else if ((wavelengths_vertint(i) .ge. 8.68e-7) .and. (wavelengths_vertint(i) .le. 8.71e-7)) then - wavelengths_index_vertint(i) = 4. - else - print*,'wavelengths_vertint of ',wavelengths_vertint(i),' is an invalid value.' - return - end if - end do + ilam870 = mie%getChannel(8.70e-7) + if(ilam870 <= 0) ilam870 = 0 ! Determine if going to do Angstrom parameter calculation ! ------------------------------------------------------- do_angstrom = .false. ! If both 470 and 870 channels provided (and not the same) then ! possibly will do Angstrom parameter calculation - if(ilam470 .ne. 0. .and. & - ilam870 .ne. 0. .and. & + if(ilam470 .ne. 0 .and. & + ilam870 .ne. 0 .and. & ilam470 .ne. ilam870) do_angstrom = .true. if( present(angstrom) .and. do_angstrom ) then @@ -3520,6 +3473,8 @@ subroutine Aero_Compute_Diags (mie_table, km, klid, nbegin, nbins, rlow, rup, ch end do endif + allocate(tau(i1:i2,j1:j2,km),source = 0.) + allocate(ssa(i1:i2,j1:j2,km),source = 0.) ! Calculate the extinction and/or scattering AOD if( present(extcoef) .or. & present(scacoef) ) then @@ -3528,27 +3483,19 @@ subroutine Aero_Compute_Diags (mie_table, km, klid, nbegin, nbins, rlow, rup, ch if( present(scacoef) ) scacoef = 0. do n = nbegin, nbins - do w = 1, size(wavelengths_profile) - do k = klid, km - do j = j1, j2 - do i = i1, i2 -! call Chem_MieQuery(mie_table, n, ilam550, & - call Chem_MieQuery(mie_table, n, wavelengths_index_profile(w), & - aerosol(i,j,k,n)*delp(i,j,k)/grav, & - rh(i,j,k), tau=tau, ssa=ssa) - -! Calculate the total ext. and scat. coefficients - if( present(extcoef) ) then - extcoef(i,j,k,w) = extcoef(i,j,k,w) + & - tau * (grav * rhoa(i,j,k) / delp(i,j,k)) - endif - if( present(scacoef) ) then - scacoef(i,j,k,w) = scacoef(i,j,k,w) + & - ssa * tau * (grav * rhoa(i,j,k) / delp(i,j,k)) - endif - enddo !i - enddo !j - enddo !k + do w = 1, size(wavelengths_profile) + call mie%Query(wavelengths_profile(w),n, & + aerosol(:,:,:,n)*delp/grav, & + rh, tau=tau, ssa=ssa, __RC__) +! Calculate the total ext. and scat. coefficients + if ( present(extcoef) ) then + extcoef(:,:,:,w) = extcoef(:,:,:,w) + & + tau * (grav * rhoa / delp) + endif + if ( present(scacoef) ) then + scacoef(:,:,:,w) = scacoef(:,:,:,w) + & + ssa * tau * (grav * rhoa / delp) + endif enddo !wavelengths_profile enddo !nbins end if !present(extcoef)... @@ -3573,72 +3520,63 @@ subroutine Aero_Compute_Diags (mie_table, km, klid, nbegin, nbins, rlow, rup, ch if( present(exttaufm) ) exttaufm = 0. if( present(scataufm) ) scataufm = 0. - do w = 1, size(wavelengths_vertint) - do n = nbegin, nbins -! do w = 1, size(wavelengths_vertint) - do k = klid, km - do j = j1, j2 - do i = i1, i2 - -! call Chem_MieQuery(mie_table, n, ilam550, & - call Chem_MieQuery(mie_table, n, wavelengths_index_vertint(w), & - aerosol(i,j,k,n)*delp(i,j,k)/grav, & - rh(i,j,k), tau=tau, ssa=ssa) - -! Integrate in the vertical - if( present(exttau) ) exttau(i,j,w) = exttau(i,j,w) + tau - if( present(stexttau) ) then - if (ple(i,j,k) .le. tropp(i,j)) then - stexttau(i,j,w) = stexttau(i,j,w) + tau - elseif(ple(i,j,k) .gt. tropp(i,j) .and. ple(i,j,k-1) .lt. tropp(i,j)) then - stexttau(i,j,w) = stexttau(i,j,w) + log(tropp(i,j)/ple(i,j,k-1))/log(ple(i,j,k)/ple(i,j,k-1))*tau - endif - endif - - if( present(exttaufm) ) then - if( NO3nFlag_ ) then - exttaufm(i,j,w) = exttaufm(i,j,w) + tau - else - exttaufm(i,j,w) = exttaufm(i,j,w) + tau * fPMfm(n) - end if - end if + do w = 1, size(wavelengths_vertint) + do n = nbegin, nbins + call mie%Query(wavelengths_vertint(w), n, & + aerosol(:,:,:,n)*delp/grav, & + rh, tau=tau, ssa=ssa, __RC__) + do k = klid, km +! Integrate in the vertical + if( present(exttau) ) exttau(:,:,w) = exttau(:,:,w) + tau(:,:,k) + if( present(stexttau) ) then + where (ple(:,:,k) .le. tropp) + stexttau(:,:,w) = stexttau(:,:,w) + tau(:,:,k) + elsewhere(ple(:,:,k) .gt. tropp .and. ple(:,:,k-1) .lt. tropp) + stexttau(:,:,w) = stexttau(:,:,w) + log(tropp/ple(:,:,k-1))/log(ple(:,:,k)/ple(:,:,k-1))*tau(:,:,k) + endwhere + endif - if( present(exttau25) ) then - if( NO3nFlag_ ) then - exttau25(i,j,w) = exttau25(i,j,w) + tau - else - exttau25(i,j,w) = exttau25(i,j,w) + tau * fPM25(n) - end if + if( present(exttaufm) ) then + if( NO3nFlag_ ) then + exttaufm(:,:,w) = exttaufm(:,:,w) + tau(:,:,k) + else + exttaufm(:,:,w) = exttaufm(:,:,w) + tau(:,:,k) * fPMfm(n) end if + end if - if( present(scatau) ) scatau(i,j,w) = scatau(i,j,w) + tau*ssa - if( present(stscatau) ) then - if (ple(i,j,k) .le. tropp(i,j)) then - stscatau(i,j,w) = stscatau(i,j,w) + tau*ssa - elseif(ple(i,j,k) .gt. tropp(i,j) .and. ple(i,j,k-1) .lt. tropp(i,j)) then - stscatau(i,j,w) = stscatau(i,j,w) + log(tropp(i,j)/ple(i,j,k-1))/log(ple(i,j,k)/ple(i,j,k-1))*tau*ssa - endif - endif - if( present(scataufm) ) then - if( NO3nFlag_ ) then - scataufm(i,j,w) = scataufm(i,j,w) + tau * ssa - else - scataufm(i,j,w) = scataufm(i,j,w) + tau * ssa * fPMfm(n) - end if + if( present(exttau25) ) then + if( NO3nFlag_ ) then + exttau25(:,:,w) = exttau25(:,:,w) + tau(:,:,k) + else + exttau25(:,:,w) = exttau25(:,:,w) + tau(:,:,k) * fPM25(n) end if - - if( present(scatau25) ) then - if( NO3nFlag_ ) then - scatau25(i,j,w) = scatau25(i,j,w) + tau * ssa - else - scatau25(i,j,w) = scatau25(i,j,w) + tau * ssa * fPM25(n) - end if + end if + + if( present(scatau) ) scatau(:,:,w) = scatau(:,:,w) + tau(:,:,k)*ssa(:,:,k) + if( present(stscatau) ) then + where (ple(:,:,k) .le. tropp) + stscatau(:,:,w) = stscatau(:,:,w) + tau(:,:,k)*ssa(:,:,k) + elsewhere(ple(:,:,k) .gt. tropp .and. ple(:,:,k-1) .lt. tropp) + stscatau(:,:,w) = stscatau(:,:,w) + log(tropp/ple(:,:,k-1))/log(ple(:,:,k)/ple(:,:,k-1))*tau(:,:,k)*ssa(:,:,k) + endwhere + endif + if( present(scataufm) ) then + if( NO3nFlag_ ) then + scataufm(:,:,w) = scataufm(:,:,w) + tau(:,:,k) * ssa(:,:,k) + else + scataufm(:,:,w) = scataufm(:,:,w) + tau(:,:,k) * ssa(:,:,k) * fPMfm(n) end if + end if - enddo !i - enddo !j - enddo !k - enddo !wavelengths_vertint + if( present(scatau25) ) then + if( NO3nFlag_ ) then + scatau25(:,:,w) = scatau25(:,:,w) + tau(:,:,k) * ssa(:,:,k) + else + scatau25(:,:,w) = scatau25(:,:,w) + tau(:,:,k) * ssa(:,:,k) * fPM25(n) + end if + end if + enddo !k + enddo !wavelengths_vertint enddo !nbins endif !present(exttau)... @@ -3652,30 +3590,27 @@ subroutine Aero_Compute_Diags (mie_table, km, klid, nbegin, nbins, rlow, rup, ch do n = nbegin, nbins -! Select the name for species - do k = klid, km - do j = j1, j2 - do i = i1, i2 - call Chem_MieQuery(mie_table, n, ilam470, & - aerosol(i,j,k,n)*delp(i,j,k)/grav, & - rh(i,j,k), tau=tau) - tau470(i,j) = tau470(i,j) + tau - - call Chem_MieQuery(mie_table, n, ilam870, & - aerosol(i,j,k,n)*delp(i,j,k)/grav, & - rh(i,j,k), tau=tau) - tau870(i,j) = tau870(i,j) + tau - enddo +! Select the name for species + call mie%Query(4.70E-7, n, & + aerosol(:,:,:,n)*delp/grav, & + rh, tau=tau, __RC__) + do k = klid, km + tau470 = tau470 + tau(:,:,k) enddo - enddo + call mie%Query(8.70E-7, n, & + aerosol(:,:,:,n)*delp/grav, & + rh, tau=tau, __RC__) + do k = klid, km + tau870 = tau870 + tau(:,:,k) + enddo enddo ! nbins angstrom(i1:i2,j1:j2) = & -log(tau470(i1:i2,j1:j2)/tau870(i1:i2,j1:j2)) / & log(470./870.) endif - + deallocate(tau,ssa) __RETURN__(__SUCCESS__) end subroutine Aero_Compute_Diags !==================================================================== @@ -4271,7 +4206,7 @@ end subroutine hoppelCorrection ! !INTERFACE: ! - subroutine CAEmission (mie_table, km, nbins, cdt, grav, prefix, ratPOM, eAircraftfuel, aircraft_fuel_src, & + subroutine CAEmission (mie, km, nbins, cdt, grav, prefix, ratPOM, eAircraftfuel, aircraft_fuel_src, & aviation_lto_src, aviation_cds_src, aviation_crs_src, & fHydrophobic, pblh, tmpu, rhoa, rh, aerosolPhilic, aerosolPhobic, & delp, aviation_layers, & @@ -4283,7 +4218,7 @@ subroutine CAEmission (mie_table, km, nbins, cdt, grav, prefix, ratPOM, eAircraf implicit NONE ! !INPUT PARAMETERS: - type(Chem_Mie), intent(in) :: mie_table ! mie table + type(GOCART2G_Mie), intent(in) :: mie ! mie table integer, intent(in) :: km ! total model levels integer, intent(in) :: nbins ! number of aerosol size bins real, intent(in) :: cdt ! chemistry model time-step [sec] @@ -4354,11 +4289,13 @@ subroutine CAEmission (mie_table, km, nbins, cdt, grav, prefix, ratPOM, eAircraf real, dimension(:,:), allocatable :: f_bb_ ! scaling factor for BB emissions based on maximum allowed exttau real, dimension(:,:), allocatable :: exttau_bb_ ! increment of exttau due to BB during the current time step - real, allocatable, dimension(:,:,:,:) :: qa_bb_ ! increment of qa due to BB during the current time step (nbins,i1:i2,j1:j2:km) + real, allocatable, dimension(:,:,:,:) :: qa_bb_ ! increment of qa due to BB during the current time step (nbins,i1:i2,j1:j2:km) + ! W.Jiang note, changed to (i1:i2,j1:j2,km,nbins) for efficiency real :: cutoff_bb_exttau - integer :: nch, idx - real :: ilam550 - real :: tau, ssa + integer :: idx + integer :: ilam550, status + real :: wavelength550 + real, dimension(:,:,:), allocatable :: tau character(len=255) :: qname real, parameter :: max_bb_exttau = 30.0 @@ -4367,7 +4304,6 @@ subroutine CAEmission (mie_table, km, nbins, cdt, grav, prefix, ratPOM, eAircraf ! Source function terms for SOA from Anthropogenic VOCs real :: srcSOAanthro = 0.0 - ! Initialize local variables ! -------------------------- i2 = size(rhoa,1) @@ -4471,7 +4407,7 @@ subroutine CAEmission (mie_table, km, nbins, cdt, grav, prefix, ratPOM, eAircraf ! Limit biomass burning emissions ! ------------------------------- - allocate(qa_bb_(nbins,i1:i2,j1:j2,km)) + allocate(qa_bb_(i1:i2,j1:j2,km,nbins)) qa_bb_ = 0.0 p0 = ps @@ -4512,44 +4448,31 @@ subroutine CAEmission (mie_table, km, nbins, cdt, grav, prefix, ratPOM, eAircraf ! ----------------------------------- factor = cdt * grav / delp(:,:,k) - qa_bb_(1,:,:,k) = factor * srcHydrophobic - qa_bb_(2,:,:,k) = factor * srcHydrophilic + qa_bb_(:,:,k,1) = factor * srcHydrophobic + qa_bb_(:,:,k,2) = factor * srcHydrophilic end do K_LOOP_BB - nch = mie_table%nch - ! Get the wavelength indices ! -------------------------- ! Must provide ilam550 for AOT calculation - ilam550 = 1. - if(nch .gt. 1) then - do i = 1, nch - if ( mie_table%channels(i) .ge. 5.49e-7 .and. & - mie_table%channels(i) .le. 5.51e-7) ilam550 = i - enddo - endif + ilam550 = mie%getChannel(5.50e-7) + if (ilam550 <=0) ilam550 = 1 + wavelength550 = mie%getWavelength(ilam550, __RC__) ! Calculate the extinction and/or scattering AOD exttau_bb_(i1:i2,j1:j2) = 0.0 - + allocate(tau(i1:i2,j1:j2,km), source = 0.) do n = 1, nbins ! Select the name for species and the index - do k = 1, km - do j = j1, j2 - do i = i1, i2 - call Chem_MieQuery(mie_table, n, ilam550, & - qa_bb_(n,i,j,k)*delp(i,j,k)/grav, & - rh(i,j,k), tau=tau, ssa=ssa) - + call mie%Query(wavelength550, n, & + qa_bb_(:,:,:,n)*delp(:,:,:)/grav, & + rh, tau=tau, __RC__) + do k = 1, km ! Integrate in the vertical - exttau_bb_(i,j) = exttau_bb_(i,j) + tau - - enddo - enddo - enddo - + exttau_bb_(:,:) = exttau_bb_(:,:) + tau(:,:,k) + enddo enddo ! nbins f_bb_ = 1.0 @@ -4563,7 +4486,7 @@ subroutine CAEmission (mie_table, km, nbins, cdt, grav, prefix, ratPOM, eAircraf enddo enddo - deallocate(qa_bb_) + deallocate(qa_bb_, tau) ! Now update the tracer mixing ratios with the aerosol sources ! ------------------------------------------------------------ @@ -6762,7 +6685,7 @@ end subroutine SU_Wet_Removal !BOP ! !IROUTINE: SU_Compute_Diags - subroutine SU_Compute_Diags ( km, klid, rmed, sigma, rhop, grav, pi, nSO4, mie_table, channels, & + subroutine SU_Compute_Diags ( km, klid, rmed, sigma, rhop, grav, pi, nSO4, mie, & wavelengths_profile, wavelengths_vertint, & tmpu, rhoa, delp, ple, tropp,rh, u, v, & DMS, SO2, SO4, MSA, & @@ -6785,8 +6708,7 @@ subroutine SU_Compute_Diags ( km, klid, rmed, sigma, rhop, grav, pi, nSO4, mie_t real, intent(in) :: grav ! gravity [m/sec] real, intent(in) :: pi ! pi constant integer, intent(in) :: nSO4 ! index of SO4 relative to other internal variables - type(Chem_Mie), intent(in) :: mie_table ! mie table - real, dimension(:), intent(in) :: channels + type(GOCART2G_Mie), intent(in) :: mie ! mie table real, dimension(:), intent(in) :: wavelengths_profile real, dimension(:), intent(in) :: wavelengths_vertint real, pointer, dimension(:,:,:), intent(in) :: tmpu ! temperature [K] @@ -6837,11 +6759,10 @@ subroutine SU_Compute_Diags ( km, klid, rmed, sigma, rhop, grav, pi, nSO4, mie_t ! 29july2020, E.Sherman - refactored for process library ! !Local Variables - integer :: i, j, k, w, i1=1, j1=1, i2, j2, nch - real, allocatable, dimension(:) :: wavelengths_index_profile, wavelengths_index_vertint - real :: tau, ssa + integer :: i, j, k, w, i1=1, j1=1, i2, j2, status + real, dimension(:,:,:), allocatable :: tau, ssa real, dimension(:,:), allocatable :: tau470, tau870 - real :: ilam550, ilam470, ilam870 + integer :: ilam470, ilam870 logical :: do_angstrom real :: rh_, gf, rwet, svol @@ -6849,8 +6770,6 @@ subroutine SU_Compute_Diags ( km, klid, rmed, sigma, rhop, grav, pi, nSO4, mie_t !EOP !------------------------------------------------------------------------- ! Begin - - nch = size(channels) j2 = ubound(tmpu, 2) i2 = ubound(tmpu, 1) @@ -6858,65 +6777,20 @@ subroutine SU_Compute_Diags ( km, klid, rmed, sigma, rhop, grav, pi, nSO4, mie_t ! Get the wavelength indices ! -------------------------- -! Must provide ilam550 for AOT calculation - ilam550 = 1. - ilam470 = 0. - ilam870 = 0. - if(nch .gt. 1) then - do i = 1, nch - if ( channels(i) .ge. 5.49e-7 .and. & - channels(i) .le. 5.51e-7) ilam550 = i - if ( channels(i) .ge. 4.69e-7 .and. & - channels(i) .le. 4.71e-7) ilam470 = i - if ( channels(i) .ge. 8.69e-7 .and. & - channels(i) .le. 8.71e-7) ilam870 = i - enddo - endif - allocate(wavelengths_index_profile(size(wavelengths_profile))) - allocate(wavelengths_index_vertint(size(wavelengths_vertint))) - wavelengths_index_profile = 0. - wavelengths_index_vertint = 0. - - ! Channel values are 4.7e-7 5.5e-7 6.7e-7 8.7e-7 [meter]. Their indices are 1,2,3,4 respectively. - do i = 1, size(wavelengths_profile) - if ((wavelengths_profile(i) .ge. 5.49e-7) .and. (wavelengths_profile(i) .le. 5.51e-7)) then - wavelengths_index_profile(i) = 2. - else if ((wavelengths_profile(i) .ge. 4.69e-7) .and. (wavelengths_profile(i) .le. 4.71e-7)) then - wavelengths_index_profile(i) = 1. - else if ((wavelengths_profile(i) .ge. 6.69e-7) .and. (wavelengths_profile(i) .le. 6.71e-7)) then - wavelengths_index_profile(i) = 3. - else if ((wavelengths_profile(i) .ge. 8.69e-7) .and. (wavelengths_profile(i) .le. 8.71e-7)) then - wavelengths_index_profile(i) = 4. - else - print*,'wavelengths_profile of ',wavelengths_profile(i),' is an invalid value.' - return - end if - end do + ilam470 = mie%getChannel(4.70e-7) + if(ilam470 <= 0) ilam470 = 0 - ! Channel values are 4.7e-7 5.5e-7 6.7e-7 8.7e-7 [meter]. Their indices are 1,2,3,4 respectively. - do i = 1, size(wavelengths_vertint) - if ((wavelengths_vertint(i) .ge. 5.49e-7) .and. (wavelengths_vertint(i) .le. 5.51e-7)) then - wavelengths_index_vertint(i) = 2. - else if ((wavelengths_vertint(i) .ge. 4.69e-7) .and. (wavelengths_vertint(i) .le. 4.71e-7)) then - wavelengths_index_vertint(i) = 1. - else if ((wavelengths_vertint(i) .ge. 6.69e-7) .and. (wavelengths_vertint(i) .le. 6.71e-7)) then - wavelengths_index_vertint(i) = 3. - else if ((wavelengths_vertint(i) .ge. 8.69e-7) .and. (wavelengths_vertint(i) .le. 8.71e-7)) then - wavelengths_index_vertint(i) = 4. - else - print*,'wavelengths_profile of ',wavelengths_profile(i),' is an invalid value.' - return - end if - end do + ilam870 = mie%getChannel(8.70e-7) + if(ilam870 <= 0) ilam870 = 0 ! Determine if going to do Angstrom parameter calculation ! ------------------------------------------------------- do_angstrom = .false. ! If both 470 and 870 channels provided (and not the same) then ! possibly will do Angstrom parameter calculation - if(ilam470 .ne. 0. .and. & - ilam870 .ne. 0. .and. & + if(ilam470 .ne. 0 .and. & + ilam870 .ne. 0 .and. & ilam470 .ne. ilam870) do_angstrom = .true. @@ -7017,32 +6891,27 @@ subroutine SU_Compute_Diags ( km, klid, rmed, sigma, rhop, grav, pi, nSO4, mie_t endif ! Calculate the extinction and/or scattering AOD + allocate(tau(i1:i2,j1:j2,km), source = 0.) + allocate(ssa(i1:i2,j1:j2,km), source = 0.) if( associated(extcoef) .or. associated(scacoef) ) then if (associated(extcoef)) extcoef = 0. if (associated(scacoef)) scacoef = 0. do w = 1, size(wavelengths_profile) - do k = klid, km - do j = j1, j2 - do i = i1, i2 - call Chem_MieQuery(mie_table, 1, wavelengths_index_profile(w), & ! Only SO4 exists in the MieTable, so its index is 1 - SO4(i,j,k)*delp(i,j,k)/grav, & - rh(i,j,k), tau=tau, ssa=ssa) + call mie%Query(wavelengths_profile(w), 1, & ! Only SO4 exists in the MieTable, so its index is 1 + SO4*delp/grav, rh, & + tau=tau, ssa=ssa, __RC__) ! Calculate the total ext. and scat. coefficients - if( associated(extcoef) ) then - extcoef(i,j,k,w) = extcoef(i,j,k,w) + & - tau * (grav * rhoa(i,j,k) / delp(i,j,k)) - endif - if( associated(scacoef) ) then - scacoef(i,j,k,w) = scacoef(i,j,k,w) + & - ssa * tau * (grav * rhoa(i,j,k) / delp(i,j,k)) - endif - - enddo - enddo - enddo + if( associated(extcoef) ) then + extcoef(:,:,:,w) = extcoef(:,:,:,w) + & + tau * (grav * rhoa / delp) + endif + if( associated(scacoef) ) then + scacoef(:,:,:,w) = scacoef(:,:,:,w) + & + ssa * tau * (grav * rhoa / delp) + endif enddo endif @@ -7055,42 +6924,36 @@ subroutine SU_Compute_Diags ( km, klid, rmed, sigma, rhop, grav, pi, nSO4, mie_t if (associated(stscatau)) stscatau = 0. do w = 1, size(wavelengths_vertint) - do k = klid, km - do j = j1, j2 - do i = i1, i2 - call Chem_MieQuery(mie_table, 1, wavelengths_index_vertint(w), & ! Only SO4 exists in the MieTable, so its index is 1 - SO4(i,j,k)*delp(i,j,k)/grav, & - rh(i,j,k), tau=tau, ssa=ssa) - -! Integrate in the vertical - if( associated(exttau) ) then - exttau(i,j,w) = exttau(i,j,w) + tau - endif - - if(associated(stexttau) ) then - if (ple(i,j,k) .le. tropp(i,j)) then - stexttau(i,j,w) = stexttau(i,j,w) + tau - elseif(ple(i,j,k) .gt. tropp(i,j) .and. ple(i,j,k-1) .lt. tropp(i,j)) then - stexttau(i,j,w) = stexttau(i,j,w) + log(tropp(i,j)/ple(i,j,k-1))/log(ple(i,j,k)/ple(i,j,k-1))*tau - endif - endif + call mie%Query(wavelengths_vertint(w), 1, & ! Only SO4 exists in the MieTable, so its index is 1 + SO4*delp/grav, rh, & + tau=tau, ssa=ssa, __RC__) - if( associated(scatau) ) then - scatau(i,j,w) = scatau(i,j,w) + tau*ssa - endif - - if( associated(stscatau) ) then - if (ple(i,j,k) .le. tropp(i,j)) then - stscatau(i,j,w) = stscatau(i,j,w) + tau*ssa - elseif(ple(i,j,k) .gt. tropp(i,j) .and. ple(i,j,k-1) .lt. tropp(i,j)) then - stscatau(i,j,w) = stscatau(i,j,w) + log(tropp(i,j)/ple(i,j,k-1))/log(ple(i,j,k)/ple(i,j,k-1))*tau*ssa - endif + do k = klid, km +! Integrate in the vertical + if ( associated(exttau) ) then + exttau(:,:,w) = exttau(:,:,w) + tau(:,:,k) + endif + + if (associated(stexttau) ) then + where (ple(:,:,k) .le. tropp) + stexttau(:,:,w) = stexttau(:,:,w) + tau(:,:,k) + elsewhere(ple(:,:,k-1) .lt. tropp) + stexttau(:,:,w) = stexttau(:,:,w) + log(tropp/ple(:,:,k-1))/log(ple(:,:,k)/ple(:,:,k-1))*tau(:,:,k) + endwhere + endif + + if ( associated(scatau) ) then + scatau(:,:,w) = scatau(:,:,w) + tau(:,:,k)*ssa(:,:,k) + endif + + if ( associated(stscatau) ) then + where (ple(:,:,k) .le. tropp) + stscatau(:,:,w) = stscatau(:,:,w) + tau(:,:,k)*ssa(:,:,k) + elsewhere(ple(:,:,k-1) .lt. tropp) + stscatau(:,:,w) = stscatau(:,:,w) + log(tropp/ple(:,:,k-1))/log(ple(:,:,k)/ple(:,:,k-1))*tau(:,:,k)*ssa(:,:,k) + endwhere endif - - enddo - enddo - enddo enddo endif @@ -7102,21 +6965,19 @@ subroutine SU_Compute_Diags ( km, klid, rmed, sigma, rhop, grav, pi, nSO4, mie_t tau470(i1:i2,j1:j2) = tiny(1.0) tau870(i1:i2,j1:j2) = tiny(1.0) - do k = klid, km - do j = j1, j2 - do i = i1, i2 - - call Chem_MieQuery(mie_table, 1, ilam470, & ! Only SO4 exists in the MieTable, so its index is 1 - SO4(i,j,k)*delp(i,j,k)/grav, rh(i,j,k), tau=tau) - tau470(i,j) = tau470(i,j) + tau - - call Chem_MieQuery(mie_table, 1, ilam870, & - SO4(i,j,k)*delp(i,j,k)/grav,rh(i,j,k), tau=tau) - tau870(i,j) = tau870(i,j) + tau + call mie%Query(4.70E-7, 1, & ! Only SO4 exists in the MieTable, so its index is 1 + SO4*delp/grav, rh, & + tau=tau, __RC__) + do k = klid, km + tau470 = tau470 + tau(:,:,k) + enddo - enddo - enddo - enddo + call mie%Query(8.70E-7, 1, & + SO4*delp/grav,rh, & + tau=tau, __RC__) + do k = klid, km + tau870 = tau870 + tau(:,:,k) + enddo ! enddo ! nbins angstrom(i1:i2,j1:j2) = & @@ -7150,7 +7011,7 @@ subroutine SU_Compute_Diags ( km, klid, rmed, sigma, rhop, grav, pi, nSO4, mie_t enddo endif endif - + deallocate(tau,ssa) __RETURN__(__SUCCESS__) end subroutine SU_Compute_Diags diff --git a/Process_Library/MieQuery.H b/Process_Library/MieQuery.H new file mode 100644 index 00000000..033687dd --- /dev/null +++ b/Process_Library/MieQuery.H @@ -0,0 +1,256 @@ +#define IDENTITY_(a) a + +#if (RANK_ == 1) +# define SUFFIX_ _1d +# define __DIMS__ : +#elif (RANK_ == 2) +# define SUFFIX_ _2d +# define __DIMS__ :,: +#elif (RANK_ == 3) +# define SUFFIX_ _3d +# define __DIMS__ :,:,: +#endif + +#define __RHINTERP2__(table,rh,ah,bin,sh) reshape((table(rh,bin) * (1-ah) + table(rh+1,bin)*ah), shape(sh)) +#define __RHINTERP3__(table,rh,ah,ch,bin,sh) reshape((table(rh,ch,bin) * (1-ah) + table(rh+1,ch,bin)*ah), shape(sh)) +#define __RHINTERP5__(table,rh,ah,ch,bin,k,l,sh) reshape((table(rh,ch,bin,k,l) * (1-ah) + table(rh+1,ch,bin,k,l)*ah), shape(sh)) + +#define SUBCHANNEL_ IDENTITY_(QueryByChannel)IDENTITY_(SUFFIX_) +#define SUBWAVELENGTH_ IDENTITY_(QueryByWavelength)IDENTITY_(SUFFIX_) + +!------------------------------------------------------------------------- +!BOP +! +! !IROUTINE: GOCART2G_Query --- Return Tau, SSA, etc +! +! +! !INTERFACE: +! + + subroutine SUBCHANNEL_ ( this, channel , bin, q_mass, rh, & + tau, ssa, gasym, bext, bsca, bbck, & + reff, pmom, p11, p22, gf, rhop, rhod, & + vol, area, refr, refi, rc ) + +! !INPUT PARAMETERS: + + class (GOCART2G_Mie), intent(in ) :: this + integer, intent(in ) :: channel ! wave length [m] + integer, intent(in ) :: bin ! bin number + real, intent(in ) :: q_mass(__DIMS__) ! aerosol mass [kg/m2], + real, intent(in ) :: rh (__DIMS__) ! relative himidity + +! !OUTPUT PARAMETERS: + + real, optional, intent(out) :: tau (__DIMS__) ! aerol extinction optical depth + real, optional, intent(out) :: ssa (__DIMS__) ! single scattering albedo + real, optional, intent(out) :: gasym (__DIMS__) ! asymmetry parameter + real, optional, intent(out) :: bext (__DIMS__) ! mass extinction efficiency [m2 (kg dry mass)-1] + real, optional, intent(out) :: bsca (__DIMS__) ! mass scattering efficiency [m2 (kg dry mass)-1] + real, optional, intent(out) :: bbck (__DIMS__) ! mass backscatter efficiency [m2 (kg dry mass)-1] + real, optional, intent(out) :: reff (__DIMS__) ! effective radius (micron) + real, optional, intent(out) :: pmom (__DIMS__,:,:) ! moments of phase function + real, optional, intent(out) :: p11 (__DIMS__) ! P11 phase function at backscatter + real, optional, intent(out) :: p22 (__DIMS__) ! P22 phase function at backscatter + real, optional, intent(out) :: gf (__DIMS__) ! Growth factor (ratio of wet to dry radius) + real, optional, intent(out) :: rhop (__DIMS__) ! Wet particle density [kg m-3] + real, optional, intent(out) :: rhod (__DIMS__) ! Dry particle density [kg m-3] + real, optional, intent(out) :: vol (__DIMS__) ! Wet particle volume [m3 kg-1] + real, optional, intent(out) :: area (__DIMS__) ! Wet particle cross section [m2 kg-1] + real, optional, intent(out) :: refr (__DIMS__) ! Wet particle real part of ref. index + real, optional, intent(out) :: refi (__DIMS__) ! Wet particle imag. part of ref. index + integer, optional, intent(out) :: rc ! error code + +! !DESCRIPTION: +! +! Returns requested parameters from the Mie tables, as a function +! of species, relative humidity, and channel +! +! Notes: Needs some checking, and I still force an interpolation step + +! +! !REVISION HISTORY: +! +! 03Mar2022 Tom Clune refactoring +! 23Mar2005 Colarco +! 11Jul2005 da Silva Standardization. +! +!EOP +!------------------------------------------------------------------------- + + integer :: status + + integer, allocatable, dimension(:) :: irh, qrh + real, allocatable, dimension(:) :: rh_, arh + real, allocatable, dimension(__DIMS__) :: bext_, bsca_ + integer :: i, j + + character(len=*), parameter :: Iam = 'Query' + + if ( present(rc) ) rc = 0 + + rh_ = reshape(min(max(rh,0.),0.99),[size(rh)]) ! no (super) saturation + qrh = max(1,nint((rh_+0.001)*1000.)) ! quantized RH + arh = this%rha(qrh) ! interpolation weight in table + irh = this%rhi(qrh) ! index in table for each RH vcalue + where (irh==this%nrh ) + irh = this%nrh-1 + arh = 1.0 + end where + + if(present(bext) .or. present(tau) .or. present(ssa) ) then + bext_ = __RHINTERP3__(this%bext, irh, arh, channel, bin, bext) + endif + + if(present(bsca) .or. present(ssa) ) then + bsca_ = __RHINTERP3__(this%bsca, irh, arh, channel, bin, bsca) + endif + + if(present(bbck)) then + bbck = __RHINTERP3__(this%bbck, irh, arh, channel, bin, bbck) + endif + + if(present(gasym)) then + gasym = __RHINTERP3__(this%g, irh, arh, channel, bin, gasym) + endif + + if(present(rEff)) then + rEff = 1.E6* __RHINTERP2__(this%rEff, irh, arh, bin, rEff) + endif + + if(present(p11)) then + p11 = __RHINTERP3__(this%p11, irh, arh, channel, bin, p11) + endif + + if(present(p22)) then + p22 = __RHINTERP3__(this%p22, irh, arh, channel, bin, p22) + endif + + if(present(gf)) then + gf = __RHINTERP2__(this%gf, irh, arh, bin, gf) + endif + + if(present(rhod)) then + rhod = __RHINTERP2__(this%rhod, irh, arh, bin, rhod) + endif + + if(present(vol)) then + vol = __RHINTERP2__(this%vol, irh, arh, bin, vol) + endif + + if(present(area)) then + area = __RHINTERP2__(this%area, irh, arh, bin, area) + endif + + if(present(refr)) then + refr = __RHINTERP3__(this%refr, irh, arh, channel, bin, refr) + endif + + if(present(refi)) then + refi = __RHINTERP3__(this%refi, irh, arh, channel, bin, refi) + endif + + if (present(pmom)) then + do j = 1, size(this%pmom,5) + do i = 1, size(this%pmom,4) + pmom(__DIMS__,i,j) = __RHINTERP5__(this%pmom, irh, arh, channel, bin, i, j, rh) + enddo + enddo + endif + + if(present(tau )) tau = bext_ * q_mass + if(present(ssa )) ssa = bsca_/bext_ + if(present(bext )) bext = bext_ + if(present(bsca )) bsca = bsca_ + + end subroutine SUBCHANNEL_ + + +!------------------------------------------------------------------------- +!BOP +! +! !IROUTINE: GOCART2G_Query --- Return Tau, SSA, etc +! +! +! !INTERFACE: +! + + subroutine SUBWAVELENGTH_ ( this, wavelength, bin, q_mass, rh, & + tau, ssa, gasym, bext, bsca, bbck, & + reff, pmom, p11, p22, gf, rhop, rhod, & + vol, area, refr, refi, rc ) + +! !INPUT PARAMETERS: + + class (GOCART2G_Mie), intent(in ) :: this + real, intent(in ) :: wavelength ! wave length [m] + integer, intent(in ) :: bin ! bin number + real, intent(in ) :: q_mass(__DIMS__) ! aerosol mass [kg/m2], + real, intent(in ) :: rh (__DIMS__) ! relative himidity + +! !OUTPUT PARAMETERS: + + real, optional, intent(out) :: tau (__DIMS__) ! aerol extinction optical depth + real, optional, intent(out) :: ssa (__DIMS__) ! single scattering albedo + real, optional, intent(out) :: gasym (__DIMS__) ! asymmetry parameter + real, optional, intent(out) :: bext (__DIMS__) ! mass extinction efficiency [m2 (kg dry mass)-1] + real, optional, intent(out) :: bsca (__DIMS__) ! mass scattering efficiency [m2 (kg dry mass)-1] + real, optional, intent(out) :: bbck (__DIMS__) ! mass backscatter efficiency [m2 (kg dry mass)-1] + real, optional, intent(out) :: reff (__DIMS__) ! effective radius (micron) + real, optional, intent(out) :: pmom (__DIMS__,:,:) ! moments of phase function + real, optional, intent(out) :: p11 (__DIMS__) ! P11 phase function at backscatter + real, optional, intent(out) :: p22 (__DIMS__) ! P22 phase function at backscatter + real, optional, intent(out) :: gf (__DIMS__) ! Growth factor (ratio of wet to dry radius) + real, optional, intent(out) :: rhop (__DIMS__) ! Wet particle density [kg m-3] + real, optional, intent(out) :: rhod (__DIMS__) ! Dry particle density [kg m-3] + real, optional, intent(out) :: vol (__DIMS__) ! Wet particle volume [m3 kg-1] + real, optional, intent(out) :: area (__DIMS__) ! Wet particle cross section [m2 kg-1] + real, optional, intent(out) :: refr (__DIMS__) ! Wet particle real part of ref. index + real, optional, intent(out) :: refi (__DIMS__) ! Wet particle imag. part of ref. index + integer, optional, intent(out) :: rc ! error code + +! !DESCRIPTION: +! +! Returns requested parameters from the Mie tables, as a function +! of species, relative humidity, and channel +! +! Notes: Needs some checking, and I still force an interpolation step + +! +! !REVISION HISTORY: +! +! 03Mar2022 Tom Clune refactoring +! 23Mar2005 Colarco +! 11Jul2005 da Silva Standardization. +! +!EOP +!------------------------------------------------------------------------- + + integer :: status + integer :: channel + + character(len=*), parameter :: Iam = 'QueryByWavelength' + + if ( present(rc) ) rc = 0 + + channel = this%getChannel(wavelength, rc=rc) + if (present(rc)) then + if (rc /=0) return + endif + + call this%Query(channel, bin, q_mass, rh, & + tau, ssa, gasym, bext, bsca, bbck, & + reff, pmom, p11, p22, gf, rhop, rhod, & + vol, area, refr, refi, rc ) + + end subroutine SUBWAVELENGTH_ + +#undef SUBWAVELENGTH_ +#undef SUBCHANNEL_ +#undef __RHINTERP2__ +#undef __RHINTERP3__ +#undef __RHINTERP5__ +#undef __DIMS__ +#undef SUFFIX_ +#undef IDENTITY_ + diff --git a/components.yaml b/components.yaml index 106747ab..263629ed 100644 --- a/components.yaml +++ b/components.yaml @@ -5,34 +5,34 @@ GOCART: env: local: ./env@ remote: ../ESMA_env.git - tag: v3.4.0 + tag: v3.13.0 develop: main cmake: local: ./cmake@ remote: ../ESMA_cmake.git - tag: v3.6.2 + tag: v3.12.0 develop: develop ecbuild: local: ./cmake@/ecbuild@ remote: ../ecbuild.git - tag: geos/v1.0.6 + tag: geos/v1.2.0 HEMCO: local: ./ESMF/HEMCO_GridComp@ remote: ../HEMCO.git - branch: geos/v2.2.1 + branch: geos/v2.2.3 GMAO_Shared: local: ./ESMF/Shared/GMAO_Shared@ remote: ../GMAO_Shared.git - tag: v1.4.10 + tag: v1.5.3 sparse: ./config/GMAO_Shared.sparse develop: main MAPL: local: ./ESMF/Shared/MAPL@ remote: ../MAPL.git - tag: v2.8.6 + tag: v2.19.0 develop: develop