UPP  001
 All Data Structures Files Functions Pages
grib2_module.f
1  module grib2_module
2 !------------------------------------------------------------------------
3 !
4 ! This module generates grib2 messages and writes out the messages in
5 ! parallel.
6 !
7 ! program log:
8 ! March, 2010 Jun Wang Initial code
9 ! Jan, 2012 Jun Wang post available fields with grib2 description
10 ! are defined in xml file
11 ! March, 2015 Lin Gan Replace XML file with flat file implementation
12 ! with parameter marshalling
13 !------------------------------------------------------------------------
14  use xml_perl_data, only: param_t,paramset_t
15 !
16  implicit none
17  private
18 ! ------------------------------------------------------------------------
19 !
20 !--- general grib2 info provided by post control file
21 ! type param_t
22 ! integer :: post_avblfldidx=-9999
23 ! character(len=80) :: shortname=''
24 ! character(len=300) :: longname=''
25 ! character(len=30) :: pdstmpl=''
26 ! integer :: mass_windpoint=1
27 ! character(len=30) :: pname=''
28 ! character(len=10) :: table_info=''
29 ! character(len=20) :: stats_proc=''
30 ! character(len=80) :: fixed_sfc1_type=''
31 ! integer, dimension(:), pointer :: scale_fact_fixed_sfc1 => null()
32 ! real, dimension(:), pointer :: level => null()
33 ! character(len=80) :: fixed_sfc2_type=''
34 ! integer, dimension(:), pointer :: scale_fact_fixed_sfc2 => null()
35 ! real, dimension(:), pointer :: level2 => null()
36 ! character(len=80) :: aerosol_type=''
37 ! character(len=80) :: typ_intvl_size=''
38 ! integer :: scale_fact_1st_size=0
39 ! real :: scale_val_1st_size=0.
40 ! integer :: scale_fact_2nd_size=0
41 ! real :: scale_val_2nd_size=0.
42 ! character(len=80) :: typ_intvl_wvlen=''
43 ! integer :: scale_fact_1st_wvlen=0
44 ! real :: scale_val_1st_wvlen=0.
45 ! integer :: scale_fact_2nd_wvlen=0
46 ! real :: scale_val_2nd_wvlen=0.
47 ! real, dimension(:), pointer :: scale => null()
48 ! integer :: stat_miss_val=0
49 ! integer :: leng_time_range_prev=0
50 ! integer :: time_inc_betwn_succ_fld=0
51 ! character(len=80) :: type_of_time_inc=''
52 ! character(len=20) :: stat_unit_time_key_succ=''
53 ! character(len=20) :: bit_map_flag=''
54 ! integer :: perturb_num=0
55 ! integer :: num_ens_fcst=10
56 ! end type param_t
57 !
58 ! type paramset_t
59 ! character(len=6) :: datset=''
60 ! integer :: grid_num=255
61 ! character(len=20) :: sub_center=''
62 ! character(len=20) :: version_no=''
63 ! character(len=20) :: local_table_vers_no=''
64 ! character(len=20) :: sigreftime=''
65 ! character(len=20) :: prod_status=''
66 ! character(len=20) :: data_type=''
67 ! character(len=20) :: gen_proc_type=''
68 ! character(len=30) :: time_range_unit=''
69 ! character(len=50) :: orig_center=''
70 ! character(len=30) :: gen_proc=''
71 ! character(len=20) :: packing_method=''
72 ! character(len=20) :: field_datatype=''
73 ! character(len=20) :: comprs_type=''
74 ! character(len=50) :: type_ens_fcst=''
75 ! character(len=50) :: type_derived_fcst=''
76 ! type(param_t), dimension(:), pointer :: param => null()
77 ! end type paramset_t
78  type(paramset_t),save :: pset
79 !
80 !--- grib2 info related to a specific data file
81  integer nrecout
82  integer num_pset
83  integer isec,hrs_obs_cutoff,min_obs_cutoff
84  integer sec_intvl,stat_miss_val,time_inc_betwn_succ_fld
85  integer perturb_num,num_ens_fcst
86  character*80 type_of_time_inc,stat_unit_time_key_succ
87  logical*1,allocatable :: bmap(:)
88  integer ibm
89  integer,allocatable :: mg(:)
90 !
91  integer,parameter :: max_bytes=1000*1300000
92  integer,parameter :: max_numbit=16
93  integer,parameter :: lugi=650
94  character*255 fl_nametbl,fl_gdss3
95  logical :: first_grbtbl
96 !
97  public num_pset,pset,nrecout,gribit2,grib_info_init,first_grbtbl,grib_info_finalize,read_grib2_head,read_grib2_sngle
98  real(8), EXTERNAL :: timef
99 !-------------------------------------------------------------------------------------
100 !
101  contains
102 !
103 !-------------------------------------------------------------------------------------
104  subroutine grib_info_init()
105 !
106 !--- initialize general grib2 information and
107 !
108  implicit none
109 !
110 ! logical,intent(in) :: first_grbtbl
111 !
112 !-- local variables
113  integer ierr
114  character(len=80) outfile
115  character(len=10) envvar
116 !
117 !-- set up pset
118 !
119 !-- 1. pset is set up at READCNTRL_xml.f
120 !-- initialize items of pset that are not set in xml file
121 !
122  if(pset%grid_num==0) &
123  pset%grid_num=218
124  if(trim(pset%sub_center)=='') &
125  pset%sub_center="ncep_emc"
126  if(trim(pset%version_no)=='') &
127  pset%version_no="v2003"
128  if(trim(pset%local_table_vers_no)=='') &
129  pset%local_table_vers_no="local_table_no"
130  if(trim(pset%sigreftime)=='') &
131  pset%sigreftime="fcst"
132  if(trim(pset%prod_status)=='') &
133  pset%prod_status="oper_test"
134  if(trim(pset%data_type)=='') &
135  pset%data_type="fcst"
136  if(trim(pset%orig_center)=='') &
137  pset%orig_center="nws_ncep"
138  if(trim(pset%time_range_unit)=='') &
139  pset%time_range_unit="hour"
140  if(trim(pset%gen_proc_type)=='') &
141  pset%gen_proc_type="fcst"
142  if(trim(pset%gen_proc)=='') &
143  pset%gen_proc="gfs_avn"
144  if(trim(pset%packing_method)=='') &
145  pset%packing_method="jpeg"
146  if(trim(pset%field_datatype)=='') &
147  pset%field_datatype="flting_pnt"
148  if(trim(pset%comprs_type)=='') &
149  pset%comprs_type="lossless"
150  if(trim(pset%type_ens_fcst)=='') &
151  pset%type_ens_fcst="pos_pert_fcst"
152  if(trim(pset%type_derived_fcst)=='') &
153  pset%type_derived_fcst="unweighted_mean_all_mem"
154 !
155 !-- set up other grib2_info
156 !
157  isec=0
158  hrs_obs_cutoff=0 ! applies to only obs
159  min_obs_cutoff=0 ! applies to only obs
160  sec_intvl=0
161  stat_miss_val=0
162  type_of_time_inc='same_start_time_fcst_fcst_time_inc'
163  stat_unit_time_key_succ='missing'
164  time_inc_betwn_succ_fld=0
165 !
166 !-- open fld name tble
167 !
168  if(first_grbtbl) then
169  fl_nametbl='params_grib2_tbl_new'
170  call open_and_read_4dot2( fl_nametbl, ierr )
171  if ( ierr /= 0 ) then
172  print*, 'Couldnt open table file - return code was ',ierr
173  call mpi_abort()
174  endif
175  first_grbtbl=.false.
176  endif
177 !
178 !--
179 !
180  end subroutine grib_info_init
181 !-------------------------------------------------------------------------------------
182 !-------------------------------------------------------------------------------------
183 !
184  subroutine grib_info_finalize
185 !
186 !--- finalize grib2 information and close file
187 !
188  implicit none
189 !
190 !---
191  integer ierr
192  call close_4dot2(ierr)
193 !
194  end subroutine grib_info_finalize
195 !-------------------------------------------------------------------------------------
196 !-------------------------------------------------------------------------------------
197  subroutine gribit2(post_fname)
198 !
199 !-------
200  use ctlblk_mod, only : im,jm,im_jm,num_procs,me,jsta,jend,ifhr,sdat,ihrst,imin, &
201  mpi_comm_comp,ntlfld,fld_info,datapd,icnt,idsp
202  implicit none
203 !
204  include 'mpif.h'
205 !
206 ! real,intent(in) :: data(im,1:jend-jsta+1,ntlfld)
207  character(255),intent(in) :: post_fname
208 !
209 !------- local variables
210  integer i,j,k,n,nm,nprm,nlvl,fldlvl1,fldlvl2,cstart,cgrblen,ierr
211  integer nf,nfpe,nmod
212  integer fh, clength,lunout
213  integer idisc,icatg,iparm,itblinfo,ntrange,leng_time_range_stat
214  integer,allocatable :: nfld_pe(:),snfld_pe(:),enfld_pe(:)
215  integer(4),allocatable :: isdsp(:),iscnt(:),ircnt(:),irdsp(:)
216  integer status(mpi_status_size)
217  integer(kind=MPI_OFFSET_KIND) idisp
218  integer,allocatable :: jsta_pe(:),jend_pe(:)
219  integer,allocatable :: grbmsglen(:)
220  real,allocatable :: datafld(:,:)
221  real,allocatable :: datafldtmp(:)
222  logical, parameter :: debugprint = .false.
223 !
224  character(1), dimension(:), allocatable :: cgrib
225 !
226 !
227 !---------------- code starts here --------------------------
228 !
229  allocate(cgrib(max_bytes))
230 !
231 !******* part 1 resitribute data ********
232 !
233 !--- calculate # of fields on each processor
234 !
235  nf=ntlfld/num_procs
236  nfpe=nf+1
237  nmod=mod(ntlfld,num_procs)
238 ! print *,'ntlfld=',ntlfld,'nf=',nf,'nmod=',nmod
239  allocate(snfld_pe(num_procs),enfld_pe(num_procs),nfld_pe(num_procs))
240  do n=1,num_procs
241  if(n-1<nmod ) then
242  snfld_pe(n)=nfpe*(n-1)+1
243  enfld_pe(n)=snfld_pe(n)+nfpe-1
244  nfld_pe(n)=nfpe
245  else
246  snfld_pe(n)=nfpe*nmod+nf*(n-1-nmod)+1
247  enfld_pe(n)=snfld_pe(n)+nf-1
248  nfld_pe(n)=nf
249  endif
250  enddo
251 ! print *,'in gribit2,ntlfld=',ntlfld,'nf=',nf,'myfld=',snfld_pe(me+1),enfld_pe(me+1)
252 !
253 !--- reditribute data from partial domain data with all fields
254 !--- to whole domain data but partial fields
255 !
256  allocate(jsta_pe(num_procs),jend_pe(num_procs))
257  call mpi_allgather(jsta,1,mpi_integer,jsta_pe,1, &
258  mpi_integer,mpi_comm_comp,ierr)
259  call mpi_allgather(jend,1,mpi_integer,jend_pe,1, &
260  mpi_integer,mpi_comm_comp,ierr)
261 ! print *,'in gribit2,jsta_pe=',jsta_pe,'jend_pe=',jend_pe
262 !
263 !--- end part1
264 !
265 !********************* generate grib2 message and write out data ****
266 !
267  allocate(bmap(im_jm))
268  allocate(mg(im_jm))
269 !
270 !--- sequatial write if the number of fields to write is small
271 !
272  if(minval(nfld_pe)<1.or.num_procs==1) then
273 !
274 !-- collect data to pe 0
275  allocate(datafld(im_jm,ntlfld) )
276  if(num_procs==1) then
277  datafld=reshape(datapd,(/im_jm,ntlfld/))
278  else
279  do i=1,ntlfld
280  call mpi_gatherv(datapd(:,:,i),icnt(me),mpi_real, &
281  datafld(:,i),icnt,idsp,mpi_real,0,mpi_comm_comp,ierr)
282  enddo
283  endif
284 !
285 !-- pe 0 create grib2 message and write to the file
286  if(me==0) then
287 !
288  lunout=601
289  call baopenw(lunout,trim(post_fname),ierr)
290  print*,'write_grib2: opened ',lunout, &
291  'for grib2 data ',trim(post_fname), &
292  'return code is ',ierr
293 !
294  do i=1,ntlfld
295  nprm=fld_info(i)%ifld
296  nlvl=fld_info(i)%lvl
297  fldlvl1=fld_info(i+snfld_pe(me+1)-1)%lvl1
298  fldlvl2=fld_info(i+snfld_pe(me+1)-1)%lvl2
299  if(trim(pset%param(nprm)%table_info)=='NCEP') then
300  itblinfo=1
301  else
302  itblinfo=0
303  endif
304 ! print *,'i=',i,'nprm=',fld_info(i)%ifld,'pname=',trim(pset%param(nprm)%pname), &
305 ! 'lev_type=',trim(pset%param(nprm)%fixed_sfc1_type),'itblinfo=',itblinfo, &
306 ! 'nlvl=',nlvl,'lvl1=',fldlvl1,'lvl2=',fldlvl2, &
307 ! 'shortname=',trim(pset%param(nprm)%shortname)
308  call search_for_4dot2_entry( &
309  pset%param(nprm)%pname, &
310  itblinfo, &
311  idisc, icatg, iparm, ierr)
312  if(ierr==0) then
313  write(6,'(3(A,I4),A,A)') ' discipline ',idisc, &
314  ' category ',icatg, &
315  ' parameter ',iparm, &
316  ' for var ',trim(pset%param(nprm)%pname)
317 
318  call gengrb2msg(idisc,icatg, iparm,nprm,nlvl,fldlvl1,fldlvl2, &
319  fld_info(i)%ntrange,fld_info(i)%tinvstat,datafld(:,i), &
320  cgrib,clength)
321 ! print *,'finished gengrb2msg field=',i,'ntlfld=',ntlfld,'clength=',clength
322  call wryte(lunout, clength, cgrib)
323  else
324  print *,'WRONG, could not find ',trim(pset%param(nprm)%pname), &
325  " in WMO and NCEP table!, ierr=", ierr
326  call mpi_abort()
327  endif
328  enddo
329 !
330  call baclose(lunout,ierr)
331  print *,'finish one grib file'
332  endif ! if(me==0)
333 !
334 !for more fields, use parallel i/o
335  else
336 !
337 ! print *,'in grib2,num_procs=',num_procs
338  allocate(iscnt(num_procs),isdsp(num_procs))
339  allocate(ircnt(num_procs),irdsp(num_procs))
340  isdsp(1)=0
341  do n=1,num_procs
342  iscnt(n)=(jend_pe(me+1)-jsta_pe(me+1)+1)*im*nfld_pe(n)
343  if(n<num_procs)isdsp(n+1)=isdsp(n)+iscnt(n)
344  enddo
345 !
346  irdsp(1)=0
347  do n=1,num_procs
348  ircnt(n)=(jend_pe(n)-jsta_pe(n)+1)*im*nfld_pe(me+1)
349  if(n<num_procs)irdsp(n+1)=irdsp(n)+ircnt(n)
350  enddo
351 ! print *,'in grib2,iscnt=',iscnt(1:num_procs),'ircnt=',ircnt(1:num_procs), &
352 ! 'nfld_pe=',nfld_pe(me+1)
353  allocate(datafldtmp(im_jm*nfld_pe(me+1)) )
354  allocate(datafld(im_jm,nfld_pe(me+1)) )
355 !
356  call mpi_alltoallv(datapd,iscnt,isdsp,mpi_real, &
357  datafldtmp,ircnt,irdsp,mpi_real,mpi_comm_comp,ierr)
358 !
359 !--- re-arrange the data
360  datafld=0.
361  nm=0
362  do n=1,num_procs
363  do k=1,nfld_pe(me+1)
364  do j=jsta_pe(n),jend_pe(n)
365  do i=1,im
366  nm=nm+1
367  datafld((j-1)*im+i,k)=datafldtmp(nm)
368  enddo
369  enddo
370  enddo
371  enddo
372  deallocate(datafldtmp)
373 !
374 !-- now each process has several full domain fields, start to create grib2 message.
375 !
376 ! print *,'nfld',nfld_pe(me+1),'snfld=',snfld_pe(me+1)
377 ! print *,'nprm=', &
378 ! fld_info(snfld_pe(me+1):snfld_pe(me+1)+nfld_pe(me+1)-1)%ifld
379 ! print *,'pname=',pset%param(5)%pname
380  cstart=1
381  do i=1,nfld_pe(me+1)
382  nprm=fld_info(i+snfld_pe(me+1)-1)%ifld
383  nlvl=fld_info(i+snfld_pe(me+1)-1)%lvl
384  fldlvl1=fld_info(i+snfld_pe(me+1)-1)%lvl1
385  fldlvl2=fld_info(i+snfld_pe(me+1)-1)%lvl2
386  ntrange=fld_info(i+snfld_pe(me+1)-1)%ntrange
387  leng_time_range_stat=fld_info(i+snfld_pe(me+1)-1)%tinvstat
388  if(trim(pset%param(nprm)%table_info)=='NCEP') then
389  itblinfo=1
390  else
391  itblinfo=0
392  endif
393 ! print *,'i=',i,'nprm=',nprm,'pname=',trim(pset%param(nprm)%pname), &
394 ! 'lev_type=',trim(pset%param(nprm)%fixed_sfc1_type),'itblinfo=',itblinfo, &
395 ! 'nlvl=',nlvl,'ntrange=',ntrange,'leng_time_range_stat=', &
396 ! leng_time_range_stat,'fldlvl1=',fldlvl1,'fldlvl2=',fldlvl2,'cfld=',i+snfld_pe(me+1)-1
397  call search_for_4dot2_entry( &
398  pset%param(nprm)%pname, &
399  itblinfo, &
400  idisc, icatg, iparm, ierr)
401  if(ierr==0) then
402  if(debugprint) then
403  write(6,'(3(A,I4),A,A)') ' discipline ',idisc, &
404  ' category ',icatg, &
405  ' parameter ',iparm, &
406  ' for var ',trim(pset%param(nprm)%pname)
407  endif
408 !
409 !--- generate grib2 message ---
410 !
411  call gengrb2msg(idisc,icatg, iparm,nprm,nlvl,fldlvl1,fldlvl2,ntrange, &
412  leng_time_range_stat,datafld(:,i),cgrib(cstart),clength)
413  cstart=cstart+clength
414 !
415  else
416  if(debugprint) then
417  print *,'WRONG, could not find ',trim(pset%param(nprm)%pname), &
418  " in WMO and NCEP table!"
419  endif
420 !!! call mpi_abort()
421  endif
422 !
423  enddo
424  cgrblen=cstart-1
425 ! print *,'after collect all data,cgrblen=',cgrblen
426 !
427 !******* write out grib2 message using MPI I/O *******************
428 !
429 !--- open file that will store grib2 messages
430 !
431  call mpi_barrier(mpi_comm_comp,ierr)
432 !
433 ! print *,'bf mpi_file_open,fname=',trim(post_fname)
434  call mpi_file_open(mpi_comm_comp,trim(post_fname), &
435  mpi_mode_create+mpi_mode_wronly,mpi_info_null,fh,ierr)
436 ! print *,'af mpi_file_open,ierr=',ierr
437 !
438 !--- broadcast message size
439  allocate(grbmsglen(num_procs))
440  call mpi_allgather(cgrblen,1,mpi_integer,grbmsglen,1,mpi_integer, &
441  mpi_comm_comp,ierr)
442 ! print *,'after gather gribmsg length=',grbmsglen(1:num_procs)
443 !
444 !--- setup start point
445  idisp=0
446  do n=1,me
447  idisp=idisp+grbmsglen(n)
448  enddo
449 !
450  call mpi_file_write_at(fh,idisp,cgrib,cgrblen,mpi_character,status,ierr)
451 !
452  call mpi_file_close(fh,ierr)
453 !
454 !--- deallocate arrays
455 !
456  deallocate(grbmsglen)
457  deallocate(iscnt,isdsp,ircnt,irdsp)
458 !
459  endif
460 !
461  deallocate(datafld,bmap,mg)
462  deallocate(nfld_pe,snfld_pe,enfld_pe,jsta_pe,jend_pe)
463  deallocate(cgrib)
464 !
465  end subroutine gribit2
466 !
467 !----------------------------------------------------------------------------------------
468 !----------------------------------------------------------------------------------------
469 !
470  subroutine gengrb2msg(idisc,icatg, iparm,nprm,nlvl,fldlvl1,fldlvl2,ntrange,tinvstat, &
471  datafld1,cgrib,lengrib)
472 !
473 !----------------------------------------------------------------------------------------
474 !
475  use ctlblk_mod, only : im,jm,im_jm,ifhr,idat,sdat,ihrst,ifmin,imin,fld_info,spval, &
476  vtimeunits,modelname
477  use gridspec_mod, only: maptype
478  use grib2_all_tables_module, only: g2sec0,g2sec1, &
479  g2sec4_temp0,g2sec4_temp8,g2sec4_temp44,g2sec4_temp48, &
480  g2sec5_temp0,g2sec5_temp2,g2sec5_temp3,g2sec5_temp40, &
481  get_g2_sec5packingmethod
482  !use gdtsec3, only: getgdtnum
483  implicit none
484 !
485  integer,intent(in) :: idisc,icatg, iparm,nprm,fldlvl1,fldlvl2,ntrange,tinvstat
486  integer,intent(inout) :: nlvl
487  real,dimension(:),intent(in) :: datafld1
488  character(1),intent(inout) :: cgrib(max_bytes)
489  integer, intent(inout) :: lengrib
490 !
491  integer, parameter :: igdsmaxlen=200
492 !
493  integer, parameter :: ipdstmplenmax=100
494  integer, parameter :: ipdstmp4_0len=15
495  integer, parameter :: ipdstmp4_1len=18
496  integer, parameter :: ipdstmp4_8len=29
497  integer, parameter :: ipdstmp4_11len=32
498  integer, parameter :: ipdstmp4_12len=31
499  integer, parameter :: ipdstmp4_44len=21
500  integer, parameter :: ipdstmp4_48len=26
501 !
502  integer, parameter :: idrstmplenmax=50
503  integer, parameter :: idrstmp5_0len=5
504  integer, parameter :: idrstmp5_2len=16
505  integer, parameter :: idrstmp5_3len=18
506  integer, parameter :: idrstmp5_40len=7
507 !
508  integer :: mxbit
509  integer listsec0(2) ! Length of section 0 octets 7 & 8
510  integer listsec1(13) ! Length of section 1 from octets 6-21
511  integer ipdstmpllen ! Length of general Section 4 PDS Template
512  integer ipdstmpl(ipdstmplenmax) ! Length of Section 4 PDS Template 4.48
513  integer idrstmplen
514  integer idrstmpl(idrstmplenmax) ! Length of Section 5 PDS Template 5.40
515  integer igds(5) ! Length of section 3 GDS Octet 6-14
516  integer igdstmplen
517  integer igdtlen,igdtn
518  integer idefnum
519  integer ideflist(100)
520  integer idrsnum,numcoord,ipdsnum
521  integer scaled_val_fixed_sfc2,scale_fct_fixed_sfc1
522  integer scaled_val_fixed_sfc1,scale_fct_fixed_sfc2
523  character(80) fixed_sfc2_type
524  integer idec_scl,ibin_scl,ibmap,inumbits
525  real fldscl
526  integer igdstmpl(igdsmaxlen)
527  integer lat1,lon1,lat2,lon2,lad,ds1
528  real(4) coordlist(1)
529  logical ldfgrd
530 !
531  integer ierr,ifhrorig,ihr_start
532  integer gefs1,gefs2,gefs3,gefs_status
533  character(len=4) cdum
534  integer perturb_num,num_ens_fcst,e1_type
535 !
536 !----------------------------------------------------------------------------------------
537 ! Find out if the Post is being run for the GEFS model
538 ! Check if gen_proc is gefs
539  gefs_status=0
540  if(trim(pset%gen_proc)=='gefs') then
541  call getenv('e1',cdum)
542  read(cdum,'(I4)',iostat=gefs_status)gefs1
543  e1_type=gefs1
544 
545  if(gefs_status /= 0) print *, &
546  "GEFS Run: Could not read e1 envir. var, User needs to set in script"
547 
548  call getenv('e2',cdum)
549  read(cdum,'(I4)',iostat=gefs_status)gefs2
550  perturb_num=gefs2
551 
552  if(gefs_status /= 0) print *, &
553  "GEFS Run: Could not read e2 envir. var, User needs to set in script"
554 
555  !set default number of ens forecasts to 10 for GEFS
556  !num_ens_fcst=10
557  call getenv('e3',cdum)
558  read(cdum,'(I4)',iostat=gefs_status)gefs3
559  num_ens_fcst=gefs3
560 
561  if(gefs_status /= 0) print *, &
562  "GEFS Run: Could not read e3 envir. var, User needs to set in script"
563 
564  print*,'GEFS env var ',e1_type,perturb_num,num_ens_fcst
565 
566  ! Set pdstmpl to tmpl4_1 or tmpl4_11
567  print *, "Processing for GEFS and default setting is tmpl4_1 and tmpl4_11"
568  if (trim(pset%param(nprm)%pdstmpl)=='tmpl4_0') then
569  pset%param(nprm)%pdstmpl='tmpl4_1'
570  elseif (trim(pset%param(nprm)%pdstmpl)=='tmpl4_8') then
571  pset%param(nprm)%pdstmpl='tmpl4_11'
572  endif
573  endif
574 !
575 !----------------------------------------------------------------------------------------
576 ! Feed input keys for GRIB2 Section 0 and 1 and get outputs from arrays listsec0 and listsec1
577 !
578  call g2sec0(idisc,listsec0)
579 !
580 !----------------------------------------------------------------------------------------
581 !GRIB2 - SECTION 1
582 ! IDENTIFICATION SECTION
583 !Octet No. Content
584 !1-4 Length of the section in octets (21 or N)
585 !5 Number of the section (1)
586 !6-7 Identification of originating/generating center (See Table 0 {GRIB1}) ! keys_sec1(1)
587 !8-9 Identification of originating/generating subcenter (See Table C) ! keys_sec1(2)
588 !10 GRIB master tables version number (currently 2) (See Table 1.0) (See note 1 below) ! keys_sec1(3)
589 !11 Version number of GRIB local tables used to augment Master Tables (see Table 1.1) ! keys_sec1(4)
590 !12 Significance of reference time (See Table 1.2) ! keys_sec1(5)
591 !13-14 Year (4 digits) ! keys_sec1(6)
592 !15 Month ! keys_sec1(7)
593 !16 Day ! keys_sec1(8)
594 !17 Hour ! keys_sec1(9)
595 !18 Minute ! keys_sec1(10)
596 !19 Second ! keys_sec1(11)
597 !20 Production Status of Processed data in the GRIB message (See Table 1.3) ! keys_sec1(12)
598 !21 Type of processed data in this GRIB message (See Table 1.4) ! keys_sec1(13)
599 !22-N Reserved
600 !
601  call g2sec1(pset%orig_center,pset%sub_center, &
602  pset%version_no,pset%local_table_vers_no,&
603  pset%sigreftime,nint(sdat(3)),nint(sdat(1)),nint(sdat(2)),ihrst,imin, &
604  isec,pset%prod_status,pset%data_type,listsec1)
605 !jw : set sect1(2) to 0 to compare with cnvgrb grib file
606 ! For GEFS runs we need to set the section 1 values for Grib2
607  if(trim(pset%gen_proc)=='gefs') then
608  listsec1(2)=2
609 ! Settings below for control (1 or 2) vs perturbed (3 or 4) ensemble forecast
610  if(e1_type==1.or.e1_type==2) then
611  listsec1(13)=3
612  elseif(e1_type==3.or.e1_type==4) then
613  listsec1(13)=4
614  endif
615  print *, "After g2sec1 call we need to set listsec1(2) = ",listsec1(2)
616  print *, "After g2sec1 call we need to set listsec1(13) = ",listsec1(13)
617  else
618  listsec1(2)=0
619  endif
620 !
621  call gribcreate(cgrib,max_bytes,listsec0,listsec1,ierr)
622 !
623 !----------------------------------------------------------------------------------------
624 ! Packup Grid Definition Section (section 3) and add to GRIB2 message
625 !
626 ! Define all the above derived data types and the input values will be available
627 ! through fortran modules
628 !
629  igdtlen=19
630  ldfgrd=(maptype==203.and.(trim(pset%param(nprm)%pname)=='ugrd'.or. &
631  trim(pset%param(nprm)%pname)=='vgrd'))
632  call getgds(ldfgrd,igdsmaxlen,igdtlen,igds,igdstmpl)
633  idefnum=1
634  ideflist=0 !Used if igds(3) /= 0. Dummy array otherwise
635 !
636  call addgrid(cgrib,max_bytes,igds,igdstmpl,igdtlen,ideflist,idefnum,ierr)
637 !
638 !----------------------------------------------------------------------------------------
639 ! Packup sections 4 through 7 for a given field and add them to a GRIB2 message which are
640 ! Product Defintion Section, Data Representation Section, Bit-Map Section and Data Section
641 ! respectively
642 !
643 !GRIB2 - TEMPLATE 4.0
644 !Product definition template analysis or forecast at a horizontal level or in a
645 !horizontal layer at a point in time
646 !Revised 09/21/2007
647 !Octet Contents
648 !10 Parameter category (see Code table 4.1)
649 !11 Parameter number (see Code table 4.2)
650 !12 Type of generating process (see Code table 4.3)
651 !13 Background generating process identifier (defined by originating centre)
652 !14 Analysis or forecast generating process identified (see Code ON388 Table A)
653 !15-16 Hours of observational data cutoff after reference time (see Note)
654 !17 Minutes of observational data cutoff after reference time (see Note)
655 !18 Indicator of unit of time range (see Code table 4.4)
656 !19-22 Forecast time in units defined by octet 18
657 !23 Type of first fixed surface (see Code table 4.5)
658 !24 Scale factor of first fixed surface
659 !25-28 Scaled value of first fixed surface
660 !29 Type of second fixed surfaced (see Code table 4.5)
661 !30 Scale factor of second fixed surface
662 !31-34 Scaled value of second fixed surfaces
663 !Notes: Hours greater than 65534 will be coded as 65534
664 !
665 !PRODUCT TEMPLATE 4. 0 : 3 5 2 0 96 0 0 1 12 100 0 100 255 0 0
666 ! TEXT: HGT 1 mb valid at 12 hr after 2009110500:00:00
667 !
668  coordlist=0
669  numcoord=0
670 ! print *,'size(level)=',size(pset%param(nprm)%level),'nlvl=',nlvl, &
671 ! 'lev_type=',trim(pset%param(nprm)%fixed_sfc1_type),'fldlvl1=', &
672 ! fldlvl1,'fldlvl2=',fldlvl2
673 !lvl is shown in ctl file
674  if(fldlvl1==0.and.fldlvl2==0) then
675 
676  if(size(pset%param(nprm)%level)>1.and.size(pset%param(nprm)%level)>=nlvl) then
677  scaled_val_fixed_sfc1=nint(pset%param(nprm)%level(nlvl))
678  else if(size(pset%param(nprm)%level)==1) then
679  scaled_val_fixed_sfc1=nint(pset%param(nprm)%level(1))
680  else
681  scaled_val_fixed_sfc1=0
682  endif
683  if(size(pset%param(nprm)%scale_fact_fixed_sfc1)>1.and. &
684  size(pset%param(nprm)%scale_fact_fixed_sfc1)>=nlvl) then
685  scale_fct_fixed_sfc1=pset%param(nprm)%scale_fact_fixed_sfc1(nlvl)
686  else if(size(pset%param(nprm)%scale_fact_fixed_sfc1)==1) then
687  scale_fct_fixed_sfc1=pset%param(nprm)%scale_fact_fixed_sfc1(1)
688  else
689  scale_fct_fixed_sfc1=0
690  endif
691 !
692 !for hygrid dpes level is decided in the code, not from ctl file
693  else
694  scaled_val_fixed_sfc1=fldlvl1
695  scale_fct_fixed_sfc1=0
696  scaled_val_fixed_sfc2=fldlvl2
697  scale_fct_fixed_sfc2=0
698  endif
699 
700  fixed_sfc2_type=pset%param(nprm)%fixed_sfc2_type
701  if(size(pset%param(nprm)%level2)>1.and.size(pset%param(nprm)%level2)<nlvl) then
702  fixed_sfc2_type=''
703  endif
704 
705  if (associated(pset%param(nprm)%level2)) then
706  if(size(pset%param(nprm)%level2)>1.and.size(pset%param(nprm)%level2)>=nlvl) then
707  scaled_val_fixed_sfc2=nint(pset%param(nprm)%level2(nlvl))
708  else if(size(pset%param(nprm)%level2)==1) then
709  scaled_val_fixed_sfc2=nint(pset%param(nprm)%level2(1))
710  else
711  scaled_val_fixed_sfc2=0
712  endif
713  else
714  scaled_val_fixed_sfc2=0
715  end if
716 
717  if(size(pset%param(nprm)%scale_fact_fixed_sfc2)>1 .and. &
718  size(pset%param(nprm)%scale_fact_fixed_sfc2)>=nlvl) then
719  scale_fct_fixed_sfc2=pset%param(nprm)%scale_fact_fixed_sfc2(nlvl)
720  else if(size(pset%param(nprm)%scale_fact_fixed_sfc2)==1) then
721  scale_fct_fixed_sfc2=pset%param(nprm)%scale_fact_fixed_sfc2(1)
722  else
723  scale_fct_fixed_sfc2=0
724  endif
725 
726  ihr_start = ifhr-tinvstat
727  if(modelname=='RAPR'.and.vtimeunits=='FMIN') then
728  ifhrorig = ifhr
729  ifhr = ifhr*60 + ifmin
730  ihr_start = max(0,ifhr-tinvstat)
731  else
732  if(ifmin > 0.)then ! change time range unit to minute
733  pset%time_range_unit="minute"
734  ifhrorig = ifhr
735  ifhr = ifhr*60 + ifmin
736  ihr_start = max(0,ifhr-tinvstat*60)
737  end if
738  end if
739 ! print *,'bf g2sec4_temp0,ipdstmpl=',trim(pset%param(nprm)%pdstmpl),'fixed_sfc_type=', &
740 ! pset%param(nprm)%fixed_sfc1_type,'scale_fct_fixed_sfc1=', &
741 ! scale_fct_fixed_sfc1,'scaled_val_fixed_sfc1=',scaled_val_fixed_sfc1, &
742 ! 'sfc2_type=',trim(pset%param(nprm)%fixed_sfc2_type),scale_fct_fixed_sfc2, &
743 ! scaled_val_fixed_sfc2
744 
745  if(trim(pset%param(nprm)%pdstmpl)=='tmpl4_0') then
746  ipdsnum=0
747  ipdstmpllen=ipdstmp4_0len
748  call g2sec4_temp0(icatg,iparm,pset%gen_proc_type, &
749  pset%gen_proc,hrs_obs_cutoff,min_obs_cutoff, &
750  pset%time_range_unit,ifhr, &
751  pset%param(nprm)%fixed_sfc1_type, &
752  scale_fct_fixed_sfc1, &
753  scaled_val_fixed_sfc1, &
754  fixed_sfc2_type, &
755  scale_fct_fixed_sfc2, &
756  scaled_val_fixed_sfc2, &
757  ipdstmpl(1:ipdstmpllen))
758 ! print *,'aft g2sec4_temp0,ipdstmpl0=',ipdstmpl(1:ipdstmp4_0len)
759  elseif(trim(pset%param(nprm)%pdstmpl)=='tmpl4_1') then
760  ipdsnum=1
761  ipdstmpllen=ipdstmp4_1len
762  call g2sec4_temp1(icatg,iparm,pset%gen_proc_type, &
763  pset%gen_proc,hrs_obs_cutoff,min_obs_cutoff, &
764  pset%time_range_unit,ifhr, &
765  pset%param(nprm)%fixed_sfc1_type, &
766  scale_fct_fixed_sfc1, &
767  scaled_val_fixed_sfc1, &
768  fixed_sfc2_type, &
769  scale_fct_fixed_sfc2, &
770  scaled_val_fixed_sfc2, &
771  pset%type_ens_fcst,perturb_num,num_ens_fcst, &
772  ipdstmpl(1:ipdstmpllen))
773 ! print *,'aft g2sec4_temp1,ipdstmpl1=',ipdstmpl(1:ipdstmp4_1len)
774 !
775  elseif(trim(pset%param(nprm)%pdstmpl)=='tmpl4_8') then
776 !
777  ipdsnum=8
778  ipdstmpllen=ipdstmp4_8len
779  call g2sec4_temp8(icatg,iparm,pset%gen_proc_type, &
780  pset%gen_proc,hrs_obs_cutoff,min_obs_cutoff, &
781  pset%time_range_unit,ihr_start, &
782  pset%param(nprm)%fixed_sfc1_type, &
783  scale_fct_fixed_sfc1, &
784  scaled_val_fixed_sfc1, &
785  pset%param(nprm)%fixed_sfc2_type, &
786  scale_fct_fixed_sfc2, &
787  scaled_val_fixed_sfc2, &
788  idat(3),idat(1),idat(2),idat(4),idat(5), &
789  sec_intvl,ntrange,stat_miss_val, &
790  pset%param(nprm)%stats_proc,type_of_time_inc, &
791  pset%time_range_unit, tinvstat, &
792  stat_unit_time_key_succ,time_inc_betwn_succ_fld, &
793  ipdstmpl(1:ipdstmpllen))
794 ! print *,'aft g2sec4_temp8,ipdstmpl8=',ipdstmpl(1:ipdstmp4_8len)
795 
796  elseif(trim(pset%param(nprm)%pdstmpl)=='tmpl4_11') then
797  ipdsnum=11
798  ipdstmpllen=ipdstmp4_11len
799  call g2sec4_temp11(icatg,iparm,pset%gen_proc_type, &
800  pset%gen_proc,hrs_obs_cutoff,min_obs_cutoff, &
801  pset%time_range_unit,ifhr-tinvstat, &
802  pset%param(nprm)%fixed_sfc1_type, &
803  scale_fct_fixed_sfc1, &
804  scaled_val_fixed_sfc1, &
805  pset%param(nprm)%fixed_sfc2_type, &
806  scale_fct_fixed_sfc2, &
807  scaled_val_fixed_sfc2, &
808  pset%type_ens_fcst,perturb_num,num_ens_fcst, &
809  idat(3),idat(1),idat(2),idat(4),idat(5), &
810  sec_intvl,ntrange,stat_miss_val, &
811  pset%param(nprm)%stats_proc,type_of_time_inc, &
812  pset%time_range_unit, tinvstat, &
813  stat_unit_time_key_succ,time_inc_betwn_succ_fld, &
814  ipdstmpl(1:ipdstmpllen))
815 ! print *,'aft g2sec4_temp11,ipdstmpl11=',ipdstmpl(1:ipdstmp4_11len)
816 
817  elseif(trim(pset%param(nprm)%pdstmpl)=='tmpl4_12') then
818  ipdsnum=12
819  ipdstmpllen=ipdstmp4_12len
820  call g2sec4_temp12(icatg,iparm,pset%gen_proc_type, &
821  pset%gen_proc,hrs_obs_cutoff,min_obs_cutoff, &
822  pset%time_range_unit,ifhr-tinvstat, &
823  pset%param(nprm)%fixed_sfc1_type, &
824  scale_fct_fixed_sfc1, &
825  scaled_val_fixed_sfc1, &
826  fixed_sfc2_type, &
827  scale_fct_fixed_sfc2, &
828  scaled_val_fixed_sfc2, &
829  pset%type_derived_fcst,num_ens_fcst, &
830  idat(3),idat(1),idat(2),idat(4),idat(5), &
831  sec_intvl,ntrange,stat_miss_val, &
832  pset%param(nprm)%stats_proc,type_of_time_inc, &
833  pset%time_range_unit, tinvstat, &
834  stat_unit_time_key_succ,time_inc_betwn_succ_fld, &
835  ipdstmpl(1:ipdstmpllen))
836 ! print *,'aft g2sec4_temp12,ipdstmpl12=',ipdstmpl(1:ipdstmp4_12len)
837 
838  elseif(trim(pset%param(nprm)%pdstmpl)=='tmpl4_44') then
839 !
840  ipdsnum=44
841  ipdstmpllen=ipdstmp4_44len
842  call g2sec4_temp44(icatg,iparm,pset%param(nprm)%aerosol_type, &
843  pset%param(nprm)%typ_intvl_size, &
844  pset%param(nprm)%scale_fact_1st_size, &
845  pset%param(nprm)%scale_val_1st_size, &
846  pset%param(nprm)%scale_fact_2nd_size, &
847  pset%param(nprm)%scale_val_2nd_size, &
848  pset%gen_proc_type, &
849  pset%gen_proc,hrs_obs_cutoff,min_obs_cutoff, &
850  pset%time_range_unit,ifhr, &
851  pset%param(nprm)%fixed_sfc1_type, &
852  scale_fct_fixed_sfc1, &
853  scaled_val_fixed_sfc1, &
854  pset%param(nprm)%fixed_sfc2_type, &
855  scale_fct_fixed_sfc2, &
856  scaled_val_fixed_sfc2, &
857  ipdstmpl(1:ipdstmpllen))
858 ! print *,'aft g2sec4_temp44,ipdstmpl44=',ipdstmpl(1:ipdstmp4_44len),'ipdsnum=',ipdsnum
859 
860  elseif(trim(pset%param(nprm)%pdstmpl)=='tmpl4_48') then
861 !
862  ipdsnum=48
863  ipdstmpllen=ipdstmp4_48len
864  call g2sec4_temp48(icatg,iparm,pset%param(nprm)%aerosol_type, &
865  pset%param(nprm)%typ_intvl_size, &
866  pset%param(nprm)%scale_fact_1st_size, &
867  pset%param(nprm)%scale_val_1st_size, &
868  pset%param(nprm)%scale_fact_2nd_size, &
869  pset%param(nprm)%scale_val_2nd_size, &
870  pset%param(nprm)%typ_intvl_wvlen, &
871  pset%param(nprm)%scale_fact_1st_wvlen, &
872  pset%param(nprm)%scale_val_1st_wvlen, &
873  pset%param(nprm)%scale_fact_2nd_wvlen, &
874  pset%param(nprm)%scale_val_2nd_wvlen, &
875  pset%gen_proc_type, &
876  pset%gen_proc,hrs_obs_cutoff,min_obs_cutoff, &
877  pset%time_range_unit,ifhr, &
878  pset%param(nprm)%fixed_sfc1_type, &
879  scale_fct_fixed_sfc1, &
880  scaled_val_fixed_sfc1, &
881  pset%param(nprm)%fixed_sfc2_type, &
882  scale_fct_fixed_sfc2, &
883  scaled_val_fixed_sfc2, &
884  ipdstmpl(1:ipdstmpllen))
885 ! print *,'aft g2sec4_temp48,name=',trim(pset%param(nprm)%shortname),&
886 ! 'ipdstmpl48=',ipdstmpl(1:ipdstmp4_48len)
887 
888  endif
889 
890  if(modelname=='RAPR'.and.vtimeunits=='FMIN') then
891  ifhr = ifhrorig
892  end if
893  if(ifmin>0.)then
894  ifhr = ifhrorig
895  end if
896 
897 
898 !
899 !----------
900 ! idrstmpl array is the output from g2sec5
901 !
902  call get_g2_sec5packingmethod(pset%packing_method,idrsnum,ierr)
903  if(maxval(datafld1)==minval(datafld1))then
904  idrsnum=0
905 ! print*,' changing to simple packing for constant fields'
906  end if
907  if(modelname=='RAPR') then
908  if((abs(maxval(datafld1)-minval(datafld1)) < 1.1) .and. (datafld1(1) > 500.0))then
909  idrsnum=0
910 ! print*,' changing to simple packing for constant fields: max-min < 0.1'
911  end if
912 
913  if(trim(pset%param(nprm)%shortname)=='UGRD_ON_SPEC_HGT_LVL_ABOVE_GRND_10m'.or.&
914  trim(pset%param(nprm)%shortname)=='VGRD_ON_SPEC_HGT_LVL_ABOVE_GRND_10m'.or.&
915  trim(pset%param(nprm)%shortname)=='VWSH_ON_SPEC_HGT_LVL_ABOVE_GRND'.or.&
916  trim(pset%param(nprm)%shortname)=='VUCSH_ON_SPEC_HGT_LVL_ABOVE_GRND_0-1km'.or.&
917  trim(pset%param(nprm)%shortname)=='VVCSH_ON_SPEC_HGT_LVL_ABOVE_GRND_0-1km'.or.&
918  trim(pset%param(nprm)%shortname)=='VUCSH_ON_SPEC_HGT_LVL_ABOVE_GRND_0-6km'.or.&
919  trim(pset%param(nprm)%shortname)=='VVCSH_ON_SPEC_HGT_LVL_ABOVE_GRND_0-6km')then
920  idrsnum=0
921 ! print*,' changing to simple packing for field: ',trim(pset%param(nprm)%shortname)
922  endif
923  endif
924 
925 ! print *,'aft g2sec5,packingmethod=',pset%packing_method,'idrsnum=',idrsnum, &
926 ! 'data=',maxval(datafld1),minval(datafld1)
927 !
928 !*** set number of bits, and binary scale
929 !
930  ibmap=255
931  bmap=.true.
932  if(any(datafld1>=spval))then
933  ibmap=0
934  where(abs(datafld1)>=spval)bmap=.false.
935  endif
936 !
937  if(size(pset%param(nprm)%level)==size(pset%param(nprm)%scale)) then
938  if(size(pset%param(nprm)%level)==1) nlvl=1
939  if(nlvl/=0) then
940  fldscl=nint(pset%param(nprm)%scale(nlvl))
941  else
942  fldscl=0
943  endif
944  else if (size(pset%param(nprm)%scale)==1) then
945  fldscl=nint(pset%param(nprm)%scale(1))
946  endif
947 !
948  mxbit=16
949  if(trim(pset%param(nprm)%pname)=='APCP' .or. &
950  trim(pset%param(nprm)%pname)=='ACPCP' .or. &
951  trim(pset%param(nprm)%pname)=='NCPCP')mxbit=22
952  if(mxbit>16)print*, 'increased MXBIT for ', pset%param(nprm)%pname,mxbit
953  call g2getbits(mxbit,ibmap,fldscl,size(datafld1),bmap,datafld1,ibin_scl,idec_scl,inumbits)
954 ! print *,'idec_scl=',idec_scl,'ibin_scl=',ibin_scl,'number_bits=',inumbits
955  if( idrsnum==40 ) then
956  idrstmplen=idrstmp5_40len
957  call g2sec5_temp40(idec_scl,ibin_scl,inumbits,pset%comprs_type,idrstmpl(1:idrstmplen))
958  elseif( idrsnum==2 ) then
959  idrstmplen=idrstmp5_2len
960  call g2sec5_temp2(idec_scl,ibin_scl,idrstmpl(1:idrstmplen))
961  elseif( idrsnum==3 ) then
962  idrstmplen=idrstmp5_3len
963  call g2sec5_temp3(idec_scl,ibin_scl,pset%order_of_sptdiff,idrstmpl(1:idrstmplen))
964  elseif( idrsnum==0 ) then
965  idrstmplen=idrstmp5_0len
966  call g2sec5_temp0(idec_scl,ibin_scl,inumbits,idrstmpl(1:idrstmplen))
967  endif
968 !
969 !----------------------------------------------------------------------------------------
970 ! Define all required inputs like ibmap, numcoord, coordlist etc externally in the module
971 ! prior to calling the addfield routine. Again hide the addfield routine from the user
972 !
973 ! print *,'before addfield, data=',maxval(datafld1),minval(datafld1),'ibmap=',ibmap, &
974 ! 'max_bytes=',max_bytes,'ipdsnum=',ipdsnum,'ipdstmpllen=',ipdstmpllen,'ipdstmpl=',ipdstmpl(1:ipdstmpllen), &
975 ! 'coordlist=',coordlist,'numcoord=',numcoord,'idrsnum=',idrsnum,'idrstmpl=',idrstmpl, &
976 ! 'idrstmplen=',idrstmplen,'im_jm=',im_jm
977 
978  call addfield(cgrib,max_bytes,ipdsnum,ipdstmpl(1:ipdstmpllen), &
979  ipdstmpllen,coordlist,numcoord,idrsnum,idrstmpl, &
980  idrstmplen ,datafld1,im_jm,ibmap,bmap,ierr)
981 !
982 !---------------------------------------------------------------------------------------
983 ! Finalize GRIB message after all grids and fields have been added by adding the end
984 ! section "7777"
985 ! Again hide the gribend routine from the user
986 !
987  call gribend(cgrib,max_bytes,lengrib,ierr)
988 !
989 !-------
990  end subroutine gengrb2msg
991 !
992 !--------------------------------------------------------------------------------------
993 !
994 ! E. JAMES: 10 JUN 2021 - Adding section to read in GRIB2 files for comparison
995 ! within UPP. Two new subroutines added below.
996 !
997  subroutine read_grib2_head(filenameG2,nx,ny,nz,rlonmin,rlatmax,rdx,rdy)
998 !
999 !--- read grib2 file head information
1000 !
1001  use grib_mod
1002  implicit none
1003  character*256,intent(in) :: filenameg2
1004  integer, intent(out) :: nx,ny,nz
1005  real, intent(out) :: rlonmin,rlatmax
1006  real*8, intent(out) :: rdx,rdy
1007 !
1008 !
1009  type(gribfield) :: gfld
1010  logical :: expand=.true.
1011  integer :: ifile
1012  character(len=1),allocatable,dimension(:) :: cgrib
1013  integer,parameter :: msk1=32000
1014  integer :: lskip, lgrib,iseek
1015  integer :: currlen
1016  integer :: icount , lengrib
1017  integer :: listsec0(3)
1018  integer :: listsec1(13)
1019  integer year, month, day, hour, minute, second, fcst
1020  integer :: numfields,numlocal,maxlocal,ierr
1021  integer :: grib_edition
1022  integer :: itot
1023 ! real :: dx,dy,lat1,lon1
1024  real :: scale_factor,scale_factor2
1025 !
1026 !
1027  integer :: nn,n,j,iret
1028  real :: fldmax,fldmin,sum
1029 !
1030 !
1031  scale_factor=1.0e6
1032  scale_factor2=1.0e3
1033  ifile=10
1034  loopfile: do nn=1,1
1035 ! write(6,*) 'read in grib2 file head', trim(filenameG2)
1036  lskip=0
1037  lgrib=0
1038  iseek=0
1039  icount=0
1040  itot=0
1041  currlen=0
1042 ! Open GRIB2 file
1043  call baopenr(ifile,trim(filenameg2),iret)
1044  if (iret.eq.0) then
1045  version: do
1046  ! Search opend file for the next GRIB2 messege (record).
1047  call skgb(ifile,iseek,msk1,lskip,lgrib)
1048  ! Check for EOF, or problem
1049  if (lgrib.eq.0) then
1050  exit
1051  endif
1052  ! Check size, if needed allocate more memory.
1053  if (lgrib.gt.currlen) then
1054  if (allocated(cgrib)) deallocate(cgrib)
1055  allocate(cgrib(lgrib))
1056  currlen=lgrib
1057  endif
1058  ! Read a given number of bytes from unblocked file.
1059  call baread(ifile,lskip,lgrib,lengrib,cgrib)
1060  if(lgrib.ne.lengrib) then
1061  write(*,*) 'ERROR, read_grib2 lgrib ne lengrib', &
1062  lgrib,lengrib
1063  stop 1234
1064  endif
1065  iseek=lskip+lgrib
1066  icount=icount+1
1067  ! Unpack GRIB2 field
1068  call gb_info(cgrib,lengrib,listsec0,listsec1, &
1069  numfields,numlocal,maxlocal,ierr)
1070  if(ierr.ne.0) then
1071  write(6,*) 'Error querying GRIB2 message',ierr
1072  stop
1073  endif
1074  itot=itot+numfields
1075  grib_edition=listsec0(2)
1076  if (grib_edition.ne.2) then
1077  exit version
1078  endif
1079 ! write(*,*) 'listsec0=',listsec0
1080 ! write(*,*) 'listsec1=',listsec1
1081 ! write(*,*) 'numfields=',numfields
1082 ! get information form grib2 file
1083  n=1
1084  call gf_getfld(cgrib,lengrib,n,.false.,expand,gfld,ierr)
1085  year =gfld%idsect(6) !(FOUR-DIGIT) YEAR OF THE DATA
1086  month =gfld%idsect(7) ! MONTH OF THE DATA
1087  day =gfld%idsect(8) ! DAY OF THE DATA
1088  hour =gfld%idsect(9) ! HOUR OF THE DATA
1089  minute=gfld%idsect(10) ! MINUTE OF THE DATA
1090  second=gfld%idsect(11) ! SECOND OF THE DATA
1091  write(*,*) 'year,month,day,hour,minute,second='
1092  write(*,*) year,month,day,hour,minute,second
1093  write(*,*) 'source center =',gfld%idsect(1)
1094  write(*,*) 'Indicator of model =',gfld%ipdtmpl(5)
1095  write(*,*) 'observation level (m)=',gfld%ipdtmpl(12)
1096  write(*,*) 'map projection=',gfld%igdtnum
1097  if (gfld%igdtnum.eq.0) then ! Lat/Lon grid aka Cylindrical
1098  ! Equidistant
1099  nx = gfld%igdtmpl(8)
1100  ny = gfld%igdtmpl(9)
1101  nz = 1
1102  rdx = gfld%igdtmpl(17)/scale_factor
1103  rdy = gfld%igdtmpl(18)/scale_factor
1104  rlatmax = gfld%igdtmpl(12)/scale_factor
1105  rlonmin = gfld%igdtmpl(13)/scale_factor
1106 ! write(*,*) 'nx,ny=',nx,ny
1107 ! write(*,*) 'dx,dy=',rdx,rdy
1108 ! write(*,*) 'lat1,lon1=',rlatmax,rlonmin
1109  else if (gfld%igdtnum.eq.1) then ! Rotated Lat Lon Grid (RRFS_NA)
1110  nx = gfld%igdtmpl(8)
1111  ny = gfld%igdtmpl(9)
1112  nz = 1
1113  rdx = gfld%igdtmpl(17)/scale_factor
1114  rdy = gfld%igdtmpl(18)/scale_factor
1115  rlatmax = gfld%igdtmpl(12)/scale_factor
1116  rlonmin = gfld%igdtmpl(13)/scale_factor
1117 ! write(*,*) 'nx,ny=',nx,ny
1118 ! write(*,*) 'dx,dy=',rdx,rdy
1119 ! write(*,*) 'lat1,lon1=',rlatmax,rlonmin
1120  else if (gfld%igdtnum.eq.30) then ! Lambert Conformal Grid (HRRR)
1121  nx = gfld%igdtmpl(8)
1122  ny = gfld%igdtmpl(9)
1123  nz = 1
1124  rdx = gfld%igdtmpl(15)/scale_factor2
1125  rdy = gfld%igdtmpl(16)/scale_factor2
1126  rlatmax = gfld%igdtmpl(10)/scale_factor
1127  rlonmin = gfld%igdtmpl(11)/scale_factor
1128 ! write(*,*) 'nx,ny=',nx,ny
1129 ! write(*,*) 'dx,dy=',rdx,rdy
1130 ! write(*,*) 'lat1,lon1=',rlatmax,rlonmin
1131  else
1132  write(*,*) 'unknown projection'
1133  stop 1235
1134  endif
1135  call gf_free(gfld)
1136  enddo version ! skgb
1137  endif
1138  CALL baclose(ifile,ierr)
1139  nullify(gfld%local)
1140  if (allocated(cgrib)) deallocate(cgrib)
1141  enddo loopfile
1142  return
1143  end subroutine read_grib2_head
1144 !
1145 !---
1146 !
1147  subroutine read_grib2_sngle(filenameG2,ntot,height,var)
1148 !
1149 !--- read grib2 files
1150 !
1151  use grib_mod
1152  implicit none
1153  character*256,intent(in) :: filenameg2
1154  integer, intent(in) :: ntot
1155  real, intent(out) :: var(ntot)
1156  integer, intent(out) :: height
1157 !
1158 !
1159  type(gribfield) :: gfld
1160  logical :: expand=.true.
1161  integer :: ifile
1162  character(len=1),allocatable,dimension(:) :: cgrib
1163  integer,parameter :: msk1=32000
1164  integer :: lskip, lgrib,iseek
1165  integer :: currlen
1166  integer :: icount , lengrib
1167  integer :: listsec0(3)
1168  integer :: listsec1(13)
1169  integer year, month, day, hour, minute, second, fcst
1170  integer :: numfields,numlocal,maxlocal,ierr
1171  integer :: grib_edition
1172  integer :: itot
1173  integer :: nx,ny
1174  real :: dx,dy,lat1,lon1,rtnum
1175  real :: ref_value,bin_scale_fac,dec_scale_fac,bit_number,field_type
1176  real :: bit_map
1177  real :: scale_factor,scale_factor2
1178 !
1179 !
1180  integer :: nn,n,j,iret
1181  real :: fldmax,fldmin,sum
1182 !
1183 !
1184  scale_factor=1.0e6
1185  scale_factor2=1.0e3
1186  ifile=12
1187  loopfile: do nn=1,1
1188 ! write(6,*) 'read mosaic in grib2 file ', trim(filenameG2)
1189  lskip=0
1190  lgrib=0
1191  iseek=0
1192  icount=0
1193  itot=0
1194  currlen=0
1195 ! Open GRIB2 file
1196  call baopenr(ifile,trim(filenameg2),iret)
1197  if (iret.eq.0) then
1198  version: do
1199  ! Search opend file for the next GRIB2 messege (record).
1200  call skgb(ifile,iseek,msk1,lskip,lgrib)
1201  ! Check for EOF, or problem
1202  if (lgrib.eq.0) then
1203  exit
1204  endif
1205  ! Check size, if needed allocate more memory.
1206  if (lgrib.gt.currlen) then
1207  if (allocated(cgrib)) deallocate(cgrib)
1208  allocate(cgrib(lgrib))
1209  currlen=lgrib
1210  endif
1211  ! Read a given number of bytes from unblocked file.
1212  call baread(ifile,lskip,lgrib,lengrib,cgrib)
1213  if(lgrib.ne.lengrib) then
1214  write(*,*) 'ERROR, read_grib2 lgrib ne lengrib', &
1215  lgrib,lengrib
1216  stop 1234
1217  endif
1218 ! write(*,*) 'lengrib=',lengrib
1219  iseek=lskip+lgrib
1220  icount=icount+1
1221  ! Unpack GRIB2 field
1222  call gb_info(cgrib,lengrib,listsec0,listsec1, &
1223  numfields,numlocal,maxlocal,ierr)
1224  if(ierr.ne.0) then
1225  write(6,*) 'Error querying GRIB2 message',ierr
1226  stop
1227  endif
1228  itot=itot+numfields
1229  grib_edition=listsec0(2)
1230  if (grib_edition.ne.2) then
1231  exit version
1232  endif
1233 ! write(*,*) 'listsec0=',listsec0
1234 ! write(*,*) 'listsec1=',listsec1
1235 ! write(*,*) 'numfields=',numfields!
1236 ! get information form grib2 file
1237  n=1
1238  call gf_getfld(cgrib,lengrib,n,.false.,expand,gfld,ierr)
1239  year =gfld%idsect(6) !(FOUR-DIGIT) YEAR OF THE DATA
1240  month =gfld%idsect(7) ! MONTH OF THE DATA
1241  day =gfld%idsect(8) ! DAY OF THE DATA
1242  hour =gfld%idsect(9) ! HOUR OF THE DATA
1243  minute=gfld%idsect(10) ! MINUTE OF THE DATA
1244  second=gfld%idsect(11) ! SECOND OF THE DATA
1245 ! write(*,*) 'year,month,day,hour,minute,second='
1246 ! write(*,*) year,month,day,hour,minute,second
1247 ! write(*,*) 'source center =',gfld%idsect(1)
1248 ! write(*,*) 'Indicator of model =',gfld%ipdtmpl(5)
1249 ! write(*,*) 'observation level (m)=',gfld%ipdtmpl(12)
1250 ! write(*,*) 'map projection=',gfld%igdtnum
1251  height=gfld%ipdtmpl(12)
1252  if (gfld%igdtnum.eq.0) then ! Lat/Lon grid aka Cylindrical
1253  ! Equidistant
1254  nx = gfld%igdtmpl(8)
1255  ny = gfld%igdtmpl(9)
1256  dx = gfld%igdtmpl(17)/scale_factor
1257  dy = gfld%igdtmpl(18)/scale_factor
1258  lat1 = gfld%igdtmpl(12)/scale_factor
1259  lon1 = gfld%igdtmpl(13)/scale_factor
1260 ! write(*,*) 'nx,ny=',nx,ny
1261 ! write(*,*) 'dx,dy=',dx,dy
1262 ! write(*,*) 'lat1,lon1=',lat1,lon1
1263  else if (gfld%igdtnum.eq.1) then ! Rotated Lat Lon Grid (RRFS_NA)
1264  nx = gfld%igdtmpl(8)
1265  ny = gfld%igdtmpl(9)
1266  dx = gfld%igdtmpl(17)/scale_factor
1267  dy = gfld%igdtmpl(18)/scale_factor
1268  lat1 = gfld%igdtmpl(12)/scale_factor
1269  lon1 = gfld%igdtmpl(13)/scale_factor
1270 ! write(*,*) 'nx,ny=',nx,ny
1271 ! write(*,*) 'dx,dy=',rdx,rdy
1272 ! write(*,*) 'lat1,lon1=',rlatmax,rlonmin
1273  else if (gfld%igdtnum.eq.30) then ! Lambert Conformal Grid (HRRR)
1274  nx = gfld%igdtmpl(8)
1275  ny = gfld%igdtmpl(9)
1276  dx = gfld%igdtmpl(15)/scale_factor2
1277  dy = gfld%igdtmpl(16)/scale_factor2
1278  lat1 = gfld%igdtmpl(10)/scale_factor
1279  lon1 = gfld%igdtmpl(11)/scale_factor
1280 ! write(*,*) 'In read_grib2_sngle:'
1281 ! write(*,*) 'nx,ny=',nx,ny
1282 ! write(*,*) 'dx,dy=',dx,dy
1283 ! write(*,*) 'lat1,lon1=',lat1,lon1
1284  rtnum = gfld%idrtnum
1285 ! write(*,*) 'rtnum=',rtnum
1286  ref_value = gfld%idrtmpl(1)
1287  bin_scale_fac = gfld%idrtmpl(2)
1288  dec_scale_fac = gfld%idrtmpl(3)
1289  bit_number = gfld%idrtmpl(4)
1290  field_type = gfld%idrtmpl(5)
1291  bit_map = gfld%ibmap
1292 ! write(*,*) 'ref_value=',ref_value
1293 ! write(*,*) 'bin_scale_fac=',bin_scale_fac
1294 ! write(*,*) 'dec_scale_fac=',dec_scale_fac
1295 ! write(*,*) 'bit_number=',bit_number
1296 ! write(*,*) 'field_type=',field_type
1297 ! write(*,*) 'bit map indicator=',bit_map
1298  else
1299  write(*,*) 'unknown projection'
1300  stop 1235
1301  endif
1302  call gf_free(gfld)
1303  ! Continue to unpack GRIB2 field.
1304  num_fields: do n = 1, numfields
1305  ! e.g. U and V would =2, otherwise its usually =1
1306  call gf_getfld(cgrib,lengrib,n,.true.,expand,gfld,ierr)
1307  if (ierr.ne.0) then
1308  write(*,*) ' ERROR extracting field gf_getfld = ',ierr
1309  cycle
1310  endif
1311 ! write(*,*) 'gfld%ndpts=',n,gfld%ndpts
1312 ! write(*,*) 'gfld%ngrdpts=',n,gfld%ngrdpts
1313 ! write(*,*) 'gfld%unpacked=',n,gfld%unpacked
1314  fldmax=gfld%fld(1)
1315  fldmin=gfld%fld(1)
1316  sum=gfld%fld(1)
1317  if(ntot .ne. gfld%ngrdpts) then
1318  write(*,*) 'Error, wrong dimension ',ntot, gfld%ngrdpts
1319  stop 1234
1320  endif
1321  do j=1,gfld%ngrdpts
1322  var(j)=gfld%fld(j)
1323  enddo
1324 ! write(*,*) 'j,first,last:',j,var(954370),var(953920)
1325 ! write(*,*) 'height,max,min',height,maxval(var),minval(var)
1326  call gf_free(gfld)
1327  enddo num_fields
1328  enddo version ! skgb
1329  endif
1330  CALL baclose(ifile,ierr)
1331  if (allocated(cgrib)) deallocate(cgrib)
1332  nullify(gfld%local)
1333  enddo loopfile
1334  return
1335  end subroutine read_grib2_sngle
1336 !
1337 !----------------------------------------------------------------------------------------
1338 !
1339  subroutine g2sec3tmpl40(nx,nY,lat1,lon1,lat2,lon2,lad,ds1,len3,igds,ifield3)
1340  implicit none
1341 !
1342  integer(4),intent(inout) :: nx,ny,lat1,lon1,lat2,lon2,lad,ds1,len3
1343  integer(4),intent(inout) :: ifield3(len3),igds(5)
1344 !
1345  nx=1152
1346  ny=576
1347  lat1=89761000 !lat of 1st grd pt in micro-deg
1348  lon1=0 !east-long of 1st grd pt in micro-deg
1349  lat2=-89761000 !lat of last grd pt in micro-deg
1350  lon2=359687000 !east-long of last grd pt in micro-deg
1351  lad=313000 !lat at which projection intersects earth
1352  ds1=288 !grid spacing in x and y
1353 !
1354  igds=0
1355  igds(2)=nx*ny
1356  igds(5)=40
1357 !
1358  ifield3(1) = 6 !Earth assumed spherical with radius of 6,371,229.0m
1359  ifield3(2) = 0
1360  ifield3(3) = 0
1361  ifield3(4) = 0
1362  ifield3(5) = 0
1363  ifield3(6) = 0
1364  ifield3(7) = 0
1365  ifield3(8) = nx
1366  ifield3(9) = ny
1367  ifield3(10) = 0
1368  ifield3(11) = 0
1369  ifield3(12) = lat1
1370  ifield3(13) = lon1
1371  ifield3(14) = 48
1372  ifield3(15) = lat2
1373  ifield3(16) = lon2
1374  ifield3(17) = lad
1375  ifield3(18) = ds1
1376  ifield3(19) = 0
1377 
1378  return
1379  end subroutine g2sec3tmpl40
1380 !
1381 !-------------------------------------------------------------------------------------
1382 !
1383  subroutine g2getbits(MXBIT,ibm,scl,len,bmap,g,ibs,ids,nbits)
1384 !$$$
1385 ! This subroutine is changed from w3 lib getbit to compute the total number of bits,
1386 ! The argument list is modified to have ibm,scl,len,bmap,g,ibs,ids,nbits
1387 !
1388 ! Progrma log:
1389 ! Jun Wang Apr, 2010
1390 !
1391 ! INPUT
1392 ! ibm: integer, bitmap flag (grib2 table 6.0)
1393 ! scl: real, significant digits,OR binary precision if < 0
1394 ! len: integer, field and bitmap length
1395 ! bmap: logical(len), bitmap (.true.: keep, bitmap (.true.: keep, .false. skip)
1396 ! fld: real(len), datafield
1397 ! OUTPUT
1398 ! ibs: integer, binary scale factor
1399 ! ids: integer, decimal scale factor
1400 ! nbits: integer, number of bits to pack
1401 !
1402  IMPLICIT NONE
1403 !
1404  INTEGER,INTENT(IN) :: mxbit,ibm,len
1405  LOGICAL*1,INTENT(IN) :: bmap(len)
1406  REAL,INTENT(IN) :: scl,g(len)
1407  INTEGER,INTENT(OUT) :: ibs,ids,nbits
1408 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1409 ! INTEGER,PARAMETER :: MXBIT=16
1410 !
1411 ! NATURAL LOGARITHM OF 2 AND 0.5 PLUS NOMINAL SAFE EPSILON
1412  real,PARAMETER :: alog2=0.69314718056,hpeps=0.500001
1413 !
1414 !local vars
1415  INTEGER :: i,i1,icnt,ipo,le,irange
1416  REAL :: ground,gmin,gmax,s,rmin,rmax,range,rr,rng2,po,rln2
1417 !
1418  DATA rln2/0.69314718/
1419 
1420 
1421 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1422 ! ROUND FIELD AND DETERMINE EXTREMES WHERE BITMAP IS ON
1423  IF(ibm == 255) THEN
1424  gmax = g(1)
1425  gmin = g(1)
1426  DO i=2,len
1427  gmax = max(gmax,g(i))
1428  gmin = min(gmin,g(i))
1429  ENDDO
1430  ELSE
1431  do i1=1,len
1432  if (bmap(i1)) exit
1433  enddo
1434 ! I1 = 1
1435 ! DO WHILE(I1 <= LEN .AND. .not. BMAP(I1))
1436 ! I1=I1+1
1437 ! ENDDO
1438  IF(i1 <= len) THEN
1439  gmax = g(i1)
1440  gmin = g(i1)
1441  DO i=i1+1,len
1442  IF(bmap(i)) THEN
1443  gmax = max(gmax,g(i))
1444  gmin = min(gmin,g(i))
1445  ENDIF
1446  ENDDO
1447  ELSE
1448  gmax = 0.
1449  gmin = 0.
1450  ENDIF
1451  ENDIF
1452 
1453 ! write(0,*)' GMIN=',GMIN,' GMAX=',GMAX
1454 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1455 ! COMPUTE NUMBER OF BITS
1456  icnt = 0
1457  ibs = 0
1458  ids = 0
1459  range = gmax - gmin
1460 ! IF ( range <= 0.00 ) THEN
1461  IF ( range <= 1.e-30 ) THEN
1462  nbits = 8
1463  return
1464  END IF
1465 !*
1466  IF ( scl == 0.0 ) THEN
1467  nbits = 8
1468  RETURN
1469  ELSE IF ( scl > 0.0 ) THEN
1470  ipo = int(alog10( range ))
1471 !jw: if range is smaller than computer precision, set nbits=8
1472  if(ipo<0.and.ipo+scl<-20) then
1473  print *,'for small range,ipo=',ipo,'ipo+scl=',ipo+scl,'scl=',scl
1474  nbits=8
1475  return
1476  endif
1477 
1478  IF ( range < 1.00 ) ipo = ipo - 1
1479  po = float(ipo) - scl + 1.
1480  ids = - int( po )
1481  rr = range * 10. ** ( -po )
1482  nbits = int( alog( rr ) / rln2 ) + 1
1483  ELSE
1484  ibs = -nint( -scl )
1485  rng2 = range * 2. ** (-ibs)
1486  nbits = int( alog( rng2 ) / rln2 ) + 1
1487  END IF
1488 ! write(0,*)'in g2getnits,ibs=',ibs,'ids=',ids,'nbits=',nbits,'range=',range
1489 !*
1490  IF(nbits <= 0) THEN
1491  nbits = 0
1492  IF(abs(gmin) >= 1.) THEN
1493  ids = -int(alog10(abs(gmin)))
1494  ELSE IF (abs(gmin) < 1.0.AND.abs(gmin) > 0.0) THEN
1495  ids = -int(alog10(abs(gmin)))+1
1496  ELSE
1497  ids = 0
1498  ENDIF
1499  ENDIF
1500  nbits = min(nbits,mxbit)
1501 ! write(0,*)'in g2getnits ibs=',ibs,'ids=',ids,'nbits=',nbits
1502 !
1503  IF ( scl > 0.0 ) THEN
1504  s=10.0 ** ids
1505  IF(ibm == 255) THEN
1506  ground = g(1)*s
1507  gmax = ground
1508  gmin = ground
1509  DO i=2,len
1510  gmax = max(gmax,g(i)*s)
1511  gmin = min(gmin,g(i)*s)
1512  ENDDO
1513  ELSE
1514  do i1=1,len
1515  if (bmap(i1)) exit
1516  enddo
1517  ! I1=1
1518  ! DO WHILE(I1<=LEN.AND..not.BMAP(I1))
1519  ! I1=I1+1
1520  ! ENDDO
1521  IF(i1 <= len) THEN
1522  ground = g(i1)*s
1523  gmax = ground
1524  gmin = ground
1525  DO i=i1+1,len
1526  IF(bmap(i)) THEN
1527  gmax = max(gmax,g(i)*s)
1528  gmin = min(gmin,g(i)*s)
1529  ENDIF
1530  ENDDO
1531  ELSE
1532  gmax = 0.
1533  gmin = 0.
1534  ENDIF
1535  ENDIF
1536 
1537  range = gmax-gmin
1538  if(gmax == gmin) then
1539  ibs = 0
1540  else
1541  ibs = nint(alog(range/(2.**nbits-0.5))/alog2+hpeps)
1542  endif
1543 !
1544  endif
1545 ! write(0,*)'in g2getnits,2ibs=',ibs,'ids=',ids,'nbits=',nbits,'range=',&
1546 ! range, 'scl=',scl,'data=',maxval(g),minval(g)
1547 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1548  RETURN
1549  END subroutine g2getbits
1550 !
1551 !-------------------------------------------------------------------------------------
1552 !
1553  subroutine getgds(ldfgrd,len3,ifield3len,igds,ifield3)
1554 !
1555 !***** set up gds kpds to call Boi's code
1556 !
1557  use ctlblk_mod, only : im,jm,gdsdegr,modelname
1558  use gridspec_mod, only: dxval,dyval,cenlat,cenlon,latstart,lonstart,latlast, &
1559  & lonlast,maptype,standlon,latstartv,cenlatv,lonstartv, &
1560  cenlonv,truelat1,truelat2,latstart_r,lonstart_r, &
1561  latlast_r,lonlast_r
1562 !
1563  implicit none
1564 !
1565  logical, intent(in) :: ldfgrd
1566  integer(4),intent(in) :: len3
1567  integer(4),intent(out) :: ifield3len
1568  integer(4),intent(inout) :: ifield3(len3),igds(5)
1569 
1570 ! print *,'in getgds, im=',im,'jm=',jm,'latstart=',latstart,'lonsstart=',lonstart,'maptyp=',maptype
1571 !
1572 !** set up igds
1573  igds(1) = 0 !Source of grid definition (see Code Table 3.0)
1574  igds(2) = im*jm !Number of grid points in the defined grid.
1575  igds(3) = 0 !Number of octets needed for each additional grid points definition
1576  igds(4) = 0 !Interpretation of list for optional points definition (Code Table 3.11)
1577 !
1578 !** define grid template 3
1579  IF(maptype == 1) THEN !LAmbert Conformal
1580  igds(5) = 30 !Lambert conformal
1581  ifield3len = 22
1582  ifield3 = 0
1583 !
1584  ifield3(1) = 6 !Earth assumed spherical with radius of 6,371,229.0m
1585  ifield3(8) = im !number of points along the x-axis
1586  ifield3(9) = jm !number of points along the y-axis
1587  ifield3(10) = latstart !latitude of first grid point
1588  ifield3(11) = lonstart !longitude of first grid point
1589  ifield3(12) = 8 !Resolution and component flags
1590 ! Jili Dong change grid to earth relative
1591  if (modelname == 'FV3R') then
1592  ifield3(12) = 0 !Resolution and component flags
1593  endif
1594 
1595  ifield3(13) = truelat1
1596  ifield3(14) = standlon !longitude of meridian parallel to y-axis along which latitude increases
1597  ifield3(15) = dxval
1598  ifield3(16) = dyval
1599  IF(truelat1>0)then
1600  ifield3(17) = 0
1601  else
1602  ifield3(17) =128 !for southern hemisphere
1603  end if
1604  ifield3(18) = 64
1605  ifield3(19) = truelat1 !first latitude from the pole at which the secant cone cuts the sphere
1606  ifield3(20) = truelat2 !second latitude from the pole at which the secant cone cuts the sphere
1607 
1608 !** Polar stereographic
1609  ELSE IF(maptype == 2)THEN !Polar stereographic
1610  igds(5) = 20
1611  ifield3len = 22
1612  ifield3 = 0
1613 !
1614  ifield3(1) = 6 !Earth assumed spherical with radius of 6,371,229.0m
1615  ifield3(8) = im !number of points along the x-axis
1616  ifield3(9) = jm !number of points along the y-axis
1617  ifield3(10) = latstart !latitude of first grid point
1618  ifield3(11) = lonstart !longitude of first grid point
1619  ifield3(12) = 8 !Resolution and component flags
1620  ifield3(13) = truelat1
1621  ifield3(14) = standlon !longitude of meridian parallel to y-axis along which latitude increases
1622  ifield3(15) = dxval
1623  ifield3(16) = dyval
1624  ifield3(17) = 0
1625  ifield3(18) = 64
1626 !
1627 !** Mercator
1628  ELSE IF(maptype==3)THEN !Mercator
1629  igds(5)=10
1630  ifield3len=22
1631  ifield3=0
1632 !
1633  ifield3(1) = 6 !Earth assumed spherical with radius of 6,371,229.0m
1634  ifield3(8) = im !number of points along the x-axis
1635  ifield3(9) = jm !number of points along the y-axis
1636  ifield3(10) = latstart !latitude of first grid point
1637  ifield3(11) = lonstart !longitude of first grid point
1638  ifield3(12) = 8 !Resolution and component flags
1639  ifield3(13) = truelat1 !latitude(s) at which the Mercator projection intersects the Earth
1640  ifield3(14) = latlast !latitude of last grid point
1641  ifield3(15) = lonlast !longitude of last grid point
1642  ifield3(16) = 64 !Scanning mode
1643  ifield3(17) = 0 !Orientation of the grid, angle between i direction on the map and Equator
1644  ifield3(18) = dxval
1645  ifield3(19) = dyval
1646 !
1647 !** ARAKAWA STAGGERED E-GRID
1648  ELSE IF(maptype == 203)THEN !ARAKAWA STAGGERED E-GRID`
1649  igds(5) = 32768
1650  ifield3len = 22
1651  ifield3 = 0
1652 !
1653  ifield3(1) = 6 !Earth assumed spherical with radius of 6,371,229.0m
1654  ifield3(8) = im !number of points along the x-axis
1655  ifield3(9) = jm !number of points along the y-axis
1656  ifield3(10) = 0 !Basic angle of the initial production domain
1657  if(.not.ldfgrd) then
1658  ifield3(11) = 0 !Subdivisions of basic angle used to define extreme lons & lats:missing
1659  else
1660  ifield3(11) = 45000000 !Subdivisions of basic angle used to define extreme lons & lats
1661  endif
1662  ifield3(12) = latstart !latitude of first grid point
1663  ifield3(13) = lonstart !longitude of first grid point
1664  ifield3(14) = 0 !Resolution and component flags
1665  ifield3(15) = cenlat !center latitude of grid point
1666  ifield3(16) = cenlon !Center longitude of grid point
1667  ifield3(17) = dxval
1668  ifield3(18) = dyval
1669  ifield3(19) = 64 !Scanning mode
1670 !
1671 !** ARAKAWA STAGGERED non-E-GRID
1672  ELSE IF(maptype == 205 .OR. maptype == 6)THEN !ARAKAWA STAGGERED E-GRID`
1673  igds(5) = 32769
1674  ifield3len = 21
1675  ifield3 = 0
1676 !
1677  ifield3(1) = 6 !Earth assumed spherical with radius of 6,371,229.0m
1678  ifield3(8) = im !number of points along the x-axis
1679  ifield3(9) = jm !number of points along the y-axis
1680  ifield3(10) = 0 !Basic angle of the initial production domain
1681  if(.not.ldfgrd) then
1682  ifield3(11) = 0 !Subdivisions of basic angle used to define extreme lons & lats:missing
1683  else
1684  ifield3(11) = 45000000 !Subdivisions of basic angle used to define extreme lons & lats
1685  endif
1686  ifield3(12) = latstart !latitude of first grid point
1687  ifield3(13) = lonstart !longitude of first grid point
1688  ifield3(14) = 56 !Resolution and component flags
1689  ifield3(15) = cenlat !center latitude of grid point
1690  ifield3(16) = cenlon !Center longitude of grid point
1691  ifield3(17) = dxval
1692  ifield3(18) = dyval
1693  ifield3(19) = 64 !Scanning mode
1694  ifield3(20) = latlast !Scanning mode
1695  ifield3(21) = lonlast !Scanning mode
1696 !
1697 !** ARAKAWA ROTATED A Grid
1698  ELSE IF(maptype == 207)THEN !ARAKAWA A-GRID`
1699  igds(5) = 1
1700  ifield3len = 23
1701  ifield3 = 0
1702 !
1703  ifield3(1) = 6 !Earth assumed spherical with radius of 6,371,229.0m
1704  ifield3(8) = im !number of points along the x-axis
1705  ifield3(9) = jm !number of points along the y-axis
1706  ifield3(10) = 0 !Basic angle of the initial production domain
1707  if(.not.ldfgrd) then
1708  ifield3(11) = 0 !Subdivisions of basic angle used to define extreme lons & lats:missing
1709  else
1710  ifield3(11) = 45000000 !Subdivisions of basic angle used to define extreme lons & lats
1711  endif
1712  ifield3(12) = latstart_r !latitude of first grid point
1713  ifield3(13) = lonstart_r !longitude of first grid point
1714  ifield3(14) = 56 !Resolution and component flags
1715 ! Jili Dong change grid to earth relative (Matt Pyle)
1716  if(modelname=='FV3R') then
1717  ifield3(14) = 48 !Resolution and component flags
1718  endif
1719 
1720  ifield3(15) = latlast_r !latitude of last grid point
1721  ifield3(16) = lonlast_r !longitude of last grid point
1722  ifield3(17) = dxval
1723  ifield3(18) = dyval
1724  ifield3(19) = 64 !Scanning mode
1725  ifield3(20) = cenlat-90.*gdsdegr !Latitude of the southern pole of projection
1726  ifield3(21) = cenlon !Longitude of the southern pole of projection
1727  ifield3(22) = 0 !Angle of rotation of projection
1728  ifield3(23) = 2 !List of number of points along each meridian or parallel
1729 
1730 !
1731 !** Gaussian grid
1732  ELSE IF(maptype == 4 ) THEN !Gaussian grid
1733  igds(5) = 40
1734  ifield3len = 19
1735  ifield3 = 0
1736 !
1737  ifield3(1) = 6 !Earth assumed spherical with radius of 6,371,229.0m
1738  ifield3(8) = im
1739  ifield3(9) = jm
1740  ifield3(10) = 0
1741  ifield3(11) = 0
1742  ifield3(12) = latstart
1743  ifield3(13) = lonstart
1744  ifield3(14) = 48
1745  ifield3(15) = latlast
1746  ifield3(16) = lonlast
1747  ifield3(17) = nint(360./(im)*1000000.)
1748  ifield3(18) = nint(jm/2.0)
1749  if( latstart < latlast ) then
1750  ifield3(19) = 64 !for SN scan
1751  else
1752  ifield3(19) = 0 !for NS scan
1753  endif
1754 !
1755 !** Latlon grid
1756  ELSE IF(maptype == 0 ) THEN
1757  igds(5) = 0
1758  ifield3len = 19
1759  ifield3 = 0
1760 !
1761  ifield3(1) = 6 !Earth assumed spherical with radius of 6,371,229.0m
1762  ifield3(8) = im
1763  ifield3(9) = jm
1764  ifield3(10) = 0
1765  ifield3(11) = 0
1766  ifield3(12) = latstart
1767  ifield3(13) = lonstart
1768  ifield3(14) = 48
1769  ifield3(15) = latlast
1770  ifield3(16) = lonlast
1771 ! ifield3(17) = NINT(360./(IM)*1.0E6)
1772  ifield3(17) = abs(lonlast-lonstart)/(im-1)
1773 ! if(mod(jm,2) == 0 ) then
1774 ! ifield3(18) = NINT(180./JM*1.0E6)
1775 ! else
1776 ! ifield3(18) = NINT(180./(JM-1)*1.0E6)
1777  ifield3(18) = abs(latlast-latstart)/(jm-1)
1778 ! endif
1779  if( latstart < latlast ) then
1780  ifield3(19) = 64 !for SN scan
1781  else
1782  ifield3(19) = 0 !for NS scan
1783  endif
1784 
1785  ENDIF
1786 
1787 ! write(0,*)'igds=',igds,'igdstempl=',ifield3(1:ifield3len)
1788  end subroutine getgds
1789 !
1790 !-------------------------------------------------------------------------------------
1791 !
1792  end module grib2_module
xml_perl_data::paramset_t
Definition: xml_perl_data.f:48
xml_perl_data::param_t
Definition: xml_perl_data.f:13
grib2_module
Definition: grib2_module.f:1
xml_perl_data
Definition: xml_perl_data.f:1