[4728] | 1 | SUBROUTINE COS_SZA(KSTART,KEND,KCOL,PGEMU,PGELAM,LDRADIATIONTIMESTEP,PMU0) |
---|
| 2 | |
---|
| 3 | !**** *COS_SZA* |
---|
| 4 | ! |
---|
| 5 | ! (C) Copyright 2015- ECMWF. |
---|
| 6 | ! |
---|
| 7 | ! This software is licensed under the terms of the Apache Licence Version 2.0 |
---|
| 8 | ! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. |
---|
| 9 | ! |
---|
| 10 | ! In applying this licence, ECMWF does not waive the privileges and immunities |
---|
| 11 | ! granted to it by virtue of its status as an intergovernmental organisation |
---|
| 12 | ! nor does it submit to any jurisdiction. |
---|
| 13 | ! |
---|
| 14 | ! Purpose. |
---|
| 15 | ! -------- |
---|
| 16 | ! Compute the cosine of the solar zenith angle. Note that this |
---|
| 17 | ! is needed for three different things: (1) as input to the |
---|
| 18 | ! radiation scheme in which it is used to compute the path |
---|
| 19 | ! length of the direct solar beam through the atmosphere, (2) |
---|
| 20 | ! every timestep to scale the solar fluxes by the incoming |
---|
| 21 | ! solar radiation at top-of-atmosphere, and (3) to compute the |
---|
| 22 | ! albedo of the ocean. For (1) we ideally want an average |
---|
| 23 | ! value for the duration of a radiation timestep while for (2) |
---|
| 24 | ! we want an average value for the duration of a model |
---|
| 25 | ! timestep. |
---|
| 26 | |
---|
| 27 | !** Interface. |
---|
| 28 | ! ---------- |
---|
| 29 | ! *CALL* *COS_SZA(...) |
---|
| 30 | |
---|
| 31 | ! Explicit arguments : |
---|
| 32 | ! ------------------ |
---|
| 33 | ! PGEMU - Sine of latitude |
---|
| 34 | ! PGELAM - Geographic longitude in radians |
---|
| 35 | ! LDRadiationTimestep - Is this for a radiation timestep? |
---|
| 36 | ! PMU0 - Output cosine of solar zenith angle |
---|
| 37 | |
---|
| 38 | ! Implicit arguments : |
---|
| 39 | ! -------------------- |
---|
| 40 | ! YRRIP%RWSOVR, RWSOVRM - Solar time for model/radiation timesteps |
---|
| 41 | ! RCODECM, RSIDECM - Sine/cosine of solar declination |
---|
| 42 | ! YRERAD%LAverageSZA - Average solar zenith angle in time interval? |
---|
| 43 | ! YRRIP%TSTEP - Model timestep in seconds |
---|
| 44 | ! YRERAD%NRADFR - Radiation frequency in timesteps |
---|
| 45 | |
---|
| 46 | ! Method. |
---|
| 47 | ! ------- |
---|
| 48 | ! Compute cosine of the solar zenith angle, mu0, from lat, lon |
---|
| 49 | ! and solar time using standard formula. If |
---|
| 50 | ! YRERAD%LAverageSZA=FALSE then this is done at a single time, |
---|
| 51 | ! which is assumed to be the mid-point of either the model or |
---|
| 52 | ! the radiation timestep. If YRERAD%LAverageSZA=TRUE then we |
---|
| 53 | ! compute the average over the model timestep exactly by first |
---|
| 54 | ! computing sunrise/sunset times. For radiation timesteps, mu0 |
---|
| 55 | ! is to be used to compute the path length of the direct solar |
---|
| 56 | ! beam through the atmosphere, and the fluxes are subsequently |
---|
| 57 | ! weighted by mu0. Therefore night-time values are not used, |
---|
| 58 | ! so we average mu0 only when the sun is above the horizon. |
---|
| 59 | |
---|
| 60 | ! Externals. |
---|
| 61 | ! ---------- |
---|
| 62 | |
---|
| 63 | ! Reference. |
---|
| 64 | ! ---------- |
---|
| 65 | ! ECMWF Research Department documentation of the IFS |
---|
| 66 | ! |
---|
| 67 | ! See also: Zhou, L., M. Zhang, Q. Bao, and Y. Liu (2015), On |
---|
| 68 | ! the incident solar radiation in CMIP5 |
---|
| 69 | ! models. Geophys. Res. Lett., 42, 1930–1935. doi: |
---|
| 70 | ! 10.1002/2015GL063239. |
---|
| 71 | |
---|
| 72 | ! Author. |
---|
| 73 | ! ------- |
---|
| 74 | ! Robin Hogan, ECMWF, May 2015 |
---|
| 75 | |
---|
| 76 | ! Modifications: |
---|
| 77 | ! -------------- |
---|
| 78 | |
---|
| 79 | USE PARKIND1 , ONLY : JPIM, JPRB |
---|
| 80 | USE YOMHOOK , ONLY : LHOOK, DR_HOOK |
---|
| 81 | USE YOMCST , ONLY : RPI, RDAY |
---|
| 82 | USE YOMRIP , ONLY : YRRIP |
---|
| 83 | USE YOERIP , ONLY : YRERIP |
---|
| 84 | USE YOERAD , ONLY : YRERAD |
---|
| 85 | USE YOMLUN , ONLY : NULOUT |
---|
| 86 | |
---|
| 87 | ! ------------------------------------------------------------------ |
---|
| 88 | |
---|
| 89 | IMPLICIT NONE |
---|
| 90 | |
---|
| 91 | INTEGER(KIND=JPIM),INTENT(IN) :: KSTART ! Start column to process |
---|
| 92 | INTEGER(KIND=JPIM),INTENT(IN) :: KEND ! Last column to process |
---|
| 93 | INTEGER(KIND=JPIM),INTENT(IN) :: KCOL ! Number of columns in arrays |
---|
| 94 | REAL(KIND=JPRB), INTENT(IN) :: PGEMU(KCOL) ! Sine of latitude |
---|
| 95 | REAL(KIND=JPRB), INTENT(IN) :: PGELAM(KCOL)! Longitude in radians |
---|
| 96 | LOGICAL, INTENT(IN) :: LDRADIATIONTIMESTEP ! Is this for a radiation timestep? |
---|
| 97 | REAL(KIND=JPRB), INTENT(OUT) :: PMU0(KCOL) ! Cosine of solar zenith angle |
---|
| 98 | |
---|
| 99 | ! Solar time at the start and end of the time interval |
---|
| 100 | REAL(KIND=JPRB) :: ZSOLARTIMESTART, ZSOLARTIMEEND |
---|
| 101 | |
---|
| 102 | ! The time of half a model/radiation timestep, in radians |
---|
| 103 | REAL(KIND=JPRB) :: ZHALFTIMESTEP |
---|
| 104 | |
---|
| 105 | ! For efficiency we precompute sin(solar declination)*sin(latitude) |
---|
| 106 | REAL(KIND=JPRB) :: ZSINDECSINLAT(KSTART:KEND) |
---|
| 107 | !...and cos(solar declination)*cos(latitude) |
---|
| 108 | REAL(KIND=JPRB) :: ZCOSDECCOSLAT(KSTART:KEND) |
---|
| 109 | ! ...and cosine of latitude |
---|
| 110 | REAL(KIND=JPRB) :: ZCOSLAT(KSTART:KEND) |
---|
| 111 | |
---|
| 112 | ! Tangent of solar declination |
---|
| 113 | REAL(KIND=JPRB) :: ZTANDEC |
---|
| 114 | |
---|
| 115 | ! Hour angles (=local solar time in radians plus pi) |
---|
| 116 | REAL(KIND=JPRB) :: ZHOURANGLESTART, ZHOURANGLEEND |
---|
| 117 | REAL(KIND=JPRB) :: ZHOURANGLESUNSET, ZCOSHOURANGLESUNSET |
---|
| 118 | |
---|
| 119 | INTEGER(KIND=JPIM) :: JCOL ! Column index |
---|
| 120 | |
---|
| 121 | REAL(KIND=JPRB) :: ZHOOK_HANDLE |
---|
| 122 | |
---|
| 123 | ! ------------------------------------------------------------------ |
---|
| 124 | IF (LHOOK) CALL DR_HOOK('COS_SZA',0,ZHOOK_HANDLE) |
---|
| 125 | |
---|
| 126 | ! An average solar zenith angle can only be computed if the solar time |
---|
| 127 | ! is centred on the time interval |
---|
| 128 | IF (YRERAD%LAVERAGESZA .AND. .NOT. YRERAD%LCENTREDTIMESZA) THEN |
---|
| 129 | WRITE(NULOUT,*) 'ERROR IN COS_SZA: LAverageSZA=TRUE but LCentredTimeSZA=FALSE' |
---|
| 130 | CALL ABOR1('COS_SZA: ABOR1 CALLED') |
---|
| 131 | ENDIF |
---|
| 132 | |
---|
| 133 | DO JCOL = KSTART,KEND |
---|
| 134 | ZCOSLAT(JCOL) = SQRT(1.0_JPRB - PGEMU(JCOL)**2) |
---|
| 135 | ENDDO |
---|
| 136 | |
---|
| 137 | IF (LDRADIATIONTIMESTEP) THEN |
---|
| 138 | ! Compute the effective cosine of solar zenith angle for a radiation |
---|
| 139 | ! timestep |
---|
| 140 | |
---|
| 141 | ! Precompute quantities that may be used more than once |
---|
| 142 | DO JCOL = KSTART,KEND |
---|
| 143 | ZSINDECSINLAT(JCOL) = YRERIP%RSIDECM * PGEMU(JCOL) |
---|
| 144 | ZCOSDECCOSLAT(JCOL) = YRERIP%RCODECM * ZCOSLAT(JCOL) |
---|
| 145 | ENDDO |
---|
| 146 | |
---|
| 147 | IF (.NOT. YRERAD%LAVERAGESZA) THEN |
---|
| 148 | ! Original method: compute the value at the centre of the |
---|
| 149 | ! radiation timestep (assuming that LCentredTimeSZA=TRUE - see |
---|
| 150 | ! updtim.F90) |
---|
| 151 | DO JCOL = KSTART,KEND |
---|
| 152 | ! It would be more efficient to do it like this... |
---|
| 153 | ! PMU0(JCOL)=MAX(0.0_JPRB, ZSinDecSinLat(JCOL) & |
---|
| 154 | ! & - ZCosDecCosLat(JCOL) * COS(YRERIP%RWSOVRM + PGELAM(JCOL))) |
---|
| 155 | ! ...but for bit reproducibility with previous cycle we do it |
---|
| 156 | ! like this: |
---|
| 157 | PMU0(JCOL) = MAX(0.0_JPRB, ZSINDECSINLAT(JCOL) & |
---|
| 158 | & - YRERIP%RCODECM*COS(YRERIP%RWSOVRM)*ZCOSLAT(JCOL)*COS(PGELAM(JCOL)) & |
---|
| 159 | & + YRERIP%RCODECM*SIN(YRERIP%RWSOVRM)*ZCOSLAT(JCOL)*SIN(PGELAM(JCOL))) |
---|
| 160 | ENDDO |
---|
| 161 | |
---|
| 162 | ELSE |
---|
| 163 | ! Compute the average MU0 for the period of the radiation |
---|
| 164 | ! timestep, excluding times when the sun is below the horizon |
---|
| 165 | |
---|
| 166 | ! First compute the sine and cosine of the times of the start and |
---|
| 167 | ! end of the radiation timestep |
---|
| 168 | ZHALFTIMESTEP = YRRIP%TSTEP * REAL(YRERAD%NRADFR) * RPI / RDAY |
---|
| 169 | ZSOLARTIMESTART = YRERIP%RWSOVRM - ZHALFTIMESTEP |
---|
| 170 | ZSOLARTIMEEND = YRERIP%RWSOVRM + ZHALFTIMESTEP |
---|
| 171 | |
---|
| 172 | ! Compute tangent of solar declination, with check in case someone |
---|
| 173 | ! simulates a planet completely tipped over |
---|
| 174 | ZTANDEC = YRERIP%RSIDECM / MAX(YRERIP%RCODECM, 1.0E-12) |
---|
| 175 | |
---|
| 176 | DO JCOL = KSTART,KEND |
---|
| 177 | ! Sunrise equation: cos(hour angle at sunset) = |
---|
| 178 | ! -tan(declination)*tan(latitude) |
---|
| 179 | ZCOSHOURANGLESUNSET = -ZTANDEC * PGEMU(JCOL) & |
---|
| 180 | & / MAX(ZCOSLAT(JCOL), 1.0E-12) |
---|
| 181 | IF (ZCOSHOURANGLESUNSET > 1.0) THEN |
---|
| 182 | ! Perpetual darkness |
---|
| 183 | PMU0(JCOL) = 0.0_JPRB |
---|
| 184 | ELSE |
---|
| 185 | ! Compute hour angle at start and end of time interval, |
---|
| 186 | ! ensuring that the hour angle of the centre of the time |
---|
| 187 | ! window is in the range -PI to +PI (equivalent to ensuring |
---|
| 188 | ! that local solar time = solar time + longitude is in the |
---|
| 189 | ! range 0 to 2PI) |
---|
| 190 | IF (YRERIP%RWSOVRM + PGELAM(JCOL) < 2.0_JPRB*RPI) THEN |
---|
| 191 | ZHOURANGLESTART = ZSOLARTIMESTART + PGELAM(JCOL) - RPI |
---|
| 192 | ZHOURANGLEEND = ZSOLARTIMEEND + PGELAM(JCOL) - RPI |
---|
| 193 | ELSE |
---|
| 194 | ZHOURANGLESTART = ZSOLARTIMESTART + PGELAM(JCOL) - 3.0_JPRB*RPI |
---|
| 195 | ZHOURANGLEEND = ZSOLARTIMEEND + PGELAM(JCOL) - 3.0_JPRB*RPI |
---|
| 196 | ENDIF |
---|
| 197 | |
---|
| 198 | IF (ZCOSHOURANGLESUNSET >= -1.0) THEN |
---|
| 199 | ! Not perpetual daylight or perpetual darkness, so we need |
---|
| 200 | ! to check for sunrise or sunset lying within the time |
---|
| 201 | ! interval |
---|
| 202 | ZHOURANGLESUNSET = ACOS(ZCOSHOURANGLESUNSET) |
---|
| 203 | IF (ZHOURANGLEEND <= -ZHOURANGLESUNSET & |
---|
| 204 | & .OR. ZHOURANGLESTART >= ZHOURANGLESUNSET) THEN |
---|
| 205 | ! The time interval is either completely before sunrise or |
---|
| 206 | ! completely after sunset |
---|
| 207 | PMU0(JCOL) = 0.0_JPRB |
---|
| 208 | CYCLE |
---|
| 209 | ENDIF |
---|
| 210 | |
---|
| 211 | ! Bound the start and end hour angles by sunrise and sunset |
---|
| 212 | ZHOURANGLESTART = MAX(-ZHOURANGLESUNSET, & |
---|
| 213 | & MIN(ZHOURANGLESTART, ZHOURANGLESUNSET)) |
---|
| 214 | ZHOURANGLEEND = MAX(-ZHOURANGLESUNSET, & |
---|
| 215 | & MIN(ZHOURANGLEEND, ZHOURANGLESUNSET)) |
---|
| 216 | ENDIF |
---|
| 217 | |
---|
| 218 | IF (ZHOURANGLEEND - ZHOURANGLESTART > 1.0E-8) THEN |
---|
| 219 | ! Compute average MU0 in the interval ZHourAngleStart to |
---|
| 220 | ! ZHourAngleEnd |
---|
| 221 | PMU0(JCOL) = ZSINDECSINLAT(JCOL) & |
---|
| 222 | & + (ZCOSDECCOSLAT(JCOL) & |
---|
| 223 | & * (SIN(ZHOURANGLEEND) - SIN(ZHOURANGLESTART))) & |
---|
| 224 | & / (ZHOURANGLEEND - ZHOURANGLESTART) |
---|
| 225 | |
---|
| 226 | ! Just in case... |
---|
| 227 | IF (PMU0(JCOL) < 0.0_JPRB) THEN |
---|
| 228 | PMU0(JCOL) = 0.0_JPRB |
---|
| 229 | ENDIF |
---|
| 230 | ELSE |
---|
| 231 | ! Too close to sunrise/sunset for a reliable calculation |
---|
| 232 | PMU0(JCOL) = 0.0_JPRB |
---|
| 233 | ENDIF |
---|
| 234 | |
---|
| 235 | ENDIF |
---|
| 236 | ENDDO |
---|
| 237 | ENDIF |
---|
| 238 | |
---|
| 239 | ELSE |
---|
| 240 | ! Compute the cosine of solar zenith angle for a model timestep |
---|
| 241 | |
---|
| 242 | ! Precompute quantities that may be used more than once |
---|
| 243 | DO JCOL = KSTART,KEND |
---|
| 244 | ZSINDECSINLAT(JCOL) = YRRIP%RSIDEC * PGEMU(JCOL) |
---|
| 245 | ZCOSDECCOSLAT(JCOL) = YRRIP%RCODEC * ZCOSLAT(JCOL) |
---|
| 246 | ENDDO |
---|
| 247 | |
---|
| 248 | IF (.NOT. YRERAD%LAVERAGESZA) THEN |
---|
| 249 | ! Original method: compute the value at the centre of the |
---|
| 250 | ! model timestep |
---|
| 251 | DO JCOL = KSTART,KEND |
---|
| 252 | ! It would be more efficient to do it like this... |
---|
| 253 | ! PMU0(JCOL) = MAX(0.0_JPRB, ZSinDecSinLat(JCOL) & |
---|
| 254 | ! & - ZCosDecCosLat(JCOL)*COS(YRRIP%RWSOVR + PGELAM(JCOL))) |
---|
| 255 | ! ...but for bit reproducibility with previous cycle we do it |
---|
| 256 | ! like this: |
---|
| 257 | PMU0(JCOL) = MAX(0.0_JPRB, ZSINDECSINLAT(JCOL) & |
---|
| 258 | & - YRRIP%RCODEC*COS(YRRIP%RWSOVR)*ZCOSLAT(JCOL)*COS(PGELAM(JCOL)) & |
---|
| 259 | & + YRRIP%RCODEC*SIN(YRRIP%RWSOVR)*ZCOSLAT(JCOL)*SIN(PGELAM(JCOL))) |
---|
| 260 | ENDDO |
---|
| 261 | |
---|
| 262 | ELSE |
---|
| 263 | ! Compute the average MU0 for the period of the model timestep |
---|
| 264 | |
---|
| 265 | ! First compute the sine and cosine of the times of the start and |
---|
| 266 | ! end of the model timestep |
---|
| 267 | ZHALFTIMESTEP = YRRIP%TSTEP * RPI / RDAY |
---|
| 268 | ZSOLARTIMESTART = YRRIP%RWSOVR - ZHALFTIMESTEP |
---|
| 269 | ZSOLARTIMEEND = YRRIP%RWSOVR + ZHALFTIMESTEP |
---|
| 270 | |
---|
| 271 | ! Compute tangent of solar declination, with check in case someone |
---|
| 272 | ! simulates a planet completely tipped over |
---|
| 273 | ZTANDEC = YRRIP%RSIDEC / MAX(YRRIP%RCODEC, 1.0E-12) |
---|
| 274 | |
---|
| 275 | DO JCOL = KSTART,KEND |
---|
| 276 | ! Sunrise equation: cos(hour angle at sunset) = |
---|
| 277 | ! -tan(declination)*tan(latitude) |
---|
| 278 | ZCOSHOURANGLESUNSET = -ZTANDEC * PGEMU(JCOL) & |
---|
| 279 | & / MAX(ZCOSLAT(JCOL), 1.0E-12) |
---|
| 280 | IF (ZCOSHOURANGLESUNSET > 1.0) THEN |
---|
| 281 | ! Perpetual darkness |
---|
| 282 | PMU0(JCOL) = 0.0_JPRB |
---|
| 283 | ELSE |
---|
| 284 | ! Compute hour angle at start and end of time interval, |
---|
| 285 | ! ensuring that the hour angle of the centre of the time |
---|
| 286 | ! window is in the range -PI to +PI (equivalent to ensuring |
---|
| 287 | ! that local solar time = solar time + longitude is in the |
---|
| 288 | ! range 0 to 2PI) |
---|
| 289 | IF (YRRIP%RWSOVR + PGELAM(JCOL) < 2.0_JPRB*RPI) THEN |
---|
| 290 | ZHOURANGLESTART = ZSOLARTIMESTART + PGELAM(JCOL) - RPI |
---|
| 291 | ZHOURANGLEEND = ZSOLARTIMEEND + PGELAM(JCOL) - RPI |
---|
| 292 | ELSE |
---|
| 293 | ZHOURANGLESTART = ZSOLARTIMESTART + PGELAM(JCOL) - 3.0_JPRB*RPI |
---|
| 294 | ZHOURANGLEEND = ZSOLARTIMEEND + PGELAM(JCOL) - 3.0_JPRB*RPI |
---|
| 295 | ENDIF |
---|
| 296 | |
---|
| 297 | IF (ZCOSHOURANGLESUNSET >= -1.0) THEN |
---|
| 298 | ! Not perpetual daylight or perpetual darkness, so we need |
---|
| 299 | ! to check for sunrise or sunset lying within the time |
---|
| 300 | ! interval |
---|
| 301 | ZHOURANGLESUNSET = ACOS(ZCOSHOURANGLESUNSET) |
---|
| 302 | IF (ZHOURANGLEEND <= -ZHOURANGLESUNSET & |
---|
| 303 | & .OR. ZHOURANGLESTART >= ZHOURANGLESUNSET) THEN |
---|
| 304 | ! The time interval is either completely before sunrise or |
---|
| 305 | ! completely after sunset |
---|
| 306 | PMU0(JCOL) = 0.0_JPRB |
---|
| 307 | CYCLE |
---|
| 308 | ENDIF |
---|
| 309 | |
---|
| 310 | ! Bound the start and end hour angles by sunrise and sunset |
---|
| 311 | ZHOURANGLESTART = MAX(-ZHOURANGLESUNSET, & |
---|
| 312 | & MIN(ZHOURANGLESTART, ZHOURANGLESUNSET)) |
---|
| 313 | ZHOURANGLEEND = MAX(-ZHOURANGLESUNSET, & |
---|
| 314 | & MIN(ZHOURANGLEEND, ZHOURANGLESUNSET)) |
---|
| 315 | ENDIF |
---|
| 316 | |
---|
| 317 | ! Compute average MU0 in the model timestep, although the |
---|
| 318 | ! numerator considers only the time from ZHourAngleStart to |
---|
| 319 | ! ZHourAngleEnd that the sun is above the horizon |
---|
| 320 | PMU0(JCOL) = (ZSINDECSINLAT(JCOL) * (ZHOURANGLEEND-ZHOURANGLESTART) & |
---|
| 321 | & + ZCOSDECCOSLAT(JCOL)*(SIN(ZHOURANGLEEND)-SIN(ZHOURANGLESTART))) & |
---|
| 322 | & / (2.0_JPRB * ZHALFTIMESTEP) |
---|
| 323 | |
---|
| 324 | ! This shouldn't ever result in negative values, but just in |
---|
| 325 | ! case |
---|
| 326 | IF (PMU0(JCOL) < 0.0_JPRB) THEN |
---|
| 327 | PMU0(JCOL) = 0.0_JPRB |
---|
| 328 | ENDIF |
---|
| 329 | |
---|
| 330 | ENDIF |
---|
| 331 | ENDDO |
---|
| 332 | ENDIF |
---|
| 333 | |
---|
| 334 | ENDIF |
---|
| 335 | |
---|
| 336 | |
---|
| 337 | ! ------------------------------------------------------------------ |
---|
| 338 | IF (LHOOK) CALL DR_HOOK('COS_SZA',1,ZHOOK_HANDLE) |
---|
| 339 | END SUBROUTINE COS_SZA |
---|