SLiMv2.1: SLiMv2.1/mkspot.f90 Source File

00001 
00002 
00003 program makespot
00004   !copyright Robert Cameron, 2008
00005   !IDL original version makespot.pro by RHC & some stuff by HS
00006   implicit none
00007   character*80 version
00008   parameter(version="v2.2a")
00009   ! UPDATE EACH TIME THE IS MODIFIED
00010   ! version of this code
00011  
00012   !;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
00013   ! PURPOSE: to calculate a magnetic field, velocity profile
00014   ! for a solar like atmosphere given the below parameters
00015   !
00016   ! IMPORTANT: There is a long list of possible variables. These
00017   !            variables can be changed in the section immediately
00018   !            following these comments.
00019   !
00020   ! INPUT:  what - 'plage': implements plage with the sunspot
00021   !                'fast': sets magnetic field equal to zero and
00022   !                implements the fast sound speed in it's place
00023   !                'cspert': sets magnetic field to zero and implements
00024   !                a sound speed perturbation.
00025   !                'spot'
00026   !                If left blank then uses quiet sun
00027   !
00028   !RHC   ---I don't totally like the above concept.
00029   !
00030 
00031   ! OUTPUT: FITS files of magnetic field (bx, by, bz), flows (ux, uy,
00032   !         uz), rho, pressure, gamma_1
00033   !         !WARNING! This program assumes there is a directory ../spot
00034   !
00035   !
00036   ! MODIFICATION HISTORY
00037   !
00038 
00039   integer nx,ny,nz,i,j,k,ztop,i1
00040   parameter(nx=200,ny=1,nz=1100)
00041   character*80 what,fn
00042   real dx,dy,dz,tx,xc,yc
00043   real flux,b_0,alpha,h,dh,oc
00044   real cs_factor,plage_field
00045   real depth_r
00046   real bx_out(nx,ny,nz),by_out(nx,ny,nz),bz_out(nx,ny,nz)
00047   real p_out(nx,ny,nz),rho_out(nx,ny,nz),gamma_out(nx,ny,nz)
00048   real cs_out(nx,ny,nz),g_out(nx,ny,nz)
00049   real ux_out(nx,ny,nz),uy_out(nx,ny,nz),uz_out(nx,ny,nz)
00050 
00051 
00052   real rr(2481),rho(2481),P(2481),Cs(2481),g_rhc(2481)
00053   real rho_s(nz),P_s(nz),Cs_s(nz),g_s(nz)
00054   real r_r(2481),r_rho(2481),r_g_rhc(2481),msun_r
00055   real bz(nz),dbzdz(nz),d2bzdz2(nz),pb(nx,ny,nz),p_st(nx,ny,nz)
00056 
00057   real rho_st(nx,ny,nz),dp_st(nx,ny),dp(nx,ny)
00058   real m3(nx,ny,nz),mz3(nx,ny,nz)
00059   real u3(nx,ny,nz),ux3(nx,ny,nz),uy3(nx,ny,nz)
00060   real b2(nx,ny,nz),delta_cs2(nx,ny,nz),oordrurdz3(nx,ny,nz)
00061   real cs_o(nx,ny,nz),cs_new(nx,ny,nz)
00062 
00063   real x,y,r,z,hor_r,dr
00064   real sigma,epsilon,z_fac,fac,f,f2,hwhm
00065   real dPdz(nx,ny,nz),p1(nx,ny,nz)
00066 
00067   real G_Univ,pi,g
00068 
00069   real zz(nz),bzt(nz),frac_radius(nz)
00070   real total_b,r_sun
00071 
00072   integer ns,k1
00073   interface
00074      function vr(r,z)
00075        implicit none
00076        real  vr,r,z
00077      end function vr
00078   end interface
00079   interface
00080      function reverse(a)
00081        implicit none
00082        integer nx,ny,nz
00083        parameter (nx=200,ny=1,nz=1100)
00084        real  reverse(nz),a(nz)
00085      end function reverse
00086   end interface
00087 
00088   interface
00089      function smooth(a,n)
00090        implicit none
00091        integer nx,ny,nz
00092        parameter (nx=200,ny=1,nz=1100)
00093        real  smooth(nx,ny,nz),a(nx,ny,nz)
00094        integer n
00095      end function smooth
00096   end interface
00097   interface
00098      function interpol(a,b,c)
00099        implicit none
00100        integer nx,ny,nz
00101        parameter (nx=200,ny=1,nz=1100)
00102        real  interpol(nz),a(2481),b(2481),c(nz)
00103      end function interpol
00104   end interface
00105   what="spot"
00106   !;;;;      BEGIN VARIABLE PARAMETERS      ;;;;;;;;;;;;;;;
00107   cs_factor=0.96   ! used to alter the sound speed profile of model S
00108   plage_field=00. ! gauss
00109 
00110   !-----misc
00111   tx=145.741e8  ! total extent in the x direction - cm
00112   yc=0          !in cm  
00113   xc=40e8
00114   ztop=100
00115 
00116   !----Spot parameters
00117   ! SUNSPOT
00118   flux=4.41e22                ! gauss cm^2 for a radius of 15 Mm
00119   b_0=3000.                   ! gauss
00120   alpha=6.25e8                ! controls inclination/spread of field
00121   h=-4.e9                 ! disconnection depth - cm
00122   dh=0.5e8                    ! disconnection range - cm
00123 
00124   !basic physics
00125   pi=3.1415926535
00126   G_Univ= 6.6725985e-11*1e6/1e3  !m3/(kg s2) 
00127   g=-2.74e4   ! gravity -BUT IT IS VARIABLE IN THE WAVE PROPAGATION
00128   R_sun=6.96e10
00129 
00130 
00131   dx=tx/float(nx)  ! cm/pixel
00132   dy=dx*2          ! cm/pixel
00133   dz=12.5e8/500.   ! cm/pixel
00134 
00135 
00136   oc=1. !control for Va saturation
00137   !;;;END VARIABLES;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
00138 
00139   !============START CODE PROPER==============
00140   !---read in Model S data-----------
00141 
00142   open(11,file="Tables/R_SSM.dat")
00143   read(11,*) rr
00144   close(11)
00145   rr=rr*100                    
00146 
00147   open(11,file="Tables/rho_SSM.dat")
00148   read(11,*) rho                 !g/cm^3   
00149   close(11)                        
00150 
00151   open(11,file="Tables/Cs_SSM.dat")
00152   read(11,*) Cs
00153   close(11)
00154   Cs=Cs*cs_factor                ! cm/s
00155 
00156   open(11,file="Tables/P_SSM.dat")
00157   read(11,*) P
00158   close(11)                      ! g cm/s /cm^2  
00159 
00160 
00161   !Implementing a variable gravity
00162   g_rhc(1)=0 
00163   do i=1,2481  
00164      r_rho(i)=rho(2481+1-i)
00165      r_r(i)=r_sun-rr(2481+1-i)
00166   enddo
00167   msun_r=0
00168   do i=2,2480
00169      dr=(r_r(i+1)-r_r(i-1))/2.
00170      msun_r=msun_r+r_r(i)**2*4*pi*dr*r_rho(i)
00171      g_rhc(i)=msun_r*G_univ/r_r(i)**2
00172   enddo
00173   g_rhc(2481)=msun_r*G_univ/r_r(2481)**2
00174 
00175   ! depth
00176   do i=1,nz 
00177      zz(i)=(i-ztop)*dz   ! cm
00178   enddo
00179 
00180   ! interpolate the Model S onto our depth profile
00181   rho_s=exp(interpol(alog(rho),rr,zz))
00182   Cs_s=(interpol((cs),rr,zz))
00183   P_s=exp(interpol(alog(p),rr,zz))
00184   do i=1,2481
00185      r_g_rhc(i)=g_rhc(2481+1-i)
00186   enddo   
00187   g_s=interpol(r_g_rhc,rr,zz)
00188 
00189   do i=1,nx 
00190      do j=1,ny 
00191         g_out(i,j,:)=reverse(g_s)
00192      enddo
00193   enddo
00194 
00195 !  write(*,*) rho_s
00196 
00197   do i=1,ztop+4 
00198      i1=ztop+4-i
00199      rho_s(i1)=rho_s(ztop+5)*exp(-(i+1)/18.5*exp((i+1)/100.))
00200      ! next is strange but seems to work (try different things?)
00201      cs_s(i1)=cs_s(i1)*exp(-i/50.)+(1-exp(-i/50.))*cs(ztop+5)
00202      p_s(i1)=p_s(ztop+5)*exp(-(i+1)/8.5*exp((i+1)/100.))          
00203   enddo
00204 
00205   !===============ADD SPOT========================
00206   !--first magnetic field--
00207 
00208   do i=1,nz 
00209      z=zz(nz-1-i)
00210      z_fac=exp((-z-h)/dh)/(1+exp((-z-h)/dh))
00211      bz(i)=b_0*exp(z/alpha)*z_fac
00212   enddo
00213 
00214   ! derivatives of magnetic field
00215 
00216   ! first derivative
00217   dbzdz(1)=(bz(2)-bz(1))/(zz(2)-zz(1)) 
00218   dbzdz(2:nz-1)=(bz(3:nz)-bz(1:nz-2))/(zz(3:nz)-zz(1:nz-2)) 
00219   dbzdz(nz)=(bz(nz)-bz(nz-1))/(zz(nz)-zz(nz-1)) 
00220 
00221   ! second derivative
00222   d2bzdz2(2:nz-1)=(bz(1:nz-2)-2.*bz(2:nz-1)+bz(3:nz)) &
00223                      &         /(zz(2:nz-1)-zz(3:nz))**2
00224   d2bzdz2(1)=2.*d2bzdz2(1)-d2bzdz2(2)
00225   d2bzdz2(nz)=d2bzdz2(nz-1)
00226 
00227 
00228   ! calculate bx, by and bz
00229   ! the three different loops here (k=0,nz-1) ensure the 
00230   ! periodicity of the box
00231 
00232   do k=1,nz   
00233      do i=1,nx 
00234         do j=1,ny 
00235            x= (i-1)*dx-xc      ! 40e8
00236            y= (j-1)*dy-yc
00237            r=sqrt(x**2+y**2)
00238            f=exp(-bz(k)*r**2/(flux/pi))
00239            f2=exp(-2*bz(k)*r**2/(flux/pi))
00240            bz_out(i,j,k)=bz(k)*f           
00241            bx_out(i,j,k)=-x/2.*dbzdz(k)*f
00242            by_out(i,j,k)=-y/2.*dbzdz(k)*f
00243         enddo
00244      enddo
00245   enddo
00246 
00247   do k=1,nz  
00248      do i=1,nx 
00249         do j=1,ny 
00250            x= (i-1)*dx-xc+tx    !40e8+tx
00251            y= (j-1)*dy-yc
00252            r=sqrt(x**2+y**2)
00253            f=exp(-bz(k)*r**2/(flux/pi))
00254            f2=exp(-2*bz(k)*r**2/(flux/pi))
00255            bz_out(i,j,k)=bz_out(i,j,k)+bz(k)*f           
00256            bx_out(i,j,k)=bx_out(i,j,k)-x/2.*dbzdz(k)*f
00257            by_out(i,j,k)=by_out(i,j,k)-y/2.*dbzdz(k)*f
00258         enddo
00259      enddo
00260   enddo
00261 
00262   do k=1,nz 
00263      do i=1,nx 
00264         do j=1,ny 
00265            x= (i-1)*dx-xc-tx   !40e8-tx
00266            y= (j-1)*dy-yc
00267            r=sqrt(x**2+y**2)
00268            f=exp(-bz(k)*r**2/(flux/pi))
00269            f2=exp(-2*bz(k)*r**2/(flux/pi))
00270            bz_out(i,j,k)=bz_out(i,j,k)+bz(k)*f        
00271            bx_out(i,j,k)=bx_out(i,j,k)-x/2.*dbzdz(k)*f
00272            by_out(i,j,k)=by_out(i,j,k)-y/2.*dbzdz(k)*f
00273         enddo
00274      enddo
00275   enddo
00276 
00277   !;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
00278   ! this ensures the flux at each height remains the same. 492 is a
00279   ! random height to compare with. Possibly unnecessary, but a good
00280   ! check.
00281   ! CAUTION: this will only work in 3D and if done in 2D with a sunspot
00282   ! will make the magnetic field gaussian settle at a higher value than
00283   ! zero
00284 
00285   !2 dimensional code
00286   !r_n=fltarr(200,560)
00287   !for i=0,199 do r_n(i,*)=abs(i-99.5)
00288   !bzt=total(reform((bz_out)*(r_n)),1)/mean(r_n)
00289   !3d
00290 
00291   !if ny gt 1 and what ne 'qs' then begin
00292   bzt=0
00293   do k=1,nz
00294      do i=1,nx
00295         do j=1,ny
00296            bzt(k)=bzt(k)+bz_out(i,j,k)
00297         enddo
00298      enddo
00299   enddo
00300 
00301   do k=1,nz 
00302      bz_out(:,:,k)=bz_out(:,:,k)+(bzt(1000)-bzt(k))/float(nx*ny)
00303   enddo
00304   !endif
00305 
00306   !;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
00307 
00308   ! magnetic energy density
00309   B2=(bx_out*bx_out+by_out*by_out+bz_out*bz_out)/(8.*pi)
00310 
00311   ! scaling factor for the Lorentz force above the surface in the
00312   ! low-beta region
00313   pB=B2
00314   do k=1,nz 
00315      pB(:,:,k)=pB(:,:,k)/p_s(nz-1-k)
00316   enddo
00317   do k=nz-ztop+10,nz-1 
00318      pB(:,:,k)=pB(:,:,k)*exp((k-(nz-ztop+10))/20.)
00319   enddo
00320 
00321   !-----Sunspot's Pressure--
00322   do k=1,nz 
00323      do i=1,nx 
00324         do j=1,ny 
00325            x= i*dx-xc
00326            y= j*dy-yc
00327            fac=1-(pB(i,j,k)/oc)/(1+(pB(i,j,k)/oc))
00328            r=sqrt(x**2+y**2)
00329            f2=exp(-2*bz(k)*r**2/(flux/pi))*fac
00330            p_out(i,j,k)=p_s(nz-1-k)+f2/(64.*pi)* &
00331             &    (-8.*Bz(k)**2+2.*d2Bzdz2(k)*flux/pi &
00332             &    -dbzdz(k)**2/bz(k)*(flux/pi+2.*bz(k)*r**2))
00333            p_out(i,j,k)=max(p_out(i,j,k),p_s(nz-1-k)/10.)       
00334         enddo
00335      enddo
00336   enddo
00337 
00338   !----sunspot's density------
00339   dPdz(:,:,1)=(p_out(:,:,2)-p_out(:,:,1))/(zz(2)-zz(1)) 
00340   dPdz(:,:,2:nz-1)=(p_out(:,:,3:nz)-p_out(:,:,1:nz-2))/(zz(3)-zz(1)) 
00341   dPdz(:,:,nz)=(p_out(:,:,nz)-p_out(:,:,nz-1))/(zz(nz)-zz(nz-1)) 
00342 
00343   rho_st=rho_out
00344   do k=1,nz 
00345      do i=1,nx
00346         do j=1,ny
00347            x= i*dx-xc
00348            y= j*dy-yc
00349            r=sqrt(x**2+y**2)
00350            f2=exp(-2*bz(k)*r**2/(flux/pi))
00351            fac=1-(pb(i,j,k)/oc)/(1+(pb(i,j,k)/oc))
00352            rho_st(i,j,k)= 1./(-g_out(i,j,k))*(dPdz(i,j,k)  -         &  
00353    &            1./(4.*pi)*(r/2.) *dBzdz(k)*f2*fac*    &  
00354    &  (r*(-d2Bzdz2(k)/2.+2*bz(k)*pi/flux) +r**3/2.*dbzdz(k)**2*pi/flux))
00355         enddo
00356      enddo
00357      rho_st(:,:,1)=rho_st(:,:,2)**2/rho_st(:,:,3) 
00358                                        !extrapolation in log space
00359   enddo
00360 
00361   rho_out=rho_st
00362 !  write(*,*) rho_st
00363   rho_out(:,:,3)=rho_out(:,:,4)*rho_out(:,:,4)/rho_out(:,:,5)
00364   rho_out(:,:,2)=rho_out(:,:,3)*rho_out(:,:,3)/rho_out(:,:,4)
00365   rho_out(:,:,1)=rho_out(:,:,2)*rho_out(:,:,2)/rho_out(:,:,3)
00366 
00367   !CHECKED THAT DENSITY IS CORRECTLY REPRODUCED !OK SEPT 2007!
00368 
00369   !------sunspot's soundspeed-------
00370 
00371   do k=1,nz 
00372      k1=nz-1-k 
00373      cs_out(:,:,k)=sqrt(Cs_s(k1)**2)
00374   enddo
00375 
00376   !===============STABILIZE========================
00377   !-------------Keep rho and Cs--------------------
00378   !--reconstruct P  ---------
00379 
00380   p1=p_out
00381   P_st=P_out*1.! Only the first two of these values are kept, 
00382                ! and go directly into P_st, 
00383                 !so can change this  arbitrarily
00384 
00385   do i=nz-2,1,-1 
00386      dp=(P_out(:,:,i+1)-P_out(:,:,i-1))/2.  ! delta pressure
00387      dp_st=+Cs_out(:,:,i)**2*(rho_out(:,:,i+1)-rho_out(:,:,i-1))/2.*0.98
00388      !delta_rho*cs - what is 0.98 for? 
00389      !value of gamma? no.....stumps for epsilon?
00390      dp=max(dp, dp_st)
00391      P_st(:,:,i-1)=P_st(:,:,i)+dp
00392   enddo
00393 
00394   p_st(:,:,1)=p_st(:,:,2)*(p_st(:,:,2)/p_st(:,:,3))
00395   p_out=-p_st+minval(p1)-minval(-p_st)
00396   !p_out(*,*,nz-1)=p_out(*,*,nz-2)*(p_out(*,*,nz-2)/p_out(*,*,nz-3))
00397 
00398   !===========Construct Gamma====================
00399 
00400   ns=155   ! a place in the relatively quiet model - far form the spot
00401   P_out=P_out
00402   do k=1,nz
00403      gamma_out(:,:,k)=cs_out(:,:,k)**2*rho_out(ns,1,k)/P_out(ns,1,k)
00404   enddo
00405 
00406   !================================================
00407   !smooth
00408   p_out=exp(smooth(alog(P_out),2))
00409   rho_out=exp(smooth(alog(rho_out),2))
00410   gamma_out=smooth(gamma_out,2)
00411 
00412   cs_o=sqrt(gamma_out*P_out/rho_out)
00413 
00414 !!!!!!switch off magnetic field and set soundspeed=fast mode speed!!!!!!
00415   if (what .eq. 'fast') then 
00416      write(*,*) "Switching off magnetic field"
00417      write(*,*) "setting sound speed = fast mode speed"
00418 
00419      !Next piece of code increases P_out so as to make Cs=fast mode speed
00420      !then switches off magnetic field: the idea is to have "equivalent
00421      !speed"
00422 
00423      b2=bx_out**2+by_out**2+bz_out**2
00424 
00425      do k=1010,1100 
00426         b2(:,:,k)=b2(:,:,k)*(1./(1+exp((k-1020.)/4.)))
00427      enddo
00428 
00429      cs_new=max(sqrt(b2/(4*pi*rho_out)+cs_o**2),(5*cs_o))
00430 
00431      P_out=(cs_new/cs_o)**2*P_out
00432 
00433      ! set the magnetic field to zero
00434      Bx_out=Bx_out*1.e-6
00435      By_out=By_out*1.e-6
00436      Bz_out=Bz_out*1.e-6
00437   endif
00438 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00439 !!!!!!insert sound speed perturbation!!!!!!!!!!!!!!
00440   if (what .eq. 'cspert') then 
00441      ! the only thing that's written out is gamma1 which uses cs, 
00442      ! so if you modify  gamma by multiplying by two then you 
00443      !modify cs by sqrt(2)
00444 
00445      !     INSERT CSPERT HERE
00446 
00447      ! set the magnetic field to zero
00448      Bx_out=Bx_out*1.e-6
00449      By_out=By_out*1.e-6
00450      Bz_out=Bz_out*1.e-6
00451 
00452      ! insert a small sound speed perturbation
00453      epsilon=0.1
00454      depth_r=10e8                ! cm
00455      hor_r=20e8                  ! cm
00456 
00457      do k=1,nz  
00458         do i=1,nx 
00459            do j=1,ny 
00460               x= (i-1)*dx-xc          
00461               y= (j-1)*dy-yc
00462               r=sqrt(x**2+y**2)
00463               delta_cs2(i,j,k)=cs_out(i,j,k)**2*epsilon* &
00464          &         exp(-r**2/hor_r**2)*exp(-zz(nz-1-k)**2/depth_r**2)
00465            enddo              
00466         enddo                 
00467      enddo                     
00468  
00469      cs_out = sqrt(cs_out**2 + delta_cs2)
00470 
00471   endif
00472 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00473 
00474   !add in the moat flow
00475 
00476   do i=1,nx 
00477      do j=1,ny  
00478         do k=1,nz 
00479            x=(i-1)*dx
00480            y=(j-1)*dy
00481            r=sqrt((x-xc)**2+(y-yc)**2) 
00482            z=(k-(nz-ztop))*dz 
00483            if((r .gt. 15.e8) .and. (r .lt. 34.3e8) .and. (z .lt. 0)) then
00484               m3(i,j,k)=vr(r,z) 
00485               oordrurdz3(i,j,k)=1./r*((r+0.1e8)*vr(r+0.1e8,z) -&
00486                    &  (r-0.1e8)*vr(r-0.1e8,z))/(0.2e8)
00487            endif
00488         enddo
00489      enddo
00490   enddo
00491 
00492   u3=m3/sqrt(rho_out)
00493   oordrurdz3=oordrurdz3/sqrt(rho_out)
00494 
00495   do k=nz-ztop,0,-1 
00496      mz3(:,:,k)=mz3(:,:,k+1)+oordrurdz3(:,:,k)
00497   enddo
00498 
00499   u3=u3/maxval(u3)*250.
00500   do i=1,nx 
00501      do j=1,ny  
00502         x=(i-1)*dx
00503         y=(j-1)*dy
00504         r=sqrt((x-xc)**2+(y-yc)**2) 
00505         ux3(i,j,:)=u3(i,j,:)*(x-xc)/r 
00506         uy3(i,j,:)=u3(i,j,:)*(y-yc)/r 
00507      enddo
00508   enddo
00509 
00510   !set moat flow
00511 
00512 
00513   !  GOTO,SKIPSHIFTING
00514   ! this was originally to create the sunspot in the centre of the array
00515   ! and then shift the array to the correct orientation. Probably useful
00516   ! to keep the periodicity of the box. may use this in the future. Then
00517   ! must change xc=0 at the beginning!!!
00518   !shift into position
00519   !  bx_out=shift(bx_out(*,*,*),-(tx/2.-40e8)/tx*200.,0,0)
00520   !  by_out=shift(by_out(*,*,*),-(tx/2.-40e8)/tx*200.,0,0)
00521   !  bz_out=shift(bz_out(*,*,*),-(tx/2.-40e8)/tx*200.,0,0)
00522   !  jx_out=shift(jx_out(*,*,*),-(tx/2.-40e8)/tx*200.,0,0)
00523   !  jy_out=shift(jy_out(*,*,*),-(tx/2.-40e8)/tx*200.,0,0)
00524   !  jz_out=shift(jz_out(*,*,*),-(tx/2.-40e8)/tx*200.,0,0)
00525   !  p_out=shift(p_out(*,*,*),-(tx/2.-40e8)/tx*200.,0,0)
00526   !  rho_out=shift(rho_out(*,*,*),-(tx/2.-40e8)/tx*200.,0,0)
00527   !  gamma_out=shift(gamma_out(*,*,*),-(tx/2.-40e8)/tx*200.,0,0)
00528   !  SKIPSHIFTING:
00529 
00530   ! PLAGE
00531   if (what .eq. 'plage') then
00532      write(*,*) "Inserting plage"
00533      frac_radius=(6.959e10-reverse(zz))/6.959e10
00534      do i=1,nz 
00535         bz_out(:,:,i)=bz_out(:,:,i)+plage_field/frac_radius(i)**2
00536      enddo
00537      b_0=plage_field
00538      write(*,*) "Plage field has a value of ",plage_field," gauss"
00539   endif
00540 
00541   ! RADIUS - HWHM at ztop
00542   if (what .eq. 'spot') then 
00543      total_b=sum(sqrt(bz_out(:,:,nz-ztop-2)**2 & 
00544                 &   + bx_out(:,:,nz-ztop-2)**2 &
00545                 &   + by_out(:,:,nz-ztop-2)**2 ))
00546      sigma=sqrt(flux*pi/(4.*sum(bz_out(:,:,nz-ztop-2))))
00547      hwhm=sqrt(2*alog(2.))*sigma
00548   else
00549      hwhm=0.
00550   endif
00551 
00552   ! OUTPUT THE RESULTS!
00553   fn="spot/gamma1_wavesim.fits"
00554   call fw_cube(fn,gamma_out,nx,ny,nz,dx,dy,dz,flux,b_0,alpha,&
00555                hwhm,dh,h,what,version)
00556   fn="spot/bx1_wavesim.fits"
00557   call fw_cube(fn,bx_out,nx,ny,nz,dx,dy,dz,flux,b_0,alpha,&
00558               hwhm,dh,h,what,version)
00559   fn="spot/by1_wavesim.fits"
00560   call fw_cube(fn,by_out,nx,ny,nz,dx,dy,dz,flux,b_0,alpha,&
00561                hwhm,dh,h,what,version)
00562   fn="spot/bz1_wavesim.fits"
00563   call fw_cube(fn,bz_out,nx,ny,nz,dx,dy,dz,flux,b_0,alpha,&
00564                hwhm,dh,h,what,version)
00565   fn="spot/pressure_wavesim.fits"
00566   call fw_cube(fn,p_out,nx,ny,nz,dx,dy,dz,flux,b_0,alpha,&
00567                hwhm,dh,h,what,version)
00568   fn="spot/rho_wavesim.fits"
00569   call fw_cube(fn,rho_out,nx,ny,nz,dx,dy,dz,flux,b_0,alpha,&
00570                hwhm,dh,h,what,version)
00571   fn="spot/ux1_wavesim.fits"
00572   call fw_cube(fn,ux_out,nx,ny,nz,dx,dy,dz,flux,b_0,alpha,&
00573                hwhm,dh,h,what,version)
00574   fn="spot/uy1_wavesim.fits"
00575   call fw_cube(fn,uy_out,nx,ny,nz,dx,dy,dz,flux,b_0,alpha,&
00576                hwhm,dh,h,what,version)
00577   fn="spot/uz1_wavesim.fits"
00578   call fw_cube(fn,uz_out,nx,ny,nz,dx,dy,dz,flux,b_0,alpha,&
00579                hwhm,dh,h,what,version)
00580   fn="spot/g_wavesim.fits"
00581   call fw_cube(fn,g_out,nx,ny,nz,dx,dy,dz,flux,b_0,alpha,&
00582                hwhm,dh,h,what,version)
00583   WRITE(*,*) "SUCCESS"
00584 end program makespot
00585 
00586 subroutine fw_cube(fn,data,nx,ny,nz,dx,dy,dz,flux,b_0,alpha, &
00587                       hwhm,dh,h,what,version)
00588   implicit none
00589   character*80 fn,record,infn
00590   character*30 errmsg
00591   real data(nx,ny,nz)
00592   real dx,dy,dz,time
00593   integer status,unit,blocksize,bitpix,naxis,naxes(3)
00594   integer i,j,group,fpixel,nelements,nx,ny,nz,junk
00595   integer n,nkeys,nspace
00596   logical simple,extend
00597   character*80 what, version
00598   real flux,b_0,alpha,hwhm,dh,h
00599 
00600   group=1
00601   fpixel=1
00602   simple=.true.
00603   status=0
00604   bitpix=-32
00605   naxis=3
00606   naxes(1)=nx
00607   naxes(2)=ny
00608   naxes(3)=nz
00609   extend=.true.
00610   status=0
00611   blocksize=1
00612   nelements=naxes(1)*naxes(2)*naxes(3)
00613   call ftgiou(unit,status)
00614   call ftinit(unit,fn,blocksize,status)
00615   !     write the required header keywords                                 
00616   call ftphpr(unit,simple,bitpix,naxis,naxes,0,1,extend,status)
00617   call ftpkye(unit,'DX',dx,6,'[cm] dx',status)
00618   call ftpkye(unit,'DY',dx,6,'[cm] dy',status)
00619   call ftpkye(unit,'DZ',dx,6,'[cm] dz',status)
00620   call ftpcom(unit,'**MKSPOT.F90 PARAMETERS',status)
00621   call ftpkys(unit,'FEATURE',what,'feature in simulation',status)
00622   call ftpkys(unit,'MKSPOT',version,'version of mkspot.f90',status)
00623   call ftpkye(unit,'FLUX',flux,6,'[gauss cm^2] flux',status)
00624   call ftpkye(unit,'B0',b_0,6,'[gauss] spot strength',status)
00625   call ftpkye(unit,'alpha',alpha,6,'controls spread of field',status)
00626   call ftpkye(unit,'Radius',hwhm,6,'[cm]HWHM of sunspot at surface',status)
00627   call ftpkye(unit,'DH',dh,6,'[cm]disconnection depth',status)
00628   call ftpkye(unit,'HZ',h,6,'[cm]disconnection range',status)
00629   call ftpkye(unit,'HX',0.,6,'[cm]disconnected width',status)
00630   call ftpkye(unit,'VMAX',0.,6,'[cm/s]maximum flow',status)
00631 
00632   !     write the array to the FITS file                                      
00633   call ftppre(unit,group,fpixel,nelements,data,status)
00634 
00635   !     close the file and free the unit number       
00636   call ftclos(unit, status)
00637   call ftfiou(unit, status)
00638   if (status .ne. 0) then
00639      call ftgerr(status,errmsg)
00640      write(*,*) "err: ", errmsg
00641      stop
00642   endif
00643 end subroutine fw_cube
00644 
00645 
00646 
00647 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00648 function vr(r,z)
00649   real vr,r,z
00650   vr=0
00651   return
00652 end function vr
00653 
00654 function smooth(a,n)
00655   integer nx,ny,nz
00656   parameter (nx=200,ny=1,nz=1100)
00657   real smooth(nx,ny,nz),a(nx,ny,nz)
00658   smooth=a
00659   !if (ny .gt. 1) then
00660   !smooth(2:nx-1,2:ny-1,2:nz-1)= &
00661   !                  (a(1:nx-2,2:ny-1,2:nz-1)+a(3:nx,2:ny-1,2:nz-1) &
00662   !                  +a(2:nx-1,1:ny-2,2:nz-1)+a(2:nx-1,3:ny,2:nz-1) &
00663   !                  +a(2:nx-1,2:ny-1,1:nz-2)+a(2:nx-1,2:ny-1,3:nz) &     
00664   !                  +a(2:nx-1,2:ny-1,2:nz-1))/7.
00665   !endif
00666   !smooth(1,:,:)=smooth(2,:,:)
00667   !smooth(nx,:,:)=smooth(nx,:,:)
00668   return
00669 end function smooth
00670 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00671 function reverse(a)
00672   integer nx,ny,nz
00673   parameter (nx=200,ny=1,nz=1100)
00674   real reverse(nz),a(nz)
00675   do i=1,nz
00676      reverse(i)=a(nz+1-i)
00677   enddo   
00678   return
00679 end function reverse
00680 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00681 
00682 
00683 
00684 function interpol(y_old,x_old,x_new)
00685   implicit none
00686 
00687   integer nx,ny,nz
00688   parameter (nx=200,ny=1,nz=1100)
00689   real interpol(nz),x_old(2481),y_old(2481),X_new(nz)
00690   integer i
00691 
00692   interface
00693      function interpol_one(y_old,x_old,xin)
00694        implicit none
00695        integer nx,ny,nz
00696        parameter (nx=200,ny=1,nz=1100)
00697        real  interpol_one,x_old(2481),y_old(2481),xin
00698      end function interpol_one
00699   end interface
00700 
00701   !  write(*,*) x_old
00702   do i=1,nz
00703      interpol(i)=interpol_one(y_old,x_old,x_new(i))
00704   enddo
00705   return
00706 end function interpol
00707 
00708 
00709 function interpol_one(y_old,x_old,x_new)
00710   implicit none
00711   integer nx,ny,nz
00712   parameter (nx=200,ny=1,nz=1100)
00713   real interpol_one,x_old(2481),y_old(2481),x_new
00714   integer loc1,loc2
00715   real temp_x_old(2481)
00716   real x1,x2,fac
00717 !linear
00718   temp_x_old=x_old
00719   loc1=minloc(abs(x_new-temp_x_old),dim=1)
00720   temp_x_old(loc1)=-9e19 
00721   loc2=minloc(abs(x_new-temp_x_old),dim=1)
00722   fac=(x_new-x_old(loc2))/(x_old(loc1)-x_old(loc2))
00723   interpol_one=y_old(loc1)*fac+y_old(loc2)*(1-fac)
00724 
00725 !  write(*,*) x_new,y_old(loc1),y_old(loc2),interpol_one,fac
00726 
00727 
00728 !  write(*,*) interpol_one,exp(interpol_one)
00729 ! cubic spline interpolation
00730 !  
00731   return
00732 end function interpol_one
00733 
00734 
00735