pro map_fv3_input, indate, datpath

;  Code to take raw gfs extracts and map into input file expected for forcing
;  the diurnal warming model

;  indate:  Date for which model run is desired, yymmdd
;  datpath:  path to root directory for gfs model inputs
;            (/satellite/data/gaw/skinsst/gfs/input)

;  FV3 version is a bit more complex in that taux and tauy are not
;  provided in the model output and must be derived from ustar and 
;  the wind speed components

;  Initializations

;  Original gfs grid characteristics
;  Grid is 0-360 in lon, 90 - -90 in lat
;  In revised, grid now has centers from 0 - 359.5, 90 - -89.5
;     Is a slight 'bug' where data from -89.75 to -90 is put in last bin
;     but these data not used

xorig = 720
yorig = 360

;  Define desired model grid

ntimes = 9     ; 4 times per day x 2 days + 1 ending

minlat = -60.0 ; still use this notation, but now understand offset
maxlat = 60.0
minlon = 0.0
maxlon = 360.0

res = 0.5

xnew = (maxlon - minlon)/res
;ynew = (maxlat - minlat)/res
ynew = (maxlat - minlat)/res + 1    ; for centers at 0, 60, and -60

;  Compute subgrid boundaries for model extract

istart = fix((minlon - 0.0) / res)
iend = fix((maxlon - 0.0) / res) - 1
jstart = fix((90.0 - maxlat) / res)     ; grid runs N-S
;jend = fix((90.0 - minlat) / res) - 1
jend = fix((90.0 - minlat) / res)

;  Initialize grid arrays

taux = fltarr(xnew,ynew,ntimes)
tauy = fltarr(xnew,ynew,ntimes)
Qe = fltarr(xnew,ynew,ntimes)
Qh = fltarr(xnew,ynew,ntimes)
Qlw = fltarr(xnew,ynew,ntimes)
Qsol = fltarr(xnew,ynew,ntimes)
sst = fltarr(xnew,ynew,ntimes)
ta = fltarr(xnew,ynew,ntimes)
qa = fltarr(xnew,ynew,ntimes)

Wper = fltarr(xnew,ynew,ntimes)
swh = fltarr(xnew,ynew,ntimes)
Wdir = fltarr(xnew,ynew,ntimes)

;  Additional arrays for new derivation of taux and tauy

p = fltarr(xnew,ynew)
ustar = fltarr(xnew,ynew)
uspd = fltarr(xnew,ynew)
vspd = fltarr(xnew,ynew)
ta_raw = fltarr(xnew,ynew)
qa_raw = fltarr(xnew,ynew)

lmask = fltarr(xnew,ynew)

fdate = fltarr(ntimes)

dum = -999.

taux(*,*,*) = dum
tauy(*,*,*) = dum
Qe(*,*,*) = dum
Qh(*,*,*) = dum
Qlw(*,*,*) = dum
Qsol(*,*,*) = dum
sst(*,*,*) = dum
ta(*,*,*) = dum
qa(*,*,*) = dum

Wper(*,*,*) = dum
swh(*,*,*) = dum
Wdir(*,*,*) = dum

;  Construct date strings for day before and day after

datestr = strtrim(string(indate),1)

yrstr = strmid(datestr,0,4)
monstr = strmid(datestr,4,2)
daystr = strmid(datestr,6,2)

yr=fix(yrstr)
mon=fix(monstr)
day=fix(daystr)

jday=julianday(day,mon,yr)

;  First data comes from day before

sjday = jday - 1

if (sjday eq 0) then begin
   syr = yr - 1
   smon = 12
   sday = 31
endif else begin
   calendarday,sjday,yr,sday,smon,leap
   syr = yr
endelse

sdatestr = string(syr,smon,sday,format='(i4,i2.2,i2.2)')
sdate = fix(sdatestr)

;  Last data comes from 00z on day after

ejday = jday + 1

leap=0
IF (((yr mod 4) eq 0 and (yr mod 100) ne 0) or (yr mod 400) eq 0) THEN leap = 1

if ((leap eq 0 and ejday eq 366) or (leap eq 1 and ejday eq 367)) then begin
   eyr = yr + 1
   emon = 1
   eday = 1
endif else begin
   calendarday,ejday,yr,eday,emon,eleap
   eyr = yr
endelse

edatestr = string(eyr,emon,eday,format='(i4,i2.2,i2.2)')
edate = fix(edatestr)


;  Enter loop through model times

for n = 1, ntimes do begin

;  Construct file root for given time

   CASE n OF

   1: BEGIN
      fdatestr = sdatestr
      init = '00'
      fdate(n-1) = sjday
      END

   2: BEGIN
      fdatestr = sdatestr
      init = '06'
      fdate(n-1) = sjday + 0.25
      END

   3: BEGIN
      fdatestr = sdatestr
      init = '12'
      fdate(n-1) = sjday + 0.5
      END

   4: BEGIN
      fdatestr = sdatestr
      init = '18'
      fdate(n-1) = sjday + 0.75
      END

   5: BEGIN
      fdatestr = datestr
      init = '00'
      fdate(n-1) = jday
      END

   6: BEGIN
      fdatestr = datestr
      init = '06'
      fdate(n-1) = jday + 0.25
      END

   7: BEGIN
      fdatestr = datestr
      init = '12'
      fdate(n-1) = jday + 0.5
      END

   8: BEGIN
      fdatestr = datestr
      init = '18'
      fdate(n-1) = jday + 0.75
      END

   9: BEGIN
      fdatestr = edatestr
      init = '00'
      fdate(n-1) = ejday
      END

   ENDCASE

   froot = datpath + '/' + fdatestr + '/gfs.t' + init + 'z.f00_'

;  Initialize temporary arrays for the current time

   p(*,*) = dum
   ustar(*,*) = dum
   uspd(*,*) = dum
   vspd(*,*) = dum
   ta_raw(*,*) = dum
   qa_raw(*,*) = dum

;  Sequence through each of the input product fields

;  Wind stress now computed at the end of all variables
;;  Wind stress

;   ;fname = froot + 'uflx.grd'
;   fname = froot + 'uflx.grd.gz'

;   ;openr,1,fname
;   openr,1,fname,/compress
;   dat=fltarr(xorig,yorig)
;   readu,1,dat
;   close,1

;   taux(*,*,n-1) = dat(istart:iend,jstart:jend)

;   ;fname = froot + 'vflx.grd'
;   fname = froot + 'vflx.grd.gz'

;   ;openr,1,fname
;   openr,1,fname,/compress
;   dat=fltarr(xorig,yorig)
;   readu,1,dat
;   close,1

;   tauy(*,*,n-1) = dat(istart:iend,jstart:jend)

;  Floating point arrays first

;  Ustar

   ;fname = froot + 'fricv.grd'
   fname = froot + 'fricv.grd.gz'

   ;openr,1,fname
   openr,1,fname,/compress
   dat=fltarr(xorig,yorig)
   readu,1,dat
   close,1

   ustar(*,*) = dat(istart:iend,jstart:jend)

;  Surface Pressure

   ;fname = froot + 'sp.grd'
   fname = froot + 'sp.grd.gz'

   ;openr,1,fname
   openr,1,fname,/compress
   dat=fltarr(xorig,yorig)
   readu,1,dat
   close,1

   p(*,*) = dat(istart:iend,jstart:jend)

;  Continue with those quantities stored as integer

;  Latent heat flux
   
   ;fname = froot + 'lhtfl.grd'
   fname = froot + 'lhtfl.grd.gz'
   scale = 0.1

   ;openr,1,fname
   openr,1,fname,/compress
   dat=intarr(xorig,yorig)
   readu,1,dat
   close,1

   Qe(*,*,n-1) = dat(istart:iend,jstart:jend) * scale

;  Sensible heat flux
   
   ;fname = froot + 'shtfl.grd'
   fname = froot + 'shtfl.grd.gz'
   scale = 0.1

   ;openr,1,fname
   openr,1,fname,/compress
   dat=intarr(xorig,yorig)
   readu,1,dat
   close,1

   Qh(*,*,n-1) = dat(istart:iend,jstart:jend) * scale

;  Longwave heat flux

   ;fname1 = froot + 'ulwrf.grd'
   ;fname2 = froot + 'dlwrf.grd'
   fname1 = froot + 'ulwrf.grd.gz'
   fname2 = froot + 'dlwrf.grd.gz'
   scale = 0.1

   ;openr,1,fname1
   openr,1,fname1,/compress
   dat1=intarr(xorig,yorig)
   readu,1,dat1
   close,1

   ;openr,1,fname2
   openr,1,fname2,/compress
   dat2=intarr(xorig,yorig)
   readu,1,dat2
   close,1

   Qlw(*,*,n-1) = (dat1(istart:iend,jstart:jend) - dat2(istart:iend,jstart:jend)) * scale

;  Solar heat flux

   ;fname1 = froot + 'uswrf.grd'
   ;fname2 = froot + 'dswrf.grd'
   fname1 = froot + 'uswrf.grd.gz'
   fname2 = froot + 'dswrf.grd.gz'
   scale = 0.1

   ;openr,1,fname1
   openr,1,fname1,/compress
   dat1=intarr(xorig,yorig)
   readu,1,dat1
   close,1

   ;openr,1,fname2
   openr,1,fname2,/compress
   dat2=intarr(xorig,yorig)
   readu,1,dat2
   close,1

   Qsol(*,*,n-1) = (dat1(istart:iend,jstart:jend) - dat2(istart:iend,jstart:jend)) * scale

;  sst
   
   ;fname = froot + 't.grd'
   fname = froot + 't.grd.gz'
   scale = 0.01

   ;openr,1,fname
   openr,1,fname,/compress
   dat=intarr(xorig,yorig)
   readu,1,dat
   close,1

   sst(*,*,n-1) = dat(istart:iend,jstart:jend) * scale - 273.15

;  Ta
   
   ;fname = froot + '2t.grd'
   fname = froot + '2t.grd.gz'
   scale = 0.01

   ;openr,1,fname
   openr,1,fname,/compress
   dat=intarr(xorig,yorig)
   readu,1,dat
   close,1

   ta(*,*,n-1) = dat(istart:iend,jstart:jend) * scale - 273.15
   ta_raw(*,*) = dat(istart:iend,jstart:jend) * scale

;  qa
   
   ;fname = froot + 'q.grd'
   fname = froot + 'q.grd.gz'
   scale = 0.1

   ;openr,1,fname
   openr,1,fname,/compress
   dat=intarr(xorig,yorig)
   readu,1,dat
   close,1

   qa(*,*,n-1) = dat(istart:iend,jstart:jend) * scale
   qa_raw(*,*) = dat(istart:iend,jstart:jend) / 10000.0

;  Significant wave height
   
   ;fname = froot + 'htsgw.grd'
   fname = froot + 'htsgw.grd.gz'
   scale = 0.01

   ;openr,1,fname
   openr,1,fname,/compress
   dat=intarr(xorig,yorig)
   readu,1,dat
   close,1

   swh(*,*,n-1) = dat(istart:iend,jstart:jend) * scale

;  Primary wave period
   
   ;fname = froot + 'perpw.grd'
   fname = froot + 'perpw.grd.gz'
   scale = 0.01

   ;openr,1,fname
   openr,1,fname,/compress
   dat=intarr(xorig,yorig)
   readu,1,dat
   close,1

   Wper(*,*,n-1) = dat(istart:iend,jstart:jend) * scale

;  Primary wave direction
   
   ;fname = froot + 'dirpw.grd'
   fname = froot + 'dirpw.grd.gz'
   scale = 0.1

   ;openr,1,fname
   openr,1,fname,/compress
   dat=intarr(xorig,yorig)
   readu,1,dat
   close,1

   Wdir(*,*,n-1) = dat(istart:iend,jstart:jend) * scale

;  U wind speed
   
   ;fname = froot + '10u.grd'
   fname = froot + '10u.grd.gz'
   scale = 0.01

   ;openr,1,fname
   openr,1,fname,/compress
   dat=intarr(xorig,yorig)
   readu,1,dat
   close,1

   uspd(*,*) = dat(istart:iend,jstart:jend) * scale

;  V wind speed
   
   ;fname = froot + '10v.grd'
   fname = froot + '10v.grd.gz'
   scale = 0.01

   ;openr,1,fname
   openr,1,fname,/compress
   dat=intarr(xorig,yorig)
   readu,1,dat
   close,1

   vspd(*,*) = dat(istart:iend,jstart:jend) * scale

;  Now, for FV3, derive taux and tauy for the current time

   rho = p / ((1.0 + 0.61*qa_raw) * ta_raw * 287.058)
   tau = rho * ustar^2

   wdspd = sqrt(uspd^2 + vspd^2)

   wdspd(where(wdspd lt 0.00001)) = 0.00001
   taux(*,*,n-1) = tau * uspd/wdspd
   tauy(*,*,n-1) = tau * vspd/wdspd

endfor     ; End loop through model times

;  Read in and map static land mask file

fname = datpath + '/gfs_lsm.grd.gz'
scale = 0.01

openr,1,fname,/compress
dat=intarr(xorig,yorig)
readu,1,dat
close,1

lmask(*,*) = dat(istart:iend,jstart:jend) * scale

;  Open files for storing sst initialization and model forcing data

outfluxf = datpath + '/' + datestr + '_flx.dat'
openw,1,outfluxf

outsstf = datpath + '/' + datestr + '_sst.dat'
openw,2,outsstf

;  All data now stored in output arrays.  Enter loop through grid points to 
;  output in desired format

for j = 0, ynew-1 do begin
   for i = 0, xnew-1 do begin

      ;  Check grid point for ocean coverage

;      if (lmask(i,j) lt 0.05) then begin
      if (lmask(i,j) lt 0.05 and sst(i,j,0) gt -2.0) then begin

         ;  Compute lat and lon for grid point

         ;alat = maxlat - j*res - res/2.
         ;alon = minlon + i*res + res/2.
         alat = maxlat - j*res
         alon = minlon + i*res

         ;  Write sst and initialization data

         printf,2,indate,alat,alon,sst(i,j,0),i,j, $
                  format='(i8,2(1x,f7.2),1x,f6.2,2(1x,i4))'

         ;  Write forcing data for corresponding grid points

         for n = 0, ntimes-1 do begin

            printf,1,fdate(n),taux(i,j,n),tauy(i,j,n),Qsol(i,j,n),Qlw(i,j,n), $
                  Qe(i,j,n),Qh(i,j,n),swh(i,j,n),Wper(i,j,n),Wdir(i,j,n), $
		  sst(i,j,n),ta(i,j,n),qa(i,j,n), $
		  format='(f6.2,2(1x,f10.7),1x,f8.2,3(1x,f7.2),2(1x,f7.2),1x,f6.1,2(1x,f6.2),1x,f6.1)'
		  ;format='(f6.2,2(1x,f10.7),1x,f8.2,3(1x,f7.2),2(1x,f7.2),1x,f6.1)'
		  ;format='(f6.2,2(1x,f10.7),4(1x,f7.2),2(1x,f7.2),1x,f6.1)'
		  ;format='(f6.2,2(1x,f10.7),4(1x,f9.4),2(1x,f9.4))'

         endfor

      endif     ; End of if was ocean grid point

   endfor
endfor

end