phys/module_ra_eclipse.F

   1 module module_ra_eclipse
   2
   3 contains
   4
   5 !--------------------------------------------
   6 ! Solar eclipse prediction moduel
   7 !
   8 ! Alex Montornes and Bernat Codina. University of Barcelona. 2015
   9 !
  10 ! Based on Montornes A., Codina B., Zack J., Sola, Y.: Implementation of the
  11 ! Bessel's method for solar eclipses prediction in the WRF-ARW model. 2015
  12 !
  13 ! This key point of this subroutine is the evaluation of the degree of obscuration,
  14 ! i.e. the part of the solar disk that is hidden by the Moon. This variable is
  15 ! introduced in the SW schemes and used as a correction of the incoming radiation.
  16 ! In other words, if obscur is the part of the solar disk hidden by the Moon, then
  17 ! (1-obscur) is the part of the solar disk emitting radiation.
  18 !
  19 ! Note that if obscur = 0., then we recover the standard code
  20 !
  21 !----------------------------
  22 ! History
  23 !
  24 !       amontornes      2015/09         First implementation
  25 !
  26 !----------------------------
  27 ! MAIN ROUTINE
  28 !
  29 subroutine solar_eclipse(ims,ime,jms,jme,its,ite,jts,jte,    &
  30                           julian,gmt,year,xtime,radt,        &
  31                           degrad,xlon,xlat,Obscur,mask,      &
  32                           elat_track,elon_track,sw_eclipse   )
  33 !-------------------------------------------------------------
  34   implicit none
  35 !
  36 !--------------
  37 !----------
  38 ! VARIABLES
  39 !----------
  40 !--------------
  41 !
  42 !----------
  43 ! Input/output vars
  44 !----------
  45 ! Grid indeces
  46   integer,                             intent(in)    :: ims,ime,jms,jme, &
  47                                                         its,ite,jts,jte
  48 ! Temporal variables
  49   real,                                intent(in)    :: gmt,degrad,julian,xtime,radt
  50   integer,                             intent(in)    :: year
  51 ! Geographical information
  52   real,    dimension(ims:ime,jms:jme), intent(in)    :: xlat,xlon
  53   real,    dimension(ims:ime,jms:jme), intent(inout) :: Obscur
  54 ! Eclipse variables
  55   integer, dimension(ims:ime,jms:jme), intent(inout) :: mask
  56   real,                                intent(inout) :: elat_track,elon_track
  57 ! Namelist option 1 - enabled / 0 - disabled
  58   integer,                             intent(in)    :: sw_eclipse
  59
  60 !----------
  61 ! Local vars
  62 !----------
  63 !
  64 ! Constants
  65   real,    parameter       :: pi      = 3.1415926535897932384626433
  66   real,    parameter       :: epsil   = 0.0818192 ! Earth excentricity
  67 !
  68 ! Besseliam elements from file eclipse_besselian_elements.dat
  69 !       X, Y         :: coordinates of the eclipse's axis in the fundamental plane
  70 !       d, mu        :: declination and azimuth (deg)
  71 !       l1, l2       :: radii of the penumbra and umbra regions in the fundamental
  72 !                       plane
  73 !       tanf1, tanf2 :: angle of the shadow's conus
  74 !       t0, min, max :: time when the Besselian elements are valid
  75 !       Dt           :: time correction
  76   integer, parameter       :: NBessel = 4 ! Dimension
  77   real, dimension(NBessel) :: X, Y, d, l1, l2, mu
  78   real                     :: tanf1, tanf2, t0, tmin, tmax, Dt
  79 !
  80 ! d and mu in radians
  81   real                     :: rmu, rd
  82 !
  83 ! Besselian elements at the current time t1
  84   real                 :: cx, cy, cd, cl1, cl2, cmu
  85 !
  86 ! Current time t1, and hourly angle H
  87   real                 :: t1, &
  88                            H
  89 !
  90 ! Latitude and longitude in radians
  91   real                 :: rlat, rlon
  92 !
  93 !
  94 ! Geocentric coordinates
  95 !----------
  96   real                 :: cos_lat1, sin_lat1 ! cos/sin of the corrected latitude
  97   real                 :: xi, eta1, zeta1      ! Geocentric coordinates at the observation point
  98 !
  99 ! Eclipse variables
 100   real                 :: DELTA, DELTA2      ! Distance of the observation point
 101                                              ! to the eclipse axis
 102   real                 :: LL1, LL2           ! Penumbra and Umbra radius
 103                                              ! from the point observer's plane
 104   logical              :: eclipse_exist
 105 !
 106 ! Loop indeces
 107   integer              :: i, j
 108 !
 109 ! Internal variables
 110   real                 :: A, B
 111   real                 :: da,eot,xt24
 112   integer              :: julday
 113 !
 114 ! Eclipse path
 115   real :: E, E_1, tan_rd1, rd1, rho1, cy1, tan_gamma, cgamma, sin_beta, &
 116           tan_C, C, cm, tan_lat1, tan_lat, sin_H,                       &
 117           cos_H, tan_H, beta, lon, tan_rd2, rd2, rho2
 118 !
 119 ! Missing value
 120   real :: MISSING
 121
 122
 123 !--------------
 124 !----------
 125 ! CODE CUSTOMIZATION
 126 !----------
 127 !--------------
 128 !
 129   MISSING = 1E30
 130
 131 !--------------
 132 !----------
 133 ! MAIN
 134 !----------
 135 !--------------
 136 !
 137 !----------
 138 ! Check if the eclipse is enabled by the namelist.input
 139 ! If eclipse is not enabled -> set eclipse variables to zero and nothing to do
 140   if(sw_eclipse .eq. 0) then
 141      Obscur(:,:) = 0.
 142      mask(:,:)   = 0
 143      elat_track  = MISSING
 144      elon_track  = MISSING
 145      return
 146   endif
 147 !
 148 !----------
 149 ! Set current time
 150   xt24   = mod(xtime+radt,1440.)
 151   t1     = gmt+xt24/60
 152   julday = int(julian)+1
 153   if(t1 .ge. 24.) then
 154         t1     = t1 - 24.
 155   endif
 156 !
 157 !----------
 158 ! Check if the eclipse exists for the current simulation date
 159   eclipse_exist = .false.
 160 !
 161 ! and load the besselian elements
 162 ! Note: eclipse_besselian_elements.dat should be in the running directory!
 163   call load_besselian_elements(t1, julday, year, X, Y, d, l1, l2, mu, &
 164                                tanf1, tanf2, t0, tmin, tmax, NBessel, &
 165                                eclipse_exist,Dt)
 166 !
 167 !----------
 168 ! If the eclipse does not exist, then we set all the output vars to zero and finish
 169   if(.not. eclipse_exist) then
 170      Obscur(:,:) = 0.
 171      mask(:,:)   = 0
 172      elat_track  = MISSING
 173      elon_track  = MISSING
 174      return
 175   endif
 176 !
 177 !----------
 178 ! Compute the besselian elements at t1
 179   call compute_besselian_t(t1, NBessel, X, Y, d, l1, l2, mu, t0, Dt, &
 180                                cx, cy, cd, cl1, cl2, cmu)
 181 !
 182 !----------
 183 ! Convert the beselian elements from deg to radians
 184   rmu = cmu * degrad
 185   rd  = cd  * degrad
 186
 187 !---------------------------
 188 ! STEP 1: Compute the eclipse track. This is just a diagnosis for the user
 189 !
 190 !-- Initialization
 191   elat_track=MISSING
 192   elon_track=MISSING
 193 !
 194 !-- Compute the excentricity correction
 195   E         = sqrt(1 - epsil**2)
 196   E_1       = 1/E
 197 !-- Get d1
 198   tan_rd1   = E_1 * tan(rd)
 199   rd1       = atan(tan_rd1)
 200 !-- Get rho1
 201   rho1      = sin(rd)/sin(rd1)
 202 !-- Correct cy1
 203   cy1       = cy / rho1
 204 !-- Get d2
 205   tan_rd2   = E * tan(rd)
 206   rd2       = atan(tan_rd2)
 207 !-- Get rho2
 208   rho2      = cos(rd)/cos(rd2)
 209 !-- Compute gamma
 210   tan_gamma = cx / cy1
 211   cgamma    = atan(tan_gamma)
 212 !-- Compute beta
 213   sin_beta  = cx/sin(cgamma)
 214 !
 215 !-- If |sin(beta)|>1 then the axis of the eclipse is tangent to the Earth's
 216 !   surface or directly it does not crosses the Earth.
 217   if(abs(sin_beta) .lt. 1) then
 218         beta      = asin(sin_beta)
 219         !-- Compute C
 220         tan_C     = cy1 / cos(beta)
 221         C         = atan(tan_C)
 222         !-- Compute c
 223         cm        = cy1 / sin(C)
 224         !-- Compute lat1
 225         sin_lat1  = cm * sin(C + rd1)
 226         cos_lat1  = sqrt(1 - sin_lat1**2)
 227         tan_lat1  = sin_lat1/cos_lat1
 228         !-- Convert to geographic latitude
 229         tan_lat   = E_1 * tan_lat1
 230         elat_track= atan(tan_lat)/degrad
 231         !-- Compute the hourly angle
 232         sin_H     = cx / cos_lat1
 233         cos_H     = cm * cos(C+rd1)/cos_lat1
 234         tan_H     = sin_H/cos_H
 235         H         = atan(tan_H)
 236
 237         if(cos_H .lt. 0) then
 238                 if(sin_H .ge. 0) then
 239                         H = pi - abs(H)
 240                 else
 241                         H = pi + abs(H)
 242                 end if
 243         end if
 244         !-- Compute the longitude
 245         elon_track= (rmu - H)/degrad
 246         !-- Convert to geopgraphic longitude
 247         if(elon_track .gt. 360) then
 248                 elon_track = elon_track -360
 249         endif
 250         if(elon_track .gt. 180) then
 251                 elon_track = elon_track -360
 252         endif
 253         if(elon_track .le. -180) then
 254                 elon_track = elon_track +360
 255         endif
 256         !-- Correct the Earth's movement
 257         elon_track = -elon_track + 1.002738 * 15*Dt/3600
 258   endif
 259
 260 !---------------------------
 261 ! STEP 2: Compute the obscuricity for each grid-point. This step it is important
 262 !         in order to reduce the incoming radiation
 263 ! Loop over latitudes
 264   do j=jts,jte
 265      ! Loop over longitudes
 266      do i=its,ite
 267         !--Correct the Earth's movement
 268         lon = -xlon(i,j) - 1.002738 * 15*Dt/3600
 269         if(lon .lt. 0) then
 270             lon = lon + 360
 271         endif
 272         if(lon .gt. 360) then
 273             lon = lon -360
 274         endif
 275
 276         !-- Convert degrees to radians
 277         rlat = xlat(i,j)*degrad
 278         rlon = lon*degrad
 279
 280         !-- Hourly angle
 281         H    = rmu - rlon
 282
 283         !-- Adjust the Earth's sphericity
 284         A    = sqrt( 1 - ( epsil * sin(rlat) )**2 )
 285
 286         !-- Compute lat1
 287         cos_lat1 = cos(rlat) / A
 288         sin_lat1 = sin(rlat) * E / A
 289
 290         !-- Geocentric coordinates at the observation point
 291         xi     = cos_lat1 * sin(H)
 292         eta1   = sin_lat1 * cos(rd1) * E - cos_lat1 * sin(rd1) * cos(H)
 293         zeta1  = sin_lat1 * sin(rd2) * E + cos_lat1 * cos(rd2) * cos(H)
 294
 295         !-- Correct cy (sembla que no es necessari)
 296         cy1  = cy/rho1
 297
 298         !-- Distance to the shadow axis
 299         DELTA2 = (xi - cx)**2 + (eta1  - cy1)**2;
 300         DELTA  = sqrt(DELTA2)
 301
 302         !-- Penumbra and Umbra cone's circumference over the Earth's surface
 303         LL1    = (cl1 - zeta1 * tanf1) ! Penumbra
 304         LL2    = (cl2 - zeta1 * tanf2) ! Umbra
 305
 306         ! Check if we are inside the eclipse
 307         if(DELTA .gt. LL1) then ! Without obscuration
 308               Obscur(i,j) = 0.
 309               mask(i,j)   = 0
 310         elseif(DELTA .le. LL1 .and. DELTA .gt. abs(LL2)) then   ! Penumbra
 311               Obscur(i,j) = (LL1 - DELTA) / (LL1 + LL2)         ! Partial eclipse
 312               mask(i,j)   = 1
 313         elseif(DELTA .le. abs(LL2)) then                        ! Umbra
 314                 if(LL2 .lt. 0) then
 315                       Obscur(i,j) = 1.                          ! Total eclipse
 316                       mask(i,j)   = 2
 317                 else
 318                       Obscur(i,j) = (LL1 - DELTA) / (LL1 + LL2) ! Annular eclipse
 319                       mask(i,j)   = 3
 320                 endif
 321         endif
 322
 323         ! Check that the obscuration is not greater than 1
 324         if(Obscur(i,j) .gt. 1) then
 325               Obscur(i,j) = 1.
 326         endif
 327
 328      enddo
 329      ! End loop over longitudes
 330   enddo
 331 ! End loop over longitudes
 332 end subroutine
 333 !
 334 !----------------------------
 335 ! INTERNAL SUBROUTINES
 336 !       load_besselian_elements --> load the besselian elements from eclipse_besselian_elements.dat
 337 !       compute_besselian_t     --> compute the besselian elements for the current time
 338 !
 339 !----------------------------
 340 subroutine load_besselian_elements(t, julday, year, X, Y, d, l1, l2, mu,  &
 341                                    tanf1, tanf2, t0, tmin, tmax, NBessel, &
 342                                    eclipse_exist,Dt)
 343   IMPLICIT NONE
 344 !
 345 !--------------
 346 !----------
 347 ! VARIABLES
 348 !----------
 349 !--------------
 350 !
 351 !----------
 352 ! Input/output vars
 353 !----------
 354   integer,                             intent(in)    :: julday, year, NBessel
 355   real,                                intent(in)    :: t
 356   logical,                             intent(inout) :: eclipse_exist
 357   real, dimension(NBessel),            intent(inout) :: X, Y, d, l1, l2, mu
 358   real,                                intent(inout) :: tanf1, tanf2, t0, tmin, tmax, Dt
 359 !----------
 360 ! Local vars
 361 !----------
 362   integer :: ryear, rjulday
 363   logical :: file_exist
 364
 365 !----------
 366 ! Initialization
 367 !----------
 368   eclipse_exist = .false.
 369
 370 !----------
 371 ! Load the information from eclipse_besselian_elements.dat
 372 !----------
 373   inquire( file="eclipse_besselian_elements.dat", exist=file_exist)
 374   if (file_exist) then
 375          open (20,FILE='eclipse_besselian_elements.dat',status="old",action="read")
 376   else
 377         call wrf_error_fatal('load_besselian_elements: eclipse_besselian_elements.dat not found')
 378         return
 379   end if
 380
 381   do
 382      ! Load the list of eclipse features
 383      read(20,*,end=1)ryear,rjulday,t0,tmin,tmax,  &
 384                          X(1), X(2), X(3), X(4),  &
 385                          Y(1), Y(2), Y(3), Y(4),  &
 386                          d(1), d(2), d(3), d(4),  &
 387                         l1(1),l1(2),l1(3),l1(4),  &
 388                         l2(1),l2(2),l2(3),l2(4),  &
 389                         mu(1),mu(2),mu(3),mu(4),  &
 390                         tanf1,tanf2,Dt
 391      ! Check if the current date matches with the loaded episode
 392      if(ryear .eq. year .and. rjulday .eq. julday .and. tmin .le. t .and. t .le. tmax) then
 393              eclipse_exist = .true.
 394              close(20)
 395              return
 396      else
 397      endif
 398   end do
 399   1 close(20)
 400   return
 401 end subroutine
 402 !----------------------------
 403 !----------------------------
 404 !----------------------------
 405 subroutine compute_besselian_t(t1, NBessel, X, Y, d, l1, l2, mu, t0, Dt, &
 406                                cx, cy, cd, cl1, cl2, cmu)
 407   IMPLICIT NONE
 408   integer,                             intent(in)    :: NBessel
 409   real,                                intent(in)    :: t1
 410   real, dimension(NBessel),            intent(in)    :: X, Y, d, l1, l2, mu
 411   real,                                intent(in)    :: t0, Dt
 412   real,                                intent(inout) :: cx, cy, cd, cl1, cl2, cmu
 413
 414   real                                               :: t
 415   integer                                            :: i
 416
 417   cx  = 0.
 418   cy  = 0.
 419   cd  = 0.
 420   cl1 = 0.
 421   cl2 = 0.
 422   cmu = 0.
 423
 424   t = (t1 + DT/3600) - t0
 425
 426   do i=1,NBessel
 427       cx  = cx  + X(i)  * t**(i-1)
 428       cy  = cy  + Y(i)  * t**(i-1)
 429       cd  = cd  + d(i)  * t**(i-1)
 430       cl1 = cl1 + l1(i) * t**(i-1)
 431       cl2 = cl2 + l2(i) * t**(i-1)
 432       cmu = cmu + mu(i) * t**(i-1)
 433   enddo
 434
 435 end subroutine
 436 !----------------------------
 437 !----------------------------
 438
 439 end module module_ra_eclipse