SUBROUTINE NGMSLP2 C$$$ SUBPROGRAM DOCUMENTATION BLOCK C . . . C SUBPROGRAM: NGMSLP2 NMC SEA LEVEL PRESSURE REDUCTION C PRGRMMR: TREADON ORG: W/NP2 DATE: 93-02-02 C C ABSTRACT: C C THIS ROUTINE COMPUTES SEA LEVEL PRESSURE USING THE C HYDROSTATIC EQUATION WITH THE SHUELL CORRECTION. THE C FOLLOWING IS BASED ON DOCUMENTATION IN SUBROUTINE C OUTHYDRO OF THE NGM: C C THE FUNDAMENTAL HYDROSTATIC EQUATION IS C D(HEIGHT) C --------- = TAU = VIRTUAL TEMPERATURE * (RGAS/GRAVITY) C D (Z) C WHERE C Z = MINUS LOG OF PRESSURE (-LN(P)). C C SEA-LEVEL PRESSURE IS COMPUTED FROM THE FORMULA C PRESS(MSL) = PRESS(GROUND) * EXP( F) C WHERE C F = HEIGHT OF GROUND / MEAN TAU C MEAN TAU = ( TAU(GRND) + TAU(SL) ) / 2 C C IN THE NGM TAU(GRND) AND TAU(SL) ARE FIRST SET USING A C 6.5DEG/KM LAPSE RATE FROM SIGMA LEVEL K=1. THIS IS MODIFIED C BY A CORRECTION BASED ON THE CRITICAL TAU OF THE SHUELL C CORRECTION: C TAUCR=(RGASD/GRAVITY) * 290.66 C C 1) WHERE ONLY TAU(SL) EXCEEDS TAUCR, CHANGE TAU(SL) TO TAUCR. C C 2) WHERE BOTH TAU(SL) AND TAU(GRND) EXCEED TAUCR, C CHANGE TAU(SL) TO TAUCR-CONST*(TAU(GRND)-TAUCR )**2 C WHERE CONST = .005 (GRAVITY/RGASD) C C THE AVERAGE OF TAU(SL) AND TAU(GRND) IS THEN USED TOGETHER C WITH THE GROUND HEIGHT AND PRESSURE TO DERIVE THE PRESSURE C AT SEA LEVEL. C C EXTRAPOLATING FROM TEMPERATURES IN THE FIRST ATMOSPHERIC ETA C LAYER LEAD TO EXCESSIVELY NOISY UNDERGROUND TEMPERATURE FIELDS. C AFTER EXPERIMENTATION WE OPTED TO USE A MEAN TEMPERATURE FROM C THE LOWEST DP=30MB ABOVE THE GROUND. FROM THIS MEAN TEMPERATURE C WE EXTRAPOLATE TO GET THE GROUND AND SEA LEVEL TEMPERATURES. C USING A TRUE MASS WEIGHTED LAYER MEAN INSTEAD OF THE SIMPLE C ARITHMETIC MEAN CODED SHOWED LITTLE QUALITATIVE DIFFERENCE. C C HEIGHT OF THE 1000MB SURFACE IS COMPUTED FROM THE MSL PRESSURE C FIELD USING THE FORMULA: C C P(MSL) - P(1000MB) = MEAN DENSITY * GRAVITY * HGT(1000MBS) C C WHERE P(MSL) IS THE SEA LEVEL PRESSURE FIELD WE HAVE JUST C COMPUTED. C C BELOW GROUND VIRTUAL TEMPERATURES ARE COMPUTED ASSUMING A CONSTANT C 6.5DEG/KM LAPSE RATE. TO CONVERT VIRTUAL TEMPERATURE TO A DRY C TEMPERATURE WE USE THE DP=30MB MEAN SPECIFIC HUMIDITY. C C LINE COMMENTED OUT WITH 'CX' ARE LEFT OVER FROM EARLIER C EXPERIMENTATION. THEY REMAIN TO REMIND READERS WHAT C WAS TRIED. C . C C PROGRAM HISTORY LOG: C 93-02-02 RUSS TREADON C 98-06-08 T BLACK - CONVERSION FROM 1-D TO 2-D C 00-01-04 JIM TUCCILLO - MPI VERSION C C USAGE: CALL NGMSLP2 C INPUT ARGUMENT LIST: C NONE C C OUTPUT ARGUMENT LIST: C NONE C C OUTPUT FILES: C NONE C C SUBPROGRAMS CALLED: C UTILITIES: C NONE C LIBRARY: C COMMON - VRBLS C EXTRA C MAPOT C LOOPS C DYNAMD C MASKS C PVRBLS C INDX C C ATTRIBUTES: C LANGUAGE: FORTRAN C MACHINE : CRAY C-90 C$$$ C C C INCLUDE PARAMETERS. SET LOCAL PARAMETERS. INCLUDE "parmeta" INCLUDE "params" PARAMETER (GAMMA=6.5/1000.,ZSL=0.0) PARAMETER (TAUCR=RD*GI*290.66,CONST=0.005*G/RD) PARAMETER (GORD=G/RD,DP=60.E2) PARAMETER (ISMTHT=4,ISMTHQ=2) C C DECLARE VARIABLES REAL TBAR(IM,JM),QBAR(IM,JM),ALPBAR(IM,JM) REAL TSFC(IM,JM),QSFC(IM,JM),EGRID(IM,JM) C C INCLUDE COMMON BLOCKS INCLUDE "VRBLS.comm" INCLUDE "EXTRA.comm" INCLUDE "MAPOT.comm" INCLUDE "LOOPS.comm" INCLUDE "DYNAMD.comm" INCLUDE "MASKS.comm" INCLUDE "PVRBLS.comm" INCLUDE "INDX.comm" INCLUDE "CTLBLK.comm" C C********************************************************************** C START NGMSLP HERE. C C LOOP OVER HORIZONTAL GRID. C !$omp parallel do !$omp& private(alpm,alps,l,llmh,nlev,pm,psfc,ptop,qsum,thsfc,tsum) doout30: DO J=JSTA,JEND doin30: DO I=1,IM LLMH = LMH(I,J) C C LOCATE TOP OF LAYER OVER WHICH TO COMPUTE MEAN FIELDS. PSFC = PD(I,J)+PT PTOP = PSFC-DP QSFC(I,J) = (1.-SM(I,J))*QS(I,J)+SM(I,J)*QZ0(I,J) QSFC(I,J) = AMAX1(H1M12,QSFC(I,J)) THSFC = (1.-SM(I,J))*THS(I,J)+SM(I,J)*THZ0(I,J) TSFC(I,J) = THSFC*(P1000/PSFC)**(-CAPA) C C COMPUTE MEAN FIELDS. NLEV = 1 ALPS = D50*(ALPINT(I,J,LLMH)+ALPINT(I,J,LLMH+1)) TSUM = T(I,J,LLMH) QSUM = Q(I,J,LLMH) TBAR(I,J) = TSUM QBAR(I,J) = QSUM ALPBAR(I,J) = ALPS IF (LLMH.EQ.LM) CYCLE doin30 do10: DO L = LLMH-1,1,-1 ALPM = D50*(ALPINT(I,J,L)+ALPINT(I,J,L+1)) PM = EXP(ALPM) IF (PM.LT.PTOP) EXIT do10 NLEV = NLEV + 1 ALPS = ALPS + ALPM TSUM = TSUM + T(I,J,L) QSUM = QSUM + Q(I,J,L) END DO do10 RNLEV = 1./NLEV ALPBAR(I,J) = ALPS*RNLEV TBAR(I,J) = TSUM*RNLEV QBAR(I,J) = QSUM*RNLEV END DO doin30 END DO doout30 C C SMOOTH LAYER MEAN TEMPERATURE AND SPECIFIC HUMIDITY. C BOUND SPECIFIC HUMIDITY SO THAT IT IS NON-NEGATIVE. C CALL P2FILT(ISMTHT,HBM2,TBAR) CALL P2FILT(ISMTHT,HBM2,TSFC) CALL P2FILT(ISMTHQ,HBM2,QBAR) CALL P2FILT(ISMTHQ,HBM2,QSFC) CALL BOUNDL(QBAR,H1M12,H99999,IM,JM) CALL BOUNDL(QSFC,H1M12,H99999,IM,JM) C C LOOP OVER HORIZONTAL GRID. C !$omp parallel do !$omp& private(alpavg,alpsfc,llmh,pavg,psfc,ptop,qavg,rhoavg,rrhog, !$omp& tau,tauavg,tausfc,tausl,tavg,tvrbar,tvrsfc,tvrsl, !$omp& tvrt,tvrtal,zbar,zl,zsfc) doout50: DO J=JSTA,JEND doin50: DO I=1,IM LLMH = LMH(I,J) ZSFC = FIS(I,J)*GI PSFC = PD(I,J)+PT ALPSFC = ALOG(PSFC) PTOP = PSFC-DP SLP(I,J) = PSFC CX IF (LLMH.EQ.LM) GOTO 50 C C COMPUTE LAYER MEAN TAU (VIRTUAL TEMP*RD/G). TVRT = TBAR(I,J)*(H1+D608*QBAR(I,J)) TAU = TVRT*RD*GI C C COMPUTE HEIGHT OF MEAN ATMOSPHERIC LAYER. CX QSFC = (1.-SM(I,J))*QS(I,J)+SM(I,J)*QZ0(I,J) CX QSFC = AMAX1(H1M12,QSFC) QAVG = D50*(QSFC(I,J)+QBAR(I,J)) CX THSFC = (1.-SM(I,J))*THS(I,J)+SM(I,J)*THZ0(I,J) CX TSFC = THSFC*(P1000/PSFC)**(-CAPA) TAVG = D50*(TSFC(I,J)+TBAR(I,J)) TVRBAR = TAVG*(H1+D608*QAVG) ZBAR = RD*GI*TVRBAR*(ALPSFC-ALPBAR(I,J)) + ZSFC C C COMPUTE TAU AT THE GROUND (Z=ZSFC) AND SEA LEVEL (Z=0) C ASSUMING A CONSTANT LAPSE RATE OF GAMMA=6.5DEG/KM. TVRSFC = TVRT + (ZBAR-ZSFC)*GAMMA TAUSFC = TVRSFC*RD*GI TVRSL = TVRT + (ZBAR- ZSL)*GAMMA TAUSL = TVRSL*RD*GI C C IF NEED BE APPLY SHEULL CORRECTION. IF ((TAUSL.GT.TAUCR).AND.(TAUSFC.LE.TAUCR)) THEN TAUSL=TAUCR ELSEIF ((TAUSL.GT.TAUCR).AND.(TAUSFC.GT.TAUCR)) THEN TAUSL = TAUCR-CONST*(TAUSFC-TAUCR)**2 ENDIF C C COMPUTE MEAN TAU. TAUAVG = D50*(TAUSL+TAUSFC) C C COMPUTE SEA LEVEL PRESSURE. IF (LLMH.LT.LM) SLP(I,J) = PSFC*EXP(ZSFC/TAUAVG) C C COMPUTE 1000MB HEIGHTS. ALPAVG = D50*(ALOG(PSFC)+ALOG(SLP(I,J))) PAVG = EXP(ALPAVG) RHOAVG = PAVG*GI/TAUAVG RRHOG = H1/(RHOAVG*G) Z1000(I,J) = (SLP(I,J)-P1000)*RRHOG C C COMPUTE TEMPERATURES ON ETA LEVELS BELOW GROUND. IF ((ZSFC-ZSL).EQ.0.) CYCLE doin50 CX TLAPSE = GORD*(TAUSFC-TAUSL)/(ZSFC-ZSL) IF (LLMH.EQ.LM) CYCLE doin50 do40: DO L = LLMH+1,LM ZL = GI*DFL(L) TVRTL = TVRT + (ZBAR-ZL)*GAMMA T(I,J,L) = TVRTL/(H1+D608*QBAR(I,J)) END DO do40 C C MOVE TO NEXT HORIZONTAL GRIDPOINT. END DO doin50 END DO doout50 C C WE NOW HAVE THE SHEULL SEA LEVEL PRESSURE FIELD C IN ARRAY SLP. APPLY A FIVE POINT SMOOTHER. C !$omp parallel do DO J=JSTA,JEND DO I=1,IM EGRID(I,J)=SLP(I,J) ENDDO ENDDO C CALL EXCH(SLP) do90: DO KS=1,KSLPD do80: DO J=JSTA_M2,JEND_M2 IEND=IM-1-MOD(J+1,2) DO I=2,IEND IF (FIS(I,J).GT.H1000) GO TO 70 IF (FIS(I+IHW(J),J-1).GT.H1000)GO TO 70 IF (FIS(I+IHE(J),J-1).GT.H1000)GO TO 70 IF (FIS(I+IHW(J),J+1).GT.H1000)GO TO 70 IF (FIS(I+IHE(J),J+1).GT.H1000)GO TO 70 IF (KS.GT.1) CYCLE do80 70 EGRID(I,J)=D125* 1 (SLP(I+IHW(J),J-1)+SLP(I+IHE(J),J-1) 1 +SLP(I+IHW(J),J+1)+SLP(I+IHE(J),J+1) 2 +H4*SLP(I,J)) ENDDO END DO do80 DO J=JSTA,JEND DO I=1,IM SLP(I,J) = EGRID(I,J) ENDDO ENDDO END DO do90 C C END OF ROUTINE. C RETURN END