c==============================================================================
C	Model of firn gravitational and thermal fractionation
C       driven by seasonal temperature fluctuations or glacial terminations
C------------------------------------------------------------------------------
C     Explanation:
C       This is finite difference scheme for the time-dependent 
c       diffusion equation, based on Alley and Koci (1990),
c
C     References:
c	Alley, R.B., and Koci, B.R., Recent warming in central Greenland? Ann. 
c	  Glaciol. 14: (1990).
c       Bender, M.L., Sowers, T., Barnola, J.-M., and Chappellaz, J., 1994.
c         Changes in the O2/N2 ratio of the atmosphere during recent decades 
c         reflected in the composition of air in the firn at Vostok Station,
c         Antarctica, Geophysical Research Letters 21:189-192.
c       Grew, K.E., and Ibbs, T.L., 1952.  Thermal diffusion in gases, 
c         Cambridge University Press.
c       Reid, R.C., Prausnitz, J.M., and Sherwood, T.K., 1977. The properties
c         of Gases and Liquids, Third Edition, New York: McGraw-Hill.
c       Schwander, J., Stauffer, B., and Sigg, A., 1988. Air mixing in firn and
c         the age of the air at pore close-off, Annals of Glaciology 10:141-145
c       Schwander, J., Barnola, J.-M., Andrie, C., Leuenberger, M., Ludin, A.,
c         Raynaud, D., and Stauffer, B., 1993. The Age of the Air in the Firn 
c         and the Ice at Summit, Greenland, J. Geophys. Res. 98:2831-2838.
c       Sowers, T., Bender, M., and Raynaud, D., 1989. Elemental and isotopic
c         composition of occluded O2 and N2 in polar ice, JGR 94:5137-5150.
c       Yen, Y.-C., 1981. Review of thermal properties of snow, ice, and sea
c         ice, CRREL Report 81-10, US Army Corps of Engineers, Cold Regions 
c         Research and Engineering Laboratory, Hanover, N.H.
c       G. WELLER and P. SCHWERDTFEGER, 1970, Thermal properties and heat 
c         transfer processes of the snow of the central Antarctic plateau,
c         International Association of Hydrological Sciences Publication
c       Cuffey K.M., et al, 1994, Calibration of the $\delta^18$O isotopic 
c         paleothermometer for central Greenland, using borehole temperatures,
c         Journal of Glaciology
c       Severinghaus J.P., et al, 1998, Timing of abrupt climate change at the
c         end of the Younger Dryas interval from thermally fractionated gases
c         in polar ice , Nature
c       Orsi, A. J., B. D. Cornuelle, and J. P. Severinghaus (2012), Little Ice
c         Age cold interval in West Antarctica: Evidence from borehole temperature
c         at the West Antarctic Ice Sheet (WAIS) Divide, Geophys. Res. Lett., 39, 
c         L09710, doi:10.1029/2012GL051260. 
c
c     Accompanying files:
c     'inputdata_MegaWAIS_11wsc.txt'  input file for parameters
c     'inputfiles/indfreq5.txt' parameters for the fourier basis functions
c     'ofile' is the input temperature history
c     'fileiniT' is the input temperature profile

c
c     Units: length          meters
c            concentration   permil with respect to atmosphere
c            time            s  (except sometimes I use years, as noted)
c            mass            kg (except mass difference which is in g/mol)
cc
c     Variables (vector quantities positive downwards):
c dkt(i)=temperature gradient of diffusivity (m2/s/C)
c dro(i)=depth gradient of density (kg/m3/m)
c ro(i)=density (kg/m3)
c sg(i)=load (kg/m2)
c t(i)=temperature (C)
c tnew(i)=temperatures at new time step
cc ta(ko)=time of temperature history record (a)
cc tt(ko)=temperature history at lock-in depth (C)
cc ts(ko)=temperature history at surface (C)
c z(i)=depth below surface (m)
c
c aa=annual average temperature at surface
c a=flow-law parameter (Pa-3s-2)
c bdot=accumulation rate of ice (m/s)
c cp=heat capacity (J/kg/C)
c dct=gradient of c with t (J/kg/C2)
c dkro=gradient of diffusivity with density (m2/s/kg/m3)
c dz=depth step (m)
c geotherm= initial temperature profile (C/m)
c g=gravitational acceleration (m/s2)
c id=number of grid points in array
c jprint=test whether output this loop
c ntstep=number of time steps
c ntprint= number of times during experiment to print depth profile output
c qs=strain heating (C/s)
c rh=ice thickness (m)
c rhoi=ice density (kg/m3)
c rk=thermal diffusivity of firn (m2/s)
c spa=seconds per annum
c time=point in experiment (s)
c tmax=experiment length; read in as a and recalc as s
c tstep=time step; read in as steps/year and recalc as length in s
c tst: also time step in steps/year
c vz=vertical velocity (m/s)
c w1,w2,w3,w4,w5,w6=computational intermediates
c
c tinit(tlevels) = save the initial temperature profile t(Z)
c aainit: save the initial surface temperature.
c zi(tlevels) = depth of the initial temperature profile
c tprofile, tprofile2 and tprofileini are used to store data for file .42 and 43
c critdepth1, first0, first1, first2 used in the parameterisation of the density profile
c critdepth2, second0, second1, second2  used in the parameterisation of the density profile

c -- indexes ---
c jt : time step index
c bb: basis function index
c nfin and fin are intermediates for the file naming to work

c -- file naming ---
c there is an automatic file naming for all the basis functions. 
c This code uses Fourier series, so file names will have sin or cos, 
c followed by p, followed by 10000+ the index of the frequency number, bb.
c firn_sin_p_10000.32-43 is the output of the model without any added frequencies,
c that can be used a a simple model output, and as the baseline for all the sin/cos
c perturbations.
c .32 has an output every 10 steps
c .35 is the time series of surface temp with the annual average
c .33 is the temperature profile.
c .42 and .43 are a series of temperature profile, for diagnostic
c look at the bottom of the code for details.
c - all files are included in a folder given as the input 'ofolder'.
c
c --- to run the code, install gfortran and run in the command line:
c      gfortran  -O3 firn_BT.f -o firn_BTg
c     ./firn_BTg
c
c --------------------------------------------------------------------------
      program firn
      implicit real*8(a-h,o-z)
      integer level1, tlevels, tlevels1, dimlevels, layers
      parameter ( level1=501, tlevels = 621, tlevels1 = tlevels+1)
      PARAMETER (dimlevels = 65, layers = 500)
      real dkt(tlevels), dro(tlevels), ro(tlevels), sg(tlevels)
      real  t(tlevels), tnew(tlevels), z(tlevels), tinit(tlevels)
      real  deriv(tlevels), c(tlevels),diff(tlevels),s(tlevels)
      real tprofile(tlevels1, 1002), zi(tlevels)
      real tprofile2(tlevels1, 1002), tprofileini(tlevels1, 1002)

      integer jt, tst, agetst, bb
      integer duration
      real ny, amp, aaa, pi, to
      integer age, indname, acind, age2
       character*80 fin, fin32, fin42, fin33,fin34, fin35, fin43
       character*80 fin49, fileiniT
      real indf(402), indcos(402) 
      real indsign(402)
      real act(230), acr(230), acstart
      integer charnum, nfin, bmax, ntst, itst, period, t42freq
      character*80 ofolder, ofile
      real dz, dz2, aainit, bdot

c
c initialize
      open (26,file='inputfiles/inputdata_MegaWAIS_11wsc.txt'
     $          ,status='old')
      read (26,*) dz, dz2, bdot, rh, aa, geotherm
      read (26,*) duration,tstep,ntprint,pressure,ustart,uscale
      read (26,*) critdepth1, first0, first1, first2
      read (26,*) critdepth2, second0, second1, second2
      read(26,'(A80)') ofolder
      read(26,'(A80)') ofile
      read(26,'(A80)') fileiniT
      read(26,'(I3)') bmax
      close(26)

      print*, 'I read the file'
c calculate real ice density (Bader)
      rhoice = 916.5 - 0.14438*aa - 0.00015175*aa**2
c      print*, 'rhoice = ', rhoice
      spa = 3600.*24.*365.25
      tmax = duration*spa
      tst=tstep
      tstep = spa/tstep
      ko = 1
      g = 9.82
      R = 8.314
      sg(1) = 0.
      dct = 7.755
      ntstep = tmax/tstep
      kprint = ntstep/(ntprint-1)     
      timestep = tstep
      levels = dimlevels - 1
      idl = tlevels - 1
      aaK = aa + 273.15
      pi = 3.1415927
c get the folder name for output files
c ------ get the file depending on sin/cos -----
      CALL system('mkdir ' // ofolder )
c find length of folder string ignoring trailing blanks:
      DO 15, I = LEN(ofolder), 1, -1 
         IF(ofolder(I:I) .NE. ' ') GO TO 20 
15    CONTINUE 
20    charnum = I 
      nfin=charnum+17
      fin = ofolder(1:charnum)


c calculate depths and densities in open firn column using quadratic fit to data:
       rhoclose = second0 + second1*critdepth2 + second2*critdepth2**2	
      do 40 i=1,level1
      z(i) = dz*(i-1)
      if (z(i).lt.critdepth1) ro(i)=first0+ first1*z(i)+ first2*z(i)**2
      if (z(i).ge.critdepth1) ro(i)=second0+second1*z(i)+second2*z(i)**2		
      if (z(i).ge.critdepth2) ro(i) = rhoice-(rhoice-rhoclose)
     $  *exp(-(z(i)-critdepth2)/(rhoice-rhoclose)*(second1+2*second2
     $     *critdepth2))
  					
  40	continue
        do 41 i=level1+1, tlevels
           z(i) = z(level1)+ dz2*(i-level1)
            ro(i) = rhoice-(rhoice-rhoclose)
     $  *exp(-(z(i)-critdepth2)/(rhoice-rhoclose)*(second1+2*second2
     $     *critdepth2))
 41        continue
           print*, 'dz=', dz, 'until', z(level1),', then dz=', dz2, 
     $          'until', z(tlevels)
c        print*, ro(1), ro(10), ro(30), ro(100), ro(300)
c ****************************************************
c calculate density gradients
        do 50 i=2,level1-1
         dro(i)=(ro(i+1)-ro(i-1))/(2.*dz)
	 sg(i)=sg(i-1)+ro(i)*dz
	 dkt(i)=1.084d-16*ro(i)**3-1.88d-13*ro(i)**2
     $  +8.291d-11*ro(i)-1.569d-8
  50   continue
       i=level1
       dro(i)=(ro(i+1)-ro(i-1))/(dz+dz2)
	 sg(i)=sg(i-1)+ro(i)*dz
	 dkt(i)=1.084d-16*ro(i)**3-1.88d-13*ro(i)**2
     $  +8.291d-11*ro(i)-1.569d-8

        do 51 i=level1+1,tlevels-1
         dro(i)=(ro(i+1)-ro(i-1))/(2.*dz2)
	 sg(i)=sg(i-1)+ro(i)*dz2
	 dkt(i)=1.084d-16*ro(i)**3-1.88d-13*ro(i)**2
     $  +8.291d-11*ro(i)-1.569d-8
 51   continue

        sg(tlevels)=sg(tlevels-1)+ro(tlevels)*dz
	dkt(tlevels)=1.084d-16*ro(tlevels)**3-1.88d-13*ro(tlevels)**2
     $  +8.291d-11*ro(tlevels)-1.569d-8

        fin = ofolder(1:charnum) // "/Initialprofile.txt"
        open(52, file=fin, status='replace')
        do 53 i=1,tlevels
        write(52, 521) z(i), ro(i), dro(i), sg(i), dkt(i)
 521    format(1x, f10.5, 1x, f10.5, 1x, g10.3, 1x, g10.3, 1x, g10.3)
 53   continue
      close(52)
c-------------------------------------------------------------
c ----  initialize temperatures with stationary profile ------
        open(27, file = fileiniT, status = 'old')
     	do 58 i=1,tlevels
           read(27,*), zi(i), t(i)
 58     continue
        close(27)
c        print*, 'I calculated density and inital T'
c ---- find a stationary  temperature profile:
c$$$              do 951 jt = 1,5000*tst
c$$$                 time = jt*tstep
c$$$                 ta=time/spa
c$$$              t(1)=aa+10.*(cos(2.*3.1415926*(ta-int(ta)))+
c$$$     $         0.3*cos(4.*3.1415926*(ta-int(ta))))
c$$$
c$$$      do 304 i = 2, level1-1
c$$$       dt = (t(i+1)-t(i-1))/2/dz
c$$$       d2t = (t(i+1)-2.*t(i) + t(i-1))/dz/dz
c$$$       rr = ro(i)
c$$$       rk = (1-0.00882*(t(i)+30))*(-1.229d-14*rr**3
c$$$     $  + 2.1312d-11*rr**2-9.4d-9*rr + 1.779d-6)
c$$$       dkro = (1-0.00882*(t(i)+30))*(-3.687d-14*rr**2
c$$$     $  + 4.2624d-11*rr - 9.4d-9)
c$$$       a = 4.26d-13*exp(-7217./(273.15+t(i)))      
c$$$       cp = 2096.*(1+0.0037*t(i))
c$$$       qs = 2.*(bdot/rh)**(4./3.)/(rr*cp*a**(1./3.))
c$$$       vz = rhoice*bdot/rr*(1-sg(i)/rhoice/rh)
c$$$       a1 = (rk/rr+dkro)*dro(i)-vz    
c$$$       a2 = dkt(i) + rk*dct/cp
c$$$       a3 = sg(i)*bdot*rhoice/rr**3*dro(i)*g/cp + qs
c$$$       tnew(i) = t(i) + tstep*(rk*d2t+a1*dt+a2*dt*dt+a3)
c$$$ 304   continue
c$$$
c$$$      i=level1
c$$$      dt = (t(i+1)-t(i-1))/(dz+dz2)
c$$$       d2t = ((t(i+1)-t(i))/dz2-(t(i)-t(i-1))/dz)/((dz+dz2)/2)
c$$$       rr = ro(i)
c$$$       rk = (1-0.00882*(t(i)+30))*(-1.229d-14*rr**3
c$$$     $  + 2.1312d-11*rr**2-9.4d-9*rr + 1.779d-6)
c$$$       dkro = (1-0.00882*(t(i)+30))*(-3.687d-14*rr**2
c$$$     $  + 4.2624d-11*rr - 9.4d-9)
c$$$       a = 4.26d-13*exp(-7217./(273.15+t(i)))      
c$$$       cp = 2096.*(1+0.0037*t(i))
c$$$       qs = 2.*(bdot/rh)**(4./3.)/(rr*cp*a**(1./3.))
c$$$       vz = rhoice*bdot/rr*(1-sg(i)/rhoice/rh)
c$$$       a1 = (rk/rr+dkro)*dro(i)-vz    
c$$$       a2 = dkt(i) + rk*dct/cp
c$$$       a3 = sg(i)*bdot*rhoice/rr**3*dro(i)*g/cp + qs
c$$$       tnew(i) = t(i) + tstep*(rk*d2t+a1*dt+a2*dt*dt+a3)
c$$$
c$$$      do 303 i = level1+1, tlevels-1
c$$$       dt = (t(i+1)-t(i-1))/2/dz2
c$$$       d2t = (t(i+1)-2.*t(i) + t(i-1))/dz2/dz2
c$$$       rr = ro(i)
c$$$       rk = (1-0.00882*(t(i)+30))*(-1.229d-14*rr**3
c$$$     $  + 2.1312d-11*rr**2-9.4d-9*rr + 1.779d-6)
c$$$       dkro = (1-0.00882*(t(i)+30))*(-3.687d-14*rr**2
c$$$     $  + 4.2624d-11*rr - 9.4d-9)
c$$$       a = 4.26d-13*exp(-7217./(273.15+t(i)))      
c$$$       cp = 2096.*(1+0.0037*t(i))
c$$$       qs = 2.*(bdot/rh)**(4./3.)/(rr*cp*a**(1./3.))
c$$$       vz = rhoice*bdot/rr*(1-sg(i)/rhoice/rh)
c$$$       a1 = (rk/rr+dkro)*dro(i)-vz    
c$$$       a2 = dkt(i) + rk*dct/cp
c$$$       a3 = sg(i)*bdot*rhoice/rr**3*dro(i)*g/cp + qs
c$$$       tnew(i) = t(i) + tstep*(rk*d2t+a1*dt+a2*dt*dt+a3)
c$$$ 303   continue
c$$$       if(mod(ta, 500.).eq.0) write(*,'(1X,A3,F7.1,1x,A2,1x,
c$$$     $      G12.5,1x,A4,1x,G12.5,1x,A6,1x,3G12.5)') 
c$$$     $       'ta', ta,'aa', aa, 't(1)', t(1),
c$$$     $      't(100)', t(100), t(tlevels-1), t(tlevels)
c$$$c assign new temperature and concentration
c$$$       do 551 i = 2, tlevels-2
c$$$        t(i) = tnew(i)
c$$$ 551  continue
c$$$ 951  continue

c save the stationary profile:
        fin = ofolder(1:charnum) // "/StationaryTprofile.txt"
      open(54, file=fin, status='replace')
        do 59 i=1,tlevels
        write(54, 522) z(i), t(i)
 522    format(1x, f10.5, 1x, f10.5)
        tinit(i)=t(i)
 59   continue
      close(54)
      aainit=aa

c       add a number of the frequency we want to calculate on the prompt
c     length of indfreq0 is 204, indfreq2 and 4 are 202, indfreq3 is 42
c	open(28,file='indfreq0.txt', status='old') 

	open(28,file='inputfiles/indfreq5.txt', status='old') 
	do 56 i=1,402
		 read(28,*) indf(i), indcos(i), indsign(i)
 56           continue 
        close(28)

c	============================================================
c	============================================================
c	loop over all the harmonics we are considering
c	indf(bb) index is for the frequency number 
c	indcos(bb)is 1 for  sin or 2 for cos
c	indsign(bb) is for plus or minus
c              bb=1
	do 1001 bb =1,bmax
c              aa=aaa
c ****************************************************
c -- to save on run time, we can change itst to start the run later
c in this configuration this section is not useful.
c - figure out the length of this run:
        ntst=ntstep
        itst=duration
        if (indf(bb).gt.1) then
c     get 20 periods: length=itst
c        itst=min(int(4000*20/indf(bb)),duration)
c        print* ,'itst', itst
c time of departure: duration-length, agetst used for the initial profile
        agetst = int((duration-itst)/10)
c        print*, 'agetst', agetst
c     get a power of 10 for the accurate initial profile
c        itst=(int(duration/10)-agetst)*10
c     number of timesteps
        ntst=itst*tst
c        print*,'duration=', itst, ntst
        end if

c  re-initialise temperature profile:
              aa=aainit
              if(bb.ge.2) then
              do 65 i=1,tlevels-1
                 t(i)=tprofileini(i+1, agetst+2)
 65              continue
c                 print*, 'reinitialised temp for agetst', 
c     $                tprofileini(1,agetst+2)
              end if

c      ==============================================================
c ----  insert history of temp ----- 
        open(29, file=ofile , status='old')

c re-initisalise the time of file reading
        to=0.   
c     start at the right point for duration=itst
c 66        CONTINUE
          read(29,*), to, aa
c           IF (to.lt.duration-itst) GOTO 66
        acind=1
c        print*, 'read temperature history', ofile
c        PRINT*, 'start date', to, '+', itst, '=', duration
c       ========================================================
c       main timestepping loop
c       --------------------------------------------------------
c
c 		Top boundary condition simulation
c
c           do 1000 jt = 1,ntstep
           do 1000 jt = 1,ntst
       time = jt*tstep
       if(jt.eq.1) print*, 'Im in the loop'
c*********************************************************
c  calculate surface temperature over time (top boundary condition, model forcing)
      ta=time/spa
	  day=(ta-aint(ta))*365.25
c read anual temp from file
c          if(jt.gt.100*tst.and.int(jt/tst)*tst.eq.jt) then
c          if(int(jt/tst)*tst.eq.jt) then
          if(ta.ge.to+1) then
             read(29,*), to, aa
          end if


c     --- seasonal cycle ---
c     -either 1. approximate generic seasonal cycle (day 0 is Jan 0)
c$$$	   t(1) = aa + 13.09
c$$$	  if (day.gt.30..and.day.le.95.)
c$$$     $          t(1) = (aa+13.09) - 19.7/65.*(day-30)
c$$$	  if (day.gt.95..and.day.le.260.) t(1) = aa - 6.61
c$$$	  if (day.gt.260..and.day.le.350.) 
c$$$     $         t(1)=(aa-6.61)+19.7/90.*(day-260)
c     -or 2. sinusoidal seasonal cycle:
      season=10.*(cos(2.*3.1415926*(ta-int(ta)))+
     $         0.3*cos(4.*3.1415926*(ta-int(ta))))


c     --- perturbation to temperature --- 
      amp=1
c  either 1. triangles
c      ny=25.
c      if (ta.ge.100..and.ta.lt.(100+ny)) t(1) = t(1)+amp*(ta-100)/ny
c      if (ta.ge.(100+ny).and.ta.lt.(100+2*ny)) 
c     $     t(1) = t(1)+amp*(1-(ta-100-ny)/ny)

c or 2. add a step function for 300 years:
c          if(ta.le.100) aaa = aa
c          if (ta.gt.100.and.ta.le.300) aaa = aa+1
c          if (ta.gt.300) aaa = aa

c or 3. sin/cos
c change indsign(bb) to make it positive when it gets to the end:
c this is not necessary, but it simplifies the plotting later
c make sure that the matlab file takes this into account.
      if(mod(indf(bb),2.).eq.0) then
         indsign(bb)=1
         else
           indsign(bb)=-1
            end if

c now add the sin/cos perturbation
	if (indcos(bb).eq.1) then
	  t(1)=aa+season+indsign(bb)*amp*sin(2.*pi*(indf(bb)/4000.)*ta)	   
	  else 
	  t(1)=aa+season+indsign(bb)*amp*cos(2.*pi*(indf(bb)/4000.)*ta)
	  if (jt.eq.10) print *, "i changed t1 this cos time"
         if (jt.eq.10) 
     $         print *, indsign(bb)*amp*cos(2*pi*(indf(bb)/4000.)*ta)
	end if   
  	 if (jt.eq.10*tst) then
           write(*,'(1X,A3,F4.1,1x,A5,1x,G11.4,1x,G11.4,1x,G11.4)') 
     $           'ta=',ta, " arg=",(indf(bb)/4000.)*ta, 
     $           indsign(bb)*amp,
     $           indsign(bb)*amp*sin(2.*pi*(indf(bb)/4000.)*ta) 
         end if

c$$$ 	if (indcos(bb).eq.1) then
c$$$	  t(1)=t(1) - indsign(bb)*amp*sin(2.*pi*(1./indf(bb))*ta)	   
c$$$	  else 
c$$$	  t(1)=t(1) - indsign(bb)*amp*cos(2.*pi*(1./indf(bb))*ta)
c$$$	  if (jt.eq.10) print *, "i changed t1 this cos time"
c$$$         if (jt.eq.10) 
c$$$     $         print *, indsign(bb)*amp*cos(2*pi*(1./indf(bb))*ta)
c$$$	end if   
c$$$  	 if (jt.eq.10*tst) then
c$$$           write(*,'(1X,A3,F4.1,1x,A5,1x,G10.4,1x,G10.4,1x,G10.4)') 
c$$$     $           'ta=',ta, " arg=",(1./indf(bb))*ta, 
c$$$     $           indsign(bb)*amp,
c$$$     $           indsign(bb)*amp*sin(2.*pi*(1./indf(bb))*ta) 
c$$$         end if


c      --- Temperature model (courtesy of Richard Alley; Alley and Koci, 1989): ---
c       if (int(jt/100000)*100000.eq.jt) print*, jt,ta,  t(1)
      do 301 i = 2, level1-1
       dt = (t(i+1)-t(i-1))/2/dz
       d2t = (t(i+1)-2.*t(i) + t(i-1))/dz/dz
       rr = ro(i)
       rk = (1-0.00882*(t(i)+30))*(-1.229d-14*rr**3
     $  + 2.1312d-11*rr**2-9.4d-9*rr + 1.779d-6)
       dkro = (1-0.00882*(t(i)+30))*(-3.687d-14*rr**2
     $  + 4.2624d-11*rr - 9.4d-9)
       a = 4.26d-13*exp(-7217./(273.15+t(i)))      
       cp = 2096.*(1+0.0037*t(i))
       qs = 2.*(bdot/rh)**(4./3.)/(rr*cp*a**(1./3.))
       vz = rhoice*bdot/rr*(1-sg(i)/rhoice/rh)
       a1 = (rk/rr+dkro)*dro(i)-vz    
       a2 = dkt(i) + rk*dct/cp
       a3 = sg(i)*bdot*rhoice/rr**3*dro(i)*g/cp + qs
       tnew(i) = t(i) + tstep*(rk*d2t+a1*dt+a2*dt*dt+a3)

  301 continue

      i=level1
      dt = (t(i+1)-t(i-1))/(dz+dz2)
       d2t = ((t(i+1)-t(i))/dz2-(t(i)-t(i-1))/dz)/((dz+dz2)/2)
       rr = ro(i)
       rk = (1-0.00882*(t(i)+30))*(-1.229d-14*rr**3
     $  + 2.1312d-11*rr**2-9.4d-9*rr + 1.779d-6)
       dkro = (1-0.00882*(t(i)+30))*(-3.687d-14*rr**2
     $  + 4.2624d-11*rr - 9.4d-9)
       a = 4.26d-13*exp(-7217./(273.15+t(i)))      
       cp = 2096.*(1+0.0037*t(i))
       qs = 2.*(bdot/rh)**(4./3.)/(rr*cp*a**(1./3.))
       vz = rhoice*bdot/rr*(1-sg(i)/rhoice/rh)
       a1 = (rk/rr+dkro)*dro(i)-vz    
       a2 = dkt(i) + rk*dct/cp
       a3 = sg(i)*bdot*rhoice/rr**3*dro(i)*g/cp + qs
       tnew(i) = t(i) + tstep*(rk*d2t+a1*dt+a2*dt*dt+a3)

      do 302 i = level1+1, tlevels-1
       dt = (t(i+1)-t(i-1))/2/dz2
       d2t = (t(i+1)-2.*t(i) + t(i-1))/dz2/dz2
       rr = ro(i)
       rk = (1-0.00882*(t(i)+30))*(-1.229d-14*rr**3
     $  + 2.1312d-11*rr**2-9.4d-9*rr + 1.779d-6)
       dkro = (1-0.00882*(t(i)+30))*(-3.687d-14*rr**2
     $  + 4.2624d-11*rr - 9.4d-9)
       a = 4.26d-13*exp(-7217./(273.15+t(i)))      
       cp = 2096.*(1+0.0037*t(i))
       qs = 2.*(bdot/rh)**(4./3.)/(rr*cp*a**(1./3.))
       vz = rhoice*bdot/rr*(1-sg(i)/rhoice/rh)
       a1 = (rk/rr+dkro)*dro(i)-vz    
       a2 = dkt(i) + rk*dct/cp
       a3 = sg(i)*bdot*rhoice/rr**3*dro(i)*g/cp + qs
       tnew(i) = t(i) + tstep*(rk*d2t+a1*dt+a2*dt*dt+a3)

 302  continue

c     =======================================================
c     assign new temperatures and new concentrations
c       -----------------------------------------------------
c		Temperatures:
c
       do 550 i = 2, tlevels-2

        t(i) = tnew(i)

  550  continue

c      ==============================================================
c       Output sampling routines for the open firn: 
c ------ get the file depending on sin/cos -----
       if (jt.eq.1) then
           if (indcos(bb).eq.2)
     $       fin = ofolder(1:charnum) // "/firn_cos_p_10000"
           if (indcos(bb).eq.1)
     $       fin = ofolder(1:charnum) // "/firn_sin_p_10000"
	indname = 10000 + int(indf(bb))	  
        write(fin(nfin-4:nfin), '(I5)') indname
           fin32 = fin(1:nfin) // ".32"
           fin33 = fin(1:nfin) // ".33"
           fin34 = fin(1:nfin) // ".34"
           fin35 = fin(1:nfin) // ".35"
           fin42 = fin(1:nfin) // ".42"
           fin43 = fin(1:nfin) // ".43"
           fin49 = fin(1:nfin) // ".log"
           write(*, '(A40)') fin32
c       if(jt.eq.1) then
      open(32, file=fin32, status = 'replace')
      open(33, file=fin33, status = 'replace')
      open(34, file=fin34, status = 'replace')
      open(35, file=fin35, status = 'replace')
      open(42, file=fin42, status = 'replace')
      open(43, file=fin43, status = 'replace')
      open(49, file=fin49, status = 'replace')
c      print*, fin32, bb, '/204'
      end if
c--- end of getting file name and opening the file, if(jt.eq.1)... ---

c  .32: write a time series at the surface and lock-in depth:
	  if (ta.gt.samplestart) then
	    if (int(jt/20)*20.eq.jt) then	
c             if(ta.le.samplestart+10) then
	  write (32,720) ta,t(1),t(40),
     $              t(100),t(200),t(300),t(tlevels)
 720   format(2x,f10.4,2x,f10.5,2x,2x,f10.5,2x,
     $         f10.5,2x,f10.5,2x,f10.5,2x,f10.5)      
        end if

	  end if
c

c .33: write a depth profile at the end of the run, 
c write a depth profile at the end of the run, on dec 31, 
          if (jt.eq.(ntst-tst)) then
      do 900 jj=1,tlevels
        write (33,800) z(jj), t(jj), ro(jj)
 800    format (2x,f7.2,2x,f10.5,2x,f10.4,2x,f10.5,2x,f10.7,2x,f10.5)
 900    continue
        end if

c .34: write a depth profile just before the end of the run,
c       if (jt.eq.(ntstep-2*tst)) then
c        if(ta.eq.int((duration-1)/indf(bb))*indf(bb)) then
        if(ta.eq.int((itst-1)/indf(bb))*indf(bb)) then

        do 901 jj=1,tlevels
          write (34,801) z(jj), t(jj), ro(jj)
 801      format (2x,f7.2,2x,f10.5,2x,f10.4,2x,f10.5,2x,f10.7,2x,f10.5)
 901    continue
       end if

c .35: write the annual average surface temperature	
	    if (int(ta).eq.ta) then		
	  write (35,720) ta,aa,t(40),
     $              t(100),t(200),t(300),t(tlevels)
        end if


c     put interesting profiles into a matrix .42
        if (jt.eq.1) then
           tprofile(1,1)=0
           tprofile(1,2) = ta
           tprofile2(1,1)=0
           tprofile2(1,2) = ta
           do 310 jj = 1, tlevels
              tprofile (jj+1,1) = z(jj)
              tprofile (jj+1,2) = t(jj)
              tprofile2 (jj+1,1) = z(jj)
              tprofile2 (jj+1,2) = t(jj)

c              tprofile (jj+1,2) = tav(jj)
 310       continue
           age = 1
           age2=1
        end if

c decide what time of year to print: 
c Actually, with the cos, 0=dec 21 pretty much...
c        age = int(ta/10.)
c        if(int(ta*10)/10.eq.ta) then
c           if (int(ta/100)*100.eq.ta) then
        t42freq=10
        age = int(ta/t42freq)
           if (int(ta/1000)*1000.eq.ta) then
              write(*,'(A4,1x,I4,1x,G11.4,1x,G11.4,1x,2G15.8)')
     $             'age',age, ta,jt,t(1), t(100)
           end if

        if(age*t42freq.eq.ta) then
           tprofile(1, age+2) = ta
           do 311 jj = 1, tlevels
              Tprofile (jj+1,age+2) = t(jj)
c              tprofile (jj+1,age+2) = tav(jj)
311        continue
c          print*, 'checked file 42: fill', age
         end if

         if (indf(bb).ge.1) then
            period=int(2000/indf(bb))
            if(int(ta/period)*period.eq.ta) then
              tprofile2(1, age2+2) = ta
           do 313 jj = 1, tlevels
              Tprofile2 (jj+1,age2+2) = t(jj)
c              tprofile (jj+1,age+2) = tav(jj)
 313       continue
              age2=age2+1
c          print*, 'checked file 43: fill', age2
            end if
         if(jt.eq.1) print*,'period=', period
         end if
 
         if(jt.eq.(ntst-1)) then
            print*, 'Ill write the file soon', age, duration, ta
            do 312 jj = 1,tlevels
                write(42,'(1X,200G18.8)') (tprofile(jj,nn),nn=1,
     $              min(age+2,100))
 312        continue
            do 314 jj = 1,tlevels
                write(43,'(1X,200G18.8)') (tprofile2(jj,nn),nn=1,age2+2)
 314         continue

             if(bb.eq.1) then
                do 315 jj=1,tlevels+1
                   do 316 nn=1, age+2
                      tprofileini(jj,nn)=tprofile(jj,nn)
 316               continue
 315            continue
            end if

	    write(49,*), 'output of firnJ11long03b4.f'
	    write(49,*), 'inputdata_MegaWAIIS_11wsc.txt'
	    write(49,*), fileiniT
c            write(49,*), 'output of constant profile'
           write(49,*), 'output of ', ofile
	    write(49,*), 'dz = ', dz
	    write(49,*), 'aa = ', aaa
	    write(49,*), 'duration = ', itst, duration,  ' years'
	    write(49,*), 'tstep = ', tst, ' times per year'
	    write(49,*), 'amp = ', amp, 'degC amplitude'
         write(49,*),'sin/cos output of the WSD borehole with '
         write(49,*),'2 cos real seasonal cycle'
         end if
c          print*, 'checked log'


c 			go to 1200
	

 1000   CONTINUE
        close(29)
 1001    continue

c       ------------------------------------------------------
c       end of main time stepping loop
c       =======================================================

 1200	end
