From e5bfb051d7476afdc73b00a220d4e1d795fb8860 Mon Sep 17 00:00:00 2001 From: EdwardColon-NOAA Date: Tue, 2 Mar 2021 10:48:23 -0500 Subject: [PATCH 01/14] Update params_grib2_tbl_new --- parm/params_grib2_tbl_new | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/parm/params_grib2_tbl_new b/parm/params_grib2_tbl_new index e84973b74..1cd40828f 100755 --- a/parm/params_grib2_tbl_new +++ b/parm/params_grib2_tbl_new @@ -290,11 +290,14 @@ 3 5 4 0 EBSSTSTD 0 19 30 0 EDPARM 0 7 204 1 EFHL + 0 3 222 0 EFSH + 0 2 233 0 ESHR 0 7 9 0 EHLX 4 2 1 0 ELCDEN 4 0 1 0 ELECTMP 10 3 194 1 ELEV 0 19 238 1 ELLINX + 0 0 205 0 ELMELT 0 191 197 1 ELONN 0 191 193 1 ELON 0 1 211 1 EMNP @@ -1086,6 +1089,7 @@ 0 1 119 0 UCSCLW 0 0 28 0 UCTMP 0 3 29 0 UDRATE + 0 2 234 0 UESH 0 2 17 0 UFLX 0 2 2 0 UGRD 0 2 23 0 UGUST @@ -1135,6 +1139,7 @@ 4 1 1 0 VEL1 4 1 2 0 VEL2 4 1 3 0 VEL3 + 0 2 235 0 VESH 0 2 18 0 VFLX 0 6 48 0 VFRCICE 0 6 49 0 VFRCIW From 2dad12ebf8fdf60feb4d3985939556cd8b974a18 Mon Sep 17 00:00:00 2001 From: EdwardColon-NOAA Date: Tue, 2 Mar 2021 10:59:24 -0500 Subject: [PATCH 02/14] Update post_avblflds.xml --- parm/post_avblflds.xml | 57 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 57 insertions(+) diff --git a/parm/post_avblflds.xml b/parm/post_avblflds.xml index 92dc1414b..d54abc9db 100755 --- a/parm/post_avblflds.xml +++ b/parm/post_avblflds.xml @@ -6989,5 +6989,62 @@ 5.0 + + 979 + EFSH_ON_LFC + EFSH + level_free_convection + 4.0 + + + + 980 + EFSH_ON_EQUIL_LVL + EFSH + equil_lvl + 4.0 + + + + 981 + ELMELT_ON_LFC + ELMELT + level_free_convection + 6.0 + + + + 982 + ELMELT_ON_EQUIL_LVL + ELMELT + equil_lvl + 6.0 + + + + 983 + UESH_ON_LFC + UESH + level_free_convection + 6.0 + + + + 984 + VESH_ON_LFC + VESH + level_free_convection + 6.0 + + + + 985 + ESHR_ON_LFC + ESHR + level_free_convection + 6.0 + + + From 8827d9cd43e3dae7d654dab420b47ab8c5583c4a Mon Sep 17 00:00:00 2001 From: "Edward.Colon" Date: Thu, 5 Aug 2021 23:24:30 +0000 Subject: [PATCH 03/14] This commit includes changes to the computation of the effective layer top and bottom. --- sorc/ncep_post.fd/ALLOCATE_ALL.f | 16 +- sorc/ncep_post.fd/DEALLOCATE.f | 1 + sorc/ncep_post.fd/MISCLN.f | 573 +++++++++++++++++++++++++++---- sorc/ncep_post.fd/TTBLEX.f | 129 ++++++- sorc/ncep_post.fd/UPP_PHYSICS.f | 464 ++++++++++++++++++++++++- sorc/ncep_post.fd/VRBLS2D_mod.f | 5 +- 6 files changed, 1105 insertions(+), 83 deletions(-) diff --git a/sorc/ncep_post.fd/ALLOCATE_ALL.f b/sorc/ncep_post.fd/ALLOCATE_ALL.f index 322ccb8ad..57c6b8033 100644 --- a/sorc/ncep_post.fd/ALLOCATE_ALL.f +++ b/sorc/ncep_post.fd/ALLOCATE_ALL.f @@ -79,7 +79,7 @@ SUBROUTINE ALLOCATE_ALL() !$omp parallel do private(i,j,l) do l=1,lm do j=jsta_2l,jend_2u - do i=1,im + do i=1,lm u(i,j,l)=0. v(i,j,l)=0. t(i,j,l)=spval @@ -107,7 +107,7 @@ SUBROUTINE ALLOCATE_ALL() !$omp parallel do private(i,j,l) do l=1,lp1 do j=jsta_2l,jend_2u - do i=1,im + do i=1,lm pint(i,j,l)=spval alpint(i,j,l)=spval zint(i,j,l)=spval @@ -147,7 +147,7 @@ SUBROUTINE ALLOCATE_ALL() !$omp parallel do private(i,j,l) do l=1,lm do j=jsta_2l,jend_2u - do i=1,im + do i=1,lm cwm(i,j,l)=spval F_ice(i,j,l)=spval F_rain(i,j,l)=spval @@ -195,7 +195,7 @@ SUBROUTINE ALLOCATE_ALL() !$omp parallel do private(i,j,l) do l=1,lm do j=jsta_2l,jend_2u - do i=1,im + do i=1,lm NRAIN(i,j,l)=spval radius_cloud(i,j,l)=spval radius_ice(i,j,l)=spval @@ -669,6 +669,7 @@ SUBROUTINE ALLOCATE_ALL() allocate(z500(im,jsta_2l:jend_2u)) allocate(z700(im,jsta_2l:jend_2u)) allocate(teql(im,jsta_2l:jend_2u)) + allocate(ieql(im,jsta_2l:jend_2u)) allocate(cfracl(im,jsta_2l:jend_2u)) allocate(cfracm(im,jsta_2l:jend_2u)) allocate(cfrach(im,jsta_2l:jend_2u)) @@ -694,6 +695,7 @@ SUBROUTINE ALLOCATE_ALL() t700(i,j)=spval z700(i,j)=spval teql(i,j)=spval + ieql(i,j)=spval cfracl(i,j)=spval cfracm(i,j)=spval cfrach(i,j)=spval @@ -1229,7 +1231,7 @@ SUBROUTINE ALLOCATE_ALL() !Initialization !$omp parallel do private(i,j) do j=jsta_2l,jend_2u - do i=1,im + do i=1,lm dusmass(i,j)=spval ducmass(i,j)=spval dusmass25(i,j)=spval @@ -1271,7 +1273,7 @@ SUBROUTINE ALLOCATE_ALL() !Initialization !$omp parallel do private(i,j) do j=jsta_2l,jend_2u - do i=1,im + do i=1,lm acswupt(i,j)=spval swdnt(i,j)=spval acswdnt(i,j)=spval @@ -1285,7 +1287,7 @@ SUBROUTINE ALLOCATE_ALL() !Initialization !$omp parallel do private(i,j) do j=jsta_2l,jend_2u - do i=1,im + do i=1,lm ddvdx(i,j)=spval ddudy(i,j)=spval uuavg(i,j)=spval diff --git a/sorc/ncep_post.fd/DEALLOCATE.f b/sorc/ncep_post.fd/DEALLOCATE.f index bb80f4496..1d84e1068 100644 --- a/sorc/ncep_post.fd/DEALLOCATE.f +++ b/sorc/ncep_post.fd/DEALLOCATE.f @@ -237,6 +237,7 @@ SUBROUTINE DE_ALLOCATE deallocate(z500) deallocate(z700) deallocate(teql) + deallocate(ieql) deallocate(cfracl) deallocate(cfracm) deallocate(cfrach) diff --git a/sorc/ncep_post.fd/MISCLN.f b/sorc/ncep_post.fd/MISCLN.f index 3538cd00e..4ecbf4569 100644 --- a/sorc/ncep_post.fd/MISCLN.f +++ b/sorc/ncep_post.fd/MISCLN.f @@ -86,7 +86,7 @@ SUBROUTINE MISCLN use vrbls3d, only: pmid, uh, vh, t, zmid, zint, pint, alpint, q, omga use vrbls3d, only: catedr,mwt,gtg use vrbls2d, only: pblh, cprate, fis, T500, T700, Z500, Z700,& - teql + teql,ieql use masks, only: lmh use params_mod, only: d00, d50, h99999, h100, h1, h1m12, pq0, a2, a3, a4, & rhmin, rgamog, tfrz, small, g @@ -95,7 +95,8 @@ SUBROUTINE MISCLN jsta_2l, jend_2u, MODELNAME, SUBMODELNAME use rqstfld_mod, only: iget, lvls, id, iavblfld, lvlsxml use grib2_module, only: pset - use upp_physics, only: FPVSNEW,CALRH_PW,CALCAPE,CALCAPE2,TVIRTUAL + use upp_physics, only: FPVSNEW,CALRH_PW,CALCAPE,CALCAPE2,TVIRTUAL, & + CALCAPE1D use gridspec_mod, only: gridtype !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none @@ -114,6 +115,7 @@ SUBROUTINE MISCLN real,PARAMETER :: D2000=2000 real,PARAMETER :: HCONST=42000000. real,PARAMETER :: K2C=273.16 + REAL,PARAMETER :: DM9999=-9999.0 ! ! DECLARE VARIABLES. @@ -129,7 +131,8 @@ SUBROUTINE MISCLN RH1D, EGRID1, EGRID2, EGRID3, EGRID4, & EGRID5, EGRID6, EGRID7, EGRID8, & MLCAPE,MLCIN,MLLCL,MUCAPE,MUCIN,MUMIXR, & - FREEZELVL,MUQ1D,SLCL + FREEZELVL,MUQ1D,SLCL,THE,MAXTHE + integer,dimension(im,jsta:jend) :: MAXTHEPOS real, dimension(:,:,:),allocatable :: OMGBND, PWTBND, QCNVBND, & PBND, TBND, QBND, & UBND, VBND, RHBND, & @@ -150,12 +153,13 @@ SUBROUTINE MISCLN MAXWP, MAXWZ, MAXWU, MAXWV, & MAXWT INTEGER,dimension(:,:),allocatable :: LLOW, LUPP + INTEGER,dimension(:,:),allocatable :: LLOW1, LUPP1 REAL, dimension(:,:),allocatable :: CANGLE,ESHR,UVECT,VVECT,& EFFUST,EFFVST,FSHR,HTSFC,& ESRH ! integer I,J,jj,L,ITYPE,ISVALUE,LBND,ILVL,IFD,ITYPEFDLVL(NFD), & - iget1, iget2, iget3, LLMH + iget1, iget2, iget3, LLMH,imax,jmax,lmax real DPBND,PKL1,PKU1,FAC1,FAC2,PL,TL,QL,QSAT,RHL,TVRL,TVRBLO, & ES1,ES2,QS1,QS2,RH1,RH2,ZSF,DEPTH(2),work1,work2,work3, & SCINtmp,MUCAPEtmp,MUCINtmp,MLLCLtmp,ESHRtmp,MLCAPEtmp,STP,& @@ -169,11 +173,34 @@ SUBROUTINE MISCLN integer ISTART,ISTOP,JSTART,JSTOP,MIDCAL real dummy(IM,jsta:jend) integer idummy(IM,jsta:jend) +! NEW VARIABLES USED FOR EFFECTIVE LAYER + INTEGER,dimension(:,:),allocatable :: EL_BASE, EL_TOPS + LOGICAL,dimension(:,:),allocatable :: FOUND_BASE, FOUND_TOPS + INTEGER,dimension(:,:),allocatable :: L_THETAE_MAX + INTEGER,dimension(:,:),allocatable :: CAPE9, CINS9 + CHARACTER(LEN=5) :: IM_CH, JSTA_CH, JEND_CH, ME_CH + CHARACTER(LEN=60) :: EFFL_FNAME + CHARACTER(LEN=60) :: EFFL_FNAME2 + INTEGER :: IREC, IUNIT + INTEGER :: IREC2, IUNIT2 + LOGICAL :: debugprint + INTEGER :: EL_SCHEME + INTEGER :: LLL + INTEGER :: LLCL_PAR, LEQL_PAR + REAL :: LMASK, PSFC, CAPE_PAR, CINS_PAR, LPAR0 + REAL, DIMENSION(4) :: PARCEL0 + REAL, DIMENSION(:), ALLOCATABLE :: TPAR_B, TPAR_T + REAL, DIMENSION(:), ALLOCATABLE :: TPAR_TMP + REAL, DIMENSION(:), ALLOCATABLE :: P_AMB, T_AMB, Q_AMB, ZINT_AMB + REAL, DIMENSION(:,:,:), ALLOCATABLE :: TPAR_BASE, TPAR_TOPS ! !**************************************************************************** ! START MISCLN HERE. ! + debugprint = .TRUE. + EL_SCHEME = 1 + allocate(USHR1(IM,jsta_2l:jend_2u),VSHR1(IM,jsta_2l:jend_2u), & USHR6(IM,jsta_2l:jend_2u),VSHR6(IM,jsta_2l:jend_2u)) allocate(UST(IM,jsta_2l:jend_2u),VST(IM,jsta_2l:jend_2u), & @@ -1401,7 +1428,9 @@ SUBROUTINE MISCLN DO J=JSTA,JEND DO I=1,IM GRID1(I,J)=Z1D(I,J) - IF (SUBMODELNAME == 'RTMA') FREEZELVL(I,J)=GRID1(I,J) + IF (SUBMODELNAME == 'RTMA' .OR. MODELNAME == 'FV3R') THEN + FREEZELVL(I,J)=GRID1(I,J) + ENDIF ENDDO ENDDO CALL BOUND (GRID1,D00,H99999) @@ -1490,7 +1519,7 @@ SUBROUTINE MISCLN END IF ! HIGHEST FREEZING LEVEL RELATIVE HUMIDITY - IF (IGET(350)>0)THEN + IF (IGET(350)>0)THEN GRID1=spval !$omp parallel do private(i,j) DO J=JSTA,JEND @@ -1563,7 +1592,7 @@ SUBROUTINE MISCLN END IF ! HIGHEST -10C ISOTHERM RELATIVE HUMIDITY - IF (IGET(777)>0)THEN + IF (IGET(777)>0)THEN GRID1=spval !$omp parallel do private(i,j) DO J=JSTA,JEND @@ -2115,7 +2144,6 @@ SUBROUTINE MISCLN !$omp parallel do private(i,j) DO J=JSTA,JEND DO I=1,IM - IF (EGRID1(I,J) > EGRID2(I,J)) THEN EGRID2(I,J) = EGRID1(I,J) LB2(I,J) = LVLBND(I,J,LBND) @@ -2237,7 +2265,7 @@ SUBROUTINE MISCLN endif ENDIF IF (IGET(110)>0) THEN - GRID1=spval + GRID1=spval !$omp parallel do private(i,j) DO J=JSTA,JEND DO I=1,IM @@ -2793,8 +2821,8 @@ SUBROUTINE MISCLN ! ! SIGMA 0.85000-1.00000 MOISTURE CONVERGENCE. IF (IGET(103)>0) THEN + GRID1=spval ! CONVERT TO DIVERGENCE FOR GRIB - GRID1=spval !$omp parallel do private(i,j) DO J=JSTA,JEND DO I=1,IM @@ -2910,7 +2938,8 @@ SUBROUTINE MISCLN EGRID1(I,J) = LOG(PMID(I,J,LM)/EGRID2(I,J)) & / LOG(PMID(I,J,LM)/PMID(I,J,LM-1)) - IF (MODELNAME == 'GFS' .OR. MODELNAME == 'FV3R') THEN + IF (MODELNAME == 'GFS' .OR. MODELNAME == 'FV3R' & + .OR. MODELNAME == 'FV3R') THEN EGRID1(I,J) = LOG(PMID(I,J,LM)/EGRID2(I,J)) & / max(1.e-6,LOG(PMID(I,J,LM)/PMID(I,J,LM-1))) EGRID1(I,J) =max(-10.0,min(EGRID1(I,J), 10.0)) @@ -2987,7 +3016,7 @@ SUBROUTINE MISCLN QS2 = CON_EPS*ES2/(PMID(I,J,LM-1)+CON_EPSM1*ES2) RH2 = Q(I,J,LM-1)/QS2 GRID1(I,J) = (RH1+(RH2-RH1)*EGRID1(I,J))*100. - ENDIF + ENDIF ENDDO ENDDO CALL BOUND(GRID1,D00,H100) @@ -3059,7 +3088,7 @@ SUBROUTINE MISCLN DO J=JSTA,JEND DO I=1,IM IF(OMGA(I,J,LM)0) THEN + IF (IGET(582)>0) THEN ! dong add missing value for cape GRID1=spval !$omp parallel do private(i,j) @@ -3129,7 +3154,9 @@ SUBROUTINE MISCLN DO I=1,IM IF(T1D(I,J) < spval) THEN GRID1(I,J) = EGRID1(I,J) - IF (SUBMODELNAME == 'RTMA') MLCAPE(I,J)=GRID1(I,J) + IF (SUBMODELNAME == 'RTMA' .OR. MODELNAME == 'FV3R') THEN + MLCAPE(I,J)=GRID1(I,J) + ENDIF ENDIF ENDDO ENDDO @@ -3164,7 +3191,9 @@ SUBROUTINE MISCLN DO I=1,IM IF(T1D(I,J) < spval) THEN GRID1(I,J) = - GRID1(I,J) - IF (SUBMODELNAME == 'RTMA') MLCIN(I,J) = GRID1(I,J) + IF (SUBMODELNAME == 'RTMA' .OR. MODELNAME == 'FV3R') THEN + MLCIN(I,J) = GRID1(I,J) + ENDIF ENDIF ENDDO ENDDO @@ -3194,7 +3223,9 @@ SUBROUTINE MISCLN DO J=JSTA,JEND DO I=1,IM IF(T1D(I,J) < spval) GRID1(I,J)=EGRID2(I,J) - IF (SUBMODELNAME == 'RTMA') MLLCL(I,J) = GRID1(I,J) + IF (SUBMODELNAME == 'RTMA' .OR. MODELNAME == 'FV3R') THEN + MLLCL(I,J) = GRID1(I,J) + ENDIF ENDDO ENDDO ! @@ -3252,7 +3283,9 @@ SUBROUTINE MISCLN DPBND = 300.E2 CALL CALCAPE(ITYPE,DPBND,P1D,T1D,Q1D,LB2,EGRID1, & EGRID2,EGRID3,EGRID4,EGRID5) - IF (SUBMODELNAME == 'RTMA') MUMIXR(I,J) = Q1D(I,J) +! IF (SUBMODELNAME == 'RTMA' .OR. MODELNAME == 'FV3R') THEN +! MUMIXR(I,J) = Q1D(I,J) +! ENDIF IF (IGET(584)>0) THEN ! dong add missing value to cin GRID1 = spval @@ -3261,12 +3294,14 @@ SUBROUTINE MISCLN DO I=1,IM IF(T1D(I,J) < spval) THEN GRID1(I,J) = EGRID1(I,J) - IF (SUBMODELNAME == 'RTMA') MUCAPE(I,J) = GRID1(I,J) + IF (SUBMODELNAME == 'RTMA' .OR. MODELNAME == 'FV3R') THEN + MUCAPE(I,J) = GRID1(I,J) + ENDIF ENDIF ENDDO ENDDO CALL BOUND(GRID1,D00,H99999) -! IF (SUBMODELNAME == 'RTMA') THEN +! IF (SUBMODELNAME == 'RTMA' .OR. MODELNAME == 'FV3R') THEN ! CALL BOUND(MUCAPE,D00,H99999) ! ENDIF if(grib=='grib2') then @@ -3298,7 +3333,8 @@ SUBROUTINE MISCLN DO I=1,IM IF(T1D(I,J) < spval) THEN GRID1(I,J) = - GRID1(I,J) - IF (SUBMODELNAME == 'RTMA') THEN + IF (SUBMODELNAME == 'RTMA' .OR. & + MODELNAME == 'FV3R') THEN MUCAPE(I,J) = GRID1(I,J) MUQ1D(I,J) = Q1D(I,J) ENDIF @@ -3342,7 +3378,6 @@ SUBROUTINE MISCLN enddo endif ENDIF - !Equilibrium Temperature IF (IGET(982)>0) THEN DO J=JSTA,JEND @@ -3394,7 +3429,7 @@ SUBROUTINE MISCLN ! GENERAL THUNDER PARAMETER ??? 458 ??? IF (IGET(444)>0) THEN - GRID1 = spval + GRID1 = spval !$omp parallel do private(i,j) DO J=JSTA,JEND DO I=1,IM @@ -3422,11 +3457,361 @@ SUBROUTINE MISCLN endif ENDIF ENDIF + +! +! --- Effective (inflow) Layer (EL) +! + IF ( EL_SCHEME > 0 ) THEN + ALLOCATE(EL_BASE(IM,JSTA_2L:JEND_2U)) + ALLOCATE(EL_TOPS(IM,JSTA_2L:JEND_2U)) + ALLOCATE(FOUND_BASE(IM,JSTA_2L:JEND_2U)) + ALLOCATE(FOUND_TOPS(IM,JSTA_2L:JEND_2U)) +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=1,IM + EL_BASE(I,J) = LM + EL_TOPS(I,J) = LM + FOUND_BASE(I,J) = .FALSE. + FOUND_TOPS(I,J) = .FALSE. + ENDDO + ENDDO + END IF +! + IF ( EL_SCHEME == 1 ) THEN + ITYPE = 2 + DPBND = 0. + + DO L = LM, 1, -1 + +! SET AIR PARCELS FOR LEVEL L +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=1,IM + EGRID1(I,J) = -H99999 + EGRID2(I,J) = -H99999 + IDUMMY(I,J) = 0 + P1D(I,J) = PMID(I,J,L) + T1D(I,J) = T(I,J,L) + Q1D(I,J) = Q(I,J,L) + ENDDO + ENDDO + +!--- CALCULATE CAPE/CIN FOR ALL AIR PARCELS on LEVEL L + IF (debugprint) WRITE(1000+ME,'(1x,A,I2.2,2x,A,I3)') & + 'NEW EL_SCHEME:', EL_SCHEME, & + ' CALCULATING CAPE/CINS ON LEVEL:',L + CALL CALCAPE(ITYPE,DPBND,P1D,T1D,Q1D,IDUMMY,EGRID1, & + EGRID2,EGRID3,EGRID4,EGRID5) + +!--- CHECK CAPE/CIN OF EACH AIR PARCELS WITH EL CRITERIA +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=1,IM + IF ( .NOT. FOUND_BASE(I,J) ) THEN + IF ( EGRID1(I,J) >= 100. .AND. EGRID2(I,J) >= -250. ) THEN + EL_BASE(I,J) = L + FOUND_BASE(I,J) = .TRUE. + ELSE + EL_BASE(I,J) = LM + FOUND_BASE(I,J) = .FALSE. + END IF + ELSE + IF ( .NOT. FOUND_TOPS(I,J) ) THEN + IF ( EGRID1(I,J) < 100. .OR. EGRID2(I,J) < -250. ) THEN + EL_TOPS(I,J) = L + 1 + FOUND_TOPS(I,J) = .TRUE. + ELSE + EL_TOPS(I,J) = LM + FOUND_TOPS(I,J) = .FALSE. + END IF + END IF + END IF + ENDDO + ENDDO + + END DO ! L +! + ELSE IF ( EL_SCHEME == 2 ) THEN + +!--- SEARCH FOR EL ALONG EACH PROFILE/COLUMN + ALLOCATE(L_THETAE_MAX(IM,JSTA:JEND)) +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=1,IM + EGRID1(I,J) = -H99999 + EGRID2(I,J) = -H99999 + L_THETAE_MAX(I,J) = -9999 + ENDDO + ENDDO + +!------ SEARCH FOR PARCEL WITH MAX THETA-E OF EACH PROFILE/COLUMN + DO L = LM,1,-1 +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=1,IM + EGRID1(I,J) = -H99999 + P1D(I,J) = PMID(I,J,L) + T1D(I,J) = T(I,J,L) + Q1D(I,J) = Q(I,J,L) + ENDDO + ENDDO +! CALL CALTHTE(PMID(1,jsta,L),T(1,jsta,L),Q(1,jsta,L),EGRID1) + CALL CALTHTE(P1D,T1D,Q1D,EGRID1) +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=1,IM + IF (EGRID1(I,J) > EGRID2(I,J)) THEN + EGRID2(I,J) = EGRID1(I,J) + L_THETAE_MAX(I,J) = L + END IF + END DO + END DO + + END DO ! L + +!--- SET UP THE PARCELS WIH MAX THETA-E + ALLOCATE(CAPE9(IM,JSTA:JEND)) + ALLOCATE(CINS9(IM,JSTA:JEND)) +!$omp parallel do private(i,j,LLL) + DO J=JSTA,JEND + DO I=1,IM + LLL = L_THETAE_MAX(I,J) + P1D(I,J) = PMID(I,J,LLL) + T1D(I,J) = T(I,J,LLL) + Q1D(I,J) = Q(I,J,LLL) + ENDDO + ENDDO +!--- COMPUTE CAPE/CIN OF PARCEL WITH MAX THETA-E + ITYPE = 2 + DPBND = 0. +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=1,IM + EGRID1(I,J) = -H99999 + EGRID2(I,J) = -H99999 + CAPE9(I,J) = -H99999 + CINS9(I,J) = -H99999 + IDUMMY(I,J) = 0 + ENDDO + ENDDO + CALL CALCAPE(ITYPE,DPBND,P1D,T1D,Q1D,IDUMMY,EGRID1, & + EGRID2,EGRID3,EGRID4,EGRID5) +!--- SANITY CHECK IF CAPE/CIN OF MAX THETA-E PARCEL OF THIS COLUMN +!$omp parallel do private(i,j,LLL) + DO J=JSTA,JEND + DO I=1,IM + CAPE9(I,J) = EGRID1(I,J) + CINS9(I,J) = EGRID2(I,j) + ENDDO + ENDDO + + ALLOCATE(TPAR_B(LM), TPAR_T(LM)) + ALLOCATE(TPAR_TMP(LM)) + ALLOCATE(P_AMB(LM), T_AMB(LM), Q_AMB(LM), ZINT_AMB(LM)) + IF (debugprint) THEN + ALLOCATE(TPAR_BASE(IM,JSTA:JEND,LM)) + ALLOCATE(TPAR_TOPS(IM,JSTA:JEND,LM)) +!$omp parallel do private(i,j,l) + DO L=1,LM + DO J=JSTA,JEND + DO I=1,IM + TPAR_BASE(I,J,L) = DM9999 + TPAR_TOPS(I,J,L) = DM9999 + ENDDO + ENDDO + ENDDO + END IF + +!$OMP PARALLEL DO PRIVATE(i,j,l,P_AMB,T_AMB,Q_AMB,ZINT_AMB,LMASK,PSFC, & +!$OMP& TPAR_B,TPAR_T,LPAR0,CAPE_PAR,CINS_PAR, & +!$OMP& PARCEL0,LLCL_PAR,LEQL_PAR,TPAR_TMP) + DO J=JSTA,JEND + DO I=1,IM + +!-----------GET THE AMBIENT PROFILE + DO L=1,LM + P_AMB(L) = PMID(I,J,L) + T_AMB(L) = T(I,J,L) + Q_AMB(L) = Q(I,J,L) + ZINT_AMB(L) = ZINT(I,J,L) + END DO + LMASK = NINT(LMH(I,J)) + PSFC = PMID(I,J,NINT(LMH(I,J))) + EL_BASE(I,J)=LM + EL_TOPS(I,J)=LM + FOUND_BASE(I,J) = .FALSE. + FOUND_TOPS(I,J) = .FALSE. + TPAR_B(1:LM) = DM9999 + TPAR_T(1:LM) = DM9999 + TPAR_TMP(1:LM) = DM9999 + +!---------- SANITY CHECK +!---------- USING CAPE/CINS OF AIR PARCEL WITH MAX THETA-E FOR SANITY CHECK FIRST +!---------- CAPE4/CINS4 ARE CAPE/CINS CALCULATED FROM CODE ABOVE FOR AIR PARCEL +!---------- WITH MAX THETA-E. + IF (CAPE9(I,J) >= 100. .OR. CINS9(I,J) >= -250.) THEN !EL FOR THIS COLUMN +! EL SHOULD EXIST ALONG THIS COLUMN/PROFILE. +!-------------NEED TO SEARCH FROM BOTTOM TO TOP, COMPUTE CAPE/CINS OF AIR PARCEL +! INITIATING AT EACH LEVEL UNTIL THE CRITERIA MEETS. + + VLOOP: DO L=LM,1,-1 +!---------------CALCULATE CAPE/CIN OF AIR PARCEL INITIALIZED AT LEVEL L + IF (.NOT. FOUND_BASE(I,J)) THEN !SEARCH FOR BASE FIRST + LPAR0 = FLOAT(L) + CAPE_PAR = D00 + CINS_PAR = D00 + PARCEL0(1) = LPAR0 + PARCEL0(2:4) = D00 + LLCL_PAR = 1 + LEQL_PAR = LM + TPAR_B(1:LM) = DM9999 + TPAR_T(1:LM) = DM9999 + + WRITE(1000+me,*)"CALCAPE2_1D_base: parcel@I J L:",I,J,L + + CALL CALCAPE1D(P_AMB,T_AMB,Q_AMB,ZINT_AMB,LPAR0, & + PSFC,LMASK,CAPE_PAR,CINS_PAR,TPAR_B, & + PARCEL0,LLCL_PAR,LEQL_PAR) + IF (CAPE_PAR >= 100. .AND. CINS_PAR >= -250.) THEN ! BASE OF EL + EL_BASE(I,J) = L + FOUND_BASE(I,J) = .TRUE. + ELSE + TPAR_B(1:LM) = DM9999 ! RESET TPAR IF NOT FOUND BASE + ENDIF + ELSE !SEARCH FOR TOP AFTER BASE IS FOUND + IF (.NOT. FOUND_TOPS(I,J)) THEN !SEARCH FOR BASE FIRST + LPAR0 = FLOAT(L) + CAPE_PAR = D00 + CINS_PAR = D00 + TPAR_T(1:LM) = DM9999 + PARCEL0(1) = LPAR0 + PARCEL0(2:4) = D00 + LLCL_PAR = 1 + LEQL_PAR = LM + + WRITE(1000+me,*)"CALCAPE2_1D_top: parcel@I J L:",I,J,L + + CALL CALCAPE1D(P_AMB,T_AMB,Q_AMB,ZINT_AMB,LPAR0, & + PSFC,LMASK,CAPE_PAR,CINS_PAR,TPAR_T, & + PARCEL0,LLCL_PAR,LEQL_PAR) + IF (CAPE_PAR < 100. .OR. CINS_PAR < -250.) THEN ! TOP OF EL + EL_TOPS(I,J) = L+1 + FOUND_TOPS(I,J) = .TRUE. + IF (EL_TOPS(I,J) == EL_BASE(I,J)) THEN + TPAR_T(1:LM) = TPAR_B(1:LM) + ELSE IF (EL_TOPS(I,J) < EL_BASE(I,J) ) THEN + TPAR_T(1:LM) = TPAR_TMP(1:LM) + ELSE + WRITE(0,'(1x,A,A)') "TOP OF EFFECTIVE LAYER IS", & + " LOWER THAN BASE. WRONG! ABORT ..." + STOP 9 + END IF + EXIT VLOOP + ELSE + TPAR_TMP(1:LM) = TPAR_T(1:LM) + TPAR_T(1:LM) = DM9999 ! RESET TPAR IF NOT FOUND TOP + ENDIF + ENDIF ! FOUND_TOPS + + ENDIF ! FOUND_BASE OR NOT + + ENDDO VLOOP + + IF (debugprint) THEN + DO L=1,LM + TPAR_BASE(I,J,L) = TPAR_B(L) + TPAR_TOPS(I,J,L) = TPAR_T(L) + END DO + END IF + + ENDIF ! IF PASSING SANITY CHECK + + IF ( FOUND_BASE(I,J) /= FOUND_TOPS(I,J) ) THEN + WRITE(0,'(1x,A,A,A,I6,1x,I6)') "BASE & TOP OF ", & + " EFFECTIVE LAYER ARE NOT FOUND TOGETHER. WRONG! ", & + " ABORT! ABORT! ... AT GRID POINT: ", I, J + STOP 10 + END IF + + ENDDO ! J + ENDDO ! I + + IF(ALLOCATED(L_THETAE_MAX)) DEALLOCATE(L_THETAE_MAX) + IF(ALLOCATED(CAPE9)) DEALLOCATE(CAPE9) + IF(ALLOCATED(CINS9)) DEALLOCATE(CINS9) + IF(ALLOCATED(TPAR_B)) DEALLOCATE(TPAR_B) + IF(ALLOCATED(TPAR_T)) DEALLOCATE(TPAR_T) + IF(ALLOCATED(TPAR_TMP)) DEALLOCATE(TPAR_TMP) + IF(ALLOCATED(P_AMB)) DEALLOCATE(P_AMB) + IF(ALLOCATED(T_AMB)) DEALLOCATE(T_AMB) + IF(ALLOCATED(Q_AMB)) DEALLOCATE(Q_AMB) + IF(ALLOCATED(ZINT_AMB)) DEALLOCATE(ZINT_AMB) + + END IF ! EL_SCHEME :1 OR 2 + + IF (ALLOCATED(FOUND_BASE)) DEALLOCATE(FOUND_BASE) + IF (ALLOCATED(FOUND_TOPS)) DEALLOCATE(FOUND_TOPS) + + IF (debugprint .AND. EL_SCHEME > 0) THEN + WRITE(IM_CH,'(I5.5)') IM + WRITE(JSTA_CH,'(I5.5)') JSTA + WRITE(JEND_CH,'(I5.5)') JEND + EFFL_FNAME="EFFL_NEW_"//IM_CH//"_"//JSTA_CH//"_"//JEND_CH & + //".dat" + EFFL_FNAME2="EFFL_NEW_LVLS_"//IM_CH//"_"//JSTA_CH//"_"//JEND_CH & + //".dat" + IUNIT=10000+JSTA + IUNIT2=20000+JSTA + IREC=0 + IREC2=0 + OPEN(IUNIT,FILE=TRIM(ADJUSTL(EFFL_FNAME)),FORM='FORMATTED') + IF (EL_SCHEME == 2) THEN + OPEN(IUNIT2,FILE=TRIM(ADJUSTL(EFFL_FNAME2)),FORM='FORMATTED') + END IF +! OPEN(IUNIT,FILE=TRIM(ADJUSTL(EFFL_FNAME)),FORM='UNFORMATTED', & +! ACCESS='DIRECT',RECL=4*6) + DO J=JSTA,JEND + DO I=1,IM + IREC = IREC + 1 + IREC2 = IREC2 + 1 +! WRITE(IUNIT,'(1x,I6,2x,I6,2x,I6,2x,I6)')I,J,EL_BASE(I,J),EL_TOPS(I,J) + WRITE(IUNIT,'(1x,I6,2x,I6,2(2x,I6,2x,F12.3))') I, J, & + EL_BASE(I,J),PMID(I,J,EL_BASE(I,J)), & + EL_TOPS(I,J),PMID(I,J,EL_TOPS(I,J)) +! WRITE(IUNIT,REC=IREC) I, J, & +! EL_BASE(I,J),PMID(I,J,EL_BASE(I,J)), & +! EL_TOPS(I,J),PMID(I,J,EL_TOPS(I,J)) + IF (EL_SCHEME == 2) THEN + WRITE(IUNIT2,'(1x,I6,2x,I6,2(2x,I6,2x,F12.3))') I, J, & + EL_BASE(I,J),PMID(I,J,EL_BASE(I,J)), & + EL_TOPS(I,J),PMID(I,J,EL_TOPS(I,J)) + DO L=LM,1,-1 + IREC2=IREC2+1 + WRITE(IUNIT2,'(1x,I4,5(2x,F12.3))') & + LM+1-L, PMID(I,J,L), ZINT(I,J,L), T(I,J,L), & + TPAR_BASE(I,J,L), TPAR_TOPS(I,J,L) +! WRITE(IUNIT2,REC=IREC2) & +! LM+1-L, PMID(I,J,L), ZINT(I,J,L), T(I,J,L), & +! TPAR_BASE(I,J,L), TPAR_TOPS(I,J,L) + END DO + END IF + END DO + ENDDO + CLOSE(IUNIT) + + IF (EL_SCHEME == 2) THEN + CLOSE(IUNIT2) + END IF + + ENDIF + + IF(ALLOCATED(TPAR_BASE)) DEALLOCATE(TPAR_BASE) + IF(ALLOCATED(TPAR_TOPS)) DEALLOCATE(TPAR_TOPS) ! ! EXPAND HRRR CAPE/CIN RELATED VARIABLES ! ! CAPE AND CINS 0-3KM, FOLLOW ML PROCEDURE WITH HEIGHT 0-3KM - +! FIELD1=.FALSE. FIELD2=.FALSE. ! @@ -3444,8 +3829,12 @@ SUBROUTINE MISCLN FIELD2=.TRUE. ENDIF ! +! IF(FIELD1)ITYPE=2 +! IF(FIELD2)ITYPE=2 + IF(FIELD1.OR.FIELD2)THEN ITYPE = 2 + ! !$omp parallel do private(i,j) DO J=JSTA,JEND @@ -3469,12 +3858,12 @@ SUBROUTINE MISCLN Q1D(I,J) = (QBND(I,J,1) + QBND(I,J,2) + QBND(I,J,3))/3 ENDDO ENDDO -! + DPBND = 0. CALL CALCAPE2(ITYPE,DPBND,P1D,T1D,Q1D,LB2, & EGRID1,EGRID2,EGRID3,EGRID4,EGRID5, & EGRID6,EGRID7,EGRID8) -! + ! CAPE1, CINS2, LFC3, ESRHL4,ESRHH5, ! DCAPE6,DGLD7, ESP8) ! @@ -3538,7 +3927,7 @@ SUBROUTINE MISCLN ! LFC HEIGHT IF (IGET(952)>0) THEN - GRID1=spval + GRID1=spval !$omp parallel do private(i,j) DO J=JSTA,JEND DO I=1,IM @@ -3560,12 +3949,14 @@ SUBROUTINE MISCLN endif ENDIF !952 + ! EFFECTIVE STORM RELATIVE HELICITY AND STORM MOTION. allocate(UST(IM,jsta_2l:jend_2u),VST(IM,jsta_2l:jend_2u), & HELI(IM,jsta_2l:jend_2u,2)) allocate(LLOW(IM,jsta_2l:jend_2u),LUPP(IM,jsta_2l:jend_2u), & CANGLE(IM,jsta_2l:jend_2u)) +! allocate(LLOW1(IM,jsta_2l:jend_2u),LUPP1(IM,jsta_2l:jend_2u)) iget1 = IGET(953) iget2 = -1 @@ -3585,6 +3976,38 @@ SUBROUTINE MISCLN LUPP(I,J) = INT(EGRID5(I,J)) ENDDO ENDDO +!--- OUTPUT EL BASE & TOP COMPUTED BY OLD SCHEME + IF (debugprint) THEN + WRITE(IM_CH,'(I5.5)') IM + WRITE(JSTA_CH,'(I5.5)') JSTA + WRITE(JEND_CH,'(I5.5)') JEND + EFFL_FNAME="EFFL_OLD_"//IM_CH//"_"//JSTA_CH//"_"//JEND_CH & + //".dat" + IUNIT=10000+JSTA + IREC=0 + OPEN(IUNIT,FILE=TRIM(ADJUSTL(EFFL_FNAME)),FORM='FORMATTED') + DO J=JSTA,JEND + DO I=1,IM + IREC = IREC + 1 +! WRITE(IUNIT,'(1x,I6,2x,I6,2x,I6,2x,I6)')I,J,LLOW(I,J),LUPP(I,J) + WRITE(IUNIT,'(1x,I6,2x,I6,2(2x,I6,2x,F12.3))') I, J, & + LLOW(I,J),PMID(I,J,LLOW(I,J)), & + LUPP(I,J),PMID(I,J,LUPP(I,J)) + END DO + ENDDO + CLOSE(IUNIT) + ENDIF + +!--- IF USSING EL BASE & TOP COMPUTED BY NEW SCHEME FOR THE RELATED VARIABLES + IF ( EL_SCHEME > 0 ) THEN +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=1,IM + LLOW(I,J) = EL_BASE(I,J) + LUPP(I,J) = EL_TOPS(I,J) + ENDDO + ENDDO + END IF ! CALL CALHEL(DEPTH,UST,VST,HELI,USHR1,VSHR1,USHR6,VSHR6) CALL CALHEL2(LLOW,LUPP,DEPTH,UST,VST,HELI,CANGLE) @@ -3613,7 +4036,8 @@ SUBROUTINE MISCLN ENDIF !953 - IF (SUBMODELNAME == 'RTMA') THEN !Start RTMA block + + IF (SUBMODELNAME == 'RTMA' .OR. MODELNAME=='FV3R') THEN !Start RTMA block !EL field allocation @@ -3622,6 +4046,40 @@ SUBROUTINE MISCLN allocate(EFFUST(IM,jsta_2l:jend_2u),EFFVST(IM,jsta_2l:jend_2u),& ESRH(IM,jsta_2l:jend_2u)) +! + DO J=JSTA,JEND + DO I=1,IM + MAXTHE(I,J)=-H99999 + THE(I,J)=-H99999 + MAXTHEPOS(I,J)=0 + MUQ1D(I,J) = 0. + ENDDO + ENDDO + DO L=LM,1,-1 + + DO J=JSTA,JEND + DO I=1,IM + EGRID1(I,J) = -H99999 + P1D(I,J)=PMID(I,J,L) + T1D(I,J)=T(I,J,L) + Q1D(I,J)=Q(I,J,L) + ENDDO + ENDDO + CALL CALTHTE(P1D,T1D,Q1D,EGRID1) + DO J=JSTA,JEND + DO I=1,IM + THE(I,J)=EGRID1(I,J) + IF(THE(I,J)>=MAXTHE(I,J))THEN + MAXTHE(I,J)=THE(I,J) + MAXTHEPOS(I,J)=L + MUQ1D(I,J) = Q(I,J,L) ! save the Q of air parcel with max theta-e (MU Parcel) + ENDIF + ENDDO + ENDDO + + ENDDO + +! !get surface height IF(gridtype == 'E')THEN JVN = 1 @@ -3725,14 +4183,14 @@ SUBROUTINE MISCLN DO J=JSTA,JEND DO I=1,IM IF(LLOW(I,J)0) THEN GRID1(I,J)=STP @@ -3934,7 +4392,7 @@ SUBROUTINE MISCLN !Fixed Layer Tornado Parameter IF (IGET(990)>0) THEN - DO J=JSTA,JEND + DO J=JSTA,JEND DO I=1,IM LLMH = NINT(LMH(I,J)) P1D(I,J) = PMID(I,J,LLMH) @@ -4077,6 +4535,7 @@ SUBROUTINE MISCLN CALL CALCAPE2(ITYPE,DPBND,P1D,T1D,Q1D,LB2, & EGRID1,EGRID2,EGRID3,EGRID4,EGRID5, & EGRID6,EGRID7,EGRID8) + GRID1=spval DO J=JSTA,JEND DO I=1,IM @@ -4126,9 +4585,6 @@ SUBROUTINE MISCLN GRID1=spval DO J=JSTA,JEND DO I=1,IM - IF(T700(I,J) < spval .and. T500(I,J) < spval .and.& - Z700(I,J) < spval .and. Z500(I,J) < spval .and.& - MUCAPE(I,J) < spval .and. MUQ1D(I,J) < spval .and. FSHR(I,J) < spval) THEN LAPSE=-((T700(I,J)-T500(I,J))/((Z700(I,J)-Z500(I,J)))) SHIP=(MUCAPE(I,J)*D1000*MUQ1D(I,J)*LAPSE*(T500(I,J)-K2C)*FSHR(I,J))/HCONST IF (MUCAPE(I,J)<1300.)THEN @@ -4141,7 +4597,6 @@ SUBROUTINE MISCLN SHIP=SHIP*(FREEZELVL(I,J)/2400.) ENDIF GRID1(I,J)=SHIP - ENDIF ENDDO ENDDO if(grib=='grib2') then @@ -4168,7 +4623,7 @@ SUBROUTINE MISCLN !$omp parallel do private(i,j) DO J=JSTA,JEND DO I=1,IM - IF(T1D(I,J) < spval ) GRID1(I,J) = CANGLE(I,J) + IF(T1D(I,J) < spval ) GRID1(I,J) = CANGLE(I,J) ! IF(EGRID1(I,J)<100. .OR. EGRID2(I,J)>-250.) THEN ! GRID1(I,J) = 0. ! ENDIF @@ -4240,21 +4695,21 @@ SUBROUTINE MISCLN ! Downdraft CAPE - ITYPE = 1 - ! DO J=JSTA,JEND - ! DO I=1,IM - ! LB2(I,J) = (LVLBND(I,J,1) + LVLBND(I,J,2) + & - ! LVLBND(I,J,3))/3 - ! P1D(I,J) = (PBND(I,J,1) + PBND(I,J,2) + PBND(I,J,3))/3 - ! T1D(I,J) = (TBND(I,J,1) + TBND(I,J,2) + TBND(I,J,3))/3 - ! Q1D(I,J) = (QBND(I,J,1) + QBND(I,J,2) + QBND(I,J,3))/3 - ! ENDDO - ! ENDDO - - DPBND = 400.E2 - ! CALL CALCAPE2(ITYPE,DPBND,P1D,T1D,Q1D,LB2, & - ! EGRID1,EGRID2,EGRID3,EGRID4,EGRID5, & - ! EGRID6,EGRID7,EGRID8) +! ITYPE = 1 +! DO J=JSTA,JEND +! DO I=1,IM +! LB2(I,J) = (LVLBND(I,J,1) + LVLBND(I,J,2) + & +! LVLBND(I,J,3))/3 +! P1D(I,J) = (PBND(I,J,1) + PBND(I,J,2) + PBND(I,J,3))/3 +! T1D(I,J) = (TBND(I,J,1) + TBND(I,J,2) + TBND(I,J,3))/3 +! Q1D(I,J) = (QBND(I,J,1) + QBND(I,J,2) + QBND(I,J,3))/3 +! ENDDO +! ENDDO + +! DPBND = 400.E2 +! CALL CALCAPE2(ITYPE,DPBND,P1D,T1D,Q1D,LB2, & +! EGRID1,EGRID2,EGRID3,EGRID4,EGRID5, & +! EGRID6,EGRID7,EGRID8) IF (IGET(954)>0) THEN GRID1 = spval @@ -4289,6 +4744,8 @@ SUBROUTINE MISCLN if (allocated(heli)) deallocate(heli) if (allocated(llow)) deallocate(llow) if (allocated(lupp)) deallocate(lupp) + if (allocated(llow1)) deallocate(llow1) + if (allocated(lupp1)) deallocate(lupp1) if (allocated(cangle))deallocate(cangle) if (allocated(effust))deallocate(effust) if (allocated(effvst))deallocate(effvst) diff --git a/sorc/ncep_post.fd/TTBLEX.f b/sorc/ncep_post.fd/TTBLEX.f index 09bd3423b..0e11d5da5 100644 --- a/sorc/ncep_post.fd/TTBLEX.f +++ b/sorc/ncep_post.fd/TTBLEX.f @@ -112,5 +112,132 @@ SUBROUTINE TTBLEX(TREF,TTBL,ITB,JTB,KARR,PMIDL & ENDDO ! RETURN - END + END SUBROUTINE TTBLEX !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& +! +!------------------------------------------------------------------------------------- +! + SUBROUTINE TTBLEX1D(TREF,TTBL,ITB,JTB,KARR,PMIDL, & + PL,QQ,PP,RDP,THE0,STHE,RDTHE,THESP, & + IPTB,ITHTB) +!FPP$ NOCONCUR R +!$$$ SUBPROGRAM DOCUMENTATION BLOCK +! . . . +! SUBPROGRAM: TTBLEX COMPUTES T ALONG A MOIST ADIABAT +! PRGRMMR: BLACK ORG: W/NP2 DATE: ??-??-?? +! +! ABSTRACT: +! THIS ROUTINE COMPUTES THE TEMPERATURE ALONG A MOIST +! ADIABAT GIVEN THE SATURATION POTENTIAL TEMPERATURE +! AND THE PRESSURE +! . +! +! PROGRAM HISTORY LOG: +! ??-??-?? T BLACK - ORIGINATOR +! 98-06-12 T BLACK - CONVERSION FROM 1-D TO 2-D +! 00-01-04 JIM TUCCILLO - MPI VERSION +! 01-10-22 H CHUANG - MODIFIED TO PROCESS HYBRID MODEL OUTPUT +! 02-01-15 MIKE BALDWIN - WRF VERSION +! +! OUTPUT FILES: +! NONE +! +! SUBPROGRAMS CALLED: +! UTILITIES: +! NONE +! +! ATTRIBUTES: +! LANGUAGE: FORTRAN +!---------------------------------------------------------------------- + use ctlblk_mod, only: jsta, jend, im, jsta_2l, jend_2u, me +!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + implicit none +!---------------------------------------------------------------------- + + integer,intent(in) :: ITB,JTB + integer,intent(in) :: KARR + real,dimension(JTB,ITB),intent(in) :: TTBL + real, intent(in) :: PMIDL + real, intent(out) :: TREF +! real,dimension(IM,jsta:jend),intent(out) :: QQ,PP + real,intent(out) :: QQ,PP + real, intent(in) :: THESP + real,dimension(ITB), intent(in) :: THE0,STHE +! integer,dimension(IM,jsta:jend),intent(out) :: IPTB,ITHTB + integer,intent(out) :: IPTB,ITHTB + real,intent(in) :: PL,RDP,RDTHE + +! + integer I,J,ITH,IP,IPTBK + real PK,TPK,T00K,T10K,T01K,T11K,BTHE00K,STHE00K,BTHK,STHK, & + TTHK,BTHE10K,STHE10K +! real QQ, PP +! INTEGER IPTB, ITHTB +!----------------------------------------------------------------------- + +!!$omp parallel do & +!!$omp& private(i,j,bthe00k,bthe10k,bthk,ip,iptbk,ith,pk,sthe00k,sthe10k,& +!!$omp& sthk,t00k,t01k,t10k,t11k,tpk,tthk) +! DO J=JSTA,JEND +! DO I=1,IM + IF(KARR > 0) THEN +!--------------SCALING PRESSURE & TT TABLE INDEX------------------------ + PK = PMIDL + TPK = (PK-PL)*RDP + QQ = TPK-AINT(TPK) + IPTB = INT(TPK) + 1 +!--------------KEEPING INDICES WITHIN THE TABLE------------------------- + IF(IPTB < 1) THEN + IPTB = 1 + QQ = 0. + ENDIF +! + IF(IPTB >= ITB) THEN + IPTB = ITB-1 + QQ = 0. + ENDIF +!--------------BASE AND SCALING FACTOR FOR THE-------------------------- + IPTBK = IPTB + BTHE00K = THE0(IPTBK) + STHE00K = STHE(IPTBK) + BTHE10K = THE0(IPTBK+1) + STHE10K = STHE(IPTBK+1) +!--------------SCALING THE & TT TABLE INDEX----------------------------- + BTHK = (BTHE10K-BTHE00K)*QQ+BTHE00K + STHK = (STHE10K-STHE00K)*QQ+STHE00K + TTHK = (THESP-BTHK)/STHK*RDTHE + PP = TTHK-AINT(TTHK) +! write(1000+me,*)' i=',i,' j=',j,' tthk=',tthk,' thesp=',thesp(i,j) & +! , ' bthk=',bthk,' sthk=',sthk,' rdthe=',rdthe + + ITHTB = INT(TTHK)+1 +!--------------KEEPING INDICES WITHIN THE TABLE------------------------- + IF(ITHTB < 1) THEN + ITHTB = 1 + PP = 0. + ENDIF +! + IF(ITHTB >= JTB) THEN + ITHTB = JTB-1 + PP = 0. + ENDIF +!--------------TEMPERATURE AT FOUR SURROUNDING TT TABLE PTS.------------ + ITH = ITHTB + IP = IPTB + T00K = TTBL(ITH ,IP ) + T10K = TTBL(ITH+1,IP ) + T01K = TTBL(ITH ,IP+1) + T11K = TTBL(ITH+1,IP+1) +!--------------PARCEL TEMPERATURE------------------------------------- + TREF = (T00K+(T10K-T00K)*PP+(T01K-T00K)*QQ & + + (T00K-T10K-T01K+T11K)*PP*QQ) + ENDIF +! ENDDO +! ENDDO +! + RETURN + END SUBROUTINE TTBLEX1D +!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& +! +!------------------------------------------------------------------------------------- +! diff --git a/sorc/ncep_post.fd/UPP_PHYSICS.f b/sorc/ncep_post.fd/UPP_PHYSICS.f index b9ca9995d..684de251a 100644 --- a/sorc/ncep_post.fd/UPP_PHYSICS.f +++ b/sorc/ncep_post.fd/UPP_PHYSICS.f @@ -38,7 +38,7 @@ module upp_physics private - public :: CALCAPE, CALCAPE2 + public :: CALCAPE, CALCAPE2, CALCAPE1D public :: CALRH public :: CALRH_GFS, CALRH_GSD, CALRH_NAM public :: CALRH_PW @@ -51,7 +51,7 @@ module upp_physics ! SUBROUTINE CALRH(P1,T1,Q1,RH) - use ctlblk_mod, only: im, jsta, jend, MODELNAME + use ctlblk_mod, only: im, jsta, jend, MODELNAME, spval implicit none REAL,dimension(IM,jsta:jend),intent(in) :: P1,T1 @@ -137,7 +137,7 @@ SUBROUTINE CALRH_NAM(P1,T1,Q1,RH) ! DO J=JSTA,JEND DO I=1,IM - IF (T1(I,J) < SPVAL) THEN + IF (T1(I,J) < spval) THEN IF (ABS(P1(I,J)) >= 1) THEN QC = PQ0/P1(I,J)*EXP(A2*(T1(I,J)-A3)/(T1(I,J)-A4)) ! @@ -156,7 +156,7 @@ SUBROUTINE CALRH_NAM(P1,T1,Q1,RH) ! ENDIF ELSE - RH(I,J) = SPVAL + RH(I,J) = spval ENDIF ENDDO ENDDO @@ -245,7 +245,7 @@ SUBROUTINE CALRH_GFS(P1,T1,Q1,RH) !$omp parallel do private(i,j,es,qc) DO J=JSTA,JEND DO I=1,IM - IF (T1(I,J) < SPVAL .AND. P1(I,J) < SPVAL.AND.Q1(I,J)/=SPVAL) THEN + IF (T1(I,J) < spval .AND. P1(I,J) < spval.AND.Q1(I,J)/=spval) THEN ! IF (ABS(P1(I,J)) > 1.0) THEN ! IF (P1(I,J) > 1.0) THEN IF (P1(I,J) >= 1.0) THEN @@ -269,7 +269,7 @@ SUBROUTINE CALRH_GFS(P1,T1,Q1,RH) ENDIF ELSE - RH(I,J) = SPVAL + RH(I,J) = spval ENDIF ENDDO ENDDO @@ -295,7 +295,7 @@ SUBROUTINE CALRH_GSD(P1,T1,Q1,RHB) DO J=JSTA,JEND DO I=1,IM - IF (T1(I,J) < SPVAL .AND. P1(I,J) < SPVAL .AND. Q1(I,J) < SPVAL) THEN + IF (T1(I,J) < spval .AND. P1(I,J) < spval .AND. Q1(I,J) < spval) THEN ! - compute relative humidity Tx=T1(I,J)-273.15 POL = 0.99999683 + TX*(-0.90826951E-02 + & @@ -309,7 +309,7 @@ SUBROUTINE CALRH_GSD(P1,T1,Q1,RHB) E = P1(I,J)/100.*Q1(I,J)/(0.62197+Q1(I,J)*0.37803) RHB(I,J) = MIN(1.,E/ES) ELSE - RHB(I,J) = SPVAL + RHB(I,J) = spval ENDIF ENDDO ENDDO @@ -344,8 +344,8 @@ SUBROUTINE CALRH_PW(RHPW) k=lm-l+1 DO J=JSTA,JEND DO I=1,IM - if(t(i,j,k) THETAA, ADD TO THE CAPE SUM. +! (B) IF THETAP < THETAA, ADD TO THE CINS SUM. +! (7) ARE WE AT EQUILIBRIUM LEVEL? +! (A) IF YES, STOP THE SUMMATION. +! (B) IF NO, CONTIUNUE THE SUMMATION. +! (8) ENFORCE LIMITS ON CAPE AND CINS (I.E. NO NEGATIVE CAPE) +! +! PROGRAM HISTORY LOG: +! 93-02-10 RUSS TREADON +! 93-06-19 RUSS TREADON - GENERALIZED ROUTINE TO ALLOW FOR +! TYPE 2 CAPE/CINS CALCULATIONS. +! 94-09-23 MIKE BALDWIN - MODIFIED TO USE LOOK UP TABLES +! INSTEAD OF COMPLICATED EQUATIONS. +! 94-10-13 MIKE BALDWIN - MODIFIED TO CONTINUE CAPE/CINS CALC +! UP TO AT HIGHEST BUOYANT LAYER. +! 98-06-12 T BLACK - CONVERSION FROM 1-D TO 2-D +! 98-08-18 T BLACK - COMPUTE APE INTERNALLY +! 00-01-04 JIM TUCCILLO - MPI VERSION +! 02-01-15 MIKE BALDWIN - WRF VERSION +! 03-08-24 G MANIKIN - ADDED LEVEL OF PARCEL BEING LIFTED +! AS OUTPUT FROM THE ROUTINE AND ADDED +! THE DEPTH OVER WHICH ONE SEARCHES FOR +! THE MOST UNSTABLE PARCEL AS INPUT +! 10-09-09 G MANIKIN - CHANGED COMPUTATION TO USE VIRTUAL TEMP +! - ADDED EQ LVL HGHT AND THUNDER PARAMETER +! 15-xx-xx S MOORTHI - optimization and threading +! +! USAGE: CALL CALCAPE(ITYPE,DPBND,P1D,T1D,Q1D,L1D,CAPE, +! CINS,PPARC) +! INPUT ARGUMENT LIST: +! ITYPE - INTEGER FLAG SPECIFYING HOW PARCEL TO LIFT IS +! IDENTIFIED. SEE COMMENTS ABOVE. +! DPBND - DEPTH OVER WHICH ONE SEARCHES FOR MOST UNSTABLE PARCEL +! P1D - ARRAY OF PRESSURE OF PARCELS TO LIFT. +! T1D - ARRAY OF TEMPERATURE OF PARCELS TO LIFT. +! Q1D - ARRAY OF SPECIFIC HUMIDITY OF PARCELS TO LIFT. +! L1D - ARRAY OF MODEL LEVEL OF PARCELS TO LIFT. +! +! OUTPUT ARGUMENT LIST: +! CAPE - CONVECTIVE AVAILABLE POTENTIAL ENERGY (J/KG) +! CINS - CONVECTIVE INHIBITION (J/KG) +! PPARC - PRESSURE LEVEL OF PARCEL LIFTED WHEN ONE SEARCHES +! OVER A PARTICULAR DEPTH TO COMPUTE CAPE/CIN +! +! OUTPUT FILES: +! STDOUT - RUN TIME STANDARD OUT. +! +! SUBPROGRAMS CALLED: +! UTILITIES: +! BOUND - BOUND (CLIP) DATA BETWEEN UPPER AND LOWER LIMTS. +! TTBLEX - LOOKUP TABLE ROUTINE TO GET T FROM THETAE AND P +! +! LIBRARY: +! COMMON - +! +! ATTRIBUTES: +! LANGUAGE: FORTRAN 90 +! MACHINE : CRAY C-90 +!$$$ +! +! use vrbls3d, only: pmid, t, q, zint +! use vrbls2d, only: teql +! use masks, only: lmh + + use params_mod, only: d00, h1m12, h99999, h10e5, capa, elocp, eps, & + oneps, g + use lookup_mod, only: thl, rdth, jtb, qs0, sqs, rdq, itb, ptbl, & + plq, ttbl, pl, rdp, the0, sthe, rdthe, ttblq, & + itbq, jtbq, rdpq, the0q, stheq, rdtheq + use ctlblk_mod, only: jsta_2l, jend_2u, lm, jsta, jend, im, me, spval +! +!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + implicit none +! +! INCLUDE/SET PARAMETERS. CONSTANTS ARE FROM BOLTON (MWR, 1980). + real,PARAMETER :: ISMTHP=2,ISMTHT=2,ISMTHQ=2 +! +! DECLARE PASSED VARIABLES +! + real, intent(in) :: LPAR0 + real, intent(in) :: PSFC + real, intent(in) :: LMASK + real, dimension(1:LM), intent(in) :: P_A, T_A, Q_A, ZINT_A + + real, intent(out) :: CAPE, CINS + real, dimension(1:LM), intent(out) :: TPAR + real, dimension(4), intent(inout) :: PARCEL0 + integer, intent(out) :: LLCL, LEQL +! +! DECLARE LOCAL VARIABLES +! +! integer, dimension(im,jsta:jend) :: IPTB, ITHTB, PARCEL, KLRES, KHRES, LCL, IDX + integer :: IPTB, ITHTB, KLRES, KHRES +! +! real, dimension(im,jsta:jend) :: THESP, PSP, CAPE20, QQ, PP, THUND + real :: THESP, PSP, QQ, PP + + LOGICAL THUNDER(IM,jsta:jend), NEEDTHUN + real PSFCK,PKL,TBTK,QBTK,APEBTK,TTHBTK,TTHK,APESPK,TPSPK, & + BQS00K,SQS00K,BQS10K,SQS10K,BQK,SQK,TQK,PRESK,GDZKL,THETAP, & + THETAA,P00K,P10K,P01K,P11K,TTHESK,ESATP,QSATP,TVP,TV + real DPBND + integer LP0 +! real,external :: fpvsnew + integer I,J,L,KNUML,KNUMH,LBEG,LEND,IQ, KB,ITTBK + integer PTYPE, PARCEL_L + +! integer I,J,L,KNUML,KNUMH,LBEG,LEND,IQ,IT,LMHK, KB,ITTBK +! +!************************************************************** +! START CALCAPE1D HERE. +! +! COMPUTE CAPE/CINS +! +! WHICH IS: THE SUM FROM THE LCL TO THE EQ LEVEL OF +! G * (LN(THETAP) - LN(THETAA)) * DZ +! +! (POSITIVE AREA FOR CAPE, NEGATIVE FOR CINS) +! +! WHERE: +! THETAP IS THE PARCEL THETA +! THETAA IS THE AMBIENT THETA +! DZ IS THE THICKNESS OF THE LAYER +! +! USING LCL AS LEVEL DIRECTLY BELOW SATURATION POINT +! AND EQ LEVEL IS THE HIGHEST POSITIVELY BUOYANT LEVEL. +! +! IEQL = EQ LEVEL +! P_thetaemax - real pressure of theta-e max parcel (Pa) +! +! INITIALIZE CAPE AND CINS ARRAYS +! + CAPE = D00 + CINS = D00 + PARCEL0(1:4) = D00 + LLCL = 0 + LEQL = LM + THESP = D00 + PSP = D00 +! + DO L=1,LM +! TPAR(L) = D00 + TPAR(L) = -9999.0 + ENDDO +! +! BASED ON LP0, CHOOSE THE WAY TO SELECT THE AIR PARCEL: +! + PTYPE = 0 + DPBND = 300.00 !DEFAULT GIVEN DEPTH + IF (LPAR0 <= -100.0) THEN +!-------FOR PTYPE=1--FIND MAXIMUM THETA E LAYER IN LOWEST DPBND ABOVE GROUND------- + PTYPE = 1 + DPBND = ABS(LPAR0) !GIVEN DEPTH + LP0 = -1 !DUMMY SINCE PARCEL IS IN SEARCHING. + ELSE IF ( NINT(LPAR0) == -9 ) THEN !GIVEN PARCEL +!-------FOR PTYPE=2--AIR PARCEL IS GIVEN BY PARCEL0------- + PTYPE = 2 + LP0 = 1 + ELSE IF ( NINT(LPAR0) == -1 ) THEN !SURFACE BASED PARCEL + PTYPE = 0 + LP0 = LM + ELSE IF ( NINT(LPAR0) == 0 ) THEN !IF LPAR0 IS 0, SET IT TO MODEL TOP + PTYPE = 0 + LP0 = 1 + ELSE IF ( NINT(LPAR0) >= LM ) THEN !IF LPAR0 IS BEYOND MODEL BASE, SET IT TO BASE LEVEL + PTYPE = 0 + LP0 = LM + ELSE IF ( NINT(LPAR0) >=1 .AND. NINT(LPAR0) <=LM ) THEN !LPAR0 IS IN GOOD RANGE + PTYPE = 0 + LP0 = NINT(LPAR0) + ELSE !INVALID VALUE OF LPAR0 + WRITE(0,*) "CALCAPE1D: INVALID VALUE PASSED TO INITIAL LEVEL OF AIR PARCEL" + STOP 1 + END IF + PARCEL_L = LP0 !PARCEL_L TO STORE LEVEL OF FOUNDED PARCEL(PTYPE=1) +! +!-------FOR ITYPE=1--FIND MAXIMUM THETA E LAYER IN LOWEST DPBND ABOVE GROUND------- +!--------------TRIAL MAXIMUM BUOYANCY LEVEL VARIABLES------------------- +!-------FOR ITYPE=2--FIND THETA E LAYER OF GIVEN T1D, Q1D, P1D--------------------- + + PSFCK = PSFC + DO KB=1,LM + PKL = P_A(KB) + IF ( (PTYPE == 0 .AND. KB == LP0) .OR. & + (PTYPE == 2 .AND. KB == LP0) .OR. & + (PTYPE == 1 .AND. (PKL >= PSFCK-DPBND .AND. PKL <= PSFCK)) ) THEN + + IF (PTYPE == 2 ) THEN !USING INFO OF GIVEN PARCEL + PKL = PARCEL0(2) + TBTK = PARCEL0(3) + QBTK = MAX(0.0, MIN(MAX(H1M12,PARCEL0(4)),H99999)) + APEBTK = (H10E5/PKL)**CAPA + ELSE !PARCEL IS CHOSEN FROM AMBIENT ATMOS. + PKL = P_A(KB) + TBTK = T_A(KB) + QBTK = max(0.0, Q_A(KB)) + APEBTK = (H10E5/PKL)**CAPA + ENDIF + +!----------Breogan Gomez - 2009-02-06 +! To prevent QBTK to be less than 0 which leads to a unrealistic value of PRESK +! and a floating invalid. + +! if(QBTK < 0) QBTK = 0 + +!--------------SCALING POTENTIAL TEMPERATURE & TABLE INDEX-------------- + TTHBTK = TBTK*APEBTK + TTHK = (TTHBTK-THL)*RDTH + QQ = TTHK - AINT(TTHK) + ITTBK = INT(TTHK) + 1 +!--------------KEEPING INDICES WITHIN THE TABLE------------------------- + IF(ITTBK < 1) THEN + ITTBK = 1 + QQ = D00 + ENDIF + IF(ITTBK >= JTB) THEN + ITTBK = JTB-1 + QQ = D00 + ENDIF +!--------------BASE AND SCALING FACTOR FOR SPEC. HUMIDITY--------------- + BQS00K = QS0(ITTBK) + SQS00K = SQS(ITTBK) + BQS10K = QS0(ITTBK+1) + SQS10K = SQS(ITTBK+1) +!--------------SCALING SPEC. HUMIDITY & TABLE INDEX--------------------- + BQK = (BQS10K-BQS00K)*QQ + BQS00K + SQK = (SQS10K-SQS00K)*QQ + SQS00K + TQK = (QBTK-BQK)/SQK*RDQ + PP = TQK-AINT(TQK) + IQ = INT(TQK)+1 +!--------------KEEPING INDICES WITHIN THE TABLE------------------------- + IF(IQ < 1) THEN + IQ = 1 + PP = D00 + ENDIF + IF(IQ >= ITB) THEN + IQ = ITB-1 + PP = D00 + ENDIF +!--------------SATURATION PRESSURE AT FOUR SURROUNDING TABLE PTS.------- + P00K = PTBL(IQ ,ITTBK ) + P10K = PTBL(IQ+1,ITTBK ) + P01K = PTBL(IQ ,ITTBK+1) + P11K = PTBL(IQ+1,ITTBK+1) +!--------------SATURATION POINT VARIABLES AT THE BOTTOM----------------- + TPSPK = P00K + (P10K-P00K)*PP + (P01K-P00K)*QQ & + + (P00K-P10K-P01K+P11K)*PP*QQ + +!!from WPP::tgs APESPK=(H10E5/TPSPK)**CAPA + if (TPSPK > 1.0e-3) then + APESPK = (max(0.,H10E5/ TPSPK))**CAPA + else + APESPK = 0.0 + endif + + TTHESK = TTHBTK * EXP(ELOCP*QBTK*APESPK/TTHBTK) +!--------------CHECK FOR MAXIMUM THETA E-------------------------------- + IF(TTHESK > THESP) THEN + PSP = TPSPK + THESP = TTHESK + PARCEL_L = KB ! MARK THE LEVEL OF MAX THETA-E + ENDIF + END IF ! KB == LP0 or PTYPE==1 (DPBND) + ENDDO ! KB loop + +!----FIND THE PRESSURE OF THE PARCEL THAT WAS LIFTED +!-----SAVE THE INITIAL PARCEL AND OUTPUT + IF (PTYPE == 0 .OR. PTYPE == 1) THEN + PARCEL0(1) = FLOAT(PARCEL_L) + PARCEL0(2) = P_A(PARCEL_L) + PARCEL0(3) = T_A(PARCEL_L) + PARCEL0(4) = Q_A(PARCEL_L) + ELSE IF (PTYPE == 1) THEN + PARCEL0(1) = FLOAT(LP0) + PARCEL0(2) = P_A(LP0) + PARCEL0(3) = T_A(LP0) + PARCEL0(4) = Q_A(LP0) + END IF +! +!-----CHOOSE LAYER DIRECTLY BELOW PSP AS LCL AND------------------------ +!-----ENSURE THAT THE LCL IS ABOVE GROUND.------------------------------ +!-------(IN SOME RARE CASES FOR ITYPE=2, IT IS NOT)--------------------- +!-----Note: L=1,LM --> from top to bottom + DO L=1,LM + IF (P_A(L) < PSP) LLCL = L+1 + ENDDO + IF (LLCL > NINT(LMASK)) LLCL = NINT(LMASK) +!----------------------------------------------------------------------- +!---------FIND TEMP OF PARCEL LIFTED ALONG MOIST ADIABAT (TPAR)--------- +!----------------------------------------------------------------------- +!-----Note: loop L=LM,1,-1 --> from bottom to top + DO L=LM,1,-1 + KLRES = 0 + KHRES = 0 +!--------------SCALING PRESSURE & TT TABLE INDEX------------------------ + IF(L <= LLCL) THEN + IF(P_A(L) < PLQ)THEN + KLRES = 1 + ELSE + KHRES = 1 + ENDIF + ENDIF +!*** +!*** COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE 0) THEN + CALL TTBLEX1D(TPAR(L),TTBL,ITB,JTB,KLRES,P_A(L), & + PL,QQ,PP,RDP,THE0,STHE,RDTHE,THESP, & + IPTB,ITHTB) + ENDIF +!*** +!*** COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE>PLQ +!** + IF(KHRES > 0) THEN + CALL TTBLEX1D(TPAR(L),TTBLQ,ITBQ,JTBQ,KHRES,P_A(L), & + PLQ,QQ,PP,RDPQ,THE0Q,STHEQ,RDTHEQ,THESP, & + IPTB,ITHTB) + ENDIF + +!------------SEARCH FOR EQ LEVEL---------------------------------------- + IF(KHRES > 0) THEN + IF(TPAR(L) > T_A(L)) LEQL = L + ENDIF +! + IF(KLRES > 0) THEN + IF(TPAR(L) > T_A(L) .AND. & + P_A(L)>100.) LEQL = L + ENDIF +!----------------------------------------------------------------------- + ENDDO ! end of do l=lm,1,-1 loop +!------------COMPUTE CAPE AND CINS-------------------------------------- + LBEG = 1000 + LEND = 0 + LBEG = MIN(LEQL,LBEG) + LEND = MAX(LLCL,LEND) +! +!-----Note: LOOP for CAPE: L=LEQL, LLCL --> FROM HIGH TO LOW + DO L=LBEG,LEND + + IF(L >= LEQL .AND. L <= LLCL) THEN + PRESK = P_A(L) + GDZKL = (ZINT_A(L)-ZINT_A(L+1)) * G + ESATP = min(FPVSNEW(TPAR(L)),PRESK) + QSATP = EPS*ESATP/(PRESK-ESATP*ONEPS) +! TVP = TPAR(I,J,L)*(1+0.608*QSATP) + TVP = TVIRTUAL(TPAR(L),QSATP) + THETAP = TVP*(H10E5/PRESK)**CAPA +! TV = T(I,J,L)*(1+0.608*Q(I,J,L)) + TV = TVIRTUAL(T_A(L),Q_A(L)) + THETAA = TV*(H10E5/PRESK)**CAPA + IF(THETAP < THETAA) THEN + CINS = CINS + (LOG(THETAP)-LOG(THETAA))*GDZKL + ELSEIF(THETAP > THETAA) THEN + CAPE = CAPE + (LOG(THETAP)-LOG(THETAA))*GDZKL + ENDIF + ENDIF + + ENDDO +! +! ENFORCE LOWER LIMIT OF 0.0 ON CAPE AND UPPER +! LIMIT OF 0.0 ON CINS. +! + CAPE = MAX(D00,CAPE) + CINS = MIN(CINS,D00) +! + RETURN +! + END SUBROUTINE CALCAPE1D +! +!------------------------------------------------------------------------------------- ! elemental function TVIRTUAL(T,Q) ! diff --git a/sorc/ncep_post.fd/VRBLS2D_mod.f b/sorc/ncep_post.fd/VRBLS2D_mod.f index 1067d1432..11c542b1d 100644 --- a/sorc/ncep_post.fd/VRBLS2D_mod.f +++ b/sorc/ncep_post.fd/VRBLS2D_mod.f @@ -28,7 +28,7 @@ module vrbls2d ,PBLH(:,:),PBLHGUST(:,:),HBOTD(:,:),HTOPD(:,:),HBOTS(:,:),HTOPS(:,:) & ,CLDEFI(:,:),ALBASE(:,:),SI(:,:),LSPA(:,:) & ,RSWINC(:,:),VIS(:,:),PD(:,:),MXSNAL(:,:),MIXHT(:,:) & - ,SNONC(:,:),EPSR(:,:),RSWTOA(:,:),TEQL(:,:) & + ,SNONC(:,:),EPSR(:,:),RSWTOA(:,:),TEQL(:,:) & ! HWRF additions ,MDLTAUX(:,:),MDLTAUY(:,:),CD10(:,:),CH10(:,:) & ,ACSWUPT(:,:),SWDNT(:,:),ACSWDNT(:,:) & @@ -82,7 +82,8 @@ module vrbls2d ,avgesnow(:,:),avgpotevp(:,:),avgprec_cont(:,:),avgcprate_cont(:,:)& ,ti(:,:),aod550(:,:),du_aod550(:,:),ss_aod550(:,:),su_aod550(:,:) & ,bc_aod550(:,:),oc_aod550(:,:) - integer, allocatable :: IVGTYP(:,:),ISLTYP(:,:),ISLOPE(:,:) + integer, allocatable :: IVGTYP(:,:),ISLTYP(:,:),ISLOPE(:,:) & + ,IEQL(:,:) ! Add 2d aerosol diagnosis fields for GOCART (NGAC) real, allocatable :: & DUSMASS(:,:),DUCMASS(:,:),DUSMASS25(:,:),DUCMASS25(:,:) & From de51871700a10b131455ddd4f7b6f164818613da Mon Sep 17 00:00:00 2001 From: "Edward.Colon" Date: Thu, 5 Aug 2021 23:43:29 +0000 Subject: [PATCH 04/14] This commit includes changes to the computation of the effective layer top and bottom. --- sorc/ncep_post.fd/MISCLN.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sorc/ncep_post.fd/MISCLN.f b/sorc/ncep_post.fd/MISCLN.f index 4ecbf4569..61caf5967 100644 --- a/sorc/ncep_post.fd/MISCLN.f +++ b/sorc/ncep_post.fd/MISCLN.f @@ -3726,7 +3726,7 @@ SUBROUTINE MISCLN ENDIF ! IF PASSING SANITY CHECK - IF ( FOUND_BASE(I,J) /= FOUND_TOPS(I,J) ) THEN + IF ( FOUND_BASE(I,J) .neqv. FOUND_TOPS(I,J) ) THEN WRITE(0,'(1x,A,A,A,I6,1x,I6)') "BASE & TOP OF ", & " EFFECTIVE LAYER ARE NOT FOUND TOGETHER. WRONG! ", & " ABORT! ABORT! ... AT GRID POINT: ", I, J From 5b399d4805a43f6701a449d43fa32411ddb36d4e Mon Sep 17 00:00:00 2001 From: EdwardColon-NOAA Date: Mon, 16 Aug 2021 11:16:15 -0400 Subject: [PATCH 05/14] Update ALLOCATE_ALL.f i loop indices corrected. --- sorc/ncep_post.fd/ALLOCATE_ALL.f | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/sorc/ncep_post.fd/ALLOCATE_ALL.f b/sorc/ncep_post.fd/ALLOCATE_ALL.f index 57c6b8033..7b84283f9 100644 --- a/sorc/ncep_post.fd/ALLOCATE_ALL.f +++ b/sorc/ncep_post.fd/ALLOCATE_ALL.f @@ -79,7 +79,7 @@ SUBROUTINE ALLOCATE_ALL() !$omp parallel do private(i,j,l) do l=1,lm do j=jsta_2l,jend_2u - do i=1,lm + do i=1,im u(i,j,l)=0. v(i,j,l)=0. t(i,j,l)=spval @@ -107,7 +107,7 @@ SUBROUTINE ALLOCATE_ALL() !$omp parallel do private(i,j,l) do l=1,lp1 do j=jsta_2l,jend_2u - do i=1,lm + do i=1,im pint(i,j,l)=spval alpint(i,j,l)=spval zint(i,j,l)=spval @@ -147,7 +147,7 @@ SUBROUTINE ALLOCATE_ALL() !$omp parallel do private(i,j,l) do l=1,lm do j=jsta_2l,jend_2u - do i=1,lm + do i=1,im cwm(i,j,l)=spval F_ice(i,j,l)=spval F_rain(i,j,l)=spval @@ -195,7 +195,7 @@ SUBROUTINE ALLOCATE_ALL() !$omp parallel do private(i,j,l) do l=1,lm do j=jsta_2l,jend_2u - do i=1,lm + do i=1,im NRAIN(i,j,l)=spval radius_cloud(i,j,l)=spval radius_ice(i,j,l)=spval From 85519598043dfebf2178707f2cb943b31591b13d Mon Sep 17 00:00:00 2001 From: EdwardColon-NOAA Date: Mon, 16 Aug 2021 11:32:37 -0400 Subject: [PATCH 06/14] Update MISCLN.f --- sorc/ncep_post.fd/MISCLN.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sorc/ncep_post.fd/MISCLN.f b/sorc/ncep_post.fd/MISCLN.f index 61caf5967..3b089d36a 100644 --- a/sorc/ncep_post.fd/MISCLN.f +++ b/sorc/ncep_post.fd/MISCLN.f @@ -198,7 +198,7 @@ SUBROUTINE MISCLN !**************************************************************************** ! START MISCLN HERE. ! - debugprint = .TRUE. + debugprint = .FALSE. EL_SCHEME = 1 allocate(USHR1(IM,jsta_2l:jend_2u),VSHR1(IM,jsta_2l:jend_2u), & From 662bfedfa6bafc3e2e50683b238adb0f19d0d8b8 Mon Sep 17 00:00:00 2001 From: EdwardColon-NOAA Date: Mon, 16 Aug 2021 11:53:34 -0400 Subject: [PATCH 07/14] Update MISCLN.f Removed EL_SCHEME flag and the EL_SCHEME 2 as it was intended as a method of validating the output of EL_SCHEME effective layer top and bottom calculations. The HWT/SFE confirmed that the EL_SCHEME 1 computations were viable making EL_SCHEME 2 superfluous. The debugprint flag was also set to false as it currently not needed in the execution of the MISCLN subprogram, however, it may be switch on at any point it becomes necessary to trace a problem. --- sorc/ncep_post.fd/MISCLN.f | 270 ++----------------------------------- 1 file changed, 12 insertions(+), 258 deletions(-) diff --git a/sorc/ncep_post.fd/MISCLN.f b/sorc/ncep_post.fd/MISCLN.f index 3b089d36a..cf02fd640 100644 --- a/sorc/ncep_post.fd/MISCLN.f +++ b/sorc/ncep_post.fd/MISCLN.f @@ -184,7 +184,6 @@ SUBROUTINE MISCLN INTEGER :: IREC, IUNIT INTEGER :: IREC2, IUNIT2 LOGICAL :: debugprint - INTEGER :: EL_SCHEME INTEGER :: LLL INTEGER :: LLCL_PAR, LEQL_PAR REAL :: LMASK, PSFC, CAPE_PAR, CINS_PAR, LPAR0 @@ -199,7 +198,7 @@ SUBROUTINE MISCLN ! START MISCLN HERE. ! debugprint = .FALSE. - EL_SCHEME = 1 + allocate(USHR1(IM,jsta_2l:jend_2u),VSHR1(IM,jsta_2l:jend_2u), & USHR6(IM,jsta_2l:jend_2u),VSHR6(IM,jsta_2l:jend_2u)) @@ -3461,7 +3460,7 @@ SUBROUTINE MISCLN ! ! --- Effective (inflow) Layer (EL) ! - IF ( EL_SCHEME > 0 ) THEN + ALLOCATE(EL_BASE(IM,JSTA_2L:JEND_2U)) ALLOCATE(EL_TOPS(IM,JSTA_2L:JEND_2U)) ALLOCATE(FOUND_BASE(IM,JSTA_2L:JEND_2U)) @@ -3475,9 +3474,9 @@ SUBROUTINE MISCLN FOUND_TOPS(I,J) = .FALSE. ENDDO ENDDO - END IF + ! - IF ( EL_SCHEME == 1 ) THEN + ITYPE = 2 DPBND = 0. @@ -3497,8 +3496,7 @@ SUBROUTINE MISCLN ENDDO !--- CALCULATE CAPE/CIN FOR ALL AIR PARCELS on LEVEL L - IF (debugprint) WRITE(1000+ME,'(1x,A,I2.2,2x,A,I3)') & - 'NEW EL_SCHEME:', EL_SCHEME, & + IF (debugprint) WRITE(1000+ME,'(A,I3)') & ' CALCULATING CAPE/CINS ON LEVEL:',L CALL CALCAPE(ITYPE,DPBND,P1D,T1D,Q1D,IDUMMY,EGRID1, & EGRID2,EGRID3,EGRID4,EGRID5) @@ -3530,229 +3528,12 @@ SUBROUTINE MISCLN ENDDO END DO ! L -! - ELSE IF ( EL_SCHEME == 2 ) THEN - -!--- SEARCH FOR EL ALONG EACH PROFILE/COLUMN - ALLOCATE(L_THETAE_MAX(IM,JSTA:JEND)) -!$omp parallel do private(i,j) - DO J=JSTA,JEND - DO I=1,IM - EGRID1(I,J) = -H99999 - EGRID2(I,J) = -H99999 - L_THETAE_MAX(I,J) = -9999 - ENDDO - ENDDO - -!------ SEARCH FOR PARCEL WITH MAX THETA-E OF EACH PROFILE/COLUMN - DO L = LM,1,-1 -!$omp parallel do private(i,j) - DO J=JSTA,JEND - DO I=1,IM - EGRID1(I,J) = -H99999 - P1D(I,J) = PMID(I,J,L) - T1D(I,J) = T(I,J,L) - Q1D(I,J) = Q(I,J,L) - ENDDO - ENDDO -! CALL CALTHTE(PMID(1,jsta,L),T(1,jsta,L),Q(1,jsta,L),EGRID1) - CALL CALTHTE(P1D,T1D,Q1D,EGRID1) -!$omp parallel do private(i,j) - DO J=JSTA,JEND - DO I=1,IM - IF (EGRID1(I,J) > EGRID2(I,J)) THEN - EGRID2(I,J) = EGRID1(I,J) - L_THETAE_MAX(I,J) = L - END IF - END DO - END DO - - END DO ! L - -!--- SET UP THE PARCELS WIH MAX THETA-E - ALLOCATE(CAPE9(IM,JSTA:JEND)) - ALLOCATE(CINS9(IM,JSTA:JEND)) -!$omp parallel do private(i,j,LLL) - DO J=JSTA,JEND - DO I=1,IM - LLL = L_THETAE_MAX(I,J) - P1D(I,J) = PMID(I,J,LLL) - T1D(I,J) = T(I,J,LLL) - Q1D(I,J) = Q(I,J,LLL) - ENDDO - ENDDO -!--- COMPUTE CAPE/CIN OF PARCEL WITH MAX THETA-E - ITYPE = 2 - DPBND = 0. -!$omp parallel do private(i,j) - DO J=JSTA,JEND - DO I=1,IM - EGRID1(I,J) = -H99999 - EGRID2(I,J) = -H99999 - CAPE9(I,J) = -H99999 - CINS9(I,J) = -H99999 - IDUMMY(I,J) = 0 - ENDDO - ENDDO - CALL CALCAPE(ITYPE,DPBND,P1D,T1D,Q1D,IDUMMY,EGRID1, & - EGRID2,EGRID3,EGRID4,EGRID5) -!--- SANITY CHECK IF CAPE/CIN OF MAX THETA-E PARCEL OF THIS COLUMN -!$omp parallel do private(i,j,LLL) - DO J=JSTA,JEND - DO I=1,IM - CAPE9(I,J) = EGRID1(I,J) - CINS9(I,J) = EGRID2(I,j) - ENDDO - ENDDO - - ALLOCATE(TPAR_B(LM), TPAR_T(LM)) - ALLOCATE(TPAR_TMP(LM)) - ALLOCATE(P_AMB(LM), T_AMB(LM), Q_AMB(LM), ZINT_AMB(LM)) - IF (debugprint) THEN - ALLOCATE(TPAR_BASE(IM,JSTA:JEND,LM)) - ALLOCATE(TPAR_TOPS(IM,JSTA:JEND,LM)) -!$omp parallel do private(i,j,l) - DO L=1,LM - DO J=JSTA,JEND - DO I=1,IM - TPAR_BASE(I,J,L) = DM9999 - TPAR_TOPS(I,J,L) = DM9999 - ENDDO - ENDDO - ENDDO - END IF - -!$OMP PARALLEL DO PRIVATE(i,j,l,P_AMB,T_AMB,Q_AMB,ZINT_AMB,LMASK,PSFC, & -!$OMP& TPAR_B,TPAR_T,LPAR0,CAPE_PAR,CINS_PAR, & -!$OMP& PARCEL0,LLCL_PAR,LEQL_PAR,TPAR_TMP) - DO J=JSTA,JEND - DO I=1,IM - -!-----------GET THE AMBIENT PROFILE - DO L=1,LM - P_AMB(L) = PMID(I,J,L) - T_AMB(L) = T(I,J,L) - Q_AMB(L) = Q(I,J,L) - ZINT_AMB(L) = ZINT(I,J,L) - END DO - LMASK = NINT(LMH(I,J)) - PSFC = PMID(I,J,NINT(LMH(I,J))) - EL_BASE(I,J)=LM - EL_TOPS(I,J)=LM - FOUND_BASE(I,J) = .FALSE. - FOUND_TOPS(I,J) = .FALSE. - TPAR_B(1:LM) = DM9999 - TPAR_T(1:LM) = DM9999 - TPAR_TMP(1:LM) = DM9999 - -!---------- SANITY CHECK -!---------- USING CAPE/CINS OF AIR PARCEL WITH MAX THETA-E FOR SANITY CHECK FIRST -!---------- CAPE4/CINS4 ARE CAPE/CINS CALCULATED FROM CODE ABOVE FOR AIR PARCEL -!---------- WITH MAX THETA-E. - IF (CAPE9(I,J) >= 100. .OR. CINS9(I,J) >= -250.) THEN !EL FOR THIS COLUMN -! EL SHOULD EXIST ALONG THIS COLUMN/PROFILE. -!-------------NEED TO SEARCH FROM BOTTOM TO TOP, COMPUTE CAPE/CINS OF AIR PARCEL -! INITIATING AT EACH LEVEL UNTIL THE CRITERIA MEETS. - - VLOOP: DO L=LM,1,-1 -!---------------CALCULATE CAPE/CIN OF AIR PARCEL INITIALIZED AT LEVEL L - IF (.NOT. FOUND_BASE(I,J)) THEN !SEARCH FOR BASE FIRST - LPAR0 = FLOAT(L) - CAPE_PAR = D00 - CINS_PAR = D00 - PARCEL0(1) = LPAR0 - PARCEL0(2:4) = D00 - LLCL_PAR = 1 - LEQL_PAR = LM - TPAR_B(1:LM) = DM9999 - TPAR_T(1:LM) = DM9999 - - WRITE(1000+me,*)"CALCAPE2_1D_base: parcel@I J L:",I,J,L - - CALL CALCAPE1D(P_AMB,T_AMB,Q_AMB,ZINT_AMB,LPAR0, & - PSFC,LMASK,CAPE_PAR,CINS_PAR,TPAR_B, & - PARCEL0,LLCL_PAR,LEQL_PAR) - IF (CAPE_PAR >= 100. .AND. CINS_PAR >= -250.) THEN ! BASE OF EL - EL_BASE(I,J) = L - FOUND_BASE(I,J) = .TRUE. - ELSE - TPAR_B(1:LM) = DM9999 ! RESET TPAR IF NOT FOUND BASE - ENDIF - ELSE !SEARCH FOR TOP AFTER BASE IS FOUND - IF (.NOT. FOUND_TOPS(I,J)) THEN !SEARCH FOR BASE FIRST - LPAR0 = FLOAT(L) - CAPE_PAR = D00 - CINS_PAR = D00 - TPAR_T(1:LM) = DM9999 - PARCEL0(1) = LPAR0 - PARCEL0(2:4) = D00 - LLCL_PAR = 1 - LEQL_PAR = LM - - WRITE(1000+me,*)"CALCAPE2_1D_top: parcel@I J L:",I,J,L - - CALL CALCAPE1D(P_AMB,T_AMB,Q_AMB,ZINT_AMB,LPAR0, & - PSFC,LMASK,CAPE_PAR,CINS_PAR,TPAR_T, & - PARCEL0,LLCL_PAR,LEQL_PAR) - IF (CAPE_PAR < 100. .OR. CINS_PAR < -250.) THEN ! TOP OF EL - EL_TOPS(I,J) = L+1 - FOUND_TOPS(I,J) = .TRUE. - IF (EL_TOPS(I,J) == EL_BASE(I,J)) THEN - TPAR_T(1:LM) = TPAR_B(1:LM) - ELSE IF (EL_TOPS(I,J) < EL_BASE(I,J) ) THEN - TPAR_T(1:LM) = TPAR_TMP(1:LM) - ELSE - WRITE(0,'(1x,A,A)') "TOP OF EFFECTIVE LAYER IS", & - " LOWER THAN BASE. WRONG! ABORT ..." - STOP 9 - END IF - EXIT VLOOP - ELSE - TPAR_TMP(1:LM) = TPAR_T(1:LM) - TPAR_T(1:LM) = DM9999 ! RESET TPAR IF NOT FOUND TOP - ENDIF - ENDIF ! FOUND_TOPS - - ENDIF ! FOUND_BASE OR NOT - - ENDDO VLOOP - - IF (debugprint) THEN - DO L=1,LM - TPAR_BASE(I,J,L) = TPAR_B(L) - TPAR_TOPS(I,J,L) = TPAR_T(L) - END DO - END IF - - ENDIF ! IF PASSING SANITY CHECK - - IF ( FOUND_BASE(I,J) .neqv. FOUND_TOPS(I,J) ) THEN - WRITE(0,'(1x,A,A,A,I6,1x,I6)') "BASE & TOP OF ", & - " EFFECTIVE LAYER ARE NOT FOUND TOGETHER. WRONG! ", & - " ABORT! ABORT! ... AT GRID POINT: ", I, J - STOP 10 - END IF - - ENDDO ! J - ENDDO ! I - - IF(ALLOCATED(L_THETAE_MAX)) DEALLOCATE(L_THETAE_MAX) - IF(ALLOCATED(CAPE9)) DEALLOCATE(CAPE9) - IF(ALLOCATED(CINS9)) DEALLOCATE(CINS9) - IF(ALLOCATED(TPAR_B)) DEALLOCATE(TPAR_B) - IF(ALLOCATED(TPAR_T)) DEALLOCATE(TPAR_T) - IF(ALLOCATED(TPAR_TMP)) DEALLOCATE(TPAR_TMP) - IF(ALLOCATED(P_AMB)) DEALLOCATE(P_AMB) - IF(ALLOCATED(T_AMB)) DEALLOCATE(T_AMB) - IF(ALLOCATED(Q_AMB)) DEALLOCATE(Q_AMB) - IF(ALLOCATED(ZINT_AMB)) DEALLOCATE(ZINT_AMB) - - END IF ! EL_SCHEME :1 OR 2 + IF (ALLOCATED(FOUND_BASE)) DEALLOCATE(FOUND_BASE) IF (ALLOCATED(FOUND_TOPS)) DEALLOCATE(FOUND_TOPS) - IF (debugprint .AND. EL_SCHEME > 0) THEN + IF (debugprint) THEN WRITE(IM_CH,'(I5.5)') IM WRITE(JSTA_CH,'(I5.5)') JSTA WRITE(JEND_CH,'(I5.5)') JEND @@ -3765,44 +3546,18 @@ SUBROUTINE MISCLN IREC=0 IREC2=0 OPEN(IUNIT,FILE=TRIM(ADJUSTL(EFFL_FNAME)),FORM='FORMATTED') - IF (EL_SCHEME == 2) THEN - OPEN(IUNIT2,FILE=TRIM(ADJUSTL(EFFL_FNAME2)),FORM='FORMATTED') - END IF -! OPEN(IUNIT,FILE=TRIM(ADJUSTL(EFFL_FNAME)),FORM='UNFORMATTED', & -! ACCESS='DIRECT',RECL=4*6) + + DO J=JSTA,JEND DO I=1,IM IREC = IREC + 1 - IREC2 = IREC2 + 1 -! WRITE(IUNIT,'(1x,I6,2x,I6,2x,I6,2x,I6)')I,J,EL_BASE(I,J),EL_TOPS(I,J) - WRITE(IUNIT,'(1x,I6,2x,I6,2(2x,I6,2x,F12.3))') I, J, & + IREC2 = IREC2 + 1 + WRITE(IUNIT,'(1x,I6,2x,I6,2(2x,I6,2x,F12.3))') I, J, & EL_BASE(I,J),PMID(I,J,EL_BASE(I,J)), & EL_TOPS(I,J),PMID(I,J,EL_TOPS(I,J)) -! WRITE(IUNIT,REC=IREC) I, J, & -! EL_BASE(I,J),PMID(I,J,EL_BASE(I,J)), & -! EL_TOPS(I,J),PMID(I,J,EL_TOPS(I,J)) - IF (EL_SCHEME == 2) THEN - WRITE(IUNIT2,'(1x,I6,2x,I6,2(2x,I6,2x,F12.3))') I, J, & - EL_BASE(I,J),PMID(I,J,EL_BASE(I,J)), & - EL_TOPS(I,J),PMID(I,J,EL_TOPS(I,J)) - DO L=LM,1,-1 - IREC2=IREC2+1 - WRITE(IUNIT2,'(1x,I4,5(2x,F12.3))') & - LM+1-L, PMID(I,J,L), ZINT(I,J,L), T(I,J,L), & - TPAR_BASE(I,J,L), TPAR_TOPS(I,J,L) -! WRITE(IUNIT2,REC=IREC2) & -! LM+1-L, PMID(I,J,L), ZINT(I,J,L), T(I,J,L), & -! TPAR_BASE(I,J,L), TPAR_TOPS(I,J,L) - END DO - END IF END DO ENDDO CLOSE(IUNIT) - - IF (EL_SCHEME == 2) THEN - CLOSE(IUNIT2) - END IF - ENDIF IF(ALLOCATED(TPAR_BASE)) DEALLOCATE(TPAR_BASE) @@ -3999,7 +3754,6 @@ SUBROUTINE MISCLN ENDIF !--- IF USSING EL BASE & TOP COMPUTED BY NEW SCHEME FOR THE RELATED VARIABLES - IF ( EL_SCHEME > 0 ) THEN !$omp parallel do private(i,j) DO J=JSTA,JEND DO I=1,IM @@ -4007,7 +3761,7 @@ SUBROUTINE MISCLN LUPP(I,J) = EL_TOPS(I,J) ENDDO ENDDO - END IF + ! CALL CALHEL(DEPTH,UST,VST,HELI,USHR1,VSHR1,USHR6,VSHR6) CALL CALHEL2(LLOW,LUPP,DEPTH,UST,VST,HELI,CANGLE) From d8788f10204b4147228667b25abcb3f486df995d Mon Sep 17 00:00:00 2001 From: EdwardColon-NOAA Date: Mon, 16 Aug 2021 12:13:54 -0400 Subject: [PATCH 08/14] Update UPP_PHYSICS.f Removed the definition of CALCAPE1D which was invoked in MISCLN.f as the CAPE/CIN computational routine used in EL_SCHEME 2. It is now a defunct option. --- sorc/ncep_post.fd/UPP_PHYSICS.f | 433 +------------------------------- 1 file changed, 1 insertion(+), 432 deletions(-) diff --git a/sorc/ncep_post.fd/UPP_PHYSICS.f b/sorc/ncep_post.fd/UPP_PHYSICS.f index 684de251a..4487d6364 100644 --- a/sorc/ncep_post.fd/UPP_PHYSICS.f +++ b/sorc/ncep_post.fd/UPP_PHYSICS.f @@ -38,7 +38,7 @@ module upp_physics private - public :: CALCAPE, CALCAPE2, CALCAPE1D + public :: CALCAPE, CALCAPE2 public :: CALRH public :: CALRH_GFS, CALRH_GSD, CALRH_NAM public :: CALRH_PW @@ -1765,437 +1765,6 @@ END SUBROUTINE CALCAPE2 ! !------------------------------------------------------------------------------------- ! -! - - SUBROUTINE CALCAPE1D(P_A,T_A,Q_A,ZINT_A,LPAR0,PSFC,LMASK,CAPE,CINS, & - TPAR,PARCEL0,LLCL,LEQL) -!$$$ SUBPROGRAM DOCUMENTATION BLOCK -! . . . -! SUBPROGRAM: CALCAPE1D COMPUTES CAPE AND CINS ALONG 1-D VERTICAL PROFILE -! PRGRMMR: ZHAO ORG: IMSG@NCEP/EMC DATE: 21-04-xx -! BASED ON SUBROUTINE CALCAPE BY TREADON -! -! ABSTRACT: -! -! THIS ROUTINE COMPUTES CAPE AND CINS GIVEN TEMPERATURE, -! PRESSURE, AND SPECIFIC HUMIDTY. IN "STORM AND CLOUD -! DYNAMICS" (1989, ACADEMIC PRESS) COTTON AND ANTHES DEFINE -! CAPE (EQUATION 9.16, P501) AS -! -! EL -! CAPE = SUM G * LN(THETAP/THETAA) DZ -! LCL -! -! WHERE, -! EL = EQUILIBRIUM LEVEL, -! LCL = LIFTING CONDENSTATION LEVEL, -! G = GRAVITATIONAL ACCELERATION, -! THETAP = LIFTED PARCEL POTENTIAL TEMPERATURE, -! THETAA = AMBIENT POTENTIAL TEMPERATURE. -! -! NOTE THAT THE INTEGRAND LN(THETAP/THETAA) APPROXIMATELY -! EQUALS (THETAP-THETAA)/THETAA. THIS RATIO IS OFTEN USED -! IN THE DEFINITION OF CAPE/CINS. -! -! TWO TYPES OF CAPE/CINS CAN BE COMPUTED BY THIS ROUTINE. THE -! SUMMATION PROCESS IS THE SAME FOR BOTH CASES. WHAT DIFFERS -! IS THE DEFINITION OF THE PARCEL TO LIFT. FOR ITYPE=1 THE -! PARCEL WITH THE WARMEST THETA-E IN A DPBND PASCAL LAYER ABOVE -! THE MODEL SURFACE IS LIFTED. THE ARRAYS P1D, T1D, AND Q1D -! ARE NOT USED. FOR ITYPE=2 THE ARRAYS P1D, T1D, AND Q1D -! DEFINE THE PARCEL TO LIFT IN EACH COLUMN. BOTH TYPES OF -! CAPE/CINS MAY BE COMPUTED IN A SINGLE EXECUTION OF THE POST -! PROCESSOR. -! -! THIS ALGORITHM PROCEEDS AS FOLLOWS. -! FOR EACH COLUMN, -! (1) INITIALIZE RUNNING CAPE AND CINS SUM TO 0.0 -! (2) COMPUTE TEMPERATURE AND PRESSURE AT THE LCL USING -! LOOK UP TABLE (PTBL). USE EITHER PARCEL THAT GIVES -! MAX THETAE IN LOWEST DPBND ABOVE GROUND (ITYPE=1) -! OR GIVEN PARCEL FROM T1D,Q1D,...(ITYPE=2). -! (3) COMPUTE THE TEMP OF A PARCEL LIFTED FROM THE LCL. -! WE KNOW THAT THE PARCEL'S -! EQUIVALENT POTENTIAL TEMPERATURE (THESP) REMAINS -! CONSTANT THROUGH THIS PROCESS. WE CAN -! COMPUTE TPAR USING THIS KNOWLEDGE USING LOOK -! UP TABLE (SUBROUTINE TTBLEX). -! (4) FIND THE EQUILIBRIUM LEVEL. THIS IS DEFINED AS THE -! HIGHEST POSITIVELY BUOYANT LAYER. -! (IF THERE IS NO POSITIVELY BUOYANT LAYER, CAPE/CINS -! WILL BE ZERO) -! (5) COMPUTE CAPE/CINS. -! (A) COMPUTE THETAP. WE KNOW TPAR AND P. -! (B) COMPUTE THETAA. WE KNOW T AND P. -! (6) ADD G*(THETAP-THETAA)*DZ TO THE RUNNING CAPE OR CINS SUM. -! (A) IF THETAP > THETAA, ADD TO THE CAPE SUM. -! (B) IF THETAP < THETAA, ADD TO THE CINS SUM. -! (7) ARE WE AT EQUILIBRIUM LEVEL? -! (A) IF YES, STOP THE SUMMATION. -! (B) IF NO, CONTIUNUE THE SUMMATION. -! (8) ENFORCE LIMITS ON CAPE AND CINS (I.E. NO NEGATIVE CAPE) -! -! PROGRAM HISTORY LOG: -! 93-02-10 RUSS TREADON -! 93-06-19 RUSS TREADON - GENERALIZED ROUTINE TO ALLOW FOR -! TYPE 2 CAPE/CINS CALCULATIONS. -! 94-09-23 MIKE BALDWIN - MODIFIED TO USE LOOK UP TABLES -! INSTEAD OF COMPLICATED EQUATIONS. -! 94-10-13 MIKE BALDWIN - MODIFIED TO CONTINUE CAPE/CINS CALC -! UP TO AT HIGHEST BUOYANT LAYER. -! 98-06-12 T BLACK - CONVERSION FROM 1-D TO 2-D -! 98-08-18 T BLACK - COMPUTE APE INTERNALLY -! 00-01-04 JIM TUCCILLO - MPI VERSION -! 02-01-15 MIKE BALDWIN - WRF VERSION -! 03-08-24 G MANIKIN - ADDED LEVEL OF PARCEL BEING LIFTED -! AS OUTPUT FROM THE ROUTINE AND ADDED -! THE DEPTH OVER WHICH ONE SEARCHES FOR -! THE MOST UNSTABLE PARCEL AS INPUT -! 10-09-09 G MANIKIN - CHANGED COMPUTATION TO USE VIRTUAL TEMP -! - ADDED EQ LVL HGHT AND THUNDER PARAMETER -! 15-xx-xx S MOORTHI - optimization and threading -! -! USAGE: CALL CALCAPE(ITYPE,DPBND,P1D,T1D,Q1D,L1D,CAPE, -! CINS,PPARC) -! INPUT ARGUMENT LIST: -! ITYPE - INTEGER FLAG SPECIFYING HOW PARCEL TO LIFT IS -! IDENTIFIED. SEE COMMENTS ABOVE. -! DPBND - DEPTH OVER WHICH ONE SEARCHES FOR MOST UNSTABLE PARCEL -! P1D - ARRAY OF PRESSURE OF PARCELS TO LIFT. -! T1D - ARRAY OF TEMPERATURE OF PARCELS TO LIFT. -! Q1D - ARRAY OF SPECIFIC HUMIDITY OF PARCELS TO LIFT. -! L1D - ARRAY OF MODEL LEVEL OF PARCELS TO LIFT. -! -! OUTPUT ARGUMENT LIST: -! CAPE - CONVECTIVE AVAILABLE POTENTIAL ENERGY (J/KG) -! CINS - CONVECTIVE INHIBITION (J/KG) -! PPARC - PRESSURE LEVEL OF PARCEL LIFTED WHEN ONE SEARCHES -! OVER A PARTICULAR DEPTH TO COMPUTE CAPE/CIN -! -! OUTPUT FILES: -! STDOUT - RUN TIME STANDARD OUT. -! -! SUBPROGRAMS CALLED: -! UTILITIES: -! BOUND - BOUND (CLIP) DATA BETWEEN UPPER AND LOWER LIMTS. -! TTBLEX - LOOKUP TABLE ROUTINE TO GET T FROM THETAE AND P -! -! LIBRARY: -! COMMON - -! -! ATTRIBUTES: -! LANGUAGE: FORTRAN 90 -! MACHINE : CRAY C-90 -!$$$ -! -! use vrbls3d, only: pmid, t, q, zint -! use vrbls2d, only: teql -! use masks, only: lmh - - use params_mod, only: d00, h1m12, h99999, h10e5, capa, elocp, eps, & - oneps, g - use lookup_mod, only: thl, rdth, jtb, qs0, sqs, rdq, itb, ptbl, & - plq, ttbl, pl, rdp, the0, sthe, rdthe, ttblq, & - itbq, jtbq, rdpq, the0q, stheq, rdtheq - use ctlblk_mod, only: jsta_2l, jend_2u, lm, jsta, jend, im, me, spval -! -!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none -! -! INCLUDE/SET PARAMETERS. CONSTANTS ARE FROM BOLTON (MWR, 1980). - real,PARAMETER :: ISMTHP=2,ISMTHT=2,ISMTHQ=2 -! -! DECLARE PASSED VARIABLES -! - real, intent(in) :: LPAR0 - real, intent(in) :: PSFC - real, intent(in) :: LMASK - real, dimension(1:LM), intent(in) :: P_A, T_A, Q_A, ZINT_A - - real, intent(out) :: CAPE, CINS - real, dimension(1:LM), intent(out) :: TPAR - real, dimension(4), intent(inout) :: PARCEL0 - integer, intent(out) :: LLCL, LEQL -! -! DECLARE LOCAL VARIABLES -! -! integer, dimension(im,jsta:jend) :: IPTB, ITHTB, PARCEL, KLRES, KHRES, LCL, IDX - integer :: IPTB, ITHTB, KLRES, KHRES -! -! real, dimension(im,jsta:jend) :: THESP, PSP, CAPE20, QQ, PP, THUND - real :: THESP, PSP, QQ, PP - - LOGICAL THUNDER(IM,jsta:jend), NEEDTHUN - real PSFCK,PKL,TBTK,QBTK,APEBTK,TTHBTK,TTHK,APESPK,TPSPK, & - BQS00K,SQS00K,BQS10K,SQS10K,BQK,SQK,TQK,PRESK,GDZKL,THETAP, & - THETAA,P00K,P10K,P01K,P11K,TTHESK,ESATP,QSATP,TVP,TV - real DPBND - integer LP0 -! real,external :: fpvsnew - integer I,J,L,KNUML,KNUMH,LBEG,LEND,IQ, KB,ITTBK - integer PTYPE, PARCEL_L - -! integer I,J,L,KNUML,KNUMH,LBEG,LEND,IQ,IT,LMHK, KB,ITTBK -! -!************************************************************** -! START CALCAPE1D HERE. -! -! COMPUTE CAPE/CINS -! -! WHICH IS: THE SUM FROM THE LCL TO THE EQ LEVEL OF -! G * (LN(THETAP) - LN(THETAA)) * DZ -! -! (POSITIVE AREA FOR CAPE, NEGATIVE FOR CINS) -! -! WHERE: -! THETAP IS THE PARCEL THETA -! THETAA IS THE AMBIENT THETA -! DZ IS THE THICKNESS OF THE LAYER -! -! USING LCL AS LEVEL DIRECTLY BELOW SATURATION POINT -! AND EQ LEVEL IS THE HIGHEST POSITIVELY BUOYANT LEVEL. -! -! IEQL = EQ LEVEL -! P_thetaemax - real pressure of theta-e max parcel (Pa) -! -! INITIALIZE CAPE AND CINS ARRAYS -! - CAPE = D00 - CINS = D00 - PARCEL0(1:4) = D00 - LLCL = 0 - LEQL = LM - THESP = D00 - PSP = D00 -! - DO L=1,LM -! TPAR(L) = D00 - TPAR(L) = -9999.0 - ENDDO -! -! BASED ON LP0, CHOOSE THE WAY TO SELECT THE AIR PARCEL: -! - PTYPE = 0 - DPBND = 300.00 !DEFAULT GIVEN DEPTH - IF (LPAR0 <= -100.0) THEN -!-------FOR PTYPE=1--FIND MAXIMUM THETA E LAYER IN LOWEST DPBND ABOVE GROUND------- - PTYPE = 1 - DPBND = ABS(LPAR0) !GIVEN DEPTH - LP0 = -1 !DUMMY SINCE PARCEL IS IN SEARCHING. - ELSE IF ( NINT(LPAR0) == -9 ) THEN !GIVEN PARCEL -!-------FOR PTYPE=2--AIR PARCEL IS GIVEN BY PARCEL0------- - PTYPE = 2 - LP0 = 1 - ELSE IF ( NINT(LPAR0) == -1 ) THEN !SURFACE BASED PARCEL - PTYPE = 0 - LP0 = LM - ELSE IF ( NINT(LPAR0) == 0 ) THEN !IF LPAR0 IS 0, SET IT TO MODEL TOP - PTYPE = 0 - LP0 = 1 - ELSE IF ( NINT(LPAR0) >= LM ) THEN !IF LPAR0 IS BEYOND MODEL BASE, SET IT TO BASE LEVEL - PTYPE = 0 - LP0 = LM - ELSE IF ( NINT(LPAR0) >=1 .AND. NINT(LPAR0) <=LM ) THEN !LPAR0 IS IN GOOD RANGE - PTYPE = 0 - LP0 = NINT(LPAR0) - ELSE !INVALID VALUE OF LPAR0 - WRITE(0,*) "CALCAPE1D: INVALID VALUE PASSED TO INITIAL LEVEL OF AIR PARCEL" - STOP 1 - END IF - PARCEL_L = LP0 !PARCEL_L TO STORE LEVEL OF FOUNDED PARCEL(PTYPE=1) -! -!-------FOR ITYPE=1--FIND MAXIMUM THETA E LAYER IN LOWEST DPBND ABOVE GROUND------- -!--------------TRIAL MAXIMUM BUOYANCY LEVEL VARIABLES------------------- -!-------FOR ITYPE=2--FIND THETA E LAYER OF GIVEN T1D, Q1D, P1D--------------------- - - PSFCK = PSFC - DO KB=1,LM - PKL = P_A(KB) - IF ( (PTYPE == 0 .AND. KB == LP0) .OR. & - (PTYPE == 2 .AND. KB == LP0) .OR. & - (PTYPE == 1 .AND. (PKL >= PSFCK-DPBND .AND. PKL <= PSFCK)) ) THEN - - IF (PTYPE == 2 ) THEN !USING INFO OF GIVEN PARCEL - PKL = PARCEL0(2) - TBTK = PARCEL0(3) - QBTK = MAX(0.0, MIN(MAX(H1M12,PARCEL0(4)),H99999)) - APEBTK = (H10E5/PKL)**CAPA - ELSE !PARCEL IS CHOSEN FROM AMBIENT ATMOS. - PKL = P_A(KB) - TBTK = T_A(KB) - QBTK = max(0.0, Q_A(KB)) - APEBTK = (H10E5/PKL)**CAPA - ENDIF - -!----------Breogan Gomez - 2009-02-06 -! To prevent QBTK to be less than 0 which leads to a unrealistic value of PRESK -! and a floating invalid. - -! if(QBTK < 0) QBTK = 0 - -!--------------SCALING POTENTIAL TEMPERATURE & TABLE INDEX-------------- - TTHBTK = TBTK*APEBTK - TTHK = (TTHBTK-THL)*RDTH - QQ = TTHK - AINT(TTHK) - ITTBK = INT(TTHK) + 1 -!--------------KEEPING INDICES WITHIN THE TABLE------------------------- - IF(ITTBK < 1) THEN - ITTBK = 1 - QQ = D00 - ENDIF - IF(ITTBK >= JTB) THEN - ITTBK = JTB-1 - QQ = D00 - ENDIF -!--------------BASE AND SCALING FACTOR FOR SPEC. HUMIDITY--------------- - BQS00K = QS0(ITTBK) - SQS00K = SQS(ITTBK) - BQS10K = QS0(ITTBK+1) - SQS10K = SQS(ITTBK+1) -!--------------SCALING SPEC. HUMIDITY & TABLE INDEX--------------------- - BQK = (BQS10K-BQS00K)*QQ + BQS00K - SQK = (SQS10K-SQS00K)*QQ + SQS00K - TQK = (QBTK-BQK)/SQK*RDQ - PP = TQK-AINT(TQK) - IQ = INT(TQK)+1 -!--------------KEEPING INDICES WITHIN THE TABLE------------------------- - IF(IQ < 1) THEN - IQ = 1 - PP = D00 - ENDIF - IF(IQ >= ITB) THEN - IQ = ITB-1 - PP = D00 - ENDIF -!--------------SATURATION PRESSURE AT FOUR SURROUNDING TABLE PTS.------- - P00K = PTBL(IQ ,ITTBK ) - P10K = PTBL(IQ+1,ITTBK ) - P01K = PTBL(IQ ,ITTBK+1) - P11K = PTBL(IQ+1,ITTBK+1) -!--------------SATURATION POINT VARIABLES AT THE BOTTOM----------------- - TPSPK = P00K + (P10K-P00K)*PP + (P01K-P00K)*QQ & - + (P00K-P10K-P01K+P11K)*PP*QQ - -!!from WPP::tgs APESPK=(H10E5/TPSPK)**CAPA - if (TPSPK > 1.0e-3) then - APESPK = (max(0.,H10E5/ TPSPK))**CAPA - else - APESPK = 0.0 - endif - - TTHESK = TTHBTK * EXP(ELOCP*QBTK*APESPK/TTHBTK) -!--------------CHECK FOR MAXIMUM THETA E-------------------------------- - IF(TTHESK > THESP) THEN - PSP = TPSPK - THESP = TTHESK - PARCEL_L = KB ! MARK THE LEVEL OF MAX THETA-E - ENDIF - END IF ! KB == LP0 or PTYPE==1 (DPBND) - ENDDO ! KB loop - -!----FIND THE PRESSURE OF THE PARCEL THAT WAS LIFTED -!-----SAVE THE INITIAL PARCEL AND OUTPUT - IF (PTYPE == 0 .OR. PTYPE == 1) THEN - PARCEL0(1) = FLOAT(PARCEL_L) - PARCEL0(2) = P_A(PARCEL_L) - PARCEL0(3) = T_A(PARCEL_L) - PARCEL0(4) = Q_A(PARCEL_L) - ELSE IF (PTYPE == 1) THEN - PARCEL0(1) = FLOAT(LP0) - PARCEL0(2) = P_A(LP0) - PARCEL0(3) = T_A(LP0) - PARCEL0(4) = Q_A(LP0) - END IF -! -!-----CHOOSE LAYER DIRECTLY BELOW PSP AS LCL AND------------------------ -!-----ENSURE THAT THE LCL IS ABOVE GROUND.------------------------------ -!-------(IN SOME RARE CASES FOR ITYPE=2, IT IS NOT)--------------------- -!-----Note: L=1,LM --> from top to bottom - DO L=1,LM - IF (P_A(L) < PSP) LLCL = L+1 - ENDDO - IF (LLCL > NINT(LMASK)) LLCL = NINT(LMASK) -!----------------------------------------------------------------------- -!---------FIND TEMP OF PARCEL LIFTED ALONG MOIST ADIABAT (TPAR)--------- -!----------------------------------------------------------------------- -!-----Note: loop L=LM,1,-1 --> from bottom to top - DO L=LM,1,-1 - KLRES = 0 - KHRES = 0 -!--------------SCALING PRESSURE & TT TABLE INDEX------------------------ - IF(L <= LLCL) THEN - IF(P_A(L) < PLQ)THEN - KLRES = 1 - ELSE - KHRES = 1 - ENDIF - ENDIF -!*** -!*** COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE 0) THEN - CALL TTBLEX1D(TPAR(L),TTBL,ITB,JTB,KLRES,P_A(L), & - PL,QQ,PP,RDP,THE0,STHE,RDTHE,THESP, & - IPTB,ITHTB) - ENDIF -!*** -!*** COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE>PLQ -!** - IF(KHRES > 0) THEN - CALL TTBLEX1D(TPAR(L),TTBLQ,ITBQ,JTBQ,KHRES,P_A(L), & - PLQ,QQ,PP,RDPQ,THE0Q,STHEQ,RDTHEQ,THESP, & - IPTB,ITHTB) - ENDIF - -!------------SEARCH FOR EQ LEVEL---------------------------------------- - IF(KHRES > 0) THEN - IF(TPAR(L) > T_A(L)) LEQL = L - ENDIF -! - IF(KLRES > 0) THEN - IF(TPAR(L) > T_A(L) .AND. & - P_A(L)>100.) LEQL = L - ENDIF -!----------------------------------------------------------------------- - ENDDO ! end of do l=lm,1,-1 loop -!------------COMPUTE CAPE AND CINS-------------------------------------- - LBEG = 1000 - LEND = 0 - LBEG = MIN(LEQL,LBEG) - LEND = MAX(LLCL,LEND) -! -!-----Note: LOOP for CAPE: L=LEQL, LLCL --> FROM HIGH TO LOW - DO L=LBEG,LEND - - IF(L >= LEQL .AND. L <= LLCL) THEN - PRESK = P_A(L) - GDZKL = (ZINT_A(L)-ZINT_A(L+1)) * G - ESATP = min(FPVSNEW(TPAR(L)),PRESK) - QSATP = EPS*ESATP/(PRESK-ESATP*ONEPS) -! TVP = TPAR(I,J,L)*(1+0.608*QSATP) - TVP = TVIRTUAL(TPAR(L),QSATP) - THETAP = TVP*(H10E5/PRESK)**CAPA -! TV = T(I,J,L)*(1+0.608*Q(I,J,L)) - TV = TVIRTUAL(T_A(L),Q_A(L)) - THETAA = TV*(H10E5/PRESK)**CAPA - IF(THETAP < THETAA) THEN - CINS = CINS + (LOG(THETAP)-LOG(THETAA))*GDZKL - ELSEIF(THETAP > THETAA) THEN - CAPE = CAPE + (LOG(THETAP)-LOG(THETAA))*GDZKL - ENDIF - ENDIF - - ENDDO -! -! ENFORCE LOWER LIMIT OF 0.0 ON CAPE AND UPPER -! LIMIT OF 0.0 ON CINS. -! - CAPE = MAX(D00,CAPE) - CINS = MIN(CINS,D00) -! - RETURN -! - END SUBROUTINE CALCAPE1D ! !------------------------------------------------------------------------------------- ! From e1bb6a2fd1deaa4c5869f40554ac34fc847ccbfe Mon Sep 17 00:00:00 2001 From: EdwardColon-NOAA Date: Mon, 16 Aug 2021 15:51:49 -0400 Subject: [PATCH 09/14] Update MISCLN.f --- sorc/ncep_post.fd/MISCLN.f | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/sorc/ncep_post.fd/MISCLN.f b/sorc/ncep_post.fd/MISCLN.f index cf02fd640..0cb596bdb 100644 --- a/sorc/ncep_post.fd/MISCLN.f +++ b/sorc/ncep_post.fd/MISCLN.f @@ -2937,9 +2937,8 @@ SUBROUTINE MISCLN EGRID1(I,J) = LOG(PMID(I,J,LM)/EGRID2(I,J)) & / LOG(PMID(I,J,LM)/PMID(I,J,LM-1)) - IF (MODELNAME == 'GFS' .OR. MODELNAME == 'FV3R' & - .OR. MODELNAME == 'FV3R') THEN - EGRID1(I,J) = LOG(PMID(I,J,LM)/EGRID2(I,J)) & + IF (MODELNAME == 'GFS' .OR. MODELNAME == 'FV3R') THEN + EGRID1(I,J) = LOG(PMID(I,J,LM)/EGRID2(I,J)) & / max(1.e-6,LOG(PMID(I,J,LM)/PMID(I,J,LM-1))) EGRID1(I,J) =max(-10.0,min(EGRID1(I,J), 10.0)) IF ( ABS(PMID(I,J,LM)-PMID(I,J,LM-1)) < 0.5 ) THEN From 18b60ba8e3ac65d2c641b5b557ea84408928bb71 Mon Sep 17 00:00:00 2001 From: EdwardColon-NOAA Date: Mon, 16 Aug 2021 17:12:49 -0400 Subject: [PATCH 10/14] Update MISCLN.f --- sorc/ncep_post.fd/MISCLN.f | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/sorc/ncep_post.fd/MISCLN.f b/sorc/ncep_post.fd/MISCLN.f index 0cb596bdb..b00579605 100644 --- a/sorc/ncep_post.fd/MISCLN.f +++ b/sorc/ncep_post.fd/MISCLN.f @@ -95,8 +95,7 @@ SUBROUTINE MISCLN jsta_2l, jend_2u, MODELNAME, SUBMODELNAME use rqstfld_mod, only: iget, lvls, id, iavblfld, lvlsxml use grib2_module, only: pset - use upp_physics, only: FPVSNEW,CALRH_PW,CALCAPE,CALCAPE2,TVIRTUAL, & - CALCAPE1D + use upp_physics, only: FPVSNEW,CALRH_PW,CALCAPE,CALCAPE2,TVIRTUAL use gridspec_mod, only: gridtype !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none From 2d3fee79eef9a1c60704b11b8e2e2a45ca1b136a Mon Sep 17 00:00:00 2001 From: "Edward.Colon" Date: Fri, 20 Aug 2021 02:35:41 +0000 Subject: [PATCH 11/14] Added conditional block that only invokes the effective layer top and bottom correction for RTMA. --- sorc/ncep_post.fd/MISCLN.f | 40 +++++++++++++++++++++----------------- 1 file changed, 22 insertions(+), 18 deletions(-) diff --git a/sorc/ncep_post.fd/MISCLN.f b/sorc/ncep_post.fd/MISCLN.f index b00579605..08e9c7750 100644 --- a/sorc/ncep_post.fd/MISCLN.f +++ b/sorc/ncep_post.fd/MISCLN.f @@ -152,7 +152,6 @@ SUBROUTINE MISCLN MAXWP, MAXWZ, MAXWU, MAXWV, & MAXWT INTEGER,dimension(:,:),allocatable :: LLOW, LUPP - INTEGER,dimension(:,:),allocatable :: LLOW1, LUPP1 REAL, dimension(:,:),allocatable :: CANGLE,ESHR,UVECT,VVECT,& EFFUST,EFFVST,FSHR,HTSFC,& ESRH @@ -3455,6 +3454,9 @@ SUBROUTINE MISCLN ENDIF ENDIF + + IF (SUBMODELNAME == 'RTMA')THEN + ! ! --- Effective (inflow) Layer (EL) ! @@ -3560,6 +3562,8 @@ SUBROUTINE MISCLN IF(ALLOCATED(TPAR_BASE)) DEALLOCATE(TPAR_BASE) IF(ALLOCATED(TPAR_TOPS)) DEALLOCATE(TPAR_TOPS) + + ENDIF ! ! EXPAND HRRR CAPE/CIN RELATED VARIABLES ! @@ -3709,7 +3713,6 @@ SUBROUTINE MISCLN HELI(IM,jsta_2l:jend_2u,2)) allocate(LLOW(IM,jsta_2l:jend_2u),LUPP(IM,jsta_2l:jend_2u), & CANGLE(IM,jsta_2l:jend_2u)) -! allocate(LLOW1(IM,jsta_2l:jend_2u),LUPP1(IM,jsta_2l:jend_2u)) iget1 = IGET(953) iget2 = -1 @@ -3722,13 +3725,25 @@ SUBROUTINE MISCLN IF (iget1 > 0 .OR. IGET(162) > 0 .OR. IGET(953) > 0) THEN DEPTH(1) = 3000.0 DEPTH(2) = 1000.0 + IF (SUBMODELNAME == 'RTMA') THEN +!--- IF USSING EL BASE & TOP COMPUTED BY NEW SCHEME FOR THE +!RELATED VARIABLES !$omp parallel do private(i,j) - DO J=JSTA,JEND - DO I=1,IM - LLOW(I,J) = INT(EGRID4(I,J)) - LUPP(I,J) = INT(EGRID5(I,J)) + DO J=JSTA,JEND + DO I=1,IM + LLOW(I,J) = EL_BASE(I,J) + LUPP(I,J) = EL_TOPS(I,J) + ENDDO ENDDO - ENDDO + ELSE +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=1,IM + LLOW(I,J) = INT(EGRID4(I,J)) + LUPP(I,J) = INT(EGRID5(I,J)) + ENDDO + ENDDO + ENDIF !--- OUTPUT EL BASE & TOP COMPUTED BY OLD SCHEME IF (debugprint) THEN WRITE(IM_CH,'(I5.5)') IM @@ -3751,15 +3766,6 @@ SUBROUTINE MISCLN CLOSE(IUNIT) ENDIF -!--- IF USSING EL BASE & TOP COMPUTED BY NEW SCHEME FOR THE RELATED VARIABLES -!$omp parallel do private(i,j) - DO J=JSTA,JEND - DO I=1,IM - LLOW(I,J) = EL_BASE(I,J) - LUPP(I,J) = EL_TOPS(I,J) - ENDDO - ENDDO - ! CALL CALHEL(DEPTH,UST,VST,HELI,USHR1,VSHR1,USHR6,VSHR6) CALL CALHEL2(LLOW,LUPP,DEPTH,UST,VST,HELI,CANGLE) @@ -4496,8 +4502,6 @@ SUBROUTINE MISCLN if (allocated(heli)) deallocate(heli) if (allocated(llow)) deallocate(llow) if (allocated(lupp)) deallocate(lupp) - if (allocated(llow1)) deallocate(llow1) - if (allocated(lupp1)) deallocate(lupp1) if (allocated(cangle))deallocate(cangle) if (allocated(effust))deallocate(effust) if (allocated(effvst))deallocate(effvst) From 3ac68ca903b4e9eb14a3d87f7cdfd005268ef01d Mon Sep 17 00:00:00 2001 From: "Edward.Colon" Date: Sun, 22 Aug 2021 00:53:21 +0000 Subject: [PATCH 12/14] Addressed remaining bug linked to the use of the FV3R first guess. --- sorc/ncep_post.fd/ALLOCATE_ALL.f | 6 +++--- sorc/ncep_post.fd/INITPOST.F | 4 ++-- sorc/ncep_post.fd/MISCLN.f | 19 +++++++++---------- 3 files changed, 14 insertions(+), 15 deletions(-) diff --git a/sorc/ncep_post.fd/ALLOCATE_ALL.f b/sorc/ncep_post.fd/ALLOCATE_ALL.f index 653adb1ef..a0acb04c6 100644 --- a/sorc/ncep_post.fd/ALLOCATE_ALL.f +++ b/sorc/ncep_post.fd/ALLOCATE_ALL.f @@ -1231,7 +1231,7 @@ SUBROUTINE ALLOCATE_ALL() !Initialization !$omp parallel do private(i,j) do j=jsta_2l,jend_2u - do i=1,lm + do i=1,im dusmass(i,j)=spval ducmass(i,j)=spval dusmass25(i,j)=spval @@ -1273,7 +1273,7 @@ SUBROUTINE ALLOCATE_ALL() !Initialization !$omp parallel do private(i,j) do j=jsta_2l,jend_2u - do i=1,lm + do i=1,im acswupt(i,j)=spval swdnt(i,j)=spval acswdnt(i,j)=spval @@ -1287,7 +1287,7 @@ SUBROUTINE ALLOCATE_ALL() !Initialization !$omp parallel do private(i,j) do j=jsta_2l,jend_2u - do i=1,lm + do i=1,im ddvdx(i,j)=spval ddudy(i,j)=spval uuavg(i,j)=spval diff --git a/sorc/ncep_post.fd/INITPOST.F b/sorc/ncep_post.fd/INITPOST.F index 8a8991bc8..87cd31af5 100644 --- a/sorc/ncep_post.fd/INITPOST.F +++ b/sorc/ncep_post.fd/INITPOST.F @@ -1781,7 +1781,7 @@ SUBROUTINE INITPOST IM,1,JM,1,IM,JS,JE,1) do j = jsta_2l, jend_2u do i = 1, im - IF(SUBMODELNAME == 'RTMA')THEN !use 1st level of unstaggered U for U10 + IF(SUBMODELNAME == 'RTMA' .and. MODELNAME == 'RAPR')THEN !use 1st level of unstaggered U for U10 U10 ( i, j ) = uh ( i, j, lm ) ELSE U10 ( i, j ) = dummy( i, j ) @@ -1793,7 +1793,7 @@ SUBROUTINE INITPOST IM,1,JM,1,IM,JS,JE,1) do j = jsta_2l, jend_2u do i = 1, im - IF( SUBMODELNAME == 'RTMA')THEN!use 1st level of unstaggered V for V10 + IF( SUBMODELNAME == 'RTMA' .and. MODELNAME == 'RAPR')THEN!use 1st level of unstaggered V for V10 V10 ( i, j ) = vh ( i, j, lm ) ELSE V10 ( i, j ) = dummy( i, j ) diff --git a/sorc/ncep_post.fd/MISCLN.f b/sorc/ncep_post.fd/MISCLN.f index 08e9c7750..90a3f58aa 100644 --- a/sorc/ncep_post.fd/MISCLN.f +++ b/sorc/ncep_post.fd/MISCLN.f @@ -1425,7 +1425,7 @@ SUBROUTINE MISCLN DO J=JSTA,JEND DO I=1,IM GRID1(I,J)=Z1D(I,J) - IF (SUBMODELNAME == 'RTMA' .OR. MODELNAME == 'FV3R') THEN + IF (SUBMODELNAME == 'RTMA') THEN FREEZELVL(I,J)=GRID1(I,J) ENDIF ENDDO @@ -3150,7 +3150,7 @@ SUBROUTINE MISCLN DO I=1,IM IF(T1D(I,J) < spval) THEN GRID1(I,J) = EGRID1(I,J) - IF (SUBMODELNAME == 'RTMA' .OR. MODELNAME == 'FV3R') THEN + IF (SUBMODELNAME == 'RTMA') THEN MLCAPE(I,J)=GRID1(I,J) ENDIF ENDIF @@ -3187,7 +3187,7 @@ SUBROUTINE MISCLN DO I=1,IM IF(T1D(I,J) < spval) THEN GRID1(I,J) = - GRID1(I,J) - IF (SUBMODELNAME == 'RTMA' .OR. MODELNAME == 'FV3R') THEN + IF (SUBMODELNAME == 'RTMA') THEN MLCIN(I,J) = GRID1(I,J) ENDIF ENDIF @@ -3219,7 +3219,7 @@ SUBROUTINE MISCLN DO J=JSTA,JEND DO I=1,IM IF(T1D(I,J) < spval) GRID1(I,J)=EGRID2(I,J) - IF (SUBMODELNAME == 'RTMA' .OR. MODELNAME == 'FV3R') THEN + IF (SUBMODELNAME == 'RTMA') THEN MLLCL(I,J) = GRID1(I,J) ENDIF ENDDO @@ -3279,7 +3279,7 @@ SUBROUTINE MISCLN DPBND = 300.E2 CALL CALCAPE(ITYPE,DPBND,P1D,T1D,Q1D,LB2,EGRID1, & EGRID2,EGRID3,EGRID4,EGRID5) -! IF (SUBMODELNAME == 'RTMA' .OR. MODELNAME == 'FV3R') THEN +! IF (SUBMODELNAME == 'RTMA') THEN ! MUMIXR(I,J) = Q1D(I,J) ! ENDIF IF (IGET(584)>0) THEN @@ -3290,14 +3290,14 @@ SUBROUTINE MISCLN DO I=1,IM IF(T1D(I,J) < spval) THEN GRID1(I,J) = EGRID1(I,J) - IF (SUBMODELNAME == 'RTMA' .OR. MODELNAME == 'FV3R') THEN + IF (SUBMODELNAME == 'RTMA') THEN MUCAPE(I,J) = GRID1(I,J) ENDIF ENDIF ENDDO ENDDO CALL BOUND(GRID1,D00,H99999) -! IF (SUBMODELNAME == 'RTMA' .OR. MODELNAME == 'FV3R') THEN +! IF (SUBMODELNAME == 'RTMA') THEN ! CALL BOUND(MUCAPE,D00,H99999) ! ENDIF if(grib=='grib2') then @@ -3329,8 +3329,7 @@ SUBROUTINE MISCLN DO I=1,IM IF(T1D(I,J) < spval) THEN GRID1(I,J) = - GRID1(I,J) - IF (SUBMODELNAME == 'RTMA' .OR. & - MODELNAME == 'FV3R') THEN + IF (SUBMODELNAME == 'RTMA')THEN MUCAPE(I,J) = GRID1(I,J) MUQ1D(I,J) = Q1D(I,J) ENDIF @@ -3795,7 +3794,7 @@ SUBROUTINE MISCLN ENDIF !953 - IF (SUBMODELNAME == 'RTMA' .OR. MODELNAME=='FV3R') THEN !Start RTMA block + IF (SUBMODELNAME == 'RTMA') THEN !Start RTMA block !EL field allocation From 574b0d8824cfc6d0eabf690e0cec500f39095986 Mon Sep 17 00:00:00 2001 From: "Edward.Colon" Date: Tue, 31 Aug 2021 14:55:47 +0000 Subject: [PATCH 13/14] minor code changes made to ALLOCATE_ALL.f, MISCNL.f and UPP_PHYSICS.f to address segmentation fault occurring during the rtma regression tests --- sorc/ncep_post.fd/ALLOCATE_ALL.f | 2 +- sorc/ncep_post.fd/MISCLN.f | 41 ++++++++++++-------------------- sorc/ncep_post.fd/UPP_PHYSICS.f | 2 +- 3 files changed, 17 insertions(+), 28 deletions(-) diff --git a/sorc/ncep_post.fd/ALLOCATE_ALL.f b/sorc/ncep_post.fd/ALLOCATE_ALL.f index a0acb04c6..f444e518f 100644 --- a/sorc/ncep_post.fd/ALLOCATE_ALL.f +++ b/sorc/ncep_post.fd/ALLOCATE_ALL.f @@ -695,7 +695,7 @@ SUBROUTINE ALLOCATE_ALL() t700(i,j)=spval z700(i,j)=spval teql(i,j)=spval - ieql(i,j)=spval + ieql(i,j)=0 cfracl(i,j)=spval cfracm(i,j)=spval cfrach(i,j)=spval diff --git a/sorc/ncep_post.fd/MISCLN.f b/sorc/ncep_post.fd/MISCLN.f index 90a3f58aa..f051ae865 100644 --- a/sorc/ncep_post.fd/MISCLN.f +++ b/sorc/ncep_post.fd/MISCLN.f @@ -3148,12 +3148,10 @@ SUBROUTINE MISCLN !$omp parallel do private(i,j) DO J=JSTA,JEND DO I=1,IM - IF(T1D(I,J) < spval) THEN - GRID1(I,J) = EGRID1(I,J) - IF (SUBMODELNAME == 'RTMA') THEN - MLCAPE(I,J)=GRID1(I,J) - ENDIF - ENDIF + IF(T1D(I,J) < spval) THEN + GRID1(I,J) = EGRID1(I,J) + IF (SUBMODELNAME == 'RTMA')MLCAPE(I,J)=GRID1(I,J) + ENDIF ENDDO ENDDO CALL BOUND(GRID1,D00,H99999) @@ -3185,12 +3183,10 @@ SUBROUTINE MISCLN !$omp parallel do private(i,j) DO J=JSTA,JEND DO I=1,IM - IF(T1D(I,J) < spval) THEN - GRID1(I,J) = - GRID1(I,J) - IF (SUBMODELNAME == 'RTMA') THEN - MLCIN(I,J) = GRID1(I,J) - ENDIF - ENDIF + IF(T1D(I,J) < spval) THEN + GRID1(I,J) = - GRID1(I,J) + IF (SUBMODELNAME == 'RTMA') MLCIN(I,J)=GRID1(I,J) + ENDIF ENDDO ENDDO ! @@ -3219,9 +3215,7 @@ SUBROUTINE MISCLN DO J=JSTA,JEND DO I=1,IM IF(T1D(I,J) < spval) GRID1(I,J)=EGRID2(I,J) - IF (SUBMODELNAME == 'RTMA') THEN - MLLCL(I,J) = GRID1(I,J) - ENDIF + IF (SUBMODELNAME == 'RTMA') MLLCL(I,J) = GRID1(I,J) ENDDO ENDDO ! @@ -3279,21 +3273,17 @@ SUBROUTINE MISCLN DPBND = 300.E2 CALL CALCAPE(ITYPE,DPBND,P1D,T1D,Q1D,LB2,EGRID1, & EGRID2,EGRID3,EGRID4,EGRID5) -! IF (SUBMODELNAME == 'RTMA') THEN -! MUMIXR(I,J) = Q1D(I,J) -! ENDIF + IF (SUBMODELNAME == 'RTMA') MUMIXR(I,J) = Q1D(I,J) IF (IGET(584)>0) THEN ! dong add missing value to cin GRID1 = spval !$omp parallel do private(i,j) DO J=JSTA,JEND DO I=1,IM - IF(T1D(I,J) < spval) THEN - GRID1(I,J) = EGRID1(I,J) - IF (SUBMODELNAME == 'RTMA') THEN - MUCAPE(I,J) = GRID1(I,J) - ENDIF - ENDIF + IF(T1D(I,J) < spval) THEN + GRID1(I,J) = EGRID1(I,J) + IF (SUBMODELNAME == 'RTMA') MUCAPE(I,J)=GRID1(I,J) + ENDIF ENDDO ENDDO CALL BOUND(GRID1,D00,H99999) @@ -3939,8 +3929,7 @@ SUBROUTINE MISCLN GRID1=spval DO J=JSTA,JEND DO I=1,IM - IF(LLOW(I,J) Date: Wed, 1 Sep 2021 22:53:03 +0000 Subject: [PATCH 14/14] Clean up before merging. --- sorc/ncep_post.fd/MISCLN.f | 3 + sorc/ncep_post.fd/TTBLEX.f | 124 -------------------------------- sorc/ncep_post.fd/UPP_PHYSICS.f | 2 + 3 files changed, 5 insertions(+), 124 deletions(-) diff --git a/sorc/ncep_post.fd/MISCLN.f b/sorc/ncep_post.fd/MISCLN.f index f051ae865..c333ad586 100644 --- a/sorc/ncep_post.fd/MISCLN.f +++ b/sorc/ncep_post.fd/MISCLN.f @@ -45,6 +45,9 @@ !! 20-11-10 J Meng - USE UPP_PHYSICS MODULE !! 21-03-25 E Colon - 3D-RTMA-specific SPC fields added as output !! 21-04-01 J Meng - computation on defined points only +!! 21-09-01 E Colon - Correction to the effective layer top and +!! bottoma calculation which is only employed +!! for RTMA usage. !! !! USAGE: CALL MISCLN !! INPUT ARGUMENT LIST: diff --git a/sorc/ncep_post.fd/TTBLEX.f b/sorc/ncep_post.fd/TTBLEX.f index 0e11d5da5..21748a6f4 100644 --- a/sorc/ncep_post.fd/TTBLEX.f +++ b/sorc/ncep_post.fd/TTBLEX.f @@ -116,128 +116,4 @@ END SUBROUTINE TTBLEX !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& ! !------------------------------------------------------------------------------------- -! - SUBROUTINE TTBLEX1D(TREF,TTBL,ITB,JTB,KARR,PMIDL, & - PL,QQ,PP,RDP,THE0,STHE,RDTHE,THESP, & - IPTB,ITHTB) -!FPP$ NOCONCUR R -!$$$ SUBPROGRAM DOCUMENTATION BLOCK -! . . . -! SUBPROGRAM: TTBLEX COMPUTES T ALONG A MOIST ADIABAT -! PRGRMMR: BLACK ORG: W/NP2 DATE: ??-??-?? -! -! ABSTRACT: -! THIS ROUTINE COMPUTES THE TEMPERATURE ALONG A MOIST -! ADIABAT GIVEN THE SATURATION POTENTIAL TEMPERATURE -! AND THE PRESSURE -! . -! -! PROGRAM HISTORY LOG: -! ??-??-?? T BLACK - ORIGINATOR -! 98-06-12 T BLACK - CONVERSION FROM 1-D TO 2-D -! 00-01-04 JIM TUCCILLO - MPI VERSION -! 01-10-22 H CHUANG - MODIFIED TO PROCESS HYBRID MODEL OUTPUT -! 02-01-15 MIKE BALDWIN - WRF VERSION -! -! OUTPUT FILES: -! NONE -! -! SUBPROGRAMS CALLED: -! UTILITIES: -! NONE -! -! ATTRIBUTES: -! LANGUAGE: FORTRAN -!---------------------------------------------------------------------- - use ctlblk_mod, only: jsta, jend, im, jsta_2l, jend_2u, me -!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none -!---------------------------------------------------------------------- - - integer,intent(in) :: ITB,JTB - integer,intent(in) :: KARR - real,dimension(JTB,ITB),intent(in) :: TTBL - real, intent(in) :: PMIDL - real, intent(out) :: TREF -! real,dimension(IM,jsta:jend),intent(out) :: QQ,PP - real,intent(out) :: QQ,PP - real, intent(in) :: THESP - real,dimension(ITB), intent(in) :: THE0,STHE -! integer,dimension(IM,jsta:jend),intent(out) :: IPTB,ITHTB - integer,intent(out) :: IPTB,ITHTB - real,intent(in) :: PL,RDP,RDTHE - -! - integer I,J,ITH,IP,IPTBK - real PK,TPK,T00K,T10K,T01K,T11K,BTHE00K,STHE00K,BTHK,STHK, & - TTHK,BTHE10K,STHE10K -! real QQ, PP -! INTEGER IPTB, ITHTB -!----------------------------------------------------------------------- - -!!$omp parallel do & -!!$omp& private(i,j,bthe00k,bthe10k,bthk,ip,iptbk,ith,pk,sthe00k,sthe10k,& -!!$omp& sthk,t00k,t01k,t10k,t11k,tpk,tthk) -! DO J=JSTA,JEND -! DO I=1,IM - IF(KARR > 0) THEN -!--------------SCALING PRESSURE & TT TABLE INDEX------------------------ - PK = PMIDL - TPK = (PK-PL)*RDP - QQ = TPK-AINT(TPK) - IPTB = INT(TPK) + 1 -!--------------KEEPING INDICES WITHIN THE TABLE------------------------- - IF(IPTB < 1) THEN - IPTB = 1 - QQ = 0. - ENDIF -! - IF(IPTB >= ITB) THEN - IPTB = ITB-1 - QQ = 0. - ENDIF -!--------------BASE AND SCALING FACTOR FOR THE-------------------------- - IPTBK = IPTB - BTHE00K = THE0(IPTBK) - STHE00K = STHE(IPTBK) - BTHE10K = THE0(IPTBK+1) - STHE10K = STHE(IPTBK+1) -!--------------SCALING THE & TT TABLE INDEX----------------------------- - BTHK = (BTHE10K-BTHE00K)*QQ+BTHE00K - STHK = (STHE10K-STHE00K)*QQ+STHE00K - TTHK = (THESP-BTHK)/STHK*RDTHE - PP = TTHK-AINT(TTHK) -! write(1000+me,*)' i=',i,' j=',j,' tthk=',tthk,' thesp=',thesp(i,j) & -! , ' bthk=',bthk,' sthk=',sthk,' rdthe=',rdthe - - ITHTB = INT(TTHK)+1 -!--------------KEEPING INDICES WITHIN THE TABLE------------------------- - IF(ITHTB < 1) THEN - ITHTB = 1 - PP = 0. - ENDIF -! - IF(ITHTB >= JTB) THEN - ITHTB = JTB-1 - PP = 0. - ENDIF -!--------------TEMPERATURE AT FOUR SURROUNDING TT TABLE PTS.------------ - ITH = ITHTB - IP = IPTB - T00K = TTBL(ITH ,IP ) - T10K = TTBL(ITH+1,IP ) - T01K = TTBL(ITH ,IP+1) - T11K = TTBL(ITH+1,IP+1) -!--------------PARCEL TEMPERATURE------------------------------------- - TREF = (T00K+(T10K-T00K)*PP+(T01K-T00K)*QQ & - + (T00K-T10K-T01K+T11K)*PP*QQ) - ENDIF -! ENDDO -! ENDDO -! - RETURN - END SUBROUTINE TTBLEX1D -!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& -! -!------------------------------------------------------------------------------------- ! diff --git a/sorc/ncep_post.fd/UPP_PHYSICS.f b/sorc/ncep_post.fd/UPP_PHYSICS.f index 0b499dc7b..1e7ad92c5 100644 --- a/sorc/ncep_post.fd/UPP_PHYSICS.f +++ b/sorc/ncep_post.fd/UPP_PHYSICS.f @@ -574,6 +574,7 @@ SUBROUTINE CALCAPE(ITYPE,DPBND,P1D,T1D,Q1D,L1D,CAPE, & ! 10-09-09 G MANIKIN - CHANGED COMPUTATION TO USE VIRTUAL TEMP ! - ADDED EQ LVL HGHT AND THUNDER PARAMETER ! 15-xx-xx S MOORTHI - optimization and threading +! 21-09-01 E COLON - equivalent level height index for RTMA ! ! USAGE: CALL CALCAPE(ITYPE,DPBND,P1D,T1D,Q1D,L1D,CAPE, ! CINS,PPARC) @@ -1077,6 +1078,7 @@ SUBROUTINE CALCAPE2(ITYPE,DPBND,P1D,T1D,Q1D,L1D, & ! 19-09-03 J MENG - MODIFIED TO ADD 0-3KM CAPE/CINS, LFC, ! EFFECTIVE HELICITY, DOWNDRAFT CAPE, ! DENDRITIC GROWTH LAYER DEPTH, ESP +! 21-09-01 E COLON - equivalent level height index for RTMA ! ! USAGE: CALL CALCAPE2(ITYPE,DPBND,P1D,T1D,Q1D,L1D, & ! CAPE,CINS,LFC,ESRHL,ESRHH, &