Commit 3f3c6f2e authored by ph290's avatar ph290
Browse files

update

parent 3a16a08e
import os
import subprocess
import shutil
import glob
from math import cos, asin, sqrt
import multiprocessing as mp
from functools import partial
import uuid
##################################################
# you may need to change things here #
##################################################
base_directory = '/gpfs/ts0/home/ph290/s2p3_rv2.0/'
num_procs = mp.cpu_count() # this will use all available processors. Note that on a multi-node machine it can only use the processors on one node
# num_procs = 1
# output_directory = '/gpfs/ts0/projects/Research_Project-148395/s2p3_rv2.0/output/era5/test/' #where you want the output to go
output_directory = '/gpfs/ts0/projects/Research_Project-148395/s2p3_rv2.0/output/radiation/test_drag/' #where you want the output to go
#output_file_name = 'global_tropics_era5_test'
#meterological_file_name = 'meterological_data'
#domain_file_name = 's12_m2_s2_n2_h_map_test.dat'
#nutrient_file_name = 'initial_nitrate_test.dat'
#executable_file_name = 's2p3_rv2.0'
output_file_name = 'drag_test'
meterological_file_name = 'meterological_data'
domain_file_name = 's12_m2_s2_n2_h_map_yucatan.dat'
nutrient_file_name = 'initial_nitrate_yucatan.dat'
executable_file_name = 's2p3_rv2.0'
# met_data_location = '/gpfs/ts0/projects/Research_Project-148395/s2p3_rv2.0/era5/test/' # The location containing the tar.gz met files (in the format met_data_year.tar.gz)
met_data_location = '/gpfs/ts0/projects/Research_Project-148395/s2p3_rv2.0/yucatan/' # The location containing the tar.gz met files (in the format met_data_year.tar.gz)
met_data_temporary_location = base_directory+'met/spatial_data/' # The location that met data for each year will be un tar.gziped into
start_year = 1997
#having to restaer becase of isca restart
end_year = 2001 # same as start year resuls in a 1 year run
depth_min = 10 # NOTE that these numbers MUST be the same as those used in the scripts used to produce the meterology and nutrient files, otherwse data will not be taken for teh correct lats/lons and/or the script will fail
depth_max = 50
write_error_output = False
parallel_processing = True
#######################################################
# Variables to output from model #
# =1 means output, =0 means do not output #
# the columns in the output will be left to right in #
# the order that the variables are here top-to-bottom #
#######################################################
include_depth_output=1
include_temp_surface_output=1
include_temp_bottom_output=1
include_chlorophyll_surface_output=1
include_phyto_biomass_surface_output=0
include_phyto_biomass_bottom_output=0
include_PAR_surface_output=1 # daily mean surface of PAR W m-2
include_PAR_bottom_output=1 # daily mean bottom of PAR W m-2
include_windspeed_output=0 # windspeed
include_stressx_output=0 # x component of surface wind drag
include_stressy_output=0 # y component of surface wind drag
include_Etide_output=0 # Mixing power in the tidal currents (assumes constant mixing efficiency 0.003)
include_Ewind_output=0 # Mixing power in the wind (assumes constant mixing efficiency 0.023 and slippage factor=0.025)
include_u_mean_surface_output=0 # daily mean surface u in cm/s !!
include_u_mean_bottom_output=0 # daily mean bottom u in cm/s !!
include_grow1_mean_surface_output=0 # daily mean growth rate d-1 surface
include_grow1_mean_bottom_output=0 # daily mean growth rate d-1 bottom
include_uptake1_mean_surface_output=0 # daily mean DIN uptake rate mmol DIN (mg C)-1 d-1 surface
include_uptake1_mean_bottom_output=0 # daily mean DIN uptake rate mmol DIN (mg C)-1 d-1 bottom
include_tpn1_output=0 # total water column net production / mg C m-2 d-1
include_tpg1_output=0 # total water column gross production / mg C m-2 hd-1
include_speed3_output=0 # depth-mean current speed
##################################################
# functions used by the script #
##################################################
def distance(lat1, lon1, lat2, lon2):
p = 0.017453292519943295
a = 0.5 - cos((lat2-lat1)*p)/2 + cos(lat1*p)*cos(lat2*p) * (1-cos((lon2-lon1)*p)) / 2
return 12742 * asin(sqrt(a))
def closest(data, lat1,lon1):
return min(data, key=lambda p: distance(lat1,lon1,p[0],p[1]))
# def return_domain_lon(filename,i):
# f=open(filename)
# lines=f.readlines()
# return filter(lambda a: a != '', lines[i+1].split(' '))[0:2]
#
#
f=open(base_directory+'domain/'+domain_file_name)
lines=f.readlines()
f2=open(base_directory+'domain/'+nutrient_file_name)
lines2=f2.readlines()
lat_domain=[]
lon_domain=[]
alldepth=[]
smaj1=[]
smin1=[]
smaj2=[]
smin2=[]
smaj3=[]
smin3=[]
smaj4=[]
smin4=[]
smaj5=[]
smin5=[]
woa_nutrient=[]
counter = 0
for i,line in enumerate(lines[1::]):
depth = float(line[77:84])
if ((depth >= depth_min) & (depth <= depth_max) & (depth > 0.0)):
lon_domain.append(line[0:8])
lat_domain.append(line[8:16])
alldepth.append(line[77:84])
smaj1.append(line[16:22])
smin1.append(line[22:28])
smaj2.append(line[28:34])
smin2.append(line[34:40])
smaj3.append(line[40:46])
smin3.append(line[46:52])
smaj4.append(line[52:58])
smin4.append(line[58:64])
smaj5.append(line[64:70])
smin5.append(line[70:76])
woa_nutrient.append(lines2[counter][16:22])
counter += counter
def run_model(domain_file_name,lats_lons,year,start_year,unique_job_id,met_data_temporary_location,lon_domain,lat_domain,smaj1,smin1,smaj2,smin2,smaj3,smin3,smaj4,smin4,smaj5,smin5,woa_nutrient,alldepth,include_depth_output,include_temp_surface_output,include_temp_bottom_output,include_chlorophyll_surface_output,include_phyto_biomass_surface_output,include_phyto_biomass_bottom_output,include_PAR_surface_output,include_PAR_bottom_output,include_windspeed_output,include_stressx_output,include_stressy_output,include_Etide_output,include_Ewind_output,include_u_mean_surface_output,include_u_mean_bottom_output,include_grow1_mean_surface_output,include_grow1_mean_bottom_output,include_uptake1_mean_surface_output,include_uptake1_mean_bottom_output,include_tpn1_output,include_tpg1_output,include_speed3_output,i):
#modifying so that the fortran code looks for the correct met file, rather than us having to copy it into the working dorectory
# lon,lat = return_domain_lon(base_directory+'domain/'+domain_file_name,i)
lon_domain_tmp = float(lon_domain[i])
if lon_domain_tmp < 0.0:
lon_domain_tmp = 360.0+lon_domain_tmp
# forcing_lat_lon = closest(lats_lons, float(lat_domain[i]),float(lon_domain[i]))
# print i
# print str(forcing_lat_lon[0])
# print str(forcing_lat_lon[1])
# print str(float(lat_domain[i]))
# print str(lon_domain_tmp)
run_command = """./"""+executable_file_name+""" << EOF
"""+str(start_year)+"""
"""+str(year)+"""
"""+str(float(lat_domain[i]))+"""
"""+str(lon_domain_tmp)+"""
../domain/"""+domain_file_name+"""
../domain/"""+nutrient_file_name+"""
"""+unique_job_id+"""
"""+met_data_temporary_location+"""
map
"""+str(i+1)+"""
"""+str(smaj1[i])+"""
"""+str(smin1[i])+"""
"""+str(smaj2[i])+"""
"""+str(smin2[i])+"""
"""+str(smaj3[i])+"""
"""+str(smin3[i])+"""
"""+str(smaj4[i])+"""
"""+str(smin4[i])+"""
"""+str(smaj5[i])+"""
"""+str(smin5[i])+"""
"""+str(woa_nutrient[i])+"""
"""+str(alldepth[i])+"""
"""+str(include_depth_output)+"""
"""+str(include_temp_surface_output)+"""
"""+str(include_temp_bottom_output)+"""
"""+str(include_chlorophyll_surface_output)+"""
"""+str(include_phyto_biomass_surface_output)+"""
"""+str(include_phyto_biomass_bottom_output)+"""
"""+str(include_PAR_surface_output)+"""
"""+str(include_PAR_bottom_output)+"""
"""+str(include_windspeed_output)+"""
"""+str(include_stressx_output)+"""
"""+str(include_stressy_output)+"""
"""+str(include_Etide_output)+"""
"""+str(include_Ewind_output)+"""
"""+str(include_u_mean_surface_output)+"""
"""+str(include_u_mean_bottom_output)+"""
"""+str(include_grow1_mean_surface_output)+"""
"""+str(include_grow1_mean_bottom_output)+"""
"""+str(include_uptake1_mean_surface_output)+"""
"""+str(include_uptake1_mean_bottom_output)+"""
"""+str(include_tpn1_output)+"""
"""+str(include_tpg1_output)+"""
"""+str(include_speed3_output)+"""
"""+str(start_year)+"""
"""+str(year)+"""
"""+str(float(lat_domain[i]))+"""
"""+str(lon_domain_tmp)+"""
../domain/"""+domain_file_name+"""
../domain/"""+nutrient_file_name+"""
"""+unique_job_id+"""
"""+met_data_temporary_location+"""
map
"""+str(i+1)+"""
"""+str(smaj1[i])+"""
"""+str(smin1[i])+"""
"""+str(smaj2[i])+"""
"""+str(smin2[i])+"""
"""+str(smaj3[i])+"""
"""+str(smin3[i])+"""
"""+str(smaj4[i])+"""
"""+str(smin4[i])+"""
"""+str(smaj5[i])+"""
"""+str(smin5[i])+"""
"""+str(woa_nutrient[i])+"""
"""+str(alldepth[i])+"""
"""+str(include_depth_output)+"""
"""+str(include_temp_surface_output)+"""
"""+str(include_temp_bottom_output)+"""
"""+str(include_chlorophyll_surface_output)+"""
"""+str(include_phyto_biomass_surface_output)+"""
"""+str(include_phyto_biomass_bottom_output)+"""
"""+str(include_PAR_surface_output)+"""
"""+str(include_PAR_bottom_output)+"""
"""+str(include_windspeed_output)+"""
"""+str(include_stressx_output)+"""
"""+str(include_stressy_output)+"""
"""+str(include_Etide_output)+"""
"""+str(include_Ewind_output)+"""
"""+str(include_u_mean_surface_output)+"""
"""+str(include_u_mean_bottom_output)+"""
"""+str(include_grow1_mean_surface_output)+"""
"""+str(include_grow1_mean_bottom_output)+"""
"""+str(include_uptake1_mean_surface_output)+"""
"""+str(include_uptake1_mean_bottom_output)+"""
"""+str(include_tpn1_output)+"""
"""+str(include_tpg1_output)+"""
"""+str(include_speed3_output)+"""
EOF"""
# print run_command
proc = subprocess.Popen([run_command], shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
out, err = proc.communicate()
# return out
return out,err
##################################################
# main program #
##################################################
unique_job_id = str(uuid.uuid4())
num_lines = sum(1 for line in open(base_directory+'domain/'+domain_file_name)) - 1
# num_lines = 10
print 'initial unzipping of met data to extrcat lat and lon data'
subprocess.call('tar -C '+met_data_temporary_location+' -zxf '+met_data_location+'met_data_'+str(start_year)+'.tar.gz', shell=True)
files = glob.glob(met_data_temporary_location+'*_'+str(start_year)+'.dat')
w, h = 2, len(files) ;
lats_lons = [[0 for x in range(w)] for y in range(h)]
for i,file in enumerate(files):
tmp = file.split('lat')[-1].split('.dat')[0].split('lon')
lats_lons[i][0] = float(tmp[0])
lats_lons[i][1] = float(tmp[1].split('_')[0])
# ##testing
# year=2006
# pool = mp.Pool(processes=num_procs)
# func = partial(run_model, domain_file_name, lats_lons, year, start_year, unique_job_id, met_data_temporary_location,lon_domain,lat_domain,smaj1,smin1,smaj2,smin2,smaj3,smin3,smaj4,smin4,smaj5,smin5,woa_nutrient,alldepth,include_depth_output,include_temp_surface_output,include_temp_bottom_output,include_chlorophyll_surface_output,include_phyto_biomass_surface_output,include_phyto_biomass_bottom_output,include_PAR_surface_output,include_PAR_bottom_output,include_windspeed_output,include_stressx_output,include_stressy_output,include_Etide_output,include_Ewind_output,include_u_mean_surface_output,include_u_mean_bottom_output,include_grow1_mean_surface_output,include_grow1_mean_bottom_output,include_uptake1_mean_surface_output,include_uptake1_mean_bottom_output,include_tpn1_output,include_tpg1_output,include_speed3_output)
# results, errors = zip(*pool.map(func, range(1)))
# print results[0].split('\n')[0]
# print results[0].split('\n')[1]
# print results[0].split('\n')[2]
# # print errors[0].split('\n')[0]
# pool.close()
print 'looping through years'
for year in range(start_year,end_year+1):
print year
#clean up and prexisting met files
try:
files_to_delete = glob.glob(met_data_temporary_location+'*.dat')
[os.remove(f) for f in files_to_delete]
except:
print 'no met files to clean up'
subprocess.call('tar -C '+met_data_temporary_location+' -zxf '+met_data_location+'met_data_'+str(year)+'.tar.gz', shell=True)
try:
shutil.move(output_directory+output_file_name+'_'+str(year), output_directory+output_file_name+'_'+str(year)+'_previous')
except:
print 'no previous output file to move'
if parallel_processing:
pool = mp.Pool(processes=num_procs)
func = partial(run_model, domain_file_name, lats_lons, year, start_year, unique_job_id, met_data_temporary_location,lon_domain,lat_domain,smaj1,smin1,smaj2,smin2,smaj3,smin3,smaj4,smin4,smaj5,smin5,woa_nutrient,alldepth,include_depth_output,include_temp_surface_output,include_temp_bottom_output,include_chlorophyll_surface_output,include_phyto_biomass_surface_output,include_phyto_biomass_bottom_output,include_PAR_surface_output,include_PAR_bottom_output,include_windspeed_output,include_stressx_output,include_stressy_output,include_Etide_output,include_Ewind_output,include_u_mean_surface_output,include_u_mean_bottom_output,include_grow1_mean_surface_output,include_grow1_mean_bottom_output,include_uptake1_mean_surface_output,include_uptake1_mean_bottom_output,include_tpn1_output,include_tpg1_output,include_speed3_output)
# results,errors = pool.map(func, range(num_lines))
results, errors = zip(*pool.map(func, range(len(lat_domain))))
# results = pool.map(func, range(num_lines))
with open(output_directory+output_file_name+'_'+str(year),'w') as fout:
for result in results:
fout.write(result)
if write_error_output:
with open(output_directory+output_file_name+'_error_'+str(year),'w') as fout:
for error in errors:
fout.write(error)
pool.close()
if not parallel_processing:
# non parallel version
with open(output_directory+output_file_name+'_'+str(year),'w') as fout:
for i in range(len(lat_domain)):
out,err = run_model(domain_file_name, lats_lons, year, start_year, unique_job_id, met_data_temporary_location,lon_domain,lat_domain,smaj1,smin1,smaj2,smin2,smaj3,smin3,smaj4,smin4,smaj5,smin5,woa_nutrient,alldepth,include_depth_output,include_temp_surface_output,include_temp_bottom_output,include_chlorophyll_surface_output,include_phyto_biomass_surface_output,include_phyto_biomass_bottom_output,include_PAR_surface_output,include_PAR_bottom_output,include_windspeed_output,include_stressx_output,include_stressy_output,include_Etide_output,include_Ewind_output,include_u_mean_surface_output,include_u_mean_bottom_output,include_grow1_mean_surface_output,include_grow1_mean_bottom_output,include_uptake1_mean_surface_output,include_uptake1_mean_bottom_output,include_tpn1_output,include_tpg1_output,include_speed3_output,i)
fout.write(result)
#clean up and leftover met files
try:
files_to_delete = glob.glob(met_data_temporary_location+'*.dat')
[os.remove(f) for f in files_to_delete]
except:
print 'no met files to clean up'
remove_files = glob.glob(base_directory+'main/*'+unique_job_id+'*')
try:
remove_files.remove(base_directory+'main/restart'+unique_job_id+'.dat')
except:
pass
for remove_file in remove_files:
os.remove(remove_file)
# for year in range(start_year,end_year+1):
# pool = mp.Pool(processes=num_procs)
# func = partial(run_model2, domain_file_name, lats_lons,year,start_year,unique_job_id, met_data_temporary_location)
# # results,errors = pool.map(func, range(num_lines))
# results = pool.map(func, range(num_lines))
# with open(output_directory+output_file_name+'_'+str(year)+'_error','w') as fout:
# for result in results:
# fout.write(result)
#!/bin/sh
#PBS -V # export all environment variables to the batch job.
#PBS -d . # set working directory to .
#PBS -q ptq # submit to the Parallel queue
#PBS -l walltime=00:00:20:00 # Maximum wall time for the job. days: hours: minutes:seconds
#PBS -A Research_Project-148395 # research project to submit under.
#PBS -l nodes=1:ppn=16 # or nodes=number of nodes required. ppn=number of processors per node
#PBS -m e -M ph290@exeter.ac.uk # email me at job completion
python2.7 run_map_parallel_radiation.py
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment