!Fortran code for P.B. Parks diamond shell ablation model

!gives ablation rate in atoms/sec

! Tpell is plasmas temperature in eV  

! npell is plasma density in m^-3

! rsphere is pellet radius

! mz is mass of Carbon in amu (=12.0) 

          REAL(r8), PARAMETER ::C0=8.146777e-9_r8, ZstarPlus1C = 2.86_r8, &

     &    IstC = 60.0_r8, xiexp = 0.601_r8, lamdaa = 0.0933979540623963,     &

     &    lamdab = -0.7127242270013098, lamdac = -0.2437544205933372,   &

     &    lamdad = -0.8534855445478313, fugCG = 0.777686  

          REAL(r8), PARAMETER :: onethird=1.0_r8/3.0_r8, AvN=6.0221409e+23

          REAL(r8), PARAMETER :: twothird=2.0_r8/3.0_r8, fourthird=4.0_r8/3.0_r8

          Albedo = 23.920538030089528 * log(1.0_r8 + 0.20137080524063228_r8 * ZstarPlus1C)

          flelectro = exp(-1.936)

          fL = (1.0_r8 - Albedo / 100.0_r8) * flelectro

          Tcut = MAX(Tpell,30.0_r8)

          loglamCSlow = log(2.0 * Tcut / IstC * sqrt(2.71828 * 2.0))

          BLamdaq = 1.0_r8 / (REAL(zimp2) * loglamCSlow) * (4.0_r8 / (2.5_r8 + 2.2_r8 * sqrt(ZstarPlus1C)))

          Gpr = (C0* ((mz)**twothird) * (gamm1**onethird) *((fL*npell/1e6_r8)**onethird) *   &

     &         (rsphere**fourthird) * (Tpell**(5.5_r8*onethird)) * (BLamdaq**twothird) )

          av = 10.420403555938629 * ((Tcut / 2000.0)**lamdaa)

          bv = 0.6879779829877795 * ((Tcut / 2000.0)**lamdab)

          cv = 1.5870910225610804 * ((Tcut / 2000.0)**lamdac)

          dv = 2.9695640286641840 * ((Tcut / 2000.0)**lamdad)

          CG = ( fugCG * av * log(1.0_r8 + bv * ((npell/1e20_r8)**twothird)  *        &

     &        (rsphere**twothird) ) / log(cv + dv * ((npell/1e20_r8)**twothird)  *        &

     &         (rsphere**twothird) ) )

          rabl = (AvN/(mz)) * xiexp * CG * Gpr    !ablation rate in atoms/sec