!$Id: chemmain.F90 125 2009-04-21 10:22:57Z acosce $ !! ========================================================================= !! INCA - INteraction with Chemistry and Aerosols !! !! Copyright Laboratoire des Sciences du Climat et de l'Environnement (LSCE) !! Unite mixte CEA-CNRS-UVSQ !! !! Contributors to this INCA subroutine: !! !! Didier Hauglustaine, LSCE, hauglustaine@cea.fr !! Stacy Walters, NCAR, stacy@ucar.edu !! C. Textor, LSCE !! !! Anne Cozic, LSCE, anne.cozic@cea.fr !! Yann Meurdesoif, LSCE, yann.meurdesoif@cea.fr !! !! This software is a computer program whose purpose is to simulate the !! atmospheric gas phase and aerosol composition. The model is designed to be !! used within a transport model or a general circulation model. This version !! of INCA was designed to be coupled to the LMDz GCM. LMDz-INCA accounts !! for emissions, transport (resolved and sub-grid scale), photochemical !! transformations, and scavenging (dry deposition and washout) of chemical !! species and aerosols interactively in the GCM. Several versions of the INCA !! model are currently used depending on the envisaged applications with the !! chemistry-climate model. !! !! This software is governed by the CeCILL license under French law and !! abiding by the rules of distribution of free software. You can use, !! modify and/ or redistribute the software under the terms of the CeCILL !! license as circulated by CEA, CNRS and INRIA at the following URL !! "http://www.cecill.info". !! !! As a counterpart to the access to the source code and rights to copy, !! modify and redistribute granted by the license, users are provided only !! with a limited warranty and the software's author, the holder of the !! economic rights, and the successive licensors have only limited !! liability. !! !! In this respect, the user's attention is drawn to the risks associated !! with loading, using, modifying and/or developing or reproducing the !! software by the user in light of its specific status of free software, !! that may mean that it is complicated to manipulate, and that also !! therefore means that it is reserved for developers and experienced !! professionals having in-depth computer knowledge. Users are therefore !! encouraged to load and test the software's suitability as regards their !! requirements in conditions enabling the security of their systems and/or !! data to be ensured and, more generally, to use and operate it in the !! same conditions as regards security. !! !! The fact that you are presently reading this means that you have had !! knowledge of the CeCILL license and that you accept its terms. !! ========================================================================= #include SUBROUTINE CHEMMAIN( & mmr , & nstep , & calday , & ncdate , & !not used ncsec , & !not used lat , & delt , & ps , & pmid , & pdel , & area , & oro , & tsurf , & ! not used albs , & zma , & phis , & cldfr , & cldfr_st , & cldfr_cv , & cldtop , & cldbot , & cwat , & flxrst , & flxrcv , & flxsst , & flxscv , & flxupd , & flxmass_w , & tfld , & sh , & ql , & ! variable en prevision de la chimie strato rh , & nx , & ny , & source) !----------------------------------------------------------------------- ! ! INCA -- INteractions with Chemistry in the Atmosphere ! ! Advances the volumetric mixing ratio forward one time step via a ! combination of explicit, EBI, QSSA, fully implicit, and/or rodas ! algorithms. ! ! Didier Hauglustaine and Stacy Walters, 2000. !----------------------------------------------------------------------- USE CHEM_MODS, ONLY : nas, invariants, hrates, hrates_cv, o3_prod, o3_loss, o3_st_flx, & so2_p_h2soh, so2_p_dmsoh, so2_p_dmsno3, so2_p_dmsooh, asmsam_p_dmsooh, & dmso_p_dmsoh, asso4m_p_so2oh, wet3d_so2, wet3d_dms, wet3d_dmso, wet3d_hno3, wet3d_nh3, & wet3d_noy, wet3d_h2o2, wet3d_hono, hno3_p_g, hno3_p_a, hno3_l_g, hno3_l_a, & nh3_l_g, nh3_l_a, csno3_p_a1, csno3_p_a2, cino3_p_a, & wet3d_asso4m, wet3d_asnh4m, wet3d_asno3m, wet3d_csno3m, wet3d_cino3m, & hono_p_g, hono_p_a, hono_l_g, hono_l_a, & asno3m_p_nh3hno3, asnh4m_p_nh3hno3, hno3_p_nh3hno3, nh3_p_nh3hno3, & ASAPp1a_p, ASAPp2a_p, ASARp1a_p, ASARp2a_p, pom_p_g, & co2_basprod, co2_nmhcprod, co2_radicalprod, & hno3_prod, hno3_loss, co_prod, co_loss, ch4_loss, n2o_loss, nadv_mass, & no_daytime, day_cnt , adv_mass, h2o, h2oc, extfrc, wetloss USE SPECIES_NAMES USE CHEM_TRACNM USE CHEM_CONS, ONLY : gravit, uma ! USE MASS_DIAGS USE AIRPLANE_SRC, ONLY : itrop USE SFLX USE TRANSPORT_CONTROLS USE INCA_DIM USE MOD_INCA_PARA USE MOD_GRID_INCA, ONLY : PLON_GLO USE RATE_INDEX_MOD USE INCA_WRITE_FIELD_P USE XIOS_INCA USE OXYDANT_COM USE PARAM_CHEM USE LIGHTNING #ifdef STRAT USE HETCHEM USE LGLIVED_MOD #endif IMPLICIT NONE !----------------------------------------------------------------------- ! ... Dummy arguments !----------------------------------------------------------------------- INTEGER, INTENT(in) :: nstep ! time index INTEGER, INTENT(in) :: lat ! latitude index (always 1 and place holder in LMDz) INTEGER, INTENT(in) :: nx, ny ! dyn grid INTEGER, INTENT(in) :: ncdate ! date at end present time step INTEGER, INTENT(in) :: ncsec ! seconds relative to ncdate INTEGER, INTENT(in) :: cldtop(PLON) ! cloud top level ( 1 ... PLEV ) INTEGER, INTENT(in) :: cldbot(PLON) ! cloud bot level ( 1 ... PLEV ) REAL, INTENT(in) :: calday ! time of year in days for midpoint time REAL, INTENT(in) :: delt ! timestep in seconds REAL, INTENT(inout) :: mmr(PLON,PLEV,PCNST) ! xported species ( mmr ) REAL, INTENT(in) :: ps(PLON) ! surface press ( pascals ) REAL, INTENT(in) :: pmid(PLON,PLEV) ! midpoint press ( pascals ) REAL, INTENT(in) :: pdel(PLON,PLEV) ! delta press across midpoints REAL, INTENT(in) :: area(PLON) ! grid cell area REAL, INTENT(in) :: oro(PLON) ! surface orography flag REAL, INTENT(in) :: tsurf(PLON) ! surface temperature REAL, INTENT(in) :: zma(PLON,PLEV) ! abs geopot height at midpoints ( m ) REAL, INTENT(in) :: phis(PLON) ! surf geopot REAL, INTENT(in) :: cldfr(PLON,PLEV) ! cloud fraction (all clouds) REAL, INTENT(in) :: cldfr_cv(PLON,PLEV) ! cloud fraction (convective clouds) REAL, INTENT(in) :: cldfr_st(PLON,PLEV) ! cloud fraction (large scale clouds) REAL, INTENT(in) :: cwat(PLON,PLEV) ! total cloud water (kg/kg) REAL, INTENT(in) :: flxrst(PLON,PLEVP) !liquid water flux (stratiform) kgH2O/m2/s REAL, INTENT(in) :: flxrcv(PLON,PLEVP) !solid water flux (stratiform) kgH2O/m2/s REAL, INTENT(in) :: flxsst(PLON,PLEVP) !liquid water flux (convection) kgH2O/m2/s REAL, INTENT(in) :: flxscv(PLON,PLEVP) !solid water flux (convection) kgH2O/m2/s REAL, INTENT(in) :: flxmass_w(PLON,PLEV) ! vertical mass flux REAL, INTENT(in) :: flxupd(PLON,PLEV) ! entrainment flux kgAIR/m2/s REAL, INTENT(in) :: tfld(PLON,PLEV) ! midpoint temperature REAL, INTENT(in) :: sh(PLON,PLEV) ! specific humidity ( kg/kg ) REAL, INTENT(in) :: ql(PLON,PLEV) ! eau liquide REAL, INTENT(in) :: rh(PLON,PLEV) ! relative humidity REAL, INTENT(in) :: albs(PLON) ! surface albedo REAL, INTENT(out) :: source(PLON,nbtrac) !----------------------------------------------------------------------- ! ... Local variables !----------------------------------------------------------------------- INTEGER :: i, k, m, n, it,j INTEGER :: klev # if GRPCNT != 0 REAL :: group_ratios(PLON,PLEV,GRPCNT) # endif # if NCOL != 0 REAL :: col_dens(PLON,PLEV,NCOL) # else REAL :: col_dens(1) # endif # if RELCNT != 0 REAL :: rel_ratios(PLON,PLEV,RELCNT) # endif # if HETCNT != 0 REAL, ALLOCATABLE, SAVE :: het_rates(:,:,:) REAL, ALLOCATABLE, SAVE :: het_rates_st(:,:,:) REAL, ALLOCATABLE, SAVE :: het_rates_cv(:,:,:) !$OMP THREADPRIVATE(het_rates,het_rates_st,het_rates_cv) # else REAL :: het_rates(1) = 0.0 # endif REAL :: vmr(PLON,PLEV,PCNST) ! xported species ( vmr ) REAL :: reaction_rates(PLON,PLEV,RXNCNT) REAL, DIMENSION(PLON,PLEV) :: h2ovmr ! water vapor volume mixing ratio REAL, DIMENSION(PLON,PLEV) :: mbar ! mean wet atmospheric mass ( amu ) REAL, DIMENSION(PLON,PLEV) :: zmid ! midpoint geopotential in km REAL, DIMENSION(PCNST) :: vprodj ! volume production REAL, DIMENSION(PCNST) :: vlossj ! volume loss REAL, DIMENSION(PLON) :: zen_angle REAL, DIMENSION(PLON) :: loc_angle REAL :: sunon(PLON) REAL :: sunoff(PLON) REAL :: delt_inverse REAL :: zelev(PLON) REAL :: zalt(PLON,PLEV) LOGICAL :: polar_night(PLON) LOGICAL :: polar_day(PLON) LOGICAL :: zangtz(PLON) REAL :: tfld_lim(PLON,PLEV) REAL :: albs_lim(PLON) REAL :: tfld_glo(PLON_GLO,PLEV) REAL :: pmid_glo(PLON_GLO,PLEV) REAL, DIMENSION(PLON,PLEV,RXNCNT) :: reacflux LOGICAL,SAVE :: first = .TRUE. !$OMP THREADPRIVATE(first) !----------------------------------------------------------------------- ! ... Function interface !----------------------------------------------------------------------- REAL :: TSECND ! Exemple utilisation de writefied pour le debuggage ! DO i=1,PCNST ! CALL writefield_p('q_'//tracnam(i), mmr(:,:,i), PLEV) ! ENDDO reaction_rates(:,:,:) = 0.0 IF (first) THEN ALLOCATE(het_rates(PLON,PLEV,HETCNT)) ALLOCATE(het_rates_st(PLON,PLEV,HETCNT)) ALLOCATE(het_rates_cv(PLON,PLEV,HETCNT)) het_rates(:,:,:) = 0.0 het_rates_st(:,:,:) = 0.0 het_rates_cv(:,:,:) = hrates_cv(:,:,:) first=.FALSE. ENDIF delt_inverse = 1. / delt !----------------------------------------------------------------------- ! ... T and albedo adjustment for photorates !----------------------------------------------------------------------- tfld_lim(:,:) = MAX(tfld(:,:),226.) albs_lim(:) = MIN(albs(:),0.5) # if PHTCNT != 0 !----------------------------------------------------------------------- ! ... Calculate parameters for diurnal geometry !----------------------------------------------------------------------- CALL DIURNAL_GEOM( & lat, calday, polar_night, polar_day, & sunon, sunoff, loc_angle, zen_angle, & zangtz) # endif !----------------------------------------------------------------------- ! ... Initialize xform between mass and volume mixing ratios !----------------------------------------------------------------------- CALL INTI_MR_XFORM( sh, mbar ) !----------------------------------------------------------------------- ! ... Xform from mmr to vmr !----------------------------------------------------------------------- CALL MMR2VMR( vmr, mmr, mbar ) !----------------------------------------------------------------------- ! ... Xform water vapor from mmr to vmr and adjust in stratosphere !----------------------------------------------------------------------- #ifdef STRAT CALL ADJH2O( h2o, sh, ql, h2oc, mbar, vmr ) #else CALL ADJH2O( h2ovmr, sh,ql, h2oc, mbar, vmr ) #endif !----------------------------------------------------------------------- ! ... Xform geopotential height from m to km !----------------------------------------------------------------------- zelev(:) = 1.e-3 * phis(:PLON) / gravit zmid(:,:) = 1.e-3 * zma(:PLON,:) / gravit DO i = 1, PLON DO k = 1, PLEV zalt(i,k) = zmid(i,k) + zelev(i) END DO END DO # if NFS != 0 !----------------------------------------------------------------------- ! ... Set the "invariants" !----------------------------------------------------------------------- CALL SETINV(invariants, tfld, h2ovmr, pmid) # endif #ifndef DUSS # if NCOL != 0 !----------------------------------------------------------------------- ! ... Set the column densities at the upper boundary !----------------------------------------------------------------------- CALL SET_UB_COL(col_dens, vmr, invariants, pdel) # endif #endif !----------------------------------------------------------------------- ! ... Set rates for "tabular" and user specified reactions !----------------------------------------------------------------------- # if SETRXNCNT != 0 CALL SETRXT( reaction_rates, tfld ) # endif #ifndef DUSS # if USRRXNCNT != 0 CALL USRRXT( & reaction_rates(:,:,:) , & tfld , & invariants , & rh , & ps , & pmid , & mmr , & invariants(1,1,INDEXM) , & lat , & zen_angle , & vmr ) # endif #endif # if GASCNT != 0 # ifdef RADJFLAG ! changement d'unite des taux de reactions des reactions a 2 reactants ! ils sont calcules dans usrrxt en s-1/(mol.cm-3) et le code dans ! exp_sol et imp_sol attends des s-1 ! ce n'est pas le cas des reactions a 1 reactant dont les taux ! sont directement calcule en s-1 par usrrxt (so4aerrxt etc.) CALL ADJRXT(reaction_rates, invariants, invariants(1,1,INDEXM)) # endif # endif # if PHTCNT != 0 !----------------------------------------------------------------------- ! ... Compute the photolysis rates !----------------------------------------------------------------------- # if NCOL != 0 !----------------------------------------------------------------------- ! ... Set the column densities !----------------------------------------------------------------------- CALL SETCOL(col_dens, vmr, invariants, pdel) # endif !----------------------------------------------------------------------- ! ... Calculate the photodissociation rates !----------------------------------------------------------------------- CALL PHOTO(reaction_rates, pmid, pdel, & tfld, zalt, col_dens, zen_angle, & albs, cwat, cldfr, sunon, sunoff, & polar_night, polar_day) !----------------------------------------------------------------------- ! ... Adjust the photodissociation rates !----------------------------------------------------------------------- CALL PHTADJ(reaction_rates, invariants, invariants(1,1,INDEXM)) # endif # if HETCNT != 0 !----------------------------------------------------------------------- ! ... Compute the heterogeneous rates at time = t(n+1) !----------------------------------------------------------------------- ! ... For stratiform precipitation CALL SETHET( & het_rates_st ,& pmid ,& pdel ,& lat ,& zmid ,& tfld ,& delt ,& invariants(1,1,INDEXM) ,& flxrst ,& flxsst ,& flxupd ,& cldtop ,& cldbot ,& cldfr_st ,& 1 ,& vmr ) ! ... For convective precipitation CALL SETHET( & het_rates_cv ,& pmid ,& pdel ,& lat ,& zmid ,& tfld ,& delt ,& invariants(1,1,INDEXM) ,& flxrcv ,& flxscv ,& flxupd ,& cldtop ,& cldbot ,& cldfr_cv ,& 2 ,& vmr ) ! ... Total removal rate het_rates = het_rates_st + het_rates_cv hrates = het_rates hrates_cv = het_rates_cv # endif # if EXTCNT != 0 !----------------------------------------------------------------------- ! ... Compute the extraneous forcings !----------------------------------------------------------------------- !Initialize for this timeslice DO m = 1,EXTCNT extfrc(:,:,m) = 0. END DO #ifndef GES CALL SETEXT_LGHT( & calday, extfrc, lat, zmid, & area, cldtop, vmr, tfld, & col_dens, zen_angle, delt, & invariants) #ifdef NMHC CALL SETEXT_O3L( & calday, extfrc, lat, & zmid, area, cldtop, & vmr, tfld, col_dens, & zen_angle, delt, & invariants) CALL SETEXT_BBG( & calday, extfrc, lat, & zmid, area, cldtop, & vmr, tfld, col_dens, & zen_angle, delt, pdel, & invariants) IF (flag_plane .GE. 1) THEN CALL SETEXT_AIR( & calday, extfrc, lat, & zmid, area, cldtop, & vmr, tfld, col_dens, & zen_angle, delt, & invariants) ENDIF #endif #else DO m = 1,EXTCNT extfrc(:,:,m) = 0. END DO #endif # endif # if GRPCNT != 0 !----------------------------------------------------------------------- ! ... Set the group ratios !----------------------------------------------------------------------- CALL SET_GRP_RATIOS( & group_ratios, reaction_rates(:,:,:), & reaction_rates, vmr, mmr, nas, mbar, invariants) !----------------------------------------------------------------------- ! ... Vertical mass flux !----------------------------------------------------------------------- DO i = 1, PLON !klev = itrop(i) klev = 11 o3_st_flx(:)=flxmass_w(:,klev)*mmr(:,klev,id_OX)/uma/nadv_mass(id_O3) END DO # endif # if RELCNT != 0 !----------------------------------------------------------------------- ! ... Set the relationship ratios !----------------------------------------------------------------------- # if GRPCNT != 0 CALL SET_REL_RATIOS(& rel_ratios , PLON, PLEV, & RELCNT, group_ratios, & GRPCNT, reaction_rates(1,1,PHTCNTP1), & GASCNT, reaction_rates, PHTCNT, vmr, & PCNST, invariants, NFS) #else CALL SET_REL_RATIOS(& rel_ratios , PLON, PLEV, & RELCNT, reaction_rates(1,1,PHTCNTP1), & GASCNT, reaction_rates, PHTCNT, vmr, & PCNST, invariants, NFS) #endif #endif #if RELCNT != 0 || GRPCNT != 0 !----------------------------------------------------------------------- ! ... Modify the reaction rate of any reaction ! with group member or proportional reactant(s) !----------------------------------------------------------------------- #if RELCNT != 0 CALL RXT_MOD(reaction_rates, het_rates, rel_ratios, group_ratios) #else CALL RXT_MOD(reaction_rates, het_rates, group_ratios) #endif #endif !======================================================================= ! ... Call the class solution algorithms !======================================================================= #if CLSCNT1 != 0 !----------------------------------------------------------------------- ! ... Solve for "explicit" species !----------------------------------------------------------------------- CALL EXP_SOL(vmr, reaction_rates, het_rates, extfrc, & #ifdef NMHC co_prod, co_loss, ch4_loss, n2o_loss, invariants(:,:,INDEXM),& #endif nstep, delt) #endif #ifndef DUSS # if CLSCNT4 != 0 !----------------------------------------------------------------------- ! ... Solve for Implicit species !----------------------------------------------------------------------- CALL IMP_SOL( & vmr, & reaction_rates, & het_rates, & extfrc, & o3_prod, & o3_loss, & #ifdef NMHC co2_basprod, & co2_nmhcprod, & co2_radicalprod, & hno3_prod, & hno3_loss, & #endif #ifdef AER so2_p_h2soh, & so2_p_dmsoh, & so2_p_dmsno3, & so2_p_dmsooh, & asmsam_p_dmsooh, & dmso_p_dmsoh, & asso4m_p_so2oh, & wet3d_so2, & wet3d_dms, & wet3d_dmso, & #ifdef NMHC pom_p_g, & hno3_p_g, & hno3_p_a, & hno3_l_g, & hno3_l_a, & hono_p_g, & hono_p_a, & hono_l_g, & hono_l_a, & nh3_l_g, & nh3_l_a, & csno3_p_a1, & csno3_p_a2, & cino3_p_a, & wet3d_h2o2, & wet3d_hono, & #endif wet3d_hno3, & wet3d_nh3, & wet3d_noy, & wet3d_asso4m, & wet3d_asnh4m, & wet3d_asno3m, & wet3d_csno3m, & wet3d_cino3m, & #endif pdel, & #ifdef AERONLY invariants, & #endif invariants(:,:,INDEXM), & nstep, & delt_inverse, & lat, wetloss ) # endif #endif #ifdef STRAT !----------------------------------------------------------------------- ! ... Sedimentation in the stratosphere !----------------------------------------------------------------------- CALL STRSEDCALC( invariants(:,:,INDEXM), pdel, vmr, delt,h2o, h2oc ) #endif DO i = 1, PLON #ifdef AER #if !defined(AERONLY) wflux(i,id_ASSO4M) = SUM(wet3d_asso4m(i,:)) wflux(i,id_ASNH4M) = SUM(wet3d_asnh4m(i,:)) wflux(i,id_ASNO3M) = SUM(wet3d_asno3m(i,:)) wflux(i,id_CSNO3M) = SUM(wet3d_csno3m(i,:)) wflux(i,id_CINO3M) = SUM(wet3d_cino3m(i,:)) #endif wfluxso2(i) = SUM(wet3d_so2(i,:)) wfluxdms(i) = SUM(wet3d_dms(i,:)) wfluxdmso(i) = SUM(wet3d_dmso(i,:)) wfluxnoy(i) = SUM(wet3d_noy(i,:)) wfluxnh3(i) = SUM(wet3d_nh3(i,:)) wfluxhno3(i) = SUM(wet3d_hno3(i,:)) #endif ENDDO #if !defined(DUSS) && defined(AER) !DH Nitrates formation CALL AERTHERM( & delt , & tfld , & rh , & pmid , & invariants(:,:,INDEXM) , & asno3m_p_nh3hno3 , & asnh4m_p_nh3hno3 , & hno3_p_nh3hno3 , & nh3_p_nh3hno3 , & mmr , & vmr ) #endif !----------------------------------------------------------------------- ! ... Check for negative values and reset to zero !----------------------------------------------------------------------- CALL NEGTRC( 'After chemistry ', vmr ) IF (calcul_flux) THEN !----------------------------------------------------------------------- ! ... Compute the flux through each reaction !----------------------------------------------------------------------- CALL reac_flx(vmr, reaction_rates, invariants, reacflux) CALL xios_inca_change_context("inca") DO n=1, RXNCNT CALL xios_inca_send_field("flux_"//trim(reacname(n)), reacflux(:,:,n)) ENDDO CALL xios_inca_change_context("LMDZ") ENDIF !----------------------------------------------------------------------- ! ... Xform from vmr to mmr !----------------------------------------------------------------------- CALL VMR2MMR( vmr, mmr, & # if GRPCNT != 0 nas, & group_ratios, & # endif mbar ) #ifndef DUSS #ifdef AER CALL GASTOPARTICLE(delt,mmr,mbar,invariants(:,:,INDEXM)) #endif #endif #if defined(NMHC) && defined(AER) CALL soa_main(mmr, & mbar, & pmid, & tfld, & ASAPp1a_p, & ASAPp2a_p, & ASARp1a_p, & ASARp2a_p, & delt) #endif !DH Update the drydep flux with latest surface mixing ratio CALL DRYDEP (pmid, & tfld, & mmr) ! calcul des sources DO it=1,nbtrac DO k = 1, PLON source(k,it) = eflux(k,it)-dflux(k,it) ENDDO ENDDO flux_source(:,:) = source(:,:) #if defined(NMHC) && defined(AER) CALL xios_inca_change_context("inca") ! CALL xios_inca_send_field("O3_surf",vmr(:,1,id_O3)) ! Define surface concentration of O3 --- CALL xios_inca_send_field("SO2_surf",vmr(:,1,id_SO2)) ! Define surface concentration of SO2 --- CALL xios_inca_send_field("NO2_surf",vmr(:,1,id_NO2)) ! Define surface concentration of NO2 --- CALL xios_inca_send_field("NH3_surf",vmr(:,1,id_NH3)) ! Define surface concentration of NH3 --- CALL xios_inca_send_field("CH2O_surf",vmr(:,1,id_CH2O)) ! Define surface concentration of CH2O --- CALL xios_inca_send_field("OH_surf",vmr(:,1,id_OH)) ! Define surface concentration of OH --- CALL xios_inca_change_context("LMDZ") ! #endif END SUBROUTINE CHEMMAIN SUBROUTINE endrun !---------------------------------------------------------------------- ! ... Abort the model !---------------------------------------------------------------------- IMPLICIT NONE EXTERNAL abort CALL abort END SUBROUTINE endrun