text:
      subroutine carbon
      !! This code simulates organic C, N, and P cycling in the soil
      !! The code was developed by Armen R. Kemanian and Stefan Julich
      !! It has been adapted from Kemanian and Stockle (2010) (European Journal of Agronomy 32:22-29)
      !! and crafted to accomodate to SWAT conventions
      !! Plant residues and manure residues are decomposed separately
      !! For convenience, the denitrification subroutine is called from here
      !! March 2009: testing has been minimal and further adjustments are expected
      !! manuscript describing this subroutine to be submitted to Ecological Modelling (September, 2010)
      !! use with caution and report anomalous results to akemanian@psu.edu, and jeff.arnold@ars.usda.edu, stefan.julich@tudor.lu

!!!!!!!
!      cmup_kgh       kg/ha    current soil carbon for first soil layer
!      cmtot_kgh      kg/ha    current soil carbon integrated - aggregating
!                                 all soil layers

      !! local variables
      !! cx = saturated soil carbon concentration (%) (Hassink and Whitmore, 1997)
      !! decf = decomposition factor
      !! net_N = nitrogen net mineralization
      !! net_P = phosphorus net mineralization
      !! rnet_N, rnet_P, mnet_N, mnet_P for residue and manure respectively
      !! sol_cdec = carbon decomposition (CO2)
      !! tilf = tillage factor
      !! resc_hum = humified carbon residue
      !! manc_hum = humified carbon manure

      use parm

      INTERFACE

	real*8 Function fwf(fc,wc,pwp)
      real*8, intent (in) :: fc, wc, pwp
      end function

      real*8 Function fof(void,por)
        real*8, intent (in) :: void, por
      end function

      real*8 Function ftilf(tillage, wc, sat)
        real*8, intent (in out) :: tillage
        real*8, intent (in) :: wc, sat
      End function

      real*8 Function fcx(pclay)
        real*8, intent (in) :: pclay
      End function

	real*8 Function fsol_cdec(pcarbon, cx, cfdec, tilf, csf, sol_cmass)
        real*8, intent (in) :: pcarbon, cx, cfdec, tilf, csf, sol_cmass
	End function

      real*8 Function fCNnew(yy1,yy2,CNpool,yy5)
        real*8, intent(in) :: yy1, yy2, CNpool
        real*8, intent(in) :: yy5
      End function

	real*8 Function fhc(pclay, pcarbon, cx)
        real*8, intent(in) :: pclay, pcarbon, cx
	End function

	real*8 Function fnetmin(poold, R1, R2, hc, dummy, poolm, xinorg, cc1)
        real*8, intent(in) :: R1, R2, hc, poolm, xinorg, cc1
        real*8, intent(in out) :: poold, dummy
      End function

      END INTERFACE

    !! private variables
      real*8 :: cx, decf, rhc, mhc, sol_cdec, tilf
      real*8 :: resc_hum, manc_hum
      real*8 :: xx, xx1, xx2, xx3, xx4, csf
      real*8 :: rdc, mdc, wdn, cdg, sut
      real*8 :: CNsoil, CPsoil, NPsoil
      real*8 :: CNres, CPres, CNman, CPman, rCNnew, mCNnew
      real*8 :: sol_thick, sol_mass, sol_cmass, sol_nmass
      real*8 :: net_N, net_P, rnet_N, rnet_P, mnet_N, mnet_P
      real*8 :: wc, fc, wf, of, void
      real*8 :: sat, ffres1, ffres2, ffman1, ffman2

    !! mass balance variables
      real*8 :: sum_c_i, sum_n_i, sum_p_i
      real*8 :: sum_c_f, sum_n_f, sum_p_f
      real*8 :: bal_c, bal_n, bal_p

    !! output variables for soil profile
      real*8 :: cmass_pro, nmass_pro, sol_orgp_pro
      real*8 :: sol_rsd_pro, sol_fon_pro, sol_fop_pro
      real*8 :: sol_mc_pro, sol_mn_pro, sol_mp_pro
      real*8 :: sol_no3_pro, sol_solp_pro, sol_nh3_pro
      real*8 :: sol_cdec_pro, wdn_pro, net_N_pro, solN_net_min
      real*8 :: solN_net_min_pro

      integer :: j, k, kk

      j = 0; wdn = 0
      j = ihru

    !! initialize
      cmass_pro = 0.
      nmass_pro = 0.
      sol_orgp_pro = 0.
      sol_rsd_pro = 0.
      sol_fon_pro = 0.
      sol_fop_pro = 0.
      sol_mc_pro = 0.
      sol_mn_pro = 0.
      sol_mp_pro = 0.
      sol_no3_pro = 0.
      sol_nh3_pro = 0.
      sol_solp_pro = 0.
      sol_cdec_pro = 0.
      wdn_pro = 0.
      net_N_pro = 0.
      solN_net_min = 0.
      solN_net_min_pro=0.

!!    zero new carbon variables for output.hru
      cmup_kgh(j) = 0.
      cmtot_kgh(j) = 0.
!!    zero new carbon variables for output.hru

      if (sol_cbn(1,j) == 0.) return

      do k = 1, sol_nly(j)

      sum_c_i = 0.
      sum_n_i = 0.
      sum_p_i = 0.

      rdc = 0.
      mdc = 0.
      rhc  = 0.
      mhc = 0.
      sol_cdec = 0.
      resc_hum = 0.
      manc_hum = 0.
      rnet_N = 0.
      rnet_P = 0.
      mnet_N = 0.
      mnet_P = 0.
      net_N = 0.
      net_P = 0.

      rCNnew = 0.
      mCNnew = 0.

      wc = 0.
      fc = 0.
      wf = 0.
      of = 0.
      void = 0.
      cdg = 0.
      sut = 0.
      csf = 0.

      ffres = 0.
      ffres1 = 0.
      ffres2 = 0.

      ffman = 0.
      ffman1 = 0.
      ffman2 = 0.

	if (k == 1) then
	    sol_thick = sol_z(k,j)
	else
	    sol_thick = sol_z(k,j) - sol_z(k-1,j)
	end if

      !! soil carbon and nitrogen mass
      sol_mass = (sol_thick / 1000.) * 10000. * sol_bd(k,j)
     &     * 1000. * (1- sol_rock(k,j) / 100.)
      sol_cmass = sol_mass * (sol_cbn(k,j) / 100.)
      sol_nmass = sol_mass * (sol_n(k,j) / 100.)

	!! sum initial C,N,P for mass balance
	sum_c_i = sol_cmass + 0.43 * sol_rsd (k,j) + sol_mc(k,j)
	sum_n_i = sol_nmass + sol_no3(k,j) + sol_nh3(k,j) + sol_fon(k,j)
     &	+ sol_mn(k,j)
	sum_p_i = sol_orgp(k,j) + sol_solp(k,j) + sol_fop(k,j)
     &	+ sol_mp(k,j)

      kk = k
      if (k == 1) kk = 2

	if (sol_tmp(kk,j) > 0. .and. sol_st(kk,j) > 0.) then
		!! microbial processes if temp > 0 C
		!! microbial processes if water > pwp

		!!compute soil water factor - sut
		fc = sol_fc(kk,j) + sol_wpmm(kk,j)  ! units mm
		wc = sol_st(kk,j) + sol_wpmm(kk,j)  ! units mm
		sat = sol_ul(kk,j) + sol_wpmm(kk,j) ! units mm
		void = sol_por(kk,j) * (1. - wc / sat)

        wf = fwf(fc,wc,sol_wpmm(kk,j))
        of = fof(void,sol_por(kk,j))
        sut = wf * of

        !! compute soil temperature factor - cdg
        cdg = fcgd(sol_tmp(kk,j))

		!! compute combined factor
		xx = 0.
!!		xx = sqrt(cdg * sut)
		xx = (cdg * sut) ** cf(j)
		if (xx < 0.) xx = 0.
		if (xx > 1.) xx = 1.
		csf = xx

		!! call denitrification (to use void and cdg factor)
		wdn = 0.
		if (cdg > 0. .and. void <= 0.1) then
			call ndenit(k,j,cdg,wdn,void)
		end if
		wshd_dnit = wshd_dnit + wdn * hru_dafr(j)
		wdntl = wdntl + wdn

		!! calculate soil carbon 'decomposition'
		!! tillage factor
		tilf = ftilf(tillagef(k,j), wc, sat)
		!! saturated carbon concentration (%)
		cx = fcx(sol_clay(k,j))
		!! soil carbon decomposition
		sol_cdec=fsol_cdec(sol_cbn(k,j),cx,cfdec(j),
     & 		tilf,csf,sol_cmass)

        !! residue and manure decomposition and N and P mineralization
        CNsoil = sol_cbn(k,j) / sol_n(k,j)
        if (sol_orgp(k,j) < .01) then
          NPsoil = 25.
        else
          NPsoil = (sol_mass * sol_n(k,j) / 100.)/ sol_orgp(k,j)
        end if
        sol_orgp(k,j) = dmin1(sol_orgp(k,j), 25.)
        CPsoil = CNsoil * NPsoil

        if (sol_rsd(k,j) > 0.00001) then

            ! This IF STATEMENT preserved in case residues with different properties are considered
            !if (idplt(nro(j),icr(j),j) > 0) then
            !    decr = rsdco_pl(idplt(nro(j),icr(j),j)) * csf
            !else
            !    decr = 0.05 *  csf
            !end if

            !! residue decomposition
            rdc = 0.05 * csf * sol_rsd(k,j)

			!! humification factor
			rhc = fhc(sol_clay(k,j),sol_cbn(k,j),cx) * cfh(j)

            if (sol_fon(k,j) < .01) then
              CNres = 15.
            else
              CNres = 0.43 * sol_rsd(k,j) / sol_fon(k,j)
            end if
            CNres = dmax1(CNres, 15.)

            if (sol_fop(k,j) < .01) then
              CPres = 400.
            else
              CPres = 0.43 * sol_rsd(k,j) / sol_fop(k,j)
            end if
            CPres = dmax1(CPres, 400.)

            !! CN of new organic matter (humified residue)
            rCNnew = fCNnew(sol_no3(k,j),sol_mass,CNres, 1.1D+02)   ! 110.)

			!! N mineralization / immobilization
			xx1 = sol_rsd(k,j)
			xx2 = sol_no3(k,j) + sol_nh3(k,j)
			rnet_N = fnetmin(rdc,CNres,rCNnew,rhc,ffres1,xx1,xx2, 4.3D-01)

            !! P mineralization / immobilization
            xx2 = sol_solp(k,j)
            rnet_P = fnetmin(rdc,CPres,CPsoil,rhc,ffres2,xx1,xx2, 4.3D-01)

            !! check if P limits N mineralization and re-adjust
            if (ffres2 < ffres1) then
                  rnet_N = rdc * 0.43 * (1. / CNres - rhc / rCNnew)
                ffres = ffres2
            else
                ffres = ffres1
            end if
        end if

        ! manure decomposition and N and P mineralization
        if (sol_mc(k,j) > 0.00001 .and. sol_mp(k,j) > 0.00001) then

            !! residue decomposition
            !! decomposition rate about 1/2 of residue
            mdc = 0.025 * csf * sol_mc(k,j)

            !! humification factor
			mhc = 1.6 * fhc(sol_clay(k,j),sol_cbn(k,j),cx)

			CNman = sol_mc(k,j) / sol_mn(k,j)
			CPman = sol_mc(k,j) / sol_mp(k,j)

			!! CN of new organic matter (humified manure)
			mCNnew = fCNnew(sol_no3(k,j),sol_mass,CNman, 5.5D+01)

			!! N mineralization / immobilization
			xx1 = sol_mc(k,j)
			xx2 = sol_no3(k,j) + sol_nh3(k,j)
			mnet_N = fnetmin(mdc,CNman,mCNnew,mhc,ffman1,xx1,xx2,1.0D+00)

            !! P mineralization / immobilization
            xx2 = sol_solp(k,j)
            mnet_P = fnetmin(mdc,CPman,CPsoil,mhc,ffman2,xx1,xx2,1.0D+00)

            !! check if P or N limit mineralization and re-adjust
            if (ffman2 < ffman1) then
                  mnet_N = mdc * (1. / CNman - mhc / mCNnew)
                ffman = ffman1
            else
                ffman = ffman2
            end if
        end if

        !! check if sufficient mineral N for both residue and manure decomposition
		if ((sol_no3(k,j) + sol_nh3(k,j)) > 0.) then
			xx = (rnet_N + mnet_N) / (sol_no3(k,j) + sol_nh3(k,j))
            if (xx < -1.) then
                rdc = -rdc / xx
                rnet_N = -rnet_N / xx
                rnet_P = -rnet_P / xx
                ffres = -ffres / xx

                mdc = -mdc / xx
                mnet_N = -mnet_N / xx
                mnet_P = -mnet_P / xx
                ffman = -ffman / xx
            end if
        end if

        !! check if sufficient mineral P for both residue and manure decomposition
        if (sol_solp(k,j) > 0.) then
            xx = (rnet_P + mnet_P) / sol_solp(k,j)
            if (xx < -1.) then
                rdc = -rdc / xx
                rnet_N = -rnet_N / xx
                rnet_P = -rnet_P / xx
                ffres = -ffres / xx

                mdc = -mdc / xx
                mnet_N = -mnet_N / xx
                mnet_P = -mnet_P / xx
                ffman = -ffman / xx
            end if
        end if

        resc_hum = rhc * rdc * 0.43
        manc_hum = mhc * mdc
        net_N = rnet_N + mnet_N
        net_P = rnet_P + mnet_P
        if (resc_hum == 0.) rCNnew = 1000.
        if (manc_hum == 0.) mCNnew = 1000.

		!! C N P pools update
		sol_cmass = sol_cmass + resc_hum + manc_hum - sol_cdec
 		sol_cbn(k,j) = 100. * sol_cmass / sol_mass

		sol_nmass = sol_nmass - sol_cdec / CNsoil
     &       	  + resc_hum / rCNnew + manc_hum / mCNnew

!! april 2010 output net_min
		solN_net_min=-(- sol_cdec / CNsoil
     &             + resc_hum / rCNnew + manc_hum / mCNnew)
 		sol_n(k,j) = 100. * sol_nmass / sol_mass
		sol_orgn(k,j) = sol_nmass ! for output
        sol_orgp(k,j) = sol_orgp(k,j) - sol_cdec / CPsoil
     &         + (resc_hum + manc_hum) / CPsoil

                if (ffres > 1.) ffres = 1.

        sol_rsd(k,j) = sol_rsd(k,j) * (1. - ffres)
        sol_fon(k,j) = sol_fon(k,j) * (1. - ffres)
        sol_fop(k,j) = sol_fop(k,j) * (1. - ffres)

                if (ffman > 1.) ffman = 1.

        sol_mc(k,j) = sol_mc(k,j) * (1. - ffman)
        sol_mn(k,j) = sol_mn(k,j) * (1. - ffman)
        sol_mp(k,j) = sol_mp(k,j) * (1. - ffman)

		!sol_no3(k,j) = sol_no3(k,j) + net_N +
     &	!	sol_cdec * (1. / CNsoil)

		! add positive n-mineralization to ammonia pool in the layer
		if (rnet_N>0.) sol_nh3(k,j) = sol_nh3(k,j) + rnet_N
		if (mnet_N>0.) sol_nh3(k,j) = sol_nh3(k,j) + mnet_N

		if (rnet_N<0.) then
			if (abs(rnet_N) < sol_nh3(k,j)) then
				sol_nh3(k,j) = sol_nh3(k,j) + rnet_N
			else
				xx4 = sol_nh3(k,j) + rnet_N
				sol_nh3(k,j) = 0.
				sol_no3(k,j) = sol_no3(k,j) + xx4
			end if
		end if

		if (mnet_N<0.) then
			if (abs(mnet_N) < sol_nh3(k,j)) then
				sol_nh3(k,j) = sol_nh3(k,j) + mnet_N
			else
				xx4 = sol_nh3(k,j) + mnet_N
				sol_nh3(k,j) = 0.
				sol_no3(k,j) = sol_no3(k,j) + xx4
			end if
		end if

		sol_nh3(k,j) = sol_nh3(k,j) + sol_cdec * (1. / CNsoil)

		sol_solp(k,j) = sol_solp(k,j) + net_P +
     &		sol_cdec * (1. / CPsoil)

          wshd_rmn = wshd_rmn + net_N * hru_dafr(j)
		wshd_rmp = wshd_rmp + net_P * hru_dafr(j)

		If (sol_rsd(k,j) < 1e-10) sol_rsd(k,j) = 1e-10
		If (sol_fon(k,j) < 1e-11) sol_fon(k,j) = 1e-11
		If (sol_fop(k,j) < 1e-12) sol_fop(k,j) = 1e-12
		If (sol_mc(k,j) < 1e-10) sol_mc(k,j) = 1e-10
		If (sol_mn(k,j) < 1e-11) sol_mn(k,j) = 1e-11
		If (sol_mp(k,j) < 1e-12) sol_mp(k,j) = 1e-12
		If (sol_no3(k,j) < 1e-12) sol_no3(k,j) = 1e-12

	end if

	!! balance file cswat_balance
	sum_c_f = sol_cmass + 0.43 * sol_rsd(k,j) + sol_cdec
     &	+ (1. - rhc) * rdc * 0.43 + sol_mc(k,j) + (1. - mhc) * mdc
	sum_n_f = sol_nmass + sol_no3(k,j) + sol_nh3(k,j) + sol_fon(k,j)
     &	+ sol_mn(k,j) + wdn
	sum_p_f = sol_orgp(k,j) + sol_solp(k,j) + sol_fop(k,j)
     &	+ sol_mp(k,j)

	bal_c = sum_c_i - sum_c_f
	bal_n = sum_n_i - sum_n_f
	bal_p = sum_p_i - sum_p_f

	!! writing daily output by layer for testing purposes of the routine SJ and AK 2010
	!!if (i==365) then
!!	   write (98,9000) iyr, i, k, j, sol_cmass, sol_cbn(k,j),
!!     &	sol_nmass, sol_n(k,j), sol_orgp(k,j), sol_rsd(k,j),
!!     &	sol_fon(k,j), sol_fop(k,j), sol_solp(k,j), sol_mc(k,j),
!!     &	sol_mn(k,j), sol_mp(k,j), sol_no3(k,j),
!!     &	sol_cmass/sol_nmass, sol_nmass/sol_orgp(k,j), sol_nh3(k,j),
!!     &    tilf, sol_cdec, wdn, net_N
!!      	 write (99,9001) iyr, i, k, j, bal_c, sum_c_i, sum_c_f,
!!     &     bal_n,sum_n_i, sum_n_f, bal_p, sum_p_i, sum_p_f
	!!end if

!!    carbon outputs for .hru file
      if (k == 1) cmup_kgh(j) = sol_cmass
      cmtot_kgh(j) = cmtot_kgh(j) + sol_cmass
!!    carbon outputs for .hru file

	cmass_pro = cmass_pro + sol_cmass
	nmass_pro = nmass_pro + sol_nmass
	sol_rsd_pro = sol_rsd_pro + sol_rsd(k,j)
	sol_cdec_pro = sol_cdec_pro + sol_cdec
	sol_fon_pro = sol_fon_pro + sol_fon(k,j)
	sol_no3_pro = sol_no3_pro + sol_no3(k,j)
	sol_nh3_pro = sol_nh3_pro + sol_nh3(k,j)
	sol_orgp_pro = sol_orgp_pro + sol_orgp(k,j)
	sol_fop_pro = sol_fop_pro + sol_fop(k,j)
	sol_solp_pro = sol_solp_pro + sol_solp(k,j)
	sol_mc_pro = sol_mc_pro + sol_mc(k,j)
	sol_mn_pro = sol_mn_pro + sol_mn(k,j)
	sol_mp_pro = sol_mp_pro + sol_mp(k,j)
	wdn_pro = wdn_pro + wdn
	net_N_pro = net_N_pro + net_N
	solN_net_min_pro = solN_net_min_pro + solN_net_min

      end do

	!! writing daily profile output

        write (100,9002) iyr, i, j, cmass_pro, sol_rsd_pro, sol_mc_pro

!9000  format(i4,';',i3,';',i1,';',i4,20(';',f10.3))
!9001  format(i4,';',i3,';',i1,';',i4,10(';',f10.3))
9002  format(i4,';',i3,';',i4,';',f11.3,';',f11.3,';',f11.3)

      return
      end subroutine

	!! LOCAL FUNCTIONS
	real*8 Function fwf(fc,wc,pwp)
          real*8, intent (in) :: fc, wc, pwp
		xx2 = 0.
		if (wc <= pwp) then
			xx2 = 0.4 * wc / pwp
		else if (wc <= fc) then
			xx2 = 0.4 + 0.6 * (wc - pwp)/(fc - pwp)
		else
			xx2 = 1.
		end if

        !!fwf = (1. + (1. - xx2) / (1. - 0.75)) * (xx2 / 1.) ** (1./ (1. - 0.75))
        fwf = (1. + (1. - xx2) / 0.25) * (xx2) ** 4.
      End function

      real*8 Function fof(void,por)
        real*8, intent (in) :: void, por
        xx3 = 0.
        if (void >= 0.1) then
            xx3 = 0.2 + 0.8 * (void - 0.1) / (por - 0.1)
        else
            xx3 = 0.2 * void / 0.1
        end if
        fof = 0.5 + 0.5 * xx3 / (xx3 + Exp(-20. * xx3))
      End function

	real*8 Function fcgd(xx)
        real*8, intent (in) :: xx
		tn = -5.
	  top = 35.
		tx = 50.
		qq = (tn - top)/(top - tx)
	  fcgd = ((xx-tn)**qq)*(tx-xx)/(((top-tn)**qq)*(tx-top))
		if (fcgd < 0.) fcgd = 0.
	End function

      real*8 Function ftilf(tillage, wc, sat)
        real*8, intent (in out) :: tillage
        real*8, intent (in) :: wc, sat
        !! tillage factor effect on decomposition
        !! tillage factor returns to baseline (=1) based on WC
        tillage = tillage * (1. - 0.02 * wc/sat)
        if (tillage < 0.) tillage = 0.
        ftilf = 1. + tillage
      End function

      real*8 Function fcx(pclay)
        real*8, intent (in) :: pclay
        !! saturated soil carbon concentration (%) from Hassink and Whitmore 1997
        fcx = 2.11 + 0.0375 * pclay
      End function

	real*8 Function fsol_cdec(pcarbon, cx, cfdec, tilf, csf, sol_cmass)
        real*8, intent (in) :: pcarbon, cx, cfdec, tilf, csf, sol_cmass
		!! decomposition adjustment by current SOC
		decf = (pcarbon / cx) ** 0.5
		! if (decf > 1.) decf = 1.
		!! maximum soil carbon decomposition = 0.045 kg C / kg C per year
		fsol_cdec = cfdec / 365. * decf * tilf * csf * sol_cmass
	End function

      real*8 Function fCNnew(yy1,yy2,CNpool,yy5)
        real*8, intent(in) :: yy1,yy2,CNpool
        real*8, intent(in) :: yy5
      !! CN ratio of newly formed organic matter
      !! based on CN or decomposing residue and nitrate in soil
      !! the same approach used for crop residues and manure
      !! CNpool = the CN of decomposing pool
      !! yy1 = the layer nitrate mass
      !! yy2 = the layer soil mass
      !! yy3 = nitrate concentration g/g
      !! yy4 = correction factor based on CN of decomposing pool
      !! yy5 = input-dependent constant to correct CN ratio

        yy3 = yy1 / yy2
        yy4 = 5.5 * (1. - 1. / (1. + (CNpool / yy5)**3.))
        fCNnew = 8.5 + yy4 * (0.5 + 0.5 / (1. + (yy3 / 8e-6)**3.))
      End function

	real*8 Function fhc(pclay, pcarbon, cx)
        real*8, intent(in) :: pclay, pcarbon, cx
	!! maximum and actual humification factor
	!! hx = maximum humification factor
	!! hf = humification adjustment factor
	!! pclay = %clay
	!! pcarbon = %carbon
	!! cx = saturated soil carbon, %

	real*8 :: hx, hf

		hx = 0.09 + 0.09 * (1. - Exp(-5.5 * pclay / 100.))
		!! humification adjustment by current SOC
		if (pcarbon > cx) then
			hf = 0.
		else
			hf = 1. - (pcarbon / cx) ** 6.
		end if
		fhc = hx * hf
	End function

	real*8 Function fnetmin(poold, R1, R2, hc, dummy, poolm, xinorg, cc1)
	!! This function computes net mineralization
	!! R1 = CN or CP ratio of decomposing pool
	!! R2 = CN or CP ratio of receiving pool
	!! hc = humification rate
	!! dummy = fraction of pool being decomposed
	!! poolm = current mass of pool being decomposed
	!! poold = mass of pool decomposed
	!! xinorg = mass of NO3 or P in solution
	!! xx = net mineralization of N or P
	!! cc1 = pool's carbon fraction
        real*8, intent(in) :: R1, R2, hc, poolm, xinorg, cc1
        real*8, intent(in out) :: poold, dummy

		xx = 0.
		xx = poold * cc1 * (1. / R1 - hc / R2)

		if (xx > 0.) then
			dummy = poold / poolm
		else if (abs(xx)< xinorg) then
			!! net mineralization is positive or
			!! immobilization not consuming all mineral N or P
			dummy = poold / poolm
		else
			!! immobilization, consuming all mineral N or P
			xx = -xinorg
			poold = xx / cc1 * 1. / (1. / R1 - hc / R2)
			dummy = poold / poolm
    		end if

        fnetmin = xx
      End function

max_stars_repo_name: GISWAT/erosion-sediment

max_issues_repo_name: GISWAT/erosion-sediment

max_forks_repo_name: GISWAT/erosion-sediment


! { dg-do compile }
! { dg-options "-fcoarray=lib -fdump-tree-original" }
!
! Allocate/deallocate with libcaf.
!

 subroutine test()
 integer(4), allocatable :: xx[:], yy(:)[:]
 integer :: stat
 character(len=200) :: errmsg
 allocate(xx[*], stat=stat, errmsg=errmsg)
 allocate(yy(2)[*], stat=stat, errmsg=errmsg)
 deallocate(xx,yy,stat=stat, errmsg=errmsg)
 end

! { dg-final { scan-tree-dump-times "_gfortran_caf_register \\(4, 1, &xx.token, \\(void \\*\\) &xx, &stat.., &errmsg, 200\\);" 1 "original" } }
! { dg-final { scan-tree-dump-times "_gfortran_caf_register \\(8, 1, &yy.token, \\(void \\*\\) &yy, &stat.., &errmsg, 200\\);" 1 "original" } }
! { dg-final { scan-tree-dump-times "_gfortran_caf_deregister .&xx.token, 0, &stat.., &errmsg, 200.;" 1 "original" } }
! { dg-final { scan-tree-dump-times "_gfortran_caf_deregister .&yy.token, 0, &stat.., &errmsg, 200.;" 1 "original" } }
! { dg-final { scan-tree-dump-times "_gfortran_caf_deregister .&yy.token, 0, 0B, 0B, 0.;" 1 "original" } }
! { dg-final { scan-tree-dump-times "_gfortran_caf_deregister .&xx.token, 0, 0B, 0B, 0.;" 1 "original" } }


text:
Describe Users/NikZ here.
.

20091025 14:33:56 nbsp Welcome to the Wiki. Are you Andre Lee, or are you just posting to his public figure page in the first person? Users/JasonAller

20091025 14:38:58 nbsp It generally pays off to take the time to add links to text when you are pasting it into the wiki. For instance wrap the word sustainability in square brackets and quotes like {{{sustainability}}} and it will render on the page as sustainability. Users/JasonAller

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
PROGRAM pg123
!-------------------------------------------------------------------------
!      program p12.3 three dimensional analysis of Laplace's equation
!      using 8-node bricks, DEC, preconditioned conjugate gradient solver
!      diagonal preconditioner; parallel; externally generated model
!-------------------------------------------------------------------------
!USE mpi_wrapper  !remove comment for serial compilation
 USE precision; USE global_variables; USE mp_interface; USE input
 USE output; USE loading; USE timing; USE maths; USE gather_scatter
 USE new_library; USE common_mod; USE geometry_mod; IMPLICIT NONE
!neq,ntot  are now global variables - not declared
 INTEGER, PARAMETER::ndim=3,nodof=1
 INTEGER::nod,nn,nr,nip,i,j,k,m,iters,limit,iel,partitioner,meshgen,     &
   node_end,node_start,nodes_pp,loaded_freedoms,fixed_freedoms,          &
   fixed_freedoms_pp,fixed_freedoms_start,nlen,nres,is,it,               &
   loaded_freedoms_pp,loaded_freedoms_start,nels,ndof,npes_pp,n_smplx,ne
 REAL(iwp),PARAMETER::zero=0.0_iwp,penalty=1.e20_iwp
 REAL(iwp)::kx,ky,kz,ka,det,tol,up,alpha,beta,q; CHARACTER(LEN=6)::ch
 CHARACTER(LEN=15)::element; CHARACTER(LEN=50)::argv
 LOGICAL::converged=.false.
 REAL(iwp),ALLOCATABLE::points(:,:),weights(:),eld_pp(:,:),kay(:),       &
   fun(:),jac(:,:),der(:,:),deriv(:,:),col(:,:),row(:,:),                &
   kcx(:,:),kcy(:,:),kcz(:,:),diag_precon_pp(:),p_pp(:),r_pp(:),x_pp(:), &
   xnew_pp(:),u_pp(:),pmul_pp(:,:),utemp_pp(:,:),d_pp(:),val(:,:),       &
   diag_precon_tmp(:,:),store_pp(:),storkc_pp(:,:,:),eld(:),timest(:),   &
   val_f(:),g_coord_pp(:,:,:),ptl_pp(:)
 INTEGER,ALLOCATABLE::rest(:,:),g_num_pp(:,:),g_g_pp(:,:),no(:),         &
   no_pp(:),no_f_pp(:),no_pp_temp(:),sense(:),node(:),smplx(:,:)
!--------------------------input and initialisation-----------------------
 ALLOCATE(timest(20)); timest=zero; timest(1)=elap_time()
 CALL find_pe_procs(numpe,npes); CALL getname(argv,nlen)
 CALL read_p123(argv,numpe,element,fixed_freedoms,kx,ky,kz,limit,        &
   loaded_freedoms,meshgen,nels,nip,nn,nod,nr,nres,partitioner,tol)
 CALL calc_nels_pp(argv,nels,npes,numpe,partitioner,nels_pp)
 ndof=nod*nodof; ntot=ndof
 ALLOCATE(g_num_pp(nod,nels_pp),g_coord_pp(nod,ndim,nels_pp))
 g_num_pp=0; g_coord_pp=zero
 CALL read_g_num_pp(argv,iel_start,nn,npes,numpe,g_num_pp)
 IF(meshgen == 2) CALL abaqus2sg(element,g_num_pp)
 CALL read_g_coord_pp(argv,g_num_pp,nn,npes,numpe,g_coord_pp)
 IF (nr>0) THEN; ALLOCATE(rest(nr,nodof+1)); rest=0
 CALL read_rest(argv,numpe,rest); END IF
 ALLOCATE (points(nip,ndim),jac(ndim,ndim),storkc_pp(ntot,ntot,nels_pp), &
   deriv(ndim,nod),kcx(ntot,ntot),weights(nip),der(ndim,nod),kay(ndim),  &
   pmul_pp(ntot,nels_pp),utemp_pp(ntot,nels_pp),col(ntot,1),eld(ntot),   &
   g_g_pp(ntot,nels_pp),kcy(ntot,ntot),row(1,ntot),kcz(ntot,ntot))
!----------  find the steering array and equations per process -----------
 timest(2)=elap_time(); g_g_pp=0; neq=0
 IF(nr>0) THEN; CALL rearrange_2(rest)
   elements_0: DO iel=1, nels_pp
     CALL find_g4(g_num_pp(:,iel),g_g_pp(:,iel),rest)
   END DO elements_0
 ELSE
   g_g_pp=g_num_pp  !When nr = 0, g_num_pp and g_g_pp are identical
 END IF
 neq=MAXVAL(g_g_pp); neq=max_p(neq); CALL calc_neq_pp
 CALL calc_npes_pp(npes,npes_pp); CALL make_ggl(npes_pp,npes,g_g_pp)
 DO i=1,neq_pp;IF(nres==ieq_start+i-1)THEN;it=numpe;is=i;END IF;END DO
 IF(numpe==it)THEN
   OPEN(11,FILE=argv(1:nlen)//'.res',STATUS='REPLACE',ACTION='WRITE')
   WRITE(11,'(A,I5,A)')"This job ran on ", npes,"  processes"
   WRITE(11,'(A,3(I7,A))')"There are ",nn," nodes",nr,                   &
     " restrained and   ",neq," equations"
   WRITE(11,'(A,F10.4)')"Time after setup is ",elap_time()-timest(1)
 END IF
 ALLOCATE(p_pp(neq_pp),r_pp(neq_pp),x_pp(neq_pp),xnew_pp(neq_pp),        &
   u_pp(neq_pp),diag_precon_pp(neq_pp),d_pp(neq_pp))
 r_pp=zero; p_pp=zero; x_pp=zero; xnew_pp=zero; diag_precon_pp=zero
!-------------- element stiffness using DEC and storage ----------------
 storkc_pp=zero;  dim_cmplx = 3;  dim_embbd = 3;  k=dim_cmplx+1
 ALLOCATE(num_pelm_pp(1)); num_pelm_pp(1)=nod;  extra_pelm=0;  glb_offset=0
 selectcase(element)
 case('hexahedron');   ne = 6;  glb_num_elm = (/ nod, 19, 18, ne /);
 case('tetrahedron');  ne = 1;  glb_num_elm = (/ nod, 19, 18, ne /);
 end SELECT
 num_elm=glb_num_elm;  ALLOCATE(lcl_complex(k));
 ALLOCATE(&
   lcl_complex(1)%centers(nod,dim_embbd),&
   lcl_complex(k)%orientation(ne),&
   lcl_complex(k)%node_indx(ne,k))
 lcl_complex(k)%orientation=0;  indx_offset = 0

 kay(1)=kx; kay(2)=ky; kay(3)=kz;
 elements_1: DO iel=1,nels_pp
   !- split elements into simplices
   CALL elm2smplx(num_elm(dim_cmplx+1),lcl_complex(dim_cmplx+1)%node_indx,&
      g_coord_pp(:,:,iel),element,nod)

   !- Recursively compute element (co-)boundaries
   DO k=dim_cmplx+1,2,-1;  CALL calc_bndry_cobndry(k);  END DO

   !- setup connectivity
   DO k=1,dim_cmplx+1
     lcl_complex(k)%lcl_node_indx = lcl_complex(k)%node_indx
   END DO

   !- calc circumcenter
   lcl_complex(1)%centers = g_coord_pp(:,:,iel)
   DO k=2,dim_cmplx+1; CALL calc_circumcenters(k); END DO

   !- calc primal edge and dual area
   CALL calc_prml_unsgnd_vlm(2); CALL calc_dual_vlm(2);  CALL calc_prml_dir()

   kcx=zero
   DO i=1,nod
     DO j=1,lcl_complex(1)%num_cobndry(i)
       k = lcl_complex(1)%cobndry(i)%indx(j)
       IF (ABS(lcl_complex(2)%dual_volume(k) / &
          max(lcl_complex(2)%prml_volume(k),small))<smalls) CYCLE
       m = lcl_complex(2)%bndry(k)%indx(1)
       IF (i==m) m = lcl_complex(2)%bndry(k)%indx(2)
       ka = DOT_PRODUCT(kay,abs(lcl_complex(2)%prml_dir(k,:)))
       kcx(i,m) = kcx(i,m) - ka*lcl_complex(2)%dual_volume(k) / &
          max(lcl_complex(2)%prml_volume(k),small)
       kcx(i,i) = kcx(i,i) + ka*lcl_complex(2)%dual_volume(k) / &
          max(lcl_complex(2)%prml_volume(k),small)
     END DO
   END DO;  storkc_pp(:,:,iel)=kcx

   !- Deallocate variables
   DO k=1,dim_cmplx+1;
     IF (ALLOCATED(lcl_complex(k)%recv_indx))   DEALLOCATE(lcl_complex(k)%recv_indx)
     IF (ALLOCATED(lcl_complex(k)%num_bndry))   DEALLOCATE(lcl_complex(k)%num_bndry)
     IF (ALLOCATED(lcl_complex(k)%bndry))       DEALLOCATE(lcl_complex(k)%bndry)
     IF (ALLOCATED(lcl_complex(k)%num_cobndry)) DEALLOCATE(lcl_complex(k)%num_cobndry)
     IF (ALLOCATED(lcl_complex(k)%cobndry))     DEALLOCATE(lcl_complex(k)%cobndry)
     IF (ALLOCATED(lcl_complex(k)%node_indx))   DEALLOCATE(lcl_complex(k)%node_indx)
   END DO
 END DO elements_1
!------------------ build the diagonal preconditioner --------------------
 ALLOCATE(diag_precon_tmp(ntot,nels_pp)); diag_precon_tmp=zero
 elements_1a: DO iel=1,nels_pp
   DO i=1,ndof
     diag_precon_tmp(i,iel)=diag_precon_tmp(i,iel)+storkc_pp(i,i,iel)
   END DO
 END DO elements_1a; CALL scatter(diag_precon_pp,diag_precon_tmp)
 DEALLOCATE(diag_precon_tmp)
!------------- read in fixed freedoms and assign to equations ------------
 IF(fixed_freedoms>0) THEN
   ALLOCATE(node(fixed_freedoms),no_pp_temp(fixed_freedoms),             &
     val_f(fixed_freedoms),no(fixed_freedoms),sense(fixed_freedoms))
   node=0; no=0; no_pp_temp=0; sense=0; val_f=zero
   CALL read_fixed(argv,numpe,node,sense,val_f)
   CALL find_no2(g_g_pp,g_num_pp,node,sense,no)
   CALL reindex(ieq_start,no,no_pp_temp,fixed_freedoms_pp,               &
     fixed_freedoms_start,neq_pp)
   ALLOCATE(no_f_pp(fixed_freedoms_pp),store_pp(fixed_freedoms_pp))
   no_f_pp=0; store_pp=zero; no_f_pp=no_pp_temp(1:fixed_freedoms_pp)
   DEALLOCATE(node,no,sense,no_pp_temp)
 END IF
 IF(fixed_freedoms==0) fixed_freedoms_pp=0; IF(nr>0) DEALLOCATE(rest)
!------------ read in loaded freedoms and get starting r_pp --------------
 IF(loaded_freedoms>0) THEN
   ALLOCATE(node(loaded_freedoms),val(nodof,loaded_freedoms),            &
     no_pp_temp(loaded_freedoms)); val=zero; node=0
   CALL read_loads(argv,numpe,node,val)
   CALL reindex(ieq_start,node,no_pp_temp,loaded_freedoms_pp,            &
     loaded_freedoms_start,neq_pp); ALLOCATE(no_pp(loaded_freedoms_pp))
   no_pp=no_pp_temp(1:loaded_freedoms_pp); DEALLOCATE(no_pp_temp)
   DO i = 1, loaded_freedoms_pp
     r_pp(no_pp(i)-ieq_start+1) = val(1,loaded_freedoms_start+i-1)
   END DO; q=SUM_P(r_pp); DEALLOCATE(node,val)
 END IF
!-------------------------- invert preconditioner ------------------------
 IF(fixed_freedoms_pp>0) THEN
   DO i=1,fixed_freedoms_pp; j=no_f_pp(i)-ieq_start+1
     diag_precon_pp(j)=diag_precon_pp(j)+penalty
     store_pp(i)=diag_precon_pp(j)
   END DO
 END IF;  diag_precon_pp = 1._iwp/diag_precon_pp
!---------------- initialise preconditioned conjugate gradient -----------
 IF(fixed_freedoms_pp>0) THEN
   DO i=1,fixed_freedoms_pp
     j=no_f_pp(i)-ieq_start+1; k=fixed_freedoms_start+i-1
     r_pp(j)=store_pp(i)*val_f(k)
   END DO
 END IF;  d_pp=diag_precon_pp*r_pp; p_pp=d_pp; x_pp=zero
!--------------- preconditioned conjugate gradient iterations ------------
 iters=0; timest(3)=elap_time()
 iterations: DO
   iters=iters+1; u_pp=zero; pmul_pp=zero; utemp_pp=zero
   CALL gather(p_pp,pmul_pp)
   elements_2 : DO iel = 1, nels_pp
     utemp_pp(:,iel) = MATMUL(storkc_pp(:,:,iel),pmul_pp(:,iel))
   END DO elements_2 ; CALL scatter(u_pp,utemp_pp)
   IF(fixed_freedoms_pp>0) THEN
     DO i=1,fixed_freedoms_pp
       j=no_f_pp(i)-ieq_start+1; u_pp(j)=p_pp(j)*store_pp(i)
     END DO
   END IF
   up=DOT_PRODUCT_P(r_pp,d_pp); alpha=up/DOT_PRODUCT_P(p_pp,u_pp)
   xnew_pp=x_pp+p_pp*alpha; r_pp=r_pp-u_pp*alpha
   d_pp=diag_precon_pp*r_pp; beta=DOT_PRODUCT_P(r_pp,d_pp)/up
   p_pp=d_pp+p_pp*beta; CALL checon_par(xnew_pp,tol,converged,x_pp)
   IF(converged .OR. iters==limit) EXIT
 END DO iterations
 timest(4)=elap_time()
 IF(numpe==it)THEN
   WRITE(11,'(A,I5)')"The number of iterations to convergence was ",iters
   WRITE(11,'(A,E12.4)')"The total load is ",q
   WRITE(11,'(A)')   "The potentials are:"
   WRITE(11,'(A)') " Freedom       Potential"
   DO i=1,4
     WRITE(11,'(I8,A,E12.4)') nres+i-1, "     ", xnew_pp(is+i-1)
   END DO
   WRITE(11,'(A,F10.4)')"Integration, preconditioning and loading took ",&
     timest(3)-timest(2)
   WRITE(11,'(A,F10.4)') "Time spent in the solver was ",                &
     timest(4)-timest(3)
 END IF
!------------------- output potentials for ParaView ----------------------
 IF(numpe==1)THEN; WRITE(ch,'(I6.6)') numpe
   OPEN(12,file=argv(1:nlen)//".ensi.NDPTL-"//ch,status='replace',       &
     action='write')
   WRITE(12,'(A)') "Alya Ensight Gold --- Scalar per-node variable file"
   WRITE(12,'(A/A/A)') "part", "    1","coordinates"
 END IF
 CALL calc_nodes_pp(nn,npes,numpe,node_end,node_start,nodes_pp)
 ALLOCATE(ptl_pp(nodes_pp*ndim))
 ptl_pp=zero; utemp_pp=zero; CALL gather(xnew_pp(1:),utemp_pp)
 CALL scatter_nodes(npes,nn,nels_pp,g_num_pp,nod,nodof,nodes_pp,         &
   node_start,node_end,utemp_pp,ptl_pp,1)
 CALL dismsh_ensi_p(12,1,nodes_pp,npes,numpe,1,ptl_pp)
 IF(numpe==it)                                                           &
   WRITE(11,'(A,F10.4)') "This analysis took ",elap_time()-timest(1)
 IF(numpe==it) CLOSE(11); IF(numpe==1) CLOSE(12); CALL SHUTDOWN()
END PROGRAM pg123

max_stars_repo_name: ParaFEM/ParaGEMS

max_issues_repo_name: ParaFEM/ParaGEMS

max_forks_repo_name: ParaFEM/ParaGEMS


text:
      subroutine tltrm(kptbl)
      use tfstk
      use tmacro
      use sad_main
      use tparastat
      implicit none
      integer*4 i
      integer*4 kptbl(np0,6)
      call tclrpara
      do i=1,np0
c         write(*,*)'tltrm ',i,np0,kptbl(i,3)
c     $       ilist(1,abs(kptbl(i,3)))
        if(kptbl(i,3) .ne. 0)then
          call tfree(int8(abs(kptbl(i,3))))
          kptbl(i,3)=0
        endif
      enddo
      return
      end

max_stars_repo_name: noboruatkek/SAD

max_issues_repo_name: noboruatkek/SAD

max_forks_repo_name: noboruatkek/SAD


text:
c Identical to vegas.for, except for the fact that a larger number of
c integration variables (58 rather than 10) is allowed here
c
c
c-
c Integration-histogramming package using vegas.
c Calls:
c      subroutine delete(fname)         ! deletes the file fname
c      character * (*) fname
c
c      subroutine time(ctime)           ! returns a string with the time
c      character * 8 ctime
c
c      subroutine date(cdate)           ! returns a string with the date
c      character * 9 cdate
c

      subroutine strnum(string,num)
c- writes the number num on the string string starting at the blank
c- following the last non-blank character
      character * (*) string
      character * 20 tmp
      l = len(string)
      write(tmp,'(i15)')num
      j=1
      dowhile(tmp(j:j).eq.' ')
        j=j+1
      enddo
      ipos = istrl(string)
      ito = ipos+1+(15-j)
      if(ito.gt.l) then
         write(*,*)'error, string too short'
         write(*,*) string
         stop
      endif
      string(ipos+1:ito)=tmp(j:)
      end

      function istrl(string)
c returns the position of the last non-blank character in string
      character * (*) string
      i = len(string)
      dowhile(i.gt.0.and.string(i:i).eq.' ')
         i=i-1
      enddo
      istrl = i
      end

      subroutine strcat(str1,str2,str)
c concatenates str1 and str2 into str. Ignores trailing blanks of str1,str2
      character *(*) str1,str2,str
      l1=istrl(str1)
      l2=istrl(str2)
      l =len(str)
      if(l.lt.l1+l2) then
          write(*,*) 'error: l1+l2>l in strcat'
          write(*,*) 'l1=',l1,' str1=',str1
          write(*,*) 'l2=',l2,' str2=',str2
          write(*,*) 'l=',l
          stop
      endif
      if(l1.ne.0) str(1:l1)=str1(1:l1)
      if(l2.ne.0) str(l1+1:l1+l2)=str2(1:l2)
      if(l1+l2+1.le.l) str(l1+l2+1:l)= ' '
      end

      block data skiphistos
c If integrate will be allowed to save and restore histograms
c (the default mode of operation till Nov 2013), then set:
c  skiph=.false.
c here. Otherwise, if skiph=.true., no information on the histograms
c will be saved in the .sv1 and .sv2 by the routine integrate()
      implicit none
      logical skiph
      common/cskiph/skiph
      data skiph/.true./
      end

      subroutine
     # integrate(init,sig,str,n,inew,idim,icalls,avgi,sd,chi2a,isave)
      implicit real * 8 (a-h,o-z)
      external init,sig
      common/bveg1/xl(58),xu(58),acc,ndim,ncall,itmx,nprn
      common/bveg2/xi(50,58),si,si2,swgt,schi,ndo,it
      common/newver/newver
      character * (*) str
      character * 7 statf,newver
      character * 70 fname
      logical isave
c- if inew = 0 ignore
      if(inew.eq.0) return
c- add str to list of save files
      if(isave) then
         fname = str
         call addfil(fname)
      endif
c- initialize histograms
      call init
c
      if(inew.ge.2) then
         call resume(str,init,'IN')
c also fixes the value of it
         call vegas4(sig,avgi,sd,chi2a)
      endif
c
      if(inew.eq.10) then
c- save the file '.sv2' in ASCII format
         call strcat(str,'.sv2',fname)
c- ndim = idim for upward compatibility with old restart files
         ndim = idim
         call savea(fname,newver)
         return
      endif
c
      if(inew.eq.11) then
c- save the file '.sv1' in standard format
         call strcat(str,'.sv1',fname)
         call save(fname,newver)
         return
      endif
c
c Initialize parameters for Vegas
      do j=1,58
          xl(j) = 0
          xu(j) = 1
      enddo
      acc = -1
      ndim = idim
      ncall = icalls
      itmx = 1
c
c for negative n abs(n) is the number of iterations to reach,
c for positive n, n is the number of iterations to perform
      if(n.lt.0) then
         if(inew.eq.1.or.inew.eq.3) then
            nsteps = abs(n)
         elseif(inew.eq.2) then
            nsteps = abs(n) - it
         endif
      else
         nsteps = n
      endif
c create new versions of save file
cRF
c$$$      statf = newver
      statf = 'UNKNOWN'
c----------------------------------------------------------
      do j=1,nsteps
         if(j.eq.1) then
            if(inew.eq.1) then
                  call vegas(sig,avgi,sd,chi2a)
            elseif(inew.eq.2) then
                  call vegas2(sig,avgi,sd,chi2a)
            elseif(inew.eq.3) then
                  call init
                  call vegas1(sig,avgi,sd,chi2a)
            else
                  write(*,*) 'irregular inew flag ',inew
                  stop
            endif
         else
            call vegas3(sig,avgi,sd,chi2a)
         endif
         call accum
         if(isave) then
            call strcat(str,'.sv2',fname)
            call save(fname,statf)
            call strcat(str,'.sv1',fname)
            call save(fname,statf)
            call strcat(str,'.sv2',fname)
            call delete(fname)
            statf = 'UNKNOWN'
         endif
      enddo
      return
      entry integrms1
c # of iterations message
      write(*,*)
      write(*,*)'if you enter a number of iterations n>0'
      write(*,*)'the program will execute n iterations'
      write(*,*)'if n<0, it will execute iterations until'
      write(*,*)'total number of iterations performed = abs(n)'
      write(*,*)'(it makes a difference only if in restart mode)'
      write(*,*)
      return
      entry integrms2
c inew values
      write(*,*)
      write(*,*) 'enter 0 to exclude,1 for new run, 2 to restart'
      write(*,*) '3 to restart keeping the grid but reset histo''s '
      write(*,*) '10 to produce ASCII save files from standard ones'
      write(*,*) '11 to produce standard save files from ASCII ones'
      write(*,*) '(10/11 used to transport save files across machines)'
      write(*,*)
      return
      entry nopr(nnn)
      nprn = nnn
      end

      subroutine resume(string,init,fl)
      character * (*) string, fl *2
      character * 70 fname
      if(fl.eq.'IN') call init
      call strcat(string,'.sv1',fname)
      call restart(fname,itmp)
      if(itmp.eq.1.and.fl.eq.'NO') goto 99
      if(itmp.eq.0) return
c- retry with spare save file
      if(fl.eq.'IN') call init
      call strcat(string,'.sv2',fname)
      call restart(fname,itmp)
      if(itmp.eq.1.and.fl.eq.'NO') goto 99
      if(itmp.eq.0) return
c- retry with ASCII spare save file
      if(fl.eq.'IN') call init
      call strcat(string,'.sv2',fname)
      call restara(fname,itmp)
      if(itmp.eq.1.and.fl.eq.'NO') goto 99
      if(itmp.eq.0) return
c- give up.
      write(*,*) 'unusable restart file'
      stop
99    continue
c In this case the accumulated values in the histograms common block
c could have been altered in the attempt to read the save file.
c Abort. (if itmp=2 this is not the case)
      write(*,*) 'corrupted file',fname,'.'
      write(*,'(a,/)')
     # 'The attempt to read the corrupted file may have corrupted',
     # 'the histogram common block. Delete the corrupted file, or',
     # 'rerun one more iteration, or run with nit=11.'
      stop
      end

      subroutine savea(name,statf)
      include "dbook.inc"
      PARAMETER (NMB=NPLOTS)
      character * 80 runstr,string
      character * 3 books(nmb)
      real * 8 xl,xu,acc
      common/bveg1/xl(58),xu(58),acc,ndim,ncall,itmx,nprn
      real * 8 xi,si,si2,swgt,schi
      common/bveg2/xi(50,58),si,si2,swgt,schi,ndo,it

      common/runstr/runstr
      character * (*) name, statf
      logical skiph
      common/cskiph/skiph
c
c Salva gli istogrammi in uso (marcati 'YES' in book(j))
c e tutte le informazioni per vegas, il seme del numero
c casuale, etc.
c Marchia la fine del file con i numeri acos(-1.),exp(1.),log(10.),
c
      open(unit=97,file=name,status=statf)
101   format(a)
102   format(25a)
103   format(6d13.7)
104   format(6i12)
105   format(4d20.14)
c version of save file
      string = 'SAVEA,VER 8-3-92'
      write(97,101) string
c runstring
      write(97,101)  runstr
c booked histograms
      if(.not.skiph)then
      write(97,104)  nmb
      write(97,102)  (book(j),j=1,nmb)
c Mbook block
      do j=1,nmb
         if(book(j).eq.'YES') then
            write(97,103)
     &      (hist(j,k),k=1,nbin(j)),uscore(j),oscore(j)
            write(97,104)
     &      (ihis(j,k),k=1,nbin(j)),iuscore(j),ioscore(j),ient(j)
         endif
      enddo
      endif
c Vegas block
      write(97,104)
     & ndo,it,num1,num2,ndim
      write(97,105)
     & ((xi(j,k),j=1,ndo),k=1,ndim),si,si2,swgt,schi

      a=acos(-1.)
      b=exp(1.)
      c=log(10.)
c end mark
      write(97,103)
     & a,b,c
      close(97)
      return
c
c Loads the histograms found in the (ascii) save file into the common block.
c Previous values are overwritten.
c Does not modify histograms not present in the file.
c If all goes well, it riturns iflag= 0.
c If there are errors (e.g. IO errors, the mark at the end of the file
c is not found, etc.) and in the input operation the common blocks
c may have been corrupted, returns iflag=1.
c If there are errors, but no common block has been altered, returns iflag=2.
c
      entry restara(name,iflag)
      open(unit=97,file=name,status='OLD',err=2)
c version of save file
      string = ' '
      read(97,101,err=2,end=2) string
      if(string.ne.'SAVEA,VER 8-3-92') goto 2
c run string
      string = ' '
      read(97,101,err=2,end=2) string
c marked histograms
      if(.not.skiph)then
      read(unit=97,fmt=104,err=2,end=2) nhs
      if(nhs.gt.nmb) then
         write(*,*) 'Error: save files had nmb=',nhs
         write(*,*) 'Running program has nmb=',nmb,'<',nhs
         write(*,*) 'Make nmb >= ',nhs,' in running program'
         stop
      endif
      read(unit=97,fmt=102,err=1,end=1) (books(j),j=1,nhs)
c Mbook block
      do j=1,nhs
         if(books(j).eq.'YES') then
            read(unit=97,fmt=103,err=1,end=1)
     &      (hist(j,k),k=1,nbin(j)),uscore(j),oscore(j)
            read(unit=97,fmt=104,err=1,end=1)
     &      (ihis(j,k),k=1,nbin(j)),iuscore(j),ioscore(j),ient(j)
         endif
      enddo
      endif
c vegas block
      read(unit=97,fmt=104,err=1,end=1)
     & ndo,it,num1,num2,ndim
      read(unit=97,fmt=105,err=1,end=1)
     & ((xi(J,k),j=1,ndo),k=1,ndim),si,si2,swgt,schi
c end mark
      read(unit=97,fmt=103,err=1,end=1)
     & a0,b0,c0
c
      tiny = .1e-5
      a = abs(a0-acos(-1.))
      b = abs(b0-exp(1.))
      c = abs(c0-log(10.))
      if( a.gt.tiny .or. b.gt.tiny .or. c.gt.tiny ) goto 1
      close(unit=97,iostat=ios)
      iflag = 0
      if(string.ne.runstr) then
         write(6,*)
     #  'You are using save files belonging to a different run,'
         write(6,*) string
         write(6,*) 'instead of'
         write(6,*) runstr
      endif
      write(6,*) string
      return
 1    continue
      iflag = 1
      goto 3
 2    continue
      iflag = 2
      goto 3
 3    continue
      close(unit=97,iostat=ios)
      end

      subroutine save(name,statf)
      include "dbook.inc"
      PARAMETER (NMB=NPLOTS)
      character * 80 runstr,string
      character * 3 books(nmb)
      real * 8 xl,xu,acc
      common/bveg1/xl(58),xu(58),acc,ndim,ncall,itmx,nprn
      real * 8 xi,si,si2,swgt,schi
      common/bveg2/xi(50,58),si,si2,swgt,schi,ndo,it

      common/runstr/runstr
      character * (*) name, statf
      logical skiph
      common/cskiph/skiph
c
c Salva gli istogrammi in uso (marcati 'YES' in book(j))
c e tutte le informazioni per vegas, il seme del numero
c casuale, etc.
c Marchia la fine del file con i numeri acos(-1.),exp(1.),log(10.),
c
      open(unit=97,file=name,form='UNFORMATTED',status=statf)
c version of save program
      string = 'SAVE,VER 8-3-92'
      write(97) string
c run string
      write(97)  runstr
c Mbook Block
      if(.not.skiph)then
      write(97) nmb
      write(97)  (book(j),j=1,nmb)
      do j=1,nmb
         if(book(j).eq.'YES') then
           write(97)
     &     (hist(j,k),k=1,nbin(j))
     &     ,uscore(j),oscore(j)
     &     ,(ihis(j,k),k=1,nbin(j))
     &     ,iuscore(j),ioscore(j)
     &     ,ient(j)
         endif
      enddo
      endif
c Vegas Block
      write(97)
     & ndo,it,num1,num2,ndim
      write(97)
     & ((xi(J,k),j=1,ndo),k=1,ndim),si,si2,swgt,schi
c
      a=acos(-1.)
      b=exp(1.)
      c=log(10.)
      write(97) a,b,c
      close(97)
      return
c
c Loads the histograms found in the file into the common block.
c Previous values are overwritten.
c Does not modify histograms not present in the file.
c If all goes well, it riturns iflag= 0.
c If there are errors (e.g. IO errors, the mark at the end of the file
c is not found, etc.) and in the input operation the common blocks
c may have been corrupted, returns iflag=1.
c If there are errors, but no common block has been altered, returns iflag=2.
c
      entry restart(name,iflag)
      open(unit=97,file=name,form='unformatted',status='OLD',err=2)
c version of save file
      string = ' '
      read(unit=97,err=2,end=2) string
      if(string.ne.'SAVE,VER 8-3-92') goto 2
c run string
      string = ' '
      read(unit=97,err=2,end=2) string
c Mbook block
      if(.not.skiph)then
      read(97,err=2,end=2) nhs
      if(nhs.gt.nmb) then
         write(*,*) 'Error: save files had nmb=',nhs
         write(*,*) 'Running program has nmb=',nmb,'<',nhs
         write(*,*) 'Make nmb >= ',nhs,' in running program'
         stop
      endif
      read(97,err=1,end=1) (books(j),j=1,nhs)
      do j=1,nhs
         if(books(j).eq.'YES') then
            read(97,err=1,end=1)
     &      (hist(j,k),k=1,nbin(j))
     &      ,uscore(j),oscore(j)
     &      ,(ihis(j,k),k=1,nbin(j))
     &      ,iuscore(j),ioscore(j)
     &      ,ient(j)
         endif
      enddo
      endif
c Vegas Block
      read(97,err=1,end=1)
     & ndo,it,num1,num2,ndim
      read(97,err=1,end=1)
     & ((xi(J,k),j=1,ndo),k=1,ndim),si,si2,swgt,schi
c
      a0 = 0
      b0 = 0
      c0 = 0
      read(97,err=1,end=1)   a0,b0,c0
      a=acos(-1.)
      b=exp(1.)
      c=log(10.)
      if( a.ne.a0 .or. b.ne.b0 .or. c.ne.c0 ) goto 1
      close(97)
      iflag = 0
      if(string.ne.runstr) then
         write(6,*)
     #   'You are using save files belonging to a different run,'
         write(6,*) string
         write(6,*) 'instead of'
         write(6,*) runstr
      endif
      write(6,*) string
      return
 1    continue
      iflag = 1
      goto 3
 2    continue
      iflag = 2
      goto 3
 3    continue
      close(unit=97,iostat=ios)
      end

       block data vegas0
       implicit  real*8 (a-h,o-z)
       common/bveg1/xl(58),xu(58),acc,ndim,ncall,itmx,nprn
       common/bveg2/xi(50,58),si,si2,swgt,schi,ndo,it
       data ncall/10000/,itmx/10/,nprn/ 1000/,acc/-1.d0/,
     1   xl/58*1.d-3/,
     2   xu/58*1.d0/
       data xi/2900*1.d0/
       end
C
C
C      NCALL IS THE NUMBER OF CALLS TO VEGAS.
C      NPRN >  0 VEGAS PRINTS THE RESULTS OF EACH ITERATION.
C      NPRN 0 VEGAS PRINTS NOTHING.
C      NPRN < 0 VEGAS PRINTS ALL.
C      XL(I) IS LOWER INTEGRATION LIMIT ON I TH AXIS.
C      XU(I) IS UPPER INTEGRATION LIMIT ON I THE AXIS.
c
         subroutine vegas(fxn,avgi,sd,chi2a)
c
c   routine performs n dim monte carlo inte
c written by p lepage
c
         implicit real*8 (a-h,o-z)
         implicit integer*4 (i-n)
         common/bveg1/xl(58),xu(58),acc,ndim,ncall,itmx,nprn
         common/bveg2/xi(50,58),si,si2,swgt,schi,ndo,it
         dimension d(50,58),di(50,58),xin(50),r(50),dx(58),dt(58),
     1   x(58),kg(58),ia(58)
       dimension RAND(58)
         data ndmx/50/,alph/1.5d0/,one/1.d0/,mds/1/
         ndo=1
         do 1 j=1,ndim
 1       xi(1,j)=one
c
         entry vegas1(fxn,avgi,sd,chi2a)
c       initialises  cumulative  variables but not grid
         it=0
         si=0.
         si2=si
         swgt=si
         schi=si
c
         entry vegas2(fxn,avgi,sd,chi2a)
c        no initialisation
         nd=ndmx
         ng=1
         if(mds.eq.0)go to 2
         ng=(ncall/2.)**(1./ndim)
         mds=1
         if((2*ng-ndmx).lt.0)go to 2
         mds=-1
         npg=ng/ndmx+1
         nd=ng/npg
         ng=npg*nd
 2       k=ng**ndim
         npg=ncall/k
         if(npg.lt.2)npg=2
         calls=npg*k
         dxg=one/ng
         dv2g=(calls*dxg**ndim)**2/npg/npg/(npg-one)
         xnd=nd
         ndm=nd-1
         dxg=dxg*xnd
         xjac=one/calls
         do 3 j=1,ndim
         dx(j)=xu(j)-xl(j)
 3       xjac=xjac*dx(j)
c
c    rebin preserving bin density
c
         if(nd.eq.ndo)go to 8
         rc=ndo/xnd
         do 7 J=1,ndim
         k=0
         xn=0.
         dr=xn
         i=k
 4       k=k+1
         dr=dr+one
         xo=xn
         xn=xi(k,j)
 5       if(rc.gt.dr)go to 4
         i=i+1
         dr=dr-rc
         xin(i)=xn-(xn-xo)*dr
         if(i.lt.ndm)go to 5
         do 6 i=1,ndm
 6       xi(i,j)=xin(i)
 7       xi(nd,j)=one
         ndo=nd
c
 8       if(nprn.ne.0)write(6,200)ndim,calls,it,itmx,acc
     1   ,mds,nd,(xl(j),xu(j),j=1,ndim)
c
         entry vegas3(fxn,avgi,sd,chi2a)
c         main integration loop
 9       it=it+1
          ti=0.
         tsi=ti
         do 10 j=1,ndim
         kg(j)=1
         do 10 i=1,nd
         d(i,j)=ti
 10      di(i,j)=ti
c
 11      fb=0.
         f2b=fb
         k=0
 12      k=k+1
       call randa(ndim,rand)
         wgt=xjac
         do 15 j=1,ndim
         xn=(kg(j)-rand(j))*dxg+one
         ia(j)=xn
         if(ia(j).gt.1)go to 13
         xo=xi(ia(j),j)
         rc=(xn-ia(j))*xo
         go to 14
13       xO=xi(ia(j),j)-xi(ia(j)-1,j)
         rc=xi(ia(j)-1,j)+(xn-ia(j))*xo
 14      x(j)=xl(j)+rc*dx(j)
 15      wgt=wgt*xo*xnd
c
         f=wgt
         f=f*fxn(x,wgt)
         f2=f*f
         fb=fb+f
         f2b=f2b+f2
         do 16 j=1,ndim
         di(ia(j),j)=di(ia(j),j)+f
 16      if(mds.ge.0)d(ia(j),J)=d(ia(j),J)+f2
         if(k.lt.npg) go to 12
c
888    FORMAT(1X,'F',G14.6,'F2',G14.6,'FB',G14.6,'F2B',G14.6)
         f2b= sqrt(f2b*      NPG)
         f2b=(f2b-fb)*(f2b+fb)
1661   FORMAT(1X,'F2B',G14.6,'NPG',  I10)
         ti=ti+fb
         tsi=tsi+f2b
33     FORMAT(1X,'TSI',G14.6,'F2B',G14.6)
         if(mds.ge.0)go to 18
         do 17 j=1,ndim
 17      d(ia(j),j)=d(ia(j),j)+f2b
 18      k=ndim
 19      kg(k)=mod(kg(k),ng)+1
         if(kg(k).ne.1)go to 11
         k=k-1
         if(k.gt.0)go to 19
c
c final results for this iteration
c
        tsi=tsi*dv2g
        ti2=ti*ti
88     format(1x,'tsi',g14.6)
       if(tsi.eq.0.d0)then
         rpp=1.d15
       else
         rpp=abs(ti2/tsi)
       endif
       if(rpp.lt.1.d14) then
        wgt=ti2/tsi
        si=si+ti*wgt
        si2=si2+ti2
        swgt=swgt+wgt
        schi=schi+ti2*wgt
995    FORMAT(1X,'SWGT',G14.6,'SI2',G14.6)
        avgi=si/swgt
        sd=swgt*it/si2
        chi2a=sd*(schi/swgt-avgi*avgi)/(it-.999)
        sd=dsqrt(one/sd)
       else
        write(*,*) '***VEGAS WARNING***: zero error!'
        write(*,*) ' we guess that integral is exact '
        avgi=ti
        si=ti
        si2=0.d0
        sd=0
        chi2a=-1
        return
       endif
c
        if(nprn.eq.0)go to 21
        tsi=dsqrt(tsi)
        write(6,201)it,ti,tsi,avgi,sd,chi2a
cRF: also update the grids for the integer sum
        call regrid_MC_integer
        if(nprn.ge.0)go to 21
        do 20 j=1,ndim
 20     write(6,202) j,(xi(i,j),di(i,j),d(i,j),i=1,nd)
c
c      refine grid
c
 21     do 23 j=1,ndim
        xo=d(1,j)
        xn=d(2,j)
        d(1,j)=(xo+xn)/2.
        dt(j)=d(1,j)
        do 22 i=2,ndm
        d(i,j)=xo+xn
        xo=xn
        xn=d(i+1,j)
        d(i,j)=(d(i,j)+xn)/3.
 22     dt(j)=dt(j)+d(i,j)
        d(nd,j)=(xn+xo)/2.
 23     dt(j)=dt(j)+d(nd,j)
c
        do 28 j=1,ndim
        rc=0.
        do 24 i=1,nd
        r(i)=0.
        if(d(i,j).le.0.)go to 24
        xo=dt(j)/d(i,j)
        r(i)=((xo-one)/xo/dlog(xo))**alph
 24     rc=rc+r(i)
        rc=rc/xnd
        k=0
        xn=0.
        dr=xn
        i=k
 25     k=k+1
        dr=dr+r(k)
        xo=xn
        xn=xi(k,j)
 26     if(rc.gt.dr)go to 25
        i=i+1
        dr=dr-rc
        xin(i)=xn-(xn-xo)*dr/r(k)
        if(i.lt.ndm)go to 26
        do 27 i=1,ndm
 27     xi(i,j)=xin(i)
 28     xi(nd,j)=one
c
        if(it.lt.itmx.and.acc*dabs(avgi).lt.sd)go to 9
 200    format(1X,'0input parameters for vegas:  ndim=',i3,
     1  '   ncall=',f8.0/28x,'  it=',i5,'    itmx=',i5/28x,
     2  '  acc=',g9.3/28x,'  mds=',i3,'     nd=',i4/28x,
     3  '  (xl,xu)=',(t40,'( ',g12.6,' , ',g12.6,' )'))
 201    format(///' integration by vegas' / '0iteration no.',i5,
     1  ':  integral=',g14.8/21x,'std dev =',g10.4 /
     2  ' accumulated results:   integral=',g14.8/
     3  24x,'std dev =',g10.4 / 24x,'chi**2 per it''n =',g10.4)
 202    format(1X,'0data for axis',i2,/,' ',6x,'x',7x,'  delt i ',
     1  2x,'conv','ce   ',11x,'x',7x,'  delt i ',2x,'conv','ce  '
     2  ,11x,'x',7x,'   delt i ',2x,'conv','CE  ',/,
     3  (1X,' ',3g12.4,5x,3g12.4,5x,3g12.4))
        return
        entry vegas4(fxn,avgi,sd,chi2a)
        if(si2.eq.0.d0)then
          avgi=si
          sd=0
          chi2a=-1
        else
          avgi=si/swgt
          sd=swgt*it/si2
          chi2a=sd*(schi/swgt-avgi*avgi)/(it-.999)
          sd=dsqrt(one/sd)
        endif
        if(nprn.ne.0)then
           write(6,201)it,0.d0,0.d0,avgi,sd,chi2a
        endif
        return
        end

c        subroutine save(ndim)
c        implicit real*8 (a-h,o-z)
c       implicit integer*4 (i-n)
c        common/bveg2/xi(50,58),si,si2,swgt,schi,ndo,it
c
c      stores vegas data   (unit 7) for later initialisation
c
c        write(7,200) ndo,it,si,si2,swgt,schi,
c     1       ((xi(i,j),i=1,ndo),j=1,ndim)
c        return
c        entry restr(ndim)
c
c         enters initialisation data for vegas
c
c        read(7,200) ndo,it,si,si2,swgt,schi,
c     1    ((xi(i,j),i= 1,ndo),j=1,ndim)
c 200    format(2i8,4z16/(5z16))
c        return
c        end

      subroutine randa(n,rand)
      implicit double precision (a-h,o-z)
      common/caso/caso(5)
      dimension rand(58)
      real*8 ran2
      external ran2
      do 1 i=1,n
      rand(i)=ran2()
1     continue
      do 2 i=1,5
      caso(i)=ran2()
2     continue
      return
      end

      subroutine sysdep
c use newver='NEW' for vaxes, 'UNKNOWN' in other machines,
      character * 7 newver
      common/newver/newver
      newver = 'UNKNOWN'
      end

      subroutine delete(fname)
      character * 80 str
      character * (*) fname
      l = len(fname)
      k = 1
      dowhile(fname(k:k).eq.' '.and.k.lt.l)
         k = k+1
      enddo
      dowhile(fname(l:l).eq.' '.and.l.gt.k)
         l = l-1
      enddo
      if(l-k.gt.70) then
         write(*,*) 'delete: filename > 70 chars not allowed'
         stop
      endif
      if(l.eq.k) then
         write(*,*) 'delete: void filename'
         stop
      endif
      str(1:) = 'rm -f '
      str(7:) = fname(k:l)
      call system(str)
      end

c
c
c A simple integration wrapper
c
c
      subroutine xinteg(sig,indim,initn,incalls,avg,stddev)
      implicit real * 8 (a-h,o-z)
      common/bveg1/xl(58),xu(58),acc,ndim,ncall,itmx,nprn
      common/bveg2/xi(50,58),si,si2,swgt,schi,ndo,it
      common/cwww/w1max,w1vgs,w1evt
      common/ciweight/iweight
      external sig
c
      do j=1,58
        xl(j)=0.d0
        xu(j)=1.d0
      enddo
c Set nprn equal to zero to print nothing
      nprn=1
      acc=-1.d0
      ndim=indim
      ncall=incalls
      itmx=1
      do j=1,initn
        if(j.eq.1)then
          call vegas(sig,avgi,sd,chi2a)
        else
          call vegas3(sig,avgi,sd,chi2a)
        endif
      enddo
      avg=avgi
      stddev=sd
      return
      end

ext: for

max_stars_repo_name: JaySandesara/madminer

max_issues_repo_name: JaySandesara/madminer

max_forks_repo_name: JaySandesara/madminer


text:
!> \file rhoinsitu.f90
!! \BRIEF
!> Module with rhoinsitu subroutine - compute in situ density from S, Tis, P
MODULE mrhoinsitu
CONTAINS
!>     Compute in situ density from salinity (psu), in situ temperature (C), & pressure (db).
!!     This subroutine is needed because rho is a function (using scalars not arrays)
SUBROUTINE rhoinsitu(salt, tempis, pdbar, N, rhois)

  !     Purpose:
  !     Compute in situ density from salinity (psu), in situ temperature (C), & pressure (db)
  !     Needed because rho is a function

  USE msingledouble
  USE mrho
  IMPLICIT NONE

  INTEGER :: N

! INPUT variables
  ! salt   = salinity [psu]
  ! tempis = in situ temperature [C]
  ! pdbar  = pressure [db]

  !> salinity [psu]
  REAL(kind=rx), INTENT(in), DIMENSION(N) :: salt
  !> in situ temperature [C]
  REAL(kind=rx), INTENT(in), DIMENSION(N) :: tempis
  !> pressure [db]
  REAL(kind=rx), INTENT(in), DIMENSION(N) :: pdbar
!!!f2py optional , depend(salt) :: n=len(salt)
!f2py depend(salt) :: n=len(salt)

! OUTPUT variables:
  ! rhois  = in situ density

  !> in situ density [kg/m3]
  REAL(kind=rx), INTENT(out), DIMENSION(N) :: rhois

! Local variables
  INTEGER :: i

! REAL(kind=rx) ::  rho
! EXTERNAL rho

  DO i = 1,N
     rhois(i) = rho(salt(i), tempis(i), pdbar(i)/10.)
  END DO

  RETURN
END SUBROUTINE rhoinsitu
END MODULE mrhoinsitu

max_stars_repo_name: tomaslovato/mocsy

max_issues_repo_name: tomaslovato/mocsy

max_forks_repo_name: tomaslovato/mocsy


text:
!-----------------------------------------------------------------------
!     Function:      db_midpoint
!     Authors:       S. Lazerson (lazerson@pppl.gov)
!     Date:          03/02/2012
!     Description:   This subroutine calculates the flux due to a line
!                    segment using midpoint formula.
!-----------------------------------------------------------------------
      REAL FUNCTION db_midpoint(x0,y0,z0,x1,y1,z1,cg,ldb)
!-----------------------------------------------------------------------
!     Libraries
!-----------------------------------------------------------------------
      USE stel_kinds, ONLY: rprec
      USE virtual_casing_mod
      USE biotsavart
      USE bsc_T
!-----------------------------------------------------------------------
!     Local Variables
!          ier            Error Flag
!          iunit          File ID Number
!-----------------------------------------------------------------------
      INTEGER, intent(in) :: CG
      real(rprec)   , intent(in) :: x0,y0,z0,x1,y1,z1
      LOGICAL, intent(in), optional :: ldb
      LOGICAL     :: ldb_local
      INTEGER     :: ier
      real(rprec) :: dx,dy,dz,xp,yp,zp,bx,by,bz
      REAL(rprec) :: xvec(3),bvec(3)
!-----------------------------------------------------------------------
!     Begin Function
!-----------------------------------------------------------------------
      dx = x1-x0
      dy = y1-y0
      dz = z1-z0
      xp = .5*(x0+x1)
      yp = .5*(y0+y1)
      zp = .5*(z0+z1)
      dl = 1.0
      ldb_local=.false.
      IF (PRESENT(ldb)) ldb_local = ldb
      IF (cg > 0) THEN
         xvec(1)=xp
         xvec(2)=yp
         xvec(3)=zp
         CALL bsc_b(coil_group(cg),xvec,bvec)
         dflux_midpoint = bvec(1)*dx+bvec(2)*dy+bvec(3)*dz
      ELSE
         ier = 1
         CALL bfield_vc(xp,yp,zp,bx,by,bz,ier)
         !CALL bfield_virtual_casing_adapt(xp,yp,zp,bx,by,bz,ier)
         !CALL bfield_volint_adapt(xp,yp,zp,bx,by,bz,ier)
         db_midpoint = bx*dx+by*dy+bz*dz
      END IF
      IF (ldb_local) db_midpoint = db_midpoint/SQRT(dx*dx+dy*dy+dz*dz)
      RETURN
!-----------------------------------------------------------------------
!     End Function
!-----------------------------------------------------------------------
      END FUNCTION db_midpoint

max_stars_repo_name: joseluisvelasco/STELLOPT

max_issues_repo_name: joseluisvelasco/STELLOPT

max_forks_repo_name: joseluisvelasco/STELLOPT


! PR fortran/49103
! { dg-do run }
  integer :: a(2), b(2), i, j
  open (10, status='scratch')
  do j = 1, 2
    a = (/ 0, 0 /)
    b = (/ 1, 1 /)
    do i = 1, 2
      write (10, *) a
      write (10, *) b
    end do
  end do
  rewind (10)
  do i = 0, 7
    read (10, *) a
    if (any (a .ne. mod (i, 2))) STOP 1
  end do
  close (10)
end


submodule (parent) container
  implicit none

contains

    module procedure init
      p%mother = mother
      p%father = father
    end procedure init

    module subroutine harmonize(p)
      type(parent_type), intent(inout) :: p
      real :: avg

      avg = 0.5 * (p%father + p%mother)
      p%father = avg
      p%mother = avg
    end subroutine harmonize

    module function parent_weight(p) result(w)
      type(parent_type), intent(in) :: p
      real :: w

      w = p%mother**2 + p%father**2
    end function parent_weight

end submodule container


text:
      SUBROUTINE DPOFA(A,LDA,N,INFO)
      INTEGER LDA,N,INFO
      DOUBLE PRECISION A(LDA,1)
C
C     DPOFA FACTORS A DOUBLE PRECISION SYMMETRIC POSITIVE DEFINITE
C     MATRIX.
C
C     DPOFA IS USUALLY CALLED BY DPOCO, BUT IT CAN BE CALLED
C     DIRECTLY WITH A SAVING IN TIME IF  RCOND  IS NOT NEEDED.
C     (TIME FOR DPOCO) = (1 + 18/N)*(TIME FOR DPOFA) .
C
C     ON ENTRY
C
C        A       DOUBLE PRECISION(LDA, N)
C                THE SYMMETRIC MATRIX TO BE FACTORED.  ONLY THE
C                DIAGONAL AND UPPER TRIANGLE ARE USED.
C
C        LDA     INTEGER
C                THE LEADING DIMENSION OF THE ARRAY  A .
C
C        N       INTEGER
C                THE ORDER OF THE MATRIX  A .
C
C     ON RETURN
C
C        A       AN UPPER TRIANGULAR MATRIX  R  SO THAT  A = TRANS(R)*R
C                WHERE  TRANS(R)  IS THE TRANSPOSE.
C                THE STRICT LOWER TRIANGLE IS UNALTERED.
C                IF  INFO .NE. 0 , THE FACTORIZATION IS NOT COMPLETE.
C
C        INFO    INTEGER
C                = 0  FOR NORMAL RETURN.
C                = K  SIGNALS AN ERROR CONDITION.  THE LEADING MINOR
C                     OF ORDER  K  IS NOT POSITIVE DEFINITE.
C
C     LINPACK.  THIS VERSION DATED 08/14/78 .
C     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
C
C     SUBROUTINES AND FUNCTIONS
C
C     BLAS DDOT
C     FORTRAN DSQRT
C
C     INTERNAL VARIABLES
C
      DOUBLE PRECISION DDOT,T
      DOUBLE PRECISION S
      INTEGER J,JM1,K
C     BEGIN BLOCK WITH ...EXITS TO 40
C
C
         DO 30 J = 1, N
            INFO = J
            S = 0.0D0
            JM1 = J - 1
            IF (JM1 .LT. 1) GO TO 20
            DO 10 K = 1, JM1
               T = A(K,J) - DDOT(K-1,A(1,K),1,A(1,J),1)
               T = T/A(K,K)
               A(K,J) = T
               S = S + T*T
   10       CONTINUE
   20       CONTINUE
            S = A(J,J) - S
C     ......EXIT
            IF (S .LE. 0.0D0) GO TO 40
            A(J,J) = DSQRT(S)
   30    CONTINUE
         INFO = 0
   40 CONTINUE
      RETURN
      END


text:
program main
    use omp_lib
    implicit none

    integer, parameter :: N = 200
    double precision, dimension(N) :: a, b, c
    double precision :: output(N, N, N)
    integer :: i, j, k
    integer :: bounds(6)
    double precision :: t
    integer :: P(2, 2)

    do i = 1, N
        a(i) = i - 1
        b(i) = i - 1
        c(i) = i - 1
    end do

    bounds = (/ 1, 2, 1, 3, 2, 3 /)
    !output(:, :, :) = spread(spread(a(1, 1, :), 1, N) * spread(b(1, :, 1), 2, N), 1, N) &
    !                          * spread(spread(c(:, 1, 1), 2, N), 3, N)
    !output(:, :, :) = spread(spread(a, 1, N), 2, N) * spread(spread(b, 1, N), 3, N) &
    !                          * spread(spread(c, 2, N), 3, N)
    !output(:, :, :) = spread(spread(a, 1, N) * spread(b, 2, N), 1, N) &
    !                          * spread(spread(c, 2, N), 3, N)
    !do i = 1, N
    !    do j = 1, N
    !        do k = 1, N
    !            if (output(i, j, k) .ne. (i-1) * (j-1) * (k-1)) then
    !                print *, "not equal", i-1, j-1, k-1, output(i, j, k), &
    !                                      (i-1) * (j-1) * (k-1)
    !                stop
    !            end if
    !        end do
    !    end do
    !end do

    !do i = 1, 100
    !    output(:, :, :) = spread(spread(a, 1, N) * spread(b, 2, N), 1, N) &
    !                          * spread(spread(c, 2, N), 3, N)
    !end do

    t = omp_get_wtime()
    do i = 1, 10
        output(:, :, :) = spread(spread(a, bound(1), N), bound(2), N) * &
                                 spread(spread(b, bound(3), N), bound(4), N) * &
                                 spread(spread(c, bound(5), N), bound(6), N)
    end do
    t = omp_get_wtime() - t
    print *, t

contains
    integer function bound(i)
        implicit  none
        integer :: i
        bound = bounds(i)
    end function
end program main

max_stars_repo_name: markflorisson/minivect

max_issues_repo_name: markflorisson/minivect

max_forks_repo_name: markflorisson/minivect


text:
      subroutine EDGETOCELL(
     &           edgeData
     &           ,iedgeDatalo0
     &           ,iedgeDatahi0
     &           ,cellData
     &           ,icellDatalo0
     &           ,icellDatahi0
     &           ,icellBoxlo0
     &           ,icellBoxhi0
     &           ,dir
     &           )
      implicit none
      integer*8 ch_flops
      COMMON/ch_timer/ ch_flops
      integer CHF_ID(0:5,0:5)
      data CHF_ID/ 1,0,0,0,0,0 ,0,1,0,0,0,0 ,0,0,1,0,0,0 ,0,0,0,1,0,0 ,0
     &,0,0,0,1,0 ,0,0,0,0,0,1 /
      integer iedgeDatalo0
      integer iedgeDatahi0
      REAL*8 edgeData(
     &           iedgeDatalo0:iedgeDatahi0)
      integer icellDatalo0
      integer icellDatahi0
      REAL*8 cellData(
     &           icellDatalo0:icellDatahi0)
      integer icellBoxlo0
      integer icellBoxhi0
      integer dir
      integer i
      integer ii
      do i = icellBoxlo0,icellBoxhi0
      ii = i+CHF_ID(0,dir)
      cellData(i) = (0.500d0)*(
     &                    edgeData(i)
     &                   +edgeData(ii))
      enddo
      return
      end
      subroutine EDGETOINCREMENTCELL(
     &           edgeData
     &           ,iedgeDatalo0
     &           ,iedgeDatahi0
     &           ,cellData
     &           ,icellDatalo0
     &           ,icellDatahi0
     &           ,icellBoxlo0
     &           ,icellBoxhi0
     &           ,dir
     &           )
      implicit none
      integer*8 ch_flops
      COMMON/ch_timer/ ch_flops
      integer CHF_ID(0:5,0:5)
      data CHF_ID/ 1,0,0,0,0,0 ,0,1,0,0,0,0 ,0,0,1,0,0,0 ,0,0,0,1,0,0 ,0
     &,0,0,0,1,0 ,0,0,0,0,0,1 /
      integer iedgeDatalo0
      integer iedgeDatahi0
      REAL*8 edgeData(
     &           iedgeDatalo0:iedgeDatahi0)
      integer icellDatalo0
      integer icellDatahi0
      REAL*8 cellData(
     &           icellDatalo0:icellDatahi0)
      integer icellBoxlo0
      integer icellBoxhi0
      integer dir
      integer i0
      integer ii0
      ii0=CHF_ID(0, dir)
      do i0 = icellBoxlo0,icellBoxhi0
         cellData(i0) = cellData(i0) + (0.500d0)*(
     &      edgeData(i0) + edgeData(i0+ii0))
      enddo
      return
      end
      subroutine EDGETOCELLMAX(
     &           edgeData
     &           ,iedgeDatalo0
     &           ,iedgeDatahi0
     &           ,cellData
     &           ,icellDatalo0
     &           ,icellDatahi0
     &           ,icellBoxlo0
     &           ,icellBoxhi0
     &           ,dir
     &           )
      implicit none
      integer*8 ch_flops
      COMMON/ch_timer/ ch_flops
      integer CHF_ID(0:5,0:5)
      data CHF_ID/ 1,0,0,0,0,0 ,0,1,0,0,0,0 ,0,0,1,0,0,0 ,0,0,0,1,0,0 ,0
     &,0,0,0,1,0 ,0,0,0,0,0,1 /
      integer iedgeDatalo0
      integer iedgeDatahi0
      REAL*8 edgeData(
     &           iedgeDatalo0:iedgeDatahi0)
      integer icellDatalo0
      integer icellDatahi0
      REAL*8 cellData(
     &           icellDatalo0:icellDatahi0)
      integer icellBoxlo0
      integer icellBoxhi0
      integer dir
      integer i
      integer ii
      do i = icellBoxlo0,icellBoxhi0
      ii = i+CHF_ID(0,dir)
      cellData(i) = max(
     &                    edgeData(i),
     &                    edgeData(ii))
      enddo
      return
      end


text:
Prostitution is the act of accepting money in exchange for engaging in sex sexual activities. It is currently illegal in Davis and the rest of California. However, illicit activities do take place in Davis.  Some students pay for their Move to D1 steep tuition costs, rent and other living expenses by providing sexual services for money, primarily to people from out of town.

Some contemporary feminists see little difference between the role a wife plays and that of a sex worker paid for hire.

History
Prostitution has a long and noble legacy, one that began hundreds of thousands of years ago among our prehistoric ancestors. Indeed, sex was not only one of the oldest commodities, it is also the profession on which both monogamous marriage and human societies were originally founded.

Where to find it
The nearest locality that has legal prostitution is a mere few hours away by car near Reno, Nevada.  The Old Bridge Ranch and Wild Horse Canyon Resort are just over the border, and several other legal brothels lie further on in the state. Prostitution is legal in all Nevada counties excluding Clark County including Las Vegas, Washoe County including Reno, Carson City, Douglas County, Lincoln County and Eureka County.  This basically means, all the inhabited places.

The nearest place to reliably find local ladies of the night is http://www.sacbee.com/content/politics/local_government/story/13386995p14228570c.html Stockton Blvd, Broadway, and along the American and Sacramento rivers in Sacramento.

Sacramento City Council is currently debating whether pictures of convicted johns and male (but not female) prostitutes should be put onto a city website to shame them, and they will make their decision in September.
   Does anyone know what happened with this?  Its allegedly past September now.

http://sfbay.craigslist.org/ Craigslist is also about as good a place as any to find a prostitute.

Qualms and Debates
The laws surrounding prostitution in California creates a moral dilemma for some, but one that is more open to debate in the Davis Politics political climate of Davis.  Some feel that paying for sex cheapens the whole experience of sex, often removing even the semblance of love from the encounter, however, others make the argument that it is actually an empowering experience for the sex worker, as the sex worker is actually receiving a fair payment for services performed  in contrast to some people who put out all the time without receiving getting anything in return (a common scenario of unfulfilling, longstanding relationships, such as marriage).  These beliefs and qualms about exchanging money for physical pleasure are influenced by social norms (many of which are directly imported from earlier, more misogynistic days) and Spiritual Organizations religious beliefs.

  We should have a hotline for anonymous counseling about harm reduction and alternatives for people considering selling themselves in some way. No one should have to go this route because of economic necessity.

    I am available to provide counseling  feel free to email me at MailTo(starbaby5 AT hotmail DOT com) and we can chat and I can help you find another way out if you would like that, or just give you support and a place to vent. I graduated awhile ago and no longer live in Davis, so nothing you say to me will get around anywhere.   Users/CristinaDeptula

      Beyond that, Cristina sounds a bit sexist and homophobic, considering she ruled out male prostitutes in their entirety, not to mention lesbian prostitutes. What, Christina, theyre not good enough to be considered part of a moral dilemma? Or is that stacking a moral dilemma on a moral dilemma?  SS

Bonuses

If youre lonely, dont know many people in town, kind of a loner but need some good (or maybe just some) sex, prostitution should be a pretty consistent bet.

Other Drawbacks

Prostitutes may have STDs (although anyone else you sleep with may as well). This is mainly an issue when sex workers are not adequately screened for STDs and STIs, in an unregulated environment such as California. Many sex workers regularly get tested.

Local Prostitution References
Prostitution sometimes occurs in Davis. Here are local stories and allegations of prostitution in Davis.

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. (WWWlink)

The day before Valentines Day, 2006, the following ad was spotted on Craigslist:

Attractive Indian Male college student for hire!  21

I am a 21 year old Indian male attending UC Davis. I am currently experiencing financial difficulties and have come to this as a last resort. I would be more than happy to accompany a woman and ensure her complete and total satisfaction. All I ask is $60/30 minutes, $100/hour, and $80 for each additional hour up to a total of 8 hours.

Keep in mind that you are paying for my time and my time only. Anything else is purely between 2 consentual adults.
I look forward to hearing from you!

Attached were two pictures, presumably enough for enterprising students to figure out the identity and shame the enterprising strumpet in public, 18thcentury style.

Allegations of Prostitution in Local Massage Parlors
  On January 22, 2009, the owner of Spring 380 was arrested on charges of prostitution. An article from the Davis Enterprise reporting the event is on Massage Services.
  It only took the Bee nine days to pick up on this, they make references to three (3) other potential massage parlors in town.