PROGRAM GFSGRIDS C C DECEMBER 2016 SAMPLATSKY MDL MOS-2000 C APRIL 2017 SAMPLATSKY MODIFIED SOME DOCUMENTATION C APRIL 2018 HUANG REDIFINED PLAIN TO BE 32 CHARACTER TO C BE CONSISTENT WITH OTHER SUBROUTINES. C ALSO INITIALIZED PLAIN AS 32 BLANK C CHARACTERS. C SEPTEMBER 2020 GHIRARDELLI INITIALIZED IOSTAT TO PASS check C all CHECK C DECEMBER 2022 SAMPLATSKY ADDED LOGIC TO PRODUCE GRIDS C NEEDED TO PERFORM MRMS QC AT C HH:30 AND HH:00 VALID TIMES. C C PURPOSE C READ TDLPACK GRID FILE CONTAINING GFS PW AND LI AT 3-HOUR C PROJECTIONS. PROGRAM WILL TIME INTERPOLATE TO A 15 MINUTE C TIME BETWEEN PROJECTIONS, BASED ON THE DATE READ IN FROM C THE NCEPDATE FILE. THIS PROGRAM IS SET UP TO INTERPOLATE C TO 15 MIN BEFORE, AND 15 MIN AFTER THE TOP OF THE HOUR C SPECIFIED BY THE INGEST DATE. THE GRID THAT IS OUTPUT WILL C BE THE SAME AS THE INGESTED GRID. C C NOTE THIS ROUTINE IS SET FOR A VERY SPECIFIC FUNCTION, AND C EVEN MINOR CHANGES TO THE PROCESS WOULD HAVE SIGNIFICANT C IMPACT ON THIS CODE. C C C VARIABLES C CONVERTX = CHARACTER VARIABLE HOLDING THE KEYWORD USED TO OPEN C OUTPUT FILE WITH THE CORRECT ENDIAN. C (CHARACTER*20) (INTERNAL). C C C VARIABLES C NUMVAR = NUMBER OF INPUT AND OUTPUT MRMS VARIABLES. C (PARAMETER) C IPACK(N) = ARRAY CONTAINING PACKED DATA IN TDLPACK GRID C INPUT FILE (N=1,ND1). (INTERNAL) C NID = MAXIMUM NUMBER OF VARIABLE IDS TO BE PROCESSED. C DIMENSION OF ID_OLD( , ), ID_NEW( , ), JTYPE( ). C ID_OLD(J,K) = VARIABLE IDS ON INPUT FILE (J=1,4; K=1,NID). C (INPUT) C ID_NEW(J,K) = VARIABLE IDS ON OUTPUT FILES (J=1,4; K=1,NID). C (INPUT) C NBYTES = NUMBER OF BYTES IN RECORD TO FOLLOW. C IS0( ) = MOS-2000 GRIB SECTION 0 (8 BYTES) C 1 - "TDLP" C 2 - RECORD LENGTH IN BYTES C IS1( ) = MOS-2000 GRIB SECTION 1 C 9-12 TDL 4-WORD ID C IS2( ) = MOS-2000 GRIB SECTION 2 C 2 - MAP PROJECTION C 3 - NX NUMBER OF GRID POINTS IN X DIRECTION C 4 - NY NUMBER OF GRID POINTS IN Y DIRECTION C 5 - LATITUDE OF LOWER LEFT CORNER, DEGREE*10000 C 6 - LONGITUDE OF LOWER LEFT CORNER, DEGREE*10000 C 7 - GRID ORIENTATION, DEGREE*10000 C 8 - GRID LENGTH, MILLIMETERS C 9 - LATITUDE AT WHICH MESH APPLIES, C DEGREES*10000; ALSO TANGENCY LATITUDE C IS4( ) = MOS-2000 GRIB SECTION 4 C L3264B = WORD LENGTH OF MACHINE YOU ARE RUNNING ON C 32 OR 64. C L3264W = 64/L3264B C ND7 = DIMENSION OF IS0( ), IS2( ), AND IS4( )(=60). C ND1 = DIMENSION OF IPACK( ) AND DATA( ) - THE C PACKED AND UNPACKED DATA, RESPECTIVELY. C ND1 = DIMENSION OF CCALL( ), THE "STATION CALL C LETTERS". C KREC = COUNTER FOR ALL PACKED DATA RECORDS READ FROM C INPUT FILE. (INTERNAL) C ISTRT = FIRST UNPACKED "DATA" RECORD TO PRINT OUT. C (INPUT) C IFIN = LAST UNPACKED "DATA" RECORD TO PRINT OUT. C (INPUT) C JTYPE(N) = THE TYPE OF OUTPUT FOR THE CERTAIN ID (N=1,NID). C = 1 VECTOR FILE, = 2 GRID FILE, = 3 BOTH C C COMMON BLOCKS: NONE C C C SUBPROGRAMS CALLED: C UNPACK, CONTNG, QCANL, WRVTGD C C23456789112345678921234567893123456789412345678951234567896123456789712 C PARAMETER (ND1=3000000,ND7=60,L3264B=32,L3264W=64/L3264B) C DIMENSION GFSDATA(8,ND1),DATA(ND1) C DIMENSION IS0(ND7),IS1(ND7),IS2(ND7),IS4(ND7),IPACK(ND1), 1 IWORK(ND1),ISTA(ND1),IPLAIN(L3264W,4),ID(4),NBYTES(2) C CHARACTER*1 STA1(8),XTMP CHARACTER*32 PLAIN CHARACTER*8 CCALL(ND1),STA8(1) CHARACTER*60 INFILE,OUTFILE,DATFILE CHARACTER*255 JUNK C EQUIVALENCE (IPLAIN,PLAIN) EQUIVALENCE (STA8,STA1) CINTEL CHARACTER(LEN=20) :: CONVERTX C CONVERTX='BIG_ENDIAN' CINTEL DATA KFILDI/5/,KFILDO/6/,KFILVT/99/ DATA PLAIN/' '/ C NBYWD = L3264B/8 C CALL TIMPR(KFILDO,KFILDO,'START GFSGRIDS ') C C READ CONTROL FILE INFO C READ(KFILDI,100) JUNK 100 FORMAT(A255) READ(KFILDI,105) ID1PW READ(KFILDI,105) ID1LI READ(KFILDI,105) MESH 105 FORMAT(I10) C READ(KFILDI,110) KFILDT,DATFILE READ(KFILDI,110) IJUNK,JUNK READ(KFILDI,110) KFILIN,INFILE READ(KFILDI,110) IJUNK,JUNK READ(KFILDI,110) KFILOUT,OUTFILE 110 FORMAT(I3,4X,A60) C C INITIALIZE IOSTAT IOSTAT=0 C WRITE(KFILDO,120) ID1PW,ID1LI,MESH,KFILDT,DATFILE,KFILIN,INFILE, 1 KFILOUT,OUTFILE 120 FORMAT(/,' 1ST WORD OF PW ID: ',I10.9, 1 /,' 1ST WORD OF LI ID: ',I10.9, 2 /,' NOMINAL MESH OF GRID:',I10,/, 2 /,' UNIT AND NAME OF DATE FILE: ',I2,2X,A60, 3 /,' UNIT AND NAME OF INPUT FILE: ',I2,2X,A60, 4 /,' UNIT AND NAME OF OUTPUT FILE: ',I2,2X,A60) C C OPEN DATE FILE WITH CONVERT= SPECIFIER C OPEN(UNIT=KFILDT,FORM='FORMATTED',STATUS='OLD', 1 IOSTAT=IOS,ERR=150) 150 IF (IOSTAT.NE.0) THEN WRITE(KFILDO,155)KFILDT,IOS 155 FORMAT(/,' ****TROUBLE OPENING FILE ON UNIT NO.',I3, 1 '. IOSTAT =',I5,' STOP 150') CALL W3TAGE('LMP_GFSGRIDS') STOP 150 ENDIF C C READ AND PRINT THE DATE TO BE PROCESSED C CALL GET_NCEPDATE(KFILDT,IYR,IMO,IDA,IHR,KDATE,IER) C IF(IER.NE.0)THEN WRITE(KFILDO,160) 160 FORMAT(/' ****ERROR: CAN NOT READ NCEP DATE FILE - ', 1 'FATAL ERROR IN GFSGRIDS. STOP 160.') CALL W3TAGE('LMP_GFSGRIDS') STOP 160 ENDIF C WRITE(KFILDO,170) KDATE 170 FORMAT(/,'DATE-HOUR TO WRITE TO OUTPUT FILE:',2X,I10) C C OPEN TDLPACK GRID INPUT AND OUTPUT FILES (BIG-ENDIAN). C OPEN(UNIT=KFILIN,FORM='UNFORMATTED', 1 STATUS='OLD',CONVERT=CONVERTX,IOSTAT=IOS,ERR=200) C 200 IF (IOSTAT.NE.0) THEN WRITE(KFILDO,210) KFILIN,IOS 210 FORMAT(/,' **** TROUBLE OPENING INPUT FILE UNIT NO.',I3, 1 '. IOSTAT=',I5,' STOP 200') CALL W3TAGE('LMP_GFSGRIDS') STOP 200 END IF C OPEN(UNIT=KFILOUT,FORM='UNFORMATTED', 1 STATUS='NEW',CONVERT=CONVERTX,IOSTAT=IOS,ERR=250) C 250 IF (IOSTAT.NE.0) THEN WRITE(KFILDO,260) KFILOUT,IOS 260 FORMAT(/,' **** TROUBLE OPENING OUTPUT FILE UNIT NO.',I3, 1 '. IOSTAT=',I5,' STOP 250') CALL W3TAGE('LMP_GFSGRIDS') STOP 250 END IF C C BASED ON IHR AS DEFINED ABOVE, SET GFS DATE. EACH HOUR C IS "OFFSET" FROM THE GFS, WHICH RUNS 4 TIMES PER DAY, AND C USES THE LAMP CONCEPT WHERE THE SHORTEST OFFSET IS 4 HRS, C AND MAX IS 9 HRS. C IF (MOD(IHR,6).EQ.0) IOFFSET=6 ! 00, 06, 12, 18 UTC IF (MOD(IHR,6).EQ.1) IOFFSET=7 ! 01, 07, 13, 19 UTC IF (MOD(IHR,6).EQ.2) IOFFSET=8 ! 02, 08, 14, 20 UTC IF (MOD(IHR,6).EQ.3) IOFFSET=9 ! 03, 09, 15, 21 UTC IF (MOD(IHR,6).EQ.4) IOFFSET=4 ! 04, 10, 16, 22 UTC IF (MOD(IHR,6).EQ.5) IOFFSET=5 ! 05, 11, 17, 23 UTC C CALL UPDAT(KDATE,-IOFFSET,IDATEGFS) C C THE INPUT GFS DATA FILE WILL CONTAIN 8 RECORDS FOR A SINGLE C "IDATEGFS" ... READ IN AND SAVE THE ENTIRE FILE. HOWEVER C THE PROGRAM WILL BE STOPPED IF THE DATE DOES NOT MATCH C IDATEGFS, OR THE FIRST WORD OF THE ID DOES NOT MATCH THE C SPECIFIED ID1PW OR ID1LI. C C THE ORDER THE DATA EXIST IN THE FILE WILL BE BY PROJECTION C FOR ID1PW, FOLLOWED BY PROJECTION FOR ID1LI. C DO K=1,8 C 300 READ (UNIT=KFILIN,END=900,ERR=900) 1 (NBYTES(J),J=1,L3264W), 2 (IPACK(J),J=1,NBYTES(L3264W)/NBYWD) C C IF THE DATE FROM THE FILE DOES NOT MATCH IDATEGFS, THIS IS C A FATAL ERROR. THE WRONG INPUT DATA FILE WAS SUPPLIED. C IF (IPACK(5).NE.IDATEGFS) THEN WRITE(KFILDO,310) IPACK(5),IDATEGFS 310 FORMAT(/,' **** THE DATE READ FROM THE INPUT FILE OF',I12, 1 ' DOES NOT MATCH EXPECTED DATE OF',I12,'.', 2 /,' THIS IS A FATAL ERROR. STOP 310') CALL W3TAGE('LMP_GFSGRIDS') STOP 310 END IF C C IF THE FIRST WORD OF THE ID DOES NOT MATCH ID1PW OR ID1LI, C THE INGEST FILE IS NOT VALID. C IF ((IPACK(6).NE.ID1PW).AND.(IPACK(6).NE.ID1LI)) THEN WRITE(KFILDO,320) IPACK(6),ID1PW,ID1LI 320 FORMAT(/,' **** READ ID(1) OF',I10.9,' WHEN EXPECTING TO', 1 ' FIND',I10.9,' OR',I10.9,'. ', 2 ' FATAL ERROR, STOP 320') CALL W3TAGE('LMP_GFSGRIDS') STOP 320 END IF C C LAST CHECK, MAKE SURE THE PROJECTIONS ARE RIGHT. THEY C SHOULD BE 3,6,9,12. IF THEY AREN'T READ IN PROPERLY, C ASSUME AN INCORRECT INGEST AND STOP. C IF ((IPACK(8).LT.3.OR.IPACK(8).GT.12).OR. 1 (MOD(IPACK(8),3).NE.0)) THEN WRITE(KFILDO,330) IPACK(8),IPROJ 330 FORMAT(/,' **** READ IN PROJECTION OF ',I3,' WHEN 3, 6, 9,', 1 ' OR 12 WAS EXPECTED. FATAL ERROR, STOP 330') CALL W3TAGE('LMP_GFSGRIDS') STOP 330 END IF C C ALL CHECKS PASSED, UNPACK AND SAVE THE RECORD. C CALL UNPACK(KFILDO,IPACK,IWORK,DATA,ND1,IS0,IS1,IS2, 1 IS4,ND7,MISSP,MISSS,2,L3264B,IER) C IF (IER.NE.0) THEN WRITE(KFILDO,340) IER,(IPACK(J),J=5,9) 340 FORMAT(/,' **** IER=',I5,' WHEN TRYING TO UNPACK RECORD', 1 I12,4I10.9,' ... FATAL ERROR, STOP 340') CALL W3TAGE('LMP_GFSGRIDS') STOP 340 END IF C C TRANSFER DATA TO GFSDATA, FOR LATER USE. C DO N=1,ND1 GFSDATA(K,N)=DATA(N) END DO C END DO ! DO K=1,8 C C AT THIS POINT, ALL DATA ARE READ IN AND SAVED. DOUBLE LOOP C BELOW WILL DO ALL REMAINING PROCESSING. C DO IMIN=1,4 DO IVAR=1,2 C C BASED ON IHR, DETERMINE WHICH SAVED RECORDS WILL BE USED C ALONG WITH THE PROPER WEIGHTS TO USE FOR THE TIME C INTERPOLATION. C C THE LOGIC BELOW IS CUMBERSOME, NOT READILY SURE OF BETTER C WAY... C C 00,06,12,18 C 400 IF (MOD(IHR,6).EQ.0) THEN IF (IMIN.LE.2) THEN IF (IVAR.EQ.1) THEN N1=1 ! 3-HR PROJ N2=3 ! 6-HR PROJ ELSE N1=2 ! 3-HR PROJ N2=4 ! 6-HR PROJ END IF ELSE IF (IVAR.EQ.1) THEN N1=3 ! 6-HR PROJ N2=5 ! 9-HR PROJ ELSE N1=4 ! 6-HR PROJ N2=6 ! 9-HR PROJ END IF END IF IF (IMIN.EQ.1) THEN WT1=30.0/180.0 WT2=150.0/180.0 ELSE IF (IMIN.EQ.2) THEN WT1=15.0/180.0 WT2=165.0/180.0 ELSE IF (IMIN.EQ.3) THEN WT1=180.0/180.0 WT2=0.0/180.0 ELSE IF (IMIN.EQ.4) THEN WT1=165.0/180.0 WT2=15.0/180.0 END IF C C 01,07,13,19 C ELSE IF (MOD(IHR,6).EQ.1) THEN IF (IMIN.LE.2) THEN IF (IVAR.EQ.1) THEN N1=3 ! 6-HR PROJ N2=5 ! 9-HR PROJ ELSE N1=4 ! 6-HR PROJ N2=6 ! 9-HR PROJ END IF ELSE IF (IVAR.EQ.1) THEN N1=3 ! 6-HR PROJ N2=5 ! 9-HR PROJ ELSE N1=4 ! 6-HR PROJ N2=6 ! 9-HR PROJ END IF END IF IF (IMIN.EQ.1) THEN WT1=150.0/180.0 WT2=30.0/180.0 ELSE IF (IMIN.EQ.2) THEN WT1=135.0/180.0 WT2=45.0/180.0 ELSE IF (IMIN.EQ.3) THEN WT1=120.0/180.0 WT2=60.0/180.0 ELSE IF (IMIN.EQ.4) THEN WT1=105.0/180.0 WT2=75.0/180.0 END IF C C 02,08,14,20 C ELSE IF (MOD(IHR,6).EQ.2) THEN IF (IMIN.LE.2) THEN IF (IVAR.EQ.1) THEN N1=3 ! 6-HR PROJ N2=5 ! 9-HR PROJ ELSE N1=4 ! 6-HR PROJ N2=6 ! 9-HR PROJ END IF ELSE IF (IVAR.EQ.1) THEN N1=3 ! 6-HR PROJ N2=5 ! 9-HR PROJ ELSE N1=4 ! 6-HR PROJ N2=6 ! 9-HR PROJ END IF END IF IF (IMIN.EQ.1) THEN WT1=90.0/180.0 WT2=90.0/180.0 ELSE IF (IMIN.EQ.2) THEN WT1=75.0/180.0 WT2=105.0/180.0 ELSE IF (IMIN.EQ.3) THEN WT1=60.0/180.0 WT2=120.0/180.0 ELSE IF (IMIN.EQ.4) THEN WT1=45.0/180.0 WT2=135.0/180.0 END IF C C 03,09,15,21 C ELSE IF (MOD(IHR,6).EQ.3) THEN IF (IMIN.LE.2) THEN IF (IVAR.EQ.1) THEN N1=3 ! 6-HR PROJ N2=5 ! 9-HR PROJ ELSE N1=4 ! 6-HR PROJ N2=6 ! 9-HR PROJ END IF ELSE IF (IVAR.EQ.1) THEN N1=5 ! 9-HR PROJ N2=7 ! 12-HR PROJ ELSE N1=6 ! 9-HR PROJ N2=8 ! 12-HR PROJ END IF END IF IF (IMIN.EQ.1) THEN WT1=30.0/180.0 WT2=150.0/180.0 ELSE IF (IMIN.EQ.2) THEN WT1=15.0/180.0 WT2=165.0/180.0 ELSE IF (IMIN.EQ.3) THEN WT1=180.0/180.0 WT2=0.0/180.0 ELSE IF (IMIN.EQ.4) THEN WT1=165.0/180.0 WT2=15.0/180.0 END IF C C 04,10,16,22 C ELSE IF (MOD(IHR,6).EQ.4) THEN IF (IMIN.LE.2) THEN IF (IVAR.EQ.1) THEN N1=1 ! 3-HR PROJ N2=3 ! 6-HR PROJ ELSE N1=2 ! 3-HR PROJ N2=4 ! 6-HR PROJ END IF ELSE IF (IVAR.EQ.1) THEN N1=1 ! 3-HR PROJ N2=3 ! 6-HR PROJ ELSE N1=2 ! 3-HR PROJ N2=4 ! 6-HR PROJ END IF END IF IF (IMIN.EQ.1) THEN WT1=150.0/180.0 WT2=30.0/180.0 ELSE IF (IMIN.EQ.2) THEN WT1=135.0/180.0 WT2=45.0/180.0 ELSE IF (IMIN.EQ.3) THEN WT1=120.0/180.0 WT2=60.0/180.0 ELSE IF (IMIN.EQ.4) THEN WT1=105.0/180.0 WT2=75.0/180.0 END IF C C 05,11,17,23 C ELSE IF (MOD(IHR,6).EQ.5) THEN IF (IMIN.LE.2) THEN IF (IVAR.EQ.1) THEN N1=1 ! 3-HR PROJ N2=3 ! 6-HR PROJ ELSE N1=2 ! 3-HR PROJ N2=4 ! 6-HR PROJ END IF ELSE IF (IVAR.EQ.1) THEN N1=1 ! 3-HR PROJ N2=3 ! 6-HR PROJ ELSE N1=2 ! 3-HR PROJ N2=4 ! 6-HR PROJ END IF END IF IF (IMIN.EQ.1) THEN WT1=90.0/180.0 WT2=90.0/180.0 ELSE IF (IMIN.EQ.2) THEN WT1=75.0/180.0 WT2=105.0/180.0 ELSE IF (IMIN.EQ.3) THEN WT1=60.0/180.0 WT2=120.0/180.0 ELSE IF (IMIN.EQ.4) THEN WT1=45.0/180.0 WT2=135.0/180.0 END IF C END IF ! IF-ELSE CHAIN ON IHR C C WITH THE RIGHT PROJECTIONS AND WEIGHTS DETERMINED, C PERFORM TIME INTERPOLATION. IF EITHER NEEDED GRID C CONTAINS A MISSING VALUE, THEN A MISSING VALUE WILL C BE OUTPUT. C DO N=1,ND1 IF(GFSDATA(N1,N).LT.9998.5.AND.GFSDATA(N2,N).LT.9998.5) THEN DATA(N)=(WT1*GFSDATA(N1,N))+(WT2*GFSDATA(N2,N)) ELSE DATA(N)=9999.0 END IF END DO C C SET UP OUTPUT PARAMETERS. C C NDATE WILL BE THE OUTPUT DATE. IF IMIN=1, THEN NDATE C MUST BE 1HR LESS THAN KDATE. C 500 IF (IMIN.LE.2) THEN CALL UPDAT(KDATE,-1,NDATE) ELSE NDATE=KDATE END IF C C SET ID1 BASED ON IVAR. C IF (IVAR.EQ.1) THEN ID(1)=ID1PW ELSE ID(1)=ID1LI END IF C C SET ID2 BASED ON IMIN. THE LAST 2 DIGITS OF THE 2ND C WORD WILL HOLD THE MINUTES. C IF (IMIN.EQ.1) THEN ID(2)=30 ELSE IF (IMIN.EQ.2) THEN ID(2)=44 ELSE IF (IMIN.EQ.3) THEN ID(2)=00 ELSE IF (IMIN.EQ.4) THEN ID(2)=14 END IF C C ID3,4 WILL BE 0. C ID(3)=0 ID(4)=0 C C SETTING PARAMETERS OF OUTPUT GRID DOMAIN C ALAT=IS2(5)*1.0/10000. ALON=IS2(6)*1.0/10000. ORIENT=IS2(7)*1.0/10000. XMESHL=IS2(8)*1.0/1000. XLAT=IS2(9)*1.0/10000. NPROJ=IS2(2) NX=IS2(3) NY=IS2(4) NSTA=NX*NY ITYPE=2 C C C CALL ROUTINE TO PACK AND WRITE EACH QC'D VARIABLE IN C TDLPACK GRID FORMAT. C CALL WRVTGD(KFILDO,KFILVT,KFILOUT,ID,NDATE,CCALL,DATA, 1 NSTA,ALAT,ALON,ORIENT,MESH,XLAT,NX,NY,ND7, 2 IS1(17),IS1(18),ITYPE,L3264B,L3264W,ISTOP, 3 IPLAIN,PLAIN,IERV,IERG,NPROJ) C END DO END DO C C A SUCCESSFUL COMPLETION C CALL TIMPR(KFILDO,KFILDO,'END GFSGRIDS ') STOP C C ERRORS IN READING INGEST DATA C 900 WRITE(KFILDO,910) KFILIN 910 FORMAT(/,' **** ERRORS WERE ENCOUNTERED TRIED TO READ KFILIN.', 1 ' STOP 300') CALL W3TAGE('LMP_GFSGRIDS') STOP 300 C END