! $Id: Kr85_mod.f,v 1.1 2009/06/09 21:51:52 daven Exp $ MODULE Kr85_MOD ! !****************************************************************************** ! Module Kr85_MOD contains routines and variables for the Kr85 radionuclide ! simulation. (jsw, bmy, 8/21/03, 11/6/08) ! ! Module Variables: ! ============================================================================ ! (1 ) N_SOURCES (INTEGER) : Maximum number of Kr85 point sources ! (2 ) N_YEARS (INTEGER) : Maximum number of years for Kr85 emissions ! (3 ) SMALLNUM (REAL*8 ) : A small number, used to prevent underflow ! ! Module Routines: ! ============================================================================ ! (1 ) GET_SOURCE_IJ : Returns (I,J) location of each Kr85 point source ! (2 ) GET_EMITTED_Kr85 : Returns Kr85 emission from a point src in [kg] ! (3 ) EMISSKr85 : Adds Kr85 emissions into the tracer array ! (4 ) CHEMKr85 : Performs radioactive (1st-order) loss for Kr85 ! ! GEOS-CHEM modules referenced by biomass_mod.f ! ============================================================================ ! (1 ) diag_mod.f : Module containing GEOS-CHEM diagnostic arrays ! (2 ) error_mod.f : Module containing I/O error and NaN check routines ! (3 ) time_mod.f : Module containing routines for computing time & date ! (4 ) tracer_mod.f : Module containing GEOS-CHEM tracer array STT etc. ! ! References: ! ============================================================================ ! (1 ) Jacob, D.J., M.J. Prather, S.C. Wofsy, M.B. McElroy, "Atmospheric ! distribution of 85Kr simulated with a general circulation model", ! JGR, 92(D6), pp. 6614-6626, June 20, 1987. ! ! NOTES: ! (1 ) Now references "tracer_mod.f" (bmy, 7/20/04) ! (2 ) Modifications for GEOS-5 nested grids (yxw, dan, bmy, 11/6/08) !****************************************************************************** ! IMPLICIT NONE !================================================================= ! MODULE PRIVATE DECLARATIONS -- keep certain internal variables ! and routines from being seen outside "Kr85_mod.f" !================================================================= ! PRIVATE module routines PRIVATE :: GET_SOURCE_IJ PRIVATE :: GET_EMITTED_KR85 ! PRIVATE module variables PRIVATE :: N_SOURCES, N_YEARS, SMALLNUM !================================================================= ! MODULE VARIABLES !================================================================= INTEGER, PARAMETER :: N_SOURCES = 8 INTEGER, PARAMETER :: N_YEARS = 6 REAL*8, PARAMETER :: SMALLNUM = 1d-20 !================================================================= ! MODULE ROUTINES -- follow below the "CONTAINS" statement !================================================================= CONTAINS !------------------------------------------------------------------------------ SUBROUTINE GET_SOURCE_IJ( N_SOURCE, I, J ) ! !****************************************************************************** ! Subroutine GET_SOURCE_IJ returns the (I,J) grid box indices for each ! Kr85 point source. For now we have hardwired this, since there are only ! a few stations. Worry about making this more general at some future time. ! (bmy, 8/21/03, 11/6/08) ! ! Arguments as Input: ! ============================================================================ ! (1 ) N_SOURCE (INTEGER) : Number of Kr85 point source (1-8) ! ! Arguments as Output ! ============================================================================ ! (1-2) I, J (INTEGER) : Lon & lat indices for the N_SOURCEth Kr85 source ! ! NOTES: ! (1 ) Updated for 0.5 x 0.666 nested grids (yxw, dan, bmy, 11/6/08) !****************************************************************************** ! ! References to F90 modules USE ERROR_MOD, ONLY : ERROR_STOP # include "define.h" ! Switches ! Arguments INTEGER, INTENT(IN) :: N_SOURCE INTEGER, INTENT(OUT) :: I, J ! Local variables INTEGER :: GRID ! Station (I,J) arrays 4x5 2x25 1x1 INTEGER :: I1(3) = (/ 21, 41, 101 /) INTEGER :: J1(3) = (/ 32, 63, 105 /) !--- INTEGER :: I2(3) = (/ 14, 27, 66 /) INTEGER :: J2(3) = (/ 35, 69, 137 /) !--- INTEGER :: I3(3) = (/ 36, 71, 176 /) INTEGER :: J3(3) = (/ 37, 73, 145 /) !--- INTEGER :: I4(3) = (/ 38, 75, 186 /) INTEGER :: J4(3) = (/ 35, 69, 137 /) !--- INTEGER :: I5(3) = (/ 37, 73, 181 /) INTEGER :: J5(3) = (/ 36, 71, 141 /) !--- INTEGER :: I6(3) = (/ 39, 77, 191 /) INTEGER :: J6(3) = (/ 36, 71, 141 /) !--- INTEGER :: I7(3) = (/ 65, 129, 321 /) INTEGER :: J7(3) = (/ 32, 63, 125 /) !--- INTEGER :: I8(3) = (/ 49, 97, 241 /) INTEGER :: J8(3) = (/ 37, 73, 145 /) !================================================================= ! GET_SOURCE_IJ begins here! !================================================================= ! Select flag for grid type #if defined( GRID4x5 ) GRID = 1 #elif defined( GRID2x25 ) GRID = 2 #elif defined( GRID1x1 ) GRID = 3 #elif defined( GRID05x0666 ) GRID = 3 !(dan ) #endif ! Select proper (I,J) for each station SELECT CASE( N_SOURCE ) CASE( 1 ) I = I1(GRID) J = J1(GRID) CASE( 2 ) I = I2(GRID) J = J2(GRID) CASE( 3 ) I = I3(GRID) J = J3(GRID) CASE( 4 ) I = I4(GRID) J = J4(GRID) CASE( 5 ) I = I5(GRID) J = J5(GRID) CASE( 6 ) I = I6(GRID) J = J6(GRID) CASE( 7 ) I = I7(GRID) J = J7(GRID) CASE( 8 ) I = I8(GRID) J = J8(GRID) CASE DEFAULT CALL ERROR_STOP( 'N_SOURCE must be between 1-8!', & 'GET_SOURCE_IJ (Kr85_mod.f)' ) END SELECT ! Return to calling program END SUBROUTINE GET_SOURCE_IJ !------------------------------------------------------------------------------ FUNCTION GET_EMITTED_Kr85( N_SOURCE, YEARCOUNT ) RESULT( Kr85 ) ! !****************************************************************************** ! Subroutine GET_EMITTED_Kr85 returns the amount of Kr85 emitted from a ! given point source ! ! Arguments as Input: ! ============================================================================ ! (1 ) N_SOURCE (INTEGER) : Kr85 point source index (1-N_SOURCES) ! (2 ) YEARCOUNT (INTEGER) : Year of Kr85 emissions to use (1-N_YEARS) ! ! NOTES: !****************************************************************************** ! ! References to F90 modules USE ERROR_MOD, ONLY : ERROR_STOP USE TIME_MOD, ONLY : GET_TS_EMIS ! Arguments INTEGER, INTENT(IN) :: N_SOURCE, YEARCOUNT ! Local variables REAL*8 :: Kr85, DTSRCE REAL*8, PARAMETER :: SEC_PER_YR = 86400d0 * 365.25d0 ! Kr85 point src emissions ! Units: MCi/year Year1 Year2 Year3 Year4 Year5 Year6 REAL*8 :: S1(6) = (/ 0.48d0,0.46d0,0.47d0,0.54d0,0.56d0,0.69d0 /) REAL*8 :: S2(6) = (/ 0.10d0,0.00d0,0.09d0,0.06d0,0.01d0,0.00d0 /) REAL*8 :: S3(6) = (/ 0.70d0,0.94d0,0.84d0,1.40d0,1.19d0,1.13d0 /) REAL*8 :: S4(6) = (/ 0.31d0,0.28d0,0.54d0,0.31d0,0.31d0,0.31d0 /) REAL*8 :: S5(6) = (/ 0.79d0,0.64d0,0.83d0,0.91d0,1.27d0,1.95d0 /) REAL*8 :: S6(6) = (/ 0.00d0,0.05d0,0.03d0,0.00d0,0.00d0,0.08d0 /) REAL*8 :: S7(6) = (/ 0.06d0,0.00d0,0.28d0,0.11d0,0.19d0,0.09d0 /) REAL*8 :: S8(6) = (/ 3.56d0,3.77d0,3.19d0,3.07d0,3.00d0,2.40d0 /) !================================================================= ! GET_EMITTED_Kr85 begins here! !================================================================= ! Error check year IF ( YEARCOUNT < 1 .or. YEARCOUNT > 6 ) THEN CALL ERROR_STOP( 'YEARCOUNT must be between 1-6!', & 'GET_EMITTED_KR85 (Kr85_mod.f)' ) ENDIF ! Return Kr85 for the given point source & year SELECT CASE( N_SOURCE ) CASE( 1 ) Kr85 = S1(YEARCOUNT) CASE( 2 ) Kr85 = S2(YEARCOUNT) CASE( 3 ) Kr85 = S3(YEARCOUNT) CASE( 4 ) Kr85 = S4(YEARCOUNT) CASE( 5 ) Kr85 = S5(YEARCOUNT) CASE( 6 ) Kr85 = S6(YEARCOUNT) CASE( 7 ) Kr85 = S7(YEARCOUNT) CASE( 8 ) Kr85 = S8(YEARCOUNT) CASE DEFAULT CALL ERROR_STOP( 'N_SOURCE must be between 1-8!', & 'GET_SOURCE_IJ (Kr85_mod.f)' ) END SELECT ! Emission timestep [s] DTSRCE = GET_TS_EMIS() * 60d0 ! Convert from [MCi/yr] to [kg/emission timestep] ! 1 kg of Kr85 is equivalent to 2.55 MCi (cf Jacob et al 1987) Kr85 = Kr85 * 2.55d0 * ( DTSRCE / SEC_PER_YR ) ! Return to calling program END FUNCTION GET_EMITTED_Kr85 !------------------------------------------------------------------------------ SUBROUTINE EMISSKr85 ! !****************************************************************************** ! Subroutine EMISSKr85 places Kr85 emissions from point sources (e.g. nuclear ! reprocessing plants) into the tracer array. (jsw, bmy, 8/21/03, 7/20/04) ! ! NOTES: ! (1 ) Now references STT from "tracer_mod.f" (bmy, 7/20/04) !****************************************************************************** ! ! References to F90 modules !-------------------------------------------------------------- ! Prior to 12/7/04: ! Need to reassign the diagnostic number !USE DIAG_MOD, ONLY : AD03 !-------------------------------------------------------------- USE TIME_MOD, ONLY : GET_TS_EMIS, GET_YEAR USE TRACER_MOD, ONLY : STT # include "CMN_SIZE" ! Size parameters # include "CMN_DIAG" ! Diagnostics # include "CMN_O3" ! FSCALYR ! Local Variables LOGICAL, SAVE :: FIRSTEMISS = .TRUE. INTEGER, SAVE :: YEARCOUNT, LASTYEAR INTEGER :: I, J, N REAL*8 :: Kr85_KG, DTSRCE !================================================================= ! EMISSKr85 begins here! !================================================================= ! Emission timestep [s] DTSRCE = GET_TS_EMIS() * 60d0 ! First-time initialization IF ( FIRSTEMISS ) THEN YEARCOUNT = FSCALYR LASTYEAR = GET_YEAR() FIRSTEMISS = .FALSE. WRITE( 6, 100 ) YEARCOUNT 100 FORMAT( ' - EMISSKr85: Using Kr85 emissions from year ',i3) ENDIF ! If it's a new year, increment YEARCOUNT IF ( GET_YEAR() /= LASTYEAR ) THEN YEARCOUNT = YEARCOUNT + 1 LASTYEAR = GET_YEAR() WRITE( 6, 100 ) YEARCOUNT ENDIF !================================================================= ! Add Kr85 emissions [kg] to the STT tracer array ! NOTE: Assumes a global (not a window!) simulation !================================================================= DO N = 1, N_SOURCES ! Get (I,J) for each Kr85 point source CALL GET_SOURCE_IJ( N, I, J ) ! Get emitted Kr85 from each point source [kg] Kr85_KG = GET_EMITTED_Kr85( N, YEARCOUNT ) ! Add Kr85 into STT array STT(I,J,1,1) = STT(I,J,1,1) + Kr85_KG !-------------------------------------------------------------- ! Prior to 12/7/04: ! Need to reassign the diagnostic number (bmy, 12/7/04) !! Archive emitted Kr85 for ND04 diagnostic [kg] !IF ( ND03 > 0 ) THEN ! AD03(I,J,1,1) = AD03(I,J,1,1) + Kr85_KG !ENDIF !-------------------------------------------------------------- ENDDO ! Return to calling program END SUBROUTINE EMISSKr85 !------------------------------------------------------------------------------ SUBROUTINE CHEMKr85 ! !****************************************************************************** ! Subroutine CHEMKr85 applies first-order loss to the Kr85 tracer. ! (jsw, bmy, 8/21/03, 7/20/04) ! ! NOTES: ! (1 ) Now references STT from "tracer_mod.f" (bmy, 7/20/04) !****************************************************************************** ! ! References to F90 modules !------------------------------------------ ! Prior to 12/7/04: ! Need to reassign diagnostic number !USE DIAG_MOD, ONLY : AD03 !------------------------------------------ USE TIME_MOD, ONLY : GET_TS_CHEM USE TRACER_MOD, ONLY : STT # include "CMN_SIZE" ! Size parameters # include "CMN_DIAG" ! ND03 ! Local variables INTEGER :: I, J, L REAL*8 :: DTCHEM, KRATE, LOSS_FACTOR, Kr_LOST !================================================================= ! CHEMKr85 begins here! !================================================================= ! Chemistry timestep [s] DTCHEM = GET_TS_CHEM() * 60d0 ! The decay for 85Kr is calculated by: dC/dt = -kC ! where k = 1/15.52yr = 2.042E-9 s^-1 KRATE = 2.042d-9 ! Multiplication factor to compute tracer lost LOSS_FACTOR = 1d0 - EXP( -2.042d-9 * DTCHEM ) !================================================================= ! Apply radioactive decay to Kr85 tracer !================================================================= !$OMP PARALLEL DO !$OMP+DEFAULT( SHARED ) !$OMP+PRIVATE( I, J, L, Kr_LOST ) DO L = 1, LLPAR DO J = 1, JJPAR DO I = 1, IIPAR ! Amount of Kr85 lost to radioactive decay [kg] Kr_LOST = STT(I,J,L,1) * LOSS_FACTOR ! Prevent underflow IF ( Kr_LOST < SMALLNUM ) Kr_LOST = 0d0 ! Subtract Kr85 lost from the tracer array STT(I,J,L,1) = STT(I,J,L,1) - Kr_LOST !------------------------------------------------------------- ! Prior to 12/7/04: ! Need to reassign the diagnostic number !! Archive Kr85 lost by decay [kg] in ND04 diagnostic !IF ( ND03 > 0 ) THEN ! AD03(I,J,L,2) = AD03(I,J,L,2) + Kr_LOST !ENDIF !------------------------------------------------------------- ENDDO ENDDO ENDDO !$OMP END PARALLEL DO ! Return to calling program END SUBROUTINE CHEMKr85 !------------------------------------------------------------------------------ END MODULE Kr85_MOD