C**********************************************************************
C program UV Index
C**********************************************************************
C program will read fixed fields
C get date , cycle and fcst hour
C determine Solar Zenith Angles
C Read ozone field
C From LUT (ozone, sza) determine Erythemal Dose Rates
C Determine which GADS aerosol fields to use
C Read GEFS_aer AOD and SSA
C Call Erythemal_UV inputs:AOD,SSA,SZA output:Ratio
C Adjust field for Elevation
C Read SNOD (snow depth)
C Read ALBDO (albedo) files
C Determine where snow depth is greater than 0.02 (~1 inch) then use
C the albdo
C Adjust field for Albedo (Snow)
C Read UVBCLR and UVBCLD fields
C Determine 2nd 3 hr UVBCLR AND UVBCLD
C Determine 2nd 3 hr UV Transmissivity
C Adjust field for UV Trans (Clouds)
C
C OUTPUT FILES:
C fort.052 - UV Index clrsky field (GRIB)
C
C SUBPROGRAMS CALLED: (LIST ALL CALLED FROM ANYWHERE IN CODES)
C UNIQUE: - ROUTINES THAT ACCOMPANY SOURCE FOR COMPILE
C grib - grib UV Index field
C LIBRARY:
C W3LIB - link to /nwprod/w3lib
C
C ATTRIBUTES:
C LANGUAGE: FORTRAN
C MACHINE: Linux
C
C Changes: 2009-02-25 C.Long Changed IFCST = 24
C to IFCST = FCSTHR
C 2012-12-03 C.Long Changed open calls to read symbolic linked
C file names
C 2016-04-13 C.Long Changed to read hourly files
C 2018-03-01 C.Long Change to use SNOD and ALBDO
C 2020-02-13 C.Long Made code ready to use GEFS AOD and SSA
C 2021-07-09 Hai-Tien Lee
C add UVIUSC text and grib output
C debug - avoid negative ALBADJ=A0 when ALB=0
C 2022-01-24 Hai-Tien Lee
C calculate UVI with 0.25deg fields for both Input and Output
C$$$
C-----------------------------------------------------------
C
PARAMETER (IGRD=1440, JGRD=721, NXY=IGRD*JGRD)
PARAMETER (MXBIT=32,LENPDS=28,LENGDS=32)
PARAMETER (MXSIZE=30+LENPDS+LENGDS+NXY*(MXBIT+1)/8)
C
INTEGER YYYYMMDD, EDATE, YYYY, MM, DD, DOY, DPY, UTC, CYCL,
> FCSTHR, YY, XHR, YHR
C
REAL OZONE(igrd,jgrd), SZA(igrd,jgrd), GLAT,
> EDRLUT(18,27), ELEV(igrd,jgrd), UV(NXY),
> UVBCLD(igrd,jgrd), UVBCLR(igrd,jgrd), ALBEDO(igrd,jgrd),
> UVBCLDx(igrd,jgrd), UVBCLRx(igrd,jgrd), ALBEDOx(igrd,jgrd),
> UVBCLDy(igrd,jgrd), UVBCLRy(igrd,jgrd), ALBEDOy(igrd,jgrd),
> ELADJ(igrd,jgrd), EL(igrd,jgrd),
> ALBSNO(igrd,jgrd), SNOD(igrd,jgrd)
C
REAL ERY1(igrd,jgrd), ERY2(igrd,jgrd), ERY3(igrd,jgrd),
> ERY4(igrd,jgrd), ERY5(igrd,jgrd), ERYOUT(igrd,jgrd)
REAL OPT(igrd,jgrd), SSA(igrd,jgrd)
REAL OPTS(igrd,jgrd), SSAS(igrd,jgrd)
REAL OPTW(igrd,jgrd), SSAW(igrd,jgrd)
REAL MON(12), MONL(12)
REAL UVTRANS(igrd,jgrd), AERO(igrd,jgrd)
REAL UVIUCS(igrd,jgrd),UVC(NXY)
C
REAL PI, D2R, R2D, FZ, FC, SOLDEC
C
INTEGER IHR, IFCST, IDSCALE, IDMODEL, IDPARA, IDLEVEL, ITOT
C
CHARACTER*1 KBUF(MXSIZE)
C
DATA MON/0,31,59,90,120,151,181,212,243,273,304,334/
DATA MONL/0,31,60,91,121,152,182,213,244,275,305,335/
C
COMMON /COMGRIB/YY,MM,DD,IHR,IFCST,IDSCALE,IDMODEL,IDPARA,
& IDLEVEL,ITOT
C
C-----------------------------------------------------------------------
C
PI = 3.1415927
D2R = PI/180.0
R2D = 180.0/PI
C
C...Read Input
C
READ(5,*) YYYYMMDD
READ(5,*) CYCL
READ(5,*) FCSTHR
READ(5,*) XHR
READ(5,*) YHR
C
WRITE(6,910) YYYYMMDD, CYCL, FCSTHR, XHR, YHR
C
YYYY = INT(YYYYMMDD/10000)
YY = MOD(YYYY,100)
MM = INT((YYYYMMDD-YYYY*10000)/100)
DD = YYYYMMDD-YYYY*10000-MM*100
C
C...Read Fixed Fields
C
C...Read Erythemal(oz,sza) LUT
C...FILE='ery_lut_fix'
READ(10,*) EDRLUT
C write(6,911) ((edrlut(i,j),i=1,18),j=1,27)
WRITE(6,FMT='("Erythemal LUT read")')
C
C...Read Elevation
C
C...FILE='orograp_fix'
READ(11,907) ELEV
C
WRITE(6,FMT='("Elevation data read")')
C
C...Read Input Fields (1,1) is 90North, 0East
C
C...FILE='albdox_in'
READ(22,*) nx1,ny1
READ(22,*) ALBEDOx
WRITE(6,FMT='("Albedo_x Data Read",i4,i4)') nx1,ny1
C
C...FILE='uvclrx_in'
READ(23,*) nx2, ny2
READ(23,*) UVBCLRx
WRITE(6,FMT='("UVBCLR_x Data Read",i4,i4)') nx2,ny2
C
C...FILE='uvcldx_in'
READ(24,*) nx3, ny3
READ(24,*) UVBCLDx
WRITE(6,FMT='("UVBCLD_x Data Read",i4,i4)') nx3,ny3
C
C...FILE='uvozoy_in'
READ(30,*) nx0, ny0
READ(30,*) OZONE
WRITE(6,FMT='("Ozone Data Read",i4,i4)') nx0, ny0
C
C...FILE='snod_in'
READ(31,*) nx1,ny1
READ(31,*) SNOD
WRITE(6,FMT='("Snow Depth Read",i4,i4)') nx1,ny1
C
C...FILE='albdoy_in'
READ(32,*) nx1,ny1
READ(32,*) ALBEDOy
WRITE(6,FMT='("Albedo_y Read",i4,i4)') nx1,ny1
C
C...FILE='uvclry_in'
READ(33,*) nx2, ny2
READ(33,*) UVBCLRy
WRITE(6,FMT='("UVBCLR_y Data Read",i4,i4)') nx2,ny2
C
C...FILE='uvcldy_in'
READ(34,*) nx3, ny3
READ(34,*) UVBCLDy
WRITE(6,FMT='("UVBCLD_y Data Read",i4,i4)') nx3,ny3
C
C Read GEFS Optical Thickness and SSA
C
C...FILE='opt_in'
READ(36,*) nx3, ny3
READ(36,*) OPT
WRITE(6,FMT='("AOD Data Read",i4,i4)') nx3,ny3
C
C...FILE='ssa_in'
READ(37,*) nx3, ny3
READ(37,*) SSA
WRITE(6,FMT='("SSA Data Read",i4,i4)') nx3,ny3
C--------------------------------------------------------------
C
C...Create X hour file from Y hour files
C
C...h6 = 6*6h-5*5h
C...h5 = 5*5h-4*4h
C...h4 = 4*4h-3*3h
C...h3 = 3*3h-2*2h
C...h2 = 2*2h-1*1h
C...h1 = 1h
C
IF (XHR .GT. 0) THEN
DO 50 J = 1, JGRD
DO 49 I = 1, IGRD
IF (SNOD(I,J) .LT. 0.01 .OR. SNOD(I,J) .GT. 5.0) THEN
SNOD(I,J)=-1.0
ENDIF
C...ALBEDO averages are accumulated over the 6 hour forecast period
ALBx = ALBEDOx(I,J)
ALBy = ALBEDOy(I,J)
IF (ALBx .GE. 0.0) THEN
ALBEDO(I,J) = FLOAT(YHR)*ALBy - FLOAT(XHR)*ALBx
ELSE
ALBEDO(I,J) = ALBy
ENDIF
IF ((SNOD(I,J) .GT. 0.02 .AND. SNOD(I,J).LT.5) .and.
>(ALBEDO(I,J).GT.0.0 .AND. ALBEDO(I,J).LE.100)) THEN
ALBSNO(I,J) = ALBEDO(I,J)
ELSE
ALBSNO(I,J) = -1.0
ENDIF
C...DUVB and CDUVB averages are accumulated over the 6 hour forecast period
CLDX = UVBCLDX(I,J)
CLDY = UVBCLDY(I,J)
IF (CLDX .GE. 0.0) THEN
UVBCLD(I,J) = FLOAT(YHR)*CLDY - FLOAT(XHR)*CLDX
ELSE
UVBCLD(I,J) = CLDY
ENDIF
C
CLRX = UVBCLRX(I,J)
CLRY = UVBCLRY(I,J)
IF (CLRX .GE. 0.0) THEN
UVBCLR(I,J) = FLOAT(YHR)*CLRY - FLOAT(XHR)*CLRX
ELSE
UVBCLR(I,J) = CLRY
ENDIF
49 CONTINUE
50 CONTINUE
ELSE
DO 52 J = 1, JGRD
DO 51 I = 1, IGRD
IF (SNOD(I,J) .LT. 0.01 .OR. SNOD(I,J).GT.5.0) THEN
SNOD(I,J)=-1.0
ENDIF
C
ALBEDO(I,J) = ALBEDOy(I,J)
C
C...Check to see if snow depth and albedo values are reasonable
C
IF ((SNOD(I,J) .GT. 0.02 .AND. SNOD(I,J).LT.5) .AND.
> (ALBEDO(I,J).GT.0 .AND. ALBEDO(I,J).LE.100)) THEN
ALBSNO(I,J) = ALBEDO(I,J)
ELSE
ALBSNO(I,J) = -1.0
ENDIF
C
UVBCLD(I,J) = UVBCLDY(I,J)
UVBCLR(I,J) = UVBCLRY(I,J)
51 CONTINUE
52 CONTINUE
ENDIF
C...Determine DOY and Days Per Year
C
IF (MOD(YYYY,4).EQ.0) THEN
DOY = MONL(MM) + DD
DPY = 366
ELSE
DOY = MON(MM) + DD
DPY = 365
ENDIF
WRITE(6,919) YYYY, MM, DD, DOY
C
C...from the cycle time and forecast hour
C...determine the day and UTC time of the product
C
FC = FLOAT(CYCL)
FZ = FLOAT(FCSTHR)
UTC = MOD((CYCL+FCSTHR),24)
C
DOY = DOY + (CYCL+FCSTHR)/24
WRITE(6,920) DOY, CYCL, FCSTHR, UTC
C
C...determine Solar Declination
C
SOLDEC = 23.45*SIN(2.0*PI*(DOY-80.0)/DPY)
c. WRITE(6,921) SOLDEC
SOLDEC = D2R*SOLDEC
C
C...determine Solar Zenith Angle
C
C GLON = 360.0/igrd
C GLAT = 180.0/(jgrd-1)
GLON = 0.25
GLAT = 0.25
C
DO 10 IL = 1, jgrd
RLT = 90.0 - (IL-1)*GLAT
RLAT = D2R*RLT
DO 9 IJ = 1, igrd
RLON = (IJ-1)*GLON
HRR = (UTC-12.0)*15.0 + RLON
HOUR = D2R*HRR
COSZ = SIN(SOLDEC)*SIN(RLAT)+COS(SOLDEC)*COS(RLAT)*COS(HOUR)
Z = R2D*ACOS(COSZ)
c...only work with SZA smaller than 85 degrees
IF (Z.GE.0.0 .AND. Z.LE.85.0) THEN
SZA(IJ,IL) = Z
ELSE
SZA(IJ,IL) = -1.0
ENDIF
9 CONTINUE
c. WRITE(6,923)il,rlt, MAXVAL(SZA(1:igrd,il)),
c. > MINVAL(SZA(1:igrd,il), MASK = SZA(1:igrd,il) .GT. 0)
10 CONTINUE
C
C...DETERMIN EARTH-SUN DISTANCE RATIO
C
CALL SUNDIS(YYYYMMDD, ESRAT)
WRITE(6,908) ESRAT
C
C...DETERMINE SEA LEVEL, NO AEROSOL, NO SNOW, NO CLOUDS EDR
C
DO 20 LTT = 1, jgrd
DO 19 LNN = 1, igrd
IF (SZA(LNN,LTT).GT.0.0 .AND. SZA(LNN,LTT).LT.85.0) THEN
OZ = 1+(OZONE(LNN,LTT)-80.0)/20.0
SZ = 1+SZA(LNN,LTT)/5.0
CALL W3FT01(SZ, OZ, EDRLUT, EOUT, 18, 27, 0, 1)
ERY1(LNN,LTT) = EOUT*ESRAT
C
C if (lnn.eq.576) write(6,940)glat(ltt),ozone(lnn,ltt),
C > oz, sza(lnn,ltt), sz, eout
C
c. if(lnn.eq.576) write(6,945) glat(ltt),sza(lnn,ltt),
c. > ozone(lnn,ltt),albedo(lnn,ltt),uvbclr(lnn,ltt),uvbcld(lnn,ltt)
ENDIF
19 CONTINUE
c. WRITE(6,924)ltt,glat(ltt), MAXVAL(ERY1(1:igrd,ltt)),
c. > MINVAL(ERY1(1:igrd,ltt))
20 CONTINUE
C
C...ADJUST Arrays FOR ELEVATION
C...Array operations
C
EL = ELEV/1000.0
ELADJ = -0.0009 + 5.4457*EL - 0.0414*EL**2
ERY2 = ERY1*(1.0+ELADJ/100.0)
c. DO 61 LTT = 1, JGRD
c. WRITE(6,926)ltt,glat(ltt), MAXVAL(ERY2(1:igrd,ltt)),
c. > MINVAL(ERY2(1:igrd,ltt))
c. 61 CONTINUE
C Assign UVIUCS as ERY2 (astro, ozone, elevation adjusted UVI)
UVIUCS=ERY2
C
C...Adjust for AOD and SSA
C...Compute Erythemal Ratio With vs W/O Aerosols
C
DO 60 LTT = 1, jgrd
LTX = JGRD - LTT + 1
DO 59 LNN = 1, igrd
IF (SZA(LNN,LTT).GT.0.0 .AND. SZA(LNN,LTT).LT.85.0) THEN
CALL ERYTH(OPT(LNN,LTT),SSA(LNN,LTT),SZA(LNN,LTT),RAT)
ELSE
RAT = 0.0
ENDIF
C
AERO(LNN,LTT) = RAT*100.0
ERY3(LNN,LTT) = ERY2(LNN,LTT)*RAT
59 CONTINUE
C WRITE(6,925)ltt,glat(ltt), MAXVAL(ERY3(1:igrd,ltt)),
C > MINVAL(ERY3(1:igrd,ltt))
60 CONTINUE
C
C...ADJUST FOR ALBEDO (SNOW)
C
ERY4 = ERY3
C
A0 = -1.32278
A1 = 0.43339
A2 = 0.00053
A3 = 0.00002
C
C...Adjust for snow if snow dept is greater than 0.02m
C
DO 80 LTT = 1, JGRD
N = 0
NN = 0
DO 79 LNN = 1, IGRD
ALB = ALBSNO(LNN,LTT)
C
IF (SZA(LNN,LTT).GT.0 .AND. SZA(LNN,LTT).LT.85) THEN
C 2021709 Hai-Tien Lee debug - avoid negative ALBADJ=A0 when ALB=0
if (ALB .gt. 0.) then
ALBADJ = A3*ALB**3 + A2*ALB**2 + A1*ALB + A0
else
ALBADJ=0.
endif
ERY4(LNN,LTT) = ERY3(LNN,LTT)*(1.0 + ALBADJ/100.0)
UVIUCS(LNN,LTT) = ERY2(LNN,LTT)*(1.0 + ALBADJ/100.0)
N = N + 1
NN = NN + 1
ENDIF
79 CONTINUE
c. WRITE(6,927)ltt,glat(ltt), MAXVAL(ERY4(1:igrd,ltt)),
c. > MINVAL(ERY4(1:igrd,ltt)), N, NN
80 CONTINUE
C
C...Determine UV Trans
C
ERY5 = ERY4
c
DO 90 LTT = 1, JGRD
DO 89 LNN = 1, IGRD
C23456789012345678901234567893123456789412345678951234567896123456789712
IF(UVBCLD(LNN,LTT).GT.0.0.AND.UVBCLR(LNN,LTT).GT.0.0)THEN
UVT = UVBCLD(LNN,LTT)/UVBCLR(LNN,LTT)
IF (UVT .GT. 1.0) UVT = 1.0
ERY5(LNN,LTT) = ERY4(LNN,LTT)*UVT
ELSE
UVT = -0.01
ENDIF
UVTRANS(LNN,LTT) = UVT*100.0
89 CONTINUE
c. WRITE(6,928)j,glat(j), MAXVAL(ERY5(1:igrd,j)),
c. > MINVAL(ERY5(1:igrd,j))
90 CONTINUE
C
C...changed from ery1 to ery2 August 25, 2005
C
WRITE(51,909) YYYYMMDD, CYCL, FCSTHR
WRITE(51,902) ERY2
C
WRITE(52,909) YYYYMMDD, CYCL, FCSTHR
WRITE(52,902) ERY5
C
WRITE(53,909) YYYYMMDD, CYCL, FCSTHR
WRITE(53,902) UVTRANS
C
WRITE(54,909) YYYYMMDD, CYCL, FCSTHR
WRITE(54,902) SZA
C
WRITE(55,909) YYYYMMDD, CYCL, FCSTHR
WRITE(55,902) AERO
C
WRITE(56,909) YYYYMMDD, CYCL, FCSTHR
WRITE(56,902) ALBEDO
C
WRITE(57,909) YYYYMMDD, CYCL, FCSTHR
WRITE(57,902) ALBSNO
C
WRITE(58,909) YYYYMMDD, CYCL, FCSTHR
WRITE(58,902) OZONE
C
WRITE(59,909) YYYYMMDD, CYCL, FCSTHR
WRITE(59,902) SNOD
C add UVIUCS output to LUN=60
C WRITE(60,909) YYYYMMDD, CYCL, FCSTHR
C WRITE(60,902) UVIUCS
C Experimental output
C WRITE(41,909) YYYYMMDD, CYCL, FCSTHR
C WRITE(41,902) ERY1
C WRITE(42,909) YYYYMMDD, CYCL, FCSTHR
C WRITE(42,902) ERY3
C WRITE(43,909) YYYYMMDD, CYCL, FCSTHR
C WRITE(43,902) ERY4
C
C...GRIB PARAMETERS
C IYY, IMM, IDD ARE PROVIDED BELOW
IHR=12
IFCST=FCSTHR
IDSCALE=0
IDMODEL=2
IDPARA=206
IDLEVEL=1
C...MODEL ID= 2, UV INDEX
C...PARAMETER ID=206, UV INDEX (unitless)
C...TYPE OF LEVEL= 1, SURFACE
C...SCALING FACTOR
C
C...Grib UV field as a product (Convert unit to Grib standard: W/m**2)
C
C 20230620 Hai-Tien Lee
C Note that the UV Index is unitless. It is obtained by scaling
C the erythemal-weighted UV radiative flux by 25 mW/m^2.
C This program, however, encodes UVI at flux unit "mW/m^2".
UVMAX = MAXVAL(ERY5)
UVMIN = MINVAL(ERY5)
WRITE(6,913) UVMAX, UVMIN
do j = 1, jgrd
do i = 1, igrd
k = (j-1)*igrd + i
uv(k) = ery5(i,j)*1.e-3
enddo
enddo
C
WRITE(6,912) MXSIZE
CALL GRIB(uv, kbuf)
C
open(62,status='replace',form='unformatted',access='direct',
> recl=itot)
WRITE(62,rec=1) (KBUF(i),i=1,ITOT)
C
C Add UVIUSC grib output
do j = 1, jgrd
do i = 1, igrd
k = (j-1)*igrd + i
UVC(k) = UVIUCS(i,j)*1.e-3
enddo
enddo
C
WRITE(6,912) MXSIZE
CALL GRIB(UVC, KBUF)
C
OPEN(63, status='replace', form='unformatted',access='direct',
> recl=itot)
WRITE(63,rec=1) (KBUF(i),i=1,ITOT)
C...FORMAT CARDS
C
900 FORMAT(1x,I8)
901 FORMAT(18F8.2)
902 FORMAT(10F8.2)
903 FORMAT(10F8.4)
904 FORMAT(8F9.4)
905 FORMAT(I4,2I2)
906 FORMAT(48I1)
907 FORMAT(10F8.1)
908 FORMAT(1X,'Earth-Sun Distance Ratio =', f9.6)
909 FORMAT(I8,1X,I2,1X,I3)
910 FORMAT(' YYYYMMDD =', I8,/,' CYCL =', I3,/,' FCST HR =', I3,
> /,'1st HOUR OF 6=',I3,/,'2nd Hour of 6=',i3)
911 FORMAT(18f5.0)
912 FORMAT(1X,'MXSIZE=',I8)
913 FORMAT(1X,'ERY5 max =',F8.2,' ERY min =',F8.2)
919 FORMAT(1X,'YYYY =', I5,', MM =', I3,', DD =', I3,', DOY =', I4)
920 FORMAT(1X,'DOY =', I4, ' CYCL=', I4,' FCST HR=', I4,' UTC=', I4)
921 FORMAT(1X,'SOLDEC (DEG)=', F8.2)
922 FORMAT(1X,'GLON =', F8.5)
923 FORMAT(1X,'IL=', i4,' RLAT=',F8.2,' MAX SZA =', F8.3,
> ' MIN SZA =', F8.3)
924 FORMAT(1X,'IL=', i4,' RLAT=',F8.2,' MAX ERY1 =', F8.3,
> ' MIN ERY1 =', F8.3)
925 FORMAT(1X,'IL=', i4,' RLAT=',F8.2,' MAX ERY2 =', F8.3,
> ' MIN ERY2 =', F8.3)
926 FORMAT(1X,'IL=', i4,' RLAT=',F8.2,' MAX ERY3 =', F8.3,
> ' MIN ERY3 =', F8.3)
927 FORMAT(1X,'IL=', i4,' RLAT=',F8.2,' MAX ERY4 =', F8.3,
> ' MIN ERY4 =', F8.3, ' N=',I4,' NN =', I4)
928 FORMAT(1X,'IL=', i4,' RLAT=',F8.2,' MAX ERY5 =', F8.3,
> ' MIN ERY5 =', F8.3)
940 FORMAT(1X,'LAT=',F6.2,' OZONE=',F8.2,' OZ =', F6.2,' SZA =',
> F8.2,' SZ =', F6.2, ' EOUT=',F8.2)
945 FORMAT(1x,'LAT=', F7.2,' SZA=', F6.2,' OZONE=', F7.2,
> ' ALBEDO=', F6.2,' UV CLR=', F5.2,' UV CLD=', F5.2)
950 FORMAT(2I4,4F8.4)
C
END