diff --git a/bin/calc_montg1.py b/bin/calc_montg1.py index f93f0c8d..1096666a 100755 --- a/bin/calc_montg1.py +++ b/bin/calc_montg1.py @@ -1,17 +1,17 @@ #!/usr/bin/env python -import argparse -import matplotlib -matplotlib.use('Agg') -import matplotlib.pyplot -import abfile.abfile as abf -import numpy import modeltools.hycom -import logging -import re -import os +import subprocess +import argparse +import datetime import os.path - +import sys +import re +import logging +import abfile +import pickle +import numpy # Set up logger +##_loglevel=logging.DEBUG _loglevel=logging.INFO logger = logging.getLogger(__name__) logger.setLevel(_loglevel) @@ -22,23 +22,21 @@ logger.addHandler(ch) logger.propagate=False -# Note that this part can be easily implemented into the python code -# that aims to generate nesting [ab] files. I may not have time to manage it -# but evrything has been provided here or other earlier developed NERSC scripts. -# By Mostafa Bakhoday-Paskyabi -# July 4th 2019 (almost at last day of my work at NERSC) -# -# ../bin/calc_montg1.py ../nest/015/archv.2016_001_00.a ./ -# ../bin/calc_montg1.py YOUR_INPUT_AFILE OUTPUT_PATH -# It saves with same name as a-file to outpput path +""" +usage: +python ../bin/topaz6_calc_montg1.py archv_file restart_file output_path +python ../bin/topaz6_calc_montg1.py ../nest/070/archv.2022_004_00.a ./data/TP6restart.2021_365_00_0000.a ./ -cline1="HYCOM archive files from climatology generated from" -cline2="NEMO Relaxation fields on native grid, MBP2019" -cline3="Expt 01.5 nhybrd=50 nsigma= 0 ds00= 1.00 dp00= 1.00 dp00x= 450.0 dp00f=1.150" +python /cluster/home/achoth/pyScripts_AO_Betzy/ValidationFreq/Alf_calc_montg1.py /cluster/work/users/achoth/TP5a0.06/nest/080_NewMontg/archv.1997_001_00.b /cluster/work/users/achoth/TP5a0.06/expt_08.1/data/restart.1997_001_00_0000.b /cluster/work/users/achoth/TP5a0.06/nest/test_dir/ + +""" +### AO +# Get component of montgomery potential at surface that depends on pbavg def montg1_pb(thstar,p) : kdm=thstar.shape[0] + #AO thref=_constants.thref thref=1e-3 # Part that depends on pbavg : montgpb=-thstar[kdm-1,:,:]*thref**2 @@ -49,149 +47,179 @@ def montg1_pb(thstar,p) : - # Get component of montgomery potential at surface that does not depend on pbavg def montg1_no_pb(psikk,thkk,thstar,p) : kdm=thstar.shape[0] + ##AO thref=_constants.thref thref=1e-3 montgc=psikk+(p[kdm,:,:]*(thkk-thstar[kdm-1,:,:]))*thref**2 for k in reversed(range(kdm-1)) : - montgc [:,:]=montgc [:,:]+ p[k+1,:,:]*(thstar[k+1,:,:]-thstar[k,:])*thref**2 + montgc [:,:]=montgc [:,:]+ p[k+1,:,:]*(thstar[k+1,:,:]-thstar[k,:])*thref**2 return montgc -def main(archv_files,opath,header_line1=cline1,header_line2=cline2,header_line3=cline3 ) : - - +### AO + + +def approx_montg1(thstar,p,srfhgt,idm,jdm,rstr_file): + # estimate montg1--------------------------------------------- + # Read input variables from lowest layer of a restart file + thref=1e-3 + #psikk_file="./TP6restart.2021_041_00_0000.a" + psikk_file=rstr_file#[0]#[-30:] + logger.info("----->rstr_file :%s"%psikk_file) + m=re.match( "^(.*)(\.[ab])", psikk_file) + if m : psikk_file=m.group(1) + logger.info("----->rstr_file :%s"%psikk_file) + rfile=abfile.ABFileRestart(psikk_file,"r",idm=idm,jdm=jdm) + psikk=rfile.read_field("psikk",0) + thkk=rfile.read_field("thkk",0) + rfile.close() + #print(( "psikk.min(),psikk.max()=",psikk.min(),psikk.max()) + #print( "thkk.min(),thkkk.max()=",thkk.min(),thkk.max() + # + #montg1pb = modeltools.hycom._montg_tools.montg1_pb(thstar,p) + #montg1c = modeltools.hycom._montg_tools.montg1_no_pb(psikk,thkk,thstar,p) + montg1pb = montg1_pb(thstar,p) + montg1c = montg1_no_pb(psikk,thkk,thstar,p) + print( "montg1c.min(),montg1c.max()=",montg1c.min(),montg1c.max()) + print( "montg1pb.min(),montg1pb.max()=",montg1pb.min(),montg1pb.max() ) + pbavg = (srfhgt - montg1c) / (montg1pb+thref) + print( "pbavg.min(),pbavg.max()=",pbavg.min(),pbavg.max() ) + montg1 = montg1pb*pbavg + montg1c + logger.info("Estimated montg1 ") + print( " min max montg1=",montg1.min(),montg1.max() ) + # end estimate montg1--------------------------------------------- + return montg1 + + + +def main(archv_files,restart_file,out_path) : +#def main(dtlist,archv_files,restart_file,archv_basedir="") : +# main([],args.archive_file,archv_basedir=args.archv_basedir) + + # Set paths using REGION.src environment (must be available at python invocation) + # Open blkdat files. Get some properties + header_line1="HYCOM nested archive files" + header_line2="Compute montg1 from Nemo ssh" + header_line3="Expt 01.7 nhybrd=50 nsigma= 0 ds00= 1.00 dp00= 1.00 dp00x= 450.0 dp00f=1.150" + # bp=modeltools.hycom.BlkdatParser("blkdat.input") idm = bp["idm"] jdm = bp["jdm"] kdm = bp["kdm"] thflag = bp["thflag"] - thbase = numpy.float(bp["thbase"]) + thbase = bp["thbase"] kapref = bp["kapref"] iversn = bp["iversn"] - iexpt = bp["iexpt"] yrflag = bp["yrflag"] - #thref=modeltools.hycom.thref + iexpt = bp["iexpt"] thref=1e-3 - if kapref == -1 : - kapnum = 2 - msg="Only kapref>=0 is implemented for now" - logger.error(msg) - raise ValueError(msg) - else : - kapnum = 1 - - if kapnum > 1 : - msg="Only kapnum=1 is implemented for now" - logger.error(msg) - raise ValueError(msg) - - - # hycom sigma and kappa, written in python. NB: sigver is not used here. - # Modify to use other equations of state. For now we assume sigver is: - # 1 (7-term eqs referenced to 0 bar) - # 2 (7-term eqs referenced to 2000 bar) - if thflag == 0 : - sigver=1 - else : - sigver=2 - sig = modeltools.hycom.Sigma(thflag) - if kapref > 0 : kappa = modeltools.hycom.Kappa(kapref,thflag*1000.0e4) # - - - - # Loop through archive files + # First read nemo ssh from the meain files + #archvf_00=archv_files[0][-19:-4]+"00.a" + for archv_file in archv_files : - archv_file=archv_file[:-2] - logger.debug("Processing %s"%(archv_file)) - arcfile=abf.ABFileArchv(archv_file,"r") - - temp = numpy.ma.zeros((jdm,idm)) # Only needed when calculating density - saln = numpy.ma.zeros((jdm,idm)) # Only needed when calculating density + archvf_00=archv_file + logger.info("----->Processing archv file :%s"%archvf_00) + logger.info("----->Restart file _pthkk :%s"%restart_file) + logger.info("Daily mean nemo file:%s"%archvf_00) + if os.path.isfile(archvf_00): + logger.info("file %s OK:"%archvf_00) + else: + logger.info("Ooobs file does not exist %s:"%archvf_00) + print( "" ) + ##logger.debug("Reading nemo_srfhgt from %s"%(archvf_00)) + arcfile_in=abfile.ABFileArchv(archvf_00,"r") + nemo_srfhgt=arcfile_in.read_field("srfhgt",0) + #nemo_montg1=arcfile_in.read_field("montg1",thflag) + # --------------------------------------------------- + logger.info("compute thstar and p for montg1") + sig = modeltools.hycom.Sigma(thflag) + # Read all layers .. (TODO: If there are memory problems, read and estimate sequentially) + temp = numpy.ma.zeros((jdm,idm)) # Only needed when calculating density + saln = numpy.ma.zeros((jdm,idm)) # Only needed when calculating density th3d =numpy.ma.zeros((kdm,jdm,idm)) thstar=numpy.ma.zeros((kdm,jdm,idm)) dp =numpy.ma.zeros((jdm,idm)) p =numpy.ma.zeros((kdm+1,jdm,idm)) - pup =numpy.ma.zeros((jdm,idm)) - montg=numpy.ma.zeros((jdm,idm)) + logger.info("---------------") + logger.info("-----Reading layers to get thstar and p----------") for k in range(kdm) : - logger.info("Reading layer %d from %s"%(k,archv_file)) - temp =arcfile.read_field("temp",k+1) - saln =arcfile.read_field("salin",k+1) - dp [:,:]=arcfile.read_field("thknss",k+1) - th3d [k ,:,:]=sig.sig(temp,saln) - thbase + ##logger.debug("Reading layer %d from %s"%(k,archvf_00)) + temp =arcfile_in.read_field("temp",k+1) + saln =arcfile_in.read_field("salin",k+1) + #dp [k ,:,:]=arcfile.read_field("thknss",k+1) + dp [:,:]=arcfile_in.read_field("thknss",k+1) + th3d [k ,:,:]=sig.sig(temp,saln)-float(thbase) + #dens_sig=sig.sig(temp,saln) p [k+1,:,:]= p[k,:,:] + dp[:,:] - if k>0 : thstarup =numpy.ma.copy(thstar[k ,:,:]) thstar[k ,:,:]=numpy.ma.copy(th3d [k ,:,:]) if kapref > 0 : - thstar[k ,:,:]=thstar [k ,:,:] + kappa.kappaf(temp[:,:], - saln[:,:], - th3d[k,:,:]+thbase, - p[k,:,:]) + thstar[k ,:,:]=thstar [k ,:,:] + kappa.kappaf( + temp[:,:], saln[:,:], th3d[k,:,:]+float(thbase), p[k,:,:]) elif kapref < 0 : msg="Only kapref>=0 is implemented for now" logger.error(msg) raise ValueError(msg) - if k > 0 : - montg = montg - p[k,:,:] * (thstar[k ,:,:] - thstarup) * thref**2 - - thkk = thstar[k ,:,:] - psikk = montg - - # This part of montg1 does nto depend on pbavg - montg1c = montg1_pb(thstar,p) - # This part depends linearly on pbavg - montg1pb = montg1_no_pb(psikk,thkk,thstar,p) - # Barotropic pressure. NB: srfhgt is in geopotential height : - srfhgt=arcfile.read_field("srfhgt",0) - pbavg = (srfhgt - montg1c) / (montg1pb+thref) - montg1 = montg1pb*pbavg + montg1c - logger.info("Estimated montg1 ") - - # Apply mask - #montg1=numpy.ma.masked_where(srfhgt.mask,montg1) - #montg1[~msk]=2**100 - #montg1 = numpy.ma.masked_where(srfhgt<2.0**99,montg1) - msk=numpy.where(srfhgt>2.0**99,1,0) - - - # Open new archive file, and write montg1 to it. - fnameout = opath+str(archv_file[-17:]) - arcfile_out=abf.ABFileArchv(fnameout,"w",iversn=iversn,yrflag=yrflag,iexpt=iexpt,mask=False,cline1=header_line1,cline2=header_line2,cline3=header_line3) - - fields=arcfile.get_fields() - for key in sorted(fields.keys()) : - fieldname = fields[key]["field"] - time_step = fields[key]["step"] - model_day = fields[key]["day"] - k = fields[key]["k"] - dens = fields[key]["dens"] - fld =arcfile.read_field(fieldname,k) - - if fieldname == "montg1" : - logger.info("Writing field %10s at level %3d to %s (modified)"%(fieldname,k,fnameout)) - arcfile_out.write_field(montg1,None,fieldname,time_step,model_day,k,dens) - else : - arcfile_out.write_field(fld ,None,fieldname,time_step,model_day,k,dens) - - logger.info("Finished writing to %s"%fnameout) + # --------------------------------------------------- + #logger.info("Adding nemo ssh to the tide elevation:%s"%archv_file_in_string) + logger.info("---------------------------------") + archv_file_out = out_path+str(archvf_00[-19:]) + logger.info("----->Adding montg1 archv_file_out:%s"%archv_file_out) + logger.info("Output file %s"%(archv_file_out)) + arcfile_out=abfile.ABFileArchv(archv_file_out,"w",iversn=iversn,yrflag=yrflag,iexpt=iexpt,mask=False,cline1=header_line1,cline2=header_line2,cline3=header_line3) + # + ###AA regr = open("montg_regress_tide.pckl",'rb') + ###AA a,b,nemomeandt = pickle.load(regr) + ###AA tide_montg1 = (nemo_srfhgt - nemomeandt) * a + b + estimated_montg1=approx_montg1(thstar,p,nemo_srfhgt,idm,jdm,restart_file) + #---# + #read all keys and write updated version + for keys in sorted(arcfile_in.fields.keys()) : + fieldname = arcfile_in.fields[keys]["field"] + time_step = arcfile_in.fields[keys]["step"] + model_day = arcfile_in.fields[keys]["day"] + k = arcfile_in.fields[keys]["k"] + dens = arcfile_in.fields[keys]["dens"] + fld = arcfile_in.read_field(fieldname,k) + #logger.info("Current field %10s at level %3d, and dens %s"%(fieldname,k,dens)) + if fieldname == "montg1" : + logger.info("Writing field %10s at level %3d to %s (modified)"%(fieldname,k,archv_file_out)) + arcfile_out.write_field(estimated_montg1,None,fieldname,time_step,model_day,k,dens) + #elif fieldname == "srfhgt" : + # logger.info("Writing field %10s at level %3d to %s (modified)"%(fieldname,k,archv_file_out)) + # #arcfile_out.write_field(tide_srfhgt,None,fieldname,time_step,model_day,k,dens) + else : + arcfile_out.write_field(fld ,None,fieldname,time_step,model_day,k,dens) + + arcfile_in.close() arcfile_out.close() - arcfile.close() + + # Remove archive files lying around + #logger.info("Finished - files placed in location %s"%os.path.abspath("./")) + + if __name__ == "__main__" : + basecmd=os.path.basename(sys.argv[0]) - parser = argparse.ArgumentParser(description='Calculate montgomery potential in 1st layer (used in archv files)') - parser.add_argument('--header_line1', type=str, help="First header line, applicable only for archm & archv",default=cline1) - parser.add_argument('--header_line2', type=str, help="Montgomery with or without barotropic pressure",default=cline2) - parser.add_argument('--header_line3', type=str, help="Montgomery with or without barotropic pressure",default=cline3) + class DateTimeParseAction(argparse.Action) : + def __call__(self, parser, args, values, option_string=None): + tmp = datetime.datetime.strptime( values, "%Y-%m-%dT%H:%M:%S") + setattr(args, self.dest, tmp) + + + parser = argparse.ArgumentParser(description='') + #parser.add_argument("--archv-basedir",type = str, help = "",default="") + parser.add_argument('archive_file', type=str, nargs="+") + parser.add_argument('restart_file', type=str, help="Restart file produced by model run") + parser.add_argument('out_path', type=str, help='Output files path') + #parser_opt2.set_defaults(subparser_name="2") - parser.add_argument('archive_file', type=str, help='Files to add montg1 to',nargs="+") - parser.add_argument('opath', type=str, help='Output files path') args = parser.parse_args() - main(args.archive_file,args.opath,header_line1=args.header_line1,header_line2=args.header_line2,header_line3=args.header_line3) + + main(args.archive_file,args.restart_file,args.out_path) diff --git a/bin/expt_preprocess_calc_Montg.sh b/bin/expt_preprocess_calc_Montg.sh new file mode 100755 index 00000000..f1a8118c --- /dev/null +++ b/bin/expt_preprocess_calc_Montg.sh @@ -0,0 +1,777 @@ +#!/bin/bash +set -x + +# Initialize environment (sets Scratch dir ($S), Data dir $D ++ ) +if [ -f EXPT.src ] ; then + export BASEDIR=$(cd .. && pwd) +else + echo "Could not find EXPT.src. This script must be run in expt dir" + exit 1 +fi +export BINDIR=$(cd $(dirname $0) && pwd)/ +source ../REGION.src || { echo "Could not source ../REGION.src "; exit 1; } +source ./EXPT.src || { echo "Could not source EXPT.src "; exit 1; } +source $BINDIR/common_functions.sh || { echo "Could not source common_functions.sh "; exit 1; } + +# KAL - new input start and end date times in ISO format +if [ $# -lt 2 ] ; then + tellerror " Need start and stop times as input" + exit 1 +else + starttime=$1 + endtime=$2 + initstr="" + if [ $# -eq 3 ] ; then + initstr="$3" + fi +fi +echo "Start time is $starttime" +echo "Stop time is $endtime" + +# Parse starttime and endtime +if [[ $starttime =~ ([0-9]{4})-([0-9]{2})-([0-9]{2})T([0-9]{2}):([0-9]{2}):([0-9]{2}) ]] ; then + start_year=${BASH_REMATCH[1]} + echo ${start_year} + start_month=${BASH_REMATCH[2]} + start_day=${BASH_REMATCH[3]} + start_hour=${BASH_REMATCH[4]} + start_min=${BASH_REMATCH[5]} + start_sec=${BASH_REMATCH[6]} + start_hsec=$(echo ${start_min}\*60+${start_sec} | bc ) + start_hsec=$(echo 000${start_hsec} | tail -c5) + start_dsec=$(echo ${start_hour}\*3600+${start_hsec} | bc ) + start_dsec=$(echo 0000${start_dsec} | tail -c6) + start_oday=$(date -u -d "${start_year}-${start_month}-${start_day} 00:00:00 UTC" +%j) +else + tellerror "start time not in righ format" ; exit 1 +fi +if [[ $endtime =~ ([0-9]{4})-([0-9]{2})-([0-9]{2})T([0-9]{2}):([0-9]{2}):([0-9]{2}) ]] ; then + end_year=${BASH_REMATCH[1]} + end_month=${BASH_REMATCH[2]} + end_day=${BASH_REMATCH[3]} + end_hour=${BASH_REMATCH[4]} + end_min=${BASH_REMATCH[5]} + end_sec=${BASH_REMATCH[6]} + end_hsec=$(echo ${end_min}\*60+${end_sec} | bc ) + end_hsec=$(echo 0000${end_hsec} | tail -c5) + end_dsec=$(echo ${end_hour}\*3600+${end_hsec} | bc ) + end_dsec=$(echo 0000${end_dsec} | tail -c6) + end_oday=$(date -u -d "$end_year-$end_month-$end_day 00:00:00 UTC" +%j) +else + tellerror "end time not in righ format" ; exit 1 +fi +#echo hsec $start_hsec $end_hsec +#echo $starttime +#echo dsec $start_dsec $end_dsec +#echo oday $start_oday $end_oday +#exit + +# Init error counter (global var used by function tellerror) +numerr=0 + +# Create data and scratch dirs. Enter scratch dir +if [ ! -d $D ] ; then + mkdir -p $D || { echo "Could not create data dir : $S " ; exit 1 ; } +fi +if [ ! -d $S ] ; then + mkdir -p $S || { echo "Could not create data dir : $S " ; exit 1 ; } +fi +cd $S || { tellerror "no scratch dir $S" ; exit 1 ;} + + + +# --- fetch some things from blkdat.input +cp $P/blkdat.input blkdat.input || tellerror "No blkdat.input file" +export LBFLAG=`grep "'lbflag' =" blkdat.input | awk '{printf("%03d", $1)}'` +export EB=`grep "'iexpt ' =" blkdat.input | awk '{printf("%03d", $1)}'` +export PRIVER=`grep "'priver' =" blkdat.input | awk '{printf("%1d", $1)}'` +export NTRACR=`grep "'ntracr' =" blkdat.input | awk '{printf("%03d", $1)}'` +export YRFLAG=`grep "'yrflag' =" blkdat.input | awk '{printf("%1d", $1)}'` +export JERLV=`grep "'jerlv0' =" blkdat.input | awk '{printf("%1d", $1)}'` +export SSSRLX=`grep "'sssflg' =" blkdat.input | awk '{printf("%1d", $1)}'` +export SSTRLX=`grep "'sstflg' =" blkdat.input | awk '{printf("%1d", $1)}'` +export RLX=`grep "'relax ' =" blkdat.input | awk '{printf("%1d", $1)}'` +export TRCRLX=`grep "'trcrlx' =" blkdat.input | awk '{printf("%1d", $1)}'` +export THKDF4=`grep "'thkdf4' =" blkdat.input | awk '{printf("%f", $1)}'` +export VELDF4=`grep "'veldf4' =" blkdat.input | awk '{printf("%f", $1)}'` +export KAPREF=`grep "'kapref' =" blkdat.input | awk '{printf("%f", $1)}'` +export VSIGMA=`grep "'vsigma' =" blkdat.input | awk '{printf("%1d", $1)}'` +export FLXOFF=`grep "'flxoff' =" blkdat.input | awk '{printf("%1d", $1)}'` +export STDFLG=`grep "'stdflg' =" blkdat.input | awk '{printf("%1d", $1)}'` +export BNSTFQ=$(blkdat_get blkdat.input bnstfq) +export NESTFQ=$(blkdat_get blkdat.input nestfq) +export THKDF2=$(blkdat_get blkdat.input thkdf2) +export THFLAG=$(blkdat_get blkdat.input thflag) +export BACLIN=$(blkdat_get blkdat.input baclin) +export BATROP=$(blkdat_get blkdat.input batrop) +export CPLIFQ=$(blkdat_get blkdat.input cplifq) +export ICEFLG=$(blkdat_get blkdat.input iceflg) +export MOMTYP=$(blkdat_get blkdat.input momtyp) +export VISCO2=$(blkdat_get blkdat.input visco2) +export VELDF2=$(blkdat_get blkdat.input veldf2) +export IDM=$(blkdat_get blkdat.input idm) +export JDM=$(blkdat_get blkdat.input jdm) +export LWFLAG=`grep "'lwflag' =" blkdat.input | awk '{printf("%1d", $1)}'` + +restarti=$(blkdat_get_string blkdat.input nmrsti "restart_in") + +# Add period to restart file name if not present... +if [ $restarti != "restart_in" -a ${restarti: -1} != "." ] ;then + restarti="${restarti}." +fi + + +# +# --- Set up time limits +# +cmd="$BINDIR/hycom_limits.py $starttime $endtime $initstr" +echo "*Setting up HYCOM time limits : $cmd " +eval $cmd || tellerror "$cmd failed" +if [ $ICEFLG -eq 2 ] ; then + cmd="$BINDIR/cice_limits.py $initstr $starttime $endtime $NMPI $P/ice_in" + echo "*Setting up CICE time limits : $cmd " + eval $cmd || tellerror "$cmd failed" +fi + +# --- fetch some things from the ice model +if [ $ICEFLG -eq 2 ] ; then + #export icedt=`egrep "^[ ,]*dt[ ]*=" ../ice_in | sed "s/.*=//" | sed "s/,.*//"` + export icedt=$($HYCOM_PYTHON_ROUTINES/namelist_extract.py ice_in setup_nml dt) + export ice_restart_dir=$($HYCOM_PYTHON_ROUTINES/namelist_extract.py ice_in setup_nml restart_dir) + export ice_restart_file=$($HYCOM_PYTHON_ROUTINES/namelist_extract.py ice_in setup_nml restart_file) + export ice_restart_pointer_file=$($HYCOM_PYTHON_ROUTINES/namelist_extract.py ice_in setup_nml pointer_file) +fi + + +# --- check that iexpt from blkdat.input agrees with E from this script. +[ "$EB" == "$E" ] || tellerror " blkdat.input iexpt ${EB} different from this experiment ${E} set in EXPT.src" +echo "Fetched from blkdat.input:" +echo "--------------------------" +echo "EB is $EB " +echo "PRIVER is $PRIVER" +echo "YRFLAG is $YRFLAG" +echo "JERLV is $JERLV " +echo "SSS is $SSSRLX" +echo "SST is $SSTRLX" +echo "BNSTFQ is $BNSTFQ" +echo "NESTFQ is $NESTFQ" +echo "FLXOFF is $FLXOFF" +echo "STDFLG is $STDFLG" +echo "--------------------" + +# Check baroclinic time step +if [ $YRFLAG -ge 1 ] ; then + tmp2=21600. +else + tmp2=86400. +fi +tmp=$(echo $tmp2"%"$BACLIN | bc) +testbc=$(echo $tmp"=="0 | bc ) +if [ $testbc -ne 1 ] ; then + tellerror "$tmp2 seconds is not divisible with baclin=$BACLIN". +fi + +# Check barotropic time step +tmp=$(echo $BACLIN"/"$BATROP | bc -l) +testbt=$(echo $tmp"%2==0" | bc -l ) +if [ $testbt -ne 1 ] ; then + tellerror "($BACLIN / $BATROP ) %2 not zero ". +fi + + +# Check coupling time step +if [ $ICEFLG -eq 2 ] ; then + # Check that icedt = cplifq*baclin + test1=$(echo ${CPLIFQ}'<'0.0 | bc -l) + + # cplifq negative = number of time steps between coupling + if [ $test1 -eq 1 ] ; then + dtcpl=$(echo "-1.*"$CPLIFQ"*"$BACLIN | bc -l) + # cplifq positive = fraction of days + else + dtcpl=$(echo $CPLIFQ"*"86400. | bc -l) + fi + echo "Coupling time step=$dtcpl" + + # Check ice model dt against $dtcpl. dtcpl must be a multiple of icedt + test2=$(echo "scale=0;"$dtcpl"%"$icedt | bc -l) + test2=$(echo ${test2}'=='0.0 | bc -l) + if [ $test2 -eq 0 ] ; then + tellerror "coupling time step $dtcpl not an integer multiple of ice model time step $icedt" + fi +fi + + +# TODO: Limitation for now. Note that newest hycom can initialize with yrflag=3 (realistic forcing) +if [ $YRFLAG -ne 3 ] ; then + tellerror "must use yrflag=3 in blkdat.input" + exit 1 +fi + +# Check momtyp. veldf2 and visco2 must be zero in this case +if [ $MOMTYP -eq 4 ] ; then + test1=$(echo ${VISCO2}'!='0.0 | bc -l) + test2=$(echo ${VELDF2}'!='0.0 | bc -l) + if [ $test1 -ne 0 -o $test2 -ne 0 ] ; then + tellerror "momtyp=4; visco2 and veldf2 must be zero" + fi +fi + + + +touch ports.input tracer.input +rm ports.input tracer.input +[ -s $P/tracer.input ] && cp $P/tracer.input tracer.input + + +# Get init flag from start time +if [ "$initstr" == "--init" ] ;then + init=1 +else + init=0 +fi +tstart=$(cat limits | tr -s " " | sed "s/^[ ]*//" | cut -d " " -f1) +tstop=$(cat limits | tr -s " " | sed "s/^[ ]*//" | cut -d " " -f2) +echo "Fetched from limits:" +echo "--------------------" +echo "init is $init" +echo "tstart is $tstart" +echo "tstop is $tstop" +echo "--------------------" + +#C +#C --- turn on detailed debugging. +#C +#touch PIPE_DEBUG + +# +# --- pget, pput "copy" files between scratch and permanent storage. +# --- Can both be cp if the permanent filesystem is mounted locally. +# --- KAL - we use pget=pput=cp +export pget=/bin/cp +export pput=/bin/cp +export plink='ln -sf' + +echo "Initialization complete - now copying necessary files to scratch area" +echo + +# +# +# --- Set up MPI partition for this segment +# +# +#echo "we go to copy the damn file",$NMPI +[ "$NMPI" == "" ] && NMPI=0 # Init if unset +if [ $NMPI != 0 ] ; then + export NPATCH=`echo $NMPI | awk '{printf("%04d", $1)}'` + /bin/rm -f patch.input + echo /bin/cp $BASEDIR/topo/partit/depth_${R}_${T}.${NPATCH} patch.input + /bin/cp $BASEDIR/topo/partit/depth_${R}_${T}.${NPATCH} patch.input || \ + tellerror "Could not get patch.input ($BASEDIR/topo/partit/depth_${R}_${T}.${NPATCH})" +fi + +# +# --- input files from file server. +# +echo "**Retrieving grid and topography files" +${pget} $BASEDIR/topo/regional.grid.a regional.grid.a || tellerror "no grid file regional.grid.a" +${pget} $BASEDIR/topo/regional.grid.b regional.grid.b || tellerror "no grid file regional.grid.a" +${pget} $BASEDIR/topo/depth_${R}_${T}.a regional.depth.a || tellerror "no topo file depth_${R}_${T}.a" +${pget} $BASEDIR/topo/depth_${R}_${T}.b regional.depth.b || tellerror "no topo file depth_${R}_${T}.b" +${pget} $BASEDIR/topo/kmt_${R}_${T}.nc cice_kmt.nc || tellerror "no kmt file $BASEDIR/topo/kmt_${R}_${T}.nc " +${pget} $BASEDIR/topo/cice_grid.nc cice_grid.nc || tellerror "no cice grid file $BASEDIR/topo/cice_grid.nc " + + +if [ "$SSTRLX" -eq 3 ] ; then + [ -f $CLMDIR/seatmp.a ] || tellerror "File $CLMDIR/seatmp.a does not exist" + [ -f $CLMDIR/seatmp.b ] || tellerror "File $CLMDIR/seatmp.b does not exist" + ln -sf $CLMDIR/seatmp.a forcing.seatmp.a || tellwarn "Could not link $CLMDIR/seatmp.a" + ln -sf $CLMDIR/seatmp.b forcing.seatmp.b || tellwarn "Could not link $CLMDIR/seatmp.b" +fi + +# TODO: Not handled yet +if [ 0 -eq 1 ] ; then + [ -f $CLMDIR/surtmp.a ] || tellerror "File $CLMDIR/surtmp.a does not exist" + [ -f $CLMDIR/surtmp.b ] || tellerror "File $CLMDIR/surtmp.b does not exist" + ln -sf $CLMDIR/surtmp.a forcing.surtmp.a || tellwarn "Could not link $CLMDIR/surtmp.a" + ln -sf $CLMDIR/surtmp.b forcing.surtmp.b || tellwarn "Could not link $CLMDIR/surtmp.b" +fi + +# --- +# --- Pre-prepared forcing option - for now just sets up links +# --- TODO: To verify, we have to go into file and check timings +# --- +echo "**Setting up pre-prepared synoptic forcing from force/synoptic/$E" +DIR=$BASEDIR/force/synoptic/$E/ +echo $DIR + +declare -a arr=("radflx" "shwflx" "vapmix" "airtmp" "precip" "mslprs" "wndewd" "wndnwd" "dewpt") +#for i in tauewd taunwd wndspd radflx shwflx vapmix \ +# airtmp precip uwind vwind clouds relhum slp ; do + +#for i in radflx shwflx vapmix \ +# airtmp precip mslprs \ +# wndewd wndnwd ; do +for i in "${arr[@]}"; do + echo "|--> $i" + [ -f $DIR/$i.a ] || tellerror "File $DIR/$i.a does not exist" + [ -f $DIR/$i.b ] || tellerror "File $DIR/$i.b does not exist" + + # Check range of file against start and stop times + if [ -f $DIR/$i.a -a -f $DIR/$i.b ] + then + [ ! -s $DIR/${i}.a ] && tellerror "$DIR/$i.a: File is empty" + [ ! -s $DIR/${i}.b ] && tellerror "$DIR/$i.b: File is empty" + ln -sf $DIR/${i}.a forcing.${i}.a || tellerror "Could not fetch $DIR/$i.a" + ln -sf $DIR/${i}.b forcing.${i}.b || tellerror "Could not fetch $DIR/$i.b" + + frcstart=$(head -n 6 forcing.$i.b | tail -n1 | sed "s/.*=//" | sed "s/^[ ]*//" | cut -d " " -f 1) + frcstop=$(tail -n 1 forcing.$i.b | sed "s/.*=//" | sed "s/^[ ]*//" | cut -d " " -f 1) + test1=$(echo ${tstart#-}'>='$frcstart | bc -l) + test2=$(echo $tstop '<='$frcstop | bc -l) + [ $test1 -eq 1 ] || tellerror "File $S/forcing.$i.b: forcing starts after model starts" + [ $test2 -eq 1 ] || tellerror "File $S/forcing.$i.b: forcing stops before model stops" + fi + + +done + + +# MOSTAFA: END +# +# --- time-invarent heat flux offset +# +#TODO handle offsets +#setenv OF "" +#setenv OF "_413" +if [ "$CLMDIR" != "" ] ; then + if [ "$OF" != "" ]; then + touch forcing.offlux.a forcing.offlux.b + if [ ! -s forcing.offlux.a ] ; then + ${pget} $CLMDIR/offlux${OF}.a forcing.offlux.a & + fi + if [ ! -s forcing.offlux.b ] ; then + ${pget} $CLMDIR/offlux${OF}.b forcing.offlux.b & + fi + fi +fi + +# MOSTAFA: BEGIN +# For time-invariant offlux +# copy flux off set files if flxoff=1 +echo "FLXOFF = $FLXOFF" +if [ $FLXOFF -eq 1 ] ; then + echo "====================================================" + echo " -------flux off set true: copy flux off set files-" + cp ${D}/../../relax/${E}/offlux.a forcing.offlux.a || tellerror "Could not get offlux .a file" + cp ${D}/../../relax/${E}/offlux.b forcing.offlux.b || tellerror "Could not get offlux .b file" + echo "====================================================" + else + echo "fLxoff=F: No attempt to use flux offset correction" +fi + + +# MOSTAFA: END + +# +# --- river forcing +# --- KAL: rivers are experiment-dependent +# +if [ $PRIVER -eq 0 ] ; then +echo "**No river forcing. Set the priver to 1 to add river forcing" +fi +if [ $PRIVER -eq 1 ] ; then + echo "**Setting up river forcing from priver" + cp $BASEDIR/force/rivers/$E/rivers.a forcing.rivers.a || tellerror "Could not get river .a file" + cp $BASEDIR/force/rivers/$E/rivers.b forcing.rivers.b || tellerror "Could not get river .b file" + if [ $NTRACR -ne 0 ] ; then + echo "**Setting up bio river forcing" + cp $BASEDIR/force/rivers/$E/ECO_no3.a rivers.ECO_no3.a || tellwarn "Could not get NO3 river .a file" + cp $BASEDIR/force/rivers/$E/ECO_no3.b rivers.ECO_no3.b || tellwarn "Could not get NO3 river .b file" + cp $BASEDIR/force/rivers/$E/ECO_sil.a rivers.ECO_sil.a || tellwarn "Could not get SIL river .a file" + cp $BASEDIR/force/rivers/$E/ECO_sil.b rivers.ECO_sil.b || tellwarn "Could not get SIL river .b file" + cp $BASEDIR/force/rivers/$E/ECO_pho.a rivers.ECO_pho.a || tellwarn "Could not get PHO river .a file" + cp $BASEDIR/force/rivers/$E/ECO_pho.b rivers.ECO_pho.b || tellwarn "Could not get PHO river .b file" + fi +fi + + +# +# --- kpar forcing +# +if [ $JERLV -eq 0 ] ; then + echo "**Setting up kpar forcing" + ln -sf $BASEDIR/force/seawifs/kpar.a forcing.kpar.a || tellerror "Could not get kpar.a file" + ln -sf $BASEDIR/force/seawifs/kpar.b forcing.kpar.b || tellerror "Could not get kpar.b file" +fi + + +# +# --- relaxation +# +if [ $SSSRLX -eq 1 -o $SSTRLX -eq 1 -o $RLX -eq 1 -o $init -eq 1 ] ; then + echo "**Setting up relaxation" + for i in saln temp intf ; do + j=$(echo $i | head -c3) + [ ! -f $BASEDIR/relax/${E}/relax_$j.a ] && tellerror "$BASEDIR/relax/${E}/relax_$j.a does not exist" + [ ! -f $BASEDIR/relax/${E}/relax_$j.b ] && tellerror "$BASEDIR/relax/${E}/relax_$j.b does not exist" + ln -sf $BASEDIR/relax/${E}/relax_$j.a relax.$i.a || tellerror "Could not get relax.$i.a" + ln -sf $BASEDIR/relax/${E}/relax_$j.b relax.$i.b || tellerror "Could not get relax.$i.b" + done +fi +if [ $RLX -eq 1 ] ; then + echo "**Setting up relaxation masks" + [ ! -f $BASEDIR/relax/${E}/relax_rmu.a ] && tellerror "$BASEDIR/relax/${E}/relax_rmu.a does not exist" + [ ! -f $BASEDIR/relax/${E}/relax_rmu.b ] && tellerror "$BASEDIR/relax/${E}/relax_rmu.b does not exist" + ln -sf $BASEDIR/relax/${E}/relax_rmu.a relax.rmu.a || tellerror "Could not get relax.rmu.a" + ln -sf $BASEDIR/relax/${E}/relax_rmu.b relax.rmu.b || tellerror "Could not get relax.rmu.b" +fi +# +# --- tracer relaxation +# +if [ $TRCRLX -ne 0 -o $NTRACR -eq -1 ] ; then + echo "**Setting up tracer relaxation" + for i in ECO_no3 ECO_pho ECO_sil ECO_oxy CO2_dic CO2_alk; do + j=$(echo $i | head -c7) + [ ! -f $BASEDIR/relax/${E}/relax.$j.a ] && tellerror "$BASEDIR/relax/${E}/relax.$j.a does not exist" + [ ! -f $BASEDIR/relax/${E}/relax.$j.b ] && tellerror "$BASEDIR/relax/${E}/relax.$j.b does not exist" + ln -sf $BASEDIR/relax/${E}/relax.$j.a relax.$i.a || tellerror "Could not get relax.$i.a" + ln -sf $BASEDIR/relax/${E}/relax.$j.b relax.$i.b || tellerror "Could not get relax.$i.b" + done + echo "**Setting up tracer relaxation masks" + [ ! -f $BASEDIR/relax/${E}/relax_rmu.a ] && tellerror "$BASEDIR/relax/${E}/relax_rmutr.a does not exist" + [ ! -f $BASEDIR/relax/${E}/relax_rmu.b ] && tellerror "$BASEDIR/relax/${E}/relax_rmutr.b does not exist" + ln -sf $BASEDIR/relax/${E}/relax_rmu.a relax.rmutr.a || tellerror "Could not get relax.rmutr.a" + ln -sf $BASEDIR/relax/${E}/relax_rmu.b relax.rmutr.b || tellerror "Could not get relax.rmutr.b" + + [ ! -f $INPUTDIR/co2_annmean_gl.txt ] && tellerror "$INPUTDIR/co2_annmean_gl.txt does not exist" + ln -sf $INPUTDIR/co2_annmean_gl.txt co2_annmean_gl.txt || tellerror "Could not get co2_annmean_gl.txt" +fi +# +# - thermobaric reference state? +# +echo "**Setting up various forcing files" +[ -f tbaric.a ] && rm tbaric.a +[ -f tbaric.b ] && rm tbaric.b +if [ ${KAPREF:0:1} == "-" ] ;then + ${pget} $BASEDIR/topo/tbaric.a tbaric.a || tellerror "Could not get tbaric.a" + ${pget} $BASEDIR/topo/tbaric.b tbaric.b || tellerror "Could not get tbaric.b" +fi + +# +# Spatially varying isopycnal densities +# +if [ ${VSIGMA} -ne 0 ] ;then + ${pget} $BASEDIR/relax/${E}/iso_sigma.a iso.sigma.a || tellerror "Could not get iso.sigma.a" + ${pget} $BASEDIR/relax/${E}/iso_sigma.b iso.sigma.b || tellerror "Could not get iso.sigma.b" +fi + + +# Thickness diffusion +testthkdf2=$(echo $THKDF2'<'0.0 | bc -l) +testthkdf4=$(echo $THKDF4'<'0.0 | bc -l) +[ -f thkdf4.a ] && rm thkdf4.a +[ -f thkdf4.b ] && rm thkdf4.b +if [ ${testthkdf4} -eq 1 ] ; then + ${pget} ${D}/../../relax/${E}/thkdf4.a thkdf4.a || tellerror "Could not get thkdf4.a" + ${pget} ${D}/../../relax/${E}/thkdf4.b thkdf4.b || tellerror "Could not get thkdf4.b" +fi + +[ -f thkdf2.a ] && rm thkdf2.a +[ -f thkdf2.b ] && rm thkdf2.b +if [ ${testthkdf2} -eq 1 ] ; then + ${pget} ${D}/../../relax/${E}/thkdf2.a thkdf2.a || tellerror "Could not get thkdf2.a" + ${pget} ${D}/../../relax/${E}/thkdf2.b thkdf2.b || tellerror "Could not get thkdf2.b" +fi +testveldf4=$(echo $VELDF4'<'0.0 | bc -l) +if [ ${testveldf4} -eq 1 ] ; then + ${pget} ${D}/../../relax/${E}/veldf4.a veldf4.a || tellerror "Could not get thkdf2.a" + ${pget} ${D}/../../relax/${E}/veldf4.b veldf4.b || tellerror "Could not get thkdf2.b" +fi + + +# TODO Limited set of tests for now. +# Link in nest dir if nesting activated +#tmp=$(echo $BNSTFQ'>'0.0 | bc -l) +#tmp2=$(echo $NESTFQ'>'0.0 | bc -l) +tmp=$(echo $BNSTFQ'!='0.0 | bc -l) +tmp2=$(echo $NESTFQ'!='0.0 | bc -l) +if [ $tmp -eq 1 -o $tmp2 -eq 1 ] ; then + nestdir=$BASEDIR/nest/$E + echo "Nesting input from $nestdir" + ls nest + if [ -d $nestdir ] ; then + [ -e nest ] && rm nest + ln -s $nestdir nest + else + tellerror "Nesting dir $nest does not exist" + fi +fi + +#export waveSDIR=/work/shared/nersc/msc/STOKES/Globww3/tmp +#export waveSDIR=/work/shared/nersc/msc/STOKES/Globww3 +echo "====================================================" +echo "STDFLG = $STDFLG" +if [ $STDFLG -eq 1 ] ; then + echo "====================================================" + echo " -------Setting up Wave Stokes forcing----" + for foo in forcing.stokex.a forcing.stokex.b forcing.stokey.a forcing.stokey.b\ + forcing.transx.a forcing.transx.b forcing.transy.a forcing.transy.b\ + forcing.tauwx.a forcing.tauwx.b forcing.tauwy.a forcing.tauwy.b \ + forcing.twomx.a forcing.twomx.b forcing.twomy.a forcing.twomy.b \ + forcing.t01.a forcing.t01.b ; do + + echo "|--> linking to $waveSDIR/${foo}" + ln -sf $waveSDIR/${foo} . || tellerror "Could not fetch stokes forcing" + done + echo "====================================================" + else + echo "STD=F: No attempt to link to the Stokes/Wave forcing in SCRATCH" +fi +#read -t 5 +# Need ports.input file in these cases +if [ $tmp -eq 1 -a $LBFLAG -ne 2 -a $LBFLAG -ne 4 ] ; then + tellerror "Must have lbflag = 2 or 4 when bnstfq <> 0.0 " +elif [ $tmp -eq 1 -a $LBFLAG -eq 1 ] ; then + # Port flow - file must be present in experiment dir + cp $P/ports.input . || tellerror "Could not get port file ${P}/ports.input for port flow" +elif [ $tmp -eq 1 -a $LBFLAG -eq 2 ] ; then + # Nest flow - use file in experiment dir if present. Otherwise look in nest dir + if [ -f $P/ports.nest ] ; then + echo "Using file $P/ports.nest for nesting: $P/ports.nest -> ./ports.input" + cp $P/ports.nest ports.input || tellerror "Could not get port file ${P}/ports.nest for nest flow" + elif [ -f $nestdir/ports.nest ] ; then + echo "Using file $nestdir/ports.nest for nesting: $P/ports.nest -> ./ports.input" + cp $nestdir/ports.nest ports.input || tellerror "Could not get port file ${nestdir}/ports.nest for nest flow" + else + tellerror "Could not get port file ports.nest in $P or ${nestdir} for nest flow" + fi +fi + +# Need nest rmu in this case: +if [ $tmp2 -eq 1 ] ; then + # Nest relaxation - use file in experiment dir if present. Otherwise look in nest dir +# if [ -f $P/rmu_nest.a -a -f $P/rmu_nest.a ] ; then +# echo "Using file $P/rmu_nest.[ab] for nesting relaxation: $P/rmu_nest.[ab] -> ./rmu.[ab]" +# cp $P/rmu_nest.a rmu.a || tellerror "Could not get port file ${P}/rmu_nest.a for nest relax" +# cp $P/rmu_nest.b rmu.b || tellerror "Could not get port file ${P}/rmu_nest.b for nest relax" +# elif [ -f $nestdir/rmu_nest.a -a -f $nestdir/rmu_nest.a ] ; then +# echo "Using file $nestdir/rmu_nest.[ab] for nesting: $nestdir/rmu_nest.[ab] -> ./rmu.[ab]" +# cp $nestdir/rmu_nest.a rmu.a || tellerror "Could not get port file ${nestdir}/rmu_nest.a for nest relax" +# cp $nestdir/rmu_nest.b rmu.b || tellerror "Could not get port file ${nestdir}/rmu_nest.b for nest relax" +# fi + if [ -f $nestdir/rmu.a -a -f $nestdir/rmu.b ] ; then + echo "Using file $nestdir/rmu.[ab] for nesting" + else + tellerror "Could not find files $nestdir/rmu.[ab] for nest relaxation" + fi +fi + + +## Retrieve nersc tidal data set if active +#if [ "$tideflag" == "T" ] ; then +# ${pget} ${BASEDIR}/tides_nersc/$E/${tidechoice}obc_elev.dat . || tellerror "Could not get tidal data ${tidechoice}obc_elev.dat " +#fi + + + + +echo + +# +# --- move old restart files to KEEP, typically from batch system rerun. +# +touch restart.dummy.b restart.dummy.a +for f in restart* ; do + /bin/mv $f KEEP/$f +done + +# +# --- try to get HYCOM restart from various areas +# +if [ $init -eq 1 ] ; then + echo "No restart needed" +else + + #HYCOM restart + filename=${restarti}${start_year}_${start_oday}_${start_hour}_${start_hsec} + echo $D/${filename}_mem001.a + + # Try to fetch restart from data dir $D + if [ -f $D/${filename}.a -a -f $D/${filename}.b ] ; then + echo "using HYCOM restart files ${filename}.[ab] from data dir $D" + cp $D/${filename}.a . + cp $D/${filename}.b . + + if [ $tmp -eq 1 -o $tmp2 -eq 1 ] ; then + echo "debug: $(pwd)" + ### "Calculate and Write Montgometry potential into the nesting files" + echo "filename="${filename} + for yy in `seq ${start_year} ${end_year}`; do + y=$(( $yy % 4 )) + if [ $y -eq 0 ] + then + echo "$yy is Leap Year!" + end_day=366 + else + end_day=365 + echo "$yy is not a Leap Year!" + fi + for dn in `seq -w ${start_oday} ${end_day}`; do + python ../calc_montg1.py /cluster/work/users/achoth/TP5a0.06/nest/080_NewMontg/archv.${yy}_${dn}_00.b /cluster/work/users/achoth/TP5a0.06/expt_08.1/data/${filename}.b ${nestdir}/ + done + done + echo " Nesting Files Modified Successfully " + fi + + elif [ -f $D/${filename}_mem001.a -a -f $D/${filename}_mem001.b ]; then + echo "using HYCOM restart files ${filename}_mem???.[ab] from data dir $D" + for f in ${plink} $D/${filename}_mem*.? ; do + ${plink} $f . + done + + else + tellerror "Could not find HYCOM restart file ${filename}.[ab] in $D" + fi + + #CICE restart + if [ $ICEFLG -eq 2 ] ; then + filenameice="${ice_restart_dir}/${ice_restart_file}.${start_year}-${start_month}-${start_day}-${start_dsec}" + + # Try to fetch restart from data dir $D + if [ -f $D/${filenameice}.nc ] ; then + echo "using CICE restart file ${filenameice}.nc from data dir $D" + cp $D/${filenameice}.nc ${filenameice}.nc + echo ${filenameice}.nc > ${ice_restart_pointer_file} + + elif [ -f $D/${filenameice}_mem001.nc ]; then + echo "using CICE restart file ${filenameice}_mem???.nc from data dir $D" + for f in $D/${filenameice}_mem*.nc ; do + ${plink} $f cice/. + done + echo ${filenameice}_mem000.nc > ${ice_restart_pointer_file} + + else + tellerror "Could not find CICE restart file ${filenameice} in $D" + fi + fi + + +fi + + + +# +# --- model executable. One executable to rule them all (if it wasnt for CICE, that is...) +# +if [ $SIGVER -eq 1 ] ; then + TERMS=7 + MYTHFLAG=0 +elif [ $SIGVER -eq 2 ] ; then + TERMS=7 + MYTHFLAG=2 +elif [ $SIGVER -eq 3 ] ; then + TERMS=9 + MYTHFLAG=0 +elif [ $SIGVER -eq 4 ] ; then + TERMS=9 + MYTHFLAG=2 +elif [ $SIGVER -eq 5 ] ; then + TERMS=17 + MYTHFLAG=0 +elif [ $SIGVER -eq 6 ] ; then + TERMS=17 + MYTHFLAG=2 +else + echo "SIGVER = $SIGVER" + echo "So far only 7 term eq of state is supported (SIGVER=1 or 2) is supported" + exit 1 +fi +TERMS2=$(echo 0$TERMS | tail -c3) +echo "SIGVER = $SIGVER .There are $TERMS terms in equation of state" + +# Set up rel path and stmt fnc +compdir=$(source_dir $V $TERMS $THFLAG) +compdir=$P/build/${compdir} +if [ $ICEFLG -eq 2 ] ; then + echo "Retrieving hycom_cice from $compdir" + /bin/cp $compdir/hycom_cice . || tellerror "Could not get hycom_cice executable at " +elif [ $ICEFLG -eq 0 ] ; then + echo "Retrieving hycom_oasis from $compdir" + /bin/cp $compdir/hycom_oasis . || tellerror "Could not get hycom_oasis executable at " +fi + + +# +# --- summary printout +# +[ ! -d old ] && mkdir -p old +touch summary_out +/bin/mv summary_out old/summary_old + + +# +# --- heat transport output +# +touch flxdp_out.a flxdp_out.b +/bin/mv flxdp_out.a old/flxdp_old.a +/bin/mv flxdp_out.b old/flxdp_old.b +# +touch ovrtn_out +/bin/mv ovrtn_out old/ovrtn_old + + +# +# --- clean up old archive files, typically from batch system rerun. +# +[ ! -d KEEP ] && mkdir KEEP +touch archv.dummy.b archv.dummy.a +for f in arch* ; do + /bin/mv $f KEEP/ +done + + +# +# --- let all file copies complete. +# +wait +# +# --- zero file length means no rivers. +# +if [ ! -s forcing.rivers.a ] ; then + /bin/rm forcing.rivers.[ab] + echo "problem with rivers" + ls -l *rivers* +fi + + +#C +#C --- Nesting input archive files for next segment. +#C +#if (-e ./nest) then +# cd ./nest +# touch archv_${NB}.tar +# if (-z archv_${NB}.tar) then +# ${pget} ${D}/nest/archv_${NB}.tar archv_${NB}.tar & +# endif +# cd .. +#endif +#C +#chmod ug+x hycom +#/bin/ls -laFq +#C +#if (-e ./nest) then +# ls -laFq nest +#endif + + +if [ $numerr -eq 0 ] ; then + echo "No fatal errors. Ok to start model set up in $S" +else + echo "Some fatal errors occured. See above" +fi + +# Tell where stuff ended up +exit $numerr # Fails if any fatal errors occured + + diff --git a/bin/topaz_archv2restart_TP5.csh b/bin/topaz_archv2restart_TP5.csh new file mode 100755 index 00000000..e65c4e90 --- /dev/null +++ b/bin/topaz_archv2restart_TP5.csh @@ -0,0 +1,71 @@ +#!/bin/csh -x +# +set echo +# +# --- Form a HYCOM restart file from a HYCOM archive file. +# +set in_year = $1 +set in_day = $2 +set in_day = `printf '%03d' $in_day` + +set out_year = $1 +set out_day = $2 +set out_day = `printf '%03d' $out_day` + +set dir = `pwd` +set parentdir = `dirname $dir` +set foldername = `basename $parentdir` + +setenv RID = $foldername #TP5a0.06 +source $pwd/EXPT.src + +setenv R /cluster/work/users/achoth//${RID} +setenv D /cluster/work/users/achoth/${RID}/expt_${X}/data +setenv O /cluster/work/users/achoth/${RID}/expt_${X}/data/NewRestart + + +# +#mkdir -p ${O} +cd ${D} +# +touch regional.depth.a regional.depth.b +if (-z regional.depth.a) then + /bin/rm regional.depth.a + /bin/ln -s ${R}/topo/depth_${RID}_{T}.a regional.depth.a +endif +if (-z regional.depth.b) then + /bin/rm regional.depth.b + /bin/ln -s ${R}/topo/depth_${RID}_{T}.b regional.depth.b +endif +# +touch regional.grid.a regional.grid.b +if (-z regional.grid.a) then + /bin/rm regional.grid.a + /bin/ln -s ${R}/topo/regional.grid.a . +endif +if (-z regional.grid.b) then + /bin/rm regional.grid.b + /bin/ln -s ${R}/topo/regional.grid.b . +endif +# +# --- input archive file +# --- input restart file: (tempelate) can be any restart file from any date +# --- output restart file +# +###AO ####/bin/rm -f ${O}/restart.2022_270_00_0000.a ${O}/restart.2022_270_00_0000.b +#aprun -n 1 ~/HYCOM-tools/archive/src/archv2restart <