;+
; NAME: noahSoilsWTexture
;
; PURPOSE: Runs noah for all soils of a given texture specified in a pair of file.
;          Looks through the database of SHPs generated by MOSCEM (runMoscemSHP.pro)
;          and compHT2HK.pro, and uses the general.txt file from UNSODA to pick out texture
;
; CATEGORY: Batch running noah
;
;
; CALLING SEQUENCE: noahSoilsWTexture, texture, datafilename, texturefilename
;
; INPUTS: texture = a string identifying the texture you wish to run noah for
;                   e.g. "sandy loam" "clay" "silt clay loam" see the routine textIndex
;                   for a complete list of textures
;         datafilename = filename of compHT2HK.pro outputfile
;         texturefilename = filename of "general.txt" file from UNSODA database
;
;
; OPTIONAL INPUTS: <none>
;
; KEYWORD PARAMETERS: <none>
;
; OUTPUTS: noah output files named "out_texture_soilNumber"
;            where : texture = the inputtexture name (without spaces)
;                    soilNumber= the soil index number in the UNSODA database
;
; OPTIONAL OUTPUTS: <none>
;
; COMMON BLOCKS: <none>
;
; SIDE EFFECTS: file SOILPARM.TBL in current directory is overwritten
;
; RESTRICTIONS: Current directory must contain Noah control files
;               IHOPstyp must be set to 1
;
; PROCEDURE: Get an array of [soilNumber,Texture] from the texture file.
;            Read through the datafile, for each soil
;              If it matches the texture requested, and SHPs are reasonable
;                Write SOILPARM.TBL file
;                Run noah
;                Rename output file
;
; EXAMPLE:
;
; MODIFICATION HISTORY:
;           8/2/2004 edg - Original
;-


FUNCTION getTextures, textureFile
  openr, un, /get, textureFile
  line=''
  output=strarr(2)
  i=0
  WHILE NOT eof(un) DO BEGIN
     readf, un, line
     data=strsplit(line, '~', /extract, /preserve_null)
     texture=(strsplit(data[3], '"', /extract))[0]
     output=[[output],[[data[0], texture] ]]
     i++
  ENDWHILE
  close, un
  free_lun, un
  output=output[*,1:--i]
  return, output
END

;; Writes a one line SOILPARM.TBL file for input to the Noah model.
;; IHOPstyp must be set to soil number 1
PRO writeSoilFile, data, index, vg=vg, cam=cam
  data=float(data)
  b=data[5]
  DRYSMC=data[9]
  MAXSMC=data[8]
  IF keyword_set(cam) THEN $
    SATPSI=data[6]/100.         ; convert cm to m
  IF keyword_set(vg) THEN $
    SATPSI=data[6]*100.         ; convert 1/cm to 1/m
  Ks=data[7]/8640000. ; convert cm/day to m/s

  IF Ks EQ 0 THEN BEGIN
     print, "ERROR, ks = 0"
     print, "data[7]=",data[7]
  ENDIF

  IF keyword_set(vg) THEN tmpb=2*(1.0-1.0/b) ELSE tmpb=b

  REFSMC=MAXSMC*(5.79E-9/Ks)^(1.0/(2.0*(tmpb+3)))
  REFSMC+=(MAXSMC-REFSMC)/3.0
  Ds=tmpb*Ks*(SATPSI/MAXSMC)

  ;; the table comes from the NOAH default SOILPARM.TBL file
  ;; index picks the correct soil type from the table
  F11=([-0.472,-1.044,-0.569,0.162,0.162,-0.327,-1.491,-1.118,-1.297,-3.209,-1.916,-2.138])[index]
  QTZ=([ 0.92 , 0.82 , 0.60 ,0.25 ,0.10 , 0.40 , 0.60 , 0.10 , 0.35 , 0.52 , 0.10 , 0.25 ])[index]

  openw, oun, /get, 'SOILPARM.TBL'
  printf, oun, 'Soil Parameters'
  printf, oun, 'STAS'
  printf, oun, '1,1   BB      DRYSMC      F11     MAXSMC   REFSMC   SATPSI  SATDK       SATDW     WLTSMC  QTZ    '
  printf, oun, ' 1,'+strjoin(strcompress([b,DRYSMC,F11,MAXSMC,REFSMC, $
                                          SATPSI,Ks,Ds,DRYSMC,QTZ]),',')
  close, oun
  free_lun, oun
END

;; run the noah model and rename the output file
PRO runNoah, i, texture
  print, "running noah "+texture+strcompress(i)
  spawn, "./noah >>outputfile"
  spawn, "mv out.* out"+"_"+strcompress(texture, /remove_all)+"_"+strcompress(fix(i), /remove_all)
END

;; for now this index is used into the arrays F11 and QTZ in writeSoilFile
FUNCTION textIndex, texture
  CASE strlowcase(texture) OF 
     "sand" :            index=0
     "loamy sand" :      index=1
     "sandy loam" :      index=2
     "silt loam" :       index=3
     "silt" :            index=4
     "loam" :            index=5
     "sandy clay loam" : index=6
     "silty clay loam" : index=7
     "clay loam" :       index=8
     "sandy clay" :      index=9
     "silty clay" :      index=10
     "clay" :            index=11
  ENDCASE
  return, index
END



PRO noahSoilsWTexture, texture, dataFile, textureFile, vg=vg, cam=cam, soilNum=soilNum
  print, ""
  print, texture
  print, ""

  IF NOT keyword_set(vg) THEN cam=1

  textData=getTextures(textureFile)
  textureIndex=textIndex(texture)
  line=''
  openr, un, /get, dataFile
  IF NOT keyword_set(soilNum) THEN BEGIN 
     openw, oun, /get, strcompress(texture, /remove_all)+".data"
     soilNum = -9999
  ENDIF

  readf, un, line ; read in the column headings and discard

  WHILE NOT eof(un) DO BEGIN
     readf, un, line
     data=float(strsplit(line, /extract))

     IF soilnum EQ -9999 OR soilNum EQ data[0] THEN BEGIN 
        IF (data[1] LT 999.9 AND data[1] NE 0 AND data[4] GT 0.25 AND $
            data[6] LT 10000000000 AND data[7] LT 10000000000 AND data[7] GT 0) THEN BEGIN
           
           dex=where(textData[0,*] EQ data[0])
           IF dex[0] NE -1 THEN BEGIN
              IF textData[1,dex] EQ texture OR $
                strcompress(textData[1,dex], /remove_all) $
                EQ strcompress(texture, /remove_all) $
                THEN BEGIN
                 IF soilNum NE -9999 THEN $
                   printf, oun, data, format='(20F15.5)'
                 writeSoilFile, data, textureIndex[0], cam=cam, vg=vg
                 runNoah, data[0], texture
              ENDIF  ; (if it is the texture we are currently looking at)
           ENDIF  ; (if we can find this soil in both databases)
        ENDIF  ; (if data is valid)
     ENDIF  ; (if soilnum eq data)   
  ENDWHILE

  IF soilNum NE -9999 THEN oun=un
  close, un, oun
  free_lun, un, oun
END