PROGRAM Setup_plasma USE rcm_mod_subs ! ! Code to set up plasma parameters and initial distribution of the ! inner edges and grid-based plasma for RCM. ! Must run RCM grid setup and RCM B-field setup first to use this code. ! IMPLICIT NONE INTEGER (iprec) :: type_of_plasmasheet, ie, i, & Lrec, j, k, n_pt, nn, kc, option_eeta, option_edges,& i_deepest, plasma_type, i_last, ii, n_t, istat, t,tsize, ikk, iss REAL (rprec), DIMENSION (iesize, jsize) :: N, N_c, N_h, N_3, kT, kT_c, kT_h, kT_3, & E0, E0_c, E0_h, E0_3, E_avg, E_avg_c, E_avg_h, E_avg_3, & U, P, eta_an, pvg_an, & eta_numeric, eta_numeric_grid, & pvg_numeric, pvg_numeric_grid, & kT_numeric, N_numeric, p_numeric, U_numeric, E_avg_numeric, etac_ REAL (rprec), DIMENSION (isize, jsize, ksize) :: eta_init REAL (rprec), DIMENSION (iesize) :: eta_numeric_edge, pvg_numeric_edge REAL (rprec) :: y (iesize), x_i(iesize),kappa (iesize),kappa_c(iesize,jsize),& kappa_h(iesize,jsize), kappa_3(iesize,jsize),kappaave(jsize,391,10,2),& kappa_c_e(391,jsize), kappa_h_e(391,jsize), kappa_c_i(391,jsize), & kappa_h_i(391,jsize), & kappa_c_o(391,jsize), kappa_h_o(391,jsize), kappa_3_o(391,jsize), k_c (iesize), & T_ps(jsize),T_pse_c(391,jsize), T_pse_c_ind(jsize), T_pse_h_ind(jsize), & T_pse_h(391,jsize),T_psi_c(391,jsize), T_psi_h(391,jsize),T_psi_c_ind(jsize),& T_psi_h_ind(jsize), & N_ps(jsize),N_pse_c(391,jsize),N_pse_c_ind(jsize),N_pse_h_ind(jsize),& N_pse_h(391,jsize), N_psi_c(391,jsize),N_psi_h(391,jsize),N_psi_c_ind(jsize),& N_psi_h_ind (jsize), & n_sw, V_sw, delta_E1, & lambda_max (iesize), lambda_min (iesize), lambda_c (iesize),& eta_tmp, p_tmp, & Solve_for_y, Fdist_kappa_dflux, Fdist_3_kappa_dflux, & p_chan(ksize), pc_chan(kcsize,jsize), & vm_20, bi_20, bj_20,vm_bnd(jsize), bi_bnd(jsize), bj_bnd(jsize), & x_pt, y_pt, x_1, y_1, x_2, y_2, tmp,& bi_initial_edges (jsize, ksize), & bi_initial_grid (jsize,kcsize), Kp_for_composition, r_pp, fudgeAE(391) REAL (rprec), ALLOCATABLE :: n_sw_inp(:), v_sw_inp(:), imf_bz_inp(:) REAL (rprec), ALLOCATABLE :: delta_E(:), energy (:) CHARACTER (LEN=1) :: answer CHARACTER (LEN=6), DIMENSION(ksize) :: remark_edge='-------' CHARACTER (LEN=6), DIMENSION(kcsize) :: remark_grid='-------' CHARACTER (LEN=12) :: char_date, char_time CHARACTER (LEN=50) :: char_label, title_string=' ' CHARACTER (LEN=80) :: form_string CHARACTER (LEN=80) :: filename_sw='inputs.dat' ! file with solar wind inputs LOGICAL :: FD, L_bc_same_as_ic ! ! WRITE (*,'(//A)') 'SETUP OF PLASMA INITIAL AND BOUNDARY CONDITIONS FOR THE RCM' WRITE (*,'(//,A)') '---------------------------------------------------------' ! IF (ksize /= kcsize) STOP 'ksize/=kcsize' ! ! CALL Rcm_read_params () ibnd_type=1 ! CALL Read_kp () CALL Read_grid () CALL Read_bfield () ! read mark-times for B-field DO i = 1, SIZE(ibtime) CALL Get_bfield (ibtime, ibtime(i), .TRUE.) ! retrieve B-field for T write (*,'(A,F6.2)') ' GET_BFIELD: Fstoff=', fstoff CALL Get_boundary (boundary, bndloc) CALL Get_bfield_at_minus20_0_0 (bi_20,bj_20, vm_20) !!! Carefull!! at 6.6 Re!!! WRITE (*,'(A,I7.7, TR5,A,F5.2)') ' TIME=',ibtime(i), ' VM(-6.6,0,0)=',vm_20 END DO WRITE (*,'(/A/)') 'IF BOUNDARY CONDITION WILL BE TIME-INDEPENDENT,' WRITE (*,'(A)') 'CHOOSE VALUE OF VM(-6.6) TO USE WITH B-FIELD AT T=0' WRITE (*,'(/,T2,A)', ADVANCE='NO') 'ENTER VM (-6.6,0,0): ' READ (*,*) vm_20 WRITE (*,'(//A)') 'IF BOUNDARY CONDITION WILL BE TIME-INDEPDENT,' WRITE (*,'(A)') 'CHOOSE VALUE OF KP FOR THE ION COMPOSITION,' WRITE (*,'(A)') '[OTHERWISE, IT WILL VARY BASED ON KP IN RCMKP_INP FILE,' WRITE (*,'(A)') ' TO ALLOW THAT, ENTER -1.0 FOR KP]' WRITE (*,'(/A)',ADVANCE = 'NO') 'ENTER KP FOR ION COMPOSITION:_' READ (*,*) Kp_for_composition ! ! OPEN (LUN_1, FILE=FILENAME_SW, STATUS='OLD') READ (LUN_1,*) READ (LUN_1,*) n_t = 0 do_loop_nt: DO READ (LUN_1,'(A80)',END=101) form_string n_t = n_t + 1 END DO do_loop_nt 101 REWIND (LUN_1) print*,'kcsize,n_t',kcsize,n_t ALLOCATE (itime_etac(n_t), etac_inp (kcsize,n_t,j), n_sw_inp(n_t), v_sw_inp(n_t),imf_bz_inp(n_t)) ! READ (LUN_1,*) READ (LUN_1,*) DO ii = 1, n_t READ (LUN_1,*,IOSTAT=istat) itime_etac(ii), & n_sw_inp(ii), v_sw_inp(ii),tmp,tmp,tmp,tmp,tmp,imf_bz_inp(ii) END DO CLOSE (LUN_1) form_string = '' ! ! 2. Set up plasma parameters: ! ! ! Option 1: ! Average plasma sheet temperature is from Wolf writeup based on Huang ! and Frank 1986 GRL paper. Later, Huang and Frank published a correction ! stating that their Temperatures were really energies (1.5*T). We correct ! for this here by dividing the Wolf expressions by 1.5. ! ! Option 2: ! Plasma sheet n and T are from Borovsky et al 1998 (jgr), based on the solar ! wind n and V values delayed by about 2 hours. ! ! Option 3: Tsyganenko_Mugai 2003 statistical values ! ! Option 4: Choose N_ps and T_ps manually. ! WRITE (*,'(///,A,//)') 'THIS IS SET UP OF ALAM VALUES FOR THE RCM.' WRITE (*,'(A)') 'IT NEEDS SOME VALUES FOR ION PLASMA SHEET N AND T AT -20 RE' WRITE (*,'(A)') 'FOR THE PURPOSES OF CHOOSING ALAM VALUES' WRITE (*,'(A)') '(DURING RUNS, N AND T CAN VARY REGARDLESS & &OF THE VALUES USED HERE)' WRITE (*,'(A)') ' THERE ARE 4 WAYS TO CHOOSE N AND T:' WRITE (*,'(A)') ' OPTION 1: USE Kp-BASED STATISTICAL FORMULA [MSM]' WRITE (*,'(A)') ' OPTION 2: USE BOROVSKY ET AL SOLAR WIND-BASED FORMULA' WRITE (*,'(A)') ' OPTION 3: Tsyganenko_mukai 2003' WRITE (*,'(A)') ' OPTION 4: SPECIFY YOUR OWN VALUES' WRITE (*,'(A)') ' IN ANY CASE, VALUES USED SHOULD PROBABLY BE STORM-TIME VALUES' WRITE (*,'(A)') ' IF YOU PLAN TO USE TIME-DEPENDENT PLASM B.C., ONLY OPTION 3 works' WRITE (*,'(A)',ADVANCE='NO') ' ENTER OPTION DESIRED: ' READ (*,*) type_of_plasmasheet IF (type_of_plasmasheet < 1 .OR. type_of_plasmasheet > 4) STOP 'ILLEGAL OPTION' ! IF (type_of_plasmasheet == 2) THEN WRITE (*,'(A)',ADVANCE='NO') & 'ENTER 2-HOUR DELAYED SOLAR WIND N [CM-3] AND V [KM/S]: ' READ (*,*) n_sw, v_sw N_ps = 0.0785*n_sw**0.62 T_ps = (-3.65+0.0190*V_sw)*1000.0 ! rescale from -20 Re (Borovsky) to -13 Re based on Paterson et al: N_ps = N_ps *1.0 T_ps = T_ps *1.27 ELSE IF (type_of_plasmasheet == 1) THEN WRITE (*,'(T2,A,F4.2,A)',ADVANCE='NO') 'ENTER Kp [DEFAULT IS 1.74]: ' READ (*,*) kp !defines the electron number density which is N_ps=N(1) N_ps = 0.25*(0.4*(kp-1.0)+0.5*(5.0-kp)) ! in cm-3 T_ps = ((1.0E+3)*(0.2885*(kp-1)+0.15625*(5-kp)))*7.8/1.5/0.5/1.5 ! in eV ELSE IF (type_of_plasmasheet == 3) THEN WRITE (*,'(A)',ADVANCE='NO') & 'ENTER SOLAR WIND N [CM-3]. V [KM/S], IMF_Bz [nT]: ' READ (*,*) n_sw, v_sw, imf_bz CALL Central_plasma_sheet_tm2003 (-13.0, 0.0, n_sw, v_sw, & imf_bz, t_ps, n_ps) ELSE IF (type_of_plasmasheet == 4) THEN ! OPEN (41, FILE='inputs_T_n_kappa.dat') tsize=289 !!!!!!!!!!!!!! open Chih-ping files !!!!!!!!!! OPEN (41, FILE='kappa-8hr-IMF_for-RCM-setup-plasma.txt') ! generate in /data/arc1/RCM-MURI/idlcode/create-kappa-N-T-SVM_OMNI_for-RCM-setup-plasma_20130317-18.pro READ (41,*)((((kappaave(j,t,ikk,iss),iss=1,2),ikk=1,10),t=1,tsize),j=1,jsize) ! j:MLT, t:time, ikk:kappa para, iss:i+ or e- ! print*,'kappa=',kappaave(26,30,2,1),kappaave(30,130,10,2) ! print*,'kappa1=',kappaave(27,100,1:10,2) print*,'kappa1=',kappaave(j,1,1,1),kappaave(j,1,4,1) CLOSE (41) print*,'after openining 8hr average boundary condition.txt' OPEN (43, FILE='fudge-AE_for-RCM-setup-plasma.txt') ! generate in /RCM_theory/geotail/create-fudge-AE.pro READ (43,*)(fudgeAE(t),t=1,tsize) ! t:time CLOSE (43) !!!!! If we want initial conditions!!!!!! OPEN (51, FILE='eta_initial.dat') READ (51,*) eta_init CLOSE (51) print*,'after openining eta_initial.dat' ! print*, eta_init(:,27,39) ELSE STOP 'ILLEGAV VALUE OF TYPE_OF_PLASMASHEET' END IF ! ! for the initial conditions, read in the initial magnetic field: CALL Get_bfield (ibtime, ibtime(1), .TRUE.) ! retrieve B-field for T CALL Get_boundary (boundary, bndloc) WRITE (*,'(A)',ADVANCE='NO') 'USE TIME-INDEPENDENT B.C.(==INIT. COND)? [T/F]: ' READ (*,*) L_bc_same_as_ic IF (L_BC_SAME_as_ic) THEN label%intg(6)=0 FD = .TRUE. CALL Write_array ('rcmetac_inp', 1, Label, ARRAY_2D = etac, & SETUP = FD, ASCI=asci_flag) ELSE WRITE (*,*) 'TIME-DEPENDENT PLASMA BOUNDARY CONDITION:' DO ii = 1, n_t ! begin time loop ! DO ii = 1,2 DO ie = 1,iesize ! print*,'time loop=',ii !!!!!! Set up energy channels:!!!!!! ! ! IF (kmin (ie) >= kmax(ie)) THEN WRITE (*,*) 'species IE=', ie, 'NEEDS AT LEAST 2 CHANNELS' STOP 'CHANGE AND RE-RUN' END IF ! k_c (ie) = (kmax(ie) + kmin (ie) ) / 2.0 !kmax, kmin defined in rcm_include.h ! lambda_max (ie) = SIGN(1.0, ie-1.5) * 7.5 * kT(ie,27) / vm_bnd(27) ! lambda_min (ie) = SIGN(1.0, ie-1.5) * 0.01 * (1.5*kT(ie,27) / vm_bnd(27)) ! ! lambda_max (ie) = SIGN(1.0, ie-1.5) * 7.5 * kT(ie,27) /vm_20 ! lambda_min (ie) = SIGN(1.0, ie-1.5) * 0.1 * (1.5*kT(ie,27) /vm_20) lambda_max (ie) = SIGN(1.0, ie-1.5) * 1000000. / vm_20 !sign(..)=-1 for e-,+1 for protons lambda_min (ie) = SIGN(1.0, ie-1.5) * 10. / vm_20 ! print*, 'ie=',ie, lambda_max(ie), lambda_min(ie) y (ie) = Solve_for_y (lambda_max(ie)/lambda_min(ie)) x_i (ie) = 2.0*y(ie) / (kmax(ie) - kmin (ie)) lambda_c(ie) = lambda_min(ie) * EXP(x_i(ie)*((kmax(ie)-kmin(ie))/2.0-1.0)) ! DO k = kmin (ie), kmax (ie) IF (k < k_c(ie)) THEN alam (k) = EXP (x_i(ie)*(k-k_c(ie))) ELSE alam (k) = 1.0 + x_i(ie) * (k-k_c(ie)) END IF alam (k) = alam (k) * lambda_c(ie) alamc(k) = alam (k) ! if(ie.eq.2.and.k.le.120)print*, 'alam(k)',k,alam(k) END DO ! ! DO k = kmin (ie), kmax (ie) ikflav (k) = ie ikflavc(k) = ie IF (ie == 1) THEN ! fudge (k) = 1.0/3.0 fudge (k) = fudgeAE(ii) IF (k == 1) fudge(k) = 0.0 !plasmapause !!! ELSE fudge (k) = 0.0 END IF fudgec (k) = fudge (k) WRITE (remark_edge(k),'(F5.3)') fudge(k) WRITE (remark_grid(k),'(F5.3)') fudgec(k) END DO ! WRITE (*,'(A, A2, A, F7.1, A, F17.1, A, I4.4)') & 'SPECIES ', species_char(ie), ' ALAM_MIN=', lambda_min(ie), & ' ALAM_MAX=', lambda_max(ie), ' NUM. OF CHANNELS=', kmax(ie)-kmin(ie)+1 ! ALLOCATE (energy(kmin(ie):kmax(ie)), delta_e (kmin(ie):kmax(ie))) energy = ABS(alam(kmin(ie):kmax(ie))) * vm_20 !print*,'energy=', energy CALL Get_delta_E (kmin(ie), kmax(ie), energy, delta_E ) ! ! Modify the width of the highest energy channel: !! do j=1,jsize if (ie.eq.1) then kappa_h(ie,j)=kappaave(j,1,1,1) kappa_c(ie,j)=kappaave(j,1,4,1) end if if (ie.eq.2) then kappa_h(ie,j)=kappaave(j,1,1,2) kappa_c(ie,j)=kappaave(j,1,4,2) end if if (ie.eq.3) then kappa_h(ie,j)=6. kappa_c(ie,j)=6. end if end do kappa_h(1,27)=6. !CHANGE I MADE kappa_h(2,27)=6. !CHANGE I MADE delta_E(kmax(ie)) = energy(kmax(ie))**(kappa_h(ie,27)-0.5) / & (kappa_h(ie,27)-1.5) / & (0.5*(energy(kmax(ie))+energy(kmax(ie)-1)))**(kappa_h(ie,27)-1.5) if(ii.eq.1)then print*,'term1=', energy(kmax(ie))**(kappa_h(ie,27)-0.5) print*,'term2=', (kappa_h(ie,27)-1.5) print*,'term3=',(0.5*(energy(kmax(ie))/energy(kmax(ie)-1)))**(kappa_h(ie,27)-1.5) print*,'Emax=',delta_E(kmax(ie)),kappa_h(ie,27),ie,energy(kmax(ie)) endif ! !!!!!!!! End set up energy channels !!!!!!! ! n_sw = n_sw_inp (ii) ! v_sw = v_sw_inp (ii) ! imf_bz = imf_bz_inp(ii) ! CALL Central_plasma_sheet_tm2003 (-13.0, 0.0, n_sw, v_sw, & ! imf_bz, t_ps, n_ps) CALL Get_bfield (ibtime, ibtime(ii),.true.) ! retrieve B-field for T IF (ibtime(ii) /= itime_etac(ii) ) STOP 'ibtime and itime_etac differ!' ! although don't have to do this, jtime can be itime_etac, and it will interpolate! ! CALL Get_boundary (boundary, bndloc) ! CALL Get_bfield_at_minus13_0_0 (bi_13,bj_13, vm_13) ! Kp_for_composition = Get_kp (ikptime, kpinput, jtime=itime_etac(ii)) ! DO j=1,jsize ! vm_bnd(j) = Gntrp (vm, bndloc(j), REAL(j,rprec), 0) kappa_h(1,j)=kappaave(j,ii,1,2) kappa_c(1,j)=kappaave(j,ii,4,2) kappa_3(1,j)=kappaave(j,ii,8,2) kappa_h(2,j)=kappaave(j,ii,1,1) kappa_c(2,j)=kappaave(j,ii,4,1) kappa_3(2,j)=kappaave(j,ii,8,1) kappa_h(3,j)=6. kappa_c(3,j)=6. kappa_3(3,j)=6. N_c (1,j) = kappaave(j,ii,5,2) N_h (1,j) = kappaave(j,ii,2,2) N_3 (1,j) = kappaave(j,ii,9,2) N_c (2,j) = kappaave(j,ii,5,1) N_h (2,j) = kappaave(j,ii,2,1) N_3 (2,j) = kappaave(j,ii,9,1) N_c (3,j) = 0. N_h (3,j) = 0. N_3 (3,j) = 0. kT_c (1,j) = kappaave(j,ii,6,2)*1000. kT_h (1,j) = kappaave(j,ii,3,2)*1000. kT_3 (1,j) = kappaave(j,ii,10,2)*1000. kT_c (2,j) = kappaave(j,ii,6,1)*1000. kT_h (2,j) = kappaave(j,ii,3,1)*1000. kT_3 (2,j) = kappaave(j,ii,10,1)*1000. kT_c (3,j) = 5000. kT_h (3,j) = 5000. kT_3 (3,j) = 5000. kT(1,j) = (N_c (1,j)*kT_c (1,j) + N_h (1,j)*kT_h (1,j) + N_3 (1,j)*kT_3 (1,j))/(N_c (1,j)+ N_h (1,j) + N_3 (1,j) ) kT(2,j) = (N_c (2,j)*kT_c (2,j) + N_h (2,j)*kT_h (2,j) + N_3 (2,j)*kT_3 (2,j))/(N_c (2,j)+ N_h (2,j) + N_3 (2,j)) kT(3,j) = (N_c (3,j)*kT_c (3,j) + N_h (3,j)*kT_h (3,j) + N_3 (3,j)*kT_3 (3,j))/(N_c (3,j)+ N_h (3,j) + N_3 (3,j)) E_avg_c (1,j) = 1.5 * kT_c (1,j) E_avg_h (1,j) = 1.5 * kT_h (1,j) E_avg_3 (1,j) = 1.5 * kT_3 (1,j) E_avg_c (2,j) = 1.5 * kT_c (2,j) E_avg_h (2,j) = 1.5 * kT_h (2,j) E_avg_3 (2,j) = 1.5 * kT_3 (2,j) E_avg_c (3,j) = 1.5 * kT_c (3,j) E_avg_h (3,j) = 1.5 * kT_h (3,j) E_avg_3 (3,j) = 1.5 * kT_3 (3,j) kappa_c_o(ii,j)=6. kappa_h_o(ii,j)=6. kappa_3_o(ii,j)=6. ! E0_c (1,j) = E_avg_c (1,j) / (1.5 * kappa_c_e(ii,j) / (kappa_c_e(ii,j)-1.5)) ! E0_c (2,j) = E_avg_c (2,j) / (1.5 * kappa_c_i(ii,j) / (kappa_c_i(ii,j)-1.5)) ! E0_c (3,j) = E_avg_c (3,j) / (1.5 * kappa_c_o(ii,j) / (kappa_c_o(ii,j)-1.5)) ! E0_h (1,j) = E_avg_h (1,j) / (1.5 * kappa_h_e(ii,j) / (kappa_h_e(ii,j)-1.5)) ! E0_h (2,j) = E_avg_h (2,j) / (1.5 * kappa_h_i(ii,j) / (kappa_h_i(ii,j)-1.5)) ! E0_h (3,j) = E_avg_h (3,j) / (1.5 * kappa_h_o(ii,j) / (kappa_h_o(ii,j)-1.5)) E0_c (1,j) = E_avg_c (1,j) / (1.5 * kappa_c(1,j) / (kappa_c(1,j)-1.5)) E0_c (2,j) = E_avg_c (2,j) / (1.5 * kappa_c(2,j) / (kappa_c(2,j)-1.5)) E0_c (3,j) = E_avg_c (3,j) / (1.5 * kappa_c(3,j) / (kappa_c(3,j)-1.5)) E0_h (1,j) = E_avg_h (1,j) / (1.5 * kappa_h(1,j) / (kappa_h(1,j)-1.5)) E0_h (2,j) = E_avg_h (2,j) / (1.5 * kappa_h(2,j) / (kappa_h(2,j)-1.5)) E0_h (3,j) = E_avg_h (3,j) / (1.5 * kappa_h(3,j) / (kappa_h(3,j)-1.5)) E0_3 (1,j) = E_avg_3 (1,j) / (1.5 * kappa_3(1,j) / (kappa_3(1,j)-1.5)) E0_3 (2,j) = E_avg_3 (2,j) / (1.5 * kappa_3(2,j) / (kappa_3(2,j)-1.5)) E0_3 (3,j) = E_avg_3 (3,j) / (1.5 * kappa_3(3,j) / (kappa_3(3,j)-1.5)) if(j.eq.27)print*,'kT_C=', kT_c (1,j), kT_h (1,j),E_avg_c (1,j),E_avg_h (1,j),E0_c (1,j),E0_h (1,j) END DO ! J loop DO k = kmin(ie), kmax (ie) DO j = 1, jsize ! print*,'k=',k,'j=',j,'before flux' energy(k)=ABS(alam(k))*vm_bnd(j) delta_E1=(delta_E(k)/vm_20)*vm_bnd(j) IF (Fdist_3_kappa_dflux (N_c(ie,j)*1.0E+6, N_h(ie,j)*1.0E+6, N_3(ie,j)*1.0E+6, & xmass(ie), E0_c(ie,j)*1.6E-19, E0_h(ie,j)*1.6E-19, E0_3(ie,j)*1.6E-19, & energy(k)*1.6E-19, kappa_c(ie,j), kappa_h(ie,j), kappa_3(ie,j)) < 3.0*TINY(1.0)) THEN etac (k,j) = zero ELSE etac(k,j)=Fdist_3_kappa_dflux (N_c(ie,j)*1.0E+6, N_h(ie,j)*1.0E+6, N_3(ie,j)*1.0E+6, & xmass(ie), E0_c(ie,j)*1.6E-19, E0_h(ie,j)*1.6E-19, E0_3(ie,j)*1.6E-19, & energy(k)*1.6E-19, kappa_c(ie,j), kappa_h(ie,j), kappa_3(ie,j))*4.0*pi*& SQRT(xmass(ie) / 2.0 / (energy(k)*1.6E-19)) * & (delta_E1*1.6E-19) ! number density in SI units END IF ! ! Convert to ETA in program units ([ETA]=1/Weber=1/([B][L][L])): ! etac (k,j) = etac (k,j) * (vm_bnd(j)**(-1.5) / 1.0E-9 * Re * 1.0E+3) ! if(k.eq.2)print*,'after flux' if(etac(k,j).ge.0..and.etac(k,j).le.1.e30)then etac (k,j) = etac (k,j) else etac (k,j) = 0. endif ! if(ii.eq.2.and.j.eq.39.and.k.ge.81.and.k.le.160) print*,'Nc',N_c(ie,j),N_h(ie,j),& ! E0_c(ie,j),E0_h(ie,j),kappa_c(ie,j), kappa_h(ie,j) ! if(ii.eq.2.and.j.eq.39.and.k.ge.81.and.k.le.160) print*,'k=',k, '1st etac(k,j)',etac(k,j) END DO ! j loop END DO ! k loop DEALLOCATE (energy, delta_e ) END DO !end ie loop print*,'end ie loop' ! eta = 0.0 ! ! ! ! Set EETA: if (i,j) is occupied, EETA = ETA, else 0. ! ! DO k = 1, ksize bi_initial_edges (:,k) = bndloc (:) + 0.02 END DO option_eeta = 2 IF (option_eeta /= 1 .AND. option_eeta /= 2) STOP 'TROUBLE' ! IF (option_eeta == 1) THEN DO kc = 1, kcsize bi_initial_grid (:,kc) = bndloc (:) END DO ELSE DO j = 1, jsize i_deepest = MAXVAL (bi_initial_edges(j,:)) DO kc = 1, kcsize bi_initial_grid (j,kc) = REAL(i_deepest) END DO END DO END IF ! DO kc = 1, kcsize DO j = 1, jsize DO i = isize, 1, -1 IF (REAL(i) > bi_initial_grid(j,kc)) THEN ! eeta (i,j,kc) = zero ! no initial condition eeta (i,j,kc) = eta_init(i,j,kc) ELSE eeta (i,j,kc) = etac(kc,j) END IF END DO END DO END DO print*,'1' ! ! ! Here is a module to set up a plasmasphere. We re-design ! the electron channel with k=1 to be alam=0. Will need to ! rescale the remaining electron ETA channels to make up ! for loss of eta, because ETA in the plasmasphere channel ! has nothing to do with plasma sheet electrons. ! WILL NEED INITIAL B-FIELD HERE, MAKE SURE IT IS AVAILABLE!!! ! ! print*,etac(2,49) alam (1) = 0.0 alamc(1) = 0.0 tmp = SUM (etac(1:kmax(1),:))/SUM(etac(2:kmax(1),:)) ! print*,'tmp=',tmp etac (2:kmax(1),:) = etac(2:kmax(1),:)*tmp eta (2:kmax(1)) = eta (2:kmax(1))*tmp eta (1) = -1.0E+34 etac(1,:) = 0.0 !since ETAC is used for boundary condition, keep it 0 for PP print*,'2' ! this shoud have been done before (above): fudge(1) = 0.0 fudgec(1) = 0.0 print*,'21',ii ! print*, etac(2,49) ! ! etac_inp (:,ii,:) = etac (:,:) !2017-09-27, this cause trouble in gfortran, command it out since it doesn't seem to be used print*,'22',ii,itime_etac(ii) label%intg(6)=itime_etac(ii) print*,'23' IF (ii == 1) THEN FD = .TRUE. ELSE FD = .FALSE. END IF print*,'24' CALL Write_array ('rcmetac_inp', ii, & Label, ARRAY_2D = etac, SETUP = FD, ASCI=asci_flag) ! ! Now determine position of the plasmapause: ! r_pp = 5.6 - 0.46*Get_kp (ikptime, kpinput, jtime=0) kc = 1 DO j = 1, jsize DO i = 1, isize IF ( rmin(i,j) <= r_pp) THEN eeta(i,j,kc) = 6.38E+21 * vm(i,j)**(-1.5)* & Plasmasphere_den_CA92_sat (rmin(i,j), doy, sunspot_number) ELSE IF (rmin(i,j) < 6.6) THEN eeta(i,j,kc) = 0.8*(2.26E+21*rmin(i,j)**(-0.5) + & 1.92E+17*rmin(i,j)**4*(1.0-EXP(0.1*(2-rmin(i,j))))) ELSE eeta(i,j,kc) = 0.8*(2.26E+21*6.6**(-0.5) + & 1.92E+17*6.6**4*(1.0-EXP(0.1*(2-6.6)))) END IF IF (eeta(i,j,kc) > eta(1)) THEN eta (1) = eeta (i,j,kc) END IF END DO END DO eta (1) = -eta(1) ! "inverted" edge !!!!!!!!!!!!!!!1 ! IF (eta (1) < -1.0E+30) STOP 'ETA(1) TROUBLE' ! ! ! ! Set up inner edges: ! n_pt = 0 ! point counter DO k = 1, ksize mpoint(k) = n_pt + 1 Nn=0 DO j = jwrap, jsize-2, 2 n_pt = n_pt + 1 Nn = Nn + 1 itrack(n_pt) = n_pt IF (k /= 1) THEN bi(n_pt) = bi_initial_edges(j,k) bj(n_pt) = REAL(j,rprec) etab(n_pt) = eta(k) ! ELSE ! PLASMAPAUSE LOCALTION, DERIVED FROM GRID-BASED K=1 CHANNEL: ! bi (n_pt) = -999.0 bj (n_pt) = REAL(j,rprec) do_find_pp: DO i = isize, 2, -1 IF (rmin(i,j) < r_pp .AND. rmin(i-1,j) >= r_pp) THEN bi (n_pt) = REAL(i,rprec) - 0.5 EXIT do_find_pp END IF END DO do_find_pp IF (bi(n_pt) < 0.0) STOP 'NEGATIVE BI' ! END IF x_pt = Gntrp (xmin, bi(n_pt), bj(n_pt), 0) y_pt = Gntrp (ymin, bi(n_pt), bj(n_pt), 0) IF (k == 1) THEN write (*,'(2(TR1,i4),2(TR2,F5.2),2(TR2,F8.2))') & k, n_pt, bi(n_pt), bj(n_pt), x_pt, y_pt END IF END DO npoint(k) = Nn END DO mpoint(ksize+1) = n_pt+1 ! ! NOTE: itrack(nptmax) should be set to a large non-zero ! value, since it will be used by the RCM function ! NCODE to assign unique itrack values to new points ! added during a run (St) ! ! ! ! ! 7. Compute pressure for energy channels: ! p = (2/3)*alam*eta*vm**(5/2)*2.51*e-34 [dynes/cm**2] ! DO k = 1, ksize p_chan(k) = ABS(alam(k))*eta(k)*vm_bnd(27)**2.5*(2./3.) p_chan(k) = p_chan(k) * 2.51E-34 END DO ! DO k = 1, kcsize DO j = 1, jsize pc_chan(k,j) = ABS(alamc(k))*etac(k,j)*vm_bnd(j)**2.5*(2./3.) pc_chan(k,j) = pc_chan(k,j) * 2.51E-34 END DO END DO p_chan = p_chan / 10.0 / 1.0E-9 pc_chan= pc_chan/ 10.0 / 1.0E-9 ! convert to nPa ! print*, 'pc_chan=',pc_chan ! ! ! !|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! ! Compute plasma properties for each species on edges and on grid: ! ! ** U is energy density in erg/cm-3 ! ** Eavg is average particle density in eV: = U/density ! ** kT is (2/3)*Eavg, in eV ! ** pressure p = 2/3*U, in erg/cm-3 ! eta_numeric_edge (:) = zero eta_numeric_grid (:,:) = zero pvg_numeric_edge (:) = zero pvg_numeric_grid (:,:) = zero ! DO ie = 1, iesize DO k = kmin(ie), kmax(ie) DO j = 1, jsize IF (k /= 1) THEN ! avoid adding plasmasphere eta_numeric_edge(ie) = eta_numeric_edge(ie) + eta (k) eta_numeric_grid(ie,j) = eta_numeric_grid(ie,j) + etac (k,j) ! if(ii.eq.30.and.j.eq.39.and.k.ge.81.and.k.le.160) print*,'k=',k, '2nd etac(k,j)',etac(k,j) pvg_numeric_edge(ie) = pvg_numeric_edge(ie)+2.0/3.0*ABS(alam (k))*eta (k) pvg_numeric_grid(ie,j) = pvg_numeric_grid(ie,j)+2.0/3.0*ABS(alamc(k))*etac(k,j) END IF END DO END DO END DO ! DO ie = 1, iesize DO k = kmin(ie), kmax(ie) DO j = 1, jsize IF (k /= 1) THEN eta_numeric(ie,j) = eta_numeric_grid(ie,j) + eta_numeric_edge(j) pvg_numeric(ie,j) = pvg_numeric_grid(ie,j) + pvg_numeric_edge(j) END IF END DO END DO END DO ! DO ie = 1, iesize DO k = kmin(ie), kmax(ie) DO j = 1, jsize IF (k/=1) THEN eta_numeric_edge(ie) = eta_numeric_edge(ie) / eta_numeric(ie,27) * 100.0 eta_numeric_grid(ie,j) = eta_numeric_grid(ie,j) / eta_numeric(ie,j) * 100.0 pvg_numeric_edge(ie) = pvg_numeric_edge(ie) / pvg_numeric(ie,27) * 100.0 pvg_numeric_grid(ie,j) = pvg_numeric_grid(ie,j) / pvg_numeric(ie,j) * 100.0 END IF END DO END DO END DO ! DO ie = 1, iesize DO j= 1, jsize N_numeric(ie,j) = eta_numeric(ie,j) * (vm_bnd(j)**1.5)*1.E-9/Re/1.E+3/1.E+6 U_numeric(ie,j) = 1.5*2.51E-34*(vm_bnd(j)**2.5)* pvg_numeric(ie,j) kT_numeric(ie,j) = vm_bnd(j) * pvg_numeric(ie,j) / eta_numeric(ie,j) ! in eV E_avg_numeric(ie,j) = 1.5 * kT_numeric(ie,j) P_numeric(ie,j) = (two/three)*U_numeric(ie,j) /10.0 /1.0E-9 ! in nPa if(ii.ge.200)print*,'ie=',ie,' j=',j,' n=',n_numeric(ie,j),kT_numeric(ie,j) END DO END DO ! ! WRITE (*,*) 'PLASMA INITIAL CONDITIONS, T=0:' CALL Print_plasma_summary_info (6) ! to screen print*,'3' ! ! !............................................................ !!1!!!1!!!! start writing initial conditions!!!! if (ii.eq.1) then WRITE (*,'(T2,A)',ADVANCE='NO') & 'PRESS ENTER TO WRITE FILES, CTRL-C TO QUIT WITHOUT WRITING:' READ (*,*) ! ! ! ! Write plasma information: ! CALL Write_plasma () ! ! ! ! Write initial plasma conditions: ! Lrec = 1 label%intg (1) = Lrec CALL Write_array ('rcmbi_inp', Lrec, Label, ARRAY_1D=bi, & SETUP = .TRUE., ASCI=asci_flag) CALL Write_array ('rcmbj_inp', Lrec, Label, ARRAY_1D=bj, & SETUP = .TRUE., ASCI=asci_flag) CALL Write_array ('rcmetab_inp', Lrec, Label, ARRAY_1D=etab, & SETUP = .TRUE., ASCI=asci_flag) CALL Write_array ('rcmitrack_inp', Lrec, Label, ARRAY_1D=itrack,& SETUP = .TRUE., ASCI=asci_flag) CALL Write_array ('rcmmpoint_inp', Lrec, Label, ARRAY_1D=mpoint,& SETUP = .TRUE., ASCI=asci_flag) CALL Write_array ('rcmnpoint_inp', Lrec, Label, ARRAY_1D=npoint,& SETUP = .TRUE., ASCI=asci_flag) CALL Write_array ('rcmeeta_inp', Lrec, Label, ARRAY_3D = eeta,& SETUP = .TRUE., ASCI=asci_flag) print*,'4' ! ! !============================================================= !=== Print out a summary of the plasma setup: =============== ! OPEN (UNIT = LUN_1, FILE = 'setup_plasma_1.list', STATUS = 'REPLACE') CALL date_and_time (char_date, char_time) ! ! A. Write a header: ! WRITE (LUN_1,901) WRITE (LUN_1,901) WRITE (LUN_1,'(/,T16,A)') 'SUMMARY OF RCM INITIAL PLASMA SETUP' WRITE (LUN_1,'(T16,A6,A4,A1,A2,A1,A2,TR5,A6,A2,A1,A2,A1,A2,/)') & 'DATE: ', char_date(1:4),'/',char_date(5:6),'/',char_date(7:8),& 'TIME: ', char_time(1:2),':',char_time(3:4),':',char_time(5:6) ! ! ! B. Write out inner edges info: ! WRITE (LUN_1,'(A)') ' 1. PLASMA DESCRIBED by INNER EDGES:' WRITE (LUN_1,901) WRITE (LUN_1, 903, ADVANCE = 'NO') 903 FORMAT (T3,'K', T9,'ALAM(K)', T19,'W @boundRe', T31,'ETA(K)', & T43,'IE', T49,'% ETA(IE)', T60, '% P(IE)', T70,' Fudge') WRITE (LUN_1, 902) 902 FORMAT (T6,'|',T17,'|',T28,'|',T39,'|',T47,'|', T58,'|', T68,'|', T76, '|') WRITE (LUN_1,901) ! ! print*, 'eta_numeric=',eta_numeric(ikflav(k),j) DO k = 1, ksize DO j = 1, jsize ! print*, 'eta_numeric=',eta_numeric(ikflav(k),j) p_tmp = p_chan (k) / p_numeric(ikflav(k),j)*100.0 eta_tmp = eta (k) / eta_numeric(ikflav(k),j)*100.0 WRITE (LUN_1, 904, ADVANCE='NO') & k, alam(k), alam(k)*vm_bnd(j), eta(k), ikflav(k), & eta_tmp, p_tmp, remark_edge(k) WRITE (LUN_1, 902) END DO END DO WRITE (LUN_1, '(T28,ES10.2)', ADVANCE='NO') SUM(eta) WRITE (LUN_1, 902) 904 FORMAT (I3, T8,F8.1, T19,F8.1, T28,ES10.2, T42,I2, T49,F6.1, T60,F6.1, T70,A) WRITE (LUN_1, 901) ! ! ! C. Write out grid plasma info: ! WRITE (LUN_1,'(A)') ' 2. PLASMA DESCRIBED ON THE GRID:' WRITE (LUN_1,901) WRITE (LUN_1, 903, ADVANCE = 'NO') WRITE (LUN_1, 902) WRITE (LUN_1,901) ! DO k = 1, kcsize DO j = 1, jsize p_tmp = pc_chan (k,j) / p_numeric(ikflav(k),j)*100.0 eta_tmp = etac (k,j) / eta_numeric(ikflav(k),j)*100.0 ! print*,'p_tmp=',p_tmp ! print*,'eta_tmp=',eta_tmp ! print*,'etac=', etac WRITE (LUN_1, 904, ADVANCE= 'NO') & k, alamc (k), alamc (k)*vm_bnd(j), etac (k,j), ikflavc(k), & eta_tmp, p_tmp, remark_grid (k) WRITE (LUN_1, 902) END DO END DO WRITE (LUN_1, '(T28,ES10.2)' ) SUM(etac) WRITE (LUN_1, 901) ! WRITE (LUN_1,*) 'PLASMA INITIAL CONDITIONS, T=0:' CALL Print_plasma_summary_info (LUN_1) ! ! OPEN (UNIT = LUN_2, FILE = 'setup_plasma_2.list', STATUS = 'REPLACE') WRITE (LUN_2,901) WRITE (LUN_2,901) 901 FORMAT & (' ---------------------------------------------------------------------------') ! ! ! ! ! Print out plasma arrays to text file 'setup_plasma.list': ! Currently, print out only the plasmasphere/plasmapause channel: ! CALL Outbij (K_INIT = 1, K_FINAL = 1 , NCOL = 80, NTP = LUN_2) ! DO kc = 1, 1 CALL Outp (R = eeta(:,:,kc), IBEG = 1, IEND = isize, IINC = 1, & JBEG = jwrap, JEND = jsize-1, JINC = 1, & XSCALE = zero, ILABEL = label%intg, & TITLE = title_string, NTP = LUN_2, NCOL = 80) END DO ! CLOSE (UNIT = LUN_2) ! ! end if !!!!!!! end writing initial conditions!!!!! END DO ! time loop ! END IF WRITE (LUN_1,'(//,T10,A//)') 'END OF RCM INTIAL PLASMA SETUP SUMMARY' CLOSE (LUN_1) ! STOP ! CONTAINS ! ! SUBROUTINE Print_plasma_summary_info (L) INTEGER (iprec) :: L ! DO ie = 1, iesize DO ie = 1, 2 print*,'time=',ii,' test1' WRITE (L,'(T2,A)') '--------------------------------------------------------------------' WRITE (L,'(T2,A,A)') ' '//species_char(ie), ' SETUP INFORMATION: ' WRITE (L,'(T2,A)') '--------------------------------------------------------------------' WRITE (L,'(T2,A)') ' PARAM, UNITS | ACTUAL (NUMERIC) | ANALYTIC ' WRITE (L,'(T2,A)') '--------------------------------------------------------------------' WRITE (L, 801) 'Kp ', '|', Kp, '|', Kp WRITE (L, 801) 'Kp_for_com', '|', Kp_for_composition, '|', Kp_for_composition WRITE (L, 801) 'TOTAL ETA, Wb-1', '|', eta_numeric (ie,33), '|', eta_an(ie,33) WRITE (L, 803) ' -EDGE-BASED, %', '|', eta_numeric_edge(ie),'|', 'N/A' WRITE (L, 803) ' -GRID-BASED, %', '|', eta_numeric_grid(ie,33),'|', 'N/A' WRITE (L, 801) 'N(-20,0,0), cm-3','|', n_numeric(ie,39), '|', N_c(ie,39)+N_h(ie,39) WRITE (L, 801) 'N(-20,0,0), cm-3','|', n_numeric(ie,27), '|', N_c(ie,27)+N_h(ie,27) WRITE (L, 801) 'TOTAL Pv**gamma ','|', pvg_numeric(ie,33), '|'!, Pvg_an (ie,33) WRITE (L, 803) ' -EDGE-BASED, %', '|', pvg_numeric_edge(ie),'|', 'N/A' WRITE (L, 803) ' -GRID-BASED, %', '|', pvg_numeric_grid(ie,33),'|', 'N/A' WRITE (L, 801) '(-20,0,0), eV','|', E_avg_numeric(ie,27), '|', E_avg_c(ie,27) WRITE (L, 801) 'kT (-20,0,0), eV','|', kT_numeric(ie,27), '|', kT(ie,27) write (L, 801) 'kT (-20,0,0), eV','|', kT_numeric(ie,39), '|', kT(ie,39) WRITE (L, 801) 'kT (-20,0,0), eV','|', kT_numeric(ie,15), '|', kT(ie,15) WRITE (L, 801) 'P(-20,0,0), nPa ','|', P_numeric(ie,27), '|'!, P(ie,27) WRITE (L,'(T2,A)') '--------------------------------------------------------------------' 801 FORMAT (T2,A, T21,A1, T30,ES9.2, T46,A1, T53,ES9.2) 803 FORMAT (T2,A, T21,A1, T33, F6.2, T46,A1, T59,A) DO j=15,15 if (ie.eq.2)then endif END DO END DO ! RETURN END SUBROUTINE Print_plasma_summary_info ! ! END PROGRAM Setup_plasma