PROGRAM HYCOM_DIURNAL IMPLICIT NONE C C hycom_diurnal - Usage: C hycom_diurnal latitude wstart wend whrinc C outputs daily to hrly diurnal scale factor every whrinc hours: C from wind day wstart C to wind day wend C C Standard Output: day vs scale factor suitable for plotting C C this version for "serial" Unix systems. C C Alan J. Wallcraf, NRL, March 2012. C INTEGER IOS,IForce_File_Number,i,j,k,ipt INTEGER IARGC,itest,jtest INTEGER NARG,n2pad CHARACTER*240 CARG C INTEGER IDM,JDM,IDAY,IHOUR,IHR,IYEAR,ILAT REAL*8 alat,wstart,wstop,whrinc,TT,timeref, dum1,dum2 c double precision dtime_diurnl real*8 day365 c real*8, dimension (0:24,-91:91) :: & diurnl ! hourly vs latitude shortwave scale factor table C C READ ARGUMENTS. C NARG = IARGC() C IF (NARG.EQ.4) THEN CALL GETARG(1,CARG) READ(CARG,*) alat CALL GETARG(2,CARG) READ(CARG,*) wstart CALL GETARG(3,CARG) READ(CARG,*) wstop CALL GETARG(4,CARG) READ(CARG,*) whrinc ELSE WRITE(6,*) + 'Usage: hycom_diurnal latitude wstart wend whrinc' CALL EXIT(1) ENDIF C ILAT = nint(alat) write(6,'(a,i4)') '# time vs diurnal scale factor for lat = ',ilat C C TIME LOOP C TT=wstart call forday(TT,3, iyear,iday,ihour) day365 = mod(iday+364,365) call thermf_diurnal(diurnl, day365) dtime_diurnl = TT DO write(6,'(f12.4,f10.3)') TT,diurnl(ihour,ilat) C TT=TT+whrinc/24.d0 IF(TT.GT.wstop) EXIT call forday(TT,3, iyear,iday,ihour) if (TT-dtime_diurnl.gt.1.0) then day365 = mod(iday+364,365) call thermf_diurnal(diurnl, day365) endif ENDDO CALL EXIT(0) END subroutine forday(dtime,yrflag, iyear,iday,ihour) implicit none c real*8 dtime integer yrflag, iyear,iday,ihour c c --- converts model day to "calendar" date (year,ordinal-day,hour). c real*8 dtim1,day integer iyr,nleap c if (yrflag.eq.0) then c --- 360 days per model year, starting Jan 16 iyear = int((dtime+15.001d0)/360.d0) + 1 iday = mod( dtime+15.001d0 ,360.d0) + 1 ihour = (mod( dtime+15.001d0 ,360.d0) + 1.d0 - iday)*24.d0 c elseif (yrflag.eq.1) then c --- 366 days per model year, starting Jan 16 iyear = int((dtime+15.001d0)/366.d0) + 1 iday = mod( dtime+15.001d0 ,366.d0) + 1 ihour = (mod( dtime+15.001d0 ,366.d0) + 1.d0 - iday)*24.d0 c elseif (yrflag.eq.2) then c --- 366 days per model year, starting Jan 01 iyear = int((dtime+ 0.001d0)/366.d0) + 1 iday = mod( dtime+ 0.001d0 ,366.d0) + 1 ihour = (mod( dtime+ 0.001d0 ,366.d0) + 1.d0 - iday)*24.d0 c elseif (yrflag.eq.3) then c --- model day is calendar days since 01/01/1901 iyr = (dtime-1.d0)/365.25d0 nleap = iyr/4 dtim1 = 365.d0*iyr + nleap + 1.d0 day = dtime - dtim1 + 1.d0 if (dtim1.gt.dtime) then iyr = iyr - 1 elseif (day.ge.367.d0) then iyr = iyr + 1 elseif (day.ge.366.d0 .and. mod(iyr,4).ne.3) then iyr = iyr + 1 endif nleap = iyr/4 dtim1 = 365.d0*iyr + nleap + 1.d0 c iyear = 1901 + iyr iday = dtime - dtim1 + 1.001d0 ihour = (dtime - dtim1 + 1.001d0 - iday)*24.d0 c endif return end subroutine thermf_diurnal(diurnal, date) implicit none c real*8 diurnal(0:24,-91:91),date c c --- Calculate a table of latitude vs hourly scale factors c --- for the distribution of daily averaged solar radiation c --- the clear sky insolation formula of Lumb (1964) is used with c --- correction for the seasonally varying earth-sun distance. c --- According to reed (1977) the lumb formula gives values in close c --- agreement with the daily mean values of the seckel and beaudry c --- (1973) formulae derived from data in the smithsonian c --- meteorological tables --- (list, 1958). c c --- Lumb, F. E., 1964: The influence of cloud on hourly amounts of c --- total solar radiation at sea surface.Quart. J. Roy. Meteor. Soc. c --- 90, pp43-56. c c --- date = julian type real date - 1.0 (range 0. to 365.), c --- where 00z jan 1 = 0.0. c c --- Base on "QRLUMB" created 2-4-81 by Paul J Martin. NORDA Code 322. c real*8, parameter :: pi = 3.14159265 real*8, parameter :: raddeg = pi/180.0 c integer lat,ihr real*8 sindec,cosdec,alatrd,fd,ourang,sinalt,ri,qsum real*8 sum c c calc sin and cosin of the declination angle of the sun. call declin(date,sindec,cosdec) c c loop through latitudes do lat= -90,90 c calc latitude of site in radians. alatrd = lat*raddeg c c loop through hours sum = 0.0 do ihr= 0,23 c calc hour angle of the sun (the angular distance of the sun c from the site, measured to the west) in radians. fd = real(ihr)/24.0 ourang = (fd-0.5)*2.0*pi c calc sine of solar altitude. sinalt = sin(alatrd)*sindec+cos(alatrd)*cosdec*cos(ourang) c c calc clear-sky solar insolation from lumb formula. if (sinalt.le.0.0) then diurnal(ihr,lat) = 0.0 else ri=1.00002+.01671*cos(0.01720242*(date-2.1)) diurnal(ihr,lat) = 2793.0*ri*ri*sinalt*(.61+.20*sinalt) endif sum = sum + diurnal(ihr,lat) enddo !ihr if (sum.gt.0.0) then c rescale so that sum is 24.0 (daily average to diurnal factor) qsum = 24.0/sum do ihr= 0,23 diurnal(ihr,lat) = diurnal(ihr,lat)*qsum enddo !ihr endif diurnal(24,lat) = diurnal(0,lat) !copy for table lookup enddo !lat do ihr= 0,24 diurnal(ihr,-91) = diurnal(ihr,-90) !copy for table lookup diurnal(ihr, 91) = diurnal(ihr, 90) !copy for table lookup enddo !ihr return c contains subroutine declin(date,sindec,cosdec) implicit none c real*8 date,sindec,cosdec c c subroutine to calc the sin and cosin of the solar declination angle c as a function of the date. c date = julian type real date - 1.0 (range 0. to 365.), where 00z c jan 1 = 0.0. c sindec = returned sin of the declination angle. c cosdec = returned cosin of the declination angle. c formula is from fnoc pe model. c created 10-7-81. paul j martin. norda code 322. c real a c a=date sindec=.39785*sin(4.88578+.0172*a+.03342*sin(.0172*a)- & .001388*cos(.0172*a)+.000348*sin(.0344*a)-.000028*cos(.0344*a)) cosdec=sqrt(1.-sindec*sindec) return end subroutine declin end subroutine thermf_diurnal