From 057363277fc0820ee7d4828ccbd23c4308b9992b Mon Sep 17 00:00:00 2001 From: Colin Lee Date: Fri, 4 Dec 2020 11:59:09 -0800 Subject: [PATCH 1/4] Forward and reverse models compile and run --- src/Core/hco_calc_mod.F90 | 582 +++++++++++++++++- src/Core/hco_diagn_mod.F90 | 145 +++++ src/Core/hco_state_mod.F90 | 7 + src/Core/hco_types_mod.F90 | 3 + src/Core/hcoio_diagn_mod.F90 | 32 + src/Core/hcoio_write_esmf_mod.F90 | 8 + .../Shared/hco_interface_common.F90 | 509 +++++++++++++++ 7 files changed, 1283 insertions(+), 3 deletions(-) diff --git a/src/Core/hco_calc_mod.F90 b/src/Core/hco_calc_mod.F90 index 5539d8ea..f07e7877 100644 --- a/src/Core/hco_calc_mod.F90 +++ b/src/Core/hco_calc_mod.F90 @@ -79,6 +79,9 @@ MODULE HCO_Calc_Mod PUBLIC :: HCO_CheckDepv PUBLIC :: HCO_EvalFld PUBLIC :: HCO_MaskFld +#ifdef ADJOINT + PUBLIC :: GET_CURRENT_EMISSIONS_ADJ +#endif ! ! !PRIVATE MEMBER FUNCTIONS: ! @@ -420,8 +423,17 @@ SUBROUTINE HCO_CalcEmis( HcoState, UseConc, RC ) IF ( Diagn_AutoFillLevelDefined(HcoState%Diagn,3) .AND. DoDiagn ) THEN CALL Diagn_Update( HcoState, ExtNr=ExtNr, & Cat=PrevCat, Hier=-1, HcoID=PrevSpc, & - AutoFill=1, Array3D=CatFlx, COL=-1, RC=RC ) + AutoFill=1, Array3D=CatFlx, COL=HcoState%Diagn%HcoDiagnIDDefault, RC=RC ) IF ( RC /= HCO_SUCCESS ) RETURN +#ifdef ADJOINT + IF (HcoState%IsAdjoint) THEN + CALL Diagn_Update( HcoState, ExtNr=ExtNr, & + Cat=PrevCat, Hier=-1, HcoID=PrevSpc, & + AutoFill=1, Array3D=CatFlx, & + COL=HcoState%Diagn%HcoDiagnIDAdjoint, RC=RC ) + IF ( RC /= HCO_SUCCESS ) RETURN + ENDIF +#endif ENDIF ! Reset CatFlx array and the previously used hierarchy @@ -462,8 +474,17 @@ SUBROUTINE HCO_CalcEmis( HcoState, UseConc, RC ) IF ( Diagn_AutoFillLevelDefined(HcoState%Diagn,2) .AND. DoDiagn ) THEN CALL Diagn_Update( HcoState, ExtNr=ExtNr, & Cat=-1, Hier=-1, HcoID=PrevSpc, & - AutoFill=1,Array3D=SpcFlx, COL=-1, RC=RC ) + AutoFill=1,Array3D=SpcFlx, COL=HcoState%Diagn%HcoDiagnIDDefault, RC=RC ) IF ( RC /= HCO_SUCCESS ) RETURN +#ifdef ADJOINT + IF (HcoState%IsAdjoint) THEN + CALL Diagn_Update( HcoState, ExtNr=ExtNr, & + Cat=-1, Hier=-1, HcoID=PrevSpc, & + AutoFill=1,Array3D=SpcFlx, & + COL=HcoState%Diagn%HcoDiagnIDAdjoint, RC=RC ) + IF ( RC /= HCO_SUCCESS ) RETURN + ENDIF +#endif ENDIF ! Reset arrays and previous hierarchy. @@ -622,8 +643,22 @@ SUBROUTINE HCO_CalcEmis( HcoState, UseConc, RC ) Cat=ThisCat,Hier=ThisHir, HcoID=ThisSpc, & !AutoFill=1, Array3D=TmpFlx, PosOnly=.TRUE.,& AutoFill=1, Array3D=TmpFlx, & - COL=-1, RC=RC ) + COL=HcoState%Diagn%HcoDiagnIDDefault, RC=RC ) IF ( RC /= HCO_SUCCESS ) RETURN +#ifdef ADJOINT + IF (HcoState%IsAdjoint) THEN + ! I don't know why I chose collection=-1 instead of + ! collection=HcoState%Diagn%HcoDiagnIDAdjoint like in the other + ! parts of the adjoint code here, but it's what worked in the + ! old repo so I'm keeping it for now. May need to change + CALL Diagn_Update( HcoState, ExtNr=ExtNr, & + Cat=ThisCat,Hier=ThisHir, HcoID=ThisSpc, & + AutoFill=1, Array3D=TmpFlx, & + COL=-1, RC=RC ) + IF ( RC /= HCO_SUCCESS ) RETURN + ENDIF + +#endif ENDIF ! Update previously used species, category and hierarchy @@ -2608,5 +2643,546 @@ SUBROUTINE GetDilFact ( HcoState, EmisL1, EmisL1Unit, & RC = HCO_SUCCESS END SUBROUTINE GetDilFact +#ifdef ADJOINT +!BOP +! +! !IROUTINE: Get_Current_Emissions +! +! !DESCRIPTION: Subroutine Get\_Current\_Emissions calculates the current +! emissions for the specified emission container. +! This subroutine is only called by HCO\_CalcEmis and for base emission +! containers, i.e. containers of type 1. +!\\ +!\\ +! !INTERFACE: +! + SUBROUTINE Get_Current_Emissions_Adj( HcoState, BaseDct, & + nI, nJ, nL, OUTARR_3D, MASK, RC, UseLL ) +! +! !USES: +! + USE HCO_State_Mod, ONLY : HCO_State + USE HCO_tIdx_MOD, ONLY : tIDx_GetIndx + USE HCO_FileData_Mod, ONLY : FileData_ArrIsDefined +! +! !INPUT PARAMETERS: +! + INTEGER, INTENT(IN) :: nI ! # of lons + INTEGER, INTENT(IN) :: nJ ! # of lats + INTEGER, INTENT(IN) :: nL ! # of levs +! +! !INPUT/OUTPUT PARAMETERS: +! + + TYPE(HCO_State), POINTER :: HcoState ! HEMCO state object + TYPE(DataCont), POINTER :: BaseDct ! base emission + ! container + REAL(hp), INTENT(INOUT) :: OUTARR_3D(nI,nJ,nL) ! output array + REAL(hp), INTENT(INOUT) :: MASK (nI,nJ,nL) ! mask array + INTEGER, INTENT(INOUT) :: RC +! +! !OUTPUT PARAMETERS: +! + INTEGER, INTENT( OUT), OPTIONAL :: UseLL +! +! !REMARKS: +! This routine uses multiple loops over all grid boxes (base emissions +! and scale factors use separate loops). In an OMP environment, this approach +! seems to be faster than using only one single loop (but repeated calls to +! point to containers, etc.). The alternative approach is used in routine +! Get\_Current\_Emissions\_B at the end of this module and may be employed +! on request. +! +! !REVISION HISTORY: +! 25 Aug 2012 - C. Keller - Initial Version +! 09 Nov 2012 - C. Keller - MASK update. Masks are now treated +! separately so that multiple masks can be +! added. +! 06 Jun 2014 - R. Yantosca - Cosmetic changes in ProTeX headers +! 07 Sep 2014 - C. Keller - Mask update. Now set mask to zero as soon as +! on of the applied masks is zero. +! 03 Dec 2014 - C. Keller - Now calculate time slice index on-the-fly. +! 29 Dec 2014 - C. Keller - Added scale factor masks. +! 02 Mar 2015 - C. Keller - Now check for missing values. Missing values are +! excluded from emission calculation. +! 26 Oct 2016 - R. Yantosca - Don't nullify local ptrs in declaration stmts +! 11 May 2017 - C. Keller - Added universal scaling +!EOP +!------------------------------------------------------------------------------ +!BOC +! +! !LOCAL VARIABLES: +! + ! Pointers + TYPE(DataCont), POINTER :: ScalDct + TYPE(DataCont), POINTER :: MaskDct + TYPE(DataCont), POINTER :: LevDct1 + TYPE(DataCont), POINTER :: LevDct2 + + ! Scalars + REAL(sp) :: TMPVAL, MaskScale + REAL(hp) :: DilFact + REAL(hp) :: ScalFact + INTEGER :: tIDx, IDX + INTEGER :: I, J, L, N + INTEGER :: LowLL, UppLL, ScalLL, TmpLL + INTEGER :: ERROR + INTEGER :: TotLL, nnLL + CHARACTER(LEN=255) :: MSG, LOC + LOGICAL :: NegScalExist + LOGICAL :: MaskFractions + + ! testing only + INTEGER :: IX, IY + + !================================================================= + ! GET_CURRENT_EMISSIONS begins here + !================================================================= + + ! Initialize + ScalDct => NULL() + MaskDct => NULL() + + ! Enter + CALL HCO_ENTER(HcoState%Config%Err,'GET_CURRENT_EMISSIONS', RC ) + IF(RC /= HCO_SUCCESS) RETURN + + ! testing only: + IX = 3 !-1 + IY = 8 !-1 + + ! Check if container contains data + IF ( .NOT. FileData_ArrIsDefined(BaseDct%Dta) ) THEN + MSG = 'Array not defined: ' // TRIM(BaseDct%cName) + CALL HCO_ERROR( HcoState%Config%Err, MSG, RC ) + RETURN + ENDIF + + ! Initialize mask. By default, assume that we use all grid boxes. + MASK(:,:,:) = 1.0_hp + MaskFractions = HcoState%Options%MaskFractions + + ! Verbose + IF ( HCO_IsVerb(HcoState%Config%Err,2) ) THEN + WRITE(MSG,*) 'Evaluate field ', TRIM(BaseDct%cName) + CALL HCO_MSG(HcoState%Config%Err,MSG,SEP1=' ') + ENDIF + + ! ---------------------------------------------------------------- + ! Set base emissions + ! ---------------------------------------------------------------- + + ! Initialize ERROR. Will be set to 1 if error occurs below + ERROR = 0 + + ! Initialize variables to compute average vertical level index + totLL = 0 + nnLL = 0 + + ! Check for level index containers + IF ( BaseDct%levScalID1 > 0 ) THEN + CALL Pnt2DataCont( HcoState, BaseDct%levScalID1, LevDct1, RC ) + IF ( RC /= HCO_SUCCESS ) RETURN + ELSE + LevDct1 => NULL() + ENDIF + IF ( BaseDct%levScalID2 > 0 ) THEN + CALL Pnt2DataCont( HcoState, BaseDct%levScalID2, LevDct2, RC ) + IF ( RC /= HCO_SUCCESS ) RETURN + ELSE + LevDct2 => NULL() + ENDIF + + ! Loop over all latitudes and longitudes +!$OMP PARALLEL DO & +!$OMP DEFAULT( SHARED ) & +!$OMP PRIVATE( I, J, L, tIdx, TMPVAL, DilFact, LowLL, UppLL ) & +!$OMP SCHEDULE( DYNAMIC ) + DO J = 1, nJ + DO I = 1, nI + + ! Get current time index for this container and at this location + tIDx = tIDx_GetIndx( HcoState, BaseDct%Dta, I, J ) + IF ( tIDx < 1 ) THEN + WRITE(MSG,*) 'Cannot get time slice index at location ',I,J,& + ': ', TRIM(BaseDct%cName), tIDx + ERROR = 1 + EXIT + ENDIF + + ! Get lower and upper vertical index + CALL GetVertIndx ( HcoState, BaseDct, LevDct1, LevDct2, & + I, J, LowLL, UppLL, RC ) + IF ( RC /= HCO_SUCCESS ) THEN + WRITE(MSG,*) 'Error getting vertical index at location ',I,J,& + ': ', TRIM(BaseDct%cName) + ERROR = 1 ! Will cause error + EXIT + ENDIF + + ! average upper level + totLL = totLL + UppLL + nnLL = nnLL + 1 + + ! Loop over all levels + DO L = LowLL, UppLL + + ! Get base value. Use uniform value if scalar field. + IF ( BaseDct%Dta%SpaceDim == 1 ) THEN + TMPVAL = BaseDct%Dta%V2(tIDx)%Val(1,1) + ELSEIF ( BaseDct%Dta%SpaceDim == 2 ) THEN + TMPVAL = BaseDct%Dta%V2(tIDx)%Val(I,J) + ELSE + TMPVAL = BaseDct%Dta%V3(tIDx)%Val(I,J,L) + ENDIF + + ! If it's a missing value, mask box as unused and set value to + ! zero +#if defined( ESMF_ ) + ! SDE 2017-01-07: Temporary kludge. MAPL ExtData sets missing + ! data to 1e15, but HEMCO uses a different value! + IF ( ( TMPVAL == HCO_MISSVAL ) .or. ( TMPVAL > 1.0e+14 ) ) THEN +#else + IF ( TMPVAL == HCO_MISSVAL ) THEN +#endif + MASK(I,J,:) = 0.0_hp + OUTARR_3D(I,J,L) = 0.0_hp + + ! Pass base value to output array + ELSE + + ! Get dilution factor. Never dilute 3D emissions. + IF ( BaseDct%Dta%SpaceDim == 3 ) THEN + DilFact = 1.0 + + ! 2D dilution factor + ELSE + CALL GetDilFact ( HcoState, BaseDct%Dta%EmisL1, & + BaseDct%Dta%EmisL1Unit, BaseDct%Dta%EmisL2, & + BaseDct%Dta%EmisL2Unit, I, J, L, LowLL, & + UppLL, DilFact, RC ) + IF ( RC /= HCO_SUCCESS ) THEN + WRITE(MSG,*) 'Error getting dilution factor at ',I,J,& + ': ', TRIM(BaseDct%cName) + ERROR = 1 + EXIT + ENDIF + ENDIF + + ! Scale base emission by dilution factor + OUTARR_3D(I,J,L) = DilFact * TMPVAL + ENDIF + ENDDO !L + + ENDDO !I + ENDDO !J +!$OMP END PARALLEL DO + + ! Check for error + IF ( ERROR == 1 ) THEN + CALL HCO_ERROR( HcoState%Config%Err, MSG, RC ) + RETURN + ENDIF + + ! ---------------------------------------------------------------- + ! Apply scale factors + ! The container IDs of all scale factors associated with this base + ! container are stored in vector Scal_cID. + ! ---------------------------------------------------------------- + + ! Loop over scale factors + IF ( BaseDct%nScalID > 0 ) THEN + + DO N = 1, BaseDct%nScalID + + ! Get the scale factor container ID for the current slot + IDX = BaseDct%Scal_cID(N) + + ! Point to data container with the given container ID + CALL Pnt2DataCont( HcoState, IDX, ScalDct, RC ) + IF ( RC /= HCO_SUCCESS ) RETURN + + ! Sanity check: scale field cannot be a base field + IF ( (ScalDct%DctType == HCO_DCTTYPE_BASE) ) THEN + MSG = 'Wrong scale field type: ' // TRIM(ScalDct%cName) + CALL HCO_ERROR( HcoState%Config%Err, MSG, RC ) + RETURN + ENDIF + + ! Skip this scale factor if no data defined. This is possible + ! if scale factors are only defined for a given time range and + ! the simulation datetime is outside of this range. + IF ( .NOT. FileData_ArrIsDefined(ScalDct%Dta) ) THEN + IF ( HCO_IsVerb(HcoState%Config%Err,2) ) THEN + MSG = 'Skip scale factor '//TRIM(ScalDct%cName)// & + ' because it is not defined for this datetime.' + CALL HCO_MSG(HcoState%Config%Err,MSG) + ENDIF + CYCLE + ENDIF + + ! Verbose mode + IF ( HCO_IsVerb(HcoState%Config%Err,2) ) THEN + MSG = 'Applying scale factor ' // TRIM(ScalDct%cName) + CALL HCO_MSG(HcoState%Config%Err,MSG) + ENDIF + + ! Get vertical extension of this scale factor array. + IF( (ScalDct%Dta%SpaceDim<=2) ) THEN + ScalLL = 1 + ELSE + ScalLL = SIZE(ScalDct%Dta%V3(1)%Val,3) + ENDIF + + ! Check if there is a mask field associated with this scale + ! factor. In this case, get a pointer to the corresponding + ! mask field and evaluate scale factors only inside the mask + ! region. + IF ( ASSOCIATED(ScalDct%Scal_cID) ) THEN + CALL Pnt2DataCont( HcoState, ScalDct%Scal_cID(1), MaskDct, RC ) + IF ( RC /= HCO_SUCCESS ) RETURN + + ! Must be mask field + IF ( MaskDct%DctType /= HCO_DCTTYPE_MASK ) THEN + MSG = 'Invalid mask for scale factor: '//TRIM(ScalDct%cName) + MSG = TRIM(MSG) // '; mask: '//TRIM(MaskDct%cName) + CALL HCO_ERROR( HcoState%Config%Err, MSG, RC ) + RETURN + ENDIF + ENDIF + + ! Reinitialize error flag. Will be set to 1 or 2 if error occurs, + ! and to -1 if negative scale factor is ignored. + ERROR = 0 + + ! Loop over all latitudes and longitudes +!$OMP PARALLEL DO & +!$OMP DEFAULT( SHARED ) & +!$OMP PRIVATE( I, J, tIdx, TMPVAL, L, LowLL, UppLL, tmpLL, MaskScale ) & +!$OMP SCHEDULE( DYNAMIC ) + DO J = 1, nJ + DO I = 1, nI + + ! ------------------------------------------------------------ + ! If there is a mask associated with this scale factors, check + ! if this grid box is within or outside of the mask region. + ! Values that partially fall into the mask region are either + ! treated as binary (100% inside or outside), or partially + ! (using the real grid area fractions), depending on the + ! HEMCO options. + ! ------------------------------------------------------------ + + ! Default mask scaling is 1.0 (no mask applied) + MaskScale = 1.0_sp + + ! If there is a mask applied to this scale factor ... + IF ( ASSOCIATED(MaskDct) ) THEN + CALL GetMaskVal ( MaskDct, I, J, & + MaskScale, MaskFractions, RC ) + IF ( RC /= HCO_SUCCESS ) THEN + ERROR = 4 + EXIT + ENDIF + ENDIF + + ! We can skip this grid box if mask is completely zero + IF ( MaskScale <= 0.0_sp ) CYCLE + + ! Get current time index for this container and at this location + tIDx = tIDx_GetIndx( HcoState, ScalDct%Dta, I, J ) + IF ( tIDx < 1 ) THEN + WRITE(*,*) 'Cannot get time slice index at location ',I,J,& + ': ', TRIM(ScalDct%cName), tIDx + ERROR = 3 + EXIT + ENDIF + + ! Check if this is a mask. If so, add mask values to the MASK + ! array. For now, we assume masks to be binary, i.e. 0 or 1. + ! We may want to change that in future to also support values + ! in between. This is especially important when regridding + ! high resolution masks onto coarser grids! + ! ------------------------------------------------------------ + IF ( ScalDct%DctType == HCO_DCTTYPE_MASK ) THEN + + ! Get mask value + CALL GetMaskVal ( ScalDct, I, J, & + TMPVAL, MaskFractions, RC ) + IF ( RC /= HCO_SUCCESS ) THEN + ERROR = 4 + EXIT + ENDIF + + ! Pass to output mask + MASK(I,J,:) = MASK(I,J,:) * TMPVAL + + ! testing only + IF ( HCO_IsVerb(HcoState%Config%Err,2) .AND. I==1 .AND. J==1 ) THEN + write(MSG,*) 'Mask field ', TRIM(ScalDct%cName), & + ' found and added to temporary mask.' + CALL HCO_MSG(HcoState%Config%Err,MSG) + ENDIF + + ! Advance to next grid box + CYCLE + ENDIF! DctType=MASK + + ! ------------------------------------------------------------ + ! For non-mask fields, apply scale factors to all levels + ! of the base field individually. If the scale factor + ! field has more than one vertical level, use the + ! vertical level closest to the corresponding vertical + ! level of the base emission field + ! ------------------------------------------------------------ + + ! Get lower and upper vertical index + CALL GetVertIndx ( HcoState, BaseDct, & + LevDct1, LevDct2, I, J, LowLL, UppLL, RC ) + IF ( RC /= HCO_SUCCESS ) THEN + ERROR = 1 ! Will cause error + EXIT + ENDIF + + ! Loop over all vertical levels of the base field + DO L = LowLL,UppLL + ! If the vertical level exceeds the number of available + ! scale factor levels, use the highest available level. + IF ( L > ScalLL ) THEN + TmpLL = ScalLL + ! Otherwise use the same vertical level index. + ELSE + TmpLL = L + ENDIF + + ! Get scale factor for this grid box. Use same uniform + ! value if it's a scalar field + IF ( ScalDct%Dta%SpaceDim == 1 ) THEN + TMPVAL = ScalDct%Dta%V2(tidx)%Val(1,1) + ELSEIF ( ScalDct%Dta%SpaceDim == 2 ) THEN + TMPVAL = ScalDct%Dta%V2(tidx)%Val(I,J) + ELSE + TMPVAL = ScalDct%Dta%V3(tidx)%Val(I,J,TmpLL) + ENDIF + + ! Set missing value to one + IF ( TMPVAL == HCO_MISSVAL ) TMPVAL = 1.0_sp + + ! Eventually apply mask scaling + IF ( MaskScale /= 1.0_sp ) THEN + TMPVAL = TMPVAL * MaskScale + ENDIF + + ! For negative scale factor, proceed according to the + ! negative value setting specified in the configuration + ! file (NegFlag = 2: use this value): + IF ( TMPVAL < 0.0_sp .AND. HcoState%Options%NegFlag /= 2 ) THEN + + ! NegFlag = 1: ignore and show warning + IF ( HcoState%Options%NegFlag == 1 ) THEN + ERROR = -1 ! Will prompt warning + CYCLE + + ! Return w/ error otherwise + ELSE + WRITE(*,*) 'Negative scale factor at ',I,J,TmpLL,tidx,& + ': ', TRIM(ScalDct%cName), TMPVAL + ERROR = 1 ! Will cause error + EXIT + ENDIF + ENDIF + + ! ------------------------------------------------------- + ! Apply scale factor in accordance to field operator + ! ------------------------------------------------------- + + ! Oper 1: multiply + IF ( ScalDct%Oper == 1 ) THEN + OUTARR_3D(I,J,L) = OUTARR_3D(I,J,L) * TMPVAL + + ! Oper -1: divide + ELSEIF ( ScalDct%Oper == -1 ) THEN + ! Ignore zeros to avoid NaN + IF ( TMPVAL /= 0.0_sp ) THEN + OUTARR_3D(I,J,L) = OUTARR_3D(I,J,L) / TMPVAL + ENDIF + + ! Oper 2: square + ELSEIF ( ScalDct%Oper == 2 ) THEN + OUTARR_3D(I,J,L) = OUTARR_3D(I,J,L) * TMPVAL * TMPVAL + + ! Return w/ error otherwise (Oper 3 is only allowed for masks!) + ELSE + WRITE(*,*) 'Illegal operator for ', TRIM(ScalDct%cName), ScalDct%Oper + ERROR = 2 ! Will cause error + EXIT + ENDIF + + ENDDO !LL + + ! Verbose mode + if ( HCO_IsVerb(HcoState%Config%Err,3) .and. i == ix .and. j == iy ) then + write(MSG,*) 'Scale field ', TRIM(ScalDct%cName) + CALL HCO_MSG(HcoState%Config%Err,MSG) + write(MSG,*) 'Time slice: ', tIdx + CALL HCO_MSG(HcoState%Config%Err,MSG) + write(MSG,*) 'IX, IY: ', IX, IY + CALL HCO_MSG(HcoState%Config%Err,MSG) + write(MSG,*) 'Scale factor (IX,IY,L1): ', TMPVAL + CALL HCO_MSG(HcoState%Config%Err,MSG) + write(MSG,*) 'Mathematical operation : ', ScalDct%Oper + CALL HCO_MSG(HcoState%Config%Err,MSG) +! write(lun,*) 'Updt (IX,IY,L1): ', OUTARR_3D(IX,IY,1) + endif + + ENDDO !I + ENDDO !J +!$OMP END PARALLEL DO + + ! error check + IF ( ERROR > 0 ) THEN + IF ( ERROR == 1 ) THEN + MSG = 'Negative scale factor found (aborted): ' // TRIM(ScalDct%cName) + ELSEIF ( ERROR == 2 ) THEN + MSG = 'Illegal mathematical operator for scale factor: ' // TRIM(ScalDct%cName) + ELSEIF ( ERROR == 3 ) THEN + MSG = 'Encountered negative time index for scale factor: ' // TRIM(ScalDct%cName) + ELSEIF ( ERROR == 3 ) THEN + MSG = 'Mask error in ' // TRIM(ScalDct%cName) + ELSE + MSG = 'Error when applying scale factor: ' // TRIM(ScalDct%cName) + ENDIF + ScalDct => NULL() + CALL HCO_ERROR( HcoState%Config%Err, MSG, RC ) + RETURN + ENDIF + + ! eventually prompt warning for negative values + IF ( ERROR == -1 ) THEN + MSG = 'Negative scale factor found (ignored): ' // TRIM(ScalDct%cName) + CALL HCO_WARNING( HcoState%Config%Err, MSG, RC ) + ENDIF + + ! Free pointer + MaskDct => NULL() + + ENDDO ! N + ENDIF ! N > 0 + + ! Update optional variables + IF ( PRESENT(UseLL) ) THEN + UseLL = 1 + IF ( nnLL > 0 ) UseLL = NINT(REAL(TotLL,4)/REAL(nnLL,4)) + ENDIF + + ! Weight output emissions by mask + OUTARR_3D = OUTARR_3D * MASK + + ! Cleanup and leave w/ success + ScalDct => NULL() + CALL HCO_LEAVE ( HcoState%Config%Err, RC ) + + END SUBROUTINE Get_Current_Emissions_Adj +!EOC +#endif !EOC END MODULE HCO_Calc_Mod diff --git a/src/Core/hco_diagn_mod.F90 b/src/Core/hco_diagn_mod.F90 index 737eeac7..375eaa13 100644 --- a/src/Core/hco_diagn_mod.F90 +++ b/src/Core/hco_diagn_mod.F90 @@ -113,6 +113,7 @@ MODULE HCO_Diagn_Mod USE HCO_Arr_Mod USE HCO_Clock_Mod USE HCO_State_Mod, ONLY : HCO_State + USE MAPL_CommsMod, ONLY : MAPL_am_I_Root IMPLICIT NONE PRIVATE @@ -449,6 +450,71 @@ SUBROUTINE HcoDiagn_Init( HcoState, RC ) ! Pass this collection ID to fixed variable for easy further ! reference to this collection HcoState%Diagn%HcoDiagnIDRestart = CollectionID +#ifdef ADJOINT + IF ( HcoState%isAdjoint ) THEN + ! ------------------------------------------------------------------ + ! Default diagnostics + ! ------------------------------------------------------------------ + CALL DiagnCollection_GetDefaultDelta ( HcoState, & + deltaYMD, deltaHMS, RC ) + IF ( RC /= HCO_SUCCESS ) RETURN + + ! Try to get prefix from configuration file + CALL GetExtOpt ( HcoState%Config, CoreNr, 'DiagnPrefix', & + OptValChar=DiagnPrefix, FOUND=FOUND, RC=RC ) + IF ( RC /= HCO_SUCCESS ) RETURN + IF ( .NOT. FOUND ) THEN +#if defined( MODEL_GEOS ) + DiagnPrefix = 'HEMCO_Diagnostics.$YYYY$MM$DD$HH$MN.nc' +#else + DiagnPrefix = 'HEMCO_diagnostics' +#endif + ENDIF + + ! Output time stamp location + CALL GetExtOpt ( HcoState%Config, CoreNr, 'DiagnTimeStamp', & + OptValChar=OutTimeStampChar, FOUND=FOUND, RC=RC ) + IF ( RC /= HCO_SUCCESS ) RETURN + IF ( .NOT. FOUND ) THEN + OutTimeStamp = HcoDiagnStart + ELSE + CALL TRANLC( OutTimeStampChar ) + IF ( TRIM(OutTimeStampChar) == 'start' ) THEN + OutTimeStamp = HcoDiagnStart + + ELSEIF ( TRIM(OutTimeStampChar) == 'mid' ) THEN + OutTimeStamp = HcoDiagnMid + + ELSEIF ( TRIM(OutTimeStampChar) == 'end' ) THEN + OutTimeStamp = HcoDiagnEnd + + ELSE + WRITE(MSG,*) 'Unrecognized output time stamp location: ', & + TRIM(OutTimeStampChar), ' - will use default (start)' + CALL HCO_WARNING(HcoState%Config%Err,MSG,RC,THISLOC=LOC,WARNLEV=1) + OutTimeStamp = HcoDiagnStart + ENDIF + ENDIF + + CALL DiagnCollection_Create( HcoState%Diagn, & + NX = HcoState%NX, & + NY = HcoState%NY, & + NZ = HcoState%NZ, & + TS = HcoState%TS_EMIS, & + AM2 = HcoState%Grid%AREA_M2%Val, & + COL = CollectionID, & + PREFIX = TRIM(DiagnPrefix), & + deltaYMD = deltaYMD, & + deltaHMS = deltaHMS, & + OutTimeStamp = OutTimeStamp, & + RC = RC ) + IF ( RC /= HCO_SUCCESS ) RETURN + + ! Pass this collection ID to fixed variable for easy further + ! reference to this collection + HcoState%Diagn%HcoDiagnIDAdjoint = CollectionID + endif +#endif ! ------------------------------------------------------------------ ! Manual diagnostics @@ -600,6 +666,26 @@ SUBROUTINE Diagn_DefineFromConfig( HcoState, RC ) HcoID = HCO_GetHcoID( TRIM(SpcName), HcoState ) IF ( HcoID <= 0 ) CYCLE +#ifdef ADJOINT + if ( cName(1:6) == 'SFEmis' ) then + ! ------------------------------------------------------------------ + ! Add it to the HEMCO diagnostics collection + ! ------------------------------------------------------------------ + CALL Diagn_Create( HcoState, & + cName = cName, & + long_name = lName, & + HcoID = HcoID, & + ExtNr = ExtNr, & + Cat = Cat, & + Hier = Hier, & + SpaceDim = SpaceDim, & + OutUnit = OutUnit, & + OutOper = 'CumulSum', & + AutoFill = 1, & + COL = HcoState%Diagn%HcoDiagnIDAdjoint, & + RC = RC ) + else +#endif ! ------------------------------------------------------------------ ! Add it to the HEMCO diagnostics collection ! ------------------------------------------------------------------ @@ -615,6 +701,9 @@ SUBROUTINE Diagn_DefineFromConfig( HcoState, RC ) AutoFill = 1, & COL = HcoState%Diagn%HcoDiagnIDDefault, & RC = RC ) +#ifdef ADJOINT + endif +#endif IF ( RC /= HCO_SUCCESS ) RETURN ENDDO @@ -1035,6 +1124,12 @@ SUBROUTINE Diagn_Create( HcoState, cName, & ThisDiagn%AreaScal = 1.0_hp / Scal ENDIF + IF (HCO_IsVerb(HcoState%Config%Err,3)) THEN + WRITE(MSG, *) ' ThisDiagn%AreaScal = ', ThisDiagn%AreaScal + CALL HCO_MSG( HcoState%Config%Err, MSG) + WRITE(MSG, *) ' ThisDiagn%MassScal = ', ThisDiagn%MassScal + CALL HCO_MSG( HcoState%Config%Err, MSG) + ENDIF !---------------------------------------------------------------- ! Determine the normalization factors applied to the diagnostics ! before they are written out. Diagnostics are always stored @@ -1861,6 +1956,15 @@ SUBROUTINE Diagn_UpdateDriver( HcoState, cID, cName, & CYCLE ENDIF + IF (HCO_IsVerb(HcoState%Config%Err, 3)) THEN + WRITE(MSG,*) 'ThisDiagn%cName: ', trim(ThisDiagn%cName) + CALL HCO_MSG(HcoState%Config%Err, MSG) + WRITE(MSG,*) 'ThisDiagn%AvgFlag: ', ThisDiagn%AvgFlag + CALL HCO_MSG(HcoState%Config%Err, MSG) + WRITE(MSG,*) 'ThisDiagn%SpaceDim: ', ThisDiagn%SpaceDim + CALL HCO_MSG(HcoState%Config%Err, MSG) + ENDIF + ! Increase counter CNT = CNT + 1 @@ -1909,6 +2013,7 @@ SUBROUTINE Diagn_UpdateDriver( HcoState, cID, cName, & IF ( AS /= 0 ) THEN CALL HCO_ERROR( HcoState%Config%Err,& 'Allocation error Arr3D', RC, THISLOC=LOC ) + if (MAPL_am_I_Root()) WRITE(*,*) 'Allocation error Arr3D_HP' RETURN ENDIF Arr3D = Array3D_HP @@ -1917,6 +2022,7 @@ SUBROUTINE Diagn_UpdateDriver( HcoState, cID, cName, & IF ( AS /= 0 ) THEN CALL HCO_ERROR( HcoState%Config%Err,& 'Allocation error Arr3D', RC, THISLOC=LOC ) + if (MAPL_am_I_Root()) WRITE(*,*) 'Allocation error Arr3D' RETURN ENDIF Arr3D = Array3D @@ -1930,6 +2036,7 @@ SUBROUTINE Diagn_UpdateDriver( HcoState, cID, cName, & IF ( AS /= 0 ) THEN CALL HCO_ERROR( HcoState%Config%Err,& 'Allocation error Arr2D', RC, THISLOC=LOC ) + if (MAPL_am_I_Root()) WRITE(*,*) 'Allocation error Arr2D_HP' RETURN ENDIF Arr2D = Array2D_HP @@ -1938,6 +2045,7 @@ SUBROUTINE Diagn_UpdateDriver( HcoState, cID, cName, & IF ( AS /= 0 ) THEN CALL HCO_ERROR( HcoState%Config%Err,& 'Allocation error Arr2D', RC, THISLOC=LOC ) + if (MAPL_am_I_Root()) WRITE(*,*) 'Allocation error Arr2D' RETURN ENDIF Arr2D = Array2D @@ -2034,6 +2142,7 @@ SUBROUTINE Diagn_UpdateDriver( HcoState, cID, cName, & ENDIF ELSE MSG = 'No array passed for updating ' // TRIM(ThisDiagn%cName) + if (MAPL_am_I_Root()) WRITE(*,*) MSG CALL HCO_ERROR ( HcoState%Config%Err, MSG, RC, THISLOC=LOC ) RETURN ENDIF @@ -2284,6 +2393,7 @@ SUBROUTINE Diagn_Get( HcoState, & ! Get collection number CALL DiagnCollection_DefineID( HcoState%Diagn, PS, RC, COL=COL, & ThisColl=ThisColl, HcoState=HcoState ) + IF ( RC /= HCO_SUCCESS .and. MAPL_am_I_Root()) WRITE(*,*) 'Failed to get collection number ', col IF ( RC /= HCO_SUCCESS ) RETURN ! Set AutoFill flag @@ -2366,6 +2476,7 @@ SUBROUTINE Diagn_Get( HcoState, & ! Before returning container, make sure its data is ready for output. IF ( ASSOCIATED (DgnCont ) ) THEN CALL DiagnCont_PrepareOutput ( HcoState, DgnCont, RC ) + IF ( RC /= HCO_SUCCESS .and. MAPL_am_I_Root() ) WRITE(*,*) 'DiagnCont_PrepareOutput returned RC: ', RC, ' for Col: ', DgnCont%CollectionID IF ( RC /= HCO_SUCCESS ) RETURN FLAG = HCO_SUCCESS @@ -2978,12 +3089,14 @@ SUBROUTINE DiagnCont_PrepareOutput( HcoState, DgnCont, RC ) !----------------------------------------------------------------------- CALL DiagnCollection_Find( HcoState%Diagn, DgnCont%CollectionID, & FOUND, RC, ThisColl=ThisColl ) + IF ( RC /= HCO_SUCCESS .and. MAPL_am_I_Root() ) WRITE(*,*) 'DiagnCollection_Find failed' IF ( RC /= HCO_SUCCESS ) RETURN ! This should never happen IF ( .NOT. FOUND .OR. .NOT. ASSOCIATED(ThisColl) ) THEN WRITE(MSG,*) 'Diagnostics ', TRIM(DgnCont%cName), ' has invalid ', & 'collection ID of ', DgnCont%CollectionID + if (MAPL_am_I_Root()) WRITE(*,*) MSG CALL HCO_ERROR( HcoState%Config%Err, MSG, RC, THISLOC=LOC ) RETURN ENDIF @@ -2997,6 +3110,7 @@ SUBROUTINE DiagnCont_PrepareOutput( HcoState, DgnCont, RC ) IF ( DgnCont%SpaceDim == 2 ) THEN CALL HCO_ArrAssert( DgnCont%Arr2D, ThisColl%NX, & ThisColl%NY, RC ) + IF ( RC /= HCO_SUCCESS .and. MAPL_am_I_Root() ) WRITE(*,*) 'HCO_ArrAssert 2d failed' IF ( RC /= HCO_SUCCESS ) RETURN ! Make sure it's zero @@ -3005,6 +3119,7 @@ SUBROUTINE DiagnCont_PrepareOutput( HcoState, DgnCont, RC ) ELSEIF ( DgnCont%SpaceDim == 3 ) THEN CALL HCO_ArrAssert( DgnCont%Arr3D, ThisColl%NX, & ThisColl%NY, ThisColl%NZ, RC ) + IF ( RC /= HCO_SUCCESS .and. MAPL_am_I_Root() ) WRITE(*,*) 'HCO_ArrAssert 3d failed' IF ( RC /= HCO_SUCCESS ) RETURN ! Make sure it's zero @@ -3051,6 +3166,7 @@ SUBROUTINE DiagnCont_PrepareOutput( HcoState, DgnCont, RC ) ! Get current month and year CALL HcoClock_Get( HcoState%Clock, cYYYY=YYYY, cMM=MM, RC=RC ) + IF ( RC /= HCO_SUCCESS .and. MAPL_am_I_Root() ) WRITE(*,*) 'HcoClock_Get failed' IF ( RC /= HCO_SUCCESS ) RETURN ! Days per year @@ -3084,6 +3200,7 @@ SUBROUTINE DiagnCont_PrepareOutput( HcoState, DgnCont, RC ) ELSE WRITE(MSG,*) 'Illegal time averaging of ', DgnCont%TimeAvg, & ' for diagnostics ', TRIM(DgnCont%cName) + if (MAPL_am_I_Root()) WRITE(*,*) TRIM(MSG) CALL HCO_ERROR( HcoState%Config%Err, MSG, RC, THISLOC=LOC ) RETURN ENDIF @@ -3093,6 +3210,7 @@ SUBROUTINE DiagnCont_PrepareOutput( HcoState, DgnCont, RC ) ! Error trap IF ( norm1 <= 0.0_hp ) THEN MSG = 'Illegal normalization factor: ' // TRIM(DgnCont%cName) + if (MAPL_am_I_Root()) WRITE(*,*) TRIM(MSG) CALL HCO_ERROR( HcoState%Config%Err, MSG, RC, THISLOC=LOC ) RETURN ENDIF @@ -3106,6 +3224,16 @@ SUBROUTINE DiagnCont_PrepareOutput( HcoState, DgnCont, RC ) ! For 3D: IF ( DgnCont%SpaceDim == 3 ) THEN + IF (MAPL_am_I_Root()) THEN + WRITE(*,*) ' ', TRIM(DgnCont%cName), ' 3D', & + ' Scaling: ', totScal + IF ( DgnCont%AreaFlag == 0) THEN + IF ( ASSOCIATED(DgnCont%Arr2D) ) THEN + WRITE(*,*) ' ', TRIM(DgnCont%cName), & + ' Appying area scaling factor ' + ENDIF + ENDIF + ENDIF DO J = 1, ThisColl%NY DO I = 1, ThisColl%NX @@ -3127,6 +3255,16 @@ SUBROUTINE DiagnCont_PrepareOutput( HcoState, DgnCont, RC ) ! For 2D: ELSEIF ( DgnCont%SpaceDim == 2 ) THEN + IF (MAPL_am_I_Root()) THEN + WRITE(*,*) ' ', TRIM(DgnCont%cName), ' 2D', & + ' Scaling: ', totScal + IF ( DgnCont%AreaFlag == 0) THEN + IF ( ASSOCIATED(DgnCont%Arr2D) ) THEN + WRITE(*,*) ' ', TRIM(DgnCont%cName), & + ' Appying area scaling factor ' + ENDIF + ENDIF + ENDIF DO J = 1, ThisColl%NY DO I = 1, ThisColl%NX @@ -3149,6 +3287,10 @@ SUBROUTINE DiagnCont_PrepareOutput( HcoState, DgnCont, RC ) ! For 1D: ELSE + IF (MAPL_am_I_Root()) THEN + WRITE(*,*) ' ', TRIM(DgnCont%cName), ' 1D', & + ' Scaling: ', totScal + ENDIF DgnCont%Scalar = DgnCont%Scalar * totscal ENDIF @@ -4636,6 +4778,9 @@ SUBROUTINE DiagnBundle_Init ( Diagn ) Diagn%HcoDiagnIDDefault = -999 Diagn%HcoDiagnIDRestart = -999 Diagn%HcoDiagnIDManual = -999 +#ifdef ADJOINT + Diagn%HcoDiagnIDAdjoint = -999 +#endif Diagn%nnCollections = 0 ENDIF diff --git a/src/Core/hco_state_mod.F90 b/src/Core/hco_state_mod.F90 index e7c9b730..7773b9e0 100644 --- a/src/Core/hco_state_mod.F90 +++ b/src/Core/hco_state_mod.F90 @@ -119,6 +119,9 @@ MODULE HCO_State_Mod TYPE(ESMF_GridComp), POINTER :: GridComp TYPE(ESMF_State), POINTER :: IMPORT TYPE(ESMF_State), POINTER :: EXPORT +#endif +#ifdef ADJOINT + LOGICAL :: isAdjoint #endif END TYPE HCO_State ! @@ -322,6 +325,10 @@ SUBROUTINE HcoState_Init( HcoState, HcoConfig, nSpecies, RC ) HcoState%TS_CHEM = 0.0_sp HcoState%TS_DYN = 0.0_sp +#ifdef ADJOINT + HcoState%isAdjoint = .false. +#endif + ! Nullify temporary array. This array may be used as temporary ! place to write emissions into. HcoState%Buffer3D => NULL() diff --git a/src/Core/hco_types_mod.F90 b/src/Core/hco_types_mod.F90 index a9a4faa1..960d1f8c 100644 --- a/src/Core/hco_types_mod.F90 +++ b/src/Core/hco_types_mod.F90 @@ -534,6 +534,9 @@ MODULE HCO_TYPES_MOD INTEGER :: HcoDiagnIDDefault = -999 INTEGER :: HcoDiagnIDRestart = -999 INTEGER :: HcoDiagnIDManual = -999 +#ifdef ADJOINT + INTEGER :: HcoDiagnIDAdjoint = -999 +#endif INTEGER :: nnCollections = 0 END TYPE DiagnBundle ! diff --git a/src/Core/hcoio_diagn_mod.F90 b/src/Core/hcoio_diagn_mod.F90 index 396afe9c..057956ad 100644 --- a/src/Core/hcoio_diagn_mod.F90 +++ b/src/Core/hcoio_diagn_mod.F90 @@ -23,6 +23,7 @@ MODULE HCOIO_DIAGN_MOD ! !USES: ! USE HCO_ERROR_MOD + USE MAPL_CommsMod, ONLY : MAPL_am_I_Root IMPLICIT NONE PRIVATE @@ -102,6 +103,9 @@ SUBROUTINE HcoDiagn_Write( HcoState, Restart, RC ) ! INTEGER :: I, COL CHARACTER(LEN=255) :: MSG, LOC +#ifdef ADJOINT + INTEGER :: MaxIdx +#endif !================================================================= ! HcoDiagn_Write begins here! @@ -112,6 +116,7 @@ SUBROUTINE HcoDiagn_Write( HcoState, Restart, RC ) ! To write restart (enforced) IF ( RESTART ) THEN + IF (MAPL_am_I_Root()) WRITE(*,*) "Doing HCOIO_DIAGN_WRITEOUT for restart" CALL HCOIO_DIAGN_WRITEOUT ( HcoState, & ForceWrite = .TRUE., & UsePrevTime = .FALSE., & @@ -119,10 +124,12 @@ SUBROUTINE HcoDiagn_Write( HcoState, Restart, RC ) RC = RC ) IF( RC /= HCO_SUCCESS) RETURN + IF (MAPL_am_I_Root()) WRITE(*,*) "HcoClock_SetLast" ! Set last flag for use below CALL HcoClock_SetLast ( HcoState%Clock, .TRUE., RC ) IF( RC /= HCO_SUCCESS) RETURN + IF (MAPL_am_I_Root()) WRITE(*,*) "Doing HCOIO_DIAGN_WRITEOUT default column" CALL HCOIO_DIAGN_WRITEOUT ( HcoState, & ForceWrite = .FALSE., & UsePrevTime = .FALSE., & @@ -130,6 +137,19 @@ SUBROUTINE HcoDiagn_Write( HcoState, Restart, RC ) RC = RC ) IF( RC /= HCO_SUCCESS) RETURN +#ifdef ADJOINT + IF (HcoState%isAdjoint) THEN + if (MAPL_am_I_Root()) WRITE(*,*) "Doing HCOIO_DIAGN_WRITEOUT adjoint" + CALL HCOIO_DIAGN_WRITEOUT ( HcoState, & + ForceWrite = .FALSE., & + UsePrevTime = .FALSE., & + COL = HcoState%Diagn%HcoDiagnIDAdjoint, & + RC = RC ) + IF( RC /= HCO_SUCCESS) RETURN + ENDIF +#endif + + if (MAPL_am_I_Root()) WRITE(*,*) "HcoClock_SetLast" ! Reset IsLast flag. This is to ensure that the last flag is not ! carried over (ckeller, 11/1/16). CALL HcoClock_SetLast ( HcoState%Clock, .FALSE., RC ) @@ -140,7 +160,13 @@ SUBROUTINE HcoDiagn_Write( HcoState, Restart, RC ) ! Loop over all collections that shall be written out. ! HCOIO_DIAGN_WRITEOUT will determine whether it is time to ! write a collection or not. +#ifndef ADJOINT DO I = 1, 3 +#else + MaxIdx = 3 + IF (HcoState%isAdjoint) MaxIdx = 4 + DO I = 1, MaxIdx +#endif ! Define collection ID SELECT CASE ( I ) @@ -150,6 +176,10 @@ SUBROUTINE HcoDiagn_Write( HcoState, Restart, RC ) COL = HcoState%Diagn%HcoDiagnIDRestart CASE ( 3 ) COL = HcoState%Diagn%HcoDiagnIDManual +#ifdef ADJOINT + CASE ( 4 ) + COL = HcoState%Diagn%HcoDiagnIDAdjoint +#endif END SELECT #if !defined ( ESMF_ ) @@ -163,6 +193,7 @@ SUBROUTINE HcoDiagn_Write( HcoState, Restart, RC ) ! HCO_RestartWrite. (ckeller, 10/9/17) IF ( I == 2 ) CYCLE #endif + if (MAPL_am_I_Root()) WRITE(*,*) "Doing HCOIO_DIAGN_WRITEOUT, I = ", I ! Restart file CALL HCOIO_DIAGN_WRITEOUT ( HcoState, & @@ -170,6 +201,7 @@ SUBROUTINE HcoDiagn_Write( HcoState, Restart, RC ) UsePrevTime = .FALSE., & COL = COL, & RC = RC ) + if (MAPL_am_I_Root()) WRITE(*,*) "Back from HCOIO_DIAGN_WRITEOUT, RC = ", RC IF(RC /= HCO_SUCCESS) RETURN ENDDO ENDIF diff --git a/src/Core/hcoio_write_esmf_mod.F90 b/src/Core/hcoio_write_esmf_mod.F90 index a429c529..2865eb12 100644 --- a/src/Core/hcoio_write_esmf_mod.F90 +++ b/src/Core/hcoio_write_esmf_mod.F90 @@ -151,10 +151,18 @@ SUBROUTINE HCOIO_WRITE_ESMF ( HcoState, RC, OnlyIfFirst, COL ) ThisDiagn => NULL() DO WHILE ( .TRUE. ) + IF (MAPL_am_I_Root()) WRITE(*,*) "Getting next diagnostic in list" + ! Get next diagnostics in list. This will return the next ! diagnostics container that contains content to be written ! out on this time step. CALL Diagn_Get ( HcoState, EOI, ThisDiagn, FLAG, RC, COL=PS ) + + IF (MAPL_am_I_Root()) THEN + IF ( RC == HCO_SUCCESS .and. FLAG == HCO_SUCCESS ) WRITE(*,*) "Got it! Name: ", TRIM(ThisDiagn%cName) + IF ( RC /= HCO_SUCCESS) WRITE(*,*) "Fail! RC: ", RC, " Flag: ", FLAG + ENDIF + IF ( RC /= HCO_SUCCESS ) RETURN IF ( FLAG /= HCO_SUCCESS ) EXIT diff --git a/src/Interfaces/Shared/hco_interface_common.F90 b/src/Interfaces/Shared/hco_interface_common.F90 index 38718a9e..0aa77e29 100644 --- a/src/Interfaces/Shared/hco_interface_common.F90 +++ b/src/Interfaces/Shared/hco_interface_common.F90 @@ -28,6 +28,9 @@ MODULE HCO_Interface_Common PUBLIC :: SetHcoTime PUBLIC :: GetHcoVal PUBLIC :: GetHcoDiagn +#ifdef FALSE + PUBLIC :: CALC_EMS_SF_ADJ +#endif ! ! !REMARKS: ! These utilities were mostly migrated from GEOS-Chem HCO_Interface_Mod. @@ -330,5 +333,511 @@ SUBROUTINE GetHcoDiagn ( HcoState, ExtState, DiagnName, & RC = HCO_SUCCESS END SUBROUTINE GetHcoDiagn + +#ifdef FALSE ! the adjoint-relevant code in here has been moved to + ! Hco_CalcEmis in #ifdef ADJOINT blocks so I'm disabling + ! this for now +! +! !DESCRIPTION: Subroutine CALC_EMS_SF_ADJ stores adjoint sensitvities to +! emissions categories +!\\ +!\\ +! !INTERFACE +! + SUBROUTINE CALC_EMS_SF_ADJ(SpcId, DT, Input_Opt, State_Chm, State_Diag, RC) +! +! !USES: +! + USE Input_Opt_Mod, ONLY : OptInput + USE State_Chm_Mod, ONLY : ChmState + USE State_Diag_Mod, ONLY : DgnState + USE HCO_STATE_MOD, ONLY : HCO_State + USE HCO_ARR_MOD, ONLY : HCO_ArrAssert + USE HCO_DATACONT_MOD, ONLY : ListCont_NextCont + USE HCO_FILEDATA_MOD, ONLY : FileData_ArrIsDefined + USE HCO_Calc_Mod, ONLY : GET_CURRENT_EMISSIONS_ADJ + USE HCO_Scale_Mod, ONLY : HCO_ScaleArr + USE HCO_Diagn_Mod + USE HCO_Types_Mod + + +! !INPUT PARAMETERS: +! + INTEGER, INTENT(IN ) :: spcId ! Which species + TYPE(OptInput), INTENT(IN ) :: Input_Opt ! Input opts + REAL(fp), INTENT(IN ) :: DT ! Time step [s] +! +! !INPUT/OUTPUT PARAMETERS: +! + TYPE(ChmState), INTENT(INOUT) :: State_Chm ! Chemistry state + TYPE(DgnState), INTENT(INOUT) :: State_Diag ! Diags State + INTEGER, INTENT(INOUT) :: RC ! Success/Failure +! +! !REVISION HISTORY: +! 18 May 2020 - C. Keller - Initial version +!EOP +!------------------------------------------------------------------------------ +!BOC +! +! !LOCAL VARIABLES: +! + ! Working pointers: list and data container + TYPE(ListCont), POINTER :: Lct + TYPE(DataCont), POINTER :: Dct + + ! Temporary emissions arrays + REAL(hp), TARGET :: SpcFlx( HcoState%NX, & + HcoState%NY, & + HcoState%NZ ) + REAL(hp), TARGET :: CatFlx( HcoState%NX, & + HcoState%NY, & + HcoState%NZ ) + REAL(hp), TARGET :: TmpFlx( HcoState%NX, & + HcoState%NY, & + HcoState%NZ ) + REAL(hp) :: Mask ( HcoState%NX, & + HcoState%NY, & + HcoState%NZ ) + REAL(hp) :: HirFlx( HcoState%NX, & + HcoState%NY, & + HcoState%NZ ) + REAL(hp) :: HirMsk( HcoState%NX, & + HcoState%NY, & + HcoState%NZ ) + + ! Integers + INTEGER :: ThisSpc, PrevSpc ! current and previous species ID + INTEGER :: ThisCat, PrevCat ! current and previous category + INTEGER :: ThisHir, PrevHir ! current and previous hierarchy + INTEGER :: SpcMin, SpcMax ! range of species to be considered + INTEGER :: CatMin, CatMax ! range of categories to be considered + INTEGER :: ExtNr ! Extension Nr to be used + INTEGER :: nI, nJ, nL + INTEGER :: nnSpec, FLAG + + LOGICAL :: Found, DoDiagn, EOL, UpdateCat, UseConc + + ! For error handling & verbose mode + CHARACTER(LEN=255) :: MSG + + ! testing / debugging + integer :: ix,iy + + !================================================================= + ! CALC_EMS_SF_ADJ begins here! + !================================================================= + ! Initialize + Lct => NULL() + Dct => NULL() + + + !----------------------------------------------------------------- + ! Initialize variables + !----------------------------------------------------------------- + + ! Initialize + SpcFlx(:,:,:) = 0.0_hp + CatFlx(:,:,:) = 0.0_hp + HirFlx(:,:,:) = 0.0_hp + HirMsk(:,:,:) = 0.0_hp + PrevSpc = -1 + PrevHir = -1 + PrevCat = -1 + nnSpec = 0 + + ! Pass emission grid dimensions + nI = HcoState%NX + nJ = HcoState%NY + nL = HcoState%NZ + + ! Pass calculation options + SpcMin = HcoState%Options%SpcMin !Lower species ID + SpcMax = HcoState%Options%SpcMax !Upper species ID + CatMin = HcoState%Options%CatMin !Lower emission category + CatMax = HcoState%Options%CatMax !Upper emission category + ExtNr = HcoState%Options%ExtNr !Extension number + DoDiagn = HcoState%Options%AutoFillDiagn !Write AutoFill diagnostics? + + ! I am not sure what do about units here, so set this to false for now + UseConc = .false. + + ! Enter routine + CALL HCO_ENTER (HcoState%Config%Err,'CALC_EMS_SF_ADJ (HCO_INTERFACE_MOD.F90)', RC ) + IF(RC /= HCO_SUCCESS) RETURN + + ! Verbose mode + IF ( HCO_IsVerb(HcoState%Config%Err,2) ) THEN + WRITE (MSG, *) 'Run HEMCO calculation w/ following options:' + CALL HCO_MSG ( HcoState%Config%Err, MSG ) + WRITE (MSG, "(A20,I5)") 'Extension number:', ExtNr + CALL HCO_MSG ( HcoState%Config%Err, MSG ) + WRITE (MSG, "(A20,I5,I5)") 'Tracer range:', SpcMin, SpcMax + CALL HCO_MSG ( HcoState%Config%Err, MSG ) + WRITE (MSG, "(A20,I5,I5)") 'Category range:', CatMin, CatMax + CALL HCO_MSG ( HcoState%Config%Err, MSG ) + WRITE (MSG, *) 'Auto diagnostics: ', DoDiagn + CALL HCO_MSG ( HcoState%Config%Err, MSG ) + ENDIF + + !================================================================= + ! Walk through all containers of EmisList and determine the + ! emissions for all containers that qualify for calculation. + ! The containers in EmisList are sorted by species, category and + ! hierarchy. This enables a straightforward, piece-by-piece + ! assembly of the final emission array (start with lowest + ! hierarchy emissions, then overwrite piece-by-piece with higher + ! hierarchy values). + !================================================================= + + ! Point to the head of the emissions linked list + EOL = .FALSE. ! End of list + Lct => NULL() + CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG ) + + ! Do until end of EmisList (==> loop over all emission containers) + DO + ! Have we reached the end of the list? + IF ( FLAG /= HCO_SUCCESS ) THEN + EOL = .TRUE. + ELSE + EOL = .FALSE. + ENDIF + + ! ------------------------------------------------------------ + ! Select container and update all working variables & arrays. + ! ------------------------------------------------------------ + IF ( .NOT. EOL ) THEN + + ! Dct is the current data container + Dct => Lct%Dct + if ( HCO_IsVerb(HcoState%Config%Err,2) ) THEN + WRITE(MSG,*) ' CALC_EMS_SF_ADJ processing container ', trim(Dct%cName) + CALL HCO_MSG(HcoState%Config%Err, MSG) + endif + + ! Check if this is a base field + IF ( Dct%DctType /= HCO_DCTTYPE_BASE ) THEN + CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG ) + CYCLE + ENDIF + + ! Sanity check: Make sure this container holds data. + ! 'Empty' containers are possible if the simulation time + ! is outside of the specified data time range and time + ! slice cycling is deactivated (CycleFlag > 1). + IF( .NOT. FileData_ArrIsDefined(Lct%Dct%Dta) ) THEN + CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG ) + CYCLE + ENDIF + + ! Check if this is the specified extension number + IF ( Dct%ExtNr /= ExtNr ) THEN + CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG ) + CYCLE + ENDIF + + ! Advance to next container if the species ID is outside + ! the specified species range (SpcMin - SpcMax). Consider + ! all species above SpcMin if SpcMax is negative! + IF( ( Dct%HcoID < SpcMin ) .OR. & + ( (Dct%HcoID > SpcMax) .AND. (SpcMax > 0) ) ) THEN + CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG ) + CYCLE + ENDIF + + ! Advance to next emission field if the emission category of + ! the current container is outside of the specified species + ! range (CatMin - CatMax). Consider all categories above CatMin + ! if CatMax is negative! + IF( ( Dct%Cat < CatMin ) .OR. & + ( (Dct%Cat > CatMax) .AND. (CatMax > 0) ) ) THEN + CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG ) + CYCLE + ENDIF + + ! Check if this container holds data in the desired unit format, + ! i.e. concentration data if UseConc is enabled, emission data + ! otherwise. + IF ( UseConc .NEQV. Dct%Dta%IsConc ) THEN + CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG ) + CYCLE + ENDIF + + ! Check if this is our species + IF ( Dct%HcoID /= SpcId ) THEN + CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG ) + CYCLE + ENDIF + + ! Update working variables + ThisSpc = Dct%HcoID + ThisCat = Dct%Cat + ThisHir = Dct%Hier + + ! If end of list, use dummy values for ThisSpc, ThisCat and ThisHir + ! to make sure that emissions are added to HEMCO in the section + ! below! + ELSE + ThisSpc = -1 + ThisCat = -1 + ThisHir = -1 + ENDIF + + !-------------------------------------------------------------------- + ! Before computing emissions of current data container make sure that + ! emissions of previous container are properly archived. + !-------------------------------------------------------------------- + + ! Add emissions on hierarchy level to the category flux array. Do + ! this only if this is a new species, a new category or a new + ! hierarchy level. + ! Note: no need to add to diagnostics because hierarchy level + ! diagnostics are filled right after computing the emissions of + ! a given data container (towards the end of the DO loop). + IF ( (ThisHir /= PrevHir) .OR. & + (ThisSpc /= PrevSpc) .OR. & + (ThisCat /= PrevCat) ) THEN + + ! Add hierarchy level emissions to category array over the + ! covered regions. + CatFlx = ( (1.0_hp - HirMsk) * CatFlx ) + HirFlx + + ! Reset + HirFlx = 0.0_hp + HirMsk = 0.0_hp + ENDIF + + !-------------------------------------------------------------------- + ! If this is a new species or category, pass the previously collected + ! emissions to the species array. Update diagnostics at category level. + ! Skip this step for first species, i.e. if PrevSpc is still -1. + !-------------------------------------------------------------------- + UpdateCat = .FALSE. + IF ( ThisCat /= PrevCat ) UpdateCat = .TRUE. + IF ( ThisSpc /= PrevSpc ) UpdateCat = .TRUE. + IF ( PrevCat <= 0 .OR. PrevSpc <= 0 ) UpdateCat = .FALSE. + IF ( UpdateCat ) THEN + + ! CatFlx holds the emissions for this category. Pass this to + ! the species array SpcFlx. + SpcFlx(:,:,:) = SpcFlx(:,:,:) + CatFlx(:,:,:) + + IF (Input_Opt%IS_FD_SPOT_THIS_PET) THEN + WRITE(*,*) PrevCat, ' CatFlx(IFD,JFD) = ', & + SUM(CatFlx(Input_Opt%IFD,Input_Opt%JFD,:)) + ENDIF + + ! verbose + IF ( HCO_IsVerb(HcoState%Config%Err,3) ) THEN + WRITE(MSG,*) 'Added category emissions to species array: ' + CALL HCO_MSG(HcoState%Config%Err,MSG) + WRITE(MSG,*) 'Species : ', PrevSpc + CALL HCO_MSG(HcoState%Config%Err,MSG) + WRITE(MSG,*) 'Category : ', PrevCat + CALL HCO_MSG(HcoState%Config%Err,MSG) + WRITE(MSG,*) 'Cat. emissions: ', SUM(CatFlx) + CALL HCO_MSG(HcoState%Config%Err,MSG) + WRITE(MSG,*) 'Spc. emissions: ', SUM(SpcFlx) + CALL HCO_MSG(HcoState%Config%Err,MSG) + ENDIF + CatFlx(:,:,:) = DT * State_Chm%SpeciesAdj(:,:,:, SpcId) * CatFlx(:,:,:) + ! verbose + IF ( HCO_IsVerb(HcoState%Config%Err,3) ) THEN + WRITE(MSG,*) 'Cat. emissions: ', SUM(CatFlx) + CALL HCO_MSG(HcoState%Config%Err,MSG) + ENDIF + IF (Input_Opt%IS_FD_SPOT_THIS_PET) THEN + WRITE(*,*) PrevCat, ' CatFlxAdj(IFD,JFD) = ', & + SUM(CatFlx(Input_Opt%IFD,Input_Opt%JFD,:)) + ENDIF + + ! Add category emissions to diagnostics at category level + ! (only if defined in the diagnostics list). + IF ( Diagn_AutoFillLevelDefined(HcoState%Diagn,3) .AND. DoDiagn ) THEN + CALL Diagn_Update( HcoState, ExtNr=ExtNr, & + Cat=PrevCat, Hier=-1, HcoID=PrevSpc, & + AutoFill=1, Array3D=CatFlx, & + COL=HcoState%Diagn%HcoDiagnIDAdjoint, RC=RC ) + IF ( RC /= HCO_SUCCESS ) RETURN + ENDIF + + ! Reset CatFlx array and the previously used hierarchy + ! ==> Emission hierarchies are only important within the + ! same category, hence always start over at lowest hierarchy + ! when entering a new category. + CatFlx(:,:,:) = 0.0_hp + PrevHir = -1 + ENDIF + + !-------------------------------------------------------------------- + ! If this is a new species, pass previously calculated emissions + ! to the final emissions array in HcoState. + ! Update diagnostics at extension number level. + ! Don't do before first emission calculation, i.e. if PrevSpc + ! is still the initialized value of -1! + !-------------------------------------------------------------------- + IF ( ThisSpc /= PrevSpc .AND. PrevSpc > 0 ) THEN + + IF (Input_Opt%IS_FD_SPOT_THIS_PET) THEN + WRITE(*,*) PrevCat, ' SpcFlx(IFD,JFD) = ', & + SUM(SpcFlx(Input_Opt%IFD,Input_Opt%JFD,:)) + ENDIF + ! testing only + IF ( HCO_IsVerb(HcoState%Config%Err,3) ) THEN + WRITE(MSG,*) 'Added total emissions to output array: ' + CALL HCO_MSG(HcoState%Config%Err,MSG) + WRITE(MSG,*) 'Species: ', PrevSpc + CALL HCO_MSG(HcoState%Config%Err,MSG) + WRITE(MSG,*) 'SpcFlx : ', SUM(SpcFlx) + CALL HCO_MSG(HcoState%Config%Err,MSG) + ENDIF + + ! Add to diagnostics at extension number level. + ! The same diagnostics may be updated multiple times during + ! the same time step, continuously adding emissions to it. + IF ( Diagn_AutoFillLevelDefined(HcoState%Diagn,2) .AND. DoDiagn ) THEN + CALL Diagn_Update(HcoState, ExtNr=ExtNr, & + Cat=-1, Hier=-1, HcoID=PrevSpc, & + AutoFill=1,Array3D=SpcFlx, & + COL=HcoState%Diagn%HcoDiagnIDAdjoint, RC=RC ) + IF ( RC /= HCO_SUCCESS ) RETURN + ENDIF + + ! Reset arrays and previous hierarchy. + SpcFlx(:,:,:) = 0.0_hp + PrevCat = -1 + PrevHir = -1 + ENDIF + + !-------------------------------------------------------------------- + ! Exit DO loop here if end of list + !-------------------------------------------------------------------- + IF ( EOL ) EXIT + + !-------------------------------------------------------------------- + ! Update/archive information on species level if needed + !-------------------------------------------------------------------- + IF ( ThisSpc /= PrevSpc .AND. ThisSpc > 0 ) THEN + + ! Update number of species for which emissions have been + ! calculated. + nnSpec = nnSpec + 1 + + ! verbose mode + IF ( HCO_IsVerb(HcoState%Config%Err,2) ) THEN + WRITE(MSG,*) 'Calculating emissions for species ', & + TRIM(HcoState%Spc(ThisSpc)%SpcName) + CALL HCO_MSG( HcoState%Config%Err, MSG, SEP1='-', SEP2='-' ) + ENDIF + ENDIF + + ! verbose mode + IF ( HCO_IsVerb(HcoState%Config%Err,2) ) THEN + WRITE(MSG,*) 'Dct has ', & + Dct%Dta%SpaceDim, ' dimensions.' + CALL HCO_MSG( HcoState%Config%Err, MSG, SEP1='-', SEP2='-' ) + ENDIF + !-------------------------------------------------------------------- + ! Get current emissions and write into TmpFlx array. The array Mask + ! denotes all valid grid boxes for this inventory. + !-------------------------------------------------------------------- + TmpFlx(:,:,:) = 0.0_hp + CALL GET_CURRENT_EMISSIONS_ADJ( HcoState, Dct, & + nI, nJ, nL, TmpFlx, Mask, RC ) + IF ( RC /= HCO_SUCCESS ) RETURN + + ! Eventually add universal scale factor + CALL HCO_ScaleArr( HcoState, ThisSpc, TmpFlx, RC ) + IF ( RC /= HCO_SUCCESS ) RETURN + + ! Check for negative values according to the corresponding setting + ! in the configuration file: 2 means allow negative values, 1 means + ! set to zero and prompt a warning, else return with error. + IF ( HcoState%Options%NegFlag /= 2 ) THEN + + IF ( ANY(TmpFlx < 0.0_hp) ) THEN + + ! Set to zero and prompt warning + IF ( HcoState%Options%NegFlag == 1 ) THEN + WHERE ( TmpFlx < 0.0_hp ) TmpFlx = 0.0_hp + MSG = 'Negative emissions set to zero: '// TRIM(Dct%cName) + CALL HCO_WARNING( HcoState%Config%Err, MSG, RC ) + + ! Return with error + ELSE + MSG = 'Negative emissions in: '// TRIM(Dct%cName) // '. ' // & + 'To allow negatives, edit settings in the configuration file.' + CALL HCO_ERROR( HcoState%Config%Err, MSG, RC ) + RETURN + ENDIF + ENDIF + ENDIF + + ! ------------------------------------------------------------ + ! Collect all emissions of the same category (and species) on + ! the hierarchy level into array HirFlx. HirMsk contains the + ! combined covered region. That is, if there are two regional + ! inventories with the same hierarchy HirMsk will cover both + ! of these regions. + ! The specified field hierarchies determine whether the + ! temporary emissions are added (if hierarchy is the same + ! as the previously used hierarchy), or if they overwrite the + ! previous values in HirFlx (if hierarchy is higher than the + ! previous hierarchy). + ! ------------------------------------------------------------ + + ! Add emissions to the hierarchy array HirFlx if this hierarchy + ! is the same as previous hierarchy + IF ( ThisHir == PrevHir ) THEN + HirFlx = HirFlx + TmpFlx + HirMsk = HirMsk + Mask + + ! Make sure mask values do not exceed 1.0 + WHERE(HirMsk > 1.0 ) HirMsk = 1.0 + + ! If hierarchy is larger than those of the previously used + ! fields, overwrite HirFlx with new values. + ELSE + + HirFlx = TmpFlx + HirMsk = Mask + + ENDIF + + ! Update diagnostics at hierarchy level. Make sure that only + ! positive values are used. + ! The same diagnostics may be updated multiple times + ! during the same time step, continuously adding + ! emissions to it. + ! Now remove PosOnly flag. TmpFlx is initialized to zero, so it's + ! ok to keep negative values (ckeller, 7/12/15). + IF ( Diagn_AutoFillLevelDefined(HcoState%Diagn,4) .AND. DoDiagn ) THEN + CALL Diagn_Update( HcoState, ExtNr=ExtNr, & + Cat=ThisCat,Hier=ThisHir, HcoID=ThisSpc, & + !AutoFill=1, Array3D=TmpFlx, PosOnly=.TRUE.,& + AutoFill=1, Array3D=TmpFlx, & + COL=-1, RC=RC ) + IF ( RC /= HCO_SUCCESS ) RETURN + ENDIF + + ! Update previously used species, category and hierarchy + PrevSpc = ThisSpc + PrevCat = ThisCat + PrevHir = ThisHir + + ! Advance to next emission container + CALL ListCont_NextCont( HcoState%EmisList, Lct, FLAG ) + + ENDDO ! Loop over EmisList + + ! Make sure internal pointers are nullified + Lct => NULL() + Dct => NULL() + + ! Leave w/ success + CALL HCO_LEAVE ( HcoState%Config%Err, RC ) + + END SUBROUTINE CALC_EMS_SF_ADJ +#endif !EOC END MODULE HCO_Interface_Common From 63532e9d38b7be99afb84eeb45a4a0c9b779c45c Mon Sep 17 00:00:00 2001 From: Colin Lee Date: Sun, 4 Apr 2021 11:05:40 -0700 Subject: [PATCH 2/4] Removed logging that stopped GCClassic from compiling --- src/Core/hco_diagn_mod.F90 | 39 ------------------------------------ src/Core/hcoio_diagn_mod.F90 | 8 -------- 2 files changed, 47 deletions(-) diff --git a/src/Core/hco_diagn_mod.F90 b/src/Core/hco_diagn_mod.F90 index 375eaa13..0652ee3f 100644 --- a/src/Core/hco_diagn_mod.F90 +++ b/src/Core/hco_diagn_mod.F90 @@ -113,7 +113,6 @@ MODULE HCO_Diagn_Mod USE HCO_Arr_Mod USE HCO_Clock_Mod USE HCO_State_Mod, ONLY : HCO_State - USE MAPL_CommsMod, ONLY : MAPL_am_I_Root IMPLICIT NONE PRIVATE @@ -2013,7 +2012,6 @@ SUBROUTINE Diagn_UpdateDriver( HcoState, cID, cName, & IF ( AS /= 0 ) THEN CALL HCO_ERROR( HcoState%Config%Err,& 'Allocation error Arr3D', RC, THISLOC=LOC ) - if (MAPL_am_I_Root()) WRITE(*,*) 'Allocation error Arr3D_HP' RETURN ENDIF Arr3D = Array3D_HP @@ -2022,7 +2020,6 @@ SUBROUTINE Diagn_UpdateDriver( HcoState, cID, cName, & IF ( AS /= 0 ) THEN CALL HCO_ERROR( HcoState%Config%Err,& 'Allocation error Arr3D', RC, THISLOC=LOC ) - if (MAPL_am_I_Root()) WRITE(*,*) 'Allocation error Arr3D' RETURN ENDIF Arr3D = Array3D @@ -2036,7 +2033,6 @@ SUBROUTINE Diagn_UpdateDriver( HcoState, cID, cName, & IF ( AS /= 0 ) THEN CALL HCO_ERROR( HcoState%Config%Err,& 'Allocation error Arr2D', RC, THISLOC=LOC ) - if (MAPL_am_I_Root()) WRITE(*,*) 'Allocation error Arr2D_HP' RETURN ENDIF Arr2D = Array2D_HP @@ -2045,7 +2041,6 @@ SUBROUTINE Diagn_UpdateDriver( HcoState, cID, cName, & IF ( AS /= 0 ) THEN CALL HCO_ERROR( HcoState%Config%Err,& 'Allocation error Arr2D', RC, THISLOC=LOC ) - if (MAPL_am_I_Root()) WRITE(*,*) 'Allocation error Arr2D' RETURN ENDIF Arr2D = Array2D @@ -2142,7 +2137,6 @@ SUBROUTINE Diagn_UpdateDriver( HcoState, cID, cName, & ENDIF ELSE MSG = 'No array passed for updating ' // TRIM(ThisDiagn%cName) - if (MAPL_am_I_Root()) WRITE(*,*) MSG CALL HCO_ERROR ( HcoState%Config%Err, MSG, RC, THISLOC=LOC ) RETURN ENDIF @@ -2393,7 +2387,6 @@ SUBROUTINE Diagn_Get( HcoState, & ! Get collection number CALL DiagnCollection_DefineID( HcoState%Diagn, PS, RC, COL=COL, & ThisColl=ThisColl, HcoState=HcoState ) - IF ( RC /= HCO_SUCCESS .and. MAPL_am_I_Root()) WRITE(*,*) 'Failed to get collection number ', col IF ( RC /= HCO_SUCCESS ) RETURN ! Set AutoFill flag @@ -2476,7 +2469,6 @@ SUBROUTINE Diagn_Get( HcoState, & ! Before returning container, make sure its data is ready for output. IF ( ASSOCIATED (DgnCont ) ) THEN CALL DiagnCont_PrepareOutput ( HcoState, DgnCont, RC ) - IF ( RC /= HCO_SUCCESS .and. MAPL_am_I_Root() ) WRITE(*,*) 'DiagnCont_PrepareOutput returned RC: ', RC, ' for Col: ', DgnCont%CollectionID IF ( RC /= HCO_SUCCESS ) RETURN FLAG = HCO_SUCCESS @@ -3089,14 +3081,12 @@ SUBROUTINE DiagnCont_PrepareOutput( HcoState, DgnCont, RC ) !----------------------------------------------------------------------- CALL DiagnCollection_Find( HcoState%Diagn, DgnCont%CollectionID, & FOUND, RC, ThisColl=ThisColl ) - IF ( RC /= HCO_SUCCESS .and. MAPL_am_I_Root() ) WRITE(*,*) 'DiagnCollection_Find failed' IF ( RC /= HCO_SUCCESS ) RETURN ! This should never happen IF ( .NOT. FOUND .OR. .NOT. ASSOCIATED(ThisColl) ) THEN WRITE(MSG,*) 'Diagnostics ', TRIM(DgnCont%cName), ' has invalid ', & 'collection ID of ', DgnCont%CollectionID - if (MAPL_am_I_Root()) WRITE(*,*) MSG CALL HCO_ERROR( HcoState%Config%Err, MSG, RC, THISLOC=LOC ) RETURN ENDIF @@ -3110,7 +3100,6 @@ SUBROUTINE DiagnCont_PrepareOutput( HcoState, DgnCont, RC ) IF ( DgnCont%SpaceDim == 2 ) THEN CALL HCO_ArrAssert( DgnCont%Arr2D, ThisColl%NX, & ThisColl%NY, RC ) - IF ( RC /= HCO_SUCCESS .and. MAPL_am_I_Root() ) WRITE(*,*) 'HCO_ArrAssert 2d failed' IF ( RC /= HCO_SUCCESS ) RETURN ! Make sure it's zero @@ -3119,7 +3108,6 @@ SUBROUTINE DiagnCont_PrepareOutput( HcoState, DgnCont, RC ) ELSEIF ( DgnCont%SpaceDim == 3 ) THEN CALL HCO_ArrAssert( DgnCont%Arr3D, ThisColl%NX, & ThisColl%NY, ThisColl%NZ, RC ) - IF ( RC /= HCO_SUCCESS .and. MAPL_am_I_Root() ) WRITE(*,*) 'HCO_ArrAssert 3d failed' IF ( RC /= HCO_SUCCESS ) RETURN ! Make sure it's zero @@ -3166,7 +3154,6 @@ SUBROUTINE DiagnCont_PrepareOutput( HcoState, DgnCont, RC ) ! Get current month and year CALL HcoClock_Get( HcoState%Clock, cYYYY=YYYY, cMM=MM, RC=RC ) - IF ( RC /= HCO_SUCCESS .and. MAPL_am_I_Root() ) WRITE(*,*) 'HcoClock_Get failed' IF ( RC /= HCO_SUCCESS ) RETURN ! Days per year @@ -3200,7 +3187,6 @@ SUBROUTINE DiagnCont_PrepareOutput( HcoState, DgnCont, RC ) ELSE WRITE(MSG,*) 'Illegal time averaging of ', DgnCont%TimeAvg, & ' for diagnostics ', TRIM(DgnCont%cName) - if (MAPL_am_I_Root()) WRITE(*,*) TRIM(MSG) CALL HCO_ERROR( HcoState%Config%Err, MSG, RC, THISLOC=LOC ) RETURN ENDIF @@ -3210,7 +3196,6 @@ SUBROUTINE DiagnCont_PrepareOutput( HcoState, DgnCont, RC ) ! Error trap IF ( norm1 <= 0.0_hp ) THEN MSG = 'Illegal normalization factor: ' // TRIM(DgnCont%cName) - if (MAPL_am_I_Root()) WRITE(*,*) TRIM(MSG) CALL HCO_ERROR( HcoState%Config%Err, MSG, RC, THISLOC=LOC ) RETURN ENDIF @@ -3224,16 +3209,6 @@ SUBROUTINE DiagnCont_PrepareOutput( HcoState, DgnCont, RC ) ! For 3D: IF ( DgnCont%SpaceDim == 3 ) THEN - IF (MAPL_am_I_Root()) THEN - WRITE(*,*) ' ', TRIM(DgnCont%cName), ' 3D', & - ' Scaling: ', totScal - IF ( DgnCont%AreaFlag == 0) THEN - IF ( ASSOCIATED(DgnCont%Arr2D) ) THEN - WRITE(*,*) ' ', TRIM(DgnCont%cName), & - ' Appying area scaling factor ' - ENDIF - ENDIF - ENDIF DO J = 1, ThisColl%NY DO I = 1, ThisColl%NX @@ -3255,16 +3230,6 @@ SUBROUTINE DiagnCont_PrepareOutput( HcoState, DgnCont, RC ) ! For 2D: ELSEIF ( DgnCont%SpaceDim == 2 ) THEN - IF (MAPL_am_I_Root()) THEN - WRITE(*,*) ' ', TRIM(DgnCont%cName), ' 2D', & - ' Scaling: ', totScal - IF ( DgnCont%AreaFlag == 0) THEN - IF ( ASSOCIATED(DgnCont%Arr2D) ) THEN - WRITE(*,*) ' ', TRIM(DgnCont%cName), & - ' Appying area scaling factor ' - ENDIF - ENDIF - ENDIF DO J = 1, ThisColl%NY DO I = 1, ThisColl%NX @@ -3287,10 +3252,6 @@ SUBROUTINE DiagnCont_PrepareOutput( HcoState, DgnCont, RC ) ! For 1D: ELSE - IF (MAPL_am_I_Root()) THEN - WRITE(*,*) ' ', TRIM(DgnCont%cName), ' 1D', & - ' Scaling: ', totScal - ENDIF DgnCont%Scalar = DgnCont%Scalar * totscal ENDIF diff --git a/src/Core/hcoio_diagn_mod.F90 b/src/Core/hcoio_diagn_mod.F90 index 057956ad..f3986502 100644 --- a/src/Core/hcoio_diagn_mod.F90 +++ b/src/Core/hcoio_diagn_mod.F90 @@ -23,7 +23,6 @@ MODULE HCOIO_DIAGN_MOD ! !USES: ! USE HCO_ERROR_MOD - USE MAPL_CommsMod, ONLY : MAPL_am_I_Root IMPLICIT NONE PRIVATE @@ -116,7 +115,6 @@ SUBROUTINE HcoDiagn_Write( HcoState, Restart, RC ) ! To write restart (enforced) IF ( RESTART ) THEN - IF (MAPL_am_I_Root()) WRITE(*,*) "Doing HCOIO_DIAGN_WRITEOUT for restart" CALL HCOIO_DIAGN_WRITEOUT ( HcoState, & ForceWrite = .TRUE., & UsePrevTime = .FALSE., & @@ -124,12 +122,10 @@ SUBROUTINE HcoDiagn_Write( HcoState, Restart, RC ) RC = RC ) IF( RC /= HCO_SUCCESS) RETURN - IF (MAPL_am_I_Root()) WRITE(*,*) "HcoClock_SetLast" ! Set last flag for use below CALL HcoClock_SetLast ( HcoState%Clock, .TRUE., RC ) IF( RC /= HCO_SUCCESS) RETURN - IF (MAPL_am_I_Root()) WRITE(*,*) "Doing HCOIO_DIAGN_WRITEOUT default column" CALL HCOIO_DIAGN_WRITEOUT ( HcoState, & ForceWrite = .FALSE., & UsePrevTime = .FALSE., & @@ -139,7 +135,6 @@ SUBROUTINE HcoDiagn_Write( HcoState, Restart, RC ) #ifdef ADJOINT IF (HcoState%isAdjoint) THEN - if (MAPL_am_I_Root()) WRITE(*,*) "Doing HCOIO_DIAGN_WRITEOUT adjoint" CALL HCOIO_DIAGN_WRITEOUT ( HcoState, & ForceWrite = .FALSE., & UsePrevTime = .FALSE., & @@ -149,7 +144,6 @@ SUBROUTINE HcoDiagn_Write( HcoState, Restart, RC ) ENDIF #endif - if (MAPL_am_I_Root()) WRITE(*,*) "HcoClock_SetLast" ! Reset IsLast flag. This is to ensure that the last flag is not ! carried over (ckeller, 11/1/16). CALL HcoClock_SetLast ( HcoState%Clock, .FALSE., RC ) @@ -193,7 +187,6 @@ SUBROUTINE HcoDiagn_Write( HcoState, Restart, RC ) ! HCO_RestartWrite. (ckeller, 10/9/17) IF ( I == 2 ) CYCLE #endif - if (MAPL_am_I_Root()) WRITE(*,*) "Doing HCOIO_DIAGN_WRITEOUT, I = ", I ! Restart file CALL HCOIO_DIAGN_WRITEOUT ( HcoState, & @@ -201,7 +194,6 @@ SUBROUTINE HcoDiagn_Write( HcoState, Restart, RC ) UsePrevTime = .FALSE., & COL = COL, & RC = RC ) - if (MAPL_am_I_Root()) WRITE(*,*) "Back from HCOIO_DIAGN_WRITEOUT, RC = ", RC IF(RC /= HCO_SUCCESS) RETURN ENDDO ENDIF From 9517b513e321436ead6ae318e453d9968c07ea73 Mon Sep 17 00:00:00 2001 From: Colin Lee Date: Tue, 6 Apr 2021 15:33:07 -0700 Subject: [PATCH 3/4] Removed ifdef FALSE code --- .../Shared/hco_interface_common.F90 | 508 ------------------ 1 file changed, 508 deletions(-) diff --git a/src/Interfaces/Shared/hco_interface_common.F90 b/src/Interfaces/Shared/hco_interface_common.F90 index 0aa77e29..c49d8083 100644 --- a/src/Interfaces/Shared/hco_interface_common.F90 +++ b/src/Interfaces/Shared/hco_interface_common.F90 @@ -28,9 +28,6 @@ MODULE HCO_Interface_Common PUBLIC :: SetHcoTime PUBLIC :: GetHcoVal PUBLIC :: GetHcoDiagn -#ifdef FALSE - PUBLIC :: CALC_EMS_SF_ADJ -#endif ! ! !REMARKS: ! These utilities were mostly migrated from GEOS-Chem HCO_Interface_Mod. @@ -334,510 +331,5 @@ SUBROUTINE GetHcoDiagn ( HcoState, ExtState, DiagnName, & END SUBROUTINE GetHcoDiagn -#ifdef FALSE ! the adjoint-relevant code in here has been moved to - ! Hco_CalcEmis in #ifdef ADJOINT blocks so I'm disabling - ! this for now -! -! !DESCRIPTION: Subroutine CALC_EMS_SF_ADJ stores adjoint sensitvities to -! emissions categories -!\\ -!\\ -! !INTERFACE -! - SUBROUTINE CALC_EMS_SF_ADJ(SpcId, DT, Input_Opt, State_Chm, State_Diag, RC) -! -! !USES: -! - USE Input_Opt_Mod, ONLY : OptInput - USE State_Chm_Mod, ONLY : ChmState - USE State_Diag_Mod, ONLY : DgnState - USE HCO_STATE_MOD, ONLY : HCO_State - USE HCO_ARR_MOD, ONLY : HCO_ArrAssert - USE HCO_DATACONT_MOD, ONLY : ListCont_NextCont - USE HCO_FILEDATA_MOD, ONLY : FileData_ArrIsDefined - USE HCO_Calc_Mod, ONLY : GET_CURRENT_EMISSIONS_ADJ - USE HCO_Scale_Mod, ONLY : HCO_ScaleArr - USE HCO_Diagn_Mod - USE HCO_Types_Mod - - -! !INPUT PARAMETERS: -! - INTEGER, INTENT(IN ) :: spcId ! Which species - TYPE(OptInput), INTENT(IN ) :: Input_Opt ! Input opts - REAL(fp), INTENT(IN ) :: DT ! Time step [s] -! -! !INPUT/OUTPUT PARAMETERS: -! - TYPE(ChmState), INTENT(INOUT) :: State_Chm ! Chemistry state - TYPE(DgnState), INTENT(INOUT) :: State_Diag ! Diags State - INTEGER, INTENT(INOUT) :: RC ! Success/Failure -! -! !REVISION HISTORY: -! 18 May 2020 - C. Keller - Initial version -!EOP -!------------------------------------------------------------------------------ -!BOC -! -! !LOCAL VARIABLES: -! - ! Working pointers: list and data container - TYPE(ListCont), POINTER :: Lct - TYPE(DataCont), POINTER :: Dct - - ! Temporary emissions arrays - REAL(hp), TARGET :: SpcFlx( HcoState%NX, & - HcoState%NY, & - HcoState%NZ ) - REAL(hp), TARGET :: CatFlx( HcoState%NX, & - HcoState%NY, & - HcoState%NZ ) - REAL(hp), TARGET :: TmpFlx( HcoState%NX, & - HcoState%NY, & - HcoState%NZ ) - REAL(hp) :: Mask ( HcoState%NX, & - HcoState%NY, & - HcoState%NZ ) - REAL(hp) :: HirFlx( HcoState%NX, & - HcoState%NY, & - HcoState%NZ ) - REAL(hp) :: HirMsk( HcoState%NX, & - HcoState%NY, & - HcoState%NZ ) - - ! Integers - INTEGER :: ThisSpc, PrevSpc ! current and previous species ID - INTEGER :: ThisCat, PrevCat ! current and previous category - INTEGER :: ThisHir, PrevHir ! current and previous hierarchy - INTEGER :: SpcMin, SpcMax ! range of species to be considered - INTEGER :: CatMin, CatMax ! range of categories to be considered - INTEGER :: ExtNr ! Extension Nr to be used - INTEGER :: nI, nJ, nL - INTEGER :: nnSpec, FLAG - - LOGICAL :: Found, DoDiagn, EOL, UpdateCat, UseConc - - ! For error handling & verbose mode - CHARACTER(LEN=255) :: MSG - - ! testing / debugging - integer :: ix,iy - - !================================================================= - ! CALC_EMS_SF_ADJ begins here! - !================================================================= - ! Initialize - Lct => NULL() - Dct => NULL() - - - !----------------------------------------------------------------- - ! Initialize variables - !----------------------------------------------------------------- - - ! Initialize - SpcFlx(:,:,:) = 0.0_hp - CatFlx(:,:,:) = 0.0_hp - HirFlx(:,:,:) = 0.0_hp - HirMsk(:,:,:) = 0.0_hp - PrevSpc = -1 - PrevHir = -1 - PrevCat = -1 - nnSpec = 0 - - ! Pass emission grid dimensions - nI = HcoState%NX - nJ = HcoState%NY - nL = HcoState%NZ - - ! Pass calculation options - SpcMin = HcoState%Options%SpcMin !Lower species ID - SpcMax = HcoState%Options%SpcMax !Upper species ID - CatMin = HcoState%Options%CatMin !Lower emission category - CatMax = HcoState%Options%CatMax !Upper emission category - ExtNr = HcoState%Options%ExtNr !Extension number - DoDiagn = HcoState%Options%AutoFillDiagn !Write AutoFill diagnostics? - - ! I am not sure what do about units here, so set this to false for now - UseConc = .false. - - ! Enter routine - CALL HCO_ENTER (HcoState%Config%Err,'CALC_EMS_SF_ADJ (HCO_INTERFACE_MOD.F90)', RC ) - IF(RC /= HCO_SUCCESS) RETURN - - ! Verbose mode - IF ( HCO_IsVerb(HcoState%Config%Err,2) ) THEN - WRITE (MSG, *) 'Run HEMCO calculation w/ following options:' - CALL HCO_MSG ( HcoState%Config%Err, MSG ) - WRITE (MSG, "(A20,I5)") 'Extension number:', ExtNr - CALL HCO_MSG ( HcoState%Config%Err, MSG ) - WRITE (MSG, "(A20,I5,I5)") 'Tracer range:', SpcMin, SpcMax - CALL HCO_MSG ( HcoState%Config%Err, MSG ) - WRITE (MSG, "(A20,I5,I5)") 'Category range:', CatMin, CatMax - CALL HCO_MSG ( HcoState%Config%Err, MSG ) - WRITE (MSG, *) 'Auto diagnostics: ', DoDiagn - CALL HCO_MSG ( HcoState%Config%Err, MSG ) - ENDIF - - !================================================================= - ! Walk through all containers of EmisList and determine the - ! emissions for all containers that qualify for calculation. - ! The containers in EmisList are sorted by species, category and - ! hierarchy. This enables a straightforward, piece-by-piece - ! assembly of the final emission array (start with lowest - ! hierarchy emissions, then overwrite piece-by-piece with higher - ! hierarchy values). - !================================================================= - - ! Point to the head of the emissions linked list - EOL = .FALSE. ! End of list - Lct => NULL() - CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG ) - - ! Do until end of EmisList (==> loop over all emission containers) - DO - ! Have we reached the end of the list? - IF ( FLAG /= HCO_SUCCESS ) THEN - EOL = .TRUE. - ELSE - EOL = .FALSE. - ENDIF - - ! ------------------------------------------------------------ - ! Select container and update all working variables & arrays. - ! ------------------------------------------------------------ - IF ( .NOT. EOL ) THEN - - ! Dct is the current data container - Dct => Lct%Dct - if ( HCO_IsVerb(HcoState%Config%Err,2) ) THEN - WRITE(MSG,*) ' CALC_EMS_SF_ADJ processing container ', trim(Dct%cName) - CALL HCO_MSG(HcoState%Config%Err, MSG) - endif - - ! Check if this is a base field - IF ( Dct%DctType /= HCO_DCTTYPE_BASE ) THEN - CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG ) - CYCLE - ENDIF - - ! Sanity check: Make sure this container holds data. - ! 'Empty' containers are possible if the simulation time - ! is outside of the specified data time range and time - ! slice cycling is deactivated (CycleFlag > 1). - IF( .NOT. FileData_ArrIsDefined(Lct%Dct%Dta) ) THEN - CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG ) - CYCLE - ENDIF - - ! Check if this is the specified extension number - IF ( Dct%ExtNr /= ExtNr ) THEN - CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG ) - CYCLE - ENDIF - - ! Advance to next container if the species ID is outside - ! the specified species range (SpcMin - SpcMax). Consider - ! all species above SpcMin if SpcMax is negative! - IF( ( Dct%HcoID < SpcMin ) .OR. & - ( (Dct%HcoID > SpcMax) .AND. (SpcMax > 0) ) ) THEN - CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG ) - CYCLE - ENDIF - - ! Advance to next emission field if the emission category of - ! the current container is outside of the specified species - ! range (CatMin - CatMax). Consider all categories above CatMin - ! if CatMax is negative! - IF( ( Dct%Cat < CatMin ) .OR. & - ( (Dct%Cat > CatMax) .AND. (CatMax > 0) ) ) THEN - CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG ) - CYCLE - ENDIF - - ! Check if this container holds data in the desired unit format, - ! i.e. concentration data if UseConc is enabled, emission data - ! otherwise. - IF ( UseConc .NEQV. Dct%Dta%IsConc ) THEN - CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG ) - CYCLE - ENDIF - - ! Check if this is our species - IF ( Dct%HcoID /= SpcId ) THEN - CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG ) - CYCLE - ENDIF - - ! Update working variables - ThisSpc = Dct%HcoID - ThisCat = Dct%Cat - ThisHir = Dct%Hier - - ! If end of list, use dummy values for ThisSpc, ThisCat and ThisHir - ! to make sure that emissions are added to HEMCO in the section - ! below! - ELSE - ThisSpc = -1 - ThisCat = -1 - ThisHir = -1 - ENDIF - - !-------------------------------------------------------------------- - ! Before computing emissions of current data container make sure that - ! emissions of previous container are properly archived. - !-------------------------------------------------------------------- - - ! Add emissions on hierarchy level to the category flux array. Do - ! this only if this is a new species, a new category or a new - ! hierarchy level. - ! Note: no need to add to diagnostics because hierarchy level - ! diagnostics are filled right after computing the emissions of - ! a given data container (towards the end of the DO loop). - IF ( (ThisHir /= PrevHir) .OR. & - (ThisSpc /= PrevSpc) .OR. & - (ThisCat /= PrevCat) ) THEN - - ! Add hierarchy level emissions to category array over the - ! covered regions. - CatFlx = ( (1.0_hp - HirMsk) * CatFlx ) + HirFlx - - ! Reset - HirFlx = 0.0_hp - HirMsk = 0.0_hp - ENDIF - - !-------------------------------------------------------------------- - ! If this is a new species or category, pass the previously collected - ! emissions to the species array. Update diagnostics at category level. - ! Skip this step for first species, i.e. if PrevSpc is still -1. - !-------------------------------------------------------------------- - UpdateCat = .FALSE. - IF ( ThisCat /= PrevCat ) UpdateCat = .TRUE. - IF ( ThisSpc /= PrevSpc ) UpdateCat = .TRUE. - IF ( PrevCat <= 0 .OR. PrevSpc <= 0 ) UpdateCat = .FALSE. - IF ( UpdateCat ) THEN - - ! CatFlx holds the emissions for this category. Pass this to - ! the species array SpcFlx. - SpcFlx(:,:,:) = SpcFlx(:,:,:) + CatFlx(:,:,:) - - IF (Input_Opt%IS_FD_SPOT_THIS_PET) THEN - WRITE(*,*) PrevCat, ' CatFlx(IFD,JFD) = ', & - SUM(CatFlx(Input_Opt%IFD,Input_Opt%JFD,:)) - ENDIF - - ! verbose - IF ( HCO_IsVerb(HcoState%Config%Err,3) ) THEN - WRITE(MSG,*) 'Added category emissions to species array: ' - CALL HCO_MSG(HcoState%Config%Err,MSG) - WRITE(MSG,*) 'Species : ', PrevSpc - CALL HCO_MSG(HcoState%Config%Err,MSG) - WRITE(MSG,*) 'Category : ', PrevCat - CALL HCO_MSG(HcoState%Config%Err,MSG) - WRITE(MSG,*) 'Cat. emissions: ', SUM(CatFlx) - CALL HCO_MSG(HcoState%Config%Err,MSG) - WRITE(MSG,*) 'Spc. emissions: ', SUM(SpcFlx) - CALL HCO_MSG(HcoState%Config%Err,MSG) - ENDIF - CatFlx(:,:,:) = DT * State_Chm%SpeciesAdj(:,:,:, SpcId) * CatFlx(:,:,:) - ! verbose - IF ( HCO_IsVerb(HcoState%Config%Err,3) ) THEN - WRITE(MSG,*) 'Cat. emissions: ', SUM(CatFlx) - CALL HCO_MSG(HcoState%Config%Err,MSG) - ENDIF - IF (Input_Opt%IS_FD_SPOT_THIS_PET) THEN - WRITE(*,*) PrevCat, ' CatFlxAdj(IFD,JFD) = ', & - SUM(CatFlx(Input_Opt%IFD,Input_Opt%JFD,:)) - ENDIF - - ! Add category emissions to diagnostics at category level - ! (only if defined in the diagnostics list). - IF ( Diagn_AutoFillLevelDefined(HcoState%Diagn,3) .AND. DoDiagn ) THEN - CALL Diagn_Update( HcoState, ExtNr=ExtNr, & - Cat=PrevCat, Hier=-1, HcoID=PrevSpc, & - AutoFill=1, Array3D=CatFlx, & - COL=HcoState%Diagn%HcoDiagnIDAdjoint, RC=RC ) - IF ( RC /= HCO_SUCCESS ) RETURN - ENDIF - - ! Reset CatFlx array and the previously used hierarchy - ! ==> Emission hierarchies are only important within the - ! same category, hence always start over at lowest hierarchy - ! when entering a new category. - CatFlx(:,:,:) = 0.0_hp - PrevHir = -1 - ENDIF - - !-------------------------------------------------------------------- - ! If this is a new species, pass previously calculated emissions - ! to the final emissions array in HcoState. - ! Update diagnostics at extension number level. - ! Don't do before first emission calculation, i.e. if PrevSpc - ! is still the initialized value of -1! - !-------------------------------------------------------------------- - IF ( ThisSpc /= PrevSpc .AND. PrevSpc > 0 ) THEN - - IF (Input_Opt%IS_FD_SPOT_THIS_PET) THEN - WRITE(*,*) PrevCat, ' SpcFlx(IFD,JFD) = ', & - SUM(SpcFlx(Input_Opt%IFD,Input_Opt%JFD,:)) - ENDIF - ! testing only - IF ( HCO_IsVerb(HcoState%Config%Err,3) ) THEN - WRITE(MSG,*) 'Added total emissions to output array: ' - CALL HCO_MSG(HcoState%Config%Err,MSG) - WRITE(MSG,*) 'Species: ', PrevSpc - CALL HCO_MSG(HcoState%Config%Err,MSG) - WRITE(MSG,*) 'SpcFlx : ', SUM(SpcFlx) - CALL HCO_MSG(HcoState%Config%Err,MSG) - ENDIF - - ! Add to diagnostics at extension number level. - ! The same diagnostics may be updated multiple times during - ! the same time step, continuously adding emissions to it. - IF ( Diagn_AutoFillLevelDefined(HcoState%Diagn,2) .AND. DoDiagn ) THEN - CALL Diagn_Update(HcoState, ExtNr=ExtNr, & - Cat=-1, Hier=-1, HcoID=PrevSpc, & - AutoFill=1,Array3D=SpcFlx, & - COL=HcoState%Diagn%HcoDiagnIDAdjoint, RC=RC ) - IF ( RC /= HCO_SUCCESS ) RETURN - ENDIF - - ! Reset arrays and previous hierarchy. - SpcFlx(:,:,:) = 0.0_hp - PrevCat = -1 - PrevHir = -1 - ENDIF - - !-------------------------------------------------------------------- - ! Exit DO loop here if end of list - !-------------------------------------------------------------------- - IF ( EOL ) EXIT - - !-------------------------------------------------------------------- - ! Update/archive information on species level if needed - !-------------------------------------------------------------------- - IF ( ThisSpc /= PrevSpc .AND. ThisSpc > 0 ) THEN - - ! Update number of species for which emissions have been - ! calculated. - nnSpec = nnSpec + 1 - - ! verbose mode - IF ( HCO_IsVerb(HcoState%Config%Err,2) ) THEN - WRITE(MSG,*) 'Calculating emissions for species ', & - TRIM(HcoState%Spc(ThisSpc)%SpcName) - CALL HCO_MSG( HcoState%Config%Err, MSG, SEP1='-', SEP2='-' ) - ENDIF - ENDIF - - ! verbose mode - IF ( HCO_IsVerb(HcoState%Config%Err,2) ) THEN - WRITE(MSG,*) 'Dct has ', & - Dct%Dta%SpaceDim, ' dimensions.' - CALL HCO_MSG( HcoState%Config%Err, MSG, SEP1='-', SEP2='-' ) - ENDIF - !-------------------------------------------------------------------- - ! Get current emissions and write into TmpFlx array. The array Mask - ! denotes all valid grid boxes for this inventory. - !-------------------------------------------------------------------- - TmpFlx(:,:,:) = 0.0_hp - CALL GET_CURRENT_EMISSIONS_ADJ( HcoState, Dct, & - nI, nJ, nL, TmpFlx, Mask, RC ) - IF ( RC /= HCO_SUCCESS ) RETURN - - ! Eventually add universal scale factor - CALL HCO_ScaleArr( HcoState, ThisSpc, TmpFlx, RC ) - IF ( RC /= HCO_SUCCESS ) RETURN - - ! Check for negative values according to the corresponding setting - ! in the configuration file: 2 means allow negative values, 1 means - ! set to zero and prompt a warning, else return with error. - IF ( HcoState%Options%NegFlag /= 2 ) THEN - - IF ( ANY(TmpFlx < 0.0_hp) ) THEN - - ! Set to zero and prompt warning - IF ( HcoState%Options%NegFlag == 1 ) THEN - WHERE ( TmpFlx < 0.0_hp ) TmpFlx = 0.0_hp - MSG = 'Negative emissions set to zero: '// TRIM(Dct%cName) - CALL HCO_WARNING( HcoState%Config%Err, MSG, RC ) - - ! Return with error - ELSE - MSG = 'Negative emissions in: '// TRIM(Dct%cName) // '. ' // & - 'To allow negatives, edit settings in the configuration file.' - CALL HCO_ERROR( HcoState%Config%Err, MSG, RC ) - RETURN - ENDIF - ENDIF - ENDIF - - ! ------------------------------------------------------------ - ! Collect all emissions of the same category (and species) on - ! the hierarchy level into array HirFlx. HirMsk contains the - ! combined covered region. That is, if there are two regional - ! inventories with the same hierarchy HirMsk will cover both - ! of these regions. - ! The specified field hierarchies determine whether the - ! temporary emissions are added (if hierarchy is the same - ! as the previously used hierarchy), or if they overwrite the - ! previous values in HirFlx (if hierarchy is higher than the - ! previous hierarchy). - ! ------------------------------------------------------------ - - ! Add emissions to the hierarchy array HirFlx if this hierarchy - ! is the same as previous hierarchy - IF ( ThisHir == PrevHir ) THEN - HirFlx = HirFlx + TmpFlx - HirMsk = HirMsk + Mask - - ! Make sure mask values do not exceed 1.0 - WHERE(HirMsk > 1.0 ) HirMsk = 1.0 - - ! If hierarchy is larger than those of the previously used - ! fields, overwrite HirFlx with new values. - ELSE - - HirFlx = TmpFlx - HirMsk = Mask - - ENDIF - - ! Update diagnostics at hierarchy level. Make sure that only - ! positive values are used. - ! The same diagnostics may be updated multiple times - ! during the same time step, continuously adding - ! emissions to it. - ! Now remove PosOnly flag. TmpFlx is initialized to zero, so it's - ! ok to keep negative values (ckeller, 7/12/15). - IF ( Diagn_AutoFillLevelDefined(HcoState%Diagn,4) .AND. DoDiagn ) THEN - CALL Diagn_Update( HcoState, ExtNr=ExtNr, & - Cat=ThisCat,Hier=ThisHir, HcoID=ThisSpc, & - !AutoFill=1, Array3D=TmpFlx, PosOnly=.TRUE.,& - AutoFill=1, Array3D=TmpFlx, & - COL=-1, RC=RC ) - IF ( RC /= HCO_SUCCESS ) RETURN - ENDIF - - ! Update previously used species, category and hierarchy - PrevSpc = ThisSpc - PrevCat = ThisCat - PrevHir = ThisHir - - ! Advance to next emission container - CALL ListCont_NextCont( HcoState%EmisList, Lct, FLAG ) - - ENDDO ! Loop over EmisList - - ! Make sure internal pointers are nullified - Lct => NULL() - Dct => NULL() - - ! Leave w/ success - CALL HCO_LEAVE ( HcoState%Config%Err, RC ) - - END SUBROUTINE CALC_EMS_SF_ADJ -#endif !EOC END MODULE HCO_Interface_Common From b2d3a5b583178348b9f1698b1f8a8fb2b619cc1b Mon Sep 17 00:00:00 2001 From: Colin Lee Date: Wed, 7 Apr 2021 09:53:47 -0700 Subject: [PATCH 4/4] Removed dead code and commented out debug prints --- src/Core/hcoio_write_esmf_mod.F90 | 10 +- .../Shared/hco_interface_common.F90 | 508 ------------------ 2 files changed, 5 insertions(+), 513 deletions(-) diff --git a/src/Core/hcoio_write_esmf_mod.F90 b/src/Core/hcoio_write_esmf_mod.F90 index 2865eb12..74df14d6 100644 --- a/src/Core/hcoio_write_esmf_mod.F90 +++ b/src/Core/hcoio_write_esmf_mod.F90 @@ -151,17 +151,17 @@ SUBROUTINE HCOIO_WRITE_ESMF ( HcoState, RC, OnlyIfFirst, COL ) ThisDiagn => NULL() DO WHILE ( .TRUE. ) - IF (MAPL_am_I_Root()) WRITE(*,*) "Getting next diagnostic in list" + ! IF (MAPL_am_I_Root()) WRITE(*,*) "Getting next diagnostic in list" ! Get next diagnostics in list. This will return the next ! diagnostics container that contains content to be written ! out on this time step. CALL Diagn_Get ( HcoState, EOI, ThisDiagn, FLAG, RC, COL=PS ) - IF (MAPL_am_I_Root()) THEN - IF ( RC == HCO_SUCCESS .and. FLAG == HCO_SUCCESS ) WRITE(*,*) "Got it! Name: ", TRIM(ThisDiagn%cName) - IF ( RC /= HCO_SUCCESS) WRITE(*,*) "Fail! RC: ", RC, " Flag: ", FLAG - ENDIF + ! IF (MAPL_am_I_Root()) THEN + ! IF ( RC == HCO_SUCCESS .and. FLAG == HCO_SUCCESS ) WRITE(*,*) "Got it! Name: ", TRIM(ThisDiagn%cName) + ! IF ( RC /= HCO_SUCCESS) WRITE(*,*) "Fail! RC: ", RC, " Flag: ", FLAG + ! ENDIF IF ( RC /= HCO_SUCCESS ) RETURN IF ( FLAG /= HCO_SUCCESS ) EXIT diff --git a/src/Interfaces/Shared/hco_interface_common.F90 b/src/Interfaces/Shared/hco_interface_common.F90 index 0aa77e29..c49d8083 100644 --- a/src/Interfaces/Shared/hco_interface_common.F90 +++ b/src/Interfaces/Shared/hco_interface_common.F90 @@ -28,9 +28,6 @@ MODULE HCO_Interface_Common PUBLIC :: SetHcoTime PUBLIC :: GetHcoVal PUBLIC :: GetHcoDiagn -#ifdef FALSE - PUBLIC :: CALC_EMS_SF_ADJ -#endif ! ! !REMARKS: ! These utilities were mostly migrated from GEOS-Chem HCO_Interface_Mod. @@ -334,510 +331,5 @@ SUBROUTINE GetHcoDiagn ( HcoState, ExtState, DiagnName, & END SUBROUTINE GetHcoDiagn -#ifdef FALSE ! the adjoint-relevant code in here has been moved to - ! Hco_CalcEmis in #ifdef ADJOINT blocks so I'm disabling - ! this for now -! -! !DESCRIPTION: Subroutine CALC_EMS_SF_ADJ stores adjoint sensitvities to -! emissions categories -!\\ -!\\ -! !INTERFACE -! - SUBROUTINE CALC_EMS_SF_ADJ(SpcId, DT, Input_Opt, State_Chm, State_Diag, RC) -! -! !USES: -! - USE Input_Opt_Mod, ONLY : OptInput - USE State_Chm_Mod, ONLY : ChmState - USE State_Diag_Mod, ONLY : DgnState - USE HCO_STATE_MOD, ONLY : HCO_State - USE HCO_ARR_MOD, ONLY : HCO_ArrAssert - USE HCO_DATACONT_MOD, ONLY : ListCont_NextCont - USE HCO_FILEDATA_MOD, ONLY : FileData_ArrIsDefined - USE HCO_Calc_Mod, ONLY : GET_CURRENT_EMISSIONS_ADJ - USE HCO_Scale_Mod, ONLY : HCO_ScaleArr - USE HCO_Diagn_Mod - USE HCO_Types_Mod - - -! !INPUT PARAMETERS: -! - INTEGER, INTENT(IN ) :: spcId ! Which species - TYPE(OptInput), INTENT(IN ) :: Input_Opt ! Input opts - REAL(fp), INTENT(IN ) :: DT ! Time step [s] -! -! !INPUT/OUTPUT PARAMETERS: -! - TYPE(ChmState), INTENT(INOUT) :: State_Chm ! Chemistry state - TYPE(DgnState), INTENT(INOUT) :: State_Diag ! Diags State - INTEGER, INTENT(INOUT) :: RC ! Success/Failure -! -! !REVISION HISTORY: -! 18 May 2020 - C. Keller - Initial version -!EOP -!------------------------------------------------------------------------------ -!BOC -! -! !LOCAL VARIABLES: -! - ! Working pointers: list and data container - TYPE(ListCont), POINTER :: Lct - TYPE(DataCont), POINTER :: Dct - - ! Temporary emissions arrays - REAL(hp), TARGET :: SpcFlx( HcoState%NX, & - HcoState%NY, & - HcoState%NZ ) - REAL(hp), TARGET :: CatFlx( HcoState%NX, & - HcoState%NY, & - HcoState%NZ ) - REAL(hp), TARGET :: TmpFlx( HcoState%NX, & - HcoState%NY, & - HcoState%NZ ) - REAL(hp) :: Mask ( HcoState%NX, & - HcoState%NY, & - HcoState%NZ ) - REAL(hp) :: HirFlx( HcoState%NX, & - HcoState%NY, & - HcoState%NZ ) - REAL(hp) :: HirMsk( HcoState%NX, & - HcoState%NY, & - HcoState%NZ ) - - ! Integers - INTEGER :: ThisSpc, PrevSpc ! current and previous species ID - INTEGER :: ThisCat, PrevCat ! current and previous category - INTEGER :: ThisHir, PrevHir ! current and previous hierarchy - INTEGER :: SpcMin, SpcMax ! range of species to be considered - INTEGER :: CatMin, CatMax ! range of categories to be considered - INTEGER :: ExtNr ! Extension Nr to be used - INTEGER :: nI, nJ, nL - INTEGER :: nnSpec, FLAG - - LOGICAL :: Found, DoDiagn, EOL, UpdateCat, UseConc - - ! For error handling & verbose mode - CHARACTER(LEN=255) :: MSG - - ! testing / debugging - integer :: ix,iy - - !================================================================= - ! CALC_EMS_SF_ADJ begins here! - !================================================================= - ! Initialize - Lct => NULL() - Dct => NULL() - - - !----------------------------------------------------------------- - ! Initialize variables - !----------------------------------------------------------------- - - ! Initialize - SpcFlx(:,:,:) = 0.0_hp - CatFlx(:,:,:) = 0.0_hp - HirFlx(:,:,:) = 0.0_hp - HirMsk(:,:,:) = 0.0_hp - PrevSpc = -1 - PrevHir = -1 - PrevCat = -1 - nnSpec = 0 - - ! Pass emission grid dimensions - nI = HcoState%NX - nJ = HcoState%NY - nL = HcoState%NZ - - ! Pass calculation options - SpcMin = HcoState%Options%SpcMin !Lower species ID - SpcMax = HcoState%Options%SpcMax !Upper species ID - CatMin = HcoState%Options%CatMin !Lower emission category - CatMax = HcoState%Options%CatMax !Upper emission category - ExtNr = HcoState%Options%ExtNr !Extension number - DoDiagn = HcoState%Options%AutoFillDiagn !Write AutoFill diagnostics? - - ! I am not sure what do about units here, so set this to false for now - UseConc = .false. - - ! Enter routine - CALL HCO_ENTER (HcoState%Config%Err,'CALC_EMS_SF_ADJ (HCO_INTERFACE_MOD.F90)', RC ) - IF(RC /= HCO_SUCCESS) RETURN - - ! Verbose mode - IF ( HCO_IsVerb(HcoState%Config%Err,2) ) THEN - WRITE (MSG, *) 'Run HEMCO calculation w/ following options:' - CALL HCO_MSG ( HcoState%Config%Err, MSG ) - WRITE (MSG, "(A20,I5)") 'Extension number:', ExtNr - CALL HCO_MSG ( HcoState%Config%Err, MSG ) - WRITE (MSG, "(A20,I5,I5)") 'Tracer range:', SpcMin, SpcMax - CALL HCO_MSG ( HcoState%Config%Err, MSG ) - WRITE (MSG, "(A20,I5,I5)") 'Category range:', CatMin, CatMax - CALL HCO_MSG ( HcoState%Config%Err, MSG ) - WRITE (MSG, *) 'Auto diagnostics: ', DoDiagn - CALL HCO_MSG ( HcoState%Config%Err, MSG ) - ENDIF - - !================================================================= - ! Walk through all containers of EmisList and determine the - ! emissions for all containers that qualify for calculation. - ! The containers in EmisList are sorted by species, category and - ! hierarchy. This enables a straightforward, piece-by-piece - ! assembly of the final emission array (start with lowest - ! hierarchy emissions, then overwrite piece-by-piece with higher - ! hierarchy values). - !================================================================= - - ! Point to the head of the emissions linked list - EOL = .FALSE. ! End of list - Lct => NULL() - CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG ) - - ! Do until end of EmisList (==> loop over all emission containers) - DO - ! Have we reached the end of the list? - IF ( FLAG /= HCO_SUCCESS ) THEN - EOL = .TRUE. - ELSE - EOL = .FALSE. - ENDIF - - ! ------------------------------------------------------------ - ! Select container and update all working variables & arrays. - ! ------------------------------------------------------------ - IF ( .NOT. EOL ) THEN - - ! Dct is the current data container - Dct => Lct%Dct - if ( HCO_IsVerb(HcoState%Config%Err,2) ) THEN - WRITE(MSG,*) ' CALC_EMS_SF_ADJ processing container ', trim(Dct%cName) - CALL HCO_MSG(HcoState%Config%Err, MSG) - endif - - ! Check if this is a base field - IF ( Dct%DctType /= HCO_DCTTYPE_BASE ) THEN - CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG ) - CYCLE - ENDIF - - ! Sanity check: Make sure this container holds data. - ! 'Empty' containers are possible if the simulation time - ! is outside of the specified data time range and time - ! slice cycling is deactivated (CycleFlag > 1). - IF( .NOT. FileData_ArrIsDefined(Lct%Dct%Dta) ) THEN - CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG ) - CYCLE - ENDIF - - ! Check if this is the specified extension number - IF ( Dct%ExtNr /= ExtNr ) THEN - CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG ) - CYCLE - ENDIF - - ! Advance to next container if the species ID is outside - ! the specified species range (SpcMin - SpcMax). Consider - ! all species above SpcMin if SpcMax is negative! - IF( ( Dct%HcoID < SpcMin ) .OR. & - ( (Dct%HcoID > SpcMax) .AND. (SpcMax > 0) ) ) THEN - CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG ) - CYCLE - ENDIF - - ! Advance to next emission field if the emission category of - ! the current container is outside of the specified species - ! range (CatMin - CatMax). Consider all categories above CatMin - ! if CatMax is negative! - IF( ( Dct%Cat < CatMin ) .OR. & - ( (Dct%Cat > CatMax) .AND. (CatMax > 0) ) ) THEN - CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG ) - CYCLE - ENDIF - - ! Check if this container holds data in the desired unit format, - ! i.e. concentration data if UseConc is enabled, emission data - ! otherwise. - IF ( UseConc .NEQV. Dct%Dta%IsConc ) THEN - CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG ) - CYCLE - ENDIF - - ! Check if this is our species - IF ( Dct%HcoID /= SpcId ) THEN - CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG ) - CYCLE - ENDIF - - ! Update working variables - ThisSpc = Dct%HcoID - ThisCat = Dct%Cat - ThisHir = Dct%Hier - - ! If end of list, use dummy values for ThisSpc, ThisCat and ThisHir - ! to make sure that emissions are added to HEMCO in the section - ! below! - ELSE - ThisSpc = -1 - ThisCat = -1 - ThisHir = -1 - ENDIF - - !-------------------------------------------------------------------- - ! Before computing emissions of current data container make sure that - ! emissions of previous container are properly archived. - !-------------------------------------------------------------------- - - ! Add emissions on hierarchy level to the category flux array. Do - ! this only if this is a new species, a new category or a new - ! hierarchy level. - ! Note: no need to add to diagnostics because hierarchy level - ! diagnostics are filled right after computing the emissions of - ! a given data container (towards the end of the DO loop). - IF ( (ThisHir /= PrevHir) .OR. & - (ThisSpc /= PrevSpc) .OR. & - (ThisCat /= PrevCat) ) THEN - - ! Add hierarchy level emissions to category array over the - ! covered regions. - CatFlx = ( (1.0_hp - HirMsk) * CatFlx ) + HirFlx - - ! Reset - HirFlx = 0.0_hp - HirMsk = 0.0_hp - ENDIF - - !-------------------------------------------------------------------- - ! If this is a new species or category, pass the previously collected - ! emissions to the species array. Update diagnostics at category level. - ! Skip this step for first species, i.e. if PrevSpc is still -1. - !-------------------------------------------------------------------- - UpdateCat = .FALSE. - IF ( ThisCat /= PrevCat ) UpdateCat = .TRUE. - IF ( ThisSpc /= PrevSpc ) UpdateCat = .TRUE. - IF ( PrevCat <= 0 .OR. PrevSpc <= 0 ) UpdateCat = .FALSE. - IF ( UpdateCat ) THEN - - ! CatFlx holds the emissions for this category. Pass this to - ! the species array SpcFlx. - SpcFlx(:,:,:) = SpcFlx(:,:,:) + CatFlx(:,:,:) - - IF (Input_Opt%IS_FD_SPOT_THIS_PET) THEN - WRITE(*,*) PrevCat, ' CatFlx(IFD,JFD) = ', & - SUM(CatFlx(Input_Opt%IFD,Input_Opt%JFD,:)) - ENDIF - - ! verbose - IF ( HCO_IsVerb(HcoState%Config%Err,3) ) THEN - WRITE(MSG,*) 'Added category emissions to species array: ' - CALL HCO_MSG(HcoState%Config%Err,MSG) - WRITE(MSG,*) 'Species : ', PrevSpc - CALL HCO_MSG(HcoState%Config%Err,MSG) - WRITE(MSG,*) 'Category : ', PrevCat - CALL HCO_MSG(HcoState%Config%Err,MSG) - WRITE(MSG,*) 'Cat. emissions: ', SUM(CatFlx) - CALL HCO_MSG(HcoState%Config%Err,MSG) - WRITE(MSG,*) 'Spc. emissions: ', SUM(SpcFlx) - CALL HCO_MSG(HcoState%Config%Err,MSG) - ENDIF - CatFlx(:,:,:) = DT * State_Chm%SpeciesAdj(:,:,:, SpcId) * CatFlx(:,:,:) - ! verbose - IF ( HCO_IsVerb(HcoState%Config%Err,3) ) THEN - WRITE(MSG,*) 'Cat. emissions: ', SUM(CatFlx) - CALL HCO_MSG(HcoState%Config%Err,MSG) - ENDIF - IF (Input_Opt%IS_FD_SPOT_THIS_PET) THEN - WRITE(*,*) PrevCat, ' CatFlxAdj(IFD,JFD) = ', & - SUM(CatFlx(Input_Opt%IFD,Input_Opt%JFD,:)) - ENDIF - - ! Add category emissions to diagnostics at category level - ! (only if defined in the diagnostics list). - IF ( Diagn_AutoFillLevelDefined(HcoState%Diagn,3) .AND. DoDiagn ) THEN - CALL Diagn_Update( HcoState, ExtNr=ExtNr, & - Cat=PrevCat, Hier=-1, HcoID=PrevSpc, & - AutoFill=1, Array3D=CatFlx, & - COL=HcoState%Diagn%HcoDiagnIDAdjoint, RC=RC ) - IF ( RC /= HCO_SUCCESS ) RETURN - ENDIF - - ! Reset CatFlx array and the previously used hierarchy - ! ==> Emission hierarchies are only important within the - ! same category, hence always start over at lowest hierarchy - ! when entering a new category. - CatFlx(:,:,:) = 0.0_hp - PrevHir = -1 - ENDIF - - !-------------------------------------------------------------------- - ! If this is a new species, pass previously calculated emissions - ! to the final emissions array in HcoState. - ! Update diagnostics at extension number level. - ! Don't do before first emission calculation, i.e. if PrevSpc - ! is still the initialized value of -1! - !-------------------------------------------------------------------- - IF ( ThisSpc /= PrevSpc .AND. PrevSpc > 0 ) THEN - - IF (Input_Opt%IS_FD_SPOT_THIS_PET) THEN - WRITE(*,*) PrevCat, ' SpcFlx(IFD,JFD) = ', & - SUM(SpcFlx(Input_Opt%IFD,Input_Opt%JFD,:)) - ENDIF - ! testing only - IF ( HCO_IsVerb(HcoState%Config%Err,3) ) THEN - WRITE(MSG,*) 'Added total emissions to output array: ' - CALL HCO_MSG(HcoState%Config%Err,MSG) - WRITE(MSG,*) 'Species: ', PrevSpc - CALL HCO_MSG(HcoState%Config%Err,MSG) - WRITE(MSG,*) 'SpcFlx : ', SUM(SpcFlx) - CALL HCO_MSG(HcoState%Config%Err,MSG) - ENDIF - - ! Add to diagnostics at extension number level. - ! The same diagnostics may be updated multiple times during - ! the same time step, continuously adding emissions to it. - IF ( Diagn_AutoFillLevelDefined(HcoState%Diagn,2) .AND. DoDiagn ) THEN - CALL Diagn_Update(HcoState, ExtNr=ExtNr, & - Cat=-1, Hier=-1, HcoID=PrevSpc, & - AutoFill=1,Array3D=SpcFlx, & - COL=HcoState%Diagn%HcoDiagnIDAdjoint, RC=RC ) - IF ( RC /= HCO_SUCCESS ) RETURN - ENDIF - - ! Reset arrays and previous hierarchy. - SpcFlx(:,:,:) = 0.0_hp - PrevCat = -1 - PrevHir = -1 - ENDIF - - !-------------------------------------------------------------------- - ! Exit DO loop here if end of list - !-------------------------------------------------------------------- - IF ( EOL ) EXIT - - !-------------------------------------------------------------------- - ! Update/archive information on species level if needed - !-------------------------------------------------------------------- - IF ( ThisSpc /= PrevSpc .AND. ThisSpc > 0 ) THEN - - ! Update number of species for which emissions have been - ! calculated. - nnSpec = nnSpec + 1 - - ! verbose mode - IF ( HCO_IsVerb(HcoState%Config%Err,2) ) THEN - WRITE(MSG,*) 'Calculating emissions for species ', & - TRIM(HcoState%Spc(ThisSpc)%SpcName) - CALL HCO_MSG( HcoState%Config%Err, MSG, SEP1='-', SEP2='-' ) - ENDIF - ENDIF - - ! verbose mode - IF ( HCO_IsVerb(HcoState%Config%Err,2) ) THEN - WRITE(MSG,*) 'Dct has ', & - Dct%Dta%SpaceDim, ' dimensions.' - CALL HCO_MSG( HcoState%Config%Err, MSG, SEP1='-', SEP2='-' ) - ENDIF - !-------------------------------------------------------------------- - ! Get current emissions and write into TmpFlx array. The array Mask - ! denotes all valid grid boxes for this inventory. - !-------------------------------------------------------------------- - TmpFlx(:,:,:) = 0.0_hp - CALL GET_CURRENT_EMISSIONS_ADJ( HcoState, Dct, & - nI, nJ, nL, TmpFlx, Mask, RC ) - IF ( RC /= HCO_SUCCESS ) RETURN - - ! Eventually add universal scale factor - CALL HCO_ScaleArr( HcoState, ThisSpc, TmpFlx, RC ) - IF ( RC /= HCO_SUCCESS ) RETURN - - ! Check for negative values according to the corresponding setting - ! in the configuration file: 2 means allow negative values, 1 means - ! set to zero and prompt a warning, else return with error. - IF ( HcoState%Options%NegFlag /= 2 ) THEN - - IF ( ANY(TmpFlx < 0.0_hp) ) THEN - - ! Set to zero and prompt warning - IF ( HcoState%Options%NegFlag == 1 ) THEN - WHERE ( TmpFlx < 0.0_hp ) TmpFlx = 0.0_hp - MSG = 'Negative emissions set to zero: '// TRIM(Dct%cName) - CALL HCO_WARNING( HcoState%Config%Err, MSG, RC ) - - ! Return with error - ELSE - MSG = 'Negative emissions in: '// TRIM(Dct%cName) // '. ' // & - 'To allow negatives, edit settings in the configuration file.' - CALL HCO_ERROR( HcoState%Config%Err, MSG, RC ) - RETURN - ENDIF - ENDIF - ENDIF - - ! ------------------------------------------------------------ - ! Collect all emissions of the same category (and species) on - ! the hierarchy level into array HirFlx. HirMsk contains the - ! combined covered region. That is, if there are two regional - ! inventories with the same hierarchy HirMsk will cover both - ! of these regions. - ! The specified field hierarchies determine whether the - ! temporary emissions are added (if hierarchy is the same - ! as the previously used hierarchy), or if they overwrite the - ! previous values in HirFlx (if hierarchy is higher than the - ! previous hierarchy). - ! ------------------------------------------------------------ - - ! Add emissions to the hierarchy array HirFlx if this hierarchy - ! is the same as previous hierarchy - IF ( ThisHir == PrevHir ) THEN - HirFlx = HirFlx + TmpFlx - HirMsk = HirMsk + Mask - - ! Make sure mask values do not exceed 1.0 - WHERE(HirMsk > 1.0 ) HirMsk = 1.0 - - ! If hierarchy is larger than those of the previously used - ! fields, overwrite HirFlx with new values. - ELSE - - HirFlx = TmpFlx - HirMsk = Mask - - ENDIF - - ! Update diagnostics at hierarchy level. Make sure that only - ! positive values are used. - ! The same diagnostics may be updated multiple times - ! during the same time step, continuously adding - ! emissions to it. - ! Now remove PosOnly flag. TmpFlx is initialized to zero, so it's - ! ok to keep negative values (ckeller, 7/12/15). - IF ( Diagn_AutoFillLevelDefined(HcoState%Diagn,4) .AND. DoDiagn ) THEN - CALL Diagn_Update( HcoState, ExtNr=ExtNr, & - Cat=ThisCat,Hier=ThisHir, HcoID=ThisSpc, & - !AutoFill=1, Array3D=TmpFlx, PosOnly=.TRUE.,& - AutoFill=1, Array3D=TmpFlx, & - COL=-1, RC=RC ) - IF ( RC /= HCO_SUCCESS ) RETURN - ENDIF - - ! Update previously used species, category and hierarchy - PrevSpc = ThisSpc - PrevCat = ThisCat - PrevHir = ThisHir - - ! Advance to next emission container - CALL ListCont_NextCont( HcoState%EmisList, Lct, FLAG ) - - ENDDO ! Loop over EmisList - - ! Make sure internal pointers are nullified - Lct => NULL() - Dct => NULL() - - ! Leave w/ success - CALL HCO_LEAVE ( HcoState%Config%Err, RC ) - - END SUBROUTINE CALC_EMS_SF_ADJ -#endif !EOC END MODULE HCO_Interface_Common