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