diff --git a/kaipy/gamera/msphViz.py b/kaipy/gamera/msphViz.py index 26a93377..4b4c2705 100644 --- a/kaipy/gamera/msphViz.py +++ b/kaipy/gamera/msphViz.py @@ -23,11 +23,12 @@ eMax = 5.0 #Max current for contours #Default pressure colorbar vP = kv.genNorm(vMin=1.0e-2,vMax=10.0,doLog=True) -szStrs = ['small','std','big','dm'] +szStrs = ['small','std','big','bigger','fullD','dm'] szBds = {} szBds["std"] = [-40.0 ,20.0,2.0] szBds["big"] = [-100.0,20.0,2.0] szBds["bigger"] = [-200.0,25.0,2.0] +szBds["fullD"] = [-300.0,30.0,3.0] # full domain for double res szBds["small"] = [-10.0 , 5.0,2.0] szBds["dm"] = [-30.0 ,10.0,40.0/15.0] @@ -48,6 +49,209 @@ def GetSizeBds(args): return xyBds +#Plot absolute error in the requested, or given, equatorial field +def PlotEqErrAbs(gsphP,gsphO,nStp,xyBds,Ax,fieldNames,AxCB=None,doClear=True,doDeco=True,vMin=1e-9,vMax=1e-4,doLog=True,doVerb=True): + #specify two gsph objects as gsphPredicted (gsphP) and gsphObserved (gsphO) + # + normAbs = kv.genNorm(vMin,vMax,doLog=doLog) + cmapAbs = "PRGn" + if doLog: + cmapAbs = "magma" + + if (AxCB is not None): + #Add the colorbar to AxCB + AxCB.clear() + + cbString = "" + for fn in fieldNames: + cbString = cbString + "'" + fn + "', " + cbString = cbString + " Absolute Error" + kv.genCB(AxCB,normAbs,cbString,cM=cmapAbs) + #Now do main plotting + if (doClear): + Ax.clear() + dataAbs = None + for fn in fieldNames: + dataP = gsphP.EggSlice(fn,nStp,doEq=True,doVerb=doVerb) + dataO = gsphO.EggSlice(fn,nStp,doEq=True,doVerb=doVerb) + if dataAbs is None: + dataAbs = np.square(dataO - dataP) + else: + dataAbs = dataAbs + np.square(dataO - dataP) + dataAbs = np.sqrt(dataAbs) + Ax.pcolormesh(gsphP.xxi,gsphP.yyi,dataAbs,cmap=cmapAbs,norm=normAbs) + + kv.SetAx(xyBds,Ax) + + if (doDeco): + kv.addEarth2D(ax=Ax) + Ax.set_xlabel('SM-X [Re]') + Ax.set_ylabel('SM-Y [Re]') + return dataAbs + +#Plot relative error in the requested, or given, equatorial field +def PlotEqErrRel(gsphP,gsphO,nStp,xyBds,Ax,fieldNames,AxCB=None,doClear=True,doDeco=True,vMin=1e-16,vMax=1,doLog=True,doVerb=True): + #specify two gsph objects as gsphPredicted (gsphP) and gsphObserved (gsphO) + # + normRel = kv.genNorm(vMin,vMax,doLog=doLog) + cmapRel = "PRGn" + if doLog: + cmapRel = "magma" + + if (AxCB is not None): + #Add the colorbar to AxCB + AxCB.clear() + cbString = "" + for fn in fieldNames: + cbString = cbString + "'" + fn + "', " + cbString = cbString + " Relative Error" + kv.genCB(AxCB,normRel,cbString,cM=cmapRel) + #Now do main plotting + if (doClear): + Ax.clear() + dataAbs = None + dataBase = None + for fn in fieldNames: + dataP = gsphP.EggSlice(fn,nStp,doEq=True,doVerb=doVerb) + dataO = gsphO.EggSlice(fn,nStp,doEq=True,doVerb=doVerb) + if dataAbs is None: + dataBase = np.square(dataP) + dataAbs = np.square(dataO - dataP) + else: + dataBase = dataBase + np.square(dataP) + dataAbs = dataAbs + np.square(dataO - dataP) + dataBase = np.sqrt(dataBase) + dataAbs = np.sqrt(dataAbs) + # math breaks when Base field is exactly 0 + dataBase[dataBase == 0] = np.finfo(np.float32).tiny + dataRel = np.absolute(dataAbs/dataBase) + Ax.pcolormesh(gsphP.xxi,gsphP.yyi,dataRel,cmap=cmapRel,norm=normRel) + + kv.SetAx(xyBds,Ax) + + if (doDeco): + kv.addEarth2D(ax=Ax) + Ax.set_xlabel('SM-X [Re]') + Ax.set_ylabel('SM-Y [Re]') + return dataRel + +#Plot absolute error along the requested logical axis +def PlotLogicalErrAbs(gsphP,gsphO,nStp,Ax,fieldNames,meanAxis,AxCB=None,doClear=True,doDeco=True,vMin=1e-16,vMax=1,doLog=True,doVerb=True): + #specify two gsph objects as gsphPredicted (gsphP) and gsphObserved (gsphO) + # + normAbs = kv.genNorm(vMin,vMax,doLog=doLog) + cmapAbs = "PRGn" + if doLog: + cmapAbs = "magma" + + if (AxCB is not None): + #Add the colorbar to AxCB + AxCB.clear() + cbString = "" + for fn in fieldNames: + cbString = cbString + "'" + fn + "', " + cbString = cbString + " Absolute Error along " + if meanAxis == 0: + cbString = cbString + "I axis" + elif meanAxis == 1: + cbString = cbString + "J axis" + elif meanAxis == 2: + cbString = cbString + "K axis" + kv.genCB(AxCB,normAbs,cbString,cM=cmapAbs) + #Now do main plotting + if (doClear): + Ax.clear() + dataAbs = CalcTotalErrAbs(gsphP,gsphO,nStp,fieldNames,doVerb=doVerb,meanAxis=meanAxis) + dataAbs = np.transpose(dataAbs) # transpose to put I/J on the horizontal axis + Ax.pcolormesh(dataAbs,cmap=cmapAbs,norm=normAbs) + + if (doDeco): + if meanAxis == 0: + Ax.set_xlabel('J Indices') + Ax.set_ylabel('K Indices') + elif meanAxis == 1: + Ax.set_xlabel('I Indices') + Ax.set_ylabel('K Indices') + elif meanAxis == 2: + Ax.set_xlabel('I Indices') + Ax.set_ylabel('J Indices') + return dataAbs + +#Plot relative error along the requested logical axis +def PlotLogicalErrRel(gsphP,gsphO,nStp,Ax,fieldNames,meanAxis,AxCB=None,doClear=True,doDeco=True,vMin=1e-16,vMax=1,doLog=True,doVerb=True): + #specify two gsph objects as gsphPredicted (gsphP) and gsphObserved (gsphO) + # + normRel = kv.genNorm(vMin,vMax,doLog=doLog) + cmapRel = "PRGn" + if doLog: + cmapRel = "magma" + + if (AxCB is not None): + #Add the colorbar to AxCB + AxCB.clear() + cbString = "" + for fn in fieldNames: + cbString = cbString + "'" + fn + "', " + cbString = cbString + " Relative Error along " + if meanAxis == 0: + cbString = cbString + "I axis" + elif meanAxis == 1: + cbString = cbString + "J axis" + elif meanAxis == 2: + cbString = cbString + "K axis" + kv.genCB(AxCB,normAbs,cbString,cM=cmapAbs) + #Now do main plotting + if (doClear): + Ax.clear() + dataRel = CalcTotalErrRel(gsphP,gsphO,nStp,fieldNames,doVerb=doVerb,meanAxis=meanAxis) + dataRel = np.transpose(dataRel) # transpose to put I/J on the horizontal axis + Ax.pcolormesh(dataRel,cmap=cmapRel,norm=normRel) + + if (doDeco): + if meanAxis == 0: + Ax.set_xlabel('J Indices') + Ax.set_ylabel('K Indices') + elif meanAxis == 1: + Ax.set_xlabel('I Indices') + Ax.set_ylabel('K Indices') + elif meanAxis == 2: + Ax.set_xlabel('I Indices') + Ax.set_ylabel('J Indices') + return dataRel + +#Calculate total cumulative absolute error between two cases +def CalcTotalErrAbs(gsphP,gsphO,nStp,fieldNames,doVerb=True,meanAxis=None): + dataAbs = None + for fn in fieldNames: + dataP = gsphP.GetVar(fn,nStp,doVerb=doVerb) + dataO = gsphO.GetVar(fn,nStp,doVerb=doVerb) + if dataAbs is None: + dataAbs = np.square(dataO - dataP) + else: + dataAbs = dataAbs + np.square(dataO - dataP) + dataAbs = np.sqrt(dataAbs) + return np.mean(dataAbs,axis=meanAxis) + +#Calculate total cumulative relative error between two cases +def CalcTotalErrRel(gsphP,gsphO,nStp,fieldNames,doVerb=True,meanAxis=None): + dataAbs = None + dataBase = None + for fn in fieldNames: + dataP = gsphP.GetVar(fn,nStp,doVerb=doVerb) + dataO = gsphO.GetVar(fn,nStp,doVerb=doVerb) + if dataAbs is None: + dataBase = np.square(dataP) + dataAbs = np.square(dataO - dataP) + else: + dataBase = dataBase + np.square(dataP) + dataAbs = dataAbs + np.square(dataO - dataP) + dataBase = np.sqrt(dataBase) + dataAbs = np.sqrt(dataAbs) + # math breaks when Base field is exactly 0 + dataBase[dataBase == 0] = np.finfo(np.float32).tiny + dataRel = np.absolute(dataAbs/dataBase) + return np.mean(dataRel,axis=meanAxis) + #Plot equatorial field def PlotEqB(gsph,nStp,xyBds,Ax,AxCB=None,doClear=True,doDeco=True,doBz=False): vBZ = kv.genNorm(dbMax) diff --git a/scripts/quicklook/gamerrVid.py b/scripts/quicklook/gamerrVid.py new file mode 100755 index 00000000..b2b082bb --- /dev/null +++ b/scripts/quicklook/gamerrVid.py @@ -0,0 +1,211 @@ +#!/usr/bin/env python +#Make video of error between two Gamera cases +import argparse +from argparse import RawTextHelpFormatter +import matplotlib as mpl +mpl.use('Agg') +import matplotlib.pyplot as plt +import kaipy.kaiViz as kv +import matplotlib.gridspec as gridspec +import numpy as np +import numpy as np +import kaipy.gamera.msphViz as mviz +import kaipy.gamera.magsphere as msph +import kaipy.gamera.rcmpp as rcmpp +from alive_progress import alive_bar +import kaipy.kdefs as kdefs +from mpl_toolkits.axes_grid1.inset_locator import inset_axes +import os +import errno + +cLW = 0.25 + +if __name__ == "__main__": + #Defaults + fdir1 = os.getcwd() + ftag1 = "msphere" + fdir2 = os.getcwd() + ftag2 = "msphere" + oDir = "vid2D" + ts = 0 #[min] + te = 200 #[min] + dt = 60.0 #[sec] + Nblk = 1 #Number of blocks + nID = 1 #Block ID of this job + noMPI = False # Don't add MPI tiling + noLog = False + fieldNames = "Bx, By, Bz" + doVerb = False + + MainS = """Creates simple multi-panel figure for Gamera magnetosphere run + Left Panel - Residual vertical magnetic field + Right Panel - Pressure (or density) and hemispherical insets + """ + + parser = argparse.ArgumentParser(description=MainS, formatter_class=RawTextHelpFormatter) + parser.add_argument('-d1',type=str,metavar="directory",default=fdir1,help="Directory to read first dataset from (default: %(default)s)") + parser.add_argument('-id1',type=str,metavar="runid",default=ftag1,help="RunID of first dataset (default: %(default)s)") + parser.add_argument('-d2',type=str,metavar="directory",default=fdir2,help="Directory to read second dataset from (default: %(default)s)") + parser.add_argument('-id2',type=str,metavar="runid",default=ftag2,help="RunID of second dataset (default: %(default)s)") + parser.add_argument('-o',type=str,metavar="directory",default=oDir,help="Subdirectory to write to (default: %(default)s)") + parser.add_argument('-ts' ,type=int,metavar="tStart",default=ts,help="Starting time [min] (default: %(default)s)") + parser.add_argument('-te' ,type=int,metavar="tEnd" ,default=te,help="Ending time [min] (default: %(default)s)") + parser.add_argument('-dt' ,type=int,metavar="dt" ,default=dt,help="Cadence [sec] (default: %(default)s)") + parser.add_argument('-Nblk' ,type=int,metavar="Nblk",default=Nblk,help="Number of job blocks (default: %(default)s)") + parser.add_argument('-nID' ,type=int,metavar="nID" ,default=nID,help="Block ID of this job [1-Nblk] (default: %(default)s)") + parser.add_argument('-f',type=str,metavar="fieldnames",default=fieldNames,help="Comma-separated fields to plot (default: %(default)s)") + parser.add_argument('-linear',action='store_true', default=noLog,help="Plot linear line plot instead of logarithmic (default: %(default)s)") + parser.add_argument('-v',action='store_true', default=doVerb,help="Do verbose output (default: %(default)s)") + #parser.add_argument('-nompi', action='store_true', default=noMPI,help="Don't show MPI boundaries (default: %(default)s)") + + mviz.AddSizeArgs(parser) + + #Finalize parsing + args = parser.parse_args() + fdir1 = args.d1 + ftag1 = args.id1 + fdir2 = args.d2 + ftag2 = args.id2 + ts = args.ts + te = args.te + dt = args.dt + oSub = args.o + Nblk = args.Nblk + nID = args.nID + fieldNames = args.f + noLog = args.linear + doVerb = args.v + #noMPI = args.noMPI + + fnList = [item.strip() for item in fieldNames.split(',')] + + #Setup timing info + tOut = np.arange(ts*60.0,te*60.0,dt) + Nt = len(tOut) + vO = np.arange(0,Nt) + + print("Writing %d outputs between minutes %d and %d"%(Nt,ts,te)) + if (Nblk>1): + #Figure out work bounds + dI = (Nt//Nblk) + i0 = (nID-1)*dI + i1 = i0+dI + if (nID == Nblk): + i1 = Nt #Make sure we get last bit + print("\tBlock #%d: %d to %d"%(nID,i0,i1)) + else: + i0 = 0 + i1 = Nt + + #Setup output directory + oDir = os.getcwd() + "/" + oSub + print("Writing output to %s"%(oDir)) + + #Check/create directory if necessary + if (not os.path.exists(oDir)): + try: + print("Creating directory %s"%(oDir)) + os.makedirs(oDir) + except OSError as exc: + if exc.errno == errno.EEXIST and os.path.isdir(oDir): + pass + else: + raise + + #Get domain size + xyBds = mviz.GetSizeBds(args) + + #--------- + #Figure parameters + figSz = (12,7.5) + + + #====== + #Init data + gsph1 = msph.GamsphPipe(fdir1,ftag1) + gsph2 = msph.GamsphPipe(fdir2,ftag2) + + #====== + #Setup figure + fig = plt.figure(figsize=figSz) + gs = gridspec.GridSpec(5,2,height_ratios=[20,5,1,5,9],hspace=0.025) + + + AxTL = fig.add_subplot(gs[0,0]) + AxTR = fig.add_subplot(gs[0,1]) + AxB = fig.add_subplot(gs[-1,0:2]) + AxB2 = AxB.twinx() # second plot on bottom axis + + AxCT = fig.add_subplot(gs[2,0:2]) + + print(fig.axes) + + errTimes = [] + errListRel = [] + errListAbs = [] + relColor = "tab:blue" + absColor = "tab:orange" + + #Loop over sub-range + titstr = "Comparing '%s' to '%s'"%(fdir1,fdir2) + with alive_bar(i1-i0,title=titstr.ljust(kdefs.barLab),length=kdefs.barLen,disable=doVerb) as bar: + for i in range(i0,i1): + #Convert time (in seconds) to Step # + nStp = np.abs(gsph1.T-tOut[i]).argmin()+gsph1.s0 + if doVerb: + print("Minute = %5.2f / Step = %d"%(tOut[i]/60.0,nStp)) + npl = vO[i] + + AxTL.clear() + AxTR.clear() + AxB.clear() + AxB2.clear() + + #plot upper left msph error + mviz.PlotEqErrRel(gsph1,gsph2,nStp,xyBds,AxTL,fnList,AxCB=AxCT,doVerb=doVerb) + AxTL.set_title("Equatorial Slice of Relative Error") + + #plot upper right k-axis error + mviz.PlotLogicalErrRel(gsph1,gsph2,nStp,AxTR,fnList,2,doVerb=doVerb) + AxTR.set_title("Per-Cell Relative Error along K-Axis") + if (not noMPI): + #plot I-MPI decomp on logical plot + if(gsph2.Ri > 1): + for im in range(gsph2.Ri): + i0 = im*gsph2.dNi + AxTR.plot([i0, i0],[0, gsph2.Nj],"deepskyblue",linewidth=0.25,alpha=0.5) + #plot J-MPI decomp on logical plot + if (gsph2.Rj>1): + for jm in range(1,gsph2.Rj): + j0 = jm*gsph2.dNj + AxTR.plot([0, gsph2.Ni],[j0, j0],"deepskyblue",linewidth=0.25,alpha=0.5) + + #plot bottom line plot + errTimes.append(tOut[i]/60.0) + errListRel.append(mviz.CalcTotalErrRel(gsph1,gsph2,nStp,fnList,doVerb=doVerb)) + errListAbs.append(mviz.CalcTotalErrAbs(gsph1,gsph2,nStp,fnList,doVerb=doVerb)) + if noLog: + AxB.plot(errTimes, errListRel,color=relColor) + AxB2.plot(errTimes, errListAbs,color=absColor) + else: + AxB.semilogy(errTimes, errListRel,color=relColor) + AxB2.semilogy(errTimes, errListAbs,color=absColor) + AxB.set_xlabel('Time (min)') + AxB.set_ylabel('Per-Cell Mean Relative Error',color=relColor) + AxB.tick_params(axis='y',which='both',colors=relColor,left=True,right=True,labelleft=True,labelright=False) + AxB2.set_ylabel('Per-Cell Mean Absolute Error',color=absColor) + #AxB2.yaxis.tick_right() + AxB2.tick_params(axis='y',which='both',colors=absColor,left=True,right=True,labelleft=False,labelright=True) + AxB.set_title("'" + fieldNames + "' Per-Cell Error Over Time") + + gsph1.AddTime(nStp,AxTL,xy=[0.025,0.84],fs="x-large") + + #Add MPI decomp + if (not noMPI): + mviz.PlotMPI(gsph2,AxTL) + + fOut = oDir+"/vid.%04d.png"%(npl) + kv.savePic(fOut,bLenX=45) + + bar() + diff --git a/scripts/quicklook/gamerrpic.py b/scripts/quicklook/gamerrpic.py new file mode 100755 index 00000000..b00ff361 --- /dev/null +++ b/scripts/quicklook/gamerrpic.py @@ -0,0 +1,99 @@ +#!/usr/bin/env python +#Make video of error between two Gamera cases +import argparse +from argparse import RawTextHelpFormatter +import matplotlib as mpl +mpl.use('Agg') +import matplotlib.pyplot as plt +import kaipy.kaiViz as kv +import matplotlib.gridspec as gridspec +import numpy as np +import numpy as np +import kaipy.gamera.msphViz as mviz +import kaipy.gamera.magsphere as msph +import kaipy.gamera.rcmpp as rcmpp +from mpl_toolkits.axes_grid1.inset_locator import inset_axes +import os +import errno + +cLW = 0.25 + +if __name__ == "__main__": + #Defaults + fdir1 = os.getcwd() + ftag1 = "msphere" + fdir2 = os.getcwd() + ftag2 = "msphere" + nStp=1 + fieldNames = "Bx, By, Bz" + doMPI = False #[Add MPI tiling] + noMPI = False + + MainS = """Creates simple multi-panel figure for Gamera magnetosphere run + Left Panel - Residual vertical magnetic field + Right Panel - Pressure (or density) and hemispherical insets + """ + + parser = argparse.ArgumentParser(description=MainS, formatter_class=RawTextHelpFormatter) + parser.add_argument('-d1',type=str,metavar="directory1",default=fdir1,help="Directory to read first dataset from (default: %(default)s)") + parser.add_argument('-id1',type=str,metavar="runid1",default=ftag1,help="RunID of first dataset (default: %(default)s)") + parser.add_argument('-d2',type=str,metavar="directory2",default=fdir2,help="Directory to read second dataset from (default: %(default)s)") + parser.add_argument('-id2',type=str,metavar="runid2",default=ftag2,help="RunID of second dataset (default: %(default)s)") + parser.add_argument('-n',type=int,metavar="nStp",default=nStp,help="Step number to plot (default: %(default)s)") + parser.add_argument('-f',type=str,metavar="fieldnames",default=fieldNames,help="Comma-separated fields to plot (default: %(default)s)") + #parser.add_argument('-nompi', action='store_true', default=noMPI,help="Don't show MPI boundaries (default: %(default)s)") + + mviz.AddSizeArgs(parser) + + #Finalize parsing + args = parser.parse_args() + fdir1 = args.d1 + ftag1 = args.id1 + fdir2 = args.d2 + ftag2 = args.id2 + nStp = args.n + fieldNames = args.f + oName = "gamErrPic.png" + + #Get domain size + xyBds = mviz.GetSizeBds(args) + + #--------- + #Figure parameters + figSz = (12,7.5) + + + #====== + #Init data + gsph1 = msph.GamsphPipe(fdir1,ftag1) + gsph2 = msph.GamsphPipe(fdir2,ftag2) + + #====== + #Setup figure + fig = plt.figure(figsize=figSz) + gs = gridspec.GridSpec(2,2,height_ratios=[20,1],hspace=0.025) + + + AxL = fig.add_subplot(gs[0,0]) + AxR = fig.add_subplot(gs[0,1]) + + AxCL = fig.add_subplot(gs[-1,0]) + AxCR = fig.add_subplot(gs[-1,1]) + + fnList = [item.strip() for item in fieldNames.split(',')] + + AxL.clear() + AxR.clear() + + mviz.PlotEqErrRel(gsph1,gsph2,nStp,xyBds,AxL,fnList,AxCB=AxCL) + mviz.PlotEqErrAbs(gsph1,gsph2,nStp,xyBds,AxR,fnList,AxCB=AxCR) + + gsph1.AddTime(nStp,AxL,xy=[0.025,0.89],fs="x-large") + + #Add MPI decomp + if (doMPI): + mviz.PlotMPI(gsph2,AxL) + mviz.PlotMPI(gsph2,AxR) + + kv.savePic(oName,bLenX=45) + diff --git a/scripts/quicklook/swpic.py b/scripts/quicklook/swpic.py index 913f270a..572af55b 100755 --- a/scripts/quicklook/swpic.py +++ b/scripts/quicklook/swpic.py @@ -48,15 +48,28 @@ if __name__ == "__main__": fOut = swtag.split('.')[0]+'.'+imgtype # pulling UT variable for plotting - t0Fmt = "%Y-%m-%d %H:%M:%S" + t0Fmts = ["%Y-%m-%d %H:%M:%S","%Y-%m-%dT%H:%M:%S.%f"] utfmt='%H:%M \n%Y-%m-%d' UTall = kh5.PullVar(swIn,"UT") + #Identify the correct time format + t0Fmt = None + for tfmt in t0Fmts: + try: + datetime.datetime.strptime(UTall[1].decode('utf-8'),tfmt) + t0Fmt = tfmt + break # datetime parse succeeded + except ValueError: + pass # datetime parse failed + if t0Fmt is None: + print("Time format in bcwind.h5 did not match any expected format.") + sys.exit() + utall = [] for n in range(len(UTall)): utall.append(datetime.datetime.strptime(UTall[n].decode('utf-8'),t0Fmt)) - + # pulling the solar wind values from the table varlist = kh5.getRootVars(swIn) D = kh5.PullVar(swIn,"D") @@ -74,6 +87,6 @@ if __name__ == "__main__": else: pltInterp = 0*D doEps = False - swBCplots.swQuickPlot(UTall,D,Temp,Vx,Vy,Vz,Bx,By,Bz,SYMH,pltInterp,fOut,doEps=doEps,doTrim=doTrim) + swBCplots.swQuickPlot(UTall,D,Temp,Vx,Vy,Vz,Bx,By,Bz,SYMH,pltInterp,fOut,doEps=doEps,doTrim=doTrim,t0fmt=t0Fmt) diff --git a/src/chimp/streamline.F90 b/src/chimp/streamline.F90 index c0c1b251..60e0046a 100644 --- a/src/chimp/streamline.F90 +++ b/src/chimp/streamline.F90 @@ -706,7 +706,9 @@ module streamline if (inDom) Np = Np + 1 enddo - if (MaxFL == Np) then + ! check if exceeded tube bounds + if(Np > MaxFL) then + Np = MaxFL !$OMP CRITICAL write(*,*) ANSIRED write(*,*) ""