! Test and compare PET between Thornthwaite and Hargreaves

Use NetCDF

Integer, parameter :: Nmon=12
Real, Allocatable :: Tx(:,:), Tm(:,:), T2m(:,:), R(:,:), Lat(:)
Real, Allocatable :: PET_Th(:,:),PET_Ha(:,:),PET(:), Tx0(:), Tm0(:),T20(:)
Integer :: Nst, Yr1, Yr2, Ntime, imon
Real :: xlat
Character (Len=255) :: Fname
Logical :: Status

Fname="/home/tanpv/Data/OBS/T_R_Monthly_for_Dr/R_Tx_Tm_T2m_1980-2020_mon_3D_VN.nc"
Call Read_Data (Fname)

Allocate (PET_Th(Ntime,Nst),PET_Ha(Ntime,Nst))
Allocate (Tx0(0:Ntime-1),Tm0(0:Ntime-1),T20(0:Ntime-1),PET(0:Ntime-1))

Do ist=1,Nst
   xlat=Lat(ist)
   T20(0:Ntime-1)=T2m(0:Ntime-1,ist)
   Tx0(0:Ntime-1)=Tx (0:Ntime-1,ist)
   Tm0(0:Ntime-1)=Tm (0:Ntime-1,ist)
   Call Thornthwaite (T20,Ntime,xlat,PET)
   PET_Th(1:Ntime,ist)=PET(0:Ntime-1)
   Call Hargreaves (Tx0,Tm0,Ntime,xlat,PET)
   PET_Ha(1:Ntime,ist)=PET(0:Ntime-1)
Enddo
Fname="PET_Th_Ha.nc"
Call Write_PET (Fname)

CONTAINS
!----------------------------------------------------------------------------------------
Subroutine Write_PET (Fname)
  Character (Len=*) :: Fname
  Real :: Miss,Stn(Nst),xtime(Ntime)
  Integer :: ist, it

  Miss=-99.0
  Do ist=1,Nst; Stn(ist)=ist; Enddo
  Do it=1,Ntime; xtime(it)=it; Enddo

  Call wr_ncdf_2d (Fname,xtime,Stn,Ntime,Nst,PET_Th,"th","PET-Thornthwaite","mm/mon","month","stn",Miss)
  Call wr_ncdf_2d_apd (Fname,PET_Ha,"ha",Ntime,Nst,"PET-Hargreaves","mm/mon","month","stn",Miss)
End Subroutine Write_PET
!----------------------------------------------------------------------------------------
Subroutine Read_Data (Fname)
  Character (Len=*) :: Fname
  Integer :: Nyrs,Nmn,iyr,imon,ist,it
  Real, Allocatable :: X1(:,:,:),X2(:,:,:),X3(:,:,:),X4(:,:,:),xyear(:),xmon(:),xStn(:)

  Call rd_ncdf_3d (Fname,xStn,xyear,xmon,Nst,Nyrs,Nmn,X1,"Tx")
  Deallocate (xStn,xyear,xmon)
  Call rd_ncdf_3d (Fname,xStn,xyear,xmon,Nst,Nyrs,Nmn,X2,"Tm")
  Deallocate (xStn,xyear,xmon)
  Call rd_ncdf_3d (Fname,xStn,xyear,xmon,Nst,Nyrs,Nmn,X3,"T2m")
  Deallocate (xStn,xyear,xmon)
  Call rd_ncdf_3d (Fname,xStn,xyear,xmon,Nst,Nyrs,Nmn,X4,"R")
  Yr1=Int(xyear(1)); Yr2=Int(xyear(Nyrs))
  Ntime=Nyrs*Nmn
  Deallocate (xStn,xyear,xmon)
  Call rd_ncdf_1d (Fname,xStn,Nst,Lat,"Lat")
  Allocate (Tx(0:Ntime-1,Nst),Tm(0:Ntime-1,Nst),T2m(0:Ntime-1,Nst),R(0:Ntime-1,Nst))
  Do ist=1,Nst
     it=-1
     Do iyr=1,Nyrs
        Do imon=1,Nmn
           it=it+1
           Tx(it,ist)=X1(ist,iyr,imon)
           Tm(it,ist)=X2(ist,iyr,imon)
           T2m(it,ist)=X3(ist,iyr,imon)
           R(it,ist)=X4(ist,iyr,imon)
        Enddo
     Enddo
  Enddo
End Subroutine Read_Data
!----------------------------------------------------------------------------------------
! thornthwaite.F90:
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

!-------------------------------------------------------------------
Subroutine Hargreaves (Tx, Tm, n, lat, etpSeries)
  ! ETo = 0.0023 * RA * (Tave + 17.8) * TD**0.5
  ! ETo = Krs * RA * (Tave + HT) * TD**HE
  Integer, Parameter :: Nmon=12, Nlat=51
  Real, Parameter :: Krs=0.0023, HT=17.8, HE=0.5
  real, parameter :: days (Nmon) = (/31,28,31,30,31,30,31,31,30,31,30,31/)
  integer :: n
  real :: Tx(0:n-1), Tm(0:n-1)
  real :: Tave (0:n-1), TD(0:n-1), lat, etpSeries (0:n-1)
  Real :: RA0(Nlat,Nmon), Lat0(Nlat), RA(Nmon)
  integer :: month, i, k

  DATA Lat0 &
  /-50.0,-48.0,-46.0,-44.0,-42.0,-40.0,-38.0,-36.0,-34.0,-32.0,-30.0,-28.0,-26.0,-24.0,-22.0,-20.0,-18.0,-16.0,-14.0,-12.0,-10.0, -8.0, -6.0, -4.0, -2.0,  0.0, &
     2.0,  4.0,  6.0,  8.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0, 22.0, 24.0, 26.0, 28.0, 30.0, 32.0, 34.0, 36.0, 38.0, 40.0, 42.0, 44.0, 46.0, 48.0, 50.0/
  DATA RA0 &
  /17.5, 17.6, 17.7, 17.8, 17.8, 17.9, 17.9, 17.9, 17.8, 17.8, 17.8, 17.7, 17.6, 17.5, 17.4, 17.3, 17.1, 16.9, 16.7, 16.6, 16.4, 16.1, 15.8, 15.5, 15.3, 15.0, &
   14.7, 14.3, 13.9, 13.6, 13.2, 12.8, 12.4, 12.0, 11.6, 11.2, 10.7, 10.2,  9.8,  9.3,  8.8,  8.3,  7.9,  7.4,  6.9,  6.4,  5.9,  5.3,  4.9,  4.3,  3.8, 14.7, &
   14.9, 15.1, 15.3, 15.5, 15.7, 15.8, 16.0, 16.1, 16.2, 16.4, 16.4, 16.4, 16.5, 16.5, 16.5, 16.5, 16.4, 16.4, 16.3, 16.3, 16.1, 16.0, 15.8, 15.7, 15.5, 15.3, &
   15.0, 14.8, 14.5, 14.2, 13.9, 13.6, 13.3, 13.0, 12.7, 12.3, 11.9, 11.5, 11.1, 10.7, 10.2,  9.8,  9.4,  9.0,  8.6,  8.1,  7.6,  7.1,  6.6,  6.1, 10.9, 11.2, &
   11.5, 11.9, 12.2, 12.5, 12.8, 13.2, 13.5, 13.8, 14.0, 14.3, 14.4, 14.6, 14.8, 15.0, 15.1, 15.2, 15.3, 15.4, 15.5, 15.5, 15.6, 15.6, 15.7, 15.7, 15.6, 15.5, &
   15.4, 15.3, 15.3, 15.1, 14.9, 14.7, 14.6, 14.4, 14.2, 13.9, 13.7, 13.4, 13.1, 12.8, 12.4, 12.1, 11.8, 11.4, 11.0, 10.6, 10.2,  9.8,  9.4,  7.0,  7.5,  7.9, &
    8.4,  8.8,  9.2,  9.6, 10.1, 10.5, 10.9, 11.3, 11.6, 12.0, 12.3, 12.6, 13.0, 13.2, 13.5, 13.7, 14.0, 14.2, 14.4, 14.7, 14.9, 15.1, 15.3, 15.3, 15.5, 15.4, &
   15.6, 15.7, 15.7, 15.7, 15.6, 15.6, 15.6, 15.5, 15.4, 15.3, 15.3, 15.2, 15.0, 14.8, 14.7, 14.5, 14.3, 14.0, 13.7, 13.3, 13.0, 12.7,  4.2,  4.7,  5.2,  5.7, &
    6.1,  6.6,  7.1,  7.5,  8.0,  8.5,  8.9,  9.3,  9.7, 10.2, 10.6, 11.0, 11.4, 11.7, 12.1, 12.5, 12.8, 13.1, 13.4, 13.8, 14.1, 14.4, 14.6, 14.9, 15.1, 15.3, &
   15.5, 15.7, 15.8, 16.0, 16.1, 16.3, 16.3, 16.4, 16.4, 16.5, 16.5, 16.5, 16.5, 16.4, 16.4, 16.4, 16.2, 16.1, 16.0, 15.9, 15.8,  3.1,  3.5,  4.0,  4.4,  4.9, &
    5.3,  5.8,  6.3,  6.8,  7.3,  7.8,  8.2,  8.7,  9.1,  9.6, 10.0, 10.4, 10.8, 11.2, 11.6, 12.0, 12.4, 12.8, 13.2, 13.5, 13.9, 14.2, 14.4, 14.7, 15.0, 15.3, &
   15.5, 15.7, 15.9, 16.1, 16.4, 16.4, 16.6, 16.7, 16.8, 17.0, 17.0, 17.1, 17.2, 17.2, 17.3, 17.3, 17.2, 17.2, 17.2, 17.1,  3.5,  4.0,  4.4,  4.9,  5.4,  5.9, &
    6.3,  6.8,  7.2,  7.7,  8.1,  8.6,  9.1,  9.5, 10.0, 10.4, 10.8, 11.2, 11.6, 12.0, 12.4, 12.7, 13.1, 13.4, 13.7, 14.1, 14.3, 14.6, 14.9, 15.1, 15.3, 15.5, &
   15.7, 15.9, 16.1, 16.3, 16.4, 16.5, 16.6, 16.7, 16.8, 16.8, 16.8, 16.7, 16.7, 16.7, 16.7, 16.6, 16.6, 16.5, 16.4,  5.5,  6.0,  6.5,  6.9,  7.4,  7.9,  8.3, &
    8.8,  9.2,  9.6, 10.1, 10.4, 10.9, 11.2, 11.6, 12.0, 12.3, 12.6, 12.9, 13.2, 13.5, 13.7, 14.0, 14.3, 14.5, 14.8, 14.9, 15.1, 15.2, 15.4, 15.5, 15.6, 15.7, &
   15.7, 15.8, 15.9, 15.8, 15.8, 15.7, 15.7, 15.7, 15.6, 15.5, 15.4, 15.3, 15.2, 15.0, 14.7, 14.5, 14.3, 14.1,  8.9,  9.3,  9.7, 10.2, 10.6, 11.0, 11.4, 11.7, &
   12.0, 12.4, 12.7, 13.0, 13.2, 13.4, 13.7, 13.9, 14.1, 14.3, 14.5, 14.7, 14.8, 14.9, 15.0, 15.1, 15.2, 15.3, 15.3, 15.3, 15.3, 15.3, 15.3, 15.2, 15.1, 15.0, &
   14.9, 14.8, 14.6, 14.5, 14.3, 14.1, 13.9, 13.6, 13.4, 13.1, 12.8, 12.5, 12.2, 11.9, 11.5, 11.2, 10.9, 12.9, 13.2, 13.4, 13.7, 14.0, 14.2, 14.4, 14.6, 14.9, &
   15.1, 15.3, 15.4, 15.5, 15.6, 15.7, 15.8, 15.8, 15.8, 15.8, 15.8, 15.9, 15.8, 15.7, 15.6, 15.5, 15.4, 15.3, 15.1, 15.0, 14.8, 14.7, 14.4, 14.1, 13.9, 13.6, &
   13.3, 13.0, 12.6, 12.3, 12.0, 11.6, 11.2, 10.8, 10.6, 10.0,  9.6,  9.1,  8.7,  8.3,  7.8,  7.4, 16.5, 16.6, 16.7, 16.7, 16.8, 16.9, 17.0, 17.0, 17.1, 17.2, &
   17.3, 17.2, 17.2, 17.1, 17.0, 17.0, 16.8, 16.7, 16.5, 16.4, 16.2, 16.0, 15.8, 15.5, 15.3, 15.1, 14.8, 14.5, 14.2, 13.9, 13.6, 13.3, 12.8, 12.4, 12.0, 11.6, &
   11.1, 10.7, 10.3,  9.9,  9.5,  9.0,  8.5,  8.0,  7.5,  7.0,  6.5,  6.0,  5.5,  5.0,  4.5, 18.2, 18.2, 18.3, 18.3, 18.3, 18.3, 18.3, 18.2, 18.2, 18.1, 18.1, &
   17.9, 17.8, 17.7, 17.5, 17.4, 17.1, 16.8, 16.6, 16.5, 16.2, 16.0, 15.7, 15.4, 15.1, 14.8, 14.4, 14.1, 13.7, 13.3, 12.9, 12.5, 12.0, 11.6, 11.1, 10.7, 10.2, &
    9.7,  9.3,  8.8,  8.3,  7.8,  7.2,  6.6,  6.1,  5.7,  5.2,  4.7,  4.3,  3.7,  3.2 /
   !-------------------------------------------------
   ! Calculate RA = f(lat,mon)
   Do k=1,Nmon
      Do i=1,Nlat-1
         if (lat>Lat0(i).AND.lat<=Lat0(i+1)) exit
      Enddo
      RA(k)=RA0(i,k)+(RA0(i+1,k)-RA0(i,k))/(Lat0(i+1)-Lat0(i))*(lat-Lat0(i))
   Enddo
   ! Calculate PET
   k=0
   Do i=0,n-1
      if (Tx(i)>0.AND.Tm(i)>0) then
         Tave(i)=(Tx(i)+Tm(i))/2.0; TD(i)=Tx(i)-Tm(i)
      endif
      k=mod(i+1,Nmon); if (k==0) k=Nmon
      etpSeries(i) = Krs * RA(k) * (Tave(i) + HT) * TD(i)**HE  ! mm/day
      etpSeries(i) = etpSeries(i) * days(k)                    ! mm/month
   Enddo
End Subroutine
!-------------------------------------------------------------------

  Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90"
  Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90"
END