Gust calculation
342 views
Skip to first unread message
Will H
Feb 24, 2024, 12:17:15 PM2/24/24
to wrfpython-talk
Has anyone ever figured out how to calculate wind gusts from wrf output?
Leo
Feb 24, 2024, 2:26:57 PM2/24/24
to wrfpython-talk, hathew...@gmail.com
I would also be interested in knowing how to calculate wind gusts using the output of the WRF model.
Will H
Feb 24, 2024, 9:34:07 PM2/24/24
to Leo, wrfpython-talk
I have a function but I can't get it to work properly.Â
It's the NAM/GFS Fortran code
Leonardo Jiménez
Feb 25, 2024, 1:15:00 AM2/25/24
to Will H, wrfpython-talk
I found some equations in the article ( https://repositorio.aemet.es/bitstream/20.500.11765/5238/1/P11-trabajo%20Vindel%20et%20al.pdf ).Â
If you can share the equation that you found that uses NAM/GFS Fortran code.Â
Marco Miani
Feb 25, 2024, 3:58:37 AM2/25/24
to Leonardo Jiménez, Will H, wrfpython-talk
Hi all.Â
I start admitting my ignorance with this, but still I'd like to ask my question:Â
Does it really make any sense to calculate wind gusts, using HOURLY wind data?Â
Wind gusts, typically, are calculated relying on records sampled at minutes intervals. Instead, in our cases (WRF and paper posted here), deal with hourly wind data (model outputs).Â
Just a thought....
--
You received this message because you are subscribed to the Google Groups "wrfpython-talk" group.
To unsubscribe from this group and stop receiving emails from it, send an email to wrfpython-tal...@ucar.edu.
To view this discussion on the web visit https://groups.google.com/a/ucar.edu/d/msgid/wrfpython-talk/CADkHQmbPsRAzVj6HVW4adLXfHma-AMrAkWTEaArKgou8%2B06rfg%40mail.gmail.com.
Will H
Feb 25, 2024, 5:43:50 AM2/25/24
to wrfpython-talk, Marco Miani, Will H, wrfpython-talk, leonardoj...@gmail.com
Well for my case I sample the WRF data every 15 mins. But the thought I have is that by calculating gusts you could see the trend for energy production (wind farms) or SFC fronts. Â
I'm also making a python meteogram using wrf data and the gusts would be part of that as well for a specific location using lat long. Â
As for the codes I can share them later. The are written in Fortran though.Â
Leonardo Jiménez
Feb 25, 2024, 7:06:14 AM2/25/24
to wrfpython-talk, Marco Miani
I have some knowledge of the Fortran language and I am interested in the topic.
Reading the code is fine.
Will H
Feb 25, 2024, 7:55:04 AM2/25/24
to wrfpython-talk, leonardoj...@gmail.com, wrfpython-talk, Marco Miani
FROTRAN CODE 1
   !
   !  Within PBL depth, calculate excess wind speed over surface speed
   !  at each level. Multiply this excess by a coefficient (f(z)) that
   !  decreases with height from 1.0 to 0.5 at 1 km height, and is 0.5          <-  code does not follow this statement
   !  for any height > 1 km. Add the maximum weighted wind excess back  Â
   !  to the surface wind. [gust = vsfc + max (f(z)*(v(k)-vsfc) ]
   !
   !  One primary difference between the calculation below and the others is
   !  that the individual components of the wind at each level are used rather
   !  then a simple calculation of the speed.
   !
   !$OMP PARALLEL DO  &
   !$OMP PRIVATE (i,j)
   DO j = jts,min(jte,jde-1)
   DO i = its,min(ite,ide-1)
     !  The following code is adapted from the RUC and appears to be similar
     !  to the simpler NAM/GFS calculation
     !
     gust10m(i,j) = 0.
     speed = sqrt(u10(i,j)**2 + v10(i,j)**2)
     IF (kpbl(i,j) .gt. 1) THEN
      DO k=kpbl(i,j),1,-1
        zagl = 0.50*(z2(i,k+1,j)+z2(i,k,j)) - ht(i,j)  ! Z at mid point
        wght = max (0.5,(1-(zagl/2000.)))
        u0   = wght * (u2(i,k,j)-u10(i,j))
        v0   = wght * (v2(i,k,j)-v10(i,j))
        uu   = u0+u10(i,j)
        vv   = v0+v10(i,j)
        speed = max (speed,sqrt(uu*uu+vv*vv))
      END DO
      k   = kpbl(i,j)
      wght = 1.0 - min(0.5,pblh(i,j)/2000.)
      u0  = u10(i,j) + wght*(u2(i,k,j)-u10(i,j))
      v0  = v10(i,j) + wght*(v2(i,k,j)-v10(i,j))
      speed = sqrt(u10(i,j)**2 + v10(i,j)**2)
      speed = max (speed,sqrt(u0*u0+v0*v0))
     END IF
     gust10m(i,j) = speed
     g10_max(i,j) = max (g10_max(i,j),speed)
   END DO
   END DO
   !$OMP END PARALLEL DO
FORTRAN CODE 2:
SUBROUTINE CALGUST(LPBL,ZPBL,GUST)
C$$$ SUBPROGRAM DOCUMENTATION BLOCK
C . . .
C SUBPROGRAM: CALGUST COMPUTE MAX WIND LEVEL
C PRGRMMR: MANIKIN ORG: W/NP2 DATE: 97-03-04
C
C ABSTRACT:
C THIS ROUTINE COMPUTES SURFACE WIND GUST BY MIXING
C DOWN MOMENTUM FROM THE LEVEL AT THE HEIGHT OF THE PBL
C
C
C PROGRAM HISTORY LOG:
C 03-10-15 GEOFF MANIKIN
C 05-03-09 H CHUANG - WRF VERSION
C 05-06-30 R ROZUMALSKI - DYNAMIC MEMORY ALLOCATION AND SMP
C THREAD-SAFE VERSION
C
C USAGE: CALL CALGUST(GUST)
C INPUT ARGUMENT LIST:
C NONE
C
C OUTPUT ARGUMENT LIST:
C GUST - SPEED OF THE MAXIMUM SFC WIND GUST
C
C OUTPUT FILES:
C NONE
C
C SUBPROGRAMS CALLED:
C UTILITIES:
C H2V
C
C LIBRARY:
C COMMON -
C LOOPS
C OPTIONS
C MASKS
C INDX
C
C ATTRIBUTES:
C LANGUAGE: FORTRAN 90
C MACHINE : CRAY C-90
C$$$
C
C
use vrbls3d
use vrbls2d
C
C INCLUDE ETA GRID DIMENSIONS. SET/DERIVE PARAMETERS.
C
!
INCLUDE "params"
C
INCLUDE "CTLBLK.comm"
C
C DECLARE VARIABLES.
C
INTEGER :: LPBL(IM,JM)
REAL :: GUST(IM,JM)
REAL ZPBL(IM,jsta_2l:jend_2u)
C
C
C*****************************************************************************
C START CALMXW HERE.
C
C LOOP OVER THE GRID.
C
DO J=JSTA,JEND
DO I=1,IM
! GUST(I,J) = SPVAL
GUST(I,J) = 0.
ENDDO
ENDDO
C
C ASSUME THAT U AND V HAVE UPDATED HALOS
C
!$omp parallel do
!$omp& private(ie,iw,mxww,u0,v0,wind)
DO 20 J=JSTA_M,JEND_M
DO 20 I=2,IM-1
L=LPBL(I,J)
IF(MODELNAME .EQ. 'NMM')THEN
IE=I+MOD(J+1,2)
IW=I+MOD(J+1,2)-1
USFC=D25*(U10(I,J-1)+U10(IW,J)+
X U10(IE,J)+U10(I,J+1))
VSFC=D25*(V10(I,J-1)+V10(IW,J)+
X V10(IE,J)+V10(I,J+1))
SFCWIND=SQRT(USFC**2 + VSFC**2)
U0 = D25*(U(I,J-1,L)+U(IW,J,L)+
X U(IE,J,L)+U(I,J+1,L))
V0 = D25*(V(I,J-1,L)+V(IW,J,L)+
X V(IE,J,L)+V(I,J+1,L))
WIND=SQRT(U0**2 + V0**2)
ELSE IF(MODELNAME .EQ. 'NCAR')THEN
USFC=U10(I,J)
VSFC=V10(I,J)
SFCWIND=SQRT(USFC**2 + VSFC**2)
U0=U(I,J,L)
V0=V(I,J,L)
WIND=SQRT(U0**2 + V0**2)
END IF
DELWIND=WIND - SFCWIND
ZSFC=FIS(I,J)*GI
DELWIND=DELWIND*(1.0-AMIN1(0.5,ZPBL(I,J)/2000.))
GUST(I,J)=SFCWIND+DELWIND
10 CONTINUE
20 CONTINUE
C END OF ROUTINE.
C
RETURN
END
   !  Within PBL depth, calculate excess wind speed over surface speed
   !  at each level. Multiply this excess by a coefficient (f(z)) that
   !  decreases with height from 1.0 to 0.5 at 1 km height, and is 0.5          <-  code does not follow this statement
   !  for any height > 1 km. Add the maximum weighted wind excess back  Â
   !  to the surface wind. [gust = vsfc + max (f(z)*(v(k)-vsfc) ]
   !
   !  One primary difference between the calculation below and the others is
   !  that the individual components of the wind at each level are used rather
   !  then a simple calculation of the speed.
   !
   !$OMP PARALLEL DO  &
   !$OMP PRIVATE (i,j)
   DO j = jts,min(jte,jde-1)
   DO i = its,min(ite,ide-1)
     !  The following code is adapted from the RUC and appears to be similar
     !  to the simpler NAM/GFS calculation
     !
     gust10m(i,j) = 0.
     speed = sqrt(u10(i,j)**2 + v10(i,j)**2)
     IF (kpbl(i,j) .gt. 1) THEN
      DO k=kpbl(i,j),1,-1
        zagl = 0.50*(z2(i,k+1,j)+z2(i,k,j)) - ht(i,j)  ! Z at mid point
        wght = max (0.5,(1-(zagl/2000.)))
        u0   = wght * (u2(i,k,j)-u10(i,j))
        v0   = wght * (v2(i,k,j)-v10(i,j))
        uu   = u0+u10(i,j)
        vv   = v0+v10(i,j)
        speed = max (speed,sqrt(uu*uu+vv*vv))
      END DO
      k   = kpbl(i,j)
      wght = 1.0 - min(0.5,pblh(i,j)/2000.)
      u0  = u10(i,j) + wght*(u2(i,k,j)-u10(i,j))
      v0  = v10(i,j) + wght*(v2(i,k,j)-v10(i,j))
      speed = sqrt(u10(i,j)**2 + v10(i,j)**2)
      speed = max (speed,sqrt(u0*u0+v0*v0))
     END IF
     gust10m(i,j) = speed
     g10_max(i,j) = max (g10_max(i,j),speed)
   END DO
   END DO
   !$OMP END PARALLEL DO
FORTRAN CODE 2:
SUBROUTINE CALGUST(LPBL,ZPBL,GUST)
C$$$ SUBPROGRAM DOCUMENTATION BLOCK
C . . .
C SUBPROGRAM: CALGUST COMPUTE MAX WIND LEVEL
C PRGRMMR: MANIKIN ORG: W/NP2 DATE: 97-03-04
C
C ABSTRACT:
C THIS ROUTINE COMPUTES SURFACE WIND GUST BY MIXING
C DOWN MOMENTUM FROM THE LEVEL AT THE HEIGHT OF THE PBL
C
C
C PROGRAM HISTORY LOG:
C 03-10-15 GEOFF MANIKIN
C 05-03-09 H CHUANG - WRF VERSION
C 05-06-30 R ROZUMALSKI - DYNAMIC MEMORY ALLOCATION AND SMP
C THREAD-SAFE VERSION
C
C USAGE: CALL CALGUST(GUST)
C INPUT ARGUMENT LIST:
C NONE
C
C OUTPUT ARGUMENT LIST:
C GUST - SPEED OF THE MAXIMUM SFC WIND GUST
C
C OUTPUT FILES:
C NONE
C
C SUBPROGRAMS CALLED:
C UTILITIES:
C H2V
C
C LIBRARY:
C COMMON -
C LOOPS
C OPTIONS
C MASKS
C INDX
C
C ATTRIBUTES:
C LANGUAGE: FORTRAN 90
C MACHINE : CRAY C-90
C$$$
C
C
use vrbls3d
use vrbls2d
C
C INCLUDE ETA GRID DIMENSIONS. SET/DERIVE PARAMETERS.
C
!
INCLUDE "params"
C
INCLUDE "CTLBLK.comm"
C
C DECLARE VARIABLES.
C
INTEGER :: LPBL(IM,JM)
REAL :: GUST(IM,JM)
REAL ZPBL(IM,jsta_2l:jend_2u)
C
C
C*****************************************************************************
C START CALMXW HERE.
C
C LOOP OVER THE GRID.
C
DO J=JSTA,JEND
DO I=1,IM
! GUST(I,J) = SPVAL
GUST(I,J) = 0.
ENDDO
ENDDO
C
C ASSUME THAT U AND V HAVE UPDATED HALOS
C
!$omp parallel do
!$omp& private(ie,iw,mxww,u0,v0,wind)
DO 20 J=JSTA_M,JEND_M
DO 20 I=2,IM-1
L=LPBL(I,J)
IF(MODELNAME .EQ. 'NMM')THEN
IE=I+MOD(J+1,2)
IW=I+MOD(J+1,2)-1
USFC=D25*(U10(I,J-1)+U10(IW,J)+
X U10(IE,J)+U10(I,J+1))
VSFC=D25*(V10(I,J-1)+V10(IW,J)+
X V10(IE,J)+V10(I,J+1))
SFCWIND=SQRT(USFC**2 + VSFC**2)
U0 = D25*(U(I,J-1,L)+U(IW,J,L)+
X U(IE,J,L)+U(I,J+1,L))
V0 = D25*(V(I,J-1,L)+V(IW,J,L)+
X V(IE,J,L)+V(I,J+1,L))
WIND=SQRT(U0**2 + V0**2)
ELSE IF(MODELNAME .EQ. 'NCAR')THEN
USFC=U10(I,J)
VSFC=V10(I,J)
SFCWIND=SQRT(USFC**2 + VSFC**2)
U0=U(I,J,L)
V0=V(I,J,L)
WIND=SQRT(U0**2 + V0**2)
END IF
DELWIND=WIND - SFCWIND
ZSFC=FIS(I,J)*GI
DELWIND=DELWIND*(1.0-AMIN1(0.5,ZPBL(I,J)/2000.))
GUST(I,J)=SFCWIND+DELWIND
10 CONTINUE
20 CONTINUE
C END OF ROUTINE.
C
RETURN
END
Leonardo Jiménez
Feb 26, 2024, 10:54:09 AM2/26/24
to Matias Ezequiel Suarez, Will H, wrfpython-talk, Marco Miani
It's always interesting to read your Python code.
El lunes, 26 de febrero de 2024, Matias Ezequiel Suarez <suarezm...@gmail.com> escribió:
El lunes, 26 de febrero de 2024, Matias Ezequiel Suarez <suarezm...@gmail.com> escribió:
Hello Will, here in Córdoba, Argentina, we are using the following parameterization to calculate wind gusts from WRF output:<img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAA60AAANzCAIAAADN82nXAAAAAXNSR0IArs4c6QAAAANzQklUCAgI2+FP4AAAIABJREFUeJzsnXmcFcW5939PVZ9lzuzDwDDDvim7KCouiIoLIMaYG3dNXKLXuJt4XRITs1yNMVf0xkTjlhu97lEiaoxBg6igrLIzgDIgMDAw+z5n6a7n/aO6+/Q5M2wm8d73Ut9PPmQ8p091bV391K+eeoqYGQaDwWAwGAwGwyGG+J/OgMFgMBgMBoPB8D+AsYMNBoPBYDAYDIcixg42GAwGg8FgMByKGDvYYDAYDAaDwXAoYuxgg8FgMBgMBsOhiLGDDQaDwWAwGAyHIsYONhgMBoPBYDAcihg72GAwGAwGg8FwKGLsYIPBYDAYDAbDoYixgw0Gg8FgMBgMhyLGDjYYDAaDwWAwHIoYO9hgMBgMBoPBcChi7GCDwWAwGAwGw6GIsYMNBoPBYDAYDIcixg42GAwGg8FgMByKGDvYYDAYDAaDwXAoYuxgg8FgMBgMBsOhiLGDDQaDwWAwGAyHIsYONhgMBoPBYDAcihg72GAwGAwGg8FwKGLsYIPBYDAYDAbDoYixgw0Gg8FgMBgMhyLGDjYYDAaDwWAwHIoYO9hgMBgMBoPBcChi7GCDwWAwGAwGw6GIsYMNBoPBYDAYDIcixg42GAwGg8FgMByKGDvYYDAYDAaDwXAoYuxgg8FgMBgMBsOhiLGDDQaDwWAwGAyHIsYONhgMBoPBYDAcilhfzW0YAED7vYyZaL9X/cNggMEUyNc/+t4KQE+TjR4+D5b9S9XD3u4FHHD9/51wt9xw4G8Z+PfvugszEXVvOw0F62HfxT7wb/9B1XfwzbrX/sPMRPuvS++O6XR6Ksq+es7fxVfT7VwOthQH3E8OnK+0vPunx+wwI7MPBvrGP3H4/TKtw0wAiPggf5tNtyLv/UoAAP0va0eDwfBP5X+XHvxVGsHuHf9pox0H0P+Z9W/wSmSWfd/1wMxKqcB/pj/MSnnfiRzglQf1QwYzQykoVoqVw+ww26xsVolkqvsvu6e2n/S9VzVD1ySYwf7/e9V+0KX6Sjio7s3dQLdelPX533lHw/84f3/X3W8K++gRwXn4P+nuBwUR/UM6sHkIDAbDPviK9ODgQOSZffoVLogoY/QkVqyYiEAKxGCCACAABlsgZi0TQHxJEzatTFBg4PasKwBwAAYIIG+icOB3YoAZSVIKShEIpNNhnTy5fzMUsSKiCIRvjDMUAIJQSjlCAZCQhLSex+AESIEhSNuc/qtCEZgQYiUhBZhASXaISP8tWOir0mXf+8shIAt592X2VRm3otK2u2JmkGRwEmQrahdgoEWJxkZsb8DWrbu37qhbs2ZNpHHpSy/+ureTtKRFbr0qIrd20z3AS5NYgEgR+/lIgRUpxRKAre9Nut5ItxQBkhwB4WrPBPIUaKVU9muVXM0p+KG2Kt3e59VV9+bXVaTgACB2+zAzC9Fd48/o5/t9rzOUYsVkOVBJQJIEmECsVy8IBNgEBktvMkCAbl2LlSBBYAIRbADMMvu+7ABgtx0D7atFM7c+3dx3y63NYLB1YJ1nH4+NAgDW+nS3u+h5nf6cHWYmYaGn6WKmyij8a7y7E4O7jxLpqRRzsK2zekIgfTAzkwAgFPxa6oFMndVvd1/dzEo/Y8mhh9oOpqBzwwRyvN/6fZuClUPE6U8UEen27d45u5OZH32f7J9kPCyArmH3cVDKu0WG+svMjjv6KADECPRS+NdkJi4AOKwEiWAeDnrEP2Blt7sSbERhg+EQ4auwgzNfPIA73ukPCdljq3uNglJMgfGRfYtZv4j+Iat4fsaCQqPDYGYBEoIOfE0tUF6AwMxKv/q6SyTaxmIIEswsqGcHCf99G3QAYLDjOBCkU/arlEkxs2ISUjAze9XFpCuwh/T3XSeBi3vKW6YLh9JaLLB8+fJf/Pa/urq6dndIpVQn9QVgI8rM508ejWw7Zt95IDC77yRmAEqXUSmlVPWetkQi0dmVYGYrHMrJySkuDOfn5xMpJmZFQgiJdN3v1wiAp6zv8wLemzmbVWmBDxkH2Vddrw9mPRuQJBSr4AVKKACOcryL3RsKZghI0rPHnnOVPe3sqYB7+crtRfsuy5d4JLPmohl/B6Yu+8558Ifph6Knp3dvOdy3V5L7AP4d400PPaRbeffdrwAwWC8FMYGIyLUUkXWN/tLP7n7Xl5A5G9G/7WnikTGYM4GQrm39iGWNGMEfetO2HuYn3cYcAiBI+D/p/q33nz2Mz+kPv8TwbTAYDjG+Gj3YNR0ybYj08B2UE+IQikQb4AD1CdTXIyxRUoLSCAkgn0BACA6BqNsL4MAQSkFkugX6llw7YIOaCAqUaEI4jCF5ad3lQFAAC16zXb3yyh+bwhWxWEzJCAApZW5ubmtrKxEJ5di2HWnfUVpaeuu3pliEaEBBiRMYcGABiAIEWCAAKcBx6PWFmz799NM2WZGbm5siGY1GbVvZth3mrq6urlJZc9P1l4aIFKNJQSkUC4QEirySMjNxD2/07HdeUJPzrdJ0nQV+qBNjBUCQzM2JDRoytLKyson7kaBmWQbA4gSAXmVlTkY9ZZChDJEAYBMTURxgoJGovRPvfFK/evXqtes27N69W6ikn9uEiABIyahlWQP6VvTv339MuVVWVnbltGERy7IApVRPZrDWxRXA3t86A3trW/2t41sJnq3Julp8bTvo0cgAExSlpyPeCkM3j0mvOyqlmoR8/Kk/bu/qbVmWjJUwsxWJWpbV1dEJQCIFQLINoDVuA8jPz+/Tp88xI0uHDhG9AQnEXDVc55AC7S4CRo/wL/DMPApq2Ah0DNcyAagnzS/AfjxBmVmn53UoV7NM/5ZE9pUeROQADJbc3bZTfp79r+QBDBFpu40BZiV6MIiZGcQEctxacrXMgAK6r/S9ulXMrEgCkG6yTqDsSNdAWuF25+daee0CKaBJwFEghZBEIYQEwlAgSHfdJqhXB2qSHWRqq15/E/5DnVWh/lOf+UAobwWCWNkkhG4v5c/V/bGCASBBYKJ2wAG1MwDkESygMKBnZ0mvvnRAUN6iBAd6CDNLeE8rUXpNIzB78oroZp+AHr1+M9YGs77t4Qk1GAz/F/mn28HehLxnZSLLPtYraKlU6q33lr/11ltrtzYBEKyEEJNG97344otPOXKElK5q+w+c5nsvG7ZBL78855X3PqmtrRXxxM033zxo5nFCHNytmLm+vn7evHm11BuATSH9uf/WFMohooJU7bhx4/CtKcHUGZy0qbp6Vzgnr6ur67CBvYUQARUWVVVV7777bjPKANjCYmYhLKVUmLsA9M9tuu66Sz98f9mcOXNWfbErlUqN6d/rZz/7WWG/wozV3P3KJMxBKc6nu4mQtUg8duzofx83OuXg579fNWfOHNdKcAAgGo0CEEL0uL+tRxQrm8WuXbW/efGdjz76qNEpBQCyCgoKJh4xpk+fPpZltbW1bfhi5+bNm5MO27a9ffv2bdu2rXN2x2Kxq8/6lRZWhRDdLe+9cbCKZg/r7D1elnZvoH33Xp3CJ598snZ3mIiSlENEjrZQGUQUtINTXu8CkMd1paWl37v0jBkzpmrzK7ueA823t/ISEWfKz3tTBA+cffxkH0llLeUHu9w/UOPzNftgboLjEnyF3lue8LXMA1lk0ANLZvLo/p97W2pw80DEzIro/fcXPvPm+5s3b3ZSiUGDBj384+8NHFDuP1Dd60VvJKCenuXuOQn25L1dHhwBhBCO47jDI6G7dqtxlPO3D1bMmTNn/Re7ksnk4f1Lf/rTn46tKEqPq93UcT/zvj7iNUSPF+6vXPu8i8FgOMQ5oKXGfxw971vXGk8HyGF88Dl++cuHK2uTABRFASgWACx0Abj+a+Ouu3bGAIABi/W68MFGIAjutVdKKe1N2g5SwJurcNddP2gMVwDok9h+yy23XDejQpKIHHDq2j/YBlIp5+Oq+meeeebt9TaADlEKuF7Rxanqs88++65vjuzbtzjP9ekEgDiggDdWJO68886OUC8Arz9y/eHDe/VxU2YACVA8kVxWVfvhhx/+7q0tgKuofePo3CuuuGLIkF47drScc+0sAG2ylIj6JD+/8sor77jwSEGIBPwdu6vyGaXw4gz4Hoqu7zILMDMFNDyvIfXbWrFKkXAYf9uCG264vS46kplDTjuAey4bd8lFp5WrhBBCQFtv3feqKwApCAAtQFcSD722+fnnn09SsVKqCNv69u37k389c9Kko3KQFpdTkB0dnX+ct/7555/fEu9HRDmpxlGjRv3pwUtDUsYUa1fR9CvcK6T+/72YCNnemT21tWIwIHxPR9XNxHcABZWEUIADEJCrNf69eiLaDLYRcpTakxQrVqz43q8/bGtra7HKAECliGhYeNudd945ONr2xRdfvPL+2k2bNu22hgJgCjNzkV03ffr0WTedYAkUuE6/Ac0sU/TusXS+bygR6b+1rzZ7RiIRATaA4EQ6reEB+9OD0/cN+lhnabq+aQjPYYmyZLtgHWbM7IIaYbaR1/2+wa98n1e40zckvH8dRhiQBN33Mlf2RUYizEwqK30/swEXrHRZuvtke+OyAmBDMPBJNa666vra0GEAwtwB4O5LjvzWZWf1YlvnJ1AWVz1lZld1pfTdg/Wf2fcy9gP4ng8UXLsLzOj0rCAJMJAEFCAVLOHWjy57A7BlW+OM7z4shOigQgClyW1XXXXVHReNFRARUPq+nO2Lr58vAQlAKZuIuo/2e/fiVWDO8oPPvmQvPz7AnmwwGP5v8FU85Ps1tRnMzAr485/fu+22H9bU1GT9UGsPQojZs2e/9fZaZz8WzH7IWmnVL1cFNDa33n//L4lI34t9DiZxd7sRIRyW40aV3XffnYMHDw7eVGszkycf37dvsQj4wAFgQClesmSJEEJKSUTz5s1jzwIG3F130Uh43Oj+11576YQJE3SCxxxzzE9/etvwob0kIR6P6+rSKQBIpVLdayAz4kR2EbOUqgModcY1RIjF3M97SBw9JJ7VlIq5sYVvu+3u559/3tPCafTo0Y8/ftcJxx9lCUghA/9DQX7sX8495vHHf11WVqbTrKioCMp15G1l29sdeyzUvsseaDjPCTjzpcpgBaWUSjn45JOVV1xx/VNPv+gote8K1YlIIXKjOP74o0477TQ/t7pE5eXlxx3ba/z4weecc8qv//OmSy+9NFhMAH/9619feul1N3Sd5z4cNNTIo+cMdP+8myfrl2a/2iS6PaHs2lzkRwWBu9af/vsAn9J0VSDrFgGBPPCVYsWArdTCT9ZdeeXNTz/9rOqWVPC36d61j/S9rhjI/L76pM4DM3d1MQDLsvQABU/uzVyR6IGsVA+q/nv4NjANcRQ7bv2sueKKG5966g/BpQRdsHg8rodTKaVlWUSUSqX2llW/NoL/AhBiXzPSvfHVqjwGg+H/S76iya43uomsbd3aszIO6oR45u3P7vnNHMduDYdSJw3uPGOkGprfUOJUCYIgdFFBJ/KbrYEP/fffmhS1A8wEiKA9dzD5cZgdhgNSSSAJqgW+98sXv+jMa5DlCQ4nOMxCspDiIL0iiBTDsQAL6AWUC/X9y84qTu0BFBE7bDGFGVj+6SrAIYCgiJXOTxzoJHpv2cY2WdQuSttF6dsfrmsDHJDtVheHwRFwLpAHVZRLIXTmqV3fv3pGCdlFULnAYRV5OTJuqbYoOiPcEeH2M046yvKEO1fkgfCWGpXOmCcBQn8C6CyxnhUAIAjSEpjQqozXc0h5oRUAdoiVBCwgL4IQkkKxBVKkHDiC9W+CTooiS+CxIWyIZqBe0aU/eulvVdFmWdYsy8JOc3E09dgPzh8eQRG4EBxmOwJH/y8PnAeuAMYVOz+4Zmaxs0Nwom/vQtd1UfiOi771p7TY0+3NqtIhAnTuyBNs3XgNe7ka2teTKaAuJYAEaJsSry2pO+uOty66971Ve3K6iiYoCB0WJO0QCaUDI8CrnxBUiJ08oACqf+/8MHcpkPKcZUtyEQFyVTIfdjnUbZccfc4Jo4rsuhRFUhRpsfo3y36/fuXjze1wvJAaXtMzw2E4zE6gofdmRerphySSetVb15troEAoZJmPIIDhtXN2DWXjCpPsCpluNbtzUsVQILd/KnIUOWAtSiuwA8UEgmCQ24GFCNyLBVjoz3swnli7pwSLrDzpFLqMiqEY9cBuJZ5dan/t9pcu//c/fVLb28rv35U5P2fKNjO9+umWvreXlP3L2IH33AX7AEOx7qLsELNgSFYWeEQ/ypedYac9yp05dmPMaTp98pFRQJJFLIgFQXrldWuDiIiU39ZexJFA/fuCdLduQEQCkiCC7ej2ATCBmoA6Rc8tUV+77eVv3/vGwtoyWTywPXOaEQGGDSjKlYkId+QgEVFdUW6dNuXIECwJynzivI103rxXr+EoZQMqc9wIZNJ7lgP173iO5tqV3f82UA/ej3tUkr2Pe76jwWD4P8ZXsU8uSw6B52/nw8C6dZWPPPJYbm7u5VdfdeaZIwtCANDGmDt3871PvBfUD9ra2jZu/GzcmMO6J3uAZGk/DNiO/eJrn65evZrlaByAvrJflLtvTECI444rLi4uru5009RRfhcuXKiuGCOFzHqn7tq1u6amRgjXvm9oaFixsm3QkfnkZpUJpNytRaisrAQwfPjwwQMgvNX5WCx2ww03PP744+22XVBQcPW3rh462OrRlA8u4wZLDe+lmKUK97i+7H3i7dEh12T3PPncnwghtH8wUbZ3bEZtk2sH/fa3T1VV7Q7ql1dccUVpabH06yHDUxnQzgZSnnhcaX5+fmNLvKCg4EBa6supm26rZf40a9cUAzU19Rfe/EBra2uDNVxKySmOaZ18f4mTdo6UEp5fdbDz6Lu7vuYQELjyyhkfffQRvIYTQiSTyY8+2jxh+rBAE8hgf/ar/aBqQLc1Mr0UMi7wYwKwFjH3Gq5r3/f102dvn5n/q6ynknsKa4DMXyFgXQUN1L10Zl3VqG1OXnvttZs7+wEQIsrMubm5e8twsDKDD3V2+l4wByLy91y632ReExSJ9apRbgw33HDDPU/OT6VShYWF559//pAh/Shw9yz8p9h76A6uof1c+WlzwAOewXvqkt/97nc/7xxAREQhZo7FYt7PWAfDARCLxa677rrHHnusM5UqLCy8+rKrhwyS1MPtMnqj3+568ar7MoLX7un/3Hd/9vuJP2IYDAYDvrr4wb6uBi/CgPcVM4eJKvr0PmNs0Q9+8IM+JVEhWDAA2ERDzhzK9rRHH320lkYys4MIgN21HaPHQClbStnzjP5AMgOkEALQBKzd0PL4c/PCOcOQ1LaFQyBFSgfqOqitFQShlWpmFpQSzPkUnnHSqM1/3cGKO6iYiOIo3LYnvrkmNaBfKEJenE5CO/CXD9d2yWIASRJCiE5Z+vrfPp125CkSyHdNQEXEHcDu+lRLQ60kOvfUCflAWO/gYUSIbjhr8Hen3ZdwnHA4HGIlAJGxk4YBViAdi5QBSgcwdu0YRVKXhJDtkecbIq65AHcLEQgOBHl+o5FcKOEoQUqpMAiKw+EQAHYgLe3wp4Kr8wzFxDakAhZtxnPvfNYRGU5EwkkJEqVWzcUzeucoRwoJV4PTUeEEEQjM4BAkgHzgxPED3vtwde9caengppn51zvnAXfXv/5Eur4xIiDjCQASaXvX09l03ABiwIZkPwwqSRAsHa0ZCkA7RGPK2tZeBFHUFi4goqhj2SSVAAOWthX8PhNsAgJDsLcOQbDBKV1dkgUAKbRkpZg5THAcZ0TfyNjBsY+qkwAciihQXOSu3LQ9MWOYBHJIuuUisiH1LQjp9mUdkDqNIL1X340OQA48zdbtP3A8L17yBhFfh4T2+Cct9SIUeIC6Wz9MwpdItV+4N2nTTqhgSkuRpO1+cseQFISC0na58LyGoJ9fQd5pZOn5iX7SFOl6EORWPCwIJrY8s0uBUwSW2B0PV7cXNof7AwhzFzMnhEgCSv+WAMACKaXI2w7qeXDsPX1PEdZl96skCQHy4wFLHSgBrjapr6MS4LKZE8+bfqRSKhQKSaIQB8L3AkRwIPw+qcjdK0m+lR3oycKNNKzgCcUciN/guLVKBEj3edGtIhhwQMzYlYhs7yhsjvRjZl0/Nkkn/VuAKAYoxvUzBl175i/ith2JRMKsAkHZ/TlShs+07j8K0p1UEQQYYH8sYmZJQinlubu7C10Est00OaPggEPs9hMSftzl/dJ9JmMwGP6P8T+/7qO14bKy3rNm3du7NFcIIUCSSBIJQAoxc+bQnJwcIpLSHQbTkQcOXq/NWv5TzM0tuPfeeyORyI033qg/1L5ozNzZ2fnlSpR2JyAi4Mwzzwzajjr9RYsWcTfJecGCBf5lOqsLFixob2/3L2KwUoqBlStX6jo55ZRT/JVN8na1WFJGw2HyQkdlueS6siLYYVbMSkEvBLvps3L0/5SjvPXc4MsgqM2QtzuHmRXYZsUMpdi2e3adzDKjs2DAcdRTTz0f/AkzH3XUUZFQeqk0KAp6QRjSt9MO2aFQKK0taS9ZsGKv4MrRf7h2vwIz9LcK0MfgqUAGMvuMcpgdpRzlKKUc5djKcWvMcbIOt8uqc3g+nToDGWXfy5pJXl5edzEMgTe0dt0ZPnx41gXNzc0Z9Q9vJd496s9r356eoazysttP0gvQDjsOO7or+uGNtS+y8lRrfY3jKKXSq949oruf0nnzTkVM1zMr9/NAUwJw4ChW7l2Uo49p4IA069ePW//ujdhRipkd5fj10L3+0z7fmRKv49WPr82n9xKAdUXtN/10xwAzw3HYL7v3ITvKcZSje5SuWwCCRNiyIuGwa8UGvMbd+gcrVkqxYvdMR2Y4SvuBub3a/Tzoae1lRrFyFDuOX6tKZXosEIhAzPCjT2Yl4nh/67IodkekkJQ5kYjeEZGeEGQ6Riul3F7qtbsK9lX2fVdIR00OvgJch2nWrcy6P+gxzSu7Yys7OKYZDAYDvjI92CVTCdYwc5QIgMOKiAS0KQJmDgnJzDESo4b02vVZXJDIU/WFhYUnjh+cB3hROA8yC0TMbIOYuYOQYrpv1rO7GuM/+MnFeXlJHek2yREAEpLtLzFcpuNRECQz5xGOGFo0tE+0pqYmTkXM3GUVMPO8xWtmnndaCSzSrpZA1R6s39qqIv0AgAUzOmQJFN74cOM5M4/OYaVtI5tVAli4vNJi54gxYwcWUQjQAU4TAAPaeO8CFGNjZUN1dfVlZx4hCCEQgC4STIgDANoInXEsXrgpmUyed9Y4Atph7a7DvE93t7W1DSnLmTSp/6AoJJCnPWZdAVXvAAOBkgCD4kSKURtHZWXtph0dbW1tNfG8ZjHEfUu5vhIAoCQzlOU6MqqANgYALUBtc3LB+p2wSlMsiShXtQE4ZtzEKCDdSL3pVnf37LsWPBGhGDh/6lFjyvMnHtXHAgNIgBjUqWuGiPWrWginGckkD+xDDCQEoDU5rWcTpVJoaHCGlchoGHlEShcXaAM5RE2MrV+kNlXVtrW1OSxisVhhSX55edFh+SgpplxIBtqAdlHQFYoxs5bAOkRJp+xVTUIAIi4dBxU5CFuIej1TlwIgHR+NGQQhWctibv0B2gcXgkGAIqV1spKCXMF15J7Ax5IRFq63aAIAqB1QwLY2rFq1e2tNezQaPXJ4r9GjY30IAtDuGrpmkwATJQEC2gFF1NCFefPW5ajWc885QQC2wtvLOrdt2zasb/Skk0ZUkHBX9wmtEA6wqxWrVtVs3tMejUYnjCgePSqvN0gCuQGHB6356b7aSVBAa4pWrKit2bLysounh0FMVJ3A4sW7tu9py8/PP+3ofmW9UeL/FmiA/KwqsXz9nng8Prwi99hj+/WPQAI5rpbvdw04RF4UWzQ6tHFjc+WWho6OjsLCwiFDBh83EhZQ5MYmFwCnAMXUCSStAhuCwQIhImoPlzcCNol4HDkJ5OWhRIJc8Rs6Sm6jg32kr1VhhwR7kSjaCUpScwJ/+1tlDC3nzDwBIIfor0vjW7duHVERnjJlRCn0ugAUI0FQQBdjw4bWxu3rZk4/oZiVINFJcBT+unTH4YcPzCsCAWGWzHAIJEh3Xz1fySVBQBGzIAgIJk6BdI91iHZ1YOVKt4dM0D0EJIBYYP0kQeQQOhgJma/j+dgUUUp1hMsbGUmS8Thy4sjPR64E4I02AAOb1zVVV1dfNG2cJFiBdUKbwELoNmoh1OyhVZW7mpubUzbl5eWVlcZGjeo1IA8WkA8mkF57AUkFThAxEAcxURMQT9EnC6qSyeSFM8YB6CDaU4f5y+ra2toGlcUmTeo/OAd+bxTepomsyBgGg+EQ4au1g/cOMwt9bpbnE0ZEjrfa5TiOUlKvb952222xWM7fs1LFzFrBsoE5c95etGjR1KlTTzyxYtmyenjuaF9ObEbmOpqnXkAIcfrpp+vQB0IIxQ4RVVZWdsYBLySbYrVgQWVWCtpqf+edd86eebTvlkdEKQfLli0Dck477TQtB2oNaXdd2+7duzfvaqyurl7z+RefffZZV0d8woQJl5w+njyXPAfo6Ois3FJbWVm5aF3VqlWrnASmTZtmzxhXU1P78BNzFi9e3G5VKKWiqjkcDt999Yyvf30aMkVlAkGrRyQSicQ77y+fO3fu0k27HMeJi1IAzdRbezl3qw0KppP+yrts0aJFwbLrfwcNGsTp+2Yckpy1cElARUVZeUUZBezltesqX5/7UUdHR22HnUgkGlvbmpububX9vPPOu+W731y3bsObcz90HKc9hdbW1ub2rvr6+pbmVmZ+7alZgwf1dQ/nAxzFKdCrr855Zs7HjY2NNvIBKEgALBxmHhRufPPN5wCk7NSf/rzyww8/1H6Tbm4dPPbYYy8/uiGZTMYcBeC9Pz5RXJSHgHUYlLfJDZXAPXbFLIVYRwVhd1meAPTt25e8Qz5sRy1fveW1116bv7omlUqlZBGAXNUwePDgX9x60ciRQ5jS/rhMZCtOvL1SAAAgAElEQVSntYM3b968umrb0qVLl1fuSaVS35w+2VYnNDWkfvzjHy/blgSQoxqPOeaYJ39+vRBEDMdRy1ZXzZ49+6MVu5i5nfIA5KJ+yJAh991y4ciRw4JeqgpQSjW0JLZv3768cuu6deuWrNmWSqUmHN73wgunJ+P83HPP/eGtFalUypb5AJ4XdRdeeOGNl5wmhLCZPv10zS+ffnPbtm1dKAWQw80FBQUP3nnJ0Ucf4al96frR99pdn3zxxRff/HBFZ2dnHIV+BQ4t6br55punHXuYFK4DScrGnDlvP7PIFdSFEMTEzL/5zW9++1sKx3cIIXJTXbNnzy4ucH0LDjB9r4mhlGrtUFVVVSs3Vy9duvTTyp3JZPL8mZNtPqGhIXXPPfd8usUGkKNqjz322Ed/fkN9fXP17tbq6uotu+s2bty49vPqzs7OE8YPnjH9BACOchTJ+vrme+65h4gciufm5oacBAAbkogcz3slHo+HU22DBw9+7cmH4I0bishRatmqqtmzZ3+4cpdSKkEFzJzLjUOGDLnv1vNHjhzmPqqB+nnWqx8iAkNK+cgjj/zmNxRO7QKQl0q8+uqrbamWPXv2bNnVsmPHjlWbt2zatMluT02YMOHCM8emvXSZ4Y32cUfMmfPmS+8s3rVrVxJFABghABKdQohTJ/S7/PLLJ44olSLDI85hdHR0bthat379+gWVVatWrbITmD59eurMcbt31z385OuLFy/uoDIiCjl6TJt27rkzvN5otGGD4ZDmK7aDtS4oAv8GdweDiGwoELS/YieBCdsbsPKzplyoWCz2uzsumXRsjmc6ZkiDBwhDMMEmwcD6HZj19LuD++bdecsFUaDAaRFQBFIEBjNJ5+AP6yBidl0rXS/GKOCAZ0we/fp/NzXzELDe5kYd6P3BsrbBk/MsUBJIsnh//ocMMbRvZOfOnSnZF0BC5ABYVbVjUw3KyrWCajkClRvR2pWTF+o446T+USgCbJDNuOiW/2hoaKgPDYV7ytq4XqEdLVZZSmuLIIfVt+54as2aNfXhIQCSKkQ0IRJq3xMe8dvXdzzzzDPtqRDJw+MUg0SzVQ7gh4+/35Q/8TtTS7WC4rdiC8hm+v3cmmeeeaaruTEcDn/nm2cOGzZs0/a6t99+u7WTAYeFZAAq3c3I26TCrBucFBhCepGDseGzmgTF0kawShBRWWm+nwQzEUnvzDMnsGzNRFoa5/TdwCnQp6vXvPBeJYBGayAQs6kIGNBHbrMLD7OB9Zs+e23ucmZuiAwFkBS9ifrkWLUARE5+Am481BZYnUn6zs/nrF69mlSksHjkrFtmjBgxYsuO2mefffbjTTYIBRUFnYAFrFz/xUNPzCYix+rv27WOiNk2N/NAJrbE7oKCgvyCWMBhRf+f4/tNEqCPxQ5sNRO6VORdBwjFyibUNrXqs8oEJwmUw23HjR+qY7uuqsaDDz61YtNOADaVUoi6ZIFSqtmqqN7Jl/3w97/85b3TDicBROCA0Ax5/iU3rm+tANBpFTDnh5CnpGorGrNsN26/44nGxlin1R9AiWMt3VDbKgjAxmo8+OBTyz/bCcAWZUQUpxiAZirfWY1v//CZ++//9zMPI0mIQjG4EfKqa36wpCYPQKcoZA4RjYCF4UUj5q7Hgw/Mampq6sBAWO5T0Mx9fvnKukTe4OnTD3vk4b8sWLCg3Sqh8OFxFSYioBxxXP3TF3/72yNOHUzC04N17IBdwAcLtt/z6OttbW3FUBddcsWkI4e9+uqrf/m0jYhWtIQv//mbs245e/q0MX1AAD7Z0PTAk7N3RkaxKCBFAGwhmLnDLgYQJkGg3sX10UKkAALqgPkfbfvRo290dHQUKufi/aXfDFx02U1rm8sAdIpi5twQDeMQt+WPXbYT/3bXk42Nsc5wLwDFNhZtrK1nXHrzjz9v7c3MXbKEKJd5GEJICBUHHKRA6CBZWd3SJPsTUUQ1JTpQmpMCEI8O6+zs3JMqJyJiW0jR2/ls5LHnsCJIdIEU0dKdmDXr6RUbtwOw0ZuJE1YhMzejYmc1Lvnhcw888NPpI0gCIZCCWlLZ+OATr+4MjQWKmIkZUliOUm2qBECEJRH1Km62CnD+t37V2NhYZw0D4IgIML7Yqm4J9bXB8DybU0REqAWt29B+0wMvNzQ0hBXLUP/bv33SqFGjFiyrmj179m41EMCrq+Wr33/h1kuO+9alkwe49jMcpktuf2Lt2rX1chCAlIwAR0Ssjjpr+G/nbH/22Wfb7BBZhydkLjNTqIKZ73p8fmPBMXpMy9Hu4z1sntMBN4xCbDD8H+d/hR7M3fYiaM87G2LXrrof3fe6Umr8hPG3337xUcWu5fv36MH6Xl1diZ/+9GHbtm+//Ye5sYyvfPfTAzkvKrssQa87b+cHEQ0eiL59+1bXZ7gbLly48PITpzlENlBdvefzzz+3LOsnP7nyppv+o8N2E9HXv/nmghOvniyl6zO6aNFSABMmTMjNgR8LiYi+853vPPfccw0N2e6kwT0oAwYMWL9+fYZnIeO9996bn2oPhUKRSCSRSOga0N8y81NPPXXJ8d+L5UQ9PwdSzEmH7r131uyl9ZZlEXDvvfdOHUdEOB74+tenf+vul7dt2xYUg/3UsvOWoYBix44depU+2AT+Pn3FSpLlrhNkKc3dkoW3d+qSSy46bOpFzz03Z/aCKgQ8R3VMiW984+vHz/h6ZWXtXY/MbmxsDBxeleEGrYDnn39p1apNlmUplbr00kuPO65CEopLB4wd+6Pzr3tq586dZWVl+vr8/Px/+Zd/2bFjx9xVu9P+l8zDhg0bUZ4LoMCp8SMc97iX3/cPhqd/d3dg9dtdAatXrwZGwOu9BQUFxx7bG0BlZdV373g8lUqVlA4ePXq0iOSuX79+S22HLp0Qoqur6/77/+PUp26NhNyj6VIKU6ZM+WzuNsdxXOdX2yaixsbGf/u3B5qaE8FsHHXUUQrYtKnquu8/nkqlSnoPHjVqlIwWrF+/fntdu86PEKKzs/MXv/jVqU99n8LpMaesrIx2dwohtLunbta1a9cuXbo0J6yi0Wh7V1ZgEPrDH/4we3ZBU32Lf33QG9hO2A8//OSUh64mKbw4ZaSU+mDB9vvuu68rMkQIcdKJJ11x+ZQoMHbsrSuvebqmpkY7rP76178ePfqJ3gNAIN18n+yOLV++HIEBYejQof37949wG4BxZV3uYTyMDz764he/+EU8MoyIpkyZst/0HcbkyZM3zv3Ctm23hpUNoLGx8Y47ftXo1bAu3THHHAPgmmuueeiFxXv27Amuk2T2B2zZskX/6sorr5w+fWJ5AZRCnY1rr/1VbV36eY9GoxdddCrp7aWEysrPv3vnk6lUqqR04KhRo2SkYN26dV/UdZJHV1fXffc9dPrTNwvL0rWq62dRTYFfP9pbevjw4f369Ys4rQDG9LWFgB6R6hsyxtVgq8Hdy8ifrm64++6762mwlBIKN9544/nT+lrA8LEDR44cedO//7d//XPPPcfMd108WXprXAMHDty4caMkCUBJqR3Z33333ZDdoce0ZDLp31pHX37yyScvOf77uTlRGAyGQ5uvwA5WSiki7aqYcbuAKaO0OAhwAtJmfN6CNWt2vvnJ1oULF7ZxPln9V1c1PvrsitsuOGpgf+QAACLufucDxfUUBBSoFvjPJ9/ZvrPqrLPPOmFsLwluBcW7unSuJKfJ+HF6HW/v58TqVXJ3qVABkGBmziUx46Rxla9vVY6yrUJmbhf58z/d2ODMiAokgbc+WGFZzqRJRx9ZjgtPHvDsu7sAdMoCArWJ0tl/W3HH5SflShCQUvjo4w8U2dNPPSKqX0KsLCIJPu/0waccecuUKx91HMe2ypnZYZlf2IsUhEQIDKK7b5rxjdNG/sudLwJwZAEz56CxtLjwd3d+c+TIYQJYtaryhv94r6WlpVmWA2gPlbW18bzlVaecNCZCDoHaIVNM33/4g/eXtHSK3oLFjCPUCeOpBDYAG1avAr79qnPvvvvuOJXo9iWCrk1FaasdUMwKEORp5w7Q3N5FKARnvDUtIcnVQT0Nyd0TSPB+G2hoCrQVIgATRpbj8nNP/NOHm5nZCUeZWZFQJKIMEhgUQ78JJfMmT5wzZ06n6gMgjBDSjS8AalR48S8r2+VAQSKC+ra4EqS9DFmF1L9ecPKsWbNKCgZLIMqpcUNL7xjcd/v2EfNvfApAF+cQUYGz5+pvXnXpiblCiDArIopmeF27cwzPgPMCJigifyM/2bq8DChYBOqCcAh/WZjY1V4oJBFRrr0HwJ3XnNVHIg689d4SO5UYefjh//Xgv4ZCIanAfML9Lyx56aWXaiNjALRaFS11PH/F1uMnHaaNgl4Cd1574XFniBtvvLkt3FcpxUIB+GD+u0RUaCVTqRQ7UggRVQ2XnHuyA/xl7pKUnTh85OH/Netay7IsdpiPuf+FxS+//HKtGKuYWmRZSz3mr9xy/KTDcsDMXEj49c9vmfxpy89+9rMmqzczay/jHOeLO//t5otOHQRg4Yaun//85xs6KwDEZSERJewOCpX88vtTx48fn2Tx5z//+TevbyCi9lAZM7fKPss2N31WG+9XHuvFLEg0Arvqknc9+kZneKiEUKBvTp+YD+QBDtRF50z73e9+12T1AVBnl/7nC3OfuHO6RTh6sDXhX89aWCOuXvJOMlzCrCIqxcw3njf+5JNHlerz5JSQQAuwqzZ+x6NvdIaHWqADTL+EcNd3LzjmDHnTTTc3h8oBgEIA5r3/VyllYSiZSqUcJwwgx265aOYpZRIXTB03cuL4b3/75hYqAyBUgsF5hXmONzylgI1VW3Od5unTp3/vGxOEYAfkAC+8X7+rtsOWxcwcVk0ALvzaseX5sNgmok7I199b4iRTow4//A+/ulLvLmU+7oHnl7388su1oZEAWkRJa739/qdbjp90WCkcATpqSHj8tWcvrBFXL3u7kwoBDnMCwM3njp0yZVSpcABIEBGfd9rQkyfcOvU7jyWTyYQs0Gt9efmFWibQp/o1QdY12dfd+2IXhiZlHhENi24//8y+faAA5EGcc2yfP04csXz58iY5DEAdhj388qJTJ40eOaJXb7Ag/tHNM86devj5d73MzClRwOAoN/Yu6f3onefpMW316sqbH3ivpaWlMVwOoFWUtbTxB59unTJ5VC5JLVn4a0reYG+UYIPhkOCreNT3ZjL2CAN1dQ2XXXbzfffdt2DBAl+Ti8fj8+bNu+qqq/7zkT/Fkyn27IYeXSf3nb7jOAs/rn777bfLysquu+4bniaHYGq+FhL4ZfpG+y6R7w+akRpwwgkn+J/om8bj8SVLVjqAozB//nxmnjFjhiVw9tln+6lpoaWzs3P+/I/0trLq6prq6upQKHTiiSe4qZF7VoQASkoKevXqFZRIdX701mmllEUYNWpIQUGBr+IAGD169JjRw0ICUuCoo0b7oTP8ycCmTZv0/m69Q27JkhXvv/++lvqUUieccIKEG8NYAFLQkCFRv5hpxc5LYR+tk0wm0yqs90dW4A5fD1Nu0AO9Jdxx2PEiOSgFOAqOckstgWjU1bP9U7hyc3P11hgBhCwrHA7vPWO8Ywc6OjrIi2b6yiuvrFxZaXsG63HHHQagqKhIl9FxnJCAH+EkeBqWJaUgkkKKbqE8ssuot8Z7OnpWReni24yPP17+0EMPBYt23nnnnX76ZJ3uqlWrAFx//fWRcEgQQhKWRVdeftyoUaN06+ufrF69mgP3lVIcNgIFBQVK6XksCSEGDhz46qsPvPHGw/fff78OyjFhwoQjxhSJwF0sy7IIUkhLWld9e/KoUaM4EGhC38V9QICQJY+fVFJcXOxPeACMHDnyjNMGC0FS0LgxsRtvvNGfduoMn3TSSWecMamib06/8sg113xz5syZ/lf6snXr1rHXagy88MILuv/o9IcMGeI3yuGHDwhW7KJFi9rbO7RHtXdkW9pP3W0+bwODlG5Yru7p897T7+jodLeZKnXYCBQVFeka1mkOHjz4tdfcGh40aJCu4fFj83QBexWhpKQkS1L1IWDTpk2WZV155dnhkLQkEZBKpZ5//vngaBCLxS644AIReLhWr15NRNdff10sGg1JSUQhKa+4/LhRo0YhsJKzZs0aTod0gD86Zi2ahaSbN0FCH2RTUlJYXFycKQBnDK0MPPPMM11dXbqYSqmhQ4da7sY1koAlpR4SlVKO4wZaefrpp5Xn3yyBMWOG6uUdvzL9MS0kMPHI0TfddBMCAxq8+Ov+45b1h8FgOET4p+vBXrDXtDOArwEDrL34XJORSZBwgLfent9pC1CElc0QlkXM3IU8SOykI594d+emLU88+OCNwy0iIESKAWLBzEyKyN9rr++l1VnFYAUBoB5U06h+MuvxLlnw0C0Xl4Rcy8whaRMpWAAcHb+CpA3RRUIADpQQwoJDIOlOHvZ6+rw/xOsVfP2f+cDEw3oPKhH19fVxFIGRQBGAd5ZsHDv5yM2b8fnOZK88+6RjevUGSobERg/I3bZtWxflCqJOUUJEL85dcfyMKWHgb0uqBOO4Y44tC7MESVfMcAQQAggshC2EnXGUF7EgElBCIAJIqOKITDTFKUREFFIJi+M5DEmIgB3lzDjeesSpbrXKAUpRBMCuPY0OwExKcbPEE3/8oDXUF0DY6RTAsP7lIUBC65QAIBmSU64B5AVUBiDBAoo5XUv6NAEFd5k/PzeHG9P9xCHJ4Ka2rv4o1nu6wYqZW2G9/fbHu+yyUCiUpAiAnPxCpZTWhkPMtm2Hu7aHw+GLTx9pCSEAC47QR99JOMrxjuDSe8/hMApjUYsdkRFjznWuYHAqBYAcCitWbVZZa4IvufutqVM7vnfJMf0qMKgArz1134BeeXlASEgAEpBKESt4jg1SpSy2QwpSQFvcWfHsgsaBIJHyJnvp+QxbAFKguMJ37vkvAOt3qrq6OoekEuHedmUsFrv7yjNnzhybDwVGFOLmy875bOvnp4wkvYpiAyk7FbVCJx49fOHnG8BIiSiAXXVNtpeHCMDgEKgwxrIJpEBIARjbL6ciir6shoyPnT7r4qeffvqCc8/pCzBwy6Uzq6q2njaSwlqb5JRSKiajkyccvmDTBgBJygGwq67FdtcEEIYDIAKZF0p5z4vDzDlIRYB8TukVpBmTyvJVCxTscK5SKoREa8OumBpp6VDSzOeecuSCt57vkH0ApFSIiPY0tvmnI9YpvP7B2i45GkBYNYPAeaFWz2+4pAxK2pJsAEnKTdpYsvqLKZPH5GpTWEHq2HDMYEcAFjshN3ICHFbMXK/kG/Pd9C3VzGDOC7V7p5Z1T3/Rqq1TJo+JcVIKlgjl5UB26DUBG4xx/aPlEZSDh43PnTbrwqeffvqb3/h6GRDiFAQsyDBSwUGG3TgrJEgUA8OK6LSjRo0sgj4spBH08pyl21oEZCnYJqJC3nXpuZcOKiS9ugKggHDjZV/bsnnrGYdDu8UkdStQ6ISJwxZ8ptsuCsKO2taEXq6ACgEO6Wfcce1O1mtfjgRyBAgslC2EiEOwcsJkW5yESs+BlVJSSr3HdHccc+avT1ijAAhOMqvB/fqEoGMGIwQw84kTwjGua5BDwEhwEYBP1q2vTaBXhASJEADlFEVlV1NCSBDBcrpCKp7jqJAUETCDzzpR/Abb20SF4zgpkiDsqmtK6YjgRHok9559HV3xy5zkbDAY/r/jn24HB4eSbAFjL6PMxIkTfz7mzEQisXpL0wcffPBFfZKZdbgD7Uu3cePG3/3unf+4YUZQiSAiaOXGW3bvPrNnZofo/vvv7+zsnDFjxrhxY/zAaDZni456fLcdW0rJrHeDfLkzmQCAACnl5MmT58yZE8zb4sWLHb547tx5AM444wx9jCkLMXPmzMcee0yXV1fUpk2bvtiB4QPw8ccfAzjllFP09MKXWFzHA0/s1On7sZBB1D1jwfmJfw0RCUJOTo5vSCultFqjtZ+6BlRWViJylJ9gJBJBJkJoRU30qP8Gbb6MHgKUlJTwdn2eqvKvrK+vByr863XB//jHP67dYzGztuQc7R0IJiJLKQB5zu7y8vJLzrive6MEBdru3/r6JXkxOgBo11//qD995fvvv79u3tOnnHLKdRedPnBg7+heThHOFMN6roesECUMFiAFam9v9yeQvjZJhGg0+tFHH7XJ/gBEOHTYYYfNPPkb06aNGRxGiCDc9zqmTp18CrnasFYihRD7iAfoH1tIGbVN3jQGggQIuTm4+aZ/DXn+J6eddvJpp50sAoa7ECLlZJcxq6r9ULLZn/t1BUgpc3Jy4vG4vwign1MtlhORbhcKhHlpbW31te3Vq1tTqRSH0nn4xS8ejUajUae9q6trQ30eANu2iUgPhU888cSUyY94OYHfNPAsa7+FdYJr17amUikOpx1eDyR9QULvgPQrx3f40eUCISeHbrzxGl9+RuZ9/bpizzOKCA8++Ev37HOQAjc24ZVXXgGGa3UWQEFBwXnnHasXj7xOhalTp5x26hThHbSnO4Dj3StLKNWtptwMd29Vb2Et88RH/1t/QAt8i8rK1lQqpSPn6Iuj0SgDyp0xEYCcCPr27bujOV1jSqk1axpHHVMcrA2/foIrUbpQAu6YxuyeVKfHtIzMm1MzDIZDj6/CPxiAp5vagZtmnEDLpBUFygUfP37A0SyAvIsm93YuHfa719f8/ve/r8FwIopTDgjNsvyP89beeukZvUqsCAMAk2uaCBLQ68lQgUEQYMSJbEW/fbVy4YaOVqti+Zb2q37yVn5+fgg2gKTM27ZtW1JGmRkkmbnV6vvS3E+XrtqYSCSi8W2JROLNx28PWyLPHS73WmBtNBOISfl2eRiSganHjZ37+nMEBtiWYWbeGR/69mL87cOPlYx8/bRjcoEwO0R0zqn9X3qysQkDiWDDAtBKA1/568ZvXTRyyWfNxaH45OPLwlAE0vEn0i8bFoKVYCWhj22zmR2VIQ7rPViO9y6GZLYA4Z6+lgRxCLJXYV5VoxacREhaiVRCQZ9WhRXrW5LUy0EI+nSu9KtOkVaFGaEwFHunYXO6mykSgLCEbhPhvpCZBTEzJGFgRR9r1U4CJSiXiGwVJaI1G7ecePJ499QwAEAMeGzWTzfvatq2bdsrCzavW7euUY4gfV4aczFqY7HYnZdPmzhxYti1mXQmlRAEm4QQIEWC2bUO3XrxFnx9o8AGdFwsVOSLU44e+Pqq2ihRO/UCiXZZTETbVey/F9S/vuB3Z5999r3fOSoiUZTRHdzoKFnvV/98XfcBQKaru2IiYjduiWByJzMKKRJksR0hPHr3t5uu/0ZjWyIajZaWRCzLknCInIh3wrb/ZOnYyR2AYtR10fLlO/60sGrp0mUOFRGRZNtxnMyVDQWQgCD3pAgWpIhggQUQZpuIIn68F/dUQk6AbKANQjHqu8Ty5Vte/3DLsmXLlCghIoGkTshXagUIzJIgPW0SelXA9ft3XB8bcEEskuho1dNRkBNPdJBikpBQAIryRIjcSZqyCZ59pkgwsKumVXFMQRLI4RwA7y/ZCSBf7QEAxPOAfvltvXr1KiuI9evXb/KEk/NcNZeFf6Rf1i5eCDAUkcP0RU1jkqIKEgzJOQDmL64BkMc6/U4//T6FsX79+p10xMl53mxBAGn7Gu4iPwEWEgAiHILr2sGBxzc9g+X0mdVCEEWVoxfDALSDbKYHfr+wOl5qyyiBcu1aADdcOLl/BBanx0btFJQk2EAXwIyGuFy2bNucj7YuXbrMRgkIxAkigtBPh/RPrBYSIKUPtMuY0bmrPYLZtUBd7ytBKtOZTavyW6obUhTT58DpnW9WKET6HD4woBgsIPr2LqaGpPTOAbWRu6chrmfvApBETEq/CIhIQbEXbpygGEywCvMLqIm1XszM8WRS+b2d9eOfbnHjH2EwHCL8r4gXEURbjZJ0TFwIS15w/pHR6A0/emIue7FRmdm27dWrV089dSICQiYAxUpA9minMrB7d91zzz2nx+uqqqpq1QXA4hSAhMgFwCLDw2zPnj1tNa3MXMi7R48enSlh7meY7H4MG4Dx42PRaJRVxsvg4YefDHV0DBo0aPjwgcJ7wxXm4fjjj69a1OlKkkoBmDt3bkXvCIBJkyZFQxlZ8PUz/+UdFH2DPnx+hpRS7gHWmbqUH8U5XRbPqU5fXVtbG/yWiBoaGoByX0dkwLICMnMgh+jmEetrrlodHDNmDP6yU7e1Lx1t2LCBca4r5Hgv3D69S/J7l0w4YtjR50y95ppf1e927Ve93/zUU08966xjLAnfIkxXC7lN7B4J5uVM/+n7TAczzGACbr311g+//4fGxsbuIhmD33zzTVX15/vvv4etjL7R49KEL7tmX+l1m2y72WsCjQCYUVycX1iSDy9GGCHb4VgxlIJNqKmp/3DRikWLFi1Zv5OZIwOOraioaNzRqZSSlJai09pkZu/NEKpd2ZKgT5IjMOAoJwWrpqZ+/qIVixcvXrp+h1IqWjGpoqKiflci+ORSug+4Jm8sFvMbWreLQFq0IyLLsrLyAG9iw8xEwpeKKeDLywAz2trasprgjjvuqKgoKo9xLEaxGKJR16KPMgS5kQF1DvWqAFG2nee1FAB0dHRktewdd9xRUVHQNxexGILpRwAB5AVzLgKbQf1297xogv3Bvx7o9lh5dSVIEJHtOESkhNywoWrevHnBjlRWVnbOOdNEoBH15w6rlBK7d9d/sHjlxx9/vHzDTsdxcvodV1FRUVvdpVNwu0fghwzWUbyDdQJkjRwZXwXvG9SDu7q6/LbTf3R0aC9tTndFoKioiLkz+HNd+dTtplnPr0/QidxgMBg0X4EdHNSZ9nY74Wtigh1/DVHHgi0DXXHWYbPnLq2qqqqXuQAcEU0pVdPS0QXkMaSwdJzUFAGA1v9c/0vWGgkUGIy1azaG7CYiKmF9wlF6XHZY5/N3vOsAACAASURBVLNAZ4KZJRxmjnAHgIjTMay8OMfbt6/fvoH8pzVvrXEKfe4XA1pBYQbZApRHcsrRw7ct7hRCOJTLzHGKplKiXKW+dsbUQlISwmIBoBC44IyJbyx8A0BLqC+EaKHi1iQ99tQLjNjMkybkpc9Xg18WEDitMOlPJVg6jiO8PVv6UyVICdJWiQKxVxwJgiJIRCIRN1yB5+KtSyQIDtiXgW1EwNi4tXryKWMBS0t7TIo8CROeWu/AO6mY0r0iIItKBhcAJx0xII8bAMQpF4DDEUFi1Wdf7GlHWZ4gIEJg5hA7Qght5eaDjh83aOWu3SBiYUHIEHf1652fJ1iCiG0iUpDKe/259gT0Ad5EgANmYikEAY42JvW50O5vBAFFQEHv8Jxff+fHP/7xJ1sTRNQqezOzTREQWqkcwPtrN73w7rbrzxrkv2mZHb06Ac+1QLhHxumk0/9ow42EFoOVVhz/H3tvHmZVcbyPv9Xn3GVmmJ1lWGSGRWUxCuIGahRNhImaH7iDKAmSGAGJAiafRPGbkEhiSDAJ7oDERDCK4o4iCCgoGpVFUZaALAPMBrMyc7dzun5/9D095y4zDKAI4b4PD8+ds3RXV1VXV9ep7paQ6phYHVlnBrFFgIkwgUzXbtzkfOmGc+JaBeHVN9c+sXjzrl27fLKBiH489ILi4uKCUzs+//zKzc98BgAwAdZL+pJ751FNJgBqd22COu/XJKASeG3x548u3rRr1y4vNxDRmKEXFBcXdzqt03PPvbvxmQ2KQmYWHO3qEkwUXcfGsesHENVpUspmEKnIPUj3MqHiyjpqD0TXAFiuSKcqQ53BztGkDhvAWV2pR3e0BwMsGAQS0XOh1SrPaLuckqObdkAYTVH76OlsBkiqs4KdvmMDOLOQe3SH2usgvnwiwzn711GPGOfe+dMgUpnrULFU26maOEZxwIIAySzI0THDtFnuB37zyKuNnnbMLKQFINvaO+XHP84VUkCtXoueZrcfeGXx+ife2Lpr1y6/rBdC/HjIwOLi4o6nd3zuuRUb//kZgZxDa1QrpEq5gSQCWJIEk1CHH0bPtVbZDCL6pwHHVaXoDo+kjINqshDIyswU2gmGIFB9fYCjsla+twTBUYDo5MEAt0nzx4TK3cxkG2xDCHZ4SIS0dD9qIKUkNB3Eo7W86dVon0oFhFNI4aTA8RcPjrFH0UiAaRhXXHHFY489pm8JIXROKoMlyJbyy6176uvr+/XulpGR7h7KVbIaAeeff/6cORcTwSIIASmbnqkOY+nS9f98bY07lnDeeef97MYL6+utLKu0a9dThBOPaTmiwM5GXu74sd7P4eKLL/7nmiXuWKxqzmWX9RLRABhUq88664zs7JW1tbX6GfXD4/Gcc84ASkYJx8Zf1W9dAvSo7gQXY7jNDB1U01eaPGD1YRMMZGdnqysqb9IwjA8//HDMj4e6pRaJRGlmJz84Lg7URE801hXd8iI7O3PQoEEqBxoOqVLKN974/Iwb+gpqily6CU5LS4vuw+rOyaam/OmYi7FZj0ykwk76O4D7Yeg4MWBZdts8869//cO8t8qffvrp2pCLEufh995772fFowynsbrVMSF5p4Y4CXLMlKNprhU3ZitihBAEsm1bubBqOqE5wETM2LE78re//a3S38fj8SCEQYMGjR9f7PGgwUmO5JhjTWI+C7ipEq6SyWFXdNUi46ud4b/97W9l/t6GYZBNgwYNmjix2CA0JNvog5uURMku5hkhBCAZEBTdANidq90UjIyGJ8ktIB1UhlMFgdq2bevmGzOvWrWqZ/erov0xRied3Fw9KXc4k6i6qo8LEodbfjQjggQ7H941w1WjtGSdGWOTAsQ8T8TMdXV1EogeU0dNQfTX3tiyfft2eNUnLGLmPn36XHxxp6hfCdi2HbHZMIzteyJ///vfyz19VBUDBw68885i08BBoLGxMd4+NMkuBu5oruuZaHg4ztcHUFtbK6WEs3l227Ztm47wtCGE2LdvHzWVAyGEojlO3F26dNHeamItUKvxhJ5yO8nuRIYw4tYAtPw7hRRS+B/GMdsiUQJSIpmtAsJABBQmCrs/NDNLSAnpAbyM3j06mrLRBAxmnzzoteu/07OjH7CEiIA21+OmSY8Nmfzv637zxkW3/PGpFZX1oAao+IsQIGJkAJ1zzYs6yws7yUs6yu8W8CUd5aWd+OJO8rud+YJuOLtz2CODXg4pGjLk/sI8fP9U+4ZzzeHndzmnk8qgPYR9JDgrVWJ5a8AQEGnAxQOKMtGQbtcpi+yzQ2a48YLvFJySCy/Dw9GMwDQgy7SLLz4rXdYYJA2SNgubhZ/rzj+zW7s0meFsl6bayNDb6OqhVwKSiZmYo5EzIRkSsBkHaw8KbqJQ+yK6deo0jab1ZIIkYEAabPfq1t4j6wVJsGUbPtvwrd9trfgMtaB6dbIdRNCCFHqWJQBR29BgQ0WS48+sFkIQJNj2AWmECTcPaRMp9dpBnwzZRBHmOrPgyRdW7g2KGsAmkkQSLMEGbAGLAEFMbBsUdVkMkODoGnaVKCgAASZY4AiYwGTDV1nVECYRBmzAhrDYlORlggSzMCQJ5VVFSARBf5yz+ML/72eBBqR7MOmqDmuennjviL6nezb5ZLWfayTDlmwhrbTyoM2sAoZgSZAU/TjAtvAELI4IWIAFYUGouhzVUefBRfOVo3MH6Ci+2jCZAASDQd2nBDFBOomNQm3UTUQ2EIFctX5Xg8gPitxGZDNE+w4dvR5pqri9bQEy+rEE0peWIQEmg2EwSAX3JKncS63PQkVBBWAyGQwbsIB3P99dZ+Y3iLx65ERrAXtVtoZlq7YQgxhpGZnS0Qf1TwISktW2HWQ4eRMxYWn31E5zQ5ViqLCxECxBEAaIZDTn2wOY4H69svyoN9g22LaQZlP6P1/6dGs5QqCwk5ZsQURAISAABFhEEA39qv0QlLYoMgKRiAWATZaGYmO/03P83GBwxOBICP6ISP/HS59ubr78IFOERFzrmJWKRh3HqD5G95xoarUN2AIsuGkKHV0XYQmwhGAyaoGddXLG029XmadYwmMJT4Zd0UZW/nLMFe1UfJpRATyyaNXZV9++eEPlexv21It2B828g2ZehDi/YwcPSY8KkFg2yGJEiAHJGemZBBgw1FckEirKGw3AW2zYMEO2bUc5J5psIMfkBVlS2sxqLCBIwbJ/77Q0qjFk2OSIDZ8k/9Zd1QclbJDtmo8dbIgua/Mh5JGBLLPmnD5ksDTAFmCB6+sbgKgmMERjIERRfRLKQgYagyqd2EmIktKZ+Sh7wmwDEmyTk9KcQgop/M/j2G0Vzs6KsUT7wsD6DRsDwbB0xSB1DFUZQtOMjgcq+tWxY8cePbqq+I0t7QcffHrTpk2GYRiGEQgEZsyYUVFR666GiFRU2BBC/RNEphCGs5JaOGGGmJASs2EYBiCiO/QeOkKQ1Ho2ha+AjIz0Pn36xFwnuvzyywkxswQCTGEUF1+s01U1Bg8ebBhJFtqrtzTx+kV37FNzQ0dWFCKRiL5i29Eje1UgtuktJzjXowdlZGTod9XFmTNnVVXXRzialylltCgdd1G7pMat43YT5nAJ3YvEtddeq27p3XMbGxsfffQJC5CSJUebILmpjTHhKBfNuq60tLQ4hnz00UfhSMRiloAtpQqcu7eDkA71zGxZViQSWbFiJQEGoU2G/5aRl86d+/Bpp50GV8S3Xbt2ToQPfr9fc1WVVFpaqsqUkDbbQHKlUqXFRfe1LMLhsJKFLW3Nz6SKV1lZqV5U/sm2bduklMrzjYmIMwshZDReGxOFjdclQMlU18iMyspKXY6Uctu2bUoS7iJckcLYTxaxU18tNWaWLMmlq1HfhaMBQsUBy7ZUHNR9101zu3x07txZdzQAoVBoypSZm7fstCQkS5ulhLSlHbbwyqtLpk9/UH8M8fsRR1h5eTnr7wMMZujy2Vm11mz5EbzyytvTpz+oyrRlTP5SXC+TToKHu5fp1F63dKK8cvw2CTz++ON6P2N18cILLzzjjNMMJxTREMSCBQs8Hk+fPu0rKytVgUqmWnaKOD0NjorAJTUAaWnxH5TKysqkyhrXlEd3Ym46L1D9Sc7HDiLKzkS/fv3c9ioSiXzxRVg6PFFKW1ZWFo2LMxPR4MGDPWZMSr1b9ETktmns2oREKwMzExAnCF1UKh6cQgonCY6RHywliAyCZGmpOTdDSthhIMS8amvg1nse+uvCNfuARhghFd0EGEQwJGAT9pbWSpheGfLYwVy7dMqPfpAFmQmEgPK68DvrtteZeXWUWUeZtUa7GtH23Y82BpviysQQKiLlgfRAmmybbJuQJqSPba9tmYBHGz7JkGywNFiqLS5NSMP5px5xInaajfoTNjdlczogIhXv9ANpkBefc5qfawVsg2QGV2QbVUMGnd4GIBEhstRhYSbYBJ/ZDb1P8Zl2o0cGDI4IGc4Sdd8f2N4P9jlZd3r8k4AF2CSCFlnwOFEQDwmfFIYNhCm6MNwShg2PHd0wFJYwI2SEiMNAhMwwGRZgCVOvmwcLC8JCNDM4h3DD0LPbcG2GrAGTtLnG6Lxlvzls8r+e/Y9dYmMvY9ma6oOcr+KdNjwWm6U13MAIQYRjdcMZb4QQpgfwAAXAfbddcdlp/vbhbX4O+GRjg9GmwWjz9PLyvy3au19QFSEIEYYRIW9EeG3gYMiyYdgwBIhtaZNhkxFRx6+RATLSgIIcs41h+WQAxAxZbbZbsXHvL578aEUJfViGJ5eUvLz040bKJRaQFEFGSKZtqZBVjCAQIWrT7rQg5T3y3HuldlTkecCpufjD/w3LlOU+CnsR8uPARed08xJUsDDda5sy6OGQ4ludkf/SyvVf1IgSRhWL/+63IpIsqNxWGZUXM9TnESAIBEkcaAiGhUdK1nHZWtvbCATJFxJ+G6YN02GjSiZu0r38vFwCS4INrjfavv/FvoUr91QCtYz9YZ8NIcmwySuFb0tZsBaoh2gAhUEhIguoqg8qvVLfHGpD4TAQEaYlPBIModbpc9vcXMEsmARTQLT78PN9z64sKQXqGPtDPpZesM8SpiXM7ZXhOqCRjABEBBQGRYBGK0q/hJAQQZsswIIhyRMBLBL1wYhFpgpdSxgSRkRtrixMKTwNEiGRFnVxDI9NhoRXOhG+TODW4QMzZXmmLFdSqDJO2VqdNvSuhdfdt2zqgq/+uqT2wZf23Dnrw8GjfnPvo2/mFHQOASpS7vWBRPS7R5jSwpT23PIt/63FXkF7CNsrag8SMoBbhl+YKSuzeL9yT2uMU7ZVpw29+7nrpr4dU/4t0+597K2cDqeEAIYp2YgABw4GwEQQ6kp90IoAEXgt8tnskfCoGLkFw4IhgYglVeTYZpIQFpkRwCbTghEi0cD03lYsfHdXgHIClGOwFNLOMg/c9dPLPQQC6g1UEWbOXlkRyD63f580D9rnZpgcEgzBaED+h5/ve/bdfXuBOmB/2GNJQ8Kj7MO2ynAtUC88derITGa/AUPWgy2wZVGaLdKfX7F1Wy1KydgLsbWioY5EGAhDBCyyyGtFP+Z4bXhsIhsqVZmygfE3X5FplWXZ5YawwKGDZvv5b6w5wFQDHAQ1CqOiAbsqLAPSJM6z9nb2VP181FB1WCUz24BNRhhmhDwqfG+Thw2fqsUmYZP6YiYksQ0JMhgC5OXo96t4lzduwpxCCin8D+OY5gfribiUEsKJ9RJ16JAGYMGCBb169bpqQHsZ9ZsBcnJ/bfu9996DEwsZPHjwpd/tql149QU/znJFIpGYqqGWkoA5Jt2w6TegzqB3x4OPDConr7m7gsSAAWfMm9eU73veeee1aZMe9WgpSqGTR0hDhw5dPe+TaA4oUd8+fdukxQQRyVm+zWCLEQ6H6+rq4Nqqs7y83LJsMg1Vpg3U1NS7EyGIqKamxrZtMojBhmFYjNraWqKOcFpSXV1tsRPNBW666aZ5Kx+rra1Vtau03fLyst/+9rd/lTu9Xm+5Vaj3D1YfUt944419+/ad7i29795fItpY0jS4I3YEeD3mH//4k9///tlX15arh1UMae7cub7yottvv1WqrcUAADZRWVmZYpHNbJomIvj4448j/1+h6fM66YxkmuYZZ5yxfv16d7xnyZIlixcv9sgQEUW8hZpvzGwYxq9//evc3Nybr+g76uYb27VrJ6Wsqqr6/e8XPPyL67xer0pjUIfQKc3s0KHD1VcP0OLPzs7Mysqqr6/XGlVWVjZ27M979uwpq3bs2bPng9efEyKaa8uuXGGVrGwLYVn2unXrouFhZiGImffu3bunFJ3aSyGEZVt6O4VohiiauHrOOb2efNKJVVsRIvrzn/+8aFH3/PYFn376KVGR5sOGDRumTVt4ir3t5z//ece8NGauDyAYDKq5kto1uaSkJGiBBTMxu7JmY2uBYRh//vOfX3qpR15bVUtP/eS6deumTfN3sbdPnDixc16GlFwboJqaGqYu+lNMZWVlJJpWAtUrGxsb2ZUBXF1dHYlEhNfDgJSyulroMKHSor1794YtQETTh4cO/d5TK6u3bNnCzgcHIjIMY8OGDZvWrQYgyJZS+tFw7rnnjhp1s+5ZuZnIysraG2jq0eXl5WPG/O6cbqK+vr5u95a33nxJAEOGXDZv5YHNmzdHpeAuf+0HrvID55577i23jNCdriGoU1yiUigpKQlZUHulN1kMtckfZNiG7nGaHrUVtMr7D4fDDz30pOopUkYPM7/q6qvatoUKMIdtWrZs1RtvvMHMl19+ucfAueeeMXt29JOCIQxmnjFjxqJFPdq2z1+7di3QLdrBhVi/fv20af5T7G0TJ07smt8GQF42cnJyShpZ86esrOzHP552bnejvr6+dvd/Fy9+kQiRSKSmpgZoOl6xoqLCsiwyPUqgAvjOd06/5pprXnrpJW0N3nvvvQ1DCvr1O91DkFK+/voXSnwqS/iXv/xlhw5+4eSR20B1dZ3ed4KIYOPAgQOWZQmPR8klwqirq1MWTBCYubq6OiKj3KPYvOdUSDiFFE4eGL/5zW++4SrUuXHKrADOOESkUlalgLT9YvV7H20/2Pbt1Z/5cvI7nNpOEEKECCgE2g8s/2jv488sluTJs7Z9//wef/3lTTmG4QUbIBtIT/O8uPgjq7HeFn4PIul2lYeDP7/1svZtc/22jA6STq4eKJqaR0TRjYYJDBwksW791vc/2ytgg8gkOyNSeVqPouILuhtEphN7BcBOtgYla2P023J0CI75TSp/E0jPo1dffKnBzvRwOCey5ye3XH9210yTyCSDYBAJEBkUXWPTuW32v19e6pMNfo747IMTR1zYt2fXDEhy0ijUONgIUVkbnvP65r/NeWn3/hAJj5CWySEyqbq+4fP1m8tqzdNO61QfFv94eePDsxeV1UsSpkGWkAFBXFVdV7JjT4BzOhflrfho859mf7h9dxUZpomIAduAVVUbemf1jkhjbWGv7nkC2T70Ov20j955Mcw+LwfAZMA2yJYsTLuhoGPXUaNGrvt4tTTSPRQxOSA45LEb9+4ty/E0Fv9giN9hiuZh1JUBAJiAAc4x6OqL+3bKMUrWvRW2hF8etIwMSebazQeef2tjTdgfSusY8dD+ML3xftVzr7xjsU9AemSjKYNpsqayfO+p7dNPPbWbDwQmg2ACuW0L3luyKCza+CkMJgH2sgSbWdg/9tabzjrrzC/WvU/C60XYJxs8iPgi+yHlDcUXFnU9pS7sWfLmqw1m569KDrz10eYGX5cAZ35ZEnro0RfKK6ozI3sL841/TLu9e34bH6RJJAATsjFCX65dbYsMP4UlDCLRaHt2VgTD1Xu+O/iK4ovPJCLTyQRVn0AkuIaM15f951/L9856+s21m/dI4TNheRE2ZchDEYu9L7/x/rot+7ZWcKCmolNhR390LkcgwSQZbIBMICubdn65btu+Wj8CLDyC2ARVHGisq9p7as+epeXVXooIEiYsgHbuqazavWnUj26qrKfnX139p3mrSqtCTOwREUOGTURqGuw3V2zetbt0R2XkjFMLJAkPWABZ2bTji7U7yuo8sl6SxxCAbVdWNdbs333aaT3LKg94KSiITFgGaNeeyv27N906ZkRdBE89u+TBOSurD0oQPAh7yDI4XF8f2ri1XkbqO/bosuCNzY8//eq2sgZheEhKk8MkqKq2YcMn6xplVpfT2j/72vsPPbWiui5MJD0IexA2ENlZUrlmfVmgtqJ73565hDTi7w7s+8UHrx6oDfv5oE0eU0hJYCGYhEVmOjeYxBOu7X/fXTfnemwfkZdZMDMJhA5+su6/6dzAwmtwhEVa0PJUllVV1uLK7/b+7sWDMoEMwRed32fTmtdU+VJ4DZKSiMlgMizyZHC9SXLCtWdNvevmTI/0Eu2oo2dfWf2n2asqq20SwkOWQWGBcE0DL16+eefuvTv3W2ee1oEBASFZ7Kqnf7/8wR+eWr2nMmgSeSliwDJh1dXbazdW1x8oK+zVgwWef2nlSyu2SviY2GfYBod9htyxo2Ll6t2LV216+z8lT7+yfsHrH3uNUHq6f/pdw9uYZrts7Phyw67SBp/dAOERAEl7f3VjbfXenj17llfUmBwSZHggDUm7S/ZXl3x+649vzjbYIDCIQvUfr/9vmjzI5PHAtskTsj1lZRWVdfL7gwf2OvPsf7722cNzX9q9PwwyPWSbMghQdV3DF5//t7TGOKNnAUxKY+kHXzSgE9Xs2Lh5h8+qiZhZTObKDz4Le7uGPTlvr9n9xD+eC9six9rXOdN+/Nc3XXFu53RIL1EQVHMwMvu1Lx+eu6isVhqGB2wbHDKIqmsPlny1O8DZHbvlLv9o84NPrtlWUi2JTISJQwLh6trgO+9/FW6s79qrWyaRSn9zvrFpu51CCin8jyN5WuHXCvc5GupciahHzNE/UAZj/rPL/vDPD4goC6Xdu3f/8ZUD+/btm5Hura6ufnnlZy+++GIVTvF6veNGDh5x0wWdVMgQDCAAsoEX/hP43e9+Vy7ziaiNvf/KK6/80/gLTUKO2qEoemqDWnzDDNZLr/X2nBVkPvPsst8/s5aIguQjovzQjuLi4ocmfs8U5ItGv1SkNhoTbaGNzrkJMb+hdihj3kfi/0174oWP6gF0xK5XX53V3VDfvOO3rKoDLIkb/+/Zzz//3GafYRir5v8qP8fIVyVH6RcAGkDrNmwe9cvHAdQaHQDY8AEwKAAgL1IK4K3n/3bw4MFrx/w/ADVmBwAWeQB4EGTmvPDurl27zp879c6J972/zWDmRiOLiFgSM/tRD6DA+u+LL84vzASAcmDLfwO/eOilnTt32pQOwMv1nTt3vv2Giy6//LTyWtx885Qa0QmAgUiHDh0G9+vy3e9+9/u9vX4vshO0hF37VCio9IkDQN1BLFhRunTp0i++KiUi0w4BSJfVAHyoAlAnOgFoQD4RCVgFBQUDT88799xzr7+gQ3oa2gDMUCnPpcDbSzdNfXJZY2OjDR8zp3P9eeedd8/YoaecgkVLK2fOnFkv2jFzYUHOgAEDrujf9txzu3TyQAA7GjFlyp8+2RkCYOIgAD+HAQQoIz09ffiQC0ePvuRUDwTB79DPwB7gscfenLv4CymlBYOIPNyQlpZ217UDRo4Y3AFSQJiwOPrFQi1JktVkXnn1yBK7q5QyYOQCsMkDtdjHkVearCaiPtn1z//7yewY7qlTAwwGHwBV1chR9z+7ffv2ELUB0NbTWFxcfMuNAxobcfO4GeFwOEIZAPyo7969+x9/9v0+fdp/sGbn/fffX+HpDiBoZLKzLsuLBmbOs/b6/f41L//dFOQFM3MViaoaOeq+BV999VVYZDJzvrdx6NCht9x4TiCAET970LbtMPwA0mVdjx49fv+zK87o22Hrtvo77rijwtMTQEC00R8E1BLS07Mr589/9PIf3mVZVpXRkYgi7AdgIgCgrbU7Jyfn+X8/ePXVo8u5JzOHjQxmVjuU+eVBAAXWliVLnssXIKAKCITx4DMbX3/99aqwX/dZkyNE9INBfUeNuvL8U2AQ0iEJ5GEJoIqMoI3fzFn7yiuvHBS5AISMEFGBr+HGG2+8+4ZeBkGdwNx8+RYRXTmo16hRV6vyVTrTGx/umjp1aoUZbTsAEiyl9HIjgHx7T5TDRF6AgaXrKu+5555SbxGAoEwzDMMGE1GaPEhE7cObFy16NjcTEyb+ZsU2g5kDRnpUFaRM44iWnQU/EflQM2zYsN/f/l0TCABVNfKW+57fvn27Si/J8zQUFxffMmLAwYO4ZdyMSCQSRjqAdLuuR48eD/zskjPP6KL2yKgGwhbun7v2lVdeaaBcAAJhAO399TfeeONPrzvjyy/3/XTyH4ioRhQQkdJhtV97rrWXmVc9PyM3J90PZuYGEpbEkk345z+fXfP5biJKlzX6+RCldejQ4friC3/4w37d/TAIfoDBQdCu3WXFP50OoM7oCCBCPmb2IsjM7a2Szp07PzPv/okT73v/vwJAg5GlLUyarAPQPrxt0aIFp2WBGWC1bUU0OzrlB6eQwsmAY+AHN4cm37ERCEUiP/6/R7/44os6kQ+ALCYiMkJElGbVpKenX3VJv1tvvfWUHGEIw8sg52hQSSRZBkiUlOx7652PmPnsM04dMOA76SSJyIx+6orf6zdmqQQRM1eTMX/+K3/75+sAPD4rIyMj35RDhgyZMPKHhiDh8oMdxJyHh2YyIdypF/pKgOjLTf/dtL3cNM0svxx86SC/8oAT9kGzAGb8Z+NXX375pTfNn5mZ+YNLzjWE8MTmrhFRBAhHIhW1IWYmj2mapmq7LcDMdkOYiNrlt5FSVtYc9Hg8bKr9TU0islkEg0EPR0zTzM/JrK6uDUjh9/vVmXyKRTZTIBDwk9WubZ6fbWa2hJBSNgpz69avtm8vFUL07tm1sLCzz5BEFLDpz3/+a2bbzj16nuWgCQAAIABJREFU9Ojfp1vbtrke2ETkYRZCmFFuNemAGnkMtcjOdRSI2sY3RIKZK2vtzz77bMuuPeXl5TU1VZFIJBIKtmnTJrNNVk5OTmHnzkVFRUWF+VlZmT62iUhFSQWaZkEWyLbt/QH67LPPKg/UtWnTpv93euTlZftIElFJ5cG33nqra2GPPn36dMj1CCEMdUozs+KwbdvVAWzfvn1f6YGGhoZQJNimTZvTu3ft2bNHhscSJEzJhoj5qh0EbJvLakMbNmwo3V/j9/u7dsw766zvtDEtQcKM7sAQ0y+YOUBG5f4q5TuaJhmGYQoAIAYz2ySCwaA6P6KNIfNys7wuPzj6VQJqHZdhS66X9OGHH/93x75OnTpdfEHftLQ0UxAz79xXt2rVqpr6YH5+/nn9uhcWFqZTWAgRjJg1NTUhKTIyMkA29KmBkizL4kjYNM2OuT4BYUCoPhhXy0Xn98nIyBBgADv21q5evbq6viE/P//8fr2LijqnwRJChG1RVVVje3wej8dQuUAEAIFGC4BphnNzc6r31wGQHsM0TQ8EM1ugQCDgIRiGkZeTUVFxQJp+wzB80Vx3wczBsC2l9FOkbX6ux1lKqDytUCi89vOdu3btqjl40O/3F3Vqe84552T7pBDCG90dTxLIgAAQYQkgyKK6uvY/n22vrKxM83s6dux4Xr+epmmmQ33TV+eZEzMaBA5ZvhJSnU01NTVhGH6/31AHLAsGIG1h27YdCno8no45GWq7EAYfjFB1dXWITI/H4zVJP8/SCAaDPkTats0zgXA4XB+RACQbkUgkbCsTR8FgUJ3tLC27oaEhO9Po1q1rBktBwgJsKRtZrFnz8TYluwv6pqenCzAR7dhbvWrVqqrGUG5u7qCzorIzhKH6rw1IySFQdXXtfz7bVllZmeb3d+rUacBZ3Twejylh23ZlXQCAQcLj8ahEMAlm5nAwAqBzXoYQZKp1mSQABIlt266ss7/44os9+yoaGxvT/P68vLyz+nTv2LG9n6QQwnSd++3YvYBpmvCQEELtr8xMoVDIYGkYRm5uRl3dwVCEVS6TlFIIU+lSY2NjOrht25w051ulthWiFauiU0ghhf8BHGs/2JUE2eQDBQAJlANr1x1477OS7du3B+sbDcNok2UWFRWdfWrHAQNOyfdCAJlqIwXHD0b0vAAOgSSgMoLVORo+MIGMmEitWjZnu3dy1c2vJWPPvrqSg96CAn9mFgjIkBACuWqx3pH6wXD53PpKAGBnH3sVH/I3U4I6DUE9r3Im0xUHYjgJKI/ZKTPqU7pK8LquqCfV/wZAzvM+V6vCat26iyrp5i3bACxScxjSHDABctViO+X4AXZS0T1RucBVaqwf7CTqqdap6yqJMgJIJ86qKtXnLKt6hZOvomrxucZLNbZJkM0ccC1+9zj/ExACbKdAN8cMlcEDkg4lroyUaLtUXR6HGC2dMCAduUScuhjwR/0tSuoHB8lgIOLKthGOvKSThaOa74tqexM90XwZUrngQnFPbXBGTru0lKXTZF+0FRaBbERrh5MM5G6vEZWprZIwiMhyZKRqEU4bTZXg6+gb6TNu2BYkbJDWWMNNv0MhOa2OuPTHdn7ru2HnNzndhKNtiX41crYE5kYQu/QfTa1WEldSbvKDbYdLik7p1OK2ME3axWgkHLJ8YhBRwKFct1S6mh/VKIagqKQCDEEIOi2FIxR3S03HVuhDiexYlrq5SkAa2DkTO9o62yknTkNCLsvjZdsgQ325krH8Yaf3ma5ywo4WkUOMdKm0z2UxJINIUcIREDv8d3q04qHeADEaNbAA6WhIxMU9OHqlBaR3heBYVVF89jl16XULKS84hRROEhwjP1i7BU3em8t5VL5sBJKZGaaU0tlUVRKRYBiGYTIDMKEWxNiAs4+/62Qyt88rEF9Lorfqbrnt7POj6YQgAqlzpVr/gcxdZrOecdyaDL2wST8Qd7yC6/nYKIVMeEbF56zos4Bg135waMoDURvXxxEm1K678XteNdVFrNKt1YJxJYWmU22FszGqpjaGnyrhjg9v9YmWNTMzGcysjqiwY8Lh0ZYCYFgARFM2ODEbukalV3Y08hddiiml9DgHXDGzTWBnlyVFra6JwZZav+iSgmqa6Tpf0N08O1bWcHboa44bKo4bpdY518OtD3FdichoeTu/uNoTP1C4byklUW+xEx9lMMmmSlmQ8uB1S1WZzs6/7uwgVZHbz5euWgSDpbSY2SDBzLaAIKHOl5AEAileqWh09MhlIZjZJKXnIBCkJOeQNaUD+jeRewPmaKOStr2JPkTrRWx/tJ0/pJRqi2b9neGwygdAIMmWKlafaK0rgpMrryjXOw1L1knkSNAHQ23+yGBIm9RphAATxUmfEO1BzDHxTqVjSoKONBPUQz0pw4YwAFNFTJveksTMQm2AzTYAAZJSsiHi2EKS2Vn2mkz/hWQpnQVwTRRytN/FPe/4yjYz22AABgw0zbAlEUkV32V15qJM4J5ahKpGliY644xwCimk8L+KY7FfhNugxA3qbghEz1hiImdvB/UNXRCaAi8tOO7RW4cyXk0luJbnqyE1hs6oh/o1zxP0yHHU5UCtBydnvwggOXO0kypZqv/RzFCtT4FqDUgfhepsGaHWpMcczEvQbjElW4LdglvWXKXuk7pclMdT5S4zqcqR65U4Ypr+jA7+Ta6h2oU6kbDE44jhii0lbQgO1Wq3NsLVd9w/Wng9+i5inCE92CchNenFw9FVtQBVVxczl2veqzAMQ23+IVzNTKRNnZqm+O8409D7wLSWwqP2bOIaFSff1qtxc13VVXJTJDvuxSRzmGTq1/QWSLKzENMpJ/H0bP00ABl7smDUvLi4Hfe62kiGY8sXQljuvq/mThSTkOZ2Ot365m41MyP2SXfDHZKdJrj4SsmobUFG7lspJziFFE4SfIv5wc0hurZM/QbgXmfW8obHiRHfo8eJGRU4NK++ljKPBXO+CaF+W7V8Azgu9PMb5t5x0cYUomiNbTkK+3NMemJKo1JIIQWNY3eeXOtBsfjWifl2CTiecQTMOf7mXSkcd4hTkuYC9imkcGRIWfUUUkhB43iKB5+wMblvBXEfFk8IJFLbmnzK5t499jgeaDhJkGJ1Cm6k9CGFFFL4hnDi+cEpg5gCDseB/lrq+hZV7hjUnupTxx4pnqeQQgopHA84nvzgWJyo40Tro9qt2VriyEr+lvBNiCzJLh9uD/goeHLcs/PY4SgF1xwnT9QunMK3h1SvTCGFFI4xjsf8YIXmNhZI4bjFsRHZN5o1fnKq2TfEz/8lJ/gkUYxvt5knCZNTSCGF4wrHbzw4hRRSSAFIBQlTSCGFFFL4pnD8xoNTSCGFFFJIIYUUUkjhm8OxOEcjhRRSSOHIkYoEp5BCCimk8M0gFQ9OIYUUUkghhRRSSOFkRMoPTiGFFFJIIYUUUkjhZETKD04hhRRSSCGFFFJI4WREyg9OIYUUUkghhRRSSOFkRMoPTiGFFFJIIYUUUkjhZETKD04KCchvm4b/ETDzMd2jml0H9XHsoX0ppHBCobmOo3oUxz7mPHx4tutYd88UUjhqHM8aG0fb8UxqChopP/h4x4nekb7R499SOAIcVxp1WMQcV5QfAzTXceKuqz+PrJelumcKSXEc9jVN0vGssUn7ZgrHOU5eP5iZAclsJ0ZEbAgr2Z6lnPDvG4KUTRGdE6UjxbEx6Z/HoW094eDmYVVV1dG8fgxwSKEnVe+TRFVUM6WU27dv3717d9IHmrtCyXZV1nyrq6tL+norqcJJI4LjBy0Yz2NTY9xdIqqpqTl+dICZT5ShUIGZj8A+p/Ct4OQ9R4OIlCubtHfRt7p3vxAn3vzkxDJSJyKklKWlpVu3bl23bt2HH364atWq733ve//6178OqxCi+KPU1Z/HRnwtDGZxtzSd6sfRhDyPN+jmlJSU/PGPf3z22WdramoAFBUVPfroo8XFxWieUW5uxKGsrGz9+vWrVr375ptvMtO6deuOOEJ8BG+lcJSIU/7EB6SUX+O4kNjd1A9Hi1a9+eabzLx27dqvq8aTBJFI5LPPPnvjjTeef/75jh07Ll269NumKIVD4+T1g9EU+bCJyB0aP5ZOaGLtLT58wsyJv4nPQ61qPjXz+wTH8uXLL7/88riLAwcOdP+p3NuWG53Iw29Uo5LWlUhDC07wN03hMQYDTADzvHnzJkyYEAgE4LR3586dd911V3FxseJGc35PlIfqNwDguhuuf+ONN0KBIKKcxB13jDs8qhLmQv9LPD8RkSj9r1ci7tKUvl133XVvvvlmY2Ojvj5u3DgcN4PO8UADACf5Xqjvye6x+1e/+tWSJUs3b94cDDYyMyCGDv3BccK9FFrGyesHKwVl5pdffnn58uWA+MEPfqCCMccSh9VJTvIedTI3f8CAAa+99tq6detmz55dUlKiLhYUFHy7VB0uVKdLHIP1Fe2Q/Q/LmkD33Xfv9OnTx48ff+mll27YsGHOnDllZWVENHr0aDg8aX3wb/To0b169fpg9fsrV65UV/Ly8g6bqv9dhp+IEELETU6+OQGpkn/0ox/17t179erVbi06llrRepfx+HQue/To0bPnjgMHDpSU7FJXcnJyvl2SUmgl4j+SniRggMFqbjfk+0OXLVsG4KqrrnrttddaXwK+gYBjY2Pw5UUvVVRUDL92WGFh4dddfAonPF5//fWrr75a/V68ePHhztzKysqeeuqpvXv3FhcXX3nllfgmx9fW9xE9sDFzIBBYtGhRRUXFNddcU1hYeBwOeEcDBuY+NfenY38yfvz4WbNmqYv19fU/uf2n69evX792nd/vl1ImnQkk5acOHAaDwezs7LAVBvCnP/7pnnvuOSzCGhsbX3rppfLy8muvvbawsDDO1TiWyTMnJxI1X10/xglLoVAoKysrHA4DmDFjxuTJk4+Z0JN6t42NjYsWLaqo2H/NNdcUFXU9NpQcDerq6rKzswEQ0YwZMyZNmpTqNcc/Tt54MAACVVRWvPPOO+rPpUuXNjY2pqenf1vTTWb07t17z+4SBr+9bMnixYtxvM59E5EaKb9pKE1QRlbB7/cfbiH9+/cvKysD8PLLL+/duxfftoIlRoh79+6tFo0tWbJkyZIlOMFVK469JXtK7rzzTgDjx4/Xn6QyMzOfXfBsxIp4TQ8Oc3mAftjn8xUUFOzesxuAYRiHS17v3r1LSkqYedmyZYsXLz5BuX1CQ2u+EoFbc772Tqp1D7GfX6JatHs3ACHEsVSDJBM/5l69epWUlABi2bJlixe/fsyIOWJkZWW1b9++oqIi0bKlcNzixFuPdZRgSAkbzMQAxPPPv8DO0utQKPTKK6+gdSNu0vXaR4lwOKIMEIGiS00ZFBuv16u5W1/s4Yb8m3v+kAv/v6E+7643GAzu2rXrayy5hSX5rUEwGEy60r81VR/uK4q9mZmZcVcOC2pJFhHl5uYe7rsamvhDqERiH2FAsptpcU0IhUL6Vm1trX7miFXrG+0pAJz9epvdtTeO8j9O/0MoEMzJyenVq1dM9jOzcoJboEcFzAOBgLsLKA1mgIi0+8sEiVYZCuUMKYmohw8cOJD0sa+9d5/onyK/XvrD4bDW/Kqqqm+OOapkVV1SsR4/q7RDoZDOATtw4MDR7H9yLJGWlganZ7W8Kcfh4kTvMsctjheNP5agqFUnAP/4xz88Ho/H61W33Kvvj/3OQT6f56qrrmIwgBtuuKGFJxM/WaJ5Fzku+fKQSDrgfYtRQ13v3r17L7300j/96U9fS7HNcaP1zVT0PPjgg0dQewu1tKx4rTGsSa+rYkeMGKGqHjNmzCEpaQ5Hk7O4b9++RKZp1fX7/SpbA8C1117b+g54xNJ0v9j69MTW33X/KSUWLFigolytr8L9e9++fYMHD26uC/h8vrgrrZrSE/n9/quuukr9ef3117fS9B2lefy27MnXZdW/Xvp9Pp8WwXXXXdea8lsvI7URp47+7t2795JLLlF9MHFO6/f748aXr9eZaz3catnygNgCjoCNR6/YyglWCU6tp6Q1JR9lCSkkxUmXF6H1iJm3bNmy9pNPr7vuumAw+PrrrwNYsmRJeXl5+/ZtiYiYcHSfpY7ge+7Chc+9s3RZfn7+BYMGRl8X5F6jmliaTqwk155KLTwZ16JWNu3YZ4nFEfbUU0/dfffddXV1idsmHDFadkaTfqfTFxPpOcqpgn695UKSxnGjSaJOAmnSEpR6zJ0792c/+1lRUVH79u0PSckRt6g5zX9qXoIQGU1PEjHwwgsvLFu2LD8//4ILLiAGmCESaEhIlY3R7fib8e1KpPCwNqWKbVf8W3Gtdte4bt1aFeRu27Zta94iIkWlujhv3ryJEycePHjQ3QXcL6pAVPR6Yiy+GWmq6wsXLnznnXfy8vIGDhzo9o2o+Z3aKHZbj68Fh6VyR1Z7K+f58cH4Fms5mr6vW7Fw4UKl+fH7wCTkarfGULjJdgtx3rx5d999d21trVuL3FW4tSixlqSUfC1I2geff/75d955R1mDuMda1ucWKop7IFGNj2ys12+p1EocTnpS66lt4eIRlHM0r7dmmPgWw2eHhZPOD9YgomeeeQbAjTfeqP1gZn7xxRfvuON29YT74aSFtKAcLbzVAnw+35VXXokEM6c6v3tMjRuf4ixjUgV1v5UYUY4rvLkGttD2xNJab5LijLv+U1u9WbNmqdMBEm3lke2s2frgX1KBPvzww3V1deRa19+CJhyy5FbyPA4ZGRl8mHsLADjvvPNafqBlLWoNyPkmGNcoxTTECZGZXZW6g0BuxKgKx3TP1lOV+Fv9aHlkjSMg6Z+tqX3Dhg3qd+J+Di11NDAYRDRr1qzGxsZDiJsBSp443pxHq664I/Fu5rR+at0aHFLnD8uzaU3tibUcMprQ+sLdrx9ZF0aC5sdNQhIpOYJhxf3KrFmzamtrtRYlrULR4Pf7W5hptL6NrRxEkmp1WlqasgaaLXH2NrHwQ0rhsNT4kBobN54ys+59rVy/0bKne2San3RMb6XIWu7piYQdLvHHIU7CvAgBCLVXxL+efiY7M+fqK3943TXXK5Vl5meeeQYwAKHSG1v+RNKaHqUcgsMhkDj6kn5LCGG6C4lTslaW35xVjXOLm/vtRgs9yl1aUgrdTUval+JGo7gu7V4oph52G1B2woGqjkS+cGz2SHOsa46wxLZkZWUlLeqQJqA5e02uoy4SyXPvxePxeMidhUZNUdDEF1u4kpQJiTJqubS4i0mVSj/QJEQCOwFsRX4Me8GqRfHOgaBD9s1DIlGaLYuMYo8gOQITrxYmHtbARlFLFKVTSsnM7i7g1nNNXro/LWlxrXdcEl9JatkSn2+htEQ9b46wQ3al1qi3frcFTT7cPtsC3D3xKMtRaQytafVhwX1MaZwWJYVKkGil59TC3RZKaL1tOYKBtYXB65DVuQuJ0x+3iN3ecOLzeiLRMpozCM31r6TDUBw3WjYyRyyvFu4epW38dnES+sFRrFixcs+ePcOHD/f5vH6/b+jQoer6mjVrdu8u0TJNNNCHa4xaaUoSq3NrrZ7PJfYBNK/fh7SkcX9u3759zZo1h0tnYjmHROLzrT9KuuXt5RnMCccEtjCFOCwknQfrj1+HOy1BrElNrKgFOjUTmrOVSctMnAK5f7RcwhEb9KQPuIXYHFVxFLqZfwTeRiL9usDWy+sIXBP3AxUVFeqK0pnW1Ktb6q66Oe4lamArm9Yad7CFaUMr5w/6sUPamdYrUiKFbnz11Vdr1qxpjYIdAZJq1OFyPin014nE8uNqP2R17gfiim3lV/ukIw5aZONhOZqJdDYn3KRj8ZGZqVY6ba3Xn2bnsUc3i+DYJGO3+6FbwbEBaT2StiCgo5dXIo5g+Dt+cNL5wTp2Mn/+s0zyhpuuZ2Imvvnmm+GMN/Pn/0vrSZL+z2Bbrnh35Zixt40cOVJd3Lp165133tm9e3chRNeuXadMmaIW5h+WSySl3LRp0+y5c66/8QZFTwwkIGnFihW33XbbT3/6U/XKzp07J0+e3K1bNyIqLCy87777gsGg7i3BYPCxxx4bNGhQWlpaRkbGZZdd9vbbb8e0hQhAKBR64YUXrr766p49ey5YsEDd2rx58/jx47t16yaEKCgouO2227Zv394c8aqcurq62bNnX3PNNR07dhRCZGRk9O7d+84774w74pVcH3bD4fCbb745cuTI22+/XV1ZsGDB6aefTkRnn3325s2biSgUCr377rt79uxR0lm5cuXMmTN/97vfTZ06dcKECU899ZR7Nrxnz56pU6f27t3bMIwup3QZP2F8ZWVlc+6LlHLlypW33ablKLdu3Txx4l3duvUwDOratcvkyffU1tbHvRsMBt999929e/cqA5RIT8uCjuNbXV3dnDlzFNMMw0hPT9dMi6PWDffCl0AgMHfu3O9///vZ2dl+v/+MM86YOnWq3mzB/XpNTc369esXLFgwefLkXr16uZ9xP7lmzZrx48f37NnT6/VmZGQMGjRIZRDFDR5xVJWUlEyfPv3iiy/Oyckhotzc3Msuu+zRRx9Ve5FqpikhrlixQjHtvvunTph4Z5RpDEj+76bN82bPue6660aMGKHEuvHLL6b84p777rtv+vTpM2fOnDlz5gMPPHD//fdP+cU9alcEy7ImT5583333TZs2bebMmTMfmjn1/qn3/urXaz/5VBP50UcfqUb5fL6MjIyBAweqRuFInaHa2to5c+Zce+3wjh07GAZlZKQpwa1fv94tJl14WVnFZ59tVH9+svbTvzw0c+ZfH5r+xz/cd//UTVs2N1sNEbu4p64tX778L3/5i1a5uU/NdYkwZoeZF1988eqrr87JyRFC9O7de/r06co+uGtQTN68efOcOXOuv/76ESNGIJn38/77748cObKwsNDn8xFRfn7+gAEDRo4c+emnnx6SV1+Xndm9e3dT7+7SZfz48eXl5boKPUcKBoMvvPDCD3/4wx49esRVVFRUZBhGx44dx44dqytqzkrrjuD3+906k+iBhUKhODv27LPPKjvWr1+/zZubl68DKeXmzZvnzp2rRZColqtXr1Yi8Hq9Qoj8/Pxzzz1Xi+CQQb5ELdKG684775w7d27iK+p/xczc3FzDMJQWqU7dQl2tsSGJdSE2rK54MmfOnGuvv27kqJsRb4IEkVix4t3bbrvtJz/5ifrGu2vXrsmTJ3fr1oPIKCwsvPfee4PBJlJDochjjz0xaNCg9PR0NSC+9dZbbjL0uNmy/hQUdBo79qfNK2ryrWPcaKWNVT+WL1+uBn11ZdeuXZMmTVJdRg/6+vlwOPzoo4+qQT89Pf2yyy5Tu0+2gEPKyz1LUR5InOejiZk0aZLKfDvBwCcZJLNkDjZGsjPz8vPz1XdGZg4EAunp6cqYnn766dGHnbsaDz/88A+GFudm56iPuMOHDw8EApMmTfI6O05onHrqqWqrl1Zi0KBB0UIIEDR8+HB13bZtKeW0adOKh/wgOzP6TXz06NHBYDBpveedd14gEGDmhQsXdu2aZOPxJ598Ujdt6dKlgwcPdi+MuPfee0tLS7WWu+H1ehctWtQc/X/+85/Vh7azzjqruLj4jDPO0C8S0bhx48LhsPv5GTNmfO9739Mr3MeOHcvM48ePd9f4yCOP7Ny5MzGQCWfYI6IHHnhASimlLC0tHTNmDBlCbwgCQRDUpUuXnTt3uqUppXz44YeLi4t1joEjx7s8HiOaOROFOPXU07UcpZQ7duxoORn3gQceaL3Q//KXvygaFNP69u3rLmrcuHGhUChOCaurq/UDq1at+u1vf6tboXkCoEOHDhs3btRvPf744+np6YgdXKurq+PoKS8vV0miRUVF119/ff/+/fXDI0aMsG07jo3qR1lZ2ejRo1XJubm57tgtEV199dWJTGuSKQGCFNMGXTDQ7/UZgOEIRdXyxpuLuxYVulunSjiz31mWtJm5oaFhwIABTbslEEDoVli07O2lijyVYlhUVHTddde5G3XTTTepRrUGSs204Bxt/84PfjC0b9/eiYLTr9g2Dxp0kdKrplZTtO0Q9NXOHc1WyvzVzh0th89/P/0BRVb//v2JQIS5c+d+8sknZ599duJbl156qaJNY+DAgW4zoi2PG5MmTQLQpUuXGTNmLFy4cNq0aXqd5fLlyw/JusOyMz6f76WXXoorQfXuxB1tO3furHq3rujSSy9175tx3333lZaWKs8yjhV+v785g1ZeXq50plu3bnE6ozqC1gRlx3QKgbJjEyZMcBP58MMPt8yf1ougc+fOWgQFBQWqLe+8845bOd3QF3fs2BHH57jepA2XVhulRf369Ut8fvDgwXFa5GadtiGJ3U3Rk5TUlnhCGHZNDE9+97sHhgwpzsrKITIA/OhHPwoFwpPvnuL1eolI9zUiOv/8gaFQhJlfeOGFrl2LlG13t2X27Nm6WKU/7rzeZvRHEBk+ny+p/vTvf5Zm4BHwRzNn2rRpQ4YMyczMVNQ2N+gT0fnnn5846LuNsLuNR0bPI488UlxcrHNphg0bFggEJk+efPSez/GAk84PVnjmn/8SoHHjxrEzVkkpR4wYQaTcIKxdu5ZdI5/ynpl50qRJHo9H69lFF13Uv3//9u3bjxo16le/+tWIESPcXUgHtJozUhpSyvnz56v+o9Q3zhTeeOONPo9fwFB3Lxg08PTevbJzc35825jf/va3Y8eO1VmqyhbfcsstAnR2v/7jxo0bO3Zsjx49dMk+n2/Hjh2Kntdee02tTdZGYcCAAWon8BEjRowfP764uNixLADg9Xo/+eQTjjVk4XB4+PDhAPLy8j744AN9fePGjXqRr26R5uQdd9yR5vPrMW38nRMm3zMlbpDbtWvX5q1b7p486ZabR/m9PkXnmWeeOWzYsOEOVq9ezbbcunlLVlaW3++fNWvW8uXLFy1aNHz4cG3sBg8g0Xn4AAAgAElEQVQeHBWBLVkysz1p0l3uDqzk2K5du1GjRv36178eMWKEeyjVcmTmrVu3Tpo0adSoUV6vVyVvnnnmmcNdWLVqVWtsfQtMc0skcUR0+8FerzcvL2/06NF//tOM6b9/oLi42G3f27dvX7Jrt2L3Bx98MGHChDg/u7a6hl00lpWV9ejRg4imTJmiiZ8/f75+/ne/+y1zjNcopVy8eHFWTnZ2bs7s2bMbGhqYuaam5uG/z+rYIXrm85AhQ2KYBgA466zvDB8+XMtx1apVqi41Bidt+xdffFFUVKSJKS0tjeNMIBBQRrx9QYd1G9ari6WlpUr5J0+erP14V6PEtGm/b0FMiQiFQkq1cnNzP/jgA80oLbgYbZdS0TB+/Phbb/1RQUEndbegoEA1XHGg5Rq3bNkyadKkm2++2cW9s9SLw4YNGzZsmOIeM+sxbODA84XA4MGDp02bNmPGDDfnATzwwO+UHBV58+fPHzBggL47bNgwLVz1vzr3zuv1uj3O0tJSNQFLFEQiWrAzN910k7YzmgaPx/Pxxx+7OZCZmRnXu3Vpl156aZRg5ldee9VtcxBr0CZMmBBXkTZobpSWlnbv3p2AeyZPsaStBKzigqrGadOm6YfHjZvgtvkTJkyYPHmymwDDMHbt2tUsayRLy44TQaJKaBHs2NE0ZWqlCJQclRbF9sGz4gyXet6lRQOBJi1yayCA6dOnc4KVU92NiKLdTTLLGBviZl1zpCponiiea7VUj914441er1eABIiAQRcM7H16r9zsnDFjxkybNs09IAIYO3bsqFGjAPTv33/cuHG33XabIlLd9fv9mquvv/56CwNiK/THZrYVA4ko0Q+O54/T0qT8ueGG67xeM+qTCFww6Pxevfrk5OTFtZGIyBC3/WTszbeMgqB+Z/e/Y/y4234ytluP7rpYv9/vdF7bbcNbpEcAQtNz9913uxt+4YUX9uvXT3k+asR094KRN93Mh57ptEoNjg1OUj/4sksHC9Dq1avdHH/ttdfUPE+pBcf6wbYj2M2bN7dr1053lalTp6qpmMKOHTvc5yGrHdFbCduOnH12v6ROgJRy51e7Ck9x/ADCpCmTGwKNmvrdu3d37txZ19ujR4+Vy1doXbRtW4UTFKZMmeIu/JNPPnF3mL/85S/uuf727dvdc8SLLroojmxlYpAsMhQIBC666KIoyUR///vftR8spdyx/auC9h3U3VMKu+bm5y1cuPDTTz+dNGmSEOL2229nZpulZGZbnt0vSoMyLjFRSVveOX4CEubfM2bM0GRv3Lgx6gfbUkpLSitOju6UEmbesWOH2+uKk6OUsn///iqkl2js3ErVXJdumWkXXnihrnrWrFnuu24/eMyYMYFAwOXf86effupWg5E3jYizR8qBUKg+UOW+e9FFFxHRoEGD4uiZNm2aoxve6uqYif6bb74phPD4vB9/+omeTzIz27Jq/4Hi4uK8vLxNmzbp57UizZ0bE5/Q79q2rR0CdxdQNnrhwoWa+IULFyZyVY1hc56aKx3OK/VTe4HFNuo3AADh96fX1NQkFpUIRd6oUaOISAixfPlyd1xQSqkFpzQqTnDMPH78nerumDFjWlOjLlnBHahL+mT//v3VRKioqOvatTHu3aZNm3QEt6CgvdsPZmbbts8+++xEy6PqVRGmAQMGcGy/GzduHIDm4oKJaNnObNu2LdHOqOpUeFW3Wl2M692sbUWyitxdO9GgxcUplM4MumAgy6Z5OzNPmzZNsdfv91dX1zqs423bthUURGd9Xbt2zctrsmNEpOxYs5DMttQiiNN8TVXXrl2JKE4EzKxE4G7dIaHju+6MMje0G1dUVLRu3Tr3LaVFSk8KCgoS343vbg77XDbEr2xpa4IFLavljh07Ck/papBQX1Z+MeWeQEOjLnbXrl3aEhJRjx49VqxY4S659QPizJkzWx4QXa1o8oOT9lNlYzV/3KFfXZ36TCellNLasWN7NIYtAIEpU37R0BDQ6qraqL98du/Z450Vy6XDdUvaakqmviap6IaUltsPbt48TlN+sJsebUP0iBnn+SgtBUAQVfvjPzYeEsfe/dU4Gf1g9YUosRvr1AgAHTt2TJRKVMOkVPucwwlFSMcE2yyllMuXL9e6knS0bg5SytGjRyutTQwJSMu+4/afmSATdMUVV7jfYmab5VP/mKe/ur69bGli+bp/Knvqhrbj9957b0K9XFNVG3UZAQLWr1/rTCvt1157Tb2Y2JdUd925c6eeKRYUFITDQXc/1O4ggI8//khfX7t2rfZObNtmW57Rpy8RhIgaFz3mMTPbslth0em9e9ksObafn3baaYpm96fJqHGWfP21TXJMjP0rOSrEyVFK2bdvXxUPbs4paQGaaRdccEHSBxTT9HjjTilReecKL774opsk9WPjxo1q4k5EhmGUlZU18cm2H3jgAR0zdjv3C56ZrwIg//73v/V0Qd1qaGhQEiSixx9/XDOqtLQ0JyePyLjjjjtcrGHtPYRCIV2FutKnTx9nhJjtOJFN8Qmbpc1y9OjRagVPYhewbfuUU05RJVx11VWJsSgAHTp0CIZDSpALnpmvcg/+/e9/s2s2K6VUjVLxpMcff7yJxBaRqO3ueR0z79ixQ23iQUQFBQWRSMTNmR+PHqM4P3r06JiJXOugw/lJVS46NxMAQSULMbOUunz76afnKQ+eiNw5M8zMthxz62jlUlwzbDg3vcU6AbegoMCdQyKlXLBgQW5ubiuJV0i0M24+VFdXt2vXTofi1q+PBvULCwt79eoV9zxLPq3n6YrVjzzyiPuulLJDhw6qHLdB06pbU1PXvn2B/nq+YcMGrRsLFiwAgQSiOuMUKKVsbGz0+XxEBpHxxBNPuNvltmPuAOHatWtbFQqRzLb80Y9+lDQI4hZBnLbMnz//cEUQ1SLCnKdmJ9U9PUz84Q9/SLz7j3/8Q7c0TotibUhMZ2psbPT5PEqwjz/+qJ63uxkQT4rNbPPo0aNVXUkGRCl1NvaQIUNUmXpqKqV86h/z9P45y5YtS2zLWWdFP8AOGNBfSksXK6VUikqEX//6/zimk9rMdk1VdYd27d2KKtm22VJT5di8CFvPOVWclYg0f7RCNjQ06C+Qysa20EY35s2bB0R3lVm6dKkmVY9xKn6kJ1HsxBTYFYeOk5eU3NAQ8PnSAAFBjz3RRI/2fAYPHpzIzxUrVjg8ES+8EJ8xEmc9El9v5d1vAifXOjlmBqAWBJSXlxNRenp6bm5uXl6eSipX8xsAZWVlS5YsUTxCQnZ/RkaG8la7dOmSWMvgwYN16GXr1q04zF0mpEySZa9e93g8ajOETp066evODIyGDRsGx/92H46qq9YfE+MS/Nm1+Zfa2dRNLQlkZ2fdf//9uqKlS5fquzNnzlQ/Eo+3UM937dpVZwGWlZXphXqqCiGE4uTQoVecc845+t3+/fvrVCQhBIgaGhpYbQah1p67NpMC0axZsx566CFC09EnCjq4mPQ0Zi1H5V1R7NqvwYMHK+8fwJYtW9xkM3NDQ0NigUmRKHrNtO9973tJXyksLBw5cqRbFZOW5t75SHFbOei33nqruiilfOedd/QzQgj3Vj7uLJQnn3xSXbnkkkvcd5k5LS1t8ODBijnvvvuudqNnzJihnPKkWZ4AvF5v3KkfjY2NboLj0mA0bNtO2l4i+slPfqJ+v/766+7QODtfridOnOj1eNWlJ554Qt3VjdLlqEYBYPC7776blIxEHFLbi4qKNPPLysriVuEgQRm41ZuwAjh48OChH2aAkJaW5mQHNklZ2Qc1xii71AJt+i3TNHVzJk6c6H5gxIgR0ePfWwd27dWld1DWFRFRdnb2/fffr8nQduaRRx556KGH4p4HoCwGu5JfdRdWFVHs+eFa5bKzM6dOnao13G2U4nTGXZ3f77/sssuUMV+5cqW7aXrvhSFDhrgzHPr379+qA8ydZWFJb2oRlJeXu0UAYOTIka0XgeKMuw+2jMTUTziDCAAiUlZRy+vJJ59UZjmuuwFIS0u77LLL1INxrGuOzkOCiHSOop5fuXcXVgqvkPS08GuvvVb92LZtm1YGpSRZWVnqivuboaYtOzt76tSpuhz3gNjCXi6zZ89WJcTZWADp6enKpBDRe++99/+z9+bxURVZ//Cp25B0AmRTQgJmQSRs8wohLqCIoKigOBIIKouDEnCJ8MwIzPs8A+gzLjgOIcwMmwsEZx4WlWCQYVOWAAKJMwNB2RMEkqAkAUk6CUl3J+k6vz/OvZXqe7s7DQkQoL8fP9i5t24tp06dOnWq6hy5LGMbZRjb2FBVQEAcOXIk/SkmfVFDkvngQjxCYKD5kUceoUrI4rFNmzb0Q9jaRWdRu3QzptyViqKI2VOmp5wGL98RULPg1oqjQVJSXIRkjFmtVpvVSu8QkdZVHBAAV61aNXToULmbwdA9rRQTIDCmqmUMGOmpUVFR5COJtCXm3jcTAiAgQ2CMIQPUSkBGM1pDzcHEkIFDrQuj6jTMIgAhQcH+/v42mw0YVFdXGadYOiiJzjZFylzwN1kXVKEA5MZSAQZJSUlTp04FUICxvLyTdIraYrHs3LmTPoyLi3NH9qSkJOFFYefO3U899bRKE4BWfq2RgwKmjh3v8NJ7ier+Xas6ACDgk8Of0gSZOv0TwW+77Taiu9VaDcDVSxJaVpxz4KgAMzGFIciR2EgRi46OJncTYvLQ1IvLGKi6xDLR6EamMQ0iJiUlffrpp/T222+/NYaWYJrbDbmj6cm4ceOWLVtGr44cOSLypEaL33RzDhBsVvvuPd8yReGcb968mfxXg8IsFkt5+UWLxXL4h0MUumLt2rW//e1v77//fm0ccQAm38ci4qouCwQz6xxCM0Q9/Tiom3+efHgxxiZMmEBLMkT8/PPP6Yg/ACCD/1u5wmw2T548maLQ1disu/d8S8Vv2bIFEZFxAKi0VJWVlVkslYcOHQFQGMCXX36Zk5NDJl5jGDaBiooK0XEeAiOPGjVKcPvu3btFFAAGDNHBDD6bLksV1iVGY2QcEigGaYOgtAsKCQgIIE4uL7/olJvC6gHJ+QYyIN/MTLv42LVb3Mn8fEBYvHjxDz/8MHfuXDp/clk1p8aqLOcmLCJjbPTo0VOnTqU/8/LyqAgR4EOH228PQ+2CvyAIfaKGmAEICwl1SbqkpJH/9V9TABAUdiI/D8EBAFarXZv12ebNm0HjxoqKivLycovFcujQIRNTEPHLL7/87rvvxFlkk8lEpmWyjFwuZQREh8o5hIWFxcXFnTx5EhEXL178/fffX04X0HgnBkGikDoGG6Y2B1PvljWOoKCgwMDAmpoaMY8wxoCDzWbfvedbBGTMtGXLN4gkQRWLxWKxlJWXlx86dIwxEwBfuzZz6tScBx54QG41SIqRxpMNShIDJy8oAHqH2QAgwr8KgoQGh/j7+dvtdgCwXKpyACiSsEWAO+9Sr81UVFTpwrVq/KOEBN2m0lENZWVCRGAwMmnU1N/+F0jrAbqeLVeJc85plnZiLdi8eTMyQEQFWEVFhcViKSsr++H7owxaA9RnZGRMmTKlf//+THI4KKQcCL0AkTEWEhLi7+9fW1vLhWlGNQ0BA4YM6ZSw6CwhmK1Wu1C4idWpFAt1WFnZoUOHGGPI2Jq1a/9ryuv9+/cnDlHApCiKiSnAERRkDDRjFIMGzYdXV1eBJN+EjGrgba5WRSfh3bPeVcStpQcDwPbt2wsKCrp06XLkyBF1v16asyvKLTExMRVVlQCQmZn58ccf0x1nMXIQUGGsXbt2nHOjQUsMJGEPFj5EPHQwA0aVYMDatWvnofKebZCMsfDw8KKiIuNz+iHWakbI6r7Lqnbo0CEqKurs2Z9lNfr48eMiQdu2bV1WCQDuv/9+8YS89ohSAgICSPMgE6A384dxjtct5UGbDmneorHn0rYRFBSEahiHhhwaCgKkw3CI2BRfMDpNVyZamzZt3LWXdE0q/aeffvKQuTEHmp7p1aVLl+ghZWU2mxkDxhTOubA0HDt2DDQnX8nJyU5rDMSgoLYR4ZEDBw6M6Rzbu3dv2kw8efLkhQsXKAEpN041cV406qwpApxzRdGbB+QLLkZER0c/8cQTZCBfuXJlSkoKtevYsWOHDh16acKL7du3pxFNzqro7cSJEwHUYz2kAgQFhURERAwcODA2Nvb/6/0rcSbSQ9FEJXDW54wgbqeshI8qxpiLmC7OxGkUYiLRfStPNp5zCA8PLygooN86zpEv2oI0pgBgzpw5z44eTfnv3bv3gQceSExMfPPNN+VTkl7CyazgCh06dIiOjiYhJtv7dZkwxsrLyg4dOkRVlUc306D7Sm4RY0wTaGdJUJAJQ4xNRExOTpa/Jc6MjIx86KE44hlx0Ba0nSVErKurA6lHvKQJpaMucFrVaJgzZ87o0aPpreiCt956q0+fPpercwvKi1Jc0gpdme7E7Ca4SODYsWMi54kTJ2q+w0i35kFBQRERHQcOHBgbG3333XfTytllzXVPdNJA94mVbFiumkApO3ToQLykaylp3jRNywQRH8qMilLMTvE8IiIiKiqqqKgIEWUHlHJNFEUR62qZPsnJycBUUxA9bNeuXWRE1MCBA2M733H33XfLI4uWHGCAKEi00WWaDh06uHxF1zYa6qOlp4eiv6I6q/2FiIypM6auPnKThY5BM6bor3Pnzr3wwgs1NTW1tbW1tbXV1dUMlZqaGlutdfXq1U8++aTI4YrXkE3BLacHL/noQ1DYxIkTG643MtVaAwBBIcGjn3s2fdkyBRSr1ZqZmSkf/GKqztoADggMOCCFtwJt1UguNpGBamIRXetsJxM/UeQHXJahMrsxBF6PFFyKkZdQZ25xyUDyw1Z+rTWVW62zk70ZgDHGASljRGRMshAwbN/h9p9++knWVKqrrUJ26HR0eRUYFhYSEOBvtdoBwFZrR2nk6LYC5VYzufJMM1O5HiEcEYU9w2KxbNm0ed26df/cuMFut6sKd5u2CIrTx4xC96kMQMqwKNEEgAjqmS0GyJwoeVkDVZdY3pd0t7BhjIWFhdFBHdDuwbhbdIk/RVlmszkkJIQOWFMTxORqNgciAmJDwCoELLNcpASB5jbr168PvS0kMDCwXbt2bdq0CQoJFnwizxDkPpmKttvt/v7+TvVxQx/1cwTgorYNH3FwuLTGyoMGAZInT9r2zVYAyMnJyc/Pj4uLY4z949O/o4OrpkQGCCD2iwMDA9evXx8aGtqmTUC7du3atg0Su/O6JYpniM5CRLG6MIIOWVEvy3ZKITq8UVgvF9QERVGjYMo6n5xA+wNQGtqonVAy1o2+GjkqafFHH019LUW8Wrdu3bp165KTkxcsWCDOYHhfT5DUC5A4GRGB4e3ht509exa0sw2axUGdfSsqKjZ/vSUzM3PD+o3Ci21gYKBxGQzEObpBL5XYvn17UiCC2rZjoDCA8vIKAAWABwYGblj/z5CQkMC2bWggCO+EmvXLibwOh0PhjEGDW7fLoglCI3N/UlLSRx99lJKSIgQmdcGkSZP+9re/eViVAVK3Us0bVk1aW6hcvU9DnWroNm/RiQqUWdRN+UBzwPr168PCQgIDA9sGOZHOAwWMM6MoXe175lQfmh3EASrUJl9G2irniqIgADOpTVOcG0n5tFJMIIk1TpJfXsMwBMXJvitvwQn+adu2LdOyl2ciq9WqsTtYylT6BAS0Wb9+fWhocGBgYHC7oLZt2waFtNPqqFt4OACA1yMDBYHT5puBSIZ41GpfY4NMcIWysl8oAfVXSFhomzZ6Vgdw0m3U7gYkLQW0fQY5W3NggNAxQOVtAIDdu3fv3LmTAoEi44jIQFGYgow7bSdeJ5PwraIHE1sUFxevX78eAGTtFuTRxdiLL76YvmwZIiLg8uXLdSl1UAcq6NeamtbboAR7XuUIjcQYJ0xqg/O0bcjMaEUA91zlcv9Xl4NL3Qu08xUAIDsEJd89umlA/BkWFnbuXAljTPZHBtK5Ebk41KzjRj3PAy5cuLB27dp169bt2LGjbWCbZ599dt26dStWrPjss8880EGura4mAogoph8jTTyvPYxQ10ja2V/jV0LahoaGkjol7syBM2GtVqu7gsQyr3Pnzp6bz7RjdoyxGmvNQw891MrP1KBSuEoPzkd4jx07Fh8fL9tL6BCFsSyjLUrXdvHKQ3eTS38yFn7yySfz5s2rra39xz/+kZCQIJvoqFGoOS0xm81ij1jXFpfVMEKQlDEmri7JH4rP5Y5rKAihurqaMcY5r6ura1QmuIQ7shiNWODKKhkaGlpQUACSvijrjsb8xecvv/zy/Qn3TJ48WQ7vkp6enpubu3v3bs+7WC4hM4D8A6GhCeTRid7+8suFL7744quvvsrKymrTru2oUaO++uqrFStWfP7557pMvCxd/vOuu+6iH2KHhHimdevW5IXW+Vt9JjorwBV0q05aGrNCxFdeeeWee+555ZVX5MAly5YtO3DgwK5duzzvougqLDJ3mUzHRcbPVS5y7kRBOpIhfn6tABpCpnsFVycPPYsCHUj3BRGQzxBY1Lm0hjGi6zKddBK/5cPHINHnrrvuEjnIH4plMMj0qakZMGCAv39raJiF3cZSNS5Kja1xSR+3rQZ1Qd6qVSvUjoYOGDDAz+wPrqgt6wkiT476kHXeGBQGDx7MUPH39zcH+gNASFBoYGBgz1/1EOeer2DgNBdulXtyRN9PP/0UHXzQwIeNASZQm/IffPDBO7t04aCGCnMRM4ZDfa16zq++vh5UO7FTQaBZgMSpJiaUV/cdTZHqVPMMMKYlpzEH0vC7LHtSQ2IuZhmosFh0FaEZmqFibAuA6gmFntAFIwDo1ClS5O8uphFpkGTX5Jz36tGTiWEDoAgjmRQHi2HDc7UJAKAwYUekjMVRV0T28suvRkZGvv7668XFxStXrjz/y4Vly9OHDh3q5+cHzASgBLcLknLkiA5AqK9V1ZE6Rz0ZEVFrOP1mCMg5ICiaUm60lBib7Hkw33HHHeIT2UGPkEFC2opdvx49eoiCxL+ISHJWyCC1XI7o4BfLy4iEwm+dy6oScdvfFk69j8B37s4yRiuQa0iIjIwUz8Wp2YZqqPuBXIh4lZ4N2SriOdeMYQqYFDA5UGJTdRxJ4wvAr1XrF16awBUAgNWrV6ODb1j/zwsXLpAPKVGN22+/Xfy5e/duzU6j6DrOONO4g7gUi5JnJXWESlMgIgqTP3WcqDpF/QBtgmwuJdgpKyKoYavXAeAAqKmpYYwBYk1Npe4rhbvInBKYAFoBS0hIyM3NXbNmjVgGA8DBgwf/+Mc/XlYrPEN4XAJJzkye/EqHDpFTp/723LmSFStWnS+58Gn634c9MczsF0Ctk7VAoyrvrqyLFy9SA9UrQWJXlwEw3Ll7l86W7C43jYwcwWF86w1Q88onPzGaFRISEvbv35+RkaHrgrfffttt1k4zTsN4BAQFRFxlrh3ndTsc5ElHsLfNZlNlpiZDEJFkCHny0tbSHBsLsYYMUE9trm30iSo7b5C6Mr3LzVVAPXiAiAxBzDuyLBL5WywWBYBJyzD6H3C3s63YcXrkkUdcjmX5w9tvD1dUfYDv3r1TlZNydaXmaTU00YlqaQHvrcrr1I9aIRUVFQoZdzlE3NaB+osD7vq24fYzY0zmE5qXxdu6ujraMa7nDlBc7BswdDILCuk9duzYrKysHTu3b/56U2Zm5pdffrn87+mLlizUCW2XdL4GuCX0YNFbdEHyhRde8Jz+ueeeE1/pwvYSxDrP5b62h0WSB6FMUq+qqsqlsksPxT6gOyWMTqe5q49nOEkc52NhAHDhwgVyxB0XFyduh8TGxor7Lt988404hKDLh3Mu9pEfeughY+bG5bjLWtErI82PHz++bNkySrlnz56xY8cKP1+UBUgXuuUihKKpM+eAZkswVslosPfGcYRsYomNjRWX5bds2SJclMiJaV6sqlJDOj/00ENg6EehN+v7l7Hjx4/X1dUBY927d7/77rvlT2QThUCPHj2EnV7cl3cJEoiI2LVrV3Eo8+OPP3Z31d15GmsQmu6IhoYDNi4NJMnJydStxcXFmzdv/vjjj0NCQmS3FYgoGoWIH330kZGpLlfsyh23detWl5SkzIVjB3KCe1mluMyQfojZyOWVf0SkidmFUUdrKTngAwC6PCRDSB4dLJbKH3/8URSRlJSUn58v3J0CwJdffnlFzXJRfwA4f/482RpJziDi0aPHyU8caLGdhZVdNJOMW3JjG12vUkGyQEPE7t27i8w//vhjI0G84RkvBb7xK9ouEH+KhlgsFmGOIa+dui7IzMz0shQZtCiSRRP9Fsey5UHdsMwGEHJeHgI9eznJEJcNd57TGiEm1cclW4pZgyZEd2WBNiG6FZUSGJ1O8WCmcoaOUeVX4vqgfCahZ88e/v7+JANJxrrjDd1z41DVgYgg87nLxQxNxOoThO49e8gC3GU1jBOTzWaTRZBu1hb94qF1ulq1ENwSejARPTMzs7Cw0MSUUYkjXaTR1kwouRpBxPRPl9U5ap06VT0ypIjgc9r6iYNh1etuQnIJREa3VnWJVcZiaLXXcEAOWM/r0NXSUFU3ERgquoFh5EvyoCTSkKTjvF5XNP1ek/ElNZYc0IpM6PYGgGKxVK5evVoUJMudH344XFtbj4hdunQZNGiQGJCg3qhFZHCpptpofRE5MLr3rTBUWFX1JVQ7QT3RSPdwGUJocIh8tgkBLtVUM4bicKH2XFEtzopJbqxs2yN7PGr2ASP1QkND6SCUh6OiulaAJiw0okFlZeVnn32mM2zQv4cOHSIhTkQTnSh3jbsrFJmZmQwBOL755pu6BOrUxYBORSMiMGAm1fbGQPlq3foNG9YL+xAgGXZR7FQAgKIofn5+jz76KFlYf/zxxzlz5jSQt6FEBUApKyuz2+0MQAFGvYOAlZcqyK4AiN8MvH8AACAASURBVGIfQOUZBzI3IlRwaa8ePe+/T718+cc//nH7ju3jx4+Xj9zQJCTsfF999dWmTZvkTNR0XL0257o8A4T7zPLycnLTZsQPP/wgOo78zemaANL1WW8gciBfLu4+FxYmZNrdcGl0m9Qj72TlAob6E3EOV2ZUBHgtJYV8MwniK4ry5ptvvvHGG5Sy0bbIwke0RdaxZFkhtGrh6Fqz5UNwcHBIcAhIzFVVXSnnbFwwMABe75DzJ0sbIq5duxYAFEWRWZcxHDz4YUBADuvWrd+4cbMLpRZBHFETJdbX19eDox4c8hIFm+PcC2PstddTfj3iGbXRikL/vvnmmyIShDfshIgACqKiruUYVF4SPoXU6FG62nqTLclMcmgp7Pfr1q3buHkTp70J4BTmAdQdPwddq9FNX5SPk+YNJoYmZOr6GaV9IVFPoYg7HA65f0UaWlQjoqySOpnIVSjqZRHpCBxRhv6SupsjOjjnxD8Amq9u4BwctGMp7PpWqxU1fgEFBj0yGAV9Nm4UZdO8rxajHbkVLaUbNQBAm89yGwnEcjrmV4/xApKYRYkbEREUAIaPPvooPfnqq6/Wb/hnQ56cjPOqfxFwv8xwWR+iJ9Gtpem77nBL6MGE999/HwDuv//+4MaO7fft25diiTHGzp8/v2LFCicZwcDDMlT3/NKlSy51UCOYZvGVlSFjAoKwHOiSiT91GhLVX77WCnSbVTJzUnpx/k+udlVVFUVvmjhxolAFKM1rr70m/pw1a5bwISAXJETGBx98IIqjH6REonbCRJ65XdAIABCzs7N1z2hdi4CWCosIiE2vhEIg3ygXSwgRn8J1WRIoYrChMgAA+/bta/RzGYwxl0TTQcROI6IxDcKGDc6uJwTKLl5ctGgRACQmJgoTqdFC0E6LyUkPX375ZZFg9OjRdKi64SvGGGNHjx4V3o4QcdKkSaCZGd566y3SJ3Rz/44dO3r16mW81p2TkyNl7LQMEJseAkYtijEmSs/NzWXAXn31VeOQEc6GASApKUlorlK74Nix4w8++CA9cWfVFpB38WbOnEm+EXVYu3Yt6WSC22VQJRstyB3oc6KeB751JxZ09puGFIgU8sP4NjIy8ujRozk5OTp7z7PPPksJunbt6n39RRHyOWNRVlVV1dy5c8FZzojz9BaLJVc7Gkv5VFRU0LfigDtIgov+lAsCTWMTAu3FF18UBRErCp5hjI0ePZrW9rp6HtMGgngiOFwoK+BKhrsDFa0z7Im2REZGHjt2LCcnRyesRM296QJtvmjoAqMg1VVbNmp4YDaxPp88ebIYzvJwk3Hs2DFVhriyGhovddXW1oJhppP7lyCMETrzkMjf5QaUPCHKyzw5DfGPbra9dOmSPCFqrjzV+2TCVYLO26agD2MsKSlJyFi5IcePH5dlrHjOGFMPNTm3Tk7j5J1dvVHHdD5SG1rB2KRJk0Tpzz333OrPVquFOtdHBDclFvVS7HjvX79FAG8NpKWlUXudgrWKIHDaE64GeHG8PjVFDOnY2FgKLaglwkcGDQYAxpgxyDBlQs4dAWDMmDFcg8uK6eLoDB48mJarDz40wPiByHb48OFycfSjvr5eMOK8efO0NtY6sJZjLWLdp/9YTnZOxqDkfLEDtQg6iH36xtO3ZrOZAimJzAsLC++//37G2Lhx41yGhJk2bZo4edm3b1/yJCBw9OhRus78+uuvqzXSYpVxxEGPqCaEQY88LOrjXISDIv2MGTMGtJGWmppaY7NyxLy8vMzMTOEilzpr27ZtDoejsLBw7rxUZiJTh6lfv35FRUXl5eWlpaUiT5Xaiot+JMYgJ47Uj1qQNbUrx4wZIwwLqampVms1oiMvL0+O8WZoSwOEOQ0AEhISyBGHSHnkyBEi2pQpU3T5iJABABASEiJig9O/5eXlZJUZMmSIiLYq89i8efMURWGghAaHaW1xIDo4OgY98rAmIgEU1jfh3tR589esWZOZmTl//l+HDXsKAJKSkuQIWzp7Z1xc3OzZsz///PPMzMzU1NT+D/YDBUaMGCHq8NyY58nqzxibO3cu1ZA6UTRQOHUXkXVdDB+OFeUWWsgxKVSSkdTCRiVInZqampGR8eWXX86fP3/YsGEAMGrUKJmGniGHY01ISPjppyIRiQoRjR0n5zlw4ED6sGvXru4q7AE0BAipqalEvZMnTxL1OOdB2tqG4jYbY3SZzYFk+3j//fd1mQtCDRgwQB7mCxcuBAYJ995DEbxFVcU6TcRyaxScIt4BGOUMIhYWFvbr148xNn78eFEBzvnRw0fELk3nmDt3bMuqr68vKipKTf2zyaTy3n333VdYWKiNbkQpKJqfn98XX3yBDXG2HAVFZ2gjWy5IRwqhrwBA3759Bc+kpaURzyQlJbn8xGWcLS8hdwGXsGDRQlDYPfeoXSDSe9MFRh54/vnnhSCdN2+e1WrliCfyG8agOGytcpEh8KHYeJkzZw490WbS2kGPPCRURsZYQsK9c+fOy8jIyMzMnD/vL08OdZIhHiA4TaaJ7hVKE+LTTz9NMflkznE4HBSniTGWlpZmpIkcG6+0tJTiWVLmKv8w8DP7E/9o3zoKis5QuePHj6dSiAac1yNHW41VXIlLTk6WqedA/vDgQfKiNCHh3tTUNKJPWtpfdDKWoKrFDJ4c/hQ6OO1vUFV4vYPaqPbm/DSH1ttc++/Tf/ydzmoDwPnz5xEdyOu5GoO2/uHBA+WT3Pf0vTctdf7atWupPkOHPkkqu6gM8bmiKHI0cpnHiDIMlLHPj5Pb7qWIu164yfXg3NzcVatWySZMGt6rVq36/vvvafTKvWi1Wvfs2TNvfmrHOyLlCf7OO+9MS0vLysrKysqaNzfVxNTNKdpWy87OJvur3W7Pzs5OS0sTa8qgoKB33nknOzub4j0aGYJKt1qt2dnZFL9evWPE4A+zZm7btq2mpsZms+3cuXP27NmiPv7+/qmpqZQtIlZUVGzbto0sZITY2NilS5d+//33HOscWFtYdCp9+cfde3YTevCwp4Z+uW5tYWEhanqwvA6Oi4ubNGnStGnTEhMT/f39AwMDSY64nDYcDsdzz42hnTUACAkJmTFjRmZm5rZt22bPnk2BeWbOnCmSc17vcDhyc3Pfe38OM2k2aQYvvzo5IyMjNzdXDubOeT3pGeTrW0gQUJh/gJkxRqHhhwwZYrQs9ukbP3DQw4yZFKWVkM5Lly5FdGRn7yWNkKitKMp7772XnZ1NuoXdbt+Xkz1vfpoQMUFBQXPefudf+7JF3M7NmzcLPVhbQviBIVS91gpD0FCHg46hEwTRtm7d+tZbb9FUJBGtAbLSL2aIWbNmpaampqSkhIeHm83mt99+W54PZA6nBaECptjozlqtHIgOB9aXnC9Ww0OQZFRaSbfKFMZMr7zyit1ul6VbcXGxuLIjz3+iT8e9MJZISnXY/PUW3QV8coUxY8YMMQRkPpw1a9a2bduqq6vlgcO5Grj5scceo2QLFiwwEopw/vx5dzEvqJ4vv/yy4DeX7O2u4+jzkJCgGTOmydxu7DiSKvL+OwCMHTt2xYoV3333ncztnrFlyxZdE8hhy4wZM2j4C8IGBgYuXLhwX062rdaOiJWVldu3bxfrVQCIiIhYunRpbm6u1WrNycl59913ZbL/4Q9/ILIjBU5TGK0Vs7OzqRdycnJox0xekHsDnY8kd3LGmeL8sUfVsIsMFAXUHfz4+N4PPfSgjiBLly6lj4QeTImpoOnTp48Y+Uxr/1YBAQHz5893V8nS0tJu3bqJb42cQwMBEUmOzZkzRyzJ6K1RjnkAdYGO82fOnLl161bqAnFjj7qAvsrOzvamC4x6sDggBNqY9Q8wg8IEF4m3xEXZ2dnUkMrKym3btsnrQMFF2jRaV3rhnBhujDGyQWgSQVHA9Morr9TW1jaqFRlpQpMIsaXdbs/Kypo1axY0yF7zvLmp2Xv3kXymCTE5OVk2ZtGESDQpKCpMT0+XJcPw4cO/XJdZUFRIFVAZlc6MMBYXFzd58uRp06YR/whGleSSo7z84pZNm58b/azgGbPZPGfOnP3795MZlSOWnC+lQjVpaVKUVlobG2QsItrtdqdJn4F/gDktdd532TkkpiotFdu3bhPbYgAQ0zn2k2XObfx0ebce3UFhNM8++eST6zLXFhacRlUvry85XxzXvRsoTO4jIfMVpZWQ+dnZ2ampqYLPheZD6yhbrar5iHyC24UIzUdMmi0WN7MerDsGIPhJvfzfcPm9IQL4ihUrBFeJ5DpBKCtbYpjt2bMHHXzvt3ucvnXeX/CgB9NpIdJQdZI3IyNjw4YNxJfATLpdm6NHjyLinDlzQJGO+GujV4Sep6BWRpn+3HPPUQKjS3xqWkRExOzZbxUXlzZK7aVL041eOABg+PDhBw78hyisNplzinFgpKcCTAEmh62XDXXySgAUUFqxlJQUq9XK0VH0U2HnLneKE8ZdunRJT0/niNt2bCdZQ88TExMrKio0k6qqE8hk2bNnD0f8du8ed1Og6EdEnDX7LflkkaIA1UfUXN/XBlmwdOlS90Q7gK40s+PHjwcHB4uwljKCgoKmTJly5swZY4nix/x5f1HAxEDpHtfDKaY8IkcsLy+fNGlSw7l3jTiDBw/Oyspy2QSLxfLixAk0pmRi9urVKzNzreh3URZ1ImOMuFRpZXrt9RSrtXrDhvU0NkHRD7k1azO4jhk4IsfFixdTAtmabgQ1SlEUylzUUDTqCvDJJ59ERUXp1l3Ucfv379dRSRu/rvH99997P0k4DQEARVFSUlIEHUAnrxjs2JmFiAsXLlbAiZtJRwFQMr5aCyYQNqEGnmewZu0XiI7KSktERIR45efnRzfJQkJCUlNTvVk5yEwoyxmZYSIiImbPnl1cXOzgst5G+zaOoqKCO+/qLBp2511dli1PR47bt+6Qu2DEyMTyCgt9qdODpYLCZ8+eSU4ePcBisUxOftnEWnnmGRFuXZQiy+fDhw83ThzEDZs2uhQ1AJCRkYGIlZaKiPAOohR/f3/qguDg4Hnz5nnuAqMejIizZ89mEkBhr73uxEU6ZGVlcc4XLlzo4h0DUJit1s4RaSYtL784eXKyLENk0umkojvO37BhA5O8CeloIgaU4HZBQJoQ33vvPV2/UD1vCw0jcjz11FPOEl6dC8SE2KdPH5U0mi1Ax6jG+k+fPl3np1jkT0EriD6Wil8mTX4RWoFOZj7yyCMya8kHiJ2aCkDzI53zlEuhNMZJX+YuBdjzzz6nFeJAdJRZLiZPnuiyv3bu3EltlDchyf8SMyl06vvbvbs553v37jX2lGh7y7cHN79H9xaFgoKC1q1bBwYG0vlXxhgHR1VVFa/H1q1ba8qEGvkGEe12O7l0DQgIoN0furLjQE6Hfux2e3BwMBmxKEO73V5TUxMVFWViSl1d3bmSYkQMDQ3lXPU1bbPZbDYbrd2NIOrbrNbS0lJFofA5qnJWXl5OmijnvLS01N8caDabmRrrmNGt3piYGMZYZWVlmaXcz88v0BwAAMigtrbWWl0TGBhIIXPKysoqKyv9/f3J6T0i2mw2q9UaGhpKV5d69uxJh01TUlJ69OhRXl4eHh5+3333kThAF44dAbSzTfSWjnt99913ubm5JSUlAQEBsbGxjzzySEREBLm/YazBYwPnnDyQh4WFce2spM1ms1ttABAdG+NuYjh58uQ333xz/vz5qJg7HnvssZioWMYYAkdAu71u1apV50tKExISHnvsMfIliYA5+77LysoKCPAfMmSI8C9bUFBgNgcKn7IAUFtbS/3ITEpdXV3JuWIACAkJ8dCPCJCf/+O2rV+XlpbGxEQ9/vjjUVFONUc3F2V0z3Nycg4cOFBaWioTjcanPK1q1Fa/LSsry87OLiwsLCkpCQ0N7d2790MPPeQybJ6Mv6T9dfqM6QzY0KFDN27ZiKBGZiFaKcAQ8fwvpTt27DhzqoAxFh0d/eijj0ZGRupKd2oO8JLSkp07dp0+fRoRIyMj77333t69e8vxWtXLKwCIePLkya1btxaXlsTExDzxxBPRUdGADrvdXlxaoiiKurmPCgBUVFQgYmSnjn6t/USpnHOFKYBYeanq0KFDJpNJbI96wPnz57dnbSsoKADOoqOjhwwZ0nS/lTk5OQcO/Ke0tDQgoI3E7XrYbLaSkhISC/SEMUYDEACio6Pd+boX4FL0yvz8/K1bt5aWlsbExDz++OPR0dGUf1BQkMlkQgdnjHFAm81mq7V37NjRv7VfTY3tQmlJ69atA9uq4dPr6uqqq62KotwWEVZSUtJaad2uXTtxPbGiooKDg76lkXvw4MFDhw4VFRXV1taGhob26tXr4YcfbghF5DW6d+9OuqOQM+3bt7///vuF2kre25jqSq/hFLXVbluxYsUv5y8mJCQMefwxhQKhI+zN3rtr1y6z2fzoo4/2iY8X/NyjR48TJ04wxlJSUrp3715eXt6hQ4d77703Pr43AHh1MQbhwvlftmZ9XVBQwFDR8Yzol8LCQmFxoIckKADAncx3LgSsNuv5klLB+QSLxaIoSkREhJ+fHw2eg9+77QIPDCzGtdNDbQwSFz32xOPRUdF2m624uDg4OJhsfiR/7Ha7zWbr2LGjn59fdXX1hQsXWrdu3bZtW2psbW1ttbVGUZSoqCgAYFp0HkQsvXB+x44dhWeKACAqKuqxxx4j0nk51ux2e8m5YkVR2gUHiTzpqGvHjh055yUlJWazOdAcgEgOStXpMiYmBhErKyvLy8vNZnNAQABV1V5XW11d3cY/oENEBDB1QpQSKDabzWarCQ0NpQ3MHj165OXlAbKU11K69YizWCw0IRoNRqJF5eXlFRVVrVu3bteuDQ1YmmftdjvVCtS5GwGg5OLFbdu2nT3tRB+Uzj3ToDabzWoXMwCAinILyWR10i8rE4qKaGOgOYBIXV5eXlFRIXJgjFmtVrvVFhoaGhRCgYS44JALFy5s27at8EwRY0yuj+isgoICPz+/wMBA9Xq6gna7vabmUnR0tAKmurq6n3/+WVEUEaKImiC1/TKOy1973OR68GVB7q2W33PeA6WLou5adOedd9Kab/Hixa+++qpuYnZHDW+Emlw6SFd3r1j5MGYr/9Bla/xT/L7cmrvTAnU5t2SeQUS60IaIKSkpdJ3uciusa6kH+rtLaczwqhLNHavcKPnfTOjcuTN5m1q8ePFrr71mJFfTSUc5uCzoinPjzoHEmlK9Rsu62qXcuHCnqciUcr1QdzUpeBZHnTt3pkODV8Y/TdEfro0eQnq5G9GtXyXKCW5K5rxV4sl5A3ejQveqZarI7tYz8raIhzrLxgxhO2+WOUmuifiz6asvUT3jj0a/Au+6T5fmBl0xyjRhjJ07d44a4o0N1QPcEVB+7nLuAedecEnVRvvR3fLG5Ti9Bkqq4MDmzfbmgyxnXMqlpqyUjEObMSYKaorcvjbqqTeCuqVBN67h6lOpoWjhj1zThI2LcCYdPBAwPtQJSTkHmVG9hzdGFncprw31GGPCDMokVzCNDpYbiDm9h08Pdg0PIqll8kETayX0SFlVbVQuN1qoSwHU9Nq6y82bzrpioeZZvrdMrgCDuBcO+elmYRMzBK/7VDf3eP7qivmq6fW8AtwoIqIlwEgZb9QmL5Vgl5lfVibuinapTvkArihz9ajE9H96GneXVY1mMfpcFrwv62rWShF56wppdOVwk+EW8h/sgwcIDfhqW7tv3LF049ZcABG///57AHjqqadcnmT1wYdrA3EroLm2WdwNT36lrpp9uKXggX+u407gTTDp3BDw6cE+AEjBC4STDd8IvDkgC/GsrCzq31mzZt2gxzx8uKEhAr8LOdPoHcEmFnRZoft88IEgM6pvKrzp4dODfQBwDkXjw80EWYgvWbIEAMaMGSPig/jgwzWAWHRdvShTunWdKKgZ7yT4cBPDHf+4S+DDzQSfHuwD1NTUiOiLLkP1+kTATYBdu3atW7cuIiLib3/7Gzj3qa9/fbiqYJoPKc9yxh285E+RjAQarfSEq3Lfwu8mAUr/NR9k9iD+oZ0KH//cCjD98Y9/vN51aBG4Kb2ByHB3Xz4/P/9//ud/fvjhB3p4/Pjxffv2/fDDD6dOnerSpQt5Sr65KXMTQ/T1qVOnnnrqKYfDsWHDhh49esC1utrigw8EkjOHDh2iP4Wc+fHHH++66y6SM+BGDnvDn0JMnThxggSaKGjv3r1CoFEEPh9uElwFuXXixIn//u//PnToEC2rjh8/vmfPnmvJPy3TIdXNDZ//4FsaXbp0ofBjNIvQnRLy6M45//rrr8m71k2/SLiJgYj//ve/n3nmmYqKioyMDAoy5IMP1xJ33XUXxVhhkqe8qyFn7rrrLuERRYAKQsSvv/66X79+Tcnfh5sbPv65NeHTg29prFy50mw2R0ZGtm/fPiQkJCQkpNGAZD7ccDh48ODcuXN///vf9+3b93rXxYdbESRnOnbsePvtt4eGhgYHB18lOSMLtKtakA83JXz8c2vCpwf74IMPPvhw7eDbX/LBBx9aDnx6sA8++OCDDz744IMPtyJ8/iJ88MEHH3zwwQcffLgV4dODb134tgJ88MGHWwc+iedDU+Djn5sVvnMRPvjggw8++OCDDz7civDZg33wwQcffPDBBx98uBXh04N9aCI4AL/edWgymjs6kQ/Nj5u3j67ZppwviKAAIt68DHXzo0Vx7xVUpkXVv1mgxvi7Mdvl04N98AEAEW7YMezDjY5r5kTMF0RQ4BZv/nVHE4Wth+679mL8cnmpZfoNJLpt2bKlvLzc+/SIKBNc166MjAybzQYtfgV+A+nBN4Xd0YcWBgTggKAgME4R9dTh3cJsRbLQQXQgOlqgNLlieENtZOAArlHg5ml7s0NHHB+5dBAEYVJcXkQHAPcR6mpDUFhWmCiOqTGlB9b10FMi56vUm2gwmiCiyyYYq0Etuo5KsDrfuZES//8bv//1k09nrv2SBK14Bc5ihBpLsSFFPHPn0YSIePDgwWeffbZXr155eXlyDi1wDXAD6cE+tGi0fAXFXd0YMJAkVAscpaDJDuaM612pawdZgN5SDb8C6PQAYhVibw+z9a0Dl2NHTOrS9H/Na3YLQKaw/FD3xIP50KhLuRTsut5sRggpJOv0inLDqFI03+nWIXV1daNHj57/1/m/++3vkpOTndJLUw89EY11t+SgxPHx8StWrDh9+nS/fv3+9a9/QQsW3TdM5wEoN1Rtb1TQetEoPlB75fJtWVnZ1V6FNx0uByGpwAAKYyZKU1ZWBs6r25aHm20seENtxlHhAACnTp06e/bs1a/UDYqGfTPdKQh5tiYmv5WBEgCAJAD9qTFYC5VjNw3ELhBjTN4SQmxYvIEXypPLqUcsbJq/3oZyjZBtxsZhCC1jlhQkdzgcSSNHrcvMnJg8MXV+KigMFOY6qXhgUHz1mSMCwPjx49PS0ioqKoYMGXLo0KFmb0Jzwec3zQcnIAACMtRzNjb8wKKCwoqKigsXLuTl5R07diQnJycv72RVVVULN9cJsSj/AADGWF1d3eHDhzdt2vTFF19ERERs3779OtfVB2cg4k9FZ+fOnbty9SqLxQIAsbGxH3744dChQ6931VoaOADQxgb9jYiFheqAzc/PP3r0aE5OTl5eXlVV1Q1kxGp2GDfoi4qK/vznP69evbqiogIAoqOjP/zww2HDhl23Kt68UNVEBgwEmwIAANOfnb2+pwhuPmhkVn8joAIsOTn50+XLn3j8ic2bNzPGkBkmcS96RxjgjQuYsWPHfvbZZ1FRUYcOHQoJCXF5NuY6A33wwRmcczom6+odnv7xlJGLhgwZosvhWlT0iiC3jnM+c+bMPn36BPibxYicNm3a9a2hF1DPB986WLZsWUCbQJ3RuGvXrte7XleCqzo6OK/nvF7m8NOnCxgz6aYc3YC9VaGOI8758uXLAwMDdWItLi7uetfwZgNH5Ijo4OjgXJPFHJFz7nA4Gh0anHPxrZSdlrXDaea62tOQPMpcltUC50GO6NBI5kDuQL5gwSIAJSIioqysjHPkvIGoRmKK/kKt1fR8z549L/5mQvQdUdXV1Si9pd8VFRWdOnUCgBEjRqAh55aAVs2mUPtwswA9LsFjY2OPHDly4sSJtWvXfv755/Swffv2Xn5+3aE7NNa5c+cuXbqUXywrOltET0JCQq5T1RoBGuzZNz2opbNmzfrTn/70asprgwcPPvzDoaVLl5aWliLiiy++eL0r2OJgPBMZGxtz+PDh48ePZmZmfvbZZ/RQHrC3GtCwW00MlpKSMmjQoB9++CE9Pb2kpAQAJkyYcN1qefNiWfqyLRs2BQcHIwOz2dzKrzUA8HpHTU2N3W5/4IEHpk6dCgZxh+434rdv375kyZK2gW1MJlNAm0BFURCxurqac15VVZWenh4WFtaM9RcVY4wtWbJk+/btJpOJFlFBQTR3uD2Fb7fbL126VFVVlZGRYTabm7FW3oNp5gQG7OixozNmzGCMpaenh4SEMOb2WLyuC8QgWrdu3XvvvZebmxtoDpg+fbp8dFikDwoKmjt37rhx47766qvPP//8ueeea3Hz1zXXvBtwfdcEVHp1dfWKFSvmzZt35swZ3duWtmRpEZDX4ogdOnQAAMZg0qRJeL2JdmVFkx2isvISYyZmUoBBWlqaIdXVsr9erhWBc6yutq5YsWrevPlnzhQ2mv7GBdktOOdLly4FgClTpqjPOa+oqBgzZkxcXJzVahUPr19NLxskc9LS0nQyp5lQx7FWtlTJ7yIiImjAJie/dGMRrRmha7jMYPTKYrGMGTOme/fugsF8aEasWLFi8ODBkrlBAVAYY61bt36gX/9VK1YaP3E4HOXl5eXl5S4zzMrKGjZsWGRkJGUndKwuXbq89NJLV7UTMzMzx4wZExsbK7fFg7pF54PpbMB1gZM1F7HfA/0ZY4mJ+VdynAAAIABJREFUiS7ToM4k7Gx8z/4uJyEhgdqVkpJSVlbmuegBAwYAQGRk5KVLl5qnMc2HFnouQmeQb3aRLTKMioqiYTN06FDdq1tqnvC2sRx5fcPuVXx8PE2rU6dOxeutB18xqM7h4RHAABR2LfXgRmtlRFRUDGMmAGXYsKeucZWuKnTtpT27wsJCf39/ADhx4oQuvd1u9/B5o8+vBrwsKyoqiiaPYcOGIaLD0bysVYdY5+6dGLBTpqTciEPVHa64LYWFhWSWO378uMiKcrPZbE3J2QcP4Jzb7fbx48eruiMzPfHEExUVFbKqJTbiEXH27Nk0ZMjKa8yNfixatIhm89jY2AMHDly79iAuXLhQ6MFdunTZv39/bm5ubm7uwYMH6d/s7Ow5c+ZQKwYPHnwt6yZDnELhnP/fyhWgMEVRBPPLaRr+NJyIsNlsr6a8BgpjjPXq1Ss7O1t32tBl0Rs2bKDemTNnztVoWlNwne9JoMEKb7fbz549a7xiaUzpJWw2W1FRka44yt9ms509exYRGWPkOxqdL5xeWYk3Iph0XcElodVoMQzQ1HAFh6YQRGjVqhU0raeu4CtvnEC5zJZYX/xJzQkMFLtUNCi48CeKyBCvIjPoKkkc65L97Hb72bOFiA7GsKzsF10OIp8rHiy6WjUlH12VLhcKMIYw94M/19rtYWFh3bp10yXw8/MDw4g2Qr4Q2cQqufvWHR8aEyMiyRwaKeS0obkvq7VCNLksGhHFgPXzMzdRvqFhk/Rqw10pRFIvqyHkGxFk7ty5NpstLCyse/fulECQhRisUf9cV1D5Jo6sqw2X1XNZYbvdXlhY6CGBOzDG/Pz8Xn55kvonQFJSUrt27UD1G8EROAMGCIwxq9X64YcfUsqffvrJpc87+tGhQwdi8u3bt/ft29f7+jQRnPOkpCTQxvLdd9+dkJAQHx8fHx/fR0P//v1nzpyZkpICAF26dLlmddOBPPMwBOD49v++xRB/PeKZbt27OwAczmkaPmGMA9KFRgTMO5l/d5/eH3/4EXCcNGnSgQMH+vXrJ/vBcCdYhg8ffueddwJAWloaLTJbzihoWfeFf/7554cffviDDz4QTxqd5xrNcNCgQX/+859dZmI2m4cPH06/iY/BeeLUoeV0W9NxxW1h0hgRJ5zatm3bkMC7nnKpjF5WDYUC4SGla19p7rq44eqrUzVk14lXA3Lm586dGzRo0J/+9CeXKf39/YljEVFwrMjBOGc3hWOb2OrLdR7iMuXq1asBgJRgl23RezVx5yJaqkyzd6XxSBy68ffEGCOZQ8w/atSo5q2JDEEKMs9Qb4oB26ZNmybmTK3jmkd9+a1R52sWyemy42QJfwU9u3LlStAYTFdQ8/KJkDmyPGmZE4qXA//cuXMDBw6cO3eueHJZzUFE0op0ZTn9hQgAH3744cWLFynZ+fPn5cS6JegXX3wBADNnzrzGiqYYWVQfD3Qg9yN0POl6grHMzMzTp88gwmuvveY5LSKqLvYBd+/efd999+Xn5wPA0qVLP/nkE39/f53V0kNWL730EgCUlZWtXr26KcbNZkcL0oOXL1/eo0ePf/3rX8HBweKhZ2J5puPy5ct79uypy1D34dq1a//5z3/u27dvxowZ0Ni0eitYiN15cmXSfwLajHgZElBe2DQYoZtmW3I5qRhnX/lJQ+magahNmzYMgSEw5AxAnFq7Np1OFfv0009pCLi72IGIGRkZGzZsyM7OnjFjhpFWYmK4Sjpfs8Od5QkBcr8/aKmsQIDbb7/daGbwRi1udhg5zaWI8FyNNWvWUA/+/ve/v0qVlPlW+8EBuFSxy4id5iGlcS2qG9e6yjQvhIS/grutjLHc3NyKigoGcHvYbUZZofsBzcFdOlI0I1lQgvyk6dnSD52sJkn173//W1DenTz3gDZt2jEtwktFRUWDQREUBgqdUqutr5s7L1X9gMFP536mzDl5fVcUUVZBQcG6devCwsKmT59+ja3vUj8qDEyBgYHuSuzVqxcAiKPM1w0MFn+4BBE6der0+JDHGIAJwKS91O0JU
Matias Ezequiel Suarez
Feb 26, 2024, 11:22:09 AM2/26/24
to Leonardo Jiménez, Will H, wrfpython-talk, Marco Miani
These are the python lines that calculate the wind gust from a wrfout.
frame = wrf.ALL_TIMES
wspd_wdir10 = wrf.getvar(wrfout, 'wspd_wdir10', timeidx=frame)
wspd_wdir = wrf.getvar(wrfout, 'wspd_wdir', timeidx=frame)
pblh = wrf.getvar(wrfout, 'PBLH', timeidx=frame)
height = wrf.getvar(wrfout, 'height_agl', timeidx=frame)
u10 = wrf.getvar(wrfout, 'U10', timeidx=frame)
v10 = wrf.getvar(wrfout, 'V10', timeidx=frame)
wrfout.close()
wspd_pblh = wrf.interplevel(wspd_wdir, height, pblh)Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â # 1.
g_pblh = np.where(pblh <= 1000, pblh, 1000)Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â # 2.Â
G = (wspd_wdir10[0] + (wspd_pblh[0] - wspd_wdir10[0])*(1-(g_pblh/2000)))*3.6Â Â # 3. gust in km/h
I have attached the complete script to this email.
_________________________________________________
MatÃas
Ian Dragaud
Feb 27, 2024, 4:03:56 AM2/27/24
to wrfpython-talk
Hi,
Adding the diagnostic max wind variable to the WRF model output is an interesting option.
It was also proposed here: https://forum.mmm.ucar.edu/threads/gust_wind-implementation.9917/
The WRF model User’s Guides details how to use it. Look for output_diagnostics.
Best,
Ian
Adding the diagnostic max wind variable to the WRF model output is an interesting option.
It was also proposed here: https://forum.mmm.ucar.edu/threads/gust_wind-implementation.9917/
The WRF model User’s Guides details how to use it. Look for output_diagnostics.
Best,
Ian
--
You received this message because you are subscribed to the Google Groups "wrfpython-talk" group.
To unsubscribe from this group and stop receiving emails from it, send an email to wrfpython-tal...@ucar.edu.
To view this discussion on the web visit https://groups.google.com/a/ucar.edu/d/msgid/wrfpython-talk/CAE2oGYS5sgAjYo8wuzWyuzryo%2B6fXHw79aF2rHx20CtyVK0Ncg%40mail.gmail.com.
Will H
Feb 27, 2024, 7:13:46 AM2/27/24
to wrfpython-talk, Marco Miani, Matias Ezequiel Suarez, leonardoj...@gmail.com, Ian Dragaud
First of all a bit thank you to Matias for figuring this out ahead of time. I have been trying to get this to work for months now and I see where my mistake was made.
I have this issue raised in the wrf forum
https://forum.mmm.ucar.edu/threads/found-this-burried-deep-in-the-forum-did-wrf-arw-ever-get-the-gust-calculation.15331/
If everyone can show their support for it then perhaps it will get implemented? Add some comments about why this is important to you and how you would use the gust values in WRF perhaps?
I have this issue raised in the wrf forum
https://forum.mmm.ucar.edu/threads/found-this-burried-deep-in-the-forum-did-wrf-arw-ever-get-the-gust-calculation.15331/
If everyone can show their support for it then perhaps it will get implemented? Add some comments about why this is important to you and how you would use the gust values in WRF perhaps?
Matias Ezequiel Suarez
Feb 27, 2024, 10:35:49 AM2/27/24
to Will H, wrfpython-talk, leonardoj...@gmail.com, Marco Miani
Hello Will, here in Córdoba, Argentina, we are using the following parameterization to calculate wind gusts from WRF output:
If you are interested, I can share the Python lines to perform this calculation.
To view this discussion on the web visit https://groups.google.com/a/ucar.edu/d/msgid/wrfpython-talk/85ce9bb5-6f64-4ecf-a2f0-8d9042e66e57n%40ucar.edu.
Will H
Feb 28, 2024, 12:07:17 PM2/28/24
to wrfpython-talk, Will H, Marco Miani, suarezm...@gmail.com, leonardoj...@gmail.com, iandr...@lamma.ufrj.br
https://forum.mmm.ucar.edu/threads/found-this-burried-deep-in-the-forum-did-wrf-arw-ever-get-the-gust-calculation.15331/#post-39728
Might be worth adding your comments as to why WRF needs to add this?
Might be worth adding your comments as to why WRF needs to add this?
Will H
Jul 25, 2024, 7:58:56 AM7/25/24
to wrfpython-talk, Will H, Marco Miani, suarezm...@gmail.com, leonardoj...@gmail.com, iandr...@lamma.ufrj.br, Gabriel Cassol
Has anyone done this in GrADS?
Will H
Nov 8, 2024, 5:58:59 AM11/8/24
to wrfpython-talk, Will H, Marco Miani, suarezm...@gmail.com, leonardoj...@gmail.com, iandr...@lamma.ufrj.br, gabriel...@gmail.com
I made an error here's the source material and the new code.
Code:
SUBROUTINE CALGUST(LPBL,ZPBL,GUST)
C$$$ SUBPROGRAM DOCUMENTATION BLOCK
C . . .
C SUBPROGRAM: CALGUST COMPUTE MAX WIND LEVEL
C PRGRMMR: MANIKIN ORG: W/NP2 DATE: 97-03-04
C
C ABSTRACT:
C THIS ROUTINE COMPUTES SURFACE WIND GUST BY MIXING
C DOWN MOMENTUM FROM THE LEVEL AT THE HEIGHT OF THE PBL
C
C
C PROGRAM HISTORY LOG:
C 03-10-15 GEOFF MANIKIN
C 05-03-09 H CHUANG - WRF VERSION
C 05-06-30 R ROZUMALSKI - DYNAMIC MEMORY ALLOCATION AND SMP
C THREAD-SAFE VERSION
C
C USAGE: CALL CALGUST(GUST)
C INPUT ARGUMENT LIST:
C NONE
C
C OUTPUT ARGUMENT LIST:
C GUST - SPEED OF THE MAXIMUM SFC WIND GUST
C
C OUTPUT FILES:
C NONE
C
C SUBPROGRAMS CALLED:
C UTILITIES:
C H2V
C
C LIBRARY:
C COMMON -
C LOOPS
C OPTIONS
C MASKS
C INDX
C
C ATTRIBUTES:
C LANGUAGE: FORTRAN 90
C MACHINE : CRAY C-90
C$$$
C
C
use vrbls3d
use vrbls2d
C
C INCLUDE ETA GRID DIMENSIONS. SET/DERIVE PARAMETERS.
C
!
INCLUDE "params"
C
INCLUDE "CTLBLK.comm"
C
C DECLARE VARIABLES.
C
INTEGER :: LPBL(IM,JM)
REAL :: GUST(IM,JM)
REAL ZPBL(IM,jsta_2l:jend_2u)
C
C
C*****************************************************************************
C START CALMXW HERE.
C
C LOOP OVER THE GRID.
C
DO J=JSTA,JEND
DO I=1,IM
! GUST(I,J) = SPVAL
GUST(I,J) = 0.
ENDDO
ENDDO
C
C ASSUME THAT U AND V HAVE UPDATED HALOS
C
!$omp parallel do
!$omp& private(ie,iw,mxww,u0,v0,wind)
DO 20 J=JSTA_M,JEND_M
DO 20 I=2,IM-1
L=LPBL(I,J)
IF(MODELNAME .EQ. 'NMM')THEN
IE=I+MOD(J+1,2)
IW=I+MOD(J+1,2)-1
USFC=D25*(U10(I,J-1)+U10(IW,J)+
X U10(IE,J)+U10(I,J+1))
VSFC=D25*(V10(I,J-1)+V10(IW,J)+
X V10(IE,J)+V10(I,J+1))
SFCWIND=SQRT(USFC**2 + VSFC**2)
U0 = D25*(U(I,J-1,L)+U(IW,J,L)+
X U(IE,J,L)+U(I,J+1,L))
V0 = D25*(V(I,J-1,L)+V(IW,J,L)+
X V(IE,J,L)+V(I,J+1,L))
WIND=SQRT(U0**2 + V0**2)
ELSE IF(MODELNAME .EQ. 'NCAR')THEN
USFC=U10(I,J)
VSFC=V10(I,J)
SFCWIND=SQRT(USFC**2 + VSFC**2)
U0=U(I,J,L)
V0=V(I,J,L)
WIND=SQRT(U0**2 + V0**2)
END IF
DELWIND=WIND - SFCWIND
ZSFC=FIS(I,J)*GI
DELWIND=DELWIND*(1.0-AMIN1(0.5,ZPBL(I,J)/2000.))
GUST(I,J)=SFCWIND+DELWIND
10 CONTINUE
20 CONTINUE
C END OF ROUTINE.
C
RETURN
END
Key References​
Python:
# Calculate the 10-meter wind speed
wspd_wdir10 = np.sqrt((u10**2) + (v10**2))
wspd_wdir = np.sqrt((u**2) + (v**2))
pblh = wrf.getvar(ncfile, 'PBLH', meta=True)
height = wrf.getvar(ncfile, 'height_agl', meta=True, units="m")
# Initialize the maximum gust speed array (same shape as surface wind speed)
surface_gust_speed = np.zeros_like(wspd_wdir10) # Assuming wspd_wdir10 is the surface wind speed
# Iterate over each vertical level
for k in range(u.shape[0]):
# Calculate the wind speed at level k
wind_speed_k = np.sqrt(u[k]**2 + v[k]**2)
height_k = np.minimum(height[k], 1000) # Cap height at 1000 meters
# Calculate the weight based on capped height
weight_k = (1 - (height_k / 2000)) # For heights > 1000m, weight_k will be 0.5
# Adjust the wind speed by the weighted difference from the surface
surface_gust_speed = wspd_wdir10 + ((wind_speed_k - wspd_wdir10) * weight_k)
# Apply Gaussian smoothing and convert to mph
gusts = gaussian_filter(surface_gust_speed * 2.23694, sigma=1) # Convert m/s to mph
Will H
Nov 8, 2024, 10:20:33 AM11/8/24
to wrfpython-talk, Will H, Marco Miani, suarezm...@gmail.com, leonardoj...@gmail.com, iandr...@lamma.ufrj.br, gabriel...@gmail.com
Python:
  # Load the dataset
  ncfile = Dataset(ncfile_path)
  # Extract atmospheric variables
  # Fetch sea level pressure data from the netCDF file, which may be used in other analyses or conditions
  slp = wrf.getvar(ncfile, "slp", meta=True) Â
  # Fetch the 10-meter wind components (U10 for x-direction, V10 for y-direction)
  # These represent the horizontal wind speeds at 10 meters above ground level
  u10 = wrf.getvar(ncfile, "U10", meta=True)
  v10 = wrf.getvar(ncfile, "V10", meta=True)
  # Fetch the upper-level wind components (ua for x-direction, va for y-direction) at multiple vertical levels
  # These represent the horizontal wind speeds at various heights in the atmosphere
  u = wrf.getvar(ncfile, "ua", meta=True, units="m s-1")
  v = wrf.getvar(ncfile, "va", meta=True, units="m s-1")
  # Convert 10-meter wind components from meters per second to knots
  # This is done by multiplying by the conversion factor 1.94384449
  u10_knots = u10 * 1.94384449
  v10_knots = v10 * 1.94384449
  # Calculate the total 10-meter wind speed by combining u10 and v10 components using the Pythagorean theorem
  wspd_wdir10 = np.sqrt((u10**2) + (v10**2))
 Â
  wspd_wdir = np.sqrt((u**2) + (v**2))
  # PBLH represents the top height of the atmospheric boundary layer
  pblh = wrf.getvar(ncfile, 'PBLH', meta=True)
  # This provides information on the altitude of each vertical level
  height = wrf.getvar(ncfile, 'height_agl', meta=True, units="m")
  # This corresponds to vPBL in the wind gust parameterization formula, using interpolation
  # According to RUC20, the capped height should be 1000m for heights above 1000m [20:20†source]
  capped_pblh = np.minimum(pblh, 1000)
  wspd_pblh = wrf.interplevel(wspd_wdir, height, capped_pblh)
  # Cap the height at the minimum of either the planetary boundary layer height (PBLH) or 1000 meters
  # According to the RUC20 parameterization used by NOAA, the capped height should be 1000m for heights above 1000m [20:20†source]
  height_k = np.minimum(pblh, 1000)
  # Calculate the weight based on the capped height
  # This weight decreases linearly with height, reducing the influence of wind at higher levels
  weight_k = (1 - (height_k / 2000))  # For heights ≥1000m, weight_k becomes 0.5 [20:20†source]
  # Adjust the surface gust speed based on the parameterization formula
  # This combines the surface wind speed (10-meter level) with the weighted difference between
  # the wind speed at the top of the PBL and the surface wind speed
  surface_gust_speed = wspd_wdir10 + ((wspd_pblh - wspd_wdir10) * weight_k)
  # Apply Gaussian smoothing to the resulting gust speed field to smooth out abrupt changes in gust speed values
  # Convert the smoothed gust speeds from meters per second to miles per hour by multiplying by 2.23694
  gusts = gaussian_filter(surface_gust_speed * 2.23694, sigma=1)  # Convert m/s to mph
  # gusts as a function of the 10-meter wind speed and the wind at the top of the planetary boundary layer.
  # Reference: Vindel et al. (2020), Eq. (6), which uses the form W_G = v10 + (vPBL - v10) * (1 - hPBL/2000) [20:20†source].
Will H
Nov 21, 2024, 6:12:21 AM11/21/24
to wrfpython-talk, Will H, Marco Miani, suarezm...@gmail.com, leonardoj...@gmail.com, iandr...@lamma.ufrj.br, gabriel...@gmail.com
solve the problem with help from someone on my linkedin
If pblh < surface wind then it gives a nan. So I had to adjust the script for nan calculations
Reply all
Reply to author
Forward
0 new messages