Changeset 900 in lmdz_wrf

Timestamp:

Jun 19, 2016, 6:09:01 PM (9 years ago)

Author:

lfita

Message:

Adding `pinterp': Function to vertically interpolate using subroutines from the p_interp.F90 NCAR program

Location:

trunk/tools

Files:

• module_ForInterpolate.F90 (modified) (3 diffs)
• nc_var.py (modified) (2 diffs)
• nc_var_tools.py (modified) (2 diffs)

trunk/tools/module_ForInterpolate.F90

@@ -923 +923 @@
 
 
- SUBROUTINE interp (data_out, data_in, pres_field, interp_levels, psfc, ter, tk, qv, ix, iy, iz, it, &
-                     num_metgrid_levels, LINLOG, extrapolate, GEOPT, MISSING)
+ SUBROUTINE interp (data_in, pres_field, interp_levels, psfc, ter, tk, qv, LINLOG, extrapolate,       &
+                     GEOPT, MISSING, data_out, ix, iy, iz, it, num_metgrid_levels)
 ! Interpolation subroutine from the p_interp.F90 NCAR program
 !  Program to read wrfout data and interpolate to pressure levels
@@ -930 +930 @@
 !  November 2007 - Cindy Bruyere
 !
-     INTEGER                                          :: ix, iy, iz, it
-     INTEGER                                          :: num_metgrid_levels, LINLOG
-     REAL, DIMENSION(ix, iy, num_metgrid_levels, it)  :: data_out
-     REAL, DIMENSION(ix, iy, iz, it)                  :: data_in, pres_field, tk, qv
-     REAL, DIMENSION(ix, iy, it)                      :: psfc
-     REAL, DIMENSION(ix, iy)                          :: ter
-     REAL, DIMENSION(num_metgrid_levels)              :: interp_levels
-
-     INTEGER                                          :: i, j, itt, k, kk, kin
-     REAL, DIMENSION(num_metgrid_levels)              :: data_out1D
-     REAL, DIMENSION(iz)                              :: data_in1D, pres_field1D
-     INTEGER                                          :: extrapolate
-     REAL                                             :: MISSING
-     REAL, DIMENSION(ix, iy, num_metgrid_levels, it)  :: N
-     REAL                                             :: sumA, sumN, AVE_geopt
-     LOGICAL                                          :: GEOPT
+     INTEGER, INTENT(IN)                                 :: ix, iy, iz, it
+     INTEGER, INTENT(IN)                                 :: num_metgrid_levels, LINLOG
+     REAL, DIMENSION(ix,iy,iz,it), INTENT(IN)            :: data_in, pres_field, tk, qv
+     REAL, DIMENSION(ix,iy,it), INTENT(IN)               :: psfc
+     REAL, DIMENSION(ix,iy), INTENT(IN)                  :: ter
+     REAL, DIMENSION(num_metgrid_levels), INTENT(IN)     :: interp_levels
+     INTEGER, INTENT(IN)                                 :: extrapolate
+     REAL, INTENT(IN)                                    :: MISSING
+     LOGICAL, INTENT(IN)                                 :: GEOPT
+     REAL, DIMENSION(ix,iy,num_metgrid_levels,it),                                                    &
+       INTENT(OUT)                                       :: data_out
+
+! Local
+     INTEGER                                             :: i, j, itt, k, kk, kin
+     REAL, DIMENSION(num_metgrid_levels)                 :: data_out1D
+     REAL, DIMENSION(iz)                                 :: data_in1D, pres_field1D
+     REAL, DIMENSION(ix, iy, num_metgrid_levels, it)     :: N
+     REAL                                                :: sumA, sumN, AVE_geopt
 
 !!!!!!! Variables
@@ -1172 +1174 @@
  END SUBROUTINE int1D
 
+ FUNCTION virtual (tmp,rmix)
+!      This function returns virtual temperature in K, given temperature
+!      in K and mixing ratio in kg/kg.
+
+     real                              :: tmp, rmix, virtual
+
+     virtual=tmp*(0.622+rmix)/(0.622*(1.+rmix))
+
+ END FUNCTION virtual
+
 END MODULE module_ForInterpolate

trunk/tools/nc_var.py

@@ -41 +41 @@
   'ncreplace', 'ncstepdiff', 'netcdf_concatenation', 'netcdf_fold_concatenation',    \
   'Partialmap_Entiremap', 'Partialmap_EntiremapFor', 'Partialmap_EntiremapForExact', \
-  'remapnn',                                                                         \
+  'pinterp', 'remapnn',                                                              \
   'seasmean', 'sellonlatbox', 'sellonlatlevbox', 'selvar', 'setvar_asciivalues',     \
   'sorttimesmat', 'spacemean', 'SpatialWeightedMean', 'statcompare_files',           \
@@ -218 +218 @@
 elif oper == 'opersvarsfiles':
     ncvar.compute_opersvarsfiles(opts.values, opts.varname)
+elif oper == 'pinterp':
+    ncvar.pinterp(opts.values, opts.ncfile, opts.varname)
 elif oper == 'remapnn':
     ncvar.remapnn(opts.values, opts.ncfile, opts.varname)

trunk/tools/nc_var_tools.py

@@ -81 +81 @@
 # Partialmap_EntiremapFor: Function to transform from a partial global map (e.g.: only land points) to an entire one usinf Fortran code
 # Partialmap_EntiremapForExact: Function to transform from a partial global map (e.g.: only land points) to an entire one using Fortran code with exact location
+# pinterp: Function to vertically interpolate using subroutines from the p_interp.F90 NCAR program
 # put_variable_slice: Function to add a variable from a netcdf object to a new one following a slice
 # remapnn: Function to remap to the nearest neightbor a variable using projection from another file
@@ -15793 +15794 @@
     return
 
+def pinterp(values, ncfile, variables):
+    """ Function to vertically interpolate using subroutines from the p_interp.F90 NCAR program
+      Using fortran codes: module_generic.F90
+      foudre: f2py -m module_ForInterpolate --f90exec=/usr/bin/gfortran-4.7 -c module_generic.F90 module_ForInterpolate.F90  >& run_f2py.log
+      ciclad: f2py --f90flags="-fPIC" -L/opt/canopy-1.3.0/Canopy_64bit/System/lib/ -L/usr/lib64/ -L/opt/canopy-1.3.0/Canopy_64bit/System/lib/ -m module_ForInterpolate -c module_generic.F90 module_ForInterpolate.F90 
+      values= [interplevs],[linloginterp],[extrap]
+        [interplevs]: ':' separated list of pressure values to interpolate [hPa] (decreassing in hPa)
+        [linloginterp]: kind of linear interpolation
+          1: linear interp in pressure
+          2: linear interp in ln(pressure)
+        [extrap]: whether to set to missing value below/above model ground and top (0), or extrapolate (1)
+      ncfile= WRF file to use
+      variables = ',' separated list of names of variables to interpolate ('all', fo all 4D-atmospheric variables)
+        'WRFght': for WRF geopotential height
+        'WRFrh': for WRF relative humidity
+        'WRFt': for WRF temperature
+    """
+    import module_ForInterpolate as fin
+    fname = 'pinterp'
+
+    if values == 'h':
+        print fname + '_____________________________________________________________'
+        print pinterp.__doc__
+        quit()
+
+    arguments = '[interplevs],[linloginterp],[extrap]'
+    gen.check_arguments(fname, values, arguments, ',')
+
+    interplevs = np.array(values.split(',')[0].split(':'), dtype=np.float64)
+    linloginterp = np.int32(values.split(',')[1])
+    extrap = np.int32(values.split(',')[2])
+
+    ofile = 'pinterp.nc'
+
+    onc = NetCDFFile(ncfile, 'r')
+
+    # Variables to interpolate
+    WRFdims = ['Time', 'bottom_top', 'south_north', 'west_east']
+    newWRFdims = ['Time', 'pres', 'south_north', 'west_east']
+    notCHK = ['WRFght', 'WRFt', 'WRFrh']
+
+    varns = gen.str_list(variables, ',')
+    if variables == 'all':
+        varns = []
+        for vn in onc.variables:
+            ovn = onc.variables[vn]
+            coincdims = set(WRFdims) & set(ovn.dimensions)
+           
+            if len(WRFdims) == len(ovn.dimensions) and len(coincdims) == len(WRFdims):
+                varns.append(vn)
+        varns = varns + notCHK
+        
+    for vn in varns:
+        if not gen.searchInlist(onc.variables, vn) and                               \
+          not gen.searchInlist(notCHK, vn):
+            print errormsg
+            print '  ' + fname + "': WRF file '" + ncfile +                          \
+              "' does not have variable '" + vn + "' !!"
+            quit(-1)
+
+    # looking for WRF required variables
+    WRFvarrequired = ['P', 'PB', 'PSFC', 'PH', 'PHB', 'HGT', 'T', 'QVAPOR', 'XLONG', \
+      'XLAT', 'Times' , 'P_TOP']
+    for var in WRFvarrequired:
+        if not gen.searchInlist(onc.variables, var):
+            print errormsg
+            print '  ' + fname + ": file '" + ncfile + "' does not have variable '" +\
+              var + "' !!"
+            quit(-1)
+
+# Computing necessary variables
+    # pressure
+    ovar1 = onc.variables['P']
+    ovar2 = onc.variables['PB']
+
+    pres = (ovar1[:] + ovar2[:]).astype('float64')
+
+    dimx = pres.shape[3]
+    dimy = pres.shape[2]
+    dimz = pres.shape[1]
+    dimt = pres.shape[0]
+
+    # sfc pressure
+    ovar1 = onc.variables['PSFC']
+    psfc = ovar1[:].astype('float64')
+
+    # geopotential height
+    ovar1 = onc.variables['PH']
+    ovar2 = onc.variables['PHB']
+    geop0 = ovar1[:] + ovar2[:]
+    unstg = list(geop0.shape)
+    unstg[1] = unstg[1] - 1
+    geop = np.zeros(tuple(unstg), dtype=np.float)
+    geop = 0.5*(geop0[:,1:dimz,:,:] + geop0[:,0:dimz-1,:,:]).astype('float64')
+
+    # terrain height
+    ovar1 = onc.variables['HGT']
+    hgt = ovar1[0,:,:].astype('float64')
+
+    # water vapour micing ratio
+    ovar1 = onc.variables['QVAPOR']
+    qv = ovar1[:].astype('float64')
+
+    # temperature
+    Rd = 287.04
+    Cp = 7.*Rd/2.
+    RCP = Rd/Cp
+    p0 = 100000.
+    ovar1 = onc.variables['T']
+    temp = ((ovar1[:]+300.)*(pres[:]/p0)**RCP).astype('float64')
+
+    onewnc = NetCDFFile(ofile, 'w')
+# Creation of dimensions
+    newdim = onewnc.createDimension('west_east', dimx)
+    newdim = onewnc.createDimension('south_north', dimy)
+    newdim = onewnc.createDimension('pres', len(interplevs))
+    newdim = onewnc.createDimension('Time', None)
+
+# Creation of variable dimensions
+    WRFvardims = ['XLONG', 'XLAT', 'Times']
+    for var in WRFvardims:
+        ovar = onc.variables[var]
+        varinf = variable_inf(ovar)
+        for vd in varinf.dimns:
+            if not gen.searchInlist(onewnc.dimensions, vd):
+                newdim = onewnc.createDimension(vd, len(onc.dimensions[vd]))
+        newvar = onewnc.createVariable(var, nctype(varinf.dtype),                    \
+          tuple(varinf.dimns))
+        newvar[:] = ovar[:]
+        for attrn in ovar.ncattrs():
+            attrv = ovar.getncattr(attrn)
+            attr = set_attribute(newvar, attrn, attrv)    
+
+    # Creation of pressure variable dimension
+    newvar = onewnc.createVariable('pres', 'f8', ('pres'))
+    newvar[:] = interplevs * 100.
+    basicvardef(newvar, 'pressure', 'Pressure', 'Pa')
+    attr = set_attribute(newvar, 'positive', 'down')
+    attr = set_attribute(newvar, 'axis', 'Z')
+
+    print '  ' + fname + ' interpolating at pressure levels _______' 
+    print '     ', interplevs
+    for vn in varns:
+        print '  ' + fname + ": interpolaion of '" + vn + "' ..."
+        newvarattr = {}
+        varin = None
+
+        if gen.searchInlist(notCHK, vn):
+            if vn == 'WRFght': 
+                varin = geop
+                isgeopt = True
+                varattrs = gen.variables_values('WRFght')
+                CFvn = varattrs[0]
+                newvarattr['standard_name'] = varattrs[1]
+                newvarattr['long_name'] = varattrs[4].replace('|',' ')
+                newvarattr['units'] = varattrs[5]
+            elif vn == 'WRFrh': 
+                # relative humidity
+                ovar1 = 10.*0.6112*exp(17.67*(temp-273.16)/(temp-29.65))
+                ovar2 = 0.622*ovar1/(0.01*pres-(1.-0.622)*ovar1)
+                varin = 100.*qv/ovar2
+                isgeop = False
+                varattrs = gen.variables_values('WRFrh')
+                CFvn = varattrs[0]
+                newvarattr['standard_name'] = varattrs[1]
+                newvarattr['long_name'] = varattrs[4].replace('|',' ')
+                newvarattr['units'] = varattrs[5]
+            elif vn == 'WRFt':
+                varin = temp
+                isgeop = False
+                varattrs = gen.variables_values('WRFt')
+                CFvn = varattrs[0]
+                newvarattr['standard_name'] = varattrs[1]
+                newvarattr['long_name'] = varattrs[4].replace('|',' ')
+                newvarattr['units'] = varattrs[5]
+        elif not gen.searchInlist(WRFvarrequired, vn) or vn == 'QVAPOR':
+            ovarin = onc.variables[vn]
+            varinf = variable_inf(ovarin)
+            varin = ovarin[:]
+            isgeop = False
+            varattrs = gen.variables_values(vn)
+            CFvn = varattrs[0]
+            newvarattr['standard_name'] = varattrs[1]
+            newvarattr['long_name'] = varattrs[4].replace('|',' ')
+            newvarattr['units'] = varattrs[5]
+            for attrn in ovarin.ncattrs():
+                attrv = ovarin.getncattr(attrn)
+                newvarattr[attrn] = attrv
+
+        if varin is not None:
+            varint = varin.transpose()
+            prest = pres.transpose()
+            psfct = psfc.transpose()
+            hgtt = hgt.transpose()
+            tempt = temp.transpose()
+            qvt = qv.transpose()
+            varinterp = fin.module_forinterpolate.interp( data_in=varint,            \
+              pres_field=prest, interp_levels=interplevs, psfc=psfct,                \
+              ter=hgtt, tk=tempt, qv=qvt, linlog=linloginterp,                       \
+              extrapolate=extrap, geopt=isgeop, missing=gen.fillValueF)
+
+            newvar = onewnc.createVariable(CFvn, 'f4', tuple(newWRFdims),            \
+              fill_value=gen.fillValueF)
+            newvar[:] = varinterp.transpose()
+
+            for attrn in newvarattr.keys():
+                attr = set_attribute(newvar, attrn, newvarattr[attrn])
+            onewnc.sync()
+
+    # Global attributes
+
+    onewnc.setncattr('author', 'L. Fita')
+    newattr = set_attributek(onewnc, 'institution', unicode('Laboratoire de M' +     \
+      unichr(233) + 't' + unichr(233) + 'orologie Dynamique'), 'U')
+    onewnc.setncattr('university', 'Pierre Marie Curie - Jussieu')
+    onewnc.setncattr('center', 'Centre National de Recherches Scientifiques')
+    onewnc.setncattr('city', 'Paris')
+    onewnc.setncattr('country', 'France')
+    onewnc.setncattr('script', 'nc_var_tools.py')
+    onewnc.setncattr('function', 'pinterp')
+    onewnc.setncattr('version', '1.0')
+    onewnc.setncattr('oringal', 'Interpolation using p_interp.F90 subroutines')
+    newattr = set_atributek(onewnc, 'original_subroutines_author', unicode('Cindy ' +\
+      'Bruy' + unichr(232) + 're'), 'U')
+    onewnc.setncattr('original_subroutines_date', 'November 2007')
+
+    for attrn in onc.ncattrs():
+        attrv = onc.getncattr(attrn)
+        attr = set_attribute(onewnc, attrn, attrv)    
+
+    onewnc.sync()
+    onewnc.close()
+
+    print fname + ": successfull written of '" + ofile + "' !!"
+
+    return
+
 #quit()