text:
 	SUBROUTINE RA_WIND  ( irpntr, sknt, drct, gust, ipwbef,
     +			      ipwaft, iret )
C************************************************************************
C* RA_WIND								*
C*									*
C* This subroutine finds the wind information in an airways report.	*
C* The first numeric field of length 4 is assumed to be the wind	*
C* field.  Pointers to the first field before and after the wind	*
C* field are returned.							*
C*									*
C* RA_WIND  ( IRPNTR, SKNT, DRCT, GUST, IPWBEF, IPWAFT, IRET )		*
C*									*
C* Input parameters:							*
C*	IRPNTR		INTEGER		First field after header	*
C*									*
C* Output parameters:							*
C*	SKNT		REAL		Wind speed in knots		*
C*	DRCT		REAL		Wind direction in degrees	*
C*	GUST		REAL		Wind gust as reported		*
C*	IPWBEF		INTEGER		Field before wind		*
C*	IPWAFT		INTEGER		Field after wind		*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C*					 -3 = no wind group found    	*
C**									*
C* Log:									*
C* B. Doty/RDS		 9/87	GEMPAK 4				*
C* I. Graffman/RDS	 3/88						*
C* M. desJardins/GSFC	 9/89	GEMPAK5					*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
	INCLUDE		'racmn.cmn'
C*
	LOGICAL		good
C----------------------------------------------------------------------
	iret   = 0
	ipwbef = 0
	ipwaft = 0
	sknt   = RMISSD
	drct   = RMISSD
	gust   = RMISSD
C
C*	Find first four digit number after the pointer field.
C
	good  = .false.
	iwind = irpntr
	DO WHILE  ( ( .not. good ) .and. ( iwind .le. nfield ) )
	    IF  ( ( iftype (iwind) .eq. 2 ) .and.
     +		  ( ifsize (iwind) .eq. 4 ) )  THEN
		good = .true.
	      ELSE
		iwind = iwind + 1
	    END IF
	END DO
C
C*	If there are more than 4 slashes already, this is not the wind.
C
	IF  ( good )  THEN
	    islash = 0
	    DO  i = irpntr, iwind
		IF  ( iftype (i) .eq. 3 )  islash = islash + 1
	    END DO
	    IF  ( islash .gt. 4 )  good = .false.
	END IF
C
C*	If found, get the wind.
C
	IF  ( good )  THEN
	    ival   = ifintg (iwind)
	    idrctn = ival / 100
	    ispeed = MOD ( ival, 100 )
C
C*	    Check numbers for validity.
C
	     IF  ( idrctn .gt. 36 )  THEN
		IF  (( idrctn .lt. 51 ) .or. ( idrctn .gt. 86 ))  THEN
		    good = .false.
		  ELSE
		    ispeed = ispeed + 100
		    idrctn = idrctn - 50
		END IF
	      ELSE
		IF  ( idrctn .lt. 0 )  good = .false.
		IF  ( idrctn .eq. 36 )  idrctn = 0
	    END IF
	    IF  ((ispeed .eq. 0) .and. (idrctn .ne. 0))  good = .false.
	END IF
C
C*	Check for estimated winds and look for gusts.
C
	IF  ( good )  THEN
	    ipwbef = iwind - 1
C
C*	    Move pointer for special case of estimated winds.
C
	    IF  ( ( iftype (ipwbef) .eq. 1 ) .and.
     +		  ( cfield (ipwbef) .eq. 'E' ) )  THEN
		ipwbef = ipwbef - 1
	    END IF
	    ipwaft = iwind + 1
	    sknt   = FLOAT (ispeed)
	    drct   = FLOAT (idrctn) * 10.
C
C*	    Gusts are reported as a G followed by a number.
C
	    IF  ( ipwaft .lt. nfield )  THEN
		IF  ( ( cfield (ipwaft) .eq. 'G' ) .or.
     +		      ( cfield (ipwaft) .eq. '+' ) )  THEN
		    ipwaft = ipwaft + 1
		    IF  ( iftype (ipwaft) .eq. 2 )  THEN
			gust   = ifintg (ipwaft)
			ipwaft = ipwaft + 1
		    END IF
		END IF
	    END IF
	END IF
C
C*	If wind was not found, check for pattern:  #/#//
C*	Assume that the two slashes indicate missing wind.
C
	IF  ( .not. good )  THEN
	    iwind = irpntr + 3
	    DO WHILE  ( ( .not. good ) .and. ( iwind .lt. nfield ) )
		IF  ( ( iftype (iwind)   .eq. 3 ) .and.
     +		      ( iftype (iwind+1) .eq. 3 ) .and.
     +		      ( iftype (iwind-1) .eq. 2 ) .and.
     +		      ( iftype (iwind-2) .eq. 3 ) .and.
     +		      ( iftype (iwind-3) .eq. 2 ) )  THEN
		    ipwbef = iwind
		    ipwaft = iwind + 1
		    good   = .true.
		  ELSE
		    iwind  = iwind + 1
		END IF
	    END DO
	END IF
C
C*	If not valid, set the return code.
C
	IF  ( .not. good )  iret = -3
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
The APS (Anaerobic Phased Solids) Digester is quite possibly the coolest thing ever.

Essentially, food waste, green waste and any other organic (as in containing fixed carbon) material goes in one end, and electricity (generated from methane produced by bacteria), compost and fiber products come out the other.

The technology is being developed by Dr. Ruihong Zhang a researcher for the department of Biological and Agricultural Engineering, and http://www.onsitepowersystems.com Onsite Power Systems, who have licensed the technology. Zhangs digester is better than previous designs because hers does not require the solids to be ground up (which costs energy) before being fed into the digester.

Onsite power has a good deal of information on their http://www.onsitepowersystems.com/DigesterOverview.htm website. You can also read the http://www.californiaaggie.com/news/2004/10/27/SciTech/Solid.Waste.Digester.Turns.Garbage.Into.Electricity1318557.shtml Aggie article but it has no paragraphs.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


!c  topocorrect - approximate topo correction
!c  Reference :
!http://earth-info.nga.mil/GandG/publications/tr8350.2/tr8350.2-a/Appendix.pdf

    subroutine topo(demAccessor, dopAccessor, slrngAccessor)
      use topozeroState
      use topozeroMethods
      use geometryModule
      use orbitModule
      use linalg3Module
      use fortranUtils, ONLY : getPI

      implicit none
      include 'omp_lib.h'
      integer*8 demAccessor, dopAccessor, slrngAccessor
      integer lineFile, stat
      !integer width, length
      real*8, allocatable ::lat(:),lon(:),z(:),zsch(:)
      real*4, allocatable :: losang(:), incang(:), elevang(:)
      real*4, allocatable :: distance(:)
      real*8, allocatable :: rho(:), dopline(:)
      real*4, allocatable :: dem(:,:), demline(:)

      integer*1, allocatable :: mask(:), omask(:)
      real*8, allocatable :: orng(:), ctrack(:)
      real*8, allocatable :: oview(:)
      real*8 ctrackmin, ctrackmax, dctrack
      real*8 sch(3),xyz(3),llh(3),delta(3)

      real*8 tline, rng, dopfact
      real*8 llh_prev(3), xyz_prev(3)
      real*8 xyzsat(3), velsat(3)
      real*8 schsat(3), llhsat(3)
      real*8 ltpsat(3), ltpvel(3)
      real*8 enu(3)
      real*8 n_img(3), n_img_enu(3), n_trg_enu(3)
      real*8 that(3), chat(3), nhat(3), vhat(3)
      real*8 enumat(3,3), xyz2enu(3,3)

!!      real*8 xyz2(3), vxyz2(3)

      integer, allocatable :: converge(:)
      integer totalconv, owidth, ofactor

      real*8 height, rcurv, vmag
      real*8 aa, bb,cc
      real*8 r2d,refhgt,hnadir
      integer pixel
      integer nearrangeflag
      real*8 beta, alpha, gamm
      real*8 costheta,sintheta,cosalpha
      real*8 arg,rminoraxis,rlatg,st,ct
      real*8 fraclat, fraclon
      real*4 z1,z2,demlat,demlon
      real*4 demmax
      real*8 cospsi
      integer line,iter,ind
      integer idemlat,idemlon,i_type,i,j!,i_cnt1,i_cnt2,i_loff,i_el,i_sl

      !!!Variables for cropped DEM
      integer udemwidth, udemlength
      integer ustartx, uendx
      integer ustarty, uendy
      double precision umin_lon, umax_lon
      double precision umin_lat, umax_lat
      double precision ufirstlat, ufirstlon
      double precision hgts(2)

      real*8 pi
      integer,parameter :: b1=1
      integer,parameter :: b2=1
      integer binarysearch

      !!Geometry objects
      type(ellipsoidType) :: elp
      type(pegType) :: peg
      type(pegtransType) :: ptm

      procedure(intpTemplate), pointer :: intp_dem => null()
      procedure(interpolateOrbit_f), pointer :: intp_orbit => null()

        !!!Set up DEM interpolation method
        if (method.eq.SINC_METHOD) then
          intp_dem => intp_sinc
        else if (method.eq.BILINEAR_METHOD) then
          intp_dem => intp_bilinear
        else if (method.eq.BICUBIC_METHOD) then
          intp_dem => intp_bicubic
        else if (method.eq.NEAREST_METHOD) then
          intp_dem => intp_nearest
        else if (method.eq.AKIMA_METHOD) then
          intp_dem => intp_akima
        else if (method.eq.BIQUINTIC_METHOD) then
            intp_dem => intp_biquintic
        else
          print *, 'Undefined interpolation method.'
          stop
        endif
        call prepareMethods(method)

        !!!Set up orbit interpolation method
        if (orbitmethod .eq. HERMITE_METHOD) then
            intp_orbit => interpolateWGS84Orbit_f

            if(orbit%nVectors .lt. 4) then
                print *, 'Need atleast 4 state vectors for using hermite polynomial interpolation'
                stop
            endif
            print *, 'Orbit interpolation method: hermite'
        else if (orbitmethod .eq. SCH_METHOD) then
            intp_orbit => interpolateSCHOrbit_f

            if(orbit%nVectors .lt. 4) then
                print *, 'Need atleast 4 state vectors for using SCH interpolation'
                stop
            endif
            print *, 'Orbit interpolation method: sch'
        else if (orbitmethod .eq. LEGENDRE_METHOD) then
            intp_orbit => interpolateLegendreOrbit_f

            if(orbit%nVectors .lt. 9) then
                print *, 'Need atleast 9 state vectors for using legendre polynomial interpolation'
                stop
            endif
            print *, 'Orbit interpolation method: legendre'
        else
            print *, 'Undefined orbit interpolation method.'
            stop
        endif

        ofactor = 2
        owidth = ofactor*width + 1
        pi = getPI()
        hgts(1) = MIN_H
        hgts(2) = MAX_H

        lineFile  = 0

        totalconv = 0

        height = 0.0d0
        min_lat = 10000.
        max_lat = -10000.
        min_lon = 10000.
        max_lon = -10000.

!$omp parallel
        if(omp_get_thread_num().eq.1) then
            write(6,*), 'Max threads used: ', omp_get_num_threads()
        end if
!$omp end parallel

        if ((slrngAccessor.eq.0).and.(r0.eq.0.0d0)) then
            print *, 'Both the slant range accessor and starting range are zero'
            stop
        endif

!c  allocate variable arrays
        allocate (lat(width))
        allocate (lon(width))
        allocate (z(width))
        allocate (zsch(width))
        allocate (rho(width))
        allocate (dopline(width))
        allocate (distance(width))
        allocate (losang(2*width))
        allocate (incang(2*width))
        allocate (elevang(width))
        allocate (converge(width)) !!PSA

        if (maskAccessor.gt.0) then
            allocate (omask(owidth))
            allocate (orng(owidth))
            allocate (mask(width))
            allocate (ctrack(owidth))
            allocate (oview(owidth))
        endif

!c  some constants
        refhgt=0
        r2d=180.d0/pi
        elp%r_a = major
        elp%r_e2 = eccentricitySquared

        !!!PSA - Keep track of near range issues
        nearrangeflag = 0

        !!!Determining the bbox of interest
        !!!For detailed explanation of steps - see main loop below
        line=1
        !!!Doppler for geometry (not carrier) is const / range variant only
        call getLine(dopAccessor, dopline, line)
        call getLine(slrngAccessor, rho, line)

        !!!First line
        do line=1,2
            tline = t0 + (line-1) * NAzlooks * (length-1.0d0)/prf
!!            stat = interpolateWGS84Orbit_f(orbit, tline, xyzsat, velsat)
            stat = intp_orbit(orbit, tline, xyzsat, velsat)
            if (stat.ne.0) then
                print *, 'Error getting statevector for bounds computation'
                exit
            endif
            vmag = norm(velsat)
            call unitvec(velsat, vhat)
            i_type = XYZ_2_LLH
            call latlon(elp, xyzsat, llhsat, i_type)
            height = llhsat(3)
            call tcnbasis(xyzsat, velsat, elp, that, chat, nhat)

            peg%r_lat = llhsat(1)
            peg%r_lon = llhsat(2)
            peg%r_hdg = peghdg
            call radar_to_xyz(elp, peg, ptm)
            rcurv = ptm%r_radcur

            do ind=1,2
                pixel = (ind-1)*(width-1) + 1
!                rng=r0 + (pixel-1) * Nrnglooks *rspace
                rng = rho(pixel)
                dopfact = (0.5d0 * wvl * dopline(pixel)/vmag) * rng

                do iter=1,2

                    !!PSA - SWOT specific near range check
                    !!If slant range vector doesn't hit ground, pick nadir point
                    if (rng .le. (llhsat(3)-hgts(iter)+1.0d0)) then
                        llh = llhsat
                        print *, 'Possible near nadir imaging'
                        nearrangeflag = 1
                    else
                        zsch(pixel) = hgts(iter)
                        aa =  height + rcurv
                        bb = rcurv + zsch(pixel)
                        costheta = 0.5*((aa/rng) + (rng/aa) - (bb/aa)*(bb/rng))
                        sintheta = sqrt(1.0d0 - costheta*costheta)
                        gamm = costheta * rng
                        alpha  = (dopfact - gamm * dot(nhat,vhat)) / dot(vhat,that)
                        beta  = -ilrl * sqrt(rng*rng*sintheta*sintheta - alpha*alpha)
                        delta = gamm * nhat + alpha * that + beta * chat
                        xyz = xyzsat + delta
                        i_type=XYZ_2_LLH
                        call latlon(elp,xyz,llh,i_type)
                    endif

                    min_lat = min(min_lat, llh(1)*r2d)
                    max_lat = max(max_lat, llh(1)*r2d)
                    min_lon = min(min_lon, llh(2)*r2d)
                    max_lon = max(max_lon, llh(2)*r2d)
                end do
            end do
        end do

        !!!Account for margins
        min_lon = min_lon - MARGIN
        max_lon = max_lon + MARGIN
        min_lat = min_lat - MARGIN
        max_lat = max_lat + MARGIN

        print *,'DEM parameters:'
        print *,'Dimensions: ',idemwidth,idemlength
        print *,'Top Left: ',firstlon,firstlat
        print *,'Spacing: ',deltalon,deltalat
        print *, 'Lon: ', firstlon, firstlon+(idemwidth-1)*deltalon
        print *, 'Lat: ', firstlat+(idemlength-1)*deltalat, firstlat

        print *, ' '
        print *, 'Estimated DEM bounds needed for global height range: '
        print *, 'Lon: ', min_lon, max_lon
        print *, 'Lat: ', min_lat, max_lat

        !!!!Compare with what has been provided as input
        umin_lon = max(min_lon, firstlon)
        if (min_lon .lt. firstlon) then
            print *, 'Warning: west limit may be insufficient for global height range'
        endif

        umax_lon = min(max_lon, firstlon + (idemwidth-1)*deltalon)
        if (max_lon .gt. (firstlon + (idemwidth-1)*deltalon)) then
            print *, 'Warning: east limit may be insufficient for global height range'
        endif

        umax_lat = min(max_lat, firstlat)
        if (max_lat .gt. firstlat) then
            print *, 'Warning: north limit may be insufficient for global height range'
        endif

        umin_lat = max(min_lat, firstlat + (idemlength-1)*deltalat)
        if (min_lat .lt. (firstlat + (idemlength-1)*deltalat)) then
            print *, 'Warning: south limit may be insufficient for global height range'
        endif

        !!!!Usable part of the DEM limits
        ustartx = int((umin_lon - firstlon)/deltalon)+1
        if (ustartx .lt. 1) ustartx = 1

        uendx = int((umax_lon-firstlon)/deltalon + 0.5d0)+1
        if (uendx.gt.idemwidth) uendx = idemwidth

        ustarty = int((umax_lat-firstlat)/deltalat)+1
        if (ustarty.lt.1) ustarty=1

        uendy = int((umin_lat-firstlat)/deltalat + 0.5) + 1
        if (uendy.gt.idemlength) ustarty=idemlength

        ufirstlon = firstlon + deltalon * (ustartx-1)
        ufirstlat = firstlat + deltalat * (ustarty-1)

        udemwidth = uendx - ustartx + 1
        udemlength = uendy - ustarty + 1

        print *, ' '
        print *, 'Actual DEM bounds used: '
        print *,'Dimensions: ',udemwidth,udemlength
        print *,'Top Left: ',ufirstlon,ufirstlat
        print *,'Spacing: ',deltalon,deltalat
        print *, 'Lon: ', ufirstlon, ufirstlon + deltalon*(udemwidth-1)
        print *, 'Lat: ', ufirstlat + deltalat * (udemlength-1), ufirstlat
        print *, 'Lines: ', ustarty, uendy
        print *, 'Pixels: ', ustartx, uendx

!c  allocate dem array
        allocate (dem(udemwidth,udemlength))
        allocate (demline(idemwidth))

        !!!Read the useful part of the DEM
        do j=1,udemlength
            lineFile = j + ustarty - 1
!            print *, 'Line: ', lineFile
            call getLine(demAccessor,demline,lineFile)
            dem(:,j) = demline(ustartx:uendx)
        enddo

        demmax = maxval(dem)
        print *, 'Max DEM height: ', demmax

        print *, 'Primary iterations: ', numiter
        print *, 'Secondary iterations: ', extraiter
        print *, 'Distance threshold : ', thresh

        !!Initialize range values
!!        do pixel=1,width
!!           rho(pixel) = r0 + rspace*(pixel-1)*Nrnglooks
!!        enddo

        height = 0.0d0
        min_lat = 10000.
        max_lat = -10000.
        min_lon = 10000.
        max_lon = -10000.

      !!!File for debugging
!!      open(31, file='distance',access='direct',recl=4*width,form='unformatted')

      do line=1, length         !c For each line

         !!!!Set up the geometry
         !!Step 1: Get satellite position
         !!Get time
         tline = t0 + Nazlooks*(line - 1.0d0)/prf
         !!Get state vector

!!         stat = interpolateLegendreOrbit_f(orbit, tline, xyz2, vxyz2)
!!         print *, 'Line: ', line
!!         print *, tline, xyz2, vxyz2
!!         stat = interpolateWGS84Orbit_f(orbit, tline, xyzsat, velsat)
         stat = intp_orbit(orbit, tline, xyzsat, velsat)
!!         print *, tline, xyzsat, velsat

         call unitvec(velsat, vhat)
         vmag = norm(velsat)
         !!vhat is unit vector along velocity
         !!vmag is the magnitude of the velocity

         !!Step 2: Get local radius of curvature along heading
         !!Convert satellite position to lat lon
         i_type = XYZ_2_LLH
         call latlon(elp, xyzsat, llhsat, i_type)
         height = llhsat(3)

!!         print *, 'Sat pos: ', line
!!         print *, llhsat(1)*r2d, llhsat(2)*r2d, llhsat(3)

         !!Step 3: Get TCN basis using satellite basis
         call tcnbasis(xyzsat, velsat, elp, that, chat, nhat)
         !!that is along local tangent to the planet
         !!chat is along the cross track direction
         !!nhat is along the local normal

         !!Step 4: Get Doppler information for the line
         !! For native doppler, this corresponds to doppler polynomial
         !! For zero doppler, its a constant zero polynomial
         call getLineSequential(dopAccessor, dopline, i_type)
!!         print *, 'VEL:', velsat
!!         print *, 'TCN:', xyzsat
!!         print *, that
!!         print *, chat
!!         print *, nhat
!!         print *, vhat
         !!Get the slant range
         call getLineSequential(slrngAccessor, rho, i_type)

         !!Step 4: Set up SCH basis right below the satellite
         peg%r_lat = llhsat(1)
         peg%r_lon = llhsat(2)
         peg%r_hdg = peghdg
         hnadir = 0.0d0

!!         print *, 'Heading: ', peghdg
         call radar_to_xyz(elp, peg, ptm)
         rcurv = ptm%r_radcur
         converge = 0
         z = 0.
         zsch = 0.

         if (mod(line,1000).eq.1) then
             print *, 'Processing line: ', line, vmag
             print *, 'Dopplers: ', dopline(1), dopline(width/2), dopline(width)
         endif

        !!Initialize lat,lon to middle of input DEM
        lat(:) = ufirstlat + 0.5d0*deltalat*udemlength
        lon(:) = ufirstlon + 0.05d0*deltalon*udemwidth

        !!PSA - SWOT specific near range check
        !!Computing nadir height
        if (nearrangeflag .ne. 0) then
            demlat=(llhsat(1)*r2d-ufirstlat)/deltalat+1
            demlon=(llhsat(2)*r2d-ufirstlon)/deltalon+1
            if(demlat.lt.1)demlat=1
            if(demlat.gt.udemlength-1)demlat=udemlength-1
            if(demlon.lt.1)demlon=1
            if(demlon.gt.udemwidth-1)demlon=udemwidth-1

            !!!!! This whole part can be put into a function
            idemlat=int(demlat)
            idemlon=int(demlon)
            fraclat=demlat-idemlat
            fraclon=demlon-idemlon
            hnadir = intp_dem(dem,idemlon,idemlat,fraclon,fraclat,udemwidth,udemlength)
        endif

         !!!!Start the iterations
         do iter=1,numiter+extraiter+1

            !$omp parallel do private(pixel,sch,beta,alpha,gamm) &
            !$omp private(i_type,llh,idemlat,idemlon,xyz,arg) &
            !$omp private(z1,z2,fraclat,fraclon,demlat,demlon) &
            !$omp private(llh_prev,xyz_prev,aa,bb,cc, rng) &
            !$omp private(costheta,sintheta,delta,dopfact)&
            !$omp shared(ufirstlat,ufirstlon,deltalat,deltalon)&
            !$omp shared(xyzsat,that,chat,nhat,vhat,peg,ptm)&
            !$omp shared(length,width,Nazlooks,height,r2d,dem) &
            !$omp shared(rcurv,rho,elp,lat,lon,z,zsch,line)&
            !$omp shared(extraiter,ilrl,iter,dopline,vmag,hnadir) &
            !$omp shared(distance,converge,thresh,numiter)&
            !$omp shared(udemwidth,udemlength,totalconv,wvl)
            do pixel=1,width
               rng = rho(pixel)
               dopfact = (0.5d0 * wvl * dopline(pixel)/vmag) * rng

               !!PSA - Check for near range issues
!!                if (nearrangeflag .ne. 0) then
!!                    if (rng .le. (llhsat(2)-hnadir+1.0d0)) then
!!                endif

               !! If pixel hasnt converged
               if(converge(pixel).eq.0) then

                   !!!!Use previous llh in degrees and meters
                   llh_prev(1) = lat(pixel)/r2d
                   llh_prev(2) = lon(pixel)/r2d
                   llh_prev(3) = z(pixel)

!!                  print *, 'ITER: ', iter
!!                  print *, 'PREV: ', lat(pixel), lon(pixel), z(pixel)

                  !!!!Solve for new position at height zsch
                  aa =  height + rcurv
                  bb = rcurv + zsch(pixel)

!!                  print *, aa, bb, rng
                  !!!!Normalize reasonably to avoid overflow
                  costheta = 0.5*((aa/rng) + (rng/aa) - (bb/aa)*(bb/rng))
                  sintheta = sqrt(1.0d0 - costheta*costheta)

!!                  print *, costheta, sintheta
                  !!Components along unit vectors

                  !!Vector from satellite to point on ground can be written as
                  !! vec(dr) = alpha * vec(that) + beta * vec(chat) + gamma *
                  !! vec(nhat)
                  gamm = costheta * rng
                  alpha  = (dopfact - gamm * dot(nhat,vhat)) / dot(vhat,that)
                  beta  = -ilrl * sqrt(rng*rng*sintheta*sintheta - alpha*alpha)
!!                  print *, alpha, beta, gamm

                  !!! xyz position of target
                  delta = gamm * nhat + alpha * that + beta * chat
                  xyz = xyzsat + delta

                  i_type=XYZ_2_LLH
                  call latlon(elp,xyz,llh,i_type)

!!                  print *, 'NOW:', llh(1)*r2d, llh(2)*r2d, llh(3)
                  !c  convert lat, lon, hgt to xyz coordinates
                  lat(pixel)=llh(1)*r2d
                  lon(pixel)=llh(2)*r2d
                  demlat=(lat(pixel)-ufirstlat)/deltalat+1
                  demlon=(lon(pixel)-ufirstlon)/deltalon+1
                  if(demlat.lt.1)demlat=1
                  if(demlat.gt.udemlength-1)demlat=udemlength-1
                  if(demlon.lt.1)demlon=1
                  if(demlon.gt.udemwidth-1)demlon=udemwidth-1

                  !!!!! This whole part can be put into a function
                  idemlat=int(demlat)
                  idemlon=int(demlon)
                  fraclat=demlat-idemlat
                  fraclon=demlon-idemlon
!!!                  z1=dem(idemlon,idemlat)*(1-fraclon)+dem(idemlon+1,idemlat)*fraclon
!!!                  z2=dem(idemlon,idemlat+1)*(1-fraclon)+dem(idemlon+1,idemlat+1)*fraclon
                  !!!Can change this to Akima
!!!                  z(pixel)=z1*(1-fraclat)+z2*fraclat

                  z(pixel) = intp_dem(dem,idemlon,idemlat,fraclon,fraclat,udemwidth,udemlength)
                  !!!!!! This whole part can be put into a function

                  if(z(pixel).lt.-500.0)z(pixel)=-500.0

                  ! given llh, where h = z(pixel,line) in WGS84, get the SCH height
                  llh(1) = lat(pixel)/r2d
                  llh(2) = lon(pixel)/r2d
                  llh(3) = z(pixel)

!!                  print *, 'UPDATED: ', lat(pixel), lon(pixel), z(pixel)
                  i_type = LLH_2_XYZ
                  call latlon(elp,xyz,llh,i_type)

                  i_type = XYZ_2_SCH
                  call convert_sch_to_xyz(ptm,sch,xyz,i_type)
                  ! print *, 'after = ', sch
!!                  print *, 'ZSCH:' , zsch(pixel), sch(3)
                  zsch(pixel) = sch(3)

                  !!!!Absolute distance
                  distance(pixel) = sqrt((xyz(1)-xyzsat(1))**2 +(xyz(2)-xyzsat(2))**2 + (xyz(3)-xyzsat(3))**2) - rng
!!                  print *, 'DIST: ', distance(pixel)
                  if(abs(distance(pixel)).le.thresh) then
                     zsch(pixel) = sch(3)
                     converge(pixel) = 1
                     totalconv = totalconv+1

                 else if(iter.gt.(numiter+1)) then

                    i_type=LLH_2_XYZ
                    call latlon(elp, xyz_prev,llh_prev,i_type)

                     xyz(1) = 0.5d0*(xyz_prev(1)+xyz(1))
                     xyz(2) = 0.5d0*(xyz_prev(2)+xyz(2))
                     xyz(3) = 0.5d0*(xyz_prev(3)+xyz(3))

                     !!!!Repopulate lat,lon,z
                     i_type=XYZ_2_LLH
                     call latlon(elp,xyz,llh,i_type)
                     lat(pixel) = llh(1)*r2d
                     lon(pixel) = llh(2)*r2d
                     z(pixel) = llh(3)

                     i_type=XYZ_2_SCH
                     call convert_sch_to_xyz(ptm,sch,xyz,i_type)
                     zsch(pixel) = sch(3)
                     !!!!Absolute distance
                     distance(pixel) = sqrt((xyz(1)-xyzsat(1))**2 +(xyz(2)-xyzsat(2))**2 + (xyz(3)-xyzsat(3))**2) - rng
                 endif
             endif

            end do
            !$omp end parallel do

         end do

         !!!!Final computation.
         !!!! The output points are exactly at range pixel
         !!!!distance from the satellite

         !$omp parallel do private(pixel,cosalpha) &
         !$omp private(xyz,llh,delta,rng,i_type,sch,aa,bb) &
         !$omp private(n_img,n_img_enu,n_trg_enu,cospsi) &
         !$omp private(alpha,beta,gamm,costheta,sintheta,dopfact) &
         !$omp private(enumat,xyz2enu,enu) &
         !$omp private(demlat,demlon,idemlat,idemlon,fraclat,fraclon)&
         !$omp shared(zsch,line,rcurv,rho,height,losang,width,velsat) &
         !$omp shared(peghdg,r2d,ilrl,lat,lon,z,xyzsat,distance,incang)&
         !$omp shared(elp,ptm,that,chat,vhat,nhat,vmag,dopline,wvl,dem)&
         !$omp shared(udemwidth,udemlength,ufirstlat,ufirstlon)&
         !$omp shared(deltalat,deltalon,elevang)
         do pixel=1,width

            rng = rho(pixel)
            dopfact = (0.5d0 * wvl * dopline(pixel)/vmag) * rng

            !!!!Solve for new position at height zsch
            aa =  height + rcurv
            bb = rcurv + zsch(pixel)

            costheta = 0.5*((aa/rng) + (rng/aa) - (bb/aa)*(bb/rng))
            sintheta = sqrt(1.0d0 - costheta*costheta)

            gamm = costheta * rng
            alpha  = (dopfact -gamm * dot(nhat,vhat)) / dot(vhat,that)
            beta  = -ilrl * sqrt(rng*rng*sintheta*sintheta - alpha*alpha)

            !!! xyz position of target
            delta = gamm * nhat + alpha * that + beta * chat
            xyz = xyzsat + delta

            i_type=XYZ_2_LLH
            call latlon(elp,xyz,llh,i_type)

            !!!!Copy into output arrays
            lat(pixel) = llh(1)*r2d
            lon(pixel) = llh(2)*r2d
            z(pixel)   = llh(3)

!!            distance(pixel) = ((xyz(1)-xyzsat(1))* velsat(1)+(xyz(2)-xyzsat(2))*velsat(2) + (xyz(3)-xyzsat(3))*velsat(3)) - dopfact * vmag
            distance(pixel) = sqrt((xyz(1)-xyzsat(1))**2 + (xyz(2)-xyzsat(2))**2 + (xyz(3)-xyzsat(3))**2) - rng

            !!!Computations in ENU coordinates around target
            call enubasis(llh(1), llh(2), enumat)
            xyz2enu = transpose(enumat)
            enu = matmul(xyz2enu,delta)

            cosalpha = abs(enu(3)) / norm(enu)

            !!!!LOS vectors
            losang(2*pixel-1) = acos(cosalpha)*r2d
            losang(2*pixel) = (atan2(-enu(2), -enu(1))-0.5*pi)*r2d
            elevang(pixel) = acos(costheta)*r2d

            !!!ctrack gets stored in zsch
            zsch(pixel) = rng * sintheta

            !!!!Get local incidence angle

            demlat=(lat(pixel)-ufirstlat)/deltalat+1
            demlon=(lon(pixel)-ufirstlon)/deltalon+1
            if(demlat.lt.2)demlat=2
            if(demlat.gt.udemlength-1)demlat=udemlength-1
            if(demlon.lt.2)demlon=2
            if(demlon.gt.udemwidth-1)demlon=udemwidth-1

            !!!!! This whole part can be put into a function
            idemlat=int(demlat)
            idemlon=int(demlon)
            fraclat=demlat-idemlat
            fraclon=demlon-idemlon

            !!!Slopex
            aa = intp_dem(dem,idemlon-1,idemlat,fraclon,fraclat,udemwidth,udemlength)
            bb  = intp_dem(dem,idemlon+1,idemlat,fraclon,fraclat,udemwidth,udemlength)
            gamm = lat(pixel)/r2d
            alpha = (bb-aa)* r2d / (2.0d0 * reast(elp, gamm) * deltalon)

            !!!Slopey
            aa = intp_dem(dem,idemlon,idemlat-1,fraclon,fraclat,udemwidth,udemlength)
            bb  = intp_dem(dem,idemlon,idemlat+1,fraclon,fraclat,udemwidth,udemlength)
            beta = (bb-aa)*r2d/(2.0d0 * rnorth(elp,gamm)*deltalat)

            enu = enu / norm(enu)
            costheta = (enu(1)*alpha + enu(2)*beta-enu(3))/sqrt(1.0d0+alpha*alpha+beta*beta)
            incang(2*pixel) = acos(costheta)*r2d

            !!!! Calculate psi angle between image plane and local slope
            call cross(delta, velsat, n_img)
            call unitvec(n_img, n_img)
            n_img_enu = matmul(xyz2enu, -ilrl*n_img)
            n_trg_enu = [-alpha, -beta, 1.0d0]
            cospsi = dot(n_trg_enu, n_img_enu) / (norm(n_trg_enu)*norm(n_img_enu))
            incang(2*pixel-1) = acos( cospsi )*r2d

            !!! Temporary hack needed by dense baseline in the
            !!! derivative computation. Todo: create two new layers
            !incang(2*pixel-1) = alpha                                                                                !incang(2*pixel) = beta

         end do
         !$omp end parallel do

         !c Maybe add hmin and hmax?
         min_lat = min(minval(lat), min_lat)
         max_lat = max(maxval(lat), max_lat)
         min_lon = min(minval(lon), min_lon)
         max_lon = max(maxval(lon), max_lon)
!!         write(31,rec=line)(distance(j),j=1,width)
         call setLineSequential(latAccessor, lat)
         call setLineSequential(lonAccessor, lon)
         call setLineSequential(heightAccessor, z)
         if(losAccessor.gt.0) then
             call setLineSequential(losAccessor,losang)
         endif

         if (incAccessor.gt.0) then
             call setLineSequential(incAccessor, incang)
         endif

         if (maskAccessor.gt.0) then
             ctrackmin = minval(zsch) - demmax
             ctrackmax = maxval(zsch) + demmax
             dctrack = (ctrackmax-ctrackmin)/(owidth-1.0d0)

             !!!Sort lat / lon by ctrack
             call InsertionSort(zsch, lat, lon, width)

             !!!Interpolate heights to regular ctrack grid

             !$omp parallel do private(pixel,llh,xyz,rng,aa,bb,i_type)&
             !$omp private(demlat,demlon,idemlat,idemlon,fraclat,fraclon)&
             !$omp shared(ctrackmin,ctrackmax,dctrack,dem,r2d)&
             !$omp shared(orng,owidth,lat,lon,ufirstlat,ufirstlon)&
             !$omp shared(deltalat,deltalon,udemlength,udemwidth)&
             !$omp shared(xyzsat,elp,ctrack,oview,nhat)
             do pixel=1,owidth
                aa = ctrackmin + (pixel-1)*dctrack
                ctrack(pixel) = aa
                i_type = binarysearch(zsch, width, aa)
!!                print *, line, pixel, aa, i_type
                if (i_type.eq.width) i_type = width-1
                if (i_type.eq.0) i_type=1

                !!!!Simple bi-linear interpolation
                fraclat = (aa - zsch(i_type)) / (zsch(i_type+1) - zsch(i_type))
                demlat = lat(i_type) + fraclat*(lat(i_type+1)-lat(i_type))
                demlon = lon(i_type) + fraclat*(lon(i_type+1)-lon(i_type))

                llh(1) = demlat/r2d
                llh(2) = demlon/r2d

                demlat=(demlat-ufirstlat)/deltalat+1
                demlon=(demlon-ufirstlon)/deltalon+1
                if(demlat.lt.2)demlat=2
                if(demlat.gt.udemlength-1)demlat=udemlength-1
                if(demlon.lt.2)demlon=2
                if(demlon.gt.udemwidth-1)demlon=udemwidth-1

                !!!!! This whole part can be put into a function
                idemlat=int(demlat)
                idemlon=int(demlon)
                fraclat=demlat-idemlat
                fraclon=demlon-idemlon
                llh(3) = intp_dem(dem,idemlon,idemlat,fraclon,fraclat,udemwidth,udemlength)
                i_type=LLH_2_XYZ
                call latlon(elp,xyz,llh,i_type)

                xyz = xyz - xyzsat
                bb = norm(xyz)
                orng(pixel) = bb
                aa = abs(sum(nhat*xyz))
                oview(pixel) = acos(aa/bb)*r2d
            end do
            !$omp end parallel do

            !!!Again sort in terms of slant range
            call InsertionSort(orng, ctrack, oview, owidth)

            mask = 0
            omask = 0
            aa = elevang(1)
            do pixel=2,width
                bb=elevang(pixel)
                if (bb.le.aa) then
                    mask(pixel) = 1
                else
                    aa = bb
                endif
            end do

            aa = elevang(width)
            do pixel=width-1,1,-1
                bb = elevang(pixel)
                if (bb.ge.aa) then
                    mask(pixel) = 1
                else
                    aa = bb
                endif
            end do

            !!!!If we wanted to work with shadow
            !!!!in cross track sorted coords
            !aa = oview(1)
            !do pixel=2,owidth
                !bb = oview(pixel)
                !if (bb.le.aa) then
                    !omask(pixel) = 1
                !else
                    !aa = bb
                !endif
            !enddo

            !aa = oview(width)
            !do pixel=width-1,1,-1
                !bb = oview(pixel)
                !if (bb.ge.aa) then
                    !omask(pixel) = 1
                !else
                    !aa = bb
                !endif
            !end do

            aa = ctrack(1)
            do pixel=2,width
                bb = ctrack(pixel)
                if ((bb.le.aa).and.(omask(pixel).lt.2)) then
                    omask(pixel) = omask(pixel) + 2
                else
                    aa = bb
                endif
            end do

            aa = ctrack(owidth)
            do pixel=owidth-1,1,-1
                bb = ctrack(pixel)
                if ((bb.ge.aa).and.(omask(pixel).lt.2)) then
                    omask(pixel) = omask(pixel) + 2
                else
                    aa = bb
                endif
            end do

            do pixel=1, owidth
                if (omask(pixel).gt.0) then
!!                    idemlat = nint((orng(pixel) - r0)/ (rspace * Nrnglooks))+1
                    idemlat = binarysearch(rho, width, orng(pixel))
                    if ((idemlat.ge.1) .and. (idemlat.le.width)) then
                        if (mask(idemlat) .lt. omask(pixel)) then
                            mask(idemlat) = mask(idemlat) + omask(pixel)
                        endif
                    endif
                endif
            enddo

            !!!!!If using shadow from ctrack coords
            !do pixel=1, owidth
                !if (omask(pixel).gt.0) then
                    !idemlat = nint((orng(pixel) - r0)/ (rspace * Nrnglooks))+1
                    !if ((idemlat.ge.1) .and. (idemlat.le.width)) then
                            !mask(idemlat) = omask(pixel)
                    !endif
                !endif
            !enddo

            call setLineSequential(MaskAccessor, mask)
        endif
      end do

      print *, 'Total convergence:', totalconv, ' out of ', width*length
      call unprepareMethods(method)
!!       close(31)

       if (maskAccessor.gt.0) then
          deallocate(omask)
          deallocate(orng)
          deallocate(mask)
          deallocate(ctrack)
          deallocate(oview)
       endif

       deallocate (demline)
       deallocate (converge)
       deallocate (distance)
       deallocate (lat)
       deallocate (lon)
       deallocate (z)
       deallocate (zsch)
       deallocate (rho)
       deallocate (dem)
       deallocate (losang)
       deallocate (incang)
       deallocate (elevang)
      end

      SUBROUTINE InsertionSort(a,b,c,num)
        REAL*8, DIMENSION(num) :: a,b,c
        REAL*8 :: tempa,tempb,tempc
        INTEGER :: i, j, num

        DO i = 2, num
            j = i - 1
            tempa = a(i)
            tempb = b(i)
            tempc = c(i)
            DO WHILE (j>=1 .AND. a(j)>tempa)
                a(j+1) = a(j)
                b(j+1) = b(j)
                c(j+1) = c(j)
                j = j - 1
            END DO
            a(j+1) = tempa
            b(j+1) = tempb
            c(j+1) = tempc
        END DO
     END SUBROUTINE InsertionSort

    function binarysearch(array, length, val)
        implicit none
        integer :: length
        real*8, dimension(length) :: array
        real*8 :: val

        integer :: binarysearch, ind

        integer :: left, middle, right

        left = 1
        right = length
        do
            if (left > right) then
                exit
            endif
            middle = nint((left+right) / 2.0)

            if (left .eq. (right-1)) then
                binarySearch = left
                return
            elseif (array(middle).le.val) then
                left = middle
            elseif (array(middle).gt.val) then
                right = middle
            end if
        end do

        binarysearch = left
    end function binarysearch


text:
C
      SUBROUTINE FPLATE (N, X, Y, F, RPAR, IPAR)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION Y(N), F(N)
      COMMON/TRANS/NX,NXM1,NY,NYM1,NDEMI,NACHS1,NACHS2,NDUMMY,
     &       OMEGA,STIFFN,DELX,USH4,FAC,WEIGHT
      DO 1 I=1,NX
      DO 1 J=1,NY
      K=I+NX*(J-1)
C -------- SECOND DERIVATIVE ----
      F(K)=Y(K+NDEMI)
C ------ CENTRAL POINT---
      UC=16.D0*Y(K)
      IF(I.GT.1)THEN
         UC=UC+Y(K)
         UC=UC-8.D0*Y(K-1)
      END IF
      IF(I.LT.NX)THEN
         UC=UC+Y(K)
         UC=UC-8.D0*Y(K+1)
      END IF
      IF(J.GT.1)THEN
         UC=UC+Y(K)
         UC=UC-8.D0*Y(K-NX)
      END IF
      IF(J.LT.NY)THEN
         UC=UC+Y(K)
         UC=UC-8.D0*Y(K+NX)
      END IF
      IF(I.GT.1 .AND.J.GT.1 )UC=UC+2.D0*Y(K-NX-1)
      IF(I.LT.NX.AND.J.GT.1 )UC=UC+2.D0*Y(K-NX+1)
      IF(I.GT.1 .AND.J.LT.NY)UC=UC+2.D0*Y(K+NX-1)
      IF(I.LT.NX.AND.J.LT.NY)UC=UC+2.D0*Y(K+NX+1)
      IF(I.GT.2)UC=UC+Y(K-2)
      IF(I.LT.NXM1)UC=UC+Y(K+2)
      IF(J.GT.2)UC=UC+Y(K-2*NX)
      IF(J.LT.NYM1)UC=UC+Y(K+2*NX)
      IF(J.EQ.NACHS1.OR.J.EQ.NACHS2)THEN
         XI=I*DELX
         FORCE=EXP(-5.D0*(X-XI-2.D0)**2)+EXP(-5.D0*(X-XI-5.D0)**2)
      ELSE
         FORCE=0.D0
      END IF
      F(K+NDEMI)=-OMEGA*Y(K+NDEMI)-FAC*UC+FORCE*WEIGHT
  1   CONTINUE
      RETURN
      END
c
      SUBROUTINE JPLATE(N,X,Y,DFY,LDFY,RPAR,IPAR)
C -------- JACOBIAN FOR PLATE PROBLEM --------
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION Y(N),DFY(LDFY,N)
      COMMON/TRANS/NX,NXM1,NY,NYM1,NDEMI,NACHS1,NACHS2,NDUMMY,
     &       OMEGA,STIFFN,DELX,USH4,FAC,WEIGHT
      DO 1 I=1,N
      DO 1 J=1,N
  1   DFY(I,J)=0.D0
C -------- LA BOUCLE -------
      DO 2 I=1,NX
      DO 2 J=1,NY
      K=I+NX*(J-1)
C -------- DERIVEE DEUXIEME ----
      DFY(K,K+NDEMI)=1.D0
C ------ POINT CENTRAL ---
      DFY(K+NDEMI,K)=-FAC*16.D0
      IF(I.GT.1)THEN
         DFY(K+NDEMI,K)=DFY(K+NDEMI,K)-FAC
         DFY(K+NDEMI,K-1)=FAC*8.D0
      END IF
      IF(I.LT.NX)THEN
         DFY(K+NDEMI,K)=DFY(K+NDEMI,K)-FAC
         DFY(K+NDEMI,K+1)=FAC*8.D0
      END IF
      IF(J.GT.1)THEN
         DFY(K+NDEMI,K)=DFY(K+NDEMI,K)-FAC
         DFY(K+NDEMI,K-NX)=FAC*8.D0
      END IF
      IF(J.LT.NY)THEN
         DFY(K+NDEMI,K)=DFY(K+NDEMI,K)-FAC
         DFY(K+NDEMI,K+NX)=FAC*8.D0
      END IF
      IF(I.GT.1 .AND.J.GT.1 )DFY(K+NDEMI,K-NX-1)=-FAC*2.D0
      IF(I.LT.NX.AND.J.GT.1 )DFY(K+NDEMI,K-NX+1)=-FAC*2.D0
      IF(I.GT.1 .AND.J.LT.NY)DFY(K+NDEMI,K+NX-1)=-FAC*2.D0
      IF(I.LT.NX.AND.J.LT.NY)DFY(K+NDEMI,K+NX+1)=-FAC*2.D0
      IF(I.GT.2)DFY(K+NDEMI,K-2)=-FAC
      IF(I.LT.NXM1)DFY(K+NDEMI,K+2)=-FAC
      IF(J.GT.2)DFY(K+NDEMI,K-2*NX)=-FAC
      IF(J.LT.NYM1)DFY(K+NDEMI,K+2*NX)=-FAC
      DFY(K+NDEMI,K+NDEMI)= -OMEGA
  2   CONTINUE
      RETURN
      END
c
      SUBROUTINE JPLATS(N,X,Y,DFY,LDFY,RPAR,IPAR)
C -------- JACOBIAN FOR PLATE PROBLEM, IN THE CASE WHERE THE OPTION
C -------- IWORK(9)=N/2 IS USED --------
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION Y(N),DFY(LDFY,N)
      COMMON/TRANS/NX,NXM1,NY,NYM1,NDEMI,NACHS1,NACHS2,NDUMMY,
     &       OMEGA,STIFFN,DELX,USH4,FAC,WEIGHT
C------ METTRE A ZERO -------
      DO 1 I=1,NDEMI
      DO 1 J=1,N
  1   DFY(I,J)=0.D0
C -------- LA BOUCLE -------
      DO 2 I=1,NX
      DO 2 J=1,NY
      K=I+NX*(J-1)
C ------ POINT CENTRAL ---
      DFY(K,K)=-FAC*16.D0
      IF(I.GT.1)THEN
         DFY(K,K)=DFY(K,K)-FAC
         DFY(K,K-1)=FAC*8.D0
      END IF
      IF(I.LT.NX)THEN
         DFY(K,K)=DFY(K,K)-FAC
         DFY(K,K+1)=FAC*8.D0
      END IF
      IF(J.GT.1)THEN
         DFY(K,K)=DFY(K,K)-FAC
         DFY(K,K-NX)=FAC*8.D0
      END IF
      IF(J.LT.NY)THEN
         DFY(K,K)=DFY(K,K)-FAC
         DFY(K,K+NX)=FAC*8.D0
      END IF
      IF(I.GT.1 .AND.J.GT.1 )DFY(K,K-NX-1)=-FAC*2.D0
      IF(I.LT.NX.AND.J.GT.1 )DFY(K,K-NX+1)=-FAC*2.D0
      IF(I.GT.1 .AND.J.LT.NY)DFY(K,K+NX-1)=-FAC*2.D0
      IF(I.LT.NX.AND.J.LT.NY)DFY(K,K+NX+1)=-FAC*2.D0
      IF(I.GT.2)DFY(K,K-2)=-FAC
      IF(I.LT.NXM1)DFY(K,K+2)=-FAC
      IF(J.GT.2)DFY(K,K-2*NX)=-FAC
      IF(J.LT.NYM1)DFY(K,K+2*NX)=-FAC
      DFY(K,K+NDEMI)= -OMEGA
  2   CONTINUE
      RETURN
      END
c
      SUBROUTINE JPLATSB(N,X,Y,DFY,LDFY,RPAR,IPAR)
C -------- JACOBIAN FOR PLATE PROBLEM, IN THE CASE WHERE THE OPTION
C -------- IWORK(9)=N/2 AND THE BANDED STRUCTURE IS USED --------
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION Y(N),DFY(LDFY,N)
      COMMON/TRANS/NX,NXM1,NY,NYM1,NDEMI,NACHS1,NACHS2,NDUMMY,
     &       OMEGA,STIFFN,DELX,USH4,FAC,WEIGHT
C------ METTRE A ZERO -------
      DO 1 I=1,LDFY
      DO 1 J=1,N
  1   DFY(I,J)=0.D0
      MU=2*NX+1
      FAC2=FAC*2.0D0
      FAC8=FAC*8.0D0
      FAC16=FAC*16.0D0
C -------- LA BOUCLE -------
      DO 2 I=1,NX
      DO 2 J=1,NY
      K=I+NX*(J-1)
C ------ POINT CENTRAL ---
      DFY(MU,K)=-FAC16
      IF(I.GT.1)THEN
         DFY(MU,K)=DFY(MU,K)-FAC
         DFY(MU+1,K-1)=FAC8
      END IF
      IF(I.LT.NX)THEN
         DFY(MU,K)=DFY(MU,K)-FAC
         DFY(MU-1,K+1)=FAC8
      END IF
      IF(J.GT.1)THEN
         DFY(MU,K)=DFY(MU,K)-FAC
         DFY(MU+NX,K-NX)=FAC8
      END IF
      IF(J.LT.NY)THEN
         DFY(MU,K)=DFY(MU,K)-FAC
         DFY(MU-NX,K+NX)=FAC8
      END IF
      IF(I.GT.1 .AND.J.GT.1 )DFY(MU+NX+1,K-NX-1)=-FAC2
      IF(I.LT.NX.AND.J.GT.1 )DFY(MU+NX-1,K-NX+1)=-FAC2
      IF(I.GT.1 .AND.J.LT.NY)DFY(MU-NX+1,K+NX-1)=-FAC2
      IF(I.LT.NX.AND.J.LT.NY)DFY(MU-NX-1,K+NX+1)=-FAC2
      IF(I.GT.2)DFY(MU+2,K-2)=-FAC
      IF(I.LT.NXM1)DFY(MU-2,K+2)=-FAC
      IF(J.GT.2)DFY(MU+2*NX,K-2*NX)=-FAC
      IF(J.LT.NYM1)DFY(MU-2*NX,K+2*NX)=-FAC
      DFY(MU,K+NDEMI)= -OMEGA
  2   CONTINUE
      RETURN
      END

      SUBROUTINE XPLATE(N,X,Y,FX,RPAR,IPAR)
C -------- DERIVATIVE OF F WITH RESPECT TO X --------
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION Y(N),FX(N)
      COMMON/TRANS/NX,NXM1,NY,NYM1,NDEMI,NACHS1,NACHS2,NDUMMY,
     &       OMEGA,STIFFN,DELX,USH4,FAC,WEIGHT
C -------- METTRE A ZERO --------
      DO 1 I=1,N
   1  FX(I)=0.D0
C -------- LA BOUCLE -------
      DO 4 I=1,NX
      DO 4 J=1,NY
      K=I+NX*(J-1)
      IF(J.EQ.NACHS1.OR.J.EQ.NACHS2)THEN
         XI=I*DELX
         FX(K+NDEMI)=WEIGHT*(
     &             -10.D0*(X-XI-2.D0)*EXP(-5.D0*(X-XI-2.D0)**2)
     &             -10.D0*(X-XI-5.D0)*EXP(-5.D0*(X-XI-5.D0)**2))
      END IF
   4  CONTINUE
      RETURN
      END

max_stars_repo_name: cpmech/hwode

max_issues_repo_name: cpmech/hwode

max_forks_repo_name: cpmech/hwode


text:
      SUBROUTINE FLTRAN (A,B,C, MAT)

C  Given a set of attitude numbers for a plane, FLTRAN returns a
C  4 X 4 transformation matrix which rotates the plane so as to be
C  parallel to the z=0 plane.  This is useful for visualizing a
C  planar 3D primitive.

      REAL     A,B,C, MAT(4,4), ZZ,UU, VECL, SATAN2

      ZZ = 0.0
      UU = 1.0

      CALL EBLTM3 (ZZ,ZZ,ZZ, ZZ,ZZ,ZZ,
     1     -SATAN2(-B,C), -ASIN(A/VECL(A,B,C)),ZZ, UU,UU,UU, MAT)

      END


text:
!==============================================================================!
  subroutine Load_Geo(grid, this_proc)
!------------------------------------------------------------------------------!
!   Reads:  name.geo                                                           !
!----------------------------------[Modules]-----------------------------------!
  use all_mod, only: delta, WallDs
  use Grid_Mod
!------------------------------------------------------------------------------!
  implicit none
!---------------------------------[Arguments]----------------------------------!
  type(Grid_Type) :: grid
  integer         :: this_proc
!-----------------------------------[Locals]-----------------------------------!
  integer           :: c, s
  character(len=80) :: name_in
!==============================================================================!

  !----------------------------!
  !     Read the file with     !
  !   geometrical quantities   !
  !----------------------------!
  call Name_File(this_proc, name_in, '.geo', len_trim('.geo'))
  open(9, file=name_in, FORM='unformatted')
  if(this_proc < 2) write(*,*) '# Now reading the file:', name_in

  read(9) (grid % xc(c), c = 1, grid % n_cells)
  read(9) (grid % yc(c), c = 1, grid % n_cells)
  read(9) (grid % zc(c), c = 1, grid % n_cells)

  read(9) (grid % xc(c), c=-1,-grid % n_bnd_cells,-1)
  read(9) (grid % yc(c), c=-1,-grid % n_bnd_cells,-1)
  read(9) (grid % zc(c), c=-1,-grid % n_bnd_cells,-1)

  read(9) (grid % vol(c), c = 1, grid % n_cells)
  read(9) (delta(c),      c = 1, grid % n_cells)
  read(9) (WallDs(c),     c = 1, grid % n_cells)

  read(9) (grid % sx(s), s = 1, grid % n_faces)
  read(9) (grid % sy(s), s = 1, grid % n_faces)
  read(9) (grid % sz(s), s = 1, grid % n_faces)

  read(9) (grid % dx(s), s = 1, grid % n_faces)
  read(9) (grid % dy(s), s = 1, grid % n_faces)
  read(9) (grid % dz(s), s = 1, grid % n_faces)

  read(9) (grid % f(s), s = 1, grid % n_faces)

  read(9) (grid % xf(s), s = 1, grid % n_faces)
  read(9) (grid % yf(s), s = 1, grid % n_faces)
  read(9) (grid % zf(s), s = 1, grid % n_faces)

  close(9)

  end subroutine

max_stars_repo_name: mhadziabdic/T-FlowS

max_issues_repo_name: mhadziabdic/T-FlowS

max_forks_repo_name: mhadziabdic/T-FlowS


text:
CDECK RATNUC
      SUBROUTINE  RATNUC
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      SAVE
C     ------------------------------------------------------------------
      COMMON  /VECTA/  VA(401),SA(401),PA(401),PEA(401),TA(401)
      COMMON  /VECTA/  RA(401),TAUA(401),FRA(401),FCA(401)
      COMMON  /VECTA/  QNUC(401),RATX(401)
      COMMON  /VECTA/  CPA(401),DADA(401)
      COMMON  /VECTA/  CC(401),UA(401),QT(401),EA(401)
      COMMON  /VECTA/  FCT(401),TAUC(401)
      COMMON  /VECTA/  PVS(401)
      COMMON  /VECTA/  OPACA(401),DOPDVA(401),DOPDTA(401)
      COMMON  /VECTAI/ KONV(401)
C     ------------------------------------------------------------------
      COMMON  /COMPOS/  XXX,YYY,ZZZ
      COMMON  /COMPOS/  XXA(401),YYA(401),ZZA(401)
      COMMON  /COMPOS/  EFF(401),DRODEF(401)
      COMMON  /COMPOSI/ JPOINT
C     ------------------------------------------------------------------
      COMMON  /COMPROG/  FACMIX,FLUM,PEXTER,TAUEX,PSTART
      COMMON  /COMPROI/  NP,NPP1
C     ------------------------------------------------------------------
      DO 1 I=2,NPP1
         QNUC(I)=0.0D0
         RATX(I)=0.0D0
         X=XXA(I)
         Y=YYA(I)
         ZCN=ZZA(I)
         EPP=0.D0
         ECN=0.D0
         T=TA(I)/1.D+6
         T13=T**(1.D0/3.D0)
         RO=1.D0/VA(I)
         IF(T.GE.1.0D0) THEN
            IF(X.GE.1.D-20) THEN
               SQTR=SQRT(RO/T**3)
               W=1.22D16*(X/(1.D0+X))*(1.D0+SQTR)/SQRT(T13)
     1              *EXP(-102.6D0/T13)
               ALPHA=5.48D17*(Y*0.25D0/X)**2*EXP(-100.D0/T13)
C     changed to avoid data dependencies
C               IF(ALPHA.LE.1.D5) GAMA=SQRT(ALPHA*(2.D0+ALPHA))-ALPHA
C               IF(ALPHA.GT.1.D5) GAMA=1.D0-0.5D0/ALPHA
C     1              *(1.D0-1.D0/ALPHA)
               IF(ALPHA.LE.1.D5) THEN
                  GAMA=SQRT(ALPHA*(2.D0+ALPHA))-ALPHA
                  ELSE
                     GAMA=1.D0-0.5D0/ALPHA *(1.D0-1.D0/ALPHA)
               ENDIF
               PSI=1.D0+GAMA*(0.959D0+0.47D0*W)/(1.D0+W)
               G11=1.D0+0.0012D0*T13+0.0078D0*T13**2+0.0006D0*T
               F11=1.D0+0.25D0*SQTR
               EPP=2.206D6*RO*X**2/T13**2*EXP(-33.81D0/T13)
     1              *F11*G11*PSI
               IF(T.GE.6.D0) THEN
                  IF(ZCN.GE.1.D-10) THEN
                     F14=1.D0+1.75D0*SQTR
                     G14=1.D0+0.0027D0*T13-0.0037D0*T13**2-0.0007D0*T
                     ECN=7.94D27*RO*X*ZCN/T13**2*EXP(-152.313D0/T13)
     1                    *F14*G14
                  ENDIF
               ENDIF
 2             CONTINUE
               E=EPP+ECN
               RAT=EPP/6256.D0+ECN/6035.D0
               RAT=RAT/1.D15
               QNUC(I)=E
               RATX(I)=RAT
            ENDIF
         ENDIF
    1 CONTINUE
      QNUC(1)=0.0D0
      RATX(1)=0.0D0
      RETURN
      END


module stppm2_5mod

!$$$ module documentation block
!           .      .    .                                       .
! module:   stppm2_5mod    module for stppm2_5 and its tangent linear stppm2_5_tl
!  pgrmmr:
!
! abstract: module for stppm2_5 and its tangent linear stppm2_5_tl
!
! program history log:
!   2010-10-15  pagowski - initial code for pm2_5
!   2016-05-18  guo     - replaced ob_type with polymorphic obsNode through type casting
!
! subroutines included:
!   sub stppm2_5
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block
  use m_obsNode  , only: obsNode
  use m_pm2_5Node, only: pm2_5Node
  use m_pm2_5Node, only: pm2_5Node_typecast
  use m_pm2_5Node, only: pm2_5Node_nextcast
  implicit none

  private
  public stppm2_5

contains

  subroutine stppm2_5(pm2_5head,rval,sval,out,sges,nstep)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    stppm2_5        calcuate penalty and stepsize from q
!                            with addition of nonlinear qc.
!   prgmmr: derber           org: np23                date: 1991-02-26
!
! abstract: calculate penalty and contribution to stepsize from pm2_5
!           using nonlinear qc.
!
! program history log:
!   2010-10-15  pagowski - initial code for pm2_5
!
!   input argument list:
!     pm2_5head
!     rpm2_5       - search direction for pm2_5
!     spm2_5       - analysis increment for pm2_5
!     sges     - stepsize estimates (nstep)
!     nstep    - number of stepsize estimates (== 0 means use outer iteration values)
!
!   output argument list:
!     out(1:nstep)   - contribution of penalty from pm2_5 sges(1:nstep)
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
    use kinds, only: r_kind,i_kind,r_quad
    use qcmod, only: nlnqc_iter,varqc_iter
    use constants, only: half,one,two,tiny_r_kind,cg_term,zero_quad,max_varname_length,zero
    use gsi_bundlemod, only: gsi_bundle
    use gsi_bundlemod, only: gsi_bundlegetpointer
    use gridmod, only: cmaq_regional,wrf_mass_regional
    USE chemmod, ONLY: d_2_5,s_2_5,nh4_mfac,oc_mfac,laeroana_gocart
    implicit none

! declare passed variables
    class(obsNode), pointer             ,intent(in   ) :: pm2_5head
    integer(i_kind)                     ,intent(in   ) :: nstep
    real(r_quad),dimension(max(1,nstep)),intent(inout) :: out
    type(gsi_bundle)                    ,intent(in   ) :: rval,sval
    real(r_kind),dimension(max(1,nstep)),intent(in   ) :: sges

! declare local variables
    integer(i_kind) ier,istatus
    integer(i_kind) j1,j2,j3,j4,j5,j6,j7,j8,kk
    real(r_kind) cg_pm2_5,val,val2,wgross,wnotgross,pm2_5_pg
    real(r_kind),dimension(max(1,nstep))::pen
    real(r_kind) w1,w2,w3,w4,w5,w6,w7,w8,qq
    real(r_kind),pointer,dimension(:):: rpm2_5,spm2_5
    type(pm2_5Node), pointer :: pm2_5ptr

    character(len=max_varname_length) :: aeroname

    out=zero_quad

!   If no pm2_5 data return
    if(.not. associated(pm2_5head))return

    if (cmaq_regional .or. (wrf_mass_regional .and. .not. laeroana_gocart)) then

       ier=0
       call gsi_bundlegetpointer(sval,'pm2_5',spm2_5,istatus);ier=ier+istatus
       call gsi_bundlegetpointer(rval,'pm2_5',rpm2_5,istatus);ier=ier+istatus

       if(ier/=0) return

       !pm2_5ptr => pm2_5head
       pm2_5ptr => pm2_5Node_typecast(pm2_5head)
       do while (associated(pm2_5ptr))
          if(pm2_5ptr%luse)then
             if(nstep > 0)then
                j1=pm2_5ptr%ij(1)
                j2=pm2_5ptr%ij(2)
                j3=pm2_5ptr%ij(3)
                j4=pm2_5ptr%ij(4)
                j5=pm2_5ptr%ij(5)
                j6=pm2_5ptr%ij(6)
                j7=pm2_5ptr%ij(7)
                j8=pm2_5ptr%ij(8)
                w1=pm2_5ptr%wij(1)
                w2=pm2_5ptr%wij(2)
                w3=pm2_5ptr%wij(3)
                w4=pm2_5ptr%wij(4)
                w5=pm2_5ptr%wij(5)
                w6=pm2_5ptr%wij(6)
                w7=pm2_5ptr%wij(7)
                w8=pm2_5ptr%wij(8)

                val= w1* rpm2_5(j1)+w2* rpm2_5(j2)+w3* rpm2_5(j3)+w4* rpm2_5(j4)+ &
                     w5* rpm2_5(j5)+w6* rpm2_5(j6)+w7* rpm2_5(j7)+w8* rpm2_5(j8)
                val2=w1* spm2_5(j1)+w2* spm2_5(j2)+w3* spm2_5(j3)+w4* spm2_5(j4)+ &
                     w5* spm2_5(j5)+w6* spm2_5(j6)+w7* spm2_5(j7)+w8* spm2_5(j8)-pm2_5ptr%res

                do kk=1,nstep
                   qq=val2+sges(kk)*val
                   pen(kk)=qq*qq*pm2_5ptr%err2
                end do
             else
                pen(1)=pm2_5ptr%res*pm2_5ptr%res*pm2_5ptr%err2
             end if

!  modify penalty term if nonlinear qc

             if (nlnqc_iter .and. pm2_5ptr%pg > tiny_r_kind .and. &
                  pm2_5ptr%b  > tiny_r_kind) then
                pm2_5_pg=pm2_5ptr%pg*varqc_iter
                cg_pm2_5=cg_term/pm2_5ptr%b
                wnotgross= one-pm2_5_pg
                wgross = pm2_5_pg*cg_pm2_5/wnotgross
                do kk=1,max(1,nstep)
                   pen(kk)= -two*log((exp(-half*pen(kk))+wgross)/(one+wgross))
                end do
             endif

             out(1) = out(1)+pen(1)*pm2_5ptr%raterr2
             do kk=2,nstep
                out(kk) = out(kk)+(pen(kk)-pen(1))*pm2_5ptr%raterr2
             end do
          end if

          !pm2_5ptr => pm2_5ptr%llpoint
          pm2_5ptr => pm2_5Node_nextcast(pm2_5ptr)

       end do

    endif

    if (wrf_mass_regional .and. laeroana_gocart) then

       !pm2_5ptr => pm2_5head
       pm2_5ptr => pm2_5Node_typecast(pm2_5head)

       do while (associated(pm2_5ptr))
          if(pm2_5ptr%luse)then
             if(nstep > 0)then
                j1=pm2_5ptr%ij(1)
                j2=pm2_5ptr%ij(2)
                j3=pm2_5ptr%ij(3)
                j4=pm2_5ptr%ij(4)
                j5=pm2_5ptr%ij(5)
                j6=pm2_5ptr%ij(6)
                j7=pm2_5ptr%ij(7)
                j8=pm2_5ptr%ij(8)
                w1=pm2_5ptr%wij(1)
                w2=pm2_5ptr%wij(2)
                w3=pm2_5ptr%wij(3)
                w4=pm2_5ptr%wij(4)
                w5=pm2_5ptr%wij(5)
                w6=pm2_5ptr%wij(6)
                w7=pm2_5ptr%wij(7)
                w8=pm2_5ptr%wij(8)

                istatus=0

                val2=-pm2_5ptr%res
                val=zero

                aeroname='bc1'

                call gsi_bundlegetpointer(sval,trim(aeroname),spm2_5,istatus)
                if(istatus /= 0) then
                   write(6,*) 'error gsi_bundlegetpointer in stppm2_5 for ',&
                        aeroname
                   call stop2(456)
                endif

                call gsi_bundlegetpointer(rval,trim(aeroname),rpm2_5,istatus)

                if(istatus /= 0) then
                   write(6,*) 'error gsi_bundlegetpointer in stppm2_5 for ',&
                        aeroname
                   call stop2(457)
                endif

                val= w1* rpm2_5(j1)+w2* rpm2_5(j2)+w3* rpm2_5(j3)+w4* rpm2_5(j4)+ &
                     w5* rpm2_5(j5)+w6* rpm2_5(j6)+w7* rpm2_5(j7)+w8* rpm2_5(j8)+val
                val2=w1* spm2_5(j1)+w2* spm2_5(j2)+w3* spm2_5(j3)+w4* spm2_5(j4)+ &
                     w5* spm2_5(j5)+w6* spm2_5(j6)+w7* spm2_5(j7)+w8* spm2_5(j8)+val2

                nullify(spm2_5,rpm2_5)

                aeroname='bc2'

                call gsi_bundlegetpointer(sval,trim(aeroname),spm2_5,istatus)
                if(istatus /= 0) then
                   write(6,*) 'error gsi_bundlegetpointer in stppm2_5 for ',&
                        aeroname
                   call stop2(456)
                endif

                call gsi_bundlegetpointer(rval,trim(aeroname),rpm2_5,istatus)

                if(istatus /= 0) then
                   write(6,*) 'error gsi_bundlegetpointer in stppm2_5 for ',&
                        aeroname
                   call stop2(457)
                endif

                val= w1* rpm2_5(j1)+w2* rpm2_5(j2)+w3* rpm2_5(j3)+w4* rpm2_5(j4)+ &
                     w5* rpm2_5(j5)+w6* rpm2_5(j6)+w7* rpm2_5(j7)+w8* rpm2_5(j8)+val
                val2=w1* spm2_5(j1)+w2* spm2_5(j2)+w3* spm2_5(j3)+w4* spm2_5(j4)+ &
                     w5* spm2_5(j5)+w6* spm2_5(j6)+w7* spm2_5(j7)+w8* spm2_5(j8)+val2

                nullify(spm2_5,rpm2_5)

                aeroname='sulf'

                call gsi_bundlegetpointer(sval,trim(aeroname),spm2_5,istatus)
                if(istatus /= 0) then
                   write(6,*) 'error gsi_bundlegetpointer in stppm2_5 for ',&
                        aeroname
                   call stop2(456)
                endif

                call gsi_bundlegetpointer(rval,trim(aeroname),rpm2_5,istatus)

                if(istatus /= 0) then
                   write(6,*) 'error gsi_bundlegetpointer in stppm2_5 for ',&
                        aeroname
                   call stop2(457)
                endif

                val= w1* rpm2_5(j1)+w2* rpm2_5(j2)+w3* rpm2_5(j3)+w4* rpm2_5(j4)+ &
                     w5* rpm2_5(j5)+w6* rpm2_5(j6)+w7* rpm2_5(j7)+w8* rpm2_5(j8)+val
                val2=(w1* spm2_5(j1)+w2* spm2_5(j2)+w3* spm2_5(j3)+w4* spm2_5(j4)+ &
                     w5* spm2_5(j5)+w6* spm2_5(j6)+w7* spm2_5(j7)+w8* spm2_5(j8)) * nh4_mfac +val2

                nullify(spm2_5,rpm2_5)

                aeroname='p25'

                call gsi_bundlegetpointer(sval,trim(aeroname),spm2_5,istatus)
                if(istatus /= 0) then
                   write(6,*) 'error gsi_bundlegetpointer in stppm2_5 for ',&
                        aeroname
                   call stop2(456)
                endif

                call gsi_bundlegetpointer(rval,trim(aeroname),rpm2_5,istatus)

                if(istatus /= 0) then
                   write(6,*) 'error gsi_bundlegetpointer in stppm2_5 for ',&
                        aeroname
                   call stop2(457)
                endif

                val= w1* rpm2_5(j1)+w2* rpm2_5(j2)+w3* rpm2_5(j3)+w4* rpm2_5(j4)+ &
                     w5* rpm2_5(j5)+w6* rpm2_5(j6)+w7* rpm2_5(j7)+w8* rpm2_5(j8)+val
                val2=w1* spm2_5(j1)+w2* spm2_5(j2)+w3* spm2_5(j3)+w4* spm2_5(j4)+ &
                     w5* spm2_5(j5)+w6* spm2_5(j6)+w7* spm2_5(j7)+w8* spm2_5(j8)+val2

                nullify(spm2_5,rpm2_5)

                aeroname='oc1'

                call gsi_bundlegetpointer(sval,trim(aeroname),spm2_5,istatus)
                if(istatus /= 0) then
                   write(6,*) 'error gsi_bundlegetpointer in stppm2_5 for ',&
                        aeroname
                   call stop2(456)
                endif

                call gsi_bundlegetpointer(rval,trim(aeroname),rpm2_5,istatus)

                if(istatus /= 0) then
                   write(6,*) 'error gsi_bundlegetpointer in stppm2_5 for ',&
                        aeroname
                   call stop2(457)
                endif

                val= w1* rpm2_5(j1)+w2* rpm2_5(j2)+w3* rpm2_5(j3)+w4* rpm2_5(j4)+ &
                     w5* rpm2_5(j5)+w6* rpm2_5(j6)+w7* rpm2_5(j7)+w8* rpm2_5(j8)+val
                val2=(w1* spm2_5(j1)+w2* spm2_5(j2)+w3* spm2_5(j3)+w4* spm2_5(j4)+ &
                     w5* spm2_5(j5)+w6* spm2_5(j6)+w7* spm2_5(j7)+w8* spm2_5(j8)) * oc_mfac +val2

                nullify(spm2_5,rpm2_5)

                aeroname='oc2'

                call gsi_bundlegetpointer(sval,trim(aeroname),spm2_5,istatus)
                if(istatus /= 0) then
                   write(6,*) 'error gsi_bundlegetpointer in stppm2_5 for ',&
                        aeroname
                   call stop2(456)
                endif

                call gsi_bundlegetpointer(rval,trim(aeroname),rpm2_5,istatus)

                if(istatus /= 0) then
                   write(6,*) 'error gsi_bundlegetpointer in stppm2_5 for ',&
                        aeroname
                   call stop2(457)
                endif

                val= w1* rpm2_5(j1)+w2* rpm2_5(j2)+w3* rpm2_5(j3)+w4* rpm2_5(j4)+ &
                     w5* rpm2_5(j5)+w6* rpm2_5(j6)+w7* rpm2_5(j7)+w8* rpm2_5(j8)+val
                val2=(w1* spm2_5(j1)+w2* spm2_5(j2)+w3* spm2_5(j3)+w4* spm2_5(j4)+ &
                     w5* spm2_5(j5)+w6* spm2_5(j6)+w7* spm2_5(j7)+w8* spm2_5(j8)) * oc_mfac +val2

                nullify(spm2_5,rpm2_5)

                aeroname='seas1'

                call gsi_bundlegetpointer(sval,trim(aeroname),spm2_5,istatus)
                if(istatus /= 0) then
                   write(6,*) 'error gsi_bundlegetpointer in stppm2_5 for ',&
                        aeroname
                   call stop2(456)
                endif

                call gsi_bundlegetpointer(rval,trim(aeroname),rpm2_5,istatus)

                if(istatus /= 0) then
                   write(6,*) 'error gsi_bundlegetpointer in stppm2_5 for ',&
                        aeroname
                   call stop2(457)
                endif

                val= w1* rpm2_5(j1)+w2* rpm2_5(j2)+w3* rpm2_5(j3)+w4* rpm2_5(j4)+ &
                     w5* rpm2_5(j5)+w6* rpm2_5(j6)+w7* rpm2_5(j7)+w8* rpm2_5(j8)+val
                val2=w1* spm2_5(j1)+w2* spm2_5(j2)+w3* spm2_5(j3)+w4* spm2_5(j4)+ &
                     w5* spm2_5(j5)+w6* spm2_5(j6)+w7* spm2_5(j7)+w8* spm2_5(j8)+val2

                nullify(spm2_5,rpm2_5)

                aeroname='seas2'

                call gsi_bundlegetpointer(sval,trim(aeroname),spm2_5,istatus)
                if(istatus /= 0) then
                   write(6,*) 'error gsi_bundlegetpointer in stppm2_5 for ',&
                        aeroname
                   call stop2(456)
                endif

                call gsi_bundlegetpointer(rval,trim(aeroname),rpm2_5,istatus)

                if(istatus /= 0) then
                   write(6,*) 'error gsi_bundlegetpointer in stppm2_5 for ',&
                        aeroname
                   call stop2(457)
                endif

                val= w1* rpm2_5(j1)+w2* rpm2_5(j2)+w3* rpm2_5(j3)+w4* rpm2_5(j4)+ &
                     w5* rpm2_5(j5)+w6* rpm2_5(j6)+w7* rpm2_5(j7)+w8* rpm2_5(j8)+val
                val2=(w1* spm2_5(j1)+w2* spm2_5(j2)+w3* spm2_5(j3)+w4* spm2_5(j4)+ &
                     w5* spm2_5(j5)+w6* spm2_5(j6)+w7* spm2_5(j7)+w8* spm2_5(j8)) * s_2_5 +val2

                nullify(spm2_5,rpm2_5)

                aeroname='dust1'

                call gsi_bundlegetpointer(sval,trim(aeroname),spm2_5,istatus)
                if(istatus /= 0) then
                   write(6,*) 'error gsi_bundlegetpointer in stppm2_5 for ',&
                        aeroname
                   call stop2(456)
                endif

                call gsi_bundlegetpointer(rval,trim(aeroname),rpm2_5,istatus)

                if(istatus /= 0) then
                   write(6,*) 'error gsi_bundlegetpointer in stppm2_5 for ',&
                        aeroname
                   call stop2(457)
                endif

                val= w1* rpm2_5(j1)+w2* rpm2_5(j2)+w3* rpm2_5(j3)+w4* rpm2_5(j4)+ &
                     w5* rpm2_5(j5)+w6* rpm2_5(j6)+w7* rpm2_5(j7)+w8* rpm2_5(j8)+val
                val2=w1* spm2_5(j1)+w2* spm2_5(j2)+w3* spm2_5(j3)+w4* spm2_5(j4)+ &
                     w5* spm2_5(j5)+w6* spm2_5(j6)+w7* spm2_5(j7)+w8* spm2_5(j8)+val2

                nullify(spm2_5,rpm2_5)

                aeroname='dust2'

                call gsi_bundlegetpointer(sval,trim(aeroname),spm2_5,istatus)
                if(istatus /= 0) then
                   write(6,*) 'error gsi_bundlegetpointer in stppm2_5 for ',&
                        aeroname
                   call stop2(456)
                endif

                call gsi_bundlegetpointer(rval,trim(aeroname),rpm2_5,istatus)

                if(istatus /= 0) then
                   write(6,*) 'error gsi_bundlegetpointer in stppm2_5 for ',&
                        aeroname
                   call stop2(457)
                endif

                val= w1* rpm2_5(j1)+w2* rpm2_5(j2)+w3* rpm2_5(j3)+w4* rpm2_5(j4)+ &
                     w5* rpm2_5(j5)+w6* rpm2_5(j6)+w7* rpm2_5(j7)+w8* rpm2_5(j8)+val
                val2=(w1* spm2_5(j1)+w2* spm2_5(j2)+w3* spm2_5(j3)+w4* spm2_5(j4)+ &
                     w5* spm2_5(j5)+w6* spm2_5(j6)+w7* spm2_5(j7)+w8* spm2_5(j8)) * d_2_5 +val2

                nullify(spm2_5,rpm2_5)

                do kk=1,nstep
                   qq=val2+sges(kk)*val
                   pen(kk)=qq*qq*pm2_5ptr%err2
                end do
             else
                pen(1)=pm2_5ptr%res*pm2_5ptr%res*pm2_5ptr%err2
             end if

!  modify penalty term if nonlinear qc

             if (nlnqc_iter .and. pm2_5ptr%pg > tiny_r_kind .and. &
                  pm2_5ptr%b  > tiny_r_kind) then
                pm2_5_pg=pm2_5ptr%pg*varqc_iter
                cg_pm2_5=cg_term/pm2_5ptr%b
                wnotgross= one-pm2_5_pg
                wgross = pm2_5_pg*cg_pm2_5/wnotgross
                do kk=1,max(1,nstep)
                   pen(kk)= -two*log((exp(-half*pen(kk))+wgross)/(one+wgross))
                end do
             endif

             out(1) = out(1)+pen(1)*pm2_5ptr%raterr2
             do kk=2,nstep
                out(kk) = out(kk)+(pen(kk)-pen(1))*pm2_5ptr%raterr2
             end do
          end if

          !pm2_5ptr => pm2_5ptr%llpoint
          pm2_5ptr => pm2_5Node_nextcast(pm2_5ptr)

       end do

    endif

    return
  end subroutine stppm2_5

end module stppm2_5mod


text:
Pictures, videos, audio, or texts about pornography or what some people call the good stuff.

 Where to find it in Davis
  X49er VideoX
  Aggie Liquor
  XBordersX
  Fast and Easy Mart Fast and Easy Mart (Fast and Sleazy)
  Illusions
  Newsbeat
  Quick Shop Market
  7Eleven
  Internet Service Providers The Internet
XPublicXPubic screenings
In the past, UC Davis has been host to free public screenings of pornographic movies. http://theaggie.org/article/3419 For example, in 2009, ASUCDs Entertainment Council filled 123 Sciences Lecture Hall SciLec with viewers eager to see Pirates II: Stagnettis Revenge on the biggest screen of the biggest lecture hall on campus.

 Outside of Davis
  Frys Electronics (in Sacramento)
  Goldies Adult Superstore (in Sacramento  201 N 12th St, (916) 4475860) is pretty much the biggest adult store around.
  Video Clearance Center (http://www.google.com/local?hlen&lr&c2coff1&clientfirefoxa&rlsorg.mozilla:enUS:official&qvideo+clearance+center&nearSacramento,+CA&saX&oilocalr Several locations in Sacramento) may have the widest variety of cheap porn.
  see Computer and Internet Internet Service Providers

 Pornography In The Dorms

Pornography is particularly annoying when everyone in the dorms simultaneously realizes that cybernudies are available at 10mb/s.  When spring fever hits, its possible one giddy person will hog up all the bandwidth and affect everyone else. Its truly disturbing to sneak into a roommates room and see the weird fudge hes downloading. Wild. Taking a picture of his screen with a digicam provides good laughs and possible blackmailing later on! For more nakedness, see: Sex.

 Production

Understandably, production is kept fairly quiet.  Any local noted pornographers, feel free to add yourselves.

  In 2004 http://suicidegirls.com SuicideGirls featured a photo of a young woman on Stevenson Bridge.  Her name on the site was http://suicidegirls.com/girls/Absinthe/ Absinthe.  She was a Davis native, attending Davis High School and graduating in 2004.
   This has been bothering me for a long time.  I knew the girl in question and I saw the pictures on SuicideGirls.  The pictures were not porn, they were erotica.  Theres a pretty big difference.  She was/is a good kid and should not have to bear the accusation of being a porn model for doing something a little wild in her youth.
   Welcome to Davis, I recognize Absinthe.... I told her she could do modeling in like Jr. High. I also know of several mature Davisites and another 04 alumni who became active in the porn industry. Good times Users/StevenDaubert

  Davis professor Darrell Hamamoto created a 50minute pornographic film entitled Skin on Skin and a shorter film, Yellowcaust: A Patriot Act, consisting of clips from the first. They mix intensely sexually explicit scenes with a deep message, mainly... a juxtaposition of the subjugation of Asians against graphic sex. There is a documentary about Skin on Skin titled http://www.mastersofthepillow.com Masters of the Pillow. It was featured on a segment of The Daily Show and is available to watch http://www.thedailyshow.com/watch/wednovember192003/theysohorny here.

  In 2011, genderqueer porn star http://jizlee.com Jiz Lee gave a talk on campus for GASCs Generation Sex Week. They spoke on their involvement with the queer porn movement in San Francisco, visibility for marginalized communities in porn including queer community queer women and trans people of color, ethical pornography, and porn as art and activism.

  In 2011 A former student of UC Davis appeared at the Moonlite Bunny Ranch, a legal brothel in Nevada turn reality show. Dani James appeared in a soft core porn directed and produced by the same staff as their long time running HBO series Cathouse. Dani starred with a fellow nameless UCD student in the soft core porn style light bondage scene. (http://www.bunnyranch.com/showgal.php?gcontent/BR/Bunnies/danijames/1393/2_1&s27 link)

In Peets (downtown) the other day, I saw some weird guy peddling porn...just something randomUsers/ArianeMetz

Hrmm I dont like the implication that porn is worse than erotica or that a good person wouldnt be associated with porn... (see note about SuicideGirls). Its fine if people wanna differentiate between erotica and porn, but not if were making value judgements about porn/erotica and the people involved in them.  Users/JenniferGiang


text:
!
! sum program
!
! Example of a do directive
!
! 2007,2008
! Creation: A. Muller, 2007
! Modification: F. Silber-Chaussumier

      program sum
      implicit none
      integer N
      parameter (N=10)
      integer i,a(N,2),b(N),c(N)

      do 5 i=1,N
         a(i,1)=0
         a(i,2)=0
         b(i)=0
         c(i)=-1
 5    continue
!$omp PARALLEL DO
      do 10 i=1,N
         a(i,1)=i
         a(i,2)=2*i
         b(i)=i*10
         c(i)=0
 10   continue

!$omp PARALLEL DO
      do 20 i=2,N
         c(i)=a(i-1,1)+a(i,1)+b(i)
 20   continue

!$omp PARALLEL DO
      do 30 i=1,N/2
         b(i)=i
         b(N+1-i)=i
 30   continue

      print *,a
      print *,b
      print *,c
      end


text:
!
! Copyright (C) 2001-2007 Quantum ESPRESSO group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!
!-----------------------------------------------------------------------
SUBROUTINE atomic_wfc (ik, wfcatom)
  !-----------------------------------------------------------------------
  !
  ! This routine computes the superposition of atomic wavefunctions
  ! for k-point "ik" - output in "wfcatom"
  !
  USE kinds,      ONLY : DP
  USE constants,  ONLY : tpi, fpi, pi
  USE cell_base,  ONLY : omega, tpiba
  USE ions_base,  ONLY : nat, ntyp => nsp, ityp, tau
  USE basis,      ONLY : natomwfc
  USE gvect,      ONLY : mill, eigts1, eigts2, eigts3, g
  USE klist,      ONLY : xk
  USE wvfct,      ONLY : npwx, npw, nbnd, igk
  USE us,         ONLY : tab_at, dq
  USE uspp_param, ONLY : upf
  USE spin_orb,   ONLY : lspinorb, rot_ylm, fcoef, lmaxx, domag, &
                         starting_spin_angle
  !
  implicit none
  !
  integer, parameter :: npol=1 !substitute for noncollin_module%npol
  integer, intent(in) :: ik
  complex(DP), intent(out) :: wfcatom (npwx, npol, natomwfc)
  !
  integer :: n_starting_wfc, lmax_wfc, nt, l, nb, na, m, lm, ig, iig, &
             i0, i1, i2, i3, nwfcm
  real(DP), allocatable :: qg(:), ylm (:,:), chiq (:,:,:), gk (:,:)
  complex(DP), allocatable :: sk (:), aux(:)
  complex(DP) :: kphase, lphase
  real(DP) :: arg, px, ux, vx, wx

  call start_clock ('atomic_wfc')

  ! calculate max angular momentum required in wavefunctions
  lmax_wfc = 0
  do nt = 1, ntyp
     lmax_wfc = MAX ( lmax_wfc, MAXVAL (upf(nt)%lchi(1:upf(nt)%nwfc) ) )
  enddo
  !
  nwfcm = MAXVAL ( upf(1:ntyp)%nwfc )
  !
  allocate ( ylm (npw,(lmax_wfc+1)**2), chiq(npw,nwfcm,ntyp), &
             sk(npw), gk(3,npw), qg(npw) )
  !
  do ig = 1, npw
     gk (1,ig) = xk(1, ik) + g(1, igk(ig) )
     gk (2,ig) = xk(2, ik) + g(2, igk(ig) )
     gk (3,ig) = xk(3, ik) + g(3, igk(ig) )
     qg(ig) = gk(1, ig)**2 +  gk(2, ig)**2 + gk(3, ig)**2
  enddo
  !
  !  ylm = spherical harmonics
  !
  call ylmr2 ((lmax_wfc+1)**2, npw, gk, qg, ylm)
  !
  ! set now q=|k+G| in atomic units
  !
  do ig = 1, npw
     qg(ig) = sqrt(qg(ig))*tpiba
  enddo
  !
  n_starting_wfc = 0
  !
  ! chiq = radial fourier transform of atomic orbitals chi
  !
  do nt = 1, ntyp
     do nb = 1, upf(nt)%nwfc
        if ( upf(nt)%oc (nb) >= 0.d0) then
           do ig = 1, npw
              px = qg (ig) / dq - int (qg (ig) / dq)
              ux = 1.d0 - px
              vx = 2.d0 - px
              wx = 3.d0 - px
              i0 = INT( qg (ig) / dq ) + 1
              i1 = i0 + 1
              i2 = i0 + 2
              i3 = i0 + 3
              chiq (ig, nb, nt) = &
                     tab_at (i0, nb, nt) * ux * vx * wx / 6.d0 + &
                     tab_at (i1, nb, nt) * px * vx * wx / 2.d0 - &
                     tab_at (i2, nb, nt) * px * ux * wx / 2.d0 + &
                     tab_at (i3, nb, nt) * px * ux * vx / 6.d0
           enddo
        endif
     enddo
  enddo

  deallocate (qg, gk)
  allocate ( aux(npw) )
  !
  wfcatom(:,:,:) = (0.0_dp, 0.0_dp)
  !
  do na = 1, nat
     arg = (xk(1,ik)*tau(1,na) + xk(2,ik)*tau(2,na) + xk(3,ik)*tau(3,na)) * tpi
     kphase = CMPLX(cos (arg), - sin (arg) ,kind=DP)
     !
     !     sk is the structure factor
     !
     do ig = 1, npw
        iig = igk (ig)
        sk (ig) = kphase * eigts1 (mill (1,iig), na) * &
                           eigts2 (mill (2,iig), na) * &
                           eigts3 (mill (3,iig), na)
     enddo
     !
     nt = ityp (na)
     do nb = 1, upf(nt)%nwfc
        if (upf(nt)%oc(nb) >= 0.d0) then
           l = upf(nt)%lchi(nb)
           lphase = (0.d0,1.d0)**l
           !
           !  the factor i^l MUST BE PRESENT in order to produce
           !  wavefunctions for k=0 that are real in real space
           !
              !
              call atomic_wfc___ ( )
              !
           !
        END IF
        !
     END DO
     !
  END DO

  if (n_starting_wfc /= natomwfc) call errore ('atomic_wfc', &
       'internal error: some wfcs were lost ', 1)

  deallocate(aux, sk, chiq, ylm)

  call stop_clock ('atomic_wfc')
  return

CONTAINS

   SUBROUTINE atomic_wfc___( )
   !
   ! ... LSDA or nonmagnetic case
   !
   DO m = 1, 2 * l + 1
      lm = l**2 + m
      n_starting_wfc = n_starting_wfc + 1
      if (n_starting_wfc > natomwfc) call errore &
         ('atomic_wfc___', 'internal error: too many wfcs', 1)
      !
      DO ig = 1, npw
         wfcatom (ig, 1, n_starting_wfc) = lphase * &
            sk (ig) * ylm (ig, lm) * chiq (ig, nb, nt)
      ENDDO
      !
   END DO
   !
   END SUBROUTINE atomic_wfc___
   !
END SUBROUTINE atomic_wfc

max_stars_repo_name: utdsimmons/ohpc

max_issues_repo_name: utdsimmons/ohpc

max_forks_repo_name: utdsimmons/ohpc