MODULE CABLE_WEATHERGENERATOR IMPLICIT NONE ! Parameters INTEGER, PARAMETER,PRIVATE :: sp = 8 REAL(sp),PARAMETER,PRIVATE :: & Pi = 3.14159265 ,& ! Pi PiBy2 = 1.57079632 ,& ! Pi/2 SecDay = 86400.0 ,& ! Seconds/day SolarConst = 1370*SecDay/1e6 ,& ! Solar constant [MJ/m2/day] epsilon = 0.736 ,& SBoltz = 5.67e-8 ! Stefan-Boltzmann constant [W/m2/K4] ! Global variables TYPE WEATHER_GENERATOR_TYPE ! general INTEGER :: np INTEGER :: ndtime ! number of subdiurnal steps in 24 hr REAL :: delT ! subdiurnal timestep (in s) REAL(sp),DIMENSION(:),ALLOCATABLE :: LatDeg ! Daily input REAL(sp),DIMENSION(:),ALLOCATABLE :: & WindDay ,& TempMinDay ,& ! Daily minimum air temp [degC] TempMaxDay ,& ! Daily maximum air temp [degC] TempMinDayNext ,& ! Daily minimum air temp tomorrow [degC] TempMaxDayPrev ,& ! Daily maximum air temp yesterday [degC] SolarMJDay ! Total daily down-ward irradiation [MJ/m2/d] (swdown) ! Daily constants REAL(sp) :: DecRad ! Declination in radians REAL(sp),DIMENSION(:),ALLOCATABLE :: & WindDark ,& ! Wind [m/s] during night hrs WindLite ,& ! Wind [m/s] during daylight hrs SolarNorm ,& ! (daily solar)/(solar const): geometry LatRad ,& ! Latitude in radians DayLength ,& TimeSunsetPrev ,& TimeSunrise ,& TimeMaxTemp ,& TimeSunset ,& TempSunsetPrev ,& TempSunset ,& TempNightRate ,& TempNightRatePrev ,& TempRangeDay ,& TempRangeAft ! Additional input to Subdiurnal met calculations ! Daily Met REAL(sp),DIMENSION(:),ALLOCATABLE :: & PrecipDay ,& ! 24hr-total precipitation [m/day] (current day = 24h from 0900 on previous day) SnowDay ,& ! 24hr-total precipitation [m/day] (current day = 24h from 0900 on previous day) !CLNx PrecipDayNext ,& ! 24hr-total precipitation [m/day] (next day = 24 h from 0900 on current day) !CLNx VapPmbDay ,& ! 24hr-av water vapour pressure [mb] !CLNx VapPmb0900 ,& ! 0900 water vapour pressure [mb] !CLNx VapPmb1500 ,& ! 1500 water vapour pressure [mb] !CLNx VapPmb1500Prev ,& ! 1500 (prev day) water vapour pressure [mb] !CLNx VapPmb0900Next ,& ! 0900(next day) water vapour pressure [mb] PmbDay ! 24hr-av pressure [mb] ! Solar and Temperature Params ! Hourly Met outgoing REAL(sp),DIMENSION(:),ALLOCATABLE :: & PhiSd ,& ! downward solar irradiance [W/m2] PhiLd ,& ! down longwave irradiance [W/m2] Precip ,& ! precip [mm/h] Snow ,& ! precip [mm/h] Wind ,& ! wind [m/s] Temp ,& ! temp [degC] VapPmb ,& ! vapour pressure [mb] Pmb ,& ! pressure [mb] coszen ! cos(theta) END TYPE WEATHER_GENERATOR_TYPE CONTAINS !=============================================================================== SUBROUTINE WGEN_INIT( WG, np, latitude, dels ) IMPLICIT NONE TYPE(WEATHER_GENERATOR_TYPE) :: WG INTEGER, INTENT(IN) :: np REAL, INTENT(IN) :: latitude(np), dels WG%np = np WG%delT = dels WG%ndtime = NINT ( REAL(SecDay)/dels ) ALLOCATE ( WG%LatDeg (np) ) ALLOCATE ( WG%WindDay (np) ) ALLOCATE ( WG%TempMinDay (np) ) ! Daily minimum air temp [degC] ALLOCATE ( WG%TempMaxDay (np) ) ! Daily maximum air temp [degC] ALLOCATE ( WG%TempMinDayNext (np) ) ! Daily minimum air temp tomorrow [degC] ALLOCATE ( WG%TempMaxDayPrev (np) ) ! Daily maximum air temp yesterday [degC] ALLOCATE ( WG%SolarMJDay (np) ) ALLOCATE ( WG%WindDark (np) ) ! Wind [m/s] during night hrs ALLOCATE ( WG%WindLite (np) ) ! Wind [m/s] during daylight hrs ALLOCATE ( WG%SolarNorm (np) ) ! (daily solar)/(solar const): geometry ALLOCATE ( WG%LatRad (np) ) ! Latitude in radians ALLOCATE ( WG%DayLength (np) ) ALLOCATE ( WG%TimeSunsetPrev (np) ) ALLOCATE ( WG%TimeSunrise (np) ) ALLOCATE ( WG%TimeMaxTemp (np) ) ALLOCATE ( WG%TimeSunset (np) ) ALLOCATE ( WG%TempSunsetPrev (np) ) ALLOCATE ( WG%TempSunset (np) ) ALLOCATE ( WG%TempNightRate (np) ) ALLOCATE ( WG%TempNightRatePrev (np) ) ALLOCATE ( WG%TempRangeDay (np) ) ALLOCATE ( WG%TempRangeAft (np) ) ALLOCATE ( WG%PrecipDay (np) ) ! 24hr-total precipitation [m/day] (current day = 24h from 0900 on previous day) ALLOCATE ( WG%SnowDay (np) ) ! [m/d] ALLOCATE ( WG%PmbDay (np) ) ! 24hr-av pressure [mb] ALLOCATE ( WG%PhiSd (np) ) ! downward solar irradiance [W/m2] ALLOCATE ( WG%PhiLd (np) ) ! down longwave irradiance [W/m2] ALLOCATE ( WG%Precip (np) ) ! precip [mm/h] ALLOCATE ( WG%Snow (np) ) ! precip [mm/h] ALLOCATE ( WG%Wind (np) ) ! wind [m/s] ALLOCATE ( WG%Temp (np) ) ! temp [degC] ALLOCATE ( WG%VapPmb (np) ) ! vapour pressure [mb] ALLOCATE ( WG%Pmb (np) ) ! pressure [mb] ALLOCATE ( WG%coszen (np) ) ! cos(theta) WG%LatDeg(:) = latitude(:) END SUBROUTINE WGEN_INIT SUBROUTINE WGEN_DAILY_CONSTANTS( WG, np, YearDay ) !------------------------------------------------------------------------------- ! Routine for calculating solar and met parameters that are constant over the ! the course of a 24hr day. These have been split off from version 01 of ! SubDiurnalMet which no longer calculates a day's worth of subdiurnal met in ! one call. Subdiurnal met for ntime times is now calculated using ntime calls, ! using the day-constants calculated here. !------------------------------------------------------------------------------- IMPLICIT NONE TYPE(WEATHER_GENERATOR_TYPE) :: WG INTEGER, INTENT(IN) :: np, YearDay ! Local variables REAL(sp) :: YearRad REAL(sp),DIMENSION(np) :: LatDeg1 REAL(sp),DIMENSION(np) :: TanTan REAL(sp),DIMENSION(np) :: HDLRad REAL(sp),DIMENSION(np) :: TimeSunriseNext REAL(sp) :: RatioWindLiteDark ! ------------------------- ! Downward solar irradiance ! ------------------------- LatDeg1 = SIGN(MIN(ABS(WG%LatDeg),89.9), WG%LatDeg) ! avoid singularity at pole WG%LatRad = LatDeg1*Pi/180.0 ! latitude in radians YearRad = 2.0*Pi*(YearDay-1)/365.0 ! day of year in radians !YearRadPrev = 2.0*Pi*(YearDay-2)/365.0 ! previous day of year in radians ! (Ok for YD - 2 = -1) ! DecRad = Declination in radians (+23.5 deg on 22 June, -23.5 deg on 22 Dec): WG%DecRad = 0.006918 - 0.399912*COS(YearRad) + 0.070257*SIN(YearRad) & - 0.006758*COS(2.0*YearRad) + 0.000907*SIN(2.0*YearRad) & - 0.002697*COS(3.0*YearRad) + 0.001480*SIN(3.0*YearRad) ! Paltridge and Platt eq [3.7] ! Daylength: HDLRad = Half Day Length in radians (dawn:noon = noon:dusk): TanTan = -TAN(WG%LatRad)*TAN(WG%DecRad) WHERE (TanTan .LE. -1.0) HDLRad = Pi ! polar summer: sun never sets ELSEWHERE (TanTan .GE. 1.0) HDLRad = 0.0 ! polar winter: sun never rises ELSEWHERE HDLRad = ACOS(TanTan) ! Paltridge and Platt eq [3.21] END WHERE ! (HDLRad = their capital PI) WG%DayLength = 24.0*2.0*HDLRad / (2.0*Pi) ! Daylength (dawn:dusk) in hours ! Daily solar irradiance without atmosphere, normalised by solar constant, ! with both energy fluxes in MJ/m2/day, calculated from solar geometry: WG%SolarNorm = & ! Paltridge and Platt eq [3.22] (HDLRad*SIN(WG%LatRad)*SIN(WG%DecRad) + & COS(WG%LatRad)*COS(WG%DecRad)*SIN(HDLRad)) / Pi ! ---- ! Wind ! ---- RatioWindLiteDark = 3.0 ! (daytime wind) / (nighttime wind) WG%WindDark = WG%WindDay / & ( (24.0-WG%DayLength)/24.0 + RatioWindLiteDark*WG%DayLength/24.0 ) WG%WindLite = WG%WindDay / & ( (1.0/RatioWindLiteDark)*(24.0-WG%DayLength)/24.0 + WG%DayLength/24.0 ) ! ----------- ! Temperature ! ----------- ! These are parameters required for the calculation of temperature according to ! Cesaraccio et al 2001 for sunrise-to-sunrise-next-day. Because we are calculating temp for ! midnight-to-midnight, we need to calculate midnight-to-sunrise temps using data for the ! previous day (-24h), hence the extra parameters denoted by *'s below which are not ! mentioned per se in Cesaraccio. Cesaraccio symbology for incoming met data is included ! here as comments for completeness: ! Sym in Cesaraccio et al 2001 ! TempMinDay Tn ! TempMaxDay Tx ! TempMinDayNext Tp ! TempMaxDayPrev !$write(71, "( 1000e16.6)") WG%DecRad !$ write(72, "( 1000e16.6)") TAN(WG%DecRad) !$ write(73, "( 1000e16.6)") TAN(WG%LatRad) !$ write(74, "( 1000e16.6)") WG%DayLength WHERE ( WG%DayLength .LT. 0.01 ) ! Polar night WG%TimeSunrise = 9. ! Hn WG%TimeSunset = 16. ! Ho WG%TimeMaxTemp = 13. ! Hx ELSEWHERE ( WG%DayLength .GT. 23.0 ) ! Polar day WG%TimeSunrise = 1. ! Hn WG%TimeSunset = 23. ! Ho WG%TimeMaxTemp = 13. ! Hx ELSEWHERE WG%TimeSunrise = (ACOS(MIN(TAN(WG%LatRad)*TAN(WG%DecRad),0.9999)))*12./Pi ! Hn WG%TimeSunset = WG%TimeSunrise + WG%DayLength ! Ho WG%TimeMaxTemp = WG%TimeSunset - MIN(4.,( WG%DayLength * 0.4)) ! Hx END WHERE WG%TimeSunsetPrev = WG%TimeSunset - 24. ! * Ho-24h (a negative hour) TimeSunriseNext = WG%TimeSunrise + 24. ! Hp WG%TempSunset = WG%TempMaxDay - & (0.39 * (WG%TempMaxDay - WG%TempMinDayNext)) ! To WG%TempSunsetPrev = WG%TempMaxDayPrev - & (0.39 * (WG%TempMaxDayPrev - WG%TempMinDay)) ! * To-24h WG%TempRangeDay = WG%TempMaxDay - WG%TempMinDay ! alpha = Tx-Tn WG%TempRangeAft = WG%TempMaxDay - WG%TempSunset ! R = Tx-To WG%TempNightRate = (WG%TempMinDayNext - WG%TempSunset)/ & SQRT(TimeSunriseNext-WG%TimeSunset) ! b = (Tp-To)/sqrt(Hp-Ho) WG%TempNightRatePrev = (WG%TempMinDay - WG%TempSunsetPrev)/ & SQRT(WG%TimeSunrise-WG%TimeSunsetPrev) ! * b-24h = (Tn-(To-24h))/sqrt(Hn-(Ho-24h)) END SUBROUTINE WGEN_DAILY_CONSTANTS !============================================================================== SUBROUTINE WGEN_SUBDIURNAL_MET(WG, np, itime) !------------------------------------------------------------------------------- ! Routine for downscaling daily met data to ntime subdiurnal values per day, ! at the instants ((it/ntime)*2400, it=1,ntime). ! ! ALGORITHMS: ! * Downward solar irradiance: ! * Fit a daily irradiance time series based on solar geometry, normalised to ! observed daily total SolarMJDay. Reference: Paltridge and Platt (1976). ! * Downward longwave irradiance: ! * Computed from downscaled air temperature using Swinbank (1963) formula. ! * Precipitation: ! * Assume steady through 24 hours at average value (needs improvement?) ! * Wind: ! * Set RatioWindLiteDark = (daytime wind) / (nighttime wind), a predetermined ! value, typically 3. Then calculate wind by day and wind by night (both ! steady, but different) to make the average work out right. ! * Temperature: ! * Fit a sine wave from sunrise (TempMinDay) to peak at TempMaxDay, another ! from TempMaxDay to sunset, TempMinDay, and a square-root function from ! sunset to sunrise the next day (Cesaraccio et al 2001). ! * Water Vapour Pressure, Air Pressure: ! * Hold constant. Return rel humidity at TempMinDay as a diagnostic on obs. ! ! HISTORY: ! * 04-sep-2003 (MRR): 01 Written and tested in Program SubDiurnalWeather ! * 30-nov-2007 (PRB): 02 Vectorise with deferred shape arrays, adding np dimension. ! * 11-dec-2007 (PRB): 03 Switch to Cesaraccio et al 2001 algorithm for temperature ! * 28/02/2012 (VH) : 04 cahnge precip from uniform distribution to evenly distributed over ! the periods 0600:0700; 0700:0800; 1500:1600; 1800:1900 !------------------------------------------------------------------------------- IMPLICIT NONE TYPE(WEATHER_GENERATOR_TYPE) :: WG INTEGER, INTENT(IN) :: np, itime ! Local variables REAL(sp) :: TimeNoon , test1, test2, adjust_fac(np) REAL(sp) :: TimeRad REAL(sp) :: rntime ! Real version of ntime REAL(sp) :: ritime ! Real version of current time REAL(sp),DIMENSION(np):: PhiLd_Swinbank ! down longwave irradiance [W/m2] !------------------------------------------------------------------------------- ritime = REAL(itime) * WG%delT/3600. ! Convert the current time to real rntime = REAL(WG%ndtime) * WG%delT/3600. ! Convert ntime to real ! Instantaneous downward hemispheric solar irradiance PhiSd TimeNoon = ritime/rntime - 0.5 ! Time in day frac (-0.5 to 0.5, zero at noon) TimeRad = 2.0*Pi*TimeNoon ! Time in day frac (-Pi to Pi, zero at noon) WHERE (ritime >= WG%TimeSunrise .AND. ritime <= WG%TimeSunset& .AND.WG%SolarNorm.GT.1.e-3 ) ! Sun is up WG%PhiSd = MAX((WG%SolarMJDay/WG%SolarNorm) & ! PhiSd [MJ/m2/day] * ( SIN(WG%DecRad)*SIN(WG%LatRad) + COS(WG%DecRad)*COS(WG%LatRad) & * COS(TimeRad) ), 0.0) ! fix to avoid negative PhiSD vh 13/05/08 ! Paltridge and Platt eq [3.4] WG%coszen = ( SIN(WG%DecRad)*SIN(WG%LatRad) + COS(WG%DecRad)*COS(WG%LatRad)*COS(TimeRad) ) ELSEWHERE ! sun is down WG%PhiSd = 0.0 WG%coszen = 0.0 END WHERE WG%PhiSd = WG%PhiSd*1e6/SecDay ! Convert PhiSd: [MJ/m2/day] to [W/m2] ! ------------- ! Precipitation ! ------------- !Precip = PrecipDay*1000./rntime ! convert from m/d to mm/h !Precip = (PrecipDay*1000.*9./24. + PrecipDayNext*1000.*15./24.)/rntime !Precip = (PrecipDay*1000.*9./24. + PrecipDayNext*1000.*15./24.)/24. ! hourly precip [mm/h] IF ( ABS(ritime-REAL(INT(ritime))) .GT. 1e-7) THEN WRITE(*,*) "Only works for integer hourly timestep! tstep = ", ritime STOP "cable_weathergenerator.F90!" ENDIF IF ((ritime >= 15. .AND. ritime < 16.).OR.(ritime >= 18. .AND. ritime < 19.)) THEN WG%Precip = WG%PrecipDay *1000./2. WG%Snow = WG%SnowDay*1000./2. ELSE WG%Precip = 0. WG%Snow = 0. ENDIF ! ---- ! Wind ! ---- WHERE (ritime >= WG%TimeSunrise .AND. ritime <= WG%TimeSunset) ! Sun is up WG%Wind = WG%WindLite ELSEWHERE ! Sun is down WG%Wind = WG%WindDark END WHERE ! ----------- ! Temperature ! ----------- ! Calculate temperature according to Cesaraccio et al 2001, including midnight to ! sunrise period using previous days info, and ignoring the period from after the ! following midnight to sunrise the next day, normally calculated by Cesaraccio. WHERE ( ritime <= WG%TimeSunrise ) ! Midnight to sunrise WG%Temp = WG%TempSunsetPrev + WG%TempNightRatePrev * SQRT(ritime - WG%TimeSunsetPrev) ELSEWHERE (ritime > WG%TimeSunrise .AND. ritime <= WG%TimeMaxTemp) ! Sunrise to time of maximum temperature WG%Temp = WG%TempMinDay + & WG% TempRangeDay *SIN (((ritime-WG%TimeSunrise)/ & (WG%TimeMaxTemp-WG%TimeSunrise))*PiBy2) ELSEWHERE (ritime > WG%TimeMaxTemp .AND. ritime <= WG%TimeSunset) ! Time of maximum temperature to sunset WG%Temp = WG%TempSunset + & WG%TempRangeAft * SIN(PiBy2 + ((ritime-WG%TimeMaxTemp)/4.*PiBy2)) ELSEWHERE (ritime > WG%TimeSunset ) ! Sunset to midnight WG%Temp = WG%TempSunset + WG%TempNightRate * SQRT(ritime - WG%TimeSunset) END WHERE ! ----------------------------------- ! Water Vapour Pressure, Air Pressure ! ----------------------------------- !CLN VapPmb = VapPmbDay !CLN Pmb = PmbDay !CLN !CLN IF (ritime <= 9.) THEN !CLN ! before 9am !CLN VapPmb = VapPmb1500Prev + (VapPmb0900 - VapPmb1500Prev) * (9. + ritime)/18. !CLN ELSEIF (ritime > 9 .AND. ritime <= 15.) THEN !CLN ! between 9am and 15:00 !CLN VapPmb = VapPmb0900 + (VapPmb1500 - VapPmb0900) * (ritime - 9.)/(15.-9.) !CLN ELSEIF (ritime > 15.) THEN !CLN ! after 15:00 !CLN VapPmb = VapPmb1500 + (VapPmb0900Next - VapPmb1500) * (ritime - 15.)/18. !CLN END IF ! ---------------------------- ! Downward longwave irradiance ! ---------------------------- PhiLd_Swinbank = 335.97 * (((WG%Temp + 273.16) / 293.0)**6) ! [W/m2] (Swinbank 1963) ! ------------------------------- ! Alternate longwave formulation ! ---------------------------- WG%PhiLd = epsilon * SBoltz * (WG%Temp + 273.16)**4 ! [W/m2] (Brutsaert) WHERE (WG%PhiSd.GT.50.0) adjust_fac = ((1.17)**(WG%SolarNorm))/1.17 ELSEWHERE adjust_fac = 0.9 ENDWHERE WG%PhiLd = WG%PhiLd /adjust_fac * (1.0 + WG%PhiSd/8000.) ! adjustment (formulation from Gab Abramowitz) WHERE ((WG%PhiLd.GT.500.00).OR.(WG%PhiLd.LT.100.00)) WG%PhiLd = PhiLd_Swinbank ENDWHERE IF (ANY((WG%PhiLd.GT.750.00).OR.(WG%PhiLd.LT.100.00))) THEN !write(*,*) 'PhiLD out of range' ENDIF END SUBROUTINE WGEN_SUBDIURNAL_MET END MODULE CABLE_WEATHERGENERATOR