!//=========================================================================
!// Oslo CTM3
!//=========================================================================
!// Based on UCI CTM core p-7.1 (1/2013).
!//
!// Amund Sovde Haslerud, September 2017
!//=========================================================================
!// Utilities for running the CTM, but more specific to Oslo CTM3.
!//=========================================================================
module utilities_oslo
  !// ----------------------------------------------------------------------
  !// Utilities for the Oslo chemistry.
  !// Try to avoid using modules here, this module is supposed to be
  !// compiled before other modules, as they may need the routines here.
  !//
  !// Contains:
  !//   subroutine gotoZC_IJ
  !//   subroutine backfromZC_IJ
  !//   subroutine ZC_MASS2CONC
  !//   subroutine ZC_CONC2MASS
  !//   real(r8) function troe
  !//   real(r8) function rate3B
  !//   subroutine get_chmcycles
  !//   subroutine us76_atmosphere
  !//   real(r8) function h2o_sat
  !//   subroutine source_e90
  !//   subroutine decay_e90
  !//   subroutine init_e90
  !//   subroutine tpause_e90
  !//   subroutine tpauseb_e90
  !//   subroutine tpauseb_o3
  !//   subroutine SZA_PN
  !//   subroutine stringUpCase
  !//   subroutine landfrac2mosaic
  !//   subroutine GROWSEASON
  !//   subroutine MAPPED_GROWSEASON
  !//
  !// Stefanie Falk, March - August 2018
  !//   Moved GROWSEASON to utilities.
  !//   Created mapping of daily preprocessed growseason.
  !// Amund Sovde Haslerud, September 2017
  !//   Some updates to use different tropopause definitions.
  !// Ole Amund Sovde, November 2015, October 2008
  !//
  !// ----------------------------------------------------------------------
  use cmn_precision, only: r8
  !// ----------------------------------------------------------------------
  implicit none
  !// ----------------------------------------------------------------------
  character(len=*), parameter, private :: f90file = 'utilities_oslo.f90'
  !// ----------------------------------------------------------------------
  public
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
contains
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine gotoZC_IJ(ZC_LOCAL,BTT,XSTT,I,J,II,JJ,L_START,L_END, &
       IPAR,JPAR,LPAR,IDBLK,JDBLK,MPBLKIB,MPBLKIE,MP, &
       NPAR, NOTRPAR, TRACER_ID_MAX, chem_idx, Xchem_idx)
    !// --------------------------------------------------------------------
    !// Putting BTT (transported tracers) and XSTT (non-transported tracers)
    !// into ZC_LOCAL, to avoid striding in chemistry integration routine.
    !// ONLY moves the column between L_START and L_END, so that this
    !// routine is applicable to both the troposphere and the stratosphere.
    !//
    !// Can consider only putting tropospheric components into ZC_LOCAL for
    !// tropospheric chemistry (stratospheric chemistry will also need some
    !// tropospheric components), but I doubt it will be any faster.
    !//
    !// Ole Amund Sovde, October 2008
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Input
    integer, intent(in) :: IPAR, JPAR, LPAR, IDBLK, JDBLK, MP, &
                           NPAR, NOTRPAR, TRACER_ID_MAX
    real(r8), intent(in), dimension(LPAR,NPAR,IDBLK,JDBLK) :: BTT
    real(r8), intent(in), dimension(LPAR,NOTRPAR,IPAR,JPAR):: XSTT
    integer, intent(in) :: I,J, II,JJ, L_START,L_END
    integer,dimension(IDBLK*JDBLK),intent(in) :: MPBLKIB, MPBLKIE
    integer, intent(in) :: chem_idx(NPAR), Xchem_idx(NOTRPAR)
    !// Output
    real(r8), intent(out), dimension(TRACER_ID_MAX,LPAR) :: ZC_LOCAL

    !// Locals
    integer :: TRACER_ID, N, L, K
    logical, parameter :: LOCALDEBUG = .false.
    !// --------------------------------------------------------------------
    character(len=*), parameter :: subr = 'gotoZC_IJ'
    !// --------------------------------------------------------------------

    !// Initialize (may want to remove it for efficiency)
    ZC_LOCAL(:,:) = 0._r8

    !// Transported species
    do N = 1,NPAR
       TRACER_ID = chem_idx(N)
       !// When tracer list was read, it was checked that TRACER_ID
       !// should be ok.
       do L = L_START, L_END
          ZC_LOCAL(TRACER_ID,L) = BTT(L,N,II,JJ)
          if (LOCALDEBUG) then
             !// Should not be necessary - rather test on BTT earlier.
             if (ZC_LOCAL(TRACER_ID,L).ne.ZC_LOCAL(TRACER_ID,L)) then
                write(6,'(a,5i5)') f90file//':'//subr// &
                     ': Tracer is NAN:',TRACER_ID,L,II,JJ,N,L_START, L_END
                do K = 1, LPAR
                   write(6,'(a,i5,es16.6)') '  L, BTT(L)',K,BTT(K,N,II,JJ)
                end do
                stop
             end if
          end if
       end do
    end do

    !// Non-transported species
    do N = 1,NOTRPAR
       TRACER_ID = Xchem_idx(N)
       !// When tracer list was read, it was checked that TRACER_ID
       !// should be ok.
       do L = L_START, L_END
          ZC_LOCAL(TRACER_ID,L) = XSTT(L,N,I,J)
       end do
    end do

    !// --------------------------------------------------------------------
  end subroutine gotoZC_IJ
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine backfromZC_IJ(ZC_LOCAL,BTT,XSTT,I,J,II,JJ,L_START, &
       L_END,IPAR,JPAR,LPAR,IDBLK,JDBLK,MPBLKIB,MPBLKIE,MP, &
       NPAR, NOTRPAR, TRACER_ID_MAX, chem_idx, Xchem_idx)
    !// --------------------------------------------------------------------
    !// Putting ZC_LOCAL back into BTT (transported tracers) and XSTT
    !// (non-transported tracers). See gotoZC_IJ.
    !// ONLY moves the column between L_START and L_END, so that this
    !// routine is applicable to both the troposphere and the stratosphere.
    !//
    !// Ole Amund Sovde, November 2014, October 2008
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Input
    integer, intent(in) :: IPAR, JPAR, LPAR, IDBLK, JDBLK, MP, &
                           NPAR, NOTRPAR, TRACER_ID_MAX
    real(r8), intent(in), dimension(TRACER_ID_MAX,LPAR) :: ZC_LOCAL
    integer, intent(in) :: I,J, II,JJ, L_START,L_END
    integer,dimension(IDBLK*JDBLK),intent(in) :: MPBLKIB, MPBLKIE
    integer, intent(in) :: chem_idx(NPAR), Xchem_idx(NOTRPAR)
    !// Output
    real(r8), intent(out), dimension(LPAR,NPAR,IDBLK,JDBLK) :: BTT
    real(r8), intent(out), dimension(LPAR,NOTRPAR,IPAR,JPAR):: XSTT

    !// Locals
    integer :: TRACER_ID, N, L
    !// --------------------------------------------------------------------
    character(len=*), parameter :: subr = 'backfromZC_IJ'
    !// --------------------------------------------------------------------

    !// Transported species back to BTT (no striding in BTT)
    do N = 1,NPAR
       TRACER_ID = chem_idx(N)
       !// When tracer list was read, it was checked that TRACER_ID
       !// should be ok.
       do L = L_START, L_END
          BTT(L,N,II,JJ) = ZC_LOCAL(TRACER_ID,L)
       end do
    end do

    !// Non-transported species
    do N = 1,NOTRPAR
       TRACER_ID = Xchem_idx(N)
       !// When tracer list was read, it was checked that TRACER_ID
       !// should be ok.
       do L = L_START, L_END
          XSTT(L,N,I,J) = ZC_LOCAL(TRACER_ID,L)
       end do
    end do

    !// --------------------------------------------------------------------
  end subroutine backfromZC_IJ
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine ZC_MASS2CONC(ZC_LOCAL,L_START,L_END,LPAR, &
          NPAR,NOTRPAR,TRACER_ID_MAX, trsp_idx, Xtrsp_idx, &
          TMASS,XTMASS,DV,AVOGNR)
    !// --------------------------------------------------------------------
    !// Converting ZC_LOCAL from mass to concentration.
    !// ONLY converts the column between L_START and L_END, so that this
    !// routine is applicable to both the troposphere and the stratosphere.
    !//
    !// The conversion uses molecular weight (MW [g/mol]), volume (DV [m3]),
    !// and Avogadros number (Na):
    !// ZC(molec/cm3) = ZC(kg) * 1/(MW * kg/(1e3g)) * Na / (DV * 1.e6cm3/m3)
    !//               = ZC(kg) * 1e3/MW * Na / DV / 1e6
    !//               = ZC(kg) * 1e-3/MW * Na / DV
    !//
    !// Ole Amund Sovde, October 2008, March 2012
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Input
    integer, intent(in) :: LPAR, NPAR, NOTRPAR, TRACER_ID_MAX
    real(r8), intent(in) :: DV(LPAR) !// Volume of box
    integer, intent(in) :: L_START, L_END
    integer, intent(in) :: trsp_idx(TRACER_ID_MAX), Xtrsp_idx(TRACER_ID_MAX)
    real(r8), intent(in) :: TMASS(NPAR), XTMASS(NOTRPAR)
    real(r8), intent(in) :: AVOGNR

    !// Input/Output
    real(r8), intent(inout), dimension(TRACER_ID_MAX,LPAR) :: ZC_LOCAL

    !// Locals
    integer :: TRACER_ID, N, XN, L
    real(r8) :: RDUM
    logical, parameter :: LOCALDEBUG = .false.
    !// --------------------------------------------------------------------
    character(len=*), parameter :: subr = 'ZC_MASS2CONC'
    !// --------------------------------------------------------------------

    !// Loop through chemical species
    !// While this produce some striding in ZC_LOCAL, it reduces the number
    !// of IF statements. Old routine checked on species for all
    !// tracers and levels.
    do TRACER_ID = 1, TRACER_ID_MAX

       N  = trsp_idx(TRACER_ID)  !// Transported index, if transported
       XN = Xtrsp_idx(TRACER_ID) !// Non-transported index, if not transported

       if (N .gt. 0) then
          !// Transported species
          RDUM = 1.e-3_r8 * AVOGNR / TMASS(N)   !// See above for comments
       else if (XN .gt. 0) then
          !// Non-transported species
          RDUM = 1.e-3_r8 * AVOGNR / XTMASS(XN) !// See above for comments
       else
          !// Species not in use, go to next species
          cycle
       end if

       !// Loop through partial column
       do L = L_START,L_END
          !// kg * 1d-3molec/kg / cm3 = molec/cm3
          ZC_LOCAL(TRACER_ID,L) = ZC_LOCAL(TRACER_ID,L) * RDUM / DV(L)
          if (LOCALDEBUG) then
             if (ZC_LOCAL(TRACER_ID,L) .ne. ZC_LOCAL(TRACER_ID,L)) then
                write(6,'(a,4i5,f6.1,es16.4)') f90file//':'//subr// &
                     ': ', TRACER_ID, L, N, XN, RDUM, DV(L)
                stop
             end if
          end if !// if (LOCALDEBUG) then

       end do !// do L = L_START, L_END
    end do !// do TRACER_ID = 1, TRACER_ID_MAX

    !// --------------------------------------------------------------------
  end subroutine ZC_MASS2CONC
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine ZC_CONC2MASS(ZC_LOCAL,L_START,L_END,LPAR, &
           NPAR,NOTRPAR,TRACER_ID_MAX, trsp_idx, Xtrsp_idx, &
           TMASS,XTMASS,DV,AVOGNR)
    !// --------------------------------------------------------------------
    !// Converting ZC_LOCAL from concentration to mass.
    !// ONLY converts the column between L_START and L_END, so that this
    !// routine is applicable to both the troposphere and the stratosphere.
    !//
    !// The conversion uses molecular weight (MW [g/mol]), volume (DV [m3]),
    !// and Avogadros number (Na):
    !// ZC(molec/cm3) = ZC(kg) * 1/(MW * kg/(1e3g) * Na / (DV * 1.e6cm3/m3)
    !//               = ZC(kg) * 1e3/MW * Na / DV / 1e6
    !//               = ZC(kg) * 1e-3/MW * Na / DV
    !//
    !// ZC(kg) = ZC(molec/cm3) * 1e3 * DV * MW / Na
    !//
    !// Ole Amund Sovde, October 2008
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Input
    integer, intent(in) :: LPAR, NPAR, NOTRPAR, TRACER_ID_MAX
    real(r8), intent(in) :: DV(LPAR) !// Volume of box
    integer, intent(in) :: L_START,L_END
    integer, intent(in) :: trsp_idx(TRACER_ID_MAX), Xtrsp_idx(TRACER_ID_MAX)
    real(r8), intent(in) :: TMASS(NPAR), XTMASS(NOTRPAR)
    real(r8), intent(in) :: AVOGNR
    !// Input/Output
    real(r8), intent(inout), dimension(TRACER_ID_MAX,LPAR) :: ZC_LOCAL

    !// Locals
    integer :: TRACER_ID, N, XN, L
    real(r8) :: RDUM
    !// --------------------------------------------------------------------
    character(len=*), parameter :: subr = 'ZC_CONC2MASS'
    !// --------------------------------------------------------------------

    !// Loop through chemical species
    !// While this produce some striding in ZC_LOCAL, it reduces the number
    !// of IF statements. Old routine checked on species for all
    !// tracers and levels.
    do TRACER_ID = 1, TRACER_ID_MAX

       N  = trsp_idx(TRACER_ID)  !// Transported index, if transported
       XN = Xtrsp_idx(TRACER_ID) !// Non-transported index, if not transported

       if (N .gt. 0) then
          !// Transported species
          RDUM = 1000._r8 * TMASS(N) / AVOGNR   !// See above for comments
       else if (XN .gt. 0) then
          !// Non-transported species
          RDUM = 1000._r8 * XTMASS(XN) / AVOGNR !// See above for comments
       else
          !// Species not in use, go to next species
          cycle
       end if

       !// Loop through partial column
       do L = L_START, L_END

          !// kg * 1d-3molec/kg / cm3 = molec/cm3
          ZC_LOCAL(TRACER_ID,L) = ZC_LOCAL(TRACER_ID,L) * RDUM * DV(L)

       end do !// do L = L_START, L_END

    end do !// do TRACER_ID = 1, TRACER_ID_MAX

    !// --------------------------------------------------------------------
  end subroutine ZC_CONC2MASS
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  real(r8) function TROE (KZERO,KINF,FC)
    !// --------------------------------------------------------------------
    !// Calculates ratecoeff. for three body reactions according to
    !// expressions developed by troe. see f.ex. atkinson et al. 1989
    !//
    !// This is the old tropospheric calculation of 3-body reactions,
    !// and is included for historical reasons. It is replaced by
    !// rate3B, which is more flexible.
    !// Amund Sovde, December 2014
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// I/O parameters
    real(r8), intent(in) :: KZERO,KINF,FC
    !// Local parameters
    real(r8) :: RHLP1,RHLP2,RHLP3
    !// --------------------------------------------------------------------
    character(len=*), parameter :: subr = 'TROE'
    !// --------------------------------------------------------------------
    write(6,'(a)') f90file//':'//subr//': Not to be used! Use rate3B instead'
    stop
    RHLP1 = KZERO / KINF
    RHLP2 = KZERO * KINF / (KZERO + KINF)
    RHLP3 = 1._r8 / (1._r8 + log10(RHLP1) * log10(RHLP1))
    TROE  = RHLP2 * (FC**RHLP3)
    !// --------------------------------------------------------------------
  end function TROE
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  real(r8) function rate3B(INUM, TZ300,AIRDENS,kzero300,Nzero, &
                         kinf300,Minf, FC,ACT)
    !// --------------------------------------------------------------------
    !// Three-body reaction rates. This is very much the same routine
    !// as the stratospheric code BODY3_BOX, but written a bit more
    !// cleanly.
    !// INUM:    Reaction number for debug use
    !// TZ300:   T/300
    !// AIRDENS: molec/cm3
    !// kzero300: k_o^300 in JPL
    !// kzero300: k_o^300 in JPL
    !// Nzero:    n in JPL
    !// Minf:     m in JPL
    !// FC:     constant, usually 0.6d0
    !// ACT:    0: Regular calculation (k_f)
    !//         1: Calculated as an activation channel (k_f^ca)
    !//
    !// Ole Amund Sovde, December 2014
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Input
    real(r8), intent(in)  :: &
         TZ300, AIRDENS, kzero300, Nzero, kinf300, Minf, FC
    integer, intent(in) :: INUM, ACT
    !// Local variables
    real(r8) :: kzero, kinf, B, D, RX
    !// --------------------------------------------------------------------
    character(len=*), parameter :: subr = 'rate3B'
    !// --------------------------------------------------------------------

    if (kinf300 .eq. 0._r8) then
       !// No high-pressure limit is given; reaction should
       !// only apply for low pressure.
       RX = kzero300 * TZ300**(-Nzero)
    else
       !// Assume both low and high-pressure limit is given
       kzero = kzero300 * TZ300**(-Nzero)
       kinf  = kinf300  * TZ300**(-Minf)
       B = kzero * AIRDENS / kinf
       D = 1._r8 / ( 1._r8 + log10(B) * log10(B) )
       RX = kzero / ( 1._r8 + B ) * (FC**D)
    end if

    if (ACT .eq. 0) then
       !// Standard 3-body reaction; multiply by air density.
       RATE3B = RX * AIRDENS
    else if (ACT .eq. 1) then
       !// Activation channel; skip multiplication by AIRDENS
       !// Check the JPL reactions if you are in doubt.
       RATE3B = RX
    else
       write(6,'(a, i7)') f90file//':'//subr// &
            ': ACT must be 0 or 1! INUM:', INUM
       stop
    end if

    !// --------------------------------------------------------------------
  end function rate3B
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine get_chmcycles(CHMCYCLES, NRCHEM, NROPSM)
    !// --------------------------------------------------------------------
    !// Internal UiO loop for chemistry/boundary layer mixing
    !// Number of cycles per NOPS when NROPSM==3.
    !//
    !// Ole Amund Sovde, November 2009
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Input
    integer, intent(in) :: NRCHEM, NROPSM
    !// Output
    integer, intent(out) :: CHMCYCLES
    !// --------------------------------------------------------------------
    character(len=*), parameter :: subr = 'get_chmcycles'
    !// --------------------------------------------------------------------

    if (NRCHEM .eq. 1) then
       CHMCYCLES = 4
       !// Simple test to take NROPSM into account,
       if (NROPSM .ge. 6 .and. NROPSM .lt. 12) CHMCYCLES = 2
       if (NROPSM .ge. 12) CHMCYCLES = 1
    else if (NRCHEM .eq. 2) then
       CHMCYCLES = 2
       !// Simple test to take NROPSM into account
       if (NROPSM .ge. 6) CHMCYCLES = 1
    else
       CHMCYCLES = 1
    end if

    !// --------------------------------------------------------------------
  end subroutine get_chmcycles
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine US76_Atmosphere(alt, sigma, delta, theta)
    !// --------------------------------------------------------------------
    ! PURPOSE: Compute the properties of the 1976 standard atmosphere to 86 km.
    ! AUTHOR: Ralph Carmichael, Public Domain Aeronautical Software
    ! NOTE: If alt > 86, the values returned will not be correct, but they will
    !   not be too far removed from the correct values for density.
    !   The reference document does not use the terms pressure and temperature
    !   above 86 km.
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// A R G U M E N T S
    !// --------------------------------------------------------------------
    REAL(r8),INTENT(IN)::  alt    ! geometric altitude, km.
    REAL(r8),INTENT(OUT):: sigma  ! density/sea-level standard density
    REAL(r8),INTENT(OUT):: delta  ! pressure/sea-level standard pressure
    REAL(r8),INTENT(OUT):: theta  ! temperature/sea-level standard temperature
    !// --------------------------------------------------------------------
    !// L O C A L   C O N S T A N T S
    !// --------------------------------------------------------------------
    REAL(r8),PARAMETER:: REARTH = 6369.0_r8 ! radius of the Earth (km)
    REAL(r8),PARAMETER:: GMR = 34.163195_r8 ! hydrostatic constant
    INTEGER,PARAMETER:: NTAB = 8 ! number of entries in the defining tables
    !// --------------------------------------------------------------------
    !// L O C A L   V A R I A B L E S
    !// --------------------------------------------------------------------
    INTEGER :: i,j,k         ! counters
    REAL(r8) :: h            ! geopotential altitude (km)
    REAL(r8) :: tgrad, tbase ! temperature gradient and base temp of this layer
    REAL(r8) :: tlocal       ! local temperature
    REAL(r8) :: deltah       ! height above base of this layer
    !// --------------------------------------------------------------------
    !//  L O C A L   A R R A Y S   ( 1 9 7 6   S T D.  A T M O S P H E R E )
    !// --------------------------------------------------------------------
    REAL(r8),DIMENSION(NTAB),PARAMETER:: htab = &
         (/   0.0_r8,  11.0_r8,   20.0_r8,   32.0_r8,   47.0_r8, &
             51.0_r8,  71.0_r8,   84.852_r8/)
    REAL(r8),DIMENSION(NTAB),PARAMETER:: ttab = &
         (/288.15_r8, 216.65_r8, 216.65_r8, 228.65_r8, 270.65_r8, &
           270.65_r8, 214.65_r8, 186.946_r8/)
    REAL(r8),DIMENSION(NTAB),PARAMETER:: ptab = &
         (/1.0_r8, 2.233611e-1_r8, 5.403295e-2_r8, 8.5666784e-3_r8, &
         1.0945601e-3_r8, 6.6063531e-4_r8, 3.9046834e-5_r8, 3.68501e-6_r8/)
    REAL(r8),DIMENSION(NTAB),PARAMETER:: gtab = &
         (/-6.5_r8, 0.0_r8, 1.0_r8, 2.8_r8, 0.0_r8, -2.8_r8, -2.0_r8, 0.0_r8/)
    !// --------------------------------------------------------------------
    character(len=*), parameter :: subr = 'US76_Atmosphere'
    !// --------------------------------------------------------------------

    ! convert geometric to geopotential altitude
    h = alt * REARTH / (alt + REARTH)

    i = 1
    j = NTAB                    ! setting up for binary search
    DO
       k = (i + j) / 2              ! integer division
       IF (h .lt. htab(k)) THEN
          j = k
       ELSE
          i = k
       END IF
       IF (j .le. i+1) EXIT
    END DO

    tgrad = gtab(i)             ! i will be in 1...NTAB-1
    tbase = ttab(i)
    deltah = h - htab(i)
    tlocal = tbase + tgrad * deltah
    theta = tlocal / ttab(1)      ! temperature ratio

    IF (tgrad .eq. 0._r8) THEN    ! pressure ratio
       delta = ptab(i) * EXP(-GMR * deltah / tbase)
    ELSE
       delta = ptab(i) * (tbase / tlocal)**(GMR / tgrad)
    END IF

    sigma = delta / theta         ! density ratio

    !// --------------------------------------------------------------------
  end subroutine US76_Atmosphere
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  real(r8) function h2o_sat(T)
    !// --------------------------------------------------------------------
    !// Returns saturated value of H2O [molec/cm3] at temperature T.
    !// For T<273.15K: From Marti and Mauersberger, GRL Vol 20, No 5,
    !// 363-366, 1993.
    !// --------------------------------------------------------------------
    use cmn_parameters, only: KBOLTZ, AVOGNR, TK_0C
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Input
    real(r8), intent(in) :: T
    !// Locals
    real(r8) :: AMP0, E, TC
    !// --------------------------------------------------------------------
    AMP0 = KBOLTZ * T ! mb * cm3/molec
    if (T .lt. TK_0C) then
       H2O_SAT = 10._r8**(12.537_r8 - 2663.5_r8/T) / (AMP0 * 100._r8)
    else
       TC = T - TK_0C
       E  = 6.1_r8 * 10._r8**(7.45_r8 * TC / (235._r8 + TC))
       H2O_SAT = 217.e-6_r8 * E / T / 18._r8 * AVOGNR
    end if
    !// --------------------------------------------------------------------
  end function h2o_sat
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine source_e90(BTT,DTSRCE,MP)
    !// --------------------------------------------------------------------
    !// Adds source to E90 tracer.
    !//
    !// Ole Amund Sovde, February 2012
    !// --------------------------------------------------------------------
    use cmn_size, only: LPAR, NPAR, IDBLK, JDBLK
    use cmn_ctm, only: MPBLK, MPBLKJB, MPBLKJE, MPBLKIB, MPBLKIE, &
         AREAXY
    use cmn_chem, only: Ne90
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Input/Output
    integer, intent(in) :: MP
    real(r8),  intent(inout) ::  BTT(LPAR,NPAR,IDBLK,JDBLK)
    real(r8),  intent(in) :: DTSRCE

    !// Locals
    real(r8) :: ZTOTAREA    !// Total area
    real(r8) :: ADDTC       !// Emitted mass
    integer :: I,II,J,JJ

    !// Global emission rate kg/s
    real(r8), parameter :: g_rate = 1.0605e12_r8/31536000_r8
    !// --------------------------------------------------------------------

    !// Check for E90-tracer
    if (Ne90 .le. 0) return

    !// Increase source at surface
    ZTOTAREA = 1._r8 / sum(AREAXY)

    do J = MPBLKJB(MP), MPBLKJE(MP)
       JJ   = J - MPBLKJB(MP) + 1
       do I = MPBLKIB(MP), MPBLKIE(MP)
          II   = I - MPBLKIB(MP) + 1
          ADDTC = g_rate * AREAXY(I,J) * ZTOTAREA * DTSRCE
          BTT(1,Ne90,II,JJ) = BTT(1,Ne90,II,JJ) + ADDTC
       end do !// do I = MPBLKIB(MP), MPBLKIE(MP)
    end do !// do J = MPBLKJB(MP), MPBLKJE(MP)

    !// --------------------------------------------------------------------
  end subroutine source_e90
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine decay_e90 (BTT,BXT,BXX,BYT,BYY,BZT,BZZ,BXY,BXZ,BYZ,DTCHEM,MP)
    !// --------------------------------------------------------------------
    !// Provides a simple e-fold decay of species throughout the model
    !// domain. Applies for the E90 tracer, i.e. uses 90 days e-fold decay.
    !// Take from UCI routine DECAY.
    !// --------------------------------------------------------------------
    use cmn_precision, only: rMom
    use cmn_size, only: LPAR, NPAR, IDBLK, JDBLK
    use cmn_ctm, only: AREAXY
    use cmn_chem, only: Ne90
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Input/output
    real(r8), intent(inout), dimension(LPAR,NPAR,IDBLK,JDBLK):: BTT
    real(rMom), intent(inout), dimension(LPAR,NPAR,IDBLK,JDBLK):: &
         BXT, BXX, BYT, BYY, BZT, BZZ, BXY, BXZ, BYZ
    real(r8), intent(in)::  DTCHEM
    integer, intent(in):: MP

    !// Locals
    integer :: II,JJ,L
    real(r8)  :: F1L
    !// e-fold decay over DECTIM days
    real(r8), parameter :: DECTIM = 90._r8
      !// ------------------------------------------------------------------

    !// Check for E90-tracer
    if (Ne90 .le. 0) return

    if (DECTIM .gt. 0._r8) then
       F1L = exp(-DTCHEM / (86400._r8 * DECTIM))
       do JJ = 1, JDBLK
          do II = 1, IDBLK
             do L = 1, LPAR
                BTT(L,Ne90,II,JJ) = BTT(L,Ne90,II,JJ) * F1L
                BZT(L,Ne90,II,JJ) = BZT(L,Ne90,II,JJ) * F1L
                BZZ(L,Ne90,II,JJ) = BZZ(L,Ne90,II,JJ) * F1L
                BXZ(L,Ne90,II,JJ) = BXZ(L,Ne90,II,JJ) * F1L
                BYZ(L,Ne90,II,JJ) = BYZ(L,Ne90,II,JJ) * F1L
                BXT(L,Ne90,II,JJ) = BXT(L,Ne90,II,JJ) * F1L
                BXX(L,Ne90,II,JJ) = BXX(L,Ne90,II,JJ) * F1L
                BYT(L,Ne90,II,JJ) = BYT(L,Ne90,II,JJ) * F1L
                BYY(L,Ne90,II,JJ) = BYY(L,Ne90,II,JJ) * F1L
                BXY(L,Ne90,II,JJ) = BXY(L,Ne90,II,JJ) * F1L
             end do
          end do
       end do
    end if

    !// --------------------------------------------------------------------
  end subroutine decay_e90
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine init_e90()
    !// --------------------------------------------------------------------
    !// Initialize e90. Set constant in troposphere, assume e-fold above.
    !// If E90 needs to be initialised, LMTROP is first set from PVU.
    !// This is because LMTROP should be calculated after E90 tropopause
    !// is calculated (in physics_oslo.f90), which would fail if E90 was not
    !// set.
    !//
    !// Amund Sovde Haslerud, September 2017
    !//   Added possible initialisation of LMTROP.
    !// Amund Sovde, February 2012
    !// --------------------------------------------------------------------
    use cmn_size, only: IPAR, JPAR, LPAR
    use cmn_ctm, only: STT, AIR
    use cmn_chem, only: Ne90, TMASSMIX2MOLMIX, TMOLMIX2MASSMIX
    use cmn_oslo, only: LMTROP
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Locals
    integer :: I, J, L
    real(r8) :: rtmp
    real(r8), parameter :: init_ppbv = 110.e-9_r8
    !// --------------------------------------------------------------------
    character(len=*), parameter :: subr = 'init_e90'
    !// --------------------------------------------------------------------

    !// Check for E90-tracer
    if (Ne90 .le. 0) return

    !// Is perhaps E90 already set?
    rtmp = maxval(STT(:,:,:,Ne90) / AIR(:,:,:)) * TMASSMIX2MOLMIX(Ne90)
    if (rtmp .gt. 1.e-12_r8) return

    if (maxval(LMTROP) .eq. 0) then
       write(6,'(a)') f90file//':'//subr// &
           ': Why is notLMTROP initialised from PVU before E90 '// &
           'initialisation?'
       stop
    end if

    do J = 1, JPAR
       do I = 1, IPAR
          !// Set 110 ppb up to LMTROP-2
          do L = 1, LMTROP(I,J)-2
             STT(I,J,L,Ne90) = init_ppbv * AIR(I,J,L) * TMOLMIX2MASSMIX(NE90)
          end do
          !// Exponential above
          do L = LMTROP(I,J)-1, LPAR
             STT(I,J,L,Ne90) = init_ppbv * AIR(I,J,L) * TMOLMIX2MASSMIX(NE90) &
                  * exp(-real(L-LMTROP(I,J), r8) / real(LPAR-LMTROP(I,J), r8))
          end do
       end do
    end do

    write(6,'(a)') f90file//':'//subr//': Initialized E90 tracer'

    !// --------------------------------------------------------------------
  end subroutine init_e90
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine tpause_e90()
    !// --------------------------------------------------------------------
    !// Initialize LSTRATAIS_E90, a 3D logical array where stratopsheric
    !// air is T and tropospheric air is F, based on the E90 tracer.
    !//
    !// Amund Sovde Haslerud, August 2017
    !//   Changed names.
    !// Ole Amund Sovde, February 2012
    !// --------------------------------------------------------------------
    use cmn_size, only: IPAR, JPAR, LPAR
    use cmn_ctm, only: STT, AIR
    use cmn_chem, only: Ne90, TMOLMIX2MASSMIX, &
         E90VMR_TP, LSTRATAIR_E90, LPAUZTOP
    use cmn_oslo, only: LMTROP
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Locals
    integer :: I,J,L
    real(r8) :: tpval
    !// --------------------------------------------------------------------
    character(len=*), parameter :: subr = 'tpause_e90'
    !// --------------------------------------------------------------------

    !// Check for E90-tracer
    if (Ne90 .le. 0) return

    !// E90 vmr value for tropopause, converted to mmr
    tpval = E90VMR_TP * TMOLMIX2MASSMIX(Ne90)

    !// Initialize LSTRATAIR_E90
    do L = 1, LPAR
       do J = 1, JPAR
          do I = 1, IPAR
             LSTRATAIR_E90(L,I,J) = STT(I,J,L,Ne90) .lt. tpval*AIR(I,J,L)
          end do
       end do
    end do

    !// Fint max height level of tropospheric air
    do J = 1, JPAR
       do I = 1, IPAR
          do L = 1, LPAR
             if (.not.LSTRATAIR_E90(L,I,J)) LPAUZTOP(I,J) = L
          end do
       end do
    end do

    write(6,'(a)') f90file//':'//subr//': Initialized LSTRATAIR_E90'

    !// --------------------------------------------------------------------
  end subroutine tpause_e90
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine tpauseb_e90(BTT,AIRB,LSTRATAIR_E90B,MP)
    !// --------------------------------------------------------------------
    !// Define the stratosphere based on E90 tracer.
    !// LSTRATAIR_E90(L,I,J) is a 3D logical, .true. = stratospheric
    !// called within the IJ-blocks.
    !// MPSPLIT splits it into private array LSTRATAIR_E90B.
    !//
    !// Amund Sovde Haslerud, August 2017
    !//   Changed names.
    !// Ole Amund Sovde, February 2012
    !// --..----------------------------------------------------------------
    use cmn_size, only: LPAR, NPAR, IDBLK, JDBLK
    use cmn_ctm, only: MPBLKIB, MPBLKIE, MPBLKJB, MPBLKJE
    use cmn_chem, only: Ne90, TMOLMIX2MASSMIX, E90VMR_TP, LPAUZTOP
    use cmn_oslo, only: LMTROP
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Input/output
    real(r8), intent(in) ::  BTT(LPAR,NPAR,IDBLK,JDBLK)
    real(r8), intent(in) ::  AIRB(LPAR,IDBLK,JDBLK)
    integer, intent(in) ::  MP
    logical, intent(out) :: LSTRATAIR_E90B(LPAR,IDBLK,JDBLK)

    !// Locals
    real(r8) :: tpval
    integer :: I, J, L, II, JJ
    !// ..------------------------------------------------------------------
    character(len=*), parameter :: subr = 'toauseb_e90'
    !// --------------------------------------------------------------------

    !// Check for E90-tracer
    if (Ne90 .le. 0) return

    !// E90 vmr value for tropopause, converted to mmr
    tpval = E90VMR_TP * TMOLMIX2MASSMIX(Ne90)

    do J = MPBLKJB(MP), MPBLKJE(MP)
       JJ    = J - MPBLKJB(MP) + 1
       do I = MPBLKIB(MP), MPBLKIE(MP)
          II    = I - MPBLKIB(MP) + 1
          do L = 1, LPAR
             LSTRATAIR_E90B(L,II,JJ) = &
                  BTT(L,Ne90,II,JJ) .lt. tpval * AIRB(L,II,JJ)
             !// Fint max height level of tropospheric air
             if (.not.LSTRATAIR_E90B(L,II,JJ)) LPAUZTOP(I,J) = L
          end do
       end do
    end do

    !// --------------------------------------------------------------------
  end subroutine tpauseb_e90
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine tpauseb_o3(BTT,AIRB,LSTRATAIR_O3,MP)
    !// --------------------------------------------------------------------
    !// Define the stratosphere based on O3 isopleths.
    !// Note that LSTRATAIR_O3(L,I,J) is a 3D logical, .true. = stratospheric
    !// called within the IJ-blocks.
    !//
    !// Amund Sovde Haslerud, August 2017
    !//   Changed names.
    !// Ole Amund Sovde, February 2012
    !// --------------------------------------------------------------------
    use cmn_size, only: IPAR, JPAR, LPAR, NPAR, IDBLK, JDBLK
    use cmn_ctm, only: MPBLKIB, MPBLKIE, MPBLKJB, MPBLKJE
    use cmn_chem, only: TMOLMIX2MASSMIX, O3iso1, O3iso2
    use cmn_met, only: ZOFLE
    use cmn_oslo, only: trsp_idx
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Input/output
    real(r8), intent(in) ::  BTT(LPAR,NPAR,IDBLK,JDBLK)
    real(r8), intent(in) ::  AIRB(LPAR,IDBLK,JDBLK)
    integer, intent(in) ::  MP
    logical, intent(out) :: LSTRATAIR_O3(LPAR,IPAR,JPAR,2)

    !// Locals
    real(r8) :: tpval
    integer :: I, J, L, II, JJ
    !// ------------------------------------------------------------------
    character(len=*), parameter :: subr = 'toauseb_o3'
    !// --------------------------------------------------------------------

    !// Check for O3
    if (trsp_idx(1) .le. 0) return

    !// Select O3 isopleth 1 as chemical tropopause (mmr)
    tpval = O3iso1 * TMOLMIX2MASSMIX(trsp_idx(1))

    do J = MPBLKJB(MP), MPBLKJE(MP)
       JJ    = J - MPBLKJB(MP) + 1
       do I = MPBLKIB(MP), MPBLKIE(MP)
          II    = I - MPBLKIB(MP) + 1
          do L = 1, LPAR
             LSTRATAIR_O3(L,I,J,1) = &
                  (ZOFLE(L,I,J) .gt. 4000._r8) .and. &
                   BTT(L,trsp_idx(1),II,JJ) .ge. tpval*AIRB(L,II,JJ)
          end do
          !// Loop down from top and set possible strat boxes
          do L = LPAR, 1, -1
             if (LSTRATAIR_O3(L,I,J,1)) then
                exit
             else
                LSTRATAIR_O3(L,I,J,1) = .true.
             end if
          end do
          !// Repeat to set whole tropopause to false
          !// Will mess up tropopause folds; limit to tropics!
          !if (abs(YDGRD(J)).lt. 30._r8) then
          !   do L = LPAR, 1, -1
          !      if (.not.LSTRATAIR_O3(L,I,J,1)) exit
          !   end do
          !   LSTRATAIR_O3(1:L,I,J,1) = .false.
          !end if

       end do
    end do

    !// Select O3 isopleth 2 as chemical tropopause
    tpval = O3iso2 * TMOLMIX2MASSMIX(trsp_idx(1))

    do J = MPBLKJB(MP), MPBLKJE(MP)
       JJ    = J - MPBLKJB(MP) + 1
       do I = MPBLKIB(MP), MPBLKIE(MP)
          II    = I - MPBLKIB(MP) + 1
          do L = 1, LPAR
             LSTRATAIR_O3(L,I,J,2) = &
                  (ZOFLE(L,I,J) .gt. 4000._r8) .and. &
                  BTT(L,trsp_idx(1),II,JJ) .ge. tpval*AIRB(L,II,JJ)
          end do
          !// Loop down from top and set possible strat boxes
          do L = LPAR, 1, -1
             if (LSTRATAIR_O3(L,I,J,2)) then
                exit
             else
                LSTRATAIR_O3(L,I,J,2) = .true.
             end if
          end do
          !// Repeat to set whole tropopause to false
          !// Will mess up tropopause folds; limit to tropics!
          !if (abs(YDGRD(J)).lt. 30._r8) then
          !   do L = LPAR, 1, -1
          !      if (.not.LSTRATAIR_O3(L,I,J,2)) exit
          !   end do
          !   LSTRATAIR_O3(1:L,I,J,2) = .false.
          !end if

       end do
    end do

    !// --------------------------------------------------------------------
  end subroutine tpauseb_o3
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine SZA_PN(GMTIME,NDAY,YGRDJ,XGRDI, SZA,PN)
    !// --------------------------------------------------------------------
    !// Calculate solar zenith angle (SZA) and whether it is polar
    !// night (PN) or not.
    !//
    !// Based on Version FJX6.6.
    !// GMTIME = UT for when J-values are wanted
    !//   (for implicit solver this is at the end of the time step)
    !// NDAY   = integer day of the year (used for solar lat and declin)
    !// YGRDJ  = laitude (radians) for grid (I,J)
    !// XGDRI  = longitude (radians) for grid (I,J)
    !//
    !// SZA = solar zenith angle in degrees
    !// COSSZA = U0 = cos(SZA)
    !// PN = polar night (1) not polar night (0)
    !// --------------------------------------------------------------------
    use cmn_parameters, only: CPI
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Input
    real(r8), intent(in) ::   GMTIME,YGRDJ,XGRDI
    integer, intent(in) ::  NDAY
    !// Output
    real(r8), intent(out) ::  SZA
    integer, intent(out)::  PN
    !// Locals
    real(r8) :: PI, PI180, LOCT, minSZA
    real(r8) :: SINDEC, SOLDEK, COSDEC, SINLAT, SOLLAT, COSLAT,COSZ,COSSZA
    integer :: I
    !// --------------------------------------------------------------------

    PI     = CPI
    PI180  = PI/180._r8
    SINDEC = 0.3978_r8 * sin(0.9863_r8 * (real(NDAY, r8) - 80._r8) * PI180)
    SOLDEK = asin(SINDEC)
    COSDEC = cos(SOLDEK)
    SINLAT = sin(YGRDJ)
    SOLLAT = asin(SINLAT)
    COSLAT = cos(SOLLAT)

    !// Find if we have polar night:
    !// Loop through 24 hours and check for minimum SZA
    minSZA = 100._r8
    do I = 1, 24
       LOCT   = ((real(I, r8) * 15._r8) - 180._r8) * PI180 + XGRDI
       COSSZA = COSDEC * COSLAT * cos(LOCT) + SINDEC * SINLAT
       minSZA    = min(minSZA, acos(COSSZA) / PI180)
    end do
    if (minSZA .ge. 90._r8) then
       PN = 1
    else
       PN = 0
    end if

    !// Find SZA
    LOCT   = (((GMTIME) * 15._r8) - 180._r8) * PI180 + XGRDI
    COSSZA = COSDEC * COSLAT * cos(LOCT) + SINDEC * SINLAT
    SZA    = acos(COSSZA) / PI180

    !// --------------------------------------------------------------------
  end subroutine SZA_PN
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine stringUpCase(UCSTR,INSTR)
    !// --------------------------------------------------------------------
    !// Put a string into upper case letters (English alphabet).
    !//
    !// Ole Amund Sovde, August 2009
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    !// Input
    character(len=*),intent(in) :: INSTR
    !// Output
    character(len=*),intent(out) :: UCSTR

    !// Locals
    character(len=*),parameter :: LC='abcdefghijklmnopqrstuvwxyz'
    character(len=*),parameter :: UC='ABCDEFGHIJKLMNOPQRSTUVWXYZ'
    integer :: i, n,strlen
    !// --------------------------------------------------------------------

    UCSTR = INSTR
    strlen = min(len(INSTR),len(UCSTR))

    do i = 1, strlen
       !// Find location of letter in lower case constant string
       n = index( LC, INSTR(i:i))
       !// If current substring is a lower case letter, make it upper case,
       !// otherwise keep the input letter.
       if ( n .ne. 0 ) then
          UCSTR(i:i) = UC(n:n)
       else
          UCSTR(i:i) = INSTR(i:i)
       end if
    end do

    !// --------------------------------------------------------------------
  end subroutine stringUpCase
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine landfrac2mosaic(LFmosaic,nLFmosaic,LFin,nLFin, lat, LANDUSE_IDX)
    !// --------------------------------------------------------------------
    !// Convert land type fractions to EMEP categories, to be used
    !// in dry deposition scheme.
    !// Used two times in drydeposition_oslo.f90 and one time in
    !// bcoc_oslo.f90.
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    integer, intent(in) :: nLFmosaic, nLFin, LANDUSE_IDX
    real(r8), intent(in) :: lat
    real(r8), dimension(nLFin), intent(in) :: LFin
    real(r8), dimension(nLFmosaic), intent(out) :: LFmosaic
    !// --------------------------------------------------------------------
    character(len=*), parameter :: subr = 'landfrac2mosaic'
    !// --------------------------------------------------------------------

    LFmosaic(:) = 0._r8

    !// Will set land categories, and then estimate ocean, since there
    !// may be slight inconsistencies between 1-PLAND and 1 - land
    !// categories. There should not be, but I have seen that
    !// 1 - land categories can be a tiny negative number.

    !// Vegetation type fractions (Fraction Land-use)
    if (LANDUSE_IDX .eq. 2) then
       !// ----------------------------------------------------------------
       !// ISLSCP2 MODIS land fraction and type data
       !// Note: Type 0 represents water and is not included in LS_FRAC.
       !//       It has been used to calculate land fraction PLAND.
       !//  0=Water Bodies                      1=Evergreen Needleleaf Forests
       !//  2=Evergreen Broadleaf Forests       3=Deciduous Needleleaf Forests
       !//  4=Deciduous Broadleaf Forests       5=Mixed Forests
       !//  6=Closed Shrublands                 7=Open Shrublands
       !//  8=Woody Savannas                    9=Savannas
       !// 10=Grasslands                       11=Permanent Wetlands
       !// 12=Croplands                        13=Urban and Built-Up
       !// 14=Cropland/Natural Vegetation Mosaic 15=Permanent Snow and Ice
       !// 16=Barren or Sparsely Vegetated       17=Unclassified
       !// ----------------------------------------------------------------
       !// mOSaic categories (reduced from Simpson et al., 2012)
       !//  1. Needleleaftree (temperated/boreal)
       !//  2. Deciduoustree (temperated/boral)
       !//  3. Needleleaftree (mediterranean)
       !//  4. Broadleaftree (mdeiterranean)
       !//  5. Crops <a,b,c>
       !//  6. Moorland (savanna++)
       !//  7. Grassland
       !//  8. Scrubs (med.)
       !//  9. Wetlands
       !//  10. Tundra
       !//  11. Desert
       !//  12. Water
       !//  13. Urban
       !//  14. Ice/Snow
       !// MODIS -> mOSaic
       !// 1: 1(90-45S, 45-90N),3, 0.5*5
       !// 2: 4, 0.5*5
       !// 3: 1 (45S-45N)
       !// 4: 2
       !// 5: 12,14
       !// 6: - [6,7,8,9]
       !// 7: 10
       !// 8: 6,7,8,9 [-]
       !// 9: 11
       !// 10: 16(90S-60S,60N-90N)
       !// 11: 16(60S-60N), 17
       !// 12: 1.d0 - PLAND(I,J)
       !// 13: 13
       !// 14: 15 (but will be treated seperately under snow/ice)
       !// ----------------------------------------------------------------
       LFmosaic(1) = LFin(3)+0.5_r8*LFin(5)
       LFmosaic(2) = LFin(4)+0.5_r8*LFin(5)
       LFmosaic(3) = 0._r8
       if (lat.lt.-45 .or. lat.gt.45) then
          LFmosaic(1) = LFmosaic(1)+LFin(1)
       else
          LFmosaic(3) = LFmosaic(3)+LFin(1)
       end if
       LFmosaic(4) = LFin(2)
       LFmosaic(5) = LFin(12) + LFin(14)
       !LFmosaic(6) = 0._r8
       LFmosaic(6) = sum(LFin(6:9))
       LFmosaic(7) = LFin(10)
       !LFmosaic(8) = sum(LFin(6:9))
       LFmosaic(8) = 0._r8
       LFmosaic(9) = LFin(11)

       if (lat.lt.-60._r8 .or. lat.gt.60._r8) then
          LFmosaic(10) = LFin(16)
          LFmosaic(11) = 0._r8
       else
          LFmosaic(10) = 0._r8
          LFmosaic(11) = LFin(16)
       end if
       LFmosaic(11) = LFmosaic(11) + LFin(17)
       LFmosaic(13) = LFin(13)
       LFmosaic(14) = LFin(15)
       !// Set Category "ocean" at the end
       LFmosaic(12) = max(0._r8, 1._r8 - (sum(LFmosaic(1:11)) + sum(LFmosaic(13:))))

    else if (LANDUSE_IDX .eq. 3) then

       !// --------------------------------------------------------------------
       !// mOSaic categories (reduced from Simpson et al., 2012)
       !//  1. Needleleaftree (temperated/boreal)
       !//  2. Deciduoustree (temperated/boral)
       !//  3. Needleleaftree (mediterranean)
       !//  4. Broadleaftree (mediterranean)
       !//  5. Crops <a,b,c>
       !//  6. Moorland (savanna++)
       !//  7. Grassland
       !//  8. Scrubs (med.)
       !//  9. Wetlands
       !//  10. Tundra
       !//  11. Desert
       !//  12. Water
       !//  13. Urban
       !//  14. Ice/Snow
       !// CLM categories
       !// CLM  CTM indices
       !//  1    17    Barren land
       !//  2     1    Needleaf evergreen temperate tree
       !//  3     2    Needleaf evergreen boreal tree
       !//  4     3    Needleaf deciduous boreal tree
       !//  5     4    Broadleaf evergreen tropical tree
       !//  6     5    Broadleaf evergreen temperate tree
       !//  7     6    Broadleaf deciduous tropical tree
       !//  8     7    Broadleaf deciduous temperate tree
       !//  9     8    Broadleaf deciduous boreal tree
       !// 10     9    Broadleaf evergreen temperate shrub
       !// 11    10    Broadleaf deciduous temperate shrub
       !// 12    11    Broadleaf deciduous boreal shrub
       !// 13    12    Arctic C3 grass (cold)
       !// 14    13    C3 grass (cool)
       !// 15    14    C4 grass (warm)
       !// 16    15    Crop1
       !// 17    16    Crop2
       !// CLM -> mOSaic
       !// 1: 1 (90-45S, 45-90N) 2,3
       !// 2: 5,7,8
       !// 3: 1 (45S-45N)
       !// 4: 4,6
       !// 5: 15,16
       !// 6: - [14]
       !// 7: 13,14 [12,13]
       !// 8: 9,10
       !// 9: -
       !// 10: 11,12 [11]
       !// 11: 17
       !// 12: Ocean (1-PLAND)
       !// 13: -
       !// 14: -
       !// --------------------------------------------------------------------
       LFmosaic(1) = LFin(2)+ LFin(3)
       LFmosaic(2) = LFin(5) + sum(LFin(7:8))
       LFmosaic(3) = 0._r8
       if ((lat.lt.-45) .or. (lat.gt.45)) then
          LFmosaic(1) = LFmosaic(1)+LFin(1)
       else
          LFmosaic(3) = LFmosaic(3)+LFin(1)
       end if
       LFmosaic(4) = LFin(4) + LFin(6)
       LFmosaic(5) = sum(LFin(15:16))
       !LFmosaic(6) = 0._r8
       LFmosaic(6) = LFin(14)
       !LFmosaic(7) = LFin(13) + LFin(14)
       LFmosaic(7) = sum(LFin(12:13))
       LFmosaic(8) = LFin(9) + LFin(10)
       LFmosaic(9) = 0._r8
       !LFmosaic(10) = LFin(11) + LFin(12)
       LFmosaic(10) = LFin(11)
       if ((lat.gt.-60) .and. (lat.lt.60)) then
          LFmosaic(11) = LFin(17)
          LFmosaic(14) = 0._r8
       else
          LFmosaic(14) = LFin(17)
          LFmosaic(11) = 0._r8
       end if
       LFmosaic(13) = 0._r8

       !// Ocean may not be fully compatible with 1-PLAND:
       LFmosaic(12) = max(0._r8, 1._r8 - (sum(LFmosaic(1:11)) + sum(LFmosaic(13:))))

    else
       write(6,'(a,i5)') f90file//':'//subr// &
            ': unknown LANDUSE_IDX ',LANDUSE_IDX
       stop
    end if

    !// --------------------------------------------------------------------
  end subroutine landfrac2mosaic
  !// ----------------------------------------------------------------------
!// ----------------------------------------------------------------------
  subroutine set_vegetation_height(start_height,lat, out_height)
    !// --------------------------------------------------------------------
    !// MODIFY forest vegetation height dependent on latitude.
    !//
    !// Based on the modification polwards from 60deg North
    !// in Simpson et al. 2012.
    !// Extended southward to the tropics and usage on southern hemesphere.
    !//
    !// Additional ESTIMATION:
    !// Average Tropical forest is about 2-times higher (e.g. 40m) than
    !// midlatitude forest. Decrease height linearly to 20m at 40deg.
    !// Constant height thereafter till 60deg.
    !// Valide for both, northern and southern, hemisphere.
    !//
    !// Not sure if this is realistic, though.
    !//
    !// Used two times in drydeposition_oslo.f90 and one time in
    !// bcoc_oslo.f90.
    !//
    !// Stefanie Falk, March 2018
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    real(r8), intent(in) :: start_height, lat
    real(r8), intent(out) :: out_height
    !// Locals
    real(r8) :: rate
    !// --------------------------------------------------------------------
    character(len=*), parameter :: subr = 'get_vegetation_height'
    !// --------------------------------------------------------------------
    rate = 0.05_r8

    if (abs(lat) .le. 10._r8) then
       out_height = start_height*2._r8
    else if (abs(lat) .le. 40._r8) then
       out_height = (1._r8-(abs(lat)-10._r8)*rate/3._r8)*start_height*2._r8
    else if (abs(lat) .le. 60._r8) then
       !// Forest vegetation height
       out_height = start_height
    else if (abs(lat) .lt. 74._r8) then
       !// From 60-74 latitude, linearly decreasing with rate
       out_height = (1._r8-(abs(lat)-60._r8)*rate)*start_height
    else
       !// High-latitude forest
       out_height = start_height*3._r8/10._r8
    end if

    if (out_height .lt. 0._r8) then
       write(6,'(a,i5)') f90file//':'//subr// &
            ': There is something wrong with your vegetation height!'
       write(6,'(a)') 'Start (m): '
       write(6,*) start_height
       write(6,'(a)') 'Stop (m): '
       write(6,*), out_height
       stop
    end if
    !// --------------------------------------------------------------------
  end subroutine set_vegetation_height
  !// ----------------------------------------------------------------------
  !// ----------------------------------------------------------------------
  subroutine MAPPED_GROWSEASON ( JDAY, I, J, GDAY, GLEN )
    !// --------------------------------------------------------------------
    !  DESCRIPTION
    !     This function uses a preprocessed map of grow season GDAY and GLEN.
    !     The preprocessed data was created using the 5 days/5 degC
    !     criterum to estimate the begin and end of growing season
    !     from the OpenIFS meteorological data of the same year
    !     for latitudes between 45 and 85deg north.
    !     All others are basically the same as computed with GROWSEASON.
    !
    !// Stefanie Falk, August 2018
    !//   Based on GROWSEASON (below).
    !//   Enable usage of preprocessed GDAY/GLEN
    !// --------------------------------------------------------------------
    use cmn_sfc, only: GDAY_MAP, GLEN_MAP
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    integer, intent(in)  :: JDAY
    integer, intent(in)  :: I, J
    integer, intent(out) :: GDAY, GLEN
    !// --------------------------------------------------------------------
    character(len=*), parameter :: subr = 'MAPPED_GROWSEASON'
    !// --------------------------------------------------------------------
    GLEN = GLEN_MAP(I,J)
    GDAY = GDAY_MAP(I,J,JDAY)

  end subroutine MAPPED_GROWSEASON
  !// ----------------------------------------------------------------------

  !// ----------------------------------------------------------------------
  subroutine GROWSEASON ( JDAY, LAT, GDAY, GLEN )
    !// --------------------------------------------------------------------
    !  DESCRIPTION
    !    This function computes the day of the growing season
    !    corresponding to the given day of year.
    !
    !  HISTORY:
    !    07/21/11 : Imported from SMOKE-BEIS v3.14 and modified  (Tan)
    !               Variation of growing season depends on latitude
    !               (Guenther)
    !//
    !// Amund Sovde Haslerud, October 2017
    !//   Based on MEGANv2.10 code, modified for Oslo CTM3
    !// --------------------------------------------------------------------
    use cmn_ctm, only: LYEAR
    !// --------------------------------------------------------------------
    implicit none
    !// --------------------------------------------------------------------
    integer, intent(in)  :: JDAY
    real(r8), intent(in) :: LAT
    integer, intent(out) :: GDAY, GLEN
    !// Locals
    integer  :: &
         GSEASON_START, GSEASON_END, &
         MAY31,NOV1,DEC31
    !// --------------------------------------------------------------------
    character(len=*), parameter :: subr = 'GROWSEASON'
    !// --------------------------------------------------------------------

    !// Days of year for important dates, taking leap year into account
    if (LYEAR) then
       MAY31 = 152
       NOV1  = 306
       DEC31 = 366
    else
       MAY31 = 151
       NOV1  = 305
       DEC31 = 365
    end if

    if ( LAT .le. 23.0_r8 .and. LAT .ge. -23.0_r8 ) then
       ! tropical regions, year round
       GSEASON_START = 1
       GSEASON_END = DEC31

       GDAY = JDAY - GSEASON_START + 1
       GLEN = GSEASON_END - GSEASON_START + 1

    else if ( LAT .lt. -23.0_r8 ) then
       ! southern hemisphere
       if ( LAT .lt. -60.0_r8 ) then
          ! antarctic start = 0 end = 0, no growing
          GDAY = 0
          GLEN = 0
       else
          ! southern hemisphere temperate, NOV, DEC, JAN-MAY
          if (JDAY .ge. NOV1 .AND. JDAY .le. DEC31 ) THEN
             GSEASON_START = NOV1
             GSEASON_END   = DEC31

             GDAY = JDAY - GSEASON_START + 1
             !// Total length of growing season (Nov-May)
             GLEN = DEC31-NOV1 + 1 + MAY31
          else if (JDAY .ge. 1 .and. JDAY .le. MAY31) then
             GSEASON_START = 1
             GSEASON_END   = MAY31
             !// Add Nov+Dec
             GDAY = JDAY - GSEASON_START + 1 + 61
             GLEN = 61 + GSEASON_END - GSEASON_START + 1
             !// Total length of growing season (Nov-May)
             GLEN = DEC31-NOV1 + 1 + MAY31
          else
             GDAY = 0
             GLEN = 0
          end if

       end if

    else if ( LAT .gt. 23._r8 ) then

       ! northern hemisphere
       if ( LAT .gt. 65.0_r8 ) then
          ! arctic start = 0 end = 0, no growing season
          GDAY = 0
          GLEN = 0
       else
          ! northern hemisphere temperate
          ! start= (lat-23)*4.5            189
          ! end = 365 -((lat-23)*3.3)      226

          GSEASON_START = INT( (LAT - 23.0_r8) * 4.5_r8 )
          GSEASON_END   = DEC31 - INT( (LAT - 23.0_r8) * 3.3_r8 )

          if (JDAY .ge. GSEASON_START .and. JDAY .le. GSEASON_END) then
             GDAY = JDAY - GSEASON_START + 1
          else
             GDAY = 0
          end if
          GLEN = GSEASON_END - GSEASON_START + 1
       end if
    else
       write(6,'(a,es16.6)') f90file//':'//subr//': Invalid LAT ',LAT
       stop
    end if

    !// --------------------------------------------------------------------
  end subroutine GROWSEASON
  !// ----------------------------------------------------------------------
  !// ----------------------------------------------------------------------
end module utilities_oslo
!//=========================================================================


text:
program test_basis
use basis, only: get_cartesian_shell, getstart, get_2ints_size
use utils, only: assert
implicit none

integer :: l, m, n
integer, allocatable, dimension(:) :: lpower, mpower, npower
l = 0
allocate(lpower((l+1)*(l+2)/2), mpower((l+1)*(l+2)/2), npower((l+1)*(l+2)/2))
call get_cartesian_shell(l, lpower, mpower, npower)
call assert(all(lpower == [0]))
call assert(all(mpower == [0]))
call assert(all(npower == [0]))
deallocate(lpower, mpower, npower)

l = 1
allocate(lpower((l+1)*(l+2)/2), mpower((l+1)*(l+2)/2), npower((l+1)*(l+2)/2))
call get_cartesian_shell(l, lpower, mpower, npower)
call assert(all(lpower == [1, 0, 0]))
call assert(all(mpower == [0, 1, 0]))
call assert(all(npower == [0, 0, 1]))
deallocate(lpower, mpower, npower)

l = 2
allocate(lpower((l+1)*(l+2)/2), mpower((l+1)*(l+2)/2), npower((l+1)*(l+2)/2))
call get_cartesian_shell(l, lpower, mpower, npower)
call assert(all(lpower == [2, 1, 1, 0, 0, 0]))
call assert(all(mpower == [0, 1, 0, 2, 1, 0]))
call assert(all(npower == [0, 0, 1, 0, 1, 2]))
deallocate(lpower, mpower, npower)

l = 3
allocate(lpower((l+1)*(l+2)/2), mpower((l+1)*(l+2)/2), npower((l+1)*(l+2)/2))
call get_cartesian_shell(l, lpower, mpower, npower)
call assert(all(lpower == [3, 2, 2, 1, 1, 1, 0, 0, 0, 0]))
call assert(all(mpower == [0, 1, 0, 2, 1, 0, 3, 2, 1, 0]))
call assert(all(npower == [0, 0, 1, 0, 1, 2, 0, 1, 2, 3]))
deallocate(lpower, mpower, npower)

l = 4
allocate(lpower((l+1)*(l+2)/2), mpower((l+1)*(l+2)/2), npower((l+1)*(l+2)/2))
call get_cartesian_shell(l, lpower, mpower, npower)
call assert(all(lpower == [4, 3, 3, 2, 2, 2, 1, 1, 1, 1, 0, 0, 0, 0, 0]))
call assert(all(mpower == [0, 1, 0, 2, 1, 0, 3, 2, 1, 0, 4, 3, 2, 1, 0]))
call assert(all(npower == [0, 0, 1, 0, 1, 2, 0, 1, 2, 3, 0, 1, 2, 3, 4]))
deallocate(lpower, mpower, npower)

! -------------------------------------------------------
call assert(all(getstart([6, 3, 3, 1, 1, 1]) == [1, 7, 10, 13, 14, 15]))
call assert(all(getstart([6, 6, 6, 3, 3, 1, 1, 1]) == &
    [1, 7, 13, 19, 22, 25, 26, 27]))

! -------------------------------------------------------
n = 15
m = n*16/2
call assert(get_2ints_size(n) == m*(m+1)/2)

end program

max_stars_repo_name: certik/hfsolver

max_issues_repo_name: certik/hfsolver

max_forks_repo_name: certik/hfsolver


text:
      subroutine chksto ( flux  , subs  , stoch , nstoc ,
     j                    itemid, nitem )

      integer      nstoc , nitem
      character*10 flux(nstoc),subs(nstoc),itemid(nitem)
      real         stoch(nstoc)

c     Subroutine to check tables R6-R7-R8

c     Check if the same flux (velocity, dispersion)
c     if it occurs more than once, has a consistent set of
c     stoichiometry rules

c     After the check, the doubles are removed

c     flux   array of items (fluxes, velocity, dispersion)
c     subs   array of affected substances
c     stoch  stoichionetry constants
c     nstoc  nr of stochi lines

c     fluxi  local flux index array
c     subsi  local substance index array
c     repet  indicates if a line is a repeated line

      integer      istoc  , iitem , istoc2, iitem2, isubs , irepet,
     j             nstoc2
      integer      fluxi(nstoc), subsi(nstoc), repet(nstoc)
      logical      done
      integer      isflux(nitem), isrepe(nitem)
      real         effect(nitem), effec2(nitem)

c     zero local arrays

      write ( 11, * ) 'chksto'
      do iitem = 1,nitem
          isflux(iitem) = 0
          isrepe(iitem) = 0
      enddo

c     index table on items

      do istoc = 1,nstoc
          repet(istoc) = 1
          call zoek (flux(istoc),nitem,itemid,10,iitem)
          if (iitem.le.0) then
              write ( 11, * ) flux(istoc)
              stop 'BUG CHKSTO 001'
          endif
          fluxi(istoc) = iitem
          isflux(iitem) = 1
          call zoek (subs(istoc),nitem,itemid,10,iitem)
          if (iitem.le.0) stop 'BUG CHKSTO 002'
          subsi(istoc) = iitem
      enddo
      write ( 11, * ) 'stochi table indexed'

c     check occurence of repeated lines

      do istoc = 1,nstoc
c         check table processed so far, upward
          do istoc2 = istoc-1,1,-1
              if ( fluxi(istoc2) .eq. fluxi(istoc) .and.
     j             subsi(istoc2) .eq. subsi(istoc) ) then
                  repet(istoc) = repet(istoc2) + 1
                  isrepe(fluxi(istoc)) = repet(istoc)
                  goto 10
              endif
          enddo
   10     continue
      enddo

c     check if multiple occurences are consistent
c     if a flux occurs in more than one process, it should affect
c     the same substances in the same way

      do iitem = 1,nitem
c         Only fluxes:
          if ( isflux(iitem) .eq. 1 ) then
              if ( isrepe(iitem) .gt. 1 )
     j        write (11,*) ' flux ',itemid(iitem),isrepe(iitem)
c             Zero effect on substances
              do iitem2 = 1,nitem
                  effect(iitem2) = 0.0
              enddo
c             Compute effect on substances first occurrence
              do istoc = 1,nstoc
                  if ( fluxi(istoc) .eq. iitem ) then
                      if ( repet(istoc) .eq. 1 ) then
                          isubs = subsi(istoc)
                          effect(isubs) =  stoch(istoc)
                      endif
                  endif
              enddo
              do irepet = 2,isrepe(iitem)
              write (11,*) ' check ',irepet
c             Zero effect on substances
              do iitem2 = 1,nitem
                  effec2(iitem2) = 0.0
              enddo
c             Compute effect on substances following occurences
              do istoc = 1,nstoc
                  if ( fluxi(istoc) .eq. iitem ) then
                      if ( repet(istoc) .eq. irepet ) then
                          isubs = subsi(istoc)
                          effec2(isubs) =  stoch(istoc)
                      endif
                  endif
              enddo
c             Check
              do iitem2 = 1,nitem
                  if ( effec2(iitem2) .ne. effect(iitem2) ) then
                      write (*,*) ' Flux/vel/disp ',itemid(iitem),
     j               ' not consistently defined'
                      write (*,*) ' Effect on substance ',
     j                itemid(iitem2),effec2(iitem2),effect(iitem2)
                      stop
                  endif
              enddo
              enddo
          endif
      enddo

c     Clear multiple occurences

  190 nstoc2 = nstoc
      do istoc = 1,nstoc2
          if ( repet(istoc) .gt. 1 ) then
              do istoc2 = istoc+1,nstoc2
                  flux (istoc2-1) = flux (istoc2)
                  fluxi(istoc2-1) = fluxi(istoc2)
                  subs (istoc2-1) = subs (istoc2)
                  subsi(istoc2-1) = subsi(istoc2)
                  repet(istoc2-1) = repet(istoc2)
                  stoch(istoc2-1) = stoch(istoc2)
              enddo
              nstoc = nstoc - 1
              goto 190
          endif
      enddo

      return
      end

ext: f

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
*----------------------------------------------------------------------*
      subroutine occ_contr(occ,ierr,contr,occ_vtx_in,njoined)
*----------------------------------------------------------------------*
*     get resulting occupation of contr
*     occ_vtx is obtained from occvtx4contr with mode==1
*----------------------------------------------------------------------*
      implicit none

      include 'opdim.h'
      include 'stdunit.h'
      include 'def_contraction.h'
      include 'ifc_operators.h'

      integer, parameter ::
     &     ntest = 00

      type(contraction), intent(in), target ::
     &     contr
      integer, intent(in) ::
     &     occ_vtx_in(ngastp,2,contr%nvtx), njoined
      integer, intent(out) ::
     &     ierr,
     &     occ(ngastp,2,njoined)

      integer ::
     &     nvtx, narc, nins, idx1, idx2, iarc, type, type_last
      type(cntr_arc), pointer ::
     &     arc(:)
      integer, pointer ::
     &     occ_vtx(:,:,:), occ_cnt(:,:)

      if (ntest.ge.100)
     &     call write_title(lulog,wst_dbg_subr,'occ_contr')
      call quit(1,'occ_contr','obsolete, use occ_contr2')

      nvtx = contr%nvtx

      if (ntest.ge.100) then
        write(lulog,*) 'occupations on input:'
        call wrt_occ_n(lulog,occ_vtx_in,nvtx)
      end if

      ierr = 0

      ! get a copy of occ_vtx
      allocate(occ_vtx(ngastp,2,nvtx))
      occ_vtx = occ_vtx_in

      narc = contr%narc
      arc => contr%arc
      ! remove all contracted indices
      do iarc = 1, narc
        idx1 = arc(iarc)%link(1)
        idx2 = arc(iarc)%link(2)
        occ_cnt => arc(iarc)%occ_cnt
        occ_vtx(1:ngastp,1:2,idx1) = occ_vtx(1:ngastp,1:2,idx1)
     &       - occ_cnt
        occ_vtx(1:ngastp,1:2,idx2) = occ_vtx(1:ngastp,1:2,idx2)
     &       - iocc_dagger(occ_cnt)
      end do

      if (ntest.ge.100) then
        write(lulog,*) 'remainder occupations:'
        call wrt_occ_n(lulog,occ_vtx,nvtx)
      end if

      ! loop over remainder-occupations
      type_last = -1
      idx2 = 0
      do idx1 = 1, nvtx
        ! get type of remainder: 0, ex, dx, general
        if (iocc_zero(occ_vtx(1:ngastp,1:2,idx1))) then
          type = 0
        else if (iocc_zero(iocc_xdn(2,occ_vtx(1:ngastp,1:2,idx1)))) then
          type = 1
        else if (iocc_zero(iocc_xdn(1,occ_vtx(1:ngastp,1:2,idx1)))) then
          type = 2
        else
          type = 3
        end if
        !  0: remove -- do nothing
        !  ex/dx: if previous vertex was of same type, add up
        if (type.eq.1.or.type.eq.2) then
          if (type.eq.type_last) then
            occ_vtx(1:ngastp,1:2,idx2) = occ_vtx(1:ngastp,1:2,idx2) +
     &                                   occ_vtx(1:ngastp,1:2,idx1)
          else
            idx2 = idx2+1
            if (idx2.lt.idx1)
     &           occ_vtx(1:ngastp,1:2,idx2) = occ_vtx(1:ngastp,1:2,idx1)
          end if
          type_last = type
        else if (type.eq.3) then
          !  general: remains as is
          idx2 = idx2+1
          if (idx2.lt.idx1)
     &         occ_vtx(1:ngastp,1:2,idx2) = occ_vtx(1:ngastp,1:2,idx1)
          type_last = type
        end if
      end do

      if (ntest.ge.100) then
        write(lulog,*) 'after joining vertices (remaining:',idx2,')'
        call wrt_occ_n(lulog,occ_vtx,idx2)
      end if

      ! exception: density storage of (pure DX|pure EX)
      if (idx2.eq.2.and.njoined.eq.1) then
        if (iocc_nonzero(iocc_xdn(1,occ_vtx(1:ngastp,1:2,1))).or.
     &      iocc_nonzero(iocc_xdn(2,occ_vtx(1:ngastp,1:2,2)))) then
          call wrt_occ_n(lulog,occ_vtx,2)
          call quit(1,'occ_contr',
     &         'cannot express supervertex as density')
        end if
        occ(1:ngastp,1:2,1) = occ_vtx(1:ngastp,1:2,1) +
     &                    occ_vtx(1:ngastp,1:2,2)
      else if (idx2.gt.njoined) then
        ierr = -1-idx2
c        write(lulog,*) 'expected super-vertices: ',njoined
c        write(lulog,*) 'but found the following final occupation:'
c        call wrt_occ_n(lulog,occ_vtx,idx2)
c        call quit(1,'occ_contr',
c     &         'resulting supervertex is incompatible with expected
c     &       number of joined vertices')
      else if (idx2.lt.njoined) then
        ! try to fix: insert 0-occ after DX, before EX
c dbg
c        print *,'fixing 0-occ'
c dbg
        nvtx = idx2
        nins = njoined-nvtx
        idx1 = 1
        do idx2 = 1, nvtx
          ! no de-excitation follows?
          if (iocc_zero(iocc_xdn(2,occ_vtx(1:ngastp,1:2,idx2)))) then
            occ(1:ngastp,1:2,idx1:idx1-1+nins) = 0
            occ(1:ngastp,1:2,idx1+nins) = occ_vtx(1:ngastp,1:2,idx2)
            idx1 = idx1+nins+1
          else
            occ(1:ngastp,1:2,idx1) = occ_vtx(1:ngastp,1:2,idx2)
            idx1 = idx1+1
          end if
        end do
        if (idx1.le.njoined) then
          occ(1:ngastp,1:2,idx1:njoined) = 0
        end if
      else
        occ(1:ngastp,1:2,1:njoined) = occ_vtx(1:ngastp,1:2,1:njoined)
      end if

      deallocate(occ_vtx)

      return
      end

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public


text:
      SUBROUTINE POLIMD
C
C     This FORTRAN program produces a polarization vector
C     plot from an image of the polarization and the
C     polarization position angle on a requested graphics
C     device. Optionally the polarization and position
C     angle errors may be supplied to overplot the error
C     extent. Various labelling options may also be
C     exercized.
C
C     The linear polarization and polarization position angle
C     files are read. The linear polarization and
C     polarization position angle error files can also
C     be read. The binning of the data is required to plot on
C     the correct grid (such as for overlaying on an image).
C     A global shift to the position angles can be specified.
C
C     Written by  J. R. Walsh, ST-ECF, ESO (jwalsh@eso.org)  December 1999
C     Version 1.00
C
      IMPLICIT NONE

      INTEGER STAT
      INTEGER STAT1
      INTEGER STAT2
      INTEGER STAT3
      INTEGER STATP

      INTEGER IMDSCR1
      INTEGER IMDSCR2
      INTEGER IMDSCR3
      INTEGER IMDSCR4

      INTEGER NAXISI
      INTEGER DIMEN(7)
      INTEGER DTYPE

      INTEGER P1
      INTEGER P2
      INTEGER Q1
      INTEGER Q2
      INTEGER PE1
      INTEGER PE2
      INTEGER QE1
      INTEGER QE2
      INTEGER NELEM

      INTEGER LPOL
      INTEGER POPA
      INTEGER LPOLER
      INTEGER POPAER

      INTEGER XYBIN
      INTEGER XYSAM

      INTEGER GRDSCR

      REAL POSCAL
      REAL POLOW
      REAL POHIGH
      REAL PAOFF
      REAL VPXY(4)

      REAL MEMR(1)

      CHARACTER*32 INPOL
      CHARACTER*32 INPA
      CHARACTER*32 INPOER
      CHARACTER*32 INPAER

      CHARACTER*16 DEVICE

      CHARACTER*16 XLAB
      CHARACTER*16 YLAB
      CHARACTER*64 LABEL

      CHARACTER*132 OUTEXT

      LOGICAL LOER
      LOGICAL LLAB
      LOGICAL LERASE

      COMMON/MEM/MEMR
C
C     Get name of the linear polarization data set name
C
100   CALL UCLGST('inpol',INPOL,STAT)
      CALL UIMOPN(INPOL,1,IMDSCR1,STAT1)
      IF (STAT1.NE.0) THEN
        WRITE(OUTEXT,105) INPOL
105     FORMAT(' Error opening data file ',A32)
        CALL UMSPUT(OUTEXT,3,0,STAT)
        GO TO 990
      ENDIF
      CALL UIMGID(IMDSCR1,DTYPE,NAXISI,DIMEN,STAT2)
      IF (STAT2.NE.0) THEN
        WRITE(OUTEXT,107) INPOL
107     FORMAT(' Error reading data file header ',A32)
        CALL UMSPUT(OUTEXT,3,0,STAT)
        GO TO 990
      ENDIF
      IF (NAXISI.NE.2) THEN
        WRITE(OUTEXT,109)
109     FORMAT(' Data file is not 2-dimensional')
        CALL UMSPUT(OUTEXT,3,0,STAT)
        GO TO 990
      ENDIF
      P1=DIMEN(1)
      P2=DIMEN(2)
C
C     Allocate dynamic memory for the 2-D input image
C     of dimensions P1 by P2 (real number)
C
      NELEM=P1*P2
      CALL UDMGET(NELEM,6,LPOL,STAT)
      IF (STAT.NE.0) THEN
        WRITE(OUTEXT,111)
111     FORMAT(' Unable to assign memory for internal 2-D array')
        CALL UMSPUT(OUTEXT,3,0,STAT)
        GO TO 990
      ENDIF
C
C     Read the input image to LPOL
C
      CALL UIGS2R(IMDSCR1,1,P1,1,P2,MEMR(LPOL),STAT3)
      IF (STAT3.NE.0) THEN
        WRITE(OUTEXT,113) INPOL
113     FORMAT(' Error reading image into internal 2-D array ',A32)
        CALL UMSPUT(OUTEXT,3,0,STAT)
        GO TO 990
      ENDIF
C
C     Get name of the polarization position angle data set
C     name
C
120   CALL UCLGST('inpopa',INPA,STAT)
      CALL UIMOPN(INPA,1,IMDSCR2,STAT1)
      IF (STAT1.NE.0) THEN
        WRITE(OUTEXT,105) INPA
        CALL UMSPUT(OUTEXT,3,0,STAT)
        GO TO 990
      ENDIF
      CALL UIMGID(IMDSCR2,DTYPE,NAXISI,DIMEN,STAT2)
      IF (STAT2.NE.0) THEN
        WRITE(OUTEXT,107) INPA
        CALL UMSPUT(OUTEXT,3,0,STAT)
        GO TO 990
      ENDIF
      IF (NAXISI.NE.2) THEN
        WRITE(OUTEXT,109)
        CALL UMSPUT(OUTEXT,3,0,STAT)
        GO TO 990
      ENDIF
      Q1=DIMEN(1)
      Q2=DIMEN(2)
C
C     Check input dimensions of polarization position angle
C     image same as polarization image
C
      IF (Q1.NE.P1.OR.Q2.NE.P2) THEN
        WRITE(OUTEXT,125)
125     FORMAT(' Dimensions of input intensity and polarization
     : position angle images not equal')
        CALL UMSPUT(OUTEXT,3,0,STAT)
        GO TO 990
      ENDIF
C
C     Allocate dynamic memory for the 2-D input image
C     of dimensions Q1 by Q2 (real number)
C
      NELEM=Q1*Q2
      CALL UDMGET(NELEM,6,POPA,STAT)
      IF (STAT.NE.0) THEN
        WRITE(OUTEXT,111)
        CALL UMSPUT(OUTEXT,3,0,STAT)
        GO TO 990
      ENDIF
C
C     Read the input image to POPA
C
      CALL UIGS2R(IMDSCR2,1,Q1,1,Q2,MEMR(POPA),STAT3)
      IF (STAT3.NE.0) THEN
        WRITE(OUTEXT,113) INPA
        CALL UMSPUT(OUTEXT,3,0,STAT)
        GO TO 990
      ENDIF
C
C     Close the input polarization and position angle images
C
      CALL UIMCLO(IMDSCR1,STAT)
      CALL UIMCLO(IMDSCR2,STAT)
C
C     Get the switch to determine if imagers of the
C     polarization error and position angle error are
C     to be input
C
      CALL UCLGSB('poi_err',LOER,STAT)
      IF (.NOT.LOER) THEN
C
C       Need to allocate dynamic memory for input arrays of
C       errors even if not used
C
        PE1=P1
        PE2=P2
        QE1=Q1
        QE2=Q2
      ENDIF
      IF (LOER) THEN
C
C       Get name of the linear polarization error name
C
200     CALL UCLGST('inpoer',INPOER,STAT)
        CALL UIMOPN(INPOER,1,IMDSCR3,STAT1)
        IF (STAT1.NE.0) THEN
          WRITE(OUTEXT,105) INPOER
          CALL UMSPUT(OUTEXT,3,0,STAT)
          GO TO 990
        ENDIF
        CALL UIMGID(IMDSCR3,DTYPE,NAXISI,DIMEN,STAT2)
        IF (STAT2.NE.0) THEN
          WRITE(OUTEXT,107) INPOER
          CALL UMSPUT(OUTEXT,3,0,STAT)
          GO TO 990
        ENDIF
        IF (NAXISI.NE.2) THEN
          WRITE(OUTEXT,109)
          CALL UMSPUT(OUTEXT,3,0,STAT)
          GO TO 990
        ENDIF
        PE1=DIMEN(1)
        PE2=DIMEN(2)
      ENDIF
C
C     Allocate dynamic memory for the 2-D input image
C     of dimensions PE1 by PE2 (real number)
C
      NELEM=PE1*PE2
      CALL UDMGET(NELEM,6,LPOLER,STAT)
      IF (STAT.NE.0) THEN
        WRITE(OUTEXT,111)
        CALL UMSPUT(OUTEXT,3,0,STAT)
        GO TO 990
      ENDIF
C
C     Read the input image to LPOLER
C
      IF (LOER) THEN
        CALL UIGS2R(IMDSCR3,1,PE1,1,PE2,MEMR(LPOLER),STAT3)
        IF (STAT3.NE.0) THEN
          WRITE(OUTEXT,113) INPOL
          CALL UMSPUT(OUTEXT,3,0,STAT)
          GO TO 990
        ENDIF
      ENDIF
C
C     Get name of the polarization position angle error
C     name
C
220   IF (LOER) THEN
        CALL UCLGST('inpaer',INPAER,STAT)
        CALL UIMOPN(INPAER,1,IMDSCR4,STAT1)
        IF (STAT1.NE.0) THEN
          WRITE(OUTEXT,105) INPAER
          CALL UMSPUT(OUTEXT,3,0,STAT)
          GO TO 990
        ENDIF
        CALL UIMGID(IMDSCR4,DTYPE,NAXISI,DIMEN,STAT2)
        IF (STAT2.NE.0) THEN
          WRITE(OUTEXT,107) INPAER
          CALL UMSPUT(OUTEXT,3,0,STAT)
          GO TO 990
        ENDIF
        IF (NAXISI.NE.2) THEN
          WRITE(OUTEXT,109)
          CALL UMSPUT(OUTEXT,3,0,STAT)
          GO TO 990
        ENDIF
        QE1=DIMEN(1)
        QE2=DIMEN(2)
C
C       Check input dimensions of polarization position angle
C       error image same as polarization error image
C
        IF (QE1.NE.PE1.OR.QE2.NE.PE2) THEN
          WRITE(OUTEXT,125)
          CALL UMSPUT(OUTEXT,3,0,STAT)
          GO TO 990
        ENDIF
      ENDIF
C
C     Allocate dynamic memory for the 2-D input image
C     of dimensions QE1 by QE2 (real number)
C
      NELEM=QE1*QE2
      CALL UDMGET(NELEM,6,POPAER,STAT)
      IF (STAT.NE.0) THEN
        WRITE(OUTEXT,111)
        CALL UMSPUT(OUTEXT,3,0,STAT)
        GO TO 990
      ENDIF
C
C     Read the input image to POPAER
C
      IF (LOER) THEN
        CALL UIGS2R(IMDSCR4,1,Q1,1,Q2,MEMR(POPAER),STAT3)
        IF (STAT3.NE.0) THEN
          WRITE(OUTEXT,113) INPAER
          CALL UMSPUT(OUTEXT,3,0,STAT)
          GO TO 990
        ENDIF
      ENDIF
C
C     Close the input polarization and position angle error
C     images
C
      IF (LOER) THEN
        CALL UIMCLO(IMDSCR3,STAT)
        CALL UIMCLO(IMDSCR4,STAT)
      ENDIF
C
C     Get the pixel binning factor (same in X and Y) used to
C     form the polarization and position angle images
C
      CALL UCLGSI('xybin',XYBIN,STAT2)
      IF (STAT2.NE.0) THEN
        WRITE(OUTEXT,311)
311     FORMAT(' Failed to get binning factor')
        CALL UMSPUT(OUTEXT,3,0,STAT)
        GO TO 990
      ENDIF
C
C     Get the scalar to multiply the polarization values
C     to produce useful vector lengths
C
      CALL UCLGSR('poscal',POSCAL,STAT2)
      IF (STAT2.NE.0) THEN
        WRITE(OUTEXT,315)
315     FORMAT(' Failed to get scalar value')
        CALL UMSPUT(OUTEXT,3,0,STAT)
        GO TO 990
      ENDIF
C
C     Get the low cut-off in polarization below which not to
C     plot the vectors
C
      CALL UCLGSR('polow',POLOW,STAT2)
      IF (STAT2.NE.0) THEN
        WRITE(OUTEXT,317)
317     FORMAT(' Failed to get cut-off polarization')
        CALL UMSPUT(OUTEXT,3,0,STAT)
        GO TO 990
      ENDIF
C
C     Get the upper cut-off in polarization above which not to
C     plot the vectors
C
      CALL UCLGSR('pohigh',POHIGH,STAT2)
      IF (STAT2.NE.0) THEN
        WRITE(OUTEXT,319)
319     FORMAT(' Failed to get cut-off polarization')
        CALL UMSPUT(OUTEXT,3,0,STAT)
        GO TO 990
      ENDIF
C
C     Get the value of the offset angle to add to all
C     the position angles
C
      CALL UCLGSR('pacorr',PAOFF,STAT2)
      IF (STAT2.NE.0) THEN
        WRITE(OUTEXT,321)
321     FORMAT(' Failed to get position angle offset')
        CALL UMSPUT(OUTEXT,3,0,STAT)
        GO TO 990
      ENDIF
C
C     Get the number of values to sample in X and Y when
C     plotting for clarity of presentation (expressed as
C     number of values to skip in parameter file)
C
      CALL UCLGSI('xysam',XYSAM,STAT2)
      IF (STAT2.NE.0) THEN
        WRITE(OUTEXT,323)
323     FORMAT(' Failed to get spatial sampling step')
        CALL UMSPUT(OUTEXT,3,0,STAT)
        GO TO 990
      ENDIF
      XYSAM=XYSAM+1
C
C     Get the switch to determine if the vector map is
C     to be fully labelled
C
      CALL UCLGSB('polab',LLAB,STAT)
C
C     Get a label, an X-label and a Y-label for the plot
C     if applicable
C
      IF (LLAB) THEN
        CALL UCLGST('title',LABEL,STAT)
        CALL UCLGST('xlabel',XLAB,STAT)
        CALL UCLGST('ylabel',YLAB,STAT)
      ENDIF
C
C     Get the name of the graphics device
C
      CALL UCLGST('device',DEVICE,STAT)
C
C     Get the 4 values of the viewport X minimum
C     and maximum, Y minumum and maximum values for
C     rescaling the plot
C
       CALL UCLGSR('left',VPXY(1),STAT2)
       CALL UCLGSR('right',VPXY(2),STAT2)
       CALL UCLGSR('bottom',VPXY(3),STAT2)
       CALL UCLGSR('top',VPXY(4),STAT2)

      IF (STAT2.NE.0) THEN
        WRITE(OUTEXT,333)
333     FORMAT(' Failed to get array of viewport limits')
        CALL UMSPUT(OUTEXT,3,0,STAT)
        GO TO 990
      ENDIF
C
C     Get the switch to determine whether to erase the
C     graphics device
C
      CALL UCLGSB('erase',LERASE,STAT)
      IF (DEVICE(:7).EQ.'stdplot') THEN
        LERASE=.FALSE.
      ENDIF
C
C     Attempt to open the graphics station. If noerase
C     then open in append mode
C
400   IF (LERASE) THEN
        CALL UGOPEN(5,DEVICE,GRDSCR,STAT)
      ENDIF
      IF (.NOT.LERASE) THEN
        CALL UGOPEN(4,DEVICE,GRDSCR,STAT)
      ENDIF
      IF (STAT.NE.0) THEN
        CALL UMSPUT(' Failed to open graphics station',3,0,STAT)
        GO TO 990
      ENDIF
C
C     Clear the screen if required
C
      IF (LERASE) THEN
        CALL UGCLRS(GRDSCR,STAT)
        IF (STAT.NE.0) THEN
          CALL UMSPUT(' Error clearing graphics station',3,0,STAT)
          GO TO 990
        ENDIF
      ENDIF
C
C     Call the subroutine to plot the polarization vector map
C     on the graphics station
C
      CALL PLOTPOPA(P1,P2,MEMR(LPOL),MEMR(POPA),LOER,
     :              MEMR(LPOLER),MEMR(POPAER),GRDSCR,
     :              XYBIN,POSCAL,POLOW,POHIGH,PAOFF,
     :              XYSAM,LLAB,XLAB,YLAB,LABEL,VPXY,STATP)
C
C     Free the dynamic memory allocated which is no longer
C     required
C
      CALL UDMFRE(LPOL,6,STAT)
      CALL UDMFRE(POPA,6,STAT)
      CALL UDMFRE(LPOLER,6,STAT)
      CALL UDMFRE(POPAER,6,STAT)
C
C     Close the graphics device
C
      CALL UGCLOS(GRDSCR,STAT)

      IF (STATP.EQ.0) THEN
        GO TO 999
      ENDIF

990   WRITE(OUTEXT,991)
991   FORMAT(' Program failed. No vector plot')
      CALL UMSPUT(OUTEXT,1,0,STAT)

999   END

      SUBROUTINE PLOTPOPA(N1,N2,POL,PA,LOER,POLER,PAER,
     :                    GRDSCR,XYBIN,POSCAL,POLOW,POHIGH,
     :                    PAOFF,XYSAM,LLAB,XLAB,YLAB,LABEL,
     :                    VPXY,STATO)
C
C     This subroutine plots using the F77VOS graphics routines
C     the polarization vector field given by the array of
C     polarization POL and position angle PA. If LOER is TRUE
C     then the vectors formed with the polarization error and
C     the position angle error included are overplotted in
C     dashed lines. If the data were binned (XYBIN) to form
C     the polarization then the original grid of the data is
C     reproduced. The polarization vector length is scaled
C     by POSCAL. Vectors are plotted for all polarization
C     values greater than POLOW. A global offset to the
C     position angle is given by PAOFF. Every XYSAM
C     polarization vector can be plotted to reduce crowding.
C     Full labelling is specified by POLAB TRUE
C
C     Input parameters:
C
C       N1 X dimension of input arrays
C       N2 Y dimension of input arrays
C       POL(N1,N2) input array of linear polarization (%)
C       PA(N1,N2) input array of position angle (degrees)
C       LOER logical flag to indicate whether to plot error vectors
C       POLER(N1,N2) input array of linear polarization error (%)
C       PAER(N1,N2) input array of position angle error (degrees)
C       GRDSCR graphics device descriptor
C       XYBIN integer binning factor for X and Y used for polarization images
C       POSCAL scalar to apply to polarization values to plot as pixel lengths
C       POLOW low limit of polarization to plot
C       POHIGH high limit of polarization to plot
C       PAOFF additive correction to apply to PA
C       XYSAM factor to determine plotting every XYSAM vector (in X and Y)
C       LLAB logical flag to indicate whether to label and annotate plot
C       XLAB label for the X axis
C       YLAB label for the Y axis
C       LABEL text to label the plot
C       VWXY(4) array of the viewport limits X min and max, Y min and max
C       STATO status return - returned non-zero if plotting problem occurred
C
      IMPLICIT NONE
      INTEGER N1
      INTEGER N2
      INTEGER GRDSCR
      INTEGER XYBIN
      INTEGER XYSAM
      INTEGER STATO

      REAL POL(N1,N2)
      REAL PA(N1,N2)
      REAL POLER(N1,N2)
      REAL PAER(N1,N2)
      REAL POSCAL
      REAL POLOW
      REAL POHIGH
      REAL PAOFF
      REAL VPXY(4)

      CHARACTER*16 XLAB
      CHARACTER*16 YLAB
      CHARACTER*64 LABEL

      LOGICAL LOER
      LOGICAL LLAB
C
C     Local variables
C
      INTEGER I,J
      INTEGER STAT
      INTEGER STAT2
      INTEGER NPT
      INTEGER LINE
      INTEGER GPLTYP

      REAL TYPSIZ
      REAL RESIZ
      REAL VPMINX
      REAL VPMAXX
      REAL VPMINY
      REAL VPMAXY
      REAL WMINX
      REAL WMAXX
      REAL WMINY
      REAL WMAXY
      REAL ANG
      REAL XLFRAC
      REAL YLFRAC

      REAL SINRD
      REAL COSRD

      REAL POSX(2)
      REAL POSY(2)

      CHARACTER*3 PCLAB
      CHARACTER*80 OUTEXT

      STATO=0
C
C     Set the viewport limits of the graphics device
C
      VPMINX=VPXY(1)
      VPMINY=VPXY(3)
      VPMAXX=VPXY(2)
      VPMAXY=VPXY(4)
      CALL UGPVPT(GRDSCR,VPMINX,VPMAXX,VPMINY,VPMAXY,STAT2)
      IF (STAT2.NE.0) THEN
        WRITE(OUTEXT,21)
21      FORMAT(' Error setting viewport limits of graphics window')
        CALL UMSPUT(OUTEXT,3,0,STAT)
        STATO=1
        GO TO 990
      ENDIF
C
C     Set the World Coordinates of the data window
C
      WMINX=1.0
      WMAXX=REAL(N1*XYBIN)
      WMINY=1.0
      WMAXY=REAL(N2*XYBIN)
      CALL UGPWND(GRDSCR,WMINX,WMAXX,WMINY,WMAXY,STAT2)
      IF (STAT2.NE.0) THEN
        WRITE(OUTEXT,31)
31      FORMAT(' Error setting World Coordinates of graphics window')
        CALL UMSPUT(OUTEXT,3,0,STAT)
        STATO=1
        GO TO 990
      ENDIF
C
C     Get the character height and reset in proportion to the
C     viewport size
C
      CALL UGGGPR(GRDSCR,39,TYPSIZ,STAT)
      RESIZ=ABS(VPXY(2)-VPXY(1) + VPXY(4)-VPXY(3))
      TYPSIZ=RESIZ*TYPSIZ
      CALL UGPGPR(GRDSCR,39,TYPSIZ,STAT)
C
C     Plot the axis box and label
C
      IF (LLAB) THEN
        CALL UGLABL(GRDSCR,LABEL,XLAB,YLAB,STAT2)
        IF (STAT2.NE.0) THEN
          WRITE(OUTEXT,41)
41        FORMAT(' Error displaying and labelling plot frame')
          CALL UMSPUT(OUTEXT,3,0,STAT)
          STATO=1
          GO TO 990
        ENDIF
      ENDIF
C
C     Plot the polarization vectors on the unbinned grid at
C     the required sampling with a solid line
C
      NPT=2
      LINE=1 ! Solid line
      CALL UGPGPI(GRDSCR,GPLTYP,LINE,STAT2)
      DO I=XYSAM,N1-XYSAM,XYSAM
        DO J=XYSAM,N2-XYSAM,XYSAM
          IF (POL(I,J).GT.POLOW.AND.POL(I,J).LT.POHIGH) THEN
            ANG=-PA(I,J) + PAOFF
            POSX(1)=REAL(I*XYBIN) - POSCAL*0.5*POL(I,J)*SINRD(ANG)
            POSY(1)=REAL(J*XYBIN) - POSCAL*0.5*POL(I,J)*COSRD(ANG)
            POSX(2)=REAL(I*XYBIN) + POSCAL*0.5*POL(I,J)*SINRD(ANG)
            POSY(2)=REAL(J*XYBIN) + POSCAL*0.5*POL(I,J)*COSRD(ANG)
            CALL UGLINE(GRDSCR,POSX,POSY,NPT,STAT2)
          ENDIF
        ENDDO
      ENDDO
C
C     If error plotting required plot the vector with the
C     polarization - error at position angle - error and
C     the vector with the polarization + error at position
C     angle + error
C
      IF (LOER) THEN
        LINE=2 ! Dashed line
        CALL UGPGPI(GRDSCR,GPLTYP,LINE,STAT2)
        DO I=XYSAM,N1-XYSAM,XYSAM
          DO J=XYSAM,N2-XYSAM,XYSAM
            IF (POL(I,J).GT.POLOW.AND.POL(I,J).LT.POHIGH) THEN
              ANG=-PA(I,J) - PAER(I,J) + PAOFF
              POSX(1)=REAL(I*XYBIN) -
     :               POSCAL*0.5*ABS(POL(I,J)-POLER(I,J))*SINRD(ANG)
              POSY(1)=REAL(J*XYBIN) -
     :               POSCAL*0.5*ABS(POL(I,J)-POLER(I,J))*COSRD(ANG)
              POSX(2)=REAL(I*XYBIN) +
     :                POSCAL*0.5*ABS(POL(I,J)-POLER(I,J))*SINRD(ANG)
              POSY(2)=REAL(J*XYBIN) +
     :                POSCAL*0.5*ABS(POL(I,J)-POLER(I,J))*COSRD(ANG)
              CALL UGLINE(GRDSCR,POSX,POSY,NPT,STAT2)
            ENDIF
          ENDDO
        ENDDO
        DO I=XYSAM,N1-XYSAM,XYSAM
          DO J=XYSAM,N2-XYSAM,XYSAM
            IF (POL(I,J).GT.POLOW.AND.POL(I,J).LT.POHIGH) THEN
              ANG=-PA(I,J) + PAER(I,J) + PAOFF
              POSX(1)=REAL(I*XYBIN) -
     :               POSCAL*0.5*ABS(POL(I,J)+POLER(I,J))*SINRD(ANG)
              POSY(1)=REAL(J*XYBIN) -
     :               POSCAL*0.5*ABS(POL(I,J)+POLER(I,J))*COSRD(ANG)
              POSX(2)=REAL(I*XYBIN) +
     :                POSCAL*0.5*ABS(POL(I,J)+POLER(I,J))*SINRD(ANG)
              POSY(2)=REAL(J*XYBIN) +
     :                POSCAL*0.5*ABS(POL(I,J)+POLER(I,J))*COSRD(ANG)
              CALL UGLINE(GRDSCR,POSX,POSY,NPT,STAT2)
            ENDIF
          ENDDO
        ENDDO
      ENDIF
C
C     Plot a marker to indicate the polarization vector length
C     as 50% in the lower right corner
C
      LINE=1 ! Solid line
      PCLAB='50%'
      CALL UGPGPI(GRDSCR,GPLTYP,LINE,STAT2)
      XLFRAC=0.85
      YLFRAC=0.15
      POSX(1)=XLFRAC*REAL(N1*XYBIN) - 50.0*0.50*POSCAL
      POSX(2)=XLFRAC*REAL(N1*XYBIN) + 50.0*0.50*POSCAL
      POSY(1)=YLFRAC*REAL(N2*XYBIN)
      POSY(2)=YLFRAC*REAL(N2*XYBIN)
      CALL UGLINE(GRDSCR,POSX,POSY,NPT,STAT2)
      POSX(1)=XLFRAC*REAL(N1*XYBIN) - 50.0*0.50*POSCAL
      POSY(1)=1.04*YLFRAC*REAL(N2*XYBIN)
      CALL UGTEXT(GRDSCR,POSX(1),POSY(1),PCLAB,STAT2)

990   CONTINUE

999   END

max_stars_repo_name: iraf-community/stsdas


! ##################################################################################################################################
! 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.

      SUBROUTINE GET_SETID ( CARD, SETID )

! Gets SET ID from CASE CONTROL cards:  LOAD, METHOD, MPC, NLPARM, SPC, TEMP

      USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG, DOUBLE
      USE IOUNT1, ONLY                :  WRT_ERR, WRT_LOG, ERR, F04, F06
      USE SCONTR, ONLY                :  CC_ENTRY_LEN, FATAL_ERR, BLNK_SUB_NAM
      USE TIMDAT, ONLY                :  TSEC
      USE SUBR_BEGEND_LEVELS, ONLY    :  GET_SETID_BEGEND

      USE GET_SETID_USE_IFs

      IMPLICIT NONE

      CHARACTER(LEN=LEN(BLNK_SUB_NAM)):: SUBR_NAME = 'GET_SETID'
      CHARACTER(LEN=*), INTENT(IN)    :: CARD              ! A Case Control card (can be modified by subr CSHIFT, called herein)
      CHARACTER(LEN=LEN(CARD)+1)      :: CARD1             ! CARD shifted to begin in col after "=" sign
      CHARACTER(LEN=LEN(CARD)+1)      :: ERRTOK            ! An output from subr STOKEN, called herein
      CHARACTER( 3*BYTE)              :: EXCEPT            ! An input/output to/from subr STOKEN, called herein
      CHARACTER( 3*BYTE)              :: THRU              ! An inputoutput to/from subr STOKEN, called herein
      CHARACTER( 8*BYTE)              :: TOKEN(3)          ! An output from subr STOKEN, called herein
      CHARACTER( 8*BYTE)              :: TOKTYP(3)         ! An output from subr STOKEN, called herein

      INTEGER(LONG), INTENT(OUT)      :: SETID             ! Set ID read from CARD after '=', if CARD contains an integer here.
!                                                            SETID is set to -1 if 'ALL' is found after '=' & to 0 if 'NONE' found
      INTEGER(LONG)                   :: ECOL              ! Column on CARD where '=' is located
      INTEGER(LONG)                   :: IERR              ! An output from subr CSHIFT, called herein
      INTEGER(LONG)                   :: IERROR            ! An output from subr STOKEN, called herein
      INTEGER(LONG)                   :: IOCHK             ! IOSTAT error value from READ
      INTEGER(LONG)                   :: ISTART            ! An input to subr STOKEN, called herein
      INTEGER(LONG)                   :: NTOKEN            ! An output from subr STOKEN, called herein
      INTEGER(LONG)                   :: TOKLEN            ! An input to subr STOKEN, called herein
      INTEGER(LONG), PARAMETER        :: SUBR_BEGEND = GET_SETID_BEGEND

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9001) SUBR_NAME,TSEC
 9001    FORMAT(1X,A,' BEGN ',F10.3)
      ENDIF

! **********************************************************************************************************************************
! Get SETID

      SETID = 0
      CALL CSHIFT ( CARD, '=', CARD1, ECOL, IERR )
      IF (IERR /= 0) THEN
         FATAL_ERR = FATAL_ERR + 1
         WRITE(ERR,1703)
         WRITE(F06,1703)
      ELSE
         ISTART = 1
         THRU   = 'OFF'
         EXCEPT = 'OFF'
         TOKLEN = CC_ENTRY_LEN
         CALL STOKEN ( SUBR_NAME, CARD1, ISTART, TOKLEN, NTOKEN, IERROR, TOKTYP, TOKEN, ERRTOK, THRU, EXCEPT )
         IF (NTOKEN > 1) THEN
            FATAL_ERR = FATAL_ERR + 1
            WRITE(ERR,1704) NTOKEN
            WRITE(F06,1704) NTOKEN
         ELSE
            IF ((ISTART <= TOKLEN) .AND. (CARD1(ISTART:ISTART) /= '$')) THEN
               FATAL_ERR = FATAL_ERR + 1
               WRITE(ERR,1704) NTOKEN
               WRITE(F06,1704) NTOKEN
            ELSE
               IF (TOKTYP(1) == 'INTEGER ') THEN
                  READ(TOKEN(1),'(I8)',IOSTAT=IOCHK) SETID
                  IF ((IOCHK < 0) .OR. (IOCHK > 0)) THEN
                     FATAL_ERR = FATAL_ERR + 1
                     WRITE(ERR,1705)
                     WRITE(F06,1705)
                  ENDIF
                  IF (SETID <= 0) THEN
                     FATAL_ERR = FATAL_ERR + 1
                     WRITE(ERR,1706)
                     WRITE(F06,1706)
                  ENDIF
               ELSE
                  FATAL_ERR = FATAL_ERR + 1
                  WRITE(ERR,1707)
                  WRITE(F06,1707)
               ENDIF
            ENDIF
         ENDIF
      ENDIF

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9002) SUBR_NAME,TSEC
 9002    FORMAT(1X,A,' END  ',F10.3)
      ENDIF

      RETURN

! **********************************************************************************************************************************
 1703 FORMAT(' *ERROR  1703: CANNOT FIND EQUAL SIGN (=) ON ABOVE CASE CONTROL ENTRY ')

 1704 FORMAT(' *ERROR  1704: ABOVE CASE CONTROL ENTRY MUST HAVE ONLY ONE INTEGER SET NUMBER FOLLOWING THE "=" SIGN.'               &
                    ,/,14X,' ABOVE ENTRY DOES NOT MEET THIS REQUIREMENT (NTOKEN = ', I8, ')')

 1705 FORMAT(' *ERROR  1705: ERROR READING SET ID ON PREVIOUS CASE CONTROL ENTRY')

 1706 FORMAT(' *ERROR  1706: ZERO OR NEGATIVE SET ID NOT ALLOWED ON PREVIOUS CASE CONTROL ENTRY')

 1707 FORMAT(' *ERROR  1707: SET ID MUST BE AN INTEGER OF <= 8 DIGITS')

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

      END SUBROUTINE GET_SETID