text:
     subroutine get_clr_src_dir(solalt,viewalt,od_g_msl,od_g_vert,od_a_vert,ext_haero,htmsl,ssa,agv,ao,aav,aod_ref,redp_lvl,scale_ht_a,ags,aas,ic,idebug,istart,iend,srcdir,sumi_g,sumi_a,opac_slant,nsteps,ds,tausum_a)

     use mem_namelist, ONLY: earth_radius
     use constants_laps, ONLY: R, GRAV

     include 'trigd.inc'
     include 'rad.inc'

!    Input idebug for select angles

     trans(od) = exp(-min(od,80.))
     opac(od) = 1.0 - trans(od)
     scurve(x) = (-0.5 * cos(x*3.14159265)) + 0.5
     real*8 dsdst,dradius_start,daltray,dcurvat2
     dcurvat2(dsdst,dradius_start,daltray) = &
              sqrt(dsdst**2 + dradius_start**2 &
       - (2D0*dsdst*dradius_start*cosd(90D0+daltray))) - dradius_start

     real tausum_a(istart:iend)
     real*8 sbar,htbar

     real agv      ! gas   relative to zenith value for observer ht/viewalt
     real ao       ! ozone relative to zenith at sea level pressure/viewalt
     real aav      ! aero  relative to zenith value for observer ht/viewalt
     real ags      ! gas   relative to zenith value for observer ht/solalt
     real aos      ! ozone relative to zenith value for observer ht/solalt
     real aas      ! aero  relative to zenith value for observer ht/solalt

     od_o_msl = (o3_du/300.) * ext_o(ic)
     od_o_vert = od_o_msl * patm_o3(htmsl)

!    htmin_view = htminf(htmsl,viewalt,earth_radius)
!    patm_htmsl = ztopsa(htmin_view) / 1013.25
     patm_htmsl = ztopsa(htmsl) / 1013.25
     rearg = 2.0 + 0.667 * patm_htmsl

!    Calculate src term of gas+aerosol along ray normalized by TOA radiance
     od_g_slant = od_g_vert * agv
     od_o_slant = od_o_msl  * ao
     od_a_slant = od_a_vert * aav
     od_slant = od_g_slant + od_o_slant + od_a_slant
     od_slant = max(od_slant,1e-30) ! QC check
     r_e = 6371.e3

     scale_ht_g = R * ztotsa(htmsl) / GRAV
     alpha_sfc_g = od_g_vert / scale_ht_g

     scale_ht_g = R * ztotsa(0.) / GRAV
     alpha_ref_g = od_g_msl  / scale_ht_g

     alpha_sfc_a = od_a_vert / scale_ht_a
     alpha_ref_a = aod_ref / scale_ht_a

     opac_vert = opac(od_g_vert + od_o_vert + od_a_vert)
     if(opac_vert .eq. 0)then ! very small od
         opac_vert = od_g_vert + od_o_vert + od_a_vert
     endif

     opac_slant = opac(od_slant)
     if(opac_slant .eq. 0)then ! very small od
         opac_slant = od_slant
     endif

!    Integral [0 to x] a * exp(-b*x) = a * (1. - exp*(-b*x)) / b
!    ds = 25.

     if(idebug .eq. 1)then
         write(6,*)
         write(6,*)'solalt/viewalt = ',solalt,viewalt
         write(6,*)'ags/aas = ',ags,aas
         write(6,*)'agv/aav = ',agv,aav
         write(6,*)'od_g_vert = ',od_g_vert
         write(6,*)'scale_ht_g = ',scale_ht_g
         write(6,*)'od_o_msl/vert = ',od_o_msl,od_o_vert
         write(6,*)'od_o_vert = ',od_o_vert
         write(6,*)'od_a_vert = ',od_a_vert
         write(6,*)'opac_vert = ',opac_vert
         write(6,*)'od_g_slant = ',od_g_slant
         write(6,*)'od_o_slant,ao = ',od_o_slant,ao
!        if(viewalt .ge. 0.)then
           write(6,*)'aorel = ',od_o_slant/od_o_vert
!        endif
         write(6,*)'od_a_slant = ',od_a_slant
         write(6,*)'opac_slant = ',opac_slant
         write(6,*)'od_slant = ',od_slant
         write(6,*)'ds/nsteps = ',ds,nsteps
     endif ! i

     tausum = 0.
     tausum_a(:) = 0.
     sumi = 0. ; sumi_g = 0.; sumi_a = 0.
     frac_opac = 0.

     httopill = 100000.

     if(idebug .eq. 1)then
       write(6,*)'           sbar         htbar       dtaug    dtauo   tausum od_sol_a  od_solar   rad       di     sumi_g   sumi_a opac_curr frac_opac sumi_mn sumi_ext'
     endif

     opac_curr = 0.

     if(solalt .ge. 0.)then
         frac_iso_a = 1.0
     else
         frac_iso_a = 0.5
     endif

     do i = istart,iend
         sbar = (float(i)-0.5) * ds
         if(htmsl .le. 1000e3)then ! low altitude
           htbar = sbar * sind(viewalt)
           xybar = sbar * cosd(viewalt)
!          htbar = htbar + curvat(xybar,earth_radius*(4./3.)) ! Earth Curvature
           htbar = htbar + xybar**2 / (rearg * r_e) ! Earth Curvature
         else                      ! high altitude
           radius_start_eff = earth_radius + htmsl
!          Is this showing up negative at nadir?
           htbar = dcurvat2(sbar,dble(radius_start_eff),dble(viewalt))
         endif
         htbar_msl = htbar+htmsl

         zapp = 90. - solalt
         od_solar_slant_o3 = od_o_msl * patm_o3(htbar_msl) * airmasso(zapp,htbar_msl)

         patm_bar = ztopsa(htbar_msl) / 1013.25
         scale_ht_g = R * ztotsa(htbar_msl) / GRAV
         alphabar_g = (od_g_msl * patm_bar) / scale_ht_g

!        od_solar_slant_g = od_g_vert * ags * exp(-htbar/scale_ht_g)
         od_solar_slant_g = od_g_msl  * ags * exp(-htbar_msl/scale_ht_g)
!        od_solar_slant_a = od_a_vert * aas * exp(-htbar/scale_ht_a)
         od_solar_slant_a = aod_ref   * aas * exp(-(htbar_msl-redp_lvl)/scale_ht_a) * frac_iso_a
         od_solar_slant = od_solar_slant_g &
                        + od_solar_slant_o3 + od_solar_slant_a

         alphabar_o = alpha_o3(od_o_msl,htbar_msl)
!        alphabar_a = alpha_sfc_a * exp(-htbar/scale_ht_a)
         alphabar_a = alpha_ref_a * exp(-(htbar_msl-redp_lvl)/scale_ht_a)

!        AOD stratospheric = .012 (nighttime - layer extends above observer)
         if(htbar_msl .gt. h1_ha .and. htbar_msl .le. h2_ha .and. solalt .lt. 333. .and. aod_ha .gt. 0.)then
             alphabar_a = alphabar_a + alpha_ha * ext_haero
         elseif(htbar_msl .gt. h2_ha .and. htbar_msl .lt. h3_ha)then
             frac_ha = (h3_ha - htbar_msl) / (h3_ha - h2_ha)
             alphabar_a = alphabar_a + alpha_ha * ext_haero * scurve(frac_ha)
         endif

         dtau = ds * (alphabar_g + alphabar_o + alphabar_a)
         tausum = tausum + dtau
!        if(tausum .lt. 1.0)distod1 = sbar
         tausum_a(i) = tausum

         if(htbar .gt. httopill .or. htbar .lt. -(htmsl+500.) .or. frac_opac .gt. 11.4)then
           if(idebug .eq. 1)then
!            write(6,*)' criteria ',htbar,httopill,-(htmsl+500.),frac_opac
             write(6,11)i,sbar,htbar,alphabar_g*ds,alphabar_o*ds,tausum,od_solar_slant_a,od_solar_slant,rad,di,sumi_g,sumi_a,opac_curr,frac_opac,sumi_mean,sumi_extrap
           endif
           goto 900
         endif

         opac_last = opac_curr
         opac_curr = opac(tausum)
         if(opac_curr .eq. 0.)opac_curr = tausum ! very small od
         do = opac_curr - opac_last

!        Experiment with in tandem with 'skyglow_phys'
!        rad = 1.0
         rad = trans(od_solar_slant)
         radg = max(rad,.0001) ! secondary scattering
         rada = max(rad,.00001) ! secondary scattering
         rad = max(rad,.0001) ! secondary scattering

!        di = (dtau * rad) * exp(-tausum)
         if(dtau .gt. 0.)then
             ssa_eff = (alphabar_g + ssa * alphabar_a) / (alphabar_g + alphabar_a)
         else
             ssa_eff = 1.0
         endif

         alphabar = alphabar_g + alphabar_o + alphabar_a
         if(alphabar .gt. 0.)then
             di_g = do * radg       * (alphabar_g / alphabar)
             di_a = do * rada * ssa * (alphabar_a / alphabar)
         else
             di_g = 0.
             di_a = 0.
         endif
         di = di_g + di_a

         sumi_g = sumi_g + di_g
         sumi_a = sumi_a + di_a
         sumi = sumi + di

         frac_opac = opac_curr / opac_slant

         sumi_mean = sumi / opac_curr
!        sumi_extrap = (sumi + (1.-frac_opac) * rad) / opac_slant
         sumi_extrap = (sumi_mean * frac_opac) + (1.-frac_opac) * rad

         if(idebug .eq. 1)then
           if(i .eq. 10 .or. i .eq. 100 .or. i .eq. 1000 .or. i .eq. 10000 .or. i .eq. iend)then
             write(6,11)i,sbar,htbar,alphabar_g*ds,alphabar_o*ds,tausum,od_solar_slant_a,od_solar_slant,rad,di,sumi_g,sumi_a,opac_curr,frac_opac,sumi_mean,sumi_extrap
11           format(i8,f11.0,f13.1,2f10.7,f8.4,f9.6,f9.4,2f9.5,6f9.4)
           endif
         endif
     enddo ! i

900  continue

     srcdir = sumi_extrap

     return
     end

     subroutine get_clr_src_dir_low(solalt,solazi,viewalt,viewazi &      ! I
           ,od_g_msl,od_g_vert,od_o_msl,od_o_vert                 &      ! I
           ,od_a_vert,ext_haero,htmsl                             &      ! I
           ,ssa,agv,ao,aav,aod_ref,redp_lvl,scale_ht_a &                 ! I
           ,ags_in,aas_in,ags_a,aas_a,isolalt_lo,isolalt_hi,del_solalt & ! I
           ,ic,idebug,refdist_solalt,solalt_ref &                        ! I
           ,srcdir,sumi_g,sumi_a,opac_slant,nsteps,ds,tausum_a)

     use mem_namelist, ONLY: earth_radius
     use constants_laps, ONLY: R, GRAV

     include 'trigd.inc'
     include 'rad.inc'

     real*8 opac_last,opac_curr,tausum,trans,opac,od
     real logint

     trans(od) = exp(-min(od,80d0))
     opac(od) = 1d0 - trans(od)
     angdif(X,Y)=MOD(X-Y+540.,360.)-180.
     scurve(x) = (-0.5 * cos(x*3.14159265)) + 0.5
     logint(y1,y2,f) = exp(log(y1)*(1.-f) + log(y2)*f)
     curvat2(sdst,radius_start,altray) = &
             sqrt(sdst**2 + radius_start**2 &
               - (2.*sdst*radius_start*cosd(90.+altray))) - radius_start
     real*8 dsdst,dradius_start,daltray,dcurvat2
     dcurvat2(dsdst,dradius_start,daltray) = &
              sqrt(dsdst**2 + dradius_start**2 &
       - (2D0*dsdst*dradius_start*cosd(90D0+daltray))) - dradius_start
     expbl(x) = exp(max(x,-80.))

     real tausum_a(nsteps)

     real agv      ! gas   relative to zenith value for observer ht/viewalt
     real ao       ! ozone relative to zenith at sea level pressure/viewalt
     real aav      ! aero  relative to zenith value for observer ht/viewalt
     real ah       ! haero relative to zenith value for observer ht/viewalt
     real ags_in   ! gas   relative to zenith value for observer ht/solalt
     real aas_in   ! aero  relative to zenith value for observer ht/solalt
     real ags_a(isolalt_lo:isolalt_hi) ! gas rel to zen val fr obs ht/salt
     real aas_a(isolalt_lo:isolalt_hi) ! aero rel to zen val fr obs ht/salt

!    return ! timing

!    od_o_msl = (o3_du/300.) * ext_o(ic)
     od_o_vert = od_o_msl * patm_o3(htmsl)

!    htmin_view = htminf(htmsl,viewalt,earth_radius)
!    patm_htmsl = ztopsa(htmin_view) / 1013.25
     patm_htmsl = ztopsa(htmsl) / 1013.25
     rearg = 2.0 + 0.667 * patm_htmsl
     radius_earth_eff = earth_radius * (1. + (0.33333 * patm_htmsl))
     radius_start_eff = radius_earth_eff + htmsl

!    High aerosols (just during low sun)
     if(solalt .lt. 333.)then
         od_h_vert = aod_ha * patm_ah(htmsl,h1_ha,h2_ha)
     else
         od_h_vert = 0.
     endif
     if(patm_ah(htmsl,h1_ha,h2_ha) .eq. 1.)then
         ah = am_homo_lyr(90.-viewalt,h1_ha-htmsl,h2_ha-htmsl,earth_radius)
     elseif(patm_ah(htmsl,h1_ha,h2_ha) .eq. 0.)then ! above layer
         ah = 0.
     else ! in layer
         ah = am_homo_lyr(90.-viewalt,0.,h2_ha-htmsl,earth_radius)
     endif

!    Calculate src term of gas+aerosol along ray normalized by TOA radiance
     od_g_slant = od_g_vert * agv
     od_o_slant = od_o_msl  * ao
     od_a_slant = od_a_vert * aav
     od_h_slant = od_h_vert * ah
     od_slant = od_g_slant + od_o_slant + od_a_slant + od_h_slant
     od_slant = max(od_slant,1e-30) ! QC check
     r_e = 6371.e3

     scale_ht_g = R * ztotsa(htmsl) / GRAV
     alpha_sfc_g = od_g_vert / scale_ht_g

     scale_ht_g = R * ztotsa(0.) / GRAV
     alpha_ref_g = od_g_msl  / scale_ht_g

!    scale_ht_a = 1500. ! Input this?
!    alpha_sfc_a = od_a_vert / scale_ht_a
     alpha_ref_a = aod_ref / scale_ht_a

     opac_vert = opac(dble(od_g_vert + od_o_vert + od_a_vert))
     if(opac_vert .eq. 0)then ! very small od
         opac_vert = od_g_vert + od_o_vert + od_a_vert
     endif

     opac_slant = opac(dble(od_slant))
     if(opac_slant .eq. 0)then ! very small od
         opac_slant = od_slant
     endif

!    Integral [0 to x] a * exp(-b*x) = a * (1. - exp*(-b*x)) / b
     ds = 100. / max(sind(viewalt),0.1)
     istart = max(1,nint((htmsl-100000.)/ds))

     dsolalt_ref = solalt_ref - solalt

     if(idebug .eq. 1)then
         write(6,*)
         write(6,*)'subroutine get_clr_src_dir_low'
         write(6,*)'solalt/ref/solazi = ',solalt,solalt_ref,solazi
         write(6,*)'viewalt/viewazi = ',viewalt,viewazi
         write(6,*)'ags/aas = ',ags_in,aas_in
         write(6,*)'agv/aav = ',agv,aav
         write(6,*)'od_g_vert/od_g_msl = ',od_g_vert,od_g_msl
         write(6,*)'scale_ht_g = ',scale_ht_g
         write(6,*)'alpha_sfc_g/alpha_ref_g = ',alpha_sfc_g,alpha_ref_g
         write(6,*)'od_o_msl/vert = ',od_o_msl,od_o_vert
         write(6,*)'od_a_vert = ',od_a_vert
         write(6,*)'alpha_ha / ext_haero = ',alpha_ha,ext_haero
         write(6,*)'opac_vert = ',opac_vert
         write(6,*)'od_g_slant = ',od_g_slant
         write(6,*)'od_o_slant/ao = ',od_o_slant,ao
         write(6,*)'od_a_slant = ',od_a_slant
         write(6,*)'od_h_vert/slant = ',od_h_vert,od_h_slant
         write(6,*)'opac_slant = ',opac_slant
         write(6,*)'od_slant = ',od_slant
         write(6,*)'istart/ds = ',istart,ds
     endif ! i

!    ags_a = ags_in; aas_a = aas_in ! test
     tausum = 0.
     sumi = 0. ; sumi_g = 0.; sumi_a = 0.
     frac_opac = 0.

     solaltill = max(solalt,0.)
     httopill = max(100000.-sind(solaltill)*60000.,40000.)
     htbotill = -(htmsl + 500.)

     if(idebug .eq. 1)then
       write(6,*)'           sbar    htbar_msl   dtau_g   dtau_a    tausum  dsolalt      ags     aas   od_sol_g  od_solar  rad        di      sumi_g    sumi_a opac_curr frac_opac sumi_mn sumi_ext'
     endif

     opac_curr = 0.

     dsolalt_dxy = cosd(angdif(solazi,viewazi)) / 110000.

     if(solalt .ge. 0.)then
         frac_iso_a = 1.0
     else
         frac_iso_a = 0.5
     endif

     iend = max(istart+50000,nsteps)

     do i = istart,iend ! nsteps
         sbar = (float(i)-0.5) * ds

!        htbar = sbar * sind(viewalt)            ! Works at low altitude
!        htbar = htbar + sbar**2 / (rearg * r_e) ! Works at LEO
!        htbar = curvat2(sbar,radius_start_eff,viewalt) ! Works at GEO
!        Works from moon
         htbar = dcurvat2(dble(sbar),dble(radius_start_eff),dble(viewalt))

         htbar_msl = htbar+htmsl
         xybar = (sbar - refdist_solalt) * cosd(viewalt)

         patm_bar = ztopsa(htbar_msl) / 1013.25
         scale_ht_g = R * ztotsa(htbar_msl) / GRAV
         alphabar_g = (od_g_msl * patm_bar) / scale_ht_g

         alphabar_o = alpha_o3(od_o_msl,htbar_msl)
!        alphabar_a = alpha_sfc_a * exp(-htbar/scale_ht_a)
         alphabar_a = alpha_ref_a * expbl(-(htbar_msl-redp_lvl)/scale_ht_a)

!        AOD stratospheric = .012 (nighttime - layer extends above observer)
         if(htbar_msl .gt. h1_ha .and. htbar_msl .le. h2_ha .and. solalt .lt. 333. .and. aod_ha .gt. 0.)then
             alphabar_a = alphabar_a + alpha_ha * ext_haero
         elseif(htbar_msl .gt. h2_ha .and. htbar_msl .lt. h3_ha)then
             frac_ha = (h3_ha - htbar_msl) / (h3_ha - h2_ha)
             alphabar_a = alphabar_a + alpha_ha * ext_haero * scurve(frac_ha)
         endif

         dtau = ds * (alphabar_g + alphabar_o + alphabar_a)
         tausum = tausum + dtau
!        if(tausum .lt. 1.0)distod1 = sbar
!        tausum_a(i) = tausum

         if(htbar .gt. httopill .or. htbar .lt. htbotill .or. frac_opac .gt. 11.4 .or. (sbar.gt.refdist_solalt .and. htbar_msl.gt.130e3))then ! efficiency
!          tausum_a(min(i+1,nsteps):nsteps) = tausum
           goto 900
         endif

         opac_last = opac_curr
         opac_curr = opac(tausum)
         do = opac_curr - opac_last

!        Update ags,aas
         dsolalt = dsolalt_ref + dsolalt_dxy * xybar
         dsolalt = min(max(dsolalt,-90.),+90.)
         solalt_step = solalt + dsolalt
         dsolaltb = &
            max(min(dsolalt/del_solalt,float(isolalt_hi)/del_solalt-.0001),float(isolalt_lo)/del_solalt)
         isolaltl = floor(dsolaltb); isolalth = isolaltl + 1
         fsolalt = dsolaltb - float(isolaltl)
!        ags = ags_a(isolaltl) * (1.-fsolalt) + ags_a(isolalth) * fsolalt
!        aas = aas_a(isolaltl) * (1.-fsolalt) + aas_a(isolalth) * fsolalt
         if(isolaltl .ge. isolalt_lo .and. isolaltl .le. isolalt_hi)then
           ags = logint(ags_a(isolaltl),ags_a(isolalth),fsolalt)
           aas = logint(aas_a(isolaltl),aas_a(isolalth),fsolalt)
         else
           write(6,*)' ERROR: with isolalt',isolaltl,isolalt_lo,dsolalt,del_solalt,dsolaltb,float(isolalt_lo)/del_solalt
           stop
         endif

!        Considering smoothing out with solar disk
         horz_dep = horz_depf(htbar_msl,earth_radius)
         if(solalt_step + 0.5 .lt. -horz_dep)then
            sol_occ = 0.0 ! invisible
         else
            sol_occ = 1.0 ! visible
         endif

         if(solalt_step .ge. 0.)then
           zapp = 90. - solalt_step
           od_solar_slant_o3 = od_o_msl * patm_o3(htbar_msl) * airmasso(zapp,htbar_msl)
         elseif(sol_occ .gt. 0.)then
           zappi = 90. + solalt_step
           patm = exp(-(htbar_msl/scale_ht_g))
           ztruei = zappi + refractd_app(abs(solalt_step) ,patm)
           patmo  = patm_o3(htbar_msl)
           htmin_ray = htminf(htbar_msl,solalt_step,earth_radius)
           patm2o = patm_o3(htmin_ray)
           ztrue0 = 90. + refractd_app(0.,patm2o)
!          alttrue0 = -refractd_app(0.,patm2o)
!          ao1 = airmasso(ztruei,htmin_ray) * patm2o ! ztruei from htmin
           ao2 = airmasso(ztruei,htbar_msl) * patmo  ! upward from observer
           ao3 = airmasso(ztrue0,htmin_ray) * patm2o ! 0 deg from htmin
           ao_calc = (ao3-ao2) + ao3
           od_solar_slant_o3 = od_o_msl * ao_calc * sol_occ
         else
           od_solar_slant_o3 = 999.
         endif

         if(sol_occ .gt. 0.)then

           if(htmsl .le. 150e3 .and. .false.)then
!                        M         L              S
             od_solar_slant = &
                    od_g_vert * ags * exp(-htbar/scale_ht_g)  &
                  + od_solar_slant_o3 &
                  + od_a_vert * aas * exp(-htbar/scale_ht_a) * frac_iso_a
           else ! presently testing
             od_solar_slant_g = od_g_msl * ags * exp(-htbar_msl/scale_ht_g)
             od_solar_slant_a = aod_ref  * aas * exp(-(htbar_msl-redp_lvl)/scale_ht_a) * frac_iso_a
             od_solar_slant = od_solar_slant_g &
                            + od_solar_slant_o3 + od_solar_slant_a
           endif

           rad = trans(dble(od_solar_slant))
           radg = max(rad,.0001) ! secondary scattering
           rada = max(rad,.00001) ! secondary scattering
           rad = max(rad,.0001) ! secondary scattering

!          di = (dtau * rad) * exp(-tausum)
           if(dtau .gt. 0.)then
             ssa_eff = (alphabar_g + ssa * alphabar_a) / (alphabar_g + alphabar_a)
           else
             ssa_eff = 1.0
           endif

           alphabar = alphabar_g + alphabar_o + alphabar_a
           if(alphabar .gt. 0.)then
             di_g = do * radg       * (alphabar_g / alphabar)
             di_a = do * rada * ssa * (alphabar_a / alphabar)
           else
             di_g = 0.
             di_a = 0.
           endif
           di = di_g + di_a

         else ! sol_occ = 0
           od_solar_slant = 999.
           rad = 0. ; di_g = 0. ; di_a = 0. ; di = 0.

         endif ! sol_occ

         sumi_g = sumi_g + di_g
         sumi_a = sumi_a + di_a
         sumi = sumi + di

         frac_opac = opac_curr / opac_slant

         sumi_mean = sumi / opac(tausum)
!        sumi_extrap = (sumi + (1.-frac_opac) * rad) / opac_slant
         sumi_extrap = (sumi_mean * frac_opac) + (1.-frac_opac) * rad

         if(idebug .eq. 1)then
!          if(i .le. 10 .or. ((i .eq. (i/25)*25) .and. i .le. 400) .or. i .eq. 500 .or. i .eq. 1000 .or. i .eq. 1500 .or. i .eq. 2000. .or. i .eq. 2500 .or. i .eq. 10000)then
           if((i .le. 10 .or. i .eq. (i/25)*25) .and. htbar_msl .gt. 0. .and. htbar_msl .le. 110e3)then
             if(sol_occ.gt.0.)then
               write(6,12)horz_dep,solalt_step,htmin_ray,sol_occ ! ,ao2,ao3,ao
             endif
12           format(33x,'   horz_dep/solalt/htmin/sol_occ = ',2f8.2,f9.0,f8.3,3f9.4)
             write(6,11)i,sbar,htbar_msl,alphabar_g*ds,alphabar_a*ds,tausum,dsolalt,ags,aas,od_solar_slant_g,od_solar_slant,rad,di,sumi_g,sumi_a,opac_curr,frac_opac,sumi_mean,sumi_extrap
11           format(i6,f11.0,f12.1,2e10.3,f8.4,f9.4,f10.2,4f9.4,e9.2,2f11.8,2f9.4,2f9.4)
           endif
         endif
     enddo ! i

900  continue

     if(idebug .eq. 1)then
       write(6,11)i,sbar,htbar_msl,alphabar_g*ds,alphabar_a*ds,tausum,dsolalt,ags,aas,od_solar_slant_g,od_solar_slant,rad,di,sumi_g,sumi_a,opac_curr,frac_opac,sumi_mean,sumi_extrap
       if(sol_occ.gt.0.)write(6,12)horz_dep,solalt_step,htmin_ray,sol_occ ! ,ao2,ao3,ao
     endif

     srcdir = sumi_extrap

     return
     end

     subroutine get_clr_src_dir_topo(solalt,solazi,viewalt,viewazi &   ! I
           ,emis_ang,r_missing_data &                                  ! I
           ,od_g_msl,od_g_vert,od_o_msl,od_a_vert,htmsl,dist_to_topo & ! I
           ,ssa,agv,aav,aod_ref,redp_lvl,scale_ht_a &                  ! I
           ,ags_a,aas_a,isolalt_lo,isolalt_hi,del_solalt &             ! I
           ,ic,idebug,nsteps,refdist_solalt,solalt_ref &               ! I
           ,sumi_g,sumi_a,opac_slant,tausum_a)                         ! O

     use mem_namelist, ONLY: earth_radius
     use constants_laps, ONLY: R, GRAV
     include 'trigd.inc'
     include 'rad.inc'

     trans(od) = exp(-min(od,80.))
     opac(od) = 1.0 - trans(od)
     angdif(X,Y)=MOD(X-Y+540.,360.)-180.
     curvat(hdst,radius) = hdst**2 / (2. * radius)
     curvat2(sdst,radius_start,altray) = &
             sqrt(sdst**2 + radius_start**2 &
               - (2.*sdst*radius_start*cosd(90.+altray))) - radius_start
     real*8 dsdst,dradius_start,daltray,dcurvat2
     real*8 send,sbar,htbar,dist_to_topo,ds
     dcurvat2(dsdst,dradius_start,daltray) = &
              sqrt(dsdst**2 + dradius_start**2 &
       - (2D0*dsdst*dradius_start*cosd(90D0+daltray))) - dradius_start
     expbl(x) = exp(max(x,-80.))

     real tausum_a(nsteps)
     real ags_a(isolalt_lo:isolalt_hi) ! gas rel to zen val fr obs ht/salt
     real aas_a(isolalt_lo:isolalt_hi) ! aero rel to zen val fr obs ht/salt

!    Calculate src term of gas+aerosol along ray normalized by TOA radiance
     od_g_slant = od_g_vert * agv
     od_a_slant = od_a_vert * aav
     od_slant = od_g_slant + od_a_slant
     r_e = 6371.e3
     patm_htmsl = ztopsa(htmsl) / 1013.25
     radius_earth_eff = earth_radius * (1. + (0.33333 * patm_htmsl))
     radius_start_eff = radius_earth_eff + htmsl

     scale_ht_g = R * ztotsa(htmsl) / GRAV
!    alpha_sfc_g = od_g_vert / scale_ht_g

!    scale_ht_g = 8000.
     scale_ht_g = R * ztotsa(0.) / GRAV
     alpha_ref_g = od_g_msl  / scale_ht_g

!    alpha_sfc_a = od_a_vert / scale_ht_a
     alpha_ref_a = aod_ref / scale_ht_a

     opac_vert = opac(od_g_vert + od_a_vert)
     if(opac_vert .eq. 0)then ! very small od
         opac_vert = od_g_vert + od_a_vert
     endif
     opac_slant = opac(od_g_slant + od_a_slant)
     if(opac_slant .eq. 0)then ! very small od
         opac_slant = od_g_slant + od_a_slant
     endif

!    Integral [0 to x] a * exp(-b*x) = a * (1. - exp*(-b*x)) / b
     ds = min(dist_to_topo/10d0,1000d0)
     if(dist_to_topo .le. 10000d0)then
       nsteps_topo = nint(dist_to_topo/ds) ! max(10,int(htmsl/1000.))
     else
       nsteps_topo = int(dist_to_topo/ds) + 1
     endif
     tausum = 0.
     sumi = 0.; sumi_g = 0.; sumi_a = 0.
     httopill = max(100000.-sind(solalt)*60000.,40000.)
     htbotill = -(htmsl + 500.)
     httopo = 0.
     istart = max(1,nint((htmsl-100000d0)/ds))

     dsolalt_ref = solalt_ref - solalt

     if(idebug .eq. 1)then
         write(6,*)
         write(6,*)'subroutine get_clr_src_dir_topo'
         write(6,*)'solalt/ref/solazi = ',solalt,solalt_ref,solazi
         write(6,*)'viewalt/viewazi/emis = ',viewalt,viewazi,emis_ang
         write(6,*)'agv/aav = ',agv,aav
         write(6,*)'od_g_vert = ',od_g_vert
         write(6,*)'scale_ht_g = ',scale_ht_g
         write(6,*)'alpha_sfc_g/alpha_ref_g = ',alpha_sfc_g,alpha_ref_g
         write(6,*)'od_o_msl = ',od_o_msl
         write(6,*)'od_a_vert = ',od_a_vert
         write(6,*)'opac_vert = ',opac_vert
         write(6,*)'aod_ref/redp_lvl/scale_ht_a',aod_ref,redp_lvl,scale_ht_a
         write(6,*)'od_g_slant = ',od_g_slant
         write(6,*)'od_a_slant = ',od_a_slant
         write(6,*)'opac_slant = ',opac_slant
         write(6,*)'od_slant = ',od_slant
         write(6,*)'dist_to_topo = ',dist_to_topo
         write(6,*)'istart/nsteps_topo/ds = ',istart,nsteps_topo,ds
         write(6,*)'htbotill/httopill',htbotill,httopill
     endif ! i

     if(idebug .eq. 1)then ! substitute radg/rada?
       write(6,*)'         sbar      htbar      htmsl     dtau    tausum  dsolalt alphbr_g alphbr_a  od_solar   rad       di     sumi_g   sumi_a  opac_curr frac_opac sumi_mn sumi_ext'
     endif

     opac_curr = 0.

     dsolalt_dxy = cosd(angdif(solazi,viewazi)) / 110000.

     iwrite_slant = 0

     do i = istart,nsteps_topo
       sbar = (dble(i)-0.5d0) * ds
       send =  dble(i)        * ds
!      htbar = sbar * sind(viewalt)
!      xybar = sbar * cosd(viewalt)
!      htbar = htbar + curvat(xybar,earth_radius*(4./3.)) ! Earth Curvature
!      htbar = curvat2(sbar,radius_start_eff,viewalt) ! From GEO
       htbar = dcurvat2(sbar,dble(radius_start_eff),dble(viewalt)) !Fm Moon
       htbar_msl_last = htbar_msl
       htbar_msl = htbar+htmsl
       ht_msl_last = ht_msl
       ht_msl = 1.5 * htbar_msl - 0.5 * htbar_msl_last
       if(htbar_msl .lt. httopill)then
         if(send .lt. dist_to_topo)then
           frac_step = 1.0
         else
           frac_step = 1d0 - (send - dist_to_topo)/ds
         endif
         if(emis_ang .ne. r_missing_data)then
           xybar = (sbar - refdist_solalt) * cosd(emis_ang)
         else
           xybar = (sbar - refdist_solalt) * cosd(viewalt)
         endif

!        Update ags,aas
         dsolalt = dsolalt_ref + dsolalt_dxy * xybar
         dsolalt = min(max(dsolalt,-90.),+90.)
         solalt_step = solalt + dsolalt
         dsolaltb = &
            max(min(dsolalt/del_solalt,float(isolalt_hi)/del_solalt-.0001),float(isolalt_lo)/del_solalt)
         isolaltl = floor(dsolaltb); isolalth = isolaltl + 1
         fsolalt = dsolaltb - float(isolaltl)
         ags = ags_a(isolaltl) * (1.-fsolalt) + ags_a(isolalth) * fsolalt
         aas = aas_a(isolaltl) * (1.-fsolalt) + aas_a(isolalth) * fsolalt

!        Considering smoothing out with solar disk
         horz_dep = horz_depf(htbar_msl,earth_radius)
         if(solalt_step + 0.5 .lt. -horz_dep)then
            sol_occ = 0.0 ! invisible
         else
            sol_occ = 1.0 ! visible
         endif

         if(solalt_step .ge. 0.)then
           zapp = 90. - solalt_step
           od_solar_slant_o3 = od_o_msl * patm_o3(htbar_msl) * airmasso(zapp,htbar_msl)
         elseif(sol_occ .gt. 0.)then
           zappi = 90. + solalt_step
           patm = exp(-(htbar_msl/scale_ht_g))
           ztruei = zappi + refractd_app(abs(solalt_step) ,patm)
           patmo  = patm_o3(htbar_msl)
           htmin_ray = htminf(htbar_msl,solalt_step,earth_radius)
           patm2o = patm_o3(htmin_ray)
           ztrue0 = 90. + refractd_app(0.,patm2o)
           ao2 = airmasso(ztruei,htbar_msl) * patmo  ! upward from observer
           ao3 = airmasso(ztrue0,htmin_ray) * patm2o ! 0 deg from htmin
           ao_calc = (ao3-ao2) + ao3
           od_solar_slant_o3 = od_o_msl * ao_calc * sol_occ
         else ! more accurate calculation technique was imported
           od_solar_slant_o3 = 0.
         endif

         if(i .eq. nsteps_topo .and. idebug .eq. 1)then
             write(6,*)'aas terms',isolaltl,aas_a(isolaltl),isolalth,aas_a(isolalth),fsolalt
             write(6,*)'od_solar_slant terms',od_g_vert,ags,htbar,scale_ht_g,exp(-htbar/scale_ht_g),aod_ref,aas,htbar,scale_ht_a,expbl(-(htbar_msl-redp_lvl)/scale_ht_a)
             iwrite_slant = 1
         endif

         if(sol_occ .gt. 0.)then
           od_solar_slant = od_g_msl * ags * expbl(-htbar_msl/scale_ht_g) &
                    + od_solar_slant_o3 &
                    + aod_ref * aas * expbl(-(htbar_msl-redp_lvl)/scale_ht_a)
         else
           od_solar_slant = 999.
         endif

!        alphabar_g = alpha_sfc_g * exp(-htbar/scale_ht_g)
!        alphabar_g = alpha_ref_g * exp(-htbar_msl/scale_ht_g)
         patm_bar = ztopsa(htbar_msl) / 1013.25
         scale_ht_g = R * ztotsa(htbar_msl) / GRAV
         alphabar_g = (od_g_msl * patm_bar) / scale_ht_g

!        alphabar_a = alpha_sfc_a * exp(-htbar/scale_ht_a)
         alphabar_a = alpha_ref_a * expbl(-(htbar_msl-redp_lvl)/scale_ht_a)
         dtau = ds * (alphabar_g + alphabar_a) * frac_step
         tausum = tausum + dtau
!        if(tausum .lt. 1.0)distod1 = sbar
!        tausum_a(i) = tausum

         opac_last = opac_curr
         opac_curr = opac(tausum)
         do = opac_curr - opac_last

!        rad = 1.0
         rad = trans(od_solar_slant)
         radg = rad ! max(rad,.0001) ! secondary scattering
         rada = rad ! max(rad,.00001) ! secondary scattering
!        rad = max(rad,.0001) ! secondary scattering

!        di = (dtau * rad) * exp(-tausum)
         if(dtau .gt. 0.)then
             ssa_eff = (alphabar_g + ssa * alphabar_a) / (alphabar_g + alphabar_a)
         else
             ssa_eff = 1.0
         endif

         di   = do * rad  * ssa_eff
         di_g = do * radg       * alphabar_g / (alphabar_g + alphabar_a)
         di_a = do * rada * ssa * alphabar_a / (alphabar_g + alphabar_a)

         sumi   = sumi   + di
         sumi_g = sumi_g + di_g
         sumi_a = sumi_a + di_a

         frac_opac = opac_curr / opac_slant

         sumi_mean = sumi / opac(tausum)
!        sumi_extrap = (sumi + (1.-frac_opac) * rad) / opac_slant
         sumi_extrap = (sumi_mean * frac_opac) + (1.-frac_opac) * rad

         if(idebug .eq. 1)then
!          if(i .le. 10)then
             write(6,11)i,sbar,htbar,htbar_msl,dtau,tausum,dsolalt,alphabar_g,alphabar_a,od_solar_slant,rad,di,sumi_g,sumi_a,opac_curr,frac_opac,sumi_mean,sumi_extrap
11           format(i6,f11.0,f11.0,f9.0,f9.4,f9.6,f9.4,2f9.6,9f9.4)
!          endif
         endif
       endif ! htbar < httopill

       if(htbar_msl .lt. httopo .or. (htbar_msl .gt. httopill .and. viewalt .gt. 0.))then
         if(idebug .eq. 1)then
           write(6,*)' exit loop:',htbar_msl,ht_msl,httopo,frac_step
         endif
         goto 900
       endif
     enddo ! i

     if(idebug .eq. 1)then
       write(6,*)' end of loop:',htbar_msl,ht_msl,httopo,frac_step
     endif

900  continue

     return
     end

max_stars_repo_name: maxinye/laps-mirror

max_stars_repo_licenses: Intel Unlicense OLDAP-2.2.1 NetCDF

max_issues_repo_name: longwosion/laps-mirror

max_issues_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense

max_forks_repo_name: longwosion/laps-mirror

max_forks_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense


text:
/*****************************************************************************
 * wb_dma_test_base.f
 *****************************************************************************/
+tool.infact

+UVM_TESTNAME=wb_dma_transfer_test
+uvm_set_type_override=wb_dma_transfer_seq,wb_dma_infact_single_transfer_seq
+uvm_set_type_override=wb_dma_config_seq,wb_dma_infact_pairwise_config_seq

ext: f


text:
      program MainAngle
!     ------Programm for calculation main angle and ellipticity-------
      intejer j
      real ff,q,d1,d2,ksi,delta,u,z,k1,k2,k3
      complex e1,e2,e3
      parameter (pi=3.1415926)
      common ff,e1,e2,e3,q,d1,d2
      write(*,*) 'Input E1, E2, E3 /*.,-*/'
      read(*,*)e1,e2,e3
      write(*,*) 'Input length of wave'
      read(*,*)q
      open(1,file='EXIT.dat')
      write(*,*) 'Input thicknessis d1 and d2'
      read(*,*)d1,d2
 1    write(*,*) 'Input maximal angle of incidence'
      read(*,*)W
      write(*,*) 'Input minimal angle of incidence'
      read(*,*)V
      write(*,6)
      write(1,6)
      ff=W
      call FZU(z,u,delta,ksi)
      k1=u
      ff=V
      call FZU(z,u,delta,ksi)
      k2=u
      g=k1*k2
      if(f.lt.0) goto 2
      write(*,*) 'Root does not exist'
      goto 1
 2    ff=(W+V)/2
      call FZU(z,u,delta,ksi)
      k3=u
      if(abs(k3).le.0.0001) goto 4
      g=k1*k3
      if(g.lt.0) goto 3
      W=ff
      goto 2
 3    V=ff
      goto 2
 4    write(*,5) ff, u, z, delta, ksi
      write(1,5) ff, u, z, delta, ksi
      write(*,*) 'Continue? j=1, j=0 - Stop'
      read(*,*)j
      if(j.eq.1)goto 1
      if(j.eq.0)goto 7
 5    format(2x,f7.4,4(2x,f12.7))
 6    format(6x,'F',14x,'COS',14x,'TAN',14x,'DELTA',12x,'KSI')
 7    stop
      end program MainAngle
!     --------------Calculation of ellipsometric parameters-------------
      SUBROUTINE FZU(z,u,delta,ksi)
      real ff,q,d1,d2,ksi,delta,u,z
      complex dl1,dl2,e1,e2,e3,cf1,cf2,cf3
      complex r01p,r12p,r23p,r01s,r12s,r23s
      complex r,rs,rp,ex1,ex2,ex3
      parameter (pi=3.1415926)
      common ff,e1,e2,e3,q,d1,d2
      f=ff*pi/180
      cf1=CSQRT(1-sin(f)**2/(e1*e1))
      cf2=CSQRT(1-sin(f)**2/(e2*e2))
      cf3=CSQRT(1-sin(f)**2/(e3*e3))
      r01p=(e1*cos(f)-cf1)/(e1*cos(f)+cf1)
      r12p=(e2*cf1-e1*cf2)/(e2*cf1+e1*cf2)
      r23p=(e3*cf2-e2*cf3)/(e3*cf2+e2*cf3)
      r01s=(cos(f)-e1*cf1)/(cos(f)+e1*cf1)
      r12s=(e1*cf1-e2*cf2)/(e1*cf1+e2*cf2)
      r23s=(e2*cf1-e2*cf2)/(e1*cf1+e2*cf2)
      dl1=2*pi/q*d1*CSQRT(e1*e1-sin(f)**2)
      ex1=CEXP(2*(0.,-1.)*dl1)
      dl2=2*pi/q*d2*CSQRT(e2*e2-sin(f)**2)
      ex2=CEXP(2*(0.,-1.)*dl2)
      ex3=CEXP(2*(0.,-1.)*(dl1+dl2))
!     ---------------------------------------------------------------
      rp=r01p+r12p*ex1+r01p*r12p*r23p*ex2+r23p*ex3
      rp=rp/(1+r01p*r12p*ex1+r12p*r23p*ex2+r01p*r23p*ex3)
      rs=r01s+r12s*ex1+r01s*r12s*r23s*ex2+r23s*ex3
      rs=rs/(1+r01s+r12s*ex1+r12s*r23s*ex2+r01s*r23s*ex3)
      r=rp/rs
      z=CABS(r)
      u=REAL(r)/z
      ksi=ATAN(z)*180/pi
      delta=ACOS(u)*180/pi
      ff=f*180/pi
      end SUBROUTINE FZU

max_stars_repo_name: RomaniukVadim/ellipsometric-calculator

max_issues_repo_name: RomaniukVadim/ellipsometric-calculator

max_forks_repo_name: RomaniukVadim/ellipsometric-calculator


text:
*> \brief \b CGGHD3
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download CGGHD3 + dependencies
*>
*> [TGZ]
*>
*> [ZIP]
*>
*> [TXT]
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*        SUBROUTINE CGGHD3( COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q,
*       $                   LDQ, Z, LDZ, WORK, LWORK, INFO )
*
*        .. Scalar Arguments ..
*        CHARACTER          COMPQ, COMPZ
*        INTEGER            IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, N, LWORK
*        ..
*        .. Array Arguments ..
*        COMPLEX            A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
*       $                   Z( LDZ, * ), WORK( * )
*        ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*>
*> CGGHD3 reduces a pair of complex matrices (A,B) to generalized upper
*> Hessenberg form using unitary transformations, where A is a
*> general matrix and B is upper triangular.  The form of the
*> generalized eigenvalue problem is
*>    A*x = lambda*B*x,
*> and B is typically made upper triangular by computing its QR
*> factorization and moving the unitary matrix Q to the left side
*> of the equation.
*>
*> This subroutine simultaneously reduces A to a Hessenberg matrix H:
*>    Q**H*A*Z = H
*> and transforms B to another upper triangular matrix T:
*>    Q**H*B*Z = T
*> in order to reduce the problem to its standard form
*>    H*y = lambda*T*y
*> where y = Z**H*x.
*>
*> The unitary matrices Q and Z are determined as products of Givens
*> rotations.  They may either be formed explicitly, or they may be
*> postmultiplied into input matrices Q1 and Z1, so that
*>
*>      Q1 * A * Z1**H = (Q1*Q) * H * (Z1*Z)**H
*>
*>      Q1 * B * Z1**H = (Q1*Q) * T * (Z1*Z)**H
*>
*> If Q1 is the unitary matrix from the QR factorization of B in the
*> original equation A*x = lambda*B*x, then CGGHD3 reduces the original
*> problem to generalized Hessenberg form.
*>
*> This is a blocked variant of CGGHRD, using matrix-matrix
*> multiplications for parts of the computation to enhance performance.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] COMPQ
*> \verbatim
*>          COMPQ is CHARACTER*1
*>          = 'N': do not compute Q;
*>          = 'I': Q is initialized to the unit matrix, and the
*>                 unitary matrix Q is returned;
*>          = 'V': Q must contain a unitary matrix Q1 on entry,
*>                 and the product Q1*Q is returned.
*> \endverbatim
*>
*> \param[in] COMPZ
*> \verbatim
*>          COMPZ is CHARACTER*1
*>          = 'N': do not compute Z;
*>          = 'I': Z is initialized to the unit matrix, and the
*>                 unitary matrix Z is returned;
*>          = 'V': Z must contain a unitary matrix Z1 on entry,
*>                 and the product Z1*Z is returned.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The order of the matrices A and B.  N >= 0.
*> \endverbatim
*>
*> \param[in] ILO
*> \verbatim
*>          ILO is INTEGER
*> \endverbatim
*>
*> \param[in] IHI
*> \verbatim
*>          IHI is INTEGER
*>
*>          ILO and IHI mark the rows and columns of A which are to be
*>          reduced.  It is assumed that A is already upper triangular
*>          in rows and columns 1:ILO-1 and IHI+1:N.  ILO and IHI are
*>          normally set by a previous call to CGGBAL; otherwise they
*>          should be set to 1 and N respectively.
*>          1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
*> \endverbatim
*>
*> \param[in,out] A
*> \verbatim
*>          A is COMPLEX array, dimension (LDA, N)
*>          On entry, the N-by-N general matrix to be reduced.
*>          On exit, the upper triangle and the first subdiagonal of A
*>          are overwritten with the upper Hessenberg matrix H, and the
*>          rest is set to zero.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER
*>          The leading dimension of the array A.  LDA >= max(1,N).
*> \endverbatim
*>
*> \param[in,out] B
*> \verbatim
*>          B is COMPLEX array, dimension (LDB, N)
*>          On entry, the N-by-N upper triangular matrix B.
*>          On exit, the upper triangular matrix T = Q**H B Z.  The
*>          elements below the diagonal are set to zero.
*> \endverbatim
*>
*> \param[in] LDB
*> \verbatim
*>          LDB is INTEGER
*>          The leading dimension of the array B.  LDB >= max(1,N).
*> \endverbatim
*>
*> \param[in,out] Q
*> \verbatim
*>          Q is COMPLEX array, dimension (LDQ, N)
*>          On entry, if COMPQ = 'V', the unitary matrix Q1, typically
*>          from the QR factorization of B.
*>          On exit, if COMPQ='I', the unitary matrix Q, and if
*>          COMPQ = 'V', the product Q1*Q.
*>          Not referenced if COMPQ='N'.
*> \endverbatim
*>
*> \param[in] LDQ
*> \verbatim
*>          LDQ is INTEGER
*>          The leading dimension of the array Q.
*>          LDQ >= N if COMPQ='V' or 'I'; LDQ >= 1 otherwise.
*> \endverbatim
*>
*> \param[in,out] Z
*> \verbatim
*>          Z is COMPLEX array, dimension (LDZ, N)
*>          On entry, if COMPZ = 'V', the unitary matrix Z1.
*>          On exit, if COMPZ='I', the unitary matrix Z, and if
*>          COMPZ = 'V', the product Z1*Z.
*>          Not referenced if COMPZ='N'.
*> \endverbatim
*>
*> \param[in] LDZ
*> \verbatim
*>          LDZ is INTEGER
*>          The leading dimension of the array Z.
*>          LDZ >= N if COMPZ='V' or 'I'; LDZ >= 1 otherwise.
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is COMPLEX array, dimension (LWORK)
*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
*> \endverbatim
*>
*> \param[in]  LWORK
*> \verbatim
*>          LWORK is INTEGER
*>          The length of the array WORK.  LWORK >= 1.
*>          For optimum performance LWORK >= 6*N*NB, where NB is the
*>          optimal blocksize.
*>
*>          If LWORK = -1, then a workspace query is assumed; the routine
*>          only calculates the optimal size of the WORK array, returns
*>          this value as the first entry of the WORK array, and no error
*>          message related to LWORK is issued by XERBLA.
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>          = 0:  successful exit.
*>          < 0:  if INFO = -i, the i-th argument had an illegal value.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date January 2015
*
*> \ingroup complexOTHERcomputational
*
*> \par Further Details:
*  =====================
*>
*> \verbatim
*>
*>  This routine reduces A to Hessenberg form and maintains B in
*>  using a blocked variant of Moler and Stewart's original algorithm,
*>  as described by Kagstrom, Kressner, Quintana-Orti, and Quintana-Orti
*>  (BIT 2008).
*> \endverbatim
*>
*  =====================================================================
      SUBROUTINE CGGHD3( COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q,
     $                   LDQ, Z, LDZ, WORK, LWORK, INFO )
*
*  -- LAPACK computational routine (version 3.6.1) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     January 2015
*
*
      IMPLICIT NONE
*
*     .. Scalar Arguments ..
      CHARACTER          COMPQ, COMPZ
      INTEGER            IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, N, LWORK
*     ..
*     .. Array Arguments ..
      COMPLEX            A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
     $                   Z( LDZ, * ), WORK( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      COMPLEX            CONE, CZERO
      PARAMETER          ( CONE = ( 1.0E+0, 0.0E+0 ),
     $                     CZERO = ( 0.0E+0, 0.0E+0 ) )
*     ..
*     .. Local Scalars ..
      LOGICAL            BLK22, INITQ, INITZ, LQUERY, WANTQ, WANTZ
      CHARACTER*1        COMPQ2, COMPZ2
      INTEGER            COLA, I, IERR, J, J0, JCOL, JJ, JROW, K,
     $                   KACC22, LEN, LWKOPT, N2NB, NB, NBLST, NBMIN,
     $                   NH, NNB, NX, PPW, PPWO, PW, TOP, TOPQ
      REAL               C
      COMPLEX            C1, C2, CTEMP, S, S1, S2, TEMP, TEMP1, TEMP2,
     $                   TEMP3
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            ILAENV
      EXTERNAL           ILAENV, LSAME
*     ..
*     .. External Subroutines ..
      EXTERNAL           CGGHRD, CLARTG, CLASET, CUNM22, CROT, XERBLA
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          REAL, CMPLX, CONJG, MAX
*     ..
*     .. Executable Statements ..
*
*     Decode and test the input parameters.
*
      INFO = 0
      NB = ILAENV( 1, 'CGGHD3', ' ', N, ILO, IHI, -1 )
      LWKOPT = MAX( 6*N*NB, 1 )
      WORK( 1 ) = CMPLX( LWKOPT )
      INITQ = LSAME( COMPQ, 'I' )
      WANTQ = INITQ .OR. LSAME( COMPQ, 'V' )
      INITZ = LSAME( COMPZ, 'I' )
      WANTZ = INITZ .OR. LSAME( COMPZ, 'V' )
      LQUERY = ( LWORK.EQ.-1 )
*
      IF( .NOT.LSAME( COMPQ, 'N' ) .AND. .NOT.WANTQ ) THEN
         INFO = -1
      ELSE IF( .NOT.LSAME( COMPZ, 'N' ) .AND. .NOT.WANTZ ) THEN
         INFO = -2
      ELSE IF( N.LT.0 ) THEN
         INFO = -3
      ELSE IF( ILO.LT.1 ) THEN
         INFO = -4
      ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN
         INFO = -5
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -7
      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
         INFO = -9
      ELSE IF( ( WANTQ .AND. LDQ.LT.N ) .OR. LDQ.LT.1 ) THEN
         INFO = -11
      ELSE IF( ( WANTZ .AND. LDZ.LT.N ) .OR. LDZ.LT.1 ) THEN
         INFO = -13
      ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN
         INFO = -15
      END IF
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'CGGHD3', -INFO )
         RETURN
      ELSE IF( LQUERY ) THEN
         RETURN
      END IF
*
*     Initialize Q and Z if desired.
*
      IF( INITQ )
     $   CALL CLASET( 'All', N, N, CZERO, CONE, Q, LDQ )
      IF( INITZ )
     $   CALL CLASET( 'All', N, N, CZERO, CONE, Z, LDZ )
*
*     Zero out lower triangle of B.
*
      IF( N.GT.1 )
     $   CALL CLASET( 'Lower', N-1, N-1, CZERO, CZERO, B(2, 1), LDB )
*
*     Quick return if possible
*
      NH = IHI - ILO + 1
      IF( NH.LE.1 ) THEN
         WORK( 1 ) = CONE
         RETURN
      END IF
*
*     Determine the blocksize.
*
      NBMIN = ILAENV( 2, 'CGGHD3', ' ', N, ILO, IHI, -1 )
      IF( NB.GT.1 .AND. NB.LT.NH ) THEN
*
*        Determine when to use unblocked instead of blocked code.
*
         NX = MAX( NB, ILAENV( 3, 'CGGHD3', ' ', N, ILO, IHI, -1 ) )
         IF( NX.LT.NH ) THEN
*
*           Determine if workspace is large enough for blocked code.
*
            IF( LWORK.LT.LWKOPT ) THEN
*
*              Not enough workspace to use optimal NB:  determine the
*              minimum value of NB, and reduce NB or force use of
*              unblocked code.
*
               NBMIN = MAX( 2, ILAENV( 2, 'CGGHD3', ' ', N, ILO, IHI,
     $                      -1 ) )
               IF( LWORK.GE.6*N*NBMIN ) THEN
                  NB = LWORK / ( 6*N )
               ELSE
                  NB = 1
               END IF
            END IF
         END IF
      END IF
*
      IF( NB.LT.NBMIN .OR. NB.GE.NH ) THEN
*
*        Use unblocked code below
*
         JCOL = ILO
*
      ELSE
*
*        Use blocked code
*
         KACC22 = ILAENV( 16, 'CGGHD3', ' ', N, ILO, IHI, -1 )
         BLK22 = KACC22.EQ.2
         DO JCOL = ILO, IHI-2, NB
            NNB = MIN( NB, IHI-JCOL-1 )
*
*           Initialize small unitary factors that will hold the
*           accumulated Givens rotations in workspace.
*           N2NB   denotes the number of 2*NNB-by-2*NNB factors
*           NBLST  denotes the (possibly smaller) order of the last
*                  factor.
*
            N2NB = ( IHI-JCOL-1 ) / NNB - 1
            NBLST = IHI - JCOL - N2NB*NNB
            CALL CLASET( 'All', NBLST, NBLST, CZERO, CONE, WORK, NBLST )
            PW = NBLST * NBLST + 1
            DO I = 1, N2NB
               CALL CLASET( 'All', 2*NNB, 2*NNB, CZERO, CONE,
     $                      WORK( PW ), 2*NNB )
               PW = PW + 4*NNB*NNB
            END DO
*
*           Reduce columns JCOL:JCOL+NNB-1 of A to Hessenberg form.
*
            DO J = JCOL, JCOL+NNB-1
*
*              Reduce Jth column of A. Store cosines and sines in Jth
*              column of A and B, respectively.
*
               DO I = IHI, J+2, -1
                  TEMP = A( I-1, J )
                  CALL CLARTG( TEMP, A( I, J ), C, S, A( I-1, J ) )
                  A( I, J ) = CMPLX( C )
                  B( I, J ) = S
               END DO
*
*              Accumulate Givens rotations into workspace array.
*
               PPW  = ( NBLST + 1 )*( NBLST - 2 ) - J + JCOL + 1
               LEN  = 2 + J - JCOL
               JROW = J + N2NB*NNB + 2
               DO I = IHI, JROW, -1
                  CTEMP = A( I, J )
                  S = B( I, J )
                  DO JJ = PPW, PPW+LEN-1
                     TEMP = WORK( JJ + NBLST )
                     WORK( JJ + NBLST ) = CTEMP*TEMP - S*WORK( JJ )
                     WORK( JJ ) = CONJG( S )*TEMP + CTEMP*WORK( JJ )
                  END DO
                  LEN = LEN + 1
                  PPW = PPW - NBLST - 1
               END DO
*
               PPWO = NBLST*NBLST + ( NNB+J-JCOL-1 )*2*NNB + NNB
               J0 = JROW - NNB
               DO JROW = J0, J+2, -NNB
                  PPW = PPWO
                  LEN  = 2 + J - JCOL
                  DO I = JROW+NNB-1, JROW, -1
                     CTEMP = A( I, J )
                     S = B( I, J )
                     DO JJ = PPW, PPW+LEN-1
                        TEMP = WORK( JJ + 2*NNB )
                        WORK( JJ + 2*NNB ) = CTEMP*TEMP - S*WORK( JJ )
                        WORK( JJ ) = CONJG( S )*TEMP + CTEMP*WORK( JJ )
                     END DO
                     LEN = LEN + 1
                     PPW = PPW - 2*NNB - 1
                  END DO
                  PPWO = PPWO + 4*NNB*NNB
               END DO
*
*              TOP denotes the number of top rows in A and B that will
*              not be updated during the next steps.
*
               IF( JCOL.LE.2 ) THEN
                  TOP = 0
               ELSE
                  TOP = JCOL
               END IF
*
*              Propagate transformations through B and replace stored
*              left sines/cosines by right sines/cosines.
*
               DO JJ = N, J+1, -1
*
*                 Update JJth column of B.
*
                  DO I = MIN( JJ+1, IHI ), J+2, -1
                     CTEMP = A( I, J )
                     S = B( I, J )
                     TEMP = B( I, JJ )
                     B( I, JJ ) = CTEMP*TEMP - CONJG( S )*B( I-1, JJ )
                     B( I-1, JJ ) = S*TEMP + CTEMP*B( I-1, JJ )
                  END DO
*
*                 Annihilate B( JJ+1, JJ ).
*
                  IF( JJ.LT.IHI ) THEN
                     TEMP = B( JJ+1, JJ+1 )
                     CALL CLARTG( TEMP, B( JJ+1, JJ ), C, S,
     $                            B( JJ+1, JJ+1 ) )
                     B( JJ+1, JJ ) = CZERO
                     CALL CROT( JJ-TOP, B( TOP+1, JJ+1 ), 1,
     $                          B( TOP+1, JJ ), 1, C, S )
                     A( JJ+1, J ) = CMPLX( C )
                     B( JJ+1, J ) = -CONJG( S )
                  END IF
               END DO
*
*              Update A by transformations from right.
*
               JJ = MOD( IHI-J-1, 3 )
               DO I = IHI-J-3, JJ+1, -3
                  CTEMP = A( J+1+I, J )
                  S = -B( J+1+I, J )
                  C1 = A( J+2+I, J )
                  S1 = -B( J+2+I, J )
                  C2 = A( J+3+I, J )
                  S2 = -B( J+3+I, J )
*
                  DO K = TOP+1, IHI
                     TEMP = A( K, J+I  )
                     TEMP1 = A( K, J+I+1 )
                     TEMP2 = A( K, J+I+2 )
                     TEMP3 = A( K, J+I+3 )
                     A( K, J+I+3 ) = C2*TEMP3 + CONJG( S2 )*TEMP2
                     TEMP2 = -S2*TEMP3 + C2*TEMP2
                     A( K, J+I+2 ) = C1*TEMP2 + CONJG( S1 )*TEMP1
                     TEMP1 = -S1*TEMP2 + C1*TEMP1
                     A( K, J+I+1 ) = CTEMP*TEMP1 + CONJG( S )*TEMP
                     A( K, J+I ) = -S*TEMP1 + CTEMP*TEMP
                  END DO
               END DO
*
               IF( JJ.GT.0 ) THEN
                  DO I = JJ, 1, -1
                     C = DBLE( A( J+1+I, J ) )
                     CALL CROT( IHI-TOP, A( TOP+1, J+I+1 ), 1,
     $                          A( TOP+1, J+I ), 1, C,
     $                          -CONJG( B( J+1+I, J ) ) )
                  END DO
               END IF
*
*              Update (J+1)th column of A by transformations from left.
*
               IF ( J .LT. JCOL + NNB - 1 ) THEN
                  LEN  = 1 + J - JCOL
*
*                 Multiply with the trailing accumulated unitary
*                 matrix, which takes the form
*
*                        [  U11  U12  ]
*                    U = [            ],
*                        [  U21  U22  ]
*
*                 where U21 is a LEN-by-LEN matrix and U12 is lower
*                 triangular.
*
                  JROW = IHI - NBLST + 1
                  CALL CGEMV( 'Conjugate', NBLST, LEN, CONE, WORK,
     $                        NBLST, A( JROW, J+1 ), 1, CZERO,
     $                        WORK( PW ), 1 )
                  PPW = PW + LEN
                  DO I = JROW, JROW+NBLST-LEN-1
                     WORK( PPW ) = A( I, J+1 )
                     PPW = PPW + 1
                  END DO
                  CALL CTRMV( 'Lower', 'Conjugate', 'Non-unit',
     $                        NBLST-LEN, WORK( LEN*NBLST + 1 ), NBLST,
     $                        WORK( PW+LEN ), 1 )
                  CALL CGEMV( 'Conjugate', LEN, NBLST-LEN, CONE,
     $                        WORK( (LEN+1)*NBLST - LEN + 1 ), NBLST,
     $                        A( JROW+NBLST-LEN, J+1 ), 1, CONE,
     $                        WORK( PW+LEN ), 1 )
                  PPW = PW
                  DO I = JROW, JROW+NBLST-1
                     A( I, J+1 ) = WORK( PPW )
                     PPW = PPW + 1
                  END DO
*
*                 Multiply with the other accumulated unitary
*                 matrices, which take the form
*
*                        [  U11  U12   0  ]
*                        [                ]
*                    U = [  U21  U22   0  ],
*                        [                ]
*                        [   0    0    I  ]
*
*                 where I denotes the (NNB-LEN)-by-(NNB-LEN) identity
*                 matrix, U21 is a LEN-by-LEN upper triangular matrix
*                 and U12 is an NNB-by-NNB lower triangular matrix.
*
                  PPWO = 1 + NBLST*NBLST
                  J0 = JROW - NNB
                  DO JROW = J0, JCOL+1, -NNB
                     PPW = PW + LEN
                     DO I = JROW, JROW+NNB-1
                        WORK( PPW ) = A( I, J+1 )
                        PPW = PPW + 1
                     END DO
                     PPW = PW
                     DO I = JROW+NNB, JROW+NNB+LEN-1
                        WORK( PPW ) = A( I, J+1 )
                        PPW = PPW + 1
                     END DO
                     CALL CTRMV( 'Upper', 'Conjugate', 'Non-unit', LEN,
     $                           WORK( PPWO + NNB ), 2*NNB, WORK( PW ),
     $                           1 )
                     CALL CTRMV( 'Lower', 'Conjugate', 'Non-unit', NNB,
     $                           WORK( PPWO + 2*LEN*NNB ),
     $                           2*NNB, WORK( PW + LEN ), 1 )
                     CALL CGEMV( 'Conjugate', NNB, LEN, CONE,
     $                           WORK( PPWO ), 2*NNB, A( JROW, J+1 ), 1,
     $                           CONE, WORK( PW ), 1 )
                     CALL CGEMV( 'Conjugate', LEN, NNB, CONE,
     $                           WORK( PPWO + 2*LEN*NNB + NNB ), 2*NNB,
     $                           A( JROW+NNB, J+1 ), 1, CONE,
     $                           WORK( PW+LEN ), 1 )
                     PPW = PW
                     DO I = JROW, JROW+LEN+NNB-1
                        A( I, J+1 ) = WORK( PPW )
                        PPW = PPW + 1
                     END DO
                     PPWO = PPWO + 4*NNB*NNB
                  END DO
               END IF
            END DO
*
*           Apply accumulated unitary matrices to A.
*
            COLA = N - JCOL - NNB + 1
            J = IHI - NBLST + 1
            CALL CGEMM( 'Conjugate', 'No Transpose', NBLST,
     $                  COLA, NBLST, CONE, WORK, NBLST,
     $                  A( J, JCOL+NNB ), LDA, CZERO, WORK( PW ),
     $                  NBLST )
            CALL CLACPY( 'All', NBLST, COLA, WORK( PW ), NBLST,
     $                   A( J, JCOL+NNB ), LDA )
            PPWO = NBLST*NBLST + 1
            J0 = J - NNB
            DO J = J0, JCOL+1, -NNB
               IF ( BLK22 ) THEN
*
*                 Exploit the structure of
*
*                        [  U11  U12  ]
*                    U = [            ]
*                        [  U21  U22  ],
*
*                 where all blocks are NNB-by-NNB, U21 is upper
*                 triangular and U12 is lower triangular.
*
                  CALL CUNM22( 'Left', 'Conjugate', 2*NNB, COLA, NNB,
     $                         NNB, WORK( PPWO ), 2*NNB,
     $                         A( J, JCOL+NNB ), LDA, WORK( PW ),
     $                         LWORK-PW+1, IERR )
               ELSE
*
*                 Ignore the structure of U.
*
                  CALL CGEMM( 'Conjugate', 'No Transpose', 2*NNB,
     $                        COLA, 2*NNB, CONE, WORK( PPWO ), 2*NNB,
     $                        A( J, JCOL+NNB ), LDA, CZERO, WORK( PW ),
     $                        2*NNB )
                  CALL CLACPY( 'All', 2*NNB, COLA, WORK( PW ), 2*NNB,
     $                         A( J, JCOL+NNB ), LDA )
               END IF
               PPWO = PPWO + 4*NNB*NNB
            END DO
*
*           Apply accumulated unitary matrices to Q.
*
            IF( WANTQ ) THEN
               J = IHI - NBLST + 1
               IF ( INITQ ) THEN
                  TOPQ = MAX( 2, J - JCOL + 1 )
                  NH  = IHI - TOPQ + 1
               ELSE
                  TOPQ = 1
                  NH = N
               END IF
               CALL CGEMM( 'No Transpose', 'No Transpose', NH,
     $                     NBLST, NBLST, CONE, Q( TOPQ, J ), LDQ,
     $                     WORK, NBLST, CZERO, WORK( PW ), NH )
               CALL CLACPY( 'All', NH, NBLST, WORK( PW ), NH,
     $                      Q( TOPQ, J ), LDQ )
               PPWO = NBLST*NBLST + 1
               J0 = J - NNB
               DO J = J0, JCOL+1, -NNB
                  IF ( INITQ ) THEN
                     TOPQ = MAX( 2, J - JCOL + 1 )
                     NH  = IHI - TOPQ + 1
                  END IF
                  IF ( BLK22 ) THEN
*
*                    Exploit the structure of U.
*
                     CALL CUNM22( 'Right', 'No Transpose', NH, 2*NNB,
     $                            NNB, NNB, WORK( PPWO ), 2*NNB,
     $                            Q( TOPQ, J ), LDQ, WORK( PW ),
     $                            LWORK-PW+1, IERR )
                  ELSE
*
*                    Ignore the structure of U.
*
                     CALL CGEMM( 'No Transpose', 'No Transpose', NH,
     $                           2*NNB, 2*NNB, CONE, Q( TOPQ, J ), LDQ,
     $                           WORK( PPWO ), 2*NNB, CZERO, WORK( PW ),
     $                           NH )
                     CALL CLACPY( 'All', NH, 2*NNB, WORK( PW ), NH,
     $                            Q( TOPQ, J ), LDQ )
                  END IF
                  PPWO = PPWO + 4*NNB*NNB
               END DO
            END IF
*
*           Accumulate right Givens rotations if required.
*
            IF ( WANTZ .OR. TOP.GT.0 ) THEN
*
*              Initialize small unitary factors that will hold the
*              accumulated Givens rotations in workspace.
*
               CALL CLASET( 'All', NBLST, NBLST, CZERO, CONE, WORK,
     $                      NBLST )
               PW = NBLST * NBLST + 1
               DO I = 1, N2NB
                  CALL CLASET( 'All', 2*NNB, 2*NNB, CZERO, CONE,
     $                         WORK( PW ), 2*NNB )
                  PW = PW + 4*NNB*NNB
               END DO
*
*              Accumulate Givens rotations into workspace array.
*
               DO J = JCOL, JCOL+NNB-1
                  PPW  = ( NBLST + 1 )*( NBLST - 2 ) - J + JCOL + 1
                  LEN  = 2 + J - JCOL
                  JROW = J + N2NB*NNB + 2
                  DO I = IHI, JROW, -1
                     CTEMP = A( I, J )
                     A( I, J ) = CZERO
                     S = B( I, J )
                     B( I, J ) = CZERO
                     DO JJ = PPW, PPW+LEN-1
                        TEMP = WORK( JJ + NBLST )
                        WORK( JJ + NBLST ) = CTEMP*TEMP -
     $                                       CONJG( S )*WORK( JJ )
                        WORK( JJ ) = S*TEMP + CTEMP*WORK( JJ )
                     END DO
                     LEN = LEN + 1
                     PPW = PPW - NBLST - 1
                  END DO
*
                  PPWO = NBLST*NBLST + ( NNB+J-JCOL-1 )*2*NNB + NNB
                  J0 = JROW - NNB
                  DO JROW = J0, J+2, -NNB
                     PPW = PPWO
                     LEN  = 2 + J - JCOL
                     DO I = JROW+NNB-1, JROW, -1
                        CTEMP = A( I, J )
                        A( I, J ) = CZERO
                        S = B( I, J )
                        B( I, J ) = CZERO
                        DO JJ = PPW, PPW+LEN-1
                           TEMP = WORK( JJ + 2*NNB )
                           WORK( JJ + 2*NNB ) = CTEMP*TEMP -
     $                                          CONJG( S )*WORK( JJ )
                           WORK( JJ ) = S*TEMP + CTEMP*WORK( JJ )
                        END DO
                        LEN = LEN + 1
                        PPW = PPW - 2*NNB - 1
                     END DO
                     PPWO = PPWO + 4*NNB*NNB
                  END DO
               END DO
            ELSE
*
               CALL CLASET( 'Lower', IHI - JCOL - 1, NNB, CZERO, CZERO,
     $                      A( JCOL + 2, JCOL ), LDA )
               CALL CLASET( 'Lower', IHI - JCOL - 1, NNB, CZERO, CZERO,
     $                      B( JCOL + 2, JCOL ), LDB )
            END IF
*
*           Apply accumulated unitary matrices to A and B.
*
            IF ( TOP.GT.0 ) THEN
               J = IHI - NBLST + 1
               CALL CGEMM( 'No Transpose', 'No Transpose', TOP,
     $                     NBLST, NBLST, CONE, A( 1, J ), LDA,
     $                     WORK, NBLST, CZERO, WORK( PW ), TOP )
               CALL CLACPY( 'All', TOP, NBLST, WORK( PW ), TOP,
     $                      A( 1, J ), LDA )
               PPWO = NBLST*NBLST + 1
               J0 = J - NNB
               DO J = J0, JCOL+1, -NNB
                  IF ( BLK22 ) THEN
*
*                    Exploit the structure of U.
*
                     CALL CUNM22( 'Right', 'No Transpose', TOP, 2*NNB,
     $                            NNB, NNB, WORK( PPWO ), 2*NNB,
     $                            A( 1, J ), LDA, WORK( PW ),
     $                            LWORK-PW+1, IERR )
                  ELSE
*
*                    Ignore the structure of U.
*
                     CALL CGEMM( 'No Transpose', 'No Transpose', TOP,
     $                           2*NNB, 2*NNB, CONE, A( 1, J ), LDA,
     $                           WORK( PPWO ), 2*NNB, CZERO,
     $                           WORK( PW ), TOP )
                     CALL CLACPY( 'All', TOP, 2*NNB, WORK( PW ), TOP,
     $                            A( 1, J ), LDA )
                  END IF
                  PPWO = PPWO + 4*NNB*NNB
               END DO
*
               J = IHI - NBLST + 1
               CALL CGEMM( 'No Transpose', 'No Transpose', TOP,
     $                     NBLST, NBLST, CONE, B( 1, J ), LDB,
     $                     WORK, NBLST, CZERO, WORK( PW ), TOP )
               CALL CLACPY( 'All', TOP, NBLST, WORK( PW ), TOP,
     $                      B( 1, J ), LDB )
               PPWO = NBLST*NBLST + 1
               J0 = J - NNB
               DO J = J0, JCOL+1, -NNB
                  IF ( BLK22 ) THEN
*
*                    Exploit the structure of U.
*
                     CALL CUNM22( 'Right', 'No Transpose', TOP, 2*NNB,
     $                            NNB, NNB, WORK( PPWO ), 2*NNB,
     $                            B( 1, J ), LDB, WORK( PW ),
     $                            LWORK-PW+1, IERR )
                  ELSE
*
*                    Ignore the structure of U.
*
                     CALL CGEMM( 'No Transpose', 'No Transpose', TOP,
     $                           2*NNB, 2*NNB, CONE, B( 1, J ), LDB,
     $                           WORK( PPWO ), 2*NNB, CZERO,
     $                           WORK( PW ), TOP )
                     CALL CLACPY( 'All', TOP, 2*NNB, WORK( PW ), TOP,
     $                            B( 1, J ), LDB )
                  END IF
                  PPWO = PPWO + 4*NNB*NNB
               END DO
            END IF
*
*           Apply accumulated unitary matrices to Z.
*
            IF( WANTZ ) THEN
               J = IHI - NBLST + 1
               IF ( INITQ ) THEN
                  TOPQ = MAX( 2, J - JCOL + 1 )
                  NH  = IHI - TOPQ + 1
               ELSE
                  TOPQ = 1
                  NH = N
               END IF
               CALL CGEMM( 'No Transpose', 'No Transpose', NH,
     $                     NBLST, NBLST, CONE, Z( TOPQ, J ), LDZ,
     $                     WORK, NBLST, CZERO, WORK( PW ), NH )
               CALL CLACPY( 'All', NH, NBLST, WORK( PW ), NH,
     $                      Z( TOPQ, J ), LDZ )
               PPWO = NBLST*NBLST + 1
               J0 = J - NNB
               DO J = J0, JCOL+1, -NNB
                     IF ( INITQ ) THEN
                     TOPQ = MAX( 2, J - JCOL + 1 )
                     NH  = IHI - TOPQ + 1
                  END IF
                  IF ( BLK22 ) THEN
*
*                    Exploit the structure of U.
*
                     CALL CUNM22( 'Right', 'No Transpose', NH, 2*NNB,
     $                            NNB, NNB, WORK( PPWO ), 2*NNB,
     $                            Z( TOPQ, J ), LDZ, WORK( PW ),
     $                            LWORK-PW+1, IERR )
                  ELSE
*
*                    Ignore the structure of U.
*
                     CALL CGEMM( 'No Transpose', 'No Transpose', NH,
     $                           2*NNB, 2*NNB, CONE, Z( TOPQ, J ), LDZ,
     $                           WORK( PPWO ), 2*NNB, CZERO, WORK( PW ),
     $                           NH )
                     CALL CLACPY( 'All', NH, 2*NNB, WORK( PW ), NH,
     $                            Z( TOPQ, J ), LDZ )
                  END IF
                  PPWO = PPWO + 4*NNB*NNB
               END DO
            END IF
         END DO
      END IF
*
*     Use unblocked code to reduce the rest of the matrix
*     Avoid re-initialization of modified Q and Z.
*
      COMPQ2 = COMPQ
      COMPZ2 = COMPZ
      IF ( JCOL.NE.ILO ) THEN
         IF ( WANTQ )
     $      COMPQ2 = 'V'
         IF ( WANTZ )
     $      COMPZ2 = 'V'
      END IF
*
      IF ( JCOL.LT.IHI )
     $   CALL CGGHRD( COMPQ2, COMPZ2, N, JCOL, IHI, A, LDA, B, LDB, Q,
     $                LDQ, Z, LDZ, IERR )
      WORK( 1 ) = CMPLX( LWKOPT )
*
      RETURN
*
*     End of CGGHD3
*
      END

max_stars_repo_name: lingochamp/OpenBLAS

max_issues_repo_name: lingochamp/OpenBLAS

max_forks_repo_name: lingochamp/OpenBLAS


text:
PROGRAM q2
    IMPLICIT NONE
    REAL:: a, b, AreaTR, AreaSR, Coeff
    INTEGER::  N      ! N is number of segment
    WRITE(*, fmt = '(/A)', ADVANCE = 'NO') "Enter the number of segment: "
    READ(*, *) N
    WRITE(*, fmt = '(/A)') "Enter the value of a & b"
    READ(*, *) a, b

    CALL TrapezoidalRule(a, b, AreaTR, N)
    WRITE(*, *) "Area using Trapezoidal Rule is: ", AreaTR
    Coeff = sqrt(1./AreaTR)
    WRITE(*, *) "Coefficient A from Trapezoidal is: ", Coeff

    CALL SimpsonsRule(a, b, AreaSR, N)
    WRITE(*, *) "Area using Simpsons Rule is: ", AreaSR
    Coeff = sqrt(1./AreaSR)
    WRITE(*, *) "Coefficient A from Simpsons is: ", Coeff

    CONTAINS
    SUBROUTINE TrapezoidalRule(a, b, Area, N)
        IMPLICIT NONE
        REAL:: a, b, Area, Xi, h
        INTEGER:: N, i
        ! N is number of segment
        IF(N .LT. 1) THEN
            Area = 0.0
            RETURN
        ENDIF
        h = (b-a)/N
        Area = (f(a) + f(b)) / 2.0
        DO i = 1, N-1
            Xi = a + i * h
            Area = Area + f(Xi)
        ENDDO
        Area = h * Area
    END SUBROUTINE TrapezoidalRule

    SUBROUTINE SimpsonsRule(a, b, Area, N)
        IMPLICIT NONE
        REAL:: a, b, Area, h
        INTEGER:: N, i
        ! N is number of segment
        IF(N .LT. 2) THEN
            Area = 0.0
            RETURN
        ENDIF
        h = (b-a)/N
        Area = f(a) + f(b)
        DO i = 1, N/2
            Area = Area + 4 * f(a + (2*i - 1) * h)
        ENDDO
        DO i = 1, (N+1)/2 - 1
            Area = Area + 2 * f(a + 2 * i * h)
        ENDDO
        Area = (h/3) * Area
    END SUBROUTINE SimpsonsRule

    REAL FUNCTION f(x)
        IMPLICIT NONE
        REAL:: x
        f = EXP(-x**2)**2
        return
    END FUNCTION f

END PROGRAM q2

max_stars_repo_path: Numerical Integration/q2.f90

max_issues_repo_path: Numerical Integration/q2.f90

max_forks_repo_path: Numerical Integration/q2.f90


text:
!=========================================================================================
! Purpose: Read precalculated binary background error corvariance file
!          in GSI and modify it as you want.
! Version: 0.1
! Author: Feng Zhu, feng.zhu@wisc.edu
! Copyright: BSD
!=========================================================================================

program bereader
    implicit none
    !-------------------------------------------------------------------------------------
    ! settings
    !-------------------------------------------------------------------------------------
    integer, parameter :: i_kind = selected_int_kind(8) ! 1=byte, 4=short, 8=long, 16=llong
    integer, parameter :: r_kind = selected_real_kind(6) ! 6=single, 15=double, 20=quad

    ! for modifying
    ! real(r_kind) :: coeff = 1
    real(r_kind) :: coeff = .5
    !-------------------------------------------------------------------------------------

    character(len=*), parameter :: file_dir = '/data/users/fzhu/Tools/comGSI_v3.3_large_nmsgmax/fix/Big_Endian/'
    character(len=*), parameter :: file_name = 'nam_nmmstat_na.gcv'
    ! character(len=*), parameter :: file_name = 'nam_glb_berror.f77.gcv'
    character(len=*), parameter :: file_path = file_dir//file_name

    integer(i_kind) :: msig, mlat

    integer(i_kind) :: inerr = 22

    !-------------------------------------------------------------------------------------
    ! for reading 'berror_bal_reg'
    !-------------------------------------------------------------------------------------
    real(r_kind), dimension(:), allocatable ::  clat_avn, sigma_avn
    real(r_kind), dimension(:, :), allocatable :: wgv_avn ,bv_avn
    real(r_kind), dimension(:, :, :), allocatable :: agv_avn
    !-------------------------------------------------------------------------------------

    !-------------------------------------------------------------------------------------
    ! for reading 'berror_bal_reg'
    !-------------------------------------------------------------------------------------
    character(len=5) :: varshort
    character(len=256) :: var
    integer(i_kind) :: isig, istat
    real(r_kind),dimension(:, :), allocatable ::  corqq_avn
    real(r_kind), dimension(:, :), allocatable :: corz_avn, hwll_avn, vztdq_avn

    !-------------------------------------------------------------------------------------

    integer k

    call system('cp '//file_path//' ./input') ! backup the orginal file in the input dir
    call system('rm -f ./output/'//file_name) ! clean up the output dir

    !=====================================================================================
    ! Read 'berror_stats'.
    !=====================================================================================
    write(*, *)
    write(*, *) '-------------------'
    write(*, *) 'START: BE reader...'
    write(*, *) '-------------------'
    write(*, *)
    write(*, *) 'BE file path: ', file_path
    write(*, *)

    ! berror_get_dims OUT: msig, mlat
    open(99, file='./output/'//file_name, form='unformatted', status='replace', convert='big_endian')
    open(inerr, file='./input/'//file_name, form='unformatted', status='old', convert='big_endian')
    rewind inerr
    read(inerr) msig, mlat
    write(*, *) 'msig: ', msig
    write(*, *) 'mlat: ', mlat
    write(*, *)
    write(99) msig, mlat

    ! berror_read_bal_reg IN: msig, mlat; OUT: agvi, wgvi, bvi
    allocate(clat_avn(mlat))
    allocate(sigma_avn(1:msig))
    allocate(agv_avn(0:mlat+1, 1:msig, 1:msig))
    allocate(bv_avn(0:mlat+1, 1:msig), wgv_avn(0:mlat+1, 1:msig))

    read(inerr) clat_avn, (sigma_avn(k), k=1,msig)
    read(inerr) agv_avn, bv_avn, wgv_avn
    write(*, *) 'Reading bal done.'
    write(*, *)
    write(99) clat_avn, (sigma_avn(k), k=1,msig)
    write(99) agv_avn, bv_avn, wgv_avn

    ! berror_read_wgt_reg IN: msig, mlat; OUT: corz, corp, hwll, hwllp, vz, rlsig
    read: do
        read(inerr, iostat=istat) varshort, isig
        write(*, *) 'varshort: ', varshort
        write(*, *) 'isig: ', isig
        var = varshort
        if (istat /= 0) exit
        write(99, iostat=istat) varshort, isig
        allocate(corz_avn(1:mlat, 1:isig))
        allocate(hwll_avn(0:mlat+1, 1:isig))
        allocate(vztdq_avn(1:isig, 0:mlat+1))
        if (var/='q') then
            read(inerr) corz_avn
            write(99) corz_avn*coeff
        else
            allocate ( corqq_avn(1:mlat,1:isig) )
            read(inerr) corz_avn, corqq_avn
            write(99) corz_avn*coeff, corqq_avn*coeff
        end if

        read(inerr) hwll_avn
        write(99) hwll_avn

        if (isig>1) then
            read(inerr) vztdq_avn
            write(99) vztdq_avn
        end if

        ! write(*, *) 'corz_avn: ', corz_avn

        deallocate ( corz_avn )
        deallocate ( hwll_avn )
        deallocate ( vztdq_avn )
        if (var=='q') deallocate ( corqq_avn )
    enddo read
    close(inerr)
    close(99)

    deallocate(clat_avn,sigma_avn)
    deallocate(agv_avn,bv_avn,wgv_avn)

    write(*, *)
    write(*, *) 'Reading wgt done.'
    write(*, *)

    write(*, *) '-------------------------------'
    write(*, *) 'DONE: Reading BE file complete!'
    write(*, *) '-------------------------------'
    write(*, *)
    close(inerr)
    !=====================================================================================

end program bereader

max_stars_repo_name: lyricorpse/BeReader

max_issues_repo_name: lyricorpse/BeReader

max_forks_repo_name: lyricorpse/BeReader


text:
program test
  integer  :: a(10) ,I = 0,j

  a = (/ (I, I = 1, 10) /)
  print *,I
  print *,a
end program test

max_stars_repo_name: clementval/fc

max_issues_repo_name: clementval/fc

max_forks_repo_name: clementval/fc


text:
C Copyright(C) 1999-2020 National Technology & Engineering Solutions
C of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains certain rights in this software.
C
C See packages/seacas/LICENSE for details

C=======================================================================
      SUBROUTINE PTSNRM( PLAPTS, INPT,  CUTPT, CUTNRM , IERR)
C=======================================================================

C   --*** PTRNRM *** (MESH) Turns three point plane to point-normal form
C   --   Written by Ray J. Meyers 29 May, 1990
C   --
C   -- Given three non-colinear points, returns a description of the plane
C   -- consisting of a point on the plane and a normal vector to the plane.
C   --
C   --Parameters:
C   --   PLAPTS - IN - the 3 points defining the cutting plane;
C   --   INPT - IN - a point in the mesh, ie, on the opposite side
C                    from the cutting plane normal
C   --   CUTPT - OUT - a point on the cutting plane;
C   --   CUTNRM - OUT - the normal to the cutting plane;
C   --   IERR - OUT - non-zero if the input points are colinear (no plane
C                     defined)

      REAL PLAPTS(3,3)
      REAL INPT(3)
      REAL CUTPT(3)
      REAL CUTNRM(3)
      INTEGER IERR
      REAL VEC1(3), VEC2(3)
      REAL TOL, DOT

C DEFINE TWO VECTORS IN THE PLANE : VEC1 = P1-P2, VEC2 = P3-P2

      VEC1(1) = PLAPTS(1,1) - PLAPTS(1,2)
      VEC1(2) = PLAPTS(2,1) - PLAPTS(2,2)
      VEC1(3) = PLAPTS(3,1) - PLAPTS(3,2)
      VEC2(1) = PLAPTS(1,3) - PLAPTS(1,2)
      VEC2(2) = PLAPTS(2,3) - PLAPTS(2,2)
      VEC2(3) = PLAPTS(3,3) - PLAPTS(3,2)

C DEFINE A TOLERANCE BASED ON THE LENGTH OF  VEC1

      TOL = 1E-06 * AMAX1 ( VEC1(1), VEC1(2), VEC1(3) )

C TAKE THE CROSS PRODUCT OF VEC1 AND VEC2 AS THE NORMAL OF THE PLANE AND
C NORMALIZE

      CUTNRM(1) = VEC1(2)*VEC2(3) - VEC1(3)*VEC2(2)
      CUTNRM(2) = VEC1(3)*VEC2(1) - VEC1(1)*VEC2(3)
      CUTNRM(3) = VEC1(1)*VEC2(2) - VEC1(2)*VEC2(1)

      DIST = SQRT( CUTNRM(1)*CUTNRM(1) + CUTNRM(2)*CUTNRM(2) +
     $             CUTNRM(3)*CUTNRM(3) )
      CUTNRM(1) = CUTNRM(1)/DIST
      CUTNRM(2) = CUTNRM(2)/DIST
      CUTNRM(3) = CUTNRM(3)/DIST

C IF THE NORMAL IS (0,0,0), THE ORIGINAL POINTS WERE COLINEAR

      IF( ABS(CUTNRM(1)) .LT. TOL .AND. ABS(CUTNRM(2)) .LT. TOL
     $    .AND. ABS(CUTNRM(3)) .LT. TOL) THEN
          CALL PRTERR ('CMDERR', 'Points do not define a plane')
          IERR = 1
      ELSE
          IERR = 0

C USE THE SECOND INPUT POINT AS THE CHOSEN CUTPT

          CUTPT(1) = PLAPTS(1,2)
          CUTPT(2) = PLAPTS(2,2)
          CUTPT(3) = PLAPTS(3,2)

C CHECK THE INPT TO SEE IF THE NORMAL IS THE CORRECT SENSE

          VEC1(1) = INPT(1) - CUTPT(1)
          VEC1(2) = INPT(2) - CUTPT(2)
          VEC1(3) = INPT(3) - CUTPT(3)
          DOT = VEC1(1)*CUTNRM(1) + VEC1(2)*CUTNRM(2) +
     $          VEC1(3)*CUTNRM(3)
          IF(DOT .GT. 0) THEN
             CUTNRM(1) = -CUTNRM(1)
             CUTNRM(2) = -CUTNRM(2)
             CUTNRM(3) = -CUTNRM(3)
          END IF
      END IF

      RETURN
      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text: Ryan Gonzales,18, resides and works in Davis. Now done with high school he has started doing important things with his life. He works as a production specialist at Davis Community Television and has been working there since the end of 2005. He plays drums in a hardcore punk band currently called No Incentive. With some friends Ryan does a radio show called Better than Silence on KDRT 95.7 fm thursdays and fridays.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
SUBROUTINE CALFQC(ISTL_,IS2TL_,MVAR,MO,CON,CON1,IT)

  ! **  SUBROUTINE CALFQC CALCULATES MASS SOURCES AND SINKS ASSOCIATED
  ! **  WITH CONSTANT AND TIME SERIES INFLOWS AND OUTFLOWS; CONTROL
  ! **  STRUCTURE INFLOWS AND OUTLOWS; WITHDRAWAL AND RETURN STRUCTURE
  ! **  OUTFLOWS; AND  EMBEDED CHANNEL INFLOWS AND OUTFLOWS
  !
  ! CHANGE RECORD
  ! DATE MODIFIED     BY               DESCRIPTION
  !----------------------------------------------------------------------------------!
  !    2015-11       PAUL M. CRAIG     ELIMINATED DUPLICATE 3TL/2TL FQC COMPUTATIONS
  !                                    ADDED THE NEW HYDRAULIC CONTROL STRUCTURES WITH
  !                                    RETURN FLOWS OPTIONALLY PUT IN DIFFRENT LAYERS
  !    2015-06       PAUL M. CRAIG     IMPLEMENTED SIGMA-Z (SGZ) IN EE7.3
  !    2011-03       Paul M. Craig     Rewritten to F90
  !----------------------------------------------------------------------------------!

  USE GLOBAL
  !  USE OMP_LIB
  IMPLICIT NONE

  INTEGER, INTENT(IN) :: ISTL_, IS2TL_, MVAR, MO, IT

  REAL    :: CON(LCM,KCM), CON1(LCM,KCM)
  REAL    :: FQOUT(NQCTLM), QTOUT(NQCTLM)

  INTEGER :: M,L,K,ID,JD,KD,NWR,IU,JU,KU,LU,NS
  INTEGER :: LD,NMD,NJP,LJP,KTMP,NTMP
  INTEGER :: NQSTMP, NCSTMP, NCTL
  INTEGER :: LMDCHHT, LMDCHUT, LMDCHVT

  REAL    :: QWRABS,RPORTS,QVJPTMP,QCJPTMP,QVJPENT,CONUP
  REAL    :: RQWD, QUKTMP, QVKTMP, TMPVAL, DOSAT

  LOGICAL :: USENQSIJ

  ! *** VARIABLES USED FOR ANTI-DIFFUSION CALCULATIONS
  ! ***
  ! *** FQCPAD  - MASS LOADINGS INTO MODEL DOMAIN FROM ALL BC TYPES FOR EACH CELL AND LAYER  (OUTFLOWS EXCLUDED)
  ! *** QSUMPAD - TOTAL FLOWS INTO MODEL DOMAIN FROM ALL BC TYPES FOR EACH CELL AND LAYER  (OUTFLOWS EXCLUDED)

  M=MO

  ! *** SELECTIVE ZEROING
  IF( KC > 1 )THEN
    IF( NGWSER > 0 .OR. ISGWIT/=0 )THEN
      DO L=1,LC
        FQC(L,KSZ(L),IT)=0.
        CONQ(L,KSZ(L),IT)=0.
      ENDDO
    ENDIF

    ! *** ZERO EVAP/RAINFALL
    IF( MVAR == 2 )THEN
      DO L=1,LC
        FQC(L,KC,IT)=0.
        CONQ(L,KC,IT)=0.
      ENDDO
      IF( ISADAC(MVAR) >= 2 )THEN
        DO L=1,LC
          FQCPAD(L,KC,IT)=0.
        ENDDO
      ENDIF
      IF( ISADAC(MVAR) > 0 )THEN
        DO L=1,LC
          QSUMPAD(L,KC,IT)=0.
        ENDDO
      ENDIF
    ENDIF

    ! *** ZERO ALL DEFINED BCs
    DO NS=1,NBCS
      L=LBCS(NS)
      DO K=1,KC
        FQC(L,K,IT)=0.
        CONQ(L,K,IT)=0.
        FQCPAD(L,K,IT)=0
        QSUMPAD(L,K,IT)=0.
      ENDDO
    ENDDO

  ELSE
    ! *** KC=1 RESET ALL
    FQC(1:LCM,KC,IT)=0.
    CONQ(1:LCM,KC,IT)=0.
    IF( ISADAC(MVAR) >= 2 )THEN
      FQCPAD(1:LCM,1:KCM,IT)=0.
    ENDIF
    QSUMPAD(1:LCM,1:KCM,IT)=0.
  ENDIF

  ! **  CORRECT CONCENTRATION FOR WATER SURFACE ELEVATION DATA ASSIMILATION
  !IF(ISWSEDA > 0 )THEN
  !  DO K=1,KC
  !    DO L=1,LC
  !      FQC(L,K)=FQC(L,K)+QWSEDA(L,K)*CONQ(L,K)
  !      QSUMPAD(L,K)=QSUMPAD(L,K)+MAX(QWSEDA(L,K),0.)
  !     QSUMNAD(L,K)=QSUMNAD(L,K)+MIN(QWSEDA(L,K),0.)
  !    ENDDO
  !  ENDDO
  !ENDIF

  ! *** HANDLE LEGACY WATER QUALITY MASS LOADING BC TYPE, I.E. IWQPSL = 0 OR IWQPSL = 1
  USENQSIJ = .TRUE.
  IF( MVAR == 8 .AND. IWQPSL /= 2 ) USENQSIJ = .FALSE.

  ! **************************************************************************************
  ! ***     INITIALIZE VOLUMETRIC SOURCE-SINK FLUXES AND AUXILLARY VARIABLES
  ! **************************************************************************************

  ! *** CONQ INITIALIZATION: 3TL STANDARD TIME STEP
  IF( ISTL_ == 3 )THEN
    ! *** INITIALIZE BOTTOM LAYER FOR GW INTERACTIONS
    IF( NGWSER > 0 .OR. ISGWIT /= 0 )THEN
      DO L=1,LC
        CONQ(L,KSZ(L),IT) = CON(L,KSZ(L))
      ENDDO
    ENDIF

    ! *** INITIALIZE TOP LAYER FOR EVAP/RAINFALL
    IF( MVAR == 2 )THEN
      DO L=1,LC
        CONQ(L,KC,IT)=CON(L,KC)
      ENDDO
    ENDIF

    ! *** INITIALIZE ALL DEFINED BCs
    DO NS=1,NBCS
      L=LBCS(NS)
      DO K=KSZ(L),KC
        CONQ(L,K,IT)=CON(L,K)
      ENDDO
    ENDDO
  ENDIF  ! *** END OF CONQ INITIALIZATION: 3TL STANDARD TIME STEP

  ! *** CONQ INITIALIZATION: 3TL CORRECTION STEP
  IF( ISTL_ == 2 .AND. IS2TL_ == 0 )THEN
    ! *** INITIALIZE BOTTOM LAYER FOR GW INTERACTIONS
    IF( NGWSER > 0 .OR. ISGWIT /= 0 )THEN
      DO L=1,LC
        CONQ(L,KSZ(L),IT) = 0.5*(CON(L,KSZ(L))+CON1(L,KSZ(L)))
      ENDDO
    ENDIF

    ! *** INITIALIZE TOP LAYER FOR EVAP/RAINFALL
    IF( MVAR == 2 )THEN
      DO L=1,LC
        CONQ(L,KC,IT) = 0.5*(CON(L,KC)+CON1(L,KC))
      ENDDO
    ENDIF

    ! *** INITIALIZE ALL DEFINED BCs
    DO NS=1,NBCS
      L=LBCS(NS)
      DO K=KSZ(L),KC
        CONQ(L,K,IT) = 0.5*(CON(L,K)+CON1(L,K))
      ENDDO
    ENDDO
  ENDIF   ! *** END OF CONQ INITIALIZATION: 3TL CORRECTION STEP

  ! *** CONQ INITIALIZATION: 2TL STANDARD TIME STEP
  IF( IS2TL_ == 1 )THEN
    ! *** INITIALIZE BOTTOM LAYER FOR GW INTERACTIONS
    IF( NGWSER > 0 .OR. ISGWIT /= 0 )THEN
      DO L=1,LC
        CONQ(L,KSZ(L),IT) = 0.5*(3.*CON(L,KSZ(L))-CON1(L,KSZ(L)))
      ENDDO
    ENDIF

    ! *** INITIALIZE TOP LAYER FOR EVAP/RAINFALL
    IF( MVAR == 2 )THEN
      DO L=1,LC
        CONQ(L,KC,IT) = 0.5*(3.*CON(L,KC)-CON1(L,KC))
      ENDDO
    ENDIF

    ! *** INITIALIZE ALL DEFINED BCs
    DO NS=1,NBCS
      L=LBCS(NS)
      DO K=KSZ(L),KC
        CONQ(L,K,IT) = 0.5*(3.*CON(L,K)-CON1(L,K))
      ENDDO
    ENDDO
  ENDIF  ! *** END OF CONQ INITIALIZATION: 2TL STANDARD TIME STEP

  ! *** GET MASS OUT FOR CONTROL STRUCTURES
  IF( NQCTL > 0 )THEN
    DO NCTL=1,NQCTL
      IU=IQCTLU(NCTL)
      JU=JQCTLU(NCTL)
      LU=LIJ(IU,JU)
      FQOUT(NCTL) = 0.
      QTOUT(NCTL) = 0.
      DO K=KSZ(LU),KC
        QTOUT(NCTL) = QTOUT(NCTL) + QCTLT(K,NCTL,1)
        FQOUT(NCTL) = FQOUT(NCTL) + QCTLT(K,NCTL,1)*CONQ(LU,K,IT)
      ENDDO
    ENDDO
  ENDIF

  IF( MVAR == 8 )THEN
    M=4+NTOX+NSED+NSND+MO
  ENDIF

  ! *********************************************************************!
  ! *** STANDARD VOLUMETRICS SOURCE SINK LOCATIONS
  ! *** EFDCPlus - 2015-11 - ELIMINATED DUPLICATE MASS LOADINGS SECTIONS FOR THE DIFFERENT TIME STEP OPTIONS

  ! *** FLOW BOUNDARY CELLS
  IF( USENQSIJ )THEN
    DO NS=1,NQSIJ
      L=LQS(NS)
      NQSTMP=NQSERQ(NS)
      NCSTMP=NCSERQ(NS,M)
      DO K=KSZ(L),KC
        ! ***                             |----------------  IN  -----------------|   |-----------  OUT  ----------------|
        FQC(L,K,IT)     = FQC(L,K,IT)     + MAX(QSS(K,NS),     0.) *CQS(K,NS,M)       + MIN(QSS(K,NS),0.)*CONQ(L,K,IT)      &
                                          + MAX(QSERCELL(K,NS),0.) *CSERT(K,NCSTMP,M) + MIN(QSERCELL(K,NS),0.)*CONQ(L,K,IT)

        FQCPAD(L,K,IT)  = FQCPAD(L,K,IT)  + MAX(QSS(K,NS),0.)*CQS(K,NS,M) + MAX(QSERCELL(K,NS),0.)*CSERT(K,NCSTMP,M)
        QSUMPAD(L,K,IT) = QSUMPAD(L,K,IT) + MAX(QSS(K,NS),0.)             + MAX(QSERCELL(K,NS),0.)
        !QSUMNAD(L,K,IT)=QSUMNAD(L,K,IT)+MIN(QSS(K,NS),0.)+MIN(QSERCELL(K,NS),0.)
      ENDDO
    ENDDO
  ENDIF

  ! ***  JET-PLUME VOLUMETRICS SOURCE SINK LOCATIONS
  IF( NQJPIJ > 0 )THEN
    DO NJP=1,NQJPIJ
      IF( ICALJP(NJP) == 1 )THEN
        ! *** DISCHARGE ONLY USING QSER AND STANDARD CONCENTRATION SERIES
        RPORTS = FLOAT(NPORTJP(NJP))
        LJP = LIJ(IQJP(NJP),JQJP(NJP))
        KTMP = KEFFJP(NJP)
        ! ***  QVJPTMP=TIME SERIES DISCHARGE FROM JET-PLUME
        QVJPTMP = 0.
        DO K=KSZ(LJP),KC
          QVJPTMP = QVJPTMP+QSERT(K,NQSERJP(NJP))
        ENDDO

        ! **QCJPTMP=ENTRAINMENT FLUX
        QCJPTMP=0.
        QVJPENT=0.
        ! *** REMOVE ENTRAINMENT FLUX AND CALCULATE TOTAL ENTRAIMENT FLUX
        DO K=KSZ(LJP),KC
          FQC(LJP,K,IT)  = FQC(LJP,K,IT) - QJPENT(K,NJP)*CONQ(LJP,K,IT)*RPORTS
          QCJPTMP        = QCJPTMP       + QJPENT(K,NJP)*CONQ(LJP,K,IT)
          QVJPENT        = QVJPENT       + QJPENT(K,NJP)
          !QSUMNAD(LJP,K,IT)=QSUMNAD(LJP,K,IT)-RPORTS*QJPENT(K,NJP)
        ENDDO

        ! *** PLACE JET FLUX AND ENTRAINMENT FLUX IS EFFECTIVE LAYER
        FQC(LJP,KTMP,IT)     =     FQC(LJP,KTMP,IT) + RPORTS*QCJPTMP + RPORTS*QQCJP(NJP)*CQCJP(1,NJP,M) + RPORTS*QVJPTMP*CSERT(1,NCSERJP(NJP,M),M)
        FQCPAD(LJP,KTMP,IT)  =  FQCPAD(LJP,KTMP,IT) + RPORTS*QCJPTMP + RPORTS*QQCJP(NJP)*CQCJP(1,NJP,M) + RPORTS*QVJPTMP*CSERT(1,NCSERJP(NJP,M),M)
        QSUMPAD(LJP,KTMP,IT) = QSUMPAD(LJP,KTMP,IT) + RPORTS*QVJPENT + RPORTS*QQCJP(NJP)                + RPORTS*QVJPTMP
      ENDIF

      IF( ICALJP(NJP) == 2 )THEN
        ! *** WITHDRAWAL/RETURN TYPE USING W/R SERIES
        RPORTS=FLOAT(NPORTJP(NJP))
        LJP=LIJ(IQJP(NJP),JQJP(NJP))
        KTMP=KEFFJP(NJP)
        NS=NQWRSERJP(NJP)
        LU=LIJ(IUPCJP(NJP),JUPCJP(NJP))
        KU=KUPCJP(NJP)
        CONUP=CONQ(LU,KU,IT)
        QCJPTMP=0.
        QVJPENT=0.

        ! *** REMOVE ENTRAIMENT FLUX AND CALCULATE TOTAL ENTRAINMENT
        DO K=KSZ(LJP),KC
          FQC(LJP,K,IT)=FQC(LJP,K,IT)-RPORTS*QJPENT(K,NJP)*CONQ(LJP,K,IT)
          QCJPTMP=QCJPTMP+QJPENT(K,NJP)*CONQ(LJP,K,IT)
          QVJPENT=QVJPENT+QJPENT(K,NJP)
          !QSUMNAD(LJP,K,IT)=QSUMNAD(LJP,K,IT)-RPORTS*QJPENT(K,NJP)
        ENDDO

        ! *** PLACE ENTRAINMENT, CONSTANT AND TIME SERIES FLUXES IN EFFECTIVE CELL
        FQC(LJP,KTMP,IT)     =     FQC(LJP,KTMP,IT) + RPORTS*QCJPTMP+RPORTS*QWRCJP(NJP)*(CWRCJP(NJP,M)+CONUP)+RPORTS*QWRSERT(NS)*(CQWRSERT(NS,M)+CONUP)
        FQCPAD(LJP,KTMP,IT)  =  FQCPAD(LJP,KTMP,IT) + RPORTS*QCJPTMP+RPORTS*QWRCJP(NJP)*(CWRCJP(NJP,M)+CONUP)+RPORTS*QWRSERT(NS)*(CQWRSERT(NS,M)+CONUP)
        QSUMPAD(LJP,KTMP,IT) = QSUMPAD(LJP,KTMP,IT) + RPORTS*QVJPENT+RPORTS*QWRCJP(NJP)+RPORTS*QWRSERT(NS)

        ! *** REMOVAL WITHDRAWAL FROM UPSTREAM CELL
        FQC(LU,KU,IT)=FQC(LU,KU,IT)-RPORTS*QWRCJP(NJP)*CONUP-RPORTS*QWRSERT(NS)*CONUP
        !QSUMNAD(LU,KU,IT)=QSUMNAD(LU,KU,IT)-RPORTS*QWRCJP(NJP)-RPORTS*QWRSERT(NS)
      ENDIF
    ENDDO
  ENDIF

  ! *** CONTROL STRUCTURES
  DO NCTL=1,NQCTL
    RQWD=1.
    IU=IQCTLU(NCTL)
    JU=JQCTLU(NCTL)
    LU=LIJ(IU,JU)
    ID=IQCTLD(NCTL)
    JD=JQCTLD(NCTL)
    IF( ID == 0 .AND. JD == 0 )THEN
      LD=LC
      RQWD=0.
    ELSE
      LD=LIJ(ID,JD)
    ENDIF
    DO K=KSZ(LU),KC
      FQC(LU,K,IT) = FQC(LU,K,IT) - QCTLT(K,NCTL,1)*CONQ(LU,K,IT)
      IF( QTOUT(NCTL) > 1.E-6 )THEN
        ! *** INTO THE DOMAIN ONLY
        TMPVAL = RQWD*QCTLT(K,NCTL,2)/QTOUT(NCTL)*FQOUT(NCTL)
        FQC(LD,K,IT)     = FQC(LD,K,IT)    + TMPVAL
        FQCPAD(LD,K,IT)  = FQCPAD(LD,K,IT) + TMPVAL
        QSUMPAD(LD,K,IT) = QSUMPAD(LD,K,IT)+ RQWD*QCTLT(K,NCTL,2)
        !QSUMNAD(L,K,IT)=QSUMNAD(L,K,IT)-QCTLT(K,NCTL,1)
      ENDIF
    ENDDO
  ENDDO

  ! *** WITHDRAWAL CONCENTRATION RISE RETURN
  DO NWR=1,NQWR
    ! *** Handle +/- Flows for Withdrawal/Return Structures
    NQSTMP=NQWRSERQ(NWR)
    IF( QWRSERT(NQSTMP) >= 0. )THEN
      !! *** Original Withdrawal/Return
      IU=IQWRU(NWR)
      JU=JQWRU(NWR)
      KU=KQWRU(NWR)
      ID=IQWRD(NWR)
      JD=JQWRD(NWR)
      KD=KQWRD(NWR)
    ELSE
      ! *** Reverse Flow Withdrawal/Return
      ID=IQWRU(NWR)
      JD=JQWRU(NWR)
      KD=KQWRU(NWR)
      IU=IQWRD(NWR)
      JU=JQWRD(NWR)
      KU=KQWRD(NWR)
      QWR(NWR)=0.  ! *** Only allow time variable flows when using -W/R
    ENDIF
    QWRABS = ABS(QWRSERT(NQSTMP))
    LU=LIJ(IU,JU)
    LD=LIJ(ID,JD)
    NCSTMP=NQWRSERQ(NWR)

    FQC(LU,KU,IT)=   FQC(LU,KU,IT)   -(QWR(NWR)+QWRABS)*CONQ(LU,KU,IT)
    FQC(LD,KD,IT)=   FQC(LD,KD,IT)   +QWR(NWR)*(CONQ(LU,KU,IT)+CQWR(NWR,M))+QWRABS*(CONQ(LU,KU,IT)+CQWRSERT(NCSTMP,M))
    FQCPAD(LD,KD,IT)=FQCPAD(LD,KD,IT)+QWR(NWR)*(CONQ(LU,KU,IT)+CQWR(NWR,M))+QWRABS*(CONQ(LU,KU,IT)+CQWRSERT(NCSTMP,M))
    QSUMPAD(LD,KD,IT)=QSUMPAD(LD,KD,IT)+QWR(NWR)+QWRABS
    !QSUMNAD(LU,KU,IT)=QSUMNAD(LU,KU,IT)-(QWR(NWR)+QWRSERT(NQSTMP))
  ENDDO

  ! *** SUBGRID SCALE CHANNEL EXCHANGE
  IF( MDCHH >= 1 )THEN
    DO K=1,KC
      DO NMD=1,MDCHH
        LMDCHHT=LMDCHH(NMD)
        LMDCHUT=LMDCHU(NMD)
        LMDCHVT=LMDCHV(NMD)
        IF( MDCHTYP(NMD) == 1 )THEN
          QUKTMP=QCHANU(NMD)*DZC(L,K)
          QVKTMP=0.
        ENDIF
        IF( MDCHTYP(NMD) == 2 )THEN
          QVKTMP=QCHANV(NMD)*DZC(L,K)
          QUKTMP=0.
        ENDIF
        IF( MDCHTYP(NMD) == 3 )THEN
          QUKTMP=QCHANU(NMD)*DZC(L,K)
          QVKTMP=QCHANV(NMD)*DZC(L,K)
        ENDIF
        FQC(LMDCHHT,K,IT)=FQC(LMDCHHT,K,IT)+MAX(QUKTMP,0.)*CONQ(LMDCHUT,K,IT)+MIN(QUKTMP,0.)*CONQ(LMDCHHT,K,IT)  &
                                            +MAX(QVKTMP,0.)*CONQ(LMDCHVT,K,IT)+MIN(QVKTMP,0.)*CONQ(LMDCHHT,K,IT)
        FQC(LMDCHUT,K,IT)=FQC(LMDCHUT,K,IT)-MAX(QUKTMP,0.)*CONQ(LMDCHUT,K,IT)-MIN(QUKTMP,0.)*CONQ(LMDCHHT,K,IT)
        FQC(LMDCHVT,K,IT)=FQC(LMDCHVT,K,IT)-MAX(QVKTMP,0.)*CONQ(LMDCHVT,K,IT)-MIN(QVKTMP,0.)*CONQ(LMDCHHT,K,IT)
      ENDDO
    ENDDO
  ENDIF

  ! *** GROUNDWATER FLUXES
  IF( ISGWIT > 0. )THEN
    IF( ISTRAN(5) > 0 .AND. MVAR == 5 )THEN
      DO NC=1,NTMP
        DO L=2,LA
          CONGW(L,NC) = GWCSERT(NGWSL(L),NC)
        END DO
      END DO
    ENDIF

    ! *** COMPUTE IN/OUT FLUXES BUT BYPASS SEDIMENT LOADS
    IF( MVAR /= 6 .AND. MVAR /= 7 )THEN
      DO L=2,LA
        IF( QGW(L) < 0. )THEN        ! *** 2018-10-24 PMC CHANGED SIGN CONVENTION FOR SEEPAGE +(IN), -(OUT)
          IF( H1P(L) > HDRY )THEN
            FQC(L,KSZ(L),IT) = FQC(L,KSZ(L),IT) + QGW(L)*CONQ(L,KSZ(L),IT)
          ENDIF
        ELSE
          FQC(L,KSZ(L),IT) = FQC(L,KSZ(L),IT) + QGW(L)*GWCSERT(NGWSL(L),M)
        ENDIF
      ENDDO
    ENDIF
  ENDIF

  ! *** ICE MELT FOR WATER QUALITY CONSTITUENTS
  IF( MVAR == 8 .AND. ISICE > 2 )THEN
    IF( MO == 19 )THEN
      DOSAT = 14.   ! *** FRESH WATER D.O. SATURATION FOR NEAR FREEZING WATER (G/M3)
      DO L=2,LA
        IF( ICERATE(L) > 0. )THEN
          FQC(L,KC,IT)     = FQC(L,KC,IT)     + ICERATE(L)*DOSAT
          !FQCPAD(L,KC,IT)  = FQCPAD(L,KC,IT)  + ICERATE(L)*DOSAT
          !QSUMPAD(L,KC,IT) = QSUMPAD(L,KC,IT) + ICERATE(L)
        ENDIF
        ICERATE(L) = 0.
      ENDDO
    ENDIF
  ENDIF

  ! *** TEMPERATURE ADJUSTMENTS FOR RAINFALL & EVAPORATION
  IF( M == 2 )THEN
    IF( ISTOPT(2) == 0 .OR. ISTOPT(2) == 3 )THEN
      ! ***
      DO L=2,LA
        FQC(L,KC,IT)=FQC(L,KC,IT)+RAINT(L)*TEMO*DXYP(L)
      ENDDO
    ELSE    ! ISTOPT(2) == 1 .OR. ISTOPT(2) == 2 .OR. ISTOPT(2) == 4  .OR. ISTOPT(2) == 5
      DO L=2,LA
        FQC(L,KC,IT)     = FQC(L,KC,IT)     + RAINT(L)*TATMT(L)*DXYP(L)
        FQCPAD(L,KC,IT)  = FQCPAD(L,KC,IT)  + RAINT(L)*TATMT(L)*DXYP(L)
        QSUMPAD(L,KC,IT) = QSUMPAD(L,KC,IT) + RAINT(L)*DXYP(L)
      ENDDO
    ENDIF
  ENDIF

  IF( M == 2 )THEN
    IF( ISTOPT(2) == 0 )THEN
      DO L=2,LA
        FQC(L,KC,IT) = FQC(L,KC,IT) - EVAPSW(L)*CONQ(L,KC,IT)
      ENDDO
    ENDIF
  ENDIF

  RETURN

END

max_stars_repo_name: dsi-llc/EFDCPlus

max_forks_repo_name: dsi-llc/EFDCPlus


text:
      SUBROUTINE INDXA (I,IR,IDXA,NA)
      COMMON /CBLKT/  NPP        ,K          ,EPS        ,CNV        ,
     1                NM         ,NCMPLX     ,IK
      NA = 2**IR
      IDXA = I-NA+1
      IF (I-NM) 102,102,101
  101 NA = 0
  102 RETURN
      END

max_stars_repo_name: sbwhitecap/clocc-hg

max_issues_repo_name: sbwhitecap/clocc-hg

max_forks_repo_name: sbwhitecap/clocc-hg


text:
      SUBROUTINE free_mem_boozer
      USE booz_params
      USE booz_persistent
      IMPLICIT NONE
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER :: istat1=0, istat2=0, istat3=0
C-----------------------------------------------
      IF (ALLOCATED(jlist)) DEALLOCATE(jlist)
      IF (ALLOCATED(lsurf_boz)) DEALLOCATE(lsurf_boz)
      !IF (ALLOCATED(jlist)) DEALLOCATE (jlist, lsurf_boz)

      IF (ALLOCATED(bsubumnc)) DEALLOCATE(bsubumnc)
      IF (ALLOCATED(bsubvmnc)) DEALLOCATE(bsubvmnc)
      IF (ALLOCATED(bmodmnc)) DEALLOCATE(bmodmnc)
      IF (ALLOCATED(rmnc)) DEALLOCATE(rmnc)
      IF (ALLOCATED(zmns)) DEALLOCATE(zmns)
      IF (ALLOCATED(lmns)) DEALLOCATE(lmns)
      IF (ALLOCATED(xm)) DEALLOCATE(xm)
      IF (ALLOCATED(xn)) DEALLOCATE(xn)
      IF (ALLOCATED(xm_nyq)) DEALLOCATE(xm_nyq)
      IF (ALLOCATED(xn_nyq)) DEALLOCATE(xn_nyq)
      IF (ALLOCATED(hiota)) DEALLOCATE(hiota)
      IF (ALLOCATED(phip)) DEALLOCATE(phip)
      IF (ALLOCATED(gpsi)) DEALLOCATE(gpsi)
      IF (ALLOCATED(ipsi)) DEALLOCATE(ipsi)
      IF (ALLOCATED(pmns)) DEALLOCATE(pmns)
      IF (ALLOCATED(beta_vol)) DEALLOCATE(beta_vol)
      IF (ALLOCATED(pres)) DEALLOCATE(pres)
      IF (ALLOCATED(phi)) DEALLOCATE(phi)
      IF (ALLOCATED(buco)) DEALLOCATE(buco)
      IF (ALLOCATED(bvco)) DEALLOCATE(bvco)
      IF (ALLOCATED(rmncb)) DEALLOCATE(rmncb)
      IF (ALLOCATED(zmnsb)) DEALLOCATE(zmnsb)
      IF (ALLOCATED(pmnsb)) DEALLOCATE(pmnsb)
      IF (ALLOCATED(gmncb)) DEALLOCATE(gmncb)
      IF (ALLOCATED(bmncb)) DEALLOCATE(bmncb)
      IF (ALLOCATED(bmod_b)) DEALLOCATE(bmod_b)
!      IF (ALLOCATED(bsubumnc)) DEALLOCATE(
!     1    bsubumnc, bsubvmnc, bmodmnc, rmnc, zmns, lmns,
!     2    xm, xn, xm_nyq, xn_nyq, hiota, phip, gpsi,
!     3    ipsi, pmns, beta_vol, pres, phi, buco, bvco,
!     3    rmncb, zmnsb, pmnsb, gmncb, bmncb, bmod_b, stat=istat1 )

      IF (ALLOCATED(cosm_b)) DEALLOCATE(cosm_b)
      IF (ALLOCATED(sinm_b)) DEALLOCATE(sinm_b)
      IF (ALLOCATED(cosn_b)) DEALLOCATE(cosn_b)
      IF (ALLOCATED(sinn_b)) DEALLOCATE(sinn_b)
      IF (ALLOCATED(cosm_nyq)) DEALLOCATE(cosm_nyq)
      IF (ALLOCATED(sinm_nyq)) DEALLOCATE(sinm_nyq)
      IF (ALLOCATED(cosn_nyq)) DEALLOCATE(cosn_nyq)
      IF (ALLOCATED(sinn_nyq)) DEALLOCATE(sinn_nyq)
      IF (ALLOCATED(sfull)) DEALLOCATE(sfull)
      IF (ALLOCATED(scl)) DEALLOCATE(scl)
      IF (ALLOCATED(xmb)) DEALLOCATE(xmb)
      IF (ALLOCATED(xnb)) DEALLOCATE(xnb)
      IF (ALLOCATED(thgrd)) DEALLOCATE(thgrd)
      IF (ALLOCATED(ztgrd)) DEALLOCATE(ztgrd)

!      IF (ALLOCATED(cosm_b)) DEALLOCATE(
!     1    cosm_b, sinm_b, cosn_b, sinn_b,
!     2    cosm_nyq, sinm_nyq, cosn_nyq, sinn_nyq,
!     2    sfull, scl, xmb, xnb, thgrd, ztgrd, stat=istat2)

      IF (ALLOCATED(bsubumns)) DEALLOCATE(bsubumns)
      IF (ALLOCATED(bsubvmns)) DEALLOCATE(bsubvmns)
      IF (ALLOCATED(bmodmns)) DEALLOCATE(bmodmns)
      IF (ALLOCATED(rmns)) DEALLOCATE(rmns)
      IF (ALLOCATED(zmnc)) DEALLOCATE(zmnc)
      IF (ALLOCATED(lmnc)) DEALLOCATE(lmnc)
      IF (ALLOCATED(pmnc)) DEALLOCATE(pmnc)
      IF (ALLOCATED(rmnsb)) DEALLOCATE(rmnsb)
      IF (ALLOCATED(zmncb)) DEALLOCATE(zmncb)
      IF (ALLOCATED(pmncb)) DEALLOCATE(pmncb)
      IF (ALLOCATED(gmnsb)) DEALLOCATE(gmnsb)
      IF (ALLOCATED(bmnsb)) DEALLOCATE(bmnsb)
!      IF (ALLOCATED(bsubumns)) DEALLOCATE(
!     1    bsubumns, bsubvmns, bmodmns, rmns, zmnc, lmnc, pmnc,
!     2    rmnsb, zmncb, pmncb, gmnsb, bmnsb, stat=istat3)

      IF (istat1 .ne.0 .or. istat2 .ne. 0)
     1  PRINT *,' Deallocation error in Free_mem_boozer'

      END SUBROUTINE free_mem_boozer

max_stars_repo_name: joseluisvelasco/STELLOPT

max_issues_repo_name: joseluisvelasco/STELLOPT

max_forks_repo_name: joseluisvelasco/STELLOPT


text:
      MODULE m_felim
c.........................................................felim
c energy limits setup and checking
c
      CONTAINS
      SUBROUTINE felim(
     >                 mrad,lll,zz,nqn,vv,rc,
     <                 elim)
C     -----------------
C--->  FIND  E-LIMIT
C     -----------------

      IMPLICIT NONE
C     ..
C     .. Scalar Arguments ..
      INTEGER, INTENT (IN) :: mrad,lll,nqn
      REAL   , INTENT (IN) :: zz
      REAL   , INTENT (OUT):: elim
C     ..
C     .. Array Arguments ..
      REAL   , INTENT (IN) :: rc(mrad),vv(mrad)
C     ..
C     .. Local Scalars ..
      REAL val
      INTEGER n
C     ..
      IF (lll.EQ.0) THEN
         elim = -2.0*zz*zz/ (1.50*nqn*nqn)
      ELSE
         elim = vv(1) + lll/rc(1)**2
         DO 10 n = 2,mrad
            val = vv(n) + lll/rc(n)**2
            IF (val.LE.elim) elim = val
   10    CONTINUE
      END IF
      END SUBROUTINE felim
      END MODULE m_felim

max_stars_repo_name: MRedies/FLEUR

max_issues_repo_name: MRedies/FLEUR

max_forks_repo_name: MRedies/FLEUR