/* file:  ComputeTemperaturePrecipitation.cpp
   Compute the temperature and precipitation at time t from 
   the data at present day and delta18O
 */

#include "./elements.h"

void ComputeDelta18oTemperaturePrecipitation(IssmDouble Delta18oSurfacePresent, IssmDouble Delta18oSurfaceLgm, IssmDouble Delta18oSurfaceTime,
				     IssmDouble Delta18oPresent, IssmDouble Delta18oLgm, IssmDouble Delta18oTime, 
				     IssmDouble* PrecipitationsPresentday,
				     IssmDouble* TemperaturesLgm, IssmDouble* TemperaturesPresentday, 
				     IssmDouble* monthlytemperaturesout, IssmDouble* monthlyprecout){

  IssmDouble monthlytemperaturestmp[12],monthlyprectmp[12];
  IssmDouble delta18oLapseRate=-6.2*pow(10.,-3);
  IssmDouble glacialindex; // used to vary present day temperature
  int imonth=0;
  // printf("Delta18oTime %f\n",Delta18oTime);
  // printf("Delta18oPresent %f\n",Delta18oPresent);
  // printf("Delta18oLgm %f\n",Delta18oLgm);

  glacialindex = (Delta18oTime-Delta18oPresent-delta18oLapseRate*(Delta18oSurfaceTime-Delta18oSurfacePresent))
    /(Delta18oLgm-Delta18oPresent-delta18oLapseRate*(Delta18oSurfaceLgm-Delta18oSurfacePresent));
  // printf(" monthlytemperatureseeeeee %f\n",TemperaturesPresentday[1]);
  // printf(" monthlypreceeeeeeee %f\n",TemperaturesLgm[1]);
  // printf(" glacialindex %f\n",glacialindex);

  for (int imonth = 0; imonth<12; imonth++){  
    monthlytemperaturestmp[imonth] = TemperaturesPresentday[imonth];//glacialindex*TemperaturesLgm[imonth] + (1-glacialindex)*TemperaturesPresentday[imonth];
    //monthlyprectmp[imonth] = 1.5*pow(2,((monthlytemperaturestmp[imonth]-273.15-0)/10)); //equation from rob's paper
    monthlyprectmp[imonth] = PrecipitationsPresentday[imonth];

    /*Assign output pointer*/
    *(monthlytemperaturesout+imonth) = monthlytemperaturestmp[imonth];
    *(monthlyprecout+imonth) = monthlyprectmp[imonth];
  }
  // printf(" tempera %f\n",monthlytemperaturestmp[1]);
}
