diff --git a/aqm_files.cmake b/aqm_files.cmake index 168b6ba4..67709fdc 100644 --- a/aqm_files.cmake +++ b/aqm_files.cmake @@ -87,7 +87,6 @@ set(VDIFF "${CCTM_ROOT}/vdiff/acm2") set(localCCTM "src/model/src") list(APPEND aqm_CCTM_files ${AERO}/AERO_DATA.F - ${AERO}/aero_depv.F ${AERO}/aero_driver.F ${AERO}/AERO_EMIS.F ${AERO}/AEROMET_DATA.F @@ -236,8 +235,12 @@ list(APPEND aqm_CCTM_files ${localCCTM}/o3totcol.f ${localCCTM}/vdiffacmx.F ${localCCTM}/PTMAP.F + ${localCCTM}/PT3D_DATA_MOD.F ${localCCTM}/PT3D_DEFN.F + ${localCCTM}/PT3D_FIRE_DEFN.F + ${localCCTM}/PT3D_STKS_DEFN.F ${localCCTM}/ASX_DATA_MOD.F ${localCCTM}/DUST_EMIS.F ${localCCTM}/AERO_PHOTDATA.F + ${localCCTM}/aero_depv.F ) diff --git a/examples/aqm.rc b/examples/aqm.rc index 637196e1..13246deb 100644 --- a/examples/aqm.rc +++ b/examples/aqm.rc @@ -128,6 +128,11 @@ myemis_frequency: hourly # supported schemes are: sofiev, none # #myemis_plume_rise: sofiev +# +# Optionally, provide an empirical weight for the top plume layer. +# If this option is absent, or a value < 0.0 is provided, fire +# emissions will be distributed linearly among the plume layers +#myemis_plume_top_fraction: 0.8 # list of emission species # input emissions can be manipulated by adding multiple diff --git a/src/aqm_cap.F90 b/src/aqm_cap.F90 index dfec4a49..15a32313 100644 --- a/src/aqm_cap.F90 +++ b/src/aqm_cap.F90 @@ -8,11 +8,12 @@ module AQM use NUOPC_Model, inheritModel => SetServices use aqm_comp_mod + use aqm_const_mod, only: rad_to_deg implicit none ! -- import fields - integer, parameter :: importFieldCount = 35 + integer, parameter :: importFieldCount = 36 character(len=*), dimension(importFieldCount), parameter :: & importFieldNames = (/ & "canopy_moisture_storage ", & @@ -30,6 +31,7 @@ module AQM "inst_net_sw_flx ", & "inst_pbl_height ", & "inst_pres_height_surface ", & + "inst_pres_interface ", & "inst_pres_levels ", & "inst_rainfall_amount ", & "inst_sensi_heat_flx ", & @@ -61,9 +63,6 @@ module AQM private - real(ESMF_KIND_R8), parameter :: pi = 3.1415926535897931 - real(ESMF_KIND_R8), parameter :: rad2deg = 180./pi - public SetServices !----------------------------------------------------------------------------- @@ -264,7 +263,7 @@ subroutine DataInitialize(model, rc) real(ESMF_KIND_R8) :: dts type(ESMF_Time) :: startTime type(ESMF_TimeInterval) :: TimeStep - type(ESMF_CoordSys_Flag) :: aqmGridCoordSys + type(ESMF_CoordSys_Flag) :: coordSys character(len=ESMF_MAXSTR) :: msgString, name @@ -413,7 +412,7 @@ subroutine DataInitialize(model, rc) end if ! -- get local coordinate arrays - call ESMF_GridGet(grid, coordSys=aqmGridCoordSys, rc=rc) + call ESMF_GridGet(grid, coordSys=coordSys, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__)) & @@ -427,13 +426,29 @@ subroutine DataInitialize(model, rc) file=__FILE__)) & return ! bail out - if (aqmGridCoordSys == ESMF_COORDSYS_SPH_RAD) coord = coord * rad2deg - - call aqm_model_domain_coord_set(item, coord, de=localDe, rc=rc) - - if (aqm_rc_check(rc)) then + if (coordSys == ESMF_COORDSYS_SPH_RAD) then + call aqm_model_domain_coord_set(item, coord, scale=rad_to_deg, de=localDe, rc=rc) + if (aqm_rc_check(rc)) then + call ESMF_LogSetError(ESMF_RC_INTNRL_BAD, & + msg="Failed to set coordinates for air quality model", & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc) + return ! bail out + end if + else if (coordSys == ESMF_COORDSYS_SPH_DEG) then + call aqm_model_domain_coord_set(item, coord, de=localDe, rc=rc) + if (aqm_rc_check(rc)) then + call ESMF_LogSetError(ESMF_RC_INTNRL_BAD, & + msg="Failed to set coordinates for air quality model", & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc) + return ! bail out + end if + else call ESMF_LogSetError(ESMF_RC_INTNRL_BAD, & - msg="Failed to set coordinates for air quality model", & + msg="Unsupported coordinate system - Failed to set coordinates for air quality model", & line=__LINE__, & file=__FILE__, & rcToReturn=rc) @@ -442,6 +457,7 @@ subroutine DataInitialize(model, rc) end do end do + deallocate(minIndexPDe, maxIndexPDe, minIndexPTile, maxIndexPTile, & computationalLBound, computationalUBound, & deToTileMap, localDeToDeMap, stat=localrc) diff --git a/src/aqm_comp_mod.F90 b/src/aqm_comp_mod.F90 index a1c937db..c2553167 100644 --- a/src/aqm_comp_mod.F90 +++ b/src/aqm_comp_mod.F90 @@ -477,6 +477,12 @@ subroutine aqm_comp_import(state, fieldNames, rc) line=__LINE__, & file=__FILE__)) & return ! bail + case ("inst_pres_interface") + call ESMF_FieldGet(field, localDe=localDe, farrayPtr=stateIn % pri, rc=rc) + if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__)) & + return ! bail case ("inst_pres_levels") call ESMF_FieldGet(field, localDe=localDe, farrayPtr=stateIn % prl, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, & diff --git a/src/drv/cmaq_mod.F90 b/src/drv/cmaq_mod.F90 index 193bb285..ac476109 100644 --- a/src/drv/cmaq_mod.F90 +++ b/src/drv/cmaq_mod.F90 @@ -443,11 +443,12 @@ subroutine cmaq_emis_init(rc) type(aqm_internal_emis_type), pointer :: em ! -- local parameters - character(len=*), parameter :: etype(4) = (/ & + character(len=*), parameter :: etype(5) = (/ & "anthropogenic", & "biogenic ", & "fengsha ", & - "gbbepx " & + "gbbepx ", & + "point-source " & /) ! -- begin @@ -465,7 +466,7 @@ subroutine cmaq_emis_init(rc) if (associated(em)) then select case (trim(etype(item))) - case ("anthropogenic", "gbbepx") + case ("anthropogenic", "gbbepx", "point-source") ! -- check if internal emissions reference table was allocated if (aqm_rc_test(.not.associated(em % table), & diff --git a/src/io/aqmio/aqmio.F90 b/src/io/aqmio/aqmio.F90 index fc937e00..7abcf587 100644 --- a/src/io/aqmio/aqmio.F90 +++ b/src/io/aqmio/aqmio.F90 @@ -39,6 +39,7 @@ module AQMIO public :: AQMIO_Close public :: AQMIO_Read public :: AQMIO_ReadTimes + public :: AQMIO_DataRead public :: AQMIO_Sync public :: AQMIO_Write @@ -46,16 +47,17 @@ module AQMIO !------------------------------------------------------------------------------ - function AQMIO_Create(grid, vm, rc) + function AQMIO_Create(grid, vm, allpes, rc) type(ESMF_Grid), intent(in) :: grid type(ESMF_VM), intent(in), optional :: vm + logical, intent(in), optional :: allpes integer, intent(out), optional :: rc type(ESMF_GridComp) :: AQMIO_Create ! -- local variables integer :: localrc - integer :: i, localDe, localDeCount, localpe, peCount, npe + integer :: i, iope, localDe, localDeCount, localpe, peCount, npe integer :: tile, tileCount integer, dimension(:), allocatable :: localTile, tileToPet, pes, recvpes type(ESMF_GridComp) :: IOComp, taskComp @@ -256,6 +258,11 @@ function AQMIO_Create(grid, vm, rc) file=__FILE__, & rcToReturn=rc)) return ! bail out + iope = 0 + if (present(allpes)) then + if (allpes) iope = localpe + end if + ! -- flag PET if local I/O must be performed do localDe = 0, localDeCount - 1 call ESMF_GridCompGet(is % IO % IOLayout(localDe) % taskComp, vm=tasksVM, rc=localrc) @@ -268,7 +275,7 @@ function AQMIO_Create(grid, vm, rc) line=__LINE__, & file=__FILE__, & rcToReturn=rc)) return ! bail out - is % IO % IOLayout(localDe) % localIOflag = (localpe == 0) + is % IO % IOLayout(localDe) % localIOflag = (localpe == iope) end do AQMIO_Create = IOComp @@ -1177,7 +1184,7 @@ subroutine AQMIO_FieldRead(IO, field, & ! -- local variables integer :: localrc - integer :: ilen, jlen, lbuf, lde, localpe, rank + integer :: ilen, jlen, lbuf, lde, rank integer :: varId, ncStatus, ndims, xtype integer :: kmin, kmax, klen, uid integer, dimension(3) :: elb, eub @@ -1265,13 +1272,7 @@ subroutine AQMIO_FieldRead(IO, field, & file=__FILE__, & rcToReturn=rc)) return ! bail out - call ESMF_VMGet(vm, localPet=localpe, rc=localrc) - if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & - line=__LINE__, & - file=__FILE__, & - rcToReturn=rc)) return ! bail out - - if (localpe == 0) then + if (IO % IOLayout(lde) % localIOflag) then if (IO % IOLayout(lde) % ncid > 0) then #if HAVE_NETCDF dataSetName = "NetCDF data set" @@ -1634,6 +1635,251 @@ subroutine AQMIO_FieldRead(IO, field, & end subroutine AQMIO_FieldRead +!------------------------------------------------------------------------------ + + subroutine AQMIO_DataRead(IOComp, fArray, variableName, timeSlice, localDe, rc) + type(ESMF_GridComp), intent(inout) :: IOComp + real(ESMF_KIND_R4), pointer :: fArray(:) + character(len=*), intent(in) :: variableName + integer, intent(in), optional :: timeSlice + integer, intent(in), optional :: localDe + integer, intent(out), optional :: rc + + ! -- local variables + integer :: localrc + integer :: ilen, lde, rank + integer :: varId, ncStatus, ndims, xtype + integer :: uid + integer :: ibuf(1) + integer, dimension(:), allocatable :: dimids + integer, dimension(:), allocatable :: elemCount + integer, dimension(:), allocatable :: elemStart + real(ESMF_KIND_R4), dimension(:), allocatable :: buf + real(ESMF_KIND_R4), dimension(:), pointer :: fp + character(len=ESMF_MAXSTR) :: dataSetName + type(ESMF_VM) :: vm + type(ioWrapper) :: is + type(ioData), pointer :: IO + + ! -- begin + if (present(rc)) rc = ESMF_SUCCESS + + if (.not.ESMF_GridCompIsPetLocal(IOComp)) return + + call ESMF_GridCompGetInternalState(IOComp, is, localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) return ! bail out + + if (.not.associated(is % IO)) return + if (.not.associated(is % IO % IOLayout)) return + + IO => is % IO + + lde = 0 + if (present(localDe)) lde = localDe + + call ESMF_GridCompGet(IO % IOLayout(lde) % taskComp, vm=vm, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) return ! bail out + + if (IO % IOLayout(lde) % localIOflag) then + if (IO % IOLayout(lde) % ncid > 0) then +#if HAVE_NETCDF + dataSetName = "NetCDF data set" + + ncStatus = nf90_inquire(IO % IOLayout(lde) % ncid, unlimitedDimId=uid) + if (ESMF_LogFoundNetCDFError(ncerrToCheck=ncStatus, & + msg="No unlimited dimension (time) in "//trim(dataSetName), & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) return ! bail out + + ncStatus = nf90_inq_varid(IO % IOLayout(lde) % ncid, trim(variableName), varId) + if (ESMF_LogFoundNetCDFError(ncerrToCheck=ncStatus, & + msg="Variable "//trim(variableName)//" not defined in "//trim(dataSetName), & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) return ! bail out + + ncStatus = nf90_inquire_variable(IO % IOLayout(lde) % ncid, varId, & + xtype=xtype, ndims=ndims) + if (ESMF_LogFoundNetCDFError(ncerrToCheck=ncStatus, & + msg="Error inquiring variable "//trim(variableName)//" in "//trim(dataSetName), & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) return ! bail out + + call AQMIO_VariableCheckType(variableName, xtype, ESMF_TYPEKIND_R4, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) return ! bail out + + allocate(elemStart(ndims), elemCount(ndims), stat=localrc) + if (ESMF_LogFoundAllocError(statusToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) return ! bail out + elemStart = 1 + elemCount = 1 + + allocate(dimids(ndims), stat=localrc) + if (ESMF_LogFoundAllocError(statusToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) return ! bail out + + ncStatus = nf90_inquire_variable(IO % IOLayout(lde) % ncid, varId, dimIds=dimids) + if (ESMF_LogFoundNetCDFError(ncerrToCheck=ncStatus, & + msg="Error inquiring variable "//trim(variableName)//" in "//trim(dataSetName), & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) return ! bail out + + if (uid == -1) then + if (present(timeSlice)) then + if (timeSlice == 1) then + call ESMF_LogWrite("No time record found in "//trim(dataSetName) & + // " - proceed only for first time step", & + ESMF_LOGMSG_WARNING, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) return ! bail out + else + call ESMF_LogSetError(ESMF_RC_NOT_FOUND, & + msg="No time record found in "//dataSetName, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc) + return ! bail out + end if + end if + else + if (dimids(ndims) == uid) then + if (present(timeSlice)) elemStart(ndims) = timeSlice + ndims = ndims - 1 + else + if (present(timeSlice)) then + call ESMF_LogSetError(ESMF_RC_NOT_FOUND, & + msg="No time record found for variable "// variableName, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc) + return ! bail out + end if + end if + end if + + rank = 1 + + if (rank /= ndims) localrc = ESMF_RC_ARG_INCOMP + if (ESMF_LogFoundError(rcToCheck=localrc, & + msg="Variable rank incompatible with netCDF variable "//trim(variableName), & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) return ! bail out + + ! -- allocate array according to dimension on file + ncStatus = nf90_inquire_dimension(IO % IOLayout(lde) % ncid, dimids(ndims), len=ibuf(1)) + if (ESMF_LogFoundNetCDFError(ncerrToCheck=ncStatus, & + msg="Error inquiring dimension for "//trim(variableName)//" in "//trim(dataSetName), & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) return ! bail out + end if + end if + + call ESMF_VMBroadcast(vm, ibuf, 1, 0, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) return ! bail out + + ilen = ibuf(1) + + if (associated(fArray)) then + ilen = min(ilen,size(fArray)) + else + allocate(fp(ilen), stat=localrc) + if (ESMF_LogFoundAllocError(statusToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) return ! bail out + fArray => fp + end if + + fArray = 0._ESMF_KIND_R4 + + if (IO % IOLayout(lde) % localIOflag) then + if (IO % IOLayout(lde) % ncid > 0) then + + deallocate(dimids, stat=localrc) + if (ESMF_LogFoundDeallocError(statusToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) return ! bail out + + elemCount(1) = ilen + + allocate(buf(ilen), stat=localrc) + if (ESMF_LogFoundAllocError(statusToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) return ! bail out + + ncStatus = nf90_get_var(IO % IOLayout(lde) % ncid, varId, fArray, & + start=elemStart, count=elemCount) + if (ESMF_LogFoundNetCDFError(ncerrToCheck=ncStatus, & + msg="Error reading field "//trim(variableName)//" from "//trim(dataSetName), & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) return ! bail out + + deallocate(elemStart, elemCount, stat=localrc) + if (ESMF_LogFoundDeallocError(statusToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) return ! bail out + + end if + end if + + call ESMF_VMBroadcast(vm, fArray, ilen, 0, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) return ! bail out + + if (IO % IOLayout(lde) % localIOflag) then + if (IO % IOLayout(lde) % ncid > 0) then + ! -- nothing else to do +#else + call ESMF_LogSetError(rcToCheck=ESMF_RC_LIB_NOT_PRESENT, & + msg="- netCDF support is unavailable", & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc) + return +#endif + else if (IO % IOLayout(lde) % iounit > 0) then + + call ESMF_LogSetError(rcToCheck=ESMF_RC_LIB_NOT_PRESENT, & + msg="- binary format is not supported", & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc) + return + + end if + end if + + end subroutine AQMIO_DataRead + !------------------------------------------------------------------------------ subroutine AQMIO_TimesRead(IO, variableName, timesList, localDe, rc) @@ -1828,7 +2074,7 @@ subroutine AQMIO_FieldWrite(IO, field, & ! -- local variables integer :: localrc - integer :: ilen, jlen, lbuf, lde, localpe + integer :: ilen, jlen, lbuf, lde integer :: varId, ncStatus integer :: ndims, rank integer :: kmin, kmax, klen @@ -1882,12 +2128,6 @@ subroutine AQMIO_FieldWrite(IO, field, & file=__FILE__, & rcToReturn=rc)) return ! bail out - call ESMF_VMGet(vm, localPet=localpe, rc=localrc) - if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & - line=__LINE__, & - file=__FILE__, & - rcToReturn=rc)) return ! bail out - if (typekind == ESMF_TYPEKIND_I4) then allocate(buf_i4(minIndexPTile(1):maxIndexPTile(1), & @@ -2074,7 +2314,7 @@ subroutine AQMIO_FieldWrite(IO, field, & return ! bail out end if - if (localpe == 0) then + if (IO % IOLayout(lde) % localIOflag) then if (IO % IOLayout(lde) % ncid > 0) then #if HAVE_NETCDF dataSetName = "NetCDF data set" diff --git a/src/model/Makefile.am b/src/model/Makefile.am index c5fed5d3..cff39671 100644 --- a/src/model/Makefile.am +++ b/src/model/Makefile.am @@ -163,7 +163,6 @@ libCCTM_a_SOURCES += \ PHOT = $(CCTM)/phot/inline libPHOT = $(PHOT)/$(libCCTM)- libCCTM_a_SOURCES += \ - $(PHOT)/AERO_PHOTDATA.F \ $(PHOT)/CLOUD_OPTICS.F \ $(PHOT)/complex_number_module.F90 \ $(PHOT)/CSQY_DATA.F \ @@ -244,13 +243,16 @@ libCCTM_a_SOURCES += \ # local replacements for CCTM files liblocalCCTM = $(localCCTM)/$(libCCTM)- libCCTM_a_SOURCES += \ - $(localCCTM)/o3totcol.f \ - $(localCCTM)/vdiffacmx.F \ - $(localCCTM)/PTMAP.F \ - $(localCCTM)/PT3D_DEFN.F \ - $(localCCTM)/ASX_DATA_MOD.F \ - $(localCCTM)/DUST_EMIS.F - + ${localCCTM}/o3totcol.f \ + ${localCCTM}/vdiffacmx.F \ + ${localCCTM}/PTMAP.F \ + ${localCCTM}/PT3D_DATA_MOD.F \ + ${localCCTM}/PT3D_DEFN.F \ + ${localCCTM}/PT3D_FIRE_DEFN.F \ + ${localCCTM}/PT3D_STKS_DEFN.F \ + ${localCCTM}/ASX_DATA_MOD.F \ + ${localCCTM}/DUST_EMIS.F \ + ${localCCTM}/AERO_PHOTDATA.F libCCTM_a_CPPFLAGS = -DSUBST_FILES_ID=\"FILES_CTM.EXT\" @@ -272,7 +274,8 @@ libCCTM_a_FFLAGS = $(CCTM_FFLAGS) libCCTM_a_FFLAGS += \ -I $(ICL)/filenames -I $(ICL)/const -I $(ICL)/emctrl -I$(ICL)/mpi \ -I $(top_builddir)/src/io/ioapi -I $(top_builddir)/src/shr \ - -I $(top_builddir)/src/io/aqmio + -I $(top_builddir)/src/io/aqmio \ + -I $(AERO) libCCTM_a_FCFLAGS = $(CCTM_FCFLAGS) libCCTM_a_FCFLAGS += -I $(top_builddir)/src/io/ioapi @@ -546,12 +549,6 @@ $(libPA)pa_update.$(OBJEXT) : $(ICL)/const/CONST.EXT $(ICL)/emctrl/EMISPRM.EXT $ $(libGRID)PAGRD_DEFN.$(OBJEXT) $(libUTIL)UTILIO_DEFN.$(OBJEXT) # phot -$(libPHOT)AERO_PHOTDATA.$(OBJEXT) : $(ICL)/const/CONST.EXT \ - $(libAERO)AERO_DATA.$(OBJEXT) $(libAERO)AEROMET_DATA.$(OBJEXT) \ - $(libSPCS)CGRID_SPCS.$(OBJEXT) $(libPHOT)twoway_rrtmg_aero_optics.$(OBJEXT) \ - $(libPHOT)CSQY_DATA.$(OBJEXT) $(libGRID)GRID_CONF.$(OBJEXT) \ - $(libMECHS)RXNS_DATA_MODULE.$(OBJEXT) $(libAERO)SOA_DEFN.$(OBJEXT) \ - $(libUTIL)UTILIO_DEFN.$(OBJEXT) $(libPHOT)CLOUD_OPTICS.$(OBJEXT) : \ $(libPHOT)CSQY_DATA.$(OBJEXT) $(libGRID)GRID_CONF.$(OBJEXT) \ $(libUTIL)UTILIO_DEFN.$(OBJEXT) @@ -563,7 +560,7 @@ $(libPHOT)opphot.$(OBJEXT) : $(ICL)/filenames/FILES_CTM.EXT \ $(libPHOT)PHOT_MOD.$(OBJEXT) $(libMECHS)RXNS_DATA_MODULE.$(OBJEXT) \ $(libUTIL)UTILIO_DEFN.$(OBJEXT) $(libPHOT)phot.$(OBJEXT) : $(ICL)/filenames/FILES_CTM.EXT \ - $(libAERO)AERO_DATA.$(OBJEXT) $(libPHOT)AERO_PHOTDATA.$(OBJEXT) \ + $(libAERO)AERO_DATA.$(OBJEXT) $(liblocalCCTM)AERO_PHOTDATA.$(OBJEXT) \ $(libSPCS)CGRID_SPCS.$(OBJEXT) $(libPHOT)CLOUD_OPTICS.$(OBJEXT) \ $(libSTENEX)noop_modules.$(OBJEXT) $(libGRID)PCGRID_DEFN.$(OBJEXT) \ $(libPHOT)PHOT_MET_DATA.$(OBJEXT) $(libPHOT)PHOT_MOD.$(OBJEXT) \ @@ -573,7 +570,7 @@ $(libPHOT)PHOT_MET_DATA.$(OBJEXT) : $(ICL)/const/CONST.EXT $(ICL)/filenames/FILE $(libPHOT)CLOUD_OPTICS.$(OBJEXT) $(libGRID)GRID_CONF.$(OBJEXT) \ $(libUTIL)UTILIO_DEFN.$(OBJEXT) $(libPHOT)PHOT_MOD.$(OBJEXT) : $(ICL)/const/CONST.EXT \ - $(libPHOT)AERO_PHOTDATA.$(OBJEXT) $(libPHOT)CLOUD_OPTICS.$(OBJEXT) \ + $(liblocalCCTM)AERO_PHOTDATA.$(OBJEXT) $(libPHOT)CLOUD_OPTICS.$(OBJEXT) \ $(libPHOT)CSQY_DATA.$(OBJEXT) $(libMECHS)RXNS_DATA_MODULE.$(OBJEXT) \ $(libUTIL)UTILIO_DEFN.$(OBJEXT) $(libPHOT)PHOTOLYSIS_ALBEDO.$(OBJEXT) : $(ICL)/filenames/FILES_CTM.EXT \ @@ -676,10 +673,19 @@ $(liblocalCCTM)vdiffacmx.$(OBJEXT) : $(ICL)/filenames/FILES_CTM.EXT \ $(libVDIFF)VDIFF_MAP.$(OBJEXT) $(liblocalCCTM)PTMAP.$(OBJEXT) : $(libAERO)AERO_DATA.$(OBJEXT) \ $(libSPCS)CGRID_SPCS.$(OBJEXT) $(libUTIL)UTILIO_DEFN.$(OBJEXT) -$(liblocalCCTM)PT3D_DEFN.$(OBJEXT) : $(libAERO)AERO_DATA.$(OBJEXT) \ +$(liblocalCCTM)PT3D_DEFN.$(OBJEXT) : $(libGRID)GRID_CONF.$(OBJEXT) \ + $(liblocalCCTM)PTMAP.$(OBJEXT) $(liblocalCCTM)PT3D_FIRE_DEFN.$(OBJEXT) \ + $(liblocalCCTM)PT3D_STKS_DEFN.$(OBJEXT) +$(liblocalCCTM)PT3D_FIRE_DEFN.$(OBJEXT) : $(libAERO)AERO_DATA.$(OBJEXT) \ $(libGRID)GRID_CONF.$(OBJEXT) $(libSPCS)CGRID_SPCS.$(OBJEXT) \ + $(liblocalCCTM)PT3D_DATA_MOD.$(OBJEXT) \ $(liblocalCCTM)PTMAP.$(OBJEXT) $(libMECHS)RXNS_DATA_MODULE.$(OBJEXT) \ $(libEMIS)STK_EMIS.$(OBJEXT) $(libUTIL)UTILIO_DEFN.$(OBJEXT) +$(liblocalCCTM)PT3D_STKS_DEFN.$(OBJEXT) : $(libAERO)AERO_DATA.$(OBJEXT) \ + $(liblocalCCTM)ASX_DATA_MOD.$(OBJEXT) \ + $(libGRID)GRID_CONF.$(OBJEXT) $(libSPCS)CGRID_SPCS.$(OBJEXT) \ + $(liblocalCCTM)PT3D_DATA_MOD.$(OBJEXT) $(liblocalCCTM)PTMAP.$(OBJEXT) \ + $(libEMIS)STK_EMIS.$(OBJEXT) $(libUTIL)UTILIO_DEFN.$(OBJEXT) $(liblocalCCTM)ASX_DATA_MOD.$(OBJEXT) : $(ICL)/const/CONST.EXT $(ICL)/filenames/FILES_CTM.EXT $(ICL)/mpi/PE_COMM.EXT \ $(libDEPV)DEPVVARS.$(OBJEXT) $(libGRID)GRID_CONF.$(OBJEXT) \ $(libDEPV)LSM_MOD.$(OBJEXT) $(libSTENEX)noop_modules.$(OBJEXT) \ @@ -688,3 +694,9 @@ $(liblocalCCTM)DUST_EMIS.$(OBJEXT) : $(ICL)/const/CONST.EXT $(ICL)/filenames/FIL $(libAERO)AERO_DATA.$(OBJEXT) $(liblocalCCTM)ASX_DATA_MOD.$(OBJEXT) \ $(libGRID)GRID_CONF.$(OBJEXT) $(libGRID)HGRD_DEFN.$(OBJEXT) \ $(libEMIS)LUS_DEFN.$(OBJEXT) $(libUTIL)UTILIO_DEFN.$(OBJEXT) +$(liblocalCCTM)AERO_PHOTDATA.$(OBJEXT) : $(ICL)/const/CONST.EXT \ + $(libAERO)AERO_DATA.$(OBJEXT) $(libAERO)AEROMET_DATA.$(OBJEXT) \ + $(libSPCS)CGRID_SPCS.$(OBJEXT) $(libPHOT)twoway_rrtmg_aero_optics.$(OBJEXT) \ + $(libPHOT)CSQY_DATA.$(OBJEXT) $(libGRID)GRID_CONF.$(OBJEXT) \ + $(libMECHS)RXNS_DATA_MODULE.$(OBJEXT) $(libAERO)SOA_DEFN.$(OBJEXT) \ + $(libUTIL)UTILIO_DEFN.$(OBJEXT) diff --git a/src/model/Makefile.in b/src/model/Makefile.in index 4177970a..078af39b 100644 --- a/src/model/Makefile.in +++ b/src/model/Makefile.in @@ -182,7 +182,6 @@ am_libCCTM_a_OBJECTS = $(AERO)/libCCTM_a-AERO_DATA.$(OBJEXT) \ $(MECHS)/libCCTM_a-RXNS_FUNC_MODULE.$(OBJEXT) \ $(PA)/libCCTM_a-PA_DEFN.$(OBJEXT) \ $(PA)/libCCTM_a-pa_update.$(OBJEXT) \ - $(PHOT)/libCCTM_a-AERO_PHOTDATA.$(OBJEXT) \ $(PHOT)/libCCTM_a-CLOUD_OPTICS.$(OBJEXT) \ $(PHOT)/libCCTM_a-complex_number_module.$(OBJEXT) \ $(PHOT)/libCCTM_a-CSQY_DATA.$(OBJEXT) \ @@ -232,12 +231,16 @@ am_libCCTM_a_OBJECTS = $(AERO)/libCCTM_a-AERO_DATA.$(OBJEXT) \ $(VDIFF)/libCCTM_a-VDIFF_DIAG.$(OBJEXT) \ $(VDIFF)/libCCTM_a-VDIFF_MAP.$(OBJEXT) \ $(VDIFF)/libCCTM_a-vdiffproc.$(OBJEXT) \ - $(localCCTM)/libCCTM_a-o3totcol.$(OBJEXT) \ - $(localCCTM)/libCCTM_a-vdiffacmx.$(OBJEXT) \ - $(localCCTM)/libCCTM_a-PTMAP.$(OBJEXT) \ - $(localCCTM)/libCCTM_a-PT3D_DEFN.$(OBJEXT) \ - $(localCCTM)/libCCTM_a-ASX_DATA_MOD.$(OBJEXT) \ - $(localCCTM)/libCCTM_a-DUST_EMIS.$(OBJEXT) + ${localCCTM}/libCCTM_a-o3totcol.$(OBJEXT) \ + ${localCCTM}/libCCTM_a-vdiffacmx.$(OBJEXT) \ + ${localCCTM}/libCCTM_a-PTMAP.$(OBJEXT) \ + ${localCCTM}/libCCTM_a-PT3D_DATA_MOD.$(OBJEXT) \ + ${localCCTM}/libCCTM_a-PT3D_DEFN.$(OBJEXT) \ + ${localCCTM}/libCCTM_a-PT3D_FIRE_DEFN.$(OBJEXT) \ + ${localCCTM}/libCCTM_a-PT3D_STKS_DEFN.$(OBJEXT) \ + ${localCCTM}/libCCTM_a-ASX_DATA_MOD.$(OBJEXT) \ + ${localCCTM}/libCCTM_a-DUST_EMIS.$(OBJEXT) \ + ${localCCTM}/libCCTM_a-AERO_PHOTDATA.$(OBJEXT) libCCTM_a_OBJECTS = $(am_libCCTM_a_OBJECTS) AM_V_P = $(am__v_P_@AM_V@) am__v_P_ = $(am__v_P_@AM_DEFAULT_V@) @@ -483,12 +486,12 @@ libCCTM_a_SOURCES = $(AERO)/AERO_DATA.F $(AERO)/aero_depv.F \ $(ICL)/emctrl/EMISPRM.EXT $(ICL)/filenames/FILES_CTM.EXT \ $(ICL)/mpi/PE_COMM.EXT $(INIT)/initscen.F $(INIT)/load_cgrid.F \ $(MECHS)/RXNS_DATA_MODULE.F90 $(MECHS)/RXNS_FUNC_MODULE.F90 \ - $(PA)/PA_DEFN.F $(PA)/pa_update.F $(PHOT)/AERO_PHOTDATA.F \ - $(PHOT)/CLOUD_OPTICS.F $(PHOT)/complex_number_module.F90 \ - $(PHOT)/CSQY_DATA.F $(PHOT)/OMI_1979_to_2015.dat \ - $(PHOT)/opphot.F $(PHOT)/phot.F $(PHOT)/PHOT_MET_DATA.F \ - $(PHOT)/PHOT_MOD.F $(PHOT)/PHOTOLYSIS_ALBEDO.F \ - $(PHOT)/PHOT_OPTICS.dat $(PHOT)/SEAS_STRAT_O3_MIN.F \ + $(PA)/PA_DEFN.F $(PA)/pa_update.F $(PHOT)/CLOUD_OPTICS.F \ + $(PHOT)/complex_number_module.F90 $(PHOT)/CSQY_DATA.F \ + $(PHOT)/OMI_1979_to_2015.dat $(PHOT)/opphot.F $(PHOT)/phot.F \ + $(PHOT)/PHOT_MET_DATA.F $(PHOT)/PHOT_MOD.F \ + $(PHOT)/PHOTOLYSIS_ALBEDO.F $(PHOT)/PHOT_OPTICS.dat \ + $(PHOT)/SEAS_STRAT_O3_MIN.F \ $(PHOT)/twoway_rrtmg_aero_optics.F90 $(PLRISE)/delta_zs.f \ $(PLRISE)/fire_plmris.F $(PLRISE)/openlayout.F \ $(PLRISE)/oppt3d_diag.F $(PLRISE)/plmris.F $(PLRISE)/plsprd.f \ @@ -508,10 +511,12 @@ libCCTM_a_SOURCES = $(AERO)/AERO_DATA.F $(AERO)/aero_depv.F \ $(VDIFF)/opddep.F $(VDIFF)/opddep_fst.F $(VDIFF)/opddep_mos.F \ $(VDIFF)/rddepv.F $(VDIFF)/SEDIMENTATION.F $(VDIFF)/tri.F \ $(VDIFF)/VDIFF_DIAG.F $(VDIFF)/VDIFF_MAP.F \ - $(VDIFF)/vdiffproc.F $(localCCTM)/o3totcol.f \ - $(localCCTM)/vdiffacmx.F $(localCCTM)/PTMAP.F \ - $(localCCTM)/PT3D_DEFN.F $(localCCTM)/ASX_DATA_MOD.F \ - $(localCCTM)/DUST_EMIS.F + $(VDIFF)/vdiffproc.F ${localCCTM}/o3totcol.f \ + ${localCCTM}/vdiffacmx.F ${localCCTM}/PTMAP.F \ + ${localCCTM}/PT3D_DATA_MOD.F ${localCCTM}/PT3D_DEFN.F \ + ${localCCTM}/PT3D_FIRE_DEFN.F ${localCCTM}/PT3D_STKS_DEFN.F \ + ${localCCTM}/ASX_DATA_MOD.F ${localCCTM}/DUST_EMIS.F \ + ${localCCTM}/AERO_PHOTDATA.F # local version of CCTM source files localCCTM = $(builddir)/src @@ -598,8 +603,8 @@ libCCTM_a_CPPFLAGS = -DSUBST_FILES_ID=\"FILES_CTM.EXT\" \ -D_AQM_ libCCTM_a_FFLAGS = $(CCTM_FFLAGS) -I $(ICL)/filenames -I $(ICL)/const \ -I $(ICL)/emctrl -I$(ICL)/mpi -I $(top_builddir)/src/io/ioapi \ - -I $(top_builddir)/src/shr -I $(top_builddir)/src/io/aqmio \ - $(ESMF_F90COMPILEOPTS) $(ESMF_F90COMPILEPATHS) \ + -I $(top_builddir)/src/shr -I $(top_builddir)/src/io/aqmio -I \ + $(AERO) $(ESMF_F90COMPILEOPTS) $(ESMF_F90COMPILEPATHS) \ $(ESMF_F90COMPILEFREECPP) $(ESMF_F90COMPILECPPFLAGS) libCCTM_a_FCFLAGS = $(CCTM_FCFLAGS) -I $(top_builddir)/src/io/ioapi all: all-am @@ -873,8 +878,6 @@ $(PHOT)/$(am__dirstamp): $(PHOT)/$(DEPDIR)/$(am__dirstamp): @$(MKDIR_P) $(PHOT)/$(DEPDIR) @: > $(PHOT)/$(DEPDIR)/$(am__dirstamp) -$(PHOT)/libCCTM_a-AERO_PHOTDATA.$(OBJEXT): $(PHOT)/$(am__dirstamp) \ - $(PHOT)/$(DEPDIR)/$(am__dirstamp) $(PHOT)/libCCTM_a-CLOUD_OPTICS.$(OBJEXT): $(PHOT)/$(am__dirstamp) \ $(PHOT)/$(DEPDIR)/$(am__dirstamp) $(PHOT)/libCCTM_a-complex_number_module.$(OBJEXT): \ @@ -1003,29 +1006,41 @@ $(VDIFF)/libCCTM_a-VDIFF_MAP.$(OBJEXT): $(VDIFF)/$(am__dirstamp) \ $(VDIFF)/$(DEPDIR)/$(am__dirstamp) $(VDIFF)/libCCTM_a-vdiffproc.$(OBJEXT): $(VDIFF)/$(am__dirstamp) \ $(VDIFF)/$(DEPDIR)/$(am__dirstamp) -$(localCCTM)/$(am__dirstamp): - @$(MKDIR_P) $(localCCTM) - @: > $(localCCTM)/$(am__dirstamp) -$(localCCTM)/$(DEPDIR)/$(am__dirstamp): - @$(MKDIR_P) $(localCCTM)/$(DEPDIR) - @: > $(localCCTM)/$(DEPDIR)/$(am__dirstamp) -$(localCCTM)/libCCTM_a-o3totcol.$(OBJEXT): \ - $(localCCTM)/$(am__dirstamp) \ - $(localCCTM)/$(DEPDIR)/$(am__dirstamp) -$(localCCTM)/libCCTM_a-vdiffacmx.$(OBJEXT): \ - $(localCCTM)/$(am__dirstamp) \ - $(localCCTM)/$(DEPDIR)/$(am__dirstamp) -$(localCCTM)/libCCTM_a-PTMAP.$(OBJEXT): $(localCCTM)/$(am__dirstamp) \ - $(localCCTM)/$(DEPDIR)/$(am__dirstamp) -$(localCCTM)/libCCTM_a-PT3D_DEFN.$(OBJEXT): \ - $(localCCTM)/$(am__dirstamp) \ - $(localCCTM)/$(DEPDIR)/$(am__dirstamp) -$(localCCTM)/libCCTM_a-ASX_DATA_MOD.$(OBJEXT): \ - $(localCCTM)/$(am__dirstamp) \ - $(localCCTM)/$(DEPDIR)/$(am__dirstamp) -$(localCCTM)/libCCTM_a-DUST_EMIS.$(OBJEXT): \ - $(localCCTM)/$(am__dirstamp) \ - $(localCCTM)/$(DEPDIR)/$(am__dirstamp) +${localCCTM}/$(am__dirstamp): + @$(MKDIR_P) ${localCCTM} + @: > ${localCCTM}/$(am__dirstamp) +${localCCTM}/$(DEPDIR)/$(am__dirstamp): + @$(MKDIR_P) ${localCCTM}/$(DEPDIR) + @: > ${localCCTM}/$(DEPDIR)/$(am__dirstamp) +${localCCTM}/libCCTM_a-o3totcol.$(OBJEXT): \ + ${localCCTM}/$(am__dirstamp) \ + ${localCCTM}/$(DEPDIR)/$(am__dirstamp) +${localCCTM}/libCCTM_a-vdiffacmx.$(OBJEXT): \ + ${localCCTM}/$(am__dirstamp) \ + ${localCCTM}/$(DEPDIR)/$(am__dirstamp) +${localCCTM}/libCCTM_a-PTMAP.$(OBJEXT): ${localCCTM}/$(am__dirstamp) \ + ${localCCTM}/$(DEPDIR)/$(am__dirstamp) +${localCCTM}/libCCTM_a-PT3D_DATA_MOD.$(OBJEXT): \ + ${localCCTM}/$(am__dirstamp) \ + ${localCCTM}/$(DEPDIR)/$(am__dirstamp) +${localCCTM}/libCCTM_a-PT3D_DEFN.$(OBJEXT): \ + ${localCCTM}/$(am__dirstamp) \ + ${localCCTM}/$(DEPDIR)/$(am__dirstamp) +${localCCTM}/libCCTM_a-PT3D_FIRE_DEFN.$(OBJEXT): \ + ${localCCTM}/$(am__dirstamp) \ + ${localCCTM}/$(DEPDIR)/$(am__dirstamp) +${localCCTM}/libCCTM_a-PT3D_STKS_DEFN.$(OBJEXT): \ + ${localCCTM}/$(am__dirstamp) \ + ${localCCTM}/$(DEPDIR)/$(am__dirstamp) +${localCCTM}/libCCTM_a-ASX_DATA_MOD.$(OBJEXT): \ + ${localCCTM}/$(am__dirstamp) \ + ${localCCTM}/$(DEPDIR)/$(am__dirstamp) +${localCCTM}/libCCTM_a-DUST_EMIS.$(OBJEXT): \ + ${localCCTM}/$(am__dirstamp) \ + ${localCCTM}/$(DEPDIR)/$(am__dirstamp) +${localCCTM}/libCCTM_a-AERO_PHOTDATA.$(OBJEXT): \ + ${localCCTM}/$(am__dirstamp) \ + ${localCCTM}/$(DEPDIR)/$(am__dirstamp) libCCTM.a: $(libCCTM_a_OBJECTS) $(libCCTM_a_DEPENDENCIES) $(EXTRA_libCCTM_a_DEPENDENCIES) $(AM_V_at)-rm -f libCCTM.a @@ -1050,7 +1065,7 @@ mostlyclean-compile: -rm -f $(STENEX)/*.$(OBJEXT) -rm -f $(UTIL)/*.$(OBJEXT) -rm -f $(VDIFF)/*.$(OBJEXT) - -rm -f $(localCCTM)/*.$(OBJEXT) + -rm -f ${localCCTM}/*.$(OBJEXT) distclean-compile: -rm -f *.tab.c @@ -1499,12 +1514,6 @@ $(PA)/libCCTM_a-pa_update.o: $(PA)/pa_update.F $(PA)/libCCTM_a-pa_update.obj: $(PA)/pa_update.F $(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o $(PA)/libCCTM_a-pa_update.obj `if test -f '$(PA)/pa_update.F'; then $(CYGPATH_W) '$(PA)/pa_update.F'; else $(CYGPATH_W) '$(srcdir)/$(PA)/pa_update.F'; fi` -$(PHOT)/libCCTM_a-AERO_PHOTDATA.o: $(PHOT)/AERO_PHOTDATA.F - $(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o $(PHOT)/libCCTM_a-AERO_PHOTDATA.o `test -f '$(PHOT)/AERO_PHOTDATA.F' || echo '$(srcdir)/'`$(PHOT)/AERO_PHOTDATA.F - -$(PHOT)/libCCTM_a-AERO_PHOTDATA.obj: $(PHOT)/AERO_PHOTDATA.F - $(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o $(PHOT)/libCCTM_a-AERO_PHOTDATA.obj `if test -f '$(PHOT)/AERO_PHOTDATA.F'; then $(CYGPATH_W) '$(PHOT)/AERO_PHOTDATA.F'; else $(CYGPATH_W) '$(srcdir)/$(PHOT)/AERO_PHOTDATA.F'; fi` - $(PHOT)/libCCTM_a-CLOUD_OPTICS.o: $(PHOT)/CLOUD_OPTICS.F $(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o $(PHOT)/libCCTM_a-CLOUD_OPTICS.o `test -f '$(PHOT)/CLOUD_OPTICS.F' || echo '$(srcdir)/'`$(PHOT)/CLOUD_OPTICS.F @@ -1679,35 +1688,59 @@ $(VDIFF)/libCCTM_a-vdiffproc.o: $(VDIFF)/vdiffproc.F $(VDIFF)/libCCTM_a-vdiffproc.obj: $(VDIFF)/vdiffproc.F $(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o $(VDIFF)/libCCTM_a-vdiffproc.obj `if test -f '$(VDIFF)/vdiffproc.F'; then $(CYGPATH_W) '$(VDIFF)/vdiffproc.F'; else $(CYGPATH_W) '$(srcdir)/$(VDIFF)/vdiffproc.F'; fi` -$(localCCTM)/libCCTM_a-vdiffacmx.o: $(localCCTM)/vdiffacmx.F - $(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o $(localCCTM)/libCCTM_a-vdiffacmx.o `test -f '$(localCCTM)/vdiffacmx.F' || echo '$(srcdir)/'`$(localCCTM)/vdiffacmx.F +${localCCTM}/libCCTM_a-vdiffacmx.o: ${localCCTM}/vdiffacmx.F + $(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o ${localCCTM}/libCCTM_a-vdiffacmx.o `test -f '${localCCTM}/vdiffacmx.F' || echo '$(srcdir)/'`${localCCTM}/vdiffacmx.F + +${localCCTM}/libCCTM_a-vdiffacmx.obj: ${localCCTM}/vdiffacmx.F + $(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o ${localCCTM}/libCCTM_a-vdiffacmx.obj `if test -f '${localCCTM}/vdiffacmx.F'; then $(CYGPATH_W) '${localCCTM}/vdiffacmx.F'; else $(CYGPATH_W) '$(srcdir)/${localCCTM}/vdiffacmx.F'; fi` + +${localCCTM}/libCCTM_a-PTMAP.o: ${localCCTM}/PTMAP.F + $(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o ${localCCTM}/libCCTM_a-PTMAP.o `test -f '${localCCTM}/PTMAP.F' || echo '$(srcdir)/'`${localCCTM}/PTMAP.F + +${localCCTM}/libCCTM_a-PTMAP.obj: ${localCCTM}/PTMAP.F + $(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o ${localCCTM}/libCCTM_a-PTMAP.obj `if test -f '${localCCTM}/PTMAP.F'; then $(CYGPATH_W) '${localCCTM}/PTMAP.F'; else $(CYGPATH_W) '$(srcdir)/${localCCTM}/PTMAP.F'; fi` -$(localCCTM)/libCCTM_a-vdiffacmx.obj: $(localCCTM)/vdiffacmx.F - $(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o $(localCCTM)/libCCTM_a-vdiffacmx.obj `if test -f '$(localCCTM)/vdiffacmx.F'; then $(CYGPATH_W) '$(localCCTM)/vdiffacmx.F'; else $(CYGPATH_W) '$(srcdir)/$(localCCTM)/vdiffacmx.F'; fi` +${localCCTM}/libCCTM_a-PT3D_DATA_MOD.o: ${localCCTM}/PT3D_DATA_MOD.F + $(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o ${localCCTM}/libCCTM_a-PT3D_DATA_MOD.o `test -f '${localCCTM}/PT3D_DATA_MOD.F' || echo '$(srcdir)/'`${localCCTM}/PT3D_DATA_MOD.F -$(localCCTM)/libCCTM_a-PTMAP.o: $(localCCTM)/PTMAP.F - $(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o $(localCCTM)/libCCTM_a-PTMAP.o `test -f '$(localCCTM)/PTMAP.F' || echo '$(srcdir)/'`$(localCCTM)/PTMAP.F +${localCCTM}/libCCTM_a-PT3D_DATA_MOD.obj: ${localCCTM}/PT3D_DATA_MOD.F + $(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o ${localCCTM}/libCCTM_a-PT3D_DATA_MOD.obj `if test -f '${localCCTM}/PT3D_DATA_MOD.F'; then $(CYGPATH_W) '${localCCTM}/PT3D_DATA_MOD.F'; else $(CYGPATH_W) '$(srcdir)/${localCCTM}/PT3D_DATA_MOD.F'; fi` -$(localCCTM)/libCCTM_a-PTMAP.obj: $(localCCTM)/PTMAP.F - $(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o $(localCCTM)/libCCTM_a-PTMAP.obj `if test -f '$(localCCTM)/PTMAP.F'; then $(CYGPATH_W) '$(localCCTM)/PTMAP.F'; else $(CYGPATH_W) '$(srcdir)/$(localCCTM)/PTMAP.F'; fi` +${localCCTM}/libCCTM_a-PT3D_DEFN.o: ${localCCTM}/PT3D_DEFN.F + $(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o ${localCCTM}/libCCTM_a-PT3D_DEFN.o `test -f '${localCCTM}/PT3D_DEFN.F' || echo '$(srcdir)/'`${localCCTM}/PT3D_DEFN.F -$(localCCTM)/libCCTM_a-PT3D_DEFN.o: $(localCCTM)/PT3D_DEFN.F - $(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o $(localCCTM)/libCCTM_a-PT3D_DEFN.o `test -f '$(localCCTM)/PT3D_DEFN.F' || echo '$(srcdir)/'`$(localCCTM)/PT3D_DEFN.F +${localCCTM}/libCCTM_a-PT3D_DEFN.obj: ${localCCTM}/PT3D_DEFN.F + $(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o ${localCCTM}/libCCTM_a-PT3D_DEFN.obj `if test -f '${localCCTM}/PT3D_DEFN.F'; then $(CYGPATH_W) '${localCCTM}/PT3D_DEFN.F'; else $(CYGPATH_W) '$(srcdir)/${localCCTM}/PT3D_DEFN.F'; fi` -$(localCCTM)/libCCTM_a-PT3D_DEFN.obj: $(localCCTM)/PT3D_DEFN.F - $(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o $(localCCTM)/libCCTM_a-PT3D_DEFN.obj `if test -f '$(localCCTM)/PT3D_DEFN.F'; then $(CYGPATH_W) '$(localCCTM)/PT3D_DEFN.F'; else $(CYGPATH_W) '$(srcdir)/$(localCCTM)/PT3D_DEFN.F'; fi` +${localCCTM}/libCCTM_a-PT3D_FIRE_DEFN.o: ${localCCTM}/PT3D_FIRE_DEFN.F + $(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o ${localCCTM}/libCCTM_a-PT3D_FIRE_DEFN.o `test -f '${localCCTM}/PT3D_FIRE_DEFN.F' || echo '$(srcdir)/'`${localCCTM}/PT3D_FIRE_DEFN.F -$(localCCTM)/libCCTM_a-ASX_DATA_MOD.o: $(localCCTM)/ASX_DATA_MOD.F - $(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o $(localCCTM)/libCCTM_a-ASX_DATA_MOD.o `test -f '$(localCCTM)/ASX_DATA_MOD.F' || echo '$(srcdir)/'`$(localCCTM)/ASX_DATA_MOD.F +${localCCTM}/libCCTM_a-PT3D_FIRE_DEFN.obj: ${localCCTM}/PT3D_FIRE_DEFN.F + $(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o ${localCCTM}/libCCTM_a-PT3D_FIRE_DEFN.obj `if test -f '${localCCTM}/PT3D_FIRE_DEFN.F'; then $(CYGPATH_W) '${localCCTM}/PT3D_FIRE_DEFN.F'; else $(CYGPATH_W) '$(srcdir)/${localCCTM}/PT3D_FIRE_DEFN.F'; fi` -$(localCCTM)/libCCTM_a-ASX_DATA_MOD.obj: $(localCCTM)/ASX_DATA_MOD.F - $(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o $(localCCTM)/libCCTM_a-ASX_DATA_MOD.obj `if test -f '$(localCCTM)/ASX_DATA_MOD.F'; then $(CYGPATH_W) '$(localCCTM)/ASX_DATA_MOD.F'; else $(CYGPATH_W) '$(srcdir)/$(localCCTM)/ASX_DATA_MOD.F'; fi` +${localCCTM}/libCCTM_a-PT3D_STKS_DEFN.o: ${localCCTM}/PT3D_STKS_DEFN.F + $(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o ${localCCTM}/libCCTM_a-PT3D_STKS_DEFN.o `test -f '${localCCTM}/PT3D_STKS_DEFN.F' || echo '$(srcdir)/'`${localCCTM}/PT3D_STKS_DEFN.F -$(localCCTM)/libCCTM_a-DUST_EMIS.o: $(localCCTM)/DUST_EMIS.F - $(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o $(localCCTM)/libCCTM_a-DUST_EMIS.o `test -f '$(localCCTM)/DUST_EMIS.F' || echo '$(srcdir)/'`$(localCCTM)/DUST_EMIS.F +${localCCTM}/libCCTM_a-PT3D_STKS_DEFN.obj: ${localCCTM}/PT3D_STKS_DEFN.F + $(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o ${localCCTM}/libCCTM_a-PT3D_STKS_DEFN.obj `if test -f '${localCCTM}/PT3D_STKS_DEFN.F'; then $(CYGPATH_W) '${localCCTM}/PT3D_STKS_DEFN.F'; else $(CYGPATH_W) '$(srcdir)/${localCCTM}/PT3D_STKS_DEFN.F'; fi` -$(localCCTM)/libCCTM_a-DUST_EMIS.obj: $(localCCTM)/DUST_EMIS.F - $(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o $(localCCTM)/libCCTM_a-DUST_EMIS.obj `if test -f '$(localCCTM)/DUST_EMIS.F'; then $(CYGPATH_W) '$(localCCTM)/DUST_EMIS.F'; else $(CYGPATH_W) '$(srcdir)/$(localCCTM)/DUST_EMIS.F'; fi` +${localCCTM}/libCCTM_a-ASX_DATA_MOD.o: ${localCCTM}/ASX_DATA_MOD.F + $(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o ${localCCTM}/libCCTM_a-ASX_DATA_MOD.o `test -f '${localCCTM}/ASX_DATA_MOD.F' || echo '$(srcdir)/'`${localCCTM}/ASX_DATA_MOD.F + +${localCCTM}/libCCTM_a-ASX_DATA_MOD.obj: ${localCCTM}/ASX_DATA_MOD.F + $(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o ${localCCTM}/libCCTM_a-ASX_DATA_MOD.obj `if test -f '${localCCTM}/ASX_DATA_MOD.F'; then $(CYGPATH_W) '${localCCTM}/ASX_DATA_MOD.F'; else $(CYGPATH_W) '$(srcdir)/${localCCTM}/ASX_DATA_MOD.F'; fi` + +${localCCTM}/libCCTM_a-DUST_EMIS.o: ${localCCTM}/DUST_EMIS.F + $(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o ${localCCTM}/libCCTM_a-DUST_EMIS.o `test -f '${localCCTM}/DUST_EMIS.F' || echo '$(srcdir)/'`${localCCTM}/DUST_EMIS.F + +${localCCTM}/libCCTM_a-DUST_EMIS.obj: ${localCCTM}/DUST_EMIS.F + $(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o ${localCCTM}/libCCTM_a-DUST_EMIS.obj `if test -f '${localCCTM}/DUST_EMIS.F'; then $(CYGPATH_W) '${localCCTM}/DUST_EMIS.F'; else $(CYGPATH_W) '$(srcdir)/${localCCTM}/DUST_EMIS.F'; fi` + +${localCCTM}/libCCTM_a-AERO_PHOTDATA.o: ${localCCTM}/AERO_PHOTDATA.F + $(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o ${localCCTM}/libCCTM_a-AERO_PHOTDATA.o `test -f '${localCCTM}/AERO_PHOTDATA.F' || echo '$(srcdir)/'`${localCCTM}/AERO_PHOTDATA.F + +${localCCTM}/libCCTM_a-AERO_PHOTDATA.obj: ${localCCTM}/AERO_PHOTDATA.F + $(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o ${localCCTM}/libCCTM_a-AERO_PHOTDATA.obj `if test -f '${localCCTM}/AERO_PHOTDATA.F'; then $(CYGPATH_W) '${localCCTM}/AERO_PHOTDATA.F'; else $(CYGPATH_W) '$(srcdir)/${localCCTM}/AERO_PHOTDATA.F'; fi` .F.f: $(F77COMPILE) -F $< @@ -1909,11 +1942,11 @@ $(UTIL)/libCCTM_a-get_envlist.o: $(UTIL)/get_envlist.f $(UTIL)/libCCTM_a-get_envlist.obj: $(UTIL)/get_envlist.f $(AM_V_F77)$(F77) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o $(UTIL)/libCCTM_a-get_envlist.obj `if test -f '$(UTIL)/get_envlist.f'; then $(CYGPATH_W) '$(UTIL)/get_envlist.f'; else $(CYGPATH_W) '$(srcdir)/$(UTIL)/get_envlist.f'; fi` -$(localCCTM)/libCCTM_a-o3totcol.o: $(localCCTM)/o3totcol.f - $(AM_V_F77)$(F77) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o $(localCCTM)/libCCTM_a-o3totcol.o `test -f '$(localCCTM)/o3totcol.f' || echo '$(srcdir)/'`$(localCCTM)/o3totcol.f +${localCCTM}/libCCTM_a-o3totcol.o: ${localCCTM}/o3totcol.f + $(AM_V_F77)$(F77) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o ${localCCTM}/libCCTM_a-o3totcol.o `test -f '${localCCTM}/o3totcol.f' || echo '$(srcdir)/'`${localCCTM}/o3totcol.f -$(localCCTM)/libCCTM_a-o3totcol.obj: $(localCCTM)/o3totcol.f - $(AM_V_F77)$(F77) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o $(localCCTM)/libCCTM_a-o3totcol.obj `if test -f '$(localCCTM)/o3totcol.f'; then $(CYGPATH_W) '$(localCCTM)/o3totcol.f'; else $(CYGPATH_W) '$(srcdir)/$(localCCTM)/o3totcol.f'; fi` +${localCCTM}/libCCTM_a-o3totcol.obj: ${localCCTM}/o3totcol.f + $(AM_V_F77)$(F77) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o ${localCCTM}/libCCTM_a-o3totcol.obj `if test -f '${localCCTM}/o3totcol.f'; then $(CYGPATH_W) '${localCCTM}/o3totcol.f'; else $(CYGPATH_W) '$(srcdir)/${localCCTM}/o3totcol.f'; fi` ID: $(am__tagged_files) $(am__define_uniq_tagged_files); mkid -fID $$unique @@ -2027,6 +2060,8 @@ clean-generic: distclean-generic: -test -z "$(CONFIG_CLEAN_FILES)" || rm -f $(CONFIG_CLEAN_FILES) -test . = "$(srcdir)" || test -z "$(CONFIG_CLEAN_VPATH_FILES)" || rm -f $(CONFIG_CLEAN_VPATH_FILES) + -rm -f ${localCCTM}/$(DEPDIR)/$(am__dirstamp) + -rm -f ${localCCTM}/$(am__dirstamp) -test -z "$(AERO)/$(DEPDIR)/$(am__dirstamp)" || rm -f $(AERO)/$(DEPDIR)/$(am__dirstamp) -test -z "$(AERO)/$(am__dirstamp)" || rm -f $(AERO)/$(am__dirstamp) -test -z "$(BIOG)/$(DEPDIR)/$(am__dirstamp)" || rm -f $(BIOG)/$(DEPDIR)/$(am__dirstamp) @@ -2059,8 +2094,6 @@ distclean-generic: -test -z "$(UTIL)/$(am__dirstamp)" || rm -f $(UTIL)/$(am__dirstamp) -test -z "$(VDIFF)/$(DEPDIR)/$(am__dirstamp)" || rm -f $(VDIFF)/$(DEPDIR)/$(am__dirstamp) -test -z "$(VDIFF)/$(am__dirstamp)" || rm -f $(VDIFF)/$(am__dirstamp) - -test -z "$(localCCTM)/$(DEPDIR)/$(am__dirstamp)" || rm -f $(localCCTM)/$(DEPDIR)/$(am__dirstamp) - -test -z "$(localCCTM)/$(am__dirstamp)" || rm -f $(localCCTM)/$(am__dirstamp) maintainer-clean-generic: @echo "This command is intended for maintainers to use" @@ -2423,12 +2456,6 @@ $(libPA)pa_update.$(OBJEXT) : $(ICL)/const/CONST.EXT $(ICL)/emctrl/EMISPRM.EXT $ $(libGRID)PAGRD_DEFN.$(OBJEXT) $(libUTIL)UTILIO_DEFN.$(OBJEXT) # phot -$(libPHOT)AERO_PHOTDATA.$(OBJEXT) : $(ICL)/const/CONST.EXT \ - $(libAERO)AERO_DATA.$(OBJEXT) $(libAERO)AEROMET_DATA.$(OBJEXT) \ - $(libSPCS)CGRID_SPCS.$(OBJEXT) $(libPHOT)twoway_rrtmg_aero_optics.$(OBJEXT) \ - $(libPHOT)CSQY_DATA.$(OBJEXT) $(libGRID)GRID_CONF.$(OBJEXT) \ - $(libMECHS)RXNS_DATA_MODULE.$(OBJEXT) $(libAERO)SOA_DEFN.$(OBJEXT) \ - $(libUTIL)UTILIO_DEFN.$(OBJEXT) $(libPHOT)CLOUD_OPTICS.$(OBJEXT) : \ $(libPHOT)CSQY_DATA.$(OBJEXT) $(libGRID)GRID_CONF.$(OBJEXT) \ $(libUTIL)UTILIO_DEFN.$(OBJEXT) @@ -2440,7 +2467,7 @@ $(libPHOT)opphot.$(OBJEXT) : $(ICL)/filenames/FILES_CTM.EXT \ $(libPHOT)PHOT_MOD.$(OBJEXT) $(libMECHS)RXNS_DATA_MODULE.$(OBJEXT) \ $(libUTIL)UTILIO_DEFN.$(OBJEXT) $(libPHOT)phot.$(OBJEXT) : $(ICL)/filenames/FILES_CTM.EXT \ - $(libAERO)AERO_DATA.$(OBJEXT) $(libPHOT)AERO_PHOTDATA.$(OBJEXT) \ + $(libAERO)AERO_DATA.$(OBJEXT) $(liblocalCCTM)AERO_PHOTDATA.$(OBJEXT) \ $(libSPCS)CGRID_SPCS.$(OBJEXT) $(libPHOT)CLOUD_OPTICS.$(OBJEXT) \ $(libSTENEX)noop_modules.$(OBJEXT) $(libGRID)PCGRID_DEFN.$(OBJEXT) \ $(libPHOT)PHOT_MET_DATA.$(OBJEXT) $(libPHOT)PHOT_MOD.$(OBJEXT) \ @@ -2450,7 +2477,7 @@ $(libPHOT)PHOT_MET_DATA.$(OBJEXT) : $(ICL)/const/CONST.EXT $(ICL)/filenames/FILE $(libPHOT)CLOUD_OPTICS.$(OBJEXT) $(libGRID)GRID_CONF.$(OBJEXT) \ $(libUTIL)UTILIO_DEFN.$(OBJEXT) $(libPHOT)PHOT_MOD.$(OBJEXT) : $(ICL)/const/CONST.EXT \ - $(libPHOT)AERO_PHOTDATA.$(OBJEXT) $(libPHOT)CLOUD_OPTICS.$(OBJEXT) \ + $(liblocalCCTM)AERO_PHOTDATA.$(OBJEXT) $(libPHOT)CLOUD_OPTICS.$(OBJEXT) \ $(libPHOT)CSQY_DATA.$(OBJEXT) $(libMECHS)RXNS_DATA_MODULE.$(OBJEXT) \ $(libUTIL)UTILIO_DEFN.$(OBJEXT) $(libPHOT)PHOTOLYSIS_ALBEDO.$(OBJEXT) : $(ICL)/filenames/FILES_CTM.EXT \ @@ -2553,10 +2580,19 @@ $(liblocalCCTM)vdiffacmx.$(OBJEXT) : $(ICL)/filenames/FILES_CTM.EXT \ $(libVDIFF)VDIFF_MAP.$(OBJEXT) $(liblocalCCTM)PTMAP.$(OBJEXT) : $(libAERO)AERO_DATA.$(OBJEXT) \ $(libSPCS)CGRID_SPCS.$(OBJEXT) $(libUTIL)UTILIO_DEFN.$(OBJEXT) -$(liblocalCCTM)PT3D_DEFN.$(OBJEXT) : $(libAERO)AERO_DATA.$(OBJEXT) \ +$(liblocalCCTM)PT3D_DEFN.$(OBJEXT) : $(libGRID)GRID_CONF.$(OBJEXT) \ + $(liblocalCCTM)PTMAP.$(OBJEXT) $(liblocalCCTM)PT3D_FIRE_DEFN.$(OBJEXT) \ + $(liblocalCCTM)PT3D_STKS_DEFN.$(OBJEXT) +$(liblocalCCTM)PT3D_FIRE_DEFN.$(OBJEXT) : $(libAERO)AERO_DATA.$(OBJEXT) \ $(libGRID)GRID_CONF.$(OBJEXT) $(libSPCS)CGRID_SPCS.$(OBJEXT) \ + $(liblocalCCTM)PT3D_DATA_MOD.$(OBJEXT) \ $(liblocalCCTM)PTMAP.$(OBJEXT) $(libMECHS)RXNS_DATA_MODULE.$(OBJEXT) \ $(libEMIS)STK_EMIS.$(OBJEXT) $(libUTIL)UTILIO_DEFN.$(OBJEXT) +$(liblocalCCTM)PT3D_STKS_DEFN.$(OBJEXT) : $(libAERO)AERO_DATA.$(OBJEXT) \ + $(liblocalCCTM)ASX_DATA_MOD.$(OBJEXT) \ + $(libGRID)GRID_CONF.$(OBJEXT) $(libSPCS)CGRID_SPCS.$(OBJEXT) \ + $(liblocalCCTM)PT3D_DATA_MOD.$(OBJEXT) $(liblocalCCTM)PTMAP.$(OBJEXT) \ + $(libEMIS)STK_EMIS.$(OBJEXT) $(libUTIL)UTILIO_DEFN.$(OBJEXT) $(liblocalCCTM)ASX_DATA_MOD.$(OBJEXT) : $(ICL)/const/CONST.EXT $(ICL)/filenames/FILES_CTM.EXT $(ICL)/mpi/PE_COMM.EXT \ $(libDEPV)DEPVVARS.$(OBJEXT) $(libGRID)GRID_CONF.$(OBJEXT) \ $(libDEPV)LSM_MOD.$(OBJEXT) $(libSTENEX)noop_modules.$(OBJEXT) \ @@ -2565,6 +2601,12 @@ $(liblocalCCTM)DUST_EMIS.$(OBJEXT) : $(ICL)/const/CONST.EXT $(ICL)/filenames/FIL $(libAERO)AERO_DATA.$(OBJEXT) $(liblocalCCTM)ASX_DATA_MOD.$(OBJEXT) \ $(libGRID)GRID_CONF.$(OBJEXT) $(libGRID)HGRD_DEFN.$(OBJEXT) \ $(libEMIS)LUS_DEFN.$(OBJEXT) $(libUTIL)UTILIO_DEFN.$(OBJEXT) +$(liblocalCCTM)AERO_PHOTDATA.$(OBJEXT) : $(ICL)/const/CONST.EXT \ + $(libAERO)AERO_DATA.$(OBJEXT) $(libAERO)AEROMET_DATA.$(OBJEXT) \ + $(libSPCS)CGRID_SPCS.$(OBJEXT) $(libPHOT)twoway_rrtmg_aero_optics.$(OBJEXT) \ + $(libPHOT)CSQY_DATA.$(OBJEXT) $(libGRID)GRID_CONF.$(OBJEXT) \ + $(libMECHS)RXNS_DATA_MODULE.$(OBJEXT) $(libAERO)SOA_DEFN.$(OBJEXT) \ + $(libUTIL)UTILIO_DEFN.$(OBJEXT) # Tell versions [3.59,3.63) of GNU make to not export all variables. # Otherwise a system limit (for SysV at least) may be exceeded. diff --git a/src/model/src/ASX_DATA_MOD.F b/src/model/src/ASX_DATA_MOD.F index 223d12b9..16b30957 100755 --- a/src/model/src/ASX_DATA_MOD.F +++ b/src/model/src/ASX_DATA_MOD.F @@ -133,6 +133,7 @@ Module ASX_DATA_MOD !> 3-D meteorological fields: Real, Allocatable :: KZMIN ( :,:,: ) ! minimum Kz [m**2/s] Real, Allocatable :: PRES ( :,:,: ) ! layer 1 pressure [Pa] + Real, Allocatable :: PRESF ( :,:,: ) ! full layer pressure [Pa] Real, Allocatable :: QV ( :,:,: ) ! water vapor mixing ratio Real, Allocatable :: QC ( :,:,: ) ! cloud water mixing ratio Real, Allocatable :: THETAV ( :,:,: ) ! potential temp @@ -144,6 +145,8 @@ Module ASX_DATA_MOD Real, Allocatable :: RJACM ( :,:,: ) ! reciprocal mid-layer Jacobian Real, Allocatable :: RJACF ( :,:,: ) ! reciprocal full-layer Jacobian Real, Allocatable :: RRHOJ ( :,:,: ) ! reciprocal density X Jacobian + Real, Allocatable :: UWINDA ( :,:,: ) ! [m/s] + Real, Allocatable :: VWINDA ( :,:,: ) ! [m/s] End Type MET_Type Type :: GRID_Type @@ -496,10 +499,13 @@ Subroutine INIT_MET ( JDATE, JTIME, MOSAIC, ABFLUX, HGBIDI ) & Met_Data%CONVCT ( NCOLS,NROWS ), & Met_Data%PBL ( NCOLS,NROWS ), & Met_Data%NACL_EMIS( NCOLS,NROWS ), + & Met_Data%UWINDA ( NCOLS,NROWS,NLAYS ), + & Met_Data%VWINDA ( NCOLS,NROWS,NLAYS ), & Met_Data%UWIND ( NCOLS+1,NROWS+1,NLAYS ), & Met_Data%VWIND ( NCOLS+1,NROWS+1,NLAYS ), & Met_Data%KZMIN ( NCOLS,NROWS,NLAYS ), & Met_Data%PRES ( NCOLS,NROWS,NLAYS ), + & Met_Data%PRESF ( NCOLS,NROWS,1:NLAYS+1 ), & Met_Data%QV ( NCOLS,NROWS,NLAYS ), & Met_Data%QC ( NCOLS,NROWS,NLAYS ), & Met_Data%THETAV ( NCOLS,NROWS,NLAYS ), @@ -963,6 +969,14 @@ Subroutine GET_MET ( JDATE, JTIME, TSTEP, MOSAIC, ABFLUX, HGBIDI ) Call M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT1 ) End If + VNAME = 'PRESF' + If ( .Not. INTERPX( MET_CRO_3D, VNAME, PNAME, + & STRTCOLMC3,ENDCOLMC3, STRTROWMC3,ENDROWMC3, 1,NLAYS+1, + & JDATE, JTIME, Met_Data%PRESF ) ) Then + XMSG = MSG1 // TRIM( VNAME ) // ' from ' // MET_CRO_3D + Call M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT1 ) + End If + VNAME = 'ZF' If ( .Not. INTERPX( MET_CRO_3D, VNAME, PNAME, & STRTCOLMC3,ENDCOLMC3, STRTROWMC3,ENDROWMC3, 1,NLAYS, @@ -1035,6 +1049,22 @@ Subroutine GET_MET ( JDATE, JTIME, TSTEP, MOSAIC, ABFLUX, HGBIDI ) CALL M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT1 ) END IF + VNAME = 'UWINDA' + If ( .Not. INTERPX( MET_CRO_2D, VNAME, PNAME, + & STRTCOLMC2,ENDCOLMC2, STRTROWMC2,ENDROWMC2, 1,NLAYS, + & JDATE, JTIME, Met_Data%UWINDA ) ) Then + XMSG = MSG1 // TRIM( VNAME ) // ' from ' // MET_CRO_2D + Call M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT1 ) + End If + + VNAME = 'VWINDA' + If ( .Not. INTERPX( MET_CRO_2D, VNAME, PNAME, + & STRTCOLMC2,ENDCOLMC2, STRTROWMC2,ENDROWMC2, 1,NLAYS, + & JDATE, JTIME, Met_Data%VWINDA ) ) Then + XMSG = MSG1 // TRIM( VNAME ) // ' from ' // MET_CRO_2D + Call M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT1 ) + End If + C-------------------------------- MET_CRO_2D -------------------------------- C Vegetation and surface vars VNAME = 'LAI' diff --git a/src/model/src/PT3D_DATA_MOD.F b/src/model/src/PT3D_DATA_MOD.F new file mode 100644 index 00000000..086b81c6 --- /dev/null +++ b/src/model/src/PT3D_DATA_MOD.F @@ -0,0 +1,14 @@ + MODULE PT3D_DATA_MOD + + IMPLICIT NONE + + LOGICAL :: PT3DEMIS ! flag in-lining plume rise + INTEGER :: NPTGRPS ! no. pt src input file groups + + REAL, ALLOCATABLE :: VDEMIS_PT( :,:,:,: ) ! 3D pt src non-PM emis + REAL, ALLOCATABLE :: VDEMIS_PT_FIRE( :,:,:,: ) ! 3D pt src non-PM emis + REAL, ALLOCATABLE :: PMEMIS_PT( :,:,:,: ) ! 3D pt src PM emis + + CHARACTER( * ), PARAMETER :: CTM_PT3DEMIS = 'CTM_PT3DEMIS' ! env var for in-line 3d pt src emis + + END MODULE PT3D_DATA_MOD diff --git a/src/model/src/PT3D_DEFN.F b/src/model/src/PT3D_DEFN.F index 64e90fcf..0a7163db 100644 --- a/src/model/src/PT3D_DEFN.F +++ b/src/model/src/PT3D_DEFN.F @@ -1,94 +1,50 @@ MODULE PT3D_DEFN + USE PT3D_DATA_MOD + USE PT3D_FIRE_DEFN + USE PT3D_STKS_DEFN + IMPLICIT NONE - LOGICAL, SAVE :: PT3DEMIS ! flag in-lining plume rise - INTEGER, SAVE :: NPTGRPS ! no. pt src input file groups - REAL, ALLOCATABLE, SAVE :: VDEMIS_PT( :,:,:,: ) ! 3D pt src non-PM emis - REAL, ALLOCATABLE, SAVE :: VDEMIS_PT_FIRE( :,:,:,: ) ! 3D pt src non-PM emis - REAL, ALLOCATABLE, SAVE :: PMEMIS_PT( :,:,:,: ) ! 3D pt src PM emis - REAL, ALLOCATABLE, SAVE :: VFRAC( :,:,: ) ! vertical fraction - REAL, ALLOCATABLE, SAVE :: BUFFER( : ) ! emission buffer + PRIVATE PUBLIC PT3DEMIS, NPTGRPS, VDEMIS_PT, VDEMIS_PT_FIRE, PMEMIS_PT, & PT3D_INIT, GET_PT3D_EMIS - PRIVATE - -C Stack emissions for sources within domain - INTEGER, SAVE :: EMLYRS ! no. of emis layers - INTEGER, SAVE :: PM_EMLYRS ! no. of emis layers for PM - -C Gas-phase Namelist emissions factors - REAL, ALLOCATABLE, SAVE :: EMIS_FAC( : ) ! per species scaling factor - - CHARACTER( 240 ) :: XMSG = ' ' - -C Emission type - CHARACTER( * ), PARAMETER :: ETYPE = 'gbbepx' - -C Species name for Fire Radiative Power - CHARACTER( * ), PARAMETER :: EMFRP = 'FRP' CONTAINS -C======================================================================= - FUNCTION PT3D_INIT ( N_SPC_EMIS, EMLAYS, JDATE, JTIME, TSTEP ) - & RESULT ( SUCCESS ) + & RESULT ( SUCCESS ) - USE AQM_EMIS_MOD - USE GRID_CONF ! horizontal & vertical domain specifications - USE CGRID_SPCS ! CGRID mechanism species - USE STK_EMIS, ONLY : STKSPC ! hourly point source emissions - USE PTMAP ! defines pt src species mapping to VDEMIS* arrays + USE GRID_CONF ! horizontal & vertical domain specifications + USE PTMAP + USE STK_EMIS, ONLY : STKSPC ! hourly point source emissions USE UTILIO_DEFN IMPLICIT NONE -C Includes: -C INCLUDE SUBST_FILES_ID ! file name parameters (for CTM_PT3D_DIAG) - -C Arguments: INTEGER, INTENT( IN ) :: N_SPC_EMIS ! total no. of model emissions species INTEGER, INTENT( IN ) :: EMLAYS ! number of emissions layers INTEGER, INTENT( IN ) :: JDATE ! Julian date (YYYYDDD) INTEGER, INTENT( IN ) :: JTIME ! time (HHMMSS) INTEGER, INTENT( IN ) :: TSTEP ! output time step - LOGICAL SUCCESS - -C External Functions: - INTEGER, EXTERNAL :: SETUP_LOGDEV - -C Parameters: - -C Local Variables: - CHARACTER( 16 ), SAVE :: CTM_PT3DEMIS = 'CTM_PT3DEMIS ' ! env var for in-line - ! 3d pt src emis - CHARACTER( 16 ) :: PNAME = 'PT3D_INIT ' ! procedure name - CHARACTER( 16 ) :: VNAME ! variable name buffer - CHARACTER( 16 ), SAVE, ALLOCATABLE :: STKGNAME( : ) ! stack groups file name - - INTEGER IOS ! i/o and allocate memory status - - INTEGER LOGDEV - - INTEGER IDX - INTEGER N, NSPC, NSPC1, NSPC2, NSPC3 - INTEGER S, S_OFFSET, V + ! -- local variables + INTEGER :: IOS + INTEGER :: N, S, V INTEGER, ALLOCATABLE :: MAP( : ) + CHARACTER( 240 ) :: XMSG = ' ' + CHARACTER( 16 ) :: PNAME = 'PT3D_INIT ' ! procedure name - TYPE( AQM_INTERNAL_EMIS_TYPE ), POINTER :: EM + LOGICAL :: SUCCESS -C----------------------------------------------------------------------- + ! -- begin SUCCESS = .TRUE. - LOGDEV = SETUP_LOGDEV() - -C In-line 3D point source emissions? + ! -- In-line 3D point source emissions? - PT3DEMIS = ENVYN( 'CTM_PT3DEMIS', + PT3DEMIS = ENVYN( CTM_PT3DEMIS, & 'Flag for in-line 3d point source emissions', & .FALSE., IOS ) @@ -99,242 +55,67 @@ FUNCTION PT3D_INIT ( N_SPC_EMIS, EMLAYS, JDATE, JTIME, TSTEP ) RETURN END IF -C check if emissions are being provided - - EM => AQM_EMIS_GET( ETYPE ) - IF ( .NOT.ASSOCIATED( EM ) ) RETURN - -C get number of different file groups (sectors) - - NPTGRPS = 0 - -C create point source internal array + ! -- merge PM maps from fire and point-source emissions - ALLOCATE( STKSPC( NPTGRPS ), STAT = IOS ) - CALL CHECKMEM( IOS, 'STKSPC', PNAME ) + IF ( .NOT. PTMAP_INIT( ) ) THEN + XMSG = 'Could not merge point source mappings' + CALL M3WARN ( PNAME, JDATE, JTIME, XMSG ) + SUCCESS = .FALSE.; RETURN + END IF -C set number of emissions layers depending on whether plumerise is on + ! -- initialize emission types - SELECT CASE ( TRIM( EM % PLUMERISE ) ) - CASE ("sofiev") - EMLYRS = NLAYS - PM_EMLYRS = NLAYS - CASE DEFAULT - EMLYRS = 1 - PM_EMLYRS = 1 - END SELECT - -C get point source emission mapping + IF ( .NOT. PT3D_FIRE_INIT ( N_SPC_EMIS, EMLAYS, JDATE, JTIME ) ) THEN + XMSG = 'Could not initialize fire emissions' + CALL M3WARN ( PNAME, JDATE, JTIME, XMSG ) + SUCCESS = .FALSE.; RETURN + END IF - IF ( .NOT. PTMAP_INIT( ) ) THEN - XMSG = 'Could not get point source mappings' + IF ( .NOT. PT3D_STKS_INIT ( ) ) THEN + XMSG = 'Could not initialize point-source emissions' CALL M3WARN ( PNAME, JDATE, JTIME, XMSG ) SUCCESS = .FALSE.; RETURN END IF -C allocate emission arrays + ! -- allocate emission arrays - ALLOCATE ( VDEMIS_PT( NCOLS,NROWS,EMLYRS,N_SPC_PTEM ), STAT = IOS ) + ALLOCATE ( VDEMIS_PT( NCOLS, NROWS, EMLAYS, N_SPC_PTEM ), STAT = IOS ) CALL CHECKMEM( IOS, 'VDEMIS_PT', PNAME ) VDEMIS_PT = 0.0 ! array assignment - ALLOCATE ( VDEMIS_PT_FIRE( NCOLS,NROWS,EMLYRS,N_SPC_PTEM ), STAT = IOS ) + ALLOCATE ( VDEMIS_PT_FIRE( NCOLS, NROWS, EMLAYS, N_SPC_PTEM ), STAT = IOS ) CALL CHECKMEM( IOS, 'VDEMIS_PT_FIRE', PNAME ) VDEMIS_PT_FIRE = 0.0 ! array assignment - ALLOCATE ( PMEMIS_PT( NCOLS,NROWS,EMLYRS,N_SPC_PTPM ), STAT = IOS ) + ALLOCATE ( PMEMIS_PT( NCOLS, NROWS, EMLAYS, N_SPC_PTPM ), STAT = IOS ) CALL CHECKMEM( IOS, 'PMEMIS_PT', PNAME ) PMEMIS_PT = 0.0 ! array - ALLOCATE ( VFRAC( NCOLS,NROWS,EMLYRS ), STAT = IOS ) - CALL CHECKMEM( IOS, 'VFRAC', PNAME ) - VFRAC = 1.0 ! array + ! -- get number of different file groups (sectors) - ALLOCATE ( BUFFER( NCOLS * NROWS ), STAT = IOS ) - CALL CHECKMEM( IOS, 'BUFFER', PNAME ) - BUFFER = 0.0 ! array + NPTGRPS = 0 - SUCCESS = .TRUE.; RETURN + ! -- create point source internal array + + ALLOCATE( STKSPC( NPTGRPS ), STAT = IOS ) + CALL CHECKMEM( IOS, 'STKSPC', PNAME ) END FUNCTION PT3D_INIT -C======================================================================= SUBROUTINE GET_PT3D_EMIS ( JDATE, JTIME, TSTEP ) -! Revision History. -! Aug 12, 15 D. Wong: added code to handle parallel I/O implementation - -C----------------------------------------------------------------------- - -C Time step part of laypoint - - USE AQM_EMIS_MOD - USE AQM_FIRES_MOD - USE AQM_RC_MOD - USE RXNS_DATA, ONLY : MECHNAME !Get Chemical Mechanism Name - USE GRID_CONF ! horizontal & vertical domain specifications - USE CGRID_SPCS ! CGRID mechanism species - USE AERO_DATA, ONLY : PMEM_MAP_NAME - USE PTMAP ! defines pt src species mapping to VDEMIS* arrays - USE UTILIO_DEFN - - IMPLICIT NONE - -C Includes: -C INCLUDE SUBST_CONST ! physical and mathematical constants -C INCLUDE SUBST_FILES_ID ! file name parameters (for CTM_PT3D_DIAG) - -C Arguments: INTEGER, INTENT( IN ) :: JDATE, JTIME INTEGER, INTENT( IN ) :: TSTEP( 3 ) -C Parameters: - -C External functions: - INTEGER, EXTERNAL :: SETUP_LOGDEV - -C Local variables: - CHARACTER( 16 ) :: PNAME = 'GET_PT3D_EMIS ' ! procedure name - CHARACTER( 16 ) :: VNAME ! variable name buffer - - INTEGER IOS ! i/o and allocate memory status - INTEGER L, S, V ! counters - INTEGER C, R, K, N - INTEGER LOCALRC - - LOGICAL :: IS_NOT_NVPOA, SAVE_POC - - LOGICAL, SAVE :: FIRSTIME = .TRUE. - INTEGER, SAVE :: LOGDEV - - TYPE( AQM_INTERNAL_EMIS_TYPE ), POINTER :: EM - -C----------------------------------------------------------------------- - - IF ( FIRSTIME ) THEN - FIRSTIME = .FALSE. - LOGDEV = SETUP_LOGDEV() - END IF - - EM => AQM_EMIS_GET( ETYPE ) - IF ( .NOT.ASSOCIATED( EM ) ) RETURN - -C For each time step, compute the layer fractions... - - WRITE( XMSG,'(A, I7.6)' ) - & 'Calculating emissions point source layer fractions for', JTIME - WRITE( LOGDEV,* ) ' ' - CALL M3MSG2( XMSG ) - C ... initialize emission arrays ... VDEMIS_PT = 0.0 ! array assignment VDEMIS_PT_FIRE = 0.0 ! array assignment PMEMIS_PT = 0.0 ! array assignment -C ... initialize vertical fraction arrays ... -C ... fire emissions are added to surface only by default ... - - VFRAC = 0.0 - VFRAC(:,:,1) = 1.0 - -C Retrieve fire emissions and distribute according to plume-rise algorithm - -C ... plumerise ... - - SELECT CASE ( TRIM( EM % PLUMERISE ) ) - CASE ( "sofiev" ) - N = INDEX1( EMFRP, SIZE( EM % SPECIES ), EM % SPECIES ) - IF ( N > 0 ) THEN -C read in frp - BUFFER = 0.0 - CALL AQM_EMIS_READ( ETYPE, EMFRP, BUFFER, RC=LOCALRC ) - IF ( AQM_RC_CHECK( LOCALRC, - & MSG="failure while reading frp from " // - & TRIM( ETYPE ) // " emissions", - & FILE=__FILE__, LINE=__LINE__ ) ) RETURN - CALL AQM_PLUME_SOFIEV( EM, BUFFER, VFRAC, RC=LOCALRC ) - IF ( AQM_RC_CHECK( LOCALRC, - & MSG="failed to compute plume rise", - & FILE=__FILE__, LINE=__LINE__ ) ) RETURN - ELSE - CALL M3WARN( PNAME, JDATE, JTIME, - & EMFRP // "species not found. " // - & "Adding fire emissions to surface only" ) - END IF - CASE ( "none" ) -C no plume rise - CASE DEFAULT -C plume rise is disabled by default - END SELECT - -C ... gas species ... - - IS_NOT_NVPOA = ( INDEX( MECHNAME, 'NVPOA' ) .EQ. 0 ) - - DO S = 1, N_GSPC_EMIS - N = PTEM_MAP( S ) - IF ( N .GT. 0 ) THEN - BUFFER = 0.0 - CALL AQM_EMIS_READ( ETYPE, EM % TABLE( N, 1 ), BUFFER, RC=LOCALRC ) - IF ( AQM_RC_CHECK( LOCALRC, MSG="Failure while reading " // - & TRIM( EM % TABLE( N, 1 ) ) // " from " // TRIM( ETYPE ) // " emissions", - & FILE=__FILE__, LINE=__LINE__)) RETURN -C Read Non-Carbon Organic Matter too if POC is Requested - SAVE_POC = .FALSE. - IF ( IS_NOT_NVPOA .AND. EM % TABLE( N, 1 ) .EQ. 'POC' ) THEN - CALL AQM_EMIS_READ( ETYPE, 'PNCOM', BUFFER, RC=LOCALRC ) - IF ( AQM_RC_CHECK( LOCALRC, - & MSG="Failure while reading PNCOM" // - & " from " // TRIM( ETYPE ) // " emissions", - & FILE=__FILE__, LINE=__LINE__)) RETURN - SAVE_POC = IS_NOT_NVPOA - END IF - DO L = 1, EMLYRS - K = 0 - DO R = 1, MY_NROWS - DO C = 1, MY_NCOLS - K = K + 1 - VDEMIS_PT( C,R,L,N ) = VFRAC( C,R,L ) * BUFFER( K ) - END DO - END DO - END DO - IF ( SAVE_POC ) THEN - DO L = 1, EMLYRS - K = 0 - DO R = 1, MY_NROWS - DO C = 1, MY_NCOLS - K = K + 1 - VDEMIS_PT_FIRE( C,R,L,N ) = VDEMIS_PT( C,R,L,N ) - END DO - END DO - END DO - END IF - END IF - END DO - -C ... aerosol species ... - - DO S = 1, N_SPC_PTPM - V = PTPM_MAP( S ) - BUFFER = 0.0 - CALL AQM_EMIS_READ( ETYPE, PMEM_MAP_NAME( V ), BUFFER, RC=LOCALRC ) - IF ( AQM_RC_CHECK( LOCALRC, MSG="Failure while reading " // - & TRIM( PMEM_MAP_NAME( V ) ) // " from " // TRIM( ETYPE ) // " emissions", - & FILE=__FILE__, LINE=__LINE__)) RETURN - DO L = 1, PM_EMLYRS - K = 0 - DO R = 1, MY_NROWS - DO C = 1, MY_NCOLS - K = K + 1 - PMEMIS_PT( C,R,L,S ) = VFRAC( C,R,L ) * BUFFER( K ) - END DO - END DO - END DO - END DO - - RETURN + CALL GET_PT3D_FIRE_EMIS ( JDATE, JTIME ) + CALL GET_PT3D_STKS_EMIS ( JDATE, JTIME ) END SUBROUTINE GET_PT3D_EMIS diff --git a/src/model/src/PT3D_FIRE_DEFN.F b/src/model/src/PT3D_FIRE_DEFN.F new file mode 100644 index 00000000..65fb0e6f --- /dev/null +++ b/src/model/src/PT3D_FIRE_DEFN.F @@ -0,0 +1,294 @@ + MODULE PT3D_FIRE_DEFN + + USE PT3D_DATA_MOD + + IMPLICIT NONE + + REAL, ALLOCATABLE :: VFRAC( :,:,: ) ! vertical fraction + REAL, ALLOCATABLE :: BUFFER( : ) ! emission buffer + +C Species names from input file used for point source non-PM emissions mapping + INTEGER, ALLOCATABLE :: PTEM_FIRE_MAP( : ) + INTEGER, ALLOCATABLE :: PTPM_FIRE_MAP( : ) + +C Emission layers for sources within domain + INTEGER :: EMLYRS ! no. of emis layers + INTEGER :: PM_EMLYRS ! no. of emis layers for PM + +C Emission counters + INTEGER :: N_GSPC_FIRE_EMIS ! number of gas species in diagnostic file + INTEGER :: N_SPC_FIRE_PTPM + + CHARACTER( 240 ) :: XMSG = ' ' + +C Emission type + CHARACTER( * ), PARAMETER :: ETYPE = 'gbbepx' + +C Species name for Fire Radiative Power + CHARACTER( * ), PARAMETER :: EMFRP = 'FRP' + + PRIVATE + + PUBLIC PT3D_FIRE_INIT, GET_PT3D_FIRE_EMIS + + CONTAINS + +C======================================================================= + + FUNCTION PT3D_FIRE_INIT ( N_SPC_EMIS, EMLAYS, JDATE, JTIME ) + & RESULT ( SUCCESS ) + + USE AQM_EMIS_MOD + USE GRID_CONF ! horizontal & vertical domain specifications + USE CGRID_SPCS ! CGRID mechanism species + USE PTMAP ! defines pt src species mapping to VDEMIS* arrays + USE UTILIO_DEFN + + IMPLICIT NONE + +C Includes: +C None + +C Arguments: + INTEGER, INTENT( IN ) :: N_SPC_EMIS ! total no. of model emissions species + INTEGER, INTENT( IN ) :: EMLAYS ! number of emissions layers + INTEGER, INTENT( IN ) :: JDATE ! Julian date (YYYYDDD) + INTEGER, INTENT( IN ) :: JTIME ! time (HHMMSS) + + LOGICAL SUCCESS + +C Parameters: + +C Local Variables: + CHARACTER( 16 ) :: PNAME = 'PT3D_INIT ' ! procedure name + CHARACTER( 16 ) :: VNAME ! variable name buffer + CHARACTER( 16 ), SAVE, ALLOCATABLE :: STKGNAME( : ) ! stack groups file name + + INTEGER IOS ! i/o and allocate memory status + + + INTEGER IDX + INTEGER N, NSPC, NSPC1, NSPC2, NSPC3 + INTEGER S, S_OFFSET, V + + INTEGER, ALLOCATABLE :: MAP( : ) + INTEGER, ALLOCATABLE :: SPC_PTEM_FIRE_FAC( : ) + INTEGER, ALLOCATABLE :: SPC_PTEM_FIRE_MAP( : ) + + TYPE( AQM_INTERNAL_EMIS_TYPE ), POINTER :: EM + +C----------------------------------------------------------------------- + + SUCCESS = .TRUE. + +C check if emissions are being provided + + EM => AQM_EMIS_GET( ETYPE ) + IF ( .NOT.ASSOCIATED( EM ) ) RETURN + +C set number of emissions layers depending on whether plumerise is on + + CALL AQM_EMIS_DESC( ETYPE, NLAYS=EMLYRS ) + PM_EMLYRS = EMLYRS + +C get point source emission mapping + + IF (.NOT.PTMAP_TYPE_INIT( EM, + & N_GSPC_FIRE_EMIS, N_SPC_FIRE_PTPM, + & PTEM_FIRE_MAP, PTPM_FIRE_MAP, + & SPC_PTEM_FIRE_FAC, SPC_PTEM_FIRE_MAP, + & "FIRE" ) ) THEN + XMSG = 'Could not get point source mappings' + CALL M3WARN ( PNAME, JDATE, JTIME, XMSG ) + SUCCESS = .FALSE.; RETURN + END IF + + DEALLOCATE( SPC_PTEM_FIRE_FAC, SPC_PTEM_FIRE_MAP ) + + ALLOCATE ( VFRAC( NCOLS,NROWS,EMLYRS ), STAT = IOS ) + CALL CHECKMEM( IOS, 'VFRAC', PNAME ) + VFRAC = 1.0 ! array + + ALLOCATE ( BUFFER( NCOLS * NROWS ), STAT = IOS ) + CALL CHECKMEM( IOS, 'BUFFER', PNAME ) + BUFFER = 0.0 ! array + + SUCCESS = .TRUE.; RETURN + + END FUNCTION PT3D_FIRE_INIT + +C======================================================================= + + SUBROUTINE GET_PT3D_FIRE_EMIS ( JDATE, JTIME ) + +! Revision History. +! Aug 12, 15 D. Wong: added code to handle parallel I/O implementation + +C----------------------------------------------------------------------- + +C Time step part of laypoint + + USE AQM_EMIS_MOD + USE AQM_FIRES_MOD + USE AQM_RC_MOD + USE RXNS_DATA, ONLY : MECHNAME !Get Chemical Mechanism Name + USE GRID_CONF ! horizontal & vertical domain specifications + USE CGRID_SPCS ! CGRID mechanism species + USE AERO_DATA, ONLY : PMEM_MAP_NAME + USE PTMAP ! defines pt src species mapping to VDEMIS* arrays + USE UTILIO_DEFN + + IMPLICIT NONE + +C Arguments: + INTEGER, INTENT( IN ) :: JDATE, JTIME + +C Parameters: + +C External functions: + INTEGER, EXTERNAL :: SETUP_LOGDEV + +C Local variables: + CHARACTER( 16 ) :: PNAME = 'GET_PT3D_EMIS ' ! procedure name + CHARACTER( 16 ) :: VNAME ! variable name buffer + + INTEGER IOS ! i/o and allocate memory status + INTEGER L, S, V ! counters + INTEGER C, R, K, M, N + INTEGER LOCALRC + + LOGICAL :: IS_NOT_NVPOA, SAVE_POC + + LOGICAL, SAVE :: FIRSTIME = .TRUE. + INTEGER, SAVE :: LOGDEV + + TYPE( AQM_INTERNAL_EMIS_TYPE ), POINTER :: EM + +C----------------------------------------------------------------------- + + IF ( FIRSTIME ) THEN + FIRSTIME = .FALSE. + LOGDEV = SETUP_LOGDEV() + END IF + + EM => AQM_EMIS_GET( ETYPE ) + IF ( .NOT.ASSOCIATED( EM ) ) RETURN + +C For each time step, compute the layer fractions... + + WRITE( XMSG,'(A, I7.6)' ) + & 'Calculating emissions point source layer fractions for', JTIME + WRITE( LOGDEV,* ) ' ' + CALL M3MSG2( XMSG ) + +C ... initialize vertical fraction arrays ... +C ... fire emissions are added to surface only by default ... + + VFRAC = 0.0 + VFRAC(:,:,1) = 1.0 + +C Retrieve fire emissions and distribute according to plume-rise algorithm + +C ... plumerise ... + + SELECT CASE ( TRIM( EM % PLUMERISE ) ) + CASE ( "sofiev" ) + N = INDEX1( EMFRP, SIZE( EM % SPECIES ), EM % SPECIES ) + IF ( N > 0 ) THEN +C read in frp + BUFFER = 0.0 + CALL AQM_EMIS_READ( ETYPE, EMFRP, BUFFER, RC=LOCALRC ) + IF ( AQM_RC_CHECK( LOCALRC, + & MSG="failure while reading frp from " // + & TRIM( ETYPE ) // " emissions", + & FILE=__FILE__, LINE=__LINE__ ) ) RETURN + CALL AQM_PLUME_SOFIEV( EM, BUFFER, VFRAC, RC=LOCALRC ) + IF ( AQM_RC_CHECK( LOCALRC, + & MSG="failed to compute plume rise", + & FILE=__FILE__, LINE=__LINE__ ) ) RETURN + ELSE + CALL M3WARN( PNAME, JDATE, JTIME, + & EMFRP // "species not found. " // + & "Adding fire emissions to surface only" ) + END IF + CASE ( "none" ) +C no plume rise + CASE DEFAULT +C plume rise is disabled by default + END SELECT + +C ... gas species ... + + IS_NOT_NVPOA = ( INDEX( MECHNAME, 'NVPOA' ) .EQ. 0 ) + + DO S = 1, N_GSPC_FIRE_EMIS + M = PTEM_FIRE_MAP( S ) + IF ( M .GT. 0 ) THEN + N = PTEM_MAP( S ) +c IF ( N .GT. 0 ) THEN + BUFFER = 0.0 + CALL AQM_EMIS_READ( ETYPE, EM % TABLE( M, 1 ), + & BUFFER, RC=LOCALRC ) + IF ( AQM_RC_CHECK( LOCALRC, MSG="Failure while reading " + & // TRIM( EM % TABLE( M, 1 ) ) // " from " + & // TRIM( ETYPE ) // " emissions", + & FILE=__FILE__, LINE=__LINE__)) RETURN +C Read Non-Carbon Organic Matter too if POC is Requested + SAVE_POC = .FALSE. + IF ( IS_NOT_NVPOA .AND. EM % TABLE( M, 1 ) .EQ. 'POC' ) THEN + CALL AQM_EMIS_READ( ETYPE, 'PNCOM', BUFFER, RC=LOCALRC ) + IF ( AQM_RC_CHECK( LOCALRC, + & MSG="Failure while reading PNCOM" // + & " from " // TRIM( ETYPE ) // " emissions", + & FILE=__FILE__, LINE=__LINE__)) RETURN + SAVE_POC = IS_NOT_NVPOA + END IF + DO L = 1, EMLYRS + K = 0 + DO R = 1, MY_NROWS + DO C = 1, MY_NCOLS + K = K + 1 + VDEMIS_PT( C,R,L,N ) = VFRAC( C,R,L ) * BUFFER( K ) + END DO + END DO + END DO + IF ( SAVE_POC ) THEN + DO L = 1, EMLYRS + K = 0 + DO R = 1, MY_NROWS + DO C = 1, MY_NCOLS + K = K + 1 + VDEMIS_PT_FIRE( C,R,L,N ) = VDEMIS_PT( C,R,L,N ) + END DO + END DO + END DO + END IF + END IF + END DO + +C ... aerosol species ... + + DO S = 1, N_SPC_FIRE_PTPM + N = PTPM_FIRE_MAP( S ) + V = PTPM_MAP( N ) + BUFFER = 0.0 + CALL AQM_EMIS_READ( ETYPE, PMEM_MAP_NAME( V ), BUFFER, RC=LOCALRC ) + IF ( AQM_RC_CHECK( LOCALRC, MSG="Failure while reading " // + & TRIM( PMEM_MAP_NAME( V ) ) // " from " // TRIM( ETYPE ) // " emissions", + & FILE=__FILE__, LINE=__LINE__)) RETURN + DO L = 1, PM_EMLYRS + K = 0 + DO R = 1, MY_NROWS + DO C = 1, MY_NCOLS + K = K + 1 + PMEMIS_PT( C,R,L,N ) = VFRAC( C,R,L ) * BUFFER( K ) + END DO + END DO + END DO + END DO + + RETURN + + END SUBROUTINE GET_PT3D_FIRE_EMIS + + END MODULE PT3D_FIRE_DEFN diff --git a/src/model/src/PT3D_STKS_DEFN.F b/src/model/src/PT3D_STKS_DEFN.F new file mode 100644 index 00000000..3a0cf2a5 --- /dev/null +++ b/src/model/src/PT3D_STKS_DEFN.F @@ -0,0 +1,501 @@ + MODULE PT3D_STKS_DEFN + + USE GRID_CONF ! horizontal & vertical domain specifications + USE PTMAP + USE AQM_EMIS_MOD, ONLY : AQM_EMIS_DESC, AQM_EMIS_GET, + & AQM_EMIS_READ, AQM_INTERNAL_EMIS_TYPE + + IMPLICIT NONE + + INTEGER :: N_GSPC_STKS_EMIS + INTEGER :: N_SPC_PT3DEM + INTEGER :: N_SPC_STKS_PTPM + + INTEGER, ALLOCATABLE :: SPC_PTEM_STKS_FAC( : ) + INTEGER, ALLOCATABLE :: SPC_PTEM_STKS_MAP( : ) + INTEGER, ALLOCATABLE :: PTEM_STKS_MAP( : ) + INTEGER, ALLOCATABLE :: PTPM_STKS_MAP( : ) + + REAL :: CNVTP ! intermediate combined conv. factor + +C Emission type + CHARACTER( * ), PARAMETER :: ETYPE = 'point-source' + + PRIVATE + + PUBLIC :: N_SPC_STKS_PTPM, PTPM_STKS_MAP + PUBLIC :: GET_PT3D_STKS_EMIS, PT3D_STKS_INIT + + CONTAINS + + SUBROUTINE GET_PT3D_STKS_EMIS( JDATE, JTIME ) + + USE CGRID_SPCS ! CGRID mechanism species + USE UTILIO_DEFN + USE AQM_RC_MOD + USE AERO_DATA, ONLY : PMEM_MAP_NAME + USE ASX_DATA_MOD, ONLY: MET_DATA, GRID_DATA, CONVPA + USE PT3D_DATA_MOD + + IMPLICIT NONE + +C Includes: + INCLUDE SUBST_FILES_ID ! file name parameters + INCLUDE SUBST_CONST + +C Arguments: + INTEGER, INTENT( IN ) :: JDATE ! Julian date (YYYYDDD) + INTEGER, INTENT( IN ) :: JTIME ! time (HHMMSS) + +C External Functions: + INTEGER, EXTERNAL :: SETUP_LOGDEV + +C Local Variables: + INTEGER :: LOGDEV + INTEGER :: C, R, L, M, N, S, V, K, I ! loop induction variables + INTEGER :: LBOT ! layer containing plume bottom + INTEGER :: LTOP ! layer containing plume top + INTEGER :: LPBL ! first L: ZF(L) above mixing layer - ONLY for REPORT + INTEGER :: LSTK ! first L: ZF(L) > STKHT + INTEGER :: IJ, IS + INTEGER :: NSRC, NP + INTEGER :: IOS, LOCALRC + INTEGER, POINTER :: IP( : ) + INTEGER, POINTER :: JP( : ) + INTEGER, POINTER :: IJMAP( : ) + REAL :: MV ! mininum LFRAC + REAL :: PSFC ! surface pressure [Pa] + REAL :: TSTK ! temperature at top of stack [K] + REAL :: TSUM ! tmp layer frac sum for renormalizing + REAL :: WSTK ! wind speed at top of stack [m/s] + REAL :: ZBOT ! plume bottom elevation [m] + REAL :: ZTOP ! plume top elevation [m] + REAL :: ZDIFF ! ZTOP - ZBOT + REAL :: ZF0, ZF1 + REAL :: DDZ ! 1 / ZDIFF + REAL :: ZPLM ! plume centerline height above stack [m] + REAL :: USTMP ! temp storage for ustar [m/s] + REAL :: HFLX ! converted heat flux + REAL, ALLOCATABLE :: BUFFER( : ) + REAL, ALLOCATABLE :: DTHDZ( : ) + REAL, ALLOCATABLE :: PRESF( : ) + REAL, ALLOCATABLE :: WSPD( : ) + REAL, ALLOCATABLE :: ZZF( : ) + REAL, ALLOCATABLE :: ZSTK( : ) + REAL, ALLOCATABLE :: DDZF( : ) + REAL, ALLOCATABLE :: TFRAC( : ) + REAL, ALLOCATABLE :: VFRAC( :, : ) + +C Parameters: + REAL, PARAMETER :: CMLMR = 1.0E+06 ! ppmV/Molar Mixing Ratio + REAL, PARAMETER :: CNVPA2MB = 1.0E-2 ! convert Pa to mb + REAL, PARAMETER :: USTARMIN = 0.1 ! Min valid value for USTAR + CHARACTER( 8 ) :: CINT ! integer to character buffer for Cwarning messages + CHARACTER( 16 ) :: PNAME = 'GET_PT3D_EMIS' + CHARACTER( 120 ) :: XMSG = ' ' + CHARACTER( 10 ), PARAMETER :: BLANK10 = ' ' + + LOGICAL, SAVE :: FIRSTIME = .TRUE. + INTEGER, SAVE :: EMLYRS + + TYPE( AQM_INTERNAL_EMIS_TYPE ), POINTER :: EM + +C----------------------------------------------------------------------- + + INTERFACE + SUBROUTINE PREPLM( FIREFLG, EMLAYS, HMIX, HTS, PSFC, TS, DDZF, QV, + & TA, UW, VW, ZH, ZF, PRES, LSTK, LPBL, TSTK, + & WSTK, DTHDZ, WSPD ) + LOGICAL, INTENT( IN ) :: FIREFLG + INTEGER, INTENT( IN ) :: EMLAYS + REAL, INTENT( IN ) :: HMIX + REAL, INTENT( IN ) :: HTS + REAL, INTENT( IN ) :: PSFC + REAL, INTENT( IN ) :: TS + REAL, INTENT( IN ) :: DDZF( : ) + REAL, INTENT( IN ) :: QV ( : ) + REAL, INTENT( IN ) :: TA ( : ) + REAL, INTENT( IN ) :: UW ( : ) + REAL, INTENT( IN ) :: VW ( : ) + REAL, INTENT( IN ) :: ZH ( : ) + REAL, INTENT( IN ) :: ZF ( : ) + REAL, INTENT( IN ) :: PRES( 0: ) + INTEGER, INTENT( OUT ) :: LSTK + INTEGER, INTENT( OUT ) :: LPBL + REAL, INTENT( OUT ) :: TSTK + REAL, INTENT( OUT ) :: WSTK + REAL, INTENT( OUT ) :: DTHDZ( : ) + REAL, INTENT( OUT ) :: WSPD ( : ) + END SUBROUTINE PREPLM + + SUBROUTINE PLMRIS( EMLAYS, LSTK, HFX, HMIX, + & STKDM, STKHT, STKTK, STKVE, + & TSTK, USTAR, DTHDZ, TA, WSPD, + & ZF, ZH, ZSTK, WSTK, ZPLM ) + INTEGER, INTENT( IN ) :: EMLAYS + INTEGER, INTENT( IN ) :: LSTK + REAL, INTENT( IN ) :: HFX + REAL, INTENT( IN ) :: HMIX + REAL, INTENT( IN ) :: STKDM + REAL, INTENT( IN ) :: STKHT + REAL, INTENT( IN ) :: STKTK + REAL, INTENT( IN ) :: STKVE + REAL, INTENT( IN ) :: TSTK + REAL, INTENT( IN ) :: USTAR + REAL, INTENT( IN ) :: DTHDZ( : ) + REAL, INTENT( IN ) :: TA ( : ) + REAL, INTENT( IN ) :: WSPD ( : ) + REAL, INTENT( IN ) :: ZF ( 0: ) + REAL, INTENT( IN ) :: ZH ( : ) + REAL, INTENT( IN ) :: ZSTK ( : ) + REAL, INTENT( INOUT ) :: WSTK + REAL, INTENT( OUT ) :: ZPLM + END SUBROUTINE PLMRIS + + SUBROUTINE PLSPRD( DTHDZ, ZF, KZ, CEFSTK, PLTOP, PLBOT ) + REAL, INTENT ( IN ) :: DTHDZ( : ) + REAL, INTENT ( IN ) :: ZF( 0: ) + INTEGER, INTENT ( IN ) :: KZ + REAL, INTENT ( IN ) :: CEFSTK + REAL, INTENT( OUT ) :: PLTOP + REAL, INTENT( OUT ) :: PLBOT + END SUBROUTINE PLSPRD + + END INTERFACE + +C----------------------------------------------------------------------- + + EM => AQM_EMIS_GET( ETYPE ) + IF ( .NOT.ASSOCIATED( EM ) ) RETURN + + IF (EM % COUNT == 0) RETURN + + IF ( FIRSTIME ) THEN + + LOGDEV = SETUP_LOGDEV() + +C set number of emissions layers depending on whether plumerise is on + + CALL AQM_EMIS_DESC( ETYPE, NLAYS=EMLYRS ) + + FIRSTIME = .FALSE. + END IF + +C Allocate Buffer space for Reading Emissions + NSRC = SIZE( EM % IJMAP ) + + ALLOCATE ( BUFFER( SIZE( EM % LAT ) ), STAT = IOS ) + CALL CHECKMEM( IOS, 'BUFFER', PNAME ) + + ALLOCATE ( VFRAC( NSRC, EMLYRS ), STAT = IOS ) + CALL CHECKMEM( IOS, 'TFRAC', PNAME ) + + ALLOCATE ( TFRAC( EMLYRS ), STAT = IOS ) + CALL CHECKMEM( IOS, 'TFRAC', PNAME ) + + IF ( EMLYRS .GT. 1 ) THEN + ALLOCATE ( DTHDZ( EMLYRS ), STAT = IOS ) + CALL CHECKMEM( IOS, 'DTHDZ', PNAME ) + + ALLOCATE ( PRESF( 0:EMLYRS ), STAT = IOS ) + CALL CHECKMEM( IOS, 'PRESF', PNAME ) + + ALLOCATE ( WSPD( EMLYRS ), STAT = IOS ) + CALL CHECKMEM( IOS, 'WSPD', PNAME ) + + ALLOCATE ( ZZF( 0:EMLYRS ), STAT = IOS ) + CALL CHECKMEM( IOS, 'ZZF', PNAME ) + + ALLOCATE ( ZSTK( EMLYRS ), STAT = IOS ) + CALL CHECKMEM( IOS, 'ZSTK', PNAME ) + + ALLOCATE ( DDZF( EMLYRS ), STAT = IOS ) + CALL CHECKMEM( IOS, 'DDZF', PNAME ) + END IF + + ! -- loop over mapped emissions only + + DO IS = 1, NSRC + IJ = EM % IJMAP( IS ) + C = EM % IP( IJ ) + R = EM % JP( IJ ) + + IF ( EMLYRS .GT. 1 ) THEN +C Loop through sources and compute plume rise + + ZZF( 0 ) = 0.0 + ZZF( 1:EMLYRS ) = Met_Data % ZF( C,R,1:EMLYRS ) + +C Set surface pressure (convert to mb from Pa) + PSFC = CNVPA2MB * Met_Data % PRSFC( C,R ) + + PRESF( 0:EMLYRS ) = CNVPA2MB * Met_Data % PRESF( C,R,: ) + +C Compute vertical gradient quantities + ZF0 = Met_Data % ZF( C,R,1 ) + ZSTK( 1 ) = ZF0 - EM % STKHT( IJ ) + DDZF( 1 ) = 1.0 / ZF0 + + DO L = 2, EMLYRS + ZF1 = Met_Data % ZF( C,R,L ) + ZSTK( L ) = ZF1 - EM % STKHT( IJ ) + DDZF( L ) = 1.0 / ( ZF1 - ZF0 ) + ZF0 = ZF1 + END DO + +C Compute derived met vars needed before layer assignments + CALL PREPLM( .FALSE., EMLYRS, + & Met_Data % PBL ( C,R ), EM % STKHT( IJ ), PSFC, + & Met_Data % TEMP2 ( C,R ), DDZF, + & Met_Data % QV ( C,R,: ), Met_Data % TA ( C,R,: ), + & Met_Data % UWINDA( C,R,: ), Met_Data % VWINDA ( C,R,: ), + & Met_Data % ZH ( C,R,: ), Met_Data % ZF ( C,R,: ), + & PRESF, LSTK, LPBL, TSTK, WSTK, + & DTHDZ, WSPD ) + +C Trap USTAR at a minimum realistic value + USTMP = MAX( Met_Data % USTAR( C,R ), USTARMIN ) + +C Convert heat flux (watts/m2 to m K /s ) + HFLX = Met_Data % HFX( C,R ) / ( CPD * Met_Data % DENS1( C,R ) ) + + CALL PLMRIS( EMLYRS, LSTK, HFLX, Met_Data % PBL( C,R ), + & EM % STKDM( IJ ), EM % STKHT( IJ ), + & EM % STKTK( IJ ), EM % STKVE( IJ ), + & TSTK, USTMP, + & DTHDZ, Met_Data % TA( C,R,: ), + & WSPD, ZZF, + & Met_Data % ZH( C,R,: ), ZSTK, + & WSTK, ZPLM ) + +C Determine the bottom and top heights of the plume. +C Default Turner approach. Plume thickness = amount of plume rise +C Plume rise DH = ZPLM minus the stack height STKHT + ZTOP = EM % STKHT( IJ ) + & + 1.5 * ( ZPLM - EM % STKHT( IJ ) ) + ZBOT = EM % STKHT( IJ ) + & + 0.5 * ( ZPLM - EM % STKHT( IJ ) ) + +C Set up for computing plume fractions, assuming uniform distribution in pressure +C (~mass concentration -- minor hydrostatic assumption) from bottom to top. + + IF ( ZTOP .LT. EM % STKHT( IJ ) ) THEN + WRITE( CINT,'( I8 )' ) S + WRITE( XMSG,94010 ) 'ERROR: Top of plume is less than ' + & // 'top of stack for source:' // CINT + CALL M3MESG( XMSG ) + WRITE( LOGDEV,* ) ' Zbot: ', ZBOT, ' Ztop: ', ZTOP + WRITE( LOGDEV,* ) ' Stack Top: ', EM % STKHT( IJ ), + & ' Plume Top: ', ZPLM + CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT2 ) + END IF + +C Allocate plume to layers (compute layer plume fractions) + +C Compute LBOT, LTOP such that +C ZZF( LBOT-1 ) <= ZBOT < ZZF( LBOT ) and +C ZZF( LTOP-1 ) <= ZTOP < ZZF( LTOP ) + + DO L = 1, EMLYRS - 1 +c IF ( ZBOT .LE. Met_Data % ZF( C,R,L ) ) THEN + IF ( ZBOT .LE. ZZF( L ) ) THEN + LBOT = L + GO TO 122 + ELSE + TFRAC( L ) = 0.0 ! fractions below plume + END IF + END DO + LBOT = EMLYRS ! fallback + +122 CONTINUE ! loop exit: bottom found at LBOT + + IF ( ZTOP .LE. ZZF( LBOT ) ) THEN ! plume in this layer + + TFRAC( LBOT ) = 1.0 + LTOP = LBOT + + DO L = LBOT + 1, EMLYRS ! fractions above plume + TFRAC( L ) = 0.0 + END DO + + ELSE IF ( LBOT .EQ. EMLYRS ) THEN ! plume above top layer + + TFRAC( LBOT ) = 1.0 + + DO L = 1, EMLYRS - 1 ! fractions below plume + TFRAC( L ) = 0.0 + END DO + + ELSE ! plume crosses layers + + DO L = LBOT + 1, EMLYRS +c IF ( ZTOP .LE. Met_Data % ZF( C,R,L ) ) THEN + IF ( ZTOP .LE. ZZF( L ) ) THEN + LTOP = L + GO TO 126 + END IF + END DO + LTOP = EMLYRS ! fallback + +126 CONTINUE + + ZDIFF = ZTOP - ZBOT + IF ( ZDIFF .GT. 0.0 ) THEN + + DDZ = 1.0 / ZDIFF +c TFRAC( LBOT ) = DDZ * ( Met_Data % ZF( C,R,LBOT ) - ZBOT ) +c TFRAC( LTOP ) = DDZ * ( ZTOP - Met_Data % ZF( C,R,LTOP-1 ) ) + TFRAC( LBOT ) = DDZ * ( ZZF( LBOT ) - ZBOT ) + TFRAC( LTOP ) = DDZ * ( ZTOP - ZZF( LTOP-1 ) ) + + ELSE ! ZDIFF .le. 0 + + WRITE( CINT,'( I8 )' ) S + WRITE( XMSG,94020 ) + & 'Infinitely small plume created for source:,' + & // CINT // CRLF() // BLANK10 + & // 'All emissions put in first layer.' + CALL M3WARN( PNAME, JDATE, JTIME, XMSG ) + LBOT = 1; LTOP = 1 + TFRAC( LBOT ) = 1.0 + + END IF + + DO L = LBOT + 1, LTOP - 1 ! layers in plume + TFRAC( L ) = DDZ * ( Met_Data % ZF( C,R,L ) - Met_Data % ZF( C,R,L-1 ) ) + END DO + + DO L = LTOP + 1, EMLYRS ! fractions above plume + TFRAC( L ) = 0.0 + END DO + + END IF + +C If layer fractions are negative, put in the first layer + MV = MINVAL( TFRAC( 1:EMLYRS ) ) + IF ( MV .LT. 0.0 ) THEN + + WRITE( CINT,'( I8 )' ) S + WRITE( XMSG,94010 ) + & 'WARNING: One or more negative plume ' + & // 'fractions found for source:' // CINT + & // CRLF() // BLANK10 // 'Plume reset to ' + & // 'put all emissions in surface layer.' + CALL M3MESG( XMSG ) + + TFRAC( 1 ) = 1.0 + TFRAC( 2:EMLYRS ) = 0.0 + + END IF + + ELSE + + TFRAC( 1:EMLYRS ) = 1.0 + + END IF + + VFRAC( IS, : ) = TFRAC + + END DO + + DEALLOCATE ( TFRAC ) + + ! -- (1) non-PM emissions + + DO S = 1, N_GSPC_STKS_EMIS + M = PTEM_STKS_MAP( S ) + IF ( M .LT. 1 ) CYCLE + N = PTEM_MAP( S ) + + BUFFER = 0.0 + CALL AQM_EMIS_READ( ETYPE, EM % TABLE( M, 1 ), BUFFER, RC=LOCALRC) + IF ( AQM_RC_CHECK( LOCALRC, + & MSG="Failure while reading " // + & TRIM( EM % TABLE( M, 1 ) ) // " from " // + & TRIM( ETYPE ) // " emissions", + & FILE=__FILE__, LINE=__LINE__ ) ) RETURN + IF ( EM % TABLE( M, 1 ) .EQ. 'POC' ) THEN + CALL AQM_EMIS_READ( ETYPE, 'PNCOM', BUFFER, RC=LOCALRC) + IF ( AQM_RC_CHECK( LOCALRC, + & MSG="Failure while reading PNCOM emissions from" // + & TRIM( ETYPE ) // " emissions", + & FILE=__FILE__, LINE=__LINE__ ) ) RETURN + END IF + + ! -- add emissions + DO L = 1, EMLYRS + DO IS = 1, NSRC + IJ = EM % IJMAP( IS ) + C = EM % IP( IJ ) + R = EM % JP( IJ ) + VDEMIS_PT( C,R,L,N ) = VDEMIS_PT( C,R,L,N ) + & + VFRAC( IS, L ) * BUFFER( IJ ) + END DO + END DO + END DO + + ! -- (2) PM emissions + + DO S = 1, N_SPC_STKS_PTPM + N = PTPM_STKS_MAP( S ) + V = PTPM_MAP( N ) + BUFFER = 0.0 + CALL AQM_EMIS_READ( ETYPE, PMEM_MAP_NAME( V ), BUFFER, RC=LOCALRC ) + IF ( AQM_RC_CHECK( LOCALRC, MSG="Failure while reading " // + & TRIM( PMEM_MAP_NAME( V ) ) // " from " // TRIM( ETYPE ) // " emissions", + & FILE=__FILE__, LINE=__LINE__)) RETURN + DO L = 1, EMLYRS + DO IS = 1, NSRC + IJ = EM % IJMAP( IS ) + C = EM % IP( IJ ) + R = EM % JP( IJ ) + PMEMIS_PT( C,R,L,N ) = PMEMIS_PT( C,R,L,N ) + & + VFRAC( IS, L ) * BUFFER( IJ ) + END DO + END DO + END DO + + ! -- free up memory + + DEALLOCATE ( BUFFER, VFRAC ) + + IF ( EMLYRS .GT. 1 ) DEALLOCATE ( DDZF, DTHDZ, PRESF, WSPD, ZSTK, ZZF ) + +C------------------ FORMAT STATEMENTS ------------------------------ + +94010 FORMAT( 12( A, :, I8, :, 1X ) ) +94020 FORMAT( 10( A, :, I7, :, 1X ) ) + + + END SUBROUTINE GET_PT3D_STKS_EMIS + + + FUNCTION PT3D_STKS_INIT( ) RESULT ( SUCCESS ) + + USE AQM_RC_MOD, ONLY: AQM_RC_TEST + + TYPE( AQM_INTERNAL_EMIS_TYPE ), POINTER :: EM + + LOGICAL :: SUCCESS + + SUCCESS = .TRUE. + + EM => AQM_EMIS_GET( ETYPE ) + IF ( .NOT.ASSOCIATED( EM ) ) RETURN + + IF (EM % COUNT == 0) RETURN + + SUCCESS = PTMAP_TYPE_INIT( EM, + & N_GSPC_STKS_EMIS, N_SPC_STKS_PTPM, + & PTEM_STKS_MAP, PTPM_STKS_MAP, + & SPC_PTEM_STKS_FAC, SPC_PTEM_STKS_MAP, + & "STKS" ) + + IF ( AQM_RC_TEST( .NOT. SUCCESS, + & MSG="Failure initializing mapping for" // + & TRIM( ETYPE ) // " emissions", + & FILE=__FILE__, LINE=__LINE__ ) ) RETURN + + END FUNCTION PT3D_STKS_INIT + + END MODULE PT3D_STKS_DEFN diff --git a/src/model/src/PTMAP.F b/src/model/src/PTMAP.F index c7e9a7ec..8d651a60 100644 --- a/src/model/src/PTMAP.F +++ b/src/model/src/PTMAP.F @@ -16,53 +16,89 @@ MODULE PTMAP INTEGER, SAVE :: N_SPC_PTEM = 0 ! merged no. of unique species INTEGER, SAVE :: N_SPC_PTPM = 0 ! merged no. of unique species for PM -C Emission type - CHARACTER( * ), PARAMETER :: ETYPE = 'gbbepx' - PRIVATE PUBLIC :: N_GSPC_EMIS, N_SPC_PTEM, N_SPC_PTPM PUBLIC :: SPC_PTEM_FAC, SPC_PTEM_MAP PUBLIC :: PTEM_MAP, PTPM_MAP + PUBLIC :: PTMAP_INIT + PUBLIC :: PTMAP_TYPE_INIT CONTAINS FUNCTION PTMAP_INIT ( ) RESULT ( SUCCESS ) - USE CGRID_SPCS ! CGRID mechanism species - USE AERO_DATA, ONLY : N_EMIS_PM, PMEM_MAP_NAME ! defines aerosol species - USE AQM_EMIS_MOD, ONLY : AQM_EMIS_GET, AQM_INTERNAL_EMIS_TYPE + USE CGRID_SPCS ! CGRID mechanism species + USE AERO_DATA, ONLY : N_EMIS_PM, PMEM_MAP_NAME ! defines aerosol species + USE AQM_EMIS_MOD, ONLY : AQM_EMIS_GET, AQM_EMIS_ISPRESENT, + & AQM_INTERNAL_EMIS_TYPE USE UTILIO_DEFN IMPLICIT NONE LOGICAL SUCCESS - INTEGER, EXTERNAL :: SETUP_LOGDEV + INTEGER I, IDX, IOS, IPT, NPT + INTEGER N, N_GAS_EMIS, NSPC, NSPC1, NSPC2, NSPC3 + INTEGER S, S_OFFSET, V - CHARACTER( 16 ) :: PNAME = 'PTMAP_INIT' ! procedure name + INTEGER, ALLOCATABLE :: MAP( : ) + CHARACTER( 16 ), ALLOCATABLE :: EMSPC( : ) - INTEGER IOS, LOGDEV + TYPE PT_EMIS_TYPE + TYPE( AQM_INTERNAL_EMIS_TYPE ), POINTER :: EM + END TYPE PT_EMIS_TYPE - INTEGER IDX - INTEGER N, N_GAS_EMIS, NSPC, NSPC1, NSPC2, NSPC3 - INTEGER S, S_OFFSET, V + TYPE( PT_EMIS_TYPE ) :: PT( 2 ) - INTEGER, ALLOCATABLE :: MAP( : ) + CHARACTER( * ), PARAMETER :: PNAME = 'PTMAP_INIT' ! procedure name - TYPE( AQM_INTERNAL_EMIS_TYPE ), POINTER :: EM + CHARACTER( * ), PARAMETER :: ETYPES( 2 ) = + & (/ 'gbbepx ', 'point-source' /) C----------------------------------------------------------------------- - LOGDEV = SETUP_LOGDEV() - SUCCESS = .TRUE. C check if emissions are being provided - EM => AQM_EMIS_GET( ETYPE ) - IF ( .NOT.ASSOCIATED( EM ) ) RETURN + NPT = SIZE( ETYPES ) + + I = 0 + NSPC = 0 + DO IPT = 1, NPT + IF ( AQM_EMIS_ISPRESENT( ETYPES( IPT ) ) ) THEN + I = I + 1 + PT( I ) % EM => AQM_EMIS_GET( ETYPES( IPT ) ) + NSPC = NSPC + SIZE( PT( I ) % EM % TABLE, DIM=1 ) + END IF + END DO + + IF ( I == 0 ) RETURN + + NPT = I + + ALLOCATE ( EMSPC( NSPC ), STAT = IOS ) + CALL CHECKMEM( IOS, 'EMSPC', PNAME ) + MAP = 0 + + EMSPC = "" + + ! -- merge emission tables and remove duplicates + + I = 0 + DO IPT = 1, NPT + DO N = 1, SIZE( PT( IPT ) % EM % TABLE, DIM=1 ) + IDX = INDEX1( TRIM( PT( IPT ) % EM % TABLE( N, 1 ) ), NSPC, EMSPC ) + IF ( IDX .LE. 0 ) THEN + I = I + 1 + EMSPC( I ) = PT( IPT ) % EM % TABLE( N, 1 ) + END IF + END DO + END DO + + NSPC = I C compute emission offsets @@ -72,8 +108,6 @@ FUNCTION PTMAP_INIT ( ) RESULT ( SUCCESS ) C create auxiliary arrays mapping fire emission species to CMAQ gas and aerosol species - NSPC = SIZE( EM % TABLE, DIM=1 ) - ALLOCATE( MAP( NSPC ), STAT = IOS ) CALL CHECKMEM( IOS, 'MAP', PNAME ) MAP = 0 @@ -98,11 +132,14 @@ FUNCTION PTMAP_INIT ( ) RESULT ( SUCCESS ) CALL CHECKMEM( IOS, 'PTEM_MAP', PNAME ) PTEM_MAP = -1 + N = 0 + S_OFFSET = 0 DO S = 1, N_GC_EMIS - IDX = INDEX1( GC_EMIS( S ), NSPC, EM % TABLE( 1, 1 ) ) + IDX = INDEX1( GC_EMIS( S ), NSPC, EMSPC ) IF ( IDX .GT. 0 ) THEN - PTEM_MAP ( S ) = IDX + N = N + 1 + PTEM_MAP ( S ) = N SPC_PTEM_MAP( S ) = S SPC_PTEM_FAC( S ) = GC_EMIS_FAC( S ) END IF @@ -110,10 +147,11 @@ FUNCTION PTMAP_INIT ( ) RESULT ( SUCCESS ) S_OFFSET = N_GC_EMIS DO S = 1, N_NR_EMIS - IDX = INDEX1( NR_EMIS( S ), NSPC, EM % TABLE( 1, 1 ) ) + IDX = INDEX1( NR_EMIS( S ), NSPC, EMSPC ) IF ( IDX .GT. 0 ) THEN + N = N + 1 V = S + S_OFFSET - PTEM_MAP ( V ) = IDX + PTEM_MAP ( V ) = N SPC_PTEM_MAP( V ) = S + NSPC2 SPC_PTEM_FAC( V ) = NR_EMIS_FAC( S ) END IF @@ -121,37 +159,183 @@ FUNCTION PTMAP_INIT ( ) RESULT ( SUCCESS ) S_OFFSET = S_OFFSET + N_NR_EMIS DO S = 1, N_TR_EMIS - IDX = INDEX1( TR_EMIS( S ), NSPC, EM % TABLE( 1, 1 ) ) + IDX = INDEX1( TR_EMIS( S ), NSPC, EMSPC ) IF ( IDX .GT. 0 ) THEN + N = N + 1 V = S + S_OFFSET - PTEM_MAP ( V ) = IDX + PTEM_MAP ( V ) = N SPC_PTEM_MAP( V ) = S + NSPC3 SPC_PTEM_FAC( V ) = TR_EMIS_FAC( S ) END IF END DO + N_SPC_PTEM = N N_GSPC_EMIS = N_GAS_EMIS C ... aerosol species ... N = 0 DO S = 1, NSPC - IDX = INDEX1( EM % TABLE( S, 1 ), N_EMIS_PM, PMEM_MAP_NAME ) + IDX = INDEX1( EMSPC( S ), N_EMIS_PM, PMEM_MAP_NAME ) IF ( IDX .GT. 0 ) THEN N = N + 1 MAP( N ) = IDX END IF END DO - N_SPC_PTEM = NSPC N_SPC_PTPM = N ALLOCATE( PTPM_MAP( N_SPC_PTPM ), STAT = IOS ) CALL CHECKMEM( IOS, 'PTPM_MAP', PNAME ) PTPM_MAP = MAP( 1:N_SPC_PTPM ) - DEALLOCATE( MAP ) + DEALLOCATE( EMSPC, MAP ) END FUNCTION PTMAP_INIT + + FUNCTION PTMAP_TYPE_INIT( EM, + & N_GSPC_TYPE_EMIS, N_SPC_TYPE_PTPM, + & PTEM_TYPE_MAP, PTPM_TYPE_MAP, + & SPC_PTEM_TYPE_FAC, SPC_PTEM_TYPE_MAP, + & ELABEL ) + & RESULT ( SUCCESS ) + + + USE CGRID_SPCS ! CGRID mechanism species + USE AERO_DATA, ONLY : N_EMIS_PM, PMEM_MAP_NAME ! defines aerosol species + USE AQM_EMIS_MOD, ONLY : AQM_INTERNAL_EMIS_TYPE + USE UTILIO_DEFN + + IMPLICIT NONE + + ! -- arguments + + TYPE( AQM_INTERNAL_EMIS_TYPE ), POINTER :: EM + INTEGER, INTENT( OUT ) :: N_GSPC_TYPE_EMIS + INTEGER, INTENT( OUT ) :: N_SPC_TYPE_PTPM + INTEGER, ALLOCATABLE, INTENT( OUT ) :: PTEM_TYPE_MAP( : ) + INTEGER, ALLOCATABLE, INTENT( OUT ) :: PTPM_TYPE_MAP( : ) + INTEGER, ALLOCATABLE, INTENT( OUT ) :: SPC_PTEM_TYPE_FAC( : ) + INTEGER, ALLOCATABLE, INTENT( OUT ) :: SPC_PTEM_TYPE_MAP( : ) + CHARACTER( * ), INTENT( IN ) :: ELABEL + + ! -- local variables + + LOGICAL SUCCESS + + INTEGER IOS + + INTEGER IDX + INTEGER N, N_GAS_EMIS, NSPC, NSPC1, NSPC2, NSPC3 + INTEGER S, S_OFFSET, V + + INTEGER, ALLOCATABLE :: MAP( : ) + + CHARACTER( * ), PARAMETER :: PNAME = 'PTEM_TYPE_MAP' ! procedure name + +C----------------------------------------------------------------------- + + SUCCESS = .TRUE. + + N_GSPC_TYPE_EMIS = 0 + N_SPC_TYPE_PTPM = 0 + + IF ( .NOT.ASSOCIATED( EM ) ) RETURN + +C compute emission offsets + + NSPC1 = N_GC_EMIS + NSPC2 = NSPC1 + N_AE_EMIS + NSPC3 = NSPC2 + N_NR_EMIS + +C create auxiliary arrays mapping fire emission species to CMAQ gas and aerosol species + + NSPC = SIZE( EM % TABLE, DIM=1 ) + + ALLOCATE( MAP( NSPC ), STAT = IOS ) + CALL CHECKMEM( IOS, 'MAP: ' // TRIM( ELABEL ), PNAME ) + MAP = 0 + +C ... gas species ... + + NSPC1 = N_GC_EMIS + NSPC2 = NSPC1 + N_AE_EMIS + NSPC3 = NSPC2 + N_NR_EMIS + + N_GAS_EMIS = N_GC_EMIS + N_NR_EMIS + N_TR_EMIS + + ALLOCATE( SPC_PTEM_TYPE_FAC( N_GAS_EMIS ), STAT = IOS ) + CALL CHECKMEM( IOS, 'SPC_PTEM_TYPE_FAC: ' // TRIM( ELABEL ), PNAME ) + SPC_PTEM_TYPE_FAC = 1.0 + + ALLOCATE( SPC_PTEM_TYPE_MAP( N_GAS_EMIS ), STAT = IOS ) + CALL CHECKMEM( IOS, 'SPC_PTEM_TYPE_MAP: ' // TRIM( ELABEL ), PNAME ) + SPC_PTEM_TYPE_MAP = -1 + + ALLOCATE( PTEM_TYPE_MAP( N_GAS_EMIS ), STAT = IOS ) + CALL CHECKMEM( IOS, 'PTEM_TYPE_MAP: ' // TRIM( ELABEL ), PNAME ) + PTEM_TYPE_MAP = -1 + + N = 0 + + S_OFFSET = 0 + DO S = 1, N_GC_EMIS + IDX = INDEX1( GC_EMIS( S ), NSPC, EM % TABLE( 1, 1 ) ) + IF ( IDX .GT. 0 ) THEN + N = N + 1 + PTEM_TYPE_MAP ( S ) = IDX + SPC_PTEM_TYPE_MAP( S ) = S + SPC_PTEM_TYPE_FAC( S ) = GC_EMIS_FAC( S ) + END IF + END DO + + S_OFFSET = N_GC_EMIS + DO S = 1, N_NR_EMIS + IDX = INDEX1( NR_EMIS( S ), NSPC, EM % TABLE( 1, 1 ) ) + IF ( IDX .GT. 0 ) THEN + N = N + 1 + V = S + S_OFFSET + PTEM_TYPE_MAP ( V ) = IDX + SPC_PTEM_TYPE_MAP( V ) = S + NSPC2 + SPC_PTEM_TYPE_FAC( V ) = NR_EMIS_FAC( S ) + END IF + END DO + + S_OFFSET = S_OFFSET + N_NR_EMIS + DO S = 1, N_TR_EMIS + IDX = INDEX1( TR_EMIS( S ), NSPC, EM % TABLE( 1, 1 ) ) + IF ( IDX .GT. 0 ) THEN + N = N + 1 + V = S + S_OFFSET + PTEM_TYPE_MAP ( V ) = IDX + SPC_PTEM_TYPE_MAP( V ) = S + NSPC3 + SPC_PTEM_TYPE_FAC( V ) = TR_EMIS_FAC( S ) + END IF + END DO + + N_GSPC_TYPE_EMIS = N + +C ... aerosol species ... + + N = 0 + DO S = 1, N_SPC_PTPM + V = PTPM_MAP( S ) + IDX = INDEX1( PMEM_MAP_NAME( V ), NSPC, EM % TABLE( 1, 1 ) ) + IF ( IDX .GT. 0 ) THEN + N = N + 1 + MAP( N ) = S + END IF + END DO + + N_SPC_TYPE_PTPM = N + + ALLOCATE( PTPM_TYPE_MAP( N_SPC_TYPE_PTPM ), STAT = IOS ) + CALL CHECKMEM( IOS, 'PTPM_TYPE_MAP: '// TRIM( ELABEL ), PNAME ) + PTPM_TYPE_MAP = MAP( 1:N_SPC_TYPE_PTPM ) + + DEALLOCATE( MAP ) + + END FUNCTION PTMAP_TYPE_INIT + END MODULE PTMAP diff --git a/src/model/src/aero_depv.F b/src/model/src/aero_depv.F new file mode 100644 index 00000000..f0caefa5 --- /dev/null +++ b/src/model/src/aero_depv.F @@ -0,0 +1,894 @@ + +!------------------------------------------------------------------------! +! The Community Multiscale Air Quality (CMAQ) system software is in ! +! continuous development by various groups and is based on information ! +! from these groups: Federal Government employees, contractors working ! +! within a United States Government contract, and non-Federal sources ! +! including research institutions. These groups give the Government ! +! permission to use, prepare derivative works of, and distribute copies ! +! of their work in the CMAQ system to the public and to permit others ! +! to do so. The United States Environmental Protection Agency ! +! therefore grants similar permission to use the CMAQ system software, ! +! but users are requested to provide copies of derivative works or ! +! products designed to operate in the CMAQ system to the United States ! +! Government without restrictions as to use by others. Software ! +! that is used with the CMAQ system but distributed under the GNU ! +! General Public License or the GNU Lesser General Public License is ! +! subject to their copyright restrictions. ! +!------------------------------------------------------------------------! + + +C RCS file, release, date & time of last delta, author, state, [and locker] +C $Header: /project/yoj/arc/CCTM/src/aero/aero5/aero_depv.F,v 1.12 2012/01/19 13:12:14 yoj Exp $ + +C what(1) key, module and SID; SCCS file; date and time of last delta: +C @(#)aero_depv.F 1.3 /project/mod3/CMAQ/src/ae_depv/aero_depv/SCCS/s.aero_depv.F 18 Jun 1997 12:55:48 + +C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + SUBROUTINE AERO_DEPV ( CGRID, JDATE, JTIME, TSTEP, MOSAIC, VDEP_AE ) + +C----------------------------------------------------------------------- +C aerosol dry deposition routine +C written 4/9/97 by Dr. Francis S. Binkowski +C uses code from modpar and vdvg from the aerosol module. +C This routine uses a single block to hold information +C for the lowest layer. +C NOTES: This version assumes that RA is available on the met file. +c Array structure for vector optimization +C 26 Apr 97 Jeff - many mods +C 13 Dec 97 Jeff - expect uncoupled CGRID, concs as micro-g/m**3, #/m**3 +C +C 1/11/99 David Wong at LM - change NUMCELLS to CELLNUM in the loop index +C FSB 3/17/99 changed to accommodate surface area/second moment and +C encapsulated the actual drydep calculation into a subroutine which +C is attached to this code +C Jeff - Dec 00 - move CGRID_MAP into f90 module +C FSB 12/11/2000. Logic added to allow deposition of particles at their +C "wet" diameters; that is, accounting for the water on the particles. +C This is done by adjusting the third and second moments for the +C presence of water assuming that the geometric standard deviations +C are not changed by this process. This appears to be a very good +C assumption. +C 30 Aug 01 J.Young: Dyn alloc; Use HGRD_DEFN +C Jan 03 J.Young: Change CGRID dimensions, eliminate re-allocations +C 6 Mar 03 J.Young: eliminate a lot of allocate/deallocates +C 7 Aug 03 S.Roselle: updated code for loading the min aero conc array +C 17 Dec 03 S.Roselle: Adjust 2nd and 3rd moments to include SOA, +C without affecting the geometric standard deviations. +C 31 Jan 05 J.Young: dyn alloc - establish both horizontal & vertical +C domain specifications in one module +C 07 Jun 05 P.Bhave: Added code to handle new species in the AE4 +C mechanism: ANAI, ANAJ, ANAK, ACLI, ACLJ, ACLK, ASO4K, AH2OK, +C and ANO3K; look for ASEAS only when using AE3 mechanism +C 30 Jan 08 S.Napelenok & P.Bhave: Added code to handle new SOA species +C in AE5; defined DRY aerosol to include nonvolatile SOA spcs +C 14 Apr 08 J.Kelly: Added code to handle new species ANH4K and SRFCOR. +C Also added code to handle variable coarse mode standard deviation +C in AE5 (no longer fixed at 2.2). +C 08 Sep 08 P.Bhave: Backward compatibility with AE4 mechanisms +C standardized names of all coarse-mode variables +C 19 Apr 10 S.Howard: aero re-engineering for modularity +C 23 Apr 10 J.Young: replace chem mechanism include files with namelists +C 10 Mar 11 S.Howard: Renamed met_data to aeromet_data +C 25 Mar 11 S.Roselle: Replaced I/O API include files with UTILIO_DEFN +C 20 May 11 D.Schwede: Modified for mosaic +C 31 Aug 11 J.Bash: Moved shared mosaic variables to MOSAIC_MOD +C 27 Sep 11 David Wong: replaced all run time dynamic arrays with allocatable +C arrays to avoid run time memory issue +C 08 Jun 12 J.Young: remove full character blank padding for GNU Fortran (GCC) 4.1.2 +C 07 Jul 14 B.Hutzell: replaced mechanism include file(s) with fortran module +C 07 Nov 14 J.Bash: Updated for the ASX_DATA_MOD shared data module. +C May 16 B. Murphy, H. Pye: Updated treatment of aerosol moments +C----------------------------------------------------------------------- + + USE GRID_CONF ! horizontal & vertical domain specifications + USE RXNS_DATA ! chemical mechanism data + USE CGRID_SPCS ! CGRID mechanism species + USE UTILIO_DEFN + USE AERO_DATA ! aero variable data + USE AEROMET_DATA ! Includes CONST.EXT + USE Mosaic_Mod, Only: ADEPVJ ! Shared mosaic variables + Use LSM_Mod, Only: N_LUFRAC + USE ASX_DATA_MOD, Only: Met_Data + + IMPLICIT NONE + +C Includes: + + INCLUDE SUBST_FILES_ID ! file name parameters + +C Arguments + REAL, POINTER :: CGRID( :,:,:,: ) + INTEGER, INTENT( IN ) :: JDATE ! current model date , coded YYYYDDD + INTEGER, INTENT( IN ) :: JTIME ! current model time , coded HHMMSS + INTEGER, INTENT( IN ) :: TSTEP ! model time step, coded HHMMSS + LOGICAL, INTENT( IN ) :: MOSAIC ! use mosaic option + REAL, INTENT( OUT ) :: VDEP_AE( :,:,: ) ! surrogate deposition velocities [ m s**-1 ] + +C Parameters + REAL, PARAMETER :: T0 = 288.15 ! [ K ] ! starting standard surface temp. + REAL, PARAMETER :: TWO3 = 2.0 / 3.0 + INTEGER, PARAMETER :: N_AE_DEP_SPC = 9 ! no. of surrogates for aerosol dry dep velocities + +C set up species dimension and indices for deposition velocity internal array VDEP + INTEGER, PARAMETER :: VDNATK = 1, ! Aitken mode number + & VDNACC = 2, ! accumulation mode number + & VDNCOR = 3, ! coarse mode number + & VDMATK = 4, ! Aitken mode mass + & VDMACC = 5, ! accumulation mode mass + & VDMCOR = 6, ! coarse mode mass + & VDSATK = 7, ! Aitken mode surface area + & VDSACC = 8, ! accumulation mode surface area + & VDSCOR = 9 ! coarse mode surface area + +C Local variables: + + CHARACTER( 16 ) :: VDAE_NAME( N_AE_DEP_SPC )! dep vel surrogate name table + DATA VDAE_NAME( 1 ) / 'VNUMATKN' / + DATA VDAE_NAME( 2 ) / 'VNUMACC ' / + DATA VDAE_NAME( 3 ) / 'VNUMCOR ' / + DATA VDAE_NAME( 4 ) / 'VMASSI ' / + DATA VDAE_NAME( 5 ) / 'VMASSJ ' / + DATA VDAE_NAME( 6 ) / 'VMASSC ' / + DATA VDAE_NAME( 7 ) / 'VSRFATKN' / + DATA VDAE_NAME( 8 ) / 'VSRFACC ' / + DATA VDAE_NAME( 9 ) / 'VSRFCOR ' / + + INTEGER, ALLOCATABLE, SAVE :: DEPV_SUR( : ) ! pointer to surrogate + +C Meteorological variables + + CHARACTER( 16 ), SAVE :: AE_VRSN ! Aerosol version name + + INTEGER, SAVE :: NCELLS ! number of cells per layer + + REAL, ALLOCATABLE, SAVE :: XXLSGAT( :,: ) ! log of standard deviation + REAL, ALLOCATABLE, SAVE :: XXLSGAC( :,: ) + REAL, ALLOCATABLE, SAVE :: XXLSGCO( :,: ) + + REAL, ALLOCATABLE, SAVE :: DGATK( :,: ) ! geometric mean diameter + REAL, ALLOCATABLE, SAVE :: DGACC( :,: ) + REAL, ALLOCATABLE, SAVE :: DGCOR( :,: ) + + REAL, ALLOCATABLE, SAVE :: PDENSAT( :,: ) ! particle density + REAL, ALLOCATABLE, SAVE :: PDENSAC( :,: ) + REAL, ALLOCATABLE, SAVE :: PDENSCO( :,: ) + + REAL, ALLOCATABLE, SAVE :: XLM( :,: ) ! mean free path [ m ] + REAL, ALLOCATABLE, SAVE :: AMU( :,: ) ! dynamic viscosity [ kg m**-1 s**-1 ] + + REAL, ALLOCATABLE, SAVE :: VDEP( :,:,: ) ! deposition velocity [ m/s ] + REAL, ALLOCATABLE, SAVE :: VDEPJ( :,:,:,: ) ! deposition velocity [ m/s ] + + REAL M3_WET, M3SUBT, M3_DRY + REAL M2_WET, M2_DRY + + LOGICAL, SAVE :: FIRSTIME = .TRUE. + INTEGER, SAVE :: LOGDEV ! unit number for the log file + CHARACTER( 16 ), SAVE :: PNAME = 'AERO_DEPV' + CHARACTER( 16 ) :: VNAME ! varable name + CHARACTER( 96 ) :: XMSG = ' ' + + INTEGER C, R, V, N, J ! loop counters + INTEGER SPC, S ! species loop counter + INTEGER ALLOCSTAT + + INTERFACE + SUBROUTINE GETDEP_V ( XLM, AMU, DGATK, DGACC, DGCOR, + & XXLSGAT, XXLSGAC, XXLSGCO, + & PDENSAT, PDENSAC, PDENSCO, + & VDEP, VDEPJ, MOSAIC ) + REAL, INTENT( IN ) :: XLM( :,: ) ! atmospheric mean free path [ m ] + REAL, INTENT( IN ) :: AMU( :,: ) ! atmospheric dynamic viscosity [ kg/(m s) ] + REAL, INTENT( IN ) :: DGATK( :,: ) ! nuclei mode geometric mean diameter [ m ] + REAL, INTENT( IN ) :: DGACC( :,: ) ! accumulation geometric mean diameter [ m ] + REAL, INTENT( IN ) :: DGCOR( :,: ) ! coarse mode geometric mean diameter [ m ] + REAL, INTENT( IN ) :: XXLSGAT( :,: ) ! Aitken mode + REAL, INTENT( IN ) :: XXLSGAC( :,: ) ! accumulation mode + REAL, INTENT( IN ) :: XXLSGCO( :,: ) ! coarse mode + REAL, INTENT( IN ) :: PDENSAT( :,: ) ! average particle density in nuclei mode + REAL, INTENT( IN ) :: PDENSAC( :,: ) ! average particle density in accumulation mode + REAL, INTENT( IN ) :: PDENSCO( :,: ) ! average particle density in coarse mode + REAL, INTENT( OUT ) :: VDEP( :,:,: ) ! deposition velocity [ m/s ] + REAL, INTENT( OUT ) :: VDEPJ( :,:,:,: ) ! deposition velocity [ m/s ] for each land use category + LOGICAL, INTENT( IN ) :: MOSAIC + END SUBROUTINE GETDEP_V + END INTERFACE + +C----------------------------------------------------------------------- + + IF ( FIRSTIME ) THEN + FIRSTIME = .FALSE. + LOGDEV = INIT3() + + NCELLS = NCOLS * NROWS + +C Allocate arrays + ALLOCATE( XXLSGAT( NCOLS,NROWS ), + & XXLSGAC( NCOLS,NROWS ), + & XXLSGCO( NCOLS,NROWS ), STAT = ALLOCSTAT ) + IF ( ALLOCSTAT .NE. 0 ) THEN + XMSG = 'Failure allocating XXLSGAT, XXLSGAC or XXLSGCO' + CALL M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT1 ) + END IF + + ALLOCATE( DGATK( NCOLS,NROWS ), + & DGACC( NCOLS,NROWS ), + & DGCOR( NCOLS,NROWS ), STAT = ALLOCSTAT ) + IF ( ALLOCSTAT .NE. 0 ) THEN + XMSG = 'Failure allocating DGATK, DGACC or DGCOR' + CALL M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT1 ) + END IF + + ALLOCATE( PDENSAT( NCOLS,NROWS ), + & PDENSAC( NCOLS,NROWS ), + & PDENSCO( NCOLS,NROWS ), STAT = ALLOCSTAT ) + IF ( ALLOCSTAT .NE. 0 ) THEN + XMSG = 'Failure allocating PDENSAT, PDENSAC or PDENSCO' + CALL M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT1 ) + END IF + + ALLOCATE( XLM( NCOLS,NROWS ), + & AMU( NCOLS,NROWS ), STAT = ALLOCSTAT ) + IF ( ALLOCSTAT .NE. 0 ) THEN + XMSG = 'Failure allocating XLM or AMU' + CALL M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT1 ) + END IF + + ALLOCATE ( VDEP( NCOLS,NROWS,N_AE_DEP_SPC ), STAT = ALLOCSTAT ) + IF ( ALLOCSTAT .NE. 0 ) THEN + XMSG = 'Failure allocating VDEP' + CALL M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT1 ) + END IF + + ALLOCATE ( VDEPJ( N_LUFRAC,NCOLS,NROWS,N_AE_DEP_SPC ), STAT = ALLOCSTAT ) + IF ( ALLOCSTAT .NE. 0 ) THEN + XMSG = 'Failure allocating VDEPJ' + CALL M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT1 ) + END IF + IF ( MOSAIC ) THEN + ALLOCATE ( ADEPVJ( N_LUFRAC,N_AE_DEPV,NCOLS,NROWS ), STAT = ALLOCSTAT ) + IF ( ALLOCSTAT .NE. 0 ) THEN + XMSG = 'Failure allocating ADEV' + CALL M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT1 ) + END IF + END IF + + ALLOCATE ( DEPV_SUR( N_AE_DEPV ), STAT = ALLOCSTAT ) + IF ( ALLOCSTAT .NE. 0 ) THEN + XMSG = 'Failure allocating DEPV_SUR' + CALL M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT1 ) + END IF + +C Set the dep vel surrogate pointers + DO V = 1, N_AE_DEPV + N = INDEX1( AE_DEPV( V ), N_AE_DEP_SPC, VDAE_NAME ) + IF ( N .NE. 0 ) THEN + DEPV_SUR( V ) = N + ELSE + XMSG = 'Could not find ' // AE_DEPV( V ) // ' in aerosol' // + & ' surrogate table. >>> Dep vel set to zero <<< ' + CALL M3WARN( PNAME, JDATE, JTIME, XMSG ) + DEPV_SUR( V ) = 0 + END IF + END DO + + END IF ! FIRSTIME + + IF ( N_AE_SPC .LE. 0 ) RETURN + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +C Put the grid cell physical data in the block arrays +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + + DO R = 1, MY_NROWS + DO C = 1, MY_NCOLS + +C *** Set meteorological data for the grid cell. + AIRDENS = Met_Data%DENS1( C,R ) + AIRTEMP = Met_Data%TEMP2( C,R ) + AIRPRES = Met_Data%PRSFC( C,R ) + +C *** extract grid cell concentrations of aero species from CGRID +C into aerospc_conc in aero_data module +C Also determines second moment from surface area and adds wet +C species + CALL EXTRACT_AERO( CGRID( C,R,1,: ), .TRUE. ) + +C *** Calculate geometric mean diameters and standard deviations of the +C "wet" size distribution + CALL GETPAR( .FALSE. ) + +C Save getpar values to arrays + XXLSGAT( C,R ) = AEROMODE_LNSG( 1 ) + XXLSGAC( C,R ) = AEROMODE_LNSG( 2 ) + XXLSGCO( C,R ) = AEROMODE_LNSG( 3 ) + + DGATK( C,R ) = AEROMODE_DIAM( 1 ) + DGACC( C,R ) = AEROMODE_DIAM( 2 ) + DGCOR( C,R ) = AEROMODE_DIAM( 3 ) + + PDENSAT( C,R ) = AEROMODE_DENS( 1 ) + PDENSAC( C,R ) = AEROMODE_DENS( 2 ) + PDENSCO( C,R ) = AEROMODE_DENS( 3 ) + +C Calculate mean free path [ m ]: + XLM( C,R ) = 6.6328E-8 * STDATMPA * AIRTEMP / ( T0 * AIRPRES ) + +C *** Calcualte dynamic viscosity [ kg m**-1 s**-1 ]: + AMU( C,R ) = 1.458E-6 * AIRTEMP * SQRT( AIRTEMP ) + & / ( AIRTEMP + 110.4 ) + + END DO ! Column LOOP + END DO ! Row LOOP + +C *** get dry deposition velocities: + + CALL GETDEP_V ( XLM, AMU, DGATK, DGACC, DGCOR, + & XXLSGAT, XXLSGAC, XXLSGCO, + & PDENSAT, PDENSAC, PDENSCO, + & VDEP, VDEPJ, MOSAIC ) + +C Return dry deposition velocities for aerosols (first layer only). + + DO R = 1, NROWS + DO C = 1, NCOLS + DO V = 1, N_AE_DEPV + IF ( DEPV_SUR( V ) .GT. 0 ) THEN + VDEP_AE( V,C,R ) = VDEP( C,R,DEPV_SUR( V ) ) + ELSE + VDEP_AE( V,C,R ) = 0.0 + END IF + END DO + END DO + END DO + + IF ( MOSAIC ) THEN + DO R = 1, NROWS + DO C = 1, NCOLS + DO V = 1, N_AE_DEPV + DO J = 1, N_LUFRAC + IF ( DEPV_SUR( V ) .GT. 0 ) THEN + ADEPVJ( J,V,C,R ) = VDEPJ( J,C,R,DEPV_SUR( V ) ) + ELSE + ADEPVJ( J,V,C,R ) = 0.0 + END IF + END DO + END DO + END DO + END DO + END IF + + RETURN + END SUBROUTINE AERO_DEPV + +C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + SUBROUTINE GETDEP_V ( XLM, AMU, DGATK, DGACC, DGCOR, + & XXLSGAT, XXLSGAC, XXLSGCO, + & PDENSAT, PDENSAC, PDENSCO, + & VDEP, VDEPJ, MOSAIC ) + +C *** Calculate deposition velocity for Aitken, accumulation, and +C coarse modes. +C Reference: +C Binkowski F. S., and U. Shankar, The regional particulate +C model 1. Model description and preliminary results. +C J. Geophys. Res., 100, D12, 26191-26209, 1995. + +C May 05 D.Schwede: added impaction term to coarse mode dry deposition +C 25 May 05 J.Pleim: Updated dry dep velocity calculation for aerosols +C to Venkatram and Pleim (1999) +C 20 Jul 05 J.Pleim: Changed impaction term using modal integration of +C Stokes**2 / 400 (Giorgi, 1986, JGR) +C 14 Apr 08 J.Kelly: Added code to calculate deposition velocity of +C coarse surface area and to account for variable +C standard deviation of the coarse mode. +C 08 Sep 08 P.Bhave: Backward compatibility with AE4 mechanisms +C standardized names of all coarse-mode variables +C----------------------------------------------------------------------- + + USE LSM_MOD, Only: n_lufrac ! to get n_lufrac + USE AEROMET_DATA ! Includes CONST.EXT + USE ASX_DATA_MOD, Only: Met_Data, Mosaic_Data, Grid_Data + + IMPLICIT NONE + +C *** input arguments + +C atmospheric properties + REAL, INTENT( IN ) :: XLM( :,: ) ! atmospheric mean free path [ m ] + REAL, INTENT( IN ) :: AMU( :,: ) ! atmospheric dynamic viscosity [ kg/(m s) ] + +C aerosol properties: + +C modal diameters: [ m ] + REAL, INTENT( IN ) :: DGATK( :,: ) ! nuclei mode geometric mean diameter [ m ] + REAL, INTENT( IN ) :: DGACC( :,: ) ! accumulation geometric mean diameter [ m ] + REAL, INTENT( IN ) :: DGCOR( :,: ) ! coarse mode geometric mean diameter [ m ] + +C log of modal geometric standard deviations + REAL, INTENT( IN ) :: XXLSGAT( :,: ) ! Aitken mode + REAL, INTENT( IN ) :: XXLSGAC( :,: ) ! accumulation mode + REAL, INTENT( IN ) :: XXLSGCO( :,: ) ! coarse mode + +C average modal particle densities [ kg/m**3 ] + REAL, INTENT( IN ) :: PDENSAT( :,: ) ! average particle density in nuclei mode + REAL, INTENT( IN ) :: PDENSAC( :,: ) ! average particle density in accumulation mode + REAL, INTENT( IN ) :: PDENSCO( :,: ) ! average particle density in coarse mode + +C Mosaic variables + LOGICAL, INTENT( IN ) :: MOSAIC + +C *** output arguments + + REAL, INTENT( OUT ) :: VDEP( :,:,: ) ! deposition velocity [ m/s ] + REAL, INTENT( OUT ) :: VDEPJ( :,:,:,: ) ! deposition velocity [ m/s ] for each land use category + +C *** array indices hardcoded to match SUBROUTINE AERO_DEPV + INTEGER, PARAMETER :: VDNATK = 1, ! Aitken mode number + & VDNACC = 2, ! accumulation mode number + & VDNCOR = 3, ! coarse mode number + & VDMATK = 4, ! Aitken mode mass + & VDMACC = 5, ! accumulation mode mass + & VDMCOR = 6, ! coarse mode mass + & VDSATK = 7, ! Aitken mode surface area + & VDSACC = 8, ! accumulation mode surface area + & VDSCOR = 9 ! coarse mode surface area + +C modal Knudsen numbers + REAL KNATK ! Aitken mode Knudsen number + REAL KNACC ! accumulation " + REAL KNCOR ! coarse mode + +C modal particle diffusivities for number, 2nd, and 3rd moment, or mass: + REAL DCHAT0N, DCHAT0A, DCHAT0C + REAL DCHAT2N, DCHAT2A, DCHAT2C + REAL DCHAT3N, DCHAT3A, DCHAT3C + +C modal sedimentation velocities for number, 2nd, and 3rd moment, or mass: + REAL VGHAT0N, VGHAT0A, VGHAT0C + REAL VGHAT2N, VGHAT2A, VGHAT2C + REAL VGHAT3N, VGHAT3A, VGHAT3C + + INTEGER NCELL, J, C, R + + REAL DCONST1, DCONST1N, DCONST1A, DCONST1C + REAL DCONST2, DCONST3N, DCONST3A, DCONST3C + REAL SC0N, SC0A, SC0C ! Schmidt numbers for number + REAL SC2N, SC2A, SC2C ! Schmidt numbers for 2ND MOMENT + REAL SC3N, SC3A, SC3C ! Schmidt numbers for 3rd moment + REAL(8) STOKEN, STOKEA, STOKEC ! Stokes numbers for each mode + REAL RD0N, RD0A, RD0C ! canopy resistance for number + REAL RD2N, RD2A, RD2C ! canopy resistance for 2nd moment + REAL RD3N, RD3A, RD3C ! canopy resisteance for 3rd moment + REAL UTSCALE ! scratch function of USTAR and WSTAR + REAL NU ! kinematic viscosity [ m**2 s**-1 ] + REAL USTFAC ! scratch function of USTAR, NU, and GRAV + REAL TWOXLM ! 2 X atmospheric mean free path + + REAL, PARAMETER :: BHAT = 1.246 ! Constant from Cunningham slip correction + REAL, PARAMETER :: THREEPI = 3.0 * PI + REAL, PARAMETER :: TWO3 = 2.0 / 3.0 + +C Scalar variables for VARIABLE standard deviations. + + REAL L2SGAT, L2SGAC ! see usage + REAL L2SGCO + + REAL EAT1 ! Aitken mode exp( log^2( sigmag )/8 ) + REAL EAC1 ! accumulation mode exp( log^2( sigmag )/8 ) + REAL ECO1 ! coarse mode exp( log^2( sigmag )/8 ) + + REAL ESAT04 ! Aitken " **4 + REAL ESAC04 ! accumulation " + REAL ESCO04 ! coarse " + + REAL ESAT08 ! Aitken " **8 + REAL ESAC08 ! accumulation " + REAL ESCO08 ! coarse " + + REAL ESAT12 ! Aitken " **12 + REAL ESAC12 ! accumulation " + REAL ESCO12 ! coarse " + + REAL ESAT16 ! Aitken " **16 + REAL ESAC16 ! accumulation " + REAL ESCO16 ! coarse " + + REAL ESAT20 ! Aitken " **20 + REAL ESAC20 ! accumulation " + REAL ESCO20 ! coarse " + + REAL ESAT28 ! Aitken " **28 + REAL ESAC28 ! accumulation " + REAL ESCO28 ! coarse " + + REAL ESAT32 ! Aitken " **32 + REAL ESAC32 ! accumulation " + REAL ESCO32 ! coarse " + + REAL ESAT36 ! Aitken " **36 + REAL ESAC36 ! accumulation " + REAL ESCO36 ! coarse " + + REAL ESAT48 ! Aitken " **48 + REAL ESAC48 ! accumulation " + REAL ESCO48 ! coarse " + + REAL ESAT64 ! Aitken " **64 + REAL ESAC64 ! accumulation " + REAL ESCO64 ! coarse " + + REAL ESAT128 ! Aitken " **128 + REAL ESAC128 ! accumulation " + REAL(8) ESCO128 ! coarse " + + REAL ESAT160 ! Aitken " **160 + REAL ESAC160 ! accumulation " + REAL ESCO160 ! coarse " + + REAL ESATM12 ! Aitken " **(-12) + REAL ESACM12 ! accumulation " + REAL ESCOM12 ! coarse " + + REAL ESATM16 ! Aitken " **(-16) + REAL ESACM16 ! accumulation " + REAL ESCOM16 ! coarse " + + REAL ESATM20 ! Aitken " **(-20) + REAL ESACM20 ! accumulation " + REAL ESCOM20 ! coarse " + + REAL ESATM32 ! Aitken " **(-32) + REAL ESACM32 ! accumulation " + REAL ESCOM32 ! coarse " + + REAL(8) EIM ! Impaction efficiency + +C----------------------------------------------------------------------- + + VDEP = 0.0 ! array assignment + IF ( MOSAIC ) THEN + VDEPJ = 0.0 ! array assignment + END IF + + DO R = 1, SIZE( Met_Data%TEMP2, 2 ) + DO C = 1, SIZE( Met_Data%TEMP2, 1 ) + +C *** Calculate Knudsen numbers + + TWOXLM = XLM( C,R ) + XLM( C,R ) + KNATK = TWOXLM / DGATK( C,R ) + KNACC = TWOXLM / DGACC( C,R ) + KNCOR = TWOXLM / DGCOR( C,R ) + +C *** Calculate functions of variable standard deviation. + + L2SGAT = XXLSGAT( C,R ) * XXLSGAT( C,R ) + L2SGAC = XXLSGAC( C,R ) * XXLSGAC( C,R ) + L2SGCO = XXLSGCO( C,R ) * XXLSGCO( C,R ) + + EAT1 = EXP( 0.125 * L2SGAT ) + EAC1 = EXP( 0.125 * L2SGAC ) + ECO1 = EXP( 0.125 * L2SGCO ) + + ESAT04 = EAT1 ** 4 + ESAC04 = EAC1 ** 4 + ESCO04 = ECO1 ** 4 + + ESAT08 = ESAT04 * ESAT04 + ESAC08 = ESAC04 * ESAC04 + ESCO08 = ESCO04 * ESCO04 + + ESAT12 = ESAT04 * ESAT08 + ESAC12 = ESAC04 * ESAC08 + ESCO12 = ESCO04 * ESCO08 + + ESAT16 = ESAT08 * ESAT08 + ESAC16 = ESAC08 * ESAC08 + ESCO16 = ESCO08 * ESCO08 + + ESAT20 = ESAT16 * ESAT04 + ESAC20 = ESAC16 * ESAC04 + ESCO20 = ESCO16 * ESCO04 + + ESAT28 = ESAT20 * ESAT08 + ESAC28 = ESAC20 * ESAC08 + ESCO28 = ESCO20 * ESCO08 + + ESAT32 = ESAT16 * ESAT16 + ESAC32 = ESAC16 * ESAC16 + ESCO32 = ESCO16 * ESCO16 + + ESAT36 = ESAT16 * ESAT20 + ESAC36 = ESAC16 * ESAC20 + ESCO36 = ESCO16 * ESCO20 + + ESAT48 = ESAT36 * ESAT12 + ESAC48 = ESAC36 * ESAC12 + ESCO48 = ESCO36 * ESCO12 + + ESAT64 = ESAT32 * ESAT32 + ESAC64 = ESAC32 * ESAC32 + ESCO64 = ESCO32 * ESCO32 + + ESAT128 = ESAT64 * ESAT64 + ESAC128 = ESAC64 * ESAC64 + ESCO128 = ESCO64 * ESCO64 + + ESAT160 = ESAT128* ESAT32 + ESAC160 = ESAC128* ESAC32 + ESCO160 = ESCO128* ESCO32 + +C *** calculate inverses: + + ESATM12 = 1.0 / ESAT12 + ESACM12 = 1.0 / ESAC12 + ESCOM12 = 1.0 / ESCO12 + + ESATM16 = 1.0 / ESAT16 + ESACM16 = 1.0 / ESAC16 + ESCOM16 = 1.0 / ESCO16 + + ESATM20 = 1.0 / ESAT20 + ESACM20 = 1.0 / ESAC20 + ESCOM20 = 1.0 / ESCO20 + + ESATM32 = 1.0 / ESAT32 + ESACM32 = 1.0 / ESAC32 + ESCOM32 = 1.0 / ESCO32 + + DCONST1 = BOLTZMANN * Met_Data%TEMP2( C,R ) / ( THREEPI * AMU( C,R ) ) + DCONST1N = DCONST1 / DGATK( C,R ) + DCONST1A = DCONST1 / DGACC( C,R ) + DCONST1C = DCONST1 / DGCOR( C,R ) + DCONST2 = GRAV / ( 18.0 * AMU( C,R ) ) + DCONST3N = DCONST2 * PDENSAT( C,R ) * DGATK( C,R ) * DGATK( C,R ) + DCONST3A = DCONST2 * PDENSAC( C,R ) * DGACC( C,R ) * DGACC( C,R ) + DCONST3C = DCONST2 * PDENSCO( C,R ) * DGCOR( C,R ) * DGCOR( C,R ) + +C i-mode + DCHAT0N = DCONST1N * ( ESAT04 + BHAT * KNATK * ESAT16 ) + DCHAT2N = DCONST1N * ( ESATM12 + BHAT * KNATK * ESATM16 ) + DCHAT3N = DCONST1N * ( ESATM20 + BHAT * KNATK * ESATM32 ) + VGHAT0N = DCONST3N * ( ESAT16 + BHAT * KNATK * ESAT04 ) + VGHAT2N = DCONST3N * ( ESAT48 + BHAT * KNATK * ESAT20 ) + VGHAT3N = DCONST3N * ( ESAT64 + BHAT * KNATK * ESAT28 ) + +C j-mode + DCHAT0A = DCONST1A * ( ESAC04 + BHAT * KNACC * ESAC16 ) + DCHAT2A = DCONST1A * ( ESACM12 + BHAT * KNACC * ESACM16 ) + DCHAT3A = DCONST1A * ( ESACM20 + BHAT * KNACC * ESACM32 ) + VGHAT0A = DCONST3A * ( ESAC16 + BHAT * KNACC * ESAC04 ) + VGHAT2A = DCONST3A * ( ESAC48 + BHAT * KNACC * ESAC20 ) + VGHAT3A = DCONST3A * ( ESAC64 + BHAT * KNACC * ESAC28 ) + +C coarse mode + DCHAT0C = DCONST1C * ( ESCO04 + BHAT * KNCOR * ESCO16 ) + DCHAT2C = DCONST1C * ( ESCOM12 + BHAT * KNCOR * ESCOM16 ) + DCHAT3C = DCONST1C * ( ESCOM20 + BHAT * KNCOR * ESCOM32 ) + VGHAT0C = DCONST3C * ( ESCO16 + BHAT * KNCOR * ESCO04 ) + VGHAT2C = DCONST3C * ( ESCO48 + BHAT * KNCOR * ESCO20 ) + VGHAT3C = DCONST3C * ( ESCO64 + BHAT * KNCOR * ESCO28 ) + +C now calculate the deposition velocities + + NU = AMU( C,R ) / Met_Data%DENS1( C,R ) + USTFAC = Met_Data%USTAR( C,R ) ** 2 / ( GRAV * NU ) + STOKEN = DCONST3N * USTFAC + STOKEA = DCONST3A * USTFAC + STOKEC = DCONST3C * USTFAC + UTSCALE = Met_Data%USTAR( C,R ) + & + 0.24 * Met_Data%WSTAR( C,R )**2 + & / Met_Data%USTAR( C,R ) + +C first do 0th moment for the deposition of number + +C Aitken mode + SC0N = NU / DCHAT0N + EIM = STOKEN ** 2 / 400.0 * ESAT64 + EIM = MIN( EIM, 1.0 ) + RD0N = 1.0 / ( UTSCALE * ( SC0N ** ( -TWO3 ) + EIM ) ) + + VDEP( C,R,VDNATK ) = VGHAT0N + & / ( 1.0 - EXP( -VGHAT0N * ( Met_Data%RA( C,R ) + RD0N ) ) ) + +C accumulation mode + SC0A = NU / DCHAT0A + EIM = STOKEA ** 2 / 400.0 * ESAC64 + EIM = MIN( EIM, 1.0 ) + RD0A = 1.0 / ( UTSCALE * ( SC0A ** ( -TWO3 ) + EIM ) ) + + VDEP( C,R,VDNACC ) = VGHAT0A + & / ( 1.0 - EXP( -VGHAT0A * ( Met_Data%RA( C,R ) + RD0A ) ) ) + +C coarse mode + SC0C = NU / DCHAT0C + EIM = STOKEC ** 2 / 400.0 * ESCO64 + EIM = MIN( EIM, 1.0 ) + RD0C = 1.0 / ( UTSCALE * ( SC0C ** ( -TWO3 ) + EIM ) ) + + VDEP( C,R,VDNCOR ) = VGHAT0C + & / ( 1.0 - EXP( -VGHAT0C * ( Met_Data%RA( C,R ) + RD0C ) ) ) + +C now do 2nd moment for the deposition of surface area + +C Aitken mode + SC2N = NU / DCHAT2N + EIM = STOKEN ** 2 / 400.0 * ESAT128 + EIM = MIN( EIM, 1.0 ) + RD2N = 1.0 / ( UTSCALE * ( SC2N ** ( -TWO3 ) + EIM ) ) + + VDEP( C,R,VDSATK ) = VGHAT2N + & / ( 1.0 - EXP( -VGHAT2N * ( Met_Data%RA( C,R ) + RD2N ) ) ) + +C accumulation mode + SC2A = NU / DCHAT2A + EIM = STOKEA ** 2 / 400.0 * ESAC128 + EIM = MIN( EIM, 1.0 ) + RD2A = 1.0 / ( UTSCALE * ( SC2A ** ( -TWO3 ) + EIM ) ) + + VDEP( C,R,VDSACC ) = VGHAT2A + & / ( 1.0 - EXP( -VGHAT2A * ( Met_Data%RA( C,R ) + RD2A ) ) ) + +C coarse mode + SC2C = NU / DCHAT2C + EIM = STOKEC ** 2 / 400.0 * ESCO128 + EIM = MIN( EIM, 1.0 ) + RD2C = 1.0 / ( UTSCALE * ( SC2C ** ( -TWO3 ) + EIM ) ) + + VDEP( C,R,VDSCOR ) = VGHAT2C + & / ( 1.0 - EXP( -VGHAT2C * ( Met_Data%RA( C,R ) + RD2C ) ) ) + +C now do 3rd moment for the deposition of mass + +C Aitken mode + SC3N = NU / DCHAT3N + EIM = STOKEN ** 2 / 400.0 * ESAT160 + EIM = MIN( EIM, 1.0 ) + RD3N = 1.0 / ( UTSCALE * ( SC3N ** ( -TWO3 ) + EIM ) ) + + VDEP( C,R,VDMATK ) = VGHAT3N + & / ( 1.0 - EXP( -VGHAT3N * ( Met_Data%RA( C,R ) + RD3N ) ) ) + +C accumulation mode + SC3A = NU / DCHAT3A + EIM = STOKEA ** 2 / 400.0 * ESAC160 + EIM = MIN( EIM, 1.0 ) + RD3A = 1.0 / ( UTSCALE * ( SC3A ** ( -TWO3 ) + EIM ) ) + + VDEP( C,R,VDMACC ) = VGHAT3A + & / ( 1.0 - EXP( -VGHAT3A * ( Met_Data%RA( C,R ) + RD3A ) ) ) + +C coarse mode + SC3C = NU / DCHAT3C + EIM = STOKEC ** 2 / 400.0 * ESCO160 + EIM = MIN( EIM, 1.0 ) + RD3C = 1.0 / ( UTSCALE * ( SC3C ** ( -TWO3 ) + EIM ) ) + + VDEP( C,R,VDMCOR ) = VGHAT3C + & / ( 1.0 - EXP( -VGHAT3C * ( Met_Data%RA( C,R ) + RD3C ) ) ) + +C Do mosaic calculations - essentially a repeat of the above, using the mosaic vars + IF ( MOSAIC ) THEN + + DO J = 1, N_LUFRAC + IF ( Mosaic_Data%RA( C,R,J ) .GT. 0.0 .And. Grid_Data%LUFRAC( C,R,J ) .GT. 0.0 ) THEN + +C now calculate the deposition velocities + + NU = AMU( C,R ) / Met_Data%DENS1( C,R ) + USTFAC = Mosaic_Data%USTAR( C,R,J )**2 / ( GRAV * NU ) + STOKEN = DCONST3N * USTFAC + STOKEA = DCONST3A * USTFAC + STOKEC = DCONST3C * USTFAC + UTSCALE = Mosaic_Data%USTAR( C,R,J ) + & + 0.24 * Met_Data%WSTAR( C,R )**2 + & / Mosaic_Data%USTAR( C,R,J ) + +C first do 0th moment for the deposition of number + +C Aitken mode + + SC0N = NU / DCHAT0N + EIM = STOKEN ** 2 / 400.0 * ESAT64 + EIM = MIN( EIM, 1.0 ) + RD0N = 1.0 / ( UTSCALE * ( SC0N ** ( -TWO3 ) + EIM ) ) + + VDEPJ( J,C,R,VDNATK ) = VGHAT0N + & / ( 1.0 - EXP(-VGHAT0N * ( Mosaic_Data%RA( C,R,J ) + RD0N ) ) ) + +C accumulation mode + + SC0A = NU / DCHAT0A + EIM = STOKEA ** 2 / 400.0 * ESAC64 + EIM = MIN( EIM, 1.0 ) + RD0A = 1.0 / ( UTSCALE * ( SC0A ** ( -TWO3 ) + EIM ) ) + + VDEPJ( J,C,R,VDNACC ) = VGHAT0A + & / ( 1.0 - EXP(-VGHAT0A * ( Mosaic_Data%RA( C,R,J ) + RD0A ) ) ) + +C coarse mode + + SC0C = NU / DCHAT0C + EIM = STOKEC ** 2 / 400.0 * ESCO64 + EIM = MIN( EIM, 1.0 ) + RD0C = 1.0 / ( UTSCALE * ( SC0C ** ( -TWO3 ) + EIM ) ) + + VDEPJ( J,C,R,VDNCOR ) = VGHAT0C + & / ( 1.0 - EXP(-VGHAT0C * ( Mosaic_Data%RA( C,R,J ) + RD0C ) ) ) + +C now do 2nd moment for the deposition of surface area + +C Aitken mode + + SC2N = NU / DCHAT2N + EIM = STOKEN ** 2 / 400.0 * ESAT128 + EIM = MIN( EIM, 1.0 ) + RD2N = 1.0 / ( UTSCALE * ( SC2N ** ( -TWO3 ) + EIM ) ) + + VDEPJ( J,C,R,VDSATK ) = VGHAT2N + & / ( 1.0 - EXP(-VGHAT2N * ( Mosaic_Data%RA( C,R,J ) + RD2N ) ) ) + +C accumulation mode + + SC2A = NU / DCHAT2A + EIM = STOKEA ** 2 / 400.0 * ESAC128 + EIM = MIN( EIM, 1.0 ) + RD2A = 1.0 / ( UTSCALE * ( SC2A ** ( -TWO3 ) + EIM ) ) + + VDEPJ( J,C,R,VDSACC ) = VGHAT2A + & / ( 1.0 - EXP(-VGHAT2A * ( Mosaic_Data%RA( C,R,J ) + RD2A ) ) ) + +C coarse mode + + SC2C = NU / DCHAT2C + EIM = STOKEC ** 2 / 400.0 * ESCO128 + EIM = MIN( EIM, 1.0 ) + RD2C = 1.0 / ( UTSCALE * ( SC2C ** ( -TWO3 ) + EIM ) ) + + VDEPJ( J,C,R,VDSCOR ) = VGHAT2C + & / ( 1.0 - EXP(-VGHAT2C * ( Mosaic_Data%RA( C,R,J ) + RD2C ) ) ) + +C now do 3rd moment for the deposition of mass + +C Aitken mode + + SC3N = NU / DCHAT3N + EIM = STOKEN ** 2 / 400.0 * ESAT160 + EIM = MIN( EIM, 1.0 ) + RD3N = 1.0 / ( UTSCALE * ( SC3N ** ( -TWO3 ) + EIM ) ) + + VDEPJ( J,C,R,VDMATK ) = VGHAT3N + & / ( 1.0 - EXP(-VGHAT3N * ( Mosaic_Data%RA( C,R,J ) + RD3N ) ) ) + +C accumulation mode + + SC3A = NU / DCHAT3A + EIM = STOKEA ** 2 / 400.0 * ESAC160 + EIM = MIN( EIM, 1.0 ) + RD3A = 1.0 / ( UTSCALE * ( SC3A ** ( -TWO3 ) + EIM ) ) + + VDEPJ( J,C,R,VDMACC ) = VGHAT3A + & / ( 1.0 - EXP(-VGHAT3A * ( Mosaic_Data%RA( C,R,J ) + RD3A ) ) ) + +C coarse mode + + SC3C = NU / DCHAT3C + EIM = STOKEC ** 2 / 400.0 * ESCO160 + EIM = MIN( EIM, 1.0 ) + RD3C = 1.0 / ( UTSCALE * ( SC3C ** ( -TWO3 ) + EIM ) ) + + VDEPJ( J,C,R,VDMCOR ) = VGHAT3C + & / ( 1.0 - EXP(-VGHAT3C * ( Mosaic_Data%RA( C,R,J ) + RD3C ) ) ) + + END IF ! RA > 0 + + END DO ! n_lufrac + + END IF ! mosaic + END DO ! end loop on C + END DO ! end loop on R + + RETURN + END SUBROUTINE GETDEP_V diff --git a/src/model/src/vdiffacmx.F b/src/model/src/vdiffacmx.F index 3c12cfa7..1cdb48ef 100644 --- a/src/model/src/vdiffacmx.F +++ b/src/model/src/vdiffacmx.F @@ -56,6 +56,7 @@ SUBROUTINE VDIFFACMX ( DTSEC, SEDDY, DDEP, ICMP, DDEPJ, DDEPJ_FST, CNGRD ) USE VDIFF_MAP USE UTILIO_DEFN USE BIDI_MOD +C USE PT3D_EMIS_DEFN IMPLICIT NONE @@ -64,19 +65,19 @@ SUBROUTINE VDIFFACMX ( DTSEC, SEDDY, DDEP, ICMP, DDEPJ, DDEPJ_FST, CNGRD ) CHARACTER( 120 ) :: XMSG = ' ' C Arguments: - REAL, INTENT( IN ) :: DTSEC ! model time step in seconds + REAL*8, INTENT( IN ) :: DTSEC ! model time step in seconds C--- SEDDY is strictly an input, but it gets modified here - REAL, INTENT( INOUT ) :: SEDDY ( :,:,: ) ! flipped EDDYV - REAL, INTENT( INOUT ) :: DDEP ( :,:,: ) ! ddep accumulator - REAL, INTENT( INOUT ) :: ICMP ( :,:,: ) ! component flux accumlator - REAL, INTENT( INOUT ), OPTIONAL :: DDEPJ ( :,:,:,: ) ! ddep for mosaic - REAL, INTENT( INOUT ), OPTIONAL :: DDEPJ_FST( :,:,:,: ) ! ddep for stomtal/cuticular pathway + REAL*8, INTENT( INOUT ) :: SEDDY ( :,:,: ) ! flipped EDDYV + REAL*8, INTENT( INOUT ) :: DDEP ( :,:,: ) ! ddep accumulator + REAL*8, INTENT( INOUT ) :: ICMP ( :,:,: ) ! component flux accumlator + REAL*8, INTENT( INOUT ), OPTIONAL :: DDEPJ ( :,:,:,: ) ! ddep for mosaic + REAL*8, INTENT( INOUT ), OPTIONAL :: DDEPJ_FST( :,:,:,: ) ! ddep for stomtal/cuticular pathway REAL, INTENT( INOUT ) :: CNGRD ( :,:,:,: ) ! cgrid replacement C Parameters: C explicit, THETA = 0, implicit, THETA = 1 ! Crank-Nicholson: THETA = 0.5 - REAL, PARAMETER :: THETA = 0.5, + REAL*8, PARAMETER :: THETA = 0.5, & THBAR = 1.0 - THETA C External Functions: None @@ -87,26 +88,26 @@ SUBROUTINE VDIFFACMX ( DTSEC, SEDDY, DDEP, ICMP, DDEPJ, DDEPJ_FST, CNGRD ) LOGICAL, SAVE :: FIRSTIME = .TRUE. - REAL, ALLOCATABLE, SAVE :: DD_FAC ( : ) ! combined subexpression - REAL, ALLOCATABLE, SAVE :: DDBF ( : ) ! secondary DDEP - REAl, ALLOCATABLE, SAVE :: CMPF ( : ) ! intermediate CMP - REAL, ALLOCATABLE, SAVE :: CONC ( :,: ) ! secondary CGRID expression - REAL, ALLOCATABLE, SAVE :: EMIS ( :,: ) ! emissions subexpression - REAL DTDENS1 ! DT * layer 1 air density + REAL*8, ALLOCATABLE, SAVE :: DD_FAC ( : ) ! combined subexpression + REAL*8, ALLOCATABLE, SAVE :: DDBF ( : ) ! secondary DDEP + REAL*8, ALLOCATABLE, SAVE :: CMPF ( : ) ! intermediate CMP + REAL*8, ALLOCATABLE, SAVE :: CONC ( :,: ) ! secondary CGRID expression + REAL*8, ALLOCATABLE, SAVE :: EMIS ( :,: ) ! emissions subexpression + REAL*8 DTDENS1 ! DT * layer 1 air density C ACM Local Variables - REAL DFACP, DFACQ - REAL RP, RQ - REAL, ALLOCATABLE, SAVE :: DEPVCR ( : ) ! dep vel in one cell - REAL, ALLOCATABLE, SAVE :: EFAC1 ( : ) - REAL, ALLOCATABLE, SAVE :: EFAC2 ( : ) - REAL, ALLOCATABLE, SAVE :: POL ( : ) ! prodn/lossrate = PLDV/DEPV - REAL PLDV_HONO ! PLDV for HONO - REAL DEPV_NO2 ! dep vel of NO2 - REAL DEPV_HNO3 ! dep vel of HNO3 + REAL*8 DFACP, DFACQ + REAL*8 RP, RQ + REAL*8, ALLOCATABLE, SAVE :: DEPVCR ( : ) ! dep vel in one cell + REAL*8, ALLOCATABLE, SAVE :: EFAC1 ( : ) + REAL*8, ALLOCATABLE, SAVE :: EFAC2 ( : ) + REAL*8, ALLOCATABLE, SAVE :: POL ( : ) ! prodn/lossrate = PLDV/DEPV + REAL*8 PLDV_HONO ! PLDV for HONO + REAL*8 DEPV_NO2 ! dep vel of NO2 + REAL*8 DEPV_HNO3 ! dep vel of HNO3 INTEGER, SAVE :: NO2_HIT, HONO_HIT, HNO3_HIT, NO2_MAP, HNO3_MAP INTEGER, SAVE :: NH3_HIT - REAL DTS + REAL*8 DTS INTEGER, SAVE :: LOGDEV INTEGER ASTAT @@ -167,6 +168,7 @@ SUBROUTINE VDIFFACMX ( DTSEC, SEDDY, DDEP, ICMP, DDEPJ, DDEPJ_FST, CNGRD ) END IF ! if Firstime +C CALL GET_PT3D_EMIS () C ------------------------------------------- Row, Col LOOPS ----------- DTS = DTSEC @@ -220,7 +222,11 @@ SUBROUTINE VDIFFACMX ( DTSEC, SEDDY, DDEP, ICMP, DDEPJ, DDEPJ_FST, CNGRD ) IF ( V .EQ. HNO3_HIT ) THEN S = HNO3_MAP CONC( S,1 ) = POL( V ) + ( CONC( S,1 ) - POL( V ) ) * EFAC1( V ) - DEPV_HNO3 = DEPVCR( V ) + PLDV_HONO / CONC( NO2_MAP,1 ) + IF (CONC( NO2_MAP,1 ) .NE. 0) THEN + DEPV_HNO3 = DEPVCR( V ) + PLDV_HONO / CONC( NO2_MAP,1 ) + ELSE + DEPV_HNO3 = DEPVCR( V ) + END IF DD_FAC( V ) = DTDENS1 * DD_CONV( V ) * DEPV_HNO3 DDBF( V ) = DDBF( V ) + THETA * DD_FAC( V ) * CONC( S,1 ) @@ -231,7 +237,11 @@ SUBROUTINE VDIFFACMX ( DTSEC, SEDDY, DDEP, ICMP, DDEPJ, DDEPJ_FST, CNGRD ) C reduce the NO2 conc. in the atmosphere without affecting the depositional loss. ELSE IF ( V .EQ. NO2_HIT ) THEN S = NO2_MAP - DEPV_NO2 = DEPVCR( V ) + 2.0 * PLDV_HONO / CONC( S,1 ) + IF (CONC( S,1 ) .NE. 0) THEN + DEPV_NO2 = DEPVCR( V ) + 2.0 * PLDV_HONO / CONC( S,1 ) + ELSE + DEPV_NO2 = DEPVCR( V ) + END IF EFAC1 ( V ) = EXP( -DEPV_NO2 * RP ) EFAC2 ( V ) = EXP( -DEPV_NO2 * RQ ) POL ( V ) = PLDV( V,C,R ) / DEPV_NO2 diff --git a/src/shr/aqm_const_mod.F90 b/src/shr/aqm_const_mod.F90 index 80099ae4..979689a7 100644 --- a/src/shr/aqm_const_mod.F90 +++ b/src/shr/aqm_const_mod.F90 @@ -10,6 +10,10 @@ module aqm_const_mod real(AQM_KIND_R8), parameter :: one = 1.0_AQM_KIND_R8 real(AQM_KIND_R8), parameter :: half = 0.5_AQM_KIND_R8 + real(AQM_KIND_R8), parameter :: aqm_pi = 3.14159265358979323846_AQM_KIND_R8 + real(AQM_KIND_R8), parameter :: deg_to_rad = aqm_pi / 180._AQM_KIND_R8 + real(AQM_KIND_R8), parameter :: rad_to_deg = 180._AQM_KIND_R8 / aqm_pi + ! -- include CMAQ constants include SUBST_CONST diff --git a/src/shr/aqm_emis_mod.F90 b/src/shr/aqm_emis_mod.F90 index 17bc086c..e4325770 100644 --- a/src/shr/aqm_emis_mod.F90 +++ b/src/shr/aqm_emis_mod.F90 @@ -7,7 +7,8 @@ module aqm_emis_mod use aqm_types_mod use aqm_tools_mod use aqm_internal_mod - use aqm_model_mod, only : aqm_model_get, aqm_state_type + use aqm_model_mod, only : aqm_model_get, aqm_model_domain_get, aqm_state_type + use aqm_const_mod, only : deg_to_rad, rad_to_deg implicit none @@ -196,12 +197,16 @@ subroutine aqm_emis_src_create(config, em, rc) em(item) % type = "" em(item) % path = "" em(item) % file = "" - em(item) % frequency = "" - em(item) % format = "netcdf" - em(item) % iomode = "read" - em(item) % iofmt = AQMIO_FMT_NETCDF - em(item) % irec = 0 + em(item) % frequency = "" + em(item) % format = "netcdf" + em(item) % iomode = "read" + em(item) % iofmt = AQMIO_FMT_NETCDF + em(item) % irec = 0 + em(item) % count = 0 + em(item) % layers = 1 em(item) % scalefactor = 1.0 + em(item) % topfraction = -1.0 + em(item) % gridded = .true. em(item) % sync = .false. em(item) % verbose = .false. em(item) % logprefix = "" @@ -209,12 +214,28 @@ subroutine aqm_emis_src_create(config, em, rc) em(item) % plumerise = "" em(item) % specfile = "" em(item) % specprofile = "" + em(item) % latname = "" + em(item) % lonname = "" + em(item) % stkdmname = "" + em(item) % stkhtname = "" + em(item) % stktkname = "" + em(item) % stkvename = "" nullify(em(item) % species) nullify(em(item) % sources) nullify(em(item) % units) nullify(em(item) % factors) nullify(em(item) % fields) + nullify(em(item) % lat) + nullify(em(item) % lon) + nullify(em(item) % stkdm) + nullify(em(item) % stkht) + nullify(em(item) % stktk) + nullify(em(item) % stkve) + nullify(em(item) % rates) nullify(em(item) % dens_flag) + nullify(em(item) % ip) + nullify(em(item) % jp) + nullify(em(item) % ijmap) nullify(em(item) % table) call ESMF_ConfigGetAttribute(config, emName, rc=localrc) @@ -586,6 +607,178 @@ subroutine aqm_emis_src_init(model, em, rc) rcToReturn=rc)) & return ! bail out end if + if (trim(em % plumerise) /= "none") then + call ESMF_ConfigGetAttribute(config, em % topfraction, & + label=trim(em % name)//"_plume_top_fraction:", default=-1.0, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + if (em % verbose) then + write(msgString,'(g20.8)') em % topfraction + call ESMF_LogWrite(trim(em % logprefix)//": "//pName & + //": plume_top_fraction: "//adjustl(msgString), ESMF_LOGMSG_INFO, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + end if + if (em % topfraction > 1.0) then + call ESMF_LogSetError(ESMF_RC_NOT_VALID, & + msg="plume top fraction must not exceed 1.0", & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc) + return ! bail out + end if + end if + case ("point-source") + em % gridded = .false. + ! -- get plumerise type + call ESMF_ConfigGetAttribute(config, value, & + label=trim(em % name)//"_plume_rise:", default="default", rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + em % plumerise = ESMF_UtilStringLowerCase(value, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + if (em % verbose) then + call ESMF_LogWrite(trim(em % logprefix)//": "//pName & + //": plume_rise: "//trim(em % plumerise), ESMF_LOGMSG_INFO, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + end if + ! -- get emission layers + call ESMF_ConfigGetAttribute(config, em % layers, & + label=trim(em % name)//"_layers:", default=1, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + if (em % layers <= 0) then + call ESMF_LogSetError(ESMF_RC_NOT_VALID, & + msg="- layers must be greater than 0", & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc) + return ! bail out + end if + if (em % verbose) then + write(msgString, '(a,": ",a,": layers: ",i0)') trim(em % logprefix), & + trim(pName), em % layers + call ESMF_LogWrite(trim(msgString), ESMF_LOGMSG_INFO, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + end if + ! -- get coordinate labels + call ESMF_ConfigFindLabel(config, trim(em % name) // "_latlon_names:", rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + call ESMF_ConfigGetAttribute(config, em % latname, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + call ESMF_ConfigGetAttribute(config, em % lonname, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + if (em % verbose) then + call ESMF_LogWrite(trim(em % logprefix)//": "//pName & + //": latlon_names: "//trim(em % latname)//" "//trim(em % lonname), & + ESMF_LOGMSG_INFO, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + end if + ! -- get stack parameters + call ESMF_ConfigGetAttribute(config, em % stkdmname, & + label=trim(em % name)//"_stack_diameter:", default='STKDM', rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + if (em % verbose) then + call ESMF_LogWrite(trim(em % logprefix)//": "//pName & + //": stack_diameter: "//trim(em % stkdmname), ESMF_LOGMSG_INFO, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + end if + call ESMF_ConfigGetAttribute(config, em % stkhtname, & + label=trim(em % name)//"_stack_height:", default='STKHT', rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + if (em % verbose) then + call ESMF_LogWrite(trim(em % logprefix)//": "//pName & + //": stack_height: "//trim(em % stkhtname), ESMF_LOGMSG_INFO, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + end if + call ESMF_ConfigGetAttribute(config, em % stktkname, & + label=trim(em % name)//"_stack_temperature:", default='STKTK', rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + if (em % verbose) then + call ESMF_LogWrite(trim(em % logprefix)//": "//pName & + //": stack_temperature: "//trim(em % stktkname), ESMF_LOGMSG_INFO, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + end if + call ESMF_ConfigGetAttribute(config, em % stkvename, & + label=trim(em % name)//"_stack_velocity:", default='STKVE', rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + if (em % verbose) then + call ESMF_LogWrite(trim(em % logprefix)//": "//pName & + //": stack_velocity: "//trim(em % stkvename), ESMF_LOGMSG_INFO, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + end if case ("product") em % iomode = "create" readFactors = .false. @@ -774,13 +967,12 @@ subroutine aqm_emis_src_init(model, em, rc) end if end if - ! -- create fields + ! -- create arrays allocate(em % species(fieldCount), & em % factors(fieldCount), & em % sources(fieldCount), & em % units(fieldCount), & - em % dens_flag(fieldCount), & - em % fields(fieldCount), stat=stat) + em % dens_flag(fieldCount), stat=stat) if (ESMF_LogFoundAllocError(statusToCheck=stat, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & file=__FILE__, & @@ -804,13 +996,6 @@ subroutine aqm_emis_src_init(model, em, rc) em % factors(fieldCount) = tmpFactorList(item) em % sources(fieldCount) = tmpSourceList(item) em % units(fieldCount) = tmpUnitsList(item) - em % fields(fieldCount) = ESMF_FieldCreate(grid, ESMF_TYPEKIND_R4, & - name=em % sources(fieldCount), rc=localrc) - if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & - line=__LINE__, & - file=__FILE__, & - rcToReturn=rc)) & - return ! bail out end if end do @@ -821,6 +1006,39 @@ subroutine aqm_emis_src_init(model, em, rc) rcToReturn=rc)) & return + ! -- create fields if needed + if (em % gridded) then + allocate(em % fields(fieldCount), stat=stat) + if (ESMF_LogFoundAllocError(statusToCheck=stat, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return + do item = 1, fieldCount + em % fields(item) = ESMF_FieldCreate(grid, ESMF_TYPEKIND_R4, & + name=em % sources(item), rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + end do + else + allocate(em % rates(fieldCount), stat=stat) + if (ESMF_LogFoundAllocError(statusToCheck=stat, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return + ! -- initialize ungridded emissions (point sources) + call aqm_emis_pts_init(model, em, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + end if + end if ! -- build unique list of species @@ -913,6 +1131,15 @@ subroutine aqm_emis_finalize(model, rc) em => is % wrap % emis(item) + if (associated(em % sources)) then + deallocate(em % sources, stat=stat) + if (ESMF_LogFoundDeallocError(statusToCheck=stat, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + nullify(em % sources) + end if if (associated(em % species)) then deallocate(em % species, stat=stat) if (ESMF_LogFoundDeallocError(statusToCheck=stat, msg=ESMF_LOGERR_PASSTHRU, & @@ -983,6 +1210,107 @@ subroutine aqm_emis_finalize(model, rc) return ! bail out nullify(em % fields) end if + if (associated(em % lat)) then + deallocate(em % lat, stat=stat) + if (ESMF_LogFoundDeallocError(statusToCheck=stat, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + nullify(em % lat) + end if + if (associated(em % lon)) then + deallocate(em % lon, stat=stat) + if (ESMF_LogFoundDeallocError(statusToCheck=stat, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + nullify(em % lon) + end if + if (associated(em % stkdm)) then + deallocate(em % stkdm, stat=stat) + if (ESMF_LogFoundDeallocError(statusToCheck=stat, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + nullify(em % stkdm) + end if + if (associated(em % stkht)) then + deallocate(em % stkht, stat=stat) + if (ESMF_LogFoundDeallocError(statusToCheck=stat, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + nullify(em % stkht) + end if + if (associated(em % stktk)) then + deallocate(em % stktk, stat=stat) + if (ESMF_LogFoundDeallocError(statusToCheck=stat, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + nullify(em % stktk) + end if + if (associated(em % stkve)) then + deallocate(em % stkve, stat=stat) + if (ESMF_LogFoundDeallocError(statusToCheck=stat, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + nullify(em % stkve) + end if + if (associated(em % ip)) then + deallocate(em % ip, stat=stat) + if (ESMF_LogFoundDeallocError(statusToCheck=stat, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + nullify(em % ip) + end if + if (associated(em % jp)) then + deallocate(em % jp, stat=stat) + if (ESMF_LogFoundDeallocError(statusToCheck=stat, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + nullify(em % jp) + end if + if (associated(em % ijmap)) then + deallocate(em % ijmap, stat=stat) + if (ESMF_LogFoundDeallocError(statusToCheck=stat, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + nullify(em % ijmap) + end if + if (associated(em % rates)) then + do n = 1, size(em % rates) + if (associated(em % rates(n) % values)) then + deallocate(em % rates(n) % values, stat=stat) + if (ESMF_LogFoundDeallocError(statusToCheck=stat, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return + nullify(em % rates(n) % values) + end if + end do + deallocate(em % rates, stat=stat) + if (ESMF_LogFoundDeallocError(statusToCheck=stat, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + nullify(em % rates) + end if call ESMF_AlarmDestroy(em % alarm, rc=localrc) if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & @@ -1028,6 +1356,7 @@ subroutine aqm_emis_finalize(model, rc) end subroutine aqm_emis_finalize + subroutine aqm_emis_update(model, rc) type(ESMF_GridComp) :: model integer, optional, intent(out) :: rc @@ -1115,12 +1444,23 @@ subroutine aqm_emis_update(model, rc) end do else em % irec = em % irec + 1 - call AQMIO_Read(em % IO, em % fields, timeSlice=em % irec, rc=localrc) - if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & - line=__LINE__, & - file=__FILE__, & - rcToReturn=rc)) & - return ! bail out + if (em % gridded) then + call AQMIO_Read(em % IO, em % fields, timeSlice=em % irec, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + else + do n = 1, size(em % sources) + call AQMIO_DataRead(em % IO, em % rates(n) % values, em % sources(n), timeSlice=em % irec, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + end do + end if end if call ESMF_AlarmRingerOff(em % alarm, rc=localrc) if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & @@ -1134,6 +1474,7 @@ subroutine aqm_emis_update(model, rc) end subroutine aqm_emis_update + function aqm_emis_get(etype) result (ep) character(len=*), intent(in) :: etype type(aqm_internal_emis_type), pointer :: ep @@ -1153,6 +1494,7 @@ function aqm_emis_get(etype) result (ep) end function aqm_emis_get + function aqm_emis_data_get() result (ep) type(aqm_internal_emis_type), pointer :: ep(:) @@ -1162,6 +1504,7 @@ function aqm_emis_data_get() result (ep) end function aqm_emis_data_get + logical function aqm_emis_ispresent(etype) character(len=*), intent(in) :: etype @@ -1169,38 +1512,55 @@ logical function aqm_emis_ispresent(etype) end function aqm_emis_ispresent + subroutine aqm_emis_desc( etype, nlays, nvars, vnames, units ) - character(len=*), intent(in) :: etype - integer, intent(out) :: nlays - integer, intent(out) :: nvars - character(len=16), intent(out) :: vnames(:) - character(len=16), intent(out) :: units(:) + character(len=*), intent(in) :: etype + integer, optional, intent(out) :: nlays + integer, optional, intent(out) :: nvars + character(len=16), optional, intent(out) :: vnames(:) + character(len=16), optional, intent(out) :: units(:) ! -- local variables - integer :: item + integer :: localrc + integer :: item, nsrc type(aqm_internal_emis_type), pointer :: em ! -- begin - nlays = 0 - nvars = 0 - vnames = "" - units = "" - ! -- get emission data nullify(em) em => aqm_emis_get(etype) if (associated(em)) then - nlays = 1 - nvars = size( em % table, dim=1 ) - vnames( 1:nvars ) = em % table( 1:nvars, 1 ) - units ( 1:nvars ) = em % table( 1:nvars, 2 ) + if (present(nlays)) then + ! -- get number of vertical levels from model + call aqm_model_domain_get(nl=nlays, rc=localrc) + if (aqm_rc_check(localrc, msg="Failure to retrieve model layers", & + file=__FILE__, line=__LINE__)) return + select case (trim(em % plumerise)) + case ("briggs","default") + nlays = min(nlays, em % layers) + case ("sofiev") + ! -- use full vertical column for fire emissions + case default + nlays = 1 + end select + end if + nsrc = size( em % table, dim=1 ) + if (present(nvars)) nvars = nsrc + if (present(vnames)) vnames( 1:nsrc ) = em % table( 1:nsrc, 1 ) + if (present(units)) units ( 1:nsrc ) = em % table( 1:nsrc, 2 ) + else + if (present(nlays)) nlays = 0 + if (present(nvars)) nvars = 0 + if (present(vnames)) vnames = "" + if (present(units)) units = "" end if end subroutine aqm_emis_desc - subroutine aqm_emis_read(etype, spcname, buffer, localDe, rc) - character(len=*), intent(in) :: etype + + subroutine aqm_emis_grd_read(em, spcname, buffer, localDe, rc) + type(aqm_internal_emis_type) :: em character(len=*), intent(in) :: spcname real, intent(inout) :: buffer(*) integer, optional, intent(in) :: localDe @@ -1213,7 +1573,6 @@ subroutine aqm_emis_read(etype, spcname, buffer, localDe, rc) character(len=ESMF_MAXSTR) :: msgString real(ESMF_KIND_R4), pointer :: fptr(:,:) type(aqm_state_type), pointer :: stateIn - type(aqm_internal_emis_type), pointer :: em ! -- begin if (present(rc)) rc = AQM_RC_SUCCESS @@ -1221,9 +1580,6 @@ subroutine aqm_emis_read(etype, spcname, buffer, localDe, rc) ! -- NOTE: input emissions need to be converted to surface densities here ! -- since grid cell area is set to 1 internally. - em => aqm_emis_get(etype) - ! -- bail out if no emissions available - if (.not.associated(em)) return ! -- bail out if this is a product if (trim(em % iomode) /= "read") return @@ -1244,12 +1600,6 @@ subroutine aqm_emis_read(etype, spcname, buffer, localDe, rc) return ! bail out end if - if (trim(em % type) == "fengsha") then - ! -- ensure fengsha input variables are not normalized by area like - ! -- emissions conversions below - em % dens_flag(item) = 1 - end if - select case (em % dens_flag(item)) case (:-1) ! -- this case indicates that input emissions are provided as totals/cell @@ -1299,8 +1649,108 @@ subroutine aqm_emis_read(etype, spcname, buffer, localDe, rc) return end select if (em % verbose) then - write(msgString, '(a12,": read",12x,": ",a16,": min/max = ",2g20.8)') & - em % logprefix, spcname, minval(fptr), maxval(fptr) + write(msgString, '(a,": read: ",a16,": min/max = ",2g20.8)') & + trim(em % logprefix), spcname, minval(fptr), maxval(fptr) + call ESMF_LogWrite(trim(msgString), ESMF_LOGMSG_INFO, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + end if + end if + end do + + end subroutine aqm_emis_grd_read + + + subroutine aqm_emis_pts_read(em, spcname, buffer, ip, jp, ijmap, localDe, rc) + type(aqm_internal_emis_type) :: em + character(len=*), intent(in) :: spcname + real, intent(inout) :: buffer(*) + integer, optional, pointer :: ip(:) + integer, optional, pointer :: jp(:) + integer, optional, pointer :: ijmap(:) + integer, optional, intent(in) :: localDe + integer, optional, intent(out) :: rc + + ! -- local variables + integer :: localrc + integer :: item, i, j, m, n + character(len=ESMF_MAXSTR) :: msgString + real(ESMF_KIND_R4) :: em_min, em_max + type(aqm_state_type), pointer :: stateIn + + ! -- begin + if (present(rc)) rc = AQM_RC_SUCCESS + + ! -- NOTE: input emissions need to be converted to surface densities here + ! -- since grid cell area is set to 1 internally. + + ! -- bail out if this is a product + if (trim(em % iomode) /= "read") return + + if (em % count == 0) return + + call aqm_model_get(stateIn=stateIn, rc=localrc) + if (aqm_rc_check(localrc, msg="Failure to retrieve model input state", & + file=__FILE__, line=__LINE__, rc=rc)) return + + ! -- return grid locations and map + if (present(ip)) ip => em % ip + if (present(jp)) jp => em % jp + if (present(ijmap)) ijmap => em % ijmap + + do item = 1, size(em % species) + if (trim(spcname) == trim(em % species(item))) then + + select case (em % dens_flag(item)) + case (:-1) + ! -- this case indicates that input emissions are provided as totals/cell + ! -- while surface densities are required and should never occur in CMAQ + do m = 1, size(em % ijmap) + n = em % ijmap(m) + i = em % ip(n) + j = em % jp(n) + buffer(n) = buffer(n) & + + em % factors(item) * em % rates(item) % values(n) / stateIn % area(i,j) & + / stateIn % area(i,j) + end do + case (0) + ! -- emissions are totals over each grid cell + do m = 1, size(em % ijmap) + n = em % ijmap(m) + i = em % ip(n) + j = em % jp(n) + buffer(n) = buffer(n) & + + em % factors(item) * em % rates(item) % values(n) / stateIn % area(i,j) + end do + case (1:) + ! -- emissions are already provided as surface densities, no need to normalize + do m = 1, size(em % ijmap) + n = em % ijmap(m) + buffer(n) = buffer(n) & + + em % factors(item) * em % rates(item) % values(n) + end do + case default + ! -- this case should never occur + call ESMF_LogSetError(ESMF_RC_INTNRL_BAD, & + msg="uncategorized emission source", & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc) + return + end select + if (em % verbose) then + em_min = huge(0._ESMF_KIND_R4) + em_max = -em_min + do m = 1, size(em % ijmap) + n = em % ijmap(m) + em_min = min( em_min, em % rates(item) % values(n) ) + em_max = max( em_max, em % rates(item) % values(n) ) + end do + write(msgString, '(a,": read: ",a16,": min/max = ",2g20.8)') & + trim(em % logprefix), spcname, em_min, em_max call ESMF_LogWrite(trim(msgString), ESMF_LOGMSG_INFO, rc=localrc) if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & line=__LINE__, & @@ -1311,6 +1761,493 @@ subroutine aqm_emis_read(etype, spcname, buffer, localDe, rc) end if end do + end subroutine aqm_emis_pts_read + + + subroutine aqm_emis_read(etype, spcname, buffer, ip, jp, ijmap, localDe, rc) + character(len=*), intent(in) :: etype + character(len=*), intent(in) :: spcname + real, intent(inout) :: buffer(*) + integer, optional, pointer :: ip(:) + integer, optional, pointer :: jp(:) + integer, optional, pointer :: ijmap(:) + integer, optional, intent(in) :: localDe + integer, optional, intent(out) :: rc + + ! -- local variables + type(aqm_internal_emis_type), pointer :: em + + ! -- begin + if (present(rc)) rc = AQM_RC_SUCCESS + + ! -- retrieve emissions + em => aqm_emis_get(etype) + + ! -- bail out if no emissions available + if (.not.associated(em)) return + + ! -- bail out if this is a product + if (trim(em % iomode) /= "read") return + + ! -- select emission reader + if (em % gridded) then + call aqm_emis_grd_read(em, spcname, buffer, localDe=localDe, rc=rc) + else + call aqm_emis_pts_read(em, spcname, buffer, ip=ip, jp=jp, ijmap=ijmap, localDe=localDe, rc=rc) + end if + end subroutine aqm_emis_read - + + + subroutine aqm_emis_pts_map(model, em, rc) +! --------------------------------------------------------------------------------- +! +! This subroutine is based on source code included in NASA's MAPL library: +! +! https://github.com/GEOS-ESM/MAPL +! +! MAPL is licensed under the Apache license version 2.0. See notice below: +! +! --------------------------------------------------------------------------------- +! +! NASA Docket No. GSC-15,354-1, and identified as "GEOS-5 GCM Modeling Software +! +! Copyright (C) 2008 United States Government as represented by the Administrator +! of the National Aeronautics and Space Administration. All Rights Reserved. +! +! Licensed under the Apache License, Version 2.0 (the "License"); you may not use +! this file except in compliance with the License. You may obtain a copy of the +! License at +! +! http://www.apache.org/licenses/LICENSE-2.0 +! +! Unless required by applicable law or agreed to in writing, software distributed +! under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR +! CONDITIONS OF ANY KIND, either express or implied. See the License for the +! specific language governing permissions and limitations under the License. +! +! --------------------------------------------------------------------------------- + + type(ESMF_GridComp) :: model + type(aqm_internal_emis_type) :: em + integer, optional, intent(out) :: rc + + ! -- local variables + integer :: localrc, stat + integer :: localDeCount + integer :: ic, jc, ijcount + integer :: m, n + integer, dimension(2) :: ec, tc + + real(AQM_KIND_R8), allocatable :: lon(:,:) + real(AQM_KIND_R8), allocatable :: lat(:,:) + real(AQM_KIND_R8), allocatable :: lonc(:,:) + real(AQM_KIND_R8), allocatable :: latc(:,:) + real(AQM_KIND_R8), allocatable :: lonp(:) + real(AQM_KIND_R8), allocatable :: latp(:) + + real(ESMF_KIND_R8), pointer :: fp(:,:) + real(ESMF_KIND_R8), pointer :: coord(:,:) + + type(ESMF_Grid) :: grid + type(ESMF_Field) :: field + type(ESMF_RouteHandle) :: rh + type(ESMF_CoordSys_Flag) :: coordSys + type(ESMF_Index_Flag) :: indexflag + + real(AQM_KIND_R8) :: emin, emax, fmin, fmax + character(len=ESMF_MAXSTR) :: msg + + character(len=*), parameter :: pName = "init: source" + + ! -- begin + if (present(rc)) rc = ESMF_SUCCESS + + nullify(fp) + + ! -- bail out if this is not point source + if (trim(em % type) /= "point-source") return + + ! -- bail out if this is a product + if (trim(em % iomode) /= "read") return + + ! -- get model grid + call ESMF_GridCompGet(model, grid=grid, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + + ! -- create work field to compute grid coordinate halos + field = ESMF_FieldCreate(grid, ESMF_TYPEKIND_R8, staggerLoc=ESMF_STAGGERLOC_CORNER, & + totalLWidth=[1,1], totalUWidth=[1,1], rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + + call ESMF_FieldHaloStore(field, rh, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + + ! -- compute halos only on active PETs + call ESMF_GridGet(grid, localDeCount=localDeCount, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + + if (localDeCount > 0) then + + ! -- get coordinate system + call ESMF_GridGet(grid, coordSys=coordSys, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + + ! -- (1) center coordinates + call ESMF_GridGet(grid, staggerloc=ESMF_STAGGERLOC_CENTER, & + localDe=0, exclusiveCount=ec, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + + ic = ec(1) + jc = ec(2) + + call ESMF_GridGetCoord(grid, coordDim=1, staggerloc=ESMF_STAGGERLOC_CENTER, & + fArrayPtr=coord, totalCount=tc, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + + allocate(lon(ic, jc), lat(ic, jc), lonc(ic+1, jc+1), latc(ic+1, jc+1), stat=stat) + if (ESMF_LogFoundAllocError(statusToCheck=stat, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return + + lon = coord + + call ESMF_GridGetCoord(grid, coordDim=2, staggerloc=ESMF_STAGGERLOC_CENTER, & + fArrayPtr=coord, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + + lat = coord + + ! -- (2) corner coordinates + call ESMF_GridGetCoord(grid, coordDim=1, staggerloc=ESMF_STAGGERLOC_CORNER, & + fArrayPtr=coord, totalCount=tc, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + + call ESMF_FieldGet(field, farrayPtr=fp, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + + fp(1:tc(1),1:tc(2)) = coord + + end if + + call ESMF_FieldHalo(field, rh, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + + if (localDeCount > 0) then + + lonc = fp(1:ec(1)+1,1:ec(2)+1) + + call ESMF_GridGetCoord(grid, coordDim=2, staggerloc=ESMF_STAGGERLOC_CORNER, & + fArrayPtr=coord, totalCount=tc, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + + fp(1:tc(1),1:tc(2)) = coord + + end if + + call ESMF_FieldHalo(field, rh, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + + if (localDeCount > 0) latc = fp(1:ec(1)+1,1:ec(2)+1) + + call ESMF_FieldDestroy(field, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + + call ESMF_FieldHaloRelease(rh, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + + if (localDeCount < 1) return + + if (em % count == 0) return + + if (coordSys == ESMF_COORDSYS_SPH_DEG) then + lon = lon * deg_to_rad + lat = lat * deg_to_rad + lonc = lonc * deg_to_rad + latc = latc * deg_to_rad + else if (coordSys == ESMF_COORDSYS_SPH_RAD) then + ! -- nothing to do + else + ! -- unsupported coordinate system + call ESMF_LogSetError(ESMF_RC_NOT_VALID, & + msg="coordinate system unsupported by emission mapping method.", & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc) + return + end if + + allocate(lonp(em % count), latp(em % count), stat=stat) + if (ESMF_LogFoundAllocError(statusToCheck=stat, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return + + ! -- emission coordinates are expected in degrees + lonp = em % lon * deg_to_rad + latp = em % lat * deg_to_rad + + allocate(em % ip(em % count), em % jp(em % count), stat=stat) + if (ESMF_LogFoundAllocError(statusToCheck=stat, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return + + ! -- map point-source locations to grid + call aqm_gridloc_get(lon, lat, lonc, latc, lonp, latp, em % ip, em % jp, ijcount) + + if (ijcount > 0) then + allocate(em % ijmap(ijcount), stat=stat) + if (ESMF_LogFoundAllocError(statusToCheck=stat, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return + m = 0 + do n = 1, em % count + if (em % ip(n) > 0 .and. em % jp(n) > 0) then + m = m + 1 + em % ijmap(m) = n + end if + end do + em % count = m + else + ! -- no local sources, disable emissions + em % count = 0 + end if + + if (em % verbose) then + + ! -- log coordinates + lon = lon * rad_to_deg + lat = lat * rad_to_deg + lonc = lonc * rad_to_deg + latc = latc * rad_to_deg + lonp = lonp * rad_to_deg + latp = latp * rad_to_deg + + if (coordSys == ESMF_COORDSYS_SPH_DEG) then + call ESMF_LogWrite(trim(em % logprefix) // ": " // trim(pName) // & + ": grid coord. in degrees", ESMF_LOGMSG_INFO, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + else if (coordSys == ESMF_COORDSYS_SPH_RAD) then + call ESMF_LogWrite(trim(em % logprefix) // ": " // trim(pName) // & + ": grid coord. in radians", ESMF_LOGMSG_INFO, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + end if + + where(lon < 0.) lon = lon + 360. + where(lonc < 0.) lonc = lonc + 360. + where(lonp < 0.) lonp = lonp + 360. + + write( msg, '(": mapping: centers: min/max: ",4g20.5)') & + minval(lon), maxval(lon), minval(lat), maxval(lat) + call ESMF_LogWrite(trim(em % logprefix) // ": " // trim(pName) // & + msg, ESMF_LOGMSG_INFO, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + write( msg, '(": mapping: corners: min/max: ",4g20.5)') & + minval(lonc), maxval(lonc), minval(latc), maxval(latc) + call ESMF_LogWrite(trim(em % logprefix) // ": " // trim(pName) // & + msg, ESMF_LOGMSG_INFO, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + + emin = huge(0.0) + fmin = huge(0.0) + emax = -emin + fmax = -fmin + do n = 1, em % count + m = em % ijmap(n) + emin = min(emin,latp(m)) + emax = max(emax,latp(m)) + fmin = min(fmin,lonp(m)) + fmax = max(fmax,lonp(m)) + end do + write( msg, '(": mapping: pt-srcs: min/max: ",4g20.5)') fmin, fmax, emin, emax + call ESMF_LogWrite(trim(em % logprefix) // ": " // trim(pName) // & + msg, ESMF_LOGMSG_INFO, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + + end if + + ! -- free up memory + deallocate(lat, latc, latp, lon, lonc, lonp, stat=stat) + if (ESMF_LogFoundDeallocError(statusToCheck=stat, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return + + end subroutine aqm_emis_pts_map + + + subroutine aqm_emis_pts_init(model, em, rc) + type(ESMF_GridComp) :: model + type(aqm_internal_emis_type) :: em + integer, optional, intent(out) :: rc + + ! -- local variables + integer :: localrc, item + character(len=ESMF_MAXSTR) :: msgString + + character(len=*), parameter :: pName = "init: source" + + ! -- begin + if (present(rc)) rc = ESMF_SUCCESS + + ! -- bail out if gridded emissions + if (em % gridded) return + + ! -- initialize pointers & counters + do item = 1, size(em % sources) + nullify(em % rates(item) % values) + end do + + ! -- read stack locations + call AQMIO_DataRead(em % IO, em % lat, em % latname, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + call AQMIO_DataRead(em % IO, em % lon, em % lonname, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + + ! -- read stack parameters + call AQMIO_DataRead(em % IO, em % stkdm, em % stkdmname, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + call AQMIO_DataRead(em % IO, em % stkht, em % stkhtname, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + call AQMIO_DataRead(em % IO, em % stktk, em % stktkname, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + call AQMIO_DataRead(em % IO, em % stkve, em % stkvename, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + + ! -- set initial (maximum) count + em % count = size(em % lat) + + ! -- map point source emissions to grid + call aqm_emis_pts_map(model, em, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + + if (em % verbose) then + write(msgString, '("mapped to ",i0," local grid points")') em % count + call ESMF_LogWrite(trim(em % logprefix)//": "//pName & + //": "//trim(msgString), ESMF_LOGMSG_INFO, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + end if + + end subroutine aqm_emis_pts_init + end module aqm_emis_mod diff --git a/src/shr/aqm_fires_mod.F90 b/src/shr/aqm_fires_mod.F90 index 785950fd..405d74c7 100644 --- a/src/shr/aqm_fires_mod.F90 +++ b/src/shr/aqm_fires_mod.F90 @@ -25,6 +25,7 @@ subroutine aqm_plume_sofiev(em, frp, profile, rc) integer :: c, r, l integer :: lev0, lev1 real :: Hp, pblh, th0, th1, dz + real :: tfrac, w real(AQM_KIND_R8) :: hbl real(AQM_KIND_R8), pointer :: phi(:) type(aqm_state_type), pointer :: state @@ -44,16 +45,23 @@ subroutine aqm_plume_sofiev(em, frp, profile, rc) ! -- get model info call aqm_model_get(stateIn=state, rc=localrc) if (aqm_rc_check(localrc, msg="Failed to retrieve model state", & - file=__FILE__, line=__LINE__)) return + file=__FILE__, line=__LINE__, rc=rc)) return ! -- get domain info call aqm_model_domain_get(ids=is, ide=ie, jds=js, jde=je, nl=nl, rc=localrc) if (aqm_rc_check(localrc, msg="Failed to retrieve grid coordinates", & - file=__FILE__, line=__LINE__)) return + file=__FILE__, line=__LINE__, rc=rc)) return nx = ie - is + 1 ny = je - js + 1 + ! -- compute layer empirical weights + if (aqm_rc_test((em % topfraction > 1.0), & + msg="Plume top fraction must be between 0.0 and 1.0", & + file=__FILE__, line=__LINE__, rc=rc)) return + + w = min( 1.0 - em % topfraction, 1.0 ) + ! -- select free-troposphere vertical level k = 0 do r = 1, ny @@ -64,11 +72,11 @@ subroutine aqm_plume_sofiev(em, frp, profile, rc) lev0 = minloc(phi, 1, mask=phi >= hbl) if (aqm_rc_test((phi(lev0) < hbl), & msg="Could not find first free-troposphere layer", & - file=__FILE__, line=__LINE__)) return + file=__FILE__, line=__LINE__,rc=rc)) return lev1 = lev0 + 1 if (aqm_rc_test((lev1 > nl), & msg="Not enough vertical levels", & - file=__FILE__, line=__LINE__)) return + file=__FILE__, line=__LINE__,rc=rc)) return dz = onebg * ( phi(lev1) - phi(lev0) ) th0 = state % temp(c,r,lev0) * (p_ref / state % prl(c,r,lev0)) ** rcp @@ -84,10 +92,13 @@ subroutine aqm_plume_sofiev(em, frp, profile, rc) if (lev1 > lev0) then dz = 0.0 - do l = lev0, lev1 - profile(c,r,l) = ( phi(l) - dz ) / phi(lev1) + tfrac = 0.0 + do l = lev0, lev1-1 + profile(c,r,l) = w * ( phi(l) - dz ) / phi(lev1) dz = phi(l) + tfrac = tfrac + profile(c,r,l) end do + profile(c,r,lev1) = 1.0 - tfrac else profile(c,r,lev0) = 1.0 end if diff --git a/src/shr/aqm_internal_mod.F90 b/src/shr/aqm_internal_mod.F90 index 82465396..3b4a1108 100644 --- a/src/shr/aqm_internal_mod.F90 +++ b/src/shr/aqm_internal_mod.F90 @@ -4,6 +4,10 @@ module aqm_internal_mod implicit none + type aqm_internal_rate_type + real(ESMF_KIND_R4), dimension(:), pointer :: values => null() + end type aqm_internal_rate_type + type aqm_internal_emis_type character(len=ESMF_MAXSTR) :: name character(len=ESMF_MAXSTR) :: type @@ -15,10 +19,20 @@ module aqm_internal_mod character(len=ESMF_MAXSTR) :: plumerise character(len=ESMF_MAXSTR) :: specfile character(len=ESMF_MAXSTR) :: specprofile + character(len=ESMF_MAXSTR) :: latname + character(len=ESMF_MAXSTR) :: lonname + character(len=ESMF_MAXSTR) :: stkdmname + character(len=ESMF_MAXSTR) :: stkhtname + character(len=ESMF_MAXSTR) :: stktkname + character(len=ESMF_MAXSTR) :: stkvename character(len=6) :: period + logical :: gridded logical :: sync logical :: verbose real :: scalefactor + real :: topfraction + integer(ESMF_KIND_I4) :: layers + integer :: count integer :: irec integer :: iofmt character(len=ESMF_MAXSTR) :: iomode @@ -28,8 +42,18 @@ module aqm_internal_mod character(len=ESMF_MAXSTR), dimension(:), pointer :: species => null() character(len=ESMF_MAXSTR), dimension(:), pointer :: units => null() integer, dimension(:), pointer :: dens_flag => null() + integer, dimension(:), pointer :: ip => null() + integer, dimension(:), pointer :: jp => null() + integer, dimension(:), pointer :: ijmap => null() + real(ESMF_KIND_R4), dimension(:), pointer :: lat => null() + real(ESMF_KIND_R4), dimension(:), pointer :: lon => null() + real(ESMF_KIND_R4), dimension(:), pointer :: stkdm => null() + real(ESMF_KIND_R4), dimension(:), pointer :: stkht => null() + real(ESMF_KIND_R4), dimension(:), pointer :: stktk => null() + real(ESMF_KIND_R4), dimension(:), pointer :: stkve => null() real(ESMF_KIND_R4), dimension(:), pointer :: factors => null() type(ESMF_Field), dimension(:), pointer :: fields => null() + type(aqm_internal_rate_type), dimension(:), pointer :: rates => null() character(len=ESMF_MAXSTR), dimension(:,:), pointer :: table => null() end type diff --git a/src/shr/aqm_methods.F90 b/src/shr/aqm_methods.F90 index 443512d6..8501c2f5 100644 --- a/src/shr/aqm_methods.F90 +++ b/src/shr/aqm_methods.F90 @@ -91,6 +91,7 @@ LOGICAL FUNCTION DESC3( FNAME ) integer :: localrc integer :: is, ie, js, je + integer :: EMLAYS type(aqm_config_type), pointer :: config ! -- begin @@ -127,7 +128,16 @@ LOGICAL FUNCTION DESC3( FNAME ) ELSE IF ( TRIM( FNAME ) .EQ. TRIM( EMIS_1 ) ) THEN - call aqm_emis_desc("anthropogenic", NLAYS3D, NVARS3D, VNAME3D, UNITS3D) + NLAYS3D = 0 + + call aqm_emis_desc("gbbepx", NLAYS=EMLAYS) + NLAYS3D = MAX(EMLAYS, NLAYS3D) + + call aqm_emis_desc("point-source", NLAYS=EMLAYS) + NLAYS3D = MAX(EMLAYS, NLAYS3D) + + call aqm_emis_desc("anthropogenic", NLAYS=EMLAYS, NVARS=NVARS3D, VNAMES=VNAME3D, UNITS=UNITS3D) + NLAYS3D = MAX(EMLAYS, NLAYS3D) ELSE IF ( TRIM( FNAME ) .EQ. TRIM( GRID_DOT_2D ) ) THEN NVARS3D = 1 @@ -208,7 +218,7 @@ LOGICAL FUNCTION DESC3( FNAME ) GDNAM3D = 'Cubed-Sphere' VGTYP3D = VGSGPN3 ! non-hydrostatic sigma-P - VGTOP3D = 20. + VGTOP3D = 20. * 101.325 ! -- actual sigma levels are not required for AQM since transport ! -- is performed in the atmosphere. Bogus sigma levels are set ! -- to satisfy d(sigma) = 1 @@ -216,22 +226,24 @@ LOGICAL FUNCTION DESC3( FNAME ) VGLVS3D( is ) = DBLE(NLAYS3D + 1 - is) END DO - NVARS3D = 11 + NVARS3D = 14 VNAME3D( 1:NVARS3D ) = & (/ 'JACOBF ', 'JACOBM ', & 'DENSA_J ', 'TA ', & 'QV ', 'QC ', & 'PRES ', 'DENS ', & + 'UWINDA ', 'VWINDA ', & 'ZH ', 'ZF ', & - 'CFRAC_3D ' & + 'CFRAC_3D ', 'PRESF ' & /) UNITS3D( 1:NVARS3D ) = & (/ 'M ', 'M ', & 'KG/M**2 ', 'K ', & 'KG/KG ', 'KG/KG ', & 'Pa ', 'KG/M**3 ', & + 'M/S ', 'M/S ', & 'M ', 'M ', & - 'FRACTION ' & + 'FRACTION ', 'Pa ' & /) call aqm_model_get(config=config, rc=localrc) @@ -353,7 +365,8 @@ logical function envyn(name, description, defaultval, status) case ('CTM_PHOTODIAG') envyn = config % ctm_photodiag case ('CTM_PT3DEMIS') - envyn = aqm_emis_ispresent("gbbepx") + envyn = aqm_emis_ispresent("gbbepx") .or. & + aqm_emis_ispresent("point-source") case ('CTM_GRAV_SETL') envyn = .false. case ('CTM_WBDUST_FENGSHA') @@ -887,6 +900,8 @@ logical function interpx( fname, vname, pname, & end do case ("PRES") p3d => stateIn % prl + case ("PRESF") + p3d => stateIn % pri case ("CFRAC_3D") p3d => stateIn % cldfl case ("PV") @@ -917,6 +932,10 @@ logical function interpx( fname, vname, pname, & p3d => stateIn % tr(:,:,:,config % species % p_atm_qg) set_non_neg = .true. end if + case ("UWINDA") + p3d => stateIn % uwind + case ("VWINDA") + p3d => stateIn % vwind case ("ZF") k = 0 do l = lay0, lay1 diff --git a/src/shr/aqm_model_mod.F90 b/src/shr/aqm_model_mod.F90 index 21b00372..095e62ed 100644 --- a/src/shr/aqm_model_mod.F90 +++ b/src/shr/aqm_model_mod.F90 @@ -86,15 +86,29 @@ subroutine aqm_model_destroy(rc) if (associated(aqm_model(de) % config % species)) then deallocate(aqm_model(de) % config % species, stat=localrc) if (aqm_rc_test((localrc /= 0), & - msg="Failure to allocate model species memory", & + msg="Failure to deallocate model species memory", & file=__FILE__, line=__LINE__, rc=rc)) return nullify(aqm_model(de) % config % species) end if deallocate(aqm_model(de) % config, stat=localrc) if (aqm_rc_test((localrc /= 0), & - msg="Failure to allocate model config memory", & + msg="Failure to deallocate model config memory", & file=__FILE__, line=__LINE__, rc=rc)) return nullify(aqm_model(de) % config) + if (associated(aqm_model(de) % domain % lon)) then + deallocate(aqm_model(de) % domain % lon, stat=localrc) + if (aqm_rc_test((localrc /= 0), & + msg="Failure to deallocate model domain longitudes memory", & + file=__FILE__, line=__LINE__, rc=rc)) return + nullify(aqm_model(de) % domain % lon) + end if + if (associated(aqm_model(de) % domain % lat)) then + deallocate(aqm_model(de) % domain % lat, stat=localrc) + if (aqm_rc_test((localrc /= 0), & + msg="Failure to deallocate model domain latitudes memory", & + file=__FILE__, line=__LINE__, rc=rc)) return + nullify(aqm_model(de) % domain % lat) + end if end if deallocate(aqm_model, stat=localrc) if (aqm_rc_test((localrc /= 0), msg="Failure to allocate model memory", & @@ -195,15 +209,17 @@ subroutine aqm_model_domain_set(minIndexPDe, maxIndexPDe, minIndexPTile, maxInde end subroutine aqm_model_domain_set - subroutine aqm_model_domain_coord_set(coordDim, coord, de, rc) + subroutine aqm_model_domain_coord_set(coordDim, coord, scale, de, rc) - integer, intent(in) :: coordDim - real(AQM_KIND_R8), pointer :: coord(:,:) - integer, optional, intent(in) :: de - integer, optional, intent(out) :: rc + integer, intent(in) :: coordDim + real(AQM_KIND_R8), intent(in) :: coord(:,:) + real(AQM_KIND_R8), optional, intent(in) :: scale + integer, optional, intent(in) :: de + integer, optional, intent(out) :: rc !-- local variables - integer :: deCount, localrc + integer :: deCount, localrc + real(AQM_KIND_R8) :: fscale type(aqm_model_type), pointer :: model !-- begin @@ -215,11 +231,20 @@ subroutine aqm_model_domain_coord_set(coordDim, coord, de, rc) if (deCount < 1) return + fscale = 1._AQM_KIND_R8 + if (present(scale)) fscale = scale + select case (coordDim) case(1) - model % domain % lon => coord + allocate(model % domain % lon, source=fscale * coord, stat=localrc) + if (aqm_rc_test((localrc /= 0), & + msg="Failure to allocate model config memory", & + file=__FILE__, line=__LINE__, rc=rc)) return case(2) - model % domain % lat => coord + allocate(model % domain % lat, source=fscale * coord, stat=localrc) + if (aqm_rc_test((localrc /= 0), & + msg="Failure to allocate model config memory", & + file=__FILE__, line=__LINE__, rc=rc)) return case default call aqm_rc_set(AQM_RC_FAILURE, & msg="coordDim can only be 1 or 2.", & diff --git a/src/shr/aqm_prod_mod.F90 b/src/shr/aqm_prod_mod.F90 index ab844479..ce4b5172 100644 --- a/src/shr/aqm_prod_mod.F90 +++ b/src/shr/aqm_prod_mod.F90 @@ -205,7 +205,7 @@ subroutine aqm_prod_init(model, rc) else if (btest(verbosity,8)) then - write(msgString,'(a,": ",a,": types : none")') trim(name), & + write(msgString,'(a,": ",a,": ",a,": types: none")') trim(name), & trim(rName), trim(pName) call ESMF_LogWrite(msgString, ESMF_LOGMSG_INFO, rc=localrc) if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & diff --git a/src/shr/aqm_state_mod.F90 b/src/shr/aqm_state_mod.F90 index fc2c194f..9dc762fc 100644 --- a/src/shr/aqm_state_mod.F90 +++ b/src/shr/aqm_state_mod.F90 @@ -36,6 +36,7 @@ module aqm_state_mod real(AQM_KIND_R8), dimension(:,:,:), pointer :: dkt => null() real(AQM_KIND_R8), dimension(:,:,:), pointer :: phii => null() real(AQM_KIND_R8), dimension(:,:,:), pointer :: phil => null() + real(AQM_KIND_R8), dimension(:,:,:), pointer :: pri => null() real(AQM_KIND_R8), dimension(:,:,:), pointer :: prl => null() real(AQM_KIND_R8), dimension(:,:,:), pointer :: smois => null() real(AQM_KIND_R8), dimension(:,:,:), pointer :: stemp => null() diff --git a/src/shr/aqm_tools_mod.F90 b/src/shr/aqm_tools_mod.F90 index 3c829c7a..71203b73 100644 --- a/src/shr/aqm_tools_mod.F90 +++ b/src/shr/aqm_tools_mod.F90 @@ -1,12 +1,14 @@ module aqm_tools_mod - use aqm_types_mod, only : AQM_KIND_R4 + use aqm_types_mod, only : AQM_KIND_R4, & + AQM_KIND_R8 implicit none private public :: aqm_units_conv + public :: aqm_gridloc_get contains @@ -120,4 +122,212 @@ subroutine get_index(src, dst, list, isrc, idst, psrc, pdst) end subroutine get_index +! --------------------------------------------------------------------------------- +! grid mapping (spherical geometry) +! --------------------------------------------------------------------------------- +! +! The following source code has been adapted from portions of NASA's MAPL library: +! +! https://github.com/GEOS-ESM/MAPL +! +! MAPL is licensed under the Apache license version 2.0. See notice below: +! +! --------------------------------------------------------------------------------- +! +! NASA Docket No. GSC-15,354-1, and identified as "GEOS-5 GCM Modeling Software +! +! Copyright (C) 2008 United States Government as represented by the Administrator +! of the National Aeronautics and Space Administration. All Rights Reserved. +! +! Licensed under the Apache License, Version 2.0 (the "License"); you may not use +! this file except in compliance with the License. You may obtain a copy of the +! License at +! +! http://www.apache.org/licenses/LICENSE-2.0 +! +! Unless required by applicable law or agreed to in writing, software distributed +! under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR +! CONDITIONS OF ANY KIND, either express or implied. See the License for the +! specific language governing permissions and limitations under the License. +! +! --------------------------------------------------------------------------------- + + subroutine aqm_gridloc_get(center_lons,center_lats,corner_lons,corner_lats,lons,lats,ii,jj,ijcount,rc) + real(AQM_KIND_R8), intent(in) :: center_lats(:,:),center_lons(:,:) + real(AQM_KIND_R8), intent(in) :: corner_lats(:,:),corner_lons(:,:) + real(AQM_KIND_R8), intent(in) :: lons(:),lats(:) + integer, intent(out) :: ii(:),jj(:) + integer, intent(out) :: ijcount + integer, optional, intent(out) :: rc + + ! -- local variables + integer :: npts,i,n,niter,im,jm,ilb,jlb,iub,jub,ifound,jfound + integer :: lold,uold,lnew,unew + logical :: in_region,in_sub_region + + ! -- begin + ii = -1 + jj = -1 + ijcount = 0 + + npts = size(lats) + if (npts /= size(lons)) return + if (npts /= size(ii)) return + if (npts /= size(jj)) return + + im=size(corner_lons,1)-1 + jm=size(corner_lons,2)-1 + niter = max(im,jm) + + do i=1,npts + ifound=-1 + jfound=-1 + ilb=1 + iub=im + jlb=1 + jub=jm + in_region = point_in_polygon([lons(i),lats(i)],[center_lons(ilb,jlb),center_lats(ilb,jlb)], & + [corner_lons(ilb,jlb),corner_lats(ilb,jlb)], & + [corner_lons(iub+1,jlb),corner_lats(iub+1,jlb)], & + [corner_lons(iub+1,jub+1),corner_lats(iub+1,jub+1)], & + [corner_lons(ilb,jub+1),corner_lats(ilb,jub+1)]) + if (in_region) then + ! bisect first dimension + lnew=ilb + unew=iub + do n = 1,niter + lold=lnew + uold=unew + unew=lold+(uold-lold)/2 + in_sub_region = point_in_polygon([lons(i),lats(i)], [center_lons(lnew,jlb),center_lats(lnew,jlb)], & + [corner_lons(lnew,jlb),corner_lats(lnew,jlb)], & + [corner_lons(unew+1,jlb),corner_lats(unew+1,jlb)], & + [corner_lons(unew+1,jub+1),corner_lats(unew+1,jub+1)], & + [corner_lons(lnew,jub+1),corner_lats(lnew,jub+1)]) + if (in_sub_region) then + lnew=lold + unew=unew + else + lnew=unew+1 + unew=uold + end if + if (unew==lnew) then + ifound=unew + exit + end if + enddo + ! bisect 2nd dimension + lnew=jlb + unew=jub + do n = 1,niter + lold=lnew + uold=unew + unew=lold+(uold-lold)/2 + in_sub_region = point_in_polygon([lons(i),lats(i)], [center_lons(ifound,lnew),center_lats(ifound,lnew)] , & + [corner_lons(ifound,lnew),corner_lats(ifound,lnew)], & + [corner_lons(ifound+1,lnew),corner_lats(ifound+1,lnew)], & + [corner_lons(ifound+1,unew+1),corner_lats(ifound+1,unew+1)], & + [corner_lons(ifound,unew+1),corner_lats(ifound,unew+1)]) + if (in_sub_region) then + lnew=lold + unew=unew + else + lnew=unew+1 + unew=uold + end if + if (unew==lnew) then + jfound=unew + exit + end if + enddo + end if + ii(i)=ifound + jj(i)=jfound + if (ifound > 0 .and. jfound > 0) ijcount = ijcount + 1 + enddo + + end subroutine aqm_gridloc_get + + function point_in_polygon(p0,pinside,a1,a2,a3,a4) result(in_poly) + real(AQM_KIND_R8), intent(in) :: p0(2),pinside(2),a1(2),a2(2),a3(2),a4(2) + logical :: in_poly + + real(AQM_KIND_R8) :: p1c(3),p2c(3),a1c(3),a2c(3),a3c(3),a4c(3) + logical :: intersect(4) + p1c=convert_to_cart(p0) + p2c=convert_to_cart(pinside) + a1c=convert_to_cart(a1) + a2c=convert_to_cart(a2) + a3c=convert_to_cart(a3) + a4c=convert_to_cart(a4) + + intersect(1) = lines_intersect(p1c,p2c,a1c,a2c) + intersect(2) = lines_intersect(p1c,p2c,a2c,a3c) + intersect(3) = lines_intersect(p1c,p2c,a3c,a4c) + intersect(4) = lines_intersect(p1c,p2c,a4c,a1c) + if (mod(count(intersect),2)==0) then + in_poly=.true. + else + in_poly=.false. + end if + + end function point_in_polygon + + function lines_intersect(b0,b1,a0,a1) result(intersect) + real(AQM_KIND_R8), intent(in) :: b0(3),b1(3),a0(3),a1(3) + logical :: intersect + real(AQM_KIND_R8) :: p(3),q(3),t(3) + real(AQM_KIND_R8) :: s1,s2,s3,s4 + logical :: signs(4) + + intersect=.false. + q=cross_prod(b0,b1) + p=cross_prod(a0,a1) + t=normal_vect(cross_prod(p,q)) + + s1=dot_product(cross_prod(a0,p),t) + s2=dot_product(cross_prod(a1,p),t) + s3=dot_product(cross_prod(b0,q),t) + s4=dot_product(cross_prod(b1,q),t) + + signs(1) = -s1 <0.d0 + signs(2) = s2 <0.d0 + signs(3) = -s3 < 0.d0 + signs(4) = s4 < 0.d0 + + intersect = ((count(signs)==0) .or. (count(signs)==4)) + + end function lines_intersect + + function normal_vect(vin) result(vout) + real(AQM_KIND_R8), intent(in) :: vin(3) + real(AQM_KIND_R8) :: vout(3) + vout=vin/sqrt(vin(1)*vin(1)+vin(2)*vin(2)+vin(3)*vin(3)) + + end function normal_vect + + function cross_prod(v1,v2) result(vout) + real(AQM_KIND_R8), intent(in) :: v1(3),v2(3) + real(AQM_KIND_R8) :: vout(3) + vout(1)=v1(2)*v2(3)-v1(3)*v2(2) + vout(2)=v1(3)*v2(1)-v1(1)*v2(3) + vout(3)=v1(1)*v2(2)-v1(2)*v2(1) + end function cross_prod + + function convert_to_cart(v) result(xyz) + real(AQM_KIND_R8), intent(in) :: v(2) + real(AQM_KIND_R8) :: xyz(3) + + xyz(1)=cos(v(2))*cos(v(1)) + xyz(2)=cos(v(2))*sin(v(1)) + xyz(3)=sin(v(2)) + + end function convert_to_cart + + function vect_mag(v) result(mag) + real(AQM_KIND_R8), intent(in) :: v(3) + real :: mag + mag = sqrt(v(1)*v(1)+v(2)*v(2)+v(3)*v(3)) + end function vect_mag + end module aqm_tools_mod