/* Simples Energiebilanzmodell (EBM)

aus: v.Storch/Guess/Heimann "Klimamodellierung"

S. Guess, 1997
 
Modifikationen: 
   - atmosphaerische Transmissivitaet = 0.225
   - temperaturabhaengige albedo als stueckweise lineare funktion
     mit Eis-albedo von 0.5 fuer T<273.15 und a=0.1 fuer T>303.15
   - Schalter fuer Ein/Ausschaltung der nicht-linearen Feedbacks
     und des stochastischen Antriebs fuer die atm. Transmissivitaet
     
mh, Thu, Jan 11, 2001 */

#include <stdio.h>
#include <stdlib.h>  
#include <math.h>

/* Maximal-Anzahl Zeitschritte: */ 
#define NPT 40000

float zufall01(long *idum) ;
float zufallgauss(long *idum) ;

main( void ) 
{
   /* DEKLARATIONEN: */ 
      /* Bodennahe Temperatur abhaengig von Zeit : */
      double temp[NPT+1] , t_end ;
      /* Strahlungsfluesse f : */
      double f_kw_in , f_lw_out , f_netto ; 
      /* Kurzwellige Albedo und langwellige Transmissivitaet : */
      double albedo , transmiss ;
      /* Stochastische Stoerung */ 
      double stoer[NPT+1] ;

      int i , i_dt_tage , n_jahre , n_zeitschritte ;
      int j_albedo_feedback , j_stoer_transmiss ; 
      long *p_stoer ;
      long i_stoer ;  
      double ckl , ckd , ckt , cko ;
      double albedo_0 , transmiss_0 , f_kw_in_0 ;
      double albedo_min, albedo_max, temp_min, temp_max ;
      double sigma , stdabw_stoer ;
      double dt , ck , sim_time ; 
      double pi; 

      /* Ausgabe-File : */ 
      FILE *checkoutputfile1; 
      char szoutputfilename1[FILENAME_MAX], *psz_strchr_out ;  

   /* Schalter und Kontrollwerte: */

     /* Schalter fuer Temperatur-Albedo-Ruekkopplung:1=ein,0=aus:*/ 
     j_albedo_feedback = 1 ; 
     /* Definition der Temperatur-Albedo-Ruekkopplung */
     temp_min = 273.15 ;
     temp_max = 303.15 ;
     albedo_max = 0.5 ;
     albedo_min = 0.1 ;
     
     /* Schalter fuer Stoerungen in Transmissivitaet:1=ein,0=aus:*/ 
     j_stoer_transmiss = 0 ; 
     /* Standardabweichung der Stoerungen in der Transmissivitaet: */
     stdabw_stoer = 0.03 ; 

     /* Anfangswert fuer Temperatur [K]: */
     temp[0] = 290.0 ;

     /* Zeitschrittweite in Tagen : */
     i_dt_tage = 10 ;

     /* Anzahl Simulationsjahre : */
     n_jahre = 1000 ; 
     n_zeitschritte = (int)((365./(float)(i_dt_tage)) * n_jahre ) ; 


   /* KONSTANTEN: */  

     /* Stefan-Boltzmann-Konstante in W/(m**2 K**4): */
     sigma = 5.67E-8 ;
     /* Pi: */ 
     pi = 4.0 * atan(1.0) ; 

     /* Waermekapazitaet:   Nur Atmosphaere:*/
     ckl = 1.0E7 ; 
     /* Waermekapazitaet: Ozean.Deckschicht 50m (real ca.70m):*/
     ckd = 2.0E8 ;
     /* Waermekapazitaet: Oberer Ozean: Schicht 250m (real ca.360m)*/
     ckt = 1.0E9 ;
     /* Waermekapazitaet: Ozean, 3900m, 70% der Erdoberflaeche */
     cko = 1.6E10 ;
    /* Auswahl des Kompartiments der Waermekapazitaet: */
     ck = ckd ;

     /* Albedo Mittelwert constant : */ 
     albedo_0 = 0.30 ;
     
     /* Transmissivitaet der Atmosphaere (ist so getunt, dass
        bei Standardwerten der Parameter (u.a. Albedo=0.3)
        eine Gleichgewichtstemperatur
        von 288.13K herauskommt: */
     transmiss_0 = 0.224889 ;
                                          
     /* Kurzwellige Einstrahlung (Solarkonstante) in W/m**2: */ 
     f_kw_in_0 = 342.0 ;

     /* Zeitschrittweite in Grundeinheit Sekunden: */
     dt = 60*60*24 * (float) i_dt_tage ;
 
     /* Zufallszahlengenerator initialisieren: */
     i_stoer = -14 ;
     p_stoer = &i_stoer ;

  /* AUSGABE-FILE OEFFNEN: */ 

     checkoutputfile1 = fopen( "ebmout.dat" , "w" );
     if( checkoutputfile1 == NULL)
     { 
      puts("Fehler: ebmout.dat kann nicht geschrieben werden \n");
      exit(EXIT_FAILURE);
     } ;
  

  /* VORWAERTS-INTEGRATION: ZEITSCHRITTE DURCHLAUFEN: */

     for( i = 0; i <= n_zeitschritte ; i = i+1 )
      {

    /* Simulationszeit (Jahre) */

	 sim_time = dt * (float) i /(365*24*60*60) ;
					  
					  
    /* Normalverteilte Pseudo-Zufallszahlen: */ 
         stoer[i] = stdabw_stoer * (double) zufallgauss( &i_stoer ) ;

         /* Parameter: */

         if (j_albedo_feedback < 1 ) {
             albedo = albedo_0;
         } else {                         
             albedo = albedo_min * (temp[i]-temp_min)/(temp_max-temp_min)
                    + albedo_max * (temp_max-temp[i])/(temp_max-temp_min);
         if ( temp[i]<temp_min ) {albedo=albedo_max;};
         if ( temp[i]>temp_max ) {albedo=albedo_min;};
         };
         
         transmiss = transmiss_0 
                + (float) j_stoer_transmiss * stoer[i] * transmiss_0 ;

         /* Strahlungsfluesse: */ 
         f_kw_in = f_kw_in_0 ; 

	 f_lw_out = sigma*temp[i]*temp[i]*temp[i]*temp[i] ;

         f_netto = f_kw_in*(1 - albedo) - f_lw_out*(1.0+transmiss)/2.0 ; 
 
         /* Temperatur fuer naechsten Zeitschritt: */ 
         temp[i+1] = temp[i]   + dt * (1/(ck))* f_netto ; 


         /* AUSGABE: */ 
         fprintf( checkoutputfile1 , "   %f   %f  %f  %f  \n", 
             sim_time, temp[i], albedo, transmiss); 
      }     
    
  /* ENDE: */ 
 
    i = i-1 ;
    fclose(checkoutputfile1);

}
/************************ ENDE MAIN ********************/


/* Gleichverteilte Zufallszahlen zwischen 0 und 1: *****/
float zufall01( long * i_stoer )
{  long lauf, lauf2;
   static long l1=0;
   static long l2[2836];
   float resultat;
   if ( *i_stoer <= 0 || !l1 ) 
   {  if ( ((-1)*(*i_stoer)) < 1 )
      { *i_stoer = 1; 
      }
      else 
      { *i_stoer = (-1)*(*i_stoer); 
      }; 
      for( lauf=2836+7 ; lauf >= 0 ; lauf-- ) 
      {  lauf2=(*i_stoer)/127773;
         *i_stoer=16807*(*i_stoer-lauf2*127773) - 2836*lauf2;
         if( *i_stoer < 0 )
         { *i_stoer += 2147483647; 
         }; 
         if( lauf < 2836 ) 
         { l2[lauf] = *i_stoer;
         }; 
      };
      l1=l2[0];
   }; 
   lauf2 = (*i_stoer)/127773;
   *i_stoer = 16807 * (*i_stoer-lauf2*127773) - 2836 * lauf2;
   if( *i_stoer < 0 ) 
   { *i_stoer += 2147483647;
   }; 
   lauf = l1/(1+(2147483647-1)/2836);
   l1 = l2[lauf];
   l2[lauf] = *i_stoer;
   if ( (resultat=(1.0/2147483647)*l1) > (1.0 - 1.2e-7) )
   { return (1.0 - 1.2e-7); 
   }
   else 
   { return resultat;
   }; 
}

/* Gauss(0;1)verteilte Zufallszahlen aus zufall01: *****/
#include <math.h>
float zufallgauss(long *i_stoer)
 {
  static long longintval=0;
  static float resultat;
  float fwurzel,betrag2; 
  float wert1,wert2;
  float zufall01(long *i_stoer);
  if  (longintval == 0) 
  {  do 
     {  wert1 = 2.0*zufall01( i_stoer ) - 1.0;
        wert2 = 2.0*zufall01( i_stoer ) - 1.0;
        betrag2=wert1*wert1+wert2*wert2;
     } 
     while (betrag2 >= 1.0 || betrag2 == 0.0);
     fwurzel=sqrt(-2.0*log(betrag2)/betrag2);
     resultat=wert1*fwurzel;
     longintval=1;
     return wert2*fwurzel;
  } 
  else 
  {  longintval=0;
     return resultat;
  };
}