Merge pull request #97 from TerribleNews/adjoint

HEMCO updates for GCHP adjoint from @TerribleNews
geoschem · Oct 4, 2021 · 792f28a · 792f28a
2 parents 2e5b703 + d288de2
commit 792f28a
Showing 1 changed file with 70 additions and 5 deletions.
diff --git a/src/Core/hco_calc_mod.F90 b/src/Core/hco_calc_mod.F90
@@ -895,6 +895,8 @@ SUBROUTINE Get_Current_Emissions( HcoState, BaseDct,   nI,   nJ,           &
           RC = HCO_FAIL
           RETURN
        ENDIF
+    ELSE
+       LevDct1_Unit = -1
     ENDIF
 
     ! Get the units of LevDct2 (if it exists)
@@ -905,6 +907,8 @@ SUBROUTINE Get_Current_Emissions( HcoState, BaseDct,   nI,   nJ,           &
           CALL HCO_ERROR ( MSG, RC, THISLOC=LOC )
           RETURN
        ENDIF
+    ELSE
+       LevDct2_Unit = -1
     ENDIF
 
     ! Throw an error if boxheight is missing and the units are in meters
@@ -2786,6 +2790,13 @@ SUBROUTINE Get_Current_Emissions_Adj( HcoState, BaseDct,   &
     CHARACTER(LEN=255)      :: MSG, LOC
     LOGICAL                 :: NegScalExist
     LOGICAL                 :: MaskFractions
+    LOGICAL                 :: isLevDct1
+    LOGICAL                 :: isLevDct2
+    LOGICAL                 :: isMaskDct
+    LOGICAL                 :: isPblHt
+    LOGICAL                 :: isBoxHt
+    INTEGER                 :: LevDct1_Unit
+    INTEGER                 :: LevDct2_Unit
 
     ! testing only
     INTEGER                 :: IX, IY
@@ -2797,9 +2808,10 @@ SUBROUTINE Get_Current_Emissions_Adj( HcoState, BaseDct,   &
     ! Initialize
     ScalDct => NULL()
     MaskDct => NULL()
+    LOC     = 'GET_CURRENT_EMISSIONS_ADJ (hco_calc_mod.F90)'
 
     ! Enter
-    CALL HCO_ENTER(HcoState%Config%Err,'GET_CURRENT_EMISSIONS', RC )
+    CALL HCO_ENTER(HcoState%Config%Err,'GET_CURRENT_EMISSIONS_ADJ', RC )
     IF(RC /= HCO_SUCCESS) RETURN
 
     ! testing only:
@@ -2848,6 +2860,55 @@ SUBROUTINE Get_Current_Emissions_Adj( HcoState, BaseDct,   &
        LevDct2 => NULL()
     ENDIF
 
+    ! Test whether LevDct1 and LevDct2 are associated
+    isLevDct1 = ASSOCIATED( LevDct1 )
+    isLevDct2 = ASSOCIATED( LevDct2 )
+
+    ! Get the units of LevDct1 (if it exists)
+    IF ( isLevDct1 ) THEN
+       LevDct1_Unit = GetEmisLUnit( HcoState, LevDct1 )
+       IF ( LevDct1_Unit < 0 ) THEN
+          MSG = 'LevDct1 units are not defined!'
+          CALL HCO_ERROR ( HcoState%Config%Err, MSG, RC, THISLOC=LOC )
+          RC = HCO_FAIL
+          RETURN
+       ENDIF
+    ELSE
+       LevDct1_Unit = -1
+    ENDIF
+
+    ! Get the units of LevDct2 (if it exists)
+    IF ( isLevDct2 ) THEN
+       LevDct2_Unit = GetEmisLUnit( HcoState, LevDct2 )
+       IF ( LevDct2_Unit < 0 ) THEN
+          MSG = 'LevDct2_Units are not defined!'
+          CALL HCO_ERROR ( HcoState%Config%Err, MSG, RC, THISLOC=LOC )
+          RETURN
+       ENDIF
+    ELSE
+       LevDct2_Unit = -1
+    ENDIF
+
+    ! Throw an error if boxheight is missing and the units are in meters
+    IF ( LevDct1_Unit == HCO_EMISL_M  .or.                                  &
+         LevDct2_Unit == HCO_EMISL_M ) THEN
+       IF ( .NOT. ASSOCIATED(HcoState%Grid%BXHEIGHT_M%Val) ) THEN
+          MSG = 'Boxheight (in meters) is missing in HEMCO state'
+          CALL HCO_ERROR ( HcoState%Config%Err, MSG, RC, THISLOC=LOC )
+          RETURN
+       ENDIF
+    ENDIF
+
+    ! Throw an error if boxheight is missing and the units are in PBL frac
+    IF ( LevDct1_Unit == HCO_EMISL_PBL  .or.                                &
+         LevDct2_Unit == HCO_EMISL_PBL ) THEN
+       IF ( .NOT. ASSOCIATED(HcoState%Grid%PBLHEIGHT%Val) ) THEN
+          MSG = 'Boundary layer height is missing in HEMCO state'
+          CALL HCO_ERROR ( HcoState%Config%Err, MSG, RC, THISLOC=LOC )
+          RETURN
+       ENDIF
+    ENDIF
+
     ! Loop over all latitudes and longitudes
 !$OMP PARALLEL DO                                                      &
 !$OMP DEFAULT( SHARED )                                                &
@@ -2866,8 +2927,10 @@ SUBROUTINE Get_Current_Emissions_Adj( HcoState, BaseDct,   &
        ENDIF
 
        ! Get lower and upper vertical index
-       CALL GetVertIndx ( HcoState, BaseDct, LevDct1, LevDct2, &
-                          I, J, LowLL, UppLL, RC )
+       CALL GetVertIndx ( HcoState,     BaseDct,   isLevDct1, LevDct1,       &
+                          LevDct1_Unit, isLevDct2, LevDct2,   LevDct2_Unit,  &
+                          I,            J,         LowLL,     UppLL,         &
+                          RC                                                )
        IF ( RC /= HCO_SUCCESS ) THEN
           WRITE(MSG,*) 'Error getting vertical index at location ',I,J,&
                        ': ', TRIM(BaseDct%cName)
@@ -3091,8 +3154,10 @@ SUBROUTINE Get_Current_Emissions_Adj( HcoState, BaseDct,   &
           ! ------------------------------------------------------------ 
 
           ! Get lower and upper vertical index
-          CALL GetVertIndx ( HcoState, BaseDct, &
-                             LevDct1, LevDct2, I, J, LowLL, UppLL, RC )
+          CALL GetVertIndx( HcoState, BaseDct,       isLevDct1,              &
+                            LevDct1,  LevDct1_Unit,  isLevDct2,              &
+                            LevDct2,  LevDct2_Unit,  I,                      &
+                            J,        LowLL,         UppLL,      RC         )
           IF ( RC /= HCO_SUCCESS ) THEN
              ERROR = 1 ! Will cause error
              EXIT