SUBROUTINE CALWXT_BOURG(T,Q,PINT,LMHK,LM,PPT,ZINT,PTYPE) C C FILE: BOURG.f C WRITTEN: 06 JULY 1999, MICHAEL BALDWIN C REVISIONS: C 20 SEP 1999 - BALDWIN: MAKE MORE CONSISTENT C WITH BOURGOUIN (1992) C 24 MAR 2006 - MANIKIN: ADDED TO WRF SOUNDING POST C C ROUTINE TO COMPUTE PRECIPITATION TYPE USING A DECISION TREE C APPROACH THAT USES THE SO-CALLED "ENERGY METHOD" OF C BOURGOUIN OF AES (CANADA) 1992 C C *** NOTE: THIS SCHEME HAS SOME "IN-BETWEEN" SCENARIOS WHICH C IT ADDRESSES WITH RANDOM NUMBER GENERATION. AS A RESULT, C RESULTS MAY NOT BE REPRODUCABLE. *** C C LIST OF VARIABLES NEEDED C INPUT: C T,TD,P,PINT,LMH,LM,PPT C C T - Mid layer temp (K) C pressures in log P) C Q - Mid layer spec humidity C PINT - Interfacial pressure (Pa) C LMHK - Number of layers C LM - MAX Number of layers C PPT - Precip (m) C ZINT - Interfacial geopotential (m) C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C NOTE: VERTICAL ORDER OF ARRAYS MUST BE LAYER 1 = TOP C ---- . C . C . C LAYER LMH = BOTTOM C (JUST LIKE IN THE ETA MODEL) C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C INTERNAL: C C C OUTPUT: C PTYPE - INSTANTANEOUS WEATHER TYPE. C ACTS LIKE A 4 BIT BINARY C 1111 = RAIN/FREEZING RAIN/ICE PELLETS/SNOW C WHERE THE ONE'S DIGIT IS FOR SNOW C THE TWO'S DIGIT IS FOR ICE PELLETS C THE FOUR'S DIGIT IS FOR FREEZING RAIN C AND THE EIGHT'S DIGIT IS FOR RAIN C C C------------------------------------------------------------- C IN OTHER WORDS... C C PTYPE=1 SNOW C PTYPE=2 ICE PELLETS/MIX WITH ICE PELLETS C PTYPE=4 FREEZING RAIN/MIX WITH FREEZING RAIN C PTYPE=8 RAIN C------------------------------------------------------------- C C INITIALIZE WEATHER TYPE ARRAY TO ZERO (IE, OFF). C WE DO THIS SINCE WE WANT PTYPE TO REPRESENT THE C INSTANTANEOUS WEATHER TYPE ON RETURN. REAL T(LM),Q(LM),PINT(LM+1),ZINT(LM) INTEGER PTYPE PARAMETER(PTHRES=0.02,G=9.80665) C PTYPE=0 IF (PPT.LE.PTHRES) RETURN T1=RTC() PSFCK=PINT(LMHK+1) C C SKIP THIS POINT IF NO PRECIP THIS TIME STEP C C FIND THE DEPTH OF THE WARM LAYER BASED AT THE SURFACE C THIS WILL BE THE CUT OFF POINT BETWEEN COMPUTING C THE SURFACE BASED WARM AIR AND THE WARM AIR ALOFT C C C LOWEST LAYER T C tlmhk = T(LMHK) iwrml = lmhk + 1 IF (tlmhk.ge.273.15) THEN DO l = lmhk, 2, -1 IF (t(l).ge.273.15.and.t(l-1).lt.273.15.and. + iwrml.eq.lmhk+1) iwrml = l END DO END IF C C NOW FIND THE HIGHEST ABOVE FREEZING LEVEL C lhiwrm = lmhk + 1 DO l = lmhk, 1, -1 IF (t(l).ge.273.15) lhiwrm = l END DO C ENERGY VARIABLES C SURFW IS THE POSITIVE ENERGY BETWEEN THE GROUND C AND THE FIRST SUB-FREEZING LAYER ABOVE GROUND C AREANE IS THE NEGATIVE ENERGY BETWEEN THE GROUND C AND THE HIGHEST LAYER ABOVE GROUND C THAT IS ABOVE FREEZING C AREAPE IS THE POSITIVE ENERGY "ALOFT" C WHICH IS THE WARM ENERGY NOT BASED AT THE GROUND C (THE TOTAL WARM ENERGY = SURFW + AREAPE) C C PINTK1 IS THE PRESSURE AT THE BOTTOM OF THE LAYER C PINTK2 IS THE PRESSURE AT THE TOP OF THE LAYER C DZKL IS THE THICKNESS OF THE LAYER C IFRZL IS A FLAG THAT TELLS US IF WE HAVE HIT C A BELOW FREEZING LAYER C pintk1 = psfck ifrzl = 0 areane = 0.0 areape = 0.0 surfw = 0.0 DO 20 l = lmhk, 1, -1 IF (ifrzl.eq.0.and.t(l).le.273.15) ifrzl = 1 pintk2=pint(L) DZKL=ZINT(L)-ZINT(L+1) c dzkl = t(i,j,l) * (q(i,j,l)*0.608+1.0) * rog * c + alog(pintk1/pintk2) area1 = alog(t(l)/273.15) * g * dzkl IF (t(l).ge.273.15) THEN IF (l.lt.iwrml) areape = areape + area1 IF (l.ge.iwrml) surfw = surfw + area1 ELSE IF (l.gt.lhiwrm) areane = areane + abs(area1) END IF pintk1 = pintk2 20 CONTINUE C C DECISION TREE TIME C IF (areape.lt.2.0) THEN C VERY LITTLE OR NO POSITIVE ENERGY ALOFT, CHECK FOR C POSITIVE ENERGY JUST ABOVE THE SURFACE TO DETERMINE RAIN VS. SNOW IF (surfw.lt.5.6) THEN C NOT ENOUGH POSITIVE ENERGY JUST ABOVE THE SURFACE C SNOW = 1 PTYPE = 1 goto 800 ELSE IF (surfw.gt.13.2) THEN C ENOUGH POSITIVE ENERGY JUST ABOVE THE SURFACE C RAIN = 8 PTYPE = 8 goto 800 ELSE C TRANSITION ZONE, ASSUME EQUALLY LIKELY RAIN/SNOW C PICKING A RANDOM NUMBER, IF <=0.5 SNOW t2=rtc() ta=t2-t1 call srand(ta) r1 = rand() IF (r1.le.0.5) THEN C SNOW = 1 PTYPE = 1 goto 800 ELSE C RAIN = 8 PTYPE = 8 goto 800 END IF END IF END IF C IF (areape.ge.2.0) THEN C SOME POSITIVE ENERGY ALOFT, CHECK FOR ENOUGH NEGATIVE ENERGY C TO FREEZE AND MAKE ICE PELLETS TO DETERMINE IP VS. ZR IF (areane.gt.66.0+0.66*areape) THEN C ENOUGH NEGATIVE AREA TO MAKE IP, C NOW NEED TO CHECK IF THERE IS ENOUGH POSITIVE ENERGY C JUST ABOVE THE SURFACE TO MELT IP TO MAKE RAIN IF (surfw.lt.5.6) THEN C NOT ENOUGH ENERGY AT THE SURFACE TO MELT IP C ICE PELLETS = 2 PTYPE = 2 goto 800 ELSE IF (surfw.gt.13.2) THEN C ENOUGH ENERGY AT THE SURFACE TO MELT IP C RAIN = 8 PTYPE = 8 goto 800 ELSE C TRANSITION ZONE, ASSUME EQUALLY LIKELY IP/RAIN C PICKING A RANDOM NUMBER, IF <=0.5 IP t2=rtc() ta=t2-t1 call srand(ta) r1 = rand() IF (r1.le.0.5) THEN C ICE PELLETS = 2 PTYPE = 2 goto 800 ELSE C RAIN = 8 PTYPE = 8 goto 800 END IF END IF ELSE IF (areane.lt.46.0+0.66*areape) THEN C NOT ENOUGH NEGATIVE ENERGY TO REFREEZE, CHECK SURFACE TEMP C TO DETERMINE RAIN VS. ZR IF (tlmhk.lt.273.15) THEN C FREEZING RAIN = 4 PTYPE = 4 goto 800 ELSE C RAIN = 8 PTYPE = 8 goto 800 END IF ELSE C TRANSITION ZONE, ASSUME EQUALLY LIKELY IP/ZR C PICKING A RANDOM NUMBER, IF <=0.5 IP t2=rtc() ta=t2-t1 call srand(ta) r1 = rand() IF (r1.le.0.5) THEN C STILL NEED TO CHECK POSITIVE ENERGY C JUST ABOVE THE SURFACE TO MELT IP VS. RAIN IF (surfw.lt.5.6) THEN C ICE PELLETS = 2 PTYPE = 2 goto 800 ELSE IF (surfw.gt.13.2) THEN C RAIN = 8 PTYPE = 8 goto 800 ELSE C TRANSITION ZONE, ASSUME EQUALLY LIKELY IP/RAIN C PICKING A RANDOM NUMBER, IF <=0.5 IP t2=rtc() ta=t2-t1 call srand(ta) r2 = rand() IF (r2.le.0.5) THEN C ICE PELLETS = 2 PTYPE = 2 goto 800 ELSE C RAIN = 8 PTYPE = 8 goto 800 END IF END IF ELSE C NOT ENOUGH NEGATIVE ENERGY TO REFREEZE, CHECK SURFACE TEMP C TO DETERMINE RAIN VS. ZR IF (tlmhk.lt.273.15) THEN C FREEZING RAIN = 4 PTYPE = 4 goto 800 ELSE C RAIN = 8 PTYPE = 8 goto 800 END IF END IF END IF END IF 800 CONTINUE RETURN END C C