diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 7b2eb2da..e80f514e 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -17,3 +17,8 @@ add_subdirectory(overgridid.fd) add_subdirectory(webtitle.fd) add_subdirectory(ocnicepost.fd) + +add_subdirectory(ensadd.fd) +add_subdirectory(ensppf.fd) +add_subdirectory(ensstat.fd) +add_subdirectory(wave_stat.fd) diff --git a/src/ensadd.fd/CMakeLists.txt b/src/ensadd.fd/CMakeLists.txt new file mode 100644 index 00000000..02bcd2f2 --- /dev/null +++ b/src/ensadd.fd/CMakeLists.txt @@ -0,0 +1,12 @@ +list(APPEND fortran_src + ENSADD.f90 + printinfr.f90 +) + +set(exe_name ensadd.x) +add_executable(${exe_name} ${fortran_src}) +target_link_libraries(${exe_name} PRIVATE bacio::bacio_4 + w3emc::w3emc_4 + g2::g2_d) + +install(TARGETS ${exe_name} RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}) diff --git a/src/ensadd.fd/ENSADD.f90 b/src/ensadd.fd/ENSADD.f90 new file mode 100644 index 00000000..d7fddf6d --- /dev/null +++ b/src/ensadd.fd/ENSADD.f90 @@ -0,0 +1,187 @@ +!$$$ MAIN PROGRAM DOCUMENTATION BLOCK +! +! MAIN PROGRAM: ENSADD_G2 +! PRGMMR: ZHU ORG: NP23 DATE: 1999-08-31 +! +! ABSTRACT: THIS PROGRAM WILL EXTEND PDS MESSAGE WHICH WILL +! INCLUDE ENS(5) MESSAGE +! +! PROGRAM HISTORY LOG: +! 96-10-?? MARK IREDELL - Originator +! 97-03-17 YUEJIAN ZHU - Added DOCBLOACK +! 99-07-26 YUEJIAN ZHU - Modified to IBM-SP +! 14-10-06 BO Cui - Modified to Decode/Incode GRIB2 Data +! +! USAGE: +! +! INPUT FILES: +! UNIT 11 GRIB FILE +! UNIT 5 READ *, IENST,IENSI (For ensemble message) +! +! OUTPUT FILES: +! UNIT 51 GRIB FILE +! +! SUBPROGRAMS CALLED: +! GETGB2 -- W3LIB ROUTINE +! PUTGB2 -- W3LIB ROUTINE +! +! ATTRIBUTES: +! LANGUAGE: FORTRAN +! +!$$$ + +program ensadd_g2 + +use grib_mod +use params + +implicit none + +type(gribfield) :: gfldo + +integer :: currlen=0 +logical :: unpack=.true. +logical :: expand=.false. + +integer,dimension(200) :: kids,kpdt,kgdt +integer kskp,kdisc,kpdtn,kgdtn,i + +integer,dimension(200) :: jids,jpdt,jgdt,iids,ipdt,igdt +integer jskp,jdisc,jpdtn,jgdtn,idisc,ipdtn,igdtn + +integer temp(200) + +integer ienst,iensi,lpgb,lpgi,lpge,icount,iretb,ireti,irete,iret,jret,ipdtnum_out +character*255 cpgb,cpge +namelist /namin/ ienst,iensi,cpgb,cpge + +CALL W3TAGB('ENSADD',1999,0243,0068,'NP23') + +read (5,namin) +lpgb=len_trim(cpgb) +lpge=len_trim(cpge) + +print *, cpgb(1:lpgb),' ',cpge(1:lpge) +print *, ' ' + +call baopenr(11,cpgb(1:lpgb),iretb) +call baopen (51,cpge(1:lpge),irete) + +! loop over variables + +kids=-9999;kpdt=-9999; kgdt=-9999 +kdisc=-1; kpdtn=-1; kgdtn=-1 + +icount=0 +kskp=0 +do + + print *, '----- Read Ensemble Forecast ------' + print *, ' ' + + call getgb2(11,0,kskp,kdisc,kids,kpdtn,kpdt,kgdtn,kgdt,unpack,kskp,gfldo,iret) + + if(iret.ne.0) then + if(iret.eq.99 ) exit + print *,' getgb2 error = ',iret + cycle + !call errexit(17) + endif + + icount=icount+1 + + if(iret.eq.0) then + call printinfr(gfldo,icount) + else + print*, 'there is no fcst for ens member' + endif + + ! gfldo%idsect(2): Identification of originating + ! gfldo%idsect(2)=2: NCEP Ensemble Products + ! Original GFS has gfldo%idsect(2)=0, change it to be 2 + + gfldo%idsect(2)=2 + gfldo%idsect(13)=3 + + ! when product difinition template 4.0 change to 4.1 + ! when product difinition template 4.8 change to 4.11 + ! ipdtlen aslo change, need do modification for output + +! print*, '1 gfldo%ipdtmpl=',gfldo%ipdtmpl + + if(gfldo%ipdtnum.eq.0) then + + temp=-9999 + + temp(1:gfldo%ipdtlen)=gfldo%ipdtmpl(1:gfldo%ipdtlen) + if(gfldo%ipdtmpl(3).eq.2) then + temp(3)=4 + endif + + deallocate (gfldo%ipdtmpl) + + gfldo%ipdtnum=1 + if(gfldo%ipdtnum.eq.1) gfldo%ipdtlen=18 + if(gfldo%ipdtnum.eq.1) allocate (gfldo%ipdtmpl(gfldo%ipdtlen)) + + gfldo%ipdtmpl(1:15)=temp(1:15) + gfldo%ipdtmpl(16)=ienst + gfldo%ipdtmpl(17)=iensi + gfldo%ipdtmpl(18)=10 + + endif + + if(gfldo%ipdtnum.eq.8) then + + temp=-9999 + + temp(1:gfldo%ipdtlen)=gfldo%ipdtmpl(1:gfldo%ipdtlen) + if(gfldo%ipdtmpl(3).eq.2) then + temp(3)=4 + endif + + deallocate (gfldo%ipdtmpl) + + gfldo%ipdtnum=11 + if(gfldo%ipdtnum.eq.11) gfldo%ipdtlen=32 + if(gfldo%ipdtnum.eq.11) allocate (gfldo%ipdtmpl(gfldo%ipdtlen)) + + gfldo%ipdtmpl(1:15)=temp(1:15) + gfldo%ipdtmpl(16)=ienst + gfldo%ipdtmpl(17)=iensi + gfldo%ipdtmpl(18)=10 + gfldo%ipdtmpl(19:32)=temp(16:29) + + endif + +! print*, '2 gfldo%ipdtmpl=',gfldo%ipdtmpl + +! print*, 'gfldo%ipdtnum=',gfldo%ipdtnum +! print*, 'gfldo%ipdtlen=',gfldo%ipdtlen + + print *, '----- Write Ensemble Forecast ------' + print *, ' ' + + call printinfr(gfldo,icount) + call putgb2(51,gfldo,jret) + + ! end of probability forecast calculation + + 200 continue + + call gf_free(gfldo) + +enddo + +! end of ivar loop + +! close files + +call baclose(11,iretb) +call baclose(51,irete) + +CALL W3TAGE('ENSADD') + +stop +end + diff --git a/src/ensadd.fd/printinfr.f90 b/src/ensadd.fd/printinfr.f90 new file mode 100644 index 00000000..3727b52b --- /dev/null +++ b/src/ensadd.fd/printinfr.f90 @@ -0,0 +1,82 @@ +subroutine printinfr(gfld,ivar) + +! SUBPROGRAM: printinfr +! +! PRGMMR: Bo Cui DATE: 2013-06-11 +! +! USAGE: print grib2 data information +! +! INPUT: gfld,ivar + +use grib_mod +use params + +implicit none + +type(gribfield) :: gfld + +integer,dimension(200) :: jids,jpdt,jgdt,ipdt,igdt +integer :: currlen=0 +logical :: unpack=.true. + +integer kf,j,ivar,i +real fldmin,fldmax + +!kf=gfld%ngrdpts +kf=gfld%ndpts + +fldmin=gfld%fld(1) +fldmax=gfld%fld(1) + +do j=2,kf + if (gfld%fld(j).gt.fldmax) fldmax=gfld%fld(j) + if (gfld%fld(j).lt.fldmin) fldmin=gfld%fld(j) +enddo + +! print out + +! gfld%ipdtnum 1/11: ens. fcst or control, high reslution +! gfld%ipdtnum 2/12: ens. average fcst +! gfld%ipdtnum 0/8: cdas reanalysis + +! gfld%ipdtnum 11: ens. fcst or control in a continuous or non-continuous time interval +! gfld%ipdtnum 12: derived fcst based on ens. members in a continuous/non-continuous time interval +! gfld%ipdtnum 8: Rstatistically processed values in a continuous/non-continuous time interval + +if(gfld%ipdtnum.eq.11.or.gfld%ipdtnum.eq.12) then + write(6,100) + write(6,102) ivar,gfld%ipdtnum,(gfld%ipdtmpl(i),i=1,3),gfld%ipdtmpl(10),gfld%ipdtmpl(12), & + (gfld%idsect(i),i=6,9),gfld%ipdtmpl(9),gfld%ipdtmpl(30), & + (gfld%ipdtmpl(i),i=16,17), & + kf,fldmax,fldmin,gfld%fld(8601) + write(6,*) + +elseif(gfld%ipdtnum.eq.8.or.gfld%ipdtnum.eq.0) then + write(6,300) + write(6,302) ivar,gfld%ipdtnum,(gfld%ipdtmpl(i),i=1,3),gfld%ipdtmpl(10),gfld%ipdtmpl(12), & + (gfld%idsect(i),i=6,9),gfld%ipdtmpl(9), & + kf,fldmax,fldmin,gfld%fld(8601) + write(6,*) + +else + write(6,200) + write(6,202) ivar,gfld%ipdtnum,(gfld%ipdtmpl(i),i=1,3),gfld%ipdtmpl(10),gfld%ipdtmpl(12), & + (gfld%idsect(i),i=6,9),gfld%ipdtmpl(9),(gfld%ipdtmpl(i),i=16,17), & + kf,fldmax,fldmin,gfld%fld(8601) + write(6,*) +endif + +100 format(' REC PDTN PD1 PD2 PD3 PD10 PD12 YEAR MN DY HR FHR TR ', & + 'E16 E17 LEN MAX MIN EXAMPLE') +102 format(i4,i5,4i4,i8,i6,3i3,2i4,2i4,i8,3f11.2) + +200 format(' REC PDTN PD1 PD2 PD3 PD10 PD12 YEAR MN DY HR FHR ', & + 'E16 E17 LEN MAX MIN EXAMPLE') +202 format(i4,i5,4i4,i8,i6,3i3,3i4,i8,3f11.2) + +300 format(' REC PDTN PD1 PD2 PD3 PD10 PD12 YEAR MN DY HR FHR ', & + ' LEN MAX MIN EXAMPLE') +302 format(i4,i5,4i4,i8,i6,3i3,i4,8x,i8,3f11.2) + +return +end diff --git a/src/ensppf.fd/CMakeLists.txt b/src/ensppf.fd/CMakeLists.txt new file mode 100644 index 00000000..d65d6730 --- /dev/null +++ b/src/ensppf.fd/CMakeLists.txt @@ -0,0 +1,13 @@ +list(APPEND fortran_src + ENSPPF.f90 + printinfr.f90 + init_parm.f90 +) + +set(exe_name ensppf.x) +add_executable(${exe_name} ${fortran_src}) +target_link_libraries(${exe_name} PRIVATE bacio::bacio_4 + w3emc::w3emc_4 + g2::g2_d) + +install(TARGETS ${exe_name} RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}) diff --git a/src/ensppf.fd/ENSPPF.f90 b/src/ensppf.fd/ENSPPF.f90 new file mode 100644 index 00000000..4e295701 --- /dev/null +++ b/src/ensppf.fd/ENSPPF.f90 @@ -0,0 +1,350 @@ +!$$$ MAIN PROGRAM DOCUMENTATION BLOCK +! +! MAIN PROGRAM: ENSPPF ENNSEMBLE PRECIP. PROBABILITY +! PRGMMR: YUEJIAN ZHU ORG:NP23 DATE: 97-03-17 +! +! ABSTRACT: THIS PROGRAM WILL CALCULATE ENSEMBLE BASED PRECIP. +! PROBABILITY FORECAST (PPF) +! +! PROGRAM HISTORY LOG: +! 97-03-17 YUEJIAN ZHU (WD20YZ) +! 99-07-26 YUEJIAN ZHU (WX20YZ) MODITY TO IBM-SP +! 00-04-17 YUEJIAN ZHU (WX20YZ) INCREASE MEMBERS FROM 17 TO 23 +! 03-09-09 YUEJIAN ZHU (WX20YZ) REDUCE MEMBERS FROM 23 TO 11 +! TO MAKE COMPARRABLE RESOLUTIONS AND SKILLS +! 04-02-10 YUEJIAN ZHU (WX20YZ) RESOLVE A BUG WHICH PRODUCE A +! INCORRECT DATA MESSAGE +! 06-01-25 YUEJIAN ZHU (WX20YZ) CHANGE FOR NEW IMPLEMANTED ENSEMBLE +! PRODUCTION (14+1+1 MEMBERS AND ONE SIDE ENSMBLE) +! 14-11-06 BO Cui: Modify for grib2 encode/decode +! Change for full ensemble member, +! production(20+1 members and one side ensemble) +! 18-05-10 Xianwu Xue : Add 'npert' to calculate 'mem' +! +! +! USAGE: +! +! INPUT FILES: +! UNIT 11 PRECIPITATION GRIB FILE ( 144*73 ) +! +! OUTPUT FILES: +! UNIT 51 PPF GRIB FILE ( 144*73 ) +! +! SUBPROGRAMS CALLED: +! GETGB2 -- W3LIB ROUTINE +! PUTGB2 -- W3LIB ROUTINE +! +! ATTRIBUTES: +! LANGUAGE: FORTRAN +! +program ppf + +use grib_mod +use params + +implicit none + +type(gribfield) :: gfld,gfldo + +integer :: currlen=0 +logical :: unpack=.true. +logical :: expand=.false. + +integer,dimension(200) :: kids,kpdt,kgdt +integer kskp,kdisc,kpdtn,kgdtn,i + +integer,dimension(200) :: jids,jpdt,jgdt,iids,ipdt,igdt +integer jskp,jdisc,jpdtn,jgdtn,idisc,ipdtn,igdtn +common /param/jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt + +integer len,mem,icyc,ijd,ncnt,n,l,icnt,m,ii,k +integer lpgb,lpge,iretb,irete,imem,iret +integer lpgi,ireti,jj +integer temp(200),ipdt9,P2,iunit2 +real bb,cc + +parameter(len=60) + +integer :: npert = 20 + +real, allocatable :: ff(:,:),pp(:,:),ss(:,:),aa(:),gg(:,:),hh(:,:) + +character*255 cpgb,cpgi,cpge + +namelist /namin/ cpgb,cpge,npert + +real rk(9) +data rk/0.254,1.00,2.54,5.00,6.35,10.00,12.7,25.4,50.8/ + +CALL W3TAGB('ENSPPF',2000,0110,0073,'NP20 ') + +read (5,namin,end=1020) +!write (6, namin) +!print *, npert + +mem = npert + 1 + +lpgb=len_trim(cpgb) +lpge=len_trim(cpge) +print *, cpgb(1:lpgb) +print *, cpge(1:lpge) +call baopenr(11,cpgb(1:lpgb),iretb) +call baopenw(51,cpge(1:lpge),irete) + +print*,' ' + +! find grib message + +iids=-9999;ipdt=-9999; igdt=-9999 +idisc=-1; ipdtn=-1; igdtn=-1 +call init_parm(ipdtn,ipdt,igdtn,igdt,idisc,iids) +call getgb2(11,0,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt,unpack,jskp,gfld,iret) + +if(iret.eq.0) then + + ijd=gfld%ngrdpts + + ! check Indicator of unit of time ranga + ! gfld%ipdtmpl(8): indicator of unit of time range + ! 11: 6 hours; 1: hour; 10: 3 hour; 12: 12 hours + + if(gfld%ipdtmpl(8).eq.1) then + print *, ' Unit of time range is 1 hour '; print *, ' ' + elseif(gfld%ipdtmpl(8).eq.11) then + print *, ' Unit of time range is 6 hours '; print *, ' ' + else + print *, 'Unit of time range is not 1 hour or 6 hours '; print *, ' ' + endif + + ! check gfld%ipdtmpl(30) for PDT number 4.11 + ! length of the time range in units defined by the previous octet + + if(gfld%ipdtmpl(30).eq.6) then + print *, ' Length of time range is 6 hours '; print *, ' ' + else + print *, ' Length of time range is not 6 hours '; print *, ' ' + endif + + call gf_free(gfld) + + allocate(ff(ijd,mem),pp(ijd,mem),ss(ijd,mem),aa(ijd)) + allocate(gg(ijd,mem),hh(ijd,mem)) + + ncnt=0 + + do n = 1, len !### len=steps + + ! Part I: get ctl + 20 ensemble members precipitation data + + icnt=0 + ff=0.0 + do m=1, mem + + iids=-9999;ipdt=-9999; igdt=-9999 + idisc=-1; ipdtn=-1; igdtn=-1 + + ! read and process input ensemble member + + ! read in control member for m=1 + + if(m.eq.1) then + ipdt(16)=1 ! type of ensemble forecast + ipdt(17)=0 ! perturbation number + else + ipdt(16)=3 ! type of ensemble forecast + ipdt(17)=m-1 ! perturbation number + endif + + ! ipdt(8)=1 ! time unit: 1 hour - kpds(13) in grib1 + + ! gfld%ipdtmpl(9): forecast time in units + + ipdt(9)=int((n-1)*6) + + ! gfld%ipdtmpl(29): indicator of unit of time ranga for process + ! 11: 6 hours; 1: hour; 10: 3 hour; 12: 12 hours + + ipdt(29)=1 ! time unit: 1 for 1 hour + + ! check gfld%ipdtmpl(30) for PDT number 4.11 + ! length of the time range in units defined by the previous octet + + ipdt(30)=6 ! length of the time range is 6 hours + + ipdtn=11; igdtn=-1 + + call init_parm(ipdtn,ipdt,igdtn,igdt,idisc,iids) + call getgb2(11,0,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt,unpack,jskp,gfld,iret) + + if(iret.ne.0) print *, 'there is no varibale for member ',m + + if(iret.eq.0) then + icnt=icnt + 1 + print *, '----- Input Data for Current Time ------' + call printinfr(gfld,m) + do ii=1,ijd + ff(ii,icnt)=gfld%fld(ii) + enddo + else + ncnt=ncnt+1 + if(ncnt.le.1) then + print *,' n=',n,' iret=',iret + endif + endif ! end of iret.eq.0 + + if(icnt.eq.1) then + if (n.ge.4.and.mod(n,2).eq.0) then + gfldo=gfld + else + call gf_free(gfld) + endif + else + call gf_free(gfld) + endif + + enddo !### for m = 1, mem + ! + ! + ! PART II: to calculate the probability scores + ! + ! skip n=1-3 and other odd steps, when you calculate 24-hour interval + + if (n.ge.4.and.mod(n,2).eq.0) then + + ! change grib2 pdt message for new ensemble products + + gfldo%idsect(2)=2 ! Identification of originating/generating subcenter + ! 2: NCEP Ensemble Products + + gfldo%idsect(13)=5 ! Type of processed data in this GRIB message + ! 5: Control and Perturbed Forecast Products + + temp=-9999 + +! print *, 'gfldo%ipdtlen=',gfldo%ipdtlen +! print *, 'gfldo%ipdtmpl=',gfldo%ipdtmpl + + temp(1:gfldo%ipdtlen)=gfldo%ipdtmpl(1:gfldo%ipdtlen) + + temp(3)=5 ! Type of generating process + ! 5: Probability Forecast + + deallocate (gfldo%ipdtmpl) + + gfldo%ipdtnum=9 ! Probability forecasts from ensemble + if(gfldo%ipdtnum.eq.9) gfldo%ipdtlen=36 + if(gfldo%ipdtnum.eq.9) allocate (gfldo%ipdtmpl(gfldo%ipdtlen)) + + gfldo%ipdtmpl(1:15)=temp(1:15) + + gfldo%ipdtmpl(1)=1 ! Parameter category : 1 Moisture + gfldo%ipdtmpl(2)=8 ! Parameter number : 8 Total Precipitation(APCP) + +! unit of time range(time interval) is 6 hours for gefs +! for ipdt 4.11, unit of time range: gfldo%ipdtmpl(29) +! forecast time p1: gfldo%ipdtmpl(9) +! this is for ipdt 4.11 only + + if(temp(8).eq.1) then + if(temp(29).eq.1) then + iunit2=1 + P2=abs(temp(30)*iunit2) + ipdt9=int(temp(9)+P2) + gfldo%ipdtmpl(9)=int(ipdt9-24) ! Forecast time in units + endif + endif + + gfldo%ipdtmpl(16)=0 ! Forecast probability number + gfldo%ipdtmpl(17)=mem ! Total number of forecast probabilities + gfldo%ipdtmpl(18)=1 ! Probability Type + ! 1: Probability of event above upper limit + + gfldo%ipdtmpl(19)=0 ! Scale factor of lower limit + gfldo%ipdtmpl(20)=0 ! Scaled value of lower limit + gfldo%ipdtmpl(21)=3 ! Scale factor of upper limit + + ! gfldo%ipdtmpl(22) will be set below + + gfldo%ipdtmpl(23:36)=temp(19:32) + + ! gfldo%ipdtmpl(34): Length of the time range over which processing done + + if(gfldo%ipdtmpl(8).eq.1) then + gfldo%ipdtmpl(34)=24 ! Length of the time range + endif + + ! start to calculate the probability scores + + do k = 1, 9 + aa=0.0 + do ii = 1, ijd + do m = 1, icnt + bb=(ff(ii,m)+pp(ii,m)+gg(ii,m)+hh(ii,m)) + if (bb.ge.rk(k)) then + aa(ii) = aa(ii) + 1.0 + endif + enddo + enddo + do ii = 1, ijd + aa(ii) = aa(ii)*100.0/float(icnt) + if (aa(ii).ge.99.0) then + aa(ii) = 100.0 + endif + enddo + + ! + ! testing print + ! + ! if (n.eq.2.and.k.eq.2) then + ! write(*,999) (aa(ii),ii=1001,1100) + ! endif + !999 format (10f8.1) + + print *, '----- Output for Current Time ------' + + ! gfldo%ipdtmpl(22): Scaled value of upper limit + + gfldo%ipdtmpl(22)=rk(k)*(10**gfldo%ipdtmpl(21)) + + gfldo%fld(1:ijd)=aa(1:ijd) + + call printinfr(gfldo,k) + call putgb2(51,gfldo,iret) + + enddo !### for k = 1, 9 + + call gf_free(gfldo) + + endif + + if (icnt.gt.0) then + do ii = 1, ijd + do jj = 1, mem + hh(ii,jj)=gg(ii,jj) + gg(ii,jj)=pp(ii,jj) + pp(ii,jj)=ff(ii,jj) + enddo + enddo + endif + + enddo !### for n = 1, len + + call baclose(11,iretb) + call baclose(51,irete) + + CALL W3TAGE('ENSPPF') + + stop +else + print*,' getgbeh, cannot get maxgrd ' +endif + +1020 Continue + +print *,'Wrong Data Input, Output or Wrong Message Input' + +stop + +end + + diff --git a/src/ensppf.fd/init_parm.f90 b/src/ensppf.fd/init_parm.f90 new file mode 100644 index 00000000..69cc2bdb --- /dev/null +++ b/src/ensppf.fd/init_parm.f90 @@ -0,0 +1,57 @@ +subroutine init_parm(ipdtn,ipdt,igdtn,igdt,idisc,iids) + +! SUBPROGRAM: init_parm +! +! PRGMMR: Bo Cui DATE: 2013-12-01 +! +! ABSTRACT: This subroutine returns the Grid Definition, and +! Product Definition for a given data field. +! +! PROGRAM HISTORY LOG: +! 2013-12-01 Bo Cui +! +! USAGE: call init_parm(ipdtn,ipdt,igdtn,igdt) +! +! INPUT: ipdtn,ipdt,igdtn,igdt +! OUTPUT: jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt +! +! REMARKS: None +! +! ATTRIBUTES: +! LANGUAGE: Fortran 90 +! MACHINE: IBM SP +! +!$$$ + +implicit none + +integer,dimension(200) :: jids,jpdt,jgdt,iids,ipdt,igdt +integer jskp,jdisc,jpdtn,jgdtn,idisc,ipdtn,igdtn +common /param/jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt + +jskp=0 + +jdisc=-1 +jpdtn=ipdtn +jgdtn=igdtn + +jids=-9999 +jpdt=-9999 +jgdt=-9999 + +! input product defining values + +jids=iids +jpdt=ipdt +jgdt=igdt + +!jpdt(1)=ipdt(1) +!jpdt(2)=ipdt(2) +!jpdt(10)=ipdt(10) +!jpdt(11)=ipdt(11) +!jpdt(12)=ipdt(12) + +return +end + + diff --git a/src/ensppf.fd/printinfr.f90 b/src/ensppf.fd/printinfr.f90 new file mode 100644 index 00000000..877da8df --- /dev/null +++ b/src/ensppf.fd/printinfr.f90 @@ -0,0 +1,96 @@ +subroutine printinfr(gfld,ivar) + +! SUBPROGRAM: printinfr +! +! PRGMMR: Bo Cui DATE: 2013-06-11 +! +! PROGRAM HISTORY LOG: +! 14-11-04 Bo Cui: Add information printing for ipdt 4.9 +! +! USAGE: print grib2 data information +! +! INPUT: gfld,ivar + +use grib_mod +use params + +implicit none + +type(gribfield) :: gfld + +integer,dimension(200) :: jids,jpdt,jgdt,ipdt,igdt +integer :: currlen=0 +logical :: unpack=.true. + +integer kf,j,ivar,i +real fldmin,fldmax + +kf=gfld%ndpts + +fldmin=gfld%fld(1) +fldmax=gfld%fld(1) + +do j=2,kf + if (gfld%fld(j).gt.fldmax) fldmax=gfld%fld(j) + if (gfld%fld(j).lt.fldmin) fldmin=gfld%fld(j) +enddo + +! print out + +! gfld%ipdtnum 1/11: ens. fcst or control, high reslution +! gfld%ipdtnum 2/12: ens. average fcst +! gfld%ipdtnum 0/8: cdas reanalysis + +! gfld%ipdtnum 11: ens. fcst or control in a continuous or non-continuous time interval +! gfld%ipdtnum 12: derived fcst based on ens. members in a continuous/non-continuous time interval +! gfld%ipdtnum 8: Rstatistically processed values in a continuous/non-continuous time interval + +! gfld%ipdtnum 9: Probability forecasts in a continuous or non-continuous time interval + +if(gfld%ipdtnum.eq.11.or.gfld%ipdtnum.eq.12) then + write(6,100) + write(6,102) ivar,gfld%ipdtnum,(gfld%ipdtmpl(i),i=1,3),gfld%ipdtmpl(10),gfld%ipdtmpl(12), & + (gfld%idsect(i),i=6,9),gfld%ipdtmpl(9),gfld%ipdtmpl(30), & + (gfld%ipdtmpl(i),i=16,17), & + kf,fldmax,fldmin,gfld%fld(8601) + write(6,*) + +elseif(gfld%ipdtnum.eq.9) then + write(6,400) + write(6,102) ivar,gfld%ipdtnum,(gfld%ipdtmpl(i),i=1,3),gfld%ipdtmpl(10),gfld%ipdtmpl(12), & + (gfld%idsect(i),i=6,9),gfld%ipdtmpl(9),gfld%ipdtmpl(34), & + (gfld%ipdtmpl(i),i=17,18), & + kf,fldmax,fldmin,gfld%fld(8601) + write(6,*) + +elseif(gfld%ipdtnum.eq.8.or.gfld%ipdtnum.eq.0) then + write(6,300) + write(6,302) ivar,gfld%ipdtnum,(gfld%ipdtmpl(i),i=1,3),gfld%ipdtmpl(10),gfld%ipdtmpl(12), & + (gfld%idsect(i),i=6,9),gfld%ipdtmpl(9), & + kf,fldmax,fldmin,gfld%fld(8601) + write(6,*) + +else + write(6,200) + write(6,202) ivar,gfld%ipdtnum,(gfld%ipdtmpl(i),i=1,3),gfld%ipdtmpl(10),gfld%ipdtmpl(12), & + (gfld%idsect(i),i=6,9),gfld%ipdtmpl(9),(gfld%ipdtmpl(i),i=16,17), & + kf,fldmax,fldmin,gfld%fld(8601) + write(6,*) +endif + +100 format(' REC PDTN PD1 PD2 PD3 PD10 PD12 YEAR MN DY HR FHR TR ', & + 'E16 E17 LEN MAX MIN EXAMPLE') +102 format(i4,i5,4i4,i8,i6,3i3,2i4,2i4,i8,3f11.2) + +200 format(' REC PDTN PD1 PD2 PD3 PD10 PD12 YEAR MN DY HR FHR ', & + 'E16 E17 LEN MAX MIN EXAMPLE') +202 format(i4,i5,4i4,i8,i6,3i3,3i4,i8,3f11.2) + +300 format(' REC PDTN PD1 PD2 PD3 PD10 PD12 YEAR MN DY HR FHR ', & + ' LEN MAX MIN EXAMPLE') +400 format(' REC PDTN PD1 PD2 PD3 PD10 PD12 YEAR MN DY HR FHR TR ', & + 'E17 E18 LEN MAX MIN EXAMPLE') +302 format(i4,i5,4i4,i8,i6,3i3,i4,8x,i8,3f11.2) + +return +end diff --git a/src/ensstat.fd/CMakeLists.txt b/src/ensstat.fd/CMakeLists.txt new file mode 100644 index 00000000..aa815d72 --- /dev/null +++ b/src/ensstat.fd/CMakeLists.txt @@ -0,0 +1,17 @@ +list(APPEND fortran_src + ensstat.f90 + init_parm.f90 + printinfr.f90 + gtbits.f90 + isrchne.f90 + EPDF.f + SORTMM.f +) + +set(exe_name ensstat.x) +add_executable(${exe_name} ${fortran_src}) +target_link_libraries(${exe_name} PRIVATE bacio::bacio_4 + w3emc::w3emc_4 + g2::g2_d) + +install(TARGETS ${exe_name} RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}) diff --git a/src/ensstat.fd/EPDF.f b/src/ensstat.fd/EPDF.f new file mode 100644 index 00000000..ab0a1afa --- /dev/null +++ b/src/ensstat.fd/EPDF.f @@ -0,0 +1,365 @@ +c +c FUNCTION EPDF(EF,EW,N,Q,ICTL) +c Prgmmr: Yuejian Zhu Org: np23 Date: 2007-0801 +c +c function epdf to build up Probability Density Function by using picewise +c linear approximation. +c +c parameters: +c input: +c ef -- forecasts for n ensemble members +c ew -- forecast weights for n ensemble members +c n -- ensemble members +c qq -- forecast value (specified) when ictl = 1 +c -- probability value (0-1.0, specified) when ictl = -1 +c -- 1.0 for weighted mean when ictl = 0 +c -- -1.0 for unweighted mean when ictl = 0 +c -- 2.0 for weighted spread when ictl = 0 +c -- -2.0 for unweighted spread when ictl = 0 +c ictl -- 1 to calculate probability from forecast value +c -1 to calculate forecast value from giving probability +c 0 to calculate special requests +c +c Fortran 77 on IBMSP +c + function epdf(ef,ew,n,qq,ictl) + dimension ef(n),ew(n) + dimension q(n),u(n),w(n),qa(n+1),d(n,3) +C--------+---------+---------+---------+---------+---------+---------+---------+ +c +c normalized the weights in case of input weights are not normalized +c +c print *, "ew=",(ew(i),i=1,n) +c print *, "ef=",(ef(i),i=1,n) + ew1 = 0.0 + do i = 1, n + ew1 = ew1 + ew(i) + enddo + do i = 1, n + ew(i) = ew(i)/ew1 + enddo +c + if (ictl.eq.0) then +c +c for special requests: +c if qq=1.0, to calculate weighted mean +c if qq=-1.0, to calculate unweighted mean +c if qq=2.0, to calculate weighted spread +c if qq=-2.0, to calculate unweighted spread +c + wm = 0.0 + uwm = 0.0 + do i = 1, n + wm = wm + ef(i)*ew(i) + uwm = uwm + ef(i)/float(n) + enddo + epdf = 0.0 + if (qq.eq.1.0) then + epdf = wm + else if (qq.eq.-1.0) then + epdf = uwm + else if (qq.eq.2.0) then + ssprd = 0.0 + do i = 1, n + ssprd = ssprd + (ef(i)-wm)*ew(i)*(ef(i)-wm)*ew(i) + enddo + epdf = sqrt(ssprd*float(n*n)/float(n-1)) + else if (qq.eq.-2.0) then + ssprd = 0.0 + do i = 1, n + ssprd = ssprd + (ef(i)-uwm)*(ef(i)-uwm) + enddo + epdf = sqrt(ssprd/float(n-1)) +c return + else + epdf=-9999.99 +c return + endif + return + endif +c +c This sorting program will: +c keep original data order +c get new order (low to high) +c new order data with original index +c + do i = 1, n + d(i,1) = ef(i) + enddo + call sortmm(d,n,3,1) +c write (*,'("ORG DATA SET",13f5.1)') (d(i,1),i=1,n) +c write (*,'("NEW DATA SET",13f5.1)') (d(i,2),i=1,n) +c write (*,'("N DATA INDEX",13f5.1)') (d(i,3),i=1,n) + + do i = 1, n + ef(i) = d(i,2) + enddo +c +c consider statistical the same values if their difference is less wqual to +c + cr = (ef(n) - ef(1))*0.001 +c +c to eliminate the duplication, give extra weight +c m will be new dimension after this process (m<=n) +c + m = 1 + q(1) = ef(1) + u(1) = ew(int(d(1,3))) + do i = 2, n + if ((ef(i)-ef(i-1)).gt.cr) then + m = m + 1 + q(m) = ef(i) + u(m) = ew(int(d(i,3))) + else + u(m) = u(m) + ew(int(d(i,3))) + endif + enddo +c +c for extreme case, all ensemble values are the same +c + if (m.eq.1) then + if (ictl.eq.1) then + if (qq.lt.q(1)) then + epdf=0.0 + return + else if (qq.gt.q(m)) then + epdf=100.0 + return + else + epdf=50.0 + return + endif + else if (ictl.eq.-1) then + epdf=q(1) + return + else if (ictl.eq.0) then + if (abs(qq).eq.1) then + epdf=q(1) + return + else if (abs(qq).eq.2) then + epdf=0.0 + return + else + epdf=-9999.99 + return + endif + else + epdf=-9999.99 + return + endif + endif + + do i = 1, n + ef(i) = d(i,1) + enddo + +c write (*,'("ENS WEIGHTS:",13f5.3)') (u(i),i=1,m) +c +c to calculate ensemble segment weights +c + w(1) = u(1) / (q(2) - q(1)) + w(m) = u(m) / (q(m) - q(m-1)) + do i = 2, m-1 +ccc find a bug to to calculate segment weight (08/07/2007) +c w(i) = (u(i+1) + u(i)) / (q(i+1) - q(i-1)) +ccc Yuejian suggested weight +c w(i) = (u(i+1) + 2.0*u(i) + u(i-1)) / (q(i+1) - q(i-1)) / 2.0 +ccc Keith suggested weight + w(i) = 2.0*u(i) / (q(i+1) - q(i-1)) + enddo + wsum=0.0 + do i = 1, m + wsum = wsum + w(i) + enddo +c +c normalized PDF weights for each value +c +c do i = 1, m +c w(i) = w(i)/wsum +c enddo +c write (*,'("PDF WEIGHTS:",13f5.3)') (w(i),i=1,m) +c +c calculate left and right bounds by approximation (considering tails) +c +c modfied qlt and qrt approximation for extreme cases (small spread) +c 08/17/2007 -Yuejian Zhu + if (3*m.le.2*n) then + qdelta=(q(m)-q(1))/float(n-1) + qltd=2.0/(w(1)*float(n+1)) + qrtd=2.0/(w(m)*float(n+1)) + if (qltd.gt.2.0*qdelta) then + qltd=2.0*qdelta + endif + if (qrtd.gt.2.0*qdelta) then + qrtd=2.0*qdelta + endif + qlt=q(1) - qltd + qrt=q(m) + qrtd + else +c +c find bug to calculate qlt and qrt +c 08/22/2007 -Yuejian Zhu + qlt=q(1) - 1.0/(w(1)*float(n+1)) + qrt=q(m) + 1.0/(w(m)*float(n+1)) +c qlt=q(1) - 2.0/(w(1)*float(n+1)) +c qrt=q(m) + 2.0/(w(m)*float(n+1)) + endif +c +c normalized area of each trapezoidal +c + qa(1) = (q(1) - qlt)*w(1)/2.0 + qun = qa(1) + do i = 2, m + qa(i) = (w(i-1)+w(i))*(q(i)-q(i-1))/2.0 + qun = qun + qa(i) + enddo + qa(m+1) = (qrt - q(m))*w(m)/2.0 + qun = qun + qa(m+1) + qa(1) = qa(1)/qun + do i = 2, m+1 + qa(i) = qa(i-1) + qa(i)/qun + enddo + +c write (*,'("CDF AREAS :",14f5.3)') (qa(i),i=1,m+1) +c print *, "left bound is ",qlt,": right bound is ",qrt + + if (ictl.eq.1) then +c +c to find out the position of give value qq +c + if (qq.le.qlt) then + epdf=0.0 + return + endif + if (qq.ge.qrt) then + epdf=1.0 + return + endif + + k=m + do i = 1, m + if (qq.lt.q(i)) then + k=i-1 + exit + endif + enddo +c 101 continue + + if (k.eq.0) then + ww = w(1)*(qq - qlt)/(q(1)-qlt) + fta = ww*(qq - qlt)/(2.0*qun) + epdf = fta + return + else if (k.eq.m) then +c +c cumulative trapezoid area (CTA) +c + cta = qa(m) +c +c fractional trapezoid area +c + ww = w(m)*(qrt - qq)/(qrt-q(m)) + fta = ww*(qrt - qq)/(2.0*qun) +c +c the area is the probability (0.0-1.0) of given value qq +c + epdf = cta + fta + return + else +c +c cumulative trapezoid area (CTA) +c + cta = qa(k) +c +c fractional trapezoid area +c + ww = w(k) + (w(k+1)-w(k))*(qq - q(k))/(q(k+1) - q(k)) + fta = (ww + w(k))*(qq - q(k))/(2.0*qun) +c +c the area is the probability (0.0-1.0) of given value qq +c + epdf = cta + fta + return + + endif + + elseif (ictl.eq.-1) then + + if (qq.le.0.0) then + epdf=qlt + return + endif + if (qq.ge.1.0) then + epdf=qrt + return + endif + + k = m + do i = 1, m + if (qq.lt.qa(i)) then + k=i-1 + exit + endif + enddo +c 102 continue + +c +c to solve the quadratic equation ( ax2 + bx + c = 0 ) +c + if (k.eq.0) then + fta = qq*qun + a = w(1)/(q(1)-qlt) + b = 0.0 + c = -2.0*fta + if (a.eq.0.0) then + epdf=-9999.99 + return + else + epdf = qlt + sqrt(-c/a) + return + endif + else if (k.eq.m) then + fta = (qq - qa(m))*qun + a = -w(m)/(qrt-q(m)) + b = 2.0*w(m) + c = -2.0*fta + b2m4ac = b*b - 4.0*a*c + if ( b2m4ac.lt.0.0) then + print *, "There is no real solution for this problem" + epdf=-9999.99 + return + endif + if (a.eq.0.0) then + epdf = q(m) - c/b + return + else + epdf = q(m) + (sqrt(b2m4ac) - b)/(2.0*a) + return + endif + else + fta = (qq - qa(k))*qun + a = (w(k+1)-w(k))/(q(k+1)-q(k)) + b = 2.0*w(k) + c = -2.0*fta + b2m4ac = b*b - 4.0*a*c + if ( b2m4ac.lt.0.0) then + print *, "There is no real solution for this problem" + epdf=-9999.99 + return + endif + if (a.eq.0.0) then + epdf = q(k) - c/b + return + else + epdf = q(k) + (sqrt(b2m4ac) - b)/(2.0*a) + return + endif + endif + else + epdf=-9999.99 + return + endif + + return + end + diff --git a/src/ensstat.fd/SORTMM.f b/src/ensstat.fd/SORTMM.f new file mode 100644 index 00000000..0086ed19 --- /dev/null +++ b/src/ensstat.fd/SORTMM.f @@ -0,0 +1,51 @@ + subroutine sortmm(a,n,nc,k) +CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC +C C +C USAGE: SORT ONE DIMENSION DATA WITH LENGTH N C +C CODE : F77 on IBMSP --- Yuejian Zhu (07/08/99) C +C C +C INPUT: array a(n,nc) C +C n -- input data length C +C nc -- variable second dimension (using nc=3) C +C k -- input data location a(*,k) C +C C +C OUTPUT: re-order array a(n,nc) C +C a(n,1) -- original data (if k=1) C +C a(n,2) -- new output after sorting (low -> high) C +C a(n,3) -- new output at original order index C +C C +CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC +C +c--------+----------+----------+----------+----------+----------+----------+-- + dimension a(n,nc),b(n,nc),js(n) + do i1 = 1, n + iless = 0 + imore = 0 + ieq = 0 + aa=a(i1,k) + do i2 = 1, n + bb=a(i2,k) + if ( aa.lt.bb ) iless = iless + 1 + if ( aa.gt.bb ) imore = imore + 1 + if ( aa.eq.bb ) then + ieq = ieq + 1 + js(ieq) = i2 + endif + enddo + if ( ieq.eq.1) then + b(imore+1,2)=aa + b(imore+1,1)=i1 + else + do i3 = 1, ieq + b(imore+i3,2)=aa + b(imore+i3,1)=js(i3) + enddo + endif + enddo + do jj= 1, n + a(jj,3) = b(jj,1) + a(jj,2) = b(jj,2) + enddo + return + end + diff --git a/src/ensstat.fd/ensstat.f90 b/src/ensstat.fd/ensstat.f90 new file mode 100644 index 00000000..77815920 --- /dev/null +++ b/src/ensstat.fd/ensstat.f90 @@ -0,0 +1,391 @@ +program ens_avgspr_g2 +! +! main program: ens_avgspr_g2 +! +! prgmmr: Bo Cui org: np/wx20 date: 2013-10-01 +! +! abstract: calculate ensemble mean & spread +! +! usage: +! +! input file: ncep/cmc/fnmoc ensemble forecast +! +! output file: ensemble mean and spread + +! programs called: +! baopenr grib i/o +! baopenw grib i/o +! baclose grib i/o +! getgb2 grib reader +! putgb2 grib writer +! init_parm define grid definition and product definition +! printinfr print grib2 data information +! change_template4 change data values for specified Product Definition Template + +! exit states: +! cond = 0 - successful run +! cond = 1 - I/O abort +! +! attributes: +! language: fortran 90 +! modified by: +! Xianwu Xue 05/11/2018 +! added 'navg_min' from namelist to determine the minimum members +!$$$ + +use grib_mod +use params + +!implicit none + +type(gribfield) :: gfld,gfldo +integer :: currlen=0 +logical :: unpack=.true. +logical :: expand=.false. + +integer,dimension(200) :: kids,kpdt,kgdt,temp +integer kskp,kdisc,kpdtn,kgdtn +integer firstfile + +integer,dimension(200) :: jids,jpdt,jgdt,iids,ipdt,igdt +integer jskp,jdisc,jpdtn,jgdtn,idisc,ipdtn,igdtn +common /param/jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt + +integer nmemd,nmvar,nvar,ivar,im,imem,n,inum +parameter (nmemd=62,nmvar=50) + +real, allocatable :: fgrid(:,:),fst(:),ens_avg(:),ens_spr(:) +real weight(nmemd) + +integer maxgrd,iret,jret,icount,i,ipdtnum_out + +integer ipd1,ipd2,ipd3,ipd10,ipd11,ipd12,ipdn + +integer iunit,lfipg(nmemd),icfipg(nmemd) +integer nfiles,nenspost,iskip(nmemd),tfiles,ifile +integer lfopg1,lfopg2 +integer icfopg1,icfopg2 + +character*255 cfipg(nmemd) +character*255 cfopg1,cfopg2 + +real gmin,gmax +integer nbit + +integer :: navg_min = 10 + +namelist /namens/nfiles,nenspost,cfipg,iskip,cfopg1,cfopg2,navg_min + +read (5,namens) +!write (6,namens) +!print *, navg_min + +print *, 'Input variables include ' + +! stop this program if there is no enough files put in + +print *, ' '; print *, 'Input files size ', nfiles + +if(nfiles.gt.2) then + + ! set the fort.* of intput file, open forecast files + + print *, ' ' + print *, 'Input files include ' + + iunit=9 + + tfiles=nfiles + + do ifile=1,nfiles + iunit=iunit+1 + icfipg(ifile)=iunit + lfipg(ifile)=len_trim(cfipg(ifile)) + print '(a4,i3,a98)', 'fort.',icfipg(ifile), cfipg(ifile)(1:lfipg(ifile)) + call baopenr(icfipg(ifile),cfipg(ifile)(1:lfipg(ifile)),iret) + if ( iret .ne. 0 ) then + print *,'there is no NAEFS forecast, ifile,iret = ',cfipg(ifile)(1:lfipg(ifile)),iret + tfiles=nfiles-1 + iskip(ifile)=1 + endif + enddo + + ! set the fort.* of output file + + print *, ' ' + print *, 'Output files include ' + + iunit=iunit+1 + icfopg1=iunit + lfopg1=len_trim(cfopg1) + call baopenwa(icfopg1,cfopg1(1:lfopg1),iret) + print *, 'fort.',icfopg1, cfopg1(1:lfopg1) + if(iret.ne.0) then + print *,'there is no output ensemble average = ',cfopg1(1:lfopg1),iret + endif + + iunit=iunit+1 + icfopg2=iunit + lfopg2=len_trim(cfopg2) + call baopenwa(icfopg2,cfopg2(1:lfopg2),iret) + print *, 'fort.',icfopg2, cfopg2(1:lfopg2) + if(iret.ne.0) then + print *,'there is no output ensemble spread = ',cfopg2(1:lfopg2),iret + endif + + ! find grib message, maxgrd: number of grid points in the defined grid + + !do ifile=1,tfiles + do ifile=1,nfiles + if(iskip(ifile).eq.0) then + iids=-9999;ipdt=-9999; igdt=-9999 + idisc=-1; ipdtn=-1; igdtn=-1 + call init_parm(ipdtn,ipdt,igdtn,igdt,idisc,iids) + call getgb2(icfipg(ifile),0,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt,unpack,jskp,gfld,iret) + maxgrd=gfld%ngrdpts + call gf_free(gfld) + firstfile=ifile + if(iret.eq.0) exit + endif + enddo + + if(iret.eq.0) then + + allocate(fgrid(maxgrd,tfiles),fst(tfiles),ens_avg(maxgrd),ens_spr(maxgrd)) + + print *, ' ' + + ! loop over variables + + kids=-9999;kpdt=-9999; kgdt=-9999 + kdisc=-1; kpdtn=-1; kgdtn=-1 + + icount=0 + kskp=0 + do + + inum=0 + fgrid=-9999.9999 + + print *, '----- Start Read Ensemble Forecast ------' + print *, ' ' + + call getgb2(icfipg(firstfile),0,kskp,kdisc,kids,kpdtn,kpdt,kgdtn,kgdt,unpack,kskp,gfldo,iret) + + if(iret.ne.0) then + if(iret.eq.99 ) exit + print *,' getgb2 error = ',iret + cycle + !call errexit(17) + endif + + icount=icount+1 + + if(iret.eq.0) then + inum=inum+1 + call printinfr(gfldo,icount) + fgrid(1:maxgrd,inum)=gfldo%fld(1:maxgrd) + else + print*, 'there is no fcst for ens member 1' + endif + + ! save parameter message for other members + + ipd1=gfldo%ipdtmpl(1) + ipd2=gfldo%ipdtmpl(2) + ipd3=gfldo%ipdtmpl(3) + ipd10=gfldo%ipdtmpl(10) + ipd11=gfldo%ipdtmpl(11) + ipd12=gfldo%ipdtmpl(12) + ipdn=gfldo%ipdtnum + + ! loop over NAEFS members, get operational ensemble forecast + + print *, '----- Ensemble Forecast for Member ------' + print *, ' ' + + ! do imem=2,nfiles + do imem=firstfile+1,nfiles + + if(iskip(ifile).eq.0) then + + iids=-9999;ipdt=-9999; igdt=-9999 + idisc=-1; ipdtn=-1; igdtn=-1 + ipdt(1)=ipd1; ipdt(2)=ipd2; ipdt(10)=ipd10; ipdt(11)=ipd11; ipdt(12)=ipd12 + igdtn=-1; ipdtn=ipdn + call init_parm(ipdtn,ipdt,igdtn,igdt,idisc,iids) + call getgb2(icfipg(imem),0,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt,unpack,jskp,gfld,iret) + + ! print GEFS data message + + if(iret.eq.0) then + inum=inum+1 + call printinfr(gfld,icount) + fgrid(1:maxgrd,inum)=gfld%fld(1:maxgrd) + else + print*, 'there is no fcst for ens member',imem + endif + + call gf_free(gfld) + + endif + + enddo ! end of imem loop + + ! end of imem loop, calculate ensemble mean and spread + + print *, ' '; print *,' variable has member',inum; print *, ' ' + if(inum.gt.navg_min) then + + print *, ' '; print *, ' Combined Ensemble Data Example at Point 8601 ' + write (*,'(10f8.1)') (fgrid(8601,i),i=1,inum) + print *, ' ' + + do n=1,maxgrd + + fst(1:inum)=fgrid(n,1:inum) + + do i=1,inum + weight(i)=1/float(inum) + enddo + + ens_avg(n)=epdf(fst,weight,inum,1.0,0) + ens_spr(n)=epdf(fst,weight,inum,2.0,0) + + enddo + + print *, 'ens_avg(8601)= ',ens_avg(8601) + print *, 'ens_spr(8601)= ',ens_spr(8601) + + print *, ' ' + print *, '----- Output ensemble average and spread for Current Time ------' + print *, ' ' + + ! fnmoc tmax and tmin have message different from ncep, modify them as ncep + + if(gfldo%ipdtnum.eq.11) then + if(gfldo%ipdtmpl(1).eq.0.and.gfldo%ipdtmpl(2).eq.0.and.gfldo%ipdtmpl(27).eq.2) then + gfldo%ipdtmpl(1)=0 + gfldo%ipdtmpl(2)=4 + endif + if(gfldo%ipdtmpl(1).eq.0.and.gfldo%ipdtmpl(2).eq.0.and.gfldo%ipdtmpl(27).eq.3) then + gfldo%ipdtmpl(1)=0 + gfldo%ipdtmpl(2)=5 + endif + endif + + ! when product difinition template 4.1/4.11 chenge to 4.2/4.12 + ! ipdtlen aslo change, need do modification for output + ! code table 4.0, 2=derived forecast + + if(gfldo%ipdtnum.eq.1) then + + temp=-9999 + temp(1:gfldo%ipdtlen)=gfldo%ipdtmpl(1:gfldo%ipdtlen) + + deallocate (gfldo%ipdtmpl) + + gfldo%ipdtnum=2 + if(gfldo%ipdtnum.eq.2) gfldo%ipdtlen=17 + if(gfldo%ipdtnum.eq.2) allocate (gfldo%ipdtmpl(gfldo%ipdtlen)) + + gfldo%ipdtmpl(1:15)=temp(1:15) + gfldo%ipdtmpl(17)=temp(18) + + elseif(gfldo%ipdtnum.eq.11) then + + temp=-9999 + temp(1:gfldo%ipdtlen)=gfldo%ipdtmpl(1:gfldo%ipdtlen) + + deallocate (gfldo%ipdtmpl) + + gfldo%ipdtnum=12 + if(gfldo%ipdtnum.eq.12) gfldo%ipdtlen=31 + if(gfldo%ipdtnum.eq.12) allocate (gfldo%ipdtmpl(gfldo%ipdtlen)) + + gfldo%ipdtmpl(1:15)=temp(1:15) + gfldo%ipdtmpl(17:31)=temp(18:32) + + endif + + ! if(gfldo%ipdtnum.eq.1) ipdtnum_out=2 + ! if(gfldo%ipdtnum.eq.11) ipdtnum_out=12 + ! call change_template4(gfldo%ipdtnum,ipdtnum_out,gfldo%ipdtmpl,gfldo%ipdtlen) + + ! extensions for ensemble mean + + gfldo%ipdtmpl(16)=0 ! code table 4.7, 0=unweighted mean of all Members + gfldo%ipdtmpl(17)=inum ! template 4.2, number of forecast in the ensemble + + ! get the number of bits + ! gfldo%idrtmpl(3) : GRIB2 DRT 5.40 decimal scale factor + + write(6,*) 'gfldo%idrtmpl(3)=',gfldo%idrtmpl(3) + + call gtbits(0,gfldo%idrtmpl(3),maxgrd,0,ens_avg,gmin,gmax,nbit) + + ! gfldo%idrtmpl(4) : GRIB2 DRT 5.40 number of bits + + gfldo%idrtmpl(4)=nbit + + gfldo%fld(1:maxgrd)=ens_avg(1:maxgrd) + + print *, '----- Ensemble Average for Current Time ------' + call putgb2(icfopg1,gfldo,jret) + call printinfr(gfldo,icount) + + ! extensions for ensemble spread + + gfldo%ipdtmpl(16)=2 ! code table 4.7, 2=standard deviation w.r.t cluster mean + gfldo%ipdtmpl(17)=inum ! template 4.2, number of forecast in the ensemble + + ! get the number of bits + ! gfldo%idrtmpl(3) : GRIB2 DRT 5.40 decimal scale factor + + call gtbits(0,gfldo%idrtmpl(3),maxgrd,0,ens_spr,gmin,gmax,nbit) + + ! gfldo%idrtmpl(4) : GRIB2 DRT 5.40 number of bits + + gfldo%idrtmpl(4)=nbit + + gfldo%fld(1:maxgrd)=ens_spr(1:maxgrd) + + print *, '----- Ensemble Spread for Current Time ------' + call putgb2(icfopg2,gfldo,jret) + call printinfr(gfldo,icount) + + ! end of probability forecast calculation + + endif + + call gf_free(gfldo) + + enddo + + ! end of ivar loop + + ! close files + + do ifile=1,nfiles + call baclose(icfipg(ifile),iret) + enddo + + call baclose(icfopg1,iret) + call baclose(icfopg2,iret) + + print *,'Probability Calculation Successfully Complete' + + stop + + else !if(iret.ne.0) then + print *,'there is no maxgrd information' + endif !if(iret.eq.0) then +endif !if(nfiles.gt.2) then + +print *, 'There is not Enough Files Input, Stop!' +call errmsg('There is not Enough Files Input, Stop!') +!call errexit(1) + +stop +end + diff --git a/src/ensstat.fd/gtbits.f90 b/src/ensstat.fd/gtbits.f90 new file mode 100644 index 00000000..1fa7b22f --- /dev/null +++ b/src/ensstat.fd/gtbits.f90 @@ -0,0 +1,98 @@ + SUBROUTINE GTBITS(IBM,IDS,LEN,MG,G,GMIN,GMAX,NBIT) +!$$$ SUBPROGRAM DOCUMENTATION BLOCK +! +! SUBPROGRAM: GTBITS COMPUTE NUMBER OF BITS AND ROUND FIELD. +! PRGMMR: IREDELL ORG: W/NMC23 DATE: 92-10-31 +! +! ABSTRACT: THE NUMBER OF BITS REQUIRED TO PACK A GIVEN FIELD +! AT A PARTICULAR DECIMAL SCALING IS COMPUTED USING THE FIELD RANGE. +! THE FIELD IS ROUNDED OFF TO THE DECIMAL SCALING FOR PACKING. +! THE MINIMUM AND MAXIMUM ROUNDED FIELD VALUES ARE ALSO RETURNED. +! GRIB BITMAP MASKING FOR VALID DATA IS OPTIONALLY USED. +! +! PROGRAM HISTORY LOG: +! 92-10-31 IREDELL +! +! USAGE: CALL GTBITS(IBM,IDS,LEN,MG,G,GMIN,GMAX,NBIT) +! INPUT ARGUMENT LIST: +! IBM - INTEGER BITMAP FLAG (=0 FOR NO BITMAP) +! IDS - INTEGER DECIMAL SCALING +! (E.G. IDS=3 TO ROUND FIELD TO NEAREST MILLI-VALUE) +! LEN - INTEGER LENGTH OF THE FIELD AND BITMAP +! MG - INTEGER (LEN) BITMAP IF IBM=1 (0 TO SKIP, 1 TO KEEP) +! G - REAL (LEN) FIELD +! +! OUTPUT ARGUMENT LIST: +! GROUND - REAL (LEN) FIELD ROUNDED TO DECIMAL SCALING +! (SET TO ZERO WHERE BITMAP IS 0 IF IBM=1) +! GMIN - REAL MINIMUM VALID ROUNDED FIELD VALUE +! GMAX - REAL MAXIMUM VALID ROUNDED FIELD VALUE +! NBIT - INTEGER NUMBER OF BITS TO PACK +! +! SUBPROGRAMS CALLED: +! ISRCHNE - FIND FIRST VALUE IN AN ARRAY NOT EQUAL TO TARGET VALUE +! +! ATTRIBUTES: +! LANGUAGE: CRAY FORTRAN +! +!$$$ +! DIMENSION MG(LEN),G(LEN),GROUND(LEN) + +! implicit none + + INTEGER MG(LEN),IBM,IDS,LEN,NBIT + INTEGER I,I1 + REAL G(LEN),GROUND(LEN),GMAX,GMIN,DS +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +! ROUND FIELD AND DETERMINE EXTREMES WHERE BITMAP IS ON + DS=10.**IDS + IF(IBM.EQ.0) THEN + GROUND(1)=NINT(G(1)*DS)/DS + GMAX=GROUND(1) + GMIN=GROUND(1) + DO I=2,LEN +! print *, ' I G(I)=',G(I) + GROUND(I)=NINT(G(I)*DS)/DS +! print *, ' I GROUND(I) =',GROUND(I) + GMAX=MAX(GMAX,GROUND(I)) +! print *, ' I GMAX=',GMAX + GMIN=MIN(GMIN,GROUND(I)) +! print *, ' I GMIN=',GMIN +! print *, ' I ground,gmin,gmax=',I,GROUND(I),GMAX,GMIN + ENDDO + ELSE + I1=ISRCHNE(LEN,MG,1,0) + IF(I1.GT.0.AND.I1.LE.LEN) THEN + DO I=1,I1-1 + GROUND(I)=0. + ENDDO + GROUND(I1)=NINT(G(I1)*DS)/DS + GMAX=GROUND(I1) + GMIN=GROUND(I1) + DO I=I1+1,LEN + IF(MG(I).NE.0) THEN + GROUND(I)=NINT(G(I)*DS)/DS + GMAX=MAX(GMAX,GROUND(I)) + GMIN=MIN(GMIN,GROUND(I)) + ELSE + GROUND(I)=0. + ENDIF + ENDDO + ELSE + DO I=1,LEN + GROUND(I)=0. + ENDDO + GMAX=0. + GMIN=0. + ENDIF + ENDIF +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +! COMPUTE NUMBER OF BITS + NBIT=LOG((GMAX-GMIN)*DS+0.9)/LOG(2.)+1. +! +! GET REAL FIELD ROUNDED TO DECIMAL SCALING + +! G(1:LEN)=GROUND(1:LEN) +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + RETURN + END diff --git a/src/ensstat.fd/init_parm.f90 b/src/ensstat.fd/init_parm.f90 new file mode 100644 index 00000000..a079e641 --- /dev/null +++ b/src/ensstat.fd/init_parm.f90 @@ -0,0 +1,57 @@ +subroutine init_parm(ipdtn,ipdt,igdtn,igdt,idisc,iids) + +! SUBPROGRAM: init_parm +! +! PRGMMR: Bo Cui DATE: 2013-06-11 +! +! ABSTRACT: This subroutine returns the Grid Definition, and +! Product Definition for a given data field. +! +! PROGRAM HISTORY LOG: +! 2013-02-27 Bo Cui +! +! USAGE: call init_parm(ipdtn,ipdt,igdtn,igdt) +! +! INPUT: ipdtn,ipdt,igdtn,igdt +! OUTPUT: jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt +! +! REMARKS: None +! +! ATTRIBUTES: +! LANGUAGE: Fortran 90 +! MACHINE: IBM SP +! +!$$$ + +implicit none + +integer,dimension(200) :: jids,jpdt,jgdt,iids,ipdt,igdt +integer jskp,jdisc,jpdtn,jgdtn,idisc,ipdtn,igdtn +common /param/jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt + +jskp=0 + +jdisc=-1 +jpdtn=ipdtn +jgdtn=igdtn + +jids=-9999 +jpdt=-9999 +jgdt=-9999 + +! input product defining values + +jids=iids +jpdt=ipdt +jgdt=igdt + +!jpdt(1)=ipdt(1) +!jpdt(2)=ipdt(2) +!jpdt(10)=ipdt(10) +!jpdt(11)=ipdt(11) +!jpdt(12)=ipdt(12) + +return +end + + diff --git a/src/ensstat.fd/isrchne.f90 b/src/ensstat.fd/isrchne.f90 new file mode 100644 index 00000000..89cdeb87 --- /dev/null +++ b/src/ensstat.fd/isrchne.f90 @@ -0,0 +1,45 @@ + FUNCTION ISRCHNE(N,X,INCX,TARGET) +!$$$ SUBPROGRAM DOCUMENTATION BLOCK +! . . . . +!SUBPROGRAM: ISRCHNE Searches vector given a target +! PRGMMR: gilbert ORG: W/NP11 DATE: 99-02-11 +! +!ABSTRACT: Searches a vector for the first element +! not equal to a target +! +!PROGRAM HISTORY LOG: +! 99-02-11 Gilbert +! +!USAGE: index=ISRCHNE(n, x, incx, target) +! INPUT ARGUMENT LIST: +! n - Number of elements to be searched +! x - Real or integer array of dimension (n-1) * |incx| + 1. +! Array x contains the vector to be searched. +! incx - Increment between elements of the searched array. +! target - Value for which to search in the array. +! +! OUTPUT VALUE +! index - Index of the first element equal or not equal to target. If +! target is not found, n+1 is returned. If n <= 0, 0 is +! returned. +! +!REMARKS: This code and documentation was taken directly from the +! man page for routine ISRCHNE on a CRAY UNICOS system. +! +!ATTRIBUTES: +! LANGUAGE: Fortran +! +!$$$ + INTEGER X(*), TARGET + J=1 + ISRCHNE=0 + IF(N.LE.0) RETURN + IF(INCX.LT.0) J=1-(N-1)*INCX + DO 100 I=1,N + IF(X(J).NE.TARGET) EXIT + J=J+INCX + 100 CONTINUE + 200 ISRCHNE=I + RETURN + END + diff --git a/src/ensstat.fd/printinfr.f90 b/src/ensstat.fd/printinfr.f90 new file mode 100644 index 00000000..3727b52b --- /dev/null +++ b/src/ensstat.fd/printinfr.f90 @@ -0,0 +1,82 @@ +subroutine printinfr(gfld,ivar) + +! SUBPROGRAM: printinfr +! +! PRGMMR: Bo Cui DATE: 2013-06-11 +! +! USAGE: print grib2 data information +! +! INPUT: gfld,ivar + +use grib_mod +use params + +implicit none + +type(gribfield) :: gfld + +integer,dimension(200) :: jids,jpdt,jgdt,ipdt,igdt +integer :: currlen=0 +logical :: unpack=.true. + +integer kf,j,ivar,i +real fldmin,fldmax + +!kf=gfld%ngrdpts +kf=gfld%ndpts + +fldmin=gfld%fld(1) +fldmax=gfld%fld(1) + +do j=2,kf + if (gfld%fld(j).gt.fldmax) fldmax=gfld%fld(j) + if (gfld%fld(j).lt.fldmin) fldmin=gfld%fld(j) +enddo + +! print out + +! gfld%ipdtnum 1/11: ens. fcst or control, high reslution +! gfld%ipdtnum 2/12: ens. average fcst +! gfld%ipdtnum 0/8: cdas reanalysis + +! gfld%ipdtnum 11: ens. fcst or control in a continuous or non-continuous time interval +! gfld%ipdtnum 12: derived fcst based on ens. members in a continuous/non-continuous time interval +! gfld%ipdtnum 8: Rstatistically processed values in a continuous/non-continuous time interval + +if(gfld%ipdtnum.eq.11.or.gfld%ipdtnum.eq.12) then + write(6,100) + write(6,102) ivar,gfld%ipdtnum,(gfld%ipdtmpl(i),i=1,3),gfld%ipdtmpl(10),gfld%ipdtmpl(12), & + (gfld%idsect(i),i=6,9),gfld%ipdtmpl(9),gfld%ipdtmpl(30), & + (gfld%ipdtmpl(i),i=16,17), & + kf,fldmax,fldmin,gfld%fld(8601) + write(6,*) + +elseif(gfld%ipdtnum.eq.8.or.gfld%ipdtnum.eq.0) then + write(6,300) + write(6,302) ivar,gfld%ipdtnum,(gfld%ipdtmpl(i),i=1,3),gfld%ipdtmpl(10),gfld%ipdtmpl(12), & + (gfld%idsect(i),i=6,9),gfld%ipdtmpl(9), & + kf,fldmax,fldmin,gfld%fld(8601) + write(6,*) + +else + write(6,200) + write(6,202) ivar,gfld%ipdtnum,(gfld%ipdtmpl(i),i=1,3),gfld%ipdtmpl(10),gfld%ipdtmpl(12), & + (gfld%idsect(i),i=6,9),gfld%ipdtmpl(9),(gfld%ipdtmpl(i),i=16,17), & + kf,fldmax,fldmin,gfld%fld(8601) + write(6,*) +endif + +100 format(' REC PDTN PD1 PD2 PD3 PD10 PD12 YEAR MN DY HR FHR TR ', & + 'E16 E17 LEN MAX MIN EXAMPLE') +102 format(i4,i5,4i4,i8,i6,3i3,2i4,2i4,i8,3f11.2) + +200 format(' REC PDTN PD1 PD2 PD3 PD10 PD12 YEAR MN DY HR FHR ', & + 'E16 E17 LEN MAX MIN EXAMPLE') +202 format(i4,i5,4i4,i8,i6,3i3,3i4,i8,3f11.2) + +300 format(' REC PDTN PD1 PD2 PD3 PD10 PD12 YEAR MN DY HR FHR ', & + ' LEN MAX MIN EXAMPLE') +302 format(i4,i5,4i4,i8,i6,3i3,i4,8x,i8,3f11.2) + +return +end diff --git a/src/wave_stat.fd/CMakeLists.txt b/src/wave_stat.fd/CMakeLists.txt new file mode 100644 index 00000000..541c40fd --- /dev/null +++ b/src/wave_stat.fd/CMakeLists.txt @@ -0,0 +1,12 @@ +list(APPEND fortran_src + grbit2.f90 + wave_stat.f90 +) + +set(exe_name wave_stat.x) +add_executable(${exe_name} ${fortran_src}) +target_link_libraries(${exe_name} PRIVATE bacio::bacio_4 + w3emc::w3emc_4 + g2::g2_d) + +install(TARGETS ${exe_name} RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}) diff --git a/src/wave_stat.fd/grbit2.f90 b/src/wave_stat.fd/grbit2.f90 new file mode 100644 index 00000000..fa6631c3 --- /dev/null +++ b/src/wave_stat.fd/grbit2.f90 @@ -0,0 +1,263 @@ +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +subroutine grbit2 (IFL2,nx,ny,lcgrib,rdlon,rdlat,bmp,fld,ierr,ymdc,fhr,& + parcode,ensid,nprb,prbid,prblv,nensm) +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +! +! Generates grib2 output +! +! Vera Gerald, NCEP June 2011 +! +! Changes: +! - Add parcode to read in discipline, category and parameter number +! +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +! Inputs +! +! parcode - three integers defining discipline, category and param number +! from grib2 tables +! ensid - 1, 2 or 3 ; where: 1=mean , 2=spread and 3=probability +! nprb - Total number of probability levels +! prbid - Index of probability value +! prblv - Value of probability threshold +! nensm - number of ensemble members +! IFL2 - Unit number for output file +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +! + CHARACTER(len=1) :: cgrib(lcgrib),cin(lcgrib) + real :: FLD(nx*ny) + real :: coordlist(1) + real :: rdlon(3),rdlat(3) + integer :: pds(25),gds(22),ensid,prbid,prblv,nensm + integer :: listsec0(3),tmpln + integer :: listsec1(13) + integer :: parcode(3) + integer :: igds(5),igdstmpl(200),ipdstmpl(200) + integer :: ymdc,ymd,ym,y4,mm,dd,cc,fhr,yy,cen + integer :: nx, ny, nxny, ibfl + integer :: ideflist,idefnum + integer :: idrstmpl(200) + integer :: intlon,intlat + Logical*1 :: bmp(nx*ny) + character*11 :: envvar + character(len=80):: g2file +! +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +! +! Initialize date parameters +! + cc = mod(ymdc,100) + ymd = ymdc/100 + dd = mod(ymd,100) + ym = ymd/100 + mm = mod(ym,100) + y4 = ym/100 +! +! Create grib message +! +! Initialize new GRIB2 message and pack +! +! Section 0 (Indicator Section) +! + listsec0 = 00 +! + listsec0(1) = parcode(1) ! Discipline-GRIB Master Table Number (see Code Table 0.0) + listsec0(2) = 2 ! GRIB Edition Number (currently 2) + listsec0(3) = 0 +! +! Section 1 (Identification Section) +! + listsec1(1) = 7 ! 7 Id of orginating centre (Common Code Table C-1) + listsec1(2) = 0 !"EMC"! Id of orginating sub-centre (local table)/Table C/on388 + listsec1(3) = 2 ! GRIB Master Tables Version Number (Code Table 1.0) + listsec1(4) = 1 !per Brent! GRIB Local Tables Version Number (Code Table 1. + listsec1(5) = 1 ! Significance of Reference Time (Code Table 1.2) + listsec1(6) = y4 ! Reference Time - Year (4 digits) + listsec1(7) = mm ! Reference Time - Month + listsec1(8) = dd ! Reference Time - Day + listsec1(9) = cc ! Reference Time - Hour(cycle:0,6,12,18) + listsec1(10) = 0 ! Reference Time - Minute + listsec1(11) = 0 ! Reference Time - Second + listsec1(12) = 0 ! Production status of data (Code Table 1.3) + listsec1(13) = 1 ! Type of processed data (Code Table 1.4) + +! + call gribcreate(cgrib,lcgrib,listsec0,listsec1,ierr) + if (ierr.ne.0) then + write(6,*) ' ERROR creating new GRIB2 field = ',ierr + return + endif +! +! Section 3 (Grid Definition Section) +! + + nxny = nx * ny ! NI # rows(lat) & NJ # rows long +! + igds = 00 +! + igds(1) = 0 !Source of grid definition (see Code Table 3.0) + igds(2) = nxny ! num of grid points + igds(3) = 0 !Number of octets needed for each additional grid points dfn + igds(4) = 0 !Interpretation of list for optional points definition (Code Table 3.11) + igds(5) = 0 !Grid Definition Template Number (Code Table 3.1) +! + idefnum = 0 + ideflist=0 !Used if igds(3) .ne. 0. Dummy array otherwise +! + igdstmpl = 00 +! +! Define geographical ranges from arguments +! + iloni=INT(rdlon(1)*1.E6) + ilone=INT(rdlon(2)*1.E6) + idlon=INT(rdlon(3)*1.E6) + ilati=INT(rdlat(1)*1.E6) + ilate=INT(rdlat(2)*1.E6) + idlat=INT(rdlat(3)*1.E6) + + igdstmpl(1) = 6 !Earth assumed spherical with radius of 6,371,229.0m + igdstmpl(2) = 0 ! + igdstmpl(3) = 0 + igdstmpl(4) = 0 + igdstmpl(5) = 0 + igdstmpl(6) = 0 + igdstmpl(7) = 0 + igdstmpl(8) = nx ! num points along parallel + igdstmpl(9) = ny ! num points along meridian + igdstmpl(10) = 0 + igdstmpl(11) = 0 + igdstmpl(12) = ilati ! lat of first grid point + igdstmpl(13) = iloni ! long of first grid point + igdstmpl(14) = 48 ! res and comp flags + igdstmpl(15) = ilate ! lat of last grid point + igdstmpl(16) = ilone ! long of last grid point + igdstmpl(17) = idlat ! Increment of lat + igdstmpl(18) = idlon ! Increment of long + igdstmpl(19) = 0 ! scanning mode + igdstmpl(21) = 0 + igdstmpl(22) = 0 + +! + call addgrid(cgrib,lcgrib,igds,igdstmpl,200,ideflist, & + idefnum,ierr) + if (ierr.ne.0) then + write(6,*) ' ERROR adding GRIB2 grid = ',ierr + return + endif +! +! Section 4 - Product Definition Section +! +! Define production template number (4.tmpln) +! + if(ensid.eq.3)then + ipdsnum = 5 ! ensemble for probal template(4.5) + else + ipdsnum = 2 ! ensemble spread or mean template(4.2) + endif +! +! Parameter category (see Code Table 4.1) +! + ipdstmpl = 00 +! + ipdstmpl(1) = parcode(2) ! Category +! +! Get parameter number (see Code Table 4.2) +! + ipdstmpl(2) = parcode(3) ! Parameter number +! +! Type of generating process: analysis or forecast(see code Table 4.3) +! For Wave Model/Forecast fields +! + ipdstmpl(3) = 4 ! ensemb + ipdstmpl(4) = 0 !background generating process identifier +! (defined by originating Center) + ipdstmpl(5) = 10 ! :analysis or forecast generating process identifier +! (defined by originating Center) + ipdstmpl(6) = 0 ! hours of observational data cutoff after reference time + ipdstmpl(7) = 0 ! minutes of observational data cutoff after reference time + ipdstmpl(8) = 1 !indicator of unit of time range (see Code Table 4.4) + ipdstmpl(9) = fhr !forecast time in units defined by pdstmpl(8) + ipdstmpl(10) = 1 ! type of level (see Code Table 4.5) 1st level + ipdstmpl(11) = 0 ! scale factor of pdstmpl(10) + ipdstmpl(12) = 0 ! scaled value of pdstmpl(10) + ipdstmpl(13) = 255 ! type of level (See Code Table 4.5) 2nd level + ipdstmpl(14) = 0 ! scale factor of ipdstmpl(13) + ipdstmpl(15) = 0 ! scaled value of ipdstmpl(13) +! +! Choose proper parameters for given ensid (mean, spread or probab) +! + if (ensid.eq.1) then ! ensemble mean + ipdstmpl(16) = 0 ! + ipdstmpl(17) = nensm ! Number of esemble members + elseif (ensid.eq.2) then ! ensemble spread + ipdstmpl(16) = 4 + ipdstmpl(17) = nensm ! Number of ensemble members + elseif (ensid.eq.3) then ! ensemble probability + ipdstmpl(16) = prbid ! Forecast probability number + ipdstmpl(17) = 8 ! Total number of forecast probabilities + ipdstmpl(18) = 1 ! Probability type + ipdstmpl(19) = 0 ! Scale factor of lower limit + ipdstmpl(20) = 0 ! Scaled value of lower limit + ipdstmpl(21) = 2 ! Scale factor of upper limit + ipdstmpl(22) = prblv ! Scaled value of upper limit + endif +! + numcoord=0 + coordlist= 0.0 !needed for hybrid vertical coordinate +! set bitmap flag + ibfl = 192 ! GDS/BMS FLAG (RIGHT ADJ COPY OF OCTET 8) + if (btest(ibfl,6)) then + ibmap=0 + else + ibmap=255 ! Bitmap indicator ( see Code Table 6.0 ) + bmp=.true. + endif + idrsnum=40 ! Data Rep Template Number ( see Code Table 5.0 )/Simple packing +! + idrstmpl=0 +! +!******************************************************************************** +! idrstmpl(1): reference value(R) (IEEE 32-bit floating-point value) * +! idrstmpl(2): binary scale factor (E) * +! idrstmpl(3): decimal scale factor (D) * +! idrstmpl(4): number of bits used for each packed value for simple packing * +! or for each group reference value for complex packing or * +! spatial differencing * +! idrstmpl(5): type of original field values (See Code Table 5.1) * +!******************************************************************************** +! + idrstmpl(1) = 0 + idrstmpl(2) = 0 + + idrstmpl(3) = 2 ! binary scale factor (E) + idrstmpl(4) = 0 + idrstmpl(5) = 0 +! + call addfield(cgrib,lcgrib,ipdsnum,ipdstmpl,200, & + coordlist,numcoord,idrsnum,idrstmpl,200,fld, & + nxny,ibmap,bmp,ierr) + if (ierr.ne.0) then + write(6,*) ' ERROR adding GRIB2 field = ',ierr + return + endif +! +! + print*,'grid2 PDS ',ipdsnum,idrsnum,ipdstmpl(1:25) +! + print*,'grib2 GDS',igds(1:5),igdstmpl(1:22) +! +! +! End GRIB2 field +! + call gribend(cgrib,lcgrib,lengrib,ierr) + if (ierr.ne.0) then + write(6,*) ' ERROR ending new GRIB2 message = ',ierr + return + endif + print *,' writing ',lengrib,' bytes...' + call wryte(ifl2,lengrib,cgrib) +! +!.. encode next record +! + return +! + end diff --git a/src/wave_stat.fd/wave_stat.f90 b/src/wave_stat.fd/wave_stat.f90 new file mode 100644 index 00000000..fa122de5 --- /dev/null +++ b/src/wave_stat.fd/wave_stat.f90 @@ -0,0 +1,348 @@ +!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -+ + + ++ +! multiwave_gwes_stats +! +! Origin: H.S. Chen, 2004 +! Upgrades: D. Cao, 2008 +! NFCENS upgrade: program nfcombwave_ensemble, Y.Y. Chao, Jun 2011 +! GWES upgrade: J-H. Alves, Jan 2014 +! +! Latest changes (GWES upgrade) +! - Became multiwave_gwes_stats +! - Resolution increased to 0.5 degrees +! - Allow output grid to be different than input grid (use +! allocatable arrays to get resolutions from input grib2 data) +! - Buoy output using wgrib2 to extract directly from input data +! - Parallelization of stats using mpiserial acroos parameters and +! statistics +! +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - --+ + + ++ + module ens_common +! + implicit none +! + integer,parameter :: iu11=11,iu12=12 + integer,parameter :: iu51=51,iu52=52,iu53=53,iu54=54,iu59=59 + real,parameter :: defv=1.0e10,tolr=1.0e-8,zero=0.0 +! + character :: parname*5, statstype*6 + integer :: kpds(25),kgds(22) + integer :: ymdc,fhr,nnme,nnsc + integer :: parcode(3) + integer :: istop + real,allocatable :: dumy1d(:) +! Grid + real :: rdlon(3),rdlat(3) + integer :: nlola(2) +! + character(len=13) :: fname + logical*1,allocatable :: lbms2d(:,:),lbms1d(:) + real,allocatable :: fsum(:,:),fmean(:,:),fspre(:,:) + character,allocatable :: id_member(:)*2 + character(len=7),allocatable :: id_file(:) + real,allocatable :: f(:,:,:),fprob(:,:,:),scale(:) +! + real :: gmean,gspre + real,allocatable :: g(:),gprob(:) +! +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - --+ + + ++ + endmodule ens_common +! +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - --+ + + ++ +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - --+ + + ++ + PROGRAM gwes_stats +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -+ + + ++ +!$$$ MAIN PROGRAM DOCUMENTATION BLOCK +! . . . . +! MAIN PROGRAM: ensemble conduct ensemble statistics. +! PRGMMR: H.S. CHEN ORG: W/MP21 DATE: 04-05-20 +! +! ABSTRACT: conduct ensemble statistics. +! +! PROGRAM HISTORY LOG: +! 04-05-20 H.S. Chen Origination and implement at NCEP. +! +! USAGE: +!** Interface. +! ---------- +! The program runs independently and must run after the +! completion of all individual ensemble POST runs +! +! Input File: +! ------------ +! gwes.stats.inp: Provides date, parameter definitions, ensemble member +! names and numbers, probability level and input data +! file names (data_NN) +! Output File: +! ------------- +! mean_out - grib2 file containing combined ensemble mean +! spread_out - grib2 file containing combined ensemble spread +! prob_out - grib2 file containing combined ensemble probabilities of +! Hs exceeding a given value +! +! Method. +! ------- +! Interpolation using wgrib2 +! +! Externals. +! ---------- +! wgrib2 +! +! ATTRIBUTES: +! LANGUAGE: FORTRAN 90 +! +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - --+ + + ++ +! + use ens_common +! + integer :: ierr,i,j,k,l,m,n,i1,i2,j1,j2,iret,leve +! +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - --+ + + ++ +! CALL W3TAGB('nww2ec ',0097,0027,0075,'NP21 ') +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - --+ + + ++ +! + open (unit=iu59,file='test_out') +! + read(*,*) ymdc,fhr,parname,parcode + write(iu59,'(i10.10,2x,i3.3,2x,a5,i3,1x,i3,1x,i3)') ymdc,fhr,parname,parcode +! read(*,*) statstype +! write(iu59,'(a6)') statstype + read(*,*) nnme + write(iu59,'(i4)') nnme + allocate( id_member(nnme) ) + read(*,*) id_member(:) + write(iu59,'(20(1x,a2))') id_member(:) + read(*,*) nnsc + write(iu59,'(i4)') nnsc +! + allocate( scale(nnsc) ) + scale(:) = 0.0 + read(*,*) scale(:) + write(iu59,'(10f8.3)') scale(:) +! + read(*,*) nlola + write(iu59,'(i6,1x,i6)') nlola + nx=nlola(1) + ny=nlola(2) + nxy=nx*ny + lcgrib=4*nxy + allocate (dumy1d(nxy),lbms1d(nxy),lbms2d(nx,ny),fsum(nx,ny), & + fmean(nx,ny),fspre(nx,ny)) +! + read(*,*) rdlon + write(iu59,'(3f9.4)') rdlon +! + read(*,*) rdlat + write(iu59,'(3f9.4)') rdlat +! + allocate( f(nx,ny,nnme) ) + allocate( fprob(nx,ny,nnsc) ) +! + allocate( g(nnme) ) + allocate( gprob(nnsc) ) +! + allocate( id_file(nnme) ) + read(*,*) id_file(:) + write(iu59,'(10(1x,a7))') id_file(:) +! +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -+ + + ++ +! + f(:,:,:) = 0.0 + do m=1,nnme +! + fname=id_file(m) + write(iu59,*) fname +! + open (iu11,file=fname,access='SEQUENTIAL',form='UNFORMATTED') +! + write(iu59,'("*** read f(m), m =",i4)') m + read (iu11,end=12,err=91,iostat=ierr) & + ((f(i,j,m),i=1,nx),j=ny,1,-1) +! + where( f(:,:,m).ge.defv ) f(:,:,m) = zero +! + 12 continue + close(iu11) + enddo +! +! Calculate ensemble statistics. +! + call ens_statistics +! +! Output. +! +! Set up bitmap, lbms2d(:,:) and lbms1d(:), for missing data indicator; +! true for data, false for no data. +! + lbms2d(:,:) = .true. + where( fmean(:,:).le.zero ) lbms2d(:,:) = .false. +! + ij=1 + do j=1,ny + do i=1,nx + lbms1d(ij)=lbms2d(i,j) + ij=ij+1 + enddo + enddo +! +! ensemble mean. +! +! if (statstype .eq. 'avrage') then + call baopen(iu52,'mean_out',iret) + write(*,'("** After baopen for mean, iret = ",i5)') iret + dumy1d(:)=0.0 + ij=1 + do j=1,ny + do i=1,nx + dumy1d(ij)=fmean(i,j) + ij=ij+1 + enddo + enddo +! + call grbit2 (iu52,nx,ny,lcgrib,rdlon,rdlat,lbms1d,dumy1d,iret, & + ymdc,fhr,parcode,1,0,0,0,nnme) + if( iret.ne.0 ) then + write(*,'("**** ERR in grbit2, ymdc,fhr,parcode,parname: ")') + write(*,'(1x,i10.10,1x,i3.3,i5,1x,a5)') ymdc,fhr,parcode,parname + endif +! endif +! +! ensemble spread. +! +! if (statstype .eq. 'spread') then + call baopen(iu53,'spread_out',iret) + write(*,'("** After baopen for spread, iret = ",i5)') iret + dumy1d(:)=0.0 + ij=1 + do j=1,ny + do i=1,nx + dumy1d(ij)=fspre(i,j) + ij=ij+1 + enddo + enddo +! + call grbit2 (iu53,nx,ny,lcgrib,rdlon,rdlat,lbms1d,dumy1d,iret, & + ymdc,fhr,parcode,2,0,0,0,nnme) + if( iret.ne.0 ) then + write(*,'("**** ERR in grbit2, ymdc,fhr,parcode,parname: ")') + write(*,'(1x,i10.10,1x,i3.3,i5,1x,a5)') ymdc,fhr,parcode,parname + endif +! endif +! +! ensemble probability. +! +! if (statstype .eq. 'probab') then + call baopen(iu54,'prob_out',iret) + write(*,'("** After baopen for prob, iret = ",i5)') iret + do l=1,nnsc +! + dumy1d(:)=0.0 + ij=1 + do j=1,ny + do i=1,nx + dumy1d(ij)=fprob(i,j,l) + ij=ij+1 + enddo + enddo +! + leve = nint(100.*scale(l)) + call grbit2 (iu54,nx,ny,lcgrib,rdlon,rdlat,lbms1d,dumy1d, & + iret,ymdc,fhr,parcode,3,nnsc,l,leve,nnme) + if( iret.ne.0 ) then + write(*,'("** ERR in grbit2, ymdc,fhr,parcode,parname,leve:")') + write(*,'(1x,i10.10,1x,i3.3,i5,1x,a5,i7)') & + ymdc,fhr,parcode,parname,leve + endif + enddo +! endif +! +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -+ + + ++ +! + 91 continue +! +! deallocate. +! + deallocate( id_member ) + deallocate( id_file,scale,fprob,f ) + deallocate( g,gprob ) + deallocate (dumy1d,lbms1d,lbms2d,fsum,fmean,fspre) +! + close(iu51) +! + call baclose(iu52,iret) + write(*,'("** After baclose for mean, iret= ",i5)') iret + call baclose(iu53,iret) + write(*,'("** After baclose for spread, iret= ",i5)') iret + call baclose(iu54,iret) + write(*,'("** After baclose for prob, iret= ",i5)') iret + close(iu59) +! + STOP + END +! +! + subroutine ens_statistics +!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -+ + + ++ +!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -+ + + ++ +! Calculate ensemble statistics. +! input : nx = number of x points. +! ny = number of y points. +! me = number of ensemble members. +! nnsc = number of scales for probability. +! f(nx,ny,me) = field array. +! fsum(nx,ny) = only used for calculating sum. +! fscale(nnsc) = scale array. +! statstype = type of statstics to compute (mean, spread, +! prob) - feature added to increase paralellism +! output: fmean(nx,ny) = ensemble mean array. +! fspre(nx,ny) = ensemble spread array. +! fprob(nx,ny,nnsc) = ensemble probability. +!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -+ + + ++ +! + use ens_common +! +! local +! + integer :: m,n,i,j,k +! +!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -+ + + ++ +! +! mean. +! +! if(statstype .eq. 'avrage')then +! write(*,*)statstype + fsum(:,:) = 0.0 + do m=1,nnme + fsum(:,:) = fsum(:,:) + f(:,:,m) + enddo + fmean(:,:) = fsum(:,:)/real(nnme) +! endif +! +! spread - deviation. +! +! if(statstype .eq. 'spread')then +! write(*,*)statstype + fspre(:,:) = 0.0 + do m=1,nnme + fspre(:,:) = fspre(:,:) + (f(:,:,m)-fmean(:,:))**2 + enddo + fspre(:,:) = sqrt( fspre(:,:)/real(nnme) ) +! endif +! +! probability. +! +! if(statstype .eq. 'probab')then +! write(*,*)statstype + fprob(:,:,:) = 0.0 + do n=1,nnsc + do m=1,nnme + where( f(:,:,m).ge.scale(n) ) fprob(:,:,n)=fprob(:,:,n)+1. + enddo + fprob(:,:,n) = fprob(:,:,n)/real(nnme) + enddo +! endif +! +!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -+ + + ++ +! + return + end +!