! =================================================================== ! block data init_swind ! =================================================================== ! ! Initialization of common blocks for solar wind - joule heating sub. ! ------------------------------------------------------------------- ! implicit none ! ----- Constants ------------------------------------------------------ ! ---------------------------------------------------------------------- ! model parameters ! ---------------------------------------------------------------------- integer n_lon ! number of longitude beams integer n_lat ! number of latitude beams integer LM ! number of layers parameter (n_lon=96, n_lat=50) PARAMETER (LM=39) ! ----- Commons --------------------------------------------------! !-----------------------------------------------------------------! ! EE : intensity of electric field (n_lon/n_lat) : V/m ! ni : concent. of atm. ions (n_lon/n_lat/LM) : mo/m^3 ! qi : ionization rate (n_lon/n_lat/LM) : 1/m^3/s ! nuu : frequency of e-neut. collision (n_lon/n_lat/LM) : 1/sec ! QJOUL : Joule heating (n_lon/n_lat/LM) : K/sec real*4 EE (n_lon,n_lat) real*4 ni (n_lon,n_lat,LM) real*4 qi (n_lon,n_lat,LM) real*4 nuu (n_lon,n_lat,LM) real*4 Qjoul (n_lon,n_lat,LM) common / qqJOUL / EE, ni, qi, nuu, Qjoul DATA EE /nEE*0./ DATA qi /nnX*0./ DATA ni /nnX*0./ DATA nuu /nnX*0./ DATA Qjoul /nnX*0./ end ! ======================================================================== ! SUBROUTINE SW_JOULq(nnYEAR, nday, tofday) ! ======================================================================== ! ! Calculation of the stratospheric Joule heating rate induced by the Solar ! wind input. !--------------------------------------------------------------------------! ! INPUT: ! ====== ! ------------------- arguments of subroutine ! nnYEAR - number of year to choose solar wind parameters ! from 1979 to 2002 ! ! nday - number of the year day (1-365) ! tofday - time of day in hours (1 - 24) ! ! ------------------- Commons ! T (K) - 3D temperature field (n_lont,n_lat,LM) ! prest (hPa) - 3D pressure filed (n_lont,n_lat,LM) ! YM (grad) - latitude array (n_lat) ! XM (grad) - longitudinal array (n_lon) ! ! ! Bz(nT) - vertical component of IMF :: (SWyr, SWda, SWm) ! Vsw(km/s) - velocity of solar wind :: (SWyr, SWda, SWm) ! nsw(mo/m^3) - ion concentration of Swind:: (SWyr, SWda, SWm) ! ------------------------------------------------------------------------ ! ! These are 3 fields over period 1979-2002 are input data for sub. ! e.g. Bz( nnYEAR, nday, tofday) ! 1979:1999 1:365 1-4 ! 00:00 06:00 12:00 18:00 GMT ! ! OUTPUT ! ====== ! ! ------------------- Commons ! ! EE : intensity of electric field (n_lon/n_lat) : V/m ! ni : concent. of atm. ions (n_lon/n_lat/LM) : mo/m^3 ! qi : ionization rate (n_lon/n_lat/LM) : 1/m^3/s ! nuu : frequency of e-neut. collision (n_lon/n_lat/LM) : 1/sec ! QJOUL : Joule heating (n_lon/n_lat/LM) : K/sec !-----------------------------------------------------------------! ! ! ======================================================================== ! implicit none !--------------------------------------------------------------7-2 ! Input data !--------------------------------------------------------------7-2 ! ----- Arguments --------------------------------------------- ! integer nnYEAR, nday real tofday ! ----- Commons ----------------------------------------------- ! ! ---------------------------------------------------------------------- ! model parameters ! ---------------------------------------------------------------------- integer n_lon ! number of longitude beams integer n_lat ! number of latitude beams integer LM ! number of layers integer LM1 ! number of boundares parameter (n_lon=96, n_lat=50) PARAMETER (LM=39, LM1=LM+1) ! ----- Commons ------------------------------------------------------ ! -------------------------------------------------------------------- ! ! Solar wind parameters and input common blocks ! ! -------------------------------------------------------------------- ! integer SWm, SWda, SWyr parameter ( SWm = 4, SWda = 365, SWyr = 21+3) !-----------------------------------------------------------------! ! Bz : vertical component of IMF : nT ! Vsw : velocity of solar wind : km/s ! nsw : ion concentration of Swind : mo/m^3 real*4 Bz (SWyr,SWda,SWm) real*4 Vsw (SWyr,SWda,SWm) real*4 nsw (SWyr,SWda,SWm) common / SWind / Bz, Vsw, nsw ! T (K) - 3D temperature field (n_lont,n_lat,LM) ! prest (hPa) - 3D pressure filed (n_lont,n_lat,LM) ! YM (grad) - latitude array (n_lat) ! XM (grad) - longitudinal array (n_lon) real, dimension (n_lon,n_lat,LM) :: T, PREST real, dimension (n_lon) :: XM real, dimension (n_lat) :: YM common /xINPUTS/ T, PREST, YM !--------------------------------------------------------------7-2 ! Output data !--------------------------------------------------------------7-2 !-----------------------------------------------------------------! ! EE : intensity of electric field (n_lon/n_lat) : V/m ! ni : concent. of atm. ions (n_lon/n_lat/LM) : mo/m^3 ! qi : ionization rate (n_lon/n_lat/LM) : 1/m^3/s ! nuu : frequency of e-neut. collision (n_lon/n_lat/LM) : 1/sec ! QJOUL : Joule heating (n_lon/n_lat/LM) : K/sec real*4 EE (n_lon,n_lat) real*4 ni (n_lon,n_lat,LM) real*4 qi (n_lon,n_lat,LM) real*4 nuu (n_lon,n_lat,LM) real*4 Qjoul (n_lon,n_lat,LM) common / qqJOUL / EE, ni, qi, nuu, Qjoul ! ----- Local variables --------------------------------------- ! ! ----- Constants ------------------------------------------------------------- integer n11min, n11max parameter (n11min = 10, n11max = 14 ) integer I, J, K, ila integer t1, t2, iyear1, iyear2, inday1, inday2 x, isolar, iyear , iprest real xt , plevel real r0 , Teta, alfa, Dp real e , me , mp , Mi real na , ro , Tem , cp x, X1 , X2 , xBz , xVsw, xnsw, xRe x, xEE , xBB , xni , xqi, xnuu, xQjoul x, xnaSW, naSW(n_lon,n_lat,LM), tofday1 x, lambda x, localT, argum, PI, mm, Aee, Bee x, Acc, Bcc, xEEe, xEEc, localTT real AA, BBmax, BBmin real Aeff real y11min(n11min), y11max(n11max) logical aanii, ptest, iclock ! ----- parameters of routine --------------------- ! parameter ( x mp = 1.673e-27 ! proton mass (kg) x, me = 9.110e-31 ! electron mass (kg) x, e = 1.602e-19 ! electron charge (q) x, Mi = 400. ! mean ion mass number (Mion = Mi * mp) x, cp = 1.005e3 ! specific heat capacity at constant pressure (J*kg/K) x, AA = 1.740e-18 ! constant for latitude function of ionization rate x, BBmax = 2.840e-17 ! -"- for minimum of 11yr-solar activity x, BBmin = 1.930e-17 ! -"- for maximum of 11yr-solar activity x ) data PI/3.14159592/ data y11min /1984,1985,1986,1987,1988,1994,1995,1996,1997,1998/ data y11max /1979,1980,1981,1982,1983,1989,1990,1991,1992,1993 > ,1999,2000,2001,2002/ aanii = .TRUE. ptest = .FALSE. ! .TRUE. iclock = .TRUE. data iprest/1/ c ------------------ What is nnYEAR ---> 11-max or 11-min year. ---- c c--------------------------------------------------------------------c if( nnYEAR .eq. 2000) then write(*,*) 'Sorry, the year 2000 data are not available' write(*,*) 'Please, use the 1999 or 2001 years' stop endif isolar = 2 do iyear = 1,n11max if (nnYEAR .eq. y11max(iyear)) isolar = 1 enddo do iyear = 1,n11min if (nnYEAR .eq. y11min(iyear)) isolar = 0 enddo if( isolar .eq. 2) then write(*,*) 'Sorry, you are wrong with nnYEAR parameter' stop endif c ----- Bz, Vsw and nsw for time moment - nnYEAR, nday, tofday ----- c c--------------------------------------------------------------------c c-------- time moment if (tofday .eq. 24.0) then tofday1 = 23.999 else tofday1 = tofday endif t1 = int(tofday1) / 6 + 1 t2 = t1 + 1 xt = tofday1 - float(t1-1) * 6. if (t2 .eq. 5 .and. nday .ne. 365) then iyear1 = nnYEAR iyear2 = nnYEAR inday1 = nday inday2 = nday+1 t2 = 1 else iyear1 = nnYEAR iyear2 = nnYEAR inday1 = nday inday2 = nday t2 = t1 + 1 endif if (t2 .eq. 5 .and. nday .eq. 365) then if(nnYEAR .eq. 2002) then iyear1 = nnYEAR iyear2 = nnYEAR inday1 = nday inday2 = nday t2 = t1 else iyear1 = nnYEAR iyear2 = nnYEAR+1 inday1 = nday inday2 = 1 t2 = 1 endif endif iyear1 = iyear1 - 1979 + 1 iyear2 = iyear2 - 1979 + 1 c-------- Bz X1 = Bz(iyear1,inday1,t1) X2 = Bz(iyear2,inday2,t2) xBz = X1 + ( X2 - X1 ) / 6. * xt c-------- Vsw X1 = Vsw(iyear1,inday1,t1) X2 = Vsw(iyear2,inday2,t2) xVsw = X1 + ( X2 - X1 ) / 6. * xt c-------- nsw X1 = nsw(iyear1,inday1,t1) X2 = nsw(iyear2,inday2,t2) xnsw = X1 + ( X2 - X1 ) / 6. * xt if(ptest) then ! ----------------------- test printing----start write(2,*) ' ' write(2,*) ' ############################' write(2,*) ' year1 :', iyear1 write(2,*) ' year2 :', iyear2 write(2,*) ' nday1 :', inday1 write(2,*) ' nday2 :', inday2 write(2,*) ' time1 :', t1 write(2,*) ' time2 :', t2 write(2,*) ' Universal time (hour)=tofday1 :', tofday1 write(2,*) ' xt :', xt write(2,*) ' Bz1 = ', Bz(iyear1,inday1,t1) write(2,*) ' Bz2 = ', Bz(iyear2,inday2,t2) write(2,*) ' Bz = ', xBz write(2,*) ' Vsw1 = ', Vsw(iyear1,inday1,t1) write(2,*) ' Vsw2 = ', Vsw(iyear2,inday2,t2) write(2,*) ' Vsw = ', xVsw write(2,*) ' nsw1 = ', nsw(iyear1,inday1,t1) write(2,*) ' nsw2 = ', nsw(iyear2,inday2,t2) write(2,*) ' nsw = ', xnsw write(2,*) '*************************' write(2,*) ' ' write(2,*) ' ' write(2,*) ' ' ****** pause !---------------------------- test printing end endif *----------------------------- pressure of solar wind Dp= 1.7e-6 * xNsw * xVsw * xVsw * --------------------------- checking 1 if( Dp .eq. 0.) then write(500,*) ' Dp = 0', Dp , tofday1, nday write(500,*) ' Nsw ',xNsw write(500,*) ' Vsw ',xVsw write(*,*) ' Dp = 0', Dp , tofday1, nday write(*,*) ' Nsw ',xNsw write(*,*) ' Vsw ',xVsw Dp = 1. endif if( Dp .lt. 0.) then write(500,*) ' Dp < 0', Dp , tofday1, nday write(500,*) ' Nsw ',xNsw write(500,*) ' Vsw ',xVsw write(*,*) ' Dp < 0', Dp , tofday1, nday write(*,*) ' Nsw ',xNsw write(*,*) ' Vsw ',xVsw Dp = 1. endif if( Dp .gt. 100.) then write(500,*) ' Dp very large', Dp , tofday1, nday write(500,*) ' Nsw ',xNsw write(500,*) ' Vsw ',xVsw write(*,*) ' Dp very large', Dp , tofday1, nday write(*,*) ' Nsw ',xNsw write(*,*) ' Vsw ',xVsw Dp = 1. endif write(*,*) 'check1 is ok', Dp * --------------------------- check 1 do I = 1, n_lon do J = 1, n_lat c --- calculation of the magnetopause radius Re ----------- c c----------------------------------------------------------c if(xBz .ge. 0.) then r0 = (11.4 + 0.013 * xBz) * Dp**(-1./6.6) else r0 = (11.4 + 0.140 * xBz) * Dp**(-1./6.6) endif alfa = (0.58 - 0.010 * xBz)* > (1. + 0.010 * Dp) * ------------------ Longitudinal dependency of EE * ------------------------------------------------- * sun localT = tofday + (XM(I)/15.) localTT = localT if (localT .gt. 24.) localT = localT - 24. localT = localT - 12. ! localT = 0. <=> noon if(iclock) then *-------------- clock angle of longitudinal, mm = PI ! rate of cos decreasing at argum from 0 to +-1 argum = localT / 12. ! argum = [-1,1] if( localT .ge. -6. and. localT .lt. 6) then !################################################################# c ------------ calculation of the E(V/m) ----------------- c c----------------------------------------------------------c Teta = cos(mm*argum) xRe = r0 * (2./(1.+ Teta)) ** alfa ! distance to magnetopause ! in radius of Earth xEE = (527.01 - 31.5 * xRe) * 1.e-3 ! V/m if(xEE .lt. 0.) xEE = 0. ! physical correction !################################################################# else !################################################################# Teta = 0. ! Teta = cos(pi/2) xRe = r0 * (2./(1.+ Teta)) ** alfa ! distance to magnetopause xEEc = (527.01 - 31.5 * xRe) * 1.e-3 ! V/m if(xEEc .lt. 0.) xEEc = 0. ! physical correction xEE = xEEc * 0.4 !################################################################# endif EE(I,J) = xEE *----------------- end of iclock endif if(ptest) then if(I .eq. 5 .and. J .eq. 47) then write(2,*) '----------latitude ', Ym(J) write(2,*) '----------longitude', Xm(I) write(2,*) ' ' write(2,*) ' Universal time (hour) :', tofday write(2,*) ' Local time (hour) :', localTT write(2,*) ' Local time - 12 (hour) :', localT write(2,*) ' argum = (Local time - 12)/12 :', argum write(2,*) ' mm = pi :', mm write(2,*) ' ' write(2,*) ' ' write(2,*) ' ' write(2,*) '----clock angle=argum*mm ', argum*mm write(2,*) '----cos(clock angle) ', Teta write(2,*) '----nsw ', xnsw write(2,*) '----Vsw ', xVsw write(2,*) '----Dp ', Dp write(2,*) '----alfa ', alfa write(2,*) '----r0 ', r0 write(2,*) '--(2/1+cos(clock angle))**alfa ', (2./(1.+Teta))**alfa write(2,*) ' ' write(2,*) ' ' write(2,*) '--------- A ', AA write(2,*) '--------- B ', BBmax write(2,*) '------ sin ', sin(Ym(J)*.017453) write(2,*) '----sin**4 ', sin(Ym(J)*.017453)**4 write(2,*) ' ' write(2,*) ' ' write(2,*) '- mi ', mi write(2,*) '- mp ', mp write(2,*) '- me ', me write(2,*) '- me*me ', me*me write(2,*) '- mi*mp/me/me', mi*mp/me/me write(2,*) '- e ', e write(2,*) '-(eE)**2 ', e*e*xEE*xEE write(2,*) '- cp ', cp write(2,*) '------------------------------------------' write(2,*) ' Bz(nT) Re(radE) EE(V/m) nuu(1/s) ni(1/m^3) qi(1/m^3/s) HQ(K/day) Tem(K) ro(kg/m^3) na(1/m^3) Aeff plev(hPa) ' endif endif do K = 1, LM plevel = prest(I,J,K) if ( plevel .le. 100. .and. plevel .ge.7) then *######################################################## from 100 to 7 hPa c -------- Temperature (K) Tem = T(I,J,K) c -------- density of neutral atmosphere (kg/m^3) ro = prest(I,J,K) * 100. / 287. / Tem c -------- concentration of neutral particles (mol/m^3) xnaSW = prest(I,J,K) * 100./287./Tem / 28.96 * 6.02e26 naSW(I,J,K) = xnaSW c -------- ionization rate if ( isolar .eq. 1) then xBB = BBmin elseif( isolar .eq. 0) then xBB = BBmax endif xqi = (AA + xBB * sin(YM(J)*0.017453)**4 ) * xnaSW * 5. ! qi(I,J,K) = xqi !-------------- Aeff dependency on altitude (Makarova & Shirochkov) if(plevel .ge. 50.) Aeff = 1.000e-09 if(plevel .le. 50. .and. plevel .ge. 10.) then Aeff = (5.e-9 - 1.e-9)/(10.-50.)*(plevel-50.)+1.e-9 endif if(plevel .le. 10. .and. plevel .ge. 3.) then Aeff = (1.e-8 - 5.e-9)/(3.- 10.)*(plevel-10.)+5.e-9 endif if(plevel .le. 3.) then Aeff = 1.e-8 endif c -------- ion concentration (ion/m^3) xni = sqrt ( xqi / Aeff ) ni(I,J,K) = xni c -------- frequency of e-neutral collisions (1/sec) xnuu = 4.4e-15 * (Tem/300.) * xnaSW nuu(I,J,K) = xnuu if(ptest) then ! ----------------------- test printing----start write(3,3) xnuu, Tem, Tem/300., xnaSW, xqi, Aeff, plevel !---------------------------- test printing end endif c--------- calculation of Joule Heating rate(K/day) if(aanii) then xQJOUL = xni * e*e * xEE**2 * (mi*mp/me/me/xnuu) ! W/m^3 xQJOUL = xQJOUL / cp / ro * 8.64e4 ! K/day QJOUL(I,J,K) = xQJOUL ! K/day endif if(ptest) then if(I .eq. 5 .and. J .eq. 47 .and. plevel .ge. 3 .and. plevel .le. 70) then ! ----------------------- test printing----start write(2,1) xBz, xRe, xEE, xnuu, xni, xqi, xQjoul,Tem,ro,xnaSW,Aeff,plevel !---------------------------- test printing end endif endif *######################################################### from 100 to 7 hPa endif enddo ! K enddo ! J enddo ! I 1 format(13(1pe11.3)) 3 format(7(1pe12.4)) return end