JHUAPL/kaiju
1
0
Fork 0
mirror of https://github.com/JHUAPL/kaiju.git synced 2026-01-09 18:57:53 -05:00
Code Issues Packages Projects Releases 1 Wiki Activity

First wave of edits. Partially tested. Removed build folder

Browse Source
This commit is contained in:
Anthony
2020-04-17 09:45:06 -06:00
parent b86fce9004
commit 6102a7a065
11 changed files with 63 additions and 35 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
.gitignore vendored
View File

@@ -1,5 +1,6 @@
# User dependent stuff
cmake/user.cmake
build
# Compiled Object files
*.slo
@@ -39,4 +40,4 @@ cmake/user.cmake
*.h5part
# Mac nonsense
*.DS_STORE
*.DS_STORE

4
src/base/types/volttypes.F90
View File

@@ -135,10 +135,12 @@ module volttypes
contains
! null default subroutines for inner mag base type
subroutine baseInit(imag,iXML,isRestart,vApp)
subroutine baseInit(imag,iXML,isRestart,Rp_m,iono_m,vApp)
class(innerMagBase_T), intent(inout) :: imag
type(XML_Input_T), intent(in) :: iXML
logical, intent(in) :: isRestart
real(rp), intent(in) :: Rp_m
real(rp), intent(in) :: iono_m
type(voltApp_T), intent(inout) :: vApp
end subroutine

1
src/rcm/modules.f90
View File

@@ -144,6 +144,7 @@ MODULE ionosphere_exchange
rm%nLat_ion = isize
rm%nLon_ion = jsize-jwrap+1
rm%planet_radius = 6380.e3 ! Default to Earth radius
ALLOCATE( rm%gcolat(rm%nLat_ion) )
ALLOCATE( rm%glong(rm%nLon_ion) )

3
src/rcm/rcm_mhd.f90
View File

@@ -122,7 +122,8 @@ subroutine rcm_mhd(mhdtime,mhdtimedt,RM,iflag,iXML)
! Set up RCM ionospheric grid:
!call Grid_torcm (75.0_rprec, 15.0_rprec, 0.0_rprec, radius_earth_m, radius_iono_m) ! set up RCM ionospheric grid here
call Grid_torcm (HighLatBD,LowLatBD, 0.0_rprec, radius_earth_m, radius_iono_m) ! set up RCM ionospheric grid here
!call Grid_torcm (HighLatBD,LowLatBD, 0.0_rprec, radius_earth_m, radius_iono_m) ! set up RCM ionospheric grid here
call Grid_torcm (HighLatBD,LowLatBD, 0.0_rprec, RM%planet_radius, RM%iono_radius) ! set up RCM ionospheric grid here
! Setup Ionosphere intermediate Grid by equating it to the RCM grid, without angular overlap:
call setupIon(RM)

2
src/rcm/rcm_mhd_interfaces.f90
View File

@@ -13,6 +13,8 @@ module rcm_mhd_interfaces
type rcm_mhd_T
integer(iprec) :: nLat_ion
integer(iprec) :: nLon_ion
real(rprec) :: planet_radius !Planet radius in meters
real(rprec) :: iono_radius !Ionosphee radius in meters
real(rprec),allocatable :: gcolat(:) !> RCM Latitude grid points
real(rprec),allocatable :: glong(:) !> RCM Longitude grid points

15
src/rcm/tomhd.f90
View File

@@ -10,7 +10,7 @@
Read_grid, Read_plasma,Get_boundary, &
xmass, densrcm,denspsph,imin_j,rcmdir, &
eflux,eavg,ie_el
USE constants, ONLY : mass_proton,radius_earth_m,nt,ev,pressure_factor,density_factor
USE constants, ONLY : mass_proton,nt,ev!,radius_earth_m,pressure_factor,density_factor
USE rice_housekeeping_module
Use rcm_mhd_interfaces
!
@@ -92,6 +92,9 @@
LOGICAL,PARAMETER :: avoid_boundaries = .false.
INTEGER(iprec) :: im,ipl,jm,jpl,km,kpl
!Replacing constant density_factor and pressure_factor to allow for other planet radii
!real(rprec) :: pressure_factor = 2./3.*ev/RM%planet_radius*nt
!real(rprec) :: density_factor = nt/RM%planet_radius
! INCLUDE 'rcmdir.h'
ierr = 0
@@ -210,14 +213,15 @@
IF (vm(i,j) < 0.0) CYCLE
DO k = 1, kcsize
pressrcm(i,j) = pressrcm(i,j) + &
pressure_factor*ABS(alamc(k))*eeta_avg(i,j,k)*vm(i,j)**2.5 ! in pascals
!pressure_factor*ABS(alamc(k))*eeta_avg(i,j,k)*vm(i,j)**2.5 ! in pascals
2./3.*ev/RM%planet_radius*nt*ABS(alamc(k))*eeta_avg(i,j,k)*vm(i,j)**2.5 ! in pascals
! normalize everything to the mass_proton, otherwise answer is below
! floating point minimum answer and gets zero in ples/m^3
! FIXME: This version is mass weighted, not sure why.
if(alamc(k) >0.0)then ! only add the ion contribution
densrcm(i,j) = densrcm(i,j) + &
density_factor/mass_proton*xmass(ikflavc(k))*eeta_avg(i,j,k)*vm(i,j)**1.5
!density_factor/mass_proton*xmass(ikflavc(k))*eeta_avg(i,j,k)*vm(i,j)**1.5
nt/RM%planet_radius/mass_proton*xmass(ikflavc(k))*eeta_avg(i,j,k)*vm(i,j)**1.5
end if
END DO
if (use_plasmasphere) then
@@ -244,7 +248,8 @@
RM%fac = birk (:,jwrap:jdim)
RM%toMHD = .false.
dRad = ellBdry%dRadMHD*radius_earth_m
!dRad = ellBdry%dRadMHD*radius_earth_m
dRad = ellBdry%dRadMHD*RM%planet_radius
do j=jwrap,jdim
jp = j-jwrap+1

24
src/rcm/torcm.f90
View File

@@ -192,7 +192,7 @@
CALL Read_alam (kcsize, alamc, ikflavc, fudgec, almdel, almmax, almmin, iesize, ierr)
IF (ierr < 0) RETURN
CALL Press2eta ! this populates EETA_NEW array
CALL Press2eta(RM%planet_radius) ! this populates EETA_NEW array
! Here is the second IF promised at the beginning.
! It must come after calls to calc_ftv and press2eta.
@@ -328,7 +328,7 @@
use rcm_mhd_interfaces
IMPLICIT NONE
type(rcm_mhd_T),intent(in) :: RM
type(rcm_mhd_T), intent(in) :: RM
REAL(rprec), INTENT (IN) :: big_vm
INTEGER(iprec), INTENT (OUT) :: ierr
!
@@ -411,15 +411,15 @@
den (:, 1) = den (:,jsize-2)
den (:, 2) = den (:,jsize-1)
xmin (:,jwrap:jsize) = RM%x_bmin (:, :,1)/radius_earth_m
xmin (:,jwrap:jsize) = RM%x_bmin (:, :,1)/RM%planet_radius
xmin (:, 1) = xmin (:,jsize-2)
xmin (:, 2) = xmin (:,jsize-1)
ymin (:,jwrap:jsize) = RM%x_bmin (:, :,2)/radius_earth_m
ymin (:,jwrap:jsize) = RM%x_bmin (:, :,2)/RM%planet_radius
ymin (:, 1) = ymin (:,jsize-2)
ymin (:, 2) = ymin (:,jsize-1)
zmin (:,jwrap:jsize) = RM%x_bmin (:, :,3)/radius_earth_m
zmin (:,jwrap:jsize) = RM%x_bmin (:, :,3)/RM%planet_radius
zmin (:, 1) = zmin (:,jsize-2)
zmin (:, 2) = zmin (:,jsize-1)
@@ -519,7 +519,7 @@
!====================================================================
!
SUBROUTINE Press2eta
SUBROUTINE Press2eta(planet_radius)
! ====================================================================
!
@@ -571,11 +571,12 @@
IMPLICIT NONE
real(rprec):: xmin,xmax
real(rprec):: ptemp,eta_correction
real(rprec), intent(in) ::planet_radius
!
! real(rprec):: Erf
! EXTERNAL Erf
!
real(rprec) ::trans = radius_earth_m/nt
!real(rprec) ::trans = radius_earth_m/nt
real(rprec),parameter :: eps=1.0e-30
real(rprec):: sqrtpi,factor,fac0,fac1,t
integer(iprec) :: i,j,k
@@ -609,7 +610,8 @@
STOP 'ILLEGAL IKFLAVC(K) IN PRESS2ETA'
END IF
if (t > 0.)then
fac0 = trans*den(i,j)/((vm(i,j))**1.5)
!fac0 = trans*den(i,j)/((vm(i,j))**1.5)
fac0 = planet_radius/nt*den(i,j)/((vm(i,j))**1.5)
xmax = SQRT(ev*ABS(almmax(k))*vm(i,j)/boltz/t)
xmin = SQRT(ev*ABS(almmin(k))*vm(i,j)/boltz/t)
fac1 = (Erf(xmax)-Erf(xmin)) -2.0/sqrtpi* &
@@ -631,7 +633,8 @@
ptemp = 0.0
DO k=1,kdim
ptemp = ptemp + &
pressure_factor*ABS(alamc(k))*eeta_new(i,j,k)*vm(i,j)**2.5
!pressure_factor*ABS(alamc(k))*eeta_new(i,j,k)*vm(i,j)**2.5
2./3.*ev/planet_radius*nt*ABS(alamc(k))*eeta_new(i,j,k)*vm(i,j)**2.5
END DO
eta_correction = 0.0
IF (ptemp/= 0)eta_correction = press(i,j)/ptemp
@@ -642,7 +645,8 @@
ptemp = 0.0
DO k=1,kdim
ptemp = ptemp + &
pressure_factor*ABS(alamc(k))*eeta_new(i,j,k)*vm(i,j)**2.5
!pressure_factor*ABS(alamc(k))*eeta_new(i,j,k)*vm(i,j)**2.5
2./3.*ev/planet_radius*nt*ABS(alamc(k))*eeta_new(i,j,k)*vm(i,j)**2.5
END DO
IF (ABS(ptemp-press(i,j)) > 0.01*press(i,j))then
Write(6,*)' Warning, pressures do not match at i,j =',i,j

11
src/voltron/innermagsphere.F90
View File

@@ -29,12 +29,12 @@ module innermagsphere
contains
!Figure out which inner magnetosphere model we're using and initialize it
subroutine InitInnerMag(vApp,isRestart,iXML)
subroutine InitInnerMag(vApp,gApp,iXML)
type(voltApp_T) , intent(inout) :: vApp
logical, intent(in) :: isRestart
type(gamApp_T), intent(in) :: gApp
type(XML_Input_T), intent(inout) :: iXML
character(len=strLen) :: imStr
real(rp) :: iono_rad_m
if (.not. vApp%doDeep) return !Why are you even here?
@@ -55,7 +55,10 @@ module innermagsphere
stop
end select
call vApp%imagApp%doInit(iXML,isRestart,vApp)
!iono_rad_m = gApp%Model%Units%Rion*gApp%Model%Units%gx0 !Convert Rp to m
!Assume for now that ionosphere is always 1.01*planet radius
iono_rad_m = 1.01*gApp%Model%Units%gx0 ! m
call vApp%imagApp%doInit(iXML,gApp%Model%isRestart,gApp%Model%Units%gx0,iono_rad_m,vApp)
end subroutine InitInnerMag

29
src/voltron/rcmimag.F90
View File

@@ -30,6 +30,7 @@ module rcmimag
integer, parameter :: MAXRCMIOVAR = 30
character(len=strLen), private :: h5File
real(rp), private :: Rp_cgs
!Information taken from MHD flux tubes
!TODO: Figure out RCM boundaries
@@ -79,12 +80,13 @@ module rcmimag
contains
!Initialize RCM inner magnetosphere model
subroutine initRCM(imag,iXML,isRestart,vApp)
subroutine initRCM(imag,iXML,isRestart,Rp_m,iono_m,vApp)
class(rcmIMAG_T), intent(inout) :: imag
type(XML_Input_T), intent(in) :: iXML
logical, intent(in) :: isRestart
real(rp), intent(in) :: Rp_m
real(rp), intent(in) :: iono_m
type(voltApp_T), intent(inout) :: vApp
character(len=strLen) :: RunID,RCMH5
logical :: fExist
@@ -93,7 +95,6 @@ module rcmimag
t0 => vApp%time, &
dtCpl => vApp%DeepDT, &
nRes => vApp%IO%nRes)
call iXML%Set_Val(RunID,"/gamera/sim/runid","sim")
RCMApp%rcm_runid = trim(RunID)
@@ -124,6 +125,12 @@ module rcmimag
call iXML%Set_Val(SmoothOp%nIter,"imag/nIter",4)
call iXML%Set_Val(SmoothOp%nRad ,"imag/nRad" ,8)
Rp_cgs = Rp_m*1.0e+2
!Store planet radius in rcmCpl so that units stay consistent between models
RCMApp%planet_radius = Rp_m
RCMApp%iono_radius = iono_m
write(*,*) 'RCM planet radius = ',RCMApp%planet_radius
write(*,*) 'RCM ionosphere radius = ',RCMApp%iono_radius
end associate
end subroutine initRCM
@@ -145,7 +152,7 @@ module rcmimag
associate(RCMApp => imag%rcmCpl)
!Lazily grabbing rDeep here, convert to RCM units
rEqMin = vApp%rDeep*Re_cgs*1.0e-2 !Re=>meters
rEqMin = vApp%rDeep*Rp_cgs*1.0e-2 !Rp=>meters
llBC = vApp%mhd2chmp%lowlatBC
@@ -221,14 +228,14 @@ module rcmimag
!Find maximum extent of closed field region
maxRad = maxval(norm2(RCMApp%X_bmin,dim=3),mask=vApp%imag2mix%isClosed)
maxRad = maxRad/(Re_cgs*1.0e-2)
maxRad = maxRad/(Rp_cgs*1.0e-2)
vApp%rTrc = 1.25*maxRad
end associate
contains
!Calculate Alfven bounce timescale
!D = #/m3, B = T, L = Re
!D = #/m3, B = T, L = Rp
function AlfvenBounce(D,B,L) result(dTb)
real(rp), intent(in) :: D,B,L
real(rp) :: dTb
@@ -242,7 +249,7 @@ module rcmimag
nCC = D*rcmNScl !Get n in #/cc
bNT = B*1.0e+9 !Convert B to nT
Va = 22.0*bNT/sqrt(nCC) !km/s, from NRL plasma formulary
dTb = (L*Re_km)/Va
dTb = (L*Rp_cgs*(1.0e-5))/Va !km
end function AlfvenBounce
end subroutine AdvanceRCM
@@ -435,10 +442,10 @@ module rcmimag
call genStream(ebModel,ebState,x0,t,bTrc)
!Get diagnostics from field line
!Minimal surface (bEq in Re, bMin in EB)
!Minimal surface (bEq in Rp, bMin in EB)
call FLEq(ebModel,bTrc,bEq,bMin)
bMin = bMin*oBScl*1.0e-9 !EB=>Tesla
bEq = bEq*Re_cgs*1.0e-2 !Re=>meters
bEq = bEq*Rp_cgs*1.0e-2 !Rp=>meters
!Plasma quantities
!dvB = Flux-tube volume (Re/EB)
@@ -498,8 +505,8 @@ module rcmimag
colat = PI/2 - lat
L = 1.0/(sin(colat)**2.0)
ijTube%Vol = 32./35.*L**4.0/mdipole
ijTube%X_bmin(XDIR) = L*cos(lon)*Re_cgs*1.0e-2 !Re=>meters
ijTube%X_bmin(YDIR) = L*sin(lon)*Re_cgs*1.0e-2 !Re=>meters
ijTube%X_bmin(XDIR) = L*cos(lon)*Rp_cgs*1.0e-2 !Rp=>meters
ijTube%X_bmin(YDIR) = L*sin(lon)*Rp_cgs*1.0e-2 !Rp=>meters
ijTube%X_bmin(ZDIR) = 0.0
ijTube%bmin = mdipole/L**3.0

4
src/voltron/sstimag.F90
View File

@@ -34,10 +34,12 @@ module sstimag
contains
!Initialize EQ Map data
subroutine initSST(imag,iXML,isRestart,vApp)
subroutine initSST(imag,iXML,isRestart,Rp_m,iono_m,vApp)
class(eqData_T), intent(inout) :: imag
type(XML_Input_T), intent(in) :: iXML
logical, intent(in) :: isRestart !Do you even care?
real(rp), intent(in) :: Rp_m !Still don't care
real(rp), intent(in) :: iono_m !Really don't care
type(voltApp_T), intent(inout) :: vApp
character(len=strLen) :: eqFile

2
src/voltron/voltapp.F90
View File

@@ -139,7 +139,7 @@ module voltapp
!Set first deep coupling (defaulting to 0)
call xmlInp%Set_Val(vApp%DeepT, "coupling/tDeep", 0.0_rp)
!Initialize deep coupling type/inner magnetosphere model
call InitInnerMag(vApp,gApp%Model%isRestart,xmlInp)
call InitInnerMag(vApp,gApp,xmlInp)
endif
if(present(optFilename)) then