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