1204 lines
36 KiB
Fortran
1204 lines
36 KiB
Fortran
!------------------------------------------------------------------------------
|
|
! Matthew Cooper - Dalhousie University !
|
|
!------------------------------------------------------------------------------
|
|
!BOP
|
|
!
|
|
! !MODULE: regrid_a2a_mod.F90
|
|
!
|
|
! !DESCRIPTION: Module REGRID\_A2A\_MOD uses an algorithm adapted from MAP\_A2A
|
|
! code to regrid from one horizonatal grid to another.
|
|
!\\
|
|
!\\
|
|
! !INTERFACE:
|
|
!
|
|
MODULE REGRID_A2A_MOD
|
|
!
|
|
! !USES:
|
|
!
|
|
IMPLICIT NONE
|
|
PRIVATE
|
|
!
|
|
! !PRIVATE MEMBER FUNCTIONS:
|
|
!
|
|
PRIVATE :: XMAP
|
|
PRIVATE :: YMAP
|
|
PRIVATE :: READ_INPUT_GRID
|
|
!
|
|
! !PUBLIC MEMBER FUNCTIONS:
|
|
!
|
|
PUBLIC :: DO_REGRID_A2A
|
|
PUBLIC :: DO_REGRID_DKH
|
|
PUBLIC :: MAP_A2A
|
|
!
|
|
! !REVISION HISTORY:
|
|
! 13 Mar 2012 - M. Cooper - Initial version
|
|
! 03 Apr 2012 - M. Payer - Now use functions GET_AREA_CM2(I,J,L),
|
|
! GET_YEDGE(I,J,L) and GET_YSIN(I,J,L) from the
|
|
! new grid_mod.F90
|
|
! 22 May 2012 - L. Murray - Implemented several bug fixes
|
|
! 23 Aug 2012 - R. Yantosca - Add capability for starting from hi-res grids
|
|
! (generic 0.5x0.5, generic 0.25x0.25, etc.)
|
|
! 23 Aug 2012 - R. Yantosca - Add subroutine READ_INPUT_GRID, which reads the
|
|
! grid parameters (lon & lat edges) w/ netCDF
|
|
! 27 Aug 2012 - R. Yantosca - Now parallelize key DO loops
|
|
!EOP
|
|
!------------------------------------------------------------------------------
|
|
!BOC
|
|
CONTAINS
|
|
!EOC
|
|
!------------------------------------------------------------------------------
|
|
! Harvard University Atmospheric Chemistry Modeling Group !
|
|
!------------------------------------------------------------------------------
|
|
!BOP
|
|
!
|
|
! !IROUTINE: do_regrid_a2a
|
|
!
|
|
! !DESCRIPTION: Subroutine DO\_REGRID\_A2A regrids 2-D data in the
|
|
! horizontal direction.
|
|
!\\
|
|
!\\
|
|
! !INTERFACE:
|
|
!
|
|
SUBROUTINE DO_REGRID_A2A( FILENAME, IM, JM, INGRID, OUTGRID, IS_MASS, &
|
|
netCDF )
|
|
!
|
|
! !USES:
|
|
!
|
|
USE GRID_MOD, ONLY : GET_XEDGE
|
|
USE GRID_MOD, ONLY : GET_YSIN
|
|
USE GRID_MOD, ONLY : GET_AREA_CM2
|
|
USE GRID_MOD, ONLY : GET_IJ
|
|
USE FILE_MOD, ONLY : IOERROR
|
|
USE inquireMod, ONLY : findFreeLUN
|
|
|
|
# include "CMN_SIZE" ! Size parameters
|
|
# include "CMN_GCTM" ! Size parameters
|
|
!
|
|
! !INPUT PARAMETERS:
|
|
!
|
|
! Name of file with lon and lat edge information on the INPUT GRID
|
|
CHARACTER(LEN=*), INTENT(IN) :: FILENAME
|
|
|
|
! Number of lon centers and lat centers on the INPUT GRID
|
|
INTEGER, INTENT(IN) :: IM
|
|
INTEGER, INTENT(IN) :: JM
|
|
|
|
! Data array on the input grid
|
|
REAL*8, INTENT(IN) :: INGRID(IM,JM)
|
|
|
|
! IS_MASS=0 if data is units of concentration (molec/cm2/s, unitless, etc.)
|
|
! IS_MASS=1 if data is units of mass (kg/yr, etc.)
|
|
INTEGER, INTENT(IN) :: IS_MASS
|
|
|
|
! Read from netCDF file? (needed for debugging, will disappear later)
|
|
LOGICAL, OPTIONAL,INTENT(IN) :: netCDF
|
|
!
|
|
! !OUTPUT PARAMETERS:
|
|
!
|
|
! Data array on the OUTPUT GRID
|
|
REAL*8, INTENT(OUT) :: OUTGRID(IIPAR,JJPAR)
|
|
!
|
|
! !REVISION HISTORY:
|
|
|
|
! 13 Mar 2012 - M. Cooper - Initial version
|
|
! 22 May 2012 - L. Murray - Bug fix: INSIN should be allocated w/ JM+1.
|
|
! 22 May 2012 - R. Yantosca - Updated comments, cosmetic changes
|
|
! 25 May 2012 - R. Yantosca - Bug fix: declare the INGRID argument as
|
|
! INTENT(IN) to preserve the values of INGRID
|
|
! in the calling routine
|
|
! 06 Aug 2012 - R. Yantosca - Now make IU_REGRID a local variable
|
|
! 06 Aug 2012 - R. Yantosca - Move calls to findFreeLUN out of DEVEL block
|
|
! 23 Aug 2012 - R. Yantosca - Now use f10.4 format for hi-res grids
|
|
! 23 Aug 2012 - R. Yantosca - Now can read grid info from netCDF files
|
|
! 27 Aug 2012 - R. Yantosca - Add parallel DO loops
|
|
!EOP
|
|
!------------------------------------------------------------------------------
|
|
!BOC
|
|
!
|
|
! !LOCAL VARIABLES:
|
|
!
|
|
! Scalars
|
|
INTEGER :: I, J
|
|
INTEGER :: IOS, M
|
|
INTEGER :: IU_REGRID
|
|
REAL*8 :: INAREA, RLAT
|
|
CHARACTER(LEN=15) :: HEADER1
|
|
CHARACTER(LEN=20) :: FMT_LAT, FMT_LON, FMT_LEN
|
|
LOGICAL :: USE_NETCDF
|
|
|
|
! Arrays
|
|
REAL*8 :: INLON (IM +1) ! Lon edges on INPUT GRID
|
|
REAL*8 :: INSIN (JM +1) ! SIN( lat edges ) on INPUT GRID
|
|
REAL*8 :: OUTLON (IIPAR+1) ! Lon edges on OUTPUT GRID
|
|
REAL*8 :: OUTSIN (JJPAR+1) ! SIN( lat edges ) on OUTPUT GRID
|
|
REAL*8 :: IN_GRID(IM,JM ) ! Shadow variable for INGRID
|
|
|
|
! dkh debug
|
|
REAL*8 :: total
|
|
REAL*4 :: LON
|
|
REAL*4 :: LAT
|
|
INTEGER :: IIJJ(2)
|
|
|
|
!======================================================================
|
|
! Initialization
|
|
!
|
|
! NOTE: In the near future ASCII input will be replaced by netCDF!
|
|
!======================================================================
|
|
|
|
! Save value of netCDF to shadow variable
|
|
IF ( PRESENT( netCDF ) ) THEN
|
|
USE_netCDF = netCDF
|
|
ELSE
|
|
USE_netCDF = .FALSE.
|
|
ENDIF
|
|
|
|
! Longitude edges on the OUTPUT GRID
|
|
! NOTE: May have to make OUTLON a 2-D array later for the GI model
|
|
DO I = 1, IIPAR+1
|
|
OUTLON(I) = GET_XEDGE( I )
|
|
ENDDO
|
|
|
|
! SIN( lat edges ) on the OUTPUT GRID
|
|
! NOTE: May have to make OUTSIN a 2-D array later for the GI model
|
|
DO J = 1, JJPAR+1
|
|
OUTSIN(J) = GET_YSIN( 1, J, 1 )
|
|
ENDDO
|
|
|
|
! Read the input grid specifications
|
|
IF ( USE_netCDF ) THEN
|
|
|
|
!------------------------------------------
|
|
! %%% FROM NETCDF FILE %%%
|
|
!------------------------------------------
|
|
|
|
! Read the grid specifications from a netCDF file
|
|
CALL READ_INPUT_GRID( IM, JM, FILENAME, INLON, INSIN )
|
|
|
|
ELSE
|
|
|
|
!------------------------------------------
|
|
! %%% FROM ASCII FILE %%%
|
|
!
|
|
! NOTE: Deprecated, will be removed later.
|
|
!------------------------------------------
|
|
|
|
! Find a free file LUN
|
|
IU_REGRID = findFreeLUN()
|
|
|
|
! Open file containing lon & lat edges on the INPUT GRID
|
|
OPEN( IU_REGRID, FILE=TRIM( FILENAME ), STATUS='OLD', IOSTAT=IOS )
|
|
IF ( IOS /= 0 ) CALL IOERROR( IOS, IU_REGRID, 'latlonread' )
|
|
|
|
! Create the approprate FORMAT strings
|
|
WRITE(FMT_LEN,*) IM+1
|
|
|
|
! NOTE: If the resolution of the grid is high enough, we have
|
|
! to allow for an extra digit in the input file. This will
|
|
! become obsolete once we migrate to netCDF format (bmy, 8/23/12)
|
|
IF ( IM > 1000 ) THEN
|
|
FMT_LON='(' // TRIM ( FMT_LEN ) // 'F10.4)' ! For hi-res grids
|
|
ELSE
|
|
FMT_LON='(' // TRIM ( FMT_LEN ) // 'F9.3)' ! For all other grids
|
|
ENDIF
|
|
|
|
WRITE(FMT_LEN,*) JM
|
|
FMT_LAT='(' // TRIM ( FMT_LEN ) // 'F15.10)'
|
|
|
|
! Read lon edges & SIN( lat edges ) on the INPUT GRID
|
|
READ( IU_REGRID, '(A15)',IOSTAT=IOS ) HEADER1
|
|
READ( IU_REGRID,FMT_LON,IOSTAT=IOS ) ( INLON(M), M=1,IM+1 )
|
|
READ( IU_REGRID,FMT_LAT,IOSTAT=IOS ) ( INSIN(M), M=1,JM+1 )
|
|
|
|
! Close file
|
|
CLOSE( IU_REGRID )
|
|
|
|
ENDIF
|
|
|
|
|
|
|
|
!======================================================================
|
|
! Regridding
|
|
!======================================================================
|
|
|
|
! Copy the input argument INGRID to a local shadow variable,
|
|
! so that we can preserve the value of INGRID in the calling routine
|
|
IN_GRID = INGRID
|
|
|
|
! Convert input to per area units if necessary
|
|
IF ( IS_MASS == 1 ) THEN
|
|
|
|
!$OMP PARALLEL DO &
|
|
!$OMP DEFAULT( SHARED ) &
|
|
!$OMP PRIVATE( I, J, RLAT, INAREA )
|
|
DO J = 1, JM
|
|
RLAT = INSIN(J+1) - INSIN(J)
|
|
INAREA = ( 2d0 * PI * Re * RLAT * 1d4 * Re ) / DBLE( IM )
|
|
DO I = 1, IM
|
|
IN_GRID(I,J) = IN_GRID(I,J) / INAREA
|
|
ENDDO
|
|
ENDDO
|
|
!$OMP END PARALLEL DO
|
|
|
|
ENDIF
|
|
|
|
! dkh debug
|
|
! display total ingrid
|
|
IF ( IS_MASS == 3 ) THEN
|
|
total = 0d0
|
|
DO J = 1, JM
|
|
RLAT = INSIN(J+1) - INSIN(J)
|
|
INAREA = ( 2d0 * PI * Re * RLAT * 1d4 * Re ) / DBLE( IM )
|
|
DO I = 1, IM
|
|
total = total + IN_GRID(I,J) * INAREA * 1d-4 ! 1d-4 because INAREA is cm2, but HTAP is per m2
|
|
ENDDO
|
|
ENDDO
|
|
! and convert from kg per s to kg per month (July)
|
|
print*, ' HTAP CO sum in = ', total * 60d0 * 60d0 * 24d0 * 31d0
|
|
ENDIF
|
|
|
|
! DO J = 1, JM
|
|
! RLAT = INSIN(J+1) - INSIN(J)
|
|
! INAREA = ( 2d0 * PI * Re * RLAT * 1d4 * Re ) / DBLE( IM )
|
|
! DO I = 1, IM
|
|
! LAT = -179.95d0 + (JM - 1 ) * 0.1d0
|
|
! LON = - 89.95d0 + (IM - 1 ) * 0.1d0
|
|
! IIJJ = GET_IJ(LON, LAT)
|
|
! OUTGRID(IIJJ(1),IIJJ(2)) = OUTGRID(IIJJ(1),IIJJ(2)) &
|
|
! + IN_GRID(I,J) * INAREA * 1d-4 ! 1d-4 because INAREA is cm2, but HTAP is per m2
|
|
! ENDDO
|
|
! ENDDO
|
|
! DO J = 1, JJPAR
|
|
! DO I = 1, IIPAR
|
|
! OUTGRID(I,J) = OUTGRID(I,J) / GET_AREA_CM2( J ) * 1d4 ! cm2 to m2
|
|
! ENDDO
|
|
! ENDDO
|
|
|
|
! Call MAP_A2A to do the regridding
|
|
CALL MAP_A2A( IM, JM, INLON, INSIN, IN_GRID, &
|
|
IIPAR, JJPAR, OUTLON, OUTSIN, OUTGRID, 0, 0 )
|
|
|
|
! Convert back from "per area" if necessary
|
|
IF ( IS_MASS == 1 ) THEN
|
|
|
|
!$OMP PARALLEL DO &
|
|
!$OMP DEFAULT( SHARED ) &
|
|
!$OMP PRIVATE( I, J )
|
|
DO J = 1, JJPAR
|
|
DO I = 1, IIPAR
|
|
OUTGRID(I,J) = OUTGRID(I,J) * GET_AREA_CM2( J )
|
|
ENDDO
|
|
ENDDO
|
|
!$OMP END PARALLEL DO
|
|
|
|
ENDIF
|
|
|
|
! dkh debug
|
|
! display total ingrid
|
|
IF ( IS_MASS == 3 ) THEN
|
|
total = 0d0
|
|
DO J = 1, JJPAR
|
|
DO I = 1, IIPAR
|
|
total = total + OUTGRID(I,J) * GET_AREA_CM2( J ) * 1d-4 ! cm2 to m2
|
|
ENDDO
|
|
ENDDO
|
|
print*, ' HTAP CO sum in = ', total * 60d0 * 60d0 * 24d0 * 31d0
|
|
ENDIF
|
|
END SUBROUTINE DO_REGRID_A2A
|
|
!EOC
|
|
!------------------------------------------------------------------------------
|
|
!BOP
|
|
!
|
|
! !IROUTINE: do_regrid_dkh
|
|
!
|
|
! !DESCRIPTION: Subroutine DO\_REGRID\_DKH regrids 2-D data in the
|
|
! horizontal direction.
|
|
!\\
|
|
!\\
|
|
! !INTERFACE:
|
|
!
|
|
SUBROUTINE DO_REGRID_DKH( FILENAME, IM, JM, INGRID, OUTGRID, IS_MASS, &
|
|
netCDF )
|
|
!
|
|
! !USES:
|
|
!
|
|
USE GRID_MOD, ONLY : GET_XEDGE
|
|
USE GRID_MOD, ONLY : GET_YSIN
|
|
USE GRID_MOD, ONLY : GET_AREA_CM2
|
|
USE GRID_MOD, ONLY : GET_IJ
|
|
USE GRID_MOD, ONLY : GET_IJ_GLOBAL
|
|
USE FILE_MOD, ONLY : IOERROR
|
|
USE inquireMod, ONLY : findFreeLUN
|
|
USE GRID_MOD, ONLY : GET_XOFFSET, GET_YOFFSET
|
|
|
|
# include "CMN_SIZE" ! Size parameters
|
|
# include "CMN_GCTM" ! Size parameters
|
|
!
|
|
! !INPUT PARAMETERS:
|
|
!
|
|
! Name of file with lon and lat edge information on the INPUT GRID
|
|
CHARACTER(LEN=*), INTENT(IN) :: FILENAME
|
|
|
|
! Number of lon centers and lat centers on the INPUT GRID
|
|
INTEGER, INTENT(IN) :: IM
|
|
INTEGER, INTENT(IN) :: JM
|
|
|
|
! Data array on the input grid
|
|
REAL*8, INTENT(IN) :: INGRID(IM,JM)
|
|
|
|
! IS_MASS=0 if data is units of concentration (molec/cm2/s, unitless, etc.)
|
|
! IS_MASS=1 if data is units of mass (kg/yr, etc.)
|
|
INTEGER, INTENT(IN) :: IS_MASS
|
|
|
|
! Read from netCDF file? (needed for debugging, will disappear later)
|
|
LOGICAL, OPTIONAL,INTENT(IN) :: netCDF
|
|
!
|
|
! !OUTPUT PARAMETERS:
|
|
!
|
|
! Data array on the OUTPUT GRID
|
|
REAL*8, INTENT(OUT) :: OUTGRID(IIPAR,JJPAR)
|
|
!
|
|
! !REVISION HISTORY:
|
|
|
|
! 17 Nov 2013 - D. Henze - Initial version, based on A2A
|
|
!EOP
|
|
!------------------------------------------------------------------------------
|
|
!BOC
|
|
!
|
|
! !LOCAL VARIABLES:
|
|
!
|
|
! Scalars
|
|
INTEGER :: I, J
|
|
INTEGER :: IOS, M
|
|
INTEGER :: IU_REGRID
|
|
REAL*8 :: INAREA, RLAT
|
|
CHARACTER(LEN=15) :: HEADER1
|
|
CHARACTER(LEN=20) :: FMT_LAT, FMT_LON, FMT_LEN
|
|
LOGICAL :: USE_NETCDF
|
|
|
|
! Arrays
|
|
REAL*8 :: INLON (IM +1) ! Lon edges on INPUT GRID
|
|
REAL*8 :: INSIN (JM +1) ! SIN( lat edges ) on INPUT GRID
|
|
REAL*8 :: OUTLON (IIPAR+1) ! Lon edges on OUTPUT GRID
|
|
REAL*8 :: OUTSIN (JJPAR+1) ! SIN( lat edges ) on OUTPUT GRID
|
|
REAL*8 :: IN_GRID(IM,JM ) ! Shadow variable for INGRID
|
|
|
|
! (yd 2015/6/19 Global offset for nested)
|
|
INTEGER :: P, Q
|
|
REAL*8 :: P0, Q0
|
|
REAL*8 :: OUTGRID_GLOB(IIPAR_L,JJPAR_L)
|
|
!(quzhen 2016/2/1)
|
|
INTEGER :: LON_COUNT
|
|
INTEGER :: IIJJ_PREV !lon index of previous longitude
|
|
REAL*8 :: RATIO
|
|
|
|
! dkh debug
|
|
REAL*8 :: total
|
|
REAL*4 :: LON
|
|
REAL*4 :: LAT
|
|
INTEGER :: IIJJ(2)
|
|
|
|
!======================================================================
|
|
! Initialization
|
|
!
|
|
! NOTE: In the near future ASCII input will be replaced by netCDF!
|
|
!======================================================================
|
|
|
|
! Save value of netCDF to shadow variable
|
|
IF ( PRESENT( netCDF ) ) THEN
|
|
USE_netCDF = netCDF
|
|
ELSE
|
|
USE_netCDF = .FALSE.
|
|
ENDIF
|
|
|
|
! Longitude edges on the OUTPUT GRID
|
|
! NOTE: May have to make OUTLON a 2-D array later for the GI model
|
|
DO I = 1, IIPAR+1
|
|
OUTLON(I) = GET_XEDGE( I )
|
|
ENDDO
|
|
|
|
! SIN( lat edges ) on the OUTPUT GRID
|
|
! NOTE: May have to make OUTSIN a 2-D array later for the GI model
|
|
DO J = 1, JJPAR+1
|
|
OUTSIN(J) = GET_YSIN( 1, J, 1 )
|
|
ENDDO
|
|
|
|
! Read the input grid specifications
|
|
IF ( USE_netCDF ) THEN
|
|
|
|
!------------------------------------------
|
|
! %%% FROM NETCDF FILE %%%
|
|
!------------------------------------------
|
|
|
|
! Read the grid specifications from a netCDF file
|
|
CALL READ_INPUT_GRID( IM, JM, FILENAME, INLON, INSIN )
|
|
|
|
ELSE
|
|
|
|
!------------------------------------------
|
|
! %%% FROM ASCII FILE %%%
|
|
!
|
|
! NOTE: Deprecated, will be removed later.
|
|
!------------------------------------------
|
|
|
|
! Find a free file LUN
|
|
IU_REGRID = findFreeLUN()
|
|
|
|
! Open file containing lon & lat edges on the INPUT GRID
|
|
OPEN( IU_REGRID, FILE=TRIM( FILENAME ), STATUS='OLD', IOSTAT=IOS )
|
|
IF ( IOS /= 0 ) CALL IOERROR( IOS, IU_REGRID, 'latlonread' )
|
|
|
|
! Create the approprate FORMAT strings
|
|
WRITE(FMT_LEN,*) IM+1
|
|
|
|
! NOTE: If the resolution of the grid is high enough, we have
|
|
! to allow for an extra digit in the input file. This will
|
|
! become obsolete once we migrate to netCDF format (bmy, 8/23/12)
|
|
IF ( IM > 1000 ) THEN
|
|
FMT_LON='(' // TRIM ( FMT_LEN ) // 'F10.4)' ! For hi-res grids
|
|
ELSE
|
|
FMT_LON='(' // TRIM ( FMT_LEN ) // 'F9.3)' ! For all other grids
|
|
ENDIF
|
|
|
|
WRITE(FMT_LEN,*) JM
|
|
FMT_LAT='(' // TRIM ( FMT_LEN ) // 'F15.10)'
|
|
|
|
! Read lon edges & SIN( lat edges ) on the INPUT GRID
|
|
READ( IU_REGRID, '(A15)',IOSTAT=IOS ) HEADER1
|
|
READ( IU_REGRID,FMT_LON,IOSTAT=IOS ) ( INLON(M), M=1,IM+1 )
|
|
READ( IU_REGRID,FMT_LAT,IOSTAT=IOS ) ( INSIN(M), M=1,JM+1 )
|
|
|
|
! Close file
|
|
CLOSE( IU_REGRID )
|
|
|
|
ENDIF
|
|
|
|
|
|
|
|
!======================================================================
|
|
! Regridding
|
|
!======================================================================
|
|
|
|
! Copy the input argument INGRID to a local shadow variable,
|
|
! so that we can preserve the value of INGRID in the calling routine
|
|
IN_GRID = INGRID
|
|
|
|
! Convert input to per area units if necessary
|
|
IF ( IS_MASS == 1 ) THEN
|
|
|
|
!$OMP PARALLEL DO &
|
|
!$OMP DEFAULT( SHARED ) &
|
|
!$OMP PRIVATE( I, J, RLAT, INAREA )
|
|
DO J = 1, JM
|
|
RLAT = INSIN(J+1) - INSIN(J)
|
|
INAREA = ( 2d0 * PI * Re * RLAT * 1d4 * Re ) / DBLE( IM )
|
|
DO I = 1, IM
|
|
IN_GRID(I,J) = IN_GRID(I,J) / INAREA
|
|
ENDDO
|
|
ENDDO
|
|
!$OMP END PARALLEL DO
|
|
|
|
ENDIF
|
|
|
|
DO J = 1, JM
|
|
RLAT = INSIN(J+1) - INSIN(J)
|
|
INAREA = ( 2d0 * PI * Re * RLAT * 1d4 * Re ) / DBLE( IM )
|
|
DO I = 1, IM
|
|
LON = 0.05d0 + (I - 1 ) * 0.1d0
|
|
IF ( LON > 180d0 ) LON = LON - 360d0
|
|
LAT = - 89.95d0 + (J - 1 ) * 0.1d0
|
|
#if .not. defined(NESTED_CH)
|
|
IIJJ = GET_IJ(LON, LAT)
|
|
#elif defined(NESTED_CH)
|
|
!(quzhen 2014/12/31)
|
|
IIJJ = GET_IJ_GLOBAL(LON, LAT)
|
|
RATIO = 1d0
|
|
|
|
!(quzhen 2016/2/1 ONLY DO THIS IN NESTED_CH DOMAIN)
|
|
IF (LON .GT. 69.65 .AND. LON .LT. 150.35) THEN !in domain
|
|
IF (MOD(IIJJ(1), 3) .EQ. 1) THEN
|
|
! lower boundary
|
|
IF (IIJJ(1) .NE. IIJJ_PREV) THEN ! a new grid
|
|
RATIO = (0.75 - 2./3.)/0.1
|
|
LON_COUNT = 0
|
|
ENDIF
|
|
LON_COUNT = LON_COUNT + 1
|
|
IF(LON_COUNT .EQ. 6) THEN
|
|
! upper boundary
|
|
RATIO = (1./3. - 0.25)/0.1
|
|
ENDIF
|
|
ENDIF
|
|
|
|
IF (MOD(IIJJ(1), 3) .EQ. 2) THEN ! 2nd pattern
|
|
! lower boundary
|
|
IF (IIJJ(1) .NE. IIJJ_PREV) THEN
|
|
RATIO = (0.35 - 1./3.)/0.1
|
|
LON_COUNT = 0
|
|
ENDIF
|
|
LON_COUNT = LON_COUNT + 1
|
|
IF(LON_COUNT .EQ. 7) THEN
|
|
! upper boundary
|
|
RATIO = 0.05/0.1
|
|
ENDIF
|
|
ENDIF
|
|
|
|
IF (MOD(IIJJ(1), 3) .EQ. 0) THEN ! 3rd pattern
|
|
! lower boundary
|
|
IF (IIJJ(1) .NE. IIJJ_PREV) THEN
|
|
!INAREA = INAREA * 0.05/0.1
|
|
RATIO = 0.05/0.1
|
|
LON_COUNT = 0
|
|
ENDIF
|
|
LON_COUNT = LON_COUNT + 1
|
|
!print*,'quzhen 3rd pattern,lon_count',lon_count, ratio
|
|
IF(LON_COUNT .EQ. 7) THEN
|
|
! upper boundary
|
|
!INAREA = INAREA * (2./3. - 0.65)/0.1
|
|
RATIO = (2./3. - 0.65)/0.1
|
|
ENDIF
|
|
ENDIF
|
|
|
|
ENDIF
|
|
IIJJ_PREV = IIJJ(1)
|
|
|
|
!(quzhen 2015/1/23 regrid globally first)
|
|
OUTGRID_GLOB(IIJJ(1),IIJJ(2)) = OUTGRID_GLOB(IIJJ(1),IIJJ(2)) &
|
|
+ IN_GRID(I,J) * INAREA * 1d-4 * RATIO ! 1d-4 because INAREA is cm2, but HTAP is per m2
|
|
#endif
|
|
|
|
#if .not. defined(NESTED_CH)
|
|
OUTGRID(IIJJ(1),IIJJ(2)) = OUTGRID(IIJJ(1),IIJJ(2)) &
|
|
+ IN_GRID(I,J) * INAREA * 1d-4 ! 1d-4 because INAREA is cm2, but HTAP is per m2
|
|
#endif
|
|
|
|
ENDDO
|
|
ENDDO
|
|
|
|
#if (defined (NESTED_NA) || defined (NESTED_CH) || defined (NESTED_SD) )
|
|
!(quzhen 2015/1/23 select data in the window area)
|
|
P0 = GET_XOFFSET( GLOBAL=.TRUE. )
|
|
Q0 = GET_YOFFSET( GLOBAL=.TRUE. )
|
|
|
|
DO P = 1, IIPAR
|
|
DO Q = 1, JJPAR
|
|
OUTGRID( P, Q ) = OUTGRID_GLOB ( P+P0, Q+Q0)
|
|
ENDDO
|
|
ENDDO
|
|
#endif
|
|
|
|
DO J = 1, JJPAR
|
|
DO I = 1, IIPAR
|
|
OUTGRID(I,J) = OUTGRID(I,J) / GET_AREA_CM2( J ) * 1d4 ! cm2 to m2
|
|
ENDDO
|
|
ENDDO
|
|
|
|
! Convert back from "per area" if necessary
|
|
IF ( IS_MASS == 1 ) THEN
|
|
|
|
!$OMP PARALLEL DO &
|
|
!$OMP DEFAULT( SHARED ) &
|
|
!$OMP PRIVATE( I, J )
|
|
DO J = 1, JJPAR
|
|
DO I = 1, IIPAR
|
|
OUTGRID(I,J) = OUTGRID(I,J) * GET_AREA_CM2( J )
|
|
ENDDO
|
|
ENDDO
|
|
!$OMP END PARALLEL DO
|
|
|
|
ENDIF
|
|
|
|
END SUBROUTINE DO_REGRID_DKH
|
|
!EOC
|
|
!------------------------------------------------------------------------------
|
|
! Harvard University Atmospheric Chemistry Modeling Group !
|
|
!------------------------------------------------------------------------------
|
|
!BOP
|
|
!
|
|
! !IROUTINE: map_a2a
|
|
!
|
|
! !DESCRIPTION: Subroutine MAP\_A2A is a horizontal arbitrary grid to arbitrary
|
|
! grid conservative high-order mapping regridding routine by S-J Lin.
|
|
!\\
|
|
!\\
|
|
! !INTERFACE:
|
|
!
|
|
! (1 ) INLON (REAL*8 ) : Longitude edges of input grid
|
|
! (2 ) INSIN (REAL*8 ) : Sine of input grid latitude edges
|
|
! (3 ) INGRID (REAL*8 ) : Data array to be regridded
|
|
|
|
SUBROUTINE map_a2a( im, jm, lon1, sin1, q1, &
|
|
in, jn, lon2, sin2, q2, ig, iv)
|
|
!
|
|
! !INPUT PARAMETERS:
|
|
!
|
|
! Longitude and Latitude dimensions of INPUT grid
|
|
INTEGER, INTENT(IN) :: im, jm
|
|
|
|
! Longitude and Latitude dimensions of OUTPUT grid
|
|
INTEGER, INTENT(IN) :: in, jn
|
|
|
|
! IG=0: pole to pole;
|
|
! IG=1 J=1 is half-dy north of south pole
|
|
INTEGER, INTENT(IN) :: ig
|
|
|
|
! IV=0: Regrid scalar quantity
|
|
! IV=1: Regrid vector quantity
|
|
INTEGER, INTENT(IN) :: iv
|
|
|
|
! Longitude edges (degrees) of INPUT and OUTPUT grids
|
|
REAL*8, INTENT(IN) :: lon1(im+1), lon2(in+1)
|
|
|
|
! Sine of Latitude Edges (radians) of INPUT and OUTPUT grids
|
|
REAL*8, INTENT(IN) :: sin1(jm+1), sin2(jn+1)
|
|
|
|
! Quantity on INPUT grid
|
|
REAL*8, INTENT(IN) :: q1(im,jm)
|
|
!
|
|
! !OUTPUT PARAMETERS:
|
|
!
|
|
! Regridded quantity on OUTPUT grid
|
|
REAL*8, INTENT(OUT) :: q2(in,jn)
|
|
!
|
|
! !REVISION HISTORY:
|
|
! (1) Original subroutine by S-J Lin. Converted to F90 freeform format
|
|
! and inserted into "Geos3RegridModule" by Bob Yantosca (9/21/00)
|
|
! (2) Added F90 type declarations to be consistent w/ TypeModule.f90.
|
|
! Also updated comments. (bmy, 9/21/00)
|
|
! 21 Sep 2000 - R. Yantosca - Initial version
|
|
! 27 Aug 2012 - R. Yantosca - Add parallel DO loops
|
|
!EOP
|
|
!------------------------------------------------------------------------------
|
|
!BOC
|
|
!
|
|
! !LOCAL VARIABLES:
|
|
!
|
|
INTEGER :: i,j,k
|
|
REAL*8 :: qtmp(in,jm)
|
|
|
|
!===================================================================
|
|
! E-W regridding
|
|
!===================================================================
|
|
IF ( im .eq. in ) THEN
|
|
|
|
! Don't call XMAP if both grids have the same # of longitudes
|
|
! but save the input data in the QTMP array
|
|
!$OMP PARALLEL DO &
|
|
!$OMP DEFAULT( SHARED ) &
|
|
!$OMP PRIVATE( I, J )
|
|
DO j=1,jm-ig
|
|
DO i=1,im
|
|
qtmp(i,j+ig) = q1(i,j+ig)
|
|
ENDDO
|
|
ENDDO
|
|
!$OMP END PARALLEL DO
|
|
|
|
ELSE
|
|
|
|
! Otherwise, call XMAP to regrid in the E-W direction
|
|
CALL xmap(im, jm-ig, lon1, q1(1,1+ig),in, lon2, qtmp(1,1+ig) )
|
|
|
|
ENDIF
|
|
|
|
!===================================================================
|
|
! N-S regridding
|
|
!===================================================================
|
|
IF ( jm .eq. jn ) THEN
|
|
|
|
! Don't call XMAP if both grids have the same # of longitudes,
|
|
! but assign the value of QTMP to the output Q2 array
|
|
!$OMP PARALLEL DO &
|
|
!$OMP DEFAULT( SHARED ) &
|
|
!$OMP PRIVATE( I, J )
|
|
DO j=1,jm-ig
|
|
DO i=1,in
|
|
q2(i,j+ig) = qtmp(i,j+ig)
|
|
ENDDO
|
|
ENDDO
|
|
!$OMP END PARALLEL DO
|
|
|
|
ELSE
|
|
|
|
! Otherwise, call YMAP to regrid in the N-S direction
|
|
CALL ymap(in, jm, sin1, qtmp(1,1+ig), jn, sin2, q2(1,1+ig), ig, iv)
|
|
|
|
ENDIF
|
|
|
|
END SUBROUTINE map_a2a
|
|
!EOC
|
|
!------------------------------------------------------------------------------
|
|
! Prasad Kasibhatla - Duke University !
|
|
!------------------------------------------------------------------------------
|
|
!BOP
|
|
!
|
|
! !IROUTINE: ymap
|
|
!
|
|
! !DESCRIPTION: Routine to perform area preserving mapping in N-S from an
|
|
! arbitrary resolution to another.
|
|
!\\
|
|
!\\
|
|
! !INTERFACE:
|
|
!
|
|
SUBROUTINE ymap(im, jm, sin1, q1, jn, sin2, q2, ig, iv)
|
|
!
|
|
! !INPUT PARAMETERS:
|
|
!
|
|
|
|
! original E-W dimension
|
|
INTEGER, INTENT(IN) :: im
|
|
|
|
! original N-S dimension
|
|
INTEGER, INTENT(IN) :: jm
|
|
|
|
! Target N-S dimension
|
|
INTEGER, INTENT(IN) :: jn
|
|
|
|
! IG=0: scalars from SP to NP (D-grid v-wind is also IG=0)
|
|
! IG=1: D-grid u-wind
|
|
INTEGER, INTENT(IN) :: ig
|
|
|
|
! IV=0: scalar;
|
|
! IV=1: vector
|
|
INTEGER, INTENT(IN) :: iv
|
|
|
|
! Original southern edge of the cell sin(lat1)
|
|
REAL*8, INTENT(IN) :: sin1(jm+1-ig)
|
|
|
|
! Original data at center of the cell
|
|
REAL*8, INTENT(IN) :: q1(im,jm)
|
|
|
|
! Target cell's southern edge sin(lat2)
|
|
REAL*8, INTENT(IN) :: sin2(jn+1-ig)
|
|
!
|
|
! !OUTPUT PARAMETERS:
|
|
!
|
|
! Mapped data at the target resolution
|
|
REAL*8, INTENT(OUT) :: q2(im,jn)
|
|
!
|
|
! !REMARKS:
|
|
!
|
|
! sin1 (1) = -1 must be south pole; sin1(jm+1)=1 must be N pole.
|
|
!
|
|
! sin1(1) < sin1(2) < sin1(3) < ... < sin1(jm) < sin1(jm+1)
|
|
! sin2(1) < sin2(2) < sin2(3) < ... < sin2(jn) < sin2(jn+1)!
|
|
!
|
|
! !AUTHOR:
|
|
! Developer: Prasad Kasibhatla
|
|
! March 6, 2012
|
|
!
|
|
! !REVISION HISTORY
|
|
! 06 Mar 2012 - P. Kasibhatla - Initial version
|
|
! 27 Aug 2012 - R. Yantosca - Added parallel DO loops
|
|
! 27 Aug 2012 - R. Yantosca - Change REAL*4 variables to REAL*8 to better
|
|
! ensure numerical stability
|
|
!EOP
|
|
!------------------------------------------------------------------------------
|
|
!BOC
|
|
!
|
|
! !LOCAL VARIABLES:
|
|
!
|
|
INTEGER :: i, j0, m, mm, j
|
|
REAL*8 :: dy1(jm)
|
|
REAL*8 :: dy
|
|
!------------------------------------------------------------------------------
|
|
! Prior to 8/27/12:
|
|
! Change REAL*4 to REAL*8, to eliminate numerical noise (bmy, 8/27/12)
|
|
! REAL*4 :: qsum, sum
|
|
!------------------------------------------------------------------------------
|
|
REAL*8 :: qsum, sum
|
|
|
|
! YMAP begins here!
|
|
do j=1,jm-ig
|
|
dy1(j) = sin1(j+1) - sin1(j)
|
|
enddo
|
|
|
|
!===============================================================
|
|
! Area preserving mapping
|
|
!===============================================================
|
|
|
|
!$OMP PARALLEL DO &
|
|
!$OMP DEFAULT( SHARED ) &
|
|
!$OMP PRIVATE( I, J0, J, M, QSUM, MM, DY )
|
|
do 1000 i=1,im
|
|
j0 = 1
|
|
do 555 j=1,jn-ig
|
|
do 100 m=j0,jm-ig
|
|
|
|
!=========================================================
|
|
! locate the southern edge: sin2(i)
|
|
!=========================================================
|
|
if(sin2(j) .ge. sin1(m) .and. sin2(j) .le. sin1(m+1)) then
|
|
|
|
if(sin2(j+1) .le. sin1(m+1)) then
|
|
|
|
! entire new cell is within the original cell
|
|
q2(i,j)=q1(i,m)
|
|
j0 = m
|
|
goto 555
|
|
else
|
|
|
|
! South most fractional area
|
|
qsum=(sin1(m+1)-sin2(j))*q1(i,m)
|
|
|
|
do mm=m+1,jm-ig
|
|
|
|
! locate the northern edge: sin2(j+1)
|
|
if(sin2(j+1) .gt. sin1(mm+1) ) then
|
|
|
|
! Whole layer
|
|
qsum = qsum + dy1(mm)*q1(i,mm)
|
|
else
|
|
|
|
! North most fractional area
|
|
dy = sin2(j+1)-sin1(mm)
|
|
qsum=qsum+dy*q1(i,mm)
|
|
j0 = mm
|
|
goto 123
|
|
endif
|
|
enddo
|
|
goto 123
|
|
endif
|
|
endif
|
|
100 continue
|
|
123 q2(i,j) = qsum / ( sin2(j+1) - sin2(j) )
|
|
555 continue
|
|
1000 continue
|
|
!$OMP END PARALLEL DO
|
|
|
|
!===================================================================
|
|
! Final processing for poles
|
|
!===================================================================
|
|
if ( ig .eq. 0 .and. iv .eq. 0 ) then
|
|
|
|
!------------------------------------------------------------------------------
|
|
! Prior to 8/27/12:
|
|
! Change REAL*4 to REAL*8, to eliminate numerical noise (bmy, 8/27/12)
|
|
! ! South pole
|
|
! sum = 0.
|
|
! do i=1,im
|
|
! sum = sum + q2(i,1)
|
|
! enddo
|
|
!
|
|
! sum = sum / float(im)
|
|
! do i=1,im
|
|
! q2(i,1) = sum
|
|
! enddo
|
|
!
|
|
! ! North pole:
|
|
! sum = 0.
|
|
! do i=1,im
|
|
! sum = sum + q2(i,jn)
|
|
! enddo
|
|
!
|
|
! sum = sum / float(im)
|
|
! do i=1,im
|
|
! q2(i,jn) = sum
|
|
! enddo
|
|
!------------------------------------------------------------------------------
|
|
! South pole
|
|
sum = 0.d0
|
|
do i=1,im
|
|
sum = sum + q2(i,1)
|
|
enddo
|
|
|
|
sum = sum / DBLE( im )
|
|
do i=1,im
|
|
q2(i,1) = sum
|
|
enddo
|
|
|
|
! North pole:
|
|
sum = 0.d0
|
|
do i=1,im
|
|
sum = sum + q2(i,jn)
|
|
enddo
|
|
|
|
sum = sum / DBLE( im )
|
|
do i=1,im
|
|
q2(i,jn) = sum
|
|
enddo
|
|
|
|
endif
|
|
|
|
END SUBROUTINE YMAP
|
|
!EOC
|
|
!------------------------------------------------------------------------------
|
|
! Prasad Kasibhatla - Duke University !
|
|
!------------------------------------------------------------------------------
|
|
!BOP
|
|
!
|
|
! !IROUTINE: xmap
|
|
!
|
|
! !DESCRIPTION: Routine to perform area preserving mapping in E-W from an
|
|
! arbitrary resolution to another.
|
|
! Periodic domain will be assumed, i.e., the eastern wall bounding cell
|
|
! im is lon1(im+1) = lon1(1); Note the equal sign is true geographysically.
|
|
!\\
|
|
!\\
|
|
! !INTERFACE:
|
|
!
|
|
SUBROUTINE xmap(im, jm, lon1, q1, in, lon2, q2)
|
|
!
|
|
! !INPUT PARAMETERS:
|
|
!
|
|
! Original E-W dimension
|
|
INTEGER, INTENT(IN) :: im
|
|
|
|
! Target E-W dimension
|
|
INTEGER, INTENT(IN) :: in
|
|
|
|
! Original N-S dimension
|
|
INTEGER, INTENT(IN) :: jm
|
|
|
|
! Original western edge of the cell
|
|
REAL*8, INTENT(IN) :: lon1(im+1)
|
|
|
|
! Original data at center of the cell
|
|
REAL*8, INTENT(IN) :: q1(im,jm)
|
|
|
|
! Target cell's western edge
|
|
REAL*8, INTENT(IN) :: lon2(in+1)
|
|
!
|
|
! !OUTPUT PARAMETERS:
|
|
!
|
|
! Mapped data at the target resolution
|
|
REAL*8, INTENT(OUT) :: q2(in,jm)
|
|
!
|
|
! !REMARKS:
|
|
! lon1(1) < lon1(2) < lon1(3) < ... < lon1(im) < lon1(im+1)
|
|
! lon2(1) < lon2(2) < lon2(3) < ... < lon2(in) < lon2(in+1)
|
|
!
|
|
! !AUTHOR:
|
|
! Developer: Prasad Kasibhatla
|
|
! March 6, 2012
|
|
!
|
|
! !REVISION HISTORY
|
|
! 06 Mar 2012 - P. Kasibhatla - Initial version
|
|
! 27 Aug 2012 - R. Yantosca - Added parallel DO loops
|
|
! 27 Aug 2012 - R. Yantosca - Change REAL*4 variables to REAL*8 to better
|
|
! ensure numerical stability
|
|
!EOP
|
|
!------------------------------------------------------------------------------
|
|
!BOC
|
|
!
|
|
! !LOCAL VARIABLES:
|
|
!
|
|
INTEGER :: i1, i2, i, i0, m, mm, j
|
|
REAL*8 :: qtmp(-im:im+im)
|
|
REAL*8 :: x1(-im:im+im+1)
|
|
REAL*8 :: dx1(-im:im+im)
|
|
REAL*8 :: dx
|
|
!------------------------------------------------------------------------------
|
|
! Prior to 8/27/12:
|
|
! Change REAL*4 to REAL*8, to eliminate numerical noise (bmy, 8/27/12)
|
|
! REAL*4 :: qsum
|
|
!------------------------------------------------------------------------------
|
|
REAL*8 :: qsum
|
|
LOGICAL :: found
|
|
|
|
! XMAP begins here!
|
|
do i=1,im+1
|
|
x1(i) = lon1(i)
|
|
enddo
|
|
|
|
do i=1,im
|
|
dx1(i) = x1(i+1) - x1(i)
|
|
enddo
|
|
|
|
!===================================================================
|
|
! check to see if ghosting is necessary
|
|
! Western edge:
|
|
!===================================================================
|
|
found = .false.
|
|
i1 = 1
|
|
do while ( .not. found )
|
|
if( lon2(1) .ge. x1(i1) ) then
|
|
found = .true.
|
|
else
|
|
i1 = i1 - 1
|
|
if (i1 .lt. -im) then
|
|
write(6,*) 'failed in xmap western edge '
|
|
stop
|
|
else
|
|
x1(i1) = x1(i1+1) - dx1(im+i1)
|
|
dx1(i1) = dx1(im+i1)
|
|
endif
|
|
endif
|
|
enddo
|
|
|
|
!===================================================================
|
|
! Eastern edge:
|
|
!===================================================================
|
|
found = .false.
|
|
i2 = im+1
|
|
do while ( .not. found )
|
|
if( lon2(in+1) .le. x1(i2) ) then
|
|
found = .true.
|
|
else
|
|
i2 = i2 + 1
|
|
if (i2 .gt. 2*im) then
|
|
write(6,*) 'failed in xmap eastern edge'
|
|
stop
|
|
else
|
|
dx1(i2-1) = dx1(i2-1-im)
|
|
x1(i2) = x1(i2-1) + dx1(i2-1)
|
|
endif
|
|
endif
|
|
enddo
|
|
|
|
!$OMP PARALLEL DO &
|
|
!$OMP DEFAULT( SHARED ) &
|
|
!$OMP PRIVATE( J, QTMP, I, I0, M, QSUM, MM, DX )
|
|
do 1000 j=1,jm
|
|
|
|
!=================================================================
|
|
! Area preserving mapping
|
|
!================================================================
|
|
|
|
qtmp(0)=q1(im,j)
|
|
do i=1,im
|
|
qtmp(i)=q1(i,j)
|
|
enddo
|
|
qtmp(im+1)=q1(1,j)
|
|
|
|
! check to see if ghosting is necessary
|
|
! Western edge
|
|
if ( i1 .le. 0 ) then
|
|
do i=i1,0
|
|
qtmp(i) = qtmp(im+i)
|
|
enddo
|
|
endif
|
|
|
|
! Eastern edge:
|
|
if ( i2 .gt. im+1 ) then
|
|
do i=im+1,i2-1
|
|
qtmp(i) = qtmp(i-im)
|
|
enddo
|
|
endif
|
|
|
|
i0 = i1
|
|
|
|
do 555 i=1,in
|
|
do 100 m=i0,i2-1
|
|
|
|
!=============================================================
|
|
! locate the western edge: lon2(i)
|
|
!=============================================================
|
|
if(lon2(i) .ge. x1(m) .and. lon2(i) .le. x1(m+1)) then
|
|
|
|
if(lon2(i+1) .le. x1(m+1)) then
|
|
|
|
! entire new grid is within the original grid
|
|
q2(i,j)=qtmp(m)
|
|
i0 = m
|
|
goto 555
|
|
else
|
|
|
|
! Left most fractional area
|
|
qsum=(x1(m+1)-lon2(i))*qtmp(m)
|
|
do mm=m+1,i2-1
|
|
|
|
! locate the eastern edge: lon2(i+1)
|
|
if(lon2(i+1) .gt. x1(mm+1) ) then
|
|
|
|
! Whole layer
|
|
qsum = qsum + dx1(mm)*qtmp(mm)
|
|
|
|
else
|
|
! Right most fractional area
|
|
dx = lon2(i+1)-x1(mm)
|
|
qsum=qsum+dx*qtmp(mm)
|
|
i0 = mm
|
|
goto 123
|
|
endif
|
|
enddo
|
|
goto 123
|
|
endif
|
|
endif
|
|
100 continue
|
|
123 q2(i,j) = qsum / ( lon2(i+1) - lon2(i) )
|
|
555 continue
|
|
1000 continue
|
|
!$OMP END PARALLEL DO
|
|
|
|
END SUBROUTINE xmap
|
|
!EOC
|
|
!------------------------------------------------------------------------------
|
|
! Harvard University Atmospheric Chemistry Modeling Group !
|
|
!------------------------------------------------------------------------------
|
|
!BOP
|
|
!
|
|
! !IROUTINE: read_input_grid
|
|
!
|
|
! !DESCRIPTION: Routine to read variables and attributes from a netCDF
|
|
! file. This routine was automatically generated by the Perl script
|
|
! NcdfUtilities/perl/ncCodeRead.
|
|
!\\
|
|
!\\
|
|
! !INTERFACE:
|
|
!
|
|
SUBROUTINE READ_INPUT_GRID( IM, JM, fileName, lon_edges, lat_sines )
|
|
!
|
|
! !USES:
|
|
!
|
|
! Modules for netCDF read
|
|
USE m_netcdf_io_open
|
|
USE m_netcdf_io_get_dimlen
|
|
USE m_netcdf_io_read
|
|
USE m_netcdf_io_readattr
|
|
USE m_netcdf_io_close
|
|
|
|
IMPLICIT NONE
|
|
|
|
# include "netcdf.inc"
|
|
!
|
|
! !INPUT PARAMETERS:
|
|
!
|
|
INTEGER, INTENT(IN) :: IM ! # of longitudes
|
|
INTEGER, INTENT(IN) :: JM ! # of latitudes
|
|
CHARACTER(LEN=*), INTENT(IN) :: fileName ! File w/ grid info
|
|
!
|
|
! !OUTPUT PARAMETERS:
|
|
!
|
|
REAL*8, INTENT(OUT) :: lon_edges(IM+1) ! Lon edges [degrees]
|
|
REAL*8, INTENT(OUT) :: lat_sines(JM+1) ! SIN( latitude edges )
|
|
!
|
|
! !REMARKS:
|
|
! Created with the ncCodeRead script of the NcdfUtilities package,
|
|
! with subsequent hand-editing.
|
|
!
|
|
! !REVISION HISTORY:
|
|
! 23 Aug 2012 - R. Yantosca - Initial version
|
|
!EOP
|
|
!------------------------------------------------------------------------------
|
|
!BOC
|
|
!
|
|
! !LOCAL VARIABLES:
|
|
!
|
|
! Scalars
|
|
INTEGER :: fId ! netCDF file ID
|
|
|
|
! Arrays
|
|
INTEGER :: st1d(1), ct1d(1) ! netCDF start & count
|
|
|
|
!======================================================================
|
|
! Read data from file
|
|
!======================================================================
|
|
|
|
! Open file for reading
|
|
CALL Ncop_Rd( fId, TRIM( fileName ) )
|
|
|
|
! Read lon_edges from file
|
|
st1d = (/ 1 /)
|
|
ct1d = (/ IM+1 /)
|
|
CALL NcRd( lon_edges, fId, "lon_edges", st1d, ct1d )
|
|
|
|
! Read lat_sines from file
|
|
st1d = (/ 1 /)
|
|
ct1d = (/ JM+1 /)
|
|
CALL NcRd( lat_sines, fId, "lat_sines", st1d, ct1d )
|
|
|
|
! Close netCDF file
|
|
CALL NcCl( fId )
|
|
|
|
END SUBROUTINE READ_INPUT_GRID
|
|
!EOC
|
|
END MODULE REGRID_A2A_MOD
|