!*****************************************************************!
!*                4C (FORSEE) Simulation Model                   *!
!*                                                               *!
!*                                                               *!
!*          contains:                                            *!
!*       s_t_prof      generates initial soil temp. profile      *!
!*       BTFOUR        TRICOF use to develope soil-surface-temp. *!
!*                                                               *!
!*                  Copyright (C) 1996-2018                      *!
!*     Potsdam Institute for Climate Impact Reserach (PIK)       *!
!*          Authors and contributors see AUTHOR file             *!
!*  This file is part of 4C and is licensed under BSD-2-Clause   *!
!*                   See LICENSE file or under:                  *!
!*     http://www.https://opensource.org/licenses/BSD-2-Clause   *!
!*                           Contact:                            *!
!*       https://gitlab.pik-potsdam.de/foresee/4C                *!
!*                                                               *!
!*****************************************************************!
SUBROUTINE s_t_prof

! Generation of initial soil temperature profile

use data_par
use data_soil
use data_soil_t
use data_simul

implicit none

integer i, ia, ie, k
real ath, dth, rfr, sn, uhf, u, vk, vh, fourterm  
real tcsu, hcsu

    ia = 1
    ie = 365
    TQ = 10.

      CALL BTFOUR
      tcsu = 0.
      hcsu = 0.
      rfr = 2. * pi / 365.  ! radial frequency; Radialfrequenz
      UHF=2.*pi/(IE-IA+1)
      u = uhf * it   

! calculation of temperature profile commonly from day 1 (it=1) set in data_soil;
! Temperaturprofil berechnen, standardmaessig it=1 (1. Tag) in data_soil gesetzt
      do i = 1, nlay
           tcsu = tcsu + t_cond(i)*thick(i)
           hcsu = hcsu + h_cap(i)*thick(i)   
           ath = tcsu / hcsu  ! for a weighted mean both values are divided by the depth (i), thus they cancel each other; fuer gewichtetes Mittel beide Werte durch depth(i) teilen ==> weggekuerzt
           DTH=SQRT(2*ATH/RFR)
           VH=mid(I)/DTH
           fourterm = 0.
           do k = 1, nk
              VK=VH*SQRT(K+0.)
              SN=FTA(K)*EXP(-VK)*SIN(U*K+FTO(K)-VK)
              fourterm = fourterm + SN
           enddo
          temps(i) = TQ + fourterm
          if (flag_dayout .eq. 3) write (2244, *) i, temps(i), mid(i), ath, dth, fourterm
      enddo

END  subroutine s_t_prof
 
!******************************************************************************

SUBROUTINE BTFOUR

! using TRICOF for a Fourier series development for ground surface temperature; 
! Fourierreihenentwicklung fuer Boden-Oberflaechen-Temperatur unter Nutzung von TRICOF

use data_climate
use data_par
use data_soil
use data_soil_t
use data_simul

implicit none

integer i, n, nt, nts, nf, nend, naf, no, ne, lf
real a0
real, dimension(184):: FA,FB

! set amount of auxiliary points NF for transformation; 
! Anzahl der Stuetzstellen NF fuer Transformation festlegen

      nend = 365
      naf = 1

      NT=NEND-NAF+1
      NTS=1
      NF=(NT+NTS-1)/NTS
      N=(NF-1)/2
      IF((2*N-NF+1) .LT. 0) THEN 
          NF=NF-1
          NT=(NF*NTS)-NTS+1
          NEND=NAF+NT-1
      ENDIF
      NE=1+(NF/2)
      NO=NE-2
      NK=NO
      
! calculation of auxiliary points; Stuetzstellen berechnen

      LF=NAF
      DO I=1,NF
          airtemp   = tp(lf,1)
          airtemp_1 = tp(lf-1,1)
          airtemp_2 = tp(lf-2,1)
          rad       = rd(lf,1)
          iday = lf
          call surf_t
          if (lf .eq. 1) temps(1) = temps_surf
          Four_sp(i) = temps_surf
          LF=LF+NTS
      ENDDO

! Fourier transformation; FOURIERTRANSFORMATION
      CALL TRICOF(Four_sp,NF,FA,NE,FB,NO,1)
      A0 = FA(1) / 2.
      TQ = A0

! coefficient to transform solution; Koeffizienten fuer Loesung transformieren 
      DO I=1,NK
        FTA(I) = SQRT(FA(I+1)*FA(I+1) + FB(I)*FB(I))
        FTA(I) = FTA(I) * SIGN(1.,FB(I))
        if(FB(I).eq. 0.) then
		    FTO(I) = pi/2.
		else
	        FTO(I) = ATAN(FA(I+1)/FB(I))
		end if
        FTO(I) = FTO(I) - (NEND+NAF)*PI*I/(NEND-NAF)
      ENDDO

END SUBROUTINE BTFOUR
 
!******************************************************************************