From 18916a26dc8cf766009f82f293f97a76ecb70495 Mon Sep 17 00:00:00 2001 From: Wei Li Date: Fri, 20 Oct 2023 19:06:15 -0500 Subject: [PATCH 1/5] add localized files --- src/model/src/PT3D_DEFN.F | 106 +++---- src/model/src/PT3D_FIRE_DEFN.F | 223 ++++--------- src/model/src/PT3D_STKS_DEFN.F | 114 ++----- src/model/src/RUNTIME_VARS.F | 29 +- src/model/src/centralized_io_module.F | 49 ++- src/model/src/get_env_mod.f90 | 438 ++++++++++++++++++++++++++ 6 files changed, 628 insertions(+), 331 deletions(-) create mode 100644 src/model/src/get_env_mod.f90 diff --git a/src/model/src/PT3D_DEFN.F b/src/model/src/PT3D_DEFN.F index 0a7163d..39ed225 100644 --- a/src/model/src/PT3D_DEFN.F +++ b/src/model/src/PT3D_DEFN.F @@ -1,6 +1,6 @@ MODULE PT3D_DEFN - USE PT3D_DATA_MOD + !USE PT3D_DATA_MOD USE PT3D_FIRE_DEFN USE PT3D_STKS_DEFN @@ -8,31 +8,34 @@ MODULE PT3D_DEFN PRIVATE - PUBLIC PT3DEMIS, NPTGRPS, VDEMIS_PT, VDEMIS_PT_FIRE, PMEMIS_PT, - & PT3D_INIT, GET_PT3D_EMIS + PUBLIC PT3D_INIT, GET_PT3D_EMIS !, VDEMIS_PT + + CHARACTER( 16 ), ALLOCATABLE, SAVE :: STKENAME( : ) CONTAINS - FUNCTION PT3D_INIT ( N_SPC_EMIS, EMLAYS, JDATE, JTIME, TSTEP ) + FUNCTION PT3D_INIT ( JDATE, JTIME, TSTEP ) & RESULT ( SUCCESS ) USE GRID_CONF ! horizontal & vertical domain specifications - USE PTMAP - USE STK_EMIS, ONLY : STKSPC ! hourly point source emissions + USE STK_EMIS + USE DESID_VARS USE UTILIO_DEFN + USE PTBILIN, ONLY: NPTGRPS + USE RUNTIME_VARS IMPLICIT NONE - INTEGER, INTENT( IN ) :: N_SPC_EMIS ! total no. of model emissions species - INTEGER, INTENT( IN ) :: EMLAYS ! number of emissions layers +C INTEGER, INTENT( IN ) :: N_SPC_EMIS ! total no. of model emissions species +C INTEGER, INTENT( IN ) :: EMLAYS ! number of emissions layers INTEGER, INTENT( IN ) :: JDATE ! Julian date (YYYYDDD) INTEGER, INTENT( IN ) :: JTIME ! time (HHMMSS) INTEGER, INTENT( IN ) :: TSTEP ! output time step ! -- local variables INTEGER :: IOS - INTEGER :: N, S, V - INTEGER, ALLOCATABLE :: MAP( : ) + INTEGER :: N, ISRM, V +C INTEGER, ALLOCATABLE :: MAP( : ) CHARACTER( 240 ) :: XMSG = ' ' CHARACTER( 16 ) :: PNAME = 'PT3D_INIT ' ! procedure name @@ -44,78 +47,59 @@ FUNCTION PT3D_INIT ( N_SPC_EMIS, EMLAYS, JDATE, JTIME, TSTEP ) ! -- In-line 3D point source emissions? - PT3DEMIS = ENVYN( CTM_PT3DEMIS, - & 'Flag for in-line 3d point source emissions', - & .FALSE., IOS ) - IF ( PT3DEMIS ) THEN + IF ( NPTGRPS .GT. 0 ) THEN XMSG = 'Using in-line 3d point source emissions option' CALL M3MSG2( XMSG ) ELSE RETURN END IF - ! -- merge PM maps from fire and point-source emissions - IF ( .NOT. PTMAP_INIT( ) ) THEN - XMSG = 'Could not merge point source mappings' - CALL M3WARN ( PNAME, JDATE, JTIME, XMSG ) - SUCCESS = .FALSE.; RETURN - END IF +C Initialize stack emissions data + ALLOCATE ( STKENAME( NPTGRPS ), STAT = IOS ) ! stk emis files array + CALL CHECKMEM( IOS, 'STKENAME', PNAME ) + STKENAME = ' ' ! array - ! -- initialize emission types + ! Retrive Point Source Filenames from Emissions Filename Vector + DO N = 1, NPTGRPS + ISRM = MAP_PTtoISRM( N ) + STKENAME( N ) = DESID_STREAM_NAME( ISRM ) + END DO - IF ( .NOT. PT3D_FIRE_INIT ( N_SPC_EMIS, EMLAYS, JDATE, JTIME ) ) THEN - XMSG = 'Could not initialize fire emissions' + ! Initialize Point Source Emissions + IF ( .NOT. STK_EMIS_INIT( STKENAME, JDATE, JTIME ) ) THEN + XMSG = 'Could not initialize stack parameters' CALL M3WARN ( PNAME, JDATE, JTIME, XMSG ) SUCCESS = .FALSE.; RETURN END IF - IF ( .NOT. PT3D_STKS_INIT ( ) ) THEN - XMSG = 'Could not initialize point-source emissions' - CALL M3WARN ( PNAME, JDATE, JTIME, XMSG ) - SUCCESS = .FALSE.; RETURN - END IF - - ! -- allocate emission arrays - - ALLOCATE ( VDEMIS_PT( NCOLS, NROWS, EMLAYS, N_SPC_PTEM ), STAT = IOS ) - CALL CHECKMEM( IOS, 'VDEMIS_PT', PNAME ) - VDEMIS_PT = 0.0 ! array assignment - - ALLOCATE ( VDEMIS_PT_FIRE( NCOLS, NROWS, EMLAYS, N_SPC_PTEM ), STAT = IOS ) - CALL CHECKMEM( IOS, 'VDEMIS_PT_FIRE', PNAME ) - VDEMIS_PT_FIRE = 0.0 ! array assignment - - ALLOCATE ( PMEMIS_PT( NCOLS, NROWS, EMLAYS, N_SPC_PTPM ), STAT = IOS ) - CALL CHECKMEM( IOS, 'PMEMIS_PT', PNAME ) - PMEMIS_PT = 0.0 ! array - - ! -- get number of different file groups (sectors) - - NPTGRPS = 0 - - ! -- create point source internal array - - ALLOCATE( STKSPC( NPTGRPS ), STAT = IOS ) - CALL CHECKMEM( IOS, 'STKSPC', PNAME ) END FUNCTION PT3D_INIT - SUBROUTINE GET_PT3D_EMIS ( JDATE, JTIME, TSTEP ) + SUBROUTINE GET_PT3D_EMIS ( JDATE, JTIME, TSTEP,EMVAR_PT, ISRM, + & VDEMIS_PT, PTLAYS,L_DESID_DIAG) INTEGER, INTENT( IN ) :: JDATE, JTIME INTEGER, INTENT( IN ) :: TSTEP( 3 ) - -C ... initialize emission arrays ... - - VDEMIS_PT = 0.0 ! array assignment - VDEMIS_PT_FIRE = 0.0 ! array assignment - PMEMIS_PT = 0.0 ! array assignment - - CALL GET_PT3D_FIRE_EMIS ( JDATE, JTIME ) - CALL GET_PT3D_STKS_EMIS ( JDATE, JTIME ) + CHARACTER(16), INTENT( IN ) :: EMVAR_PT( : ) + INTEGER, INTENT( IN ) :: ISRM + + INTEGER, INTENT( OUT) :: PTLAYS + REAL, INTENT(INOUT) :: VDEMIS_PT( :,:,:,: ) + +! Local variables: + LOGICAL :: L_DESID_DIAG ! not really used + + + IF (ISRM == 2) !For Fire + & CALL GET_PT3D_FIRE_EMIS ( JDATE, JTIME, EMVAR_PT, ISRM, + & VDEMIS_PT, PTLAYS ) + + IF (ISRM == 3) !For point-source + & CALL GET_PT3D_STKS_EMIS ( JDATE, JTIME, EMVAR_PT, ISRM, + & VDEMIS_PT, PTLAYS ) END SUBROUTINE GET_PT3D_EMIS diff --git a/src/model/src/PT3D_FIRE_DEFN.F b/src/model/src/PT3D_FIRE_DEFN.F index 65fb0e6..8306838 100644 --- a/src/model/src/PT3D_FIRE_DEFN.F +++ b/src/model/src/PT3D_FIRE_DEFN.F @@ -1,23 +1,13 @@ MODULE PT3D_FIRE_DEFN - USE PT3D_DATA_MOD + ! USE PT3D_DATA_MOD IMPLICIT NONE - REAL, ALLOCATABLE :: VFRAC( :,:,: ) ! vertical fraction - REAL, ALLOCATABLE :: BUFFER( : ) ! emission buffer - -C Species names from input file used for point source non-PM emissions mapping - INTEGER, ALLOCATABLE :: PTEM_FIRE_MAP( : ) - INTEGER, ALLOCATABLE :: PTPM_FIRE_MAP( : ) C Emission layers for sources within domain INTEGER :: EMLYRS ! no. of emis layers - INTEGER :: PM_EMLYRS ! no. of emis layers for PM -C Emission counters - INTEGER :: N_GSPC_FIRE_EMIS ! number of gas species in diagnostic file - INTEGER :: N_SPC_FIRE_PTPM CHARACTER( 240 ) :: XMSG = ' ' @@ -29,100 +19,18 @@ MODULE PT3D_FIRE_DEFN PRIVATE - PUBLIC PT3D_FIRE_INIT, GET_PT3D_FIRE_EMIS + PUBLIC GET_PT3D_FIRE_EMIS CONTAINS C======================================================================= - FUNCTION PT3D_FIRE_INIT ( N_SPC_EMIS, EMLAYS, JDATE, JTIME ) - & RESULT ( SUCCESS ) - - USE AQM_EMIS_MOD - USE GRID_CONF ! horizontal & vertical domain specifications - USE CGRID_SPCS ! CGRID mechanism species - USE PTMAP ! defines pt src species mapping to VDEMIS* arrays - USE UTILIO_DEFN - - IMPLICIT NONE - -C Includes: -C None - -C Arguments: - INTEGER, INTENT( IN ) :: N_SPC_EMIS ! total no. of model emissions species - INTEGER, INTENT( IN ) :: EMLAYS ! number of emissions layers - INTEGER, INTENT( IN ) :: JDATE ! Julian date (YYYYDDD) - INTEGER, INTENT( IN ) :: JTIME ! time (HHMMSS) - - LOGICAL SUCCESS - -C Parameters: - -C Local Variables: - CHARACTER( 16 ) :: PNAME = 'PT3D_INIT ' ! procedure name - CHARACTER( 16 ) :: VNAME ! variable name buffer - CHARACTER( 16 ), SAVE, ALLOCATABLE :: STKGNAME( : ) ! stack groups file name - - INTEGER IOS ! i/o and allocate memory status - - - INTEGER IDX - INTEGER N, NSPC, NSPC1, NSPC2, NSPC3 - INTEGER S, S_OFFSET, V - - INTEGER, ALLOCATABLE :: MAP( : ) - INTEGER, ALLOCATABLE :: SPC_PTEM_FIRE_FAC( : ) - INTEGER, ALLOCATABLE :: SPC_PTEM_FIRE_MAP( : ) - - TYPE( AQM_INTERNAL_EMIS_TYPE ), POINTER :: EM - -C----------------------------------------------------------------------- - - SUCCESS = .TRUE. - -C check if emissions are being provided - - EM => AQM_EMIS_GET( ETYPE ) - IF ( .NOT.ASSOCIATED( EM ) ) RETURN - -C set number of emissions layers depending on whether plumerise is on - - CALL AQM_EMIS_DESC( ETYPE, NLAYS=EMLYRS ) - PM_EMLYRS = EMLYRS - -C get point source emission mapping - - IF (.NOT.PTMAP_TYPE_INIT( EM, - & N_GSPC_FIRE_EMIS, N_SPC_FIRE_PTPM, - & PTEM_FIRE_MAP, PTPM_FIRE_MAP, - & SPC_PTEM_FIRE_FAC, SPC_PTEM_FIRE_MAP, - & "FIRE" ) ) THEN - XMSG = 'Could not get point source mappings' - CALL M3WARN ( PNAME, JDATE, JTIME, XMSG ) - SUCCESS = .FALSE.; RETURN - END IF - - DEALLOCATE( SPC_PTEM_FIRE_FAC, SPC_PTEM_FIRE_MAP ) - ALLOCATE ( VFRAC( NCOLS,NROWS,EMLYRS ), STAT = IOS ) - CALL CHECKMEM( IOS, 'VFRAC', PNAME ) - VFRAC = 1.0 ! array - - ALLOCATE ( BUFFER( NCOLS * NROWS ), STAT = IOS ) - CALL CHECKMEM( IOS, 'BUFFER', PNAME ) - BUFFER = 0.0 ! array - - SUCCESS = .TRUE.; RETURN - - END FUNCTION PT3D_FIRE_INIT - -C======================================================================= - - SUBROUTINE GET_PT3D_FIRE_EMIS ( JDATE, JTIME ) + SUBROUTINE GET_PT3D_FIRE_EMIS ( JDATE, JTIME, EMVAR_PT, ISRM, + & VDEMIS_PT, PTLAYS ) ! Revision History. -! Aug 12, 15 D. Wong: added code to handle parallel I/O implementation +! Oct 20, 23 Wei Li: change to use DESID mapping C----------------------------------------------------------------------- @@ -134,19 +42,22 @@ SUBROUTINE GET_PT3D_FIRE_EMIS ( JDATE, JTIME ) USE RXNS_DATA, ONLY : MECHNAME !Get Chemical Mechanism Name USE GRID_CONF ! horizontal & vertical domain specifications USE CGRID_SPCS ! CGRID mechanism species - USE AERO_DATA, ONLY : PMEM_MAP_NAME - USE PTMAP ! defines pt src species mapping to VDEMIS* arrays + USE DESID_VARS, ONLY : DESID_N_ISTR USE UTILIO_DEFN + USE RUNTIME_VARS, only: LOGDEV IMPLICIT NONE C Arguments: INTEGER, INTENT( IN ) :: JDATE, JTIME - + CHARACTER(16), INTENT( IN ) :: EMVAR_PT( : ) + INTEGER, INTENT( IN ) :: ISRM + INTEGER, INTENT( OUT) :: PTLAYS + REAL, INTENT(INOUT) :: VDEMIS_PT( :,:,:,: ) C Parameters: C External functions: - INTEGER, EXTERNAL :: SETUP_LOGDEV +C INTEGER, EXTERNAL :: SETUP_LOGDEV C Local variables: CHARACTER( 16 ) :: PNAME = 'GET_PT3D_EMIS ' ! procedure name @@ -154,33 +65,47 @@ SUBROUTINE GET_PT3D_FIRE_EMIS ( JDATE, JTIME ) INTEGER IOS ! i/o and allocate memory status INTEGER L, S, V ! counters - INTEGER C, R, K, M, N + INTEGER C, R, K, M, N, NSPC INTEGER LOCALRC + REAL, ALLOCATABLE :: VFRAC( :,:,: ) ! vertical fraction + REAL, ALLOCATABLE :: BUFFER( : ) ! emission buffe + LOGICAL :: IS_NOT_NVPOA, SAVE_POC LOGICAL, SAVE :: FIRSTIME = .TRUE. - INTEGER, SAVE :: LOGDEV +C INTEGER, SAVE :: LOGDEV TYPE( AQM_INTERNAL_EMIS_TYPE ), POINTER :: EM C----------------------------------------------------------------------- - IF ( FIRSTIME ) THEN - FIRSTIME = .FALSE. - LOGDEV = SETUP_LOGDEV() - END IF +C IF ( FIRSTIME ) THEN +C FIRSTIME = .FALSE. +C LOGDEV = SETUP_LOGDEV() +C END IF + NULLIFY(EM) EM => AQM_EMIS_GET( ETYPE ) IF ( .NOT.ASSOCIATED( EM ) ) RETURN - + !number of species in this emission file + NSPC = SIZE( EM % TABLE, DIM=1 ) + CALL AQM_EMIS_DESC( ETYPE, NLAYS=EMLYRS ) + PTLAYS = EMLYRS C For each time step, compute the layer fractions... WRITE( XMSG,'(A, I7.6)' ) - & 'Calculating emissions point source layer fractions for', JTIME + & 'Calculating emissions point source layer fractions for',JTIME WRITE( LOGDEV,* ) ' ' CALL M3MSG2( XMSG ) + ALLOCATE ( VFRAC( NCOLS,NROWS,EMLYRS ), STAT = IOS ) + CALL CHECKMEM( IOS, 'VFRAC', PNAME ) + + ALLOCATE ( BUFFER( NCOLS * NROWS ), STAT = IOS ) + CALL CHECKMEM( IOS, 'BUFFER', PNAME ) + BUFFER = 0.0 ! array + C ... initialize vertical fraction arrays ... C ... fire emissions are added to surface only by default ... @@ -221,11 +146,12 @@ SUBROUTINE GET_PT3D_FIRE_EMIS ( JDATE, JTIME ) IS_NOT_NVPOA = ( INDEX( MECHNAME, 'NVPOA' ) .EQ. 0 ) - DO S = 1, N_GSPC_FIRE_EMIS - M = PTEM_FIRE_MAP( S ) + DO S = 1, DESID_N_ISTR + VNAME = EMVAR_PT( S ) + IF ( VNAME .EQ. '' ) CYCLE + M = INDEX1( VNAME, NSPC, EM % TABLE( :, 1 ) ) + IF ( M .GT. 0 ) THEN - N = PTEM_MAP( S ) -c IF ( N .GT. 0 ) THEN BUFFER = 0.0 CALL AQM_EMIS_READ( ETYPE, EM % TABLE( M, 1 ), & BUFFER, RC=LOCALRC ) @@ -234,58 +160,39 @@ SUBROUTINE GET_PT3D_FIRE_EMIS ( JDATE, JTIME ) & // TRIM( ETYPE ) // " emissions", & FILE=__FILE__, LINE=__LINE__)) RETURN C Read Non-Carbon Organic Matter too if POC is Requested - SAVE_POC = .FALSE. - IF ( IS_NOT_NVPOA .AND. EM % TABLE( M, 1 ) .EQ. 'POC' ) THEN - CALL AQM_EMIS_READ( ETYPE, 'PNCOM', BUFFER, RC=LOCALRC ) - IF ( AQM_RC_CHECK( LOCALRC, - & MSG="Failure while reading PNCOM" // - & " from " // TRIM( ETYPE ) // " emissions", - & FILE=__FILE__, LINE=__LINE__)) RETURN - SAVE_POC = IS_NOT_NVPOA - END IF - DO L = 1, EMLYRS +C SAVE_POC = .FALSE. +C IF ( IS_NOT_NVPOA .AND. EM % TABLE( M, 1 ) .EQ. 'POC' ) THEN +C CALL AQM_EMIS_READ( ETYPE, 'PNCOM', BUFFER, RC=LOCALRC ) +C IF ( AQM_RC_CHECK( LOCALRC, +C & MSG="Failure while reading PNCOM" // +C & " from " // TRIM( ETYPE ) // " emissions", +C & FILE=__FILE__, LINE=__LINE__)) RETURN +C SAVE_POC = IS_NOT_NVPOA +C END IF + DO L = 1, PTLAYS K = 0 - DO R = 1, MY_NROWS - DO C = 1, MY_NCOLS + DO R = 1, NROWS + DO C = 1, NCOLS K = K + 1 - VDEMIS_PT( C,R,L,N ) = VFRAC( C,R,L ) * BUFFER( K ) + VDEMIS_PT(S,L,C,R) = VFRAC( C,R,L ) * BUFFER( K ) END DO END DO END DO - IF ( SAVE_POC ) THEN - DO L = 1, EMLYRS - K = 0 - DO R = 1, MY_NROWS - DO C = 1, MY_NCOLS - K = K + 1 - VDEMIS_PT_FIRE( C,R,L,N ) = VDEMIS_PT( C,R,L,N ) - END DO - END DO - END DO - END IF + ! IF ( SAVE_POC ) THEN !comment out this for now + ! DO L = 1, PTLAYS + ! K = 0 + ! DO R = 1, NROWS + ! DO C = 1, NCOLS + ! K = K + 1 + ! VDEMIS_PT(S,L,C,R) = VDEMIS_PT(S,L,C,R) + VDEMIS_PT( S,L,C,R ) + ! END DO + ! END DO + ! END DO + ! END IF END IF END DO - -C ... aerosol species ... - - DO S = 1, N_SPC_FIRE_PTPM - N = PTPM_FIRE_MAP( S ) - V = PTPM_MAP( N ) - BUFFER = 0.0 - CALL AQM_EMIS_READ( ETYPE, PMEM_MAP_NAME( V ), BUFFER, RC=LOCALRC ) - IF ( AQM_RC_CHECK( LOCALRC, MSG="Failure while reading " // - & TRIM( PMEM_MAP_NAME( V ) ) // " from " // TRIM( ETYPE ) // " emissions", - & FILE=__FILE__, LINE=__LINE__)) RETURN - DO L = 1, PM_EMLYRS - K = 0 - DO R = 1, MY_NROWS - DO C = 1, MY_NCOLS - K = K + 1 - PMEMIS_PT( C,R,L,N ) = VFRAC( C,R,L ) * BUFFER( K ) - END DO - END DO - END DO - END DO + ! DEALLOCATE temporary variables + DEALLOCATE ( BUFFER, VFRAC ) RETURN diff --git a/src/model/src/PT3D_STKS_DEFN.F b/src/model/src/PT3D_STKS_DEFN.F index 3a0cf2a..a27ff2f 100644 --- a/src/model/src/PT3D_STKS_DEFN.F +++ b/src/model/src/PT3D_STKS_DEFN.F @@ -1,20 +1,12 @@ MODULE PT3D_STKS_DEFN USE GRID_CONF ! horizontal & vertical domain specifications - USE PTMAP + !USE PT3D_DATA_MOD USE AQM_EMIS_MOD, ONLY : AQM_EMIS_DESC, AQM_EMIS_GET, & AQM_EMIS_READ, AQM_INTERNAL_EMIS_TYPE IMPLICIT NONE - INTEGER :: N_GSPC_STKS_EMIS - INTEGER :: N_SPC_PT3DEM - INTEGER :: N_SPC_STKS_PTPM - - INTEGER, ALLOCATABLE :: SPC_PTEM_STKS_FAC( : ) - INTEGER, ALLOCATABLE :: SPC_PTEM_STKS_MAP( : ) - INTEGER, ALLOCATABLE :: PTEM_STKS_MAP( : ) - INTEGER, ALLOCATABLE :: PTPM_STKS_MAP( : ) REAL :: CNVTP ! intermediate combined conv. factor @@ -23,19 +15,19 @@ MODULE PT3D_STKS_DEFN PRIVATE - PUBLIC :: N_SPC_STKS_PTPM, PTPM_STKS_MAP - PUBLIC :: GET_PT3D_STKS_EMIS, PT3D_STKS_INIT + PUBLIC :: GET_PT3D_STKS_EMIS !, PT3D_STKS_INIT CONTAINS - SUBROUTINE GET_PT3D_STKS_EMIS( JDATE, JTIME ) + SUBROUTINE GET_PT3D_STKS_EMIS( JDATE, JTIME, EMVAR_PT, ISRM, + & VDEMIS_PT, PTLAYS ) USE CGRID_SPCS ! CGRID mechanism species USE UTILIO_DEFN USE AQM_RC_MOD - USE AERO_DATA, ONLY : PMEM_MAP_NAME USE ASX_DATA_MOD, ONLY: MET_DATA, GRID_DATA, CONVPA - USE PT3D_DATA_MOD + USE DESID_VARS, ONLY : DESID_N_ISTR + USE RUNTIME_VARS, only: LOGDEV IMPLICIT NONE @@ -46,19 +38,23 @@ SUBROUTINE GET_PT3D_STKS_EMIS( JDATE, JTIME ) C Arguments: INTEGER, INTENT( IN ) :: JDATE ! Julian date (YYYYDDD) INTEGER, INTENT( IN ) :: JTIME ! time (HHMMSS) + CHARACTER(16), INTENT( IN ) :: EMVAR_PT( : ) + INTEGER, INTENT( IN ) :: ISRM + INTEGER, INTENT( OUT) :: PTLAYS + REAL, INTENT(INOUT) :: VDEMIS_PT( :,:,:,: ) C External Functions: - INTEGER, EXTERNAL :: SETUP_LOGDEV +C INTEGER, EXTERNAL :: SETUP_LOGDEV C Local Variables: - INTEGER :: LOGDEV +C INTEGER :: LOGDEV INTEGER :: C, R, L, M, N, S, V, K, I ! loop induction variables INTEGER :: LBOT ! layer containing plume bottom INTEGER :: LTOP ! layer containing plume top INTEGER :: LPBL ! first L: ZF(L) above mixing layer - ONLY for REPORT INTEGER :: LSTK ! first L: ZF(L) > STKHT INTEGER :: IJ, IS - INTEGER :: NSRC, NP + INTEGER :: NSRC, NP, NSPC INTEGER :: IOS, LOCALRC INTEGER, POINTER :: IP( : ) INTEGER, POINTER :: JP( : ) @@ -92,6 +88,7 @@ SUBROUTINE GET_PT3D_STKS_EMIS( JDATE, JTIME ) REAL, PARAMETER :: USTARMIN = 0.1 ! Min valid value for USTAR CHARACTER( 8 ) :: CINT ! integer to character buffer for Cwarning messages CHARACTER( 16 ) :: PNAME = 'GET_PT3D_EMIS' + CHARACTER( 16 ) :: VNAME ! variable name buffer CHARACTER( 120 ) :: XMSG = ' ' CHARACTER( 10 ), PARAMETER :: BLANK10 = ' ' @@ -165,6 +162,7 @@ END SUBROUTINE PLSPRD C----------------------------------------------------------------------- + NULLIFY(EM) EM => AQM_EMIS_GET( ETYPE ) IF ( .NOT.ASSOCIATED( EM ) ) RETURN @@ -172,11 +170,12 @@ END SUBROUTINE PLSPRD IF ( FIRSTIME ) THEN - LOGDEV = SETUP_LOGDEV() +C LOGDEV = SETUP_LOGDEV() C set number of emissions layers depending on whether plumerise is on CALL AQM_EMIS_DESC( ETYPE, NLAYS=EMLYRS ) + PTLAYS = EMLYRS FIRSTIME = .FALSE. END IF @@ -400,13 +399,14 @@ END SUBROUTINE PLSPRD DEALLOCATE ( TFRAC ) - ! -- (1) non-PM emissions - - DO S = 1, N_GSPC_STKS_EMIS - M = PTEM_STKS_MAP( S ) + ! use DESID_N_ISTR instead of PTMAP + NSPC = SIZE( EM % TABLE, DIM=1 ) + DO S = 1, DESID_N_ISTR + VNAME = EMVAR_PT( S ) + IF ( VNAME .EQ. '' ) CYCLE + M = INDEX1( VNAME, NSPC, EM % TABLE( :, 1 ) ) IF ( M .LT. 1 ) CYCLE - N = PTEM_MAP( S ) - +C N = PTEM_MAP( S ) BUFFER = 0.0 CALL AQM_EMIS_READ( ETYPE, EM % TABLE( M, 1 ), BUFFER, RC=LOCALRC) IF ( AQM_RC_CHECK( LOCALRC, @@ -414,46 +414,26 @@ END SUBROUTINE PLSPRD & TRIM( EM % TABLE( M, 1 ) ) // " from " // & TRIM( ETYPE ) // " emissions", & FILE=__FILE__, LINE=__LINE__ ) ) RETURN - IF ( EM % TABLE( M, 1 ) .EQ. 'POC' ) THEN - CALL AQM_EMIS_READ( ETYPE, 'PNCOM', BUFFER, RC=LOCALRC) - IF ( AQM_RC_CHECK( LOCALRC, - & MSG="Failure while reading PNCOM emissions from" // - & TRIM( ETYPE ) // " emissions", - & FILE=__FILE__, LINE=__LINE__ ) ) RETURN - END IF +C IF ( EM % TABLE( M, 1 ) .EQ. 'POC' ) THEN +C CALL AQM_EMIS_READ( ETYPE, 'PNCOM', BUFFER, RC=LOCALRC) +C IF ( AQM_RC_CHECK( LOCALRC, +C & MSG="Failure while reading PNCOM emissions from" // +C & TRIM( ETYPE ) // " emissions", +C & FILE=__FILE__, LINE=__LINE__ ) ) RETURN +C END IF ! -- add emissions - DO L = 1, EMLYRS + DO L = 1, PTLAYS DO IS = 1, NSRC IJ = EM % IJMAP( IS ) C = EM % IP( IJ ) R = EM % JP( IJ ) - VDEMIS_PT( C,R,L,N ) = VDEMIS_PT( C,R,L,N ) + VDEMIS_PT( S,L,C,R ) = VDEMIS_PT( S,L,C,R ) & + VFRAC( IS, L ) * BUFFER( IJ ) END DO END DO END DO - ! -- (2) PM emissions - - DO S = 1, N_SPC_STKS_PTPM - N = PTPM_STKS_MAP( S ) - V = PTPM_MAP( N ) - BUFFER = 0.0 - CALL AQM_EMIS_READ( ETYPE, PMEM_MAP_NAME( V ), BUFFER, RC=LOCALRC ) - IF ( AQM_RC_CHECK( LOCALRC, MSG="Failure while reading " // - & TRIM( PMEM_MAP_NAME( V ) ) // " from " // TRIM( ETYPE ) // " emissions", - & FILE=__FILE__, LINE=__LINE__)) RETURN - DO L = 1, EMLYRS - DO IS = 1, NSRC - IJ = EM % IJMAP( IS ) - C = EM % IP( IJ ) - R = EM % JP( IJ ) - PMEMIS_PT( C,R,L,N ) = PMEMIS_PT( C,R,L,N ) - & + VFRAC( IS, L ) * BUFFER( IJ ) - END DO - END DO - END DO ! -- free up memory @@ -470,32 +450,4 @@ END SUBROUTINE PLSPRD END SUBROUTINE GET_PT3D_STKS_EMIS - FUNCTION PT3D_STKS_INIT( ) RESULT ( SUCCESS ) - - USE AQM_RC_MOD, ONLY: AQM_RC_TEST - - TYPE( AQM_INTERNAL_EMIS_TYPE ), POINTER :: EM - - LOGICAL :: SUCCESS - - SUCCESS = .TRUE. - - EM => AQM_EMIS_GET( ETYPE ) - IF ( .NOT.ASSOCIATED( EM ) ) RETURN - - IF (EM % COUNT == 0) RETURN - - SUCCESS = PTMAP_TYPE_INIT( EM, - & N_GSPC_STKS_EMIS, N_SPC_STKS_PTPM, - & PTEM_STKS_MAP, PTPM_STKS_MAP, - & SPC_PTEM_STKS_FAC, SPC_PTEM_STKS_MAP, - & "STKS" ) - - IF ( AQM_RC_TEST( .NOT. SUCCESS, - & MSG="Failure initializing mapping for" // - & TRIM( ETYPE ) // " emissions", - & FILE=__FILE__, LINE=__LINE__ ) ) RETURN - - END FUNCTION PT3D_STKS_INIT - END MODULE PT3D_STKS_DEFN diff --git a/src/model/src/RUNTIME_VARS.F b/src/model/src/RUNTIME_VARS.F index eee3d62..264b803 100644 --- a/src/model/src/RUNTIME_VARS.F +++ b/src/model/src/RUNTIME_VARS.F @@ -73,7 +73,7 @@ MODULE RUNTIME_VARS PUBLIC - INTEGER :: OUTDEV = 6 ! File Unit for Standard Output + INTEGER :: OUTDEV = -1 ! File Unit for Standard Output (AQM change from 6) INTEGER :: LOGDEV = -1 ! File Unit for Ascii Log File INTEGER :: TOTPE = 1 ! Number of Total Processors INTEGER :: NPROCS = 1 ! Number of Total Processors @@ -97,7 +97,7 @@ MODULE RUNTIME_VARS ! this is for MPAS LOGICAL :: ncd_64bit_offset = .FALSE. - INTEGER :: cell_num !beiming tang + !----------------------------------------------------------------------------------- !>> Define Environment Variables for Controlling CMAQ Processes !----------------------------------------------------------------------------------- @@ -341,7 +341,7 @@ SUBROUTINE INIT_ENV_VARS( JDATE, JTIME ) CALL MPI_COMM_SIZE ( MPI_COMM_WORLD, TOTPE, ERROR ) #else MYPE = 0 - TOTPE = 0 + TOTPE = 1 !(AQM change from 0) #endif IF ( MYPE .EQ. 0 ) VARDEV = OUTDEV @@ -349,8 +349,8 @@ SUBROUTINE INIT_ENV_VARS( JDATE, JTIME ) ! Get Simulation Scenario Name to Label Log Files, etc CALL GET_ENV( APPL_NAME, 'CTM_APPL', APPL_NAME, VARDEV ) - ! Start I/O-API and set up log file(s) - CALL SETUP_LOGDEV() + ! Start I/O-API and set up log file(s) + !CALL SETUP_LOGDEV() !comment out Otherwise will crash (AQM) IF ( MYPE .EQ. 0 ) THEN CALL LOG_HEADING( OUTDEV, "Environment Variable Report" ) CALL LOG_SUBHEADING( OUTDEV, "Grid and High-Level Model Parameters" ) @@ -609,7 +609,7 @@ SUBROUTINE INIT_ENV_VARS( JDATE, JTIME ) CALL GET_ENV ( NLAYS_DIAG, 'NLAYS_PHOTDIAG', NLAYS_DIAG, VARDEV ) ! Get Desired Wavelengths for Diagnostic Output - CALL GET_ENV( 'NWAVE_PHOTDIAG', NWAVE, WAVE_ENV, VARDEV ) !this could be wrong, why char before NWAVE, beiming??? + CALL GET_ENV( 'NWAVE_PHOTDIAG', NWAVE, WAVE_ENV, VARDEV ) END IF ! Get flag to use core-shell mixing model for aerosol optical properties @@ -925,17 +925,18 @@ SUBROUTINE INIT_ENV_VARS( JDATE, JTIME ) CALL M3EXIT( 'INIT_ENV_VARS', JDATE, JTIME, XMSG, EXIT_STATUS ) END IF - IF ( OCEAN_CHEM ) THEN - IF ( (INDEX( MECHNAME, 'CB6R5M_AE7_AQ') .GT. 0 ) .OR. - & (INDEX( MECHNAME, 'CB6R5_AE7_AQ' ) .GT. 0) ) then - USE_MARINE_GAS_EMISSION = .TRUE. - ENDIF - ENDIF + ! comment out since we do not have DMS emissions (AQM) + !IF ( OCEAN_CHEM ) THEN + ! IF ( (INDEX( MECHNAME, 'CB6R5M_AE7_AQ') .GT. 0 ) .OR. + !& (INDEX( MECHNAME, 'CB6R5_AE7_AQ' ) .GT. 0) ) then + ! USE_MARINE_GAS_EMISSION = .TRUE. + ! ENDIF + !ENDIF ! for MPAS #ifdef mpas - call get_env(ncd_64bit_offset, 'ncd_64bit_offset', .false., vardev) - call get_env(cell_num, 'cell_num', 1, vardev) + call get_env (ncd_64bit_offset, 'ncd_64bit_offset', .false., vardev) + call get_env( cell_num, 'cell_num', 1, vardev) #endif #ifdef twoway diff --git a/src/model/src/centralized_io_module.F b/src/model/src/centralized_io_module.F index 6121fbc..c0605e6 100644 --- a/src/model/src/centralized_io_module.F +++ b/src/model/src/centralized_io_module.F @@ -317,7 +317,7 @@ MODULE CENTRALIZED_IO_MODULE & ,retrieve_ocean_data_mpas #else & ,boundary_files_setup, - & retrieve_grid_cro_2d_data, + & retrieve_grid_cro_2d_data, !public this function (AQM) & retrieve_grid_dot_2d_data, & retrieve_ocean_data #endif @@ -560,6 +560,7 @@ subroutine gridded_files_setup file_xcell(f_met) = xcell3d file_ycell(f_met) = ycell3d + IF (INDEX1( 'TSEASFC', NVARS3D, VNAME3D ) .gt. 0) then TSEASFC_AVAIL = .true. adj = 0 @@ -775,8 +776,9 @@ subroutine gridded_files_setup ! emission file, could be one or multiple layer - call desid_read_namelist() - call desid_init_regions() + !called somewhere else (AQM) + !call desid_read_namelist() + !call desid_init_regions() allocate (cio_emis_file_name(N_FILE_GR), & cio_emis_file_loc(N_FILE_GR), @@ -1347,7 +1349,7 @@ subroutine retrieve_time_dep_gridded_data (jdate, jtime, vname) Call M3EXIT( PNAME, 0, 0, XMSG, XSTAT1 ) END IF #else - If ( .Not. XTRACT3( MET_CRO_2D, 'SLTYP', + If ( .Not. XTRACT3( MET_CRO_2D, 'SLTYP', & 1, 1, STRTROWMC2, ENDROWMC2, STRTCOLMC2, ENDCOLMC2, & jdate, jtime, SOILCAT_A ) ) THEN XMSG = ' Error interpolating variable SLTYP from ' // MET_CRO_2D @@ -1524,6 +1526,7 @@ subroutine retrieve_time_dep_gridded_data (jdate, jtime, vname) XMSG = 'Could not extract ' // cio_emis_file_name(fnum) // ' file' CALL M3EXIT ( PNAME, data_jdate, data_jtime, XMSG, XSTAT1 ) END IF + #endif else if (cio_grid_var_name(v,2) == 'e3d') then @@ -1586,7 +1589,7 @@ subroutine retrieve_time_dep_gridded_data (jdate, jtime, vname) end if data_jdate = loc_jdate(f_emis(fnum)) data_jtime = loc_jtime(f_emis(fnum)) - + IF ( .NOT. XTRACT3( cio_emis_file_name(fnum), loc_vname, & 1, cio_emis_file_layer(fnum), & cio_emis_file_startrow(fnum), cio_emis_file_endrow(fnum), @@ -1596,6 +1599,7 @@ subroutine retrieve_time_dep_gridded_data (jdate, jtime, vname) CALL M3EXIT ( PNAME, data_jdate, data_jtime, XMSG, XSTAT1 ) END IF + #endif else if (cio_grid_var_name(v,2) == 'ic') then @@ -1650,7 +1654,7 @@ subroutine retrieve_time_dep_gridded_data (jdate, jtime, vname) cio_grid_data_tstamp(1, buf_loc, v) = data_jdate cio_grid_data_tstamp(2, buf_loc, v) = data_jtime - end do + end do !end v #ifndef mpas ! assign TEMPG to TSEASFC when TSEASFC is not available in the input file @@ -4693,7 +4697,7 @@ end subroutine megan_setup subroutine centralized_io_init (in_ncols) use lsm_mod, only: n_lufrac, init_lsm - USE UTILIO_DEFN !, only : m3exit !beiming tang + USE UTILIO_DEFN, only : m3exit USE RUNTIME_VARS, only: log_heading, logdev #ifdef mpas @@ -4801,9 +4805,9 @@ subroutine centralized_io_init (in_ncols) call gridded_files_setup - call boundary_files_setup + !call boundary_files_setup !(AQM) - call stack_files_setup + !call stack_files_setup !(AQM) if (BIOGEMIS_BEIS) then call biogemis_setup @@ -4841,7 +4845,6 @@ subroutine centralized_io_init (in_ncols) call soilinp_setup end if - call retrieve_grid_cro_2d_data call retrieve_grid_dot_2d_data @@ -4860,9 +4863,9 @@ subroutine centralized_io_init (in_ncols) #ifdef mpas call retrieve_stack_data_mpas (cio_model_sdate, cio_model_stime) #else - call retrieve_boundary_data (cio_model_sdate, cio_model_stime) + !call retrieve_boundary_data (cio_model_sdate, cio_model_stime) !(AQM) - call retrieve_stack_data (cio_model_sdate, cio_model_stime) + !call retrieve_stack_data (cio_model_sdate, cio_model_stime) !(AQM) #endif end subroutine centralized_io_init @@ -6282,6 +6285,7 @@ subroutine r_interpolate_var_2d (vname, date, time, data, loc_size_spatial = size_c3d / nlays end if + if ((cio_grid_data_tstamp(1, loc_tail, var_loc) .lt. date) .or. & ((cio_grid_data_tstamp(2, loc_tail, var_loc) .lt. time) .and. & (cio_grid_data_tstamp(1, loc_tail, var_loc) .eq. date))) then @@ -6290,12 +6294,20 @@ subroutine r_interpolate_var_2d (vname, date, time, data, loc_jtime = cio_grid_data_tstamp(2, loc_tail, var_loc) CALL NEXTIME ( loc_jdate, loc_jtime, loc_tstep ) + call retrieve_time_dep_gridded_data (loc_jdate,loc_jtime, vname) - call retrieve_time_dep_gridded_data (loc_jdate, loc_jtime, vname) loc_head = head_grid(var_loc) loc_tail = tail_grid(var_loc) end if + !add to log + !write(-1,*)'loc_head and tail:',loc_head,loc_tail,var_loc,loc_tstep,time2sec(loc_tstep) + !write(-1,*)'cio_grid_data_tstamp(1, 2,var_loc):',cio_grid_data_tstamp(1, 2, var_loc) + !write(-1,*)'cio_grid_data_tstamp(2,2,var_loc):',cio_grid_data_tstamp(2, 2, var_loc) + !write(-1,*)'cio_grid_data_tstamp(1, loc_head, var_loc)',cio_grid_data_tstamp(1, loc_head, var_loc) + !write(-1,*)'cio_grid_data_tstamp(2, loc_head,var_loc)',cio_grid_data_tstamp(2, loc_head, var_loc) + !write(-1,*)'cio_grid_data_tstamp(2, loc_tail,var_loc)',cio_grid_data_tstamp(2, loc_tail, var_loc) + if ((cio_grid_data_tstamp(1, 2, var_loc) .eq. date) .and. & (cio_grid_data_tstamp(2, 2, var_loc) .eq. time)) then count = count + 1 @@ -6353,8 +6365,10 @@ subroutine r_interpolate_var_2d (vname, date, time, data, store_beg_ind = cio_grid_data_inx(1,2,var_loc) store_end_ind = cio_grid_data_inx(2,2,var_loc) - cio_grid_data(store_beg_ind:store_end_ind) = cio_grid_data(head_beg_ind:head_end_ind) * ratio1 - & + cio_grid_data(tail_beg_ind:tail_end_ind) * ratio2 + + cio_grid_data(store_beg_ind:store_end_ind) = ! cio_grid_data(head_beg_ind:head_end_ind) !* ratio1 + + & cio_grid_data(tail_beg_ind:tail_end_ind) !* ratio2 + end if adj_lvl = 0 @@ -6465,6 +6479,7 @@ subroutine i_interpolate_var_2d (vname, date, time, data) loc_jdate = cio_grid_data_tstamp(1, loc_tail, var_loc) loc_jtime = cio_grid_data_tstamp(2, loc_tail, var_loc) CALL NEXTIME ( loc_jdate, loc_jtime, loc_tstep ) + call retrieve_time_dep_gridded_data (loc_jdate, loc_jtime, vname) loc_head = head_grid(var_loc) loc_tail = tail_grid(var_loc) @@ -6887,8 +6902,8 @@ subroutine r_interpolate_var_3d (vname, date, time, data, fname) lcount = lcount + 1 end if - cio_grid_data(store_beg_ind:store_end_ind) = cio_grid_data(head_beg_ind:head_end_ind) * ratio1 - & + cio_grid_data(tail_beg_ind:tail_end_ind) * ratio2 + cio_grid_data(store_beg_ind:store_end_ind) = ! cio_grid_data(head_beg_ind:head_end_ind) !* ratio1 + + & cio_grid_data(tail_beg_ind:tail_end_ind) !* ratio2 end if end if diff --git a/src/model/src/get_env_mod.f90 b/src/model/src/get_env_mod.f90 new file mode 100644 index 0000000..6eca3a7 --- /dev/null +++ b/src/model/src/get_env_mod.f90 @@ -0,0 +1,438 @@ +!------------------------------------------------------------------------! +! 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. ! +!------------------------------------------------------------------------! + + module get_env_module + +! Function: get environment variables + +! Revision History: +! 2010 D.Wong: initial implementation +! 2 Feb 2010 D.Wong: provided an optional outputing device option, +! absorbed get_envlist function + + implicit none + + integer, parameter :: max_str_len = 10000 + + character (max_str_len) :: loc_str + + interface get_env + module procedure get_env_int, & + get_env_float, & + get_env_double, & + get_env_char, & + get_env_logical, & + get_envlist + end interface + + contains + +! -------------------------------------------------------------------------------- + subroutine get_env_int (env_value, env_var, default_env_value, logdev) + + integer, intent(out) :: env_value + character (*), intent(in) :: env_var + integer, intent(in) :: default_env_value + integer, intent(in), optional :: logdev + + integer :: loc_logdev + logical :: default, regular + + call getenv (env_var, loc_str) + + if (present(logdev)) then + loc_logdev = logdev + else + loc_logdev = 6 + end if + + regular = .false. + default = .false. + + !change the if.esle here + if ( trim(env_var) == 'N_EMIS_GR') then + env_value = 1 + regular = .true. + else if ( trim(env_var) == 'N_EMIS_PT') then + env_value = 2 + regular = .true. + else + env_value = default_env_value + default = .true. + end if + + if ( loc_logdev .gt. 0 ) then + if (default) then + write( loc_logdev, '(A16,2x,A,2x,i10, 1x, a9)' ), env_var,'|', env_value, '(default)' + else if (regular) then + write( loc_logdev, '(A16,2x,A,2x,i10)' ), env_var,'|', env_value + end if + end if + + end subroutine get_env_int + +! -------------------------------------------------------------------------------- + subroutine get_env_float (env_value, env_var, default_env_value, logdev) + + real, intent(out) :: env_value + character (*), intent(in) :: env_var + real, intent(in) :: default_env_value + integer, intent(in), optional :: logdev + + integer :: loc_logdev + logical :: default, regular + + call getenv (env_var, loc_str) + + if (present(logdev)) then + loc_logdev = logdev + else + loc_logdev = 6 + end if + + regular = .false. + default = .false. + + if (len(trim(loc_str)) == 0) then + env_value = default_env_value + default = .true. + else + read (loc_str, *) env_value + regular = .true. + end if + + if ( loc_logdev .gt. 0 ) then + if (default) then + write( loc_logdev, '(A16,2x,A,2x,e10.3, 1x, a9)' ), env_var,'|', env_value, '(default)' + else if (regular) then + write( loc_logdev, '(A16,2x,A,2x,e10.3)' ), env_var,'|', env_value + end if + end if + + end subroutine get_env_float + +! -------------------------------------------------------------------------------- + subroutine get_env_double (env_value, env_var, default_env_value, logdev) + + real (8), intent(out) :: env_value + character (*), intent(in) :: env_var + real, intent(in) :: default_env_value + integer, intent(in), optional :: logdev + + integer :: loc_logdev + logical :: default, regular + + call getenv (env_var, loc_str) + + if (present(logdev)) then + loc_logdev = logdev + else + loc_logdev = 6 + end if + + regular = .false. + default = .false. + + if (len(trim(loc_str)) == 0) then + env_value = default_env_value + default = .true. + else + read (loc_str, *) env_value + regular = .true. + end if + + if ( loc_logdev .gt. 0 ) then + if (default) then + write( loc_logdev, '(A16,2x,A,2x,e10.3, 1x, a9)' ), env_var,'|', env_value, '(default)' + else if (regular) then + write( loc_logdev, '(A16,2x,A,2x,e10.3)' ), env_var,'|', env_value + end if + end if + + end subroutine get_env_double + +! -------------------------------------------------------------------------------- + subroutine get_env_char (env_value, env_var, default_env_value, logdev) + + character (*), intent(out) :: env_value + character (*), intent(in) :: env_var + character (*), intent(in) :: default_env_value + integer, intent(in), optional :: logdev + + integer :: loc_logdev, length, STATUS !add STATUS for ENVSTR in aqm_methods + logical :: default, regular + character (50) :: myfmt + + call getenv (env_var, loc_str) + + if (present(logdev)) then + loc_logdev = logdev + else + loc_logdev = 6 + end if + + regular = .false. + default = .false. + + !change the if.esle here + if ( trim(env_var) == 'MISC_CTRL_NML') then + call nameval ('MISC_CTRL', env_value) + regular = .true. + else if ( trim(env_var) == 'DESID_CTRL_NML') then + call nameval ('DESID_CTRL', env_value) + regular = .true. + else if ( trim(env_var) == 'DESID_CHEM_CTRL_NML') then + call nameval ('DESID_CHEM_CTRL', env_value) + regular = .true. + else if ( trim(env_var) == 'gc_matrix_nml') then + call nameval ('gc_matrix_nml', env_value) + regular = .true. + else if ( trim(env_var) == 'ae_matrix_nml') then + call nameval ('ae_matrix_nml', env_value) + regular = .true. + else if ( trim(env_var) == 'nr_matrix_nml') then + call nameval ('nr_matrix_nml', env_value) + regular = .true. + else if ( trim(env_var) == 'tr_matrix_nml') then + call nameval ('tr_matrix_nml', env_value) + regular = .true. + else if ( trim(env_var) == 'GR_EMIS_LAB_001') then + env_value = 'GRIDDED_EMIS' + regular = .true. + else if ( trim(env_var) == 'STK_EMIS_LAB_001') then + env_value = 'PT_FIRES' + regular = .true. + else if ( trim(env_var) == 'STK_EMIS_LAB_002') then + env_value = 'PT_OTHER' + regular = .true. + else if ( trim(env_var) == 'GRID_NAME') then + !add spaces for 'GRID_NAME =STRTEMP(1:16)' in RUNTIME_VAR + env_value = 'Cubed-Sphere ' + regular = .true. + else if ( trim(env_var) == 'BIOG_SPRO') then + CALL ENVSTR('BIOG_SPRO', 'Speciation profile for biogenics', & + default_env_value, env_value, STATUS) + regular = .true. + else + env_value = default_env_value + default = .true. + end if + + if ( loc_logdev .gt. 0 ) then + length = len_trim(env_value) + if (default) then + if (length .eq. 0) then + write( loc_logdev, '(A16, 2x, A, 13x, a9)') env_var, '|', '(default)' + else + write (myfmt, '(a18, i3.3, a9)') '(A16, 2x, A, 2x, A', length, ', 1x, a9)' + write( loc_logdev, myfmt) env_var, '|', env_value, '(default)' + end if + else if (regular) then + write (myfmt, '(a18, i3.3, a1)') '(A16, 2x, A, 2x, A', length, ')' + write( loc_logdev, myfmt) env_var,'|', env_value + end if + end if + + end subroutine get_env_char + +! -------------------------------------------------------------------------------- + subroutine get_env_logical (env_value, env_var, default_env_value, logdev) + + logical, intent(out) :: env_value + character (*), intent(in) :: env_var + logical, intent(in) :: default_env_value + integer, intent(in), optional :: logdev + + integer :: length, status !add status for envyn in aqm_methods + integer :: loc_logdev + logical :: default, regular + LOGICAL, EXTERNAL :: ENVYN !add for envyn + + call getenv (env_var, loc_str) + + if (present(logdev)) then + loc_logdev = logdev + else + loc_logdev = 6 + end if + + length = len(trim(loc_str)) + regular = .false. + default = .false. + + !change the if.esle here + if (trim(env_var) == 'CTM_EMISCHK') then + env_value = .FALSE. + regular = .true. + else if (trim(env_var) == 'CTM_GRAV_SETL') then + env_value = ENVYN('CTM_GRAV_SETL','GRAV_SETL setting', & + default_env_value,STATUS) + regular = .true. + else if (trim(env_var) == 'CTM_WB_DUST') then + env_value = ENVYN('CTM_WB_DUST','In-line dust emissions', & + default_env_value,STATUS) + regular = .true. + else if (length <= 0) then + env_value = default_env_value + default = .true. + else if ((length == 1) .and. ((loc_str(1:1) .eq. 'Y') .or. & + (loc_str(1:1) .eq. 'y') .or. & + (loc_str(1:1) .eq. 'T') .or. & + (loc_str(1:1) .eq. 't'))) then + env_value = .true. + regular = .true. + else if ((length == 1) .and. ((loc_str(1:1) .eq. 'N') .or. & + (loc_str(1:1) .eq. 'n') .or. & + (loc_str(1:1) .eq. 'F') .or. & + (loc_str(1:1) .eq. 'f'))) then + env_value = .false. + regular = .true. + else if ((trim(loc_str) == '.TRUE.') .or. & + (trim(loc_str) == '.true.') .or. & + (trim(loc_str) == '.True.') .or. & + (trim(loc_str) == 'TRUE') .or. & + (trim(loc_str) == 'true') .or. & + (trim(loc_str) == 'True') .or. & + (trim(loc_str) == 'YES') .or. & + (trim(loc_str) == 'yes') .or. & + (trim(loc_str) == 'Yes')) then + env_value = .true. + regular = .true. + else if ((trim(loc_str) == '.FALSE.') .or. & + (trim(loc_str) == '.false.') .or. & + (trim(loc_str) == '.False.') .or. & + (trim(loc_str) == 'FALSE') .or. & + (trim(loc_str) == 'false') .or. & + (trim(loc_str) == 'False') .or. & + (trim(loc_str) == 'NO') .or. & + (trim(loc_str) == 'no') .or. & + (trim(loc_str) == 'No')) then + env_value = .false. + regular = .true. + else + write (loc_logdev, *) ' Note: Variable ', trim(env_var), ' improperly formatted' + env_value = default_env_value + default = .true. + end if + + if ( loc_logdev .gt. 0 ) then + if (default) then + write( loc_logdev, '(A16,2x,A,10x,L, 1x, a9)' ), env_var,'|', env_value, '(default)' + else if (regular) then + write( loc_logdev, '(A16,2x,A,10x,L)' ), env_var,'|', env_value + end if + end if + + end subroutine get_env_logical + +! -------------------------------------------------------------------------------- + subroutine get_envlist ( env_var, nvars, val_list, logdev ) + +! get a list env var (quoted string of items delimited by white space, +! commas or semi-colons) and parse out the items into variables. Two data +! types: character strings and integers (still represented as strings in +! the env var vaules). +! Examples: +! 1) setenv AVG_CONC_VARS "O3 NO NO2" +! 2) setenv AVG_CONC_LAYS "2 5" < start at two, end at 5 +! 3) setenv NPCOLSXNPROWS "4 3" +! 4) setenv BCOL_ECOL "3 8" +! 5) setenv BROW_EROW "2 10" +! 6) setenv BLAY_ELAY "1 5" + +! In example (1), not only parse out the named items "O3", "NO" and "NO2", +! but also obtain the count on the number of items (=3). + +! Revision: 2013/02/11 David Wong: increased the max env var length from 256 to 1000 +! 13 Dec 2013 J.Young: 1000 breaks BUFLEN in IOAPI's envgets.c. Change to 512. +! 17 Jun 2016 J.Young: IOAPI's envgets.c BUFLEN has been increased to 10000. +! 20 Jun 2016 J.Young: Forget IOAPI's envgets.c: use Fortran GETENV + + character( * ), intent ( in ) :: env_var + integer, intent ( out ) :: nvars + character( 16 ), intent ( out ) :: val_list( : ) + integer, intent(in), optional :: logdev + + integer :: max_len + character( 16 ) :: pname = 'GET_ENVLIST' + character( 16*size( val_list ) ) :: e_val + character( 1 ) :: chr + character( 96 ) :: xmsg + + integer :: jp( 16*size( val_list ) ), kp( 16*size( val_list ) ), status + integer ip, v + + integer :: loc_logdev + + if (present(logdev)) then + loc_logdev = logdev + else + loc_logdev = 6 + end if + + max_len = 16 * size( val_list ) + + call get_env( e_val, env_var, ' ', loc_logdev ) + + if ( env_var .eq. " " ) then + xmsg = 'Environment variable ' // env_var // ' not set' + call m3warn( pname, 0, 0, xmsg ) + nvars = 0 + return + end if + + nvars = 1 + + ip = 0 + +101 continue + ip = ip + 1 + if ( ip .gt. max_len ) go to 301 + chr = e_val( ip:ip ) + if ( chr .eq. ' ' .or. ichar ( chr ) .eq. 09 ) go to 101 + jp( nvars ) = ip ! 1st char + +201 continue + ip = ip + 1 + if ( ip .gt. max_len ) then + xmsg = 'Environment variable value too long' + call m3exit( pname, 0, 0, xmsg, 2 ) + end if + chr = e_val( ip:ip ) + if ( chr .ne. ' ' .and. & + chr .ne. ',' .and. & + chr .ne. ';' .or. & + ichar ( chr ) .eq. 09 ) then ! 09 = horizontal tab + go to 201 + else + kp( nvars ) = ip - 1 ! last char in this item + nvars = nvars + 1 + end if + + go to 101 + +301 continue + nvars = nvars - 1 + + do v = 1, nvars + val_list( v ) = e_val( jp( v ):kp( v ) ) + end do + + end subroutine get_envlist + + end module get_env_module From a419a75c184bd13f68184ea6f5709db631a81c8f Mon Sep 17 00:00:00 2001 From: Wei Li Date: Fri, 20 Oct 2023 19:18:07 -0500 Subject: [PATCH 2/5] Add other modified files --- CMakeLists.txt | 2 + aqm_files.cmake | 40 ++- src/aqm_comp_mod.F90 | 21 +- src/drv/cmaq_mod.F90 | 142 +++++----- src/drv/cmaq_model_mod.F90 | 19 +- src/io/ioapi/IODECL3.EXT | 319 ++++++++++++++++++++++ src/io/ioapi/daymon.F | 102 +++++++ src/io/ioapi/m3err.F | 59 ++++ src/io/ioapi/m3utilio.F90 | 53 +++- src/model/src/PTMAP.F | 341 ----------------------- src/shr/aqm_config_mod.F90 | 51 ++++ src/shr/aqm_emis_mod.F90 | 27 +- src/shr/aqm_methods.F90 | 538 ++++++++++++++++++++++++++++++++----- 13 files changed, 1217 insertions(+), 497 deletions(-) create mode 100644 src/io/ioapi/IODECL3.EXT create mode 100644 src/io/ioapi/daymon.F create mode 100644 src/io/ioapi/m3err.F delete mode 100644 src/model/src/PTMAP.F diff --git a/CMakeLists.txt b/CMakeLists.txt index 0682a02..7df5792 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -93,6 +93,7 @@ target_compile_definitions(CCTM PUBLIC SUBST_FILES_ID="FILES_CTM.EXT" SUBST_PE_COMM="PE_COMM.EXT" SUBST_COMM=NOOP_COMM SUBST_BARRIER=NOOP_BARRIER + SUBST_HI_LO_BND_PE=NOOP_HI_LO_BND_PE SUBST_SUBGRID_INDEX=NOOP_SUBGRID_INDEX AQCHEM=DUMMY_AQCHEM CONVCLD_ACM=DUMMY_CONVCLD_ACM @@ -105,6 +106,7 @@ target_compile_definitions(CCTM PUBLIC SUBST_FILES_ID="FILES_CTM.EXT" WR_INIT=DUMMY_WR_INIT verbose_aero verbose_gas + verbose_cio # mpas _AQM_) diff --git a/aqm_files.cmake b/aqm_files.cmake index 394512d..127f003 100644 --- a/aqm_files.cmake +++ b/aqm_files.cmake @@ -36,6 +36,7 @@ list(APPEND aqm_aqmio_files list(APPEND aqm_ioapi_files src/io/ioapi/FDESC3.EXT src/io/ioapi/PARMS3.EXT + src/io/ioapi/IODECL3.EXT src/io/ioapi/crlf.F src/io/ioapi/currec.f src/io/ioapi/currstep.f @@ -59,16 +60,19 @@ list(APPEND aqm_ioapi_files src/io/ioapi/upcase.f src/io/ioapi/wkday.F src/io/ioapi/yr2day.F + src/io/ioapi/daymon.F src/io/ioapi/m3exit.F90 src/io/ioapi/m3mesg.F90 src/io/ioapi/m3msg2.F90 src/io/ioapi/m3warn.F90 + src/io/ioapi/m3err.F src/io/ioapi/m3utilio.F90 ) set(CCTM_ROOT "src/model/CMAQ/CCTM/src") set(AERO "${CCTM_ROOT}/aero/aero6") set(BIOG "${CCTM_ROOT}/biog/beis4") +set(MEGAN "${CCTM_ROOT}/biog/megan3") set(CLOUD "${CCTM_ROOT}/cloud/acm_ae6") set(DEPV "${CCTM_ROOT}/depv/m3dry") set(EMIS "${CCTM_ROOT}/emis/emis") @@ -86,12 +90,14 @@ set(STENEX "${CCTM_ROOT}/STENEX/noop") set(UTIL "${CCTM_ROOT}/util/util") set(VDIFF "${CCTM_ROOT}/vdiff/acm2_m3dry") set(DRIV "${CCTM_ROOT}/driver") +set(CIO "${CCTM_ROOT}/cio") set(localCCTM "src/model/src") list(APPEND aqm_CCTM_files ${AERO}/AERO_DATA.F ${AERO}/aero_driver.F ${AERO}/aero_nml_modes.F ${AERO}/AEROMET_DATA.F + ${AERO}/AERO_EMIS.F ${AERO}/AEROSOL_CHEMISTRY.F ${AERO}/aero_subs.F ${AERO}/coags.f @@ -111,6 +117,13 @@ list(APPEND aqm_CCTM_files ${BIOG}/parsline.f ${BIOG}/tmpbeis.F ${BIOG}/wrdaymsg.f + ${MEGAN}/MEGAN_DEFN.F + ${MEGAN}/megan_gspro.F + ${MEGAN}/megan_hrno_mod.F + ${MEGAN}/megan_fx_mod.f90 + ${MEGAN}/BDSNP_MOD.F + ${MEGAN}/MAP_CV2CB6_AE7.EXT + ${MEGAN}/SPC_CB6_AE7.EXT ${CLOUD}/hlconst.F ${CLOUD}/cldproc_acm.F ${CLOUD}/getalpha.F @@ -144,12 +157,14 @@ list(APPEND aqm_CCTM_files ${EMIS}/UDTYPES.F ${EMIS}/biog_emis_param_module.F ${EMIS}/CMAQ_Control_DESID.nml - #${EMIS}/desid_module.F ${EMIS}/desid_param_module.F ${EMIS}/desid_util.F ${EMIS}/desid_vars.F + ${EMIS}/desid_module.F ${EMIS}/lus_data_module.F ${EMIS}/stack_group_data_module.F + # ${EMIS}/PT3D_DEFN.F + ${EMIS}/PTMET.F ${GAS}/../../reactive_tracers/DEGRADE_PARAMETERS.F ${GAS}/../../reactive_tracers/DEGRADE_ROUTINES.F ${GAS}/hrdata_mod.F @@ -219,10 +234,11 @@ list(APPEND aqm_CCTM_files ${STENEX}/noop_util_module.f ${UTIL}/findex.f ${UTIL}/log_header.F - ${UTIL}/get_env_mod.f90 + #${UTIL}/get_env_mod.f90 ${UTIL}/setup_logdev.F ${UTIL}/subhdomain.F ${UTIL}/UTILIO_DEFN.F + #${UTIL}/RUNTIME_VARS.F ${UTIL}/util_family_module.F ${UTIL}/CMAQ_Control_Misc.nml ${DRIV}/ELMO_PROC.F @@ -237,22 +253,22 @@ list(APPEND aqm_CCTM_files ${VDIFF}/VDIFF_DATA.F ${VDIFF}/VDIFF_DIAG.F ${VDIFF}/VDIFF_MAP.F - #${VDIFF}/vdiffacmx.F - #${VDIFF}/vdiffproc.F - ${VDIFF}/../../biog/megan3/BDSNP_MOD.F + ${VDIFF}/vdiffproc.F + #${CIO}/centralized_io_module.F ${localCCTM}/o3totcol.f - ${localCCTM}/AERO_EMIS.F - ${localCCTM}/RUNTIME_VARS.F + #${localCCTM}/AERO_EMIS.F #${localCCTM}/PTMAP.F #${localCCTM}/PT3D_DATA_MOD.F - #${localCCTM}/PT3D_DEFN.F - #${localCCTM}/PT3D_FIRE_DEFN.F - #${localCCTM}/PT3D_STKS_DEFN.F + ${localCCTM}/PT3D_DEFN.F + ${localCCTM}/PT3D_FIRE_DEFN.F + ${localCCTM}/PT3D_STKS_DEFN.F + ${localCCTM}/vdiffacmx.F ${localCCTM}/ASX_DATA_MOD.F ${localCCTM}/DUST_EMIS.F ${localCCTM}/AERO_PHOTDATA.F - ${localCCTM}/vdiffacmx.F - #${localCCTM}/phot.F + ${localCCTM}/phot.F + ${localCCTM}/RUNTIME_VARS.F + ${localCCTM}/get_env_mod.f90 ${localCCTM}/centralized_io_module.F ${localCCTM}/centralized_io_util_module.F ) diff --git a/src/aqm_comp_mod.F90 b/src/aqm_comp_mod.F90 index ecc81c5..62f4f8d 100644 --- a/src/aqm_comp_mod.F90 +++ b/src/aqm_comp_mod.F90 @@ -12,6 +12,10 @@ module aqm_comp_mod use aqm_prod_mod use aqm_internal_mod use cmaq_model_mod + USE CENTRALIZED_IO_MODULE !(CMAQ) + USE ELMO_PROC !(CMAQ) + USE RUNTIME_VARS, ONLY: STDATE, STTIME !(CMAQ) + implicit none @@ -147,6 +151,7 @@ subroutine aqm_comp_advance(model, rc) type(ESMF_Time) :: currTime type(ESMF_TimeInterval) :: timeStep type(aqm_config_type), pointer :: config => null() + logical, save :: first_step = .true. !(add new) ! -- begin rc = ESMF_SUCCESS @@ -176,8 +181,11 @@ subroutine aqm_comp_advance(model, rc) return ! bail out ! -- set model internal timestep vector (HHMMSS) - tstep( 1 ) = config % ctm_tstep ! TSTEP(1) = local output step - tstep( 2 ) = tstep( 1 ) ! TSTEP(2) = sciproc sync. step (chem) + !tstep( 1 ) = config % ctm_tstep ! TSTEP(1) = local output step + !tstep( 2 ) = tstep( 1 ) ! TSTEP(2) = sciproc sync. step (chem) + !test different tstep(1) and tstep(2) + tstep( 1 ) = 1 * 10000 ! TSTEP(1) = local output step + tstep( 2 ) = config % ctm_tstep ! TSTEP(2) = sciproc sync. step (chem) tstep( 3 ) = tstep( 2 ) ! TSTEP(3) = twoway model time step call ESMF_TimeGet(currTime, yy=yy, mm=mm, dd=dd, h=h, m=m, s=s, & @@ -198,6 +206,15 @@ subroutine aqm_comp_advance(model, rc) file=__FILE__)) & return ! bail out + !Initiliaze IO Arrays and Open Files (CMAQ) + if (first_step) then + first_step = .false. + !update STDATE(default is 1995192) and STTIME here for restart + STDATE = jdate + STTIME = jtime + call centralized_io_init + endif + ! -- advance CMAQ call cmaq_model_advance(jdate, jtime, tstep, rc=localrc) if (aqm_rc_check(localrc, file=__FILE__, line=__LINE__)) then diff --git a/src/drv/cmaq_mod.F90 b/src/drv/cmaq_mod.F90 index ac47610..ff3b705 100644 --- a/src/drv/cmaq_mod.F90 +++ b/src/drv/cmaq_mod.F90 @@ -14,16 +14,19 @@ module cmaq_mod use cgrid_spcs use AERO_DATA, only : aerolist, n_aerolist, & - aerospc_mw, n_emis_pm, & - map_pmemis, pmem_map, & - pmem_map_name, pmem_units + aerospc_mw use M3UTILIO, only : M3MESG use UTILIO_DEFN, only : INDEX1, INIT3, MXVARS3 + use DESID_VARS !(CMAQ) + USE CENTRALIZED_IO_MODULE !(CMAQ) + USE ELMO_PROC, ONLY : ELMO_DRIVER, WRITE_ELMO, MAP_ELMO !(CMAQ) implicit none integer :: cmaq_logdev + integer, save :: my_ncols + integer, save :: my_nrows ! -- pointer to CMAQ concentration array real, pointer :: CGRID(:,:,:,:) => null() @@ -65,6 +68,7 @@ subroutine cmaq_species_read(nspecies, rc) ! -- read from namelist CGRID gas chem, aerosol, non-reactive, ! -- and tracer species definitions ! -- This is done only on DE 0 and shared with other DEs on this PET + nspecies = 0 if (aqm_rc_test(.not.cgrid_spcs_init(), & msg="Error in CGRID_SPCS:CGRID_SPCS_INIT", & @@ -75,7 +79,10 @@ subroutine cmaq_species_read(nspecies, rc) file=__FILE__, line=__LINE__, rc=rc)) return nspecies = n_gc_trns + n_ae_trns + n_nr_trns - + + !add log + write(-1,* ) ' GC/AE/NR Species Number: ', n_gc_trns,'/',n_ae_trns,'/',n_nr_trns + end subroutine cmaq_species_read subroutine cmaq_init(rc) @@ -98,6 +105,10 @@ subroutine cmaq_init(rc) msg="Failure defining horizontal domain", & file=__FILE__, line=__LINE__, rc=rc)) return + !define 'my_ncols' here + my_ncols = NCOLS + my_nrows = NROWS + ! -- set I/O flag IO_PE_INCLUSIVE = ( MYPE .EQ. 0 ) @@ -107,7 +118,8 @@ subroutine cmaq_init(rc) ! -- Set up horizontal domain and calculate processor-to-subdomain maps for ! -- process analysis, if required IF ( LIPR .OR. LIRR ) THEN - IF (aqm_rc_test( .NOT. PAGRD_INIT( NPROCS, MYPE ), & +! IF (aqm_rc_test( .NOT. PAGRD_INIT( NPROCS, MYPE ), & + IF (aqm_rc_test( .NOT. PAGRD_INIT( MYPE ), & msg="Failure defining PA domain configuration", & FILE=__FILE__, LINE=__LINE__, rc=rc)) RETURN END IF @@ -119,6 +131,12 @@ subroutine cmaq_init(rc) CGRID => PCGRID( 1:MY_NCOLS,1:MY_NROWS,:,: ) ! required for PinG + !Initiliaze IO Arrays and Open Files (CMAQ) + call desid_read_namelist() + call desid_init_regions() + !Initialize ELMO Arrays and Maps + call map_elmo + end subroutine cmaq_init subroutine cmaq_advance(jdate, jtime, tstep, run_aero, run_rescld, rc) @@ -136,6 +154,12 @@ SUBROUTINE VDIFF ( CGRID, JDATE, JTIME, TSTEP ) INTEGER :: JDATE, JTIME INTEGER :: TSTEP( 3 ) END SUBROUTINE VDIFF + !add PHOT + SUBROUTINE PHOT ( CGRID, JDATE, JTIME, TSTEP ) + REAL, POINTER :: CGRID( :,:,:,: ) + INTEGER, INTENT( IN ) :: JDATE, JTIME + INTEGER, INTENT( IN ) :: TSTEP( : ) + END SUBROUTINE PHOT SUBROUTINE CLDPROC ( CGRID, JDATE, JTIME, TSTEP ) REAL, POINTER :: CGRID( :,:,:,: ) INTEGER :: JDATE, JTIME @@ -160,6 +184,8 @@ END SUBROUTINE AERO CALL VDIFF ( CGRID, JDATE, JTIME, TSTEP ) if (run_rescld) then + !add PHOT before CLDPROC + CALL PHOT ( CGRID, JDATE, JTIME, TSTEP ) CALL CLDPROC ( CGRID, JDATE, JTIME, TSTEP ) end if @@ -169,6 +195,11 @@ END SUBROUTINE AERO CALL AERO ( CGRID, JDATE, JTIME, TSTEP ) end if + !add ELMO here + !ALways let INIT_TIME=.TRUE. here to output each small time step + CALL ELMO_DRIVER( CGRID, JDATE, JTIME, TSTEP, INIT_TIME=.TRUE.) + CALL WRITE_ELMO( JDATE, JTIME, TSTEP, INIT_TIME=.TRUE.) + end subroutine cmaq_advance subroutine cmaq_import(tracers, prl, phii, temp, start_index, rc) @@ -438,6 +469,9 @@ subroutine cmaq_emis_init(rc) integer :: localrc, stat integer :: item integer :: ltable, n, spc + integer :: IRULE + integer, save :: N_RULE + Character( 16 ) :: pmem_units !units for PM emissions for all species integer, allocatable :: umap(:) real(AQM_KIND_R4) :: ucnv type(aqm_internal_emis_type), pointer :: em @@ -476,30 +510,38 @@ subroutine cmaq_emis_init(rc) ! -- add internal units to emissions reference table ltable = size(em % table, dim=1) - ! -- define mapping of CMAQ aerosol and related emission species - call map_pmemis() ! -- set destination units for PM emissions for all species pmem_units = "G/S" ! -- set internal units for all species - ! -- (a) gas species - do n = 1, n_gc_emis - spc = index1( gc_emis( n ), ltable, em % table(:,1) ) - if (spc > 0) em % table(spc,2) = "MOL/S" - end do - ! -- (b) non reactive - do n = 1, n_nr_emis - spc = index1( nr_emis( n ), ltable, em % table(:,1) ) - if (spc > 0) em % table(spc,2) = "MOL/S" + ! --- use DESID_EMVAR_TABLE from DESID_VARS.F + ! Find Total Number of Rules + N_RULE = 0 + DO IRULE = 1,SIZE( DESID_RULES_NML ) + IF( DESID_RULES_NML( IRULE )%SPEC .EQ. '' ) EXIT + N_RULE = IRULE + END DO + do n = 1, DESID_N_EMVAR_TABLE + spc = index1(DESID_EMVAR_TABLE( n )%NAME, ltable, em % table(:,1) ) + if (spc > 0) then + do IRULE = 1,N_RULE + CALL UPCASE( DESID_RULES_NML( IRULE )%EMVAR ) + if (DESID_EMVAR_TABLE( n )%NAME .EQ. DESID_RULES_NML( IRULE)%EMVAR ) then + CALL UPCASE( DESID_RULES_NML( IRULE )%PHASE ) + if (DESID_RULES_NML( IRULE )%PHASE .EQ. 'GAS') then + em % table(spc,2) = "MOL/S" + else + em % table(spc,2) = pmem_units + endif + endif + end do + + endif end do + spc = index1( "NH3_FERT", ltable, em % table(:,1) ) - if (spc > 0) em % table(spc,2) = "MOL/S" - ! -- (c) aerosols - do n = 1, n_emis_pm - spc = index1( pmem_map_name( n ), ltable, em % table(:,1) ) - if (spc > 0) em % table(spc,2) = pmem_units - end do + if (spc > 0) em % table(spc,2) = "MOL/S" ! -- perform unit conversion for input species, if needed ! -- (a) map input species to internal species @@ -520,44 +562,13 @@ subroutine cmaq_emis_init(rc) end do ! -- (b) perform unit conversion for input species - ! --- 1. gas species - do n = 1, size(em % species) - if (umap(n) > 0) then - spc = index1( em % species(n), n_gc_emis, gc_emis ) - if (spc > 0) then - ucnv = aqm_units_conv( em % units(n), em % table(umap(n),2), gc_molwt(gc_emis_map(spc)), em % dens_flag(n) ) - if (aqm_rc_test(ucnv == 0._AQM_KIND_R4, & - msg=trim(em % species(n))//": invalid input units ("//trim(em % units(n))//")", & - file=__FILE__, line=__LINE__, rc=rc)) return - em % factors(n) = ucnv * em % factors(n) - umap(n) = 0 - end if - end if - end do - ! --- 2. non reactive - do n = 1, size(em % species) - if (umap(n) > 0) then - if ( trim(em % species(n)) == "NH3_FERT" ) then - spc = index1( "NH3", n_nr_emis, nr_emis ) - else - spc = index1( em % species(n), n_nr_emis, nr_emis ) - end if - if (spc > 0) then - ucnv = aqm_units_conv( em % units(n), em % table(umap(n),2), nr_molwt(nr_emis_map(spc)), em % dens_flag(n) ) - if (aqm_rc_test(ucnv == 0._AQM_KIND_R4, & - msg=trim(em % species(n))//": invalid input units ("//trim(em % units(n))//")", & - file=__FILE__, line=__LINE__, rc=rc)) return - em % factors(n) = ucnv * em % factors(n) - umap(n) = 0 - end if - end if - end do - ! --- 3. aerosols - do n = 1, size(em % species) + ! --- use DESID_EMVAR_TABLE from DESID_VARS.F; + do n = 1, size(em % species) if (umap(n) > 0) then - spc = index1( em % species(n), n_emis_pm, pmem_map_name ) + spc = index1(em % species(n), DESID_N_EMVAR_TABLE, DESID_EMVAR_TABLE( : )%NAME ) if (spc > 0) then - ucnv = aqm_units_conv( em % units(n), em % table(umap(n),2), aerospc_mw(pmem_map(spc)), em % dens_flag(n) ) + ucnv = aqm_units_conv( em % units(n), em % table(umap(n),2), DESID_EMVAR_TABLE( spc )%MW, em % dens_flag(n) ) + if (aqm_rc_test(ucnv == 0._AQM_KIND_R4, & msg=trim(em % species(n))//": invalid input units ("//trim(em % units(n))//")", & file=__FILE__, line=__LINE__, rc=rc)) return @@ -567,6 +578,7 @@ subroutine cmaq_emis_init(rc) end if end do + ! -- (c) free up memory deallocate(umap, stat=stat) if (aqm_rc_test(stat /= 0, & @@ -824,20 +836,20 @@ subroutine cmaq_prod_pm25( pm25, cgrid, frac, idx, nlays_in) integer :: i, ibeg, iend, imod, mode, spc integer :: c, r, l - ! -- local parameters + ! -- local parameters (updated species for CB6r5_aero7) character(len=*), parameter :: pm25_species(*) = & - (/ "ASO4I ", "ANO3I ", "ANH4I ", "ANAI ", "ACLI ", "AECI ", "AOTHRI ", & ! I-mode (Atken) - "ALVPO1I", "ASVPO1I", "ASVPO2I", "ALVOO1I", "ALVOO2I", "ASVOO1I", "ASVOO2I", & + (/ "ASO4I ", "ANO3I ", "ANH4I ", "ANAI ", "ACLI ", "AECI ", "AOTHRI ", "APOCI ", & ! I-mode (Atken) + "ALVPO1I", "ASVPO1I", "ASVPO2I", "ALVOO1I", "ALVOO2I", "ASVOO1I", "ASVOO2I", "APNCOMI", & "ASO4J ", "ANO3J ", "ANH4J ", "ANAJ ", "ACLJ ", "AECJ ", "AOTHRJ ", & ! J-mode (accum) "ALVPO1J", "ASVPO1J", "ASVPO2J", "ASVPO3J", "AIVPO1J", & - "AXYL1J ", "AXYL2J ", "AXYL3J ", "ATOL1J ", "ATOL2J ", "ATOL3J ", "ABNZ1J ", "ABNZ2J ", & - "ABNZ3J ", "AISO1J ", "AISO2J ", "AISO3J ", "ATRP1J ", "ATRP2J ", "ASQTJ ", "AALK1J ", & - "AALK2J ", "APAH1J ", "APAH2J ", "APAH3J ", "AORGCJ ", "AOLGBJ ", "AOLGAJ ", & + "AMTNO3J", "AMTHYDJ", "APOCJ ", "APNCOMJ", "AAVB1J ", "AAVB2J ", "AAVB3J ", "AAVB4J ", & + "AMT1J ", "AISO1J ", "AISO2J ", "AISO3J ", "AMT2J ", "AMT3J ", "ASQTJ ", "AMT4J ", & + "AMT5J ", "AMT6J ", "AORGCJ ", "AOLGBJ ", "AOLGAJ ", & "ALVOO1J", "ALVOO2J", "ASVOO1J", "ASVOO2J", "ASVOO3J", "APCSOJ ", & "AFEJ ", "ASIJ ", "ATIJ ", "ACAJ ", "AMGJ ", "AMNJ ", "AALJ ", "AKJ ", & "ASOIL ", "ACORS ", "ASEACAT", "ACLK ", "ASO4K ", "ANO3K ", "ANH4K " /) ! K-mode (coarse) - integer, parameter :: nspc(3) = (/ 14, 49, 7 /) + integer, parameter :: nspc(3) = (/ 16, 47, 7 /) ! -- begin pm25 = 0. diff --git a/src/drv/cmaq_model_mod.F90 b/src/drv/cmaq_model_mod.F90 index 305e7de..24fb0c5 100644 --- a/src/drv/cmaq_model_mod.F90 +++ b/src/drv/cmaq_model_mod.F90 @@ -4,7 +4,8 @@ module cmaq_model_mod use aqm_types_mod use aqm_model_mod use cmaq_mod - + use RUNTIME_VARS !(AQM) + implicit none private @@ -40,6 +41,9 @@ subroutine cmaq_model_init(rc) ! -- initialize CMAQ ! -- NOTE: CMAQ can only run on 1DE/PET domain decomposition (DE 0) + ! Initialize all runscript environmental variables (AQM) + CALL INIT_ENV_VARS( 0, 0 ) + ! -- initialize species from namelists on DE 0 call cmaq_species_read(numSpecies, rc=localrc) if (aqm_rc_check(localrc, msg="Failed to initialize CMAQ species", & @@ -52,6 +56,11 @@ subroutine cmaq_model_init(rc) call aqm_model_domain_get(nt=nt, rc=localrc) if (aqm_rc_check(localrc, msg="Failed to retrieve model domain on local DE", & file=__FILE__, line=__LINE__, rc=rc)) return + + !add log AQM + write( logdev,* ) ' NT number of tracer is: ', nt + write( logdev,* ) ' p_aqm_beg/ndiag is:', config % species %p_aqm_beg,'/',config % species % ndiag + if (aqm_rc_test((config % species % p_aqm_beg + numSpecies - 1 > nt), & msg="Coupling tracer fields cannot hold all the required species", & file=__FILE__, line=__LINE__, rc=rc)) return @@ -117,10 +126,10 @@ subroutine cmaq_model_advance(jdate, jtime, tstep, rc) ! -- import advected species mixing ratios if (config % init_conc .and. first_step) then call cmaq_conc_init(jdate, jtime, tstep, rc=localrc) - if (aqm_rc_check(localrc, msg="Failed to initialize concentrations", & - file=__FILE__, line=__LINE__, rc=rc)) return - first_step = .false. - if (config % verbose) call cmaq_conc_log(trim(config % name) // ": init") + if (aqm_rc_check(localrc, msg="Failed to initialize concentrations", & + file=__FILE__, line=__LINE__, rc=rc)) return + first_step = .false. + if (config % verbose) call cmaq_conc_log(trim(config % name) // ": init") else call cmaq_import(stateIn % tr, stateIn % prl, stateIn % phii, stateIn % temp, config % species % p_aqm_beg) if (config % verbose) call cmaq_conc_log(trim(config % name) // ": import") diff --git a/src/io/ioapi/IODECL3.EXT b/src/io/ioapi/IODECL3.EXT new file mode 100644 index 0000000..377c767 --- /dev/null +++ b/src/io/ioapi/IODECL3.EXT @@ -0,0 +1,319 @@ + +!......................................................................... +! Version "@(#)$Header$" +! EDSS/Models-3 I/O API. Copyright (C) 1992-2002 MCNC +! Distributed under the GNU LESSER GENERAL PUBLIC LICENSE version 2.1 +! See file "LGPL.txt" for conditions of use. +!.................................................................... +! INCLUDE FILE IODECL3.EXT +! +! +! DO NOT EDIT !! +! +! The EDSS/Models-3 I/O API depends in an essential manner +! upon the contents of this INCLUDE file. ANY CHANGES are +! likely to result in very obscure, difficult-to-diagnose +! bugs caused by an inconsistency between standard "libioapi.a" +! object-libraries and whatever code is compiled with the +! resulting modified INCLUDE-file. +! +! By making any changes to this INCLUDE file, the user +! explicitly agrees that in the case any assistance is +! required of MCNC or of the I/O API author, Carlie J. Coats, Jr. +! as a result of such changes, THE USER AND/OR HIS PROJECT OR +! CONTRACT AGREES TO REIMBURSE MCNC AND/OR THE I/O API AUTHOR, +! CARLIE J. COATS, JR., AT A RATE TRIPLE THE NORMAL CONTRACT +! RATE FOR THE SERVICES REQUIRED. +! +! CONTAINS: declarations and usage comments for the Models-3 (M3) +! Interprocess Communication Applications Programming +! Interface (API) +! +! DEPENDENT UPON: consistency with the API itself. +! +! RELATED FILES: PARM3.EXT, FDESC3.EXT +! +! REVISION HISTORY: +! prototype 3/1992 by Carlie J. Coats, Jr., MCNC Environmental +! Programs +! +! Modified 2/2002 by CJC: updated dates, license, compatibility +! with both free and fixed Fortran 9x source forms +! +!.................................................................... + + LOGICAL CHECK3 ! is JDATE:JTIME available for FNAME? + LOGICAL CLOSE3 ! close FNAME + LOGICAL DESC3 ! Puts M3 file descriptions into FDESC3.EXT + LOGICAL FILCHK3 ! check file type and dimensions + INTEGER INIT3 ! Initializes M3 API and returns unit for log + LOGICAL SHUT3 ! Shuts down API + LOGICAL OPEN3 ! opens an M3 file + LOGICAL READ3 ! read M3 file for variable,layer,timestep + LOGICAL WRITE3 ! write timestep to M3 file + LOGICAL XTRACT3 ! extract window from timestep in a M3 file + LOGICAL INTERP3 ! do time interpolation from a M3 file + LOGICAL DDTVAR3 ! do time derivative from M3 file + + LOGICAL INTERPX ! time interpolation from a window + ! extraction from an M3 gridded file +!! LOGICAL PINTERPB ! parallel time interpolation from an + ! M3 boundary file + + LOGICAL INQATT3 ! inquire attributes in M3 file + LOGICAL RDATT3 ! read numeric attributes by name from M3 file + LOGICAL WRATT3 ! add new numeric attributes " + LOGICAL RDATTC ! read CHAR attributes " + LOGICAL WRATTC ! add new CHAR attributes " + + LOGICAL SYNC3 ! flushes file to disk, etc. + + EXTERNAL CHECK3 , CLOSE3, DESC3 , FILCHK3, INIT3 , + & SHUT3 , OPEN3 , READ3 , WRITE3 , XTRACT3, + & INTERP3, DDTVAR3, INQATT3, RDATT3 , WRATT3 , + & RDATTC , WRATTC, SYNC3, INTERPX ! , PINTERPB + +!....................................................................... +!.................. API FUNCTION USAGE AND EXAMPLES .................. +!....... +!....... In the examples below, names (FILENAME, PROGNAME, VARNAME) +!....... should be CHARACTER*16, STATUS and RDFLAG are LOGICAL, dates +!....... are INTEGER, coding the Julian date as YYYYDDD, times are +!....... INTEGER, coding the time as HHMMSS, and LOGDEV is the FORTRAN +!....... INTEGER unit number for the program's log file; and layer, +!....... row, and column specifications use INTEGER FORTRAN array +!....... index conventions (in particular, they are based at 1, not +!....... based at 0, as in C). +!....... Parameter values for "ALL...", for grid and file type IDs, +!....... and for API dimensioning values are given in PARMS3.EXT; +!....... file descriptions are passed via commons BDESC3 and CDESC3 +!....... in file FDESC3.EXT. +!....... +!....... CHECK3(): check whether timestep JDATE:JTIME is available +!....... for variable VNAME in file FILENAME. +!....... FORTRAN usage is: +!....... +!....... STATUS = CHECK3 ( FILENAME, VNAME, JDATE, JTIME ) +!....... IF ( .NOT. STATUS ) THEN +!....... ... (data-record not available in file FNAME) +!....... END IF +!....... +!....... CLOSE3(): check whether timestep JDATE:JTIME is available +!....... for variable VNAME in file FILENAME. +!....... FORTRAN usage is: +!....... +!....... STATUS = CLOSE3 ( FILENAME ) +!....... IF ( .NOT. STATUS ) THEN +!....... ... could not flush file to disk successfully, +!....... or else file not currently open. +!....... END IF +!....... +!....... DESC3(): return description of file FILENAME to the user +!....... in commons BDESC3 and CDESC3, file FDESC3.EXT. +!....... FORTRAN usage is: +!....... +!....... STATUS = DESC3 ( FILENAME ) +!....... IF ( .NOT. STATUS ) THEN +!....... ... (file not yet opened) +!....... END IF +!....... ... +!....... (Now common FDESC3 (file FDESC3.EXT) contains the descriptive +!....... information for this file.) +!....... +!....... FILCHK3(): check whether file type and dimensions for file +!....... FILENAME match the type and dimensions supplied by the user. +!....... FORTRAN usage is: +!....... +!....... STATUS = FILCHK3 ( FILENAME, FTYPE, NCOLS, NROWS, NLAYS, NTHIK ) +!....... IF ( .NOT. STATUS ) THEN +!....... ... (file type and dimensions do not match +!....... the supplied FTYPE, NCOLS, NROWS, NLAYS, NTHIK) +!....... END IF +!....... ... +!....... +!....... INIT3(): set up the M3 API, open the program's log file, and +!....... return the unit FORTRAN number for log file. May be called +!....... multiple times (in which case, it always returns the log-file's +!....... unit number). Note that block data INITBLK3.FOR must also be +!....... linked in. +!....... FORTRAN usage is: +!....... +!....... LOGDEV = INIT3 ( ) +!....... IF ( LOGDEV .LT. 0 ) THEN +!....... ... (can't proceed: probably can't open the log. +!....... Stop the program) +!....... END IF +!....... +!....... OPEN3(): open file FILENAME from program PROGNAME, with +!....... requested read-write/old-new status. For files opened for WRITE, +!....... record program-name and other history info in their headers. +!....... May be called multiple times for the same file (in which case, +!....... it returns true unless the request is for READ-WRITE status +!....... for a file already opened READ-ONLY). Legal statuses are: +!....... FSREAD3: "old read-only" +!....... FSRDWR3: "old read-write" +!....... FSNEW3: "new (read-write)" +!....... FSUNKN3: "unknown (read_write)" +!....... FORTRAN usage is: +!....... +!....... STATUS = OPEN3 ( FILENAME, FSTATUS, PROGNAME ) +!....... IF ( .NOT. STATUS ) THEN +!....... ... (process the error) +!....... END IF +!....... +!....... READ3(): read data from FILENAME for timestep JDATE:JTIME, +!....... variable VNAME, layer LAY, into location ARRAY. +!....... If VNAME==ALLVARS3=='ALL ', reads all variables; +!....... if LAY==ALLAYS3==-1, reads all layers. +!....... Offers random access to the data by filename, date&time, variable, +!....... and layer. For DICTIONARY files, logical name for file being +!....... requested maps into the VNAME argument. For time-independent +!....... files (including DICTIONARY files), JDATE and JTIME are ignored. +!....... FORTRAN usage is: +!....... +!....... STATUS = READ3 ( FILENAME, VNAME, LAY, JDATE, JTIME, ARRAY ) +!....... IF ( .NOT. STATUS ) THEN +!....... ... (read failed -- process this error.) +!....... END IF +!....... +!....... SHUT3(): Flushes and closes down all M3 files currently open. +!....... Must be called before program termination; if it returns FALSE +!....... the run must be considered suspect. +!....... FORTRAN usage is: +!....... +!....... STATUS = SHUT3 ( ) +!....... IF ( .NOT. STATUS ) THEN +!....... ... (Flush of files to disk probably didn't work; +!....... look at netCDF error messages) +!....... END IF +!....... +!....... WRITE3(): write data from ARRAY to file FILENAME for timestep +!....... JDATE:JTIME. For GRIDDED, BUONDARY, and CUSTOM files, VNAME +!....... must be a variable found in the file, or else ALLVARS3=='ALL' +!....... to write all variables from ARRAY. For other file types, +!....... VNAME _must_ be ALLVARS3. +!....... FORTRAN usage is: +!....... +!....... STATUS = WRITE3 ( FILENAME, VNAME, JDATE, JTIME, ARRAY ) +!....... IF ( .NOT. STATUS ) THEN +!....... ... (write failed -- process this error.) +!....... END IF +!....... +!....... XTRACT3(): read/extract gridded data into location ARRAY +!....... from FILENAME for time step JDATE:JTIME, variable VNAME +!....... and the data window defined by +!....... LOLAY <= layer <= HILAY, +!....... LOROW <= row <= HIROW, +!....... LOCOL <= column <= HICOL +!....... FORTRAN usage is: +!....... +!....... STATUS = XTRACT3 ( FILENAME, VNAME, +!....... & LOLAY, HILAY, +!....... & LOROW, HIROW, +!....... & LOCOL, HICOL, +!....... & JDATE, JTIME, ARRAY ) +!....... IF ( .NOT. STATUS ) THEN +!....... ... (extract failed -- process this error.) +!....... END IF +!....... +!....... INTERP3(): read/interpolate gridded, boundary, or custom data +!....... into location ARRAY from FILENAME for time JDATE:JTIME, variable +!....... VNAME, and all layers. Note use of ASIZE = transaction size = +!....... size of ARRAY, for error-checking. +!....... FORTRAN usage is: +!....... +!....... STATUS = INTERPX ( FILENAME, VNAME, CALLER, JDATE, JTIME, +!....... & ASIZE, ARRAY ) +!....... IF ( .NOT. STATUS ) THEN +!....... ... (interpolate failed -- process this error.) +!....... END IF +!....... +!....... INTERPX(): read/interpolate/window gridded, boundary, or custom +!....... data into location ARRAY from FILENAME for time JDATE:JTIME, +!....... variable VNAME, and all layers. +!....... FORTRAN usage is: +!....... +!....... STATUS = INTERPX ( FILENAME, VNAME, CALLER, +!....... & COL0, COL1, ROW0, ROW1, LAY0, LAY1, +!....... & JDATE, JTIME, ARRAY ) +!....... IF ( .NOT. STATUS ) THEN +!....... ... (windowed interpolate failed -- process this error.) +!....... END IF +!....... +!....... DDTVAR3(): read and calculate mean time derivative (per second) +!....... for gridded, boundary, or custom data. Put result into location +!....... ARRAY from FILENAME for time JDATE:JTIME, variable VNAME, and all +!....... layers. Note use of ASIZE = transaction size = size of ARRAY, +!....... for error-checking. Note d/dt( time-independent )==0.0 +!....... FORTRAN usage is: +!....... +!....... STATUS = DDTVAR3 ( FILENAME, VNAME, JDATE, JTIME, +!....... & ASIZE, ARRAY ) +!....... IF ( .NOT. STATUS ) THEN +!....... ... (operation failed -- process this error.) +!....... END IF +!....... +!....... INQATT(): inquire how many attributes there are for a +!....... particular file and variable (or for the file globally, +!....... if the variable-name ALLVAR3 is used)), and what the +!....... names, types, and array-dimensions of these attributes are. +!....... FORTRAN usage is: +!....... +!....... STATUS = INQATT3( FNAME, VNAME, MXATTS, +!....... & NATTS, ANAMES, ATYPES, ASIZES ) +!....... IF ( .NOT. STATUS ) THEN +!....... ... (operation failed -- process this error.) +!....... END IF +!....... +!....... RDATT3(): Reads an INTEGER, REAL, or DOUBLE attribute by name +!....... for a specified file and variable into a user-specified array. +!....... If variable name is ALLVAR3, reads the file-global attribute. +!....... FORTRAN usage is: +!....... +!....... STATUS = RDATT3( FNAME, VNAME, ANAME, ATYPE, AMAX, +!....... & ASIZE, AVAL ) +!....... IF ( .NOT. STATUS ) THEN +!....... ... (operation failed -- process this error.) +!....... END IF +!....... +!....... WRATT3(): Writes an INTEGER, REAL, or DOUBLE attribute by name +!....... for a specified file and variable. If variable name is ALLVAR3, +!....... reads the file-global attribute. +!....... +!....... STATUS = WRATT3( FNAME, VNAME, +!....... & ANAME, ATYPE, AMAX, AVAL ) +!....... IF ( .NOT. STATUS ) THEN +!....... ... (operation failed -- process this error.) +!....... END IF +!....... +!....... RDATTC(): Reads a CHARACTER string attribute by name +!....... for a specified file and variable into a user-specified array. +!....... If variable name is ALLVAR3, reads the file-global attribute. +!....... FORTRAN usage is: +!....... +!....... STATUS = RDATTC( FNAME, VNAME, ANAME, CVAL ) +!....... IF ( .NOT. STATUS ) THEN +!....... ... (operation failed -- process this error.) +!....... END IF +!....... +!....... WRATT3(): Writes a CHARACTER string attribute by name +!....... for a specified file and variable. If variable name is ALLVAR3, +!....... reads the file-global attribute. +!....... +!....... STATUS = WRATTC( FNAME, VNAME, ANAME, CVAL ) +!....... IF ( .NOT. STATUS ) THEN +!....... ... (operation failed -- process this error.) +!....... END IF +!....... +!....... SYNC3(): Synchronize FILENAME with disk (flush output; +!....... re-read header and invalidate data-buffers for input. +!....... FORTRAN usage is: +!....... +!....... STATUS = SYNC3 ( FILENAME ) +!....... IF ( .NOT. STATUS ) THEN +!....... ... (file not yet opened, or disk-synch failed) +!....... END IF +!....... ... +!....... +!................ end IODECL3.EXT .................................... + diff --git a/src/io/ioapi/daymon.F b/src/io/ioapi/daymon.F new file mode 100644 index 0000000..0cc905f --- /dev/null +++ b/src/io/ioapi/daymon.F @@ -0,0 +1,102 @@ + + SUBROUTINE DAYMON( JDATE, MNTH, MDAY ) + +C*********************************************************************** +C Version "$Id: daymon.F 219 2015-08-17 18:05:54Z coats $" +C EDSS/Models-3 I/O API. +C Copyright (C) 1992-2002 MCNC and Carlie J. Coats, Jr., and +C (C) 2003-2013 Baron Advanced Meteorological Systems +C Distributed under the GNU LESSER GENERAL PUBLIC LICENSE version 2.1 +C See file "LGPL.txt" for conditions of use. +C......................................................................... +C function body starts at line 49 +C +C FUNCTION: +C +C This routine determines the month and day of the month +C for the Julian date YYYYDDD that is input +C +C REVISION HISTORY: +C +C 3/1995 Adapted for Models-3/EDSS from ROM GREG.FOR by CJC +C +C 2/2002 Unification by CJC with global-climate DAYMON, which +C uses a 360-day "year" +C +C Version 1/2007 by CJC: handle negative JDATEs correctly +C +C Modified 03/2010 by CJC: F9x changes for I/O API v3.1 +C +C Version 5/2013 by CJC: handle standard-year cases +C*********************************************************************** + + + IMPLICIT NONE + +C........... ARGUMENTS and their descriptions: + + INTEGER, INTENT(IN ) :: JDATE ! Julian date, format YYYYDDD = 1000*Year + Day + INTEGER, INTENT( OUT) :: MNTH ! month (1...12) + INTEGER, INTENT( OUT) :: MDAY ! day-of-month (1...28,29,30,31) + + +C........... SCRATCH LOCAL VARIABLES: + + INTEGER IBIAS, IDATE, YEAR, DAY, L, J + + +C*********************************************************************** +C begin body of subroutine DAYMON + + IF ( JDATE .GT. -1000 ) THEN + IDATE = JDATE + IBIAS = 0 + ELSE + YEAR = -JDATE + YEAR = YEAR / 1000 + 1 + IBIAS = 2800000 * YEAR + IDATE = JDATE + IBIAS + END IF + +#ifdef IO_360 + DAY = MOD( IDATE, 1000 ) - 1 + MNTH = DAY / 30 + 1 + MDAY = MOD( DAY , 30 ) + 1 + RETURN +#endif + +#ifdef IO_365 + YEAR = JDATE / 1000 + DAY = MOD( IDATE, 1000 ) + J = MOD( DAY + 305, 365 ) + J = MOD( J, 153 ) / 61 + ( J / 153 ) * 2 + J + + MNTH = MOD( J / 31 + 2, 12 ) + 1 + MDAY = MOD( J, 31 ) + 1 + RETURN +#endif + + YEAR = JDATE / 1000 + DAY = MOD( IDATE, 1000 ) + IF ( YEAR .LE. 2 ) THEN !! "standard-year data + L = 365 + ELSE IF ( MOD( YEAR, 400 ) .EQ. 0 ) THEN + L = 366 + ELSE IF ( MOD( YEAR, 100 ) .EQ. 0 ) THEN + L = 365 + ELSE IF ( MOD( YEAR, 4 ) .EQ. 0 ) THEN + L = 366 + ELSE + L = 365 + END IF + + J = MOD( DAY + 305, L ) + J = MOD( J, 153 ) / 61 + ( J / 153 ) * 2 + J + + MNTH = MOD( J / 31 + 2, 12 ) + 1 + MDAY = MOD( J, 31 ) + 1 + + RETURN + + END SUBROUTINE DAYMON + diff --git a/src/io/ioapi/m3err.F b/src/io/ioapi/m3err.F new file mode 100644 index 0000000..3d653d0 --- /dev/null +++ b/src/io/ioapi/m3err.F @@ -0,0 +1,59 @@ + + SUBROUTINE M3ERR ( CALLER, JDATE, JTIME, MSGTXT, FATAL ) + +C*********************************************************************** +C Version "$Id: m3err.F 219 2015-08-17 18:05:54Z coats $" +C EDSS/Models-3 I/O API. +C Copyright (C) 1992-2002 MCNC and Carlie J. Coats, Jr., +C (C) 2003-2010 by Baron Advanced Meteorological Systems. +C Distributed under the GNU LESSER GENERAL PUBLIC LICENSE version 2.1 +C See file "LGPL.txt" for conditions of use. +C......................................................................... +C +C DEPRECATED!! +C Use M3EXIT(), instead. +C +C subroutine body starts at line 51 +C +C FUNCTION: Generate simple error messages for Models-3 core; +C terminate program execution via exit( 2 ) iff FATAL +C +C PRECONDITIONS REQUIRED: +C JDATE:JTIME represented as YYYYDDD:HHMMSS +C +C SUBROUTINES AND FUNCTIONS CALLED: DT2STR, INIT3, SHUT3 +C +C REVISION HISTORY: +C prototype 5/92 by CJC +C Revised 8/96 to close currently-open POSIX-OK Fortran units. +C Modified 1/97 by CJC to trim trailing blanks from MSGTXT +C Modified 10/98 by CJC: rewritten in terms of m3exit(), m3warn() +C*********************************************************************** + + IMPLICIT NONE + +C........... INCLUDES: + + INCLUDE 'IODECL3.EXT' + + +C........... ARGUMENTS and their descriptions: + + CHARACTER*(*), INTENT(IN ) :: CALLER ! name of the caller + INTEGER , INTENT(IN ) :: JDATE, JTIME ! model date&time for the error + CHARACTER*(*), INTENT(IN ) :: MSGTXT ! error message + LOGICAL , INTENT(IN ) :: FATAL ! terminate program iff TRUE + + +C*********************************************************************** +C begin body of subroutine M3ERR + + IF ( FATAL ) THEN + CALL M3EXIT( CALLER, JDATE, JTIME, MSGTXT, 2 ) + ELSE ! not endflag + CALL M3WARN( CALLER, JDATE, JTIME, MSGTXT ) + END IF + RETURN + + END SUBROUTINE M3ERR + diff --git a/src/io/ioapi/m3utilio.F90 b/src/io/ioapi/m3utilio.F90 index 398a0eb..5cfe18d 100644 --- a/src/io/ioapi/m3utilio.F90 +++ b/src/io/ioapi/m3utilio.F90 @@ -46,6 +46,14 @@ LOGICAL FUNCTION CURRSTEP ( JDATE, JTIME, & END FUNCTION CURRSTEP END INTERFACE + INTERFACE + SUBROUTINE DAYMON( JDATE, MNTH, MDAY ) + INTEGER, INTENT(IN ) :: JDATE ! Julian date, format YYYYDDD = 1000*Year + Day + INTEGER, INTENT( OUT) :: MNTH ! month (1...12) + INTEGER, INTENT( OUT) :: MDAY ! day-of-month (1...28,29,30,31) + END SUBROUTINE DAYMON + END INTERFACE + INTERFACE ! get file-description for FNAME LOGICAL FUNCTION DESC3( FNAME ) CHARACTER*(*), INTENT(IN ) :: FNAME ! file name @@ -58,13 +66,36 @@ CHARACTER*10 FUNCTION HHMMSS ( JTIME ) END FUNCTION HHMMSS END INTERFACE - INTERFACE - INTEGER FUNCTION INDEX1( NAME, N, NLIST ) +! INTERFACE +! INTEGER FUNCTION INDEX1( NAME, N, NLIST ) +! CHARACTER*(*), INTENT(IN ) :: NAME ! Character string being searched for +! INTEGER , INTENT(IN ) :: N ! Length of array to be searched +! CHARACTER*(*), INTENT(IN ) :: NLIST(*) ! array to be searched +! END FUNCTION INDEX1 +! END INTERFACE + + INTERFACE INDEXKEY + + INTEGER FUNCTION INDEX1( NAME, N, NLIST ) CHARACTER*(*), INTENT(IN ) :: NAME ! Character string being searched for INTEGER , INTENT(IN ) :: N ! Length of array to be searched CHARACTER*(*), INTENT(IN ) :: NLIST(*) ! array to be searched - END FUNCTION INDEX1 - END INTERFACE + END FUNCTION INDEX1 + + INTEGER FUNCTION INDEXINT1( IKEY, N, NLIST ) + INTEGER, INTENT(IN ) :: IKEY ! integer being searched for + INTEGER, INTENT(IN ) :: N ! Length of array to be searched + INTEGER, INTENT(IN ) :: NLIST(*) ! array to be searched + END FUNCTION INDEXINT1 + + INTEGER FUNCTION INDEXL1( LKEY, N, NLIST ) + INTEGER*8, INTENT(IN ) :: LKEY ! integer being searched for + INTEGER, INTENT(IN ) :: N ! Length of array to be searched + INTEGER*8, INTENT(IN ) :: NLIST(*) ! array to be searched + END FUNCTION INDEXL1 + + END INTERFACE !! indexkey + INTERFACE INTEGER FUNCTION JUNIT() @@ -114,6 +145,15 @@ SUBROUTINE M3MSG2( MESSAGE ) END SUBROUTINE M3MSG2 END INTERFACE + INTERFACE + SUBROUTINE M3EXIT( CALLER, JDATE, JTIME, MSGTXT, ISTAT ) + CHARACTER*(*), INTENT(IN ) :: CALLER ! name of the caller + INTEGER , INTENT(IN ) :: JDATE, JTIME ! model date&time for the error + CHARACTER*(*), INTENT(IN ) :: MSGTXT ! error message + INTEGER , INTENT(IN ) :: ISTAT ! exit status for program + END SUBROUTINE M3EXIT + END INTERFACE + INTERFACE READ3 MODULE PROCEDURE READ3_INTEGER MODULE PROCEDURE READ3_REAL @@ -172,7 +212,12 @@ LOGICAL FUNCTION OPEN3( FNAME, FSTATUS, PGNAME ) CHARACTER(LEN=*), INTENT(IN) :: FNAME INTEGER, INTENT(IN) :: FSTATUS CHARACTER(LEN=*), INTENT(IN) :: PGNAME + include SUBST_FILES_ID + OPEN3 = .TRUE. + IF (trim(fname) == trim(LUFRAC_CRO)) THEN + OPEN3 = .FALSE. + END IF END FUNCTION OPEN3 LOGICAL FUNCTION CLOSE3( FNAME ) diff --git a/src/model/src/PTMAP.F b/src/model/src/PTMAP.F deleted file mode 100644 index 8d651a6..0000000 --- a/src/model/src/PTMAP.F +++ /dev/null @@ -1,341 +0,0 @@ - MODULE PTMAP - - IMPLICIT NONE - - INTEGER, SAVE :: N_GSPC_EMIS = 0 ! number of gas species in diagnostic file - -C Species names from input file used for point source non-PM emissions mapping - INTEGER, ALLOCATABLE, SAVE :: PTEM_MAP( : ) - INTEGER, ALLOCATABLE, SAVE :: SPC_PTEM_MAP( : ) - REAL, ALLOCATABLE, SAVE :: SPC_PTEM_FAC( : ) - -C Species names from input file used for point source PM emissions mapping - INTEGER, ALLOCATABLE, SAVE :: PTPM_MAP( : ) - -C Mapping for point source non-PM emissions (maps only GC, NR and TR) - INTEGER, SAVE :: N_SPC_PTEM = 0 ! merged no. of unique species - INTEGER, SAVE :: N_SPC_PTPM = 0 ! merged no. of unique species for PM - - PRIVATE - - PUBLIC :: N_GSPC_EMIS, N_SPC_PTEM, N_SPC_PTPM - PUBLIC :: SPC_PTEM_FAC, SPC_PTEM_MAP - PUBLIC :: PTEM_MAP, PTPM_MAP - - PUBLIC :: PTMAP_INIT - PUBLIC :: PTMAP_TYPE_INIT - - CONTAINS - - FUNCTION PTMAP_INIT ( ) RESULT ( SUCCESS ) - - USE CGRID_SPCS ! CGRID mechanism species - USE AERO_DATA, ONLY : N_EMIS_PM, PMEM_MAP_NAME ! defines aerosol species - USE AQM_EMIS_MOD, ONLY : AQM_EMIS_GET, AQM_EMIS_ISPRESENT, - & AQM_INTERNAL_EMIS_TYPE - USE UTILIO_DEFN - - IMPLICIT NONE - - LOGICAL SUCCESS - - INTEGER I, IDX, IOS, IPT, NPT - INTEGER N, N_GAS_EMIS, NSPC, NSPC1, NSPC2, NSPC3 - INTEGER S, S_OFFSET, V - - INTEGER, ALLOCATABLE :: MAP( : ) - CHARACTER( 16 ), ALLOCATABLE :: EMSPC( : ) - - TYPE PT_EMIS_TYPE - TYPE( AQM_INTERNAL_EMIS_TYPE ), POINTER :: EM - END TYPE PT_EMIS_TYPE - - TYPE( PT_EMIS_TYPE ) :: PT( 2 ) - - CHARACTER( * ), PARAMETER :: PNAME = 'PTMAP_INIT' ! procedure name - - CHARACTER( * ), PARAMETER :: ETYPES( 2 ) = - & (/ 'gbbepx ', 'point-source' /) - -C----------------------------------------------------------------------- - - SUCCESS = .TRUE. - -C check if emissions are being provided - - NPT = SIZE( ETYPES ) - - I = 0 - NSPC = 0 - DO IPT = 1, NPT - IF ( AQM_EMIS_ISPRESENT( ETYPES( IPT ) ) ) THEN - I = I + 1 - PT( I ) % EM => AQM_EMIS_GET( ETYPES( IPT ) ) - NSPC = NSPC + SIZE( PT( I ) % EM % TABLE, DIM=1 ) - END IF - END DO - - IF ( I == 0 ) RETURN - - NPT = I - - ALLOCATE ( EMSPC( NSPC ), STAT = IOS ) - CALL CHECKMEM( IOS, 'EMSPC', PNAME ) - MAP = 0 - - EMSPC = "" - - ! -- merge emission tables and remove duplicates - - I = 0 - DO IPT = 1, NPT - DO N = 1, SIZE( PT( IPT ) % EM % TABLE, DIM=1 ) - IDX = INDEX1( TRIM( PT( IPT ) % EM % TABLE( N, 1 ) ), NSPC, EMSPC ) - IF ( IDX .LE. 0 ) THEN - I = I + 1 - EMSPC( I ) = PT( IPT ) % EM % TABLE( N, 1 ) - END IF - END DO - END DO - - NSPC = I - -C compute emission offsets - - NSPC1 = N_GC_EMIS - NSPC2 = NSPC1 + N_AE_EMIS - NSPC3 = NSPC2 + N_NR_EMIS - -C create auxiliary arrays mapping fire emission species to CMAQ gas and aerosol species - - ALLOCATE( MAP( NSPC ), STAT = IOS ) - CALL CHECKMEM( IOS, 'MAP', PNAME ) - MAP = 0 - -C ... gas species ... - - NSPC1 = N_GC_EMIS - NSPC2 = NSPC1 + N_AE_EMIS - NSPC3 = NSPC2 + N_NR_EMIS - - N_GAS_EMIS = N_GC_EMIS + N_NR_EMIS + N_TR_EMIS - - ALLOCATE( SPC_PTEM_FAC( N_GAS_EMIS ), STAT = IOS ) - CALL CHECKMEM( IOS, 'SPC_PTEM_FAC', PNAME ) - SPC_PTEM_FAC = 1.0 - - ALLOCATE( SPC_PTEM_MAP( N_GAS_EMIS ), STAT = IOS ) - CALL CHECKMEM( IOS, 'SPC_PTEM_MAP', PNAME ) - SPC_PTEM_MAP = -1 - - ALLOCATE( PTEM_MAP( N_GAS_EMIS ), STAT = IOS ) - CALL CHECKMEM( IOS, 'PTEM_MAP', PNAME ) - PTEM_MAP = -1 - - N = 0 - - S_OFFSET = 0 - DO S = 1, N_GC_EMIS - IDX = INDEX1( GC_EMIS( S ), NSPC, EMSPC ) - IF ( IDX .GT. 0 ) THEN - N = N + 1 - PTEM_MAP ( S ) = N - SPC_PTEM_MAP( S ) = S - SPC_PTEM_FAC( S ) = GC_EMIS_FAC( S ) - END IF - END DO - - S_OFFSET = N_GC_EMIS - DO S = 1, N_NR_EMIS - IDX = INDEX1( NR_EMIS( S ), NSPC, EMSPC ) - IF ( IDX .GT. 0 ) THEN - N = N + 1 - V = S + S_OFFSET - PTEM_MAP ( V ) = N - SPC_PTEM_MAP( V ) = S + NSPC2 - SPC_PTEM_FAC( V ) = NR_EMIS_FAC( S ) - END IF - END DO - - S_OFFSET = S_OFFSET + N_NR_EMIS - DO S = 1, N_TR_EMIS - IDX = INDEX1( TR_EMIS( S ), NSPC, EMSPC ) - IF ( IDX .GT. 0 ) THEN - N = N + 1 - V = S + S_OFFSET - PTEM_MAP ( V ) = N - SPC_PTEM_MAP( V ) = S + NSPC3 - SPC_PTEM_FAC( V ) = TR_EMIS_FAC( S ) - END IF - END DO - - N_SPC_PTEM = N - N_GSPC_EMIS = N_GAS_EMIS - -C ... aerosol species ... - - N = 0 - DO S = 1, NSPC - IDX = INDEX1( EMSPC( S ), N_EMIS_PM, PMEM_MAP_NAME ) - IF ( IDX .GT. 0 ) THEN - N = N + 1 - MAP( N ) = IDX - END IF - END DO - - N_SPC_PTPM = N - - ALLOCATE( PTPM_MAP( N_SPC_PTPM ), STAT = IOS ) - CALL CHECKMEM( IOS, 'PTPM_MAP', PNAME ) - PTPM_MAP = MAP( 1:N_SPC_PTPM ) - - DEALLOCATE( EMSPC, MAP ) - - END FUNCTION PTMAP_INIT - - - FUNCTION PTMAP_TYPE_INIT( EM, - & N_GSPC_TYPE_EMIS, N_SPC_TYPE_PTPM, - & PTEM_TYPE_MAP, PTPM_TYPE_MAP, - & SPC_PTEM_TYPE_FAC, SPC_PTEM_TYPE_MAP, - & ELABEL ) - & RESULT ( SUCCESS ) - - - USE CGRID_SPCS ! CGRID mechanism species - USE AERO_DATA, ONLY : N_EMIS_PM, PMEM_MAP_NAME ! defines aerosol species - USE AQM_EMIS_MOD, ONLY : AQM_INTERNAL_EMIS_TYPE - USE UTILIO_DEFN - - IMPLICIT NONE - - ! -- arguments - - TYPE( AQM_INTERNAL_EMIS_TYPE ), POINTER :: EM - INTEGER, INTENT( OUT ) :: N_GSPC_TYPE_EMIS - INTEGER, INTENT( OUT ) :: N_SPC_TYPE_PTPM - INTEGER, ALLOCATABLE, INTENT( OUT ) :: PTEM_TYPE_MAP( : ) - INTEGER, ALLOCATABLE, INTENT( OUT ) :: PTPM_TYPE_MAP( : ) - INTEGER, ALLOCATABLE, INTENT( OUT ) :: SPC_PTEM_TYPE_FAC( : ) - INTEGER, ALLOCATABLE, INTENT( OUT ) :: SPC_PTEM_TYPE_MAP( : ) - CHARACTER( * ), INTENT( IN ) :: ELABEL - - ! -- local variables - - LOGICAL SUCCESS - - INTEGER IOS - - INTEGER IDX - INTEGER N, N_GAS_EMIS, NSPC, NSPC1, NSPC2, NSPC3 - INTEGER S, S_OFFSET, V - - INTEGER, ALLOCATABLE :: MAP( : ) - - CHARACTER( * ), PARAMETER :: PNAME = 'PTEM_TYPE_MAP' ! procedure name - -C----------------------------------------------------------------------- - - SUCCESS = .TRUE. - - N_GSPC_TYPE_EMIS = 0 - N_SPC_TYPE_PTPM = 0 - - IF ( .NOT.ASSOCIATED( EM ) ) RETURN - -C compute emission offsets - - NSPC1 = N_GC_EMIS - NSPC2 = NSPC1 + N_AE_EMIS - NSPC3 = NSPC2 + N_NR_EMIS - -C create auxiliary arrays mapping fire emission species to CMAQ gas and aerosol species - - NSPC = SIZE( EM % TABLE, DIM=1 ) - - ALLOCATE( MAP( NSPC ), STAT = IOS ) - CALL CHECKMEM( IOS, 'MAP: ' // TRIM( ELABEL ), PNAME ) - MAP = 0 - -C ... gas species ... - - NSPC1 = N_GC_EMIS - NSPC2 = NSPC1 + N_AE_EMIS - NSPC3 = NSPC2 + N_NR_EMIS - - N_GAS_EMIS = N_GC_EMIS + N_NR_EMIS + N_TR_EMIS - - ALLOCATE( SPC_PTEM_TYPE_FAC( N_GAS_EMIS ), STAT = IOS ) - CALL CHECKMEM( IOS, 'SPC_PTEM_TYPE_FAC: ' // TRIM( ELABEL ), PNAME ) - SPC_PTEM_TYPE_FAC = 1.0 - - ALLOCATE( SPC_PTEM_TYPE_MAP( N_GAS_EMIS ), STAT = IOS ) - CALL CHECKMEM( IOS, 'SPC_PTEM_TYPE_MAP: ' // TRIM( ELABEL ), PNAME ) - SPC_PTEM_TYPE_MAP = -1 - - ALLOCATE( PTEM_TYPE_MAP( N_GAS_EMIS ), STAT = IOS ) - CALL CHECKMEM( IOS, 'PTEM_TYPE_MAP: ' // TRIM( ELABEL ), PNAME ) - PTEM_TYPE_MAP = -1 - - N = 0 - - S_OFFSET = 0 - DO S = 1, N_GC_EMIS - IDX = INDEX1( GC_EMIS( S ), NSPC, EM % TABLE( 1, 1 ) ) - IF ( IDX .GT. 0 ) THEN - N = N + 1 - PTEM_TYPE_MAP ( S ) = IDX - SPC_PTEM_TYPE_MAP( S ) = S - SPC_PTEM_TYPE_FAC( S ) = GC_EMIS_FAC( S ) - END IF - END DO - - S_OFFSET = N_GC_EMIS - DO S = 1, N_NR_EMIS - IDX = INDEX1( NR_EMIS( S ), NSPC, EM % TABLE( 1, 1 ) ) - IF ( IDX .GT. 0 ) THEN - N = N + 1 - V = S + S_OFFSET - PTEM_TYPE_MAP ( V ) = IDX - SPC_PTEM_TYPE_MAP( V ) = S + NSPC2 - SPC_PTEM_TYPE_FAC( V ) = NR_EMIS_FAC( S ) - END IF - END DO - - S_OFFSET = S_OFFSET + N_NR_EMIS - DO S = 1, N_TR_EMIS - IDX = INDEX1( TR_EMIS( S ), NSPC, EM % TABLE( 1, 1 ) ) - IF ( IDX .GT. 0 ) THEN - N = N + 1 - V = S + S_OFFSET - PTEM_TYPE_MAP ( V ) = IDX - SPC_PTEM_TYPE_MAP( V ) = S + NSPC3 - SPC_PTEM_TYPE_FAC( V ) = TR_EMIS_FAC( S ) - END IF - END DO - - N_GSPC_TYPE_EMIS = N - -C ... aerosol species ... - - N = 0 - DO S = 1, N_SPC_PTPM - V = PTPM_MAP( S ) - IDX = INDEX1( PMEM_MAP_NAME( V ), NSPC, EM % TABLE( 1, 1 ) ) - IF ( IDX .GT. 0 ) THEN - N = N + 1 - MAP( N ) = S - END IF - END DO - - N_SPC_TYPE_PTPM = N - - ALLOCATE( PTPM_TYPE_MAP( N_SPC_TYPE_PTPM ), STAT = IOS ) - CALL CHECKMEM( IOS, 'PTPM_TYPE_MAP: '// TRIM( ELABEL ), PNAME ) - PTPM_TYPE_MAP = MAP( 1:N_SPC_TYPE_PTPM ) - - DEALLOCATE( MAP ) - - END FUNCTION PTMAP_TYPE_INIT - - END MODULE PTMAP diff --git a/src/shr/aqm_config_mod.F90 b/src/shr/aqm_config_mod.F90 index 4147b7d..3aca5ce 100644 --- a/src/shr/aqm_config_mod.F90 +++ b/src/shr/aqm_config_mod.F90 @@ -23,6 +23,9 @@ module aqm_config_mod character(len=AQM_MAXSTR) :: csqy_data = "" character(len=AQM_MAXSTR) :: optics_data = "" character(len=AQM_MAXSTR) :: omi = "" + character(len=AQM_MAXSTR) :: desid_chem_ctrl = "" + character(len=AQM_MAXSTR) :: desid_ctrl = "" + character(len=AQM_MAXSTR) :: misc_ctrl = "" character(len=AQM_MAXSTR) :: mp_map = "" character(len=AQM_MAXSTR) :: ctm_stdout = "" integer :: dy_map_beg = 0 @@ -140,6 +143,33 @@ subroutine aqm_config_read(model, config, rc) rcToReturn=rc)) & return ! bail out + ! -- DESID_CHEM_CTRL + call ESMF_ConfigGetAttribute(cf, config % desid_chem_ctrl, & + label="desid_chem_ctrl:", rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + + ! -- DESID_CTRL + call ESMF_ConfigGetAttribute(cf, config % desid_ctrl, & + label="desid_ctrl:", rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + + ! -- MISC_CTRL for ELMO + call ESMF_ConfigGetAttribute(cf, config % misc_ctrl, & + label="misc_ctrl:", rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + ! -- read run settings call ESMF_ConfigGetAttribute(cf, config % run_aero, & label="run_aerosol:", default=.true., rc=localrc) @@ -538,6 +568,27 @@ subroutine aqm_config_log(config, name, rc) file=__FILE__, & rcToReturn=rc)) & return ! bail out + call ESMF_LogWrite(trim(name) // ": config: read: desid_chem_ctrl: " & + // config % desid_chem_ctrl, ESMF_LOGMSG_INFO, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + call ESMF_LogWrite(trim(name) // ": config: read: desid_ctrl: " & + // config % desid_ctrl, ESMF_LOGMSG_INFO, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out + call ESMF_LogWrite(trim(name) // ": config: read: misc_ctrl: " & + // config % misc_ctrl, ESMF_LOGMSG_INFO, rc=localrc) + if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, & + line=__LINE__, & + file=__FILE__, & + rcToReturn=rc)) & + return ! bail out if (config % ctm_aod) then call ESMF_LogWrite(trim(name) // ": config: read: ctm_aod: true", & ESMF_LOGMSG_INFO, rc=localrc) diff --git a/src/shr/aqm_emis_mod.F90 b/src/shr/aqm_emis_mod.F90 index 711264d..7227aed 100644 --- a/src/shr/aqm_emis_mod.F90 +++ b/src/shr/aqm_emis_mod.F90 @@ -1512,18 +1512,20 @@ logical function aqm_emis_ispresent(etype) end function aqm_emis_ispresent - - subroutine aqm_emis_desc( etype, nlays, nvars, vnames, units ) + !Add number of points for fire and point source + subroutine aqm_emis_desc( etype, nlays, nvars, vnames, units, npoints ) character(len=*), intent(in) :: etype integer, optional, intent(out) :: nlays integer, optional, intent(out) :: nvars character(len=16), optional, intent(out) :: vnames(:) character(len=16), optional, intent(out) :: units(:) + integer, optional, intent(out) :: npoints ! -- local variables integer :: localrc integer :: item, nsrc type(aqm_internal_emis_type), pointer :: em + type(aqm_state_type), pointer :: stateIn ! -- begin ! -- get emission data @@ -1549,11 +1551,26 @@ subroutine aqm_emis_desc( etype, nlays, nvars, vnames, units ) if (present(nvars)) nvars = nsrc if (present(vnames)) vnames( 1:nsrc ) = em % table( 1:nsrc, 1 ) if (present(units)) units ( 1:nsrc ) = em % table( 1:nsrc, 2 ) + !add npoints here ;treat grid data as points + if (present(npoints)) then + if (etype == 'gbbepx') then + nullify(stateIn) + 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 + npoints = size (stateIn % area) + else if (etype == 'point-source') then + npoints = size (em % ijmap) + else + npoints = 0 + end if + end if else if (present(nlays)) nlays = 0 if (present(nvars)) nvars = 0 if (present(vnames)) vnames = "" if (present(units)) units = "" + if (present(npoints)) npoints = 0 end if end subroutine aqm_emis_desc @@ -1617,7 +1634,7 @@ subroutine aqm_emis_grd_read(em, spcname, buffer, localDe, rc) if (abs(fptr(i,j)) < emAccept) then buffer(k) = buffer(k) & + em % factors(item) * fptr(i,j) / stateIn % area(i,j) & - / stateIn % area(i,j) + / stateIn % area(i,j) end if end do end do @@ -1720,7 +1737,7 @@ subroutine aqm_emis_pts_read(em, spcname, buffer, ip, jp, ijmap, localDe, rc) j = em % jp(n) buffer(n) = buffer(n) & + em % factors(item) * em % rates(item) % values(n) / stateIn % area(i,j) & - / stateIn % area(i,j) + / stateIn % area(i,j) end do case (0) ! -- emissions are totals over each grid cell @@ -1736,7 +1753,7 @@ subroutine aqm_emis_pts_read(em, spcname, buffer, ip, jp, ijmap, localDe, rc) do m = 1, size(em % ijmap) n = em % ijmap(m) buffer(n) = buffer(n) & - + em % factors(item) * em % rates(item) % values(n) + + em % factors(item) * em % rates(item) % values(n) end do case default ! -- this case should never occur diff --git a/src/shr/aqm_methods.F90 b/src/shr/aqm_methods.F90 index 7412e7d..49c1278 100644 --- a/src/shr/aqm_methods.F90 +++ b/src/shr/aqm_methods.F90 @@ -86,7 +86,6 @@ LOGICAL FUNCTION DESC3( FNAME ) IMPLICIT NONE CHARACTER(LEN=*), INTENT(IN) :: FNAME - CHARACTER(LEN=len(FNAME)) :: FNAME_TRIM !(Wei Li) INCLUDE SUBST_FILES_ID @@ -107,8 +106,7 @@ LOGICAL FUNCTION DESC3( FNAME ) STIME3D = 0 TSTEP3D = 0 - FNAME_TRIM = TRIM(FNAME_TRIM) - !!Replace INIT_GASC,AERO,NONR,TRAC to INIT_CONC_1 (Wei Li) + !!Replace INIT_GASC,AERO,NONR,TRAC to INIT_CONC_1 IF ( (TRIM(FNAME) .EQ. TRIM(INIT_CONC_1)) ) THEN ! -- Input initial background values for the following species @@ -126,21 +124,31 @@ LOGICAL FUNCTION DESC3( FNAME ) call aqm_emis_desc("biogenic", NLAYS3D, NVARS3D, VNAME3D, UNITS3D) -! EMIS_1 is not used anymore. Change to other env variables. (Wei Li) - ELSE IF ( ( (FNAME_TRIM(1:8) .EQ. 'GR_EMIS_') .AND. (len(FNAME_TRIM) .EQ. 11 )) .OR. & - ( (FNAME_TRIM(1:9) .EQ. 'STK_EMIS_').AND. (len(FNAME_TRIM) .EQ. 12 )) ) THEN - +! EMIS_1 is not used anymore. Change to other env variables. + ELSE IF ( TRIM( FNAME ) .EQ. 'GR_EMIS_001' ) THEN NLAYS3D = 0 - - call aqm_emis_desc("gbbepx", NLAYS=EMLAYS) + call aqm_emis_desc("anthropogenic", NLAYS=EMLAYS, NVARS=NVARS3D, VNAMES=VNAME3D, UNITS=UNITS3D) NLAYS3D = MAX(EMLAYS, NLAYS3D) + !This is to add the TSTEP3D for ratio calculation in interpolate_var function + call aqm_model_get(config=config, rc=localrc) + if (aqm_rc_check(localrc, msg="Failure to retrieve model input state", & + file=__FILE__, line=__LINE__)) return + + SDATE3D = config % ctm_stdate + STIME3D = config % ctm_sttime + TSTEP3D = config % ctm_tstep - call aqm_emis_desc("point-source", NLAYS=EMLAYS) + ELSE IF ( TRIM( FNAME ) .EQ. 'STK_EMIS_001' ) THEN !fire stream + NLAYS3D = 0 + call aqm_emis_desc("gbbepx", NLAYS=EMLAYS,NVARS=NVARS3D, VNAMES=VNAME3D, UNITS=UNITS3D, NPOINTS=NROWS3D) NLAYS3D = MAX(EMLAYS, NLAYS3D) - call aqm_emis_desc("anthropogenic", NLAYS=EMLAYS, NVARS=NVARS3D, VNAMES=VNAME3D, UNITS=UNITS3D) + ELSE IF ( TRIM( FNAME ) .EQ. 'STK_EMIS_002' ) THEN !STKS stream + NLAYS3D = 0 + call aqm_emis_desc("point-source", NLAYS=EMLAYS,NVARS=NVARS3D, VNAMES=VNAME3D, UNITS=UNITS3D,NPOINTS=NROWS3D) NLAYS3D = MAX(EMLAYS, NLAYS3D) + ELSE IF ( TRIM( FNAME ) .EQ. TRIM( GRID_DOT_2D ) ) THEN NVARS3D = 1 VNAME3D( 1:NVARS3D ) = & @@ -221,7 +229,7 @@ LOGICAL FUNCTION DESC3( FNAME ) SDATE3D = config % ctm_stdate STIME3D = config % ctm_sttime TSTEP3D = config % ctm_tstep - + ELSE IF ( TRIM( FNAME ) .EQ. TRIM( MET_CRO_3D ) ) THEN CALL aqm_model_domain_get(nl=NLAYS3D, rc=localrc) @@ -366,14 +374,14 @@ logical function envyn(name, description, defaultval, status) envyn = .false. em => aqm_emis_get("biogenic") if (associated(em)) envyn = (trim(em % period) == "summer") - ! case ('CTM_AOD') - ! envyn = config % ctm_aod + case ('CTM_AOD') + envyn = config % ctm_aod case ('CTM_BIOGEMIS') envyn = aqm_emis_ispresent("biogenic") case ('CTM_DEPVFILE') envyn = config % ctm_depvfile - ! case ('CTM_PMDIAG') - ! envyn = config % ctm_pmdiag + case ('CTM_PMDIAG') + envyn = config % ctm_pmdiag case ('CTM_PHOTODIAG') envyn = config % ctm_photodiag case ('CTM_PT3DEMIS') @@ -569,8 +577,11 @@ subroutine nameval(name, eqname) integer :: deCount, localrc type(aqm_config_type), pointer :: config type(aqm_internal_emis_type), pointer :: em + !INTEGER, EXTERNAL :: SETUP_LOGDEV + !INTEGER, SAVE :: LOGDEV ! -- begin + !LOGDEV = SETUP_LOGDEV() eqname = "" nullify(config) @@ -591,6 +602,12 @@ subroutine nameval(name, eqname) eqname = config % tr_matrix_nml case ('CSQY_DATA') eqname = config % csqy_data + case ('MISC_CTRL') + eqname = config % misc_ctrl + case ('DESID_CTRL') + eqname = config % desid_ctrl + case ('DESID_CHEM_CTRL') + eqname = config % desid_chem_ctrl case ('GSPRO') nullify(em) em => aqm_emis_get("biogenic") @@ -621,7 +638,6 @@ logical function interpx( fname, vname, pname, & implicit none character(len=*), intent(in) :: fname, vname, pname - CHARACTER(LEN=len(fname)) :: FNAME_TRIM !(Wei Li) integer, intent(in) :: col0, col1, row0, row1, lay0, lay1 integer, intent(in) :: jdate, jtime real, intent(out) :: buffer(*) @@ -646,7 +662,6 @@ logical function interpx( fname, vname, pname, & ! -- begin interpx = .false. - FNAME_TRIM = TRIM(fname) !(Wei Li) lbuf = (col1-col0+1) * (row1-row0+1) * (lay1-lay0+1) buffer(1:lbuf) = 0. @@ -861,14 +876,15 @@ logical function interpx( fname, vname, pname, & return end select -! EMIS_1 is not used anymore. Change to other env variables. (Wei Li) - else if ( ( (FNAME_TRIM(1:8) .EQ. 'GR_EMIS_') .AND. (len(FNAME_TRIM) .EQ. 11 )) ) then +! EMIS_1 is not used anymore. Change to other env variables. + else if ( trim(fname) .EQ. 'GR_EMIS_001') then ! -- read in emissions call aqm_emis_read("anthropogenic", vname, buffer, rc=localrc) if (aqm_rc_test((localrc /= 0), & msg="Failure to read emissions for " // vname, & file=__FILE__, line=__LINE__)) return + else if (trim(fname) == trim(MET_CRO_3D)) then call aqm_model_get(config=config, stateIn=stateIn, rc=localrc) @@ -1050,12 +1066,15 @@ LOGICAL FUNCTION XTRACT3 ( FNAME, VNAME, & LAY0, LAY1, ROW0, ROW1, COL0, COL1, & JDATE, JTIME, BUFFER ) - use aqm_types_mod, only : AQM_KIND_R4 - use aqm_model_mod, only : aqm_state_type, aqm_model_get + use aqm_types_mod, only : AQM_KIND_R4, AQM_KIND_R8, AQM_MAXSTR + use aqm_model_mod, only : aqm_config_type,aqm_state_type, & + aqm_model_get, aqm_model_domain_get use aqm_rc_mod, only : aqm_rc_check, aqm_rc_test use aqm_const_mod, only : con_mr2ppm_o3, thrs_p_strato use aqm_emis_mod, only : aqm_emis_read use aqm_config_mod + use aqm_const_mod, only : eps1, grav, onebg, rdgas + USE M3UTILIO, ONLY : M3MESG implicit none @@ -1074,20 +1093,32 @@ LOGICAL FUNCTION XTRACT3 ( FNAME, VNAME, & ! -- local variables integer :: localrc - integer :: c, r, l, k, lbuf, lu_index + integer :: c, r, l, k, n, nl, lbuf, lu_index type(aqm_config_type), pointer :: config type(aqm_state_type), pointer :: stateIn + !add from interpx + logical :: set_non_neg + character(len=16) :: varname + character(len=AQM_MAXSTR) :: msgString + real(AQM_KIND_R8), dimension(:,:), pointer :: lat, lon + real(AQM_KIND_R8), dimension(:,:), pointer :: p2d + real(AQM_KIND_R8), dimension(:,:,:), pointer :: p3d + include SUBST_FILES_ID + logical, parameter :: debug = .true. - ! -- begin + ! -- begin + nullify(p2d) + nullify(p3d) nullify(config) nullify(stateIn) + set_non_neg = .false. lbuf = (LAY1-LAY0+1)*(ROW1-ROW0+1)*(COL1-COL0+1) BUFFER(1:lbuf) = 0. - XTRACT3 = .TRUE. + XTRACT3 = .FALSE. IF (TRIM(FNAME) == TRIM(GRID_CRO_2D)) THEN @@ -1095,12 +1126,30 @@ LOGICAL FUNCTION XTRACT3 ( FNAME, VNAME, & if (aqm_rc_check(localrc, msg="Failure to retrieve model input state", & file=__FILE__, line=__LINE__)) return + call aqm_model_domain_get(lon=lon, lat=lat, rc=localrc) + if (aqm_rc_check(localrc, msg="Failure to retrieve grid coordinates", & + file=__FILE__, line=__LINE__)) return + +! if (vname(1:7) == 'LUFRAC_') then +! if (aqm_rc_test((LAY0.NE.1).OR.(LAY1.NE.1), & +! msg=TRIM(VNAME)//" is 2D. LAY0 and LAY1 must be 1", & +! file=__FILE__, line=__LINE__)) return +! lu_index = 0 +! read(vname(8:), *, iostat=localrc) lu_index +! if (aqm_rc_test(localrc /= 0, msg="Failure to identify LU_INDEX", & +! file=__FILE__, line=__LINE__)) return +! k = 0 +! do r = row0, row1 +! do c = col0, col1 +! k = k + 1 +! if (int(stateIn % stype(c,r)) == lu_index) buffer(k) = 1.0 +! end do +! end do +! end if + if (vname(1:7) == 'LUFRAC_') then - if (aqm_rc_test((LAY0.NE.1).OR.(LAY1.NE.1), & - msg=TRIM(VNAME)//" is 2D. LAY0 and LAY1 must be 1", & - file=__FILE__, line=__LINE__)) return lu_index = 0 - read(vname(8:), *, iostat=localrc) lu_index + read(vname(8:9), *, iostat=localrc) lu_index if (aqm_rc_test(localrc /= 0, msg="Failure to identify LU_INDEX", & file=__FILE__, line=__LINE__)) return k = 0 @@ -1110,8 +1159,327 @@ LOGICAL FUNCTION XTRACT3 ( FNAME, VNAME, & if (int(stateIn % stype(c,r)) == lu_index) buffer(k) = 1.0 end do end do + else + select case (trim(vname)) + case ('HT') + p2d => stateIn % ht + case ('LAT') + p2d => lat + case ('LON') + k = 0 + do r = row0, row1 + do c = col0, col1 + k = k + 1 + if (lon(c,r) > 180.) then + buffer(k) = lon(c,r) - 360. + else + buffer(k) = lon(c,r) + end if + end do + end do + case ('LWMASK') + k = 0 + do r = row0, row1 + do c = col0, col1 + k = k + 1 + buffer(k) = stateIn % slmsk(c,r) + if (nint(buffer(k)) == 2) buffer(k) = 0. ! set sea ice points as water + end do + end do + case ('MSFX2') + buffer(1:lbuf) = 1. + case ('PURB') + case default + return + end select end if + ELSE IF (trim(fname) == trim(MET_CRO_2D)) THEN + + 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 + + select case (trim(vname)) + case ("HFX") + p2d => stateIn % hfx + case ("LAI") + p2d => stateIn % xlai + case ("LH") + p2d => stateIn % lh + case ("PRSFC") + p2d => stateIn % psfc + case ("PBL") + p2d => stateIn % hpbl + case ("Q2") + p2d => stateIn % q2m + case ("RADYNI") + p2d => stateIn % cmm + case ("RSTOMI") + k = 0 + do r = row0, row1 + do c = col0, col1 + k = k + 1 + if ( stateIn % rc(c,r) /= 0.0 ) buffer(k) = 1.0 / stateIn % rc(c,r) + end do + end do + case ("RA") + k = 0 + do r = row0, row1 + do c = col0, col1 + k = k + 1 + buffer(k) = sqrt(stateIn % uwind(c,r,1) * stateIn % uwind(c,r,1) + & + stateIn % vwind(c,r,1) * stateIn % vwind(c,r,1)) / & + ( stateIn % ustar(c,r) * stateIn % ustar(c,r) ) + end do + end do + case ("RS") + p2d => stateIn % rc + case ("RC") + k = 0 + do r = row0, row1 + do c = col0, col1 + k = k + 1 + buffer(k) = 100. * stateIn % rainc(c,r) + end do + end do + case ("RGRND") + p2d => stateIn % rgrnd + case ("RN") + k = 0 + do r = row0, row1 + do c = col0, col1 + k = k + 1 + buffer(k) = 100. * (stateIn % rain(c,r) - stateIn % rainc(c,r)) + end do + end do + set_non_neg = .true. + case ("SEAICE") + p2d => stateIn % fice + case ("SLTYP") + p2d => stateIn % stype + case ("SNOCOV") + p2d => stateIn % sncov + case ("SOIM1") + p2d => stateIn % smois(:,:,1) + case ("SOIM2") + p2d => stateIn % smois(:,:,2) + case ("SOIT1") + p2d => stateIn % stemp(:,:,1) + case ("SOIT2") + p2d => stateIn % stemp(:,:,2) + case ("TEMPG") + p2d => stateIn % tsfc + case ("TEMP2") + p2d => stateIn % t2m + case ("USTAR") + p2d => stateIn % ustar + case ("VEG") + p2d => stateIn % vfrac + case ("WR") + p2d => stateIn % wr + case ("WSPD10") + k = 0 + do r = row0, row1 + do c = col0, col1 + k = k + 1 + buffer(k) = sqrt(stateIn % u10m(c,r) * stateIn % u10m(c,r) & + + stateIn % v10m(c,r) * stateIn % v10m(c,r)) + end do + end do + case ("ZRUF") + k = 0 + do r = row0, row1 + do c = col0, col1 + k = k + 1 + buffer(k) = 0.01 * stateIn % zorl(c,r) + end do + end do + case ("CLAYF","DRAG","SANDF","UTHR") + ! -- 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 + + ELSE IF (trim(fname) == trim(OCEAN_1)) THEN + select case (trim(vname)) + case ("OPEN") + 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 + ! -- set to complement to land mask + k = 0 + do r = row0, row1 + do c = col0, col1 + k = k + 1 + buffer(k) = 1.0 - stateIn % slmsk(c,r) + if (nint(stateIn % slmsk(c,r)) == 2) buffer(k) = 1.0 ! set sea ice points as water + end do + end do + case ("SURF") + ! -- zero + case default + return + end select + + ! EMIS_1 is not used anymore. Change to other env variables. + ELSE IF ( trim(fname) .EQ. 'GR_EMIS_001') then + ! -- read in emissions + call aqm_emis_read("anthropogenic", vname, buffer, rc=localrc) + if (aqm_rc_test((localrc /= 0), & + msg="Failure to read emissions for " // vname, & + file=__FILE__, line=__LINE__)) return + + ELSE IF (trim(fname) == trim(MET_CRO_3D)) THEN + + 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 + + select case (trim(vname)) + case ("JACOBF") + call aqm_model_domain_get(nl=nl, rc=localrc) + if (aqm_rc_check(localrc, msg="Failure to retrieve model coordinates", & + file=__FILE__, line=__LINE__)) return + k = 0 + do l = lay0, lay1 + n = min(l + 1, nl) + do r = row0, row1 + do c = col0, col1 + k = k + 1 + buffer(k) = onebg * (stateIn % phil(c,r,n) - stateIn % phil(c,r,n-1)) + end do + end do + end do + case ("JACOBM") + k = 0 + do l = lay0, lay1 + do r = row0, row1 + do c = col0, col1 + k = k + 1 + buffer(k) = onebg * (stateIn % phii(c,r,l+1) - stateIn % phii(c,r,l)) + end do + end do + end do + case ("DENS") + k = 0 + do l = lay0, lay1 + do r = row0, row1 + do c = col0, col1 + k = k + 1 + buffer(k) = stateIn % prl(c,r,l) / ( rdgas * stateIn % temp(c,r,l) ) + end do + end do + end do + case ("DENSA_J") + k = 0 + do l = lay0, lay1 + do r = row0, row1 + do c = col0, col1 + k = k + 1 + ! -- rho + buffer(k) = stateIn % prl(c,r,l) / ( rdgas * stateIn % temp(c,r,l) ) + ! -- Jacobian + buffer(k) = buffer(k) & + * onebg * (stateIn % phii(c,r,l+1) - stateIn % phii(c,r,l)) + end do + end do + end do + case ("PRES") + p3d => stateIn % prl + case ("PRESF") + p3d => stateIn % pri + case ("CFRAC_3D") + p3d => stateIn % cldfl + case ("PV") + buffer(1:lbuf) = 1.0 + case ("QV") + p3d => stateIn % tr(:,:,:,config % species % p_atm_qv) + set_non_neg = .true. + case ("QC") + p3d => stateIn % tr(:,:,:,config % species % p_atm_qc) + set_non_neg = .true. + case ("QR") + if (config % species % p_atm_qr > 0) then + p3d => stateIn % tr(:,:,:,config % species % p_atm_qr) + set_non_neg = .true. + end if + case ("QI") + if (config % species % p_atm_qi > 0) then + p3d => stateIn % tr(:,:,:,config % species % p_atm_qi) + set_non_neg = .true. + end if + case ("QS") + if (config % species % p_atm_qs > 0) then + p3d => stateIn % tr(:,:,:,config % species % p_atm_qs) + set_non_neg = .true. + end if + case ("QG") + if (config % species % p_atm_qg > 0) then + p3d => stateIn % tr(:,:,:,config % species % p_atm_qg) + set_non_neg = .true. + end if + case ("UWINDA") + p3d => stateIn % uwind + case ("VWINDA") + p3d => stateIn % vwind + case ("ZF") + k = 0 + do l = lay0, lay1 + do r = row0, row1 + do c = col0, col1 + k = k + 1 + buffer(k) = onebg * stateIn % phii(c,r,l+1) + end do + end do + end do + set_non_neg = .true. + case ("ZH") + k = 0 + do l = lay0, lay1 + do r = row0, row1 + do c = col0, col1 + k = k + 1 + buffer(k) = onebg * stateIn % phil(c,r,l) + end do + end do + end do + set_non_neg = .true. + case ("TA") + p3d => stateIn % temp + case default + ! set to 0 + end select + + else if (trim(fname) == trim(MET_DOT_3D)) then + + select case (trim(vname)) + case ("UWINDC") + ! u-wind is on C grid, while imported wind component are on A grid + ! this needs to be fixed + ! set to 0 for now + case ("VWINDC") + ! set to 0 for now + end select + ELSE IF (TRIM(FNAME) .EQ. 'MODIS_FPAR') THEN IF (TRIM(VNAME) .EQ. 'MODIS_FPAR_T') THEN @@ -1133,7 +1501,7 @@ LOGICAL FUNCTION XTRACT3 ( FNAME, VNAME, & end do END IF - !!Replace INIT_GASC,AERO,NONR,TRAC to INIT_CONC_1 (Wei Li) + !!Replace INIT_GASC,AERO,NONR,TRAC to INIT_CONC_1 ELSE IF ( TRIM(FNAME) .EQ. TRIM(INIT_CONC_1) ) THEN ! -- initialize gas-phase species (ppmV) @@ -1171,6 +1539,43 @@ LOGICAL FUNCTION XTRACT3 ( FNAME, VNAME, & END IF + !!interpx + if (associated(p2d)) then + k = 0 + do r = row0, row1 + do c = col0, col1 + k = k + 1 + buffer(k) = p2d(c,r) + end do + end do + else if (associated(p3d)) then + k = 0 + do l = lay0, lay1 + do r = row0, row1 + do c = col0, col1 + k = k + 1 + buffer(k) = p3d(c,r,l) + end do + end do + end do + end if + + if (set_non_neg) buffer(1:lbuf) = max( 0., buffer(1:lbuf) ) + + XTRACT3 = .TRUE. + + call aqm_model_get(config=config, rc=localrc) + if (aqm_rc_check(localrc, msg="Failure to retrieve model configuration", & + file=__FILE__, line=__LINE__)) return + + if (config % verbose) then + varname = vname + write(msgString, '(a,": interpx: ",a16,": ",a16,": min/max = ",2g20.8)') & + trim(config % name), fname, varname, & + minval(buffer(1:lbuf)), maxval(buffer(1:lbuf)) + call m3mesg(msgString) + end if + END FUNCTION XTRACT3 @@ -1231,24 +1636,29 @@ LOGICAL FUNCTION WRITE3_REAL2D( FNAME, VNAME, JDATE, JTIME, BUFFER ) type(aqm_state_type), pointer :: stateOut WRITE3_REAL2D = .TRUE. -!CTM_AOD_1 seems to be removed. (Wei Li) -! IF ( TRIM( FNAME ) .EQ. TRIM( CTM_AOD_1 ) ) THEN -! +!move to WRITE3_REAL4D below since we specify all model layers in CMAQ_Control_Misc.nml. +! IF ( TRIM( FNAME ) .EQ. TRIM( CTM_ELMO_1 ) ) THEN +! IF ( TRIM( FNAME ) .EQ. TRIM( CTM_DEPV_DIAG ) ) THEN !test depv +! IF ( TRIM( FNAME ) .EQ. TRIM( CTM_DRY_DEP_1 ) ) THEN !test depv +! IF ( TRIM( FNAME ) .EQ. TRIM( CTM_DUST_EMIS_1 ) ) THEN !test emis dust ! WRITE3_REAL2D = .FALSE. -! + ! IF ( TRIM( VNAME ) .EQ. TRIM( ALLVAR3 ) ) THEN -! +! IF ( TRIM( VNAME ) .EQ. 'VMASSJ' ) THEN +! IF ( TRIM( VNAME ) .EQ. 'ASOIL' ) THEN !DDEP +! IF ( TRIM( VNAME ) .EQ. 'PMCOARSE_SOIL' ) THEN !emission + ! 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 % aod = BUFFER -! + ! END IF -! + ! WRITE3_REAL2D = .TRUE. -! + ! END IF END FUNCTION WRITE3_REAL2D @@ -1274,32 +1684,34 @@ LOGICAL FUNCTION WRITE3_REAL4D( FNAME, VNAME, JDATE, JTIME, BUFFER ) type(aqm_state_type), pointer :: stateOut type(aqm_config_type), pointer :: config - integer, parameter :: p_pm25at = 23 + integer, parameter :: p_pm25at = 1 !(change from 23) WRITE3_REAL4D = .TRUE. -!CTM_PMDIAG_1 seems to be removed. (Wei Li) -! IF ( TRIM( FNAME ) .EQ. TRIM( CTM_PMDIAG_1 ) ) THEN -! -! WRITE3_REAL4D = .FALSE. -! -! IF ( TRIM( VNAME ) .EQ. TRIM( ALLVAR3 ) ) THEN -! -! nullify(config) -! nullify(stateOut) -! call aqm_model_get(config=config, stateOut=stateOut, rc=localrc) -! if (aqm_rc_check(localrc, msg="Failure to retrieve model output state", & -! file=__FILE__, line=__LINE__)) return -! -! do s = 0, config % species % ndiag - 2 -! stateOut % tr(:,:,:,config % species % p_diag_beg + s) = & -! buffer(:,:,:,p_pm25at + s) -! end do -! -! END IF -! -! WRITE3_REAL4D = .TRUE. -! -! END IF +!CTM_PMDIAG_1 seems to be removed. Use CTM_ELMO_1. + IF ( TRIM( FNAME ) .EQ. TRIM( CTM_ELMO_1 ) ) THEN + + WRITE3_REAL4D = .FALSE. + + IF ( TRIM( VNAME ) .EQ. TRIM( ALLVAR3 ) ) THEN + + nullify(config) + nullify(stateOut) + call aqm_model_get(config=config, stateOut=stateOut, rc=localrc) + if (aqm_rc_check(localrc, msg="Failure to retrieve model output state", & + file=__FILE__, line=__LINE__)) return + + do s = 0, config % species % ndiag - 2 + stateOut % tr(:,:,:,config % species % p_diag_beg + s) = & + buffer(:,:,:,p_pm25at + s) + end do + ! add AOD here; point to the 4th species in ELMO_INST + !stateOut % aod = BUFFER(:,:,1,4) + + END IF + + WRITE3_REAL4D = .TRUE. + + END IF END FUNCTION WRITE3_REAL4D From 2306c94c6b465408c67384eabcadda355af70670 Mon Sep 17 00:00:00 2001 From: Wei Li Date: Fri, 20 Oct 2023 19:19:33 -0500 Subject: [PATCH 3/5] Add other modified files --- src/model/src/AERO_EMIS.F | 579 -------------------------------------- 1 file changed, 579 deletions(-) delete mode 100644 src/model/src/AERO_EMIS.F diff --git a/src/model/src/AERO_EMIS.F b/src/model/src/AERO_EMIS.F deleted file mode 100644 index 0088145..0000000 --- a/src/model/src/AERO_EMIS.F +++ /dev/null @@ -1,579 +0,0 @@ - -!------------------------------------------------------------------------! -! The Community Multiscale Air Quality (CMAQ) system software is in ! -! continuous development by various groups and is based on information ! -! from these groups: Federal Government employees, contractors working ! -! within a United States Government contract, and non-Federal sources ! -! including research institutions. These groups give the Government ! -! permission to use, prepare derivative works of, and distribute copies ! -! of their work in the CMAQ system to the public and to permit others ! -! to do so. The United States Environmental Protection Agency ! -! therefore grants similar permission to use the CMAQ system software, ! -! but users are requested to provide copies of derivative works or ! -! products designed to operate in the CMAQ system to the United States ! -! Government without restrictions as to use by others. Software ! -! that is used with the CMAQ system but distributed under the GNU ! -! General Public License or the GNU Lesser General Public License is ! -! subject to their copyright restrictions. ! -!------------------------------------------------------------------------! - -C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: - MODULE AERO_EMIS - -C Emissions data and code required for the modal aerosol module in CMAQ -C Based on original codes by Dr. Francis S. Binkowski and J. Young - -C Dependent Upon: NONE - -C Revision History: - -C 30 Aug 01 J.Young: dyn alloc - Use HGRD_DEFN -C 09 Oct 03 J.Gipson: added MW array for AE emis species to module contents -C 31 Jan 05 J.Young: dyn alloc - establish both horizontal & vertical -C domain specifications in one module, GRID_CONF -C 26 Apr 05 P.Bhave: removed code supporting the "old type" of emission -C files that had unspeciated PM10 and PM2.5 only -C removed need for 'AERO_SPC.EXT' by declaring the -C required variables locally -C 13 Jun 05 P.Bhave: added vars needed for sea-salt emission processing -C inherit N_AE_EMIS,AE_EMIS,AE_EMIS_MAP from AE_EMIS.EXT -C moved RHO* parameters from RDEMIS_AE to this module -C for use by SSEMIS routine -C 24 Aug 07 J.Young: Modified to enable in-line plume rise calculation for -C 3D pt source emissions. Distinguish between PM (primary, -C unspeciated, file data) and AE (model speciated). Re- -C named RDEMIS_AE to GET_AERO_EMIS. -C 11 Apr 08 J.Kelly: added code to emit coarse surface area -C 4 Jan 10 J.Young: restructure; eliminate ref to older AERO versions -C 21 Feb 10 J.Young: move sea salt emissions to its own module (SSEMIS) -C 23 Apr 10 J.Young: replace include files with mechanism namelists -C 30 Apr 10 J.Young: update to use aero_reeng by Steve Howard, Prakash Bhave, -C Jeff Young, and Sergey Napelenok -C 23 Jul 10 D.Wong: remove CLOSE3 and BARRIER -C 24 Feb 11 J.Young: Reorganized module with initialization and timestepping -C procedures -C 25 Feb 11 J.Young: add windblown dust module -C 25 Mar 11 S.Roselle: replaced I/O API include files with UTILIO_DEFN -C 11 May 11 D.Wong: incorporated twoway model implementation -C 18 Aug 11 David Wong: In the merge inline point source PM species calculation, -C arrays EMBUFF and PMEMIS_PT have incorrect index values -C 17 Apr 13 J.Young: replace "SPFC ASO4" (found by Havala Pye) with "SPFC_ASO4" -C 07 Nov 14 J.Bash: Updated for the ASX_DATA_MOD shared data module. -C----------------------------------------------------------------------- - - USE AERO_DATA, ONLY: DESID_N_AERO_REF, N_MODE - USE DESID_VARS, ONLY: DESID_LAYS, DESID_STREAM_AERO, DESID_N_SRM, CELLVOL - - IMPLICIT NONE - SAVE -C aerosol emissions: [ppmv/s] for mass & number spcs, [m2/mol/s] for surface area spcs - PUBLIC DESID_SIZE_DIST, AERO_EMIS_INIT, DESID_INIT_SIZE_DIST, - & MAP_ISTRtoAERO, MAP_ISTRtoMODE, MAP_NUMtoISTR, MAP_SRFtoISTR, - & MAP_ISTRtoNUM, MAP_ISTRtoSRF, MAP_ISTRtoSD, DESID_STREAM_AERO, - & SD_SPLIT - PRIVATE - -C Variables for converting mass emissions rate to number emissions rate - REAL :: FACNUM( DESID_N_AERO_REF,N_MODE ) - -C Variables for converting mass emissions rate to 2nd moment emissions rate - REAL :: FACSRF( DESID_N_AERO_REF,N_MODE ) - -C Variables for Saving split factors between emission modes - REAL, ALLOCATABLE :: SD_SPLIT( :,: ) - -C Emission rate of all aerosol species interpolated to current time - INTEGER, ALLOCATABLE :: MAP_ISTRtoAERO( : ) - INTEGER, ALLOCATABLE :: MAP_ISTRtoMODE( : ) - INTEGER, ALLOCATABLE :: MAP_NUMtoISTR ( : ) - INTEGER, ALLOCATABLE :: MAP_SRFtoISTR ( : ) - INTEGER, ALLOCATABLE :: MAP_ISTRtoNUM ( : ) - INTEGER, ALLOCATABLE :: MAP_ISTRtoSRF ( : ) - INTEGER, ALLOCATABLE :: MAP_ISTRtoSD ( :,: ) - INTEGER, ALLOCATABLE :: MAP_AEROtoDIFF( :,: ) ! indices of aero species to CGRID - -C Miscellaneous variables - CHARACTER( 200 ) :: XMSG = ' ' - - CONTAINS - -C----------------------------------------------------------------------- - FUNCTION AERO_EMIS_INIT ( JDATE, JTIME, TSTEP ) RESULT ( SUCCESS) - -C Revision History: - -C 30 Aug 01 J.Young: dynamic allocation - Use INTERPX -C 29 Jul 03 P.Bhave: added compatibility with emission files that contain -C PM10, PEC, POA, PNO3, PSO4, and PMF, but do not -C contain PMC -C 20 Aug 03 J.Young: return aero emissions in molar mixing ratio, ppm units -C 09 Oct 03 J.Gipson: added MW array for AE emis species to module contents -C 01 Sep 04 P.Bhave: changed MW for primary organics from 120 to 220 g/mol, -C to match MWPOA in subroutine ORGAER3. -C 31 Jan 05 J.Young: dyn alloc - removed HGRD_ID, VGRID_ID, and COORD_ID -C include files because those parameters are now -C inherited from the GRID_CONF module -C 26 Apr 05 P.Bhave: removed code supporting the "old type" of emission -C files that had unspeciated PM10 and PM2.5 only -C removed need for 'AERO_CONST.EXT' by declaring the -C required variables locally -C simplified the CONVM, CONVN, CONVS calculations -C updated and enhanced in-line documentation -C 03 May 05 P.Bhave: fixed bug in the H2SO4 unit conversion, initially -C identified by Jinyou Liang of CARB -C 13 Jun 05 P.Bhave: calculate sea-salt emissions; execute if MECHNAME = AE4 -C read input fields from new OCEAN_1 file -C read extra input fields from MET_CRO_2D and MET_CRO_3D -C write diagnostic sea-salt emission file -C added TSTEP to call vector for diagnostic output file -C inherit MWs from AE_SPC.EXT instead of hardcoding -C find pointers to CGRID indices instead of hardcoding -C 08 Mar 07 P.Bhave& added capability for emission files that contain -C S.Roselle: POC or POA -C 30 Jan 08 P.Bhave: added compatibility with AE5 mechanisms -C 23 Mar 08 J.Young: modifications to allow for in-line point source emissions -C 11 Apr 08 J.Kelly: added code to emit coarse surface area -C 09 Sep 08 P.Bhave: backward compatibility with AE4 mechanisms -C 20 Feb 10 J.Young: move ssemis out to its own F90 module -C 24 Feb 11 J.Young: add windblown dust emissions option -C 25 Mar 11 S.Roselle: Replaced I/O API include files with UTILIO_DEFN -C 07 Jul 14 B.Hutzell: replaced mechanism include file(s) with fortran module -C 17 Sep 14 K.Fahey: Changed geometric mean diameter and geometric -C standard deviation of emitted particles according to -C Elleman and Covert (2010) -C 15 Apr 16 J.Young: Use aerosol factors from the AERO_DATA module's named constants; -C Moved K.Fahey's mods to geometric mean diameter and standard -C deviation to the AERO_DATA module - -C References: -C CRC76, "CRC Handbook of Chemistry and Physics (76th Ed)", -C CRC Press, 1995 -C Elleman & Covert, "Aerosol size distribution modeling with the Community -C Multiscale Air Quality modeling system in the Pacific -C Northwest: 3. Size distribution of particles emitted -C into a mesoscale model", J. Geophys. Res., Vol 115, -C No D3, doi:10.1029/2009JD012401, 2010 -C Hobbs, P.V. "Basic Physical Chemistry for the Atmospheric Sciences", -C Cambridge Univ. Press, 206 pp, 1995. -C Snyder, J.P. "Map Projections-A Working Manual", U.S. Geological Survey -C Paper 1395 U.S.GPO, Washington, DC, 1987. -C Binkowski & Roselle Models-3 Community Multiscale Air Quality (CMAQ) -C model aerosol component 1: Model Description. -C J. Geophys. Res., Vol 108, No D6, 4183 -C doi:10.1029/2001JD001409, 2003 -C----------------------------------------------------------------------- - - USE AERO_DATA, ONLY: DESID_AERO_REF, N_AEROSPC, AEROSPC, - & AERO_MISSING, MAP_AERO - USE GRID_CONF, ONLY: GDTYP_GD, XCELL_GD, YCELL_GD, YORIG_GD, GL_NROWS, X3FACE_GD - USE DUST_EMIS, ONLY: DUST_EMIS_INIT - USE DESID_VARS, ONLY: MAP_ISTRtoEMVAR - USE PRECURSOR_DATA, ONLY: MAP_PRECURSOR - USE RUNTIME_VARS, ONLY: OCEAN_CHEM, WB_DUST - USE SSEMIS, ONLY: SSEMIS_INIT - USE UTILIO_DEFN !(Wei Li), ONLY: INDEX1, M3EXIT, LATGRD3, XSTAT1, XSTAT2 - USE VDIFF_MAP, ONLY : N_SPC_DIFF, DIFF_SPC - - INCLUDE SUBST_CONST ! physical and mathematical constants - INCLUDE SUBST_FILES_ID ! file name parameters - -C Arguments: - - INTEGER, INTENT( IN ) :: JDATE ! current model date, coded YYYYDDD - INTEGER, INTENT( IN ) :: JTIME ! current model time, coded HHMMSS - INTEGER, INTENT( IN ) :: TSTEP ! time step vector (HHMMSS) - ! TSTEP(1) = local output step - LOGICAL SUCCESS - -C External Functions: - INTEGER, EXTERNAL :: FINDEX ! looks up number in table. - -C Local Variables: - REAL DGV, SG, SPLIT_ACCUM - -C Domain decomposition info from emission and meteorology files - INTEGER GXOFF, GYOFF ! origin offset - -C Miscellaneous variables - INTEGER STATUS ! ENV..., ALLOCATE status - CHARACTER( 16 ), SAVE :: PNAME = 'AERO_EMIS_INIT ' - CHARACTER( 16 ) :: VNAME ! temp var for species names - CHARACTER( 50 ) :: VARDESC ! variable for reading environ. variables - INTEGER L, N, S, V, IAERO, ISRM, ! Loop indices - & IEM, IDIFF, ISPC - -C ---------------------------------------------------------------------- - - SUCCESS = .TRUE. - -C *** Map data modules - CALL MAP_AERO() - CALL MAP_PRECURSOR() - -C *** set up for sea-spray emission processing - IF ( OCEAN_CHEM ) THEN - IF ( .NOT. SSEMIS_INIT( JDATE, JTIME, TSTEP ) ) THEN - XMSG = 'Failure initializing sea-spray emission processing' - CALL M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT2 ) - END IF - END IF - -C *** set up for dust emission processing - IF ( WB_DUST ) THEN - IF ( .NOT. DUST_EMIS_INIT( JDATE, JTIME, TSTEP ) ) THEN - XMSG = 'Failure initializing dust emission processing' - CALL M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT2 ) - END IF - END IF - -C *** Set up emissions size distribution arrays - ! Calculate factors for converting 3rd moment emission rates into - ! number and surface area emission rates. See Equation 7b of - ! Binkowski & Roselle (2003) - DO IEM = 1,DESID_N_AERO_REF - DO N = 1, N_MODE - DGV = DESID_AERO_REF( IEM )%DGVEM( N ) - SG = DESID_AERO_REF( IEM )%SGEM ( N ) - - IF ( DESID_AERO_REF( IEM )%SPLIT( N ) .GT. 0.0 ) THEN - FACNUM( IEM,N ) = EXP( 4.5 * LOG( SG ) ** 2 ) / DGV ** 3 - FACSRF( IEM,N ) = PI * EXP( 0.5 * LOG( SG ) ** 2 ) / DGV - ELSE - FACNUM( IEM,N ) = 0.0 - FACSRF( IEM,N ) = 0.0 - END IF - END DO - - END DO - - ! Map the Modal-Dependent Names to Transported Species - ALLOCATE ( MAP_AEROtoDIFF( N_AEROSPC, N_MODE ) ) - DO ISPC = 1,N_AEROSPC - DO N = 1,N_MODE - MAP_AEROtoDIFF( ISPC, N ) = INDEX1( AEROSPC( ISPC )%name( N ), - & N_SPC_DIFF, DIFF_SPC ) - END DO - END DO - - - ! Modify the reference emissions splits based on what transported - ! aerosol species are actually available. For example, if the aerosol - ! namelist only includes the accumulation mode (J) but not the - ! Aitken mode (I) for a particular species, then the split for - ! Aitken mode should be added to the Accumulation mode. Save - ! these scale factors as a function of transported species and - ! mode. - ALLOCATE( SD_SPLIT( N_SPC_DIFF, DESID_N_AERO_REF ) ) - SD_SPLIT = 0.0 - DO IEM = 1,DESID_N_AERO_REF - ! For the Fine Mode Reference Distribution, lump Aitken - ! with Accumulation mode if Aitken Mode does not exist - IF ( DESID_AERO_REF( IEM )%NAME .EQ. 'FINE_REF' ) THEN - DO ISPC = 1,N_AEROSPC - SPLIT_ACCUM = 0.0 - DO N = 1,N_MODE-1 - IF ( AERO_MISSING( ISPC,N ) ) THEN - SPLIT_ACCUM = SPLIT_ACCUM + DESID_AERO_REF( IEM )%SPLIT( N ) - ELSE - SD_SPLIT( MAP_AEROtoDIFF( ISPC,N ),IEM ) = - & SD_SPLIT( MAP_AEROtoDIFF( ISPC,N ),IEM ) + - & DESID_AERO_REF( IEM )%SPLIT( N ) + SPLIT_ACCUM - SPLIT_ACCUM = 0.0 - END IF - END DO - END DO - ELSE - ! Arbitrary Distribution -> Apply factor to species - ! if it exists in each mode - DO ISPC = 1, N_AEROSPC - DO N = 1, N_MODE - IF ( .NOT. AERO_MISSING( ISPC,N ) ) THEN - SD_SPLIT( MAP_AEROtoDIFF( ISPC,N ),IEM ) = - & DESID_AERO_REF( IEM )%SPLIT( N ) - END IF - END DO - END DO - END IF - END DO - - ALLOCATE ( MAP_NUMtoISTR ( N_MODE ), - & MAP_SRFtoISTR ( N_MODE ), STAT = STATUS ) - CALL CHECKMEM( STATUS, 'MAP_NUMtoEM', PNAME ) - CALL CHECKMEM( STATUS, 'MAP_SRFtoEM', PNAME ) - - END FUNCTION AERO_EMIS_INIT - -C----------------------------------------------------------------------- - - SUBROUTINE DESID_INIT_SIZE_DIST ( JDATE, JTIME ) - -C EM_SD_INIT initializes the structures that map modes and streams to -C reference modes including splits, diameters, and standard deviations. - -C----------------------------------------------------------------------- - USE AERO_DATA, ONLY: DESID_AERO_REF, DESID_N_AERO_REF - USE DESID_VARS, ONLY: DESID_SD_NML - USE DESID_UTIL, ONLY: DESID_GET_RULE_STREAMS - USE UTILIO_DEFN, ONLY: INDEX1, XSTAT1 - - IMPLICIT NONE - - INTEGER, INTENT( IN ) :: JDATE ! current model date, coded YYYYDDD - INTEGER, INTENT( IN ) :: JTIME ! current model time, coded HHMMSS - INTEGER ISRM - - INTEGER :: N_SD_RULE - INTEGER :: N_SD( DESID_N_SRM ) - CHARACTER( 16 ) :: SD_NAME( DESID_N_SRM, 10 ) - INTEGER :: SD( DESID_N_SRM, 10 ) - LOGICAL :: RULE_STREAM( DESID_N_SRM ) - CHARACTER( 16 ) :: CSUR - CHARACTER( 16 ), SAVE :: PNAME = 'EM_SD_INIT ' - CHARACTER( 20 ) :: DESID_AERO_REF_CAPS( DESID_N_AERO_REF ) - - INTEGER IRULE, ISUR, N, NLEN, ISD, IM, IEM, NRULE - LOGICAL :: LREMOVE - - ! Find Total Number of Size Distribution Registries - N_SD_RULE = 0 - DO IRULE = 1,SIZE( DESID_SD_NML ) - IF ( DESID_SD_NML( IRULE )%STREAM .EQ. '' ) EXIT - N_SD_RULE = IRULE - END DO - - ! First Load all of the Streams with the Default FINE, COARSE, and - ! AERO Mode references - SD = 0 - SD_NAME = '' - - ! Capitalize EM_AERO_REF(:)%NAME - DO IM = 1,DESID_N_AERO_REF - DESID_AERO_REF_CAPS( IM ) = DESID_AERO_REF( IM )%NAME - CALL UPCASE( DESID_AERO_REF_CAPS( IM ) ) - ENDDO - - DO ISRM = 1,DESID_N_SRM - N_SD( ISRM ) = 2 - SD_NAME( ISRM,1 ) = 'FINE' - SD( ISRM,1 ) = INDEX1( 'FINE_REF', DESID_N_AERO_REF, DESID_AERO_REF_CAPS( : ) ) - SD_NAME( ISRM,2 ) = 'COARSE' - SD( ISRM,2 ) = INDEX1( 'COARSE_REF', DESID_N_AERO_REF, DESID_AERO_REF_CAPS( : ) ) - END DO - - ! Now Modify those defaults or add new modes to desired streams - DO IRULE = 1, N_SD_RULE - ! Expand Size Distribution Rule to All Streams if Requested - LREMOVE = .FALSE. - IF ( DESID_SD_NML( IRULE )%STREAM .EQ. '' ) CYCLE - CALL DESID_GET_RULE_STREAMS( DESID_SD_NML( IRULE )%STREAM, IRULE, - & RULE_STREAM, LREMOVE ) - IF ( LREMOVE ) CYCLE - - ! Loop through streams, set defaults, and build map array - DO ISRM = 1, DESID_N_SRM - IF ( RULE_STREAM( ISRM ) ) THEN - ! This Stream is Being Modified by a Size Distribution - ! rule - CALL UPCASE( DESID_SD_NML( IRULE )%MODE_REF ) - IF ( DESID_SD_NML( IRULE )%MODE .EQ. 'FINE' ) THEN - ! Overwrite the FINE mode. All fine particle species - ! will go to this mode by default - SD( ISRM,1 ) = INDEX1( DESID_SD_NML( IRULE )%MODE_REF, - & DESID_N_AERO_REF, DESID_AERO_REF_CAPS( : ) ) - IF ( SD( ISRM,1 ) .EQ. 0 ) THEN - WRITE( XMSG,'(A,A,A,/,A,I2,A)' ), '*** Reference Aerosol Mode (', - & DESID_SD_NML( IRULE )%MODE_REF, 'Specified in Emissions Size ', - & 'Dist Rule ',IRULE,' does not exist in AERO_DATA.' - CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 ) - END IF - - ELSEIF ( DESID_SD_NML( IRULE )%MODE .EQ. 'COARSE' ) THEN - ! Overwrite the COARSE mode. All coarse particle - ! species will go to this mode by default - SD( ISRM,2 ) = INDEX1( DESID_SD_NML( IRULE )%MODE_REF, - & DESID_N_AERO_REF, DESID_AERO_REF_CAPS( : ) ) - IF ( SD( ISRM,2 ) .EQ. 0 ) THEN - WRITE( XMSG,'(A,A,A,/,A,I2,A)' ), '*** Reference Aerosol Mode (', - & DESID_SD_NML( IRULE )%MODE_REF, 'Specified in Emissions Size ', - & 'Dist Rule ',IRULE,' does not exist in AERO_DATA.' - CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 ) - END IF - - ELSE - ! Add a New Available Mode. For example, add a mode - ! just for BC, call it PUREBC, and make sure the AEC - ! for this stream is pointing to this mode. Also make - ! sure you set AEC for FINE mode aerosol to 0.0 if - ! you have default mapping turned on. - N_SD( ISRM ) = N_SD( ISRM ) + 1 - SD_NAME( ISRM,N_SD( ISRM ) ) = DESID_SD_NML( IRULE )%MODE - SD( ISRM,N_SD( ISRM ) ) = INDEX1( DESID_SD_NML( IRULE )%MODE_REF, - & DESID_N_AERO_REF, DESID_AERO_REF_CAPS( : ) ) - IF ( SD( ISRM,N_SD( ISRM )) .EQ. 0 ) THEN - WRITE( XMSG,'(A,A,A,/,A,I2,A)' ), '*** Reference Aerosol Mode (', - & DESID_SD_NML( IRULE )%MODE_REF, 'Specified in Emissions Size ', - & 'Dist Rule ',IRULE,' does not exist in AERO_DATA.' - CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 ) - END IF - - END IF - END IF - END DO - END DO - - ! Finally, transfer this data to a global variable which - ! captures and organizes the modes of each stream - ALLOCATE( DESID_STREAM_AERO( DESID_N_SRM ) ) - DO ISRM = 1,DESID_N_SRM - N = N_SD( ISRM ) - DESID_STREAM_AERO( ISRM )%LEN = N + 1 - ALLOCATE( DESID_STREAM_AERO( ISRM )%NAME( N+1 ) ) - ALLOCATE( DESID_STREAM_AERO( ISRM )%REF( N+1 ) ) - ALLOCATE( DESID_STREAM_AERO( ISRM )%FACNUM( N+1,N_MODE ) ) - ALLOCATE( DESID_STREAM_AERO( ISRM )%FACSRF( N+1,N_MODE ) ) - - DESID_STREAM_AERO( ISRM )%NAME( 2:N+1 ) = SD_NAME( ISRM,1:N ) - DESID_STREAM_AERO( ISRM )%REF( 2:N+1 ) = SD( ISRM,1:N ) - DESID_STREAM_AERO( ISRM )%NAME( 1 ) = 'GAS' - DESID_STREAM_AERO( ISRM )%REF( 1 ) = 0 - - ! Map Factors for Converting Aerosol Mass to Number and - ! Surface Area to each Emission Stream - DESID_STREAM_AERO( ISRM )%FACNUM( :,: ) = 0.0 - DESID_STREAM_AERO( ISRM )%FACSRF( :,: ) = 0.0 - DO ISD = 2,N+1 - IEM = DESID_STREAM_AERO( ISRM )%REF( ISD ) - DO IM = 1,N_MODE - DESID_STREAM_AERO( ISRM )%FACNUM( ISD,IM ) = FACNUM( IEM,IM ) - DESID_STREAM_AERO( ISRM )%FACSRF( ISD,IM ) = FACSRF( IEM,IM ) - END DO - END DO - END DO - - END SUBROUTINE DESID_INIT_SIZE_DIST - - -C----------------------------------------------------------------------- - - SUBROUTINE DESID_SIZE_DIST ( ISRM, VDEMIS, NL ) - -C EMISS_SIZE_DIST distributes bulk aerosol emissions into size space -C using parameters precompiled in the AERO_DATA module. -C -C Revision History: - -C 16 AUG 17 BMURPHY: Created -C -C ---------------------------------------------------------------------- - - USE AERO_DATA, ONLY: AEROSPC, N_AEROSPC, AEROSPC_MWINV - USE AEROMET_DATA, ONLY: F6DPI - USE ASX_DATA_MOD, ONLY: MET_DATA - USE DESID_VARS, ONLY: DESID_N_ISTR, IDUSTSRM, ISEASRM - USE GRID_CONF, ONLY: NCOLS, NROWS - USE SSEMIS, ONLY: SEA_FACTNUM, SEA_FACTSRF - - INTEGER, INTENT( IN ) :: ISRM, NL - REAL, INTENT( INOUT ) :: VDEMIS ( :,:,:,: ) - - INTEGER :: N, S, IAERO, IM, ISD, ISTR ! Looping Variables - INTEGER :: ROW, COL, LAY, N_SD, INUM, ISRF - REAL :: FACNUM, FACSRF, MW_FAC - REAL, ALLOCATABLE :: EMISM3( :,:,:,:,: ) - REAL, ALLOCATABLE, SAVE :: GSFAC( :,:,: ) - REAL, ALLOCATABLE, SAVE :: DENS_FAC( : ) - REAL, PARAMETER :: F6DPIM9 = 1.0E-9 * F6DPI ! 1.0E-9 = Kg/ug - LOGICAL, SAVE :: FIRST_TIME = .TRUE. - -C *** Initialize Variables - IF ( FIRST_TIME ) THEN - FIRST_TIME = .FALSE. - ALLOCATE( GSFAC ( DESID_LAYS,NCOLS,NROWS ) ) - - ALLOCATE( DENS_FAC( N_AEROSPC ) ) - DO IAERO = 1,N_AEROSPC - DENS_FAC( IAERO ) = F6DPIM9 / AEROSPC( IAERO )%DENSITY - END DO - - END IF - - N_SD = DESID_STREAM_AERO( ISRM )%LEN - ALLOCATE( EMISM3( DESID_LAYS,NCOLS,NROWS,N_MODE,N_SD ) ) - EMISM3 = 0.0 - -C *** Calculate scaling factor for converting mass emissions into [ug/m3/s] -C note: RJACM converts grid heights from sigma coordinates to meters -C Also calculate scaling factors for converting to molar-mixing-ratio units - DO LAY = 1,NL - GSFAC( LAY,:,: ) = Met_Data%RJACM( :,:,LAY ) / CELLVOL( :,:,LAY ) ![ug/s] to [ug/m3/s] - END DO - -C *** Apply Aerosol Size Distribution - DO ISTR = 1, DESID_N_ISTR - ! Find which Size Distribution or Phase this emissions species belongs - ! to for this stream. If the value is a 0, then there are no emissions - ! for this species from this stream. If it is a 1, then this species is - ! a gas and the following aerosol conversions should be skipped. - ISD = MAP_ISTRtoSD( ISTR,ISRM ) - IF ( ISD .LE. 1 ) CYCLE - - ! Look up Aerosol Species and Mode of Interest - IAERO = MAP_ISTRtoAERO( ISTR ) !This maps to the CMAQ aerosol - ! species so we can retrieve density - IM = MAP_ISTRtoMODE( ISTR ) !This maps to the internal CMAQ modes - ! (ie. I, J, and K) - !DENS_FAC = F6DPIM9 / AEROSPC( IAERO )%DENSITY - - ! Convert Aerosol from [g/s] to [ug/m3/s] for all streams - ! except Dust and Sea Spray. For those streams, convert - ! [g/m3/s] to [ug/m3/s] - VDEMIS( ISTR,1:NL,:,: ) = VDEMIS( ISTR,1:NL,:,: ) * 1.0E6 - IF ( ISRM .NE. ISEASRM .AND. ISRM .NE. IDUSTSRM ) THEN - VDEMIS( ISTR,1:NL,:,: ) = VDEMIS( ISTR,1:NL,:,: ) * GSFAC( 1:NL,:,: ) - END IF - - ! Sum Total Volume of Mode N [m3/m3/s] - IF ( .NOT. AEROSPC( IAERO )%TRACER ) - & EMISM3( 1:NL,:,:,IM,ISD ) = EMISM3( 1:NL,:,:,IM,ISD ) + - & VDEMIS( ISTR,1:NL,:,: ) * DENS_FAC( IAERO ) - - ! Convert Mass Emission Rates from [ug/m3/s] to [umol/m3/s] - VDEMIS( ISTR,1:NL,:,: ) = VDEMIS( ISTR,1:NL,:,: ) * AEROSPC_MWINV( IAERO ) - END DO - -C *** Calculate the number emissions rate for each mode [1/m3/s], using -C Equation 7b of Binkowski & Roselle (2003). -C Calculate the surface area emissions rate for the fine modes [m2/m3/s], -C using Equation 7c of Binkowski & Roselle (2003). Multiplying by PI -C converts 2nd moment to surface area. - - DO ISD = 2, N_SD ! Skip the Index for the Gas Phase - IF ( ISRM .EQ. ISEASRM ) THEN - ! Apply Spatially-Dependent Number and Surface Area Scale Factors - DO IM = 1, N_MODE - INUM = MAP_NUMtoISTR(IM) - VDEMIS( INUM,1,:,: ) = VDEMIS( INUM,1,:,: ) - & + EMISM3( 1,:,:,IM,ISD ) * SEA_FACTNUM( IM,:,: ) - - ISRF = MAP_SRFtoISTR(IM) - VDEMIS( ISRF,1,:,: ) = VDEMIS( ISRF,1,:,: ) - & + EMISM3( 1,:,:,IM,ISD ) * SEA_FACTSRF( IM,:,: ) - END DO - ELSE - ! Apply Homogeneous Scale Factors Consistent with this Stream - DO IM = 1, N_MODE - INUM = MAP_NUMtoISTR(IM) - FACNUM = DESID_STREAM_AERO( ISRM )%FACNUM( ISD,IM ) - VDEMIS( INUM,1:NL,:,: ) = VDEMIS( INUM,1:NL,:,: ) + EMISM3( 1:NL,:,:,IM,ISD ) * FACNUM - - ISRF = MAP_SRFtoISTR(IM) - FACSRF = DESID_STREAM_AERO( ISRM )%FACSRF( ISD,IM ) - VDEMIS( ISRF,1:NL,:,: ) = VDEMIS( ISRF,1:NL,:,: ) + EMISM3( 1:NL,:,:,IM,ISD ) * FACSRF - END DO - END IF - END DO - - END SUBROUTINE DESID_SIZE_DIST - - END MODULE AERO_EMIS - From 6351dd4ed226464290cd2fd4d2c1eb50454bb68e Mon Sep 17 00:00:00 2001 From: lwcugb <35088762+lwcugb@users.noreply.github.com> Date: Fri, 20 Oct 2023 20:51:34 -0400 Subject: [PATCH 4/5] Update .gitmodules point to CMAQ 5.4+ --- .gitmodules | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.gitmodules b/.gitmodules index b13c486..6031757 100644 --- a/.gitmodules +++ b/.gitmodules @@ -1,4 +1,4 @@ [submodule "src/model/CMAQ"] path = src/model/CMAQ url = https://github.com/USEPA/CMAQ - branch = 5.2.1 + branch = 5.4+ From 5e768a662e332ca9225f4f7e5c83d104a40a678c Mon Sep 17 00:00:00 2001 From: lwcugb <35088762+lwcugb@users.noreply.github.com> Date: Sat, 21 Oct 2023 01:47:45 -0400 Subject: [PATCH 5/5] Update aqm_methods.F90 Forget to turn on AOD --- src/shr/aqm_methods.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/shr/aqm_methods.F90 b/src/shr/aqm_methods.F90 index 49c1278..0098d07 100644 --- a/src/shr/aqm_methods.F90 +++ b/src/shr/aqm_methods.F90 @@ -1705,7 +1705,7 @@ LOGICAL FUNCTION WRITE3_REAL4D( FNAME, VNAME, JDATE, JTIME, BUFFER ) buffer(:,:,:,p_pm25at + s) end do ! add AOD here; point to the 4th species in ELMO_INST - !stateOut % aod = BUFFER(:,:,1,4) + stateOut % aod = BUFFER(:,:,1,4) END IF