Subroutine thornthwaite (tempSeries, n, lat, etpSeries)
  !! Calculates the potential evapotranspiration from a series of monthly
  !! temperature data, following the method of Thornthwaite.
  integer :: n
  real :: tempSeries (0:n-1), lat, etpSeries (0:n-1)

  real :: T, J, J2, J3, c, N_up, omega, K_up(13), tanLatMonth
  integer :: month, i, k

  real, parameter ::  pi = 3.14159265358979
  !! Tangens of the average solar declination angle for each month of the year
  real, parameter ::  tanDelta (0:11) = (/-0.37012566, -0.23853358, -0.04679872,  0.16321764, &
                                           0.32930908,  0.40677729,  0.37477410,  0.23906300, &
                                           0.04044485, -0.16905776, -0.33306377, -0.40743608/)
  real, parameter ::  days (0:11) = (/31,28,31,30,31,30,31,31,30,31,30,31/)

  real :: tanLat
  tanLat = tan(0.0174532925*lat)
  
  print*,"tanLat:", tanLat  
  !! -----------------------------------------------------------------------------------------
  !! Compute J: annual temperature efficiency index
  J = 0.
  do month = 1, 12
    T = 0. ; k = 0
    do i = month - 1, n - 1, 12
      T = T + tempSeries (i)
      k = k + 1
    enddo
    !! T: monthly average temperature
    T = T / real(k)
    !! print*,"CHECK:", J !! NXTTNX
    if (T > 0.) J = J + (T / 5.)**1.514
  enddo

  !! Compute c exponent
  J2 = J*J
  J3 = J2*J
  c = 0.000000675*J3 - 0.0000771*J2 + 0.01792*J + 0.49239

  !! print*,"CHECK:", c !! NXTTNX
  !! Compute K: monthly correction factor, depending on latitude
  do month = 1, 12
    tanLatMonth = tanLat*tanDelta(month-1)
    if ( (tanLatMonth .gt. -1.) .and. (tanLatMonth .lt. 1.) ) omega = acos(-tanLatMonth)
    if ( tanLatMonth .le. -1. ) omega = 0.
    if ( tanLatMonth .ge. 1. )  omega = pi
    N_up = 24.*omega / pi
    K_up(month) = N_up/12. * days(month-1)/30.
  enddo

  !! Compute potential evapotranspiration series
  month = 1
  do i = 0, n - 1
    if ( (tempSeries(i) .gt. 0.) .and. (J .gt. 0.) ) then
      etpSeries(i) = K_up(month) * 16. * (10*tempSeries(i)/J)**c
    else
      etpSeries(i) = 0.
    endif
    
    if (month .lt. 12) then
      month = month + 1
    else 
      month = 1
    endif
  enddo
                      
End Subroutine thornthwaite