SUBROUTINE CRHTAB(RHCL,IER)
C---------------------------------------------------------------------
C.. CLD-RH RELATIONS OBTAINED FROM MITCHELL-HAHN PROCEDURE, HERE READ
C CLD/RH TUNING TABLES FOR DAY 0,1,...,5 AND MERGE INTO 1 FILE..
C .............K.A.C. MAR 93
C USE ONLY ONE TABLE (DAY 1) FOR ALL FCST HRS....K.A.C. FEB 94
c... 4 cld types .... KAC FEB96
c... smooth out last bunch of bins of the tables...KAC AUG97
C OUTPUT:
C RHCL - TUNING TABLES FOR ALL FORECAST DAYS
C IER - =1 IF TABLES AVAILABLE.. =-1 IF NO TABLES
C--------------------------------------------------------------------
use machine
implicit none
!!
integer mcld,nseal,ida,nbin,nlon,nlat,lon,jl,nc,lat,kcl,nsl
integer ken,icrit,isat,ibs,nb,kt,it,l,m,j,k,itim,ier,icfq,i
integer iy,im,n,kd,nbdayi,id,ld
CRH1T PARAMETER (MCLD=3,NSEAL=2,IDA=6,
cmcl3 PARAMETER (MCLD=3,NSEAL=2,IDA=1,
PARAMETER (MCLD=4,NSEAL=2,IDA=1,
2 NBIN=100,NLON=2,NLAT=4)
real (kind=kind_io8) cstem,cfrac,clsat,rhsat,binscl
real (kind=kind_io8) RHFD(NBIN,NLON,NLAT,MCLD,NSEAL)
real (kind=kind_io8) RRHFD(NBIN,NLON,NLAT,MCLD,NSEAL)
real (kind=kind_io8) RTNFFD(NBIN,NLON,NLAT,MCLD,NSEAL)
real (kind=kind_io8) RRNFFD(NBIN,NLON,NLAT,MCLD,NSEAL)
real (kind=kind_io8) RHCF(NBIN,NLON,NLAT,MCLD,NSEAL)
real (kind=kind_io8) RTNFCF(NBIN,NLON,NLAT,MCLD,NSEAL)
INTEGER KPTS(NLON,NLAT,MCLD,NSEAL)
INTEGER KKPTS(NLON,NLAT,MCLD,NSEAL)
real (kind=kind_io8) RHC(NLON,NLAT,MCLD,NSEAL)
real (kind=kind_io8) RHCL (NBIN,NLON,NLAT,MCLD,NSEAL,IDA)
real (kind=kind_io8) RHCLA(NBIN,NLON,NLAT,MCLD,NSEAL)
INTEGER ICDAYS(15),IDATE(4)
real(kind=kind_io4) fhour
real(kind=kind_io4) RHFD4(NBIN,NLON,NLAT,MCLD,NSEAL)
real(kind=kind_io4) RTNFFD4(NBIN,NLON,NLAT,MCLD,NSEAL)
C........................... BEGIN HERE ..............
IER = 1
DO 8000 ITIM=1,IDA
ICFQ = 43 + ITIM-1
REWIND ICFQ
Cmcl3 NCLDS=1,2,3 (L,M,H)..JSL=1,2 (LAND,SEA)
cmcl4 MCLD=1,2,3,4 (BL,L,M,H)
BINSCL = 1./NBIN
DO 1000 M=1,NSEAL
DO 1000 L=1,MCLD
DO 1000 K=1,NLAT
DO 1000 J=1,NLON
DO 1000 I=1,NBIN
RRHFD(I,J,K,L,M) = 0.
RRNFFD(I,J,K,L,M) = 0.
1000 CONTINUE
DO 1001 M=1,NSEAL
DO 1001 L=1,MCLD
DO 1001 K=1,NLAT
DO 1001 J=1,NLON
KKPTS(J,K,L,M) = 0
1001 CONTINUE
C.... READ THE DATA OFF THE ROTATING FILE
READ (ICFQ,ERR=998,END=999) NBDAYI,ICDAYS
DO 53 LD=1,NBDAYI
ID = ICDAYS(LD) / 10000
IM = (ICDAYS(LD)-ID*10000) / 100
IY = ICDAYS(LD)-ID*10000-IM*100
53 CONTINUE
READ (ICFQ,ERR=998,END=999) FHOUR,IDATE
csela PRINT 3003,IDATE,FHOUR,ITIM
DO 1300 KD=1,NBDAYI
READ (ICFQ) RHFD4
RHFD=RHFD4
READ (ICFQ) RTNFFD4
RTNFFD=RTNFFD4
READ (ICFQ) KPTS
DO 1002 M=1,NSEAL
DO 1002 L=1,MCLD
DO 1002 K=1,NLAT
DO 1002 J=1,NLON
DO 1002 I=1,NBIN
RRHFD(I,J,K,L,M) = RRHFD(I,J,K,L,M) + RHFD(I,J,K,L,M)
RRNFFD(I,J,K,L,M) = RRNFFD(I,J,K,L,M)+RTNFFD(I,J,K,L,M)
1002 CONTINUE
DO 1003 M=1,NSEAL
DO 1003 L=1,MCLD
DO 1003 K=1,NLAT
DO 1003 J=1,NLON
KKPTS(J,K,L,M) = KKPTS(J,K,L,M) + KPTS(J,K,L,M)
1003 CONTINUE
1300 CONTINUE
C
DO 1004 M=1,NSEAL
DO 1004 L=1,MCLD
DO 1004 K=1,NLAT
DO 1004 J=1,NLON
DO 1004 I=1,NBIN
RHCF(I,J,K,L,M) = RRHFD(I,J,K,L,M)
RTNFCF(I,J,K,L,M) = RRNFFD(I,J,K,L,M)
1004 CONTINUE
DO 1005 M=1,NSEAL
DO 1005 L=1,MCLD
DO 1005 K=1,NLAT
DO 1005 J=1,NLON
KPTS(J,K,L,M) = KKPTS(J,K,L,M)
1005 CONTINUE
C..... COMPUTE THE CUMULATIVE FREQUENCY DISTRIBUTION..
DO 200 N=1,NSEAL
DO 200 K=1,MCLD
DO 200 L=1,NLAT
DO 200 J=1,NLON
DO 190 I=2,NBIN
RHCF(I,J,L,K,N) = RHCF(I-1,J,L,K,N) + RHCF(I,J,L,K,N)
RTNFCF(I,J,L,K,N)=RTNFCF(I-1,J,L,K,N) + RTNFCF(I,J,L,K,N)
190 CONTINUE
200 CONTINUE
DO 300 N=1,NSEAL
DO 300 L=1,NLAT
DO 300 J=1,NLON
DO 300 K=1,MCLD
DO 300 I=1,NBIN
IF (KPTS(J,L,K,N).GT.0) THEN
RHCF(I,J,L,K,N) = RHCF(I,J,L,K,N) / KPTS(J,L,K,N)
RTNFCF(I,J,L,K,N) = RTNFCF(I,J,L,K,N) / KPTS(J,L,K,N)
c... cause we mix calculations of rh retune with cray and ibm words
c the last value of rhcf is close to but ne 1.0,
c so we reset it in order that the 360 loop gives compleat tabl
c... rtnfcf caused couple of problems, seems to be ever so slightly
c gt 1.0
IF (I.EQ.NBIN) THEN
RHCF(I,J,L,K,N) = 1.0
END IF
IF (RTNFCF(I,J,L,K,N).GE.1.0) THEN
RTNFCF(I,J,L,K,N) = 1.0
END IF
ELSE
RHCF(I,J,L,K,N) = -0.1
RTNFCF(I,J,L,K,N) = -0.1
END IF
300 CONTINUE
DO 255 NSL=1,NSEAL
DO 255 KCL=1,MCLD
csela PRINT 264,KCL,NSL
csela PRINT 265,((KPTS(I,L,KCL,NSL),I=1,NLON),L=1,NLAT)
255 CONTINUE
DO 360 NSL=1,NSEAL
DO 360 K=1,MCLD
DO 360 L=1,NLAT
DO 360 J=1,NLON
IF (KPTS(J,L,K,NSL).LE.0) GO TO 317
DO 320 I=1,NBIN
ICRIT = I
IF (RHCF(I,J,L,K,NSL).GE.RTNFCF(1,J,L,K,NSL)) GO TO 350
320 CONTINUE
C... NO CRITICAL RH
317 ICRIT=-1
350 RHC(J,L,K,NSL) = ICRIT * BINSCL
360 CONTINUE
csela DO 1210 NSL=1,NSEAL
csela DO 1210 K=1,MCLD
csela PRINT 1221,K,NSL
csela DO 1210 L=1,NLAT
csela PRINT 211,(RHC(J,L,K,NSL),J=1,NLON)
csela 1210 CONTINUE
DO 450 NSL=1,NSEAL
DO 450 KEN=1,MCLD
DO 450 L=1,NLAT
DO 450 JL=1,NLON
DO 400 I=1,NBIN
RHCL(I,JL,L,KEN,NSL,ITIM) = -0.1
400 CONTINUE
450 CONTINUE
DO 751 NSL=1,NSEAL
DO 751 KEN=1,MCLD
DO 751 L=1,NLAT
DO 751 JL=1,NLON
IF (KPTS(JL,L,KEN,NSL).LE.0) GO TO 751
DO 753 I=1,NBIN
DO 755 J=1,NBIN
IF (RHCF(J,JL,L,KEN,NSL).GE.RTNFCF(I,JL,L,KEN,NSL)) THEN
RHCL(I,JL,L,KEN,NSL,ITIM) = J*BINSCL
GO TO 753
END IF
755 CONTINUE
753 CONTINUE
751 CONTINUE
DO 3000 LON=1,NLON
DO 3000 LAT=1,NLAT
DO 3000 NC=1,MCLD
DO 3000 NSL=1,NSEAL
ISAT = 0
DO 67 IT=1,NBIN
CFRAC = BINSCL * (IT-1)
IF (RHCL(IT,LON,LAT,NC,NSL,ITIM).LT.0.) THEN
PRINT 1941,IT,NSL,NC,LAT,LON
STOP
END IF
IF (IT.LT.NBIN.AND.RTNFCF(IT,LON,LAT,NC,NSL).GE.1.) THEN
IF (ISAT.LE.0) THEN
ISAT = IT
RHSAT = RHCL(IT,LON,LAT,NC,NSL,ITIM)
CLSAT = CFRAC
END IF
RHCL(IT,LON,LAT,NC,NSL,ITIM) =
1 RHSAT + (1.-RHSAT)*(CFRAC-CLSAT)/(1.-CLSAT)
END IF
IF (IT.EQ.NBIN) RHCL(IT,LON,LAT,NC,NSL,ITIM) = 1.
67 CONTINUE
3000 CONTINUE
c... smooth out the table as it reaches rh=1.0, via linear interpolation
c between location of rh ge .98 and the NBIN bin (where rh=1.0)
c... previously rh=1.0 occurred for many of the latter bins in the
c table, thereby giving a cloud value of less then 1.0 for rh=1.0
nb=nbin-2
DO 4000 LON=1,NLON
DO 4000 LAT=1,NLAT
DO 4000 NC=1,MCLD
DO 4000 NSL=1,NSEAL
do 4167 it=1,nbin
RHCLA(IT,LON,LAT,NC,NSL)=RHCL(IT,LON,LAT,NC,NSL,ITIM)
4167 continue
DO 4067 IT=1,nb
ibs=it
cfrac=binscl*ibs
IF (RHCL(IT,LON,LAT,NC,NSL,ITIM).ge..98) THEN
CC Print 4011,nsl,nc,lat,lon,ibs,nbin
CC 4011 format (1h ,'nsl,nc,lat,lon,ibs,nbin=',6i4)
do 4068 kt=ibs,nbin
cstem=binscl*kt
RHCLA(kt,LON,LAT,NC,NSL) =
1 rhcl(ibs,LON,LAT,NC,NSL,ITIM)+
2 (rhcl(nbin,LON,LAT,NC,NSL,ITIM)
3 -rhcl(ibs,LON,LAT,NC,NSL,ITIM))*
3 (cstem-cfrac)/(1.-cfrac)
c if (nc.eq.2.and.lat.eq.2.and.lo.eq.1.and.nsl.eq.2) then
c print 4012,kt,cstem,cfrac,rhcl(ibs,LON,LAT,NC,NSL,ITIM),
c 1 RHCLA(kt,LON,LAT,NC,NSL)
c 4012 format(1h ,'kt,cs,cf,rhibs,rhcla=',i5,4f12.8)
c end if
4068 continue
go to 4000
end if
4067 CONTINUE
4000 CONTINUE
c... restore table data to preferred location..
DO 4200 LON=1,NLON
DO 4200 LAT=1,NLAT
DO 4200 NC=1,MCLD
DO 4200 NSL=1,NSEAL
DO 4200 IT=1,NBIN
RHCL(IT,LON,LAT,NC,NSL,ITIM)= RHCLA(IT,LON,LAT,NC,NSL)
4200 CONTINUE
8000 CONTINUE
DO 8001 KEN=1,IDA
ICFQ = 42 + KEN
REWIND ICFQ
8001 CONTINUE
RETURN
998 PRINT 988,ITIM
IER = -1
RETURN
999 PRINT 989,ITIM
IER = -1
RETURN
11 FORMAT(1H ,' from crhtab DAYS ON FILE =',I5)
51 FORMAT(1H ,' from crhtab ARCHV DATA FROM DA,MO,YR=',3I4)
202 FORMAT(1H0,' MODEL RH ',' OBS RTCLD')
203 FORMAT(2F10.2)
210 FORMAT(1H ,' NO CRIT RH FOR LAT=',I3,' AND LON BAND=',I3,
1 ' LAND(=1) SEA=',I3)
211 FORMAT(1H ,15F6.2)
264 FORMAT(1H ,' NUMBER OF GG POINTS USED IN EACH AREA..BY LATITUDE',
1 '..FOR CLOUD TYPE=',I4,'SEALAND=',I2)
265 FORMAT(1H ,15I8)
988 FORMAT(1H ,'from crhtab ERROR READING TABLES FOR TIME=',I4)
989 FORMAT(1H ,'from crhtab E.O.F READING TABLES FOR TIME=',I4)
1221 FORMAT(1H0,' CRITICAL RH FOR LON,LAT ARRAYS FOR CLD TYPE=',I3,
1 ' LAND(=1) SEA=',I3)
1941 FORMAT(1H ,' NEG RHCL FOR IT,NSL,NC,LAT,LON=',5I4,'...STOPPP..')
3003 FORMAT(5X,'...LAST DATE/TIME AND CURRENT ITIM',/,10X,
1 4I15,F7.1,I6)
END