!$Id: strat_chem_mod.f,v 1.2 2012/07/13 20:09:14 nicolas Exp $ !------------------------------------------------------------------------------ ! Harvard University Atmospheric Chemistry Modeling Group ! !------------------------------------------------------------------------------ !BOP ! ! !MODULE: strat_chem_mod ! ! !DESCRIPTION: Module STRAT\_CHEM\_MOD contains variables and routines for ! performing a simple linearized chemistry scheme for more realistic ! upper boundary conditions. Archived 3D monthly climatological production ! rates and loss frequencies are applied from the GMI combo model. ! ! In the original schem code (schem.f), only the following species ! were destroyed by photolysis in the stratosphere: ! PAN, H2O2, ACET, MEK, ALD2, RCHO, MVK, MACR, R4N2, CH2O, N2O5, HNO4, MP ! and by reaction with OH for: ! ALK4, ISOP, H2O2, ACET, MEK, ALD2, RCHO, MVK, MACR, PMN, R4N2, ! PRPE, C3H8, CH2O, C2H6, HNO4, MP ! ! The updated code includes at least all of these, and many more. The code ! is flexible enough to automatically apply the rate to any new tracers ! for future simulations that share the name in tracer_mod with the ! GMI name. (See Documentation). ! !\\ !\\ ! !INTERFACE: ! MODULE STRAT_CHEM_MOD ! ! !USES: ! IMPLICIT NONE PRIVATE ! ! !PUBLIC MEMBER FUNCTIONS: ! PUBLIC :: DO_STRAT_CHEM PUBLIC :: CLEANUP_STRAT_CHEM ! hml PUBLIC :: PROD_0 PUBLIC :: LOSS_0 PUBLIC :: PROD PUBLIC :: LOSS PUBLIC :: DTCHEM PUBLIC :: NSCHEM PUBLIC :: GET_RATES PUBLIC :: GET_RATES_INTERP PUBLIC :: Strat_TrID_GC ! ! !PRIVATE MEMBER FUNCTIONS: ! PRIVATE :: INIT_STRAT_CHEM ! ! !ADJOINT GROUP: ! ! ! !PUBLIC DATA MEMBERS: ! ! !REMARKS: ! ! References: ! ============================================================================ ! (1 ) ! !REVISION HISTORY: ! 1 Feb 2011 - L. Murray - Initial version ! 22 Oct 2011 - H.-M. Lee - Modified to implement in adjoint. ! Now we can calculte strat prod and loss sensitivity. adj32_025 !EOP !------------------------------------------------------------------------------ !BOC ! ! !PRIVATE TYPES: ! ! Scalars REAL*8 :: dTchem ! chemistry time step [s] INTEGER, PARAMETER :: NTR_GMI = 120 ! Number of species from GMI model INTEGER :: NSCHEM ! Number of species upon which to ! apply P's & k's in GEOS-Chem ! Arrays REAL*8, ALLOCATABLE :: PROD(:,:,:,:) REAL*8, ALLOCATABLE :: LOSS(:,:,:,:) REAL*8, ALLOCATABLE :: STRAT_OH(:,:,:) ! Monthly mean OH [v/v] INTEGER, SAVE :: ncID_strat_rates ! For adjoint REAL*8, ALLOCATABLE :: PROD_0(:,:,:,:) REAL*8, ALLOCATABLE :: LOSS_0(:,:,:,:) CHARACTER(LEN=16) :: GMI_TrName(NTR_GMI) ! Tracer names in GMI INTEGER :: Strat_TrID_GC(NTR_GMI) ! Maps 1:NSCHEM to STT index INTEGER :: Strat_TrID_GMI(NTR_GMI) ! Maps 1:NSCHEM to GMI index ! (At most NTR_GMI species could overlap between G-C & GMI) ! Variables used to calculate the strat-trop exchange flux REAL*8 :: TauInit ! Initial time INTEGER :: NymdInit, NhmsInit ! Initial date REAL*8 :: TpauseL_Cnt ! Tropopause counter REAL*8, ALLOCATABLE :: TpauseL(:,:) ! Tropopause level aggregator REAL*8, ALLOCATABLE :: MInit(:,:,:,:) ! Init. atm. state for STE period REAL*8, ALLOCATABLE :: Before(:,:,:) ! Init. atm. state each chem. dt REAL*8, ALLOCATABLE :: SChem_Tend(:,:,:,:) ! Stratospheric chemical tendency ! (total P - L) [kg period-1] !================================================================= ! MODULE ROUTINES -- follow below the "CONTAINS" statement !================================================================= CONTAINS !EOC !------------------------------------------------------------------------------ ! Harvard University Atmospheric Chemistry Modeling Group ! !------------------------------------------------------------------------------ !BOP ! ! !IROUTINE: DO_STRAT_CHEM ! ! !DESCRIPTION: Function DO\_STRAT\_CHEM is the driver routine for computing ! the simple linearized stratospheric chemistry scheme for a host of species ! whose prod/loss rates were determined from the GMI combo model. Ozone is ! treated using either Linoz or Synoz. ! ! !INTERFACE: ! SUBROUTINE DO_STRAT_CHEM ! ! !USES: ! USE DAO_MOD, ONLY : AD, CONVERT_UNITS USE ERROR_MOD, ONLY : DEBUG_MSG, GEOS_CHEM_STOP USE LOGICAL_MOD, ONLY : LLINOZ, LPRT USE LINOZ_MOD, ONLY : DO_LINOZ USE TIME_MOD, ONLY : GET_MONTH, TIMESTAMP_STRING USE TRACER_MOD, ONLY : ITS_A_FULLCHEM_SIM, ITS_A_TAGOX_SIM USE TRACER_MOD, ONLY : N_TRACERS, STT, TCVV, TRACER_MW_KG USE TRACERID_MOD, ONLY : IDTOX USE TROPOPAUSE_MOD, ONLY : GET_MIN_TPAUSE_LEVEL, GET_TPAUSE_LEVEL USE TROPOPAUSE_MOD, ONLY : ITS_IN_THE_TROP USE NETCDF_UTIL_MOD ! adj_group (hml, 07/25/11) USE ADJ_ARRAYS_MOD, ONLY : PROD_SF, LOSS_SF USE ADJ_ARRAYS_MOD, ONLY : ID_LOSS USE ADJ_ARRAYS_MOD, ONLY : IFD, JFD, LFD, DO_CHK_FILE USE ADJ_ARRAYS_MOD, ONLY : NSTPL USE TRACER_MOD, ONLY : STT_STRAT_TMP USE LOGICAL_ADJ_MOD,ONLY : LADJ USE LOGICAL_ADJ_MOD,ONLY : LADJ_STRAT USE CHECKPOINT_MOD, ONLY : MAKE_BEFSTRAT_CHKFILE USE TIME_MOD, ONLY : GET_NHMS USE TIME_MOD, ONLY : GET_NYMD USE TIME_MOD, ONLY : GET_TAU USE UPBDFLX_MOD, ONLY : UPBDFLX_O3, INIT_UPBDFLX # include "CMN_SIZE" ! ! !REMARKS: ! ! !REVISION HISTORY: ! 1 Feb 2011 - L. Murray - Initial version !EOP !------------------------------------------------------------------------------ !BOC ! ! !LOCAL VARIABLES: ! LOGICAL, SAVE :: FIRST = .TRUE. INTEGER, SAVE :: LASTMONTH = -999 INTEGER, SAVE :: LASTSEASON = -1 INTEGER :: I, J, L, N, LMIN INTEGER :: IORD, JORD, KORD INTEGER :: NHMS INTEGER :: NYMD INTEGER :: NN, NS, NSL REAL*8 :: TAU REAL*8 :: dt, P, k, M0 REAL*8 :: STT0(IIPAR,JJPAR,LLPAR,N_TRACERS) CHARACTER(LEN=16) :: STAMP !=============================== ! DO_STRAT_CHEM begins here! !=============================== STAMP = TIMESTAMP_STRING() WRITE( 6, 10 ) STAMP 10 FORMAT( ' - DO_STRAT_CHEM: Linearized strat chemistry at ', a ) IF ( FIRST ) THEN ! Allocate all module arrays CALL INIT_STRAT_CHEM #if defined( GEOS_3 ) ! Initialize some Synoz variables IF ( .NOT. ( LLINOZ ) ) THEN CALL GET_ORD( IORD, JORD, KORD ) CALL INIT_UPBDFLX( IORD, JORD, KORD ) ENDIF #endif ENDIF ! Get the minimum level extent of the tropopause LMIN = GET_MIN_TPAUSE_LEVEL() IF ( GET_MONTH() /= LASTMONTH ) THEN IF ( LPRT ) CALL DEBUG_MSG( '### STRAT_CHEM: at GET_RATES' ) ! Read rates for this month IF ( ITS_A_FULLCHEM_SIM() ) THEN #if defined( GRID4x5 ) || defined( GRID2x25 ) CALL GET_RATES( GET_MONTH() ) #else ! For resolutions finer than 2x2.5, nested, ! or otherwise exotic domains and resolutions CALL GET_RATES_INTERP( GET_MONTH() ) #endif ENDIF ! Save month for next iteration LASTMONTH = GET_MONTH() ENDIF ! Set first-time flag to false FIRST = .FALSE. IF ( LPRT ) CALL DEBUG_MSG( '### STRAT_CHEM: at DO_STRAT_CHEM' ) WRITE(6,*) '-----------------------------------------------------' write(6,*) ' Doing stratospheric chemistry (STRAT_CHEM_MOD) ' WRITE(6,*) '-----------------------------------------------------' !================================================================ ! Full chemistry simulations !================================================================ IF ( ITS_A_FULLCHEM_SIM() ) THEN !! Advance counter for number of times we've sampled the tropopause level !TpauseL_CNT = TpauseL_CNT + 1d0 !============================================================= ! Do chemical production and loss for non-ozone species for ! which we have explicit prod/loss rates from GMI !============================================================= !$OMP PARALLEL DO !$OMP+DEFAULT( SHARED ) !$OMP+PRIVATE( I, J, L, N, k, P, dt, M0, NS, NN ) !$OMP+SCHEDULE( DYNAMIC ) DO J = 1, JJPAR DO I = 1, IIPAR ! Add to tropopause level aggregator for later determining STE flux TpauseL(I,J) = TpauseL(I,J) + GET_TPAUSE_LEVEL(I,J) DO L = LMIN, LLPAR IF ( ITS_IN_THE_TROP( I, J, L ) ) CYCLE DO N = 1, NSCHEM ! Tracer index of active strat chem species NN = Strat_TrID_GC(N) ! Tracer index in STT ! Skip Ox; we'll always use either Linoz or Synoz ! Adj use GMI rate for Ox if LINOZ is off (hml, 10/31/11) !IF ( ITS_A_FULLCHEM_SIM() .and. NN .eq. IDTOx ) CYCLE IF ( ITS_A_FULLCHEM_SIM() .and. (NN .eq. IDTOx) .and. & LLINOZ ) CYCLE ! adj_group: make a version that applies scaling factors ! and use this if the stratosphere adjoint ID #'s are active IF ( LADJ_STRAT ) THEN DO NS = 1, NSTPL NSL = ID_LOSS(NS) ! same for ID_PROD(NS) IF ( NN .EQ. NSL ) THEN PROD(I,J,L,N) = PROD_0(I,J,L,N) & * PROD_SF(I,J,1,NS) LOSS(I,J,L,N) = LOSS_0(I,J,L,N) & * LOSS_SF(I,J,1,NS) ENDIF ENDDO ENDIF ! Check point values of STT STT_STRAT_TMP(I,J,L,NN) = STT(I,J,L,NN) dt = DTCHEM ! timestep [s] k = LOSS(I,J,L,N) ! loss freq [s-1] P = PROD(I,J,L,N) * AD(I,J,L) / TCVV(NN)! production term [kg s-1] M0 = STT(I,J,L,NN) ! initial mass [kg] ! debug test !IF ( I == IFD .and. J == JFD .and. L == LFD ) THEN ! print*, NN,' STRAT TEST fwd: k = ', k ! print*, NN,' STRAT TEST fwd: P = ', P ! print*, NN,' STRAT TEST fwd: M0= ', M0 !ENDIF ! No prod or loss at all IF ( k .eq. 0d0 .and. P .eq. 0d0 ) CYCLE ! Simple analytic solution to dM/dt = P - kM over [0,dt] IF ( k .gt. 0d0 ) THEN STT(I,J,L,NN) = M0 * exp(-k*dt) + (P/k) & * (1d0-exp(-k*dt)) ELSE STT(I,J,L,NN) = M0 + P*dt ENDIF ! Aggregate chemical tendency [kg box-1] SCHEM_TEND(I,J,L,NN) = SCHEM_TEND(I,J,L,NN) & + ( STT(I,J,L,NN) - M0 ) ENDDO ! N ENDDO ! L ENDDO ! I ENDDO ! J !$OMP END PARALLEL DO ! Make check point file IF ( DO_CHK_FILE() ) THEN NHMS = GET_NHMS() NYMD = GET_NYMD() TAU = GET_TAU() CALL MAKE_BEFSTRAT_CHKFILE( NYMD, NHMS, TAU ) ENDIF !=================================== ! Ozone !=================================== ! Make note of inital state for determining tendency later BEFORE = STT(:,:,:,IDTOX ) ! Put ozone in v/v STT(:,:,:,IDTOX ) = STT(:,:,:,IDTOX) * TCVV( IDTOX ) / AD IF ( LLINOZ ) THEN CALL DO_LINOZ ! Linoz ELSE IF ( .not. LADJ ) THEN ! Must use Linoz or strat chem Ox fluxes for the adjoint CALL UPBDFLX_O3 ! Synoz ENDIF ! Now move unit conversion into LINOZ (hml, 11/06/11) ! Put ozone back to kg STT(:,:,:,IDTOX) = STT(:,:,:,IDTOX) * AD / TCVV( IDTOX ) ! Put tendency into diagnostic array [kg box-1] SCHEM_TEND(:,:,:,IDTOX) = SCHEM_TEND(:,:,:,IDTOX) & + ( STT(:,:,:,IDTOX) - BEFORE ) !====================================================================== ! Tagged Ox simulation !====================================================================== ELSE IF ( ITS_A_TAGOX_SIM() ) THEN ! Intial conditions STT0(:,:,:,:) = STT(:,:,:,:) CALL CONVERT_UNITS( 1, N_TRACERS, TCVV, AD, STT ) ! kg -> v/v IF ( LLINOZ ) THEN CALL DO_LINOZ ! Linoz ELSE IF ( .not. LADJ ) THEN ! must use Linoz or strat chem Ox fluxes for the adjoint CALL UPBDFLX_O3 ! Synoz ENDIF CALL CONVERT_UNITS( 2, N_TRACERS, TCVV, AD, STT ) ! v/v -> kg ! Add to tropopause level aggregator for later determining STE flux TpauseL_CNT = TpauseL_CNT + 1d0 !$OMP PARALLEL DO !$OMP+DEFAULT( SHARED ) !$OMP+PRIVATE( I, J ) DO I=1,IIPAR DO J=1,JJPAR TpauseL(I,J) = TpauseL(I,J) + GET_TPAUSE_LEVEL(I,J) ENDDO ENDDO !$OMP END PARALLEL DO ! Aggregate chemical tendency [kg box-1] DO N=1,NSCHEM NN = Strat_TrID_GC(N) SCHEM_TEND(:,:,:,NN) = SCHEM_TEND(:,:,:,NN) & + ( STT(:,:,:,NN) - STT0(:,:,:,NN) ) ENDDO ELSE ! The code will need to be modified for other tagged simulations ! (e.g., CO). Simulations like CH4, CO2 with standard tracer names ! should probably just work as is with the full chemistry code above, ! but would need to be tested. WRITE( 6, * ) 'Strat chemistry needs to be activated for ' // & 'your simulation type.' !WRITE( 6, * ) 'Please see GeosCore/strat_chem_mod.F90' // & WRITE( 6, * ) 'Please see new/strat_chem_mod.f' // & 'or disable in input.geos' CALL GEOS_CHEM_STOP ENDIF END SUBROUTINE DO_STRAT_CHEM !EOC !------------------------------------------------------------------------------ ! Harvard University Atmospheric Chemistry Modeling Group ! !------------------------------------------------------------------------------ !BOP ! ! !IROUTINE: GET_RATES ! ! !DESCRIPTION: Function GET\_RATES reads from disk the chemical production ! and loss rates for the species of interest !\\ !\\ ! !INTERFACE: ! SUBROUTINE GET_RATES( THISMONTH ) ! ! !USES: ! USE BPCH2_MOD, ONLY : GET_NAME_EXT, GET_RES_EXT USE BPCH2_MOD, ONLY : GET_TAU0, READ_BPCH2 USE DIRECTORY_MOD, ONLY : DATA_DIR USE TIME_MOD, ONLY : GET_MONTH USE LOGICAL_MOD, ONLY : LLINOZ USE TRACER_MOD, ONLY : N_TRACERS, TRACER_NAME, TRACER_COEFF USE TRANSFER_MOD, ONLY : TRANSFER_3D USE NETCDF_UTIL_MOD, ONLY : NCDF_GET_VAR USE NETCDF_UTIL_MOD, ONLY : NCDF_GET_VARID USE NETCDF_UTIL_MOD, ONLY : NCDF_OPEN_FOR_READ USE NETCDF_UTIL_MOD, ONLY : NCDF_CLOSE IMPLICIT NONE # include "CMN_SIZE" ! ! !INPUT PARAMETERS: ! ! Arguments INTEGER,INTENT(IN) :: THISMONTH ! ! !REVISION HISTORY: ! 1 Feb 2011 - L. Murray - Initial version !EOP !------------------------------------------------------------------------------ !BOC ! ! !LOCAL VARIABLES: ! CHARACTER(LEN=255) :: FILENAME CHARACTER(LEN=6) :: SPNAME( NTR_GMI ) REAL*4 :: ARRAY( IIPAR, JJPAR, LGLOB ) ! Full vertical res REAL*8 :: ARRAY2( IIPAR, JJPAR, LLPAR )! Actual vertical res REAL*8 :: XTAU INTEGER :: N, M, S, F, NN, fileID INTEGER :: prodID, lossID INTEGER :: ohID !================================================================= ! GET_RATES begins here !================================================================= ! Initialize arrays LOSS = 0d0 PROD = 0d0 WRITE(6, 11 ) & ' - Getting new strat prod/loss rates for month: ', & THISMONTH 11 FORMAT( a, I2.2 ) M = THISMONTH !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! Get stratospheric OH mixing ratio [v/v] !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% FILENAME = 'strat_chem_201206/gmi.clim.OH.' // & GET_NAME_EXT() // '.' // GET_RES_EXT() // '.nc' FILENAME = TRIM( DATA_DIR ) // TRIM( FILENAME ) WRITE(6,'(a)') ' => Reading from file: ' // trim(filename) call ncdf_open_for_read( fileID, TRIM( FILENAME ) ) ohID = ncdf_get_varid( fileID, 'species' ) call ncdf_get_var( fileID, ohID, array, & (/ 1, 1, 1, m /), ! Start & (/ iipar, jjpar, lglob, 1 /) ) ! Count call ncdf_close( fileID ) ! Cast from REAL*4 to REAL*8 and resize to 1:LLPAR call transfer_3D( array, array2 ) STRAT_OH(:,:,:) = ARRAY2 DO N=1,NSCHEM NN = Strat_TrID_GMI(N) !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! Open individual species file !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% FILENAME = 'strat_chem_201206/gmi.clim.' // & TRIM( GMI_TrName(NN) ) // '.' // & GET_NAME_EXT() // '.' // GET_RES_EXT() // '.nc' FILENAME = TRIM( DATA_DIR ) // TRIM( FILENAME ) WRITE(6,'(a)') ' => Reading from file: ' // trim(filename) call ncdf_open_for_read( fileID, TRIM( FILENAME ) ) !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! Read production rate [v/v/s] !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! Get the variable IDs for the species, prod and loss rates prodID = ncdf_get_varid( fileID, 'prod' ) call ncdf_get_var( fileID, prodID, array, & (/ 1, 1, 1, m /), ! Start & (/ iipar, jjpar, lglob, 1 /) ) ! Count ! Cast from REAL*4 to REAL*8 and resize to 1:LLPAR call transfer_3D( array, array2 ) PROD(:,:,:,N) = ARRAY2 ! Save rates from file to respective arrays (hml, 09/15/11) PROD_0(:,:,:,N) = PROD(:,:,:,N) !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! Read loss frequencies [s^-1] !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% lossID = ncdf_get_varid( fileID, 'loss' ) call ncdf_get_var( fileID, lossID, array, & (/ 1, 1, 1, m /), ! Start & (/ iipar, jjpar, lglob, 1 /) ) ! Count ! Cast from REAL*4 to REAL*8 and resize to 1:LLPAR call transfer_3D( array, array2 ) LOSS(:,:,:,N) = ARRAY2 ! Save rates from file to respective arrays (hml, 09/15/11) LOSS_0(:,:,:,N) = LOSS(:,:,:,N) !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! Close species file !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% call ncdf_close( fileID ) ENDDO END SUBROUTINE GET_RATES !EOC !------------------------------------------------------------------------------ ! Harvard University Atmospheric Chemistry Modeling Group ! !------------------------------------------------------------------------------ !BOP ! ! !IROUTINE: GET_RATES_INTERP ! ! !DESCRIPTION: Function GET\_RATES\_INTERP reads from disk the chemical ! production and loss rates for the species of interest to resolutions finer ! than 2 x 2.5 (e.g., nested simluations) via simple nearest-neighbor mapping. !\\ !\\ !\\ !\\ ! !INTERFACE: ! SUBROUTINE GET_RATES_INTERP( THISMONTH ) ! ! !USES: ! USE BPCH2_MOD, ONLY : GET_NAME_EXT, GET_RES_EXT USE DIRECTORY_MOD, ONLY : DATA_DIR USE GRID_MOD, ONLY : GET_YMID, GET_XMID USE LOGICAL_MOD, ONLY : LLINOZ USE TRACER_MOD, ONLY : N_TRACERS, TRACER_NAME, TRACER_COEFF USE TRANSFER_MOD, ONLY : TRANSFER_3D USE ADJ_ARRAYS_MOD, ONLY : NSTPL USE NETCDF_UTIL_MOD, ONLY : NCDF_GET_VAR USE NETCDF_UTIL_MOD, ONLY : NCDF_GET_VARID USE NETCDF_UTIL_MOD, ONLY : NCDF_OPEN_FOR_READ USE NETCDF_UTIL_MOD, ONLY : NCDF_CLOSE # include "define.h" # include "CMN_SIZE" ! ! !INPUT PARAMETERS: ! ! Arguments INTEGER,INTENT(IN) :: THISMONTH ! ! !REVISION HISTORY: ! 1 Feb 2011 - L. Murray - Initial version !EOP !------------------------------------------------------------------------------ !BOC ! ! !LOCAL VARIABLES: ! CHARACTER(LEN=255) :: FILENAME CHARACTER(LEN=6) :: SPNAME( NTR_GMI ) REAL*4 :: ARRAY( IIPAR, JJPAR, LGLOB ) REAL*8 :: ARRAY2( IIPAR, JJPAR, LLPAR ) INTEGER :: N, M, S, F INTEGER :: NN INTEGER :: ohID INTEGER :: prodID, lossID INTEGER :: lat_varID, lon_varID REAL*4 :: XMID_COARSE(144), YMID_COARSE(91) INTEGER :: I_f2c(IIPAR), J_f2c(JJPAR) ! f2c = fine to coar map'ng INTEGER :: I, J, fileID INTEGER :: II(1), JJ(1) REAL*4 :: COLUMN( LGLOB ) !================================================================= ! GET_RATES_INTERP begins here !================================================================= ! In the original schem code, the following species were destroyed ! by photolysis in the stratosphere: ! PAN, H2O2, ACET, MEK, ALD2, RCHO, MVK, MACR, R4N2, CH2O, ! N2O5, HNO4, MP ! And by reaction with OH for: ! ALK4, ISOP, H2O2, ACET, MEK, ALD2, RCHO, MVK, MACR, PMN, R4N2, ! PRPE, C3H8, CH2O, C2H6, HNO4, MP ! The updated code includes at least all of these, and several more. ! Initialize arrays LOSS = 0d0 PROD = 0d0 ! first read in the OH file so that we can get the lat and long ! values to populate XMID_COARSE and YMID_COARSE ! Path to input data, use 2 x 2.5 file #if defined( GEOS_FP ) FILENAME = 'strat_chem_201206/gmi.clim.OH.geos5' // & '.2x25.nc' #else FILENAME = 'strat_chem_201206/gmi.clim.OH.' // GET_NAME_EXT() // & '.2x25.nc' #endif !FILENAME = TRIM( DATA_DIR_1x1 ) // TRIM( FILENAME ) !FILENAME = TRIM( DATA_DIR ) // TRIM('../GEOS_2x25/') // FILENAME = TRIM( DATA_DIR ) // TRIM('../GEOS_2x2.5/') // & TRIM( FILENAME ) WRITE(6, 11 ) & ' - Getting new strat prod/loss rates for month: ', & THISMONTH 11 FORMAT( a, I2.2 ) ! Open the netCDF file containing the rates WRITE(6,'(a)') & ' => Interpolate to resolution from file: ' & // trim(filename) call ncdf_open_for_read( fileID, TRIM( filename ) ) ! Get the lat and lon centers of the 2x2.5 GMI climatology ! WARNING MAKE 2x25 after testing !call NcRd( XMID_COARSE, fileID, 'longitude', (/1/), (/144/) ) !call NcRd( YMID_COARSE, fileID, 'latitude', (/1/), (/91/) ) lat_varID = ncdf_get_varid( fileID, 'latitude' ) lon_varID = ncdf_get_varid( fileID, 'longitude' ) !call NcRd( XMID_COARSE, fileID, 'longitude', (/1/), (/144/) ) !call NcRd( YMID_COARSE, fileID, 'latitude', (/1/), (/91/) ) call ncdf_get_var( fileID, lon_varID, XMID_COARSE, (/1/), (/144/)) call ncdf_get_var( fileID, lat_varID, YMID_COARSE, (/1/), (/91/) ) ! For each fine grid index, determine the closest coarse (2x2.5) index ! Note: This doesn't do anything special for the date line, and may ! therefore not pick the exact closest if it is on the other side. ! Note: CMN_SIZE_MOD claims in its comments that IIPAR < IGLOB, but ! in actuality, IIPAR = IGLOB and JJPAR = JGLOB, the dimensions of the nested ! region. DO I=1,IGLOB II = MINLOC( ABS( GET_XMID(I) - XMID_COARSE ) ) I_f2c(I) = II(1) !print*,'I:',I,'->',II(1) ENDDO DO J=1,JGLOB JJ = MINLOC( ABS( GET_YMID(J) - YMID_COARSE ) ) J_f2c(J) = JJ(1) !print*,'J:',J,'->',JJ(1) ENDDO ! call ncdf_close( fileID ) M = THISMONTH !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !Get Stratospheric OH concentrations [v/v] !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ohID = ncdf_get_varid( fileID, 'species' ) DO I=1,IGLOB DO J=1,JGLOB call ncdf_get_var( fileID, ohID, column, & (/ I_f2c(I), J_f2c(J), 1, m /), ! Start & (/ 1, 1, lglob, 1 /) ) ! Count array( I, J, : ) = column ENDDO ENDDO call ncdf_close( fileID ) call transfer_3D( array, array2 ) STRAT_OH(:,:,:) = ARRAY2 DO N=1,NSCHEM NN = Strat_TrID_GMI(N) !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! Open individual species file !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! Path to input data, use 2 x 2.5 file ! (lzh,02/01/2015) add geosfp - use geos-5 #if defined( GEOS_FP ) FILENAME = 'strat_chem_201206/gmi.clim.' // & TRIM( GMI_TrName(NN) ) // '.geos5' // & '.2x25.nc' #else FILENAME = 'strat_chem_201206/gmi.clim.' // & TRIM( GMI_TrName(NN) ) // '.' // GET_NAME_EXT() // & '.2x25.nc' #endif !FILENAME = TRIM( DATA_DIR ) // TRIM('../GEOS_2x25/') // FILENAME = TRIM( DATA_DIR ) // TRIM('../GEOS_2x2.5/') // & TRIM( FILENAME ) WRITE(6,'(a)') & ' => Interpolate to resolution from file: ' // & trim(filename) call ncdf_open_for_read( fileID, TRIM( FILENAME ) ) !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! Read production rate [v/v/s] !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% array = 0.0 prodID = ncdf_get_varid( fileID, 'prod' ) DO I=1,IGLOB DO J=1,JGLOB call ncdf_get_var( fileID, prodID, column, & (/ I_f2c(I), J_f2c(J), 1, m /), ! Start & (/ 1, 1, lglob, 1 /) ) ! Count array( I, J, : ) = column ENDDO ENDDO ! Cast from REAL*4 to REAL*8 and resize to 1:LLPAR if necessary call transfer_3D( array, array2 ) PROD(:,:,:,N) = ARRAY2 ! Save rates from file to respective arrays (hml, 09/15/11) PROD_0(:,:,:,N) = PROD(:,:,:,N) !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! Read loss frequencies [s^-1] !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% array = 0.0 lossID = ncdf_get_varid( fileID, 'loss' ) DO I=1,IGLOB DO J=1,JGLOB call ncdf_get_var( fileID, lossID, column, & (/ I_f2c(I), J_f2c(J), 1, m /), ! Start & (/ 1, 1, lglob, 1 /) ) ! Count array( I, J, : ) = column ENDDO ENDDO ! Cast from REAL*4 to REAL*8 and resize to 1:LLPAR if necessary call transfer_3D( array, array2 ) LOSS(:,:,:,N) = ARRAY2 ! Save rates from file to respective arrays (hml, 09/15/11) LOSS_0(:,:,:,N) = LOSS(:,:,:,N) call ncdf_close( fileID ) ENDDO !call ncdf_close( ncID_strat_rates ) END SUBROUTINE GET_RATES_INTERP !EOC !------------------------------------------------------------------------------ ! Harvard University Atmospheric Chemistry Modeling Group ! !------------------------------------------------------------------------------ !BOP ! ! !IROUTINE: calc_ste ! ! !DESCRIPTION: Subroutine CALC\_STE estimates what the stratosphere-to- ! troposphere exchange flux must have been since the last time ! it was reset !\\ !\\ ! !INTERFACE: ! SUBROUTINE CALC_STE ! ! !USES: ! USE TRACER_MOD, ONLY : STT, TRACER_MW_KG, N_TRACERS, TRACER_NAME USE TIME_MOD, ONLY : GET_TAU, GET_NYMD, GET_NHMS, EXPAND_DATE IMPLICIT NONE #include "define.h" #include "CMN_SIZE" ! ! !REVISION HISTORY: ! 28 Apr 2012 - L. Murray - Initial version !EOP !------------------------------------------------------------------------------ !BOC REAL*8 :: M1(IIPAR,JJPAR,LLPAR), M2(IIPAR,JJPAR,LLPAR) REAL*8 :: dStrat, STE, Tend, tauEnd, dt INTEGER :: N, I, J, L, NN, LTP(IIPAR,JJPAR) CHARACTER(LEN=256) :: dateStart, dateEnd !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! By simple mass balance, dStrat/dt = P - L - STE, ! where STE is the net stratosphere-to-troposphere mass exchange. ! ! Therefore, we estimate STE as ! STE = (P-L) - dStrat/dt ! ! As the tropopause is dynamic, we use the mean tropopause level during ! the period for determining initial and end stratospheric masses. ! (ltm, 04/28/2012) !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #if defined( NESTED_NA ) || defined( NESTED_CH ) || defined( NESTED_EU ) ! This method only works for a global domain. ! It could be modified for nested domains if the total mass flux across the ! boundaries during the period is taken into account. RETURN #endif ! Determine mean tropopause level for the period !$OMP PARALLEL DO !$OMP+DEFAULT( SHARED ) !$OMP+PRIVATE( I, J ) DO I = 1,IIPAR DO J = 1,JJPAR LTP(I,J) = NINT( TPauseL(I,J) / TPauseL_Cnt ) ENDDO ENDDO !$OMP END PARALLEL DO ! Period over which STE is being determined [a] tauEnd = GET_TAU() ! [h] dt = ( tauEnd - tauInit ) / 24d0 / 365.25d0 dateStart = 'YYYY-MM-DD hh:mm' CALL EXPAND_DATE(dateStart,NymdInit,NhmsInit) dateEnd = 'YYYY-MM-DD hh:mm' CALL EXPAND_DATE(dateEnd,GET_NYMD(),GET_NHMS()) ! Print to output WRITE( 6, * ) '' WRITE( 6, '(a)' ) REPEAT( '=', 79 ) WRITE( 6, '(a)' ) ' Strat-Trop Exchange' WRITE( 6, '(a)' ) REPEAT( '-', 79 ) WRITE( 6, '(a)' ) & ' Global stratosphere-to-troposphere fluxes estimated over' WRITE( 6, 100 ) TRIM(dateStart), TRIM(dateEnd) WRITE( 6, * ) '' WRITE( 6, 110 ) 'Species','[moles a-1]','* [g/mol]','= [Tg a-1]' 100 FORMAT( 2x,a16,' to ',a16 ) 110 FORMAT( 2x,a8,':',4x,a11 ,4x,a9 ,4x, a11 ) ! Loop through each species DO N=1,N_TRACERS ! Populate before (M1) and after (M2) state for the species [kg] M1 = MInit(:,:,:,N) M2 = STT(:,:,:,N) ! Zero out tropopshere and determine total change in the stratospheric ! burden of species N (dStrat) [kg] !$OMP PARALLEL DO !$OMP+DEFAULT( SHARED ) !$OMP+PRIVATE( I, J ) DO I=1,IIPAR DO J=1,JJPAR M2(I,J,1:LTP(I,J)) = 0d0 M1(I,J,1:LTP(I,J)) = 0d0 ENDDO ENDDO !$OMP END PARALLEL DO dStrat = SUM(M2)-SUM(M1) ! The total chemical tendency (P-L) over the period for species N [kg] Tend = SUM(Schem_tend(:,:,:,N)) ! Calculate flux as STE = (P-L) - dStrat/dt STE = (Tend-dStrat)/dt ! [kg a-1] ! Print to standard output WRITE(6,120) TRIM(TRACER_NAME(N)), & STE/TRACER_MW_KG(N), ! mol/a-1 & TRACER_MW_KG(N)*1d3, ! g/mol & STE*1d-9 ! Tg a-1 ENDDO 120 FORMAT( 2x,a8,':',4x,e11.3,4x,f9.1,4x,f11.4 ) WRITE( 6, * ) '' WRITE( 6, '(a)' ) REPEAT( '=', 79 ) WRITE( 6, * ) '' ! Reset variables for next STE period NymdInit = GET_NYMD() NhmsInit = GET_NHMS() TauInit = GET_TAU() TPauseL_Cnt = 0d0 TPauseL(:,:) = 0d0 SChem_tend(:,:,:,:) = 0d0 MInit(:,:,:,:) = STT(:,:,:,:) END SUBROUTINE CALC_STE !EOC !------------------------------------------------------------------------------ ! Harvard University Atmospheric Chemistry Modeling Group ! !------------------------------------------------------------------------------ !BOP ! ! !IROUTINE: init_strat_chem ! ! !DESCRIPTION: Subroutine INIT\_STRAT\_CHEM allocates all module arrays. ! It also opens the necessary rate files. !\\ !\\ ! !INTERFACE: ! SUBROUTINE INIT_STRAT_CHEM ! ! !USES: ! USE ERROR_MOD, ONLY : ALLOC_ERR USE LOGICAL_MOD, ONLY : LLINOZ USE TRACER_MOD, ONLY : ITS_A_FULLCHEM_SIM, ITS_A_TAGOX_SIM USE TRACER_MOD, ONLY : N_TRACERS, TRACER_NAME, TRACER_COEFF USE TRACER_MOD, ONLY : TRACER_COEFF, STT USE BPCH2_MOD, ONLY : GET_NAME_EXT, GET_RES_EXT USE DIRECTORY_MOD, ONLY : DATA_DIR USE TIME_MOD, ONLY : EXPAND_DATE USE TIME_MOD, ONLY : GET_TAU, GET_NYMD, GET_NHMS USE TIME_MOD, ONLY : GET_TS_CHEM USE LOGICAL_ADJ_MOD, ONLY : LADJ USE UPBDFLX_MOD, ONLY : UPBDFLX_O3, INIT_UPBDFLX # include "CMN_SIZE" ! ! !REVISION HISTORY: ! 1 Feb 2011 - L. Murray - Initial version !EOP !------------------------------------------------------------------------------ !BOC ! ! !LOCAL VARIABLES: ! CHARACTER(LEN=16) :: sname INTEGER :: AS, N, NN CHARACTER(LEN=255) :: FILENAME, FILENAMEOUT CHARACTER(LEN=6) :: SPNAME( NTR_GMI ) INTEGER :: spname_varID !================================================================= ! INIT_STRAT_CHEM begins here! !================================================================= ! Initialize counters, initial times, mapping arrays TpauseL_Cnt = 0. NSCHEM = 0 TauInit = GET_TAU() NymdInit = GET_NYMD() NhmsInit = GET_NHMS() strat_trID_GC(:) = 0 strat_trID_GMI(:) = 0 ! Initialize timestep for chemistry [s] dTchem = GET_TS_CHEM() * 60d0 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! Determine the mapping for the GMI to the GC variables based on ! tracer name, which only needs to be done once per model run. !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! List of available tracers with archived monthly climatological ! production rates, loss frequencies, and mixing ratios from the ! GMI Combo model (tracer names here are as used in GMI). GMI_TrName = (/ & 'A3O2', 'ACET', 'ACTA', 'ALD2', 'ALK4', 'ATO2', & 'B3O2', 'Br', 'BrCl', 'BrO', 'BrONO2', 'C2H6', & 'C3H8', 'CCl4', 'CF2Br2', 'CF2Cl2', 'CF2ClBr', 'CF3Br', & 'CFC113', 'CFC114', 'CFC115', 'CFCl3', 'CH2O', 'CH3Br', & 'CH3CCl3', 'CH3Cl', 'CH4', 'CO', 'Cl', 'Cl2', & 'Cl2O2', 'ClO', 'ClONO2', 'EOH', 'ETO2', 'ETP', & 'GCO3', 'GLYC', 'GLYX', 'GP', 'GPAN', 'H', & 'H2', 'H2402', 'H2O', 'H2O2', 'HAC', 'HBr', & 'HCFC141b', 'HCFC142b', 'HCFC22', 'HCOOH', 'HCl', 'HNO2', & 'HNO3', 'HNO4', 'HO2', 'HOBr', 'HOCl', 'IALD', & 'IAO2', 'IAP', 'INO2', 'INPN', 'ISN1', 'ISNP', & 'ISOP', 'KO2', 'MACR', 'MAN2', 'MAO3', 'MAOP', & 'MAP', 'MCO3', 'MEK', 'MGLY', 'MO2', 'MOH', & 'MP', 'MRO2', 'MRP', 'MVK', 'MVN2', 'N', & 'N2O', 'N2O5', 'NO', 'NO2', 'NO3', 'NOx', & 'O', 'O1D', 'O3', 'OClO', 'OH', 'Ox', & 'PAN', 'PMN', 'PO2', 'PP', 'PPN', 'PRN1', & 'PRPE', 'PRPN', 'R4N1', 'R4N2', 'R4O2', 'R4P', & 'RA3P', 'RB3P', 'RCHO', 'RCO3', 'RCOOH', 'RIO1', & 'RIO2', 'RIP', 'ROH', 'RP', 'VRO2', 'VRP' /) !===========================! ! Full chemistry simulation ! !===========================! IF ( ITS_A_FULLCHEM_SIM() ) THEN DO NN = 1, NTR_GMI sname = TRIM(GMI_TrName(NN)) DO N = 1, N_TRACERS IF ( TRIM(TRACER_NAME(N)) .eq. TRIM(sname) ) THEN IF ( LLINOZ .and. & TRIM(TRACER_NAME(N)) .eq. 'Ox' ) THEN WRITE(6,*), TRIM(TRACER_NAME(N)) // ' (via Linoz)' ! For adjoint !ELSE IF ( TRIM(TRACER_NAME(N)) .eq. 'Ox' ) THEN ! WRITE(6,*) TRIM(TRACER_NAME(N)) // ' (via Synoz)' ELSEIF ( TRIM(TRACER_NAME(N)) .eq. 'Ox' ) THEN WRITE(6,*),TRIM(TRACER_NAME(N)) // ' (Ox via GMI)' ELSE WRITE(6,*), TRIM(TRACER_NAME(N)) // & ' (via GMI rates)' ENDIF NSCHEM = NSCHEM + 1 Strat_TrID_GC(NSCHEM) = N ! Maps 1:NSCHEM to STT index Strat_TrID_GMI(NSCHEM) = NN ! Maps 1:NSCHEM to GMI_TrName index ENDIF ENDDO ENDDO ! Allocate array to hold monthly mean OH mixing ratio ALLOCATE( STRAT_OH( IIPAR, JJPAR, LLPAR ), STAT=AS ) IF ( AS /=0 ) CALL ALLOC_ERR( 'STRAT_OH' ) STRAT_OH = 0d0 !===========! ! Tagged Ox ! !===========! ELSE IF ( ITS_A_TAGOX_SIM() ) THEN IF ( LLINOZ ) THEN WRITE(6,*) 'Linoz ozone performed on: ' ELSE IF ( .not. LADJ ) THEN ! must use Linoz or strat chem Ox fluxes for the adjoint CALL UPBDFLX_O3 ! Synoz WRITE(6,*) 'Synoz ozone performed on: ' ENDIF DO N = 1, N_TRACERS IF ( TRIM(TRACER_NAME(N)) .eq. 'Ox' .or. & TRIM(TRACER_NAME(N)) .eq. 'OxStrt' ) THEN NSCHEM = NSCHEM + 1 Strat_TrID_GC(NSCHEM) = N WRITE(6,*) TRIM(TRACER_NAME(N)) ENDIF ENDDO ENDIF !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%! ! Allocate and initialize prod & loss arrays ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%! ! Allocate PROD -- array for clim. production rates [v/v/s] ALLOCATE( PROD( IIPAR, JJPAR, LLPAR, NSCHEM ), STAT=AS ) IF ( AS /= 0 ) CALL ALLOC_ERR( 'PROD' ) PROD = 0d0 ! Allocate LOSS -- array for clim. loss freq [s-1] ALLOCATE( LOSS( IIPAR, JJPAR, LLPAR, NSCHEM ), STAT=AS ) IF ( AS /= 0 ) CALL ALLOC_ERR( 'LOSS' ) LOSS = 0d0 ! For adjoint !ALLOCATE( PROD_0( IIPAR, JJPAR, LLPAR, N_TRACERS ), STAT=AS ) ALLOCATE( PROD_0( IIPAR, JJPAR, LLPAR, NSCHEM ), STAT=AS ) IF ( AS /= 0 ) CALL ALLOC_ERR( 'PROD_0' ) PROD_0 = 0d0 !ALLOCATE( LOSS_0( IIPAR, JJPAR, LLPAR, N_TRACERS ), STAT=AS ) ALLOCATE( LOSS_0( IIPAR, JJPAR, LLPAR, NSCHEM ), STAT=AS ) IF ( AS /= 0 ) CALL ALLOC_ERR( 'LOSS_0' ) LOSS_0 = 0d0 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%! ! Allocate and initialize arrays for STE calculation ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%! ! Array to hold initial state of atmosphere at the beginning ! of the period over which to estimate STE. Populate with ! initial atm. conditions from restart file [kg]. ALLOCATE( MInit( IIPAR, JJPAR, LLPAR, N_TRACERS ), STAT=AS ) IF ( AS /= 0 ) CALL ALLOC_ERR( 'MInit' ) MInit = STT ! Array to determine the mean tropopause level over the period ! for which STE is being estimated. ALLOCATE( TPAUSEL( IIPAR, JJPAR ), STAT=AS ) IF ( AS /= 0 ) CALL ALLOC_ERR( 'TPAUSEL' ) TPAUSEL = 0d0 ! Array to save chemical state before each chemistry time step [kg] ALLOCATE( BEFORE( IIPAR, JJPAR, LLPAR ), STAT=AS ) IF ( AS /= 0 ) CALL ALLOC_ERR( 'BEFORE' ) BEFORE = 0d0 ! Array to aggregate the stratospheric chemical tendency [kg period-1] ALLOCATE( SCHEM_TEND(IIPAR,JJPAR,LLPAR,N_TRACERS), STAT=AS ) IF ( AS /= 0 ) CALL ALLOC_ERR( 'SCHEM_TEND' ) SCHEM_TEND = 0d0 END SUBROUTINE INIT_STRAT_CHEM !EOC !------------------------------------------------------------------------------ ! Harvard University Atmospheric Chemistry Modeling Group ! !------------------------------------------------------------------------------ !BOP ! ! !IROUTINE: cleanup_strat_chem ! ! !DESCRIPTION: Subroutine CLEANUP\_STRAT\_CHEM deallocates all module ! arrays. !\\ !\\ ! !INTERFACE: ! SUBROUTINE CLEANUP_STRAT_CHEM ! ! !USES: ! USE NETCDF_UTIL_MOD ! !REVISION HISTORY: ! 1 Feb 2011 - L. Murray - Initial version !EOP !------------------------------------------------------------------------------ !BOC IF ( ALLOCATED( PROD ) ) DEALLOCATE( PROD ) IF ( ALLOCATED( LOSS ) ) DEALLOCATE( LOSS ) IF ( ALLOCATED( PROD_0 ) ) DEALLOCATE( PROD_0 ) IF ( ALLOCATED( LOSS_0 ) ) DEALLOCATE( LOSS_0 ) IF ( ALLOCATED( STRAT_OH ) ) DEALLOCATE( STRAT_OH ) IF ( ALLOCATED( MInit ) ) DEALLOCATE( MInit ) IF ( ALLOCATED( TPAUSEL ) ) DEALLOCATE( TPAUSEL ) IF ( ALLOCATED( BEFORE ) ) DEALLOCATE( BEFORE ) IF ( ALLOCATED( SCHEM_TEND ) ) DEALLOCATE( SCHEM_TEND ) END SUBROUTINE CLEANUP_STRAT_CHEM !EOC END MODULE STRAT_CHEM_MOD