From d172ee785a076fefedee5fff8679cf1abc11d7f9 Mon Sep 17 00:00:00 2001 From: Aaron James Ridley Date: Thu, 8 May 2025 08:32:29 -0400 Subject: [PATCH 01/12] BUG: should sort file list to process oldest files first --- srcPython/pGITM.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/srcPython/pGITM.py b/srcPython/pGITM.py index d90e8c05..48e0e592 100755 --- a/srcPython/pGITM.py +++ b/srcPython/pGITM.py @@ -381,7 +381,7 @@ def post_process_gitm(dir, doRemove, isVerbose = False): 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)) From e5ad5e2fafe0c8246b65bf4ee418e3846fa5ecc1 Mon Sep 17 00:00:00 2001 From: Aaron James Ridley Date: Thu, 8 May 2025 16:29:17 -0400 Subject: [PATCH 02/12] FEAT: work from UA and UA/data directories too! --- srcPython/post_process.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/srcPython/post_process.py b/srcPython/post_process.py index be56f573..1247c8df 100755 --- a/srcPython/post_process.py +++ b/srcPython/post_process.py @@ -426,7 +426,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() From fdf1f5ef82337ee4a9cf03fe93a30c38fdf71c16 Mon Sep 17 00:00:00 2001 From: Aaron Bukowski Date: Tue, 13 May 2025 13:34:58 +0000 Subject: [PATCH 03/12] STY: Different exit messages when not looping (in post_process.py) --- srcPython/post_process.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/srcPython/post_process.py b/srcPython/post_process.py index 1247c8df..4297baa3 100755 --- a/srcPython/post_process.py +++ b/srcPython/post_process.py @@ -515,7 +515,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 From 0828753c32a5a0cbc22caf4e4fba71ff48802866 Mon Sep 17 00:00:00 2001 From: Aaron Bukowski Date: Wed, 14 May 2025 15:38:48 +0000 Subject: [PATCH 04/12] MAINT: Don't build legacy postprocessors by default (can still be made with `make POST`) --- Makefile | 1 - src/Makefile | 2 -- 2 files changed, 3 deletions(-) diff --git a/Makefile b/Makefile index 8b66be20..de673704 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 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} From 224013caf30725c4dd2c47ff4b11c96b06ddfedd Mon Sep 17 00:00:00 2001 From: Aaron Bukowski Date: Wed, 14 May 2025 15:43:12 +0000 Subject: [PATCH 05/12] FEAT: Change default postprocessor! details: - do not copy postprocess.exe (pGITM) - do not copy pGITM.py - yes copy post_process.py - update README with new postprocess info --- Makefile | 4 +--- README.md | 21 ++++++++++----------- 2 files changed, 11 insertions(+), 14 deletions(-) diff --git a/Makefile b/Makefile index de673704..6ae6b49f 100644 --- a/Makefile +++ b/Makefile @@ -108,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 From ec4f24ae0448c157efec4eefa1437086097a191b Mon Sep 17 00:00:00 2001 From: Aaron Bukowski Date: Wed, 14 May 2025 15:43:55 +0000 Subject: [PATCH 06/12] doc: edit comments in pGITM.py - it doesn't use postGITM.exe --- srcPython/pGITM.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/srcPython/pGITM.py b/srcPython/pGITM.py index 48e0e592..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,7 +374,7 @@ 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() From 1bb9854ea48d28715fe1c9cfbe5e9b8b6871e360 Mon Sep 17 00:00:00 2001 From: Aaron Bukowski Date: Wed, 14 May 2025 15:44:20 +0000 Subject: [PATCH 07/12] doc: update docstrings & help messages in post_process.py --- srcPython/post_process.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/srcPython/post_process.py b/srcPython/post_process.py index 4297baa3..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', From eadcc7f2c48aa042d992751aa0be474479a0c91c Mon Sep 17 00:00:00 2001 From: Aaron Bukowski Date: Wed, 14 May 2025 15:57:22 +0000 Subject: [PATCH 08/12] sty: prettifying --- src/ModInputs.f90 | 2 +- src/init_altitude.f90 | 4 +- src/init_msis.Earth.f90 | 137 ++++++++++++++++----------------- src/set_inputs.f90 | 4 +- src/set_vertical_bcs.Earth.f90 | 16 ++-- 5 files changed, 81 insertions(+), 82 deletions(-) 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 From 85fd88e6f348c42ba76e06a51b93170742af71bb Mon Sep 17 00:00:00 2001 From: Aaron Bukowski Date: Thu, 15 May 2025 21:16:49 +0000 Subject: [PATCH 09/12] BUG: Formatter issues w/ MSIS --- .fprettify.rc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 From 56c84b66b3bc0aeb32ebdbbfc8cc49266cacb135 Mon Sep 17 00:00:00 2001 From: Aaron Bukowski Date: Fri, 16 May 2025 10:11:30 -0400 Subject: [PATCH 10/12] Revert "sty: prettifying" UNprettifying, oy This reverts commit eadcc7f2c48aa042d992751aa0be474479a0c91c. --- src/ModInputs.f90 | 2 +- src/init_altitude.f90 | 4 +- src/init_msis.Earth.f90 | 137 +++++++++++++++++---------------- src/set_inputs.f90 | 4 +- src/set_vertical_bcs.Earth.f90 | 16 ++-- 5 files changed, 82 insertions(+), 81 deletions(-) diff --git a/src/ModInputs.f90 b/src/ModInputs.f90 index 44b40c71..17a586f4 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 96b1abc5..cf372382 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 c650cf13..5fed7104 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,6 +554,7 @@ 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 4bcc49db..45d3e1a6 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 539dbe6d..6990e9f7 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 From 56c235c20f3036b98c900bdd1de9ddfa17be7c6d Mon Sep 17 00:00:00 2001 From: Aaron Bukowski Date: Fri, 16 May 2025 10:59:34 -0400 Subject: [PATCH 11/12] bug: force format test to use up-to-date PR --- .github/workflows/Format-Doc-Test.yml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) 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 From 1361c33804f35aee5077f5487f117fa400e3a6b2 Mon Sep 17 00:00:00 2001 From: Aaron Bukowski Date: Fri, 16 May 2025 11:05:26 -0400 Subject: [PATCH 12/12] sty: RE prettifying, yo Formatter wasn't grabbing the most code that it should have been. --- src/ModInputs.f90 | 2 +- src/init_altitude.f90 | 4 +- src/init_msis.Earth.f90 | 137 ++++++++++++++++----------------- src/set_inputs.f90 | 4 +- src/set_vertical_bcs.Earth.f90 | 16 ++-- 5 files changed, 81 insertions(+), 82 deletions(-) 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