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