Files
GEOS-Chem-adjoint-v35-note/code/paranox_mod.f
2018-08-28 00:46:26 -04:00

1010 lines
37 KiB
Fortran

!------------------------------------------------------------------------------
! Harvard University Atmospheric Chemistry Modeling Group !
!------------------------------------------------------------------------------
!BOP
!
! !MODULE: paranox_mod
!
! !DESCRIPTION: Module PARANOX\_MOD contains subroutines for reading and
! interpolating look up tables necessary for the PARANOX (PARAmeterization
! of emitted NOX) ship plume model developed by G.C.M. Vinken.
!\\
!\\
! !INTERFACE:
!
MODULE PARANOX_MOD
!
! !USES:
!
USE inquireMod, ONLY : findFreeLUN
IMPLICIT NONE
PRIVATE
!
! !PUBLIC MEMBER FUNCTIONS:
!
PUBLIC :: READ_PARANOX_LUT
!##############################################
! Prior to 5/31/13:
! Comment out, this is not used (bmy, 5/31/13)
! PUBLIC :: INTERPOLATE_LUT
!##############################################
PUBLIC :: INTERPOLATE_LUT2
PUBLIC :: FRACNOX, INTOPE
!
! !MODULE VARIABLES
!
! fracnox = look up table for fraction of NOx remaining
! for ship emissions (gvinken, 6/6/10)
! intope = look up table for integrated Ozone Production
! Efficiency for ship emiss (gvinken, 6/6/10)
REAL*4 :: fracnox(4,4,4,12,12,4,5)
REAL*4 :: intope(4,4,4,12,12,4,5)
!
! !REMARKS
! References:
! ============================================================================
! (1 ) Vinken, G.C.M., Boersma, K.F., Jacob, D.J., and Meijer, E.W.:
! Accounting for non-linear chemistry of ship plumes in the GEOS-Chem
! global chemistry transport model, Atmos. Chem. Phys., 11, 11707-11722,
! doi:10.5194/acp-11-11707-2011, 2011.
!
! !REVISION HISTORY:
! 06 Feb 2012 - M. Payer - Initial version
! 01 Mar 2012 - R. Yantosca - Use updated GET_LOCALTIME from time_mod.F
! 03 Aug 2012 - R. Yantosca - Move calls to findFreeLUN out of DEVEL block
!EOP
!------------------------------------------------------------------------------
!BOC
CONTAINS
!EOC
!------------------------------------------------------------------------------
! Harvard University Atmospheric Chemistry Modeling Group !
!------------------------------------------------------------------------------
!BOP
!
! !IROUTINE: read_paranox_lut
!
! !DESCRIPTION: Subroutine READ\_PARANOX\_LUT reads look up tables for use in
! the PARANOX ship plume model (G.C.M. Vinken)
!\\
!\\
! !INTERFACE:
!
SUBROUTINE READ_PARANOX_LUT
!
! !USES:
!
USE DIRECTORY_MOD, ONLY : DATA_DIR_1x1
USE FILE_MOD, ONLY : IOERROR
!
! !REVISION HISTORY:
! 06 Feb 2012 - M. Payer - Initial version modified from code provided by
! G.C.M. Vinken
! 01 Aug 2012 - R. Yantosca - Add reference to findFreeLUN from inqure_mod.F90
! 03 Aug 2012 - R. Yantosca - Move calls to findFreeLUN out of DEVEL block
!EOP
!------------------------------------------------------------------------------
!BOC
!
! !LOCAL VARIABLES:
INTEGER :: IOS
CHARACTER(LEN=255) :: FILENAME
INTEGER :: IU_FILE
!=================================================================
! READ_PARANOX_LUT begins here
!=================================================================
!=================================================================
! Read look up table for fraction of NOx remaining for ship
! emissions [unitless]
!=================================================================
! File name
FILENAME = TRIM( DATA_DIR_1x1 ) // 'PARANOX_201202/' //
& 'FracNOx_binary_5hrs_20gs.dat'
! Echo info
WRITE( 6, 100 ) TRIM( FILENAME )
100 FORMAT( 'READ_PARANOX_LUT: Reading ', a )
! Find a free file LUN
IU_FILE = findFreeLUN()
! Open file to read
OPEN( IU_FILE, FILE=TRIM( FILENAME ), FORM="binary",
& IOSTAT=IOS )
IF ( IOS /= 0 ) CALL IOERROR( IOS, IU_FILE, 'read_paranox_lut:1' )
! Read file
READ( IU_FILE, IOSTAT=IOS ) FRACNOX
IF ( IOS /= 0 ) CALL IOERROR( IOS, IU_FILE, 'read_paranox_lut:2' )
!PRINT*, "binary_fracnox: ", FRACNOX(1:4,1,1,2,3,4,4)
! Close file
CLOSE( IU_FILE )
!=================================================================
! Read look up table for integrated Ozone Production Efficiency
! for ship emissions [molec O3 produced / molec NOx lost]
!=================================================================
! File name
FILENAME = TRIM( DATA_DIR_1x1 ) // 'PARANOX_201202/' //
& 'IntOPE_binary_5hrs_20gs.dat'
! Echo info
WRITE( 6, 101 ) TRIM( FILENAME )
101 FORMAT( 'READ_PARANOX_LUT: Reading ', a )
! Find a free file LUN
IU_FILE = findFreeLUN()
! Open file to read
OPEN( IU_FILE, FILE=TRIM( FILENAME ), FORM="BINARY" )
IF ( IOS /= 0 ) CALL IOERROR( IOS, IU_FILE, 'read_paranox_lut:3' )
! Read file
READ( IU_FILE, IOSTAT=IOS ) INTOPE
IF ( IOS /= 0 ) CALL IOERROR( IOS, IU_FILE, 'read_paranox_lut:4' )
!PRINT*, "binary_intope: ", INTOPE(1:4,1,1,2,3,4,4)
! Close file
CLOSE( IU_FILE )
END SUBROUTINE READ_PARANOX_LUT
!EOC
!###############################################################################
! Prior to 5/31/13:
! Comment out, this subroutine is not used here (bmy, 5/31/13)
!!------------------------------------------------------------------------------
!! Harvard University Atmospheric Chemistry Modeling Group !
!!------------------------------------------------------------------------------
!!BOP
!!
!! !IROUTINE: interpolate_lut
!!
!! !DESCRIPTION: Subroutine INTERPOLATE\_LUT returns FracNOx or IntOPE from the
!! lookup table (G.C.M. Vinken, KNMI, June 2010)
!!\\
!!\\
!! !INTERFACE:
!!
! SUBROUTINE INTERPOLATE_LUT( I, J, JO1D, JNO2, fraction_nox, int_ope,
! & Input_Opt, State_Met, State_Chm )
!!
!! !USES:
!!
! USE TRACERID_MOD, ONLY : IDO3, IDTO3, IDTCO
! USE GIGC_Input_Opt_Mod, ONLY : OptInput
! USE GIGC_State_Chm_Mod, ONLY : ChmState
! USE TIME_MOD, ONLY : GET_LOCALTIME
! USE ERROR_MOD, ONLY : ERROR_STOP
! USE GIGC_State_Met_Mod, ONLY : MetState
!
! USE CMN_FJ_MOD ! Photolysis parameters
! USE CMN_SIZE_MOD
!!
!! !INPUT PARAMETERS:
!!
! INTEGER, INTENT(IN) :: I, J
! REAL*8, INTENT(IN) :: JO1D, JNO2
! TYPE(OptInput), INTENT(IN) :: OptInput ! Input options object
! TYPE(MetState), INTENT(IN) :: State_Met ! Meteorology State object
!!
!! !INPUT/OUTPUT PARAMETERS:
!!
! TYPE(ChmState), INTENT(INOUT) :: State_Chm ! Chemistry State object
!!
!! !OUTPUT PARAMETERS:
!!
! REAL, INTENT(OUT) :: fraction_nox, int_ope
!!
!! !REVISION HISTORY:
!! Jun 2010 - G.C.M. Vinken - Initial version
!! 06 Feb 2012 - M. Payer - Moved from emissions_mod.F to paranox_mod.F;
!! Added ProTeX headers
!! 15 Feb 2012 - M. Payer - Add error trap to ensure 0 < fracnox < 1.
!! 09 Nov 2012 - M. Payer - Replaced all met field arrays with State_Met
!! derived type object
!! 28 Nov 2012 - R. Yantosca - Replace SUNCOS_MID w/ State_Met%SUNCOSmid
!! 28 Nov 2012 - R. Yantosca - Replace SUNCOS_MID_5hr w/ State_Met%SUNCOSmid5
!! 14 Mar 2013 - M. Payer - Replace Ox with O3 as part of removal of
!! NOx-Ox partitioning
!! 30 May 2013 - R. Yantosca - Replace TCVV with Input_Opt%TCVV
!
!!EOP
!!------------------------------------------------------------------------------
!!BOC
!!
!! !LOCAL VARIABLES:
!!
! !=======================================================================
! ! temp : model temperature
! ! jno2 : J(NO2) value
! ! cao3 : concentration O3 in ambient air
! ! alfa0 : solar zenith angle 5 hours ago
! ! alfa5 : solar zenith angle at this time
! ! jo1d : ratio J(O1D)/J(NO2)
! ! caco : concentration CO in ambient air
! !=======================================================================
!
! INTEGER :: IJLOOP
! INTEGER, PARAMETER :: ntemp = 4
! INTEGER, PARAMETER :: njno2 = 4
! INTEGER, PARAMETER :: ncao3 = 4
! integer, PARAMETER :: nalfa0 = 12
! INTEGER, PARAMETER :: nalfa5 = 12
! INTEGER, PARAMETER :: njo1d = 4
! INTEGER, PARAMETER :: ncaco = 4
!
! REAL, DIMENSION(ntemp) :: templev
! REAL, DIMENSION(njno2) :: jno2lev
! REAL, DIMENSION(ncao3) :: cao3lev
! REAL, DIMENSION(nalfa0) :: alfa0lev
! REAL, DIMENSION(nalfa5) :: alfa5lev
! REAL, DIMENSION(njo1d) :: jo1dlev
! REAL, DIMENSION(ncaco) :: cacolev
!
! ! Temporary variable storage
! REAL :: temp_tmp, jno2_tmp, cao3_tmp
! REAL :: alfa0_tmp, alfa5_tmp, jo1d_tmp
! REAL :: caco_tmp
!
! ! Interpolation parameters
! REAL, DIMENSION(2) :: xtemp, xjno2, xcao3, xalfa0
! REAL, DIMENSION(2) :: xalfa5, xjo1d, xcaco
!
! ! For loops
! INTEGER :: itemp, ijno2, icao3, ialfa0
! INTEGER :: ialfa5, ijo1d, icaco
! INTEGER :: i0,i1,i2,i3,i4,i5,i6,i7
!
! ! array contain temp, jno2, cao3, alfa_0, alfa_5, jo1d, caco
! REAL,DIMENSION(7) :: var_array
!
! CHARACTER(LEN=255) :: MSG
!
! ! Pointers
! ! We need to define local arrays to hold corresponding values
! ! from the Chemistry State (State_Chm) object. (mpayer, 12/6/12)
! REAL*8, POINTER :: STT(:,:,:,:)
!
! !=================================================================
! ! INTERPOLATE_LUT begins here!
! !=================================================================
!
! stt => State_Chm%TRACERS
!
! ! Set the levels that were chosen in the look up table
! templev = (/ 275., 280., 285., 300. /)
! jno2lev = (/ 5.e-4, 0.0025, 0.0050, 0.012 /)
! cao3lev = (/ 5., 20., 35., 75. /)
! alfa0lev = (/ -90., -60., -45., -30.,
! $ -15., 0., 15., 30.,
! $ 45., 60., 75., 90. /)
! alfa5lev = (/ -90., -60., -45., -30.,
! $ -15., 0., 15., 30.,
! $ 45., 60., 75., 90. /)
! jo1dlev = (/ 5.e-4, 0.0015, 0.0025, 0.0055 /)
! cacolev = (/ 50., 100., 150., 1200. /)
!
!! PRINT*,"Temperature levels are: ",templev
!! PRINT*,"This is grid cell: ",I,J
!
! ! Temperature
!! PRINT*,"Temperature here is: ",State_Met%TS(I,J)
!! PRINT*,"USA: ", State_Met%TS(32,64)
!
! ! Tracer concentrations in v/v
!! PRINT*,"[O3] is: ",STT(I,J,1,IDTOX3)/ State_Met%AD(I,J,1) * TCVV(IDTO3)
!! PRINT*,"[CO] is: ",STT(I,J,1,IDTCO)/ State_Met%AD(I,J,1) * TCVV(IDTCO)
!! PRINT*,"IDTO3 is: ", IDTO3
!! PRINT*,"IDO3 is: ", IDO3
!! PRINT*,"In USA: ",STT(32,64,1,IDTO3)/State_Met%AD(32,64,1) * TCVV(IDTO3)
!
! ! SOLAR ZENITH ANGLES IN DEGREES
!! IJLOOP = ( (J-1) * IIPAR ) + I
!! PRINT*,"Local Time: ",GET_LOCALTIME(I,1,1)
!! PRINT*,"Solar Zenith Angle at this location: ",
!! $ ASIND(State_Met%SUNCOSmid(I,J))
!! IJLOOP = ( (64-1) * IIPAR ) + 32
!! PRINT*,"Local USA time: ", GET_LOCALTIME(32,1,1)
!! PRINT*,"Solar Zenith Angle at USA: ",
!! & ASIND(State_Met%SUNCOSmid(I,J))
!! PRINT*,"Solar Zenith Angle at USA - 5: ",
!! & ASIND(State_Met%SUNCOSmid5(I,J))
!
! ! Set the variables
! IJLOOP = ( (J-1) * IIPAR ) + I
! var_array(1) = State_Met%TS(I,J) ! Temperature
! var_array(2) = JNO2 ! J(NO2)
! var_array(3) = STT(I,J,1,IDTO3) / ! [O3] in ppbv
! $ State_Met%AD(I,J,1) *
! $ Input_OPt%TCVV(IDTO3) * 1.E9
! var_array(4) = ASIND(State_Met%SUNCOSmid5(I,J)) ! alfa0
! var_array(5) = ASIND(State_Met%SUNCOSmid(I,J)) ! alfa5
! var_array(6) = JO1D / JNO2 ! J(O1D)/J(NO2)
! var_array(7) = STT(I,J,1,IDTCO) / ! [CO] in ppbv
! $ State_Met%AD(I,J,1) *
! $ Input_Opt%TCVV(IDTCO) * 1.E9
!
! ! prevent NaN when jvalues are 0.
! IF (JNO2 .eq. 0.) var_array(6) = 0.
!
! ! First some error checking
! ! ########### MAYBE CHECK HERE FOR NEGATIVE VALUES?##########
!
! !
! ! Determine reference index ( itemp, ijno2, icao3, ialfa0,
! ! ialfa5, ialfa, ijo1d, icaco )
! !
! !=================================================================
! ! Find smallest temperature reference level (i) for which actual
! ! temperature is smaller, then do
! !
! ! x(1) = ( temperature_level(i+1) - actual temperature )
! ! -------------------------------------------------
! ! ( temperature_level(i+1) - temperature_level(i) )
! !
! ! then x(2) = 1.0 - x(1)
! !
! !=================================================================
!
! !----------------------
! ! Temperature:
! !----------------------
! temp_tmp = var_array(1)
!
! ! If temperature larger than largest in LUT, assign largest temp
! IF ( var_array(1) > templev( ntemp ) ) temp_tmp = templev(ntemp)
!
! ! If temp smaller, assign smallest temp level
! IF ( var_array(1) < templev(1) ) temp_tmp = templev(1)
!
! DO i0=1,ntemp-1
! itemp = i0
! IF( templev( itemp+1 ) > temp_tmp ) EXIT
! END DO
!
! xtemp(1) = ( templev( itemp+1 ) - temp_tmp ) /
! $ ( templev( itemp+1 ) - templev( itemp ) )
! xtemp(2) = 1.0 - xtemp(1)
!
! !----------------------
! ! J(NO2):
! !----------------------
! jno2_tmp = var_array(2)
!
! ! If larger than largest in LUT, assign largest level values
! IF ( var_array(2) > jno2lev( njno2 ) ) jno2_tmp = jno2lev(njno2)
!
! ! If smaller, assign smallest level value
! IF ( var_array(2) < jno2lev(1) ) jno2_tmp = jno2lev(1)
!
! DO i0=1,njno2-1
! ijno2 = i0
! IF( jno2lev( ijno2+1 ) > jno2_tmp ) EXIT
! END DO
!
! xjno2(1) = ( jno2lev( ijno2+1 ) - jno2_tmp ) /
! $ ( jno2lev( ijno2+1 ) - jno2lev( ijno2 ) )
! xjno2(2) = 1.0 - xjno2(1)
!
! !----------------------
! ! [O3]:
! !----------------------
! cao3_tmp = var_array(3)
!
! ! If larger than largest in LUT, assign largest level values
! IF ( var_array(3) > cao3lev( ncao3 ) ) cao3_tmp = cao3lev(ncao3)
!
! ! If smaller, assign smallest level value
! IF ( var_array(3) < cao3lev(1) ) cao3_tmp = cao3lev(1)
!
! DO i0=1,ncao3-1
! icao3 = i0
! IF( cao3lev( icao3+1 ) > cao3_tmp ) EXIT
! END DO
!
! xcao3(1) = ( cao3lev( icao3+1 ) - cao3_tmp ) /
! $ ( cao3lev( icao3+1 ) - cao3lev( icao3 ) )
! xcao3(2) = 1.0 - xcao3(1)
!
! !----------------------
! ! alfa0:
! !----------------------
! alfa0_tmp = var_array(4)
!
! ! If larger than largest in LUT, assign largest level values
! IF ( var_array(4) > alfa0lev( nalfa0 ) ) alfa0_tmp =
! $ alfa0lev(nalfa0)
!
! ! If smaller, assign smallest level value
! IF ( var_array(4) < alfa0lev(1) ) alfa0_tmp = alfa0lev(1)
!
! DO i0=1,nalfa0-1
! ialfa0 = i0
! IF( alfa0lev( ialfa0+1 ) > alfa0_tmp ) EXIT
! END DO
!
! xalfa0(1) = ( alfa0lev( ialfa0+1 ) - alfa0_tmp ) /
! $ ( alfa0lev( ialfa0+1 ) - alfa0lev( ialfa0 ) )
! xalfa0(2) = 1.0 - xalfa0(1)
!
! !----------------------
! ! alfa5:
! !----------------------
! alfa5_tmp = var_array(5)
!
! ! If larger than largest in LUT, assign largest level values
! IF ( var_array(5) > alfa5lev( nalfa5 ) ) alfa5_tmp =
! $ alfa5lev(nalfa5)
!
! ! If smaller, assign smallest level value
! IF ( var_array(5) < alfa5lev(1) ) alfa5_tmp = alfa5lev(1)
!
! DO i0=1,nalfa5-1
! ialfa5 = i0
! IF( alfa5lev( ialfa5+1 ) > alfa5_tmp ) EXIT
! END DO
!
! xalfa5(1) = ( alfa5lev( ialfa5+1 ) - alfa5_tmp ) /
! $ ( alfa5lev( ialfa5+1 ) - alfa5lev( ialfa5 ) )
! xalfa5(2) = 1.0 - xalfa5(1)
!
! !----------------------
! ! jo1d:
! !----------------------
! jo1d_tmp = var_array(6)
!
! ! If larger than largest in LUT, assign largest level values
! IF ( var_array(6) > jo1dlev( njo1d ) ) jo1d_tmp = jo1dlev(njo1d)
!
! ! If smaller, assign smallest level value
! IF ( var_array(6) < jo1dlev(1) ) jo1d_tmp = jo1dlev(1)
!
! DO i0=1,njo1d-1
! ijo1d = i0
! IF( jo1dlev( ijo1d+1 ) > jo1d_tmp ) EXIT
! END DO
!
! xjo1d(1) = ( jo1dlev( ijo1d+1 ) - jo1d_tmp ) /
! $ ( jo1dlev( ijo1d+1 ) - jo1dlev( ijo1d ) )
! xjo1d(2) = 1.0 - xjo1d(1)
!
! !----------------------
! ! [CO]:
! !----------------------
! caco_tmp = var_array(7)
!
! ! If larger than largest in LUT, assign largest level values
! IF ( var_array(7) > cacolev( ncaco ) ) caco_tmp = cacolev(ncaco)
!
! ! If smaller, assign smallest level value
! IF ( var_array(7) < cacolev(1) ) caco_tmp = cacolev(1)
!
! DO i0=1,ncaco-1
! icaco = i0
! IF( cacolev( icaco+1 ) > caco_tmp ) EXIT
! END DO
!
! xcaco(1) = ( cacolev( icaco+1 ) - caco_tmp ) /
! $ ( cacolev( icaco+1 ) - cacolev( icaco ) )
! xcaco(2) = 1.0 - xcaco(1)
!
!! PRINT*,"The i-values are:", itemp, ijno2, icao3, ialfa0,
!! $ ialfa5, ijo1d, icaco
!! PRINT*,"Variables are: ", var_array
!! PRINT*,"For testing, xtemp: ", xtemp
!
! !========================
! ! Linear interpolation
! !========================
!
! fraction_nox = 0.0
! int_ope = 0.0
!
! DO i1=1,2
! DO i2=1,2
! DO i3=1,2
! DO i4=1,2
! DO i5=1,2
! DO i6=1,2
! DO i7=1,2
!
! !IF ENCOUNTER -999 IN THE LUT PRINT ERROR!!
! IF ( ( fracnox( itemp+i1-1, ijno2+i2-1, icao3+i3-1,
! $ ialfa0+i4-1, ialfa5+i5-1, ijo1d+i6-1,
! $ icaco+i7-1 ) < 0. ) .or.
! $ ( fracnox( itemp+i1-1, ijno2+i2-1, icao3+i3-1,
! $ ialfa0+i4-1, ialfa5+i5-1, ijo1d+i6-1,
! $ icaco+i7-1 ) > 1. ) ) THEN
!
!! PRINT*, 'INTERPOLATE_LUT: fracnox = ,',
!! $ fracnox( itemp+i1-1, ijno2+i2-1, icao3+i3-1,
!! $ ialfa0+i4-1, ialfa5+i5-1, ijo1d+i6-1,
!! $ icaco+i7-1 )
!
! MSG = 'LUT error: Fracnox should be between 0 and 1!'
! CALL ERROR_STOP( MSG, 'INTERPOLATE_LUT ("paranox_mod.F")' )
! ENDIF
!
! !Cycle if both angles are 0
! IF ((ialfa0+i4-1 .eq. 6) .and. (ialfa5+i5-1 .eq. 6) )
! $ CYCLE
!
! ! fracnox is the array with the actual lut data
! fraction_nox = fraction_nox + xtemp(i1) * xjno2(i2) *
! $ xcao3(i3) * xalfa0(i4) * xalfa5(i5) *
! $ xjo1d(i6) * xcaco(i7) *
! $ fracnox( itemp+i1-1, ijno2+i2-1,
! $ icao3+i3-1, ialfa0+i4-1, ialfa5+i5-1,
! $ ijo1d+i6-1, icaco+i7-1 )
!
! ! intope is the array with the actual lut data
! int_ope = int_ope + xtemp(i1) * xjno2(i2) *
! $ xcao3(i3) * xalfa0(i4) * xalfa5(i5) *
! $ xjo1d(i6) * xcaco(i7) *
! $ intope( itemp+i1-1, ijno2+i2-1,
! $ icao3+i3-1, ialfa0+i4-1, ialfa5+i5-1,
! $ ijo1d+i6-1, icaco+i7-1 )
!
! END DO
! END DO
! END DO
! END DO
! END DO
! END DO
! END DO
!!
!! IF ((I .eq. 108) .and. (J .eq. 49)) THEN
!! PRINT*,"----INTERPOLATE_LUT-----"
!! PRINT*,"fraction_nox and int_OPE: ",fraction_nox,
!! $ int_ope
!! PRINT*,"Jvalues are: ", jvalues(I, J,:)
!! PRINT*,"Vars are: ", var_array
!! PRINT*,"[O3] in interpolate_lut: ", var_array(3)
!! PRINT*,"[CO] in interpolate_lut: ", var_array(7)
!! PRINT*,"The i-values are:", itemp, ijno2, icao3,
!! $ ialfa0, ialfa5, ijo1d,
!! $ icaco
!! ENDIF
!!
!! PRINT*,"fraction_nox is: ",fraction_nox
!! PRINT*,"integrated OPE: ",int_ope
!
! NULLIFY( STT )
!
! END SUBROUTINE INTERPOLATE_LUT
!!EOC
!##############################################################################
!------------------------------------------------------------------------------
! Harvard University Atmospheric Chemistry Modeling Group !
!------------------------------------------------------------------------------
!BOP
!
! !IROUTINE: interpolate_lut2
!
! !DESCRIPTION: Subroutine INTERPOLATE\_LUT2 returns FracNOx or IntOPE from
! the lookup tables (G.C.M. Vinken, KNMI, June 2010)
!\\
!\\
! !INTERFACE:
!
SUBROUTINE INTERPOLATE_LUT2( I, J, o3, no, no2, dens, JO1D, JNO2,
& fraction_nox, int_ope )
!
! !USES:
!
USE TIME_MOD, ONLY : GET_LOCALTIME
USE ERROR_MOD, ONLY : ERROR_STOP, SAFE_DIV
USE DAO_MOD, ONLY : TS, SUNCOS, SUNCOS_5hr
# include "CMN_SIZE" ! Size parameters
!
! !INPUT PARAMETERS:
!
INTEGER, INTENT(IN) :: I, J
REAL*8, INTENT(IN) :: o3, no, no2, dens, JNO2, JO1D
!
! !OUTPUT PARAMETERS:
!
REAL*4, INTENT(OUT) :: fraction_nox, int_ope
!
! !REVISION HISTORY:
! Jun 2010 - G.C.M. Vinken - Initial version
! 21 Feb 2011 - G.C.M. Vinken - Updated for NOx in LUT
! 06 Feb 2012 - M. Payer - Moved from emissions_mod.F to paranox_mod.F;
! Added ProTeX headers
! 15 Feb 2012 - M. Payer - Add error trap to ensure 0 < fracnox < 1.
! 09 Nov 2012 - M. Payer - Replaced all met field arrays with State_Met
! derived type object
! 28 Nov 2012 - R. Yantosca - Replace SUNCOS_MID w/ State_Met%SUNCOSmid
! 28 Nov 2012 - R. Yantosca - Replace SUNCOS_MID_5hr w/ State_Met%SUNCOSmid5
! 17 Jun 2013 - R. Yantosca - Bug fix: declare all REAL variables with
! REAL*4 in order to avoid numerical precision
! errors when compiling with OMP=yes.
!EOP
!------------------------------------------------------------------------------
!BOC
!
! !LOCAL VARIABLES:
!
!=======================================================================
! temp : model temperature
! jno2 : J(NO2) value
! cao3 : concentration O3 in ambient air
! alfa0 : solar zenith angle 5 hours ago
! alfa5 : solar zenith angle at this time
! jo1d : ratio J(O1D)/J(NO2)
! canox : concentration NOx in ambient air
!
! o3 : incoming o3 concentration
! no : incoming no
! no2 : incoming no2
! dens : incoming air density
!=======================================================================
INTEGER :: IJLOOP
INTEGER, PARAMETER :: ntemp = 4
INTEGER, PARAMETER :: njno2 = 4
INTEGER, PARAMETER :: ncao3 = 4
INTEGER, PARAMETER :: nalfa0 = 12
INTEGER, PARAMETER :: nalfa5 = 12
INTEGER, PARAMETER :: njo1d = 4
INTEGER, PARAMETER :: ncanox = 5
REAL*4, DIMENSION(ntemp) :: templev
REAL*4, DIMENSION(njno2) :: jno2lev
REAL*4, DIMENSION(ncao3) :: cao3lev
REAL*4, DIMENSION(nalfa0) :: alfa0lev
REAL*4, DIMENSION(nalfa5) :: alfa5lev
REAL*4, DIMENSION(njo1d) :: jo1dlev
REAL*4, DIMENSION(ncanox) :: canoxlev
! Temporary variable storage
REAL*4 :: temp_tmp, jno2_tmp, cao3_tmp
REAL*4 :: alfa0_tmp, alfa5_tmp, jo1d_tmp
REAL*4 :: canox_tmp
! Interpolation parameters
REAL*4, DIMENSION(2) :: xtemp, xjno2, xcao3, xalfa0
REAL*4, DIMENSION(2) :: xalfa5, xjo1d, xcanox
! For loops
INTEGER :: itemp, ijno2, icao3, ialfa0
INTEGER :: ialfa5, ijo1d, icanox
INTEGER :: i0,i1,i2,i3,i4,i5,i6,i7
! array contain temp, jno2, cao3, alfa_0, alfa_5, jo1d, canox
REAL*4, DIMENSION(7) :: var_array
CHARACTER(LEN=255) :: MSG
!=================================================================
! INTERPOLATE_LUT2 begins here!
!=================================================================
! Set the levels that were chosen in the look up table
templev = (/ 275., 280., 285., 310. /)
jno2lev = (/ 5.e-4, 0.0025, 0.0050, 0.012 /)
cao3lev = (/ 5., 20., 35., 75. /)
alfa0lev = (/ -90., -60., -45., -30.,
$ -15., 0., 15., 30.,
$ 45., 60., 75., 90. /)
alfa5lev = (/ -90., -60., -45., -30.,
$ -15., 0., 15., 30.,
$ 45., 60., 75., 90. /)
jo1dlev = (/ 5.e-4, 0.0015, 0.0025, 0.0055 /)
canoxlev = (/ 10., 200., 1000., 2000., 6000. /)
! PRINT*,"Temperature levels are: ",templev
! PRINT*,"This is grid cell: ",I,J
! Temperature
! PRINT*,"Temperature here is: ",TS(I,J)
! PRINT*,"USA: ",TS(32,64)
! Tracer concentrations in v/v
! PRINT*,"[O3] is: ",STT(I,J,1,IDTO3)/ State_Met%AD(I,J,1) * TCVV(IDTO3)
! PRINT*,"[CO] is: ",STT(I,J,1,IDTCO)/ State_Met%AD(I,J,1) * TCVV(IDTCO)
! PRINT*,"IDTO3 is: ", IDTO3
! PRINT*,"IDO3 is: ", IDO3
! PRINT*,"In USA: ",STT(32,64,1,IDTO3)/State_Met%AD(32,64,1) * TCVV(IDTO3)
! SOLAR ZENITH ANGLES IN DEGREES
! IJLOOP = ( (J-1) * IIPAR ) + I
! PRINT*,"Local Time: ",GET_LOCALTIME(I)
! PRINT*,"Solar Zenith Angle at this location: ",
! $ ASIND(SUNCOS(IJLOOP))
! IJLOOP = ( (64-1) * IIPAR ) + 32
! PRINT*,"Local USA time: ", GET_LOCALTIME(32)
! PRINT*,"Solar Zenith Angle at USA: ",
! & ASIND(SUNCOS(I,J))
! PRINT*,"Solar Zenith Angle at USA - 5: ",
! & ASIND(SUNCOS_5hr(IJLOOP))
! Set the variables
IJLOOP = ( (J-1) * IIPAR ) + I
var_array(1) = TS(I,J) ! Temperature
var_array(2) = JNO2 ! J(NO2), 1/s
var_array(3) = o3 / dens * 1.E9 ! [O3] in ppbv
var_array(4) = ASIND(SUNCOS(IJLOOP)) ! alfa0
var_array(5) = ASIND(SUNCOS_5hr(IJLOOP)) ! alfa5
var_array(6) = SAFE_DIV( JO1D, JNO2, 0d0 ) ! J(O1D)/J(NO2)
var_array(7) = (no + no2) / dens * 1.E12 ! [NOx] in pptv
! prevent NaN when jvalues are 0.
IF (JNO2 .eq. 0.) var_array(6) = 0.
! First some error checking
! ########### MAYBE CHECK HERE FOR NEGATIVE VALUES?##########
!
! Determine reference index ( itemp, ijno2, icao3, ialfa0,
! ialfa5, ialfa, ijo1d, icaco )
!
!========================================================================
! Find smallest temperature reference level (i) for which actual
! temperature is smaller, then do
!
! x(1) = ( temperature_level(i+1) - actual temperature )
! -------------------------------------------------
! ( temperature_level(i+1) - temperature_level(i) )
!
! then x(2) = 1.0 - x(1)
!
!========================================================================
!---------------------
! Temperature:
!---------------------
temp_tmp = var_array(1)
! If temperature larger than largest in LUT, assign largest temp
IF ( var_array(1) > templev( ntemp ) ) temp_tmp = templev(ntemp)
! If temp smaller, assign smallest temp level
IF ( var_array(1) < templev(1) ) temp_tmp = templev(1)
DO i0=1,ntemp-1
itemp = i0
IF( templev( itemp+1 ) > temp_tmp ) EXIT
END DO
xtemp(1) = ( templev( itemp+1 ) - temp_tmp ) /
$ ( templev( itemp+1 ) - templev( itemp ) )
xtemp(2) = 1.0 - xtemp(1)
!---------------------
! J(NO2):
!---------------------
jno2_tmp = var_array(2)
! If larger than largest in LUT, assign largest level values
IF ( var_array(2) > jno2lev( njno2 ) ) jno2_tmp = jno2lev(njno2)
! If smaller, assign smallest level value
IF ( var_array(2) < jno2lev(1) ) jno2_tmp = jno2lev(1)
DO i0=1,njno2-1
ijno2 = i0
IF( jno2lev( ijno2+1 ) > jno2_tmp ) EXIT
END DO
xjno2(1) = ( jno2lev( ijno2+1 ) - jno2_tmp ) /
$ ( jno2lev( ijno2+1 ) - jno2lev( ijno2 ) )
xjno2(2) = 1.0 - xjno2(1)
!---------------------
! [O3]:
!---------------------
cao3_tmp = var_array(3)
! If larger than largest in LUT, assign largest level values
IF ( var_array(3) > cao3lev( ncao3 ) ) cao3_tmp = cao3lev(ncao3)
! If smaller, assign smallest level value
IF ( var_array(3) < cao3lev(1) ) cao3_tmp = cao3lev(1)
DO i0=1,ncao3-1
icao3 = i0
IF( cao3lev( icao3+1 ) > cao3_tmp ) EXIT
END DO
xcao3(1) = ( cao3lev( icao3+1 ) - cao3_tmp ) /
$ ( cao3lev( icao3+1 ) - cao3lev( icao3 ) )
xcao3(2) = 1.0 - xcao3(1)
!---------------------
! alfa0:
!---------------------
alfa0_tmp = var_array(4)
! If larger than largest in LUT, assign largest level values
IF ( var_array(4) > alfa0lev( nalfa0 ) ) alfa0_tmp =
$ alfa0lev(nalfa0)
! If smaller, assign smallest level value
IF ( var_array(4) < alfa0lev(1) ) alfa0_tmp = alfa0lev(1)
DO i0=1,nalfa0-1
ialfa0 = i0
IF( alfa0lev( ialfa0+1 ) > alfa0_tmp ) EXIT
END DO
xalfa0(1) = ( alfa0lev( ialfa0+1 ) - alfa0_tmp ) /
$ ( alfa0lev( ialfa0+1 ) - alfa0lev( ialfa0 ) )
xalfa0(2) = 1.0 - xalfa0(1)
!---------------------
! alfa5:
!---------------------
alfa5_tmp = var_array(5)
! If larger than largest in LUT, assign largest level values
IF ( var_array(5) > alfa5lev( nalfa5 ) ) alfa5_tmp =
$ alfa5lev(nalfa5)
! If smaller, assign smallest level value
IF ( var_array(5) < alfa5lev(1) ) alfa5_tmp = alfa5lev(1)
DO i0=1,nalfa5-1
ialfa5 = i0
IF( alfa5lev( ialfa5+1 ) > alfa5_tmp ) EXIT
END DO
xalfa5(1) = ( alfa5lev( ialfa5+1 ) - alfa5_tmp ) /
$ ( alfa5lev( ialfa5+1 ) - alfa5lev( ialfa5 ) )
xalfa5(2) = 1.0 - xalfa5(1)
!---------------------
! jo1d:
!---------------------
jo1d_tmp = var_array(6)
! If larger than largest in LUT, assign largest level values
IF ( var_array(6) > jo1dlev( njo1d ) ) jo1d_tmp = jo1dlev(njo1d)
! If smaller, assign smallest level value
IF ( var_array(6) < jo1dlev(1) ) jo1d_tmp = jo1dlev(1)
DO i0=1,njo1d-1
ijo1d = i0
IF( jo1dlev( ijo1d+1 ) > jo1d_tmp ) EXIT
END DO
xjo1d(1) = ( jo1dlev( ijo1d+1 ) - jo1d_tmp ) /
$ ( jo1dlev( ijo1d+1 ) - jo1dlev( ijo1d ) )
xjo1d(2) = 1.0 - xjo1d(1)
!---------------------
! [NOx]:
!---------------------
canox_tmp = var_array(7)
! If larger than largest in LUT, assign largest level values
IF ( var_array(7) > canoxlev( ncanox ) ) canox_tmp =
$ canoxlev(ncanox)
! If smaller, assign smallest level value
IF ( var_array(7) < canoxlev(1) ) canox_tmp = canoxlev(1)
DO i0=1,ncanox-1
icanox = i0
IF( canoxlev( icanox+1 ) > canox_tmp ) EXIT
END DO
xcanox(1) = ( canoxlev( icanox+1 ) - canox_tmp ) /
$ ( canoxlev( icanox+1 ) - canoxlev( icanox ) )
xcanox(2) = 1.0 - xcanox(1)
! PRINT*,"The i-values are:", itemp, ijno2, icao3, ialfa0,
! $ ialfa5, ijo1d, icanox
! PRINT*,"Variables are: ", var_array
! PRINT*,"For testing, xtemp: ", xtemp
!======================
! Linear interpolation
!======================
fraction_nox = 0.0
int_ope = 0.0
DO i1=1,2
DO i2=1,2
DO i3=1,2
DO i4=1,2
DO i5=1,2
DO i6=1,2
DO i7=1,2
!IF ENCOUNTER -999 IN THE LUT PRINT ERROR!!
IF ( ( fracnox( itemp+i1-1, ijno2+i2-1, icao3+i3-1,
$ ialfa0+i4-1, ialfa5+i5-1, ijo1d+i6-1,
$ icanox+i7-1 ) < 0. ) .or.
$ ( fracnox( itemp+i1-1, ijno2+i2-1, icao3+i3-1,
$ ialfa0+i4-1, ialfa5+i5-1, ijo1d+i6-1,
$ icanox+i7-1 ) > 1. ) ) THEN
PRINT*, 'INTERPOLATE_LUT2: fracnox = ',
$ fracnox( itemp+i1-1, ijno2+i2-1, icao3+i3-1,
$ ialfa0+i4-1, ialfa5+i5-1, ijo1d+i6-1,
$ icanox+i7-1 )
MSG = 'LUT error: Fracnox should be between 0 and 1!'
CALL ERROR_STOP( MSG, 'INTERPOLATE_LUT2 ("paranox_mod.F")' )
ENDIF
! fracnox is the array with the actual lut data
fraction_nox = fraction_nox + xtemp(i1) * xjno2(i2) *
$ xcao3(i3) * xalfa0(i4) * xalfa5(i5) *
$ xjo1d(i6) * xcanox(i7) *
$ fracnox( itemp+i1-1, ijno2+i2-1,
$ icao3+i3-1, ialfa0+i4-1, ialfa5+i5-1,
$ ijo1d+i6-1, icanox+i7-1 )
! intope is the array with the actual lut data
int_ope = int_ope + xtemp(i1) * xjno2(i2) *
$ xcao3(i3) * xalfa0(i4) * xalfa5(i5) *
$ xjo1d(i6) * xcanox(i7) *
$ intope( itemp+i1-1, ijno2+i2-1,
$ icao3+i3-1, ialfa0+i4-1, ialfa5+i5-1,
$ ijo1d+i6-1, icanox+i7-1 )
END DO
END DO
END DO
END DO
END DO
END DO
END DO
! IF ((I .eq. 108) .and. (J .eq. 49)) THEN
! PRINT*,"-----INTERPOLATE_LUT2, for 108,49-----"
! PRINT*,"Fraction_nox and int_OPE: ", fraction_nox,
! & int_ope
! PRINT*,"Jvalues are: ", JNO2, JO1D
! PRINT*,"Vars are: ", var_array
! PRINT*,"[O3] in interpolate_lut: ", var_array(3)
! PRINT*,"[NOx] in interpolate_lut: ", var_array(7)
! PRINT*,"J(O1D)/J(NO2) : ", var_array(6)
! PRINT*,"The i-values are:", itemp, ijno2, icao3,
! $ ialfa0, ialfa5, ijo1d,
! $ icanox
! PRINT*,"Interpolation parameters: ", xtemp, xjno2, xcao3,
! $ xalfa0, xalfa5, xjo1d,
! $ xcanox
! PRINT*,"---------------------------------"
! ENDIF
!
! IF ((I .eq. 73) .and. (J .eq. 76)) THEN
! PRINT*,"-----INTERPOLATE_LUT2, for 73,76-----"
! PRINT*,"Fraction_nox and int_OPE: ", fraction_nox,
! & int_ope
! PRINT*,"Jvalues are: ", JNO2, JO1D
! PRINT*,"Vars are: ", var_array
! PRINT*,"[O3] in interpolate_lut: ", var_array(3)
! PRINT*,"[NOx] in interpolate_lut: ", var_array(7)
! PRINT*,"J(O1D)/J(NO2) : ", var_array(6)
! PRINT*,"The i-values are:", itemp, ijno2, icao3,
! $ ialfa0, ialfa5, ijo1d,
! $ icanox
! PRINT*,"Interpolation parameters: ", xtemp, xjno2, xcao3,
! $ xalfa0, xalfa5, xjo1d,
! $ xcanox
! PRINT*,"------------------------------------"
! ENDIF
END SUBROUTINE INTERPOLATE_LUT2
!EOC
END MODULE PARANOX_MOD