From 29341ff83ab98dda3e104e21b1120672d4f16373 Mon Sep 17 00:00:00 2001 From: JesseMeng-NOAA <49207433+JesseMeng-NOAA@users.noreply.github.com> Date: Mon, 28 Dec 2020 19:43:56 -0500 Subject: [PATCH] UPP refactor phase 2 from Jesse Meng (#224) * 20200820 Jesse Meng refactor branch of remove goto remove grib1 * 20200901 Jesse Meng remove redundant LAI output in SURFCE.f * 20200901 Jesse Meng continue merging refactor to develop * 20200911 Jesse Meng add precisions to selected variables * 20200911 Jesse Meng revert the precision changes * 20200924 Jesse Meng update with EMC_post/develop * 20201105 Jesse Meng merge with NOAA-EMC/EMC_post PR#196 https://github.com/NOAA-EMC/EMC_post/commit/eed90f5d8af3e1734bd489b249398774000a3eef * 20201105 Jesse Meng sync with EMC_POST/develop * 20201106 Jesse Meng add module UPP_MATH * 20201113 Jesse Meng add UPP_MATH and UPP_PHYSICS modules and updated with https://github.com/NOAA-EMC/EMC_post/pull/219 * 20201113 Jesse Meng updated CMakeLists.txt with EMC_post/develop * 20201113 Jesse Meng fix sorc/ncep_post.fd/CMakeLists.txt for added/removed files * 20201130 Jesse Meng update modules, move GRIDAVG.f inside module UPP_MATH * 20201215 Jesse Meng fix EXCH() options in CALVOR and UPP_MATH * 20201215 Jesse Meng merge with EMC_post/develop * 20201215 Jesse Meng minor clean up as suggested by reviewers * 20201218 Jesse Meng revised RH unification fv3, fv3r, nmmb, use CALRH_NAM hrrr, rap, use CALRH_GSD * 20201222 Jesse Meng convert UPP_MATH UPP_PHYSICS docblocks to doxygen format * 20201223 Jesse Meng minor revision style tweaking * 20201223 Jesse Meng fix redundance which costs nmmb runtime Co-authored-by: wx22mj Co-authored-by: wx22mj Co-authored-by: wx22mj Co-authored-by: wx22mj Co-authored-by: wx22mj Co-authored-by: wx22mj Co-authored-by: wx22mj --- sorc/ncep_post.fd/ALLOCATE_ALL.f | 6 + sorc/ncep_post.fd/BNDLYR.f | 3 +- sorc/ncep_post.fd/CALCAPE.f | 493 ----- sorc/ncep_post.fd/CALCAPE2.f | 788 -------- sorc/ncep_post.fd/CALPW.f | 3 +- sorc/ncep_post.fd/CALRH.f | 94 - sorc/ncep_post.fd/CALRH_GFS.f | 209 --- sorc/ncep_post.fd/CALRH_GSD.f | 37 - sorc/ncep_post.fd/CALRH_PW.f | 63 - sorc/ncep_post.fd/CALUPDHEL.f | 18 +- sorc/ncep_post.fd/CALVOR.f | 43 +- sorc/ncep_post.fd/CLDRAD.f | 7 +- sorc/ncep_post.fd/CMakeLists.txt | 9 +- sorc/ncep_post.fd/DEALLOCATE.f | 5 + sorc/ncep_post.fd/FDLVL.f | 3 +- sorc/ncep_post.fd/FRZLVL.f | 6 +- sorc/ncep_post.fd/FRZLVL2.f | 3 +- sorc/ncep_post.fd/INITPOST_GFS.f | 3 +- sorc/ncep_post.fd/INITPOST_GFS_NEMS.f | 3 +- sorc/ncep_post.fd/INITPOST_GFS_NEMS_MPIIO.f | 3 +- sorc/ncep_post.fd/INITPOST_GFS_NETCDF.f | 3 +- sorc/ncep_post.fd/INITPOST_GFS_NETCDF_PARA.f | 3 +- sorc/ncep_post.fd/INITPOST_GFS_SIGIO.f | 3 +- sorc/ncep_post.fd/INITPOST_NEMS.f | 1 + sorc/ncep_post.fd/INITPOST_NEMS_MPIIO.f | 1 + sorc/ncep_post.fd/INITPOST_NETCDF.f | 3 +- sorc/ncep_post.fd/LFMFLD.f | 3 +- sorc/ncep_post.fd/LFMFLD_GFS.f | 3 +- sorc/ncep_post.fd/MDL2P.f | 12 +- sorc/ncep_post.fd/MDL2STD_P.f | 14 +- sorc/ncep_post.fd/MDL2THANDPV.f | 27 +- sorc/ncep_post.fd/MDLFLD.f | 15 +- sorc/ncep_post.fd/MISCLN.f | 5 +- sorc/ncep_post.fd/OTLFT.f | 3 +- sorc/ncep_post.fd/OTLIFT.f | 3 +- sorc/ncep_post.fd/SURFCE.f | 15 +- sorc/ncep_post.fd/{GRIDAVG.f => UPP_MATH.f} | 142 +- sorc/ncep_post.fd/UPP_PHYSICS.f | 1775 ++++++++++++++++++ sorc/ncep_post.fd/makefile_module | 13 +- 39 files changed, 2029 insertions(+), 1814 deletions(-) delete mode 100644 sorc/ncep_post.fd/CALCAPE.f delete mode 100644 sorc/ncep_post.fd/CALCAPE2.f delete mode 100644 sorc/ncep_post.fd/CALRH.f delete mode 100644 sorc/ncep_post.fd/CALRH_GFS.f delete mode 100644 sorc/ncep_post.fd/CALRH_GSD.f delete mode 100644 sorc/ncep_post.fd/CALRH_PW.f rename sorc/ncep_post.fd/{GRIDAVG.f => UPP_MATH.f} (55%) create mode 100644 sorc/ncep_post.fd/UPP_PHYSICS.f diff --git a/sorc/ncep_post.fd/ALLOCATE_ALL.f b/sorc/ncep_post.fd/ALLOCATE_ALL.f index b8551b5ff..2bda4be0a 100644 --- a/sorc/ncep_post.fd/ALLOCATE_ALL.f +++ b/sorc/ncep_post.fd/ALLOCATE_ALL.f @@ -14,6 +14,7 @@ !! SU diagnostic fields !! - 19-07-24 Li(Kate) Zhang - Merge and update NGAC UPP for FV3-Chem !! - 19-11-23 Wen Meng - Add sea ice skin T +!! - 20-11-06 Jesse Meng - Add UPP_MATH module variables !! !! OUTPUT FILES: !! - STDOUT - RUN TIME STANDARD OUT. @@ -30,6 +31,7 @@ SUBROUTINE ALLOCATE_ALL() use vrbls2d use soil use masks + use upp_math, only: ddvdx, ddudy, uuavg ! !use params_mod use ctlblk_mod @@ -533,5 +535,9 @@ SUBROUTINE ALLOCATE_ALL() allocate(acswupt(im,jsta_2l:jend_2u)) allocate(swdnt(im,jsta_2l:jend_2u)) allocate(acswdnt(im,jsta_2l:jend_2u)) +! UPP_MATH MODULE DIFFERENTIAL EQUATIONS + allocate(ddvdx(im,jsta_2l:jend_2u)) + allocate(ddudy(im,jsta_2l:jend_2u)) + allocate(uuavg(im,jsta_2l:jend_2u)) ! end diff --git a/sorc/ncep_post.fd/BNDLYR.f b/sorc/ncep_post.fd/BNDLYR.f index 4c4ea74c7..371e241dc 100644 --- a/sorc/ncep_post.fd/BNDLYR.f +++ b/sorc/ncep_post.fd/BNDLYR.f @@ -30,6 +30,7 @@ !! 98-12-22 MIKE BALDWIN - BACK OUT RH OVER ICE !! 00-01-04 JIM TUCCILLO - MPI VERSION !! 02-01-15 MIKE BALDWIN - WRF VERSION +!! 20-11-10 JESSE MENG - USE UPP_PHYSICS MODULE !! !! USAGE: CALL BNDLYR(PBND,TBND,QBND,RHBND,UBND,VBND, !! WBND,OMGBND,PWTBND,QCNVBND) @@ -74,12 +75,12 @@ SUBROUTINE BNDLYR(PBND,TBND,QBND,RHBND,UBND,VBND, & jsta_m, jend_m, im, nbnd use physcons_post, only: con_rd, con_rv, con_eps, con_epsm1 use gridspec_mod, only: gridtype + use upp_physics, only: FPVSNEW ! implicit none ! ! DECLARE VARIABLES. ! - real,external :: FPVSNEW real,PARAMETER :: DPBND=30.E2 integer, dimension(IM,jsta:jend,NBND),intent(inout) :: LVLBND real, dimension(IM,jsta:jend,NBND),intent(inout) :: PBND,TBND, & diff --git a/sorc/ncep_post.fd/CALCAPE.f b/sorc/ncep_post.fd/CALCAPE.f deleted file mode 100644 index c1859e6c8..000000000 --- a/sorc/ncep_post.fd/CALCAPE.f +++ /dev/null @@ -1,493 +0,0 @@ -!> @file -! . . . -!> SUBPROGRAM: CALCAPE COMPUTES CAPE AND CINS -!! PRGRMMR: TREADON ORG: W/NP2 DATE: 93-02-10 -!! -!! 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 -!! - SUBROUTINE CALCAPE(ITYPE,DPBND,P1D,T1D,Q1D,L1D,CAPE, & - CINS,PPARC,ZEQL,THUND) - -! - use vrbls3d, only: pmid, t, q, zint - 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 - -! -!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none -! -! INCLUDE/SET PARAMETERS. CONSTANTS ARE FROM BOLTON (MWR, 1980). - real,PARAMETER :: ISMTHP=2,ISMTHT=2,ISMTHQ=2 -! -! DECLARE VARIABLES. -! - integer,intent(in) :: ITYPE - real,intent(in) :: DPBND - integer, dimension(IM,Jsta:jend),intent(in) :: L1D - real, dimension(IM,Jsta:jend),intent(in) :: P1D,T1D - real, dimension(IM,jsta:jend),intent(inout) :: Q1D,CAPE,CINS,PPARC,ZEQL -! - integer, dimension(im,jsta:jend) :: IEQL, IPTB, ITHTB, PARCEL, KLRES, KHRES, LCL, IDX -! - real, dimension(im,jsta:jend) :: THESP, PSP, CAPE20, QQ, PP, THUND - REAL, ALLOCATABLE :: TPAR(:,:,:) - - 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,external :: fpvsnew - integer I,J,L,KNUML,KNUMH,LBEG,LEND,IQ, KB,ITTBK - -! integer I,J,L,KNUML,KNUMH,LBEG,LEND,IQ,IT,LMHK, KB,ITTBK -! -!************************************************************** -! START CALCAPE HERE. -! - ALLOCATE(TPAR(IM,JSTA_2L:JEND_2U,LM)) -! -! 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 -! -!$omp parallel do - DO J=JSTA,JEND - DO I=1,IM - CAPE(I,J) = D00 - CAPE20(I,J) = D00 - CINS(I,J) = D00 - LCL(I,J) = 0 - THESP(I,J) = D00 - IEQL(I,J) = LM - PARCEL(I,J) = LM - PSP(I,J) = D00 - PPARC(I,J) = D00 - THUNDER(I,J) = .TRUE. - ENDDO - ENDDO -! -!$omp parallel do - DO L=1,LM - DO J=JSTA,JEND - DO I=1,IM - TPAR(I,J,L) = D00 - ENDDO - ENDDO - ENDDO -! -! TYPE 2 CAPE/CINS: -! NOTE THAT FOR TYPE 1 CAPE/CINS ARRAYS P1D, T1D, Q1D -! ARE DUMMY ARRAYS. -! - IF (ITYPE == 2) THEN -!$omp parallel do private(i,j) - DO J=JSTA,JEND - DO I=1,IM - Q1D(I,J) = MIN(MAX(H1M12,Q1D(I,J)),H99999) - ENDDO - ENDDO - ENDIF -!-------FOR ITYPE=1--FIND MAXIMUM THETA E LAYER IN LOWEST DPBND ABOVE GROUND------- -!-------FOR ITYPE=2--FIND THETA E LAYER OF GIVEN T1D, Q1D, P1D--------------------- -!--------------TRIAL MAXIMUM BUOYANCY LEVEL VARIABLES------------------- - - DO KB=1,LM -!hc IF (ITYPE==2.AND.KB>1) cycle - IF (ITYPE == 1 .OR. (ITYPE == 2 .AND. KB == 1)) THEN - -!$omp parallel do private(i,j,apebtk,apespk,bqk,bqs00k,bqs10k,iq,ittbk, & -!$omp & p00k,p01k,p10k,p11k,pkl,psfck,qbtk,sqk,sqs00k, & -!$omp & sqs10k,tbtk,tpspk,tqk,tthbtk,tthesk,tthk) - DO J=JSTA,JEND - DO I=1,IM - PSFCK = PMID(I,J,NINT(LMH(I,J))) - PKL = PMID(I,J,KB) - -!hc IF (ITYPE==1.AND.(PKLPSFCK)) cycle - IF (ITYPE ==2 .OR. & - (ITYPE == 1 .AND. (PKL >= PSFCK-DPBND .AND. PKL <= PSFCK)))THEN - IF (ITYPE == 1) THEN - TBTK = T(I,J,KB) - QBTK = max(0.0, Q(I,J,KB)) - APEBTK = (H10E5/PKL)**CAPA - ELSE - PKL = P1D(I,J) - TBTK = T1D(I,J) - QBTK = max(0.0, Q1D(I,J)) - 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(I,J) = TTHK - AINT(TTHK) - ITTBK = INT(TTHK) + 1 -!--------------KEEPING INDICES WITHIN THE TABLE------------------------- - IF(ITTBK < 1) THEN - ITTBK = 1 - QQ(I,J) = D00 - ENDIF - IF(ITTBK >= JTB) THEN - ITTBK = JTB-1 - QQ(I,J) = 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(I,J) + BQS00K - SQK = (SQS10K-SQS00K)*QQ(I,J) + SQS00K - TQK = (QBTK-BQK)/SQK*RDQ - PP(I,J) = TQK-AINT(TQK) - IQ = INT(TQK)+1 -!--------------KEEPING INDICES WITHIN THE TABLE------------------------- - IF(IQ < 1) THEN - IQ = 1 - PP(I,J) = D00 - ENDIF - IF(IQ >= ITB) THEN - IQ = ITB-1 - PP(I,J) = 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(I,J) + (P01K-P00K)*QQ(I,J) & - + (P00K-P10K-P01K+P11K)*PP(I,J)*QQ(I,J) - -!!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(I,J)) THEN - PSP (I,J) = TPSPK - THESP(I,J) = TTHESK - PARCEL(I,J) = KB - ENDIF - END IF - ENDDO ! I loop - ENDDO ! J loop - END IF - ENDDO ! KB loop - -!----FIND THE PRESSURE OF THE PARCEL THAT WAS LIFTED -!$omp parallel do private(i,j) - DO J=JSTA,JEND - DO I=1,IM - PPARC(I,J) = PMID(I,J,PARCEL(I,J)) - ENDDO - ENDDO -! -!-----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)--------------------- - DO L=1,LM -!$omp parallel do private(i,j) - DO J=JSTA,JEND - DO I=1,IM - IF (PMID(I,J,L) < PSP(I,J)) LCL(I,J) = L+1 - ENDDO - ENDDO - ENDDO -!$omp parallel do private(i,j) - DO J=JSTA,JEND - DO I=1,IM - IF (LCL(I,J) > NINT(LMH(I,J))) LCL(I,J) = NINT(LMH(I,J)) - IF (ITYPE > 2) THEN - IF (T(I,J,LCL(I,J)) < 263.15) THEN - THUNDER(I,J) = .FALSE. - ENDIF - ENDIF - ENDDO - ENDDO -!----------------------------------------------------------------------- -!---------FIND TEMP OF PARCEL LIFTED ALONG MOIST ADIABAT (TPAR)--------- -!----------------------------------------------------------------------- - - DO L=LM,1,-1 -!--------------SCALING PRESSURE & TT TABLE INDEX------------------------ - KNUML = 0 - KNUMH = 0 - DO J=JSTA,JEND - DO I=1,IM - KLRES(I,J) = 0 - KHRES(I,J) = 0 - IF(L <= LCL(I,J)) THEN - IF(PMID(I,J,L) < PLQ)THEN - KNUML = KNUML + 1 - KLRES(I,J) = 1 - ELSE - KNUMH = KNUMH + 1 - KHRES(I,J) = 1 - ENDIF - ENDIF - ENDDO - ENDDO -!*** -!*** COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE 0) THEN - CALL TTBLEX(TPAR(1,JSTA_2L,L),TTBL,ITB,JTB,KLRES & - , PMID(1,JSTA_2L,L),PL,QQ,PP,RDP,THE0,STHE & - , RDTHE,THESP,IPTB,ITHTB) - ENDIF -!*** -!*** COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE>PLQ -!** - IF(KNUMH > 0) THEN - CALL TTBLEX(TPAR(1,JSTA_2L,L),TTBLQ,ITBQ,JTBQ,KHRES & - , PMID(1,JSTA_2L,L),PLQ,QQ,PP,RDPQ & - ,THE0Q,STHEQ,RDTHEQ,THESP,IPTB,ITHTB) - ENDIF - -!------------SEARCH FOR EQ LEVEL---------------------------------------- -!$omp parallel do private(i,j) - DO J=JSTA,JEND - DO I=1,IM - IF(KHRES(I,J) > 0) THEN - IF(TPAR(I,J,L) > T(I,J,L)) IEQL(I,J) = L - ENDIF - ENDDO - ENDDO -! -!$omp parallel do private(i,j) - DO J=JSTA,JEND - DO I=1,IM - IF(KLRES(I,J) > 0) THEN - IF(TPAR(I,J,L) > T(I,J,L) .AND. & - PMID(I,J,L)>100.) IEQL(I,J) = L - ENDIF - ENDDO - ENDDO -!----------------------------------------------------------------------- - ENDDO ! end of do l=lm,1,-1 loop -!------------COMPUTE CAPE AND CINS-------------------------------------- - LBEG = 1000 - LEND = 0 - DO J=JSTA,JEND - DO I=1,IM - LBEG = MIN(IEQL(I,J),LBEG) - LEND = MAX(LCL(I,J),LEND) - ENDDO - ENDDO -! -!$omp parallel do private(i,j) - DO J=JSTA,JEND - DO I=1,IM - IF(T(I,J,IEQL(I,J)) > 255.65) THEN - THUNDER(I,J) = .FALSE. - ENDIF - ENDDO - ENDDO -! - DO L=LBEG,LEND - -!$omp parallel do private(i,j) - DO J=JSTA,JEND - DO I=1,IM - IDX(I,J) = 0 - IF(L >= IEQL(I,J).AND.L <= LCL(I,J)) THEN - IDX(I,J) = 1 - ENDIF - ENDDO - ENDDO -! -!$omp parallel do private(i,j,gdzkl,presk,thetaa,thetap,esatp,qsatp,tvp,tv) - DO J=JSTA,JEND - DO I=1,IM - IF(IDX(I,J) > 0) THEN - PRESK = PMID(I,J,L) - GDZKL = (ZINT(I,J,L)-ZINT(I,J,L+1)) * G - ESATP = min(FPVSNEW(TPAR(I,J,L)),PRESK) - QSATP = EPS*ESATP/(PRESK-ESATP*ONEPS) - TVP = TPAR(I,J,L)*(1+0.608*QSATP) - THETAP = TVP*(H10E5/PRESK)**CAPA - TV = T(I,J,L)*(1+0.608*Q(I,J,L)) - THETAA = TV*(H10E5/PRESK)**CAPA - IF(THETAP < THETAA) THEN - CINS(I,J) = CINS(I,J) + (LOG(THETAP)-LOG(THETAA))*GDZKL - ELSEIF(THETAP > THETAA) THEN - CAPE(I,J) = CAPE(I,J) + (LOG(THETAP)-LOG(THETAA))*GDZKL - IF (THUNDER(I,J) .AND. T(I,J,L) < 273.15 & - .AND. T(I,J,L) > 253.15) THEN - CAPE20(I,J) = CAPE20(I,J) + (LOG(THETAP)-LOG(THETAA))*GDZKL - ENDIF - ENDIF - ENDIF - ENDDO - ENDDO - ENDDO -! -! ENFORCE LOWER LIMIT OF 0.0 ON CAPE AND UPPER -! LIMIT OF 0.0 ON CINS. -! -!$omp parallel do private(i,j) - DO J=JSTA,JEND - DO I=1,IM - CAPE(I,J) = MAX(D00,CAPE(I,J)) - CINS(I,J) = MIN(CINS(I,J),D00) -! add equillibrium height - ZEQL(I,J) = ZINT(I,J,IEQL(I,J)) - IF (CAPE20(I,J) < 75.) THEN - THUNDER(I,J) = .FALSE. - ENDIF - IF (THUNDER(I,J)) THEN - THUND(I,J) = 1.0 - ELSE - THUND(I,J) = 0.0 - ENDIF - ENDDO - ENDDO -! - DEALLOCATE(TPAR) -! - RETURN - END diff --git a/sorc/ncep_post.fd/CALCAPE2.f b/sorc/ncep_post.fd/CALCAPE2.f deleted file mode 100644 index 36b478cd8..000000000 --- a/sorc/ncep_post.fd/CALCAPE2.f +++ /dev/null @@ -1,788 +0,0 @@ -!> @file -! . . . -!> COMPUTES CAPE AND CINS -!! PRGRMMR: TREADON ORG: W/NP2 DATE: 93-02-10 -!! -!! 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 -!! 19-09-03 J MENG - MODIFIED TO ADD 0-3KM CAPE/CINS, LFC, -!! EFFECTIVE HELICITY, DOWNDRAFT CAPE, -!! DENDRITIC GROWTH LAYER DEPTH, ESP -!! -!! USAGE: CALL CALCAPE2(ITYPE,DPBND,P1D,T1D,Q1D,L1D, & -!! CAPE,CINS,LFC,ESRHL,ESRHH, & -!! DCAPE,DGLD,ESP) -!! -!! 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) -!! LFC - LEVEL OF FREE CONVECTION (M) -!! ESRHL - LOWER BOUND TO ACCOUNT FOR EFFECTIVE HELICITY CALCULATION -!! ESRHH - UPPER BOUND TO ACCOUNT FOR EFFECTIVE HELICITY CALCULATION -!! DCAPE - DOWNDRAFT CAPE (J/KG) -!! DGLD - DENDRITIC GROWTH LAYER DEPTH (M) -!! ESP - ENHANCED STRETCHING POTENTIAL -!! -!! 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 - -!! - SUBROUTINE CALCAPE2(ITYPE,DPBND,P1D,T1D,Q1D,L1D, & - CAPE,CINS,LFC,ESRHL,ESRHH, & - DCAPE,DGLD,ESP) - use vrbls3d, only: pmid, t, q, zint - use vrbls2d, only: fis - use gridspec_mod, only: gridtype - use masks, only: lmh - use params_mod, only: d00, h1m12, h99999, h10e5, capa, elocp, eps, & - oneps, g, tfrz - 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, jm, me, jsta_m, jend_m - -! -!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none -! -! INCLUDE/SET PARAMETERS. CONSTANTS ARE FROM BOLTON (MWR, 1980). - real,PARAMETER :: ISMTHP=2,ISMTHT=2,ISMTHQ=2 -! -! DECLARE VARIABLES. -! - integer,intent(in) :: ITYPE - real,intent(in) :: DPBND - integer, dimension(IM,Jsta:jend),intent(in) :: L1D - real, dimension(IM,Jsta:jend),intent(in) :: P1D,T1D -! real, dimension(IM,jsta:jend),intent(inout) :: Q1D,CAPE,CINS,PPARC,ZEQL - real, dimension(IM,jsta:jend),intent(inout) :: Q1D,CAPE,CINS - real, dimension(IM,jsta:jend) :: PPARC,ZEQL - real, dimension(IM,jsta:jend),intent(inout) :: LFC,ESRHL,ESRHH - real, dimension(IM,jsta:jend),intent(inout) :: DCAPE,DGLD,ESP - integer, dimension(im,jsta:jend) ::L12,L17,L3KM -! - integer, dimension(im,jsta:jend) :: IEQL, IPTB, ITHTB, PARCEL, KLRES, KHRES, LCL, IDX -! - real, dimension(im,jsta:jend) :: THESP, PSP, CAPE20, QQ, PP, THUND - integer, dimension(im,jsta:jend) :: PARCEL2 - real, dimension(im,jsta:jend) :: THESP2,PSP2 - real, dimension(im,jsta:jend) :: CAPE4,CINS4 - REAL, ALLOCATABLE :: TPAR(:,:,:) - REAL, ALLOCATABLE :: TPAR2(:,:,:) - - 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 PRESK2, ESATP2, QSATP2, TVP2, THETAP2, TV2, THETAA2 - real,external :: fpvsnew - integer I,J,L,KNUML,KNUMH,LBEG,LEND,IQ, KB,ITTBK - integer IE,IW,JN,JS,IVE(JM),IVW(JM),JVN,JVS - integer ISTART,ISTOP,JSTART,JSTOP - real, dimension(IM,jsta:jend) :: HTSFC - -! integer I,J,L,KNUML,KNUMH,LBEG,LEND,IQ,IT,LMHK, KB,ITTBK -! -!************************************************************** -! START CALCAPE HERE. -! - ALLOCATE(TPAR(IM,JSTA_2L:JEND_2U,LM)) - ALLOCATE(TPAR2(IM,JSTA_2L:JEND_2U,LM)) -! -! 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 -! -!$omp parallel do - DO J=JSTA,JEND - DO I=1,IM - CAPE(I,J) = D00 - CAPE20(I,J) = D00 - CAPE4(I,J) = D00 - CINS(I,J) = D00 - CINS4(I,J) = D00 - LCL(I,J) = 0 - THESP(I,J) = D00 - IEQL(I,J) = LM - PARCEL(I,J) = LM - PSP(I,J) = D00 - PPARC(I,J) = D00 - THUNDER(I,J) = .TRUE. - LFC(I,J) = D00 - ESRHL(I,J) = D00 - ESRHH(I,J) = D00 - DCAPE(I,J) = D00 - DGLD(I,J) = D00 - ESP(I,J) = D00 - THESP2(I,J) = 500. - PSP2(I,J) = D00 - PARCEL2(I,J) = LM - ENDDO - ENDDO -! -!$omp parallel do - DO L=1,LM - DO J=JSTA,JEND - DO I=1,IM - TPAR(I,J,L) = D00 - TPAR2(I,J,L) = D00 - ENDDO - ENDDO - ENDDO -! -! FIND SURFACE HEIGHT -! - IF(gridtype == 'E')THEN - JVN = 1 - JVS = -1 - do J=JSTA,JEND - IVE(J) = MOD(J,2) - IVW(J) = IVE(J)-1 - enddo - ISTART = 2 - ISTOP = IM-1 - JSTART = JSTA_M - JSTOP = JEND_M - ELSE IF(gridtype == 'B')THEN - JVN = 1 - JVS = 0 - do J=JSTA,JEND - IVE(J)=1 - IVW(J)=0 - enddo - ISTART = 2 - ISTOP = IM-1 - JSTART = JSTA_M - JSTOP = JEND_M - ELSE - JVN = 0 - JVS = 0 - do J=JSTA,JEND - IVE(J) = 0 - IVW(J) = 0 - enddo - ISTART = 1 - ISTOP = IM - JSTART = JSTA - JSTOP = JEND - END IF -!!$omp parallel do private(htsfc,ie,iw) - IF(gridtype /= 'A') CALL EXCH(FIS(1:IM,JSTA:JEND)) - DO J=JSTART,JSTOP - DO I=ISTART,ISTOP - IE = I+IVE(J) - IW = I+IVW(J) - JN = J+JVN - JS = J+JVS -!mp PDSLVK=(PD(IW,J)*RES(IW,J)+PD(IE,J)*RES(IE,J)+ -!mp 1 PD(I,J+1)*RES(I,J+1)+PD(I,J-1)*RES(I,J-1))*0.25 -!mp PSFCK=AETA(LMV(I,J))*PDSLVK+PT - IF (gridtype=='B')THEN - HTSFC(I,J) = (0.25/g)*(FIS(IW,J)+FIS(IE,J)+FIS(I,JN)+FIS(IE,JN)) - ELSE - HTSFC(I,J) = (0.25/g)*(FIS(IW,J)+FIS(IE,J)+FIS(I,JN)+FIS(I,JS)) - ENDIF - ENDDO - ENDDO -! -! TYPE 2 CAPE/CINS: -! NOTE THAT FOR TYPE 1 CAPE/CINS ARRAYS P1D, T1D, Q1D -! ARE DUMMY ARRAYS. -! - IF (ITYPE == 2) THEN -!$omp parallel do private(i,j) - DO J=JSTA,JEND - DO I=1,IM - Q1D(I,J) = MIN(MAX(H1M12,Q1D(I,J)),H99999) - ENDDO - ENDDO - ENDIF -!-------FOR ITYPE=1--FIND MAXIMUM THETA E LAYER IN LOWEST DPBND ABOVE GROUND------- -!-------FOR ITYPE=2--FIND THETA E LAYER OF GIVEN T1D, Q1D, P1D--------------------- -!--------------TRIAL MAXIMUM BUOYANCY LEVEL VARIABLES------------------- - - DO KB=1,LM -!hc IF (ITYPE==2.AND.KB>1) cycle - IF (ITYPE == 1 .OR. (ITYPE == 2 .AND. KB == 1)) THEN - -!$omp parallel do private(i,j,apebtk,apespk,bqk,bqs00k,bqs10k,iq,ittbk, & -!$omp & p00k,p01k,p10k,p11k,pkl,psfck,qbtk,sqk,sqs00k, & -!$omp & sqs10k,tbtk,tpspk,tqk,tthbtk,tthesk,tthk) - DO J=JSTA,JEND - DO I=1,IM - PSFCK = PMID(I,J,NINT(LMH(I,J))) - PKL = PMID(I,J,KB) - -!hc IF (ITYPE==1.AND.(PKLPSFCK)) cycle - IF (ITYPE ==2 .OR. & - (ITYPE == 1 .AND. (PKL >= PSFCK-DPBND .AND. PKL <= PSFCK)))THEN - IF (ITYPE == 1) THEN - TBTK = T(I,J,KB) - QBTK = max(0.0, Q(I,J,KB)) - APEBTK = (H10E5/PKL)**CAPA - ELSE - PKL = P1D(I,J) - TBTK = T1D(I,J) - QBTK = max(0.0, Q1D(I,J)) - 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(I,J) = TTHK - AINT(TTHK) - ITTBK = INT(TTHK) + 1 -!--------------KEEPING INDICES WITHIN THE TABLE------------------------- - IF(ITTBK < 1) THEN - ITTBK = 1 - QQ(I,J) = D00 - ENDIF - IF(ITTBK >= JTB) THEN - ITTBK = JTB-1 - QQ(I,J) = 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(I,J) + BQS00K - SQK = (SQS10K-SQS00K)*QQ(I,J) + SQS00K - TQK = (QBTK-BQK)/SQK*RDQ - PP(I,J) = TQK-AINT(TQK) - IQ = INT(TQK)+1 -!--------------KEEPING INDICES WITHIN THE TABLE------------------------- - IF(IQ < 1) THEN - IQ = 1 - PP(I,J) = D00 - ENDIF - IF(IQ >= ITB) THEN - IQ = ITB-1 - PP(I,J) = 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(I,J) + (P01K-P00K)*QQ(I,J) & - + (P00K-P10K-P01K+P11K)*PP(I,J)*QQ(I,J) - -!!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(I,J)) THEN - PSP (I,J) = TPSPK - THESP(I,J) = TTHESK - PARCEL(I,J) = KB - ENDIF -!--------------CHECK FOR MINIMUM THETA E-------------------------------- - IF(TTHESK < THESP2(I,J)) THEN - PSP2 (I,J) = TPSPK - THESP2(I,J) = TTHESK - PARCEL2(I,J) = KB - ENDIF - END IF - ENDDO ! I loop - ENDDO ! J loop - END IF - ENDDO ! KB loop - -!----FIND THE PRESSURE OF THE PARCEL THAT WAS LIFTED -!$omp parallel do private(i,j) - DO J=JSTA,JEND - DO I=1,IM - PPARC(I,J) = PMID(I,J,PARCEL(I,J)) - ENDDO - ENDDO -! -!-----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)--------------------- - DO L=1,LM -!$omp parallel do private(i,j) - DO J=JSTA,JEND - DO I=1,IM - IF (PMID(I,J,L) < PSP(I,J)) LCL(I,J) = L+1 - ENDDO - ENDDO - ENDDO -!$omp parallel do private(i,j) - DO J=JSTA,JEND - DO I=1,IM - IF (LCL(I,J) > NINT(LMH(I,J))) LCL(I,J) = NINT(LMH(I,J)) - IF (ITYPE > 2) THEN - IF (T(I,J,LCL(I,J)) < 263.15) THEN - THUNDER(I,J) = .FALSE. - ENDIF - ENDIF - ENDDO - ENDDO -!----------------------------------------------------------------------- -!---------FIND TEMP OF PARCEL LIFTED ALONG MOIST ADIABAT (TPAR)--------- -!----------------------------------------------------------------------- - DO L=LM,1,-1 -!--------------SCALING PRESSURE & TT TABLE INDEX------------------------ - KNUML = 0 - KNUMH = 0 - DO J=JSTA,JEND - DO I=1,IM - KLRES(I,J) = 0 - KHRES(I,J) = 0 - IF(L <= LCL(I,J)) THEN - IF(PMID(I,J,L) < PLQ)THEN - KNUML = KNUML + 1 - KLRES(I,J) = 1 - ELSE - KNUMH = KNUMH + 1 - KHRES(I,J) = 1 - ENDIF - ENDIF - ENDDO - ENDDO -!*** -!*** COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE 0) THEN - CALL TTBLEX(TPAR(1,JSTA_2L,L),TTBL,ITB,JTB,KLRES & - , PMID(1,JSTA_2L,L),PL,QQ,PP,RDP,THE0,STHE & - , RDTHE,THESP,IPTB,ITHTB) - ENDIF -!*** -!*** COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE>PLQ -!** - IF(KNUMH > 0) THEN - CALL TTBLEX(TPAR(1,JSTA_2L,L),TTBLQ,ITBQ,JTBQ,KHRES & - , PMID(1,JSTA_2L,L),PLQ,QQ,PP,RDPQ & - ,THE0Q,STHEQ,RDTHEQ,THESP,IPTB,ITHTB) - ENDIF - -!------------SEARCH FOR EQ LEVEL---------------------------------------- -!$omp parallel do private(i,j) - DO J=JSTA,JEND - DO I=1,IM - IF(KHRES(I,J) > 0) THEN - IF(TPAR(I,J,L) > T(I,J,L)) IEQL(I,J) = L - ENDIF - ENDDO - ENDDO -! -!$omp parallel do private(i,j) - DO J=JSTA,JEND - DO I=1,IM - IF(KLRES(I,J) > 0) THEN - IF(TPAR(I,J,L) > T(I,J,L)) IEQL(I,J) = L - ENDIF - ENDDO - ENDDO -!----------------------------------------------------------------------- - ENDDO ! end of do l=lm,1,-1 loop -!------------COMPUTE CAPE AND CINS-------------------------------------- - LBEG = 1000 - LEND = 0 - DO J=JSTA,JEND - DO I=1,IM - LBEG = MIN(IEQL(I,J),LBEG) - LEND = MAX(LCL(I,J),LEND) - ENDDO - ENDDO -! -!$omp parallel do private(i,j) - DO J=JSTA,JEND - DO I=1,IM - IF(T(I,J,IEQL(I,J)) > 255.65) THEN - THUNDER(I,J) = .FALSE. - ENDIF - ENDDO - ENDDO -! -!reverse L order from bottom up for ESRH calculation -! - ESRHH = LCL - ESRHL = LCL -! DO L=LBEG,LEND - DO L=LEND,LBEG,-1 - -!$omp parallel do private(i,j) - DO J=JSTA,JEND - DO I=1,IM - IDX(I,J) = 0 - IF(L >= IEQL(I,J).AND.L <= LCL(I,J)) THEN - IDX(I,J) = 1 - ENDIF - ENDDO - ENDDO -! -!$omp parallel do private(i,j,gdzkl,presk,thetaa,thetap,esatp,qsatp,tvp,tv) - DO J=JSTA,JEND - DO I=1,IM - IF(IDX(I,J) > 0) THEN - PRESK = PMID(I,J,L) - GDZKL = (ZINT(I,J,L)-ZINT(I,J,L+1)) * G - ESATP = min(FPVSNEW(TPAR(I,J,L)),PRESK) - QSATP = EPS*ESATP/(PRESK-ESATP*ONEPS) - TVP = TPAR(I,J,L)*(1+0.608*QSATP) - THETAP = TVP*(H10E5/PRESK)**CAPA - TV = T(I,J,L)*(1+0.608*Q(I,J,L)) - THETAA = TV*(H10E5/PRESK)**CAPA - IF(THETAP < THETAA) THEN - CINS4(I,J) = CINS4(I,J) + (LOG(THETAP)-LOG(THETAA))*GDZKL - IF(ZINT(I,J,L)-HTSFC(I,J) <= 3000.) THEN - CINS(I,J) = CINS(I,J) + (LOG(THETAP)-LOG(THETAA))*GDZKL - ENDIF - ELSEIF(THETAP > THETAA) THEN - CAPE4(I,J) = CAPE4(I,J) + (LOG(THETAP)-LOG(THETAA))*GDZKL - IF(ZINT(I,J,L)-HTSFC(I,J) <= 3000.) THEN - CAPE(I,J) = CAPE(I,J) + (LOG(THETAP)-LOG(THETAA))*GDZKL - ENDIF - IF (THUNDER(I,J) .AND. T(I,J,L) < 273.15 & - .AND. T(I,J,L) > 253.15) THEN - CAPE20(I,J) = CAPE20(I,J) + (LOG(THETAP)-LOG(THETAA))*GDZKL - ENDIF - ENDIF - -! LFC - IF (ITYPE /= 1) THEN - PRESK2 = PMID(I,J,L+1) - ESATP2 = min(FPVSNEW(TPAR(I,J,L+1)),PRESK2) - QSATP2 = EPS*ESATP2/(PRESK2-ESATP2*ONEPS) - TVP2 = TPAR(I,J,L+1)*(1+0.608*QSATP2) - THETAP2 = TVP2*(H10E5/PRESK2)**CAPA - TV2 = T(I,J,L+1)*(1+0.608*Q(I,J,L+1)) - THETAA2 = TV2*(H10E5/PRESK2)**CAPA - IF(THETAP >= THETAA .AND. THETAP2 <= THETAA2) THEN - IF(LFC(I,J) == D00)THEN - LFC(I,J) = ZINT(I,J,L) - ENDIF - ENDIF - ENDIF -! -! ESRH/CAPE threshold check - IF(ZINT(I,J,L)-HTSFC(I,J) <= 3000.) THEN - IF(CAPE4(I,J) >= 100. .AND. CINS4(I,J) >= -250.) THEN - IF(ESRHL(I,J) == LCL(I,J)) ESRHL(I,J)=L - ENDIF - ESRHH(I,J)=L - ENDIF - - ENDIF !(IDX(I,J) > 0) - ENDDO - ENDDO - ENDDO - -!$omp parallel do private(i,j) - DO J=JSTA,JEND - DO I=1,IM - IF(ESRHH(I,J) > ESRHL(I,J)) ESRHH(I,J)=IEQL(I,J) - ENDDO - ENDDO -! -! ENFORCE LOWER LIMIT OF 0.0 ON CAPE AND UPPER -! LIMIT OF 0.0 ON CINS. -! ENFORCE LFC ABOVE LCL AND BELOW EL -! -!$omp parallel do private(i,j) - DO J=JSTA,JEND - DO I=1,IM - CAPE(I,J) = MAX(D00,CAPE(I,J)) - CINS(I,J) = MIN(CINS(I,J),D00) -! equillibrium height - ZEQL(I,J) = ZINT(I,J,IEQL(I,J)) - LFC(I,J) = MIN(LFC(I,J),ZINT(I,J,IEQL(I,J))) - LFC(I,J) = MAX(ZINT(I,J, LCL(I,J)),LFC(I,J)) - IF (CAPE20(I,J) < 75.) THEN - THUNDER(I,J) = .FALSE. - ENDIF - IF (THUNDER(I,J)) THEN - THUND(I,J) = 1.0 - ELSE - THUND(I,J) = 0.0 - ENDIF - ENDDO - ENDDO -!------------COMPUTE DCAPE-------------------------------------- -!----------------------------------------------------------------------- -!-----CHOOSE LAYER DIRECTLY BELOW PSP2 AS LCL AND------------------------ - DO L=1,LM -!$omp parallel do private(i,j) - DO J=JSTA,JEND - DO I=1,IM - IF (PMID(I,J,L) < PSP2(I,J)) LCL(I,J) = L+1 - ENDDO - ENDDO - ENDDO -!$omp parallel do private(i,j) - DO J=JSTA,JEND - DO I=1,IM - IF (LCL(I,J) > NINT(LMH(I,J))) LCL(I,J) = NINT(LMH(I,J)) - ENDDO - ENDDO - -!---------FIND TEMP OF PARCEL DESCENDED ALONG MOIST ADIABAT (TPAR)--------- -!----------------------------------------------------------------------- - - DO L=LM,1,-1 -!--------------SCALING PRESSURE & TT TABLE INDEX------------------------ - KNUML = 0 - KNUMH = 0 - DO J=JSTA,JEND - DO I=1,IM - KLRES(I,J) = 0 - KHRES(I,J) = 0 - ! IF(L <= LCL(I,J)) THEN - IF(PMID(I,J,L) < PLQ)THEN - KNUML = KNUML + 1 - KLRES(I,J) = 1 - ELSE - KNUMH = KNUMH + 1 - KHRES(I,J) = 1 - ENDIF - ! ENDIF - PSFCK = PMID(I,J,NINT(LMH(I,J))) - PKL = PMID(I,J,L) - IF(PKL >= PSFCK-DPBND) PARCEL2(I,J)=L - ENDDO - ENDDO -!*** -!*** COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE 0) THEN - CALL TTBLEX(TPAR2(1,JSTA_2L,L),TTBL,ITB,JTB,KLRES & - , PMID(1,JSTA_2L,L),PL,QQ,PP,RDP,THE0,STHE & - , RDTHE,THESP2,IPTB,ITHTB) - ENDIF -!*** -!*** COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE>PLQ -!** - IF(KNUMH > 0) THEN - CALL TTBLEX(TPAR2(1,JSTA_2L,L),TTBLQ,ITBQ,JTBQ,KHRES & - , PMID(1,JSTA_2L,L),PLQ,QQ,PP,RDPQ & - , THE0Q,STHEQ,RDTHEQ,THESP2,IPTB,ITHTB) - ENDIF - ENDDO ! end of do l=lm,1,-1 loop - - LBEG = LM - LEND = LM -!$omp parallel do private(i,j) - DO J=JSTA,JEND - DO I=1,IM - ! LBEG = MIN(PARCEL2(I,J),LBEG) - LBEG = MIN(LCL(I,J),LBEG) - ENDDO - ENDDO - - DO L=LBEG,LM -!$omp parallel do private(i,j) - DO J=JSTA,JEND - DO I=1,IM - IDX(I,J) = 0 - ! IF(L >= PARCEL2(I,J)) THEN - IF(L >= LCL(I,J)) THEN - IDX(I,J) = 1 - ENDIF - ENDDO - ENDDO -! -!$omp parallel do private(i,j,gdzkl,presk,thetaa,thetap,esatp,qsatp,tvp,tv) - DO J=JSTA,JEND - DO I=1,IM - IF(IDX(I,J) > 0) THEN - PRESK = PMID(I,J,L) - GDZKL = (ZINT(I,J,L)-ZINT(I,J,L+1)) * G - ESATP = min(FPVSNEW(TPAR2(I,J,L)),PRESK) - QSATP = EPS*ESATP/(PRESK-ESATP*ONEPS) - TVP = TPAR2(I,J,L)*(1+0.608*QSATP) - THETAP = TVP*(H10E5/PRESK)**CAPA - TV = T(I,J,L)*(1+0.608*Q(I,J,L)) - THETAA = TV*(H10E5/PRESK)**CAPA - !IF(THETAP > THETAA) THEN - IF(THETAP < THETAA) THEN - DCAPE(I,J) = DCAPE(I,J) + (LOG(THETAP)-LOG(THETAA))*GDZKL - ENDIF - ENDIF - ENDDO - ENDDO - ENDDO - -!$omp parallel do private(i,j) - DO J=JSTA,JEND - DO I=1,IM - !DCAPE(I,J) = MAX(D00,DCAPE(I,J)) - DCAPE(I,J) = MIN(D00,DCAPE(I,J)) - ENDDO - ENDDO -! -! Dendritic Growth Layer depth -! the layer with temperatures from -12 to -17 C in meters -! - L12=LM - L17=LM - DO L=LM,1,-1 -!$omp parallel do private(i,j) - DO J=JSTA,JEND - DO I=1,IM - IF(T(I,J,L) <= TFRZ-12. .AND. L12(I,J)==LM) L12(I,J)=L - IF(T(I,J,L) <= TFRZ-17. .AND. L17(I,J)==LM) L17(I,J)=L - ENDDO - ENDDO - ENDDO -!$omp parallel do private(i,j) - DO J=JSTA,JEND - DO I=1,IM - IF(L12(I,J)/=LM .AND. L17(I,J)/=LM) THEN - DGLD(I,J)=ZINT(I,J,L17(I,J))-ZINT(I,J,L12(I,J)) - DGLD(I,J)=MAX(DGLD(I,J),0.) - ENDIF - ENDDO - ENDDO -! -! Enhanced Stretching Potential -! ESP = (0-3 km MLCAPE / 50 J kg-1) * ((0-3 km lapse rate - 7.0) / 1.0 C (km-1) -! https://www.spc.noaa.gov/exper/soundings/help/params4.html -! - L3KM=LM - DO L=LM,1,-1 -!$omp parallel do private(i,j) - DO J=JSTA,JEND - DO I=1,IM - IF(ZINT(I,J,L)-HTSFC(I,J) <= 3000.) L3KM(I,J)=L - ENDDO - ENDDO - ENDDO -!$omp parallel do private(i,j) - DO J=JSTA,JEND - DO I=1,IM - ESP(I,J) = (CAPE(I,J) / 50.) * (T(I,J,LM) - T(I,J,L3KM(I,J)) - 7.0) - IF((T(I,J,LM) - T(I,J,L3KM(I,J))) < 7.0) ESP(I,J) = 0. -! IF(CAPE(I,J) < 250.) ESP(I,J) = 0. - ENDDO - ENDDO -! - DEALLOCATE(TPAR) - DEALLOCATE(TPAR2) -! - RETURN - END diff --git a/sorc/ncep_post.fd/CALPW.f b/sorc/ncep_post.fd/CALPW.f index c880080d4..c16619bed 100644 --- a/sorc/ncep_post.fd/CALPW.f +++ b/sorc/ncep_post.fd/CALPW.f @@ -33,6 +33,7 @@ !! 15-07-04 SARAH LU - CORRECT PW INTEGRATION FOR AOD (17) !! 15-07-10 SARAH LU - UPDATE TO CALCULATE ASYMETRY PARAMETER !! 19-07-25 Li(Kate) Zhang - MERGE SARHA LU's update for FV3-Chem +!! 20-11-10 JESSE MENG - USE UPP_PHYSICS MODULE !! !! USAGE: CALL CALPW(PW) !! INPUT ARGUMENT LIST: @@ -65,6 +66,7 @@ SUBROUTINE CALPW(PW,IDECID) use masks, only: htm use params_mod, only: tfrz, gi use ctlblk_mod, only: lm, jsta, jend, im + use upp_physics, only: FPVSNEW !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none ! @@ -83,7 +85,6 @@ SUBROUTINE CALPW(PW,IDECID) real,dimension(IM,jsta:jend),intent(inout) :: PW INTEGER LLMH,I,J,L REAL ALPM,DZ,PM,PWSUM,RHOAIR,DP,ES - real,external :: FPVSNEW REAL QDUM(IM,jsta:jend), PWS(IM,jsta:jend),QS(IM,jsta:jend) ! !*************************************************************** diff --git a/sorc/ncep_post.fd/CALRH.f b/sorc/ncep_post.fd/CALRH.f deleted file mode 100644 index 057f62524..000000000 --- a/sorc/ncep_post.fd/CALRH.f +++ /dev/null @@ -1,94 +0,0 @@ -!> @file -! -!> SUBPROGRAM: CALRH COMPUTES RELATIVE HUMIDITY -!! PRGRMMR: TREADON ORG: W/NP2 DATE: 92-12-22 -!! -!! ABSTRACT: -!! THIS ROUTINE COMPUTES RELATIVE HUMIDITY GIVEN PRESSURE, -!! TEMPERATURE, SPECIFIC HUMIDITY. AN UPPER AND LOWER BOUND -!! OF 100 AND 1 PERCENT RELATIVE HUMIDITY IS ENFORCED. WHEN -!! THESE BOUNDS ARE APPLIED THE PASSED SPECIFIC HUMIDITY -!! ARRAY IS ADJUSTED AS NECESSARY TO PRODUCE THE SET RELATIVE -!! HUMIDITY. -!! -!! PROGRAM HISTORY LOG: -!! ??-??-?? DENNIS DEAVEN -!! 92-12-22 RUSS TREADON - MODIFIED AS DESCRIBED ABOVE. -!! 98-06-08 T BLACK - CONVERSION FROM 1-D TO 2-D -!! 98-08-18 MIKE BALDWIN - MODIFY TO COMPUTE RH OVER ICE AS IN MODEL -!! 98-12-16 GEOFF MANIKIN - UNDO RH COMPUTATION OVER ICE -!! 00-01-04 JIM TUCCILLO - MPI VERSION -!! 02-06-11 MIKE BALDWIN - WRF VERSION -!! 06-03-19 Wen Meng - MODIFY TOP PRESSURE to 1 PA -!! -!! USAGE: CALL CALRH(P1,T1,Q1,RH) -!! INPUT ARGUMENT LIST: -!! P1 - PRESSURE (PA) -!! T1 - TEMPERATURE (K) -!! Q1 - SPECIFIC HUMIDITY (KG/KG) -!! -!! OUTPUT ARGUMENT LIST: -!! RH - RELATIVE HUMIDITY (DECIMAL FORM) -!! Q1 - ADJUSTED SPECIFIC HUMIDITY (KG/KG) -!! -!! OUTPUT FILES: -!! NONE -!! -!! SUBPROGRAMS CALLED: -!! UTILITIES: -!! LIBRARY: -!! NONE -!! -!! ATTRIBUTES: -!! LANGUAGE: FORTRAN -!! MACHINE : CRAY C-90 -!! - SUBROUTINE CALRH(P1,T1,Q1,RH) - -! - use params_mod, only: PQ0, a2, a3, a4, rhmin - use ctlblk_mod, only: jsta, jend, spval, im -!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none -! -! SET PARAMETER. -! -! DECLARE VARIABLES. -! - REAL,dimension(IM,jsta:jend),intent(in) :: P1,T1 - REAL,dimension(IM,jsta:jend),intent(inout) :: Q1 - REAL,dimension(IM,jsta:jend),intent(out) :: RH - REAL QC - integer I,J -!*************************************************************** -! -! START CALRH. -! - DO J=JSTA,JEND - DO I=1,IM - 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)) -! - RH(I,J) = Q1(I,J)/QC -! -! BOUNDS CHECK -! - IF (RH(I,J) > 1.0) THEN - RH(I,J) = 1.0 - Q1(I,J) = RH(I,J)*QC - ENDIF - IF (RH(I,J) < RHmin) THEN !use smaller RH limit for stratosphere - RH(I,J) = RHmin - Q1(I,J) = RH(I,J)*QC - ENDIF -! - ENDIF - ELSE - RH(I,J) = SPVAL - ENDIF - ENDDO - ENDDO - - RETURN - END diff --git a/sorc/ncep_post.fd/CALRH_GFS.f b/sorc/ncep_post.fd/CALRH_GFS.f deleted file mode 100644 index b6ff5da40..000000000 --- a/sorc/ncep_post.fd/CALRH_GFS.f +++ /dev/null @@ -1,209 +0,0 @@ -!> @file -! -!> SUBPROGRAM: CALRH COMPUTES RELATIVE HUMIDITY -!! PRGRMMR: TREADON ORG: W/NP2 DATE: 92-12-22 -!! -!! ABSTRACT: -!! THIS ROUTINE COMPUTES RELATIVE HUMIDITY GIVEN PRESSURE, -!! TEMPERATURE, SPECIFIC HUMIDITY. AN UPPER AND LOWER BOUND -!! OF 100 AND 1 PERCENT RELATIVE HUMIDITY IS ENFORCED. WHEN -!! THESE BOUNDS ARE APPLIED THE PASSED SPECIFIC HUMIDITY -!! ARRAY IS ADJUSTED AS NECESSARY TO PRODUCE THE SET RELATIVE -!! HUMIDITY. -!! -!! PROGRAM HISTORY LOG: -!! ??-??-?? DENNIS DEAVEN -!! 92-12-22 RUSS TREADON - MODIFIED AS DESCRIBED ABOVE. -!! 98-06-08 T BLACK - CONVERSION FROM 1-D TO 2-D -!! 98-08-18 MIKE BALDWIN - MODIFY TO COMPUTE RH OVER ICE AS IN MODEL -!! 98-12-16 GEOFF MANIKIN - UNDO RH COMPUTATION OVER ICE -!! 00-01-04 JIM TUCCILLO - MPI VERSION -!! 02-06-11 MIKE BALDWIN - WRF VERSION -!! 13-08-13 S. Moorthi - Threading -!! 06-03-19 Wen Meng - MODIFY TOP PRESSURE to 1 PA -!! -!! USAGE: CALL CALRH(P1,T1,Q1,RH) -!! INPUT ARGUMENT LIST: -!! P1 - PRESSURE (PA) -!! T1 - TEMPERATURE (K) -!! Q1 - SPECIFIC HUMIDITY (KG/KG) -!! -!! OUTPUT ARGUMENT LIST: -!! RH - RELATIVE HUMIDITY (DECIMAL FORM) -!! Q1 - ADJUSTED SPECIFIC HUMIDITY (KG/KG) -!! -!! OUTPUT FILES: -!! NONE -!! -!! SUBPROGRAMS CALLED: -!! UTILITIES: -!! LIBRARY: -!! NONE -!! -!! ATTRIBUTES: -!! LANGUAGE: FORTRAN -!! MACHINE : CRAY C-90 -!! - SUBROUTINE CALRH_GFS(P1,T1,Q1,RH) - -! - use params_mod, only: rhmin - use ctlblk_mod, only: jsta, jend, spval, im -!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none -! - real,parameter:: con_rd =2.8705e+2 ! gas constant air (J/kg/K) - real,parameter:: con_rv =4.6150e+2 ! gas constant H2O - real,parameter:: con_eps =con_rd/con_rv - real,parameter:: con_epsm1 =con_rd/con_rv-1 -! real,external::FPVSNEW - - INTERFACE - ELEMENTAL FUNCTION FPVSNEW (t) - REAL FPVSNEW - REAL, INTENT(IN) :: t - END FUNCTION FPVSNEW - END INTERFACE -! - REAL,dimension(IM,jsta:jend),intent(in) :: P1,T1 - REAL,dimension(IM,jsta:jend),intent(inout):: Q1,RH - REAL ES,QC - integer :: I,J -!*************************************************************** -! -! START CALRH. -! -!$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 (ABS(P1(I,J)) > 1.0) THEN -! IF (P1(I,J) > 1.0) THEN - IF (P1(I,J) >= 1.0) THEN - ES = MIN(FPVSNEW(T1(I,J)),P1(I,J)) - QC = CON_EPS*ES/(P1(I,J)+CON_EPSM1*ES) - -! QC=PQ0/P1(I,J)*EXP(A2*(T1(I,J)-A3)/(T1(I,J)-A4)) - - RH(I,J) = min(1.0,max(Q1(I,J)/QC,rhmin)) - q1(i,j) = rh(i,j)*qc - -! BOUNDS CHECK -! -! IF (RH(I,J) > 1.0) THEN -! RH(I,J) = 1.0 -! Q1(I,J) = RH(I,J)*QC -! ELSEIF (RH(I,J) < RHmin) THEN !use smaller RH limit for stratosphere -! RH(I,J) = RHmin -! Q1(I,J) = RH(I,J)*QC -! ENDIF - - ENDIF - ELSE - RH(I,J) = SPVAL - ENDIF - ENDDO - ENDDO - - RETURN - END - - elemental function fpvsnew(t) -!$$$ Subprogram Documentation Block -! -! Subprogram: fpvsnew Compute saturation vapor pressure -! Author: N Phillips w/NMC2X2 Date: 30 dec 82 -! -! Abstract: Compute saturation vapor pressure from the temperature. -! A linear interpolation is done between values in a lookup table -! computed in gpvs. See documentation for fpvsx for details. -! Input values outside table range are reset to table extrema. -! The interpolation accuracy is almost 6 decimal places. -! On the Cray, fpvs is about 4 times faster than exact calculation. -! This function should be expanded inline in the calling routine. -! -! Program History Log: -! 91-05-07 Iredell made into inlinable function -! 94-12-30 Iredell expand table -! 1999-03-01 Iredell f90 module -! 2001-02-26 Iredell ice phase -! -! Usage: pvs=fpvsnew(t) -! -! Input argument list: -! t Real(krealfp) temperature in Kelvin -! -! Output argument list: -! fpvsnew Real(krealfp) saturation vapor pressure in Pascals -! -! Attributes: -! Language: Fortran 90. -! -!$$$ - implicit none - integer,parameter:: nxpvs=7501 - real,parameter:: con_ttp =2.7316e+2 ! temp at H2O 3pt - real,parameter:: con_psat =6.1078e+2 ! pres at H2O 3pt - real,parameter:: con_cvap =1.8460e+3 ! spec heat H2O gas (J/kg/K) - real,parameter:: con_cliq =4.1855e+3 ! spec heat H2O liq - real,parameter:: con_hvap =2.5000e+6 ! lat heat H2O cond - real,parameter:: con_rv =4.6150e+2 ! gas constant H2O - real,parameter:: con_csol =2.1060e+3 ! spec heat H2O ice - real,parameter:: con_hfus =3.3358e+5 ! lat heat H2O fusion - real,parameter:: tliq=con_ttp - real,parameter:: tice=con_ttp-20.0 - real,parameter:: dldtl=con_cvap-con_cliq - real,parameter:: heatl=con_hvap - real,parameter:: xponal=-dldtl/con_rv - real,parameter:: xponbl=-dldtl/con_rv+heatl/(con_rv*con_ttp) - real,parameter:: dldti=con_cvap-con_csol - real,parameter:: heati=con_hvap+con_hfus - real,parameter:: xponai=-dldti/con_rv - real,parameter:: xponbi=-dldti/con_rv+heati/(con_rv*con_ttp) - real tr,w,pvl,pvi - real fpvsnew - real,intent(in):: t - integer jx - real xj,x,tbpvs(nxpvs),xp1 - real xmin,xmax,xinc,c2xpvs,c1xpvs -! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - xmin=180.0 - xmax=330.0 - xinc=(xmax-xmin)/(nxpvs-1) -! c1xpvs=1.-xmin/xinc - c2xpvs=1./xinc - c1xpvs=1.-xmin*c2xpvs -! xj=min(max(c1xpvs+c2xpvs*t,1.0),real(nxpvs,krealfp)) - xj=min(max(c1xpvs+c2xpvs*t,1.0),float(nxpvs)) - jx=min(xj,float(nxpvs)-1.0) - x=xmin+(jx-1)*xinc - - tr=con_ttp/x - if(x>=tliq) then - tbpvs(jx)=con_psat*(tr**xponal)*exp(xponbl*(1.-tr)) - elseif(x=tliq) then - tbpvs(jx+1)=con_psat*(tr**xponal)*exp(xponbl*(1.-tr)) - elseif(xp1AERFD CAN BE PROCESSED WHEN NIN=NBIN_DU +! 20-11-10 JESSE MENG - USE UPP_PHYSICS MODULE ! ! USAGE: CALL FDLVL_MASS(ITYPE,NFD,PTFD,HTFD,NIN,QIN,QTYPE,QFD) ! INPUT ARGUMENT LIST: @@ -902,6 +903,7 @@ SUBROUTINE FDLVL_MASS(ITYPE,NFD,PTFD,HTFD,NIN,QIN,QTYPE,QFD) JEND_M, IM, JM,global,MODELNAME use gridspec_mod, only: GRIDTYPE use physcons_post,only: CON_FVIRT, CON_ROG, CON_EPS, CON_EPSM1 + use upp_physics, only: FPVSNEW !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none ! @@ -910,7 +912,6 @@ SUBROUTINE FDLVL_MASS(ITYPE,NFD,PTFD,HTFD,NIN,QIN,QTYPE,QFD) ! DECLARE VARIABLES ! real,parameter:: zshul=75.,tvshul=290.66 - real,external :: fpvsnew integer,intent(in) :: ITYPE(NFD) integer,intent(in) :: NFD ! coming from calling subroutine diff --git a/sorc/ncep_post.fd/FRZLVL.f b/sorc/ncep_post.fd/FRZLVL.f index d18bcbb3d..d84bd7da4 100644 --- a/sorc/ncep_post.fd/FRZLVL.f +++ b/sorc/ncep_post.fd/FRZLVL.f @@ -33,6 +33,7 @@ !! 01-10-25 H CHUANG - MODIFIED TO PROCESS HYBRID MODEL OUTPUT !! 02-01-15 MIKE BALDWIN - WRF VERSION !! 19-10-30 Bo CUI - REMOVE "GOTO" STATEMENT +!! 20-11-10 JESSE MENG - USE UPP_PHYSICS MODULE !! !! USAGE: CALL FRZLVL(ZFRZ,RHFRZ) !! INPUT ARGUMENT LIST: @@ -71,12 +72,9 @@ SUBROUTINE FRZLVL(ZFRZ,RHFRZ,PFRZL) use params_mod, only: gi, d00, capa, d0065, tfrz, pq0, a2, a3, a4 use ctlblk_mod, only: jsta, jend, spval, lm, modelname, im use physcons_post, only: con_rd, con_rv, con_eps, con_epsm1 + use upp_physics, only: FPVSNEW implicit none - - real,external::FPVSNEW -! -! implicit none ! ! DECLARE VARIABLES. ! diff --git a/sorc/ncep_post.fd/FRZLVL2.f b/sorc/ncep_post.fd/FRZLVL2.f index 795fe2328..ff1d3ccd7 100644 --- a/sorc/ncep_post.fd/FRZLVL2.f +++ b/sorc/ncep_post.fd/FRZLVL2.f @@ -36,6 +36,7 @@ !! 02-01-15 MIKE BALDWIN - WRF VERSION !! 10-08-27 T. Smirnova - added PFRZL to the output !! 16-01-21 C. Alexander - Generalized function for any isotherm +!! 20-11-10 JESSE MENG - USE UPP_PHYSICS MODULE !! !! USAGE: CALL FRZLVL2(ISOTHERM,ZFRZ,RHFRZ,PFRZL) !! INPUT ARGUMENT LIST: @@ -68,10 +69,10 @@ SUBROUTINE FRZLVL2(ISOTHERM,ZFRZ,RHFRZ,PFRZL) use params_mod, only: gi, d00, capa, d0065, tfrz, pq0, a2, a3, a4, d50 use ctlblk_mod, only: jsta, jend, spval, lm, modelname, im use physcons_post, only: con_rd, con_rv, con_eps, con_epsm1 + use upp_physics, only: FPVSNEW implicit none - real,external::FPVSNEW !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! implicit none ! diff --git a/sorc/ncep_post.fd/INITPOST_GFS.f b/sorc/ncep_post.fd/INITPOST_GFS.f index 966e79fcb..21f91b5b4 100644 --- a/sorc/ncep_post.fd/INITPOST_GFS.f +++ b/sorc/ncep_post.fd/INITPOST_GFS.f @@ -65,6 +65,7 @@ SUBROUTINE INITPOST_GFS(NREC,iunit,iostatusFlux,iunitd3d,iostatusD3D,gfile) truelat2, psmapf, cenlat, truelat1 use rqstfld_mod, only: igds, iq, is, avbl use sfcio_module, only: sfcio_head, sfcio_data, sfcio_srohdc + use upp_physics, only: fpvsnew !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none ! @@ -84,7 +85,6 @@ SUBROUTINE INITPOST_GFS(NREC,iunit,iostatusFlux,iunitd3d,iostatusD3D,gfile) ! real,parameter:: con_eps =con_rd/con_rv ! real,parameter:: con_epsm1 =con_rd/con_rv-1 ! -! real,external::FPVSNEW ! This version of INITPOST shows how to initialize, open, read from, and ! close a NetCDF dataset. In order to change it to read an internal (binary) ! dataset, do a global replacement of _ncd_ with _int_. @@ -130,7 +130,6 @@ SUBROUTINE INITPOST_GFS(NREC,iunit,iostatusFlux,iunitd3d,iostatusD3D,gfile) integer ii,jj,js,je,jev,iyear,imn,iday,itmp,ioutcount,istatus, & I,J,L,ll,k,kf,irtn,igdout,n,Index real TSTART,TLMH,TSPH,ES, FACT,soilayert,soilayerb - real, external :: fpvsnew character*8,allocatable:: recname(:) character*16,allocatable :: reclevtyp(:) diff --git a/sorc/ncep_post.fd/INITPOST_GFS_NEMS.f b/sorc/ncep_post.fd/INITPOST_GFS_NEMS.f index 86d5bbffe..13fb5442f 100644 --- a/sorc/ncep_post.fd/INITPOST_GFS_NEMS.f +++ b/sorc/ncep_post.fd/INITPOST_GFS_NEMS.f @@ -83,6 +83,7 @@ SUBROUTINE INITPOST_GFS_NEMS(NREC,iostatusFlux,iostatusD3D, & use gridspec_mod, only: maptype, gridtype, latstart, latlast, lonstart, lonlast, cenlon, & dxval, dyval, truelat2, truelat1, psmapf, cenlat use rqstfld_mod, only: igds, avbl, iq, is + use upp_physics, only: fpvsnew ! use wrf_io_flags_mod, only: ! Do we need this? !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none @@ -102,7 +103,6 @@ SUBROUTINE INITPOST_GFS_NEMS(NREC,iostatusFlux,iostatusD3D, & ! real,parameter:: con_eps =con_rd/con_rv ! real,parameter:: con_epsm1 =con_rd/con_rv-1 ! -! real,external::FPVSNEW ! This version of INITPOST shows how to initialize, open, read from, and ! close a NetCDF dataset. In order to change it to read an internal (binary) ! dataset, do a global replacement of _ncd_ with _int_. @@ -144,7 +144,6 @@ SUBROUTINE INITPOST_GFS_NEMS(NREC,iostatusFlux,iostatusD3D, & impf,jmpf,nframed2,iunitd3d,ierr,idum,iret real TSTART,TLMH,TSPH,ES,FACT,soilayert,soilayerb,zhour,dum, & tvll,pmll,tv - real, external :: fpvsnew character*8, allocatable :: recname(:) character*16,allocatable :: reclevtyp(:) diff --git a/sorc/ncep_post.fd/INITPOST_GFS_NEMS_MPIIO.f b/sorc/ncep_post.fd/INITPOST_GFS_NEMS_MPIIO.f index 7685e8b63..0658d3683 100644 --- a/sorc/ncep_post.fd/INITPOST_GFS_NEMS_MPIIO.f +++ b/sorc/ncep_post.fd/INITPOST_GFS_NEMS_MPIIO.f @@ -96,6 +96,7 @@ SUBROUTINE INITPOST_GFS_NEMS_MPIIO(iostatusAER) dxval, dyval, truelat2, truelat1, psmapf, cenlat use rqstfld_mod, only: igds, avbl, iq, is use nemsio_module_mpi + use upp_physics, only: fpvsnew !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none ! @@ -114,7 +115,6 @@ SUBROUTINE INITPOST_GFS_NEMS_MPIIO(iostatusAER) ! real,parameter:: con_eps =con_rd/con_rv ! real,parameter:: con_epsm1 =con_rd/con_rv-1 ! -! real,external::FPVSNEW ! This version of INITPOST shows how to initialize, open, read from, and ! close a NetCDF dataset. In order to change it to read an internal (binary) ! dataset, do a global replacement of _ncd_ with _int_. @@ -160,7 +160,6 @@ SUBROUTINE INITPOST_GFS_NEMS_MPIIO(iostatusAER) impf,jmpf,nframed2,iunitd3d,ierr,idum,iret,nrec,idrt real TSTART,TLMH,TSPH,ES,FACT,soilayert,soilayerb,zhour,dum, & tvll,pmll,tv, tx1, tx2 - real, external :: fpvsnew character*16,allocatable :: recname(:) character*16,allocatable :: reclevtyp(:) diff --git a/sorc/ncep_post.fd/INITPOST_GFS_NETCDF.f b/sorc/ncep_post.fd/INITPOST_GFS_NETCDF.f index 3dcc4b2b9..ea1656250 100644 --- a/sorc/ncep_post.fd/INITPOST_GFS_NETCDF.f +++ b/sorc/ncep_post.fd/INITPOST_GFS_NETCDF.f @@ -81,6 +81,7 @@ SUBROUTINE INITPOST_GFS_NETCDF(ncid3d) dxval, dyval, truelat2, truelat1, psmapf, cenlat,lonstartv, lonlastv, cenlonv, & latstartv, latlastv, cenlatv,latstart_r,latlast_r,lonstart_r,lonlast_r use rqstfld_mod, only: igds, avbl, iq, is + use upp_physics, only: fpvsnew !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none ! @@ -104,7 +105,6 @@ SUBROUTINE INITPOST_GFS_NETCDF(ncid3d) ! real,parameter:: con_eps =con_rd/con_rv ! real,parameter:: con_epsm1 =con_rd/con_rv-1 ! -! real,external::FPVSNEW ! This version of INITPOST shows how to initialize, open, read from, and ! close a NetCDF dataset. In order to change it to read an internal (binary) ! dataset, do a global replacement of _ncd_ with _int_. @@ -149,7 +149,6 @@ SUBROUTINE INITPOST_GFS_NETCDF(ncid3d) integer ncid3d,ncid2d,varid,nhcas real TSTART,TLMH,TSPH,ES,FACT,soilayert,soilayerb,zhour,dum, & tvll,pmll,tv, tx1, tx2 - real, external :: fpvsnew character*20,allocatable :: recname(:) integer, allocatable :: reclev(:), kmsk(:,:) diff --git a/sorc/ncep_post.fd/INITPOST_GFS_NETCDF_PARA.f b/sorc/ncep_post.fd/INITPOST_GFS_NETCDF_PARA.f index 91d7b8d5d..5f7b9b864 100644 --- a/sorc/ncep_post.fd/INITPOST_GFS_NETCDF_PARA.f +++ b/sorc/ncep_post.fd/INITPOST_GFS_NETCDF_PARA.f @@ -81,6 +81,7 @@ SUBROUTINE INITPOST_GFS_NETCDF_PARA(ncid3d) dxval, dyval, truelat2, truelat1, psmapf, cenlat,lonstartv, lonlastv, cenlonv, & latstartv, latlastv, cenlatv,latstart_r,latlast_r,lonstart_r,lonlast_r use rqstfld_mod, only: igds, avbl, iq, is + use upp_physics, only: fpvsnew !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none ! @@ -104,7 +105,6 @@ SUBROUTINE INITPOST_GFS_NETCDF_PARA(ncid3d) ! real,parameter:: con_eps =con_rd/con_rv ! real,parameter:: con_epsm1 =con_rd/con_rv-1 ! -! real,external::FPVSNEW ! This version of INITPOST shows how to initialize, open, read from, and ! close a NetCDF dataset. In order to change it to read an internal (binary) ! dataset, do a global replacement of _ncd_ with _int_. @@ -149,7 +149,6 @@ SUBROUTINE INITPOST_GFS_NETCDF_PARA(ncid3d) integer ncid3d,ncid2d,varid,nhcas real TSTART,TLMH,TSPH,ES,FACT,soilayert,soilayerb,zhour,dum, & tvll,pmll,tv, tx1, tx2 - real, external :: fpvsnew character*20,allocatable :: recname(:) integer, allocatable :: reclev(:), kmsk(:,:) diff --git a/sorc/ncep_post.fd/INITPOST_GFS_SIGIO.f b/sorc/ncep_post.fd/INITPOST_GFS_SIGIO.f index bee74b5b9..36aeb74bb 100644 --- a/sorc/ncep_post.fd/INITPOST_GFS_SIGIO.f +++ b/sorc/ncep_post.fd/INITPOST_GFS_SIGIO.f @@ -88,6 +88,7 @@ SUBROUTINE INITPOST_GFS_SIGIO(lusig,iunit,iostatusFlux,iostatusD3D,idrt,sighead) use rqstfld_mod, only: IGDS, AVBL, IQ, IS use sigio_module, only: SIGIO_HEAD use sfcio_module, only: sfcio_head, sfcio_data, sfcio_srohdc + use upp_physics, only: fpvsnew ! use wrf_io_flags_mod !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none @@ -109,7 +110,6 @@ SUBROUTINE INITPOST_GFS_SIGIO(lusig,iunit,iostatusFlux,iostatusD3D,idrt,sighead) ! real,parameter:: con_eps =con_rd/con_rv ! real,parameter:: con_epsm1 =con_rd/con_rv-1 ! -! real,external::FPVSNEW ! This version of INITPOST shows how to initialize, open, read from, and ! close a NetCDF dataset. In order to change it to read an internal (binary) ! dataset, do a global replacement of _ncd_ with _int_. @@ -169,7 +169,6 @@ SUBROUTINE INITPOST_GFS_SIGIO(lusig,iunit,iostatusFlux,iostatusD3D,idrt,sighead) I,J,L,ll,k,kf,irtn,igdout,n,Index,nframe, & impf,jmpf,nframed2,iunitd3d real TSTART,TLMH,TSPH,ES, FACT,soilayert,soilayerb,zhour,dum - real, external :: fpvsnew real, allocatable:: glat1d(:),glon1d(:),qstl(:) integer ierr,idum diff --git a/sorc/ncep_post.fd/INITPOST_NEMS.f b/sorc/ncep_post.fd/INITPOST_NEMS.f index 892ec1129..701272f6c 100644 --- a/sorc/ncep_post.fd/INITPOST_NEMS.f +++ b/sorc/ncep_post.fd/INITPOST_NEMS.f @@ -66,6 +66,7 @@ SUBROUTINE INITPOST_NEMS(NREC,nfile) lonlastv, cenlonv use rqstfld_mod, only: use nemsio_module, only: nemsio_gfile, nemsio_getfilehead, nemsio_close, nemsio_getheadvar + use upp_math, only: h2u ! ! INCLUDE/SET PARAMETERS. implicit none diff --git a/sorc/ncep_post.fd/INITPOST_NEMS_MPIIO.f b/sorc/ncep_post.fd/INITPOST_NEMS_MPIIO.f index 44d403bf3..5224be6e5 100644 --- a/sorc/ncep_post.fd/INITPOST_NEMS_MPIIO.f +++ b/sorc/ncep_post.fd/INITPOST_NEMS_MPIIO.f @@ -65,6 +65,7 @@ SUBROUTINE INITPOST_NEMS_MPIIO() use rqstfld_mod, only: ! use nemsio_module, only: nemsio_gfile, nemsio_getfilehead, nemsio_close, nemsio_getheadvar use nemsio_module_mpi + use upp_math, only: h2u ! ! INCLUDE/SET PARAMETERS. implicit none diff --git a/sorc/ncep_post.fd/INITPOST_NETCDF.f b/sorc/ncep_post.fd/INITPOST_NETCDF.f index f360e9172..d8c2d1a71 100644 --- a/sorc/ncep_post.fd/INITPOST_NETCDF.f +++ b/sorc/ncep_post.fd/INITPOST_NETCDF.f @@ -82,6 +82,7 @@ SUBROUTINE INITPOST_NETCDF(ncid3d) dxval, dyval, truelat2, truelat1, psmapf, cenlat,lonstartv, lonlastv, cenlonv, & latstartv, latlastv, cenlatv,latstart_r,latlast_r,lonstart_r,lonlast_r, STANDLON use rqstfld_mod, only: igds, avbl, iq, is + use upp_physics, only: fpvsnew !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none ! @@ -105,7 +106,6 @@ SUBROUTINE INITPOST_NETCDF(ncid3d) ! real,parameter:: con_eps =con_rd/con_rv ! real,parameter:: con_epsm1 =con_rd/con_rv-1 ! -! real,external::FPVSNEW ! This version of INITPOST shows how to initialize, open, read from, and ! close a NetCDF dataset. In order to change it to read an internal (binary) ! dataset, do a global replacement of _ncd_ with _int_. @@ -151,7 +151,6 @@ SUBROUTINE INITPOST_NETCDF(ncid3d) integer ncid3d,ncid2d,varid,nhcas real TSTART,TLMH,TSPH,ES,FACT,soilayert,soilayerb,zhour,dum, & tvll,pmll,tv, tx1, tx2 - real, external :: fpvsnew character*20,allocatable :: recname(:) integer, allocatable :: reclev(:), kmsk(:,:) diff --git a/sorc/ncep_post.fd/LFMFLD.f b/sorc/ncep_post.fd/LFMFLD.f index d0446c2d6..a5d83919b 100644 --- a/sorc/ncep_post.fd/LFMFLD.f +++ b/sorc/ncep_post.fd/LFMFLD.f @@ -40,6 +40,7 @@ !! 00-01-04 JIM TUCCILLO - MPI VERSION !! 02-04-24 MIKE BALDWIN - WRF VERSION !! 19-10-30 Bo CUI - REMOVE "GOTO" STATEMENT +!! 20-11-10 JESSE MENG - USE UPP_PHYSICS MODULE !! !! !! USAGE: CALL LFMFLD(RH3310,RH6610,RH3366,PW3310) @@ -74,10 +75,10 @@ SUBROUTINE LFMFLD(RH3310,RH6610,RH3366,PW3310) use params_mod, only: d00, d50, pq0, a2, a3, a4, h1, d01, gi use ctlblk_mod, only: jsta, jend, modelname, spval, im use physcons_post, only: con_rd, con_rv, con_eps, con_epsm1 + use upp_physics, only: FPVSNEW implicit none - real,external::FPVSNEW !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! real,PARAMETER :: RHOWAT=1.E3 diff --git a/sorc/ncep_post.fd/LFMFLD_GFS.f b/sorc/ncep_post.fd/LFMFLD_GFS.f index 37df02a78..e89436e39 100644 --- a/sorc/ncep_post.fd/LFMFLD_GFS.f +++ b/sorc/ncep_post.fd/LFMFLD_GFS.f @@ -43,6 +43,7 @@ !! HAVE DIFFERENT THICKNESS AS MESO AND USE DP !! RATHER THAN DZ !! 19-10-30 Bo CUI - REMOVE "GOTO" STATEMENT +!! 20-11-10 JESSE MENG - USE UPP_PHYSICS MODULE !! !! !! USAGE: CALL LFMFLD(RH3310,RH6610,RH3366,PW3310) @@ -76,6 +77,7 @@ SUBROUTINE LFMFLD_GFS(RH4410,RH7294,RH4472,RH3310) use masks, only: lmh use params_mod, only: d00 use ctlblk_mod, only: jsta, jend, spval, im + use upp_physics, only: FPVSNEW ! implicit none ! @@ -96,7 +98,6 @@ SUBROUTINE LFMFLD_GFS(RH4410,RH7294,RH4472,RH3310) integer I,J,L,LLMH real P4410, P7294,P4472,P3310,Q4410,Q7294,Q4472,Q3310,QS4410, & QS7294,QS4472,QS3310,PS,P33,DP1,DP2,DP3,DP4 - real,external :: fpvsnew !*********************************************************************** ! START LFMFLD HERE diff --git a/sorc/ncep_post.fd/MDL2P.f b/sorc/ncep_post.fd/MDL2P.f index f8a38443f..d08e91d4d 100644 --- a/sorc/ncep_post.fd/MDL2P.f +++ b/sorc/ncep_post.fd/MDL2P.f @@ -26,6 +26,8 @@ !! 14-02-26 S Moorthi - threading datapd assignment !! 19-10-30 B CUI - REMOVE "GOTO" STATEMENT !! 20-03-25 J MENG - remove grib1 +!! 20-05-20 J MENG - CALRH unification with NAM scheme +!! 20-11-10 J MENG - USE UPP_PHYSICS MODULE !! !! USAGE: CALL MDL2P !! INPUT ARGUMENT LIST: @@ -82,6 +84,7 @@ SUBROUTINE MDL2P(iostatusD3D) imp_physics use rqstfld_mod, only: IGET, LVLS, ID, IAVBLFLD, LVLSXML use gridspec_mod, only: GRIDTYPE, MAPTYPE, DXVAL + use upp_physics, only: FPVSNEW, CALRH !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! implicit none @@ -133,7 +136,6 @@ SUBROUTINE MDL2P(iostatusD3D) integer I,J,L,LP,LL,LLMH,JJB,JJE,II,JJ,LI,IFINCR,ITD3D,ista,imois,luhi,la real fact,ALPSL,PSFC,QBLO,PNL1,TBLO,TVRL,TVRBLO,FAC,PSLPIJ, & ALPTH,AHF,PDV,QL,TVU,TVD,GAMMAS,QSAT,RHL,ZL,TL,PL,ES,part,dum1 - real,external :: fpvsnew logical log1 real dxm, tem, zero ! @@ -1259,13 +1261,7 @@ SUBROUTINE MDL2P(iostatusD3D) ENDDO ENDDO ! - IF(MODELNAME == 'GFS' .or. MODELNAME == 'FV3R')THEN - CALL CALRH_GFS(EGRID2(1,jsta),TSL(1,jsta),QSL(1,jsta),EGRID1(1,jsta)) - ELSEIF (MODELNAME == 'RAPR')THEN - CALL CALRH_GSD(EGRID2(1,jsta),TSL(1,jsta),QSL(1,jsta),EGRID1(1,jsta)) - ELSE - CALL CALRH(EGRID2(1,jsta),TSL(1,jsta),QSL(1,jsta),EGRID1(1,jsta)) - END IF + CALL CALRH(EGRID2(1,jsta),TSL(1,jsta),QSL(1,jsta),EGRID1(1,jsta)) !$omp parallel do private(i,j) DO J=JSTA,JEND diff --git a/sorc/ncep_post.fd/MDL2STD_P.f b/sorc/ncep_post.fd/MDL2STD_P.f index 40fe422b0..c1860106e 100644 --- a/sorc/ncep_post.fd/MDL2STD_P.f +++ b/sorc/ncep_post.fd/MDL2STD_P.f @@ -9,6 +9,8 @@ !! !! PROGRAM HISTORY LOG: !! 19-09-24 Y Mao - REWRITTEN FROM MISCLN.f +!! 20-05-20 J MENG - CALRH unification with NAM scheme +!! 20-11-10 J MENG - USE UPP_PHYSICS MODULE !! !! USAGE: CALL MDL2STD_P !! INPUT ARGUMENT LIST: @@ -44,6 +46,7 @@ SUBROUTINE MDL2STD_P() jsta_2l, jend_2u, MODELNAME use rqstfld_mod, only: iget, lvls, iavblfld, lvlsxml use grib2_module, only: pset + use upp_physics, only: CALRH !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! @@ -481,16 +484,7 @@ SUBROUTINE MDL2STD_P() EGRID4(1:IM,JSTA:JEND)=QFD(1:IM,JSTA:JEND,IFD,2) ! Q EGRID1 = SPVAL - IF(MODELNAME == 'GFS' .or. MODELNAME == 'FV3R')THEN - CALL CALRH_GFS(EGRID2(1,jsta),EGRID3(1,jsta),& - EGRID4(1,jsta), EGRID1(1,jsta)) - ELSEIF (MODELNAME == 'RAPR')THEN - CALL CALRH_GSD(EGRID2(1,jsta),EGRID3(1,jsta),& - EGRID4(1,jsta), EGRID1(1,jsta)) - ELSE - CALL CALRH(EGRID2(1,jsta),EGRID3(1,jsta),& - EGRID4(1,jsta), EGRID1(1,jsta)) - END IF + CALL CALRH(EGRID2(1,jsta),EGRID3(1,jsta),EGRID4(1,jsta),EGRID1(1,jsta)) !$omp parallel do private(i,j) DO J=JSTA,JEND diff --git a/sorc/ncep_post.fd/MDL2THANDPV.f b/sorc/ncep_post.fd/MDL2THANDPV.f index 7878ba96b..d561b9557 100644 --- a/sorc/ncep_post.fd/MDL2THANDPV.f +++ b/sorc/ncep_post.fd/MDL2THANDPV.f @@ -13,6 +13,8 @@ !! 14-03-06 S. Moorthi - updated for threading and some optimization !! 16-12-19 G.P. Lou - Added A-grid regional models !! 20-03-25 J MENG - remove grib1 +!! 20-03-25 J MENG - remove grib1 +!! 20-11-10 J MENG - USE UPP_MATH and UPP_PHYSICS MODULES !! !! !! USAGE: CALL MDL2THANDPV @@ -49,6 +51,8 @@ SUBROUTINE MDL2THANDPV(kth,kpv,th,pv) im, jm, jsta, jend, jsta_m, jend_m, modelname, global,gdsdegr,me use RQSTFLD_mod, only: iget, lvls, id, iavblfld, lvlsxml use gridspec_mod, only: gridtype,dyval + use upp_physics, only: FPVSNEW + use upp_math, only: DVDXDUDY, DDVDX, DDUDY, UUAVG, h2u !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none ! @@ -77,7 +81,7 @@ SUBROUTINE MDL2THANDPV(kth,kpv,th,pv) integer, dimension(im) :: iw, ie integer I,J,L,K,lp,imb2,ip1,im1,ii,jj,jmt2,ihw,ihe real DVDX,DUDY,UAVG,TPHI, es, qstl, eradi, tem - real,external :: fpvsnew + real, allocatable :: DVDXL(:,:,:), DUDYL(:,:,:), UAVGL(:,:,:) ! ! !****************************************************************************** @@ -481,8 +485,19 @@ SUBROUTINE MDL2THANDPV(kth,kpv,th,pv) ENDIF !regional models and A-grid end here !----------------------------------------------------------------- ELSE IF (GRIDTYPE == 'B')THEN + allocate(DVDXL(1:im,jsta_m:jend_m,lm)) + allocate(DUDYL(1:im,jsta_m:jend_m,lm)) + allocate(UAVGL(1:im,jsta_m:jend_m,lm)) DO L=1,LM CALL EXCH(VH(1:IM,JSTA_2L:JEND_2U,L)) + CALL DVDXDUDY(UH(:,:,L),VH(:,:,L)) + DO J=JSTA_m,Jend_m + DO I=2,im-1 + DVDXL(I,J,L) = DDVDX(I,J) + DUDYL(I,J,L) = DDUDY(I,J) + UAVGL(I,J,L) = UUAVG(I,J) + END DO + END DO END DO DO J=JSTA_m,Jend_m JMT2=JM/2+1 @@ -499,12 +514,9 @@ SUBROUTINE MDL2THANDPV(kth,kpv,th,pv) DUM1D3(L) = (T(ip1,J,L) - T(im1,J,L)) * wrk2(i,j) !dt/dx DUM1D2(L) = (PMID(I,J+1,L)-PMID(I,J-1,L)) * wrk3(i,j) !dp/dy DUM1D4(L) = (T(I,J+1,L)-T(I,J-1,L)) * wrk3(i,j) !dt/dy - DVDX = (0.5*(VH(I,J,L)+VH(I,J-1,L))-0.5*(VH(IM1,J,L) & - + VH(IM1,J-1,L)))*wrk2(i,j)*2.0 - DUDY = (0.5*(UH(I,J,L)+UH(I-1,J,L))-0.5*(UH(I,J-1,L) & - + UH(I-1,J-1,L)))*wrk3(i,j)*2.0 - UAVG = 0.25*(UH(IM1,J-1,L)+UH(IM1,J,L) & - & + UH(I,J-1,L)+UH(I,J,L)) + DVDX = DVDXL(I,J,L) + DUDY = DUDYL(I,J,L) + UAVG = UAVGL(I,J,L) ! is there a (f+tan(phi)/erad)*u term? DUM1D6(L) = DVDX - DUDY + F(I,J) + UAVG*TAN(TPHI)/ERAD !vort @@ -560,6 +572,7 @@ SUBROUTINE MDL2THANDPV(kth,kpv,th,pv) END IF END DO END DO + deallocate(DVDXL,DUDYL,UAVGL) ELSE IF (GRIDTYPE == 'E')THEN DO J=JSTA_m,Jend_m JMT2 = JM/2+1 diff --git a/sorc/ncep_post.fd/MDLFLD.f b/sorc/ncep_post.fd/MDLFLD.f index d04c9cc57..dfe0430c0 100644 --- a/sorc/ncep_post.fd/MDLFLD.f +++ b/sorc/ncep_post.fd/MDLFLD.f @@ -39,6 +39,9 @@ !! 15-11-03 S Moorthi - fix a bug in "RELATIVE HUMIDITY ON MDLSURFACES" sectio logic !! 19-10-30 Bo CUI - REMOVE "GOTO" STATEMENT !! 20-03-24 J MENG - remove grib1 +!! 20-05-20 J MENG - CALRH unification with NAM scheme +!! 20-11-10 J MENG - USE UPP_MATH MODULE +!! 20-11-10 J MENG - USE UPP_PHYSICS MODULE !! !! USAGE: CALL MDLFLD !! INPUT ARGUMENT LIST: @@ -95,6 +98,8 @@ SUBROUTINE MDLFLD me, dt, avrain, theat, ifhr, ifmin, avcnvc, lp1, im, jm use rqstfld_mod, only: iget, id, lvls, iavblfld, lvlsxml use gridspec_mod, only: gridtype,maptype,dxval + use upp_physics, only: CALRH, CALCAPE + use upp_math, only: H2U, H2V, U2H, V2H ! !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - @@ -1518,13 +1523,9 @@ SUBROUTINE MDLFLD Q1D(I,J) = Q(I,J,LL) ENDDO ENDDO - IF(MODELNAME == 'GFS')THEN - CALL CALRH_GFS(P1D(1,jsta),T1D(1,jsta),Q1D(1,jsta),EGRID4(1,jsta)) - ELSE IF (MODELNAME == 'RAPR')THEN - CALL CALRH_GSD(P1D(1,jsta),T1D(1,jsta),Q1D(1,jsta),EGRID4(1,jsta)) - ELSE - CALL CALRH(P1D(1,jsta),T1D(1,jsta),Q1D(1,jsta),EGRID4(1,jsta)) - END IF + + CALL CALRH(P1D(1,jsta),T1D(1,jsta),Q1D(1,jsta),EGRID4(1,jsta)) + !$omp parallel do private(i,j) DO J=JSTA,JEND DO I=1,IM diff --git a/sorc/ncep_post.fd/MISCLN.f b/sorc/ncep_post.fd/MISCLN.f index d57af3dde..6854dc9a5 100644 --- a/sorc/ncep_post.fd/MISCLN.f +++ b/sorc/ncep_post.fd/MISCLN.f @@ -42,6 +42,7 @@ !! levels input from control file !! 19-09-03 J Meng - ADD CAPE related variables for HRRR !! 20-03-24 J Meng - remove grib1 +!! 20-11-10 J Meng - USE UPP_PHYSICS MODULE !! !! USAGE: CALL MISCLN !! INPUT ARGUMENT LIST: @@ -77,8 +78,8 @@ !! MACHINE : CRAY C-90 !! SUBROUTINE MISCLN -! +! ! use vrbls3d, only: pmid, uh, vh, t, zmid, pint, alpint, q, omga use vrbls3d, only: catedr,mwt,gtg @@ -91,6 +92,7 @@ SUBROUTINE MISCLN jsta_2l, jend_2u, MODELNAME use rqstfld_mod, only: iget, lvls, id, iavblfld, lvlsxml use grib2_module, only: pset + use upp_physics, only: FPVSNEW, CALRH_PW, CALCAPE, CALCAPE2 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none ! @@ -149,7 +151,6 @@ SUBROUTINE MISCLN REAL, allocatable :: HTFDCTL(:) integer, allocatable :: ITYPEFDLVLCTL(:) - real,external :: fpvsnew ! !**************************************************************************** ! START MISCLN HERE. diff --git a/sorc/ncep_post.fd/OTLFT.f b/sorc/ncep_post.fd/OTLFT.f index a7416fad1..6c9a2b3ee 100644 --- a/sorc/ncep_post.fd/OTLFT.f +++ b/sorc/ncep_post.fd/OTLFT.f @@ -18,6 +18,7 @@ !! 00-01-04 JIM TUCCILLO - MPI VERSION !! 02-06-17 MIKE BALDWIN - WRF VERSION !! 11-04-12 GEOFF MANIKIN - USE VIRTUAL TEMPERATURE +!! 20-11-10 JESSE MENG - USE UPP_PHYSICS MODULE !! !! USAGE: CALL OTLFT(PBND,TBND,QBND,SLINDX) !! INPUT ARGUMENT LIST: @@ -53,12 +54,12 @@ SUBROUTINE OTLFT(PBND,TBND,QBND,SLINDX) PL, RDP, THE0, STHE, RDTHE, TTBL use ctlblk_mod, only: JSTA, JEND, IM use params_mod, only: D00, H10E5, CAPA, ELOCP, EPS, ONEPS + use upp_physics, only: FPVSNEW !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none ! ! SET LOCAL PARAMETERS. real,PARAMETER :: D8202=.820231E0 , H5E4=5.E4 , P500=50000. - real,external::FPVSNEW ! ! DECLARE VARIABLES. diff --git a/sorc/ncep_post.fd/OTLIFT.f b/sorc/ncep_post.fd/OTLIFT.f index 21acf7cf6..9e911d846 100644 --- a/sorc/ncep_post.fd/OTLIFT.f +++ b/sorc/ncep_post.fd/OTLIFT.f @@ -18,6 +18,7 @@ !! 01-10-25 H CHUANG - MODIFIED TO PROCESS HYBRID MODEL OUTPUT !! 02-06-11 MIKE BALDWIN - WRF VERSION !! 11-04-12 GEOFF MANIKIN - USE VIRTUAL TEMPERATURE +!! 20-11-10 JESSE MENG - USE UPP_PHYSICS MODULE !! !! USAGE: CALL OTLIFT(SLINDX) !! INPUT ARGUMENT LIST: @@ -48,6 +49,7 @@ SUBROUTINE OTLIFT(SLINDX) RDP, THE0, STHE, RDTHE, TTBL use ctlblk_mod, only: JSTA, JEND, IM use params_mod, only: D00,H10E5, CAPA, ELOCP, EPS, ONEPS + use upp_physics, only: FPVSNEW ! ! @@ -55,7 +57,6 @@ SUBROUTINE OTLIFT(SLINDX) ! ! SET LOCAL PARAMETERS. real,PARAMETER :: D8202=.820231E0 , H5E4=5.E4 , P500=50000. - real,external::FPVSNEW ! ! DECLARE VARIABLES. diff --git a/sorc/ncep_post.fd/SURFCE.f b/sorc/ncep_post.fd/SURFCE.f index e3d40f9d9..85945fd94 100644 --- a/sorc/ncep_post.fd/SURFCE.f +++ b/sorc/ncep_post.fd/SURFCE.f @@ -34,6 +34,8 @@ !! - 14-02-26 S Moorthi - threading datapd assignment !! - 14-11-26 S Moorthi - cleanup and some bug fix (may be?) !! - 20-03-25 J MENG - remove grib1 +!! - 20-05-20 J MENG - CALRH unification with NAM scheme +!! - 20-11-10 J MENG - USE UPP_PHYSICS MODULE !! !! USAGE: CALL SURFCE !! INPUT ARGUMENT LIST: @@ -98,6 +100,7 @@ SUBROUTINE SURFCE lp1, imp_physics, me, asrfc, tsrfc, pt, pdtop, & mpi_comm_comp, im, jm, prec_acc_dt1 use rqstfld_mod, only: iget, lvls, id, iavblfld, lvlsxml + use upp_physics, only: fpvsnew, CALRH !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none ! @@ -146,7 +149,6 @@ SUBROUTINE SURFCE RC,SFCTMP,SNCOVR,FACTRS,SOLAR, s,tk,tl,w,t2c,dlt,APE, & qv,e,dwpt,dum1,dum2,dum3,dum1s,dum3s,dum21,dum216,es - real,external :: fpvsnew !**************************************************************************** ! @@ -1443,14 +1445,9 @@ SUBROUTINE SURFCE Q1D(I,J) = QSHLTR(I,J) ENDDO ENDDO -! CALL CALRH(PSHLTR,TSHLTR,QSHLTR,EGRID1) - IF(MODELNAME == 'GFS')THEN - CALL CALRH_GFS(P1D,T1D,Q1D,EGRID1(1,jsta)) - ELSEIF(MODELNAME == 'RAPR')THEN - CALL CALRH_GSD(P1D,T1D,Q1D,EGRID1(1,jsta)) - ELSE - CALL CALRH(P1D,T1D,Q1D,EGRID1(1,jsta)) - END IF + + CALL CALRH(P1D,T1D,Q1D,EGRID1(1,jsta)) + if (allocated(q1d)) deallocate(q1d) !$omp parallel do private(i,j) DO J=JSTA,JEND diff --git a/sorc/ncep_post.fd/GRIDAVG.f b/sorc/ncep_post.fd/UPP_MATH.f similarity index 55% rename from sorc/ncep_post.fd/GRIDAVG.f rename to sorc/ncep_post.fd/UPP_MATH.f index 40b880beb..2b1ad4a75 100644 --- a/sorc/ncep_post.fd/GRIDAVG.f +++ b/sorc/ncep_post.fd/UPP_MATH.f @@ -1,3 +1,116 @@ +!> @file +! +!> SUBPROGRAM: UPP_MATH +!! @author JMENG @date 2020-05-20 +!! +!! A collection of UPP subroutines for numerical math functions calculation. +!! +!! DVDXDUDY +!! computes dudy, dvdx, uwnd +!! +!! H2U, H2V, U2H, V2H +!! interpolates variables between U, V, H, points +!! adopted from UPP subroutine GRIDAVG.f +!! +!! PROGRAM HISTORY LOG: +!! MAY 20 2020 Jesse Meng Initial code +!!------------------------------------------------------------------------ +!! + module upp_math + + use masks, only: dx, dy + use ctlblk_mod, only: im, jsta_2l, jend_2u, jsta_m, jend_m, spval + use gridspec_mod, only: gridtype +! + implicit none + + private + + public :: DDVDX, DDUDY, UUAVG + public :: dvdxdudy + public :: H2U, H2V, U2H, V2H + + REAL, ALLOCATABLE :: DDVDX(:,:) + REAL, ALLOCATABLE :: DDUDY(:,:) + REAL, ALLOCATABLE :: UUAVG(:,:) +! +!------------------------------------------------------------------------------------- +! + contains +! +!------------------------------------------------------------------------------------- + subroutine dvdxdudy(uwnd,vwnd) +! + implicit none + + REAL, dimension(im,jsta_2l:jend_2u), intent(in) :: UWND, VWND +! +!-- local variables +!-- + integer i, j + real r2dx, r2dy + INTEGER, allocatable :: IHE(:),IHW(:) +! + IF(GRIDTYPE == 'A')THEN +!$omp parallel do private(i,j,r2dx,r2dy) + DO J=JSTA_M,JEND_M + DO I=2,IM-1 + IF(VWND(I+1,J) @file +! +!> SUBPROGRAM: UPP_PHYSICS +!! @author JMENG @date 2020-05-20 +!! +!! A collection of UPP subroutines for physics variables calculation. +!! +!! CALCAPE +!! Compute CAPE/CINS and other storm related variables. +!! +!! CALCAPE2 +!! Compute additional storm related variables. +!! +!! CALRH +!! CALRH_NAM +!! CALRH_GFS +!! CALRH_GSD +!! Compute RH using various algorithms. +!! The NAM v4.1.18 ALGORITHM (CALRH_NAM) is selected as default for +!! NMMB and FV3GFS, FV3GEFS, and FV3R for the UPP 2020 unification. +!! +!! CALRH_PW +!! Algorithm use at GSD for RUC and Rapid Refresh +!! +!! FPVSNEW +!! Compute saturation vapor pressure. +!! +!! TVIRTUAL +!! Compute virtual temperature. +!! +!! PROGRAM HISTORY LOG: +!! MAY, 2020 Jesse Meng Initial code +!!------------------------------------------------------------------------------------- +!! + module upp_physics + + implicit none + + private + + public :: CALCAPE, CALCAPE2 + public :: CALRH + public :: CALRH_GFS, CALRH_GSD, CALRH_NAM + public :: CALRH_PW + public :: FPVSNEW + public :: TVIRTUAL + + contains +! +!------------------------------------------------------------------------------------- +! + SUBROUTINE CALRH(P1,T1,Q1,RH) + + use ctlblk_mod, only: im, jsta, jend, MODELNAME + implicit none + + REAL,dimension(IM,jsta:jend),intent(in) :: P1,T1 + REAL,dimension(IM,jsta:jend),intent(inout) :: Q1 + REAL,dimension(IM,jsta:jend),intent(out) :: RH + + IF(MODELNAME == 'RAPR')THEN + CALL CALRH_GSD(P1,T1,Q1,RH) + ELSE + CALL CALRH_NAM(P1,T1,Q1,RH) + END IF + + END SUBROUTINE CALRH +! +!------------------------------------------------------------------------------------- +! + SUBROUTINE CALRH_NAM(P1,T1,Q1,RH) +! SUBROUTINE CALRH(P1,T1,Q1,RH) +!$$$ SUBPROGRAM DOCUMENTATION BLOCK +! . . . +! SUBPROGRAM: CALRH COMPUTES RELATIVE HUMIDITY +! PRGRMMR: TREADON ORG: W/NP2 DATE: 92-12-22 +! +! ABSTRACT: +! THIS ROUTINE COMPUTES RELATIVE HUMIDITY GIVEN PRESSURE, +! TEMPERATURE, SPECIFIC HUMIDITY. AN UPPER AND LOWER BOUND +! OF 100 AND 1 PERCENT RELATIVE HUMIDITY IS ENFORCED. WHEN +! THESE BOUNDS ARE APPLIED THE PASSED SPECIFIC HUMIDITY +! ARRAY IS ADJUSTED AS NECESSARY TO PRODUCE THE SET RELATIVE +! HUMIDITY. +! . +! +! PROGRAM HISTORY LOG: +! ??-??-?? DENNIS DEAVEN +! 92-12-22 RUSS TREADON - MODIFIED AS DESCRIBED ABOVE. +! 98-06-08 T BLACK - CONVERSION FROM 1-D TO 2-D +! 98-08-18 MIKE BALDWIN - MODIFY TO COMPUTE RH OVER ICE AS IN MODEL +! 98-12-16 GEOFF MANIKIN - UNDO RH COMPUTATION OVER ICE +! 00-01-04 JIM TUCCILLO - MPI VERSION +! 02-06-11 MIKE BALDWIN - WRF VERSION +! 06-03-19 Wen Meng - MODIFY TOP PRESSURE to 1 PA +! +! USAGE: CALL CALRH(P1,T1,Q1,RH) +! INPUT ARGUMENT LIST: +! P1 - PRESSURE (PA) +! T1 - TEMPERATURE (K) +! Q1 - SPECIFIC HUMIDITY (KG/KG) +! +! OUTPUT ARGUMENT LIST: +! RH - RELATIVE HUMIDITY (DECIMAL FORM) +! Q1 - ADJUSTED SPECIFIC HUMIDITY (KG/KG) +! +! OUTPUT FILES: +! NONE +! +! SUBPROGRAMS CALLED: +! UTILITIES: +! LIBRARY: +! NONE +! +! ATTRIBUTES: +! LANGUAGE: FORTRAN +! MACHINE : CRAY C-90 +!$$$ +! + use params_mod, only: PQ0, a2, a3, a4, rhmin + use ctlblk_mod, only: jsta, jend, spval, im +!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + implicit none +! +! SET PARAMETER. +! +! DECLARE VARIABLES. +! + REAL,dimension(IM,jsta:jend),intent(in) :: P1,T1 + REAL,dimension(IM,jsta:jend),intent(inout) :: Q1 + REAL,dimension(IM,jsta:jend),intent(out) :: RH + REAL QC + integer I,J +!*************************************************************** +! +! START CALRH. +! + DO J=JSTA,JEND + DO I=1,IM + 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)) +! + RH(I,J) = Q1(I,J)/QC +! +! BOUNDS CHECK +! + IF (RH(I,J) > 1.0) THEN + RH(I,J) = 1.0 + Q1(I,J) = RH(I,J)*QC + ENDIF + IF (RH(I,J) < RHmin) THEN !use smaller RH limit for stratosphere + RH(I,J) = RHmin + Q1(I,J) = RH(I,J)*QC + ENDIF +! + ENDIF + ELSE + RH(I,J) = SPVAL + ENDIF + ENDDO + ENDDO +! +! +! END SUBROUTINE CALRH + END SUBROUTINE CALRH_NAM +! +!------------------------------------------------------------------------------------- +! + SUBROUTINE CALRH_GFS(P1,T1,Q1,RH) +!$$$ SUBPROGRAM DOCUMENTATION BLOCK +! . . . +! SUBPROGRAM: CALRH COMPUTES RELATIVE HUMIDITY +! PRGRMMR: TREADON ORG: W/NP2 DATE: 92-12-22 +! +! ABSTRACT: +! THIS ROUTINE COMPUTES RELATIVE HUMIDITY GIVEN PRESSURE, +! TEMPERATURE, SPECIFIC HUMIDITY. AN UPPER AND LOWER BOUND +! OF 100 AND 1 PERCENT RELATIVE HUMIDITY IS ENFORCED. WHEN +! THESE BOUNDS ARE APPLIED THE PASSED SPECIFIC HUMIDITY +! ARRAY IS ADJUSTED AS NECESSARY TO PRODUCE THE SET RELATIVE +! HUMIDITY. +! . +! +! PROGRAM HISTORY LOG: +! ??-??-?? DENNIS DEAVEN +! 92-12-22 RUSS TREADON - MODIFIED AS DESCRIBED ABOVE. +! 98-06-08 T BLACK - CONVERSION FROM 1-D TO 2-D +! 98-08-18 MIKE BALDWIN - MODIFY TO COMPUTE RH OVER ICE AS IN MODEL +! 98-12-16 GEOFF MANIKIN - UNDO RH COMPUTATION OVER ICE +! 00-01-04 JIM TUCCILLO - MPI VERSION +! 02-06-11 MIKE BALDWIN - WRF VERSION +! 13-08-13 S. Moorthi - Threading +! 06-03-19 Wen Meng - MODIFY TOP PRESSURE to 1 PA +! +! USAGE: CALL CALRH(P1,T1,Q1,RH) +! INPUT ARGUMENT LIST: +! P1 - PRESSURE (PA) +! T1 - TEMPERATURE (K) +! Q1 - SPECIFIC HUMIDITY (KG/KG) +! +! OUTPUT ARGUMENT LIST: +! RH - RELATIVE HUMIDITY (DECIMAL FORM) +! Q1 - ADJUSTED SPECIFIC HUMIDITY (KG/KG) +! +! OUTPUT FILES: +! NONE +! +! SUBPROGRAMS CALLED: +! UTILITIES: +! LIBRARY: +! NONE +! +! ATTRIBUTES: +! LANGUAGE: FORTRAN +! MACHINE : CRAY C-90 +!$$$ +! + use params_mod, only: rhmin + use ctlblk_mod, only: jsta, jend, spval, im +!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + implicit none +! + real,parameter:: con_rd =2.8705e+2 ! gas constant air (J/kg/K) + real,parameter:: con_rv =4.6150e+2 ! gas constant H2O + real,parameter:: con_eps =con_rd/con_rv + real,parameter:: con_epsm1 =con_rd/con_rv-1 +! real,external::FPVSNEW + +! INTERFACE +! ELEMENTAL FUNCTION FPVSNEW (t) +! REAL FPVSNEW +! REAL, INTENT(IN) :: t +! END FUNCTION FPVSNEW +! END INTERFACE +! + REAL,dimension(IM,jsta:jend),intent(in) :: P1,T1 + REAL,dimension(IM,jsta:jend),intent(inout):: Q1,RH + REAL ES,QC + integer :: I,J +!*************************************************************** +! +! START CALRH. +! +!$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 (ABS(P1(I,J)) > 1.0) THEN +! IF (P1(I,J) > 1.0) THEN + IF (P1(I,J) >= 1.0) THEN + ES = MIN(FPVSNEW(T1(I,J)),P1(I,J)) + QC = CON_EPS*ES/(P1(I,J)+CON_EPSM1*ES) + +! QC=PQ0/P1(I,J)*EXP(A2*(T1(I,J)-A3)/(T1(I,J)-A4)) + + RH(I,J) = min(1.0,max(Q1(I,J)/QC,rhmin)) + q1(i,j) = rh(i,j)*qc + +! BOUNDS CHECK +! +! IF (RH(I,J) > 1.0) THEN +! RH(I,J) = 1.0 +! Q1(I,J) = RH(I,J)*QC +! ELSEIF (RH(I,J) < RHmin) THEN !use smaller RH limit for stratosphere +! RH(I,J) = RHmin +! Q1(I,J) = RH(I,J)*QC +! ENDIF + + ENDIF + ELSE + RH(I,J) = SPVAL + ENDIF + ENDDO + ENDDO + + END SUBROUTINE CALRH_GFS +! +!------------------------------------------------------------------------------------- +! + SUBROUTINE CALRH_GSD(P1,T1,Q1,RHB) +! +! Algorithm use at GSD for RUC and Rapid Refresh +!------------------------------------------------------------------ +! + + use ctlblk_mod, only: jsta, jend, im + + implicit none + + integer :: j, i + real :: tx, pol, esx, es, e + real, dimension(im,jsta:jend) :: P1, T1, Q1, RHB + + + DO J=JSTA,JEND + DO I=1,IM + +! - compute relative humidity + Tx=T1(I,J)-273.15 + POL = 0.99999683 + TX*(-0.90826951E-02 + & + TX*(0.78736169E-04 + TX*(-0.61117958E-06 + & + TX*(0.43884187E-08 + TX*(-0.29883885E-10 + & + TX*(0.21874425E-12 + TX*(-0.17892321E-14 + & + TX*(0.11112018E-16 + TX*(-0.30994571E-19))))))))) + esx = 6.1078/POL**8 + + ES = esx + E = P1(I,J)/100.*Q1(I,J)/(0.62197+Q1(I,J)*0.37803) + RHB(I,J) = MIN(1.,E/ES) + + ENDDO + ENDDO + + END SUBROUTINE CALRH_GSD +! +!------------------------------------------------------------------------------------- +! + SUBROUTINE CALRH_PW(RHPW) +! +! Algorithm use at GSD for RUC and Rapid Refresh +!------------------------------------------------------------------ +! + + use vrbls3d, only: q, pmid, t + use params_mod, only: g + use ctlblk_mod, only: lm, jsta, jend, lm, im +!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + implicit none + + real,PARAMETER :: svp1=6.1153,svp2=17.67,svp3=29.65 + + REAL, dimension(im,jsta:jend):: PW, PW_SAT, RHPW + REAL deltp,sh,qv,temp,es,qs,qv_sat + integer i,j,l,k,ka,kb + + pw = 0. + pw_sat = 0. + rhpw = 0. + + DO L=1,LM + k=lm-l+1 + DO J=JSTA,JEND + DO I=1,IM +! -- use specific humidity for PW calculation + sh = q(i,j,k) + qv = sh/(1.-sh) + KA = MAX(1,K-1) + KB = MIN(LM,K+1) + +! assumes that P is in mb at this point - be careful! + DELTP = 0.5*(PMID(I,J,KB)-PMID(I,J,KA)) + PW(I,J) = PW(I,J) + sh *DELTP/G + +!Csgb -- Add more for RH w.r.t. PW-sat + + temp = T(I,J,K) +! --- use saturation mixing ratio w.r.t. water here +! for this check. + es = svp1*exp(SVP2*(Temp-273.15)/(Temp-SVP3)) +! -- get saturation specific humidity (w.r.t. total air) + qs = 0.62198*es/(pmid(i,j,k)*1.e-2-0.37802*es) +! -- get saturation mixing ratio (w.r.t. dry air) + qv_sat = qs/(1.-qs) + + pw_sat(i,j) = pw_sat(i,j) + max(sh,Qs)*DELTP/G + + if (i==120 .and. j==120 ) & + write (6,*)'pw-sat', temp, sh, qs, pmid(i,j,kb) & + ,pmid(i,j,ka),pw(i,j),pw_sat(i,j) + +!sgb - This IS RH w.r.t. PW-sat. + RHPW (i,j) = min(1.,PW(i,j) / pw_sat(i,j)) * 100. + + ENDDO + ENDDO + ENDDO + + END SUBROUTINE CALRH_PW +! +!------------------------------------------------------------------------------------- +! + elemental function fpvsnew(t) +!$$$ Subprogram Documentation Block +! +! Subprogram: fpvsnew Compute saturation vapor pressure +! Author: N Phillips w/NMC2X2 Date: 30 dec 82 +! +! Abstract: Compute saturation vapor pressure from the temperature. +! A linear interpolation is done between values in a lookup table +! computed in gpvs. See documentation for fpvsx for details. +! Input values outside table range are reset to table extrema. +! The interpolation accuracy is almost 6 decimal places. +! On the Cray, fpvs is about 4 times faster than exact calculation. +! This function should be expanded inline in the calling routine. +! +! Program History Log: +! 91-05-07 Iredell made into inlinable function +! 94-12-30 Iredell expand table +! 1999-03-01 Iredell f90 module +! 2001-02-26 Iredell ice phase +! +! Usage: pvs=fpvsnew(t) +! +! Input argument list: +! t Real(krealfp) temperature in Kelvin +! +! Output argument list: +! fpvsnew Real(krealfp) saturation vapor pressure in Pascals +! +! Attributes: +! Language: Fortran 90. +! +!$$$ + implicit none + integer,parameter:: nxpvs=7501 + real,parameter:: con_ttp =2.7316e+2 ! temp at H2O 3pt + real,parameter:: con_psat =6.1078e+2 ! pres at H2O 3pt + real,parameter:: con_cvap =1.8460e+3 ! spec heat H2O gas (J/kg/K) + real,parameter:: con_cliq =4.1855e+3 ! spec heat H2O liq + real,parameter:: con_hvap =2.5000e+6 ! lat heat H2O cond + real,parameter:: con_rv =4.6150e+2 ! gas constant H2O + real,parameter:: con_csol =2.1060e+3 ! spec heat H2O ice + real,parameter:: con_hfus =3.3358e+5 ! lat heat H2O fusion + real,parameter:: tliq=con_ttp + real,parameter:: tice=con_ttp-20.0 + real,parameter:: dldtl=con_cvap-con_cliq + real,parameter:: heatl=con_hvap + real,parameter:: xponal=-dldtl/con_rv + real,parameter:: xponbl=-dldtl/con_rv+heatl/(con_rv*con_ttp) + real,parameter:: dldti=con_cvap-con_csol + real,parameter:: heati=con_hvap+con_hfus + real,parameter:: xponai=-dldti/con_rv + real,parameter:: xponbi=-dldti/con_rv+heati/(con_rv*con_ttp) + real tr,w,pvl,pvi + real fpvsnew + real,intent(in):: t + integer jx + real xj,x,tbpvs(nxpvs),xp1 + real xmin,xmax,xinc,c2xpvs,c1xpvs +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + xmin=180.0 + xmax=330.0 + xinc=(xmax-xmin)/(nxpvs-1) +! c1xpvs=1.-xmin/xinc + c2xpvs=1./xinc + c1xpvs=1.-xmin*c2xpvs +! xj=min(max(c1xpvs+c2xpvs*t,1.0),real(nxpvs,krealfp)) + xj=min(max(c1xpvs+c2xpvs*t,1.0),float(nxpvs)) + jx=min(xj,float(nxpvs)-1.0) + x=xmin+(jx-1)*xinc + + tr=con_ttp/x + if(x>=tliq) then + tbpvs(jx)=con_psat*(tr**xponal)*exp(xponbl*(1.-tr)) + elseif(x=tliq) then + tbpvs(jx+1)=con_psat*(tr**xponal)*exp(xponbl*(1.-tr)) + elseif(xp1 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 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 +! +!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + implicit none +! +! INCLUDE/SET PARAMETERS. CONSTANTS ARE FROM BOLTON (MWR, 1980). + real,PARAMETER :: ISMTHP=2,ISMTHT=2,ISMTHQ=2 +! +! DECLARE VARIABLES. +! + integer,intent(in) :: ITYPE + real,intent(in) :: DPBND + integer, dimension(IM,Jsta:jend),intent(in) :: L1D + real, dimension(IM,Jsta:jend),intent(in) :: P1D,T1D + real, dimension(IM,jsta:jend),intent(inout) :: Q1D,CAPE,CINS,PPARC,ZEQL +! + integer, dimension(im,jsta:jend) :: IEQL, IPTB, ITHTB, PARCEL, KLRES, KHRES, LCL, IDX +! + real, dimension(im,jsta:jend) :: THESP, PSP, CAPE20, QQ, PP, THUND + REAL, ALLOCATABLE :: TPAR(:,:,:) + + 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,external :: fpvsnew + integer I,J,L,KNUML,KNUMH,LBEG,LEND,IQ, KB,ITTBK + +! integer I,J,L,KNUML,KNUMH,LBEG,LEND,IQ,IT,LMHK, KB,ITTBK +! +!************************************************************** +! START CALCAPE HERE. +! + ALLOCATE(TPAR(IM,JSTA_2L:JEND_2U,LM)) +! +! 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 +! +!$omp parallel do + DO J=JSTA,JEND + DO I=1,IM + CAPE(I,J) = D00 + CAPE20(I,J) = D00 + CINS(I,J) = D00 + LCL(I,J) = 0 + THESP(I,J) = D00 + IEQL(I,J) = LM + PARCEL(I,J) = LM + PSP(I,J) = D00 + PPARC(I,J) = D00 + THUNDER(I,J) = .TRUE. + ENDDO + ENDDO +! +!$omp parallel do + DO L=1,LM + DO J=JSTA,JEND + DO I=1,IM + TPAR(I,J,L) = D00 + ENDDO + ENDDO + ENDDO +! +! TYPE 2 CAPE/CINS: +! NOTE THAT FOR TYPE 1 CAPE/CINS ARRAYS P1D, T1D, Q1D +! ARE DUMMY ARRAYS. +! + IF (ITYPE == 2) THEN +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=1,IM + Q1D(I,J) = MIN(MAX(H1M12,Q1D(I,J)),H99999) + ENDDO + ENDDO + ENDIF +!-------FOR ITYPE=1--FIND MAXIMUM THETA E LAYER IN LOWEST DPBND ABOVE GROUND------- +!-------FOR ITYPE=2--FIND THETA E LAYER OF GIVEN T1D, Q1D, P1D--------------------- +!--------------TRIAL MAXIMUM BUOYANCY LEVEL VARIABLES------------------- + + DO KB=1,LM +!hc IF (ITYPE==2.AND.KB>1) cycle + IF (ITYPE == 1 .OR. (ITYPE == 2 .AND. KB == 1)) THEN + +!$omp parallel do private(i,j,apebtk,apespk,bqk,bqs00k,bqs10k,iq,ittbk, & +!$omp & p00k,p01k,p10k,p11k,pkl,psfck,qbtk,sqk,sqs00k, & +!$omp & sqs10k,tbtk,tpspk,tqk,tthbtk,tthesk,tthk) + DO J=JSTA,JEND + DO I=1,IM + PSFCK = PMID(I,J,NINT(LMH(I,J))) + PKL = PMID(I,J,KB) + +!hc IF (ITYPE==1.AND.(PKLPSFCK)) cycle + IF (ITYPE ==2 .OR. & + (ITYPE == 1 .AND. (PKL >= PSFCK-DPBND .AND. PKL <= PSFCK)))THEN + IF (ITYPE == 1) THEN + TBTK = T(I,J,KB) + QBTK = max(0.0, Q(I,J,KB)) + APEBTK = (H10E5/PKL)**CAPA + ELSE + PKL = P1D(I,J) + TBTK = T1D(I,J) + QBTK = max(0.0, Q1D(I,J)) + 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(I,J) = TTHK - AINT(TTHK) + ITTBK = INT(TTHK) + 1 +!--------------KEEPING INDICES WITHIN THE TABLE------------------------- + IF(ITTBK < 1) THEN + ITTBK = 1 + QQ(I,J) = D00 + ENDIF + IF(ITTBK >= JTB) THEN + ITTBK = JTB-1 + QQ(I,J) = 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(I,J) + BQS00K + SQK = (SQS10K-SQS00K)*QQ(I,J) + SQS00K + TQK = (QBTK-BQK)/SQK*RDQ + PP(I,J) = TQK-AINT(TQK) + IQ = INT(TQK)+1 +!--------------KEEPING INDICES WITHIN THE TABLE------------------------- + IF(IQ < 1) THEN + IQ = 1 + PP(I,J) = D00 + ENDIF + IF(IQ >= ITB) THEN + IQ = ITB-1 + PP(I,J) = 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(I,J) + (P01K-P00K)*QQ(I,J) & + + (P00K-P10K-P01K+P11K)*PP(I,J)*QQ(I,J) + +!!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(I,J)) THEN + PSP (I,J) = TPSPK + THESP(I,J) = TTHESK + PARCEL(I,J) = KB + ENDIF + END IF + ENDDO ! I loop + ENDDO ! J loop + END IF + ENDDO ! KB loop + +!----FIND THE PRESSURE OF THE PARCEL THAT WAS LIFTED +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=1,IM + PPARC(I,J) = PMID(I,J,PARCEL(I,J)) + ENDDO + ENDDO +! +!-----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)--------------------- + DO L=1,LM +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=1,IM + IF (PMID(I,J,L) < PSP(I,J)) LCL(I,J) = L+1 + ENDDO + ENDDO + ENDDO +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=1,IM + IF (LCL(I,J) > NINT(LMH(I,J))) LCL(I,J) = NINT(LMH(I,J)) + IF (ITYPE > 2) THEN + IF (T(I,J,LCL(I,J)) < 263.15) THEN + THUNDER(I,J) = .FALSE. + ENDIF + ENDIF + ENDDO + ENDDO +!----------------------------------------------------------------------- +!---------FIND TEMP OF PARCEL LIFTED ALONG MOIST ADIABAT (TPAR)--------- +!----------------------------------------------------------------------- + + DO L=LM,1,-1 +!--------------SCALING PRESSURE & TT TABLE INDEX------------------------ + KNUML = 0 + KNUMH = 0 + DO J=JSTA,JEND + DO I=1,IM + KLRES(I,J) = 0 + KHRES(I,J) = 0 + IF(L <= LCL(I,J)) THEN + IF(PMID(I,J,L) < PLQ)THEN + KNUML = KNUML + 1 + KLRES(I,J) = 1 + ELSE + KNUMH = KNUMH + 1 + KHRES(I,J) = 1 + ENDIF + ENDIF + ENDDO + ENDDO +!*** +!*** COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE 0) THEN + CALL TTBLEX(TPAR(1,JSTA_2L,L),TTBL,ITB,JTB,KLRES & + , PMID(1,JSTA_2L,L),PL,QQ,PP,RDP,THE0,STHE & + , RDTHE,THESP,IPTB,ITHTB) + ENDIF +!*** +!*** COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE>PLQ +!** + IF(KNUMH > 0) THEN + CALL TTBLEX(TPAR(1,JSTA_2L,L),TTBLQ,ITBQ,JTBQ,KHRES & + , PMID(1,JSTA_2L,L),PLQ,QQ,PP,RDPQ & + ,THE0Q,STHEQ,RDTHEQ,THESP,IPTB,ITHTB) + ENDIF + +!------------SEARCH FOR EQ LEVEL---------------------------------------- +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=1,IM + IF(KHRES(I,J) > 0) THEN + IF(TPAR(I,J,L) > T(I,J,L)) IEQL(I,J) = L + ENDIF + ENDDO + ENDDO +! +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=1,IM + IF(KLRES(I,J) > 0) THEN + IF(TPAR(I,J,L) > T(I,J,L) .AND. & + PMID(I,J,L)>100.) IEQL(I,J) = L + ENDIF + ENDDO + ENDDO +!----------------------------------------------------------------------- + ENDDO ! end of do l=lm,1,-1 loop +!------------COMPUTE CAPE AND CINS-------------------------------------- + LBEG = 1000 + LEND = 0 + DO J=JSTA,JEND + DO I=1,IM + LBEG = MIN(IEQL(I,J),LBEG) + LEND = MAX(LCL(I,J),LEND) + ENDDO + ENDDO +! +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=1,IM + IF(T(I,J,IEQL(I,J)) > 255.65) THEN + THUNDER(I,J) = .FALSE. + ENDIF + ENDDO + ENDDO +! + DO L=LBEG,LEND + +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=1,IM + IDX(I,J) = 0 + IF(L >= IEQL(I,J).AND.L <= LCL(I,J)) THEN + IDX(I,J) = 1 + ENDIF + ENDDO + ENDDO +! +!$omp parallel do private(i,j,gdzkl,presk,thetaa,thetap,esatp,qsatp,tvp,tv) + DO J=JSTA,JEND + DO I=1,IM + IF(IDX(I,J) > 0) THEN + PRESK = PMID(I,J,L) + GDZKL = (ZINT(I,J,L)-ZINT(I,J,L+1)) * G + ESATP = min(FPVSNEW(TPAR(I,J,L)),PRESK) + QSATP = EPS*ESATP/(PRESK-ESATP*ONEPS) +! TVP = TPAR(I,J,L)*(1+0.608*QSATP) + TVP = TVIRTUAL(TPAR(I,J,L),QSATP) + THETAP = TVP*(H10E5/PRESK)**CAPA +! TV = T(I,J,L)*(1+0.608*Q(I,J,L)) + TV = TVIRTUAL(T(I,J,L),Q(I,J,L)) + THETAA = TV*(H10E5/PRESK)**CAPA + IF(THETAP < THETAA) THEN + CINS(I,J) = CINS(I,J) + (LOG(THETAP)-LOG(THETAA))*GDZKL + ELSEIF(THETAP > THETAA) THEN + CAPE(I,J) = CAPE(I,J) + (LOG(THETAP)-LOG(THETAA))*GDZKL + IF (THUNDER(I,J) .AND. T(I,J,L) < 273.15 & + .AND. T(I,J,L) > 253.15) THEN + CAPE20(I,J) = CAPE20(I,J) + (LOG(THETAP)-LOG(THETAA))*GDZKL + ENDIF + ENDIF + ENDIF + ENDDO + ENDDO + ENDDO +! +! ENFORCE LOWER LIMIT OF 0.0 ON CAPE AND UPPER +! LIMIT OF 0.0 ON CINS. +! +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=1,IM + CAPE(I,J) = MAX(D00,CAPE(I,J)) + CINS(I,J) = MIN(CINS(I,J),D00) +! add equillibrium height + ZEQL(I,J) = ZINT(I,J,IEQL(I,J)) + IF (CAPE20(I,J) < 75.) THEN + THUNDER(I,J) = .FALSE. + ENDIF + IF (THUNDER(I,J)) THEN + THUND(I,J) = 1.0 + ELSE + THUND(I,J) = 0.0 + ENDIF + ENDDO + ENDDO +! + DEALLOCATE(TPAR) +! + END SUBROUTINE CALCAPE +! +!------------------------------------------------------------------------------------- +! + SUBROUTINE CALCAPE2(ITYPE,DPBND,P1D,T1D,Q1D,L1D, & + CAPE,CINS,LFC,ESRHL,ESRHH, & + DCAPE,DGLD,ESP) +! SUBROUTINE CALCAPE(ITYPE,DPBND,P1D,T1D,Q1D,L1D,CAPE, & +! CINS,PPARC,ZEQL,THUND) +!$$$ SUBPROGRAM DOCUMENTATION BLOCK +! . . . +! SUBPROGRAM: CALCAPE COMPUTES CAPE AND CINS +! PRGRMMR: TREADON ORG: W/NP2 DATE: 93-02-10 +! +! 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 +! 19-09-03 J MENG - MODIFIED TO ADD 0-3KM CAPE/CINS, LFC, +! EFFECTIVE HELICITY, DOWNDRAFT CAPE, +! DENDRITIC GROWTH LAYER DEPTH, ESP +! +! USAGE: CALL CALCAPE2(ITYPE,DPBND,P1D,T1D,Q1D,L1D, & +! CAPE,CINS,LFC,ESRHL,ESRHH, & +! DCAPE,DGLD,ESP) +! +! 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) +! LFC - LEVEL OF FREE CONVECTION (M) +! ESRHL - LOWER BOUND TO ACCOUNT FOR EFFECTIVE HELICITY CALCULATION +! ESRHH - UPPER BOUND TO ACCOUNT FOR EFFECTIVE HELICITY CALCULATION +! DCAPE - DOWNDRAFT CAPE (J/KG) +! DGLD - DENDRITIC GROWTH LAYER DEPTH (M) +! ESP - ENHANCED STRETCHING POTENTIAL +! +! 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: fis + use gridspec_mod, only: gridtype + use masks, only: lmh + use params_mod, only: d00, h1m12, h99999, h10e5, capa, elocp, eps, & + oneps, g, tfrz + 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, jm, me, jsta_m, jend_m +! +!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + implicit none +! +! INCLUDE/SET PARAMETERS. CONSTANTS ARE FROM BOLTON (MWR, 1980). + real,PARAMETER :: ISMTHP=2,ISMTHT=2,ISMTHQ=2 +! +! DECLARE VARIABLES. +! + integer,intent(in) :: ITYPE + real,intent(in) :: DPBND + integer, dimension(IM,Jsta:jend),intent(in) :: L1D + real, dimension(IM,Jsta:jend),intent(in) :: P1D,T1D +! real, dimension(IM,jsta:jend),intent(inout) :: Q1D,CAPE,CINS,PPARC,ZEQL + real, dimension(IM,jsta:jend),intent(inout) :: Q1D,CAPE,CINS + real, dimension(IM,jsta:jend) :: PPARC,ZEQL + real, dimension(IM,jsta:jend),intent(inout) :: LFC,ESRHL,ESRHH + real, dimension(IM,jsta:jend),intent(inout) :: DCAPE,DGLD,ESP + integer, dimension(im,jsta:jend) ::L12,L17,L3KM +! + integer, dimension(im,jsta:jend) :: IEQL, IPTB, ITHTB, PARCEL, KLRES, KHRES, LCL, IDX +! + real, dimension(im,jsta:jend) :: THESP, PSP, CAPE20, QQ, PP, THUND + integer, dimension(im,jsta:jend) :: PARCEL2 + real, dimension(im,jsta:jend) :: THESP2,PSP2 + real, dimension(im,jsta:jend) :: CAPE4,CINS4 + REAL, ALLOCATABLE :: TPAR(:,:,:) + REAL, ALLOCATABLE :: TPAR2(:,:,:) + + 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 PRESK2, ESATP2, QSATP2, TVP2, THETAP2, TV2, THETAA2 +! real,external :: fpvsnew + integer I,J,L,KNUML,KNUMH,LBEG,LEND,IQ, KB,ITTBK + integer IE,IW,JN,JS,IVE(JM),IVW(JM),JVN,JVS + integer ISTART,ISTOP,JSTART,JSTOP + real, dimension(IM,jsta:jend) :: HTSFC + +! integer I,J,L,KNUML,KNUMH,LBEG,LEND,IQ,IT,LMHK, KB,ITTBK +! +!************************************************************** +! START CALCAPE HERE. +! + ALLOCATE(TPAR(IM,JSTA_2L:JEND_2U,LM)) + ALLOCATE(TPAR2(IM,JSTA_2L:JEND_2U,LM)) +! +! 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 +! +!$omp parallel do + DO J=JSTA,JEND + DO I=1,IM + CAPE(I,J) = D00 + CAPE20(I,J) = D00 + CAPE4(I,J) = D00 + CINS(I,J) = D00 + CINS4(I,J) = D00 + LCL(I,J) = 0 + THESP(I,J) = D00 + IEQL(I,J) = LM + PARCEL(I,J) = LM + PSP(I,J) = D00 + PPARC(I,J) = D00 + THUNDER(I,J) = .TRUE. + LFC(I,J) = D00 + ESRHL(I,J) = D00 + ESRHH(I,J) = D00 + DCAPE(I,J) = D00 + DGLD(I,J) = D00 + ESP(I,J) = D00 + THESP2(I,J) = 1000. + PSP2(I,J) = D00 + PARCEL2(I,J) = LM + ENDDO + ENDDO +! +!$omp parallel do + DO L=1,LM + DO J=JSTA,JEND + DO I=1,IM + TPAR(I,J,L) = D00 + TPAR2(I,J,L) = D00 + ENDDO + ENDDO + ENDDO +! +! FIND SURFACE HEIGHT +! + IF(gridtype == 'E')THEN + JVN = 1 + JVS = -1 + do J=JSTA,JEND + IVE(J) = MOD(J,2) + IVW(J) = IVE(J)-1 + enddo + ISTART = 2 + ISTOP = IM-1 + JSTART = JSTA_M + JSTOP = JEND_M + ELSE IF(gridtype == 'B')THEN + JVN = 1 + JVS = 0 + do J=JSTA,JEND + IVE(J)=1 + IVW(J)=0 + enddo + ISTART = 2 + ISTOP = IM-1 + JSTART = JSTA_M + JSTOP = JEND_M + ELSE + JVN = 0 + JVS = 0 + do J=JSTA,JEND + IVE(J) = 0 + IVW(J) = 0 + enddo + ISTART = 1 + ISTOP = IM + JSTART = JSTA + JSTOP = JEND + END IF +!!$omp parallel do private(htsfc,ie,iw) + IF(gridtype /= 'A') CALL EXCH(FIS(1:IM,JSTA:JEND)) + DO J=JSTART,JSTOP + DO I=ISTART,ISTOP + IE = I+IVE(J) + IW = I+IVW(J) + JN = J+JVN + JS = J+JVS +!mp PDSLVK=(PD(IW,J)*RES(IW,J)+PD(IE,J)*RES(IE,J)+ +!mp 1 PD(I,J+1)*RES(I,J+1)+PD(I,J-1)*RES(I,J-1))*0.25 +!mp PSFCK=AETA(LMV(I,J))*PDSLVK+PT + IF (gridtype=='B')THEN + HTSFC(I,J) = (0.25/g)*(FIS(IW,J)+FIS(IE,J)+FIS(I,JN)+FIS(IE,JN)) + ELSE + HTSFC(I,J) = (0.25/g)*(FIS(IW,J)+FIS(IE,J)+FIS(I,JN)+FIS(I,JS)) + ENDIF + ENDDO + ENDDO +! +! TYPE 2 CAPE/CINS: +! NOTE THAT FOR TYPE 1 CAPE/CINS ARRAYS P1D, T1D, Q1D +! ARE DUMMY ARRAYS. +! + IF (ITYPE == 2) THEN +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=1,IM + Q1D(I,J) = MIN(MAX(H1M12,Q1D(I,J)),H99999) + ENDDO + ENDDO + ENDIF +!-------FOR ITYPE=1--FIND MAXIMUM THETA E LAYER IN LOWEST DPBND ABOVE GROUND------- +!-------FOR ITYPE=2--FIND THETA E LAYER OF GIVEN T1D, Q1D, P1D--------------------- +!--------------TRIAL MAXIMUM BUOYANCY LEVEL VARIABLES------------------- + + DO KB=1,LM +!hc IF (ITYPE==2.AND.KB>1) cycle + IF (ITYPE == 1 .OR. (ITYPE == 2 .AND. KB == 1)) THEN + +!$omp parallel do private(i,j,apebtk,apespk,bqk,bqs00k,bqs10k,iq,ittbk, & +!$omp & p00k,p01k,p10k,p11k,pkl,psfck,qbtk,sqk,sqs00k, & +!$omp & sqs10k,tbtk,tpspk,tqk,tthbtk,tthesk,tthk) + DO J=JSTA,JEND + DO I=1,IM + PSFCK = PMID(I,J,NINT(LMH(I,J))) + PKL = PMID(I,J,KB) + +!hc IF (ITYPE==1.AND.(PKLPSFCK)) cycle + IF (ITYPE ==2 .OR. & + (ITYPE == 1 .AND. (PKL >= PSFCK-DPBND .AND. PKL <= PSFCK)))THEN + IF (ITYPE == 1) THEN + TBTK = T(I,J,KB) + QBTK = max(0.0, Q(I,J,KB)) + APEBTK = (H10E5/PKL)**CAPA + ELSE + PKL = P1D(I,J) + TBTK = T1D(I,J) + QBTK = max(0.0, Q1D(I,J)) + 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(I,J) = TTHK - AINT(TTHK) + ITTBK = INT(TTHK) + 1 +!--------------KEEPING INDICES WITHIN THE TABLE------------------------- + IF(ITTBK < 1) THEN + ITTBK = 1 + QQ(I,J) = D00 + ENDIF + IF(ITTBK >= JTB) THEN + ITTBK = JTB-1 + QQ(I,J) = 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(I,J) + BQS00K + SQK = (SQS10K-SQS00K)*QQ(I,J) + SQS00K + TQK = (QBTK-BQK)/SQK*RDQ + PP(I,J) = TQK-AINT(TQK) + IQ = INT(TQK)+1 +!--------------KEEPING INDICES WITHIN THE TABLE------------------------- + IF(IQ < 1) THEN + IQ = 1 + PP(I,J) = D00 + ENDIF + IF(IQ >= ITB) THEN + IQ = ITB-1 + PP(I,J) = 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(I,J) + (P01K-P00K)*QQ(I,J) & + + (P00K-P10K-P01K+P11K)*PP(I,J)*QQ(I,J) + +!!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(I,J)) THEN + PSP (I,J) = TPSPK + THESP(I,J) = TTHESK + PARCEL(I,J) = KB + ENDIF +!--------------CHECK FOR MINIMUM THETA E-------------------------------- + IF(TTHESK < THESP2(I,J)) THEN + PSP2 (I,J) = TPSPK + THESP2(I,J) = TTHESK + PARCEL2(I,J) = KB + ENDIF + END IF + ENDDO ! I loop + ENDDO ! J loop + END IF + ENDDO ! KB loop + +!----FIND THE PRESSURE OF THE PARCEL THAT WAS LIFTED +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=1,IM + PPARC(I,J) = PMID(I,J,PARCEL(I,J)) + ENDDO + ENDDO +! +!-----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)--------------------- + DO L=1,LM +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=1,IM + IF (PMID(I,J,L) < PSP(I,J)) LCL(I,J) = L+1 + ENDDO + ENDDO + ENDDO +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=1,IM + IF (LCL(I,J) > NINT(LMH(I,J))) LCL(I,J) = NINT(LMH(I,J)) + IF (ITYPE > 2) THEN + IF (T(I,J,LCL(I,J)) < 263.15) THEN + THUNDER(I,J) = .FALSE. + ENDIF + ENDIF + ENDDO + ENDDO +!----------------------------------------------------------------------- +!---------FIND TEMP OF PARCEL LIFTED ALONG MOIST ADIABAT (TPAR)--------- +!----------------------------------------------------------------------- + DO L=LM,1,-1 +!--------------SCALING PRESSURE & TT TABLE INDEX------------------------ + KNUML = 0 + KNUMH = 0 + DO J=JSTA,JEND + DO I=1,IM + KLRES(I,J) = 0 + KHRES(I,J) = 0 + IF(L <= LCL(I,J)) THEN + IF(PMID(I,J,L) < PLQ)THEN + KNUML = KNUML + 1 + KLRES(I,J) = 1 + ELSE + KNUMH = KNUMH + 1 + KHRES(I,J) = 1 + ENDIF + ENDIF + ENDDO + ENDDO +!*** +!*** COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE 0) THEN + CALL TTBLEX(TPAR(1,JSTA_2L,L),TTBL,ITB,JTB,KLRES & + , PMID(1,JSTA_2L,L),PL,QQ,PP,RDP,THE0,STHE & + , RDTHE,THESP,IPTB,ITHTB) + ENDIF +!*** +!*** COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE>PLQ +!** + IF(KNUMH > 0) THEN + CALL TTBLEX(TPAR(1,JSTA_2L,L),TTBLQ,ITBQ,JTBQ,KHRES & + , PMID(1,JSTA_2L,L),PLQ,QQ,PP,RDPQ & + ,THE0Q,STHEQ,RDTHEQ,THESP,IPTB,ITHTB) + ENDIF + +!------------SEARCH FOR EQ LEVEL---------------------------------------- +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=1,IM + IF(KHRES(I,J) > 0) THEN + IF(TPAR(I,J,L) > T(I,J,L)) IEQL(I,J) = L + ENDIF + ENDDO + ENDDO +! +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=1,IM + IF(KLRES(I,J) > 0) THEN + IF(TPAR(I,J,L) > T(I,J,L)) IEQL(I,J) = L + ENDIF + ENDDO + ENDDO +!----------------------------------------------------------------------- + ENDDO ! end of do l=lm,1,-1 loop +!------------COMPUTE CAPE AND CINS-------------------------------------- + LBEG = 1000 + LEND = 0 + DO J=JSTA,JEND + DO I=1,IM + LBEG = MIN(IEQL(I,J),LBEG) + LEND = MAX(LCL(I,J),LEND) + ENDDO + ENDDO +! +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=1,IM + IF(T(I,J,IEQL(I,J)) > 255.65) THEN + THUNDER(I,J) = .FALSE. + ENDIF + ENDDO + ENDDO +! +!reverse L order from bottom up for ESRH calculation +! + ESRHH = LCL + ESRHL = LCL +! DO L=LBEG,LEND + DO L=LEND,LBEG,-1 + +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=1,IM + IDX(I,J) = 0 + IF(L >= IEQL(I,J).AND.L <= LCL(I,J)) THEN + IDX(I,J) = 1 + ENDIF + ENDDO + ENDDO +! +!$omp parallel do private(i,j,gdzkl,presk,thetaa,thetap,esatp,qsatp,tvp,tv,& +!$omp & presk2,esatp2,qsatp2,tvp2,thetap2,tv2,thetaa2) + DO J=JSTA,JEND + DO I=1,IM + IF(IDX(I,J) > 0) THEN + PRESK = PMID(I,J,L) + GDZKL = (ZINT(I,J,L)-ZINT(I,J,L+1)) * G + ESATP = min(FPVSNEW(TPAR(I,J,L)),PRESK) + QSATP = EPS*ESATP/(PRESK-ESATP*ONEPS) +! TVP = TPAR(I,J,L)*(1+0.608*QSATP) + TVP = TVIRTUAL(TPAR(I,J,L),QSATP) + THETAP = TVP*(H10E5/PRESK)**CAPA +! TV = T(I,J,L)*(1+0.608*Q(I,J,L)) + TV = TVIRTUAL(T(I,J,L),Q(I,J,L)) + THETAA = TV*(H10E5/PRESK)**CAPA + IF(THETAP < THETAA) THEN + CINS4(I,J) = CINS4(I,J) + (LOG(THETAP)-LOG(THETAA))*GDZKL + IF(ZINT(I,J,L)-HTSFC(I,J) <= 3000.) THEN + CINS(I,J) = CINS(I,J) + (LOG(THETAP)-LOG(THETAA))*GDZKL + ENDIF + ELSEIF(THETAP > THETAA) THEN + CAPE4(I,J) = CAPE4(I,J) + (LOG(THETAP)-LOG(THETAA))*GDZKL + IF(ZINT(I,J,L)-HTSFC(I,J) <= 3000.) THEN + CAPE(I,J) = CAPE(I,J) + (LOG(THETAP)-LOG(THETAA))*GDZKL + ENDIF + IF (THUNDER(I,J) .AND. T(I,J,L) < 273.15 & + .AND. T(I,J,L) > 253.15) THEN + CAPE20(I,J) = CAPE20(I,J) + (LOG(THETAP)-LOG(THETAA))*GDZKL + ENDIF + ENDIF + +! LFC + IF (ITYPE /= 1) THEN + PRESK2 = PMID(I,J,L+1) + ESATP2 = min(FPVSNEW(TPAR(I,J,L+1)),PRESK2) + QSATP2 = EPS*ESATP2/(PRESK2-ESATP2*ONEPS) +! TVP2 = TPAR(I,J,L+1)*(1+0.608*QSATP2) + TVP2 = TVIRTUAL(TPAR(I,J,L+1),QSATP2) + THETAP2 = TVP2*(H10E5/PRESK2)**CAPA +! TV2 = T(I,J,L+1)*(1+0.608*Q(I,J,L+1)) + TV2 = TVIRTUAL(T(I,J,L+1),Q(I,J,L+1)) + THETAA2 = TV2*(H10E5/PRESK2)**CAPA + IF(THETAP >= THETAA .AND. THETAP2 <= THETAA2) THEN + IF(LFC(I,J) == D00)THEN + LFC(I,J) = ZINT(I,J,L) + ENDIF + ENDIF + ENDIF +! +! ESRH/CAPE threshold check + IF(ZINT(I,J,L)-HTSFC(I,J) <= 3000.) THEN + IF(CAPE4(I,J) >= 100. .AND. CINS4(I,J) >= -250.) THEN + IF(ESRHL(I,J) == LCL(I,J)) ESRHL(I,J)=L + ENDIF + ESRHH(I,J)=L + ENDIF + + ENDIF !(IDX(I,J) > 0) + ENDDO + ENDDO + ENDDO + +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=1,IM + IF(ESRHH(I,J) > ESRHL(I,J)) ESRHH(I,J)=IEQL(I,J) + ENDDO + ENDDO +! +! ENFORCE LOWER LIMIT OF 0.0 ON CAPE AND UPPER +! LIMIT OF 0.0 ON CINS. +! ENFORCE LFC ABOVE LCL AND BELOW EL +! +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=1,IM + CAPE(I,J) = MAX(D00,CAPE(I,J)) + CINS(I,J) = MIN(CINS(I,J),D00) +! equillibrium height + ZEQL(I,J) = ZINT(I,J,IEQL(I,J)) + LFC(I,J) = MIN(LFC(I,J),ZINT(I,J,IEQL(I,J))) + LFC(I,J) = MAX(ZINT(I,J, LCL(I,J)),LFC(I,J)) + IF (CAPE20(I,J) < 75.) THEN + THUNDER(I,J) = .FALSE. + ENDIF + IF (THUNDER(I,J)) THEN + THUND(I,J) = 1.0 + ELSE + THUND(I,J) = 0.0 + ENDIF + ENDDO + ENDDO +!------------COMPUTE DCAPE-------------------------------------- +!----------------------------------------------------------------------- +!---------FIND TEMP OF PARCEL DESCENDED ALONG MOIST ADIABAT (TPAR)--------- +!----------------------------------------------------------------------- + IF (ITYPE == 1) THEN + + DO L=LM,1,-1 +!--------------SCALING PRESSURE & TT TABLE INDEX------------------------ + KNUML = 0 + KNUMH = 0 + DO J=JSTA,JEND + DO I=1,IM + KLRES(I,J) = 0 + KHRES(I,J) = 0 + PSFCK = PMID(I,J,NINT(LMH(I,J))) + PKL = PMID(I,J,L) + IF(PKL >= PSFCK-DPBND) THEN + IF(PMID(I,J,L) < PLQ)THEN + KNUML = KNUML + 1 + KLRES(I,J) = 1 + ELSE + KNUMH = KNUMH + 1 + KHRES(I,J) = 1 + ENDIF + ENDIF + ENDDO + ENDDO +!*** +!*** COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE 0) THEN + CALL TTBLEX(TPAR2(1,JSTA_2L,L),TTBL,ITB,JTB,KLRES & + , PMID(1,JSTA_2L,L),PL,QQ,PP,RDP,THE0,STHE & + , RDTHE,THESP2,IPTB,ITHTB) + ENDIF +!*** +!*** COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE>PLQ +!** + IF(KNUMH > 0) THEN + CALL TTBLEX(TPAR2(1,JSTA_2L,L),TTBLQ,ITBQ,JTBQ,KHRES & + , PMID(1,JSTA_2L,L),PLQ,QQ,PP,RDPQ & + , THE0Q,STHEQ,RDTHEQ,THESP2,IPTB,ITHTB) + ENDIF + ENDDO ! end of do l=lm,1,-1 loop + + LBEG = 1 + LEND = LM + + DO L=LBEG,LEND +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=1,IM + IDX(I,J) = 0 + IF(L >= PARCEL2(I,J).AND.L < NINT(LMH(I,J))) THEN + IDX(I,J) = 1 + ENDIF + ENDDO + ENDDO +! +!$omp parallel do private(i,j,gdzkl,presk,thetaa,thetap,esatp,qsatp,tvp,tv) + DO J=JSTA,JEND + DO I=1,IM + IF(IDX(I,J) > 0) THEN + PRESK = PMID(I,J,L) + GDZKL = (ZINT(I,J,L)-ZINT(I,J,L+1)) * G + ESATP = min(FPVSNEW(TPAR2(I,J,L)),PRESK) + QSATP = EPS*ESATP/(PRESK-ESATP*ONEPS) +! TVP = TPAR2(I,J,L)*(1+0.608*QSATP) + TVP = TVIRTUAL(TPAR2(I,J,L),QSATP) + THETAP = TVP*(H10E5/PRESK)**CAPA +! TV = T(I,J,L)*(1+0.608*Q(I,J,L)) + TV = TVIRTUAL(T(I,J,L),Q(I,J,L)) + THETAA = TV*(H10E5/PRESK)**CAPA + !IF(THETAP < THETAA) THEN + DCAPE(I,J) = DCAPE(I,J) + (LOG(THETAP)-LOG(THETAA))*GDZKL + !ENDIF + ENDIF + ENDDO + ENDDO + ENDDO + +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=1,IM + DCAPE(I,J) = MIN(D00,DCAPE(I,J)) + ENDDO + ENDDO + + ENDIF !ITYPE=1 FOR DCAPE + +! +! Dendritic Growth Layer depth +! the layer with temperatures from -12 to -17 C in meters +! + L12=LM + L17=LM + DO L=LM,1,-1 +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=1,IM + IF(T(I,J,L) <= TFRZ-12. .AND. L12(I,J)==LM) L12(I,J)=L + IF(T(I,J,L) <= TFRZ-17. .AND. L17(I,J)==LM) L17(I,J)=L + ENDDO + ENDDO + ENDDO +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=1,IM + IF(L12(I,J)/=LM .AND. L17(I,J)/=LM) THEN + DGLD(I,J)=ZINT(I,J,L17(I,J))-ZINT(I,J,L12(I,J)) + DGLD(I,J)=MAX(DGLD(I,J),0.) + ENDIF + ENDDO + ENDDO +! +! Enhanced Stretching Potential +! ESP = (0-3 km MLCAPE / 50 J kg-1) * ((0-3 km lapse rate - 7.0) / 1.0 C (km-1) +! https://www.spc.noaa.gov/exper/soundings/help/params4.html +! + L3KM=LM + DO L=LM,1,-1 +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=1,IM + IF(ZINT(I,J,L)-HTSFC(I,J) <= 3000.) L3KM(I,J)=L + ENDDO + ENDDO + ENDDO +!$omp parallel do private(i,j) + DO J=JSTA,JEND + DO I=1,IM + ESP(I,J) = (CAPE(I,J) / 50.) * (T(I,J,LM) - T(I,J,L3KM(I,J)) - 7.0) + IF((T(I,J,LM) - T(I,J,L3KM(I,J))) < 7.0) ESP(I,J) = 0. +! IF(CAPE(I,J) < 250.) ESP(I,J) = 0. + ENDDO + ENDDO +! + DEALLOCATE(TPAR) + DEALLOCATE(TPAR2) +! + END SUBROUTINE CALCAPE2 +! +!------------------------------------------------------------------------------------- +! + elemental function TVIRTUAL(T,Q) +! +! COMPUTE VIRTUAL TEMPERATURE +! + IMPLICIT NONE + REAL TVIRTUAL + REAL, INTENT(IN) :: T, Q + + TVIRTUAL = T*(1+0.608*Q) + + end function TVIRTUAL +! +!------------------------------------------------------------------------------------- +! + end module upp_physics diff --git a/sorc/ncep_post.fd/makefile_module b/sorc/ncep_post.fd/makefile_module index 5678f1ecd..5b6f2c763 100644 --- a/sorc/ncep_post.fd/makefile_module +++ b/sorc/ncep_post.fd/makefile_module @@ -72,10 +72,11 @@ OBJS = wrf_io_flags.o getVariable.o getIVariableN.o \ CMASSI.o CTLBLK.o GRIDSPEC.o LOOKUP.o PARAMR.o RHGRD.o RQSTFLD.o xml_perl_data.o \ cuparm.o params.o svptbl.o get_postfilename.o grib2_module.o \ SET_LVLSXML.o FILL_PSETFLD.o \ - BNDLYR.o BOUND.o CALCAPE.o CALDWP.o CALDRG.o CALHEL.o CALLCL.o \ - CALMCVG.o CALPOT.o CALPW.o CALRH.o CALRCH.o CALRH_GSD.o \ + UPP_MATH.o UPP_PHYSICS.o \ + BNDLYR.o BOUND.o CALDWP.o CALDRG.o CALHEL.o CALLCL.o \ + CALMCVG.o CALPOT.o CALPW.o CALRCH.o \ CALSTRM.o CALTAU.o CALTHTE.o CALVIS.o CALVIS_GSD.o CALVOR.o CALWXT.o \ - CALWXT_RAMER.o CALWXT_BOURG.o CALWXT_REVISED.o CALRH_PW.o \ + CALWXT_RAMER.o CALWXT_BOURG.o CALWXT_REVISED.o \ CALWXT_EXPLICIT.o CALWXT_DOMINANT.o \ CLDRAD.o CLMAX.o COLLECT.o COLLECT_LOC.o DEWPOINT.o \ FDLVL.o FGAMMA.o FIXED.o FRZLVL.o FRZLVL2.o \ @@ -89,14 +90,14 @@ OBJS = wrf_io_flags.o getVariable.o getIVariableN.o \ CALMICT.o MICROINIT.o GPVS.o MDL2SIGMA.o \ ETCALC.o CANRES.o CALGUST.o WETFRZLVL.o SNFRAC.o MDL2AGL.o SNFRAC_GFS.o \ AVIATION.o DEALLOCATE.o \ - CALPBL.o MDL2SIGMA2.o INITPOST_GFS.o CALRH_GFS.o LFMFLD_GFS.o \ + CALPBL.o MDL2SIGMA2.o INITPOST_GFS.o LFMFLD_GFS.o \ CALRAD_WCLOUD_newcrtm.o MDL2THANDPV.o CALPBLREGIME.o POLEAVG.o \ INITPOST_NEMS.o GETNEMSNDSCATTER.o ICAOHEIGHT.o INITPOST_GFS_NEMS.o \ - GEO_ZENITH_ANGLE.o GFIP3.o GRIDAVG.o CALUPDHEL.o INITPOST_GFS_SIGIO.o \ + GEO_ZENITH_ANGLE.o GFIP3.o CALUPDHEL.o INITPOST_GFS_SIGIO.o \ AllGETHERV_GSD.o MSFPS.o SELECT_CHANNELS.o ALLOCATE_ALL.o INITPOST_NEMS_MPIIO.o ASSIGNNEMSIOVAR.o \ INITPOST_GFS_NEMS_MPIIO.o INITPOST_NETCDF.o INITPOST_GFS_NETCDF.o INITPOST_GFS_NETCDF_PARA.o \ gtg_ctlblk.o gtg_indices.o gtg_filter.o gtg_compute.o gtg_config.o map_routines.o gtg_algo.o gtg_smoothseams.o CALVESSEL.o \ - CALHEL2.o CALCAPE2.o ETAMP_Q2F.o + CALHEL2.o ETAMP_Q2F.o .SUFFIXES: .F .f .o .f90 .c