From d9121d0dc85b196d0c4ebb3488b3d1216d9db628 Mon Sep 17 00:00:00 2001 From: vsicardi Date: Wed, 21 Nov 2018 11:33:56 +0100 Subject: [PATCH] add Oce_IC (https://earth.bsc.es/gitlab/ec-earthCPWG/Oce_IC.git) to cpg_tools --- Oce_IC/compute_density.py | 36 +++ Oce_IC/fetch_ORAS4.sh | 68 ++++++ Oce_IC/fetch_ORAS4_EC23.sh | 67 ++++++ Oce_IC/library_interp_extrap.py | 398 ++++++++++++++++++++++++++++++++ Oce_IC/prep_restart_interp.py | 179 ++++++++++++++ 5 files changed, 748 insertions(+) create mode 100644 Oce_IC/compute_density.py create mode 100755 Oce_IC/fetch_ORAS4.sh create mode 100755 Oce_IC/fetch_ORAS4_EC23.sh create mode 100644 Oce_IC/library_interp_extrap.py create mode 100644 Oce_IC/prep_restart_interp.py diff --git a/Oce_IC/compute_density.py b/Oce_IC/compute_density.py new file mode 100644 index 0000000..438bbd0 --- /dev/null +++ b/Oce_IC/compute_density.py @@ -0,0 +1,36 @@ +import numpy as np +from numpy import ma as ma +def compute_density(ptemp, psal): + """Compute the density using the Jackett and McDougall (1994) equation of state + (adapted from fortran routine eosbn2 from NEMO 3.3) + by C. Prodhomme -- 2016 + Arguments: + ptemp: numpy array of potential temperature + psal: numpy array of potential salinity + Output: + zrhop: numpy array of density + """ + + #??? + rau0 = 1035 + + zrau0r = 1.e0 / rau0 + #square root salinity + psalr = np.sqrt( np.absolute(psal)) + + + # compute volumic mass pure water at atm pressure + zr1= ( ( ( ( 6.536332e-9*ptemp-1.120083e-6 )*ptemp+1.001685e-4 )*ptemp + -9.095290e-3 )*ptemp+6.793952e-2 )*ptemp+999.842594 + # seawater volumic mass atm pressure + zr2= ( ( ( 5.3875e-9*ptemp-8.2467e-7 )*ptemp+7.6438e-5 ) *ptemp + -4.0899e-3 ) *ptemp+0.824493 + zr3= ( -1.6546e-6*ptemp+1.0227e-4 )*ptemp-5.72466e-3 + zr4= 4.8314e-4 + + # potential volumic mass (reference to the surface) + zrhop= ( zr4*psal + zr3*psalr + zr2 ) *psal + zr1 + + return(zrhop) + + diff --git a/Oce_IC/fetch_ORAS4.sh b/Oce_IC/fetch_ORAS4.sh new file mode 100755 index 0000000..4fd6fc4 --- /dev/null +++ b/Oce_IC/fetch_ORAS4.sh @@ -0,0 +1,68 @@ +#!/bin/ksh +#SBATCH -n 1 +#SBATCH -t 78:00:00 +#SBATCH -J fetchORAS4 +#SBATCH -o fetch-oras4-%j.log +set -exv +#module load netCDF-Fortran/4.2-foss-2015a +# This script download the s4 restarts from ECMWF +# and launch the interpolation of the files toward +# the target resolution +# History : Chloe Prodhomme - Initial version 2016 +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# Arguments +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +YEARI=2005 # first years of the restarts to be downloaded +YEARF=2016 # Last year to download +MEMBER_LST="3" # List of members to download +MON_LST="1231" #"0228 0331 0430 0531 0630 0731 0831 0930 1031 1130 1231" #"0131 0228 0331 0430 0531 0630 0731 0831 1031 1130 1231" # List of months to download in formate MMDD +config=1L75 #target resolution for the output grid +modelversion="Ec3.2" +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +workdir=/scratch/Earth/$USER/ORAS4/ +rm -f /scratch/Earth/$USER/ORAS4/* +mkdir -p $workdir +cd $workdir + +DATASOURCE="ec:/emos/OCEA/4/0001/" # Root of the path where to find ORAS4 on ECFS +#loop over the members +for mem in $MEMBER_LST ; do + #loop over the month + for mon in $MON_LST ; do + #loop over the years + for YEAR in `seq $YEARI $YEARF`; do + #directory for orginal oras4 files + localdir=/esnas/releases/ic/ocean/ORCA1/s4_original_file/ + #name of the output file + fileout=s4_fc${mem}_${YEAR}${mon}_restart.nc + #check if the file already on the /esnas + if [ ! -e $localdir/${fileout}.gz ]; then + #get the exact file name at ECMWF + file=$(ecaccess-file-dir ${DATASOURCE}/opa${mem}/restart/${YEAR}/0001_${YEAR}${mon}_*_restart.tar.gz) + #download the file + ecaccess-file-get ${DATASOURCE}/opa${mem}/restart/${YEAR}/$file $file + #unzip and untar + gunzip $file + tar -xvf ${file/.gz/} + #remove not needed files + rm -f *dat ${file/.gz/} + #store the file + mv 0001_${YEAR}${mon}_*_restart.nc s4_fc${mem}_${YEAR}${mon}_restart.nc + cp s4_fc${mem}_${YEAR}${mon}_restart.nc ${localdir}/$fileout + gzip ${localdir}/$fileout + else + #copy the file into the workdir + cp $localdir/${fileout}.gz . + #unzip the file + gunzip ${fileout}.gz + + fi + #Launch the interpolation toward the target resolution + python prep_restart_interp.py ${fileout} /esnas/autosubmit/con_files/mesh_mask_nemo.nemovar_O1L42.nc $config $modelversion + # remove the local files + rm -rf ${fileout}.gz + done #loop year + done #loop month +done #loop member + diff --git a/Oce_IC/fetch_ORAS4_EC23.sh b/Oce_IC/fetch_ORAS4_EC23.sh new file mode 100755 index 0000000..e9e2041 --- /dev/null +++ b/Oce_IC/fetch_ORAS4_EC23.sh @@ -0,0 +1,67 @@ +#!/bin/ksh +#SBATCH -n 1 +#SBATCH -t 78:00:00 +#SBATCH -J fetchORAS4 +#SBATCH -o fetch-oras4-%j.log +set -exv +#module load netCDF-Fortran/4.2-foss-2015a +# This script download the s4 restarts from ECMWF +# and launch the interpolation of the files toward +# the target resolution +# History : Chloe Prodhomme - Initial version 2016 +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# Arguments +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +YEARI=2015 # first years of the restarts to be downloaded +YEARF=2016 #2016 # Last year to download +MEMBER_LST="0 1 2 3 4" # List of members to download +MON_LST="1031" #"0131 0228 0331 0430 0531 0630 0731 0831 1031 1130 1231" # List of months to download in formate MMDD +config=1L42 #resolution for the output grid +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +workdir=/scratch/Earth/$USER/ORAS4/ +rm -f /scratch/Earth/$USER/ORAS4/* +mkdir -p $workdir +cd $workdir + +DATASOURCE="ec:/emos/OCEA/4/0001/" # Root of the path where to find ORAS4 on ECFS +#loop over the members +for mem in $MEMBER_LST ; do + #loop over the month + for mon in $MON_LST ; do + #loop over the years + for YEAR in `seq $YEARI $YEARF`; do + + #name of the output file + localdir=/esnas/releases/ic/ocean/ORCA1/s4_original_file/ + fileout=s4_fc${mem}_${YEAR}${mon}_restart.nc + #check if the file already on the /esnas + if [ ! -e $localdir/${fileout}.gz ]; then + #get the exact file name at ECMWF + file=$(ecaccess-file-dir ${DATASOURCE}/opa${mem}/restart/${YEAR}/0001_${YEAR}${mon}_*_restart.tar.gz) + #download the file + ecaccess-file-get ${DATASOURCE}/opa${mem}/restart/${YEAR}/$file $file + #unzip and untar + gunzip $file + tar -xvf ${file/.gz/} + #remove not needed files + rm -f *dat ${file/.gz/} + #store the file + mv 0001_${YEAR}${mon}_*_restart.nc s4_fc${mem}_${YEAR}${mon}_restart.nc + cp s4_fc${mem}_${YEAR}${mon}_restart.nc ${localdir}/$fileout + gzip ${localdir}/$fileout + else + #copy the file into the workdir + cp $localdir/${fileout}.gz . + #unzip the file + gunzip ${fileout}.gz + + fi + #Launch the interpolation toward the target resolution + python /home/Earth/cprodhom/tools/oras4_restart/prep_restart_interp.py ${fileout} /esnas/autosubmit/con_files/mesh_mask_nemo.nemovar_O1L42.nc $config "Ec2.3" + # remove the local files + rm -rf ${fileout}.gz + done #loop year + done #loop month +done #loop member + diff --git a/Oce_IC/library_interp_extrap.py b/Oce_IC/library_interp_extrap.py new file mode 100644 index 0000000..d1a628e --- /dev/null +++ b/Oce_IC/library_interp_extrap.py @@ -0,0 +1,398 @@ +import numpy as np +from numpy import ma as ma +from netCDF4 import Dataset as nc +import sys +from scipy.spatial import ckdtree +from glob import glob +import shutil +import gzip + +def getvarmask(varin): + """define if the name of the mask corresponding to each variable + input: + varin: name of the input variables + output: + mask: name of the corresponding mask variables + varlon: name of the longitude variable + varlat: name of the latitude variables""" + + if varin in ['avmu','un','ub','gru','gsu','gtu','ssu_m','u_io']: + mask='umask' ; varlon='glamu' ; varlat='gphiu' + elif varin in ['avmv','vn','vb','grv','gsv','gtv','ssv_m','v_io']: + mask='vmask' ; varlon='glamv' ; varlat='gphiv' + elif varin in ['rotb','rotn']: + mask='fmask' ; varlon='glamf' ; varlat='gphif' + elif varin in ['avt','avmt','en','dissl','hdivb','hdivn','rhop','sn','tn','sb','tb','gcx','gcxb','gcxbfs','ssh_m','sshb','sshn','sss_m','sst_m','alb_ice','sst_io','sss_io']: + mask='tmask' ; varlon='glamt' ; varlat='gphit' + else: + print("variable not found") + sys.exit(1) + return(mask, varlon, varlat) + +def read_files(fin, varin, fmin, fmout): + """This function is reading the variables needed for the interpolation + process from the input files + Arguments: + filein: Input restart type netcdf (already open) + varin: Variable to interpolate + varmask: name of the mask variable corresponding to varin + (tmask, umask, vmask, fmask) + fmin: Input maskfile type netcdf (already open) + fmout: Output mask file type netcdf (already open) + Outputs: + fieldin: Field to interpolate + depthin: 3D depth corresponding to fieldin + thickin: 3D thichness corresponding to the fieldin + thickout: 3D thickness toward which the interpolation + has to be performed + depthout: 3D depth toward which the interpolation + has to be performed + """ + + #get the mask variable corresponding to varin + varmask, varlon, varlat =getvarmask(varin) + + + #read the input mask (needed because original s4 not masked) + maskin=np.array(fmin.variables[varmask][:].squeeze()) + #generate a masked array to mask the land + maskin=ma.array(maskin>=1,mask=maskin<0.5) + #read the thickness + thickin=np.array(fmin.variables["e3t_0"][:].squeeze()) + #read the depth + depthin=np.array(fmin.variables["gdept_0"][:].squeeze()) + #read the longitude + lonin=np.array(fmin.variables[varlon][:].squeeze()) + #read the latitude + latin=np.array(fmin.variables[varlat][:].squeeze()) + + + fieldin=np.array(fin.variables[varin][:].squeeze()) + #diemnsion name + dimnames=list(fin.variables[varin].dimensions) + + + + #if the last level of the input field is fully masked () + #or almost fully masked (less than 20points) + #remove it + #check if the variable is 2D or 3D: + if len(fieldin.shape)==3: + lastlev=maskin.shape[0]-1 + if (np.sum(maskin[lastlev,:,:])<=20) or (ma.sum(maskin[lastlev,:,:]).mask): + maskin=maskin[0:lastlev,:,:] + fieldin=np.array(fieldin[0:lastlev,:,:]) + thickin=thickin[0:lastlev] + depthin=depthin[0:lastlev] + + + #read the mask field + maskout=np.array(fmout.variables[varmask][:].squeeze()) + #thickness + thickout=np.array(fmout.variables["e3t_0"][:].squeeze()) + #depth + depthout=np.array(fmout.variables["gdept_0"][:].squeeze()) + #read the longitude + lonout=np.array(fmout.variables[varlon][:].squeeze()) + #read the latitude + latout=np.array(fmout.variables[varlat][:].squeeze()) + #generate a masked array to mask the land + maskout=ma.array(maskout,mask=maskout<0.5) + + #permute the dimensions if needed + # not take into account t dimension (size 1) + if len(fieldin.shape)==3: + correctdimorder=(u'z', u'y', u'x') + else: + correctdimorder=(u'y', u'x') + + dimnames.remove(u't') + #the correctdimension order + for idim,dimname in enumerate(dimnames): + if dimname!=correctdimorder[idim]: + print("Dimension"+dimname+"permuted") + fieldin=moveaxis(fieldin, idim, + correctdimorder.index(dimname)) + + # extend the input depth in 3D + #(needed because from Nemo 3.6 both vary with time ans space) + lz0, ly0, lx0 = maskin.shape + lzout, lyout, lxout = maskout.shape + if len(depthin.shape)==1: + depthin.shape=(lz0,1,1) + depthin=np.repeat(np.repeat(depthin,ly0, axis=1), lx0, axis=2) + if len(depthout.shape)==1: + depthout.shape=(lzout,1,1) + depthout=np.repeat(np.repeat(depthout, ly0, axis=1), lx0, axis=2) + # extend the input thickness in 3D + if len(thickin.shape)==1: + thickin.shape=(lz0,1,1) + thickin=np.repeat(np.repeat(thickin,ly0, axis=1), lx0, axis=2) + if len(thickout.shape)==1: + thickout.shape=(lzout,1,1) + thickout=np.repeat(np.repeat(thickout,ly0, axis=1), lx0, axis=2) + + + #construct mask and thickness for 2D variables + if len(fieldin.shape)==2: + ly0,lx0=fieldin.shape + maskin=maskin[0,:,:] + thickin=np.zeros((ly0,lx0),'d') + depthin=np.zeros((ly0,lx0),'d') + thickout=np.zeros((ly0,lx0),'d') + depthout=np.zeros((ly0,lx0),'d') + maskout=maskout[0,:,:] + + #check that the mask corresponds to the input file + #(in s4 restarts continents are filled with 0) + + if (fieldin[np.where(maskin.mask)]!=0).any(): + print("the variable: "+varmask) + print("the input mask:") + print(fmin.filepath()) + print("does not corespond to the input file:") + print(fin.filepath()) + print("for the variable:"+varin) + sys.exit(1) + + #apply the mask + fieldin=ma.array(fieldin,mask=1-maskin) + + #check that the field is not empty + if (np.sum(fieldin)==0): + print("The field written from the file: "+f) + print("for the variable"+varin) + print("mask with the file: "+fmask) + print("seems to be empty") + print("Please check your input") + sys.exit(1) + + return(fieldin, thickin, depthin, lonin, latin, + maskout, thickout, depthout, lonout, latout) + +def interp_vert(fieldin, thick_in, thick_out): + """This function interpolates vertically a (x,y,z) field using the + conservative method. + Arguments: + fieldin: Input Field type mask array + thick_int: input thickness thickness of the level + thick_out: output thickness of the level + History : Chloe Prodhomme - Initial version adpated + from Cfutools extrap_vert function + written by V. Guemas - 2016""" + + (lz0,ly0,lx0)=fieldin.shape + (lz1,ly1,lx1)=thick_out.shape + + var0_a=np.array(fieldin) + var0=np.array(fieldin) + var0=ma.masked_array(var0,mask=None,fill_value=None) + + var0=var0.astype('d') + + mask_in=1-fieldin.mask + + var1=np.zeros((lz1,ly0,lx0),'d') + thick_jk1=np.zeros((ly0,lx0),'d') + mask_out=np.zeros((lz1,ly0,lx0),'d') + thick_bis=np.zeros((ly0,lx0),'d') + + jk0=0 + for jk1 in np.arange(lz1) : + while ( (jk0 <= (lz0-1)) & ((np.add.reduce(thick_in[0:(jk0+1)]) <= np.add.reduce(thick_out[0:(jk1+1)])) ).any()) : + thick_bis=np.multiply(np.minimum(thick_in[jk0],np.add.reduce(thick_in[0:(jk0+1)])-np.add.reduce(thick_out[0:jk1])),mask_in[jk0,:,:]) + var1[jk1,:,:]=var1[jk1,:,:]+np.multiply(var0[jk0,:,:],thick_bis) + thick_jk1=thick_jk1+thick_bis + jk0=jk0+1 + if True: #( jk0 <= (lz0-1)) : + jk0in=min(jk0,lz0-1) + thick_bis=np.multiply(np.minimum(thick_out[jk1],np.add.reduce(thick_out[0:(jk1+1)])-np.add.reduce(thick_in[0:jk0])),mask_in[jk0in,:,:]) + var1[jk1,:,:]=var1[jk1,:,:]+np.multiply(var0[jk0in,:,:],thick_bis) + thick_jk1=thick_jk1+thick_bis + mask_out[jk1,:,:]=np.where(thick_jk1>0.01,False,True) + thick_jk1=np.where(thick_jk1==0.,1.,thick_jk1) + var1[jk1,:,:]=var1[jk1,:,:]/thick_jk1 + thick_jk1[:,:]=0 + #mask the output field + fieldout=ma.array(var1, mask=mask_out) + + return(fieldout) + + +def lon_lat_to_cartesian(lon, lat, R=1): + """calculates lon, lat coordinates of a point on a sphere with + radius R + """ + lon_r = np.radians(lon) + lat_r = np.radians(lat) + x = R * np.cos(lat_r) * np.cos(lon_r) + y = R * np.cos(lat_r) * np.sin(lon_r) + z = R * np.sin(lat_r) + return x, y, z + + +def extrap(fieldin,lonin,latin,maskout): + """This function extrapolates horizontally each level of a (x,y,z) + field using the nearest neighbour method. + Argument: + fieldin: Input field (ma.array type) + lonin: 2d longitude + latin: 2d longitude + maskout: mask of the field after extrapolation + Output: + fieldout: output field (ma.array) + """ + #add a dimension of size 1 for 2d field + if len(fieldin.shape)==2: + fieldin.shape=(1,)+fieldin.shape + + #construct points to fill + Pfill=ma.array((fieldin.mask)&(maskout!=0), mask=1-maskout) + + # fieldout intitialization + fieldout=np.array(fieldin) + + + #each level has a different mask, loop over the levels + lz0=fieldin.shape[0] + + for lev in np.arange(lz0): + newvar = ma.array(fieldin[lev,:,:], + mask=fieldin.mask[lev,:,:]) + Pfilllev = Pfill[lev,:,:] + + #not masked coordinates for the present level + nomlon = lonin[np.where(newvar.mask==False)] + nomlat = latin[np.where(newvar.mask==False)] + #cartesian coordinates of not masked points + levxs, levys, levzs = lon_lat_to_cartesian(nomlon, nomlat) + levkdt = ckdtree.cKDTree(np.array([levxs, levys, levzs]).T) + #points to extrapolate at the current level + newindexes = np.where(Pfilllev==1) + latidxs = newindexes[0] + lonidxs = newindexes[1] + #print(newindexes) + exlon = lonin[latidxs, lonidxs] + exlat = latin[latidxs, lonidxs] + #print(newvar[levkdt.data[0,0], levkdt.data[0,1]]) + + #cartesian coordinates of the points to extrapolate + xt, yt, zt = lon_lat_to_cartesian(exlon, exlat) + #calculates distances and indexes of the nearest neighbours + dist, idx = levkdt.query(np.array([xt, yt, zt]).T, k=1, p=2) + + #assign to the result matrix the found values + fieldout[lev,latidxs, lonidxs] = newvar[np.where(newvar.mask==False)].flatten()[idx] + + #mask field out + fieldout=fieldout*maskout + fieldout=ma.array(fieldout,mask=1-maskout) + return(fieldout) + + + +def vertinterp_extrap(filein, varin, fmask, fmaskout, month="", modvers=""): + """read and interpolate verticaly and extrapolate horizontal the given variable + from a restart file + Arguments: + filein: Input restart file + varin: name of the variable to extrapolate + fmask: mask file corresponding to varin + fmaskout: mask to define the grid toward which the extrapolation if done + month: starting month of the restart + (needed only if there are points to fill in the caspienne sea) + modelvers: name of the model version to find the climatology to fill the caspienne + (needed only if there are points to fill in the caspienne sea) + + """ + + #read files and define variables needed for the rest of the process + (fieldin, thickin, depthin, lonin, latin, + maskout, + thickout, depthout, + lonout, latout) = read_files(filein, + varin, fmask, fmaskout) + + #check if the variable is 2D or 3D: + if len(fieldin.shape)==3: + #perform the vertical extrapolation + fieldin=interp_vert(fieldin, thickin,thickout) + + #construct points to fill + Pfill=ma.array((fieldin.mask)&(maskout!=0), mask=1-maskout) + + #find if the caspienne sea need to be filled + caspbox=np.where((lonin>45)&(lonin<56)&(latin>35)&(latin<50)) + #longitude and latitude boundary of the caspienne sea + loncasp1,loncasp2=caspbox[0].min(),caspbox[0].max() + latcasp1,latcasp2=caspbox[1].min(),caspbox[1].max() + #check is there are points to fill in the caspienne sea + if len(fieldin.shape)==2: + casptest=np.sum(Pfill[loncasp1:loncasp2,latcasp1:latcasp2], + axis=(0,1))!=0 + else: + casptest=np.sum(Pfill[:,loncasp1:loncasp2,latcasp1:latcasp2], + axis=(0,1,2))!=0 + #If there are points to be filled + if casptest: + + if ((varin=="sn")|(varin=="sb"))|((varin=="tn")|(varin=="tb")): + #Ugly solution, to be improve would need levitus file at all resolution + if modvers=="Ec3.1": + pathclim="/esnas/releases/ic/ocean/ORCA1L46/levitus_clim/" + if ((varin=="tn")|(varin=="tb")): + #name of the file + varclim="votemper" + fileclim=pathclim+"temperature_m%02d.nc"%month + + if ((varin=="sn")|(varin=="sb")): + #open climatological files + varclim="vosaline" + fileclim=pathclim+"salinity_m%02d.nc"%month + + elif modvers=="Ec2.3": + #in this case we use an old restart for the climatology + fileclim=glob("/esnas/releases/ic/ocean/ORCA1/s4/s4_fc0_1993%02d*_restart.nc.gz"%month) + print(fileclim) + shutil.copy(fileclim[0], "fileclim.nc.gz") + fileclim="fileclim.nc" + inF = gzip.open('fileclim.nc.gz', 'rb') + outF = open(fileclim, 'wb') + outF.write( inF.read() ) + inF.close() + outF.close() + varclim=varin + else: + print("no climatology is available for this model version to fill the caspienne sea:") + print(modvers) + sys.exit(1) + + #open the file + print(fileclim) + fclim=nc(fileclim, "r") + #read variable for the given month + clim=np.array(fclim.variables[varclim][:].squeeze()) + + #add the caspienne sea + fieldin[:,caspbox[0],caspbox[1]]=clim[:,caspbox[0],caspbox[1]] + else: + print(fieldin.shape) + print(caspbox[1].min()) + print(caspbox[1].max()) + print(caspbox[0].min()) + print(caspbox[0].max()) + if (len(fieldin.shape)==3): + fieldin[:,caspbox[0],caspbox[1]]=0 + elif (len(fieldin.shape)==2): + fieldin[caspbox[0],caspbox[1]]=0 + else: + print("vertinterp_extrap function can handle only 2 and 3d arrays") + sys.exit(1) + + #horizontal extrapolation + fieldout=extrap(fieldin,lonin,latin,maskout) + return(fieldout) + + + + diff --git a/Oce_IC/prep_restart_interp.py b/Oce_IC/prep_restart_interp.py new file mode 100644 index 0000000..21a7eb8 --- /dev/null +++ b/Oce_IC/prep_restart_interp.py @@ -0,0 +1,179 @@ +from library_interp_extrap import * +from compute_density import compute_density +from glob import glob +import os +import shutil +import gzip +#This fonction generate a restart interpolated verticaly on the output grid +########### +#Arguments# +########### +#filein: input restart file +#fmask: input mask corresponding to the input file +#ex: fmask="/esnas/autosubmit/con_files/mesh_mask_nemo.nemovar_O1L42.nc" +#resout: output resolution ex: resout="1L75" +#versout: corresponding model version ex:versout="Ec3.2" + + +# Input arguments +filein=sys.argv[1] +fmask=sys.argv[2] +resout=sys.argv[3] +versout=sys.argv[4] + +#go in the working directory +#workdir="/scratch/Earth/"+os.environ["USER"]+"/restart_"+resout+"_"+versout +#if (not os.path.isdir(workdir)): +# os.makedirs(workdir) +#os.chdir(workdir) + +# mask file describing the output grid +fmaskout="/esnas/autosubmit/con_files/mesh_mask_nemo."+versout+"_O"+resout+".nc" + +dirout="/esnas/releases/ic/ocean/ORCA"+resout+"/s4" + +#print("before entering in the loop") +#check if the output file is already there +if len(glob(os.path.join(dirout, filein+".gz")))==0: + #print("enter agin in the loop") + #copy the file from the /esnas into the + dirin,fileinloc=os.path.split(filein) + if dirin!="": + shutil.copy(filein, fileinloc) + + #copy the mask from the /esnas into the + dirmask,fmaskloc=os.path.split(fmask) + if dirmask!="": + shutil.copy(fmask, fmaskloc) + + #temporary output file + fileout=filein.replace(".nc","_new.nc") + + #define the variable name of longitude latitude depth + #to be copied from the mask + vardim=["nav_lon", "nav_lat", "nav_lev"] + + #variable copied with no interpolation + varnointerp=["time_counter", "rdt", "slim", "rdttra1", "kt", "ndastp", "adatrj"] + + lstfieldout=[] + + #Open the output file + fin=nc(fileinloc, "r") + #print(fmaskout) + fmout=nc(fmaskout, "r") + fmin=nc(fmaskloc, "r") + fout=nc(fileout, "w", format="NETCDF3_CLASSIC") + nlev=fmout.dimensions["z"].size + #print(fileout) + + #copy all th dimension for fin except the level + for name, dimension in fin.dimensions.iteritems(): + if name != "z": + fout.createDimension(str(name), len(dimension) if not dimension.isunlimited() else None) + #create the level dimension + fout.createDimension("z", nlev) + #print(fileout) + + #write the variable corresponding to the dimension (from the input mask) + for name in vardim: + variable=fmout.variables[name] + #print(fileout) + #print(fout) + fout.createVariable(name, variable.datatype, variable.dimensions) + fout.variables[name][:]=variable[:] + #print("after writting dimension variable") + + listdens=[] + listdens_varname=[] + #print(fin.variables) + + + for varin, variable in fin.variables.iteritems(): + ##take out the variable nav_lev, lon ,lat + if not(varin in vardim): + ##for the variable which no need interpolation copy them + if varin in varnointerp: + #print(variable.datatype) + #print(variable.dimensions) + #print(fout) + var=fout.createVariable(varin, variable.datatype, variable.dimensions) + var[:]=variable[:] + #print("after writting variable:") + #print(varin) + else: + month=int(fileinloc.split("_")[2][4:6]) + ##perform the interpolation vertical y horizontal + fieldout=vertinterp_extrap(fin, varin, fmin, fmout, + month=month, modvers=versout) + ##add the temporal dimension + shapeout=(1,)+fieldout.shape + fieldout.shape=shapeout + print(varin) + if varin in ["tn","sn"]: + print("Iam here") + listdens.append(fieldout) + listdens_varname.append(varin) + ##put 0 for mask value (only by converting from ma into np array) + fieldout=np.array(fieldout) + #print("before defining variable") + #print(varin) + #print(variable.datatype) + #print(variable.dimensions) + #print(fout) + #add a test to avoid negative values in en + if varin =="en": + fieldout[fieldout<0]=0 + print("negative values removed from en") + var=fout.createVariable(varin, variable.datatype, variable.dimensions) + var[:]=fieldout + #print("after writting variable:") + #print(varin) + #get temperature + print(listdens_varname) + ptemp=listdens[listdens_varname.index("tn")] + #get salinity + psal=listdens[listdens_varname.index("sn")] + + #print("rhop") + rhop=compute_density(ptemp, psal) + #convert into numpy array (0 instead of mask) + rhop=np.array(rhop) + #write density into the output file + #print("before writting variable:") + #print(rhop) + var=fout.createVariable("rhop", + fin.variables["tn"].datatype, + fin.variables["tn"].dimensions) + ##print("after defining rhop in the file") + + var[:]=rhop + #print("after writting rhop") + + + fin.close() + #print("after closing fin") + fout.close() + #print("after closing fout") + fmin.close() + #print("after closing fmin") + fmout.close() + #print("after closing fmout") + + #gzip the files + + #print("before gzip: "+ fileout) + with open(fileout, 'rb') as f_in, gzip.open(fileout+'.gz', 'wb') as f_out: + shutil.copyfileobj(f_in, f_out) + #print("after gzip") + #copy the output file into /esnas + shutil.copy(fileout+'.gz', os.path.join(dirout,fileinloc+'.gz')) + print("the output restart has been created and copied into:") + print(os.path.join(dirout,fileinloc+'.gz')) + + #remove the local files + os.remove(fileinloc) + os.remove(fileout+'.gz') + os.remove(fileout) + + print("local files have been removed") -- GitLab