SUBROUTINE LTGTHERMO(KFILDO,KFIL10,IDPARS,JD,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
C APRIL 2019 SHAFER MDL MOS-2000
C JANUARY 2021 SHAFER,
C SAMPLATSKY ADJUSTED SEVERAL IF STATEMENTS
C FROM A.GT.B TO A.GE.B AS THE EQ
C CONDITION WAS NEVER ACCOUNTED.
C ORIGINAL STATEMENTS WERE LEFT
C AS COMMENTED.
C
C PURPOSE
C SUBROUTINE LTGTHERMO COMPUTES VARIOUS THERMODYNAMIC/STABILITY
C PARAMETERS AT MODEL GRID POINTS FOR PREDICTION OF LIGHTNING.
C A SOUNDING IS CONSTRUCTED ABOVE EACH GRID POINT FROM AVAILABLE
C TEMPERATURE AND DEWPOINT DATA AT THE SURFACE AND AT PRESSURE
C LEVELS. THERE MUST BE AT LEAST FIVE LEVELS AVAILABLE IN ORDER
C TO PERFORM CALCULATIONS.
C
C THE FOLLOWING IDPARS(1) AND IDPARS(2) ARE ACCOMMODATED:
C
C 002 003 - TEMPERATURE AT THE LIFTING CONDENSATION LEVEL
C (DEG C) FOR THE MOST UNSTABLE PARCEL BETWEEN THE
C SURFACE UP TO SIGMA LEVEL IDPARS(7)/100.
C IDPARS(7)=0 WILL DEFAULT TO 0.7 SIGMA.
C
C 002 004 - EQUILIBRIUM LEVEL TEMPERATURE (DEC C) FOR THE
C MOST UNSTABLE PARCEL BETWEEN THE SURFACE UP TO
C SIGMA LEVEL IDPARS(7)/100. IDPARS(7)=0 WILL
C DEFAULT TO 0.7 SIGMA. IF NO EQUILIBRIUM LEVEL
C EXISTS THIS IS SET TO THE LCL TEMPERATURE FOR
C THE MOST UNSTABLE PARCEL.
C
C 007 021 - BEST LIFTED INDEX (DEG C) COMPUTED BY FINDING THE
C MOST UNSTABLE PARCEL BETWEEN THE SURFACE UP TO
C SIGMA LEVEL IDPARS(7)/100. IDPARS(7)=0 WILL DEFAULT
C TO 0.7 SIGMA.
C
C 007 104 - MOST UNSTABLE CONVECTIVE AVAILABLE POTENTIAL ENERGY
C (CAPE - J/KG) COMPUTED FOR THE LAYER IDPARS(6)<T<0C
C WHERE IDPARS(6) IS GIVEN IN KELVIN. UPPER LIMIT FOR
C FINDING MOST UNSTABLE PARCEL IS SIGMA IDPARS(7)/100.
C IDPARS(7)=0 WILL DEFAULT TO 0.7 SIGMA.
C
C IDPARS(6) = TEMPERATURE THRESHOLD (IN KELVIN) THAT
C DEFINES THE TOP OF THE LAYER FOR
C COMPUTING CAPE, E.G.,
C 0253 = -20C LEVEL
C 0000 = USE HIGHEST AVAIL LEVEL
C
C 007 105 - MOST UNSTABLE CONVECTIVE AVAILABLE POTENTIAL ENERGY
C (CAPE - J/KG) IN THE LAYER COLDER THAN A SPECIFIED
C TEMPERATURE THRESHOLD GIVEN BY IDPARS(6) IN KELVIN.
C UPPER LIMIT FOR FINDING THE MOST UNSTABLE PARCEL IS
C SIGMA IDPARS(7)/100. IDPARS(7)=0 WILL DEFAULT TO
C 0.7 SIGMA.
C
C IDPARS(6) = TEMPERATURE THRESHOLD (IN KELVIN) THAT
C DEFINES THE BOTTOM OF THE LAYER FOR
C COMPUTING CAPE, E.G.,
C 0273 = 0C LEVEL
C 0000 = DEFAULTS TO FULL SOUNDING
C
C 007 121 - CONVECTIVE INHIBITION (CIN - J/KG) FOR A PARCEL
C ORIGINATING IDPARS(6) HPA FROM THE SURFACE, IN THE
C LAYER BOUNDED BY IDPARS(7)/100 SIGMA AT THE TOP.
C IDPARS(7)=0 WILL DEFAULT TO 0.7 SIGMA.
C
C
C DATA SET USE
C KFILDO = DEFAULT UNIT NUMBER FOR OUTPUT(PRINT) FILE.
C (OUTPUT)
C KFIL10 = UNIT NUMBER OF TDL MOS-2000 FILE SYSTEM
C ACCESS.(INPUT-OUTPUT)
C
C VARIABLES
C CORE(J) = THE ARRAY TO STORE OR RETRIEVE THE DATA
C IDENTIFIED IN LSTORE(,) (J=1,ND10).
C WHEN CORE() IS FULL DATA ARE STORED ON DISK.
C (INPUT)
C DATA(J) = COMPUTED VARIABLE (J=1,ND5). (OUTPUT)
C I = LOOP CONTROL VARIABLE
C IDPARS(J) = THE PARSED, INDIVIDUAL COMPONENTS OF THE
C PREDICTOR ID CORRESPONDING TO ID() (J=1,15).
C (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
C 1 LAYER)
C J=7--LTLTLTLT (TOP OF LAYER)
C J=8--T (TRANSFORMATION)
C J=9--RR (RUN TIME OFFSET, ALWAYS + AND
C BACK 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 SEE GFETCH FOR OTHER VALUES.
C WHEN IER NE 0, DATA ARE RETURNED AS MISSING.
C (INTERNAL-OUTPUT)
C IPACK(J) = WORK ARRAY (J=1,ND5). (INTERNAL)
C IS0(J) = MOS-2000 GRIB SECTION 0 ID'S (J=1,3).
C (INTERNAL)
C IS1(J) = MOS-2000 GRIB SECTION 1 ID'S (J=1,22+).
C (INTERNAL)
C IS2(J) = MOS-2000 GRIB SECTION 2 ID'S (J=1,12).
C IS2(3) AND IS2(4) ARE USED BY THE CALLING
C PROGRAM AS THE GRID DIMENSIONS.
C (INTERNAL-OUTPUT)
C IS4(J) = MOS-2000 GRIB SECTION 4 ID'S (J=1,4).
C (INTERNAL)
C ISTAV = 0 SINCE THE DATA RETURNED ARE GRID DATA.
C (OUTPUT)
C IWORK(J) = WORK ARRAY (J=1,ND5). (INTERNAL)
C JD(J) = THE BASIC INTEGER PREDICTOR ID (J=1,4).
C THIS IS THE SAME AS ID(J), EXCEPT THAT
C THE PORTIONS PERTAINING TO PROCESSING
C ARE OMITTED:
C B = IDPARS(3),
C T = IDPARS(8),
C I = IDPARS(13),
C S = IDPARS(14),
C G = IDPARS(15), AND
C THRESH.
C ID() IS USED TO HELP IDENTIFY THE BASIC MODEL
C FIELDS AS READ FROM THE ARCHIVE. (INPUT)
C KFILDO = DEFAULT UNIT NUMBER FOR OUTPUT (PRINT) FILE.
C (INPUT)
C KFIL10 = UNIT NUMBER OF TDL MOS-2000 FILE SYSTEM ACCESS.
C (INPUT)
C L3264B = INTEGER WORD LENGTH IN BITS OF MACHINE BEING
C USED (EITHER 32 OR 64). (INPUT)
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
C DATA STORED (L=1,12) (J=1,LITEMS).
C (INPUT-OUTPUT)
C L=1,4--THE 4 ID'S 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,
C THIS IS MINUS THE RECORD NUMBER WHERE
C THE DATA START.
C L=6 --THE NUMBER OF 4-BYTE WORDS STORED.
C L=7 --2 FOR DATA PACKED IN TDL GRIB, 1 FOR NOT.
C L=8 --THE 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
C IN NGRIDC(,L) DEFINING THE
C CHARACTERISTICS OF THIS GRID.
C L=11 --THE NUMBER OF THE PREDICTOR IN THE SORTED
C LIST IN ID(,N) (N=1,NPRED) FOR WHICH
C THIS VARIABLE IS NEEDED, WHEN IT IS
C NEEDED ONLY ONCE FROM LSTORE(,).
C WHEN IT IS NEEDED MORE THAN ONCE, THE
C VALUE IS SET = 7777.
C L=12 --USED INITIALLY IN ESTABLISHING
C MSTORE(,). LATER USED AS A WAY OF
C DETERMINING WHETHER TO KEEP THIS
C VARIABLE.
C MDPARS() = PARSED ID USED IN SUBROUTINE PRSID1 FOR
C SUBROUTINE DEWPOINT
C MISSP = PRIMARY MISSING VALUE INDICATOR. RETURNED AS
C 0 WHEN DATA ARE NOT PACKED. (INTERNAL)
C MISSS = SECONDARY MISSING VALUE INDICATOR. RETURNED AS
C 0 WHEN DATA ARE NOT PACKED. (INTERNAL)
C MISTOT = TOTAL NUMBER OF TIMES A MISSING INDICATOR
C HAS BEEN ENCOUNTERED IN UNPACKING GRIDS.
C (INPUT-OUTPUT)
C NBLOCK = THE BLOCK SIZE IN WORDS OF THE MOS-2000 RANDOM
C DISK FILE. (INPUT)
C ND2X3 = DIMENSION OF SEVERAL VARIABLES. THE SIZE OF
C THE GRID IS NOT KNOWN BEFORE FDTK AND FDDP
C ARE FETCHED. ALL WORK ARRAYS ARE DIMENSIONED
C ND2X3 (INPUT)
C ND5 = DIMENSION OF IPACK( ), IWORK( ), AND
C DATA( ). (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
C BE DEALT WITH ON THIS RUN. LAST DIMENSION
C OF NGRIDC(,). (INPUT)
C NDATE = THE DATE/TIME FOR WHICH PREDICTOR IS NEEDED.
C (INPUT)
C NFETCH = INCREMENTED EACH TIME GFETCH IS ENTERED.
C IT IS A RUNNING COUNT FROM THE BEGINNING OF
C THE PROGRAM. THIS COUNT IS MAINTAINED IN
C CASE THE USER NEEDS IT(DIAGNOSTICS, ETC.).
C (OUTPUT)
C NGRIDC(L,M) = HOLDS THE GRID CHARACTERISTICS (L=1,6) FOR
C EACH GRID COMBINATION (M=1,NGRID).
C L=1--MAP PROJECTION NUMBER (3=LAMBERT, 5=
C POLAR STEREOGRAPHIC).
C L=2--GRID LENGTH IN METERS.
C L=3--LATITUDE AT WHICH THE GRID LENGTH IS
C CORRECT *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
C *1000.
C NPACK = 2 FOR TDL GRIB PACKED DATA; 1 FOR NOT PACKED
C THIS IS RETURNED FROM GFETCH. (INTERNAL)
C NSLAB = THE NUMBER OF THE SLAB IN DIR(, ,) AND
C IN NGRIDC(,) DEFINING THE CHARACTERISTICS
C OF THIS GRID. (OUTPUT)
C NTIMES = THE NUMBER OF TIMES, INCLUDING THIS ONE,
C THAT THE RECORD HAS BEEN FETCHED. THIS IS
C STORED IN LSTORE(9,). (INTERNAL)
C NWORDS = NUMBER OF WORDS RETURNED IN DATA(). THIS
C IS RETURNED FROM GFETCH (INTERNAL)
C 1 2 3 4 5 6 7 X
C
IMPLICIT NONE
C
INTEGER, PARAMETER :: MAXLEV=38 !MAX NUMBER OF PRESSURE LEVELS
C
INTEGER JD(4),IDPARS(15),LD(4),LDPARS(15)
INTEGER IS0(ND7),IS1(ND7),IS2(ND7),IS4(ND7)
INTEGER IPACK(ND5),IWORK(ND5)
INTEGER LSTORE(12,ND9)
INTEGER NGRIDC(6,ND11)
INTEGER ICCCFFF,I,J,K,L,N,Z,N500
INTEGER KMAX,NLVL
INTEGER NX,NY,ISTAV1,IER1,NLVLT,NLVLTD
INTEGER IER,ISTAV,KFILDO,KFIL10,L3264B,LITEMS,
1 MISSP,MISSS,MISTOT,NBLOCK,ND2X3,
2 ND5,ND7,ND9,ND10,ND11,
3 NDATE,NFETCH,NPACK,NSLAB,NTIMES,NWORDS
INTEGER TAVAIL(MAXLEV)
C
REAL CORE(ND10)
REAL DATA(ND5),FD1(ND5),FD2(ND5),FD3(ND5),FD4(ND5),
1 PRESS(ND5,MAXLEV),TK(ND5,MAXLEV),TDK(ND5,MAXLEV)
REAL P(MAXLEV),T(MAXLEV),D(MAXLEV)
REAL PP(MAXLEV),TT(MAXLEV),DD(MAXLEV)
REAL TTC(MAXLEV),DDC(MAXLEV),TP(MAXLEV)
REAL IP,TRPT,LCLP,PB,SATLFT,ALCL
REAL P1,TDIFF,LAPS
REAL PO,TO,DO,PT,THW,TP500,ELP
REAL PPP,TTE,TTP,PTOP,PBOT
C
IER=0
ISTAV=0
DATA=9999.
PRESS(1:ND5,1:MAXLEV)=99999.
TK(1:ND5,1:MAXLEV)=9999.
TDK(1:ND5,1:MAXLEV)=9999.
TAVAIL(1:MAXLEV)=0
C
ICCCFFF=IDPARS(1)*1000+IDPARS(2)
C
C CHECK IF VARIABLE IS ACCOMMODATED
C
IF(ICCCFFF.NE.002003.AND.ICCCFFF.NE.002004.AND.
1 ICCCFFF.NE.007021.AND.ICCCFFF.NE.007104.AND.
2 ICCCFFF.NE.007105.AND.ICCCFFF.NE.007121) THEN
WRITE(KFILDO,101)(JD(J),J=1,4)
101 FORMAT(/,' ****IDPARS(1) AND IDPARS(2) DO NOT INDICATE',
1 ' LTGTHERMO PREDICTOR. ',I9.9,2I10.9,I4.3,' NOT',
2 ' COMPUTED IN LTGTHERMO.')
IER=100
RETURN
ENDIF
C
C FETCH SURFACE PRESSURE AND STORE IT IN PRESS(,1)
C (FOR NAM FETCH 0-30MB BOUNDARY LAYER PRESSURE THEN
C ADD 15 MB TO IT TO APPROXIMATE A SURFACE PRESSURE).
C
IF(IDPARS(4).EQ.07) THEN
LD(1)=001107000+IDPARS(4)
LD(2)=970
ELSE
LD(1)=001100000+IDPARS(4)
LD(2)=0
ENDIF
LD(3)=IDPARS(9)*1000000+IDPARS(12)
LD(4)=0
CALL GFETCH(KFILDO,KFIL10,LD,7777,LSTORE,ND9,LITEMS,
1 IS0,IS1,IS2,IS4,ND7,IPACK,IWORK,PRESS(1:ND5,1),
2 ND2X3,NWORDS,NPACK,NDATE,NTIMES,CORE,ND10,
3 NBLOCK,NFETCH,NSLAB,MISSP,MISSS,L3264B,1,IER1)
IF(MISSP.NE.0)MISTOT=MISTOT+1
IF(IER1.EQ.0) THEN
NX=IS2(3)
NY=IS2(4)
IF(IDPARS(4).EQ.07) THEN
DO J=1,NX*NY
PRESS(J,1)=PRESS(J,1)+1500.
ENDDO
ENDIF
ELSE
WRITE(KFILDO,111)(JD(J),J=1,4)
111 FORMAT(/' ****VARIABLE ',I9.9,2I10.9,I4.3,
1 ' NOT COMPUTED IN LTGTHERMO. SURFACE PRESSURE COULD',
2 ' NOT BE FETCHED.')
IER=101
RETURN
ENDIF
C
C FILL PRESSURE ARRAY IN 25 HPA INCREMENTS
C
DO J=1,NX*NY
PRESS(J,1)=PRESS(J,1)/100.
ENDDO
DO K=2,MAXLEV
PRESS(1:ND5,K)=1000.-((K-2)*25.)
ENDDO
C
C FETCH TEMPERATURE AT EACH PRESSURE LEVEL
C AND STORE IT IN TK(:,K).
C
NLVLT=0
DO K=1,MAXLEV
IF(K.EQ.1) THEN
LD(1)=002001000+IDPARS(4)
LD(2)=2
ELSE
LD(1)=002000000+IDPARS(4)
LD(2)=1000.-((K-2)*25.)
ENDIF
LD(3)=IDPARS(9)*1000000+IDPARS(12)
LD(4)=0
CALL GFETCH(KFILDO,KFIL10,LD,7777,LSTORE,ND9,LITEMS,
1 IS0,IS1,IS2,IS4,ND7,IPACK,IWORK,TK(1:ND5,K),
2 ND2X3,NWORDS,NPACK,NDATE,NTIMES,CORE,ND10,
3 NBLOCK,NFETCH,NSLAB,MISSP,MISSS,L3264B,1,IER1)
IF(MISSP.NE.0)MISTOT=MISTOT+1
IF(IER1.EQ.0) THEN
TAVAIL(K)=1
NLVLT=NLVLT+1
ENDIF
IF(IS2(3).NE.NX.OR.IS2(4).NE.NY) THEN
WRITE(KFILDO,125)
125 FORMAT(/,' ****GRID CHARACTERISTICS OF INPUT GRIDS ARE',
1 ' DIFFERENT IN LTGTHERMO.')
IER=102
RETURN
ENDIF
ENDDO
C
C FETCH OR COMPUTE DEWPOINT AT EACH PRESSURE LEVEL
C AND STORE IT IN TDK(:,K). ONLY CHECK LEVELS
C THAT HAD TEMPERATURE TO SAVE TIME.
C
NLVLTD=0
DO K=1,MAXLEV
IF(TAVAIL(K).EQ.0) CYCLE
IF(K.EQ.1) THEN
LD(1)=003101000+IDPARS(4)
LD(2)=2
ELSE
LD(1)=003100000+IDPARS(4)
LD(2)=1000.-((K-2)*25.)
ENDIF
LD(3)=IDPARS(9)*1000000+IDPARS(12)
LD(4)=0
CALL GFETCH(KFILDO,KFIL10,LD,7777,LSTORE,ND9,LITEMS,
1 IS0,IS1,IS2,IS4,ND7,IPACK,IWORK,TDK(1:ND5,K),
2 ND2X3,NWORDS,NPACK,NDATE,NTIMES,CORE,ND10,
3 NBLOCK,NFETCH,NSLAB,MISSP,MISSS,L3264B,1,IER1)
IF(MISSP.NE.0)MISTOT=MISTOT+1
IF(IER1.EQ.0) THEN
NLVLTD=NLVLTD+1
CYCLE
ENDIF
CALL PRSID1(KFILDO,LD,LDPARS)
CALL DEWPT(KFILDO,KFIL10,LDPARS,LD,NDATE,
1 NGRIDC,ND11,NSLAB,IPACK,IWORK,TDK(1:ND5,K),ND5,
2 LSTORE,ND9,LITEMS,CORE,ND10,NBLOCK,
3 NFETCH,IS0,IS1,IS2,IS4,ND7,
4 FD1,FD2,FD3,FD4,ND2X3,
5 ISTAV1,L3264B,MISTOT,IER1)
IF(IER1.EQ.0)NLVLTD=NLVLTD+1
IF(MISSP.NE.0)MISTOT=MISTOT+1
IF(IS2(3).NE.NX.OR.IS2(4).NE.NY) THEN
WRITE(KFILDO,125)
IER=102
RETURN
ENDIF
ENDDO
C
C CHECK THAT WE HAVE AT LEAST 5 LEVELS AVAILABLE
C
IF(NLVLT.LT.5.OR.NLVLTD.LT.5) THEN
WRITE(KFILDO,150)(JD(J),J=1,4)
150 FORMAT(/' ****VARIABLE ',I9.9,2I10.9,I4.3,
1 ' NOT COMPUTED IN LTGTHERMO. INSUFFICIENT NUMBER',
2 ' OF PRESSURE LEVELS AVAILABLE FOR T AND TD.')
IER=103
RETURN
ENDIF
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C BEGIN MAIN LOOP TO COMPUTE REQUESTED PARAMETER
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
!$OMP PARALLEL DO DEFAULT(SHARED)
!$OMP& PRIVATE(J,K,KMAX,N,P,T,D,NLVL,PP,TT,DD,TTC,DDC,IP,N500,
!$OMP& PT,PO,THW,TO,DO,LCLP,TP,ELP,L,P1,Z,PPP,TTE,TTP,
!$OMP& TP500,PTOP,PBOT,PB)
C
DO J=1,NX*NY
C
IF(PRESS(J,1).EQ.99999..OR.
1 TK(J,1).EQ.9999..OR.TDK(J,1).EQ.9999.) CYCLE
C
C WRITE(13,*) ''
C WRITE(13,*) 'GRIDPOINT= ',J
C DO K=MAXLEV,1,-1
C WRITE(13,'(3(f7.2,2x))') PRESS(J,K),TK(J,K),
C 1 TDK(J,K)
C ENDDO
C
C DETERMINE HIGHEST AVAILABLE PRESSURE LEVEL
C
DO K=1,MAXLEV
IF(TK(J,K).NE.9999..AND.TDK(J,K).NE.9999.) KMAX=K
ENDDO
C
C CONSTRUCT SOUNDING AT THIS GRID POINT FROM
C AVAILABLE PRESSURE LEVELS
C
N=1
P(1)=PRESS(J,1)
T(1)=TK(J,1)
D(1)=TDK(J,1)
DO K=2,KMAX
IF(TK(J,K).NE.9999..AND.TDK(J,K).NE.9999..AND.
2 PRESS(J,K).LT.PRESS(J,1)) THEN
N=N+1
P(N)=PRESS(J,K)
T(N)=TK(J,K)
D(N)=TDK(J,K)
ENDIF
ENDDO
C
C WRITE(14,*) ''
C WRITE(14,*) 'GRIDPOINT= ',J
C DO K=N,1,-1
C WRITE(14,'(3(f7.2,2x))') P(K),T(K),D(K)
C ENDDO
C
C INTERPOLATE TO FILL IN MISSING PRESSURE LEVELS
C EVERY 25 HPA UP TO HIGHEST AVAILABLE LEVEL.
C
NLVL=1
PP(1)=PRESS(J,1)
TT(1)=TK(J,1)
TTC(1)=TK(J,1)-273.15
DD(1)=TDK(J,1)
DDC(1)=TDK(J,1)-273.15
C
DO K=2,KMAX
IP=1000.-((K-2)*25.)
IF(IP.LT.PP(1)) THEN
NLVL=NLVL+1
PP(NLVL)=IP
IF(PP(NLVL).EQ.500.) N500=NLVL ! LEVEL THAT CORRESPONDS TO 500 HPA
IF(TK(J,K).EQ.9999.) THEN
TT(NLVL)=TRPT(P,T,N,IP)
TTC(NLVL)=TT(NLVL)-273.15
ELSE
TT(NLVL)=TK(J,K)
TTC(NLVL)=TT(NLVL)-273.15
ENDIF
IF(TDK(J,K).EQ.9999.) THEN
DD(NLVL)=TRPT(P,D,N,IP)
DDC(NLVL)=DD(NLVL)-273.15
ELSE
DD(NLVL)=TDK(J,K)
DDC(NLVL)=DD(NLVL)-273.15
ENDIF
ENDIF
ENDDO
C
C WRITE(15,*) ''
C WRITE(15,*) 'GRIDPOINT= ',J
C DO K=NLVL,1,-1
C WRITE(15,'(3(f7.2,2x))') PP(K),TTC(K),DDC(K)
C ENDDO
C
C SET UPPER LIMIT FOR FINDING MOST UNSTABLE PARCEL
C (ALSO USED AS UPPER LIMIT FOR COMPUTING CIN)
C
IF(IDPARS(7).EQ.0) THEN
PT=PP(1)*0.7
ELSE
PT=PP(1)*IDPARS(7)/100.
IF(PT.LT.PP(NLVL)) PT=PP(NLVL)
ENDIF
C
C CALL ROUTINE MUPARCEL TO FIND THE PRESSURE OF THE
C MOST UNSTABLE PARCEL BETWEEN THE SFC AND PT. THIS IS
C THE PARCEL WITH THE WARMEST WETBULB POTENTIAL
C TEMPERATURE, THW.
C
CALL MUPARCEL(TTC,DDC,PP,PT,NLVL,PO,THW)
C
C FIND TEMP/DEWP (C) OF MOST UNSTABLE PARCEL BY
C INTERPOLATING FROM SOUNDING AT PRESSURE PO
C
TO=TRPT(P,T,N,PO)-273.15
DO=TRPT(P,D,N,PO)-273.15
C
C FIND LCL PRESSURE (MB) FOR MOST UNSTABLE PARCEL -
C COMPUTED BY FUNCTION ALCL
C
LCLP=ALCL(TO,DO,PO)
IF(LCLP.LT.PP(NLVL)) LCLP=PP(NLVL)
C
SELECT CASE(ICCCFFF)
C
CASE(002003) ! LCL TEMPERATURE
C
C FIND LCL TEMPERATURE (DEG C) BY INTERPOLATING
C FROM SOUNDING AT PRESSURE LCLP
C
IF(LCLP.LE.PP(NLVL)) THEN
DATA(J)=TTC(NLVL)
ELSE
DATA(J)=TRPT(P,T,N,LCLP)-273.15
ENDIF
C
CASE(002004) ! EQUIL LVL TEMPERATURE
C
C CALCULATE THE PARCEL TEMPERATURE TP() AT EACH
C PRESSURE LEVEL ABOVE THE LCL. THE PARCEL IS
C ASSUMED TO BE SATURATED SO IT IS LIFTED ALONG
C THE SATURATED ADIABAT, THW.
C
DO K=1,NLVL
IF(PP(K).LE.LCLP) THEN
TP(K)=SATLFT(THW,PP(K))
ELSE
TP(K)=TTC(K)
ENDIF
ENDDO
C
C START AT THE TOP PRESSURE LEVEL AND FIND
C WHERE THE PARCEL FIRST BECOMES WARMER THAN
C THE ENVIRONMENT - THIS IS THE EQUIL LEVEL.
C THEN INTERPOLATE FROM THE SOUNDING TO FIND
C THE TEMPERATURE AT THIS LEVEL.
C
IF(TP(NLVL).GT.TTC(NLVL)) THEN
ELP=PP(NLVL)
DATA(J)=TTC(NLVL) ! PARCEL WARMER THAN ENV AT TOP LEVEL
ELSE ! SET EQUIL LEVEL TO THIS PRESSURE
L=NLVL
DO K=NLVL,1,-1
IF(TP(K).LE.TTC(K)) THEN
L=K
P1=PP(K)
ELSE
EXIT
ENDIF
ENDDO
IF(L.EQ.1) THEN ! PARCEL COLDER THAN ENV ABOVE LCL
ELP=LCLP ! SET EQUIL LEVEL TO LCL PRESSURE
DATA(J)=TRPT(P,T,N,ELP)-273.15
ELSE
DO Z=1,26 ! ITERATE TO FIND EQUIL LEVEL
PPP=P1+Z-1
TTE=TRPT(PP,TTC,NLVL,PPP)
TTP=TRPT(PP,TP,NLVL,PPP)
C IF(TTP.GT.TTE) THEN
IF(TTP.GE.TTE) THEN
ELP=PPP-0.5
DATA(J)=TRPT(P,T,N,ELP)-273.15
EXIT
ENDIF
ENDDO
ENDIF
ENDIF
C
CASE(007021) ! BEST LIFTED INDEX
C
C FIND TEMPERATURE OF THE MOST UNSTABLE PARCEL
C AT 500 MB - THIS IS COMPUTED BY FUNCTION SATLFT
C
TP500=SATLFT(THW,500.)
C
C COMPUTE LIFTED INDEX BY SUBTRACTING THE PARCEL
C TEMPERATURE FROM THE ENV TEMPERATURE AT 500 MB
C
DATA(J)=TTC(N500)-TP500
C
CASE(007104) ! LAYER MUCAPE
C
C FIND THE PRESSURE THAT CORRESPONDS TO THE
C TEMPERATURE THRESHOLD GIVEN BY IDPARS(6) -
C THIS WILL DEFINE THE TOP OF THE LAYER FOR
C COMPUTING MUCAPE.
C
IF(IDPARS(6).EQ.0) THEN
PTOP=PP(NLVL) ! DEFAULT TO HIGHEST LEVEL
ELSE
L=NLVL
DO K=NLVL,1,-1 ! START AT THE TOP AND FIND LEVEL WHERE ENV
IF(TT(K).LT.IDPARS(6)) THEN ! TEMP IS WARMER THAN THRESHHOLD
L=K
P1=PP(K)
ELSE
EXIT
ENDIF
ENDDO
IF(L.EQ.1) THEN ! ENTIRE SOUNDING < T. SET PTOP TO SFC PRESSURE
PTOP=PP(1)
ELSE IF(L.EQ.NLVL) THEN ! ENTIRE SOUNDING > T. SET PTOP TO HIGHEST LVL
PTOP=PP(NLVL)
ELSE
DO Z=1,26 ! ITERATE TO FIND P AT T THRESH
PPP=P1+Z-1
TTE=TRPT(PP,TT,NLVL,PPP)
C IF(TTE.GT.IDPARS(6)) THEN
IF(TTE.GE.IDPARS(6)) THEN
PTOP=PPP-0.5
EXIT
ENDIF
ENDDO
ENDIF
ENDIF
C
C FIND THE PRESSURE THAT CORRESPONDS TO 0 C.
C THIS WILL DEFINE THE BOTTOM OF THE LAYER FOR
C COMPUTING MUCAPE.
C
L=NLVL
DO K=NLVL,1,-1 ! START AT THE TOP AND FIND LEVEL WHERE ENV
IF(TT(K).LT.273.15) THEN ! TEMP IS WARMER THAN 0C
L=K
P1=PP(K)
ELSE
EXIT
ENDIF
ENDDO
IF(L.EQ.1) THEN ! ENTIRE SOUNDING < 0C. SET PBOT TO SFC PRESSURE
PBOT=PP(1)
ELSE IF(L.EQ.NLVL) THEN ! ENTIRE SOUNDING > 0C. SET PBOT TO HIGHEST LVL
PBOT=PP(NLVL)
ELSE
DO Z=1,26 ! ITERATE TO FIND P AT 0C
PPP=P1+Z-1
TTE=TRPT(PP,TT,NLVL,PPP)
C IF(TTE.GT.273.15) THEN
IF(TTE.GE.273.15) THEN
PBOT=PPP-0.5
EXIT
ENDIF
ENDDO
ENDIF
C
C CALL ROUTINE POSAREA TO COMPUTE MUCAPE FOR
C IN THE LAYER BOUNDED BY PRESSURE PBOT AT THE
C BOTTOM AND PRESSURE PTOP AT THE TOP.
C
CALL POSAREA(TO,DO,PO,PBOT,PTOP,TTC,DDC,PP,
1 NLVL,DATA(J))
C
CASE(007105) ! LAYER MUCAPE
C
C FIND THE PRESSURE THAT CORRESPONDS TO THE
C TEMPERATURE THRESHOLD GIVEN BY IDPARS(6) -
C THIS WILL DEFINE BOTTOM OF THE LAYER FOR
C COMPUTING MUCAPE.
C
IF(IDPARS(6).EQ.0) THEN
PB=PP(1)
ELSE
L=NLVL
DO K=NLVL,1,-1 ! START AT THE TOP AND FIND LEVEL WHERE ENV
IF(TT(K).LT.IDPARS(6)) THEN ! TEMP IS COLDER THAN THRESHHOLD
L=K
P1=PP(K)
ELSE
EXIT
ENDIF
ENDDO
IF(L.EQ.1) THEN ! ENTIRE SOUNDING < T. SET PB TO SFC PRESSURE
PB=PP(1)
ELSE IF(L.EQ.NLVL) THEN ! ENTIRE SOUNDING > T. SET PB TO HIGHEST LVL
PB=PP(NLVL)
ELSE
DO Z=1,26 ! ITERATE TO FIND P AT T THRESH
PPP=P1+Z-1
TTE=TRPT(PP,TT,NLVL,PPP)
C IF(TTE.GT.IDPARS(6)) THEN
IF(TTE.GE.IDPARS(6)) THEN
PB=PPP-0.5
EXIT
ENDIF
ENDDO
ENDIF
ENDIF
C
C CALL ROUTINE POSAREA TO COMPUTE MOST UNSTABLE
C CAPE IN THE LAYER BOUNDED BY PRESSURE PB AT THE
C BOTTOM AND PRESSURE PP(NLVL) AT THE TOP.
C
CALL POSAREA(TO,DO,PO,PB,PP(NLVL),TTC,DDC,PP,
1 NLVL,DATA(J))
C
CASE(007121) ! LAYER CIN
C
C FIND TEMPERATURE AND DEWPOINT OF ORIGINATING
C PARCEL IN CELSIUS
C
PO=PP(1)-IDPARS(6) ! PRESSURE OF ORIG PARCEL
TO=TRPT(P,T,N,PO)-273.15 ! INTERP T FROM SOUNDING
DO=TRPT(P,D,N,PO)-273.15 ! INTERP TD FROM SOUNDING
C
C CALL ROUTINE NEGAREA TO COMPUTE CIN FOR
C PARCEL ORIGINATING AT PRESSURE PO IN THE
C LAYER BOUNDED BY PRESSURE PT AT THE TOP.
C
IF(PT.GT.PO) THEN
DATA(J)=0.
ELSE
CALL NEGAREA(TO,DO,PO,TTC,DDC,PP,PT,NLVL,DATA(J))
ENDIF
C
END SELECT
C
ENDDO
C
!$OMP END PARALLEL DO
C
C SET NSLAB=1 IN CASE IT GOT CHANGED EARLIER
C
NSLAB=1
C
RETURN
END