*IF DEF,A01_1A,OR,DEF,A01_1B,OR,DEF,A01_2A,OR,DEF,A01_2B,OR,DEF,A01_3A AWI3F402.7
C ******************************COPYRIGHT****************************** GTS2F400.9307
C (c) CROWN COPYRIGHT 1995, METEOROLOGICAL OFFICE, All Rights Reserved. GTS2F400.9308
C GTS2F400.9309
C Use, duplication or disclosure of this code is subject to the GTS2F400.9310
C restrictions as set forth in the contract. GTS2F400.9311
C GTS2F400.9312
C Meteorological Office GTS2F400.9313
C London Road GTS2F400.9314
C BRACKNELL GTS2F400.9315
C Berkshire UK GTS2F400.9316
C RG12 2SZ GTS2F400.9317
C GTS2F400.9318
C If no contract has been raised with this copy of the code, the use, GTS2F400.9319
C duplication or disclosure of it is strictly prohibited. Permission GTS2F400.9320
C to do so must first be obtained in writing from the Head of Numerical GTS2F400.9321
C Modelling at the above address. GTS2F400.9322
C ******************************COPYRIGHT****************************** GTS2F400.9323
C GTS2F400.9324
CLL Subroutine SOLANG ----------------------------------------------- SOLANG1A.3
CLL SOLANG1A.4
CLL Purpose : SOLANG1A.5
CLL Calculations of the earth's orbit described in the second page of SOLANG1A.6
CLL the "Calculation of incoming insolation" section of UMDP 23, i.e. SOLANG1A.7
CLL from the sin of the solar declination, the position of each point SOLANG1A.8
CLL and the time limits it calculates how much sunlight, if any, it SOLANG1A.9
CLL receives. SOLANG1A.10
CLL SOLANG1A.11
CLL Author: William Ingram SOLANG1A.12
CLL Model Modification history from model version 3.0: SOLANG1A.13
CLL version Date SOLANG1A.14
CLL SOLANG1A.15
CLL Programming standard : SOLANG1A.16
CLL Written to comply with 12/9/89 version of UMDP 4 (meteorological SOLANG1A.17
CLL standard). SOLANG1A.18
CLL Written in FORTRAN 77 with the addition of "!" comments and SOLANG1A.19
CLL underscores in variable names. SOLANG1A.20
CLL SOLANG1A.21
CLL Logical components covered : P233 SOLANG1A.22
CLL SOLANG1A.23
CLL Project task : SOLANG1A.24
CLL SOLANG1A.25
CLL External documentation: UMDP 23 SOLANG1A.26
CLL SOLANG1A.27
CLLEND ----------------------------------------------------------------- SOLANG1A.28
C*L SOLANG1A.29
SUBROUTINE SOLANG (SINDEC, T, DT, SINLAT, LONGIT, K, 5SOLANG1A.30
& LIT, COSZ) SOLANG1A.31
INTEGER!, INTENT(IN) :: SOLANG1A.32
& K ! Number of points SOLANG1A.33
REAL!, INTENT(IN) :: SOLANG1A.34
& SINDEC, ! Sin(solar declination) SOLANG1A.35
& T, DT, ! Start time (GMT) & timestep SOLANG1A.36
& SINLAT(K), ! sin(latitude) & longitude SOLANG1A.37
& LONGIT(K) ! of each point SOLANG1A.38
REAL!, INTENT(OUT) :: SOLANG1A.39
& LIT(K), ! Sunlit fraction of the timestep SOLANG1A.40
& COSZ(K) ! Mean cos(solar zenith angle) SOLANG1A.41
C ! during the sunlit fraction SOLANG1A.42
C* SOLANG1A.43
CL This routine has no dynamically allocated work areas. It calls the SOLANG1A.44
CL intrinsic functions SQRT, ACOS & SIN, but no user functions or SOLANG1A.45
CL subroutines. The only structure is a loop over all the points to be SOLANG1A.46
CL dealt with, with IF blocks nested inside to cover the various SOLANG1A.47
CL possibilities. SOLANG1A.48
INTEGER J ! Loop counter over points SOLANG1A.49
REAL TWOPI, ! 2*pi SOLANG1A.50
& S2R ! Seconds-to-radians converter SOLANG1A.51
REAL SINSIN, ! Products of the sines and of the cosines SOLANG1A.52
& COSCOS, ! of solar declination and of latitude. SOLANG1A.53
& HLD, ! Half-length of the day in radians (equal SOLANG1A.54
& ! to the hour-angle of sunset, and minus SOLANG1A.55
& COSHLD, ! the hour-angle of sunrise) & its cosine. SOLANG1A.56
& HAT, ! Local hour angle at the start time. SOLANG1A.57
& OMEGAB, ! Beginning and end of the timestep and SOLANG1A.58
& OMEGAE, ! of the period over which cosz is SOLANG1A.59
& OMEGA1, ! integrated, and sunset - all measured in SOLANG1A.60
& OMEGA2, ! radians after local sunrise, not from SOLANG1A.61
& OMEGAS, ! local noon as the true hour angle is. SOLANG1A.62
& DIFSIN, ! A difference-of-sines intermediate value SOLANG1A.63
& DIFTIM, ! and the corresponding time period SOLANG1A.64
& TRAD, DTRAD SOLANG1A.65
C ! These are the start-time and length of the timestep (T & DT) SOLANG1A.66
C ! converted to radians after midday GMT, or equivalently, hour SOLANG1A.67
C ! angle of the sun on the Greenwich meridian. SOLANG1A.68
*CALL C_PI 
SOLANG1A.69
PARAMETER ( TWOPI = 2. * PI, S2R = PI / 43200.) SOLANG1A.70
C SOLANG1A.71
TRAD = T * S2R - PI SOLANG1A.72
DTRAD = DT * S2R SOLANG1A.73
CDIR$ IVDEP SOLANG1A.74
! Fujitsu vectorization directive GRB0F405.549
!OCL NOVREC GRB0F405.550
DO 100 J = 1, K ! Loop over points SOLANG1A.75
HLD = 0. ! Logically unnecessary SOLANG1A.76
C statement without which the CRAY compiler will not vectorize this code SOLANG1A.77
SINSIN = SINDEC * SINLAT(J) SOLANG1A.78
COSCOS = SQRT( (1.-SINDEC**2) * (1.-SINLAT(J)**2) ) SOLANG1A.79
COSHLD = SINSIN / COSCOS SOLANG1A.80
IF (COSHLD.LT.-1.) THEN ! Perpetual night SOLANG1A.81
LIT(J) = 0. SOLANG1A.82
COSZ(J) = 0. SOLANG1A.83
ELSE SOLANG1A.84
HAT = LONGIT(J) + TRAD ! (3.2.2) SOLANG1A.85
IF (COSHLD.GT.1.) THEN ! Perpetual day - hour SOLANG1A.86
OMEGA1 = HAT ! angles for (3.2.3) are SOLANG1A.87
OMEGA2 = HAT + DTRAD ! start & end of timestep SOLANG1A.88
ELSE ! At this latitude some SOLANG1A.89
C points are sunlit, some not. Different ones need different treatment. SOLANG1A.90
HLD = ACOS(-COSHLD) ! (3.2.4) SOLANG1A.91
C The logic seems simplest if one takes all "times" - actually hour SOLANG1A.92
C angles - relative to sunrise (or sunset), but they must be kept in the SOLANG1A.93
C range 0 to 2pi for the tests on their orders to work. SOLANG1A.94
OMEGAB = HAT + HLD SOLANG1A.95
IF (OMEGAB.LT.0.) OMEGAB = OMEGAB + TWOPI SOLANG1A.96
IF (OMEGAB.GE.TWOPI) OMEGAB = OMEGAB - TWOPI SOLANG1A.97
IF (OMEGAB.GE.TWOPI) OMEGAB = OMEGAB - TWOPI SOLANG1A.98
C ! Line repeated - otherwise could have failure if SOLANG1A.99
C ! longitudes W are > pi rather than < 0. SOLANG1A.100
OMEGAE = OMEGAB + DTRAD SOLANG1A.101
IF (OMEGAE.GT.TWOPI) OMEGAE = OMEGAE - TWOPI SOLANG1A.102
OMEGAS = 2. * HLD SOLANG1A.103
C Now that the start-time, end-time and sunset are set in terms of hour SOLANG1A.104
C angle, can set the two hour-angles for (3.2.3). The simple cases are SOLANG1A.105
C start-to-end-of-timestep, start-to-sunset, sunrise-to-end and sunrise- SOLANG1A.106
C -to-sunset, but two other cases exist and need special treatment. SOLANG1A.107
IF (OMEGAB.LE.OMEGAS .OR. OMEGAB.LT.OMEGAE) THEN SOLANG1A.108
OMEGA1 = OMEGAB - HLD SOLANG1A.109
ELSE SOLANG1A.110
OMEGA1 = - HLD SOLANG1A.111
ENDIF SOLANG1A.112
IF (OMEGAE.LE.OMEGAS) THEN SOLANG1A.113
OMEGA2 = OMEGAE - HLD SOLANG1A.114
ELSE SOLANG1A.115
OMEGA2 = OMEGAS - HLD SOLANG1A.116
ENDIF SOLANG1A.117
IF (OMEGAE.GT.OMEGAB.AND.OMEGAB.GT.OMEGAS) OMEGA2=OMEGA1 SOLANG1A.118
C Put in an arbitrary marker for the case when the sun does not rise SOLANG1A.119
C during the timestep (though it is up elsewhere at this latitude). SOLANG1A.120
C (Cannot set COSZ & LIT within the ELSE ( COSHLD < 1 ) block SOLANG1A.121
C because 3.2.3 is done outside this block.) SOLANG1A.122
ENDIF ! This finishes the ELSE (perpetual day) block SOLANG1A.123
DIFSIN = SIN(OMEGA2) - SIN(OMEGA1) ! Begin (3.2.3) SOLANG1A.124
DIFTIM = OMEGA2 - OMEGA1 SOLANG1A.125
C Next, deal with the case where the sun sets and then rises again SOLANG1A.126
C within the timestep. There the integration has actually been done SOLANG1A.127
C backwards over the night, and the resulting negative DIFSIN and DIFTIM SOLANG1A.128
C must be combined with positive values representing the whole of the SOLANG1A.129
C timestep to get the right answer, which combines contributions from SOLANG1A.130
C the two separate daylit periods. A simple analytic expression for the SOLANG1A.131
C total sun throughout the day is used. (This could of course be used SOLANG1A.132
C alone at points where the sun rises and then sets within the timestep) SOLANG1A.133
IF (DIFTIM.LT.0.) THEN SOLANG1A.134
DIFSIN = DIFSIN + 2. * SQRT(1.-COSHLD**2) SOLANG1A.135
DIFTIM = DIFTIM + 2. * HLD SOLANG1A.136
ENDIF SOLANG1A.137
IF (DIFTIM.EQ.0.) THEN SOLANG1A.138
C Pick up the arbitrary marker for night points at a partly-lit latitude SOLANG1A.139
COSZ(J) = 0. SOLANG1A.140
LIT(J) = 0. SOLANG1A.141
ELSE SOLANG1A.142
COSZ(J) = DIFSIN*COSCOS/DIFTIM + SINSIN ! (3.2.3) SOLANG1A.143
LIT(J) = DIFTIM / DTRAD SOLANG1A.144
ENDIF SOLANG1A.145
ENDIF ! This finishes the ELSE (perpetual night) block SOLANG1A.146
100 CONTINUE SOLANG1A.147
RETURN SOLANG1A.148
END SOLANG1A.149
*ENDIF DEF,A01_1A,OR,DEF,A01_1B,OR,DEF,A01_2A,OR,DEF,A01_3A ADB1F400.383