SUBROUTINE ADJPPT
C ******************************************************************
C$$$ SUBPROGRAM DOCUMENTATION BLOCK
C . . .
C SUBPROGRAM: ADJPPT PRECIPITATION ASSIMILATION ADJUSTMENT
C PRGRMMR: BALDWIN ORG: W/NP22 DATE: 98-08-26
C
C ABSTRACT:
C ADJPPT MODIFIES THE ENVIRONMENT AND RELEASES LATENT HEAT
C TO ATTEMPT TO MATCH THE AMOUNT OF OBSERVED PRECIPITATION.
C IF THE OBSERVED PRECIPITATION IS DETERMINED TO BE CONVECTIVE IN TYPE,
C THE LATENT HEAT RELEASE AND MOISTURE CHANGE IS FORMULATED TO FOLLOW
C THE BETTS-MILLER-JANJIC SCHEME. OTHERWISE ADJUSTMENTS ARE MADE
C IN A MANNER CONSISTENT WITH THE GRID-SCALE CLOUD PHYSICS SCHEME.
C
C
C PROGRAM HISTORY LOG:
C 98-08-26 BALDWIN - ORIGINATOR
C
C USAGE: CALL ADJPPT FROM MAIN PROGRAM EBU
C
C INPUT ARGUMENT LIST:
C NONE
C
C OUTPUT ARGUMENT LIST:
C NONE
C
C OUTPUT FILES:
C NONE
C
C SUBPROGRAMS CALLED:
C
C UNIQUE:
C TTBLEX
C
C LIBRARY:
C NONE
C
C COMMON BLOCKS: CTLBLK
C LOOPS
C MASKS
C PHYS
C VRBLS
C CNVCLD
C PVRBLS
C ACMCLH
C INDX
C
C ATTRIBUTES:
C LANGUAGE: FORTRAN 90
C MACHINE : IBM SP
C$$$
C----------------------------------------------------------------------
INCLUDE "cuparm"
INCLUDE "parmeta"
INCLUDE "parm.tbl"
INCLUDE "mpp.h"
C----------------------------------------------------------------------
P A R A M E T E R
& (IMJM=IM*JM,JAM=6+2*(JM-10)
&, IMJM_LOC=IDIM2*JDIM2
&, LP1=LM+1,LM1=LM-1)
PARAMETER (ELIV=2.834E6, DLDT=2274.E0)
C----------------------------------------------------------------------
L O G I C A L
& RUN,FIRST,RESTRT,SIGMA
C----------------------------------------------------------------------
INCLUDE "CTLBLK.comm"
C-----------------------------------------------------------------------
INCLUDE "LOOPS.comm"
C-----------------------------------------------------------------------
INCLUDE "MASKS.comm"
C-----------------------------------------------------------------------
INCLUDE "PHYS.comm"
C-----------------------------------------------------------------------
INCLUDE "VRBLS.comm"
C-----------------------------------------------------------------------
INCLUDE "CNVCLD.comm"
C-----------------------------------------------------------------------
INCLUDE "PVRBLS.comm"
C-----------------------------------------------------------------------
INCLUDE "ACMCLH.comm"
C-----------------------------------------------------------------------
INCLUDE "INDX.comm"
C-----------------------------------------------------------------------
INCLUDE "PPTASM.comm"
C-----------------------------------------------------------------------
INCLUDE "CLDWTR.comm"
C-----------------------------------------------------------------------
D I M E N S I O N
& TREFK (LM),QREFK (LM),PK (LM),APEK (LM),TK (LM)
&,THSK (LM),PSK (LM),APESK (LM),QK (LM),THERK (LM)
&,THVREF(LM),THEVRF(LM),THVMOD(LM),DIFT (LM),DIFQ (LM)
&,QSATK (LM),FPK (LM),RELH (LM)
C
D I M E N S I O N
& LTOP (IDIM1:IDIM2,JDIM1:JDIM2),LBOT (IDIM1:IDIM2,JDIM1:JDIM2)
&,PTOP (IDIM1:IDIM2,JDIM1:JDIM2),PBOT (IDIM1:IDIM2,JDIM1:JDIM2)
&,IPTB (IDIM1:IDIM2,JDIM1:JDIM2),ITHTB (IDIM1:IDIM2,JDIM1:JDIM2)
&,PDSL (IDIM1:IDIM2,JDIM1:JDIM2),APEBT (IDIM1:IDIM2,JDIM1:JDIM2)
&,PRECL(IDIM1:IDIM2,JDIM1:JDIM2), QC(LM)
&,TBT (IDIM1:IDIM2,JDIM1:JDIM2),Q2BT (IDIM1:IDIM2,JDIM1:JDIM2)
&,QQ (IDIM1:IDIM2,JDIM1:JDIM2),PP (IDIM1:IDIM2,JDIM1:JDIM2)
&,PSP (IDIM1:IDIM2,JDIM1:JDIM2),THBT (IDIM1:IDIM2,JDIM1:JDIM2)
&,THESP (IDIM1:IDIM2,JDIM1:JDIM2),P (IDIM1:IDIM2,JDIM1:JDIM2)
&,ILRES (IMJM_LOC),JLRES (IMJM_LOC)
&,IHRES (IMJM_LOC),JHRES (IMJM_LOC)
&,DDATA(IDIM1:IDIM2,JDIM1:JDIM2), ADATA(IDIM1:IDIM2,JDIM1:JDIM2)
&,APE (IDIM1:IDIM2,JDIM1:JDIM2,LM)
&,TREF (IDIM1:IDIM2,JDIM1:JDIM2,LM)
C
C-----------------------------------------------------------------------
C CPREC: model convective precip at each time step.
C PREC = CPREC + APREC
C
C DDATA: Pobs at each time step
C
C If observed precip is larger than model precip, give the
C convective adjustment the first chance to make the rain. (Initially
C I wanted to keep the conv/gridscale precip invariant through ADJPPT,
C but then thought better of it - in the case of Pobs >> Pmod, the
C instability might build up too much if we don't let the conv adj
C take care of as much of it as possible).
C
C ADATA: After convective adjustment, let ADATA = DDATA - CPREC
C (i.e. ADATA is the portion of the observed precip not able to
C be accounted for by convective adjustment, and to be accounted
C for in the grid-scale portion of the adjustment scheme).
C 'CPREC' is not an actual variable in this routine.
C
C PPTSUM is used to keep track of total observed precip at grid point
C (itest,jtest) during the 3-hour assimilation period (to make sure that
C we are partitioning the hourly precip obs correctly). PPTSUM is zero
C before the first call to ADJPPT, and the accumulated value is saved
C between subsequent calls.
C
LOGICAL IGSADJ(IDIM1:IDIM2,JDIM1:JDIM2)
REAL PFACTOR(IDIM1:IDIM2,JDIM1:JDIM2)
C
C IGSADJ: flag for whether to make grid-scale adjustment
C (T - do GS adjustment. F - don't do GS adjustment)
C IGSADJ is false if
C 1) PPTDAT = 999., or
C 2) DDATA = 0., or
C 3) DDATA .LE. PREC, or
C 4) While DDATA > PREC, all (pre-adj) Pmod is convective, and
C the entire amount of DDATA is accounted for during convective
C adjustment
C IGSADJ is true if the convective adjustment does not
C produce enough convective precip to account for DDATA
C
DATA PPTSUM/0./
SAVE PPTSUM
C
C----------------------------------------------------------------------
IF (NTSD.EQ.1) THEN
TLAT=0.
RETURN
ENDIF
IF (NTSD.GT.NHEAT) RETURN
C
C-----------------------------------------------------------------------
C--------------INITIALIZE GRID-SCALE ADJUSTMENT MASK--------------------
C--------------(later the mask will be updated during conv adj)---------
C-----------------------------------------------------------------------
C
IGSADJ = .FALSE.
C
C PREPATORY CALCULATIONS
C-----------------------------------------------------------------------
DTCNVC=NCNVC*DT
RDTCNVC=1./DTCNVC
TAUK=DTCNVC/TREL
CTHRS=(0.006350/86400.)*DTCNVC
TIMES=(NTSD-1)*DT
COUNT=MOD(TIMES,3600.)
IHR=(TIMES-1.0)/3600.+1
PHYST=DTCNVC
C
IF (COUNT.GE.PHYST.OR.COUNT.EQ.0.0) THEN
FRACT=PHYST/3600.
ELSE
FRACT=COUNT/3600.
END IF
C
c
c Check to see if this is the last time step before the end. If so,
c applying the remaining fraction of precip to this time step.
c
c No, don't do that. For the 80km runs, this would mean an 25% increase
c of precip in the last time step, and it'll affect the temperature and
c the moisture fields inappriately. Better to increase the model precip
c after precip assim adjustment by 25% in the end (i.e. multiply by
c 'ENDFCTR').
c
TEND=3.0
c print*,'tend=',tend
IF (TEND*3600.-TIMES .LT. PHYST) THEN
c FEND = TEND - TIMES/3600.
ENDFCTR = (TEND*3600-TIMES+PHYST)/PHYST
ELSE
c FEND = 0.
ENDFCTR = 1.
ENDIF
c
C
IF (TIMES .LT. PHYST) THEN
ZER=1.0E-05*FRACT
ELSE
ZER=1.0E-05*PHYST/3600.
ENDIF
C
C ZER IS OUR ZERO THRESHOLD; .01 MM PER HOUR
C (CORRESPONDS TO 1 HUNDRETH OF AN INCH PER DAY)
C
SIXSIX=PHYST/3600.
C
C Under one of the scenarios (when Pobs > 0 and Pmod=0), we need to
C create a layer of precipitating cloud from scratch. We specify 3
C cloud-thicknesses based on the precipitation amount:
C
PTRES1=2.81E-03*SIXSIX
PTRES2=3.75E-04*SIXSIX
PTRES3=1.0E-03*SIXSIX
C
c print*,'mype,mtstpe=',mype,mtstpe
c print*,'itstloc,jtstloc=',itstloc,jtstloc
IF (MYPE.EQ.MTSTPE) THEN
WRITE(98,*) 'NTSD=',NTSD,' TIMES=',TIMES,' FRACT=',FRACT,
& ' ENDFCTR=', ENDFCTR
WRITE(98,*) 'IHR=', IHR,' PPTDAT=',PPTDAT(ITSTLOC,JTSTLOC,IHR)
c print*,'pptdat(itstloc,jtstloc,ihr)=',pptdat(itstloc,jtstloc,ihr)
ENDIF
C-----------------------------------------------------------------------
C FRACT IS THE FRACTION OF IHR'S PRECIP THAT WE WANT FOR
C THIS ADJUSTMENT, WE WANT (PHYST/3600-FRACT) WORTH OF IHR-1 PRECIP
C WE HAVE DATA ONLY FOR IHR=1,3
C-----------------------------------------------------------------------
C SET UP OBSERVED PRECIP FOR THIS TIMESTEP IN DDATA
C-----------------------------------------------------------------------
C
PFACTOR = 1.
C
!$omp parallel do private(pdiff,pexp,pptsum)
DO 110 J=MYJS,MYJE
DO 100 I=MYIS,MYIE
C---
R2D=57.2957795 ! 180.0/PI
GLATMIN=27.5
GLATMAX=42.5
GLATD=GLAT(I,J)*R2D
IF (GLATD.GE.GLATMAX .AND. (SM(I,J)+SICE(I,J)).GT.0.5) THEN
PPTDAT(I,J,IHR)=999.0
END IF
C---
IF (PPTDAT(I,J,IHR).GT.900.) GO TO 100
C---- rewrite 12-11 WNE
C IF (IHR.EQ.1 .OR. PPTDAT(I,J,IHR-1).GT.900.) THEN
IF (.not.(IHR.NE.1 .AND. PPTDAT(I,J,IHR-1).LE.900.)) THEN
DDATA(I,J) = PPTDAT(I,J,IHR)*FRACT
ELSE
DDATA(I,J) = PPTDAT(I,J,IHR)*FRACT
& + PPTDAT(I,J,IHR-1)*(SIXSIX-FRACT)
c & + PPTDAT(I,J,3)*FEND
ENDIF
C---
CC IF ((SM(I,J)+SICE(I,J)).GT.0.5) THEN
IF (SM(I,J).GT.0.5) THEN
IF (GLATD.GT.GLATMIN .AND. GLATD.LT.GLATMAX) THEN
AFAC=1.0-((GLATD-GLATMIN)/(GLATMAX-GLATMIN))
DDATA(I,J) = AFAC*DDATA(I,J) + (1.0-AFAC)*PREC(I,J)
END IF
END IF
C---
C
C Use the difference between Pobs and Pmod to modify RH in the cloud
C (M. Baldwin, 20 Apr 99) by a factor of
C 1.0+0.2*(exp(r)-exp(-r))/(exp(r)+exp(-r)) where: r=(Pobs-Pmod in mm)/25mm
C (Pobs and Pmod are hourly precip). The value of the factor would be
C between 0.8 and 1.2.
C
CC The moisture field at T0 seems to be too wet. So only use PFACTOR to
CC dry the atmosphere if Pobs < Pmod. Otherwise set PFACTOR to 1.
CC No, don't do that ... for now.
PDIFF = (DDATA(I,J)-PREC(I,J))/(0.025*FRACT)
CC PDIFF = AMIN1(0.,PDIFF)
PEXP = EXP(PDIFF)
PEXP = PEXP * PEXP
PFACTOR(I,J) = 1.0 + 0.2 * (PEXP-1.)/(PEXP+1.)
C
C
C If PREC > 0 (i.e. Pmod > 0), partition DDATA into 'convective'
C and 'grid scale', based on the ratio of APREC/PREC.
C
C IF Pmod = 0, first assume that all observed precip are in fact
C convective (we will try to let convective adjustment take care
C of it first. If there's any leftover DDATA un-accounted for,
C we then let grid-scale precip take care of it.
C
IF (I.EQ.ITSTLOC .AND. J.EQ.JTSTLOC .AND. MYPE.EQ.MTSTPE) THEN
c print*,'should be writing to unit 98'
PPTSUM = PPTSUM + DDATA(I,J)*ENDFCTR
WRITE(98,*) 'DDATA=',DDATA(I,J),' PREC=',PREC(I,J),
& ' APREC=', APREC(I,J), ' ZER=',ZER,' PPTSUM=', PPTSUM
WRITE(98,*) 'PFACTOR=', PFACTOR(I,J), ' PDIFF=', PDIFF,
& ' PEXP=', PEXP
c print*,'i,j,ddata(i,j),prec(i,j),pptsum=',
c * i,j,ddata(i,j),prec(i,j),pptsum
ENDIF
100 CONTINUE
110 CONTINUE
C
C Set minimum cloud depth for deep convection. This would be scaled
C by the total atmosphere depth (PSFCIJ) at this horizontal point later on.
C
PSHNEW=20000.
C
C The big loop - looping through all horizontal grid points
C In M. Baldwin's original ADJPPT.f, there was no 'big (i,j) loop'.
C I am replacing the many little loops in his code with this big loop.
C The upper and lower limits of I and J are chosen to be this way to
C be consistent with his loop limits in the buoyancy calculation.
C
!$omp parallel do
!$omp& private(adjust,ai,apekl,apekxx,apekxy,apesp,apests,
!$omp& bi,bq,bqs00k,bqs10k,climit,cratio,
!$omp& delt,deltacp,delcwm,delq,depmin,depth,depwl,detacl,
!$omp& dsp,dsp0k,dspbk,dsptk,dthem,
!$omp& efi,efinew,elv,etabot,etatop,etbig,
!$omp& factor,fi,fiw,fiwl1,
!$omp& iq,iqtb,it,itb,ittb,ittbk,ivi,
!$omp& knuml,knunh,l0,l0m1,lb,lbtk,lbm1,lcbottm,
!$omp& lmhij,lmhk,ltp1,ltpk,numlev,oldcwm,oldq,oldrh,
!$omp& p00k,p01k,p10k,p11k,petal,pk0,pkb,pkl,pkt,
!$omp& pp1,prec1,preck,precmax,presk,psfc,psfck,
!$omp& qbt,qckl,qkl,qq1,qi,qint,qw,
!$omp& ratio,rdp0t,relhum,rhfctr,sq,sqs00k,sqs10k,stabdl,
!$omp& therkx,therky,tkl,tmt0,tmt15,tpsp,tq,trefkx,
!$omp& ttemp,tth,tthbt,tthes,wfix,wmin,yltmp)
DO 910 J=MYJS,MYJE
DO 900 I=MYIS,MYIE
IF (PPTDAT(I,J,IHR).GT.900. .OR.
& DDATA(I,J).LE.ZER .AND. PREC(I,J).LE.ZER) GOTO 900
C-----------------------------------------------------------------------
C--------------PREPARATIONS---------------------------------------------
C-----------------------------------------------------------------------
THESP(I,J)=0.
PDSL (I,J)=RES(I,J)*PD(I,J)
LBOT (I,J)=LMH(I,J)
PBOT(I,J)=AETA(LBOT(I,J))*PDSL(I,J)+PT
TREF(I,J,1)=T(I,J,1)
C-----------------------------------------------------------------------
C--- CASE 1. Pobs = 0, Pmod > 0
C--------------IF OBSERVED PRECIP IS LESS THAN OR EQUAL TO ZER----------
C--------------TAKE BACK THE LATENT HEAT RELEASE------------------------
C-----------------------------------------------------------------------
IF (DDATA(I,J).LE.ZER .AND. PREC(I,J).GT.ZER) THEN
CLDEFI(I,J)=STEFI
DO 130 L=1,LM
IF (HTM(I,J,L).LT.0.5) GO TO 130
C--------- -FIND THE PRE-MODIFIED RELATIVE HUMIDITY FOR THIS POINT------
PETAL=PDSL(I,J)*AETA(L)+PT
QCKL=PQ0/PETAL
2 *EXP(A2*(T(I,J,L)-A3)/(T(I,J,L)-A4))
if (QCKL.eq.0.0) then
write(0,*)"QCKL=0.0",QCKL,I,J,L,MYPE
write(0,*)PQ0,PETAL,T(I,J,L)
write(*,*)"QCKL=0.0",QCKL,I,J,L,MYPE
write(*,*)PQ0,PETAL,T(I,J,L)
CALL MPI_FINALIZE(IERR)
STOP 8
end if
RELHUM=Q(I,J,L)/QCKL
OLDRH=RELHUM
OLDQ=Q(I,J,L)
OLDCWM=CWM(I,J,L)
C MODIFY THE TEMP AND PRECIP
T(I,J,L)=T(I,J,L)-TLAT(I,J,L)
C Reduce RH by the factor PFACTOR:
QCKL=PQ0/PETAL
2 *EXP(A2*(T(I,J,L)-A3)/(T(I,J,L)-A4))
RELHUM= RELHUM * PFACTOR(I,J)
IF (TLAT(I,J,L).GT.0.) Q(I,J,L)=RELHUM*QCKL
IF (I.EQ.ITSTLOC .AND. J.EQ.JTSTLOC .AND. MYPE.EQ.MTSTPE)
2 WRITE(98,129) L, PETAL,TLAT(I,J,L), OLDRH, PFACTOR(I,J),
2 RELHUM, OLDQ, Q(I,J,L)
129 FORMAT(I2,X,f7.0,x,f7.4,X,3(2X,F6.4), X, 2(2X,E12.6))
C If any part of model-predicted rainfall was grid-scale, decrease the
C cloud water mixing ratio ( if > WMIN) to the minimum value, WMIN:
IF (APREC(I,J).GT.0.) THEN
TTEMP=0.025*(T(I,J,L)-273.16)
WFIX=0.9814*EXP(0.01873*L)
WMIN=0.1E-3*EXP(TTEMP)*WFIX
c if(mype.eq.17) then
c print*,'in 130 loop'
c print*,'i,j,l=',i,j,l
c print*,'cwm(i,j,l),wmin=',cwm(i,j,l),wmin
c print*,'ttemp=',ttemp
c print*,'t(i,j,l)=',t(i,j,l)
c endif
CWM(I,J,L) = AMIN1(WMIN,CWM(I,J,L))
ENDIF
c
c Calculate the water vapor and cloud water/ice increments for
c 1) the entire column
c 2) sfc-700mb
c
DELQ = (Q(I,J,L)-OLDQ) * DETA(L)*PDSL(I,J)/G
DELCWM = (CWM(I,J,L)-OLDCWM) * DETA(L)*PDSL(I,J)/G
c
VAPINC(I,J)=VAPINC(I,J)+DELQ
CLDINC(I,J)=CLDINC(I,J)+DELCWM
C
IF (PETAL.GE.70000.) THEN
VAPINC7(I,J)=VAPINC7(I,J)+DELQ
CLDINC7(I,J)=CLDINC7(I,J)+DELCWM
ENDIF
130 CONTINUE
C Take back the PREC from ACPREC and CUPREC as well. For CUPREC,
C the amount taken back depends how much convective precip there
C was before the adjustment.
ACPREC(I,J)=ACPREC(I,J)-PREC(I,J)
CUPREC(I,J)=CUPREC(I,J)-(PREC(I,J)-APREC(I,J))
CUPPT(I,J)= CUPPT(I,J)-(PREC(I,J)-APREC(I,J))
PREC(I,J)=0.
APREC(I,J)=0.
GO TO 900
ENDIF
C
C CASE 2, Pmod > Pobs > 0
C
IF (DDATA(I,J).LE.PREC(I,J)) THEN
C THIS IS THE ADJUSTMENT WE DO IF WE HAD TOO MUCH PRECIP, CONVECTIVE
C OR OTHERWISE, IN THE MODEL. MULTIPLY THE LATENT HEAT
C AT EACH LEVEL BY THE FRACTION: DATA/MODEL RAINFALL
C MATCH THE RH THAT THE PROFILE HAD PRIOR TO THIS ADJUSTMENT
C-----------FIND THE PRE-MODIFIED RELATIVE HUMIDITY FOR THIS POINT------
ADJUST=DDATA(I,J)/PREC(I,J)
if (i.eq.itstloc .and. j.eq.jtstloc .and. mype.eq.mtstpe)
& write(98,*)
& 'Check for Case 2, DDATA=', DDATA(I,J),' PREC=',PREC(I,J),
& ' ADJUST=', ADJUST, ' SR=', SR(I,J)
C Compute the ratio of convective precip/total precip:
CRATIO=(PREC(I,J)-APREC(I,J))/PREC(I,J)
C
DO 140 L=1,LM
IF (HTM(I,J,L).LT.0.5) GO TO 140
PETAL=PDSL(I,J)*AETA(L)+PT
c if(i.eq.16.and.j.eq.35.and.mype.eq.13.and.ntsd.eq.26) then
c print*,'i,j,l,aeta(l),pdsl(i,j)=',i,j,l,aeta(l),pdsl(i,j)
c endif
QCKL=PQ0/PETAL
2 *EXP(A2*(T(I,J,L)-A3)/(T(I,J,L)-A4))
RELHUM=Q(I,J,L)/QCKL
OLDRH=RELHUM
OLDQ=Q(I,J,L)
OLDCWM=CWM(I,J,L)
C MODIFY THE TEMP CHANGE AND PRECIP
T(I,J,L)=T(I,J,L)+TLAT(I,J,L)*(ADJUST-1.)
CYL
CYL Assume ice process below freezing, and water process otherwise.
CYL This is to be consistent with the latent heat calculation in PRECPD.
CYL
IF (T(I,J,L).GE.273.15) THEN
ELV=ELWV
ELSE
ELV=ELIV
ENDIF
DELT=TLAT(I,J,L)*(ADJUST-1.)
C The following is the 'accum pcp' version of DELT to take care of the
C fractional time step at the end of each EDAS segment:
DELTACP=TLAT(I,J,L)*(ADJUST*ENDFCTR-1.)
PREC(I,J)=DELT*DETA(L)*PDSL(I,J)*CP/(ROW*G*ELV)+PREC(I,J)
CUPREC(I,J)=DELTACP*DETA(L)*PDSL(I,J)*CP/(ROW*G*ELV)*CRATIO
2 +CUPREC(I,J)
ACPREC(I,J)=DELTACP*DETA(L)*PDSL(I,J)*CP/(ROW*G*ELV)
2 +ACPREC(I,J)
CUPPT(I,J)=DELTACP*DETA(L)*PDSL(I,J)*CP/(ROW*G*ELV)*CRATIO
2 +CUPPT(I,J)
C Reduce RH by the factor PFACTOR:
QCKL=PQ0/PETAL
2 *EXP(A2*(T(I,J,L)-A3)/(T(I,J,L)-A4))
RELHUM= RELHUM * PFACTOR(I,J)
IF (TLAT(I,J,L).GT.0.) Q(I,J,L)=RELHUM*QCKL
IF (I.EQ.ITSTLOC .AND. J.EQ.JTSTLOC .AND. MYPE.EQ.MTSTPE)
2 WRITE(98,128) L, PETAL,TLAT(I,J,L), T(I,J,L), OLDRH,
3 RELHUM, OLDQ, Q(I,J,L)
128 FORMAT(I2,X,f7.0,x,f7.4,X,F7.2,2(2X,F6.4), X, 2(2X,E12.6))
C
C If the model had grid-scale precip prior to the adjustment, and
C the cloud water is above the minimum (WMIN) for generating rain,
C reduce the cloud water(CWM) proportionally, but keep it above WMIN,
C
IF (TLAT(I,J,L).GT.0. .and. APREC(I,J).GT.0.
2 .and. CWM(I,J,L).GT. WMIN) THEN
TTEMP=0.025*(T(I,J,L)-273.16)
WFIX=0.9814*EXP(0.01873*L)
WMIN=0.1E-3*EXP(TTEMP)*WFIX
c if(mype.eq.17) then
c print*,'in 140 loop'
c print*,'i,j,l=',i,j,l
c print*,'cwm(i,j,l),wmin=',cwm(i,j,l),wmin
c print*,'ttemp=',ttemp
c print*,'t(i,j,l)=',t(i,j,l)
c print*,'adjust=',adjust
c endif
CWM(I,J,L) = AMAX1(WMIN,CWM(I,J,L)*ADJUST)
ENDIF
c
c Calculate the water vapor and cloud water/ice increments for
c 1) the entire column
c 2) sfc-700mb
c
DELQ = (Q(I,J,L)-OLDQ) * DETA(L)*PDSL(I,J)/G
DELCWM = (CWM(I,J,L)-OLDCWM) * DETA(L)*PDSL(I,J)/G
c
VAPINC(I,J)=VAPINC(I,J)+DELQ
CLDINC(I,J)=CLDINC(I,J)+DELCWM
C
IF (PETAL.GE.70000.) THEN
VAPINC7(I,J)=VAPINC7(I,J)+DELQ
CLDINC7(I,J)=CLDINC7(I,J)+DELCWM
ENDIF
140 CONTINUE
C
C We didn't adjust APREC yet. Here we'll reduce APREC proportionally:
APREC(I,J) = PREC(I,J) * (1.-CRATIO)
GO TO 900
ENDIF
C
C-----------------------------------------------------------------------
C Case 3 ------IF WE ARE HERE, THEN Pmod < Pobs ------------------------
C--------------IF OBSERVED PRECIP IS GREATER THAN ZER-------------------
C--------------DETERMINE IF IT IS CONVECTIVE OR GRID-SCALE--------------
C-----------------------------------------------------------------------
C GO THROUGH THE BETTS/MILLER/JANJIC CLOUD SEARCH, IF THE CLOUD
C IS CONSIDERED DEEP, THEN WE HAVE CONVECTION.
C-----------------------------------------------------------------------
C--------------PADDING SPECIFIC HUMIDITY IF TOO SMALL-------------------
C RESTORE APE TO SCRATCH ARRAY
DO 150 L=1,LM
APESTS=PDSL(I,J)*AETA(L)+PT
APE(I,J,L)=(1.E5/APESTS)**CAPA
IF(Q(I,J,L).LT.EPSQ)Q(I,J,L)=HTM(I,J,L)*EPSQ
150 CONTINUE
C--------------SEARCH FOR MAXIMUM BUOYANCY LEVEL------------------------
DO 170 KB=1,LM
IF (HTM(I,J,L).LT.0.5) GO TO 170
C--------------TRIAL MAXIMUM BUOYANCY LEVEL VARIABLES-------------------
PKL=AETA(KB)*PDSL(I,J)+PT
LMHK=LMH(I,J)
PSFCK=AETA(LMHK)*PDSL(I,J)+PT
C
IF(KB.LE.LMHK .AND. PKL.GE.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
C--------------KEEPING INDICES WITHIN THE TABLE-------------------------
IF(ITTB.LT.1)THEN
ITTB=1
QQ1=0.
ENDIF
C
IF(ITTB.GE.JTB)THEN
ITTB=JTB-1
QQ1=0.
ENDIF
C--------------BASE AND SCALING FACTOR FOR SPEC. HUMIDITY---------------
ITTBK=ITTB
BQS00K=QS0(ITTBK)
SQS00K=SQS(ITTBK)
BQS10K=QS0(ITTBK+1)
SQS10K=SQS(ITTBK+1)
C--------------SCALING SPEC. HUMIDITY & 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
C--------------KEEPING INDICES WITHIN THE TABLE-------------------------
IF(IQTB.LT.1)THEN
IQTB=1
PP1=0.
ENDIF
C
IF(IQTB.GE.ITB)THEN
IQTB=ITB-1
PP1=0.
ENDIF
C--------------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)
C--------------SATURATION POINT VARIABLES AT THE BOTTOM-----------------
TPSP=P00K+(P10K-P00K)*PP1+(P01K-P00K)*QQ1
1 +(P00K-P10K-P01K+P11K)*PP1*QQ1
APESP=(1.E5/TPSP)**CAPA
TTHES=TTHBT*EXP(ELOCP*QBT*APESP/TTHBT)
C--------------CHECK FOR MAXIMUM BUOYANCY-------------------------------
IF(TTHES.GT.THESP(I,J))THEN
PSP (I,J)=TPSP
THBT (I,J)=TTHBT
THESP(I,J)=TTHES
ENDIF
ENDIF
C-----------------------------------------------------------------------
170 CONTINUE
C
C---------CHOOSE CLOUD BASE AS MODEL LEVEL JUST BELOW PSP--------------
C
DO 190 L=1,LM1
IF (HTM(I,J,L).LT.0.5) GO TO 190
AETAL=AETA(L)
P(I,J)=PDSL(I,J)*AETAL+PT
IF(P(I,J).LT.PSP(I,J).AND.P(I,J).GE.PQM)LBOT(I,J)=L+1
190 CONTINUE
C***
C*** WARNING: LBOT MUST NOT BE GT LMH(I,J)-1 IN SHALLOW CONVECTION
C*** MAKE SURE CLOUD BASE IS AT LEAST PONE ABOVE THE SURFACE
C***
LMHIJ=LMH(I,J)
PBOT(I,J)=AETA(LBOT(I,J))*PDSL(I,J)+PT
PSFCK=AETA(LMHIJ)*PDSL(I,J)+PT
C
IF(PBOT(I,J).GE.PSFCK-PONE.OR.LBOT(I,J).GE.LMHIJ)THEN
DO 200 L=1,LMHIJ-1
IF (HTM(I,J,L).LT.0.5) GO TO 200
P(I,J)=AETA(L)*PDSL(I,J)+PT
IF(P(I,J).LT.PSFCK-PONE)LBOT(I,J)=L
200 CONTINUE
PBOT(I,J)=AETA(LBOT(I,J))*PDSL(I,J)+PT
ENDIF
C--------------CLOUD TOP COMPUTATION------------------------------------
LTOP(I,J)=LBOT(I,J)
PTOP(I,J)=PBOT(I,J)
C-----------------------------------------------------------------------
c
DO 250 L=LM,1,-1
IF (HTM(I,J,L).LT.0.5) GO TO 250
c
C--------------SCALING PRESSURE & TT TABLE INDEX------------------------
KNUML=0
KNUMH=0
C
PRESK=PDSL(I,J)*AETA(L)+PT
C
IF(PRESK.LT.PLQ)THEN
KNUML=KNUML+1
ILRES(KNUML)=I
JLRES(KNUML)=J
ELSE
KNUMH=KNUMH+1
IHRES(KNUMH)=I
JHRES(KNUMH)=J
ENDIF
C***
C*** COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE<PL
C**
IF(KNUML.GT.0)THEN
CALL TTBLEX(TREF(IDIM1,JDIM1,L),TTBL,ITB,JTB,KNUML
1, ILRES,JLRES,PDSL,AETA(L),HTM(IDIM1,JDIM1,L)
2, PT,PL,QQ(IDIM1,JDIM1),PP(IDIM1,JDIM1)
3, RDP,THE0,STHE,RDTHE
4, THESP(IDIM1,JDIM1),IPTB(IDIM1,JDIM1)
5, ITHTB(IDIM1,JDIM1))
ENDIF
C***
C*** COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE>PL
C**
IF(KNUMH.GT.0)THEN
CALL TTBLEX(TREF(IDIM1,JDIM1,L),TTBLQ,ITBQ,JTBQ,KNUMH
1, IHRES,JHRES,PDSL,AETA(L),HTM(IDIM1,JDIM1,L)
2, PT,PLQ,QQ(IDIM1,JDIM1),PP(IDIM1,JDIM1)
3, RDPQ,THE0Q,STHEQ,RDTHEQ
4, THESP(IDIM1,JDIM1),IPTB(IDIM1,JDIM1)
5, ITHTB(IDIM1,JDIM1))
ENDIF
250 CONTINUE
C--------------BUOYANCY CHECK-------------------------------------------
DO 280 L=LM,1,-1
IF (HTM(I,J,L).LT.0.5) GO TO 280
IF(TREF(I,J,L).GT.T(I,J,L)-DTTOP)LTOP(I,J)=L
280 CONTINUE
C-----------------CLOUD TOP PRESSURE------------------------------------
PTOP(I,J)=AETA(LTOP(I,J))*PDSL(I,J)+PT
C--------------CLEAN UP AND GATHER DEEP CONVECTION POINTS---------------
IF ((PPTDAT(I,J,IHR).LT.900 .AND. DDATA(I,J).LE.ZER).OR.
& PPTDAT(I,J,IHR).LT.ZER) THEN
LTOP(I,J)=LBOT(I,J)
PTOP(I,J)=PBOT(I,J)
ENDIF
IF(LTOP(I,J).GT.LBOT(I,J))THEN
LTOP(I,J)=LBOT(I,J)
PTOP(I,J)=PBOT(I,J)
ENDIF
IF(HBM2(I,J).LT.0.90)THEN
LTOP(I,J)=LBOT(I,J)
PTOP(I,J)=PBOT(I,J)
ENDIF
C
C If the cloud is too shallow for convective precip, go to grid scale.
C
PSFCIJ=PD(I,J)+PT
DEPMIN=PSHNEW*PSFCIJ*1.E-5
DEPTH=PBOT(I,J)-PTOP(I,J)
if (i.eq.itstloc .and. j.eq.jtstloc .and. mype.eq.mtstpe)
& write(98,*) 'PTOP=',ptop(i,j), ' PBOT=',pbot(i,j),
& ' DEPTH=', DEPTH, ' DEPMIN=', DEPMIN
C
IF(DEPTH .LT. DEPMIN) THEN
IGSADJ(I,J) = .TRUE.
ADATA(I,J) = DDATA(I,J)
GOTO 600
ENDIF
C***********************************************************************
C************* IF CLOUD IS DEEP ENOUGH THEN ASSUME CONVECTION **********
C************* IS OBSERVED, MAKE CONVECTIVE-TYPE ADJUSTMENT ************
C***********************************************************************
C***********************************************************************
C
C ESTIMATE THE CHANGE IN EFI, BASIALLY MULTIPLYING CURRENT EFI
C BY FORECAST PREC/OBS PREC
C
C Don't worry - if we are here then DDATA > 0.
FACTOR=PREC(I,J)/DDATA(I,J)
C
C IF THERE WAS NO FORECAST PRECIP, LEAVE EFI ALONE
C
IF (PREC(I,J).LE.ZER) FACTOR=1.
EFINEW=CLDEFI(I,J)*FACTOR
EFI=CLDEFI(I,J)*FCB+EFINEW*FCC
IF (EFI.GT.1.0) EFI=1.0
IF (EFI.LT.0.2) EFI=0.2
IF (SM(I,J).LT.1.0.AND.DDATA(I,J).LT.CTHRS) EFI=1.0
CLDEFI(I,J)=EFI
C
C TAKE BACK ANY LATENT HEAT/PRECIP THAT WAS RELEASED PREVIOUSLY
C SINCE WE'LL BE ADJUSTING TO THE PROFILE THAT RELEASES HEAT
C THAT SUMS UP TO THE OBSERVED PRECIP
C
DO L=1,LM
IF (HTM(I,J,L).GT.0.5) T(I,J,L)=T(I,J,L)-TLAT(I,J,L)
ENDDO
C
C TAKE BACK THE PRECIP TOO
C
CUPREC(I,J)=CUPREC(I,J)-(PREC(I,J)-APREC(I,J))
CUPPT(I,J)= CUPPT(I,J)-(PREC(I,J)-APREC(I,J))
ACPREC(I,J)=ACPREC(I,J)-PREC(I,J)
PREC (I,J)=0.
C
LTPK=LTOP(I,J)
LBTK=LBOT(I,J)
CDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCD
CDCDCDCDCDCDC DEEP CONVECTION DCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCD
CDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCD
LB =LBTK
EFI =CLDEFI(I,J)
C--------------INITIALIZE VARIABLES IN THE CONVECTIVE COLUMN------------
C***
C*** ONE SHOULD NOTE THAT THE VALUES ASSIGNED TO THE ARRAY TREFK
C*** IN THE 410 LOOP ARE REALLY ONLY RELEVANT IN ANCHORING THE
C*** REFERENCE TEMPERATURE PROFILE AT LEVEL LB. WHEN BUILDING THE
C*** REFERENCE PROFILE FROM CLOUD BASE, THEN ASSIGNING THE
C*** AMBIENT TEMPERATURE TO TREFK IS ACCEPTABLE. HOWEVER, WHEN
C*** BUILDING THE REFERENCE PROFILE FROM SOME OTHER LEVEL (SUCH AS
C*** ONE LEVEL ABOVE THE GROUND), THEN TREFK SHOULD BE FILLED WITH
C*** THE TEMPERATURES IN TREF(I,J,L) WHICH ARE THE TEMPERATURES OF
C*** THE MOIST ADIABAT THROUGH CLOUD BASE. BY THE TIME THE LINE
C*** NUMBERED 450 HAS BEEN REACHED, TREFK ACTUALLY DOES HOLD THE
C*** REFERENCE TEMPERATURE PROFILE.
C***
DO 410 L=1,LM
IF (HTM(I,J,L).LT.0.5) GO TO 410
DIFT (L)=0.
DIFQ (L)=0.
TKL =T(I,J,L)
TK (L)=TKL
TREFK (L)=TKL
QKL =Q(I,J,L)
QK (L)=QKL
QREFK (L)=QKL
PKL =AETA(L)*PDSL(I,J)+PT
PK (L)=PKL
PSK (L)=PKL
APEKL =APE(I,J,L)
APEK (L)=APEKL
THERK (L)=TREF(I,J,L)*APEKL
410 CONTINUE
C--------------DEEP CONVECTION REFERENCE TEMPERATURE PROFILE------------
cdrun
LTP1=LTPK+1
LBM1=LB-1
PKB=PK(LB)
PKT=PK(LTPK)
C--------------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)
C
DO 420 L=LTPK,LBM1
IVI=LTPK+LBM1-L
IF(T(I,J,IVI+1).LT.TFRZ)GO TO 430
STABDL=STABD
TREFKX=((THERKY-THERKX)*STABDL
1 +TREFKX*APEKXX)/APEKXY
TREFK(IVI)=TREFKX
APEKXX=APEKXY
THERKX=THERKY
APEKXY=APEK(IVI-1)
THERKY=THERK(IVI-1)
L0=IVI
PK0=PK(L0)
420 CONTINUE
C--------------FREEZING LEVEL AT OR ABOVE THE CLOUD TOP-----------------
L0M1=L0-1
GO TO 445
C--------------TEMPERATURE REFERENCE PROFILE ABOVE FREEZING LEVEL-------
430 L0M1=L0-1
RDP0T=1./(PK0-PKT)
DTHEM=THERK(L0)-TREFK(L0)*APEK(L0)
CCCCCCCCCCCCCCCDIR$ SHORTLOOP
DO 440 L=LTPK,L0M1
TREFK(L)=(THERK(L)-(PK(L)-PKT)*DTHEM*RDP0T)/APEK(L)
440 CONTINUE
C-----------------------------------------------------------------------
C------------- ADJUST TEMP PROFILE TO MATCH OBSERVED PPT ---------------
C-----------------------------------------------------------------------
445 CONTINUE
PRECMAX=0.
C
DO L=LTPK,LB
PRECMAX=PDSL(I,J)*DETA(L)*(TREFK(L)-TK(L))*CPRLG+PRECMAX
ENDDO
if (i.eq.itstloc .and. j.eq.jtstloc .and. mype.eq.mtstpe)
& write(98,*)'PRECMAX=', PRECMAX
C
IF (PRECMAX.LE.0.) THEN
C SEND THIS TO THE GRID-SCALE ADJUSTMENT
C NOT ENOUGH POSITIVE AREA TO DO CONVECTIVE ADJUSTMENT
LTOP(I,J)=LBOT(I,J)-3
PTOP(I,J)=PBOT(I,J)-2.*PSHNEW
ADATA(I,J) = DDATA(I,J)
IGSADJ(I,J) = .TRUE.
GOTO 600
ENDIF
C
C IF THE OBS PRECIP IS GREATER THAN THE MAX POSSIBLE PRECIP
C (WHICH IS THE AMOUNT YOU'D GET BY GOING ALL THE WAY TO THE
C REF PROFILE) ONLY DO THE MAX POSSIBLE PRECIP. SET EFI TO EFIMN
C IN AN ATTEMPT TO TRY TO GET THE GRID-SCALE PRECIP TO START
C TAKING OVER
C
RATIO=DDATA(I,J)/PRECMAX
IF (RATIO.GT.1.) THEN
RATIO=1.
EFI=EFIMN
CLDEFI(I,J)=EFIMN
IGSADJ(I,J) = .TRUE.
ENDIF
C
DO L=LTPK,LB
TREFK(L)=TK(L)+RATIO*(TREFK(L)-TK(L))
ENDDO
C--------------DEEP CONVECTION REFERENCE HUMIDITY PROFILE---------------
C DEFINE DSPS
DSPBK=((EFI-EFIMN)*SLOPBS+DSPBSS)*SM(I,J)
1 +((EFI-EFIMN)*SLOPBL+DSPBSL)*(1.-SM(I,J))
DSP0K=((EFI-EFIMN)*SLOP0S+DSP0SS)*SM(I,J)
1 +((EFI-EFIMN)*SLOP0L+DSP0SL)*(1.-SM(I,J))
DSPTK=((EFI-EFIMN)*SLOPTS+DSPTSS)*SM(I,J)
1 +((EFI-EFIMN)*SLOPTL+DSPTSL)*(1.-SM(I,J))
cccccccccccccccccCDIR$ SHORTLOOP
450 CONTINUE
DEPTH=PFRZ*PSFCIJ*1.E-5
DEPWL=PKB-PK0
DO 460 L=LTPK,LB
C--------------SATURATION PRESSURE DIFFERENCE---------------------------
IF(DEPWL .GE. DEPTH) THEN
IF(L.LT.L0)THEN
DSP=((PK0-PK(L))*DSPTK+(PK(L)-PKT)*DSP0K)/(PK0-PKT)
ELSE
DSP=((PKB-PK(L))*DSP0K+(PK(L)-PK0)*DSPBK)/(PKB-PK0)
ENDIF
ELSE
DSP=DSP0K
IF(L.LT.L0)
1 DSP=((PK0-PK(L))*DSPTK+(PK(L)-PKT)*DSP0K)/(PK0-PKT)
ENDIF
C--------------HUMIDITY PROFILE-----------------------------------------
IF(PK(L).GT.PQM)THEN
PSK(L)=PK(L)+DSP
APESK(L)=(1.E5/PSK(L))**CAPA
THSK(L)=TREFK(L)*APEK(L)
QREFK(L)=PQ0/PSK(L)*EXP(A2*(THSK(L)-A3*APESK(L))
1 /(THSK(L)-A4*APESK(L)))
ELSE
QREFK(L)=Q(I,J,L)
ENDIF
460 CONTINUE
C--------------HEATING, MOISTENING, PRECIPITATION-----------------------
PRECK =0.
cccccccccccccccccccccccCDIR$ SHORTLOOP
DO 530 L=LTPK,LB
PRECK =DETA(L)*(TREFK(L)-TK(L))+PRECK
530 CONTINUE
C
C--------------UPDATE PRECIPITATION, TEMPERATURE & MOISTURE-------------
C
PREC (I,J)=PDSL(I,J)*PRECK*CPRLG+PREC (I,J)
CUPREC(I,J)=PDSL(I,J)*PRECK*CPRLG*ENDFCTR + CUPREC(I,J)
CUPPT(I,J)=PDSL(I,J)*PRECK*CPRLG*ENDFCTR + CUPPT(I,J)
ACPREC(I,J)=PDSL(I,J)*PRECK*CPRLG*ENDFCTR + ACPREC(I,J)
ADATA(I,J)=DDATA(I,J)-PDSL(I,J)*PRECK*CPRLG
APREC(I,J) = 0.
if (i.eq.itstloc .and. j.eq.jtstloc .and. mype.eq.mtstpe)
& write(98,*) 'After deep conv, PREC=',PREC(I,J),' CUPREC=',
& CUPREC(I,J),' ACPREC=',ACPREC(I,J),' ADATA=',ADATA(I,J)
C
DO 580 L=LTPK,LB
C Calculate the relative humidity before the T and Q update:
PETAL=PDSL(I,J)*AETA(L)+PT
QCKL=PQ0/PETAL
2 *EXP(A2*(T(I,J,L)-A3)/(T(I,J,L)-A4))
RELHUM=Q(I,J,L)/QCKL
C
OLDRH=RELHUM
OLDQ=Q(I,J,L)
C
T(I,J,L)=TREFK(L)
Q(I,J,L)=QREFK(L)
C
C Increase RH by factor PFACTOR, but keep it under 90%:
C
QCKL=PQ0/PETAL
2 *EXP(A2*(T(I,J,L)-A3)/(T(I,J,L)-A4))
C
RELHUM= AMIN1(0.80, RELHUM*PFACTOR(I,J))
Q(I,J,L)=RELHUM*QCKL
C
c Calculate the water vapor increment for
c 1) the entire column
c 2) sfc-700mb
c
DELQ = (Q(I,J,L)-OLDQ) * DETA(L)*PDSL(I,J)/G
VAPINC(I,J)=VAPINC(I,J)+DELQ
C
IF (PETAL.GE.70000.) THEN
VAPINC7(I,J)=VAPINC7(I,J)+DELQ
ENDIF
C tlat(l) in the following print is meaningless. Included so as to be
C consistent with the earlier prints (for Pobs < Pmod)
IF (MYPE.EQ.MTSTPE .AND. I.EQ.ITSTLOC .AND. J.EQ.JTSTLOC)
2 WRITE(98,129) L,PETAL,TLAT(I,J,L), OLDRH, PFACTOR(I,J),
3 RELHUM, OLDQ, Q(I,J,L)
580 CONTINUE
C
CDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCD
CDCDCDCDCDCDC END OF DEEP CONVECTION DCDCDCDCDCDCDCD
CDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCD
C
C-----------------------------------------------------------------------
600 CONTINUE
C-----------------------------------------------------------------------
C--------------GATHER GRID-SCALE PRECIP ADJUSTMENT POINTS---------------
C-----------------------------------------------------------------------
C
IF(.NOT.IGSADJ(I,J)) GO TO 900
C
C***********************************************************************
C
LMHK=LMH(I,J)
C
IF (APREC(I,J).GT.ZER.AND.ADATA(I,J).GT.ZER) THEN
C
C THIS IS THE ADJUSTMENT WE DO IF WE HAVE RAIN BOTH IN THE
C DATA AND IN THE MODEL, MULTIPLY THE LATENT HEAT
C AT EACH LEVEL BY THE FRACTION: DATA/MODEL RAINFALL
C WHILE NOT CHANGING DELTA Q
C THE Q CHANGE IS IMPLICIT, INCREASING Q
C AND THEN REMOVING IT VIA CONDENSATION.
C MATCH THE RH THAT THE PROFILE HAD PRIOR TO THIS ADJUSTMENT
C IF THE RATIO OF OBS PPT TO PREC(K) IS > 10, SEND IT TO THE
C PARABOLIC PROFILE PART OF THIS ROUTINE
C
C Near the top of the model, the cirrus clouds produce a not-insignificant
C amount of snow. Whether this is physically true is debatable (the ice
C crystals dropping from these cirri play a role, to be sure, but it should
C be more in the sense of a catalyst (seeder-feeder mechanism) than in terms
C of _amount_ of snow. Anyway, if we increase the TLAT and Q at these levels
C we could be in trouble. Let's only do the adjustment below 200mb.
C
ADJUST=ADATA(I,J)/APREC(I,J)
if (i.eq.itstloc .and. j.eq.jtstloc .and. mype.eq.mtstpe)
& write(98,*) 'adjust=', adjust
IF (ADJUST.LE.10.0) THEN
C-----------FIND THE PRE-MODIFIED RELATIVE HUMIDITY FOR THIS POINT------
DO 640 L=1,LMHK
IF (HTM(I,J,L).LT.0.5) GO TO 640
PETAL=PDSL(I,J)*AETA(L)+PT
IF (PETAL.LE.20000.) go to 640
QCKL=PQ0/PETAL
2 *EXP(A2*(T(I,J,L)-A3)/(T(I,J,L)-A4))
RELHUM=Q(I,J,L)/QCKL
OLDQ=Q(I,J,L)
OLDCWM=CWM(I,J,L)
C MODIFY THE TEMP CHANGE AND PRECIP
T(I,J,L)=T(I,J,L)+TLAT(I,J,L)*(ADJUST-1.)
CYL IF (T(I,J,L).GT.258.) THEN
IF (T(I,J,L).GT.258. .and. SR(I,J).LT.0.9) THEN
ELV=ELWV-DLDT*(T(I,J,L)-A3)
ELSE
ELV=ELIV
END IF
DELT=TLAT(I,J,L)*(ADJUST-1.)
C The following is the 'accum pcp' version of DELT to take care of the
C fractional time step at the end of each EDAS segment:
DELTACP=TLAT(I,J,L)*(ADJUST*ENDFCTR-1.)
PREC(I,J)=DELT*DETA(L)*PDSL(I,J)*CP/(ELV*ROW*G)
2 +PREC(I,J)
APREC(I,J)=DELTACP*DETA(L)*PDSL(I,J)*CP/(ELV*ROW*G)
2 +APREC(I,J)
ACPREC(I,J)=DELTACP*DETA(L)*PDSL(I,J)
2 *CP/(ELV*ROW*G) + ACPREC(I,J)
C SET THE RH TO BE THE SAME AS IT WAS BEFORE THE LATENT HEAT MODIFI.
QCKL=PQ0/PETAL * EXP(A2*(T(I,J,L)-A3)/(T(I,J,L)-A4))
RELHUM= RELHUM * PFACTOR(I,J)
Q(I,J,L)=RELHUM*QCKL
C
C Cloud adjustment for grid-scale precip: increase/decrease
C the cloud water mixing ratio proportionally (take care not to go below
C the minimum CWM to produce rain). If the precip is convective,
C not adjusting cloud.
C
C Cloud adjustment.
C First, calculate minimum cloud water for rain production. Note: we're
C only doing the adjustment in levels where TLAT > 0, i.e. where rain was
C produced.
C
C 2002/10/11: only do the cloud adjustment if ADJUST < 1 (i.e. to only reduce
C cloud, not to increase it. Remember that CWM is cloud water+cloud ice.
C If we just increase CWM proportionally, we could be increasing cloud ice
C by nearly 10 fold right underneath 200mb (10 fold is the cap).
C Much more prudent to not increase CWM at all, just make sure it's above
C the minimum required to form precip.
C
IF (TLAT(I,J,L).GT.0. .and. ADJUST.LT.1.) THEN
TTEMP=0.025*(T(I,J,L)-273.16)
WFIX=0.9814*EXP(0.01873*L)
WMIN=0.1E-3*EXP(TTEMP)*WFIX
c if(mype.eq.17) then
c print*,'in 640 loop'
c print*,'i,j,l=',i,j,l
c print*,'cwm(i,j,l),wmin=',cwm(i,j,l),wmin
c print*,'ttemp=',ttemp
c print*,'t(i,j,l)=',t(i,j,l)
c print*,'adjust=',adjust
c endif
CWM(I,J,L) = AMAX1(WMIN,CWM(I,J,L)*ADJUST)
ENDIF
c
c Calculate the water vapor and cloud water/ice increments for
c 1) the entire column
c 2) sfc-700mb
c
DELQ = (Q(I,J,L)-OLDQ) * DETA(L)*PDSL(I,J)/G
DELCWM = (CWM(I,J,L)-OLDCWM) * DETA(L)*PDSL(I,J)/G
c
VAPINC(I,J)=VAPINC(I,J)+DELQ
CLDINC(I,J)=CLDINC(I,J)+DELCWM
C
IF (PETAL.GE.70000.) THEN
VAPINC7(I,J)=VAPINC7(I,J)+DELQ
CLDINC7(I,J)=CLDINC7(I,J)+DELCWM
ENDIF
640 CONTINUE
ENDIF
ENDIF
C
IF ((APREC(I,J).LE.ZER.AND.ADATA(I,J).GT.ZER)
& .OR.ADJUST.GT.10.) THEN
C FINALLY, IF THE MODEL SAYS THAT NO RAIN OCCURED WHILE THE
C DATA SAYS THAT WE GOT SOMETHING,
C OR IF THE RATIO OF OBS PRECIP TO FORECASTED PRECIP IS > 10
C THE FOLLOWING OCCURS
C WE SPECIFY A PARABOLIC LATENT HEAT PROFILE
C AND COLLECT THE APPROPRIATE AMOUNT OF RAIN FROM THIS HEAT PROFILE
C WE ARE GOING TO INCREASE Q IN THE HEATED LAYERS SO THE RH WILL BE
C 80%. THIS SHOULD HELP TO GET SOME MODEL PRECIP THE NEXT TIMESTEP
C CLOUD BASE IS THE FIRST LEVEL ABOVE GROUND WHERE RH>80%
C CLOUD TOP IS THE FIRST LEVEL ABOVE CLOUD BASE WHERE RH<80%
C IF THIS CLOUD IS TOO SHALLOW, THEN SPECIFY A X MB CLOUD
C X MB ABOVE GROUND, DEPENDING UPON PPT RATE
C-----------------------------------------------------------------------
PSFC=PD(I,J)+PT
FIWL1=0.
CLIMIT=1.E-20
C
DO 650 L=1,LMHK
IF (HTM(I,J,L).LT.0.5) GO TO 650
TMT0=(T(I,J,L)-273.16)
TMT15=AMIN1(TMT0,-15.)
AI=0.008855
BI=1.
IF(TMT0.LT.-20.)THEN
AI=0.007225
BI=0.9674
ENDIF
QW=PQ0/(PDSL(I,J)*AETA(L)+PT)
1 *EXP(A2*(T(I,J,L)-A3)/(T(I,J,L)-A4))
QI=QW*(BI+AI*AMIN1(TMT0,0.))
QINT=QW*(1.-0.00032*TMT15*(TMT15+15.))
IF(TMT0.LE.-40.)QINT=QI
C-------------------ICE-WATER ID NUMBER IW------------------------------
IF(TMT0.LT.-15.)THEN
FI=Q(I,J,L)-0.75*QI
IF(FI.GT.0.0.OR.CWM(I,J,L).GT.CLIMIT) THEN
FIW=1.
ELSE
FIW=0.
ENDIF
ENDIF
IF(TMT0.LT.0.0.AND.TMT0.GE.-15.0)THEN
FIW=0.
IF(FIWL1.GT.0.0.AND.CWM(I,J,L).GT.CLIMIT)FIW=1.
ENDIF
IF(TMT0.GE.0.)THEN
FIW=0.
ENDIF
QC(L)=(1.-FIW)*QINT+FIW*QI
FIWL1=FIW
RELH(L)=Q(I,J,L)/QC(L)
650 CONTINUE
C
PBOT=0.
PTOP=0.
C
DO 660 L=LMHK,2,-1
IF (HTM(I,J,L).LT.0.5) GO TO 660
PETAL=PDSL(I,J)*AETA(L)+PT
IF (PETAL.LT.20000.) go to 660
IF(PBOT(I,J).EQ.0.0.AND.RELH(L).GE.0.80) THEN
PBOT(I,J)=PDSL(I,J)*AETA(L)+PT
PTOP(I,J)=PBOT(I,J)
ENDIF
IF(PBOT(I,J).GT.0.0.AND.PTOP(I,J).EQ.PBOT(I,J)
& .AND. RELH(L).LT.0.80) PTOP(I,J)=PDSL(I,J)*AETA(L)+PT
660 CONTINUE
C
IF (PBOT(I,J)-PTOP(I,J).LT.20000.0) THEN
C
C CLOUD SEARCH BASED UPON RH FAILED TO PRODUCE A SIGNIFICANT CLOUD
C SO SPECIFIY CLOUD FROM PRECIP RATE
C
C the following have been specified before loop 900:
C PTRES1=2.81E-03*SIXSIX
C PTRES2=3.75E-04*SIXSIX
C PTRES3=1.0E-03*SIXSIX
C
C THIS IS THE THRESHOLD VALUE FOR DETERMINING THE DEPTH OF CLOUD
C TO BE HEATED (OR WHETHER TO INCREASE RH IN Q ENHANCEMENT)
C
C CLOUD BASE IS 150 MB ABOVE SURFACE FOR ALL CLOUDS
C
PBOT(I,J)=PSFC-15000.
C
C CLOUD TOP IS AT 200 MB FOR INTENSE PRECIP
C
IF (ADATA(I,J).GE.PTRES1) PTOP(I,J)=20000.
C
C CLOUD DEPTH IS 450 MB FOR MODERATE PRECIP
C (HIGHEST CLOUD TOP ALLOWED IS AT 200 MB)
C
IF (ADATA(I,J).GE.PTRES2.AND.ADATA(I,J).LT.PTRES1) THEN
PTOP(I,J)=PBOT(I,J)-45000.
IF (PTOP(I,J).LT.20000.) PTOP(I,J)=20000.
END IF
C
C CLOUD DEPTH IS 300 MB FOR LIGHT PRECIP
C (HIGHEST CLOUD TOP ALLOWED IS AT 200 MB)
C
IF (ADATA(I,J).LT.PTRES2) THEN
PTOP(I,J)=PBOT(I,J)-30000.
IF (PTOP(I,J).LT.20000.) PTOP(I,J)=20000.
ENDIF
ENDIF
C
C FIND LAYERS JUST ABOVE PTOP AND PBOT
C
DO 670 L=1,LM
IF (HTM(I,J,L).LT.0.5) GO TO 670
PK(L)=PDSL(I,J)*AETA(L)+PT
IF (PK(L).LT.PBOT(I,J)) LCBOT=L
IF (PK(L).LT.PTOP(I,J)) LCTOP=L
670 CONTINUE
C
NUMLEV=LCBOT-LCTOP+1
PREC1=(ADATA(I,J)-APREC(I,J))/NUMLEV
DETACL=0.0
C
DO 680 L=LCTOP,LCBOT
DETACL=DETACL+DETA(L)
680 CONTINUE
C
C THIS VERSION HAS A PARABOLIC PROFILE OF PRECIP
C WHICH ALLOWS FOR A CHANGE IN LATENT HEATING WITH
C TEMPERATURE, ESPECIALLY NEAR THE FREEZING LEVEL
C
C I DEFINE THE ETATOP AND ETABOT TO BE THE INTERFACIAL
C LAYERS OF THE CLOUD OUTSIDE THE ACTUAL AETA(LCTOP)
C AND AETA(LCBOT)
C
ETATOP=AETA(LCTOP)-DETA(LCTOP)/2.0
ETABOT=AETA(LCBOT)+DETA(LCBOT)/2.0
DO 690 L=LCTOP,LCBOT
IF (T(I,J,L).GT.258.) THEN
ELV=ELWV-DLDT*(T(I,J,L)-A3)
ELSE
ELV=ELIV
END IF
OLDQ=Q(I,J,L)
ETBIG=AETA(L)*AETA(L)-(ETATOP+ETABOT)*AETA(L)+ETABOT*ETATOP
PRECL(I,J)=-6.0*PREC1*ETBIG/
& ((ETATOP-ETABOT)*(ETATOP-ETABOT))
DELT=PRECL(I,J)*G*ROW*ELV/(CP*DETA(L)*PDSL(I,J))
T(I,J,L)=DELT+T(I,J,L)
PREC(I,J)=DELT*DETA(L)*PDSL(I,J)*CP/(ELV*ROW*G)+PREC(I,J)
ACPREC(I,J)=DELT*DETA(L)*PDSL(I,J)*CP/
& (ELV*ROW*G)*ENDFCTR + ACPREC(I,J)
APREC(I,J)=DELT*DETA(L)*PDSL(I,J)*CP/
& (ELV*ROW*G)*ENDFCTR + APREC(I,J)
IF (I.EQ.ITSTLOC .AND. J.EQ.JTSTLOC .AND. MYPE.EQ.MTSTPE) THEN
yltmp = DELT*DETA(L)*PDSL(I,J)*CP/(ELV*ROW*G)
write(98,689) L, etbig,precl(i,j),yltmp,
& aprec(i,j), prec(i,j)
689 format(i2,3x,5(2x,e12.6))
endif
C
C KICK THE RH UP TO 80% IN THE HEATED LAYERS IF THE OBS PRECIP
C RATE IS MORE THAN 1.0 MM/HR
IF (ADATA(I,J).GE.PTRES3) THEN
QC(L)=HTM(I,J,L)*PQ0/(PDSL(I,J)*AETA(L)+PT)
2 *EXP(HTM(I,J,L)*A2*(T(I,J,L)-A3)/(T(I,J,L)-A4))
Q(I,J,L)=0.80*QC(L)
END IF
C
c
c Calculate the water vapor increment for
c 1) the entire column
c 2) sfc-700mb
c
DELQ = (Q(I,J,L)-OLDQ) * DETA(L)*PDSL(I,J)/G
VAPINC(I,J)=VAPINC(I,J)+DELQ
IF (PETAL.GE.70000.) THEN
VAPINC7(I,J)=VAPINC7(I,J)+DELQ
ENDIF
690 CONTINUE
CYL If below this newly added layer of cloud, the air is too dry, the
c added rain will quickly evaporate, leading to rapid cooling (could
c be 3 degs in one timestep in PRECPD), and the shock might lead to
c a blowup. So we need to moisten the air below the cloud layer too.
c We first tried setting it to 80%, minimum (80% or the original RH,
c whichever is greater), but that proved to be too much. So now
c I'm just moistening (when necessary) the three layers underneath
c the cloud base - 80% at LCBOT+1, 70% at LCBOT+2, and 60% at LCBOT+3.
c
c The cloud base shouldn't go below the ground surface.
c
LCBOTTM = MIN0(LMHK,LCBOT+3)
DO 700 L = LCBOT+1, LCBOTTM
IF (ADATA(I,J).GE.PTRES3) THEN
OLDQ=Q(I,J,L)
RHFCTR = 1. - 0.1*FLOAT(L-LCBOT+1)
QC(L)=HTM(I,J,L)*PQ0/(PDSL(I,J)*AETA(L)+PT)
2 *EXP(HTM(I,J,L)*A2*(T(I,J,L)-A3)/(T(I,J,L)-A4))
Q(I,J,L)=AMAX1(Q(I,J,L),RHFCTR*QC(L))
C
c Calculate the water vapor increment for
c 1) the entire column
c 2) sfc-700mb
c
DELQ = (Q(I,J,L)-OLDQ) * DETA(L)*PDSL(I,J)/G
VAPINC(I,J)=VAPINC(I,J)+DELQ
C
IF (PETAL.GE.70000.) THEN
VAPINC7(I,J)=VAPINC7(I,J)+DELQ
ENDIF
END IF
700 CONTINUE
C
C Cloud adjustment: for the layer of cloud we specified,
C we set cloud water mixing ratio to WMIN or the original CWM, which
C ever is greater.
C
DO 710 L = LCTOP, LCBOT
OLDCWM=CWM(I,J,L)
TTEMP=0.025*(T(I,J,L)-273.16)
WFIX=0.9814*EXP(0.01873*L)
WMIN=0.1E-3*EXP(TTEMP)*WFIX
c if(mype.eq.17) then
c print*,'in 710 loop'
c print*,'i,j,l=',i,j,l
c print*,'cwm(i,j,l),wmin=',cwm(i,j,l),wmin
c print*,'ttemp=',ttemp
c print*,'t(i,j,l)=',t(i,j,l)
c endif
CWM(I,J,L) = AMAX1(WMIN,CWM(I,J,L))
c
c Calculate the cloud water/ice increment for
c 1) the entire column
c 2) sfc-700mb
c
DELCWM = (CWM(I,J,L)-OLDCWM) * DETA(L)*PDSL(I,J)/G
CLDINC(I,J)=CLDINC(I,J)+DELCWM
C
IF (PETAL.GE.70000.) THEN
CLDINC7(I,J)=CLDINC7(I,J)+DELCWM
ENDIF
710 CONTINUE
C
ENDIF
C
C-----------------------------------------------------------------------
C***********************************************************************
C*******END OF HORIZONTAL LOOP FOR GRID-SCALE TYPE ADJUSTMENT **********
C***********************************************************************
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C--------------SAVE CLOUD TOP AND BOTTOM FOR RADIATION------------------
HTOP(I,J)=MIN(FLOAT(LTOP(I,J)),HTOP(I,J))
HBOT(I,J)=MAX(FLOAT(LBOT(I,J)),HBOT(I,J))
C***********************************************************************
C
900 CONTINUE
910 CONTINUE
C
C Zero out latent heat array to be ready for the next round of tracking/
C adjustments:
C
TLAT = 0.0
C
IF (MYPE.EQ.MTSTPE) THEN
WRITE(98,*) ' AT END OF ADJPPT, PREC=',PREC(ITSTLOC,JTSTLOC),
& ' APREC=',APREC(ITSTLOC,JTSTLOC),
& ' ACPREC=',ACPREC(ITSTLOC,JTSTLOC),
& ' CUPREC=', CUPREC(ITSTLOC,JTSTLOC)
WRITE(98,*)
ENDIF
C
RETURN
END