source: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/ifs/cos_sza.F90 @ 5218

Last change on this file since 5218 was 5159, checked in by abarral, 5 months ago

Put dimensions.h and paramet.h into modules
File size: 13.1 KB
Line 
1SUBROUTINE 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
79USE PARKIND1 , ONLY : JPIM, JPRB
80USE YOMHOOK  , ONLY : LHOOK, DR_HOOK
81USE YOMCST   , ONLY : RPI, RDAY
82USE YOMRIP   , ONLY : YRRIP
83USE YOERIP   , ONLY : YRERIP
84USE YOERAD   , ONLY : YRERAD
85USE YOMLUN   , ONLY : NULOUT
86
87!     ------------------------------------------------------------------
88
89IMPLICIT NONE
90
91INTEGER(KIND=JPIM),INTENT(IN) :: KSTART      ! Start column to process
92INTEGER(KIND=JPIM),INTENT(IN) :: KEND        ! Last column to process
93INTEGER(KIND=JPIM),INTENT(IN) :: KCOL        ! Number of columns in arrays
94REAL(KIND=JPRB),   INTENT(IN) :: PGEMU(KCOL) ! Sine of latitude
95REAL(KIND=JPRB),   INTENT(IN) :: PGELAM(KCOL)! Longitude in radians
96LOGICAL, INTENT(IN) :: LDRADIATIONTIMESTEP   ! Is this for a radiation timestep?
97REAL(KIND=JPRB),  INTENT(OUT) :: PMU0(KCOL)  ! Cosine of solar zenith angle
98
99! Solar time at the start and end of the time interval
100REAL(KIND=JPRB) :: ZSOLARTIMESTART, ZSOLARTIMEEND
101
102! The time of half a model/radiation timestep, in radians
103REAL(KIND=JPRB) :: ZHALFTIMESTEP
104
105! For efficiency we precompute sin(solar declination)*sin(latitude)
106REAL(KIND=JPRB) :: ZSINDECSINLAT(KSTART:KEND)
107!...and cos(solar declination)*cos(latitude)
108REAL(KIND=JPRB) :: ZCOSDECCOSLAT(KSTART:KEND)
109! ...and cosine of latitude
110REAL(KIND=JPRB) :: ZCOSLAT(KSTART:KEND)
111
112! Tangent of solar declination
113REAL(KIND=JPRB) :: ZTANDEC
114
115! Hour angles (=local solar time in radians plus pi)
116REAL(KIND=JPRB) :: ZHOURANGLESTART, ZHOURANGLEEND
117REAL(KIND=JPRB) :: ZHOURANGLESUNSET, ZCOSHOURANGLESUNSET
118
119INTEGER(KIND=JPIM) :: JCOL        ! Column index
120
121REAL(KIND=JPRB) :: ZHOOK_HANDLE
122
123!     ------------------------------------------------------------------
124IF (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
128IF (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')
131ENDIF
132
133DO JCOL = KSTART,KEND
134  ZCOSLAT(JCOL) = SQRT(1.0_JPRB - PGEMU(JCOL)**2)
135ENDDO
136
137IF (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
239ELSE
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
334ENDIF
335
336
337!     ------------------------------------------------------------------
338IF (LHOOK) CALL DR_HOOK('COS_SZA',1,ZHOOK_HANDLE)
339END SUBROUTINE COS_SZA
Note: See TracBrowser for help on using the repository browser.