      program thetafao
C---------------------------------------------------------------
c      version 0.7, Last Modified 11/14/2022
C      WRITTEN FOR: EU FOCUS PRZM/VFSMOD tool 
C      Written by: R. Munoz-Carpena (rmc)
C                  University of Florida, carpena@ufl.edu
C---------------------------------------------------------------
C Program to calculate soil moisture content between runoff events.
C This is a necessary step for continuous simulation of the PRZM/VFSMOD
C EU FOCUS tool. The result of the calculation will produce the initial
C moisture content of the soil (OI) for the next VFSMOD run in the time
C series. It follows FAO-56 adjusted ET calculations (Allen et al.,1998)
C based on Dr. M. Qußemada (U. Politecnica Madrid) spreadsheet calculations
C and good results in the comparison with field measured soil moisture.
c---------------------------------------------------------------
C   Input parameters
c   isoiltype (USDA, S:Sand;L:Loam;s:Silt;C:Clay):-1:user,1:S,2:LS,3:SL,4:L,5:sL,6:s,7:sCL,8:sC,9:C
c   OI(m3/m3): top soil initial water content (same as in VFSMOD *.iso file)
c   FC(m3/m3): top soil fied capacity water content (read internally or provided by user when isoiltype=-1)
c   WP(m3/m3): top soil wilting point water content (read internally or provided by user when isoiltype=-1)
c   Zr(m): maximum grass root zone depth (typical values (0.5-1.5 m)
c   pfrac[-]: fraction of easily estractable water (typical 0.6 for Bermuda grass)
c   Hm(m): height of vegetation (from VFSMOD *.igr file, H(cm)/100)
c   iFH(optional): input MET file formating flag where,
C     iFM= 0 (or not present), 8 columns,last two columns are Tmin,Tmax
C     iFM= 1, 7 columns, last column is HRmean
C     iFM= 2, 9 columns, last 3 columns are Tmin,Tmax, Kcmid for a crop other than grass
C     iFM= 3, 8 columns, last 2 columns are HRmin, Kcmid for a crop other than grass 
c   REW(mm): readily extractable water (soil dependent)
c   soil(isoiltype,1): FC, top soil field capacity (m3/m3)
c   soil(isoiltype,2): WP, top soil wilting point (m3/m3)
c   soil(isoiltype,1): top soil REW(mm) (see above)
c   TAW(mm): total available water
C---------------------------------------------------------------
c    Compiling for Win32 and Unix environments:
c       1. The i/o for these operating systems is different.
c       2. Change the comments in the finput.f program to reflect
c          your operating system.   3/9/2012
c---------------------------------------------------------------
c    CHANGES
c    v0.7, 11/14/2022. Added AET (mm) to last column of modified .MET output file and fix field 
c      length for date that was short of 1 character on the output .MET file
c    v0.6, 12/19/2019. Changed iRH into iFM (input formatting MET)
c      to modify the last columns of the MET file with user cKmid
c---------------------------------------------------------------

      implicit double precision (a-h, o-z)
      CHARACTER*75 LISFIL(6)

c------------------------------------------
c  open files and get inputs 
c------------------------------------------
      call finput(lisfil)
      call getinp(isoiltype,OI,FC,WP,Zr,pfrac,Hm,TAW,cKmid,iFM)
c--------------------------
c Calculate runoff volume by SCS method
c--------------------------
      call theta10(OI,Zr,pfrac,Hm,cKmid,WP,FC,iFM)
c--------------------------
c Output results
c--------------------------
      call  results(isoiltype,OI,FC,WP,Zr,pfrac,Hm,TAW,cKmid,iFM) 
c----------
c Output message at end of program -----------------
c----------
      WRITE(*,*)
      WRITE(*,*)'...FINISHED!...'
      WRITE(*,*)


      close(1)
      close(2)
      close(3)
      close(4)

      stop
      end


      SUBROUTINE FINPUT(LISFIL)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     Create input and output file names from a command line input string    C
C     NOTE: Maximum length of command line string = 50                       C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      CHARACTER*50 FILENM1
      CHARACTER*75 LISFIL(6)
      CHARACTER*3,SCOD(6)
      CHARACTER*1,DUMMY1
      character*200 linein
      character*1 slash
      DATA(SCOD(I),I=1,3)/'in','met','out'/

c*** Command line option to input filename
c*** Comment out the following depending for which system you compile

CWIN32*** Start of Win32 file i/o     ***
CWIN32      slash='\'
CWIN32      INARGS=NARGS()-1
CWIN32      IF (INARGS.EQ.1) THEN
CWIN32      CALL GETARG(1,FILENM1,IFSTATUS)
CWIN32*** End of Win32 file i/o       ***
cUNIX *** Start Unix file i/o         ***
CUNIX
      slash='/'
CUNIX
      INARGS=IARGC()
CUNIX
      IF (INARGS.EQ.1) THEN
CUNIX
         CALL GETARG(1,FILENM1)
cUNIX *** End of UNIX file i/o section ***
       ELSE
         WRITE(*,*)
         WRITE(*,105)
         WRITE(*,110)
         WRITE(*,130)
         WRITE(*,140)
         WRITE(*,150)
         WRITE(*,*)
         STOP
      ENDIF
c----------
c Output message at end of program -----------------
c----------
      WRITE(*,*)
      WRITE(*,*)'THETAFAO v0.7, 11/2022'
      WRITE(*,*)'R.Munoz-Carpena (UFL), carpena@ufl.edu'
      WRITE(*,*)

C------ create I/O filenames from input string -------------------------
c------    or read filenames from a project file -----------------------
c
      ilstr=index(filenm1,'.')
      if (ilstr.gt.0) then
c      *** using project file (.prj or .lis) to read filenames
c      ***  check to see if extension is .prj or .lis
          ilstr1=index(filenm1,'.prj')
          ilstr2=index(filenm1,'.lis')
          if ((ilstr1.gt.0).or.(ilstr2.gt.0)) then
c          *** fill filename array with safe names
             do 11 i=1,3
               dummy1=scod(i)
               IF(DUMMY1.EQ.'i') THEN
                   WRITE(LISFIL(I),'(5A)')
     &             'inputs',slash,'dummy','.',SCOD(I)
                  ELSE
                   WRITE(LISFIL(I),'(5A)')
     &             'output',slash,'dummy','.',SCOD(I)
               ENDIF
11           continue
c
             open(unit=99,file=filenm1,status='old')
12           read(99,'(a)',end=18) linein
             lpos=index(linein,'=')
             lstr=len(linein)
             if ((lpos.gt.0).and.(lstr.gt.0)) then
                do 14 jj=1,3
                   lpp = index(linein(1:lpos-1),scod(jj))
                   if (lpp.gt.0) lisfil(jj)=linein(lpos+1:)
14              continue
             endif
             go to 12
c           ****** done
18           continue
            else
             WRITE(*,*)
             WRITE(*,105)
             WRITE(*,110)
             WRITE(*,130)
             WRITE(*,140)
             WRITE(*,150)
             WRITE(*,*)
             STOP
          endif
        else
c     **** rafa's i/o scheme
          ILSTR=INDEX(FILENM1,' ')-1
          DO 101 I=1,3
            DUMMY1=SCOD(I)
            IF(DUMMY1.EQ.'i'.or.DUMMY1.EQ.'m') THEN
                  WRITE(LISFIL(I),'(5A)')
     &        'inputs',slash,FILENM1(:ILSTR),'.',SCOD(I)
               ELSE
                  WRITE(LISFIL(I),'(5A)')
     &        'output',slash,FILENM1(:ILSTR),'.',SCOD(I)
            ENDIF
c--- recreate .MET file in output directory
            WRITE(LISFIL(4),'(5A)')
     &        'output',slash,FILENM1(:ILSTR),'.',SCOD(2)
101      CONTINUE
      endif

      write(*,*)'...Opening files...'
      DO 102 I=1,4
         write(*,'(70A)')lisfil(i)
102   CONTINUE

C---------Open I/O files -------------------------------
      OPEN(1,FILE=LISFIL(1),STATUS='OLD')
      OPEN(2,FILE=LISFIL(2),STATUS='OLD')
      OPEN(3,FILE=LISFIL(3),STATUS='unknown')
      OPEN(4,FILE=LISFIL(4),STATUS='unknown')
      WRITE(3,220)LISFIL(3)

105   FORMAT('Name:    thetafao')
110   FORMAT(9x,'(VFSMOD-FOCUS:calculate soil moisture between events)')
130   FORMAT('Usage:   thetafao filename (max 8 characters, no ext.)')
CWIN32
140   FORMAT('Version: 0.7 for Windows -Nov 2022')
CUNIX140   FORMAT('Version: 0.7 for Unix -Nov 2022')
150   FORMAT('Author: R.Munoz-Carpena (UFL), carpena@ufl.edu')
160   FORMAT(72('-'))
220   FORMAT('File: ',A40,9x,'THETAFAO v0.7, 11/2022')

      RETURN
      END


      subroutine getinp(isoiltype,OI,FC,WP,Zr,pfrac,Hm,
     & TAW,cKmid,iFM)
C---------------------------------------------------------------
C   Input parameters
c   isoiltype (USDA, S:Sand;L:Loam;s:Silt;C:Clay):-1:user,1:S,2:LS,3:SL,4:L,5:sL,6:s,7:sCL,8:sC,9:C
c   OI(m3/m3): top soil initial water content (same as in VFSMOD *.iso file)
c   FC,soil(isoiltype,1): top soil field capacity (m3/m3)
c   WP,soil(isoiltype,2): top soil wilting point (m3/m3)
c   REW(mm): readily extractable water (soil dependent)
c   Zr(m): maximum grass root zone depth (typical values (0.5-1.5 m)
c   pfrac[-]: fraction of easily estractable water (typical 0.6 for Bermuda grass)
c   Hm(m): height of vegetation (from VFSMOD *.igr file, H(cm)/100)
c   iFH(optional): input MET file formating flag where,
C     iFM= 0 (or not present), 8 columns,last two columns are Tmin,Tmax
C     iFM= 1, 7 columns, last column is HRmean
C     iFM= 2, 9 columns, last 3 columns are Tmin,Tmax, Kcmid for a crop other than grass
C     iFM= 3, 8 columns, last 2 columns are HRmin, Kcmid for a crop other than grass 
c   REW,soil(isoiltype,1): top soil REW(mm) (see above)
c   TAW(mm): total available water
C---------------------------------------------------------------
      implicit double precision (a-h, o-z)
      dimension soil(9,3)

      data(soil(1,J),J=1,3)/0.12d0,0.045d0,4.5d0/
      data(soil(2,J),J=1,3)/0.15d0,0.065d0,6.d0/
      data(soil(3,J),J=1,3)/0.23d0,0.11d0,8.d0/
      data(soil(4,J),J=1,3)/0.25d0,0.12d0,9.d0/
      data(soil(5,J),J=1,3)/0.29d0,0.15d0,9.5d0/
      data(soil(6,J),J=1,3)/0.32d0,0.17d0,9.5d0/
      data(soil(7,J),J=1,3)/0.335d0,0.205d0,9.5d0/
      data(soil(8,J),J=1,3)/0.36d0,0.23d0,10.d0/
      data(soil(9,J),J=1,3)/0.36d0,0.22d0,10.d0/

      cKmid=1.d0
      iFM=0
c   iFM =1 or 2: cKmid is read in met file and Kc
c   calculation is bypassed
      read(1,*)isoiltype
      if(isoiltype.eq.-1) then
            backspace(1)
            read(1,*,end=10)isoiltype,OI,FC,WP
         elseif(isoiltype.le.9) then
            FC=soil(isoiltype,1)
            WP=soil(isoiltype,2)
            OI=FC
         else
         	write(*,*)'ERROR: wrong soil type selection (-1,9)'
         	STOP
      endif
      read(1,*,end=10)Zr,pfrac,Hm
      read(1,*,end=4)iFM
4     return
10    print*,'Error reading input file, please check inputs'
      end


      subroutine theta10(OI,Zr,pfrac,Hm,cKmid,WP,FC,iFM)
C---------------------------------------------------------------
C Soil moisture calculation, FAO (1998)
C MET file: the last columns can be changed to be read from the 
C   modified EUFOCUS MET files (last 2 columns Tmax,Tmin), and/or 
C    provide crop coefficient values (kcmid) for plants other than 
C    than grass using the different values of iFM: 
C  iFM= 0 (or not present), 8 columns,last two columns are Tmin,Tmax
C  iFM= 1, 7 columns, last column is HRmean
C  iFM= 2, 9 columns, last 3 columns are Tmin,Tmax, Kcmid for a crop other than grass
C  iFM= 3, 8 columns, last 2 columns are HRmin, Kcmid for a crop other than grass 
C---------------------------------------------------------------
      implicit double precision (a-h, o-z)
      character*7 datein
      character*100 dum

c-----Initial calculations before the data processing loop
c-----root zone and top layer initial depletion
      Dr=(FC-OI)*Zr*1000.d0
      TAW=(FC-WP)*Zr*1000.d0
      
c-----Start calculation loop for original .MET file and save file with additional column in
c-----output directory. The scratch file is used to read PRZM date fix format (first column 'A7')
c-----while the other inputs are are read as free format. 
      do 10 i=1,10000
        open(25,status='scratch')
        read(2,101,end=30)datein,dum
        write(25,*)dum
        rewind(25)
        SELECT CASE (iFM)
           CASE (1)
              read(25,*)prec,ETo,temp,u2,rad,RHm
           CASE (2)
              read(25,*)prec,ETo,temp,u2,rad,tempmax,tempmin,cKmid
              RHm=RHmin(tempmax,tempmin)
           CASE (3)
              read(25,*)prec,ETo,temp,u2,rad,RHm,cKmid
           CASE DEFAULT
              read(25,*)prec,ETo,temp,u2,rad,tempmax,tempmin
cdebug        write(*,'(i6,8f8.2)')i,prec,ETo,temp,u2,rad,tempmax,tempmin
              RHm=RHmin(tempmax,tempmin)
        END SELECT
        close(25)
c------change units to mm/day (P,ETo) and m/s (u2) from default in MET file ------
        precmm=prec*10.d0
        ETomm=ETo*10.d0
        u2m=u2/100.d0
		if(u2m.lt.0.d0)u2m=0.d0
c-----Calculate adjusted crop coefficient (cKb) based on daily wind and relative temperature
        cKbadj=(0.04d0*(u2m-2.d0)-0.004d0*(RHm-45.d0))*(Hm/3.d0)**0.3d0
        cKb=cKmid+cKbadj
        ETc=cKb*ETomm
c-----Calculate root zone layer depletion (Dr) from the daily soil water balance	  	
        DP= precmm-ETc-Dr
        if(DP.lt.0.d0) DP=0.d0
        Dr=Dr-precmm+ETc+DP
        if(Dr.lt.0.d0) Dr=0.d0
        if(Dr.gt.TAW)  Dr=TAW
c-----Calculate soil water stress coefficient (cKs)	
        peff=pfrac+0.04d0*(5.d0-ETc)
        RAW=TAW*peff
        if(Dr.gt.RAW)then
            cKs=(TAW-Dr)/(TAW-RAW)
          else
            cKs=1.d0
        endif
c-----Calculate the actual ET (ETa) and calculate the topsoil moisture
        ETa=cKb*cKs*ETomm
c-----Calculate the soil moisture (depth Zr)
        thetamm=1000.d0*FC*Zr-Dr
        theta=thetamm/(Zr*1000.d0)
c-----Write the new theta10 and ETa columns to the output/*.met file	  	
cdebug        write(*,200)i,precmm,ETomm,u2m,RHm,cKb,ETc,pfrac,peff,
cdebug     &     RAW,TAW,Dr,DP,cKs,ETa,theta
        if(iFM.eq.1.or.iFM.eq.3) then
c           write(4,102)datein,prec,ETo,temp,u2,rad,RHm,theta
           write(4,102)datein,prec,ETo,temp,u2,rad,RHm,theta,ETa/10.d0
         else	
c           write(4,100)datein,prec,ETo,temp,u2,rad,tempmax,tempmin,
c     &       theta
           write(4,100)datein,prec,ETo,temp,u2,rad,tempmax,tempmin,
     &       theta,ETa/10.d0
        endif
10    continue

100   format(a7,2f10.2,f10.1,f10.0,f10.1,2f10.2,f10.3,f10.3)
101   format(a7,a100)
102   format(a7,2f10.2,f10.1,f10.0,f10.1,f10.2,f10.3,f10.3)
cdebug 200   format(i6,20f10.3)

30    return
      end
 
      
      function RHmin(tempmax,tempmin)
c---------------------------------------------------------------
c--- Estimate missing RHmin based on assumption Tmin<>Tdew
c---------------------------------------------------------------
      implicit double precision (a-h, o-z)
      
      a1= 17.625d0
      b1= 243.04d0
      c1= 0.61121d0
      ea= c1*dexp(a1*tempmin/(tempmin+b1))
      es= c1*dexp(a1*tempmax/(tempmax+b1))
      RHmin=100.d00*ea/es
      if(RHmin.gt.100.d0) RHmin=100.d0

      return
      end


      subroutine results(isoiltype,OI,FC,WP,Zr,pfrac,Hm,TAW,cKmid,iFM) 
C---------------------------------------------------------------
C      Output summary of hydrology results nicely
C---------------------------------------------------------------
      implicit double precision (a-h, o-z)
      character*13 stype(9)
      data stype/'Sand','LoamySand','SandyLoam','Loam','SiltyLoam',
     1 'Silt','SiltyClayLoam','SiltyClay','Clay'/
     
      TAW=(FC-WP)*Zr*1000.d0
      write(3,*)' '
      write(3,*)'TOP SOIL MOISTURE CALCULATION FAO (1998) METHOD'
      write(3,*)' '
      write(3,*)'INPUTS'
      write(3,*)'------'
      if(isoiltype.eq.-1) then
      		write(3,200)'User'
      	else
      		write(3,200)stype(isoiltype)
      endif
      write(3,250)FC
      write(3,300)WP
      write(3,325)OI
      write(3,450)Zr
      write(3,500)pfrac      
      write(3,600)TAW
	  if (iFM.eq.2.or.iFM.eq.3) then
		write(3,655)
	   else
		write(3,650)cKmid
      endif
      write(3,700)Hm
      SELECT CASE (iFM)
         CASE (1)
             write(3,800)iFM,'last column RHmean            '
         CASE (2)
             write(3,800)iFM,'last 3 columns Tmax,Tmin,cKmid'
         CASE (3)
             write(3,800)iFM,'last 2 columns RHmean,cKmid   '
         CASE DEFAULT
             write(3,800)iFM,'last 2 columns Tmax,Tmin      '
      END SELECT
      write(3,*)' '

200   format('Soil type',34x,'=',a13)
250   format('Top soil field capacity, FC(m3/m3)',14x,'=',f9.3)
300   format('Top soil wilting point, WP(m3/m3)',15x,'=',f9.3)
325   format('Top soil initial water content, OI(m3/m3)',7x,'=',f9.3)
450   format('Maximum grass root zone depth, Zr(m)',12x,'=',f8.2)
500   format('Fraction of easily extractable water,pfrac',6x,'=',f8.2)
600   format('Total available water, TAW(mm)',18x,'=',f8.2)
650   format('Mid season crop coeff., Kcmid',19x,'='f8.2)
655   format('Mid season crop coeff., Kcmid',19x,
     & '=    variable (.MET file)')
700   format('Vegetation height, H(m)',25x,'=',f8.2)
800   format('Input format option for MET file, iFM =',i2,7x,'=',a34)

      return
      end

