program gefs_bias_g2     
!
! main program: gefs_bias_g2       
!
! prgmmr: Bo Cui           org: np/wx20        date: 2006-12-10
!                          mod np/wx20         date: 2008-01-25
!                          mod np/wx20         date: 2009-04_01
!                          mod np/wx20         date: 2011-04_01
!                          mod: Bo.Cui         date: 2015-06-01
!
! abstract: 1. update bias estimation between GDAS analysis and NCEP ensemble forecast 
!              add bias estimation calculation between NCEP and CMC analysis (CMC - NCEP)
!           2. add more variables for bias estimation (+14)
!           3. add variable ULWRF(OLR). ULWRF(OLR) is 6hr average value and ULWRF(OLR) at GEFS control
!              f00hr is almost instantaneous value. Therefore, choose 6hr forecast of previous cycle as
!              the ULWRF(OLR) analysis
!           4. add variable 2m dew point temperature and relative humility, calculate & output bias
!              calculate dew point temperature forecast from relative humidity and 2m temperature
!              then calculate the bias of 2m dew point temperature, 52 variables in total for bias estimation
!              51 variables for reanalysis difference (no 2m dpt)
!              49 variables for NCEP CMC analysis difference (no tmax,tmin and vvel)

! modification: add new variable TCDC (total cloud cover), 6 hour average
!              
! usage:
!
!   input file: grib
!     unit 11 -    : prior bias estimation                                               
!     unit 12 -    : analysis
!     unit 13 -    : analysis for Tmax, Tmin and ULWRF
!     unit 14 -    : ncep operational forecast or CDAS reanalysis
!     unit 15 -    : CFS reanalysis for variable tmax, tmin, ulwrf(surface) and ulwrf(top)
!
!   output file: grib
!     unit 51 -    : updated bias estimation pgrba file
!
!   parameters
!     fgrid -      : ensemble forecast
!     agrid -      : analysis data (GDAS)
!     rgrid -      : reanalysis data (CDAS)
!     bias  -      : bias estimation
!     dec_w -      : decay averaging weight 
!     nvar  -      : number of variables

! programs called:
!   baopenr          grib i/o
!   baopenw          grib i/o
!   baclose          grib i/o
!   getgb2           grib2 reader
!   putgb2           grib2 writer
!   init_parm        define grid definition and product definition
!   printinfr        print grib2 data information
!
! attributes:
!   language: fortran 90
!
!$$$

use grib_mod
use params
use naefs_mod

implicit none

integer      ivar,i,k,icstart,odate,ij
!parameter   (nvar=52)

real,       allocatable :: agrid(:),fgrid(:),bias(:),rgrid(:)
real,       allocatable :: rh2m(:),tmp2m(:)
real,       allocatable :: anl_ncep(:),anl_cmc(:)

real        dec_w

integer     ifile,afile,afile_m06,cfile,rfile,rfile_m06,rfile_m12,ofile

integer     maxgrd,ndata,index,j,n,iret,jret             
integer     ipd11_new(nvar),ipd12_new(nvar)              
integer     ipd11_cfs(nvar),ipd12_cfs(nvar)              

character*7 cfortnn
character*4 nens

namelist/message/icstart,nens,dec_w,odate

read(5,message,end=1020)
write(6,message)

if(nens.ne.'mdf'.and.nens.ne.'anl') then

ifile=11
afile=12
afile_m06=13
cfile=14
ofile=51

! set the fort.* of intput files

write(cfortnn,'(a5,i2)') 'fort.',afile
call baopenr(afile,cfortnn,iret)
if (iret .ne. 0) then; print*,'there is no analysis data, stop!'; endif
if (iret .ne. 0) goto 1020

write(cfortnn,'(a5,i2)') 'fort.',afile_m06
call baopenr(afile_m06,cfortnn,iret)
if (iret .ne. 0) then; print*,'there is no analysis data for Tmax and Tmin, stop!'; endif
if (iret .ne. 0) goto 1020

write(cfortnn,'(a5,i2)') 'fort.',cfile
call baopenr(cfile,cfortnn,iret)
if (iret .ne. 0) then; print*,'there is no NCEP forecast data, stop!'; endif
if (iret .ne. 0) goto 1020

if(icstart.eq.0) then
  write(cfortnn,'(a5,i2)') 'fort.',ifile
  call baopenr(ifile,cfortnn,iret)
  if (iret .ne. 0) then; print*,'there is no bias estimation data, please check!'; endif
endif

! set the fort.* of output file

write(cfortnn,'(a5,i2)') 'fort.',ofile
call baopenw(ofile,cfortnn,iret)
if (iret .ne. 0) then; print*,'there is no output bias data, stop!'; endif
if (iret .ne. 0) goto 1020

! find grib message, maxgrd: number of grid points in the defined grid

ipdt=-9999; igdt=-9999
ipdtn=-1; igdtn=-1
call init_parm(ipdtn,ipdt,igdtn,igdt)
call getgb2(cfile,0,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt,unpack,jskp,gfld,iret)
maxgrd=gfld%ngrdpts
call gf_free(gfld)

if (iret .ne. 0) then; print*,' getgbeh ,afile,index,iret =',afile,index,iret; endif
if (iret .ne. 0) goto 1020

! get NCEP ensemble ipdt message: fixed surface and its scaled value

call getipdt_g2_surface(cfile,ipd11_new,ipd12_new)

ipd11=ipd11_new
ipd12=ipd12_new

allocate (agrid(maxgrd),fgrid(maxgrd),bias(maxgrd))
allocate (rh2m(maxgrd),tmp2m(maxgrd))

do ivar = 1, nvar  

  bias=0.0            
  agrid=-9999.9999
  fgrid=-9999.9999

  ! read and process variable of input data

  ipdt=-9999
  igdt=-9999

  ipdt(1)=ipd1(ivar)
  ipdt(2)=ipd2(ivar)
  ipdt(10)=ipd10(ivar)
  ipdt(11)=ipd11(ivar)
  ipdt(12)=ipd12(ivar)

  ! get initialized bias estimation

  if(icstart.eq.1) then
    print *, '----- Cold Start for Bias Estimation -----'
    print*, '  '
    bias=0.0
  else
    print *, '----- Initialized Bias for Current Time -----'

    igdtn=-1
    if(nens.eq.'gfs') ipdtn=ipdn(ivar)
    if(nens.eq.'avg') ipdtn=ipdnm(ivar)
    call init_parm(ipdtn,ipdt,igdtn,igdt)
    call getgb2(ifile,0,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt,unpack,jskp,gfld,iret)

    if(iret.ne.0) then
      print*, 'There is no variable', jpdt(1),jpdt(2),jpdt(10),jpdt(12)
      call gf_free(gfld)
    endif

    if(iret .eq. 0) then    
      bias(1:maxgrd)=gfld%fld(1:maxgrd)
      call printinfr(gfld,ivar)
      call gf_free(gfld)
    else
      bias=0.0
    endif
  endif

  ! get ensemble forecast

  print *, '----- NCEP ensemble forecast for current Time ------'

  ! ipdn for GFS forecast; ipdnm for ensemble average forecast

  igdtn=-1
  if(nens.eq.'gfs') ipdtn=ipdn(ivar)
  if(nens.eq.'avg') ipdtn=ipdnm(ivar)

  call init_parm(ipdtn,ipdt,igdtn,igdt)
  if(ipd1(ivar).eq.0.and.ipd2(ivar).eq.6.and.ipd10(ivar).eq.103.and.ipd12(ivar).eq.2) then
    print *, ' '
    print*, 'Start to Calcualte GEFS Dew Point Temperature Forecast '; print *, ' '
    call get_dpt_g2(cfile,maxgrd,ipdtn,ipdt,igdtn,igdt,gfldo,iret)
    print *, ' '
  elseif(ipd1(ivar).eq.2.and.ipd2(ivar).eq.1.and.ipd10(ivar).eq.103.and.ipd12(ivar).eq.10) then
    print *, ' '
    print*, 'Start to Calcualte GEFS 10m Wind Speed '; print *, ' '
    call get_wspd10m(cfile,maxgrd,ipdtn,ipdt,igdtn,igdt,gfldo,iret)
  else
    call getgb2(cfile,0,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt,unpack,jskp,gfldo,iret)
  endif

  if(iret.ne.0) then
    print*, 'There is no variable', jpdt(1),jpdt(2),jpdt(10),jpdt(12)
    print *, ' '
    call gf_free(gfld)
  endif
     
  if(iret.ne.0) goto 100

  if(iret.eq.0) then
    call printinfr(gfldo,ivar)
    fgrid(1:maxgrd)=gfldo%fld(1:maxgrd)
  endif

  ! get analysis data

  print *, '----- Analysis for Current Time ------'

  igdtn=-1; ipdtn=ipdn(ivar)

  call init_parm(ipdtn,ipdt,igdtn,igdt)

  if(ipd1(ivar).eq.0.and.ipd2(ivar).eq.6.and.ipd10(ivar).eq.103.and.ipd12(ivar).eq.2) then
    print *, ' '
    print*, 'Start to Calcualte GEFS Dew Point Temperature Analysis '; print *, ' '
    call get_dpt_g2(afile,maxgrd,ipdtn,ipdt,igdtn,igdt,gfld,iret)
  elseif(ipd1(ivar).eq.2.and.ipd2(ivar).eq.1.and.ipd10(ivar).eq.103.and.ipd12(ivar).eq.10) then
    print *, ' '
    print*, 'Start to Calcualte GEFS 10m Wind Speed Analysis '; print *, ' '
    call get_wspd10m(afile,maxgrd,ipdtn,ipdt,igdtn,igdt,gfld,iret)
  elseif(ipd1(ivar).eq.0.and.ipd2(ivar).eq.4.and.ipd10(ivar).eq.103.and.ipd12(ivar).eq.2) then
    call getgb2(afile_m06,0,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt,unpack,jskp,gfld,iret)
  elseif(ipd1(ivar).eq.0.and.ipd2(ivar).eq.5.and.ipd10(ivar).eq.103.and.ipd12(ivar).eq.2) then
    call getgb2(afile_m06,0,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt,unpack,jskp,gfld,iret)
  elseif(ipd1(ivar).eq.5.and.ipd2(ivar).eq.193.and.ipd10(ivar).eq.1.and.ipd12(ivar).eq.0) then
    call getgb2(afile_m06,0,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt,unpack,jskp,gfld,iret)
  elseif(ipd1(ivar).eq.5.and.ipd2(ivar).eq.193.and.ipd10(ivar).eq.8.and.ipd12(ivar).eq.0) then
    call getgb2(afile_m06,0,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt,unpack,jskp,gfld,iret)
  elseif(ipd1(ivar).eq.6.and.ipd2(ivar).eq.1.and.ipd10(ivar).eq.10.and.ipd12(ivar).eq.0) then
    call getgb2(afile_m06,0,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt,unpack,jskp,gfld,iret)
  else
    call getgb2(afile,0,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt,unpack,jskp,gfld,iret)
  endif

  if(iret.ne.0) then
    print*, 'There is no variable', jpdt(1),jpdt(2),jpdt(10),jpdt(12)
    print *, ' '
    call gf_free(gfld)
    bias(1:maxgrd)=0.0                           
  endif

! if(iret.ne.0) goto 200
  if(iret.ne.0) goto 100

  if(iret.eq.0) then
    call printinfr(gfld,ivar)
    agrid(1:maxgrd)=gfld%fld(1:maxgrd)
  endif

  call gf_free(gfld)

  ! apply the decay average

  call decay(bias,fgrid,agrid,maxgrd,dec_w)

  ! save the first moment bias, give the bias a time one day before

  200 continue

  ! output bias estimation

  gfldo%fld(1:maxgrd)=bias(1:maxgrd)

  print *, '----- Output Bias Estimation for Current Time ------'

  ! gfldo%idrtmpl(3) : GRIB2 DRT 5.40 decimal scale factor

  if(gfldo%ipdtmpl(1).eq.3.and.gfldo%ipdtmpl(2).eq.5) then
    continue
  elseif(gfldo%ipdtmpl(1).eq.3.and.gfldo%ipdtmpl(2).eq.1) then
    gfldo%idrtmpl(3)=2
  elseif(gfldo%ipdtmpl(1).eq.3.and.gfldo%ipdtmpl(2).eq.0) then
    gfldo%idrtmpl(3)=2
  else
    gfldo%idrtmpl(3)=2
  endif

  gfldo%ipdtmpl(3)=6                !  gfldo%ipdtmpl(3)=6 GRIB2 Code Table 4.3 forecast error

  call putgb2(ofile,gfldo,jret)
  call printinfr(gfldo,ivar)

  call gf_free(gfldo)

  100 continue

! end of bias estimation 

enddo

call baclose(ifile,iret)
call baclose(afile,iret)
call baclose(afile_m06,iret)
call baclose(cfile,iret)
call baclose(ofile,iret)

print *,'Bias Estimation Successfully Complete'

endif

if(nens.eq.'mdf') then

ifile=11
afile=12
afile_m06=13
rfile=14
rfile_m06=15
ofile=51

! set the fort.* of intput files

write(cfortnn,'(a5,i2)') 'fort.',afile
call baopenr(afile,cfortnn,iret)
if (iret .ne. 0) then; print*,'there is no analysis data, stop!';  endif
if (iret .ne. 0) goto 1020

write(cfortnn,'(a5,i2)') 'fort.',afile_m06
call baopenr(afile_m06,cfortnn,iret)
if (iret .ne. 0) then; print*,'there is no analysis data for Tmax and Tmin, stop!'; endif
if (iret .ne. 0) goto 1020

write(cfortnn,'(a5,i2)') 'fort.',rfile
call baopenr(rfile,cfortnn,iret)
if (iret .ne. 0) then; print*,'there is no reanalysis data, stop!'; endif
if (iret .ne. 0) goto 1020

write(cfortnn,'(a5,i2)') 'fort.',rfile_m06
call baopenr(rfile_m06,cfortnn,iret)
if (iret .ne. 0) then; print*,'there is no reanalysis data for surface variables, stop!'; endif
if (iret .ne. 0) goto 1020

if(icstart.eq.0) then
  write(cfortnn,'(a5,i2)') 'fort.',ifile
  call baopenr(ifile,cfortnn,iret)
  if (iret .ne. 0) then; print*,'there is no bias between CDAS and GDAS, please check!'; endif
endif

! set the fort.* of output file

write(cfortnn,'(a5,i2)') 'fort.',ofile
call baopenw(ofile,cfortnn,iret)
if (iret .ne. 0) then; print*,'there is no output bias data, stop!'; endif
if (iret .ne. 0) goto 1020

! find grib message, maxgrd: number of grid points in the defined grid

ipdt=-9999; igdt=-9999
ipdtn=-1; igdtn=-1
call init_parm(ipdtn,ipdt,igdtn,igdt)
call getgb2(afile,0,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt,unpack,jskp,gfld,iret)
maxgrd=gfld%ngrdpts
call gf_free(gfld)

if (iret .ne. 0) then; print*,' getgbeh ,afile,index,iret =',afile,index,iret; endif
if (iret .ne. 0) goto 1020

! get NCEP analysis ipdt message: fixed surface and its scaled value

call getipdt_g2_surface(afile,ipd11_new,ipd12_new)

print*,' ipd11_new,ipd12_new=',ipd11_new,ipd12_new

ipd11=ipd11_new
ipd12=ipd12_new

! get CFS analysis ipdt message: fixed surface and its scaled value

call getipdt_g2_surface(rfile,ipd11_new,ipd12_new)

ipd11_cfs=ipd11_new
ipd12_cfs=ipd12_new

print*,' ipd11_new,ipd12_new=',ipd11_new,ipd12_new

allocate (agrid(maxgrd),rgrid(maxgrd),bias(maxgrd))

do ivar = 1, nvar  

  bias=0.0             
  agrid=-9999.9999
  rgrid=-9999.9999

  ! read and process variable of input data

  ipdt=-9999
  igdt=-9999

  ipdtn=ipdn(ivar)
  ipdt(1)=ipd1(ivar)
  ipdt(2)=ipd2(ivar)
  ipdt(10)=ipd10(ivar)
  ipdt(11)=ipd11(ivar)
  ipdt(12)=ipd12(ivar)

  ! get initialized bias estimation between CDAS and GDAS

  if(icstart.eq.1) then
    print *, '----- Cold Start for Bias Estimation between CDAS and GDAS -----'
    print*, '  '
    bias=0.0
  else
    print *, '----- Initialized Bias for Current Time -----'

    igdtn=-1
    ipdtn=ipdnr(ivar)
!   ipdtn=ipdn(ivar)
    call init_parm(ipdtn,ipdt,igdtn,igdt)
    call getgb2(ifile,0,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt,unpack,jskp,gfld,iret)

    if (iret.ne.0) print*, 'There is no variable', jpdt(1),jpdt(2),jpdt(10),jpdt(12)
    if (iret .eq. 0) then
      bias(1:maxgrd)=gfld%fld(1:maxgrd)
      call printinfr(gfld,ivar)
      call gf_free(gfld)
    else
      bias=0.0
    endif

  endif

  ! get analysis data GDAS

  print *, '----- Analysis for Current Time ------'

  igdtn=-1
  ipdtn=ipdn(ivar) 
  call init_parm(ipdtn,ipdt,igdtn,igdt)

  if(ipd1(ivar).eq.0.and.ipd2(ivar).eq.6.and.ipd10(ivar).eq.103.and.ipd12(ivar).eq.2) then
    print *, ' '
    print*, 'Start to Calcualte GEFS Dew Point Temperature Analysis '; print *, ' '
    call get_dpt_g2(afile,maxgrd,ipdtn,ipdt,igdtn,igdt,gfld,iret)
  elseif(ipd1(ivar).eq.2.and.ipd2(ivar).eq.1.and.ipd10(ivar).eq.103.and.ipd12(ivar).eq.10) then
    print *, ' '
    print*, 'Start to Calcualte GEFS 10m Wind Speed Analysis '; print *, ' '
    call get_wspd10m(afile,maxgrd,ipdtn,ipdt,igdtn,igdt,gfld,iret)
  elseif(ipd1(ivar).eq.0.and.ipd2(ivar).eq.4.and.ipd10(ivar).eq.103.and.ipd12(ivar).eq.2) then
    call getgb2(afile_m06,0,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt,unpack,jskp,gfld,iret)
  elseif(ipd1(ivar).eq.0.and.ipd2(ivar).eq.5.and.ipd10(ivar).eq.103.and.ipd12(ivar).eq.2) then
    call getgb2(afile_m06,0,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt,unpack,jskp,gfld,iret)
  elseif(ipd1(ivar).eq.5.and.ipd2(ivar).eq.193.and.ipd10(ivar).eq.1.and.ipd12(ivar).eq.0) then
    call getgb2(afile_m06,0,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt,unpack,jskp,gfld,iret)
  elseif(ipd1(ivar).eq.5.and.ipd2(ivar).eq.193.and.ipd10(ivar).eq.8.and.ipd12(ivar).eq.0) then
    call getgb2(afile_m06,0,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt,unpack,jskp,gfld,iret)
  else
    call getgb2(afile,0,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt,unpack,jskp,gfld,iret)
  endif

  if (iret.ne.0) then
    print*, 'There is no variable', jpdt(1),jpdt(2),jpdt(10),jpdt(12)
    call gf_free(gfld)
  endif

  if (iret.ne.0) goto 300
  if (iret.eq.0) then
    call printinfr(gfld,ivar)
    agrid(1:maxgrd)=gfld%fld(1:maxgrd)
  endif

  call gf_free(gfld)

  ! get CFS reanalysis data CDAS  

  igdtn=-1
  ipdtn=ipdnr(ivar) 
  ipdt(11)=ipd11_cfs(ivar)
  ipdt(12)=ipd12_cfs(ivar)
  call init_parm(ipdtn,ipdt,igdtn,igdt)

  print *, '----- CDAS reanalysis for current Time ------'

  if(ipd1(ivar).eq.0.and.ipd2(ivar).eq.6.and.ipd10(ivar).eq.103.and.ipd12(ivar).eq.2) then
    print *, ' '
    print*, 'Start to Calcualte CDAS Dew Point Temperature Forecast '; print *, ' '
    call get_dpt_g2(rfile,maxgrd,ipdtn,ipdt,igdtn,igdt,gfldo,iret)
  elseif(ipd1(ivar).eq.2.and.ipd2(ivar).eq.1.and.ipd10(ivar).eq.103.and.ipd12(ivar).eq.10) then
    print *, ' '
    print*, 'Start to Calcualte CFS 10m Wind Speed Analysis '; print *, ' '
    call get_wspd10m(rfile,maxgrd,ipdtn,ipdt,igdtn,igdt,gfldo,iret)
  elseif(ipd1(ivar).eq.0.and.ipd2(ivar).eq.4.and.ipd10(ivar).eq.103.and.ipd12(ivar).eq.2) then
    call getgb2(rfile_m06,0,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt,unpack,jskp,gfldo,iret)
  elseif(ipd1(ivar).eq.0.and.ipd2(ivar).eq.5.and.ipd10(ivar).eq.103.and.ipd12(ivar).eq.2) then
    call getgb2(rfile_m06,0,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt,unpack,jskp,gfldo,iret)
  elseif(ipd1(ivar).eq.5.and.ipd2(ivar).eq.193.and.ipd10(ivar).eq.1.and.ipd12(ivar).eq.0) then
    call getgb2(rfile_m06,0,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt,unpack,jskp,gfldo,iret)
  elseif(ipd1(ivar).eq.5.and.ipd2(ivar).eq.193.and.ipd10(ivar).eq.8.and.ipd12(ivar).eq.0) then
    call getgb2(rfile_m06,0,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt,unpack,jskp,gfldo,iret)
  else
    call getgb2(rfile,0,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt,unpack,jskp,gfldo,iret)
  endif

  if (iret.ne.0) then
    print*, 'There is no variable', jpdt(1),jpdt(2),jpdt(10),jpdt(12)
    call gf_free(gfldo)
  endif

  if (iret.ne.0) goto 300
  if (iret.eq.0) then
    call printinfr(gfldo,ivar)
    rgrid(1:maxgrd)=gfldo%fld(1:maxgrd)
  endif

  ! apply the decay average

  call decay(bias,rgrid,agrid,maxgrd,dec_w)

  gfldo%fld(1:maxgrd)=bias(1:maxgrd)

  print *, '----- Output Bias Estimation between CDAS and GDAS ------'

  ! gfldo%idrtmpl(3) : GRIB2 DRT 5.40 decimal scale factor

  if(gfldo%ipdtmpl(1).eq.3.and.gfldo%ipdtmpl(2).eq.5) then
    continue
  elseif(gfldo%ipdtmpl(1).eq.3.and.gfldo%ipdtmpl(2).eq.1) then
    gfldo%idrtmpl(3)=2
  elseif(gfldo%ipdtmpl(1).eq.3.and.gfldo%ipdtmpl(2).eq.0) then
    gfldo%idrtmpl(3)=2
  else
    gfldo%idrtmpl(3)=2
  endif

  gfldo%ipdtmpl(3)=7                !  gfldo%ipdtmpl(3)=7 GRIB2 Code Table 4.3 analysis error

  call putgb2(ofile,gfldo,iret)
  call printinfr(gfldo,ivar)

  call gf_free(gfldo)

  300 continue

! end of bias estimation between CDAS and GDAS 

enddo

call baclose(ifile,iret)
call baclose(afile,iret)
call baclose(afile_m06,iret)
call baclose(rfile,iret)
call baclose(rfile_m06,iret)
call baclose(ofile,iret)

print *,'Bias Estimation Between CDAS and GDAS Successfully Complete'

endif

if(nens.eq.'anl') then

ifile=11
afile=12
afile_m06=13
rfile=14
rfile_m12=15
ofile=51

! set the fort.* of intput files

write(cfortnn,'(a5,i2)') 'fort.',afile
call baopenr(afile,cfortnn,iret)
if (iret .ne. 0) then; print*,'there is no NCEP analysis data, stop!';  endif
if (iret .ne. 0) goto 1020

write(cfortnn,'(a5,i2)') 'fort.',afile_m06
call baopenr(afile_m06,cfortnn,iret)
if (iret .ne. 0) then; print*,'there is no NCEP analysis for  ULWRF!';  endif
if (iret .ne. 0) goto 1020

write(cfortnn,'(a5,i2)') 'fort.',rfile
call baopenr(rfile,cfortnn,iret)
if (iret .ne. 0) then; print*,'there is no CMC analysis data, stop!'; endif
if (iret .ne. 0) goto 1020

write(cfortnn,'(a5,i2)') 'fort.',rfile_m12
call baopenr(rfile_m12,cfortnn,iret)
if (iret .ne. 0) then; print*,'there is no CMC analysis for ULWRF, stop!'; endif
!if (iret .ne. 0) goto 1020

if(icstart.eq.0) then
  write(cfortnn,'(a5,i2)') 'fort.',ifile
  call baopenr(ifile,cfortnn,iret)
  if (iret .ne. 0) then; print*,'there is no bias between CMC and NCEP analysis, please check!'; endif
endif

! set the fort.* of output file

write(cfortnn,'(a5,i2)') 'fort.',ofile
call baopenw(ofile,cfortnn,iret)
if (iret .ne. 0) then; print*,'there is no output bias data, stop!'; endif
if (iret .ne. 0) goto 1020

! find grib message, maxgrd: number of grid points in the defined grid

ipdt=-9999; igdt=-9999
ipdtn=-1; igdtn=-1
call init_parm(ipdtn,ipdt,igdtn,igdt)
call getgb2(afile,0,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt,unpack,jskp,gfld,iret)
maxgrd=gfld%ngrdpts
call gf_free(gfld)

if (iret .ne. 0) then; print*,' getgbeh ,afile,index,iret =',afile,index,iret; endif
if (iret .ne. 0) goto 1020

! get NCEP analysis ipdt message: fixed surface and its scaled value

call getipdt_g2_surface(afile,ipd11_new,ipd12_new)

ipd11=ipd11_new
ipd12=ipd12_new

! get CMC analysis ipdt message: fixed surface and its scaled value

call getipdt_g2_surface(rfile,ipd11_new,ipd12_new)

ipd11_cmc=ipd11_new
ipd12_cmc=ipd12_new

allocate (anl_ncep(maxgrd),anl_cmc(maxgrd),bias(maxgrd))
allocate (tmp2m(maxgrd),rh2m(maxgrd))

do ivar = 1, nvar  

  bias=0.0             
  anl_cmc=-9999.9999
  anl_ncep=-9999.9999

  ! read and process variable of input data

  ipdt=-9999
  igdt=-9999

  ipdt(1)=ipd1(ivar)
  ipdt(2)=ipd2(ivar)
  ipdt(10)=ipd10(ivar)
  ipdt(11)=ipd11(ivar)
  ipdt(12)=ipd12(ivar)

  ! get initialized bias estimation between NCEP and CMC analysis

  if(icstart.eq.1) then
    print *, '----- Cold Start for Bias Estimation between NCEP and CMC analysis -----'
    print*, '  '
    bias=0.0
  else

    print *, '----- Initialized Bias for Current Time -----'

    igdtn=-1
    ipdtn=ipdn(ivar)
    call init_parm(ipdtn,ipdt,igdtn,igdt)
    call getgb2(ifile,0,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt,unpack,jskp,gfld,iret)

    if (iret.ne.0) then
      print*, 'There is no variable', jpdt(1),jpdt(2),jpdt(10),jpdt(12)
      call gf_free(gfld)
    endif

    if (iret .eq. 0) then
      bias(1:maxgrd)=gfld%fld(1:maxgrd)
      call printinfr(gfld,ivar)
      call gf_free(gfld)
    else
      bias=0.0
    endif

  endif

  ! get NCEP analysis data GDAS

  igdtn=-1; ipdtn=ipdn(ivar)
  call init_parm(ipdtn,ipdt,igdtn,igdt)

  ! there are no Tmax and Tmin

  if(ipd1(ivar).eq.0.and.ipd2(ivar).eq.4.and.ipd10(ivar).eq.103.and.ipd12(ivar).eq.2) goto 500
  if(ipd1(ivar).eq.0.and.ipd2(ivar).eq.5.and.ipd10(ivar).eq.103.and.ipd12(ivar).eq.2) goto 500

  print *, '----- NCEP Analysis for Current Cycle ------'

  if(ipd1(ivar).eq.0.and.ipd2(ivar).eq.6.and.ipd10(ivar).eq.103.and.ipd12(ivar).eq.2) then
    print *, ' '
    print*, 'Start to Calcualte NCEP Dew Point Temperature Analysis '; print *, ' '
    call get_dpt_g2(afile,maxgrd,ipdtn,ipdt,igdtn,igdt,gfldo,iret)
  elseif(ipd1(ivar).eq.2.and.ipd2(ivar).eq.1.and.ipd10(ivar).eq.103.and.ipd12(ivar).eq.10) then
    print *, ' '
    print*, 'Start to Calcualte GEFS 10m Wind Speed '; print *, ' '
    call get_wspd10m(afile,maxgrd,ipdtn,ipdt,igdtn,igdt,gfldo,iret)
  elseif(ipd1(ivar).eq.5.and.ipd2(ivar).eq.193.and.ipd10(ivar).eq.1.and.ipd12(ivar).eq.0) then
    call getgb2(afile_m06,0,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt,unpack,jskp,gfldo,iret)
  elseif(ipd1(ivar).eq.5.and.ipd2(ivar).eq.193.and.ipd10(ivar).eq.8.and.ipd12(ivar).eq.0) then
    call getgb2(afile_m06,0,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt,unpack,jskp,gfldo,iret)
  else
    call getgb2(afile,0,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt,unpack,jskp,gfldo,iret)
  endif

  if (iret.ne.0) then
    print*, 'There is no variable', jpdt(1),jpdt(2),jpdt(10),jpdt(12)
    call gf_free(gfldo)
  endif

  if (iret.ne.0) goto 500

  if (iret.eq.0) then
    call printinfr(gfldo,ivar)
    anl_ncep(1:maxgrd)=gfldo%fld(1:maxgrd)
  endif

  ! get CMC analysis data
  ! there is no two ULWRF(surface/top) and VVEL(850mb) for analysis, get from ens average
  ! ULWRF(surface): PRODUCT TEMPLATE 4. 12 :  5 193 4 0 70 0 0 11 1 1 0 0 
  ! ULWRF(top):     PRODUCT TEMPLATE 4. 12 :  5 193 4 0 70 0 0 11 1 8 0 0

  ipdt(11)=ipd11_cmc(ivar)
  ipdt(12)=ipd12_cmc(ivar)

  igdtn=-1; ipdtn=ipdn(ivar)
  call init_parm(ipdtn,ipdt,igdtn,igdt)

  print *, '----- CMC Analysis for Current Cycle ------'

  if(ipd1(ivar).eq.0.and.ipd2(ivar).eq.6.and.ipd10(ivar).eq.103.and.ipd12(ivar).eq.2) then
    print *, ' '
    print*, 'Start to Calcualte CMC Dew Point Temperature Analysis '; print *, ' '
    call get_dpt_g2(rfile,maxgrd,ipdtn,ipdt,igdtn,igdt,gfld,iret)
  elseif(ipd1(ivar).eq.2.and.ipd2(ivar).eq.1.and.ipd10(ivar).eq.103.and.ipd12(ivar).eq.10) then
    print *, ' '
    print*, 'Start to Calcualte CMC 10m Wind Speed '; print *, ' '
    call get_wspd10m(rfile,maxgrd,ipdtn,ipdt,igdtn,igdt,gfld,iret)
  elseif(ipd1(ivar).eq.5.and.ipd2(ivar).eq.193.and.ipd10(ivar).eq.1.and.ipd12(ivar).eq.0) then
    ipdtn=12
    call init_parm(ipdtn,ipdt,igdtn,igdt)
    call getgb2(rfile_m12,0,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt,unpack,jskp,gfld,iret)
  elseif(ipd1(ivar).eq.5.and.ipd2(ivar).eq.193.and.ipd10(ivar).eq.8.and.ipd12(ivar).eq.0) then
    ipdtn=12
    call init_parm(ipdtn,ipdt,igdtn,igdt)
    call getgb2(rfile_m12,0,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt,unpack,jskp,gfld,iret)
  elseif(ipd1(ivar).eq.6.and.ipd2(ivar).eq.1.and.ipd10(ivar).eq.10.and.ipd12(ivar).eq.0) then
    print *, ' '
    print*, 'Start to get CMC TCDC Analysis '; print *, ' '
    ipdtn=1 
    call init_parm(ipdtn,ipdt,igdtn,igdt)
    call getgb2(rfile,0,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt,unpack,jskp,gfld,iret)
  else
    call getgb2(rfile,0,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt,unpack,jskp,gfld,iret)
  endif

  if (iret.ne.0) then
    print*, 'There is no variable', jpdt(1),jpdt(2),jpdt(10),jpdt(12)
    call gf_free(gfld)
  endif

  if (iret.ne.0) goto 500

  if (iret.eq.0) then
    call grid_cnvncep_g2(gfld,ivar)
    anl_cmc(1:maxgrd)=gfld%fld(1:maxgrd)
  endif

  call gf_free(gfld)

  ! apply the decay average

  call decay(bias,anl_cmc,anl_ncep,maxgrd,dec_w)

  print *, '----- Output Bias Estimation between CMC and NCEP analysis ------'

  ! gfldo%idrtmpl(3) : GRIB2 DRT 5.40 decimal scale factor

  if(gfldo%ipdtmpl(1).eq.3.and.gfldo%ipdtmpl(2).eq.5) then
    continue
  elseif(gfldo%ipdtmpl(1).eq.3.and.gfldo%ipdtmpl(2).eq.1) then
    gfldo%idrtmpl(3)=2
  elseif(gfldo%ipdtmpl(1).eq.3.and.gfldo%ipdtmpl(2).eq.0) then
    gfldo%idrtmpl(3)=2
  else
    gfldo%idrtmpl(3)=2
  endif

  gfldo%fld(1:maxgrd)=bias(1:maxgrd)

  gfldo%ipdtmpl(3)=7                !  gfldo%ipdtmpl(3)=7 GRIB2 Code Table 4.3 analysis error

  call putgb2(ofile,gfldo,jret)

  call printinfr(gfldo,ivar)

  call gf_free(gfldo)

  500 continue

! end of bias estimation between CMC and NCEP analysis

enddo

call baclose(ifile,iret)
call baclose(afile,iret)
call baclose(afile_m06,iret)
call baclose(rfile,iret)
call baclose(rfile_m06,iret)
call baclose(ofile,iret)

print *,'Bias Estimation Between CMC and NCEP Analysis Successfully Complete'

endif
stop

1020  continue

stop
end

subroutine decay(aveeror,fgrid,agrid,maxgrd,dec_w)

! apply the decaying average scheme
!
!     input      
!            fgrid  ---> ensemble forecast
!            agrid  ---> analysis data
!            aveeror---> bias estimation
!            dec_w  ---> decay weight
!
!     output
!            aveeror---> updated bias estimation

implicit none

integer maxgrd,ij
real aveeror(maxgrd),fgrid(maxgrd),agrid(maxgrd)
real dec_w           

do ij=1,maxgrd
  if(fgrid(ij).gt.-99999.0.and.fgrid(ij).lt.999999.0.and.agrid(ij).gt.-99999.0.and.agrid(ij).lt.999999.0) then
      if(aveeror(ij).gt.-99999.0.and.aveeror(ij).lt.999999.0) then
        aveeror(ij)= (1-dec_w)*aveeror(ij)+dec_w*(fgrid(ij)-agrid(ij))
      else
        aveeror(ij)= dec_w*(fgrid(ij)-agrid(ij))
      endif
  else
    if(aveeror(ij).gt.-99999.0 .and.aveeror(ij).lt.999999.0) then
      aveeror(ij)= aveeror(ij)                   
    else
      aveeror(ij)= 0.0                                
    endif
  endif
enddo

return
end


subroutine cal_dewpt(dpt,tmp,rhp,maxgrd)

! calculate dew point temperature 
!
!    Compute Dew Point Temperature (Bolton 1980)
!    es = 6.112 * exp((17.67 * T)/(T + 243.5));
!    e = es * (RH/100.0);
!    Td = log(e/6.112)*243.5/(17.67-log(e/6.112));

!         where:
!           T = temperature in deg C;
!           es = saturation vapor pressure in mb;
!           e = vapor pressure in mb;
!           RH = Relative Humidity in percent;
!           Td = dew point in deg C 
!
!     parameters
!
!        input
!                  tmp  ---> temperature                   
!                  rhp  ---> relative humidity 
!
!        output
!                  dpt  ---> dew point temperature                   

implicit none

integer maxgrd,ij
real dpt(maxgrd),tmp(maxgrd),rhp(maxgrd)
real T,Td,es,e,RH

do ij=1,maxgrd
  T=tmp(ij)-273.15
  RH=rhp(ij)  
  if(RH.eq.0) then
    print *,'RH 0 value',ij,RH,e,Td,dpt(ij)
    if(ij.eq.1) then
      if(rhp(ij+1).ne.0.0) then
        RH=rhp(ij+1)
      else
        RH=0.5
      endif
    elseif(ij.eq.maxgrd) then
      if(rhp(ij-1).ne.0.0) then
        RH=rhp(ij-1)
      else
        RH=0.5
      endif
    else
      if(rhp(ij-1).ne.0.0.and.rhp(ij+1).ne.0.0) then
        RH=(rhp(ij+1)+rhp(ij-1))/2.0
      else
        RH=0.5
      endif
    endif
  endif
  es=6.112 * exp((17.67 * T)/(T + 243.5))
  e=es*(RH/100.0)
  Td=log(e/6.112)*243.5/(17.67-log(e/6.112))+273.15
  dpt(ij)=min(T+273.15,Td)                 
enddo
 
return 
end