diff --git a/CCTM/src/ddm3d/DDM3D_CHEM.F b/CCTM/src/ddm3d/DDM3D_CHEM.F index aa3fddda3..c598e524d 100644 --- a/CCTM/src/ddm3d/DDM3D_CHEM.F +++ b/CCTM/src/ddm3d/DDM3D_CHEM.F @@ -13,6 +13,7 @@ Module DDM3D_CHEM Implicit None + Real(8), Allocatable, Save :: YCHECK( : ) ! Concs used by DDM Real(8), Allocatable, Save :: YCDDM( : ) ! Concs used by DDM ! (avg of pre- and post-chem, ! or mid-chemstep concs, @@ -31,16 +32,20 @@ Module DDM3D_CHEM Real, Allocatable, Save :: PRD( : ) Real, Allocatable, Save :: PRD_RATE( : ) ! Contribution to PRD from rxn rate sens Real, Allocatable, Save :: SOLD( : ) + Integer, Allocatable, Save :: IPVT ( : ) ! an integer vector of pivot indices. + Integer, Allocatable :: SENS_INDEX( : ) ! CGRID_INDEX translated by JNEW2OLD + Integer, Allocatable :: MECH_INDEX( : ) ! CGRID_INDEX translated by JNEW2OLD C Variables used for hddm-3d - Real*8, Allocatable, Save :: SRK2 ( : ) ! rate constants + Real(8), Allocatable, Save :: SRK2 ( : ) ! rate constants LOGICAL, Allocatable, Save :: ORDER1 ( : ) ! true if order 1; else, false - Real, Allocatable, Save :: PDT2( :, : ) ! Used for 2nd order call of JAC - Real, Allocatable, Save :: SMID( :, : ) ! SENGRID in middle of timestep - Real(8), Allocatable, Save :: SMIDJAC( : ) ! SMID for the relevant 1st order - ! sensitivity parameter - Real*8, Allocatable, Save :: RK ( : ) + Real, Allocatable, Save :: PDT2( :, : ) ! Used for 2nd order call of JAC + Real, Allocatable, Save :: SMID( :, : ) ! SENGRID in middle of timestep + Real(8), Allocatable, Save :: SEND( :, : ) ! SENGRID at end of timestep + Real(8), Allocatable, Save :: SMIDJAC( : ) ! SMID for the relevant 1st order + ! sensitivity parameter + Real(8), Allocatable, Save :: RK ( : ) Integer N_EBI_MID ! the midpoint ebi step; half of N_EBI_STEPS Logical ODD_STEPS ! true if N_EBI_STEPS is odd @@ -51,27 +56,35 @@ Module DDM3D_CHEM C----------------------------------------------------------------------- Subroutine INIT_DDM3D_CHEM - Use RXNS_DATA, Only: NUMB_MECH_SPC, NRXNS + Use RXNS_DATA, Only: NUMB_MECH_SPC, NRXNS, CGRID_INDEX Use UTILIO_DEFN ! IOAPI parameters and functions declarations Implicit None Character( 16 ), Save :: PNAME = 'INIT_DDM3D_CHEM' - Integer LOGDEV - Character( 96 ) :: XMSG = ' ' - Integer ALLOCSTAT + Integer :: LOGDEV + Character( 96 ) :: XMSG = ' ' + Integer :: ALLOCSTAT + Integer :: ISPC + LOGDEV = INIT3 () + ALLOCATE( YCHECK ( NUMB_MECH_SPC), STAT = ALLOCSTAT ) + IF ( ALLOCSTAT .NE. 0 ) THEN + XMSG = 'Failure allocating YCHECK' + CALL M3EXIT ( PNAME, 0, 0, XMSG, XSTAT2 ) + END IF + ALLOCATE( YCDDM ( NUMB_MECH_SPC), STAT = ALLOCSTAT ) IF ( ALLOCSTAT .NE. 0 ) THEN - XMSG = 'Failure allocating DDM_CHECK' + XMSG = 'Failure allocating YCDDM' CALL M3EXIT ( PNAME, 0, 0, XMSG, XSTAT2 ) END IF ALLOCATE( DDM_CHECK( NUMB_MECH_SPC), STAT = ALLOCSTAT ) IF ( ALLOCSTAT .NE. 0 ) THEN - XMSG = 'Failure allocating YCDDM' + XMSG = 'Failure allocating DDM_CHECK' CALL M3EXIT ( PNAME, 0, 0, XMSG, XSTAT2 ) END IF @@ -99,12 +112,19 @@ Subroutine INIT_DDM3D_CHEM CALL M3EXIT ( PNAME, 0, 0, XMSG, XSTAT2 ) END IF - ALLOCATE( SMID( NPMAX, NUMB_MECH_SPC ), STAT = ALLOCSTAT ) + ALLOCATE( SMID( NUMB_MECH_SPC,NPMAX ), STAT = ALLOCSTAT ) IF ( ALLOCSTAT .NE. 0 ) THEN XMSG = 'Failure allocating SMID' CALL M3EXIT ( PNAME, 0, 0, XMSG, XSTAT2 ) END IF + ALLOCATE( SEND( 1,NUMB_MECH_SPC ), STAT = ALLOCSTAT ) + IF ( ALLOCSTAT .NE. 0 ) THEN + XMSG = 'Failure allocating SEND' + CALL M3EXIT ( PNAME, 0, 0, XMSG, XSTAT2 ) + END IF + SEND = 0.0D0 + ALLOCATE( PRD( NUMB_MECH_SPC ), & SOLD( NUMB_MECH_SPC ), & IPVT ( NUMB_MECH_SPC ), @@ -115,16 +135,29 @@ Subroutine INIT_DDM3D_CHEM CALL M3EXIT ( PNAME, 0, 0, XMSG, XSTAT2 ) END IF + ALLOCATE( SENS_INDEX( NUMB_MECH_SPC ), + & MECH_INDEX( NUMB_MECH_SPC ), STAT = ALLOCSTAT ) + IF ( ALLOCSTAT .NE. 0 ) THEN + XMSG = 'Failure allocating SENS_INDEX or MECH_INDEX' + CALL M3EXIT ( PNAME, 0, 0, XMSG, XSTAT2 ) + END IF + + DO ISPC = 1, NUMB_MECH_SPC + MECH_INDEX(ISPC) = CGRID_INDEX(ISPC) + SENS_INDEX(ISPC) = CGRID_INDEX(ISPC) + END DO + Return End Subroutine INIT_DDM3D_CHEM C----------------------------------------------------------------------- - Subroutine SOLVE_DDM3D_CHEM( C,R,L, CHEMSTEP) + Subroutine SOLVE_DDM3D_CHEM( C,R,L,CHEMSTEP,OLD2NEW,NEW2OLD ) Use DDM3D_DEFN, Only: SENGRID, NPMAX, NP, DATENUM, IPT, & IDATE, HIGH, IREGION, IRXN, IPARM, STARTDATE - Use RXNS_DATA, Only: NRXNS, NREACT, NPRDCT, SC, IRR, CGRID_INDEX, NUMB_MECH_SPC + Use RXNS_DATA, Only: NRXNS, NREACT, NPRDCT, SC, IRR, CGRID_INDEX, NUMB_MECH_SPC, + & RXLABEL, CHEMISTRY_SPC Use MECHANISM_FUNCTIONS @@ -134,18 +167,29 @@ Subroutine SOLVE_DDM3D_CHEM( C,R,L, CHEMSTEP) Character( 16 ), Save :: PNAME = 'SOLVE_DDM3D_CHEM' - Integer C,R,L - Real( 8 ) CHEMSTEP + Integer, Intent( In ) :: C,R,L + Real( 8 ), Intent( In ) :: CHEMSTEP + Integer, Optional, Intent( In ) :: OLD2NEW( :,: ) + Integer, Optional, Intent( In ) :: NEW2OLD( :,: ) Integer I,J,S,N + Integer JROW, JCOL Integer INFO ! see s_lu.F - Real SUM + + Real(8) SUMAT + Real(8) SUMSP + Real(8) TOTAL + Real(8) DIFF + Real(8) DIFFSP + Real IREGTEMP ! Holds relevant value of IREGION Real KSTEP ! Holds k times timestep(in min) Real RXNMULT ! Holds product of concs of reactants + Integer NRCT ! Counter over reactants Integer NPROD ! Counter over products - + Integer IREACT ! index for reaction reactant + Integer IPROD ! index for reaction produce Integer HIPARM( 2 ) ! index numbers of the 1st order sens ! parameters to which ! we're taking @@ -153,33 +197,124 @@ Subroutine SOLVE_DDM3D_CHEM( C,R,L, CHEMSTEP) Integer HITMP1 Integer HITMP2 - Character( 96 ) :: XMSG = ' ' Logical, Save :: FIRSTIME = .TRUE. + Logical, Save :: REORDER = .TRUE. ! reorder YCDDM Integer, Save :: LOGDEV - + Character( 96 ) :: XMSG = ' ' IF ( FIRSTIME ) THEN - FIRSTIME = .FALSE. - LOGDEV = INIT3() + FIRSTIME = .FALSE. + LOGDEV = INIT3() + DDM_LOG = LOGDEV + ERROR_LOG = LOGDEV CALL SET_MECHANISM( ) ! determine formulas for Mechanism Jacobain and Species Rate of Change + + + IF( PRESENT( OLD2NEW ) .AND. PRESENT( OLD2NEW ) )THEN + + REORDER = .FALSE. ! no reordering because YCDMM is already redeorderd by chemsolver + +! YCHECK = 0.0D0 +! DO I = 1,NUMB_MECH_SPC +! S = NEW2OLD( I,1 ) +! YCHECK( S ) = YCDDM( I ) +! SEND(1,I) = YCDDM( I ) +! WRITE(LOGDEV,'(A,A,ES16.4)')' ORIG: ' // CHEMISTRY_SPC(I),' = ',YCHECK(S) +! END DO +! DO I = 1,NUMB_MECH_SPC +! S = NEW2OLD( I,1 ) +! WRITE(LOGDEV,'(A,A,2(ES16.4,1X))')' SORT: ' // CHEMISTRY_SPC(S),' = ',YCDDM( I ),YCHECK(S)-YCDDM( I ) +! END DO +! CALL EVALUATE_F_JAC_MECH( YCHECK, SRK, PDT2 ) ! Evaluate Jacobian based on YCDDM and SKR values + +! overwrite reordering using chemsolver conversing maps + + JNEW2OLD(1:NUMB_MECH_SPC,1) = NEW2OLD(1:NUMB_MECH_SPC,1) + JOLD2NEW(1:NUMB_MECH_SPC,1) = OLD2NEW(1:NUMB_MECH_SPC,1) + RESET_JACOBIAN = .TRUE. + DO RXN = 1, NRXNS + DO NRCT = 1, NREACT( RXN ) + IREACT = ISPECIES_REACTION( NRCT,RXN ) + ISPECIES_REACTION( NRCT,RXN ) = OLD2NEW( IREACT,1 ) + END DO + + DO NPROD = 1, NPRDCT( RXN ) + IPROD = ISPECIES_REACTION( NPROD+3,RXN ) + ISPECIES_REACTION( NPROD+3,RXN ) = OLD2NEW( IPROD,1 ) + END DO + END DO + +! DO S = 1, NUMB_MECH_SPC +! MECH_INDEX(S) = CGRID_INDEX(JOLD2NEW(S,1)) +! SENS_INDEX(S) = CGRID_INDEX(JNEW2OLD(S,1)) +! END DO + +! ELSE +! +! DO S = 1, NUMB_MECH_SPC +! MECH_INDEX(S) = CGRID_INDEX(JOLD2NEW(S,1)) +! SENS_INDEX(S) = CGRID_INDEX(JNEW2OLD(S,1)) +! END DO +! DO I = 1,NUMB_MECH_SPC +! S = JOLD2NEW( I,1 ) +! YCHECK( S ) = YCDDM( I ) +! WRITE(LOGDEV,'(A,A,ES16.4)')' SORT: ' // CHEMISTRY_SPC(I),' = ',YCHECK(S) +! END DO +! DO I = 1,NUMB_MECH_SPC +! S = JNEW2OLD( I,1 ) +! WRITE(LOGDEV,'(A,A,2(ES16.4,1X))')' ORIGINAL: ' // CHEMISTRY_SPC(S),' = ',YCDDM( S ),YCHECK(I)-YCDDM( S ) +! END DO +! YCDDM = YCHECK + + END IF + +! DO RXN = 1, NRXNS + +! WRITE(LOGDEV,'(A,1X,I4,1X,A)')'For reaction number and label,',RXN, +! & TRIM(RXLABEL(RXN)) +! WRITE(LOGDEV,*)'Reactants' +! WRITE(LOGDEV,'(10X,40(A,", "))') +! & (CHEMISTRY_SPC(JNEW2OLD(ISPECIES_REACTION( NRCT,RXN ),1)),NRCT = 1, NREACT( RXN )) + +! WRITE(LOGDEV,*)'Products' +! WRITE(LOGDEV,'(10X,40(A,", "))') +! & (CHEMISTRY_SPC(JNEW2OLD(ISPECIES_REACTION( NPROD+3,RXN ),1)),NPROD = 1, NPRDCT( RXN )) + +! END DO +! WRITE(LOGDEV,*)'Original order: chemistry species, cgrid index' +! DO S = 1, NUMB_MECH_SPC +! WRITE(LOGDEV,'(A16,1X,I4)')CHEMISTRY_SPC(S),CGRID_INDEX(S) +! END DO +! WRITE(LOGDEV,*)'Sort ordered: chemistry species, cgrid index' +! DO S = 1, NUMB_MECH_SPC +! WRITE(LOGDEV,'(A16,1X,I4)')CHEMISTRY_SPC(JOLD2NEW(S,1)),MECH_INDEX(S) +! END DO END IF + + IF ( REORDER ) THEN ! reorder to speed-up LU decomposition of Jacobian + + DO I = 1,NUMB_MECH_SPC + S = JOLD2NEW( I,1 ) + YCHECK( S ) = YCDDM( I ) + END DO + YCDDM = YCHECK + + END IF CALL EVALUATE_F_JAC_MECH( YCDDM, SRK, PDT ) ! Evaluate Jacobian based on YCDDM and SKR values DO 433 J = 1, NUMB_MECH_SPC DO 434 I = 1, NUMB_MECH_SPC - A( I, J ) = 0.0 + A( I, J ) = 0.0 A1( I, J ) = 0.0 - A( I, J ) = -0.5 * CHEMSTEP * PDT( I, J ) - IF ( I .EQ. J ) THEN - A( I, J ) = 1.0 + A( I, J ) - END IF - A1( I, J ) = 0.5 * CHEMSTEP * PDT( I, J ) + A( I, J ) = -0.5 * CHEMSTEP * PDT( I, J ) + A1( I, J ) = 0.5 * CHEMSTEP * PDT( I, J ) IF ( I .EQ. J ) THEN + A( I, J ) = 1.0 + A( I, J ) A1( I, J ) = 1.0 + A1( I, J ) END IF 434 CONTINUE @@ -195,10 +330,9 @@ Subroutine SOLVE_DDM3D_CHEM( C,R,L, CHEMSTEP) DO S = 1, NUMB_MECH_SPC PRD( S ) = 0.0 PRD_RATE( S ) = 0.0 - SOLD(S) = SENGRID( C, R, L, NP, CGRID_INDEX(S) ) - IF (ABS(SOLD( S )) .LT. 1.e-25 ) THEN - SOLD(S) = 0. - END IF + J = JOLD2NEW(S,1) + SOLD(J) = SENGRID( C, R, L, NP,CGRID_INDEX(S)) + IF ( ABS(SOLD( J )) .LT. 1.e-25 ) SOLD(J) = 0.0 END DO C Begin code specific to reaction rate sensitivities @@ -215,24 +349,24 @@ Subroutine SOLVE_DDM3D_CHEM( C,R,L, CHEMSTEP) ! and slots 4- are for products IF ( NREACT( RXN ) .EQ. 1 ) THEN RXNMULT = KSTEP - & * YCDDM( IRR( RXN, 1 ) ) + & * YCDDM( ISPECIES_REACTION( 1,RXN ) ) ELSE IF ( NREACT( RXN ) .EQ. 2 ) THEN RXNMULT = KSTEP - & * YCDDM( IRR( RXN, 1 ) ) - & * YCDDM( IRR( RXN, 2 ) ) + & * YCDDM( ISPECIES_REACTION( 1,RXN ) ) + & * YCDDM( ISPECIES_REACTION( 2,RXN ) ) ELSE IF ( NREACT( RXN ) .EQ. 3 ) THEN RXNMULT = KSTEP - & * YCDDM( IRR( RXN, 1 ) ) - & * YCDDM( IRR( RXN, 2 ) ) - & * YCDDM( IRR( RXN, 3 ) ) + & * YCDDM( ISPECIES_REACTION( 1,RXN ) ) + & * YCDDM( ISPECIES_REACTION( 2,RXN ) ) + & * YCDDM( ISPECIES_REACTION( 3,RXN ) ) ELSE XMSG = 'NREACT out of expected bounds of 1-3.' CALL M3EXIT( PNAME, 0, 0, XMSG, XSTAT1 ) END IF DO NRCT = 1, NREACT( RXN ) ! Loop over the reactants - PRD_RATE( IRR( RXN, NRCT ) ) = PRD_RATE( IRR( RXN, NRCT ) ) - & - RXNMULT * IREGTEMP ! Subtract RXNMULT from PRD_RATE for reactant species + PRD_RATE( ISPECIES_REACTION( NRCT,RXN ) ) = PRD_RATE( ISPECIES_REACTION( NRCT,RXN ) ) + & - RXNMULT * IREGTEMP ! Subtract RXNMULT from PRD_RATE for reactant species END DO DO NPROD = 1, NPRDCT( RXN ) ! Loop over the products @@ -244,35 +378,36 @@ Subroutine SOLVE_DDM3D_CHEM( C,R,L, CHEMSTEP) ! coefficients of products ! and do not need the +3 (see ! RXDT.EXT) - PRD_RATE( IRR( RXN, 3+NPROD ) ) = PRD_RATE( IRR( RXN, 3+NPROD ) ) - & + ( RXNMULT * SC( RXN, NPROD ) * IREGTEMP ) + PRD_RATE( ISPECIES_REACTION( 3+NPROD,RXN ) ) = PRD_RATE( ISPECIES_REACTION( 3+NPROD,RXN ) ) + & + ( RXNMULT * SC( RXN, NPROD ) * IREGTEMP ) END DO END IF END DO ! RXN END IF ! RXNFLAG C End code specific to reaction rate sensitivities DO S = 1, NUMB_MECH_SPC - SUM = 0.0 + TOTAL = 0.0D0 DO J = 1, NUMB_MECH_SPC - SUM = SUM + A1( S, J ) * SOLD( J ) + TOTAL = TOTAL + A1( S, J ) * SOLD( J ) END DO C edits by A.Digar - PRD( S ) = SUM + PRD_RATE( S ) + PRD( S ) = TOTAL + PRD_RATE( S ) C end edit END DO CALL SGESL( A, NUMB_MECH_SPC, NUMB_MECH_SPC, IPVT, PRD, 0 ) DO S = 1, NUMB_MECH_SPC + J = CGRID_INDEX(JNEW2OLD(S,1)) IF ( ABS ( PRD ( S ) ) .LT. 1.e-25 ) THEN IF ( HIGH ) THEN - SMID( NP, S ) = 0.5 * SENGRID( C, R, L, NP, CGRID_INDEX(S) ) + SMID( S,NP ) = 0.5 * SOLD(S) END IF - SENGRID( C, R, L, NP, CGRID_INDEX(S) ) = 0. + SENGRID( C, R, L, NP,J ) = 0.0 ELSE IF ( HIGH ) THEN ! SMID is the average of SENGRID before and after chemistry - SMID( NP, S ) = 0.5 * ( SENGRID( C, R, L, NP, CGRID_INDEX(S) ) + PRD( S) ) + SMID( S,NP ) = 0.5 * ( SOLD(S) + PRD( S) ) END IF - SENGRID( C, R, L, NP, CGRID_INDEX(S) ) = PRD( S ) + SENGRID( C, R, L, NP,J ) = PRD( S ) END IF END DO @@ -291,15 +426,14 @@ Subroutine SOLVE_DDM3D_CHEM( C,R,L, CHEMSTEP) END DO DO S = 1, NUMB_MECH_SPC - SMIDJAC( S ) = SMID( HIPARM( 1 ), S ) + SMIDJAC( S ) = SMID( S,HIPARM( 1 ) ) END DO C Added by A.Digar DO S = 1, NUMB_MECH_SPC - PRD( S ) = 0.0 - SOLD(S) = SENGRID( C, R, L,NP,CGRID_INDEX(S) ) - IF (ABS(SOLD( S )) .LT. 1.e-25 ) THEN - SOLD(S) = 0. - END IF + PROD(S) = 0.0 + J = JOLD2NEW(S,1) + SOLD(J) = SENGRID( C, R, L, NP, CGRID_INDEX(S) ) + IF (ABS(SOLD( J )) .LT. 1.e-25 ) SOLD(J) = 0.0 END DO cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc @@ -311,6 +445,7 @@ Subroutine SOLVE_DDM3D_CHEM( C,R,L, CHEMSTEP) HITMP1 = HIPARM( 1 ) HITMP2 = HIPARM( 2 ) + DO N = 1, 2 ! loop for two 1st-order sens parameters IF ( ( IPT( HITMP1 ) .EQ. 5 ) .AND. ( RXNFLAG( HITMP1 ) ) ) THEN ! check for rate constant sens, date & time IREGTEMP = IREGION ( C, R, L, HITMP1 ) @@ -327,7 +462,7 @@ Subroutine SOLVE_DDM3D_CHEM( C,R,L, CHEMSTEP) CALL EVALUATE_F_JAC_MECH( YCDDM, RK, PDT ) ! Evaluate Jacobian based on YCDDM and RK values DO S = 1, NUMB_MECH_SPC DO J = 1, NUMB_MECH_SPC - PRD( S ) = PRD( S ) + CHEMSTEP * PDT( S,J ) * SMID( HITMP2,J ) * IREGTEMP + PRD( S ) = PRD( S ) + CHEMSTEP * PDT( S,J ) * SMID( J,HITMP2 ) * IREGTEMP END DO END DO IF ( IPT( HITMP1 ) .eq. IPT( HITMP2 ) ) THEN @@ -346,13 +481,13 @@ Subroutine SOLVE_DDM3D_CHEM( C,R,L, CHEMSTEP) CALL EVALUATE_F_JAC_MECH( SMIDJAC, SRK2, PDT2 ) ! Evaluate Jacobian based on SMIDJAC and SKR2 values DO S = 1, NUMB_MECH_SPC - SUM = 0.0 + TOTAL = 0.0 DO J = 1, NUMB_MECH_SPC - SUM = SUM + A1( S, J ) * SOLD( J ) - & + CHEMSTEP * PDT2( S,J ) * SMID( HIPARM( 2 ),J ) + TOTAL = TOTAL + A1( S, J ) * SOLD( J ) + & + CHEMSTEP * PDT2( S,J ) * SMID( J,HIPARM( 2 ) ) END DO C edits by A.Digar - PRD( S ) = SUM + PRD( S ) + PRD( S ) = TOTAL + PRD( S ) C end of edits END DO @@ -360,9 +495,9 @@ Subroutine SOLVE_DDM3D_CHEM( C,R,L, CHEMSTEP) DO S = 1, NUMB_MECH_SPC IF ( ABS ( PRD ( S ) ) .LT. 1.e-25 ) THEN - SENGRID( C, R, L, NP, CGRID_INDEX(S) ) = 0. + SENGRID( C, R, L, NP, CGRID_INDEX(JNEW2OLD(S,1)) ) = 0.0 ELSE - SENGRID( C, R, L, NP, CGRID_INDEX(S) ) = PRD( S ) + SENGRID( C, R, L, NP, CGRID_INDEX(JNEW2OLD(S,1)) ) = PRD( S ) END IF END DO diff --git a/CCTM/src/ddm3d/MECHANISM_FUNCTIONS.F b/CCTM/src/ddm3d/MECHANISM_FUNCTIONS.F index 7b2cd81db..c5862a711 100644 --- a/CCTM/src/ddm3d/MECHANISM_FUNCTIONS.F +++ b/CCTM/src/ddm3d/MECHANISM_FUNCTIONS.F @@ -17,6 +17,7 @@ MODULE MECHANISM_FUNCTIONS REAL( 8 ), ALLOCATABLE :: DYDT( : ) ! time derivative of species REAL( 8 ), ALLOCATABLE :: JACOBIAN( :,: ) ! mechanism's jacobain matrix + INTEGER, PARAMETER :: MAX_NCELLS = 1 CHARACTER( 16 ), ALLOCATABLE :: SPECIES( : ) @@ -65,6 +66,7 @@ MODULE MECHANISM_FUNCTIONS INTEGER :: DDM_LOG = 6 ! Unit number of output log INTEGER :: ERROR_LOG = 6 LOGICAL :: CHECK_MECHANISM = .FALSE. ! write out Jacobian and derivatives values + LOGICAL :: CHECK_SORTING = .FALSE. ! write out result of mechanism species sorting #ifdef verbose_ddm3d LOGICAL :: REPORT_CHART = .TRUE. ! write out species derivations and mechanism jacobian @@ -84,8 +86,13 @@ MODULE MECHANISM_FUNCTIONS INTEGER, ALLOCATABLE :: JOLD2NEW( :,: ) ! YC species map INTEGER, ALLOCATABLE :: JNEW2OLD( :,: ) ! YC species map + INTEGER, ALLOCATABLE :: ISPECIES_REACTION( :,: ) ! species index for each reaction + + INTEGER :: CHANGING_SPECIES = 0 ! mechanism species affected by reactions + INTEGER :: STEADY_SPECIES = 0 ! mechanism species not affected by reactions + + LOGICAL, ALLOCATABLE :: JACOBIAN_FILLED( :,:,: ) ! Is Jacobian position nonzero - LOGICAL, ALLOCATABLE :: JACOBIAN_FILLED( :,:,: ) LOGICAL :: LSUNLIGHT = .TRUE. REAL( 8 ), ALLOCATABLE, PRIVATE :: ATMPRES ( : ) ! Cell pressure, Atm @@ -96,8 +103,10 @@ MODULE MECHANISM_FUNCTIONS REAL( 8 ), ALLOCATABLE, PRIVATE :: RJIN ( :, : ) ! J-values for a cell LOGICAL, ALLOCATABLE, PRIVATE :: LAND ( : ) ! land_zone value for specific cell + LOGICAL :: RESET_JACOBIAN = .FALSE. + CONTAINS - SUBROUTINE SET_MECHANISM( ) + SUBROUTINE SET_MECHANISM( ) USE RXNS_DATA USE RXNS_FUNCTION @@ -118,20 +127,37 @@ SUBROUTINE SET_MECHANISM( ) C..Parameters: - INTEGER :: I, J, RXN, IP, IPNEG, IL + INTEGER :: I, J, RXN, IP, IPNEG, IL, ISPC INTEGER :: NCELL INTEGER :: IOS INTEGER :: C, L, R, S ! Loop indices INTEGER :: SPC ! array index + INTEGER :: MIN_NEW_SPECIES ! Index holder for sort routine + INTEGER :: MIN_OLD_SPECIES ! Index holder for sort routine + INTEGER :: MIN_COUNT ! Current minimum number of PD terms in sort + + INTEGER INEW, JNEW ! Index for sorted species number + INTEGER IOLD, JOLD ! Index for old species number + INTEGER IMINNEW ! Index holder for sort routine + INTEGER IMINOLD ! Index holder for sort routine + INTEGER MINVALU ! Current number of PD terms in sort + LOGICAL :: EXISTS LOGICAL :: EFLAG + LOGICAL :: SWAPPING + LOGICAL, SAVE :: INITIALIZED = .FALSE. CHARACTER( 132 ) :: MSG ! Message text + CHARACTER( 32 ) :: SWAPPED_NAME REAL :: HEIGHT REAL( 8 ) :: FACTOR ! conversion factor + INTEGER, ALLOCATABLE :: SORTED_COUNT ( : ) ! # of Jacobain partial derivations per species + INTEGER, ALLOCATABLE :: JACOBIAN_COUNT( : ) ! # of Jacobain partial derivations per species + CHARACTER( 32 ), ALLOCATABLE :: SWAPPED_NAMES( : ) + IF ( INITIALIZED ) RETURN EFLAG = .FALSE. ERROR_LOG = DDM_LOG @@ -174,17 +200,6 @@ SUBROUTINE SET_MECHANISM( ) WRITE(ERROR_LOG,*)TRIM(MSG) EFLAG = .TRUE. END IF - ALLOCATE( JOLD2NEW( NUMB_MECH_SPC,1 ), - & JNEW2OLD( NUMB_MECH_SPC,1 ), STAT = IOS ) - IF ( IOS .NE. 0 ) THEN - MSG = 'Error allocating JOLD2NEW or JNEW2OLD' - WRITE(ERROR_LOG,*)TRIM(MSG) - EFLAG = .TRUE. - END IF - DO I = 1, NUMB_MECH_SPC - JOLD2NEW( I,1 ) = I - JNEW2OLD( I,1 ) = I - END DO ALLOCATE( RKI_SAV( NRXNS,MAX_NCELLS ), STAT = IOS ) IF ( IOS .NE. 0 ) THEN MSG = 'Error allocating RKI_SAV' @@ -193,6 +208,51 @@ SUBROUTINE SET_MECHANISM( ) END IF END IF + ALLOCATE( JOLD2NEW( NUMB_MECH_SPC,1 ), + & JNEW2OLD( NUMB_MECH_SPC,1 ), STAT = IOS ) + IF ( IOS .NE. 0 ) THEN + MSG = 'Error allocating JOLD2NEW or JNEW2OLD' + WRITE(ERROR_LOG,*)TRIM(MSG) + EFLAG = .TRUE. + ELSE + MSG = 'JOLD2NEW or JNEW2OLD allocated' + WRITE(ERROR_LOG,*)TRIM(MSG) + END IF + DO I = 1, NUMB_MECH_SPC + JOLD2NEW( I,1 ) = I + JNEW2OLD( I,1 ) = I + END DO + + ALLOCATE( SORTED_COUNT( NUMB_MECH_SPC), + & JACOBIAN_COUNT( NUMB_MECH_SPC), STAT = IOS ) + IF ( IOS .NE. 0 ) THEN + MSG = 'Error allocating JACOBIAN_COUNT' + WRITE(ERROR_LOG,*)TRIM(MSG) + EFLAG = .TRUE. + END IF + JACOBIAN_COUNT = 0 + SORTED_COUNT = 0 + + ALLOCATE( SWAPPED_NAMES( NUMB_MECH_SPC), STAT = IOS ) + IF ( IOS .NE. 0 ) THEN + MSG = 'Error allocating SWAPPED_NAMES' + WRITE(ERROR_LOG,*)TRIM(MSG) + EFLAG = .TRUE. + END IF + SWAPPED_NAMES = 'BLANK' + + ALLOCATE( ISPECIES_REACTION(MXPRD+3,NRXNS), STAT = IOS ) + IF ( IOS .NE. 0 ) THEN + MSG = 'Error allocating ISPECIES_REACTION' + WRITE(ERROR_LOG,*)TRIM(MSG) + EFLAG = .TRUE. + END IF + DO RXN = 1,NRXNS + DO SPC = 1,MXPRD+3 + ISPECIES_REACTION( SPC,RXN ) = IRR( RXN,SPC ) + END DO + END DO + ALLOCATE( DYDT( NUMB_MECH_SPC ), STAT = IOS ) IF ( IOS .NE. 0 ) THEN MSG = 'Error allocating DYDT array' @@ -208,7 +268,142 @@ SUBROUTINE SET_MECHANISM( ) ! map species time derivative and jacobian array CALL CHART_IRR - + +! Set the number of Partial derivative terms in the Jacobian and +! count the number of terms for each species + DO I = 1, NUMB_MECH_SPC + DO J = 1, NUMB_MECH_SPC + IF ( JACOBIAN_FILLED( J,I,1 ) ) THEN + JACOBIAN_COUNT( J ) = JACOBIAN_COUNT( J ) + 1 + END IF + END DO + END DO + +! Sort the species, putting all with zero partial derivative +! terms at the bottom and those with fewest PD terms at top. +! Set arrays for species with zero PD terms + STEADY_SPECIES = NUMB_MECH_SPC + DO JOLD = 1, NUMB_MECH_SPC + IF ( JACOBIAN_COUNT( JOLD ) .GT. 0 ) THEN + CHANGING_SPECIES = CHANGING_SPECIES + 1 + JNEW = CHANGING_SPECIES + JNEW2OLD( JNEW,1 ) = JOLD + JOLD2NEW( JOLD,1 ) = JNEW + ELSE + JNEW2OLD( CHANGING_SPECIES,1 ) = JOLD + JOLD2NEW( JOLD,1 ) = STEADY_SPECIES + STEADY_SPECIES = STEADY_SPECIES - 1 + END IF + END DO + STEADY_SPECIES = NUMB_MECH_SPC-CHANGING_SPECIES + IF (CHECK_SORTING )THEN + WRITE( DDM_LOG,'(A)')'Initial Mechanism Species sorted based on their mechanism activity' + WRITE( DDM_LOG,'(A,2(I4,1X))')'Changeing Species, Steady State Species ',CHANGING_SPECIES, + & STEADY_SPECIES + WRITE( DDM_LOG,'(A16,1X,4(A16,1X))')'ISPECIES', 'Original','Sorted','Check','JCOUNT' + DO ISPC = 1, CHANGING_SPECIES + I = JOLD2NEW( ISPC,1 ) + J = JNEW2OLD( I,1 ) + WRITE(DDM_LOG,'(I4,13X,3(A16,1X),I16)')ISPC, + & CHEMISTRY_SPC(ISPC),CHEMISTRY_SPC(I),CHEMISTRY_SPC(J),JACOBIAN_COUNT(I) + END DO + END IF + +! Now sort by number of PD terms, fewest at position 1, most at +! the end position. + DO JNEW = 1, CHANGING_SPECIES +c Uncomment the following three lines to turn off species ordering; +c not recommended since computational efficiency reduced +! INEW2OLD( JNEW,NCS ) = JNEW +! IOLD2NEW( JNEW,NCS ) = JNEW +! IF ( JNEW .NE. 0 ) CYCLE + JOLD = JNEW2OLD( JNEW,1 ) + MINVALU = JACOBIAN_COUNT( JOLD ) + IMINOLD = JOLD + IMINNEW = JNEW + + DO INEW = JNEW + 1, CHANGING_SPECIES + IOLD = JNEW2OLD( INEW,1 ) + IF ( JACOBIAN_COUNT( IOLD ) .LT. MINVALU ) THEN + MINVALU = JACOBIAN_COUNT( IOLD ) + IMINOLD = IOLD + IMINNEW = INEW + END IF + END DO + + JNEW2OLD( IMINNEW,1 ) = JOLD + JNEW2OLD( JNEW,1 ) = IMINOLD + JOLD2NEW( JOLD,1 ) = IMINNEW + JOLD2NEW( IMINOLD,1 ) = JNEW + + END DO + DO J = 1,NUMB_MECH_SPC + JNEW = JOLD2NEW( J,1 ) + SORTED_COUNT(JNEW) = JACOBIAN_COUNT( J ) + END DO + JACOBIAN_COUNT = SORTED_COUNT +! Alternative method for sorting species for species names +! SWAPPED_NAMES = CHEMISTRY_SPC +! DO J = CHANGING_SPECIES-1,1,-1 +! SWAPPING = .FALSE. +! DO I = 1,J +! IF ( JACOBIAN_COUNT(I) .GT. JACOBIAN_COUNT(I+1) )THEN +! MINVALU = JACOBIAN_COUNT(I+1) +! JACOBIAN_COUNT(I+1) = JACOBIAN_COUNT(I) +! JACOBIAN_COUNT(I) = MINVALU +! SWAPPED_NAME = SWAPPED_NAMES(I+1) +! SWAPPED_NAMES(I+1) = SWAPPED_NAMES(I) +! SWAPPED_NAMES(I) = SWAPPED_NAME +! SWAPPING = .TRUE. +! END IF +! END DO +! IF( .NOT. SWAPPING ) EXIT +! END DO +! DO I = 1, NUMB_MECH_SPC +! WRITE(DDM_LOG,'(I4,1X,A)')I,SWAPPED_NAMES(I) +! END DO +! DO J = 1,CHANGING_SPECIES +! DO I = 1,CHANGING_SPECIES +! IF( TRIM(CHEMISTRY_SPC(J)) .EQ. TRIM( SWAPPED_NAMES(I) ) )THEN +! JOLD2NEW(J,1) = I +! JNEW2OLD(I,1) = J +! WRITE(DDM_LOG,'(I4,1X,A,2(1X,I4))')J,CHEMISTRY_SPC(J),JOLD2NEW(J,1),JNEW2OLD(I,1) +! END IF +! END DO +! END DO + IF (CHECK_SORTING )THEN +! Write sorting results + WRITE( DDM_LOG,'(A)')'Mechanism Species sorted based on their mechanism activity' + WRITE( DDM_LOG,'(A,2(I4,1X))')'Changeing Species, Steady State Species ',CHANGING_SPECIES, + & STEADY_SPECIES + WRITE( DDM_LOG,'(A16,1X,4(A16,1X))')'ISPECIES', 'Original','Sorted','Check','JCOUNT' + DO ISPC = 1, NUMB_MECH_SPC + I = JOLD2NEW( ISPC,1 ) + J = JNEW2OLD( I,1 ) + S = JNEW2OLD( ISPC,1 ) + WRITE(DDM_LOG,'(I4,13X,3(A16,1X),I16)')ISPC, + & CHEMISTRY_SPC(ISPC),CHEMISTRY_SPC(S),CHEMISTRY_SPC(J),JACOBIAN_COUNT(ISPC) + + END DO + END IF + + RESET_JACOBIAN = .TRUE. + +! Reset the ISPECIES_REACTION array using the new species order developed above. + + DO RXN = 1, NRXNS + DO I = 1, NREACT( RXN ) + J = ISPECIES_REACTION( I,RXN ) + ISPECIES_REACTION( I,RXN ) = JOLD2NEW( J,1 ) + END DO + + DO I = 1, NPRDCT( RXN ) + J = ISPECIES_REACTION( I+3,RXN ) + ISPECIES_REACTION( I+3,RXN ) = JOLD2NEW( J,1 ) + END DO + END DO + DEALLOCATE( JACOBIAN_COUNT, SORTED_COUNT ) + IF( CHECK_MECHANISM )CHECK_MECHANISM = .FALSE. 95000 FORMAT(I4,1X,A16,' = ',4(ES12.4,', ')) @@ -291,7 +486,7 @@ SUBROUTINE CHART_IRR() CALL M3EXIT ( PNAME, 0, 0, MSG, XSTAT2 ) END IF JACOBIAN_FILLED = .FALSE. - DO I = 1, NUMB_MECH_SPC + DO I = 1, NUMB_MECH_SPC ! set diagonal to to true JACOBIAN_FILLED(I,I,1:2) = .TRUE. END DO @@ -591,6 +786,7 @@ SUBROUTINE EVALUATE_F_JAC_MECH( YIN, RKI, JAC ) INTEGER NP ! Loop index over partial derivation terms INTEGER IPART ! index for partial derivation of reaction INTEGER NRK ! Reaction number + INTEGER IOS REAL( 8 ) :: CR2 ! Temporary product for 3 reactant reaction REAL( 8 ) :: FRACN ! Stoichiometric coefficient @@ -598,8 +794,40 @@ SUBROUTINE EVALUATE_F_JAC_MECH( YIN, RKI, JAC ) LOGICAL, SAVE :: INITIALIZE = .TRUE. + INTEGER, ALLOCATABLE :: JACOBIAN_COL( : ) + INTEGER, ALLOCATABLE :: JACOBIAN_ROW( : ) + CHARACTER(16), ALLOCATABLE :: JACOBIAN_COL_SPECIES( : ) + CHARACTER(16), ALLOCATABLE :: JACOBIAN_ROW_SPECIES( : ) + C*********************************************************************** + IF( RESET_JACOBIAN )THEN ! reroder row and columns of Jacobian + ! using sorted species maps + + DO NRK = 1, NRXNS + DO NP = 1, REACTION_CHART( NRK )%JACOB_OCCURANCES + JROW = REACTION_CHART( NRK )%JACOB_PARTIAL_ROW( NP ) + JCOL = REACTION_CHART( NRK )%JACOB_PARTIAL_COL( NP ) + JROW = JOLD2NEW( JROW,1 ) + JCOL = JOLD2NEW( JCOL,1 ) + REACTION_CHART( NRK )%JACOB_PARTIAL_ROW( NP ) = JROW + REACTION_CHART( NRK )%JACOB_PARTIAL_COL( NP ) = JCOL + END DO + END DO + + DO NRK = 1, NRXNS + DO NP = 1,REACTION_CHART( NRK )%NREACTANTS + JR1 = JOLD2NEW( REACTION_CHART( NRK )%REACTANT( NP ),1 ) + REACTION_CHART( NRK )%REACTANT( NP ) = JR1 + END DO + END DO + + RESET_JACOBIAN = .FALSE. + + END IF + IF( INITIALIZE )THEN + INITIALIZE = .FALSE. + END IF c...initialize Jacobian JAC( :,: ) = 0.0 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc @@ -635,14 +863,11 @@ SUBROUTINE EVALUATE_F_JAC_MECH( YIN, RKI, JAC ) cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Add Reaction's Partial Derivative to Jacobian cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - DO NP = 1, REACTION_CHART( NRK )%JACOB_OCCURANCES - + DO NP = 1, REACTION_CHART( NRK )%JACOB_OCCURANCES JROW = REACTION_CHART( NRK )%JACOB_PARTIAL_ROW( NP ) JCOL = REACTION_CHART( NRK )%JACOB_PARTIAL_COL( NP ) FRACN = REACTION_CHART( NRK )%JACOB_PARTIAL_COEFF( NP ) IPART = REACTION_CHART( NRK )%JACOB_PARTIAL_INDEX( NP ) - !JROW = JOLD2NEW( JROW,1 ) - !JCOL = JOLD2NEW( JCOL,1 ) JAC( JROW, JCOL ) = JAC( JROW, JCOL ) + REAL(FRACN * EXPLIC( IPART ) ) END DO END DO LOOP_REACTIONS @@ -718,6 +943,7 @@ SUBROUTINE EVALUATE_F_MECH( YIN, TAIR, DAIR, RKI, YDOT ) ! YDOT = 0.0D0 PROD = 0.0D0 LOSS = 0.0D0 + cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Loop over reactions to calculate dc/dt cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc diff --git a/CCTM/src/gas/ros3/rbdata_mod.F b/CCTM/src/gas/ros3/rbdata_mod.F index 23d3c56e9..efd9d23a6 100644 --- a/CCTM/src/gas/ros3/rbdata_mod.F +++ b/CCTM/src/gas/ros3/rbdata_mod.F @@ -87,6 +87,9 @@ MODULE RBDATA REAL( 8 ), ALLOCATABLE :: REVERSE_CONV( : ) ! CHEM to CGRID Species conversion factor REAL( 8 ), ALLOCATABLE :: Y( :,: ) ! Species concentrations REAL( 8 ), ALLOCATABLE :: Y_DEGRADE( :,: ) ! Concentration for degradation +#ifdef sens + REAL( 8 ), ALLOCATABLE :: YAVE( :,: ) ! Species concentrations +#endif c..Block variables INTEGER :: BLKID ! Block ID diff --git a/CCTM/src/gas/ros3/rbdriver.F b/CCTM/src/gas/ros3/rbdriver.F index 6e4cd49d3..1edeefc98 100644 --- a/CCTM/src/gas/ros3/rbdriver.F +++ b/CCTM/src/gas/ros3/rbdriver.F @@ -83,6 +83,10 @@ SUBROUTINE CHEM( CONC, JDATE, JTIME, TSTEP ) USE PA_DEFN, Only: LIRR ! Process Anaylsis control and data variable USE PA_IRR_CLT USE CENTRALIZED_IO_MODULE, ONLY : INTERPOLATE_VAR, OCEAN, SZONE +#ifdef sens + USE DDM3D_CHEM + Use DDM3D_DEFN, Only: DATENUM, STARTDATE, IPT, IDATE, HIGH, NP, NPMAX, CKTIME +#endif IMPLICIT NONE @@ -154,6 +158,8 @@ SUBROUTINE CHEM( CONC, JDATE, JTIME, TSTEP ) INTEGER ROW ! Row index INTEGER SPC ! Species loop index INTEGER VAR ! Variable number on I/O API file + + INTEGER NUMB_CELLS INTEGER, ALLOCATABLE, SAVE :: IRRCELL( : ) ! Cell No. of an IRR cell @@ -192,25 +198,6 @@ SUBROUTINE RBSOLVER ( JDATE, JTIME, CHEMSTEP, NCSP, INTEGER, INTENT( INOUT ) :: NIRRCLS INTEGER, INTENT( IN ) :: IRRCELL( : ) END SUBROUTINE RBSOLVER -! SUBROUTINE FIND_DEGRADED( JDATE, JTIME, CALL_DEGRADE ) -! INTEGER, INTENT( IN ) :: JDATE ! current model date , coded YYYYDDD -! INTEGER, INTENT( IN ) :: JTIME ! current model time , coded HHMMSS -! LOGICAL, INTENT( OUT ) :: CALL_DEGRADE ! whether to call degradation routines -! END SUBROUTINE FIND_DEGRADED -! SUBROUTINE INIT_DEGRADE_BLK( CBLK, TCELL, DCELL, PRESS_CELL, QV_CELL, PHOTO_CELL, -! & JDATE, JTIME ) -! REAL( 8 ), INTENT( IN ) :: CBLK( :,: ) ! species concentration in cell -! REAL( 8 ), INTENT( IN ) :: TCELL( : ) ! cell temperature [ k ] -! REAL( 8 ), INTENT( IN ) :: DCELL( : ) ! cell air density [ kg/m^3 ] -! REAL( 8 ), INTENT( IN ) :: PRESS_CELL( : ) ! cell Pressure [ Pa ] -! REAL( 8 ), INTENT( IN ) :: QV_CELL( : ) ! cell water vapor mass mixing ratio [ kg/kg ] -! REAL( 8 ), INTENT( IN ) :: PHOTO_CELL( :,: ) ! Photolysis table for cell [1/s] -! INTEGER, INTENT( IN ) :: JDATE ! current model date , coded YYYYDDD -! INTEGER, INTENT( IN ) :: JTIME ! current model time , coded HHMMSS -! END SUBROUTINE INIT_DEGRADE_BLK -! SUBROUTINE FINAL_DEGRADE_BLK( CBLK ) -! REAL( 8 ), INTENT( INOUT ) :: CBLK( :,: ) ! species conc in cell -! END SUBROUTINE FINAL_DEGRADE_BLK SUBROUTINE HETCHEM_UPDATE_AERO( CGRID ) REAL, POINTER :: CGRID( :,:,:,: ) ! species concentration in cell END SUBROUTINE HETCHEM_UPDATE_AERO @@ -226,12 +213,6 @@ END SUBROUTINE HETCHEM_UPDATE_AERO CALL M3EXIT( PNAME, JDATE, JTIME, MSG, XSTAT1 ) #endif -#ifdef sens - MSG = 'ERROR: Rosenbrock Chemistry Solver does not perform DDM3 calculation.' - WRITE(LOGDEV,'(A)')TRIM( MSG ) - MSG = 'Must use the EBI solver for the chemical mechanism' - CALL M3EXIT( PNAME, JDATE, JTIME, MSG, XSTAT1 ) -#endif IF ( NUMB_MECH_SPC .EQ. 0 ) THEN CALL M3MESG( '*** WARNING: Gas-Chemistry processing bypassed!' ) @@ -366,6 +347,20 @@ END SUBROUTINE HETCHEM_UPDATE_AERO ALLOCATE( Y_DEGRADE( BLKSIZE, NSPCSD ) ) +#ifdef sens + CALL INIT_DDM3D_CHEM() + +C For higher order sensitivities + IF ( HIGH ) THEN + DO IRXN = 1, NRXNS + IF( NREACT( IRXN ) .LE. 1 ) THEN + ORDER1( IRXN ) = .TRUE. + ELSE + ORDER1( IRXN ) = .FALSE. + END IF + END DO + END IF +#endif END IF ! First call ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc @@ -406,6 +401,21 @@ END SUBROUTINE HETCHEM_UPDATE_AERO CALL HETCHEM_RATES( TA, PRES, QV, CONC, DENS ) +#ifdef sens +C Set the date and hour counters used in sensitivity calls + DATENUM = 1 + JDATE - STARTDATE + +C For reaction rate sensitivities + DO NP = 1, NPMAX + IF ( IPT( NP ) .EQ. 5 ) THEN + CALL CKTIME( JDATE,JTIME,NP,RXNFLAG(NP) ) ! Rxnflag set to true if ipt=5 and time, date within bounds + IF ( IDATE( NP, DATENUM ) .NE. 1 ) RXNFLAG( NP ) = .FALSE. + ELSE + RXNFLAG( NP ) = .FALSE. + END IF + END DO +#endif + ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Set flag for reordering of cells and put cells in sequential c order initially @@ -581,17 +591,6 @@ END SUBROUTINE HETCHEM_UPDATE_AERO IF ( CALL_DEG ) THEN - -C Update degradation array with species treated by Rosenbach solver -C -C DO ISP = 1, ISCHANG( NCS ) -C ISPOLD = INEW2OLD( ISP, NCS ) -C DO NCELL = 1, NUMCELLS -C Y_DEGRADE( NCELL,ISPOLD ) = Y( NCELL,ISP ) -C END DO -C END DO - - C Update CGRID based on the degradation routines CALL FINAL_DEGRADE_BLK( Y_DEGRADE ) UPDATE_DEGRADED: DO ISP = 1, N_REACT @@ -610,6 +609,42 @@ END SUBROUTINE HETCHEM_UPDATE_AERO END DO UPDATE_DEGRADED END IF !WTH +#ifdef sens + + NUMB_CELLS = NUMCELLS + + DO NCELL = 1, NUMB_CELLS + DO IRXN = 1, NRXNS + SRK( IRXN ) = RKI( NCELL,IRXN ) + IF ( HIGH ) THEN + IF ( NREACT( IRXN ) .LE. 1 ) THEN + SRK2( IRXN ) = 0.0 + ELSE IF ( NREACT( IRXN ) .EQ. 2 ) THEN + SRK2( IRXN ) = REAL( RKI( NCELL,IRXN ), 4 ) + ELSE IF ( NREACT( IRXN ) .EQ. 3 ) THEN + SRK2( IRXN ) = REAL( RKI( NCELL,IRXN ),4 ) + ELSE + SRK2( IRXN ) = 0.0 + END IF + END IF + END DO + + DO ISP = 1,NUMB_MECH_SPC +! ISPOLD = INEW2OLD( ISP,NCS ) +! YCDDM( ISPOLD ) = YAVE(NCELL,ISP) + YCDDM( ISP ) = YAVE(NCELL,ISP) + END DO + CELLNUM = NORDCELL( OFFSET + NCELL ) + COL = CCOL( CELLNUM ) + ROW = CROW( CELLNUM ) + LEV = CLEV( CELLNUM ) + +! CALL SOLVE_DDM3D_CHEM( COL,ROW,LEV,CHEMSTEP ) + CALL SOLVE_DDM3D_CHEM( COL,ROW,LEV,CHEMSTEP,IOLD2NEW,INEW2OLD ) + END DO + NUMCELLS = NUMB_CELLS +#endif + #ifdef rbstats DO NCELL = 1, NUMCELLS diff --git a/CCTM/src/gas/ros3/rbinit.F b/CCTM/src/gas/ros3/rbinit.F index 8ee55d5a9..44e7c9060 100644 --- a/CCTM/src/gas/ros3/rbinit.F +++ b/CCTM/src/gas/ros3/rbinit.F @@ -238,6 +238,13 @@ SUBROUTINE RBINIT CALL M3EXIT( PNAME, 0, 0, XMSG, XSTAT2 ) END IF +#ifdef sens + ALLOCATE( YAVE( BLKSIZE,NUMB_MECH_SPC ), STAT = STATUS ) + IF ( STATUS .NE. 0 ) THEN + XMSG = 'ERROR allocating YAVE' + CALL M3EXIT( PNAME, 0, 0, XMSG, XSTAT2 ) + END IF +#endif ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Set convergence tolerances for each species; currently uses c one set of tolerances for all species diff --git a/CCTM/src/gas/ros3/rbsolver.F b/CCTM/src/gas/ros3/rbsolver.F index 3d536689b..4faff0ec9 100644 --- a/CCTM/src/gas/ros3/rbsolver.F +++ b/CCTM/src/gas/ros3/rbsolver.F @@ -338,6 +338,9 @@ SUBROUTINE RBDECOMP( NCSP ) END IF +#ifdef sens + YAVE = 0.0D0 +#endif DO 100 WHILE ( TNOW .LT. CHEMSTEP ) CALL RBFEVAL( NCSP, Y, YDOT ) @@ -462,6 +465,11 @@ SUBROUTINE RBDECOMP( NCSP ) DO N = 1, NUMB_MECH_SPC DO NCELL = 1, NUMCELLS +#ifdef sens + YAVE( NCELL,N ) = YAVE( NCELL,N ) + & + 0.5D0 * (Y( NCELL,N )+ MAX(YP( NCELL,N),CONMIN)) + & * DT +#endif Y( NCELL,N ) = MAX( YP( NCELL,N ), CONMIN ) END DO END DO @@ -580,6 +588,12 @@ SUBROUTINE RBDECOMP( NCSP ) 100 END DO ! end time integration loop + +#ifdef sens +! complete calculation for YAVE + YAVE(1:NUMCELLS,1:NUMB_MECH_SPC) = YAVE(1:NUMCELLS,1:NUMB_MECH_SPC) / CHEMSTEP +! YAVE = 0.5D0*(YAVE+Y) +#endif RETURN 92100 FORMAT( ' Convergence failure ', diff --git a/CCTM/src/gas/smvgear/GRVARS.F b/CCTM/src/gas/smvgear/GRVARS.F index 8bf607f49..2b737e670 100644 --- a/CCTM/src/gas/smvgear/GRVARS.F +++ b/CCTM/src/gas/smvgear/GRVARS.F @@ -311,6 +311,11 @@ MODULE GRVARS REAL( 8 ), ALLOCATABLE, SAVE :: BLKCONC( :, : ) ! Species conc. for cells in block ! in original species order (ppm) REAL( 8 ), ALLOCATABLE, SAVE :: CNEW( :, : ) ! Species conc. for cells in block +#ifdef sens + REAL( 8 ), ALLOCATABLE :: CAVEG( :,: ) ! Average species concentrations over time step + REAL( 8 ), ALLOCATABLE :: CINIT( :,: ) ! species concentrations at start of subtime step + REAL( 8 ), ALLOCATABLE :: CFINI( :,: ) ! species concentrations at start of subtime step +#endif ! in sorted species order (ppm) REAL( 8 ), ALLOCATABLE, SAVE :: EMBLK( :, : ) ! Species emissions in each cell REAL( 8 ), ALLOCATABLE, SAVE :: GLOSS( :, : ) ! dc/dt for each species (i.e., RHS) @@ -450,6 +455,16 @@ SUBROUTINE GRVARS_INIT( JDATE, JTIME ) CALL M3EXIT( PNAME, JDATE, JTIME, MSG, XSTAT2 ) END IF +#ifdef sens + ALLOCATE( CINIT( BLKSIZE,NUMB_MECH_SPC ), + & CFINI( BLKSIZE,NUMB_MECH_SPC ), + & CAVEG( BLKSIZE,NUMB_MECH_SPC ), STAT = IOS ) + IF ( IOS .NE. 0 ) THEN + MSG = 'ERROR allocating CINIT,CFINI,CAVEG' + CALL M3EXIT( PNAME, JDATE, JTIME, MSG, XSTAT2 ) + END IF +#endif + ALLOCATE ( CCOL( MXCELLS ), & CROW( MXCELLS ), & CLEV( MXCELLS ), diff --git a/CCTM/src/gas/smvgear/grdriver.F b/CCTM/src/gas/smvgear/grdriver.F index bd3132ed7..588012608 100644 --- a/CCTM/src/gas/smvgear/grdriver.F +++ b/CCTM/src/gas/smvgear/grdriver.F @@ -90,6 +90,10 @@ SUBROUTINE CHEM( CONC, JDATE, JTIME, TSTEP ) USE PA_DEFN, Only: LIRR ! Process Anaylsis control and data variable USE PA_IRR_CLT USE CENTRALIZED_IO_MODULE, ONLY : INTERPOLATE_VAR, OCEAN, SZONE +#ifdef sens + USE DDM3D_CHEM, DDM_RK => RK + Use DDM3D_DEFN, Only: DATENUM, STARTDATE, IPT, IDATE, HIGH, NP, NPMAX, CKTIME +#endif IMPLICIT NONE @@ -232,13 +236,6 @@ END SUBROUTINE HETCHEM_UPDATE_AERO CALL M3EXIT( PNAME, JDATE, JTIME, MSG, XSTAT1 ) #endif -#ifdef sens - MSG = 'ERROR: SMVGEAR Chemistry Solver does not perform DDM3D calculations.' - WRITE(LOGDEV,'(A)')TRIM( MSG ) - MSG = 'Must use the EBI solver for the chemical mechanism' - CALL M3EXIT( PNAME, JDATE, JTIME, MSG, XSTAT1 ) -#endif - IF ( NUMB_MECH_SPC .EQ. 0 ) RETURN ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc @@ -299,6 +296,20 @@ END SUBROUTINE HETCHEM_UPDATE_AERO ALLOCATE( Y_DEGRADE( BLKSIZE, NSPCSD ) ) +#ifdef sens + CALL INIT_DDM3D_CHEM() + +C For higher order sensitivities + IF ( HIGH ) THEN + DO IRXN = 1, NRXNS + IF( NREACT( IRXN ) .LE. 1 ) THEN + ORDER1( IRXN ) = .TRUE. + ELSE + ORDER1( IRXN ) = .FALSE. + END IF + END DO + END IF +#endif ENDIF ! First call #ifdef sunws @@ -348,6 +359,21 @@ END SUBROUTINE HETCHEM_UPDATE_AERO C solver finds a solution. CALL HETCHEM_RATES( TA, PRES, QV, CONC, DENS ) + +#ifdef sens +C Set the date and hour counters used in sensitivity calls + DATENUM = 1 + JDATE - STARTDATE + +C For reaction rate sensitivities + DO NP = 1, NPMAX + IF ( IPT( NP ) .EQ. 5 ) THEN + CALL CKTIME( JDATE,JTIME,NP,RXNFLAG(NP) ) ! Rxnflag set to true if ipt=5 and time, date within bounds + IF ( IDATE( NP, DATENUM ) .NE. 1 ) RXNFLAG( NP ) = .FALSE. + ELSE + RXNFLAG( NP ) = .FALSE. + END IF + END DO +#endif ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Set flag for reordering of cells and put cells in sequential @@ -433,6 +459,9 @@ END SUBROUTINE HETCHEM_UPDATE_AERO CNEW( NCELL,ISPNEW ) = REAL( MAX( CONC( COL,ROW,LEV,SPC ), CONCMIN), 8 ) END IF BLKCONC( NCELL, ISP ) = CNEW( NCELL,ISPNEW ) +#ifdef sens + CINIT( NCELL, ISP ) = CNEW( NCELL,ISPNEW ) +#endif ENDDO ENDDO @@ -694,6 +723,9 @@ END SUBROUTINE HETCHEM_UPDATE_AERO ROW = CROW( CELLNUM ) COL = CCOL( CELLNUM ) LEV = CLEV( CELLNUM ) +#ifdef sens + CFINI( NCELL, ISPOLD ) = CNEW( NCELL,ISP ) +#endif IF( CONVERT_CONC( ISPOLD ) )THEN CONC( COL,ROW,LEV,SPC ) = REAL( REVERSE_CONV( ISPOLD ) & * BLKDENS( NCELL ) * CNEW( NCELL,ISP ), 4) @@ -737,8 +769,39 @@ END SUBROUTINE HETCHEM_UPDATE_AERO IF ( LIRRBLK ) CALL PA_IRR_BLKENDC ( OFFSET, CCOL, CROW, CLEV, & NORDCELL, NIRRCLS, IRRCELL ) - ENDIF +#ifdef sens + + DO NCELL = 1, NUMCELLS + DO IRXN = 1, NRXNS + SRK( IRXN ) = RK( NCELL,IRXN ) + IF ( HIGH ) THEN + IF ( NREACT( IRXN ) .LE. 1 ) THEN + SRK2( IRXN ) = 0.0 + ELSE IF ( NREACT( IRXN ) .EQ. 2 ) THEN + SRK2( IRXN ) = REAL( RK( NCELL,IRXN ), 4 ) + ELSE IF ( NREACT( IRXN ) .EQ. 3 ) THEN + SRK2( IRXN ) = REAL( RK( NCELL,IRXN ),4 ) + ELSE + SRK2( IRXN ) = 0.0 + END IF + END IF + END DO + + DO ISP = 1,NUMB_MECH_SPC +! ISPOLD = INEW2OLD( ISP,NCS ) +! YCDDM( ISP ) = 0.5D0 * ( CINIT(NCELL,ISP)+CFINI( NCELL,ISP ) ) + YCDDM( ISP ) = CAVEG(NCELL,ISP) + END DO + CELLNUM = NORDCELL( OFFSET + NCELL ) + COL = CCOL( CELLNUM ) + ROW = CROW( CELLNUM ) + LEV = CLEV( CELLNUM ) + CALL SOLVE_DDM3D_CHEM( COL,ROW,LEV,CHEMSTEP ) +! CALL SOLVE_DDM3D_CHEM( COL,ROW,LEV,CHEMSTEP,IOLD2NEW,INEW2OLD ) + END DO +#endif + ENDIF 500 CONTINUE ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc diff --git a/CCTM/src/gas/smvgear/grsmvgear.F b/CCTM/src/gas/smvgear/grsmvgear.F index 8a2acdcf7..754489251 100644 --- a/CCTM/src/gas/smvgear/grsmvgear.F +++ b/CCTM/src/gas/smvgear/grsmvgear.F @@ -257,6 +257,10 @@ SUBROUTINE SMVGEAR( IRUN, JDATE, JTIME, CHSTEP, LIRRFLAG, END IF +#ifdef sens + CAVEG = 0.0D0 +#endif + cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Set up for debug output and initialize variables cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc @@ -443,14 +447,20 @@ SUBROUTINE SMVGEAR( IRUN, JDATE, JTIME, CHSTEP, LIRRFLAG, IF( CALL_DEG ) THEN !:WTH Initialize for degradation routines YNEW_DEGRADE = 0.0D0 DO I = 1, ISCHAN - ISPOLD = INEW2OLD( I,NCS ) - SPC = CGRID_INDEX( ISPOLD ) + SPC = CGRID_INDEX( I ) DO NCELL = 1, NUMCELLS YNEW_DEGRADE( NCELL,SPC ) = BLKCONC( NCELL, I ) END DO END DO ENDIF +#ifdef sens + DO ISPOLD = 1, ISCHAN + DO NCELL = 1, NUMCELLS + CFINI( NCELL,ISPOLD ) = BLKCONC( NCELL, ISPOLD ) + END DO + END DO +#endif cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Update coefficients of the order. NQQ is the order. ASET and c PERTST are defined in SUBROUTINE JSPARSE. Note that PERTST @@ -987,6 +997,18 @@ SUBROUTINE SMVGEAR( IRUN, JDATE, JTIME, CHSTEP, LIRRFLAG, ENDIF ENDIF +#ifdef sens + DO ISPOLD = 1, ISCHAN + ISPNEW = IOLD2NEW( ISPOLD,1 ) + DO NCELL = 1, NUMCELLS + CAVEG( NCELL, ISPOLD ) = CAVEG( NCELL, ISPOLD ) + & + 0.5D0 * ( CFINI( NCELL, ISPOLD ) + CONC( NCELL, ISPNEW ) ) + & * DELT + CFINI( NCELL, ISPOLD ) = CONC( NCELL, ISPNEW ) + ENDDO + ENDDO +#endif + IF( LCONCOUT ) THEN WRITE( IUNCOUT ) IRUN, BLKID, CELLOUT, CCOLOUT, CROWOUT, & CLEVOUT, NSTEPS, ( RUNMIN + XELAPS ), @@ -1169,6 +1191,15 @@ SUBROUTINE SMVGEAR( IRUN, JDATE, JTIME, CHSTEP, LIRRFLAG, & CNEW( NCELL, ISPC ) = CONCMIN 1180 CONTINUE +#ifdef sens + CHEMSTEP = 1.0D0 / CHEMSTEP + DO I = 1, ISCHAN + DO NCELL = 1, NUMCELLS + CAVEG( NCELL,I ) = MAX( CAVEG( NCELL,I )*CHEMSTEP,CONCMIN ) + ENDDO + ENDDO +#endif + IF( LTRCOUT ) THEN WRITE( IUNTRC, 93200 ) IRUN, N999, NSUBFUN, NPDERIV, NUMNEWT, & NCFAIL, NEFAIL, NZERO, XELAPS, ZERO,