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