运行Noahmp在create_UV.perl时出错

gaochao123 · 发表于 2020-5-3 19:32:07

登录后查看更多精彩内容~

您需要登录才可以下载或查看，没有帐号？立即注册

x

在运行create_UV.perl时出了这个错，/mnt/d/extracted_Gldas/Wind/GLDAS_Wind.2019090100.grb

Program received signal SIGSEGV: Segmentation fault - invalid memory reference.

Backtrace for this error:
#0  0x7f4e7b83b2ed in ???
#1  0x7f4e7b83a503 in ???
#2  0x7f4e7ae9ef1f in ???
#3  0x7f4e7c004532 in w3fi74_
      at /usr/local/w3lib-2.0.6/w3fi74.f:422
#4  0x7f4e7c002af7 in w3fi72_
      at /usr/local/w3lib-2.0.6/w3fi72.f:255
#5  0x7f4e7c001778 in putgb_
      at /usr/local/w3lib-2.0.6/putgb.f:197
#6  0x7f4e7c001320 in ???
#7  0x7f4e7c00145c in ???
#8  0x7f4e7ae81b96 in ???
#9  0x7f4e7c001019 in ???
#10  0xffffffffffffffff in ???

所有的风速的grb文件都是一样的错误。
我的create_UV.perl脚本：
#!/usr/bin/perl

$filename = "create_UV";
system ("mpif90 -o $filename $filename.f90 -L/usr/local/w3lib-2.0.6 -lw3");

@nums = ("00","01","02","03","04","05","06","07","08","09",
      "10","11","12","13","14","15","16","17","18","19",
   "20","21","22","23","24","25","26","27","28","29",
   "30","31");

@yrs = ("19");

$day_start = 244;
$day_end = 250;

@hrs = ("00","03","06","09","12","15","18","21");

@noleap_days = (0,31,59,90,120,151,181,212,243,273,304,334,365);
@leap_days = (0,31,60,91,121,152,182,213,244,274,305,335,366);

$data_dir = "/mnt/d/extracted_Gldas";
$results_dir = "/mnt/d/extracted_Gldas/u_v_wind";

for $yy (@yrs)
{

# This will be the jday time loop

for($julday=$day_start;$julday<=$day_end;$julday++)

{

# This little section finds the text month and day

@modays = @noleap_days;
if($yy == "00" || $yy == "04" || $yy == "08" || $yy == "12" || $yy == "16" || $yy == "20" ) {@modays = @leap_days}

for($mo=1;$mo<=12;$mo++)
  {
if($julday>$modays[$mo-1] && $julday<=$modays[$mo])
   {
   $mon = $mo;
   $day = $julday - $modays[$mo-1];
   }
  }

for $hr (@hrs)
{

  $file_in  = "$data_dir/Wind/GLDAS_Wind.20$yy$nums[$mon]$nums[$day]$hr.grb";
  $file1_out = "$results_dir/U/GLDAS_U.20$yy$nums[$mon]$nums[$day]$hr.grb";
  $file2_out = "$results_dir/V/GLDAS_V.20$yy$nums[$mon]$nums[$day]$hr.grb";

  print ("$file_in \n");
  system ("./$filename $file_in $file1_out $file2_out");
}
}
} # End of outer time loop

system("rm $filename");

w3fi74.f是在第422行：ICOMP = ISHFT(ICOMP,3)；
w3fi72.f是在第255行出错：CALL W3FI74(IGDS,ICOMP,GDS,LENGDS,NPTS,IGERR)（应该跟w3fi74.f的错误有关）
putgb.f错误：CALL W3FI72(0,FR,0,NBIT,0,IPDS,PDS,       //第195行
   &          1,255,IGDS,0,0,IBM,KF,IBDS,    //第196行
   &          KFO,GRIB,LGRIB,IRET)             //第197行

希望大家找找问题在哪里，不甚感谢
另在评论区附上我的w3fi74.f，w3fi72.f，putgb.f的所有内容

gaochao123 · 发表于 2020-5-3 19:34:22

本帖最后由 gaochao123 于 2020-5-3 19:41 编辑

这个是putgb.f
C-----------------------------------------------------------------------
   SUBROUTINE PUTGB(LUGB,KF,KPDS,KGDS,LB,F,IRET)
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM: PUTGB       PACKS AND WRITES A GRIB MESSAGE
C PRGMMR: IREDELL       ORG: W/NMC23    DATE: 94-04-01
C
C ABSTRACT: PACK AND WRITE A GRIB MESSAGE.
C THIS SUBPROGRAM IS NEARLY THE INVERSE OF GETGB.
C
C PROGRAM HISTORY LOG:
C 94-04-01  IREDELL
C 95-10-31  IREDELL    REMOVED SAVES AND PRINTS
C 09-10-15  GAYNO    INCREASED MAXBIT FROM 16 TO 32
C
C USAGE: CALL PUTGB(LUGB,KF,KPDS,KGDS,LB,F,IRET)
C INPUT ARGUMENTS:
C    LUGB       INTEGER UNIT OF THE UNBLOCKED GRIB DATA FILE
C    KF          INTEGER NUMBER OF DATA POINTS
C    KPDS       INTEGER (200) PDS PARAMETERS
C       (1) - ID OF CENTER
C       (2) - GENERATING PROCESS ID NUMBER
C       (3) - GRID DEFINITION
C       (4) - GDS/BMS FLAG (RIGHT ADJ COPY OF OCTET 8)
C       (5) - INDICATOR OF PARAMETER
C       (6) - TYPE OF LEVEL
C       (7) - HEIGHT/PRESSURE , ETC OF LEVEL
C       (8) - YEAR INCLUDING (CENTURY-1)
C       (9) - MONTH OF YEAR
C       (10)  - DAY OF MONTH
C       (11)  - HOUR OF DAY
C       (12)  - MINUTE OF HOUR
C       (13)  - INDICATOR OF FORECAST TIME UNIT
C       (14)  - TIME RANGE 1
C       (15)  - TIME RANGE 2
C       (16)  - TIME RANGE FLAG
C       (17)  - NUMBER INCLUDED IN AVERAGE
C       (18)  - VERSION NR OF GRIB SPECIFICATION
C       (19)  - VERSION NR OF PARAMETER TABLE
C       (20)  - NR MISSING FROM AVERAGE/ACCUMULATION
C       (21)  - CENTURY OF REFERENCE TIME OF DATA
C       (22)  - UNITS DECIMAL SCALE FACTOR
C       (23)  - SUBCENTER NUMBER
C       (24)  - PDS BYTE 29, FOR NMC ENSEMBLE PRODUCTS
C                128 IF FORECAST FIELD ERROR
C                64 IF BIAS CORRECTED FCST FIELD
C                32 IF SMOOTHED FIELD
C                WARNING: CAN BE COMBINATION OF MORE THAN 1
C       (25)  - PDS BYTE 30, NOT USED
C    KGDS       INTEGER (200) GDS PARAMETERS
C       (1) - DATA REPRESENTATION TYPE
C       (19)  - NUMBER OF VERTICAL COORDINATE PARAMETERS
C       (20)  - OCTET NUMBER OF THE LIST OF VERTICAL COORDINATE
C                PARAMETERS
C                OR
C                OCTET NUMBER OF THE LIST OF NUMBERS OF POINTS
C                IN EACH ROW
C                OR
C                255 IF NEITHER ARE PRESENT
C       (21)  - FOR GRIDS WITH PL, NUMBER OF POINTS IN GRID
C       (22)  - NUMBER OF WORDS IN EACH ROW
C    LATITUDE/LONGITUDE GRIDS
C       (2) - N(I) NR POINTS ON LATITUDE CIRCLE
C       (3) - N(J) NR POINTS ON LONGITUDE MERIDIAN
C       (4) - LA(1) LATITUDE OF ORIGIN
C       (5) - LO(1) LONGITUDE OF ORIGIN
C       (6) - RESOLUTION FLAG (RIGHT ADJ COPY OF OCTET 17)
C       (7) - LA(2) LATITUDE OF EXTREME POINT
C       (8) - LO(2) LONGITUDE OF EXTREME POINT
C       (9) - DI LONGITUDINAL DIRECTION OF INCREMENT
C       (10)  - DJ LATITUDINAL DIRECTION INCREMENT
C       (11)  - SCANNING MODE FLAG (RIGHT ADJ COPY OF OCTET 28)
C    GAUSSIAN  GRIDS
C       (2) - N(I) NR POINTS ON LATITUDE CIRCLE
C       (3) - N(J) NR POINTS ON LONGITUDE MERIDIAN
C       (4) - LA(1) LATITUDE OF ORIGIN
C       (5) - LO(1) LONGITUDE OF ORIGIN
C       (6) - RESOLUTION FLAG  (RIGHT ADJ COPY OF OCTET 17)
C       (7) - LA(2) LATITUDE OF EXTREME POINT
C       (8) - LO(2) LONGITUDE OF EXTREME POINT
C       (9) - DI LONGITUDINAL DIRECTION OF INCREMENT
C       (10)  - N - NR OF CIRCLES POLE TO EQUATOR
C       (11)  - SCANNING MODE FLAG (RIGHT ADJ COPY OF OCTET 28)
C       (12)  - NV - NR OF VERT COORD PARAMETERS
C       (13)  - PV - OCTET NR OF LIST OF VERT COORD PARAMETERS
C                            OR
C                PL - LOCATION OF THE LIST OF NUMBERS OF POINTS IN
C                      EACH ROW (IF NO VERT COORD PARAMETERS
C                      ARE PRESENT
C                            OR
C                255 IF NEITHER ARE PRESENT
C    POLAR STEREOGRAPHIC GRIDS
C       (2) - N(I) NR POINTS ALONG LAT CIRCLE
C       (3) - N(J) NR POINTS ALONG LON CIRCLE
C       (4) - LA(1) LATITUDE OF ORIGIN
C       (5) - LO(1) LONGITUDE OF ORIGIN
C       (6) - RESOLUTION FLAG  (RIGHT ADJ COPY OF OCTET 17)
C       (7) - LOV GRID ORIENTATION
C       (8) - DX - X DIRECTION INCREMENT
C       (9) - DY - Y DIRECTION INCREMENT
C       (10)  - PROJECTION CENTER FLAG
C       (11)  - SCANNING MODE (RIGHT ADJ COPY OF OCTET 28)
C    SPHERICAL HARMONIC COEFFICIENTS
C       (2) - J PENTAGONAL RESOLUTION PARAMETER
C       (3) - K    "       "       "
C       (4) - M    "       "       "
C       (5) - REPRESENTATION TYPE
C       (6) - COEFFICIENT STORAGE MODE
C    MERCATOR GRIDS
C       (2) - N(I) NR POINTS ON LATITUDE CIRCLE
C       (3) - N(J) NR POINTS ON LONGITUDE MERIDIAN
C       (4) - LA(1) LATITUDE OF ORIGIN
C       (5) - LO(1) LONGITUDE OF ORIGIN
C       (6) - RESOLUTION FLAG (RIGHT ADJ COPY OF OCTET 17)
C       (7) - LA(2) LATITUDE OF LAST GRID POINT
C       (8) - LO(2) LONGITUDE OF LAST GRID POINT
C       (9) - LATIT - LATITUDE OF PROJECTION INTERSECTION
C       (10)  - RESERVED
C       (11)  - SCANNING MODE FLAG (RIGHT ADJ COPY OF OCTET 28)
C       (12)  - LONGITUDINAL DIR GRID LENGTH
C       (13)  - LATITUDINAL DIR GRID LENGTH
C    LAMBERT CONFORMAL GRIDS
C       (2) - NX NR POINTS ALONG X-AXIS
C       (3) - NY NR POINTS ALONG Y-AXIS
C       (4) - LA1 LAT OF ORIGIN (LOWER LEFT)
C       (5) - LO1 LON OF ORIGIN (LOWER LEFT)
C       (6) - RESOLUTION (RIGHT ADJ COPY OF OCTET 17)
C       (7) - LOV - ORIENTATION OF GRID
C       (8) - DX - X-DIR INCREMENT
C       (9) - DY - Y-DIR INCREMENT
C       (10)  - PROJECTION CENTER FLAG
C       (11)  - SCANNING MODE FLAG (RIGHT ADJ COPY OF OCTET 28)
C       (12)  - LATIN 1 - FIRST LAT FROM POLE OF SECANT CONE INTER
C       (13)  - LATIN 2 - SECOND LAT FROM POLE OF SECANT CONE INTER
C    LB          LOGICAL*1 (KF) BITMAP IF PRESENT
C    F          REAL (KF) DATA
C OUTPUT ARGUMENTS:
C    IRET       INTEGER RETURN CODE
C                   0    ALL OK
C                   OTHER  W3FI72 GRIB PACKER RETURN CODE
C
C SUBPROGRAMS CALLED:
C R63W72       MAP W3FI63 PARAMETERS ONTO W3FI72 PARAMETERS
C GETBIT       GET NUMBER OF BITS AND ROUND DATA
C W3FI72       PACK GRIB
C WRYTE       WRITE DATA
C
C REMARKS: SUBPROGRAM CAN BE CALLED FROM A MULTIPROCESSING ENVIRONMENT.
C DO NOT ENGAGE THE SAME LOGICAL UNIT FROM MORE THAN ONE PROCESSOR.
C
C ATTRIBUTES:
C LANGUAGE: FORTRAN 77
C MACHINE:  CRAY, WORKSTATIONS
C
C$$$
   INTEGER KPDS(200),KGDS(200)
   LOGICAL*1 LB(KF)
   REAL F(KF)
   PARAMETER(MAXBIT=32)
   INTEGER IBM(KF),IPDS(200),IGDS(200),IBDS(200)
   REAL FR(KF)
   CHARACTER PDS(400),GRIB(1000+KF*(MAXBIT+1)/8)
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  GET W3FI72 PARAMETERS
   CALL R63W72(KPDS,KGDS,IPDS,IGDS)
   IBDS=0
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  COUNT VALID DATA
   KBM=KF
   IF(IPDS(7).NE.0) THEN
      KBM=0
      DO I=1,KF
      IF(LB(I)) THEN
         IBM(I)=1
         KBM=KBM+1
      ELSE
         IBM(I)=0
      ENDIF
      ENDDO
      IF(KBM.EQ.KF) IPDS(7)=0
   ENDIF
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  GET NUMBER OF BITS AND ROUND DATA
   IF(KBM.EQ.0) THEN
      DO I=1,KF
      FR(I)=0.
      ENDDO
      NBIT=0
   ELSE
      CALL GETBIT(IPDS(7),0,IPDS(25),KF,IBM,F,FR,FMIN,FMAX,NBIT)
      NBIT=MIN(NBIT,MAXBIT)
   ENDIF
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  PACK AND WRITE GRIB DATA
   CALL W3FI72(0,FR,0,NBIT,0,IPDS,PDS,
   &          1,255,IGDS,0,0,IBM,KF,IBDS,
   &          KFO,GRIB,LGRIB,IRET)
   IF(IRET.EQ.0) CALL WRYTE(LUGB,LGRIB,GRIB)
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
   RETURN
   END

gaochao123 · 发表于 2020-5-3 19:37:05

本帖最后由 gaochao123 于 2020-5-3 19:42 编辑

这个是w3fi72.f：
   SUBROUTINE W3FI72(ITYPE,FLD,IFLD,IBITL,
   &                IPFLAG,ID,PDS,
   &                IGFLAG,IGRID,IGDS,ICOMP,
   &                IBFLAG,IBMAP,IBLEN,IBDSFL,
   &                NPTS,KBUF,ITOT,JERR)
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C             .    . .                                     .
C SUBPROGRAM:  W3FI72       MAKE A COMPLETE GRIB MESSAGE
C PRGMMR: FARLEY          ORG: NMC421    DATE:94-11-22
C
C ABSTRACT: MAKES A COMPLETE GRIB MESSAGE FROM A USER SUPPLIED
C ARRAY OF FLOATING POINT OR INTEGER DATA.  THE USER HAS THE
C OPTION OF SUPPLYING THE PDS OR AN INTEGER ARRAY THAT WILL BE
C USED TO CREATE A PDS (WITH W3FI68).  THE USER MUST ALSO
C SUPPLY OTHER NECESSARY INFO; SEE USAGE SECTION BELOW.
C
C PROGRAM HISTORY LOG:
C 91-05-08  R.E.JONES
C 92-07-01  M. FARLEY ADDED GDS AND BMS LOGIC.  PLACED EXISTING
C                         LOGIC FOR BDS IN A ROUTINE.
C 92-10-02  R.E.JONES ADD ERROR EXIT FOR W3FI73
C 93-04-30  R.E.JONES REPLACE DO LOOPS TO MOVE CHARACTER DATA
C                         WITH XMOVEX, USE XSTORE TO ZERO CHARACTER
C                         ARRAY. MAKE CHANGE SO FLAT FIELD WILL PACK.
C 93-08-06  CAVANAUGH MODIFIED CALL TO W3FI75
C 93-10-26  CAVANAUGH ADDED CODE TO RESTORE INPUT FIELD TO ORIGINAL
C                         VALUES IF D-SCALE NOT 0
C 94-01-27  CAVANAUGH ADDED IGDS ARRAY IN CALL TO W3FI75 TO PROVIDE
C                         INFORMATION FOR BOUSTROPHEDONIC PROCESSING
C 94-03-03  CAVANAUGH INCREASED SIZE OF GDS ARRAY FOR THIN GRIDS
C 94-05-16  FARLEY    CLEANED UP DOCUMENTATION
C 94-11-10  FARLEY    INCREASED SIZE OF PFLD/IFLD ARRARYS FROM
C                         100K TO 260K FOR .5 DEGREE SST ANAL FIELDS
C 94-12-04  R.E.JONES CHANGE DOCUMENT FOR IPFLAG.
C 95-10-31  IREDELL    REMOVED SAVES AND PRINTS
C 98-05-19  Gilbert    Increased array dimensions to handle grids
C                         of up to 500,000 grid points.
C 95-10-31  IREDELL    GENERALIZED WORD SIZE
C 98-12-21  Gilbert    Replaced Function ICHAR with mova2i.
C 99-02-01  Gilbert    Changed the method of zeroing out array KBUF.
C                         the old method, using W3FI01 and XSTORE was
C                         incorrect with 4-byte integers and 8-byte reals.
C 2001-06-07  Gilbert    Removed calls to xmovex.
C                         changed IPFLD from integer to character.
C 10-02-19  GAYNO       FIX ALLOCATION OF ARRAY BMS
C
C USAGE:  CALL W3FI72(ITYPE,FLD,IFLD,IBITL,
C       &          IPFLAG,ID,PDS,
C       &          IGFLAG,IGRID,IGDS,ICOMP,
C       &          IBFLAG,IBMAP,IBLEN,IBDSFL,
C       &          IBDSFL,
C       &          NPTS,KBUF,ITOT,JERR)
C
C INPUT ARGUMENT LIST:
C    ITYPE - 0 = FLOATING POINT DATA SUPPLIED IN ARRAY 'FLD'
C             1 = INTEGER DATA SUPPLIED IN ARRAY 'IFLD'
C    FLD    - REAL ARRAY OF DATA (AT PROPER GRIDPOINTS) TO BE
C             CONVERTED TO GRIB FORMAT IF ITYPE=0.
C             SEE REMARKS #1 & 2.
C    IFLD    - INTEGER ARRAY OF DATA (AT PROPER GRIDPOINTS) TO BE
C             CONVERTED TO GRIB FORMAT IF ITYPE=1.
C             SEE REMARKS #1 & 2.
C    IBITL - 0 = COMPUTER COMPUTES LENGTH FOR PACKING DATA FROM
C                   POWER OF 2 (NUMBER OF BITS) BEST FIT OF DATA
C                   USING 'VARIABLE' BIT PACKER W3FI58.
C             8, 12, ETC. COMPUTER RESCALES DATA TO FIT INTO THAT
C                   'FIXED' NUMBER OF BITS USING W3FI59.
C             SEE REMARKS #3.
C
C    IPFLAG - 0 = MAKE PDS FROM USER SUPPLIED ARRAY (ID)
C             1 = USER SUPPLYING PDS
C             NOTE: IF PDS IS GREATER THAN 30, USE IPLFAG=1.
C             THE USER COULD CALL W3FI68 BEFORE HE CALLS
C             W3FI72. THIS WOULD MAKE THE FIRST 30 BYTES OF
C             THE PDS, USER THEN WOULD MAKE BYTES AFTER 30.
C    ID    - INTEGER ARRAY OF  VALUES THAT W3FI68 WILL USE
C             TO MAKE AN EDITION 1 PDS IF IPFLAG=0.  (SEE THE
C             DOCBLOCK FOR W3FI68 FOR LAYOUT OF ARRAY)
C    PDS    - CHARACTER ARRAY OF VALUES (VALID PDS SUPPLIED
C             BY USER) IF IPFLAG=1. LENGTH MAY EXCEED 28 BYTES
C             (CONTENTS OF BYTES BEYOND 28 ARE PASSED
C             THROUGH UNCHANGED).
C
C    IGFLAG - 0 = MAKE GDS BASED ON 'IGRID' VALUE.
C             1 = MAKE GDS FROM USER SUPPLIED INFO IN 'IGDS'
C                   AND 'IGRID' VALUE.
C             SEE REMARKS #4.
C    IGRID - # = GRID IDENTIFICATION (TABLE B)
C             255 = IF USER DEFINED GRID; IGDS MUST BE SUPPLIED
C                   AND IGFLAG MUST =1.
C    IGDS    - INTEGER ARRAY CONTAINING USER GDS INFO (SAME
C             FORMAT AS SUPPLIED BY W3FI71 - SEE DOCKBLOCK FOR
C             LAYOUT) IF IGFLAG=1.
C    ICOMP - RESOLUTION AND COMPONENT FLAG FOR BIT 5 OF GDS(17)
C             0 = EARTH ORIENTED WINDS
C             1 = GRID ORIENTED WINDS
C
C    IBFLAG - 0 = MAKE BIT MAP FROM USER SUPPLIED DATA
C             # = BIT MAP PREDEFINED BY CENTER
C             SEE REMARKS #5.
C    IBMAP - INTEGER ARRAY CONTAINING BIT MAP
C    IBLEN - LENGTH OF BIT MAP WILL BE USED TO VERIFY LENGTH
C             OF FIELD (ERROR IF IT DOESN'T MATCH).
C
C    IBDSFL - INTEGER ARRAY CONTAINING TABLE 11 FLAG INFO
C             BDS OCTET 4:
C             (1) 0 = GRID POINT DATA
C                   1 = SPHERICAL HARMONIC COEFFICIENTS
C             (2) 0 = SIMPLE PACKING
C                   1 = SECOND ORDER PACKING
C             (3) ... SAME VALUE AS 'ITYPE'
C                   0 = ORIGINAL DATA WERE FLOATING POINT VALUES
C                   1 = ORIGINAL DATA WERE INTEGER VALUES
C             (4) 0 = NO ADDITIONAL FLAGS AT OCTET 14
C                   1 = OCTET 14 CONTAINS FLAG BITS 5-12
C             (5) 0 = RESERVED - ALWAYS SET TO 0
C       BYTE 6 OPTION 1 NOT AVAILABLE (AS OF 5-16-93)
C             (6) 0 = SINGLE DATUM AT EACH GRID POINT
C                   1 = MATRIX OF VALUES AT EACH GRID POINT
C       BYTE 7 OPTION 0 WITH SECOND ORDER PACKING N/A (AS OF 5-16-93)
C             (7) 0 = NO SECONDARY BIT MAPS
C                   1 = SECONDARY BIT MAPS PRESENT
C             (8) 0 = SECOND ORDER VALUES HAVE CONSTANT WIDTH
C                   1 = SECOND ORDER VALUES HAVE DIFFERENT WIDTHS
C
C OUTPUT ARGUMENT LIST:
C    NPTS    - NUMBER OF GRIDPOINTS IN ARRAY FLD OR IFLD
C    KBUF    - ENTIRE GRIB MESSAGE ('GRIB' TO '7777')
C             EQUIVALENCE TO INTEGER ARRAY TO MAKE SURE IT
C             IS ON WORD BOUNARY.
C    ITOT    - TOTAL LENGTH OF GRIB MESSAGE IN BYTES
C    JERR    - = 0, COMPLETED MAKING GRIB FIELD WITHOUT ERROR
C                1, IPFLAG NOT 0 OR 1
C                2, IGFLAG NOT 0 OR 1
C                3, ERROR CONVERTING IEEE F.P. NUMBER TO IBM370 F.P.
C                4, W3FI71 ERROR/IGRID NOT DEFINED
C                5, W3FK74 ERROR/GRID REPRESENTATION TYPE NOT VALID
C                6, GRID TOO LARGE FOR PACKER DIMENSION ARRAYS
C                   SEE AUTOMATION DIVISION FOR REVISION!
C                7, LENGTH OF BIT MAP NOT EQUAL TO SIZE OF FLD/IFLD
C                8, W3FI73 ERROR, ALL VALUES IN IBMAP ARE ZERO
C
C OUTPUT FILES:
C    FT06F001 - STANDARD FORTRAN OUTPUT PRINT FILE
C
C SUBPROGRAMS CALLED:
C    LIBRARY:
C    W3LIB - W3FI58, W3FI59, W3FI68, W3FI71, W3FI73, W3FI74
C                W3FI75, W3FI76
C    FORTRAN 90 INTRINSIC - BIT_SIZE
C
C REMARKS:
C 1)  IF BIT MAP TO BE INCLUDED IN MESSAGE, NULL DATA SHOULD
C    BE INCLUDED IN FLD OR IFLD.  THIS ROUTINE WILL TAKE CARE
C    OF 'DISCARDING' ANY NULL DATA BASED ON THE BIT MAP.
C 2)  UNITS MUST BE THOSE IN GRIB DOCUMENTATION:  NMC O.N. 388
C    OR WMO PUBLICATION 306.
C 3)  IN EITHER CASE, INPUT NUMBERS WILL BE MULTIPLIED BY
C    '10 TO THE NTH' POWER FOUND IN ID(25) OR PDS(27-28),
C    THE D-SCALING FACTOR, PRIOR TO BINARY PACKING.
C 4)  ALL NMC PRODUCED GRIB FIELDS WILL HAVE A GRID DEFINITION
C    SECTION INCLUDED IN THE GRIB MESSAGE.  ID(6) WILL BE
C    SET TO '1'.
C    - GDS WILL BE BUILT BASED ON GRID NUMBER (IGRID), UNLESS
C       IGFLAG=1 (USER SUPPLYING IGDS).  USER MUST STILL SUPPLY
C       IGRID EVEN IF IGDS PROVIDED.
C 5)  IF BIT MAP USED THEN ID(7) OR PDS(8) MUST INDICATE THE
C    PRESENCE OF A BIT MAP.
C 6)  ARRAY KBUF SHOULD BE EQUIVALENCED TO AN INTEGER VALUE OR
C    ARRAY TO MAKE SURE IT IS ON A WORD BOUNDARY.
C 7)  SUBPROGRAM CAN BE CALLED FROM A MULTIPROCESSING ENVIRONMENT.
C
C ATTRIBUTES:
C LANGUAGE: FORTRAN 90
C
C$$$
C
   REAL          FLD(*)
C
   INTEGER       IBDSFL(*)
   INTEGER       IBMAP(*)
   INTEGER       ID(*)
   INTEGER       IFLD(*)
   INTEGER       IGDS(*)
   INTEGER       IB(4)
   INTEGER       NLEFT, NUMBMS
C
   CHARACTER * 1 BDS11(11)
   CHARACTER * 1 KBUF(*)
   CHARACTER * 1 PDS(*)
   CHARACTER * 1 GDS(200)
   CHARACTER(1),ALLOCATABLE:: BMS(:)
   CHARACTER(1),ALLOCATABLE:: PFLD(:)
   CHARACTER(1),ALLOCATABLE:: IPFLD(:)
   CHARACTER * 1 SEVEN
   CHARACTER * 1 ZERO
C
C
C ASCII REP OF  /'G', 'R', 'I', 'B'/
C
   DATA  IB / 71,  82,  73,  66/
C
   IER = 0
   IBERR  = 0
   JERR = 0
   IGRIBL = 8
   IPDSL  = 0
   LENGDS = 0
   LENBMS = 0
   LENBDS = 0
   ITOSS  = 0
C
C$          1.0 PRODUCT DEFINITION SECTION(PDS).
C
C SET ID(6) TO 1 ...OR... MODIFY PDS(8) ...
C    REGARDLESS OF USER SPECIFICATION...
C NMC GRIB FIELDS WILL ALWAYS HAVE A GDS
C
   IF (IPFLAG .EQ.0) THEN
      ID(6) = 1
      CALL W3FI68(ID,PDS)
   ELSE IF (IPFLAG .EQ. 1) THEN
      IF (IAND(mova2i(PDS(8)),64) .EQ. 64) THEN
C       BOTH GDS AND BMS
      PDS(8) = CHAR(192)
      ELSE IF (mova2i(PDS(8)) .EQ. 0) THEN
C       GDS ONLY
      PDS(8) = CHAR(128)
      END IF
      CONTINUE
   ELSE
C    PRINT *,' W3FI72 ERROR, IPFLAG IS NOT 0 OR 1 IPFLAG = ',IPFLAG
      JERR = 1
      GO TO 900
   END IF
C
C    GET LENGTH OF PDS
C
   IPDSL = mova2i(PDS(1)) * 65536 + mova2i(PDS(2)) * 256 +
   &       mova2i(PDS(3))
C
C$          2.0 GRID DEFINITION SECTION (GDS).
C
C    IF IGFLAG=1 THEN USER IS SUPPLYING THE IGDS INFORMATION
C
   IF (IGFLAG .EQ. 0) THEN
      CALL W3FI71(IGRID,IGDS,IGERR)
      IF (IGERR .EQ. 1) THEN
C       PRINT *,' W3FI71 ERROR, GRID TYPE NOT DEFINED...',IGRID
      JERR = 4
      GO TO 900
      END IF
   END IF
   IF (IGFLAG .EQ. 0  .OR.  IGFLAG .EQ.1) THEN
      CALL W3FI74(IGDS,ICOMP,GDS,LENGDS,NPTS,IGERR)
      IF (IGERR .EQ. 1) THEN
C       PRINT *,' W3FI74 ERROR, GRID REP TYPE NOT VALID...',IGDS(3)
      JERR = 5
      GO TO 900
      ELSE
      END IF
   ELSE
C    PRINT *,' W3FI72 ERROR, IGFLAG IS NOT 0 OR 1 IGFLAG = ',IGFLAG
      JERR = 2
      GO TO 900
   END IF
C
C$          3.0 BIT MAP SECTION (BMS).
C
C    SET ITOSS=1 IF BITMAP BEING USED.  W3FI75 WILL TOSS DATA
C    PRIOR TO PACKING.  LATER CODING WILL BE NEEDED WHEN THE
C    'PREDEFINED' GRIDS ARE FINALLY 'DEFINED'.
C
   IF (mova2i(PDS(8)) .EQ. 64 .OR.
   & mova2i(PDS(8)) .EQ. 192) THEN
      ITOSS = 1
      IF (IBFLAG .EQ. 0) THEN
      IF (IBLEN .NE. NPTS) THEN
C          PRINT *,' W3FI72 ERROR, IBLEN .NE. NPTS = ',IBLEN,NPTS
         JERR = 7
         GO TO 900
      END IF
      IF (MOD(IBLEN,16).NE.0) THEN
         NLEFT  = 16 - MOD(IBLEN,16)
      ELSE
         NLEFT  = 0
      END IF
      NUMBMS = 6 + (IBLEN+NLEFT) / 8
      ALLOCATE(BMS(NUMBMS))
      ZERO = CHAR(00)
      BMS = ZERO
      CALL W3FI73(IBFLAG,IBMAP,IBLEN,BMS,LENBMS,IER)
      IF (IER .NE. 0) THEN
C          PRINT *,' W3FI73 ERROR, IBMAP VALUES ARE ALL ZERO'
         JERR = 8
         GO TO 900
      END IF
      ELSE
C       PRINT *,' BIT MAP PREDEFINED BY CENTER, IBFLAG = ',IBFLAG
      END IF
   END IF
C
C$          4.0 BINARY DATA SECTION (BDS).
C
C$          4.1 SCALE THE DATA WITH D-SCALE FROM PDS(27-28)
C
   JSCALE = mova2i(PDS(27)) * 256 + mova2i(PDS(28))
   IF (IAND(JSCALE,32768).NE.0) THEN
      JSCALE = - IAND(JSCALE,32767)
   END IF
   SCALE  = 10.0 ** JSCALE
   IF (ITYPE .EQ. 0) THEN
      DO 410 I = 1,NPTS
      FLD(I) = FLD(I) * SCALE
  410 CONTINUE
   ELSE
      DO 411 I = 1,NPTS
      IFLD(I) = NINT(FLOAT(IFLD(I)) * SCALE)
  411 CONTINUE
   END IF
C
C$          4.2 CALL W3FI75 TO PACK DATA AND MAKE BDS.
C
   ALLOCATE(PFLD(NPTS*4))
C
   IF(IBDSFL(2).NE.0) THEN
      ALLOCATE(IPFLD(NPTS*4))
      IPFLD=char(0)
   ELSE
      ALLOCATE(IPFLD(1))
   ENDIF
C
   CALL W3FI75(IBITL,ITYPE,ITOSS,FLD,IFLD,IBMAP,IBDSFL,
   &       NPTS,BDS11,IPFLD,PFLD,LEN,LENBDS,IBERR,PDS,IGDS)
C
   IF(IBDSFL(2).NE.0) THEN
C       CALL XMOVEX(PFLD,IPFLD,NPTS*4)
      do ii = 1, NPTS*4
         PFLD(ii) = IPFLD(ii)
      enddo
   ENDIF
      DEALLOCATE(IPFLD)
C
      IF (IBERR .EQ. 1) THEN
      JERR = 3
      GO TO 900
      END IF
C          4.3 IF D-SCALE NOT 0, RESCALE INPUT FIELD TO
C                ORIGINAL VALUE
C
   IF (JSCALE.NE.0) THEN
      DSCALE = 1.0 / SCALE
      IF (ITYPE.EQ.0) THEN
            DO 412 I = 1, NPTS
               FLD(I)  = FLD(I) * DSCALE
  412       CONTINUE
      ELSE
            DO 413 I = 1, NPTS
               FLD(I)  = NINT(FLOAT(IFLD(I)) * DSCALE)
  413       CONTINUE
      END IF
   END IF
C
C$          5.0 OUTPUT SECTION.
C
C$          5.1 ZERO OUT THE OUTPUT ARRAY KBUF.
C
   ZERO = CHAR(00)
   ITOT = IGRIBL + IPDSL + LENGDS + LENBMS + LENBDS + 4
C    PRINT *,'IGRIBL  =',IGRIBL
C    PRINT *,'IPDSL =',IPDSL
C    PRINT *,'LENGDS  =',LENGDS
C    PRINT *,'LENBMS  =',LENBMS
C    PRINT *,'LENBDS  =',LENBDS
C    PRINT *,'ITOT =',ITOT
   KBUF(1:ITOT)=ZERO
C
C$          5.2 MOVE SECTION 0 - 'IS' INTO KBUF (8 BYTES).
C
   ISTART  = 0
   DO 520 I = 1,4
      KBUF(I) = CHAR(IB(I))
  520 CONTINUE
C
   KBUF(5) = CHAR(MOD(ITOT / 65536,256))
   KBUF(6) = CHAR(MOD(ITOT / 256,256))
   KBUF(7) = CHAR(MOD(ITOT       ,256))
   KBUF(8) = CHAR(1)
C
C$          5.3 MOVE SECTION 1 - 'PDS' INTO KBUF (28 BYTES).
C
   ISTART  = ISTART + IGRIBL
   IF (IPDSL.GT.0) THEN
C       CALL XMOVEX(KBUF(ISTART+1),PDS,IPDSL)
      do ii = 1, IPDSL
         KBUF(ISTART+ii) = PDS(ii)
      enddo
   ELSE
C    PRINT *,'LENGTH OF PDS LESS OR EQUAL 0, IPDSL = ',IPDSL
   END IF
C
C$          5.4 MOVE SECTION 2 - 'GDS' INTO KBUF.
C
   ISTART  = ISTART + IPDSL
   IF (LENGDS .GT. 0) THEN
C       CALL XMOVEX(KBUF(ISTART+1),GDS,LENGDS)
      do ii = 1, LENGDS
         KBUF(ISTART+ii) = GDS(ii)
      enddo
   END IF
C
C$          5.5 MOVE SECTION 3 - 'BMS' INTO KBUF.
C
   ISTART  = ISTART + LENGDS
   IF (LENBMS .GT. 0) THEN
C       CALL XMOVEX(KBUF(ISTART+1),BMS,LENBMS)
      do ii = 1, LENBMS
         KBUF(ISTART+ii) = BMS(ii)
      enddo
   END IF
C
C$          5.6 MOVE SECTION 4 - 'BDS' INTO KBUF.
C
C$                MOVE THE FIRST 11 OCTETS OF THE BDS INTO KBUF.
C
   ISTART  = ISTART + LENBMS
C    CALL XMOVEX(KBUF(ISTART+1),BDS11,11)
   do ii = 1, 11
      KBUF(ISTART+ii) = BDS11(ii)
   enddo
C
C$                MOVE THE PACKED DATA INTO THE KBUF
C
   ISTART  = ISTART + 11
   IF (LEN.GT.0) THEN
C       CALL XMOVEX(KBUF(ISTART+1),PFLD,LEN)
      do ii = 1, LEN
         KBUF(ISTART+ii) = PFLD(ii)
      enddo
   END IF
C
C$                ADD '7777' TO END OFF KBUF
C NOTE THAT THESE 4 OCTETS NOT INCLUDED IN ACTUAL SIZE OF BDS.
C
   SEVEN  = CHAR(55)
   ISTART = ITOT - 4
   DO 562 I = 1,4
      KBUF(ISTART+I) = SEVEN
562  CONTINUE
C
900  CONTINUE
   IF(ALLOCATED(BMS)) DEALLOCATE(BMS)
   IF(ALLOCATED(PFLD)) DEALLOCATE(PFLD)
   RETURN
   END

gaochao123 · 发表于 2020-5-3 19:41:03

这个是w3fi74.f
   SUBROUTINE W3FI74 (IGDS,ICOMP,GDS,LENGDS,NPTS,IGERR)
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C             .    . .                                     .
C SUBPROGRAM: W3FI74    CONSTRUCT GRID DEFINITION SECTION (GDS)
C PRGMMR: FARLEY          ORG: W/NMC42 DATE: 93-08-24
C
C ABSTRACT: THIS SUBROUTINE CONSTRUCTS A GRIB GRID DEFINITION
C SECTION.
C
C PROGRAM HISTORY LOG:
C 92-07-07  M. FARLEY ORIGINAL AUTHOR
C 92-10-16  R.E.JONES ADD CODE TO LAT/LON SECTION TO DO
C                      GAUSSIAN GRIDS.
C 93-03-29  R.E.JONES ADD SAVE STATEMENT
C 93-08-24  R.E.JONES CHANGES FOR GRIB GRIDS 37-44
C 93-09-29  R.E.JONES CHANGES FOR GAUSSIAN GRID FOR DOCUMENT
C                      CHANGE IN W3FI71.
C 94-02-15  R.E.JONES CHANGES FOR ETA MODEL GRIDS 90-93
C 95-04-20  R.E.JONES CHANGE 200 AND 201 TO 201 AND 202
C 95-10-31  IREDELL    REMOVED SAVES AND PRINTS
C 98-08-20  BALDWIN    ADD TYPE 203
C 07-03-20  VUONG    ADD TYPE 204
C 10-01-21  GAYNO    ADD GRID 205 - ROTATED LAT/LON A,B,C,D STAGGERS
C
C
C USAGE: CALL W3FI74 (IGDS, ICOMP, GDS, LENGDS, NPTS, IGERR)
C INPUT ARGUMENT LIST:
C    IGDS       - INTEGER ARRAY SUPPLIED BY W3FI71
C    ICOMP    - TABLE 7- RESOLUTION & COMPONENT FLAG (BIT 5)
C                FOR GDS(17) WIND COMPONENTS
C
C OUTPUT ARGUMENT LIST:
C    GDS    - COMPLETED GRIB GRID DEFINITION SECTION
C    LENGDS - LENGTH OF GDS
C    NPTS    - NUMBER OF POINTS IN GRID
C    IGERR    - 1, GRID REPRESENTATION TYPE NOT VALID
C
C REMARKS: SUBPROGRAM CAN BE CALLED FROM A MULTIPROCESSING ENVIRONMENT.
C
C ATTRIBUTES:
C LANGUAGE: CRAY CFT77 FORTRAN 77, IBM370 VS FORTRAN
C MACHINE:  CRAY C916-128, CRAY Y-MP8/864, CRAY Y-MP EL2/256, HDS
C
C$$$
C
   INTEGER    IGDS  (*)
C
   CHARACTER*1 GDS (*)
C
   ISUM  = 0
   IGERR = 0
C
C    PRINT *,' '
C    PRINT *,'(W3FI74-IGDS = )'
C    PRINT *,(IGDS(I),I=1,18)
C    PRINT *,' '
C
C    COMPUTE LENGTH OF GDS IN OCTETS (OCTETS 1-3)
C    LENGDS =  32 FOR LAT/LON, GNOMIC, GAUSIAN LAT/LON,
C                   POLAR STEREOGRAPHIC, SPHERICAL HARMONICS,
C                   ROTATED LAT/LON E-STAGGER
C    LENGDS =  34 ROTATED LAT/LON A,B,C,D STAGGERS
C    LENGDS =  42 FOR MERCATOR, LAMBERT, TANGENT CONE
C    LENGDS = 178 FOR MERCATOR, LAMBERT, TANGENT CONE
C
   IF (IGDS(3) .EQ. 0  .OR.  IGDS(3) .EQ. 2  .OR.
   & IGDS(3) .EQ. 4  .OR.  IGDS(3) .EQ. 5  .OR.
   & IGDS(3) .EQ. 50 .OR.  IGDS(3) .EQ. 201.OR.
   & IGDS(3) .EQ. 202.OR.  IGDS(3) .EQ. 203.OR.
   & IGDS(3) .EQ. 204 ) THEN
      LENGDS = 32
C
C    CORRECTION FOR GRIDS 37-44
C
      IF (IGDS(3).EQ.0.AND.IGDS(1).EQ.0.AND.IGDS(2).NE.
   &  255) THEN
      LENGDS = IGDS(5) * 2 + 32
      ENDIF
   ELSE IF (IGDS(3) .EQ. 1  .OR.  IGDS(3) .EQ. 3  .OR.
   &       IGDS(3) .EQ. 13) THEN
      LENGDS = 42
   ELSE IF (IGDS(3) .EQ. 205) THEN
      LENGDS = 34
   ELSE
C    PRINT *,' W3FI74 ERROR, GRID REPRESENTATION TYPE NOT VALID'
      IGERR = 1
      RETURN
   ENDIF
C
C    PUT LENGTH OF GDS SECTION IN BYTES 1,2,3
C
   GDS(1) = CHAR(MOD(LENGDS/65536,256))
   GDS(2) = CHAR(MOD(LENGDS/  256,256))
   GDS(3) = CHAR(MOD(LENGDS    ,256))
C
C    OCTET 4 = NV, NUMBER OF VERTICAL COORDINATE PARAMETERS
C    OCTET 5 = PV, PL OR 255
C    OCTET 6 = DATA REPRESENTATION TYPE (TABLE 6)
C
   GDS(4) = CHAR(IGDS(1))
   GDS(5) = CHAR(IGDS(2))
   GDS(6) = CHAR(IGDS(3))
C
C    FILL OCTET THE REST OF THE GDS BASED ON DATA REPRESENTATION
C    TYPE (TABLE 6)
C
C$$
C    PROCESS ROTATED LAT/LON A,B,C,D STAGGERS
C
   IF (IGDS(3).EQ.205) THEN
      GDS( 7) = CHAR(MOD(IGDS(4)/256,256))
      GDS( 8) = CHAR(MOD(IGDS(4) ,256))
      GDS( 9) = CHAR(MOD(IGDS(5)/256,256))
      GDS(10) = CHAR(MOD(IGDS(5) ,256))
      LATO = IGDS(6) ! LAT OF FIRST POINT
      IF (LATO .LT. 0) THEN
      LATO = -LATO
      LATO = IOR(LATO,8388608)
      ENDIF
      GDS(11) = CHAR(MOD(LATO/65536,256))
      GDS(12) = CHAR(MOD(LATO/  256,256))
      GDS(13) = CHAR(MOD(LATO    ,256))
      LONO = IGDS(7) ! LON OF FIRST POINT
      IF (LONO .LT. 0) THEN
      LONO = -LONO
      LONO = IOR(LONO,8388608)
      ENDIF
      GDS(14) = CHAR(MOD(LONO/65536,256))
      GDS(15) = CHAR(MOD(LONO/  256,256))
      GDS(16) = CHAR(MOD(LONO    ,256))
      LATEXT  = IGDS(9)  ! CENTER LAT
      IF (LATEXT .LT. 0) THEN
      LATEXT = -LATEXT
      LATEXT = IOR(LATEXT,8388608)
      ENDIF
      GDS(18) = CHAR(MOD(LATEXT/65536,256))
      GDS(19) = CHAR(MOD(LATEXT/  256,256))
      GDS(20) = CHAR(MOD(LATEXT    ,256))
      LONEXT  = IGDS(10) ! CENTER LON
      IF (LONEXT .LT. 0) THEN
      LONEXT = -LONEXT
      LONEXT = IOR(LONEXT,8388608)
      ENDIF
      GDS(21) = CHAR(MOD(LONEXT/65536,256))
      GDS(22) = CHAR(MOD(LONEXT/  256,256))
      GDS(23) = CHAR(MOD(LONEXT    ,256))
      GDS(24) = CHAR(MOD(IGDS(11)/256,256))
      GDS(25) = CHAR(MOD(IGDS(11) ,256))
      GDS(26) = CHAR(MOD(IGDS(12)/256,256))
      GDS(27) = CHAR(MOD(IGDS(12) ,256))
      GDS(28) = CHAR(IGDS(13))
      LATO = IGDS(14) ! LAT OF LAST POINT
      IF (LATO .LT. 0) THEN
      LATO = -LATO
      LATO = IOR(LATO,8388608)
      ENDIF
      GDS(29) = CHAR(MOD(LATO/65536,256))
      GDS(30) = CHAR(MOD(LATO/  256,256))
      GDS(31) = CHAR(MOD(LATO    ,256))
      LONO = IGDS(15) ! LON OF LAST POINT
      IF (LONO .LT. 0) THEN
      LONO = -LONO
      LONO = IOR(LONO,8388608)
      ENDIF
      GDS(32) = CHAR(MOD(LONO/65536,256))
      GDS(33) = CHAR(MOD(LONO/  256,256))
      GDS(34) = CHAR(MOD(LONO    ,256))
C
C    PROCESS LAT/LON GRID TYPES OR GAUSSIAN GRID OR ARAKAWA
C    STAGGERED, SEMI-STAGGERED, OR FILLED E-GRIDS
C
   ELSEIF (IGDS(3).EQ.0.OR.IGDS(3).EQ.4.OR.
   & IGDS(3).EQ.201.OR.IGDS(3).EQ.202.OR.
   & IGDS(3).EQ.203.OR.IGDS(3).EQ.204) THEN
      GDS( 7) = CHAR(MOD(IGDS(4)/256,256))
      GDS( 8) = CHAR(MOD(IGDS(4) ,256))
      GDS( 9) = CHAR(MOD(IGDS(5)/256,256))
      GDS(10) = CHAR(MOD(IGDS(5) ,256))
      LATO = IGDS(6)
      IF (LATO .LT. 0) THEN
      LATO = -LATO
      LATO = IOR(LATO,8388608)
      ENDIF
      GDS(11) = CHAR(MOD(LATO/65536,256))
      GDS(12) = CHAR(MOD(LATO/  256,256))
      GDS(13) = CHAR(MOD(LATO    ,256))
      LONO = IGDS(7)
      IF (LONO .LT. 0) THEN
      LONO = -LONO
      LONO = IOR(LONO,8388608)
      ENDIF
      GDS(14) = CHAR(MOD(LONO/65536,256))
      GDS(15) = CHAR(MOD(LONO/  256,256))
      GDS(16) = CHAR(MOD(LONO    ,256))
      LATEXT  = IGDS(9)
      IF (LATEXT .LT. 0) THEN
      LATEXT = -LATEXT
      LATEXT = IOR(LATEXT,8388608)
      ENDIF
      GDS(18) = CHAR(MOD(LATEXT/65536,256))
      GDS(19) = CHAR(MOD(LATEXT/  256,256))
      GDS(20) = CHAR(MOD(LATEXT    ,256))
      LONEXT  = IGDS(10)
      IF (LONEXT .LT. 0) THEN
      LONEXT = -LONEXT
      LONEXT = IOR(LONEXT,8388608)
      ENDIF
      GDS(21) = CHAR(MOD(LONEXT/65536,256))
      GDS(22) = CHAR(MOD(LONEXT/  256,256))
      GDS(23) = CHAR(MOD(LONEXT    ,256))
      IRES = IAND(IGDS(8),128)
      IF (IGDS(3).EQ.201.OR.IGDS(3).EQ.202.OR.
   &    IGDS(3).EQ.203.OR.IGDS(3).EQ.204) THEN
      GDS(24) = CHAR(MOD(IGDS(11)/256,256))
      GDS(25) = CHAR(MOD(IGDS(11) ,256))
      ELSE IF (IRES.EQ.0) THEN
      GDS(24) = CHAR(255)
      GDS(25) = CHAR(255)
      ELSE
      GDS(24) = CHAR(MOD(IGDS(12)/256,256))
      GDS(25) = CHAR(MOD(IGDS(12) ,256))
      END IF
      IF (IGDS(3).EQ.4) THEN
      GDS(26) = CHAR(MOD(IGDS(11)/256,256))
      GDS(27) = CHAR(MOD(IGDS(11) ,256))
      ELSE IF (IGDS(3).EQ.201.OR.IGDS(3).EQ.202.OR.
   &          IGDS(3).EQ.203.OR.IGDS(3).EQ.204)THEN
      GDS(26) = CHAR(MOD(IGDS(12)/256,256))
      GDS(27) = CHAR(MOD(IGDS(12) ,256))
      ELSE IF (IRES.EQ.0) THEN
      GDS(26) = CHAR(255)
      GDS(27) = CHAR(255)
      ELSE
      GDS(26) = CHAR(MOD(IGDS(11)/256,256))
      GDS(27) = CHAR(MOD(IGDS(11) ,256))
      END IF
      GDS(28) = CHAR(IGDS(13))
      GDS(29) = CHAR(0)
      GDS(30) = CHAR(0)
      GDS(31) = CHAR(0)
      GDS(32) = CHAR(0)
      IF (LENGDS.GT.32) THEN
      ISUM = 0
      I = 19
      DO 10 J = 33,LENGDS,2
         ISUM    = ISUM + IGDS(I)
         GDS(J) = CHAR(MOD(IGDS(I)/256,256))
         GDS(J+1) = CHAR(MOD(IGDS(I) ,256))
         I       = I + 1
10    CONTINUE
      END IF
C
C$$    PROCESS MERCATOR GRID TYPES
C
   ELSE IF (IGDS(3) .EQ. 1) THEN
      GDS( 7) = CHAR(MOD(IGDS(4)/256,256))
      GDS( 8) = CHAR(MOD(IGDS(4) ,256))
      GDS( 9) = CHAR(MOD(IGDS(5)/256,256))
      GDS(10) = CHAR(MOD(IGDS(5) ,256))
      LATO = IGDS(6)
      IF (LATO .LT. 0) THEN
      LATO = -LATO
      LATO = IOR(LATO,8388608)
      ENDIF
      GDS(11) = CHAR(MOD(LATO/65536,256))
      GDS(12) = CHAR(MOD(LATO/  256,256))
      GDS(13) = CHAR(MOD(LATO    ,256))
      LONO = IGDS(7)
      IF (LONO .LT. 0) THEN
      LONO = -LONO
      LONO = IOR(LONO,8388608)
      ENDIF
      GDS(14) = CHAR(MOD(LONO/65536,256))
      GDS(15) = CHAR(MOD(LONO/  256,256))
      GDS(16) = CHAR(MOD(LONO    ,256))
      LATEXT = IGDS(9)
      IF (LATEXT .LT. 0) THEN
      LATEXT = -LATEXT
      LATEXT = IOR(LATEXT,8388608)
      ENDIF
      GDS(18) = CHAR(MOD(LATEXT/65536,256))
      GDS(19) = CHAR(MOD(LATEXT/  256,256))
      GDS(20) = CHAR(MOD(LATEXT    ,256))
      LONEXT  = IGDS(10)
      IF (LONEXT .LT. 0) THEN
      LONEXT = -LONEXT
      LONEXT = IOR(LONEXT,8388608)
      ENDIF
      GDS(21) = CHAR(MOD(LONEXT/65536,256))
      GDS(22) = CHAR(MOD(LONEXT/  256,256))
      GDS(23) = CHAR(MOD(LONEXT    ,256))
      GDS(24) = CHAR(MOD(IGDS(13)/65536,256))
      GDS(25) = CHAR(MOD(IGDS(13)/  256,256))
      GDS(26) = CHAR(MOD(IGDS(13)    ,256))
      GDS(27) = CHAR(0)
      GDS(28) = CHAR(IGDS(14))
      GDS(29) = CHAR(MOD(IGDS(12)/65536,256))
      GDS(30) = CHAR(MOD(IGDS(12)/  256,256))
      GDS(31) = CHAR(MOD(IGDS(12)    ,256))
      GDS(32) = CHAR(MOD(IGDS(11)/65536,256))
      GDS(33) = CHAR(MOD(IGDS(11)/  256,256))
      GDS(34) = CHAR(MOD(IGDS(11)    ,256))
      GDS(35) = CHAR(0)
      GDS(36) = CHAR(0)
      GDS(37) = CHAR(0)
      GDS(38) = CHAR(0)
      GDS(39) = CHAR(0)
      GDS(40) = CHAR(0)
      GDS(41) = CHAR(0)
      GDS(42) = CHAR(0)
C$$    PROCESS LAMBERT CONFORMAL GRID TYPES
   ELSE IF (IGDS(3) .EQ. 3) THEN
      GDS( 7) = CHAR(MOD(IGDS(4)/256,256))
      GDS( 8) = CHAR(MOD(IGDS(4) ,256))
      GDS( 9) = CHAR(MOD(IGDS(5)/256,256))
      GDS(10) = CHAR(MOD(IGDS(5) ,256))
      LATO = IGDS(6)
      IF (LATO .LT. 0) THEN
      LATO = -LATO
      LATO = IOR(LATO,8388608)
      ENDIF
      GDS(11) = CHAR(MOD(LATO/65536,256))
      GDS(12) = CHAR(MOD(LATO/  256,256))
      GDS(13) = CHAR(MOD(LATO    ,256))
      LONO = IGDS(7)
      IF (LONO .LT. 0) THEN
      LONO = -LONO
      LONO = IOR(LONO,8388608)
      ENDIF
      GDS(14) = CHAR(MOD(LONO/65536,256))
      GDS(15) = CHAR(MOD(LONO/  256,256))
      GDS(16) = CHAR(MOD(LONO    ,256))
      LONM = IGDS(9)
      IF (LONM .LT. 0) THEN
      LONM = -LONM
      LONM = IOR(LONM,8388608)
      ENDIF
      GDS(18) = CHAR(MOD(LONM/65536,256))
      GDS(19) = CHAR(MOD(LONM/  256,256))
      GDS(20) = CHAR(MOD(LONM    ,256))
      GDS(21) = CHAR(MOD(IGDS(10)/65536,256))
      GDS(22) = CHAR(MOD(IGDS(10)/  256,256))
      GDS(23) = CHAR(MOD(IGDS(10)    ,256))
      GDS(24) = CHAR(MOD(IGDS(11)/65536,256))
      GDS(25) = CHAR(MOD(IGDS(11)/  256,256))
      GDS(26) = CHAR(MOD(IGDS(11)    ,256))
      GDS(27) = CHAR(IGDS(12))
      GDS(28) = CHAR(IGDS(13))
      GDS(29) = CHAR(MOD(IGDS(15)/65536,256))
      GDS(30) = CHAR(MOD(IGDS(15)/  256,256))
      GDS(31) = CHAR(MOD(IGDS(15)    ,256))
      GDS(32) = CHAR(MOD(IGDS(16)/65536,256))
      GDS(33) = CHAR(MOD(IGDS(16)/  256,256))
      GDS(34) = CHAR(MOD(IGDS(16)    ,256))
      GDS(35) = CHAR(MOD(IGDS(17)/65536,256))
      GDS(36) = CHAR(MOD(IGDS(17)/  256,256))
      GDS(37) = CHAR(MOD(IGDS(17)    ,256))
      GDS(38) = CHAR(MOD(IGDS(18)/65536,256))
      GDS(39) = CHAR(MOD(IGDS(18)/  256,256))
      GDS(40) = CHAR(MOD(IGDS(18)    ,256))
      GDS(41) = CHAR(0)
      GDS(42) = CHAR(0)
C$$    PROCESS POLAR STEREOGRAPHIC GRID TYPES
   ELSE IF (IGDS(3) .EQ. 5) THEN
      GDS( 7) = CHAR(MOD(IGDS(4)/256,256))
      GDS( 8) = CHAR(MOD(IGDS(4) ,256))
      GDS( 9) = CHAR(MOD(IGDS(5)/256,256))
      GDS(10) = CHAR(MOD(IGDS(5) ,256))
      LATO = IGDS(6)
      IF (LATO .LT. 0) THEN
      LATO = -LATO
      LATO = IOR(LATO,8388608)
      ENDIF
      GDS(11) = CHAR(MOD(LATO/65536,256))
      GDS(12) = CHAR(MOD(LATO/  256,256))
      GDS(13) = CHAR(MOD(LATO    ,256))
      LONO = IGDS(7)
      IF (LONO .LT. 0) THEN
      LONO = -LONO
      LONO = IOR(LONO,8388608)
      ENDIF
      GDS(14) = CHAR(MOD(LONO/65536,256))
      GDS(15) = CHAR(MOD(LONO/  256,256))
      GDS(16) = CHAR(MOD(LONO    ,256))
      LONM = IGDS(9)
      IF (LONM .LT. 0) THEN
      LONM = -LONM
      LONM = IOR(LONM,8388608)
      ENDIF
      GDS(18) = CHAR(MOD(LONM/65536,256))
      GDS(19) = CHAR(MOD(LONM/ 256,256))
      GDS(20) = CHAR(MOD(LONM    ,256))
      GDS(21) = CHAR(MOD(IGDS(10)/65536,256))
      GDS(22) = CHAR(MOD(IGDS(10)/  256,256))
      GDS(23) = CHAR(MOD(IGDS(10)    ,256))
      GDS(24) = CHAR(MOD(IGDS(11)/65536,256))
      GDS(25) = CHAR(MOD(IGDS(11)/  256,256))
      GDS(26) = CHAR(MOD(IGDS(11)    ,256))
      GDS(27) = CHAR(IGDS(12))
      GDS(28) = CHAR(IGDS(13))
      GDS(29) = CHAR(0)
      GDS(30) = CHAR(0)
      GDS(31) = CHAR(0)
      GDS(32) = CHAR(0)
   ENDIF
C    PRINT 10,(GDS(IG),IG=1,32)
C10    FORMAT ('  GDS= ',32(1X,Z2.2))
C
C    COMPUTE NUMBER OF POINTS IN GRID BY MULTIPLYING
C    IGDS(4) AND IGDS(5) ... NEEDED FOR PACKER
C
   IF (IGDS(3).EQ.0.AND.IGDS(1).EQ.0.AND.IGDS(2).NE.
   & 255) THEN
      NPTS = ISUM
   ELSE
      NPTS = IGDS(4) * IGDS(5)
   ENDIF
C
C    'IOR' ICOMP-BIT 5 RESOLUTION & COMPONENT FLAG FOR WINDS
C    WITH IGDS(8) INFO (REST OF RESOLUTION & COMPONENT FLAG DATA)
C
   ICOMP = ISHFT(ICOMP,3)
   GDS(17) = CHAR(IOR(IGDS(8),ICOMP))
C
   RETURN
   END

		自动登录	找回密码
密码			立即注册