SUBROUTINE CUCNVC !>-------------------------------------------------------------------------------------------------- !> SUBROUTINE CUCNVC !> !> SUBPROGRAM: CUCNVC - CONVECTIVE PRECIPITATION PARAMETERIZATION !> PROGRAMMER: JANJIC !> ORG: W/NP2 !> DATE: 93-11-02 !> !> ABSTRACT: !> CUCNVC CALCULATES THE SUB-GRID SCALE CONVECTION INCLUDING DEEP AND SHALLOW CONVECTIVE CLOUDS !> FOLLOWING THE SCHEME DESCRIBED BY JANJIC (1994) BUT WITH SIGNIFICANT MODIFICATIONS. !> IN ADDITION, THE LATENT HEAT RELEASE AND MOISTURE CHANGE DUE TO PRECIPITATING AND NON - !> PRECIPITATING CLOUDS ARE COMPUTED. !> !> PROGRAM HISTORY LOG: !> 87-09-?? JANJIC - ORIGINATOR !> 90-11-21 JANJIC - TWO SETS OF DSP PROFILES (FAST AND SLOW) REPLACE THE ORIGINAL ONE SET !> 95-03-25 BLACK - CONVERSION FROM 1-D TO 2-D IN HORIZONTAL !> 96-03-28 BLACK - ADDED EXTERNAL EDGE !> 98-11-02 BLACK - MODIFIED FOR DISTRIBUTED MEMORY !> 18-01-15 LUCCI - MODERNIZATION OF THE CODE, INCLUDING: !> * F77 TO F90/F95 !> * INDENTATION & UNIFORMIZATION CODE !> * REPLACEMENT OF COMMONS BLOCK FOR MODULES !> * DOCUMENTATION WITH DOXYGEN !> * OPENMP FUNCTIONALITY !> 26-06-19 DIEGO - INCLUDING PCC SCHEME (KOH AND FONSECA 2016) !> FOR CUMULUS-RADIATION FEEDBACK !> !> INPUT ARGUMENT LIST: !> NONE !> !> OUTPUT ARGUMENT LIST: !> NONE !> !> INPUT/OUTPUT ARGUMENT LIST: !> NONE !> !> USE MODULES: ACMCLH !> CNVCLD !> CTLBLK !> CUPARM !> DYNAM0 !> F77KINDS !> GLB_TABLE !> INDX !> LOOPS !> MAPPINGS !> MASKS !> MPPCOM !> PARMETA !> PARM_TBL !> PHYS0 !> PPTASM !> PVRBLS !> TEMPCOM !> TOPO !> VRBLS !> !> DRIVER : EBU !> NEWFLT !> !> CALLS : TTBLEX !> ZERO2 !>-------------------------------------------------------------------------------------------------- ! !----------------------------------------------------------------------------------------- ! REFERENCES: ! ! JANJIC, Z.I., 1994: THE STEP-MOUNTAIN ETA COORDINATE MODEL: ! FURTHER DEVELOPMENTS OF THE CONVECTION, VISCOUS SUBLAYER AND TURBULENCE CLOSURE SCHEMES. ! MONTHLY WEATHER REVIEW, VOL. 122, 927-945. !----------------------------------------------------------------------------------------- ! !----------------------------------------------------------------------- ! WARNING: THIS SUBROUTINE WILL NOT WORK IF (LM .LT. 12); ! MUST BE CALLED IN THE SAME STEP WITH PROFQ2 BECAUSE PROFQ DEFINES APE; !----------------------------------------------------------------------- USE ACMCLH USE CNVCLD USE CTLBLK USE CUPARM USE F77KINDS USE GLB_TABLE USE INDX USE LOOPS USE MAPPINGS USE MASKS USE MPPCOM USE PARMETA USE PARM_TBL USE PHYS0 USE PPTASM USE PVRBLS USE TEMPCOM USE TOPO USE VRBLS !-------diego---- USE RRTMCOMMS !-------diego---- ! IMPLICIT NONE ! LOGICAL(KIND=L4KIND) ::& & UNIS , UNIL , OCT90 ! NAMELIST /CUPARMDATA/ & & STABS , STABD , STABFC , DTTOP , & & RHF , EPSUP , EPSDN , EPSTH , & & PBM , PQM , PNO , PONE , PSH , & & PFRZ , PSHU , & & UNIS , UNIL , OCT90 & & FSS , EFIMN , EFMNT , FCC , & & DSPBFL , DSP0FL , DSPTFL , FSL , & & DSPBFS , DSP0FS , DSPTFS , & & TREL , EPSNTP , EFIFC , & & DSPC , EPSP , & & STEFI !------------------------ ! IMPLICIT NONE VARIABLES !------------------------ ! !------------------------------ ! INSTABILITY FOR TOO LARGE LSH !------------------------------ INTEGER(KIND=I4KIND), PARAMETER :: KSMUD = 0 INTEGER(KIND=I4KIND), PARAMETER :: NROW = 0 !-------------------------------------------------------------- ! INSTABILITY FOR TOO LARGE LSH ! ! LP1=LM+1,LM1=LM-1,LNO=1,LSH=LM/3-1,LSHU=LM/2-1,LQM=LM/5,KBM=3 ! LP1=LM+1,LM1=LM-1,LNO=1,LSH=LM/3 ,LSHU=LM/2-1,LQM=LM/5,KBM=3 ! LP1=LM+1,LM1=LM-1,LNO=3,LSH=LM/3 ,LSHU=LM/2-1,LQM=LM/5,KBM=3 ! LP1=LM+1,LM1=LM-1,LNO=2,LSH=LM/3-1,LSHU=LM/2-1,LQM=LM/5,KBM=3 ! LP1=LM+1,LM1=LM-1,LNO=2,LSH=LM/3-2,LSHU=LM/2-1,LQM=LM/5,KBM=3 !-------------------------------------------------------------- INTEGER(KIND=I4KIND), PARAMETER :: IMJM = IM * JM - JM / 2 INTEGER(KIND=I4KIND), PARAMETER :: IMJM_LOC = IDIM2 * JDIM2 ! REAL (KIND=R4KIND), DIMENSION(LM) ::& & TREFK , QREFK , PK , APEK , TK , & & THSK , PSK , APESK , QK , THERK , & & THVREF , THEVRF , THVMOD , DIFT , DIFQ , & & QSATK , FPK ! INTEGER(KIND=I4KIND), DIMENSION(LM) ::& & NTOPD , NBOTD , NTOPS , NBOTS , & & NDPTHD , NDPTHS ! INTEGER(KIND=I4KIND), DIMENSION(IDIM1:IDIM2, JDIM1:JDIM2) ::& & LTOP , LBOT , & & IPTB , ITHTB ! REAL (KIND=R4KIND), DIMENSION(IDIM1:IDIM2, JDIM1:JDIM2) ::& & PTOP , PBOT , & & PDSL , APEBT , & & TBT , Q2BT , & & QQ , PP , & & PSP , THBT , & & THESP , P , & & BTH , STH , & & T00 , T10 , & & T01 , T11 , & & WF1 , WF2 , & & WF3 , WF4 , & & PRECOL ! INTEGER(KIND=I4KIND), DIMENSION(IMJM_LOC) ::& & IBUOY , JBUOY , & & IDEEP , JDEEP , & & ISHAL , JSHAL , & & ILRES , JLRES , & & IHRES , JHRES ! REAL (KIND=R4KIND), DIMENSION(IDIM1:IDIM2, JDIM1:JDIM2, LM) ::& & APE , & & TREF , & & TMOD , & & QMOD ! REAL (KIND=R4KIND), DIMENSION(IDIM1:IDIM2, JDIM1:JDIM2) ::& & DSPB , DSP0 , & & DSPT , & & TL , QL , & & TNE , TSE , & & QNE , QSE ! REAL (KIND=R4KIND) ::& & STABS , STABD , STABFC , DTTOP , RHF , EPSUP , EPSDN , EPSTH , PBM , & & PSHU , FSS , EFIMN , EFMNT , FCC , DSPBFL , DSP0FL , DSPTFL , FSL , & & DSPBFS , DSP0FS , DSPTFS , TREL , PQM , PNO , PONE , PFRZ , PSH , & & EPSNTP , EFIFC , DSPC , EPSP , & & STEFI , FCB , DSPBSL , DSP0SL , & & DSPTSL , DSPBSS , DSP0SS , DSPTSS , & & AVGEFI , SLOPBL , SLOP0L , SLOPTL , & & SLOPBS , SLOP0S , SLOPTS , SLOPE , & & A23M4L , ELOCP , CPRLG , RCP , & & DTCNVC , RDTCNVC , TAUK , APESTS , PKL , & & PSFCK , QBT , TTHBT , TTH , QQ1 , BQS00K , & & SQS00K , BQS10K , SQS10K , BQ , SQ , TQ , PP1 , P00K , & & P10K , P01K , P11K , TPSP , & & APESP , TTHES , AETAL , PRESK , EFI , SMK , RSMK , & & PSHNEW , PSFCIJ , DEPMIN , DEPTH , & & DSPBK , DSP0K , DSPTK , TKL , QKL , APEKL , PKB , PKT , & & PK0 , TREFKX , THERKX , APEKXX , & & THERKY , APEKXY , STABDL , RDP0T , & & DTHEM , DEPWL , DSP , SUMDP , SUMDE , HCORR , TSKL , DHDT , & & THSKL , DENTPY , AVRGT , PRECK , DIFTL , DIFQL , AVRGTL , PBTK , & & PTPK , RHH , RHMAX , RHL , DRHDP , DRHEAT , FEFI , PZ0 , & & THVMKL , THTPK , TTHK , QQK , BQK , SQK , TQK , PPK , & & PART1 , PART2 , PART3 , DPMIX , SMIX , PKXXXX , SUMDT , RDPSUM , & & TCORR , PKXXXY , TRFKL , PSUM , QSUM , POTSUM , QOTSUM , OTSUM , & & DST , FPTK , DPKL , RTBAR , ROTSUM , DSTQ , DEN , DQREF , & & QRFTP , QRFKL , QNEW , DTDETA , CUTOP , CUBOT ! INTEGER(KIND=I4KIND) ::& & I , J , K , KB , NDSTN , NDSTP , LLMH , & & LMHK , IQTB , IQ , IT , NSHAL , NEGDS , NSMUD , & & LMHIJ , KNUML , KNUMH , KROW , KS , & & KOUNT , KHDEEP , N , LTPK , ITTB , ITTBK , & & LBTK , LB , LTP1 , LBM1 , & & L0 , L0M1 , ITER , LCOR , & & LQM , LTSH , LSHU , LTP2 , & & NDEEP , NDEPTH , NNEG , KHSHAL ! !---------------diego------------------------------------------------------------------------------ ! Acordingly to Koh and Fonseca 2016 (PCC Scheme) ! INTEGER :: KFLIP, TTOP, BBOT, JPRINT ! REAL (KIND=R4KIND), PARAMETER :: RD = 287.04 REAL (KIND=R4KIND), PARAMETER :: RV = 461.5 REAL (KIND=R4KIND), PARAMETER :: EPS1 = RV / RD - 1. REAL (KIND=R4KIND), PARAMETER :: CRMX = 85.0 REAL (KIND=R4KIND), PARAMETER :: CRMN = 0.14 REAL (KIND=R4KIND), PARAMETER :: PEPS = 1./2.5 REAL (KIND=R4KIND), PARAMETER :: GAM = 0.5 !-------------------------------------------------- ! These information must be at /ucl/cuparmdata.dat !-------------------------------------------------- !diego REAL (KIND=R4KIND), PARAMETER :: EFIMN = 0.20 !diego REAL (KIND=R4KIND), PARAMETER :: EFMNT = 0.70 !diego REAL (KIND=R4KIND), PARAMETER :: FCC = 5.00 !diego REAL (KIND=R4KIND), PARAMETER :: FSS = 0.60 !diego REAL (KIND=R4KIND), PARAMETER :: FSL = 0.60 !------------------------------------------------------------------ ! !diego REAL (KIND=R4KIND), DIMENSION(IDIM1:IDIM2, JDIM1:JDIM2), INTENT(INOUT) ::& REAL (KIND=R4KIND), DIMENSION(IDIM1:IDIM2, JDIM1:JDIM2) ::& & CONVCLD ! REAL (KIND=R4KIND), DIMENSION(IDIM1:IDIM2, JDIM1:JDIM2, LM) ::& & RQVSHCUTEN , RTHSHCUTEN, RHO , DZ8W , PMID , THTE ! !diego REAL (KIND=R4KIND), DIMENSION(IDIM1:IDIM2, JDIM1:JDIM2, LM), INTENT(INOUT) ::& !diego & CCLDFRA , QCCONV , QICONV ! ! LOCAL VARIABLES ! REAL (KIND=R4KIND) ::& & PAVG , PPAVG , PWCOL , DQCOL , DQCOLMIN, CUMX , QCIS , GAMMA , & & RRP , PRRT , MCOL , MPVPR , FACTL , PDQD , PRWD !diego & RRP , PRRT , MCOL , MPVPR , BBOT , TTOP , FACTL , PDQD , PRWD ! REAL (KIND=R4KIND), DIMENSION(LM) ::& & DPCOL , DQDT , DTDT , PCOL , QCOL , TCOL , DQDTS , DTDTS , & & PVPR , JPR , QCC , QIC , CCLDN , CCLD , CCLDOLD , & & T_COL , QV_COL ! REAL (KIND=R4KIND), DIMENSION(IDIM1:IDIM2, JDIM1:JDIM2) ::& & PRATEC !---------------diego------------------------------------------------------------------------------ ! OPEN(11,FILE='CUPARMDATA.dat',STATUS='OLD') REWIND 11 READ(11,CUPARMDATA) CLOSE(11) ! FCB = 1. - FCC DSPBSL = DSPBFL * FSL DSP0SL = DSP0FL * FSL DSPTSL = DSPTFL * FSL DSPBSS = DSPBFS * FSS DSP0SS = DSP0FS * FSS DSPTSS = DSPTFS * FSS ! AVGEFI = (EFIMN + 1.) * .5 ! SLOPBL = (DSPBFL-DSPBSL) / (1.-EFIMN) SLOP0L = (DSP0FL-DSP0SL) / (1.-EFIMN) SLOPTL = (DSPTFL-DSPTSL) / (1.-EFIMN) SLOPBS = (DSPBFS-DSPBSS) / (1.-EFIMN) SLOP0S = (DSP0FS-DSP0SS) / (1.-EFIMN) SLOPTS = (DSPTFS-DSPTSS) / (1.-EFIMN) ! SLOPE = (1. - EFMNT) / (1. - EFIMN) A23M4L = A2 * (A3 - A4) * ELWV ELOCP = ELIVW / CP CPRLG = CP / (ROW * G * ELWV) RCP = 1. / CP ! CALL ZERO2(DSP0) CALL ZERO2(DSPB) CALL ZERO2(DSPT) CALL ZERO2(PSP) ! AVCNVC = AVCNVC + 1. ACUTIM = ACUTIM + 1. DTCNVC = NCNVC * DT RDTCNVC = 1. / DTCNVC TAUK = DTCNVC / TREL !---------------------------------------------- ! POSSIBLE FUTURE CHANGE FOR SHALLOW CONVECTION ! TAUKSC=DTCNVC/(5.*TREL) !---------------------------------------------- ! !------------ ! DIAGNOSTICS !------------ DO K=1,LM NTOPD(K) = 0 NBOTD(K) = 0 NTOPS(K) = 0 NBOTS(K) = 0 NDPTHS(K) = 0 NDPTHD(K) = 0 END DO !------------ !PREPARATIONS !------------ ! !------- ! OPENMP !------- !$omp parallel do DO 120 J=MYJS,MYJE DO 120 I=MYIS,MYIE LBOT(I,J) = LMH(I,J) THESP(I,J) = 0. PDSL (I,J) = RES(I,J) * PD(I,J) PRECOL(I,J) = 0. TREF(I,J,1) = T(I,J,1) !----diego------------------ !diego DO K=1,LM !diego JPR(K) = 0. !diego PVPR(K) = 0. !diego QCCONV(I,J,K) = 0. !diego QICONV(I,J,K) = 0. !diego CCLDFRA(I,K,J) = 0. !diego END DO !----diego------------------- ! 120 END DO !--------------------------------------- ! PADDING SPECIFIC HUMIDITY IF TOO SMALL ! RESTORE APE TO SCRATCH ARRAY !--------------------------------------- ! !------- ! OPENMP !------- !$omp parallel do private(APESTS) DO 130 K=1,LM DO J=MYJS,MYJE DO I=MYIS,MYIE APESTS = PDSL(I,J) * AETA(K) + PT APE(I,J,K) = (1.E5 / APESTS) ** CAPA IF (Q(I,J,K) < EPSQ) Q(I,J,K) = HTM(I,J,K) * EPSQ END DO END DO 130 END DO !---------------------------------- ! SEARCH FOR MAXIMUM BUOYANCY LEVEL !---------------------------------- DO 170 KB=1,LM !--------------------------------------- ! TRIAL MAXIMUM BUOYANCY LEVEL VARIABLES !--------------------------------------- ! !------- ! OPENMP !------- !$omp parallel do !$omp private (APESP , BQ , BQS00K, BQS10K, IQ , IT , ITTB , ITTBK , IQTB , LMHK , & !$omp P00K , P01K , P10K , P11K , PKL , PP1 , PSFCK , QBT , QQ1 , SQ , & !$omp SQS00K, SQS10K, TPSP , TQ , TTH , TTHBT , TTHES) ! DO 155 J=MYJS2,MYJE2 DO 150 I=MYIS1,MYIE1 ! PKL = AETA(KB) * PDSL(I,J) + PT LMHK = LMH(I,J) PSFCK = AETA(LMHK) * PDSL(I,J) + PT !-------------------------------------------------------------------------------------------------- ! NOW SEARCHING OVER A SCALED DEPTH IN FINDING THE PARCEL WITH THE MAX THETA-E INSTEAD OF THE OLD ! 130 MB !-------------------------------------------------------------------------------------------------- IF (KB <= LMHK .AND. PKL >= 0.80*PSFCK) THEN QBT = Q(I,J,KB) TTHBT = T(I,J,KB) * APE(I,J,KB) TTH = (TTHBT-THL) * RDTH QQ1 = TTH - AINT(TTH) ITTB = INT(TTH) + 1 !--------------------------------- ! KEEPING INDICES WITHIN THE TABLE !--------------------------------- IF (ITTB < 1) THEN ITTB = 1 QQ1 = 0. END IF IF (ITTB >= JTB) THEN ITTB = JTB - 1 QQ1 = 0. END IF CONTINUE !------------------------------------------- ! BASE AND SCALING FACTOR FOR SPEC. HUMIDITY !------------------------------------------- ITTBK = ITTB BQS00K = QS0(ITTBK) SQS00K = SQS(ITTBK) BQS10K = QS0(ITTBK+1) SQS10K = SQS(ITTBK+1) !--------------------------------------- ! SCALING SPEC. HUMIDITY AND TABLE INDEX !--------------------------------------- BQ = (BQS10K-BQS00K) * QQ1 + BQS00K SQ = (SQS10K-SQS00K) * QQ1 + SQS00K TQ = (QBT-BQ) / SQ * RDQ PP1 = TQ - AINT(TQ) IQTB = INT(TQ) + 1 !--------------------------------- ! KEEPING INDICES WITHIN THE TABLE !--------------------------------- IF (IQTB < 1) THEN IQTB = 1 PP1 = 0. END IF IF (IQTB >= ITB) THEN IQTB = ITB - 1 PP1 = 0. END IF !--------------------------------------------------- ! SATURATION PRESSURE AT FOUR SURROUNDING TABLE PTS. !--------------------------------------------------- IQ = IQTB IT = ITTB P00K = PTBL(IQ ,IT ) P10K = PTBL(IQ+1,IT ) P01K = PTBL(IQ ,IT+1) P11K = PTBL(IQ+1,IT+1) !----------------------------------------- ! SATURATION POINT VARIABLES AT THE BOTTOM !----------------------------------------- TPSP = P00K + (P10K-P00K) * PP1 + (P01K-P00K) * QQ1 + (P00K-P10K-P01K+P11K) & & * PP1 * QQ1 APESP = (1.E5/TPSP) ** CAPA TTHES = TTHBT * EXP(ELOCP*QBT*APESP/TTHBT) !--------------------------- ! CHECK FOR MAXIMUM BUOYANCY !--------------------------- IF (TTHES > THESP(I,J)) THEN PSP(I,J) = TPSP THBT(I,J) = TTHBT THESP(I,J) = TTHES END IF END IF 150 END DO 155 END DO 170 END DO !------------------------------------------------ ! CHOOSE CLOUD BASE AS MODEL LEVEL JUST BELOW PSP !------------------------------------------------ DO 240 K=1,LM1 AETAL = AETA(K) !------- ! OPENMP !------- !$omp parallel do DO J=MYJS2,MYJE2 DO I=MYIS,MYIE P(I,J) = PDSL(I,J) * AETAL + PT IF (P(I,J) < PSP(I,J) .AND. P(I,J) >= PQM) LBOT(I,J) = K + 1 END DO END DO ! 240 END DO !-------------------------------------------------------------- ! WARNING: LBOT MUST NOT BE GT LMH(I,J)-1 IN SHALLOW CONVECTION ! MAKE SURE CLOUD BASE IS AT LEAST PONE ABOVE THE SURFACE !-------------------------------------------------------------- ! !------- ! OPENMP !------- !$omp parallel do private (LMHIJ, PSFCK) DO 250 J=MYJS2,MYJE2 DO 250 I=MYIS,MYIE LMHIJ = LMH(I,J) PBOT(I,J) = AETA(LBOT(I,J)) * PDSL(I,J) + PT PSFCK = AETA(LMHIJ) * PDSL(I,J) + PT IF (PBOT(I,J) >= PSFCK-PONE .OR. LBOT(I,J) >= LMHIJ) THEN ! DO K=1,LMHIJ-1 P(I,J) = AETA(K) * PDSL(I,J) + PT IF (P(I,J) < PSFCK-PONE) LBOT(I,J) = K END DO ! PBOT(I,J) = AETA(LBOT(I,J)) * PDSL(I,J) + PT END IF 250 END DO !---------------------- ! CLOUD TOP COMPUTATION !---------------------- ! !------- ! OPENMP !------- !$omp parallel do DO J=MYJS,MYJE DO I=MYIS,MYIE LTOP(I,J) = LBOT(I,J) PTOP(I,J) = PBOT(I,J) END DO END DO !------- ! OPENMP !------- !$omp parallel do !$omp private (IHRES, ILRES, IPTB, ITHTB, JHRES, JLRES, KNUMH, KNUML, PP, PRESK, QQ) ! DO 290 K=LM,1,-1 !------------------------------------ ! SCALING PRESSURE AND TT TABLE INDEX !------------------------------------ KNUML=0 KNUMH=0 ! DO 270 J=MYJS2,MYJE2 DO 270 I=MYIS1,MYIE1 PRESK = PDSL(I,J) * AETA(K) + PT IF (PRESK < PLQ) THEN KNUML = KNUML + 1 ILRES(KNUML) = I JLRES(KNUML) = J ELSE KNUMH = KNUMH + 1 IHRES(KNUMH) = I JHRES(KNUMH) = J END IF 270 END DO !--------------------------------------------------------------- ! COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE 0) THEN CALL TTBLEX(TREF(IDIM1,JDIM1,K) & & , TTBL, ITB, JTB, KNUML, ILRES, JLRES, PDSL & & , AETA(K) & & , HTM(IDIM1,JDIM1,K) & & , PT, PL & & , QQ(IDIM1,JDIM1), PP(IDIM1,JDIM1) & & , RDP, THE0, STHE, RDTHE & & , THESP(IDIM1,JDIM1), IPTB(IDIM1,JDIM1) & & , ITHTB(IDIM1,JDIM1)) END IF !--------------------------------------------------------------- ! COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE>PL !--------------------------------------------------------------- IF (KNUMH > 0) THEN CALL TTBLEX(TREF(IDIM1,JDIM1,K) & & , TTBLQ , ITBQ , JTBQ , KNUMH & & , IHRES , JHRES , PDSL & & , AETA(K) & & , HTM(IDIM1,JDIM1,K) & & , PT , PLQ & & , QQ(IDIM1,JDIM1), PP(IDIM1,JDIM1) & & , RDPQ , THE0Q , STHEQ, RDTHEQ & & , THESP(IDIM1,JDIM1), IPTB(IDIM1,JDIM1) & & , ITHTB(IDIM1,JDIM1)) END IF 290 END DO !--------------- ! BUOYANCY CHECK !--------------- DO 295 K=LM,1,-1 ! !------- ! OPENMP !------- !$omp parallel do DO J=MYJS2,MYJE2 DO I=MYIS1,MYIE1 IF (TREF(I,J,K) > T(I,J,K)-DTTOP) LTOP(I,J) = K END DO END DO 295 END DO !------------------- ! CLOUD TOP PRESSURE !------------------- ! !------- ! OPENMP !------- !$omp parallel do DO J=MYJS2,MYJE2 DO I=MYIS1,MYIE1 PTOP(I,J) = AETA(LTOP(I,J)) * PDSL(I,J) + PT END DO END DO !---------------------------------- ! DEFINE AND SMOOTH DSPS AND CLDEFI ! UNIFIED OR SEPARATE LAND/SEA CONV !---------------------------------- IF (UNIS) THEN ! !------- ! OPENMP !------- !$omp parallel do private (EFI) DO J=MYJS,MYJE DO I=MYIS,MYIE EFI = CLDEFI(I,J) DSPB(I,J) = (EFI-EFIMN) * SLOPBS + DSPBSS DSP0(I,J) = (EFI-EFIMN) * SLOP0S + DSP0SS DSPT(I,J) = (EFI-EFIMN) * SLOPTS + DSPTSS END DO END DO ELSE IF( .NOT. UNIL) THEN ! !------- ! OPENMP !------- !$omp parallel do private (EFI) DO J=MYJS,MYJE DO I=MYIS,MYIE EFI=CLDEFI(I,J) DSPB(I,J) = ((EFI-EFIMN)*SLOPBS+DSPBSS) * SM(I,J) + ((EFI-EFIMN)*SLOPBL+DSPBSL) & & * (1.-SM(I,J)) ! DSP0(I,J) = ((EFI-EFIMN)*SLOP0S+DSP0SS) * SM(I,J) + ((EFI-EFIMN)*SLOP0L+DSP0SL) & & * (1.-SM(I,J)) ! DSPT(I,J) = ((EFI-EFIMN)*SLOPTS+DSPTSS) * SM(I,J) + ((EFI-EFIMN)*SLOPTL+DSPTSL) & & * (1.-SM(I,J)) END DO END DO ELSE ! !------- ! OPENMP !------- !$omp parallel do private (EFI) DO J=MYJS,MYJE DO I=MYIS,MYIE EFI=CLDEFI(I,J) DSPB(I,J) = ((EFI-EFIMN) * SLOPBL + DSPBSL) DSP0(I,J) = ((EFI-EFIMN) * SLOP0L + DSP0SL) DSPT(I,J) = ((EFI-EFIMN) * SLOPTL + DSPTSL) END DO END DO END IF !------------------------------------------------ ! EXTENDING SEA STRUCTURES INLAND ALONG COASTLINE !------------------------------------------------ IF (NROW > 0 .AND. .NOT. UNIS .AND. .NOT. UNIL) THEN ! !------- ! OPENMP !------- !$omp parallel do DO J=MYJS,MYJE DO I=MYIS,MYIE WF1(I,J) = 0. WF2(I,J) = 0. WF3(I,J) = 0. WF4(I,J) = 0. END DO END DO ! KROW = NROW ! DO 350 KOUNT=1,KROW ! !------- ! OPENMP !------- !$omp parallel do DO 345 J=MYJS2,MYJE2 DO 345 I=MYIS1,MYIE1 WF1(I,J) = (DSPB(I+IHW(J),J-1) + DSPB(I+IHE(J),J-1) & & + DSPB(I+IHW(J),J+1) + DSPB(I+IHE(J),J+1) & & + 4. * DSPB(I,J)) * HBM2(I,J) * 0.125 ! WF2(I,J) = (DSP0(I+IHW(J),J-1) + DSP0(I+IHE(J),J-1) & & + DSP0(I+IHW(J),J+1) + DSP0(I+IHE(J),J+1) & & + 4. * DSP0(I,J)) * HBM2(I,J) * 0.125 ! WF3(I,J) = (DSPT(I+IHW(J),J-1) + DSPT(I+IHE(J),J-1) & & + DSPT(I+IHW(J),J+1) + DSPT(I+IHE(J),J+1) & & + 4. * DSPT(I,J)) * HBM2(I,J) * 0.125 ! WF4(I,J) = (CLDEFI(I+IHW(J),J-1) + CLDEFI(I+IHE(J),J-1) & & + CLDEFI(I+IHW(J),J+1) + CLDEFI(I+IHE(J),J+1) & & + 4. * CLDEFI(I,J)) * HBM2(I,J) * 0.125 345 END DO ! !------- ! OPENMP !------- !$omp parallel do private (RSMK,SMK) DO J=MYJS,MYJE DO I=MYIS,MYIE SMK = SM(I,J) RSMK = 1. - SMK DSPB(I,J) = DSPB(I,J) * SMK + WF1(I,J) * RSMK DSP0(I,J) = DSP0(I,J) * SMK + WF2(I,J) * RSMK DSPT(I,J) = DSPT(I,J) * SMK + WF3(I,J) * RSMK CLDEFI(I,J) = CLDEFI(I,J) * SMK + WF4(I,J) * RSMK END DO END DO ! 350 END DO ! END IF !----------------------------------------------- ! NITIALIZE CHANGES OF T AND Q DUE TO CONVECTION !----------------------------------------------- ! !------- ! OPENMP !------- !$omp parallel do DO 360 K=1,LM DO J=MYJS,MYJE DO I=MYIS,MYIE TMOD(I,J,K) = 0. QMOD(I,J,K) = 0. END DO END DO 360 END DO !------------------------------------------- ! CLEAN UP AND GATHER DEEP CONVECTION POINTS !------------------------------------------- ! !------- ! OPENMP !------- !$omp parallel do DO 380 J=MYJS2,MYJE2 DO 380 I=MYIS1,MYIE1 IF (LTOP(I,J) >= LBOT(I,J)) THEN LBOT(I,J) = 0 LTOP(I,J) = LBOT(I,J) PTOP(I,J) = PBOT(I,J) END IF IF (HBM2(I,J) < 0.90) THEN LBOT(I,J) = 0 LTOP(I,J) = LBOT(I,J) PTOP(I,J) = PBOT(I,J) END IF IF (PTOP(I,J) > PBOT(I,J)-PNO .OR. LTOP(I,J) > LBOT(I,J)-2) & & CLDEFI(I,J) = AVGEFI * SM(I,J) + STEFI * (1.-SM(I,J)) 380 END DO ! KHDEEP = 0 PSHNEW = 20000. DO J=MYJS2,MYJE2 DO I=MYIS1,MYIE1 PSFCIJ = PD(I,J) + PT !-------------------------------------------------------------------------------------------------- ! DEPTH OF CLOUD REQUIRED TO MAKE THE POINT A DEEP CONVECTION POINT IS NOW A SCALED VALUE OF THE ! PSFC INSTEAD OF 290 MB EVERYWHERE !-------------------------------------------------------------------------------------------------- DEPMIN = PSHNEW * PSFCIJ * 1.E-5 DEPTH = PBOT(I,J) - PTOP(I,J) IF (DEPTH >= DEPMIN) THEN KHDEEP = KHDEEP + 1 IDEEP(KHDEEP) = I JDEEP(KHDEEP) = J END IF END DO END DO ! !------------------------------------ ! HORIZONTAL LOOP FOR DEEP CONVECTION !------------------------------------ ! !------- ! OPENMP !------- !$omp parallel do !$omp private (APEK , APEKL , APEKXX, APEKXY, APESK , AVRGT , AVTGTL, DENTPY, DEPMIN, DEPTH , & !$omp DEPWL , DHDT , DIFQ , DIFQL , DIFT , DIFTL , DRHEAT, DRHDP , DSP , DSP0K , & !$omp DSPBK , DSPTK , DTHEM , EFI , FEFI , HCORR , I , J , L0 , L0M1 , & !$omp LB , LBM1 , LBTK , LCOR , LQM , LSHU , LTP1 , LTP2 , LTPK , LTSH , & !$omp PBTK , PK , PK0 , PKB , PKL , PKT , PRECK , PSFCIJ, PSK , PTHRS , & !$omp PTPK , QK , QKL , QREFK , QSATK , RDP0T , RHH , RHL , RHMAX , SUMDE , & !$omp SUMDP , THERK , THERKX, THERKY, THSK , THSKL , TK , TKL , TREFK , TREFKX, & !$omp TSKL) ! DO 600 N=1,KHDEEP ! I = IDEEP(N) J = JDEEP(N) PSFCIJ = PD(I,J) + PT LTPK = LTOP(I,J) LBTK = LBOT(I,J) !--------------- !DEEP CONVECTION !--------------- LB = LBTK EFI = CLDEFI(I,J) DSPBK = DSPB(I,J) DSP0K = DSP0(I,J) DSPTK = DSPT(I,J) !-------------------------------------------------------------------------------------------------- ! INITIALIZE VARIABLES IN THE CONVECTIVE COLUMN ! ! ONE SHOULD NOTE THAT THE VALUES ASSIGNED TO THE ARRAY TREFK IN THE 410 LOOP ARE REALLY ONLY ! RELEVANT IN ANCHORING THE REFERENCE TEMPERATURE PROFILE AT LEVEL LB. ! WHEN BUILDING THE REFERENCE PROFILE FROM CLOUD BASE, THEN ASSIGNING THE AMBIENT TEMPERATURE TO ! TREFK IS ACCEPTABLE. ! HOWEVER, WHEN BUILDING THE REFERENCE PROFILE FROM SOME OTHER LEVEL (SUCH AS ONE LEVEL ABOVE THE ! GROUND), THEN TREFK SHOULD BE FILLED WITH THE TEMPERATURES IN TREF(I,J,L) WHICH ARE THE ! TEMPERATURES OF THE MOIST ADIABAT THROUGH CLOUD BASE. ! BY THE TIME THE LINE NUMBERED 450 HAS BEEN REACHED, TREFK ACTUALLY DOES HOLD THE REFERENCE ! TEMPERATURE PROFILE. !-------------------------------------------------------------------------------------------------- DO 410 K=1,LM DIFT(K) = 0. DIFQ(K) = 0. TKL = T(I,J,K) TK(K) = TKL TREFK(K) = TKL QKL = Q(I,J,K) QK(K) = QKL QREFK(K) = QKL PKL = AETA(K) * PDSL(I,J) + PT PK(K) = PKL PSK(K) = PKL APEKL = APE(I,J,K) APEK(K) = APEKL THERK(K) = TREF(I,J,K) * APEKL 410 END DO !---------------------------------------------- ! DEEP CONVECTION REFERENCE TEMPERATURE PROFILE !---------------------------------------------- LTP1 = LTPK + 1 LBM1 = LB - 1 PKB = PK(LB) PKT = PK(LTPK) !--------------------------------------------------- ! TEMPERATURE REFERENCE PROFILE BELOW FREEZING LEVEL !--------------------------------------------------- L0 = LB ! PK0 = PK(LB) TREFKX = TREFK(LB) THERKX = THERK(LB) APEKXX = APEK(LB) THERKY = THERK(LBM1) APEKXY = APEK(LBM1) ! DO 420 K=LBM1,LTPK,-1 IF (T(I,J,K+1) < TFRZ) GO TO 430 STABDL = STABD TREFKX = ((THERKY-THERKX) * STABDL + TREFKX * APEKXX) / APEKXY TREFK(K) = TREFKX APEKXX = APEKXY THERKX = THERKY APEKXY = APEK(K-1) THERKY = THERK(K-1) L0 = K PK0 = PK(L0) 420 END DO !----------------------------------------- ! FREEZING LEVEL AT OR ABOVE THE CLOUD TOP !----------------------------------------- L0M1 = L0 - 1 GO TO 450 !--------------------------------------------------- ! TEMPERATURE REFERENCE PROFILE ABOVE FREEZING LEVEL !--------------------------------------------------- 430 L0M1 = L0 - 1 RDP0T = 1. / (PK0 - PKT) DTHEM = THERK(L0) - TREFK(L0) * APEK(L0) ! DO K=LTPK,L0M1 TREFK(K) = (THERK(K) - (PK(K) - PKT) * DTHEM * RDP0T) / APEK(K) END DO !------------------------------------------- ! DEEP CONVECTION REFERENCE HUMIDITY PROFILE !------------------------------------------- ! !-------------------------------------------------------------------------------------------------------------------------- ! REFERENCE PROFILE HAD BEEN TOO DRY IN THE CASE WHERE THE CLOUD BASE WAS CLOSE TO THE FREEZING LEVEL - THIS IS NOW CHANGED !-------------------------------------------------------------------------------------------------------------------------- 450 DEPTH = PFRZ * PSFCIJ * 1.E-5 DEPWL = PKB - PK0 DO 460 K=LTPK,LB !------------------------------- ! SATURATION PRESSURE DIFFERENCE !------------------------------- IF (DEPWL >= DEPTH) THEN IF (K < L0) THEN DSP = ((PK0 - PK(K)) * DSPTK + (PK(K) - PKT) * DSP0K) / (PK0 - PKT) ELSE DSP = ((PKB - PK(K)) * DSP0K + (PK(K) - PK0) * DSPBK) / (PKB - PK0) END IF ELSE DSP = DSP0K IF (K < L0) DSP = ((PK0 - PK(K)) * DSPTK + (PK(K) - PKT) * DSP0K) / (PK0 - PKT) END IF !----------------- ! HUMIDITY PROFILE !----------------- IF (PK(K) > PQM) THEN PSK(K) = PK(K) + DSP APESK(K) = (1.E5 / PSK(K)) ** CAPA THSK(K) = TREFK(K) * APEK(K) QREFK(K) = PQ0 / PSK(K) * EXP(A2 * (THSK(K) - A3 * APESK(K)) & & / (THSK(K) - A4 * APESK(K))) ELSE QREFK(K) = Q(I,J,K) END IF 460 END DO !------------------------------- ! ENTHALPY CONSERVATION INTEGRAL !------------------------------- DO 520 ITER=1,2 ! SUMDE = 0. SUMDP = 0. ! DO K=LTPK,LB SUMDE = ((TK(K) - TREFK(K)) * CP + (QK(K) - QREFK(K)) * ELWV) * DETA(K) + SUMDE SUMDP = SUMDP + DETA(K) END DO ! HCORR = SUMDE / (SUMDP - DETA(LTPK)) LCOR = LTPK + 1 !--------- ! FIND LQM !--------- DO K=1,LB IF (PK(K) <= PQM) LQM = K END DO !----------------------------------- ! ABOVE LQM CORRECT TEMPERATURE ONLY !----------------------------------- IF (LCOR <= LQM) THEN DO K=LCOR,LQM TREFK(K) = TREFK(K) + HCORR * RCP END DO LCOR = LQM + 1 END IF !------------------------------------------------ ! BELOW LQM CORRECT BOTH TEMPERATURE AND MOISTURE !------------------------------------------------ DO 510 K=LCOR,LB TSKL = TREFK(K) * APEK(K) / APESK(K) DHDT = QREFK(K) * A23M4L / (TSKL-A4) ** 2 + CP TREFK(K) = HCORR / DHDT + TREFK(K) THSKL = TREFK(K) * APEK(K) QREFK(K) = PQ0 / PSK(K) * EXP(A2 * (THSKL - A3 * APESK(K)) & & / (THSKL - A4 * APESK(K))) 510 END DO ! 520 END DO !----------------------------------- ! HEATING, MOISTENING, PRECIPITATION !----------------------------------- DENTPY = 0. AVRGT = 0. PRECK = 0. ! DO 530 K=LTPK,LB TKL = TK(K) DIFTL = (TREFK(K) - TKL ) * TAUK DIFQL = (QREFK(K) - QK(K)) * TAUK AVRGTL = (TKL + TKL + DIFTL) DENTPY = (DIFTL * CP + DIFQL * ELWV) * DETA(K) / AVRGTL + DENTPY AVRGT = AVRGTL * DETA(K) + AVRGT PRECK = DETA(K) * DIFTL + PRECK DIFT(K) = DIFTL DIFQ(K) = DIFQL 530 END DO ! DENTPY = DENTPY + DENTPY AVRGT = AVRGT / (SUMDP + SUMDP) !------------------------------------- ! SWAP IF ENTROPY AND/OR PRECIP .LT. 0 !------------------------------------- IF (DENTPY < EPSNTP .OR. PRECK < 0.) THEN IF (OCT90) THEN CLDEFI(I,J) = EFIMN ELSE CLDEFI(I,J) = EFIMN * SM(I,J) + STEFI * (1. - SM(I,J)) END IF !----------------------------- ! SEARCH FOR SHALLOW CLOUD TOP !----------------------------- LBTK = LBOT(I,J) LTSH = LBTK LBM1 = LBTK - 1 PBTK = PK(LBTK) !---------------------------------- ! USE NEW THRESHOLD FOR CLOUD DEPTH !---------------------------------- PSFCIJ = PD(I,J) + PT DEPMIN = PSHNEW * PSFCIJ * 1.E-5 PTPK = PBTK - DEPMIN !------------------------------------------- ! CLOUD TOP IS THE LEVEL JUST BELOW PBTK-PSH !------------------------------------------- DO K=1,LM IF (PK(K) <= PTPK) LTPK = K + 1 END DO PTPK=PK(LTPK) !----------------------------------------------- ! HIGHEST LEVEL ALLOWED IS LEVEL JUST BELOW PSHU !----------------------------------------------- IF (PTPK <= PSHU) THEN DO K=1,LM IF (PK(K) <= PSHU) LSHU = K + 1 END DO LTPK = LSHU PTPK = PK(LTPK) END IF ! LTP1 = LTPK + 1 LTP2 = LTPK + 2 ! DO K=LTPK,LBTK QSATK(K) = PQ0 / PK(K) * EXP(A2 * (TK(K) - A3) / (TK(K) - A4)) END DO ! RHH = QK(LTPK) / QSATK(LTPK) RHMAX = 0. ! DO 570 K=LTP1,LBM1 RHL = QK(K) / QSATK(K) DRHDP = (RHH-RHL) / (PK(K-1) - PK(K)) IF (DRHDP > RHMAX) THEN LTSH = K - 1 RHMAX = DRHDP END IF RHH = RHL 570 END DO ! LTOP(I,J) = LTSH !---------------------------------------- ! CLOUD MUST BE AT LEAST TWO LAYERS THICK !---------------------------------------- IF (LBOT(I,J)-LTOP(I,J) < 2) LTOP(I,J) = LBOT(I,J) - 2 ! PTOP(I,J) = PK(LTOP(I,J)) GO TO 600 END IF !-------------------------- ! DEEP CONVECTION OTHERWISE !-------------------------- DRHEAT = (PRECK * SM(I,J) + AMAX1(EPSP,PRECK) * (1. - SM(I,J))) * CP / AVRGT EFI = EFIFC * DENTPY / DRHEAT !----------------------------------- ! UNIFIED OR SEPARATE LAND/SEA CONV. !----------------------------------- IF (OCT90) THEN IF (UNIS) THEN EFI = CLDEFI(I,J) * FCB + EFI * FCC ELSE IF ( .NOT. UNIL) THEN EFI = (CLDEFI(I,J) * FCB + EFI * FCC) * SM(I,J) + 1. - SM(I,J) ELSE EFI = 1. END IF ELSE EFI = CLDEFI(I,J) * FCB + EFI * FCC END IF ! IF(EFI > 1. ) EFI = 1. IF(EFI < EFIMN) EFI = EFIMN IF(PRECK == 0.) EFI = 1. ! CLDEFI(I,J) = EFI ! FEFI = EFMNT + SLOPE * (EFI - EFIMN) FEFI = 2. - FEFI PRECK = PRECK * FEFI !----------------------------------------------- ! UPDATE PRECIPITATION, TEMPERATURE AND MOISTURE !----------------------------------------------- PRECOL(I,J) = PDSL(I,J) * PRECK * CPRLG PREC (I,J) = PDSL(I,J) * PRECK * CPRLG + PREC (I,J) CUPREC(I,J) = PDSL(I,J) * PRECK * CPRLG + CUPREC(I,J) ACPREC(I,J) = PDSL(I,J) * PRECK * CPRLG + ACPREC(I,J) CUPPT(I,J) = PDSL(I,J) * PRECK * CPRLG + CUPPT (I,J) ! DO K=LTPK,LB TMOD(I,J,K) = DIFT(K) * FEFI QMOD(I,J,K) = DIFQ(K) * FEFI TLATCU(I,J,K) = DIFT(K) * FEFI END DO ! !------------------------------------------------ ! PCC SCHEME (Acordingly to Koh and Fonseca 2016) !------------------------------------------------ ! PDQD = 0. PRWD = 0. DQCOL = 0. DQCOLMIN = 0. PWCOL = 0. CONVCLD = 0. ! DO K=1,LM T_COL(K) = 0. QV_COL(K) = 0. RHO(I,J,K) = 0. DZ8W(I,J,K) = 0. PMID(I,J,K) = 0. JPR(K) = 0. PVPR(K) = 0. QCCONV(I,J,K) = 0. QICONV(I,J,K) = 0. CCLDFRA(I,J,K) = 0. END DO ! DO K=1,LM ! !*** CONVERT FROM MIXING RATIO TO SPECIFIC HUMIDITY ! PMID(I,J,K) = PDSL(I,J) * AETA(K) + PT !HEIGHT OF MIDLAYER ETA LEVELS DZ8W(I,J,K) = PDSL(I,J) * DETA(K) + PT !DZ BETWEEN FULL LEVELS T_COL(K) = T(I,J,K) !TEMPERATURE QV_COL(K) = MAX(EPSQ, Q(I,J,K)/(1. + Q(I,J,K))) !MAXIMUM COLUMN HUMIDITY RHO(I,J,K) = DZ8W(I,J,K) / (RD * T_COL(K) * (1. + EPS1 * QV_COL(K))) !AIR DENSITY (KG/M**3) ! PRWD = QK(K) * DETA(K) + PRWD !ETA LEVELS INTERFACE SATURATION ?? DPCOL(K) = RHO(I,J,K) * G * DZ8W(I,J,K) !HYDROSTATIC EQUILIBRIUM EQUATION PDQD = (QK(K) - QREFK(K)) * DPCOL(K) + PDQD !SATURATION ADJUSTMENT ! PWCOL = PRWD / G !HUMIDITY BUOYANCY ?? DQCOL = PDQD / G !LOSS OF WATER VAPOUR MASS FROM THE !CLOUD COLUMN BY PRECIPITATION DQCOLMIN = (CRMN * TREL) / (FEFI * 86400.) !DQCOL MINIMUM VALUE ! END DO ! !------------------------------------ ! CONDITION TO MAKE PCC SCHEME WORKS! !------------------------------------ ! IF (DQCOL > DQCOLMIN) THEN ! !----------------------------------------------------------------- ! CONVECTIVE CLOUD FRACTION: BASED ON SLINGO (1987) WITH A POISSON ! VERTICAL PROFILE. PLEASE NOTE THAT THE BMJ PRECIPITATION RATE ! (PRECOL) HAS TO BE CONVERTED FROM MMS-1 TO MMDAY-1. !----------------------------------------------------------------- TTOP = 0. BBOT = 0. ! PAVG = 0. CUMX = 0. PPAVG = 0. MPVPR = 0. FACTL = 0. ! !---------------------------------------------------------------------------- !USING PRECOL(I,J)*1000 to convert to mm and 86400.0/DTCNVC TO CONVERT TO DAY !---------------------------------------------------------------------------- ! PRRT = ((PRECOL(I,J) * 1000.0)*(86400.0 / DTCNVC)) / CRMN RRP = 0.8/(LOG(CRMX/CRMN)) ! IF (PRRT < (CRMX/CRMN)) THEN CUMX = RRP*LOG(PRRT) ELSE CUMX = 0.8 END IF ! !------------------------------------------------------------------------- ! COMPUTE THE CONVECTIVE CLOUD FRACTION (CCLDFRA) AT EACH MODEL LEVEL. ! FOR N>=17 USE THE STERLING APPROXIMATION AS FOR N=17 IT GIVES A RELATIVE ! ERROR OF ~9.4x10E-8. !------------------------------------------------------------------------- ! TTOP = LTOP(I,J) BBOT = LBOT(I,J) PAVG = 1./(PEPS*3.)**2 ! DO K=1,LM IF (K.GE.TTOP.AND.K.LE.BBOT) THEN JPR(K)=(1./PEPS)*((PMID(I,J,K)-PMID(I,J,TTOP))/(PMID(I,J,BBOT)-PMID(I,J,TTOP))) ELSE JPR(K)=0.0 END IF END DO ! DO K=1,LM PVPR(K)=0. END DO ! ! "IF" WAS REMOVED FOR TESTS !diego IF (JPR(TTOP).LT.17) THEN DO K=TTOP,BBOT PVPR(K) = (PAVG)**(JPR(K))/GAMMA(JPR(K)+1.) END DO !diego ELSE !diego DO K=TTOP,BBOT !diego FACTL = JPR(K)*LOG(JPR(K))-JPR(K)+1./2.*LOG(2.*JPR(K)*ACOS(-1.))+ & !diego & LOG(1.+1./(12.*JPR(K))+1./(288.*JPR(K)**2)) !diego PVPR(K) = EXP((JPR(K))*LOG(1.*PAVG)-FACTL) !diego END DO !diego END IF ! MPVPR = MAXVAL(PVPR) DO K=TTOP,BBOT PVPR(K) = PVPR(K) / MPVPR END DO ! DO K=1,LM CCLDFRA(I,J,K) = CUMX*PVPR(K) END DO ! !--------------------------------------------------------------------------- ! COMPUTE THE CONVECTIVE CLOUD CONDENSATES (QCCONV,QICONV). PLEASE NOTE THAT ! THE EQUATION FOR QCIS IS VALID FOR PRRTs IN THE RANGE 10**(-7) TO 10**3. !--------------------------------------------------------------------------- ! QCIS=0. MCOL=0. ! DO K=LTOP(I,J),LBOT(I,J) MCOL = RHO(I,J,K)*DZ8W(I,J,K)*QV_COL(K)*CCLDFRA(I,J,K) + MCOL END DO ! CONVCLD(I,J) = PWCOL**GAM*(DQCOL-DQCOLMIN)**(1.-GAM) QCIS = CONVCLD(I,J) / MCOL ! !cloud bottom and top pressure level ! IF(MYPE.EQ.732.AND.I.EQ.1.AND.J.EQ.1)THEN !diego OPEN(1006,FILE="CUCNVC_TOPBOT.bin", FORM="UNFORMATTED", POSITION="APPEND", ACCESS="SEQUENTIAL") OPEN(1006,FILE="CUCNVC_TOPBOT.txt", FORM="FORMATTED", POSITION="APPEND") WRITE(1006,*) "NTSD=", NTSD WRITE(1006,*) "LBOT(I,J)=", LBOT(I,J) WRITE(1006,*) "LTOP(I,J)=", LTOP(I,J) !diego CLOSE(1006) END IF ! !---------------------------------------------------------------------------- ! QCCONV AND QICONV MULTIPLIED BY AN ADJUSTMENT FACTOR (SUGESTION DIEGO) ! IN ORDER TO BE COMPARABLE TO CONDENSATES FROM MYCROPHYSICS (QPI3D and QPR3D) !---------------------------------------------------------------------------- ! DO K=1,LM !diego KFLIP = LM + 1 - K IF (T_COL(K) >= TFRZ) THEN QICONV(I,J,K) = 0. QCCONV(I,J,K) = 100*((QCIS*QV_COL(K)*CCLDFRA(I,J,K))/(1.-QV_COL(K))) ELSE QICONV(I,J,K) = 100*((QCIS*QV_COL(K)*CCLDFRA(I,J,K))/(1.-QV_COL(K))) QCCONV(I,J,K) = 0. END IF END DO ! ELSE ! IF(MYPE.EQ.732.AND.I.EQ.1.AND.J.EQ.1)THEN !diego OPEN(1006,FILE="CUCNVC_TOPBOT.txt", FORM="FORMATTED", POSITION="APPEND") WRITE(1006,*) "NTSD=", NTSD WRITE(1006,*) "LBOT(I,J)=", 0. WRITE(1006,*) "LTOP(I,J)=", 0. CLOSE(1006) END IF ! !------------------------------------------------------------------------------- ! END IF !(DQCOL - DQCOLMIN) ! !---------------diego----------------------------------------------------------- !------------------ ! END OF PCC SCHEME !------------------ ! !----------------------- ! END OF DEEP CONVECTION !----------------------- ! 600 END DO ! NDEEP = 0 ! DO 620 J=MYJS2,MYJE2 DO 620 I=MYIS,MYIE LTPK = LTOP(I,J) LBTK = LBOT(I,J) LB = LMH(I,J) - 1 PSFCIJ = PD(I,J) + PT DEPMIN = PSHNEW * PSFCIJ * 1.E-5 ! IF (PTOP(I,J) < PBOT(I,J)-DEPMIN) THEN NDEEP = NDEEP + 1 NDEPTH = LB - LTPK NTOPD (LTPK ) = NTOPD (LTPK ) + 1 NBOTD (LB ) = NBOTD (LB ) + 1 NDPTHD(NDEPTH) = NDPTHD(NDEPTH) + 1 END IF 620 END DO NNEG = KHDEEP - NDEEP !--------------------------------- ! GATHER SHALLOW CONVECTION POINTS !--------------------------------- KHSHAL = 0 NDSTN = 0 NDSTP = 0 ! DO 630 J=MYJS2,MYJE2 DO 630 I=MYIS,MYIE IF (PTOP(I,J) > PBOT(I,J)-PNO .OR. LTOP(I,J) > LBOT(I,J)-2) GO TO 630 PSFCIJ = PD(I,J) + PT DEPMIN = PSHNEW * PSFCIJ * 1.E-5 IF (PTOP(I,J)+1. >= PBOT(I,J)-DEPMIN) THEN KHSHAL = KHSHAL + 1 ISHAL(KHSHAL) = I JSHAL(KHSHAL) = J END IF 630 END DO !---------------------------------------- ! HORIZONTAL LOOP FOR SHALLOW CONVECTION !---------------------------------------- ! !------- ! OPENMP !------- !$omp parallel do !$omp private(APEK, APEKL, APEKXX, APEKXY, BQK, BQS00K, BQS10K, DEN, DENTPY, DPKL, DPMIX, DQREF, & !$omp DST, DSTQ, DTDETA, FPK, FPTK, I, IQ, IT, J, LBM1, LBTK, LTP1, LTPK, OTSUM, PART1, & !$omp PART2, PART3, PK, PKL, PKXXXX, PKXXXY, POTSUM, PPK, PSUM, PTPK, PZ0, QK, QKL, QNEW, & !$omp QOTSUM, QQK, QREFK, QRFKL, QRFTP, QSATK, QSUM, RDPSUM, RTBAR, SMIX, SQK, SQS00K, & !$omp SQS10K, SUMDP, SUMDT, TCORR, THVMKL, THVREF, TK, TKL, TQK, TREFK, TREFKX, TRFKL, & !$omp TTHK) ! DO 800 N=1,KHSHAL ! I = ISHAL(N) J = JSHAL(N) !------------------- ! SHALLOW CONVECTION !------------------- PZ0 = PD(I,J) + PT LLMH = LMH(I,J) ! DO 650 K=1,LLMH TKL = T(I,J,K) TK (K) = TKL TREFK(K) = TKL QKL = Q(I,J,K) QK (K) = QKL QREFK(K) = QKL PKL = AETA(K) * PDSL(I,J) + PT PK (K) = PKL QSATK(K) = PQ0/PK(K) * EXP(A2 * (TK(K) - A3) / (TK(K) - A4)) APEKL = APE(I,J,K) !----------------------------- ! CHOOSE THE PRESSURE FUNCTION ! ! FPK (K) =ALOG(PKL) ! FPK (K) =PKL ! FPK (K) =-1./PKL !----------------------------- APEK (K) = APEKL THVMKL = TKL * APEKL * (QKL * 0.608 + 1.) THVREF(K) = THVMKL 650 END DO !------------------ ! SHALLOW CLOUD TOP !------------------ LBTK = LBOT(I,J) LBM1 = LBTK - 1 PTPK = PTOP(I,J) LTP1 = LTOP(I,J) LTPK = LTOP(I,J) - 1 ! IF (PTOP(I,J) > PBOT(I,J)-PNO .OR. LTOP(I,J) > LBOT(I,J)-2) THEN LBOT(I,J) = 0 LTOP(I,J) = LBOT(I,J) PTOP(I,J) = PBOT(I,J) GO TO 800 END IF !----------------------------------------------------- ! SCALING POTENTIAL TEMPERATURE AND TABLE INDEX AT TOP !----------------------------------------------------- THTPK = T(I,J,LTPK) * APE(I,J,LTPK) ! TTHK = (THTPK - THL) * RDTH QQK = TTHK - AINT(TTHK) IT = INT(TTHK) + 1 IF (IT < 1) THEN IT = 1 QQK = 0. END IF IF (IT >= JTB) THEN IT = JTB - 1 QQK = 0. END IF !-------------------------------------------------- ! BASE AND SCALING FACTOR FOR SPEC. HUMIDITY AT TOP !-------------------------------------------------- BQS00K = QS0(IT) SQS00K = SQS(IT) BQS10K = QS0(IT+1) SQS10K = SQS(IT+1) !---------------------------------------------- ! SCALING SPEC. HUMIDITY AND TABLE INDEX AT TOP !---------------------------------------------- BQK = (BQS10K - BQS00K) * QQK + BQS00K SQK = (SQS10K - SQS00K) * QQK + SQS00K ! TQK = (Q(I,J,LTPK) - BQK) / SQK * RDQ ! PPK = TQK - AINT(TQK) IQ = INT(TQK) + 1 IF (IQ < 1) THEN IQ = 1 PPK = 0. END IF IF (IQ >= ITB) THEN IQ = ITB - 1 PPK = 0. END IF !------------------------------------ ! CLOUD TOP SATURATION POINT PRESSURE !------------------------------------ PART1 = (PTBL(IQ+1,IT) - PTBL(IQ ,IT)) * PPK PART2 = (PTBL(IQ,IT+1) - PTBL(IQ ,IT)) * QQK PART3 = (PTBL(IQ,IT ) - PTBL(IQ+1,IT) - PTBL(IQ,IT+1) + PTBL(IQ+1,IT+1)) * PPK * QQK PTPK = PTBL(IQ,IT ) + PART1 + PART2 + PART3 ! DPMIX = PTPK - PSP(I,J) IF (ABS(DPMIX) < 3000.) DPMIX = -3000. !-------------------------- ! TEMPERATURE PROFILE SLOPE !-------------------------- SMIX = (THTPK - THBT(I,J)) / DPMIX * STABS ! TREFKX = TREFK(LBTK+1) PKXXXX = PK(LBTK+1) PKXXXY = PK(LBTK) APEKXX = APEK(LBTK+1) APEKXY = APEK(LBTK) ! DO 670 K=LBTK,LTP1,-1 TREFKX = ((PKXXXY - PKXXXX) * SMIX + TREFKX * APEKXX) / APEKXY TREFK(K) = TREFKX APEKXX = APEKXY PKXXXX = PKXXXY APEKXY = APEK(K-1) PKXXXY = PK(K-1) 670 END DO !----------------------------------------- ! TEMPERATURE REFERENCE PROFILE CORRECTION !----------------------------------------- SUMDT = 0. SUMDP = 0. ! DO K=LTP1,LBTK SUMDT = (TK(K) - TREFK(K)) * DETA(K) + SUMDT SUMDP = SUMDP + DETA(K) END DO ! RDPSUM = 1. / SUMDP FPK(LBTK) = TREFK(LBTK) ! TCORR=SUMDT*RDPSUM ! DO K=LTP1,LBTK TRFKL = TREFK(K)+TCORR TREFK(K) = TRFKL FPK (K) = TRFKL END DO !--------------------------- ! HUMIDITY PROFILE EQUATIONS !--------------------------- PSUM = 0. QSUM = 0. POTSUM = 0. QOTSUM = 0. OTSUM = 0. DST = 0. FPTK = FPK(LTP1) ! DO 700 K=LTP1,LBTK DPKL = FPK(K) - FPTK PSUM = DPKL * DETA(K) + PSUM QSUM = QK(K) * DETA(K) + QSUM RTBAR = 2. / (TREFK(K) + TK(K)) OTSUM = DETA(K) * RTBAR + OTSUM POTSUM = DPKL * RTBAR * DETA(K) + POTSUM QOTSUM = QK(K) * RTBAR * DETA(K) + QOTSUM DST = (TREFK(K) - TK(K)) * RTBAR * DETA(K) + DST 700 END DO ! PSUM = PSUM * RDPSUM QSUM = QSUM * RDPSUM ROTSUM = 1. / OTSUM POTSUM = POTSUM * ROTSUM QOTSUM = QOTSUM * ROTSUM DST = DST * ROTSUM * CP / ELWV !------------------------------- ! ENSURE POSITIVE ENTROPY CHANGE !------------------------------- IF (DST > 0.) THEN LBOT(I,J) = 0 LTOP(I,J) = LBOT(I,J) PTOP(I,J) = PBOT(I,J) GO TO 800 ELSE DSTQ = DST * EPSDN END IF !-------------------------------- ! CHECK FOR ISOTHERMAL ATMOSPHERE !-------------------------------- DEN = POTSUM - PSUM ! IF (-DEN/PSUM < 5.E-5) THEN LBOT(I,J) = 0 LTOP(I,J) = LBOT(I,J) PTOP(I,J) = PBOT(I,J) GO TO 800 !---------------------------------------- ! SLOPE OF THE REFERENCE HUMIDITY PROFILE !---------------------------------------- ELSE DQREF = (QOTSUM-DSTQ-QSUM) / DEN END IF !--------------------------------------- ! HUMIDITY DOES NOT INCREASE WITH HEIGHT !--------------------------------------- IF (DQREF < 0.) THEN LBOT(I,J) = 0 LTOP(I,J) = LBOT(I,J) PTOP(I,J) = PBOT(I,J) GO TO 800 END IF !-------------------------- ! HUMIDITY AT THE CLOUD TOP !-------------------------- QRFTP = QSUM - DQREF * PSUM !----------------- ! HUMIDITY PROFILE !----------------- DO 720 K=LTP1,LBTK QRFKL = (FPK(K) - FPTK) * DQREF + QRFTP !------------------------------------------ ! SUPERSATURATION OR NEGATIVE Q NOT ALLOWED !------------------------------------------- QNEW = (QRFKL - QK(K)) * TAUK + QK(K) ! IF (QNEW > QSATK(K)*STRESH .OR. QNEW < 0.) THEN LBOT(I,J) = 0 LTOP(I,J) = LBOT(I,J) PTOP(I,J) = PBOT(I,J) GO TO 800 END IF ! THVREF(K) = TREFK(K) * APEK(K) * (QRFKL * 0.608 + 1.) QREFK(K) = QRFKL 720 END DO !----------------------------------------------- ! ELIMINATE IMPOSSIBLE SLOPES (BETTS, DTHETA/DQ) !----------------------------------------------- DO 730 K=LTP1,LBTK DTDETA = (THVREF(K-1) - THVREF(K)) / (AETA(K) - AETA(K-1)) IF (DTDETA < EPSTH) THEN LBOT(I,J) = 0 LTOP(I,J) = LBOT(I,J) PTOP(I,J) = PBOT(I,J) GO TO 800 END IF 730 END DO !------------------ ! DIAGNOSTICS ! ! IF(DST.GT.0.)THEN ! NDSTP=NDSTP+1 ! ELSE ! NDSTN=NDSTN+1 ! ENDIF !------------------ DENTPY = 0. ! DO K=LTP1,LBTK DENTPY = ((TREFK(K) - TK(K)) * CP + (QREFK(K) - QK(K)) * ELWV) & & / (TK(K) + TREFK(K)) * DETA(K) + DENTPY END DO !-------------------------------------- ! RELAXATION TOWARDS REFERENCE PROFILES !-------------------------------------- DO 750 K=LTP1,LBTK TMOD(I,J,K) = (TREFK(K) - TK(K)) * TAUK QMOD(I,J,K) = (QREFK(K) - QK(K)) * TAUK 750 END DO !-------------------------- ! END OF SHALLOW CONVECTION !-------------------------- 800 END DO !------------ ! DIAGNOSTICS !------------ NSHAL = 0 ! DO 820 J=MYJS2,MYJE2 DO 820 I=MYIS,MYIE LTPK = LTOP(I,J) LBTK = LBOT(I,J) PTPK = PTOP(I,J) PBTK = PBOT(I,J) IF (PTPK > PBTK-PNO .OR. LTPK > LBTK-2) GO TO 820 PSFCIJ = PD(I,J) + PT DEPMIN = PSHNEW * PSFCIJ * 1.E-5 ! IF (PTPK >= PBTK-DEPMIN) THEN NSHAL = NSHAL + 1 NTOPS(LTPK) = NTOPS(LTPK) + 1 NBOTS(LBTK) = NBOTS(LBTK) + 1 NDEPTH = LBTK - LTPK NDPTHS(NDEPTH) = NDPTHS(NDEPTH) + 1 END IF 820 END DO NEGDS = KHSHAL - NSHAL !--------------------------------------------- ! SMOOTHING TEMPERATURE & HUMIDITY CORRECTIONS !--------------------------------------------- IF (KSMUD == 0) THEN ! !------- ! OPENMP !------- !$omp parallel do ! DO 900 K=1,LM !----------------------------------------- ! UPDATE THE FUNDAMENTAL PROGNOSTIC ARRAYS !----------------------------------------- DO 830 J=MYJS,MYJE DO 830 I=MYIS,MYIE T(I,J,K) = T(I,J,K) + TMOD(I,J,K) Q(I,J,K) = Q(I,J,K) + QMOD(I,J,K) !--------------------------------------------------------------------------------------------------------------- ! ACCUMULATE LATENT HEATING DUE TO CONVECTION. ! SCALE BY THE RECIPROCAL OF THE PERIOD AT WHICH THIS ROUTINE IS CALLED. THIS PERIOD IS THE CONVECTION TIMESTEP. !--------------------------------------------------------------------------------------------------------------- TCUCN(I,J,K) = TCUCN(I,J,K) + TMOD(I,J,K) * RDTCNVC 830 END DO 900 END DO ELSE ! !------- ! OPENMP !------- !$omp parallel do private (QL,QNE,QSE,TL,TNE,TSE) ! DO 910 K=1,LM ! CALL ZERO2(QL) CALL ZERO2(QNE) CALL ZERO2(QSE) CALL ZERO2(TL) CALL ZERO2(TNE) CALL ZERO2(TSE) ! DO J=MYJS,MYJE DO I=MYIS,MYIE TL(I,J) = TMOD(I,J,K) QL(I,J) = QMOD(I,J,K) END DO END DO ! NSMUD = KSMUD ! DO 870 KS=1,NSMUD ! DO J=MYJS,MYJE1 DO I=MYIS,MYIE TNE(I,J) = (TL(I+IHE(J),J+1) - TL(I,J)) * HTM(I,J,K) * HTM(I+IHE(J),J+1,K) QNE(I,J) = (QL(I+IHE(J),J+1) - QL(I,J)) * HTM(I,J,K) * HTM(I+IHE(J),J+1,K) END DO END DO ! DO J=MYJS1,MYJE DO I=MYIS,MYIE TSE(I,J) = (TL(I+IHE(J),J-1) - TL(I,J)) * HTM(I+IHE(J),J-1,K) * HTM(I,J,K) QSE(I,J) = (QL(I+IHE(J),J-1) - QL(I,J)) * HTM(I+IHE(J),J-1,K) * HTM(I,J,K) END DO END DO ! DO J=MYJS2,MYJE2 DO I=MYIS,MYIE TL(I,J) = (TNE(I,J) - TNE(I+IHW(J),J-1) + TSE(I,J) - TSE(I+IHW(J),J+1)) & & * HBM2(I,J) * 0.125 + TL(I,J) QL(I,J) = (QNE(I,J) - QNE(I+IHW(J),J-1) + QSE(I,J) - QSE(I+IHW(J),J+1)) & & * HBM2(I,J) * 0.125 + QL(I,J) END DO END DO ! 870 END DO !----------------------------------------- ! UPDATE THE FUNDAMENTAL PROGNOSTIC ARRAYS !----------------------------------------- DO J=MYJS,MYJE DO I=MYIS,MYIE T(I,J,K) = T(I,J,K) + TL(I,J) Q(I,J,K) = Q(I,J,K) + QL(I,J) END DO END DO !--------------------------------------------------------------------------------------------------------------- ! ACCUMULATE LATENT HEATING DUE TO CONVECTION. ! SCALE BY THE RECIPROCAL OF THE PERIOD AT WHICH THIS ROUTINE IS CALLED. THIS PERIOD IS THE CONVECTION TIMESTEP. !--------------------------------------------------------------------------------------------------------------- DO J=MYJS,MYJE DO I=MYIS,MYIE TCUCN(I,J,K) = TCUCN(I,J,K) + TL(I,J) * RDTCNVC END DO END DO ! 910 END DO ! END IF ! !---------------------------------------- ! SAVE CLOUD TOP AND BOTTOM FOR RADIATION !---------------------------------------- ! !------- ! OPENMP !------- !$omp parallel do ! DO J=MYJS,MYJE DO I=MYIS,MYIE IF (LTOP(I,J) > 0 .AND. LTOP(I,J) < LBOT(I,J)) THEN CUTOP = FLOAT(LTOP(I,J)) CUBOT = FLOAT(LBOT(I,J)) HTOP(I,J) = MIN(CUTOP, HTOP(I,J)) HBOT(I,J) = MAX(CUBOT, HBOT(I,J)) CNVTOP(I,J) = MIN(CUTOP,CNVTOP(I,J)) CNVBOT(I,J) = MAX(CUBOT,CNVBOT(I,J)) END IF END DO END DO ! !--------------------------------------------------------------------------------------------------------------------- ! DIAGNOSTICS ! ! WRITE(LIST,950)NTSD,NSHAL,NDEEP,NNEG,NEGDS,NDSTN,NDSTP ! ! DO 940 L=1,LM ! WRITE(LIST,952)L ! WRITE(LIST,954)NBOTS(L),NTOPS(L),NDPTHS(L),NBOTD(L),NTOPD(L),NDPTHD(L) ! 940 CONTINUE ! 950 FORMAT(' NTSD=',I3,I8,' SHALLOW, ',I4,' DEEP, ',I4,' NEG., ',I4,' NEG. SHALL.,',I4,' DST.LT.0, ',I4,' DST.GT.0') ! 952 FORMAT(' LAYER (FROM TOP),',I2) ! 954 FORMAT(' NBOTS=',I4,' NTOPS=',I4,' NDPTHS=',I4,' NBOTD=',I4,' NTOPD=',I4,' NDPTHD=',I4) !--------------------------------------------------------------------------------------------------------------------- RETURN ! END SUBROUTINE CUCNVC ! !-------------------diego----------------------------------------------- !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc FUNCTION GAMMA(X) !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !D DOUBLE PRECISION FUNCTION DGAMMA(X) ! !OBS: CHANGED KIND TO SIMPLE PRECISION (DIEGO) !---------------------------------------------------------------------- ! ! THIS ROUTINE CALCULATES THE GAMMA FUNCTION FOR A REAL ARGUMENT X. ! COMPUTATION IS BASED ON AN ALGORITHM OUTLINED IN REFERENCE 1. ! THE PROGRAM USES RATIONAL FUNCTIONS THAT APPROXIMATE THE GAMMA ! FUNCTION TO AT LEAST 20 SIGNIFICANT DECIMAL DIGITS. COEFFICIENTS ! FOR THE APPROXIMATION OVER THE INTERVAL (1,2) ARE UNPUBLISHED. ! THOSE FOR THE APPROXIMATION FOR X .GE. 12 ARE FROM REFERENCE 2. ! THE ACCURACY ACHIEVED DEPENDS ON THE ARITHMETIC SYSTEM, THE ! COMPILER, THE INTRINSIC FUNCTIONS, AND PROPER SELECTION OF THE ! MACHINE-DEPENDENT CONSTANTS. ! ! !******************************************************************* !******************************************************************* ! ! EXPLANATION OF MACHINE-DEPENDENT CONSTANTS ! ! BETA - RADIX FOR THE FLOATING-POINT REPRESENTATION ! MAXEXP - THE SMALLEST POSITIVE POWER OF BETA THAT OVERFLOWS ! XBIG - THE LARGEST ARGUMENT FOR WHICH GAMMA(X) IS REPRESENTABLE ! IN THE MACHINE, I.E., THE SOLUTION TO THE EQUATION ! GAMMA(XBIG) = BETA**MAXEXP ! XINF - THE LARGEST MACHINE REPRESENTABLE FLOATING-POINT NUMBER; ! APPROXIMATELY BETA**MAXEXP ! EPS - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT ! 1.0+EPS .GT. 1.0 ! XMININ - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT ! 1/XMININ IS MACHINE REPRESENTABLE ! ! APPROXIMATE VALUES FOR SOME IMPORTANT MACHINES ARE: ! ! BETA MAXEXP XBIG ! ! CRAY-1 (S.P.) 2 8191 966.961 ! CYBER 180/855 ! UNDER NOS (S.P.) 2 1070 177.803 ! IEEE (IBM/XT, ! SUN, ETC.) (S.P.) 2 128 35.040 ! IEEE (IBM/XT, ! SUN, ETC.) (D.P.) 2 1024 171.624 ! IBM 3033 (D.P.) 16 63 57.574 ! VAX D-FORMAT (D.P.) 2 127 34.844 ! VAX G-FORMAT (D.P.) 2 1023 171.489 ! ! XINF EPS XMININ ! ! CRAY-1 (S.P.) 5.45E+2465 7.11E-15 1.84E-2466 ! CYBER 180/855 ! UNDER NOS (S.P.) 1.26E+322 3.55E-15 3.14E-294 ! IEEE (IBM/XT, ! SUN, ETC.) (S.P.) 3.40E+38 1.19E-7 1.18E-38 ! IEEE (IBM/XT, ! SUN, ETC.) (D.P.) 1.79D+308 2.22D-16 2.23D-308 ! IBM 3033 (D.P.) 7.23D+75 2.22D-16 1.39D-76 ! VAX D-FORMAT (D.P.) 1.70D+38 1.39D-17 5.88D-39 ! VAX G-FORMAT (D.P.) 8.98D+307 1.11D-16 1.12D-308 ! !******************************************************************* !******************************************************************* ! ! ERROR RETURNS ! ! THE PROGRAM RETURNS THE VALUE XINF FOR SINGULARITIES OR ! WHEN OVERFLOW WOULD OCCUR. THE COMPUTATION IS BELIEVED ! TO BE FREE OF UNDERFLOW AND OVERFLOW. ! ! ! INTRINSIC FUNCTIONS REQUIRED ARE: ! ! INT, DBLE, EXP, LOG, REAL, SIN ! ! ! REFERENCES: AN OVERVIEW OF SOFTWARE DEVELOPMENT FOR SPECIAL ! FUNCTIONS W. J. CODY, LECTURE NOTES IN MATHEMATICS, ! 506, NUMERICAL ANALYSIS DUNDEE, 1975, G. A. WATSON ! (ED.), SPRINGER VERLAG, BERLIN, 1976. ! ! COMPUTER APPROXIMATIONS, HART, ET. AL., WILEY AND ! SONS, NEW YORK, 1968. ! ! LATEST MODIFICATION: OCTOBER 12, 1989 ! ! AUTHORS: W. J. CODY AND L. STOLTZ ! APPLIED MATHEMATICS DIVISION ! ARGONNE NATIONAL LABORATORY ! ARGONNE, IL 60439 ! !---------------------------------------------------------------------- USE F77KINDS !diego ! INTEGER I,N LOGICAL PARITY REAL(KIND=R4KIND) :: GAMMA REAL(KIND=R4KIND) & !D DOUBLE PRECISION C,CONV,EPS,FACT,HALF,ONE,P,PI,Q,RES,SQRTPI,SUM,TWELVE, & TWO,X,XBIG,XDEN,XINF,XMININ,XNUM,Y,Y1,YSQ,Z,ZERO DIMENSION C(7),P(8),Q(8) !---------------------------------------------------------------------- ! MATHEMATICAL CONSTANTS !---------------------------------------------------------------------- DATA ONE,HALF,TWELVE,TWO,ZERO/1.0E0_R4KIND,0.5E0_R4KIND,12.0E0_R4KIND,2.0E0_R4KIND,0.0E0_R4KIND/, & SQRTPI/0.9189385332046727417803297E0_R4KIND/, & PI/3.1415926535897932384626434E0_R4KIND/ !D DATA ONE,HALF,TWELVE,TWO,ZERO/1.0D0,0.5D0,12.0D0,2.0D0,0.0D0/, !D 1 SQRTPI/0.9189385332046727417803297D0/, !D 2 PI/3.1415926535897932384626434D0/ !---------------------------------------------------------------------- ! MACHINE DEPENDENT PARAMETERS !---------------------------------------------------------------------- DATA XBIG,XMININ,EPS/35.040E0_R4KIND,1.18E-38_R4KIND,1.19E-7_R4KIND/, & XINF/3.4E38_R4KIND/ !D DATA XBIG,XMININ,EPS/171.624D0,2.23D-308,2.22D-16/, !D 1 XINF/1.79D308/ !---------------------------------------------------------------------- ! NUMERATOR AND DENOMINATOR COEFFICIENTS FOR RATIONAL MINIMAX ! APPROXIMATION OVER (1,2). !---------------------------------------------------------------------- DATA P/-1.71618513886549492533811E+0_R4KIND,2.47656508055759199108314E+1_R4KIND,& -3.79804256470945635097577E+2_R4KIND,6.29331155312818442661052E+2_R4KIND,& 8.66966202790413211295064E+2_R4KIND,-3.14512729688483675254357E+4_R4KIND,& -3.61444134186911729807069E+4_R4KIND,6.64561438202405440627855E+4_R4KIND/ DATA Q/-3.08402300119738975254353E+1_R4KIND,3.15350626979604161529144E+2_R4KIND,& -1.01515636749021914166146E+3_R4KIND,-3.10777167157231109440444E+3_R4KIND,& 2.25381184209801510330112E+4_R4KIND,4.75584627752788110767815E+3_R4KIND,& -1.34659959864969306392456E+5_R4KIND,-1.15132259675553483497211E+5_R4KIND/ !D DATA P/-1.71618513886549492533811D+0,2.47656508055759199108314D+1, !D 1 -3.79804256470945635097577D+2,6.29331155312818442661052D+2, !D 2 8.66966202790413211295064D+2,-3.14512729688483675254357D+4, !D 3 -3.61444134186911729807069D+4,6.64561438202405440627855D+4/ !D DATA Q/-3.08402300119738975254353D+1,3.15350626979604161529144D+2, !D 1 -1.01515636749021914166146D+3,-3.10777167157231109440444D+3, !D 2 2.25381184209801510330112D+4,4.75584627752788110767815D+3, !D 3 -1.34659959864969306392456D+5,-1.15132259675553483497211D+5/ !---------------------------------------------------------------------- ! COEFFICIENTS FOR MINIMAX APPROXIMATION OVER (12, INF). !---------------------------------------------------------------------- DATA C/-1.910444077728E-03_R4KIND,8.4171387781295E-04_R4KIND, & -5.952379913043012E-04_R4KIND,7.93650793500350248E-04_R4KIND,& -2.777777777777681622553E-03_R4KIND,8.333333333333333331554247E-02_R4KIND,& 5.7083835261E-03_R4KIND/ !D DATA C/-1.910444077728D-03,8.4171387781295D-04, !D 1 -5.952379913043012D-04,7.93650793500350248D-04, !D 2 -2.777777777777681622553D-03,8.333333333333333331554247D-02, !D 3 5.7083835261D-03/ !---------------------------------------------------------------------- ! STATEMENT FUNCTIONS FOR CONVERSION BETWEEN INTEGER AND FLOAT !---------------------------------------------------------------------- CONV(I) = REAL(I,R4KIND) !D CONV(I) = DBLE(I) PARITY=.FALSE. FACT=ONE N=0 Y=X IF(Y.LE.ZERO)THEN !---------------------------------------------------------------------- ! ARGUMENT IS NEGATIVE !---------------------------------------------------------------------- Y=-X Y1=AINT(Y) RES=Y-Y1 IF(RES.NE.ZERO)THEN IF(Y1.NE.AINT(Y1*HALF)*TWO)PARITY=.TRUE. FACT=-PI/SIN(PI*RES) Y=Y+ONE ELSE RES=XINF GOTO 900 ENDIF ENDIF !---------------------------------------------------------------------- ! ARGUMENT IS POSITIVE !---------------------------------------------------------------------- IF(Y.LT.EPS)THEN !---------------------------------------------------------------------- ! ARGUMENT .LT. EPS !---------------------------------------------------------------------- IF(Y.GE.XMININ)THEN RES=ONE/Y ELSE RES=XINF GOTO 900 ENDIF ELSEIF(Y.LT.TWELVE)THEN Y1=Y IF(Y.LT.ONE)THEN !---------------------------------------------------------------------- ! 0.0 .LT. ARGUMENT .LT. 1.0 !---------------------------------------------------------------------- Z=Y Y=Y+ONE ELSE !---------------------------------------------------------------------- ! 1.0 .LT. ARGUMENT .LT. 12.0, REDUCE ARGUMENT IF NECESSARY !---------------------------------------------------------------------- N=INT(Y)-1 Y=Y-CONV(N) Z=Y-ONE ENDIF !---------------------------------------------------------------------- ! EVALUATE APPROXIMATION FOR 1.0 .LT. ARGUMENT .LT. 2.0 !---------------------------------------------------------------------- XNUM=ZERO XDEN=ONE DO 260 I=1,8 XNUM=(XNUM+P(I))*Z XDEN=XDEN*Z+Q(I) 260 CONTINUE RES=XNUM/XDEN+ONE IF(Y1.LT.Y)THEN !---------------------------------------------------------------------- ! ADJUST RESULT FOR CASE 0.0 .LT. ARGUMENT .LT. 1.0 !---------------------------------------------------------------------- RES=RES/Y1 ELSEIF(Y1.GT.Y)THEN !---------------------------------------------------------------------- ! ADJUST RESULT FOR CASE 2.0 .LT. ARGUMENT .LT. 12.0 !---------------------------------------------------------------------- DO 290 I=1,N RES=RES*Y Y=Y+ONE 290 CONTINUE ENDIF ELSE !---------------------------------------------------------------------- ! EVALUATE FOR ARGUMENT .GE. 12.0, !---------------------------------------------------------------------- IF(Y.LE.XBIG)THEN YSQ=Y*Y SUM=C(7) DO 350 I=1,6 SUM=SUM/YSQ+C(I) 350 CONTINUE SUM=SUM/Y-Y+SQRTPI SUM=SUM+(Y-HALF)*LOG(Y) RES=EXP(SUM) ELSE RES=XINF GOTO 900 ENDIF ENDIF !---------------------------------------------------------------------- ! FINAL ADJUSTMENTS AND RETURN !---------------------------------------------------------------------- IF(PARITY)RES=-RES IF(FACT.NE.ONE)RES=FACT/RES 900 GAMMA=RES !D900 DGAMMA = RES RETURN ! ---------- LAST LINE OF GAMMA ---------- END FUNCTION GAMMA !-------------------diego----------------------------------------------