971 lines
36 KiB
Fortran
971 lines
36 KiB
Fortran
!$Id: IASI_o3_mod.f,v 1.3 2011/02/23 00:08:48 daven Exp $
|
|
MODULE IASI_CO_OBS_MOD
|
|
|
|
IMPLICIT NONE
|
|
|
|
!mkeller
|
|
#include "CMN_SIZE"
|
|
!#include 'netcdf.inc'
|
|
|
|
!=================================================================
|
|
! MODULE VARIABLES
|
|
!=================================================================
|
|
|
|
PRIVATE
|
|
|
|
PUBLIC READ_IASI_CO_OBS
|
|
PUBLIC CALC_IASI_CO_FORCE
|
|
|
|
! Parameters
|
|
INTEGER, PARAMETER :: N_IASI_NLA = 19
|
|
INTEGER, PARAMETER :: N_IASI_NPR = 20
|
|
INTEGER, PARAMETER :: MAXIASI = 2000000
|
|
INTEGER, PARAMETER :: IASI_COL = 60!59
|
|
|
|
! Module variables
|
|
|
|
! IASI data
|
|
REAL*8, ALLOCATABLE :: IASI_LON(:)
|
|
REAL*8, ALLOCATABLE :: IASI_LAT(:)
|
|
REAL*8, ALLOCATABLE :: IASI_TIME(:)
|
|
REAL*8, ALLOCATABLE :: IASI_CO(:)
|
|
REAL*8, ALLOCATABLE :: IASI_CO_STD(:)
|
|
REAL*8, ALLOCATABLE :: IASI_ALT(:)
|
|
REAL*8, ALLOCATABLE :: IASI_CO_APR(:,:)
|
|
REAL*8, ALLOCATABLE :: IASI_CO_AVK(:,:)
|
|
REAL*8, ALLOCATABLE :: IASI_QUAL_FLAG(:)
|
|
REAL*8, ALLOCATABLE :: IASI_CLOUD_COVER(:)
|
|
REAL*8, ALLOCATABLE :: IASI_SOLAR_ZENITH(:)
|
|
REAL*8, ALLOCATABLE :: IASI_CO_LVL(:)
|
|
!REAL*8, ALLOCATABLE :: TIME_FRAC(:)
|
|
!REAL*8, ALLOCATABLE :: TEMPDATA(:,:)
|
|
|
|
! MLS grid specification
|
|
INTEGER :: N_IASI_NOB, N_IASI_DATA
|
|
|
|
|
|
! mkeller: logical flag to check whether data is available for given day
|
|
LOGICAL :: DATA_PRESENT
|
|
CONTAINS
|
|
!------------------------------------------------------------------------------
|
|
|
|
SUBROUTINE READ_IASI_CO_OBS( YYYYMMDD, N_IASI_NOB )
|
|
!
|
|
!******************************************************************************
|
|
! Subroutine READ_IASI_CO_OBS reads the file and passes back info contained
|
|
! therein. (dkh, 02/19/09)
|
|
!
|
|
! Based on READ_IASI_NH3 OBS (dkh, 04/26/10)
|
|
!
|
|
! Arguments as Input:
|
|
! ============================================================================
|
|
! (1 ) YYYYMMDD INTEGER : Current year-month-day
|
|
!
|
|
! Arguments as Output:
|
|
! ============================================================================
|
|
! (1 ) NIASI (INTEGER) : Number of IASI retrievals for current day
|
|
!
|
|
! Module variable as Output:
|
|
! ============================================================================
|
|
! (1 ) IASI (IASI_CO_OBS) : IASI retrieval for current day
|
|
!
|
|
! NOIASI:
|
|
! (1 ) Add calculation of S_OER_INV, though eventually we probably want to
|
|
! do this offline. (dkh, 05/04/10)
|
|
!******************************************************************************
|
|
!
|
|
! Reference to f90 modules
|
|
USE DIRECTORY_MOD, ONLY : DATA_DIR
|
|
USE TIME_MOD, ONLY : EXPAND_DATE
|
|
USE FILE_MOD, ONLY : IOERROR
|
|
|
|
#include "CMN_SIZE" ! size parameters
|
|
|
|
|
|
! Arguments
|
|
INTEGER, INTENT(IN) :: YYYYMMDD
|
|
|
|
! local variables
|
|
INTEGER :: FID
|
|
INTEGER :: LIASI
|
|
INTEGER :: N_IASI_NOB
|
|
INTEGER :: N, J
|
|
|
|
CHARACTER(LEN=5) :: TMP
|
|
CHARACTER(LEN=255) :: READ_FILENAME
|
|
|
|
REAL*8, PARAMETER :: FILL = -999.0D0
|
|
REAL*8, PARAMETER :: TOL = 1d-04
|
|
REAL*8 :: TMP1
|
|
REAL*8 :: DUMMY(IASI_COL)
|
|
REAL*8 :: TEMPDATA(MAXIASI, IASI_COL)
|
|
INTEGER :: I, II, III, K
|
|
INTEGER :: IU_FILE, IU_DATA, IOS
|
|
|
|
|
|
!=================================================================
|
|
! READ_IASI_CO_OBS begins here!
|
|
!=================================================================
|
|
CALL CLEANUP_IASI
|
|
! filename root
|
|
READ_FILENAME = TRIM( 'iasi_CO_LATMOS_ULB_YYYYMMDD_v20140922.txt' )
|
|
|
|
! Expand date tokens in filename
|
|
CALL EXPAND_DATE( READ_FILENAME, YYYYMMDD, 9999 )
|
|
|
|
! Construct complete filename
|
|
! READ_FILENAME = TRIM( DATA_DIR ) // TRIM( 'IASI_CO/' ) //
|
|
! & TRIM( READ_FILENAME )
|
|
READ_FILENAME = '/users/jk/15/xzhang/IASI_CO/' // TRIM( READ_FILENAME )
|
|
WRITE(6,*) ' - READ_IASI_CO_OBS: reading file: ', READ_FILENAME
|
|
IU_FILE = 67
|
|
! mkeller: check to see if file exists
|
|
INQUIRE(FILE=READ_FILENAME, EXIST = DATA_PRESENT)
|
|
|
|
IF (.NOT. DATA_PRESENT) THEN
|
|
PRINT *,"IASI file", TRIM(READ_FILENAME), " not found, "// "assuming that there is no data for this day."
|
|
RETURN
|
|
ELSE
|
|
PRINT *,"IASI file found!"
|
|
ENDIF
|
|
|
|
OPEN(IU_FILE, FILE=READ_FILENAME, IOSTAT=IOS)
|
|
N_IASI_DATA = 0
|
|
N_IASI_NOB = 0
|
|
DO
|
|
READ(IU_FILE, *, IOSTAT=IOS) DUMMY
|
|
IF (IOS /= 0) EXIT
|
|
N_IASI_DATA = N_IASI_DATA + 1
|
|
I = I + 1
|
|
TEMPDATA(I,:) = DUMMY(:)
|
|
ENDDO
|
|
CLOSE(IU_FILE)
|
|
|
|
DO I = 1, N_IASI_DATA
|
|
IF (TEMPDATA(I,3) .EQ. REAL(YYYYMMDD)) THEN
|
|
N_IASI_NOB = N_IASI_NOB + 1
|
|
ENDIF
|
|
ENDDO
|
|
|
|
ALLOCATE(IASI_LAT(N_IASI_NOB))
|
|
ALLOCATE(IASI_LON(N_IASI_NOB))
|
|
ALLOCATE(IASI_TIME(N_IASI_NOB))
|
|
ALLOCATE(IASI_SOLAR_ZENITH(N_IASI_NOB))
|
|
ALLOCATE(IASI_QUAL_FLAG(N_IASI_NOB))
|
|
ALLOCATE(IASI_CLOUD_COVER(N_IASI_NOB))
|
|
ALLOCATE(IASI_CO(N_IASI_NOB))
|
|
ALLOCATE(IASI_CO_STD(N_IASI_NOB))
|
|
ALLOCATE(IASI_CO_APR(N_IASI_NOB, N_IASI_NLA))
|
|
ALLOCATE(IASI_CO_AVK(N_IASI_NOB, N_IASI_NLA))
|
|
ALLOCATE(IASI_ALT(N_IASI_NPR))
|
|
ALLOCATE(IASI_CO_LVL(N_IASI_NOB))
|
|
|
|
DO I = 1, N_IASI_NOB
|
|
IASI_LAT(I) = TEMPDATA(I,1)
|
|
IASI_LON(I) = TEMPDATA(I,2)
|
|
IASI_TIME(I) = TEMPDATA(I,4)
|
|
IASI_SOLAR_ZENITH(I) = TEMPDATA(I,5)
|
|
IASI_QUAL_FLAG(I) = TEMPDATA(I,16)!15)
|
|
!PRINT *, "IASI_QUAL_FLAG", IASI_QUAL_FLAG(I)
|
|
IASI_CLOUD_COVER(I) = TEMPDATA(I,17)!16)
|
|
!READ(IU_FILE, *) (TEMPDATA(K), K=19,IASI_COL)
|
|
IASI_CO(I) = TEMPDATA(I,21)!20)
|
|
IASI_CO_STD(I) = TEMPDATA(I,21)*TEMPDATA(I,22)*0.5
|
|
DO J = 1, N_IASI_NLA
|
|
IASI_CO_APR(I,J) = TEMPDATA(I,22+J) !21
|
|
IASI_CO_AVK(I,J) = TEMPDATA(I,41+J) !40
|
|
!IASI_ALT(J) = REAL(J)-1d0
|
|
ENDDO
|
|
!IASI_ALT(N_IASI_NLA+1) = 60d0
|
|
ENDDO
|
|
DO J = 1, N_IASI_NLA
|
|
IASI_ALT(J) = REAL(J) -1d0
|
|
ENDDO
|
|
|
|
IASI_ALT(N_IASI_NPR) = 60d0
|
|
!PRINT *, "TEMPDATA", TEMPDATA(:)
|
|
!PRINT *, "IASI_TIME", TEMPDATA(N_IASI_NOB-10000:N_IASI_NOB,3)
|
|
!PRINT *, "IASI_CO_APR", IASI_CO_APR(1,:), "LEVEL2", IASI_CO_APR(2,:)
|
|
!PRINT *, "IASI_CO_AVK", IASI_CO_AVK(1,:), "LEVEL2", IASI_CO_AVK(2,:)
|
|
!PRINT *, "IASI_CO", IASI_CO(1:10)
|
|
!--------------------------------
|
|
! Calculate S_OER_INV
|
|
!--------------------------------
|
|
|
|
! loop over records
|
|
! Now determine how many of the levels in CO are
|
|
! 'good' and how many are just FILL.
|
|
!ALLOCATE(IASI_CO_LVL(N_IASI_NOB))
|
|
DO N = 1, N_IASI_NOB
|
|
J = 1
|
|
DO WHILE ( J .le. N_IASI_NLA )
|
|
! check if the value is good
|
|
IF (IASI_CO_APR(N,J) > FILL ) THEN
|
|
! save the number of good levels as LTES
|
|
IASI_CO_LVL(N) = N_IASI_NLA - J + 1
|
|
! and now we can exit the while loop
|
|
J = N_IASI_NLA + 1
|
|
! otherwise this level is just filler
|
|
ELSE
|
|
! so proceed to the next one up
|
|
J = J + 1
|
|
ENDIF
|
|
ENDDO
|
|
ENDDO
|
|
|
|
! Return to calling program
|
|
END SUBROUTINE READ_IASI_CO_OBS
|
|
|
|
!------------------------------------------------------------------------------
|
|
|
|
SUBROUTINE CHECK( STATUS, LOCATION )
|
|
!
|
|
!******************************************************************************
|
|
! Subroutine CHECK checks the status of calls to netCDF libraries routines
|
|
! (dkh, 02/15/09)
|
|
!
|
|
! Arguments as Input:
|
|
! ============================================================================
|
|
! (1 ) STATUS (INTEGER) : Completion status of netCDF library call
|
|
! (2 ) LOCATION (INTEGER) : Location at which netCDF library call was made
|
|
!
|
|
!
|
|
! NOIASI:
|
|
!
|
|
!******************************************************************************
|
|
!
|
|
! Reference to f90 modules
|
|
USE ERROR_MOD, ONLY : ERROR_STOP
|
|
USE NETCDF
|
|
|
|
! Arguments
|
|
INTEGER, INTENT(IN) :: STATUS
|
|
INTEGER, INTENT(IN) :: LOCATION
|
|
|
|
!=================================================================
|
|
! CHECK begins here!
|
|
!=================================================================
|
|
|
|
IF ( STATUS /= NF90_NOERR ) THEN
|
|
WRITE(6,*) TRIM( NF90_STRERROR( STATUS ) )
|
|
WRITE(6,*) 'At location = ', LOCATION
|
|
CALL ERROR_STOP('netCDF error', 'IASI_o3_mod')
|
|
ENDIF
|
|
|
|
! Return to calling program
|
|
END SUBROUTINE CHECK
|
|
|
|
!------------------------------------------------------------------------------
|
|
|
|
SUBROUTINE CALC_IASI_CO_FORCE( COST_FUNC )
|
|
!
|
|
!******************************************************************************
|
|
! Subroutine CALC_IASI_CO_FORCE calculaIASI the adjoint forcing from the IASI
|
|
! CO observations and updates the cost function. (dkh, 02/15/09)
|
|
!
|
|
!
|
|
! Arguments as Input/Output:
|
|
! ============================================================================
|
|
! (1 ) COST_FUNC (REAL*8) : Cost function [unitless]
|
|
!
|
|
!
|
|
! NOIASI:
|
|
! (1 ) Updated to GCv8 (dkh, 10/07/09)
|
|
! (2 ) Add more diagnostics. Now read and write doubled CO (dkh, 11/08/09)
|
|
! (3 ) Now use CSPEC_AFTER_CHEM and CSPEC_AFTER_CHEM_ADJ (dkh, 02/09/11)
|
|
!******************************************************************************
|
|
!
|
|
! Reference to f90 modules
|
|
USE ADJ_ARRAYS_MOD, ONLY : STT_ADJ, ADJ_FORCE
|
|
USE ADJ_ARRAYS_MOD, ONLY : N_CALC
|
|
USE ADJ_ARRAYS_MOD, ONLY : EXPAND_NAME
|
|
!USE ADJ_ARRAYS_MOD, ONLY : CO_PROF_SAV
|
|
USE ADJ_ARRAYS_MOD, ONLY : ID2C
|
|
USE CHECKPT_MOD, ONLY : CHK_STT
|
|
USE COMODE_MOD, ONLY : JLOP, VOLUME
|
|
USE COMODE_MOD, ONLY : CSPEC_AFTER_CHEM
|
|
USE COMODE_MOD, ONLY : CSPEC_AFTER_CHEM_ADJ
|
|
USE DAO_MOD, ONLY : AD, AIRDEN, AIRVOL
|
|
USE DAO_MOD, ONLY : BXHEIGHT
|
|
USE DIRECTORY_ADJ_MOD, ONLY : DIAGADJ_DIR
|
|
USE GRID_MOD, ONLY : GET_IJ, GET_XMID, GET_YMID
|
|
USE PRESSURE_MOD, ONLY : GET_PCENTER, GET_PEDGE
|
|
USE TIME_MOD, ONLY : GET_NYMD, GET_NHMS
|
|
USE TIME_MOD, ONLY : GET_TS_CHEM, GET_HOUR
|
|
USE TRACER_MOD, ONLY : XNUMOLAIR, XNUMOL
|
|
USE TRACER_MOD, ONLY : TCVV
|
|
USE TRACERID_MOD, ONLY : IDTCO
|
|
USE TROPOPAUSE_MOD, ONLY : ITS_IN_THE_TROP
|
|
|
|
|
|
# include "CMN_SIZE" ! Size params
|
|
|
|
! Arguments
|
|
REAL*8, INTENT(INOUT) :: COST_FUNC
|
|
|
|
! Local variables
|
|
INTEGER :: NTSTART, NTSTOP, NT, I_IASI
|
|
INTEGER :: IIJJ(2), I, J
|
|
INTEGER :: L, LL, LIASI, L0
|
|
INTEGER :: JLOOP
|
|
INTEGER :: GC_HOUR
|
|
REAL*8 :: GC_PRES(LLPAR)
|
|
REAL*8 :: GC_CO_NATIVE(LLPAR)
|
|
REAL*8 :: GC_CO(N_IASI_NLA), IASI_APR_WEB(N_IASI_NLA)
|
|
REAL*8 :: GC_PSURF, IASI_PSURF, GC_ALT(IIPAR, JJPAR, LLPAR+1)
|
|
REAL*8 :: CO_HAT(N_IASI_NLA)
|
|
REAL*8 :: SOBS_CO_AVK(IIPAR,JJPAR,N_IASI_NLA)
|
|
REAL*8 :: CO_PERT(N_IASI_NLA), SOBS_AVK_TOT(IIPAR,JJPAR,N_IASI_NLA)
|
|
REAL*8 :: FORCE, IASI_RATIO(N_IASI_NLA)
|
|
REAL*8 :: DIFF, CO_COL_OBS, CO_COL_HAT
|
|
REAL*8 :: IASI_PCENTER(N_IASI_NLA)
|
|
REAL*8 :: IASI_APR(N_IASI_NLA)
|
|
REAL*8 :: NEW_COST(IIPAR,JJPAR)
|
|
REAL*8 :: OLD_COST
|
|
REAL*8 :: XNUAIR, SOBS_RAND
|
|
REAL*8, SAVE :: TIME_FRAC(MAXIASI)
|
|
REAL*8 :: ALT_SURF(IIPAR,JJPAR)
|
|
REAL*8 :: SOBS_CO_STD, SOBS_CO_COL
|
|
REAL*8 :: SOBS_CO_NORM(IIPAR,JJPAR), SOBS_CO_MEAN(IIPAR,JJPAR)
|
|
REAL*8 :: SOBS_CO_TOT(IIPAR,JJPAR), SOBS_CO_AVG(IIPAR,JJPAR)
|
|
REAL*8 :: SOBS_CO(IIPAR,JJPAR), SOBS_VAR
|
|
REAL*8 :: SOBS_CO_ERR(IIPAR,JJPAR), SOBS_CO_VAR(IIPAR,JJPAR)
|
|
REAL*8 :: SOBS_ERR, SOBS_VAR_LIMIT, SOBS_CO_HAT(IIPAR,JJPAR)
|
|
REAL*8 :: SOBS_HAT(IIPAR,JJPAR)
|
|
|
|
|
|
REAL*8 :: GC_CO_NATIVE_ADJ(LLPAR)
|
|
REAL*8 :: CO_HAT_ADJ
|
|
REAL*8 :: CO_PERT_ADJ(N_IASI_NLA)
|
|
REAL*8 :: GC_CO_ADJ(N_IASI_NLA)
|
|
|
|
REAL*8 :: GC_ADJ_TEMP(IIPAR,JJPAR,LLPAR)
|
|
REAL*8 :: GC_ADJ_TEMP_COST(IIPAR,JJPAR)
|
|
! arrays needed for superobservations
|
|
LOGICAL :: SUPER_OBS = .TRUE. ! do super observations?
|
|
REAL*8 :: SOBS_COUNT(IIPAR,JJPAR)
|
|
|
|
LOGICAL, SAVE :: FIRST = .TRUE.
|
|
LOGICAL, SAVE :: SECOND = .TRUE.
|
|
INTEGER :: IU_FILE, IU_DATA, IOS
|
|
CHARACTER(LEN=255) :: FILENAME_ALT, ALT_PATH, FILE_ALT, FILENAME
|
|
|
|
|
|
!=================================================================
|
|
! CALC_IASI_CO_FORCE begins here!
|
|
!=================================================================
|
|
|
|
print*, ' - CALC_IASI_CO_FORCE '
|
|
CALL RAND
|
|
! Reset
|
|
IASI_APR(1:N_IASI_NLA) = (/10.161773,9.485584,9.1867937,8.9939589,8.8739,8.8418514,8.7545818, &
|
|
8.6420669,8.3986495,8.0001663,7.5277656,6.9161481,6.3326457,5.7343264, &
|
|
5.1686031,4.6027417,4.1054126,3.6674835,12.211041/)
|
|
XNUAIR = 28.964d-3
|
|
NEW_COST = 0D0
|
|
SOBS_COUNT = 0d0
|
|
ADJ_FORCE = 0d0
|
|
IASI_RATIO = 0d0
|
|
GC_ADJ_TEMP_COST = 0d0
|
|
SOBS_CO_NORM = 0d0
|
|
SOBS_CO_MEAN = 0d0
|
|
SOBS_CO_TOT = 0d0
|
|
SOBS_CO_AVG = 0d0
|
|
SOBS_CO = 0d0
|
|
SOBS_CO_ERR = 0d0
|
|
SOBS_CO_VAR = 0d0
|
|
|
|
GC_HOUR = GET_HOUR()
|
|
SOBS_VAR_LIMIT = 2.5e17
|
|
! Save a value of the cost function first
|
|
OLD_COST = COST_FUNC
|
|
ALT_PATH = '/users/jk/07/xzhang/met_field/'
|
|
FILE_ALT = '20000101.cn.4x5.dat'
|
|
|
|
FILENAME_ALT = TRIM(ALT_PATH) // TRIM(FILE_ALT)
|
|
OPEN(UNIT = 13, FILE = "/users/jk/07/xzhang/met_field/20000101.cn.4x5.dat", STATUS="old",ACTION="read")
|
|
READ(13,*) ALT_SURF
|
|
CLOSE(13)
|
|
PRINT *, "GET_NHMS", GET_NHMS()
|
|
IF ( SECOND ) THEN
|
|
FILENAME = 'lat_orb_iasico.NN.m'
|
|
CALL EXPAND_NAME( FILENAME, N_CALC )
|
|
FILENAME = TRIM( DIAGADJ_DIR ) // TRIM( FILENAME )
|
|
OPEN( 901, FILE=TRIM( FILENAME ), STATUS='UNKNOWN', IOSTAT=IOS, FORM='FORMATTED', ACCESS='SEQUENTIAL' )
|
|
|
|
FILENAME = 'lon_orb_iasico.NN.m'
|
|
CALL EXPAND_NAME( FILENAME, N_CALC )
|
|
FILENAME = TRIM( DIAGADJ_DIR ) // TRIM( FILENAME )
|
|
OPEN( 902, FILE=TRIM( FILENAME ), STATUS='UNKNOWN', IOSTAT=IOS, FORM='FORMATTED', ACCESS='SEQUENTIAL' )
|
|
|
|
FILENAME = 'co_chi_sq_iasico.NN.m'
|
|
CALL EXPAND_NAME( FILENAME, N_CALC )
|
|
FILENAME = TRIM( DIAGADJ_DIR ) // TRIM( FILENAME )
|
|
OPEN( 904, FILE=TRIM( FILENAME ), STATUS='UNKNOWN', IOSTAT=IOS, FORM='FORMATTED', ACCESS='SEQUENTIAL' )
|
|
|
|
FILENAME = 'diff_iasico.NN.m'
|
|
CALL EXPAND_NAME( FILENAME, N_CALC )
|
|
FILENAME = TRIM( DIAGADJ_DIR ) // TRIM( FILENAME )
|
|
OPEN( 912, FILE=TRIM( FILENAME ), STATUS='UNKNOWN', IOSTAT=IOS, FORM='FORMATTED', ACCESS='SEQUENTIAL' )
|
|
|
|
ENDIF
|
|
|
|
! Check if it is the last hour of a day
|
|
IF ( GET_NHMS() == 236000 - GET_TS_CHEM() * 100 ) THEN
|
|
|
|
! Read the IASI CO file for this day
|
|
CALL READ_IASI_CO_OBS( GET_NYMD(), N_IASI_NOB )
|
|
!PRINT *, "N_IASI", N_IASI_NOB
|
|
!PRINT *, "TIME", IASI_TIME(N_IASI_NOB-1000:N_IASI_NOB)
|
|
! TIME is YYYYMMDD.frac-of-day. Subtract date and save just time fraction
|
|
TIME_FRAC(1:N_IASI_NOB) = IASI_TIME(1:N_IASI_NOB)/240000d0
|
|
ENDIF
|
|
|
|
!IF(.NOT. DATA_PRESENT) THEN
|
|
!PRINT *,"No IASI data present for this day, nothing to do here."
|
|
!RETURN
|
|
!ENDIF
|
|
|
|
! Get the range of TES retrievals for the current hour
|
|
CALL GET_NT_RANGE( N_IASI_NOB, GET_NHMS(), TIME_FRAC, NTSTART, NTSTOP )
|
|
|
|
IF ( NTSTART == 0 .and. NTSTOP == 0 ) THEN
|
|
|
|
print*, ' No matching IASI CO obs for this hour'
|
|
RETURN
|
|
ENDIF
|
|
|
|
!print*, ' for hour range: ', GET_NHMS(), IASI_TIME(NTSTART), IASI_TIME(NTSTOP)
|
|
!print*, ' found record range: ', NTSTART, NTSTOP
|
|
|
|
|
|
DO I_IASI = NTSTART, NTSTOP, -1
|
|
IF ( (IASI_QUAL_FLAG(I_IASI) .EQ. 0d0 ) .AND. &
|
|
(IASI_CLOUD_COVER(I_IASI) .EQ. 0d0) .AND. &
|
|
( IASI_LAT(I_IASI) > -60d0 ) .AND. &
|
|
( IASI_LAT(I_IASI) < 75d0 ) .AND. &
|
|
( IASI_CO(I_IASI) > 0d0 ) ) THEN
|
|
! Get model grid coordinate indices that correspond to the observation
|
|
! For safety, initialize these up to LLTES
|
|
GC_CO = 0d0
|
|
IASI_APR_WEB= 0d0
|
|
CO_PERT = 0d0
|
|
CO_HAT = 0d0
|
|
IASI_RATIO = 0d0
|
|
|
|
LIASI = IASI_CO_LVL(I_IASI)
|
|
IIJJ = GET_IJ(REAL(IASI_LON(I_IASI),4),REAL(IASI_LAT(I_IASI),4))
|
|
I = IIJJ(1)
|
|
J = IIJJ(2)
|
|
L0 = N_IASI_NLA - LIASI
|
|
|
|
! Get GC surface pressure (mbar)
|
|
GC_PSURF = GET_PEDGE(I,J,1)
|
|
GC_ALT(I,J,1) = ALT_SURF(I,J)*1d-3
|
|
|
|
DO L = 1, LLPAR
|
|
JLOOP = JLOP(I,J,L)
|
|
GC_PRES(L) = GET_PCENTER(I,J,L)
|
|
GC_ALT(I,J,L+1) = (SUM(BXHEIGHT(I,J,1:L)) + ALT_SURF(I,J))*1d-3
|
|
GC_CO_NATIVE(L) = CHK_STT(I,J,L,IDTCO) * TCVV(IDTCO)* XNUMOLAIR* BXHEIGHT(I,J,L) *100d0 /(AIRVOL(I,J,L)* 1d6)
|
|
ENDDO
|
|
!PRINT *, "GC_ALT", GC_ALT(I,J,:)
|
|
!PRINT *, "IASI_ALT", IASI_ALT(:)
|
|
!PRINT *, "GC_CO_NATIVE", SUM(GC_CO_NATIVE(:))
|
|
CALL BIN_DATA_IASI(GC_ALT(I,J,:),IASI_ALT(L0+1:L0+LIASI),GC_CO_NATIVE(:), GC_CO, IASI_APR(L0+1:L0+LIASI),LIASI, 1)
|
|
!PRINT *, "GC_ALT", GC_ALT(I,J,:)
|
|
!PRINT *, "IASI_ALT", IASI_ALT(L0+1:L0+LIASI)
|
|
!PRINT *, " GC_CO_NATIVE", GC_CO_NATIVE(:)
|
|
!PRINT *, " GC_CO", GC_CO(:)
|
|
!--------------------------------------------------------------
|
|
! Apply IASI CO observation operator
|
|
!
|
|
! x_hat = x_a + A_k ( x_m - x_a )
|
|
!
|
|
! where
|
|
! x_hat = GC modeled column as seen by IASI [lnvmr]
|
|
! x_a = IASI apriori column [lnvmr]
|
|
! x_m = GC modeled column [lnvmr]
|
|
! A_k = IASI averaging kernel
|
|
!
|
|
! OR
|
|
!
|
|
! x_smoothed = A_k*x_m + (I-A_k)*x_a
|
|
! where
|
|
! x_smoothed = GC modeled column smoothed by IASI [vmr]
|
|
! x_a = IASI apriori partial column [vmr]
|
|
! x_m = GC modeled profile [vmr]
|
|
! A_k = IASI averaging kernel
|
|
!--------------------------------------------------------------
|
|
! x_m - x_a
|
|
! x_a + A_k * ( x_m - x_a )
|
|
!PRINT *, "GC_CO", SUM(GC_CO(:))
|
|
!PRINT *, "IASI_CO_AVK", IASI_CO_AVK(I_IASI,:)
|
|
DO L = 1, LIASI
|
|
IF ( (IASI_CO_APR(I_IASI,L+L0) > 0) .AND. &
|
|
(GC_CO(L) > 0) ) THEN
|
|
IASI_APR_WEB(L) = IASI_APR(L+L0) * AIRDEN(L,I,J) * XNUMOLAIR * 1d-14 * 1d5
|
|
!IASI_RATIO(L) = IASI_CO_APR(I_IASI,L+L0)/IASI_APR_WEB(L)
|
|
!CO_PERT(L) = IASI_RATIO(L)*GC_CO(L) - IASI_CO_APR(I_IASI,L+L0)
|
|
CO_PERT(L) = GC_CO(L) - IASI_CO_APR(I_IASI,L+L0)
|
|
CO_HAT(L) = IASI_CO_APR(I_IASI,L+L0) + IASI_CO_AVK(I_IASI,L0+L) * CO_PERT(L)
|
|
ENDIF
|
|
! actual comparison
|
|
ENDDO
|
|
CO_COL_HAT = SUM(CO_HAT(1:LIASI))
|
|
SOBS_CO_STD = (1d0 - IASI_CO_STD(I_IASI)/IASI_CO(I_IASI))**2
|
|
SOBS_COUNT(I,J) = SOBS_COUNT(I,J) + 1d0
|
|
SOBS_CO_NORM(I,J) = SOBS_CO_NORM(I,J) + SOBS_CO_STD
|
|
SOBS_CO_MEAN(I,J) = SOBS_CO_MEAN(I,J) + IASI_CO(I_IASI) * SOBS_CO_STD
|
|
SOBS_CO_TOT(I,J) = SOBS_CO_TOT(I,J) + IASI_CO(I_IASI)
|
|
SOBS_CO_HAT(I,J) = SOBS_CO_HAT(I,J) + CO_COL_HAT
|
|
SOBS_AVK_TOT(I,J,:) = SOBS_AVK_TOT(I,J,:) + IASI_CO_AVK(I_IASI,:) * SOBS_CO_STD
|
|
!PRINT *, "SOBS_CO_STD", SOBS_CO_STD
|
|
ENDIF
|
|
ENDDO
|
|
DO I = 1, IIPAR
|
|
DO J = 1, JJPAR
|
|
IF ( ( SOBS_CO_NORM(I,J) > 0d0 ) .AND. &
|
|
( SOBS_COUNT(I,J) > 0d0 ) ) THEN
|
|
SOBS_CO(I,J) = SOBS_CO_MEAN(I,J)/SOBS_CO_NORM(I,J)
|
|
SOBS_CO_AVG(I,J) = SOBS_CO_TOT(I,J)/SOBS_COUNT(I,J)
|
|
SOBS_CO_AVK(I,J,:) = SOBS_AVK_TOT(I,J,:)/SOBS_CO_NORM(I,J)
|
|
SOBS_HAT(I,J) = SOBS_CO_HAT(I,J)/SOBS_COUNT(I,J)
|
|
SOBS_CO_ERR(I,J) = SOBS_CO(I,J)*(1d0 - SQRT(SOBS_CO_NORM(I,J)/SOBS_COUNT(I,J)))
|
|
ENDIF
|
|
ENDDO
|
|
ENDDO
|
|
DO I_IASI = NTSTART, NTSTOP, -1
|
|
IF ( (IASI_QUAL_FLAG(I_IASI) .EQ. 0d0 ) .AND. &
|
|
(IASI_CLOUD_COVER(I_IASI) .EQ. 0d0) .AND. &
|
|
( IASI_LAT(I_IASI) > -60d0 ) .AND. &
|
|
( IASI_LAT(I_IASI) < 75d0 ) .AND. &
|
|
( IASI_CO(I_IASI) > 0d0 ) ) THEN
|
|
! Get model grid coordinate indices that correspond to the observation
|
|
! For safety, initialize these up to LLTES
|
|
|
|
IIJJ = GET_IJ(REAL(IASI_LON(I_IASI),4),REAL(IASI_LAT(I_IASI),4))
|
|
I = IIJJ(1)
|
|
J = IIJJ(2)
|
|
|
|
IF ( ( SOBS_CO_NORM(I,J) > 0d0 ) .AND. &
|
|
( SOBS_COUNT(I,J) > 0d0 ) ) THEN
|
|
SOBS_CO_VAR(I,J) = SOBS_CO_VAR(I,J) + (IASI_CO(I_IASI) - SOBS_CO_AVG(I,J))**2
|
|
ENDIF
|
|
ENDIF
|
|
ENDDO
|
|
|
|
DO I=1,IIPAR
|
|
DO J=1,JJPAR
|
|
FORCE = 0d0
|
|
DIFF = 0d0
|
|
CO_PERT_ADJ = 0d0
|
|
IF ( SOBS_COUNT(I,J) > 3d0 ) THEN
|
|
SOBS_VAR = SOBS_CO_VAR(I,J)/SOBS_COUNT(I,J)
|
|
SOBS_VAR = MIN(SOBS_VAR,SOBS_VAR_LIMIT)
|
|
SOBS_ERR = SQRT(SOBS_CO_ERR(I,J)**2/SOBS_COUNT(I,J)+SOBS_VAR)
|
|
CALL RANDOM_NUMBER(SOBS_RAND)
|
|
SOBS_CO_COL = SOBS_CO(I,J) !- SQRT(SOBS_VAR)*SOBS_RAND
|
|
DIFF = SOBS_HAT(I,J) - SOBS_CO_COL
|
|
FORCE = DIFF/SOBS_ERR**2
|
|
NEW_COST(I,J) = 0.5 * DIFF * FORCE
|
|
!PRINT *, "FORCE", FORCE
|
|
!PRINT *, "CO_COL_HAT", SOBS_HAT(I,J)
|
|
!PRINT *, "SOBS_CO_COL", SOBS_CO_COL
|
|
!PRINT *, "SOBS_ERR", SOBS_ERR
|
|
WRITE(912,110) (DIFF/1d12)
|
|
WRITE(902,110) (GET_XMID(I))
|
|
WRITE(901,110) (GET_YMID(J))
|
|
WRITE(904,110) (2*NEW_COST(I,J))
|
|
|
|
!DO L = 1, LIASI
|
|
! adjoint of IASI operator
|
|
!CO_PERT_ADJ(L) = SOBS_CO_AVK(I,J,L) * FORCE
|
|
!CO_PERT_ADJ(L) = CO_PERT_ADJ(L)*IASI_RATIO(L)
|
|
!ENDDO
|
|
!CALL BIN_DATA_IASI(GC_ALT(I,J,:),IASI_ALT(:),GC_CO_NATIVE_ADJ(:), CO_PERT_ADJ, IASI_APR(1+L0:L0+LIASI), LIASI, -1)
|
|
!PRINT *, "CO_PERT_ADJ", CO_PERT_ADJ(:)
|
|
!PRINT *, "GC_CO_NATIVE_ADJ(:)", GC_CO_NATIVE_ADJ(:)
|
|
DO L = 1, LLPAR
|
|
! Adjoint of unit conversion
|
|
ADJ_FORCE(I,J,L,IDTCO) = FORCE * TCVV(IDTCO) * BXHEIGHT(I,J,L) * 100d0 * XNUMOLAIR / (1d6 *AIRVOL(I,J,L))
|
|
!PRINT *, "ADJ_FORCE", ADJ_FORCE(I,J,L,IDTCO)
|
|
!PRINT *, "BXHEIGHT", BXHEIGHT(I,J,L)/AIRVOL(I,J,L)
|
|
IF ( ITS_IN_THE_TROP(I,J,L) ) THEN
|
|
STT_ADJ(I,J,L,IDTCO) = STT_ADJ(I,J,L,IDTCO) + ADJ_FORCE(I,J,L,IDTCO)
|
|
ENDIF
|
|
ENDDO
|
|
!PRINT *, "ADJ_FORCE", ADJ_FORCE(I,J,:,IDTCO)
|
|
COST_FUNC = COST_FUNC + NEW_COST(I,J)
|
|
ENDIF
|
|
ENDDO
|
|
ENDDO
|
|
|
|
110 FORMAT(F18.6,1X)
|
|
!PRINT *, "GC_ADJ_TEMP", GC_ADJ_TEMP(:,:,5)
|
|
!PRINT *, "STT_ADJ BEFORE", STT_ADJ(:,:,5,IDTCO)
|
|
!PRINT *, "GC_ADJ_TEMP", GC_ADJ_TEMP(:,:,5)/SOBS_COUNT(:,:)
|
|
!PRINT *, "STT_ADJ AFTER",STT_ADJ(:,:,5,IDTCO)
|
|
|
|
|
|
IF ( FIRST ) FIRST = .FALSE.
|
|
|
|
! Update cost function
|
|
!COST_FUNC = COST_FUNC + SUM(NEW_COST(NTSTOP:NTSTART))
|
|
|
|
print*, ' Updated value of COST_FUNC = ', COST_FUNC
|
|
print*, ' IASI CO contribution = ', COST_FUNC - OLD_COST
|
|
|
|
! Return to calling program
|
|
END SUBROUTINE CALC_IASI_CO_FORCE
|
|
|
|
!----------------------------------------------------------------------------
|
|
SUBROUTINE CLEANUP_IASI
|
|
|
|
IF(ALLOCATED(IASI_LON)) DEALLOCATE(IASI_LON)
|
|
IF(ALLOCATED(IASI_LAT)) DEALLOCATE(IASI_LAT)
|
|
IF(ALLOCATED(IASI_TIME)) DEALLOCATE(IASI_TIME)
|
|
IF(ALLOCATED(IASI_CO_AVK)) DEALLOCATE(IASI_CO_AVK)
|
|
IF(ALLOCATED(IASI_CO_APR)) DEALLOCATE(IASI_CO_APR)
|
|
IF(ALLOCATED(IASI_CO_STD)) DEALLOCATE(IASI_CO_STD)
|
|
IF(ALLOCATED(IASI_CO)) DEALLOCATE(IASI_CO)
|
|
IF(ALLOCATED(IASI_ALT)) DEALLOCATE(IASI_ALT)
|
|
IF(ALLOCATED(IASI_SOLAR_ZENITH)) DEALLOCATE(IASI_SOLAR_ZENITH)
|
|
IF(ALLOCATED(IASI_CLOUD_COVER)) DEALLOCATE(IASI_CLOUD_COVER)
|
|
IF(ALLOCATED(IASI_QUAL_FLAG)) DEALLOCATE(IASI_QUAL_FLAG)
|
|
IF(ALLOCATED(IASI_CO_LVL)) DEALLOCATE(IASI_CO_LVL)
|
|
!IF(ALLOCATED(TEMPDATA)) DEALLOCATE(TEMPDATA)
|
|
!IF(ALLOCATED(TIME_FRAC)) DEALLOCATE(TIME_FRAC)
|
|
END SUBROUTINE CLEANUP_IASI
|
|
!------------------------------------------------------------------------------
|
|
SUBROUTINE GET_NT_RANGE( N_IASI_NOB, HHMMSS, TIME_FRAC, NTSTART, NTSTOP)
|
|
!
|
|
!******************************************************************************
|
|
! Subroutine GET_NT_RANGE retuns the range of retrieval records for the
|
|
! current model hour
|
|
!
|
|
!
|
|
! Arguments as Input:
|
|
! ============================================================================
|
|
! (1 ) NTES (INTEGER) : Number of TES retrievals in this day
|
|
! (2 ) HHMMSS (INTEGER) : Current model time
|
|
! (3 ) TIME_FRAC (REAL) : Vector of times (frac-of-day) for the TES retrievals
|
|
!
|
|
! Arguments as Output:
|
|
! ============================================================================
|
|
! (1 ) NTSTART (INTEGER) : TES record number at which to start
|
|
! (1 ) NTSTOP (INTEGER) : TES record number at which to stop
|
|
!
|
|
! NOTES:
|
|
!
|
|
!******************************************************************************
|
|
!
|
|
! Reference to f90 modules
|
|
USE ERROR_MOD, ONLY : ERROR_STOP
|
|
USE TIME_MOD, ONLY : YMD_EXTRACT
|
|
|
|
! Arguments
|
|
INTEGER, INTENT(IN) :: N_IASI_NOB
|
|
INTEGER, INTENT(IN) :: HHMMSS
|
|
REAL*8, INTENT(IN) :: TIME_FRAC(N_IASI_NOB)
|
|
INTEGER, INTENT(OUT) :: NTSTART
|
|
INTEGER, INTENT(OUT) :: NTSTOP
|
|
|
|
! Local variables
|
|
INTEGER, SAVE :: NTSAVE
|
|
LOGICAL :: FOUND_ALL_RECORDS, FOUND_BAD_RECORDS
|
|
INTEGER :: NTEST
|
|
INTEGER :: HH, MM, SS
|
|
REAL*8 :: GC_HH_FRAC
|
|
REAL*8 :: H1_FRAC
|
|
|
|
!=================================================================
|
|
! GET_NT_RANGE begins here!
|
|
!=================================================================
|
|
! Initialize
|
|
FOUND_ALL_RECORDS = .FALSE.
|
|
FOUND_BAD_RECORDS = .TRUE.
|
|
NTSTART = 0
|
|
NTSTOP = 0
|
|
|
|
! set NTSAVE to NTES every time we start with a new file
|
|
IF ( HHMMSS == 230000 ) THEN
|
|
NTSAVE = N_IASI_NOB
|
|
IF ( NTSAVE > 1 ) THEN
|
|
DO WHILE (FOUND_BAD_RECORDS == .TRUE.)
|
|
IF (TIME_FRAC(NTSAVE) > 0.5d0) THEN
|
|
FOUND_BAD_RECORDS = .FALSE.
|
|
ELSE
|
|
NTSAVE = NTSAVE - 1
|
|
ENDIF
|
|
ENDDO
|
|
ENDIF
|
|
ENDIF
|
|
|
|
DO WHILE (TIME_FRAC(NTSAVE) < 0)
|
|
NTSAVE = NTSAVE -1
|
|
IF (NTSAVE == 0) EXIT
|
|
ENDDO
|
|
!print*, ' GET_NT_RANGE for ', HHMMSS
|
|
!print*, ' NTSAVE ', NTSAVE
|
|
!print*, ' N_IASI_NOB ', N_IASI_NOB
|
|
|
|
CALL YMD_EXTRACT( HHMMSS, HH, MM, SS )
|
|
|
|
|
|
! Convert HH from hour to fraction of day
|
|
GC_HH_FRAC = REAL(HH,8) / 24d0
|
|
! one hour as a fraction of day
|
|
H1_FRAC = 0d0 / 24d0
|
|
|
|
|
|
! All records have been read already
|
|
IF ( NTSAVE == 0 ) THEN
|
|
|
|
!print*, 'All records have been read already '
|
|
RETURN
|
|
! No records reached yet
|
|
ELSEIF ( TIME_FRAC(NTSAVE) + H1_FRAC < GC_HH_FRAC ) THEN
|
|
!PRINT *, "TIME_FRAC", TIME_FRAC(NTSAVE)
|
|
|
|
print*, 'No records reached yet'
|
|
RETURN
|
|
|
|
!
|
|
ELSEIF ( TIME_FRAC(NTSAVE) + H1_FRAC >= GC_HH_FRAC ) THEN
|
|
! Starting record found
|
|
NTSTART = NTSAVE
|
|
|
|
!print*, ' Starting : TIME_FRAC(NTSTART) ', TIME_FRAC(NTSTART), NTSTART
|
|
|
|
! Now search forward to find stopping record
|
|
NTEST = NTSTART
|
|
|
|
DO WHILE ( FOUND_ALL_RECORDS == .FALSE. )
|
|
|
|
! Advance to the next record
|
|
NTEST = NTEST - 1
|
|
|
|
! Stop if we reach the earliest available record
|
|
IF ( NTEST == 0 ) THEN
|
|
|
|
NTSTOP = NTEST + 1
|
|
FOUND_ALL_RECORDS = .TRUE.
|
|
|
|
!print*, ' Records found '
|
|
!print*, ' NTSTART, NTSTOP = ', NTSTART, NTSTOP
|
|
! Reset NTSAVE
|
|
NTSAVE = NTEST
|
|
|
|
! When the combined test date rounded up to the nearest
|
|
! half hour is smaller than the current model date, the
|
|
! stopping record has been passed.
|
|
ELSEIF ( TIME_FRAC(NTEST) + H1_FRAC < GC_HH_FRAC ) THEN
|
|
|
|
!print*, ' Testing : TIME_FRAC ', TIME_FRAC(NTEST), NTEST
|
|
|
|
NTSTOP = NTEST + 1
|
|
FOUND_ALL_RECORDS = .TRUE.
|
|
|
|
!print*, ' Records found '
|
|
!print*, ' NTSTART, NTSTOP = ', NTSTART, NTSTOP
|
|
|
|
! Reset NTSAVE
|
|
NTSAVE = NTEST
|
|
!ELSE
|
|
!print*, ' still looking ', NTEST
|
|
|
|
ENDIF
|
|
|
|
ENDDO
|
|
|
|
ELSE
|
|
|
|
CALL ERROR_STOP('problem', 'GET_NT_RANGE' )
|
|
|
|
ENDIF
|
|
|
|
! Return to calling program
|
|
END SUBROUTINE GET_NT_RANGE
|
|
|
|
!------------------------------------------------------------------------------------
|
|
|
|
SUBROUTINE BIN_DATA_IASI( GC_EDGE, OBS_EDGE, DATA_MODEL, DATA_IASI, OBS_IASI, LIASI, FB )
|
|
|
|
!******************************************************************************
|
|
!Based on the code from Monika. (zhe 1/19/11)
|
|
!FB = 1 for forward
|
|
!FB = -1 for adjoint
|
|
!******************************************************************************
|
|
|
|
INTEGER :: L, LL, FB
|
|
INTEGER :: LIASI, NB, LVL_CRT1, LVL_CRT2
|
|
REAL*8 :: ALT_MODEL(LLPAR), GC_EDGE(LLPAR+1), HI, LOW
|
|
REAL*8 :: DATA_MODEL(LLPAR), BIN_IASI(LLPAR,LIASI), DIFF_BIN
|
|
REAL*8 :: OBS_EDGE(LIASI+1), DATA_IASI(LIASI), OBS_IASI(LIASI)
|
|
BIN_IASI(:,:) = 0d0
|
|
LVL_CRT1 = 0
|
|
LVL_CRT2 = 0
|
|
!=================================================================
|
|
! BIN_DATA_V4 begins here!
|
|
!=================================================================
|
|
DO L = 1, LLPAR
|
|
ALT_MODEL(L) = 0.5d0*(GC_EDGE(L)+GC_EDGE(L+1))
|
|
ENDDO
|
|
IF (FB > 0) THEN
|
|
|
|
!DO L = 1, LIASI
|
|
!DO LL = 1, LLPAR
|
|
!IF ( ALT_MODEL(LL) >= OBS_EDGE(L) ) THEN
|
|
!DATA_IASI(L) = DATA_MODEL(LL)
|
|
!EXIT
|
|
!ENDIF
|
|
!ENDDO
|
|
!ENDDO
|
|
|
|
DO L = 1, LIASI
|
|
DO LL = 1, LLPAR
|
|
LOW = GC_EDGE(LL)
|
|
HI = GC_EDGE(LL+1)
|
|
IF ( GC_EDGE(LL) >= OBS_EDGE(L)) THEN
|
|
IF ( GC_EDGE(LL+1) <= OBS_EDGE(L+1) ) THEN
|
|
BIN_IASI(LL,L) = 1d0
|
|
!NB = NB + 1
|
|
ELSEIF (GC_EDGE(LL) <= OBS_EDGE(L+1)) THEN
|
|
DIFF_BIN = HI - LOW
|
|
BIN_IASI(LL,L) = ( OBS_EDGE(L+1) - LOW)/DIFF_BIN
|
|
BIN_IASI(LL,L+1) = ( HI - OBS_EDGE(L+1))/DIFF_BIN
|
|
ELSEIF (GC_EDGE(LL) > OBS_EDGE(LIASI+1)) THEN
|
|
BIN_IASI(LL,LIASI) = 1d0
|
|
ENDIF
|
|
ELSEIF (GC_EDGE(LL) < OBS_EDGE(1) ) THEN
|
|
BIN_IASI(LL,1) = 1D0
|
|
ENDIF
|
|
ENDDO
|
|
!IF (NB > 0) DATA_IASI(L) = DATA_TEM !/ NB
|
|
ENDDO
|
|
|
|
DO L = 1, LIASI
|
|
DATA_IASI(L) = 0d0
|
|
DO LL = 1, LLPAR
|
|
DATA_IASI(L) = DATA_IASI(L) + BIN_IASI(LL,L) * DATA_MODEL(LL)
|
|
ENDDO
|
|
ENDDO
|
|
DO L = 2, LIASI-1
|
|
IF (DATA_IASI(L) == 0d0) THEN
|
|
IF ( DATA_IASI(L-1) > 0d0) THEN
|
|
LVL_CRT1 = L-1
|
|
!PRINT *, "DATA_IASI1", DATA_IASI(L-1)
|
|
ENDIF
|
|
IF (DATA_IASI(L+1) > 0d0) THEN
|
|
LVL_CRT2 = L+1
|
|
!PRINT *, "DATA_IASI2", DATA_IASI(L+1)
|
|
ENDIF
|
|
IF (REAL(LVL_CRT1)*REAL(LVL_CRT2) > 0D0) THEN
|
|
DO LL = LVL_CRT1, LVL_CRT2
|
|
DATA_IASI(LL) = ((DATA_IASI(LVL_CRT1)+DATA_IASI(LVL_CRT2))/SUM(OBS_IASI(LVL_CRT1:LVL_CRT2)))*OBS_IASI(LL)
|
|
ENDDO
|
|
!PRINT *, "OBS_IASI", SUM(OBS_IASI(LVL_CRT1:LVL_CRT2))
|
|
LVL_CRT1 = 0
|
|
LVL_CRT2 = 0
|
|
ENDIF
|
|
ENDIF
|
|
ENDDO
|
|
|
|
ELSE
|
|
|
|
DATA_MODEL(:) = 0.
|
|
DO L = 1, LLPAR
|
|
DO LL = 1, LIASI
|
|
IF ( ( ALT_MODEL(L) >= OBS_EDGE(LL)) .and. ( ALT_MODEL(L) < OBS_EDGE(LL+1)) ) THEN
|
|
DATA_MODEL(L) = DATA_IASI(LL)
|
|
ENDIF
|
|
ENDDO
|
|
ENDDO
|
|
|
|
ENDIF
|
|
|
|
|
|
! Return to calling program
|
|
END SUBROUTINE BIN_DATA_IASI
|
|
|
|
|
|
|
|
|
|
FUNCTION GET_IJ_2x25( LON, LAT ) RESULT ( IIJJ )
|
|
|
|
!
|
|
!******************************************************************************
|
|
! Subroutine GET_IJ_2x25 returns I and J index from the 2 x 2.5 grid for a
|
|
! LON, LAT coord. (dkh, 11/08/09)
|
|
!
|
|
!
|
|
! Arguments as Input:
|
|
! ============================================================================
|
|
! (1 ) LON (REAL*8) : Longitude [degrees]
|
|
! (2 ) LAT (REAL*8) : Latitude [degrees]
|
|
!
|
|
! Function result
|
|
! ============================================================================
|
|
! (1 ) IIJJ(1) (INTEGER) : Long index [none]
|
|
! (2 ) IIJJ(2) (INTEGER) : Lati index [none]
|
|
!
|
|
! NOIASI:
|
|
!
|
|
!******************************************************************************
|
|
!
|
|
! Reference to f90 modules
|
|
USE ERROR_MOD, ONLY : ERROR_STOP
|
|
|
|
! Arguments
|
|
REAL*4 :: LAT, LON
|
|
|
|
! Return
|
|
INTEGER :: I, J, IIJJ(2)
|
|
|
|
! Local variables
|
|
REAL*8 :: TLON, TLAT, DLON, DLAT
|
|
REAL*8, PARAMETER :: DISIZE = 2.5d0
|
|
REAL*8, PARAMETER :: DJSIZE = 2.0d0
|
|
INTEGER, PARAMETER :: IIMAX = 144
|
|
INTEGER, PARAMETER :: JJMAX = 91
|
|
|
|
|
|
!=================================================================
|
|
! GET_IJ_2x25 begins here!
|
|
!=================================================================
|
|
|
|
TLON = 180d0 + LON + DISIZE
|
|
TLAT = 90d0 + LAT + DJSIZE
|
|
|
|
I = TLON / DISIZE
|
|
J = TLAT / DJSIZE
|
|
|
|
|
|
IF ( TLON / DISIZE - REAL(I) >= 0.5d0 ) THEN
|
|
I = I + 1
|
|
ENDIF
|
|
|
|
IF ( TLAT / DJSIZE - REAL(J) >= 0.5d0 ) THEN
|
|
J = J + 1
|
|
ENDIF
|
|
|
|
|
|
! Longitude wraps around
|
|
!IF ( I == 73 ) I = 1
|
|
IF ( I == ( IIMAX + 1 ) ) I = 1
|
|
|
|
! Check for impossible values
|
|
IF ( I > IIMAX .or. J > JJMAX .or. I < 1 .or. J < 1 ) THEN
|
|
CALL ERROR_STOP('Error finding grid box', 'GET_IJ_2x25')
|
|
ENDIF
|
|
|
|
IIJJ(1) = I
|
|
IIJJ(2) = J
|
|
|
|
! Return to calling program
|
|
END FUNCTION GET_IJ_2x25
|
|
!-------------------------------------------------------------------
|
|
|
|
END MODULE IASI_CO_OBS_MOD
|