SURFEX  V8_0
Surfex V8_0 release
 All Classes Files Functions Variables
surf_solar_shadows.F90
Go to the documentation of this file.
1 !SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
2 !SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
3 !SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
4 !SFX_LIC for details. version 1.
5 !#########################################################################
6 SUBROUTINE surf_solar_shadows (PMAP,PXHAT,PYHAT,PCOSZEN,PSINZEN,PAZIMSOL, &
7  pzs,pzs_xy,zxhat_ll,zyhat_ll,iior_ll, &
8  ijor_ll,zzs_ll,zzs_xy_ll,pdirswdt,pdirsrfswd)
9 !#########################################################################
10 !
11 !!**** * SURF_SOLAR_SHADOWS * - computes the modifications to the downwards
12 !! direct solar flux at the surface, due to
13 !! orientation and shape of this surface.
14 !!
15 !! PURPOSE
16 !! -------
17 !!
18 !!
19 !!** METHOD
20 !! ------
21 !!
22 !! EXTERNAL
23 !! --------
24 !!
25 !! IMPLICIT ARGUMENTS
26 !! ------------------
27 !!
28 !!
29 !! REFERENCE
30 !! ---------
31 !!
32 !!
33 !! AUTHOR
34 !! ------
35 !! V. Masson * Meteo-France *
36 !!
37 !! MODIFICATIONS
38 !! -------------
39 !! Original 15/01/02
40 !! V. Masson 01/03/03 add multiple wavelenghts
41 !!
42 !-------------------------------------------------------------------------------
43 !
44 !* 0. DECLARATIONS
45 ! ------------
46 !
47 !
48 !
49 USE modn_io_offline, ONLY: lshadows_other
50 
51 IMPLICIT NONE
52 !
53 !* 0.1 DECLARATIONS OF DUMMY ARGUMENTS :
54 !
55 !
56 REAL, DIMENSION(:,:), INTENT(IN) :: pmap ! map factor
57 REAL, DIMENSION(:), INTENT(IN) :: pxhat ! X coordinate
58 REAL, DIMENSION(:), INTENT(IN) :: pyhat ! Y coordinate
59 REAL, DIMENSION(:,:), INTENT(IN) :: pcoszen ! COS(zenithal solar angle)
60 REAL, DIMENSION(:,:), INTENT(IN) :: psinzen ! SIN(zenithal solar angle)
61 REAL, DIMENSION(:,:), INTENT(IN) :: pazimsol! azimuthal solar angle
62 REAL, DIMENSION(:,:), INTENT(IN) :: pzs ! (resolved) model orography
63 REAL, DIMENSION(:,:), INTENT(IN) :: pzs_xy ! orography at vort. points
64 REAL, DIMENSION(:), INTENT(IN) :: zxhat_ll ! X coordinate (all processors)
65 REAL, DIMENSION(:), INTENT(IN) :: zyhat_ll ! Y coordinate (all processors)
66 INTEGER, INTENT(IN) :: iior_ll ! position of SW corner of current processor domain
67 ! ! in the entire domain (I index along X coordinate)
68 ! ! (both including the 1 point border)
69 INTEGER, INTENT(IN) :: ijor_ll ! position of SW corner of current processor domain
70 ! ! in the entire domain (J index along Y coordinate)
71 REAL, DIMENSION(:,:), INTENT(IN) :: zzs_ll ! orography at center of grid meshes
72 ! ! (all processors)
73 REAL, DIMENSION(:,:), INTENT(IN) :: zzs_xy_ll ! orography at SW corner of grid meshes
74  ! (all processors)
75 
76 !
77 REAL, DIMENSION(:,:,:,:), INTENT(INOUT):: pdirswdt ! SW Flux received by
78  ! each subgrid triangle
79 REAL, DIMENSION(:,:,:), INTENT(OUT) :: pdirsrfswd ! SuRF. DIRect SW Flux
80 !
81 !
82 !* 0.2 DECLARATIONS OF PARAMETERS
83 !
84 REAL, PARAMETER :: xradius = 6371229. ! Radius of the Earth
85 REAL, PARAMETER :: xpi = 4.*atan(1.) ! Pi
86 INTEGER, PARAMETER :: jphext = 1 ! Number of points around the physical area of computation
87 LOGICAL :: lsphere ! tests if computations are done on the Earth sphere
88 ! ! or if the geometry of the domain is supposed flat.
89 !
90 !* 0.3 DECLARATIONS OF LOCAL VARIABLES
91 !
92 INTEGER :: iib, iie, ijb, ije
93 INTEGER :: ji, jj
94 INTEGER :: iresp
95 REAL :: zazim ! azimuthal solar angle
96 REAL :: zcosazim ! azimuthal solar angle cosine
97 REAL :: zsinazim ! azimuthal solar angle sine
98 REAL :: zcoszen ! zenithal solar angle cosine
99 REAL :: zsinzen ! zenithal solar angle sine
100 INTEGER :: jloop
101 INTEGER :: jt
102 !
103 REAL :: zx ! current X, Y and Z coordinates of
104 REAL :: zy ! the sun beam
105 REAL :: zz !
106 REAL :: za ! X, Y coordinates of a point
107 REAL :: zb ! projected according to sun dir.
108 REAL :: zdx ! current distance of the sun beam
109 REAL :: zdy ! to the next X or Y side box limit
110 REAL :: zzi ! initial position of the sun beam
111 REAL :: zzcurv ! offset due to earth curve surface
112 REAL, DIMENSION(3) :: zxt ! X, Y and Z coordinates of corners
113 REAL, DIMENSION(3) :: zyt ! of a triangle
114 REAL, DIMENSION(3) :: zzt !
115 REAL, DIMENSION(3) :: zat ! X, Y coordinates of projected corners
116 REAL, DIMENSION(3) :: zbt ! of a triangle
117 LOGICAL :: gf ! interception flag
118 !
119 !* 0.3 DECLARATIONS OF LOCAL VARIABLES FOR ENTIRE PARALLELIZED DOMAIN
120 !
121 INTEGER :: iiu_ll, iju_ll
122 INTEGER :: iib_ll, iie_ll, ijb_ll, ije_ll
123 INTEGER :: iimax_ll, ijmax_ll
124 INTEGER :: ji_ll, jj_ll ! current grid mesh
125 INTEGER :: ii_ll, ij_ll ! current intercepting grid mesh
126 REAL :: zzs_max_ll ! maximum orography in the domain
127 !
128 
129 !
130 !-------------------------------------------------------------------------------
131 !
132 !* 1. DIMENSION INITIALIZATIONS
133 ! -------------------------
134 !
135 iimax_ll=SIZE(zxhat_ll)-2*jphext
136 ijmax_ll=SIZE(zyhat_ll)-2*jphext
137 iiu_ll = iimax_ll+2*jphext
138 iju_ll = ijmax_ll+2*jphext
139 iib_ll = 1 + jphext
140 iie_ll = iiu_ll - jphext
141 ijb_ll = 1 + jphext
142 ije_ll = iju_ll - jphext
143 !
144 iib = 1 + jphext
145 iie = SIZE(pxhat) - jphext
146 ijb = 1 + jphext
147 ije = SIZE(pyhat) - jphext
148 !
149 !
150 
151 ! PRINT*,SIZE(ZZS_ll,1),SIZE(ZZS_ll,2)
152 
153 zzs_max_ll=maxval(zzs_ll)
154 zzcurv=0.
155 lsphere=any(pmap(iib:iie,ijb:ije)/=pmap(iib,ijb))
156 !
157 pdirsrfswd=0.
158 !
159 !-------------------------------------------------------------------------------
160 !
161 !* 2. LOOP ON THIS PROCESSOR GRID MESH
162 ! --------------------------------
163 !
164 IF (lshadows_other) THEN
165  DO jj=ijb,ije
166  DO ji=iib,iie
167 !
168 !* If zenithal angle greater than Pi/2, sun is down.
169 !
170  IF (pcoszen(ji,jj)<0.) cycle
171 !
172 !
173 !* If zenithal angle is vertical, there is no shadow
174 !
175  IF (pcoszen(ji,jj)==1.) cycle
176 !
177 !-------------------------------------------------------------------------------
178 !
179 !* 3. DEFINITION OF SOLAR ANGLES
180 ! --------------------------
181 !
182 !
183 !* zenithal angle between 0 and Pi/2. Equals zero for sun at zenith.
184 !
185  zcoszen =pcoszen(ji,jj)
186 !
187 !* Azimuthal angle set between 0 and 2*Pi, being equal to zero
188 ! in the East direction, Pi/2 in the North direction, etc...
189 !
190  zsinzen =psinzen(ji,jj)
191 
192 ! PRINT*,ZCOSZEN,ZSINZEN
193 
194 
195  zazim=pazimsol(ji,jj)
196 !
197 !
198 ! ZAZIM = ZAZIM - XPI/2.
199  zazim = xpi/2. - zazim
200 !
201 !* cosine and sine of azimuthal solar angle
202 !
203  zcosazim=cos(zazim)
204  zsinazim=sin(zazim)
205 !
206 !* indices of grid point relative to the entire domain
207 !
208  ji_ll = ji + iior_ll - 1
209  jj_ll = jj + ijor_ll - 1
210 !
211 !-------------------------------------------------------------------------------
212 !
213 !* 4. EXPLORATION OF INTERCEPTING OROGRAPHY FOLLOWING SUN BEAM
214 ! --------------------------------------------------------
215 !
216  DO jt=1,4
217 !
218 !* slope already in its own shadow
219 !
220  IF (all(pdirswdt(ji,jj,jt,:)==0.)) cycle
221 !
222 !* coordinates of the point in the center of the triangle
223 !
224  SELECT CASE (jt)
225  CASE (1)
226  zx=(5.*pxhat(ji)+pxhat(ji+1))/6.
227  zy=0.5*(pyhat(jj)+pyhat(jj+1))
228  zz=(pzs(ji,jj)+pzs_xy(ji,jj)+pzs_xy(ji,jj+1))/3.
229  CASE (2)
230  zx=0.5*(pxhat(ji)+pxhat(ji+1))
231  zy=(5.*pyhat(jj+1)+pyhat(jj))/6.
232  zz=(pzs(ji,jj)+pzs_xy(ji,jj+1)+pzs_xy(ji+1,jj+1))/3.
233  CASE (3)
234  zx=(5.*pxhat(ji+1)+pxhat(ji))/6.
235  zy=0.5*(pyhat(jj)+pyhat(jj+1))
236  zz=(pzs(ji,jj)+pzs_xy(ji+1,jj)+pzs_xy(ji+1,jj+1))/3.
237  CASE (4)
238  zx=0.5*(pxhat(ji)+pxhat(ji+1))
239  zy=(5.*pyhat(jj)+pyhat(jj+1))/6.
240  zz=(pzs(ji,jj)+pzs_xy(ji,jj)+pzs_xy(ji+1,jj))/3.
241  END SELECT
242 !
243 !* projection of this point according to sun direction
244 !
245  CALL proj_solar(zx,zy,zz,za,zb)
246 !
247 !* starts exploration at this point
248 !
249  ii_ll = ji_ll
250  ij_ll = jj_ll
251 !
252  zzi = zz
253 !
254 !
255 !* loop following sun beam until interception
256 !
257  DO jloop=1,2*(iiu_ll+iju_ll)
258 !
259 !* go to the next grid point encountered by the sun beam
260 !
261  IF (zsinazim>0. .AND. zcosazim>0.) THEN
262  zdy=(zyhat_ll(ij_ll+1)-zy)/zsinazim
263  zdx=(zxhat_ll(ii_ll+1)-zx)/zcosazim
264  ELSE IF (zsinazim<0. .AND. zcosazim>0.) THEN
265  zdy=(zyhat_ll(ij_ll) -zy)/zsinazim
266  zdx=(zxhat_ll(ii_ll+1)-zx)/zcosazim
267  ELSE IF (zsinazim>0. .AND. zcosazim<0.) THEN
268  zdy=(zyhat_ll(ij_ll+1)-zy)/zsinazim
269  zdx=(zxhat_ll(ii_ll) -zx)/zcosazim
270  ELSE IF (zsinazim<0. .AND. zcosazim<0.) THEN
271  zdy=(zyhat_ll(ij_ll) -zy)/zsinazim
272  zdx=(zxhat_ll(ii_ll) -zx)/zcosazim
273  ELSE IF (zsinazim==0. .AND. zcosazim<0.) THEN
274  zdx=-(zxhat_ll(ii_ll) -zx)
275  zdy=2.*zdx
276  ELSE IF (zsinazim==0. .AND. zcosazim>0.) THEN
277  zdx=(zxhat_ll(ii_ll+1)-zx)
278  zdy=2.*zdx
279  ELSE IF (zsinazim<0. .AND. zcosazim==0.) THEN
280  zdy=-(zyhat_ll(ij_ll) -zy)
281  zdx=2.*zdy
282  ELSE IF (zsinazim>0. .AND. zcosazim==0.) THEN
283  zdy=(zyhat_ll(ij_ll+1)-zy)
284  zdx=2.*zdy
285  END IF
286 !
287  IF (zdy<zdx) THEN
288  zx = zx + zdy * zcosazim
289  zy = zy + zdy * zsinazim
290  zz = zz + zdy * zcoszen / pmap(ji,jj) / zsinzen
291  IF (zsinazim>0.) THEN
292  ij_ll = ij_ll + 1
293  ELSE
294  ij_ll = ij_ll - 1
295  END IF
296  ELSE
297  zx = zx + zdx * zcosazim
298  zy = zy + zdx * zsinazim
299  zz = zz + zdx * zcoszen / pmap(ji,jj) / zsinzen
300  IF (zcosazim>0.) THEN
301  ii_ll = ii_ll + 1
302  ELSE
303  ii_ll = ii_ll - 1
304  END IF
305  END IF
306 !
307 !* effect of curvature of the earth surface
308 !
309  IF (lsphere) zzcurv = ((zz-zzi)*zsinzen/zcoszen)**2 &
310  / (2.*xradius)
311 !
312 !* sun beam goes to high
313 !
314  IF (zz > zzs_max_ll-zzcurv) EXIT
315 !
316 !* sun beam goes outside of the domain
317 !
318  IF ( ii_ll<iib_ll .OR. ii_ll>iie_ll &
319  .OR. ij_ll<ijb_ll .OR. ij_ll>ije_ll ) EXIT
320 !
321 !
322 !* treatment where the sun beam is currently located
323 !
324  IF (.NOT. ( zzs_ll(ii_ll ,ij_ll ) -zzcurv < zz .AND. &
325  zzs_xy_ll(ii_ll ,ij_ll ) -zzcurv < zz .AND. &
326  zzs_xy_ll(ii_ll+1,ij_ll ) -zzcurv < zz .AND. &
327  zzs_xy_ll(ii_ll ,ij_ll+1) -zzcurv < zz .AND. &
328  zzs_xy_ll(ii_ll+1,ij_ll+1) -zzcurv < zz ) ) THEN
329 !
330 !* exploration of the 4 triangles in the grid mesh encoutered by the sun beam
331 !
332 !* left triangle
333 !
334  zxt(1) =0.5*(zxhat_ll(ii_ll)+zxhat_ll(ii_ll+1))
335  zyt(1) =0.5*(zyhat_ll(ij_ll)+zyhat_ll(ij_ll+1))
336  zzt(1) =zzs_ll(ii_ll ,ij_ll ) - zzcurv
337  zxt(2) =zxhat_ll(ii_ll)
338  zyt(2) =zyhat_ll(ij_ll)
339  zzt(2) =zzs_xy_ll(ii_ll ,ij_ll ) - zzcurv
340  zxt(3) =zxhat_ll(ii_ll)
341  zyt(3) =zyhat_ll(ij_ll+1)
342  zzt(3) =zzs_xy_ll(ii_ll ,ij_ll+1) - zzcurv
343  CALL proj_solar(zxt(1),zyt(1),zzt(1),zat(1),zbt(1))
344  CALL proj_solar(zxt(2),zyt(2),zzt(2),zat(2),zbt(2))
345  CALL proj_solar(zxt(3),zyt(3),zzt(3),zat(3),zbt(3))
346  CALL solar_interc(zat(:),zbt(:),za,zb,gf)
347  IF (gf) THEN
348  pdirswdt(ji,jj,jt,:)=0.
349  EXIT
350  END IF
351 !
352 !* top triangle
353 !
354  zat(2) =zat(3)
355  zbt(2) =zbt(3)
356  zxt(3) =zxhat_ll(ii_ll+1)
357  zyt(3) =zyhat_ll(ij_ll+1)
358  zzt(3) =zzs_xy_ll(ii_ll+1,ij_ll+1) - zzcurv
359  CALL proj_solar(zxt(3),zyt(3),zzt(3),zat(3),zbt(3))
360  CALL solar_interc(zat(:),zbt(:),za,zb,gf)
361  IF (gf) THEN
362  pdirswdt(ji,jj,jt,:)=0.
363  EXIT
364  END IF
365 !
366 !* right triangle
367 !
368  zat(2) =zat(3)
369  zbt(2) =zbt(3)
370  zxt(3) =zxhat_ll(ii_ll+1)
371  zyt(3) =zyhat_ll(ij_ll)
372  zzt(3) =zzs_xy_ll(ii_ll+1,ij_ll ) - zzcurv
373  CALL proj_solar(zxt(3),zyt(3),zzt(3),zat(3),zbt(3))
374  CALL solar_interc(zat(:),zbt(:),za,zb,gf)
375  IF (gf) THEN
376  pdirswdt(ji,jj,jt,:)=0.
377  EXIT
378  END IF
379 !
380 !* bottom triangle
381 !
382  zat(2) =zat(3)
383  zbt(2) =zbt(3)
384  zxt(3) =zxhat_ll(ii_ll)
385  zyt(3) =zyhat_ll(ij_ll)
386  zzt(3) =zzs_xy_ll(ii_ll ,ij_ll ) - zzcurv
387  CALL proj_solar(zxt(3),zyt(3),zzt(3),zat(3),zbt(3))
388  CALL solar_interc(zat(:),zbt(:),za,zb,gf)
389  IF (gf) THEN
390  pdirswdt(ji,jj,jt,:)=0.
391  EXIT
392  END IF
393  END IF
394  END DO
395 !
396  END DO
397  END DO
398  END DO
399 !
400 
401 END IF
402 
403 !-------------------------------------------------------------------------------
404 !
405 !* 6. FLUX FOR THE GRID MESH RECEIVED OVER AN HORIZONTAL SURFACE
406 ! ----------------------------------------------------------
407 !
408 !* there is no problem of summation, because only fluxes projected
409 ! over horizontal surfaces are considered here.
410 !
411 DO jj=ijb,ije
412  DO ji=iib,iie
413  pdirsrfswd(ji,jj,:) = ( pdirswdt(ji,jj,1,:) &
414  + pdirswdt(ji,jj,2,:) &
415  + pdirswdt(ji,jj,3,:) &
416  + pdirswdt(ji,jj,4,:))&
417  * 0.25
418  END DO
419 END DO
420 !
421 !-------------------------------------------------------------------------------
422 !-------------------------------------------------------------------------------
423 !
424  CONTAINS
425 !
426 SUBROUTINE proj_solar(PX,PY,PZ,PA,PB)
427 REAL, INTENT(IN) :: px ! X coordinate (conformal plane)
428 REAL, INTENT(IN) :: py ! Y coordinate (conformal plane)
429 REAL, INTENT(IN) :: pz ! Z coordinate (conformal plane)
430 REAL, INTENT(OUT) :: pa ! X coordinate (plane perp. sun)
431 REAL, INTENT(OUT) :: pb ! Y coordinate (plane perp. sun)
432 
433 pa = zcoszen*zcosazim*px + zcoszen*zsinazim*py - zsinzen*pz
434 pb = - zsinazim*px + zcosazim*py
435 !
436 END SUBROUTINE proj_solar
437 !
438 SUBROUTINE solar_interc(PA,PB,PAD,PBD,OF)
439 REAL, DIMENSION(3), INTENT(IN) :: pa ! X coordinate
440 REAL, DIMENSION(3), INTENT(IN) :: pb ! Y coordinate
441 REAL, INTENT(IN) :: pad ! Z coordinate
442 REAL, INTENT(IN) :: pbd ! X coordinate
443 LOGICAL, INTENT(OUT) :: of ! true if interception occurs
444 !
445 REAL :: zc1, zc2, zc3 ! coefficients of the c1 A + c2 B + c3 = 0 line
446 ! ! defined by two corners of the triangle.
447 !
448 of=.false.
449 !
450 !* first side
451 !
452 zc1=(pb(1)-pb(2))
453 zc2=-(pa(1)-pa(2))
454 zc3=-pb(1)*pa(2)+pb(2)*pa(1)
455 !
456 IF ( (zc1*pad+zc2*pbd+zc3)*(zc1*pa(3)+zc2*pb(3)+zc3) <0.) RETURN
457 !
458 !* second side
459 !
460 zc1=(pb(3)-pb(2))
461 zc2=-(pa(3)-pa(2))
462 zc3=-pb(3)*pa(2)+pb(2)*pa(3)
463 !
464 IF ( (zc1*pad+zc2*pbd+zc3)*(zc1*pa(1)+zc2*pb(1)+zc3) <0.) RETURN
465 !
466 !* third side
467 !
468 zc1=(pb(3)-pb(1))
469 zc2=-(pa(3)-pa(1))
470 zc3=-pb(3)*pa(1)+pb(1)*pa(3)
471 !
472 IF ( (zc1*pad+zc2*pbd+zc3)*(zc1*pa(2)+zc2*pb(2)+zc3) <0.) RETURN
473 !
474 of=.true.
475 !
476 END SUBROUTINE solar_interc
477 !
478 END SUBROUTINE surf_solar_shadows
modn_io_offline
Definition: modn_io_offline.F90:6
proj_solar
subroutine proj_solar(PX, PY, PZ, PA, PB)
Definition: surf_solar_shadows.F90:426
surf_solar_shadows
subroutine surf_solar_shadows(PMAP, PXHAT, PYHAT, PCOSZEN, PSINZEN, PAZIMSOL, PZS, PZS_XY, ZXHAT_ll, ZYHAT_ll, IIOR_ll, IJOR_ll, ZZS_ll, ZZS_XY_ll, PDIRSWDT, PDIRSRFSWD)
Definition: surf_solar_shadows.F90:6
solar_interc
subroutine solar_interc(PA, PB, PAD, PBD, OF)
Definition: surf_solar_shadows.F90:438
Generated on Fri May 13 2016 09:31:00 for SURFEX by   doxygen 1.8.5