SUBROUTINE CYCAVG3_RAP(KFILDO,KFIL10,JD,IDPARS,NDATE, 1 NGRIDC,ND11,NSLAB,IPACK,IWORK,DATA,ND5, 2 LSTORE,ND9,LITEMS,CORE,ND10,NBLOCK,NFETCH, 3 IS0,IS1,IS2,IS4,ND7,ND2X3, 4 ISTAV,L3264B,MISTOT,IER) C C APRIL 2018 SHAFER NEW ROUTINE C FEBRUARY 2020 SHAFER MODIFIED TO ALLOW CYCLE OFFSETS C OF 3 AND 6 HOURS C C PURPOSE C TO COMPUTE A CYCLE-AVERAGED GRID FROM UP TO 3 AVAILABLE C MODEL CYCLES. THE GRIDS USED IN THE CALCULATION HAVE C THE SAME VALID TIME. THIS ROUTINE IS ENTERED WHEN C IDPARS(5)=0 AND IDPARSE(6)>0, AND CAN ACCOMMODATE ANY C GRIDDED FIELD THAT CAN BE FOUND DIRECTLY ON INPUT. C ALLOWABLE VALUES OF IDPARS(6): C C 0001 - 1-H CYCLE OFFSET (USES HH, HH-1, HH-2) C 0003 - 3-H CYCLE OFFSET (USES HH, HH-3, HH-6) C 0006 - 6-H CYCLE OFFSET (USES HH, HH-6, HH-12) C C THE MOST RECENT CYCLE IS WEIGHTED HEAVIEST IN THE C CALCULATION AND OLDEST CYCLE IS WEIGHTED LEAST. ONLY C ONE OF THE THREE NEED BE AVAILABLE, IN WHICH CASE IT C IS GIVEN 100% WEIGHT. C C NOTE: COMPUTATION OF VARIABLES IS NOT SUPPORTED, I.E. C IT IS ASSUMED THAT ANY COMPUTATIONS WILL HAVE ALREADY C BEEN DONE PRIOR TO CYCLE AVERAGING. C C C DATA SET USE C KFIL10 - UNIT NUMBER OF TDL MOS-2000 FILE SYSTEM ACCESS C (INPUT-OUTPUT) C KFILDO - DEFAULT UNIT NUMBER FOR OUTPUT (PRINT) FILE C (OUTPUT) C C VARIABLES C CORE(J) = ARRAY TO STORE OR RETRIEVE DATA IDENTIFIED IN C LSTORE( , ) (J=1,ND10). WHEN CORE( ) IS FULL, C DATA ARE STORED ON DISK. UPON RETURN TO THE C CALLING PROGRAM, THE ARRAY WILL BE IN THE SIZE C OF LX*LY (OUTPUT). C FD1(J),FD2(J),... = WORK ARRAYS (J=1,ND2X3) (INTERNAL) C DATA(J) = DATA ARRAY TO HOLD CYCLE AVERAGED GRID C (J=1,ND5) (OUTPUT) C I,J = LOOP COUNT (INTERNAL) C IDPARS(J) = PARSED, INDIVIDUAL COMPONENTS OF THE PREDICTOR C ID CORRESPONDING TO ID(J) (J=1,15) DEFINED IN C THE CALLING PROGRAM (INPUT) C J=1 -- CCC (CLASS OF VARIABLE) C J=2 -- FFF (SUBCLASS OF VARIABLE) C J=3 -- B (BINARY INDICATOR) C J=4 -- DD (DATA SOURCE, MODEL NUMBER) C J=5 -- V (VERTICAL APPLICATION) C J=6 -- LBLBLBLB (BOTTOM OF LAYER, 0 IF ONLY 1 C LAYER) C J=7 -- LTLTLTLT (TOP OF LAYER) C J=8 -- T (TRANSFORMATION) C J=9 -- RR (RUN TIME OFFSET -- PREVIOUS CYCLE. C IT IS ALWAYS A POSITIVE NUMBER AND C COUNTED BACKWARDS IN TIME) C J=10 -- OT (TIME APPLICATION) C J=11 -- OH (TIME PERIOD IN HOURS) C J=12 -- TAU (PROJECTION IN HOURS) C J=13 -- I (INTERPOLATION TYPE) C J=14 -- S (SMOOTHING INDICATOR) C J=15 -- G (GRID INDICATOR) C IER = STATUS RETURN C 0 = GOOD RETURN C 100 = VARIABLE NOT ACCOMMODATED C 103 = GRID CHARACTERISTICS ARE DIFFERENT C SEE GFETCH FOR OTHER VALUES WHEN IER.NE.0 AND C DATA ARE RETURNED AS MISSING (INTERNAL-OUTPUT) C IPACK(J) = WORK ARRAY (J=1,ND2X3) (INTERNAL) C IS0(J) = MOS-2000 GRIB SECTION 0 IDS (J=1,3) (INTERNAL) C IS1(J) = MOS-2000 GRIB SECTION 1 IDS (J=1,22+) (INTERNAL) C IS2(J) = MOS-2000 GRIB SECTION 2 IDS (J=1,12) C IS2(3) AND IS2(4) ARE USED BY THE CALLING C PROGRAM AS GRID DIMENSION (INTERNAL-OUTPUT) C IS4(J) = MOS-2000 GRIB SECTION 4 IDS (J=1,4) (INTERNAL) C ISTAV = 0 DATA RETURNED ARE GRID DATA (OUTPUT) C IWORK(J) = WORK ARRAY (J=1,ND2X3) (INTERNAL) C JD(J) = THE BASIC INTEGER PREDICTOR ID (J=1,4). THIS IS C THE SAME AS ID(J) EXCEPT THAT THE PORTIONS C PERTAINING TO PROCESSING ARE OMITTED: C B = IDPARS(3) C T = IDPARS(8) C I = IDPARS(13) C S = IDPARS(14) C G = IDPARS(15) C JD( ) IS USED TO HELP IDENTIFY THE BASIC MODEL C FIELDS AS READ FROM THE ARCHIVE (INPUT) C KFIL10 = UNIT NUMBER OF TDL MOS-2000 FILE SYSTEM ACCESS C (INPUT) C KFILDO = DEFAULT UNIT NUMBER FOR OUTPUT (PRINT) FILE C (INPUT) C L3264B = INTEGER WORD LENGTH IN BITS DEPENDING ON THE C MACHINE BEING USED (EITHER 32 OR 64) (INPUT) C LD(J) = WORK ARRAY HOLDING THE 4 ID WORDS OF WIND C DIVERGENCE BEING COMPUTED (J=1,4). (INTERNAL) C LITEMS = THE NUMBER OF ITEMS (COLUMNS) IN LSTORE( , ) C THAT HAVE BEEN USED IN THIS RUN (INPUT) C LSTORE(L,J) = THE ARRAY HOLDING INFORMATION ABOUT THE DATA C STORED (L=1,12) (J=1,LITEMS) (INPUT-OUTPUT) C L=1,4-- THE 4 IDS FOR THE DATA C L=5 -- LOCATION OF STORED DATA. WHEN IN CORE, C THIS IS THE LOCATION IN CORE( ) WHERE C THE DATA START. WHEN ON DISK, THIS IS C MINUS THE RECORD NUMBER WHERE THE DATA C START. C L=6 -- THE NUMBER OF 4-BYTE WORDS STORED C L=7 -- 2 FOR DATA PACKED IN TDL GRIB FORMAT, C 1 OTHERWISE C L=8 -- DATE/TIME OF THE DATA IN FORMAT C YYYYMMDDHH C L=9 -- NUMBER OF TIMES DATA HAVE BEEN RETRIEVED C L=10 -- NUMBER OF THE SLAB IN DIR( , ,L) AND IN C NGRIDC( ,L) DEFINING THE CHARACTERISTICS C OF THE GRID C L=11 -- NUMBER OF PREDICTORS IN THE SORTED LIST C IN IS( ,N) (N=1,NPRED) FOR WHICH THIS C VARIABLE IS NEEDED ONLY ONCE FROM C LSTORE( , ). WHEN IT IS NEEDED MORE C THAN ONCE, THE VALUE IS SET TO 7777. C L=12 -- USED INITIALLY IN ESTABLISHING C MSTORE( , ). LATER USED TO DETERMINE C WHETHER TO KEEP THIS VARIABLE C MISTOT = TOTAL NUMBER OF MISSING ITEMS ENCOUNTED IN C UNPACKING GRIDS (INPUT-OUTPUT) C NBLOCK = BLOCK SIZE IN WORDS OF THE MOS-2000 RANDOM DISK C FILE (INPUT) C ND2X3 = DIMENSION OF FD1( ), FD2( ), FD3( ), FD4( ), C FD5( ), FD6( ), FDMS( ), AND FDSINS( ) (INPUT) C ND5 = DIMENSION OF IPACK( ),IWORK( ), AND FDMD( ). C (INPUT) C ND7 = DIMENSION OF IS0( ), IS1( ), IS2( ), AND IS4( ). C NOT ALL LOCATIONS ARE USED. (INPUT) C ND9 = THE SECOND DIMENSION OF LSTORE( , ) (INPUT) C ND10 = DIMENSION OF CORE( ) (INPUT) C ND11 = MAXIMUM NUMBER OF GRID COMBINATIONS THAT CAN BE C DEALT WITH IN THIS RUN. LAST DIMENSION OF C NGRIDC( , ) (INPUT) C NDATE = DATE/TIME FOR WHICH THE PREDICTOR IS NEEDED C (INPUT) C NFETCH = INCREMENTED EACH TIME GFETCH IS ENTERED. IT IS C A RUNNING COUNT FROM THE BEGINNING OF THE MAIN C PROGRAM. THIS COUNT IS MAINTAINED IN CASE THE C USER NEEDS IT FOR DIAGNOSTICS, ETC. (OUTPUT) C NGRIDC(L,M) = HOLDING THE GRID CHARACTERISTICS (L=1,6) FOR C EACH GRID COMBINATION (M=1,NGRID) (INPUT-OUTPUT) C L=1 -- MAP PROJECTION NUMBER (3=LAMBERT, 5=POLAR C STEREOGRAPHIC) C L=2 -- GRID LENGTH IN METERS C L=3 -- LATITUDE AT WHICH GRID LENGTH IS CORRECT C *1000 C L=4 -- GRID ORIENTATION IN DEGREES *1000 C L=5 -- LATITUDE OF LL CORNER IN DEGREES *1000 C L=6 -- LONGITUDE OF LL CORNER IN DEGREES *1000 C NSLAB = THE NUMBER OF THE SLAB IN DIR( , , ) AND IN C NGRIDC( , ) DEFINING THE CHARACTERISTICS OF THE C GRID. IT IS USED TO IDENTIFY THE DATA SOURCE, C I.E., THE MODEL. (OUTPUT) C 1 2 3 4 5 6 7 X C C NONSYSTEM SUBROUTINES CALLED C NONE C IMPLICIT NONE C CHARACTER FLAG(ND5)*11 INTEGER JD(4),IDPARS(15) INTEGER IPACK(ND5),IWORK(ND5) INTEGER IS0(ND7),IS1(ND7),IS2(ND7),IS4(ND7) INTEGER LSTORE(12,ND9) INTEGER NGRIDC(6,ND11) INTEGER LD1(4),LD2(4),LD3(4) INTEGER I,IER,ISTAV,J,KFILDO,KFIL10,L3264B,LITEMS, 1 MISTOT,NBLOCK,ND2X3,ND5,ND7,ND9,ND10,ND11,NDATE, 2 NFETCH,NSLAB INTEGER NWORDS,NPACK,NTIMES,MISSP,MISSS INTEGER NX1,NY1,NX2,NY2,NX3,NY3 INTEGER NDATE1,NDATE2 INTEGER IER1,IER2,IER3,IER4,ICENTER INTEGER MODCYC,LB C REAL DATA(ND5) REAL FD1(ND2X3),FD2(ND2X3),FD3(ND2X3) REAL CORE(ND10) REAL WT3_06(3),WT3_03(3),WT3_01(3),WT3(3) REAL WT2_06(2),WT2_03(2),WT2_01(2),WT2(2) C C CYCLE WEIGHTS C MOST C RECENT -> OLDEST DATA WT3_01/ 0.50, 0.32, 0.18/ ! 3 CYCLES AVAILABLE, IDPARS(6)=1 DATA WT3_03/ 0.45, 0.34, 0.21/ ! 3 CYCLES AVAILABLE, IDPARS(6)=3 DATA WT3_06/ 0.40, 0.36, 0.24/ ! 3 CYCLES AVAILABLE, IDPARS(6)=6 DATA WT2_01/ 0.60, 0.40 / ! 2 CYCLES AVAILABLE, IDPARS(6)=1 DATA WT2_03/ 0.55, 0.45 / ! 2 CYCLES AVAILABLE, IDPARS(6)=3 DATA WT2_06/ 0.53, 0.47 / ! 2 CYCLES AVAILABLE, IDPARS(6)=6 C C INITIALIZATION C IER=0 ISTAV=0 FD1=9999. FD2=9999. FD3=9999. DATA=9999. C C CHECK FOR ALLOWABLE VALUES OF IDPARS(6) C IF(IDPARS(6).EQ.1.OR.IDPARS(6).EQ.3.OR. 1 IDPARS(6).EQ.6) THEN IF(IDPARS(6).EQ.1) WT3(1:3)=WT3_01(1:3) IF(IDPARS(6).EQ.1) WT2(1:2)=WT2_01(1:2) IF(IDPARS(6).EQ.3) WT3(1:3)=WT3_03(1:3) IF(IDPARS(6).EQ.3) WT2(1:2)=WT2_03(1:2) IF(IDPARS(6).EQ.6) WT3(1:3)=WT3_06(1:3) IF(IDPARS(6).EQ.6) WT2(1:2)=WT2_06(1:2) ELSE WRITE(KFILDO,100) IDPARS(6) 100 FORMAT(' ****CYCAVG3: IDPARS(6)=',I4,' NOT ALLOWED...' 1 ' MISSING VALUES WILL BE RETURNED.') RETURN ENDIF C C FETCH GRID FOR CYCLE NDATE AND PROJ IDPARS(12) C AND STORE IT IN FD1. IER=47 RETURNED FROM GFETCH C IS NOT FATAL, IT JUST WON'T BE USED IN THE AVERAGE. C NOTE: IDPARS(6) IN WORD 2 IS SET TO ZERO HERE. C LD1(1)=JD(1) LD1(2)=JD(2)-(IDPARS(6)*10000) LD1(3)=IDPARS(12) LD1(4)=400 C LD1(4)=IDPARS(13)*100+IDPARS(14)*10+IDPARS(15) C CALL GFETCH(KFILDO,KFIL10,LD1,7777,LSTORE,ND9,LITEMS, 1 IS0,IS1,IS2,IS4,ND7,IPACK,IWORK,FD1,ND2X3, 2 NWORDS,NPACK,NDATE,NTIMES,CORE,ND10, 3 NBLOCK,NFETCH,NSLAB,MISSP,MISSS,L3264B,1,IER1) IF(MISSP.NE.0)MISTOT=MISTOT+1 IF(IER1.NE.0.AND.IER1.NE.47) RETURN C NX1=IS2(3) NY1=IS2(4) C C FETCH GRID FOR CYCLE IDPARS(6) HOURS PREVIOUS TO NDATE C AND STORE IT IN FD2. IER=47 RETURNED FROM GFETCH IS C NOT FATAL, IT JUST WON'T BE USED IN THE AVERAGE. C NOTE: IDPARS(6) IN WORD 2 IS SET TO ZERO HERE. C LD2(1)=JD(1) LD2(2)=JD(2)-(IDPARS(6)*10000) LD2(3)=IDPARS(12)+IDPARS(6) LD2(4)=400 C LD2(4)=IDPARS(13)*100+IDPARS(14)*10+IDPARS(15) C CALL UPDAT(NDATE,-IDPARS(6),NDATE1) CALL GFETCH(KFILDO,KFIL10,LD2,7777,LSTORE,ND9,LITEMS, 1 IS0,IS1,IS2,IS4,ND7,IPACK,IWORK,FD2,ND2X3, 2 NWORDS,NPACK,NDATE1,NTIMES,CORE,ND10, 3 NBLOCK,NFETCH,NSLAB,MISSP,MISSS,L3264B,1,IER2) IF(MISSP.NE.0)MISTOT=MISTOT+1 IF(IER2.NE.0.AND.IER2.NE.47) RETURN C NX2=IS2(3) NY2=IS2(4) C C FETCH GRID FOR CYCLE 2*IDPARS(6) HOURS PREVIOUS TO NDATE C AND STORE IT IN FD3. IER=47 RETURNED FROM GFETCH IS C NOT FATAL, IT JUST WON'T BE USED IN THE AVERAGE. C NOTE: IDPARS(6) IN WORD 2 IS SET TO ZERO HERE. C LD3(1)=JD(1) LD3(2)=JD(2)-(IDPARS(6)*10000) LD3(3)=IDPARS(12)+IDPARS(6)+IDPARS(6) LD3(4)=400 C LD3(4)=IDPARS(13)*100+IDPARS(14)*10+IDPARS(15) C CALL UPDAT(NDATE1,-IDPARS(6),NDATE2) CALL GFETCH(KFILDO,KFIL10,LD3,7777,LSTORE,ND9,LITEMS, 1 IS0,IS1,IS2,IS4,ND7,IPACK,IWORK,FD3,ND2X3, 2 NWORDS,NPACK,NDATE2,NTIMES,CORE,ND10, 3 NBLOCK,NFETCH,NSLAB,MISSP,MISSS,L3264B,1,IER3) IF(MISSP.NE.0)MISTOT=MISTOT+1 IF(IER3.NE.0.AND.IER3.NE.47) RETURN C C FOR PROJECTIONS BEYOND 19 HOURS LOOK FOR EXTENDED RUN C IF(IDPARS(6) .EQ.01.AND. ! IF OFFSET=1 1 IDPARS(12).GE.20.AND. ! IF PROJ>=20 1 IER3 .EQ.47) THEN ! IF OLDEST CYC NOT THERE C MODCYC=MOD(NDATE,100) IF(MODCYC.EQ.03.OR.MODCYC.EQ.09.OR. 1 MODCYC.EQ.15.OR.MODCYC.EQ.21) LB=0 IF(MODCYC.EQ.04.OR.MODCYC.EQ.10.OR. 1 MODCYC.EQ.16.OR.MODCYC.EQ.22) LB=1 IF(MODCYC.EQ.05.OR.MODCYC.EQ.11.OR. 1 MODCYC.EQ.17.OR.MODCYC.EQ.23) LB=2 IF(MODCYC.EQ.00.OR.MODCYC.EQ.06.OR. 1 MODCYC.EQ.12.OR.MODCYC.EQ.18) LB=3 IF(MODCYC.EQ.01.OR.MODCYC.EQ.07.OR. 1 MODCYC.EQ.13.OR.MODCYC.EQ.19) LB=4 IF(MODCYC.EQ.02.OR.MODCYC.EQ.08.OR. 1 MODCYC.EQ.14.OR.MODCYC.EQ.20) LB=5 C LD3(3)=IDPARS(12)+LB CALL UPDAT(NDATE,-LB,NDATE2) CALL GFETCH(KFILDO,KFIL10,LD3,7777,LSTORE,ND9,LITEMS, 1 IS0,IS1,IS2,IS4,ND7,IPACK,IWORK,FD3,ND2X3, 2 NWORDS,NPACK,NDATE2,NTIMES,CORE,ND10, 3 NBLOCK,NFETCH,NSLAB,MISSP,MISSS,L3264B,1,IER4) IF(MISSP.NE.0)MISTOT=MISTOT+1 IF(IER4.NE.0.AND.IER4.NE.47) RETURN C ENDIF C NX3=IS2(3) NY3=IS2(4) C WRITE(KFILDO,'(1X,I10,1X,3(I9.9,1X),I10.10)') NDATE, & (LD1(J),J=1,4) WRITE(KFILDO,'(1X,I10,1X,3(I9.9,1X),I10.10)') NDATE1, & (LD2(J),J=1,4) WRITE(KFILDO,'(1X,I10,1X,3(I9.9,1X),I10.10)') NDATE2, & (LD3(J),J=1,4) C C CHECK GRID CHARACTERISTICS C IF(NX1.NE.NX2.OR.NX1.NE.NX3.OR.NY1.NE.NY2.OR.NY1.NE.NY3) THEN WRITE(KFILDO,200)(LD1(J),J=1,4),NX1,NY1, 2 (LD2(J),J=1,4),NX2,NY2, 4 (LD3(J),J=1,4),NX3,NY3 200 FORMAT(/' ****DIFFERENT GRID CHARACTERISTICS. CYCLE ', 1 'AVERAGE NOT COMPUTED IN CYCAVG3. VALUES OF ', 2 'X*Y ARE: ',3(/,2X,3I10.9,I4.3,4X,I4,' *',I4)) IER=103 RETURN ENDIF C C----------------------------------------------------- C COMPUTE CYCLE AVERAGE C----------------------------------------------------- C DO J=1,NX3*NY3 C FLAG(J)=' NO NO NO' C IF(FD1(J).NE.9999..AND.FD2(J).NE.9999..AND. 1 FD3(J).NE.9999.) THEN ! ALL CYCLES AVAILABLE DATA(J)=(FD1(J)*WT3(1))+(FD2(J)*WT3(2))+ 1 (FD3(J)*WT3(3)) FLAG(J)='YES YES YES' C ELSE IF(FD1(J).NE.9999..AND.FD2(J).NE.9999..AND. 1 FD3(J).EQ.9999.) THEN ! OLDEST CYCLE MISSING DATA(J)=(FD1(J)*WT2(1))+(FD2(J)*WT2(2)) FLAG(J)='YES YES NO' C ELSE IF(FD1(J).NE.9999..AND.FD3(J).NE.9999..AND. 1 FD2(J).EQ.9999.) THEN ! MIDDLE CYCLE MISSING DATA(J)=(FD1(J)*WT2(1))+(FD3(J)*WT2(2)) FLAG(J)='YES NO YES' C ELSE IF(FD2(J).NE.9999..AND.FD3(J).NE.9999..AND. 1 FD1(J).EQ.9999.) THEN ! NEWEST CYCLE MISSING DATA(J)=(FD2(J)*WT2(1))+(FD3(J)*WT2(2)) FLAG(J)=' NO YES YES' C ELSE IF(FD1(J).NE.9999.) THEN ! ONLY NEWEST IS AVAIL DATA(J)=FD1(J) FLAG(J)='YES NO NO' C ELSE IF(FD2(J).NE.9999.) THEN ! ONLY MIDDLE IS AVAIL DATA(J)=FD2(J) FLAG(J)=' NO YES NO' C ELSE IF(FD3(J).NE.9999.) THEN ! ONLY OLDEST IS AVAIL DATA(J)=FD3(J) FLAG(J)=' NO NO YES' C ELSE ! ALL MISSING DATA(J)=9999. ENDIF C ENDDO C C PRINT 'FLAG' STRING FOR CENTER GRIDPOINT TO C KFILDO TO LOG WHICH CYCLES WERE USED IN THE C AVERAGE CALCULATION C ICENTER=INT((NX3*NY3)/2) WRITE(KFILDO,210) NDATE,(JD(J),J=1,4),FLAG(ICENTER) 210 FORMAT(' ****CYCAVG3: ',I10,1X,3(I9.9,1X),I10.10,1X,A11) C C NSLAB WILL GET SET TO ZERO IF A GRID CANNOT BE C FOUND BY GFETCH. SET IT TO 1 HERE FOR SAFETY C BEFORE RETURNING. C NSLAB=1 C 900 RETURN END