diff --git a/README b/README
index 8b13789..d65db6d 100644
--- a/README
+++ b/README
@@ -1 +1,2 @@
-
+Modified branch to account for in-canopy effects on composition/weather
+Export 5 2D canopy data fields and 3 2D photolysis diagnostics (I. Ivanova, 02/14/2024)
diff --git a/aqm_files.cmake b/aqm_files.cmake
index 527cf3d..6727dc8 100644
--- a/aqm_files.cmake
+++ b/aqm_files.cmake
@@ -187,7 +187,6 @@ list(APPEND aqm_CCTM_files
${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
@@ -243,4 +242,6 @@ list(APPEND aqm_CCTM_files
${localCCTM}/ASX_DATA_MOD.F
${localCCTM}/DUST_EMIS.F
${localCCTM}/AERO_PHOTDATA.F
+ ${localCCTM}/phot.F
+ ${localCCTM}/centralized_io_util_module.F
)
diff --git a/examples/aqm.rc b/examples/aqm.rc
index 13246de..70240c0 100644
--- a/examples/aqm.rc
+++ b/examples/aqm.rc
@@ -34,6 +34,31 @@ omi_data: /scratch1/NCEPDEV/nems/Raffaele.Montuoro/dev/aqm/epa/data/omi_cmaq_
# - set to true for cold start
init_concentrations: true
+#
+# Inline Canopy Effects
+#
+canopy_yn: true
+
+canopy_type: canopy
+
+canopy_format: netcdf
+
+canopy_file: /scratch2/NAGAPE/arl/Patrick.C.Campbell/canopy_geofiles/gfs.t12z.geo.08.canopy_regrid.nc
+
+canopy_frequency: static
+
+canopy_species::
+ FCH 1.00000 FCH m
+ FRT 1.00000 FRT 1
+ CLU 1.00000 CLU 1
+ POPU 1.00000 POPU 10000_people/10km2
+ LAIE 1.00000 LAIE 1
+ C1R 1.00000 C1R 1
+ C2R 1.00000 C2R 1
+ C3R 1.00000 C3R 1
+ C4R 1.00000 C4R 1
+::
+
#
# Run options:
#
@@ -83,7 +108,7 @@ ctm_stdout: all
emission_sources: myemis
#
-# Emission type: anthropogenic, biogenic, gbbepx, fengsha
+# Emission type: anthropogenic, biogenic, gbbepx, fengsha, canopy
#
myemis_type: anthropogenic
diff --git a/src/aqm_cap.F90 b/src/aqm_cap.F90
index 15a3231..98c30a0 100644
--- a/src/aqm_cap.F90
+++ b/src/aqm_cap.F90
@@ -9,7 +9,7 @@ module AQM
use aqm_comp_mod
use aqm_const_mod, only: rad_to_deg
-
+
implicit none
! -- import fields
@@ -54,28 +54,36 @@ module AQM
"temperature_of_soil_layer " &
/)
! -- export fields
- integer, parameter :: exportFieldCount = 2
+ integer, parameter :: exportFieldCount = 2 + 5 + 3 !IVAI: add 5 canopy data fields add 3 photdiag arrays
character(len=*), dimension(exportFieldCount), parameter :: &
exportFieldNames = (/ &
"inst_tracer_mass_frac ", &
- "inst_tracer_diag_aod " &
+ "inst_tracer_diag_aod ", &
+ "inst_tracer_diag_coszens ", & !IVAI: photdiag
+ "inst_tracer_diag_jo3o1d ", & !IVAI: photdiag
+ "inst_tracer_diag_jno2 ", & !IVAI: photdiag
+ "inst_tracer_diag_claie ", & !IVAI: canopy via aqm_emis_read
+ "inst_tracer_diag_cfch ", & !IVAI: canopy via aqm_emis_read
+ "inst_tracer_diag_cfrt ", & !IVAI: canopy via aqm_emis_read
+ "inst_tracer_diag_cclu ", & !IVAI: canopy via aqm_emis_read
+ "inst_tracer_diag_cpopu " & !IVAI: canopy via aqm_emis_read
/)
private
public SetServices
-
+
!-----------------------------------------------------------------------------
contains
!-----------------------------------------------------------------------------
-
+
subroutine SetServices(model, rc)
type(ESMF_GridComp) :: model
integer, intent(out) :: rc
! begin
rc = ESMF_SUCCESS
-
+
! the NUOPC model component will register the generic methods
call NUOPC_CompDerive(model, inheritModel, rc=rc)
if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
@@ -120,7 +128,7 @@ subroutine SetServices(model, rc)
return ! bail out
end subroutine
-
+
!-----------------------------------------------------------------------------
subroutine InitializeP0(model, importState, exportState, clock, rc)
@@ -128,7 +136,7 @@ subroutine InitializeP0(model, importState, exportState, clock, rc)
type(ESMF_State) :: importState, exportState
type(ESMF_Clock) :: clock
integer, intent(out) :: rc
-
+
! local variables
integer :: verbosity
character(len=ESMF_MAXSTR) :: msgString, name, rcFile
@@ -195,16 +203,16 @@ subroutine InitializeP0(model, importState, exportState, clock, rc)
line=__LINE__, &
file=__FILE__)) &
return ! bail out
-
+
end subroutine
-
+
!-----------------------------------------------------------------------------
subroutine InitializeP1(model, importState, exportState, clock, rc)
type(ESMF_GridComp) :: model
type(ESMF_State) :: importState, exportState
type(ESMF_Clock) :: clock
integer, intent(out) :: rc
-
+
! begin
rc = ESMF_SUCCESS
@@ -285,7 +293,7 @@ subroutine DataInitialize(model, rc)
return ! bail out
! -- check if import fields are defined
- if (importFieldCount < 1) then
+ if (importFieldCount < 1) then
call ESMF_LogSetError(ESMF_RC_NOT_IMPL, &
msg="This component requires import fields to be defined.", &
line=__LINE__, file=__FILE__, &
@@ -294,7 +302,7 @@ subroutine DataInitialize(model, rc)
end if
! -- check if export fields are defined
- if (exportFieldCount < 1) then
+ if (exportFieldCount < 1) then
call ESMF_LogSetError(ESMF_RC_NOT_IMPL, &
msg="This component requires export fields to be defined.", &
line=__LINE__, file=__FILE__, &
@@ -417,7 +425,7 @@ subroutine DataInitialize(model, rc)
line=__LINE__, &
file=__FILE__)) &
return ! bail out
-
+
do item = 1, 2
call ESMF_GridGetCoord(grid, coordDim=item, staggerloc=ESMF_STAGGERLOC_CENTER, &
localDE=localDe, farrayPtr=coord, rc=rc)
@@ -534,7 +542,7 @@ end subroutine DataInitialize
subroutine ModelAdvance(model, rc)
type(ESMF_GridComp) :: model
integer, intent(out) :: rc
-
+
! local variables
type(ESMF_Clock) :: clock
type(ESMF_State) :: importState, exportState
@@ -548,7 +556,7 @@ subroutine ModelAdvance(model, rc)
! begin
rc = ESMF_SUCCESS
-
+
! get component's information
call NUOPC_CompGet(model, name=name, diagnostic=diagnostic, rc=rc)
if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
diff --git a/src/aqm_comp_mod.F90 b/src/aqm_comp_mod.F90
index c255316..7433755 100644
--- a/src/aqm_comp_mod.F90
+++ b/src/aqm_comp_mod.F90
@@ -325,6 +325,57 @@ subroutine aqm_comp_export(state, fieldNames, rc)
line=__LINE__, &
file=__FILE__)) &
return ! bail
+!IVAI: photdiag CTM_RJ_1 fields
+ case ("inst_tracer_diag_coszens")
+ call ESMF_FieldGet(field, localDe=localDe, farrayPtr=stateOut % coszens, rc=rc)
+ if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
+ line=__LINE__, &
+ file=__FILE__)) &
+ return ! bail
+ case ("inst_tracer_diag_jo3o1d")
+ call ESMF_FieldGet(field, localDe=localDe, farrayPtr=stateOut % jo3o1d, rc=rc)
+ if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
+ line=__LINE__, &
+ file=__FILE__)) &
+ return ! bail
+ case ("inst_tracer_diag_jno2")
+ call ESMF_FieldGet(field, localDe=localDe, farrayPtr=stateOut % jno2, rc=rc)
+ if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
+ line=__LINE__, &
+ file=__FILE__)) &
+ return ! bail
+!IVAI: canopy fields read in via 'aqm_emiss_read'
+ case ("inst_tracer_diag_claie")
+ call ESMF_FieldGet(field, localDe=localDe, farrayPtr=stateOut % claie, rc=rc)
+ if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
+ line=__LINE__, &
+ file=__FILE__)) &
+ return ! bail
+ case ("inst_tracer_diag_cfch")
+ call ESMF_FieldGet(field, localDe=localDe, farrayPtr=stateOut % cfch, rc=rc)
+ if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
+ line=__LINE__, &
+ file=__FILE__)) &
+ return ! bail
+ case ("inst_tracer_diag_cfrt")
+ call ESMF_FieldGet(field, localDe=localDe, farrayPtr=stateOut % cfrt, rc=rc)
+ if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
+ line=__LINE__, &
+ file=__FILE__)) &
+ return ! bail
+ case ("inst_tracer_diag_cclu")
+ call ESMF_FieldGet(field, localDe=localDe, farrayPtr=stateOut % cclu, rc=rc)
+ if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
+ line=__LINE__, &
+ file=__FILE__)) &
+ return ! bail
+ case ("inst_tracer_diag_cpopu")
+ call ESMF_FieldGet(field, localDe=localDe, farrayPtr=stateOut % cpopu, rc=rc)
+ if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
+ line=__LINE__, &
+ file=__FILE__)) &
+ return ! bail
+!IVAI
case default
! -- unused field
end select
diff --git a/src/model/Makefile.am b/src/model/Makefile.am
index cff3967..00d4ae9 100644
--- a/src/model/Makefile.am
+++ b/src/model/Makefile.am
@@ -168,7 +168,6 @@ libCCTM_a_SOURCES += \
$(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 \
@@ -219,7 +218,6 @@ libCCTM_a_SOURCES += \
$(UTIL)/bmatvec.F \
$(UTIL)/findex.f \
$(UTIL)/get_envlist.f \
- $(UTIL)/setup_logdev.F \
$(UTIL)/subhdomain.F \
$(UTIL)/UTILIO_DEFN.F
@@ -250,11 +248,11 @@ libCCTM_a_SOURCES += \
${localCCTM}/PT3D_DEFN.F \
${localCCTM}/PT3D_FIRE_DEFN.F \
${localCCTM}/PT3D_STKS_DEFN.F \
+ $(localCCTM)/phot.F \
${localCCTM}/ASX_DATA_MOD.F \
+ $(localCCTM)/centralized_io_util_module.F \
${localCCTM}/DUST_EMIS.F \
${localCCTM}/AERO_PHOTDATA.F
-
-
libCCTM_a_CPPFLAGS = -DSUBST_FILES_ID=\"FILES_CTM.EXT\"
libCCTM_a_CPPFLAGS += -DSUBST_CONST=\"CONST.EXT\"
libCCTM_a_CPPFLAGS += -DSUBST_EMISPRM=\"EMISPRM.EXT\"
@@ -432,6 +430,10 @@ $(libEMIS)BIOG_EMIS.$(OBJEXT) : \
$(libUTIL)UTILIO_DEFN.$(OBJEXT)
$(libEMIS)cropcal.$(OBJEXT) : \
$(libGRID)HGRD_DEFN.$(OBJEXT) $(libUTIL)UTILIO_DEFN.$(OBJEXT)
+$(libEMIS)DUST_EMIS.$(OBJEXT) : $(ICL)/const/CONST.EXT $(ICL)/filenames/FILES_CTM.EXT \
+ $(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)
$(libEMIS)EMIS_DEFN.$(OBJEXT) : $(ICL)/const/CONST.EXT $(ICL)/filenames/FILES_CTM.EXT \
$(libAERO)AERO_DATA.$(OBJEXT) $(libAERO)AERO_EMIS.$(OBJEXT) \
$(liblocalCCTM)ASX_DATA_MOD.$(OBJEXT) $(libEMIS)BEIS_DEFN.$(OBJEXT) \
diff --git a/src/model/Makefile.in b/src/model/Makefile.in
index 078af39..a787fbd 100644
--- a/src/model/Makefile.in
+++ b/src/model/Makefile.in
@@ -186,7 +186,6 @@ am_libCCTM_a_OBJECTS = $(AERO)/libCCTM_a-AERO_DATA.$(OBJEXT) \
$(PHOT)/libCCTM_a-complex_number_module.$(OBJEXT) \
$(PHOT)/libCCTM_a-CSQY_DATA.$(OBJEXT) \
$(PHOT)/libCCTM_a-opphot.$(OBJEXT) \
- $(PHOT)/libCCTM_a-phot.$(OBJEXT) \
$(PHOT)/libCCTM_a-PHOT_MET_DATA.$(OBJEXT) \
$(PHOT)/libCCTM_a-PHOT_MOD.$(OBJEXT) \
$(PHOT)/libCCTM_a-PHOTOLYSIS_ALBEDO.$(OBJEXT) \
@@ -238,7 +237,9 @@ am_libCCTM_a_OBJECTS = $(AERO)/libCCTM_a-AERO_DATA.$(OBJEXT) \
${localCCTM}/libCCTM_a-PT3D_DEFN.$(OBJEXT) \
${localCCTM}/libCCTM_a-PT3D_FIRE_DEFN.$(OBJEXT) \
${localCCTM}/libCCTM_a-PT3D_STKS_DEFN.$(OBJEXT) \
+ $(localCCTM)/libCCTM_a-phot.$(OBJEXT) \
${localCCTM}/libCCTM_a-ASX_DATA_MOD.$(OBJEXT) \
+ $(localCCTM)/libCCTM_a-centralized_io_util_module.$(OBJEXT) \
${localCCTM}/libCCTM_a-DUST_EMIS.$(OBJEXT) \
${localCCTM}/libCCTM_a-AERO_PHOTDATA.$(OBJEXT)
libCCTM_a_OBJECTS = $(am_libCCTM_a_OBJECTS)
@@ -516,6 +517,7 @@ libCCTM_a_SOURCES = $(AERO)/AERO_DATA.F $(AERO)/aero_depv.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)/phot.F $(localCCTM)/centralized_io_util_module.F \
${localCCTM}/AERO_PHOTDATA.F
# local version of CCTM source files
@@ -886,8 +888,6 @@ $(PHOT)/libCCTM_a-CSQY_DATA.$(OBJEXT): $(PHOT)/$(am__dirstamp) \
$(PHOT)/$(DEPDIR)/$(am__dirstamp)
$(PHOT)/libCCTM_a-opphot.$(OBJEXT): $(PHOT)/$(am__dirstamp) \
$(PHOT)/$(DEPDIR)/$(am__dirstamp)
-$(PHOT)/libCCTM_a-phot.$(OBJEXT): $(PHOT)/$(am__dirstamp) \
- $(PHOT)/$(DEPDIR)/$(am__dirstamp)
$(PHOT)/libCCTM_a-PHOT_MET_DATA.$(OBJEXT): $(PHOT)/$(am__dirstamp) \
$(PHOT)/$(DEPDIR)/$(am__dirstamp)
$(PHOT)/libCCTM_a-PHOT_MOD.$(OBJEXT): $(PHOT)/$(am__dirstamp) \
@@ -1024,6 +1024,9 @@ ${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-phot.$(OBJEXT): \
${localCCTM}/$(am__dirstamp) \
${localCCTM}/$(DEPDIR)/$(am__dirstamp)
${localCCTM}/libCCTM_a-PT3D_FIRE_DEFN.$(OBJEXT): \
@@ -1033,6 +1036,9 @@ ${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-centralized_io_util_module.$(OBJEXT): \
${localCCTM}/$(am__dirstamp) \
${localCCTM}/$(DEPDIR)/$(am__dirstamp)
${localCCTM}/libCCTM_a-DUST_EMIS.$(OBJEXT): \
@@ -1532,11 +1538,17 @@ $(PHOT)/libCCTM_a-opphot.o: $(PHOT)/opphot.F
$(PHOT)/libCCTM_a-opphot.obj: $(PHOT)/opphot.F
$(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o $(PHOT)/libCCTM_a-opphot.obj `if test -f '$(PHOT)/opphot.F'; then $(CYGPATH_W) '$(PHOT)/opphot.F'; else $(CYGPATH_W) '$(srcdir)/$(PHOT)/opphot.F'; fi`
-$(PHOT)/libCCTM_a-phot.o: $(PHOT)/phot.F
- $(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o $(PHOT)/libCCTM_a-phot.o `test -f '$(PHOT)/phot.F' || echo '$(srcdir)/'`$(PHOT)/phot.F
+$(localCCTM)/libCCTM_a-phot.o: $(localCCTM)/phot.F
+ $(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o $(localCCTM)/libCCTM_a-phot.o `test -f '$(localCCTM)/phot.F' || echo '$(srcdir)/'`$(localCCTM)/phot.F
+
+$(localCCTM)/libCCTM_a-phot.obj: $(localCCTM)/phot.F
+ $(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o $(localCCTM)/libCCTM_a-phot.obj `if test -f '$(localCCTM)/phot.F'; then $(CYGPATH_W) '$(localCCTM)/phot.F'; else $(CYGPATH_W) '$(srcdir)/$(localCCTM)/phot.F'; fi`
+
+$(localCCTM)/libCCTM_a-centralized_io_util_module.o: $(localCCTM)/centralized_io_util_module.F
+ $(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o $(localCCTM)/libCCTM_a-centralized_io_util_module.o `test -f '$(localCCTM)/centralized_io_util_module.F' || echo '$(srcdir)/'`$(localCCTM)/centralized_io_util_module.F
-$(PHOT)/libCCTM_a-phot.obj: $(PHOT)/phot.F
- $(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o $(PHOT)/libCCTM_a-phot.obj `if test -f '$(PHOT)/phot.F'; then $(CYGPATH_W) '$(PHOT)/phot.F'; else $(CYGPATH_W) '$(srcdir)/$(PHOT)/phot.F'; fi`
+$(localCCTM)/libCCTM_a-centralized_io_util_module.obj: $(localCCTM)/centralized_io_util_module.F
+ $(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o $(localCCTM)/libCCTM_a-centralized_io_util_module.obj `if test -f '$(localCCTM)/centralized_io_util_module.F'; then $(CYGPATH_W) '$(localCCTM)/centralized_io_util_module.F'; else $(CYGPATH_W) '$(srcdir)/$(localCCTM)/centralized_io_util_module.F'; fi`
$(PHOT)/libCCTM_a-PHOT_MET_DATA.o: $(PHOT)/PHOT_MET_DATA.F
$(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o $(PHOT)/libCCTM_a-PHOT_MET_DATA.o `test -f '$(PHOT)/PHOT_MET_DATA.F' || echo '$(srcdir)/'`$(PHOT)/PHOT_MET_DATA.F
@@ -1622,6 +1634,12 @@ $(VDIFF)/libCCTM_a-aero_sedv.o: $(VDIFF)/aero_sedv.F
$(VDIFF)/libCCTM_a-aero_sedv.obj: $(VDIFF)/aero_sedv.F
$(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o $(VDIFF)/libCCTM_a-aero_sedv.obj `if test -f '$(VDIFF)/aero_sedv.F'; then $(CYGPATH_W) '$(VDIFF)/aero_sedv.F'; else $(CYGPATH_W) '$(srcdir)/$(VDIFF)/aero_sedv.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`
+
$(VDIFF)/libCCTM_a-conv_cgrid.o: $(VDIFF)/conv_cgrid.F
$(AM_V_PPF77)$(F77) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libCCTM_a_CPPFLAGS) $(CPPFLAGS) $(libCCTM_a_FFLAGS) $(FFLAGS) -c -o $(VDIFF)/libCCTM_a-conv_cgrid.o `test -f '$(VDIFF)/conv_cgrid.F' || echo '$(srcdir)/'`$(VDIFF)/conv_cgrid.F
diff --git a/src/model/src/ASX_DATA_MOD.F b/src/model/src/ASX_DATA_MOD.F
old mode 100755
new mode 100644
index 16b3095..48d9258
--- a/src/model/src/ASX_DATA_MOD.F
+++ b/src/model/src/ASX_DATA_MOD.F
@@ -121,6 +121,17 @@ Module ASX_DATA_MOD
Real, Allocatable :: PBL ( :,: ) ! pbl height (m)
Real, Allocatable :: NACL_EMIS( :,: ) ! NACL mass emission rate of particles with d <10 um (g/m2/s)
+!> Inline Canopy Processes
+ Real, Allocatable :: FCH ( :,: ) ! Forest Canopy Height (m)
+ Real, Allocatable :: FRT ( :,: ) ! Forest Fraction
+ Real, Allocatable :: CLU ( :,: ) ! Clumping Index
+ Real, Allocatable :: POPU ( :,: ) ! Population Density (people/10km2)
+ Real, Allocatable :: LAIE ( :,: ) ! ECCC BELD3 Derived LAI (m2/m2)
+ Real, Allocatable :: C1R ( :,: ) ! cumulative LAI fraction hc to 0.75 * hc
+ Real, Allocatable :: C2R ( :,: ) ! cumulative LAI fraction hc to 0.50 * hc
+ Real, Allocatable :: C3R ( :,: ) ! cumulative LAI fraction hc to 0.35 * hc
+ Real, Allocatable :: C4R ( :,: ) ! cumulative LAI fraction hc to 0.20 * hc
+
!> FENGSHA option
Real, Allocatable :: CLAYF ( :,: ) ! Fractional Clay Content
Real, Allocatable :: SANDF ( :,: ) ! Fractional Sand Content
@@ -285,6 +296,10 @@ Module ASX_DATA_MOD
Real, Pointer, Private :: BUFF2D( :,: ) ! 2D temp var
Real, Pointer, Private :: BUFF3D( :,:,: ) ! 3D temp var
+! Canopy option control
+ CHARACTER( 20 ), SAVE :: CTM_CANOPY_SHADE = 'CTM_CANOPY_SHADE '! env var for in-line
+ LOGICAL, PUBLIC, SAVE :: CANOPY_SHADE ! flag in-lining canopy shading
+
! FENGSHA option control
CHARACTER( 18 ), SAVE :: CTM_WBDUST_FENGSHA = 'CTM_WBDUST_FENGSHA' ! env var for in-line
LOGICAL, PUBLIC, SAVE :: FENGSHA ! flag for fengsha option
@@ -630,6 +645,32 @@ Subroutine INIT_MET ( JDATE, JTIME, MOSAIC, ABFLUX, HGBIDI )
ChemMos_Data%SubName = subname
End If
+!> ccccccccccccccccccccc canopy shade option!ccccccccccccccccccccc
+ CANOPY_SHADE = ENVYN( 'CTM_CANOPY_SHADE',
+ & 'Flag for in-line canopy shading',
+ & .FALSE., IOSX )
+
+! IF ( CANOPY_SHADE ) THEN
+! XMSG = 'Using in-line canopy shading option'
+! CALL M3MSG2( XMSG )
+! END IF
+ If ( CANOPY_SHADE ) Then
+ ALLOCATE( Met_Data%FCH ( NCOLS,NROWS ),
+ & Met_Data%FRT ( NCOLS,NROWS ),
+ & Met_Data%CLU ( NCOLS,NROWS ),
+ & Met_Data%POPU ( NCOLS,NROWS ),
+ & Met_Data%LAIE ( NCOLS,NROWS ),
+ & Met_Data%C1R ( NCOLS,NROWS ),
+ & Met_Data%C2R ( NCOLS,NROWS ),
+ & Met_Data%C3R ( NCOLS,NROWS ),
+ & Met_Data%C4R ( NCOLS,NROWS ),
+ & STAT = ALLOCSTAT )
+ If ( ALLOCSTAT .Ne. 0 ) Then
+ XMSG = 'Failure allocating Canopy Shade variables'
+ Call M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ End If
+ End If
+
!> ccccccccccccccccccccc Fengsha option!ccccccccccccccccccccc
FENGSHA = ENVYN( CTM_WBDUST_FENGSHA,
& 'Flag for in-line fengsha ',
@@ -1090,6 +1131,80 @@ Subroutine GET_MET ( JDATE, JTIME, TSTEP, MOSAIC, ABFLUX, HGBIDI )
XMSG = MSG1 // TRIM( VNAME ) // ' from ' // MET_CRO_2D
Call M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
End If
+
+C Canopy vars
+ If ( CANOPY_SHADE ) Then
+ VNAME = 'FCH'
+ If ( .Not. INTERPX( MET_CRO_2D, VNAME, PNAME,
+ & STRTCOLMC2,ENDCOLMC2, STRTROWMC2,ENDROWMC2,1,1,
+ & JDATE, JTIME, Met_Data%FCH ) ) Then
+ XMSG = MSG1 // TRIM( VNAME ) // ' from ' // MET_CRO_2D
+ Call M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ End If
+
+ VNAME = 'FRT'
+ If ( .Not. INTERPX( MET_CRO_2D, VNAME, PNAME,
+ & STRTCOLMC2,ENDCOLMC2,STRTROWMC2,ENDROWMC2,1,1,
+ & JDATE, JTIME, Met_Data%FRT ) ) Then
+ XMSG = MSG1 // TRIM( VNAME ) // ' from ' // MET_CRO_2D
+ Call M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ End If
+
+ VNAME = 'CLU'
+ If ( .Not. INTERPX( MET_CRO_2D, VNAME, PNAME,
+ & STRTCOLMC2,ENDCOLMC2,STRTROWMC2,ENDROWMC2,1,1,
+ & JDATE, JTIME, Met_Data%CLU ) ) Then
+ XMSG = MSG1 // TRIM( VNAME ) // ' from ' // MET_CRO_2D
+ Call M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ End If
+
+ VNAME = 'POPU'
+ If ( .Not. INTERPX( MET_CRO_2D, VNAME, PNAME,
+ & STRTCOLMC2,ENDCOLMC2,STRTROWMC2,ENDROWMC2,1,1,
+ & JDATE, JTIME, Met_Data%POPU ) ) Then
+ XMSG = MSG1 // TRIM( VNAME ) // ' from ' // MET_CRO_2D
+ Call M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ End If
+
+ VNAME = 'LAIE'
+ If ( .Not. INTERPX( MET_CRO_2D, VNAME, PNAME,
+ & STRTCOLMC2,ENDCOLMC2,STRTROWMC2,ENDROWMC2,1,1,
+ & JDATE, JTIME, Met_Data%LAIE ) ) Then
+ Call M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ End If
+
+ VNAME = 'C1R'
+ If ( .Not. INTERPX( MET_CRO_2D, VNAME, PNAME,
+ & STRTCOLMC2,ENDCOLMC2,STRTROWMC2,ENDROWMC2,1,1,
+ & JDATE, JTIME, Met_Data%C1R ) ) Then
+ XMSG = MSG1 // TRIM( VNAME ) // ' from ' // MET_CRO_2D
+ Call M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ End If
+
+ VNAME = 'C2R'
+ If ( .Not. INTERPX( MET_CRO_2D, VNAME, PNAME,
+ & STRTCOLMC2,ENDCOLMC2,STRTROWMC2,ENDROWMC2,1,1,
+ & JDATE, JTIME, Met_Data%C2R ) ) Then
+ XMSG = MSG1 // TRIM( VNAME ) // ' from ' // MET_CRO_2D
+ Call M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ End If
+
+ VNAME = 'C3R'
+ If ( .Not. INTERPX( MET_CRO_2D, VNAME, PNAME,
+ & STRTCOLMC2,ENDCOLMC2,STRTROWMC2,ENDROWMC2,1,1,
+ & JDATE, JTIME, Met_Data%C3R ) ) Then
+ Call M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ End If
+
+ VNAME = 'C4R'
+ If ( .Not. INTERPX( MET_CRO_2D, VNAME, PNAME,
+ & STRTCOLMC2,ENDCOLMC2,STRTROWMC2,ENDROWMC2,1,1,
+ & JDATE, JTIME, Met_Data%C4R ) ) Then
+ XMSG = MSG1 // TRIM( VNAME ) // ' from ' // MET_CRO_2D
+ Call M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ End If
+ End If
+
C FENGSHA vars
If ( FENGSHA ) Then
VNAME = 'CLAYF'
@@ -1124,6 +1239,7 @@ Subroutine GET_MET ( JDATE, JTIME, TSTEP, MOSAIC, ABFLUX, HGBIDI )
Call M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
End If
End If
+
C Soil vars
VNAME = 'SOIM1'
If ( .Not. INTERPX( MET_CRO_2D, VNAME, PNAME,
diff --git a/src/model/src/centralized_io_util_module.F b/src/model/src/centralized_io_util_module.F
new file mode 100644
index 0000000..f5b0653
--- /dev/null
+++ b/src/model/src/centralized_io_util_module.F
@@ -0,0 +1,282 @@
+
+!------------------------------------------------------------------------!
+! 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. !
+!------------------------------------------------------------------------!
+
+!------------------------------------------------------------------------!
+! This module contains utility functions to support centralized I/O
+! implementation
+
+! Revision History:
+! 02/01/19, D. Wong: initial implementation
+! 08/01/19, D. Wong: modified code to work with two-way model
+! 11/20/19, F. Sidi: Modified time to sec to handle negative numbers
+!------------------------------------------------------------------------!
+
+ module centralized_io_util_module
+
+ implicit none
+
+ interface quicksort
+ module procedure quicksort1d,
+ & quicksort2d
+ end interface
+
+ contains
+
+! -------------------------------------------------------------------------
+ recursive subroutine quicksort1d (name, begin, end)
+
+ character (*), intent(out) :: name(:)
+ integer, intent(in) :: begin, end
+
+ integer :: i, j
+ character (50) :: str1, str2
+ logical :: done
+
+ str1 = name( (begin + end) / 2 )
+ i = begin
+ j = end
+ done = .false.
+ do while (.not. done)
+ do while (name(i) < str1)
+ i = i + 1
+ end do
+ do while (str1 < name(j))
+ j = j - 1
+ end do
+ if (i .ge. j) then
+ done = .true.
+ else
+ str2 = name(i)
+ name(i) = name(j)
+ name(j) = str2
+ i = i + 1
+ j = j - 1
+ end if
+ end do
+ if (begin < i-1) call quicksort(name, begin, i-1)
+ if (j+1 < end) call quicksort(name, j+1, end)
+
+ end subroutine quicksort1d
+
+! -------------------------------------------------------------------------
+ recursive subroutine quicksort2d (name, begin, end)
+
+ character (*), intent(out) :: name(:,:)
+ integer, intent(in) :: begin, end
+
+ integer :: i, j, dsize
+ character (50) :: str1, str2(3)
+ logical :: done
+
+ dsize = size(name,2)
+ str1 = name( (begin + end) / 2, 1 )
+ i = begin
+ j = end
+ done = .false.
+ do while (.not. done)
+ do while (name(i,1) < str1)
+ i = i + 1
+ end do
+ do while (str1 < name(j, 1))
+ j = j - 1
+ end do
+ if (i .ge. j) then
+ done = .true.
+ else
+ str2(1:dsize) = name(i,:)
+ name(i,:) = name(j,:)
+ name(j,:) = str2(1:dsize)
+ i = i + 1
+ j = j - 1
+ end if
+ end do
+ if (begin < i-1) call quicksort(name, begin, i-1)
+ if (j+1 < end) call quicksort(name, j+1, end)
+
+ end subroutine quicksort2d
+
+! -------------------------------------------------------------------------
+ function binary_search (name, list, n) result (loc)
+
+ character (*), intent(in) :: name, list(:)
+ integer, intent(in) :: n
+ integer :: loc
+
+ logical :: found
+ integer :: mid_loc, start_loc, end_loc
+
+ start_loc = 1
+ end_loc = n
+ found = .false.
+ do while ((start_loc .le. end_loc) .and. (.not. found))
+ mid_loc = start_loc + (end_loc - start_loc) / 2
+ if (name .lt. list(mid_loc)) then
+ end_loc = mid_loc - 1
+ else if (name .gt. list(mid_loc)) then
+ start_loc = mid_loc + 1
+ else
+ found = .true.
+ end if
+ end do
+
+ if (found) then
+ loc = mid_loc
+ else
+ loc = -1
+ end if
+
+ end function binary_search
+
+! -------------------------------------------------------------------------
+ function search (name, list, n) result (loc)
+
+ character (*), intent(in) :: name, list(:)
+ integer, intent(in) :: n
+ integer :: loc
+
+ logical :: found
+ integer :: lloc
+
+ lloc = 0
+ found = .false.
+ do while ((lloc .le. n) .and. (.not. found))
+ lloc = lloc + 1
+ if (name .eq. list(lloc)) then
+ found = .true.
+ end if
+ end do
+
+ if (found) then
+ loc = lloc
+ else
+ loc = -1
+ end if
+
+ end function search
+
+! -------------------------------------------------------------------------
+ integer function time_to_sec (time)
+
+ integer, intent(in) :: time
+ integer :: neg_time
+ integer :: time_in_sec, hr, min, sec
+
+ if (time .gt. 0) then
+ hr = time / 10000
+ min = mod(time/100, 100)
+ sec = mod(time, 100)
+ time_to_sec = hr * 3600 + min * 60 + sec
+ else
+ neg_time = abs(time)
+ hr = neg_time / 10000
+ min = mod(neg_time/100, 100)
+ sec = mod(neg_time, 100)
+ time_to_sec = -1*(hr * 3600 + min * 60 + sec)
+ end if
+
+ end function time_to_sec
+
+! -------------------------------------------------------------------------
+ integer function time_diff (time1, time2)
+
+ integer, intent(in) :: time1, time2
+
+ time_diff = time_to_sec(time1) - time_to_sec(time2)
+
+ end function time_diff
+
+!--------------------------------------------------------------------------
+ integer function next_day (jday)
+
+! This function determermins the next day for time interpolation
+ implicit none
+
+ integer, intent(in) :: jday
+ integer year, day
+
+ day = MOD(jday,1000)
+ year = INT(jday/1000)
+
+ If( day .LT. 365 ) Then
+ next_day = jday+1
+ Else
+ If( MOD(year,4) .Eq. 0 .And. MOD(year,100) .Ne. 0 ) Then
+! Leap Year
+ If( day .Eq. 365 ) Then
+ next_day = jday + 1
+ Else
+ next_day = (INT(jday/1000)+1)*1000+1
+ End If
+ Else If(MOD(year,400) .Eq. 0 ) Then
+! also a leap year, e.g. 2000 but not 2100
+ If( day .Eq. 365 ) Then
+ next_day = jday + 1
+ Else
+ next_day = (INT(jday/1000)+1)*1000+1
+ End If
+ Else
+! not a leap year
+ next_day = (INT(jday/1000)+1)*1000+1
+ End If
+ End If
+
+ end function next_day
+
+!--------------------------------------------------------------------------
+
+ function IntegrateTrapezoid(x, y)
+ !! Calculates the integral of an array y with respect to x using the trapezoid
+ !! approximation. Note that the mesh spacing of x does not have to be uniform.
+ real, intent(in) :: x(:) !! Variable x
+ real, intent(in) :: y(size(x)) !! Function y(x)
+ real :: IntegrateTrapezoid !! Integral ∫y(x)·dx
+ ! Integrate using the trapezoidal rule
+ associate(n => size(x))
+ IntegrateTrapezoid = sum((y(1+1:n-0) + y(1+0:n-1))*(x(1+1:n-0) - x(1+0:n-1)))/2
+ end associate
+ end function
+
+! ---------------------------------------------------------------------------
+
+ function interp_linear1_internal(x,y,xout) result(yout)
+ !! Interpolates for the y value at the desired x value,
+ !! given x and y values around the desired point.
+
+ implicit none
+
+ real, intent(IN) :: x(2), y(2), xout
+ real :: yout
+ real :: alph
+
+ if ( xout .lt. x(1) .or. xout .gt. x(2) ) then
+ write(*,*) "interp1: xout < x0 or xout > x1 !"
+ write(*,*) "xout = ",xout
+ write(*,*) "x0 = ",x(1)
+ write(*,*) "x1 = ",x(2)
+ stop
+ end if
+
+ alph = (xout - x(1)) / (x(2) - x(1))
+ yout = y(1) + alph*(y(2) - y(1))
+
+ return
+
+ end function interp_linear1_internal
+
+ end module centralized_io_util_module
diff --git a/src/model/src/phot.F b/src/model/src/phot.F
new file mode 100644
index 0000000..1b8fd87
--- /dev/null
+++ b/src/model/src/phot.F
@@ -0,0 +1,1441 @@
+
+!------------------------------------------------------------------------!
+! 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. !
+!------------------------------------------------------------------------!
+
+
+! RCS file, release, date & time of last delta, author, state, [and locker]
+! $Header: /project/yoj/arc/CCTM/src/phot/phot_inline/phot.F,v 1.7 2011/10/21 16:11:28 yoj Exp $
+
+! what(1) key, module and SID; SCCS file; date and time of last delta:
+! %W% %P% %G% %U%
+
+!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+ SUBROUTINE PHOT ( MDATE, MTIME, JDATE, JTIME, DTSTEP, RJ )
+
+!-----------------------------------------------------------------------
+!
+! Function: Calculates the photolysis rate constant to be used by the
+! chemical solver. It calculates these rates at each gridcell using
+! codes adapted from JPROC. Cloud correction now called within the
+! loops over MY-ROW & MY_COLS
+!
+! Preconditions: HGRD_INIT() called from PAR_INIT, which is called from
+! DRIVER
+!
+! Subroutines/Functions called: INIT3, M3EXIT, SUBHFILE, CGRID_MAP,
+! OPPHOT, LOAD_CSQY_DATA, LOAD_OPTICS_DATA, INITIALIZE_ALBEDO,
+! GET_PHOT_MET, UPDATE_SUN, GET_ALBEDO, GET_DROPLET_OPTICS,
+! GET_ICE_OPTICS, GET_AGGREGATE_OPTICS, CLEAR_HYDROMETEOR_OPTICS,
+! GET_AERO_DATA, O3TOTCOL, and NEW_OPTICS
+!
+! Revision History.
+! Started 10/08/2004 with existing PHOT and JPROC coded by
+! Dr. Francis S. Binkowski
+! Carolina Environmental Program
+! University of North Carolina at Chapel Hill
+! email: frank_binkowski@unc.edu
+! August 2005, Sarav Arunachalam, CEP, UNC-CH
+! - Minor revisions while integrating with CMAQ
+! - Error check for NPHOTS added (this version works only for SAPRC-99)
+! - Added creation of new file CTM_RJ_1 to write out RJ values
+! for O3 and NO2 (both clear sky and cloud effects), and
+! ETOT_SFC, TAU_AERO, TAU_TOT and TAUO3_TOP values for 7 wavelengths
+! June 2007, David Wong
+! -- inline with CMAQ
+! - declare RJ as assumed shape array to match with the caller routine
+! - allow PE 0 only to open the output file
+! - output species: NO2_CLOUD and O3_CLOUD with AMISS value when all cells
+! are dark and JTIME_CHK = 0
+! - output species: NO2_CLOUD and O3_CLOUD with AMISS value when CLDATT is
+! 0 and JTIME_CHK = 0
+! December 2007, Francis Binkowski
+! code has been modified to call the new on-line version that
+! has the cloud effects built in. new photolysis routine to
+! replace PHOT in CMAQ
+! January 2008, Shawn Roselle
+! - reformatted for inclusion in CMAQ
+! - added additional 3-d photolysis rate diagnostic file
+! - moved code for opening the diagnostic files to a separate subroutine
+! - moved aerosol pointer evaluation to a FORTRAN module
+! - simplified code for writing the diagnostic file
+! - changed code to call NEW_OPTICS twice, once for clear sky and
+! another time for the cloudy fraction of the grid cell. RJ's are
+! computed based on the cloud fraction weighting.
+! March 2011, Bill Hutzell
+! - enable wavelength dependent arrays to have an allocatable number
+! of wavelength bins
+! - added data structure and algorithm to compute a surface albedo that
+! depends on time and landuse catagory based on work by John Striecher
+! (AMAD/USEPA)
+! - revised writing to RJ1 file to include surface albedo
+! - moved photolysis and opacity data from CSQY module to an ASCII input
+! file
+! - added routine called LOAD_REF_DATA (inside the PHOT_MOD module) that i
+! reads this input file
+! - added call to a routine called AERO_PHOTDATA that returns opacity data
+! on the aerosol distribution
+! - revised NEW_OPTICS' arguments based on aerosol redesign in CMAQ
+! version 5.0
+! March 29, 2011 S.Roselle
+! - Replaced I/O API include files with UTILIO_DEFN
+! 07 Jul 14 B.Hutzell: replaced mechanism include file(s) with fortran module
+! 26 Sep 14 B.Hutzell: 1) moved calculation of surface albedo to its own
+! fortran module
+! 2) changed loading procedure for loading optical data;
+! two files now used
+! 3) reading and calculation of met and geo data
+! now acomplished by a fortran module
+! 4) changed description and accounting of cloud effects
+! from 2D liquid water clouds to 3D resolved and subgrid
+! clouds with multi-phases of water
+! 5) inserted calculation of aerosol optical properites via
+! fortran module to improve efficiency in radiative
+! transfer solution
+! 6) moved the O3TOTCOL routine from the PHOT_MOD to simplify
+! the NEW_OPTICS routine
+! 7) Several miscellaneous changes attempting to improve efficiency
+! June 10 15 J.Young: Modified diagnostic output timestamp to fix for other than one
+! hour time steps.
+! Aug 12, 15 D. Wong: Replaced MYPE with IO_PE_INCLUSIVE for parallel I/O implementation
+
+!----------------------------------------------------------------------
+
+C...modules
+
+ USE RXNS_DATA ! chemistry varaibles and data
+ USE CGRID_SPCS ! CGRID species number and offsets
+ USE PCGRID_DEFN ! get cgrid
+ USE UTILIO_DEFN
+ USE AERO_DATA ! describes aerosol distribution
+ USE PHOT_MOD ! photolysis in-line module - inherits CSQY_DATA module
+ USE AERO_PHOTDATA ! arrays and routines for aerosol dimensions and refractive indices
+ USE PHOTOLYSIS_ALBEDO ! surface albedo data and routines
+ USE PHOT_MET_DATA ! Met and Grid data
+ USE CLOUD_OPTICS ! data and routines for optics of cloud hydrometeors
+! USE STRATOS_O3_MINFRACS ! annual minimum fraction of ozone column density above Pressure TOP
+! USE SEAS_STRAT_O3_FRACS ! monthly minimum fraction of ozone column density above Pressure TOP
+ USE SEAS_STRAT_O3_MIN ! monthly minimum fraction of ozone column density above Pressure TOP
+!Used for canopy shade calculation
+ USE ASX_DATA_MOD, ONLY : MET_DATA !uses met data
+ USE centralized_io_util_module, ONLY: IntegrateTrapezoid, interp_linear1_internal !basic utilities
+
+#ifdef parallel
+ USE SE_MODULES ! stenex (using SE_UTIL_MODULE)
+#else
+ USE NOOP_MODULES ! stenex (using NOOP_UTIL_MODULE)
+#endif
+
+ IMPLICIT NONE
+
+!...include files
+
+ INCLUDE SUBST_FILES_ID ! file name parameters
+! INCLUDE SUBST_CONST ! physical constants--moved to PHOT_MOD.
+
+!...arguments
+
+ INTEGER, INTENT( IN ) :: MDATE ! "centered" Julian date (YYYYDDD)
+ INTEGER, INTENT( IN ) :: MTIME ! "centered" time (HHMMSS)
+ INTEGER, INTENT( IN ) :: JDATE ! current Julian date (YYYYDDD)
+ INTEGER, INTENT( IN ) :: JTIME ! current time (HHMMSS)
+ INTEGER, INTENT( IN ) :: DTSTEP( : ) ! time step vector (HHMMSS)
+
+! REAL RJ( NCOLS,NROWS,NLAYS, NPHOTAB )
+ REAL, INTENT( OUT ) :: RJ( :,:,:,: ) ! gridded J-values (1/min units)
+
+! REAL CGRID( NCOLS,NROWS,NLAYS, * ) ! Conc array
+ REAL, SAVE, POINTER :: CGRID( :,:,:,: ) ! species concentrations
+
+!...parameters
+
+ LOGICAL, PARAMETER :: CLDATT = .TRUE. ! include cloud attenuation
+
+ REAL, PARAMETER :: DENS_CONV = ( 1.0E+03 * AVO / MWAIR ) * 1.0E-06 ! convert from kg/m**3 to #/cc
+ REAL, PARAMETER :: PPM_MCM3 = 1.0E-06 ! convert from ppm to molecules / cc mol_Spec/mol_Air = ppm * 1E-06
+ REAL, PARAMETER :: PRES_CONV = 1.0 / STDATMPA ! conversion factor Pa to atm
+ REAL, PARAMETER :: ZTOA = 50.0E3 ! height of top of atmosphere [ m ] (=50km)
+ ! based a 2005 WRF model Documentation
+
+ REAL, PARAMETER :: EPSLON = 1.0E-30 ! Small number
+
+!...external functions: none
+
+!...local variables
+
+ LOGICAL, SAVE :: FIRSTIME = .TRUE. ! Flag for first call to PHOT
+ LOGICAL, SAVE :: PHOTDIAG ! Flag for PHOTDIAG file
+
+ LOGICAL, SAVE :: CALL_INIT_ALBEDO = .TRUE.
+ LOGICAL, SAVE :: CALL_GET_ALBEDO = .TRUE.
+
+ LOGICAL :: ZERO_ICE
+
+ CHARACTER( 3 ), ALLOCATABLE, SAVE :: WLTXT( : )
+ CHARACTER( 16 ) :: VARNM
+ CHARACTER( 16 ), SAVE :: PNAME = 'PHOT'
+ CHARACTER( 16 ), SAVE :: CTM_PHOTDIAG = 'CTM_PHOTDIAG'
+
+ CHARACTER( 80 ) :: VARDESC ! environment variable description
+ CHARACTER( 240 ) :: XMSG = ' '
+
+ INTEGER, SAVE :: LOGDEV
+ INTEGER, SAVE :: LGC_O3 ! pointer to O3 in CGRID
+ INTEGER, SAVE :: LGC_NO2 ! pointer to NO2 in CGRID
+ INTEGER, SAVE :: TSTEP ! output timestep in sec
+
+ INTEGER ESTAT ! status from environment var check
+ INTEGER IPHOT ! photolysis rate loop index
+ INTEGER ROW
+ INTEGER COL
+ INTEGER LEV
+ INTEGER SPC
+ INTEGER IWL
+ INTEGER L
+ INTEGER V, N, MODE
+ LOGICAL JTIME_CHK ! To check for JTIME to write RJ values
+ INTEGER ODATE ! output date
+ INTEGER OTIME ! output time
+
+ INTEGER ALLOCSTAT
+
+ INTEGER, SAVE :: TDATE
+ INTEGER, SAVE :: GXOFF, GYOFF ! global origin offset from file
+ INTEGER, SAVE :: PECOL_OFFSET ! Local Column Offset for processor
+ INTEGER, SAVE :: PEROW_OFFSET ! Local Column Offset for processor
+ INTEGER, SAVE :: TSTEP_COUNT ! counter between calls to write diagnostics
+
+! for INTERPX
+ INTEGER, SAVE :: STRTCOLGC2, ENDCOLGC2, STRTROWGC2, ENDROWGC2
+ INTEGER, SAVE :: STRTCOLMC2, ENDCOLMC2, STRTROWMC2, ENDROWMC2
+ INTEGER, SAVE :: STRTCOLMC3, ENDCOLMC3, STRTROWMC3, ENDROWMC3
+
+ REAL CURRHR ! current GMT hour
+ REAL JULIAN_DAY ! time of year [days]
+ REAL CURRHR_LST ! local standard time at each grid cell
+ REAL CTOP ! cloud top in single dimension
+ REAL CBASE ! cloud base in single dimension
+ REAL ZLEV ! height in single dimension
+ REAL ZEN ! cosine of zenith angle
+ REAL SINLAT ! sine of latitude
+ REAL COSLAT ! cosine of latitude
+ REAL RSQD ! square of soldist
+ REAL ZSFC ! surface height (msl) [ m ]
+ REAL EQT ! equation of time
+ REAL SOLDIST ! solar distance [ au ]
+ REAL SINDEC ! sine of the solar declination
+ REAL COSDEC ! cosine of the solar declination
+ REAL COSZEN ! working cosine of the solar zenith angle
+ REAL SINZEN ! working sine of the solar zenith angle
+ REAL LATCR ! local latitude
+ REAL LONCR ! local longitude
+ REAL OWATER_FRAC ! Open water fraction
+ REAL SNOW_FRAC ! Snow fractional coverage
+ REAL SEAICE_FRAC ! Sea Ice fraction
+ REAL RES_SKY_REFLECT ! reflection coefficient based on resolved sky
+ REAL RES_SKY_TRANS ! diffuse transmission coefficient based on resolved sky
+ REAL RES_SKY_TRANSD ! direct transmission coefficient based on resolved sky
+
+ REAL :: TOTAL_O3_COLUMN ! total ozone column density, DU
+
+ REAL, SAVE :: JYEAR = 0.0 ! year
+ REAL, SAVE :: JD_STRAT_O3MIN = 0.0 ! Julian day (YYYYDDD) of min fraction for stratos ozone
+
+ INTEGER, PARAMETER :: DAYS( 12 ) = (/ 0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30 /)
+ INTEGER, SAVE :: IMONTH = 0
+
+ REAL, ALLOCATABLE, SAVE :: ETOT_SFC ( : ) ! total downward irradiance at sfc [ Watts / m**2 ]
+ REAL, ALLOCATABLE, SAVE :: TAUO3_TOP( : ) ! optical depth of ozone above model domain
+ REAL, ALLOCATABLE, SAVE :: TAU_RAY ( : ) ! Rayleigh optical depth above model domain
+ REAL, ALLOCATABLE, SAVE :: TAUC_AERO( :,: ) ! aerosol optical depth at layer bottom
+ REAL, ALLOCATABLE, SAVE :: TAU_TOT ( :,: ) ! total optical depth at layer bottom
+ REAL, ALLOCATABLE, SAVE :: TAU_CLOUD( :,: ) ! cloud optical depth at layer bottom
+
+ REAL, ALLOCATABLE, SAVE :: SSA ( : ) ! aerosol single scattering albedo, column average
+
+ REAL MSCALE ! combined factor to scale ppm to Molecules / cm**3
+ ! and correct for ambient temperaure and pressure
+
+! FSB new arrays for new on-line cloud version
+
+ REAL, ALLOCATABLE, SAVE :: LWC ( : ) ! cloud liquid water content [ g/m**3 ]
+ REAL, ALLOCATABLE, SAVE :: RWC ( : ) ! rain water content [ g/m**3 ]
+ REAL, ALLOCATABLE, SAVE :: IWC ( : ) ! ice liquid water content [ g/m**3 ]
+ REAL, ALLOCATABLE, SAVE :: SWC ( : ) ! snow content [ g/m**3 ]
+ REAL, ALLOCATABLE, SAVE :: GWC ( : ) ! graupel content [ g/m**3 ]
+ REAL, ALLOCATABLE, SAVE :: CLDFRAC( : ) ! fractional cloud cover
+ REAL, ALLOCATABLE, SAVE :: BLKPRS ( : ) ! Air pressure in [ Pa ]
+ REAL, ALLOCATABLE, SAVE :: BLKTA ( : ) ! Air temperature [ K ]
+ REAL, ALLOCATABLE, SAVE :: BLKDENS( : ) ! Air density [ molecules / m**3 ]
+ REAL, ALLOCATABLE, SAVE :: BLKZH ( : ) ! layer half-height [ m ]
+ REAL, ALLOCATABLE, SAVE :: BLKO3 ( : ) ! O3 concentration [ molecules/cm**3 ]
+ REAL, ALLOCATABLE, SAVE :: BLKNO2 ( : ) ! NO2 concentration [ molecules/cm**3 ]
+ REAL, ALLOCATABLE, SAVE :: BLKZF ( : ) ! layer full-height [ m ]
+
+ REAL, ALLOCATABLE, SAVE :: BLKRJ_RES( :, : ) ! photolysis rates
+ REAL, ALLOCATABLE, SAVE :: BLKRJ_ACM( :, : ) ! photolysis rates
+
+ LOGICAL, ALLOCATABLE, SAVE :: CLOUDS( : ) ! Does layer have clouds?
+ LOGICAL :: NEW_PROFILE ! Has atmospheric temperature and density profile changed?
+ LOGICAL :: DARK ! Are this processor's cells in darkness?
+
+! Canopy in-line control
+ CHARACTER( 20 ), SAVE :: CTM_CANOPY_SHADE = 'CTM_CANOPY_SHADE ' ! env var for in-line
+ LOGICAL, SAVE :: CANOPY_SHADE ! flag in-lining canopy shading
+
+! Canopy arrays
+ REAL, ALLOCATABLE, SAVE :: RJ_CORR_C1R ( :, :) ! canopy shading correction to J-values (hc to 0.75*hc)
+ REAL, ALLOCATABLE, SAVE :: RJ_CORR_C2R ( :, :) ! canopy shading correction to J-values (hc to 0.50*hc)
+ REAL, ALLOCATABLE, SAVE :: RJ_CORR_C3R ( :, :) ! canopy shading correction to J-values (hc to 0.35*hc)
+ REAL, ALLOCATABLE, SAVE :: RJ_CORR_C4R ( :, :) ! canopy shading correction to J-values (hc to 0.20*hc)
+ REAL, ALLOCATABLE, SAVE :: RJ_CORR_BOT ( :, :) ! canopy shading correction to J-values (0.20*hc to bottom)
+ REAL, ALLOCATABLE, SAVE :: RJ_CORR ( :, :) ! total/integrated canopy shading correction to J-values
+ REAL, ALLOCATABLE, SAVE :: ZCANX ( : ) ! canopy heights[m]
+ REAL, ALLOCATABLE, SAVE :: RJ_CORRX ( : ) ! canopy height dependent photolysis attenuation factor
+ REAL :: XCAN ( 2 ) ! canopy height interpolation bounds
+ REAL :: YCAN ( 2 ) ! photolysisattenuation interpolation bounds
+ REAL ZFL, ZCAN, XCANOUT, BOTCAN ! local canopy variables
+ INTEGER :: COUNTCAN, KCAN
+ INTEGER, PARAMETER :: MAXCAN = 1000 ! Declare local maximum canopy layers
+
+!...Variables for diagnostic outputs
+
+ REAL, ALLOCATABLE, SAVE :: N_EXCEED_TROPO3( :,: ) ! Number of adjustments tropospheric ozone optical depth
+
+ REAL, ALLOCATABLE, SAVE :: TOTAL_OC( :,: ) ! total ozone column [DU]
+ REAL, ALLOCATABLE, SAVE :: TROPO_OC( :,: ) ! tropospheric ozone column [DU]
+ REAL, ALLOCATABLE, SAVE :: TROPO_O3_EXCEED( :,: ) ! Factor used to adjust tropospheric ozone optical depth
+ REAL, ALLOCATABLE, SAVE :: TRANSMIS_DIFFUSE( :,: ) ! diffuse transmission coefficient at surface
+ REAL, ALLOCATABLE, SAVE :: TRANSMIS_DIRECT( :,: ) ! direct transmission coefficient at surface
+ REAL, ALLOCATABLE, SAVE :: REFLECT_COEFF( :,: ) ! reflection coefficient at top of atmosphere
+ REAL, ALLOCATABLE, SAVE :: ETOT_SFC_WL ( :,:,: ) ! total downward irradiance at sfc [ Watts / m**2 ]
+ REAL, ALLOCATABLE, SAVE :: TAU_AERO_WL ( :,:,: ) ! total aerosol optical depth
+ REAL, ALLOCATABLE, SAVE :: TAU_CLOUD_WL( :,:,: ) ! total cloud optical depth
+ REAL, ALLOCATABLE, SAVE :: CLR_TRANSMISSION( :,: ) ! diffuse transmission coefficient of clouds
+ REAL, ALLOCATABLE, SAVE :: CLR_REFLECTION ( :,: ) ! reflection coefficient of cloud
+ REAL, ALLOCATABLE, SAVE :: CLR_TRANS_DIRECT( :,: ) ! direct transmission coefficient of clouds
+#ifdef phot_debug
+ REAL, ALLOCATABLE, SAVE :: ASY_CLOUD_WL( :,:,: ) ! columm average of cloud asymmetry factor
+ REAL, ALLOCATABLE, SAVE :: SSA_CLOUD_WL( :,:,: ) ! columm average of cloud single scattering albedo
+#endif
+ REAL, ALLOCATABLE, SAVE :: TAU_TOT_WL ( :,:,: ) ! total optical depth
+ REAL, ALLOCATABLE, SAVE :: TAUO3_TOP_WL( :,:,: ) ! optical depth of ozone above model domain
+
+ REAL, ALLOCATABLE, SAVE :: AERO_SSA ( :,:,:,: ) ! aerosol single scattering albedo
+ REAL, ALLOCATABLE, SAVE :: AERO_ASYM ( :,:,:,: ) ! aerosol asymmetry factor
+ REAL, ALLOCATABLE, SAVE :: TAU ( :,:,:,: ) ! optical depth
+ REAL, ALLOCATABLE, SAVE :: TAU_AERO ( :,:,:,: ) ! aerosol optical depth
+ REAL, ALLOCATABLE, SAVE :: ACTINIC_FX( :,:,:,: ) ! net actinic flux [watts/m**2]
+
+ INTEGER IOSX ! i/o and allocate memory status
+
+ INTERFACE
+ SUBROUTINE O3TOTCOL ( LATITUDE, LONGITUDE, JDATE, OZONE )
+ INTEGER, INTENT( IN ) :: JDATE ! Julian day of the year (yyyyddd)
+ REAL, INTENT( IN ) :: LATITUDE ! latitude of point on earth's surface
+ REAL, INTENT( IN ) :: LONGITUDE ! longitude of point on earth's surface
+ REAL, INTENT( INOUT ) :: OZONE ! total column ozone [DU]
+ END SUBROUTINE O3TOTCOL
+ END INTERFACE
+
+! ----------------------------------------------------------------------
+
+ IF ( FIRSTIME ) THEN
+
+C In-line canopy shading option? (default = false)
+
+ CANOPY_SHADE = ENVYN( 'CTM_CANOPY_SHADE',
+ & 'Flag for in-line canopy shading',
+ & .FALSE., IOSX )
+
+ IF ( CANOPY_SHADE ) THEN
+ XMSG = 'Using in-line canopy shading option'
+ CALL M3MSG2( XMSG )
+ END IF
+
+ FIRSTIME = .FALSE.
+ LOGDEV = INIT3()
+
+ TSTEP = TIME2SEC( DTSTEP( 1 ) ) ! output timestep for phot diagnostic files
+
+ CGRID => PCGRID( 1:MY_NCOLS,1:MY_NROWS,:,: )
+
+!...Get photolysis rate diagnostic file flag
+
+ PHOTDIAG = .FALSE. ! default
+ VARDESC= 'Flag for writing the photolysis rate diagnostic file'
+ PHOTDIAG = ENVYN( CTM_PHOTDIAG, VARDESC, PHOTDIAG, ESTAT )
+ IF ( ESTAT .NE. 0 ) WRITE( LOGDEV, '(5X, A)' ) VARDESC
+ IF ( ESTAT .EQ. 1 ) THEN
+ XMSG = 'Environment variable improperly formatted'
+ CALL M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT2 )
+ ELSE IF ( ESTAT .EQ. -1 ) THEN
+ XMSG =
+ & 'Environment variable set, but empty ... Using default:'
+ WRITE( LOGDEV, '(5X, A, I9)' ) XMSG, JTIME
+ ELSE IF ( ESTAT .EQ. -2 ) THEN
+ XMSG = 'Environment variable not set ... Using default:'
+ WRITE( LOGDEV, '(5X, A, I9)' ) XMSG, JTIME
+ END IF
+
+!...Get met file offsets
+
+ CALL SUBHFILE ( GRID_CRO_2D, GXOFF, GYOFF,
+ & STRTCOLGC2, ENDCOLGC2, STRTROWGC2, ENDROWGC2 )
+ CALL SUBHFILE ( MET_CRO_2D, GXOFF, GYOFF,
+ & STRTCOLMC2, ENDCOLMC2, STRTROWMC2, ENDROWMC2 )
+ CALL SUBHFILE ( MET_CRO_3D, GXOFF, GYOFF,
+ & STRTCOLMC3, ENDCOLMC3, STRTROWMC3, ENDROWMC3 )
+
+ PECOL_OFFSET = COLSD_PE( 1, MYPE+1 ) - 1
+ PEROW_OFFSET = ROWSD_PE( 1, MYPE+1 ) - 1
+
+ CALL LOAD_CSQY_DATA( )
+
+ CALL LOAD_OPTICS_DATA( )
+
+!...Allocate array needed to calculation aerosol and cloud optical properties
+
+ CALL INIT_AERO_DATA( )
+
+ CALL INIT_CLOUD_OPTICS( )
+
+!...Allocate and initialize new canopy arrays
+ IF ( CANOPY_SHADE ) THEN
+ ALLOCATE( RJ_CORRX ( MAXCAN ) )
+ ALLOCATE( ZCANX ( MAXCAN ) )
+
+ ALLOCATE( RJ_CORR_C1R ( NCOLS, NROWS ),
+ & RJ_CORR_C2R ( NCOLS, NROWS ),
+ & RJ_CORR_C3R ( NCOLS, NROWS ),
+ & RJ_CORR_C4R ( NCOLS, NROWS ),
+ & RJ_CORR_BOT ( NCOLS, NROWS ),
+ & RJ_CORR ( NCOLS, NROWS ), STAT = ALLOCSTAT )
+ IF ( ALLOCSTAT .NE. 0 ) THEN
+ XMSG = 'Failure allocating canopy photolysis rate correction arrays'
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ END IF
+
+ RJ_CORRX=0.0
+ ZCANX=0.0
+ RJ_CORR_C1R=0.0
+ RJ_CORR_C2R=0.0
+ RJ_CORR_C3R=0.0
+ RJ_CORR_C4R=0.0
+ RJ_CORR_BOT=0.0
+ RJ_CORR=0.0
+ END IF
+
+!...Initialize Surface albedo method
+
+ IF ( .NOT. INITIALIZE_ALBEDO( JDATE, JTIME, LOGDEV ) ) THEN
+ XMSG = 'Failure initializing photolysis surface albedo algorithm'
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ END IF
+
+ ALLOCATE( ETOT_SFC ( NWL ) )
+
+ ALLOCATE( LWC ( NLAYS ) )
+ ALLOCATE( RWC ( NLAYS ) )
+ ALLOCATE( IWC ( NLAYS ) )
+ ALLOCATE( SWC ( NLAYS ) )
+ ALLOCATE( GWC ( NLAYS ) )
+ ALLOCATE( BLKPRS ( NLAYS ) )
+ ALLOCATE( BLKTA ( NLAYS ) )
+ ALLOCATE( BLKDENS( NLAYS ) )
+ ALLOCATE( BLKZH ( NLAYS ) )
+ ALLOCATE( BLKO3 ( NLAYS ) )
+ ALLOCATE( BLKNO2 ( NLAYS ) )
+ ALLOCATE( BLKZF ( NLAYS+1 ) )
+ ALLOCATE( CLOUDS ( NLAYS ) )
+ ALLOCATE( CLDFRAC( NLAYS ) )
+
+ ALLOCATE( BLKRJ_RES( NLAYS,NPHOTAB ) )
+ ALLOCATE( BLKRJ_ACM( NLAYS,NPHOTAB ) )
+
+ ALLOCATE( TAUO3_TOP( NWL ) )
+ ALLOCATE( TAU_RAY ( NWL ) )
+ ALLOCATE( SSA ( NWL ) )
+
+ ALLOCATE( TAU_CLOUD( NLAYS,NWL ) )
+ ALLOCATE( TAUC_AERO( NLAYS,NWL ) )
+ ALLOCATE( TAU_TOT ( NLAYS,NWL ) )
+
+ ALLOCATE( TOTAL_OC ( NCOLS,NROWS ) )
+
+ IF ( PHOTDIAG ) THEN
+ ALLOCATE( TROPO_OC ( NCOLS,NROWS ) )
+ ALLOCATE( TROPO_O3_EXCEED( NCOLS,NROWS ) )
+ ALLOCATE( N_EXCEED_TROPO3( NCOLS,NROWS ) )
+ ALLOCATE( TRANSMIS_DIFFUSE( NCOLS,NROWS ) )
+ ALLOCATE( TRANSMIS_DIRECT ( NCOLS,NROWS ) )
+ ALLOCATE( REFLECT_COEFF ( NCOLS,NROWS ) )
+ ALLOCATE( CLR_TRANSMISSION( NCOLS,NROWS ) )
+ ALLOCATE( CLR_TRANS_DIRECT( NCOLS,NROWS ) )
+ ALLOCATE( CLR_REFLECTION ( NCOLS,NROWS ) )
+ ALLOCATE( ETOT_SFC_WL ( NCOLS,NROWS,NWL ) )
+ ALLOCATE( TAU_AERO_WL ( NCOLS,NROWS,NWL ) )
+ ALLOCATE( TAU_CLOUD_WL ( NCOLS,NROWS,NWL ) )
+#ifdef phot_debug
+ ALLOCATE( SSA_CLOUD_WL( NCOLS,NROWS,NWL ) )
+ ALLOCATE( ASY_CLOUD_WL( NCOLS,NROWS,NWL ) )
+#endif
+ ALLOCATE( TAU_TOT_WL ( NCOLS,NROWS,NWL ) )
+ ALLOCATE( TAUO3_TOP_WL( NCOLS,NROWS,NWL ) )
+
+ N_EXCEED_TROPO3 = 0.0
+ TROPO_O3_EXCEED = 0.0
+ TSTEP_COUNT = 0
+
+ DIAG_WVL( 1 ) = 1
+ DIAG_WVL( N_DIAG_WVL ) = NWL
+
+ ALLOCATE ( AERO_ASYM( NCOLS,NROWS,NLAYS,N_DIAG_WVL ), STAT = ALLOCSTAT )
+ IF ( ALLOCSTAT .NE. 0 ) THEN
+ XMSG = 'Failure allocating 3D AERO_ASYM'
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ END IF
+
+ ALLOCATE ( AERO_SSA( NCOLS,NROWS,NLAYS,N_DIAG_WVL ), STAT = ALLOCSTAT )
+ IF ( ALLOCSTAT .NE. 0 ) THEN
+ XMSG = 'Failure allocating 3D AERO_SSA'
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ END IF
+
+ ALLOCATE ( TAU_AERO( NCOLS,NROWS,NLAYS,N_DIAG_WVL ), STAT = ALLOCSTAT )
+ IF ( ALLOCSTAT .NE. 0 ) THEN
+ XMSG = 'Failure allocating 3D TAU_AERO'
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ END IF
+
+ ALLOCATE ( TAU( NCOLS,NROWS,NLAYS,N_DIAG_WVL ), STAT = ALLOCSTAT )
+ IF ( ALLOCSTAT .NE. 0 ) THEN
+ XMSG = 'Failure allocating 3D TAU'
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ END IF
+
+ ALLOCATE ( ACTINIC_FX( NCOLS,NROWS,NLAYS,NWL ), STAT = ALLOCSTAT )
+ IF ( ALLOCSTAT .NE. 0 ) THEN
+ XMSG = 'Failure allocating ACTINIC_FX'
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ END IF
+
+!...write wavelength data to a character array
+
+ ALLOCATE ( WLTXT( NWL ) )
+
+ DO IWL = 1, NWL
+ WRITE( WLTXT( IWL ),'(I3.3)' ) INT( WAVELENGTH( IWL ) )
+ END DO
+
+!...open the photolysis rate diagnostic files
+
+ ODATE = JDATE; OTIME = JTIME
+#ifndef phot_extra_tstep
+ CALL NEXTIME ( ODATE, OTIME, DTSTEP( 1 ) ) ! output timestamp ending time
+#endif
+ IF ( IO_PE_INCLUSIVE ) CALL OPPHOT ( ODATE, OTIME, DTSTEP( 1 ) )
+
+ CALL SUBST_BARRIER
+
+ END IF ! photdiag
+
+!...set pointers to species O3 and NO2 in CGRID
+
+ VARNM = 'O3'
+ LGC_O3 = INDEX1( VARNM, N_GC_SPC, GC_SPC )
+ IF ( LGC_O3 .LE. 0 ) THEN
+ XMSG = 'Could not find ' // VARNM // 'in species table'
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 )
+ END IF
+
+ VARNM = 'NO2'
+ LGC_NO2 = INDEX1( VARNM, N_GC_SPC, GC_SPC )
+ IF ( LGC_NO2 .LE. 0 ) THEN
+ XMSG = 'Could not find ' // VARNM // 'in species table'
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT3 )
+ END IF
+
+#ifdef phot_extra_tstep
+ ELSE
+ IF ( PHOTDIAG ) THEN
+ ODATE = JDATE; OTIME = JTIME
+ CALL NEXTIME ( ODATE, OTIME, DTSTEP( 2 ) ) ! sync time step
+ END IF
+#endif
+ END IF ! firstime
+
+ IF ( JD_STRAT_O3MIN .NE. JDATE ) THEN
+!...set minimum fraction of ozone column above PTOP
+
+ CALL SEASONAL_STRAT_O3( JDATE, JTIME )
+ MIN_STRATO3_FRAC = MONTH_STRAT_03_FRAC
+ MAX_TROPOO3_FRAC = MAX( 1.0 - MONTH_STRAT_03_FRAC, 0.0 )
+ WRITE( LOGDEV,*)'PHOT: MIN_STRATO3_FRAC = ',MIN_STRATO3_FRAC
+
+ JD_STRAT_O3MIN = REAL( JDATE, 4)
+ END IF
+!...initialize variables tracking whether stratosphere ozone column satisfies
+!...climatological averages.
+
+ O3_TOGGLE_AVE = 0.0
+ O3_TOGGLE_MIN = 1.0
+ N_TROPO_O3_TOGGLE = 0
+ TSTEP_COUNT = TSTEP_COUNT + 1
+
+ CALL GET_PHOT_MET( JDATE, JTIME, MDATE, MTIME )
+
+!...Get cosine of solar parameters and set DARK
+
+ CALL UPDATE_SUN( JDATE, JTIME, MDATE, MTIME )
+
+ RSQD = DIST_TO_SUN * DIST_TO_SUN
+
+ IF ( MAXVAL( COSINE_ZENITH ) .LE. 0.0 ) THEN
+ DARK = .TRUE.
+ ELSE
+ DARK = .FALSE.
+ END IF
+
+!...set surface albedos
+
+ CALL GET_ALBEDO( MDATE, MTIME, LOGDEV, COSINE_ZENITH, LAT, LON )
+
+!...SA Write COSINE_ZENITH array at the end of each output tstep
+
+ IF ( PHOTDIAG ) THEN
+#ifndef phot_extra_tstep
+ ODATE = JDATE; OTIME = JTIME
+ CALL NEXTIME ( ODATE, OTIME, DTSTEP( 2 ) ) ! sync time step
+#endif
+ JTIME_CHK = ( MOD( TIME2SEC( OTIME ), TSTEP ) .EQ. 0 )
+#ifdef parallel_io
+ IF ( .NOT. IO_PE_INCLUSIVE ) THEN
+ IF ( .NOT. OPEN3( CTM_RJ_1, FSREAD3, PNAME ) ) THEN
+ XMSG = 'Could not open ' // TRIM(CTM_RJ_1)
+ CALL M3EXIT( PNAME, ODATE, OTIME, XMSG, XSTAT1 )
+ END IF
+ IF ( .NOT. OPEN3( CTM_RJ_2, FSREAD3, PNAME ) ) THEN
+ XMSG = 'Could not open ' // TRIM(CTM_RJ_2)
+ CALL M3EXIT( PNAME, ODATE, OTIME, XMSG, XSTAT1 )
+ END IF
+ END IF
+#endif
+ ELSE
+ JTIME_CHK = .FALSE.
+ END IF
+
+
+!...If sun below horizon at all cells, zero photolysis rates & exit
+!... (assumes sun below horizon at *all* levels!)
+
+ IF ( DARK ) THEN
+
+ RJ = 0.0
+
+!...write to the log file, CTM_RJ_1 file and return
+
+ WRITE( LOGDEV, 1003 ) MYPE, JDATE, JTIME
+
+!...Initialize ETOT_SFC, TAU_AERO, TAU_TOT, TAUO3_TOP to 0.0
+
+!...Write data to output diagnostic file
+
+ TOTAL_OC = 0.0
+
+ IF ( JTIME_CHK ) THEN
+
+ TROPO_OC = 0.0
+ ETOT_SFC_WL = 0.0
+ TAUO3_TOP_WL = 0.0
+ TAU_AERO_WL = 0.0
+ TAU_CLOUD_WL = 0.0
+#ifdef phot_debug
+ SSA_CLOUD_WL = 0.0
+ ASY_CLOUD_WL = 0.0
+#endif
+ TAU_TOT_WL = 0.0
+ TAU = 0.0
+ TAU_AERO = 0.0
+ AERO_SSA = 0.0
+ AERO_ASYM = 0.0
+ ACTINIC_FX = 0.0
+
+! TROPO_O3_EXCEED = 0.0
+ TRANSMIS_DIFFUSE = 0.0
+ TRANSMIS_DIRECT = 0.0
+ REFLECT_COEFF = 0.0
+ CLR_TRANSMISSION = 0.0
+ CLR_TRANS_DIRECT = 0.0
+ CLR_REFLECTION = 0.0
+
+ END IF ! if JTIME_CHK
+
+ ELSE ! all cells not dark
+
+!...MAIN loop over all rows and columns
+ LOOP_ROWS: DO ROW = 1, MY_NROWS
+ LOOP_COLS: DO COL = 1, MY_NCOLS
+
+ PHOT_COL = COL + PECOL_OFFSET
+ PHOT_ROW = ROW + PEROW_OFFSET
+
+ COSZEN = COSINE_ZENITH( COL,ROW ) ! local cosine of solar zenith angle
+
+ IF ( COSZEN .LE. 0.0 ) THEN
+!...the cell is dark: set variables to zero and cycle
+ RJ( COL,ROW, :,: ) = 0.0
+
+ IF ( JTIME_CHK ) THEN
+ TOTAL_OC( COL,ROW ) = 0.0
+ TROPO_OC( COL,ROW ) = 0.0
+ ETOT_SFC_WL ( COL,ROW, : ) = 0.0
+ TAUO3_TOP_WL( COL,ROW, : ) = 0.0
+ TAU_AERO_WL ( COL,ROW, : ) = 0.0
+ TAU_CLOUD_WL( COL,ROW, : ) = 0.0
+#ifdef phot_debug
+ SSA_CLOUD_WL( COL,ROW, : ) = 0.0
+ ASY_CLOUD_WL( COL,ROW, : ) = 0.0
+#endif
+ TAU_TOT_WL( COL,ROW, : ) = 0.0
+ TAU ( COL,ROW, :,: ) = 0.0
+ TAU_AERO ( COL,ROW, :,: ) = 0.0
+ AERO_SSA ( COL,ROW, :,: ) = 0.0
+ AERO_ASYM ( COL,ROW, :,: ) = 0.0
+ ACTINIC_FX( COL,ROW, :,: ) = 0.0
+
+! TROPO_O3_EXCEED( COL,ROW ) = 0.0
+ TRANSMIS_DIFFUSE( COL,ROW ) = 0.0
+ TRANSMIS_DIRECT ( COL,ROW ) = 0.0
+ REFLECT_COEFF ( COL,ROW ) = 0.0
+ CLR_TRANSMISSION( COL,ROW ) = 0.0
+ CLR_TRANS_DIRECT( COL,ROW ) = 0.0
+ CLR_REFLECTION ( COL,ROW ) = 0.0
+ END IF
+
+ CYCLE LOOP_COLS
+
+ END IF
+
+!...initialize BLKRJ using F90 array operations.
+
+ BLKRJ_RES = 0.0
+ BLKRJ_ACM = 0.0
+
+!...Set height of lowest level to zero
+
+ BLKZF( 1 ) = 0.0
+
+ ZSFC = HT( COL,ROW ) ! surface height [m]
+ SINZEN = SQRT( 1.0 - COSZEN * COSZEN ) ! sine of zenith angle
+
+!...local latitude and longitude
+
+! LATCR = LAT( COL,ROW )
+! LONCR = LON( COL,ROW )
+
+!...get total ozone column based on OMI observations
+ CALL O3TOTCOL ( LAT( COL,ROW ), LON( COL,ROW ), JDATE, TOTAL_O3_COLUMN )
+
+ IF ( USE_ACM_CLOUD .OR. CLDATT ) THEN
+ OWATER_FRAC = MAX( ( 1.0 - SEAICE( COL,ROW ) ), 0.0 )
+ & * WATER_FRACTION( COL,ROW )
+ SEAICE_FRAC = SEAICE( COL,ROW ) * WATER_FRACTION( COL,ROW )
+ SNOW_FRAC = SNOCOV( COL,ROW )
+ COL_CLOUD = PHOT_COL
+ ROW_CLOUD = PHOT_ROW
+ END IF
+
+!...loop over vertical layers ambient air conditions and gas concentration
+ DO L = 1, NLAYS
+!...Fetch the grid cell ambient data at each layer.
+
+ BLKTA ( L ) = TA ( COL,ROW,L ) ! temperature [K]
+ BLKPRS ( L ) = PRES ( COL,ROW,L ) / STDATMPA ! [atmospheres]
+ BLKDENS( L ) = DENS ( COL,ROW,L ) * DENS_CONV ! [molecules / cm**3]
+ BLKZH ( L ) = ZM ( COL,ROW,L ) ! mid layer height [m]
+ BLKZF ( L+1 ) = ZFULL( COL,ROW,L ) ! full layer height [m]
+
+!...set scale factor for [ppm] -> [molecule / cm**3]
+!... To go from ppm to molecule/cc:
+!... molecule/cc = ppm * 1.0E-06 * DENS (given in molecule/cc)
+
+ MSCALE = BLKDENS( L ) * PPM_MCM3
+
+!...fetch ozone and no2 and convert to [ molecules / cm **3 ]
+!... and adjust the volume for ambient temperature and pressure.
+
+ BLKO3 ( L ) = CGRID( COL,ROW,L,LGC_O3 ) * MSCALE
+ BLKNO2( L ) = CGRID( COL,ROW,L,LGC_NO2 ) * MSCALE
+ ZLEV = BLKZF( L )
+ END DO ! loop on layers ambient conditions and gases
+
+ IF ( CLDATT .AND. CFRAC_2D( COL,ROW ) .GT. 0.0 ) THEN
+ DO L = 1, NLAYS
+
+ IF ( CFRAC_3D( COL,ROW,L ) .GT. 0.0 ) THEN
+ CLOUDS ( L ) = .TRUE.
+ CLOUD_LAYERING( L ) = .TRUE.
+ CLDFRAC( L ) = CFRAC_3D( COL,ROW,L )
+!... set hydrometeor concentrations for resolved cloud
+ MSCALE = 1.0E+3 * DENS ( COL,ROW,L )
+ IWC( L ) = MSCALE * QI( COL,ROW,L )
+ GWC( L ) = MSCALE * QG( COL,ROW,L )
+ SWC( L ) = MSCALE * QS( COL,ROW,L )
+ LWC( L ) = MSCALE * QC( COL,ROW,L )
+ RWC( L ) = MSCALE * QR( COL,ROW,L )
+ ELSE
+ CLOUDS ( L ) = .FALSE.
+ CLOUD_LAYERING( L ) = .FALSE.
+ CLDFRAC( L ) = 0.0
+ IWC( L ) = 0.0
+ GWC( L ) = 0.0
+ SWC( L ) = 0.0
+ LWC( L ) = 0.0
+ RWC( L ) = 0.0
+ END IF
+ END DO ! loop on layers clouds
+! get optical properties of resolved cloud hydrometeors
+ CALL GET_DROPLET_OPTICS( NLAYS, BLKTA, OWATER_FRAC, SEAICE_FRAC, SNOW_FRAC, LWC )
+ CALL GET_ICE_OPTICS( NLAYS, BLKTA, IWC )
+ CALL GET_AGGREGATE_OPTICS( NLAYS, RWC, SWC, GWC )
+ ELSE
+ CLOUDS = .FALSE.
+ CLOUD_LAYERING = .FALSE.
+ CLDFRAC = 0.0
+! hydrometeor concentrations
+ LWC = 0.0
+ IWC = 0.0
+ RWC = 0.0
+ SWC = 0.0
+ RWC = 0.0
+ CALL CLEAR_HYDROMETEOR_OPTICS()
+ END IF
+
+!..calculate needed aerosol properties in column
+
+! IF ( CORE_SHELL ) THEN
+ CALL GET_AERO_DATA ( COL,ROW, NLAYS, CGRID )
+! ELSE
+! CALL AERO_OPTICS_INTERNAL( COL,ROW, NLAYS, CGRID )
+! END IF
+
+! set surface albedo
+
+ FORALL ( IWL = 1:NWL )
+ ALB( IWL ) = SURFACE_ALBEDO( IWL, COL,ROW )
+ END FORALL
+!set min/max fractions of ozone column in stratosphere and troposphere
+! MIN_STRATO3_FRAC = MIN_STRAT_03_FRAC( COL, ROW )
+! MAX_TROPOO3_FRAC = MAX( 1.0 - MIN_STRAT_03_FRAC( COL, ROW ), 0.0 )
+! MIN_STRATO3_FRAC = MONTH_STRAT_03_FRAC( COL, ROW )
+! MAX_TROPOO3_FRAC = MAX( 1.0 - MONTH_STRAT_03_FRAC( COL, ROW ), 0.0 )
+!...calculate resolved-sky photolysis rates at all layers:
+
+ NEW_PROFILE = .TRUE.
+ ONLY_SOLVE_RAD = .FALSE.
+
+ CALL NEW_OPTICS ( JDATE, JTIME, NLAYS,
+ & BLKTA, BLKPRS, BLKDENS, BLKZH, BLKZF,
+ & BLKO3, BLKNO2,
+ & ZSFC, COSZEN, SINZEN, RSQD,
+ & NEW_PROFILE, CLOUDS, CLDFRAC,
+ & BLKRJ_RES, TAUC_AERO, TAU_TOT, TAUO3_TOP,
+ & TAU_RAY, SSA, TAU_CLOUD, TOTAL_O3_COLUMN )
+
+!...load diagnostic file arrays
+ IF ( PHOTDIAG .AND. .NOT. STRATO3_MINS_MET ) THEN
+ N_EXCEED_TROPO3( COL,ROW ) = N_EXCEED_TROPO3( COL,ROW ) + 1.0
+ TROPO_O3_EXCEED( COL,ROW ) = TROPO_O3_COLUMN/(MAX_TROPOO3_FRAC*TOTAL_O3_COLUMN) - 1.0
+! & + 1.0 / TROPO_O3_TOGGLE - 1.0
+ & + TROPO_O3_EXCEED( COL,ROW )
+! ELSE IF( PHOTDIAG ) THEN
+! TROPO_O3_EXCEED( COL,ROW ) = 0.0
+ END IF
+
+ IF ( JTIME_CHK ) THEN
+ TOTAL_OC( COL,ROW ) = REAL( TOTAL_O3_COLUMN )
+ TROPO_OC( COL,ROW ) = REAL( TROPO_O3_COLUMN )
+ TRANSMIS_DIFFUSE( COL,ROW ) = TRANSMISSION
+ TRANSMIS_DIRECT( COL,ROW ) = TRANS_DIRECT
+ REFLECT_COEFF( COL,ROW ) = REFLECTION
+
+
+ FORALL( IWL = 1:NWL )
+ ETOT_SFC_WL ( COL,ROW,IWL ) = IRRADIANCE( 1,IWL )
+ TAUO3_TOP_WL( COL,ROW,IWL ) = TAUO3_TOP( IWL )
+ TAU_AERO_WL ( COL,ROW,IWL ) = TAUC_AERO( 1,IWL )
+ TAU_TOT_WL ( COL,ROW,IWL ) = TAU_TOT ( 1,IWL )
+ TAU_CLOUD_WL( COL,ROW,IWL ) = TAU_CLOUD( 1,IWL )
+#ifdef phot_debug
+ SSA_CLOUD_WL( COL,ROW,IWL ) = AVE_SSA_CLD ( IWL )
+ ASY_CLOUD_WL( COL,ROW,IWL ) = AVE_ASYMM_CLD( IWL )
+#endif
+ END FORALL
+ FORALL ( LEV = 1:NLAYS, IWL = 1:NWL )
+ ACTINIC_FX( COL,ROW,LEV,IWL ) = ACTINIC_FLUX( LEV,IWL )
+ END FORALL
+
+ DO L = 1, N_DIAG_WVL
+ IWL = DIAG_WVL( L )
+ FORALL ( LEV = 1:NLAYS )
+ TAU ( COL,ROW,LEV,L ) = TAU_TOT ( LEV,IWL )
+ TAU_AERO( COL,ROW,LEV,L ) = TAUC_AERO( LEV,IWL )
+ END FORALL
+ FORALL ( LEV = 1:NLAYS, AERO_EXTI_COEF( LEV,IWL ) .GT. EPSLON )
+ AERO_SSA ( COL,ROW,LEV,L ) = AERO_SCAT_COEF( LEV,IWL )
+ & / AERO_EXTI_COEF( LEV,IWL )
+ AERO_ASYM( COL,ROW,LEV,L ) = AERO_ASYM_FAC( LEV,IWL )
+ END FORALL
+ FORALL ( LEV = 1:NLAYS, AERO_EXTI_COEF( LEV,IWL ) .LE. EPSLON )
+ AERO_SSA ( COL,ROW,LEV,L ) = 1.0
+ AERO_ASYM( COL,ROW,LEV,L ) = 0.0
+ END FORALL
+ END DO
+ END IF
+
+!Set Photolysis rates to resolved sky values
+ FORALL ( L = 1:NLAYS, IPHOT = 1:NPHOTAB )
+ RJ( COL,ROW, L,IPHOT ) = 60.0 * BLKRJ_RES( L,IPHOT )
+ END FORALL ! Loop on layers and NPHOTAB
+
+ IF ( USE_ACM_CLOUD ) THEN
+ IF ( ACM_CLOUDS( COL,ROW ) .GT. 0.0 ) THEN
+! save resolved sky reflection and transmission coefficients for possible latter use
+ RES_SKY_REFLECT = REFLECTION
+ RES_SKY_TRANS = TRANSMISSION
+ RES_SKY_TRANSD = TRANS_DIRECT
+!...find the highest layer of the sub-grid (convective) cloud
+ DO LEV = NLAYS, 1, -1
+ IF ( ACM_CFRAC( LEV, COL,ROW ) .GT. 0.0 ) EXIT
+ END DO
+!...replace the lower layers with sub-grid cloud properties
+ DO L = 1, LEV
+ SWC( L ) = 0.0
+ IF ( ACM_CFRAC( L,COL,ROW ) .GT. 0.0 ) THEN
+ CLOUDS ( L ) = .TRUE.
+ CLDFRAC( L ) = 1.0
+ MSCALE = 1.0E+3 * DENS ( COL,ROW, L )
+ LWC( L ) = MSCALE * ACM_QC( L,COL,ROW )
+ IWC( L ) = MSCALE * ACM_QI( L,COL,ROW )
+ RWC( L ) = MSCALE * ACM_QR( L,COL,ROW )
+ GWC( L ) = MSCALE * ACM_QG( L,COL,ROW )
+ ELSE
+ CLOUDS( L ) = .FALSE.
+ CLDFRAC( L ) = 0.0
+ LWC( L ) = 0.0
+ IWC( L ) = 0.0
+ RWC( L ) = 0.0
+ GWC( L ) = 0.0
+ END IF
+ CLOUD_LAYERING( L ) = .FALSE.
+ END DO
+! write(logdev,*)'ACM cloud present fraction, cloud lwc(lev),iwc(lev),rwc(1),gwc(1) = ',
+! & ACM_CLOUDS( COL,ROW ),lwc(lev),iwc(lev),rwc(1),gwc(1)
+
+! get optical properties of of subgrid cloud hydrometeors
+ CALL GET_DROPLET_OPTICS( LEV, BLKTA, OWATER_FRAC, SEAICE_FRAC, SNOW_FRAC, LWC )
+ CALL GET_ICE_OPTICS( LEV, BLKTA, IWC )
+ CALL GET_AGGREGATE_OPTICS( LEV, RWC, SWC, GWC )
+
+!...calculate the acm-cloud photolysis rates for all layers:
+ NEW_PROFILE = .FALSE.
+ CALL NEW_OPTICS ( JDATE, JTIME, NLAYS,
+ & BLKTA, BLKPRS, BLKDENS, BLKZH, BLKZF,
+ & BLKO3, BLKNO2,
+ & ZSFC, COSZEN, SINZEN, RSQD,
+ & NEW_PROFILE, CLOUDS, CLDFRAC,
+ & BLKRJ_ACM, TAUC_AERO, TAU_TOT, TAUO3_TOP,
+ & TAU_RAY, SSA, TAU_CLOUD, TOTAL_O3_COLUMN )
+
+!...load diagnostic file arrays
+!...compute a cloud-fraction weighted average of ETOT_SFC and TAU_TOT
+!... note that both TAUC_AERO and TAUO3_TOP are the same for clear and
+!... cloudy regions
+ MSCALE = MAX( 1.0 - ACM_CLOUDS( COL,ROW ), 0.0 )
+
+ IF ( JTIME_CHK ) THEN
+
+ TRANSMIS_DIRECT( COL,ROW ) = MSCALE * TRANSMIS_DIRECT( COL,ROW )
+ & + ACM_CLOUDS( COL,ROW ) * TRANS_DIRECT
+ TRANSMIS_DIFFUSE( COL,ROW ) = MSCALE * TRANSMIS_DIFFUSE( COL,ROW )
+ & + ACM_CLOUDS( COL,ROW ) * TRANSMISSION
+ REFLECT_COEFF( COL,ROW ) = MSCALE * REFLECT_COEFF( COL,ROW )
+ & + ACM_CLOUDS( COL,ROW ) * REFLECTION
+ FORALL ( IWL = 1:NWL )
+ ETOT_SFC_WL ( COL,ROW,IWL ) = MSCALE * ETOT_SFC_WL( COL,ROW,IWL )
+ & + ACM_CLOUDS( COL,ROW ) * IRRADIANCE( 1,IWL )
+ TAU_TOT_WL ( COL,ROW,IWL ) = MSCALE * TAU_TOT_WL( COL,ROW,IWL )
+ & + ACM_CLOUDS( COL,ROW ) * TAU_TOT( 1,IWL )
+ TAU_CLOUD_WL( COL,ROW,IWL ) = MSCALE * TAU_CLOUD_WL( COL,ROW,IWL )
+ & + ACM_CLOUDS( COL,ROW ) * TAU_CLOUD( 1,IWL )
+#ifdef phot_debug
+ SSA_CLOUD_WL( COL,ROW,IWL ) = MSCALE * SSA_CLOUD_WL( COL,ROW,IWL )
+ & + ACM_CLOUDS( COL,ROW ) * AVE_SSA_CLD ( IWL )
+ ASY_CLOUD_WL( COL,ROW,IWL ) = MSCALE * ASY_CLOUD_WL( COL,ROW,IWL )
+ & + ACM_CLOUDS( COL,ROW ) * AVE_ASYMM_CLD( IWL )
+#endif
+ END FORALL ! iwl
+ FORALL ( LEV = 1:NLAYS, IWL = 1:NWL )
+ ACTINIC_FX( COL,ROW,LEV,IWL ) = MSCALE * ACTINIC_FX( COL,ROW,LEV,IWL )
+ & + ACM_CLOUDS( COL,ROW ) * ACTINIC_FLUX( LEV,IWL )
+ END FORALL ! lev and iwl
+
+ DO L = 1, N_DIAG_WVL
+ IWL = DIAG_WVL( L )
+ FORALL ( LEV = 1:NLAYS)
+ TAU( COL,ROW,LEV,L ) = MSCALE * TAU( COL,ROW,LEV,L )
+ & + ACM_CLOUDS( COL,ROW ) * TAU_TOT( LEV,IWL )
+ END FORALL
+ END DO
+ END IF ! photdiag
+!Photolysis rates become a weighted average of the values from resolved and ACM skies
+ FORALL ( L = 1:NLAYS, IPHOT = 1:NPHOTAB )
+ RJ( COL,ROW, L, IPHOT ) = 60.0 * ACM_CLOUDS( COL,ROW ) * BLKRJ_ACM( L,IPHOT )
+ & + MSCALE * RJ( COL,ROW,L,IPHOT )
+ END FORALL ! Loop on layers and PHOT
+ END IF
+ END IF ! not USE_ACM_CLOUD and ACM_CLOUDS > 0
+
+!------------------------CANOPY PHOTOLYSIS CORRECTION/REDUCTION Section NOAA-ARL-------------------------------------------
+!Conditions to reduce weighted average of photolysis rates (RJ) due to canopy shading (if user-defined=true); P. C. Campbell
+!Following is based on work of ECCC in GEM-MACHv2.1: Makar et al. (2017)
+!Makar, P., Staebler, R., Akingunola, A. et al. The effects of forest canopy shading and turbulence on boundary layer ozone.
+!Nat Commun 8, 15243 (2017). https://doi.org/10.1038/ncomms15243
+
+ IF ( CANOPY_SHADE ) THEN ! compute canopy shade reduction factor (RJ_CORR)
+
+ DO LEV = 1, NLAYS !loop through model layers
+
+ IF (LEV .EQ. 1) THEN !first model layer
+ KCAN = 1
+ ELSE !check subsequent model layers
+ IF ( Met_Data%FCH( COL,ROW ) .GT. Met_Data%ZF( COL,ROW,LEV-1 )
+ & .AND. Met_Data%FCH( COL,ROW ) .LE. Met_Data%ZF( COL,ROW,LEV ) ) THEN
+ KCAN = 1
+ ELSE
+ KCAN = 0
+ END IF
+ END IF
+
+ IF (KCAN .EQ. 1) THEN !canopy could be inside model layer
+ !check for otherconditions for grid cells that do NOT have
+ !a continuos forest canopy
+ IF ( Met_Data%LAIE( COL,ROW ) .LT. 0.1
+ & .OR. Met_Data%FCH( COL,ROW ) .LT. 0.5
+! & .OR. MAX(0.0, 1.0 - Met_Data%FRT( COL,ROW)) .GT. 0.5 !IVAI: old canopy data
+ & .OR. MAX(0.0, 1.0 - Met_Data%FRT( COL,ROW)) .GT. 0.75 !IVAI: new canopy data
+ & .OR. Met_Data%POPU( COL,ROW ) .GT.10000.0
+ & .OR. (EXP(-0.5*Met_Data%LAIE( COL,ROW)*Met_Data%CLU( COL,ROW )) .GT. 0.45
+ & .AND. Met_Data%FCH(COL,ROW ) .LT. 18.0 )) THEN
+ RJ( COL,ROW, LEV, : ) = RJ( COL,ROW, LEV, :)
+ ELSE ! There is a contiguous forest canopy,apply correctoin
+ !RJ_CORR effectly represents the beam attenuation and reduces photolysis.
+ !Nilson, T. A theoretical analysis of the frequency of gaps in plant stands. Agric.
+ !Meterol. 8, 25⚌~Z~L~@~S38 (1971).
+
+!Calculate attenuation at different set cumulative LAI fractions downward through canopy (C1R, C2R, C3R, C4R data from ECCC)
+ RJ_CORR_C1R( COL,ROW ) = MAX(1.0E-10, EXP(-1.0*(0.5*(Met_Data%LAIE( COL,ROW )
+ & *Met_Data%C1R( COL,ROW ))*Met_Data%CLU( COL,ROW ))/MAX(0.05, COSZEN)))
+ RJ_CORR_C2R( COL,ROW ) = MAX(1.0E-10, EXP(-1.0*(0.5*(Met_Data%LAIE( COL,ROW )
+ & *Met_Data%C2R( COL,ROW ))*Met_Data%CLU( COL,ROW ))/MAX(0.05, COSZEN)))
+ RJ_CORR_C3R( COL,ROW ) = MAX(1.0E-10, EXP(-1.0*(0.5*(Met_Data%LAIE( COL,ROW )
+ & *Met_Data%C3R( COL,ROW ))*Met_Data%CLU( COL,ROW ))/MAX(0.05, COSZEN)))
+ RJ_CORR_C4R( COL,ROW ) = MAX(1.0E-10, EXP(-1.0*(0.5*(Met_Data%LAIE( COL,ROW )
+ & *Met_Data%C4R( COL,ROW ))*Met_Data%CLU( COL,ROW ))/MAX(0.05, COSZEN)))
+ RJ_CORR_BOT( COL,ROW ) = MAX(1.0E-10, EXP(-1.0*(0.5*Met_Data%LAIE( COL,ROW )
+ & *Met_Data%CLU( COL,ROW ))/MAX(0.05, COSZEN)))
+
+!Interpolate to get attenuation profile below canopy inside respective model layer
+ ZFL = Met_Data%ZF( COL,ROW,LEV )
+ ZCAN = ZFL ! Initialize top (m) = Bottom of model layer
+ COUNTCAN = 0 ! Initialize canopy layers
+
+ IF (LEV .EQ. 1) THEN !Find bottom in each model layer
+ BOTCAN = 0.5
+ ELSE
+ BOTCAN = Met_Data%ZF( COL,ROW,LEV-1 )
+ END IF
+
+ DO WHILE (ZCAN.GE.BOTCAN)
+ IF ( ZCAN .GT. Met_Data%FCH( COL,ROW ) ) THEN
+ COUNTCAN = COUNTCAN + 1
+ ZCANX(COUNTCAN) = ZCAN
+ RJ_CORRX (COUNTCAN) = 1.0
+ ELSE IF ( ZCAN .LE. Met_Data%FCH( COL,ROW ) .AND.
+ & ZCAN .GT. Met_Data%FCH( COL,ROW )*0.75 ) THEN
+ COUNTCAN = COUNTCAN + 1
+ XCAN(2) = Met_Data%FCH( COL,ROW )
+ YCAN(2) = 1.0
+ XCAN(1) = Met_Data%FCH( COL,ROW )*0.75
+ YCAN(1) = RJ_CORR_C1R( COL,ROW )
+ XCANOUT = ZCAN
+ ZCANX(COUNTCAN) = ZCAN
+ RJ_CORRX (COUNTCAN) = interp_linear1_internal(XCAN,YCAN,XCANOUT)
+ ELSE IF ( ZCAN .LE. Met_Data%FCH( COL,ROW )*0.75 .AND.
+ & ZCAN .GT. Met_Data%FCH( COL,ROW )*0.50 ) THEN
+ COUNTCAN = COUNTCAN + 1
+ XCAN(2) = Met_Data%FCH( COL,ROW )*0.75
+ YCAN(2) = RJ_CORR_C1R( COL,ROW )
+ XCAN(1) = Met_Data%FCH( COL,ROW )*0.50
+ YCAN(1) = RJ_CORR_C2R( COL,ROW )
+ XCANOUT = ZCAN
+ ZCANX(COUNTCAN) = ZCAN
+ RJ_CORRX (COUNTCAN) = interp_linear1_internal(XCAN,YCAN,XCANOUT)
+ ELSE IF ( ZCAN .LE. Met_Data%FCH( COL,ROW )*0.50 .AND.
+ & ZCAN .GT. Met_Data%FCH( COL,ROW )*0.35 ) THEN
+ COUNTCAN = COUNTCAN + 1
+ XCAN(2) = Met_Data%FCH( COL,ROW )*0.50
+ YCAN(2) = RJ_CORR_C2R( COL,ROW )
+ XCAN(1) = Met_Data%FCH( COL,ROW )*0.35
+ YCAN(1) = RJ_CORR_C3R( COL,ROW )
+ XCANOUT = ZCAN
+ ZCANX(COUNTCAN) = ZCAN
+ RJ_CORRX (COUNTCAN) = interp_linear1_internal(XCAN,YCAN,XCANOUT)
+ ELSE IF ( ZCAN .LE. Met_Data%FCH( COL,ROW )*0.35 .AND.
+ & ZCAN .GT. Met_Data%FCH( COL,ROW )*0.20 ) THEN
+ COUNTCAN = COUNTCAN + 1
+ XCAN(2) = Met_Data%FCH( COL,ROW )*0.35
+ YCAN(2) = RJ_CORR_C3R( COL,ROW )
+ XCAN(1) = Met_Data%FCH( COL,ROW )*0.20
+ YCAN(1) = RJ_CORR_C4R( COL,ROW )
+ XCANOUT = ZCAN
+ ZCANX(COUNTCAN) = ZCAN
+ RJ_CORRX (COUNTCAN) = interp_linear1_internal(XCAN,YCAN,XCANOUT)
+ ELSE IF ( ZCAN .LE. Met_Data%FCH( COL,ROW )*0.20 ) THEN
+ COUNTCAN = COUNTCAN + 1
+ XCAN(2) = Met_Data%FCH( COL,ROW )*0.20
+ YCAN(2) = RJ_CORR_C4R( COL,ROW )
+ XCAN(1) = 0.5
+ YCAN(1) = RJ_CORR_BOT( COL,ROW )
+ XCANOUT = ZCAN
+ ZCANX(COUNTCAN) = ZCAN
+ RJ_CORRX (COUNTCAN) = interp_linear1_internal(XCAN,YCAN,XCANOUT)
+ END IF
+ ZCAN = ZCAN-0.5 !step down in-canopy resolution of 0.5 m
+ END DO !end loop on canopy layers
+
+!Integrate to get best attenuation value to use within canopy
+ RJ_CORR( COL,ROW ) = IntegrateTrapezoid(ZCANX(COUNTCAN:1:-1),RJ_CORRX(COUNTCAN:1:-1)) /
+ & ZFL
+!Apply attenuation factors above and below canopy
+ RJ( COL,ROW, LEV, : ) = RJ( COL,ROW, LEV, : )*RJ_CORR( COL,ROW )
+!Apply attenuation value within canopy and take average above and within canopy values
+ END IF !other contiguous canopy conditions
+ END IF !canopy height could be inside model layer
+ END DO !end loop on model layers
+ END IF !canopy shade
+
+ IF ( JTIME_CHK ) THEN ! compute clear sky reflection and transmission coefficients
+ IF ( ANY( CLOUDS ) ) THEN
+ IF ( CFRAC_2D( COL,ROW ) .GT. 0.0 ) THEN ! resolved and subgrid clouds exist
+ CLOUDS = .FALSE.
+ NEW_PROFILE = .FALSE.
+ ONLY_SOLVE_RAD = .TRUE.
+ CALL NEW_OPTICS ( JDATE, JTIME, NLAYS,
+ & BLKTA, BLKPRS, BLKDENS, BLKZH, BLKZF,
+ & BLKO3, BLKNO2,
+ & ZSFC, COSZEN, SINZEN, RSQD,
+ & NEW_PROFILE, CLOUDS, CLDFRAC,
+ & BLKRJ_RES, TAUC_AERO, TAU_TOT, TAUO3_TOP,
+ & TAU_RAY, SSA, TAU_CLOUD, TOTAL_O3_COLUMN)
+ CLR_REFLECTION ( COL,ROW ) = REFLECTION
+ CLR_TRANSMISSION( COL,ROW ) = TRANSMISSION
+ CLR_TRANS_DIRECT( COL,ROW ) = TRANS_DIRECT
+ ELSE ! only subgrid clouds exist
+ CLR_REFLECTION ( COL,ROW ) = RES_SKY_REFLECT
+ CLR_TRANSMISSION( COL,ROW ) = RES_SKY_TRANS
+ CLR_TRANS_DIRECT( COL,ROW ) = RES_SKY_TRANSD
+ END IF
+ ELSE ! no cloud in vertical column
+ CLR_REFLECTION ( COL,ROW ) = REFLECTION
+ CLR_TRANSMISSION( COL,ROW ) = TRANSMISSION
+ CLR_TRANS_DIRECT( COL,ROW ) = TRANS_DIRECT
+ END IF
+ END IF
+
+ END DO LOOP_COLS
+ END DO LOOP_ROWS
+
+ END IF
+
+!...report on whether stratospheric ozone column satisfies climatological minimums
+ IF( N_TROPO_O3_TOGGLE .GT. 0 )THEN
+ O3_TOGGLE_AVE = O3_TOGGLE_AVE / REAL( N_TROPO_O3_TOGGLE )
+ WRITE( LOGDEV, 9500 )'PHOT: Exceedance of tropospheric ozone column ',
+ & 'or below top of model domains based on stratospheric column minimum ',
+ & 'at date and time; ', JDATE, JTIME, N_TROPO_O3_TOGGLE, (1.0/O3_TOGGLE_AVE - 1.0),
+ & (1.0/O3_TOGGLE_MIN - 1.0)
+ END IF
+
+!...write diagnostic data to output file at the end of every output tstep
+
+ IF ( JTIME_CHK ) THEN
+
+ VARNM = 'COSZENS'
+ IF ( .NOT. WRITE3( CTM_RJ_1, VARNM, ODATE, OTIME,
+ & COSINE_ZENITH ) ) THEN
+ XMSG = 'Error writing variable ' // VARNM
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ END IF
+
+ VARNM = 'OZONE_COLUMN'
+ IF ( .NOT. WRITE3( CTM_RJ_1, VARNM, ODATE, OTIME, TOTAL_OC ) ) THEN
+ XMSG = 'Error writing variable ' // VARNM
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ END IF
+
+ VARNM = 'TROPO_O3_COLUMN'
+ IF ( .NOT. WRITE3( CTM_RJ_1, VARNM, ODATE, OTIME, TROPO_OC ) ) THEN
+ XMSG = 'Error writing variable ' // VARNM
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ END IF
+
+ IMONTH = IMONTH + 1
+ IF( IMONTH .GT. 12 )THEN
+ IMONTH = 1
+ TDATE = 2011001
+ END IF
+ TDATE = TDATE + DAYS( IMONTH )
+! CALL SEASONAL_STRAT_O3(TDATE, JTIME )
+
+
+! VARNM = 'MIN_FRAC_STRATO3'
+! IF ( .NOT. WRITE3( CTM_RJ_1, VARNM, ODATE, OTIME, MONTH_STRAT_03_FRAC ) ) THEN
+! XMSG = 'Error writing variable ' // VARNM
+! CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+! END IF
+
+
+ VARNM = 'TRANS_DIFFUSE'
+ IF ( .NOT. WRITE3( CTM_RJ_1, VARNM, ODATE, OTIME, TRANSMIS_DIFFUSE ) ) THEN
+ XMSG = 'Error writing variable ' // VARNM
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ END IF
+
+ VARNM = 'TRANS_DIRECT'
+ IF ( .NOT. WRITE3( CTM_RJ_1, VARNM, ODATE, OTIME, TRANSMIS_DIRECT ) ) THEN
+ XMSG = 'Error writing variable ' // VARNM
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ END IF
+
+ VARNM = 'REFLECTION'
+ IF ( .NOT. WRITE3( CTM_RJ_1, VARNM, ODATE, OTIME, REFLECT_COEFF ) ) THEN
+ XMSG = 'Error writing variable ' // VARNM
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ END IF
+
+ VARNM = 'CLR_TRANS_DIF'
+ IF ( .NOT. WRITE3( CTM_RJ_1, VARNM, ODATE, OTIME, CLR_TRANSMISSION ) ) THEN
+ XMSG = 'Error writing variable ' // VARNM
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ END IF
+
+ VARNM = 'CLR_TRANS_DIR'
+ IF ( .NOT. WRITE3( CTM_RJ_1, VARNM, ODATE, OTIME, CLR_TRANS_DIRECT ) ) THEN
+ XMSG = 'Error writing variable ' // VARNM
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ END IF
+
+ VARNM = 'CLR_REFLECTION'
+ IF ( .NOT. WRITE3( CTM_RJ_1, VARNM, ODATE, OTIME, CLR_REFLECTION ) ) THEN
+ XMSG = 'Error writing variable ' // VARNM
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ END IF
+
+ VARNM = 'TROPO_O3_EXCEED'
+ TROPO_O3_EXCEED = TROPO_O3_EXCEED / REAL( MAX(1, TSTEP_COUNT) )
+ IF ( .NOT. WRITE3( CTM_RJ_1, VARNM, ODATE, OTIME, TROPO_O3_EXCEED ) ) THEN
+ XMSG = 'Error writing variable ' // VARNM
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ END IF
+ TROPO_O3_EXCEED = 0.0 ! reset sum and counter
+ TSTEP_COUNT = 0
+
+ VARNM = 'N_EXCEED_TROPO3'
+ IF ( .NOT. WRITE3( CTM_RJ_1, VARNM, ODATE, OTIME, N_EXCEED_TROPO3 ) ) THEN
+ XMSG = 'Error writing variable ' // VARNM
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ END IF
+ N_EXCEED_TROPO3 = 0.0 ! reset counter
+
+ VARNM = 'JNO2'
+ IF ( .NOT. WRITE3( CTM_RJ_1, VARNM, ODATE, OTIME, RJ( :,:,1, LNO2 ) ) ) THEN
+ XMSG = 'Error writing variable ' // VARNM
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ END IF
+
+ VARNM = 'JO3O1D'
+ IF ( .NOT. WRITE3( CTM_RJ_1, VARNM, ODATE, OTIME, RJ( :,:,1,LO3O1D ) ) ) THEN
+ XMSG = 'Error writing variable ' // VARNM
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ END IF
+
+ VARNM = 'RESOLVED_CFRAC'
+ IF ( .NOT. WRITE3( CTM_RJ_1, VARNM, ODATE, OTIME, CFRAC_2D ) ) THEN
+ XMSG = 'Error writing variable ' // VARNM
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ END IF
+
+ VARNM = 'RESOLVED_WBAR'
+ IF ( .NOT. WRITE3( CTM_RJ_1, VARNM, ODATE, OTIME, AVE_HYDROMETEORS ) ) THEN
+ XMSG = 'Error writing variable ' // VARNM
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ END IF
+
+ IF ( USE_ACM_CLOUD ) THEN
+ VARNM = 'SUBGRID_CFRAC'
+ IF ( .NOT. WRITE3( CTM_RJ_1, VARNM, ODATE, OTIME, ACM_CLOUDS ) ) THEN
+ XMSG = 'Error writing variable ' // VARNM
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ END IF
+ VARNM = 'SUBGRID_WBAR'
+ IF ( .NOT. WRITE3( CTM_RJ_1, VARNM, ODATE, OTIME, ACM_AVE_H2O ) ) THEN
+ XMSG = 'Error writing variable ' // VARNM
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ END IF
+ END IF
+
+ DO IWL = 1, NWL
+
+ VARNM = 'ETOT_SFC_W' // WLTXT( IWL )
+ IF ( .NOT. WRITE3( CTM_RJ_1, VARNM, ODATE,
+ & OTIME, ETOT_SFC_WL( :,:,IWL ) ) ) THEN
+ XMSG = 'Error writing variable ' // VARNM
+ CALL M3EXIT ( PNAME, ODATE, OTIME, XMSG, XSTAT1 )
+ END IF
+
+ VARNM = 'TAU_AERO_W' // WLTXT( IWL )
+ IF ( .NOT. WRITE3( CTM_RJ_1, VARNM, ODATE,
+ & OTIME, TAU_AERO_WL( :,:,IWL ) ) ) THEN
+ XMSG = 'Error writing variable ' // VARNM
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ END IF
+
+ VARNM = 'TAU_CLOUD_W' // WLTXT( IWL )
+ IF ( .NOT. WRITE3( CTM_RJ_1, VARNM, ODATE,
+ & OTIME, TAU_CLOUD_WL( :,:,IWL ) ) ) THEN
+ XMSG = 'Error writing variable ' // VARNM
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ END IF
+
+#ifdef phot_debug
+ VARNM = 'SSA_CLOUD_W' // WLTXT( IWL )
+ IF ( .NOT. WRITE3( CTM_RJ_1, VARNM, ODATE,
+ & OTIME, SSA_CLOUD_WL( :,:,IWL ) ) ) THEN
+ XMSG = 'Error writing variable ' // VARNM
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ END IF
+
+ VARNM = 'ASY_CLOUD_W' // WLTXT( IWL )
+ IF ( .NOT. WRITE3( CTM_RJ_1, VARNM, ODATE,
+ & OTIME, ASY_CLOUD_WL( :,:,IWL ) ) ) THEN
+ XMSG = 'Error writing variable ' // VARNM
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ END IF
+#endif
+
+ VARNM = 'TAU_TOT_W' // WLTXT( IWL )
+ IF ( .NOT. WRITE3( CTM_RJ_1, VARNM, ODATE,
+ & OTIME, TAU_TOT_WL( :,:,IWL ) ) ) THEN
+ XMSG = 'Error writing variable ' // VARNM
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ END IF
+
+ VARNM = 'TAUO3_TOP_W' // WLTXT( IWL )
+ IF ( .NOT. WRITE3( CTM_RJ_1, VARNM, ODATE,
+ & OTIME, TAUO3_TOP_WL( :,:,IWL ) ) ) THEN
+ XMSG = 'Error writing variable ' // VARNM
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ END IF
+
+ VARNM = 'ALBEDO_W' // WLTXT( IWL )
+ IF ( .NOT. WRITE3( CTM_RJ_1, VARNM, ODATE, OTIME,
+ & SURFACE_ALBEDO( IWL,:,: ) ) ) THEN
+ XMSG = 'Error writing variable ' // VARNM
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ END IF
+
+ END DO ! iwl
+
+ WRITE( LOGDEV, '( /5X, 3( A, :, 1X ), I8, ":", I6.6 )' )
+ & 'RJ Values written to', CTM_RJ_1,
+ & 'for date and time', ODATE, OTIME
+
+ DO IPHOT = 1, NPHOTAB
+ IF ( .NOT. WRITE3( CTM_RJ_2, PHOTAB( IPHOT ), ODATE,
+ & OTIME, RJ( :,:,:,IPHOT ) ) ) THEN
+ XMSG = 'Could not write ' // CTM_RJ_2 // ' file'
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ END IF
+ END DO
+
+ VARNM = 'CFRAC_3D'
+ IF ( .NOT. WRITE3( CTM_RJ_2, VARNM, ODATE, OTIME, CFRAC_3D ) ) THEN
+ XMSG = 'Could not write ' // TRIM( VARNM ) // ' to ' // CTM_RJ_2 // ' file'
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ END IF
+
+ DO IWL = 1, NWL
+ VARNM = 'ACTINIC_FX_W' // WLTXT( IWL )
+ IF ( .NOT. WRITE3( CTM_RJ_2, VARNM, ODATE, OTIME, ACTINIC_FX( :,:,:,IWL ) ) ) THEN
+ XMSG = 'Error writing variable ' // VARNM
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ END IF
+ END DO
+
+ DO L = 1, N_DIAG_WVL
+ IWL = DIAG_WVL( L )
+
+ VARNM = 'AERO_SSA_W' // WLTXT( IWL )
+ IF ( .NOT. WRITE3( CTM_RJ_2, VARNM, ODATE, OTIME, AERO_SSA( :,:,:,L ) ) ) THEN
+ XMSG = 'Error writing variable ' // VARNM
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ END IF
+
+ VARNM = 'AERO_ASYM_W' // WLTXT( IWL )
+ IF ( .NOT. WRITE3( CTM_RJ_2, VARNM, ODATE, OTIME, AERO_ASYM( :,:,:,L ) ) ) THEN
+ XMSG = 'Error writing variable ' // VARNM
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ END IF
+
+ VARNM = 'TAU_AERO_W' // WLTXT( IWL )
+ IF ( .NOT. WRITE3( CTM_RJ_2, VARNM, ODATE, OTIME, TAU_AERO( :,:,:,L ) ) ) THEN
+ XMSG = 'Error writing variable ' // VARNM
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ END IF
+
+ VARNM = 'TAU_W' // WLTXT( IWL )
+ IF ( .NOT. WRITE3( CTM_RJ_2, VARNM, ODATE, OTIME, TAU( :,:,:,L ) ) ) THEN
+ XMSG = 'Error writing variable ' // VARNM
+ CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
+ END IF
+ END DO
+
+ WRITE( LOGDEV, '( /5X, 3( A, :, 1X ), I8, ":", I6.6 )' )
+ & 'RJ and Optical Data written to', CTM_RJ_2,
+ & 'for date and time', ODATE, OTIME
+
+ END IF ! if JTIME_CHK
+
+1003 FORMAT( 8X, 'Processor ',I4.4,' is in darkness at ', I8.7, ':', I6.6,
+ & 1X, 'GMT - no photolysis')
+9500 FORMAT(3(/ A), I7, 1X, I6.6, 1X, / "Total Number: ", I9, ";Mean Value: ", F9.6,
+ & "; Max Value: ",F9.6 /)
+
+ RETURN
+ END SUBROUTINE PHOT
diff --git a/src/shr/aqm_config_mod.F90 b/src/shr/aqm_config_mod.F90
index 01aaeae..b5c9d96 100644
--- a/src/shr/aqm_config_mod.F90
+++ b/src/shr/aqm_config_mod.F90
@@ -34,7 +34,7 @@ module aqm_config_mod
logical :: biosw_yn = .false.
logical :: ctm_aod = .false.
logical :: ctm_depvfile = .false.
- logical :: ctm_photodiag = .false.
+ logical :: ctm_photdiag = .false. !IVAI
logical :: ctm_pmdiag = .false.
logical :: ctm_wb_dust = .false.
logical :: mie_optics = .false.
@@ -42,6 +42,7 @@ module aqm_config_mod
logical :: run_aero = .false.
logical :: run_rescld = .false.
logical :: verbose = .false.
+ logical :: canopy_yn = .false.
type(aqm_species_type), pointer :: species => null()
end type aqm_config_type
@@ -182,6 +183,17 @@ subroutine aqm_config_read(model, config, rc)
rcToReturn=rc)) &
return ! bail out
+!IVAI
+ ! -- read diagnostic settings
+ call ESMF_ConfigGetAttribute(cf, config % ctm_photdiag, &
+ label="ctm_photdiag:", default=.false., rc=localrc)
+ if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, &
+ line=__LINE__, &
+ file=__FILE__, &
+ rcToReturn=rc)) &
+ return ! bail out
+!IVAI
+
call ESMF_ConfigGetAttribute(cf, config % ctm_pmdiag, &
label="ctm_pmdiag:", default=.false., rc=localrc)
if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, &
@@ -197,7 +209,16 @@ subroutine aqm_config_read(model, config, rc)
file=__FILE__, &
rcToReturn=rc)) &
return ! bail out
-
+
+ ! Canopy Options
+ call ESMF_ConfigGetAttribute(cf, config % canopy_yn, &
+ label="canopy_yn:", default=.false., rc=localrc)
+ if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, &
+ line=__LINE__, &
+ file=__FILE__, &
+ rcToReturn=rc)) &
+ return ! bail out
+
call ESMF_ConfigGetAttribute(cf, value, &
label="ctm_stdout:", default="all", rc=localrc)
if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, &
@@ -232,7 +253,6 @@ subroutine aqm_config_read(model, config, rc)
! -- set other default values
config % ctm_depvfile = .false.
- config % ctm_photodiag = .false.
end subroutine aqm_config_read
@@ -545,6 +565,25 @@ subroutine aqm_config_log(config, name, rc)
rcToReturn=rc)) &
return ! bail out
end if
+!IVAI
+ if (config % ctm_photdiag) then
+ call ESMF_LogWrite(trim(name) // ": config: read: ctm_photdiag: true", &
+ ESMF_LOGMSG_INFO, rc=localrc)
+ if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, &
+ line=__LINE__, &
+ file=__FILE__, &
+ rcToReturn=rc)) &
+ return ! bail out
+ else
+ call ESMF_LogWrite(trim(name) // ": config: read: ctm_photdiag: false", &
+ 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
+!IVAI
if (config % ctm_wb_dust) then
call ESMF_LogWrite(trim(name) // ": config: read: ctm_wb_dust: true", &
ESMF_LOGMSG_INFO, rc=localrc)
@@ -562,6 +601,23 @@ subroutine aqm_config_log(config, name, rc)
rcToReturn=rc)) &
return ! bail out
end if
+ if (config % canopy_yn) then
+ call ESMF_LogWrite(trim(name) // ": config: read: canopy_yn: true", &
+ ESMF_LOGMSG_INFO, rc=localrc)
+ if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, &
+ line=__LINE__, &
+ file=__FILE__, &
+ rcToReturn=rc)) &
+ return ! bail out
+ else
+ call ESMF_LogWrite(trim(name) // ": config: read: canopy_yn: false", &
+ 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_LogWrite(trim(name) // ": config: read: ctm_stdout: " &
// config % ctm_stdout, ESMF_LOGMSG_INFO, rc=localrc)
if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, &
diff --git a/src/shr/aqm_emis_mod.F90 b/src/shr/aqm_emis_mod.F90
index e432577..711264d 100644
--- a/src/shr/aqm_emis_mod.F90
+++ b/src/shr/aqm_emis_mod.F90
@@ -1600,6 +1600,12 @@ subroutine aqm_emis_grd_read(em, spcname, buffer, localDe, rc)
return ! bail out
end if
+ if (trim(em % type) == "canopy") then
+ ! -- ensure canopy 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
diff --git a/src/shr/aqm_methods.F90 b/src/shr/aqm_methods.F90
index 8501c2f..0733d16 100644
--- a/src/shr/aqm_methods.F90
+++ b/src/shr/aqm_methods.F90
@@ -88,7 +88,7 @@ LOGICAL FUNCTION DESC3( FNAME )
CHARACTER(LEN=*), INTENT(IN) :: FNAME
INCLUDE SUBST_FILES_ID
-
+
integer :: localrc
integer :: is, ie, js, je
integer :: EMLAYS
@@ -162,7 +162,7 @@ LOGICAL FUNCTION DESC3( FNAME )
ELSE IF ( TRIM( FNAME ) .EQ. TRIM( MET_CRO_2D ) ) THEN
- NVARS3D = 35
+ NVARS3D = 44
VNAME3D( 1:NVARS3D ) = &
(/ 'PRSFC ', 'USTAR ', &
'WSTAR ', 'PBL ', &
@@ -180,6 +180,11 @@ LOGICAL FUNCTION DESC3( FNAME )
'SEAICE ', 'SOIM1 ', &
'SOIM2 ', 'SOIT1 ', &
'SOIT2 ', 'LH ', &
+ 'FCH ', 'FRT ', &
+ 'CLU ', 'POPU ', &
+ 'LAIE ', 'C1R ', &
+ 'C2R ', 'C3R ', &
+ 'C4R ', &
'CLAYF ', 'SANDF ', &
'DRAG ', 'UTHR ' /)
UNITS3D( 1:NVARS3D ) = &
@@ -199,6 +204,11 @@ LOGICAL FUNCTION DESC3( FNAME )
'FRACTION ', 'M**3/M**3 ', &
'M**3/M**3 ', 'K ', &
'K ', 'WATTS/M**2 ', &
+ 'M ', 'NO UNIT ', &
+ 'NO UNIT ', 'PEOPLE/KM**2 ', &
+ 'NO UNIT ', 'NO UNIT ', &
+ 'NO UNIT ', 'NO UNIT ', &
+ 'NO UNIT ', &
'1 ', '1 ', &
'1 ', 'M/S ' /)
@@ -362,13 +372,15 @@ logical function envyn(name, description, defaultval, status)
envyn = config % ctm_depvfile
case ('CTM_PMDIAG')
envyn = config % ctm_pmdiag
- case ('CTM_PHOTODIAG')
- envyn = config % ctm_photodiag
+ case ('CTM_PHOTDIAG')
+ envyn = config % ctm_photdiag
case ('CTM_PT3DEMIS')
envyn = aqm_emis_ispresent("gbbepx") .or. &
aqm_emis_ispresent("point-source")
case ('CTM_GRAV_SETL')
envyn = .false.
+ case ('CTM_CANOPY_SHADE')
+ envyn = config % canopy_yn
case ('CTM_WBDUST_FENGSHA')
envyn = aqm_emis_ispresent("fengsha")
case ('CTM_WB_DUST')
@@ -539,7 +551,6 @@ INTEGER FUNCTION PROMPTFFILE( PROMPT, RDONLY, FMTTED, DEFAULT, CALLER )
END FUNCTION PROMPTFFILE
-
subroutine nameval(name, eqname)
use aqm_emis_mod, only : aqm_internal_emis_type, aqm_emis_get
@@ -588,7 +599,7 @@ subroutine nameval(name, eqname)
case default
! -- nothing to do
end select
-
+
end subroutine nameval
@@ -623,6 +634,7 @@ logical function interpx( fname, vname, pname, &
real(AQM_KIND_R8), dimension(:,:,:), pointer :: p3d
type(aqm_config_type), pointer :: config
type(aqm_state_type), pointer :: stateIn
+ type(aqm_state_type), pointer :: stateOut !IVAI
! -- constants
include SUBST_FILES_ID
@@ -639,6 +651,7 @@ logical function interpx( fname, vname, pname, &
nullify(p3d)
nullify(config)
nullify(stateIn)
+ nullify(stateOut) !IVAI
set_non_neg = .false.
if (trim(fname) == trim(GRID_CRO_2D)) then
@@ -703,7 +716,7 @@ logical function interpx( fname, vname, pname, &
call aqm_model_get(stateIn=stateIn, rc=localrc)
if (aqm_rc_check(localrc, msg="Failure to retrieve model input state", &
file=__FILE__, line=__LINE__)) return
-
+
call aqm_model_get(config=config, stateIn=stateIn, rc=localrc)
if (aqm_rc_check(localrc, msg="Failure to retrieve model input state", &
file=__FILE__, line=__LINE__)) return
@@ -804,15 +817,132 @@ logical function interpx( fname, vname, pname, &
end do
end do
case ("CLAYF","DRAG","SANDF","UTHR")
- ! -- read in fengsha variables
+ ! -- fengsha variables
call aqm_emis_read("fengsha", vname, buffer, rc=localrc)
if (aqm_rc_test((localrc /= 0), &
msg="Failure to read fengsha input for " // vname, &
file=__FILE__, line=__LINE__)) return
+ case ("FCH","FRT","CLU","POPU","LAIE","C1R","C2R","C3R","C4R")
+ ! -- canopy variables
+ if (config % canopy_yn) then
+ call aqm_emis_read("canopy", vname, buffer, rc=localrc)
+ if (aqm_rc_test((localrc /= 0), &
+ msg="Failure to read canopy input for " // vname, &
+ file=__FILE__, line=__LINE__)) return
+ else
+ buffer(1:lbuf) = 0.
+ end if
case default
! return
end select
+!IVAI
+ print*, 'AQM_METHODS: FNAME= ', FNAME, VNAME !IVAI : MET_CRO_2D
+
+ IF ( TRIM( VNAME ) .EQ. TRIM('LAIE') ) THEN
+
+ print*, 'AQM_METHODS: VNAME= ', VNAME !IVAI: LAIE
+! print*, 'AQM_METHODS: LAIE = ', buffer(1:lbuf)
+
+ nullify(stateOut)
+ call aqm_model_get(stateOut=stateOut, rc=localrc)
+ if (aqm_rc_check(localrc, msg="Failure to retrieve model output state", &
+ file=__FILE__, line=__LINE__)) return
+
+ k = 0
+ do r = row0, row1
+ do c = col0, col1
+ k = k + 1
+
+ stateOut % CLAIE (c,r) = buffer(k)
+ end do
+ end do
+
+
+ END IF
+ IF ( TRIM( VNAME ) .EQ. TRIM('FCH') ) THEN
+
+ print*, 'AQM_METHODS: VNAME= ', VNAME !IVAI: FCH
+! print*, 'AQM_METHODS: FCH = ', buffer(1:lbuf)
+
+ nullify(stateOut)
+ call aqm_model_get(stateOut=stateOut, rc=localrc)
+ if (aqm_rc_check(localrc, msg="Failure to retrieve model output state", &
+ file=__FILE__, line=__LINE__)) return
+
+ k = 0
+ do r = row0, row1
+ do c = col0, col1
+ k = k + 1
+
+ stateOut % CFCH (c,r) = buffer(k)
+ end do
+ end do
+
+ END IF
+ IF ( TRIM( VNAME ) .EQ. TRIM('FRT') ) THEN
+
+ print*, 'AQM_METHODS: VNAME= ', VNAME !IVAI: FRT
+! print*, 'AQM_METHODS: FRT = ', buffer(1:lbuf)
+
+ nullify(stateOut)
+ call aqm_model_get(stateOut=stateOut, rc=localrc)
+ if (aqm_rc_check(localrc, msg="Failure to retrieve model output state", &
+ file=__FILE__, line=__LINE__)) return
+
+ k = 0
+ do r = row0, row1
+ do c = col0, col1
+ k = k + 1
+
+ stateOut % CFRT(c,r) = buffer(k)
+ end do
+ end do
+
+ END IF
+ IF ( TRIM( VNAME ) .EQ. TRIM('CLU') ) THEN
+
+ print*, 'AQM_METHODS: VNAME= ', VNAME !IVAI: CLU
+! print*, 'AQM_METHODS: CLU = ', buffer(1:lbuf)
+
+ nullify(stateOut)
+ call aqm_model_get(stateOut=stateOut, rc=localrc)
+ if (aqm_rc_check(localrc, msg="Failure to retrieve model output state", &
+ file=__FILE__, line=__LINE__)) return
+
+ k = 0
+ do r = row0, row1
+ do c = col0, col1
+ k = k + 1
+
+ stateOut % CCLU(c,r) = buffer(k)
+ end do
+ end do
+
+ END IF
+ IF ( TRIM( VNAME ) .EQ. TRIM('POPU') ) THEN
+
+ print*, 'AQM_METHODS: VNAME= ', VNAME !IVAI: POPU
+! print*, 'AQM_METHODS: POPU= ', buffer(1:lbuf)
+
+ nullify(stateOut)
+ call aqm_model_get(stateOut=stateOut, rc=localrc)
+ if (aqm_rc_check(localrc, msg="Failure to retrieve model output state", &
+ file=__FILE__, line=__LINE__)) return
+
+ k = 0
+ do r = row0, row1
+ do c = col0, col1
+ k = k + 1
+
+ stateOut % CPOPU(c,r) = buffer(k)
+ end do
+ end do
+
+ END IF
+
+!IVAI
+
else if (trim(fname) == trim(OCEAN_1)) then
select case (trim(vname))
@@ -1210,8 +1340,12 @@ LOGICAL FUNCTION WRITE3_REAL2D( FNAME, VNAME, JDATE, JTIME, BUFFER )
WRITE3_REAL2D = .FALSE.
+! print*, 'AQM_METHODS: FNAME, VNAME= ', FNAME, VNAME !IVAI: CTM_AOD_1 ALL
+
IF ( TRIM( VNAME ) .EQ. TRIM( ALLVAR3 ) ) THEN
+! print*, 'AQM_METHODS: VNAME= ', VNAME !IVAI: ADO
+
nullify(stateOut)
call aqm_model_get(stateOut=stateOut, rc=localrc)
if (aqm_rc_check(localrc, msg="Failure to retrieve model output state", &
@@ -1219,12 +1353,78 @@ LOGICAL FUNCTION WRITE3_REAL2D( FNAME, VNAME, JDATE, JTIME, BUFFER )
stateOut % aod = BUFFER
+! print*, 'AQM_METHODS: AOD pointer = ', aod !IVAI
+! print*, 'AQM_METHODS: AOD = ', stateOut % aod !IVAI
+
END IF
WRITE3_REAL2D = .TRUE.
END IF
+!IVAI
+ WRITE3_REAL2D = .TRUE.
+
+ IF ( TRIM( FNAME ) .EQ. TRIM( CTM_RJ_1 ) ) THEN
+
+ WRITE3_REAL2D = .FALSE.
+
+ print*, 'AQM_METHODS: FNAME= ', FNAME, VNAME !IVAI: JO3O1D JNO2 ... (list of 15 vars)
+
+ IF ( TRIM( VNAME ) .EQ. TRIM('COSZENS') ) THEN
+
+! print*, 'AQM_METHODS: VNAME= ', VNAME !IVAI: COSZENS
+
+ nullify(stateOut)
+ call aqm_model_get(stateOut=stateOut, rc=localrc)
+ if (aqm_rc_check(localrc, msg="Failure to retrieve model output state", &
+ file=__FILE__, line=__LINE__)) return
+
+ stateOut % coszens = BUFFER
+
+! print*, 'AQM_METHODS: COSZENS pointer = ', coszens
+! print*, 'AQM_METHODS: COSZENS = ', BUFFER
+
+ END IF
+
+ IF ( TRIM( VNAME ) .EQ. TRIM('JO3O1D') ) THEN
+
+! print*, 'AQM_METHODS: VNAME= ', VNAME !IVAI: JO3O1D
+
+ nullify(stateOut)
+ call aqm_model_get(stateOut=stateOut, rc=localrc)
+ if (aqm_rc_check(localrc, msg="Failure to retrieve model output state", &
+ file=__FILE__, line=__LINE__)) return
+
+ stateOut % JO3O1D = BUFFER
+
+! print*, 'AQM_METHODS: JO3O1D pointer = ', JO3O1D
+! print*, 'AQM_METHODS: JO3O1D = ', BUFFER
+
+ END IF
+
+ IF ( TRIM( VNAME ) .EQ. TRIM('JNO2') ) THEN
+
+! print*, 'AQM_METHODS: VNAME= ', VNAME !IVAI: JNO2
+
+ nullify(stateOut)
+ call aqm_model_get(stateOut=stateOut, rc=localrc)
+ if (aqm_rc_check(localrc, msg="Failure to retrieve model output state", &
+ file=__FILE__, line=__LINE__)) return
+
+ stateOut % JNO2 = BUFFER
+
+! print*, 'AQM_METHODS: JNO2 pointer = ', JNO2
+! print*, 'AQM_METHODS: JNO2 = ', BUFFER
+
+ END IF
+
+ WRITE3_REAL2D = .TRUE.
+
+ END IF ! CTM_RJ_1
+
+!IVAI
+
END FUNCTION WRITE3_REAL2D
LOGICAL FUNCTION WRITE3_REAL4D( FNAME, VNAME, JDATE, JTIME, BUFFER )
diff --git a/src/shr/aqm_state_mod.F90 b/src/shr/aqm_state_mod.F90
index 9dc762f..fe243f8 100644
--- a/src/shr/aqm_state_mod.F90
+++ b/src/shr/aqm_state_mod.F90
@@ -46,11 +46,25 @@ module aqm_state_mod
real(AQM_KIND_R8), dimension(:,:,:,:), pointer :: tr => null()
+!IVAI
+ ! -- canopy variables
+ real(AQM_KIND_R8), dimension(:,:), pointer :: cfch => null()
+ real(AQM_KIND_R8), dimension(:,:), pointer :: cfrt => null()
+ real(AQM_KIND_R8), dimension(:,:), pointer :: cclu => null()
+ real(AQM_KIND_R8), dimension(:,:), pointer :: cpopu => null()
+ real(AQM_KIND_R8), dimension(:,:), pointer :: claie => null()
+!IVAI
+
! -- diagnostics
real(AQM_KIND_R8), dimension(:,:), pointer :: aod => null()
+!IVAI:
+ ! -- photolysis diagnostics
+ real(AQM_KIND_R8), dimension(:,:), pointer :: coszens => null()
+ real(AQM_KIND_R8), dimension(:,:), pointer :: jo3o1d => null()
+ real(AQM_KIND_R8), dimension(:,:), pointer :: jno2 => null()
+!IVAI
end type aqm_state_type
public
-
end module aqm_state_mod