diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 1ea002b8..ab626748 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -31,3 +31,4 @@ add_subdirectory(ensadd.fd) add_subdirectory(ensppf.fd) add_subdirectory(ensstat.fd) add_subdirectory(wave_stat.fd) +add_subdirectory(gefs_6h_ave_1mem.fd) diff --git a/src/gefs_6h_ave_1mem.fd/CMakeLists.txt b/src/gefs_6h_ave_1mem.fd/CMakeLists.txt new file mode 100644 index 00000000..092d0773 --- /dev/null +++ b/src/gefs_6h_ave_1mem.fd/CMakeLists.txt @@ -0,0 +1,19 @@ +list(APPEND fortran_src + gefs_6h_ave_1mem.f90 + prlevel.f90 + printinfr.f90 +) + +if(CMAKE_Fortran_COMPILER_ID MATCHES "^(Intel)$") + set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -i4 -r8") +elseif(CMAKE_Fortran_COMPILER_ID MATCHES "^(GNU)$") + set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -fdefault-real-8") +endif() + +set(exe_name gefs_6h_ave_1mem.x) +add_executable(${exe_name} ${fortran_src}) +target_link_libraries(${exe_name} PRIVATE bacio::bacio_4 + w3emc::w3emc_d + g2::g2_d) + +install(TARGETS ${exe_name} RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}) diff --git a/src/gefs_6h_ave_1mem.fd/gefs_6h_ave_1mem.f90 b/src/gefs_6h_ave_1mem.fd/gefs_6h_ave_1mem.f90 new file mode 100755 index 00000000..2e7fd0d4 --- /dev/null +++ b/src/gefs_6h_ave_1mem.fd/gefs_6h_ave_1mem.f90 @@ -0,0 +1,473 @@ +program gefs_6h_ave_1mem +! main program: gefs_6h_ave_1mem +! Author: Hong Guan :2018-12-17 +! REF: Eric Sinsky, Wei Li, Yali Mao, Bo Cui +! Purpose: calculate accumulation or avrage (03Z,06Z) +! +! usage: daily_ave_acc.exe yyyymmdd +! +! input file: control f00 and f03 reanalysis data and ensemble forecast (master file) +! output file: accumulation + +! 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 +! +!$$$ + + use grib_mod + use params + implicit none + + type(gribfield) :: gfld,gfldo + +!type pdt_t0 +! integer :: npdt0 ! Product Definition Template Number +! integer :: icat0 ! Parameter Category by Product Discipline +! integer :: iprm0 ! Parameter Number by Product Discipline and Parameter Category +! integer :: igp0 ! Type of Generating Process +! integer :: iffs0 ! Type of first fixed surface +! integer :: isf0 ! Scale factor of first fixed surface +! integer :: isv0 !Scaled value of first fixed surface +!end type pdt_t0 +! +!! PDT parameters in the input GRIB2 file (template 4 number, category, parameter, type of level) +!type(pdt_t0), parameter :: & +!in pgrb2a +! pdt_hgt500 = pdt_t0(1, 3, 5, 4, 100, 0, 50000), & +! pdt_hgt200 = pdt_t0(1, 3, 5, 4, 100, 0, 20000), & +! pdt_ugrd200 = pdt_t0(1, 2, 2, 4, 100, 0, 20000), & +! pdt_ugrd850 = pdt_t0(1, 2, 2, 4, 100, 0, 85000), & +! pdt_vgrd200 = pdt_t0(1, 2, 3, 4, 100, 0, 20000), & +! pdt_vgrd850 = pdt_t0(1, 2, 3, 4, 100, 0, 85000), & +! pdt_t2m = pdt_t0(1, 0, 0, 4, 103, 0, 2), & +! pdt_apcp = pdt_t0(11,1, 8, 4, 1, 0, 0), & +! pdt_ulwrf_top= pdt_t0(11,5,193,4, 8, 0, 0), & +!now in pgrb2a +! pdt_tsfc = pdt_t0(1, 0, 0, 4, 1, 0, 0) + + integer :: nfield,nfield1 + parameter(nfield=42) + parameter(nfield1=39) + + type pdt_t + integer,dimension(nfield):: npdt ! Product Definition Template Number + integer,dimension(nfield):: icat ! Parameter Category by Product Discipline + integer,dimension(nfield):: iprm ! Parameter Number by Product Discipline and Parameter Category + integer,dimension(nfield):: igp ! Type of Generating Process + integer,dimension(nfield):: iffs ! Type of first fixed surface + integer,dimension(nfield):: isf ! Scale factor of first fixed surface + integer,dimension(nfield):: isv !Scaled value of first fixed surface + end type pdt_t + + type(pdt_t) :: pdt + +! data pdt%npdt/ 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8/ +! data pdt%icat/ 5, 5, 4, 4, 4, 5, 2, 2, 0, 0, 1, 1,19, 6, 6, 6, 6, 6, 0, 3, 3, 4, 4, 0, 0, 3, 3, 0, 3, 3, 0, 3, 3, 0, 6, 0, 1, 0, 1, 1, 1/!Parameter Category by Product Discipline +! data pdt%iprm/ 192,193,192,193,193,193,17,18,11,10, 7,196, 1, 1, 1, 5, 4, 3,193,194,195,194,195, 5, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0,193,16,42, 5, 8, 10, 9/!ParameterNumberbyProductDiscipline and Parameter Category +! data pdt%iffs/ 1, 1, 1, 1, 8, 8, 1, 1, 1, 1, 1, 1, 1,200,211,234,224,214, 1, 1, 1, 1, 1,103,103,233,232,233,223,222,223,213,212,213,200, 1, 1, 1, 1, 1, 1/ +! integer,dimension(nfield)::disc +! data disc/0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,0/ + + + data pdt%npdt/ 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8/ + data pdt%icat/ 5, 5, 4, 4, 4, 5, 2, 2, 0, 0, 1, 1,19, 6, 6, 6, 6, 6, 0, 3, 3, 4, 4, 0, 0, 3, 3, 0, 3, 3, 0, 3, 3, 0, 6, 0, 1, 0, 1, 1, 1, 1/!Parameter Category by Product Discipline + data pdt%iprm/ 192,193,192,193,193,193,17,18,11,10, 7,196, 1, 1, 1, 5, 4, 3,193,194,195,194,195, 5, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0,193,16,42, 5, 50, 8, 10, 9/!ParameterNumberbyProductDiscipline and Parameter Category + data pdt%iffs/ 1, 1, 1, 1, 8, 8, 1, 1, 1, 1, 1, 1, 1,200,211,234,224,214, 1, 1, 1, 1, 1,103,103,233,232,233,223,222,223,213,212,213,200, 1, 1, 1, 1, 1, 1, 1/ + integer,dimension(nfield)::disc + data disc/0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,0,0/ + + + integer :: n_time + parameter(n_time=3) !from fhr 00 to 06 + integer :: j,jpdtn,jgdtn + integer,dimension(200) :: jids,jpdt,jgdt + logical :: unpack=.true. + integer :: jdisc ! discipline#, table 0.0(met:0 hydro:1 land:2) + integer :: day,month,year,hour,fhour + character(len=256) :: datafile(n_time),outfile,outfile03,outfile006 + character(len=8) :: file_date + integer :: unit, ifid,ifid1,nx,ny,iret,jret,i,k,ifh,nfi,maxgrd,ifd,ind,ens_id + character(len=3) ::sfh + character(len=2) ::sdy + character(len=2) ::smonth + character(len=4) ::syr + character(len=8) :: pabbrev + character(len=30) :: labbrev + character(len=30) :: labbrev_short + character(len=500) :: datapath = './' + + real, allocatable :: var_save(:,:) ! (maxgrd,nfi) + real, allocatable :: var_save_acc(:) ! (maxgrd) + real, allocatable :: var_save_DSWRF(:) ! (maxgrd) + real, allocatable :: var_save_USWRF(:) ! (maxgrd) + real, allocatable :: apcp_6h(:) ! (maxgrd) + real, allocatable :: acpcp_6h(:) ! (maxgrd) + real, allocatable :: apcp_3h(:) ! (maxgrd) + real, allocatable :: acpcp_3h(:) ! (maxgrd) + real, allocatable :: dpcp(:) ! (maxgrd) + + + call GET_COMMAND_ARGUMENT(1, file_date) +! call get_environment_variable("ens_mem", ens_mem) + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!open files for output grib2 (f03 and f006) + outfile03=trim(datapath)//'gefs'//'.t00z.'//'pgrb2af003' + print *,outfile03 + call baopenwa(300,outfile03,iret) ! for add more than 1 members + outfile006=trim(datapath)//'gefs'//'.t00z.'//'pgrb2af006' + print *,outfile006 + call baopenwa(200,outfile006,iret) ! for add more than 1 members +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + unit=10 + do ifid=1,nfield1 + nfi=1 + do ifh=003,006,3 !foercast hours + + nfi=nfi+1 + unit=unit+1 + write(sfh,'(i3.3)') ifh + +!search in pgrb2a + + datafile(nfi)=trim(datapath)//'gefs.t00z.master.grb2f'//trim(sfh) + + call BAOPENR(unit, datafile(nfi), iret) + if(iret /= 0 ) then + write(*,*) "there is no GEFS forecast",datafile(nfi) + end if + write(*,*) unit, datafile(nfi), iret + + j = 0 + jids=-9999;jpdt=-9999; jgdt=-9999 + jdisc=-1; jgdtn=-1 + jpdt(1) = pdt%icat(ifid) + jpdt(2) = pdt%iprm(ifid) + jpdt(10) = pdt%iffs(ifid) + jdisc = disc(ifid) + + if (nfi /= 3 ) then + jpdtn = pdt%npdt(ifid) !template version num. + else + jpdtn =11 !template version num. + endif + + call getgb2(unit,0,j,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt,unpack,j,gfld,iret) + if (iret /= 0) then + write(*,*) "reading file iret=", iret,"j=",j,"nfi=",nfi + stop + end if + call baclose(unit, iret) + write(*,*) 'after getgb2, iret=', iret + +! assign gfldo same arribute for gfld + + maxgrd=gfld%ngrdpts + print*,maxgrd + nx = gfld%igdtmpl(8) + ny = gfld%igdtmpl(9) +! assign gfldo same arribute for gfld + + gfldo=gfld + + print*,nfi,' ',ifh + print*,maxgrd,' ',maxgrd + + if(.not. allocated(apcp_3h)) allocate(apcp_3h(maxgrd)) + if(.not. allocated(acpcp_3h)) allocate(acpcp_3h(maxgrd)) + if(.not. allocated(apcp_6h)) allocate(apcp_6h(maxgrd)) + if(.not. allocated(acpcp_6h)) allocate(acpcp_6h(maxgrd)) + if(.not. allocated(dpcp)) allocate(dpcp(maxgrd)) + if(.not. allocated(var_save)) allocate(var_save(maxgrd,3)) + + var_save(:,nfi) = gfld%fld + year = gfld%idsect(6) + month = gfld%idsect(7) + day = gfld%idsect(8) + hour = gfld%idsect(9) + + write(syr,'(i4.4)') year + write(smonth,'(i2.2)') month + write(sdy,'(i2.2)') day + +!- forecast hour + fhour=gfld%ipdtmpl(9) + ens_id=gfld%ipdtmpl(17) + pabbrev=param_get_abbrev(gfld%discipline,gfld%ipdtmpl(1),gfld%ipdtmpl(2)) + call prlevel(gfld%ipdtnum,gfld%ipdtmpl,labbrev) + if ( nfi /= 3) then + call gf_free(gfld) + endif + + 201 end do !ifh + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + labbrev_short=trim(labbrev) + print *,pabbrev,trim(labbrev_short(1:4)) + + if(.not. allocated(var_save_acc)) allocate(var_save_acc(maxgrd)) + if(.not. allocated(var_save_USWRF)) allocate(var_save_USWRF(maxgrd)) + if(.not. allocated(var_save_DSWRF)) allocate(var_save_DSWRF(maxgrd)) + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!calculate the 00-6h ave/acc + if(pabbrev=='PRATE' .or. pabbrev=='CPRAT') then + var_save_acc(:)=(var_save(:,2)+var_save(:,3)*6./3.)/2. + + if(pabbrev=='PRATE') then + apcp_6h(:)=var_save_acc(:)*6.*3600. + else + acpcp_6h(:)=var_save_acc(:)*6.*3600. + endif + print*,'For 0-6h precip. and precip. rate' + print*,'' + + else +! var_save_acc(:)=(var_save(:,2)+var_save(:,3))/2. +! rate for 6, 12, 18,... hours is accumulated values/(current time - diagnosis time) + + var_save_acc(:)=(var_save(:,2)+var_save(:,3)*6./3.)/2. + + if (pabbrev=='DSWRF' ) then + var_save_DSWRF(:)=var_save_acc(:) + endif + if (ifid==4 ) then + var_save_USWRF(:)=var_save_acc(:) + endif + + if (pabbrev=='PRES' .or. pabbrev=='TMP') then + do i=1,maxgrd + if (var_save(i,2).eq.0.0.and.var_save(i,3).eq.0.0) var_save_acc(i)=var_save(i,3) + if (var_save(i,2).eq.0.0.and.var_save(i,3).gt.0.0) var_save_acc(i)=var_save(i,3) + if (var_save(i,2).gt.0.0.and.var_save(i,3).eq.0.0) var_save_acc(i)=var_save(i,2) + if (var_save(i,2).gt.0.0.and.var_save(i,3).gt.0.0) var_save_acc(i)=(var_save(i,2)+var_save(i,3))/2. + enddo + endif + + if (pabbrev=='ALBDO') then + do i=1,maxgrd + var_save_acc(i)=0.0 + if (var_save_DSWRF(i) .gt. 0.01 .and. var_save_USWRF(i).gt.0.01) then + var_save_acc(i)=100.*var_save_USWRF(i)/var_save_DSWRF(i) + if(var_save_acc(i) .gt. 100) var_save_acc(i)=0.0 + endif + enddo + endif + + if(pabbrev=='WATR' .or. pabbrev=='TSNOWP') then + do i=1,maxgrd + var_save_acc(i)=var_save(i,2)+var_save(i,3) + enddo + endif + + if(pabbrev=='TMIN') then + do i =1,maxgrd + if(var_save(i,3).gt.var_save(i,2)) then + var_save_acc(i)=var_save(i,2) + else + var_save_acc(i)=var_save(i,3) + endif + enddo + endif + + if(pabbrev=='TMAX') then + do i =1,maxgrd + if(var_save(i,3).gt.var_save(i,2)) then + var_save_acc(i)=var_save(i,3) + else + var_save_acc(i)=var_save(i,2) + endif + enddo + endif + endif + +!output 6h variables + + gfld%fld=var_save_acc(:) + gfld%ipdtmpl(22)=6 !forecast time + gfld%ipdtmpl(30)=6 !forecast time + jret=0 + +! call putgb2(200,gfldo,jret) + call putgb2(200,gfld,jret) + write(*,*) '200 put',jret,gfld%discipline,gfld%ipdtmpl + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + write(*,'(a10,a20,4f15.2)') pabbrev,trim(labbrev),var_save(100,2),var_save(100,2),var_save(100,3),var_save_acc(100) + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!calculate the 3h ave/acc +!3h accumulation 00-03z + if(pabbrev=='PRATE' .or. pabbrev=='CPRAT') then + + var_save_acc(:)=var_save(:,2) + + print*,'For 0-3h precip. rate' + print*,'' + + if(pabbrev=='PRATE') then + apcp_3h(:)=var_save(:,2)*3600.*3 + else + acpcp_3h(:)=var_save(:,2)*3600.*3 + endif + print*,'For 0-3h precip.' + print*,'' + + else + + var_save_acc(:)=var_save(:,2) + + end if + +!output 3-h variables + + gfld%fld=var_save_acc(:) +!reassign the ipdtmpl + gfld%ipdtmpl(9)=0 !forecast time + gfld%ipdtmpl(22)=3 !forecast time + gfld%ipdtmpl(30)=3 !forecast time + jret=0 + call putgb2(300,gfld,jret) + + write(*,*) '300 put',jret,gfldo%ipdtmpl +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + call gf_free(gfldo) + + deallocate(var_save_acc) + deallocate(var_save) + if(ifid.eq.38) then + deallocate(var_save_DSWRF) + deallocate(var_save_USWRF) + endif + enddo + + +! output accumulated precipitations + + unit=100 + call BAOPENR(unit, datafile(3), iret) + if(iret /= 0 ) then + write(*,*) "there is no GEFS forecast",datafile(nfi) + end if + + do ifid=40,42 + j = 0 + jids=-9999 + jids=-9999;jpdt=-9999; jgdt=-9999 + jdisc=-1; jgdtn=-1 + + jpdtn = pdt%npdt(ifid) !template version num. + jpdtn =11 !template version num. + jpdt(1) = pdt%icat(ifid) + jpdt(2) = pdt%iprm(ifid) + jpdt(10) = pdt%iffs(ifid) +! + write(*,*) unit,datafile(3),jpdt(1),jpdt(2),jpdt(10) + call getgb2(unit,0,j,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt,unpack,j,gfld,iret) + if (iret /= 0) then + write(*,*) "reading file iret=", iret + stop + end if + + gfldo=gfld + fhour=gfld%ipdtmpl(9) + ens_id=gfld%ipdtmpl(17) + pabbrev=param_get_abbrev(gfld%discipline,gfld%ipdtmpl(1),gfld%ipdtmpl(2)) + call prlevel(gfld%ipdtnum,gfld%ipdtmpl,labbrev) + + + +! assign gfldo same arribute for gfld + if(ifid.eq.40) then + + gfld%ipdtmpl(22)=6 !forecast time + gfld%ipdtmpl(30)=6 !forecast time + gfld%fld=apcp_6h(:) + jret=0 + call putgb2(200,gfld,jret) + write(*,*) '0-6h apcp put',jret + + gfld%fld=apcp_3h(:) + gfld%ipdtmpl(22)=3 !forecast time + gfld%ipdtmpl(30)=3 !forecast time + jret=0 + call putgb2(300,gfld,jret) + write(*,*) '0-3h apcp put',jret + + endif + + if(ifid.eq.41) then + + gfld%ipdtmpl(22)=6 !forecast time + gfld%ipdtmpl(30)=6 !forecast time + gfld%fld=acpcp_6h(:) + jret=0 + call putgb2(200,gfld,jret) + write(*,*) '0-6h acpcp put',jret + + gfld%ipdtmpl(22)=3 !forecast time + gfld%ipdtmpl(30)=3 !forecast time + gfld%fld=acpcp_3h(:) + jret=0 + call putgb2(300,gfld,jret) + write(*,*) '0-3h acpcp put',jret + + endif + + if(ifid.eq.42) then + gfld%ipdtmpl(22)=6 !forecast time + gfld%ipdtmpl(30)=6 !forecast time + gfld%fld=apcp_6h(:)-acpcp_6h(:) + do i=1,maxgrd + if(gfld%fld(i).lt.0.0) then + gfld%fld(i)=0.0 + endif + enddo + jret=0 + call putgb2(200,gfld,jret) + write(*,*) '0-6h ncpcp put',jret + + gfld%ipdtmpl(22)=3 !forecast time + gfld%ipdtmpl(30)=3 !forecast time + gfld%fld=apcp_3h(:)-acpcp_3h(:) + do i=1,maxgrd + if(gfld%fld(i).lt.0.0) then + gfld%fld(i)=0.0 + endif + enddo + jret=0 + call putgb2(300,gfld,jret) + write(*,*) '0-3h ncpcp put',jret + + + endif + + enddo + + deallocate(apcp_6h) + deallocate(acpcp_6h) + deallocate(apcp_3h) + deallocate(acpcp_3h) + deallocate(dpcp) + + call baclose(100,iret) + call baclose(200,iret) + call baclose(300,iret) + + end program gefs_6h_ave_1mem diff --git a/src/gefs_6h_ave_1mem.fd/printinfr.f90 b/src/gefs_6h_ave_1mem.fd/printinfr.f90 new file mode 100755 index 00000000..877da8df --- /dev/null +++ b/src/gefs_6h_ave_1mem.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/gefs_6h_ave_1mem.fd/prlevel.f90 b/src/gefs_6h_ave_1mem.fd/prlevel.f90 new file mode 100755 index 00000000..f3b17d5f --- /dev/null +++ b/src/gefs_6h_ave_1mem.fd/prlevel.f90 @@ -0,0 +1,272 @@ +subroutine prlevel(ipdtn,ipdtmpl,labbrev) + + integer,intent(in) :: ipdtn + integer,intent(in) :: ipdtmpl(*) + character(len=40),intent(out) :: labbrev + + character(len=10) :: tmpval1,tmpval2 + + labbrev(1:40)=" " + + selectcase (ipdtn) + case(0:15) + ipos=10 + case(40:43) + ipos=11 + case(44:47) + ipos=16 + case(48) + ipos=21 + case(50:51) + ipos=10 + case(52) + ipos=13 + case(91) + ipos=10 + case default + ipos=10 + end select + + if ( ipdtmpl(ipos) .eq. 100 .and.ipdtmpl(ipos+3) .eq. 255 ) then ! Pressure Level + !write(tmpval1,*) ipdtmpl(ipos+2)/100 + call frmt(tmpval1,ipdtmpl(ipos+2),ipdtmpl(ipos+1)+2) + labbrev=trim(tmpval1)//" mb" + + elseif ( ipdtmpl(ipos) .eq. 100 .and. ipdtmpl(ipos+3) .eq. 100 ) then ! Pressure Layer + call frmt(tmpval1,ipdtmpl(ipos+2),ipdtmpl(ipos+1)+2) + call frmt(tmpval2,ipdtmpl(ipos+5),ipdtmpl(ipos+4)+2) + labbrev=trim(tmpval1)//" - "//trim(tmpval2)//" mb" + + elseif ( ipdtmpl(ipos) .eq. 101 ) then ! Mean Sea Level + labbrev(1:30)=" Mean Sea Level " + + elseif ( ipdtmpl(ipos) .eq. 102 .and. ipdtmpl(ipos+3) .eq. 255 ) then ! Altitude above MSL + call frmt(tmpval1,ipdtmpl(ipos+2),ipdtmpl(ipos+1)) + labbrev=trim(tmpval1)//" m above MSL" + + elseif ( ipdtmpl(ipos) .eq. 103 .and. ipdtmpl(ipos+3) .eq. 255 ) then ! Height above Ground + call frmt(tmpval1,ipdtmpl(ipos+2),ipdtmpl(ipos+1)) + labbrev=trim(tmpval1)//" m above ground" + + elseif ( ipdtmpl(ipos) .eq. 103 .and. ipdtmpl(ipos+3) .eq. 103 ) then ! Height above Ground + call frmt(tmpval1,ipdtmpl(ipos+2),ipdtmpl(ipos+1)) + call frmt(tmpval2,ipdtmpl(ipos+5),ipdtmpl(ipos+4)) + labbrev=trim(tmpval1)//" - "//trim(tmpval2)//" m AGL" + + elseif ( ipdtmpl(ipos) .eq. 104 .and. ipdtmpl(ipos+3) .eq. 255 ) then ! Sigma Level + call frmt(tmpval1,ipdtmpl(ipos+2),ipdtmpl(ipos+1)) + labbrev=trim(tmpval1)//" sigma" + + elseif ( ipdtmpl(ipos) .eq. 104 .and. ipdtmpl(ipos+3) .eq. 104 ) then ! Sigma Layer + call frmt(tmpval1,ipdtmpl(ipos+2),ipdtmpl(ipos+1)) + call frmt(tmpval2,ipdtmpl(ipos+5),ipdtmpl(ipos+4)) + labbrev=trim(tmpval1)//" - "//trim(tmpval2)//" sigma" + + elseif ( ipdtmpl(ipos) .eq. 105 .and. ipdtmpl(ipos+3) .eq. 255 ) then ! Hybrid Level + call frmt(tmpval1,ipdtmpl(ipos+2),ipdtmpl(ipos+1)) + labbrev=trim(tmpval1)//" hybrid lvl" + + elseif ( ipdtmpl(ipos).eq.105 .and. ipdtmpl(ipos+3).eq.105) then ! Hybrid Level + call frmt(tmpval1,ipdtmpl(ipos+2),ipdtmpl(ipos+1)) + call frmt(tmpval2,ipdtmpl(ipos+5),ipdtmpl(ipos+4)) + labbrev=trim(tmpval1)//" - "//trim(tmpval2)//" hybrid lvl" + + + elseif ( ipdtmpl(ipos) .eq. 106 .and. ipdtmpl(ipos+3) .eq. 255 ) then ! Depth Below Land Sfc + call frmt(tmpval1,ipdtmpl(ipos+2),ipdtmpl(ipos+1)) + labbrev=trim(tmpval1)//" m below land" + + elseif ( ipdtmpl(ipos).eq.106 .and. ipdtmpl(ipos+3).eq.106) then ! Depth Below Land Sfc Layer + call frmt(tmpval1,ipdtmpl(ipos+2),ipdtmpl(ipos+1)) + call frmt(tmpval2,ipdtmpl(ipos+5),ipdtmpl(ipos+4)) + labbrev=trim(tmpval1)//" - "//trim(tmpval2)//" m DBLY" + + elseif ( ipdtmpl(ipos) .eq. 107 ) then ! Isentrophic (theta) level (THEL) + labbrev(1:30)=" Isentropic level" + + elseif ( ipdtmpl(ipos).eq.108 .and. ipdtmpl(ipos+3).eq.108) then ! Press Diff from Ground Layer + !write(tmpval1,*) ipdtmpl(ipos+2)/100 + !write(tmpval2,*) ipdtmpl(ipos+5)/100 + call frmt(tmpval1,ipdtmpl(ipos+2),ipdtmpl(ipos+1)+2) + call frmt(tmpval2,ipdtmpl(ipos+5),ipdtmpl(ipos+4)+2) + labbrev=trim(tmpval1)//" - "//trim(tmpval2)//" mb SPDY" + + elseif ( ipdtmpl(ipos) .eq. 110 ) then ! Layer between two hybrid levels (HYBY) + labbrev(1:30)=" Layer bet 2-hyb lvl" + + elseif ( ipdtmpl(ipos).eq.109 .and. ipdtmpl(ipos+3).eq.255) then ! Potential Vorticity Sfc + !write(tmpval1,*) ipdtmpl(ipos+2) + call frmt(tmpval1,ipdtmpl(ipos+2),ipdtmpl(ipos+1)-6) + labbrev=trim(tmpval1)//" pv surface" + + elseif ( ipdtmpl(ipos) .eq. 111 ) then ! Eta Level (EtaL) + labbrev(1:30)=" Eta level" + elseif ( ipdtmpl(ipos) .eq. 114 ) then ! Layer between two isentropic levels (THEY) + labbrev(1:30)=" Layer bet. 2-isent." + elseif ( ipdtmpl(ipos) .eq. 117 ) then ! Mixed layer depth + labbrev(1:30)=" Mixed layer depth" + elseif ( ipdtmpl(ipos) .eq. 120 ) then ! Layer between two Eta levels (EtaY) + labbrev(1:30)=" Layer bet. 2-Eta lvl" + elseif ( ipdtmpl(ipos) .eq. 121 ) then ! Layer between two isobaric surface (IBYH) + labbrev(1:30)=" Layer bet. 2-isob." + elseif ( ipdtmpl(ipos) .eq. 125 ) then ! Specified height level above ground (HGLH) + labbrev(1:30)=" Specified height lvl" + elseif ( ipdtmpl(ipos) .eq. 126 ) then ! Isobaric Level (ISBP) + labbrev(1:30)=" Isobaric level" + + elseif ( ipdtmpl(ipos) .eq. 160 ) then ! Depth below sea level + labbrev(1:30)=" Depth below sea lvl" + + elseif ( ipdtmpl(ipos) .eq. 170 ) then ! Ionospheric D-region + labbrev(1:30)=" Ionospheric D-region lvl" + + elseif ( ipdtmpl(ipos) .eq. 1 ) then ! Ground/Water Surface + labbrev(1:30)=" Surface " + + elseif ( ipdtmpl(ipos) .eq. 2 ) then ! Cloud base level (CBL) + labbrev(1:30)=" Cloud base lvl" + elseif ( ipdtmpl(ipos) .eq. 3 ) then ! Cloud top level (CTL) + labbrev(1:30)=" Cloud top lvl" + + elseif ( ipdtmpl(ipos) .eq. 4 ) then ! Freezing Level + labbrev(1:30)=" 0 Deg Isotherm" + + elseif ( ipdtmpl(ipos) .eq. 5 ) then ! Level of adiabatic condensation lifted + labbrev(1:30)=" Level of adiabatic" ! from the surface + + elseif ( ipdtmpl(ipos) .eq. 6 ) then ! Max Wind Level + labbrev(1:30)=" Max wind lvl" + + elseif ( ipdtmpl(ipos) .eq. 7 ) then ! Tropopause + labbrev(1:30)=" Tropopause" + + elseif ( ipdtmpl(ipos) .eq. 8 ) then ! Nominal top of Atmosphere + labbrev(1:30)=" Nom. top" + + elseif ( ipdtmpl(ipos) .eq. 9 ) then ! Sea bottom + labbrev(1:30)=" Sea Bottom" + + elseif ( ipdtmpl(ipos) .eq. 10 ) then ! Entire Atmosphere (EATM) + labbrev(1:30)=" Entire Atmosphere" + + elseif ( ipdtmpl(ipos) .eq. 11 ) then ! Cumulonimbus Base + labbrev(1:30)=" Cumulonimbus Base" + + elseif ( ipdtmpl(ipos) .eq. 12 ) then ! Cumulonimbus Top + labbrev(1:30)=" Cumulonimbus Top" + + elseif ( ipdtmpl(ipos) .eq. 20 ) then ! Isothermal level + labbrev(1:30)=" Isothermal level" + + elseif ( ipdtmpl(ipos) .eq. 200 ) then ! Entire Atmosphere (EATM) + labbrev(1:30)=" Entire Atmosphere" + elseif ( ipdtmpl(ipos) .eq. 201 ) then ! Entire ocean (EOCN) + labbrev(1:30)=" Entire ocean" + elseif ( ipdtmpl(ipos) .eq. 204 ) then ! Highest tropospheric freezing level (HTFL) + labbrev(1:30)=" Highest Frz. lvl" + elseif ( ipdtmpl(ipos) .eq. 206 ) then ! Grid scale cloud bottom level (GCBL) + labbrev(1:30)=" Grid scale cloud bl" + elseif ( ipdtmpl(ipos) .eq. 207 ) then ! Grid scale cloud top level (GCTL) + labbrev(1:30)=" Grid scale cloud tl" + elseif ( ipdtmpl(ipos) .eq. 209 ) then ! Boundary layer cloud bottom level (BCBL) + labbrev(1:30)=" Boundary layer cbl" + elseif ( ipdtmpl(ipos) .eq. 210 ) then ! Boundary layer cloud top level (BCTL) + labbrev(1:30)=" Boundary layer ctl" + elseif ( ipdtmpl(ipos) .eq. 211 ) then ! Boundary layer cloud layer (BCY) + labbrev(1:30)=" Boundary layer cl" + elseif ( ipdtmpl(ipos) .eq. 212 ) then ! Low cloud bottom level (LCBL) + labbrev(1:30)=" Low cloud bot. lvl" + elseif ( ipdtmpl(ipos) .eq. 213 ) then ! Low cloud top level (LCTL) + labbrev(1:30)=" Low cloud top lvl" + elseif ( ipdtmpl(ipos) .eq. 214 ) then ! Low cloud layer (LCY) + labbrev(1:30)=" Low cloud layer" + elseif ( ipdtmpl(ipos) .eq. 215 ) then ! Cloud ceiling (CEIL) + labbrev(1:30)=" Cloud ceiling" + elseif ( ipdtmpl(ipos) .eq. 220 ) then ! Planetary Boundary Layer (PBLRI) + labbrev(1:30)=" Planetary boundary" + elseif ( ipdtmpl(ipos) .eq. 221 ) then ! Layer Between Two Hybrid Levels (HYBY) + labbrev(1:30)=" Layer 2 Hybrid lvl " + elseif ( ipdtmpl(ipos) .eq. 222 ) then ! Middle cloud bottom level (MCBL) + labbrev(1:30)=" Mid. cloud bot. lvl" + elseif ( ipdtmpl(ipos) .eq. 223 ) then ! Middle cloud top level (MCTL) + labbrev(1:30)=" Mid. cloud top lvl" + elseif ( ipdtmpl(ipos) .eq. 224 ) then ! Middle cloud layer (MCY) + labbrev(1:30)=" Middle cloud layer" + elseif ( ipdtmpl(ipos) .eq. 232 ) then ! High cloud bottom level (HCBL) + labbrev(1:30)=" High cloud bot. lvl" + elseif ( ipdtmpl(ipos) .eq. 233 ) then ! High cloud top level (HCTL) + labbrev(1:30)=" High cloud top lvl" + elseif ( ipdtmpl(ipos) .eq. 234 ) then ! High cloud layer (HCY) + labbrev(1:30)=" High cloud layer" + elseif ( ipdtmpl(ipos) .eq. 235 ) then ! Ocean isotherm level (OITL) + labbrev(1:30)=" Ocean Isotherm lvl" + elseif ( ipdtmpl(ipos) .eq. 236 ) then ! Layer between two depth below ocean sfc (OLYR) + labbrev(1:30)=" Layer 2-depth below" + elseif ( ipdtmpl(ipos) .eq. 237 ) then ! Bottom of Ocean mixed layer (OBML) + labbrev(1:30)=" Bot. Ocean mix. lyr" + elseif ( ipdtmpl(ipos) .eq. 238 ) then ! Bottom of Ocean iisothermal layer (OBIL) + labbrev(1:30)=" Bot. Ocean iso. lyr" + elseif ( ipdtmpl(ipos) .eq. 239 ) then ! Layer ocean surface and 26C ocean + labbrev(1:30)=" layer ocean sfc 26C" ! isothermal level (S26CY) + elseif ( ipdtmpl(ipos) .eq. 240 ) then ! Ocean Mixed Layer + labbrev(1:30)=" Ocean Mixed Layer" + elseif ( ipdtmpl(ipos) .eq. 241 ) then ! Ordered Sequence of Data + labbrev(1:30)=" Order Seq. Of Data" + elseif ( ipdtmpl(ipos) .eq. 242 ) then ! Convective cloud bottom level (CCBL) + labbrev(1:30)=" Con. cloud bot. lvl" + elseif ( ipdtmpl(ipos) .eq. 243 ) then ! Convective cloud top level (CCTL) + labbrev(1:30)=" Con. cloud top lvl" + elseif ( ipdtmpl(ipos) .eq. 244 ) then ! Convective cloud layer (CCY) + labbrev(1:30)=" Conv. cloud layer" + elseif ( ipdtmpl(ipos) .eq. 245 ) then ! Lowest level of the wet bulb zero (LLTW) + labbrev(1:30)=" Lowest lvl wet bulb" + elseif ( ipdtmpl(ipos) .eq. 246 ) then ! Maximum equiv. potential temp. level (MTHE) + labbrev(1:30)=" Max. equi. potential" + elseif ( ipdtmpl(ipos) .eq. 247 ) then ! Equilibrium level (EHLT) + labbrev(1:30)=" Equilibrium level" + elseif ( ipdtmpl(ipos) .eq. 248 ) then ! Shallow convective cloud bottom level (SCBL) + labbrev(1:30)=" Shallow con. cld bl" + elseif ( ipdtmpl(ipos) .eq. 249 ) then ! Shallow convective cloud top level (SCTL) + labbrev(1:30)=" Shallow con. cld tl" + elseif ( ipdtmpl(ipos) .eq. 251 ) then ! Deep convective cloud bottom level (DCBL) + labbrev(1:30)=" Deep conv. cld bl" + elseif ( ipdtmpl(ipos) .eq. 252 ) then ! Deep convective cloud top level (DCTL) + labbrev(1:30)=" Deep conv. cld tl" + elseif ( ipdtmpl(ipos) .eq. 253 ) then ! Lowest bottom level of supercooled + labbrev(1:30)=" Lowest bot. lvl sup" ! liquid water layer (LBLSW) + elseif ( ipdtmpl(ipos) .eq. 254 ) then ! Highest top level of supercooled + labbrev(1:30)=" highest top lvl sup" ! liquid water layer (HBLSW) + + else + write(labbrev,fmt='(1x,I4," (Unknown Lvl)")') ipdtmpl(ipos) + endif + + return + end + + + subroutine frmt(cval,ival,iscal) + + character(len=10),intent(out) :: cval + integer,intent(in) :: ival,iscal + + real :: rval + integer :: newscal + character(len=7) :: cformat + + if ( iscal .eq. 0 ) then + write(cval,'(I0)') ival + else + newscal=-1*iscal + rval=real(ival)*(10.0**newscal) + if ( rval .eq. real(nint(rval)) ) then + write(cval,'(1X,I0)') nint(rval) + else + write(cformat,fmt='("(f0.",I1,")")') iabs(iscal) + !write(6,'(A)') "DEGRIB2:",cformat + write(cval,fmt=cformat) rval + endif + endif + + return + end