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