text:
!***********************************************************************
subroutine dtm5(day,f,fbar,ap60,alti,hl,alat,xlon,tz,tinf,ro,d,wmm)
  !***********************************************************************
  ! CNES DTM2020 Research: F30 and Hp60/ap60
  ! ver 21/10/2020
  ! calculation of temperature and density with DTM2020_Res
  !
  !*par ** INPUT **
  !     day     = day of year [1-366]
  !     f       = f(1)=instantaneous flux at (t - 24hr)   /  f(2)=0.
  !     fbar    = fbar(1)=mean flux of last 81 days at t  /  fbar(2)=0.
  !     ap60(1) = 4hr delayed ap60 at t
  !     ap60(2) = 0hr delayed ap60 at t
  !     ap60(3) = 1hr delayed ap60 at t
  !     ap60(4) = 2hr delayed ap60 at t
  !     ap60(5) = 3hr delayed ap60 at t
  !     ap60(6) = mean of last 24 hours
  !     ap60(7) = mean of 5-6-7hr delayed at t
  !     ap60(8) = mean of 9-10-11hr delayed at t
  !     ap60(9) = mean of 14-15-16hr delayed at t
  !     ap60(10)= mean of 19-20-21hr delayed at t
  !
  !     alti   = altitude (in km) greater than 120 km
  !     hl     = local time (in radian: 0-24hr = 0-2pi)
  !     alat   = latitude (in radian)
  !     xlon   = longitude (in radian)
  !
  !*par ** OUTPUT **
  !     tz      = temperature at altitude -> alti
  !     tinf    = exospheric temperature
  !     d(1)    = partial density of atomic hydrogen (in gram/cm3)
  !     d(2)    = partial density of helium
  !     d(3)    = partial density of atomic oxygen
  !     d(4)    = partial density of molecular nitrogen
  !     d(5)    = partial density of molecular oxygen
  !     d(6)    = partial density of atomic nitrogen
  !     ro      = total density (in gram/cm3)
  !     wmm     = mean molecular mass (in gram)
  !***********************************************************************
  !
  !.. Implicit Declarations ..
  implicit none
  !
  !.. Parameters ..
  integer, parameter :: nlatm = 96
  !
  !.. Formal Arguments ..
  real, intent(in) :: alat,alti,day,hl,xlon
  real, intent(out) :: ro,tinf,tz,wmm
  real, dimension(2), intent(in) :: f,fbar
  real, dimension(7), intent(in) :: ap60
  real, dimension(6), intent(out) :: d
  !
  !.. Local Scalars ..
  integer :: i,n,latabs
  integer :: ityp = 0
  real :: asc,c,c2,c4,s,s2,dt120,dtinf, &
          dtp120,gamma,gdelaz,gdelh,gdelo,gdelt0,gdeltp,tinftz,ts,zeta
  real :: cpmg = .19081,re = 6356.77,rgas = 831.4
  real :: dzeta,dzeta2,expsz,gdelaz2,gdelhe,gdelo2,gdelt,glb,secl2,sigma, &
          sigzeta,t,t120sd,t120tt,t120tz,upapg,zlb
  real :: zlb0 = 120.,xlmg = -1.2392
  real :: t120 = 0.0,tp120=0.0,xlog = 0.0
  real :: cose = .9175,gsurf = 980.665,sine = .3978,spmg = .98163,zero = 0.
  real :: crd,dlat,dlon,gmlatd,gmlond,gmlatr,gmlonr,xl75,xl45
  !
  !.. Local Arrays ..
  integer, dimension(6) :: ma
  real, dimension(6) :: alefa,vma
  real, dimension(6) :: cc = 0.0,dbase = 0.0,fz = 0.0
  real, dimension(8) :: akp
  real, dimension(nlatm) :: daz,daz2,dh,dhe,do2,dox,dt0,dtp,dtt
  !
  !.. External Functions ..
  real, external :: bint_oe   !NB: Open Ended ap60->Hp60
  !
  !.. Common Blocks ..
  common /cons/ pi,deupi,cdr,sard
  real :: cdr,deupi,pi,sard
  common /hlocal/ hl0,ch,sh,c2h,s2h,c3h,s3h
  real :: c2h,c3h,ch,hl0,s2h,s3h,sh
  common /pardtm/ tt,h,he,o,az2,o2,az,t0,tp
  real, dimension(nlatm) :: az,az2,h,he,o,o2,t0,tp,tt
  common /plgdtm/ p10,p20,p30,p40,p50,p60,p11,p21,p31,p41,p51,p22,p32,p42,p52, &
                  p62,p33,p10mg,p20mg,p40mg,p11mg,p22mg,p31mg,p30mg,p50mg,p60mg
  real :: p10,p10mg,p11,p20,p20mg,p21,p22,p30,p31,p32,p33,p40,p40mg,p41,p42, &
          p50,p51,p52,p60,p62,p11mg,p22mg,p31mg,p30mg,p50mg,p60mg
  !
  !.. Data Declarations ..
  data alefa/-0.40,-0.38,0.,0.,0.,0./
  data ma/1,4,16,28,32,14/
  data vma/1.6606e-24,6.6423e-24,26.569e-24,46.4958e-24,53.1381e-24,23.2479e-24/
  !
  ! ... Executable Statements ...
  !
  zlb = zlb0
  !
  ro = 0.
  dtinf = 0.
  dt120 = 0.
  dtp120 = 0.
  fz(1) = 0.
  fz(2) = 0.
  fz(3) = 0.
  fz(4) = 0.
  fz(5) = 0.
  fz(6) = 0.
  !
  do i = 1,nlatm
    dtt(i) = 0.
    dh(i) = 0.
    dhe(i) = 0.
    dox(i) = 0.
    daz2(i) = 0.
    do2(i) = 0.
    daz(i) = 0.
    dt0(i) = 0.
    dtp(i) = 0.
  end do
  !
  c = sin(alat)
  c2 = c * c
  c4 = c2 * c2
  s = cos(alat)
  s2 = s * s
  p10 = c
  p20 = 1.5*c2 - 0.5
  p30 = c * (2.5*c2-1.5)
  p40 = 4.375*c4 - 3.75*c2 + 0.375
  p50 = c * (7.875*c4-8.75*c2+1.875)
  p60 = (5.5*c*p50-2.5*p40) / 3.
  p11 = s
  p21 = 3. * c * s
  p31 = s * (7.5*c2-1.5)
  p41 = c * s * (17.5*c2-7.5)
  p51 = s * (39.375*c4-26.25*c2+1.875)
  p22 = 3. * s2
  p32 = 15. * c * s2
  p42 = s2 * (52.5*c2-7.5)
  p52 = 3.*c*p42 - 2.*p32
  p62 = 2.75*c*p52 - 1.75*p42
  p33 = 15. * s * s2
  !
    crd=180./pi
    dlat=alat*crd
    dlon=xlon*crd
    call geogm(dlat,dlon,gmlatd,gmlond)
    gmlatr=gmlatd*cdr
    gmlonr=gmlond*cdr
    c = sin(gmlatr)
    c2 = c * c
    c4 = c2 * c2
    s = cos(gmlatr)
    s2 = s * s
    p10mg = c
    p20mg = 1.5*c2 - 0.5
    p30mg = c * (2.5*c2-1.5)
    p40mg = 4.375*c4 - 3.75*c2 + 0.375
    p50mg = c * (7.875*c4-8.75*c2+1.875)
    p60mg = (5.5*c*p50mg-2.5*p40mg) / 3.
    p11mg = s
    p22mg = s2*3.
    p31mg = s * (7.5*c2-1.5)
  !
  hl0 = hl
  ch = cos(hl0)
  sh = sin(hl0)
  c2h = ch*ch - sh*sh
  s2h = 2. * ch * sh
  c3h = c2h*ch - s2h*sh
  s3h = s2h*ch + c2h*sh
!
!
  latabs = int(abs(dlat))
  select case (latabs)
  case (70 : 90)
    xl75 = (90.-float(latabs))/20.
    akp(1) = bint_oe( (1-xl75)*((ap60(3)+ap60(2)+ap60(4))/3.) + xl75*((ap60(3)+ap60(4)+ap60(5))/3.) )
  case (30 : 69)
    xl45 = (69.-float(latabs))/40.
    akp(1) = bint_oe( (1-xl45)*((ap60(3)+ap60(4)+ap60(5))/3.) + xl45*((ap60(5)+ap60(4)+ap60(1))/3.) )
  case (0 : 29)
    akp(1) = bint_oe( (ap60(4)+ap60(5)+ap60(1))/3. )
  end select
  !
  !
  akp(2) = bint_oe(ap60(2)) - bint_oe(ap60(3))
  akp(3) = bint_oe(ap60(6))
  akp(4) = 0.0
  akp(5) = bint_oe(ap60(7))
  akp(6) = bint_oe(ap60(8))
  akp(7) = bint_oe(ap60(9))
  akp(8) = bint_oe(ap60(10))
  !
  !
  call gldtm_Hp(f,fbar,akp,day,tt,dtt,gdelt,1.,xlon)
  dtt(1) = 1. + gdelt
  tinf = tt(1) * dtt(1)
  call gldtm_Hp(f,fbar,akp,day,t0,dt0,gdelt0,1.,xlon)
  dt0(1) = 1. + gdelt0
  t120 = t0(1) * dt0(1)
  call gldtm_Hp(f,fbar,akp,day,tp,dtp,gdeltp,1.,xlon)
  dtp(1) = 1. + gdeltp
  tp120 = tp(1) * dtp(1)
  !-----------------------------------------------------------------------------
  sigma = tp120 / (tinf-t120)
  dzeta = (re+zlb) / (re+alti)
  zeta = (alti-zlb) * dzeta
  dzeta2 = dzeta * dzeta
  sigzeta = sigma * zeta
  expsz = exp(-sigzeta)
  tz = tinf - (tinf-t120)*expsz
  !
    call gldtm_Hp(f,fbar,akp,day,h,dh,gdelh,0.,xlon)
    dh(1)=exp(gdelh)
    dbase(1)=h(1)*dh(1)
    !
    call gldtm_Hp(f,fbar,akp,day,he,dhe,gdelhe,0.,xlon)
    dhe(1)=exp(gdelhe)
    dbase(2)=he(1)*dhe(1)
    !
    call gldtm_Hp(f,fbar,akp,day,o,dox,gdelo,1.,xlon)
    dox(1)=exp(gdelo)
    dbase(3)=o(1)*dox(1)
    !
    call gldtm_Hp(f,fbar,akp,day,az2,daz2,gdelaz2,1.,xlon)
    daz2(1)=exp(gdelaz2)
    dbase(4)=az2(1)*daz2(1)
    !
    call gldtm_Hp(f,fbar,akp,day,o2,do2,gdelo2,1.,xlon)
    do2(1)=exp(gdelo2)
    dbase(5)=o2(1)*do2(1)
    !
    call gldtm_Hp(f,fbar,akp,day,az,daz,gdelaz,1.,xlon)
    daz(1)=exp(gdelaz)
    dbase(6)=az(1)*daz(1)
    !
    glb = gsurf / (1.+zlb/re)**2
    glb = glb / (sigma*rgas*tinf)
    t120tz = t120 / tz
    xlog = LOG(t120tz)
    tinftz = tinf / tz
    t120tt = t120 / (tinf-t120)
!
    do i = 1,6
      gamma = ma(i) * glb
      upapg = 1. + alefa(i) + gamma
      fz(i) = t120tz**upapg * exp(-sigzeta*gamma)
      cc(i) = dbase(i) * fz(i)
      d(i) = cc(i) * vma(i)
      ro = ro + d(i)
    end do
      wmm=ro/(vma(1)*(cc(1)+cc(2)+cc(3)+cc(4)+cc(5)+cc(6)))
!
end subroutine dtm5
!
!***********************************************************************
subroutine gldtm_Hp(f,fbar,akp,day,a,da,gdel,ff0,xlon)
  !***********************************************************************
  ! SB 21/10/2020
  ! Corrections et nouveaux inconnus
  !*rol calcul de la fonction g(l) evoluee pour dtm_2020&Hp60
  !
  !     gdel=resultat du calcul de g(l)
  !***********************************************************************
  !
  !.. Implicit Declarations ..
  implicit none
  !
  !.. Parameters ..
  integer, parameter :: nlatm = 96
  !
  !.. Formal Arguments ..
  real, intent(in) :: day,ff0,xlon
  real, intent(out) :: gdel
  real, dimension(2), intent(in) :: f,fbar
  real, dimension(8), intent(inout) :: akp
  real, dimension(nlatm), intent(in) :: a
  real, dimension(nlatm), intent(inout) :: da
  !
  !.. Local Scalars ..
  integer :: i
  integer :: ikp = 0
  integer :: ikpm,iflux
  real :: a74,a77,a78,a88,f1f,fp,flux,f0,fp1
  real :: c2fi = 0.0
  real :: clfl,cos2te,coste,dakp,dakpm,dkp,dkpm,rsin2te,slfl,rsinte
  real :: rot = .017214206,rot2 = .034428412,roth = .261799387, &
          rots = 7.27220e-05
  !
  !.. Local Arrays ..
  real, dimension(2) :: fbm150 = 0.0,fmfb = 0.0
  !
  common /hlocal/ hl,ch,sh,c2h,s2h,c3h,s3h
  real :: c2h,c3h,ch,hl,s2h,s3h,sh
  common /plgdtm/ p10,p20,p30,p40,p50,p60,p11,p21,p31,p41,p51,p22,p32,p42,p52, &
                  p62,p33,p10mg,p20mg,p40mg,p11mg,p22mg,p31mg,p30mg,p50mg,p60mg
  real :: p10,p10mg,p11,p20,p20mg,p21,p22,p30,p31,p32,p33,p40,p40mg,p41,p42, &
          p50,p51,p52,p60,p62,p11mg,p22mg,p31mg,p30mg,p50mg,p60mg
  !
  !
  da(2) = p20
  da(3) = p40
  da(74) = p10
  da(77) = p30
  da(78) = p50
  da(79) = p60
  !
  fmfb(1) = f(1) - fbar(1)
  fmfb(2) = f(2) - fbar(2)
  fbm150(1) = fbar(1) - 150.
  fbm150(2) = fbar(2)
  da(4) = fmfb(1)
  da(6) = fbm150(1)
  da(5) = da(4) * da(4)
  da(69) = da(6) * da(6)
  da(82) = da(4) * p10
  da(83) = da(4) * p20
  da(84) = da(4) * p30
  da(85) = da(6) * p20
  da(86) = da(6) * p30
  da(87) = da(6) * p40
  !
  if(akp(1).ge.9.0 .and. akp(3).gt.7.5) akp(1)=9.0 + (akp(1)-9.0)/5.
  if(akp(4).ge.9.0 .and. akp(3).gt.7.5) akp(4)=9.0 + (akp(4)-9.0)/5.
  if(akp(5).ge.9.0 .and. akp(3).gt.7.5) akp(5)=9.0 + (akp(5)-9.0)/5.
  if(akp(6).ge.9.0 .and. akp(3).gt.7.5) akp(6)=9.0 + (akp(6)-9.0)/5.
  if(akp(7).ge.9.0 .and. akp(3).gt.7.5) akp(7)=9.0 + (akp(7)-9.0)/5.
  if(akp(8).ge.9.0 .and. akp(3).gt.7.5) akp(8)=9.0 + (akp(8)-9.0)/5.
!
  dkp  = akp(1)
  dkpm = akp(3)
  da(7) = dkp
  da(8) = p20mg * dkp
  da(60) = dkp * dkp
  da(61) = p20mg * da(60)
  da(62) = p30mg * dkp
  da(63) = p10mg * dkp
  da(67) = p60mg * dkp
  da(68) = p40mg * dkp
!
  da(64) = dkpm
  da(65) = p20mg * dkpm
  da(66) = dkpm  * dkpm
  da(73) = p20mg * da(66)
!
  flux=max(f(1),fbar(1))
  iflux=int(flux)
  select case (iflux)
  case (200:)
   da(75) = 0.333*da(60) * da(60)
   da(76) = 0.1*akp(7)*akp(7)*akp(7)*akp(7)
   da(79) = 0.1*akp(8)*akp(8)*akp(8)*akp(8)
   da(71) = 0.1*akp(5)*akp(5)*akp(5)*akp(5)
   da(72) = 0.1*akp(6)*akp(6)*akp(6)*akp(6)
  case (190:199)
   da(75) = 0.55*da(60) * da(60)
   da(76) = 0.15*akp(7)*akp(7)*akp(7)*akp(7)
   da(79) = 0.15*akp(8)*akp(8)*akp(8)*akp(8)
   da(71) = 0.15*akp(5)*akp(5)*akp(5)*akp(5)
   da(72) = 0.15*akp(6)*akp(6)*akp(6)*akp(6)
  case (180:189)
   da(75) = 0.733*da(60) * da(60)
   da(76) = 0.2*akp(7)*akp(7)*akp(7)*akp(7)
   da(79) = 0.2*akp(8)*akp(8)*akp(8)*akp(8)
   da(71) = 0.2*akp(5)*akp(5)*akp(5)*akp(5)
   da(72) = 0.2*akp(6)*akp(6)*akp(6)*akp(6)
  case (160:179)
   da(75) = da(60) * da(60)
   da(76) = 0.4*akp(7)*akp(7)*akp(7)*akp(7)
   da(79) = 0.4*akp(8)*akp(8)*akp(8)*akp(8)
   da(71) = 0.4*akp(5)*akp(5)*akp(5)*akp(5)
   da(72) = 0.4*akp(6)*akp(6)*akp(6)*akp(6)
  case (140:159)
   da(75) = da(60) * da(60)
   da(76) = 0.8*akp(7)*akp(7)*akp(7)*akp(7)
   da(79) = 0.8*akp(8)*akp(8)*akp(8)*akp(8)
   da(71) = 0.8*akp(5)*akp(5)*akp(5)*akp(5)
   da(72) = 0.8*akp(6)*akp(6)*akp(6)*akp(6)
  case default
    da(75) = da(60) * da(60)
    da(71) = akp(5)*akp(5)*akp(5)*akp(5)
    da(72) = akp(6)*akp(6)*akp(6)*akp(6)
    da(76) = akp(7)*akp(7)*akp(7)*akp(7)
    da(79) = akp(8)*akp(8)*akp(8)*akp(8)
  end select
  da(70) = akp(2)
!
  f0 = a(4)*da(4) + a(5)*da(5) + a(6)*da(6) + a(69)*da(69) + a(82)*da(82) + &
       a(83)*da(83) + a(84)*da(84) + a(85)*da(85) + a(86)*da(86) + &
       a(87)*da(87)
  f1f = 1. + f0*ff0
  !
  f0 = f0 + a(2)*da(2) + a(3)*da(3) + a(74)*da(74) + a(77)*da(77) + a(7)*da(7) + &
       a(8)*da(8) + a(60)*da(60) + a(61)*da(61) + a(68)*da(68) + a(64)*da(64) &
       + a(65)*da(65) + a(66)*da(66) + a(72)*da(72) + a(73)*da(73) + &
       a(75)*da(75) + a(76)*da(76) + a(78)*da(78) + a(79)*da(79) &
       + a(70)*da(70) + a(71)*da(71) + a(62)*da(62) + a(63)*da(63) + a(67)*da(67)
  !
  da(9) = cos(rot*(day-a(11)))
  da(10) = p20 * da(9)
  da(12) = cos(rot2*(day-a(14)))
  da(13) = p20 * da(12)
  coste = cos(rot*(day-a(18)))
  da(15) = p10 * coste
  da(16) = p30 * coste
  da(17) = da(6) * da(15)
  cos2te = cos(rot2*(day-a(20)))
  da(19) = p10 * cos2te
  da(39) = p30 * cos2te
  da(59) = da(6) * da(19)
  da(21) = p11 * ch
  da(22) = p31 * ch
  da(23) = da(6) * da(21)
  da(24) = da(21) * coste
  da(25) = p21 * ch * coste
  da(26) = p11 * sh
  da(27) = p31 * sh
  da(28) = da(6) * da(26)
  da(29) = da(26) * coste
  da(30) = p21 * sh * coste
  da(94)=p51*ch
  da(95)=p51*sh
  da(31) = p22 * c2h
  da(37) = p42 * c2h
  da(32) = p32 * c2h * coste
  da(33) = p22 * s2h
  da(38) = p42 * s2h
  da(34) = p32 * s2h * coste
  da(88) = p32 * c2h
  da(89) = p32 * s2h
  da(90) = da(6) * da(31)
  da(91) = da(6) * da(33)
  da(92) = p62 * c2h
  da(93) = p62 * s2h
  da(35) = p33 * c3h
  da(36) = p33 * s3h
  fp = a(9)*da(9) + a(10)*da(10) + a(12)*da(12) + a(13)*da(13) + a(15)*da(15) &
       + a(16)*da(16) + a(17)*da(17) + a(19)*da(19) + a(21)*da(21) + &
       a(22)*da(22) + a(23)*da(23) + a(24)*da(24) + a(25)*da(25) + &
       a(26)*da(26) + a(27)*da(27) + a(28)*da(28) + a(29)*da(29) + &
       a(30)*da(30) + a(31)*da(31) + a(32)*da(32) + a(33)*da(33) + &
       a(34)*da(34) + a(35)*da(35) + a(36)*da(36) + a(37)*da(37) + &
       a(38)*da(38) + a(39)*da(39) + a(59)*da(59) + a(88)*da(88) + a(89)*da(89) + &
       a(90)*da(90) + a(91)*da(91) + a(92)*da(92) + a(93)*da(93) + &
       a(94)*da(94) + a(95)*da(95)
  !
  da(40) = p10mg * cos2te * dkpm
  da(41) = p10mg * coste  * dkpm
  da(42) = p10mg * cos2te * dkp
  da(43) = p11mg * ch * dkp
  da(44) = p31mg * ch * dkp
  da(45) = p22mg * c2h * dkp
  da(46) = p11mg * sh * dkp
  da(47) = p31mg * sh * dkp
  da(48) = p22mg * s2h * dkp
  !
  fp = fp + a(40)*da(40) + a(41)*da(41) + a(42)*da(42) + a(43)*da(43) + &
       a(44)*da(44) + a(45)*da(45) + a(46)*da(46) + a(47)*da(47) + &
       a(48)*da(48)
    !
    clfl = cos(xlon)
    da(49) = p11 * clfl
    da(50) = p21 * clfl
    da(51) = p31 * clfl
    da(52) = p41 * clfl
    da(53) = p51 * clfl
    slfl = sin(xlon)
    da(54) = p11 * slfl
    da(55) = p21 * slfl
    da(56) = p31 * slfl
    da(57) = p41 * slfl
    da(58) = p51 * slfl
    !
    fp = fp + a(49)*da(49) + a(50)*da(50) + a(51)*da(51) + a(52)*da(52) + &
         a(53)*da(53) + a(54)*da(54) + a(55)*da(55) + a(56)*da(56) + &
         a(57)*da(57) + a(58)*da(58)
  !
  gdel = f0 + fp*f1f
  !
end subroutine gldtm_Hp
!
!***********************************************************************
subroutine lecdtm
  ! SB 21/10/2020
  !
  !reading in the DTM model coefficients in vectors, per constituent
!************************************************************************
  !.. Parameters ..
  integer, parameter :: nlatm = 96
  !
  !
  character titre*100,ct120*8,ctp120*8
  common /pardtm/ tt,h,he,o,az2,o2,az,t0,tp
  real, dimension(nlatm) :: az,az2,h,he,o,o2,t0,tp,tt
  !
  open(10,file='DTM_2020_F30_ap60',status='old')
  !
  read(iu10,530) titre
  read(iu10,540) npdtm
  write(*,600) titre
  !
  do i=1,npdtm
    read(iu10,540,end=20) ni,tt(i),dtt,h(i),dh,he(i),dhe,o(i),dox, &
            az2(i),daz2,o2(i),do2,az(i),daz,t0(i),dt0,tp(i),dtp
  end do
!
20 continue
!
530  format(a100)
540  format(i4,9(e13.6,e9.2))
600  format(a100/)
!
end subroutine lecdtm
!
!***********************************************************************
SUBROUTINE GEOGM(XLAT,XLONG,GMLAT,GMLONG)
! conversion geographic to geomagnetic latitude&longitude
!
  !.. Implicit Declarations ..
  implicit none
  !.. Formal Arguments ..
  real, intent(in) :: XLAT,XLONG
  real, intent(out) :: GMLAT,GMLONG
  !
  !.. Local Scalars ..
  real :: pi,RFAC,PLATR,PLONGR,SPL,CPL,RLAT,RLONG
  real :: SLM,CLM,PHIM1,PHIM2
  !
  pi=acos(-1.)
  if(XLONG.EQ.291.0000) XLONG=291.1
  RFAC=PI/180.
  PLATR=78.5*RFAC
  PLONGR=291.*RFAC
  SPL=SIN(PLATR)
  CPL=COS(PLATR)
  RLAT=XLAT*RFAC
  RLONG=XLONG*RFAC
  SLM=SPL*SIN(RLAT)+CPL*COS(RLAT)*COS(PLONGR-RLONG)
  CLM=SQRT(1.-SLM**2)
  PHIM1=COS(RLAT)*SIN(RLONG-PLONGR)/CLM
  PHIM2=(SPL*SLM-SIN(RLAT))/(CPL*CLM)
  GMLAT=ASIN(SLM)/RFAC
  IF(PHIM1.GE.0..AND.PHIM2.GE.0.)GO TO 100
  IF(PHIM1.GE.0..AND.PHIM2.LT.0.)GO TO 101
  IF(PHIM1.LT.0..AND.PHIM2.LT.0.)GO TO 102
  IF(PHIM1.LT.0..AND.PHIM2.GE.0.)GO TO 103
  GO TO 110
!
 100 GMLONG=ASIN(PHIM1)/RFAC
      GO TO 104
 101 GMLONG=(PI-ABS(ASIN(PHIM1)))/RFAC
      GO TO 104
 102 GMLONG=(PI+ABS(ASIN(PHIM1)))/RFAC
      GO TO 104
 103 GMLONG=(2.*PI-ABS(ASIN(PHIM1)))/RFAC
 104 RETURN
 110 PRINT*,"GEOGM: angle error"
     STOP
END SUBROUTINE GEOGM
!
!***********************************************************************
function bint_oe(ap)
  ! conversion Open Ended ap > Open Ended Kp
  !
  !.. Implicit Declarations ..
  implicit none
  !
  !.. Formal Arguments ..
  real, intent(in) :: ap
  !
  !.. Function Declaration ..
  real :: bint_oe
  !
  !.. Local Scalars ..
  integer :: i,i1
  !
  !.. Local Arrays ..
  real, dimension(38) :: aap,kp
  !
  !.. Data Declarations ..
  data (aap(i), i = 1,38)/ &
       0.,2.,3.,4.,5.,6.,7.,9.,12.,15.,18.,22.,27.,32.,39.,48.,56.,67., &
       80.,94.,111.,132.,154.,179.,207.,236.,265.,294.,324.,355.,388.,421.0, &
       456.,494.,534.,574.,617.,657./
  data (kp(i), i = 1,38)/ &
       0.,.33,.66,1.,1.33,1.66,2.,2.33,2.66,3.,3.33,3.66,4.,4.33,4.66,5.,5.33,5.66, &
       6.,6.33,6.66,7.,7.33,7.66,8.,8.33,8.66,9.,9.33,9.66,10.,10.33, &
       10.66,11.,11.33,11.66,12.,12.33/
  !
  ! ... Executable Statements ...
  !
  do i = 1,38
    if (aap(i) == ap) then
      bint_oe = kp(i)
      return
    end if
    if (aap(i) > ap) then
      i1 = i - 1
      !
      bint_oe = kp(i1) + ((kp(i)-kp(i1))/(aap(i)-aap(i1)))*(ap-aap(i1))
      return
    end if
  end do
  write (*,*) 'stop in bint_oe',ap
  stop
!
end function bint_oe

max_stars_repo_name: YoucunWang/mcm

max_stars_repo_licenses: DOC Unlicense

max_issues_repo_name: YoucunWang/mcm

max_issues_repo_licenses: DOC Unlicense

max_forks_repo_name: YoucunWang/mcm

max_forks_repo_licenses: DOC Unlicense


text:
!Crown Copyright 2012 AWE.
!
! This file is part of CloverLeaf.
!
! CloverLeaf is free software: you can redistribute it and/or modify it under
! the terms of the GNU General Public License as published by the
! Free Software Foundation, either version 3 of the License, or (at your option)
! any later version.
!
! CloverLeaf is distributed in the hope that it will be useful, but
! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
! FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
! details.
!
! You should have received a copy of the GNU General Public License along with
! CloverLeaf. If not, see http://www.gnu.org/licenses/.

!>  @brief Fortran ideal gas kernel.
!>  @author Wayne Gaudin
!>  @details Calculates the pressure and sound speed for the mesh chunk using
!>  the ideal gas equation of state, with a fixed gamma of 1.4.

MODULE ideal_gas_kernel_module

CONTAINS

SUBROUTINE ideal_gas_kernel(x_min,x_max,y_min,y_max,                &
                            density,                                &
                            energy,                                 &
                            pressure,                               &
                            soundspeed                              )

  IMPLICIT NONE

  INTEGER :: x_min,x_max,y_min,y_max
  REAL(KIND=8), DIMENSION(x_min-2:x_max+2,y_min-2:y_max+2) :: density
  REAL(KIND=8), DIMENSION(x_min-2:x_max+2,y_min-2:y_max+2) :: energy
  REAL(KIND=8), DIMENSION(x_min-2:x_max+2,y_min-2:y_max+2) :: pressure
  REAL(KIND=8), DIMENSION(x_min-2:x_max+2,y_min-2:y_max+2) :: soundspeed

  INTEGER :: j,k

  REAL(KIND=8) :: sound_speed_squared,v,pressurebyenergy,pressurebyvolume

  DO k=y_min,y_max
    DO j=x_min,x_max
      v=1.0_8/density(j,k)
      pressure(j,k)=(1.4_8-1.0_8)*density(j,k)*energy(j,k)
      pressurebyenergy=(1.4_8-1.0_8)*density(j,k)
      pressurebyvolume=-density(j,k)*pressure(j,k)
      sound_speed_squared=v*v*(pressure(j,k)*pressurebyenergy-pressurebyvolume)
      soundspeed(j,k)=SQRT(sound_speed_squared)
    ENDDO
  ENDDO

END SUBROUTINE ideal_gas_kernel

END MODULE ideal_gas_kernel_module

max_stars_repo_name: uwplse/stng

max_issues_repo_name: uwplse/stng

max_forks_repo_name: uwplse/stng


text:
SUBROUTINE remove(sel_position,sel_pool,bflg,stemf,leaff,storef,rootstoref,     &
                & rootfiberf,bcmstandstem,bcmstandleaf,bcmstandstore,           &
                & bcmflatstem,bcmflatleaf,bcmflatstore,bcmrootstorez,           &
                & bcmrootfiberz,bcmbgstemz,bczht,bcdstm,bcgrainf,bchyfg,        &
                & btmstandstem,btmstandleaf,btmstandstore,btmflatstem,          &
                & btmflatleaf,btmflatstore,btmflatrootstore,btmflatrootfiber,   &
                & btmbgstemz,btmbgleafz,btmbgstorez,btmbgrootstorez,            &
                & btmbgrootfiberz,btzht,btdstm,btgrainf,bdmstandstem,           &
                & bdmstandleaf,bdmstandstore,bdmflatstem,bdmflatleaf,           &
                & bdmflatstore,bdmflatrootstore,bdmflatrootfiber,bdmbgstemz,    &
                & bdmbgleafz,bdmbgstorez,bdmbgrootstorez,bdmbgrootfiberz,bdzht, &
                & bddstm,bdgrainf,bdhyfg,nslay,tot_mass_rem,sel_mass_left)
!
IMPLICIT NONE
!
INCLUDE 'p1werm.inc'
!
! Subroutine arguments
!
REAL :: bcdstm,bcgrainf,bcmflatleaf,bcmflatstem,bcmflatstore,bcmstandleaf,      &
      & bcmstandstem,bcmstandstore,bczht,btdstm,btgrainf,btmflatleaf,           &
      & btmflatrootfiber,btmflatrootstore,btmflatstem,btmflatstore,btmstandleaf,&
      & btmstandstem,btmstandstore,btzht,leaff,rootfiberf,rootstoref,           &
      & sel_mass_left,stemf,storef,tot_mass_rem
INTEGER :: bchyfg,bflg,nslay,sel_pool,sel_position
REAL,DIMENSION(mnsz) :: bcmbgstemz,bcmrootfiberz,bcmrootstorez,btmbgleafz,      &
                      & btmbgrootfiberz,btmbgrootstorez,btmbgstemz,btmbgstorez
REAL,DIMENSION(mnbpls) :: bddstm,bdgrainf,bdmflatleaf,bdmflatrootfiber,         &
                        & bdmflatrootstore,bdmflatstem,bdmflatstore,            &
                        & bdmstandleaf,bdmstandstem,bdmstandstore,bdzht
INTEGER,DIMENSION(mnbpls) :: bdhyfg
REAL,DIMENSION(mnsz,mnbpls) :: bdmbgleafz,bdmbgrootfiberz,bdmbgrootstorez,      &
                             & bdmbgstemz,bdmbgstorez
!
! Local variables
!
INTEGER :: idx,idy,tflg
REAL :: pool_temp1,pool_temp2,start_leaf,start_stem,start_store
REAL,DIMENSION(mnsz) :: pool_temp1z,pool_temp2z,pool_temp3z,start_rootfiber,    &
                      & start_rootstore
!
!     the following subroutine arguments are not used: bczht, bdzht,
!     btzht  jcaii  8/08/08
!
!     + + + purpose + + +
!     this subroutine performs the biomass manipulation of removing
!     biomass. the amount of each component removed is determined by
!     the fraction passed into this subroutine for each component.
!     pools are changed in the order: crop, temporary, residue and
!     locations in the order: stand with roots, flat, below ground.
!     consideration is given that if root mass or stem mass is removed,
!     then the leaves and storage portion must become flat in the same
!     pool. removal of stem mass also results in a reduction in stem
!     count. in order to avoid double accounting, proportions of standing
!     are tracked, but the final adjustment and movement is not done
!     until all removals are completed in that pool.

!     possible future enhancements
!     a)  bioflg - selects which age pools will be processed.  probably the
!     same definition as other biomass manipulation process effects use.
!     b)  xxlocflg - selects the individual mass component pools that are
!     being effected (material being removed in this case).  there would
!     likely need to be more than one of these flags, possibly one for each
!     "age" pool.  example settings could be:
!
!     crlocflg (st,yld,flt,bg,rt)	decomp1locflg (st,flt,bg,rt)
!
!     bit             val                    bit              val
!     x    st+yld+flt  0                      x    st+flt      0
!     0    yld*fract   1                      0    -           1
!     1    st*fract    2                      1    st*fract    2
!     2    fl*fract    4                      2    fl*fract    4
!     3    bg*fract    8                      3    bg*fract    8
!     4    rt*fract    16                     4    rt*fract    16
!     5    st*cutht    32                     5    st*cutht    32

!     + + + keywords + + +
!     remove, biomass manipulation

!     + + + common blocks + + +

!     + + + argument declarations + + +

!     + + + argument definitions + + +
!     sel_position - position to which percentages will be applied
!                0 - don't apply to anything
!                1 - apply to standing (and attached roots)
!                2 - apply to flat
!                3 - apply to standing (and attached roots) and flat
!                4 - apply to buried
!                5 - apply to standing (and attached roots) and buried
!                6 - apply to flat and buried
!                7 - apply to standing (and attached roots), flat and buried
!                this corresponds to the bit pattern:
!                msb(buried, flat, standing)lsb

!     sel_pool - pool to which percentages will be applied
!            0 - don't apply to anything
!            1 - apply to crop pool
!            2 - apply to temporary pool
!            3 - apply to crop and temporary pools
!            4 - apply to residue pools
!            5 - apply to crop and residue pools
!            6 - apply to temporary and residue pools
!            7 - apply to crop, temporary and residue pools
!                this corresponds to the bit pattern:
!                msb(residue, temporary, crop)lsb

!     bflg      - flag indicating what to manipulate
!       0 - all standing material is manipulated (both crop and residue)
!       1 - crop
!       2 - 1'st residue pool
!       4 - 2'nd residue pool
!       ....
!       2**n - nth residue pool

!     storef   - fraction of storage (reproductive components) removed (kg/kg)
!     leaff    - fraction of plant leaves removed (kg/kg)
!     stemf    - fraction of plant stems removed (kg/kg)
!     rootstoref - fraction of plant storage root removed (kg/kg)
!     rootfiberf - fraction of plant fibrous root removed (kg/kg)

!     bcmstandstem - crop standing stem mass (kg/m^2)
!     bcmstandleaf - crop standing leaf mass (kg/m^2)
!     bcmstandstore - crop standing storage mass (kg/m^2)
!                    (head with seed, or vegetative head (cabbage, pineapple))

!     bcmflatstem  - crop flat stem mass (kg/m^2)
!     bcmflatleaf  - crop flat leaf mass (kg/m^2)
!     bcmflatstore - crop flat storage mass (kg/m^2)

!     bcmrootstorez - crop root storage mass by soil layer (kg/m^2)
!                   (tubers (potatoes, carrots), extended leaf (onion), seeds (peanuts))
!     bcmrootfiberz - crop root fibrous mass by soil layer (kg/m^2)

!     bcmbgstemz  - crop buried stem mass by layer (kg/m^2)

!     bczht  - crop height (m)
!     bcdstm - number of crop stems per unit area (#/m^2)
!            - it is computed by taking the tillering factor
!              times the plant population density.
!     bcgrainf - internally computed grain fraction of reproductive mass

!     btmstandstem - crop standing stem mass (kg/m^2)
!     btmstandleaf - crop standing leaf mass (kg/m^2)
!     btmstandstore - crop standing storage mass (kg/m^2)
!                    (head with seed, or vegetative head (cabbage, pineapple))

!     btmflatstem  - crop flat stem mass (kg/m^2)
!     btmflatleaf  - crop flat leaf mass (kg/m^2)
!     btmflatstore - crop flat storage mass (kg/m^2)

!     btmflatrootstore - crop flat root storage mass (kg/m^2)
!                   (tubers (potatoes, carrots), extended leaf (onion), seeds (peanuts))
!     btmflatrootfiber - crop flat root fibrous mass (kg/m^2)

!     btmbgflatstemz  - crop buried stem mass by layer (kg/m^2)
!     btmbgflatleafz  - crop buried leaf mass by layer (kg/m^2)
!     btmbgflatstorez - crop buried storage mass by layer (kg/m^2)

!     btmbgrootstorez - crop root storage mass by layer (kg/m^2)
!                   (tubers (potatoes, carrots), extended leaf (onion), seeds (peanuts))
!     btmbgrootfiberz - crop root fibrous mass by layer (kg/m^2)

!     btzht  - crop height (m)
!     btdstm - number of crop stems per unit area (#/m^2)
!            - it is computed by taking the tillering factor
!              times the plant population density.
!     btgrainf - internally computed grain fraction of reproductive mass

!     bdmstandstem  - standing stem mass (kg/m^2)
!     bdmstandleaf  - standing leaf mass (kg/m^2)
!     bdmstandstore - standing storage mass (kg/m^2)

!     bdmflatstem  - flat stem mass (kg/m^2)
!     bdmflatleaf  - flat leaf mass (kg/m^2)
!     bdmflatstore - flat storage mass (kg/m^2)

!     bdmflatrootstore - flat storage root mass (kg/m^2)
!     bdmflatrootfiber - flat fibrous root mass (kg/m^2)

!     bdmbgstemz  - buried stem mass by layer (kg/m^2)
!     bdmbgleafz  - buried leaf mass by layer (kg/m^2)
!     bdmbgstorez - buried (from above ground) storage mass by layer (kg/m^2)

!     bdmbgrootstorez - buried storage root mass by layer (kg/m^2)
!     bdmbgrootfiberz - buried fibrous root mass by layer (kg/m^2)

!     bdzht  - residue height (m)
!     bddstm - number of residue stems per unit area (#/m^2)
!     bdgrainf - internally computed grain fraction of reproductive mass
!     bdhyfg - flag indicating the part of plant to apply the "grain fraction",
!              grf, to when removing that plant part for yield
!         0     grf applied to above ground storage (seeds, reproductive)
!         1     grf times growth stage factor (see growth.for) applied to above ground storage (seeds, reproductive)
!         2     grf applied to all aboveground biomass (forage)
!         3     grf applied to leaf mass (tobacco)
!         4     grf applied to stem mass (sugarcane)
!         5     grf applied to below ground storage mass (potatoes, peanuts)

!     nlay      - number of layer from which below ground biomass is removed
!     tot_mass_rem - mass of material removed by this harvest operation (kg/m^2)
!     sel_mass_left - mass of material left in pools from which mass is removed
!                     by this harvest operation (kg/m^2)

!     + + + accessed common block variable definitions + + +
!     mnbpls        - max number of decomposition pools (currently=3)
!     mnsz          - max number of soil layers

!     + + + parameters + + +

!     + + + local variables + + +

!     + + + local variable definitions + + +
!     pool_temp1 - used to substitute for non existent pools in crop,
!                  where there are no flatrootstore and flatrootfiber pools
!     pool_temp2 -  see pool_temp1
!     idx       - loop variable for soil layers
!     idy       - loop variable for decomp pools
!     tflg      - temporary flag to carry bioflag value if changes to all pools

!     + + + end specifications + + +

      !set tflg bits correctly for "all" pools if bflg=0
IF (bflg.EQ.0) THEN
  tflg = 1                          ! crop pool
  DO idy = 1,mnbpls
     tflg = tflg + 2**idy           ! decomp pools
  END DO
ELSE
  tflg = bflg
END IF

tot_mass_rem = 0.0
sel_mass_left = 0.0

!     assign crop grain fraction and representaive stem diameter values
!     to temporary pool since material may be transferred to the temporary
!     pool without a specific kill operation
btgrainf = bcgrainf

pool_temp1 = 0.0
pool_temp2 = 0.0
DO idx = 1,nslay
  pool_temp1z(idx) = 0.0
  pool_temp2z(idx) = 0.0
  pool_temp3z(idx) = 0.0
END DO

      ! crop pool
IF (btest(sel_pool,0)) THEN
          ! standing and rooted biomass
          ! set starting values
  start_store = bcmstandstore
  start_leaf = bcmstandleaf
  start_stem = bcmstandstem
  DO idx = 1,nslay
     start_rootstore(idx) = bcmrootstorez(idx)
     start_rootfiber(idx) = bcmrootfiberz(idx)
  END DO
  IF (btest(sel_position,0)) CALL rem_stand_pool(stemf,leaff,storef,rootstoref, &
    & rootfiberf,bcmstandstem,bcmstandleaf,bcmstandstore,bcmrootstorez,         &
    & bcmrootfiberz,nslay,bchyfg,bcgrainf,bcdstm,tot_mass_rem,sel_mass_left)
          ! flat biomass
  IF (btest(sel_position,1)) CALL rem_flat_pool(stemf,leaff,storef,rootstoref,  &
    & rootfiberf,bcmflatstem,bcmflatleaf,bcmflatstore,pool_temp1,pool_temp2,    &
    & bchyfg,bcgrainf,tot_mass_rem,sel_mass_left)
          ! buried biomass
              ! standing not done so root removal done here
  IF (btest(sel_position,2).AND..NOT.btest(sel_position,0))                     &
    & CALL rem_bg_pool(stemf,leaff,storef,rootstoref,rootfiberf,bcmbgstemz,     &
    & pool_temp2z,pool_temp3z,bcmrootstorez,bcmrootfiberz,nslay,bchyfg,bcgrainf,&
    & tot_mass_rem,sel_mass_left)
          ! adjust standing pools if supporting stems or roots removed
  CALL adj_stand_pool(start_stem,start_leaf,start_store,start_rootstore,        &
                    & start_rootfiber,bcmstandstem,bcmstandleaf,bcmstandstore,  &
                    & bcmrootstorez,bcmrootfiberz,bcmflatstem,bcmflatleaf,      &
                    & bcmflatstore,bcdstm,nslay)
END IF
      ! temporary pool
IF (btest(sel_pool,1)) THEN
          ! standing and rooted biomass
          ! set starting values
  start_store = btmstandstore
  start_leaf = btmstandleaf
  start_stem = btmstandstem
  DO idx = 1,nslay
     start_rootstore(idx) = btmbgrootstorez(idx)
     start_rootfiber(idx) = btmbgrootfiberz(idx)
  END DO
  IF (btest(sel_position,0)) CALL rem_stand_pool(stemf,leaff,storef,rootstoref, &
    & rootfiberf,btmstandstem,btmstandleaf,btmstandstore,btmbgrootstorez,       &
    & btmbgrootfiberz,nslay,bchyfg,btgrainf,btdstm,tot_mass_rem,sel_mass_left)
          ! flat biomass
  IF (btest(sel_position,1)) CALL rem_flat_pool(stemf,leaff,storef,rootstoref,  &
    & rootfiberf,btmflatstem,btmflatleaf,btmflatstore,btmflatrootstore,         &
    & btmflatrootfiber,bchyfg,btgrainf,tot_mass_rem,sel_mass_left)
          ! buried biomass
  IF (btest(sel_position,2)) THEN
     IF (btest(sel_position,0)) THEN
                  ! root removal already done in standing
        CALL rem_bg_pool(stemf,leaff,storef,rootstoref,rootfiberf,btmbgstemz,   &
                       & btmbgleafz,btmbgstorez,pool_temp1z,pool_temp2z,nslay,  &
                       & bchyfg,btgrainf,tot_mass_rem,sel_mass_left)
     ELSE
                  ! standing not done so do root removal here
        CALL rem_bg_pool(stemf,leaff,storef,rootstoref,rootfiberf,btmbgstemz,   &
                       & btmbgleafz,btmbgstorez,btmbgrootstorez,btmbgrootfiberz,&
                       & nslay,bchyfg,btgrainf,tot_mass_rem,sel_mass_left)
     END IF
  END IF
          ! adjust standing pools if supporting stems or roots removed
  CALL adj_stand_pool(start_stem,start_leaf,start_store,start_rootstore,        &
                    & start_rootfiber,btmstandstem,btmstandleaf,btmstandstore,  &
                    & btmbgrootstorez,btmbgrootfiberz,btmflatstem,btmflatleaf,  &
                    & btmflatstore,btdstm,nslay)
END IF
      ! residue pools
IF (btest(sel_pool,2)) THEN
  DO idy = 1,mnbpls
     IF (btest(tflg,idy)) THEN
          ! standing and rooted biomass
          ! set starting values
        start_store = bdmstandstore(idy)
        start_leaf = bdmstandleaf(idy)
        start_stem = bdmstandstem(idy)
        DO idx = 1,nslay
           start_rootstore(idx) = bdmbgrootstorez(idx,idy)
           start_rootfiber(idx) = bdmbgrootfiberz(idx,idy)
        END DO
        IF (btest(sel_position,0)) CALL rem_stand_pool(stemf,leaff,storef,      &
          & rootstoref,rootfiberf,bdmstandstem(idy),bdmstandleaf(idy),          &
          & bdmstandstore(idy),bdmbgrootstorez(1,idy),bdmbgrootfiberz(1,idy),   &
          & nslay,bdhyfg(idy),bdgrainf(idy),bddstm(idy),tot_mass_rem,           &
          & sel_mass_left)
          ! flat biomass
        IF (btest(sel_position,1)) CALL rem_flat_pool(stemf,leaff,storef,       &
          & rootstoref,rootfiberf,bdmflatstem(idy),bdmflatleaf(idy),            &
          & bdmflatstore(idy),bdmflatrootstore(idy),bdmflatrootfiber(idy),      &
          & bdhyfg(idy),bdgrainf(idy),tot_mass_rem,sel_mass_left)
          ! buried biomass
        IF (btest(sel_position,2)) THEN
           IF (btest(sel_position,0)) THEN
              ! root removal already done in standing
              CALL rem_bg_pool(stemf,leaff,storef,rootstoref,rootfiberf,        &
                             & bdmbgstemz(1,idy),bdmbgleafz(1,idy),             &
                             & bdmbgstorez(1,idy),pool_temp1z,pool_temp2z,nslay,&
                             & bdhyfg(idy),bdgrainf(idy),tot_mass_rem,          &
                             & sel_mass_left)
           ELSE
              ! standing not done so do root removal here
              CALL rem_bg_pool(stemf,leaff,storef,rootstoref,rootfiberf,        &
                             & bdmbgstemz(1,idy),bdmbgleafz(1,idy),             &
                             & bdmbgstorez(1,idy),bdmbgrootstorez(1,idy),       &
                             & bdmbgrootfiberz(1,idy),nslay,bdhyfg(idy),        &
                             & bdgrainf(idy),tot_mass_rem,sel_mass_left)
           END IF
        END IF
          ! adjust standing pools if supporting stems or roots removed
        CALL adj_stand_pool(start_stem,start_leaf,start_store,start_rootstore,  &
                          & start_rootfiber,bdmstandstem(idy),bdmstandleaf(idy),&
                          & bdmstandstore(idy),bdmbgrootstorez(1,idy),          &
                          & bdmbgrootfiberz(1,idy),bdmflatstem(idy),            &
                          & bdmflatleaf(idy),bdmflatstore(idy),bddstm(idy),     &
                          & nslay)
     END IF
  END DO
END IF

      ! check that complete crop failure shows remaining biomass
IF (tot_mass_rem+sel_mass_left.LE.0.0) sel_mass_left = bcmstandstem +           &
  & bcmstandleaf + bcmstandstore + btmstandstem + btmstandleaf + btmstandstore +&
  & btmflatstem + btmflatleaf + btmflatstore
!
END SUBROUTINE remove
!
SUBROUTINE rem_stand_pool(stemf,leaff,storef,rootstoref,rootfiberf,pool_stem,   &
                        & pool_leaf,pool_store,pool_rootstore,pool_rootfiber,   &
                        & nslay,pool_hyfg,pool_grainf,pool_dstm,tot_mass_rem,   &
                        & sel_mass_left)
!
IMPLICIT NONE
!
INCLUDE 'p1werm.inc'
!
! Subroutine arguments
!
REAL :: leaff,pool_dstm,pool_grainf,pool_leaf,pool_stem,pool_store,rootfiberf,  &
      & rootstoref,sel_mass_left,stemf,storef,tot_mass_rem
INTEGER :: nslay,pool_hyfg
REAL,DIMENSION(mnsz) :: pool_rootfiber,pool_rootstore
!
! Local variables
!
INTEGER :: idx,pool_flag
REAL :: rem_frac
!
!     + + + common blocks + + +
!
pool_flag = 0
!
rem_frac = storef
IF (pool_hyfg.LE.2) rem_frac = rem_frac*pool_grainf
CALL rem_pool(pool_store,rem_frac,pool_flag,tot_mass_rem)

rem_frac = leaff
IF (pool_hyfg.EQ.3) rem_frac = rem_frac*pool_grainf
CALL rem_pool(pool_leaf,rem_frac,pool_flag,tot_mass_rem)

rem_frac = stemf
IF (pool_hyfg.EQ.4) rem_frac = rem_frac*pool_grainf
CALL rem_pool(pool_stem,rem_frac,pool_flag,tot_mass_rem)
      ! also reduce stem count
pool_dstm = pool_dstm*(1.0-rem_frac)

rem_frac = rootstoref
IF (pool_hyfg.EQ.5) rem_frac = rem_frac*pool_grainf
DO idx = 1,nslay
  CALL rem_pool(pool_rootstore(idx),rem_frac,pool_flag,tot_mass_rem)
END DO

rem_frac = rootfiberf
DO idx = 1,nslay
  CALL rem_pool(pool_rootfiber(idx),rem_frac,pool_flag,tot_mass_rem)
END DO

      ! all but fibrous root included in harvest index
IF (pool_flag.EQ.1) THEN
  sel_mass_left = sel_mass_left + pool_store + pool_leaf + pool_stem
  DO idx = 1,nslay
     sel_mass_left = sel_mass_left + pool_rootstore(idx)
  END DO
END IF
!
END SUBROUTINE rem_stand_pool
!
SUBROUTINE rem_flat_pool(stemf,leaff,storef,rootstoref,rootfiberf,pool_stem,    &
                       & pool_leaf,pool_store,pool_rootstore,pool_rootfiber,    &
                       & pool_hyfg,pool_grainf,tot_mass_rem,sel_mass_left)
!
IMPLICIT NONE
!
! Subroutine arguments
!
REAL :: leaff,pool_grainf,pool_leaf,pool_rootfiber,pool_rootstore,pool_stem,    &
      & pool_store,rootfiberf,rootstoref,sel_mass_left,stemf,storef,tot_mass_rem
INTEGER :: pool_hyfg
!
! Local variables
!
INTEGER :: pool_flag
REAL :: rem_frac
!
!     integer idx, pool_flag
!
pool_flag = 0
!
rem_frac = storef
IF (pool_hyfg.LE.2) rem_frac = rem_frac*pool_grainf
CALL rem_pool(pool_store,rem_frac,pool_flag,tot_mass_rem)

rem_frac = leaff
IF (pool_hyfg.EQ.3) rem_frac = rem_frac*pool_grainf
CALL rem_pool(pool_leaf,rem_frac,pool_flag,tot_mass_rem)

rem_frac = stemf
IF (pool_hyfg.EQ.4) rem_frac = rem_frac*pool_grainf
CALL rem_pool(pool_stem,rem_frac,pool_flag,tot_mass_rem)

rem_frac = rootstoref
IF (pool_hyfg.EQ.5) rem_frac = rem_frac*pool_grainf
CALL rem_pool(pool_rootstore,rem_frac,pool_flag,tot_mass_rem)

rem_frac = rootfiberf
CALL rem_pool(pool_rootfiber,rem_frac,pool_flag,tot_mass_rem)

      ! all but fibrous root included in harvest index
IF (pool_flag.EQ.1) sel_mass_left = sel_mass_left + pool_store + pool_leaf +    &
                                  & pool_stem + pool_rootstore

END SUBROUTINE rem_flat_pool
!
SUBROUTINE rem_bg_pool(stemf,leaff,storef,rootstoref,rootfiberf,pool_stem,      &
                     & pool_leaf,pool_store,pool_rootstore,pool_rootfiber,nslay,&
                     & pool_hyfg,pool_grainf,tot_mass_rem,sel_mass_left)
!
IMPLICIT NONE
!
INCLUDE 'p1werm.inc'
!
! Subroutine arguments
!
REAL :: leaff,pool_grainf,rootfiberf,rootstoref,sel_mass_left,stemf,storef,     &
      & tot_mass_rem
INTEGER :: nslay,pool_hyfg
REAL,DIMENSION(mnsz) :: pool_leaf,pool_rootfiber,pool_rootstore,pool_stem,      &
                      & pool_store
!
! Local variables
!
INTEGER :: idx,pool_flag
REAL :: rem_frac
!
!     + + + common blocks + + +
!
pool_flag = 0
!
rem_frac = storef
IF (pool_hyfg.LE.2) rem_frac = rem_frac*pool_grainf
DO idx = 1,nslay
  CALL rem_pool(pool_store(idx),rem_frac,pool_flag,tot_mass_rem)
END DO

rem_frac = leaff
IF (pool_hyfg.EQ.3) rem_frac = rem_frac*pool_grainf
DO idx = 1,nslay
  CALL rem_pool(pool_leaf(idx),rem_frac,pool_flag,tot_mass_rem)
END DO

rem_frac = stemf
IF (pool_hyfg.EQ.4) rem_frac = rem_frac*pool_grainf
DO idx = 1,nslay
  CALL rem_pool(pool_stem(idx),rem_frac,pool_flag,tot_mass_rem)
END DO

rem_frac = rootstoref
IF (pool_hyfg.EQ.5) rem_frac = rem_frac*pool_grainf
DO idx = 1,nslay
  CALL rem_pool(pool_rootstore(idx),rem_frac,pool_flag,tot_mass_rem)
END DO

rem_frac = rootfiberf
DO idx = 1,nslay
  CALL rem_pool(pool_rootfiber(idx),rem_frac,pool_flag,tot_mass_rem)
END DO

      ! all but fibrous root included in harvest index
IF (pool_flag.EQ.1) THEN
  DO idx = 1,nslay
     sel_mass_left = sel_mass_left + pool_store(idx)
     sel_mass_left = sel_mass_left + pool_leaf(idx)
     sel_mass_left = sel_mass_left + pool_stem(idx)
     sel_mass_left = sel_mass_left + pool_rootstore(idx)
  END DO
END IF

END SUBROUTINE rem_bg_pool
!
SUBROUTINE rem_pool(pool_mass,pool_frac,pool_flag,tot_mass_rem)
!
IMPLICIT NONE
!
! Subroutine arguments
!
INTEGER :: pool_flag
REAL :: pool_frac,pool_mass,tot_mass_rem
!
! Local variables
!
REAL :: mass_rem
!
mass_rem = pool_mass*pool_frac
IF (mass_rem.GT.0.0) THEN
  pool_flag = 1
  pool_mass = pool_mass - mass_rem
  tot_mass_rem = tot_mass_rem + mass_rem
END IF
!
END SUBROUTINE rem_pool
!
SUBROUTINE adj_stand_pool(start_standstem,start_standleaf,start_standstore,     &
                        & start_rootstore,start_rootfiber,pool_standstem,       &
                        & pool_standleaf,pool_standstore,pool_rootstore,        &
                        & pool_rootfiber,pool_flatstem,pool_flatleaf,           &
                        & pool_flatstore,pool_dstm,nslay)
!
IMPLICIT NONE
!
INCLUDE 'p1werm.inc'
!
! Subroutine arguments
!
INTEGER :: nslay
REAL :: pool_dstm,pool_flatleaf,pool_flatstem,pool_flatstore,pool_standleaf,    &
      & pool_standstem,pool_standstore,start_standleaf,start_standstem,         &
      & start_standstore
REAL,DIMENSION(mnsz) :: pool_rootfiber,pool_rootstore,start_rootfiber,          &
                      & start_rootstore
!
! Local variables
!
INTEGER :: idx
REAL :: mov_mass,rat_leaf,rat_root,rat_rootfiber,rat_rootstore,rat_stem,        &
      & rat_store
!
!     + + + purpose + + +
!     this subroutine checks to see if a greater proportion of roots
!    (storage and fiber) have been removed than stems, and if so turns
!     the now unsupoorted stems into flat biomass. the same check is
!     then done for stems supoorting leaves and storage biomass.

!     + + + common blocks + + +

!     + + + argument declarations + + +

!     + + + argument definitions + + +
!     start_standstem - before biomass removal, crop standing stem mass (kg/m^2)
!     start_standleaf - before biomass removal, crop standing leaf mass (kg/m^2)
!     start_standstore - before biomass removal, crop standing storage mass (kg/m^2)
!                    (head with seed, or vegetative head (cabbage, pineapple))
!     start_rootstorez - before biomass removal, crop root storage mass by soil layer (kg/m^2)
!                   (tubers (potatoes, carrots), extended leaf (onion), seeds (peanuts))
!     start_rootfiberz - before biomass removal, crop root fibrous mass by soil layer (kg/m^2)

!     pool_flatstem  - pool flat stem mass (kg/m^2)
!     pool_flatleaf  - pool flat leaf mass (kg/m^2)
!     pool_flatstore - pool flat storage mass (kg/m^2)

!     pool_rootstore - pool flat root storage mass (kg/m^2)
!                   (tubers (potatoes, carrots), extended leaf (onion), seeds (peanuts))
!     pool_rootfiber - pool flat root fibrous mass (kg/m^2)

!     pool_flatstem  - pool flat stem mass (kg/m^2)
!     pool_flatleaf  - pool flat leaf mass (kg/m^2)
!     pool_flatstore - pool flat storage mass (kg/m^2)

!     pool_dstm - number of crop stems per unit area (#/m^2)
!               - it is computed by taking the tillering factor
!                 times the plant population density.

!     nslay - number of soil layers used

!     + + + local variables + + +

!     rat_root - the fraction of material remaining after removal

      ! adjust store, leaf and stem for rootstore or stem removal
IF (start_standstore.GT.0.0) THEN
  rat_store = pool_standstore/start_standstore
ELSE
  rat_store = 1.0
END IF
IF (start_standleaf.GT.0.0) THEN
  rat_leaf = pool_standleaf/start_standleaf
ELSE
  rat_leaf = 1.0
END IF
IF (start_standstem.GT.0.0) THEN
  rat_stem = pool_standstem/start_standstem
ELSE
  rat_stem = 1.0
END IF
rat_rootstore = 1.0
DO idx = 1,nslay
  IF (start_rootstore(idx).GT.0.0)                                              &
    & rat_rootstore = min(rat_rootstore,pool_rootstore(idx)/start_rootstore(idx)&
    & )
END DO
rat_rootfiber = 1.0
DO idx = 1,nslay
  IF (start_rootfiber(idx).GT.0.0)                                              &
    & rat_rootfiber = min(rat_rootfiber,pool_rootfiber(idx)/start_rootfiber(idx)&
    & )
END DO
      ! check if supporting roots removed
rat_root = min(rat_rootstore,rat_rootfiber)
IF (rat_root.LT.rat_stem) THEN
          ! reduce stem count proportionally as well
  pool_dstm = pool_dstm*(rat_root/rat_stem)
          ! move standing mass
  mov_mass = pool_standstem*(1.0-(rat_root/rat_stem))
  pool_flatstem = pool_flatstem + mov_mass
  pool_standstem = pool_standstem - mov_mass
  rat_stem = rat_root
END IF
      ! check if supporting stems removed
IF (rat_stem.LT.rat_leaf) THEN
  mov_mass = pool_standleaf*(1.0-(rat_stem/rat_leaf))
  pool_flatleaf = pool_flatleaf + mov_mass
  pool_standleaf = pool_standleaf - mov_mass
END IF
IF (rat_stem.LT.rat_store) THEN
  mov_mass = pool_standstore*(1.0-(rat_stem/rat_store))
  pool_flatstore = pool_flatstore + mov_mass
  pool_standstore = pool_standstore - mov_mass
END IF
!
END SUBROUTINE adj_stand_pool

max_stars_repo_path: Project Documents/Source Code Original/Remove.f90

max_stars_repo_name: USDA-ARS-WMSRU/upgm-standalone

max_issues_repo_name: USDA-ARS-WMSRU/upgm-standalone

max_forks_repo_name: USDA-ARS-WMSRU/upgm-standalone


text:
module spral_ssmfe_expert_ciface
  use iso_c_binding
  use spral_ssmfe_core_ciface
  use spral_ssmfe_expert
  implicit none

  type, bind(C) :: spral_ssmfe_options
     integer(C_INT) :: array_base
     integer(C_INT) :: print_level
     integer(C_INT) :: unit_error
     integer(C_INT) :: unit_warning
     integer(C_INT) :: unit_diagnostic
     integer(C_INT) :: max_iterations
     integer(C_INT) :: user_x
     integer(C_INT) :: err_est
     real(C_DOUBLE) :: abs_tol_lambda
     real(C_DOUBLE) :: rel_tol_lambda
     real(C_DOUBLE) :: abs_tol_residual
     real(C_DOUBLE) :: rel_tol_residual
     real(C_DOUBLE) :: tol_x
     real(C_DOUBLE) :: left_gap
     real(C_DOUBLE) :: right_gap
     integer(C_INT) :: extra_left
     integer(C_INT) :: extra_right
     integer(C_INT) :: max_left
     integer(C_INT) :: max_right
     logical(C_BOOL) :: minAprod
     logical(C_BOOL) :: minBprod
  end type spral_ssmfe_options

  type ssmfe_expert_ciface_keep
     type(ssmfe_expert_keep) :: keep
     type(ssmfe_rcid) :: rcid
     type(ssmfe_rciz) :: rciz
     type(ssmfe_inform) :: inform
  end type ssmfe_expert_ciface_keep

contains
  subroutine copy_options_in(coptions, foptions, cindexed)
    implicit none
    type(spral_ssmfe_options), intent(in) :: coptions
    type(ssmfe_options), intent(inout) :: foptions ! inherits defaults!
    logical, intent(out) :: cindexed

    cindexed                   = (coptions%array_base .eq. 0)
    foptions%print_level       = coptions%print_level
    foptions%unit_error        = coptions%unit_error
    foptions%unit_warning      = coptions%unit_warning
    foptions%unit_diagnostic   = coptions%unit_diagnostic
    foptions%max_iterations    = coptions%max_iterations
    foptions%user_x            = coptions%user_x
    foptions%err_est           = coptions%err_est
    foptions%abs_tol_lambda    = coptions%abs_tol_lambda
    foptions%rel_tol_lambda    = coptions%rel_tol_lambda
    foptions%abs_tol_residual  = coptions%abs_tol_residual
    foptions%rel_tol_residual  = coptions%rel_tol_residual
    foptions%tol_x             = coptions%tol_x
    foptions%left_gap          = coptions%left_gap
    foptions%right_gap         = coptions%right_gap
    foptions%extra_left        = coptions%extra_left
    foptions%extra_right       = coptions%extra_right
    foptions%max_left          = coptions%max_left
    foptions%max_right         = coptions%max_right
    foptions%minAprod          = coptions%minAprod
    foptions%minBprod          = coptions%minBprod
  end subroutine copy_options_in
end module spral_ssmfe_expert_ciface

subroutine spral_ssmfe_default_options(coptions) bind(C)
  use spral_ssmfe_expert_ciface
  implicit none

  type(spral_ssmfe_options), intent(out) :: coptions

  type(ssmfe_options) :: default_options

  coptions%array_base  = 0 ! C
  coptions%print_level       = default_options%print_level
  coptions%unit_error        = default_options%unit_error
  coptions%unit_warning      = default_options%unit_warning
  coptions%unit_diagnostic   = default_options%unit_diagnostic
  coptions%max_iterations    = default_options%max_iterations
  coptions%user_x            = default_options%user_x
  coptions%err_est           = default_options%err_est
  coptions%abs_tol_lambda    = default_options%abs_tol_lambda
  coptions%rel_tol_lambda    = default_options%rel_tol_lambda
  coptions%abs_tol_residual  = default_options%abs_tol_residual
  coptions%rel_tol_residual  = default_options%rel_tol_residual
  coptions%tol_x             = default_options%tol_x
  coptions%left_gap          = default_options%left_gap
  coptions%right_gap         = default_options%right_gap
  coptions%extra_left        = default_options%extra_left
  coptions%extra_right       = default_options%extra_right
  coptions%max_left          = default_options%max_left
  coptions%max_right         = default_options%max_right
  coptions%minAprod          = default_options%minAprod
  coptions%minBprod          = default_options%minBprod
end subroutine spral_ssmfe_default_options

subroutine spral_ssmfe_expert_standard_double(crci, left, mep, lambda, m, rr, &
     ind, ckeep, coptions, cinform) bind(C)
  use spral_ssmfe_expert_ciface
  implicit none

  type(spral_ssmfe_rcid), intent(inout) :: crci
  integer(C_INT), value :: left
  integer(C_INT), value :: mep
  real(C_DOUBLE), dimension(mep), intent(inout) :: lambda
  integer(C_INT), value :: m
  real(C_DOUBLE), dimension(2*m, 2*m, 3), intent(inout) :: rr
  integer(C_INT), dimension(m), intent(out) :: ind
  type(C_PTR), intent(inout) :: ckeep
  type(spral_ssmfe_options), intent(in) :: coptions
  type(spral_ssmfe_inform), intent(inout) :: cinform

  logical :: cindexed
  type(ssmfe_expert_ciface_keep), pointer :: fcikeep
  type(ssmfe_options) :: foptions

  ! Copy options in first to find out whether we use Fortran or C indexing
  call copy_options_in(coptions, foptions, cindexed)

  ! Translate arguments
  if (c_associated(ckeep)) then
     call c_f_pointer(ckeep, fcikeep)
  else
     allocate(fcikeep)
     ckeep = c_loc(fcikeep)
  end if
  if (crci%job .eq. 0) fcikeep%rcid%job = 0
  if ((fcikeep%rcid%job .eq. 999) .and. (fcikeep%rcid%k .gt. 0)) then
     fcikeep%rcid%i = crci%i
     fcikeep%rcid%j = crci%j
  end if

  ! Call Fortran routine
  call ssmfe_standard(fcikeep%rcid, left, mep, lambda, m, rr, ind, &
       fcikeep%keep, foptions, fcikeep%inform)

  ! Copy arguments out
  call copy_rci_out(fcikeep%rcid, crci, cindexed)
  call copy_inform_out(fcikeep%inform, cinform)
  if ((crci%job .eq.11) .and. cindexed) &
       ind(1:crci%nx) = ind(1:crci%nx) - 1
end subroutine spral_ssmfe_expert_standard_double

subroutine spral_ssmfe_expert_standard_double_complex(crci, left, mep, lambda, &
     m, rr, ind, ckeep, coptions, cinform) bind(C)
  use spral_ssmfe_expert_ciface
  implicit none

  type(spral_ssmfe_rciz), intent(inout) :: crci
  integer(C_INT), value :: left
  integer(C_INT), value :: mep
  real(C_DOUBLE), dimension(mep), intent(inout) :: lambda
  integer(C_INT), value :: m
  complex(C_DOUBLE_COMPLEX), dimension(2*m, 2*m, 3), intent(inout) :: rr
  integer(C_INT), dimension(m), intent(out) :: ind
  type(C_PTR), intent(inout) :: ckeep
  type(spral_ssmfe_options), intent(in) :: coptions
  type(spral_ssmfe_inform), intent(inout) :: cinform

  logical :: cindexed
  type(ssmfe_expert_ciface_keep), pointer :: fcikeep
  type(ssmfe_options) :: foptions

  ! Copy options in first to find out whether we use Fortran or C indexing
  call copy_options_in(coptions, foptions, cindexed)

  ! Translate arguments
  if (c_associated(ckeep)) then
     call c_f_pointer(ckeep, fcikeep)
  else
     allocate(fcikeep)
     ckeep = c_loc(fcikeep)
  end if
  if (crci%job .eq. 0) fcikeep%rciz%job = 0
  if ((fcikeep%rciz%job .eq. 999) .and. (fcikeep%rciz%k .gt. 0)) then
     fcikeep%rciz%i = crci%i
     fcikeep%rciz%j = crci%j
  end if

  ! Call Fortran routine
  call ssmfe_standard(fcikeep%rciz, left, mep, lambda, m, rr, ind, &
       fcikeep%keep, foptions, fcikeep%inform)

  ! Copy arguments out
  call copy_rci_out(fcikeep%rciz, crci, cindexed)
  call copy_inform_out(fcikeep%inform, cinform)
  if ((crci%job .eq. 11) .and. cindexed) &
       ind(1:crci%nx) = ind(1:crci%nx) - 1
end subroutine spral_ssmfe_expert_standard_double_complex

subroutine spral_ssmfe_expert_standard_shift_double(crci, sigma, left, right, &
     mep, lambda, m, rr, ind, ckeep, coptions, cinform) bind(C)
  use spral_ssmfe_expert_ciface
  implicit none

  type(spral_ssmfe_rcid), intent(inout) :: crci
  real(C_DOUBLE), value :: sigma
  integer(C_INT), value :: left
  integer(C_INT), value :: right
  integer(C_INT), value :: mep
  real(C_DOUBLE), dimension(mep), intent(inout) :: lambda
  integer(C_INT), value :: m
  real(C_DOUBLE), dimension(2*m, 2*m, 3), intent(inout) :: rr
  integer(C_INT), dimension(m), intent(out) :: ind
  type(C_PTR), intent(inout) :: ckeep
  type(spral_ssmfe_options), intent(in) :: coptions
  type(spral_ssmfe_inform), intent(inout) :: cinform

  logical :: cindexed
  type(ssmfe_expert_ciface_keep), pointer :: fcikeep
  type(ssmfe_options) :: foptions

  ! Copy options in first to find out whether we use Fortran or C indexing
  call copy_options_in(coptions, foptions, cindexed)

  ! Translate arguments
  if (c_associated(ckeep)) then
     call c_f_pointer(ckeep, fcikeep)
  else
     allocate(fcikeep)
     ckeep = c_loc(fcikeep)
  end if
  if (crci%job .eq. 0) fcikeep%rcid%job = 0
  if ((fcikeep%rcid%job .eq. 999) .and. (fcikeep%rcid%k .gt. 0)) then
     fcikeep%rcid%i = crci%i
     fcikeep%rcid%j = crci%j
  end if

  ! Call Fortran routine
  call ssmfe_standard_shift(fcikeep%rcid, sigma, left, right, mep, lambda, &
       m, rr, ind, fcikeep%keep, foptions, fcikeep%inform)

  ! Copy arguments out
  call copy_rci_out(fcikeep%rcid, crci, cindexed)
  call copy_inform_out(fcikeep%inform, cinform)
  if ((crci%job .eq. 11) .and. cindexed) &
       ind(1:crci%nx) = ind(1:crci%nx) - 1
end subroutine spral_ssmfe_expert_standard_shift_double

subroutine spral_ssmfe_expert_standard_shift_double_complex(crci, sigma, left, &
     right, mep, lambda, m, rr, ind, ckeep, coptions, cinform) bind(C)
  use spral_ssmfe_expert_ciface
  implicit none

  type(spral_ssmfe_rciz), intent(inout) :: crci
  real(C_DOUBLE), value :: sigma
  integer(C_INT), value :: left
  integer(C_INT), value :: right
  integer(C_INT), value :: mep
  real(C_DOUBLE), dimension(mep), intent(inout) :: lambda
  integer(C_INT), value :: m
  complex(C_DOUBLE_COMPLEX), dimension(2*m, 2*m, 3), intent(inout) :: rr
  integer(C_INT), dimension(m), intent(out) :: ind
  type(C_PTR), intent(inout) :: ckeep
  type(spral_ssmfe_options), intent(in) :: coptions
  type(spral_ssmfe_inform), intent(inout) :: cinform

  logical :: cindexed
  type(ssmfe_expert_ciface_keep), pointer :: fcikeep
  type(ssmfe_options) :: foptions

  ! Copy options in first to find out whether we use Fortran or C indexing
  call copy_options_in(coptions, foptions, cindexed)

  ! Translate arguments
  if (c_associated(ckeep)) then
     call c_f_pointer(ckeep, fcikeep)
  else
     allocate(fcikeep)
     ckeep = c_loc(fcikeep)
  end if
  if (crci%job .eq. 0) fcikeep%rciz%job = 0
  if ((fcikeep%rciz%job .eq. 999) .and. (fcikeep%rciz%k .gt. 0)) then
     fcikeep%rciz%i = crci%i
     fcikeep%rciz%j = crci%j
  end if

  ! Call Fortran routine
  call ssmfe_standard_shift(fcikeep%rciz, sigma, left, right, mep, lambda, &
       m, rr, ind, fcikeep%keep, foptions, fcikeep%inform)

  ! Copy arguments out
  call copy_rci_out(fcikeep%rciz, crci, cindexed)
  call copy_inform_out(fcikeep%inform, cinform)
  if ((crci%job .eq. 11) .and. cindexed) &
       ind(1:crci%nx) = ind(1:crci%nx) - 1
end subroutine spral_ssmfe_expert_standard_shift_double_complex

subroutine spral_ssmfe_expert_generalized_double(crci, left, mep, lambda, m, &
     rr, ind, ckeep, coptions, cinform) bind(C)
  use spral_ssmfe_expert_ciface
  implicit none

  type(spral_ssmfe_rcid), intent(inout) :: crci
  integer(C_INT), value :: left
  integer(C_INT), value :: mep
  real(C_DOUBLE), dimension(mep), intent(inout) :: lambda
  integer(C_INT), value :: m
  real(C_DOUBLE), dimension(2*m, 2*m, 3), intent(inout) :: rr
  integer(C_INT), dimension(m), intent(out) :: ind
  type(C_PTR), intent(inout) :: ckeep
  type(spral_ssmfe_options), intent(in) :: coptions
  type(spral_ssmfe_inform), intent(inout) :: cinform

  logical :: cindexed
  type(ssmfe_expert_ciface_keep), pointer :: fcikeep
  type(ssmfe_options) :: foptions

  ! Copy options in first to find out whether we use Fortran or C indexing
  call copy_options_in(coptions, foptions, cindexed)

  ! Translate arguments
  if (c_associated(ckeep)) then
     call c_f_pointer(ckeep, fcikeep)
  else
     allocate(fcikeep)
     ckeep = c_loc(fcikeep)
  end if
  if (crci%job .eq. 0) fcikeep%rcid%job = 0
  if ((fcikeep%rcid%job .eq. 999) .and. (fcikeep%rcid%k .gt. 0)) then
     fcikeep%rcid%i = crci%i
     fcikeep%rcid%j = crci%j
  end if

  ! Call Fortran routine
  call ssmfe_generalized(fcikeep%rcid, left, mep, lambda, m, rr, ind, &
       fcikeep%keep, foptions, fcikeep%inform)

  ! Copy arguments out
  call copy_rci_out(fcikeep%rcid, crci, cindexed)
  call copy_inform_out(fcikeep%inform, cinform)
  if ((crci%job .eq. 11) .and. cindexed) &
       ind(1:crci%nx) = ind(1:crci%nx) - 1
end subroutine spral_ssmfe_expert_generalized_double

subroutine spral_ssmfe_expert_generalized_double_complex(crci, left, mep, &
     lambda, m, rr, ind, ckeep, coptions, cinform) bind(C)
  use spral_ssmfe_expert_ciface
  implicit none

  type(spral_ssmfe_rciz), intent(inout) :: crci
  integer(C_INT), value :: left
  integer(C_INT), value :: mep
  real(C_DOUBLE), dimension(mep), intent(inout) :: lambda
  integer(C_INT), value :: m
  complex(C_DOUBLE_COMPLEX), dimension(2*m, 2*m, 3), intent(inout) :: rr
  integer(C_INT), dimension(m), intent(out) :: ind
  type(C_PTR), intent(inout) :: ckeep
  type(spral_ssmfe_options), intent(in) :: coptions
  type(spral_ssmfe_inform), intent(inout) :: cinform

  logical :: cindexed
  type(ssmfe_expert_ciface_keep), pointer :: fcikeep
  type(ssmfe_options) :: foptions

  ! Copy options in first to find out whether we use Fortran or C indexing
  call copy_options_in(coptions, foptions, cindexed)

  ! Translate arguments
  if (c_associated(ckeep)) then
     call c_f_pointer(ckeep, fcikeep)
  else
     allocate(fcikeep)
     ckeep = c_loc(fcikeep)
  end if
  if (crci%job .eq. 0) fcikeep%rciz%job = 0
  if ((fcikeep%rciz%job .eq. 999) .and. (fcikeep%rciz%k .gt. 0)) then
     fcikeep%rciz%i = crci%i
     fcikeep%rciz%j = crci%j
  end if

  ! Call Fortran routine
  call ssmfe_generalized(fcikeep%rciz, left, mep, lambda, m, rr, ind, &
       fcikeep%keep, foptions, fcikeep%inform)

  ! Copy arguments out
  call copy_rci_out(fcikeep%rciz, crci, cindexed)
  call copy_inform_out(fcikeep%inform, cinform)
  if ((crci%job .eq. 11) .and. cindexed) &
       ind(1:crci%nx) = ind(1:crci%nx) - 1
end subroutine spral_ssmfe_expert_generalized_double_complex

subroutine spral_ssmfe_expert_generalized_shift_double(crci, sigma, left, &
     right, mep, lambda, m, rr, ind, ckeep, coptions, cinform) bind(C)
  use spral_ssmfe_expert_ciface
  implicit none

  type(spral_ssmfe_rcid), intent(inout) :: crci
  real(C_DOUBLE), value :: sigma
  integer(C_INT), value :: left
  integer(C_INT), value :: right
  integer(C_INT), value :: mep
  real(C_DOUBLE), dimension(mep), intent(inout) :: lambda
  integer(C_INT), value :: m
  real(C_DOUBLE), dimension(2*m, 2*m, 3), intent(inout) :: rr
  integer(C_INT), dimension(m), intent(out) :: ind
  type(C_PTR), intent(inout) :: ckeep
  type(spral_ssmfe_options), intent(in) :: coptions
  type(spral_ssmfe_inform), intent(inout) :: cinform

  logical :: cindexed
  type(ssmfe_expert_ciface_keep), pointer :: fcikeep
  type(ssmfe_options) :: foptions

  ! Copy options in first to find out whether we use Fortran or C indexing
  call copy_options_in(coptions, foptions, cindexed)

  ! Translate arguments
  if (c_associated(ckeep)) then
     call c_f_pointer(ckeep, fcikeep)
  else
     allocate(fcikeep)
     ckeep = c_loc(fcikeep)
  end if
  if (crci%job .eq. 0) fcikeep%rcid%job = 0
  if ((fcikeep%rcid%job .eq. 999) .and. (fcikeep%rcid%k .gt. 0)) then
     fcikeep%rcid%i = crci%i
     fcikeep%rcid%j = crci%j
  end if

  ! Call Fortran routine
  call ssmfe_generalized_shift(fcikeep%rcid, sigma, left, right, mep, lambda, &
       m, rr, ind, fcikeep%keep, foptions, fcikeep%inform)

  ! Copy arguments out
  call copy_rci_out(fcikeep%rcid, crci, cindexed)
  call copy_inform_out(fcikeep%inform, cinform)
  if ((crci%job .eq. 11) .and. cindexed) &
       ind(1:crci%nx) = ind(1:crci%nx) - 1
end subroutine spral_ssmfe_expert_generalized_shift_double

subroutine spral_ssmfe_expert_generalized_shift_double_complex(crci, sigma, &
     left, right, mep, lambda, m, rr, ind, ckeep, coptions, cinform) bind(C)
  use spral_ssmfe_expert_ciface
  implicit none

  type(spral_ssmfe_rciz), intent(inout) :: crci
  real(C_DOUBLE), value :: sigma
  integer(C_INT), value :: left
  integer(C_INT), value :: right
  integer(C_INT), value :: mep
  real(C_DOUBLE), dimension(mep), intent(inout) :: lambda
  integer(C_INT), value :: m
  complex(C_DOUBLE_COMPLEX), dimension(2*m, 2*m, 3), intent(inout) :: rr
  integer(C_INT), dimension(m), intent(out) :: ind
  type(C_PTR), intent(inout) :: ckeep
  type(spral_ssmfe_options), intent(in) :: coptions
  type(spral_ssmfe_inform), intent(inout) :: cinform

  logical :: cindexed
  type(ssmfe_expert_ciface_keep), pointer :: fcikeep
  type(ssmfe_options) :: foptions

  ! Copy options in first to find out whether we use Fortran or C indexing
  call copy_options_in(coptions, foptions, cindexed)

  ! Translate arguments
  if (c_associated(ckeep)) then
     call c_f_pointer(ckeep, fcikeep)
  else
     allocate(fcikeep)
     ckeep = c_loc(fcikeep)
  endif
  if (crci%job .eq. 0) fcikeep%rciz%job = 0
  if ((fcikeep%rciz%job .eq. 999) .and. (fcikeep%rciz%k .gt. 0)) then
     fcikeep%rciz%i = crci%i
     fcikeep%rciz%j = crci%j
  end if

  ! Call Fortran routine
  call ssmfe_generalized_shift(fcikeep%rciz, sigma, left, right, mep, lambda, &
       m, rr, ind, fcikeep%keep, foptions, fcikeep%inform)

  ! Copy arguments out
  call copy_rci_out(fcikeep%rciz, crci, cindexed)
  call copy_inform_out(fcikeep%inform, cinform)
  if ((crci%job .eq. 11) .and. cindexed) &
       ind(1:crci%nx) = ind(1:crci%nx) - 1
end subroutine spral_ssmfe_expert_generalized_shift_double_complex

subroutine spral_ssmfe_expert_buckling_double(crci, sigma, left, right, mep, &
     lambda, m, rr, ind, ckeep, coptions, cinform) bind(C)
  use spral_ssmfe_expert_ciface
  implicit none

  type(spral_ssmfe_rcid), intent(inout) :: crci
  real(C_DOUBLE), value :: sigma
  integer(C_INT), value :: left
  integer(C_INT), value :: right
  integer(C_INT), value :: mep
  real(C_DOUBLE), dimension(mep), intent(inout) :: lambda
  integer(C_INT), value :: m
  real(C_DOUBLE), dimension(2*m, 2*m, 3), intent(inout) :: rr
  integer(C_INT), dimension(m), intent(out) :: ind
  type(C_PTR), intent(inout) :: ckeep
  type(spral_ssmfe_options), intent(in) :: coptions
  type(spral_ssmfe_inform), intent(inout) :: cinform

  logical :: cindexed
  type(ssmfe_expert_ciface_keep), pointer :: fcikeep
  type(ssmfe_options) :: foptions

  ! Copy options in first to find out whether we use Fortran or C indexing
  call copy_options_in(coptions, foptions, cindexed)

  ! Translate arguments
  if (c_associated(ckeep)) then
     call c_f_pointer(ckeep, fcikeep)
  else
     allocate(fcikeep)
     ckeep = c_loc(fcikeep)
  end if
  if (crci%job .eq. 0) fcikeep%rcid%job = 0
  if ((fcikeep%rcid%job .eq. 999) .and. (fcikeep%rcid%k .gt. 0)) then
     fcikeep%rcid%i = crci%i
     fcikeep%rcid%j = crci%j
  end if

  ! Call Fortran routine
  call ssmfe_buckling(fcikeep%rcid, sigma, left, right, mep, lambda, &
       m, rr, ind, fcikeep%keep, foptions, fcikeep%inform)

  ! Copy arguments out
  call copy_rci_out(fcikeep%rcid, crci, cindexed)
  call copy_inform_out(fcikeep%inform, cinform)
  if ((crci%job .eq. 11) .and. cindexed) &
       ind(1:crci%nx) = ind(1:crci%nx) - 1
end subroutine spral_ssmfe_expert_buckling_double

subroutine spral_ssmfe_expert_buckling_double_complex(crci, sigma, left, &
     right, mep, lambda, m, rr, ind, ckeep, coptions, cinform) bind(C)
  use spral_ssmfe_expert_ciface
  implicit none

  type(spral_ssmfe_rciz), intent(inout) :: crci
  real(C_DOUBLE), value :: sigma
  integer(C_INT), value :: left
  integer(C_INT), value :: right
  integer(C_INT), value :: mep
  real(C_DOUBLE), dimension(mep), intent(inout) :: lambda
  integer(C_INT), value :: m
  complex(C_DOUBLE_COMPLEX), dimension(2*m, 2*m, 3), intent(inout) :: rr
  integer(C_INT), dimension(m), intent(out) :: ind
  type(C_PTR), intent(inout) :: ckeep
  type(spral_ssmfe_options), intent(in) :: coptions
  type(spral_ssmfe_inform), intent(inout) :: cinform

  logical :: cindexed
  type(ssmfe_expert_ciface_keep), pointer :: fcikeep
  type(ssmfe_options) :: foptions

  ! Copy options in first to find out whether we use Fortran or C indexing
  call copy_options_in(coptions, foptions, cindexed)

  ! Translate arguments
  if (c_associated(ckeep)) then
     call c_f_pointer(ckeep, fcikeep)
  else
     allocate(fcikeep)
     ckeep = c_loc(fcikeep)
  end if
  if (crci%job .eq. 0) fcikeep%rciz%job = 0
  if ((fcikeep%rciz%job .eq. 999) .and. (fcikeep%rciz%k .gt. 0)) then
     fcikeep%rciz%i = crci%i
     fcikeep%rciz%j = crci%j
  end if

  ! Call Fortran routine
  call ssmfe_buckling(fcikeep%rciz, sigma, left, right, mep, lambda, &
       m, rr, ind, fcikeep%keep, foptions, fcikeep%inform)

  ! Copy arguments out
  call copy_rci_out(fcikeep%rciz, crci, cindexed)
  call copy_inform_out(fcikeep%inform, cinform)
  if ((crci%job .eq. 11) .and. cindexed) &
       ind(1:crci%nx) = ind(1:crci%nx) - 1
end subroutine spral_ssmfe_expert_buckling_double_complex

subroutine spral_ssmfe_expert_free(ckeep, cinform) bind(C)
  use spral_ssmfe_expert_ciface
  implicit none

  type(C_PTR), intent(inout) :: ckeep
  type(spral_ssmfe_inform), intent(inout) :: cinform

  type(ssmfe_expert_ciface_keep), pointer :: fcikeep

  ! Nullify pointer components of cinform
  cinform%converged       = C_NULL_PTR
  cinform%residual_norms  = C_NULL_PTR
  cinform%err_lambda      = C_NULL_PTR
  cinform%err_x           = C_NULL_PTR

  ! Check ckeep is not null
  if(.not. c_associated(ckeep)) return

  ! Associate fkeep
  call c_f_pointer(ckeep, fcikeep)

  ! Call Fortran cleanup
  call ssmfe_free(fcikeep%keep, fcikeep%inform)

  ! Free fcikeep
  deallocate(fcikeep)
  ckeep = C_NULL_PTR
end subroutine spral_ssmfe_expert_free

max_stars_repo_name: mjacobse/spral

max_issues_repo_name: mjacobse/spral

max_forks_repo_name: mjacobse/spral


program test_string

use, intrinsic:: iso_c_binding, only: c_null_char
use string_utils

implicit none (type, external)

call test_lowercase()
print *,'PASSED: HDF5 character'
call test_strip_null()
print *,'PASSED: null strip'

contains

subroutine test_lowercase()

character(*), parameter :: hello = 'HeLl0 Th3rE !>? '
  !! Fortran 2003 allocatable string

if (.not.(toLower(hello)=='hell0 th3re !>? ')) error stop 'error: lowercase conversion'

if (.not.(trim(toLower(hello))=='hell0 th3re !>?')) error stop 'Allocatable lowercase conversion error'

end subroutine test_lowercase

subroutine test_strip_null()
character(*), parameter :: hello = 'HeLl0 Th3rE !>? '

if (.not.strip_trailing_null(hello // c_null_char) == hello) error stop 'problem stripping trailing null'

end subroutine test_strip_null

end program


text:
program dv_alpha_relax

  use confort
  use globals
  use settings
  use iso_fortran_env, only: sp => real32, dp => real64
  use ieee_arithmetic, only: ieee_is_normal, ieee_is_nan
  use fileunits
  use relaxation
  use ss73solution, only: apxdisk, apx_estim, apx_refin
  use grid, only: space_linear, space_linlog

  !----------------------------------------------------------------------------!
  implicit none

  type(config) :: cfg
  integer :: model, errno
  integer :: ny = 3, i, globiter = 0
  integer, dimension(2) :: niter = [ 48, 12 ]
  real(dp), allocatable, target :: x(:), x0(:), Y(:), dY(:), M(:,:)
  real(dp), pointer, dimension(:) :: y_rho, y_temp, y_frad, tau, y_trad
  real(dp) :: rhoc, Tc, Hdisk
  logical :: user_ff, user_bf
  integer, dimension(6) :: C_

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

  ngrid = 576
  htop = 5

  !----------------------------------------------------------------------------!
  ! read command line

  parse_command_line: block

    integer :: i
    character(2**8) :: arg

    call rdargvgl

    do i = 1, command_argument_count()
      call get_command_argument(i, arg)
      select case (arg)

      ! use P_rad in alpha prescription?
      case ("-prad-alpha", "-alpha-prad")
        use_prad_in_alpha = .TRUE.
      case ("-no-prad-alpha", "-no-alpha-prad")
        use_prad_in_alpha = .FALSE.

      end select
    end do

  end block parse_command_line

  !----------------------------------------------------------------------------!
  ! read config

  call mincf_read(cfg)
  call rdconfgl(cfg)
  call rdconf(cfg)
  call mincf_free(cfg)

  user_ff = use_opacity_ff
  user_bf = use_opacity_bf

  !----------------------------------------------------------------------------!
  ! calculate the global parameters
  call cylinder(mbh, mdot, radius, rschw, omega, facc, teff, zscale)

  ! get the model number
  model = mrx_number( 'D', .FALSE., .FALSE. )
  call mrx_sel_nvar(model, ny)
  call mrx_sel_hash(model, C_)

  allocate( x(ngrid), x0(ngrid), Y(ny*ngrid), dY(ny*ngrid), M(ny*ngrid,ny*ngrid), tau(ngrid) )

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

  call apx_estim(mbh, mdot, radius, alpha, rhoc, Tc, Hdisk)
  call apx_refin(mbh, mdot, radius, alpha, rhoc, Tc, Hdisk)

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

  forall (i = 1:ngrid)
    x(i) = space_linlog(i,ngrid,htop) * Hdisk
    x0(i) = space_linlog(i,ngrid,htop) / htop
  end forall

  y_rho   => Y(C_(1)::ny)
  y_temp  => Y(C_(2)::ny)
  y_trad  => Y(C_(3)::ny)
  y_frad  => Y(C_(4)::ny)

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

  forall (i = 1:ngrid)
    y_frad(i) = x0(i) * facc
    y_temp(i) = (1 - x0(i)) * (Tc - 0.841 * Teff) + 0.841 * Teff
    y_rho(i) =  rhoc * (exp(-0.5*(x(i)/Hdisk)**2) + 1e-8)
  end forall

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

  use_opacity_ff = .false.
  use_opacity_bf = .false.

  call mrx_relax(model, x, y, niter(1), 0d0, errno)

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

  use_opacity_ff = user_ff
  use_opacity_bf = user_bf

  if ( use_opacity_ff .or. use_opacity_bf ) then

    write (uerr,'(" --- ",A)') 'bf+ff opacity is ON'
    call mrx_relax(model, x, y, niter(2), 0.5d0, errno)

  end if

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

  write (uerr,'(" --- ",A)') 'complete'

  open(33, file = trim(outfn) // '.dat', action = 'write')
  call wrdat(33)
  close(33)

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

  open(33, file = trim(outfn) // '.txt', &
      action = 'write')
  call wrpar(33)
  call wpar_gl(33)
  close(33)

  open(34, file = trim(outfn) // ".col", action = 'write')
  call wrcol(34)
  close(34)

  deallocate(x,x0,Y,M,dY,tau)

contains

  subroutine wrdat(u)
    integer, intent(in) :: u
    integer :: i
    real(dp) :: pgas,prad
    call mktaues(x,y_rho,y_temp,tau)
    do i = 1,ngrid
      pgas =  cgs_k_over_mh / miu * y_rho(i) * y_temp(i)
      prad = cgs_a / 3 * y_trad(i)**4
      write (u,'(I7,*(ES14.5E3))') i, x(i), x(i) / zscale, tau(i), &
      y_rho(i), y_temp(i), y_frad(i),   &
      fksct(y_rho(i), y_temp(i)), fkabs(y_rho(i), y_temp(i)), &
      fkcnd(y_rho(i), y_temp(i)), pgas, prad
    end do
  end subroutine

  subroutine wrpar(u)
    integer, intent(in) :: u
    write (u, fmpari) "model", model
    write (u, fmpari) "niter", globiter
    write (u, fmpare) "radius", radius
    write (u, fmpare) "zscale", fzscale(mbh,mdot,radius)
    write (u, fmparec) "rho_0", rhoc, "Central density"
    write (u, fmparec) "temp_0", Tc, "Central gas temperature"
    write (u, fmparec) "Hdisk", Hdisk, "Disk height [cm]"
    write (u, fmparfc) "alpha", alpha, "Alpha parameter"
    write (u, fmparl) "has_corona", .FALSE.
    write (u, fmparl) "has_magnetic", .FALSE.
    write (u, fmparl) "has_conduction", .FALSE.
  end subroutine

  subroutine wrcol(u)
    integer, intent(in) :: u
    write(u, fmcol) 'i', 'i4'
    write(u, fmcol) 'z', 'f4'
    write(u, fmcol) 'h', 'f4'
    write(u, fmcol) 'tau', 'f4'
    write(u, fmcol) 'rho', 'f4'
    write(u, fmcol) 'temp', 'f4'
    write(u, fmcol) 'frad', 'f4'
    write(u, fmcol) 'ksct', 'f4'
    write(u, fmcol) 'kabs', 'f4'
    write(u, fmcol) 'kcnd', 'f4'
    write(u, fmcol) 'pgas', 'f4'
    write(u, fmcol) 'prad', 'f4'
  end subroutine

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

  subroutine rdconf(cfg)
    integer :: errno
    type(config), intent(inout) :: cfg
    character(len=2048) :: buf

    call mincf_get(cfg, "alpha", buf, errno)
    if ( iand(errno, mincf_not_found) .ne. 0 )  then
      error stop "Magnetic alpha-parameter (key: alpha) is REQUIRED!"
    end if
    read (buf,*) alpha

    call mincf_get(cfg, "radius", buf, errno)
    if ( iand(errno, mincf_not_found) .ne. 0 )  then
      error stop "Radius (key: radius) is REQUIRED!"
    end if
    read (buf,*) radius

  end subroutine

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

  subroutine mktaues(z, rho, T, tau)
    real(dp), intent(in), dimension(:) :: z, rho,T
    real(dp), intent(out), dimension(:) :: tau
    real(dp) :: dz, rhom, Tm
    integer :: i
    tau(size(tau)) = 1e-12
    do i = size(tau)-1, 1, -1
      rhom = (rho(i+1) + rho(i)) / 2
      Tm = (T(i+1) + T(i)) / 2
      dz = z(i) - z(i+1)
      tau(i) = tau(i+1) - kapes0(abuX,abuZ) * rhom * dz
    end do
  end subroutine

end program dv_alpha_relax

max_stars_repo_name: gronki/pydiskvert

max_issues_repo_name: gronki/pydiskvert

max_forks_repo_name: gronki/pydiskvert


! Mikayla Webber
! 4670 Numerical Analysis
! Homework IVP Due 12.6.17

program testeuler
  implicit none
  integer :: i, neq, nsteps
  double precision, allocatable, dimension(:,:) :: archive
  double precision :: tin, tout, dt, a, b, x, y, xactual, yactual
  double precision, allocatable, dimension(:) :: yin, yout, t, f1, f2, u

  neq = 1
  nsteps = 50
  allocate(yin(neq), yout(neq), t(0:nsteps), u(neq), f1(neq), f2(neq), archive(neq, 0:nsteps))

  a = 1.0d0
  b = 2.0d0
  archive(:,0) = (/ 1.0d0, 0.0d0, 0.0d0, 1.0d0 /)
  dt = (b - a) / dble(nsteps)

  do i = 0, nsteps
    t(i) = (a + dble(i) * dt)
  end do

  archive(:,0) = (/ 1.0d0, 0.0d0, 0.0d0, 1.0d0 /)
  yin = archive(:,0)

  do i = 0, (nsteps - 1)
    tin = t(i)
    tout = t(i + 1)
    yin = archive(:,i)
    call computef(neq, tin, yin, f1)
    yout = (yin + (tout - tin) * (f1 + f2) / 2.0d0)
    archive(:,i+1) = yout
  end do

  print*, 'creating "out1"'
  open(unit=8, file='q1out1', status='replace')

  do i = 0, nsteps
    write(8,*) t(i), archive(1,i)
  end do

  close(8)

  print*, 'writing pcerror'
  open(unit=8, file='pcerror', status='replace')

  do i = 0, nsteps
    x = archive(1,i)
    y = archive(2,i)
    xactual = cos(t(i))
    yactual = t(i) / (1 + log(t(i)))
    write(8,*) i, sqrt((x - xactual)**2 + (y - yactual)**2)
  end do

  deallocate(yin, yout, archive, t, f1, f2)
  stop
end program testeuler

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

subroutine computef(neq, t, y, f)
  implicit none
  integer :: neq
  double precision :: f(neq)
  double precision, intent(in) :: t, y(neq)

  f(1) = (((1 + y(1)) / t) - (y(1) / t)**2)

  return
end


text:
!> 计算人工粘度的子程序。详见 Monaghan (1992), Hernquist 和 Katz (1989) 或第 4 章中的论述。
!>     subroutine to calculate the artificial viscosity (monaghan, 1992)
subroutine art_visc(ntotal, hsml, mass, x, vx, niac, rho, c, pair_i, pair_j, w, dwdx, dvxdt, dedt)

    use sph_kind, only: rk
    use parameter
    implicit none

    !> 在模拟中所使用的粒子总数
    !> number of particles in simulation
    integer, intent(in) :: ntotal
    !> 光滑长度
    !> smoothing length
    real(rk), intent(in) :: hsml(maxn)
    !> 粒子的质量
    !> particle masses
    real(rk), intent(in) :: mass(maxn)
    !> 粒子的坐标
    !> particle coordinates
    real(rk), intent(in) :: x(dim, maxn)
    !> 粒子的速度
    !> particle velocities
    real(rk), intent(in) :: vx(dim, maxn)
    !> 相互作用对的数目
    !> number of interaction pairs
    integer, intent(in) :: niac
    !> 密度
    !> density
    real(rk), intent(in) :: rho(maxn)
    !> 温度
    !> temperature
    real(rk), intent(in) :: c(maxn)
    !> 相互作用对的第一个粒子
    !> first partner of interaction pair
    integer, intent(in) :: pair_i(max_interaction)
    !> 相互作用对的第二个粒子
    !> second partner of interaction pair
    integer, intent(in) :: pair_j(max_interaction)
    !> 相互作用对的核函数
    !> kernel for all interaction pairs
    real(rk), intent(in) :: w(max_interaction)
    !> 相互作用对的核函数的导数
    !> derivative of kernel with respect to x, y and z
    real(rk), intent(in) :: dwdx(dim, max_interaction)
    !> 相互作用对的加速度
    !> acceleration with respect to x, y and z
    real(rk), intent(out) :: dvxdt(dim, maxn)
    !> 改变特定内部能量
    !> change of specific internal energy
    real(rk), intent(out) :: dedt(maxn)

    integer :: i, j, k, d
    real(rk) :: dx, dvx(dim), alpha, beta, etq, piv, muv, vr, rr, h, mc, mrho, mhsml

    !     parameter for the artificial viscosity:
    !     shear viscosity
    parameter(alpha=1._rk)

    !     bulk viscosity
    parameter(beta=1._rk)

    !     parameter to avoid singularities
    parameter(etq=0.1_rk)

    do i = 1, ntotal
        do d = 1, dim
            dvxdt(d, i) = 0._rk
        end do
        dedt(i) = 0._rk
    end do

    !     calculate sph sum for artificial viscosity

    do k = 1, niac
        i = pair_i(k)
        j = pair_j(k)
        mhsml = (hsml(i) + hsml(j))/2._rk
        vr = 0._rk
        rr = 0._rk
        do d = 1, dim
            dvx(d) = vx(d, i) - vx(d, j)
            dx = x(d, i) - x(d, j)
            vr = vr + dvx(d)*dx
            rr = rr + dx*dx
        end do

        !     artificial viscous force only if v_ij * r_ij < 0

        if (vr < 0._rk) then

            !     calculate muv_ij = hsml v_ij * r_ij / ( r_ij^2 + hsml^2 etq^2 )

            muv = mhsml*vr/(rr + mhsml*mhsml*etq*etq)

            !     calculate piv_ij = (-alpha muv_ij c_ij + beta muv_ij^2) / rho_ij

            mc = 0.5_rk*(c(i) + c(j))
            mrho = 0.5_rk*(rho(i) + rho(j))
            piv = (beta*muv - alpha*mc)*muv/mrho

            !     calculate sph sum for artificial viscous force

            do d = 1, dim
                h = -piv*dwdx(d, k)
                dvxdt(d, i) = dvxdt(d, i) + mass(j)*h
                dvxdt(d, j) = dvxdt(d, j) - mass(i)*h
                dedt(i) = dedt(i) - mass(j)*dvx(d)*h
                dedt(j) = dedt(j) - mass(i)*dvx(d)*h
            end do
        end if
    end do

    !     change of specific internal energy:

    do i = 1, ntotal
        dedt(i) = 0.5_rk*dedt(i)
    end do

end subroutine art_visc

max_stars_repo_name: zoziha/SPH

max_issues_repo_name: zoziha/SPH

max_forks_repo_name: zoziha/SPH


text:
      program main
c
c ... array declarations.
c
      implicit double precision (a-h, o-z)
      parameter (nxx = 20)
      parameter (nlen = nxx*nxx)
      parameter (nwk = 10*nlen)
      parameter (inwk = 3*nlen)
      dimension coef(nlen,3), rhs(nlen), u(nlen), wksp(nwk), ubar(1),
     a     rparm(16)
      integer jcoef(3), iwksp(inwk), iparm(25), p(1), ip(1)
      external cg, bic2
c
      ndim = nlen
      mdim = 3
      nw = nwk
      inw = inwk
      zeta = 1.0d-5
      nout = 6
c
c ... generate coef, jcoef, and rhs.
c
      nx = nxx
      ny = nx
      n = nx*ny
      h = 1.0d0/dble(nx + 1)
      maxnz = 3
      do i = 1,n
         coef(i,1) = 6.0d0
         coef(i,2) = -1.0d0
         coef(i,3) = -2.0d0
         rhs(i) = 0.0d0
      enddo
      k = 0
      do j = 1,ny
         y = dble(j)*h
         do i = 1,nx
            x = dble(i)*h
            k = k + 1
            if (j .eq. 1) then
               rhs(k) = rhs(k) + 2.0d0
            endif
            if (j .eq. ny) then
               rhs(k) = rhs(k) + 2.0d0*(1.0d0 + x)
               coef(k,3) = 0.0d0
            endif
            if (i .eq. 1) then
               rhs(k) = rhs(k) + 1.0d0
            endif
            if (i .eq. nx) then
               rhs(k) = rhs(k) + 1.0d0 + y
               coef(k,2) = 0.0d0
            endif
         enddo
      enddo
      jcoef(1) = 0
      jcoef(2) = 1
      jcoef(3) = nx
      call dfault (iparm,rparm)
c
c ... now, reset some default values.
c
      iparm(3) = 3
      iparm(4) = nout
      iparm(18) = 1
      iparm(19) = nx
      rparm(1) = zeta
c
c ... generate an initial guess for u and call nspcg.
c
      call vfill (n,u,0.0d0)
c
      call nspcg (bic2,cg,ndim,mdim,n,maxnz,coef,jcoef,p,ip,
     a            u,ubar,rhs,wksp,iwksp,nw,inw,iparm,rparm,ier)
      write (nout,35) nw, inw
 35   format (/1x,'nw  = ',i8/
     a         1x,'inw = ',i8)
      stop
      end

ext: f