Skip to content

Commit

Permalink
Merge branch 'dev' into bugfix/bmy/interp
Browse files Browse the repository at this point in the history
  • Loading branch information
@yantosca
yantosca committed Apr 29, 2021
2 parents d224b1e + 3eeaf58 commit 00148d0
Show file tree
Hide file tree
Showing 11 changed files with 1,139 additions and 194 deletions.
894 changes: 770 additions & 124 deletions src/Core/hco_calc_mod.F90
View file

Large diffs are not rendered by default.

106 changes: 106 additions & 0 deletions src/Core/hco_diagn_mod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -430,6 +430,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
Expand Down Expand Up @@ -581,6 +646,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
! ------------------------------------------------------------------
Expand All @@ -596,6 +681,9 @@ SUBROUTINE Diagn_DefineFromConfig( HcoState, RC )
AutoFill = 1, &
COL = HcoState%Diagn%HcoDiagnIDDefault, &
RC = RC )
#ifdef ADJOINT
endif
#endif
IF ( RC /= HCO_SUCCESS ) RETURN

ENDDO
Expand Down Expand Up @@ -964,6 +1052,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
Expand Down Expand Up @@ -1780,6 +1874,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

Expand Down Expand Up @@ -4556,6 +4659,9 @@ SUBROUTINE DiagnBundle_Init ( Diagn )
Diagn%HcoDiagnIDDefault = -999
Diagn%HcoDiagnIDRestart = -999
Diagn%HcoDiagnIDManual = -999
#ifdef ADJOINT
Diagn%HcoDiagnIDAdjoint = -999
#endif
Diagn%nnCollections = 0
ENDIF

Expand Down
38 changes: 26 additions & 12 deletions src/Core/hco_geotools_mod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -854,8 +854,15 @@ SUBROUTINE HCO_CalcVertGrid ( HcoState, PSFC, ZSFC, TK, BXHEIGHT, PEDGE, RC )

ELSEIF ( EVAL_TK ) THEN
ALLOCATE(TmpTK(HcoState%NX,HcoState%NY,HcoState%NZ))
CALL HCO_EvalFld ( HcoState, 'TK', TmpTK, RC, FOUND=FoundTK )
CALL HCO_EvalFld( HcoState, 'TK', TmpTK, RC, FOUND=FoundTK )
IF ( RC /= HCO_SUCCESS ) RETURN

! TK is sometimes listed as TMPU, so look for that too (bmy, 3/5/21)
IF ( .not. FoundTK ) THEN
CALL HCO_EvalFld( HcoState, 'TMPU', TmpTK, RC, FOUND=FoundTK )
IF ( RC /= HCO_SUCCESS ) RETURN
ENDIF

EVAL_TK = FoundTK
IF ( FoundTK ) ThisTK => TmpTk

Expand Down Expand Up @@ -901,8 +908,17 @@ SUBROUTINE HCO_CalcVertGrid ( HcoState, PSFC, ZSFC, TK, BXHEIGHT, PEDGE, RC )

! Otherwise, try to read from HEMCO configuration file
ELSEIF ( EVAL_PSFC ) THEN
CALL HCO_EvalFld ( HcoState, 'PSFC', HcoState%Grid%PSFC%Val, RC, FOUND=FoundPSFC )
CALL HCO_EvalFld( HcoState, 'PSFC', HcoState%Grid%PSFC%Val, RC, &
FOUND=FoundPSFC )
IF ( RC /= HCO_SUCCESS ) RETURN

! PSFC is sometimes listed as PS, so look for that too (bmy, 3/4/21)
IF ( .not. FoundPSFC ) THEN
CALL HCO_EvalFld( HcoState, 'PS', HcoState%Grid%PSFC%Val, RC, &
FOUND=FoundPSFC )
IF ( RC /= HCO_SUCCESS ) RETURN
ENDIF

EVAL_PSFC = FoundPSFC

! Verbose
Expand Down Expand Up @@ -1114,10 +1130,9 @@ SUBROUTINE HCO_CalcVertGrid ( HcoState, PSFC, ZSFC, TK, BXHEIGHT, PEDGE, RC )

! Set PEDGE
IF ( .NOT. FoundPEDGE ) THEN
!$OMP PARALLEL DO &
!$OMP DEFAULT( SHARED ) &
!$OMP PRIVATE( I, J, L ) &
!$OMP SCHEDULE( DYNAMIC )
!$OMP PARALLEL DO &
!$OMP DEFAULT( SHARED ) &
!$OMP PRIVATE( I, J, L )
DO L = 1, HcoState%NZ+1
DO J = 1, HcoState%NY
DO I = 1, HcoState%NX
Expand All @@ -1128,7 +1143,7 @@ SUBROUTINE HCO_CalcVertGrid ( HcoState, PSFC, ZSFC, TK, BXHEIGHT, PEDGE, RC )
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
!$OMP END PARALLEL DO
FoundPEDGE = .TRUE.

! Verbose
Expand All @@ -1152,10 +1167,9 @@ SUBROUTINE HCO_CalcVertGrid ( HcoState, PSFC, ZSFC, TK, BXHEIGHT, PEDGE, RC )
HcoState%NY, HcoState%NZ, RC )
IF ( RC /= HCO_SUCCESS ) RETURN

!$OMP PARALLEL DO &
!$OMP DEFAULT( SHARED ) &
!$OMP PRIVATE( I, J, L, P1, P2 ) &
!$OMP SCHEDULE( DYNAMIC )
!$OMP PARALLEL DO &
!$OMP DEFAULT( SHARED ) &
!$OMP PRIVATE( I, J, L, P1, P2 )
DO L = 1, HcoState%NZ
DO J = 1, HcoState%NY
DO I = 1, HcoState%NX
Expand Down Expand Up @@ -1188,7 +1202,7 @@ SUBROUTINE HCO_CalcVertGrid ( HcoState, PSFC, ZSFC, TK, BXHEIGHT, PEDGE, RC )
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO
!$OMP END PARALLEL DO

IF ( ERRZSFC ) THEN
MSG = 'Cannot calculate surface geopotential heights - at least one ' // &
Expand Down
117 changes: 116 additions & 1 deletion src/Core/hco_interp_mod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -97,6 +97,27 @@ MODULE HCO_Interp_Mod
0.066000_hp, 0.047585_hp, 0.032700_hp, &
0.020000_hp, 0.010000_hp /)

! AP parameter of native 102-layer GISS grid
REAL(hp), TARGET :: E102_EDGE_NATIVE(103) = (/ &
0.0000000, 2.7871507, 5.5743014, 8.3614521, 11.1486028, 13.9357536, &
16.7229043, 19.5100550, 22.2972057, 25.0843564, 27.8715071, 30.6586578, &
33.4458085, 36.2329593, 39.0201100, 41.8087123, 44.6089278, 47.4534183, &
50.4082336, 53.5662786, 57.0095710, 60.7533531, 64.7323011, 68.8549615, &
73.0567364, 77.2969797, 81.5364973, 85.7346430, 89.8565776, 93.8754457, &
97.7709243, 101.5277712, 105.1350991, 108.5878272, 111.8859556, 115.0302100, &
118.0249453, 120.8854039, 123.6326345, 126.2811535, 128.8360417, 131.2987506, &
133.6736353, 135.9708571, 138.2013035, 140.3700552, 142.4814670, 144.5457005, &
146.5692881, 148.5464231, 150.4712991, 152.3497225, 154.1875000, 144.5468750, &
135.1875000, 126.0781250, 117.1914062, 108.5859375, 100.3671875, 92.5898438, &
85.2265625, 78.2226562, 71.5546875, 65.2226562, 59.2226562, 53.5546875, &
48.2226562, 43.2226562, 38.5546875, 34.2226562, 30.2226562, 26.5507812, &
23.1875000, 20.0781250, 17.1896562, 14.5684375, 12.2865742, 10.3573086, &
8.7353750, 7.3664922, 6.2100156, 5.2343633, 4.4119297, 3.7186797, &
3.1341479, 2.6404328, 2.2207877, 1.8587369, 1.5477125, 1.2782115, &
1.0427319, 0.8367716, 0.6514691, 0.4772511, 0.3168814, 0.1785988, &
0.1000000, 0.0560000, 0.0320000, 0.0180000, 0.0100000, 0.0050000, &
0.0020000 /)

CONTAINS
!EOC
!------------------------------------------------------------------------------
Expand Down Expand Up @@ -520,6 +541,26 @@ SUBROUTINE ModelLev_Check( HcoState, nLev, IsModelLev, RC )
nlev <= 47 ) THEN
IsModelLev = .TRUE.
ENDIF

! Full GISS 102-layer grid
ELSEIF ( nz == 102 ) THEN
IF ( nlev <= 103 ) THEN
IsModelLev = .TRUE.
ENDIF

! Full GISS 40-layer grid
ELSEIF ( nz == 40 ) THEN
IF ( nlev <= 41 ) THEN
IsModelLev = .TRUE.
ENDIF

! Reduced GISS 74-layer grid
ELSEIF ( nz == 74 ) THEN
IF ( nlev == 102 .OR. &
nlev == 103 .OR. &
nlev <= 74 ) THEN
IsModelLev = .TRUE.
ENDIF
ENDIF

END SUBROUTINE ModelLev_Check
Expand Down Expand Up @@ -821,6 +862,78 @@ SUBROUTINE ModelLev_Interpolate( HcoState, REGR_4D, Lct, RC )

! Done!
DONE = .TRUE.

!----------------------------------------------------------------
! Reduced GISS levels
!----------------------------------------------------------------
ELSEIF ( nz == 74 ) THEN

! Determine number of output levels. If input data is on the
! native grid, we collapse them onto the reduced GISS grid.
! In all other cases, we assume the input data is already on
! the reduced levels and mappings occurs 1:1.
IF ( nlev == 102 ) THEN
nout = nz
NL = 60
ELSEIF ( nlev == 103 ) THEN
nz = nz + 1
nout = nz
NL = 61
ELSE
nout = nlev
NL = nout
ENDIF

! Make sure output array is allocated
CALL FileData_ArrCheck( HcoState%Config, Lct%Dct%Dta, nx, ny, nout, nt, RC )

! Do for every time slice
DO T = 1, nt

! Levels that are passed level-by-level.
DO L = 1, NL
Lct%Dct%Dta%V3(T)%Val(:,:,L) = REGR_4D(:,:,L,T)
ENDDO !L

! If needed, collapse from native GEOS-5 onto reduced GEOS-5
IF ( nlev == 102 .OR. nlev == 103 ) THEN

! Add one level offset if these are edges
IF ( nlev == 103 ) THEN
OS = 1
ELSE
OS = 0
ENDIF

! Collapse two levels (e.g. levels 61-62 into level 61):
CALL COLLAPSE( Lct, REGR_4D, 61+OS, 61+OS, 2, T, 22 )
CALL COLLAPSE( Lct, REGR_4D, 62+OS, 63+OS, 2, T, 22 )
CALL COLLAPSE( Lct, REGR_4D, 63+OS, 65+OS, 2, T, 22 )
CALL COLLAPSE( Lct, REGR_4D, 64+OS, 67+OS, 2, T, 22 )
CALL COLLAPSE( Lct, REGR_4D, 65+OS, 69+OS, 2, T, 22 )
CALL COLLAPSE( Lct, REGR_4D, 66+OS, 71+OS, 2, T, 22 )
CALL COLLAPSE( Lct, REGR_4D, 67+OS, 73+OS, 2, T, 22 )
! Collapse four levels:
CALL COLLAPSE( Lct, REGR_4D, 68+OS, 75+OS, 4, T, 22 )
CALL COLLAPSE( Lct, REGR_4D, 69+OS, 79+OS, 4, T, 22 )
CALL COLLAPSE( Lct, REGR_4D, 70+OS, 83+OS, 4, T, 22 )
CALL COLLAPSE( Lct, REGR_4D, 71+OS, 87+OS, 4, T, 22 )
CALL COLLAPSE( Lct, REGR_4D, 72+OS, 91+OS, 4, T, 22 )
CALL COLLAPSE( Lct, REGR_4D, 73+OS, 95+OS, 4, T, 22 )
CALL COLLAPSE( Lct, REGR_4D, 74+OS, 99+OS, 4, T, 22 )

ENDIF
ENDDO ! T

! Verbose
IF ( HCO_IsVerb(HcoState%Config%Err, 3) ) THEN
WRITE(MSG,*) 'Mapped ', nlev, ' levels onto reduced GISS levels.'
CALL HCO_MSG(HcoState%Config%Err,MSG)
ENDIF

! Done!
DONE = .TRUE.

ENDIF

ENDIF ! Vertical regridding required
Expand Down Expand Up @@ -1117,7 +1230,7 @@ SUBROUTINE COLLAPSE ( Lct, REGR_4D, OutLev, InLev1, NLEV, T, MET )
INTEGER, INTENT(IN) :: InLev1
INTEGER, INTENT(IN) :: NLEV
INTEGER, INTENT(IN) :: T
INTEGER, INTENT(IN) :: MET ! 4=GEOS-4, else GEOS-5
INTEGER, INTENT(IN) :: MET ! 4=GEOS-4, 22=GISS E2.2, else GEOS-5
!
! !INPUT/OUTPUT PARAMETERS:
!
Expand Down Expand Up @@ -1154,6 +1267,8 @@ SUBROUTINE COLLAPSE ( Lct, REGR_4D, OutLev, InLev1, NLEV, T, MET )
! Get pointer to grid edges on the native input grid
IF ( Met == 4 ) THEN
EDG => G4_EDGE_NATIVE(InLev1:TOPLEV)
ELSE IF ( Met == 22 ) THEN
EDG => E102_EDGE_NATIVE(InLev1:TOPLEV)
ELSE
EDG => G5_EDGE_NATIVE(InLev1:TOPLEV)
ENDIF
Expand Down
Loading

0 comments on commit 00148d0

Please sign in to comment.