Merge pull request #183 from sergeynk/inst_pvo3

Enable DDM-3D and ISAM calculations and output for Potential Vorticity option

CCTM/src/ddm3d/opasens.F

@@ -32,6 +32,7 @@ SUBROUTINE OPASENS ( JDATE, JTIME, TSTEP )
       USE GRID_CONF             ! horizontal & vertical domain specifications
       USE AVG_CONC              ! integral average CONC
       USE UTILIO_DEFN
+      USE RUNTIME_VARS
 #ifdef parallel
       USE SE_MODULES            ! stenex (using SE_UTIL_MODULE)
 #else
@@ -73,6 +74,17 @@ SUBROUTINE OPASENS ( JDATE, JTIME, TSTEP )
       INTEGER VGTYP_RF
       REAL VGTOP_RF
 !-----------------------------------------------------------------------
+! Define N_ASENS_VARS when a subset of species is requested and it is
+! not defined elsewhere
+
+      IF ( N_ASENS_VARS .EQ. 0 ) THEN ! this happens when 'ALL' keyword
+                                      !is not specified in the run script
+         IF ( W_VEL ) THEN
+            N_ASENS_VARS = N_ACONC_VARS - 1
+         ELSE
+            N_ASENS_VARS = N_ACONC_VARS
+         END IF
+      END IF
 
 ! Change output date/time to starting date/time - e.g. timestamp 1995196:090000
 ! represents data computed from time 1995196:090000 to 1995196:100000

CCTM/src/isam/SA_DEFN.F

@@ -1087,6 +1087,7 @@ SUBROUTINE GET_NSPC_SA ()
         ! Add Ozone if Requested
         IF ( INDEX( TAGCLASSES,'OZONE' ) .NE. 0 .OR. L_OZONE ) THEN
 
+          L_OZONE   = .TRUE.      
           L_NO3     = .TRUE.
           L_VOC     = .TRUE.
           FOUND_SPECIES = .FALSE.

CCTM/src/isam/op_sa.F

@@ -185,8 +185,10 @@ SUBROUTINE OP_SA ( JDATE, JTIME, TSTEP, NSTEPS, FILE_NAME )
         VGTOP3D = VGTOP_GD
 
 
-        DO L = ISAM_BLEV, ISAM_ELEV !1, NLAYS3D + 1
-            VGLVS3D( L ) = VGLVS_GD( L )
+        K = 0
+        DO L = ISAM_BLEV, ISAM_ELEV
+           K = 1 + K
+           VGLVS3D( K ) = VGLVS_GD( L )
         END DO
 
 !       GDNAM3D = GDNAME_GD
@@ -217,8 +219,10 @@ SUBROUTINE OP_SA ( JDATE, JTIME, TSTEP, NSTEPS, FILE_NAME )
 C
          ! print *, S, N_SPCTAG,ISAM_BLEV, ISAM_ELEV
           DO S = 1, N_SPCTAG
+            K = 0
             DO L = ISAM_BLEV, ISAM_ELEV
-              ISAM_BUFF( :,:,L,S ) = ISAM( :,:,L,S_SPCTAG(S),T_SPCTAG(S) )
+              K = K + 1
+              ISAM_BUFF( :,:,K,S ) = ISAM( :,:,L,S_SPCTAG(S),T_SPCTAG(S) )
             ENDDO
           END DO ! loop over total tags
           IF ( .NOT. WRITE3( FILE_NAME, ALLVAR3, JDATE, JTIME, 

CCTM/src/isam/wr_avg_sa.F

@@ -309,7 +309,7 @@ SUBROUTINE WR_AVG_SA ( JDATE, JTIME, TSTEP )
          L = 0
          DO K = AISAM_BLEV, AISAM_ELEV
           L = L + 1
-          BUF4( :,:,L,SPC ) = AISAM( :,:,K, S_SPCTAG(SPC), T_SPCTAG(SPC) )
+          BUF4( :,:,L,SPC ) = AISAM( :,:,L, S_SPCTAG(SPC), T_SPCTAG(SPC) )
          ENDDO !
          IF ( .NOT. WRITE3( SA_ACONC_1, VNAM_SPCTAG( SPC ),
      &      MDATE, MTIME, BUF4( :,:,:,SPC ) ) ) THEN

CCTM/src/isam/wr_sa.F

@@ -130,10 +130,12 @@ SUBROUTINE WR_SA ( JDATE, JTIME, TSTEP, NSTEPS )
       ISAM_BUFF = 0.0
 !Ckrt  print*,'S,VNAME3D(S),S_SPCTAG(S),T_SPCTAG(S)'
       DO S = 1, N_SPCTAG
+       K = 0
        DO L = ISAM_BLEV, ISAM_ELEV
+        K = K + 1
         DO R = 1, NROWS
          DO C = 1, NCOLS
-          ISAM_BUFF( C,R,L,S ) = ISAM( C,R,L,S_SPCTAG(S),T_SPCTAG(S) )
+          ISAM_BUFF( C,R,K,S ) = ISAM( C,R,L,S_SPCTAG(S),T_SPCTAG(S) )
          ENDDO
         ENDDO
        ENDDO

CCTM/src/pv_o3/pvo3.F

@@ -42,7 +42,11 @@ SUBROUTINE PVO3( CGRID, JDATE, JTIME )
       USE SA_DEFN, ONLY: ISAM, N_SPCTAG, NTAG_SA, ITAG, OTHRTAG,
      &                   ISAM_PVO3_MAP, ISAM_SPEC, L_OZONE,
      &                   ISAMRGN_NUM, ISAMRGN_MAP, NSPC_SA
-      USE em_param_module, ONLY: EM_REG_FAC, N_EM_RGN
+      USE desid_param_module, ONLY: DESID_REG_FAC, DESID_N_REG
+#endif
+
+#ifdef sens
+      USE DDM3D_DEFN, ONLY : NP, NPMAX, SENGRID, IPT
 #endif
 
       IMPLICIT NONE
@@ -204,25 +208,25 @@ SUBROUTINE PVO3( CGRID, JDATE, JTIME )
                  ISAM( C,R,L,S_O3,: ) = 1E-30
                  PVO3_STORE = 1.0
                  IF ( SUM( ISAM_PVO3_MAP(:) ) .EQ. 0 ) THEN  ! no PV tracking specified, toss it into OTHER
-                   ISAM( C,R,L,S_O3,OTHRTAG ) = PVO3_STORE
+                   ISAM( C,R,L,S_O3,OTHRTAG ) = CGRID( C,R,L,VO3 )
                  ELSE IF ( SUM( ISAM_PVO3_MAP(:) ) .GE. 1 ) THEN ! at least 1 PVO3 tag
                    DO ITAG = 1, NTAG_SA-3
                      IF( ISAM_PVO3_MAP(ITAG) .EQ. 1 ) THEN ! PV for this tag
                         IF ( ISAMRGN_NUM( ITAG ) .LT. 1 ) THEN ! full domain
-                          ISAM( C,R,L,S_O3,ITAG ) = PVO3_STORE
-                          PVO3_STORE = PVO3_STORE - ISAM( C,R,L,S_O3,ITAG )
+                          ISAM( C,R,L,S_O3,ITAG ) = CGRID( C,R,L,VO3 )
+                          PVO3_STORE = 0.0
                         ELSE ! split into regions and other
                           DO RGN = 1, ISAMRGN_NUM( ITAG )
                             ISAM( C,R,L,S_O3,ITAG ) = ISAM( C,R,L,S_O3,ITAG )
      &                            + CGRID( C,R,L,VO3 )
-     &                            * EM_REG_FAC(C,R,ISAMRGN_MAP(ITAG,RGN))
+     &                            * DESID_REG_FAC(C,R,ISAMRGN_MAP(ITAG,RGN))
                             PVO3_STORE = PVO3_STORE
-     &                            - EM_REG_FAC(C,R,ISAMRGN_MAP(ITAG,RGN))
+     &                            - DESID_REG_FAC(C,R,ISAMRGN_MAP(ITAG,RGN))
                           END DO
                         END IF 
                      END IF
                    END DO
-                   IF ( PVO3_STORE .LT. 1E-6 ) PVO3_STORE = 0.0 ! round off precision error
+                   IF ( PVO3_STORE .LT. 1E-6 .AND. PVO3_STORE .GT. 0.0 ) PVO3_STORE = 0.0 ! round off precision error
                    IF ( PVO3_STORE .GE. 0.0 ) THEN ! put the rest into OTHER tag
                      ISAM( C,R,L,S_O3,OTHRTAG ) = PVO3_STORE * CGRID( C,R,L,VO3 )
                    ELSE
@@ -236,6 +240,15 @@ SUBROUTINE PVO3( CGRID, JDATE, JTIME )
                  END IF
                END IF
 #endif
+
+#ifdef sens
+               DO NP = 1, NPMAX
+                  IF ( IPT( NP ) .EQ. 6 ) THEN
+                     SENGRID( C, R, L, NP, VO3 ) = CGRID( C, R, L, VO3 )
+                  END IF
+               END DO
+#endif
+
                IF ( VO3T .GT. 0 ) THEN
                   CGRID( C,R,L,VO3T ) = CGRID( C,R,L,VO3 )
                END IF

DOCS/Users_Guide/CMAQ_UG_ch10_HDDM-3D.md

@@ -111,8 +111,8 @@ Several sensitivities can be calculated in one simulation. In the example below,
 
     2NX
      HIGH
-     EMISNOX
-     EMISNOX
+     ENX
+     ENX
 
     RT1
      RATE
@@ -121,6 +121,14 @@ Several sensitivities can be calculated in one simulation. In the example below,
 
     END
 
+Example 5
+It is possible to calculate the sensitivity to ozone incursions at the top of the simulated volume if the base model is compiled with potential vorticity module enabled.
+
+    PO3
+     PVO3
+     SPECIES
+      O3  
+
 CMAQ-DDM-3D is flexible in the number of files that the code can handle and also allows for inline emissions streams as well.  Depending on the application and model settings, the following inline streams may be available for sensitivity calculation:
 
 |Stream Keyword | Description|
@@ -132,6 +140,7 @@ CMAQ-DDM-3D is flexible in the number of files that the code can handle and also
 |ASEA|Sea Spray Aerosol Emissions|
 |DUST|Wind-Blown Dust Emissions|
 
+
 ## 10.3.2.1 DDM-3D Control File Format
 
 For each sensitivity:

DOCS/Users_Guide/CMAQ_UG_ch11_ISAM.md

@@ -48,6 +48,10 @@ This version is available as a zip file from the following address:
 
 https://www.cmascenter.org/ioapi/download/ioapi-3.2-large-2020220.tar.gz
 
+**A note about solver configuration**
+
+ISAM is currently only implemented for the EBI solver. For some CMAQ applications, the Rosenbrock solver is desirable and is set as the default in the sample runscripts. For example, this is the case when using the CB6R5M chemical mechanism for hemispheric simulations. In such cases, the bldit script needs to be modified to not select the Rosenbrock solver. While this may incur a performance penalty in terms of CPU time and increase the likelyhood of convergence warnings, it will allow the ISAM simulation to proceed in most cases.
+
 ## 11.3 Run Instructions
 
 To begin a CMAQ simulation with source apportionment enabled, the ISAM section of the runscript must be configured.  The additional necessary environment variables are listed in Table 11-1.
@@ -202,8 +206,10 @@ The CMAQ model allows several types of emissions that are calculated in-line or
 | MGEGM | Marine Gas Emissions |
 | LTNG | Lightning NO Emissions |
 | ASEA | Sea Spray Aerosol Emissions |
-| DUST | Wind-Blown Dust Emissions |   
+| DUST | Wind-Blown Dust Emissions |
+| PVO3 | Potential Vorticity Incursion* |
 
+*Although it is not an emission stream, it is possible to tag the ozone incursions at the top of the simulated volume if the base model is compiled with potential vorticity module enabled.
 
 #### Interpretation of 'OTH' tag
 The OTH tag (e.g.“O3_OTH” in the ISAM benchmark) represents concentrations for that species attributed to 1) all other emissions streams, 2) precursor species not included in the specified tag class(es), and 3) other processes in the model.