Adding dynamic ocean (LSG) and T42 atmosphere

.gitignore

@@ -18,3 +18,4 @@ plasim/run
 puma/bin
 puma/bld
 puma/run
+postprocessor/burn7.x

README.md

@@ -53,10 +53,10 @@ See the corresponding example below:
 Now, launch the model:
 
 ```bash
+screen
 cd /cofast/$USER/PlaSim
 source modules.load
 cd plasim/run
-screen
 ./most_plasim_run
 ```
 
@@ -69,24 +69,31 @@ Ctrl+a d
 > You can monitor the model process with
 
 ```bash
+cd /cofast/$USER/PlaSim/plasim/run
 cat plasim_diag | grep Completed
 ```
 
-Results are:
+For one year, results are:
 
 ```bash
-209M  MOST.001
-4,1M  MOST_REST.001
-249K  MOST_DIAG.001
-5,3M  MOST_OCE.001.nc
- 14M  MOST_ICE.001.nc
- 69M  ice_output
-209M  MOST.002
- 29M  ocean_output
-4,1M  MOST_REST.002
-245K  MOST_DIAG.002
-5,2M  MOST_OCE.002.nc
- 15M  MOST_ICE.002.nc
-4,1M  plasim_restart
+****************************************
+* User   time         :    244.036 sec *
+* System time         :      1.324 sec *
+* Total CPU time      :    245.360 sec *
+* Memory usage        :    120.780 MB  *
+* Page reclaims       :  25832 pages   *
+* Page faults         :   5185 pages   *
+* Disk read           : 416163 blocks  *
+* Disk write          :2726849 blocks  *
+****************************************
+* Seconds per sim year:    245         *
+* Sim years per day   :    353         *
+****************************************
+
+532K  MOST_DIAG.0001.txt
+ 18M  MOST_PLA.0001.nc
+1.3M  MOST_OCE.0001.nc
+2.2M  MOST_ICE.0001.nc
+ 26M  MOST_LSG.0001.nc
 ```
 

cca_files/PBS_plasim

@@ -0,0 +1,21 @@
+#!/bin/ksh
+#PBS -S /usr/bin/ksh
+#PBS -q ns
+
+#PBS -N plasim
+#PBS -o plasim.o
+#PBS -e plasim.e
+
+#PBS -m abe
+#PBS -M jerome.sauer@uclouvain.be
+
+#PBS -l walltime=48:00:00
+
+module load netcdf4
+module load netcdf
+module switch PrgEnv-cray PrgEnv-gnu
+module load cdo
+
+cd /scratch/ms/be/cvaj/plasim/run
+./most_plasim_run
+

cca_files/PBS_plasim_mpi

@@ -0,0 +1,25 @@
+#!/bin/ksh
+#PBS -S /usr/bin/ksh
+#PBS -q nf
+#PBS -l EC_total_tasks=16
+#PBS -l EC_hyperthreads=1
+
+#PBS -N plasim_test
+#PBS -o plasim_test.o
+#PBS -e plasim_test.e
+
+#PBS -m abe
+#PBS -M jerome.sauer@uclouvain.be
+
+#PBS -l walltime=48:00:00
+
+module switch PrgEnv-cray PrgEnv-gnu
+module load cray-snplauncher
+
+module load netcdf4
+module load netcdf
+module load cdo
+
+cd /scratch/ms/be/cvaj/plasim/run
+./most_plasim_run
+

cca_files/burner_markefile

@@ -0,0 +1,2 @@
+burn7.x:	burn7.cpp
+	c++ -O2 -o burn7.x burn7.cpp /home/ms/be/cu0e/lib/netcdf-4.2/lib/libnetcdf_c++.a -I/home/ms/be/cu0e/lib/netcdf-4.2/include/ -I/opt/cray/netcdf/4.3.0/GNU/48/include -L/home/ms/be/cu0e/lib/netcdf-4.2/lib -L/opt/cray/netcdf/4.3.0/GNU/48/lib -lm -lnetcdf

cca_files/crontab.cca.tmp

@@ -0,0 +1,3 @@
+0 2,5,8,11,14,17,20,23 * * * cd /home/ms/be/cu0e/models/PLASIM/run_scripts && ksh -i charnock0.02.sh
+0 1,4,7,10,13,16,19,22 * * * cd /home/ms/be/cu0e/models/PLASIM/run_scripts && ksh -i vonk0.42.sh
+

cca_files/crontab.ccb.tmp

@@ -0,0 +1,3 @@
+0 3,6,9,12,15,18,21 * * * cd /home/ms/be/cu0e/models/PLASIM/run_scripts && ksh -i standard.sh
+0 0 * * * rsync -a --delete /scratch/ms/be/cu0e/backup/* ecgate:/perm/ms/be/cu0e/plasim/backup/
+

cca_files/make_ensemble_experiment.py

@@ -0,0 +1,86 @@
+
+import os
+import sys
+
+# check the python version
+if float(sys.version[:3]) <= 3.6:
+	print("This script require Python 3.6 !")
+	print("Try:")
+	print("\n\tmodule load python3\n\n")
+	print("and restart this script again.")
+	print("Aborting...")
+	sys.exit(1)
+
+if len(sys.argv) < 4:
+	print('Bad arguments:', sys.argv)
+	print('Usage:')
+	print('\n\t python3 make_ensemble_experiment.py where experiment ensemble_size run_folder_name\n')
+	print('Arguments:\n')
+	print('\twhere :\t\t\tWhere the experiment ensemble folders are located.')
+	print('\texperiment :\t\tName of the experiment.')
+	print('\tensemble_size :\t\tSize of the ensemble.')
+	print('\trun_folder_name :\tOptional. Addtional name added to the "run_" folder name in plasim folder. Used to selected experiment executable. If not provided, assumes run folder is "run".')
+	sys.exit(0)
+
+mail = 'jerome.sauer@uclouvain.be'
+walltime = '48:00:00'
+
+home_dir = os.getenv("HOME")
+scratch_dir = os.getenv("SCRATCH")
+plasim_dir = home_dir + "/models/PLASIM/"
+
+where = sys.argv[1]
+experiment = sys.argv[2]
+ensemble_size = int(sys.argv[3])
+
+
+run_folder = plasim_dir+"plasim/run"
+
+if len(sys.argv) == 5:
+	run_folder_ext = sys.argv[4]
+	run_folder += '_'+run_folder_ext
+
+if not os.path.isdir(run_folder):
+	print("PLASIM run directory missing!")
+	print("Run "+plasim_dir+"/most.x first.")
+	print("Aborting...")
+	sys.exit(1)
+
+print('Creating the directory holding the experiment ensemble...')
+experiment_folder = where+"/"+experiment+"/"
+os.system('mkdir -p '+experiment_folder)
+
+print('Creating the ensemble folders...')
+for i in range(1, ensemble_size+1):
+	member_number = str(i).rjust(2, '0')
+	ensemble_member_folder = experiment_folder+'run_'+experiment+'_'+member_number
+	os.system('cp -r '+run_folder+' '+ensemble_member_folder)
+	job_name = 'plasim_'+experiment+'_'+member_number
+	f = open(ensemble_member_folder+'/PBS_'+job_name, 'w')
+	f.write("""#!/bin/ksh
+#PBS -S /usr/bin/ksh
+#PBS -q ns
+		   
+""")
+
+	f.write('#PBS -N '+job_name+'\n')
+	f.write('#PBS -o '+job_name+'.o\n')
+	f.write('#PBS -e '+job_name+'.e\n\n')
+
+	f.write('#PBS -m abe\n')
+	f.write('#PBS -M '+mail+'\n\n')
+	f.write('#PBS -l walltime='+walltime+'\n\n')
+
+	f.write("""module load netcdf4
+module load netcdf
+module switch PrgEnv-cray PrgEnv-gnu
+module load cdo
+
+""")
+	f.write('cd '+ensemble_member_folder+'\n')
+	f.write('./most_plasim_run\n')
+	f.close()
+
+print('Experiment ensemble created and available in the folder:')
+print(experiment_folder)
+

cca_files/make_mpi_ensemble_algorithm_experiment.py

@@ -0,0 +1,92 @@
+
+import os
+import sys
+
+# check the python version
+if float(sys.version[:3]) <= 3.6:
+	print("This script require Python 3. !")
+	print("Try:")
+	print("\n\tmodule load python3\n\n")
+	print("and restart this script again.")
+	print("Aborting...")
+	sys.exit(1)
+
+if len(sys.argv) < 4:
+	print('Bad arguments:', sys.argv)
+	print('Usage:')
+	print('\n\t python3 make_mpi_ensemble_algorithm_experiment.py where experiment ensemble_size run_folder_name\n')
+	print('Arguments:\n')
+	print('\twhere :\t\t\tWhere the experiment ensemble folders are located.')
+	print('\texperiment :\t\tName of the experiment.')
+	print('\tensemble_size :\t\tSize of the ensemble.')
+	print('\trun_folder_name :\tOptional. Addtional name added to the "run_" folder name in plasim folder. Used to selected experiment executable. If not provided, assumes run folder is "run".')
+	sys.exit(0)
+
+mail = 'jerome.sauer@uclouvain.be'
+walltime = '48:00:00'
+memory_per_task = '500mb'
+
+home_dir = os.getenv("HOME")
+scratch_dir = os.getenv("SCRATCH")
+plasim_dir = home_dir + "/PLASIM_ALGORITHM/"
+
+where = sys.argv[1]
+experiment = sys.argv[2]
+ensemble_size = int(sys.argv[3])
+
+run_folder = plasim_dir+"plasim/run"
+
+if len(sys.argv) == 5:
+	run_folder_ext = sys.argv[4]
+	run_folder += '_'+run_folder_ext
+
+if not os.path.isdir(run_folder):
+	print("PLASIM run directory missing!")
+	print("Run "+plasim_dir+"/most.x first.")
+	print("Aborting...")
+	sys.exit(1)
+
+print('Creating the directory holding the experiment ensemble...')
+experiment_folder = where+"/"+experiment+"/"
+os.system('mkdir -p '+experiment_folder)
+
+print('Creating the ensemble folders...')
+for i in range(1, ensemble_size+1):
+	member_number = str(i).rjust(2, '0')
+	ensemble_member_folder = experiment_folder+'run_'+experiment+'_'+member_number
+	os.system('cp -r '+run_folder+' '+ensemble_member_folder)
+	job_name = 'plasim_'+experiment+'_'+member_number
+	f = open(ensemble_member_folder+'/PBS_'+job_name, 'w')
+	f.write("""#!/bin/ksh
+#PBS -S /usr/bin/ksh
+#PBS -q nf
+#PBS -l EC_total_tasks=16
+#PBS -l EC_hyperthreads=1
+
+		   
+""")
+
+	f.write('#PBS -N '+job_name+'\n')
+	f.write('#PBS -o '+job_name+'.o\n')
+	f.write('#PBS -e '+job_name+'.e\n\n')
+
+	f.write('#PBS -m abe\n')
+	f.write('#PBS -M '+mail+'\n\n')
+	f.write('#PBS -l walltime='+walltime+'\n\n')
+	f.write('#PBS -l EC_memory_per_task='+memory_per_task+'\n\n')
+
+	f.write("""module switch PrgEnv-cray PrgEnv-gnu
+module load cray-snplauncher
+
+module load netcdf4
+module load netcdf
+module load cdo
+
+""")
+	f.write('cd '+ensemble_member_folder+'\n')
+	f.write('./most_plasim_run\n')
+	f.close()
+
+print('Experiment ensemble created and available in the folder:')
+print(experiment_folder)
+

cca_files/make_mpi_ensemble_experiment.py

@@ -0,0 +1,98 @@
+
+import os
+import sys
+
+# check the python version
+if float(sys.version[:3]) <= 3.6:
+	print("This script require Python 3. !")
+	print("Try:")
+	print("\n\tmodule load python3\n\n")
+	print("and restart this script again.")
+	print("Aborting...")
+	sys.exit(1)
+
+if len(sys.argv) < 4:
+	print('Bad arguments:', sys.argv)
+	print('Usage:')
+	print('\n\t python3 make_mpi_ensemble_experiment.py where experiment ensemble_size run_folder_name\n')
+	print('Arguments:\n')
+	print('\twhere :\t\t\tWhere the experiment ensemble folders are located.')
+	print('\texperiment :\t\tName of the experiment.')
+	print('\tensemble_size :\t\tSize of the ensemble.')
+	print('\trun_folder_name :\tOptional. Addtional name added to the "run_" folder name in plasim folder. Used to selected experiment executable. If not provided, assumes run folder is "run".')
+	sys.exit(0)
+
+#mail = 'jerome.sauer@uclouvain.be'
+#walltime = '48:00:00'
+#memory_per_task = '500mb'
+
+home_dir = os.getenv("HOME")
+scratch_dir = "/cfast/jsauer/"
+plasim_dir = home_dir + "/PLASIM/"
+
+where = sys.argv[1]
+experiment = sys.argv[2]
+ensemble_size = int(sys.argv[3])
+
+run_folder = plasim_dir+"plasim/run"
+
+if len(sys.argv) == 5:
+	run_folder_ext = sys.argv[4]
+	run_folder += '_'+run_folder_ext
+
+if not os.path.isdir(run_folder):
+	print("PLASIM run directory missing!")
+	print("Run "+plasim_dir+"/most.x first.")
+	print("Aborting...")
+	sys.exit(1)
+
+print('Creating the directory holding the experiment ensemble...')
+experiment_folder = where+"/"+experiment+"/"
+os.system('mkdir -p '+experiment_folder)
+
+print('Creating the ensemble folders...')
+for i in range(1, ensemble_size+1):
+	member_number = str(i).rjust(2, '0')
+	ensemble_member_folder = experiment_folder+'run_'+experiment+'_'+member_number
+	os.system('cp -r '+run_folder+' '+ensemble_member_folder)
+
+#print('Creating the ensemble folders...')
+#for i in range(1, ensemble_size+1):
+#	member_number = str(i).rjust(2, '0')
+#	ensemble_member_folder = experiment_folder+'run_'+experiment+'_'+member_number
+#	os.system('cp -r '+run_folder+' '+ensemble_member_folder)
+#	job_name = 'plasim_'+experiment+'_'+member_number
+#	f = open(ensemble_member_folder+'/PBS_'+job_name, 'w')
+#	f.write("""#!/bin/ksh
+##PBS -S /usr/bin/ksh
+##PBS -q nf
+##PBS -l EC_total_tasks=16
+##PBS -l EC_hyperthreads=1
+#
+#		   
+#""")
+#
+#	f.write('#PBS -N '+job_name+'\n')
+#	f.write('#PBS -o '+job_name+'.o\n')
+#	f.write('#PBS -e '+job_name+'.e\n\n')
+#
+#	f.write('#PBS -m abe\n')
+#	f.write('#PBS -M '+mail+'\n\n')
+#	f.write('#PBS -l walltime='+walltime+'\n\n')
+#	f.write('#PBS -l EC_memory_per_task='+memory_per_task+'\n\n')
+#
+#	f.write("""module switch PrgEnv-cray PrgEnv-gnu
+#module load cray-snplauncher
+#
+#module load netcdf4
+#module load netcdf
+#module load cdo
+#
+#""")
+#	f.write('cd '+ensemble_member_folder+'\n')
+#	f.write('./most_plasim_run\n')
+#	f.close()
+#
+#print('Experiment ensemble created and available in the folder:')
+#print(experiment_folder)
+

cca_files/restart_ensemble_experiment.py

@@ -0,0 +1,196 @@
+
+import os
+import sys
+import glob
+import subprocess
+
+# List of variable to save for output
+
+PLASIM_VAR = ['tas','hfls','hfss']
+OCEAN_VAR = ['heata','fldoa','sst']
+LSG_VAR = ['fluxhea', 'tbound']
+
+PLASIM_VAR = ['time','lat','lon'] + PLASIM_VAR
+OCEAN_VAR = ['time','lat','lon', 'ls'] + OCEAN_VAR
+LSG_VAR = ['time','lat','lon','lev', 'wet'] + LSG_VAR
+
+# check the python version
+if float(sys.version[:3]) <= 3.6:
+	print("This script require Python >= 3.6 !")
+	print("Try:")
+	print("\n\tmodule load python3\n\n")
+	print("and restart this script again.")
+	print("Aborting...")
+	sys.exit(1)
+	
+#check if ecfs utils are loaded
+user = os.getenv('USER')
+try:
+	dummy = subprocess.run(['els', 'ectmp:/'+user+'/'], check=True, stdout=subprocess.PIPE)
+except:
+	print("This script require the ECFS toolchain to be loaded !")
+	print("Try:")
+	print("\n\tmodule load ecfs\n\n")
+	print("and restart this script again.")
+	print("Aborting...")
+	sys.exit(1)
+
+#check if netcdf utils are loaded
+try:
+	dummy = subprocess.run(['ncdump'], check=True, stderr=subprocess.PIPE)
+except:
+	print("This script require NetCDF4 toolchain to be loaded !")
+	print("Try:")
+	print("\n\tmodule load netcdf4\n\n")
+	print("and restart this script again.")
+	print("Aborting...")
+	sys.exit(1)
+
+home_dir = os.getenv("HOME")
+scratch_dir = os.getenv("SCRATCH")
+perm_dir = os.getenv("PERM")
+plasim_dir = home_dir + "/PLASIM/"
+
+if len(sys.argv) < 5:
+	print('Bad arguments:', sys.argv)
+	print('Usage:')
+	print('\n\t python3 restart_ensemble_experiment.py where experiment ensemble_size number_of_years where_to_save\n')
+	print('Arguments:\n')
+	print('\twhere :\t\t\tWhere the experiment ensemble folders are located.')
+	print('\texperiment :\t\tName of the experiment.')
+	print('\tensemble_size :\t\tSize of the ensemble.')
+	print('\tnumber_of_years :\tNumber of years simulated by one run of the experiment.')
+	print('\twhere_to_save :\t\tOptional. Where to backup the previous run. If not provide, uses $PERM.')
+	sys.exit(0)
+
+where = sys.argv[1]
+basedir = where.split('/')[-1]
+experiment = sys.argv[2]
+ensemble_size = int(sys.argv[3])
+restart_year_to_save = sys.argv[4]
+try:
+	save = sys.argv[5]
+except:
+	save = perm_dir
+
+save_light = save + "/"+experiment+"_light/"
+save += "/"+experiment+"/"
+
+# Check if the experiment folder exists
+experiment_folder = where+"/"+experiment+"/"
+if not os.path.isdir(experiment_folder):
+	print("Experiment folder not found!")
+	print("Create and start the experiment "+experiment+" first.")
+	print("Aborting...")
+	sys.exit(1)
+
+# check if the experiment is still running
+queue = subprocess.run(['ssh', 'cca', '/opt/pbs/13.0.403.161593/bin/qstat', '-u', user], stdout=subprocess.PIPE, timeout=360)
+
+#if 'plasim_'+experiment[:3] in str(queue.stdout):
+#	print("Experiment still running on cca, no need to restart!")
+#	print("Aborting...")
+#	sys.exit(1)
+queue = subprocess.run(['ssh', 'ccb', '/opt/pbs/13.0.403.161593/bin/qstat', '-u', user], stdout=subprocess.PIPE, timeout=360)
+#if 'plasim_'+experiment[:3] in str(queue.stdout):
+#	print("Experiment still running on ccb, no need to restart!")
+#	print("Aborting...")
+#	sys.exit(1)
+	
+#print('Saving the results of the previous run...')
+#os.system('mkdir -p '+save)
+#os.system('mkdir -p '+save_light)
+#os.system('mkdir -p '+scratch_dir+'/tmp/'+experiment)
+#for i in range(1, ensemble_size+1):
+#	member_number = str(i).rjust(2, '0')
+#	save_experiment_folder = save + 'run_'+experiment+'_'+member_number
+#	save_experiment_folder_light = save_light + 'run_'+experiment+'_'+member_number
+#	os.system('mkdir -p ' + save_experiment_folder)
+#	# list all the past runs saved
+#	past_run_list = os.listdir(save_experiment_folder)
+#	if len(past_run_list) == 0:
+#		past_run_list = [0]
+#	else:
+#		past_run_list = list(map(int, past_run_list))
+#		past_run_list.sort()
+#	last_experiment_index = past_run_list[-1]
+#	new_experiment_index = last_experiment_index + 1
+#	# temporarily move previous experiment run to scratch
+#	if last_experiment_index > 0:
+#		to_ecfs = save_experiment_folder+'/'+str(last_experiment_index)
+#		to_light = save_experiment_folder_light+'/'+str(last_experiment_index)
+#		os.system('mkdir -p '+scratch_dir+'/tmp/'+experiment + '/run_'+experiment+'_'+member_number)
+#		os.system('rsync -a '+to_ecfs+' '+scratch_dir+'/tmp/'+experiment + '/run_'+experiment+'_'+member_number)
+#		os.system('rm -rf '+to_ecfs)
+#		os.system('rm -rf '+to_light)
+#	# make the directory to save the result of the last experiment
+#	dest = save_experiment_folder+'/'+str(new_experiment_index)
+#	os.system('mkdir -p '+dest)
+#	os.system('mkdir -p '+dest+'/restart/')
+#	os.system('mkdir -p '+dest+'/output/')
+#	# save the result of the last experiment
+#	ensemble_member_folder = experiment_folder+'run_'+experiment+'_'+member_number
+#	os.system('rsync -a '+ensemble_member_folder+'/output/* '+dest+'/output/')
+#	os.system('rsync -a '+ensemble_member_folder+'/restart/kleiswi '+dest+'/restart/')
+#	#os.system('rsync -a '+ensemble_member_folder+'/restart/*'+restart_year_to_save+' '+dest+'/restart/')
+#	os.system('rsync -a '+ensemble_member_folder+'/restart/* '+dest+'/restart/')
+#	# generating partial output files
+#	dest_light = save_experiment_folder_light+'/'+str(new_experiment_index)
+#	os.system('mkdir -p '+dest_light)
+#	os.system('mkdir -p '+dest_light+'/restart/')
+#	os.system('mkdir -p '+dest_light+'/output/')
+#	os.system('rsync -a '+ensemble_member_folder+'/output/*.txt '+dest_light+'/output/')
+#	os.system('rsync -a '+ensemble_member_folder+'/restart/kleiswi '+dest_light+'/restart/')
+#	#os.system('rsync -a '+ensemble_member_folder+'/restart/*'+restart_year_to_save+' '+dest_light+'/restart/')
+#	os.system('rsync -a '+ensemble_member_folder+'/restart/* '+dest_light+'/restart/')
+#	nc_list = glob.glob(ensemble_member_folder+'/output/*PLA*.nc')
+#	for infile in nc_list:
+#		filename = infile.split('/')[-1]
+#		os.system('nccopy -V '+','.join(PLASIM_VAR)+' '+infile+' '+dest_light+'/output/'+filename)
+#	nc_list = glob.glob(ensemble_member_folder+'/output/*OCE*.nc')
+#	for infile in nc_list:
+#		filename = infile.split('/')[-1]
+#		os.system('nccopy -V '+','.join(OCEAN_VAR)+' '+infile+' '+dest_light+'/output/'+filename)
+#	nc_list = glob.glob(ensemble_member_folder+'/output/*LSG*.nc')
+#	for infile in nc_list:
+#		filename = infile.split('/')[-1]
+#		os.system('nccopy -V '+','.join(LSG_VAR)+' '+infile+' '+dest_light+'/output/'+filename)
+#		
+#	
+#print('Backup of the previous run done !')
+#
+#if last_experiment_index > 0:
+#	print('Creating the tar archive of the experiment run number '+str(last_experiment_index))
+#	print('and saving it in ECFS temporary storage...')
+#	j = last_experiment_index
+#	yts = int(restart_year_to_save)
+#
+#        ##
+#	os.system('cd '+scratch_dir+'/tmp/ && tar -c -f '+scratch_dir+'/tmp/plasim_'+experiment+'_years_'+str(yts*(j-1)+1)+'to'+str(yts*j)+'.tar '+experiment)
+#        ##
+#
+#	#os.system('tar -c -f '+scratch_dir+'/tmp/plasim_'+experiment+'_years_'+str(yts*(j-1)+1)+'to'+str(yts*j)+'.tar '+scratch_dir+'/tmp/'+experiment)
+#	os.system('rm -rf '+scratch_dir+'/tmp/'+experiment+'/*') 
+#	queue = subprocess.run(['emkdir', '-p', 'ectmp:'+'/'+user+'/'+basedir+'/'+experiment], check=True)
+#	os.system('ecp -t '+scratch_dir+'/tmp/plasim_'+experiment+'_years_'+str(yts*(j-1)+1)+'to'+str(yts*j)+'.tar '+ 'ectmp:'+'/'+user+'/'+basedir+'/'+experiment+'/')
+#	print('Backup and move to ecfs done !')
+#
+print('Starting the ensemble runs...')
+for i in range(1, ensemble_size+1):
+	member_number = str(i).rjust(2, '0')
+	job_name = 'plasim_'+experiment+'_'+member_number
+	ensemble_member_folder = experiment_folder+'run_'+experiment+'_'+member_number
+	os.system('qsub '+ensemble_member_folder+'/PBS_'+job_name)
+
+print("Experiment '"+experiment+"' ensemble restarted.")
+print('Check the status with: qstat -u '+user)
+
+
+
+
+
+
+
+
+
+

cca_files/start_ensemble_experiment.py

@@ -0,0 +1,46 @@
+
+import os
+import sys
+
+# check the python version
+if float(sys.version[:3]) <= 3.6:
+	print("This script require Python 3.6 !")
+	print("Try:")
+	print("\n\tmodule load python3\n\n")
+	print("and restart this script again.")
+	print("Aborting...")
+	sys.exit(1)
+
+if len(sys.argv) < 4:
+	print('Bad arguments:', sys.argv)
+	print('Usage:')
+	print('\n\t python3 start_ensemble_experiment.py where experiment ensemble_size\n')
+	print('Arguments:\n')
+	print('\twhere :\t\t\tWhere the experiment ensemble folders are located.')
+	print('\texperiment :\t\tName of the experiment.')
+	print('\tensemble_size :\t\tSize of the ensemble.')
+	sys.exit(0)
+
+where = sys.argv[1]
+experiment = sys.argv[2]
+ensemble_size = int(sys.argv[3])
+experiment_folder = where+"/"+experiment+"/"
+
+if not os.path.isdir(experiment_folder):
+	print("Experiment folder not found!")
+	print("Create the experiment "+experiment+" first.")
+	print("Aborting...")
+	sys.exit(1)
+
+#os.system('qsub '+ensemble_member_folder+'/PBS_'+job_name)
+#job_name = 'plasim_'+experiment+'_'+member_number
+#print('Starting the ensemble runs...')
+for i in range(1, ensemble_size+1):
+	member_number = str(i).rjust(2, '0')
+	ensemble_member_folder = experiment_folder+'run_'+experiment+'_'+member_number
+    #print(ensemble_member_folder)
+    #os.system('cd '+ensemble_member_folder)
+
+print("Experiment '"+experiment+"' ensemble started.")
+#print('Check the status with: qstat -u '+os.getenv('USER'))
+

cca_files/start_ensemble_experiment_BACKUP.py

@@ -0,0 +1,44 @@
+
+import os
+import sys
+
+# check the python version
+if float(sys.version[:3]) <= 3.6:
+	print("This script require Python 3.6 !")
+	print("Try:")
+	print("\n\tmodule load python3\n\n")
+	print("and restart this script again.")
+	print("Aborting...")
+	sys.exit(1)
+
+if len(sys.argv) < 4:
+	print('Bad arguments:', sys.argv)
+	print('Usage:')
+	print('\n\t python3 start_ensemble_experiment.py where experiment ensemble_size\n')
+	print('Arguments:\n')
+	print('\twhere :\t\t\tWhere the experiment ensemble folders are located.')
+	print('\texperiment :\t\tName of the experiment.')
+	print('\tensemble_size :\t\tSize of the ensemble.')
+	sys.exit(0)
+
+where = sys.argv[1]
+experiment = sys.argv[2]
+ensemble_size = int(sys.argv[3])
+experiment_folder = where+"/"+experiment+"/"
+
+if not os.path.isdir(experiment_folder):
+	print("Experiment folder not found!")
+	print("Create the experiment "+experiment+" first.")
+	print("Aborting...")
+	sys.exit(1)
+
+print('Starting the ensemble runs...')
+for i in range(1, ensemble_size+1):
+	member_number = str(i).rjust(2, '0')
+	job_name = 'plasim_'+experiment+'_'+member_number
+	ensemble_member_folder = experiment_folder+'run_'+experiment+'_'+member_number
+	os.system('qsub '+ensemble_member_folder+'/PBS_'+job_name)
+
+print("Experiment '"+experiment+"' ensemble started.")
+print('Check the status with: qstat -u '+os.getenv('USER'))
+

cca_files/test.py

@@ -0,0 +1,111 @@
+
+import os
+import sys
+import glob
+import subprocess
+
+# List of variable to save for output
+
+PLASIM_VAR = ['tas','hfls','hfss']
+OCEAN_VAR = ['heata','fldoa','sst']
+LSG_VAR = ['fluxhea', 'tbound']
+
+PLASIM_VAR = ['time','lat','lon'] + PLASIM_VAR
+OCEAN_VAR = ['time','lat','lon'] + OCEAN_VAR
+LSG_VAR = ['time','lat','lon','lev'] + LSG_VAR
+
+# check the python version
+if '3.6' not in sys.version:
+	print("This script require Python 3.6 !")
+	print("Try:")
+	print("\n\tmodule load python3\n\n")
+	print("and restart this script again.")
+	print("Aborting...")
+	sys.exit(1)
+	
+#check if ecfs utils are loaded
+user = os.getenv('USER')
+try:
+	dummy = subprocess.run(['els', 'ectmp:/'+user+'/'], check=True, stdout=subprocess.PIPE)
+except:
+	print("This script require the ECFS toolchain to be loaded !")
+	print("Try:")
+	print("\n\tmodule load ecfs\n\n")
+	print("and restart this script again.")
+	print("Aborting...")
+	sys.exit(1)
+
+#check if netcdf utils are loaded
+try:
+	dummy = subprocess.run(['ncdump'], check=True, stderr=subprocess.PIPE)
+except:
+	print("This script require NetCDF4 toolchain to be loaded !")
+	print("Try:")
+	print("\n\tmodule load netcdf4\n\n")
+	print("and restart this script again.")
+	print("Aborting...")
+	sys.exit(1)
+
+home_dir = os.getenv("HOME")
+scratch_dir = os.getenv("SCRATCH")
+perm_dir = os.getenv("PERM")
+plasim_dir = home_dir + "/models/PLASIM/"
+
+if len(sys.argv) < 5:
+	print('Bad arguments:', sys.argv)
+	print('Usage:')
+	print('\n\t python3 restart_ensemble_experiment.py where experiment ensemble_size number_of_years where_to_save\n')
+	print('Arguments:\n')
+	print('\twhere :\t\t\tWhere the experiment ensemble folders are located.')
+	print('\texperiment :\t\tName of the experiment.')
+	print('\tensemble_size :\t\tSize of the ensemble.')
+	print('\tnumber_of_years :\tNumber of years simulated by one run of the experiment.')
+	print('\twhere_to_save :\t\tOptional. Where to backup the previous run. If not provide, uses $PERM.')
+	sys.exit(0)
+
+where = sys.argv[1]
+basedir = where.split('/')[-1]
+experiment = sys.argv[2]
+ensemble_size = int(sys.argv[3])
+restart_year_to_save = sys.argv[4]
+try:
+	save = sys.argv[5]
+except:
+	save = perm_dir
+
+save_light = save + "/"+experiment+"_light/"
+save += "/"+experiment+"/"
+
+# Check if the experiment folder exists
+experiment_folder = where+"/"+experiment+"/"
+if not os.path.isdir(experiment_folder):
+	print("Experiment folder not found!")
+	print("Create and start the experiment "+experiment+" first.")
+	print("Aborting...")
+	sys.exit(1)
+
+last_experiment_index = 1		
+	
+print('Backup of the previous run done !')
+
+if last_experiment_index > 0:
+	print('Creating the tar archive of the experiment run number '+str(last_experiment_index))
+	print('and saving it in ECFS temporary storage...')
+	j = last_experiment_index
+	yts = int(restart_year_to_save)
+	os.system('tar -c -f '+scratch_dir+'/tmp/plasim_'+experiment+'_years_'+str(yts*(j-1)+1)+'to'+str(yts*j)+'.tar '+scratch_dir+'/tmp/'+experiment)
+	os.system('rm -rf '+scratch_dir+'/tmp/'+experiment+'/*') 
+	queue = subprocess.run(['emkdir', '-p', 'ectmp:'+'/'+user+'/'+basedir+'/'+experiment], check=True)
+	os.system('ecp -t '+scratch_dir+'/tmp/plasim_'+experiment+'_years_'+str(yts*(j-1)+1)+'to'+str(yts*j)+'.tar '+ 'ectmp:'+'/'+user+'/'+basedir+'/'+experiment+'/')
+	print('Backup and move to ecfs done !')
+
+print('Starting the ensemble runs...')
+for i in range(1, ensemble_size+1):
+	member_number = str(i).rjust(2, '0')
+	job_name = 'plasim_'+experiment+'_'+member_number
+	ensemble_member_folder = experiment_folder+'run_'+experiment+'_'+member_number
+	os.system('qsub '+ensemble_member_folder+'/PBS_'+job_name)
+
+print("Experiment '"+experiment+"' ensemble restarted.")
+print('Check the status with: qstat -u '+user)
+

cca_files/test2.py

@@ -0,0 +1,78 @@
+
+import os
+import sys
+import glob
+import subprocess
+
+# List of variable to save for output
+
+PLASIM_VAR = ['tas','hfls','hfss']
+OCEAN_VAR = ['heata','fldoa','sst']
+LSG_VAR = ['fluxhea', 'tbound']
+
+PLASIM_VAR = ['time','lat','lon'] + PLASIM_VAR
+OCEAN_VAR = ['time','lat','lon'] + OCEAN_VAR
+LSG_VAR = ['time','lat','lon','lev'] + LSG_VAR
+
+# check the python version
+if '3.6' not in sys.version:
+	print("This script require Python 3.6 !")
+	print("Try:")
+	print("\n\tmodule load python3\n\n")
+	print("and restart this script again.")
+	print("Aborting...")
+	sys.exit(1)
+	
+#check if ecfs utils are loaded
+user = os.getenv('USER')
+try:
+	dummy = subprocess.run(['els', 'ectmp:/'+user+'/'], check=True, stdout=subprocess.PIPE)
+except:
+	print("This script require the ECFS toolchain to be loaded !")
+	print("Try:")
+	print("\n\tmodule load ecfs\n\n")
+	print("and restart this script again.")
+	print("Aborting...")
+	sys.exit(1)
+
+#check if netcdf utils are loaded
+try:
+	dummy = subprocess.run(['ncdump'], check=True, stderr=subprocess.PIPE)
+except:
+	print("This script require NetCDF4 toolchain to be loaded !")
+	print("Try:")
+	print("\n\tmodule load netcdf4\n\n")
+	print("and restart this script again.")
+	print("Aborting...")
+	sys.exit(1)
+
+home_dir = os.getenv("HOME")
+scratch_dir = os.getenv("SCRATCH")
+perm_dir = os.getenv("PERM")
+plasim_dir = home_dir + "/models/PLASIM/"
+
+if len(sys.argv) < 5:
+	print('Bad arguments:', sys.argv)
+	print('Usage:')
+	print('\n\t python3 restart_ensemble_experiment.py where experiment ensemble_size number_of_years where_to_save\n')
+	print('Arguments:\n')
+	print('\twhere :\t\t\tWhere the experiment ensemble folders are located.')
+	print('\texperiment :\t\tName of the experiment.')
+	print('\tensemble_size :\t\tSize of the ensemble.')
+	print('\tnumber_of_years :\tNumber of years simulated by one run of the experiment.')
+	print('\twhere_to_save :\t\tOptional. Where to backup the previous run. If not provide, uses $PERM.')
+	sys.exit(0)
+
+where = sys.argv[1]
+basedir = where.split('/')[-1]
+experiment = sys.argv[2]
+ensemble_size = int(sys.argv[3])
+restart_year_to_save = sys.argv[4]
+try:
+	save = sys.argv[5]
+except:
+	save = perm_dir
+
+save_light = save + "/"+experiment+"_light/"
+save += "/"+experiment+"/"
+

cca_files/update_ensemble_experiment.py

@@ -0,0 +1,61 @@
+
+import os
+import sys
+
+# check the python version
+if float(sys.version[:3]) <= 3.6:
+	print("This script require Python 3.6 !")
+	print("Try:")
+	print("\n\tmodule load python3\n\n")
+	print("and restart this script again.")
+	print("Aborting...")
+	sys.exit(1)
+
+if len(sys.argv) < 4:
+	print('Bad arguments:', sys.argv)
+	print('Usage:')
+	print('\n\t python3 update_ensemble_experiment.py where experiment ensemble_size run_folder_name\n')
+	print('Arguments:\n')
+	print('\twhere :\t\t\tWhere the experiment ensemble folders are located.')
+	print('\texperiment :\t\tName of the experiment.')
+	print('\tensemble_size :\t\tSize of the ensemble.')
+	print('\trun_folder_name :\t\tOptional. Addtional name added to the "run_" folder name in plasim folder. Used to selected experiment executable. If not provide, assumes run folder is "run".')
+	sys.exit(0)
+
+home_dir = os.getenv("HOME")
+scratch_dir = os.getenv("SCRATCH")
+plasim_dir = home_dir + "/models/PLASIM/"
+
+where = sys.argv[1]
+experiment = sys.argv[2]
+ensemble_size = int(sys.argv[3])
+
+
+run_folder = plasim_dir+"plasim/run"
+
+if len(sys.argv) == 5:
+	run_folder_ext = sys.argv[4]
+	run_folder += '_'+run_folder_ext
+
+if not os.path.isdir(run_folder):
+	print("PLASIM run directory missing!")
+	print("Run "+plasim_dir+"/most.x first.")
+	print("Aborting...")
+	sys.exit(1)
+
+# Check if the experiment folder exists
+experiment_folder = where+"/"+experiment+"/"
+if not os.path.isdir(experiment_folder):
+	print("Experiment folder not found!")
+	print("Aborting...")
+	sys.exit(1)
+
+print('Updating the ensemble folders namelists...')
+for i in range(1, ensemble_size+1):
+	member_number = str(i).rjust(2, '0')
+	ensemble_member_folder = experiment_folder+'run_'+experiment+'_'+member_number
+	os.system('cp -r '+run_folder+'/*_namelist '+ensemble_member_folder)
+
+print('Experiment ensemble updated and available in the folder:')
+print(experiment_folder)
+

cleanrepo

@@ -0,0 +1,14 @@
+#!/bin/bash
+
+rm -f F90_INTEGER
+rm -f F90_REAL
+rm -f cc_check.x
+rm -f f90check.x
+rm -f makefile
+rm -f most.x
+rm -f most_compiler
+rm -f most_compiler_mpi
+rm -f most_debug_options
+rm -f most_info.txt
+rm -f most_last_used.cfg
+rm -f most_precision_options

configure.sh

@@ -10,7 +10,7 @@ rm -f *.x most_* F90_*
 # put your favourite compiler in front if you have more than one
 # WARNING: Portland Group Compiler pgf is reported to be buggy
 
-for MOST_F90 in ifort gfortran xlf sunf90 nagfor f90 f95 pgf90 g95 af90 NO_F90
+for MOST_F90 in ftn ifort gfortran xlf sunf90 nagfor f90 f95 pgf90 g95 af90 NO_F90
 do
    `hash 2>/dev/null $MOST_F90`
    if [ $? = 0 ] ; then break ; fi
@@ -27,6 +27,13 @@ if [ $MOST_F90 != "NO_F90" ] ; then
       MOST_F90_OPTS="-O3"
       DEBUG_F90_OPTS="-g -C -fpe0 -traceback"
       echo > most_precision_options "MOST_PREC=-r8"
+   elif [ $MOST_F90 = "ftn" ] ; then
+      MOST_F90_OPTS="-O3 -ffpe-trap=invalid,zero,overflow -ffpe-summary=none -finit-real=snan"
+      DEBUG_F90_OPTS="-g -O0 -ffpe-trap=invalid,zero,overflow -ffpe-summary=none -fcheck=all -finit-real=snan"
+      echo > most_precision_options "MOST_PREC=-fdefault-real-8"
+      #MOST_F90_OPTS="-O"
+      #DEBUG_F90_OPTS="-g -C -fpe0 -traceback"
+      #echo > most_precision_options "MOST_PREC=-r8"
    elif [ $MOST_F90 = "nagfor" ] ; then
       MOST_F90_OPTS="-O -kind=byte"
       DEBUG_F90_OPTS="-g -C -fpe0 -traceback -kind=byte"
@@ -156,7 +163,7 @@ else
 fi
 
 #check for MPI-FORTRAN-90 compiler
-for MPI_F90 in openmpif90 mpif90 mpf90 mpxlf90 NO_F90
+for MPI_F90 in ftn openmpif90 mpif90 mpf90 mpxlf90 NO_F90
 do
    `hash 2>/dev/null $MPI_F90`
    if [ $? = 0 ] ; then break ; fi
@@ -166,6 +173,8 @@ if [ $MPI_F90 != "NO_F90" ] ; then
    F90_PATH=`which $MPI_F90`
    if [ $MPI_F90 = "openmpif90" ] ; then
       MPI_RUN="openmpiexec"
+   elif [ $MPI_F90 = "ftn" ] ; then
+      MPI_RUN="mpiexec"
    else
       MPI_RUN="mpiexec"
    fi

lsg/src/lsgmod.f90

@@ -8062,14 +8062,23 @@
 !     recommended values according to 14C simulations
       if (ken==11) then
         astar=1.80e-04  !1.90e-04
-       arange=1.01e-04 !1.00e-04
-       zstar=3000.
-       lambda=3.40e-03 !3.50e-03
+        arange=1.01e-04 !1.00e-04
+        zstar=3000.
+        lambda=3.40e-03 !3.50e-03
       else if (ken==22) then
-        astar=1.44e-04           !!FL 1.3E-4   !!ORI 1.44e-04
-        arange=0.78e-04          !!FL 0.3E-4   !!ORI 0.78e-04
-        zstar=2800.              !!FL 2800.    !!ORI 2800.
-        lambda=3.60e-03          !!FL 4.5E-3   !!ORI 3.60e-03
+! These are params from 0.8e-4 to 1.3e-4, as studied by R. Sciascia 
+! in her master thesis
+!        astar=1.0479e-4            !!FL 1.3E-4   !!ORI 1.44e-04
+!        arange=0.1673e-4         !!FL 0.3E-4   !!ORI 0.78e-04
+!        zstar=2500.              !!FL 2800.    !!ORI 2800.
+!        lambda=4.50e-03          !!FL 4.5E-3   !!ORI 3.60e-03
+! Parameters after the tuning in Angeloni et al. 2020
+        !astar=0.8462e-04
+        !arange=0.3011e-04
+        astar=0.8714e-04 !0.65405e-4 !0.6635e-4 !0.6824e-04 !0.7202e-04 !0.758e-04 !0.7958e-04 !0.8336e-04 !0.8714e-04
+        arange=0.2843e-04 !0.428625e-4 !0.42235e-4 !0.4098e-04 !0.3847e-04 !0.3596e-04 !0.3345e-04 !0.30942e-04 !0.2843e-04
+        zstar=2500.
+        lambda=4.5e-03
       end if
 
       do k=1,ken

most.c

@@ -325,7 +325,7 @@ struct ItemStruct NL_list[NL_MAX_ITEMS];
 int BigEndian;
 int Oce;
 int Ice;
-int Lsg = 1;
+int Lsg;
 int SimStart;
 int SimYears;
 int nreadsr;
@@ -342,7 +342,9 @@ int nac;
 int Preprocessed;
 int SAMindex;
 int ScreenHeight;
-int Expert = 0;
+int Expert = 1;
+int SuperExpert = 0;
+int PumaEnabled;
 int CatEnabled;
 int SamEnabled;
 int LsgEnabled;
@@ -409,6 +411,8 @@ int *ModeY;
 int *ModeM;
 int *ModeN;
 
+int x11flag = 1;
+
 Colormap colormap;
 
 XColor xcolor1,xcolor2;
@@ -618,15 +622,8 @@ void ChangeModel(int NewMo)
           if (SelLsg) SelLsg->no = 0;
           for (i=0 ; i < PLANETS ; ++i) SelPlanet[i]->no = 0;
        }
-       else
-       {
-       SelOce->no   = 0;
-       SelIce->no   = 0;
-       for (i=0 ; i < PLANETS ; ++i) SelPlanet[i]->no = 0;
-       }
- 
    }
-   if (Expert)
+   if (SuperExpert)
    {
       if (NewMo == CAT)
       {
@@ -1011,22 +1008,17 @@ void InitNamelist(void)
    NL_r(PLASIM,"radmod" ,"CO2"     ,360.0);
    NL_i(PLASIM,"plasim" ,"KICK"    ,  1);
    NL_i(PLASIM,"plasim" ,"MPSTEP"  ,  0);
-//FL0318   NL_i(PLASIM,"plasim" ,"NAQUA"   ,  0);
+   NL_i(PLASIM,"plasim" ,"NAQUA"   ,  0);
    NL_i(PLASIM,"plasim" ,"NDIAG"   ,  0);
-//FL0318   NL_i(PLASIM,"plasim" ,"NGUIDBG" ,  0);
-//FL0318   NL_i(PLASIM,"plasim" ,"NQSPEC"  ,  1);
-//FL0318   NL_i(PLASIM,"plasim" ,"NVEG"    ,  0);
+   NL_i(PLASIM,"plasim" ,"NGUIDBG" ,  0);
+   NL_i(PLASIM,"plasim" ,"NQSPEC"  ,  1);
+   NL_i(PLASIM,"plasim" ,"NVEG"    ,  0);
    NL_i(PLASIM,"plasim" ,"NWPD"    ,  1);
    NL_i(PLASIM,"plasim" ,"NPRINT"  ,  0);
-//FL0318   NL_i(PLASIM,"plasim" ,"NSYNC"   ,  1); 
-//FL0318   NL_i(PLASIM,"rainmod","NCLOUDS" ,  1);
-//FL0318   NL_i(PLASIM,"rainmod","NSTORAIN",  0);
-//FL0318   NL_r(PLASIM,"plasim" ,"SYNCSTR", 0.0);
-   if (Expert)
-   {
-     NL_i(PLASIM,"plasim" ,"NSYNC"   ,  0);
-     NL_r(PLASIM,"plasim" ,"SYNCSTR", 0.0);
-   }
+   NL_i(PLASIM,"plasim" ,"NSYNC"   ,  1);
+   NL_i(PLASIM,"rainmod","NCLOUDS" ,  1);
+   NL_i(PLASIM,"rainmod","NSTORAIN",  0);
+   NL_r(PLASIM,"plasim" ,"SYNCSTR", 0.0);
 
    // SAM
 
@@ -1111,7 +1103,7 @@ void NamelistSelector(int model)
       Sel->h    = FixFontHeight + 1;
       Sel->y    = yn;
       Sel->xt   = nlpos_x;
-      Sel->yt   = Sel->y + FixFont->ascent + 1;
+      Sel->yt   = Sel->y + FixFontAscent + 1;
       if (Sel->lt > ml) ml = Sel->lt;
       if (Sel->lt > ml) ml = Sel->lt;
 
@@ -1246,6 +1238,13 @@ void InitSelections(void)
    Sel->div  = Sel->iv   =  1;
    SelMod = Sel;
 
+   // Hide PUMA ?
+   if (!PumaEnabled)
+   {
+      Sel->no = 1;
+      Sel->lt = 0;
+   }
+
    // SAM
 
    Sel = NewSel(Sel);
@@ -1320,7 +1319,7 @@ void InitSelections(void)
    Sel->h    = FixFontHeight + 1;
    Sel->w    = FixFontHeight + 1;
    Sel->yt   = Sel->y + FixFontAscent + 1;
-   Sel->div  = Sel->iv   =  0;
+   Sel->div  = Sel->iv   =  1;
    Sel->no   = 1;
    SelOce    = Sel;
    Sel->piv  = &Oce;
@@ -1361,15 +1360,15 @@ void InitSelections(void)
    Sel->h    = 0;
    Sel->w    = 0;
    Sel->yt   = Sel->y + BigFontAscent + 1;
-   
 
    // Number of CPUs
 
    Sel = NewSel(Sel);
    InitNextSelection(Sel,dys,"Cores");
+
    Sel->h    = FixFontHeight + 1;
    Sel->w    = FixFontHeight + 1;
-   Sel->yt   = Sel->y + FixFont->ascent + 1;
+   Sel->yt   = Sel->y + FixFontAscent + 1;
    Sel->w    = 4 * FixFontWidth + 2;
    Sel->type = SEL_INT;
    Sel->edco = 4;
@@ -1388,7 +1387,7 @@ void InitSelections(void)
       InitNextSelection(Sel,dys,"Instances");
       Sel->h    = FixFontHeight + 1;
       Sel->w    = FixFontHeight + 1;
-      Sel->yt   = Sel->y + FixFont->ascent + 1;
+      Sel->yt   = Sel->y + FixFontAscent + 1;
       Sel->w    = 4 * FixFontWidth + 2;
       Sel->type = SEL_INT;
       Sel->edco = 4;
@@ -1407,16 +1406,16 @@ void InitSelections(void)
    Sel->type = SEL_TEXT;
    Sel->h    = 0;
    Sel->w    = 0;
-   Sel->yt   = Sel->y + BigFont->ascent + 1;
+   Sel->yt   = Sel->y + BigFontAscent + 1;
 
    // Horizontal resolution
 
-   if (Expert)
+   if (SuperExpert)
    {
       Sel = NewSel(Sel);
       InitNextSelection(Sel,dys,"Latitudes #1");
       Sel->h    = FixFontHeight + 1;
-      Sel->yt   = Sel->y + FixFont->ascent + 1;
+      Sel->yt   = Sel->y + FixFontAscent + 1;
       Sel->w    = 4 * FixFontWidth + 2;
       Sel->edco = 4;
       Sel->type = SEL_INT;
@@ -1438,7 +1437,7 @@ void InitSelections(void)
       Sel->type = SEL_CHECK;
       Sel->h    = FixFontHeight + 1;
       Sel->w    = FixFontHeight + 1;
-      Sel->yt   = Sel->y + FixFont->ascent + 1;
+      Sel->yt   = Sel->y + FixFontAscent + 1;
       Sel->div  = Sel->iv   =  1;
       SelRes    = Sel;
       DimText1  = Sel->text;
@@ -1478,7 +1477,7 @@ void InitSelections(void)
    Sel->type = SEL_TEXT;
    Sel->h    = 0;
    Sel->w    = 0;
-   Sel->yt   = Sel->y + BigFont->ascent + 1;
+   Sel->yt   = Sel->y + BigFontAscent + 1;
 
    // Global debug switch
 
@@ -1487,7 +1486,7 @@ void InitSelections(void)
    Sel->type = SEL_CHECK;
    Sel->h    = FixFontHeight + 1;
    Sel->w    = FixFontHeight + 1;
-   Sel->yt   = Sel->y + FixFont->ascent + 1;
+   Sel->yt   = Sel->y + FixFontAscent + 1;
    Sel->div  = Sel->iv   =  0;
    Sel->piv  = &ndebug;
 
@@ -1504,7 +1503,7 @@ void InitSelections(void)
    Sel->type = SEL_CHECK;
    Sel->h    = FixFontHeight + 1;
    Sel->w    = FixFontHeight + 1;
-   Sel->yt   = Sel->y + FixFont->ascent + 1;
+   Sel->yt   = Sel->y + FixFontAscent + 1;
    Sel->div  = Sel->iv   =  0;
    Sel->piv  = &noutput;
 
@@ -1523,8 +1522,7 @@ void InitSelections(void)
    Sel = NewSel(Sel);
    InitNextSelection(Sel,dyn,"Orography");
    Sel->type = SEL_CHECK;
-//FL0318   Sel->div  = Sel->iv   =  1;
-   Sel->div  = Sel->iv   =  0;
+   Sel->div  = Sel->iv   =  1;
    Sel->piv  = &noro;
    SelOro    = Sel;
 
@@ -1533,8 +1531,7 @@ void InitSelections(void)
    Sel = NewSel(Sel);
    InitNextSelection(Sel,dyn,"Annual cycle");
    Sel->type = SEL_CHECK;
-//FL0318   Sel->div  = Sel->iv   =  1;
-   Sel->div  = Sel->iv   =  0;
+   Sel->div  = Sel->iv   =  1;
    SelAnn = Sel;
 
    // Experiment
@@ -1544,7 +1541,7 @@ void InitSelections(void)
    Sel->type = SEL_TEXT;
    Sel->h    = 0;
    Sel->w    = 0;
-   Sel->yt   = Sel->y + BigFont->ascent + 1;
+   Sel->yt   = Sel->y + BigFontAscent + 1;
 
    // Simulation Start
 
@@ -1552,7 +1549,7 @@ void InitSelections(void)
    Sel = NewSel(Sel);
    InitNextSelection(Sel,dys,"Start year");
    Sel->h    = FixFontHeight + 1;
-   Sel->yt   = Sel->y + FixFont->ascent + 1;
+   Sel->yt   = Sel->y + FixFontAscent + 1;
    Sel->w    = 4 * FixFontWidth + 2;
    Sel->edco = 4;
    Sel->type = SEL_INT;
@@ -1565,7 +1562,7 @@ void InitSelections(void)
    Sel = NewSel(Sel);
    InitNextSelection(Sel,dyn,"Years to run");
    Sel->h    = FixFontHeight + 1;
-   Sel->yt   = Sel->y + FixFont->ascent + 1;
+   Sel->yt   = Sel->y + FixFontAscent + 1;
    Sel->w    = 4 * FixFontWidth + 2;
    Sel->type = SEL_INT;
    Sel->div  = Sel->iv   = 10;
@@ -1578,7 +1575,7 @@ void InitSelections(void)
    Sel = NewSel(Sel);
    InitNextSelection(Sel,dys,"Change [gpm]");
    Sel->h    = FixFontHeight + 1;
-   Sel->yt   = Sel->y + FixFont->ascent + 1;
+   Sel->yt   = Sel->y + FixFontAscent + 1;
    Sel->w    = 6 * FixFontWidth + 2;
    Sel->type = SEL_INT;
    Sel->div  = Sel->iv   = 0;
@@ -1751,6 +1748,7 @@ int WriteRunScript(int model)
    int porm;
    FILE *fp;
    char command[256];
+   char cwd[256];
    char run[256];
 
    strcpy(exec_nam2,exec_name); // Duplicate exec name
@@ -1789,24 +1787,28 @@ int WriteRunScript(int model)
    fputs("#!/bin/bash\n",fp);
    fprintf(fp,"# run-script generated by Most %s",ctime(&CurrentDate));
    fputs("EXP=MOST    # Name your experiment here\n",fp);
+   fputs("YEAR1=1\n",fp);
+   if (ngui) {
+      fputs("YEAR2=1\n",fp);
+    } else {
+      fprintf(fp,"YEAR2=%d\n",SimYears);
+    }
+   fputs("average=\"-m\"  # Comment this line if you do not want monthly averaging\n", fp);
+   getcwd(cwd, sizeof(cwd));
+   fprintf(fp,"srv2nc=%s/scripts/srv2nc\n", cwd);
    fprintf(fp,"[ $# == 1 ] && cd $1\n");
    fprintf(fp,"rm -f %s_restart\n",ShortModelName[model]);
    fputs("rm -f Abort_Message\n",fp);
-   fputs("YEAR=0\n",fp);
-   fprintf(fp,"YEARS=%d\n",SimYears);
-   fputs("srv2nc=../../tools/srv2nc\n",fp);
-   if (Multirun > 1) fprintf(fp,"INSTANCES=%d\n",Multirun);
-   if (ngui) fputs("# Remove '#' from 'while' and 'end' lines for restart loop\n",fp);
-   if (ngui) fputs("# ",fp); /* deactivate loop for GUI case */
-   fputs("while [ $YEAR -lt $YEARS ]\n",fp);
-   if (ngui) fputs("# ",fp); /* deactivate loop for GUI case */
-   fputs("do\n",fp);
-   fputs("   YEAR=`expr $YEAR + 1`\n",fp);
-   fputs("   DATANAME=`printf '%s.%03d' $EXP $YEAR`\n",fp);
-   fputs("   DIAGNAME=`printf '%s_DIAG.%03d' $EXP $YEAR`\n",fp);
-   fputs("   RESTNAME=`printf '%s_REST.%03d' $EXP $YEAR`\n",fp);
-   fputs("   OCENAME=`printf '%s_OCE.%03d.nc' $EXP $YEAR`\n",fp);
-   fputs("   ICENAME=`printf '%s_ICE.%03d.nc' $EXP $YEAR`\n",fp);
+   fprintf(fp,"mkdir -p output\n");
+   fprintf(fp,"mkdir -p restart\n");
+   fputs("for ((YEAR=$YEAR1; YEAR<=$YEAR2; YEAR++)); do \n",fp);
+   fputs("   DATANAME=`printf '%s_PLA.%04d.nc' $EXP $YEAR`\n",fp);
+   fputs("   OCENAME=`printf '%s_OCE.%04d.nc' $EXP $YEAR`\n",fp);
+   fputs("   ICENAME=`printf '%s_ICE.%04d.nc' $EXP $YEAR`\n",fp);
+   fputs("   LSGNAME=`printf '%s_LSG.%04d.nc' $EXP $YEAR`\n",fp);
+   fputs("   DIAGNAME=`printf '%s_DIAG.%04d.txt' $EXP $YEAR`\n",fp);
+   fputs("   RESTNAME=`printf '%s_REST.%04d' $EXP $YEAR`\n",fp);
+   fputs("   LSGRESTNAME=`printf '%s_LSGREST.%04d' $EXP $YEAR`\n",fp);
    if (porm < 2)
    {
       fprintf(fp,"   ./%s\n",exec_name);
@@ -1816,11 +1818,11 @@ int WriteRunScript(int model)
       if (Multirun > 1) 
       {    
          fprintf(fp,"   %s",mpirun);
-         fprintf(fp," -np 1 %s : -np 1 %s\n",exec_name,exec_nam2);
+         fprintf(fp," -np 1 %s : -np 1 ./%s\n",exec_name,exec_nam2);
       }    
       else 
       {    
-         fprintf(fp,"   %s -np %d %s\n",mpirun,porm,exec_name);
+         fprintf(fp,"   %s -np %d ./%s\n",mpirun,porm,exec_name);
       }    
    }
    fputs("   [ -e Abort_Message ] && exit 1\n",fp);
@@ -1840,16 +1842,18 @@ int WriteRunScript(int model)
    }
    else
    {
-      fprintf(fp,"   [ -e %s ] && mv %s $DATANAME\n",outp_name,outp_name);
-      fprintf(fp,"   [ -e %s ] && mv %s $DIAGNAME\n",diag_name,diag_name);
-      fprintf(fp,"   [ -e ocean_output ] && $srv2nc ocean_output $OCENAME\n");
-      fprintf(fp,"   [ -e ice_output ] && $srv2nc ice_output $ICENAME\n");
+      fprintf(fp,"   [ -e %s ] && mv %s output/$DIAGNAME\n",diag_name,diag_name);
+      fprintf(fp,"   [ -e %s ] && $srv2nc $average %s output/$DATANAME\n",outp_name,outp_name);
+      fprintf(fp,"   [ -e ocean_output ] && $srv2nc $average ocean_output output/$OCENAME\n");
+      fprintf(fp,"   [ -e ice_output ] && $srv2nc $average ice_output output/$ICENAME\n");
+      fprintf(fp,"   [ -e lsg_output ] && $srv2nc lsg_output output/$LSGNAME\n");
       fprintf(fp,"   [ -e %s_status ] && cp %s_status %s_restart\n",
               ShortModelName[model],ShortModelName[model],ShortModelName[model]);
-      fprintf(fp,"   [ -e %s_status ] && mv %s_status $RESTNAME\n",
+      fprintf(fp,"   [ -e %s_status ] && mv %s_status restart/$RESTNAME\n",
               ShortModelName[model],ShortModelName[model]);
+      fprintf(fp,"   [ -e kleiswi ] && cp kleiswi restart/kleiswi\n");
+      fprintf(fp,"   [ -e kleiin1 ] && cp kleiin1 restart/$LSGRESTNAME\n");
    }
-   if (ngui) fputs("# ",fp); /* deactivate loop for GUI case */
    fputs("done\n",fp);
    fclose(fp);
    sprintf(command,"chmod a+x %s",run);
@@ -2002,7 +2006,7 @@ int Build(int model)
    strcpy(script_backup,bld);
    strcat(script_backup,".bak");
    rename(bld,script_backup);
-
+ 
    if (model == PLASIM) WriteResmod(res);
 
    fp = fopen(bld,"w");
@@ -2029,7 +2033,6 @@ int Build(int model)
    if (Lsg)
    {
       fputs("cp -p ../../lsg/src/lsgmod.f90 .\n",fp);
-      fputs("cat lsgmod.f90\n",fp);
       putenv("OCEANCOUP=cpl");
    }
    else
@@ -2055,8 +2058,7 @@ int Build(int model)
       fputs(" ../../most_precision_options",fp);
    }
    fprintf(fp," make_%s > makefile\n",shomo);
-   fputs("source ../../modules.load\n",fp);
-   fputs("make\n",fp);
+   fputs("make -e\n",fp);
 
    fprintf(fp,"[ $? == 0 ] && cp %s.x ../bin/%s\n",shomo,exec_name);
 
@@ -2065,13 +2067,13 @@ int Build(int model)
    system(command);
 
    sprintf(message,"Building %s - wait a minute!",FullModelName[Model]);
-   BlueMessage(message);
+   if (x11flag) BlueMessage(message);
    strcat(bld," ");
    strcat(bld,shomo);
    strcat(bld,"/bld");
    if (system(bld))
    {
-      RedMessage("Error in build process");
+      if (x11flag) RedMessage("Error in build process");
       sleep(5);
       return 0; // error
    }
@@ -2369,11 +2371,11 @@ int CheckPlasimNamelist(void)
 
    // LSG works currently only with T21 PlaSim
    
-   if (Lsg)
-   {
-       Resolution = RES_T21;
-       Latitudes = ResLat[RES_T21];
-   }
+//   if (Lsg)
+//   {
+//       Resolution = RES_T21;
+//       Latitudes = ResLat[RES_T21];
+//  }
 
    // Check # of latitudes for correct values (FFT requirements)
 
@@ -2399,7 +2401,6 @@ int CheckPlasimNamelist(void)
    if (Cores <         1) Cores =         1;
 
    if (Latitudes % Cores != 0) Cores = 1;
-
    return 0; /* Success */
 }
 
@@ -3399,10 +3400,18 @@ void WriteNamelistFile(char *nl,  int instance)
    if (!strcmp(nl,"oceanmod"))
    {
       fprintf(fp," %-12s=%6d\n","NOCEAN",Oce);
-//FL0318
-      if (Lsg)
-      {
-         fprintf(fp," %-12s=%6d\n","NLSG"  ,Lsg);
+      fprintf(fp," %-12s=%6d\n","NLSG"  ,Lsg);
+      if ( Lsg == 0) {
+         fprintf(fp," %-12s=%6d\n","NHDIFF"  , 1);
+         fprintf(fp," %-12s=%6d\n","NSHDIFF" , 1);
+         fprintf(fp," %-12s= %5.1e\n","HDIFFK"  , 1.e5);
+         if ( Latitudes == 32) { // T21
+            fprintf(fp," %-12s= %5.1e\n","HDIFFK2"  , 1.e4);
+	 } else { // T42
+            fprintf(fp," %-12s= %5.1e\n","HDIFFK2"  , 3.e4);
+         }
+      } else {
+         fprintf(fp," %-12s=%6d\n","NHDIFF"  , 0);
       }
    }
    if (!strcmp(nl,"plasim"))
@@ -4382,6 +4391,38 @@ char *vcn[6] =
    "DirectColor"
 };
 
+void InitGUItext(void)
+{
+   FILE *xpp;
+
+   // Read name of MPI execute command
+
+   xpp = fopen("most_compiler_mpi","r"); // MPI installed ?
+   if (xpp)
+   {
+      fgets(Buffer,LINEMAX,xpp);
+      if (Buffer[strlen(Buffer)-1] == 10) Buffer[strlen(Buffer)-1] = 0;
+      if (Buffer[strlen(Buffer)-1] == 13) Buffer[strlen(Buffer)-1] = 0;
+      if (!strncmp(Buffer,"MPI_RUN=",8))  strcpy(mpirun,Buffer+8);
+      fclose(xpp);
+   }
+
+   InitSelections();
+   InitNamelist();
+   ChangeModel(CAT);
+   NamelistSelector(CAT);
+   ChangeModel(SAM);
+   NamelistSelector(SAM);
+   ChangeModel(PUMA);
+   NamelistSelector(PUMA);
+   ChangeModel(PLASIM);
+   NamelistSelector(PLASIM);
+   if (ReadSettings(cfg_file))
+   {
+      UpdateResolution();
+      UpdateSelections(&SelStart);
+   }
+}
 
 void InitGUI(void)
 {
@@ -4508,12 +4549,18 @@ void InitGUI(void)
       }
    }
 
-//FL0318   xpp = fopen("Beginner","r"); // Expert mode ?
-   xpp = fopen("Expert","r"); // Expert mode ?
+   xpp = fopen("Beginner","r"); // Expert mode ?
+   if (xpp)
+   {
+      Expert = 0;
+      fclose(xpp);
+   }
+
+   xpp = fopen("SuperExpert","r"); // Expert mode ?
    if (xpp)
    {
-//FL0318      Expert = 0;
       Expert = 1;
+      SuperExpert = 1;
       fclose(xpp);
    }
 
@@ -4524,6 +4571,13 @@ void InitGUI(void)
       fclose(xpp);
    }
 
+   xpp = fopen("puma","r"); // Puma enabled
+   if (xpp)
+   {
+      PumaEnabled = 1;
+      fclose(xpp);
+   }
+
    xpp = fopen("sam","r"); // Sam enabled
    if (xpp)
    {
@@ -4532,8 +4586,7 @@ void InitGUI(void)
    }
 
    xpp = fopen("lsg/src/lsgmod.f90","r"); // LSG there ?
-//FL0318   if (xpp)
-   if (xpp  && Expert)
+   if (xpp)
    {
       LsgEnabled = 1;
       fclose(xpp);
@@ -4572,14 +4625,14 @@ void InitGUI(void)
 
    InitNamelist();
    
-   ChangeModel(PLASIM);
-   NamelistSelector(PLASIM);
    ChangeModel(CAT);
    NamelistSelector(CAT);
    ChangeModel(SAM);
    NamelistSelector(SAM);
    ChangeModel(PUMA);
    NamelistSelector(PUMA);
+   ChangeModel(PLASIM);
+   NamelistSelector(PLASIM);
 
    if (ReadSettings(cfg_file))
    {
@@ -4922,12 +4975,31 @@ int main(int argc, char *argv[])
    {
            if (!strcmp(argv[ia],"-d")) Debug = 1;
       else if (!strcmp(argv[ia],"-i")) cfg_file[0] = 0;
-      else if (!strcmp(argv[ia],"-h")) ScreenHeight = atoi(argv[++ia]);
-      else strncpy(cfg_file,argv[1],sizeof(cfg_file));
+      else if (!strcmp(argv[ia],"-c")) x11flag = 0;
+      else if (!strcmp(argv[ia],"-s")) ScreenHeight = atoi(argv[++ia]);
+      else if (!strcmp(argv[ia],"-h")) {
+         puts("Usage: most.x [Options] [config_file]");
+         puts("Options:");
+ 	 puts("-d          Debug mode");
+ 	 puts("-i          Start from defaults");
+ 	 puts("-s <height> Set screen height to <height>");
+ 	 puts("-c          Console only run (no X11), build plasim");
+ 	 puts("-h          Print this help");
+         puts("");
+ 	 puts("The old configuration is read from  most_last_used.cfg");
+         puts("unless [config_file] is specified.");
+         return 0;
+      } else strncpy(cfg_file,argv[ia],sizeof(cfg_file));
    }
    CurrentDate = time(NULL);
    BigEndian = CheckEndianess();
-   InitGUI();
-   LoopGUI();
+
+   if (!x11flag) {
+      InitGUItext();
+      SaveExit();
+   } else {
+      InitGUI();
+      LoopGUI();
+   }
    return 0;
 }

plasim/src/calmod.f90

@@ -72,7 +72,6 @@
       else
          kmon = (kyday-1) / n_days_per_month
          kday = kyday - n_days_per_month * kmon
-         kmon = kmon+1
       endif
       return
       end

plasim/src/cpl.f90

@@ -1,6 +1,10 @@
       module cplmod
+      use resmod
 !
-      parameter(nxa=64,nya=32)
+      parameter(nxa=NLAT_ATM + NLAT_ATM)
+      parameter(nya=NLAT_ATM)  
+      !parameter(nxa=64,nya=32) !T42
+      !parameter(nxa=128,nya=64) !T42
       parameter(nxo=72,nyot=76,nyo=nyot-4)
       parameter(radea=6.371E6)
       parameter(nud=6)           ! output unit

plasim/src/fluxmod.f90

@@ -10,8 +10,6 @@
 !     global parameter
 !
 
-      parameter(vonkarman=0.4)    ! von karman const.
-
 !
 !     namelist parameter:
 !
@@ -28,6 +26,8 @@
       real :: vdiff_c    = 5.    !        "
       real :: vdiff_d    = 5.    !        "
 
+      real :: vonkarman  = 0.4   ! von karman const.
+
 !
 !     arrays
 !
@@ -57,7 +57,7 @@
       use fluxmod
 !
       namelist/fluxmod_nl/nvdiff,nshfl,nevap,nstress,ntsa                  &
-     &                ,zumin,vdiff_lamm,vdiff_b,vdiff_c,vdiff_d
+     &                ,zumin,vdiff_lamm,vdiff_b,vdiff_c,vdiff_d,vonkarman
 !
       if(mypid==NROOT) then
          open(11,file=fluxmod_namelist)

plasim/src/icemod.f90

@@ -1,9 +1,10 @@
       module icemod
+
       use resmod
 !
 !     version identifier (date)
 !
-      character(len=80) :: version = '15.02.2018'
+      character(len=80) :: version = '06.03.2013 by Larry'
 !
 !
 !     Parameter
@@ -38,7 +39,7 @@
       integer :: ntskin     = 1    ! compute skin temperature (0=clim.)
       integer :: ntspd      = 32   ! number of time steps per day
       integer :: noutput    = 1    ! master switch for output
-      integer :: nout       = 32   ! output each nout timesteps
+      integer :: nout       = 0    ! output each nout timesteps
       integer :: nfluko     = 0    ! switch for flux correction
                                    ! (0 = none, 1 = heat-budget, 2 = newtonian)
       integer :: ncpl_ice_ocean = 1! coupling intervall ice - ocean
@@ -49,7 +50,6 @@
       integer :: nentropy = 0      ! switch for entropy diagnostics
       integer :: ngui   = 0        ! switch for gui
       integer :: naout  = 0        ! no additional output fields 
-      integer :: nmaxd  = 1        ! switch for flx due to xmaxd (1=add)  
 !
       real :: taunc         =  0.  ! time scale for newtonian cooling
       real :: xmind         = 0.1  ! minimal ice thickness (m)
@@ -114,9 +114,9 @@
 !
 !     Climatological fields
 !
-      real :: xclsst(NHOR,0:13) = 273.15 ! climatological sst
-      real :: xclicec(NHOR,0:13)= -999.  ! climatological ice cover
-      real :: xcliced(NHOR,0:13)= -999.  ! climatological ice thickness
+      real :: xclsst(NHOR,0:13) =-999.! climatological sst
+      real :: xclicec(NHOR,0:13)=-999.! climatological ice cover
+      real :: xcliced(NHOR,0:13)=-999.! climatological ice thickness
       real :: xflxice(NHOR,0:13)= 0.! flux correction (W/m^2)
       real :: xclsst2(NHOR)   = 0.  ! climatological sst
       real :: xclssto(NHOR)   = 0.  ! climatological sst (t-1)
@@ -256,7 +256,7 @@
 !
       namelist/icemod_nl/nout,nfluko,nperpetual_ice,ntspd,nprint,nprhor &
      &               ,nentropy,nice,nsnow,ntskin,ncpl_ice_ocean,taunc   &
-     &               ,xmind,xmaxd,thicec,newsurf,naout,nmaxd
+     &               ,xmind,xmaxd,thicec,newsurf,naout
 !
 !     copy input parameter to icemod
 !
@@ -268,6 +268,10 @@
       solar_day = psolday
       deglat(:) = pdeglat(:)
 
+      if (nout <= 0) then
+          nout = ntspd
+      endif
+
 !     compute grids properties
 !
       if(mypid == NROOT) then
@@ -280,7 +284,7 @@
 !
 !     get process id
 !
-      call mpi_info(nproc,mypid)
+      call get_mpi_info(nproc,mypid)
 !
 !     print version number and read namelist
 !
@@ -309,7 +313,6 @@
       call mpbci(nprhor)
       call mpbci(nentropy)
       call mpbci(naout)
-      call mpbci(nmaxd)
       call mpbcr(taunc)
       call mpbcr(xmind)
       call mpbcr(xmaxd)
@@ -321,15 +324,9 @@
       taunc = solar_day * taunc
 
       if (nrestart == 0) then ! read start file
-
+       
          call read_ice_surface
 !
-!     limit ice thickness to xmaxd (if set)   
-!
-         if(xmaxd > 0.) then
-          xcliced(:,:)=AMIN1(xcliced(:,:),xmaxd)
-         endif
-!
 !        initialize
 !
          call iceget
@@ -371,7 +368,6 @@
          call mpgetgp('csnow'   ,csnow   ,NHOR, 1)
          call mpgetgp('xflxicea',xflxicea,NHOR, 1)
          call mpgetgp('xheata'  ,xheata  ,NHOR, 1)
-         call mpgetgp('xcfluxr' ,xcfluxr ,NHOR, 1)
          call mpgetgp('xofluxa' ,xofluxa ,NHOR, 1)
          call mpgetgp('xqmelta' ,xqmelta ,NHOR, 1)
          call mpgetgp('xcfluxa' ,xcfluxa ,NHOR, 1)
@@ -478,6 +474,7 @@
       xscflx(:)=0.
       xcfluxn(:)=0.
       xsndch(:)=0.
+      xcfluxr(:)=0.
 !
 !     copy input to icemod
 !
@@ -496,13 +493,6 @@
       xtauy(:)=ptauy(:)
       xust3(:)=pust3(:)
 !
-!     add (sub) residual flux due to limit of ice to xmaxd 
-!     (if switched on by nmaxd)
-!
-      if(nmaxd > 0) then 
-       xheat(:)=xheat(:)-xcfluxr(:)
-      endif
-!
 !     get climatology
 !
       call iceget
@@ -652,14 +642,19 @@
 !     update ximelt
 !
       zcflux(:)=0.
-      xcfluxr(:)=0.
       if(xmaxd >= 0.) then 
        where(xiced(:) > xmaxd)
         zcflux(:)=(xiced(:)-xmaxd)*zrhoilfdt
         xiced(:)=xmaxd
         xcfluxr(:)=xcfluxr(:)+zcflux(:)
         ximelt(:)=ximelt(:)+zcflux(:)
+!
+!       diagnose the lost ice as accumulated snow 
+!       (to make the budged from the atm. output) 
+!
+        xsndch(:)=xsndch(:)+zcflux(:)*1000./CRHOI/zrhoilfdt/xdt
        end where
+       call getiflx
       endif
 !
 !     depug print out if needed
@@ -940,7 +935,6 @@
       call mpputgp('csnow'   ,csnow   ,NHOR, 1)
       call mpputgp('xflxicea',xflxicea,NHOR, 1)
       call mpputgp('xheata'  ,xheata  ,NHOR, 1)
-      call mpputgp('xcfluxr' ,xcfluxr ,NHOR, 1)
       call mpputgp('xofluxa' ,xofluxa ,NHOR, 1)
       call mpputgp('xqmelta' ,xqmelta ,NHOR, 1)
       call mpputgp('xcfluxa' ,xcfluxa ,NHOR, 1)
@@ -1889,3 +1883,38 @@
       return
       end subroutine make_ice_thickness
 
+!     ==================
+!     SUBROUTINE GETIFLX
+!     ==================
+
+      subroutine getiflx
+      use icemod
+!
+      real :: zsum(2)
+      real :: zflx(NHOR) = 0.
+!
+      zrhoilfdt=CRHOI*CLFI/xdt
+!
+      where(xls(:) < 1.)
+       zflx(:)=AMAX1(0.,xmaxd-xiced(:))*zrhoilfdt
+      endwhere
+!
+      zsum(1)=SUM(xcfluxr(:)*xgw(:),MASK=(xls(:) < 1.))
+      zsum(2)=SUM(zflx(:)*xgw(:),MASK=(xls(:) < 1.))
+      call mpsumbcr(zsum,2)
+!
+      if(zsum(1) > 0. .and. zsum(2) > 0.) then
+       zfac=zsum(1)/zsum(2)
+       if(zfac <= 1.) then
+        where(xls(:) < 1.)
+         zflx(:)=zflx(:)*zfac
+        endwhere
+       endif
+       where(xls(:) < 1.)
+        xcfluxr(:)=xcfluxr(:)-zflx(:)
+        xcflux(:)=xcflux(:)-zflx(:)
+       endwhere
+      endif
+!
+      return
+      end subroutine getiflx      

plasim/src/mpimod.f90

@@ -88,9 +88,10 @@
       subroutine mpscin(k,n) ! scatter n integer
       use mpimod
 
-      integer :: k(n)
+      integer :: k(n),l(n)
 
-      call mpi_scatter(k,n,mpi_itype,k,n,mpi_itype,NROOT,myworld,mpinfo)
+      call mpi_scatter(k,n,mpi_itype,l,n,mpi_itype,NROOT,myworld,mpinfo)
+      k(:)=l(:)
 
       return
       end subroutine mpscin
@@ -102,9 +103,10 @@
       subroutine mpscrn(p,n) ! scatter n real
       use mpimod
 
-      real :: p(*)
-
-      call mpi_scatter(p,n,mpi_rtype,p,n,mpi_rtype,NROOT,myworld,mpinfo)
+      real :: p(*),l(n)
+     
+      call mpi_scatter(p,n,mpi_rtype,l,n,mpi_rtype,NROOT,myworld,mpinfo)
+      p(1:n)=l(:)
 
       return
       end subroutine mpscrn
@@ -116,9 +118,10 @@
       subroutine mpscdn(p,n) ! scatter n double precision
       use mpimod
 
-      real (kind=8) :: p(*)
+      real (kind=8) :: p(*),l(n)
 
-      call mpi_scatter(p,n,MPI_REAL8,p,n,MPI_REAL8,NROOT,myworld,mpinfo)
+      call mpi_scatter(p,n,MPI_REAL8,l,n,MPI_REAL8,NROOT,myworld,mpinfo)
+      p(1:n)=l(:)
 
       return
       end subroutine mpscdn
@@ -225,13 +228,13 @@
       subroutine mpgacs(pcs) ! gather cross sections
       use mpimod
 
-      real :: pcs(NLAT,NLEV)
+      real :: pcs(NLAT,NLEV),l(NLPP)
 
       do jlev = 1 , NLEV
-         call mpi_gather(pcs(:,jlev),NLPP,mpi_rtype                      &
+         l(:)=pcs(1:NLPP,jlev)
+         call mpi_gather(l,NLPP,mpi_rtype                      &
      &                  ,pcs(:,jlev),NLPP,mpi_rtype                      &
      &                  ,NROOT,myworld,mpinfo)
-
       enddo
       return
       end subroutine mpgacs
@@ -398,6 +401,7 @@
       integer :: itest = 0
       real    :: rtest = 0.0
       logical :: ltest = .true.
+      character (80) :: myympname
 
       if (kind(itest) == 8) mpi_itype = MPI_INTEGER8
       if (kind(rtest) == 8) mpi_rtype = MPI_REAL8
@@ -419,12 +423,12 @@
       endif
 
       allocate(ympname(npro)) ; ympname(:) = ' '
-      call mpi_get_processor_name(ympname(1),ilen,mpinfo)
 
-      call mpi_gather(ympname,80,MPI_CHARACTER,   &   
+      call mpi_get_processor_name(myympname,ilen,mpinfo)
+
+      call mpi_gather(myympname,80,MPI_CHARACTER,   &   
                       ympname,80,MPI_CHARACTER,   &   
                       NROOT,myworld,mpinfo)
-
       return
       end subroutine mpstart
 
@@ -562,10 +566,10 @@
 
 
 !     ===================
-!     SUBROUTINE MPI_INFO
+!     SUBROUTINE GET_MPI_INFO
 !     ===================
 
-      subroutine mpi_info(nprocess,npid)    ! get nproc and pid
+      subroutine get_mpi_info(nprocess,npid)    ! get nproc and pid
       use mpimod
 
       myworld=MPI_COMM_WORLD
@@ -574,7 +578,7 @@
       call mpi_comm_rank(myworld,npid,mpinfo)
 
       return
-      end subroutine mpi_info
+      end subroutine get_mpi_info
 
 
 !     ==================

plasim/src/mpimod_multi.f90

@@ -302,12 +302,12 @@
       end
 
 
-      subroutine mpi_info(nprocess,pid)    ! get nproc and pid
+      subroutine get_mpi_info(nprocess,pid)    ! get nproc and pid
       integer nprocess, pid
       nprocess = 1
       pid = 0
       return
-      end subroutine mpi_info
+      end subroutine get_mpi_info
 
 
       subroutine mpgetsp(yn,p,kdim,klev)

plasim/src/mpimod_stub.f90

@@ -217,12 +217,12 @@
       end
 
 
-      subroutine mpi_info(nprocess,pid)    ! get nproc and pid
+      subroutine get_mpi_info(nprocess,pid)    ! get nproc and pid
       integer nprocess, pid
       nprocess = 1
       pid = 0
       return
-      end subroutine mpi_info
+      end subroutine get_mpi_info
 
 
       subroutine mpgetsp(yn,p,kdim,klev)

plasim/src/oceanmod.f90

@@ -4,7 +4,7 @@
 !
 !     version identifier (date)
 !
-      character(len=80) :: version = '15.02.2018'
+      character(len=80) :: version = '13.10.2005 by Larry'
 !
 !     Parameter
 !
@@ -32,23 +32,25 @@
                                         ! (0=none,1=heat-budget,2=newtonian)
       integer :: ndiag            = 480 ! diagnostics each ndiag timesteps
       integer :: noutput          = 1   ! master switch for output: 0= no output
-      integer :: nout             = 32  ! afterburner output each nout timesteps
+      integer :: nout             = 0   ! afterburner output each nout timesteps
       integer :: nocean           = 1   ! compute ocean yes/no
       integer :: newsurf          = 0   ! update surface arrays at restart
       integer :: ntspd            = 32  ! ocean timesteps per day
       integer :: nperpetual_ocean = 0   ! perpetual climate conditions
       integer :: nprint           = 0   ! print debug information
       integer :: nprhor           = 0   ! gp to print debug information
-      integer :: nhdiff           = 0   ! switch for horizontal heat diffusion
+      integer :: nhdiff           = 1 !0   ! switch for horizontal heat diffusion
       integer :: nentropy         = 0   ! switch for entropy diagnostics
       integer :: nlsg             = 0   ! coupling flag to lsg   
-      integer :: naomod           = 320 ! atmos/ocean(lsg) ration    
+      integer :: naomod           = 320 !320 ! atmos/ocean(lsg) ration    
+      integer :: nshdiff          = 1 !0   ! flag to use separate hdiffk north and south 
 
 !
       real :: dlayer(NLEV_OCE)  = 50.   ! layer depth (m)
       real :: taunc             =  0.   ! newtonian cooling timescale (d)
       real :: vdiffkl(NLEV_OCE) = 1.E-4 ! vertikal diffusion coeff. [m**2/s]
-      real :: hdiffk(NLEV_OCE)  = 1.E3  ! horizontal diffusion coeff. [m**2/s]
+      real :: hdiffk(NLEV_OCE)  = 1.E5  ! horizontal diffusion coeff. [m**2/s]
+      real :: hdiffk2(NLEV_OCE) = 3.E4 !1.E4  ! horizontal diffusion coeff. SOUTH [m**2/s]
 !
 !     global integer
 !
@@ -91,7 +93,7 @@
       real :: ydsst(NHOR)  = 0.         ! heat flux from vdiff (w/m2)
       real :: yqhd(NHOR)   = 0.         ! heat flux from hdiff (w/m2)
 
-      real :: yclsst(NHOR,0:13)= 273.15 ! climatological sst (K)
+      real :: yclsst(NHOR,0:13)         ! climatological sst (K)
       real :: yfsst(NHOR,0:13) = 0.     ! flux corr. sst (W/m**2)
 
       real :: yclsst2(NHOR) = 0.        ! climatological sst (K)
@@ -144,11 +146,11 @@
 !
       namelist/oceanmod_nl/ndiag,nout,nfluko,ntspd,nocean,nprint,nprhor    &
     &                  ,nperpetual_ocean,naomod,nlsg,taunc,dlayer       &
-    &                  ,vdiffkl,newsurf,hdiffk,nentropy,nhdiff
+    &                  ,vdiffkl,newsurf,hdiffk,hdiffk2,nshdiff,nentropy,nhdiff
 !
 !     get process id
 !
-      call mpi_info(nproc,mypid)
+      call get_mpi_info(nproc,mypid)
 !
 !     compute grids properties
 !
@@ -186,6 +188,10 @@
       ngui      = kgui
       ntspd     = ktspd
       solar_day = psolday
+
+      if (nout <= 0) then
+          nout = ntspd
+      endif
 !
 !     read and print namelist and distribute it
 !
@@ -215,10 +221,12 @@
       call mpbci(naomod)
       call mpbci(nentropy)
       call mpbci(nhdiff)
+      call mpbci(nshdiff)
       call mpbcr(taunc)
       call mpbcr(solar_day)
       call mpbcrn(vdiffkl,NLEV_OCE)
       call mpbcrn(hdiffk,NLEV_OCE)
+      call mpbcrn(hdiffk2,NLEV_OCE)
       call mpbcrn(dlayer,NLEV_OCE)
 !
       do jlev=1,NLEV_OCE
@@ -1307,7 +1315,8 @@
       call mpgagp(zls,yls,1)
 !
       do jlev=1,NLEV_OCE
-       zfac=hdiffk(jlev)/plarad/plarad
+       zfac=hdiffk(jlev)/plarad/plarad   ! Diffusivity Global or NORTH
+       zfac2=hdiffk2(jlev)/plarad/plarad ! Diffusivity SOUTH
 !
        call mpgagp(zt,zsst(1,jlev),1)
 !
@@ -1351,6 +1360,7 @@
          zdty(:,0)=0.
          zdty(:,NLAT)=0.
 !
+         if(nshdiff==0) then
          do jlat=1,NLAT
           jlam=jlat-1
           do jlon=1,NLON
@@ -1364,6 +1374,36 @@
            endif
           enddo
          enddo
+         else
+! NORTH
+         do jlat=1,NLAT/2
+          jlam=jlat-1
+          do jlon=1,NLON
+           jlom=jlon-1
+           if(zls(jlon,jlat) < 1.) then
+            zdtdt(jlon,jlat)=((zdtx(jlon,jlat)-zdtx(jlom,jlat))         &
+     &                        /dlam/cphi(jlat)/cphi(jlat)               &
+     &                       +(zdty(jlon,jlat)-zdty(jlon,jlam))         &
+     &                        /dmue(jlat))                              &
+     &                      *zfac
+           endif
+          enddo
+         enddo
+! SOUTH
+         do jlat=NLAT/2+1,NLAT
+          jlam=jlat-1
+          do jlon=1,NLON
+           jlom=jlon-1
+           if(zls(jlon,jlat) < 1.) then
+            zdtdt(jlon,jlat)=((zdtx(jlon,jlat)-zdtx(jlom,jlat))         &
+     &                        /dlam/cphi(jlat)/cphi(jlat)               &
+     &                       +(zdty(jlon,jlat)-zdty(jlon,jlam))         &
+     &                        /dmue(jlat))                              &
+     &                      *zfac2
+           endif
+          enddo
+         enddo
+         endif
 !
          if(nentropy > 0) then
           where(zls(:,:) < 1.)

plasim/src/plasim.f90

@@ -40,9 +40,8 @@
 !     tdiss days on every variable at the truncation wavenumber.
 !
 !     UPDATE VERSION IDENTIFIER AFTER EACH CODE CHANGE!
-      plasimwww="https://www.mi.uni-hamburg.de/en/arbeitsgruppen"       &
-     &        //"/theoretische-meteorologie/modelle/plasim.html"
-      plasimversion = "15-Mar-2018"
+
+plasimversion = "https://github.com/Edilbert/PLASIM/ : 15-Dec-2015"
 
       call mpstart(-1)       ! -1: Start MPI   >=0 arg = MPI_COMM_WORLD
       call setfilenames
@@ -145,16 +144,7 @@
          call cpu_time(tmstart)
          write(nud,'(54("*"))')
          write(nud,'("* ",17X,"PLANET SIMULATOR",17X," *")')
-         write(nud,'("* ",50X," *")')
-         write(nud,'("* ",a50," *")') "Version: "//plasimversion
-         write(nud,'("* ",a50," *")') "from "//plasimwww(1:45)
-         if(len_trim(plasimwww) > 45) then
-          write(nud,'("* ",a50," *")') "     "//plasimwww(46:90)
-         endif
-         if(len_trim(plasimwww) > 90) then
-          write(nud,'("* ",a50," *")') "     "//plasimwww(91:)
-         endif
-         write(nud,'("* ",50X," *")')
+         write(nud,'("* ",a50," *")') plasimversion
          write(nud,'(54("*"))')
          if (mrnum > 1) then
             write(nud,'("* Instance ",i3," of ",i3,32x,"*")') &

plasim/src/plasim_dummy.f90

@@ -187,7 +187,7 @@
 !
 !     get process id
 !
-      call mpi_info(nproc,mypid)
+      call get_mpi_info(nproc,mypid)
 !
 !     read and print namelist and distribute it
 !

plasim/src/plasimmod.f90

@@ -152,7 +152,7 @@
       integer :: noutput  =  1  ! master switch for output: 0=no output
       integer :: nafter   =  0  ! write data interval: 0 = once per day
       integer :: naqua    =  0  ! 1: switch to aqua planet mode
-      integer :: nveg     =  0  ! 1: run vegetation module
+      integer :: nveg     =  1  ! 1: run vegetation module
       integer :: ncoeff   =  0  ! number of modes to print
       integer :: ndiag    =  0  ! write diagnostics interval 0 = every 10th. day
       integer :: ngui     =  0  ! 1: run with GUI
@@ -493,13 +493,12 @@
       real    :: syncstr  =  0.0 ! Coupling strength (0 .. 1)
       real    :: synctime =  0.0 ! Coupling time [days]
 
-!     ************************************************************
-!     * download page and version identifier (set in plasim.f90) *
-!     * will be printed in plasim sub *prolog*                   *
-!     ************************************************************
+!     **********************************************
+!     * version identifier (to be set in puma.f90) *
+!     * will be printed in puma sub *prolog*       *
+!     **********************************************
 
-      character(len=135):: plasimwww
-      character(len=41) :: plasimversion
+      character(len=80) :: plasimversion
 
 !     ******************************
 !     * Planet dependent variables *

plasim/src/restartmod.f90

@@ -1,9 +1,4 @@
       module restartmod
-!
-!     version identifier (date)
-!
-      character(len=80) :: version = '15.02.2018'
-
       integer, parameter :: nresdim  = 200     ! Max number of records
       integer, parameter :: nreaunit =  33     ! FORTRAN unit for reading
       integer, parameter :: nwriunit =  34     ! FORTRAN unit for writing
@@ -26,10 +21,6 @@
       character (len=*)  :: yrfile
       character (len=16) :: yn ! variable name
 
-      write(nud,'(/," ***********************************************")')
-      write(nud,'(" * RESTARTMOD ",a32," *")') trim(version)
-      write(nud,'(" ***********************************************")')
-
       inquire(file=yrfile,exist=lrestart)
       if (lrestart) then
          open(nreaunit,file=yrfile,form='unformatted')
@@ -137,9 +128,9 @@
          endif
       enddo
       if (nexcheck == 1) then
-         write(nud,*)'*** WARNING in get_restart_array ***'
+         write(nud,*)'*** Error in get_restart_array ***'
          write(nud,*)'Requested array {',yn,'} was not found'
-         write(nud,*)'Default values will be used'
+         stop
       endif
       return
       end subroutine get_restart_array

plasim/src/surfmod.f90

@@ -1,7 +1,7 @@
       module surfmod
       use pumamod
 
-      character(len=80) :: version = '15.02.2018'
+      character(len=80) :: version = '12.06.2014 by Edi'
 
       integer, parameter :: nsurdim  = 100           ! max # of variables
       integer, parameter :: nsurunit =  35           ! read unit
@@ -346,14 +346,14 @@
        open(11,file=surfmod_namelist)
        read(11,surfmod_nl)
        close(11)
-       write(nud,'(/," ***********************************************")')
-       write(nud,'(" * SURFMOD ",a35," *")') trim(version)
-       write(nud,'(" ***********************************************")')
+       write(nud,'(/,"***********************************************")')
+       write(nud,'("* SURFMOD ",a35," *")') trim(version)
+       write(nud,'("***********************************************")')
        if (naqua /= 0) then
-       write(nud,'(" * AQUA planet mode - ignoring land data       *")')
+       write(nud,'("* AQUA planet mode - ignoring land data       *")')
        endif
-       write(nud,'(" * Namelist SURFMOD_NL from <surfmod_namelist> *")')
-       write(nud,'(" ***********************************************")')
+       write(nud,'("* Namelist SURFMOD_NL from <surfmod_namelist> *")')
+       write(nud,'("***********************************************")')
        write(nud,surfmod_nl)
       endif
 
@@ -368,13 +368,6 @@
          so(:)   = 0.0         ! spectral  orography
          sp(:)   = 0.0         ! spectral  pressure
          spm(:)  = 0.0         ! spectral  pressure scattered
-!
-!    add noise
-!
-         if (mypid == NROOT) then
-          call noise
-         endif
-         call mpscsp(sp,spm,1)   
       endif
 
       if (nrestart == 0 .and. naqua == 0) then ! need to read start data

postprocessor/burn7.cpp

@@ -5542,6 +5542,8 @@ void InitAll(void)
    All[183].Init("tso"  ,"climate_deep_soil_temperature"   ,"K"        ,1); // Not standard
    All[184].Init("wsoi" ,"climate_deep_soil_wetness"       ,"1"        ,1);
    All[199].Init("vegc" ,"vegetation_cover"                ,"1"        ,1); // Not standard
+   All[201].Init("tasmax" ,"maximum_air_temperature_2m"     ,"K"    ,1); // Not standard
+   All[202].Init("tasmin" ,"minimum_air_temperature_2m"     ,"K"    ,1); // Not standard
    All[203].Init("rsut" ,"toa_outgoing_shortwave_flux"     ,"W m-2"    ,1); // Not standard
    All[204].Init("ssru" ,"surface_solar_radiation_upward"  ,"W m-2"    ,1); // Not standard
    All[205].Init("stru" ,"surface_thermal_radiation_upward","W m-2"    ,1); // Not standard

postprocessor/makefile_sav

@@ -0,0 +1,3 @@
+burn7.x:	burn7.cpp
+	#c++ -O2 -o burn7.x burn7.cpp -lm -lnetcdf_c++
+	c++ -O2 -o burn7.x burn7.cpp -lm -lnetcdf 

scripts/README_burn.txt

@@ -0,0 +1,92 @@
+*burn.sh v1.0*
+--------------
+Afterburner cmorizer
+Extracts variables from PlaSim output using CMOR conventions.
+Adjusts variable names, signs, units and computes derived variables.
+
+Usage: burn.sh [options] infile outfile var
+Options:
+        -h                Prints this help
+        -s                Extract 3D variables on sigma levels
+        -a                Afterburner variable (do not cmorize)
+        -p lev1,lev2,...  Specify pressure levels (in hPa)
+
+Default pressure levels are  20,30,50,70,100,150,200,250,300,400,500,600,700,850,925,1000 hPa
+
+Variables which can be extracted:
+-------------------------------------------------------------------------------------
+| 2D CMOR2 variables:                                                               |
+-------------------------------------------------------------------------------------
+|     Name |                                Long Name |      Units | From org. vars |
+-------------------------------------------------------------------------------------
+|       ps |                     Surface Air Pressure |         Pa |             pl |
+|      psl |                       Sea Level Pressure |         Pa |            psl |
+|      tas |             Near-Surface Air Temperature |          K |            tas |
+|   tasmin |Daily Minimum Near-Surface Air Temperature|          K |         tasmin |
+|   tasmax |Daily Maximum Near-Surface Air Temperature|          K |         tasmax |
+|       ts |                      Surface Temperature |          K |             ts |
+|      uas |               Eastward Near-Surface Wind |      m s-1 |            uas |
+|      vas |              Northward Near-Surface Wind |      m s-1 |            vas |
+|  evspsbl |                              Evaporation | kg m-2 s-1 |           evap |
+|       pr |                            Precipitation | kg m-2 s-1 |             pr |
+|      prc |                 Convective Precipitation | kg m-2 s-1 |            prc |
+|      prl |                Large Scale Precipitation | kg m-2 s-1 |            prl |
+|     prsn |                            Snowfall Flux | kg m-2 s-1 |           prsn |
+|    mrros |                           Surface Runoff | kg m-2 s-1 |           mrro |
+|     mrso |              Total Soil Moisture Content | kg m-2 s-1 |           mrso |
+|      tsl |                      Temperature of Soil |          K |            tso |
+|     rsdt |         TOA Incident Shortwave Radiation |      W m-2 |       rsut rst |
+|     rsut |         TOA Outgoing Shortwave Radiation |      W m-2 |           rsut |
+|     rlut |          TOA Outgoing Longwave Radiation |      W m-2 |           rlut |
+|     rsds |  Surface Downwelling Shortwave Radiation |      W m-2 |       ssru rss |
+|     rsus |    Surface Upwelling Shortwave Radiation |      W m-2 |           ssru |
+|     rlds |   Surface Downwelling Longwave Radiation |      W m-2 |       stru rls |
+|     rlus |     Surface Upwelling Longwave Radiation |      W m-2 |           stru |
+|     hfls |          Surface Upward Latent Heat Flux |      W m-2 |           hfls |
+|     hfss |        Surface Upward Sensible Heat Flux |      W m-2 |           hfss |
+|      clt |                     Total Cloud Fraction |          % |            clt |
+|      prw |                         Water Vapor Path |     kg m-2 |            prw |
+|     huss |           Near-Surface Specific Humidity |          1 |            hus |
+|     hurs |           Near-Surface Specific Humidity |          1 |       tas td2m |
+|     orog |                         Surface Altitude |          m |             sg |
+|      snd |                               Snow Depth |          m |            snd |
+|      snm |                        Surface Snow Melt | kg m-2 s-1 |            snm |
+|     tauu |    Surface Downward Eastward Wind Stress |         Pa |           tauu |
+|     tauv |   Surface Downward Northward Wind Stress |         Pa |           tauv |
+|      sit |                        Sea Ice Thickness |          m |            sit |
+|      sic |                            Sea Ice Cover |          % |            sic |
+-------------------------------------------------------------------------------------
+| 3D CMOR2 variables:                                                               |
+-------------------------------------------------------------------------------------
+|     Name |                                Long Name |      Units | From org. vars |
+-------------------------------------------------------------------------------------
+|      hus |                        Specific Humidity |          1 |            hus |
+|      hur |                        Relative Humidity |          % |            hur |
+|       zg |                      Geopotential Height |          1 |             zg |
+|       ta |                          Air Temperature |          K |             ta |
+|       ua |                            Eastward Wind |      m s-1 |             ua |
+|       va |                           Northward Wind |      m s-1 |             va |
+|       wa |                              Upward Wind |      m s-1 |             wa |
+|      wap |                           omega (=dp/dt) |     Pa s-1 |            wap |
+|       cl |                      Cloud Area Fraction |          % |             cl |
+|      clw |      Mass Fraction of Cloud Liquid Water |          1 |            clw |
+-------------------------------------------------------------------------------------
+| Not standard variables:                                                           |
+-------------------------------------------------------------------------------------
+|     Name |                                Long Name |      Units | From org. vars |
+-------------------------------------------------------------------------------------
+|     td2m |        Near-Surface Dewpoint Temperature |          K |           td2m |
+|      rst |              TOA Net Shortwave Radiation |      W m-2 |            rst |
+|      rls |           Surface Net Longwave Radiation |      W m-2 |            rls |
+|      alb |                           Surface Albedo |          1 |            alb |
+|      prl |                Large Scale Precipitation | kg m-2 s-1 |            prl |
+|     sndc |                Surface Snow Depth Change |      m s-1 |           sndc |
+|      stf |                          Stream Function |          m |            stf |
+|      lsm |                            Land Sea Mask |          1 |            lsm |
+|      mld |              Ocean Mixed Layer Thickness |          m |            mld |
+-------------------------------------------------------------------------------------
+Variables not in this list are extracted directly using afterburner with no change.
+Examples: zeta,d,bld etc.
+Note: uas and vas are not available in PlaSim output, derived from ua and va at sigma=1
+Note: psl (sea level pressure) is not avaiable in PlaSim output,
+      derived from huss (hus at sigma=1), tas and ps using the hypsometric formula

scripts/burn.sh

@@ -0,0 +1,213 @@
+#!/bin/bash
+## burn.sh
+## Afterburner cmorizer
+## (2018) by Jost von Hardenberg <j.vonhardenberg AT isac.cnr.it>
+VERSION=v1.0
+
+#set -ex
+BURNER=/work/users/jost/plasim/postprocessor/burn7.x
+BURNER=/Users/jost/plasim/plasimnew/postprocessor/burn7.x
+PLEVS="20,30,50,70,100,150,200,250,300,400,500,600,700,850,925,1000"
+VARS3D="wap wa zeta stf psi d zg hur ta ua va hus clw cl spd psl"
+GACC=9.80665
+HUM_A=6.116441
+HUM_m=7.591386 
+HUM_Tn=240.7263
+
+function makevar {
+	outvar=$1
+	no3d=false;
+	case $outvar in
+		"ps")      invar="ps";   standard_name="surface_air_pressure" ;  long_name="Surface Air Pressure"; units="Pa";         expr="ps=ps*100";;
+		"psl")      invar="psl";   standard_name="air_pressure_at_sea_level" ;  long_name="Sea Level Pressure"; units="Pa";     expr="psl=psl*100";; 
+		"tas")     invar="tas";  standard_name="air_temperature" ;long_name="Near-Surface Air Temperature"; units="K"; expr="tas=tas";;
+		"tasmin")     invar="tasmin";  standard_name="air_temperature" ;long_name="Daily Minimum Near-Surface Air Temperature"; units="K"; expr="tasmin=tasmin";;
+		"tasmax")     invar="tasmax";  standard_name="air_temperature" ;long_name="Daily Maximum Near-Surface Air Temperature"; units="K"; expr="tasmax=tasmax";;
+		"ts")     invar="ts";    standard_name="surface_temperature" ;long_name="Surface Temperature"; units="K"; expr="ts=ts";;
+		"uas")     invar="ua";    standard_name="eastward_wind" ;long_name="Eastward Near-Surface Wind"; units="m s-1"; expr="uas=ua"; no3d=true;;
+		"vas")     invar="va";    standard_name="northward_wind" ;long_name="Northward Near-Surface Wind"; units="m s-1"; expr="vas=va"; no3d=true;;
+		"evspsbl") invar="evap"; standard_name="water_evaporation_flux" ;long_name="Evaporation";          units="kg m-2 s-1"; expr="evspsbl=-evap*1000";;
+		"pr")     invar="pr";    standard_name="precipitation_flux" ;long_name="Precipitation"; units="kg m-2 s-1"; expr="pr=pr*1000";;
+		"prc")     invar="prc";    standard_name="convective_precipitation_flux" ;long_name="Convective Precipitation"; units="kg m-2 s-1"; expr="prc=prc*1000";;
+		"prsn")     invar="prsn";    standard_name="snowfall_flux" ;long_name="Snowfall Flux"; units="kg m-2 s-1"; expr="prsn=prsn*1000";;
+		"mrros")     invar="mrro";    standard_name="surface_runoff_flux" ;long_name="Surface Runoff"; units="kg m-2 s-1"; expr="mrros=mrro*1000";;
+		"mrso")     invar="mrso";    standard_name="soil_moisture_content" ;long_name="Total Soil Moisture Content"; units="kg m-2 s-1"; expr="mrso=mrso*1000";;
+		"tsl")     invar="tso";    standard_name="soil_temperature" ;long_name="Temperature of Soil"; units="K"; expr="tsl=tso";;
+		"rsdt")     invar="rsut rst";    standard_name="toa_incoming_shortwave_flux" ;long_name="TOA Incident Shortwave Radiation"; units="W m-2"; expr="rsdt=rst-rsut";;
+		"rsut")     invar="rsut";    standard_name="toa_outgoing_shortwave_flux" ;long_name="TOA Outgoing Shortwave Radiation"; units="W m-2"; expr="rsut=-rsut";;
+		"rlut")     invar="rlut";    standard_name="toa_outgoing_longwave_flux" ;long_name="TOA Outgoing Longwave Radiation"; units="W m-2"; expr="rlut=-rlut";;
+		"rsds")     invar="ssru rss";    standard_name="surface_downwelling_shortwave_flux_in_air" ;long_name="Surface Downwelling Shortwave Radiation"; units="W m-2"; expr="rsds=rss-ssru";;
+		"rsus")     invar="ssru";    standard_name="surface_upwelling_shortwave_flux_in_air" ;long_name="Surface Upwelling Shortwave Radiation"; units="W m-2"; expr="rsus=-ssru";;
+		"rlds")     invar="stru rls";    standard_name="surface_downwelling_longwave_flux_in_air" ;long_name="Surface Downwelling Longwave Radiation"; units="W m-2"; expr="rlds=rls-stru";;
+		"rlus")     invar="stru";    standard_name="surface_upwelling_longwave_flux_in_air" ;long_name="Surface Upwelling Longwave Radiation"; units="W m-2"; expr="rlus=-stru";;
+		"hfls")     invar="hfls";    standard_name="surface_upward_latent_heat_flux" ;long_name="Surface Upward Latent Heat Flux"; units="W m-2"; expr="hfls=-hfls";;
+		"hfss")     invar="hfss";    standard_name="surface_upward_sensible_heat_flux" ;long_name="Surface Upward Sensible Heat Flux"; units="W m-2"; expr="hfss=-hfss";;
+		"clt")     invar="clt";  standard_name="cloud_area_fraction" ;   long_name="Total Cloud Fraction"; units="%"; expr="clt=clt*100";;
+		"prw")     invar="prw";    standard_name="atmosphere_water_vapor_content" ;long_name="Water Vapor Path"; units="kg m-2"; expr="prw=prw";;
+		"huss")     invar="hus";    standard_name="specific_humidity" ;long_name="Near-Surface Specific Humidity"; units="1"; expr="huss=hus"; no3d=true;;
+		"hurs")     invar="tas td2m";    standard_name="specific_humidity" ;long_name="Near-Surface Specific Humidity"; units="1"; expr="hurs=100*10^($HUM_m*( td2m/(td2m + $HUM_Tn )-tas/(tas + $HUM_Tn )   ))	"; no3d=true;;
+		"clw")     invar="clw";    standard_name="mass_fraction_of_cloud_liquid_water_in_air" ;long_name="Mass Fraction of Cloud Liquid Water"; units="1"; expr="clw=clw";;
+		"hus")     invar="hus";    standard_name="specific_humidity" ;long_name="Specific Humidity"; units="1"; expr="hus=hus";;
+		"hur")     invar="hur";    standard_name="relative_humidity" ;long_name="Relative Humidity"; units="%"; expr="hur=hur";;
+		"zg")     invar="zg";    standard_name="geopotential_height" ;long_name="Geopotential Height"; units="1"; expr="zg=zg";;
+		"ta")     invar="ta";    standard_name="air_temperature" ;long_name="Air Temperature"; units="K"; expr="ta=ta";;
+		"ua")     invar="ua";    standard_name="eastward_wind" ;long_name="Eastward Wind"; units="m s-1"; expr="ua=ua";;
+		"va")     invar="va";    standard_name="northward_wind" ;long_name="Northward Wind"; units="m s-1"; expr="va=va";;
+		"wa")     invar="wa";    standard_name="upward_wind" ;long_name="Upward Wind"; units="m s-1"; expr="wa=wa";;
+		"wap")     invar="wap";    standard_name="lagrangian_tendency_of_air_pressure" ;long_name="omega (=dp/dt)"; units="Pa s-1"; expr="wap=wap";;
+		"cl")     invar="cl";    standard_name="cloud_area_fraction_in_atmosphere_layer" ;long_name="Cloud Area Fraction"; units="%"; expr="cl=cl*100";;
+		"orog")     invar="sg";    standard_name="surface_altitude" ;long_name="Surface Altitude"; units="m"; expr="orog=sg/$GACC";;
+		"snd")     invar="snd";    standard_name="surface_snow_thickness" ;long_name="Snow Depth"; units="m"; expr="snd=snd";;
+		"snm")     invar="snm";    standard_name="surface_snow_melt_flux" ;long_name="Surface Snow Melt"; units="kg m-2 s-1"; expr="snm=snm*1000";;
+		"tauu")     invar="tauu";    standard_name="surface_downward_eastward_stress" ;long_name="Surface Downward Eastward Wind Stress"; units="Pa"; expr="tauu=tauu";;
+		"tauv")     invar="tauv";    standard_name="surface_downward_northward_stress" ;long_name="Surface Downward Northward Wind Stress"; units="Pa"; expr="tauv=tauv";;
+
+		"td2m")     invar="td2m";    standard_name="dewpoint_temperature" ;long_name="Near-Surface Dewpoint Temperature"; units="K"; expr="td2m=td2m";;
+		"rst")     invar="rst";    standard_name="toa_net_shortwave_flux" ;long_name="TOA Net Shortwave Radiation"; units="W m-2"; expr="rst=rst";;
+		"rss")     invar="rss";    standard_name="surface_net_shortwave_flux_in_air" ;long_name="Surface Net Shortwave Radiation"; units="W m-2"; expr="rss=rss";;
+		"rls")     invar="rls";    standard_name="surface_net_longwave_flux_in_air" ;long_name="Surface Net Longwave Radiation"; units="W m-2"; expr="rls=rls";;
+		"alb")     invar="alb";    standard_name="surface_albedo" ;long_name="Surface Albedo"; units="1"; expr="alb=alb";;
+		"prl")     invar="prl";    standard_name="large_scale_precipitation_flux" ;long_name="Large Scale Precipitation"; units="kg m-2 s-1"; expr="prl=prl*1000";;
+		"sndc")     invar="sndc";    standard_name="surface_snow_depth_change" ;long_name="Surface Snow Depth Change"; units="m s-1"; expr="sndc=sndc";;
+		"stf")     invar="stf";    standard_name="stream_function" ;long_name="Stream Function"; units="m"; expr="stf=stf";;
+		"lsm")     invar="lsm";    standard_name="land_sea_mask" ;long_name="Land Sea Mask"; units="1"; expr="lsm=lsm";;
+
+		"sit")     invar="sit";    standard_name="sea_ice_thickness" ;long_name="Sea Ice Thickness"; units="m"; expr="sit=sit";;
+		"sic")     invar="sic";    standard_name="sea_ice_cover" ;long_name="Sea Ice Cover"; units="%"; expr="sic=sic*100";;
+		"mld")     invar="mld";    standard_name="ocean_mixed_layer_thickness" ;long_name="Ocean Mixed Layer Thickness"; units="m"; expr="mld=mld";;
+		*)         invar=$outvar ;    expr="none";;
+	esac
+}
+
+function help {
+HELPVARS2D="ps psl tas tasmin tasmax ts uas vas evspsbl pr prc prl prsn mrros mrso tsl rsdt rsut rlut rsds rsus rlds rlus hfls hfss clt prw huss hurs orog snd snm tauu tauv sit sic"
+HELPVARS3D="hus hur zg ta ua va wa wap cl clw"
+HELPEXTRA="td2m rst rls alb prl sndc stf lsm mld"
+
+echo "burn.sh $VERSION"
+echo "Afterburner cmorizer"
+echo "Extracts variables from PlaSim output using CMOR conventions."
+echo "Adjusts variable names, signs, units and computes derived variables."
+echo
+echo Usage: burn.sh [options] infile outfile var 
+echo Options:
+echo "        -h                Prints this help"
+echo "        -s                Extract 3D variables on sigma levels"
+echo "        -a                Afterburner variable (do not cmorize)"
+echo "        -p lev1,lev2,...  Specify pressure levels (in hPa)"
+echo
+echo "Default pressure levels are " $PLEVS "hPa"
+echo
+echo Variables which can be extracted:
+echo -------------------------------------------------------------------------------------
+echo "| 2D CMOR2 variables:                                                               |"
+echo -------------------------------------------------------------------------------------
+printf "| % 8s | % 40s | % 10s | % 14s |\n" Name "Long Name" Units "From org. vars"
+echo -------------------------------------------------------------------------------------
+for v in $HELPVARS2D; do
+	makevar $v
+	printf "| % 8s | % 40s | % 10s | % 14s |\n" $v "$long_name" "$units" "$invar" 
+done
+echo -------------------------------------------------------------------------------------
+echo "| 3D CMOR2 variables:                                                               |"
+echo -------------------------------------------------------------------------------------
+printf "| % 8s | % 40s | % 10s | % 14s |\n" Name "Long Name" Units "From org. vars"
+echo -------------------------------------------------------------------------------------
+for v in $HELPVARS3D; do
+	makevar $v
+	printf "| % 8s | % 40s | % 10s | % 14s |\n" $v "$long_name" "$units" "$invar" 
+done
+echo -------------------------------------------------------------------------------------
+echo "| Not standard variables:                                                           |"
+echo -------------------------------------------------------------------------------------
+printf "| % 8s | % 40s | % 10s | % 14s |\n" Name "Long Name" Units "From org. vars"
+echo -------------------------------------------------------------------------------------
+for v in $HELPEXTRA; do
+	makevar $v
+	printf "| % 8s | % 40s | % 10s | % 14s |\n" $v "$long_name" "$units" "$invar" 
+done
+echo -------------------------------------------------------------------------------------
+
+echo Variables not in this list are extracted directly using afterburner with no change. 
+echo Examples: zeta,d,bld etc.
+echo Note: uas and vas are not available in PlaSim output, derived from ua and va at sigma=1
+echo "Note: psl (sea level pressure) is not avaiable in PlaSim output," 
+echo "      derived from huss (hus at sigma=1), tas and ps using the hypsometric formula"
+}
+
+function fixvar {
+   local infile=$1
+
+   if [ "$expr" != "none" ] && [ $after != 1 ]; then 
+	  cdo setattribute,$outvar@standard_name="$standard_name" -setattribute,$outvar@long_name="$long_name" -setattribute,$outvar@units="$units" -expr,"$expr" $infile $2
+    else
+	  cp $infile $2
+   fi
+}
+
+function getvar {
+  local infile=$1
+  local outfile=$2
+  local code=$3
+  local lev1=1
+  rm -f burn.log
+  if [ "$no3d" == "true" ]; then
+	 lev1=10
+  fi 
+  if [[ " $VARS3D " =~ .*\ $code\ .* ]] && [ "$no3d" != "true" ]; then
+	if [ $sigma == 0 ]; then
+	   cat > burn$$.nl <<EOF
+	   code=$code,vtype=p,htype=g,mean=0,netcdf=1
+	   hpa=$PLEVS
+EOF
+        else
+	   cat > burn$$.nl <<EOF
+	   code=$code,vtype=s,htype=g,mean=0,netcdf=1
+EOF
+        fi
+  else
+	cat > burn$$.nl <<EOF2
+	code=$code,vtype=s,htype=g,mean=0,netcdf=1
+	modlev=$lev1
+EOF2
+  fi
+  $BURNER $infile $outfile < ./burn$$.nl >> burn.log
+  rm ./burn$$.nl
+}
+
+##################
+## Main 
+##################
+
+sigma=0
+after=0
+while getopts "h?sap:" opt; do
+   case "$opt" in
+        h|\?) help ;;
+        s)  sigma=1 ;;
+        a)  after=1 ;;
+        p)  PLEVS=$OPTARG;;
+   esac
+done
+shift $((OPTIND-1))
+
+if [ $# -lt 3 ] ; then
+help
+exit
+fi
+
+infile=$1 
+outfile=$2
+outvar=$3
+
+makevar $outvar
+
+rm -f temp$$.*.nc
+for var in $invar
+do
+	getvar $infile temp$$.$var.nc $var
+done
+cdo merge temp$$.*.nc temp$$.all.nc 2>> burn.log 
+
+fixvar temp$$.all.nc $outfile $outvar 2>> burn.log
+rm -f temp$$.*.nc

scripts/postmon.sh

@@ -0,0 +1,50 @@
+#!/bin/bash
+# Extracts main variables from PlaSim output and computes monthly means
+# Edit $POSTDIR and $INDIR below
+#set -ex
+
+VARS2D="ps psl tas tasmin tasmax ts uas vas evspsbl pr prc prl prsn mrros mrso tsl rsdt rsut rlut rsds rsus rlds rlus hfls hfss clt prw huss snd snm tauu tauv sit sic"
+VARS3D="hus hur zg ta ua va wa wap cl clw"
+VARSEXTRA="rst rls prl sndc stf mld"
+VARSFX="orog lsm"
+PLEVS="20,30,50,70,100,150,200,250,300,400,500,600,700,850,925,1000"
+
+BURNSH="/home/jost/work/plasim/scripts/burn.sh"
+#BURNSH="/home/jost/work/plasim/scripts/burn.sh -p $PLEVS" # uncomment if you want to specify different levels
+
+if [ $# -lt 3 ] ; then
+        echo Usage: post.sh exp year1 year2	
+	exit
+fi
+
+exp=$1
+y1=$2
+y2=$3
+fn=MOST
+
+POSTDIR="$SCRATCH/plasim/exps/$exp/post"
+INDIR="$WORK/plasim/plasim/exps/$exp"
+
+mkdir -p $POSTDIR
+y1f=$(printf "%03d" $y1)
+
+mkdir -p $POSTDIR/fx
+for var in $VARSFX
+do
+	echo "Processing $var (fixed)"
+     $BURNSH $INDIR/$fn.$y1f $POSTDIR/fx/temp$$.nc $var
+     cdo -s seltimestep,1 $POSTDIR/fx/temp$$.nc $POSTDIR/fx/${exp}_${var}.nc 
+     rm -f $POSTDIR/fx/temp$$.nc 
+done
+
+for var in $VARS2D $VARS3D $VARSEXTRA
+do
+  for (( y=$y1 ; y<=$y2; y++ )); do
+     yf=$(printf "%03d" $y)
+     echo "Processing $var year $yf"
+     mkdir -p $POSTDIR/$yf
+     $BURNSH $INDIR/$fn.$yf $POSTDIR/$yf/temp$$.nc $var
+     cdo -s monmean $POSTDIR/$yf/temp$$.nc $POSTDIR/$yf/${exp}_${yf}_${var}.nc 
+     rm -f $POSTDIR/$yf/temp$$.nc 
+  done
+done

scripts/srv2nc

@@ -0,0 +1,354 @@
+#!/bin/bash
+# PlaSim netcdf postprocessor script
+# (c) 2019 Jost von Hardenberg - ISAC-CNR
+
+#set -ex
+export LC_NUMERIC="en_US.UTF-8"
+
+function help {
+   echo "Convert PLASIM service output file to netcdf format"
+   echo "PlaSIM, ML ocean, ice model and LSG outputs are recognized"
+   echo "Usage: srv2nc [-m] [-y] [-d code1,code2,...] [-f format] infile.srv outfile.nc"
+   echo "Options:"
+   echo "  -m            perform monthly means"
+   echo "  -y            perform yearly means"
+   echo "  -p            compute sea-level pressure"
+   echo "  -d codelist   remove a comma-separated list of codes from the output"
+   echo "  -f format     output format (default is zipped netcdf4 (HDF5))"
+   echo "                Alternatives: nc, nc1, nc4, nc4c, nc5"
+   echo "  -h            print this help"
+   exit
+}
+
+fmonth=""
+fyear=""
+format="nc4 -z zip"
+while getopts "d:ypmhf:" OPTION; do
+    case $OPTION in
+    d) delcodes=$OPTARG ;;
+    m) fmonth="-monmean" ;;
+    y) fyear="-yearmean" ;;
+    p) psl=1 ;;
+    f) format=$OPTARG ;;
+    h) help ;;
+    esac
+done
+shift $((OPTIND-1))
+if [[ ! "nc nc1 nc2 nc4 nc4c nc5 nc4c -z zip nc4 -z zip" == *$format* ]]
+then
+echo "Invalid format: " $format
+echo "Choose one of: nc, nc1, nc2, nc4, nc4c, nc5, nc4c -z zip, nc4 -z zip" 
+exit 1
+fi
+
+if [ $# -lt 2 ]; then
+    help
+fi
+
+infile=$1
+outfile=$2
+
+cdo="cdo -s"
+cdonctab="cdo -t ./param$$.tab -f nc -s"
+cdonc="cdo -f nc -s"
+cdozip="cdo -f $format -s"
+cdoziptab="cdo  -t ./param$$.tab -f $format -s"
+
+cat > param$$.tab << EOT
+110 mld		Mixed Layer Depth		[m]
+129 sg		Surf. Geopotential Orography	[m2/s2]
+130 ta		Temperature			[K]
+131 ua		Zonal Wind			[m/s]
+132 va		Meridional Wind			[m/s]
+133 hus		Specific Humidity		[kg/kg]
+134 ps		Surface Pressure		[hPa]
+135 wap		Vertical Pressure Velocity	[Pa/s]
+137 wa		Vertical Wind			[m/s]
+138 zeta	Vorticity			[1/s]
+139 ts		Surface Temperature		[K]
+140 mrso	Soil Wetness			[m]
+141 snd		Snow Depth			[m]
+142 prl		Large Scale Precipitation	[m/s]
+143 prc		Convective Precipitation	[m/s]
+144 prsn	Snow Fall			[m/s]
+145 bld		Boundary Layer Dissipation	[W/m**2]
+146 hfss	Surface Sensible Heat Flux	[W/m**2]
+147 hfls	Surface Latent Heat Flux	[W/m**2]
+148 stf		Streamfunction			[m**2/s]
+149 psi		Velocity Potential		[m**2/s]
+151 psl		Mean Sea Level Pressure		[hPa]
+152 pl		Log Surface Pressure		[1]
+155 d		Divergence			[1/s]
+156 zg		Geopotential Height		[m]
+157 hur		Relative Humidity		[%]
+158 tps		Tendency of Surface Pressure	[Pa/s]
+159 u3		ustar **3			[m**3/s**3]
+160 mrro	Surface Runoff			[m/s]
+161 clw		Liquid Water Content		[kg/kg]
+162 cl		Cloud Cover			[0-1]
+163 tcc		Total Cloud Cover		[0-1]
+164 clt		Total Cloud Cover (Mean)	[0-1]
+165 uas		Eastward Wind 10m		[m/s]
+166 vas		Northward Wind 10m		[m/s]
+167 tas		Temperature at 2m		[K]
+168 td2m	Dew Point Temperature at 2m	[K]
+169 tsa		Surface Temperature Accumulated	[K]
+170 tsod	Deep Soil Temperature		[K]
+171 dsw		Deep Soil Wetness		[1]
+172 lsm		Land Sea Mask			[0-1]
+173 z0		Surface Roughness		[m]
+174 alb		Surface Albedo			[1]
+175 as		Surface Albedo			[1]
+176 rss		Surface Solar Radiation		[W/m2]
+177 rls		Surface Thermal Radiation	[W/m2]
+178 rst		Top Solar Radiation		[W/m2]
+179 rlut	Top Thermal Radiation		[W/m2]
+180 tauu	U-Stress			[Pa]
+181 tauv	V-Stress			[Pa]
+182 evap	Evaporation			[m/s]
+183 tso		Soil Temperature		[K]
+184 wsoi	Soil Wetness			[1]
+199 vegc	Vegetation Cover		[0-1]
+201 tasmax	Maximum daily 2m temperature	[K]
+202 tasmin	Minimum daily 2m temperature	[K]
+203 rsut	Top Solar Radiation Upward	[W/m2]
+204 ssru	Surface Solar Radiation Upward	[W/m2]
+205 stru	Surface Therm Radiation Upward	[W/m2]
+207 tso2	Soil Temperature Level 2	[K]
+208 tso3	Soil Temperature Level 3	[K]
+209 tso4	Soil Temperature Level 4	[K]
+210 sic		Sea Ice Cover			[0-1]
+211 sit		Sea Ice Thickness		[m]
+212 vegf	Forest Cover			[0-1]
+218 snm		Snow Melt			[m/s]
+221 sndc	Snow Depth Change		[m/s]
+229 dwmax	Field capacity			[1] 
+230 prw		Vert. Integrated Spec. Hum.	[kg/m2]
+232 glac	Glacier Cover			[0-1]
+238 tsn		Snow temperature		[K]
+259 spd		Wind Speed			[m/s]
+260 pr		Total Precipitation		[m/s]
+261 ntr		Net Top Radiation		[W/m2]
+262 nbr		Net Bottom Radiation		[W/m2]
+263 hfns	Net Heat Flux			[W/m2]
+264 wfn		Net Water Flux			[m/s]
+273 dpdx	d(ps)/dx			[Pa/m]
+273 dpdy	d(ps)/dy			[Pa/m]
+277 hlpr	Half level pressure		[Pa]
+278 flpr	Full level pressure		[Pa]
+701 heata	Flux from the atmosphere	[W/m2]
+702 ofluxa	Flux from the ocean		[W/m2]
+703 tsfluxa	Flux to warm/cool ice/snow	[W/m2]
+704 smelta	Flux for snow melt		[W/m2]
+705 imelta	Flux used for ice melt		[W/m2]
+706 cfluxa	Flux to the ocean		[W/m2]
+707 fluxca	Cond. heatflux			[W/m2]
+708 qmelta	Res flux to ice			[W/m2]
+709 xflxicea	Flux correction			[W/m2]
+710 icec	Ice cover			[0-1]
+711 iced	Ice thickness			[m]
+712 scflxa	Flux from snow -> ice conversion	[W/m2]
+713 cfluxra	Flux for limiting ice to xmaxd	[W/m2]
+714 cfluxna	Flux due to neg. ice		[W/m2]
+739 ts		Surface temperature		[K]
+741 zsnow	Snow depth			[m h2o]
+769 sst		Sea surface temperature		[K]
+772 ls		Land sea mask			[0-1]
+790 clicec2	Climatological ice cover	[0-1]
+791 cliced2	Climatological ice thickness	[m]
+792 icecc	Ice cover computed prognostically	[frac.]
+794 cpmea	Fresh water (p-e) for LSG	[m/s]
+795 croffa	Runoff for LSG			[m/s]
+796 stoia	Snow converted to ice		[m h2o]
+901 heata	Heat flux from atm/ice		[W/m2]
+902 ifluxa	Heat flux into ice		[W/m2]
+903 fssta	Flux correction			[W/m2]
+904 dssta	Vertical diffusion		[m2/s]
+905 qhda	Horizontal diffusion		[m2/s]
+906 fldoa	Flux from deep ocean (LSG)	[W/m2]
+910 icec	Sea ice cover			[0-1]
+939 sst		SST (ML ocean temperature)	[K]
+972 ls		Land-sea mask			[1]
+990 clsst	Climatological SST		[K]
+EOT
+cat > paramlsg$$.tab << EOT
+1 zeta		surface elevation		[m]
+2 t		potential temperature		[K]
+3 utot		zonal velocity component	[m/s]
+4 vtot		meridional velocity component	[m/s] 
+5 s		salinity			[0/00]
+7 w		vertical velocity component	[m/s] 
+13 sice		ice thickness			[m]
+18 fluxhea	heat flux			[W/m2] 
+27 psi		horizontal barotropic stream function	[m3/s] 
+37 ub		zonal component of barotropic velocity	[m/s] 
+38 vb		meridional component of barotropic velocity	[m/s]
+40 wet		land-sea mask for scalar points	[0-1]
+41 wetvec	land-sea mask for vector points	[0-1] 
+52 taux		zonal component of wind stress	[Pa] 
+53 tauy		meridional component of wind stress	[Pa] 
+62 tc		potential temperature		[C]
+65 fluwat	net fresh water flux (P-E)	[m/s] 
+66 convadd	potential energy dissipation due to convection	[mW/m2]
+67 flukwat	fresh water flux due to Newtonian coupling	[m/s] 
+68 flukhea	heat flux due to Newtonian coupling	[W/m2] 
+69 convad	convective adjustment events	[1]
+80 fldsst	sst differences lsg-plasim	[K]
+81 fldice	ice-differences lsg-plasim	[1]
+82 fldpme	pme-differences lsg-plasim	[m/s]
+83 fldtaux	taux-differences lsg-plasim	[Pa]
+84 fldtauy	tauy-differences lsg-plasim	[Pa]
+92 tbound	boundary value of t		[K]
+97 dQ/dt 	coupling coefficient		[W/m2K] 
+98 depp		depth in scalar-points		[m] 
+99 depv		depth in vector-points		[m] 
+EOT
+
+ysize=$($cdo -s griddes $infile | grep ysize | head -1 | cut -d= -f2|sed 's/ //g')
+case $ysize in
+    32) spgrid=t21; gpgrid=n16;;
+    48) spgrid=t31; gpgrid=n24;;
+    64) spgrid=t42; gpgrid=n32;;
+    96) spgrid=t63; gpgrid=n48;;
+    128) spgrid=t85; gpgrid=n64;;
+    160) spgrid=t106; gpgrid=n80;;
+    76)  gpgridt=lsggridt$$.txt; gpgridu=lsggridu$$.txt; flsg=1;;
+esac
+
+if [ -z "$flsg" ]; then 
+    fsp=$($cdo griddes $infile | grep "gridID 2") # Spectral or not?
+    $cdo zaxisdes  $infile > inzgrid$$.txt
+    f3d=$(grep "zaxisID 2" inzgrid$$.txt) # Vertical sigma axis or not
+else
+
+# vertical axes for LSG
+cat > zaxislsgu$$.txt << EOT
+zaxistype = depth_below_sea
+size = 22
+levels = 25 75 125 175 225 275 350 450 550 650 750 850 950 1100 1300 1500 1800 2250 2750 3500 4500 5500 
+EOT
+cat > zaxislsgc$$.txt << EOT
+zaxistype = depth_below_sea
+size = 21
+levels = 75 125 175 225 275 350 450 550 650 750 850 950 1100 1300 1500 1800 2250 2750 3500 4500 5500 
+EOT
+cat > zaxislsgw$$.txt << EOT
+zaxistype = depth_below_sea
+size = 22
+levels = 50 100 150 200 250 312.5 400 500 600 700 800 900 1025 1200 1400 1650 2025 2500 3125 4000 5000 6000
+EOT
+
+    # Create scalar grid
+    cat > $gpgridt <<EOT
+#
+# gridID 1
+#
+gridtype  = curvilinear
+gridsize  = 5472
+xsize     = 72
+ysize     = 76
+xname     = lon
+xdimname  = west_east
+xlongname = "longitude"
+xunits    = "degrees_east"
+yname     = lat
+ydimname  = south_north
+ylongname = "latitude"
+yunits    = "degrees_north"
+xvals     =
+ $(awk "BEGIN {off=-2.5; for(y=93.75; y>=-93.75; y=y-2.5) {for(x=0.; x<=357.5; x=x+5.0) {printf x-off \" \" }; printf \"\n\"; off=2.5+off}}")
+yvals     = 
+ $(awk "BEGIN {for(y=93.75; y>=-93.75; y=y-2.5) {for(i=1; i<=72; ++i) {printf y \" \" } printf \"\n\"}}")
+scanningMode = 64
+EOT
+
+    # Create vector grid
+    cat > $gpgridu <<EOT
+#
+# gridID 1
+#
+gridtype  = curvilinear
+gridsize  = 5472
+xsize     = 72
+ysize     = 76
+xname     = lon
+xdimname  = west_east
+xlongname = "longitude"
+xunits    = "degrees_east"
+yname     = lat
+ydimname  = south_north
+ylongname = "latitude"
+yunits    = "degrees_north"
+xvals     =
+ $(awk "BEGIN {off=-5.; for(y=93.75; y>=-93.75; y=y-2.5) {for(x=0.; x<=357.5; x=x+5.0) {printf x-off \" \" }; printf \"\n\"; off=2.5+off}}")
+yvals     = 
+ $(awk "BEGIN {for(y=93.75; y>=-93.75; y=y-2.5) {for(i=1; i<=72; ++i) {printf y \" \" } printf \"\n\"}}")
+scanningMode = 64
+EOT
+
+fi
+
+# If 3D Field, adjust vertical grid
+if [ ! -z "$f3d" ]; then
+nlev=$( grep size inzgrid$$.txt|tail -1|cut -d= -f2|sed 's/ //g' )
+#sigmas=$($cdo outputf,%19.16f -selindexbox,1,$nlev,1,1  -selcode,333 $infile | awk ' BEGIN{p=0.} {s=($1+p)/2.; p=$1} {printf("%19.16f\n", s)}')
+sigmas=$($cdo outputf,%19.16f -selindexbox,1,$nlev,1,1  -selcode,333 $infile)
+grep -m 1 -B 1000 -A 1 "zaxisID 2" inzgrid$$.txt >zgrid$$.txt
+cat >> zgrid$$.txt << EOT
+zaxistype = hybrid
+size      = $nlev
+levels    = $( seq 1 $nlev )
+vctsize   = $(( nlev*2 + 2 ))
+vct       = $( printf "%.s0 " $(seq 1 $((nlev + 1)) ) )
+            0 $sigmas
+EOT
+fixzgrid="-setzaxis,zgrid$$.txt"
+delcode333="-delcode,333"
+else
+fixzgrid=""
+delcode333=""
+fi
+
+if [ -z "$flsg" ]; then
+# If not LSG
+# If spectral fix grids
+if [ ! -z "$fsp" ]; then
+  $cdo splitgrid ${infile} grid_$$_
+  if [ ! -z "$psl" ]; then
+     $cdonctab $fixzgrid -setgrid,$spgrid grid_$$_02.srv spectral$$.nc
+     $cdonctab $fixzgrid $delcode333 -setgrid,$gpgrid grid_$$_01.srv gaussian$$.nc
+     $cdonctab merge spectral$$.nc gaussian$$.nc tempfile0$$.nc
+     $cdonctab geopotheight -sp2gp tempfile0$$.nc geopot$$.nc
+     $cdonctab sealevelpressure -sp2gp tempfile0$$.nc press$$.nc
+     $cdozip $fmonth $fyear -merge tempfile0$$.nc geopot$$.nc press$$.nc tempfile$$.nc
+     #$cdozip $fmonth $fyear -merge tempfile0$$.nc press$$.nc tempfile$$.nc 
+  else
+     $cdonctab $fmonth $fyear $fixzgrid -setgrid,$spgrid grid_$$_02.srv spectral$$.nc
+     $cdonctab $fmonth $fyear $fixzgrid $delcode333 -setgrid,$gpgrid grid_$$_01.srv gaussian$$.nc
+     $cdozip merge spectral$$.nc gaussian$$.nc tempfile$$.nc
+  fi
+else
+   $cdoziptab setgrid,$gpgrid $fmonth $fyear $infile tempfile$$.nc
+fi
+else
+# If LSG
+   #$cdonc copy $fmonth $fyear $infile tempinfile$$.nc
+   $cdonc splitzaxis $infile zaxis$$_
+   $cdonc setzaxis,zaxislsgu$$.txt -delcode,7 zaxis$$_01.nc fixzaxis$$_01u.nc # code 7 is w
+   $cdonc setzaxis,zaxislsgw$$.txt -selcode,7 zaxis$$_01.nc fixzaxis$$_01w.nc
+   $cdonc setzaxis,zaxislsgc$$.txt zaxis$$_03.nc fixzaxis$$_03.nc
+   $cdonc merge fixzaxis$$_01u.nc fixzaxis$$_01w.nc zaxis$$_02.nc fixzaxis$$_03.nc tempin$$.nc
+   $cdonc setpartab,paramlsg$$.tab $fmonth $fyear tempin$$.nc tempinpar$$.nc
+   uvars="taux,tauy,depv,ub,vb,utot,vtot,wetvec"
+   $cdonc setgrid,$gpgridu $fmonth $fyear -selname,$uvars tempinpar$$.nc ufields$$.nc
+   $cdonc setgrid,$gpgridt $fmonth $fyear -delname,$uvars tempinpar$$.nc tfields$$.nc
+   $cdozip merge tfields$$.nc ufields$$.nc tempfile$$.nc
+fi
+
+if [ -z "$delcodes" ]
+then
+   mv tempfile$$.nc $outfile
+else
+   $cdozip delcode,$delcodes tempfile$$.nc $outfile
+fi
+rm -f spectral$$.nc gaussian$$.nc grid_$$_01.srv grid_$$_02.srv inzgrid$$.txt zgrid$$.txt tempfile$$.nc param$$.tab lsggridt$$.txt  lsggridu$$.txt ufields$$.nc tfields$$.nc zaxislsgu$$.txt zaxislsgw$$.txt zaxislsgc$$.txt zaxis$$_* fixzaxis$$_* tempin$$.nc tempinpar$$.nc paramlsg$$.tab press$$.nc geopot$$.nc tempfile0$$.nc