!     path:      $Source$
!     author:    $Author$
!     revision:  $Revision$
!     created:   $Date$
!

      program rrtmg_lw

!----------------------------------------------------------------------------
! Copyright (c) 2002-2016, Atmospheric & Environmental Research, Inc. (AER)
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!  * Redistributions of source code must retain the above copyright
!    notice, this list of conditions and the following disclaimer.
!  * Redistributions in binary form must reproduce the above copyright
!    notice, this list of conditions and the following disclaimer in the
!    documentation and/or other materials provided with the distribution.
!  * Neither the name of Atmospheric & Environmental Research, Inc., nor
!    the names of its contributors may be used to endorse or promote products
!    derived from this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
! ARE DISCLAIMED. IN NO EVENT SHALL ATMOSPHERIC & ENVIRONMENTAL RESEARCH, INC.,
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
! THE POSSIBILITY OF SUCH DAMAGE.
!                        (http://www.rtweb.aer.com/)
!----------------------------------------------------------------------------
!
! ****************************************************************************
! *                                                                          *
! *                              RRTMG_LW                                    *
! *                                                                          *
! *                                                                          *
! *                                                                          *
! *                   a rapid radiative transfer model                       *
! *                       for the longwave region                            *
! *             for application to general circulation models                *
! *                                                                          *
! *                                                                          *
! *            Atmospheric and Environmental Research, Inc.                  *
! *                        131 Hartwell Avenue                               *
! *                        Lexington, MA 02421                               *
! *                                                                          *
! *                                                                          *
! *                           Eli J. Mlawer                                  *
! *                        Jennifer S. Delamere                              *
! *                         Michael J. Iacono                                *
! *                         Shepard A. Clough                                *
! *                                                                          *
! *                                                                          *
! *                                                                          *
! *                                                                          *
! *                                                                          *
! *                                                                          *
! *                       email:  miacono@aer.com                            *
! *                       email:  emlawer@aer.com                            *
! *                       email:  jdelamer@aer.com                           *
! *                                                                          *
! *        The authors wish to acknowledge the contributions of the          *
! *        following people:  Steven J. Taubman, Karen Cady-Pereira,         *
! *        Patrick D. Brown, Ronald E. Farren, Luke Chen, Robert Bergstrom.  *
! *                                                                          *
! ****************************************************************************

! -------- Description --------

! This program is the driver for RRTMG_LW, the AER LW radiation model for
! application to GCMs, that has been adapted from RRTM_LW for improved
! efficiency.
!
! This routine:
!    a) calls RRTMG_LW_INI to intialize data and to perform
!       the g-point interval reduction from 256 to 140.
!    b) calls READPROF to read in the atmospheric profile;
!       all layering in RRTMG is ordered from surface to toa.
!    c) calls CLDPROP to set cloud optical depth based on input
!       cloud properties, or CLDPRMC to set cloud optical depth
!       for McICA
!    d) calls SETCOEF to calculate various quantities needed for
!       the radiative transfer algorithm
!    e) calls TAUMOL to calculate gaseous optical depths for each
!       of the 16 spectral bands
!    f) calls RTRNMR (for both clear and cloudy profiles) to perform the
!       radiative transfer calculation with a maximum-random cloud
!       overlap method, or calls RTRN to use random cloud overlap,
!       or calls RTRNMC to use McICA, the Monte-Carlo Independent
!       Column Approximation, to represent sub-grid scale cloud
!       variability
!    g) writes out the calculated fluxes and cooling rates
!
! Two modes of operation are possible:
!     The mode is chosen by setting flag imca below.
!
!    1) Standard, single forward model calculation (imca = 0)
!    2) Monte Carlo Independent Column Approximation (McICA, Pincus et al.,
!       JC, 2003) method is applied to the forward model calculation (imca = 1)
!       For single column calculations, this method also requires setting flag
!       nmca below to the sample size of the Monte Carlo calculation;
!       (nmca = 200 is recommended)
!
! Two random number generators are available for use when imca = 1
!     This is chosen by setting flag irng below.
!
!    1) KISSVEC (irng = 0)
!    2) Mersenne Twister (irng = 1); the default setting
!
! Two methods of cloud property input are possible:
!     Cloud properties can be input in one of two ways (controlled by input
!     flags inflag, iceflag and liqflag; see text file rrtmg_lw_instructions
!     and subroutine rrtmg_lw_cldprop.f90 for further details):
!
!    1) Input cloud fraction and cloud optical depth directly (inflglw = 0)
!    2) Input cloud fraction and cloud physical properties (inflglw = 1 or 2);
!       cloud optical properties are calculated by cldprop or cldprmc based
!       on input settings of iceflglw and liqflglw.  Ice particle size provided
!       must be appropriately defined for the ice parameterization selected.
!
! One method of aerosol property input is possible:
!     Aerosol properties can be input in only one way (controlled by input
!     flag iaer; see text file rrtmg_lw_instructions further details):
!
!    1) Input aerosol optical depth directly by layer and spectral band (iaer=10);
!       band average optical depth at the mid-point of each spectral band.
!       RRTMG_LW currently treats only aerosol absorption;
!       scattering capability is not presently available.
!
! The optional calculation of the change in upward flux as a function of surface
! temperature is available (controlled by input flag idrv).  This can be utilized
! to approximate adjustments to the upward flux profile caused only by a change in
! surface temperature between full radiation calls.  This feature uses the pre-
! calculated derivative of the Planck function with respect to surface temperature.
!
!    1) Normal forward calculation for the input profile (idrv=0)
!    2) Normal forward calculation with optional calculation of the change
!       in upward flux as a function of surface temperature for clear sky
!       and total sky flux.  Flux partial derivatives are provided in arrays
!       dtotuflux_dt and dtotuclfl_dt for total and clear sky.  (idrv=1)
!
!
! ------- Modifications -------
!
! This version of RRTMG_LW has been modified from RRTM_LW to use a reduced
! set of g-points for application to GCMs.
!
!-- Original version (derived from RRTM_LW), reduction of g-points, other
!   revisions for use with GCMs.
!     1999: M. J. Iacono, AER, Inc.
!-- Adapted for use with NCAR/CAM.
!     May 2004: M. J. Iacono, AER, Inc.
!-- Revised to add McICA capability.
!     Nov 2005: M. J. Iacono, AER, Inc.
!-- Conversion to F90 formatting for consistency with rrtmg_sw.
!     Feb 2007: M. J. Iacono, AER, Inc.
!-- Added aerosol absorption capability.
!     Mar 2008: M. J. Iacono, AER, Inc.
!-- Added capability to adjust upward flux for a change only in surface temperature.
!     Nov 2009: M. J. Iacono, E. J. Mlawer, AER, Inc.

! --------- Modules ----------

      use parkind, only : im => kind_im, rb => kind_rb
      use parrrtm, only : mxlay, nbndlw, ngptlw, maxxsec, mxmol
      use rrlw_con, only: fluxfac, heatfac, oneminus, pi
      use rrlw_wvn, only: ng, ngb, nspa, nspb, wavenum1, wavenum2, delwave
      use rrlw_vsn
      use mcica_subcol_gen_lw, only: get_alpha, mcica_subcol_lw
      use rrtmg_lw_cldprop, only: cldprop
      use rrtmg_lw_cldprmc, only: cldprmc
      use rrtmg_lw_init, only: rrtmg_lw_ini
      use rrtmg_lw_rtrn, only: rtrn
      use rrtmg_lw_rtrnmr, only: rtrnmr
      use rrtmg_lw_rtrnmc, only: rtrnmc
      use rrtmg_lw_setcoef, only: setcoef
      use rrtmg_lw_taumol, only: taumol

      implicit none

! ------- Declarations -------

! ----- Local -----

! Control
      integer(kind=im) :: nlayers             ! total number of layers
      integer(kind=im) :: istart              ! beginning band of calculation
      integer(kind=im) :: iend                ! ending band of calculation
      integer(kind=im) :: icld                ! clear/cloud flag
      integer(kind=im) :: iflag               ! control flag
      integer(kind=im) :: iout                ! output option flag
      integer(kind=im) :: iaer                ! aerosol option flag
      integer(kind=im) :: ird                 ! input unit
      integer(kind=im) :: iwr                 ! output unit
      integer(kind=im) :: i                   ! output index
      integer(kind=im) :: ig                  ! g-point index
      integer(kind=im) :: iplon               ! column loop index
      integer(kind=im) :: ims                 ! mcica statistical loop index
      integer(kind=im) :: imca                ! flag for mcica [0=off, 1=on]
      integer(kind=im) :: nmca                ! number of mcica samples (mcica mode)
      integer(kind=im) :: irng                ! flag for random number generator
                                              ! [0=kissvec, 1=mersenne twister (default)]
      integer(kind=im) :: idrv                ! flag for calculation of dFdT, the change
                                              ! in upward flux as a function of surface
                                              ! temperature [0=off, 1=on]
      integer(kind=im) :: lev, l              ! level indices
      integer(kind=im), parameter :: ncol = 1 ! total number of columns
      character page

! Atmosphere
      real(kind=rb) :: pavel(mxlay)           ! layer pressures (mb)
      real(kind=rb) :: tavel(mxlay)           ! layer temperatures (K)
      real(kind=rb) :: dz(mxlay)              ! layer thickness (m)
      real(kind=rb) :: pz(0:mxlay)            ! level (interface) pressures (hPa, mb)
      real(kind=rb) :: tz(0:mxlay)            ! level (interface) temperatures (K)
      real(kind=rb) :: tbound                 ! surface temperature (K)
      real(kind=rb) :: dtbound                ! change in surface temperature for idrv=1 (K)
      real(kind=rb) :: coldry(mxlay)          ! dry air column density (mol/cm2)
      real(kind=rb) :: wbrodl(mxlay)          ! broadening gas column density (mol/cm2)
      real(kind=rb) :: wkl(mxmol,mxlay)       ! molecular amounts (mol/cm-2)
      real(kind=rb) :: wx(maxxsec,mxlay)      ! cross-section amounts (mol/cm-2)
      real(kind=rb) :: pwvcm                  ! precipitable water vapor (cm)
      real(kind=rb) :: semiss(nbndlw)         ! lw surface emissivity
      real(kind=rb) :: fracs(mxlay,ngptlw)    !
      real(kind=rb) :: taug(mxlay,ngptlw)     ! gaseous optical depths
      real(kind=rb) :: taut(mxlay,ngptlw)     ! gaseous + aerosol optical depths

      real(kind=rb) :: tauaer(mxlay,nbndlw)   ! aerosol optical depth
!      real(kind=rb) :: ssaaer(mxlay,nbndlw)  ! aerosol single scattering albedo
                                              ! for future expansion
                                              !   (lw aerosol scattering not yet available)
!      real(kind=rb) :: asmaer(mxlay,nbndlw)  ! aerosol asymmetry parameter
                                              ! for future expansion
                                              !   (lw aerosol scattering not yet available)

! Atmosphere - setcoef
      integer(kind=im) :: laytrop             ! tropopause layer index
      integer(kind=im) :: jp(mxlay)           !
      integer(kind=im) :: jt(mxlay)           !
      integer(kind=im) :: jt1(mxlay)          !
      real(kind=rb) :: planklay(mxlay,nbndlw)   !
      real(kind=rb) :: planklev(0:mxlay,nbndlw) !
      real(kind=rb) :: plankbnd(nbndlw)       !
      real(kind=rb) :: dplankbnd_dt(nbndlw)   !

      real(kind=rb) :: colh2o(mxlay)          ! column amount (h2o)
      real(kind=rb) :: colco2(mxlay)          ! column amount (co2)
      real(kind=rb) :: colo3(mxlay)           ! column amount (o3)
      real(kind=rb) :: coln2o(mxlay)          ! column amount (n2o)
      real(kind=rb) :: colco(mxlay)           ! column amount (co)
      real(kind=rb) :: colch4(mxlay)          ! column amount (ch4)
      real(kind=rb) :: colo2(mxlay)           ! column amount (o2)
      real(kind=rb) :: colbrd(mxlay)          ! column amount (broadening gases)

      integer(kind=im) :: indself(mxlay)
      integer(kind=im) :: indfor(mxlay)
      real(kind=rb) :: selffac(mxlay)
      real(kind=rb) :: selffrac(mxlay)
      real(kind=rb) :: forfac(mxlay)
      real(kind=rb) :: forfrac(mxlay)

      integer(kind=im) :: indminor(mxlay)
      real(kind=rb) :: minorfrac(mxlay)
      real(kind=rb) :: scaleminor(mxlay)
      real(kind=rb) :: scaleminorn2(mxlay)

      real(kind=rb) :: &                      !
                         fac00(mxlay), fac01(mxlay), &
                         fac10(mxlay), fac11(mxlay)
      real(kind=rb) :: &                      !
                         rat_h2oco2(mxlay),rat_h2oco2_1(mxlay), &
                         rat_h2oo3(mxlay),rat_h2oo3_1(mxlay), &
                         rat_h2on2o(mxlay),rat_h2on2o_1(mxlay), &
                         rat_h2och4(mxlay),rat_h2och4_1(mxlay), &
                         rat_n2oco2(mxlay),rat_n2oco2_1(mxlay), &
                         rat_o3co2(mxlay),rat_o3co2_1(mxlay)

! Atmosphere/clouds - cldprop
      integer(kind=im) :: ncbands             ! number of cloud spectral bands
      integer(kind=im) :: inflag              ! flag for cloud property method
      integer(kind=im) :: iceflag             ! flag for ice cloud properties
      integer(kind=im) :: liqflag             ! flag for liquid cloud properties

      real(kind=rb) :: cldfrac(mxlay)         ! layer cloud fraction
      real(kind=rb) :: tauc(nbndlw,mxlay)     ! in-cloud optical depth (non-delta scaled)
!      real(kind=rb) :: ssac(nbndlw,mxlay)    ! in-cloud single scattering albedo (non-delta scaled)
                                              ! for future expansion
                                              !   (lw scattering not yet available)
!      real(kind=rb) :: asmc(nbndlw,mxlay)    ! in-cloud asymmetry parameter (non-delta scaled)
                                              ! for future expansion
                                              !   (lw scattering not yet available)
      real(kind=rb) :: ciwp(mxlay)            ! in-cloud ice water path
      real(kind=rb) :: clwp(mxlay)            ! in-cloud liquid water path
      real(kind=rb) :: rei(mxlay)             ! cloud ice particle effective size (microns)
                                              ! specific definition of rei depends on setting of iceflag:
                                              ! iceflag = 0: ice effective radius, r_ec, (Ebert and Curry, 1992),
                                              !              r_ec must be >= 10.0 microns
                                              ! iceflag = 1: ice effective radius, r_ec, (Ebert and Curry, 1992),
                                              !              r_ec range is limited to 13.0 to 130.0 microns
                                              ! iceflag = 2: ice effective radius, r_k, (Key, Streamer Ref. Manual, 1996)
                                              !              r_k range is limited to 5.0 to 131.0 microns
                                              ! iceflag = 3: generalized effective size, dge, (Fu, 1996),
                                              !              dge range is limited to 5.0 to 140.0 microns
                                              !              [dge = 1.0315 * r_ec]
      real(kind=rb) :: rel(mxlay)             ! cloud liquid particle effective radius (microns)

      real(kind=rb) :: taucloud(mxlay,nbndlw) ! in-cloud optical depth; delta scaled
!      real(kind=rb) :: ssacloud(mxlay,nbndlw)! in-cloud single scattering albedo; delta scaled
                                              ! for future expansion
                                              !   (lw scattering not yet available)
!      real(kind=rb) :: asmcloud(mxlay,nbndlw)! in-cloud asymmetry parameter; delta scaled
                                              ! for future expansion
                                              !   (lw scattering not yet available)

! Atmosphere/clouds - cldprmc [mcica]
      real(kind=rb) :: cldfmc(ngptlw,mxlay)   ! cloud fraction [mcica]
      real(kind=rb) :: ciwpmc(ngptlw,mxlay)   ! in-cloud ice water path [mcica]
      real(kind=rb) :: clwpmc(ngptlw,mxlay)   ! in-cloud liquid water path [mcica]
      real(kind=rb) :: relqmc(mxlay)          ! liquid particle effective radius (microns)
      real(kind=rb) :: reicmc(mxlay)          ! ice particle effective radius (microns)
      real(kind=rb) :: taucmc(ngptlw,mxlay)   ! in-cloud optical depth [mcica]
!      real(kind=rb) :: ssacmc(ngptlw,mxlay)  ! in-cloud single scattering albedo [mcica]
                                              ! for future expansion
                                              !   (lw scattering not yet available)
!      real(kind=rb) :: asmcmc(ngptlw,mxlay)  ! in-cloud asymmetry parameter [mcica]
                                              ! for future expansion

                                              !   (lw scattering not yet available)

      real(kind=rb) :: dtotuflux_dt(0:mxlay)  ! change in upward longwave flux (w/m2/k)
                                              ! with respect to surface temperature
      real(kind=rb) :: dtotuclfl_dt(0:mxlay)  ! change in clear sky upward longwave flux (w/m2/k)
                                              ! with respect to surface temperature
      real(kind=rb) :: alpha(mxlay)           ! vertical cloud fraction correlation parameter
      real(kind=rb) :: decorr_con             ! cloud fraction decorrelation length, constant (meters)
      integer(kind=im) :: juldat              ! Julian date (day of year)
      integer(kind=im) :: idcor               ! Decorrelation length type
      real(kind=rb) :: lat                    ! latitude (degrees)

! Parameters
      real(kind=rb), parameter :: cpdair = 1.004e3_rb  ! Specific heat capacity of dry air
                                                       ! at constant pressure at 273 K
                                                       ! (J kg-1 K-1)
! Output
      real(kind=rb) :: totuflux(0:mxlay)      ! upward longwave flux (w/m2)
      real(kind=rb) :: totdflux(0:mxlay)      ! downward longwave flux (w/m2)
      real(kind=rb) :: fnet(0:mxlay)          ! net longwave flux (w/m2)
      real(kind=rb) :: htr(0:mxlay)           ! longwave heating rate (k/day)
      real(kind=rb) :: totuclfl(0:mxlay)      ! clear sky upward longwave flux (w/m2)
      real(kind=rb) :: totdclfl(0:mxlay)      ! clear sky downward longwave flux (w/m2)
      real(kind=rb) :: fnetc(0:mxlay)         ! clear sky net longwave flux (w/m2)
      real(kind=rb) :: htrc(0:mxlay)          ! clear sky longwave heating rate (k/day)
! Output (mean output for McICA calculation)
      real(kind=rb) :: uflxsum(0:mxlay)       ! upward longwave flux (w/m2)
      real(kind=rb) :: dflxsum(0:mxlay)       ! downward longwave flux (w/m2)
      real(kind=rb) :: fnetsum(0:mxlay)       ! net longwave flux (w/m2)
      real(kind=rb) :: htrsum(0:mxlay)        ! longwave heating rate (k/day)

!
! Initializations

      hvrrtm = '$Revision$'
      hvrini = 'NOT USED'
      hvrcld = 'NOT USED'
      hvrclc = 'NOT USED'
      hvrrtr = 'NOT USED'
      hvrrtx = 'NOT USED'
      hvrrtc = 'NOT USED'
      hvrset = 'NOT USED'
      hvrtau = 'NOT USED'
      hvrkg  = '$Revision$'
      hvratm = 'NOT USED'
      hvrutl = 'NOT USED'
      hvrext = 'NOT USED'

      hnamrtm = '  rrtmg_lw.1col.f90:'
      hnamini = '  rrtmg_lw_init.f90:'
      hnamcld = 'rrtmglw_cldprop.f90:'
      hnamclc = 'rrtmglw_cldprmc.f90:'
      hnamrtr = '  rrtmg_lw_rtrn.f90:'
      hnamrtx = 'rrtmg_lw_rtrnmr.f90:'
      hnamrtc = 'rrtmg_lw_rtrnmc.f90:'
      hnamset = 'rrtmglw_setcoef.f90:'
      hnamtau = 'rrtmg_lw_taumol.f90:'
      hnamkg  = '   rrtmg_lw_k_g.f90:'
      hnamatm = '           rrtatm.f:'
      hnamutl = '         util_xxx.f:'
      hnamext = '            extra.f:'

      oneminus = 1._rb - 1.e-6_rb
      pi = 2._rb * asin(1._rb)
      fluxfac = pi * 2.e4_rb                ! orig:   fluxfac = pi * 2.d4
      ird = 9
      iwr = 10
      page = char(12)

      uflxsum(0:) = 0._rb
      dflxsum(0:) = 0._rb
      fnetsum(0:) = 0._rb
      htrsum(0:) = 0._rb

! Set imca to select calculation type
!  (read by subroutine readprof from input file INPUT_RRTM):
! imca = 0, use standard forward model calculation
! imca = 1, use McICA for Monte Carlo treatment of sub-grid cloud variability

! Set irng to select random number generator for McICA (used when imca = 1)
! irng = 0, KISSVEC
! irng = 1, Mersenne Twister
!      irng = 0
      irng = 1

! Set icld to select of clear or cloud calculation and cloud overlap method
!  (read by subroutine readprof from input file INPUT_RRTM):
! icld = 0, clear only
! icld = 1, with clouds using random cloud overlap
! icld = 2, with clouds using maximum_random cloud overlap
! icld = 3, with clouds using maximum cloud overlap (McICA only)
! icld = 4, with clouds using exponential cloud overlap (McICA only)
! icld = 5, with clouds using exponential_random cloud overlap (McICA only)

! Call model and data initialization, compute lookup tables, perform
! reduction of g-points from 256 to 140 for input absorption
! coefficient data and other arrays.
!
! In a GCM this call should be placed in the model initialization
! area, since this has to be called only once.

      call rrtmg_lw_ini(cpdair)

! Open the input set of atmospheres
      open (ird,file='INPUT_RRTM',form='formatted')
! Open the output file
      open (iwr,file='OUTPUT_RRTM',form='formatted')

! This is the main longitude/column loop within rrtmg.

      do iplon = 1, ncol

! Input atmospheric profile from INPUT_RRTM.
         call readprof(ird, nlayers, iout, imca, icld, iaer, idrv, &
                       pavel, tavel, pz, tz, tbound, dz, semiss, &
                       dtbound, coldry, wkl, wbrodl, wx, pwvcm, &
                       inflag, iceflag, liqflag, cldfrac, &
                       tauc, ciwp, clwp, rei, rel, tauaer, &
                       decorr_con, juldat, idcor, lat)

         istart = 1
         iend = 16
         iflag = iout

! Set nmca to sample size for Monte Carlo calculation
         if (imca.eq.0) nmca = 1
         if (imca.eq.1) nmca = 200

! Return here for multiple band output
 1000    continue
         if (iflag .gt. 0 .and. iflag .le. 40) then
            istart = iflag
            iend = iflag
         endif

! This is the statistical sampling loop for McICA

         do ims = 1, nmca

! Call sub-colum cloud generator for McICA calculations
! Output will be written for all nmca samples.  This will be excessive if
! band output (iout=99) option is selected.

            if (imca.eq.1) then
               alpha(:) = 0.0_rb
               if (icld .eq. 4 .or. icld .eq. 5) then
                  call get_alpha(iplon, nlayers, icld, idcor, decorr_con, &
                                 dz, lat, juldat, cldfrac, alpha)
               endif
               call mcica_subcol_lw(iplon, nlayers, icld, ims, irng, pavel, &
                          cldfrac, ciwp, clwp, rei, rel, tauc, alpha, &
                          cldfmc, ciwpmc, clwpmc, reicmc, relqmc, taucmc)
            endif

!  For cloudy atmosphere, use cldprop to set cloud optical properties based on
!  input cloud physical properties.  Select method based on choices described
!  in cldprop.  Cloud fraction, water path, liquid droplet and ice particle
!  effective radius must be passed into cldprop.  Cloud fraction and cloud
!  optical depth are transferred to rrtmg_lw arrays in cldprop.

!  If McICA is requested use cloud fraction and cloud physical properties
!  generated by sub-column cloud generator above.

            if (imca.eq.0) then
               call cldprop(nlayers, inflag, iceflag, liqflag, cldfrac, tauc, &
                            ciwp, clwp, rei, rel, ncbands, taucloud)
            else
               call cldprmc(nlayers, inflag, iceflag, liqflag, cldfmc, &
                            ciwpmc, clwpmc, reicmc, relqmc, ncbands, taucmc)
            endif

! Calculate information needed by the radiative transfer routine
! that is specific to this atmosphere, especially some of the
! coefficients and indices needed to compute the optical depths
! by interpolating data from stored reference atmospheres.

            call setcoef(nlayers, istart, pavel, tavel, tz, tbound, semiss, &
                         coldry, wkl, wbrodl, &
                         laytrop, jp, jt, jt1, planklay, planklev, plankbnd, &
                         idrv, dplankbnd_dt, &
                         colh2o, colco2, colo3, coln2o, colco, colch4, colo2, &
                         colbrd, fac00, fac01, fac10, fac11, &
                         rat_h2oco2, rat_h2oco2_1, rat_h2oo3, rat_h2oo3_1, &
                         rat_h2on2o, rat_h2on2o_1, rat_h2och4, rat_h2och4_1, &
                         rat_n2oco2, rat_n2oco2_1, rat_o3co2, rat_o3co2_1, &
                         selffac, selffrac, indself, forfac, forfrac, indfor, &
                         minorfrac, scaleminor, scaleminorn2, indminor)

!  Calculate the gaseous optical depths and Planck fractions for
!  each longwave spectral band.

            call taumol(nlayers, pavel, wx, coldry, &
                        laytrop, jp, jt, jt1, planklay, planklev, plankbnd, &
                        colh2o, colco2, colo3, coln2o, colco, colch4, colo2, &
                        colbrd, fac00, fac01, fac10, fac11, &
                        rat_h2oco2, rat_h2oco2_1, rat_h2oo3, rat_h2oo3_1, &
                        rat_h2on2o, rat_h2on2o_1, rat_h2och4, rat_h2och4_1, &
                        rat_n2oco2, rat_n2oco2_1, rat_o3co2, rat_o3co2_1, &
                        selffac, selffrac, indself, forfac, forfrac, indfor, &
                        minorfrac, scaleminor, scaleminorn2, indminor, &
                        fracs, taug)

! Combine gaseous and aerosol optical depths, if aerosol active
            if (iaer .eq. 0) then
               do i = 1, nlayers
                  do ig = 1, ngptlw
                     taut(i,ig) = taug(i,ig)
                  enddo
               enddo
            elseif (iaer .eq. 10) then
               do i = 1, nlayers
                  do ig = 1, ngptlw
                     taut(i,ig) = taug(i,ig) + tauaer(i,ngb(ig))
                  enddo
               enddo
            endif

! Call the radiative transfer routine.
! Either routine can be called to do clear sky calculation.  If clouds
! are present, then select routine based on cloud overlap assumption
! to be used.  Clear sky calculation is done simultaneously.
! For McICA, only RTRN is called for clear and cloudy calculations.

            if (imca .eq. 0) then
               if (icld .eq. 1) then
                  call rtrn(nlayers, istart, iend, iout, pz, semiss, ncbands, &
                            cldfrac, taucloud, planklay, planklev, plankbnd, &
                            pwvcm, fracs, taut, &
                            totuflux, totdflux, fnet, htr, &
                            totuclfl, totdclfl, fnetc, htrc, &
                            idrv, dplankbnd_dt, dtotuflux_dt, dtotuclfl_dt )
               else
                  call rtrnmr(nlayers, istart, iend, iout, pz, semiss, ncbands, &
                              cldfrac, taucloud, planklay, planklev, plankbnd, &
                              pwvcm, fracs, taut, &
                              totuflux, totdflux, fnet, htr, &
                              totuclfl, totdclfl, fnetc, htrc, &
                              idrv, dplankbnd_dt, dtotuflux_dt, dtotuclfl_dt )
               endif
            elseif (imca .eq. 1) then
               call rtrnmc(nlayers, istart, iend, iout, pz, semiss, ncbands, &
                           cldfmc, taucmc, planklay, planklev, plankbnd, &
                           pwvcm, fracs, taut, &
                           totuflux, totdflux, fnet, htr, &
                           totuclfl, totdclfl, fnetc, htrc, &
                           idrv, dplankbnd_dt, dtotuflux_dt, dtotuclfl_dt )
            endif

! If requested (idrv=1), adjust upward fluxes and heating rate for the input
! change in surface temperature (dtbound) based on the calculated change in
! upward flux as a function of surface temperature for each layer.  Downward
! fluxes are unchanged by this adjustment.

            if (idrv .eq. 1) then
! Adjust fluxes at surface
               totuflux(0) = totuflux(0) + dtotuflux_dt(0) * dtbound
               fnet(0) = totuflux(0) - totdflux(0)
               totuclfl(0) = totuclfl(0) + dtotuclfl_dt(0) * dtbound
               fnetc(0) = totuclfl(0) - totdclfl(0)

! Adjust fluxes at model levels
               do lev = 1, nlayers
                  totuflux(lev) = totuflux(lev) + dtotuflux_dt(lev) * dtbound
                  fnet(lev) = totuflux(lev) - totdflux(lev)
                  totuclfl(lev) = totuclfl(lev) + dtotuclfl_dt(lev) * dtbound
                  fnetc(lev) = totuclfl(lev) - totdclfl(lev)
                  l = lev - 1

! Re-calculate heating rates at model layers
                  htr(l)=heatfac*(fnet(l)-fnet(lev))/(pz(l)-pz(lev))
                  htrc(l)=heatfac*(fnetc(l)-fnetc(lev))/(pz(l)-pz(lev))
               enddo

! Set heating rate to zero in top layer
               htr(nlayers) = 0.0_rb
               htrc(nlayers) = 0.0_rb
            endif

! Process output.
            if (iout .lt. 0) goto 2000

            if (imca .eq. 0) then

               write(iwr,9899)wavenum1(istart),wavenum2(iend),iplon
               write(iwr,9900)
               write(iwr,9901)

               do i = nlayers, 0, -1
                  if (pz(i) .lt. 1.e-2_rb) then
                     write(iwr,9952) i,pz(i),totuflux(i),totdflux(i),fnet(i),htr(i)
                  elseif (pz(i) .lt. 1.e-1) then
                     write(iwr,9953) i,pz(i),totuflux(i),totdflux(i),fnet(i),htr(i)
                  elseif (pz(i) .lt. 1.) then
                     write(iwr,9954) i,pz(i),totuflux(i),totdflux(i),fnet(i),htr(i)
                  elseif (pz(i) .lt. 10.) then
                     write(iwr,9955) i,pz(i),totuflux(i),totdflux(i),fnet(i),htr(i)
                  elseif (pz(i) .lt. 100.) then
                     write(iwr,9956) i,pz(i),totuflux(i),totdflux(i),fnet(i),htr(i)
                  elseif (pz(i) .lt. 1000.) then
                     write(iwr,9957) i,pz(i),totuflux(i),totdflux(i),fnet(i),htr(i)
                  else
                     write(iwr,9958) i,pz(i),totuflux(i),totdflux(i),fnet(i),htr(i)
                  endif
               enddo
               write(iwr,9903)page

            elseif (imca .eq. 1) then

               do i = nlayers, 0, -1
                  uflxsum(i) = uflxsum(i) + totuflux(i)
                  dflxsum(i) = dflxsum(i) + totdflux(i)
                  fnetsum(i) = fnetsum(i) + fnet(i)
                  htrsum(i) = htrsum(i) + htr(i)
               enddo

! Output average over samples when last sample reached.  Comment this if-check
! and swap flux output write statements below to output all McICA samples.
               if (ims .eq. nmca) then

                  do i = nlayers, 0, -1
                     uflxsum(i) = uflxsum(i)/nmca
                     dflxsum(i) = dflxsum(i)/nmca
                     fnetsum(i) = fnetsum(i)/nmca
                     htrsum(i) = htrsum(i)/nmca
                  enddo

                  write(iwr,9899) wavenum1(istart),wavenum2(iend),iplon
                  write(iwr,9900)
                  write(iwr,9901)

                  do i = nlayers, 0, -1
                     if (pz(i) .lt. 1.e-2_rb) then
                        write(iwr,9952) i,pz(i),uflxsum(i),dflxsum(i),fnetsum(i),htrsum(i)
                     elseif (pz(i) .lt. 1.e-1) then
                        write(iwr,9953) i,pz(i),uflxsum(i),dflxsum(i),fnetsum(i),htrsum(i)
                     elseif (pz(i) .lt. 1.) then
                        write(iwr,9954) i,pz(i),uflxsum(i),dflxsum(i),fnetsum(i),htrsum(i)
                     elseif (pz(i) .lt. 10.) then
                        write(iwr,9955) i,pz(i),uflxsum(i),dflxsum(i),fnetsum(i),htrsum(i)
                     elseif (pz(i) .lt. 100.) then
                        write(iwr,9956) i,pz(i),uflxsum(i),dflxsum(i),fnetsum(i),htrsum(i)
                     elseif (pz(i) .lt. 1000.) then
                        write(iwr,9957) i,pz(i),uflxsum(i),dflxsum(i),fnetsum(i),htrsum(i)
                     else
                        write(iwr,9958) i,pz(i),uflxsum(i),dflxsum(i),fnetsum(i),htrsum(i)
                     endif
                  enddo
                  write(iwr,9903)page

!               do i = nlayers, 0, -1
!                  if (pz(i) .lt. 1.e-2_rb) then
!                     write(iwr,9952) i,pz(i),totuflux(i),totdflux(i),fnet(i),htr(i)
!                  elseif (pz(i) .lt. 1.e-1) then
!                     write(iwr,9953) i,pz(i),totuflux(i),totdflux(i),fnet(i),htr(i)
!                  elseif (pz(i) .lt. 1.) then
!                     write(iwr,9954) i,pz(i),totuflux(i),totdflux(i),fnet(i),htr(i)
!                  elseif (pz(i) .lt. 10.) then
!                     write(iwr,9955) i,pz(i),totuflux(i),totdflux(i),fnet(i),htr(i)
!                  elseif (pz(i) .lt. 100.) then
!                     write(iwr,9956) i,pz(i),totuflux(i),totdflux(i),fnet(i),htr(i)
!                  elseif (pz(i) .lt. 1000.) then
!                     write(iwr,9957) i,pz(i),totuflux(i),totdflux(i),fnet(i),htr(i)
!                  else
!                     write(iwr,9958) i,pz(i),totuflux(i),totdflux(i),fnet(i),htr(i)
!                  endif
!               enddo
!               write(iwr,9903)page

               endif

            endif

 2000       continue
! End statistical loop for McICA
         enddo

         if (iout .le. 40 .or. iflag .eq. 16) goto 3500
         if (iflag .eq. 99) then
            iflag = 1
         elseif (iout .eq. 99) then
            iflag = iflag + 1
         endif
         goto 1000

 3500    continue

!
! Output module version numbers
!
         write(iwr,9910) hnamrtm,hvrrtm,hnamini,hvrini,hnamcld,hvrcld, &
           hnamclc,hvrclc,hnamrtr,hvrrtr,hnamrtx,hvrrtx,hnamrtc,hvrrtc, &
           hnamset,hvrset,hnamtau,hvrtau,hnamatm,hvratm,hnamutl,hvrutl, &
           hnamext,hvrext,hnamkg,hvrkg

 4000    continue
! 4444    continue

! End longitude/column loop
      enddo

      close(ird)
      close(iwr)

 9952 format(1x,i3,9x,f7.6,3x,f8.4,6x,f8.4,6x,f12.7,10x,f9.5)
 9953 format(1x,i3,9x,f6.5,4x,f8.4,6x,f8.4,6x,f12.7,10x,f9.5)
 9954 format(1x,i3,8x,f6.4,5x,f8.4,6x,f8.4,6x,f12.7,10x,f9.5)
 9955 format(1x,i3,7x,f6.3,6x,f8.4,6x,f8.4,6x,f12.7,10x,f9.5)
 9956 format(1x,i3,6x,f6.2,7x,f8.4,6x,f8.4,6x,f12.7,10x,f9.5)
 9957 format(1x,i3,5x,f6.1,8x,f8.4,6x,f8.4,6x,f12.7,10x,f9.5)
 9958 format(1x,i3,5x,f6.1,8x,f8.4,6x,f8.4,6x,f12.7,10x,f9.5)
 9899 format(1x,'Wavenumbers: ',f6.1,' - ',f6.1,' cm-1, ATM ',i6)
 9900 format(1x,'LEVEL    PRESSURE   UPWARD FLUX   DOWNWARD FLUX    NET FLUX       HEATING RATE')
 9901 format(1x,'            mb          W/m2          W/m2           W/m2          degree/day')
 9902 format(1x,i3,3x,f11.6,4x,1p,2(g12.6,2x),g13.6,3x,g16.9,0p)
 9903 format(a)
 9910 format('  Modules and versions used in this calculation:',/,/, &
              7(5x,a20,2x,a18,10x,a20,2x,a18,/))

      contains

!************************************************************************
      subroutine readprof(ird_in, nlayers_out, iout_out, imca, icld_out, &
          iaer_out, idrv, pavel_out, tavel_out, pz_out, tz_out, tbound_out, dz_out, &
          semiss_out, dtbound_out, coldry_out, wkl_out, wbrodl_out, wx_out, &
          pwvcm_out, inflag_out, iceflag_out, liqflag_out, cldfrac_out, &
          tauc, ciwp, clwp, rei, rel, tauaer_out, &
          decorr_con_out, juldat_out, idcor_out, lat_out)
!************************************************************************

! --------- Modules ----------

      use parkind, only : im => kind_im, rb => kind_rb
      use rrlw_con, only: pi, grav, planck, boltz, clight, avogad, alosmt, &
                          gascon, radcn1, radcn2

! Note: COMMON blocks are left in this routine for array passing with
!       rrtatm.f for reading input profiles in single column mode.
!       Scalars and arrays in subroutine call are renamed to avoid conflict
!       with COMMON blocks. Variable mxlay is the maximum possible number of
!       layers, which is used to dimension the input arrays, while nlayers
!       is the actual number of input layers.
!
! Purpose: Read in atmospheric profile.

      implicit integer(i-n), real (a-h,o-z)

! ------- Parameters -------
      parameter (mxlay = 603)
      parameter (mxmol = 39)
      parameter (nbndlw = 16)
      parameter (maxinpx = mxmol)
      parameter (maxxsec = 4)
      parameter (maxprod = mxlay*maxxsec)
      parameter (amd = 28.9660_rb)         ! Effective molecular weight of dry air (g/mol)
      parameter (amw = 18.0160_rb)         ! Molecular weight of water vapor (g/mol)

      dimension altz(0:mxlay),ixtrans(14),semis(nbndlw)

      common /consts/   pic,planckc,boltzc,clightc,avogadc,alosmtc,gasconc, &
                        radcn1c,radcn2c
      common /control/  numangs, iout, istart, iend, icld, iaer
      common /profile/  nlayers,pavel(mxlay),tavel(mxlay),pz(0:mxlay),tz(0:mxlay)
      common /surface/  tbound,ireflect,semiss(nbndlw)
      common /species/  coldry(mxlay),wkl(mxmol,mxlay),wbrodl(mxlay),colmol(mxlay),nmol
      common /ifil/     ird,ipr,ipu,idum(15)
      common /xsecctrl/ nxmol,ixindx(maxinpx)
      common /xsec/     wx(maxxsec,mxlay)
      common /pathx/    ixmax,nxmol0,ixindx0(maxinpx),wx0(maxinpx,mxlay)
      common /xrrtatm/  ixsect
      common /c_drive/  ref_lat

      common /cloudin/   inflag,clddat1(mxlay),clddat2(mxlay), &
                         iceflag,liqflag,clddat3(mxlay),clddat4(mxlay)
      common /clouddat/  ncbands,cldfrac(mxlay),taucloud(mxlay,nbndlw)
      common /aerdat/    tauaer(mxlay,nbndlw)

      character*80 form1(0:1),form2(0:1),form3(0:1)
      character*1 ctest, cdollar, cprcnt,cdum
      character*3 cinp

! Dimensions for transfer to rrtmg
      integer(kind=im), intent(in) :: ird_in              ! input file unit
      integer(kind=im), intent(out) :: nlayers_out        ! total number of layers
      integer(kind=im), intent(out) :: icld_out           ! clear/cloud flag
      integer(kind=im), intent(out) :: imca               ! McICA on/off flag (1 = use McICA)
      integer(kind=im), intent(out) :: iout_out           ! output option flag
      integer(kind=im), intent(out) :: iaer_out           ! aerosol option flag
      integer(kind=im), intent(out) :: idrv               ! Planck derivative option on/off flag
                                                          ! (1 = provide upward flux adjustment
                                                          ! for change in surface temperature

      real(kind=rb), intent(out) :: pavel_out(mxlay)      ! layer pressures (mb)
      real(kind=rb), intent(out) :: tavel_out(mxlay)      ! layer temperatures (K)
      real(kind=rb), intent(out) :: dz_out(mxlay)         ! layer thickness (m)
      real(kind=rb), intent(out) :: pz_out(0:mxlay)       ! level (interface) pressures (hPa, mb)
      real(kind=rb), intent(out) :: tz_out(0:mxlay)       ! level (interface) temperatures (K)
      real(kind=rb), intent(out) :: tbound_out            ! surface temperature (K)
      real(kind=rb), intent(out) :: dtbound_out           ! surface temperature change for idrv=1 (K)
      real(kind=rb), intent(out) :: coldry_out(mxlay)     ! dry air column density (mol/cm2)
      real(kind=rb), intent(out) :: wbrodl_out(mxlay)     ! broadening gas column density (mol/cm2)
      real(kind=rb), intent(out) :: wkl_out(mxmol,mxlay)  ! molecular amounts (mol/cm2)
      real(kind=rb), intent(out) :: wx_out(maxxsec,mxlay) ! cross-section amounts (mol/cm2)
      real(kind=rb), intent(out) :: pwvcm_out             ! precipitable water vapor (cm)
      real(kind=rb), intent(out) :: semiss_out(nbndlw)    ! lw surface emissivity
      real(kind=rb), intent(out) :: decorr_con_out        ! cloud fraction decorrelation length, constant (meters)
      real(kind=rb), intent(out) :: lat_out               ! latitude
      integer(kind=im), intent(out) :: juldat_out         ! Julian date (day of year)
      integer(kind=im), intent(out) :: idcor_out          ! Decorrelation length type

      integer(kind=im), intent(out) :: inflag_out         ! cloud property option flag
      integer(kind=im), intent(out) :: iceflag_out        ! ice cloud property flag
      integer(kind=im), intent(out) :: liqflag_out        ! liquid cloud property flag

      real(kind=rb), intent(out) :: cldfrac_out(mxlay)    ! cloud fraction
      real(kind=rb), intent(out) :: tauc(nbndlw,mxlay)    ! in-cloud optical depth
!      real(kind=rb), intent(out) :: ssac(nbndlw,mxlay)   ! in-cloud single scattering albedo
                                                          !   for future expansion
!      real(kind=rb), intent(out) :: asmc(nbndlw,mxlay)   ! in-cloud asymmetry parameter
                                                          !   for future expansion
      real(kind=rb), intent(out) :: ciwp(mxlay)           ! in-cloud ice water path
      real(kind=rb), intent(out) :: clwp(mxlay)           ! in-cloud liquid water path
      real(kind=rb), intent(out) :: rel(mxlay)            ! cloud liquid particle effective radius (microns)
      real(kind=rb), intent(out) :: rei(mxlay)            ! cloud ice particle effective size (microns)
      real(kind=rb), intent(out) :: tauaer_out(mxlay,nbndlw)  ! aerosol optical depth
!      real(kind=rb), intent(out) :: ssaaer_out(mxlay,nbndlw)  ! aerosol single scattering albedo
                                                               !   for future expansion
!      real(kind=rb), intent(out) :: asmaer_out(mxlay,nbndlw)  ! aerosol asymmetry parameter
                                                               !   for future expansion

! Local
      real(kind=rb) :: fice(mxlay)                        ! cloud ice fraction

      real(kind=rb) :: amttl                              ! moist air vertical sum (molecular amount)
      real(kind=rb) :: wvttl                              ! water vapor vertical sum (molecular amount)
      real(kind=rb) :: summol                             ! sum over non-water molecules
      real(kind=rb) :: wvsh                               ! water vapor vertical total specific humitidy
      real(kind=rb) :: pwvcm                              ! precipitable water vapor (cm)
      real(kind=rb) :: dtbound                            ! change in surface temperature for idrv=1 (K)

! Local inputs for EXP or ER cloud overlap
      integer(kind=im) :: idcor                           ! decorrelation length type
                                                          !  0 = constant
                                                          !  1 = latitude and day of year dependent
      real(kind=rb) :: decorr_con                         ! cloud fraction decorrelation length, constant (meters)
      integer(kind=im) :: juldat                          ! Julian date (day of year)
      real(kind=rb) :: lat                                ! latitude (degrees)

!

! Initializations

      data cdollar /'$'/
      data cprcnt /'%'/
      data ixtrans /0,0,0,1,2,3,0,0,0,0,0,4,0,0/

      form1(0) = '(3f10.4,a3,i2,1x,2(f7.2,f8.3,f7.2))'
      form2(0) = '(3f10.4,a3,i2,23x,(f7.2,f8.3,f7.2))'
      form3(0) = '(8e10.3)'
      form1(1) = '(g15.7,g10.4,g10.4,a3,i2,1x,2(g7.2,g8.3,g7.2))'
      form2(1) = '(g15.7,g10.4,g10.4,a3,i2,23x,(g7.2,g8.3,g7.2))'
      form3(1) = '(8g15.7)'

! Pass constants to common block names for rrtatm
      pic = pi
      planckc = planck
      boltzc = boltz
      clightc = clight
      avogadc = avogad
      alosmtc = alosmt
      gasconc = gascon
      radcn1c = radcn1
      radcn2c = radcn2

      ixmax = maxinpx

      ird = ird_in

      do l = 1, mxlay
         do ix = 1, maxxsec
           wx(ix,l) = 0.0_rb
         enddo
      enddo

! Top of read input loop
 1000 continue
      read (ird,9010,end=8800) ctest
      if (ctest .eq. cprcnt) goto 8900
      if (ctest .ne. cdollar) goto 1000

      read (ird,9011) iaer, iatm, ixsect, numangs, iout, idrv, imca, icld

!  If numangs set to -1, reset to default rt code for
!  backwards compatibility with original rrtm
      if (numangs .eq. -1) numangs = 0

!  If clouds are present, read in appropriate input file, IN_CLD_RRTM.
      if (icld .ge. 1) call readcld

!  If aerosols are present, read in appropriate input file, IN_AER_RRTM.
      if (iaer .eq. 10) call readaer

!  Read in surface information.
      read (ird,9012) tbound,iemiss,ireflect,semis(1:nbndlw)

!  Read in change in surface temperature for upward flux derivative adjustment
      dtbound = 0.0_rb
      if (idrv .eq. 1) then
         read (ird,9012) dtbound
      endif

!  Read in extra records for EXP or ER cloud overlap
      if (icld .eq. 4 .or. icld .eq. 5) then
         read (ird,9014) idcor
         if (idcor .eq. 0) then
            read (ird,9015) decorr_con
         endif
         if (idcor .eq. 1) then
            read (ird,9016) juldat, lat
            if (iatm .eq. 1) lat = ref_lat
         endif
      endif

      do iband = 1, nbndlw
         semiss(iband) = 1.0_rb
         if (iemiss .eq. 1 .and. semis(1) .ne. 0._rb) then
            semiss(iband) = semis(1)
         elseif (iemiss .eq. 2) then
            if (semis(iband) .ne. 0._rb) then
               semiss(iband) = semis(iband)
            endif
         endif
      enddo

      if (iatm .eq. 0) then
         read (ird,9013) iform,nlayers,nmol
         if (nmol.eq.0) nmol = 7
         read (ird,form1(iform)) pavel(1),tavel(1),secntk,cinp, &
              ipthak,altz(0),pz(0),tz(0),altz(1),pz(1),tz(1)
         read (ird,form3(iform)) (wkl(m,1),m=1,7), wbrodl(1)
         if(nmol .gt. 7) read (ird,form3(iform)) (wkl(m,1),m=8,nmol)

         do l = 2, nlayers
            read (ird,form2(iform)) pavel(l),tavel(l),secntk,cinp, &
                 ipthrk,altz(l),pz(l),tz(l)
            read (ird,form3(iform)) (wkl(m,l),m=1,7), wbrodl(l)
            if(nmol .gt. 7) read (ird,form3(iform)) (wkl(m,l),m=8,nmol)
         enddo

         if (ixsect .eq. 1) then
            read (ird,9300) nxmol0
            nxmol = nxmol0
            call xsident(ird)
            read (ird,9301) iformx

            do l = 1, nlayers
               read (ird,9010) cdum
               read (ird, form3(iformx)) (wx0(m,l),m=1,7),wbrodx
               if (nxmol0 .gt. 7) read (ird,form3(iformx)) (wx0(m,l),m=8,nxmol0)
            enddo
         endif
! calculate approximate layer thickness
         do l = 1,nlayers
            dz(l) = altz(l) - altz(l-1)
         enddo
      else
         ipu = 7
         ipr = 66
         open(unit=ipr,file='TAPE6',status='unknown')
         call rrtatm
         if (ixsect .eq. 1) then
            do mx = 1, nxmol0
               ixindx(mx) = ixtrans(ixindx0(mx))
            enddo
         endif
      endif
      if (tbound .lt. 0) tbound = tz(0)

!  Test for mixing ratio input.
      imix = 1
      do m = 1, nmol
         if (wkl(m,1) .gt. 1.0_rb) then
            imix = 0
            goto 3600
         endif
      enddo
 3600 continue

      if (ixsect .eq. 1) then
         imixx = 0
         if (wx0(1,1) .le. 1.0_rb) imixx = 1
      endif
      amttl = 0.0_rb
      wvttl = 0.0_rb
      do l = 1, nlayers
         summol = 0.0_rb
         do imol = 2, nmol
            summol = summol + wkl(imol,l)
         enddo
         if (imix .eq. 1) then
            coldry(l) = wbrodl(l) / (1._rb - summol)
            do imol = 1, nmol
               wkl(imol,l) = coldry(l) * wkl(imol,l)
            enddo
         else
            coldry(l) = wbrodl(l) + summol
         endif
         amttl = amttl + coldry(l)+wkl(1,l)
         wvttl = wvttl + wkl(1,l)
         if (ixsect .eq. 1) then
            do ix = 1, nxmol0
               if (ixindx(ix) .ne. 0) then
                  if (imixx .eq. 1) then
                     wx(ixindx(ix),l) = coldry(l) * wx0(ix,l) * 1.e-20_rb
                  else
                     wx(ixindx(ix),l) = wx0(ix,l) * 1.e-20_rb
                  endif
               endif
            enddo
         endif
      enddo

!  Calculate total precipitable water
      wvsh = (amw * wvttl) / (amd * amttl)
      pwvcm = wvsh * (1.e3_rb * pz(0)) / (1.e2_rb * grav)

! Pass output arrays to new variables for transfer to rrtmg through subroutine call.
      nlayers_out = nlayers
      iout_out = iout
      icld_out = icld
      iaer_out = iaer
      tbound_out = tbound
      dtbound_out = dtbound
      pwvcm_out = pwvcm
      inflag_out = inflag
      iceflag_out = iceflag
      liqflag_out = liqflag

      pz_out(0) = pz(0)
      tz_out(0) = tz(0)
      do l = 1, nlayers
         pavel_out(l) = pavel(l)
         tavel_out(l) = tavel(l)
! convert layer thickness from km to m
         dz_out(l) = dz(l) * 1000._rb
         pz_out(l) = pz(l)
         tz_out(l) = tz(l)
         coldry_out(l) = coldry(l)
         wbrodl_out(l) = wbrodl(l)
         cldfrac_out(l) = cldfrac(l)
         do m = 1, mxmol
            wkl_out(m,l) = wkl(m,l)
         enddo
         do ix = 1, maxxsec
            wx_out(ix,l) = wx(ix,l)
         enddo
      enddo
      do n = 1, nbndlw
         semiss_out(n) = semiss(n)
      enddo
      decorr_con_out = decorr_con
      juldat_out = juldat
      idcor_out = idcor
      lat_out = lat

      do l = 1, nlayers
         if (inflag.eq.0) then
            do n = 1, nbndlw
               tauc(n,l) = clddat1(l)
!               ssac(n,l) = 1._rb
!               asmc(n,l) = 1._rb
            enddo
            ciwp(l) = 0._rb
            clwp(l) = 0._rb
            fice(l) = 0._rb
            rei(l) = 0._rb
            rel(l) = 0._rb
         else
            do n = 1, nbndlw
               tauc(n,l) = 0._rb
!               ssac(n,l) = 1._rb
!               asmc(n,l) = 1._rb
            enddo
            cwp = clddat1(l)
            fice(l) = clddat2(l)
            ciwp(l) = cwp * fice(l)
            clwp(l) = cwp * (1._rb - fice(l))
            rei(l) = clddat3(l)
            rel(l) = clddat4(l)
         endif
      enddo

      do l = 1, nlayers
         do n = 1, nbndlw
            tauaer_out(l,n) = tauaer(l,n)
!            ssaaer_out(l,n) = 1._rb
!            asmaer_out(l,n) = 0._rb
         enddo
      enddo

      goto 9000

 8800 continue
 8900 if (ctest.eq.'%') stop 'END OF INPUT FILE'
 9000 continue

 9010 format (a1)
 9011 format (18x,i2,29x,i1,19x,i1,13x,i2,2x,i3,1x,i1,1x,i1,i1)
 9012 format (e10.3,1x,i1,2x,i1,16e5.3)
 9013 format (1x,i1,i3,i5)
 9014 format (8x,i2)
 9015 format (f10.3)
 9016 format (5x,i5,f10.3)
 9300 format (i5)
 9301 format (1x,i1)

      end subroutine readprof

!*************************************************************************
      subroutine readcld
!*************************************************************************

! --------- Modules ----------

      use parkind, only : im => kind_im, rb => kind_rb

! Purpose:  To read in IN_CLD_RRTM, the file that contains input
!           cloud properties.

      implicit integer(i-n), real (a-h,o-z)

! ------- Parameters -------
      parameter (mxlay=603)
      parameter (nbndlw = 16)

      common /profile/   nlayers,pavel(mxlay),tavel(mxlay),pz(0:mxlay),tz(0:mxlay)
      common /cloudin/   inflag,clddat1(mxlay),clddat2(mxlay), &
                         iceflag,liqflag,clddat3(mxlay),clddat4(mxlay)
      common /clouddat/  ncbands,cldfrac(mxlay),taucloud(mxlay,nbndlw)

      character*1 ctest, cpercent

      data cpercent /'%'/
      irdcld = 11

      open(irdcld,file='IN_CLD_RRTM',form='formatted')

! Read in cloud input option.
      read(irdcld,9050) inflag, iceflag, liqflag

      do lay = 1, nlayers
         cldfrac(lay) = 0._rb
      enddo

! Top of read input loop
 1000 continue

!  For INFLAG = 0 or 1, for each cloudy layer only LAY, FRAC, and
!  DAT1 are pertinent.  If CTEST = '%', then there are no more
!  cloudy layers to process.
      read (irdcld,9100,end=9000) ctest,lay,frac,dat1,dat2,dat3,dat4
      if (ctest .eq. cpercent) goto 9000
      cldfrac(lay) = frac
      clddat1(lay) = dat1
      clddat2(lay) = dat2
      clddat3(lay) = dat3
      clddat4(lay) = dat4
      goto 1000

 9000 continue
      close(irdcld)

 9050 format (3x,i2,4x,i1,4x,i1)
 9100 format (a1,1x,i3,5e10.5)

      end subroutine readcld

!***************************************************************************
      subroutine readaer
!***************************************************************************

      use parkind, only : im => kind_im, rb => kind_rb

! Purpose:  To read in IN_AER_RRTM, the file that contains input
!           aerosol properties.

! -------- Modules --------

      use rrlw_wvn, only : wavenum1, wavenum2

      implicit integer(i-n), real(a-h,o-z)

! ------- Parameters -------
      parameter (mxlay = 603)
      parameter (nbndlw  = 16)
!      parameter (mg = 16)
!      parameter (mxstr = 16)
!      parameter (mcmu = 32)

      real aod(mxlay),aod1(nbndlw)
      integer lay(mxlay),ivec(mxlay)

      common /control/  numangs, iout, istart, iend, icld, iaer
      common /profile/  nlayers,pavel(mxlay),tavel(mxlay),pz(0:mxlay),tz(0:mxlay)

      common /aerdat/  tauaer(mxlay,nbndlw)

      character*1 ctest, cpercent

      data cpercent /'%'/

      eps = 1.e-10_rb
      irdaer = 12
      open(irdaer,file='IN_AER_RRTM',form='formatted')

      aod(:) = 0.0_rb
      tauaer(:,:) = 0.0_rb

! Read in number of different aerosol models option.
      read (irdaer, 9010) naer
!       if (naer .gt. 4) then
!          print *, 'NAER (= ', naer, ') IS GREATER THAN 4'
!          stop
!       endif

! For each aerosol read in optical properties and layer aerosol
! optical depths.
      do ia = 1, naer
	 read (irdaer, 9011) nlay, iaod

! Input restricted to direct input of aerosol optical depths
         iaod = 1

! For this aerosol, read in layers and optical depth information.
! Store a nonzero optical depth in aod to check for double specification.
         do il = 1, nlay
            read(irdaer, 9012) lay(il), (aod1(ib), ib = 1,nbndlw)
            if (aod(lay(il)) .lt. eps) then
               if (iaod .eq. 1) then
                  do ib = 1, nbndlw
                     aod(lay(il)) = max(aod(lay(il)),aod1(ib))
                     tauaer(lay(il),ib) = aod1(ib)
                  enddo
               endif
            else
               print *,'LAYER ',lay(il),' HAS MORE THAN ONE AEROSOL TYPE'
               stop
            endif
         enddo

! End of naer loop
      enddo

 9000 continue
      close(irdaer)

 9010 format (3x, i2)
 9011 format (2x, i3, 4x, i1)
 9012 format (2x, i3, 16f7.4)

      end subroutine readaer

!**********************************************************************
      subroutine xsident(ird)
!**********************************************************************

! Purpose:  This subroutine identifies which cross-sections are to be used.

      implicit integer(i-n), real (a-h,o-z)

! ------- Parameters -------
      parameter (maxinpx=39)
      parameter (maxxsec=4)

      common /xsecctrl/ nxmol,ixindx(maxinpx)

!     nxmol     - number of cross-sections input by user
!     ixindx(i) - index of cross-section molecule corresponding to Ith
!                 cross-section specified by user
!                 = 0 -- not allowed in rrtm
!                 = 1 -- ccl4
!                 = 2 -- cfc11
!                 = 3 -- cfc12
!                 = 4 -- cfc22
!
!     xsname=names, alias=aliases of the cross-section molecules
!
      character*10 xsname(maxinpx),alias(maxxsec,4),blank

      data (alias(1,i),i=1,4)/ &
         'CCL4      ', 'CCL3F     ', 'CCL2F2    ', 'CHCLF2    '/
      data (alias(2,i),i=1,4)/ &
         ' ZZZZZZZZ ', 'CFCL3     ', 'CF2CL2    ', 'CHF2CL    '/
      data (alias(3,i),i=1,4)/ &
         ' ZZZZZZZZ ', 'CFC11     ', 'CFC12     ', 'CFC22     '/
      data (alias(4,i),i=1,4)/ &
          ' ZZZZZZZZ ', 'F11       ', 'F12       ', 'F22       '/

      data blank / '          '/

      do i = 1, nxmol
         xsname(i) = blank
      enddo

! Read in the names of the molecules

      if (nxmol.gt.7) then
         read (ird,'(7a10)') (xsname(i),i=1,7)
         read (ird,'(8a10)') (xsname(i),i=8,nxmol)
      else
         read (ird,'(7a10)') (xsname(i),i=1,nxmol)
      endif

!  Match the names read in against the names stored in alias
!  and determine the index value.
      ixmax = 4
      do i = 1, nxmol
!  Left-justify all inputed names.
         call cljust (xsname(i),10)
         ixindx(i) = 0
         do j = 1, ixmax
            if ((xsname(i).eq.alias(1,j)) .or. &
                (xsname(i).eq.alias(2,j)) .or. &
                (xsname(i).eq.alias(3,j)) .or. &
                (xsname(i).eq.alias(4,j))) then
               ixindx(i) = j
            endif
         enddo
      enddo

      end subroutine xsident

      end program rrtmg_lw


text:
      subroutine prteig(ev,ev2,u,period,nt,nf,sta,iounit,outname,
     &  id,iband,nev,nsta,ih,lstcc,stcor,declav,lhd,cvecname)

ccc   OLD text file output for multi- station eigenvector analysis
ccc   changed 1 Oct, 1996

      integer iband(2),ih(nsta)
      complex u(nt,nev)
      real ev(nt),ev2(nev),stcor(2,*),declav,tpower
      character*20 outname
      character*10 cvecname
      character*4 sta(*)
      character*80 cformat,chd
      logical lstcc,lhd
ccc   output eigenvalues,vectors

      cformat = '(1x,a3,1x,2f7.2,5(2f8.3,1x)/19x,5(2f8.3,1x))'
      chd = '(1x,4hsta.,3x,4hlat.,6x,4hlon.)'

      if(lhd) then
ccc      output title, eigenvalues
         write(iounit,601) outname,period
601      format('         array  ',a20,
     &         '   period =   ',e12.4,' seconds  ')

         write(iounit,*) '# of data vectors used = ', nf
         write(iounit,*) 'Decimation level  = ',id,'; Freqs. = ',iband

         if(nf.eq.0) return

         write(iounit,*) 'Eigenvalues (normalized to local noise):'
         write(iounit,615) (ev(k),k=1,nt)
         write(iounit,*) 'Fraction of Total Power (SNR units)'
         tev = 0.0
         do k = 1,nt
            tev = tev+ev(k)
         enddo
         write(iounit,614) (ev(k)/tev,k=1,nt)
         write(iounit,*) 'Power in PCs of Sig. Comp. (nT**2/Hz):'
         write(iounit,615) (ev2(k),k=1,nev)
         write(iounit,*) 'Fraction of Total Power (nT**2/Hz)'
         tev = 0.0
         do k = 1,nev
            tev = tev + ev2(k)
         enddo
         write(iounit,614) (ev2(k)/tev,k=1,nev)
615      format(8f10.2)
614      format(8f10.6)
      end if

ccc   write out vectors in U'
      write(iounit,*)
      do j=1,nev
         write(iounit,'(a10,i3)') cvecname,j
         write(iounit,chd)
         do i=1,nsta
            write(iounit,cformat) sta(i),stcor(1,i),stcor(2,i),
     &                  (u(k,j),k=ih(i),ih(i+1)-1)
         enddo
      enddo
      write(iounit,*) '_______________________________________',
     & '_______________________________________'
      return
      end

ext: f

max_stars_repo_name: simpeg-research/iris-mt-scratch

max_issues_repo_name: simpeg-research/iris-mt-scratch

max_forks_repo_name: simpeg-research/iris-mt-scratch


text:
!***************************************************************************************************
! xnet_timers.f90 10/18/17
! This file contains modules and subroutines for internal XNet timers
!***************************************************************************************************

Module xnet_timers
  !-------------------------------------------------------------------------------------------------
  ! This module contains the performance timers for XNet.
  !-------------------------------------------------------------------------------------------------
  Use xnet_types, Only: dp
  Implicit None
  Real(dp) :: timer_burner = 0.0  ! Total burner execution time
  Real(dp) :: timer_xnet   = 0.0  ! Total XNet execution time
  Real(dp) :: timer_setup  = 0.0  ! Data loading or preprocessing time
  Real(dp) :: timer_csect  = 0.0  ! Cross section calculation time
  Real(dp) :: timer_deriv  = 0.0  ! Derivative calculation time
  Real(dp) :: timer_jacob  = 0.0  ! Jacobian building time
  Real(dp) :: timer_decmp  = 0.0  ! LU Decomposition time
  Real(dp) :: timer_bksub  = 0.0  ! Backsubstitution time
  Real(dp) :: timer_nraph  = 0.0  ! Newton Raphson iteration timer
  Real(dp) :: timer_tstep  = 0.0  ! Time integration step timer
  Real(dp) :: timer_solve  = 0.0  ! Solution time
  Real(dp) :: timer_scrn   = 0.0  ! Screening and EOS time
  Real(dp) :: timer_eos    = 0.0  ! Screening and EOS time
  Real(dp) :: timer_output = 0.0  ! Output time
  Real(dp) :: start_timer  = 0.0  ! cpu time at the beginning of the timer block
  Real(dp) :: stop_timer   = 0.0  ! cpu time at the end of the timer block
  !$omp threadprivate(timer_burner,timer_xnet,timer_setup,timer_csect,timer_deriv,timer_jacob,timer_decmp, &
  !$omp   timer_bksub,timer_nraph,timer_tstep,timer_solve,timer_scrn,timer_eos,timer_output,start_timer,stop_timer)

Contains

  Function xnet_wtime()
    !-----------------------------------------------------------------------------------------------
    ! This function returns the wall time in a manner akin to omp_get_wtime().
    !-----------------------------------------------------------------------------------------------
    Use xnet_types, Only: i8
    Implicit None

    ! Function variable
    Real(dp) :: xnet_wtime

    ! Local variables
    Integer(i8) :: clock_read
    Integer(i8) :: clock_rate
    Integer(i8) :: clock_max

    Call system_clock(clock_read,clock_rate,clock_max)
    xnet_wtime = real(clock_read,dp) / real(clock_rate,dp)

    Return
  End Function xnet_wtime

  Subroutine reset_timers
    !-----------------------------------------------------------------------------------------------
    ! This routine resets timers for zone-independent timing.
    !-----------------------------------------------------------------------------------------------
    Implicit None
    timer_xnet  = 0.0
    timer_setup = 0.0
    timer_csect = 0.0
    timer_deriv = 0.0
    timer_jacob = 0.0
    timer_decmp = 0.0
    timer_bksub = 0.0
    timer_nraph = 0.0
    timer_tstep = 0.0
    timer_solve = 0.0
    timer_scrn  = 0.0
    timer_eos   = 0.0

    Return
  End Subroutine reset_timers

End Module xnet_timers

max_stars_repo_name: vishaltiwari/XNet

max_issues_repo_name: vishaltiwari/XNet

max_forks_repo_name: vishaltiwari/XNet


text:
      program cfgchk

c*********************************************************************
c
c     dl_poly utility program for checking a CONFIG file to ensure
c     an absence of close or zero atomic contacts
c
c     copyright daresbury laboratory 1998
c     author w.smith december 1998
c
c*********************************************************************

      implicit real*8 (a-h,o-z)
      parameter (mxatms=30000)

      character*8 name,filnam
      dimension cell(9),rcell(9)
      dimension xxx(mxatms),yyy(mxatms),zzz(mxatms)

      write(*,*)'Enter name of CONFIG file'
      read(*,*)filnam
      open(7,file=filnam)

      read(7,*)
      read(7,*)levcfg,imcon
      if(imcon.gt.0)then

        read(7,*)cell(1),cell(2),cell(3)
        read(7,*)cell(4),cell(5),cell(6)
        read(7,*)cell(7),cell(8),cell(9)
        call invert(cell,rcell,det)
        if(abs(det).lt.1.d-8)write(*,'(a)')
     x    'Error - zero determinant for cell vectors'

      endif

      do i=1,mxatms

        read(7,*,end=100)name
        read(7,*)xxx(i),yyy(i),zzz(i)
        if(levcfg.gt.0)read(7,*)vxx,vyy,vzz
        if(levcfg.gt.1)read(7,*)fxx,fyy,fzz

      enddo

  100 continue

      rchk=1.d10
      natms=i-1
      last=natms
      mlast=natms/2
      nlast=(natms-1)/2

      write(*,'(a)')'Checking pair separations.......'

      do m=1,mlast

        if(m.gt.nlast)last=mlast

        do i=1,last

          j=i+m
          if(j.gt.natms)j=j-natms

          ddx=xxx(i)-xxx(j)
          ddy=yyy(i)-yyy(j)
          ddz=zzz(i)-zzz(j)
          if(imcon.gt.0)then

            sdx=rcell(1)*ddx+rcell(4)*ddy+rcell(7)*ddz
            sdy=rcell(2)*ddx+rcell(5)*ddy+rcell(8)*ddz
            sdz=rcell(3)*ddx+rcell(6)*ddy+rcell(9)*ddz
            sdx=sdx-nint(sdx)
            sdy=sdy-nint(sdy)
            sdz=sdz-nint(sdz)
            ddx=cell(1)*sdx+cell(4)*sdy+cell(7)*sdz
            ddy=cell(2)*sdx+cell(5)*sdy+cell(8)*sdz
            ddz=cell(3)*sdx+cell(6)*sdy+cell(9)*sdz

          endif

          rsq=ddx**2+ddy**2+ddz**2

          rchk=min(rchk,rsq)

          if(rsq.lt.1.d-8)then

            write(*,'(i6,1p,3e12.4)')i,xxx(i),yyy(i),zzz(i)
            write(*,'(i6,1p,3e12.4)')j,xxx(j),yyy(j),zzz(j)
            write(*,'(1p,3e12.4)')rsq,sqrt(rsq)

          endif

        enddo

      enddo

      write(*,'(a)')'Checking complete.'
      write(*,'(a,1p,e12.4)')'Minimum separation found:',sqrt(rchk)

      end
      subroutine invert(a,b,d)

c***********************************************************************
c
c     dl_poly subroutine to invert a 3 * 3 matrix using cofactors
c
c     copyright - daresbury laboratory 1992
c     author    - w. smith       april 1992
c
c***********************************************************************

      implicit none

      real*8 a,b,d,r

      dimension a(9),b(9)

c     calculate adjoint matrix

      b(1)=a(5)*a(9)-a(6)*a(8)
      b(2)=a(3)*a(8)-a(2)*a(9)
      b(3)=a(2)*a(6)-a(3)*a(5)
      b(4)=a(6)*a(7)-a(4)*a(9)
      b(5)=a(1)*a(9)-a(3)*a(7)
      b(6)=a(3)*a(4)-a(1)*a(6)
      b(7)=a(4)*a(8)-a(5)*a(7)
      b(8)=a(2)*a(7)-a(1)*a(8)
      b(9)=a(1)*a(5)-a(2)*a(4)

c     calculate determinant

      d=a(1)*b(1)+a(4)*b(2)+a(7)*b(3)
      r=0.d0
      if(abs(d).gt.0.d0)r=1.d0/d

c     complete inverse matrix

      b(1)=r*b(1)
      b(2)=r*b(2)
      b(3)=r*b(3)
      b(4)=r*b(4)
      b(5)=r*b(5)
      b(6)=r*b(6)
      b(7)=r*b(7)
      b(8)=r*b(8)
      b(9)=r*b(9)

      return
      end


text:
      FUNCTION rofunc(b)
      INTEGER NMAX
      REAL rofunc,b,EPS
      PARAMETER (NMAX=1000,EPS=1.e-7)
CU    USES select
      INTEGER j,ndata
      REAL aa,abdev,d,sum,arr(NMAX),x(NMAX),y(NMAX),select
      COMMON /arrays/ x,y,arr,aa,abdev,ndata
      do 11 j=1,ndata
        arr(j)=y(j)-b*x(j)
11    continue
      if (mod(ndata,2).eq.0) then
        j=ndata/2
        aa=0.5*(select(j,ndata,arr)+select(j+1,ndata,arr))
      else
        aa=select((ndata+1)/2,ndata,arr)
      endif
      sum=0.
      abdev=0.
      do 12 j=1,ndata
        d=y(j)-(b*x(j)+aa)
        abdev=abdev+abs(d)
        if (y(j).ne.0.) d=d/abs(y(j))
        if (abs(d).gt.EPS) sum=sum+x(j)*sign(1.0,d)
12    continue
      rofunc=sum
      return
      END

ext: for

max_stars_repo_path: NR-Functions/Numerical Recipes- Example & Functions/Functions/rofunc.for


text:
      SUBROUTINE PSLASSQ( N, X, IX, JX, DESCX, INCX, SCALE, SUMSQ )
*
*  -- ScaLAPACK auxiliary routine (version 1.7) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     May 1, 1997
*
*     .. Scalar Arguments ..
      INTEGER            IX, INCX, JX, N
      REAL               SCALE, SUMSQ
*     ..
*     .. Array Arguments ..
      INTEGER            DESCX( * )
      REAL               X( * )
*     ..
*
*  Purpose
*  =======
*
*  PSLASSQ  returns the values  scl  and  smsq  such that
*
*     ( scl**2 )*smsq = x( 1 )**2 +...+ x( n )**2 + ( scale**2 )*sumsq,
*
*  where  x( i ) = sub( X ) = X( IX+(JX-1)*DESCX(M_)+(i-1)*INCX ).
*  The value of sumsq is assumed to be non-negative and scl returns the
*  value
*
*     scl = max( scale, abs( x( i ) ) ).
*
*  scale and sumsq must be supplied in SCALE and SUMSQ respectively.
*  SCALE and SUMSQ are overwritten by scl and ssq respectively.
*
*  The routine makes only one pass through the vector sub( X ).
*
*  Notes
*  =====
*
*  Each global data object is described by an associated description
*  vector.  This vector stores the information required to establish
*  the mapping between an object element and its corresponding process
*  and memory location.
*
*  Let A be a generic term for any 2D block cyclicly distributed array.
*  Such a global array has an associated description vector DESCA.
*  In the following comments, the character _ should be read as
*  "of the global array".
*
*  NOTATION        STORED IN      EXPLANATION
*  --------------- -------------- --------------------------------------
*  DTYPE_A(global) DESCA( DTYPE_ )The descriptor type.  In this case,
*                                 DTYPE_A = 1.
*  CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
*                                 the BLACS process grid A is distribu-
*                                 ted over. The context itself is glo-
*                                 bal, but the handle (the integer
*                                 value) may vary.
*  M_A    (global) DESCA( M_ )    The number of rows in the global
*                                 array A.
*  N_A    (global) DESCA( N_ )    The number of columns in the global
*                                 array A.
*  MB_A   (global) DESCA( MB_ )   The blocking factor used to distribute
*                                 the rows of the array.
*  NB_A   (global) DESCA( NB_ )   The blocking factor used to distribute
*                                 the columns of the array.
*  RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
*                                 row of the array A is distributed.
*  CSRC_A (global) DESCA( CSRC_ ) The process column over which the
*                                 first column of the array A is
*                                 distributed.
*  LLD_A  (local)  DESCA( LLD_ )  The leading dimension of the local
*                                 array.  LLD_A >= MAX(1,LOCr(M_A)).
*
*  Let K be the number of rows or columns of a distributed matrix,
*  and assume that its process grid has dimension p x q.
*  LOCr( K ) denotes the number of elements of K that a process
*  would receive if K were distributed over the p processes of its
*  process column.
*  Similarly, LOCc( K ) denotes the number of elements of K that a
*  process would receive if K were distributed over the q processes of
*  its process row.
*  The values of LOCr() and LOCc() may be determined via a call to the
*  ScaLAPACK tool function, NUMROC:
*          LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
*          LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
*  An upper bound for these quantities may be computed by:
*          LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
*          LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
*
*  Because vectors may be viewed as a subclass of matrices, a
*  distributed vector is considered to be a distributed matrix.
*
*  The result are only available in the scope of sub( X ), i.e if
*  sub( X ) is distributed along a process row, the correct results are
*  only available in this process row of the grid. Similarly if sub( X )
*  is distributed along a process column, the correct results are only
*  available in this process column of the grid.
*
*  Arguments
*  =========
*
*  N       (global input) INTEGER
*          The length of the distributed vector sub( X ).
*
*  X       (input) REAL
*          The vector for which a scaled sum of squares is computed.
*             x( i )  = X(IX+(JX-1)*M_X +(i-1)*INCX ), 1 <= i <= n.
*
*  IX      (global input) INTEGER
*          The row index in the global array X indicating the first
*          row of sub( X ).
*
*  JX      (global input) INTEGER
*          The column index in the global array X indicating the
*          first column of sub( X ).
*
*  DESCX   (global and local input) INTEGER array of dimension DLEN_.
*          The array descriptor for the distributed matrix X.
*
*  INCX    (global input) INTEGER
*          The global increment for the elements of X. Only two values
*          of INCX are supported in this version, namely 1 and M_X.
*          INCX must not be zero.
*
*  SCALE   (local input/local output) REAL
*          On entry, the value  scale  in the equation above.
*          On exit, SCALE is overwritten with  scl , the scaling factor
*          for the sum of squares.
*
*  SUMSQ   (local input/local output) REAL
*          On entry, the value  sumsq  in the equation above.
*          On exit, SUMSQ is overwritten with  smsq , the basic sum of
*          squares from which  scl  has been factored out.
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
     $                   LLD_, MB_, M_, NB_, N_, RSRC_
      PARAMETER          ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1,
     $                     CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6,
     $                     RSRC_ = 7, CSRC_ = 8, LLD_ = 9 )
      REAL               ZERO
      PARAMETER          ( ZERO = 0.0E+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I, ICOFF, ICTXT, IIX, IOFF, IROFF, IXCOL,
     $                   IXROW, JJX, LDX, MYCOL, MYROW, NP, NPCOL,
     $                   NPROW, NQ
      REAL               TEMP1
*     ..
*     .. Local Arrays ..
      REAL               WORK( 2 )
*     ..
*     .. External Subroutines ..
      EXTERNAL           BLACS_GRIDINFO, INFOG2L, PSTREECOMB, SCOMBSSQ
*     ..
*     .. External Functions ..
      INTEGER            NUMROC
      EXTERNAL           NUMROC
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABS, MOD
*     ..
*     .. Executable Statements ..
*
*     Get grid parameters.
*
      ICTXT = DESCX( CTXT_ )
      CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )
*
*     Figure local indexes
*
      CALL INFOG2L( IX, JX, DESCX, NPROW, NPCOL, MYROW, MYCOL, IIX, JJX,
     $              IXROW, IXCOL )
*
      LDX = DESCX( LLD_ )
      IF( INCX.EQ.DESCX( M_ ) ) THEN
*
*        X is rowwise distributed.
*
         IF( MYROW.NE.IXROW )
     $      RETURN
         ICOFF = MOD( JX, DESCX( NB_ ) )
         NQ = NUMROC( N+ICOFF, DESCX( NB_ ), MYCOL, IXCOL, NPCOL )
         IF( MYCOL.EQ.IXCOL )
     $      NQ = NQ - ICOFF
*
*        Code direct from LAPACK's SLASSQ, (save subroutine call)
*
         IF( NQ.GT.0 ) THEN
            IOFF = IIX + ( JJX - 1 ) * LDX
            DO 10 I = 1, NQ
               IF( X( IOFF ).NE.ZERO ) THEN
                  TEMP1 = ABS( X( IOFF ) )
                  IF( SCALE.LT.TEMP1 ) THEN
                     SUMSQ = 1 + SUMSQ * ( SCALE / TEMP1 )**2
                     SCALE = TEMP1
                  ELSE
                     SUMSQ = SUMSQ + ( TEMP1 / SCALE )**2
                  END IF
               END IF
               IOFF = IOFF + LDX
   10       CONTINUE
         END IF
*
*        Take local result and find global
*
         WORK( 1 ) = SCALE
         WORK( 2 ) = SUMSQ
*
         CALL PSTREECOMB( ICTXT, 'Rowwise', 2, WORK, -1, IXCOL,
     $                    SCOMBSSQ )
*
         SCALE = WORK( 1 )
         SUMSQ = WORK( 2 )
*
      ELSE IF( INCX.EQ.1 ) THEN
*
*        X is columnwise distributed.
*
         IF( MYCOL.NE.IXCOL )
     $      RETURN
         IROFF = MOD( IX, DESCX( MB_ ) )
         NP = NUMROC( N+IROFF, DESCX( MB_ ), MYROW, IXROW, NPROW )
         IF( MYROW.EQ.IXROW )
     $      NP = NP - IROFF
*
*        Code direct from LAPACK's SLASSQ, (save subroutine call)
*
         IF( NP.GT.0 ) THEN
            IOFF = IIX + ( JJX - 1 ) * LDX
            DO 20 I = 1, NP
               IF( X( IOFF ).NE.ZERO ) THEN
                  TEMP1 = ABS( X( IOFF ) )
                  IF( SCALE.LT.TEMP1 ) THEN
                     SUMSQ = 1 + SUMSQ*( SCALE / TEMP1 )**2
                     SCALE = TEMP1
                  ELSE
                     SUMSQ = SUMSQ + ( TEMP1 / SCALE )**2
                  END IF
               END IF
               IOFF = IOFF + 1
   20       CONTINUE
         END IF
*
*        Take local result and find global
*
         WORK( 1 ) = SCALE
         WORK( 2 ) = SUMSQ
*
         CALL PSTREECOMB( ICTXT, 'Columnwise', 2, WORK, -1, IXCOL,
     $                    SCOMBSSQ )
*
         SCALE = WORK( 1 )
         SUMSQ = WORK( 2 )
*
      END IF
*
      RETURN
*
*     End of PSLASSQ
*
      END

max_stars_repo_name: jschueller/scalapack

max_issues_repo_name: jschueller/scalapack

max_forks_repo_name: jschueller/scalapack


SUBROUTINE p103(input_file,output_file)
!-------------------------------------------------------------------------
! Program 10.3: Eigenvalue analysis of an elastic solid in plane strain
!               using 4-node rectangular quadrilaterals. Lanczos Method.
!               Consistent mass. Mesh numbered in y-direction.
!-------------------------------------------------------------------------
 USE main
 USE geom
 IMPLICIT NONE
 CHARACTER(len=60),INTENT(IN) :: input_file
 CHARACTER(len=60),INTENT(OUT) :: output_file
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 INTEGER::i,iel,iflag=-1,iters,jflag,k,lalfa,leig,lp=6,lx,lz,nband=0,     &
   ndim=2,ndof,neig=0,nels,neq,nip=4,nlen,nmodes,nn,nod,nodof=2,nprops=3, &
   np_types,nr,nst=3,nxe,nye
 REAL(iwp)::acc,det,el,er,zero=0.0_iwp
 CHARACTER(LEN=15)::argv,element='quadrilateral'
!----------------------------- dynamic arrays-----------------------------
 INTEGER,ALLOCATABLE::etype(:),g(:),g_g(:,:),g_num(:,:),jeig(:,:),nf(:,:),&
   nu(:),num(:)
 REAL(iwp),ALLOCATABLE::alfa(:),bee(:,:),beta(:),coord(:,:),dee(:,:),     &
   del(:),der(:,:),deriv(:,:),diag(:),ecm(:,:),eig(:),fun(:),g_coord(:,:),&
   jac(:,:),kb(:,:),km(:,:),mb(:,:),mm(:,:),points(:,:),prop(:,:),ua(:),  &
   udiag(:),va(:),v_store(:,:),weights(:),w1(:),x(:),x_coords(:),y(:,:),  &
   y_coords(:),z(:,:)
!---------------------------input and initialisation----------------------

 OPEN(10,FILE=input_file)
 OPEN(11,FILE=output_file)
 READ(10,*)nxe,nye,nod,np_types
 CALL mesh_size(element,nod,nels,nn,nxe,nye)
 ndof=nod*nodof
 ALLOCATE(nf(nodof,nn),points(nip,ndim),dee(nst,nst),g_coord(ndim,nn),    &
   coord(nod,ndim),fun(nod),jac(ndim,ndim),weights(nip),g_num(nod,nels),  &
   der(ndim,nod),deriv(ndim,nod),bee(nst,ndof),num(nod),km(ndof,ndof),    &
   g(ndof),g_g(ndof,nels),mm(ndof,ndof),ecm(ndof,ndof),                   &
   prop(nprops,np_types),x_coords(nxe+1),y_coords(nye+1),etype(nels))
 READ(10,*)prop
 etype=1
 IF(np_types>1)READ(10,*)etype
 READ(10,*)x_coords,y_coords
 nf=1
 READ(10,*)nr,(k,nf(:,k),i=1,nr)
 CALL formnf(nf)
 neq=MAXVAL(nf)
 READ(10,*)nmodes,el,er,lalfa,leig,lx,lz,acc
 ALLOCATE(eig(leig),x(lx),del(lx),nu(lx),jeig(2,leig),alfa(lalfa),        &
   beta(lalfa),z(lz,leig))
!-------- loop the elements to find nband and set up global arrays -------
 elements_1: DO iel=1,nels
   CALL geom_rect(element,iel,x_coords,y_coords,coord,num,'y')
   CALL num_to_g(num,nf,g)
   g_num(:,iel)=num
   g_coord(:,num)=TRANSPOSE(coord)
   g_g(:,iel)=g
   IF(nband<bandwidth(g))nband=bandwidth(g)
 END DO elements_1
 WRITE(11,'(A,I5,A,/,A,I5,/,A,I5)')" There are",neq," equations",         &
   " The half-bandwidth (including diagonal) is",nband+1
 ALLOCATE(kb(neq,nband+1),mb(neq,nband+1),ua(0:neq),va(0:neq),            &
   diag(0:neq),udiag(0:neq),w1(0:neq),y(0:neq,leig),v_store(0:neq,lalfa))
 kb=zero
 mb=zero
 ua=zero
 va=zero
 eig=zero
 jeig=0
 x=zero
 del=zero
 nu=0
 alfa=zero
 beta=zero
 diag=zero
 udiag=zero
 w1=zero
 y=zero
 z=zero
 CALL sample(element,points,weights)
!----------------- element stiffness integration and assembly-------------
 elements_2: DO iel=1,nels
   CALL deemat(dee,prop(1,etype(iel)),prop(2,etype(iel)))
   num=g_num(:,iel)
   coord=TRANSPOSE(g_coord(:,num))
   g=g_g(:,iel)
   km=zero
   mm=zero
   integrating_pts_1: DO i=1,nip
     CALL shape_fun(fun,points,i)
     CALL shape_der(der,points,i)
     jac=MATMUL(der,coord)
     det=determinant(jac)
     CALL invert(jac)
     deriv=MATMUL(jac,der)
     call beemat(bee,deriv)
     km=km+MATMUL(MATMUL(TRANSPOSE(bee),dee),bee)*det*weights(i)
     call ecmat(ecm,fun,ndof,nodof)
     mm=mm+ecm*det*weights(i)*prop(3,etype(iel))
   END DO integrating_pts_1
   CALL formkb(kb,km,g)
   CALL formkb(mb,mm,g)
 END DO elements_2
 CALL cholin(mb)
!------------------------------find eigenvalues---------------------------
 DO iters=1,lalfa
   CALL lancz1(neq,el,er,acc,leig,lx,lalfa,lp,iflag,ua,va,eig,jeig,neig,x,&
     del,nu,alfa,beta,v_store)
   IF(iflag==0)EXIT
   IF(iflag>1)THEN
     WRITE(11,'(A,I5)')                                                   &
       " Lancz1 is signalling failure, with iflag = ",   iflag
     STOP
   END IF
!--- iflag = 1 therefore form u + a * v  ( candidate for ebe ) -----------
   udiag=va
   CALL chobk2(mb,udiag)
   CALL banmul(kb,udiag,diag)
   CALL chobk1(mb,diag)
   ua=ua+diag
 END DO
!--- iflag = 0 therefore write out the spectrum  -------------------------
 WRITE(11,'(A,I5,A/)')" It took",iters," iterations"
 WRITE(11,'(3(A,E12.4))')" Eigenvalues in the range",el," and",er," are:"
 WRITE(11,'(6E12.4)')eig(1:neig)
!------------------- calculate the eigenvectors --------------------------
 IF(neig>10)neig=10
  CALL lancz2(neq,lalfa,lp,eig,jeig,neig,alfa,beta,lz,jflag,y,w1,z,v_store)
!-------------------if jflag is zero  calculate the eigenvectors ---------
 IF(jflag==0)THEN
   DO i=1,nmodes
     udiag(:)=y(:,i)
     CALL chobk2(mb,udiag)
     WRITE(11,'(" Eigenvector number",I4," is:")')i
     WRITE(11,'(6E12.4)')udiag(1:)/MAXVAL(ABS(udiag(1:)))
   END DO
 ELSE
! lancz2 fails
   WRITE(11,'(A,I5)')" Lancz2 is signalling failure with jflag = ",jflag
 END IF
END SUBROUTINE p103