PROGRAM SNDPST C$$$ MAIN PROGRAM DOCUMENTATION BLOCK C . . . . C MAIN PROGRAM: ETA_SNDP C PRGMMR: ROGERS ORG: NP22 DATE: 1999-09-24 C C ABSTRACT: THIS ROUTINE POSTS PROFILE DATA AND WRITES C OUTPUT IN BUFR FORMAT. THIS REPLACES CODE THAT USED TO C RUN INSIDE OF CHKOUT IN THE ETA MODEL. C C PROGRAM HISTORY LOG: C 95-07-26 MIKE BALDWIN C 96-05-07 MIKE BALDWIN - USE SMASK TO SET SOIL VARS TO MISSING C 96-11-22 MIKE BALDWIN - ADD OPTION OF DOING MONOLITHIC FILE OR C SPLIT OUT FILES OR BOTH C 97-12-16 MIKE BALDWIN - NEW MULTI-LEVEL PARAMETERS (SUCH C AS SOIL MOISTURE) C 98-07-23 ERIC ROGERS - MADE Y2K COMPLIANT C 98-09-29 MIKE BALDWIN - SET ACC/AVE VARS TO MISSING AT T=0 C 99-04-01 GEOFF MANIKIN - MAJOR CHANGES FEATURING A REMOVAL OF C THE DISTNICTION OF CLASS0 - ALL OUTPUT C IS NOW CLASS1. ALSO, THE FIELDS OF C VISIBILITY AND CLOUD BASE PRESSURE C ARE ADDED. C 99-09-03 JIM TUCCILLO - REDUCED MEMORY REQUIREMENTS BY CHANGING C THE SIZE OF PRODAT AND INTRODUCING LPRO. C ALSO, PRODAT'S STRUCTURE WAS CHANGED TO C PROVIDE STRIDE-1 ACCESS. C NOTE: THIS CODE CAN STILL BE MODIFIED TO C REDUCE MEMORY. THE CHANGES TODAY WILL C NOT AFFECT ITS FUNCTIONALITY BUT WILL C ALLOW IT TO RUN ON A WINTERHAWK NODE C WITHOUT PAGING. THE MEMORY REQUIREMENT C IS NOW AT ABOUT 260 MBs. C 00-03-10 GEOFF MANIKIN - CODE CHANGED TO BE READY FOR ETA EXTENSION C TO 60 HOURS C 00--5-15 ERIC ROGERS - PUT NSOIL AND LM1 IN INCLUDED PARAMETER C FILE C 04-11-18 BRAD FERRIER - Added Cu precip rate; separated cloud water, C GEOFF MANIKIN rain, & ice from CWM array => feed into visibility C C USAGE: C INPUT ARGUMENT LIST: C NONE C C OUTPUT FILES: C NONE C C SUBPROGRAMS CALLED: C UTILITIES: C CALWXT C CALHEL C BFRIZE C LIBRARY: C BUFRLIB C C ATTRIBUTES: C LANGUAGE: FORTRAN C MACHINE : CRAY C-90 C$$$ C C*** C*** 7/26/95 M. BALDWIN C*** C*** SOUNDING POST PROCESSOR C*** PROGRAM TO READ THE PROFILE OUTPUT FILE ON THE CRAY C*** AND PRODUCE DIAGNOSTIC QUANTITIES AND PACK INTO BUFR C*** C-------------------------------------------------------------------- C C PARMS FOR HOURLY PROFILER OUTPUT C LM - MAX NUMBER OF VERTICAL LEVELS C NPNT - MAX NUMBER OF OUTPUT TIMES C C TO HOLD ALL VARIABLES FOR ANY CLASS OF OUTPUT C (MAX NO MULTI-LAYER VARIABLES*LM + NO OF SINGLE LAYER VARS) C LCL1ML1 - NUMBER OF MULTI-LAYER VARIABLES IN CLASS 1 OUTPUT C LCL1SL1 - NUMBER OF SINGLE LAYER VARIABLES IN CLASS 1 OUTPUT C LCL1ML - NUMBER OF MULTI-LAYER VARIABLES IN PROFILM C LCL1SL - NUMBER OF SINGLE LAYER VARIABLES IN PROFILM C LCL1SOIL - MAX NUMBER OF SOIL LAYER VARIABLES FOR CLASS 0 OR 1 C NSTAT - NUMBER OF STATIONS C C NOTE: THESE NUMBERS WILL BE LARGER THAN THE NUMBERS C COMING OUT OF THE MODEL IN THE BINARY FILE SINCE C WE ARE COMPUTING SOME ADDITIONAL VARIABLES TO GO INTO C BUFR IN THIS PROGRAM. C C-------------------------------------------------------------------- INCLUDE "parm.nmm" cZhou PARAMETER (LM=60,NPNT=85,NSTAT=1400 PARAMETER (LM=LM1,NPNT=NSTP &, SPVAL=-99999.0,SMISS=1.E10 &, LCL1ML=15,LCL1SL=52,LCL1SOIL=2 &, LCL1ML1=14 &, NWORD=(LCL1ML)*LM+2*LCL1SL+NSOIL*LCL1SOIL !Binbin: NWORD is the # of data in output for each record (each station) &, ROG=287.04/9.8) PARAMETER (SNOCON=1.4594E5,RAINCON=1.1787E4) C LOGICAL LVLWSE,SEQFLG(8),NEED CHARACTER*16 SEQNM1(8), SBSET CHARACTER*80 CLIST1(8),FMTO,ASSIGN CHARACTER*8 CISTAT cwas DIMENSION FPACK(NWORD),PRODAT(NSTAT,NPNT,NWORD) DIMENSION FPACK(NWORD) REAL PRODAT(NWORD,NSTAT,NPNT) C C THE PURPOSE OF LPRO IS TO HOLD THE VALUES OF RISTAT UNTIL C THEY ARE COPIED TO FRODAT. THIS ADDITION WILL ALLOW PRODAT C TO BE A REAL(4) ARRAY ( AND SAVE A CONSIDERABLE AMOUNT OF C MEMORY) . PRODAT CAN BE FURTHER REDUCED WITH SOME MORE EFFORT. C JIM TUCCILLO C REAL(8) LPRO(NSTAT,NPNT) REAL(8) RISTAT cwas REAL(8) PRODAT(NSTAT,NPNT,NWORD),RISTAT REAL(8) FRODAT(NWORD),WORKK(NWORD) DIMENSION P(LM),T(LM),U(LM),V(LM),Q(LM),PINT(LM+1),ZINT(LM+1) REAL CWTR(LM),IMXR(LM),RAIN(LM) INTEGER IDATE(3),NP1(8),LLMH(NSTAT),NLVL(2) EQUIVALENCE (CISTAT,RISTAT) C-------------------------------------------------------------------- C C SET OUTPUT UNITS FOR CLASS 1 PROFILE FILE. C LCLAS1 - OUTPUT UNIT FOR CLASS 1 BINARY FILE C LTBCL1 - INPUT UNIT FOR CLASS 1 BUFR TABLE FILE C LUNCL1 - OUTPUT UNIT FOR CLASS 1 BUFR FILE C C-------------------------------------------------------------------- I N T E G E R & LCLAS1,LTBCL1,LUNCL1,STDOUT C-------------------------------------------------------------------- L O G I C A L & MONOL,BRKOUT C-------------------------------------------------------------------- NAMELIST /MODTOP/ ETOP NAMELIST /OPTION/ MONOL,BRKOUT D A T A & LCLAS1 / 76 / &,LTBCL1 / 32 / &,LUNCL1 / 78 / &,STDOUT / 6 / &,SEQNM1 /'HEADR','PROFILE','SURF','FLUX', & 'HYDR','D10M','SLYR','XTRA'/ &,SEQFLG /.FALSE.,.TRUE.,.FALSE.,.FALSE.,.FALSE.,.FALSE., & .TRUE.,.FALSE./ &,LVLWSE /.TRUE./ PRODAT=0. FMTO='("ln -s ${DIRD}",I5.5,".",I4.4,3I2.2,'// & '" fort.",I2.2)' C C GET MODEL TOP PRESSURE C PTOP=25.0*100.0 READ(5,MODTOP,ERR=12321) 12321 CONTINUE PTOP=ETOP*100.0 C C READ IN SWITCHES TO CONTROL WHETHER TO DO... C MONOL=.TRUE. DO MONOLITHIC FILE C BRKOUT=.TRUE. DO BREAKOUT FILES C MONOL=.TRUE. BRKOUT=.FALSE. READ(11,OPTION,ERR=12322) 12322 CONTINUE C cZhou IFCSTL=-99 cZhou JHR=0 write(*,*) 'NFCST=',NFCST C---------------------------------------------------------------------- C---READ STATION DATA-------------------------------------------------- C---------------------------------------------------------------------- LUNIT=66 LRECPR=4*(8+9+LCL1ML*LM1+LCL1SL) !Binbin: LRECPR is total bytes of one record (one record for one station) !(8+9+LCL1ML*LM1+LCL1SL) is # of one record for each station from input file write(*,*)'LRECPR=',LRECPR OPEN(UNIT=LUNIT,ACCESS='DIRECT',RECL=LRECPR,IOSTAT=IER) czhou NREC=0 czhou 33 CONTINUE cZhou NREC=NREC+1 DO 4000 JHR= 1,NFCST !Binbin: Forecast time loop write(*,*) ' ===== FCST ==== ', JHR DO 3000 NST = 1, NSTAT !Binbin: Station loop, one station on loop !Note: the order of fcst time and stations are reversed in WRF and in ETA/RSM !In ETA/RSM, reading recode is first one station store all fcst hours, so first read ! station, then read fcst time !In WRF, all stations saved in one fcst hour (i.e. one files), so first read ! fcst hour, then read station !So use DO 3000/4000 loop here NREC=(JHR-1)*NSTAT + NST READ(LUNIT,REC=NREC,ERR=999) IHRST,IDATE,IFCST,ISTAT,CISTAT, & (FPACK(N),N=1,9),(FPACK(N),N=10,FPACK(7)) cZhou IF (IFCST.GT.IFCSTL) THEN c INUMS=1 c JHR=JHR+1 c IFCSTL=IFCST c ELSE c INUMS=INUMS+1 c ENDIF INUMS=NST IYR=IDATE(3) IMON=IDATE(1) IDAY=IDATE(2) RLAT=FPACK(1)*DEG RLON=FPACK(2)*DEG ELEV=FPACK(3) LLMH(INUMS)=NINT(FPACK(4)) LMH=NINT(FPACK(4)) write(*,*) 'STATION=',INUMS,' IHRST,IDATE,IFCST,ISTAT,CISTAT=', + IHRST,IDATE,IFCST,ISTAT,CISTAT,' FPACK(1-3)=',(FPACK(K),K=1,3) c do k=1,9 c write(*,*) k, FPACK(K) c end do DO 25 L=1,LMH C REVERSE ORDER SO THAT P(1) IS THE TOP AND P(LMH) IS THE BOTTOM LV=LMH-L+1 P(LV)=FPACK(L+9) T(LV)=FPACK(L+9+LMH) U(LV)=FPACK(L+9+LMH*2) V(LV)=FPACK(L+9+LMH*3) Q(LV)=FPACK(L+9+LMH*4) C CWTR, RAIN, AND IMXR NOW IN SEPARATE ARRAYS CWTR(LV)=FPACK(L+9+LMH*6) RAIN(LV)=FPACK(L+9+LMH*13) IMXR(LV)=FPACK(L+9+LMH*14) !Binbin note: IMXR is added 25 CONTINUE c do L=1,LMH c write(*,'(I5,f10.2,3f7.2,4f10.5)') L,P(L),T(L),U(L),V(L), c + Q(L),CWTR(L),RAIN(L),IMXR(L) c end do c do J=10+15*LMH,FPACK(7) c write(*,*)J, FPACK(J) c end do C USE SEA MASK TO SET SOIL/SFC VARIABLES TO MISSING VALUES C (IF SEA) C SM =FPACK(LCL1ML*LMH+54) IF (SM.GT.0.5) THEN C SMSTAV FPACK(LCL1ML*LMH+15)=SMISS C SUBSHX FPACK(LCL1ML*LMH+21)=SMISS C SNOPCX FPACK(LCL1ML*LMH+22)=SMISS C ACSNOW, SMSTOT, SNO, ACSNOM, SSROFF, BGROFF, SOILTB DO LKJ=20,26 FPACK(LCL1ML*LMH+LKJ+9)=SMISS ENDDO C SFCEXC, VEGFRC, CMC, SMC(1:4), STC(1:4) DO LKJ=34,44 FPACK(LCL1ML*LMH+LKJ+9)=SMISS ENDDO ENDIF c Print out sfc variables c write(*,*) 'single layer variables' c write(*,*) (FPACK(i),i=LCL1ML*LMH+10,LCL1ML*LMH+LCL1SL+10) C C GET PPT FOR CALWXT C PPT =FPACK(LCL1ML*LMH+16) C COMPUTE PINT,ZINT C PINT(1)=PTOP DO L=1,LMH DP1=P(L)-PINT(L) PINT(L+1)=P(L)+DP1 ENDDO ZINT(LMH+1)=FPACK(3) DO L=LMH,1,-1 TV2=T(L)*(1.0+0.608*Q(L)) ZZ=ROG*TV2*ALOG(PINT(L+1)/PINT(L)) ZINT(L)=ZINT(L+1)+ZZ ENDDO C C CALL PRECIP TYPE SUBROUTINES. C RIME=FPACK(LCL1ML*LMH+61) SR=FPACK(LCL1ML*LMH+58) TSKIN=FPACK(LCL1ML*LMH+12) CALL CALWXT(T,Q,P,PINT,LMH,LM,PPT,IWX1) CALL CALWXT_RAMER(T,Q,P,PINT,LMH,LM,PPT,IWX2) CALL CALWXT_BOURG(T,Q,PINT,LMH,LM,PPT,ZINT,IWX3) CALL CALWXT_REVISED(T,Q,P,PINT,LMH,LM,PPT,IWX4) C WARNING: EXPLICIT ALGORITHM. UNDER 18 NOT ADMITTED C WITHOUT PARENT OR GUARDIAN CALL CALWXT_EXPLICIT(LMH,TSKIN,PPT,SR,RIME,IWX5) CALL CALWXT_DOMINANT(PPT,IWX1,IWX2,IWX3,IWX4,IWX5, * CSNO,CICE,CFZR,CRAI) C COMPUTE HELICITY AND STORM MOTION C CALL CALHEL(U,V,P,ZINT,PINT,LMH,LM,HELI,UST,VST) C C COMPUTE VISIBILITY C FIRST, EXTRACT SEA LEVEL PRESSURE SR=FPACK(LCL1ML*LMH+58) SLP=FPACK(LCL1ML*LMH+10) CPRATE=FPACK(LCL1ML*LMH+60) !--- Convective precip rate !---##### Need grid-scale contributions to QRAIN and QSNO from the profilm file!! QRAIN=RAIN(LMH) QSNO=IMXR(LMH) !-- Nearly all grid-scale ice is snow IF (CPRATE .GT. 0.) THEN RAINRATE=(1-SR)*CPRATE TERM1=(T(LMH)/SLP)**0.4167 TERM2=(T(LMH)/(P(LMH)))**0.5833 TERM3=RAINRATE**0.8333 QRAIN=QRAIN+RAINCON*TERM1*TERM2*TERM3 IF (SR .GT. 0.) THEN SNORATE=SR*CPRATE TERM1=(T(LMH)/SLP)**0.47 TERM2=(T(LMH)/(P(LMH)))**0.53 TERM3=SNORATE**0.94 QSNO=QSNO+SNOCON*TERM1*TERM2*TERM3 ENDIF ENDIF TT=T(LMH) QV=Q(LMH) QCD=CWTR(LMH) ! QICE=IMXR(LMH) !--- Nearly all of grid-scale ice is snow QICE=0. PPP=P(LMH) CALL CALVIS(QV,QCD,QRAIN,QICE,QSNO,TT,PPP,HOVI) C COMPUTE CLOUD BASE PRESSURE C FIRST, EXTRACT THE CONVECTIVE CLOUD BASE HBOT=FPACK(LCL1ML*LMH+59) CLIMIT =1.0E-06 NEED = .TRUE. CDBP = SMISS CBOT = 5000 DO L=LMH,1,-1 C GSM C START AT THE FIRST LAYER ABOVE GROUND, AND FIND THE C FIRST LAYER WITH A VALUE OF CLOUD WATER GREATER THAN C THE SIGNIFICANT LIMIT (VALUE DESIGNATED BY Q. ZHAO). C THIS LAYER WILL BE THE CLOUD BOTTOM UNLESS THE BOTTOM C OF THE CONVECTIVE CLOUD (HBOT) IS FOUND BELOW IN WHICH C CASE HBOT BECOMES THE CLOUD BASE LAYER. IF ((CWTR(L)+IMXR(L)).GT.CLIMIT.AND.NEED) THEN CBOT=L IF (HBOT.GT.CBOT) THEN CBOT = HBOT ENDIF NEED=.FALSE. ENDIF ENDDO IF (CBOT.GT.LMH) THEN CDBP=SMISS ELSE CDBP=P(INT(CBOT)) ENDIF C C C SET ACC/AVERAGED VARIABLES TO MISSING IF IFCST=0 C IF (IFCST.EQ.0) THEN DO L=1,LMH FPACK(L+9+LMH*7)=SMISS FPACK(L+9+LMH*8)=SMISS ENDDO DO JK=16,29 FPACK(LCL1ML*LMH+JK)=SMISS ENDDO DO JK=32,34 FPACK(LCL1ML*LMH+JK)=SMISS ENDDO ENDIF C C ADD 9 SINGLE LEVEL VARIABLES TO THE OUTPUT C TACK THEM ON TO THE END; WE DON'T NEED CONVECTIVE C CLOUD BASE OR RIME, THOUGH, SO WRITE OVER THOSE RECORDS C NLENF = FPACK (7) NLENP = 7 + LCL1ML1*LMH + LCL1SL1 FPACK(NLENF-2) = CSNO FPACK(NLENF-1) = CICE FPACK(NLENF) = CFZR FPACK(NLENF+1) = CRAI FPACK(NLENF+2) = UST FPACK(NLENF+3) = VST FPACK(NLENF+4) = HELI FPACK(NLENF+5) = CDBP FPACK(NLENF+6) = HOVI C C PLACE DATA INTO PRODAT IN PROPER LOCATIONS PRODAT (1,INUMS,JHR) = FLOAT(IFCST) PRODAT (2,INUMS,JHR) = FLOAT(ISTAT) C RISTAT is a REAL(8) variable by virtue of the fact that it C is equivalenced to CISTAT. Everything else stored in PRODAT C is REAL(4). We have made PRODAT REAL(4) but need a REAL(8) C array for storing RISTAT - that is what LPRO is. Farther C down in the code, we will pull values out of LPRO and store C in FRODAT ( a REAL(8) array ). cwas PRODAT (3,INUMS,JHR) = RISTAT LPRO ( INUMS,JHR) = RISTAT PRODAT (4,INUMS,JHR) = FPACK (1) PRODAT (5,INUMS,JHR) = FPACK (2) PRODAT (6,INUMS,JHR) = FPACK (3) PRODAT (7,INUMS,JHR) = 1 C FPACK HAS 9 ENTRIES BEFORE THE 1ST VERTICAL PROFILE, WHILE PRODAT HAS C ONLY 7. THIS ACCOUNTS FOR THE IJ-2 C START BUILDING PRODAT FROM FPACK, BUT STOP WHEN WE GET TO THE RAIN C ARRAY WHICH IS NOT WRITTEN TO PRODAT DO IJ = 10, (LCL1ML1-1)*LMH+9 PRODAT (IJ-2,INUMS,JHR) = FPACK (IJ) ENDDO C CONTINUE BUILDING PRODAT, BUT SKIP OVER RAIN. THIS IS ACCOMPLISHED C BY ADDING LMH (THE SIZE OF THE RAIN ARRAY) TO THE FPACK INDEXING DO IJ = (LCL1ML1-1)*LMH+10,NLENP+2 MN = IJ+LMH PRODAT (IJ-2,INUMS,JHR) = FPACK (MN) ENDDO c write(*,*) 'NREC=',NREC,' done' c GOTO 33 c write(*,*) 'STATION#', NST,' done!' 3000 CONTINUE write(*,*) 'READ FCST ', JHR, 'done !' 4000 CONTINUE 999 CONTINUE write(*,*) 'Write all of data into one big Bufr file ...' C C WRITE OUT INDIVIDUAL FILES FOR EACH STATION C IF (BRKOUT) THEN DO I=1,NSTAT NLVL(1)=LLMH(I) NLVL(2)=NSOIL C DO J=1,NFCST DO IJ = 1, NWORD FRODAT(IJ) = PRODAT (IJ,I,J) ENDDO FRODAT(3) = LPRO(I,J) ISTAT=NINT(FRODAT(2)) C IF (J.EQ.1) THEN C C INITIALIZE BUFR LISTS SO BFRHDR WILL BE CALLED THE FIRST C TIME THROUGH. C WRITE(ASSIGN,FMTO) ISTAT,IYR,IMON,IDAY,IHRST,LCLAS1 CALL SYSTEM(ASSIGN) CLIST1(1)=' ' ENDIF C C CALL BUFR-IZING ROUTINE C NSEQ = 8 SBSET = 'ETACLS1' CALL BFRIZE(LTBCL1,LCLAS1,SBSET,IYR,IMON,IDAY,IHRST 1, SEQNM1,SEQFLG,NSEQ,LVLWSE,FRODAT,NLVL,CLIST1,NP1 2, WORKK,IER) IF(IER.NE.0)WRITE(6,1080)ISTAT,IER,FRODAT(1) 1080 FORMAT(' SOME SORT OF ERROR ',2I8,F9.1) C C ENDDO C C FINISHED, CLOSE UP BUFR FILES C NSEQ = 8 CALL BFRIZE(0,LCLAS1,SBSET,IYR,IMON,IDAY,IHRST 1, SEQNM1,SEQFLG,NSEQ,LVLWSE,FRODAT,NLVL,CLIST1,NP1 2, WORKK,IER) ENDDO ENDIF IF (MONOL) THEN C C WRITE OUT ONE FILE FOR ALL STATIONS C C INITIALIZE BUFR LISTS SO BFRHDR WILL BE CALLED THE FIRST C TIME THROUGH. C CLIST1(1)=' ' DO I=1,NSTAT !loop for all stations NLVL(1)=LLMH(I) NLVL(2)=NSOIL C DO J=1,NFCST !loop for all fcst hours DO IJ = 1, NWORD FRODAT(IJ) = PRODAT (IJ,I,J) ENDDO c write(*,*) 'FRODAT =', (FRODAT (IJ),IJ=1,10) FRODAT(3) = LPRO(I,J) ISTAT=NINT(FRODAT(2)) C C CALL BUFR-IZING ROUTINE C NSEQ = 8 SBSET = 'ETACLS1' CALL BFRIZE(LTBCL1,LUNCL1,SBSET,IYR,IMON,IDAY,IHRST 1, SEQNM1,SEQFLG,NSEQ,LVLWSE,FRODAT,NLVL,CLIST1,NP1 2, WORKK,IER) IF(IER.NE.0)WRITE(6,1080)ISTAT,IER,FRODAT(1) C C ENDDO write(*,*) 'BFRIZE STATION ', I ,' done' ENDDO C C FINISHED, CLOSE UP BUFR FILES C NSEQ = 8 CALL BFRIZE(0,LUNCL1,SBSET,IYR,IMON,IDAY,IHRST 1, SEQNM1,SEQFLG,NSEQ,LVLWSE,FRODAT,NLVL,CLIST1,NP1 2, WORKK,IER) ENDIF C WRITE(STDOUT,*) ' END OF SOUNDING POST ' STOP END