Skip to content

Commit

Permalink
Add fix for the lowest level in CH4 GOSAT observation operator
Browse files Browse the repository at this point in the history 
Yuzhong Zhang wrote:

   GOSAT record has negative pressure value when the surface has high
   altitude such as Tibet. Fix it by finding the lowest valid layer L0.

Signed-off-by: Melissa Sulprizio <[email protected]>
  • Loading branch information
@msulprizio
msulprizio committed Sep 23, 2019
1 parent b2fb868 commit f397e1f
Showing 1 changed file with 47 additions and 20 deletions.
67 changes: 47 additions & 20 deletions GeosCore/gosat_ch4_mod.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -445,6 +445,11 @@ SUBROUTINE CALC_GOSAT_CH4_FORCE( Input_Opt, State_Chm, State_Grid,
REAL(fp) :: DIFF
REAL(fp) :: S_OBS

! --- zyz --- Sept 19, 2018
! fix the problem that some GOSAT record has negative pressure
! values if the surface has high altidue
INTEGER :: L0

!------------------------------------------------------------------------------
! For adjoint model only:
! ! For adjoint of observation operator
Expand Down Expand Up @@ -662,6 +667,22 @@ SUBROUTINE CALC_GOSAT_CH4_FORCE( Input_Opt, State_Chm, State_Grid,
CYCLE
ENDIF

! zyz, Oct 4 2018
! GOSAT record has negative pressure value when the surface
! has high altitude such as Tibet
! Fix it by finding the lowest valid layer L0
L0 = GOS(NT)%LGOS(1)
DO L = 1, GOS(NT)%LGOS(1)
IF (GOS(NT)%PRES(L) .gt. 0.0_fp) THEN
L0 = L
EXIT
ENDIF
ENDDO
! Only find negative pressure at 1st layer so far.
! But use this condition to prevent really invalid cases
! (zyz)
IF (L0 .gt. 3 ) CYCLE

! Skip observations outside the domain
IF ( GOS(NT)%LAT(1) < State_Grid%YMin .OR.
& GOS(NT)%LAT(1) > State_Grid%YMax .OR.
Expand All @@ -688,7 +709,8 @@ SUBROUTINE CALC_GOSAT_CH4_FORCE( Input_Opt, State_Chm, State_Grid,

! skip observations where the GOSAT surface pressure is
! less than the model
IF ( (GOS(NT)%PRES(1) - State_Met%PEDGE(I,J,1)) > 50e0 ) THEN
! Use L0 for lowest valid layer in GOSAT (zyz)
IF ( (GOS(NT)%PRES(L0) - State_Met%PEDGE(I,J,1)) > 50e0 ) THEN
print*, ' Psurf threshold not met, skipping record ', NT
CYCLE
ENDIF
Expand All @@ -706,8 +728,9 @@ SUBROUTINE CALC_GOSAT_CH4_FORCE( Input_Opt, State_Chm, State_Grid,
CH4_PRIOR(:) = 0.0_fp

! Copy variable names to make coding a bit cleaner
! Use L0 for lowest layer in GOSAT record (zyz)
LGOS = GOS(NT)%LGOS(1)
DO L = 1, LGOS
DO L = L0, LGOS
CH4_PRIOR(L) = GOS(NT)%PRIOR(L) * 1d-9 ! [ppb] --> [v/v]
ENDDO

Expand All @@ -727,7 +750,7 @@ SUBROUTINE CALC_GOSAT_CH4_FORCE( Input_Opt, State_Chm, State_Grid,
! Calculate the interpolation weight matrix
MAP(:,:) = 0.0_fp
CALL GET_INTMAP( State_Grid, GC_PRES, GC_PSURF, GOS(NT)%PRES,
& LGOS, MAP )
& L0, LGOS, MAP )

!------------------------------------------------------------------------------
! For adjoint model only:
Expand Down Expand Up @@ -762,7 +785,8 @@ SUBROUTINE CALC_GOSAT_CH4_FORCE( Input_Opt, State_Chm, State_Grid,


! Interpolate GC CH4 column to GOSAT grid
DO LL = 1, LGOS
! Use L0 for lowest valid layer for GOSAT (zyz)
DO LL = L0, LGOS
GC_CH4(LL) = 0.0_fp
DO L = 1, State_Grid%NZ
GC_CH4(LL) = GC_CH4(LL)
Expand Down Expand Up @@ -823,48 +847,49 @@ SUBROUTINE CALC_GOSAT_CH4_FORCE( Input_Opt, State_Chm, State_Grid,
!--------------------------------------------------------------

! Pressure weighting array
! Use L0 for lowest valid layer in GOSAT
h(:) = 0.0_fp
p(1:LGOS) = GOS(NT)%PRES(1:LGOS)
ak(1:LGOS) = GOS(NT)%AVG_KERNEL(1:LGOS)
prior(1:LGOS) = GOS(NT)%PRIOR(1:LGOS) * 1e-9_fp ! [ppb] --> [v/v]
p(L0:LGOS) = GOS(NT)%PRES(L0:LGOS)
ak(L0:LGOS) = GOS(NT)%AVG_KERNEL(L0:LGOS)
prior(L0:LGOS) = GOS(NT)%PRIOR(L0:LGOS) * 1e-9_fp ! [ppb] --> [v/v]

! Need to integrate from the toa to surface (ajt, 05/21/13)
IF (LGOS .GT. 1) THEN
IF(GOS(NT)%PRES(2) .LT. GOS(NT)%PRES(1)) THEN
IF(GOS(NT)%PRES(L0+1) .LT. GOS(NT)%PRES(L0)) THEN
p(1:LGOS) = p(LGOS:1:-1)
ENDIF
ENDIF

L = 1
h(L) = 1./GOS(NT)%PRES(1) * ABS(
h(L) = 1./GOS(NT)%PRES(L0) * ABS(
& ( -1e0*p(L) + ( (p(L+1)-p(L))/(LOG(p(L+1)/p(L))) ) ) )
L = LGOS
h(L) = 1./GOS(NT)%PRES(1) * ABS(
L = LGOS - L0 + 1
h(L) = 1./GOS(NT)%PRES(L0) * ABS(
& ( p(L) - ( (p(L)-p(L-1))/(LOG(p(L)/p(L-1))) ) ) )
DO L=2,LGOS-1
h(L) = 1./GOS(NT)%PRES(1) * ABS(
DO L=2,LGOS-L0
h(L) = 1./GOS(NT)%PRES(L0) * ABS(
& ( -1e0*p(L) + ( (p(L+1)-p(L))/(LOG(p(L+1)/p(L))) ) ) +
& ( p(L) - ( (p(L)-p(L-1))/(LOG(p(L)/p(L-1))) ) ) )
ENDDO

! Now return to the orientation of the other variables
IF (LGOS .GT. 1) THEN
IF(GOS(NT)%PRES(2) .LT. GOS(NT)%PRES(1)) THEN
IF(GOS(NT)%PRES(L0+1) .LT. GOS(NT)%PRES(L0)) THEN
h(1:LGOS) = h(LGOS:1:-1)
p(1:LGOS) = p(LGOS:1:-1)
ENDIF
ENDIF

! Calculate Xch4_a
XCH4a = 0.0_fp
DO L = 1,LGOS
DO L = L0,LGOS
XCH4a = XCH4a + h(L) * prior(L)
ENDDO

! Calculate Xch4_m
XCH4m = 0.0_fp
XCH4m = XCH4a
DO L = 1, LGOS
DO L = L0, LGOS
XCH4m = XCH4m + ( h(L) * ak(L) * ( GC_CH4(L) - prior(L) ) )
ENDDO
GC_XCH4 = 0.0_fp
Expand Down Expand Up @@ -1819,7 +1844,7 @@ END SUBROUTINE CALC_GOSAT_CH4_FORCE
! !INTERFACE:
!
SUBROUTINE GET_INTMAP( State_Grid, GCPCEN, GCPSURF, GOSPEDGE,
& nlev, INTMAP )
& L0, L1, INTMAP )
!
! !USES:
!
Expand All @@ -1831,7 +1856,7 @@ SUBROUTINE GET_INTMAP( State_Grid, GCPCEN, GCPSURF, GOSPEDGE,
REAL(fp), INTENT(IN) :: GCPCEN(State_Grid%NZ)
REAL(fp), INTENT(IN) :: GCPSURF
REAL(fp), INTENT(IN) :: GOSPEDGE(MAXLEV)
INTEGER, INTENT(IN) :: nlev
INTEGER, INTENT(IN) :: L0, L1
!
! !OUTPUT PARAMETERS:
!
Expand All @@ -1840,6 +1865,8 @@ SUBROUTINE GET_INTMAP( State_Grid, GCPCEN, GCPSURF, GOSPEDGE,
! !REVISION HISTORY:
! 16 Jun 2017 - M. Sulprizio- Initial version based on GOSAT CH4 observation
! operator from GC Adjoint v35j
! 04 Oct 2018 - Y. Zhang - Add an input argument L0, the lowest valid layer for GOSAT
!EOP
!------------------------------------------------------------------------------
!BOC
Expand All @@ -1858,7 +1885,7 @@ SUBROUTINE GET_INTMAP( State_Grid, GCPCEN, GCPSURF, GOSPEDGE,
INTMAP(:,:) = 0e+0_fp
! Loop over each pressure level of GOS grid
DO LTM = 1, nlev
DO LTM = L0, L1
! Find the levels from GC that bracket level LTM
DO LGC = 1, State_Grid%NZ-1
Expand All @@ -1882,7 +1909,7 @@ SUBROUTINE GET_INTMAP( State_Grid, GCPCEN, GCPSURF, GOSPEDGE,
! Correct for case where GOSAT pressure is higher than the
! highest GC pressure center. In this case, just 1:1 map.
DO LTM = 1, nlev
DO LTM = L0, L1
IF ( GOSPEDGE(LTM) > GCPCEN(1) ) THEN
INTMAP(:,LTM) = 0e+0_fp
INTMAP(1,LTM) = 1e+0_fp
Expand Down

0 comments on commit f397e1f

Please sign in to comment.