SUBROUTINE ZENITH(TIMES, DAYI, HOUR) !>-------------------------------------------------------------------------------------------------- !> SUBROUTINE ZENITH !> !> SUBPROGRAM: ZENITH - COMPUTE THE SOLAR ZENITH ANGLE !> PROGRAMMER: BLACK !> ORG: W/NMC2 !> DATE: 93-10-28 !> !> ABSTRACT: !> ZENITH CALCULATES THE COSINE OF THE SOLAR ZENITH ANGLES AT EACH POINT FOR USE IN SWRAD !> !> PROGRAM HISTORY LOG: !> 93-10-28 ????? - ORIGINATOR !> 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 !> !> INPUT ARGUMENT LIST: !> TIMES - THE FORECAST TIME IN SECONDS !> !> OUTPUT ARGUMENT LIST: !> DAYI - THE DAY OF THE YEAR !> HOUR - THE HOUR OF THE DAY !> !> INPUT/OUTPUT ARGUMENT LIST: !> NONE !> !> USE MODULES: CTLBLK !> F77KINDS !> GLB_TABLE !> MAPPINGS !> MPPCOM !> PARMETA !> PARM_TBL !> PHYS !> TEMPCOM !> TOPO !> !> DRIVER : INIT !> INITS !> RADTN !> RDTEMP !> !> CALLS : ----- !>-------------------------------------------------------------------------------------------------- USE CTLBLK USE F77KINDS USE GLB_TABLE USE MAPPINGS USE MPPCOM USE PARMETA USE PARM_TBL USE PHYS USE TEMPCOM USE TOPO ! IMPLICIT NONE ! REAL (KIND=R4KIND), PARAMETER :: D00 = 0.0E0 REAL (KIND=R4KIND), PARAMETER :: D5 = 0.5E0 REAL (KIND=R4KIND), PARAMETER :: H15 = 15.E0 REAL (KIND=R4KIND), PARAMETER :: H24 = 24.E0 REAL (KIND=R4KIND), PARAMETER :: H3600 = 3600.E0 REAL (KIND=R4KIND), PARAMETER :: CSID2 = 6.570982E-2 REAL (KIND=R4KIND), PARAMETER :: CSID3 = 1.002738E0 REAL (KIND=R4KIND), PARAMETER :: CLON2 = 9.856474E-1 REAL (KIND=R4KIND), PARAMETER :: PI = 3.1415926 REAL (KIND=R4KIND), PARAMETER :: PI2 = 2. * PI REAL (KIND=R4KIND), PARAMETER :: PIH = 0.5 * PI REAL (KIND=R4KIND), PARAMETER :: DEG2RD = 1.745329E-2 REAL (KIND=R4KIND), PARAMETER :: OBLIQ = 23.440 * DEG2RD ! #include "sp.h" ! INTEGER(KIND=I4KIND), PARAMETER :: LDA = LM + 9 INTEGER(KIND=I4KIND), PARAMETER :: LA = 13 INTEGER(KIND=I4KIND), PARAMETER :: IMJM = IM * JM - JM / 2 INTEGER(KIND=I4KIND), PARAMETER :: JAM = 6 + 2 * (JM - 10) INTEGER(KIND=I4KIND), PARAMETER :: LSCRCH = 4 * LM + 1 + LA + 1 INTEGER(KIND=I4KIND), PARAMETER :: L1 = LA + 1 INTEGER(KIND=I4KIND), PARAMETER :: L2 = LA + LM + 1 INTEGER(KIND=I4KIND), PARAMETER :: L3 = LA + 2 * LM + 1 INTEGER(KIND=I4KIND), PARAMETER :: L4 = LA + 3 * LM + 1 ! REAL (KIND=R4KIND) , INTENT(IN) ::& & TIMES ! REAL (KIND=R4KIND) , INTENT(OUT) ::& & DAYI ! REAL (KIND=R4KIND) , INTENT(OUT) ::& & HOUR ! LOGICAL(KIND=L4KIND) ::& & LEAP ! INTEGER(KIND=I4KIND) ::& & IYR , INDX , I , J , KMNTH , KNT ! INTEGER(KIND=I4KIND) ::& & LL1 , LL2 ! REAL (KIND=R4KIND) ::& & CSID1 , CLON1 , DAY , SLON , DEC , RA , SIDTIM , HRANG , SINALT , & & HRLCL ! REAL (KIND=R4KIND), DIMENSION(50) ::& & CSIDX , CLONX ! INTEGER(KIND=I4KIND), DIMENSION(12) ::& & MONTH !------------------------ ! IMPLICIT NONE VARIABLES !------------------------ INTEGER(KIND=I4KIND) ::& & N , K , NSOLAR ! DATA MONTH /31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30,31/ SAVE MONTH !-------------------------------------------------------------------------------------------------- ! READ ANNUAL VALUES OF THE GREENWICH MEAN SIDEREAL TIME FOR 0000 UTC DECEMBER 31 AND OF THE ! GEOMETRIC MEAN CELESTIAL LONGITUDE OF THE SUN FOR 0000 TDT DECEMBER 31 !-------------------------------------------------------------------------------------------------- NSOLAR = 19 OPEN (NSOLAR, FILE='./SOLAR', FORM='FORMATTED', STATUS='OLD') ! LL2 = 0 DO N=1,10 LL1 = LL2 + 1 LL2 = LL1 + 4 READ(NSOLAR,25) (CSIDX(K), K=LL1,LL2) 25 FORMAT(5F9.6) END DO ! ll2 = 0 ! DO N=1,10 LL1 = LL2 + 1 LL2 = LL1 + 4 READ(NSOLAR,30) (CLONX(K), K=LL1,LL2) 30 FORMAT(5F11.6) END DO ! CLOSE(NSOLAR) ! IYR = IDAT(3) - 1900 INDX = IYR - 85 CSID1 = CSIDX(INDX) CLON1 = CLONX(INDX) DAY = D00 LEAP = .FALSE. ! IF (MOD(IDAT(3),4) == 0) THEN MONTH(2) = 29 LEAP = .TRUE. END IF ! IF (IDAT(1) > 1) THEN KMNTH = IDAT(1) - 1 DO KNT=1,KMNTH DAY = DAY + REAL(MONTH(KNT)) END DO END IF !----------------------------------------------------------------------------------- ! CALCULATE EXACT NUMBER OF DAYS FROM BEGINNING OF YEAR TO FORECAST TIME OF INTEREST !----------------------------------------------------------------------------------- DAY = DAY + REAL(IDAT(2) - 1) + (REAL(IHRST)+ TIMES/H3600) / H24 DAYI = REAL(INT(DAY) + 1) HOUR = (DAY - DAYI) * H24 + H24 !------------------------------------------------------------------------------------ ! FIND CELESTIAL LONGITUDE OF THE SUN THEN THE SOLAR DECLINATION AND RIGHT ASCENSION. !------------------------------------------------------------------------------------ SLON = (CLON1 + CLON2 * DAY) * DEG2RD ! IF (SLON > PI2) SLON = SLON - PI2 ! DEC = ASIN(SIN(SLON) * SIN(OBLIQ)) RA = ACOS(COS(SLON) / COS(DEC)) ! IF (SLON > PI) RA = PI2 - RA !------------------------------------------------------------------ ! FIND THE GREENWICH SIDEREAL TIME THEN THE LOCAL SOLAR HOUR ANGLE. !------------------------------------------------------------------ SIDTIM = CSID1 + CSID2 * DAYI + CSID3 * HOUR SIDTIM = SIDTIM * H15 * DEG2RD ! IF (SIDTIM < D00) SIDTIM = SIDTIM + PI2 IF (SIDTIM > PI2) SIDTIM = SIDTIM - PI2 ! HRANG = SIDTIM - RA ! DO 100 J=MYJS,MYJE DO 100 I=MYIS,MYIE HRLCL = HRANG - GLON(I,J) !-------------------------------------------------------------------------------------------------- ! THE ZENITH ANGLE IS THE COMPLEMENT OF THE ALTITUDE THUS THE COSINE OF THE ZENITH ANGLE EQUALS THE ! SINE OF THE ALTITUDE. !-------------------------------------------------------------------------------------------------- SINALT = SIN(DEC) * SIN(GLAT(I,J)) + COS(DEC) * COS(HRLCL) * COS(GLAT(I,J)) ! IF (SINALT < D00) SINALT = D00 ! CZEN(I,J) = SINALT ! 100 END DO !-------------------------------------------------------------------------------------------------- ! IF THE FORECAST IS IN A DIFFERENT YEAR THAN THE START TIME, RESET DAYI TO THE PROPER DAY OF THE ! NEW YEAR (IT MUST NOT BE RESET BEFORE THE SOLAR ZENITH ANGLE IS COMPUTED). !-------------------------------------------------------------------------------------------------- IF (DAYI > 365.) THEN IF (.NOT. LEAP) THEN DAYI = DAYI - 365. ELSE IF (LEAP .AND. DAYI > 366.) THEN DAYI = DAYI - 366. END IF END IF ! RETURN ! END SUBROUTINE ZENITH