source: trunk/UTIL/PYTHON/mcs.py @ 403

#!/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"])