text:
! <tsgrid.f90 - A component of the City-scale
!                 Chemistry Transport Model EPISODE-CityChem>
!*****************************************************************************!
!*
!* EPISODE - An urban-scale air quality model
!* ==========================================
!* Copyright (C) 2018  NILU - Norwegian Institute for Air Research
!*                     Instituttveien 18
!*                     PO Box 100
!*                     NO-2027 Kjeller
!*                     Norway
!*
!*                     Contact persons: Gabriela Sousa Santos - gss@nilu.no
!*                                      Paul Hamer - pdh@nilu.no
!*
!* Unless explicitly acquired and licensed from Licensor under another license,
!* the contents of this file are subject to the Reciprocal Public License ("RPL")
!* Version 1.5, https://opensource.org/licenses/RPL-1.5 or subsequent versions as
!* allowed by the RPL, and You may not copy or use this file in either source code
!* or executable form, except in compliance with the terms and conditions of the RPL.
!*
!* All software distributed under the RPL is provided strictly on an "AS IS" basis,
!* WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESS OR IMPLIED, AND LICENSOR HEREBY
!* DISCLAIMS ALL SUCH WARRANTIES, INCLUDING WITHOUT LIMITATION, ANY WARRANTIES OF
!* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, QUIET ENJOYMENT, OR NON-INFRINGEMENT.
!* See the RPL for specific language governing rights and limitations under the RPL.
!*
!* ==========================================
!* The dispersion model EPISODE (Grønskei et. al., 1993; Larssen et al., 1994;
!* Walker et al., 1992, 1999; Slørdal et al., 2003, 2008) is an Eulerian grid model
!* with embedded subgrid models for calculations of pollutant concentrations resulting
!* from different types of sources (area-, line- and point sources). EPISODE solves the
!* time dependent advection/-diffusion equation on a 3 dimensional grid.
!* Finite difference numerical methods are applied to integrate the solution forward in time.
!* It also includes extensions as the implementation of a simplified EMEP photochemistry
!* scheme for urban areas (Walker et al. 2004) and a scheme for Secondary Organic Aerosol
!* implemented by Håvard Slørdal
!*
!* Grønskei, K.E., Walker, S.E., Gram, F. (1993) Evaluation of a model for hourly spatial
!*    concentrations distributions. Atmos. Environ., 27B, 105-120.
!* Larssen, S., Grønskei, K.E., Gram, F., Hagen, L.O., Walker, S.E. (1994) Verification of
!*    urban scale time-dependent dispersion model with sub-grid elements in Oslo, Norway.
!*    In: Air poll. modelling and its appl. X. New York, Plenum Press.
!* Slørdal, L.H., McInnes, H., Krognes, T. (2008): The Air Quality Information System AirQUIS.
!*    Info. Techn. Environ. Eng., 1, 40-47, 2008.
!* Slørdal, L.H., Walker, S.-E., Solberg, S. (2003) The Urban Air Dispersion Model EPISODE
!*    applied in AirQUIS. Technical Description. NILU TR 12/2003. ISBN 82-425-1522-0.
!* Walker, S.E., Grønskei, K.E. (1992) Spredningsberegninger for on-line overvåking i Grenland.
!*    Programbeskrivelse og brukerveiledning. Lillestrøm,
!*    Norwegian Institute for Air Research (NILU OR 55/92).
!* Walker, S.E., Slørdal, L.H., Guerreiro, C., Gram, F., Grønskei, K.E. (1999) Air pollution
!*    exposure monitoring and estimation. Part II. Model evaluation and population exposure.
!*    J. Environ. Monit, 1, 321-326.
!* Walker, S.-E., Solberg, S., Denby, B. (2003) Development and implementation of a
!*    simplified EMEP photochemistry scheme for urban areas in EPISODE. NILU TR 13/2013.
!*    ISBN 82-425-1524-7
!*****************************************************************************!

      subroutine TSGRID(I)

! *** ThIS subroutine updates, calculates and writes grid model concentrations for one timestep.
! ----------------------------------------------------------------------------------
! Based on:
! Version Episode 5.5 (May29, 2012) prepared for BB-Stand-Alone
! Original source code of EPISODE by Sam-Erik Walker (NILU)
!
! Sam-Erik Walker
! Norwegian Institute for Air Research (NILU)
! Instituttveien 18 P.O. Box 100
! N-2027 Kjeller, NORWAY
! Tel: +47 63898000 Fax: +47 63898050
! E-mail: sam-erik.walker@nilu.no
!
! ----------------------------------------------------------------------------------
! ----------------------------------------------------------------------------------
! REVISION HISTORY
!
!           2016  M. Karl: commented all emep69 SOA preprocessor definitions
!           2016  M. Karl: commented sub grid cell and nested subgrid treatment
!           2016  M. Karl: EPISODE newbc, newc and oldc files are not supported.
!           2016  M. Karl: new way to open and update BCONs
!           2017  M. Karl: restart with IC files possible
!           2017  M. Karl: c,dcdt,d2cdx2,d2cdy2,excz,bc defined as double precision
!
!           2019 Dec. Kang @ CARES, write out 3D emission data (instantanous) to .dat file
!!          here, two height (Z) are given:
!                 Z(IZ) - vertical layer height,
!                 Z_act - height after considering topograhpy
!           2020 Mar. Kang @ CARES, write out 3D emission data (including BC, instantanous)
!                Changes include:
!                1. C(NC,NX+2,NY+2,NZ+1)
!!               2. location of main grid point is changed from sw and top to central point
! ----------------------------------------------------------------------------------

      use mod_util
      use mod_main
      use mod_site
      use mod_time
      use mod_mete
      use mod_conc
      use mod_grid
!MSK start
      use mod_writenc
!MSK end

!!DEC$ IF DEFINED (emep68)
!
! *** Additional EMEP1_SOA modules:
!
!      use OrganicAerosol_ml, only : OrganicAerosol,temp_soa, itemp_soa,
!     &        xn,
!     &        Grid_SOA_Fgas,Grid_avg_mw,
!     &        Grid_M0,
!     &        aer_ug_OC,    aer_ugC_OC,    aer_ug_POM,  aer_ugC_POM,
!     &        aer_ug_WOOD,  aer_ugC_WOOD,  aer_ug_FFUEL,aer_ugC_FFUEL,
!     &        aer_ug_BGNDOC,aer_ugC_BGNDOC,
!     &        aer_ug_ASOA,  aer_ugC_ASOA,  aer_ug_BSOA, aer_ugC_BSOA,
!     &        my_last_call
!
!      use My_SOA_ml,         only : FIRST_NONVOL, LAST_NONVOL,
!     &                              FIRST_SOA, LAST_SOA, BGND_SPEC
!
!!DEC$ ENDIF

      implicit none

! *** Start subroutine declarations: **************************************************************

! *** Global variables:

      INTEGER :: I

! *** I - Indicator of pre (0) or post (1) processing

! *** Local variables:

      INTEGER :: IC
      INTEGER :: IX
      INTEGER :: IY
      INTEGER :: IZ
!MSK      LOGICAL :: ATDATE
!MSK      LOGICAL :: ATTIME

! *** IC     - Compound index
! *** IX     - X-direction index
! *** IY     - Y-direction index
! *** IZ     - Vertical layer index
! *** ATDATE - If at given date          then true else false
! *** ATTIME - If at given date and time then true else false

      real             :: CMINV
      real             :: molec2ug
      double precision :: HOR_VFLUX,HOR_MFLUX,VERT_VFLUX,VERT_MFLUX,  &
                        TOT_VFLUX,TOT_MFLUX

      integer :: AV_GRIDAI
      real :: rsum_DCDT

!MSK start
!     Variables for netCDF output
      character(len=10) :: namefield
      character(len=10) :: unitname
      character(len=23) :: validity

      integer :: Nhh_in

      integer, dimension(4,1) :: mdate

      logical :: dopacking  = .false.
      logical :: dofloat    = .false.
      logical :: domirror   = .false.

      double precision        :: field3D(NX,NY,NZ)
      !double precision        :: field3D(NY,NX,NZ)
!!=====================================================================================
!! write out main grid emission data at last timestep of the simulation to .dat file
!! added by Kang @ CARES, Dec. 23, 2019
!!=====================================================================================
      double precision, allocatable :: field3D_output(:,:,:,:)
      integer :: II
      real,   allocatable :: mxm_i(:,:)
      real,   allocatable :: mxm_j(:,:)

!$$$$$$       real,   allocatable :: mxm_i_new(:,:)
!$$$$$$       real,   allocatable :: mxm_j_new(:,:)
!$$$$$$       real,   allocatable :: Znew(:)
!$$$$$$       real,   allocatable :: Znew0(:)
! Allocate
        if (.not. allocated(mxm_i) )    allocate(mxm_i(NY,NX))
        if (.not. allocated(mxm_j) )    allocate(mxm_j(NY,NX))

! Allocate
!$$$$$$         if (.not. allocated(mxm_i_new) )    allocate(mxm_i_new(NY+2*NBCY,NX+2*NBCX))
!$$$$$$         if (.not. allocated(mxm_j_new) )    allocate(mxm_j_new(NY+2*NBCY,NX+2*NBCX))
!$$$$$$         if (.not. allocated(Znew) )    allocate(Znew(NZ+NBCZ))
!$$$$$$         if (.not. allocated(Znew0) )    allocate(Znew0(NZ))
!!=====================================================================================
!MSK end

!!DEC$ IF DEFINED (emep68)
! ***  Local EMEP1_SOA-variables
!
!      logical, parameter :: DEBUG = .true.
!      real               :: ppb   = 2.5e10  ! ????
!
!!DEC$ ENDIF

! *** End subroutine declarations. ****************************************************************

      if (averaged_output) then
        AV_GRIDAI = GRIDAI
      else
        AV_GRIDAI = 0
      endif

      if (messfe .AND. ATTIME(BDAT) .AND. ISH == 1 .AND. ITS == 1 .AND. I == 1) then
        write (messun,*)
        if (averaged_output) then
           write (messun,'(A38)') &
            'TSGRID: The CM-output is NTS-averaged.'
        else
           write (messun,'(A39)') &
            'TSGRID: The CM-output is Instantaneous.'
        endif
        write (messun,*)
      endif

! *************************************************************************************************
!_NEST_Start:
!
! *** IMETHOD: Horizontal interpolation method when calculating concentration values for the NEST
! ***          domain (Declared in: mod_grid).
!
! ***   IMETHOD equal to 1 means: BILINEAR INTERPOLATION
! ***   IMETHOD equal to 2 means: BESSEL   INTERPOLATION
! ***   IMETHOD equal to 3 means: NEAREST GRIDPOINT VALUE
!
        IMETHOD  = 1

! *** IZMETHOD: Vertical interpolation method when calculating concentration values for the NEST
! ***           domian (Declared in: mod_grid).
!
! ***   IZMETHOD equal to 0 means: NO INTERPOLATION (i.e. is used when
! ***                              we have identical vertical resolution.)
! ***   IZMETHOD equal to 1 means: LINEAR INTERPOLATION

        IZMETHOD = 1

!_NEST_End.

! *************************************************************************************************

! *** Perform pre- or postprocessing calculations:
      IF (I .EQ. 0) GOTO 100
      IF (I .EQ. 1) GOTO 200

  100 CONTINUE

! *** Preprocessing calculations:
! *** ***************************

!!DEC$ IF DEFINED (emep68)
!
!! ***   The SOA-module is only applied if NC >= 68:
!
!!_TEST      if (NC >= 68 ) then
!
!        if (ATTIME(BDAT) .AND. ISH == 1 .AND. ITS == 1) then
!
!! ***     Allocate the temperature and concentration profile arrays used
!! ***     in OrganicAerosol at the start of the simulation:
!
!! ***     The three arrays below was originally declared in: module SOA_Inputs_ml.
!          if (.NOT. allocated(temp_soa))  allocate(temp_soa(1:nz))
!          if (.NOT. allocated(itemp_soa)) allocate(itemp_soa(1:nz))
!          if (.NOT. allocated(xn))        allocate(xn(1:LAST_SOA,1:nz))
!
!! ***     The three arrays below should be available outside of
!! ***     the SOA-module.
!!          if (.NOT. allocated(Grid_SOA_Fgas))
!!     &      allocate(Grid_SOA_Fgas(FIRST_SOA:LAST_SOA,nx,ny,1:nz))
!!          if (.NOT. allocated(Grid_avg_mw))
!!     &      allocate(Grid_avg_mw(nx,ny,1:nz))
!
!          if (.NOT. allocated(Grid_M0)) allocate(Grid_M0(nx,ny,1:nz))
!
!! ***     The compound numbers below are declared in 'mod_conc.for':
!! ***     See: subroutine SendConcCompoundsData(pinc,piacmpnd,piacunit,piacout,pxacmolw,pxathalf)
!!
!! ***     The numbers below are not equal to "icmpnd(ic)", but instead: "icmpnd(ic) - 4"
!          ic_org_PM2p5  = LAST_SOA      + 1
!          ic_aerOM      = ic_org_PM2p5  + n_advect
!          ic_aerOM_C    = ic_aerOM      + 1
!          ic_aerPOM     = ic_aerOM_C    + 1
!          ic_aerPOM_C   = ic_aerPOM     + 1
!          ic_aerWOOD    = ic_aerPOM_C   + 1
!          ic_aerWOOD_C  = ic_aerWOOD    + 1
!          ic_aerFFUEL   = ic_aerWOOD_C  + 1
!          ic_aerFFUEL_C = ic_aerFFUEL   + 1
!          ic_aerBGND    = ic_aerFFUEL_C + 1
!          ic_aerBGND_C  = ic_aerBGND    + 1
!          ic_aerASOA    = ic_aerBGND_C  + 1
!          ic_aerASOA_C  = ic_aerASOA    + 1
!          ic_aerBSOA    = ic_aerASOA_C  + 1
!          ic_aerBSOA_C  = ic_aerBSOA    + 1
!
!        endif  ! if (ATTIME(BDAT) .AND. ISH == 1 .AND. ITS == 1).
!
!_TEST      endif  ! if (NC >= 68 ).
!
!!DEC$ ENDIF (emep68)

! *** UPDATE BOUNDARY CONDITIONS:
! *** ***************************
!
! *** Standard practice: Updating the BC array once every hour
! *** or whenever ITS = 1 (for instance each half-hour):
!
!MSK: new BC values
!MSK:  Note: rbcon opens the netcdf files for main grid C output

      IF (ITS .EQ. 1)  call rbcon

!MSK start debug
!debug          do IZ=1,NZ+NBCZ
!debug            print *,'tsgrid BC(Z) ', IZ, BC(1,2,2,IZ)
!debug          enddo
!MSK end debug

! *** TEST POSSIBILITIES:
!

!MKS start
!MSK: FIRST TIME STEP: USE BC(:,:,:,:) AS INITIAL CONDITION FOR C(:,:,:,:)

      IF (ATTIME(BDAT) .AND. ITS .EQ. 1) THEN

         DO IC = 1,NC
           DO IZ=1,NZ+NBCZ
             DO IY=1,NY+2*NBCY
               DO IX=1,NX+2*NBCX
                 C(IC,IX,IY,IZ) = BC(IC,IX,IY,IZ)
               ENDDO
             ENDDO
           ENDDO
         ENDDO

      ENDIF

!MSK start debug
!debug           do IZ=1,NZ+NBCZ
!debug             print *,'tsgrid C(Z) 1', IZ, C(1,2,2,IZ)
!debug           enddo
!MSK end debug
!MSK end

!MSK:
! *** Only the main grid extent border is filled with BC values
!     BC(NC, NX+2*NBCX, NY+2*NBCY, NZ+NBCZ)
!     C(NC, NX+2*NBCX, NY+2*NBCY, NZ+NBCZ)
!     C(1,1) is same as BC(1,1)

!Fill the NZ+NBCZ layer
!        DO IX=1,NX+2*NBCX
!          DO IY=1,NY+2*NBCY
!            C(IC, IX, IY, NZ+NBCZ)     = BC(IC,IX,IY,NZ+NBCZ)
!          ENDDO
!        ENDDO
!Fill the East-West borders
!        DO K=1,NZ
!          DO IY=1,NY+2*NBCY
!            C(IC,  1, IY, K)      = BC(IC,  1, IY,K)
!            C(IC, NX+2*NBCX,IY,K) = BC(IC,NX+2*NBCX,IY,K)
!          ENDDO
!        ENDDO
!Fill the North-South borders
!        DO K=1,NZ
!          DO IX=1,NX+2*NBCX
!            C(IC, IX,  1, K)      = BC(IC,IX,  1,K)
!            C(IC, IX,NY+2*NBCY,K) = BC(IC,IX, NY+2*NBCY,K)
!          ENDDO
!        ENDDO
!
!MSK end

! *** Updating the boundary conditione every timestep:
!     CALL RBCON
!
! *** The IF-test below is activated if an analytical boundary solution is
! *** specified at the beginning of the simulation, and is zero thereafter:
!     IF (ATTIME(BDAT) .AND. ITS .EQ. 1)THEN
! ***   The three lines below can be activated alone if the analytic boundary
!       solution is to be updated every timestep.
!       DO IC = 1,NC
!         CALL AN_BCON(IC,0)
!       ENDDO
!     ELSE
!       DO IC = 1,NC
!         DO IZ=1,NZ+NBCZ
!           DO IY=1-NBCY,NY+NBCY
!             DO IX=1-NBCX,NX+NBCX
!               BC(IC,IX,IY,IZ) = 0.0
!             ENDDO
!           ENDDO
!         ENDDO
!       ENDDO
!     ENDIF
! *************************************************************************************************
!
! *** SPECIFY INITIAL CONDITIONS:
! *** ***************************
! *** Specify old grid model concentrations (from a previous simultion) at the beginning
! *** of the simuation:

      IF (ATTIME(BDAT) .AND. ISH .EQ. 1 .AND. ITS .EQ. 1) THEN

!MSK start
! ***   For RESTART: get initial concentrations from IC restart file
! ***   Replaces roldc from EPISODE

               call ricon

!!!!!!!!!!!!!!!!!!! for testing 3D concentration output file by Kang @ Dec.2, 2019
!             if(restart .eq. 1) then
!                open(unit=1113,file="../OUTPUT/3D_field_new.txt")  !!
!                open(unit=1112,file="../OUTPUT/3D_field_para.txt")  !!
!              do 350 IC=1,NC
!               if (IC .eq. 3) then
!                do 351 IZ=1,NZ
!                  do 352 IX = 1,NX
!                    do 353 IY = 1,NY

                     !!  field3D(IX,IY,IZ) = C(3,IX+1,IY+1,IZ)
                  !     field3D(IY,IX,IZ) = C(IC,IX+1,IY+1,IZ)
!                      write(1113,*) C(3,IX+1,IY+1,IZ)
!  351               continue
!  352             continue
!  353           continue
!               endif  !!! end if of IC=3
!  350          continue

!                 write(1112,*)"dump out 3D concentration file, NX=",NX,", NY=", NY, "NZ=",NZ, "Nhh_in=",Nhh_in, "mdate=",mdate, ", validity=",validity, ", dopacking=",dopacking,", dofloat=", dofloat, ", domirror=",domirror

!                close(1113)
!                close(1112)
!            endif
!!!!!!!!!!!!!!!!!!! end for testing 3D concentration output file by Kang @ Dec.2, 2019

!MSK start debug
!debug           do IZ=1,NZ+NBCZ
!debug             print *,'tsgrid C(Z) 2', IZ, C(1,2,2,IZ)
!debug           enddo
!MSK end debug
!MSK end

! ***  Calculate the concentrations for the NEST grid domain. These values are calculated
! ***  for the entire NEST grid, and are stored in the array:
!
! ***       BC_nest(NC,1-nbcx_nest:nx_nest+nbcx_nest,
! ***      &           1-nbcy_nest:ny_nest+nbcy_nest,
! ***      &                       nz_nest+nbcz_nest))

!MSK start  New subgrid BC values, only wanted if there are NEWBC files
!MSK     Now commented subrid BC and IC files
!MSK        if (NEWBCFE( 1 )) then
!MSK
!MSK          CALL WNEWBC
!MSK
!MSK       endif
!MSK end

!MSK start
! *** Subgrid nests are currently not used, commented below
! ***   Store the initial BC_nest variables in the array. Needed since the values
! ***   of BC_nest will be overwritten at the end of the present hour in the post-
! ***   processing treatment within this routine.
!
! ***       IC_BC_nest(nc,1-nbcx_nest:nx_nest+nbcx_nest,
! ***      &              1-nbcy_nest:ny_nest+nbcy_nest,
! ***      &                          nz_nest+nbcy_nest)

!MSK        IC_BC_nest = BC_nest

! ***   Store the internal values of IC_BC_nest in the array:
!
! ***       IC_C_nest(nc,nx_nest,ny_nest,nz_nest)
!MSK
!MSK        do ic = 1,nc
!MSK          do iz = 1,nz_nest
!MSK            do iy = 1,ny_nest
!MSK              do ix = 1,nx_nest
!MSK                IC_C_nest(ic,ix,iy,iz) = IC_BC_nest(ic,ix,iy,iz)
!MSK              enddo
!MSK            enddo
!MSK          enddo
!MSK        enddo

!MSK ***   Zero initialisation of the not used subgrid nest fields
        BC_nest(:,:,:,:) = 0.0
        IC_C_nest(:,:,:,:)  = 0.0
        IC_BC_nest(:,:,:,:) = 0.0
!MSK end

      ENDIF

!     Run model OK?

      IF (.NOT. RUNOK) THEN

!       Set grid model concentrations to missing data:

        print *,'tsgrid not runok'

        DO IC = 1,NC

          IF (C(IC,1,1,1) .EQ. MISS) THEN

!_NEST_Start:

            DO IZ=1,NZ+NBCZ
!MSK              DO IY=1-NBCY,NY+NBCY
!MSK                DO IX=1-NBCX,NX+NBCX
              DO IY=1,NY+2*NBCY
                DO IX=1,NX+2*NBCX
                  C(IC,IX,IY,IZ) = MISS
                ENDDO
              ENDDO
            ENDDO
!_NEST_End.

          ENDIF

        END DO

!MSK !       Write new grid model concentrations
!MSK
!MSK        IF (ATDATE(EDAT) .AND. ISH .EQ. NSH .AND. ITS .EQ. NTS) THEN
!MSK
!MSK          print *,'new conc file: ', NEWCFN(1)
!MSK          if ( NEWCFE( 1 ) ) then
!MSK
!MSK              call wnewc
!MSK
!MSK          endif
!MSK
!MSK          IF (MESSFE) WRITE (MESSUN,'(A)')   &
!MSK            'TSGRID: Finished call to: WNEWC'
!MSK        ENDIF
!MSK

        RETURN

      ENDIF  ! IF (.NOT. RUNOK)

!     If missing data then grid model concentrations is set equal to
!     background concentrations or if these are also missing set equal
!     to zero.

      DO IC = 1,NC

!       If missing data then set equal to background concentration:

        IF (C(IC,1,1,1) .EQ. MISS) THEN

!_NEST_Start:
          IF (MESSFE) WRITE (MESSUN,'(A67)')   &
            'TSGRID WARNING: UNDEF values calculated, background values applied.'

          DO IZ=1,NZ+NBCZ
!MSK index null not allowed
!MSK            DO IY=1-NBCY,NY+NBCY
!MSK              DO IX=1-NBCX,NX+NBCX
            DO IY=1, NY+2*NBCY
              DO IX=1, NX+2*NBCX
                C(IC,IX,IY,IZ) = BC(IC,IX,IY,IZ)
              ENDDO
            ENDDO
          ENDDO

!_NEST_End.

        ENDIF

!       If missing data then set equal to zero:

        IF (C(IC,1,1,1) .EQ. MISS) THEN

!_NEST_Start.
          IF (MESSFE) WRITE (MESSUN,'(A73)')  &
            'TSGRID WARNING: UNDEF values found in background, minimum values applied.'

          CMINV = 0.
          IF (CUNIT(IC) .EQ. 'mol/cm3') CMINV = 1.E-20

          DO IZ=1,NZ+NBCZ
!MSK index null not allowed
!MSK            DO IY=1-NBCY,NY+NBCY
!MSK              DO IX=1-NBCX,NX+NBCX
            DO IY=1, NY+2*NBCY
              DO IX=1, NX+2*NBCX
                C(IC,IX,IY,IZ) = CMINV
              ENDDO
            ENDDO
          ENDDO
!_NEST_End.

        ENDIF

      END DO

!MSK start debug
!debug          do IZ=1,NZ+NBCZ
!debug             print *,'tsgrid C(Z) 3', IZ, C(1,2,2,IZ)
!debug          enddo
!debug          print *,'tsgrid C before adv/dif ',C(1,2,2,1)
!MSK end debug

! *** Alternate direction of advection-diffusion operators every timestep:

      IF (MOD(ITS,2) .EQ. 1) THEN

! ***   Calculate vertical advection and diffusion:
!
! ***   At the moment:
! ***   NC       = 84 Total number of species in the model.
! ***   n_nochem = 16 Number of species not included in the chemistry.
! ***                 (Thus: 84 - 16 = 68 species included in the chemistry).
! ***   n_advect =  2 Number of n_nochem species that are advected.
! ***                 (these 2 are: 'ORGPM2p5' and  'ECfi'.)

        DO IC = 1,NC - n_nochem + n_advect

          if (messfe .AND. (ITS == 1 .OR. ITS == NTS) ) then
            write(messun,'(A,I4)') 'TSGRID: Preprocessing timestep ITS = ',ITS
              rsum_DCDT     = SUM(DCDT(IC,:,:,:))
           write(messun,'(A,I2,A,E12.5)')   &
              'TSGRID: Domain sum DCDT(IC=',IC,',:,:,:) = ',rsum_DCDT
          endif

          IF (GRIDVS(IC) .EQ. 0) THEN

! ***       No vertical advection or diffusion applied.
! ***       Then the emissions are introduced directly into
! ***       the appropriate grid cells.

            DO IZ = 1,NZ
              DO IY = 1,NY
                DO IX = 1,NX
!MSK                  C(IC,IX,IY,IZ) = C(IC,IX,IY,IZ) + DCDT(IC,IX,IY,IZ)*DT
                  C(IC,IX+1,IY+1,IZ) = C(IC,IX+1,IY+1,IZ) +  &
                                       DCDT(IC,IX,IY,IZ)*DT
                END DO
              END DO
            END DO

          ELSEIF (GRIDVS(IC) .EQ. 1) THEN

! ***       Applying the original combined vertical (upstream) advection
! ***       and explicit diffusion scheme:

              call adiv(ic)

          ELSEIF (GRIDVS(IC) .EQ. 2) THEN

! ***       Vertical diffusion by the explicit diffusion scheme:
              call v_dif(ic)

! ***       Vertical advection by the upstream scheme:
              call v_adv(ic)

          ELSEIF (GRIDVS(IC) .EQ. 3) THEN

! ***       Vertical diffusion by the implicit C-N diffusion scheme:
!MSK debug start pnc check
!debug       if((ic.eq.21).or.(ic.eq.22).or.(ic.eq.23) )then
!debug         print *,'tsgrid DCDT before vert diff', IC, DCDT(IC,12,10,1),C(IC,13,11,1)
!debug       endif
!MSK debug

              call difv_cn(ic)

!MSK debug start pnc check
!debug       if((ic.eq.21).or.(ic.eq.22).or.(ic.eq.23) )then
!debug         print *,'tsgrid DCDT before vert adv', IC, DCDT(IC,12,10,1),C(IC,13,11,1)
!debug       endif
!MSK debug

! ***       Vertical advection by the upstream scheme:
              call v_adv(ic)

!MSK debug start pnc check
!debug       if((ic.eq.21).or.(ic.eq.22).or.(ic.eq.23) )then
!debug         print *,'tsgrid DCDT after vert adv', IC, DCDT(IC,12,10,1),C(IC,13,11,1)
!debug       endif
!MSK debug

          ELSEIF (GRIDVS(IC) .EQ. 4) THEN

! ***       Only applying vertical diffusion by the implicit
! ***       C-N diffusion scheme:
              call difv_cn(ic)

          ELSEIF (GRIDVS(IC) .EQ. 5) THEN

! ***       Only applying vertical advection by the upstream scheme:

! ***       This should only be applied for test purposes
! ***       when no grid emissions are included.
              call v_adv(ic)

          ENDIF

         END DO  ! IC = 1,NC - n_nochem + n_advect

      ! print *,'tsgrid C after vertical ',C(1,2,2,1)

!_LHS_SOA_May_2007_End.

! ***   Calculate horisontal 2D advection:

          call advh

      ! print *,'tsgrid C after horizontal adv ',C(1,2,2,1)

! ***   Calculate horisontal 2D diffusion:

          call difh

       ! print *,'tsgrid C after horizontal dif ',C(1,2,2,1)

! ***   Calculate grid model output concentration and deposition

!MSK *** This is the connection between Eulerian grid concentration C
!MSK *** and the local main grid, receptor point, line receptor
!MSK *** and subgrid concentrations CM, CR, CLR, CS

! ***   ... in the main grid cells:
! ***       (At present: zero deposition calculated in CSUBGM.)

        call csubgm(AV_GRIDAI,GRIDAI,GRIDAI)

! ***   ... in the sub  grid cells
! ***       (At present: zero deposition calculated in CSUBGS.)

!MSK        CALL CSUBGS(AV_GRIDAI,GRIDAI,GRIDAI)

! ***   ... in the irregular receptor points
! ***       (At present: zero deposition calculated in CSUBGR.)

        call csubgr(AV_GRIDAI,GRIDAI,GRIDAI)

! ***   ... in the line source associated receptor points
! ***       (At present: zero deposition calculated in CSUBGL.)

        call csubgl(AV_GRIDAI,GRIDAI,GRIDAI)

! ***   ... in the grid model

        !print *,'tsgrid C after csub ',C(1,2,2,1)

! ***    The routine below facilitate the inclusion of dry- and
! ***    wet-deposition as a separate process instead of being
! ***    introduced in the vertical diffusion routines or in
! ***    the chemistry scheme:

        call gmdepo

        !print *,'gmphot: after gmdepo'
        !print *,'tsgrid C after gmdepo ',C(1,2,2,1)

!       Calculate grid model radioactive decay:

        call gmradi

        !print *,'tsgrid: after gmradi'
        !print *,'tsgrid C after gmradi ',C(1,2,2,1)

! ***   Calculate grid model photochemical reactions:

!_EMEP03	  IF (MESSFE) WRITE (MESSUN,*) 'TSGRID: Next action CALL GMPHOT.'

         call gmphot

        !print *,'tsgrid: after gmphot',C(1,2,2,1)

! ***   The SOA-module is applied if NC >= 68:

!!DEC$ IF DEFINED (emep68)
!
!!_TEST        IF (NC >= 68 ) THEN
!
!! ***   SOA-MODULE is to be called here. The Chemistry calculations
!! ***   have been made in SUBROUTINE GMPHOT and immediately after
!! ***   the partitioning calculations should be performed.
!
!          DO IX = 1,NX
!            DO IY = 1,NY
!
!! ***         For the present TEST setup the variable setting below
!! ***         could have been moved outside of the IX,IY - DO-LOOP.
!! ***         However, for real applications these variables will be
!! ***         updated from their 3-D model variables.
!
!              temp_soa(1:nz)  = tair_soa(IX,IY,1:nz)      ! 298.0
!              itemp_soa(1:nz) = INT(temp_soa(1:nz)+0.5)
!
!!_SOA        [xn] is concentrations in: molecules/cm3.
!              xn(FIRST_NONVOL:LAST_SOA,1:nz) =
!     &           C(FIRST_NONVOL:LAST_SOA,IX,IY,1:nz)   ! 0.1 * ppb ?
!
!!             Update the total Organic_PM_mass_conc (ug/m3) after
!!             having performed the advection and diffusion processes:
!!             If the line below is commented out we do not transport the
!!             total organic PM mass concentration.
!
!              Grid_M0(IX,IY,1:nz) = C(ic_org_PM2p5,IX,IY,1:nz)
!
!              call OrganicAerosol(DEBUG,IX,IY)
!
!!              IF (MESSFE) WRITE (MESSUN,'(A46,I3,A,I3,A)')
!!     &          ' OrganicAerosol has been completed for I,J = ('
!!     &         ,IX,',',IY,')'
!
!              C(ic_org_PM2p5,IX,IY,1:nz)  = Grid_M0(IX,IY,1:nz)
!              C(ic_aerOM,IX,IY,1:nz)      = aer_ug_OC(1:nz)
!              C(ic_aerPOM,IX,IY,1:nz)     = aer_ug_POM(1:nz)
!              C(ic_aerWOOD,IX,IY,1:nz)    = aer_ug_WOOD(1:nz)
!              C(ic_aerFFUEL,IX,IY,1:nz)   = aer_ug_FFUEL(1:nz)
!              C(ic_aerBGND,IX,IY,1:nz)    = aer_ug_BGNDOC(1:nz)
!              C(ic_aerASOA,IX,IY,1:nz)    = aer_ug_ASOA(1:nz)
!              C(ic_aerBSOA,IX,IY,1:nz)    = aer_ug_BSOA(1:nz)
!
!              C(ic_aerOM_C,IX,IY,1:nz)    = aer_ugC_OC(1:nz)
!              C(ic_aerPOM_C,IX,IY,1:nz)   = aer_ugC_POM(1:nz)
!              C(ic_aerWOOD_C,IX,IY,1:nz)  = aer_ugC_WOOD(1:nz)
!              C(ic_aerFFUEL_C,IX,IY,1:nz) = aer_ugC_FFUEL(1:nz)
!              C(ic_aerBGND_C,IX,IY,1:nz)  = aer_ugC_BGNDOC(1:nz)
!              C(ic_aerASOA_C,IX,IY,1:nz)  = aer_ugC_ASOA(1:nz)
!              C(ic_aerBSOA_C,IX,IY,1:nz)  = aer_ugC_BSOA(1:nz)
!
!            END DO
!          END DO
!
!!_TEST        ENDIF !(NC >= 68)
!
!! ***   SOA-MODULE is ended here.
!
!!DEC$ ENDIF

!_LHS_SOA_June_2007_End.

      ELSE

! ***   Calculate grid model photochemical reactions

!_LHS_SOA_May_2007_Start:
!
! ***   Commented out since the timestep is set to 2 * DT within
! ***   SUBROUTINE GMPHOT at the end of the previous timestep.
!
! ***       CALL GMPHOT
!
!_LHS_SOA_May_2007_End.
!
!
!       Calculate grid model radioactive decay

        call gmradi

!debug      print *,'tsgrid C after gmradi alt  ',C(1,2,2,1)

!MSK *** This is the connection between Eulerian grid concentration C
!MSK *** and the local main grid, receptor point, line receptor
!MSK *** and subgrid concentrations CM, CR, CLR, CS

!       Calculate grid model output concentration and deposition

! ***   ... in the main grid cells

        call csubgm(AV_GRIDAI,GRIDAI,GRIDAI)

! ***   ... in the sub  grid cells

!MSK        CALL CSUBGS(AV_GRIDAI,GRIDAI,GRIDAI)

! ***   ... in the irregular receptor points

        call csubgr(AV_GRIDAI,GRIDAI,GRIDAI)

! ***   ... in the line source associated receptor points

        call csubgl(AV_GRIDAI,GRIDAI,GRIDAI)

!_LHS_Hourly_Averaged_Output_December_2007_End.

!       ... in the grid model

!       PRINT *,'GMDEPO ',ITS

      !print *,'tsgrid C after csub alt  ',C(1,2,2,1)

        call gmdepo

      !print *,'tsgrid C after gmdepo alt  ',C(1,2,2,1)

!       Calculate vertical advection and diffusion

        DO IC = 1,NC - n_nochem + n_advect

          if (messfe .AND. (ITS == 1 .OR. ITS == NTS) ) then
            write(messun,'(A,I4)') 'TSGRID: Preprocessing timestep ITS = ',ITS
              rsum_DCDT     = SUM(DCDT(IC,:,:,:))
            write(messun,'(A,I2,A,E12.5)')  &
              'TSGRID: Domain sum DCDT(IC=',IC,',:,:,:) = ',rsum_DCDT
          endif

!MSK start   avoid tiny concentration values
          do IZ=1,NZ+NBCZ
            do IY=1, NY+2*NBCY
              do IX=1, NX+2*NBCX
                if ( C(IC,IX,IY,IZ) .lt. 1.e-20)  C(IC,IX,IY,IZ) = 0.0
              enddo
            enddo
          enddo
!MSK end

          IF (GRIDVS(IC) .EQ. 0) THEN

!           Introduce emissions directly into the grid when
!           there is no vertical operator

            DO IZ = 1,NZ
              DO IY = 1,NY
                DO IX = 1,NX
!MSK                  C(IC,IX,IY,IZ) = C(IC,IX,IY,IZ) + DCDT(IC,IX,IY,IZ)*DT
                  C(IC,IX+1,IY+1,IZ) = C(IC,IX+1,IY+1,IZ) +    &
                                       DCDT(IC,IX,IY,IZ)*DT
                END DO
              END DO
            END DO

          ELSEIF (GRIDVS(IC) .EQ. 1) THEN
              call adiv(ic)

          ELSEIF (GRIDVS(IC) .EQ. 2) THEN
              call v_dif(ic)
              call v_adv(ic)

          ELSEIF (GRIDVS(IC) .EQ. 3) THEN

              call difv_cn(ic)
              call v_adv(ic)

          ELSEIF (GRIDVS(IC) .EQ. 4) THEN
              call difv_cn(ic)

          ELSEIF (GRIDVS(IC) .EQ. 5) THEN
! ***       This should only be applied for test purposes
! ***       when no grid emissions are included.
              call v_adv(ic)
          ENDIF

        END DO

      !print *,'tsgrid C after vertical alt  ',C(1,2,2,1)

!_LHS_SOA_May_2007_End.

!       Calculate horisontal 2D diffusion

!       PRINT *,'DIVH ',ITS
          call difh

      !print *,'tsgrid C after hor dif alt  ',C(1,2,2,1)

!       Calculate horisontal 2D advection

!       PRINT *,'*ADVH ',ITS,MOD(ITS,2)
          call advh
!       PRINT *,'*END ',ITS

      !print *,'tsgrid C after hor adv alt  ',C(1,2,2,1)

      ENDIF

!     If NOT NTS-Averaging, calculate grid model concentrations.

      IF (ITS == NTS .AND. AV_GRIDAI == 0) THEN

!       ... in the main grid cells

        call csubgm(GRIDAI,0,0)

!MSK debug
!MSK debug        print *,'tsgrid CR(NO) after csubgm alt', CR(2,10000)

!       ... in the sub  grid cells

!MSK        CALL CSUBGS(GRIDAI,0,0)

!       ... in the irregular receptor points

        call csubgr(GRIDAI,0,0)

!MSK debug
!MSK debug        print *,'tsgrid CR(NO) after csubgr alt', CR(2,10000)

!       ... in the line source associated receptor points

        call csubgl(GRIDAI,0,0)

!debug      print *,'tsgrid C after csub alt  ',C(1,2,2,1)

      ENDIF

!     Finished with preprocessing

          !print *,'tsgrid: finished preprocessing'

      RETURN

  200 CONTINUE

!     Postprocessing calculations:
!     ****************************

!     If run model not ok then return

      IF (.NOT. RUNOK) RETURN

!_NEST_SIGMA_Start: *********  START FLUX-CALCULATIONS   *****************
!
      IF (INFOFE) THEN
!       Calculate total volume, total mass and volume averaged conc.:
!!MSK             CALL V_ADV(IC)  !buggy routine
!       The variables TOTVOL(MC), TOTMASS(MC) and VOLAVG(MC) are declared
!       in: /grid/mod_grid.for

        IF (ITS .EQ. NTS) THEN

          DO IC = 1,NC

            TOTVOL(IC)  = 0.0
            TOTMASS(IC) = 0.0
            VOLAVG(IC)  = 0.0

            WESTB_VF(IC) = 0.0
            WESTB_MF(IC) = 0.0

            EASTB_VF(IC) = 0.0
            EASTB_MF(IC) = 0.0

            NORTHB_VF(IC) = 0.0
            NORTHB_MF(IC) = 0.0

            SOUTHB_VF(IC) = 0.0
            SOUTHB_MF(IC) = 0.0

            TOPB_VF(IC) = 0.0
            TOPB_MF(IC) = 0.0

            BOTB_VF(IC) = 0.0
            BOTB_MF(IC) = 0.0

!           When multiplying C(IC,IX,IY,IZ) with molec2ug we
!           converts from [molecules/cm3] to [ug/m3].

            IF (CUNIT(IC)(1:3) .EQ. 'mol')THEN
!            IF(IC <= NC - n_nochem)THEN
                molec2ug = CMOLW(IC)*1.0E+12/avogad
            ELSE
                molec2ug = 1.0
            ENDIF

!MSK *** C in the FLUX calculation has been shifted by 1 in x- and y-direction

!	      PRINT *,' molec2ug = ',molec2ug
!           Calculates the total model volum and the total model mass:
            DO IZ = 1,NZ
              DO IY = 1,NY
               DO IX = 1,NX
!_Fine_Model domain	    DO IY = 8,16
!_Fine_Model domain          DO IX = 6,16

              TOTVOL(IC)  = TOTVOL(IC) + (VOL(IZ)*DEPTHM(IX,IY)/MOD_H)

              TOTMASS(IC) = TOTMASS(IC) +  &
               (C(IC,IX+1,IY+1,IZ)*molec2ug*VOL(IZ)*DEPTHM(IX,IY)/MOD_H )

                ENDDO
               ENDDO
              ENDDO

!           Calculates the Volume fluxes and the Mass fluxes out of the
!           model boundaries:

!           The EAST/WEST-Boundaries:
               DO IZ=1,NZ
                  DO IY=1,NY
!_Fine_Model domain		  DO IY=8,16

               WESTB_VF(IC) = WESTB_VF(IC) +   &
                             U(1,IY,IZ)*DY*DZ(IZ)*DEPTHM(1,IY)/MOD_H
               WESTB_MF(IC) = WESTB_MF(IC) +    &
                             C(IC,2,IY+1,IZ)*molec2ug*  &
                             U(1,IY,IZ)*DY*DZ(IZ)*DEPTHM(1,IY)/MOD_H

!_NEST_NOTE: The sign for the turbulent flux is chosen so that this
!            transport is positive if the transport is directed in the
!            positive x-direction (eastwards):
!             WESTB_MF(IC) = WESTB_MF(IC) +
!     &          ( (D(IZ)*molec2ug*(BC(IC,0,IY,IZ) - C(IC,1,IY,IZ))/DX)*
!     &                      (DY*DZ(IZ)*DEPTHM(1,IY)/MOD_H) )

               EASTB_VF(IC) = EASTB_VF(IC) +   &
                             U(NX,IY,IZ)*DY*DZ(IZ)*DEPTHM(NX,IY)/MOD_H
               EASTB_MF(IC) = EASTB_MF(IC) +    &
                             C(IC,NX+1,IY+1,IZ)*molec2ug*   &
                             U(NX,IY,IZ)*DY*DZ(IZ)*DEPTHM(NX,IY)/MOD_H

!_NEST_NOTE: The sign for the turbulent flux is chosen so that this
!            transport is positive if the transport is directed in the
!            positive x-direction (eastwards):
!             EASTB_MF(IC) = EASTB_MF(IC) -
!     &        ( (D(IZ)*molec2ug*(BC(IC,NX+1,IY,IZ) - C(IC,NX,IY,IZ))/DX)*
!     &                   (DY*DZ(IZ)*DEPTHM(NX,IY)/MOD_H) )
              ENDDO
             ENDDO

!            The NORTH/SOUTH-Boundaries:
             DO IZ=1,NZ
               DO IX=1,NX
!_Fine_Model domain		  DO IX=6,16
               NORTHB_VF(IC) = NORTHB_VF(IC) +    &
                              V(IX,NY,IZ)*DY*DZ(IZ)*DEPTHM(IX,NY)/MOD_H
               NORTHB_MF(IC) = NORTHB_MF(IC) +    &
                              C(IC,IX+1,NY+1,IZ)*molec2ug*   &
                              V(IX,NY,IZ)*DY*DZ(IZ)*DEPTHM(IX,NY)/MOD_H

!_NEST_NOTE: The sign for the turbulent flux is chosen so that this
!            transport is positive if the transport is directed in the
!            positive y-direction (northwards):
!             NORTHB_MF(IC) = NORTHB_MF(IC) -
!     &        ( (D(IZ)*molec2ug*(BC(IC,IX,NY+1,IZ) - C(IC,IX,NY,IZ))/DY)*
!     &                (DY*DZ(IZ)*DEPTHM(IX,NY)/MOD_H) )

               SOUTHB_VF(IC) = SOUTHB_VF(IC) +    &
                              V(IX,1,IZ)*DY*DZ(IZ)*DEPTHM(IX,1)/MOD_H
               SOUTHB_MF(IC) = SOUTHB_MF(IC) +    &
                              C(IC,IX+1,2,IZ)*molec2ug*   &
                              V(IX,1,IZ)*DY*DZ(IZ)*DEPTHM(IX,1)/MOD_H

!_NEST_NOTE: The sign for the turbulent flux is chosen so that this
!            transport is positive if the transport is directed in the
!            positive y-direction (northwards):
!             SOUTHB_MF(IC) = SOUTHB_MF(IC) +
!     &        ( (D(IZ)*molec2ug*(BC(IC,IX,0,IZ) - C(IC,IX,1,IZ))/DY) *
!     &                       (DY*DZ(IZ)*DEPTHM(IX,1)/MOD_H) )
               ENDDO
             ENDDO

!           The TOP/(Bottom)-Boundaries:
            DO IX=1,NX
              DO IY=1,NY
!_Fine_Model domain		DO IX=6,16
!_Fine_Model domain		  DO IY=8,16

!_SIGMA_NOTE: W(IX,IY,NZ) must be multiplied by DEPTHM(IX,IY)/MOD_H
!             in order to get the cartesian vertical velocity. None
!             of the other terms are necessary since at the top
!             boundary sigma is equal to H_0 and the additional terms
!             therefore vanish.
               TOPB_VF(IC) = TOPB_VF(IC) +    &
                            (W(IX,IY,NZ)*DX*DY*DEPTHM(IX,IY)/MOD_H)
               TOPB_MF(IC) = TOPB_MF(IC) +    &
                            C(IC,IX+1,IY+1,NZ)*molec2ug*   &
                            (W(IX,IY,NZ)*DX*DY*DEPTHM(IX,IY)/MOD_H)

!_SIGMA_NOTE:  In /mete/cdzdt.for DZDT(I,J,NZ) is multiplied with
!              square of the stretch factor, so we do not need to
!              include this here. The sign below is chosen so that
!              upward (outward) flux is positive.
! 	       TOPB_MF(IC) = TOPB_MF(IC) +  DZDT(IX,IY,NZ) *
!     &         ( molec2ug*(C(IC,IX,IY,NZ) - BC(IC,IX,IY,NZ+1)) *
!     &                     (DX*DY) )
!             BOTB_VF(IC) = BOTB_VF(IC) + W(IX,IY,1)*DX*DY
!             BOTB_MF(IC) = BOTB_MF(IC) + C(IC,IX,IY,1)*BOTB_VF(IC)

              ENDDO
             ENDDO

!            Horizontal, vertical and total Volume and Mass Flux
!            INTO the model domain, i.e. THE INFLUX:

             HOR_VFLUX  = WESTB_VF(IC)  - EASTB_VF(IC)    &
                        + SOUTHB_VF(IC) - NORTHB_VF(IC)

             HOR_MFLUX  = WESTB_MF(IC)  - EASTB_MF(IC)    &
                        + SOUTHB_MF(IC) - NORTHB_MF(IC)

!	       VERT_VFLUX = BOTB_VF(IC) - TOPB_VF(IC)
!	       VERT_MFLUX = BOTB_MF(IC) - TOPB_MF(IC)

             VERT_VFLUX = - TOPB_VF(IC)
             VERT_MFLUX = - TOPB_MF(IC)

!           TOTAL Volume and Mass transport INTO the model domain
!           i.e. The total INFLUX:

             TOT_VFLUX = HOR_VFLUX + VERT_VFLUX
             TOT_MFLUX = HOR_MFLUX + VERT_MFLUX

!           Giving the quantities as (km3/s) and as (kg/s) if
!           both are multipled by 1.0E-09:
             EASTB_VF(IC) = EASTB_VF(IC) * 1.0E-09
!            EASTB_MF(IC) = EASTB_MF(IC) * 1.0E-06
             EASTB_MF(IC) = EASTB_MF(IC) * 1.0E-09
             HOR_VFLUX = HOR_VFLUX  * 1.0E-09
!            HOR_MFLUX = HOR_MFLUX  * 1.0E-06
             HOR_MFLUX = HOR_MFLUX  * 1.0E-09
             VERT_VFLUX = VERT_VFLUX * 1.0E-09
!            VERT_MFLUX = VERT_MFLUX * 1.0E-06
             VERT_MFLUX = VERT_MFLUX * 1.0E-09
             TOT_VFLUX = TOT_VFLUX * 1.0E-09
!            TOT_MFLUX = TOT_MFLUX * 1.0E-06
             TOT_MFLUX = TOT_MFLUX * 1.0E-09

             VOLAVG(IC) = TOTMASS(IC)/TOTVOL(IC)

!_Include the output below to the log file if necessary:
!	        IF (MESSFE) WRITE (MESSUN,2002) IC,TOTMASS(IC)
!	        IF (MESSFE) WRITE (MESSUN,2003) IC,VOLAVG(IC)
!	        IF (MESSFE) WRITE (MESSUN,2004) IC,TOTVOL(IC)

!             Convert Volume from m**3 to Km**3:
             TOTVOL(IC) = TOTVOL(IC) * 1.0E-9
!             Convert total NOx Mass from ug to Kg:
             TOTMASS(IC) = TOTMASS(IC) * 1.0E-9

             WRITE (INFOUN,2005) CMPND(IC),IC,THOUR,   &
                    TOTVOL(IC),TOTMASS(IC),VOLAVG(IC),        &
                    EASTB_MF(IC),HOR_MFLUX,VERT_MFLUX,TOT_MFLUX
!     &              EASTB_VF(IC),HOR_VFLUX,VERT_VFLUX,TOT_VFLUX

!
!			IF (MESSFE) WRITE (MESSUN,2067) D(1)
!			IF (MESSFE) WRITE (MESSUN,2068) D(2)
!			IF (MESSFE) WRITE (MESSUN,2069) D(3)

           ENDDO   ! DO IC = 1,NC

         ENDIF  ! IF (ITS .EQ. NTS)

      ENDIF  !  IF (INFOFE)

!_NEST_SIGMA_End *************  END FLUX-CALCULATIONS   *******************************************

      IF (ISH == NSH .AND. ITS == NTS) THEN

! ***  Calculate the concentrations for the NEST grid domain at the end of each hour.
! ***  These values are calculated for the entire NEST grid, and are stored in the array:
!
! ***       BC_nest(NC,1-nbcx_nest:nx_nest+nbcx_nest,
! ***      &           1-nbcy_nest:ny_nest+nbcy_nest,
! ***      &                       nz_nest+nbcz_nest))

!MSK start  New subgrid BC values, only wanted if there are NEWBC files
!MSK     Now commented subrid BC and IC files
!MSK        if (NEWBCFE( 1 )) then
!MSK
!MSK          CALL WNEWBC
!MSK
!MSK        endif
!MSK end

! ***   The export of the "BC_nest" data is moved to the "ReceiveTimeDependentOutputData"
! ***   routine which is called after the present call of "CallRunEpi(picallcnt)"

      ENDIF

!MSK start
! *** Writing out RESTART concentration values:
      if (ATDATE(EDAT) .AND. ISH == NSH .AND. ITS == NTS) then

!MSK        CALL WNEWC
!!==========================================================================
!! write out main grid emission data at last timestep of the simulation to .dat file
!! added by Kang @ CARES, Dec. 23, 2019
!!==========================================================================
        if (.not. allocated(field3D_output) )  allocate(field3D_output(NX,NY,NZ,NC))
!!==========================================================================

! Write the netCDF output file of main grid C() for RESTART
          if (IC1GRIDFE) then

              ! one time step
              Nhh_in = 1

              ! current simulation date
              mdate(1,Nhh_in) = YEAR
              mdate(2,Nhh_in) = MNTH
              mdate(3,Nhh_in) = DAYM
              mdate(4,Nhh_in) = HOUR

              unitname = "ug/m^3"

              if (averaged_output) then
                 validity = 'averaged'
              else
                 validity = 'instantaneous'
              endif

              dopacking = .false.
              dofloat   = .true.
              domirror  = .false.

              do 300 IC = 1,NC

                namefield = CMPND(IC)
                if ( trim(CMPND(IC)).eq."PM2.5") then
                  namefield = "PM25      "
                endif

                ! Do neglect the outer boundaries
                ! Read C() from IX=2, IY=2, until NZ
                do 310 IZ=1,NZ
                  do 320 IX = 1,NX
                    do 330 IY = 1,NY

                       field3D(IX,IY,IZ) = C(IC,IX+1,IY+1,IZ)
                  !     field3D(IY,IX,IZ) = C(IC,IX+1,IY+1,IZ)
!!==========================================================================
!! write out main grid emission data at last timestep of the simulation to .dat file
!! added by Kang @ CARES, Dec. 23, 2019
!!==========================================================================
                 field3D_output(IX,IY,IZ,IC) = C(IC,IX+1,IY+1,IZ)
!!==========================================================================
  330               continue
  320             continue
  310           continue

                call writeconcfield(IC1GRIDFN, namefield,unitname,field3D(1:NX,1:NY,1:NZ), NX, NY, NZ, &
                                 Nhh_in, mdate, validity, dopacking, dofloat, domirror)

!!!!!!!!!!!!!!!!!!! for testing 3D concentration output file by Kang @ Dec.2, 2019
!               if (IC .eq. 3) then
!                open(unit=1111,file="../OUTPUT/3D_field_write.txt")  !!
!                open(unit=1112,file="../OUTPUT/3D_field_para_write.txt")  !!
!                do 311 IZ=1,NZ
!                    do 331 IY = 1,NY

!                      write(1111,*) (field3D(IX,IY,IZ), IX=1,NX)
                    !  write(1111,*) (C(3,IX+1,IY+1,IZ), IX=1,NX)
!  331               continue

!  311           continue

!                write(1112,*)"write out 3D concentration file, namefield=",namefield,"unitname=",unitname, &
!                              ", NX=",NX,", NY=",NY, "NZ=",NZ, "Nhh_in=",Nhh_in, "mdate=",mdate,  &
!                              ", validity=",validity,", dopacking=",dopacking,", dofloat=", dofloat,  &
!                              ", domirror=",domirror, ", IC1GRIDFN=", IC1GRIDFN

!                close(1111)
!                close(1112)
!                endif
!!!!!!!!!!!!!!!!!!! end for testing 3D concentration output file by Kang @ Dec.2, 2019

  300          continue

          endif
!MSK end
!!===================================================================================
!! write out main grid emission data at last timestep of the simulation to .dat file
!! added by Kang @ CARES, Dec. 23, 2019
!!===================================================================================

      open(unit=1116,file="../OUTPUT/3D_instantanous_mainconc_center.dat")  !! 3D data at center point of main cell

      write(1116,2140)(CMPND(IC),IC=1,NC)

       do  IZ=1,NZ
            do  IX = 1,NX
                do  IY = 1,NY
                    mxm_i(IY,IX) = SITEX0 + (IX-1)*DX + 0.5*DX
                    mxm_j(IY,IX) = SITEY0 + (IY-1)*DY + 0.5*DY
             write(1116,2150) (mdate(II,Nhh_in),II=1,4), mxm_i(IY,IX),mxm_j(IY,IX),Z(IZ)-0.5*DZ(IZ),&
                      (field3D_output(IX,IY,IZ,IC), IC=1,NC)
                enddo
          enddo
      enddo

      close(1116)

!$$$$$$
!$$$$$$       open(unit=1113,file="../OUTPUT/3D_instantanous_mainconc.dat")  !! 3D data at SW+top point of main cell
!$$$$$$
!$$$$$$       write(1113,2140)(CMPND(IC),IC=1,NC)
!$$$$$$
!$$$$$$        do  IZ=1,NZ
!$$$$$$             do  IX = 1,NX
!$$$$$$                 do  IY = 1,NY
!$$$$$$                     mxm_i(IY,IX) = SITEX0 + (IX-1)*DX
!$$$$$$                     mxm_j(IY,IX) = SITEY0 + (IY-1)*DY
!$$$$$$              write(1113,2150) (mdate(II,Nhh_in),II=1,4), mxm_i(IY,IX),mxm_j(IY,IX),Z(IZ),&
!$$$$$$                                Z(IZ)*DEPTHM(IX,IY)/MOD_H,(field3D_output(IX,IY,IZ,IC), IC=1,NC)
!$$$$$$                 enddo
!$$$$$$           enddo
!$$$$$$       enddo

        if (allocated(field3D_output))  deallocate(field3D_output)

! Deallocate
        if (allocated(mxm_i))     deallocate(mxm_i)
        if (allocated(mxm_j))     deallocate(mxm_j)

!$$$$$$        close(1113)

   !   open(unit=1114,file="../OUTPUT/3D_instantanous_mainconc_centerBCZ.dat") !!!! 3D data at center point of main cell including (BC @ x,y,z direction)

   !   write(1114,2140)(CMPND(IC),IC=1,NC)

   !    do  IZ=1,NZ+NBCZ
   !             if (IZ .lt. NZ+NBCZ) then
   !             Znew(IZ)=Z(IZ)-0.5*DZ(IZ)
   !             else
   !             Znew(IZ)=Znew(IZ-1)+DZ(IZ-1)
   !             endif

   !         do  IX = 1,NX+2*NBCX
   !             do  IY = 1,NY+2*NBCY
   !                 mxm_i_new(IY,IX) = SITEX0 + (IX-1)*DX-0.5*DX
   !                 mxm_j_new(IY,IX) = SITEY0 + (IY-1)*DY-0.5*DY

   !          if (IX .ne. 1 .and. IX .ne. NX+2*NBCX .and. IY .ne. 1 .and. IY .ne. NY+2*NBCY ) then

   !          write(1114,2150) (mdate(II,Nhh_in),II=1,4), mxm_i_new(IY,IX),mxm_j_new(IY,IX),Znew(IZ),&
!$$$$$$                                Znew(IZ)*DEPTHM(IX-1,IY-1)/MOD_H,(C(IC,IX,IY,IZ), IC=1,NC)
   !          endif

   !          if (IX .eq. 1) then
   !            if (IY .eq. NY+2*NBCY) then
   !          write(1114,2150) (mdate(II,Nhh_in),II=1,4), mxm_i_new(IY,IX),mxm_j_new(IY,IX),Znew(IZ),&
   !                            Znew(IZ)*DEPTHM(IX,IY-2)/MOD_H,(C(IC,IX,IY,IZ), IC=1,NC)
   !            endif
   !            if (IY .eq. 1) then
   !         write(1114,2150) (mdate(II,Nhh_in),II=1,4), mxm_i_new(IY,IX),mxm_j_new(IY,IX),Znew(IZ),&
   !                            Znew(IZ)*DEPTHM(IX,IY)/MOD_H,(C(IC,IX,IY,IZ), IC=1,NC)
   !            endif
   !            if (IY .ne.1 .and. IY .ne. NY+2*NBCY) then
   !         write(1114,2150) (mdate(II,Nhh_in),II=1,4), mxm_i_new(IY,IX),mxm_j_new(IY,IX),Znew(IZ),&
   !                            Znew(IZ)*DEPTHM(IX,IY-1)/MOD_H,(C(IC,IX,IY,IZ), IC=1,NC)
   !            endif
   !          endif
   !          if (IX .eq. NX+2*NBCX) then
!$$$$$$                if (IY .eq. NY+2*NBCY) then
!$$$$$$              write(1114,2150) (mdate(II,Nhh_in),II=1,4), mxm_i_new(IY,IX),mxm_j_new(IY,IX),Znew(IZ),&
!$$$$$$                                Znew(IZ)*DEPTHM(IX-2,IY-2)/MOD_H,(C(IC,IX,IY,IZ), IC=1,NC)
!$$$$$$                endif
!$$$$$$               if (IY .eq. 1) then
!$$$$$$             write(1114,2150) (mdate(II,Nhh_in),II=1,4), mxm_i_new(IY,IX),mxm_j_new(IY,IX),Znew(IZ),&
!$$$$$$                                Znew(IZ)*DEPTHM(IX-2,IY)/MOD_H,(C(IC,IX,IY,IZ), IC=1,NC)
!$$$$$$                endif
!$$$$$$               if (IY .ne.1 .and. IY .ne. NY+2*NBCY) then
!$$$$$$             write(1114,2150) (mdate(II,Nhh_in),II=1,4), mxm_i_new(IY,IX),mxm_j_new(IY,IX),Znew(IZ),&
!$$$$$$                                Znew(IZ)*DEPTHM(IX-2,IY-1)/MOD_H,(C(IC,IX,IY,IZ), IC=1,NC)
!$$$$$$                endif
!$$$$$$              endif
!$$$$$$
!$$$$$$              if ( IX .ne. 1 .and. IX .ne. NX+2*NBCX) then
!$$$$$$                if (IY .eq. 1) then
!$$$$$$              write(1114,2150) (mdate(II,Nhh_in),II=1,4), mxm_i_new(IY,IX),mxm_j_new(IY,IX),Znew(IZ),&
!$$$$$$                                Znew(IZ)*DEPTHM(IX-1,IY)/MOD_H,(C(IC,IX,IY,IZ), IC=1,NC)
!$$$$$$                endif
!$$$$$$                if (IY .eq.  NY+2*NBCY) then
!$$$$$$              write(1114,2150) (mdate(II,Nhh_in),II=1,4), mxm_i_new(IY,IX),mxm_j_new(IY,IX),Znew(IZ),&
!$$$$$$                                Znew(IZ)*DEPTHM(IX-1,IY-2)/MOD_H,(C(IC,IX,IY,IZ), IC=1,NC)
!$$$$$$                endif
!$$$$$$              endif
!$$$$$$
!$$$$$$                 enddo
!$$$$$$           enddo
!$$$$$$       enddo
!$$$$$$
! Deallocate
!$$$$$$         if (allocated(Znew))     deallocate(Znew)
!$$$$$$        close(1114)

!$$$$$$       open(unit=1115,file="../OUTPUT/3D_instantanous_mainconc_centerBC.dat") !!!! 3D data at center point of main cell including (BC @ x,y direction)
!$$$$$$
!$$$$$$       write(1115,2140)(CMPND(IC),IC=1,NC)
!$$$$$$
!$$$$$$        do  IZ=1,NZ
!$$$$$$
!$$$$$$                 Znew0(IZ)=Z(IZ)-0.5*DZ(IZ)
!$$$$$$
!$$$$$$             do  IX = 1,NX+2*NBCX
!$$$$$$                 do  IY = 1,NY+2*NBCY
!$$$$$$                     mxm_i_new(IY,IX) = SITEX0 + (IX-1)*DX-0.5*DX
!$$$$$$                     mxm_j_new(IY,IX) = SITEY0 + (IY-1)*DY-0.5*DY
!$$$$$$
!$$$$$$              if (IX .ne. 1 .and. IX .ne. NX+2*NBCX .and. IY .ne. 1 .and. IY .ne. NY+2*NBCY ) then
!$$$$$$              write(1115,2150) (mdate(II,Nhh_in),II=1,4), mxm_i_new(IY,IX),mxm_j_new(IY,IX),Znew0(IZ),&
!$$$$$$                                Znew0(IZ)*DEPTHM(IX-1,IY-1)/MOD_H,(C(IC,IX,IY,IZ), IC=1,NC)
!$$$$$$              endif
!$$$$$$
!$$$$$$              if (IX .eq. 1) then
!$$$$$$                if (IY .eq. NY+2*NBCY) then
!$$$$$$              write(1115,2150) (mdate(II,Nhh_in),II=1,4), mxm_i_new(IY,IX),mxm_j_new(IY,IX),Znew0(IZ),&
!$$$$$$                                Znew0(IZ)*DEPTHM(IX,IY-2)/MOD_H,(C(IC,IX,IY,IZ), IC=1,NC)
!$$$$$$                endif
!$$$$$$                if (IY .eq. 1) then
!$$$$$$             write(1115,2150) (mdate(II,Nhh_in),II=1,4), mxm_i_new(IY,IX),mxm_j_new(IY,IX),Znew0(IZ),&
!$$$$$$                                Znew0(IZ)*DEPTHM(IX,IY)/MOD_H,(C(IC,IX,IY,IZ), IC=1,NC)
!$$$$$$                endif
!$$$$$$                if (IY .ne.1 .and. IY .ne. NY+2*NBCY) then
!$$$$$$             write(1115,2150) (mdate(II,Nhh_in),II=1,4), mxm_i_new(IY,IX),mxm_j_new(IY,IX),Znew0(IZ),&
!$$$$$$                                Znew0(IZ)*DEPTHM(IX,IY-1)/MOD_H,(C(IC,IX,IY,IZ), IC=1,NC)
!$$$$$$                endif
!$$$$$$              endif
!$$$$$$
!$$$$$$              if (IX .eq. NX+2*NBCX) then
!$$$$$$                if (IY .eq. NY+2*NBCY) then
!$$$$$$              write(1115,2150) (mdate(II,Nhh_in),II=1,4), mxm_i_new(IY,IX),mxm_j_new(IY,IX),Znew0(IZ),&
!$$$$$$                                Znew0(IZ)*DEPTHM(IX-2,IY-2)/MOD_H,(C(IC,IX,IY,IZ), IC=1,NC)
!$$$$$$                endif
!$$$$$$               if (IY .eq. 1) then
!$$$$$$             write(1115,2150) (mdate(II,Nhh_in),II=1,4), mxm_i_new(IY,IX),mxm_j_new(IY,IX),Znew0(IZ),&
!$$$$$$                                Znew0(IZ)*DEPTHM(IX-2,IY)/MOD_H,(C(IC,IX,IY,IZ), IC=1,NC)
!$$$$$$                endif
!$$$$$$               if (IY .ne.1 .and. IY .ne. NY+2*NBCY) then
!$$$$$$             write(1115,2150) (mdate(II,Nhh_in),II=1,4), mxm_i_new(IY,IX),mxm_j_new(IY,IX),Znew0(IZ),&
!$$$$$$                                Znew0(IZ)*DEPTHM(IX-2,IY-1)/MOD_H,(C(IC,IX,IY,IZ), IC=1,NC)
!$$$$$$                endif
!$$$$$$              endif
!$$$$$$
!$$$$$$              if ( IX .ne. 1 .and. IX .ne. NX+2*NBCX) then
!$$$$$$                if (IY .eq. 1) then
!$$$$$$              write(1115,2150) (mdate(II,Nhh_in),II=1,4), mxm_i_new(IY,IX),mxm_j_new(IY,IX),Znew0(IZ),&
!$$$$$$                                Znew0(IZ)*DEPTHM(IX-1,IY)/MOD_H,(C(IC,IX,IY,IZ), IC=1,NC)
!$$$$$$                endif
!$$$$$$               if (IY .eq.  NY+2*NBCY) then
!$$$$$$              write(1115,2150) (mdate(II,Nhh_in),II=1,4), mxm_i_new(IY,IX),mxm_j_new(IY,IX),Znew0(IZ),&
!$$$$$$                                Znew0(IZ)*DEPTHM(IX-1,IY-2)/MOD_H,(C(IC,IX,IY,IZ), IC=1,NC)
!$$$$$$                endif
!$$$$$$              endif
!$$$$$$
!$$$$$$                 enddo
!$$$$$$           enddo
!$$$$$$       enddo

! Deallocate
!$$$$$$         if (allocated(mxm_i_new))     deallocate(mxm_i_new)
!$$$$$$         if (allocated(mxm_j_new))     deallocate(mxm_j_new)
!$$$$$$         if (allocated(Znew0))     deallocate(Znew0)
!$$$$$$        close(1115)

!!==========================================================================
      endif

!!DEC$ IF DEFINED (emep68)
!      if (ATDATE(EDAT) .AND. ISH == NSH .AND. ITS == NTS) then
!!_TEST	  IF (NC >= 68 ) THEN
!
!! ***     Deallocate the allocatable arrays that were declared as
!! ***     private in OrganicAerosol:
!          my_last_call = .true.
!	    call OrganicAerosol(DEBUG,1,1)
!
!!_TEST        ENDIF
!      endif
!!DEC$ ENDIF (emep68)

!     Finished with postprocessing

      RETURN

 2000 FORMAT ('TSGRID: IC,BC(IC,15,15,1) = ',I3,E16.7)
 2002 FORMAT ('TSGRID: IC,Total Model Mass    = ',I3,E16.7)
 2003 FORMAT ('TSGRID: IC,Volme Average Conc. = ',I3,E16.7)
 2004 FORMAT ('TSGRID: IC,Total Model Volume  = ',I3,E16.7)
 2005 FORMAT ('TSGRID:',2X,A12,2X,I4,1X,F7.0,1X,F10.1,1X,E12.4,1X,  &
                       E12.5,1X,E17.6,1X,E17.6,1X,E17.6,1X,E17.6)

! 2067 FORMAT ('TSGRID: H_DIFFU(1) = ',F14.7)
! 2068 FORMAT ('TSGRID: H_DIFFU(2) = ',F14.7)
! 2069 FORMAT ('TSGRID: H_DIFFU(3) = ',F14.7)
!!==============================================================================

!!! added by Kang for outputing data to .dat file @ 2020
 2140 format ("Year      Month     Day      Hour      X(m)      Y(m)      Z(m)      ",22A10)
 2150 format (I4,3I3,3F12.3,22E16.7)
!!==============================================================================

!     End of subroutine TSGRID

      end subroutine tsgrid

max_stars_repo_name: mdhillmancmcl/TheWorldAvatar-CMCL-Fork

max_issues_repo_name: mdhillmancmcl/TheWorldAvatar-CMCL-Fork

max_forks_repo_name: mdhillmancmcl/TheWorldAvatar-CMCL-Fork


! ###############################################################################################################################
! Begin MIT license text.
! _______________________________________________________________________________________________________

! Copyright 2019 Dr William R Case, Jr (dbcase29@gmail.com)

! Permission is hereby granted, free of charge, to any person obtaining a copy of this software and
! associated documentation files (the "Software"), to deal in the Software without restriction, including
! without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to
! the following conditions:

! The above copyright notice and this permission notice shall be included in all copies or substantial
! portions of the Software and documentation.

! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
! OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
! THE SOFTWARE.
! _______________________________________________________________________________________________________

! End MIT license text.

   MODULE BD_PARAM_Interface

   INTERFACE

      SUBROUTINE BD_PARAM ( CARD )

      USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG, DOUBLE
      USE IOUNT1, ONLY                :  WRT_ERR, WRT_LOG, ERR, F04, F06

      USE SCONTR, ONLY                :  BLNK_SUB_NAM, ECHO, FATAL_ERR, IERRFL, JCARD_LEN, JF, MEPSIL, MPBARLU, NUM_USETSTR,       &
                                         WARN_ERR
      USE TIMDAT, ONLY                :  TSEC
      USE SUBR_BEGEND_LEVELS, ONLY    :  BD_PARAM_BEGEND
      USE CONSTANTS_1, ONLY           :  ZERO, ONE
      USE MACHINE_PARAMS, ONLY        :  MACH_PREC
      USE DOF_TABLES, ONLY            :  TSET_CHR_LEN

      USE PARAMS, ONLY                :  ARP_TOL         , ART_KED         , ART_ROT_KED     , ART_TRAN_KED    ,                   &
                                         ART_MASS        , ART_ROT_MASS    , ART_TRAN_MASS   , AUTOSPC         , AUTOSPC_NSET    , &
                                         AUTOSPC_RAT     , AUTOSPC_INFO    , AUTOSPC_SPCF    , BAILOUT         , CRS_CCS         , &
                                         CBMIN3          , CBMIN4          , CBMIN4T         , CHKGRDS         ,                   &
                                         CUSERIN         , CUSERIN_EID     , CUSERIN_IN4     , CUSERIN_PID     , CUSERIN_SPNT_ID , &
                                         CUSERIN_XSET    , CUSERIN_COMPTYP , DARPACK         ,                                     &
                                         DELBAN          , EIGESTL         , EIGNORM2        , ELFORCEN        , EPSERR          , &
                                         EQCHK_REF_GRID  , EQCHK_NORM      , EQCHK_OUTPUT    , EQCHK_TINY      ,                   &
                                         EPSIL           , EMP0_PAUSE      , ESP0_PAUSE      , F06_COL_START   ,                   &
                                         GRDPNT          , GRIDSEQ         , HEXAXIS         ,                                     &
                                         IORQ1M          , IORQ1S          , IORQ1B          , IORQ2B          , IORQ2T          , &
                                         ITMAX           , KLLRAT          , KOORAT          , LANCMETH        , MATSPARS        , &
                                         MEMAFAC         , MIN4TRED        , MXALLOCA        , MAXRATIO        ,                   &
                                         MEFMCORD        , MEFMLOC         , MEFMGRID        ,                                     &
                                         MPFOUT          , MXITERI         , MXITERL         , OTMSKIP         , POST            , &
                                         PBARLDEC        , PBARLSHR        , PCOMPEQ         , PCHSPC1         , PCMPTSTM        , &
                                         PRTBASIC        , PRTCONN         , PRTCORD         , PRTDISP         , PRTDLR          , &
                                         PRTDOF          , PRTFOR          , PRTHMN          , PRTGMN          , PRTGOA          , &
                                         PRTCGLTM        , PRTPHIZL        , PRTIFLTM        , PRTKXX          , PRTMXX          , &
                                         PRTOU4          , PRTPHIXA        , PRTMASS         , PRTMASSD        , PRTRMG          , &
                                         PRTSCP          , PRTPSET         , PRTTSET         , PRTUSET         ,                   &
                                         PRTSTIFD        , PRTSTIFF        , PRTUO0          ,                                     &
                                         PRTYS           , PRTQSYS         ,                                                       &
                                         Q4SURFIT        , QUADAXIS        , QUAD4TYP        , RCONDK          , RELINK3         , &
                                         SEQPRT          , SEQQUIT         , SETLKTM         , SETLKTK         , SHRFXFAC        , &
                                         SKIPMGG         , SOLLIB          , SPARSE_FLAVOR   , SPARSTOR        ,                   &
                                         SPC1QUIT        , SORT_MAX        , SPC1SID         , STR_CID                           , &
                                         SUPINFO         , SUPWARN                                                               , &
                                         THRESHK         , THRESHK_LAP     , TINY            ,                                     &
                                         TSTM_DEF        , USR_JCT         , USR_LTERM_KGG   , USR_LTERM_MGG   , WINAMEM         , &
                                         WTMASS

      IMPLICIT NONE

      CHARACTER(LEN=*), INTENT(IN)    :: CARD              ! A Bulk Data card
      CHARACTER(15*BYTE)              :: PARNAM            ! The name of a parameter

      INTEGER(LONG)                   :: LOWER             ! Lower allowable value for an integer parameter
      INTEGER(LONG)                   :: UPPER             ! Upper allowable value for an integer parameter
      INTEGER(LONG), PARAMETER        :: SUBR_BEGEND = BD_PARAM_BEGEND

      END SUBROUTINE BD_PARAM

   END INTERFACE

   END MODULE BD_PARAM_Interface


text:
      program unstr
      integer a, b

c     10 j = 3
 10   j = 3
      goto 20
c     20   continue: could be fuse with node 30
 20   continue
      continue

c     if (a .eq. 0) goto 30
 30   if (a .eq. 0) goto 30
c     if:
      if (b .eq. 0) goto 20
      goto 10

      end

ext: f


text:
This editor can edit this entry and tell us a bit about themselves by clicking the Edit icon.

20110502 18:53:24 nbsp Welcome to the Wiki. In regard to the pinball machines, Why?. Can you provide more detail, are the eating quarters, or poorly maintained, or do they shock the players, or are they harboring zombies? Users/JasonAller

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------

MODULE m_radovlp
  CONTAINS
  SUBROUTINE rad_ovlp(atoms,usdus,input,hub1inp,vr,epar, uun21,udn21,dun21,ddn21)
    !***********************************************************************
    ! calculates the overlap of the radial basis functions with different
    ! spin directions. These overlapp integrals are needed to calculate
    ! the contribution to the hamiltonian from the constant constraint
    ! B-field.
    !
    ! Philipp Kurz 2000-04
    !***********************************************************************

    USE m_int21
    USE m_radfun
    USE m_types
    IMPLICIT NONE
    TYPE(t_input),INTENT(IN)   :: input
    TYPE(t_atoms),INTENT(IN)   :: atoms
    TYPE(t_hub1inp),INTENT(IN) ::hub1inp
    TYPE(t_usdus),INTENT(INOUT):: usdus

    !     .. Array Arguments ..
    REAL,    INTENT  (IN):: epar(0:,:,:)!(0:atoms%lmaxd,atoms%ntype,input%jspins)
    REAL,    INTENT  (IN):: vr(:,0:,:,:)!(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins)
    REAL,    INTENT (OUT):: uun21(0:,:),udn21(0:,:)!(0:atoms%lmaxd,atoms%ntype)
    REAL,    INTENT (OUT):: dun21(0:,:),ddn21(0:,:)!(0:atoms%lmaxd,atoms%ntype)
    !     ..
    !     .. Local Scalars ..
    INTEGER itype,l,ispin,noded,nodeu,i
    REAL    wronk
    LOGICAL l_hia
    !     ..
    !     .. Local Arrays ..
    REAL :: f(atoms%jmtd,2,0:atoms%lmaxd,input%jspins)
    REAL :: g(atoms%jmtd,2,0:atoms%lmaxd,input%jspins)
    REAL :: vrTmp(atoms%jmtd)
    !     ..
    DO itype = 1,atoms%ntype
       DO l = 0,atoms%lmax(itype)
          l_hia=.FALSE.
          DO i = atoms%n_u+1, atoms%n_u+atoms%n_hia
             IF(atoms%lda_u(i)%atomType.EQ.itype.AND.atoms%lda_u(i)%l.EQ.l) THEN
                l_hia=.TRUE.
             ENDIF
          ENDDO
          DO ispin = 1,input%jspins

             !In the case of a spin-polarized calculation with Hubbard 1 we want to treat
             !the correlated orbitals with a non-spin-polarized basis
             IF(l_hia.AND.input%jspins.EQ.2.AND..NOT.hub1inp%l_dftspinpol) THEN
                vrTmp = (vr(:,0,iType,1) + vr(:,0,iType,2))/2.0
             ELSE
                vrTmp = vr(:,0,iType,ispin)
             ENDIF

             CALL radfun(l,iType,ispin,epar(l,iType,ispin),vrTmp,atoms,&
                 f(1,1,l,ispin),g(1,1,l,ispin),usdus,nodeu,noded,wronk)

          ENDDO
          CALL int_21(f,g,atoms,itype,l,uun21,udn21,dun21,ddn21)
       ENDDO
    ENDDO

  END SUBROUTINE rad_ovlp
END

max_stars_repo_name: MRedies/FLEUR

max_issues_repo_name: MRedies/FLEUR

max_forks_repo_name: MRedies/FLEUR


text:
      SUBROUTINE proma2(a,b,c,n1,n2,n3)
!********************************************************************
!
!***this routine evaluates a matrix product
!                                                t
!           a(i,j) = b(i,k) * c(j,k)    a = b * c
!
!********************************************************************
      IMPLICIT NONE

      INTEGER (kind=4),INTENT(IN) :: n1,n2,n3
      REAL (kind=8),INTENT(IN) :: b(1:n1,1:n3), c(1:n2,1:n3)
      REAL (kind=8),INTENT(OUT) :: a(1:n1,1:n2)

      a(1:n1,1:n2) = MATMUL(b(1:n1,1:n3),TRANSPOSE(c(1:n2,1:n3)))
      RETURN

      END SUBROUTINE proma2

max_stars_repo_name: jerebenitez/IFE-simpact-openfoam

max_issues_repo_name: jerebenitez/IFE-simpact-openfoam

max_forks_repo_name: jerebenitez/IFE-simpact-openfoam


! Gamma Probability Computation
!   Hiroyuki Okamura, Hiroshima University
!   okamu@rel.hiroshima-u.ac.jp

module gamma
  implicit none
!  double precision, parameter :: npi = 3.14159265358979324d0
  double precision, private, parameter :: log2pi = 1.83787706640934548d0
  double precision, private, parameter :: logpi = 1.14472988584940017d0
  integer, private, parameter :: CN = 8
  double precision, private, parameter :: &
    B0 = 1.0d0, B1 = -1.0d0 / 2.0d0, &
    B2 = 1.0d0 / 6.0d0, B4 = -1.0d0 / 30.0d0, &
    B6 = 1.0d0 / 42.0d0, B8 = -1.0d0 / 30.0d0, &
    B10 = 5.0d0 / 66.0d0, B12 = -691.0d0 / 2730.0d0, &
    B14 = 7.0d0 / 6.0d0, B16 = -3617.0d0 / 510.0d0
  integer, private, parameter :: factmax = 20
  double precision, private, parameter :: nfact(0:factmax) = (/ &
    1.0d0, 1.0d0, 2.0d0, 6.0d0, 24.0d0, &
    120.0d0, 720.0d0, 5040.0d0, 40320.0d0, 362880.0d0, 3628800.0d0, &
    39916800.0d0, 479001600.0d0, 6227020800.0d0, 87178291200.0d0, &
    1307674368000.0d0, 20922789888000.0d0, 355687428096000.0d0, &
    6402373705728000.0d0, 121645100408832000.0d0, &
    2432902008176640000.0d0 /)
  double precision, private, parameter :: lognfact(0:factmax) = (/ &
    0.0d0, 0.0d0, &
    0.6931471805599453d0, 1.791759469228055d0, &
    3.1780538303479458d0, 4.787491742782046d0, 6.579251212010101d0, &
    8.525161361065415d0, 10.60460290274525d0, 12.801827480081469d0, &
    15.104412573075516d0, 17.502307845873887d0, 19.987214495661885d0, &
    22.552163853123425d0, 25.19122118273868d0, 27.89927138384089d0, &
    30.671860106080672d0, 33.50507345013689d0, 36.39544520803305d0, &
    39.339884187199495d0, 42.335616460753485d0 /)
  public fact, logfact, loggamma, psi, polygamma

contains
  function fact(s)
    integer, intent(in) :: s
    double precision :: fact
    if (s <= factmax) then
       fact = nfact(s)
    else
       fact = exp(loggamma(dble(s + 1)))
    end if
  end function fact

  function logfact(s)
    integer, intent(in) :: s
    double precision :: logfact
    if (s <= factmax) then
       logfact = lognfact(s)
    else
       logfact = loggamma(dble(s + 1))
    end if
  end function logfact

  function loggamma(s)
    double precision, intent(in) :: s
    double precision :: loggamma
    double precision :: x, v, w
    x = s
    v = 1.0d0
    do while (x < CN)
       v = v * x
       x = x + 1.0d0
    end do
    w = 1.0d0 / (x * x)
    loggamma = ((((((((B16 / (16.0d0 * 15.0d0)) * w &
         + (B14 / (14.0d0 * 13.0d0))) * w &
         + (B12 / (12.0d0 * 11.0d0))) * w &
         + (B10 / (10.0d0 * 9.0d0))) * w &
         + (B8 / (8.0d0 * 7.0d0))) * w &
         + (B6 / (6.0d0 * 5.0d0))) * w &
         + (B4 / (4.0d0 * 3.0d0))) * w &
         + (B2 / (2.0d0 * 1.0d0))) / x &
         + 0.5d0 * log2pi - log(v) - x + (x - 0.5d0) * log(x)
    return
  end function loggamma

  function psi(s)
    double precision, intent(in) :: s
    double precision :: psi
    double precision :: x, v, w
    x = s
    v = 0.0d0
    do while (x < CN)
       v = v + 1.0d0 / x
       x = x + 1.0d0
    end do
    w = 1.0d0 / (x * x)
    v = v + ((((((((B16 / 16.0d0) * w &
         + (B14 / 14.0d0)) * w &
         + (B12 / 12.0d0)) * w &
         + (B10 / 10.0d0)) * w &
         + (B8 / 8.0d0)) * w &
         + (B6 / 6.0d0)) * w &
         + (B4 / 4.0d0)) * w &
         + (B2 / 2.0d0)) * w &
         + 0.5d0 / x
    psi = log(x) - v
  end function psi

  function polygamma(n, s)
    integer, intent(in) :: n
    double precision, intent(in) :: s
    double precision :: polygamma
    integer :: k
    double precision :: x, t, u, v, w
    x = s
    u = 1.0d0
    do k = 1-n, -1
       u = u * k
    end do
    v = 0.0d0
    do while (x < CN)
       v = v + 1.0d0 / x**(n+1)
       x = x + 1.0d0
    end do
    w = x * x
    t = (((((((B16 &
         * (n + 15.0d0) * (n + 14.0d0) / (16.0d0 * 15.0d0 * w) + B14) &
         * (n + 13.0d0) * (n + 12.0d0) / (14.0d0 * 13.0d0 * w) + B12) &
         * (n + 11.0d0) * (n + 10.0d0) / (12.0d0 * 11.0d0 * w) + B10) &
         * (n +  9.0d0) * (n +  8.0d0) / (10.0d0 *  9.0d0 * w) + B8) &
         * (n +  7.0d0) * (n +  6.0d0) / ( 8.0d0 *  7.0d0 * w) + B6) &
         * (n +  5.0d0) * (n +  4.0d0) / ( 6.0d0 *  5.0d0 * w) + B4) &
         * (n +  3.0d0) * (n +  2.0d0) / ( 4.0d0 *  3.0d0 * w) + B2) &
         * (n +  1.0d0) * n / (2.0d0 * 1.0d0 * w) + 0.5d0 * n / x + 1.0d0
    polygamma = u * (t / x**n + n*v)
  end function polygamma

end module gamma


text:
\ fib.f         - x4 fibonacci benchmark words
\ ------------------------------------------------------------------------

 .( loading fib.f ) cr

\ ------------------------------------------------------------------------

9000000 var itters          \ 9 million

0 var counter

\ ------------------------------------------------------------------------
\ removed from the kernel because i discourage its use and this example
\ is the only place i use it

: recurse
  last name> , ;            ; foot in self shoot

\ ------------------------------------------------------------------------
\ the recursive method - this we call once not 9 milion times :)

\ this word recurses on itself one hundred and thirteen million times
\ to calculate the 40th fib (bleh)

: fib1          ( n1 --- )
  dup 1 >
  if
    dup  1- recurse
    swap 2- recurse
    +
  then ;

\ ------------------------------------------------------------------------
\ this code is more proof that anything you can do with recursion
\   can be done better without.

1 var f1
0 var f2

: fib2
  off> f2
  1 !> f1
  1
  ?do
    f2 f1 +
    f1 !> f2
    !> f1
  loop f1 ;

\ ------------------------------------------------------------------------
\ this is my version of the above itterative method

: fib3
  0 1 rot 1
  ?do
    tuck +
  loop nip ;

\ ------------------------------------------------------------------------

: (fib4) tuck + ;
: fib4 1 1 rot 1- rep (fib4) drop ;

\ ------------------------------------------------------------------------

: fib-bench
  cr ." fib 1 " timer-reset             40 fib1 drop      .elapsed
  cr ." fib 2 " timer-reset itters 0 do 40 fib2 drop loop .elapsed
  cr ." fib 3 " timer-reset itters 0 do 40 fib3 drop loop .elapsed
  cr ." fib 4 " timer-reset itters 0 do 40 fib4 drop loop .elapsed ;

\ fib1 rubs in  4:06 on my k6-3/550
\ fib2 runs in 17:56
\ fib3 runs in  0:32

\ ========================================================================


text:
C    Copyright(C) 2008-2017 National Technology & Engineering Solutions
C    of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C    NTESS, the U.S. Government retains certain rights in this software.
C
C    Redistribution and use in source and binary forms, with or without
C    modification, are permitted provided that the following conditions are
C    met:
C
C    * Redistributions of source code must retain the above copyright
C       notice, this list of conditions and the following disclaimer.
C
C    * Redistributions in binary form must reproduce the above
C      copyright notice, this list of conditions and the following
C      disclaimer in the documentation and/or other materials provided
C      with the distribution.
C
C    * Neither the name of NTESS nor the names of its
C      contributors may be used to endorse or promote products derived
C      from this software without specific prior written permission.
C
C    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C    A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C    OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C    LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C    DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C    THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
C
      SUBROUTINE MXEXEC (MYV, MYCV, MYLOC, MYCLOC, UCLOC, COFFST,
     *   OFFSET, DPOINT, LDICT, NNAMES,
     *   VOID, LVOID, NVOIDS, FILL, FDATA, CFILL, CFDATA, CHRNUM,
     *   CHRCOL, LASTER)
C
      IMPLICIT INTEGER (A-Z)
      INCLUDE 'params.inc'
C
C     This routine satisfies deferred memory requests.
C     It will service both numeric and character
C     deferred requests in the mixed mode.  In the nonmixed mode,
C     character memory is not deferred.
C
C***********************************************************************
C
C     MYV      Internal reference vector.
               DIMENSION MYV(*)
C     MYCV     Internal reference array.
               CHARACTER MYCV(*)
C     MYLOC    Address of internal reference vector.
C     MYCLOC   Address of internal character array.
C     UCLOC    Address of user's character array.
C     COFFST   Offset between internal numeric array and user's
C              character array.
C     OFFSET   Offset address between internal and user arrays.
C     DPOINT   Dictionary pointer table.
C     LDICT    Dimension of dictionary table.
C     NNAMES   Number of names in dictionary.
               DIMENSION DPOINT(LDICT,CHRCOL,3), NNAMES(2)
C     VOID     Void table.
C     LVOID    Dimension of void table.
C     NVOIDS   Number of voids.
               DIMENSION VOID(LVOID,CHRCOL,2), NVOIDS(2)
C     FILL     Flag for data fill.
C     FDATA    Data for fill.
               LOGICAL FILL
C     CFILL    Flag for character data fill.
C     CFDATA   Data for fill.
               LOGICAL CFILL
               CHARACTER*1 CFDATA
C     CHRNUM   Number of characters per numeric storage unit
C     CHRCOL   Column for character tables.
C     LASTER   Error return.
C
C***********************************************************************
C
      LASTER = SUCESS
      MEM = 0
      DO 100 IDICT = 1, NNAMES(1)
         MEM = MEM + MIN(0, DPOINT(IDICT,1,2))
  100 CONTINUE
      IF (MEM .EQ. 0) RETURN
      MEM = - MEM
C
C
      CALL MXGET (MYLOC, MEM, VOID, LVOID, NVOIDS,
     *   CHRCOL, LASTER, VROW)
      IF (LASTER .NE. SUCESS) RETURN
C
C     Now satisfy all the deferred requests.
C
      DO 130 IDICT = 1, NNAMES(1)
         IF (DPOINT(IDICT,1,2) .LT. 0) THEN
            IF (DPOINT(IDICT,1,3) .EQ. -1) THEN
               MYV(DPOINT(IDICT,1,1) - MYLOC + 1) =
     *            VOID(VROW,1,1) + OFFSET
            ELSE
               MYV(DPOINT(IDICT,1,1) - MYLOC + 1) =
     *            (VOID(VROW,1,1) - 1) * CHRNUM + 1 + COFFST
            END IF
            VOID(VROW,1,2) = VOID(VROW,1,2) + DPOINT(IDICT,1,2)
            DPOINT(IDICT,1,2) = - DPOINT(IDICT,1,2)
            DPOINT(IDICT,1,1) = VOID(VROW,1,1)
            VOID(VROW,1,1) = VOID(VROW,1,1) + DPOINT(IDICT,1,2)
C
C           Perform data fill if appropriate.
C
            IF (FILL .AND. DPOINT(IDICT,1,3) .EQ. -1) THEN
               DO 110 I = DPOINT(IDICT,1,1),
     *            DPOINT(IDICT,1,1)+DPOINT(IDICT,1,2)-1
                  MYV(I) = FDATA
  110          CONTINUE
            ELSE IF (CFILL .AND. DPOINT(IDICT,1,3) .GT. 0) THEN
               TLOC = (DPOINT(IDICT,1,1) - 1) * CHRNUM + 1 + COFFST
     *            + UCLOC - MYCLOC
               DO 120 I = TLOC, TLOC + DPOINT(IDICT,1,3) - 1
                  MYCV(I) = CFDATA
  120          CONTINUE
            END IF
C
         END IF
  130 CONTINUE
C
      CALL VTABLE (0, 0, VOID, LVOID, NVOIDS(1), CHRCOL, LASTER)
      RETURN
      END

max_stars_repo_name: mathstuf/seacas

max_stars_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_issues_repo_name: mathstuf/seacas

max_issues_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_forks_repo_name: mathstuf/seacas

max_forks_repo_licenses: BSD-3-Clause NetCDF Zlib MIT


text:
Describe Users/CristalRangel here.

20130715 04:40:28 nbsp I moved your roommate post from the wikispot.org site...and into the Roommates Wanted page. Users/tcnsc


text:
! RUN: %python %S/../test_errors.py %s %flang -fopenacc

! Check OpenACC clause validity for the following construct and directive:
!   2.8 host_data

program openacc_host_data_validity

  implicit none

  integer, parameter :: N = 256
  real(8), dimension(N, N) :: aa, bb
  logical :: ifCondition = .TRUE.

  !ERROR: At least one of USE_DEVICE clause must appear on the HOST_DATA directive
  !$acc host_data
  !$acc end host_data

  !$acc host_data use_device(aa)
  !$acc end host_data

  !$acc host_data use_device(aa) if(.true.)
  !$acc end host_data

  !$acc host_data use_device(aa) if(ifCondition)
  !$acc end host_data

  !$acc host_data use_device(aa, bb) if_present
  !$acc end host_data

  !ERROR: At most one IF_PRESENT clause can appear on the HOST_DATA directive
  !$acc host_data use_device(aa, bb) if_present if_present
  !$acc end host_data

  !$acc host_data use_device(aa, bb) if(.true.) if_present
  !$acc end host_data

  !ERROR: At most one IF clause can appear on the HOST_DATA directive
  !$acc host_data use_device(aa, bb) if(.true.) if(ifCondition)
  !$acc end host_data

end program openacc_host_data_validity

max_stars_repo_name: mkinsner/llvm

max_issues_repo_name: mkinsner/llvm

max_forks_repo_name: mkinsner/llvm


text:
      FUNCTION SADOTB( A, B )
      REAL A(3), B(3)
C*****
C  SINGLE-PRECISION VERSION
C
C  DOT PRODUCT A . B
C*****
      SADOTB = A(1)*B(1) + A(2)*B(2) + A(3)*B(3)
      RETURN
      END


text:
program array3
implicit none
integer,parameter :: imax =11
real,dimension (imax) :: x
real                  :: average,sum
integer               ::i
print*, "enter","imax", "numbers"

sum =0.0
	do i = 1,imax
		read*, x(i)
			sum = sum + x(i)
	end do
average = sum /imax
print*, "the average is ", average
print*, "the numbers are"
print*,x

end program array3

max_stars_repo_name: akkiturb/Ilkay-s-Fortran-Class


             program demo_alert
             use M_attr, only : alert, attr, attr_mode
             implicit none
             real X
                call attr_mode(manner='plain')
                call attr_mode(manner='color')
                call alert("error", "Say you didn't!")
                call alert("warn",  "I wouldn't if I were you, Will Robinson.")
                call alert("info",  "I fixed that for you, but it was a bad idea.")
                call alert("debug", "Who knows what is happening now?.")
                call alert("???    ",  "not today you don't")
                ! call to just update the macros
                call alert()
                ! conventional call to ATTR(3f) using the ALERT(3f)-defined macros
                write(*,*)attr('<bo>The year was <g><YE></g>, the month was <g><MO></g>')
                ! optional arguments
                X=211.3
                call alert('error','allowed range of X is 0 <lt> X <lt> 100, X=<r>',X)
                ! up to twenty values are allowed of intrinsic type
                call alert('info','values are<g>',10,234.567,cmplx(11.0,22.0),123.456d0,'</g>today')
             end program demo_alert


text:
module test2019_accessibility_attribute

  contains

 subroutine subA()
!
     implicit none
! Local variables
     type node_on_node_T
     end type node_on_node_T

  end subroutine subA

 end module test2019_accessibility_attribute

max_stars_repo_name: ouankou/rose

max_issues_repo_name: WildeGeist/rose

max_forks_repo_name: WildeGeist/rose


text:
c=========================================================================
      subroutine setprob
c=========================================================================

      use geoclaw_module
      use topo_module
      use dtopo_module
      use qinit_module
      use auxinit_module

      implicit double precision (a-h,o-z)

      call set_geo          !# sets basic parameters g and coord system
      call set_tsunami      !# sets parameters specific to tsunamis
      call set_flow_grades  !# sets refinement parameters
      call set_qinit        !# specifies file(s) for initial conditions
      call set_auxinit      !# specifies file(s) for auxiliary variables
      call set_topo         !# specifies topography (bathymetry) files
      call set_dtopo        !# specifies moving topo 'dtopo' from earthquake
      call setregions       !# specifies where refinement is allowed/forced
      call setgauges        !# locations of measuring gauges
      call setfixedgrids    !# specifies output on arbitrary uniform fixed grids

      return
      end

max_issues_repo_name: che-wenchao/D-Claw

max_forks_repo_name: che-wenchao/D-Claw


text:
	SUBROUTINE GD_OPNR  ( filnam, igdfln, navsz,  rnvblk, ianlsz,
     +			      anlblk, ihdrsz, maxgrd, iret )
C************************************************************************
C* GD_OPNR								*
C*									*
C* This subroutine opens an existing GEMPAK grid file for realtime	*
C* data access.  The file is opened with shared, write access.		*
C*									*
C* GD_OPNR  ( FILNAM, IGDFLN, NAVSZ, RNVBLK, IANLSZ, ANLBLK, IHDRSZ,	*
C*            MAXGRD, IRET )						*
C*									*
C* Input parameters:							*
C*	FILNAM		CHAR*		File name 			*
C*									*
C* Output parameters:							*
C*	IGDFLN		INTEGER		File number			*
C*	NAVSZ		INTEGER		Navigation block length		*
C*	RNVBLK (NAVSZ)	REAL		Navigation block		*
C*	IANLSZ		INTEGER		Analysis block length		*
C*	ANLBLK (IANLSZ)	REAL		Analysis block			*
C*	IHDRSZ		INTEGER		Grid header length		*
C*	MAXGRD		INTEGER		Maximum number of grids		*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C*					 -2 = file cannot be opened	*
C*					 -7 = not a GEMPAK5 grid file	*
C*					 -8 = nav cannot be read	*
C*					-13 = grid header too long	*
C*					-14 = file name is blank	*
C**									*
C* Log:									*
C* M. desJardins/GSFC	 6/87						*
C* M. desJardins/GSFC	 4/90	Added error for blank file name		*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
	INCLUDE		'GMBDTA.CMN'
	INCLUDE 	'grdcmn.cmn'
C
	CHARACTER*(*) 	filnam
	REAL		rnvblk (*), anlblk (*)
C*
	LOGICAL		wrtflg, shrflg
C------------------------------------------------------------------------
C*	Check for blank name.
C
	CALL ST_LSTR  ( filnam, lenf, ier )
	IF  ( lenf .eq. 0 )  THEN
	    iret = -14
	    RETURN
	END IF
C
C*	Open the file.
C
	wrtflg = .true.
	shrflg = .true.
	CALL GD_OFIL  ( filnam, wrtflg, shrflg, igdfln, navsz, rnvblk,
     +			ianlsz, anlblk, ihdrsz, maxgrd, iret )
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
      BLOCK DATA bd_gem_grid
C*************************************************************************
C
C  Information array for GEMPAK GRIB decoder Gribtogem version 2
C
C	INFO ARRAY ORDER:
C
C	1-6	= A-F grids (NMC id 21-26)
C	7	= G grid (NMC id 50)
C	8	= NMC id 5
C	9	= NMC id 6
C	10	= NMC id 27 (NH)
C	11	= NMC id 28 (SH)
C	12	= NMC id 100
C	13	= NMC id 101
C	14	= NMC id 104
C	15	= NMC id 105
C	16-19	= ECMWF id 1-4 (NH)
C	20-23	= ECMWF id 9-12 (SH)
C	24-27	= ECMWF id 5-8 (35N-35S)
C	28	= ETA 80km Resolution id 211
C	29	= Grid 256-usually thickened NMC AVN thinned grids
C	30-33	= NMC id 61-64
C       34      = NMC id 87 US Maps/Ruc
C       35      = NMC id 212 40km ETA
C       36      = NMC id 207 Alaska grid
C       37      = NMC id 214 Alaska grid
C       38      = NMC id 2 Global 2.5
C       39      = NMC id 203 Alask AVN
C       40      = NMC id 215 20km ETA
C       41      = NMC id 94 ETA 29km
C       42      = NMC id 96 ETA 48km
C       43      = NMC id 99 ETA 15km
C       44      = NMC id 192 ETA 10km West
C       45      = NMC id 196 ETA 10km East
C       46      = NMC id 3
C       47      = NMC 201
C       48	= NMC 202
C       49      = NMC 204
C       50	= NMC 205
C       51      = NMC 236 RUC2
C       52      = NMC 213 National CONUS double res
C       53	= NMC 103 N. Hemisphere polar stereographic 105W
C       54      = NMC 255, Use gds
C
C Log
C Unknown
C P.Bruehl/Undiata	5/95	Gribtogem v.2.
C Chiz/Unidata		10/96	New grids
C***************************************************************
	INCLUDE 'gem_grid.inc'

C
C     nmc is the index into the info arrays for NMC grids
C
	DATA nmc/0, 38, 46, 0, 8, 9, 14*0, 1, 2, 3, 4, 5, 6, 10, 11,
     &         21*0, 7, 10*0, 30, 31, 32, 33, 22*0, 34, 6*0, 41, 0,
     &         42, 2*0, 43, 12, 13, 0, 53, 14, 15, 86*0, 44, 3*0, 45,
     &         4*0, 47, 48, 39, 49, 50, 0, 36, 3*0, 28, 35, 52, 37, 40,
     &         20*0, 51, 18*0, 54, 29/
C
C     ecmwf is the index into the info arrays for ECMWF grids
C
	DATA ecmwf/16, 17, 18, 19, 24, 25, 26, 27, 20, 21, 22, 23,
     &             242*0,54,0/
C
	DATA ukm/0,38,46,251*0,54,0/
C
C     subgrid is a flag to indicate part of a larger grid
C
	DATA subgrid/6*.TRUE., 9*.FALSE., 12*.TRUE.,.FALSE.,.TRUE.,
     &	 	     4*.TRUE.,21*.FALSE./
C
C     tgrid_id is the grid id to check for compatibility
C
	DATA tgrid_id/4*1, 2*5, 7, 8, 9, 10, 11, 12, 13, 14, 15, 8*16,
     &     4*20, 28, 29, 4*30, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43,
     &     44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54/
C
C     scol is the starting column for subgrids
C
	DATA scol/1, 37, 1, 37, 1, 1, 9*0, 109, 73, 37, 1,
     &     109, 73, 37, 1, 109, 73, 37, 1, 0, 9999, 1, 91, 1, 91, 21*0 /
C
C     srow is the starting row for subgrids
C
	DATA srow/37, 37, 1, 1, 19, 1, 9*0, 4*37, 8*1, 0, 9999,
     &       46, 46, 1, 1, 21*0/
C
C     polepos is the position of the pole point (if applicable)
C	No longer needed, pole expansion done by {gribto}
C
C	DATA polepos/2*1333, 2*1, 1297, 1, 21*0, 0/
C
C     ngx is the number of longitude points in a large grid
C
	DATA ngx/4*73, 2*72, 9*0, 12*145, 0, 9999, 4*181, 21*0/
C
C     ngy is the number of latitude points in a large grid
C
	DATA ngy/4*73, 2*37, 9*0, 8*73, 4*29, 0, 9999, 4*91, 21*0/
C
C     rlatb is the bottom latitude of the grid
C
	DATA rlatb/6*-90.0, 20.0, 2*7.634903, -20.84062,
     &     20.84062, 17.09924, 10.51680, -0.268327, 17.5288,
     &     8*-90.0, 4*-35.0, 12.190, -90., 4*-90., 22.8756,
     &     12.190, 2*42.085, -90., 19.132, 12.190, 9.678, -3.441,
     &     17.912, 29.366, 23.476, -90., -20.826, 7.838, -25.000,
     &     0.616, 16.28, 7.838, 22.405, -9999./
C
C     rlatt is the top latitude of the grid
C
	DATA rlatt/6*90.0, 60.0, 42.8876, 44.27935, -20.84062,
     &     20.84062, 53.504, 44.48691, 32.7473, 53.7708,
     &     8*90.0, 4*35.0, 57.290, 90.,4*90., 46.0172, 57.290,
     &     2*63.976, 90., 57.634, 57.290, 53.708, 50.026, 45.508,
     &     49.150, 43.337, 90.,-20.826,35.617,60.644,45.620,55.48,
     &     35.617, 48.920, -9999. /

C
C     rlonl is the left longitude of the grid
C
	DATA rlonl/6*0.0, -140.0, 2*-133.4429, -125.0, 125.0,
     &  -129.2958, -137.1459, -139.475, -129.296, 12*0.0, -133.459,
     &  -30.,4*0., -120.4911, -133.459, 2*-175.641, 0., 174.163,
     &  -133.459, -128.826, -148.799, -98.303, -126.316, -96.745,
     &  0.,-150.000,-141.028,110.0,-84.904,-126.14,-141.028,-121.352,
     &  -9999. /
C
C     rlonr is the right longitude of the grid
C
	DATA rlonr/4*0.0,2*355.00,-52.5,2.102730,-23.74616,55.0,-55.0,
     &  -22.37376,-16.57655,-14.5993,-22.3738,12*0.0,-49.385,-30.,4*0.,
     &  -60.8284, -49.385, 2*-93.689, -2.5, -53.66, -49.385, -35.567,
     &  -3.230, -62.002, -99.844, -72.401, -1.0,30.000,-18.576,-109.129,
     &  -15.000,-57.38,-18.577,-50.811, -9999. /
C
C     ang1, ang2, ang3 are the angles for the polar stereographic grids
C
	DATA ang1/7*0, 3*90.0, -90.0, 4*90.0, 12*0, 25.,5*0, 90.0, 25.,
     &            2*90.0, 0, 90., 25., 41., 50., 33., 40., 34., 0,
     &            2*90.0, 0, 90.,25.,90.,90., -9999. /
	DATA ang2/7*0, 2*-105.0, -80.0, 100.0, 4*-105.0, 12*0, -95.,5*0,
     &       -105.0, -95., 2*-150.0,0,-150.,-95.,-97., -111., -83.,
     &       -115., -86., 0, 2*-105., 0., -60.,-95.,2*-105., -9999. /
	DATA ang3/27*0, 25., 5*0, 0, 25., 4*0, 25., 10*0,25.,2*0,-9999. /
C
C     proj is the projection type
C
	DATA proj/7*'CED',3*'STR','SPS',4*'STR',12*'CED','LCC',5*'CED',
     &            'STR', 'LCC', 2*'STR', 'CED', 'STR', 'LCC', 6*'CED',
     &            2*'STR','MER','STR','LCC',2*'STR','xxx' /

      END


text:
MODULE PolynomialBasisMappingModule

  USE KindModule, ONLY: &
    DP
  USE ProgramHeaderModule, ONLY: &
    nDOFX, nDOFZ
  USE QuadratureModule, ONLY: &
    xG5, wG5
  USE PolynomialBasisModule_Lagrange, ONLY: &
    IndL_Q, IndLX_Q, L_E, L_X1, L_X2, L_X3
  USE PolynomialBasisModule_Legendre, ONLY: &
    IndP_Q, IndPX_Q, P_E, P_X1, P_X2, P_X3, &
    MassP, MassPX

  IMPLICIT NONE
  PRIVATE

  PUBLIC :: InitializePolynomialBasisMapping
  PUBLIC :: MapNodalToModalX
  PUBLIC :: MapNodalToModalZ
  PUBLIC :: MapNodalToModal_Fluid
  PUBLIC :: MapNodalToModal_Radiation_X
  PUBLIC :: MapNodalToModal_Radiation
  PUBLIC :: MapModalToNodalX
  PUBLIC :: MapModalToNodalZ
  PUBLIC :: MapModalToNodal_Fluid
  PUBLIC :: MapModalToNodal_Radiation_X
  PUBLIC :: MapModalToNodal_Radiation

  REAL(DP), DIMENSION(:,:), ALLOCATABLE, PUBLIC :: Kij_X, Pij_X
  REAL(DP), DIMENSION(:,:), ALLOCATABLE         :: Kij_Z, Pij_Z

CONTAINS

  SUBROUTINE InitializePolynomialBasisMapping &
    ( Nodes_E, Nodes_X1, Nodes_X2, Nodes_X3 )

    REAL(DP), INTENT(in) :: Nodes_E (:)
    REAL(DP), INTENT(in) :: Nodes_X1(:)
    REAL(DP), INTENT(in) :: Nodes_X2(:)
    REAL(DP), INTENT(in) :: Nodes_X3(:)

    INTEGER :: i, j
    INTEGER :: qE, qX1, qX2, qX3

    ALLOCATE( Kij_X(1:nDOFX,1:nDOFX), Pij_X(1:nDOFX,1:nDOFX) )

    Kij_X = 0.0_DP
    DO j = 1, nDOFX
    DO i = 1, nDOFX

      DO qX3 = 1, SIZE( xG5 )
      DO qX2 = 1, SIZE( xG5 )
      DO qX1 = 1, SIZE( xG5 )

        Kij_X(i,j) &
          = Kij_X(i,j) &
              + wG5(qX1) * wG5(qX2) * wG5(qX3) &
                  * P_X1(IndPX_Q(1,i)) % P( xG5(qX1) ) &
                  * P_X2(IndPX_Q(2,i)) % P( xG5(qX2) ) &
                  * P_X3(IndPX_Q(3,i)) % P( xG5(qX3) ) &
                  * L_X1(IndLX_Q(1,j)) % P( xG5(qX1) ) &
                  * L_X2(IndLX_Q(2,j)) % P( xG5(qX2) ) &
                  * L_X3(IndLX_Q(3,j)) % P( xG5(qX3) )

      END DO
      END DO
      END DO

      Pij_X(i,j) &
        =   P_X1(IndPX_Q(1,j)) % P( Nodes_X1(IndLX_Q(1,i)) ) &
          * P_X2(IndPX_Q(2,j)) % P( Nodes_X2(IndLX_Q(2,i)) ) &
          * P_X3(IndPX_Q(3,j)) % P( Nodes_X3(IndLX_Q(3,i)) )

    END DO
    END DO

    ALLOCATE( Kij_Z(1:nDOFZ,1:nDOFZ), Pij_Z(1:nDOFZ,1:nDOFZ) )

    Kij_Z = 0.0_DP
    DO j = 1, nDOFZ
    DO i = 1, nDOFZ

      DO qX3 = 1, SIZE( xG5 )
      DO qX2 = 1, SIZE( xG5 )
      DO qX1 = 1, SIZE( xG5 )
      DO qE  = 1, SIZE( xG5 )

        Kij_Z(i,j) &
          = Kij_Z(i,j) &
              + wG5(qE) * wG5(qX1) * wG5(qX2) * wG5(qX3) &
                  * P_E (IndP_Q(0,i)) % P( xG5(qE ) ) &
                  * P_X1(IndP_Q(1,i)) % P( xG5(qX1) ) &
                  * P_X2(IndP_Q(2,i)) % P( xG5(qX2) ) &
                  * P_X3(IndP_Q(3,i)) % P( xG5(qX3) ) &
                  * L_E (IndL_Q(0,j)) % P( xG5(qE ) ) &
                  * L_X1(IndL_Q(1,j)) % P( xG5(qX1) ) &
                  * L_X2(IndL_Q(2,j)) % P( xG5(qX2) ) &
                  * L_X3(IndL_Q(3,j)) % P( xG5(qX3) )

      END DO
      END DO
      END DO
      END DO

      Pij_Z(i,j) &
        =   P_E (IndP_Q(0,j)) % P( Nodes_E (IndL_Q(0,i)) ) &
          * P_X1(IndP_Q(1,j)) % P( Nodes_X1(IndL_Q(1,i)) ) &
          * P_X2(IndP_Q(2,j)) % P( Nodes_X2(IndL_Q(2,i)) ) &
          * P_X3(IndP_Q(3,j)) % P( Nodes_X3(IndL_Q(3,i)) )

    END DO
    END DO

  END SUBROUTINE InitializePolynomialBasisMapping

  ! --- Nodal to Modal ---

  SUBROUTINE MapNodalToModalX( uN, uM )

    REAL(DP), DIMENSION(nDOFX), INTENT(in)  :: uN
    REAL(DP), DIMENSION(nDOFX), INTENT(out) :: uM

    uM = MassPX * MATMUL( Kij_X, uN )

  END SUBROUTINE MapNodalToModalX

  SUBROUTINE MapNodalToModalZ( uN, uM )

    REAL(DP), DIMENSION(nDOFZ), INTENT(in)  :: uN
    REAL(DP), DIMENSION(nDOFZ), INTENT(out) :: uM

    uM = MassP * MATMUL( Kij_Z, uN )

  END SUBROUTINE MapNodalToModalZ

  SUBROUTINE MapNodalToModal_Fluid( uN, uM )

    REAL(DP), DIMENSION(nDOFX), INTENT(in)  :: uN
    REAL(DP), DIMENSION(nDOFX), INTENT(out) :: uM

    uM = MassPX * MATMUL( Kij_X, uN )

  END SUBROUTINE MapNodalToModal_Fluid

  SUBROUTINE MapNodalToModal_Radiation_X( uN, uM )

    REAL(DP), DIMENSION(nDOFX), INTENT(in)  :: uN
    REAL(DP), DIMENSION(nDOFX), INTENT(out) :: uM

    uM = MassPX * MATMUL( Kij_X, uN )

  END SUBROUTINE MapNodalToModal_Radiation_X

  SUBROUTINE MapNodalToModal_Radiation( uN, uM )

    REAL(DP), DIMENSION(nDOFZ), INTENT(in)  :: uN
    REAL(DP), DIMENSION(nDOFZ), INTENT(out) :: uM

    uM = MassP * MATMUL( Kij_Z, uN )

  END SUBROUTINE MapNodalToModal_Radiation

  ! --- Modal to Nodal ---

  SUBROUTINE MapModalToNodalX( uN, uM )

    REAL(DP), DIMENSION(nDOFX), INTENT(out) :: uN
    REAL(DP), DIMENSION(nDOFX), INTENT(in)  :: uM

    uN = MATMUL( Pij_X, uM )

  END SUBROUTINE MapModalToNodalX

  SUBROUTINE MapModalToNodalZ( uN, uM )

    REAL(DP), DIMENSION(nDOFZ), INTENT(out) :: uN
    REAL(DP), DIMENSION(nDOFZ), INTENT(in)  :: uM

    uN = MATMUL( Pij_Z, uM )

  END SUBROUTINE MapModalToNodalZ

  SUBROUTINE MapModalToNodal_Fluid( uN, uM )

    REAL(DP), DIMENSION(nDOFX), INTENT(out) :: uN
    REAL(DP), DIMENSION(nDOFX), INTENT(in)  :: uM

    uN = MATMUL( Pij_X, uM )

  END SUBROUTINE MapModalToNodal_Fluid

  SUBROUTINE MapModalToNodal_Radiation_X( uN, uM )

    REAL(DP), DIMENSION(nDOFX), INTENT(out) :: uN
    REAL(DP), DIMENSION(nDOFX), INTENT(in)  :: uM

    uN = MATMUL( Pij_X, uM )

  END SUBROUTINE MapModalToNodal_Radiation_X

  SUBROUTINE MapModalToNodal_Radiation( uN, uM )

    REAL(DP), DIMENSION(nDOFZ), INTENT(out) :: uN
    REAL(DP), DIMENSION(nDOFZ), INTENT(in)  :: uM

    uN = MATMUL( Pij_Z, uM )

  END SUBROUTINE MapModalToNodal_Radiation

END MODULE PolynomialBasisMappingModule

max_stars_repo_name: srichers/thornado

max_issues_repo_name: srichers/thornado

max_forks_repo_name: srichers/thornado