module grib2_module !------------------------------------------------------------------------ ! ! This module generates grib2 messages and writes out the messages in ! parallel. ! ! program log: ! March, 2010 Jun Wang Initial code ! Jan, 2012 Jun Wang post available fields with grib2 description ! are defined in xml file !------------------------------------------------------------------------ use xml_data_post_t, only : param_t,paramset_t ! implicit none private ! ------------------------------------------------------------------------ ! !--- general grib2 info provided by post control file ! type param_t ! integer :: post_avblfldidx=-9999 ! character(len=80) :: shortname='' ! character(len=300) :: longname='' ! character(len=30) :: pdstmpl='' ! integer :: mass_windpoint=1 ! character(len=30) :: pname='' ! character(len=10) :: table_info='' ! character(len=20) :: stats_proc='' ! character(len=80) :: fixed_sfc1_type='' ! integer, dimension(:), pointer :: scale_fact_fixed_sfc1 => null() ! real, dimension(:), pointer :: level => null() ! character(len=80) :: fixed_sfc2_type='' ! integer, dimension(:), pointer :: scale_fact_fixed_sfc2 => null() ! real, dimension(:), pointer :: level2 => null() ! character(len=80) :: aerosol_type='' ! character(len=80) :: typ_intvl_size='' ! integer :: scale_fact_1st_size=0 ! real :: scale_val_1st_size=0. ! integer :: scale_fact_2nd_size=0 ! real :: scale_val_2nd_size=0. ! character(len=80) :: typ_intvl_wvlen='' ! integer :: scale_fact_1st_wvlen=0 ! real :: scale_val_1st_wvlen=0. ! integer :: scale_fact_2nd_wvlen=0 ! real :: scale_val_2nd_wvlen=0. ! real, dimension(:), pointer :: scale => null() ! integer :: stat_miss_val=0 ! integer :: leng_time_range_prev=0 ! integer :: time_inc_betwn_succ_fld=0 ! character(len=80) :: type_of_time_inc='' ! character(len=20) :: stat_unit_time_key_succ='' ! character(len=20) :: bit_map_flag='' ! end type param_t ! ! type paramset_t ! character(len=6) :: datset='' ! integer :: grid_num=255 ! character(len=20) :: sub_center='' ! character(len=20) :: version_no='' ! character(len=20) :: local_table_vers_no='' ! character(len=20) :: sigreftime='' ! character(len=20) :: prod_status='' ! character(len=20) :: data_type='' ! character(len=20) :: gen_proc_type='' ! character(len=30) :: time_range_unit='' ! character(len=50) :: orig_center='' ! character(len=30) :: gen_proc='' ! character(len=20) :: packing_method='' ! character(len=20) :: field_datatype='' ! character(len=20) :: comprs_type='' ! type(param_t), dimension(:), pointer :: param => null() ! end type paramset_t type(paramset_t),save :: pset ! !--- grib2 info related to a specific data file integer nrecout integer num_pset integer isec,hrs_obs_cutoff,min_obs_cutoff integer sec_intvl,stat_miss_val,time_inc_betwn_succ_fld character*80 type_of_time_inc,stat_unit_time_key_succ logical*1,allocatable :: bmap(:) integer ibm integer,allocatable :: mg(:) ! integer,parameter :: max_bytes=1000*1300000 integer,parameter :: MAX_NUMBIT=16 integer,parameter :: lugi=650 character*255 fl_nametbl,fl_gdss3 real(8) :: stime,stime1,stime2,etime,etime1 logical :: first_grbtbl ! public num_pset,pset,nrecout,gribit2,grib_info_init,first_grbtbl,grib_info_finalize real(8), EXTERNAL :: timef !------------------------------------------------------------------------------------- ! contains ! !------------------------------------------------------------------------------------- subroutine grib_info_init() ! !--- initialize general grib2 information and ! implicit none ! ! logical,intent(in) :: first_grbtbl ! !-- local variables integer ierr character(len=80) outfile character(len=10) envvar ! !-- set up pset ! !-- 1. pset is set up at READCNTRL_xml.f !-- initialize items of pset that are not set in xml file ! if(pset%grid_num==0) & pset%grid_num=218 if(trim(pset%sub_center)=='') & pset%sub_center="ncep_emc" if(trim(pset%version_no)=='') & pset%version_no="v2003" if(trim(pset%local_table_vers_no)=='') & pset%local_table_vers_no="local_table_no" if(trim(pset%sigreftime)=='') & pset%sigreftime="fcst" if(trim(pset%prod_status)=='') & pset%prod_status="oper_test" if(trim(pset%data_type)=='') & pset%data_type="fcst" if(trim(pset%orig_center)=='') & pset%orig_center="nws_ncep" if(trim(pset%time_range_unit)=='') & pset%time_range_unit="hour" if(trim(pset%gen_proc_type)=='') & pset%gen_proc_type="fcst" if(trim(pset%gen_proc)=='') & pset%gen_proc="gfs_avn" if(trim(pset%packing_method)=='') & pset%packing_method="jpeg" if(trim(pset%field_datatype)=='') & pset%field_datatype="flting_pnt" if(trim(pset%comprs_type)=='') & pset%comprs_type="lossless" ! !-- set up other grib2_info ! isec=0 hrs_obs_cutoff=0 ! applies to only obs min_obs_cutoff=0 ! applies to only obs sec_intvl=0 stat_miss_val=0 type_of_time_inc='same_start_time_fcst_fcst_time_inc' stat_unit_time_key_succ='missing' time_inc_betwn_succ_fld=0 ! !-- open fld name tble ! if(first_grbtbl) then fl_nametbl='params_grib2_tbl_new' call open_and_read_4dot2( fl_nametbl, ierr ) if ( ierr .ne. 0 ) then print*, 'Couldnt open table file - return code was ',ierr call mpi_abort() endif first_grbtbl=.false. endif ! !-- ! end subroutine grib_info_init !------------------------------------------------------------------------------------- !------------------------------------------------------------------------------------- ! subroutine grib_info_finalize ! !--- finalize grib2 information and close file ! implicit none ! !--- integer ierr call close_4dot2(ierr) ! end subroutine grib_info_finalize !------------------------------------------------------------------------------------- !------------------------------------------------------------------------------------- subroutine gribit2(post_fname) ! !------- use ctlblk_mod, only : im,jm,im_jm,num_procs,me,jsta,jend,ifhr,sdat,ihrst,imin, & mpi_comm_comp,ntlfld,fld_info,datapd,icnt,idsp implicit none ! include 'mpif.h' ! ! real,intent(in) :: data(im,1:jend-jsta+1,ntlfld) character(255),intent(in) :: post_fname ! !------- local variables integer i,j,k,n,nm,nprm,nlvl,fldlvl1,fldlvl2,cstart,cgrblen,ierr integer nf,nfpe,nmod integer fh, clength,lunout integer idisc,icatg,iparm,itblinfo,ntrange,leng_time_range_stat integer,allocatable :: nfld_pe(:),snfld_pe(:),enfld_pe(:) integer(4),allocatable :: isdsp(:),iscnt(:),ircnt(:),irdsp(:) integer status(MPI_STATUS_SIZE) integer(kind=MPI_OFFSET_KIND) idisp integer,allocatable :: jsta_pe(:),jend_pe(:) integer,allocatable :: grbmsglen(:) real,allocatable :: datafld(:,:) real,allocatable :: datafldtmp(:) ! character(1) cgrib(max_bytes) ! ! !---------------- code starts here -------------------------- ! ! !******* part 1 resitribute data ******** ! !--- calculate # of fields on each processor ! nf=ntlfld/num_procs nfpe=nf+1 nmod=mod(ntlfld,num_procs) ! print *,'ntlfld=',ntlfld,'nf=',nf,'nmod=',nmod allocate(snfld_pe(num_procs),enfld_pe(num_procs),nfld_pe(num_procs)) do n=1,num_procs if(n-11.and.size(pset%param(nprm)%level2)>=nlvl) then scaled_val_fixed_sfc2=nint(pset%param(nprm)%level2(nlvl)) else if(size(pset%param(nprm)%level2)==1) then scaled_val_fixed_sfc2=nint(pset%param(nprm)%level2(1)) else scaled_val_fixed_sfc2=0 endif if(size(pset%param(nprm)%scale_fact_fixed_sfc2)>1 .and. & size(pset%param(nprm)%scale_fact_fixed_sfc2)>=nlvl) then scale_fct_fixed_sfc2=pset%param(nprm)%scale_fact_fixed_sfc2(nlvl) else if(size(pset%param(nprm)%scale_fact_fixed_sfc2)==1) then scale_fct_fixed_sfc2=pset%param(nprm)%scale_fact_fixed_sfc2(1) else scale_fct_fixed_sfc2=0 endif ! print *,'bf g2sec4_temp0,ipdstmpl=',trim(pset%param(nprm)%pdstmpl),'fixed_sfc_type=', & ! pset%param(nprm)%fixed_sfc1_type,'scale_fct_fixed_sfc1=', & ! scale_fct_fixed_sfc1,'scaled_val_fixed_sfc1=',scaled_val_fixed_sfc1, & ! 'sfc2_type=',trim(pset%param(nprm)%fixed_sfc2_type),scale_fct_fixed_sfc2, & ! scaled_val_fixed_sfc2 if(trim(pset%param(nprm)%pdstmpl)=='tmpl4_0') then ipdsnum=0 ipdstmpllen=ipdstmp4_0len call g2sec4_temp0(icatg,iparm,pset%gen_proc_type, & pset%gen_proc,hrs_obs_cutoff,min_obs_cutoff, & pset%time_range_unit,ifhr, & pset%param(nprm)%fixed_sfc1_type, & scale_fct_fixed_sfc1, & scaled_val_fixed_sfc1, & fixed_sfc2_type, & scale_fct_fixed_sfc2, & scaled_val_fixed_sfc2, & ipdstmpl(1:ipdstmpllen)) ! print *,'aft g2sec4_temp0,ipdstmpl0=',ipdstmpl(1:ipdstmp4_0len) elseif(trim(pset%param(nprm)%pdstmpl)=='tmpl4_8') then ! ipdsnum=8 ipdstmpllen=ipdstmp4_8len call g2sec4_temp8(icatg,iparm,pset%gen_proc_type, & pset%gen_proc,hrs_obs_cutoff,min_obs_cutoff, & pset%time_range_unit,ifhr-tinvstat, & pset%param(nprm)%fixed_sfc1_type, & scale_fct_fixed_sfc1, & scaled_val_fixed_sfc1, & pset%param(nprm)%fixed_sfc2_type, & scale_fct_fixed_sfc2, & scaled_val_fixed_sfc2, & idat(3),idat(1),idat(2),idat(4),idat(5), & sec_intvl,ntrange,stat_miss_val, & pset%param(nprm)%stats_proc,type_of_time_inc, & pset%time_range_unit, tinvstat, & stat_unit_time_key_succ,time_inc_betwn_succ_fld, & ipdstmpl(1:ipdstmpllen)) print *,'aft g2sec4_temp8,ipdstmpl8=',ipdstmpl(1:ipdstmp4_8len) elseif(trim(pset%param(nprm)%pdstmpl)=='tmpl4_44') then ! ipdsnum=44 ipdstmpllen=ipdstmp4_44len call g2sec4_temp44(icatg,iparm,pset%param(nprm)%aerosol_type, & pset%param(nprm)%typ_intvl_size, & pset%param(nprm)%scale_fact_1st_size, & pset%param(nprm)%scale_val_1st_size, & pset%param(nprm)%scale_fact_2nd_size, & pset%param(nprm)%scale_val_2nd_size, & pset%gen_proc_type, & pset%gen_proc,hrs_obs_cutoff,min_obs_cutoff, & pset%time_range_unit,ifhr, & pset%param(nprm)%fixed_sfc1_type, & scale_fct_fixed_sfc1, & scaled_val_fixed_sfc1, & pset%param(nprm)%fixed_sfc2_type, & scale_fct_fixed_sfc2, & scaled_val_fixed_sfc2, & ipdstmpl(1:ipdstmpllen)) ! print *,'aft g2sec4_temp44,ipdstmpl44=',ipdstmpl(1:ipdstmp4_44len),'ipdsnum=',ipdsnum elseif(trim(pset%param(nprm)%pdstmpl)=='tmpl4_48') then ! ipdsnum=48 ipdstmpllen=ipdstmp4_48len call g2sec4_temp48(icatg,iparm,pset%param(nprm)%aerosol_type, & pset%param(nprm)%typ_intvl_size, & pset%param(nprm)%scale_fact_1st_size, & pset%param(nprm)%scale_val_1st_size, & pset%param(nprm)%scale_fact_2nd_size, & pset%param(nprm)%scale_val_2nd_size, & pset%param(nprm)%typ_intvl_wvlen, & pset%param(nprm)%scale_fact_1st_wvlen, & pset%param(nprm)%scale_val_1st_wvlen, & pset%param(nprm)%scale_fact_2nd_wvlen, & pset%param(nprm)%scale_val_2nd_wvlen, & pset%gen_proc_type, & pset%gen_proc,hrs_obs_cutoff,min_obs_cutoff, & pset%time_range_unit,ifhr, & pset%param(nprm)%fixed_sfc1_type, & scale_fct_fixed_sfc1, & scaled_val_fixed_sfc1, & pset%param(nprm)%fixed_sfc2_type, & scale_fct_fixed_sfc2, & scaled_val_fixed_sfc2, & ipdstmpl(1:ipdstmpllen)) ! print *,'aft g2sec4_temp48,name=',trim(pset%param(nprm)%shortname),& ! 'ipdstmpl48=',ipdstmpl(1:ipdstmp4_48len) endif ! !---------- ! idrstmpl array is the output from g2sec5 ! call get_g2_sec5packingmethod(pset%packing_method,idrsnum,ierr) if(maxval(datafld1)==minval(datafld1))then idrsnum=0 print*,' changing to simple packing for constant fields' end if print *,'aft g2sec5,packingmethod=',pset%packing_method,'idrsnum=',idrsnum, & 'data=',maxval(datafld1),minval(datafld1) ! !*** set number of bits, and binary scale ! ibmap=255 bmap=.true. if(any(datafld1>=SPVAL))then ibmap=0 where(abs(datafld1)>=SPVAL)bmap=.false. endif ! if(size(pset%param(nprm)%level)==size(pset%param(nprm)%scale)) then if(size(pset%param(nprm)%level)==1) nlvl=1 if(nlvl/=0) then fldscl=nint(pset%param(nprm)%scale(nlvl)) else fldscl=0 endif else if (size(pset%param(nprm)%scale)==1) then fldscl=nint(pset%param(nprm)%scale(1)) endif ! call g2getbits(ibmap,fldscl,size(datafld1),bmap,datafld1,ibin_scl,idec_scl,inumbits) ! print *,'idec_scl=',idec_scl,'ibin_scl=',ibin_scl,'number_bits=',inumbits if( idrsnum==40 ) then idrstmplen=idrstmp5_40len call g2sec5_temp40(idec_scl,ibin_scl,inumbits,pset%comprs_type,idrstmpl(1:idrstmplen)) elseif( idrsnum==2 ) then idrstmplen=idrstmp5_2len call g2sec5_temp2(idec_scl,ibin_scl,idrstmpl(1:idrstmplen)) elseif( idrsnum==3 ) then idrstmplen=idrstmp5_3len call g2sec5_temp3(idec_scl,ibin_scl,pset%order_of_sptdiff,idrstmpl(1:idrstmplen)) elseif( idrsnum==0 ) then idrstmplen=idrstmp5_0len call g2sec5_temp0(idec_scl,ibin_scl,inumbits,idrstmpl(1:idrstmplen)) endif ! !---------------------------------------------------------------------------------------- ! Define all the required inputs like ibmap, numcoord, coordlist etc externally ! in the module prior to calling the addfield routine ! Again hide the addfield routine from the user ! ! print *,'before addfield, data=',maxval(datafld1),minval(datafld1),'ibmap=',ibmap, & ! 'max_bytes=',max_bytes,'ipdsnum=',ipdsnum,'ipdstmpllen=',ipdstmpllen,'ipdstmpl=',ipdstmpl(1:ipdstmpllen), & ! 'coordlist=',coordlist,'numcoord=',numcoord,'idrsnum=',idrsnum,'idrstmpl=',idrstmpl, & ! 'idrstmplen=',idrstmplen,'im_jm=',im_jm call addfield(cgrib,max_bytes,ipdsnum,ipdstmpl(1:ipdstmpllen), & ipdstmpllen,coordlist,numcoord,idrsnum,idrstmpl, & idrstmplen ,datafld1,im_jm,ibmap,bmap,ierr) ! !--------------------------------------------------------------------------------------- ! Finalize GRIB message after all grids and fields have been added by adding the end ! section "7777" ! Again hide the gribend routine from the user ! call gribend(cgrib,max_bytes,lengrib,ierr) ! !------- end subroutine gengrb2msg ! !-------------------------------------------------------------------------------------- ! subroutine g2sec3tmpl40(nx,nY,lat1,lon1,lat2,lon2,lad,ds1,len3,igds,ifield3) implicit none ! integer(4),intent(inout) :: nx,ny,lat1,lon1,lat2,lon2,lad,ds1,len3 integer(4),intent(inout) :: ifield3(len3),igds(5) ! nx=1152 ny=576 lat1=89761000 !lat of 1st grd pt in micro-deg lon1=0 !east-long of 1st grd pt in micro-deg lat2=-89761000 !lat of last grd pt in micro-deg lon2=359687000 !east-long of last grd pt in micro-deg lad=313000 !lat at which projection intersects earth ds1=288 !grid spacing in x and y ! igds=0 igds(2)=nx*ny igds(5)=40 ! ifield3(1) = 6 !Earth assumed spherical with radius of 6,371,229.0m ifield3(2) = 0 ifield3(3) = 0 ifield3(4) = 0 ifield3(5) = 0 ifield3(6) = 0 ifield3(7) = 0 ifield3(8) = nx ifield3(9) = ny ifield3(10) = 0 ifield3(11) = 0 ifield3(12) = lat1 ifield3(13) = lon1 ifield3(14) = 48 ifield3(15) = lat2 ifield3(16) = lon2 ifield3(17) = lad ifield3(18) = ds1 ifield3(19) = 0 return end subroutine g2sec3tmpl40 ! !------------------------------------------------------------------------------------- ! subroutine g2getbits(ibm,scl,len,bmap,g,ibs,ids,nbits) !$$$ ! This subroutine is changed from w3 lib getbit to compute the total number of bits, ! The argument list is modified to have ibm,scl,len,bmap,g,ibs,ids,nbits ! ! Progrma log: ! Jun Wang Apr, 2010 ! ! INPUT ! ibm: integer, bitmap flag (grib2 table 6.0) ! scl: real, significant digits,OR binary precision if < 0 ! len: integer, field and bitmap length ! bmap: logical(len), bitmap (.true.: keep, bitmap (.true.: keep, .false. skip) ! fld: real(len), datafield ! OUTPUT ! ibs: integer, binary scale factor ! ids: integer, decimal scale factor ! nbits: integer, number of bits to pack ! IMPLICIT NONE ! INTEGER,INTENT(IN) :: IBM,LEN LOGICAL*1,INTENT(IN) :: BMAP(LEN) REAL,INTENT(IN) :: scl,G(LEN) INTEGER,INTENT(OUT) :: IBS,IDS,NBITS ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - INTEGER,PARAMETER :: MXBIT=16 ! ! NATURAL LOGARITHM OF 2 AND 0.5 PLUS NOMINAL SAFE EPSILON real,PARAMETER :: ALOG2=0.69314718056,HPEPS=0.500001 ! !local vars INTEGER :: I,I1,icnt,ipo,le,irange REAL :: GROUND,GMIN,GMAX,s,rmin,rmax,range,rr,rng2,po,rln2 ! DATA rln2/0.69314718/ ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! ROUND FIELD AND DETERMINE EXTREMES WHERE BITMAP IS ON IF(IBM == 255) THEN GMAX = G(1) GMIN = G(1) DO I=2,LEN GMAX = MAX(GMAX,G(I)) GMIN = MIN(GMIN,G(I)) ENDDO ELSE do i1=1,len if (bmap(i1)) exit enddo ! I1 = 1 ! DO WHILE(I1 <= LEN .AND. .not. BMAP(I1)) ! I1=I1+1 ! ENDDO IF(I1 <= LEN) THEN GMAX = G(I1) GMIN = G(I1) DO I=I1+1,LEN IF(BMAP(I)) THEN GMAX = MAX(GMAX,G(I)) GMIN = MIN(GMIN,G(I)) ENDIF ENDDO ELSE GMAX = 0. GMIN = 0. ENDIF ENDIF ! write(0,*)' GMIN=',GMIN,' GMAX=',GMAX ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! COMPUTE NUMBER OF BITS icnt = 0 ibs = 0 ids = 0 range = GMAX - GMIN ! IF ( range .le. 0.00 ) THEN IF ( range .le. 1.e-30 ) THEN nbits = 8 return END IF !* IF ( scl .eq. 0.0 ) THEN nbits = 8 RETURN ELSE IF ( scl > 0.0 ) THEN ipo = INT (ALOG10 ( range )) !jw: if range is smaller than computer precision, set nbits=8 if(ipo<0.and.ipo+scl<-20) then print *,'for small range,ipo=',ipo,'ipo+scl=',ipo+scl,'scl=',scl nbits=8 return endif IF ( range .lt. 1.00 ) ipo = ipo - 1 po = float(ipo) - scl + 1. ids = - INT ( po ) rr = range * 10. ** ( -po ) nbits = INT ( ALOG ( rr ) / rln2 ) + 1 ELSE ibs = -NINT ( -scl ) rng2 = range * 2. ** (-ibs) nbits = INT ( ALOG ( rng2 ) / rln2 ) + 1 END IF ! write(0,*)'in g2getnits,ibs=',ibs,'ids=',ids,'nbits=',nbits,'range=',range !* IF(nbits <= 0) THEN nbits = 0 IF(ABS(GMIN) >= 1.) THEN ids = -int(alog10(abs(gmin))) ELSE IF (ABS(GMIN) < 1.0.AND.ABS(GMIN) > 0.0) THEN ids = -int(alog10(abs(gmin)))+1 ELSE ids = 0 ENDIF ENDIF nbits = min(nbits,MXBIT) ! write(0,*)'in g2getnits ibs=',ibs,'ids=',ids,'nbits=',nbits ! IF ( scl > 0.0 ) THEN s=10.0 ** ids IF(IBM == 255) THEN GROUND = G(1)*s GMAX = GROUND GMIN = GROUND DO I=2,LEN GMAX = MAX(GMAX,G(I)*s) GMIN = MIN(GMIN,G(I)*s) ENDDO ELSE do i1=1,len if (bmap(i1)) exit enddo ! I1=1 ! DO WHILE(I1.LE.LEN.AND..not.BMAP(I1)) ! I1=I1+1 ! ENDDO IF(I1 <= LEN) THEN GROUND = G(I1)*s GMAX = GROUND GMIN = GROUND DO I=I1+1,LEN IF(BMAP(I)) THEN GMAX = MAX(GMAX,G(I)*S) GMIN = MIN(GMIN,G(I)*S) ENDIF ENDDO ELSE GMAX = 0. GMIN = 0. ENDIF ENDIF range = GMAX-GMIN if(GMAX == GMIN) then ibs = 0 else ibs = nint(alog(range/(2.**NBITS-0.5))/ALOG2+HPEPS) endif ! endif ! write(0,*)'in g2getnits,2ibs=',ibs,'ids=',ids,'nbits=',nbits,'range=',& ! range, 'scl=',scl,'data=',maxval(g),minval(g) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - RETURN END subroutine g2getbits ! !------------------------------------------------------------------------------------- ! subroutine getgds(ldfgrd,len3,ifield3len,igds,ifield3) ! !***** set up gds kpds to call Boi's code ! use CTLBLK_mod, only : im,jm use gridspec_mod, only: DXVAL,DYVAL,CENLAT,CENLON,LATSTART,LONSTART,LATLAST, & & LONLAST,MAPTYPE,STANDLON,latstartv,cenlatv,lonstartv, & cenlonv,TRUELAT1,TRUELAT2 ! implicit none ! logical, intent(in) :: ldfgrd integer(4),intent(in) :: len3 integer(4),intent(out) :: ifield3len integer(4),intent(inout) :: ifield3(len3),igds(5) integer :: g2_latstart, g2_lonstart integer :: g2_TRUELAT1, g2_TRUELAT2 integer :: g2_STANDLON integer :: g2_DXVAL, g2_DYVAL ! print *,'in getgds, im=',im,'jm=',jm,'latstart=',latstart,'lonsstart=',lonstart,'maptyp=',maptype ! print *,'in getgds bf , im=',im,'jm=',jm,'latstart=',latstart,'lonsstart=',lonstart ! print *,'in getgds bf, TRUELAT1=',TRUELAT1,'TRUELAT2=',TRUELAT2,'STANDLON=',STANDLON ! print *,'in getgds bf, latstartv=',latstartv,'cenlatv=',cenlatv,'lonstartv=',lonstartv ! print *,'in getgds bf, DXVAL= ', DXVAL, 'DYVAL= ', DYVAL ! set to 0 g2_latstart = 0 g2_lonstart = 0 g2_TRUELAT1 = 0 g2_TRUELAT2 = 0 g2_STANDLON = 0 g2_DXVAL = 0 g2_DYVAL = 0 g2_latstart = latstart * 1000 g2_lonstart = (lonstart + 360000) * 1000 g2_TRUELAT1 = TRUELAT1 * 1000 g2_TRUELAT2 = TRUELAT2 * 1000 g2_STANDLON = (STANDLON + 360000) * 1000 g2_DXVAL = DXVAL * 1000 g2_DYVAL = DYVAL * 1000 ! print *,'in getgds af , im=',im,'jm=',jm,'latstart=',latstart,'lonsstart=',lonstart ! print *,'in getgds af, TRUELAT1=',TRUELAT1,'TRUELAT2=',TRUELAT2,'STANDLON=',STANDLON ! print *,'in getgds af, latstartv=',latstartv,'cenlatv=',cenlatv,'lonstartv=',lonstartv ! print *,'in getgds af, DXVAL= ', DXVAL, 'DYVAL= ', DYVAL ! !** set up igds igds(1) = 0 !Source of grid definition (see Code Table 3.0) igds(2) = im*jm !Number of grid points in the defined grid. igds(3) = 0 !Number of octets needed for each additional grid points definition igds(4) = 0 !Interpretation of list for optional points definition (Code Table 3.11) ! !** define grid template 3 IF(MAPTYPE == 1) THEN !LAmbert Conformal igds(5) = 30 !Lambert conformal ifield3len = 22 ifield3 = 0 ! ifield3(1) = 6 !Earth assumed spherical with radius of 6,371,229.0m ifield3(8) = im !number of points along the x-axis ifield3(9) = jm !number of points along the y-axis ifield3(10) = g2_latstart !latitude of first grid point ifield3(11) = g2_lonstart !longitude of first grid point ifield3(12) = 8 !Resolution and component flags ifield3(13) = g2_TRUELAT1 ifield3(14) = g2_STANDLON !longitude of meridian parallel to y-axis along which latitude increases ifield3(15) = g2_DXVAL ifield3(16) = g2_DYVAL IF(TRUELAT1>0)then ifield3(17) = 0 else ifield3(17) =128 !for southern hemisphere end if ifield3(18) = 64 ifield3(19) = g2_TRUELAT1 !first latitude from the pole at which the secant cone cuts the sphere ifield3(20) = g2_TRUELAT2 !second latitude from the pole at which the secant cone cuts the sphere !** Polar stereographic ELSE IF(MAPTYPE == 2)THEN !Polar stereographic igds(5) = 20 ifield3len = 22 ifield3 = 0 ! ifield3(1) = 6 !Earth assumed spherical with radius of 6,371,229.0m ifield3(8) = im !number of points along the x-axis ifield3(9) = jm !number of points along the y-axis ifield3(10) = latstart !latitude of first grid point ifield3(11) = lonstart !longitude of first grid point ifield3(12) = 8 !Resolution and component flags ifield3(13) = TRUELAT1 ifield3(14) = STANDLON !longitude of meridian parallel to y-axis along which latitude increases ifield3(15) = DXVAL ifield3(16) = DYVAL ifield3(17) = 0 ifield3(18) = 64 ! !** Mercator ELSE IF(MAPTYPE.EQ.3)THEN !Mercator igds(5)=10 ifield3len=22 ifield3=0 ! ifield3(1) = 6 !Earth assumed spherical with radius of 6,371,229.0m ifield3(8) = im !number of points along the x-axis ifield3(9) = jm !number of points along the y-axis ifield3(10) = latstart !latitude of first grid point ifield3(11) = lonstart !longitude of first grid point ifield3(12) = 8 !Resolution and component flags ifield3(13) = TRUELAT1 !latitude(s) at which the Mercator projection intersects the Earth ifield3(14) = latlast !latitude of last grid point ifield3(15) = lonlast !longitude of last grid point ifield3(16) = 64 !Scanning mode ifield3(17) = 0 !Orientation of the grid, angle between i direction on the map and Equator ifield3(18) = DXVAL ifield3(19) = DYVAL ! !** ARAKAWA STAGGERED E-GRID ELSE IF(MAPTYPE == 203)THEN !ARAKAWA STAGGERED E-GRID` igds(5) = 32768 ifield3len = 22 ifield3 = 0 ! ifield3(1) = 6 !Earth assumed spherical with radius of 6,371,229.0m ifield3(8) = im !number of points along the x-axis ifield3(9) = jm !number of points along the y-axis ifield3(10) = 0 !Basic angle of the initial production domain if(.not.ldfgrd) then ifield3(11) = 0 !Subdivisions of basic angle used to define extreme lons & lats:missing else ifield3(11) = 45000000 !Subdivisions of basic angle used to define extreme lons & lats endif ifield3(12) = latstart !latitude of first grid point ifield3(13) = lonstart !longitude of first grid point ifield3(14) = 0 !Resolution and component flags ifield3(15) = CENLAT !center latitude of grid point ifield3(16) = CENLON !Center longitude of grid point ifield3(17) = DXVAL ifield3(18) = DYVAL ifield3(19) = 64 !Scanning mode ! !** ARAKAWA STAGGERED non-E-GRID ELSE IF(MAPTYPE == 205)THEN !ARAKAWA STAGGERED E-GRID` igds(5) = 32769 ifield3len = 21 ifield3 = 0 ! ifield3(1) = 6 !Earth assumed spherical with radius of 6,371,229.0m ifield3(8) = im !number of points along the x-axis ifield3(9) = jm !number of points along the y-axis ifield3(10) = 0 !Basic angle of the initial production domain if(.not.ldfgrd) then ifield3(11) = 0 !Subdivisions of basic angle used to define extreme lons & lats:missing else ifield3(11) = 45000000 !Subdivisions of basic angle used to define extreme lons & lats endif ifield3(12) = latstart !latitude of first grid point ifield3(13) = lonstart !longitude of first grid point ifield3(14) = 56 !Resolution and component flags ifield3(15) = CENLAT !center latitude of grid point ifield3(16) = CENLON !Center longitude of grid point ifield3(17) = DXVAL ifield3(18) = DYVAL ifield3(19) = 64 !Scanning mode ifield3(20) = latlast !Scanning mode ifield3(21) = lonlast !Scanning mode ! !** Gaussian grid ELSE IF(MAPTYPE == 4 ) THEN !Gaussian grid igds(5) = 40 ifield3len = 19 ifield3 = 0 ! ifield3(1) = 6 !Earth assumed spherical with radius of 6,371,229.0m ifield3(8) = im ifield3(9) = jm ifield3(10) = 0 ifield3(11) = 0 ifield3(12) = latstart ifield3(13) = lonstart ifield3(14) = 48 ifield3(15) = latlast ifield3(16) = lonlast ifield3(17) = NINT(360./(IM)*1000000.) ifield3(18) = NINT(JM/2.0) ! !** Latlon grid ELSE IF(MAPTYPE == 0 ) THEN igds(5) = 0 ifield3len = 19 ifield3 = 0 ! ifield3(1) = 6 !Earth assumed spherical with radius of 6,371,229.0m ifield3(8) = im ifield3(9) = jm ifield3(10) = 0 ifield3(11) = 0 ifield3(12) = latstart ifield3(13) = lonstart ifield3(14) = 48 ifield3(15) = latlast ifield3(16) = lonlast ifield3(17) = NINT(180./(JM-1)*1.0E6) ifield3(18) = NINT(360./(IM)*1.0E6) ifield3(19) = 0 ENDIF write(0,*)'igds=',igds,'igdstempl=',ifield3(1:ifield3len) end subroutine getgds ! !------------------------------------------------------------------------------------- ! end module grib2_module