#!/usr/bin/env python ### A. Colaitis ## # This routine transforms a diagfi.nc file into a diagfi_MCS.nc file where # the fields are directly comparable to those contained in MCS data, i.e. # fields are re-binned at times over the ranges specified in the MCS file. ### ########################################################################################### ########################################################################################### ### What is below relate to running the file as a command line executable (very convenient) if __name__ == "__main__": import sys from optparse import OptionParser ### to be replaced by argparse from netCDF4 import Dataset from os import system from times import sol2ls import numpy as np from mymath import find_nearest from myplot import getfield from make_netcdf import make_gcm_netcdf parser = OptionParser() ############################# ### Options parser.add_option('-f', '--file', action='store',dest='file', type="string", default=None, help='[NEEDED] filename.') parser.add_option('-m', '--mfile', action='store',dest='mcsfile', type="string", default=None, help='[NEEDED] filename for MCS comparison.') ############################# ### Get options and variables (opt,args) = parser.parse_args() ############################# ### Load and check data # Files nc=Dataset(opt.file) ncmcs=Dataset(opt.mcsfile) # Dimensions lon=nc.variables["longitude"][:] lat=nc.variables["latitude"][:] alt=nc.variables["altitude"][:] time=nc.variables["Time"][:] # in fraction of sols controle=nc.variables["controle"][:] day_ini=controle[3]%669 time[:]=time[:]+day_ini nx=len(lon) ny=len(lat) nz=len(alt) nt=len(time) lstime=sol2ls(time) # MCS dtimemintmp=ncmcs.variables["dtimemin"][:,:,:] dtimemaxtmp=ncmcs.variables["dtimemax"][:,:,:] ntimemintmp=ncmcs.variables["ntimemin"][:,:,:] ntimemaxtmp=ncmcs.variables["ntimemax"][:,:,:] lonmcs=ncmcs.variables["longitude"][:] latmcs=ncmcs.variables["latitude"][:] timemcs=ncmcs.variables["time"][:]%360 # IN LS dtimemin=np.ma.masked_where(dtimemintmp < 0.,dtimemintmp) dtimemin.set_fill_value([np.NaN]) dtimemax=np.ma.masked_where(dtimemaxtmp < 0.,dtimemaxtmp) dtimemax.set_fill_value([np.NaN]) ntimemin=np.ma.masked_where(ntimemintmp < 0.,ntimemintmp) ntimemin.set_fill_value([np.NaN]) ntimemax=np.ma.masked_where(ntimemaxtmp < 0.,ntimemaxtmp) ntimemax.set_fill_value([np.NaN]) # Variables temp=getfield(nc,"temp") dimensions = np.array(temp).shape print "temp" print dimensions #aps=nc.variables["aps"][:] #bps=nc.variables["bps"][:] ############################# ### Building ############################# ### We loop over chunks of gcm data corresponding to MCS time dimension ### Bin sizes for mcs data is 5 degrees ls centered on value varday=np.zeros([len(timemcs),nz,ny,nx]) varnight=np.zeros([len(timemcs),nz,ny,nx]) vardayout=np.ma.masked_invalid(varday) varnightout=np.ma.masked_invalid(varnight) i=0 for ls in timemcs: lsstart=ls-2.5 lsstop=ls+2.5 istart=find_nearest(lstime,lsstart,strict=True) istop=find_nearest(lstime,lsstop,strict=True) if ((istart is np.NaN) or (istop is np.NaN)): vardayout[i,:,:,:]=np.NaN varnightout[i,:,:,:]=np.NaN print "Time interval skipped. Ls MCS: (",lsstart,';',lsstop,')',"// Ls diagfi: (",lstime.min(),';',lstime.max(),')' i=i+1 continue print "--->> Processing Data. Ls MCS: (",lsstart,';',lsstop,')',"// Ls diagfi: (",lstime.min(),';',lstime.max(),')' varchk=temp[istart:istop,:,:,:] varchktime=time[istart:istop] ndays=np.floor(varchktime[len(varchktime)-1])-np.floor(varchktime[0]) dtmichk=dtimemin[i,:,:] dtmachk=dtimemax[i,:,:] ntmichk=ntimemin[i,:,:] ntmachk=ntimemax[i,:,:] dtmichk.set_fill_value(np.NaN) dtmachk.set_fill_value(np.NaN) ntmichk.set_fill_value(np.NaN) ntmachk.set_fill_value(np.NaN) dtmichk=dtmichk.filled() dtmachk=dtmachk.filled() ntmichk=ntmichk.filled() ntmachk=ntmachk.filled() ### We iterate for each day in the chunk, on each grid point we find ### the closest corresponding MCS grid point and the index of the ### time in the chunk closest to the time in the closest MCS grid point. ### (yea it's complicated) vartmpnight=np.zeros([ndays,nz,ny,nx]) vartmpday=np.zeros([ndays,nz,ny,nx]) nd=0 while nd < ndays: daystart=find_nearest(varchktime-varchktime[0],nd) daystop=find_nearest(varchktime-varchktime[0],nd+1) # varchktime_lon=np.zeros([daystop-daystart+1,len(lon)]) ix=0 for x in lon: varchktime_lon = 24.*(varchktime[daystart:daystop+1]-varchktime[daystart]) + x/15. iy=0 for y in lat: niy=find_nearest(latmcs,y) nix=find_nearest(lonmcs,x) nitdtmichk=find_nearest(varchktime_lon,dtmichk[niy,nix]) nitdtmachk=find_nearest(varchktime_lon,dtmachk[niy,nix]) nitntmichk=find_nearest(varchktime_lon,ntmichk[niy,nix]) nitntmachk=find_nearest(varchktime_lon,ntmachk[niy,nix]) if ((nitdtmichk is np.NaN) or (nitdtmachk is np.NaN)): vartmpday[nd,:,iy,ix]=np.NaN elif nitdtmichk > nitdtmachk: vartmpday[nd,:,iy,ix]=(np.mean(varchk[daystart+nitdtmichk:daystop+1,:,iy,ix],axis=0)+np.mean(varchk[daystart:daystart+nitdtmachk+1,:,iy,ix],axis=0))/2. else: vartmpday[nd,:,iy,ix]=np.mean(varchk[daystart+nitdtmichk:daystart+nitdtmachk+1,:,iy,ix],axis=0) if ((nitntmichk is np.NaN) or (nitntmachk is np.NaN)): vartmpnight[nd,:,iy,ix]=np.NaN elif nitntmichk > nitntmachk: vartmpnight[nd,:,iy,ix]=(np.mean(varchk[daystart+nitntmichk:daystop+1,:,iy,ix],axis=0)+np.mean(varchk[daystart:daystart+nitntmachk+1,:,iy,ix],axis=0))/2. else: vartmpnight[nd,:,iy,ix]=np.mean(varchk[daystart+nitntmichk:daystart+nitntmachk+1,:,iy,ix],axis=0) iy=iy+1 ix=ix+1 nd=nd+1 vartmpdaymasked=np.ma.masked_invalid(vartmpday) vartmpnightmasked=np.ma.masked_invalid(vartmpnight) vardayout[i,:,:,:]=np.ma.mean(vartmpdaymasked[:,:,:,:],axis=0) varnightout[i,:,:,:]=np.ma.mean(vartmpnightmasked[:,:,:,:],axis=0) i=i+1 vardayout=np.ma.masked_invalid(vardayout) varnightout=np.ma.masked_invalid(varnightout) vardayout.set_fill_value([np.NaN]) varnightout.set_fill_value([np.NaN]) vardayncdf=vardayout.filled() varnightncdf=varnightout.filled() make_gcm_netcdf (zfilename="diagfi_MCS.nc", \ zdescription="Temperatures from diagfi reworked to match MCS format", \ zlon=lon, \ zlat=lat, \ zalt=alt, \ ztime=timemcs, \ zvariables=[vardayncdf,varnightncdf], \ znames=["dtemp","ntemp"])