diff --git a/.fprettify.rc b/.fprettify.rc index 8052f7ca..9d14c90e 100644 --- a/.fprettify.rc +++ b/.fprettify.rc @@ -9,7 +9,7 @@ i: 2 recursive # Excluded files: # Anything compiled, src/venus.f90 (it's incomplete), old fortran files, makefiles, template files -exclude: ["*.exe","*.o","Venus.f90","*.f","Makef*","*emplate.f90"] +exclude: ["*.exe","*.o","Venus.f90","*.f","Makef*","*emplate.f90","*Msis21.F90"] # Whitespace option 2 (check fprettify -h for info): whitespace: 2 diff --git a/.github/workflows/Format-Doc-Test.yml b/.github/workflows/Format-Doc-Test.yml index 42cfacfd..dd6938bb 100644 --- a/.github/workflows/Format-Doc-Test.yml +++ b/.github/workflows/Format-Doc-Test.yml @@ -26,7 +26,8 @@ jobs: pip install . - name: Clone GITM uses: actions/checkout@v4 - with: + with: + ref: ${{ github.event.pull_request.head.sha }} path: GITM - name: Format check entire GITM repository diff --git a/Makefile b/Makefile index 8b66be20..6ae6b49f 100644 --- a/Makefile +++ b/Makefile @@ -61,7 +61,6 @@ LIB: cd $(GLDIR) ; make LIBPREV=${LIBDIR}/libSphere.a LIBADD cd $(MAINDIR) ; make LIBPREV=${GITM}/${GLDIR}/libUPTOGL.a LIB cd srcInterface ; make LIBPREV=${GITM}/${MAINDIR}/libUA.a LIB - make POST nompirun: make GITM @@ -109,12 +108,10 @@ rundir: if [ ! -e "EIE/README" ]; then \ ln -s ${EMPIRICALIEDIR}/data EIE;\ fi; - cd ${RUNDIR}; rm -f ./PostGITM.exe ; ln -s ${UADIR}/src/PostProcess.exe ./PostGITM.exe cd ${RUNDIR}/UA; \ mkdir restartOUT data DataIn; \ ln -s restartOUT restartIN; \ - ln -s ${UADIR}/src/pGITM .; \ - ln -s ${UADIR}/srcPython/pGITM.py .; \ + ln -s ${UADIR}/srcPython/post_process.py .; \ ln -s ${UADIR}/srcData/* DataIn; rm -f DataIn/CVS; \ ln -s ${UADIR}/data/* DataIn; rm -f DataIn/CVS cd ${RUNDIR} ; \ diff --git a/README.md b/README.md index be572f4b..c5106b60 100644 --- a/README.md +++ b/README.md @@ -115,25 +115,24 @@ blocks with 9 x 9 cells each, so the default resolution is 180 (deg lat) / (2 \* 9) = 10 deg lat, by 360 (deg lon) / (2 \* 9) = 20 deg lon. See below for how to set the resolution. -8\. Go into the directory which contains many of the outputs: +8\. Post process the output files by running: ```shell -cd UA +./post_process.py ``` -9\. Post process the output files by running: +> This can be called from `run/` and will postprocess the files in UA/output. +> It also has functionality to copy files to a remote location, monitor the +> output folder throughout a run, and more. Run `./post_process.py --help` to +> see available options. -```shell -./pGITM -``` - -10\. Go into the output directory: +9\. Go into the output directory: ```shell cd data ``` -11\. Make some plots with an old plotter: +10\. Make some plots with an old plotter: ```shell ../../../srcPython/plot_model_results.py -var=3 -alt=120 3DALL_t021221_000500.bin @@ -142,7 +141,7 @@ cd data Then look at the png file that is created. You can use a `-h` to see how to run this code. -11b\. A more advanced plotter is available through aetherpy. This is a bit more +10b\. A more advanced plotter is available through aetherpy. This is a bit more complicated, since you need to install aetherpy. If you don't use python much, this is harder. Here is how to do this: @@ -165,7 +164,7 @@ git checkout develop python setup.py develop --user ``` -11c\. Test out the new plotter: +10c\. Test out the new plotter: ```shell cd diff --git a/src/Makefile b/src/Makefile index 1922f6f1..b79819d6 100644 --- a/src/Makefile +++ b/src/Makefile @@ -145,7 +145,6 @@ SAMIEXE = ${BINDIR}/GITMSAMI.exe GITM: DEPEND @make ${EXE} - @make POST @echo @echo "${EXE} has been created" @echo @@ -158,7 +157,6 @@ ${EXE}: ${OBJECTS_EXE} SAMI: DEPENDSAMI @make ${SAMIEXE} - @make POST @echo "${SAMIEXE} has been created" ${SAMIEXE}: ${OBJECTS_SAMI} diff --git a/src/ModInputs.f90 b/src/ModInputs.f90 index 17a586f4..44b40c71 100644 --- a/src/ModInputs.f90 +++ b/src/ModInputs.f90 @@ -285,7 +285,7 @@ module ModInputs real :: CO2ppm = 225.0 logical :: DoN4SHack = .false. - + ! Allow the user to change the planet's characteristics: real :: RotationPeriodInput = Rotation_Period real :: OmegaBodyInput = 2.0*pi/Rotation_Period diff --git a/src/init_altitude.f90 b/src/init_altitude.f90 index cf372382..96b1abc5 100644 --- a/src/init_altitude.f90 +++ b/src/init_altitude.f90 @@ -16,8 +16,8 @@ subroutine get_temperature(lon, lat, alt, t, h) !--------------------------------------------------------------------------- if (UseMsis) then - call initialize_msis_routines - call get_msis_temperature(lon, lat, alt, t, h) + call initialize_msis_routines + call get_msis_temperature(lon, lat, alt, t, h) else diff --git a/src/init_msis.Earth.f90 b/src/init_msis.Earth.f90 index 5fed7104..c650cf13 100644 --- a/src/init_msis.Earth.f90 +++ b/src/init_msis.Earth.f90 @@ -29,7 +29,7 @@ subroutine call_msis(lonDeg, latDeg, altKm, f107, f107a, densities10, temp) use ModConstants, only: Boltzmanns_Constant, AMU implicit none - + real, intent(in) :: lonDeg, latDeg, altKm, f107, f107a real, intent(out) :: densities10(10) real, intent(out) :: temp @@ -58,52 +58,52 @@ subroutine call_msis(lonDeg, latDeg, altKm, f107, f107a, densities10, temp) AP = 10 if (useMsis21) then - iyd = iJulianDay - sec = utime - alt = altKm - glat = latDeg - glong = lonDeg - stl = LST - f107a_4 = f107a - f107_4 = f107 - ap_4 = AP - ! mass is not used, but passed anyways - mass = -1 - call gtd8d(iyd,sec,alt,glat,glong,stl,f107a_4,f107_4,ap_4,mass,d,t) - temp = t(2) - ! Convert to /m3 - ! 10th density is NO now! - densities10 = d * 1e6 + iyd = iJulianDay + sec = utime + alt = altKm + glat = latDeg + glong = lonDeg + stl = LST + f107a_4 = f107a + f107_4 = f107 + ap_4 = AP + ! mass is not used, but passed anyways + mass = -1 + call gtd8d(iyd, sec, alt, glat, glong, stl, f107a_4, f107_4, ap_4, mass, d, t) + temp = t(2) + ! Convert to /m3 + ! 10th density is NO now! + densities10 = d*1e6 else - CALL GTD7(iJulianDay, utime, AltKm, LatDeg, LonDeg, LST, & - f107a, f107, AP, 48, msis_dens9, msis_temp) - temp = msis_temp(2) - densities10(1:9) = msis_dens9 - - ! Very old code: - ! ! The initial profile of [NO] is refered to: - ! ! [Charles A. Barth, AGU, 1995] - ! - ! if (geo_alt < 120.) then - ! NDensityS(iLon,iLat,iAlt,iNO_,iBlock)= & - ! max(1e14-1e10*abs((geo_alt-110.0))**3.5, 100.0) - ! !10**(-0.003*(geo_alt-105.)**2 +14+LOG10(3.)) - ! else - ! m = (1e10-3.9e13)/(200) - ! k = 1e10+(-m*300.) - ! NDensityS(iLon,iLat,iAlt,iNO_,iBlock)= & - ! MAX(k+(m*geo_alt)-(geo_alt - 120.0)**2,100.0) - ! ! MAX(10**(13.-LOG10(3.)*(geo_alt-165.)/35.),1.0) - ! endif - ! - ffactor = 6.36*log(f107) - 13.8 - no = (ffactor*1.0e13 + 8.0e13)*1.24 ! 12.4 ! 12.4 is roughly exp - ! This is obviously an approximation: - h = Boltzmanns_Constant * msis_temp(2)/ & - (9.5 * 28.0 * AMU)/1000.0 - densities10(10) = no * exp(-(altKm - 100.0)/h) + CALL GTD7(iJulianDay, utime, AltKm, LatDeg, LonDeg, LST, & + f107a, f107, AP, 48, msis_dens9, msis_temp) + temp = msis_temp(2) + densities10(1:9) = msis_dens9 + + ! Very old code: + ! ! The initial profile of [NO] is refered to: + ! ! [Charles A. Barth, AGU, 1995] + ! + ! if (geo_alt < 120.) then + ! NDensityS(iLon,iLat,iAlt,iNO_,iBlock)= & + ! max(1e14-1e10*abs((geo_alt-110.0))**3.5, 100.0) + ! !10**(-0.003*(geo_alt-105.)**2 +14+LOG10(3.)) + ! else + ! m = (1e10-3.9e13)/(200) + ! k = 1e10+(-m*300.) + ! NDensityS(iLon,iLat,iAlt,iNO_,iBlock)= & + ! MAX(k+(m*geo_alt)-(geo_alt - 120.0)**2,100.0) + ! ! MAX(10**(13.-LOG10(3.)*(geo_alt-165.)/35.),1.0) + ! endif + ! + ffactor = 6.36*log(f107) - 13.8 + no = (ffactor*1.0e13 + 8.0e13)*1.24 ! 12.4 ! 12.4 is roughly exp + ! This is obviously an approximation: + h = Boltzmanns_Constant*msis_temp(2)/ & + (9.5*28.0*AMU)/1000.0 + densities10(10) = no*exp(-(altKm - 100.0)/h) endif - + end subroutine call_msis subroutine get_msis_temperature(lon, lat, alt, t, h) @@ -163,17 +163,17 @@ subroutine get_msis_temperature(lon, lat, alt, t, h) if (RCMRFlag .and. RCMROutType == "F107") then - call call_msis(lonDeg, latDeg, altKm, f107_msis, f107a_msis, msis_dens10, msis_temp1) - msis_dens = msis_dens10(1:9) - msis_temp = msis_temp1 - !CALL GTD7(iJulianDay, utime, AltKm, LatDeg, LonDeg, LST, & - ! f107a_msis, f107_msis, AP, 48, msis_dens, msis_temp) + call call_msis(lonDeg, latDeg, altKm, f107_msis, f107a_msis, msis_dens10, msis_temp1) + msis_dens = msis_dens10(1:9) + msis_temp = msis_temp1 + !CALL GTD7(iJulianDay, utime, AltKm, LatDeg, LonDeg, LST, & + ! f107a_msis, f107_msis, AP, 48, msis_dens, msis_temp) else - call call_msis(lonDeg, latDeg, altKm, f107, f107a, msis_dens10, msis_temp1) - msis_dens = msis_dens10(1:9) - msis_temp = msis_temp1 - !call GTD7(iJulianDay, utime, AltKm, LatDeg, LonDeg, LST, & - ! F107A, F107, AP, 48, msis_dens, msis_temp) + call call_msis(lonDeg, latDeg, altKm, f107, f107a, msis_dens10, msis_temp1) + msis_dens = msis_dens10(1:9) + msis_temp = msis_temp1 + !call GTD7(iJulianDay, utime, AltKm, LatDeg, LonDeg, LST, & + ! F107A, F107, AP, 48, msis_dens, msis_temp) endif t = msis_temp(2) @@ -201,7 +201,7 @@ subroutine initialize_msis_routines use ModInputs, only: UseMsisTides, useMsis21, UseMsisOnly, sw_msis use EUA_ModMsis00, ONLY: meters, tselec use msis_init, only: msisinit - + implicit none real*4 :: sw_msis4x25(25) @@ -209,11 +209,11 @@ subroutine initialize_msis_routines ! We only need to initialize MSIS once! if (.not. isFirstTime) then - return + return else - isFirstTime = .false. + isFirstTime = .false. endif - + ! We want units of /m3 and not /cm3 call meters(.true.) @@ -238,12 +238,12 @@ subroutine initialize_msis_routines sw_msis(2) = 0 if (useMsis21) then - sw_msis4x25 = sw_msis - call msisinit(parmpath = 'UA/DataIn/LowerBCs/', switch_legacy=sw_msis4x25) + sw_msis4x25 = sw_msis + call msisinit(parmpath='UA/DataIn/LowerBCs/', switch_legacy=sw_msis4x25) else - call tselec(sw_msis) + call tselec(sw_msis) endif - + end subroutine initialize_msis_routines !-------------------------------------------------------------- @@ -299,7 +299,7 @@ subroutine init_msis endif call initialize_msis_routines - + !-------------------------------------------------------------------------- ! ! From the msis90 library: @@ -393,10 +393,10 @@ subroutine init_msis NDensityS(iLon, iLat, iAlt, iH_, iBlock) = & max(msis_dens(7), 100.0) NDensityS(iLon, iLat, iAlt, iN_4S_, iBlock) = & - max(msis_dens(8), 100.0) - NDensityS(iLon, iLat, iAlt, iNO_, iBlock) = & - max(msis_dens10(10), 100.0) - + max(msis_dens(8), 100.0) + NDensityS(iLon, iLat, iAlt, iNO_, iBlock) = & + max(msis_dens10(10), 100.0) + NDensityS(iLon, iLat, iAlt, iN_2P_, iBlock) = & NDensityS(iLon, iLat, iAlt, iN_4S_, iBlock)/10000.0 NDensityS(iLon, iLat, iAlt, iN_2D_, iBlock) = & @@ -554,7 +554,6 @@ subroutine msis_bcs(iJulianDay, UTime, Alt, LatIn, LonIn, Lst, & !---------------------------------------------------------------------------- AP_I = AP - !CALL GTD7(iJulianDay, uTime, Alt, Lat, Lon, LST, & ! F107A, F107, AP_I, 48, msis_dens, msis_temp) call call_msis(lon, lat, alt, f107, f107a, msis_dens10, msis_temp1) diff --git a/src/set_inputs.f90 b/src/set_inputs.f90 index 45d3e1a6..4bcc49db 100644 --- a/src/set_inputs.f90 +++ b/src/set_inputs.f90 @@ -374,8 +374,8 @@ subroutine set_inputs write(*, *) '#MSISOBC' write(*, *) 'UseOBCExperiment (logical)' write(*, *) 'MsisOblateFactor (real)' - endif - + endif + case ("#MSIS21") call read_in_logical(UseMsis21, iError) if (iError /= 0) then diff --git a/src/set_vertical_bcs.Earth.f90 b/src/set_vertical_bcs.Earth.f90 index 6990e9f7..539dbe6d 100644 --- a/src/set_vertical_bcs.Earth.f90 +++ b/src/set_vertical_bcs.Earth.f90 @@ -139,18 +139,18 @@ subroutine set_vertical_bcs(LogRho, LogNS, Vel_GD, Temp, LogINS, iVel, VertVel) Vel_GD(-1:0, iEast_) = 0.0 Vel_GD(-1:0, iNorth_) = 0.0 ! The rest of the BCs will just stay constant. - endif + endif - if (DoN4SHack) then + if (DoN4SHack) then !! Hack for iN_4S_ for lower altitudes: do iAlt = nAlts, 0, -1 - if (NS(iAlt - 1, iN_4S_) < 1000.0) then - NS(iAlt - 1, iN_4S_) = NS(iAlt, iN_4S_) - logNS(iAlt - 1, iN_4S_) = log(NS(iAlt - 1, iN_4S_)) - endif + if (NS(iAlt - 1, iN_4S_) < 1000.0) then + NS(iAlt - 1, iN_4S_) = NS(iAlt, iN_4S_) + logNS(iAlt - 1, iN_4S_) = log(NS(iAlt - 1, iN_4S_)) + endif enddo - endif - + endif + if (.not. DuringPerturb) then Vel_GD(-1:0, iUp_) = 0.0 VertVel(-1:0, :) = 0.0 diff --git a/srcPython/pGITM.py b/srcPython/pGITM.py index d90e8c05..236d1f84 100755 --- a/srcPython/pGITM.py +++ b/srcPython/pGITM.py @@ -1,10 +1,10 @@ #!/usr/bin/env python3 ''' -A wrapper for PostGITM.exe to handle all unconcatenated output fragments +Handle all unconcatenated output fragments produced during a GITM simulation. -This script is an alternative to pGITM (written in cshell) that works robustly -and securely on a range of unix-like environments with Python available +This script works robustly and securely on a range of +unix-like environments with Python available (i.e., any modern computer.) ''' @@ -28,7 +28,7 @@ def get_args_pgitm(): parser = argparse.ArgumentParser( - description = 'Plot Aether / GITM model results') + description = 'Post-process GITM model results') parser.add_argument('-v', \ action='store_true', default = False, \ help = 'set verbose to true') @@ -374,14 +374,14 @@ def post_process_gitm(dir, doRemove, isVerbose = False): print('Processing directory doesnt exist: ', processDir) print('Make sure to put in a valid -dir= ') exit() - # Move into the UA/data directory (required by PostGITM.exe) + # Move into the UA/data directory if (isVerbose): print('Moving into processing directory : ', processDir) cwd = os.getcwd() os.chdir(processDir) # Get list of header files to process: - headers = glob('*.header') + headers = sorted(glob('*.header')) if (isVerbose): print('Found %i files to process' % len(headers)) diff --git a/srcPython/post_process.py b/srcPython/post_process.py index be56f573..1a0deb72 100755 --- a/srcPython/post_process.py +++ b/srcPython/post_process.py @@ -16,8 +16,11 @@ def parse_args_post(): - parser = argparse.ArgumentParser(description = - 'Post process and move model results') + parser = argparse.ArgumentParser( + description = "Post process and (optionally) move model results.\n"+ + "- This functions similar to pGITM.py, but can copy files to\n"+ + " a remote location, or postprocess files as they are created.", + formatter_class=argparse.RawTextHelpFormatter) parser.add_argument('-remotefile', help = 'File that contains info. about remote system', default = 'remote') @@ -31,15 +34,15 @@ def parse_args_post(): default = 'none') parser.add_argument('-dir', - help = 'remote directory to use', + help = 'remote directory to use (default none)', default = 'none') parser.add_argument('-sleep', - help = 'how long to sleep between loops', + help = 'how long to sleep between loops in seconds, (default 300)', default = 300, type = int) parser.add_argument('-totaltime', - help = 'specify how long to run in total (in hours)', + help = 'specify how long to run in total in hours, (default 0 - only run once)', default = 0, type = int) parser.add_argument('-q', @@ -426,7 +429,14 @@ def do_loop(doTarZip, user, server, dir, IsRemote): if (doTarZip): DidWork = tar_and_zip_gitm() else: + # Default should be UA/data processDir = 'UA/data' + if (not os.path.exists(processDir)): + # Maybe we are in the UA directory? + processDir = 'data' + if (not os.path.exists(processDir)): + # Maybe we are already in the data directory??? + processDir = '.' DidWork = post_process_gitm(processDir, DoRm, isVerbose = IsVerbose) #DidWork = post_process_gitm() @@ -508,7 +518,11 @@ def load_remote_file(args): currentTime = datetime.now() dt = ((currentTime - startTime).total_seconds())/3600.0 if (dt > args.totaltime): - print(" --> Stopping due to totaltime exceeded!") + if args.totaltime == 0: + # Different exit message for non-continuous runs + print(" --> All done!") + else: + print(" --> Stopping due to totaltime exceeded!") # want to break out of loop, so set loop breaking condition: DidWork = False